HierarchSearchGCT.c 146 KB
Newer Older
1
/*
2
 *  Copyright (C) 2011-2013 Karl Wette.
3
 *  Copyright (C) 2009-2010 Holger Pletsch.
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
 *
 *  Based on HierarchicalSearch.c by
 *  Copyright (C) 2005-2008 Badri Krishnan, Alicia Sintes, Bernd Machenschalk.
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
19 20
 *  along with with program; see the file COPYING. If not, write to the
 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
21
 *  MA  02111-1307  USA
22
 *
23 24 25
 */

/*********************************************************************************/
26 27 28 29 30
/**
 * \defgroup lalapps_pulsar_GCT GCT Search Application
 * \ingroup lalapps_pulsar_Apps
 */

31 32
/**
 * \author Holger Pletsch
33
 * \file
34
 * \ingroup lalapps_pulsar_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 44
#include <lal/LineRobustStats.h>
#include <RecalcToplistStats.h>
45

Holger Pletsch's avatar
Holger Pletsch committed
46
#include "HierarchSearchGCT.h"
47

48 49 50 51
#ifdef GC_SSE2_OPT
#include <gc_hotloop_sse2.h>
#else
#define ALRealloc LALRealloc
52
#define ALFree LALFree
53 54
#endif

55
/* ---------- Defines -------------------- */
Holger Pletsch's avatar
Holger Pletsch committed
56
/* #define DIAGNOSISMODE 1 */
57
#define NUDGE	10*LAL_REAL8_EPS
58

59 60 61
#define TRUE (1==1)
#define FALSE (1==0)

62 63 64 65 66 67
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

68 69 70
/* Hooks for Einstein@Home / BOINC
   These are defined to do nothing special in the standalone case
   and will be set in boinc_extras.h if EAH_BOINC is set
71
*/
72 73
#ifdef EAH_BOINC
#include "hs_boinc_extras.h"
74 75 76
// reserves 1% of progress for the last (toplist recaclculation) step
#define SHOW_PROGRESS_RESERVE(rac,dec,count,total,freq,fband)\
        SHOW_PROGRESS(rac, dec, (count) - 0.01 * (total), total, freq, fband)
77
#else
Bernd Machenschalk's avatar
Bernd Machenschalk committed
78 79
#define GET_GCT_CHECKPOINT read_gct_checkpoint // (cptname, semiCohToplist, NULL, &count)
#define SET_GCT_CHECKPOINT write_gct_checkpoint
80
#define SHOW_PROGRESS(rac,dec,skyGridCounter,tpl_total,freq,fband)
81
#define SHOW_PROGRESS_RESERVE(rac,dec,count,total,freq,fband)
82
#define MAIN  main
83 84
char**global_argv;
int global_argc;
85
#endif /* EAH_BOINC */
86

87 88
#define BLOCKSRNGMED    101     /**< Default running median window size */
#define FSTART          100.0	/**< Default Start search frequency */
89
#define FBAND           0.0  /**< Default search band */
90 91
#define FDOT            0.0       /**< Default value of first spindown */
#define DFDOT           0.0       /**< Default range of first spindown parameter */
92
#define F2DOT           0.0       /**< Default value of second spindown */
93
#define F3DOT           0.0       /**< Default value of third spindown */
94
#define DF2DOT          0.0       /**< Default range of second spindown parameter */
95
#define DF3DOT          0.0       /**< Default range of third spindown parameter */
96
#define SKYREGION       "allsky" /**< default sky region to search over -- just a single point*/
97
#define DTERMS          8    /**< Default number of dirichlet kernel terms for calculating Fstat */
98 99

/**< Default number of dirichlet kernel terms for calculating Fstat */
100 101 102 103 104
#define MISMATCH        0.3       /**< Default for metric grid maximal mismatch value */
#define DALPHA          0.001   /**< Default resolution for isotropic or flat grids */
#define DDELTA          0.001   /**< Default resolution for isotropic or flat grids */
#define FSTATTHRESHOLD  2.6	/**< Default threshold on Fstatistic for peak selection */
#define NCAND1          10      /**< Default number of candidates to be followed up from first stage */
105 106
#define FNAMEOUT        "./HS_GCT.out"  /**< Default output file basename */
#ifndef LAL_INT4_MAX
107
#define LAL_INT4_MAX    2147483647
108 109 110 111
#endif

#define BLOCKSIZE_REALLOC 50

112 113
#define Vorb_GCT   = 2.9785e04;
#define Vspin_GCT  = 465.10;
Holger Pletsch's avatar
Holger Pletsch committed
114 115
#define REARTH_GCT = 6.378140e06;
#define C_GCT      = 299792458;
116 117 118 119 120

/* ---------- Macros -------------------- */
#define HSMAX(x,y) ( (x) > (y) ? (x) : (y) )
#define HSMIN(x,y) ( (x) < (y) ? (x) : (y) )

121 122
#define GETTIME() (uvar_outputTiming ? XLALGetCPUTime() : 0)
//#define GETTIME() (uvar_outputTiming ? XLALGetTimeOfDay() : 0)
123 124 125 126 127 128 129 130 131 132 133 134

/* ---------- Exported types ---------- */
/** useful variables for each hierarchical stage */
typedef struct {
  CHAR  *sftbasename;    /**< filename pattern for sfts */
  LIGOTimeGPS tStartGPS; /**< start and end time of stack */
  REAL8 tObs;            /**< tEndGPS - tStartGPS */
  REAL8 refTime;         /**< reference time for pulsar params */
  PulsarSpinRange spinRange_startTime; /**< freq and fdot range at start-time of observation */
  PulsarSpinRange spinRange_endTime;   /**< freq and fdot range at end-time of observation */
  PulsarSpinRange spinRange_refTime;   /**< freq and fdot range at the reference time */
  PulsarSpinRange spinRange_midTime;   /**< freq and fdot range at mid-time of observation */
135
  EphemerisData *edat;             /**< ephemeris data for XLALBarycenter */
136 137 138
  LIGOTimeGPSVector *midTstack;    /**< timestamps vector for mid time of each stack */
  LIGOTimeGPSVector *startTstack;  /**< timestamps vector for start time of each stack */
  LIGOTimeGPSVector *endTstack;    /**< timestamps vector for end time of each stack */
139
  LIGOTimeGPS minStartTimeGPS;     /**< all sft data must be after this time */
140
  LIGOTimeGPS maxStartTimeGPS;       /**< all sft timestamps must be before this GPS time */
141 142
  UINT4 blocksRngMed;              /**< blocksize for running median noise floor estimation */
  UINT4 Dterms;                    /**< size of Dirichlet kernel for Fstat calculation */
143 144
  LALStringVector* assumeSqrtSX;   /**< Assume stationary Gaussian noise with detector noise-floors sqrt{SX}" */
  BOOLEAN SignalOnly;              /**< DEPRECATED: ALTERNATIVE switch to assume Sh=1 instead of estimating noise-floors from SFTs */
145
  /* parameters describing the coherent data-segments */
146 147
  REAL8 tStack;                    /**< duration of stacks */
  UINT4 nStacks;                   /**< number of stacks */
148
  LALSegList *segmentList;         /**< parsed segment list read from user-specified input file --segmentList */
149
  BSGLSetup *BSGLsetup;           /**< pre-computed setup for line-robust statistic BSGL */
150 151
  REAL8 dFreqStack;                /**< frequency resolution of Fstat calculation */
  REAL8 df1dot;                    /**< coarse grid resolution in spindown */
152
  REAL8 df2dot;                    /**< coarse grid resolution in 2nd spindown */
153
  REAL8 df3dot;                    /**< coarse grid resolution in 3rd spindown */
154 155 156 157 158
  UINT4 extraBinsFstat;            /**< Extra Fstat frequency bins required to cover residual spindowns */
  UINT4 binsFstatSearch;	   /**< nominal number of Fstat frequency bins in search band */
  UINT4 nf1dot;			/**< number of 1st spindown Fstat bins */
  UINT4 nf2dot;			/**< number of 2nd spindown Fstat bins */
  UINT4 nf3dot;			/**< number of 3rd spindown Fstat bins */
159
  SSBprecision SSBprec;            /**< SSB transform precision */
160
  FstatMethodType Fmethod;         //!< which Fstat-method/algorithm to use
161
  BOOLEAN recalcToplistStats;	   //!< do additional analysis for all toplist candidates, output F, FXvector for postprocessing */
162
  FstatMethodType FmethodRecalc;   //!< which Fstat-method/algorithm to use for the recalc step
163 164 165 166
  REAL8 mismatch1;                 /**< 'mismatch1' user-input needed here internally ... */
  UINT4 nSFTs;                     /**< total number of SFTs */
  LALStringVector *detectorIDs;    /**< vector of detector IDs */
  REAL4 NSegmentsInvX[PULSAR_MAX_DETECTORS]; /**< effective inverse number of segments per detector (needed for correct averaging in single-IFO F calculation) */
167 168
  FstatInputVector* Fstat_in_vec;	/**< Original wide-parameter search: vector of Fstat input data structures for XLALComputeFstat(), one per stack */
  FstatInputVector* Fstat_in_vec_recalc; /**< Recalculate the toplist: Vector of Fstat input data structures for XLALComputeFstat(), one per stack */
169
  FILE *timingDetailsFP;	// file pointer to write detailed coherent-timing info into
170
  PulsarParamsVector *injectionSources; ///< Source parameters to inject: comma-separated list of file-patterns and/or direct config-strings ('{...}')
171 172
} UsefulStageVariables;

Holger Pletsch's avatar
Holger Pletsch committed
173

174 175 176 177 178 179 180 181 182 183
/**
 * Struct holding various timing measurements and relevant search parameters.
 * This is used to fit timing-models with measured times to predict search run-times
 */
typedef struct
{
  UINT4 Nseg;			///< number of semi-coherent segments
  UINT4 Ndet;			///< number of detectors
  UINT4 Tcoh;			///< length of coherent segments in seconds
  UINT4 Nsft;			///< total number of SFTs
184
  UINT4 Ncand;			///< length of toplists
185 186

  UINT4 NFreqCo;		///< total number of frequency bins computed in coarse grid (including sidebands!)
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201
  REAL8 Ncoh;			///< number of coarse-grid Fstat templates ('coherent')
  REAL8 Ninc;			///< number of fine-grid templates ('incoherent')

  const char* FstatMethodStr;	///< Fstat-method used
  const char* RecalcMethodStr;	///< Fstat-method used

  // ----------
  // extended timing model:
  // runtime = Nseg * Ndet * Ncoh * tauF + Nseg * Ninc * tauSumFine + Ninc * tauExtraStats + Ntop * tauRecalc + time_Other
  REAL8 tau_Fstat;		//< time to compute F-stat for one coarse-grid template, one detector, one segment
  REAL8 tau_SumF;		//< time to sum F-stat values of one segment for one fine-grid point
  REAL8 tau_Bayes;		//< time to compute final Bayes-factor statistics for one fine-grid point
  REAL8 tau_Recalc;		//< time to recalc one toplist-entry
  REAL8 time_Other;		//< all non-scaling leftovers in overall timeing
  // ----------
202 203 204 205

} timingInfo_t;


206
/* ------------------------ Functions -------------------------------- */
207
void SetUpSFTs( LALStatus *status, UsefulStageVariables *in );
208
void PrintFstatVec( LALStatus *status, FstatResults *in, FILE *fp, PulsarDopplerParams *thisPoint,
209
                    LIGOTimeGPS refTime, INT4 stackIndex);
210 211
void PrintCatalogInfo( LALStatus *status, const SFTCatalog *catalog, FILE *fp );
void PrintStackInfo( LALStatus *status, const SFTCatalogSequence *catalogSeq, FILE *fp );
212
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 );
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227

void UpdateSemiCohToplistsOptimTriple ( 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
                             );

228
void GetSegsPosVelAccEarthOrb( LALStatus *status, REAL8VectorSequence **posSeg,
229 230
                               REAL8VectorSequence **velSeg, REAL8VectorSequence **accSeg,
                               UsefulStageVariables *usefulparams );
231
static inline INT4 ComputeU1idx( REAL8 freq_event, REAL8 f1dot_event, REAL8 A1, REAL8 B1, REAL8 U1start, REAL8 U1winInv );
232
void ComputeU2idx( REAL8 freq_event, REAL8 f1dot_event, REAL8 A2, REAL8 B2, REAL8 U2start, REAL8 U2winInv,
233
                   INT4 *U2idx);
234
int compareCoarseGridUindex( const void *a, const void *b );
235 236
int compareFineGridNC( const void *a,const void *b );
int compareFineGridsumTwoF( const void *a,const void *b );
237

238
SFTCatalogSequence *XLALSetUpStacksFromSegmentList ( const SFTCatalog *catalog, const LALSegList *segList );
239

240 241 242 243
int XLALExtrapolateToplistPulsarSpins ( toplist_t *list,
					const LIGOTimeGPS usefulParamsRefTime,
					const LIGOTimeGPS finegridRefTime);

244
static int write_TimingInfo ( const CHAR *fname, const timingInfo_t *ti );
245

246 247
/* ---------- Global variables -------------------- */
LALStatus *global_status; /* a global pointer to MAIN()s head of the LALStatus structure */
248
char *global_column_headings_stringp;
249 250 251 252 253 254

/* ###################################  MAIN  ################################### */

int MAIN( int argc, char *argv[]) {
  LALStatus status = blank_status;

255
  /* temp loop variables: generally k loops over segments and j over SFTs in a stack */
256 257
  UINT4 k;
  UINT4 skyGridCounter; /* coarse sky position counter */
258
  UINT4 f1dotGridCounter; /* coarse f1dot position counter */
259 260

  /* GPS timestamp vectors */
261 262 263 264
  LIGOTimeGPSVector *midTstack = NULL;
  LIGOTimeGPSVector *startTstack = NULL;
  LIGOTimeGPSVector *endTstack = NULL;

265
  /* General GPS times */
266 267
  LIGOTimeGPS XLAL_INIT_DECL(refTimeGPS);
  LIGOTimeGPS XLAL_INIT_DECL(tMidGPS);
268

Holger Pletsch's avatar
Holger Pletsch committed
269
  /* GPS time used for each segment's midpoint */
270
  LIGOTimeGPS XLAL_INIT_DECL(midTstackGPS);
Holger Pletsch's avatar
Holger Pletsch committed
271
  REAL8 timeDiffSeg; /* Difference to tMidGPS (midpoint) of Tobs */
272

273
  /* pos, vel, acc at midpoint of segments */
274 275
  REAL8VectorSequence *posStack = NULL;
  REAL8VectorSequence *velStack = NULL;
276 277
  REAL8VectorSequence *accStack = NULL;

278
  /* duration of each segment */
279 280
  REAL8 tStack;

Holger Pletsch's avatar
Holger Pletsch committed
281
  /* number of segments */
Holger Pletsch's avatar
Holger Pletsch committed
282
  UINT4 nStacks;
283

284 285
  /* Total observation time */
  REAL8 tObs;
286

287
  /* SFT related stuff */
288
  static LIGOTimeGPS minStartTimeGPS, maxStartTimeGPS;
289

290
  /* some useful variables for each stage */
291
  UsefulStageVariables XLAL_INIT_DECL(usefulParams);
292

293
  /* F-statistic computation related stuff */
294 295
  FstatResults* Fstat_res = NULL;			// Pointer to Fstat results structure, will be allocated by XLALComputeFstat()
  FstatQuantities Fstat_what = FSTATQ_2F;		// Quantities to be computed by XLALComputeFstat()
296
  UINT4 binsFstat1;
297

298 299
  /* Semicoherent variables */
  static SemiCoherentParams semiCohPar;
300

301
  /* coarse grid */
302
  CoarseGrid XLAL_INIT_DECL(coarsegrid);
303 304
  REAL8 dFreqStack; /* frequency resolution of Fstat calculation */
  REAL8 df1dot;  /* coarse grid resolution in spindown */
305
  UINT4 ifdot;  /* counter for coarse-grid spindown values */
306 307
  REAL8 df2dot;  /* coarse grid resolution in 2nd spindown */
  UINT4 if2dot;  /* counter for coarse-grid 2nd spindown values */
308 309
  REAL8 df3dot;  /* coarse grid resolution in 3rd spindown */
  UINT4 if3dot;  /* counter for coarse-grid 3rd spindown values */
310

311 312
  /* fine grid */
  FineGrid finegrid;
313
  UINT4 nf1dots_fg = 1; /* number of frequency and spindown values */
314
  REAL8 gammaRefine, sigmasq;  /* refinement factor and variance */
315 316
  UINT4 nf2dots_fg=1;          /* number of second spindown values */
  REAL8 gamma2Refine, sigma4;  /* 2nd spindown refinement factor and 4th moment */
317 318 319
  UINT4 nf3dots_fg=1;          /* number of third spindown values */
  REAL8 gamma3Refine=1;  /* 3rd spindown refinement */
  
320
  /* GCT helper variables */
321
  UINT4 if1dot_fg, if2dot_fg, if3dot_fg;
322
  UINT4 ifreq;
323
  INT4  U1idx;
324
  REAL8 myf0, freq_event, f1dot_event;
325
  REAL8 dfreq_fg, df1dot_fg, freqmin_fg, f1dotmin_fg, freqband_fg;
326
  REAL8 df2dot_fg, f2dotmin_fg;
327
  REAL8 df3dot_fg, f3dotmin_fg;
Holger Pletsch's avatar
Holger Pletsch committed
328
  REAL8 u1start, u1win, u1winInv;
329
  REAL8 freq_fg, f1dot_fg, f2dot_fg, f3dot_fg, f1dot_event_fg;
330 331 332
  REAL8 A1, B1;
  // currently unused: REAL8 A2;
  REAL8 B2; /* GCT helper variables for faster calculation of u1 or u2 */
333 334
  REAL8 pos[3];
  REAL8 vel[3];
335
  // currently unused: REAL8 acc[3];
336 337
  REAL8 cosAlpha, sinAlpha, cosDelta, sinDelta;
  REAL8 nvec[3]; /* unit vector pointing to sky position */
338

339
  /* These vars are currently not used, but eventually in the future.
340 341 342
     INT4  U2idx, NumU2idx;
     REAL8 myf0max, u2start, u2end;
     REAL8 u2win, u2winInv;
343
  */
344

345
  /* fstat candidate structure for candidate toplist*/
346
  toplist_t *semiCohToplist=NULL;
347 348
  toplist_t *semiCohToplist2=NULL;	// only used for SORTBY_DUAL_F_BSGL or SORTBY_TRIPLE_BStSGLtL
  toplist_t *semiCohToplist3=NULL;	// only used for SORTBY_TRIPLE_BStSGLtL
349 350 351

  /* template and grid variables */
  static DopplerSkyScanInit scanInit;   /* init-structure for DopperScanner */
352
  DopplerSkyScanState XLAL_INIT_DECL(thisScan); /* current state of the Doppler-scan */
353
  static PulsarDopplerParams dopplerpos;               /* current search-parameters */
354
  static PulsarDopplerParams thisPoint;
355
  UINT4 oldcg=0, oldfg=0;
356 357 358 359 360 361 362

  /* temporary storage for spinrange vector */
  static PulsarSpinRange spinRange_Temp;

  /* variables for logging */
  CHAR *fnamelog=NULL;
  FILE *fpLog=NULL;
363
  CHAR *logstr=NULL;
364 365 366 367 368

  /* output candidate files and file pointers */
  CHAR *fnameSemiCohCand=NULL;
  CHAR *fnameFstatVec1=NULL;
  FILE *fpFstat1=NULL;
369

370
  /* checkpoint filename */
371
  CHAR *uvar_fnameChkPoint = NULL;
372

373
  /* user variables */
374
  BOOLEAN uvar_log = FALSE;     /* logging done if true */
375

376
  BOOLEAN uvar_printCand1 = FALSE;      /* if 1st stage candidates are to be printed */
377
  BOOLEAN uvar_printFstat1 = FALSE;
378
  BOOLEAN uvar_semiCohToplist = TRUE; /* if overall first stage candidates are to be output */
379 380 381

  LALStringVector* uvar_assumeSqrtSX = NULL;    /* Assume stationary Gaussian noise with detector noise-floors sqrt{SX}" */
  BOOLEAN uvar_SignalOnly = FALSE;              /* DEPRECATED: ALTERNATIVE switch to assume Sh=1 instead of estimating noise-floors from SFTs */
382 383

  BOOLEAN uvar_recalcToplistStats = FALSE; 	/* Do additional analysis for all toplist candidates, output F, FXvector for postprocessing */
384
  BOOLEAN uvar_loudestSegOutput = FALSE; 	/* output extra info about loudest segment; requires recalcToplistStats */
385 386

  // ----- Line robust stats parameters ----------
387 388
  BOOLEAN uvar_computeBSGL = FALSE;          	/* In Fstat loop, compute line-robust statistic (BSGL=log10BSGL) using single-IFO F-stats */
  BOOLEAN uvar_BSGLlogcorr = FALSE;		/* compute log-correction in line-robust statistic BSGL (slower) or not (faster) */
389
  REAL8   uvar_Fstar0sc = 0.0;			/* (semi-coherent) BSGL transition-scale parameter 'Fstar0sc=Nseg*Fstar0coh', see documentation for XLALCreateBSGLSetup() for details */
390
  LALStringVector *uvar_oLGX = NULL;       	/* prior per-detector line-vs-Gauss odds ratios 'oLGX', see XLALCreateBSGLSetup() for details */
391
  BOOLEAN uvar_getMaxFperSeg = FALSE;          	/* In Fstat loop, compute maximum F and FX over segments */
392 393
  // --------------------------------------------

394
  REAL8 uvar_dAlpha = DALPHA;   /* resolution for flat or isotropic grids -- coarse grid*/
395
  REAL8 uvar_dDelta = DDELTA;
396
  REAL8 uvar_f1dot = FDOT;      /* first spindown value */
397 398 399
  REAL8 uvar_f1dotBand = DFDOT; /* range of first spindown parameter */
  REAL8 uvar_f2dot = F2DOT;     /* second spindown value */
  REAL8 uvar_f2dotBand = DF2DOT; /* range of second spindown parameter */
400 401
  REAL8 uvar_f3dot = F3DOT;     /* second spindown value */
  REAL8 uvar_f3dotBand = DF3DOT; /* range of second spindown parameter */
402 403 404
  REAL8 uvar_Freq = FSTART;
  REAL8 uvar_FreqBand = FBAND;

405
  REAL8 uvar_dFreq = 0;
406
  REAL8 uvar_df1dot = 0; /* coarse grid frequency and spindown resolution */
407
  REAL8 uvar_df2dot = 0; /* coarse grid second spindown resolution */
408
  REAL8 uvar_df3dot = 0; /* coarse grid third spindown resolution */
409 410 411 412 413

  REAL8 uvar_ThrF = FSTATTHRESHOLD; /* threshold of Fstat to select peaks */
  REAL8 uvar_mismatch1 = MISMATCH; /* metric mismatch for first stage coarse grid */

  REAL8 uvar_minStartTime1 = 0;
414
  REAL8 uvar_maxStartTime1 = LAL_INT4_MAX;
415 416 417 418 419

  REAL8 uvar_refTime = 0;
  INT4 uvar_nCand1 = NCAND1; /* number of candidates to be followed up from first stage */

  INT4 uvar_blocksRngMed = BLOCKSRNGMED;
420 421 422

  REAL8 uvar_tStack = 0;
  INT4  uvar_nStacksMax = 1;
423
  CHAR *uvar_segmentList = NULL;	/**< ALTERNATIVE: file containing a pre-computed segment list of tuples (startGPS endGPS duration[h] NumSFTs) */
424

425
  INT4 uvar_Dterms = DTERMS;
426
  INT4 uvar_SSBprecision = SSBPREC_RELATIVISTIC;
427
  INT4 uvar_gammaRefine = 1;
428
  INT4 uvar_gamma2Refine = 1;
429 430 431 432
  INT4 uvar_metricType1 = LAL_PMETRIC_COH_PTOLE_ANALYTIC;
  INT4 uvar_gridType1 = GRID_METRIC;
  INT4 uvar_skyPointIndex = -1;

433 434
  CHAR *uvar_ephemEarth;	/**< Earth ephemeris file to use */
  CHAR *uvar_ephemSun;		/**< Sun ephemeris file to use */
435 436 437 438 439 440 441

  CHAR *uvar_skyRegion = NULL;
  CHAR *uvar_fnameout = NULL;
  CHAR *uvar_DataFiles1 = NULL;
  CHAR *uvar_skyGridFile=NULL;
  INT4 uvar_numSkyPartitions = 0;
  INT4 uvar_partitionIndex = 0;
442
  INT4 uvar_SortToplist = 0;
443
  BOOLEAN uvar_version = 0;
444

445
  CHAR *uvar_outputTiming = NULL;
446
  CHAR *uvar_outputTimingDetails = NULL;
447

448
  CHAR *uvar_FstatMethod = XLALStringDuplicate("DemodBest");
449
  CHAR *uvar_FstatMethodRecalc = XLALStringDuplicate("DemodBest");
450

451 452
  timingInfo_t XLAL_INIT_DECL(timing);

453 454 455

  LALStringVector *uvar_injectionSources = NULL;

456 457 458 459 460 461 462
  // timing values
  REAL8 tic_RecalcToplist, time_RecalcToplist = 0;
  REAL8 tic_Fstat, time_Fstat = 0;
  REAL8 tic_SumFine, time_SumFine = 0;
  REAL8 tic_ExtraStats, time_ExtraStats = 0;
  REAL8 tic_Start, time_Total = 0;

463 464
  global_status = &status;

465 466 467 468
#ifndef EAH_BOINC
  global_argv = argv;
  global_argc = argc;
#endif
469 470 471

#ifdef EAH_LALDEBUGLEVEL
#endif
472

473 474
  uvar_ephemEarth = XLALStringDuplicate("earth00-19-DE405.dat.gz");
  uvar_ephemSun = XLALStringDuplicate("sun00-19-DE405.dat.gz");
475 476 477 478 479 480 481 482 483 484 485 486 487 488 489

  uvar_skyRegion = LALCalloc( strlen(SKYREGION) + 1, sizeof(CHAR) );
  strcpy(uvar_skyRegion, SKYREGION);

  uvar_fnameout = LALCalloc( strlen(FNAMEOUT) + 1, sizeof(CHAR) );
  strcpy(uvar_fnameout, FNAMEOUT);

  /* set LAL error-handler */
#ifdef EAH_BOINC
  lal_errhandler = BOINC_LAL_ErrHand;
#else
  lal_errhandler = LAL_ERR_EXIT;
#endif

  /* register user input variables */
490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_log,                 "log",                 BOOLEAN,      0,   OPTIONAL,   "Write log file") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_semiCohToplist,      "semiCohToplist",      BOOLEAN,      0,   OPTIONAL,   "Print toplist of semicoherent candidates" ) == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_DataFiles1,          "DataFiles1",          STRING,       0,   REQUIRED,   "1st SFT file pattern") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_skyRegion,           "skyRegion",           STRING,       0,   OPTIONAL,   "sky-region polygon (or 'allsky')") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_numSkyPartitions,    "numSkyPartitions",    INT4,         0,   OPTIONAL,   "No. of (equi-)partitions to split skygrid into") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_partitionIndex,      "partitionIndex",      INT4,         0,   OPTIONAL,   "Index [0,numSkyPartitions-1] of sky-partition to generate") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_skyGridFile,         "skyGridFile",         STRING,       0,   OPTIONAL,   "sky-grid file") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dAlpha,              "dAlpha",              REAL8,        0,   OPTIONAL,   "Resolution for flat or isotropic coarse grid") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dDelta,              "dDelta",              REAL8,        0,   OPTIONAL,   "Resolution for flat or isotropic coarse grid") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_Freq,                "Freq",                REAL8,        'f', OPTIONAL,   "Start search frequency") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dFreq,               "dFreq",               REAL8,        0,   OPTIONAL,   "Frequency resolution (required if nonzero FreqBand)") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_FreqBand,            "FreqBand",            REAL8,        'b', OPTIONAL,   "Search frequency band") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f1dot,               "f1dot",               REAL8,        0,   OPTIONAL,   "Spindown parameter") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_df1dot,              "df1dot",              REAL8,        0,   OPTIONAL,   "Spindown resolution (required if nonzero f1dotBand)") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f1dotBand,           "f1dotBand",           REAL8,        0,   OPTIONAL,   "Spindown Range") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f2dot,               "f2dot",               REAL8,        0,   OPTIONAL,   "2nd spindown parameter") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_df2dot,              "df2dot",              REAL8,        0,   OPTIONAL,   "2nd spindown resolution (required if nonzero f2dotBand)") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f2dotBand,           "f2dotBand",           REAL8,        0,   OPTIONAL,   "2nd spindown Range") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f3dot,               "f3dot",               REAL8,        0,   OPTIONAL,   "3rd spindown parameter") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_df3dot,              "df3dot",              REAL8,        0,   OPTIONAL,   "3rd spindown resolution (required if nonzero f3dotBand)") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f3dotBand,           "f3dotBand",           REAL8,        0,   OPTIONAL,   "3rd spindown Range") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_ThrF,                "peakThrF",            REAL8,        0,   OPTIONAL,   "Fstat Threshold") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_mismatch1,           "mismatch1",           REAL8,        'm', OPTIONAL,   "1st stage mismatch") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_gridType1,           "gridType1",           INT4,         0,   OPTIONAL,   "0=flat, 1=isotropic, 2=metric, 3=file") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_metricType1,         "metricType1",         INT4,         0,   OPTIONAL,   "0=none, 1=Ptole-analytic, 2=Ptole-numeric, 3=exact") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_gammaRefine,         "gammaRefine",         INT4,         'g', OPTIONAL,   "Refinement of fine grid (default: use segment times)") == XLAL_SUCCESS, XLAL_EFUNC);
  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);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_fnameout,            "fnameout",            STRING,       'o', OPTIONAL,   "Output filename") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_fnameChkPoint,       "fnameChkPoint",       STRING,       0,   OPTIONAL,   "Checkpoint filename") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nCand1,              "nCand1",              INT4,         'n', OPTIONAL,   "No. of candidates to output") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_printCand1,          "printCand1",          BOOLEAN,      0,   OPTIONAL,   "Print 1st stage candidates") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_refTime,             "refTime",             REAL8,        0,   OPTIONAL,   "Ref. time for pulsar pars [Default: mid-time]") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_ephemEarth,          "ephemEarth",          STRING,       0,   OPTIONAL,   "Location of Earth ephemeris file") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_ephemSun,            "ephemSun",            STRING,       0,   OPTIONAL,   "Location of Sun ephemeris file") == XLAL_SUCCESS, XLAL_EFUNC);
  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);
  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);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_printFstat1,         "printFstat1",         BOOLEAN,      0,   OPTIONAL,   "Print 1st stage Fstat vectors") == XLAL_SUCCESS, XLAL_EFUNC);
527 528 529

  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);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_SignalOnly,          "SignalOnly",          BOOLEAN,      'S', DEPRECATED, "DEPRECATED ALTERNATIVE: Don't estimate noise-floors but assume sqrtSX=1 instead") == XLAL_SUCCESS, XLAL_EFUNC);
530 531 532 533 534 535

  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nStacksMax,          "nStacksMax",          INT4,         0,   OPTIONAL,   "Maximum No. of segments" ) == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_tStack,              "tStack",              REAL8,        'T', OPTIONAL,   "Duration of segments (sec)" ) == XLAL_SUCCESS, XLAL_EFUNC);
  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);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_recalcToplistStats,  "recalcToplistStats",  BOOLEAN,      0,   OPTIONAL,   "Additional analysis for toplist candidates, recalculate 2F, 2FX at finegrid") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_loudestSegOutput,    "loudestSegOutput",    BOOLEAN,      0,   OPTIONAL,   "Output extra info about loudest segment; (requires --recalcToplistStats)") == XLAL_SUCCESS, XLAL_EFUNC);
536 537

  // ----- Line robust stats parameters ----------
538 539 540 541 542 543
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_computeBSGL,         "computeBSGL",         BOOLEAN,      0,   OPTIONAL,   "Compute and output line-robust statistic (BSGL)") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_Fstar0sc,            "Fstar0sc",            REAL8,        0,   OPTIONAL,   "BSGL: semi-coh transition-scale parameter 'Fstar0sc=Nseg*Fstar0coh'" ) == XLAL_SUCCESS, XLAL_EFUNC);
  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);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_BSGLlogcorr,         "BSGLlogcorr",         BOOLEAN,      0,   DEVELOPER,  "BSGL: include log-correction terms (slower) or not (faster)") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_getMaxFperSeg,       "getMaxFperSeg",       BOOLEAN,      0,   OPTIONAL,   "Compute and output maximum F and FX over segments") == XLAL_SUCCESS, XLAL_EFUNC);
  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'") == XLAL_SUCCESS, XLAL_EFUNC);
544
  // --------------------------------------------
545

546 547
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_FstatMethod,         "FstatMethod",         STRING,       0,   OPTIONAL,   "F-statistic method to use. Available methods: %s", XLALFstatMethodHelpString() ) == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_FstatMethodRecalc,   "FstatMethodRecalc",   STRING,       0,   OPTIONAL,   "F-statistic method to use for recalc. Available methods: %s", XLALFstatMethodHelpString() ) == XLAL_SUCCESS, XLAL_EFUNC);
548
  /* developer user variables */
549 550 551 552
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_blocksRngMed,        "blocksRngMed",        INT4,         0,   DEVELOPER,  "RngMed block size") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_SSBprecision,        "SSBprecision",        INT4,         0,   DEVELOPER,  "Precision for SSB transform.") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_Dterms,              "Dterms",              INT4,         0,   DEVELOPER,  "No. of terms to keep in Dirichlet Kernel" ) == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_skyPointIndex,       "skyPointIndex",       INT4,         0,   DEVELOPER,  "Only analyze this skypoint in grid" ) == XLAL_SUCCESS, XLAL_EFUNC);
553

554 555
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_outputTiming,        "outputTiming",        STRING,       0,   DEVELOPER,  "Append timing information into this file") == XLAL_SUCCESS, XLAL_EFUNC);
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_outputTimingDetails, "outputTimingDetails", STRING,       0,   DEVELOPER,  "Append detailed F-stat timing information to this file") == XLAL_SUCCESS, XLAL_EFUNC);
556

557
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_version,             "version",             BOOLEAN,      'V', SPECIAL,    "Output version information") == XLAL_SUCCESS, XLAL_EFUNC);
558

559 560
  /* inject signals into the data being analyzed */
  XLAL_CHECK_MAIN( XLALRegisterNamedUvar ( &uvar_injectionSources, "injectionSources",      STRINGVector, 0, DEVELOPER,     "CSV list of files containing signal parameters for injection [see mfdv5]") == XLAL_SUCCESS, XLAL_EFUNC );
561

562
  /* read all command line variables */
563
  BOOLEAN should_exit = 0;
564
  XLAL_CHECK_MAIN( XLALUserVarReadAllInput(&should_exit, argc, argv) == XLAL_SUCCESS, XLAL_EFUNC);
565 566
  if (should_exit)
    return(1);
567

568 569 570 571 572 573
  /* assemble version string */
  CHAR *VCSInfoString;
  if ( (VCSInfoString = XLALGetVersionString(0)) == NULL ) {
    XLALPrintError("XLALGetVersionString(0) failed.\n");
    return HIERARCHICALSEARCH_ESUB;
  }
574

575
  LogPrintfVerbatim( LOG_DEBUG, "Code-version: %s\n", VCSInfoString );
576
  // LogPrintfVerbatim( LOG_DEBUG, "CFS Hotloop variant: %s\n", OptimisedHotloopSource );
577 578 579 580 581

  if ( uvar_version )
    {
      printf ("%s\n", VCSInfoString );
      return (0);
582
    }
583

584 585
  /* some basic sanity checks on user vars */
  if ( uvar_nStacksMax < 1) {
586
    fprintf(stderr, "Invalid number of segments!\n");
587 588 589
    return( HIERARCHICALSEARCH_EBAD );
  }

590
#ifndef EXP_NO_NUM_COUNT
591 592
  /* check that the numbercount can't exceed the data type */
  {
593 594 595
    UINT8 maxseg = 1;
    maxseg = maxseg << (8*sizeof(FINEGRID_NC_T));
    maxseg -= 1;
Adam Mercer's avatar
Adam Mercer committed
596
    if ( (UINT8)uvar_nStacksMax > maxseg) {
597
      fprintf(stderr,
598 599 600
              "Number of segments exceeds %" LAL_UINT8_FORMAT "!\n"
              "Compile without GC_SSE2_OPT to extend the available segment range\n",
              maxseg);
601 602
      return( HIERARCHICALSEARCH_EBAD );
    }
603
  }
604
#endif
605

606 607 608 609 610 611
  if ( uvar_blocksRngMed < 1 ) {
    fprintf(stderr, "Invalid Running Median block size\n");
    return( HIERARCHICALSEARCH_EBAD );
  }

  if ( uvar_ThrF < 0 ) {
Holger Pletsch's avatar
Holger Pletsch committed
612
    fprintf(stderr, "Invalid value of F-statistic threshold\n");
613 614
    return( HIERARCHICALSEARCH_EBAD );
  }
615

616
  if ( uvar_f3dotBand != 0 && ( !XLALUserVarWasSet(&uvar_gammaRefine) || !XLALUserVarWasSet(&uvar_gammaRefine) || uvar_gammaRefine != 1 || uvar_gamma2Refine != 1 )){
617 618 619
	fprintf(stderr, "Search over 3rd spindown is available only with gammaRefine AND gamma2Refine manually set to 1!\n");
	return( HIERARCHICALSEARCH_EVAL );
  }
620

621 622 623
  /* check SignalOnly and assumeSqrtSX */
  XLAL_CHECK_MAIN ( !uvar_SignalOnly || (uvar_assumeSqrtSX == NULL), XLAL_EINVAL, "Cannot pass --SignalOnly AND --assumeSqrtSX at the same time!\n");

624
  /* 2F threshold for semicoherent stage */
625 626 627
#ifndef EXP_NO_NUM_COUNT
  REAL4 TwoFthreshold = 2.0 * uvar_ThrF;
#endif
628

629 630
  if ( (uvar_SortToplist < 0) || (uvar_SortToplist >= SORTBY_LAST) ) {
    XLALPrintError ( "Invalid value %d specified for toplist sorting, must be within [0, %d]\n", uvar_SortToplist, SORTBY_LAST - 1 );
631 632
    return( HIERARCHICALSEARCH_EBAD );
  }
633 634 635
  if ( (uvar_SortToplist == SORTBY_BSGL || uvar_SortToplist == SORTBY_DUAL_F_BSGL ||
        uvar_SortToplist == SORTBY_BSGLtL || uvar_SortToplist == SORTBY_BtSGLtL ||
        uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL ) && !uvar_computeBSGL ) {
636
    fprintf(stderr, "Toplist sorting by BSGL only possible if --computeBSGL given.\n");
637 638
    return( HIERARCHICALSEARCH_EBAD );
  }
639 640 641 642 643
  if ( ( uvar_SortToplist == SORTBY_BSGLtL || uvar_SortToplist == SORTBY_BtSGLtL ||
        uvar_SortToplist == SORTBY_TRIPLE_BStSGLtL ) && !uvar_getMaxFperSeg ) {
    fprintf(stderr, "Toplist sorting by B[t]SGLtL only possible if --getMaxFperSeg given.\n");
    return( HIERARCHICALSEARCH_EBAD );
  }
644

645
  /* create toplist -- semiCohToplist has the same structure
646
     as a fstat candidate, so treat it as a fstat candidate */
647
  if ( uvar_SortToplist == SORTBY_DUAL_F_BSGL )	// special treatement of 'dual' toplists: 1st one sorted by 'F', 2nd one by 'BSGL'
648
    {
649 650 651 652
      XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist, uvar_nCand1, SORTBY_F ),
                   XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_F );
      XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist2, uvar_nCand1, SORTBY_BSGL ),
                   XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BSGL );
653
    }
654 655 656 657 658 659 660 661 662
  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'
    {
      XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist, uvar_nCand1, SORTBY_BSGL ),
                   XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BSGL );
      XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist2, uvar_nCand1, SORTBY_BSGLtL ),
                   XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BSGLtL );
      XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist3, uvar_nCand1, SORTBY_BtSGLtL ),
                   XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, SORTBY_BtSGLtL );
    }
663
  else	// 'normal' single-sorting toplist cases (sortby 'F', 'nc' or 'BSGL')
664
    {
665 666
      XLAL_CHECK ( 0 == create_gctFstat_toplist ( &semiCohToplist, uvar_nCand1, uvar_SortToplist),
                   XLAL_EFUNC, "create_gctFstat_toplist() failed for nCand=%d and sortBy=%d\n", uvar_nCand1, uvar_SortToplist );
667
    }
668

669 670 671 672 673 674 675 676 677 678 679 680 681
#ifdef EAH_BOINC
  // BOINC Apps always checkpoint, so set a default filename here
  if (uvar_fnameChkPoint == NULL) {
    CHAR*fname = "checkpoint.cpt";
    uvar_fnameChkPoint = XLALMalloc(strlen(fname)+1);
    if (uvar_fnameChkPoint == NULL) {
      fprintf(stderr, "error allocating memory [HierarchSearchGCT.c %d]\n" , __LINE__);
      return(HIERARCHICALSEARCH_EMEM);
    }
    strcpy(uvar_fnameChkPoint, fname);
  }
#endif

682
  /* write the log file */
683
  if ( uvar_log )
684 685 686 687 688 689 690 691
    {
      fnamelog = LALCalloc( strlen(uvar_fnameout) + 1 + 4, sizeof(CHAR) );
      strcpy(fnamelog, uvar_fnameout);
      strcat(fnamelog, ".log");
      /* open the log file for writing */
      if ((fpLog = fopen(fnamelog, "wb")) == NULL) {
        fprintf(stderr, "Unable to open file %s for writing\n", fnamelog);
        LALFree(fnamelog);
692
        /*exit*/
693 694 695 696
        return(HIERARCHICALSEARCH_EFILE);
      }

      /* get the log string */
697
      XLAL_CHECK_MAIN( ( logstr = XLALUserVarGetLog(UVAR_LOGFMT_CFGFILE) ) != NULL, XLAL_EFUNC);
698

Holger Pletsch's avatar
Holger Pletsch committed
699
      fprintf( fpLog, "# Log file for HierarchSearchGCT.c\n\n");
700 701 702 703
      fprintf( fpLog, "# User Input:\n");
      fprintf( fpLog, "#-------------------------------------------\n");
      fprintf( fpLog, "# cmdline: %s\n", logstr );
      LALFree(logstr);
704

705
      /* add code version ID (only useful for git-derived versions) */
706
      fprintf ( fpLog, "# version: %s\n", VCSInfoString );
707

708 709
      fclose (fpLog);
      LALFree(fnamelog);
710

711
    } /* end of logging */
712

713
  tic_Start = GETTIME();
714
  BOOLEAN printHeader = 0;
715
  if ( uvar_outputTimingDetails != NULL ) {
716 717 718 719 720 721
    FILE *tmp;
    if ( (tmp = fopen ( uvar_outputTimingDetails, "r" )) == NULL ) {
      printHeader = 1;
    } else {
      fclose (tmp );
    }
722 723
    XLAL_CHECK ( (usefulParams.timingDetailsFP = fopen ( uvar_outputTimingDetails, "wb" )) != NULL, XLAL_ESYS, "Failed to open '%s' for writing\n", uvar_outputTimingDetails );
  } // if uvar_outputTimingDetails
724

725
  /* initializations of coarse and fine grids */
726
  coarsegrid.TwoF=NULL;
727
  coarsegrid.TwoFX=NULL;
728 729 730
  coarsegrid.Uindex=NULL;
  finegrid.nc= NULL;
  finegrid.sumTwoF=NULL;
731
  finegrid.sumTwoFX=NULL;
732 733
  finegrid.maxTwoFl=NULL;
  finegrid.maxTwoFXl=NULL;
734 735
  finegrid.maxTwoFlIdx=NULL;
  finegrid.maxTwoFXlIdx=NULL;
736 737

  /* initialize ephemeris info */
738 739
  EphemerisData *edat;
  XLAL_CHECK ( (edat = XLALInitBarycenter ( uvar_ephemEarth, uvar_ephemSun )) != NULL, XLAL_EFUNC );
740 741

  XLALGPSSetREAL8(&minStartTimeGPS, uvar_minStartTime1);
742
  XLALGPSSetREAL8(&maxStartTimeGPS, uvar_maxStartTime1);
743

744 745 746 747 748 749 750 751 752 753
  /* create output files for writing if requested by user */
  if ( uvar_printCand1 )
    {
      fnameSemiCohCand = LALCalloc( strlen(uvar_fnameout) + 1, sizeof(CHAR) );
      if ( fnameSemiCohCand == NULL) {
        fprintf(stderr, "error allocating memory [HierarchSearchGCT.c %d]\n" , __LINE__);
        return(HIERARCHICALSEARCH_EMEM);
      }
      strcpy(fnameSemiCohCand, uvar_fnameout);
    }
754

755 756 757 758 759 760
  if ( uvar_printFstat1 )
    {
      const CHAR *append = "_fstatVec1.dat";
      fnameFstatVec1 = LALCalloc( strlen(uvar_fnameout) + strlen(append) + 1, sizeof(CHAR) );
      strcpy(fnameFstatVec1, uvar_fnameout);
      strcat(fnameFstatVec1, append);
761
      if ( !(fpFstat1 = fopen( fnameFstatVec1, "wb")))  {
762 763 764 765 766 767 768
        fprintf ( stderr, "Unable to open Fstat file fstatvec1.out for writing.\n");
        return (HIERARCHICALSEARCH_EFILE);
      }
    }

  /*------------ Set up stacks, detector states etc. */
  /* initialize spin range vectors */
769
  XLAL_INIT_MEM( spinRange_Temp );
770 771 772

  /* some useful first stage params */
  usefulParams.sftbasename = uvar_DataFiles1;
773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798

  /* ----- prepare generation of coherent segment list */
  if ( XLALUserVarWasSet ( &uvar_segmentList ) )
    {
      if ( XLALUserVarWasSet ( &uvar_nStacksMax ) || XLALUserVarWasSet ( &uvar_tStack ) ) {
        XLALPrintError ( "Use EITHER (--nStacksMax and --tStack) OR --segmentList to define the coherent segments!\n\n" );
        return HIERARCHICALSEARCH_EBAD;
      }
      if ( (usefulParams.segmentList = XLALReadSegmentsFromFile ( uvar_segmentList )) == NULL ) {
        XLALPrintError ("Failed to parse segment-list file '%s'. xlalErrno = %d\n", uvar_segmentList, xlalErrno );
        return HIERARCHICALSEARCH_ESUB;
      }
    }
  else /* set up maximally nStacksMax fixed-size segments of length tStack */
    {
      if ( !XLALUserVarWasSet ( &uvar_tStack ) ) {
        XLALPrintError ( "Need to set --tStack or --segmentList to define the coherent segments!\n\n" );
        return HIERARCHICALSEARCH_EBAD;
      }

      usefulParams.nStacks = uvar_nStacksMax;
      usefulParams.tStack = uvar_tStack;
      usefulParams.segmentList = NULL;
    }
  /* ----- */

799

800 801 802 803
  XLAL_INIT_MEM ( usefulParams.spinRange_startTime );
  XLAL_INIT_MEM ( usefulParams.spinRange_endTime );
  XLAL_INIT_MEM ( usefulParams.spinRange_refTime );
  XLAL_INIT_MEM ( usefulParams.spinRange_midTime );
804 805 806 807 808

  /* copy user specified spin variables at reftime  */
  /* the reference time value in spinRange_refTime will be set in SetUpSFTs() */
  usefulParams.spinRange_refTime.fkdot[0] = uvar_Freq; /* frequency */
  usefulParams.spinRange_refTime.fkdot[1] = uvar_f1dot;  /* 1st spindown */
809
  usefulParams.spinRange_refTime.fkdot[2] = uvar_f2dot;  /* 2nd spindown */
810
  usefulParams.spinRange_refTime.fkdot[3] = uvar_f3dot;  /* 3rd spindown */  
811 812
  usefulParams.spinRange_refTime.fkdotBand[0] = uvar_FreqBand; /* frequency range */
  usefulParams.spinRange_refTime.fkdotBand[1] = uvar_f1dotBand; /* spindown range */
813
  usefulParams.spinRange_refTime.fkdotBand[2] = uvar_f2dotBand; /* spindown range */
814
  usefulParams.spinRange_refTime.fkdotBand[3] = uvar_f3dotBand; /* spindown range */
815 816 817

  usefulParams.edat = edat;
  usefulParams.minStartTimeGPS = minStartTimeGPS;
818
  usefulParams.maxStartTimeGPS = maxStartTimeGPS;
819 820
  usefulParams.blocksRngMed = uvar_blocksRngMed;
  usefulParams.Dterms = uvar_Dterms;
821
  usefulParams.assumeSqrtSX = uvar_assumeSqrtSX;
Holger Pletsch's avatar
Holger Pletsch committed
822
  usefulParams.SignalOnly = uvar_SignalOnly;
823
  usefulParams.SSBprec = uvar_SSBprecision;