diff --git a/examples/transient_examples/short_transient_search_gridded.py b/examples/transient_examples/short_transient_search_gridded.py index fd9b83b2f85d9e076231042efbe1c0a8b9b495a5..64e642342f40b088a5c68c1d44c7b69e68be20e9 100644 --- a/examples/transient_examples/short_transient_search_gridded.py +++ b/examples/transient_examples/short_transient_search_gridded.py @@ -35,7 +35,8 @@ search1 = pyfstat.GridSearch( sftfilepattern=os.path.join(datadir,'*simulated_transient_signal*sft'), F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref, minStartTime=minStartTime, maxStartTime=maxStartTime, - BSGL=False) + BSGL=False, + outputTransientFstatMap=True) search1.run() search1.print_max_twoF() @@ -49,7 +50,8 @@ search2 = pyfstat.GridSearch( F0s=F0s, F1s=F1s, F2s=F2s, Alphas=Alphas, Deltas=Deltas, tref=tref, minStartTime=minStartTime, maxStartTime=maxStartTime, transientWindowType='rect', t0Band=Tspan-2*Tsft, tauBand=Tspan, - BSGL=False) + BSGL=False, + outputTransientFstatMap=True) search2.run() search2.print_max_twoF() diff --git a/pyfstat/core.py b/pyfstat/core.py index 5d06781e156ea0c700fe4561deee3f1da2a122b6..e336acc3bc8d893d1b9193b3f5f002f90b00f3c2 100755 --- a/pyfstat/core.py +++ b/pyfstat/core.py @@ -642,7 +642,6 @@ class ComputeFstat(BaseSearchClass): asini=None, period=None, ecc=None, tp=None, argp=None): """ Returns twoF or ln(BSGL) fully-coherently at a single point """ - self.PulsarDopplerParams.fkdot = np.array([F0, F1, F2, 0, 0, 0, 0]) self.PulsarDopplerParams.Alpha = Alpha self.PulsarDopplerParams.Delta = Delta @@ -681,10 +680,11 @@ class ComputeFstat(BaseSearchClass): # F-stat computation self.windowRange.tau = int(2*self.Tsft) - FS = lalpulsar.ComputeTransientFstatMap( + self.FstatMap = lalpulsar.ComputeTransientFstatMap( self.FstatResults.multiFatoms[0], self.windowRange, False) + F_mn = self.FstatMap.F_mn.data - twoF = 2*np.max(FS.F_mn.data) + twoF = 2*np.max(F_mn) if self.BSGL is False: if np.isnan(twoF): return 0 @@ -705,7 +705,7 @@ class ComputeFstat(BaseSearchClass): # to return BSGL # FIXME: should we instead compute BSGL over the whole F_mn # and return the maximum of that? - idx_maxTwoF = np.argmax(FS.F_mn.data) + idx_maxTwoF = np.argmax(F_mn) self.twoFX[0] = 2*FS0.F_mn.data[idx_maxTwoF] self.twoFX[1] = 2*FS1.F_mn.data[idx_maxTwoF] @@ -751,10 +751,11 @@ class ComputeFstat(BaseSearchClass): self.transientWindowType = 'none' self.init_computefstatistic_single_point() for tau in taus: - twoFs.append(self.get_fullycoherent_twoF( + detstat = self.get_fullycoherent_twoF( tstart=tstart, tend=tstart+tau, F0=F0, F1=F1, F2=F2, Alpha=Alpha, Delta=Delta, asini=asini, period=period, ecc=ecc, - tp=tp, argp=argp)) + tp=tp, argp=argp) + twoFs.append(detstat) return taus, np.array(twoFs) diff --git a/pyfstat/grid_based_searches.py b/pyfstat/grid_based_searches.py index f0e67a35afdc765be935eb267901f5602808e94e..46179360b742ee51627482b554afc4d544435f99 100644 --- a/pyfstat/grid_based_searches.py +++ b/pyfstat/grid_based_searches.py @@ -32,7 +32,8 @@ class GridSearch(BaseSearchClass): nsegs=1, BSGL=False, minCoverFreq=None, maxCoverFreq=None, detectors=None, SSBprec=None, injectSources=None, input_arrays=False, assumeSqrtSX=None, - transientWindowType=None, t0Band=None, tauBand=None): + transientWindowType=None, t0Band=None, tauBand=None, + outputTransientFstatMap=False): """ Parameters ---------- @@ -59,6 +60,9 @@ class GridSearch(BaseSearchClass): and tau in (2*Tsft,2*Tsft+tauBand). if =0, only compute CW Fstat with t0=minStartTime, tau=maxStartTime-minStartTime. + outputTransientFstatMap: bool + if true, write output files for (t0,tau) Fstat maps + (one file for each doppler grid point!) For all other parameters, see `pyfstat.ComputeFStat` for details """ @@ -157,8 +161,15 @@ class GridSearch(BaseSearchClass): data = [] for vals in tqdm(self.input_data): - FS = self.search.get_det_stat(*vals) - data.append(list(vals) + [FS]) + detstat = self.search.get_det_stat(*vals) + windowRange = getattr(self.search, 'windowRange', None) + FstatMap = getattr(self.search, 'FstatMap', None) + data.append(list(vals) + [detstat]) + if self.outputTransientFstatMap and self.transientWindowType: + tCWfile = os.path.splitext(self.out_file)[0]+'_tCW_%.16f_%.16f_%.16f_%.16g_%.16g.dat' % (vals[2],vals[5],vals[6],vals[3],vals[4]) # freq alpha delta f1dot f2dot + fo = lal.FileOpen(tCWfile, 'w') + lalpulsar.write_transientFstatMap_to_fp ( fo, FstatMap, windowRange, None ) + del fo # instead of lal.FileClose() which is not SWIG-exported data = np.array(data, dtype=np.float) if return_data: