Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-b246709
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
HierarchSearchGCT.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2011-2013 Karl Wette.
3 * Copyright (C) 2009-2010 Holger Pletsch.
4 *
5 * Based on HierarchicalSearch.c by
6 * Copyright (C) 2005-2008 Badri Krishnan, Alicia Sintes, Bernd Machenschalk.
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with with program; see the file COPYING. If not, write to the
20 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
21 * MA 02110-1301 USA
22 *
23 */
24
25/*********************************************************************************/
26/**
27 * \defgroup lalpulsar_bin_GCT GCT Search Application
28 * \ingroup lalpulsar_bin_Apps
29 */
30
31/**
32 * \author Holger Pletsch
33 * \file
34 * \ingroup lalpulsar_bin_GCT
35 * \brief Hierarchical semicoherent CW search code based on F-Statistic,
36 * exploiting global-correlation coordinates (Phys.Rev.Lett. 103, 181102, 2009)
37 *
38 */
39
40/* ---------- Includes -------------------- */
41#include <lal/Segments.h>
42#include <lal/LALString.h>
43#include <lal/LineRobustStats.h>
44#include <RecalcToplistStats.h>
45
46#include "HierarchSearchGCT.h"
47
48#ifdef GC_SSE2_OPT
49#include <gc_hotloop_sse2.h>
50#else
51#define ALRealloc LALRealloc
52#define ALFree LALFree
53#endif
54
55/* ---------- Defines -------------------- */
56/* #define DIAGNOSISMODE 1 */
57#define NUDGE 10*LAL_REAL8_EPS
58
59#define TRUE (1==1)
60#define FALSE (1==0)
61
62#ifdef __GNUC__
63#define UNUSED __attribute__ ((unused))
64#else
65#define UNUSED
66#endif
67
68#define GET_GCT_CHECKPOINT read_gct_checkpoint // (cptname, semiCohToplist, NULL, &count)
69#define SET_GCT_CHECKPOINT write_gct_checkpoint
70char **global_argv;
71int global_argc;
72
73#define FSTART 100.0 /**< Default Start search frequency */
74#define FBAND 0.0 /**< Default search band */
75#define FDOT 0.0 /**< Default value of first spindown */
76#define DFDOT 0.0 /**< Default range of first spindown parameter */
77#define F2DOT 0.0 /**< Default value of second spindown */
78#define F3DOT 0.0 /**< Default value of third spindown */
79#define DF2DOT 0.0 /**< Default range of second spindown parameter */
80#define DF3DOT 0.0 /**< Default range of third spindown parameter */
81#define SKYREGION "allsky" /**< default sky region to search over -- just a single point*/
82#define MISMATCH 0.3 /**< Default for metric grid maximal mismatch value */
83#define DALPHA 0.001 /**< Default resolution for isotropic or flat grids */
84#define DDELTA 0.001 /**< Default resolution for isotropic or flat grids */
85#define FSTATTHRESHOLD 2.6 /**< Default threshold on Fstatistic for peak selection */
86#define NCAND1 10 /**< Default number of candidates to be followed up from first stage */
87#define FNAMEOUT "./HS_GCT.out" /**< Default output file basename */
88#ifndef LAL_INT4_MAX
89#define LAL_INT4_MAX 2147483647
90#endif
91
92#define BLOCKSIZE_REALLOC 50
93
94#define Vorb_GCT = 2.9785e04;
95#define Vspin_GCT = 465.10;
96#define REARTH_GCT = 6.378140e06;
97#define C_GCT = 299792458;
98
99/* ---------- Macros -------------------- */
100#define HSMAX(x,y) ( (x) > (y) ? (x) : (y) )
101#define HSMIN(x,y) ( (x) < (y) ? (x) : (y) )
102
103#define GETTIME() (uvar_outputTiming ? XLALGetCPUTime() : 0)
104//#define GETTIME() (uvar_outputTiming ? XLALGetTimeOfDay() : 0)
105
106/* ---------- Exported types ---------- */
107/** useful variables for each hierarchical stage */
108typedef struct {
109 CHAR *sftbasename; /**< filename pattern for sfts */
110 LIGOTimeGPS tStartGPS; /**< start and end time of stack */
111 REAL8 tObs; /**< tEndGPS - tStartGPS */
112 REAL8 refTime; /**< reference time for pulsar params */
113 PulsarSpinRange spinRange_startTime; /**< freq and fdot range at start-time of observation */
114 PulsarSpinRange spinRange_endTime; /**< freq and fdot range at end-time of observation */
115 PulsarSpinRange spinRange_refTime; /**< freq and fdot range at the reference time */
116 PulsarSpinRange spinRange_midTime; /**< freq and fdot range at mid-time of observation */
117 EphemerisData *edat; /**< ephemeris data for XLALBarycenter */
118 LIGOTimeGPSVector *midTstack; /**< timestamps vector for mid time of each stack */
119 LIGOTimeGPSVector *startTstack; /**< timestamps vector for start time of each stack */
120 LIGOTimeGPSVector *endTstack; /**< timestamps vector for end time of each stack */
121 LIGOTimeGPS minStartTimeGPS; /**< all sft data must be after this time */
122 LIGOTimeGPS maxStartTimeGPS; /**< all sft timestamps must be before this GPS time */
123 UINT4 blocksRngMed; /**< blocksize for running median noise floor estimation */
124 UINT4 Dterms; /**< size of Dirichlet kernel for Fstat calculation */
125 UINT4 DtermsRecalc; /**< Recalc: size of Dirichlet kernel for Fstat calculation */
126 LALStringVector *assumeSqrtSX; /**< Assume stationary Gaussian noise with detector noise-floors sqrt{SX}" */
127 /* parameters describing the coherent data-segments */
128 REAL8 tStack; /**< duration of stacks */
129 UINT4 nStacks; /**< number of stacks */
130 LALSegList *segmentList; /**< parsed segment list read from user-specified input file --segmentList */
131 BSGLSetup *BSGLsetup; /**< pre-computed setup for line-robust statistic BSGL */
132 REAL8 dFreqStack; /**< frequency resolution of Fstat calculation */
133 REAL8 df1dot; /**< coarse grid resolution in spindown */
134 REAL8 df2dot; /**< coarse grid resolution in 2nd spindown */
135 REAL8 df3dot; /**< coarse grid resolution in 3rd spindown */
136 UINT4 extraBinsFstat; /**< Extra Fstat frequency bins required to cover residual spindowns */
137 UINT4 binsFstatSearch; /**< nominal number of Fstat frequency bins in search band */
138 UINT4 nf1dot; /**< number of 1st spindown Fstat bins */
139 UINT4 nf2dot; /**< number of 2nd spindown Fstat bins */
140 UINT4 nf3dot; /**< number of 3rd spindown Fstat bins */
141 int SSBprec; /**< SSB transform precision */
142 FstatMethodType Fmethod; //!< which Fstat-method/algorithm to use
143 BOOLEAN recalcToplistStats; //!< do additional analysis for all toplist candidates, output F, FXvector for postprocessing */
144 FstatMethodType FmethodRecalc; //!< which Fstat-method/algorithm to use for the recalc step
145 REAL8 mismatch1; /**< 'mismatch1' user-input needed here internally ... */
146 UINT4 nSFTs; /**< total number of SFTs */
147 LALStringVector *detectorIDs; /**< vector of detector IDs */
148 REAL4 NSegmentsInvX[PULSAR_MAX_DETECTORS]; /**< effective inverse number of segments per detector (needed for correct averaging in single-IFO F calculation) */
149 FstatInputVector *Fstat_in_vec; /**< Original wide-parameter search: vector of Fstat input data structures for XLALComputeFstat(), one per stack */
150 FstatInputVector *Fstat_in_vec_recalc; /**< Recalculate the toplist: Vector of Fstat input data structures for XLALComputeFstat(), one per stack */
151 PulsarParamsVector *injectionSources; ///< Source parameters to inject: comma-separated list of file-patterns and/or direct config-strings ('{...}')
152 BOOLEAN collectFstatTiming; ///< flag whether to collect and output F-stat timing info
153} UsefulStageVariables;
154
155
156/**
157 * Struct holding various timing measurements and relevant search parameters.
158 * This is used to fit timing-models with measured times to predict search run-times
159 */
160typedef struct {
161 UINT4 Nseg; ///< number of semi-coherent segments
162 UINT4 Ndet; ///< number of detectors
163 UINT4 Tcoh; ///< length of coherent segments in seconds
164 UINT4 Nsft; ///< total number of SFTs
165 UINT4 Ncand; ///< length of toplists
166
167 UINT4 NFreqCo; ///< total number of frequency bins computed in coarse grid (including sidebands!)
168 REAL8 Ncoh; ///< number of coarse-grid Fstat templates ('coherent')
169 REAL8 Ninc; ///< number of fine-grid templates ('incoherent')
170
171 const char *FstatMethodStr; ///< Fstat-method used
172 const char *RecalcMethodStr; ///< Fstat-method used
173
174 // ----------
175 // extended timing model:
176 // runtime = Nseg * Ndet * Ncoh * tauF + Nseg * Ninc * tauSumFine + Ninc * tauExtraStats + Ntop * tauRecalc + time_Other
177 REAL8 tau_Fstat; //< time to compute F-stat for one coarse-grid template, one detector, one segment
178 REAL8 tau_SumF; //< time to sum F-stat values of one segment for one fine-grid point
179 REAL8 tau_Bayes; //< time to compute final Bayes-factor statistics for one fine-grid point
180 REAL8 tau_Recalc; //< time to recalc one toplist-entry
181 REAL8 time_Other; //< all non-scaling leftovers in overall timeing
182 // ----------
183
184} timingInfo_t;
185
186
187/* ------------------------ Functions -------------------------------- */
188void SetUpSFTs( LALStatus *status, UsefulStageVariables *in );
189void PrintFstatVec( LALStatus *status, FstatResults *in, FILE *fp, PulsarDopplerParams *thisPoint,
190 LIGOTimeGPS refTime, INT4 stackIndex );
191void PrintCatalogInfo( LALStatus *status, const SFTCatalog *catalog, FILE *fp );
192void PrintStackInfo( LALStatus *status, const SFTCatalogSequence *catalogSeq, FILE *fp );
193void UpdateSemiCohToplists( LALStatus *status, toplist_t *list1, toplist_t *list2, toplist_t *list3, FineGrid *in, REAL8 f1dot_fg, REAL8 f2dot_fg, REAL8 f3dot_fg, UsefulStageVariables *usefulparams, REAL4 NSegmentsInv, REAL4 *NSegmentsInvX, BOOLEAN have_f3dot );
194
195void UpdateSemiCohToplistsOptimTriple( LALStatus *status,
196 SortBy_t toplist_sortby,
197 toplist_t *list1,
198 toplist_t *list2,
199 toplist_t *list3,
200 FineGrid *in,
201 REAL8 f1dot_fg,
202 REAL8 f2dot_fg,
203 REAL8 f3dot_fg,
204 UsefulStageVariables *usefulparams,
205 REAL4 NSegmentsInv,
206 REAL4 *NSegmentsInvX,
207 BOOLEAN have_f3dot
208 );
209
210void GetSegsPosVelAccEarthOrb( LALStatus *status, REAL8VectorSequence **posSeg,
211 REAL8VectorSequence **velSeg, REAL8VectorSequence **accSeg,
212 UsefulStageVariables *usefulparams );
213static inline INT4 ComputeU1idx( REAL8 freq_event, REAL8 f1dot_event, REAL8 A1, REAL8 B1, REAL8 U1start, REAL8 U1winInv );
214void ComputeU2idx( REAL8 freq_event, REAL8 f1dot_event, REAL8 A2, REAL8 B2, REAL8 U2start, REAL8 U2winInv,
215 INT4 *U2idx );
216int compareCoarseGridUindex( const void *a, const void *b );
217int compareFineGridNC( const void *a, const void *b );
218int compareFineGridsumTwoF( const void *a, const void *b );
219
220SFTCatalogSequence *XLALSetUpStacksFromSegmentList( const SFTCatalog *catalog, const LALSegList *segList );
221
222int XLALExtrapolateToplistPulsarSpins( toplist_t *list,
223 const LIGOTimeGPS usefulParamsRefTime,
224 const LIGOTimeGPS finegridRefTime );
225
226static int write_TimingInfo( const CHAR *fname, const timingInfo_t *ti );
227static inline REAL4 findLoudestTwoF( const FstatResults *in );
228
229/* ---------- Global variables -------------------- */
230LALStatus *global_status; /* a global pointer to main()s head of the LALStatus structure */
231char *global_column_headings_stringp;
232
233// XLALReadSegmentsFromFile(): applications which still must support
234// the deprecated 4-column format should set this variable to non-zero
235extern int XLALReadSegmentsFromFile_support_4column_format;
236
237/* ################################### MAIN ################################### */
238
239int main( int argc, char *argv[] )
240{
241 LALStatus XLAL_INIT_DECL( status );
242
243 /* temp loop variables: generally k loops over segments and j over SFTs in a stack */
244 UINT4 k;
245 UINT4 skyGridCounter; /* coarse sky position counter */
246 UINT4 f1dotGridCounter; /* coarse f1dot position counter */
247
248 /* GPS timestamp vectors */
249 LIGOTimeGPSVector *midTstack = NULL;
250 LIGOTimeGPSVector *startTstack = NULL;
251 LIGOTimeGPSVector *endTstack = NULL;
252
253 /* General GPS times */
254 LIGOTimeGPS XLAL_INIT_DECL( refTimeGPS );
255 LIGOTimeGPS XLAL_INIT_DECL( tMidGPS );
256
257 /* GPS time used for each segment's midpoint */
258 LIGOTimeGPS XLAL_INIT_DECL( midTstackGPS );
259 REAL8 timeDiffSeg; /* Difference to tMidGPS (midpoint) of Tobs */
260
261 /* pos, vel, acc at midpoint of segments */
262 REAL8VectorSequence *posStack = NULL;
263 REAL8VectorSequence *velStack = NULL;
264 REAL8VectorSequence *accStack = NULL;
265
266 /* duration of each segment */
267 REAL8 tStack;
268
269 /* number of segments */
270 UINT4 nStacks;
271
272 /* Total observation time */
273 REAL8 tObs;
274
275 /* SFT related stuff */
276 static LIGOTimeGPS minStartTimeGPS, maxStartTimeGPS;
277
278 /* some useful variables for each stage */
279 UsefulStageVariables XLAL_INIT_DECL( usefulParams );
280
281 /* F-statistic computation related stuff */
282 FstatResults *Fstat_res = NULL; // Pointer to Fstat results structure, will be allocated by XLALComputeFstat()
283 FstatQuantities Fstat_what = FSTATQ_2F; // Quantities to be computed by XLALComputeFstat()
284 UINT4 binsFstat1;
285
286 /* Semicoherent variables */
287 static SemiCoherentParams semiCohPar;
288
289 /* coarse grid */
290 CoarseGrid XLAL_INIT_DECL( coarsegrid );
291 REAL8 dFreqStack; /* frequency resolution of Fstat calculation */
292 REAL8 df1dot; /* coarse grid resolution in spindown */
293 UINT4 ifdot; /* counter for coarse-grid spindown values */
294 REAL8 df2dot; /* coarse grid resolution in 2nd spindown */
295 UINT4 if2dot; /* counter for coarse-grid 2nd spindown values */
296 REAL8 df3dot; /* coarse grid resolution in 3rd spindown */
297 UINT4 if3dot; /* counter for coarse-grid 3rd spindown values */
298
299 /* fine grid */
300 FineGrid finegrid;
301 UINT4 nf1dots_fg = 1; /* number of frequency and spindown values */
302 REAL8 gammaRefine, sigmasq; /* refinement factor and variance */
303 UINT4 nf2dots_fg = 1; /* number of second spindown values */
304 REAL8 gamma2Refine, sigma4; /* 2nd spindown refinement factor and 4th moment */
305 UINT4 nf3dots_fg = 1; /* number of third spindown values */
306 REAL8 gamma3Refine = 1; /* 3rd spindown refinement */
307
308 /* GCT helper variables */
309 UINT4 if1dot_fg, if2dot_fg, if3dot_fg;
310 UINT4 ifreq;
311 INT4 U1idx;
312 REAL8 myf0, freq_event, f1dot_event;
313 REAL8 dfreq_fg, df1dot_fg, freqmin_fg, f1dotmin_fg, freqband_fg;
314 REAL8 df2dot_fg, f2dotmin_fg;
315 REAL8 df3dot_fg, f3dotmin_fg;
316 REAL8 u1start, u1win, u1winInv;
317 REAL8 freq_fg, f1dot_fg, f2dot_fg, f3dot_fg, f1dot_event_fg;
318 REAL8 A1, B1;
319 // currently unused: REAL8 A2;
320 REAL8 B2; /* GCT helper variables for faster calculation of u1 or u2 */
321 REAL8 pos[3];
322 REAL8 vel[3];
323 // currently unused: REAL8 acc[3];
324 REAL8 cosAlpha, sinAlpha, cosDelta, sinDelta;
325 REAL8 nvec[3]; /* unit vector pointing to sky position */
326
327 /* These vars are currently not used, but eventually in the future.
328 INT4 U2idx, NumU2idx;
329 REAL8 myf0max, u2start, u2end;
330 REAL8 u2win, u2winInv;
331 */
332
333 /* fstat candidate structure for candidate toplist*/
334 toplist_t *semiCohToplist = NULL;
335 toplist_t *semiCohToplist2 = NULL; // only used for SORTBY_DUAL_F_BSGL, SORTBY_TRIPLE_BStSGLtL or SORTBY_F_BSGLtL_BtSGLtL
336 toplist_t *semiCohToplist3 = NULL; // only used for SORTBY_TRIPLE_BStSGLtL and SORTBY_F_BSGLtL_BtSGLtL
337
338 /* template and grid variables */
339 static DopplerSkyScanInit scanInit; /* init-structure for DopperScanner */
340 DopplerSkyScanState XLAL_INIT_DECL( thisScan ); /* current state of the Doppler-scan */
341 static PulsarDopplerParams dopplerpos; /* current search-parameters */
342 static PulsarDopplerParams thisPoint;
343 UINT4 oldcg = 0, oldfg = 0;
344
345 /* temporary storage for spinrange vector */
346 static PulsarSpinRange spinRange_Temp;
347
348 /* variables for logging */
349 CHAR *fnamelog = NULL;
350 FILE *fpLog = NULL;
351 CHAR *logstr = NULL;
352
353 /* output candidate files and file pointers */
354 CHAR *fnameSemiCohCand = NULL;
355 CHAR *fnameFstatVec1 = NULL;
356 FILE *fpFstat1 = NULL;
357
358 /* checkpoint filename */
359 CHAR *uvar_fnameChkPoint = NULL;
360
361 /* user variables */
362 BOOLEAN uvar_log = FALSE; /* logging done if true */
363
364 BOOLEAN uvar_printCand1 = FALSE; /* if 1st stage candidates are to be printed */
365 BOOLEAN uvar_printFstat1 = FALSE;
366 BOOLEAN uvar_loudestTwoFPerSeg = FALSE; // output loudest per-segment Fstat candidates
367 BOOLEAN uvar_semiCohToplist = TRUE; /* if overall first stage candidates are to be output */
368
369 LALStringVector *uvar_assumeSqrtSX = NULL; /* Assume stationary Gaussian noise with detector noise-floors sqrt{SX}" */
370
371 BOOLEAN uvar_recalcToplistStats = FALSE; /* Do additional analysis for all toplist candidates, output F, FXvector for postprocessing */
372 BOOLEAN uvar_loudestSegOutput = FALSE; /* output extra info about loudest segment; requires recalcToplistStats */
373
374 // ----- Line robust stats parameters ----------
375 BOOLEAN uvar_computeBSGL = FALSE; /* In Fstat loop, compute line-robust statistic (BSGL=log10BSGL) using single-IFO F-stats */
376 BOOLEAN uvar_BSGLlogcorr = FALSE; /* compute log-correction in line-robust statistic BSGL (slower) or not (faster) */
377 REAL8 uvar_Fstar0sc = 0.0; /* (semi-coherent) BSGL transition-scale parameter 'Fstar0sc=Nseg*Fstar0coh', see documentation for XLALCreateBSGLSetup() for details */
378 LALStringVector *uvar_oLGX = NULL; /* prior per-detector line-vs-Gauss odds ratios 'oLGX', see XLALCreateBSGLSetup() for details */
379 BOOLEAN uvar_getMaxFperSeg = FALSE; /* In Fstat loop, compute maximum F and FX over segments */
380 // --------------------------------------------
381
382 REAL8 uvar_dAlpha = DALPHA; /* resolution for flat or isotropic grids -- coarse grid*/
383 REAL8 uvar_dDelta = DDELTA;
384 REAL8 uvar_f1dot = FDOT; /* first spindown value */
385 REAL8 uvar_f1dotBand = DFDOT; /* range of first spindown parameter */
386 REAL8 uvar_f2dot = F2DOT; /* second spindown value */
387 REAL8 uvar_f2dotBand = DF2DOT; /* range of second spindown parameter */
388 REAL8 uvar_f3dot = F3DOT; /* second spindown value */
389 REAL8 uvar_f3dotBand = DF3DOT; /* range of second spindown parameter */
390 REAL8 uvar_Freq = FSTART;
391 REAL8 uvar_FreqBand = FBAND;
392
393 REAL8 uvar_dFreq = 0;
394 REAL8 uvar_df1dot = 0; /* coarse grid frequency and spindown resolution */
395 REAL8 uvar_df2dot = 0; /* coarse grid second spindown resolution */
396 REAL8 uvar_df3dot = 0; /* coarse grid third spindown resolution */
397
398 REAL8 uvar_ThrF = FSTATTHRESHOLD; /* threshold of Fstat to select peaks */
399 REAL8 uvar_mismatch1 = MISMATCH; /* metric mismatch for first stage coarse grid */
400
401 REAL8 uvar_minStartTime1 = 0;
402 REAL8 uvar_maxStartTime1 = LAL_INT4_MAX;
403
404 REAL8 uvar_refTime = 0;
405 INT4 uvar_nCand1 = NCAND1; /* number of candidates to be followed up from first stage */
406
407 INT4 uvar_blocksRngMed = FstatOptionalArgsDefaults.runningMedianWindow;
408
409 REAL8 uvar_tStack = 0;
410 INT4 uvar_nStacksMax = 1;
411 CHAR *uvar_segmentList = NULL; /**< ALTERNATIVE: file containing a pre-computed segment list of tuples (startGPS endGPS duration[h] NumSFTs) */
412
413 INT4 uvar_Dterms = FstatOptionalArgsDefaults.Dterms;
414 INT4 uvar_DtermsRecalc = FstatOptionalArgsDefaults.Dterms;
415 INT4 uvar_SSBprecision = FstatOptionalArgsDefaults.SSBprec;
416 INT4 uvar_gammaRefine = 1;
417 INT4 uvar_gamma2Refine = 1;
418 INT4 uvar_metricType1 = LAL_PMETRIC_COH_PTOLE_ANALYTIC;
419 INT4 uvar_gridType1 = GRID_METRIC;
420 INT4 uvar_skyPointIndex = -1;
421
422 CHAR *uvar_ephemEarth; /**< Earth ephemeris file to use */
423 CHAR *uvar_ephemSun; /**< Sun ephemeris file to use */
424
425 CHAR *uvar_skyRegion = NULL;
426 CHAR *uvar_fnameout = NULL;
427 CHAR *uvar_DataFiles1 = NULL;
428 CHAR *uvar_skyGridFile = NULL;
429 INT4 uvar_numSkyPartitions = 0;
430 INT4 uvar_partitionIndex = 0;
431 INT4 uvar_SortToplist = 0;
432
433 CHAR *uvar_outputTiming = NULL;
434 CHAR *uvar_outputTimingDetails = NULL;
435
436 int uvar_FstatMethod = FstatOptionalArgsDefaults.FstatMethod;
437 int uvar_FstatMethodRecalc = FstatOptionalArgsDefaults.FstatMethod;
438
439 timingInfo_t XLAL_INIT_DECL( timing );
440
441 LALStringVector *uvar_injectionSources = NULL;
442
443 // timing values
444 REAL8 tic_RecalcToplist, time_RecalcToplist = 0;
445 REAL8 tic_Fstat, time_Fstat = 0;
446 REAL8 tic_SumFine, time_SumFine = 0;
447 REAL8 tic_ExtraStats, time_ExtraStats = 0;
448 REAL8 tic_Start, time_Total = 0;
449
450 global_status = &status;
451
452 global_argv = argv;
453 global_argc = argc;
454
455 // XLALReadSegmentsFromFile(): continue to support deprecated 4-column format (startGPS endGPS duration NumSFTs, duration is ignored)
456 XLALReadSegmentsFromFile_support_4column_format = 1;
457
458 uvar_ephemEarth = XLALStringDuplicate( "earth00-40-DE405.dat.gz" );
459 uvar_ephemSun = XLALStringDuplicate( "sun00-40-DE405.dat.gz" );
460
461 uvar_skyRegion = LALCalloc( strlen( SKYREGION ) + 1, sizeof( CHAR ) );
462 strcpy( uvar_skyRegion, SKYREGION );
463
464 uvar_fnameout = LALCalloc( strlen( FNAMEOUT ) + 1, sizeof( CHAR ) );
465 strcpy( uvar_fnameout, FNAMEOUT );
466
467 /* set LAL error-handler */
468 lal_errhandler = LAL_ERR_EXIT;
469
470 /* register user input variables */
471 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_log, "log", BOOLEAN, 0, OPTIONAL, "Write log file" ) == XLAL_SUCCESS, XLAL_EFUNC );
472 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_semiCohToplist, "semiCohToplist", BOOLEAN, 0, OPTIONAL, "Print toplist of semicoherent candidates" ) == XLAL_SUCCESS, XLAL_EFUNC );
473 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_DataFiles1, "DataFiles1", STRING, 0, REQUIRED, "1st SFT file pattern. Possibilities are:\n"
474 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" ) == XLAL_SUCCESS, XLAL_EFUNC );
475 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_skyRegion, "skyRegion", STRING, 0, OPTIONAL, "sky-region polygon (or 'allsky')" ) == XLAL_SUCCESS, XLAL_EFUNC );
476 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_numSkyPartitions, "numSkyPartitions", INT4, 0, OPTIONAL, "No. of (equi-)partitions to split skygrid into" ) == XLAL_SUCCESS, XLAL_EFUNC );
477 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_partitionIndex, "partitionIndex", INT4, 0, OPTIONAL, "Index [0,numSkyPartitions-1] of sky-partition to generate" ) == XLAL_SUCCESS, XLAL_EFUNC );
478 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_skyGridFile, "skyGridFile", STRING, 0, OPTIONAL, "sky-grid file" ) == XLAL_SUCCESS, XLAL_EFUNC );
479 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dAlpha, "dAlpha", REAL8, 0, OPTIONAL, "Resolution for flat or isotropic coarse grid" ) == XLAL_SUCCESS, XLAL_EFUNC );
480 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dDelta, "dDelta", REAL8, 0, OPTIONAL, "Resolution for flat or isotropic coarse grid" ) == XLAL_SUCCESS, XLAL_EFUNC );
481 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_Freq, "Freq", REAL8, 'f', OPTIONAL, "Start search frequency" ) == XLAL_SUCCESS, XLAL_EFUNC );
482 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dFreq, "dFreq", REAL8, 0, OPTIONAL, "Frequency resolution (required if nonzero FreqBand)" ) == XLAL_SUCCESS, XLAL_EFUNC );
483 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_FreqBand, "FreqBand", REAL8, 'b', OPTIONAL, "Search frequency band" ) == XLAL_SUCCESS, XLAL_EFUNC );
484 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f1dot, "f1dot", REAL8, 0, OPTIONAL, "Spindown parameter" ) == XLAL_SUCCESS, XLAL_EFUNC );
485 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_df1dot, "df1dot", REAL8, 0, OPTIONAL, "Spindown resolution (required if nonzero f1dotBand)" ) == XLAL_SUCCESS, XLAL_EFUNC );
486 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f1dotBand, "f1dotBand", REAL8, 0, OPTIONAL, "Spindown Range" ) == XLAL_SUCCESS, XLAL_EFUNC );
487 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f2dot, "f2dot", REAL8, 0, OPTIONAL, "2nd spindown parameter" ) == XLAL_SUCCESS, XLAL_EFUNC );
488 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_df2dot, "df2dot", REAL8, 0, OPTIONAL, "2nd spindown resolution (required if nonzero f2dotBand)" ) == XLAL_SUCCESS, XLAL_EFUNC );
489 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f2dotBand, "f2dotBand", REAL8, 0, OPTIONAL, "2nd spindown Range" ) == XLAL_SUCCESS, XLAL_EFUNC );
490 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f3dot, "f3dot", REAL8, 0, OPTIONAL, "3rd spindown parameter" ) == XLAL_SUCCESS, XLAL_EFUNC );
491 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_df3dot, "df3dot", REAL8, 0, OPTIONAL, "3rd spindown resolution (required if nonzero f3dotBand)" ) == XLAL_SUCCESS, XLAL_EFUNC );
492 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f3dotBand, "f3dotBand", REAL8, 0, OPTIONAL, "3rd spindown Range" ) == XLAL_SUCCESS, XLAL_EFUNC );
493 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_ThrF, "peakThrF", REAL8, 0, OPTIONAL, "Fstat Threshold" ) == XLAL_SUCCESS, XLAL_EFUNC );
494 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_mismatch1, "mismatch1", REAL8, 'm', OPTIONAL, "1st stage mismatch" ) == XLAL_SUCCESS, XLAL_EFUNC );
495 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_gridType1, "gridType1", INT4, 0, OPTIONAL, "0=flat, 1=isotropic, 2=metric, 3=file" ) == XLAL_SUCCESS, XLAL_EFUNC );
496 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_metricType1, "metricType1", INT4, 0, OPTIONAL, "0=none, 1=Ptole-analytic, 2=Ptole-numeric, 3=exact" ) == XLAL_SUCCESS, XLAL_EFUNC );
497 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_gammaRefine, "gammaRefine", INT4, 'g', OPTIONAL, "Refinement of fine grid (default: use segment times)" ) == XLAL_SUCCESS, XLAL_EFUNC );
498 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_gamma2Refine, "gamma2Refine", INT4, 'G', OPTIONAL, "Refinement of f2dot fine grid (default: use segment times, -1=use gammaRefine)" ) == XLAL_SUCCESS, XLAL_EFUNC );
499 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_fnameout, "fnameout", STRING, 'o', OPTIONAL, "Output filename" ) == XLAL_SUCCESS, XLAL_EFUNC );
500 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_fnameChkPoint, "fnameChkPoint", STRING, 0, OPTIONAL, "Checkpoint filename" ) == XLAL_SUCCESS, XLAL_EFUNC );
501 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nCand1, "nCand1", INT4, 'n', OPTIONAL, "No. of candidates to output" ) == XLAL_SUCCESS, XLAL_EFUNC );
502 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_printCand1, "printCand1", BOOLEAN, 0, OPTIONAL, "Print 1st stage candidates" ) == XLAL_SUCCESS, XLAL_EFUNC );
503 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_refTime, "refTime", REAL8, 0, OPTIONAL, "Ref. time for pulsar pars [Default: mid-time]" ) == XLAL_SUCCESS, XLAL_EFUNC );
504 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_ephemEarth, "ephemEarth", STRING, 0, OPTIONAL, "Location of Earth ephemeris file" ) == XLAL_SUCCESS, XLAL_EFUNC );
505 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_ephemSun, "ephemSun", STRING, 0, OPTIONAL, "Location of Sun ephemeris file" ) == XLAL_SUCCESS, XLAL_EFUNC );
506 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_minStartTime1, "minStartTime1", REAL8, 0, OPTIONAL, "1st stage: Only use SFTs with timestamps starting from (including) this GPS time" ) == XLAL_SUCCESS, XLAL_EFUNC );
507 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_maxStartTime1, "maxStartTime1", REAL8, 0, OPTIONAL, "1st stage: Only use SFTs with timestamps up to (excluding) this GPS time" ) == XLAL_SUCCESS, XLAL_EFUNC );
508 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_printFstat1, "printFstat1", BOOLEAN, 0, OPTIONAL, "Print 1st stage Fstat vectors" ) == XLAL_SUCCESS, XLAL_EFUNC );
509
510 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_assumeSqrtSX, "assumeSqrtSX", STRINGVector, 0, OPTIONAL, "Don't estimate noise-floors but assume (stationary) per-IFO sqrt{SX} (if single value: use for all IFOs)" ) == XLAL_SUCCESS, XLAL_EFUNC );
511
512 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nStacksMax, "nStacksMax", INT4, 0, OPTIONAL, "Maximum No. of segments" ) == XLAL_SUCCESS, XLAL_EFUNC );
513 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_tStack, "tStack", REAL8, 'T', OPTIONAL, "Duration of segments (sec)" ) == XLAL_SUCCESS, XLAL_EFUNC );
514 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_segmentList, "segmentList", STRING, 0, OPTIONAL, "ALTERNATIVE: file containing a segment list: lines of form <startGPS endGPS duration[h] NumSFTs>" ) == XLAL_SUCCESS, XLAL_EFUNC );
515 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_recalcToplistStats, "recalcToplistStats", BOOLEAN, 0, OPTIONAL, "Additional analysis for toplist candidates, recalculate 2F, 2FX at finegrid" ) == XLAL_SUCCESS, XLAL_EFUNC );
516 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_loudestSegOutput, "loudestSegOutput", BOOLEAN, 0, OPTIONAL, "Output extra info about loudest segment; (requires --recalcToplistStats)" ) == XLAL_SUCCESS, XLAL_EFUNC );
517
518 // ----- Line robust stats parameters ----------
519 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_computeBSGL, "computeBSGL", BOOLEAN, 0, OPTIONAL, "Compute and output line-robust statistic (BSGL)" ) == XLAL_SUCCESS, XLAL_EFUNC );
520 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_Fstar0sc, "Fstar0sc", REAL8, 0, OPTIONAL, "BSGL: semi-coh transition-scale parameter 'Fstar0sc=Nseg*Fstar0coh'" ) == XLAL_SUCCESS, XLAL_EFUNC );
521 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_oLGX, "oLGX", STRINGVector, 0, OPTIONAL, "BSGL: prior per-detector line-vs-Gauss odds 'oLGX' (Defaults to oLGX=1/Ndet)" ) == XLAL_SUCCESS, XLAL_EFUNC );
522 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_BSGLlogcorr, "BSGLlogcorr", BOOLEAN, 0, DEVELOPER, "BSGL: include log-correction terms (slower) or not (faster)" ) == XLAL_SUCCESS, XLAL_EFUNC );
523 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_getMaxFperSeg, "getMaxFperSeg", BOOLEAN, 0, OPTIONAL, "Compute and output maximum F and FX over segments" ) == XLAL_SUCCESS, XLAL_EFUNC );
524 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_SortToplist, "SortToplist", INT4, 0, OPTIONAL, "Sort toplist by: 0=Fstat, 1=nc, 2=B_S/GL, 3='Fstat + B_S/GL', 4=B_S/GLtL, 5=B_tS/GLtL, 6='B_S/GL + B_S/GLtL + B_tS/GLtL', 7='Fstat + B_S/GLtL + B_tS/GLtL' " ) == XLAL_SUCCESS, XLAL_EFUNC );
525 // --------------------------------------------
526
527 XLAL_CHECK_MAIN( XLALRegisterNamedUvarAuxData( &uvar_FstatMethod, "FstatMethod", UserEnum, XLALFstatMethodChoices(), 0, OPTIONAL, "F-statistic method to use" ) == XLAL_SUCCESS, XLAL_EFUNC );
528 XLAL_CHECK_MAIN( XLALRegisterNamedUvarAuxData( &uvar_FstatMethodRecalc, "FstatMethodRecalc", UserEnum, XLALFstatMethodChoices(), 0, OPTIONAL, "F-statistic method to use for recalc" ) == XLAL_SUCCESS, XLAL_EFUNC );
529
530 /* developer user variables */
531 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_blocksRngMed, "blocksRngMed", INT4, 0, DEVELOPER, "RngMed block size" ) == XLAL_SUCCESS, XLAL_EFUNC );
532 XLAL_CHECK_MAIN( XLALRegisterNamedUvarAuxData( &uvar_SSBprecision, "SSBprecision", UserEnum, &SSBprecisionChoices, 0, DEVELOPER, "Precision for SSB transform" ) == XLAL_SUCCESS, XLAL_EFUNC );
533 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_Dterms, "Dterms", INT4, 0, DEVELOPER, "Number of kernel terms (single-sided) to use in\na) Dirichlet kernel if FstatMethod=Demod*\nb) sinc-interpolation kernel if FstatMethod=Resamp*" ) == XLAL_SUCCESS, XLAL_EFUNC );
534 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_DtermsRecalc, "DtermsRecalc", INT4, 0, DEVELOPER, "Same as 'Dterms', applies to 'Recalc' step" ) == XLAL_SUCCESS, XLAL_EFUNC );
535 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_skyPointIndex, "skyPointIndex", INT4, 0, DEVELOPER, "Only analyze this skypoint in grid" ) == XLAL_SUCCESS, XLAL_EFUNC );
536
537 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_outputTiming, "outputTiming", STRING, 0, DEVELOPER, "Append timing information into this file" ) == XLAL_SUCCESS, XLAL_EFUNC );
538 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_outputTimingDetails, "outputTimingDetails", STRING, 0, DEVELOPER, "Append detailed averaged F-stat timing information to this file" ) == XLAL_SUCCESS, XLAL_EFUNC );
539
540 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_loudestTwoFPerSeg, "loudestTwoFPerSeg", BOOLEAN, 0, DEVELOPER, "Output loudest per-segment Fstat values into file '_loudestTwoFPerSeg'" ) == XLAL_SUCCESS, XLAL_EFUNC );
541
542 /* inject signals into the data being analyzed */
543 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_injectionSources, "injectionSources", STRINGVector, 0, DEVELOPER, "%s", InjectionSourcesHelpString ) == XLAL_SUCCESS, XLAL_EFUNC );
544
545 /* read all command line variables */
546 BOOLEAN should_exit = 0;
547 XLAL_CHECK_MAIN( XLALUserVarReadAllInput( &should_exit, argc, argv, lalPulsarVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
548 if ( should_exit ) {
549 return ( 1 );
550 }
551
552 /* assemble version string */
553 CHAR *VCSInfoString;
554 if ( ( VCSInfoString = XLALVCSInfoString( lalPulsarVCSInfoList, 0, "%% " ) ) == NULL ) {
555 XLALPrintError( "XLALVCSInfoString failed.\n" );
556 return HIERARCHICALSEARCH_ESUB;
557 }
558
559 LogPrintfVerbatim( LOG_DEBUG, "Code-version: %s\n", VCSInfoString );
560 // LogPrintfVerbatim( LOG_DEBUG, "CFS Hotloop variant: %s\n", OptimisedHotloopSource );
561
562 /* some basic sanity checks on user vars */
563 if ( uvar_nStacksMax < 1 ) {
564 fprintf( stderr, "Invalid number of segments!\n" );
565 return ( HIERARCHICALSEARCH_EBAD );
566 }
567
568#ifndef EXP_NO_NUM_COUNT
569 /* check that the numbercount can't exceed the data type */
570 {
571 UINT8 maxseg = 1;
572 maxseg = maxseg << ( 8 * sizeof( FINEGRID_NC_T ) );
573 maxseg -= 1;
574 if ( ( UINT8 )uvar_nStacksMax > maxseg ) {
575 fprintf( stderr,
576 "Number of segments exceeds %" LAL_UINT8_FORMAT "!\n"
577 "Compile without GC_SSE2_OPT to extend the available segment range\n",
578 maxseg );
579 return ( HIERARCHICALSEARCH_EBAD );
580 }
581 }
582#endif
583
584 if ( uvar_blocksRngMed < 1 ) {
585 fprintf( stderr, "Invalid Running Median block size\n" );
586 return ( HIERARCHICALSEARCH_EBAD );
587 }
588
589 if ( uvar_ThrF < 0 ) {
590 fprintf( stderr, "Invalid value of F-statistic threshold\n" );
591 return ( HIERARCHICALSEARCH_EBAD );
592 }
593
594 if ( uvar_f3dotBand != 0 && ( !XLALUserVarWasSet( &uvar_gammaRefine ) || !XLALUserVarWasSet( &uvar_gammaRefine ) || uvar_gammaRefine != 1 || uvar_gamma2Refine != 1 ) ) {
595 fprintf( stderr, "Search over 3rd spindown is available only with gammaRefine AND gamma2Refine manually set to 1!\n" );
596 return ( HIERARCHICALSEARCH_EVAL );
597 }
598
599 /* 2F threshold for semicoherent stage */
600#ifndef EXP_NO_NUM_COUNT
601 REAL4 TwoFthreshold = 2.0 * uvar_ThrF;
602#endif
603
604 if ( ( uvar_SortToplist < 0 ) || ( uvar_SortToplist >= SORTBY_LAST ) ) {
605 XLALPrintError( "Invalid value %d specified for toplist sorting, must be within [0, %d]\n", uvar_SortToplist, SORTBY_LAST - 1 );
606 return ( HIERARCHICALSEARCH_EBAD );
607 }
608 if ( ( uvar_SortToplist == SORTBY_BSGL || uvar_SortToplist == SORTBY_DUAL_F_BSGL ||
609 uvar_SortToplist == SORTBY_BSGLtL || uvar_SortToplist == SORTBY_BtSGLtL ||
610 uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL ||
611 uvar_SortToplist == SORTBY_F_BSGLtL_BtSGLtL ) && !uvar_computeBSGL ) {
612 fprintf( stderr, "Toplist sorting by BSGL[tL] only possible if --computeBSGL given.\n" );
613 return ( HIERARCHICALSEARCH_EBAD );
614 }
615 if ( ( uvar_SortToplist == SORTBY_BSGLtL || uvar_SortToplist == SORTBY_BtSGLtL ||
616 uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL ||
617 uvar_SortToplist == SORTBY_F_BSGLtL_BtSGLtL ) && !uvar_getMaxFperSeg ) {
618 fprintf( stderr, "Toplist sorting by B[t]SGLtL only possible if --getMaxFperSeg given.\n" );
619 return ( HIERARCHICALSEARCH_EBAD );
620 }
621
622 /* create toplist -- semiCohToplist has the same structure
623 as a fstat candidate, so treat it as a fstat candidate */
624 if ( uvar_SortToplist == SORTBY_DUAL_F_BSGL ) { // special treatement of 'dual' toplists: 1st one sorted by 'F', 2nd one by 'BSGL'
625 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist, uvar_nCand1, SORTBY_F ),
626 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_F );
627 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist2, uvar_nCand1, SORTBY_BSGL ),
628 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BSGL );
629 } else if ( uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL ) { // special treatement of 'triple' toplists: 1st one sorted by 'B_S/GL', 2nd one by 'B_S/GLtL', 3rd by 'B_tS/GLtL'
630 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist, uvar_nCand1, SORTBY_BSGL ),
631 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BSGL );
632 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist2, uvar_nCand1, SORTBY_BSGLtL ),
633 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BSGLtL );
634 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist3, uvar_nCand1, SORTBY_BtSGLtL ),
635 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BtSGLtL );
636 } else if ( uvar_SortToplist == SORTBY_F_BSGLtL_BtSGLtL ) { // special treatement of 'triple' toplists: 1st one sorted by 'F', 2nd one by 'B_S/GLtL', 3rd by 'B_tS/GLtL'
637 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist, uvar_nCand1, SORTBY_F ),
638 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_F );
639 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist2, uvar_nCand1, SORTBY_BSGLtL ),
640 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BSGLtL );
641 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist3, uvar_nCand1, SORTBY_BtSGLtL ),
642 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BtSGLtL );
643 } else { // 'normal' single-sorting toplist cases (sortby 'F', 'nc' or 'BSGL')
644 XLAL_CHECK( 0 == create_gctFstat_toplist( &semiCohToplist, uvar_nCand1, uvar_SortToplist ),
645 XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, uvar_SortToplist );
646 }
647
648 /* write the log file */
649 if ( uvar_log ) {
650 fnamelog = LALCalloc( strlen( uvar_fnameout ) + 1 + 4, sizeof( CHAR ) );
651 strcpy( fnamelog, uvar_fnameout );
652 strcat( fnamelog, ".log" );
653 /* open the log file for writing */
654 if ( ( fpLog = fopen( fnamelog, "wb" ) ) == NULL ) {
655 fprintf( stderr, "Unable to open file %s for writing\n", fnamelog );
656 LALFree( fnamelog );
657 /*exit*/
658 return ( HIERARCHICALSEARCH_EFILE );
659 }
660
661 /* get the log string */
662 XLAL_CHECK_MAIN( ( logstr = XLALUserVarGetLog( UVAR_LOGFMT_CFGFILE ) ) != NULL, XLAL_EFUNC );
663
664 fprintf( fpLog, "# Log file for HierarchSearchGCT.c\n\n" );
665 fprintf( fpLog, "# User Input:\n" );
666 fprintf( fpLog, "#-------------------------------------------\n" );
667 fprintf( fpLog, "# cmdline: %s\n", logstr );
668 LALFree( logstr );
669
670 /* add code version ID (only useful for git-derived versions) */
671 fprintf( fpLog, "# version: %s\n", VCSInfoString );
672
673 fclose( fpLog );
674 LALFree( fnamelog );
675
676 } /* end of logging */
677
678 tic_Start = GETTIME();
679 usefulParams.collectFstatTiming = ( uvar_outputTimingDetails != NULL );
680
681 /* initializations of coarse and fine grids */
682 coarsegrid.TwoF = NULL;
683 coarsegrid.TwoFX = NULL;
684 coarsegrid.Uindex = NULL;
685 finegrid.nc = NULL;
686 finegrid.sumTwoF = NULL;
687 finegrid.sumTwoFX = NULL;
688 finegrid.maxTwoFl = NULL;
689 finegrid.maxTwoFXl = NULL;
690 finegrid.maxTwoFlIdx = NULL;
691 finegrid.maxTwoFXlIdx = NULL;
692
693 /* initialize ephemeris info */
694 EphemerisData *edat;
695 XLAL_CHECK( ( edat = XLALInitBarycenter( uvar_ephemEarth, uvar_ephemSun ) ) != NULL, XLAL_EFUNC );
696
697 XLALGPSSetREAL8( &minStartTimeGPS, uvar_minStartTime1 );
698 XLALGPSSetREAL8( &maxStartTimeGPS, uvar_maxStartTime1 );
699
700 /* create output files for writing if requested by user */
701 if ( uvar_printCand1 ) {
702 fnameSemiCohCand = LALCalloc( strlen( uvar_fnameout ) + 1, sizeof( CHAR ) );
703 if ( fnameSemiCohCand == NULL ) {
704 fprintf( stderr, "error allocating memory [HierarchSearchGCT.c %d]\n", __LINE__ );
705 return ( HIERARCHICALSEARCH_EMEM );
706 }
707 strcpy( fnameSemiCohCand, uvar_fnameout );
708 }
709
710 if ( uvar_printFstat1 ) {
711 const CHAR *append = "_fstatVec1.dat";
712 fnameFstatVec1 = LALCalloc( strlen( uvar_fnameout ) + strlen( append ) + 1, sizeof( CHAR ) );
713 strcpy( fnameFstatVec1, uvar_fnameout );
714 strcat( fnameFstatVec1, append );
715 if ( !( fpFstat1 = fopen( fnameFstatVec1, "wb" ) ) ) {
716 fprintf( stderr, "Unable to open Fstat file fstatvec1.out for writing.\n" );
717 return ( HIERARCHICALSEARCH_EFILE );
718 }
719 }
720
721 /*------------ Set up stacks, detector states etc. */
722 /* initialize spin range vectors */
723 XLAL_INIT_MEM( spinRange_Temp );
724
725 /* some useful first stage params */
726 usefulParams.sftbasename = uvar_DataFiles1;
727
728 /* ----- prepare generation of coherent segment list */
729 if ( XLALUserVarWasSet( &uvar_segmentList ) ) {
730 if ( XLALUserVarWasSet( &uvar_nStacksMax ) || XLALUserVarWasSet( &uvar_tStack ) ) {
731 XLALPrintError( "Use EITHER (--nStacksMax and --tStack) OR --segmentList to define the coherent segments!\n\n" );
732 return HIERARCHICALSEARCH_EBAD;
733 }
734 if ( ( usefulParams.segmentList = XLALReadSegmentsFromFile( uvar_segmentList ) ) == NULL ) {
735 XLALPrintError( "Failed to parse segment-list file '%s'. xlalErrno = %d\n", uvar_segmentList, xlalErrno );
736 return HIERARCHICALSEARCH_ESUB;
737 }
738 } else { /* set up maximally nStacksMax fixed-size segments of length tStack */
739 if ( !XLALUserVarWasSet( &uvar_tStack ) ) {
740 XLALPrintError( "Need to set --tStack or --segmentList to define the coherent segments!\n\n" );
741 return HIERARCHICALSEARCH_EBAD;
742 }
743
744 usefulParams.nStacks = uvar_nStacksMax;
745 usefulParams.tStack = uvar_tStack;
746 usefulParams.segmentList = NULL;
747 }
748 /* ----- */
749
750
751 XLAL_INIT_MEM( usefulParams.spinRange_startTime );
752 XLAL_INIT_MEM( usefulParams.spinRange_endTime );
753 XLAL_INIT_MEM( usefulParams.spinRange_refTime );
754 XLAL_INIT_MEM( usefulParams.spinRange_midTime );
755
756 /* copy user specified spin variables at reftime */
757 /* the reference time value in spinRange_refTime will be set in SetUpSFTs() */
758 usefulParams.spinRange_refTime.fkdot[0] = uvar_Freq; /* frequency */
759 usefulParams.spinRange_refTime.fkdot[1] = uvar_f1dot; /* 1st spindown */
760 usefulParams.spinRange_refTime.fkdot[2] = uvar_f2dot; /* 2nd spindown */
761 usefulParams.spinRange_refTime.fkdot[3] = uvar_f3dot; /* 3rd spindown */
762 usefulParams.spinRange_refTime.fkdotBand[0] = uvar_FreqBand; /* frequency range */
763 usefulParams.spinRange_refTime.fkdotBand[1] = uvar_f1dotBand; /* spindown range */
764 usefulParams.spinRange_refTime.fkdotBand[2] = uvar_f2dotBand; /* spindown range */
765 usefulParams.spinRange_refTime.fkdotBand[3] = uvar_f3dotBand; /* spindown range */
766
767 usefulParams.edat = edat;
768 usefulParams.minStartTimeGPS = minStartTimeGPS;
769 usefulParams.maxStartTimeGPS = maxStartTimeGPS;
770 usefulParams.blocksRngMed = uvar_blocksRngMed;
771 usefulParams.Dterms = uvar_Dterms;
772 usefulParams.DtermsRecalc = uvar_DtermsRecalc;
773 usefulParams.assumeSqrtSX = uvar_assumeSqrtSX;
774 usefulParams.SSBprec = uvar_SSBprecision;
775 usefulParams.Fmethod = uvar_FstatMethod;
776 usefulParams.FmethodRecalc = uvar_FstatMethodRecalc;
777 usefulParams.recalcToplistStats = uvar_recalcToplistStats;
778
779 usefulParams.mismatch1 = uvar_mismatch1;
780
781 /* set reference time for pulsar parameters */
782 if ( XLALUserVarWasSet( &uvar_refTime ) ) {
783 usefulParams.refTime = uvar_refTime;
784 } else {
785 XLALPrintWarning( "Reference time will be set to mid-time of observation time\n" );
786 usefulParams.refTime = -1;
787 }
788
789 /* set Fstat calculation frequency resolution (coarse grid) */
790 if ( XLALUserVarWasSet( &uvar_FreqBand ) ) {
791 if ( XLALUserVarWasSet( &uvar_dFreq ) ) {
792 usefulParams.dFreqStack = uvar_dFreq;
793 } else {
794 XLALPrintError( "--dFreq is required if --FreqBand is given\n" );
795 return ( HIERARCHICALSEARCH_EBAD );
796 }
797 } else {
798 usefulParams.dFreqStack = 0;
799 }
800
801 /* set Fstat spindown resolution (coarse grid) */
802 if ( XLALUserVarWasSet( &uvar_f1dotBand ) ) {
803 if ( XLALUserVarWasSet( &uvar_df1dot ) ) {
804 usefulParams.df1dot = uvar_df1dot;
805 } else {
806 XLALPrintError( "--df1dot is required if --f1dotBand is given\n" );
807 return ( HIERARCHICALSEARCH_EBAD );
808 }
809 } else {
810 usefulParams.df1dot = 0;
811 }
812
813 /* set Fstat 2nd spindown resolution (coarse grid) */
814 if ( XLALUserVarWasSet( &uvar_f2dotBand ) ) {
815 if ( XLALUserVarWasSet( &uvar_df2dot ) ) {
816 usefulParams.df2dot = uvar_df2dot;
817 } else {
818 XLALPrintError( "--df2dot is required if --f2dotBand is given\n" );
819 return ( HIERARCHICALSEARCH_EBAD );
820 }
821 } else {
822 usefulParams.df2dot = 0;
823 }
824
825 /* set Fstat 3rd spindown resolution (coarse grid) */
826 if ( XLALUserVarWasSet( &uvar_f3dotBand ) ) {
827 if ( XLALUserVarWasSet( &uvar_df3dot ) ) {
828 usefulParams.df3dot = uvar_df3dot;
829 } else {
830 XLALPrintError( "--df3dot is required if --f3dotBand is given\n" );
831 return ( HIERARCHICALSEARCH_EBAD );
832 }
833 } else {
834 usefulParams.df3dot = 0;
835 }
836
837 // read signal parameters to be injected, if requested by the user
838
839 if ( uvar_injectionSources != NULL ) {
840 XLAL_CHECK_MAIN( ( usefulParams.injectionSources = XLALPulsarParamsFromUserInput( uvar_injectionSources, NULL ) ) != NULL, XLAL_EFUNC );
841 }
842
843 /* for 1st stage: read sfts, calculate detector states */
844 LogPrintf( LOG_NORMAL, "Reading input data ... " );
845 LAL_CALL( SetUpSFTs( &status, &usefulParams ), &status );
846 LogPrintfVerbatim( LOG_NORMAL, " done.\n" );
847
848 /* some useful params computed by SetUpSFTs */
849 tStack = usefulParams.tStack;
850 tObs = usefulParams.tObs;
851 nStacks = usefulParams.nStacks;
852 midTstack = usefulParams.midTstack;
853 startTstack = usefulParams.startTstack;
854 endTstack = usefulParams.endTstack;
855 tMidGPS = usefulParams.spinRange_midTime.refTime;
856 refTimeGPS = usefulParams.spinRange_refTime.refTime;
857 fprintf( stderr, "%% --- GPS reference time = %.4f , GPS data mid time = %.4f\n",
858 XLALGPSGetREAL8( &refTimeGPS ), XLALGPSGetREAL8( &tMidGPS ) );
859
860 REAL4 *loudestTwoFPerSeg = NULL;
861 if ( uvar_loudestTwoFPerSeg ) {
862 loudestTwoFPerSeg = XLALCalloc( nStacks, sizeof( REAL4 ) );
863 } // if loudestTwoFPerSeg
864
865 /* free segment list */
866 if ( usefulParams.segmentList )
867 if ( XLALSegListClear( usefulParams.segmentList ) != XLAL_SUCCESS ) {
868 XLAL_ERROR( XLAL_EFUNC );
869 }
870 XLALFree( usefulParams.segmentList );
871 usefulParams.segmentList = NULL;
872
873 /*------- set frequency and spindown resolutions and ranges for Fstat and semicoherent steps -----*/
874
875 dFreqStack = usefulParams.dFreqStack;
876 df1dot = usefulParams.df1dot;
877 df2dot = usefulParams.df2dot;
878 df3dot = usefulParams.df3dot;
879 LogPrintf( LOG_NORMAL, "dFreqStack = %e, df1dot = %e, df2dot = %e, df3dot = %e\n", dFreqStack, df1dot, df2dot, df3dot );
880
881 /* set number of fine-grid spindowns */
882 if ( XLALUserVarWasSet( &uvar_gammaRefine ) ) {
883 gammaRefine = uvar_gammaRefine;
884 } else {
885 sigmasq = 0.0; /* second moment of segments' midpoints */
886 for ( k = 0; k < nStacks; k++ ) {
887 midTstackGPS = midTstack->data[k];
888 timeDiffSeg = XLALGPSDiff( &midTstackGPS, &tMidGPS );
889 sigmasq = sigmasq + ( timeDiffSeg * timeDiffSeg );
890 }
891 sigmasq = sigmasq / ( nStacks * tStack * tStack );
892 /* Refinement factor (approximate) */
893 gammaRefine = sqrt( 1.0 + 60 * sigmasq ); /* Eq. from PRL, page 3 */
894 }
895
896 /* set number of fine-grid 2nd spindowns */
897 if ( XLALUserVarWasSet( &uvar_gamma2Refine ) ) {
898 /* use 1st spindown refinement if user value < 0 */
899 if ( uvar_gamma2Refine < 0 ) {
900 gamma2Refine = gammaRefine;
901 } else {
902 gamma2Refine = uvar_gamma2Refine;
903 }
904 } else {
905 sigmasq = sigma4 = 0.0; /* second and 4th moment of segments' midpoints */
906 for ( k = 0; k < nStacks; k++ ) {
907 midTstackGPS = midTstack->data[k];
908 timeDiffSeg = XLALGPSDiff( &midTstackGPS, &tMidGPS );
909 sigmasq = sigmasq + ( timeDiffSeg * timeDiffSeg );
910 sigma4 = sigma4 + ( timeDiffSeg * timeDiffSeg * timeDiffSeg * timeDiffSeg );
911 }
912 sigmasq = sigmasq / ( nStacks * tStack * tStack );
913 sigma4 = sigma4 / ( nStacks * tStack * tStack * tStack * tStack );
914 /* 2nd spindown refinement factor gamma2/gamma1 (approximate)
915 see Pletsch, PRD 82 042002, 2010 */
916 gamma2Refine = sqrt( 2100.0 * ( sigma4 - sigmasq * sigmasq ) );
917 }
918
919 /**** debugging information ******/
920 /* print some debug info about spinrange */
921 LogPrintf( LOG_DETAIL, "Frequency and spindown range at refTime (%d): [%f,%f], [%e,%e], [%e,%e], [%e,%e]\n",
922 usefulParams.spinRange_refTime.refTime.gpsSeconds,
923 usefulParams.spinRange_refTime.fkdot[0],
924 usefulParams.spinRange_refTime.fkdot[0] + usefulParams.spinRange_refTime.fkdotBand[0],
925 usefulParams.spinRange_refTime.fkdot[1],
926 usefulParams.spinRange_refTime.fkdot[1] + usefulParams.spinRange_refTime.fkdotBand[1],
927 usefulParams.spinRange_refTime.fkdot[2],
928 usefulParams.spinRange_refTime.fkdot[2] + usefulParams.spinRange_refTime.fkdotBand[2],
929 usefulParams.spinRange_refTime.fkdot[3],
930 usefulParams.spinRange_refTime.fkdot[3] + usefulParams.spinRange_refTime.fkdotBand[3] );
931
932 LogPrintf( LOG_DETAIL, "Frequency and spindown range at startTime (%d): [%f,%f], [%e,%e], [%e,%e], [%e,%e]\n",
933 usefulParams.spinRange_startTime.refTime.gpsSeconds,
934 usefulParams.spinRange_startTime.fkdot[0],
935 usefulParams.spinRange_startTime.fkdot[0] + usefulParams.spinRange_startTime.fkdotBand[0],
936 usefulParams.spinRange_startTime.fkdot[1],
937 usefulParams.spinRange_startTime.fkdot[1] + usefulParams.spinRange_startTime.fkdotBand[1],
938 usefulParams.spinRange_startTime.fkdot[2],
939 usefulParams.spinRange_startTime.fkdot[2] + usefulParams.spinRange_startTime.fkdotBand[2],
940 usefulParams.spinRange_startTime.fkdot[3],
941 usefulParams.spinRange_startTime.fkdot[3] + usefulParams.spinRange_startTime.fkdotBand[3] );
942
943
944 LogPrintf( LOG_DETAIL, "Frequency and spindown range at midTime (%d): [%f,%f], [%e,%e], [%e,%e], [%e,%e]\n",
945 usefulParams.spinRange_midTime.refTime.gpsSeconds,
946 usefulParams.spinRange_midTime.fkdot[0],
947 usefulParams.spinRange_midTime.fkdot[0] + usefulParams.spinRange_midTime.fkdotBand[0],
948 usefulParams.spinRange_midTime.fkdot[1],
949 usefulParams.spinRange_midTime.fkdot[1] + usefulParams.spinRange_midTime.fkdotBand[1],
950 usefulParams.spinRange_midTime.fkdot[2],
951 usefulParams.spinRange_midTime.fkdot[2] + usefulParams.spinRange_midTime.fkdotBand[2],
952 usefulParams.spinRange_midTime.fkdot[3],
953 usefulParams.spinRange_midTime.fkdot[3] + usefulParams.spinRange_midTime.fkdotBand[3] );
954
955 LogPrintf( LOG_DETAIL, "Frequency and spindown range at endTime (%d): [%f,%f], [%e,%e], [%e,%e], [%e,%e]\n",
956 usefulParams.spinRange_endTime.refTime.gpsSeconds,
957 usefulParams.spinRange_endTime.fkdot[0],
958 usefulParams.spinRange_endTime.fkdot[0] + usefulParams.spinRange_endTime.fkdotBand[0],
959 usefulParams.spinRange_endTime.fkdot[1],
960 usefulParams.spinRange_endTime.fkdot[1] + usefulParams.spinRange_endTime.fkdotBand[1],
961 usefulParams.spinRange_endTime.fkdot[2],
962 usefulParams.spinRange_endTime.fkdot[2] + usefulParams.spinRange_endTime.fkdotBand[2],
963 usefulParams.spinRange_endTime.fkdot[3],
964 usefulParams.spinRange_endTime.fkdot[3] + usefulParams.spinRange_endTime.fkdotBand[3] );
965
966 /* print debug info about stacks */
967 fprintf( stderr, "%% --- Setup, N = %d, T = %.0f s, Tobs = %.0f s, gammaRefine = %.0f, gamma2Refine = %.0f, gamma3Refine = %.0f\n",
968 nStacks, tStack, tObs, gammaRefine, gamma2Refine, gamma3Refine );
969
970
971 /*---------- set up F-statistic calculation stuff ---------*/
972 /* set reference time for calculating Fstatistic */
973 thisPoint.refTime = tMidGPS; /* midpoint of data spanned */
974
975 /* binary orbit and higher spindowns not considered */
976 thisPoint.asini = 0 /* isolated pulsar */;
977 XLAL_INIT_MEM( thisPoint.fkdot );
978
979 /*---------- set up stuff for semi-coherent part ---------*/
980 /* set up some semiCoherent parameters */
981 semiCohPar.tsMid = midTstack;
982 semiCohPar.refTime = tMidGPS; // unused??
983
984 /* reference time for finegrid is midtime */
985 finegrid.refTime = tMidGPS;
986
987 /* allocate memory for pos vel acc vectors */
988 posStack = XLALCreateREAL8VectorSequence( nStacks, 3 );
989 velStack = XLALCreateREAL8VectorSequence( nStacks, 3 );
990 accStack = XLALCreateREAL8VectorSequence( nStacks, 3 );
991
992 /* calculate Earth orbital positions, velocities and accelerations */
993 LAL_CALL( GetSegsPosVelAccEarthOrb( &status, &posStack, &velStack, &accStack, &usefulParams ), &status );
994
995 /* semicoherent configuration parameters */
996 semiCohPar.pos = posStack;
997 semiCohPar.vel = velStack;
998 semiCohPar.acc = accStack;
999 semiCohPar.outBaseName = uvar_fnameout;
1000
1001 /* get number of detectors and detector name vector */
1002 LALStringVector *detectorIDs = usefulParams.detectorIDs;
1003 const UINT4 numDetectors = detectorIDs->length;
1004
1005 // compute single-IFO F-statistics for line-robust stats
1006 if ( uvar_computeBSGL ) {
1007 Fstat_what |= FSTATQ_2F_PER_DET;
1008
1009 /* take BSGL user vars and pre-compute the corresponding BSGLsetup */
1010 REAL4 *oLGX_p = NULL;
1011 REAL4 oLGX[PULSAR_MAX_DETECTORS];
1012 if ( uvar_oLGX != NULL ) {
1013 if ( uvar_oLGX->length != numDetectors ) {
1014 fprintf( stderr, "Invalid input: length(oLGX) = %d differs from number of detectors (%d)'\n", uvar_oLGX->length, numDetectors );
1015 return ( HIERARCHICALSEARCH_EBAD );
1016 }
1017 if ( XLALParseLinePriors( &oLGX[0], uvar_oLGX ) != XLAL_SUCCESS ) {
1018 fprintf( stderr, "Invalid input oLGX'\n" );
1019 return ( HIERARCHICALSEARCH_EBAD );
1020 }
1021 oLGX_p = &oLGX[0];
1022 } // if uvar_oLGX != NULL
1023
1024 usefulParams.BSGLsetup = XLALCreateBSGLSetup( numDetectors, uvar_Fstar0sc, oLGX_p, uvar_BSGLlogcorr, nStacks );
1025 if ( usefulParams.BSGLsetup == NULL ) {
1026 fprintf( stderr, "XLALCreateBSGLSetup() failed\n" );
1027 return ( HIERARCHICALSEARCH_EBAD );
1028 }
1029 } // if uvar_computeBSGL
1030
1031 /* assemble column headings string for output file */
1032 CHAR column_headings_string_base[256];
1033 if ( XLALUserVarWasSet( &uvar_f3dot ) ) {
1034 sprintf( column_headings_string_base, "freq alpha delta f1dot f2dot f3dot nc <2F>" );
1035 } else {
1036 sprintf( column_headings_string_base, "freq alpha delta f1dot f2dot nc <2F>" );
1037 }
1038
1039 UINT4 column_headings_string_length = sizeof( column_headings_string_base );
1040 if ( uvar_computeBSGL ) {
1041 column_headings_string_length += 10 + numDetectors * 8; /* 10 for " log10BSGL" and 8 per detector for " <2F_XY>" */
1042 }
1043 if ( uvar_getMaxFperSeg ) {
1044 column_headings_string_length += 12 + 13; /* 12 for " log10BSGLtL" and 13 for " log10BtSGLtL" */
1045 column_headings_string_length += 6 + 9 + numDetectors * ( 9 + 12 ); /* 6 for " max2F", 9 for " max2Fseg" and 9+12 per detector for " max2F_XY" and " max2F_XYseg"*/
1046 }
1047 if ( uvar_recalcToplistStats ) {
1048 column_headings_string_length += 6 + 11 + numDetectors * 9; /* 6 for " <2Fr>" and 9 per detector for " <2Fr_XY>" */
1049 if ( uvar_computeBSGL ) {
1050 column_headings_string_length += 11; /* for " log10BSGLr" */
1051 if ( uvar_getMaxFperSeg ) {
1052 column_headings_string_length += 13; /* for " log10BSGLtLr" */
1053 }
1054 }
1055 if ( XLALUserVarWasSet( &uvar_f3dot ) ) {
1056 column_headings_string_length += 1;
1057 }
1058 if ( uvar_loudestSegOutput ) {
1059 column_headings_string_length += 9 + numDetectors * 7; /* for " lseg 2Fl" and 7 per detector for " 2Fl_XY" */
1060 }
1061 }
1062 char column_headings_string[column_headings_string_length];
1063 XLAL_INIT_MEM( column_headings_string );
1064 strcat( column_headings_string, column_headings_string_base );
1065 if ( uvar_computeBSGL ) {
1066 strcat( column_headings_string, " log10BSGL" );
1067 for ( UINT4 X = 0; X < numDetectors ; X ++ ) {
1068 char headingX[9];
1069 snprintf( headingX, sizeof( headingX ), " <2F_%s>", detectorIDs->data[X] );
1070 strcat( column_headings_string, headingX );
1071 } /* for X < numDet */
1072 }
1073 if ( uvar_getMaxFperSeg ) {
1074 strcat( column_headings_string, " log10BSGLtL log10BtSGLtL" );
1075 strcat( column_headings_string, " max2F" );
1076 strcat( column_headings_string, " max2Fseg" );
1077 for ( UINT4 X = 0; X < numDetectors ; X ++ ) {
1078 char headingX[16];
1079 snprintf( headingX, sizeof( headingX ), " max2F_%s", detectorIDs->data[X] );
1080 strcat( column_headings_string, headingX );
1081 }
1082 for ( UINT4 X = 0; X < numDetectors ; X ++ ) {
1083 char headingX[16];
1084 snprintf( headingX, sizeof( headingX ), " max2F_%sseg", detectorIDs->data[X] );
1085 strcat( column_headings_string, headingX );
1086 } /* for X < numDet */
1087 }
1088 if ( uvar_recalcToplistStats ) {
1089 strcat( column_headings_string, " <2Fr>" );
1090 if ( uvar_computeBSGL ) {
1091 strcat( column_headings_string, " log10BSGLr" );
1092 }
1093 for ( UINT4 X = 0; X < numDetectors ; X ++ ) {
1094 char headingX[10];
1095 snprintf( headingX, sizeof( headingX ), " <2Fr_%s>", detectorIDs->data[X] );
1096 strcat( column_headings_string, headingX );
1097 } /* for X < numDet */
1098 if ( uvar_loudestSegOutput ) {
1099 strcat( column_headings_string, " lseg 2Fl" );
1100 for ( UINT4 X = 0; X < numDetectors ; X ++ ) {
1101 char headingX[10];
1102 snprintf( headingX, sizeof( headingX ), " 2Fl_%s", detectorIDs->data[X] );
1103 strcat( column_headings_string, headingX );
1104 }
1105 }
1106 if ( uvar_computeBSGL && uvar_getMaxFperSeg ) {
1107 strcat( column_headings_string, " log10BSGLtLr" );
1108 }
1109 }
1110 global_column_headings_stringp = column_headings_string;
1111
1112 /* get effective inverse number of segments per detector (needed for correct averaging in single-IFO F calculation) */
1113 REAL4 NSegmentsInv = 1.0 / nStacks; /* also need this for multi-detector F-stat averaging later on */
1114
1115 /*-----------Create template grid for first stage ---------------*/
1116 /* prepare initialization of DopplerSkyScanner to step through paramter space */
1117 scanInit.dAlpha = uvar_dAlpha;
1118 scanInit.dDelta = uvar_dDelta;
1119 scanInit.gridType = uvar_gridType1;
1120 scanInit.metricType = uvar_metricType1;
1121 scanInit.metricMismatch = uvar_mismatch1;
1122 scanInit.projectMetric = TRUE;
1123 scanInit.obsDuration = tStack;
1124 scanInit.obsBegin = tMidGPS;
1125 scanInit.ephemeris = edat;
1126 scanInit.skyGridFile = uvar_skyGridFile;
1127 scanInit.skyRegionString = ( CHAR * )LALCalloc( 1, strlen( uvar_skyRegion ) + 1 );
1128 if ( scanInit.skyRegionString == NULL ) {
1129 fprintf( stderr, "error allocating memory [HierarchSearchGCT.c.c %d]\n", __LINE__ );
1130 return ( HIERARCHICALSEARCH_EMEM );
1131 }
1132 strcpy( scanInit.skyRegionString, uvar_skyRegion );
1133
1134 // just use first SFTs' IFO for metric (should be irrelevant)
1135 const LALDetector *firstDetector = XLALGetSiteInfo( detectorIDs->data[0] );
1136 if ( firstDetector == NULL ) {
1137 LogPrintf( LOG_CRITICAL, "\nXLALGetSiteInfo() failed for detector '%s'\n", detectorIDs->data[0] );
1138 return HIERARCHICALSEARCH_EXLAL;
1139 }
1140 scanInit.Detector = firstDetector;
1141
1142 scanInit.numSkyPartitions = uvar_numSkyPartitions;
1143 scanInit.partitionIndex = uvar_partitionIndex;
1144
1145 scanInit.Freq = usefulParams.spinRange_midTime.fkdot[0] + usefulParams.spinRange_midTime.fkdotBand[0]; /* Only used if metric is employed */
1146
1147 /* initialize skygrid */
1148 LogPrintf( LOG_DETAIL, "Setting up coarse sky grid... \n" );
1149 LAL_CALL( InitDopplerSkyScan( &status, &thisScan, &scanInit ), &status );
1150
1151
1152 /* ----- start main calculations by going over coarse grid points --------*/
1153 skyGridCounter = 0;
1154 f1dotGridCounter = 0;
1155
1156 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1157
1158 /* "spool forward" if we found a checkpoint */
1159 {
1160 UINT4 count = 0; /* The first checkpoint should have value 1 */
1161 UINT4 skycount = 0;
1162
1163 if ( 0 > GET_GCT_CHECKPOINT( uvar_fnameChkPoint, semiCohToplist, semiCohToplist2, semiCohToplist3, &count ) ) {
1164 XLALPrintError( "%s : '%s' \n", HIERARCHICALSEARCH_MSGECHECKPT, uvar_fnameChkPoint );
1165 return ( HIERARCHICALSEARCH_ECHECKPT );
1166 }
1167
1168 if ( count ) {
1169 f1dotGridCounter = ( UINT4 )( count % usefulParams.nf1dot ); /* Checkpointing counter = i_sky * nf1dot + i_f1dot */
1170 skycount = ( UINT4 )( ( count - f1dotGridCounter ) / usefulParams.nf1dot );
1171 }
1172 fprintf( stderr, "%% --- Cpt:%d, total:%d, sky:%d/%d, f1dot:%d/%d\n",
1173 count, thisScan.numSkyGridPoints * usefulParams.nf1dot, skycount + 1, thisScan.numSkyGridPoints, f1dotGridCounter + 1, usefulParams.nf1dot );
1174
1175 for ( skyGridCounter = 0; skyGridCounter < skycount; skyGridCounter++ ) {
1176 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1177 }
1178
1179 if ( count == thisScan.numSkyGridPoints * usefulParams.nf1dot ) {
1180 thisScan.state = STATE_FINISHED;
1181 }
1182
1183 }
1184
1185 /* spool forward if uvar_skyPointIndex is set
1186 ---This probably doesn't make sense when checkpointing is turned on */
1187 if ( XLALUserVarWasSet( &uvar_skyPointIndex ) ) {
1188 UINT4 count = uvar_skyPointIndex;
1189 for ( skyGridCounter = 0; ( skyGridCounter < count ) && ( thisScan.state != STATE_FINISHED ) ; skyGridCounter++ ) {
1190 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1191 }
1192 }
1193
1194 /* ################## loop over SKY coarse-grid points ################## */
1195 while ( thisScan.state != STATE_FINISHED ) {
1196 /*------------- calculate F-Statistic for each segment --------------*/
1197
1198 /* normalize skyposition: correctly map into [0,2pi]x[-pi/2,pi/2] */
1199 thisPoint.Alpha = dopplerpos.Alpha;
1200 thisPoint.Delta = dopplerpos.Delta;
1201
1202 /* Calculate unit vector n */
1203 cosAlpha = cos( thisPoint.Alpha );
1204 sinAlpha = sin( thisPoint.Alpha );
1205 cosDelta = cos( thisPoint.Delta );
1206 sinDelta = sin( thisPoint.Delta );
1207 nvec[0] = cosAlpha * cosDelta;
1208 nvec[1] = sinAlpha * cosDelta;
1209 nvec[2] = sinDelta;
1210
1211 /* calculate number of bins for Fstat overhead due to residual spin-down */
1212 semiCohPar.extraBinsFstat = usefulParams.extraBinsFstat;
1213 binsFstat1 = usefulParams.binsFstatSearch + 2 * semiCohPar.extraBinsFstat;
1214
1215 /* ################## loop over coarse-grid F1DOT values ################## */
1216 ifdot = 0;
1217
1218 while ( ifdot < usefulParams.nf1dot ) {
1219
1220 /* if checkpoint read, spool forward */
1221 if ( f1dotGridCounter > 0 ) {
1222 ifdot = f1dotGridCounter;
1223 f1dotGridCounter = 0;
1224 }
1225
1226 /* ################## loop over coarse-grid F2DOT values ################## */
1227 if2dot = 0;
1228
1229 while ( if2dot < usefulParams.nf2dot ) {
1230
1231 /* ################## loop over coarse-grid F3DOT values ################## */
1232 if3dot = 0;
1233
1234 while ( if3dot < usefulParams.nf3dot ) {
1235
1236 /* show progress */
1237 LogPrintf( LOG_NORMAL, "Coarse grid sky:%d/%d f1dot:%d/%d f2dot:%d/%d f3dot:%d/%d\n",
1238 skyGridCounter + 1, thisScan.numSkyGridPoints, ifdot + 1, usefulParams.nf1dot, if2dot + 1, usefulParams.nf2dot, if3dot + 1, usefulParams.nf3dot );
1239
1240
1241 /* ------------- Set up coarse grid --------------------------------------*/
1242 coarsegrid.freqlength = ( UINT4 )( binsFstat1 );
1243 coarsegrid.nStacks = nStacks;
1244 coarsegrid.length = coarsegrid.freqlength * coarsegrid.nStacks;
1245 coarsegrid.numDetectors = numDetectors;
1246
1247 /* allocate memory for coarsegrid */
1248 coarsegrid.TwoF = ( REAL4 * )LALRealloc( coarsegrid.TwoF, coarsegrid.length * sizeof( REAL4 ) );
1249 if ( uvar_computeBSGL ) {
1250 coarsegrid.TwoFX = ( REAL4 * )LALRealloc( coarsegrid.TwoFX, coarsegrid.numDetectors * coarsegrid.length * sizeof( REAL4 ) );
1251 }
1252 coarsegrid.Uindex = ( UINT4 * )LALRealloc( coarsegrid.Uindex, coarsegrid.length * sizeof( UINT4 ) );
1253
1254 if ( coarsegrid.TwoF == NULL || coarsegrid.Uindex == NULL ) {
1255 fprintf( stderr, "ERROR: Memory allocation [HierarchSearchGCT.c %d]\n", __LINE__ );
1256 return ( HIERARCHICALSEARCH_EMEM );
1257 }
1258
1259 /* ------------- Set up fine grid --------------------------------------*/
1260
1261 /* frequency fine-grid borders */
1262 freqmin_fg = usefulParams.spinRange_midTime.fkdot[0];
1263 freqband_fg = usefulParams.spinRange_midTime.fkdotBand[0];
1264
1265 /* fine-grid frequency resolution */
1266 dfreq_fg = dFreqStack;
1267 UINT4 nfreqs_fg = ceil( freqband_fg / dfreq_fg ); /* number of points in frequency */
1268 if ( nfreqs_fg == 0 ) {
1269 nfreqs_fg++;
1270 }
1271
1272 /* copy frequency setup parameters to fine-grid struct */
1273 finegrid.freqmin_fg = freqmin_fg;
1274 finegrid.dfreq_fg = dfreq_fg;
1275 finegrid.freqlength = nfreqs_fg ;
1276#define ALIGN_REAL4 4 /* 16 bytes / sizeof(REAL4) = 4 */
1277 finegrid.freqlengthAL = ALIGN_REAL4 * ( ( UINT4 )ceil( 1.0 * finegrid.freqlength / ALIGN_REAL4 ) );
1278
1279 /* fine-grid f1dot resolution */
1280 if ( usefulParams.nf1dot == 1 ) {
1281 nf1dots_fg = 1;
1282 } else {
1283 nf1dots_fg = ceil( gammaRefine ); /* number of spindown fine-grid points */
1284 if ( ( nf1dots_fg % 2 ) == 0 ) { /* if even, add one (to refine symmetrically) */
1285 nf1dots_fg++;
1286 }
1287 }
1288 df1dot_fg = df1dot / nf1dots_fg; /* spindown fine-grid stepsize */
1289
1290 /* adjust f1dotmin_fg, so that f1dot finegrid is centered around coarse-grid f1dot point */
1291 f1dotmin_fg = ( usefulParams.spinRange_midTime.fkdot[1] + ifdot * df1dot ) - df1dot_fg * floor( nf1dots_fg / 2.0 );
1292
1293 /* fine-grid f2dot resolution */
1294 if ( uvar_f2dotBand == 0 ) {
1295 nf2dots_fg = 1;
1296 } else {
1297 nf2dots_fg = ceil( gamma2Refine ); /* number of 2nd spindown fine-grid points */
1298 if ( ( nf2dots_fg % 2 ) == 0 ) { /* if even, add one (to refine symmetrically) */
1299 nf2dots_fg++;
1300 }
1301 }
1302 df2dot_fg = df2dot / nf2dots_fg; /* 2nd spindown fine-grid stepsize */
1303
1304 /* adjust f2dotmin_fg, so that f2dot finegrid is centered around coarse-grid f2dot point */
1305 f2dotmin_fg = ( usefulParams.spinRange_midTime.fkdot[2] + if2dot * df2dot ) - df2dot_fg * floor( nf2dots_fg / 2.0 );
1306
1307 /* fine-grid f3dot resolution */
1308 nf3dots_fg = 1;
1309 df3dot_fg = 1; /* 3rd spindown fine-grid stepsize */
1310
1311 /* adjust f3dotmin_fg, so that f3dot finegrid is centered around coarse-grid f3dot point */
1312 f3dotmin_fg = ( usefulParams.spinRange_midTime.fkdot[3] + if3dot * df3dot ) - df3dot_fg * floor( nf3dots_fg / 2.0 );
1313
1314 /* total number of fine-grid points */
1315 finegrid.length = finegrid.freqlength;
1316
1317 if ( !oldcg ) {
1318 oldcg = coarsegrid.length;
1319 LogPrintfVerbatim( LOG_NORMAL, "%% --- CG:%d ", coarsegrid.length );
1320 }
1321 if ( !oldfg ) {
1322 oldfg = finegrid.length;
1323 LogPrintfVerbatim( LOG_NORMAL, "FG:%d f1dotmin_fg:%.13g df1dot_fg:%.13g f2dotmin_fg:%.13g df2dot_fg:%.13g f3dotmin_fg:%.13g df3dot_fg:%.13g\n",
1324 finegrid.length, f1dotmin_fg, df1dot_fg, f2dotmin_fg, df2dot_fg, f3dotmin_fg, df3dot_fg );
1325 }
1326 if ( ( coarsegrid.length != oldcg ) || ( finegrid.length != oldfg ) ) {
1327 LogPrintfVerbatim( LOG_CRITICAL, "ERROR: Grid-sizes disagree!\nPrevious CG:%d FG:%d, currently CG:%d FG:%d\n",
1328 oldcg, oldfg, coarsegrid.length, finegrid.length );
1329 return ( HIERARCHICALSEARCH_EVAL );
1330 }
1331
1332 /* number of detectors, needed for sumTwoFX array */
1333 finegrid.numDetectors = coarsegrid.numDetectors;
1334
1335 /* allocate memory for finegrid points */
1336 /* FIXME: The SSE2 optimized code relies on an identical alignment modulo 16 bytes of the
1337 arrays finegrid.nc and finegrid.sumTwoF !!
1338 MacOS enforces 16 byte alignment of memory blocks with size >= 16 bytes allocated
1339 with malloc and realloc, but e.g. for 32 bit Linux this will NOT hold.
1340 Alternatives might be using (posix_)memalign under Linux.
1341 Windows ==>???
1342 */
1343
1344 finegrid.nc = ( FINEGRID_NC_T * )ALRealloc( finegrid.nc, finegrid.length * sizeof( FINEGRID_NC_T ) );
1345 finegrid.sumTwoF = ( REAL4 * )ALRealloc( finegrid.sumTwoF, finegrid.length * sizeof( REAL4 ) );
1346 if ( uvar_getMaxFperSeg ) {
1347 finegrid.maxTwoFl = ( REAL4 * )ALRealloc( finegrid.maxTwoFl, finegrid.length * sizeof( REAL4 ) );
1348 finegrid.maxTwoFlIdx = ( UINT4 * )ALRealloc( finegrid.maxTwoFlIdx, finegrid.length * sizeof( UINT4 ) );
1349
1350 }
1351 if ( uvar_computeBSGL ) {
1352 finegrid.sumTwoFX = ( REAL4 * )ALRealloc( finegrid.sumTwoFX, finegrid.numDetectors * finegrid.freqlengthAL * sizeof( REAL4 ) );
1353 if ( uvar_getMaxFperSeg ) {
1354 finegrid.maxTwoFXl = ( REAL4 * )ALRealloc( finegrid.maxTwoFXl, finegrid.numDetectors * finegrid.freqlengthAL * sizeof( REAL4 ) );
1355 finegrid.maxTwoFXlIdx = ( UINT4 * )ALRealloc( finegrid.maxTwoFXlIdx, finegrid.numDetectors * finegrid.freqlengthAL * sizeof( UINT4 ) );
1356 }
1357 }
1358
1359 if ( finegrid.nc == NULL || finegrid.sumTwoF == NULL ) {
1360 fprintf( stderr, "ERROR: Memory allocation [HierarchSearchGCT.c %d]\n", __LINE__ );
1361 return ( HIERARCHICALSEARCH_EMEM );
1362 }
1363
1364 /* copy sky coarse-grid point to finegrid, because sky is not refined */
1365 finegrid.alpha = thisPoint.Alpha;
1366 finegrid.delta = thisPoint.Delta;
1367
1368 /* ---------- Walk through fine grid f1dot --------------- */
1369 for ( if1dot_fg = 0; if1dot_fg < nf1dots_fg; if1dot_fg++ ) {
1370
1371 /* get the 1st spindown of this fine-grid point */
1372 f1dot_fg = f1dotmin_fg + if1dot_fg * df1dot_fg;
1373
1374 /* ---------- Walk through fine grid f2dot --------------- */
1375 for ( if2dot_fg = 0; if2dot_fg < nf2dots_fg; if2dot_fg++ ) {
1376
1377 /* get the 2nd spindown of this fine-grid point */
1378 f2dot_fg = f2dotmin_fg + if2dot_fg * df2dot_fg;
1379
1380 /* ---------- Walk through fine grid f3dot --------------- */
1381 for ( if3dot_fg = 0; if3dot_fg < nf3dots_fg; if3dot_fg++ ) {
1382
1383 /* get the 3rd spindown of this fine-grid point */
1384 f3dot_fg = f3dotmin_fg + if3dot_fg * df3dot_fg;
1385
1386 /* initialize the entire finegrid ( 2F-sum and number count set to 0 ) */
1387 memset( finegrid.nc, 0, finegrid.length * sizeof( FINEGRID_NC_T ) );
1388 memset( finegrid.sumTwoF, 0, finegrid.length * sizeof( REAL4 ) );
1389 if ( uvar_getMaxFperSeg ) {
1390 memset( finegrid.maxTwoFl, 0, finegrid.length * sizeof( REAL4 ) );
1391 memset( finegrid.maxTwoFlIdx, 0, finegrid.length * sizeof( UINT4 ) );
1392 }
1393 if ( uvar_computeBSGL ) {
1394 memset( finegrid.sumTwoFX, 0, finegrid.numDetectors * finegrid.freqlengthAL * sizeof( REAL4 ) );
1395 if ( uvar_getMaxFperSeg ) {
1396 memset( finegrid.maxTwoFXl, 0, finegrid.numDetectors * finegrid.freqlengthAL * sizeof( REAL4 ) );
1397 memset( finegrid.maxTwoFXlIdx, 0, finegrid.numDetectors * finegrid.freqlengthAL * sizeof( UINT4 ) );
1398 }
1399 }
1400
1401 /* compute F-statistic values for coarse grid the first time through fine grid fdots loop */
1402 const BOOLEAN doComputeFstats = ( ( if1dot_fg == 0 ) && ( if2dot_fg == 0 ) && ( if3dot_fg == 0 ) );
1403
1404 /* #########################################################################*/
1405 /* ------------- MAIN LOOP over Segments for F-statistic -------------------*/
1406
1407 for ( k = 0; k < nStacks; k++ ) {
1408
1409 /* Get pos vel acc for this segment */
1410 pos[0] = semiCohPar.pos->data[3 * k];
1411 pos[1] = semiCohPar.pos->data[3 * k + 1];
1412 pos[2] = semiCohPar.pos->data[3 * k + 2];
1413
1414 vel[0] = semiCohPar.vel->data[3 * k];
1415 vel[1] = semiCohPar.vel->data[3 * k + 1];
1416 vel[2] = semiCohPar.vel->data[3 * k + 2];
1417
1418 /* currently unused:
1419 REAL8 acc[3];
1420 acc[0] = semiCohPar.acc->data[3*k];
1421 acc[1] = semiCohPar.acc->data[3*k + 1];
1422 acc[2] = semiCohPar.acc->data[3*k + 2];
1423 */
1424
1425 /* Midpoint in time of current segment */
1426 midTstackGPS = midTstack->data[k];
1427
1428 /* Difference in time between this segment's midpoint and Fstat reftime */
1429 timeDiffSeg = XLALGPSDiff( &midTstackGPS, &thisPoint.refTime );
1430
1431 /* ---------------------------------------------------------------------------------------- */
1432
1433 /* Compute sky position associated dot products
1434 (for global-correlation coordinates, from Eq. (1) in PRL) */
1435
1436 B1 = ( pos[0] * nvec[0] \
1437 + pos[1] * nvec[1] \
1438 + pos[2] * nvec[2] ); // This is \vec r \dot \vec n
1439
1440 /* current unused:
1441 REAL8 A2 = ( acc[0] * nvec[0] \
1442 + acc[1] * nvec[1] \
1443 + acc[2] * nvec[2] ); // This is \vec a \dot \vec n
1444 */
1445
1446 B2 = ( vel[0] * nvec[0] \
1447 + vel[1] * nvec[1] \
1448 + vel[2] * nvec[2] ); // This is \vec v \dot \vec n
1449
1450 A1 = 1.0 + B2;
1451
1452 /* Setup of the cell size (windows) in u1 */
1453 u1win = dFreqStack * A1;
1454 u1winInv = 1.0 / u1win; /* Precomputing */
1455
1456 /* Setup of the cell size (windows) in u2 */
1457 /* --- currently only u1 is needed.
1458 u2win = df1dot;
1459 u2winInv = 1.0/u2win; */
1460
1461 /* Set starting frequency for Fstat calculation */
1462 thisPoint.fkdot[0] = usefulParams.spinRange_midTime.fkdot[0] - semiCohPar.extraBinsFstat * dFreqStack;
1463
1464 /* Set spindown value for Fstat calculation */
1465 thisPoint.fkdot[1] = usefulParams.spinRange_midTime.fkdot[1] + ifdot * df1dot;
1466
1467 /* Set spindown value for Fstat calculation */
1468 thisPoint.fkdot[2] = usefulParams.spinRange_midTime.fkdot[2] + if2dot * df2dot;
1469
1470 /* Set spindown value for Fstat calculation */
1471 thisPoint.fkdot[3] = usefulParams.spinRange_midTime.fkdot[3] + if3dot * df3dot;
1472
1473 /* Frequency at the segment's midpoint for later use */
1474 f1dot_event = thisPoint.fkdot[1] + thisPoint.fkdot[2] * timeDiffSeg;
1475 myf0 = thisPoint.fkdot[0] + thisPoint.fkdot[1] * timeDiffSeg +
1476 + 0.5 * thisPoint.fkdot[2] * timeDiffSeg * timeDiffSeg;
1477
1478 /* Smallest values of u1 and u2 (to be subtracted) */
1479 u1start = myf0 * A1 + f1dot_event * B1; /* Eq. (1a) */
1480 U1idx = 0;
1481
1482 /* Holger: current code structure of loops (processing f1dot by f1dot) needs only U1 calculation.
1483 u2start = f1dot_event + myf0 * A2 + 2.0 * f1dot_event * B2;
1484 myf0max = myf0 + (Fstat_res->numFreqBins - 1) * dFreqStack;
1485 u2end = f1dot_event + myf0max * A2 + 2.0 * f1dot_event * B2;
1486 NumU2idx = ceil(fabs(u2start - u2end) * u2winInv);
1487 U2idx = 0;
1488 */
1489
1490 /* ----------------------------------------------------------------- */
1491 /************************ Compute F-Statistic ************************/
1492 if ( doComputeFstats ) { /* if first time through fine grid fdots loop */
1493
1494 tic_Fstat = GETTIME();
1495 const int retn = XLALComputeFstat( &Fstat_res, usefulParams.Fstat_in_vec->data[k], &thisPoint, binsFstat1, Fstat_what );
1496 if ( retn != XLAL_SUCCESS ) {
1497 XLALPrintError( "%s: XLALComputeFstat() failed with errno=%d\n", __func__, xlalErrno );
1498 return xlalErrno;
1499 }
1500
1501 /* Loop over coarse-grid frequency bins */
1502 for ( ifreq = 0; ifreq < Fstat_res->numFreqBins; ifreq++ ) {
1503
1504 /* go to next frequency coarse-grid point */
1505 freq_event = myf0 + ifreq * dFreqStack;
1506
1507 /* compute the global-correlation coordinate indices */
1508 U1idx = ComputeU1idx( freq_event, f1dot_event, A1, B1, u1start, u1winInv );
1509
1510 /* Holger: current code structure of loops (processing f1dot by f1dot) needs only U1 calculation.
1511 ComputeU2idx ( freq_event, f1dot_event, A2, B2, u2start, u2winInv, &U2idx);
1512 */
1513
1514 /* Check U1 index value */
1515 if ( ( INT4 )ifreq != U1idx ) {
1516 fprintf( stderr, "ERROR: Incorrect Frequency-Index!\n ----> Seg: %03d ifreq: %d cg U1: %d \n",
1517 k, ifreq, U1idx );
1518 return ( HIERARCHICALSEARCH_ECG );
1519 } else {
1520 /* Holger: current code structure of loops (processing f1dot by f1dot) needs only U1 calculation.
1521 coarsegrid.Uindex = U1idx * NumU2idx + U2idx; */
1522 coarsegrid.Uindex[CG_INDEX( coarsegrid, k, ifreq )] = U1idx;
1523 }
1524
1525 /* ============ Copy the *2F* value ============ */
1526 coarsegrid.TwoF[CG_INDEX( coarsegrid, k, ifreq )] = Fstat_res->twoF[ifreq];
1527 if ( uvar_computeBSGL ) {
1528 for ( UINT4 X = 0; X < coarsegrid.numDetectors; X++ ) {
1529 INT4 detid = -1;
1530 for ( UINT4 Y = 0; Y < Fstat_res->numDetectors; Y++ ) { /* look for matching detector ID in this segment */
1531 if ( strcmp( Fstat_res->detectorNames[Y], detectorIDs->data[X] ) == 0 ) {
1532 detid = Y;
1533 }
1534 }
1535 if ( detid == -1 ) { /* if no match found, detector X was not present in this segment, so use 2FX=0.0 */
1536 coarsegrid.TwoFX[CG_FX_INDEX( coarsegrid, X, k, ifreq )] = 0.0;
1537 } else { /* if a match was found, get the corresponding 2F value */
1538 coarsegrid.TwoFX[CG_FX_INDEX( coarsegrid, X, k, ifreq )] = Fstat_res->twoFPerDet[detid][ifreq];
1539 }
1540 } /* for X < numDetectors */
1541 } /* if ( uvar_computeBSGL ) */
1542
1543 } /* END: Loop over coarse-grid frequency bins (ifreq) */
1544
1545 /* print fstat vector if required -- mostly for debugging */
1546 if ( uvar_printFstat1 ) {
1547 LAL_CALL( PrintFstatVec( &status, Fstat_res, fpFstat1, &thisPoint, refTimeGPS, k + 1 ), &status );
1548 }
1549
1550 // if requested, keep track of loudest candidate from each segment
1551 if ( loudestTwoFPerSeg ) {
1552 REAL4 loudestFstat_k = findLoudestTwoF( Fstat_res );
1553 loudestTwoFPerSeg[k] = fmaxf( loudestTwoFPerSeg[k], loudestFstat_k );
1554 }
1555 /* --- Holger: This is not needed in U1-only case. Sort the coarse grid in Uindex --- */
1556 /* qsort(coarsegrid.list, (size_t)coarsegrid.length, sizeof(CoarseGridPoint), compareCoarseGridUindex); */
1557
1558 time_Fstat += ( GETTIME() - tic_Fstat );
1559
1560 } // if (doComputeFstats)
1561 /* -------------------- END Compute F-Statistic -------------------- */
1562
1563 /* -------------------- Map fine grid to coarse grid -------------------- */
1564 tic_SumFine = GETTIME();
1565
1566 /* get the frequency of this fine-grid point at mid point of segment */
1567 /* OLD: ifreq_fg = 0; freq_tmp = finegrid.freqmin_fg + ifreq_fg * finegrid.dfreq_fg + f1dot_tmp * timeDiffSeg; */
1568 f1dot_event_fg = f1dot_fg + f2dot_fg * timeDiffSeg;
1569 freq_fg = finegrid.freqmin_fg + f1dot_fg * timeDiffSeg +
1570 0.5 * f2dot_fg * timeDiffSeg * timeDiffSeg; /* first fine-grid frequency */
1571
1572 /* compute the global-correlation coordinate indices */
1573 U1idx = ComputeU1idx( freq_fg, f1dot_event_fg, A1, B1, u1start, u1winInv );
1574
1575 if ( U1idx < 0 ) {
1576 fprintf( stderr, "ERROR: Stepped outside the coarse grid (%d)! \n", U1idx );
1577 return ( HIERARCHICALSEARCH_ECG );
1578 }
1579
1580 if ( U1idx + finegrid.freqlength - 1 >= Fstat_res->numFreqBins ) {
1581 fprintf( stderr, "ERROR: Stepped outside the coarse grid (%d:%d:%d:%d)! \n",
1582 U1idx, finegrid.freqlength, U1idx + finegrid.freqlength - 1, Fstat_res->numFreqBins );
1583 return ( HIERARCHICALSEARCH_ECG );
1584 }
1585
1586 /* coarse grid over frequency for this stack */
1587 REAL4 *cgrid2F = coarsegrid.TwoF + CG_INDEX( coarsegrid, k, U1idx );
1588
1589 /* fine grid over frequency */
1590 REAL4 *fgrid2F = finegrid.sumTwoF + FG_INDEX( finegrid, 0 );
1591#ifndef EXP_NO_NUM_COUNT
1592 FINEGRID_NC_T *fgridnc = finegrid.nc + FG_INDEX( finegrid, 0 );
1593#endif
1594
1595#ifdef GC_SSE2_OPT
1596 if ( uvar_getMaxFperSeg ) {
1597 /* disables number count keeping */
1598 REAL4 *fgrid2Fmax = finegrid.maxTwoFl + FG_INDEX( finegrid, 0 );
1599 UINT4 *fgrid2FmaxIdx = finegrid.maxTwoFlIdx + FG_INDEX( finegrid, 0 );
1600
1601 gc_hotloop_2Fmax_tracking( fgrid2F, fgrid2Fmax, fgrid2FmaxIdx, cgrid2F, k, finegrid.freqlength );
1602 } else {
1603#ifndef EXP_NO_NUM_COUNT
1604 gc_hotloop( fgrid2F, cgrid2F, fgridnc, TwoFthreshold, finegrid.freqlength );
1605#else
1606 gc_hotloop_no_nc( fgrid2F, cgrid2F, finegrid.freqlength );
1607#endif
1608 }
1609 if ( uvar_computeBSGL ) {
1610 for ( UINT4 X = 0; X < finegrid.numDetectors; X++ ) {
1611 REAL4 *cgrid2FX = coarsegrid.TwoFX + CG_FX_INDEX( coarsegrid, X, k, U1idx );
1612 REAL4 *fgrid2FX = finegrid.sumTwoFX + FG_FX_INDEX( finegrid, X, 0 );
1613
1614 if ( uvar_getMaxFperSeg ) {
1615 REAL4 *fgrid2FXmax = finegrid.maxTwoFXl + FG_FX_INDEX( finegrid, X, 0 );
1616 UINT4 *fgrid2FXmaxIdx = finegrid.maxTwoFXlIdx + FG_FX_INDEX( finegrid, X, 0 );
1617
1618 gc_hotloop_2Fmax_tracking( fgrid2FX, fgrid2FXmax, fgrid2FXmaxIdx, cgrid2FX, k, finegrid.freqlength );
1619 } else {
1620 gc_hotloop_no_nc( fgrid2FX, cgrid2FX, finegrid.freqlength );
1621 }
1622 } /* for X */
1623 }
1624#else // GC_SSE2_OPT
1625 for ( UINT4 ifreq_fg = 0; ifreq_fg < finegrid.freqlength; ifreq_fg++ ) {
1626 fgrid2F[0] += cgrid2F[0];
1627#ifndef EXP_NO_NUM_COUNT
1628 fgridnc[0] += ( TwoFthreshold < cgrid2F[0] );
1629 fgridnc++;
1630#endif // EXP_NO_NUM_COUNT
1631 fgrid2F++;
1632 cgrid2F++;
1633 }
1634 if ( uvar_computeBSGL ) {
1635 for ( UINT4 X = 0; X < finegrid.numDetectors; X++ ) {
1636 REAL4 *cgrid2FX = coarsegrid.TwoFX + CG_FX_INDEX( coarsegrid, X, k, U1idx );
1637 REAL4 *fgrid2FX = finegrid.sumTwoFX + FG_FX_INDEX( finegrid, X, 0 );
1638 for ( UINT4 ifreq_fg = 0; ifreq_fg < finegrid.freqlength; ifreq_fg++ ) {
1639 fgrid2FX[0] += cgrid2FX[0];
1640 fgrid2FX++;
1641 cgrid2FX++;
1642 }
1643 }
1644 }
1645
1646
1647 if ( uvar_getMaxFperSeg ) {
1648 cgrid2F = coarsegrid.TwoF + CG_INDEX( coarsegrid, k, U1idx );
1649 REAL4 *fgridMax2Fl = finegrid.maxTwoFl + FG_INDEX( finegrid, 0 );
1650 UINT4 *fgrid2FmaxIdx = finegrid.maxTwoFlIdx + FG_INDEX( finegrid, 0 );
1651 int isLouder;
1652 for ( UINT4 ifreq_fg = 0; ifreq_fg < finegrid.freqlength; ifreq_fg++ ) {
1653 isLouder = ( fgridMax2Fl[0] <= cgrid2F[0] );
1654 fgridMax2Fl[0] = fmaxf( fgridMax2Fl[0], cgrid2F[0] );
1655 fgrid2FmaxIdx[0] = isLouder * k + ( 1 - isLouder ) * fgrid2FmaxIdx[0];
1656 fgridMax2Fl++;
1657 fgrid2FmaxIdx++;
1658 cgrid2F++;
1659 }
1660 for ( UINT4 X = 0; X < finegrid.numDetectors; X++ ) {
1661 REAL4 *cgrid2FX = coarsegrid.TwoFX + CG_FX_INDEX( coarsegrid, X, k, U1idx );
1662 REAL4 *fgridMax2FXl = finegrid.maxTwoFXl + FG_FX_INDEX( finegrid, X, 0 );
1663 UINT4 *fgrid2FXmaxIdx = finegrid.maxTwoFXlIdx + FG_FX_INDEX( finegrid, X, 0 );
1664
1665 for ( UINT4 ifreq_fg = 0; ifreq_fg < finegrid.freqlength; ifreq_fg++ ) {
1666 isLouder = ( fgridMax2FXl[0] <= cgrid2FX[0] );
1667 fgridMax2FXl[0] = fmaxf( fgridMax2FXl[0], cgrid2FX[0] );
1668 fgrid2FXmaxIdx[0] = isLouder * k + ( 1 - isLouder ) * fgrid2FXmaxIdx[0];
1669 fgrid2FXmaxIdx++;
1670 fgridMax2FXl++;
1671 cgrid2FX++;
1672 }
1673 }
1674 }
1675#endif // GC_SSE2_OPT
1676 time_SumFine += ( GETTIME() - tic_SumFine );
1677
1678 } /* end: ------------- MAIN LOOP over Segments --------------------*/
1679
1680 /* ############################################################### */
1681
1682 tic_ExtraStats = GETTIME();
1683
1684 if ( uvar_semiCohToplist ) {
1685 /* this is necessary here, because UpdateSemiCohToplists() might set
1686 a checkpoint that needs some information from here */
1687
1688 if ( ( uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL || uvar_SortToplist == SORTBY_F_BSGLtL_BtSGLtL ) &&
1689 finegrid.sumTwoFX && finegrid.maxTwoFXl ) {
1690 LAL_CALL( UpdateSemiCohToplistsOptimTriple( &status, uvar_SortToplist, semiCohToplist, semiCohToplist2, semiCohToplist3, &finegrid, f1dot_fg, f2dot_fg, f3dot_fg, &usefulParams, NSegmentsInv, usefulParams.NSegmentsInvX, XLALUserVarWasSet( &uvar_f3dot ) ), &status );
1691 } else {
1692 LAL_CALL( UpdateSemiCohToplists( &status, semiCohToplist, semiCohToplist2, semiCohToplist3, &finegrid, f1dot_fg, f2dot_fg, f3dot_fg, &usefulParams, NSegmentsInv, usefulParams.NSegmentsInvX, XLALUserVarWasSet( &uvar_f3dot ) ), &status );
1693 }
1694 } // if semiCohToplist
1695 time_ExtraStats += ( GETTIME() - tic_ExtraStats );
1696
1697 } /* for( if1dot_fg = 0; if1dot_fg < nf1dots_fg; if1dot_fg++ ) */
1698 } /* for( if2dot_fg = 0; if2dot_fg < nf2dots_fg; if2dot_fg++ ) */
1699 } /* for( if3dot_fg = 0; if3dot_fg < nf3dots_fg; if3dot_fg++ ) */
1700 /* ---------- END walk through fine grid fdots --------------- */
1701 if3dot++; /* Increment if3dot counter */
1702
1703 } /* ########## End of loop over coarse-grid f3dot values (if3dot) ########## */
1704 if2dot++; /* Increment if2dot counter */
1705
1706 } /* ########## End of loop over coarse-grid f2dot values (if2dot) ########## */
1707 ifdot++; /* Increment ifdot counter BEFORE SET_GCT_CHECKPOINT */
1708
1709 if ( !uvar_outputTiming ) {
1710 SET_GCT_CHECKPOINT( uvar_fnameChkPoint, semiCohToplist, semiCohToplist2, semiCohToplist3, skyGridCounter * usefulParams.nf1dot + ifdot, TRUE );
1711 }
1712
1713 } /* ########## End of loop over coarse-grid f1dot values (ifdot) ########## */
1714
1715 /* continue forward till the end if uvar_skyPointIndex is set
1716 ---This probably doesn't make sense when checkpointing is turned on */
1717 if ( XLALUserVarWasSet( &uvar_skyPointIndex ) ) {
1718 while ( thisScan.state != STATE_FINISHED ) {
1719 skyGridCounter++;
1720 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1721 }
1722 } else {
1723 skyGridCounter++;
1724
1725 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1726 }
1727
1728 } /* ######## End of while loop over 1st stage SKY coarse-grid points ############ */
1729 /*---------------------------------------------------------------------------------*/
1730
1731 /* now that we have the final toplist, translate all pulsar parameters to correct reftime */
1732 xlalErrno = 0;
1733 XLALExtrapolateToplistPulsarSpins( semiCohToplist, usefulParams.spinRange_refTime.refTime, finegrid.refTime );
1734 if ( semiCohToplist2 ) { // handle (optional) second toplist
1735 XLALExtrapolateToplistPulsarSpins( semiCohToplist2, usefulParams.spinRange_refTime.refTime, finegrid.refTime );
1736 }
1737 if ( semiCohToplist3 ) { // handle (optional) second toplist
1738 XLALExtrapolateToplistPulsarSpins( semiCohToplist3, usefulParams.spinRange_refTime.refTime, finegrid.refTime );
1739 }
1740 if ( xlalErrno != 0 ) {
1741 XLALPrintError( "%s line %d : XLALExtrapolateToplistPulsarSpins() failed with xlalErrno = %d.\n\n", __func__, __LINE__, xlalErrno );
1742 return ( HIERARCHICALSEARCH_EXLAL );
1743 }
1744 LogPrintf( LOG_NORMAL, "Finished main analysis.\n" );
1745
1746
1747 // output averaged F-stat timing model (--outputTimingDetails)
1748 if ( uvar_outputTimingDetails != NULL ) {
1749 FILE *fp;
1750 XLAL_CHECK( ( fp = fopen( uvar_outputTimingDetails, "ab" ) ) != NULL, XLAL_ESYS, "Failed to open '%s' for writing\n", uvar_outputTimingDetails );
1751 for ( size_t l = 0; l < nStacks; l++ ) {
1752 XLAL_CHECK( XLALAppendFstatTiming2File( usefulParams.Fstat_in_vec->data[l], fp, ( l == 0 ) ) == XLAL_SUCCESS, XLAL_EFUNC );
1753 }
1754 fclose( fp );
1755 }
1756
1757 /* Also compute F, FX (for line-robust statistics) for all candidates in final toplist */
1758 tic_RecalcToplist = GETTIME();
1759 if ( uvar_recalcToplistStats ) {
1760
1761 LogPrintf( LOG_NORMAL, "Recalculating statistics for the final toplist...\n" );
1762 RecalcStatsParams XLAL_INIT_DECL( recalcParams );
1763 recalcParams.listEntryTypeName = "GCTtop";
1764 if ( usefulParams.Fstat_in_vec_recalc != NULL ) {
1765 recalcParams.Fstat_in_vec = usefulParams.Fstat_in_vec_recalc;
1766 } else {
1767 recalcParams.Fstat_in_vec = usefulParams.Fstat_in_vec;
1768 }
1769 timing.RecalcMethodStr = XLALGetFstatInputMethodName( recalcParams.Fstat_in_vec->data[0] );
1770 recalcParams.detectorIDs = usefulParams.detectorIDs;
1771 recalcParams.startTstack = usefulParams.startTstack;
1772 recalcParams.refTimeGPS = refTimeGPS;
1773 recalcParams.BSGLsetup = usefulParams.BSGLsetup;
1774 recalcParams.loudestSegOutput = uvar_loudestSegOutput;
1775 recalcParams.computeBSGLtL = uvar_getMaxFperSeg;
1776 XLAL_CHECK( XLAL_SUCCESS == XLALComputeExtraStatsForToplist( semiCohToplist, &recalcParams ),
1777 HIERARCHICALSEARCH_EXLAL, "XLALComputeExtraStatsForToplist() failed with xlalErrno = %d.\n\n", xlalErrno
1778 );
1779 // also recalc optional 2nd toplist if present
1780 if ( semiCohToplist2 ) {
1781 XLAL_CHECK( XLAL_SUCCESS == XLALComputeExtraStatsForToplist( semiCohToplist2, &recalcParams ),
1782 HIERARCHICALSEARCH_EXLAL, "XLALComputeExtraStatsForToplist() failed for 2nd toplist with xlalErrno = %d.\n\n", xlalErrno
1783 );
1784 }
1785 // also recalc optional 3rd toplist if present
1786 if ( semiCohToplist3 ) {
1787 XLAL_CHECK( XLAL_SUCCESS == XLALComputeExtraStatsForToplist( semiCohToplist3, &recalcParams ),
1788 HIERARCHICALSEARCH_EXLAL, "XLALComputeExtraStatsForToplist() failed for 3rd toplist with xlalErrno = %d.\n\n", xlalErrno
1789 );
1790 }
1791
1792 LogPrintf( LOG_NORMAL, "Finished recalculating toplist statistics.\n" );
1793 } // if recalcToplist
1794 time_RecalcToplist = GETTIME() - tic_RecalcToplist;
1795
1796 if ( uvar_outputTiming ) {
1797 timing.FstatMethodStr = XLALGetFstatInputMethodName( usefulParams.Fstat_in_vec->data[0] );
1798 timing.RecalcMethodStr = ( timing.RecalcMethodStr == NULL ) ? "NONE" : timing.RecalcMethodStr;
1799 timing.Nseg = coarsegrid.nStacks;
1800 timing.Ndet = coarsegrid.numDetectors;
1801 timing.Tcoh = usefulParams.tStack;
1802 timing.Nsft = usefulParams.nSFTs;
1803 timing.Ncand = uvar_nCand1;
1804
1805 timing.NFreqCo = coarsegrid.freqlength; // includes Fstat sideband bins
1806 timing.Ncoh = thisScan.numSkyGridPoints * timing.NFreqCo * usefulParams.nf1dot * usefulParams.nf2dot;
1807
1808 REAL8 nf1dot_fine = usefulParams.nf1dot * nf1dots_fg; // 'nf1dots_fg' is the number of fine-grid points *per coarse-grid point*!
1809 REAL8 nf2dot_fine = usefulParams.nf2dot * nf2dots_fg; // 'nf1dots_fg' is the number of fine-grid points *per coarse-grid point*!
1810 timing.Ninc = thisScan.numSkyGridPoints * usefulParams.binsFstatSearch * nf1dot_fine * nf2dot_fine; // excludes F-stat sideband bins
1811
1812 timing.tau_Fstat = time_Fstat / ( timing.Nseg * timing.Ncoh * timing.Ndet );
1813 timing.tau_SumF = time_SumFine / ( timing.Nseg * timing.Ninc );
1814 timing.tau_Bayes = time_ExtraStats / timing.Ninc;
1815 timing.tau_Recalc = time_RecalcToplist / timing.Ncand;
1816 time_Total = GETTIME() - tic_Start;
1817 timing.time_Other = time_Total - ( time_Fstat + time_SumFine + time_ExtraStats + time_RecalcToplist ); // what's left
1818
1819 if ( uvar_outputTiming ) {
1820 XLAL_CHECK( write_TimingInfo( uvar_outputTiming, &timing ) == XLAL_SUCCESS, XLAL_EFUNC );
1821 }
1822 } // if uvar_outputTiming
1823
1824 // output loudest per-segment candidates
1825 if ( loudestTwoFPerSeg ) {
1826 CHAR *fname = XLALStringDuplicate( uvar_fnameout );
1827 fname = XLALStringAppend( fname, "_loudestTwoFPerSeg" );
1828 FILE *fp;
1829 if ( ( fp = fopen( fname, "wb" ) ) == NULL ) {
1830 XLALPrintError( "Unable to open '%s' for writing\n", fname );
1831 return ( HIERARCHICALSEARCH_EFILE );
1832 }
1833 for ( UINT4 l = 0; l < nStacks; l ++ ) {
1834 fprintf( fp, "%.6" LAL_REAL4_FORMAT "\n", loudestTwoFPerSeg[l] );
1835 }
1836 fclose( fp );
1837 XLALFree( fname );
1838 } // if loudestTwoFPerSeg
1839
1840 char t1_suffix[16];
1841 char t2_suffix[16];
1842 char t3_suffix[16];
1843
1844 XLAL_INIT_MEM( t1_suffix );
1845 XLAL_INIT_MEM( t2_suffix );
1846 XLAL_INIT_MEM( t3_suffix );
1847
1848 /* additional toplist(s) suffix selection */
1849
1850 switch ( uvar_SortToplist ) {
1851 case SORTBY_DUAL_F_BSGL : {
1852 strcpy( t1_suffix, "" );
1853 strcpy( t2_suffix, "-BSGL" );
1854 strcpy( t3_suffix, "" );
1855 break;
1856 };
1857 case SORTBY_TRIPLE_BStSGLtL:
1858 case SORTBY_F_BSGLtL_BtSGLtL : {
1859 strcpy( t1_suffix, "" );
1860 strcpy( t2_suffix, "-BSGLtL" );
1861 strcpy( t3_suffix, "-BtSGLtL" );
1862 break;
1863 };
1864 default : {
1865 strcpy( t1_suffix, "" );
1866 strcpy( t2_suffix, "" );
1867 strcpy( t3_suffix, "" );
1868 };
1869 }
1870
1871 LogPrintf( LOG_DEBUG, "Writing output ... " );
1872 {
1873 UINT4 newlen = strlen( uvar_fnameout ) + strlen( t1_suffix ) + 1;;
1874 CHAR *fname1;
1875 XLAL_CHECK( ( fname1 = XLALCalloc( 1, newlen ) ) != NULL, XLAL_ENOMEM, "Failed to XLALCalloc(1, %d)\n\n", newlen );
1876 sprintf( fname1, "%s%s", uvar_fnameout, t1_suffix );
1877 XLAL_CHECK( write_hfs_oputput( fname1, semiCohToplist ) != -1, XLAL_EFAILED, "write_hfs_oputput('%s', toplist) failed for 1st toplist!\n", fname1 );
1878 XLALFree( fname1 );
1879 }
1880
1881
1882
1883 /* output optional additional toplists if any */
1884
1885 if ( semiCohToplist2 ) {
1886
1887 LogPrintfVerbatim( LOG_DEBUG, "toplist2 ... " );
1888 UINT4 newlen = strlen( uvar_fnameout ) + strlen( t2_suffix ) + 1;
1889 CHAR *fname2;
1890 XLAL_CHECK( ( fname2 = XLALCalloc( 1, newlen ) ) != NULL, XLAL_ENOMEM, "Failed to XLALCalloc(1, %d)\n\n", newlen );
1891 sprintf( fname2, "%s%s", uvar_fnameout, t2_suffix );
1892 XLAL_CHECK( write_hfs_oputput( fname2, semiCohToplist2 ) != -1, XLAL_EFAILED, "write_hfs_oputput('%s', toplist2) failed for 2nd toplist!\n", fname2 );
1893 XLALFree( fname2 );
1894
1895 if ( semiCohToplist3 ) {
1896 LogPrintfVerbatim( LOG_DEBUG, "toplist3 ... " );
1897 newlen = strlen( uvar_fnameout ) + strlen( t3_suffix ) + 1;
1898 CHAR *fname3;
1899 XLAL_CHECK( ( fname3 = XLALCalloc( 1, newlen ) ) != NULL, XLAL_ENOMEM, "Failed to XLALCalloc(1, %d)\n\n", newlen );
1900 sprintf( fname3, "%s%s", uvar_fnameout, t3_suffix );
1901 XLAL_CHECK( write_hfs_oputput( fname3, semiCohToplist3 ) != -1, XLAL_EFAILED, "write_hfs_oputput('%s', toplist3) failed for 3rd toplist!\n", fname3 );
1902 XLALFree( fname3 );
1903 }
1904 }
1905 LogPrintfVerbatim( LOG_DEBUG, "done.\n" );
1906
1907 clear_gct_checkpoint( uvar_fnameChkPoint );
1908
1909 /*------------ free all remaining memory -----------*/
1910
1911 if ( uvar_printCand1 ) {
1912 LALFree( fnameSemiCohCand );
1913 }
1914
1915 if ( uvar_printFstat1 ) {
1916 fclose( fpFstat1 );
1917 LALFree( fnameFstatVec1 );
1918 }
1919
1920 if ( usefulParams.injectionSources ) {
1921 XLALDestroyPulsarParamsVector( usefulParams.injectionSources );
1922 }
1923
1924 XLALDestroyFstatInputVector( usefulParams.Fstat_in_vec );
1925 XLALDestroyFstatInputVector( usefulParams.Fstat_in_vec_recalc );
1926
1927 XLALDestroyFstatResults( Fstat_res );
1928
1929 XLALDestroyTimestampVector( startTstack );
1930 XLALDestroyTimestampVector( midTstack );
1931 XLALDestroyTimestampVector( endTstack );
1932
1933 /* free Vel/Pos/Acc vectors and ephemeris */
1934 XLALDestroyREAL8VectorSequence( posStack );
1935 XLALDestroyREAL8VectorSequence( velStack );
1936 XLALDestroyREAL8VectorSequence( accStack );
1937 XLALDestroyEphemerisData( edat );
1938
1939 /* free dopplerscan stuff */
1940 LAL_CALL( FreeDopplerSkyScan( &status, &thisScan ), &status );
1941 if ( scanInit.skyRegionString ) {
1942 LALFree( scanInit.skyRegionString );
1943 }
1944
1945 XLALDestroyStringVector( detectorIDs );
1946
1947 /* free fine grid and coarse grid */
1948 if ( finegrid.nc ) {
1949 ALFree( finegrid.nc );
1950 }
1951 if ( finegrid.sumTwoF ) {
1952 ALFree( finegrid.sumTwoF );
1953 }
1954 if ( finegrid.sumTwoFX ) {
1955 ALFree( finegrid.sumTwoFX );
1956 }
1957 if ( finegrid.maxTwoFl ) {
1958 ALFree( finegrid.maxTwoFl );
1959 }
1960 if ( finegrid.maxTwoFlIdx ) {
1961 ALFree( finegrid.maxTwoFlIdx );
1962 }
1963
1964 if ( finegrid.maxTwoFXl ) {
1965 ALFree( finegrid.maxTwoFXl );
1966 }
1967
1968 if ( finegrid.maxTwoFXlIdx ) {
1969 ALFree( finegrid.maxTwoFXlIdx );
1970 }
1971
1972 if ( coarsegrid.TwoF ) {
1973 LALFree( coarsegrid.TwoF );
1974 }
1975 if ( coarsegrid.TwoFX ) {
1976 LALFree( coarsegrid.TwoFX );
1977 }
1978 if ( coarsegrid.Uindex ) {
1979 LALFree( coarsegrid.Uindex );
1980 }
1981
1982 free_gctFstat_toplist( &semiCohToplist );
1983 if ( semiCohToplist2 ) {
1984 free_gctFstat_toplist( &semiCohToplist2 );
1985 }
1986 if ( semiCohToplist3 ) {
1987 free_gctFstat_toplist( &semiCohToplist3 );
1988 }
1989
1990 XLALDestroyBSGLSetup( usefulParams.BSGLsetup );
1991
1992 XLALDestroyUserVars();
1993
1994 XLALFree( VCSInfoString );
1995
1996 XLALFree( loudestTwoFPerSeg );
1997
1998 LALCheckMemoryLeaks();
1999
2000 return HIERARCHICALSEARCH_ENORM;
2001} /* main */
2002
2003
2004
2005
2006
2007
2008
2009/**
2010 * Set up stacks, read SFTs, calculate SFT noise weights and calculate
2011 * detector-state
2012 */
2013void SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */
2014 UsefulStageVariables *in /**< input params */
2015 )
2016{
2017
2018 SFTCatalog *catalog = NULL;
2019 static SFTConstraints constraints;
2020 REAL8 timebase, tObs, deltaFsft;
2021 UINT4 k, numSFT;
2022 LIGOTimeGPS tStartGPS, tEndGPS, refTimeGPS, tMidGPS, midTstackGPS, startTstackGPS, endTstackGPS;
2023 SFTCatalogSequence catalogSeq;
2024 REAL8 midTseg, startTseg, endTseg;
2025
2026 BOOLEAN crc_check;
2027
2028 INITSTATUS( status );
2029 ATTATCHSTATUSPTR( status );
2030
2031 /* get sft catalog */
2032 constraints.minStartTime = &( in->minStartTimeGPS );
2033 constraints.maxStartTime = &( in->maxStartTimeGPS );
2034 XLAL_CHECK_LAL( status, ( catalog = XLALSFTdataFind( in->sftbasename, &constraints ) ) != NULL, XLAL_EFUNC );
2035
2036 /* check CRC sums of SFTs */
2037 XLAL_CHECK_LAL( status, XLALCheckCRCSFTCatalog( &crc_check, catalog ) == XLAL_SUCCESS, XLAL_EFUNC );
2038 if ( !crc_check ) {
2039 LogPrintf( LOG_CRITICAL, "SFT validity check failed\n" );
2040 ABORT( status, HIERARCHICALSEARCH_ESFT, HIERARCHICALSEARCH_MSGESFT );
2041 }
2042
2043 /* set some sft parameters */
2044 deltaFsft = catalog->data[0].header.deltaF;
2045 timebase = 1.0 / deltaFsft;
2046
2047 /* get sft catalogs for each stack */
2048 if ( in->segmentList ) { /* if segment list was given by user */
2049 SFTCatalogSequence *catalogSeq_p;
2050 if ( ( catalogSeq_p = XLALSetUpStacksFromSegmentList( catalog, in->segmentList ) ) == NULL ) {
2051 XLALPrintError( "%s: XLALSetUpStacksFromSegmentList() failed to set up segments from given list.\n", __func__ );
2052 ABORT( status, HIERARCHICALSEARCH_ESUB, HIERARCHICALSEARCH_MSGESUB );
2053 }
2054 catalogSeq = ( *catalogSeq_p ); /* copy top-level struct */
2055 XLALFree( catalogSeq_p ); /* free alloc'ed top-level struct after copying (contents survive in catalogSeq!) */
2056
2057 /* we need to set tStack here:
2058 * this will be used for setting Freq,f1dot resolution on segments, therefore we use the longest segment duration
2059 */
2060 UINT4 iSeg;
2061 REAL8 maxT = 0;
2062 for ( iSeg = 0; iSeg < in->segmentList->length; iSeg++ ) {
2063 REAL8 T = XLALGPSDiff( &( in->segmentList->segs[iSeg].end ), &( in->segmentList->segs[iSeg].start ) );
2064 maxT = HSMAX( maxT, T );
2065 }
2066 in->tStack = maxT;
2067 } else { /* set up nStacks segments of fixed span tStack */
2068 TRY( SetUpStacks( status->statusPtr, &catalogSeq, in->tStack, catalog, in->nStacks ), status );
2069 }
2070
2071 /* reset number of stacks */
2072 UINT4 numSegments = catalogSeq.length;
2073 in->nStacks = numSegments;
2074
2075 /* calculate start and end times and tobs from segmented catalog*/
2076 tStartGPS = catalogSeq.data[0].data[0].header.epoch;
2077 in->tStartGPS = tStartGPS;
2078 SFTCatalog *LastSegmentCat = &( catalogSeq.data[numSegments - 1] );
2079 UINT4 numSFTsInLastSeg = LastSegmentCat->length;
2080 tEndGPS = LastSegmentCat->data[numSFTsInLastSeg - 1].header.epoch;
2081 XLALGPSAdd( &tEndGPS, timebase );
2082 tObs = XLALGPSDiff( &tEndGPS, &tStartGPS );
2083 in->tObs = tObs;
2084
2085 /* get timestamps of start, mid and end times of each stack */
2086 /* set up vector containing mid times of stacks */
2087 in->midTstack = XLALCreateTimestampVector( in->nStacks );
2088
2089 /* set up vector containing start times of stacks */
2090 in->startTstack = XLALCreateTimestampVector( in->nStacks );
2091
2092 /* set up vector containing end times of stacks */
2093 in->endTstack = XLALCreateTimestampVector( in->nStacks );
2094
2095 /* now loop over stacks and get time stamps */
2096 for ( k = 0; k < in->nStacks; k++ ) {
2097
2098 if ( catalogSeq.data[k].length == 0 ) {
2099 /* something is wrong */
2100 ABORT( status, HIERARCHICALSEARCH_EVAL, HIERARCHICALSEARCH_MSGEVAL );
2101 }
2102
2103 /* start time of stack = time of first sft in stack */
2104 in->startTstack->data[k] = catalogSeq.data[k].data[0].header.epoch;
2105
2106 /* end time of stack = time of last sft in stack */
2107 numSFT = catalogSeq.data[k].length;
2108 in->endTstack->data[k] = catalogSeq.data[k].data[numSFT - 1].header.epoch;
2109
2110 /* reference time for Fstat has to be midpoint of segment */
2111 startTstackGPS = in->startTstack->data[k];
2112 endTstackGPS = in->endTstack->data[k];
2113
2114 startTseg = XLALGPSGetREAL8( &startTstackGPS );
2115 endTseg = XLALGPSGetREAL8( &endTstackGPS );
2116 /*
2117 TRY ( LALGPStoFloat( status->statusPtr, &startTseg, &startTstackGPS ), status);
2118 TRY ( LALGPStoFloat( status->statusPtr, &endTseg, &endTstackGPS ), status);
2119 */
2120 midTseg = startTseg + ( ( endTseg - startTseg + timebase ) * 0.5 );
2121
2122 XLALGPSSetREAL8( &midTstackGPS, midTseg );
2123 /*
2124 TRY ( LALFloatToGPS( status->statusPtr, &midTstackGPS, &midTseg), status);
2125 */
2126 in->midTstack->data[k] = midTstackGPS;
2127
2128 } /* loop over k */
2129
2130
2131 /* set reference time for pulsar parameters */
2132 /* first calculate the mid time of observation time span*/
2133 {
2134 REAL8 tStart8, tEnd8, tMid8;
2135
2136 tStart8 = XLALGPSGetREAL8( &tStartGPS );
2137 tEnd8 = XLALGPSGetREAL8( &tEndGPS );
2138 tMid8 = 0.5 * ( tStart8 + tEnd8 );
2139 XLALGPSSetREAL8( &tMidGPS, tMid8 );
2140 }
2141
2142 if ( in->refTime > 0 ) {
2143 REAL8 refTime = in->refTime;
2144 XLALGPSSetREAL8( &refTimeGPS, refTime );
2145 } else { /* set refTime to exact midtime of the total observation-time spanned */
2146 refTimeGPS = tMidGPS;
2147 }
2148
2149 /* get frequency and fdot bands at start time of sfts by extrapolating from reftime */
2150 in->spinRange_refTime.refTime = refTimeGPS;
2151 XLAL_CHECK_LAL( status, XLALExtrapolatePulsarSpinRange( &in->spinRange_startTime, &in->spinRange_refTime, XLALGPSDiff( &tStartGPS, &( &in->spinRange_refTime )->refTime ) ) == XLAL_SUCCESS, XLAL_EFUNC );
2152 XLAL_CHECK_LAL( status, XLALExtrapolatePulsarSpinRange( &in->spinRange_endTime, &in->spinRange_refTime, XLALGPSDiff( &tEndGPS, &( &in->spinRange_refTime )->refTime ) ) == XLAL_SUCCESS, XLAL_EFUNC );
2153 XLAL_CHECK_LAL( status, XLALExtrapolatePulsarSpinRange( &in->spinRange_midTime, &in->spinRange_refTime, XLALGPSDiff( &tMidGPS, &( &in->spinRange_refTime )->refTime ) ) == XLAL_SUCCESS, XLAL_EFUNC );
2154
2155 /* set Fstat spindown resolution (coarse grid) */
2156 in->df1dot = HSMIN( in->df1dot, in->spinRange_midTime.fkdotBand[1] );
2157
2158 /* calculate number of bins for Fstat overhead due to residual spin-down */
2159 in->extraBinsFstat = ( UINT4 )( 0.25 * ( in->tObs * in->df1dot + in->tObs * in->tObs * in->df2dot ) / in->dFreqStack + 1e-6 ) + 1;
2160
2161 /* calculate total number of bins for Fstat */
2162 if ( in->dFreqStack == 0 ) {
2163 in->binsFstatSearch = 1;
2164 } else {
2165 in->binsFstatSearch = ( UINT4 )( in->spinRange_midTime.fkdotBand[0] / in->dFreqStack + 1e-6 ) + 1;
2166 }
2167 /* number of coarse grid spindown values */
2168 if ( in->df1dot == 0 ) {
2169 in->nf1dot = 1;
2170 } else {
2171 in->nf1dot = ( UINT4 ) ceil( in->spinRange_midTime.fkdotBand[1] / in->df1dot ) + 1;
2172 }
2173 /* number of coarse grid 2nd spindown values */
2174 if ( in->df2dot == 0 ) {
2175 in->nf2dot = 1;
2176 } else {
2177 in->nf2dot = ( UINT4 ) floor( in->spinRange_midTime.fkdotBand[2] / in->df2dot + NUDGE ) + 1;
2178 }
2179 /* number of coarse grid 3rd spindown values */
2180 if ( in->df3dot == 0 ) {
2181 in->nf3dot = 1;
2182 } else {
2183 in->nf3dot = ( UINT4 ) floor( in->spinRange_midTime.fkdotBand[3] / in->df3dot + NUDGE ) + 1;
2184 }
2185
2186 /* set wings of sfts to be read */
2187 REAL8 minCoverFreq, maxCoverFreq;
2188 REAL8 asiniMax = 0, PeriodMin = 0, maxEcc = 0;
2189 // NOTE: *must* use spin-range at *mid-time* (not reftime), which is where the GCT code sets up its
2190 // template bank. This is potentially 'wider' than the physically-requested template bank, and
2191 // can therefore also require more SFT frequency bins!
2192 // NOTE2: a second trap here is that the GCT code does not strictly respect the given bands, but can exceed them on the upside
2193 // due to 'conversative' discretization on bins, seen above. (binsFstatSearch,nf1dot,nf2dot,nf3dot)
2194 // therefore we use this 'effective' spinrange in order to be able to estimate the required number of SFT bins correctly:
2195 PulsarSpinRange spinRangeEff = in->spinRange_midTime;
2196 spinRangeEff.fkdotBand[0] = ( in->binsFstatSearch - 1 ) * in->dFreqStack;
2197 spinRangeEff.fkdotBand[1] = ( in->nf1dot - 1 ) * in->df1dot;
2198 spinRangeEff.fkdotBand[2] = ( in->nf2dot - 1 ) * in->df2dot;
2199 spinRangeEff.fkdotBand[3] = ( in->nf3dot - 1 ) * in->df3dot;
2200 XLALCWSignalCoveringBand( &minCoverFreq, &maxCoverFreq, &tStartGPS, &tEndGPS, &( spinRangeEff ), asiniMax, PeriodMin, maxEcc );
2201
2202 REAL8 freqmin = minCoverFreq - in->extraBinsFstat * in->dFreqStack;
2203 REAL8 freqmax = maxCoverFreq + in->extraBinsFstat * in->dFreqStack;
2204
2205 /* fill detector name vector with all detectors present in any data sements */
2206 if ( ( in->detectorIDs = XLALListIFOsInCatalog( catalog ) ) == NULL ) {
2207 ABORT( status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2208 }
2209 const UINT4 numDetectors = in->detectorIDs->length;
2210
2211 /* set up vector of Fstat input data structs */
2212 in->Fstat_in_vec = XLALCreateFstatInputVector( in->nStacks );
2213 if ( in->Fstat_in_vec == NULL ) {
2214 ABORT( status, HIERARCHICALSEARCH_EMEM, HIERARCHICALSEARCH_MSGEMEM );
2215 }
2216 // if using different Fstat-method for main search and recalc: set-up separate Fstat input
2217 if ( in->recalcToplistStats && ( in->Fmethod != in->FmethodRecalc ) ) {
2218 in->Fstat_in_vec_recalc = XLALCreateFstatInputVector( in->nStacks );
2219 if ( in->Fstat_in_vec == NULL ) {
2220 ABORT( status, HIERARCHICALSEARCH_EMEM, HIERARCHICALSEARCH_MSGEMEM );
2221 }
2222 }
2223
2224 in->nSFTs = 0;
2225 for ( UINT4 X = 0; X < numDetectors; X++ ) {
2226 in->NSegmentsInvX[X] = 0;
2227 }
2228
2229 FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
2230 optionalArgs.SSBprec = in->SSBprec;
2231 optionalArgs.Dterms = in->Dterms;
2232 optionalArgs.runningMedianWindow = in->blocksRngMed;
2233 optionalArgs.FstatMethod = in->Fmethod;
2234 optionalArgs.collectTiming = in->collectFstatTiming;
2235 optionalArgs.injectSources = in->injectionSources;
2236
2237 FstatOptionalArgs XLAL_INIT_DECL( optionalArgsRecalc );
2238
2239 /* loop over segments and read sfts */
2240 for ( k = 0; k < in->nStacks; k++ ) {
2241
2242 /* if flag is given, assume a PSD with sqrt(S) = 1.0 */
2243 MultiNoiseFloor s_assumeSqrtSX;
2244 if ( in->assumeSqrtSX != NULL ) {
2245 const SFTCatalog *catalog_k = &( catalogSeq.data[k] );
2246 LALStringVector *detectorIDs_k = NULL;
2247 if ( ( detectorIDs_k = XLALListIFOsInCatalog( catalog_k ) ) == NULL ) {
2248 ABORT( status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2249 }
2250 if ( XLALParseMultiNoiseFloorMapped( &s_assumeSqrtSX, detectorIDs_k, in->assumeSqrtSX, in->detectorIDs ) != XLAL_SUCCESS ) {
2251 XLALPrintError( "%s: XLALParseMultiNoiseFloorMapped() failed with errno=%d", __func__, xlalErrno );
2252 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2253 }
2254 optionalArgs.assumeSqrtSX = &s_assumeSqrtSX;
2255 XLALDestroyStringVector( detectorIDs_k );
2256 } else {
2257 optionalArgs.assumeSqrtSX = NULL;
2258 }
2259
2260 /* ----- create Fstat input data struct ----- */
2261 if ( k == 0 ) {
2262 optionalArgs.prevInput = NULL;
2263 } else {
2264 optionalArgs.prevInput = in->Fstat_in_vec->data[0]; // re-use shared workspace from first segment for all subsequent segments
2265 }
2266
2267 in->Fstat_in_vec->data[k] = XLALCreateFstatInput( &catalogSeq.data[k], freqmin, freqmax, in->dFreqStack, in->edat, &optionalArgs );
2268 if ( in->Fstat_in_vec->data[k] == NULL ) {
2269 XLALPrintError( "%s: XLALCreateFstatInput() failed with errno=%d", __func__, xlalErrno );
2270 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2271 }
2272 // ----- if recalc uses a different Fstat-method from main search, we'll setup its own Fstat setup struct
2273 if ( in->recalcToplistStats && ( in->Fstat_in_vec_recalc != NULL ) ) {
2274 optionalArgsRecalc = optionalArgs;
2275 optionalArgsRecalc.FstatMethod = in->FmethodRecalc;
2276 optionalArgsRecalc.Dterms = in->DtermsRecalc;
2277 if ( k == 0 ) {
2278 optionalArgsRecalc.prevInput = NULL;
2279 } else {
2280 optionalArgsRecalc.prevInput = in->Fstat_in_vec_recalc->data[0]; // re-use shared workspace from first segment for all subsequent segments
2281 }
2282
2283 in->Fstat_in_vec_recalc->data[k] = XLALCreateFstatInput( &catalogSeq.data[k], freqmin, freqmax, in->dFreqStack, in->edat, &optionalArgsRecalc );
2284 if ( in->Fstat_in_vec->data[k] == NULL ) {
2285 XLALPrintError( "%s: XLALCreateFstatInput() failed with errno=%d", __func__, xlalErrno );
2286 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2287 }
2288 } // if Fstat_in_vec_recalc
2289
2290 if ( k == 0 ) {
2291 LogPrintf( LOG_NORMAL, "Search FstatMethod used: '%s'\n", XLALGetFstatInputMethodName( in->Fstat_in_vec->data[0] ) );
2292 LogPrintf( LOG_NORMAL, "Recalc FstatMethod used: '%s'\n", XLALGetFstatInputMethodName( in->Fstat_in_vec_recalc ? in->Fstat_in_vec_recalc->data[0] : in->Fstat_in_vec->data[0] ) );
2293 }
2294
2295 // --------------------------------------------------
2296 /* get SFT detectors and timestamps */
2297 const MultiLALDetector *multiIFO = XLALGetFstatInputDetectors( in->Fstat_in_vec->data[k] );
2298 if ( multiIFO == NULL ) {
2299 XLALPrintError( "%s: XLALGetFstatInputDetectors() failed with errno=%d", __func__, xlalErrno );
2300 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2301 }
2302 const MultiLIGOTimeGPSVector *multiTS = XLALGetFstatInputTimestamps( in->Fstat_in_vec->data[k] );
2303 if ( multiTS == NULL ) {
2304 XLALPrintError( "%s: XLALGetFstatInputTimestamps() failed with errno=%d", __func__, xlalErrno );
2305 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2306 }
2307
2308 /* ----- get effective inverse number of segments per detector (needed for correct averaging in single-IFO F calculation) ----- */
2309 for ( UINT4 X = 0; X < numDetectors; X++ ) {
2310 /* for each detector, check if present in each segment, and save the number of segments where it is */
2311 for ( UINT4 Y = 0; Y < multiTS->length; Y++ ) {
2312 if ( strcmp( multiIFO->sites[Y].frDetector.prefix, in->detectorIDs->data[X] ) == 0 ) {
2313 in->NSegmentsInvX[X] += 1;
2314 }
2315 } /* for Y < numDetectors */
2316 } /* for X < numDetectors */
2317
2318 /* ----- print debug info about SFTs in this stack ----- */
2319 LogPrintf( LOG_DETAIL, "Segment %d ", k + 1 );
2320 for ( UINT4 j = 0; j < multiIFO->length; j++ ) {
2321 LogPrintfVerbatim( LOG_DETAIL, "%s: %d ", multiIFO->sites[j].frDetector.prefix, multiTS->data[j]->length );
2322 }
2323 LogPrintfVerbatim( LOG_DETAIL, "\n" );
2324
2325 /* ----- count the total and per-segment number of SFTs used ----- */
2326 UINT4 nSFTsInSeg = 0;
2327 for ( UINT4 X = 0; X < multiTS->length; ++X ) {
2328 nSFTsInSeg += multiTS->data[X]->length;
2329 }
2330 in->nSFTs += nSFTsInSeg;
2331
2332 /* ----- if we have a segment-list: double-check number of SFTs ----- */
2333 if ( in->segmentList ) {
2334 /* check the number of SFTs we found in this segment against the nominal value, stored in the segment list field 'id' */
2335 UINT4 nSFTsExpected = in->segmentList->segs[k].id;
2336 if ( ( nSFTsExpected > 0 ) && ( nSFTsInSeg != nSFTsExpected ) ) {
2337 XLALPrintError( "%s: Segment list seems inconsistent with data read: segment %d contains %d SFTs, should hold %d SFTs\n", __func__, k, nSFTsInSeg, nSFTsExpected );
2338 ABORT( status, HIERARCHICALSEARCH_EBAD, HIERARCHICALSEARCH_MSGEBAD );
2339 }
2340 } /* if have segmentList */
2341
2342 } /* loop over k */
2343 for ( UINT4 X = 0; X < numDetectors; X++ ) {
2344 in->NSegmentsInvX[X] = 1.0 / in->NSegmentsInvX[X]; /* now it is the inverse number */
2345 }
2346 LogPrintf( LOG_NORMAL, "Number of segments: %d, total number of SFTs in segments: %d\n", in->nStacks, in->nSFTs );
2347
2348 /* realloc if nStacks != in->nStacks */
2349 /* if ( in->nStacks > nStacks ) { */
2350
2351 /* in->midTstack->length = nStacks; */
2352 /* in->midTstack->data = (LIGOTimeGPS *)LALRealloc( in->midTstack->data, nStacks * sizeof(LIGOTimeGPS)); */
2353
2354 /* in->startTstack->length = nStacks; */
2355 /* in->startTstack->data = (LIGOTimeGPS *)LALRealloc( in->startTstack->data, nStacks * sizeof(LIGOTimeGPS)); */
2356
2357 /* stackMultiSFT->length = nStacks; */
2358 /* stackMultiSFT->data = (MultiSFTVector **)LALRealloc( stackMultiSFT->data, nStacks * sizeof(MultiSFTVector *)); */
2359
2360 /* } */
2361
2362 /* we don't need the original catalog anymore*/
2363 XLALDestroySFTCatalog( catalog );
2364
2365 /* free catalog sequence */
2366 for ( k = 0; k < in->nStacks; k++ ) {
2367 if ( catalogSeq.data[k].length > 0 ) {
2368 LALFree( catalogSeq.data[k].data );
2369 } /* end if */
2370 } /* loop over stacks */
2371 LALFree( catalogSeq.data );
2372
2373 DETATCHSTATUSPTR( status );
2374 RETURN( status );
2375
2376} /* SetUpSFTs */
2377
2378
2379
2380
2381
2382
2383/**
2384 * \brief Breaks up input sft catalog into specified number of stacks
2385 * Loops over elements of the catalog, assigns a bin index and
2386 * allocates memory to the output catalog sequence appropriately. If
2387 * there are long gaps in the data, then some of the catalogs in the
2388 * output catalog sequence may be of zero length.
2389 */
2390void SetUpStacks( LALStatus *status, /**< pointer to LALStatus structure */
2391 SFTCatalogSequence *out, /**< Output catalog of sfts -- one for each stack */
2392 REAL8 tStack, /**< Output duration of each stack */
2393 SFTCatalog *in, /**< Input sft catalog to be broken up into stacks (ordered in increasing time)*/
2394 UINT4 nStacksMax ) /**< User specified number of stacks */
2395{
2396 UINT4 j, stackCounter, length;
2397 REAL8 tStart, thisTime;
2398 REAL8 Tsft;
2399
2400 INITSTATUS( status );
2401 ATTATCHSTATUSPTR( status );
2402
2403 /* check input parameters */
2404 ASSERT( in != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2405 ASSERT( in->length > 0, status, HIERARCHICALSEARCH_EVAL, HIERARCHICALSEARCH_MSGEVAL );
2406 ASSERT( nStacksMax > 0, status, HIERARCHICALSEARCH_EVAL, HIERARCHICALSEARCH_MSGEVAL );
2407 ASSERT( in != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2408 ASSERT( tStack > 0, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2409 ASSERT( out != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2410
2411 /* set memory of output catalog sequence to maximum possible length */
2412 out->length = nStacksMax;
2413 out->data = ( SFTCatalog * )LALCalloc( 1, nStacksMax * sizeof( SFTCatalog ) );
2414 if ( out->data == NULL ) {
2415 ABORT( status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2416 }
2417
2418
2419 Tsft = 1.0 / in->data[0].header.deltaF;
2420
2421 /* get first sft timestamp */
2422 /* tStart will be start time of a given stack.
2423 This initializes tStart to the first sft time stamp as this will
2424 be the start time of the first stack */
2425 tStart = XLALGPSGetREAL8( &( in->data[0].header.epoch ) );
2426
2427 /* loop over the sfts */
2428 stackCounter = 0;
2429 for ( j = 0; j < in->length; j++ ) {
2430 /* thisTime is current sft timestamp */
2431 thisTime = XLALGPSGetREAL8( &( in->data[j].header.epoch ) );
2432
2433 /* if sft lies in stack duration then add
2434 this sft to the stack. Otherwise move
2435 on to the next stack */
2436 if ( ( thisTime - tStart + Tsft <= tStack ) ) {
2437 out->data[stackCounter].length += 1;
2438
2439 length = out->data[stackCounter].length;
2440
2441 /* realloc to increase length of catalog */
2442 out->data[stackCounter].data = ( SFTDescriptor * )LALRealloc( out->data[stackCounter].data, length * sizeof( SFTDescriptor ) );
2443 if ( out->data[stackCounter].data == NULL ) {
2444 ABORT( status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2445 }
2446
2447 out->data[stackCounter].data[length - 1] = in->data[j];
2448 } else { /* move onto the next stack */
2449 if ( stackCounter + 1 == nStacksMax ) {
2450 break;
2451 }
2452
2453 stackCounter++;
2454
2455 /* reset start time of stack */
2456 tStart = XLALGPSGetREAL8( &( in->data[j].header.epoch ) );
2457
2458 /* realloc to increase length of catalog and copy data */
2459 out->data[stackCounter].length = 1; /* first entry in new stack */
2460 out->data[stackCounter].data = ( SFTDescriptor * )LALRealloc( out->data[stackCounter].data, sizeof( SFTDescriptor ) );
2461 if ( out->data[stackCounter].data == NULL ) {
2462 ABORT( status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2463 }
2464
2465 out->data[stackCounter].data[0] = in->data[j];
2466 } /* if new stack */
2467
2468 } /* loop over sfts */
2469
2470 /* realloc catalog sequence length to actual number of stacks */
2471 out->length = stackCounter + 1;
2472 out->data = ( SFTCatalog * )LALRealloc( out->data, ( stackCounter + 1 ) * sizeof( SFTCatalog ) );
2473 if ( out->data == NULL ) {
2474 ABORT( status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2475 }
2476
2477 DETATCHSTATUSPTR( status );
2478 RETURN( status );
2479
2480} /* SetUpStacks() */
2481
2482
2483
2484
2485
2486
2487
2488/** Print some sft catalog info */
2489void PrintCatalogInfo( LALStatus *status,
2490 const SFTCatalog *catalog,
2491 FILE *fp )
2492{
2493
2494 INT4 nSFT;
2495 LIGOTimeGPS start, end;
2496
2497 INITSTATUS( status );
2498 ATTATCHSTATUSPTR( status );
2499
2500 ASSERT( fp != NULL, status, HIERARCHICALSEARCH_EFILE, HIERARCHICALSEARCH_MSGEFILE );
2501 ASSERT( catalog != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2502
2503 nSFT = catalog->length;
2504 start = catalog->data[0].header.epoch;
2505 end = catalog->data[nSFT - 1].header.epoch;
2506
2507 fprintf( fp, "## Number of SFTs: %d\n", nSFT );
2508 fprintf( fp, "## First SFT timestamp: %d %d\n", start.gpsSeconds, start.gpsNanoSeconds );
2509 fprintf( fp, "## Last SFT timestamp: %d %d\n", end.gpsSeconds, end.gpsNanoSeconds );
2510
2511 DETATCHSTATUSPTR( status );
2512 RETURN( status );
2513
2514}
2515
2516
2517
2518
2519
2520
2521/** Print some stack info from sft catalog sequence*/
2522void PrintStackInfo( LALStatus *status,
2523 const SFTCatalogSequence *catalogSeq,
2524 FILE *fp )
2525{
2526
2527 INT4 nStacks, k;
2528
2529 INITSTATUS( status );
2530 ATTATCHSTATUSPTR( status );
2531
2532 ASSERT( fp != NULL, status, HIERARCHICALSEARCH_EFILE, HIERARCHICALSEARCH_MSGEFILE );
2533 ASSERT( catalogSeq != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2534 ASSERT( catalogSeq->length > 0, status, HIERARCHICALSEARCH_EVAL, HIERARCHICALSEARCH_MSGEVAL );
2535 ASSERT( catalogSeq->data != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2536
2537 nStacks = catalogSeq->length;
2538 fprintf( fp, "## Number of stacks: %d\n", nStacks );
2539
2540 for ( k = 0; k < nStacks; k++ ) {
2541 fprintf( fp, "## Stack No. %d : \n", k + 1 );
2542 TRY( PrintCatalogInfo( status->statusPtr, catalogSeq->data + k, fp ), status );
2543 }
2544
2545 fprintf( fp, "\n\n" );
2546
2547 DETATCHSTATUSPTR( status );
2548 RETURN( status );
2549
2550}
2551
2552#ifdef __GNUC__
2553#define likely(x) __builtin_expect(!!(x), 1)
2554#define unlikely(x) __builtin_expect(!!(x), 0)
2555#else
2556#define likely(x) (x)
2557#define unlikely(x) (x)
2558#endif
2559
2560/**
2561 * Get SemiCoh candidates into toplist(s)
2562 * This function allows for inserting candidates into up to 3 toplists at once, which might be sorted differently!
2563 */
2564void UpdateSemiCohToplistsOptimTriple( LALStatus *status,
2565 SortBy_t toplists_sortby,
2566 toplist_t *list1,
2567 toplist_t *list2, //< optional (can be NULL): insert candidate into this 2nd toplist as well
2568 toplist_t *list3, //< optional (can be NULL): insert candidate into this 3rd toplist as well
2569 FineGrid *in,
2570 REAL8 f1dot_fg,
2571 REAL8 f2dot_fg,
2572 REAL8 f3dot_fg,
2573 UsefulStageVariables *usefulparams,
2574 REAL4 NSegmentsInv,
2575 REAL4 *NSegmentsInvX,
2576 BOOLEAN have_f3dot
2577 )
2578{
2579
2580 REAL8 freq_fg;
2581 UINT4 ifreq_fg;
2582 GCTtopOutputEntry line;
2583 BOOLEAN delay_compute_BSGL = ( toplists_sortby == SORTBY_F_BSGLtL_BtSGLtL );
2584
2585 INITSTATUS( status );
2586 ATTATCHSTATUSPTR( status );
2587
2588 ASSERT( list1 != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2589 ASSERT( in != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2590 ASSERT( usefulparams != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2591
2592 /* Optimized version for triple toplist cases: BSGL, BSGLtL,BtSGLtL or
2593 2F, BSGLtL,BtSGLtL
2594
2595 First compute all detection metrics by which any of the toplists is sorted (even as seconary
2596 comparison criterion).
2597
2598 Next test if the candidate is loud enough to make it into at least one of the toplists.
2599 If so, build the candidate structure.
2600 If not, just try the next fine grid entry.
2601
2602 */
2603
2604
2605 /* ---------- Walk through fine-grid and insert candidates into toplist--------------- */
2606 for ( ifreq_fg = 0; ifreq_fg < in->freqlength; ifreq_fg++ ) {
2607
2608 freq_fg = in->freqmin_fg + ifreq_fg * in->dfreq_fg;
2609
2610
2611 /* local placeholders for summed 2F value over segments, not averages yet */
2612 REAL4 sumTwoF = in->sumTwoF[ifreq_fg];
2613 REAL4 sumTwoFX[PULSAR_MAX_DETECTORS];
2614
2615 /* compute BSGL */
2616
2617 line.maxTwoFl = in->maxTwoFl[ifreq_fg];
2618 for ( UINT4 X = 0; X < in->numDetectors; X++ ) {
2619 int fg_FX_idx = FG_FX_INDEX( *in, X, ifreq_fg );
2620 sumTwoFX[X] = in->sumTwoFX[fg_FX_idx]; /* here it's still the summed 2F value over segments, not the average */
2621 line.maxTwoFXl[X] = in->maxTwoFXl[fg_FX_idx];
2622 }
2623 xlalErrno = 0;
2624
2625 /* Note: sorting a toplist by 2F takes BSGL as a secondary sorting criterion, so even if the 1st toplist is
2626 sorted by 2F, we still provide an overestimate of BSGL for a test if it can make it into the toplist in
2627 the first place */
2628
2629 if ( delay_compute_BSGL ) {
2630 line.log10BSGL = LAL_REAL4_MAX;
2631 } else {
2632 line.log10BSGL = XLALComputeBSGL( sumTwoF, sumTwoFX, usefulparams->BSGLsetup );
2633 }
2634 if ( xlalErrno != 0 ) {
2635 XLALPrintError( "%s line %d : XLALComputeBSGL() failed with xlalErrno = %d.\n\n", __func__, __LINE__, xlalErrno );
2636 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2637 }
2638 if ( unlikely( line.log10BSGL < -LAL_REAL4_MAX * 0.1 ) ) {
2639 line.log10BSGL = -LAL_REAL4_MAX * 0.1; /* avoid minimum value, needed for output checking in print_gctFstatline_to_str() */
2640 }
2641
2642 line.log10BSGLtL = XLALComputeBSGLtL( sumTwoF, sumTwoFX, line.maxTwoFXl, usefulparams->BSGLsetup );
2643 if ( unlikely( xlalErrno != 0 ) ) {
2644 XLALPrintError( "%s line %d : XLALComputeBSGLtL() failed with xlalErrno = %d.\n\n", __func__, __LINE__, xlalErrno );
2645 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2646 }
2647 if ( unlikely( line.log10BSGLtL < -LAL_REAL4_MAX * 0.1 ) ) {
2648 line.log10BSGLtL = -LAL_REAL4_MAX * 0.1; /* avoid minimum value, needed for output checking in print_gctFstatline_to_str() */
2649 }
2650
2651 line.log10BtSGLtL = XLALComputeBtSGLtL( line.maxTwoFl, sumTwoFX, line.maxTwoFXl, usefulparams->BSGLsetup );
2652 if ( unlikely( xlalErrno != 0 ) ) {
2653 XLALPrintError( "%s line %d : XLALComputeBSGLtL() failed with xlalErrno = %d.\n\n", __func__, __LINE__, xlalErrno );
2654 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2655 }
2656 if ( unlikely( line.log10BtSGLtL < -LAL_REAL4_MAX * 0.1 ) ) {
2657 line.log10BtSGLtL = -LAL_REAL4_MAX * 0.1; /* avoid minimum value, needed for output checking in print_gctFstatline_to_str() */
2658 }
2659
2660 line.Freq = freq_fg; /* NOTE: this is not the final output frequency! For performance reasons, it will only later get correctly extrapolated for the final toplist */
2661 line.Alpha = in->alpha;
2662 line.Delta = in->delta;
2663 line.F1dot = f1dot_fg;
2664 line.F2dot = f2dot_fg;
2665 line.F3dot = f3dot_fg;
2666 line.have_f3dot = have_f3dot;
2667 line.nc = in->nc[ifreq_fg];
2668
2669 /* take F-stat averages over segments */
2670 line.avTwoF = sumTwoF * NSegmentsInv; /* average multi-2F by full number of segments */
2671
2672 /* now test if this candidate makes it into any of the toplists, if not we don't need
2673 to do any more copying and computing of data from finegrid to toplist entry structure */
2674
2675
2676 int isIncludedToplists = TEST_FSTAT_TOPLIST_INCLUSION( list1, &line ) ||
2677 TEST_FSTAT_TOPLIST_INCLUSION( list2, &line ) ||
2678 TEST_FSTAT_TOPLIST_INCLUSION( list3, &line );
2679
2680
2681 if ( likely( ! isIncludedToplists ) ) {
2682 continue;
2683 }
2684
2685
2686 /* if we delayed the computation of BSGL for a 2F sorted first toplist, now is the time to actually compute it */
2687 if ( delay_compute_BSGL ) {
2688 line.log10BSGL = XLALComputeBSGL( sumTwoF, sumTwoFX, usefulparams->BSGLsetup );
2689 }
2690
2691 line.numDetectors = in->numDetectors;
2692 line.avTwoFrecalc = -1.0; /* initialise this to -1.0, so that it only gets written out by print_gctFstatline_to_str if later overwritten in recalcToplistStats step */
2693 line.log10BSGLrecalc = -LAL_REAL4_MAX; /* for now, block field with minimal value, needed for output checking in print_gctFstatline_to_str() */
2694 line.log10BSGLtLrecalc = -LAL_REAL4_MAX; /* for now, block field with minimal value, needed for output checking in print_gctFstatline_to_str() */
2695 line.loudestSeg = -1;
2696 line.twoFloudestSeg = -1.0;
2697 for ( UINT4 X = 0; X < PULSAR_MAX_DETECTORS; X++ ) { /* initialise single-IFO F-stat arrays to zero */
2698 line.twoFXloudestSeg[X] = -1.0;
2699 line.avTwoFXrecalc[X] = 0.0;
2700 }
2701
2702
2703 line.maxTwoFlSeg = in->maxTwoFlIdx[ifreq_fg];
2704 for ( UINT4 X = 0; X < in->numDetectors; X++ ) {
2705 line.avTwoFX[X] = sumTwoFX[X] * NSegmentsInvX[X]; /* average single-2F by per-IFO number of segments */
2706 line.maxTwoFXlSeg[X] = in->maxTwoFXlIdx[FG_FX_INDEX( *in, X, ifreq_fg )];
2707 }
2708
2709 insert_into_gctFstat_toplist( list1, &line );
2710 insert_into_gctFstat_toplist( list2, &line );
2711 insert_into_gctFstat_toplist( list3, &line );
2712
2713 } // for ifreq_fg
2714
2715 DETATCHSTATUSPTR( status );
2716 RETURN( status );
2717
2718} /* UpdateSemiCohToplistsOptimTriple() */
2719
2720
2721
2722
2723/**
2724 * Get SemiCoh candidates into toplist(s)
2725 * This function allows for inserting candidates into up to 3 toplists at once, which might be sorted differently!
2726 */
2727void UpdateSemiCohToplists( LALStatus *status,
2728 toplist_t *list1,
2729 toplist_t *list2, //< optional (can be NULL): insert candidate into this 2nd toplist as well
2730 toplist_t *list3, //< optional (can be NULL): insert candidate into this 3rd toplist as well
2731 FineGrid *in,
2732 REAL8 f1dot_fg,
2733 REAL8 f2dot_fg,
2734 REAL8 f3dot_fg,
2735 UsefulStageVariables *usefulparams,
2736 REAL4 NSegmentsInv,
2737 REAL4 *NSegmentsInvX,
2738 BOOLEAN have_f3dot
2739 )
2740{
2741
2742 REAL8 freq_fg;
2743 UINT4 ifreq_fg;
2744 GCTtopOutputEntry line;
2745
2746 INITSTATUS( status );
2747 ATTATCHSTATUSPTR( status );
2748
2749 ASSERT( list1 != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2750 ASSERT( in != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2751 ASSERT( usefulparams != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2752
2753 /* ---------- Walk through fine-grid and insert candidates into toplist--------------- */
2754 for ( ifreq_fg = 0; ifreq_fg < in->freqlength; ifreq_fg++ ) {
2755
2756 freq_fg = in->freqmin_fg + ifreq_fg * in->dfreq_fg;
2757
2758 line.Freq = freq_fg; /* NOTE: this is not the final output frequency! For performance reasons, it will only later get correctly extrapolated for the final toplist */
2759 line.Alpha = in->alpha;
2760 line.Delta = in->delta;
2761 line.F1dot = f1dot_fg;
2762 line.F2dot = f2dot_fg;
2763 line.F3dot = f3dot_fg;
2764 line.nc = in->nc[ifreq_fg];
2765 line.avTwoF = 0.0; /* will be set to average over segments later */
2766 line.maxTwoFl = -1.0; /* initialise this to -1.0, so that it only gets written out by print_gctFstatline_to_str if actually computed */
2767 line.maxTwoFlSeg = -1;
2768 line.log10BSGL = -LAL_REAL4_MAX; /* for now, block field with minimal value, needed for output checking in print_gctFstatline_to_str() */
2769 line.log10BSGLtL = -LAL_REAL4_MAX;
2770 line.log10BtSGLtL = -LAL_REAL4_MAX;
2771
2772 line.numDetectors = in->numDetectors;
2773 for ( UINT4 X = 0; X < PULSAR_MAX_DETECTORS; X++ ) { /* initialise single-IFO F-stat arrays to zero */
2774 line.avTwoFX[X] = 0.0;
2775 line.maxTwoFXl[X] = 0.0;
2776 line.maxTwoFXlSeg[X] = -1;
2777 line.avTwoFXrecalc[X] = 0.0;
2778 }
2779 line.avTwoFrecalc = -1.0; /* initialise this to -1.0, so that it only gets written out by print_gctFstatline_to_str if later overwritten in recalcToplistStats step */
2780 line.log10BSGLrecalc = -LAL_REAL4_MAX; /* for now, block field with minimal value, needed for output checking in print_gctFstatline_to_str() */
2781 line.log10BSGLtLrecalc = -LAL_REAL4_MAX; /* for now, block field with minimal value, needed for output checking in print_gctFstatline_to_str() */
2782 line.have_f3dot = have_f3dot;
2783 line.loudestSeg = -1;
2784 line.twoFloudestSeg = -1.0;
2785 for ( UINT4 X = 0; X < PULSAR_MAX_DETECTORS; X++ ) {
2786 line.twoFXloudestSeg[X] = -1.0;
2787 }
2788
2789 /* local placeholders for summed 2F value over segments, not averages yet */
2790 REAL4 sumTwoF = in->sumTwoF[ifreq_fg];
2791 REAL4 sumTwoFX[PULSAR_MAX_DETECTORS];
2792 if ( in->sumTwoFX ) { /* if we already have FX values from the main loop, insert these, and calculate BSGL here */
2793 for ( UINT4 X = 0; X < in->numDetectors; X++ ) {
2794 sumTwoFX[X] = in->sumTwoFX[FG_FX_INDEX( *in, X, ifreq_fg )]; /* here it's still the summed 2F value over segments, not the average */
2795 }
2796 xlalErrno = 0;
2797
2798 line.log10BSGL = XLALComputeBSGL( sumTwoF, sumTwoFX, usefulparams->BSGLsetup );
2799 if ( xlalErrno != 0 ) {
2800 XLALPrintError( "%s line %d : XLALComputeBSGL() failed with xlalErrno = %d.\n\n", __func__, __LINE__, xlalErrno );
2801 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2802 }
2803 if ( line.log10BSGL < -LAL_REAL4_MAX * 0.1 ) {
2804 line.log10BSGL = -LAL_REAL4_MAX * 0.1; /* avoid minimum value, needed for output checking in print_gctFstatline_to_str() */
2805 }
2806 } else {
2807 line.log10BSGL = -LAL_REAL4_MAX; /* in non-BSGL case, block field with minimal value, needed for output checking in print_gctFstatline_to_str() */
2808 }
2809
2810 /* take F-stat averages over segments */
2811 line.avTwoF = sumTwoF * NSegmentsInv; /* average multi-2F by full number of segments */
2812 if ( in->sumTwoFX ) {
2813 for ( UINT4 X = 0; X < in->numDetectors; X++ ) {
2814 line.avTwoFX[X] = sumTwoFX[X] * NSegmentsInvX[X]; /* average single-2F by per-IFO number of segments */
2815 }
2816 }
2817
2818 if ( in->maxTwoFXl ) { /* if we already have max-per-segment values from the main loop, insert these too */
2819 line.maxTwoFl = in->maxTwoFl[ifreq_fg];
2820 line.maxTwoFlSeg = in->maxTwoFlIdx[ifreq_fg];
2821 for ( UINT4 X = 0; X < in->numDetectors; X++ ) {
2822 line.maxTwoFXl[X] = in->maxTwoFXl[FG_FX_INDEX( *in, X, ifreq_fg )];
2823 line.maxTwoFXlSeg[X] = in->maxTwoFXlIdx[FG_FX_INDEX( *in, X, ifreq_fg )];
2824 }
2825
2826 line.log10BSGLtL = XLALComputeBSGLtL( sumTwoF, sumTwoFX, line.maxTwoFXl, usefulparams->BSGLsetup );
2827 if ( xlalErrno != 0 ) {
2828 XLALPrintError( "%s line %d : XLALComputeBSGLtL() failed with xlalErrno = %d.\n\n", __func__, __LINE__, xlalErrno );
2829 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2830 }
2831 if ( line.log10BSGLtL < -LAL_REAL4_MAX * 0.1 ) {
2832 line.log10BSGLtL = -LAL_REAL4_MAX * 0.1; /* avoid minimum value, needed for output checking in print_gctFstatline_to_str() */
2833 }
2834
2835 line.log10BtSGLtL = XLALComputeBtSGLtL( line.maxTwoFl, sumTwoFX, line.maxTwoFXl, usefulparams->BSGLsetup );
2836 if ( xlalErrno != 0 ) {
2837 XLALPrintError( "%s line %d : XLALComputeBtSGLtL() failed with xlalErrno = %d.\n\n", __func__, __LINE__, xlalErrno );
2838 ABORT( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
2839 }
2840 if ( line.log10BtSGLtL < -LAL_REAL4_MAX * 0.1 ) {
2841 line.log10BtSGLtL = -LAL_REAL4_MAX * 0.1; /* avoid minimum value, needed for output checking in print_gctFstatline_to_str() */
2842 }
2843
2844 }
2845
2846 insert_into_gctFstat_toplist( list1, &line );
2847 if ( list2 ) { // also insert candidate into (optional) second toplist
2848 insert_into_gctFstat_toplist( list2, &line );
2849 }
2850 if ( list3 ) { // also insert candidate into (optional) 3rd toplist
2851 insert_into_gctFstat_toplist( list3, &line );
2852 }
2853
2854 } // for ifreq_fg
2855
2856 DETATCHSTATUSPTR( status );
2857 RETURN( status );
2858
2859} /* UpdateSemiCohToplists() */
2860
2861
2862// simply return maximal (multi-IFO) Fstat-value over frequency-bins in in->twoF
2863static inline REAL4
2864findLoudestTwoF( const FstatResults *in )
2865{
2866 REAL4 maxTwoF = 0;
2867 for ( UINT4 k = 0; k < in->numFreqBins; k ++ ) {
2868 maxTwoF = fmaxf( maxTwoF, in->twoF[k] );
2869 } // k < numFreqBins
2870
2871 return maxTwoF;
2872} // findLoudestTwoF()
2873
2874
2875/** Print Fstat vectors */
2876void PrintFstatVec( LALStatus *status,
2877 FstatResults *in,
2878 FILE *fp,
2879 PulsarDopplerParams *thisPoint,
2880 LIGOTimeGPS refTime,
2881 INT4 stackIndex )
2882{
2883 INITSTATUS( status );
2884 ATTATCHSTATUSPTR( status );
2885
2886 fprintf( fp, "%% Fstat values from stack %d (reftime -- %d %d)\n", stackIndex, refTime.gpsSeconds, refTime.gpsNanoSeconds );
2887
2888 REAL8 alpha = thisPoint->Alpha;
2889 REAL8 delta = thisPoint->Delta;
2890
2891 PulsarSpins fkdot;
2892 memcpy( fkdot, thisPoint->fkdot, sizeof( fkdot ) );
2893
2894 UINT4 length = in->numFreqBins;
2895 REAL8 deltaF = in->dFreq;
2896
2897 REAL8 f0 = fkdot[0];
2898 for ( UINT4 k = 0; k < length; k++ ) {
2899 fkdot[0] = f0 + k * deltaF;
2900
2901 /* propagate fkdot back to reference-time */
2902 XLAL_CHECK_LAL( status, XLALExtrapolatePulsarSpins( fkdot, fkdot, XLALGPSDiff( &refTime, &thisPoint->refTime ) ) == XLAL_SUCCESS, XLAL_EFUNC );
2903
2904 fprintf( fp, "%d %.13g %.12g %.12g %.13g %.13g %.6g\n",
2905 stackIndex, fkdot[0], alpha, delta, fkdot[1], fkdot[2], in->twoF[k] );
2906 }
2907
2908 fprintf( fp, "\n" );
2909
2910 DETATCHSTATUSPTR( status );
2911 RETURN( status );
2912
2913}
2914
2915
2916
2917
2918
2919
2920
2921
2922
2923
2924
2925
2926
2927
2928
2929/**
2930 * Calculate Earth orbital position, velocity and acceleration
2931 * at midpoint of each segment
2932 */
2933void GetSegsPosVelAccEarthOrb( LALStatus *status,
2934 REAL8VectorSequence **posSeg,
2935 REAL8VectorSequence **velSeg,
2936 REAL8VectorSequence **accSeg,
2937 UsefulStageVariables *usefulparams )
2938{
2939
2940 UINT4 k, nStacks;
2941 LIGOTimeGPSVector *tsMid;
2942 vect3Dlist_t *pvaUR = NULL;
2943
2944 INITSTATUS( status );
2945 ATTATCHSTATUSPTR( status );
2946
2947 ASSERT( usefulparams != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2948 ASSERT( usefulparams->nStacks > 0, status, HIERARCHICALSEARCH_EVAL, HIERARCHICALSEARCH_MSGEVAL );
2949 ASSERT( usefulparams->midTstack != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2950 ASSERT( posSeg != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2951 ASSERT( velSeg != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2952 ASSERT( accSeg != NULL, status, HIERARCHICALSEARCH_ENULL, HIERARCHICALSEARCH_MSGENULL );
2953
2954 /* local copies */
2955 nStacks = usefulparams->nStacks;
2956 tsMid = usefulparams->midTstack;
2957
2958 /* get pos,vel,acc at midpoint of each segment*/
2959 for ( k = 0; k < nStacks; k++ ) {
2960 /* initialize velocities and positions */
2961 posSeg[0]->data[3 * k] = 0.0;
2962 posSeg[0]->data[3 * k + 1] = 0.0;
2963 posSeg[0]->data[3 * k + 2] = 0.0;
2964
2965 velSeg[0]->data[3 * k] = 0.0;
2966 velSeg[0]->data[3 * k + 1] = 0.0;
2967 velSeg[0]->data[3 * k + 2] = 0.0;
2968
2969 accSeg[0]->data[3 * k] = 0.0;
2970 accSeg[0]->data[3 * k + 1] = 0.0;
2971 accSeg[0]->data[3 * k + 2] = 0.0;
2972
2973 /* get Earth's orbital pos vel acc */
2974 if ( ( pvaUR = XLALComputeOrbitalDerivatives( 3, &tsMid->data[k], usefulparams->edat ) ) == NULL ) {
2975 LogPrintf( LOG_CRITICAL, "GetSegsPosVelAccEarthOrb(): XLALComputeOrbitalDerivatives() failed.\n" );
2976 ABORT( status, HIERARCHICALSEARCH_ESFT, HIERARCHICALSEARCH_MSGESFT );
2977 }
2978
2979 posSeg[0]->data[3 * k] = pvaUR->data[0][0];
2980 posSeg[0]->data[3 * k + 1] = pvaUR->data[0][1];
2981 posSeg[0]->data[3 * k + 2] = pvaUR->data[0][2];
2982
2983 velSeg[0]->data[3 * k] = pvaUR->data[1][0];
2984 velSeg[0]->data[3 * k + 1] = pvaUR->data[1][1];
2985 velSeg[0]->data[3 * k + 2] = pvaUR->data[1][2];
2986
2987 accSeg[0]->data[3 * k] = pvaUR->data[2][0];
2988 accSeg[0]->data[3 * k + 1] = pvaUR->data[2][1];
2989 accSeg[0]->data[3 * k + 2] = pvaUR->data[2][2];
2990
2991 XLALDestroyVect3Dlist( pvaUR );
2992
2993 } /* loop over segment -- end pos vel acc of Earth's orbital motion */
2994
2995
2996 DETATCHSTATUSPTR( status );
2997 RETURN( status );
2998
2999} /* GetSegsPosVelAccEarthOrb() */
3000
3001
3002
3003
3004
3005
3006
3007/** Calculate the U1 index for a given point in parameter space */
3008static inline INT4 ComputeU1idx( REAL8 freq_event,
3009 REAL8 f1dot_event,
3010 REAL8 A1,
3011 REAL8 B1,
3012 REAL8 U1start,
3013 REAL8 U1winInv )
3014{
3015 /* compute the index of global-correlation coordinate U1, Eq. (1) */
3016 return ( ( ( freq_event * A1 + f1dot_event * B1 ) - U1start ) * U1winInv ) + 0.5;
3017
3018} /* ComputeU1idx */
3019
3020
3021
3022
3023
3024/** Calculate the U2 index for a given point in parameter space */
3025void ComputeU2idx( REAL8 freq_event,
3026 REAL8 f1dot_event,
3027 REAL8 A2,
3028 REAL8 B2,
3029 REAL8 U2start,
3030 REAL8 U2winInv,
3031 INT4 *U2idx )
3032{
3033
3034 /* compute the index of global-correlation coordinate U2 */
3035 *U2idx = ( INT4 )( ( ( ( f1dot_event + freq_event * A2 + 2.0 * f1dot_event * B2 ) - U2start ) * U2winInv ) + 0.5 );
3036
3037 return;
3038
3039} /* ComputeU2idx */
3040
3041/**
3042 * Set up 'segmented' SFT-catalogs for given list of segments and a total SFT-catalog.
3043 *
3044 * Note: this function does not allow 'empty' segments to be returned, i.e. if there is any
3045 * segment that would contain no SFTs from the given SFT-catalog, an error is returned.
3046 * These segment-lists are 'precomputed' and therefore one can assume that empty segments
3047 * are not intended.
3048 * However, the function will not complain if some SFTs from the catalog are 'unused', i.e.
3049 * they didn't fit into any of the given segments.
3050 *
3051 * \note the input segment list must be sorted, otherwise an error is returned
3052 *
3053 */
3054SFTCatalogSequence *
3055XLALSetUpStacksFromSegmentList( const SFTCatalog *catalog, /**< complete list of SFTs read in */
3056 const LALSegList *segList /**< pre-computed list of segments to split SFTs into */
3057 )
3058{
3059 SFTCatalogSequence *stacks; /* output: segmented SFT-catalogs */
3060
3061 /* check input consistency */
3062 if ( !catalog || !segList ) {
3063 XLALPrintError( "%s: invalid NULL input\n", __func__ );
3064 XLAL_ERROR_NULL( XLAL_EINVAL );
3065 }
3066 /* check that segment list is sorted */
3067 if ( ! segList->sorted ) {
3068 XLALPrintError( "%s: input segment list must be sorted! -> Use XLALSegListSort()\n", __func__ );
3069 XLAL_ERROR_NULL( XLAL_EDOM );
3070 }
3071
3072 UINT4 numSegments = segList->length;
3073 UINT4 numSFTs = catalog->length;
3074
3075 /* set memory of output catalog sequence to maximum possible length */
3076 if ( ( stacks = XLALCalloc( 1, sizeof( *stacks ) ) ) == NULL ) {
3077 XLALPrintError( "%s: XLALCalloc(%zu) failed.\n", __func__, sizeof( *stacks ) );
3078 XLAL_ERROR_NULL( XLAL_ENOMEM );
3079 }
3080 stacks->length = numSegments;
3081 if ( ( stacks->data = XLALCalloc( stacks->length, sizeof( *stacks->data ) ) ) == NULL ) {
3082 XLALPrintError( "%s: failed to allocate segmented SFT-catalog\n", __func__ );
3083 XLAL_ERROR_NULL( XLAL_ENOMEM );
3084 }
3085
3086 /* Step through segment list:
3087 * for every segment:
3088 * - find earliest and last SFT *starting* within given segment
3089 * this ensures we use all SFTs and dont lose some that fall on segment boundaries
3090 * - copy this range of SFT-headers into the segmented SFT-catalog 'stacks'
3091 */
3092 UINT4 iSeg;
3093 INT4 iSFT0 = 0, iSFT1 = 0; /* indices of earliest and last SFT fitting into segment iSeg */
3094 for ( iSeg = 0; iSeg < numSegments; iSeg ++ ) {
3095 LALSeg *thisSeg = &( segList->segs[iSeg] );
3096
3097 /* ----- find earliest SFT fitting into this segment */
3098 iSFT0 = iSFT1; /* start from previous segment's last SFT */
3099 while ( 1 ) {
3100 int cmp = XLALCWGPSinRange( catalog->data[iSFT0].header.epoch, &thisSeg->start, &thisSeg->end );
3101
3102 if ( cmp < 0 ) { /* iSFT0 lies *before* current segment => advance */
3103 iSFT0 ++;
3104 }
3105 if ( cmp == 0 ) { /* iSFT0 lies *inside* current segment ==> stop */
3106 break;
3107 }
3108
3109 /* if no more SFTs or iSFT0 lies *past* current segment => ERROR: empty segment! */
3110 if ( cmp > 0 || iSFT0 == ( INT4 )numSFTs ) {
3111 XLALPrintError( "%s: Empty segment! No SFTs fit into segment iSeg=%d\n", __func__, iSeg );
3112 XLAL_ERROR_NULL( XLAL_EDOM );
3113 }
3114 } /* while true */
3115
3116 /* ----- find last SFT still starting within segment */
3117 iSFT1 = iSFT0;
3118 while ( 1 ) {
3119 int cmp = XLALCWGPSinRange( catalog->data[iSFT1].header.epoch, &thisSeg->start, &thisSeg->end );
3120
3121 if ( cmp < 0 ) { /* start of iSFT1 lies *before* current segment ==> something is screwed up! */
3122 XLALPrintError( "%s: start of current SFT %d lies before current segment %d ==> code seems inconsistent!\n", __func__, iSFT1, iSeg );
3123 XLAL_ERROR_NULL( XLAL_EFAILED );
3124 }
3125 if ( cmp == 0 ) { /* start of iSFT1 lies *inside* current segment ==> advance */
3126 iSFT1 ++;
3127 }
3128
3129 if ( cmp > 0 || iSFT1 == ( INT4 )numSFTs ) { /* last SFT reached or start of iSFT1 lies *past* current segment => step back once and stop */
3130 iSFT1 --;
3131 break;
3132 }
3133
3134 } /* while true */
3135
3136 INT4 numSFTsInSeg = iSFT1 - iSFT0 + 1;
3137
3138 /* ----- allocate and copy this range of SFTs into the segmented catalog */
3139 stacks->data[iSeg].length = ( UINT4 )numSFTsInSeg;
3140 UINT4 size = sizeof( *stacks->data[iSeg].data );
3141 if ( ( stacks->data[iSeg].data = XLALCalloc( numSFTsInSeg, size ) ) == NULL ) {
3142 XLALPrintError( "%s: failed to XLALCalloc(%d, %d)\n", __func__, numSFTsInSeg, size );
3143 XLAL_ERROR_NULL( XLAL_ENOMEM );
3144 }
3145
3146 INT4 iSFT;
3147 for ( iSFT = iSFT0; iSFT <= iSFT1; iSFT++ ) {
3148 stacks->data[iSeg].data[iSFT - iSFT0] = catalog->data[iSFT];
3149 }
3150
3151 } /* for iSeg < numSegments */
3152
3153 return stacks;
3154
3155} /* XLALSetUpStacksFromSegmentList() */
3156
3157
3158
3159
3160/**
3161 * XLAL function to extrapolate the pulsar spin parameters of all toplist candidates
3162 * from reftime of the input toplist ('inRefTime') to a user-specified output reftime 'outRefTime'
3163 */
3164int XLALExtrapolateToplistPulsarSpins( toplist_t *list, /**< [out/in] toplist with GCTtopOutputEntry items, 'Freq,F1dot,F2dot,F3dot' fields will be overwritten */
3165 const LIGOTimeGPS outRefTime, /**< reference time as requested for the final candidate output */
3166 const LIGOTimeGPS inRefTime /**< reference time of the input toplist */
3167 )
3168{
3169
3170 /* check input parameters */
3171 if ( !list ) {
3172 XLAL_ERROR( XLAL_EFAULT, "\nNULL pointer given instead of toplist.\n" );
3173 }
3174
3175 /* convert LIGOTimeGPS into real number difference for XLALExtrapolatePulsarSpins */
3176 REAL8 deltaTau = XLALGPSDiff( &outRefTime, &inRefTime );
3177
3178 /* check if translation to reference time of toplist is necessary */
3179 if ( deltaTau == 0 ) {
3180 return XLAL_SUCCESS; /* can skip this step if reftimes are equal */
3181 }
3182
3183 PulsarSpins XLAL_INIT_DECL( fkdot );
3184
3185 UINT4 numElements = list->elems;
3186 for ( UINT4 j = 0; j < numElements; j++ ) { /* loop over toplist */
3187 /* get fkdot of each candidate */
3188 GCTtopOutputEntry *elem = toplist_elem( list, j );
3189 fkdot[0] = elem->Freq;
3190 fkdot[1] = elem->F1dot;
3191 fkdot[2] = elem->F2dot;
3192 fkdot[3] = elem->F3dot;
3193 /* propagate fkdot to reference-time */
3194 if ( XLALExtrapolatePulsarSpins( fkdot, fkdot, deltaTau ) != XLAL_SUCCESS ) {
3195 XLALPrintError( "\n%s, line %d : XLALExtrapolatePulsarSpins() failed.\n\n", __func__, __LINE__ );
3196 XLAL_ERROR( XLAL_EFUNC );
3197 }
3198 /* write back propagated frequency to toplist */
3199 elem->Freq = fkdot[0];
3200 elem->F1dot = fkdot[1];
3201 elem->F2dot = fkdot[2];
3202 elem->F3dot = fkdot[3];
3203 }
3204
3205 return XLAL_SUCCESS;
3206
3207} /* XLALExtrapolateToplistPulsarSpins() */
3208
3209
3210/**
3211 * Function to append one timing-info line to output file.
3212 *
3213 */
3214static int
3215write_TimingInfo( const CHAR *fname, const timingInfo_t *ti )
3216{
3217 /* input sanity */
3218 if ( !fname || !ti ) {
3219 XLALPrintError( "%s: invalid NULL input 'fp' | 'ti'\n", __func__ );
3220 XLAL_ERROR( XLAL_EINVAL );
3221 }
3222
3223 FILE *fp;
3224 if ( ( fp = fopen( fname, "rb" ) ) == NULL ) {
3225 XLAL_CHECK( ( fp = fopen( fname, "wb" ) ) != NULL, XLAL_ESYS, "Failed to open new timing-file '%s' for writing\n", fname );
3226 fprintf( fp, "%%%%--------------------------------------------------------------------------------\n" );
3227 fprintf( fp, "%%%% GCT Timing model:\n" );
3228 fprintf( fp, "%%%% runtime = Nseg * Ndet * Ncoh * tau_Fstat + Nseg * Ninc * tau_SumF + Ninc * tau_Bayes + Ncan * tau_Recalc + time_Other\n" );
3229 fprintf( fp, "%%%%--------------------------------------------------------------------------------\n" );
3230 fprintf( fp, "%%%%\n" );
3231 fprintf( fp, "%2s%10s %10s %10s %10s %10s | %6s %6s %10s %6s %10s %10s %10s %10s %%%15s %15s\n",
3232 "%%", "tau_Fstat", "tau_SumF", "tau_Bayes", "tau_Recalc", "time_Other",
3233 "Nseg", "Ndet", "Tcoh[s]", "Nsft", "NFreqCo", "Ncoh", "Ninc", "Ncand", "FstatMethod", "RecalcMethod" );
3234 } else {
3235 fclose( fp );
3236 XLAL_CHECK( ( fp = fopen( fname, "ab" ) ) != NULL, XLAL_ESYS, "Failed to open existing timing-file '%s' for appending\n", fname );
3237 }
3238
3239 fprintf( fp, "%12.1e %10.1e %10.1e %10.1e %10.1e %6d %6d %10d %6d %10d %10.1e %10.1e %10d %%%15s %15s\n",
3240 ti->tau_Fstat, ti->tau_SumF, ti->tau_Bayes, ti->tau_Recalc, ti->time_Other,
3241 ti->Nseg, ti->Ndet, ti->Tcoh, ti->Nsft, ti->NFreqCo, ti->Ncoh, ti->Ninc, ti->Ncand, ti->FstatMethodStr, ti->RecalcMethodStr );
3242
3243 fclose( fp );
3244 return XLAL_SUCCESS;
3245
3246} // write_TimingInfo()
likely
#define likely(x)
Definition: ComputeFstat_Demod_ComputeFaFb.c:41
unlikely
#define unlikely(x)
Definition: ComputeFstat_internal.h:42
__func__
#define __func__
log an I/O error, i.e.
Definition: CrossCorrToplist.c:865
InitDopplerSkyScan
void InitDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan, const DopplerSkyScanInit *init)
Definition: DopplerScan.c:273
FreeDopplerSkyScan
void FreeDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan)
Definition: DopplerScan.c:324
XLALNextDopplerSkyPos
int XLALNextDopplerSkyPos(PulsarDopplerParams *pos, DopplerSkyScanState *skyScan)
NextDopplerSkyPos(): step through sky-grid return 0 = OK, -1 = ERROR.
Definition: DopplerScan.c:127
STATE_FINISHED
@ STATE_FINISHED
all templates have been read
Definition: DopplerScan.h:84
GRID_METRIC
@ GRID_METRIC
generate grid using a 2D sky-metric
Definition: DopplerScan.h:93
create_gctFstat_toplist
int create_gctFstat_toplist(toplist_t **tl, UINT8 length, SortBy_t whatToSortBy)
creates a toplist with length elements, returns -1 on error (usually out of memory),...
Definition: GCTtoplist.c:246
free_gctFstat_toplist
void free_gctFstat_toplist(toplist_t **l)
frees the space occupied by the toplist
Definition: GCTtoplist.c:270
clear_gct_checkpoint
int clear_gct_checkpoint(const char *filename)
removes a checkpoint returns 0 on success, errno on failure
Definition: GCTtoplist.c:1242
insert_into_gctFstat_toplist
int insert_into_gctFstat_toplist(toplist_t *tl, GCTtopOutputEntry *elem)
Inserts an element in to the toplist either if there is space left or the element is larger than the ...
Definition: GCTtoplist.c:281
write_hfs_oputput
int write_hfs_oputput(const char *filename, toplist_t *tl)
write the final output file:
Definition: GCTtoplist.c:1284
SortBy_t
SortBy_t
enumerate all toplist-sorting options: by F (0), number-count (1), BSGL (2), "dual" toplists F + BSGL...
Definition: GCTtoplist.h:71
SORTBY_BtSGLtL
@ SORTBY_BtSGLtL
Definition: GCTtoplist.h:77
SORTBY_LAST
@ SORTBY_LAST
Definition: GCTtoplist.h:80
SORTBY_DUAL_F_BSGL
@ SORTBY_DUAL_F_BSGL
Definition: GCTtoplist.h:75
SORTBY_BSGL
@ SORTBY_BSGL
Definition: GCTtoplist.h:74
SORTBY_F
@ SORTBY_F
Definition: GCTtoplist.h:72
SORTBY_BSGLtL
@ SORTBY_BSGLtL
Definition: GCTtoplist.h:76
SORTBY_F_BSGLtL_BtSGLtL
@ SORTBY_F_BSGLtL_BtSGLtL
Definition: GCTtoplist.h:79
SORTBY_TRIPLE_BStSGLtL
@ SORTBY_TRIPLE_BStSGLtL
Definition: GCTtoplist.h:78
toplist_elem
void * toplist_elem(toplist_t *list, size_t ind)
Definition: HeapToplist.c:185
TEST_FSTAT_TOPLIST_INCLUSION
#define TEST_FSTAT_TOPLIST_INCLUSION(list, element)
Definition: HeapToplist.h:68
XLALReadSegmentsFromFile_support_4column_format
int XLALReadSegmentsFromFile_support_4column_format
Definition: SFTtimestamps.c:38
main
int main(int argc, char *argv[])
Definition: HierarchSearchGCT.c:239
ALRealloc
#define ALRealloc
Definition: HierarchSearchGCT.c:51
compareFineGridsumTwoF
int compareFineGridsumTwoF(const void *a, const void *b)
UpdateSemiCohToplists
void UpdateSemiCohToplists(LALStatus *status, toplist_t *list1, toplist_t *list2, toplist_t *list3, FineGrid *in, REAL8 f1dot_fg, REAL8 f2dot_fg, REAL8 f3dot_fg, UsefulStageVariables *usefulparams, REAL4 NSegmentsInv, REAL4 *NSegmentsInvX, BOOLEAN have_f3dot)
Get SemiCoh candidates into toplist(s) This function allows for inserting candidates into up to 3 top...
Definition: HierarchSearchGCT.c:2727
FSTART
#define FSTART
Default Start search frequency.
Definition: HierarchSearchGCT.c:73
PrintFstatVec
void PrintFstatVec(LALStatus *status, FstatResults *in, FILE *fp, PulsarDopplerParams *thisPoint, LIGOTimeGPS refTime, INT4 stackIndex)
Print Fstat vectors.
Definition: HierarchSearchGCT.c:2876
MISMATCH
#define MISMATCH
Default for metric grid maximal mismatch value.
Definition: HierarchSearchGCT.c:82
F3DOT
#define F3DOT
Default value of third spindown.
Definition: HierarchSearchGCT.c:78
DDELTA
#define DDELTA
Default resolution for isotropic or flat grids.
Definition: HierarchSearchGCT.c:84
global_argv
char ** global_argv
Definition: HierarchSearchGCT.c:70
SetUpStacks
void SetUpStacks(LALStatus *status, SFTCatalogSequence *out, REAL8 tStack, SFTCatalog *in, UINT4 nStacksMax)
Breaks up input sft catalog into specified number of stacks Loops over elements of the catalog,...
Definition: HierarchSearchGCT.c:2390
DF3DOT
#define DF3DOT
Default range of third spindown parameter.
Definition: HierarchSearchGCT.c:80
ALFree
#define ALFree
Definition: HierarchSearchGCT.c:52
XLALSetUpStacksFromSegmentList
SFTCatalogSequence * XLALSetUpStacksFromSegmentList(const SFTCatalog *catalog, const LALSegList *segList)
Set up 'segmented' SFT-catalogs for given list of segments and a total SFT-catalog.
Definition: HierarchSearchGCT.c:3055
SET_GCT_CHECKPOINT
#define SET_GCT_CHECKPOINT
Definition: HierarchSearchGCT.c:69
ALIGN_REAL4
#define ALIGN_REAL4
DALPHA
#define DALPHA
Default resolution for isotropic or flat grids.
Definition: HierarchSearchGCT.c:83
GET_GCT_CHECKPOINT
#define GET_GCT_CHECKPOINT
Definition: HierarchSearchGCT.c:68
GetSegsPosVelAccEarthOrb
void GetSegsPosVelAccEarthOrb(LALStatus *status, REAL8VectorSequence **posSeg, REAL8VectorSequence **velSeg, REAL8VectorSequence **accSeg, UsefulStageVariables *usefulparams)
Calculate Earth orbital position, velocity and acceleration at midpoint of each segment.
Definition: HierarchSearchGCT.c:2933
LAL_INT4_MAX
#define LAL_INT4_MAX
Definition: HierarchSearchGCT.c:89
global_status
LALStatus * global_status
Definition: HierarchSearchGCT.c:230
compareFineGridNC
int compareFineGridNC(const void *a, const void *b)
findLoudestTwoF
static REAL4 findLoudestTwoF(const FstatResults *in)
Definition: HierarchSearchGCT.c:2864
NUDGE
#define NUDGE
Definition: HierarchSearchGCT.c:57
global_column_headings_stringp
char * global_column_headings_stringp
Definition: HierarchSearchGCT.c:231
UpdateSemiCohToplistsOptimTriple
void UpdateSemiCohToplistsOptimTriple(LALStatus *status, SortBy_t toplist_sortby, toplist_t *list1, toplist_t *list2, toplist_t *list3, FineGrid *in, REAL8 f1dot_fg, REAL8 f2dot_fg, REAL8 f3dot_fg, UsefulStageVariables *usefulparams, REAL4 NSegmentsInv, REAL4 *NSegmentsInvX, BOOLEAN have_f3dot)
Get SemiCoh candidates into toplist(s) This function allows for inserting candidates into up to 3 top...
Definition: HierarchSearchGCT.c:2564
compareCoarseGridUindex
int compareCoarseGridUindex(const void *a, const void *b)
SKYREGION
#define SKYREGION
default sky region to search over – just a single point
Definition: HierarchSearchGCT.c:81
F2DOT
#define F2DOT
Default value of second spindown.
Definition: HierarchSearchGCT.c:77
global_argc
int global_argc
Definition: HierarchSearchGCT.c:71
TRUE
#define TRUE
Definition: HierarchSearchGCT.c:59
FALSE
#define FALSE
Definition: HierarchSearchGCT.c:60
DF2DOT
#define DF2DOT
Default range of second spindown parameter.
Definition: HierarchSearchGCT.c:79
ComputeU1idx
static INT4 ComputeU1idx(REAL8 freq_event, REAL8 f1dot_event, REAL8 A1, REAL8 B1, REAL8 U1start, REAL8 U1winInv)
Calculate the U1 index for a given point in parameter space.
Definition: HierarchSearchGCT.c:3008
FDOT
#define FDOT
Default value of first spindown.
Definition: HierarchSearchGCT.c:75
XLALExtrapolateToplistPulsarSpins
int XLALExtrapolateToplistPulsarSpins(toplist_t *list, const LIGOTimeGPS usefulParamsRefTime, const LIGOTimeGPS finegridRefTime)
XLAL function to extrapolate the pulsar spin parameters of all toplist candidates from reftime of the...
Definition: HierarchSearchGCT.c:3164
NCAND1
#define NCAND1
Default number of candidates to be followed up from first stage.
Definition: HierarchSearchGCT.c:86
HSMAX
#define HSMAX(x, y)
Definition: HierarchSearchGCT.c:100
PrintStackInfo
void PrintStackInfo(LALStatus *status, const SFTCatalogSequence *catalogSeq, FILE *fp)
Print some stack info from sft catalog sequence.
Definition: HierarchSearchGCT.c:2522
DFDOT
#define DFDOT
Default range of first spindown parameter.
Definition: HierarchSearchGCT.c:76
FBAND
#define FBAND
Default search band.
Definition: HierarchSearchGCT.c:74
ComputeU2idx
void ComputeU2idx(REAL8 freq_event, REAL8 f1dot_event, REAL8 A2, REAL8 B2, REAL8 U2start, REAL8 U2winInv, INT4 *U2idx)
Calculate the U2 index for a given point in parameter space.
Definition: HierarchSearchGCT.c:3025
write_TimingInfo
static int write_TimingInfo(const CHAR *fname, const timingInfo_t *ti)
Function to append one timing-info line to output file.
Definition: HierarchSearchGCT.c:3215
FNAMEOUT
#define FNAMEOUT
Default output file basename.
Definition: HierarchSearchGCT.c:87
FSTATTHRESHOLD
#define FSTATTHRESHOLD
Default threshold on Fstatistic for peak selection.
Definition: HierarchSearchGCT.c:85
HSMIN
#define HSMIN(x, y)
Definition: HierarchSearchGCT.c:101
PrintCatalogInfo
void PrintCatalogInfo(LALStatus *status, const SFTCatalog *catalog, FILE *fp)
Print some sft catalog info.
Definition: HierarchSearchGCT.c:2489
GETTIME
#define GETTIME()
Definition: HierarchSearchGCT.c:103
SetUpSFTs
void SetUpSFTs(LALStatus *status, UsefulStageVariables *in)
Set up stacks, read SFTs, calculate SFT noise weights and calculate detector-state.
Definition: HierarchSearchGCT.c:2013
HIERARCHICALSEARCH_EMEM
#define HIERARCHICALSEARCH_EMEM
Definition: HierarchSearchGCT.h:89
FG_INDEX
#define FG_INDEX(fg, iFreq)
Definition: HierarchSearchGCT.h:175
HIERARCHICALSEARCH_MSGEXLAL
#define HIERARCHICALSEARCH_MSGEXLAL
Definition: HierarchSearchGCT.h:105
HIERARCHICALSEARCH_MSGENULL
#define HIERARCHICALSEARCH_MSGENULL
Definition: HierarchSearchGCT.h:99
FG_FX_INDEX
#define FG_FX_INDEX(fg, iDet, iFreq)
Definition: HierarchSearchGCT.h:183
HIERARCHICALSEARCH_ESUB
#define HIERARCHICALSEARCH_ESUB
Definition: HierarchSearchGCT.h:79
HIERARCHICALSEARCH_MSGEFILE
#define HIERARCHICALSEARCH_MSGEFILE
Definition: HierarchSearchGCT.h:98
HIERARCHICALSEARCH_ECHECKPT
#define HIERARCHICALSEARCH_ECHECKPT
Definition: HierarchSearchGCT.h:88
HIERARCHICALSEARCH_ENULL
#define HIERARCHICALSEARCH_ENULL
Definition: HierarchSearchGCT.h:83
HIERARCHICALSEARCH_MSGESUB
#define HIERARCHICALSEARCH_MSGESUB
Definition: HierarchSearchGCT.h:95
HIERARCHICALSEARCH_MSGECHECKPT
#define HIERARCHICALSEARCH_MSGECHECKPT
Definition: HierarchSearchGCT.h:102
HIERARCHICALSEARCH_MSGESFT
#define HIERARCHICALSEARCH_MSGESFT
Definition: HierarchSearchGCT.h:104
FINEGRID_NC_T
#define FINEGRID_NC_T
structure for storing fine-grid points
Definition: HierarchSearchGCT.h:151
HIERARCHICALSEARCH_MSGEBAD
#define HIERARCHICALSEARCH_MSGEBAD
Definition: HierarchSearchGCT.h:97
HIERARCHICALSEARCH_EFILE
#define HIERARCHICALSEARCH_EFILE
Definition: HierarchSearchGCT.h:82
HIERARCHICALSEARCH_EBAD
#define HIERARCHICALSEARCH_EBAD
Definition: HierarchSearchGCT.h:81
HIERARCHICALSEARCH_EVAL
#define HIERARCHICALSEARCH_EVAL
Definition: HierarchSearchGCT.h:84
HIERARCHICALSEARCH_ECG
#define HIERARCHICALSEARCH_ECG
Definition: HierarchSearchGCT.h:91
HIERARCHICALSEARCH_MSGEVAL
#define HIERARCHICALSEARCH_MSGEVAL
Definition: HierarchSearchGCT.h:100
HIERARCHICALSEARCH_ENORM
#define HIERARCHICALSEARCH_ENORM
Definition: HierarchSearchGCT.h:78
HIERARCHICALSEARCH_EXLAL
#define HIERARCHICALSEARCH_EXLAL
Definition: HierarchSearchGCT.h:92
CG_INDEX
#define CG_INDEX(cg, iStack, iFreq)
Definition: HierarchSearchGCT.h:202
HIERARCHICALSEARCH_MSGEMEM
#define HIERARCHICALSEARCH_MSGEMEM
Definition: HierarchSearchGCT.h:103
CG_FX_INDEX
#define CG_FX_INDEX(cg, iDet, iStack, iFreq)
Definition: HierarchSearchGCT.h:208
HIERARCHICALSEARCH_ESFT
#define HIERARCHICALSEARCH_ESFT
Definition: HierarchSearchGCT.h:90
lal_errhandler
lal_errhandler_t lal_errhandler
LAL_ERR_EXIT
int LAL_ERR_EXIT(LALStatus *stat, const char *func, const char *file, const int line, volatile const char *id)
LAL_CALL
#define LAL_CALL(function, statusptr)
j
int j
k
int k
LALCheckMemoryLeaks
void LALCheckMemoryLeaks(void)
LALRealloc
#define LALRealloc(p, n)
LALCalloc
#define LALCalloc(m, n)
LALFree
#define LALFree(p)
lalPulsarVCSInfoList
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
delta
static double double delta
ABORT
#define ABORT(statusptr, code, mesg)
XLAL_CHECK_LAL
#define XLAL_CHECK_LAL(sp, assertion,...)
TRY
#define TRY(func, statusptr)
ATTATCHSTATUSPTR
#define ATTATCHSTATUSPTR(statusptr)
ASSERT
#define ASSERT(assertion, statusptr, code, mesg)
DETATCHSTATUSPTR
#define DETATCHSTATUSPTR(statusptr)
INITSTATUS
#define INITSTATUS(statusptr)
RETURN
#define RETURN(statusptr)
STRING
#define STRING(a)
XLALComputeExtraStatsForToplist
int XLALComputeExtraStatsForToplist(toplist_t *list, const RecalcStatsParams *recalcParams)
XLAL function to go through a (Hough or GCT) toplist and compute line-robust statistics for each cand...
Definition: RecalcToplistStats.c:43
RecalcToplistStats.h
Functions to recompute statistics for GCT/hough toplist entries.
fprintf
#define fprintf
l
int l
e
double e
gc_hotloop_sse2.h
gc_hotloop_2Fmax_tracking
static void gc_hotloop_2Fmax_tracking(REAL4 *fgrid2F, REAL4 *fgrid2Fmax, UINT4 *fgrid2FmaxIdx, REAL4 *cgrid2F, UINT4 k, UINT4 length) __attribute__((hot))
Definition: gc_hotloop_sse2.h:43
gc_hotloop
static void gc_hotloop(REAL4 *fgrid2F, REAL4 *cgrid2F, UCHAR *fgridnc, REAL4 TwoFthreshold, UINT4 length) __attribute__((hot))
Definition: gc_hotloop_sse2.h:261
gc_hotloop_no_nc
static void gc_hotloop_no_nc(REAL4 *fgrid2F, REAL4 *cgrid2F, UINT4 length) __attribute__((hot))
Definition: gc_hotloop_sse2.h:384
XLALCreateREAL8VectorSequence
REAL8VectorSequence * XLALCreateREAL8VectorSequence(UINT4 length, UINT4 veclen)
XLALDestroyREAL8VectorSequence
void XLALDestroyREAL8VectorSequence(REAL8VectorSequence *vecseq)
XLALPulsarParamsFromUserInput
PulsarParamsVector * XLALPulsarParamsFromUserInput(const LALStringVector *UserInput, const LIGOTimeGPS *refTimeDef)
Function to determine the PulsarParamsVector input from a user-input defining CW sources.
Definition: CWMakeFakeData.c:1065
InjectionSourcesHelpString
const char *const InjectionSourcesHelpString
Definition: CWMakeFakeData.c:58
XLALDestroyPulsarParamsVector
void XLALDestroyPulsarParamsVector(PulsarParamsVector *ppvect)
Destructor for PulsarParamsVector type.
Definition: CWMakeFakeData.c:740
XLALGetFstatInputTimestamps
const MultiLIGOTimeGPSVector * XLALGetFstatInputTimestamps(const FstatInput *input)
Returns the SFT timestamps stored in a FstatInput structure.
Definition: ComputeFstat.c:701
FstatMethodType
FstatMethodType
Different methods available to compute the -statistic, falling into two broad classes:
Definition: ComputeFstat.h:111
XLALGetFstatInputDetectors
const MultiLALDetector * XLALGetFstatInputDetectors(const FstatInput *input)
Returns the detector information stored in a FstatInput structure.
Definition: ComputeFstat.c:687
XLALGetFstatInputMethodName
const CHAR * XLALGetFstatInputMethodName(const FstatInput *input)
Returns the human-readable name of the -statistic method being used by a FstatInput structure.
Definition: ComputeFstat.c:672
FstatOptionalArgsDefaults
const FstatOptionalArgs FstatOptionalArgsDefaults
Global initializer for setting FstatOptionalArgs to default values.
Definition: ComputeFstat.c:90
XLALFstatMethodChoices
const UserChoices * XLALFstatMethodChoices(void)
Return pointer to a static array of all (available) FstatMethodType choices.
Definition: ComputeFstat.c:1147
XLALDestroyFstatInputVector
void XLALDestroyFstatInputVector(FstatInputVector *inputs)
Free all memory associated with a FstatInputVector structure.
Definition: ComputeFstat.c:252
XLALAppendFstatTiming2File
int XLALAppendFstatTiming2File(const FstatInput *input, FILE *fp, BOOLEAN printHeader)
Definition: ComputeFstat.c:1225
XLALCreateFstatInputVector
FstatInputVector * XLALCreateFstatInputVector(const UINT4 length)
Create a FstatInputVector of the given length, for example for setting up F-stat searches over severa...
Definition: ComputeFstat.c:231
FstatQuantities
FstatQuantities
Bit-field of -statistic quantities which can be computed by XLALComputeFstat().
Definition: ComputeFstat.h:94
XLALCreateFstatInput
FstatInput * XLALCreateFstatInput(const SFTCatalog *SFTcatalog, const REAL8 minCoverFreq, const REAL8 maxCoverFreq, const REAL8 dFreq, const EphemerisData *ephemerides, const FstatOptionalArgs *optionalArgs)
Create a fully-setup FstatInput structure for computing the -statistic using XLALComputeFstat().
Definition: ComputeFstat.c:359
XLALComputeFstat
int XLALComputeFstat(FstatResults **Fstats, FstatInput *input, const PulsarDopplerParams *doppler, const UINT4 numFreqBins, const FstatQuantities whatToCompute)
Compute the -statistic over a band of frequencies.
Definition: ComputeFstat.c:760
XLALDestroyFstatResults
void XLALDestroyFstatResults(FstatResults *Fstats)
Free all memory associated with a FstatResults structure.
Definition: ComputeFstat.c:934
FSTATQ_2F
@ FSTATQ_2F
Compute multi-detector .
Definition: ComputeFstat.h:96
FSTATQ_2F_PER_DET
@ FSTATQ_2F_PER_DET
Compute for each detector.
Definition: ComputeFstat.h:98
XLALParseMultiNoiseFloorMapped
int XLALParseMultiNoiseFloorMapped(MultiNoiseFloor *multiNoiseFloor, const LALStringVector *multiNoiseFloorDetNames, const LALStringVector *sqrtSX, const LALStringVector *sqrtSXDetNames)
Parse string-vectors (typically input by user) of N detector noise-floors for detectors ,...
Definition: DetectorStates.c:616
XLALExtrapolatePulsarSpinRange
int XLALExtrapolatePulsarSpinRange(PulsarSpinRange *range1, const PulsarSpinRange *range0, const REAL8 dtau)
General pulsar-spin extrapolation function: given a "spin-range" (ie spins + spin-bands) range0 at ti...
Definition: ExtrapolatePulsarSpins.c:99
XLALExtrapolatePulsarSpins
int XLALExtrapolatePulsarSpins(PulsarSpins fkdot1, const PulsarSpins fkdot0, REAL8 dtau)
Extrapolate the Pulsar spin-parameters (fkdot0) from the initial reference-epoch to the new referen...
Definition: ExtrapolatePulsarSpins.c:159
XLALCWSignalCoveringBand
int XLALCWSignalCoveringBand(REAL8 *minCoverFreq, REAL8 *maxCoverFreq, const LIGOTimeGPS *time1, const LIGOTimeGPS *time2, const PulsarSpinRange *spinRange, const REAL8 binaryMaxAsini, const REAL8 binaryMinPeriod, const REAL8 binaryMaxEcc)
Determines a frequency band which covers the frequency evolution of a band of CW signals between two ...
Definition: ExtrapolatePulsarSpins.c:249
XLALInitBarycenter
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
Definition: LALInitBarycenter.c:206
XLALDestroyEphemerisData
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
Definition: LALInitBarycenter.c:322
LAL_REAL4_MAX
#define LAL_REAL4_MAX
BOOLEAN
unsigned char BOOLEAN
XLAL_INIT_MEM
#define XLAL_INIT_MEM(x)
UINT8
uint64_t UINT8
REAL8
double REAL8
XLAL_INIT_DECL
#define XLAL_INIT_DECL(var,...)
CHAR
char CHAR
UINT4
uint32_t UINT4
INT4
int32_t INT4
REAL4
float REAL4
XLALCalloc
void * XLALCalloc(size_t m, size_t n)
XLALFree
void XLALFree(void *p)
LAL_REAL4_FORMAT
#define LAL_REAL4_FORMAT
LAL_UINT8_FORMAT
#define LAL_UINT8_FORMAT
XLALStringDuplicate
char char * XLALStringDuplicate(const char *s)
XLALStringAppend
int char * XLALStringAppend(char *s, const char *append)
XLALVCSInfoString
char * XLALVCSInfoString(const LALVCSInfoList vcs_list, const int verbose, const char *prefix)
XLALComputeBSGLtL
REAL4 XLALComputeBSGLtL(const REAL4 twoF, const REAL4 twoFX[PULSAR_MAX_DETECTORS], const REAL4 maxtwoFlX[PULSAR_MAX_DETECTORS], const BSGLSetup *setup)
Single-bin wrapper of XLALVectorComputeBSGLtL(), provided for backwards compatibility.
Definition: LineRobustStats.c:343
XLALComputeBtSGLtL
REAL4 XLALComputeBtSGLtL(const REAL4 maxtwoFl, const REAL4 twoFX[PULSAR_MAX_DETECTORS], const REAL4 maxtwoFXl[PULSAR_MAX_DETECTORS], const BSGLSetup *setup)
Single-bin wrapper of XLALVectorComputeBtSGLtL(), provided for backwards compatibility.
Definition: LineRobustStats.c:413
XLALDestroyBSGLSetup
void XLALDestroyBSGLSetup(BSGLSetup *setup)
Definition: LineRobustStats.c:121
XLALCreateBSGLSetup
BSGLSetup * XLALCreateBSGLSetup(const UINT4 numDetectors, const REAL4 Fstar0sc, const REAL4 oLGX[PULSAR_MAX_DETECTORS], const BOOLEAN useLogCorrection, const UINT4 numSegments)
Definition: LineRobustStats.c:74
XLALComputeBSGL
REAL4 XLALComputeBSGL(const REAL4 twoF, const REAL4 twoFX[PULSAR_MAX_DETECTORS], const BSGLSetup *setup)
Single-bin wrapper of XLALVectorComputeBSGL(), provided for backwards compatibility.
Definition: LineRobustStats.c:218
XLALParseLinePriors
int XLALParseLinePriors(REAL4 oLGX[PULSAR_MAX_DETECTORS], const LALStringVector *oLGX_string)
Parse string-vectors (typically input by user) of N per-detector line-to-Gaussian prior ratios to a ...
Definition: LineRobustStats.c:501
LogPrintfVerbatim
void void LogPrintfVerbatim(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
LogPrintf
void LogPrintf(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
LOG_CRITICAL
LOG_CRITICAL
LOG_DEBUG
LOG_DEBUG
LOG_NORMAL
LOG_NORMAL
LOG_DETAIL
LOG_DETAIL
LAL_PMETRIC_COH_PTOLE_ANALYTIC
@ LAL_PMETRIC_COH_PTOLE_ANALYTIC
Definition: PtoleMetric.h:97
PULSAR_MAX_DETECTORS
#define PULSAR_MAX_DETECTORS
maximal number of detectors we can handle (for static arrays of detector quantities)
Definition: PulsarDataTypes.h:112
PulsarSpins
REAL8 PulsarSpins[PULSAR_MAX_SPINS]
Typedef for fixed-size array holding GW frequency and derivatives fk = d^k Freq/dt^k|(tau_ref)
Definition: PulsarDataTypes.h:127
a
static const INT4 a
XLALGetSiteInfo
const LALDetector * XLALGetSiteInfo(const CHAR *name)
Find the site geometry-information 'LALDetector' for given a detector name (or prefix).
Definition: SFTnaming.c:218
XLALDestroySFTCatalog
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
XLALCheckCRCSFTCatalog
int XLALCheckCRCSFTCatalog(BOOLEAN *crc_check, SFTCatalog *catalog)
This function reads in the SFTs in the catalog and validates their CRC64 checksums.
Definition: SFTcatalog.c:524
XLALReadSegmentsFromFile
LALSegList * XLALReadSegmentsFromFile(const char *fname)
Function to read a segment list from given filename, returns a sorted LALSegList.
Definition: SFTtimestamps.c:598
XLALCreateTimestampVector
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:47
XLALListIFOsInCatalog
LALStringVector * XLALListIFOsInCatalog(const SFTCatalog *catalog)
Return a sorted string vector listing the unique IFOs in the given catalog.
Definition: SFTcatalog.c:562
XLALSFTdataFind
SFTCatalog * XLALSFTdataFind(const CHAR *file_pattern, const SFTConstraints *constraints)
Find the list of SFTs matching the file_pattern and satisfying the given constraints,...
Definition: SFTcatalog.c:71
XLALCWGPSinRange
int XLALCWGPSinRange(const LIGOTimeGPS gps, const LIGOTimeGPS *minGPS, const LIGOTimeGPS *maxGPS)
Defines the official CW convention for whether a GPS time is 'within' a given range,...
Definition: SFTtypes.c:52
XLALDestroyTimestampVector
void XLALDestroyTimestampVector(LIGOTimeGPSVector *vect)
De-allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:69
SSBprecisionChoices
const UserChoices SSBprecisionChoices
Static array of all SSBprecision choices, for use by the UserInput module parsing routines.
Definition: SSBtimes.c:41
XLALSegListClear
int XLALSegListClear(LALSegList *seglist)
XLALDestroyStringVector
void XLALDestroyStringVector(LALStringVector *vect)
XLALComputeOrbitalDerivatives
vect3Dlist_t * XLALComputeOrbitalDerivatives(UINT4 maxorder, const LIGOTimeGPS *tGPS, const EphemerisData *edat)
Compute time-derivatives up to 'maxorder' of the Earths' orbital position vector .
Definition: UniversalDopplerMetric.c:2446
XLALDestroyVect3Dlist
void XLALDestroyVect3Dlist(vect3Dlist_t *list)
Definition: UniversalDopplerMetric.c:2548
XLALUserVarReadAllInput
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
XLALDestroyUserVars
void XLALDestroyUserVars(void)
XLALUserVarGetLog
void CHAR * XLALUserVarGetLog(UserVarLogFormat format)
XLALRegisterNamedUvar
#define XLALRegisterNamedUvar(cvar, name, type, option, category,...)
XLALUserVarWasSet
int XLALUserVarWasSet(const void *cvar)
XLALRegisterNamedUvarAuxData
#define XLALRegisterNamedUvarAuxData(cvar, name, type, cdata, option, category,...)
UVAR_LOGFMT_CFGFILE
UVAR_LOGFMT_CFGFILE
XLAL_ERROR_NULL
#define XLAL_ERROR_NULL(...)
xlalErrno
#define xlalErrno
XLAL_ERROR
#define XLAL_ERROR(...)
XLAL_CHECK
#define XLAL_CHECK(assertion,...)
XLAL_CHECK_MAIN
#define XLAL_CHECK_MAIN(assertion,...)
XLALPrintError
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLALPrintWarning
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFAULT
XLAL_EFUNC
XLAL_EFUNC
XLAL_EDOM
XLAL_EDOM
XLAL_ESYS
XLAL_ESYS
XLAL_EINVAL
XLAL_EINVAL
XLAL_EFAILED
XLAL_EFAILED
XLALGPSSetREAL8
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
XLALGPSAdd
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
XLALGPSGetREAL8
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
XLALGPSDiff
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)
count
long long count
Definition: hwinject.c:371
T
T
pos
pos
size
size
end
end
out
out
deltaF
int deltaF
lalpulsar.git_version.status
string status
Definition: git_version.py:32
lalpulsar.piecewise_model.pw_fstat.SFTCatalog
def SFTCatalog(tstart, Tdata, finputdata)
Definition: pw_fstat.py:151
lalpulsar_MakeSFTDAG.segList
segList
Definition: lalpulsar_MakeSFTDAG.py:797
lalpulsar_PiecewiseSearch.line
string line
Definition: lalpulsar_PiecewiseSearch.py:914
lalpulsar_PiecewiseSearch.Tsft
Tsft
Definition: lalpulsar_PiecewiseSearch.py:598
lalpulsar_knope_result_page.f0
f0
Definition: lalpulsar_knope_result_page.py:2463
lalpulsar_knope.default
default
Definition: lalpulsar_knope.py:52
lalpulsar_knope.fp
fp
Definition: lalpulsar_knope.py:137
make_frame_cache.start
start
Definition: make_frame_cache.py:49
test_computePSD.edat
edat
Definition: test_computePSD.py:65
test_computePSD.multiTS
multiTS
Definition: test_computePSD.py:49
test_computePSD.numSFTs
numSFTs
Definition: test_computePSD.py:68
uvar_skyRegion
CHAR * uvar_skyRegion
Definition: pulsar_crosscorr.c:85
uvar_dAlpha
REAL8 uvar_dAlpha
Definition: pulsar_crosscorr.c:61
uvar_refTime
REAL8 uvar_refTime
Definition: pulsar_crosscorr.c:64
uvar_dDelta
REAL8 uvar_dDelta
Definition: pulsar_crosscorr.c:61
uvar_ephemEarth
CHAR * uvar_ephemEarth
Earth ephemeris file to use.
Definition: pulsar_crosscorr.c:78
uvar_blocksRngMed
INT4 uvar_blocksRngMed
Definition: pulsar_crosscorr.c:55
uvar_ephemSun
CHAR * uvar_ephemSun
Sun ephemeris file to use.
Definition: pulsar_crosscorr.c:79
alpha
double alpha
numDetectors
size_t numDetectors
CoarseGrid
structure for storing coarse-grid points
Definition: HierarchSearchGCT.h:189
DopplerSkyScanInit
initialization-structure passed to InitDopplerSkyScan()
Definition: DopplerScan.h:134
DopplerSkyScanInit::ephemeris
const EphemerisData * ephemeris
ephemeris-data for "exact" metric
Definition: DopplerScan.h:145
DopplerSkyScanInit::projectMetric
BOOLEAN projectMetric
project the metric orthogonal to Freq?
Definition: DopplerScan.h:143
DopplerSkyScanInit::numSkyPartitions
UINT4 numSkyPartitions
number of (roughly) equal partitions to split sky-grid into
Definition: DopplerScan.h:147
DopplerSkyScanInit::partitionIndex
UINT4 partitionIndex
index of requested sky-grid partition: in [0, numPartitions - 1]
Definition: DopplerScan.h:148
DopplerSkyScanInit::skyGridFile
const CHAR * skyGridFile
file containing a sky-grid (list of points) if GRID_FILE
Definition: DopplerScan.h:146
DopplerSkyScanInit::obsBegin
LIGOTimeGPS obsBegin
GPS start-time of time-series.
Definition: DopplerScan.h:141
DopplerSkyScanInit::dDelta
REAL8 dDelta
sky step-sizes for GRID_FLAT and GRID_ISOTROPIC
Definition: DopplerScan.h:139
DopplerSkyScanInit::Freq
REAL8 Freq
Frequency for which to build the skyGrid.
Definition: DopplerScan.h:136
DopplerSkyScanInit::gridType
DopplerGridType gridType
which type of skygrid to generate
Definition: DopplerScan.h:137
DopplerSkyScanInit::obsDuration
REAL8 obsDuration
length of time-series in seconds
Definition: DopplerScan.h:142
DopplerSkyScanInit::metricMismatch
REAL8 metricMismatch
for GRID_METRIC
Definition: DopplerScan.h:140
DopplerSkyScanInit::metricType
LALPulsarMetricType metricType
which metric to use if GRID_METRIC
Definition: DopplerScan.h:138
DopplerSkyScanInit::skyRegionString
CHAR * skyRegionString
sky-region to search: format polygon '(a1,d1), (a2,d2), ..'
Definition: DopplerScan.h:135
DopplerSkyScanInit::Detector
const LALDetector * Detector
Current detector.
Definition: DopplerScan.h:144
DopplerSkyScanState
this structure reflects the current state of a DopplerSkyScan
Definition: DopplerScan.h:152
EphemerisData
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
Definition: LALBarycenter.h:127
FineGrid::sumTwoFX
REAL4 * sumTwoFX
sum of per-IFO 2F-values, 2D array over frequencies and detectors (of length 'freqlengthAL*numDetecto...
Definition: HierarchSearchGCT.h:164
FineGrid::alpha
REAL8 alpha
right ascension
Definition: HierarchSearchGCT.h:156
FineGrid::nc
FINEGRID_NC_T * nc
number count (1D array over frequencies, of length 'length')
Definition: HierarchSearchGCT.h:165
FineGrid::maxTwoFXlIdx
UINT4 * maxTwoFXlIdx
segment index (zero based) of corresponding entry in maxTwoFXl
Definition: HierarchSearchGCT.h:169
FineGrid::delta
REAL8 delta
declination
Definition: HierarchSearchGCT.h:157
FineGrid::dfreq_fg
REAL8 dfreq_fg
fine-grid spacing in frequency
Definition: HierarchSearchGCT.h:155
FineGrid::numDetectors
UINT4 numDetectors
number of detectors for sumTwoFX array
Definition: HierarchSearchGCT.h:161
FineGrid::freqmin_fg
REAL8 freqmin_fg
fine-grid start in frequency
Definition: HierarchSearchGCT.h:154
FineGrid::freqlength
UINT4 freqlength
number of fine-grid points in frequency
Definition: HierarchSearchGCT.h:160
FineGrid::maxTwoFXl
REAL4 * maxTwoFXl
maximum of per-IFO 2F over segments, 2D array over frequencies and detectors (of length 'freqlengthAL...
Definition: HierarchSearchGCT.h:167
FineGrid::length
UINT4 length
length of multi-IFO stats vectors 'sumTwoF', 'nc' (currently 'length'= 'freqlength')
Definition: HierarchSearchGCT.h:159
FineGrid::freqlengthAL
UINT4 freqlengthAL
"aligned" number of fine-grid points in frequency: in blocks of 16 bytes, consistent with ALAlloc() [...
Definition: HierarchSearchGCT.h:163
FineGrid::refTime
LIGOTimeGPS refTime
reference time for candidates
Definition: HierarchSearchGCT.h:158
FineGrid::maxTwoFl
REAL4 * maxTwoFl
maximum of multi-IFO 2F over segments, 1D array over fine-grid frequencies (of length 'length')
Definition: HierarchSearchGCT.h:166
FineGrid::maxTwoFlIdx
UINT4 * maxTwoFlIdx
segment index (zero based) of corresponding entry in maxTwoFl
Definition: HierarchSearchGCT.h:168
FineGrid::sumTwoF
REAL4 * sumTwoF
sum of 2F-values, 1D array over fine-grid frequencies (of length 'length')
Definition: HierarchSearchGCT.h:162
FstatInputVector
A vector of XLALComputeFstat() input data structures, for e.g.
Definition: ComputeFstat.h:82
FstatInputVector::data
FstatInput ** data
Pointer to the data array.
Definition: ComputeFstat.h:87
FstatOptionalArgs
Struct of optional 'advanced level' and (potentially method-specific) arguments to be passed to the ...
Definition: ComputeFstat.h:137
FstatOptionalArgs::injectSources
PulsarParamsVector * injectSources
Vector of parameters of CW signals to simulate and inject.
Definition: ComputeFstat.h:143
FstatOptionalArgs::runningMedianWindow
UINT4 runningMedianWindow
If SFT noise weights are calculated from the SFTs, the running median window length to use.
Definition: ComputeFstat.h:141
FstatOptionalArgs::assumeSqrtSX
MultiNoiseFloor * assumeSqrtSX
Single-sided PSD values to be used for computing SFT noise weights instead of from a running median o...
Definition: ComputeFstat.h:145
FstatOptionalArgs::prevInput
FstatInput * prevInput
An FstatInput structure from a previous call to XLALCreateFstatInput(); may contain common workspace ...
Definition: ComputeFstat.h:146
FstatOptionalArgs::collectTiming
BOOLEAN collectTiming
a flag to turn on/off the collection of F-stat-method-specific timing-data
Definition: ComputeFstat.h:147
FstatOptionalArgs::Dterms
UINT4 Dterms
Number of Dirichlet kernel terms, used by some Demod methods; see FstatMethodType.
Definition: ComputeFstat.h:140
FstatOptionalArgs::FstatMethod
FstatMethodType FstatMethod
Method to use for computing the -statistic.
Definition: ComputeFstat.h:142
FstatOptionalArgs::SSBprec
SSBprecision SSBprec
Barycentric transformation precision.
Definition: ComputeFstat.h:139
FstatResults
XLALComputeFstat() computed results structure.
Definition: ComputeFstat.h:202
FstatResults::numDetectors
UINT4 numDetectors
Number of detectors over which the were computed.
Definition: ComputeFstat.h:221
FstatResults::twoFPerDet
REAL4 * twoFPerDet[PULSAR_MAX_DETECTORS]
If whatWasComputed & FSTATQ_2F_PER_DET is true, the values computed at numFreqBins frequencies space...
Definition: ComputeFstat.h:277
FstatResults::dFreq
REAL8 dFreq
Spacing in frequency between each computed -statistic.
Definition: ComputeFstat.h:214
FstatResults::numFreqBins
UINT4 numFreqBins
Number of frequencies at which the were computed.
Definition: ComputeFstat.h:217
FstatResults::twoF
REAL4 * twoF
If whatWasComputed & FSTATQ_2F is true, the multi-detector values computed at numFreqBins frequencie...
Definition: ComputeFstat.h:242
FstatResults::detectorNames
CHAR detectorNames[PULSAR_MAX_DETECTORS][3]
Names of detectors over which the were computed.
Definition: ComputeFstat.h:225
GCTtopOutputEntry
Definition: GCTtoplist.h:42
LALDetector::frDetector
LALFrDetector frDetector
LALFrDetector::prefix
CHAR prefix[3]
LALSeg::end
LIGOTimeGPS end
LALSeg::id
INT4 id
LALSeg::start
LIGOTimeGPS start
LALSegList::length
UINT4 length
LALSegList::segs
LALSeg * segs
LIGOTimeGPS::gpsNanoSeconds
INT4 gpsNanoSeconds
LIGOTimeGPSVector
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
LIGOTimeGPSVector::data
LIGOTimeGPS * data
array of timestamps
Definition: SFTfileIO.h:193
MultiLALDetector
array of detectors definitions 'LALDetector'
Definition: DetectorStates.h:108
MultiLALDetector::length
UINT4 length
number of detectors
Definition: DetectorStates.h:109
MultiLALDetector::sites
LALDetector sites[PULSAR_MAX_DETECTORS]
array of site information
Definition: DetectorStates.h:110
MultiLIGOTimeGPSVector
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
Definition: SFTfileIO.h:198
MultiNoiseFloor
array of detector-specific 'noise floors' (ie PSD values), assumed constant over the frequency-band o...
Definition: DetectorStates.h:115
PulsarDopplerParams
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
Definition: PulsarDataTypes.h:141
PulsarDopplerParams::fkdot
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ... ] where fkdot = d^kFreq/dt^k.
Definition: PulsarDataTypes.h:145
PulsarDopplerParams::Delta
REAL8 Delta
Sky position: DEC (latitude) in equatorial coords and radians.
Definition: PulsarDataTypes.h:144
PulsarDopplerParams::refTime
LIGOTimeGPS refTime
Reference time of pulsar parameters (in SSB!)
Definition: PulsarDataTypes.h:142
PulsarDopplerParams::Alpha
REAL8 Alpha
Sky position: RA (longitude) in equatorial coords and radians.
Definition: PulsarDataTypes.h:143
PulsarDopplerParams::asini
REAL8 asini
Binary: projected, normalized orbital semi-major axis (s).
Definition: PulsarDataTypes.h:146
PulsarParamsVector
Straightforward vector type of N PulsarParams structs.
Definition: CWMakeFakeData.h:56
PulsarSpinRange
Contains a "spin-range", ie spins and corresponding bands at a given (SSB) reference GPS-time .
Definition: PulsarDataTypes.h:134
PulsarSpinRange::refTime
LIGOTimeGPS refTime
SSB reference GPS-time at which spin-range is defined.
Definition: PulsarDataTypes.h:135
PulsarSpinRange::fkdotBand
PulsarSpins fkdotBand
Vector of spin-bands , MUST be same length as fkdot.
Definition: PulsarDataTypes.h:137
RecalcStatsParams
Type containing input arguments for XLALComputeExtraStatsForToplist()
Definition: RecalcToplistStats.h:66
SFTCatalog
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238
SFTCatalog::data
SFTDescriptor * data
array of data-entries describing matched SFTs
Definition: SFTfileIO.h:243
SFTCatalog::length
UINT4 length
number of SFTs in catalog
Definition: SFTfileIO.h:242
SFTCatalogSequence
sequence of SFT catalogs – for each segment
Definition: HierarchSearchGCT.h:125
SFTCatalogSequence::data
SFTCatalog * data
the catalogs
Definition: HierarchSearchGCT.h:127
SFTCatalogSequence::length
UINT4 length
the number of stacks
Definition: HierarchSearchGCT.h:126
SFTConstraints
'Constraints' for SFT-matching: which detector, within which time-stretch and which timestamps exactl...
Definition: SFTfileIO.h:212
SFTConstraints::maxStartTime
LIGOTimeGPS * maxStartTime
only include SFTs whose epoch is < maxStartTime
Definition: SFTfileIO.h:215
SFTConstraints::minStartTime
LIGOTimeGPS * minStartTime
only include SFTs whose epoch is >= minStartTime
Definition: SFTfileIO.h:214
SFTDescriptor
A 'descriptor' of an SFT: basically containing the header-info plus an opaque description of where ex...
Definition: SFTfileIO.h:226
SFTDescriptor::header
SFTtype header
SFT-header info.
Definition: SFTfileIO.h:228
SemiCoherentParams
parameters for the semicoherent stage
Definition: HierarchSearchGCT.h:132
SemiCoherentParams::acc
REAL8VectorSequence * acc
Earth orbital acceleration for each segment (new)
Definition: HierarchSearchGCT.h:137
SemiCoherentParams::pos
REAL8VectorSequence * pos
Earth orbital position for each segment.
Definition: HierarchSearchGCT.h:135
SemiCoherentParams::tsMid
LIGOTimeGPSVector * tsMid
timestamps of mid points of segments
Definition: HierarchSearchGCT.h:133
SemiCoherentParams::outBaseName
CHAR * outBaseName
file for writing output – if chosen
Definition: HierarchSearchGCT.h:138
SemiCoherentParams::refTime
LIGOTimeGPS refTime
reference time for f, fdot definition
Definition: HierarchSearchGCT.h:134
SemiCoherentParams::vel
REAL8VectorSequence * vel
Earth orbital velocity for each segment.
Definition: HierarchSearchGCT.h:136
SemiCoherentParams::extraBinsFstat
UINT4 extraBinsFstat
Extra bins required for Fstat calculation.
Definition: HierarchSearchGCT.h:140
UsefulStageVariables
useful variables for each hierarchical stage
Definition: HierarchSearchGCT.c:108
UsefulStageVariables::mismatch1
REAL8 mismatch1
'mismatch1' user-input needed here internally ...
Definition: HierarchSearchGCT.c:145
UsefulStageVariables::startTstack
LIGOTimeGPSVector * startTstack
timestamps vector for start time of each stack
Definition: HierarchSearchGCT.c:119
UsefulStageVariables::df2dot
REAL8 df2dot
coarse grid resolution in 2nd spindown
Definition: HierarchSearchGCT.c:134
UsefulStageVariables::spinRange_refTime
PulsarSpinRange spinRange_refTime
freq and fdot range at the reference time
Definition: HierarchSearchGCT.c:115
UsefulStageVariables::endTstack
LIGOTimeGPSVector * endTstack
timestamps vector for end time of each stack
Definition: HierarchSearchGCT.c:120
UsefulStageVariables::tObs
REAL8 tObs
tEndGPS - tStartGPS
Definition: HierarchSearchGCT.c:111
UsefulStageVariables::nf1dot
UINT4 nf1dot
number of 1st spindown Fstat bins
Definition: HierarchSearchGCT.c:138
UsefulStageVariables::binsFstatSearch
UINT4 binsFstatSearch
nominal number of Fstat frequency bins in search band
Definition: HierarchSearchGCT.c:137
UsefulStageVariables::Fmethod
FstatMethodType Fmethod
which Fstat-method/algorithm to use
Definition: HierarchSearchGCT.c:142
UsefulStageVariables::tStartGPS
LIGOTimeGPS tStartGPS
start and end time of stack
Definition: HierarchSearchGCT.c:110
UsefulStageVariables::df1dot
REAL8 df1dot
coarse grid resolution in spindown
Definition: HierarchSearchGCT.c:133
UsefulStageVariables::nStacks
UINT4 nStacks
number of stacks
Definition: HierarchSearchGCT.c:129
UsefulStageVariables::spinRange_startTime
PulsarSpinRange spinRange_startTime
freq and fdot range at start-time of observation
Definition: HierarchSearchGCT.c:113
UsefulStageVariables::spinRange_midTime
PulsarSpinRange spinRange_midTime
freq and fdot range at mid-time of observation
Definition: HierarchSearchGCT.c:116
UsefulStageVariables::NSegmentsInvX
REAL4 NSegmentsInvX[PULSAR_MAX_DETECTORS]
effective inverse number of segments per detector (needed for correct averaging in single-IFO F calcu...
Definition: HierarchSearchGCT.c:148
UsefulStageVariables::Fstat_in_vec_recalc
FstatInputVector * Fstat_in_vec_recalc
Recalculate the toplist: Vector of Fstat input data structures for XLALComputeFstat(),...
Definition: HierarchSearchGCT.c:150
UsefulStageVariables::tStack
REAL8 tStack
duration of stacks
Definition: HierarchSearchGCT.c:128
UsefulStageVariables::nf2dot
UINT4 nf2dot
number of 2nd spindown Fstat bins
Definition: HierarchSearchGCT.c:139
UsefulStageVariables::dFreqStack
REAL8 dFreqStack
frequency resolution of Fstat calculation
Definition: HierarchSearchGCT.c:132
UsefulStageVariables::detectorIDs
LALStringVector * detectorIDs
vector of detector IDs
Definition: HierarchSearchGCT.c:147
UsefulStageVariables::Fstat_in_vec
FstatInputVector * Fstat_in_vec
Original wide-parameter search: vector of Fstat input data structures for XLALComputeFstat(),...
Definition: HierarchSearchGCT.c:149
UsefulStageVariables::spinRange_endTime
PulsarSpinRange spinRange_endTime
freq and fdot range at end-time of observation
Definition: HierarchSearchGCT.c:114
UsefulStageVariables::refTime
REAL8 refTime
reference time for pulsar params
Definition: HierarchSearchGCT.c:112
UsefulStageVariables::maxStartTimeGPS
LIGOTimeGPS maxStartTimeGPS
all sft timestamps must be before this GPS time
Definition: HierarchSearchGCT.c:122
UsefulStageVariables::segmentList
LALSegList * segmentList
parsed segment list read from user-specified input file –segmentList
Definition: HierarchSearchGCT.c:130
UsefulStageVariables::Dterms
UINT4 Dterms
size of Dirichlet kernel for Fstat calculation
Definition: HierarchSearchGCT.c:124
UsefulStageVariables::recalcToplistStats
BOOLEAN recalcToplistStats
do additional analysis for all toplist candidates, output F, FXvector for postprocessing *‍/
Definition: HierarchSearchGCT.c:143
UsefulStageVariables::edat
EphemerisData * edat
ephemeris data for XLALBarycenter
Definition: HierarchSearchGCT.c:117
UsefulStageVariables::minStartTimeGPS
LIGOTimeGPS minStartTimeGPS
all sft data must be after this time
Definition: HierarchSearchGCT.c:121
UsefulStageVariables::nf3dot
UINT4 nf3dot
number of 3rd spindown Fstat bins
Definition: HierarchSearchGCT.c:140
UsefulStageVariables::extraBinsFstat
UINT4 extraBinsFstat
Extra Fstat frequency bins required to cover residual spindowns.
Definition: HierarchSearchGCT.c:136
UsefulStageVariables::FmethodRecalc
FstatMethodType FmethodRecalc
which Fstat-method/algorithm to use for the recalc step
Definition: HierarchSearchGCT.c:144
UsefulStageVariables::sftbasename
CHAR * sftbasename
filename pattern for sfts
Definition: HierarchSearchGCT.c:109
UsefulStageVariables::nSFTs
UINT4 nSFTs
total number of SFTs
Definition: HierarchSearchGCT.c:146
UsefulStageVariables::midTstack
LIGOTimeGPSVector * midTstack
timestamps vector for mid time of each stack
Definition: HierarchSearchGCT.c:118
UsefulStageVariables::df3dot
REAL8 df3dot
coarse grid resolution in 3rd spindown
Definition: HierarchSearchGCT.c:135
UsefulStageVariables::SSBprec
int SSBprec
SSB transform precision.
Definition: HierarchSearchGCT.c:141
UsefulStageVariables::DtermsRecalc
UINT4 DtermsRecalc
Recalc: size of Dirichlet kernel for Fstat calculation.
Definition: HierarchSearchGCT.c:125
UsefulStageVariables::collectFstatTiming
BOOLEAN collectFstatTiming
flag whether to collect and output F-stat timing info
Definition: HierarchSearchGCT.c:152
UsefulStageVariables::assumeSqrtSX
LALStringVector * assumeSqrtSX
Assume stationary Gaussian noise with detector noise-floors sqrt{SX}".
Definition: HierarchSearchGCT.c:126
UsefulStageVariables::BSGLsetup
BSGLSetup * BSGLsetup
pre-computed setup for line-robust statistic BSGL
Definition: HierarchSearchGCT.c:131
UsefulStageVariables::blocksRngMed
UINT4 blocksRngMed
blocksize for running median noise floor estimation
Definition: HierarchSearchGCT.c:123
UsefulStageVariables::injectionSources
PulsarParamsVector * injectionSources
Source parameters to inject: comma-separated list of file-patterns and/or direct config-strings ('{....
Definition: HierarchSearchGCT.c:151
timingInfo_t
Struct holding various timing measurements and relevant search parameters.
Definition: ComputeFstatistic_v2.c:119
timingInfo_t::time_Other
REAL8 time_Other
Definition: HierarchSearchGCT.c:181
timingInfo_t::Nsft
UINT4 Nsft
total number of SFTs
Definition: HierarchSearchGCT.c:164
timingInfo_t::Ninc
REAL8 Ninc
number of fine-grid templates ('incoherent')
Definition: HierarchSearchGCT.c:169
timingInfo_t::Ncoh
REAL8 Ncoh
number of coarse-grid Fstat templates ('coherent')
Definition: HierarchSearchGCT.c:168
timingInfo_t::tau_Recalc
REAL8 tau_Recalc
Definition: HierarchSearchGCT.c:180
timingInfo_t::Tcoh
UINT4 Tcoh
length of coherent segments in seconds
Definition: HierarchSearchGCT.c:163
timingInfo_t::tau_SumF
REAL8 tau_SumF
Definition: HierarchSearchGCT.c:178
timingInfo_t::Ncand
UINT4 Ncand
length of toplists
Definition: HierarchSearchGCT.c:165
timingInfo_t::NFreqCo
UINT4 NFreqCo
total number of frequency bins computed in coarse grid (including sidebands!)
Definition: HierarchSearchGCT.c:167
timingInfo_t::tau_Fstat
REAL8 tau_Fstat
Definition: HierarchSearchGCT.c:177
timingInfo_t::Ndet
UINT4 Ndet
number of detectors
Definition: HierarchSearchGCT.c:162
timingInfo_t::FstatMethodStr
const char * FstatMethodStr
Fstat-method used.
Definition: HierarchSearchGCT.c:171
timingInfo_t::RecalcMethodStr
const char * RecalcMethodStr
Fstat-method used.
Definition: HierarchSearchGCT.c:172
timingInfo_t::tau_Bayes
REAL8 tau_Bayes
Definition: HierarchSearchGCT.c:179
timingInfo_t::Nseg
UINT4 Nseg
number of semi-coherent segments
Definition: HierarchSearchGCT.c:161
toplist_t
Definition: HeapToplist.h:35
toplist_t::elems
size_t elems
Definition: HeapToplist.h:37
vect3Dlist_t
variable-length list of 3D vectors
Definition: UniversalDopplerMetric.h:77
vect3Dlist_t::data
vect3D_t * data
array of 3D vectors
Definition: UniversalDopplerMetric.h:82