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
synthesizeTransientStats.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010 Reinhard Prix
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with with program; see the file COPYING. If not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17 * MA 02110-1301 USA
18 */
19
20/*********************************************************************************/
21/**
22 * \author R. Prix
23 * \file
24 * \ingroup lalpulsar_bin_Fstatistic
25 * \brief
26 * Generate N samples of various statistics (F-stat, B-stat,...) drawn from their
27 * respective distributions, assuming Gaussian noise, and drawing signal params from
28 * their (given) priors
29 *
30 * This is based on synthesizeBstat, and is mostly meant to be used for Monte-Carlos
31 * studies of ROC curves
32 *
33 */
34
35/*
36 *
37 * Some possible use-cases to consider
38 * - transient search BF-stat (synthesize atoms)
39 * - line-veto studies (generate line-realizations)
40 * - different B-stats from different prior models (to avoid integration)
41 *
42 */
43
44#include "config.h"
45
46#include <stdio.h>
47#include <stdbool.h>
48
49#include <gsl/gsl_rng.h>
50
51#include <lal/LALString.h>
52#include <lal/SkyCoordinates.h>
53#include <lal/AVFactories.h>
54#include <lal/LALInitBarycenter.h>
55#include <lal/UserInput.h>
56#include <lal/LogPrintf.h>
57#include <lal/ComputeFstat.h>
58#include <lal/TransientCW_utils.h>
59#include <lal/ProbabilityDensity.h>
60#include <lal/SynthesizeCWDraws.h>
61#include <lal/LALPulsarVCSInfo.h>
62
63/*---------- DEFINES ----------*/
64#define SQ(x) ((x)*(x))
65#define CUBE(x) ((x)*(x)*(x))
66#define QUAD(x) ((x)*(x)*(x)*(x))
67
68/* define min/max macros if not already defined */
69#ifndef min
70#define min(a,b) ((a)<(b)?(a):(b))
71#endif
72#ifndef max
73#define max(a,b) ((a)>(b)?(a):(b))
74#endif
75
76/* ---------- local types ---------- */
77
78/** User-variables: can be set from config-file or command-line */
79typedef struct {
80 /* amplitude parameters + ranges: 4 alternative ways to specify the h0-prior (set to < 0 to deactivate all but one!) */
81 REAL8 fixedh0Nat; /**< Alternative 1: if >=0 ==> fix the GW amplitude: h0/sqrt(Sn) */
82 REAL8 fixedSNR; /**< Alternative 2: if >=0 ==> fix the optimal SNR of the injected signals */
83 REAL8 fixedh0NatMax; /**< Alternative 3: if >=0 ==> draw GW amplitude h0 in [0, h0NatMax ]: <==> 'regularized' F-stat prior (obsolete) */
84 REAL8 fixedRhohMax; /**< Alternative 4: if >=0 ==> draw rhoh=h0*(detM)^(1/8) in [0, rhohMax]: <==> canonical F-stat prior */
85
86 REAL8 cosi; /**< cos(inclination angle). If not set: randomize within [-1,1] */
87 REAL8 psi; /**< polarization angle psi. If not set: randomize within [-pi/4,pi/4] */
88 REAL8 phi0; /**< initial GW phase phi_0. If not set: randomize within [0, 2pi] */
89 INT4 AmpPriorType; /**< enumeration of types of amplitude-priors: 0=physical, 1=canonical */
90
91 /* Doppler parameters */
92 REAL8 Alpha; /**< skyposition Alpha (RA) in radians */
93 REAL8 Delta; /**< skyposition Delta (Dec) in radians */
94
95 /* transient window ranges: for injection ... */
96 CHAR *injectWindow_type; /**< Type of transient window to inject ('none', 'rect', 'exp'). */
97 REAL8 injectWindow_t0Days; /**< Earliest start-time of transient window to inject, as offset in days from dataStartGPS */
98 REAL8 injectWindow_t0DaysBand;/**< Range of start-times of transient window to inject, in days */
99 REAL8 injectWindow_tauDays; /**< Shortest transient-window timescale to inject, in days */
100 REAL8 injectWindow_tauDaysBand; /**< Range of transient-window timescale to inject, in days */
101 /* ... and for search */
102 CHAR *searchWindow_type; /**< Type of transient window to search with ('none', 'rect', 'exp'). */
103 REAL8 searchWindow_t0Days; /**< Earliest start-time of transient window to search, as offset in days from dataStartGPS */
104 REAL8 searchWindow_t0DaysBand;/**< Range of start-times of transient window to inject, in days */
105 REAL8 searchWindow_tauDays; /**< Shortest transient-window timescale to search, in days */
106 REAL8 searchWindow_tauDaysBand; /**< Range of transient-window timescale to search, in days */
107 INT4 searchWindow_dt0; /**< Step-size for search/marginalization over transient-window start-time, in seconds [Default:Tsft] */
108 INT4 searchWindow_dtau; /**< Step-size for search/marginalization over transient-window timescale, in seconds [Default:Tsft] */
109
110 /* other parameters */
111 CHAR *IFO; /**< IFO name [deprecated] */
112 LALStringVector *IFOs; /**< list of detector-names "H1,H2,L1,.." or single detector*/
113 INT4 dataStartGPS; /**< data start-time in GPS seconds */
114 INT4 dataDuration; /**< data-span to generate */
115 LALStringVector *timestampsFiles; /**< per-detector file(s) with timestamps list */
116 INT4 TAtom; /**< Fstat atoms time baseline */
117
118 BOOLEAN computeFtotal; /**< Also compute 'total' F-statistic over the full data-span */
119 INT4 numDraws; /**< number of random 'draws' to simulate for F-stat and B-stat */
120
121 CHAR *outputStats; /**< output file to write numDraw resulting statistics into */
122 CHAR *outputAtoms; /**< output F-statistic atoms into a file with this basename */
123 CHAR *outputFstatMap; /**< output F-statistic over 2D parameter space {t0, tau} into file with this basename */
124 CHAR *outputInjParams;/**< output injection parameters into this file */
125 CHAR *outputPosteriors;/**< output posterior pdfs on t0 and tau */
126 BOOLEAN SignalOnly; /**< dont generate noise-draws: will result in non-random 'signal only' values of F and B */
127
128 BOOLEAN useFReg; /**< use 'regularized' Fstat (1/D)*e^F for marginalization, or 'standard' e^F */
129
130 CHAR *ephemEarth; /**< Earth ephemeris file to use */
131 CHAR *ephemSun; /**< Sun ephemeris file to use */
132
133 INT4 randSeed; /**< GSL random-number generator seed value to use */
135
136/**
137 * Configuration settings required for and defining a coherent pulsar search.
138 * These are 'pre-processed' settings, which have been derived from the user-input.
139 */
140typedef struct {
141 AmplitudePrior_t AmpPrior; /**< amplitude-parameter priors to draw signals from */
142
143 SkyPosition skypos; /**< (Alpha,Delta,system). Use Alpha < 0 to signal 'allsky' */
144 BOOLEAN SignalOnly; /**< dont generate noise-draws: will result in non-random 'signal only' values of F and B */
145
146 transientWindowRange_t transientSearchRange; /**< transient-window range for the search (flat priors) */
147 transientWindowRange_t transientInjectRange; /**< transient-window range for injections (flat priors) */
148
149 MultiDetectorStateSeries *multiDetStates; /**< multi-detector state series covering observation time */
150
151 gsl_rng *rng; /**< gsl random-number generator */
152 CHAR *logString; /**< logstring for file-output, containing cmdline-options + code VCS version info */
153
155
156/* ---------- local prototypes ---------- */
157int main( int argc, char *argv[] );
158
159int XLALInitUserVars( UserInput_t *uvar );
160int XLALInitCode( ConfigVariables *cfg, const UserInput_t *uvar );
161int XLALInitAmplitudePrior( AmplitudePrior_t *AmpPrior, const UserInput_t *uvar );
162
163/* exportable API */
164
165/*---------- Global variables ----------*/
166extern int vrbflg; /**< defined in lal/lib/std/LALError.c */
167
168/*----------------------------------------------------------------------*/
169/* Main Function starts here */
170/*----------------------------------------------------------------------*/
171/**
172 * MAIN function
173 * Generates samples of B-stat and F-stat according to their pdfs for given signal-params.
174 */
175int main( int argc, char *argv[] )
176{
178 ConfigVariables XLAL_INIT_DECL( cfg ); /**< various derived configuration settings */
179
180 vrbflg = 1; /* verbose error-messages */
181
182 /* turn off default GSL error handler */
183 gsl_set_error_handler_off();
184
185 /* ----- register and read all user-variables ----- */
186 if ( XLALInitUserVars( &uvar ) != XLAL_SUCCESS ) {
187 LogPrintf( LOG_CRITICAL, "%s: XLALInitUserVars() failed with errno=%d\n", __func__, xlalErrno );
188 return 1;
189 }
190
191 /* do ALL cmdline and cfgfile handling */
192 BOOLEAN should_exit = 0;
193 if ( XLALUserVarReadAllInput( &should_exit, argc, argv, lalPulsarVCSInfoList ) != XLAL_SUCCESS ) {
194 LogPrintf( LOG_CRITICAL, "%s: XLALUserVarReadAllInput() failed with errno=%d\n", __func__, xlalErrno );
195 return 1;
196 }
197 if ( should_exit ) {
198 return EXIT_FAILURE;
199 }
200
201 /* ---------- Initialize code-setup ---------- */
202 if ( XLALInitCode( &cfg, &uvar ) != XLAL_SUCCESS ) {
203 LogPrintf( LOG_CRITICAL, "%s: XLALInitCode() failed with error = %d\n", __func__, xlalErrno );
205 }
206
207 /* ----- prepare stats output ----- */
208 LALFILE *fpTransientStats = NULL;
209 if ( uvar.outputStats ) {
210 if ( ( fpTransientStats = XLALFileOpen( uvar.outputStats, "wb" ) ) == NULL ) {
211 LogPrintf( LOG_CRITICAL, "Error opening file '%s' for writing..\n\n", uvar.outputStats );
213 }
214 XLALFilePrintf( fpTransientStats, "%s", cfg.logString ); /* write search log comment */
215 if ( write_transientCandidate_to_fp( fpTransientStats, NULL, 'd' ) != XLAL_SUCCESS ) { /* write header-line comment */
217 }
218 } /* if outputStats */
219
220 /* ----- prepare injection params output ----- */
221 FILE *fpInjParams = NULL;
222 if ( uvar.outputInjParams ) {
223 if ( ( fpInjParams = fopen( uvar.outputInjParams, "wb" ) ) == NULL ) {
224 LogPrintf( LOG_CRITICAL, "Error opening file '%s' for writing..\n\n", uvar.outputInjParams );
226 }
227 fprintf( fpInjParams, "%s", cfg.logString ); /* write search log comment */
228 if ( write_InjParams_to_fp( fpInjParams, NULL, 0, 0, 0 ) != XLAL_SUCCESS ) { /* write header-line comment - options outputMmunuX and numDetectors not supported here, so pass defaults to deactivate them */
230 }
231 } /* if outputInjParams */
232
233 /* ----- main MC loop over numDraws trials ---------- */
234 multiAMBuffer_t XLAL_INIT_DECL( multiAMBuffer ); /* prepare AM-buffer */
235 INT4 i;
236
237 for ( i = 0; i < uvar.numDraws; i ++ ) {
238 InjParams_t XLAL_INIT_DECL( injParamsDrawn );
239
240 /* ----- generate signal random draws from ranges and generate Fstat atoms */
241 MultiFstatAtomVector *multiAtoms;
242 multiAtoms = XLALSynthesizeTransientAtoms( &injParamsDrawn, cfg.skypos, cfg.AmpPrior, cfg.transientInjectRange, cfg.multiDetStates, cfg.SignalOnly, &multiAMBuffer, cfg.rng, -1, NULL ); // options lineX and noise_weights not supported here, so pass defaults to deactivate them
243 if ( multiAtoms == NULL ) {
244 LogPrintf( LOG_CRITICAL, "%s: XLALSynthesizeTransientAtoms() failed with xlalErrno = %d\n", __func__, xlalErrno );
246 }
247
248 /* ----- if requested, output signal injection parameters into file */
249 if ( fpInjParams && ( write_InjParams_to_fp( fpInjParams, &injParamsDrawn, uvar.dataStartGPS, 0, 0 ) ) != XLAL_SUCCESS ) { // options outputMmunuX and numDetectors not supported here, so pass defaults to deactivate them
251 } /* if fpInjParams & failure*/
252
253
254 /* ----- add meta-info on current transient-CW candidate */
256 cand.doppler.Alpha = multiAMBuffer.skypos.longitude;
257 cand.doppler.Delta = multiAMBuffer.skypos.latitude;
258 cand.windowRange = cfg.transientSearchRange;
259
260 /* ----- if needed: compute transient-Bstat search statistic on these atoms */
261 if ( fpTransientStats || uvar.outputFstatMap || uvar.outputPosteriors ) {
262 /* compute Fstat map F_mn over {t0, tau} */
263 if ( ( cand.FstatMap = XLALComputeTransientFstatMap( multiAtoms, cand.windowRange, uvar.useFReg ) ) == NULL ) {
264 XLALPrintError( "%s: XLALComputeTransientFstatMap() failed with xlalErrno = %d.\n", __func__, xlalErrno );
266 }
267 } /* if we'll need the Fstat-map F_mn */
268
269 /* ----- if requested compute marginalized Bayes factor */
270 if ( fpTransientStats ) {
271 cand.logBstat = XLALComputeTransientBstat( cand.windowRange, cand.FstatMap );
273 if ( err ) {
274 XLALPrintError( "%s: XLALComputeTransientBstat() failed with xlalErrno = %d\n", __func__, err );
276 }
277
278 if ( uvar.SignalOnly ) {
279 cand.FstatMap->maxF += 2;
280 cand.logBstat += 2;
281 }
282
283 } /* if Bstat requested */
284
285 /* ----- if requested, compute parameter posteriors for {t0, tau} */
286 pdf1D_t *pdf_t0 = NULL;
287 pdf1D_t *pdf_tau = NULL;
288 if ( fpTransientStats || uvar.outputPosteriors ) {
289 if ( ( pdf_t0 = XLALComputeTransientPosterior_t0( cand.windowRange, cand.FstatMap ) ) == NULL ) {
290 XLALPrintError( "%s: failed to compute t0-posterior\n", __func__ );
292 }
293 if ( ( pdf_tau = XLALComputeTransientPosterior_tau( cand.windowRange, cand.FstatMap ) ) == NULL ) {
294 XLALPrintError( "%s: failed to compute tau-posterior\n", __func__ );
296 }
297 /* get maximum-posterior estimate (MP) from the modes of these pdfs */
298 cand.t0_MP = XLALFindModeOfPDF1D( pdf_t0 );
299 if ( xlalErrno ) {
300 XLALPrintError( "%s: mode-estimation failed for pdf_t0. xlalErrno = %d\n", __func__, xlalErrno );
302 }
303 cand.tau_MP = XLALFindModeOfPDF1D( pdf_tau );
304 if ( xlalErrno ) {
305 XLALPrintError( "%s: mode-estimation failed for pdf_tau. xlalErrno = %d\n", __func__, xlalErrno );
307 }
308
309 } // if posteriors required
310
311 /* ----- if requested, compute Ftotal over full data-span */
312 if ( uvar.computeFtotal ) {
313 transientFstatMap_t *FtotalMap;
314 /* prepare special window to cover all the data with one F-stat calculation == Ftotal */
316 winRangeAll.type = TRANSIENT_NONE;
317
318 BOOLEAN useFReg = false;
319 if ( ( FtotalMap = XLALComputeTransientFstatMap( multiAtoms, winRangeAll, useFReg ) ) == NULL ) {
320 XLALPrintError( "%s: XLALComputeTransientFstatMap() failed with xlalErrno = %d.\n", __func__, xlalErrno );
322 }
323
324 /* we only use twoFtotal = 2 * maxF from this single-Fstat calculation */
325 REAL8 twoFtotal = 2.0 * FtotalMap->maxF;
326 if ( uvar.SignalOnly ) {
327 twoFtotal += 4;
328 }
329
330 /* ugly hack: lacking a good container for twoFtotal, we borrow fkdot[3] for this here ;) [only used for paper-MCs] */
331 cand.doppler.fkdot[3] = twoFtotal;
332
333 /* good riddance .. */
334 XLALDestroyTransientFstatMap( FtotalMap );
335
336 } /* if computeFtotal */
337
338 /* ----- if requested, output atoms-vector into file */
339 if ( uvar.outputAtoms ) {
340
341 LALFILE *fpAtoms;
342 char *fnameAtoms;
343 UINT4 len = strlen( uvar.outputAtoms ) + 20;
344 if ( ( fnameAtoms = XLALCalloc( 1, len ) ) == NULL ) {
345 XLALPrintError( "%s: failed to XLALCalloc ( 1, %d )\n", __func__, len );
347 }
348 sprintf( fnameAtoms, "%s_%04d_of_%04d.dat", uvar.outputAtoms, i + 1, uvar.numDraws );
349
350 if ( ( fpAtoms = XLALFileOpen( fnameAtoms, "wb" ) ) == NULL ) {
351 XLALPrintError( "%s: failed to open atoms-output file '%s' for writing.\n", __func__, fnameAtoms );
353 }
354 XLALFilePrintf( fpAtoms, "%s", cfg.logString ); /* output header info */
355
356 if ( write_MultiFstatAtoms_to_fp( fpAtoms, multiAtoms ) != XLAL_SUCCESS ) {
357 XLALPrintError( "%s: failed to write atoms to output file '%s'. xlalErrno = %d\n", __func__, fnameAtoms, xlalErrno );
359 }
360
361 XLALFree( fnameAtoms );
362 XLALFileClose( fpAtoms );
363 } /* if outputAtoms */
364
365 /* ----- if requested, output Fstat-map over {t0, tau} */
366 if ( uvar.outputFstatMap ) {
367 FILE *fpFstatMap;
368 char *fnameFstatMap;
369 UINT4 len = strlen( uvar.outputFstatMap ) + 20;
370 if ( ( fnameFstatMap = XLALCalloc( 1, len ) ) == NULL ) {
371 XLALPrintError( "%s: failed to XLALCalloc ( 1, %d )\n", __func__, len );
373 }
374 sprintf( fnameFstatMap, "%s_%04d_of_%04d.dat", uvar.outputFstatMap, i + 1, uvar.numDraws );
375
376 if ( ( fpFstatMap = fopen( fnameFstatMap, "wb" ) ) == NULL ) {
377 XLALPrintError( "%s: failed to open Fstat-map output file '%s' for writing.\n", __func__, fnameFstatMap );
379 }
380 fprintf( fpFstatMap, "%s", cfg.logString ); /* output header info */
381
382 fprintf( fpFstatMap, "\nFstat_mn = \\\n" );
383 if ( XLALfprintfGSLmatrix( fpFstatMap, "%.9g", cand.FstatMap->F_mn ) != XLAL_SUCCESS ) {
384 XLALPrintError( "%s: XLALfprintfGSLmatrix() failed.\n", __func__ );
386 }
387
388 XLALFree( fnameFstatMap );
389 fclose( fpFstatMap );
390
391 } /* if outputFstatMap */
392
393 /* ----- if requested, output posterior pdfs on transient params {t0, tau} into a file */
394 if ( uvar.outputPosteriors ) {
395 FILE *fpPosteriors;
396 char *fnamePosteriors;
397 UINT4 len = strlen( uvar.outputPosteriors ) + 20;
398 if ( ( fnamePosteriors = XLALCalloc( 1, len ) ) == NULL ) {
399 XLALPrintError( "%s: failed to XLALCalloc ( 1, %d )\n", __func__, len );
401 }
402 sprintf( fnamePosteriors, "%s_%04d_of_%04d.dat", uvar.outputPosteriors, i + 1, uvar.numDraws );
403
404 if ( ( fpPosteriors = fopen( fnamePosteriors, "wb" ) ) == NULL ) {
405 XLALPrintError( "%s: failed to open posteriors-output file '%s' for writing.\n", __func__, fnamePosteriors );
407 }
408 fprintf( fpPosteriors, "%s", cfg.logString ); /* output header info */
409
410 /* write them to file, using pdf-method */
411 if ( XLALOutputPDF1D_to_fp( fpPosteriors, pdf_t0, "pdf_t0" ) != XLAL_SUCCESS ) {
412 XLALPrintError( "%s: failed to output t0-posterior to file '%s'.\n", __func__, fnamePosteriors );
414 }
415 if ( XLALOutputPDF1D_to_fp( fpPosteriors, pdf_tau, "pdf_tau" ) != XLAL_SUCCESS ) {
416 XLALPrintError( "%s: failed to output tau-posterior to file '%s'.\n", __func__, fnamePosteriors );
418 }
419
420 /* free mem, close file */
421 XLALFree( fnamePosteriors );
422 fclose( fpPosteriors );
423
424 } /* if outputPosteriors */
425
426
427 /* ----- if requested, output transient-cand statistics */
428 if ( fpTransientStats && write_transientCandidate_to_fp( fpTransientStats, &cand, 'd' ) != XLAL_SUCCESS ) {
429 XLALPrintError( "%s: write_transientCandidate_to_fp() failed.\n", __func__ );
431 }
432
433 /* ----- free Memory */
434 XLALDestroyTransientFstatMap( cand.FstatMap );
436 XLALDestroyPDF1D( pdf_t0 );
437 XLALDestroyPDF1D( pdf_tau );
438
439 } /* for i < numDraws */
440
441 /* ----- close files ----- */
442 XLALFileClose( fpTransientStats );
443 if ( fpInjParams ) {
444 fclose( fpInjParams );
445 }
446
447 /* ----- free memory ---------- */
448 XLALDestroyMultiDetectorStateSeries( cfg.multiDetStates );
449 XLALDestroyMultiAMCoeffs( multiAMBuffer.multiAM );
451 /* ----- free amplitude prior pdfs ----- */
452 XLALDestroyPDF1D( cfg.AmpPrior.pdf_h0Nat );
453 XLALDestroyPDF1D( cfg.AmpPrior.pdf_cosi );
454 XLALDestroyPDF1D( cfg.AmpPrior.pdf_psi );
455 XLALDestroyPDF1D( cfg.AmpPrior.pdf_phi0 );
456
457 if ( cfg.logString ) {
458 XLALFree( cfg.logString );
459 }
460 gsl_rng_free( cfg.rng );
461
463
464 /* did we forget anything ? (doesn't cover gsl-memory!) */
466
467 return 0;
468
469} /* main() */
470
471/**
472 * Register all our "user-variables" that can be specified from cmd-line and/or config-file.
473 * Here we set defaults for some user-variables and register them with the UserInput module.
474 */
475int
477{
478 /* set a few defaults */
479 uvar->outputStats = NULL;
480
481 uvar->Alpha = -1; /* Alpha < 0 indicates "allsky" */
482 uvar->Delta = 0;
483
484 uvar->phi0 = 0;
485 uvar->psi = 0;
486
487 uvar->dataStartGPS = 814838413; /* 1 Nov 2005, ~ start of S5 */
488 uvar->dataDuration = ( INT4 ) round( LAL_YRSID_SI ); /* 1 year of data */
489
490 uvar->ephemEarth = XLALStringDuplicate( "earth00-40-DE405.dat.gz" );
491 uvar->ephemSun = XLALStringDuplicate( "sun00-40-DE405.dat.gz" );
492
493 uvar->numDraws = 1;
494 uvar->TAtom = 1800;
495
496 uvar->computeFtotal = 0;
497 uvar->useFReg = 0;
498
499 uvar->fixedh0Nat = -1;
500 uvar->fixedSNR = -1;
501 uvar->fixedh0NatMax = -1;
502 uvar->fixedRhohMax = -1;
503
504#define DEFAULT_IFO "H1"
505 uvar->IFO = XLALMalloc( strlen( DEFAULT_IFO ) + 1 );
506 strcpy( uvar->IFO, DEFAULT_IFO );
507 if ( ( uvar->IFOs = XLALCreateStringVector( "H1", NULL ) ) == NULL ) {
508 LogPrintf( LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
510 }
511
512 /* ---------- transient window defaults ---------- */
513#define DEFAULT_TRANSIENT "rect"
514 uvar->injectWindow_type = XLALMalloc( strlen( DEFAULT_TRANSIENT ) + 1 );
515 strcpy( uvar->injectWindow_type, DEFAULT_TRANSIENT );
516 uvar->searchWindow_type = XLALMalloc( strlen( DEFAULT_TRANSIENT ) + 1 );
517 strcpy( uvar->searchWindow_type, DEFAULT_TRANSIENT );
518
519 uvar->injectWindow_tauDays = 1.0;
520 uvar->injectWindow_tauDaysBand = 13.0;
521
522 REAL8 tauMaxDays = ( uvar->injectWindow_tauDays + uvar->injectWindow_tauDaysBand );
523 /* default window-ranges are t0 in [dataStartTime, dataStartTime - 3 * tauMax] */
524 uvar->injectWindow_t0Days = 0; // offset in days from uvar->dataStartGPS
525 uvar->injectWindow_t0DaysBand = fmax( 0.0, 1.0 * uvar->dataDuration / DAY24 - TRANSIENT_EXP_EFOLDING * tauMaxDays ); /* make sure it's >= 0 */
526
527 /* search-windows by default identical to inject-windows */
528 uvar->searchWindow_t0Days = uvar->injectWindow_t0Days;
529 uvar->searchWindow_t0DaysBand = uvar->injectWindow_t0DaysBand;
530 uvar->searchWindow_tauDays = uvar->injectWindow_tauDays;
531 uvar->searchWindow_tauDaysBand = uvar->injectWindow_tauDaysBand;
532
533 uvar->searchWindow_dt0 = uvar->TAtom;
534 uvar->searchWindow_dtau = uvar->TAtom;
535
536 /* register all our user-variables */
537 /* signal Doppler parameters */
538 XLALRegisterUvarMember( Alpha, REAL8, 'a', OPTIONAL, "Sky position alpha (equatorial coordinates) in radians [Default:allsky]" );
539 XLALRegisterUvarMember( Delta, REAL8, 'd', OPTIONAL, "Sky position delta (equatorial coordinates) in radians [Default:allsky]" );
540
541 /* signal amplitude parameters */
542 XLALRegisterUvarMember( fixedh0Nat, REAL8, 0, OPTIONAL, "Alternative 1: if >=0 fix the GW amplitude: h0/sqrt(Sn)" );
543 XLALRegisterUvarMember( fixedSNR, REAL8, 0, OPTIONAL, "Alternative 2: if >=0 fix the optimal SNR of the injected signals" );
544 XLALRegisterUvarMember( fixedh0NatMax, REAL8, 0, OPTIONAL, "Alternative 3: if >=0 draw GW amplitude h0 in [0, h0NatMax ] (FReg prior)" );
545 XLALRegisterUvarMember( fixedRhohMax, REAL8, 0, OPTIONAL, "Alternative 4: if >=0 draw rhoh=h0*(detM)^(1/8) in [0, rhohMax] (canonical F-stat prior)" );
546
547 XLALRegisterUvarMember( cosi, REAL8, 'i', OPTIONAL, "cos(inclination angle). If not set: randomize within [-1,1]." );
548 XLALRegisterUvarMember( psi, REAL8, 0, OPTIONAL, "polarization angle psi. If not set: randomize within [-pi/4,pi/4]." );
549 XLALRegisterUvarMember( phi0, REAL8, 0, OPTIONAL, "initial GW phase phi_0. If not set: randomize within [0, 2pi]" );
550
551 XLALRegisterUvarMember( AmpPriorType, INT4, 0, OPTIONAL, "Enumeration of types of amplitude-priors: 0=physical, 1=canonical" );
552
553 XLALRegisterUvarMember( IFO, STRING, 0, DEPRECATED, "Detector: 'G1','L1','H1,'H2', 'V1', ... (use --IFOs instead)" );
554 XLALRegisterUvarMember( IFOs, STRINGVector, 'I', OPTIONAL, "Comma-separated list of detectors, e.g. \"H1,H2,L1,G1, ...\" " );
555 XLALRegisterUvarMember( dataStartGPS, INT4, 0, OPTIONAL, "data start-time in GPS seconds" );
556 XLALRegisterUvarMember( dataDuration, INT4, 0, OPTIONAL, "data-span to generate (in seconds)" );
557 XLALRegisterUvarMember( timestampsFiles, STRINGVector, 0, OPTIONAL, "ALTERNATIVE: file(s) to read timestamps from (file-format: lines with <seconds> [<nanoseconds>]; nanoseconds currently ignored; only 1 detector/file currently supported)" );
558
559 /* transient window ranges: for injection ... */
560 XLALRegisterUvarMember( injectWindow_type, STRING, 0, OPTIONAL, "Type of transient window to inject ('none', 'rect', 'exp')" );
561 XLALRegisterUvarMember( injectWindow_tauDays, REAL8, 0, OPTIONAL, "Shortest transient-window timescale to inject, in days" );
562 XLALRegisterUvarMember( injectWindow_tauDaysBand, REAL8, 0, OPTIONAL, "Range of transient-window timescale to inject, in days" );
563 XLALRegisterUvarMember( injectWindow_t0Days, REAL8, 0, OPTIONAL, "Earliest start-time of transient window to inject, as offset in days from dataStartGPS" );
564 XLALRegisterUvarMember( injectWindow_t0DaysBand, REAL8, 0, OPTIONAL, "Range of GPS start-time of transient window to inject, in days [Default:dataDuration-3*tauMax]" );
565 /* ... and for search */
566 XLALRegisterUvarMember( searchWindow_type, STRING, 0, OPTIONAL, "Type of transient window to search with ('none', 'rect', 'exp') [Default:injectWindow]" );
567 XLALRegisterUvarMember( searchWindow_tauDays, REAL8, 0, OPTIONAL, "Shortest transient-window timescale to search, in days [Default:injectWindow]" );
568 XLALRegisterUvarMember( searchWindow_tauDaysBand, REAL8, 0, OPTIONAL, "Range of transient-window timescale to search, in days [Default:injectWindow]" );
569 XLALRegisterUvarMember( searchWindow_t0Days, REAL8, 0, OPTIONAL, "Earliest start-time of transient window to search, as offset in days from dataStartGPS [Default:injectWindow]" );
570 XLALRegisterUvarMember( searchWindow_t0DaysBand, REAL8, 0, OPTIONAL, "Range of GPS start-time of transient window to search, in days [Default:injectWindow]" );
571
572 XLALRegisterUvarMember( searchWindow_dtau, INT4, 0, OPTIONAL, "Step-size for search/marginalization over transient-window timescale, in seconds [Default:TAtom]" );
573 XLALRegisterUvarMember( searchWindow_dt0, INT4, 0, OPTIONAL, "Step-size for search/marginalization over transient-window start-time, in seconds [Default:TAtom]" );
574
575 /* misc params */
576 XLALRegisterUvarMember( computeFtotal, BOOLEAN, 0, OPTIONAL, "Also compute 'total' F-statistic over the full data-span" );
577
578 XLALRegisterUvarMember( numDraws, INT4, 'N', OPTIONAL, "Number of random 'draws' to simulate" );
579 XLALRegisterUvarMember( randSeed, INT4, 0, OPTIONAL, "GSL random-number generator seed value to use" );
580
581 XLALRegisterUvarMember( outputStats, STRING, 'o', OPTIONAL, "Output file containing 'numDraws' random draws of stats" );
582 XLALRegisterUvarMember( outputAtoms, STRING, 0, OPTIONAL, "Output F-statistic atoms into a file with this basename" );
583 XLALRegisterUvarMember( outputFstatMap, STRING, 0, OPTIONAL, "Output F-statistic over 2D parameter space {t0, tau} into file with this basename" );
584
585 XLALRegisterUvarMember( outputInjParams, STRING, 0, OPTIONAL, "Output injection parameters into this file" );
586 XLALRegisterUvarMember( outputPosteriors, STRING, 0, OPTIONAL, "output posterior pdfs on t0 and tau (in octave format) into this file " );
587
588 XLALRegisterUvarMember( SignalOnly, BOOLEAN, 'S', OPTIONAL, "Signal only: generate pure signal without noise" );
589 XLALRegisterUvarMember( useFReg, BOOLEAN, 0, OPTIONAL, "use 'regularized' Fstat (1/D)*e^F (if TRUE) for marginalization, or 'standard' e^F (if FALSE)" );
590
591 XLALRegisterUvarMember( ephemEarth, STRING, 0, OPTIONAL, "Earth ephemeris file to use" );
592 XLALRegisterUvarMember( ephemSun, STRING, 0, OPTIONAL, "Sun ephemeris file to use" );
593
594 /* 'hidden' stuff */
595 XLALRegisterUvarMember( TAtom, INT4, 0, DEVELOPER, "Time baseline for Fstat-atoms (typically Tsft) in seconds." );
596
597
598 if ( xlalErrno ) {
599 XLALPrintError( "%s: something failed in initializing user variabels .. errno = %d.\n", __func__, xlalErrno );
601 }
602
603 return XLAL_SUCCESS;
604
605} /* XLALInitUserVars() */
606
607
608/** Initialize Fstat-code: handle user-input and set everything up. */
609int
611{
612 /* generate log-string for file-output, containing cmdline-options + code VCS version info */
613 char *vcs;
614 if ( ( vcs = XLALVCSInfoString( lalPulsarVCSInfoList, 0, "%% " ) ) == NULL ) { /* short VCS version string */
615 XLALPrintError( "%s: XLALVCSInfoString failed with errno=%d.\n", __func__, xlalErrno );
617 }
618 char *cmdline;
619 if ( ( cmdline = XLALUserVarGetLog( UVAR_LOGFMT_CMDLINE ) ) == NULL ) {
620 XLALPrintError( "%s: XLALUserVarGetLog ( UVAR_LOGFMT_CMDLINE ) failed with errno=%d.\n", __func__, xlalErrno );
622 }
623 const char fmt[] = "%%%% cmdline: %s\n%%%%\n%s%%%%\n";
624 UINT4 len = strlen( vcs ) + strlen( cmdline ) + strlen( fmt ) + 1;
625 if ( ( cfg->logString = XLALMalloc( len ) ) == NULL ) {
626 XLALPrintError( "%s: XLALMalloc ( %d ) failed.\n", __func__, len );
628 }
629 sprintf( cfg->logString, fmt, cmdline, vcs );
630 XLALFree( cmdline );
631 XLALFree( vcs );
632
633 /* trivial settings from user-input */
634 cfg->SignalOnly = uvar->SignalOnly;
635
636 /* ----- parse user-input on signal amplitude-paramters + ranges ----- */
637 /* skypos */
638 cfg->skypos.longitude = uvar->Alpha; /* Alpha < 0 indicates 'allsky' */
639 cfg->skypos.latitude = uvar->Delta;
641
642 /* ----- amplitude-params: create prior pdfs reflecting the user-input */
643 if ( XLALInitAmplitudePrior( &cfg->AmpPrior, uvar ) != XLAL_SUCCESS ) {
645 }
646
647 /* ----- initialize random-number generator ----- */
648 /* read out environment variables GSL_RNG_xxx
649 * GSL_RNG_SEED: use to set random seed: default = 0, override by using --randSeed on cmdline
650 * GSL_RNG_TYPE: type of random-number generator to use: default = 'mt19937'
651 */
652 gsl_rng_env_setup();
653 /* allow overriding the random-seed per command-line */
654 if ( XLALUserVarWasSet( &uvar->randSeed ) ) {
655 gsl_rng_default_seed = uvar->randSeed;
656 }
657 cfg->rng = gsl_rng_alloc( gsl_rng_default );
658
659 LogPrintf( LOG_DEBUG, "random-number generator type: %s\n", gsl_rng_name( cfg->rng ) );
660 LogPrintf( LOG_DEBUG, "seed = %lu\n", gsl_rng_default_seed );
661
662 /* init ephemeris-data */
663 EphemerisData *edat = XLALInitBarycenter( uvar->ephemEarth, uvar->ephemSun );
664 if ( !edat ) {
665 LogPrintf( LOG_CRITICAL, "%s: XLALInitBarycenter failed: could not load Earth ephemeris '%s' and Sun ephemeris '%s'\n", __func__, uvar->ephemEarth, uvar->ephemSun );
667 }
668
669 /* init detector info */
670 if ( XLALUserVarWasSet( &uvar->IFO ) && XLALUserVarWasSet( &uvar->IFOs ) ) {
671 XLALPrintError( "%s: Please do not use both the deprecated --IFO and the newer --IFOs option.\n", __func__ );
673 } else if ( XLALUserVarWasSet( &uvar->IFO ) ) {
674 strcpy( uvar->IFOs->data[0], uvar->IFO );
675 }
676 UINT4 numDetectors = uvar->IFOs->length;
677 MultiLALDetector multiDet;
678 XLAL_CHECK( XLALParseMultiLALDetector( &multiDet, uvar->IFOs ) == XLAL_SUCCESS, XLAL_EFUNC );
679
680 /* TIMESTAMPS: either from --startTime+duration OR --timestampsFiles */
681 BOOLEAN have_startTime = XLALUserVarWasSet( &uvar->dataStartGPS );
682 BOOLEAN have_duration = XLALUserVarWasSet( &uvar->dataDuration );
683 BOOLEAN have_timestampsFiles = XLALUserVarWasSet( &uvar->timestampsFiles );
684 XLAL_CHECK( ( have_duration && have_startTime ) || !( have_duration || have_startTime ), XLAL_EINVAL, "Need BOTH {--dataStartGPS,--dataDuration} or NONE\n" );
685 // at least one of {startTime,timestamps} required
686 XLAL_CHECK( have_timestampsFiles || have_startTime, XLAL_EINVAL, "Need either --dataStartGPS and --dataDuration, OR --timestampsFiles}\n" );
687 // don't allow timestamps + {startTime+duration}
688 XLAL_CHECK( !have_timestampsFiles || !have_startTime, XLAL_EINVAL, "--timestampsFiles incompatible with --dataStartGPS and --dataDuration}\n" );
690 UINT4 numSteps = 0;
691 INT4 dataStartGPS = LAL_INT4_MAX;
692 INT4 dataEndGPS = 0;
693 INT4 dataDuration = 0;
694 /* now handle the two mutually-exclusive cases */
695 if ( have_timestampsFiles ) {
696 /* NOTE: nanoseconds will be truncated! */
697 XLAL_CHECK( ( multiTS = XLALReadMultiTimestampsFiles( uvar->timestampsFiles ) ) != NULL, XLAL_EFUNC );
698 XLAL_CHECK( ( multiTS->length > 0 ) && ( multiTS->data != NULL ), XLAL_EINVAL, "Got empty timestamps-list from XLALReadMultiTimestampsFiles()" );
699 XLAL_CHECK( ( multiTS->length == numDetectors ), XLAL_EINVAL, "Number of timestamps files does not match number of IFOs! (%d!=%d)", multiTS->length, numDetectors );
700 for ( UINT4 X = 0; X < numDetectors; X++ ) {
701 XLAL_CHECK( ( multiTS->data[X]->length > 0 ) && ( multiTS->data[X]->data != NULL ), XLAL_EINVAL, "Got empty timestamps-list for detector %s", uvar->IFOs->data[X] );
702 numSteps = multiTS->data[X]->length;
703 multiTS->data[X]->deltaT = uvar->TAtom; // Tsft information not given by timestamps-file
704 dataStartGPS = min( dataStartGPS, multiTS->data[X]->data[0].gpsSeconds );
705 dataEndGPS = max( dataEndGPS, multiTS->data[X]->data[numSteps - 1].gpsSeconds );
706 }
707 dataDuration = uvar->TAtom + dataEndGPS - dataStartGPS; // difference of last and first timestamp, plus one extra atom length
708 } // endif have_timestampsFiles
709 else {
710 /* init timestamps vector covering observation time */
711 XLAL_CHECK( ( multiTS = XLALCreateMultiLIGOTimeGPSVector( numDetectors ) ) != NULL, XLAL_EINVAL, "XLALCreateMultiLIGOTimeGPSVector(%d) failed.", numDetectors );
712 dataStartGPS = uvar->dataStartGPS;
713 dataDuration = uvar->dataDuration;
714 numSteps = ( UINT4 ) ceil( dataDuration / uvar->TAtom );
715 XLAL_CHECK( numSteps >= 2, XLAL_EINVAL, "Need timestamps covering at least 2 atoms!" );
716 for ( UINT4 X = 0; X < numDetectors; X++ ) {
717 XLAL_CHECK( ( multiTS->data[X] = XLALCreateTimestampVector( numSteps ) ) != NULL, XLAL_EFUNC, "XLALCreateTimestampVector(%d) failed.", numSteps );
718 multiTS->data[X]->deltaT = uvar->TAtom;
719 for ( UINT4 i = 0; i < numSteps; i ++ ) {
720 UINT4 ti = dataStartGPS + i * uvar->TAtom;
721 multiTS->data[X]->data[i].gpsSeconds = ti;
722 multiTS->data[X]->data[i].gpsNanoSeconds = 0;
723 }
724 }
725 } // endif !have_noiseSFTs
726
727 /* get detector states */
728 if ( ( cfg->multiDetStates = XLALGetMultiDetectorStates( multiTS, &multiDet, edat, 0.5 * uvar->TAtom ) ) == NULL ) {
729 XLALPrintError( "%s: XLALGetMultiDetectorStates() failed.\n", __func__ );
731 }
732
733 /* get rid of all temporary memory allocated for this step */
736 multiTS = NULL;
737
738
739 /* ---------- initialize transient window ranges, for injection ... ---------- */
741 int twtype;
742 XLAL_CHECK( ( twtype = XLALParseTransientWindowName( uvar->injectWindow_type ) ) >= 0, XLAL_EFUNC );
743 InjectRange.type = twtype;
744
745 /* make sure user doesn't set window=none but sets window-parameters => indicates she didn't mean 'none' */
746 if ( InjectRange.type == TRANSIENT_NONE ) {
747 if ( XLALUserVarWasSet( &uvar->injectWindow_t0Days ) || XLALUserVarWasSet( &uvar->injectWindow_t0DaysBand ) ||
748 XLALUserVarWasSet( &uvar->injectWindow_tauDays ) || XLALUserVarWasSet( &uvar->injectWindow_tauDaysBand ) ) {
749 XLALPrintError( "%s: ERROR: injectWindow_type == NONE, but window-parameters were set! Use a different window-type!\n", __func__ );
751 }
752 }
753
754 if ( uvar->injectWindow_t0DaysBand < 0 || uvar->injectWindow_tauDaysBand < 0 ) {
755 XLALPrintError( "%s: only positive t0/tau window injection bands allowed (%f, %f)\n", __func__, uvar->injectWindow_t0DaysBand, uvar->injectWindow_tauDaysBand );
757 }
758
759 /* apply correct defaults if unset: t0=dataStart, t0Band=dataDuration-3*tauMax */
760 InjectRange.t0 = dataStartGPS + uvar->injectWindow_t0Days * DAY24;
761
762 REAL8 tauMax = ( uvar->injectWindow_tauDays + uvar->injectWindow_tauDaysBand ) * DAY24;
763 if ( XLALUserVarWasSet( &uvar->injectWindow_t0DaysBand ) ) {
764 InjectRange.t0Band = uvar->injectWindow_t0DaysBand * DAY24;
765 } else {
766 InjectRange.t0Band = fmax( 0.0, dataDuration - TRANSIENT_EXP_EFOLDING * tauMax - InjectRange.t0 ); /* make sure it's >= 0 */
767 }
768
769 InjectRange.tau = ( UINT4 )( uvar->injectWindow_tauDays * DAY24 );
770 InjectRange.tauBand = ( UINT4 )( uvar->injectWindow_tauDaysBand * DAY24 );
771
772 cfg->transientInjectRange = InjectRange;
773
774 /* ---------- ... and for search -------------------- */
776 XLAL_CHECK( ( twtype = XLALParseTransientWindowName( uvar->searchWindow_type ) ) >= 0, XLAL_EFUNC );
777 SearchRange.type = twtype;
778
779 /* apply correct defaults if unset: use inect window */
780 if ( !XLALUserVarWasSet( &uvar->searchWindow_type ) ) {
781 SearchRange.type = InjectRange.type;
782 }
783 if ( !XLALUserVarWasSet( &uvar->searchWindow_t0Days ) ) {
784 SearchRange.t0 = InjectRange.t0;
785 } else {
786 SearchRange.t0 = dataStartGPS + uvar->searchWindow_t0Days * DAY24;
787 }
788 if ( !XLALUserVarWasSet( &uvar->searchWindow_t0DaysBand ) ) {
789 SearchRange.t0Band = InjectRange.t0Band;
790 } else {
791 SearchRange.t0Band = ( UINT4 )( uvar->searchWindow_t0DaysBand * DAY24 );
792 }
793 if ( !XLALUserVarWasSet( &uvar->searchWindow_tauDays ) ) {
794 SearchRange.tau = InjectRange.tau;
795 } else {
796 SearchRange.tau = ( UINT4 )( uvar->searchWindow_tauDays * DAY24 );
797 }
798 if ( !XLALUserVarWasSet( &uvar->searchWindow_tauDaysBand ) ) {
799 SearchRange.tauBand = InjectRange.tauBand;
800 } else {
801 SearchRange.tauBand = ( UINT4 )( uvar->searchWindow_tauDaysBand * DAY24 );
802 }
803
804 if ( XLALUserVarWasSet( &uvar->searchWindow_dt0 ) ) {
805 SearchRange.dt0 = uvar->searchWindow_dt0;
806 } else {
807 SearchRange.dt0 = uvar->TAtom;
808 }
809
810 if ( XLALUserVarWasSet( &uvar->searchWindow_dtau ) ) {
811 SearchRange.dtau = uvar->searchWindow_dtau;
812 } else {
813 SearchRange.dtau = uvar->TAtom;
814 }
815
816 /* make sure user doesn't set window=none but sets window-parameters => indicates she didn't mean 'none' */
817 if ( SearchRange.type == TRANSIENT_NONE )
818 if ( XLALUserVarWasSet( &uvar->searchWindow_t0Days ) || XLALUserVarWasSet( &uvar->searchWindow_t0DaysBand ) ||
819 XLALUserVarWasSet( &uvar->searchWindow_tauDays ) || XLALUserVarWasSet( &uvar->searchWindow_tauDaysBand ) ) {
820 XLALPrintError( "%s: ERROR: searchWindow_type == NONE, but window-parameters were set! Use a different window-type!\n", __func__ );
822 }
823
824 if ( uvar->searchWindow_t0DaysBand < 0 || uvar->searchWindow_tauDaysBand < 0 ) {
825 XLALPrintError( "%s: only positive t0/tau window injection bands allowed (%f, %f)\n", __func__, uvar->searchWindow_t0DaysBand, uvar->searchWindow_tauDaysBand );
827 }
828
829 cfg->transientSearchRange = SearchRange;
830
831 return XLAL_SUCCESS;
832
833} /* XLALInitCode() */
834
835
836/** Initialize amplitude-prior pdfs from the user-input
837 */
838int
840{
841 const UINT4 AmpPriorBins = 100; // defines the binnning accuracy of our prior-pdfs
842
843 /* consistency check */
844 if ( !AmpPrior || !uvar ) {
845 XLALPrintError( "%s: invalid NULL input 'AmpPrior' or 'uvar'\n", __func__ );
847 }
848 if ( AmpPrior->pdf_h0Nat || AmpPrior->pdf_cosi || AmpPrior->pdf_psi || AmpPrior->pdf_phi0 ) {
849 XLALPrintError( "%s: AmplitudePriors must be set to NULL before calling this function!\n", __func__ );
851 }
852
853 /* first check that user only provided *one* method of determining the amplitude-prior range */
854 UINT4 numSets = 0;
855 if ( uvar->fixedh0Nat >= 0 ) {
856 numSets ++;
857 }
858 if ( uvar->fixedSNR >= 0 ) {
859 numSets ++;
860 }
861 if ( uvar->fixedh0NatMax >= 0 ) {
862 numSets ++ ;
863 }
864 if ( uvar->fixedRhohMax >= 0 ) {
865 numSets ++;
866 }
867 if ( numSets != 1 ) {
868 XLALPrintError( "%s: Specify (>=0) exactly *ONE* amplitude-prior range of {fixedh0Nat, fixedSNR, fixedh0NatMax, fixedRhohMax}\n", __func__ );
870 }
871
872 /* ===== first pass: deal with all user-supplied fixed values ==> singular priors! ===== */
873
874 /* ----- h0 ----- */
875 if ( uvar->fixedh0Nat >= 0 ) /* fix h0Nat */
876 if ( ( AmpPrior->pdf_h0Nat = XLALCreateSingularPDF1D( uvar->fixedh0Nat ) ) == NULL ) {
878 }
879
880 if ( uvar->fixedSNR >= 0 ) /* dummy-pdf, as signal will be computed with h0Nat=1 then rescaled to fixedSNR */
881 if ( ( AmpPrior->pdf_h0Nat = XLALCreateSingularPDF1D( 1.0 ) ) == NULL ) {
883 }
884
885 AmpPrior->fixedSNR = uvar->fixedSNR;
886 AmpPrior->fixRhohMax = ( uvar->fixedRhohMax >= 0 );
887
888 /* ----- cosi ----- */
889 if ( XLALUserVarWasSet( &uvar->cosi ) )
890 if ( ( AmpPrior->pdf_cosi = XLALCreateSingularPDF1D( uvar->cosi ) ) == NULL ) {
892 }
893 /* ----- psi ----- */
894 if ( XLALUserVarWasSet( &uvar->psi ) )
895 if ( ( AmpPrior->pdf_psi = XLALCreateSingularPDF1D( uvar->psi ) ) == NULL ) {
897 }
898 /* ----- phi0 ----- */
899 if ( XLALUserVarWasSet( &uvar->phi0 ) )
900 if ( ( AmpPrior->pdf_phi0 = XLALCreateSingularPDF1D( uvar->phi0 ) ) == NULL ) {
902 }
903
904
905 /* ===== second pass: deal with non-singular prior ranges, taking into account the type of priors to use */
906 REAL8 h0NatMax = 0;
907 if ( uvar->fixedh0NatMax >= 0 ) { /* draw h0Nat from [0, h0NatMax] */
908 h0NatMax = uvar->fixedh0NatMax;
909 }
910 if ( uvar->fixedRhohMax >= 0 ) { /* draw h0 from [0, rhohMax/(detM)^(1/8)] */
911 h0NatMax = uvar->fixedRhohMax; /* at first, will be rescaled by (detM)^(1/8) after the fact */
912 }
913
914 switch ( uvar->AmpPriorType ) {
916
917 /* ----- h0 ----- */ // uniform in [0, h0NatMax] : not that 'physical', but simple
918 if ( AmpPrior->pdf_h0Nat == NULL )
919 if ( ( AmpPrior->pdf_h0Nat = XLALCreateUniformPDF1D( 0, h0NatMax ) ) == NULL ) {
921 }
922 /* ----- cosi ----- */
923 if ( AmpPrior->pdf_cosi == NULL )
924 if ( ( AmpPrior->pdf_cosi = XLALCreateUniformPDF1D( -1.0, 1.0 ) ) == NULL ) {
926 }
927 /* ----- psi ----- */
928 if ( AmpPrior->pdf_psi == NULL )
929 if ( ( AmpPrior->pdf_psi = XLALCreateUniformPDF1D( -LAL_PI_4, LAL_PI_4 ) ) == NULL ) {
931 }
932 /* ----- phi0 ----- */
933 if ( AmpPrior->pdf_phi0 == NULL )
934 if ( ( AmpPrior->pdf_phi0 = XLALCreateUniformPDF1D( 0, LAL_TWOPI ) ) == NULL ) {
936 }
937
938 break;
939
941 /* ----- pdf(h0) ~ h0^3 ----- */
942 if ( AmpPrior->pdf_h0Nat == NULL ) {
943 UINT4 i;
944 pdf1D_t *pdf;
945 if ( ( pdf = XLALCreateDiscretePDF1D( 0, h0NatMax, AmpPriorBins ) ) == NULL ) {
947 }
948
949 for ( i = 0; i < pdf->probDens->length; i ++ ) {
950 REAL8 xMid = 0.5 * ( pdf->xTics->data[i] + pdf->xTics->data[i + 1] );
951 pdf->probDens->data[i] = CUBE( xMid ); // pdf(h0) ~ h0^3
952 }
953 AmpPrior->pdf_h0Nat = pdf;
954 }
955 /* ----- pdf(cosi) ~ ( 1 - cosi^2)^3 ----- */
956 if ( AmpPrior->pdf_cosi == NULL ) {
957 UINT4 i;
958 pdf1D_t *pdf;
959 if ( ( pdf = XLALCreateDiscretePDF1D( -1.0, 1.0, AmpPriorBins ) ) == NULL ) {
961 }
962
963 for ( i = 0; i < pdf->probDens->length; i ++ ) {
964 REAL8 xMid = 0.5 * ( pdf->xTics->data[i] + pdf->xTics->data[i + 1] );
965 REAL8 y = 1.0 - SQ( xMid );
966 pdf->probDens->data[i] = CUBE( y );
967 }
968 AmpPrior->pdf_cosi = pdf;
969 }
970 /* ----- psi ----- */
971 if ( AmpPrior->pdf_psi == NULL )
972 if ( ( AmpPrior->pdf_psi = XLALCreateUniformPDF1D( -LAL_PI_4, LAL_PI_4 ) ) == NULL ) {
974 }
975 /* ----- phi0 ----- */
976 if ( AmpPrior->pdf_phi0 == NULL )
977 if ( ( AmpPrior->pdf_phi0 = XLALCreateUniformPDF1D( 0, LAL_TWOPI ) ) == NULL ) {
979 }
980
981 break;
982
983 default:
984 XLALPrintError( "%s: something went wrong ... unknown priorType = %d\n", __func__, uvar->AmpPriorType );
986 break;
987
988 } // switch( uvar->AmpPriorType )
989
990 return XLAL_SUCCESS;
991
992} /* XLALInitAmplitudePrior() */
#define __func__
log an I/O error, i.e.
#define LAL_INT4_MAX
#define IFO
void LALCheckMemoryLeaks(void)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
#define STRING(a)
int XLALOutputPDF1D_to_fp(FILE *fp, const pdf1D_t *pdf, const char *name)
Function to write a pdf1D to given filepointer fp.
pdf1D_t * XLALCreateSingularPDF1D(REAL8 x0)
Creator function for a 'singular' 1D pdf, containing a single value with certainty,...
pdf1D_t * XLALCreateDiscretePDF1D(REAL8 xMin, REAL8 xMax, UINT4 numBins)
Creator function for a generic discrete 1D pdf over [xMin, xMax], discretized into numBins bins.
REAL8 XLALFindModeOfPDF1D(const pdf1D_t *pdf)
Find the 'mode' of the probabilty distribution, ie.
pdf1D_t * XLALCreateUniformPDF1D(REAL8 xMin, REAL8 xMax)
Creator function for a uniform 1D pdf over [xMin, xMax].
void XLALDestroyPDF1D(pdf1D_t *pdf)
Destructor function for 1-D pdf.
#define fprintf
void XLALDestroyMultiFstatAtomVector(MultiFstatAtomVector *multiAtoms)
Free all memory associated with a MultiFstatAtomVector structure.
Definition: ComputeFstat.c:338
int XLALParseMultiLALDetector(MultiLALDetector *detInfo, const LALStringVector *detNames)
Parse string-vectors (typically input by user) of N detector-names for detectors ,...
void XLALDestroyMultiDetectorStateSeries(MultiDetectorStateSeries *mdetStates)
Helper function to get rid of a multi-IFO DetectorStateSeries Note, this is "NULL-robust" in the sens...
MultiDetectorStateSeries * XLALGetMultiDetectorStates(const MultiLIGOTimeGPSVector *multiTS, const MultiLALDetector *multiIFO, const EphemerisData *edat, REAL8 tOffset)
Get the detector-time series for the given MultiLIGOTimeGPSVector.
int XLALFileClose(LALFILE *file)
LALFILE * XLALFileOpen(const char *path, const char *mode)
int XLALFilePrintf(LALFILE *file, const char *fmt,...)
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
void XLALDestroyMultiAMCoeffs(MultiAMCoeffs *multiAMcoef)
Destroy a MultiAMCoeffs structure.
Definition: LALComputeAM.c:469
#define LAL_YRSID_SI
#define LAL_TWOPI
#define LAL_PI_4
unsigned char BOOLEAN
double REAL8
#define XLAL_INIT_DECL(var,...)
char CHAR
uint32_t UINT4
int32_t INT4
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
void XLALFree(void *p)
char char * XLALStringDuplicate(const char *s)
char * XLALVCSInfoString(const LALVCSInfoList vcs_list, const int verbose, const char *prefix)
void void int XLALfprintfGSLmatrix(FILE *fp, const char *fmt, const gsl_matrix *gij) _LAL_GCC_VPRINTF_FORMAT_(2)
void LogPrintf(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
LOG_CRITICAL
LOG_DEBUG
@ TRANSIENT_NONE
Note: in this case the window-parameters will be ignored, and treated as rect={data},...
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:47
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
MultiLIGOTimeGPSVector * XLALCreateMultiLIGOTimeGPSVector(UINT4 numDetectors)
Simple creator function for MultiLIGOTimeGPSVector with numDetectors entries.
MultiLIGOTimeGPSVector * XLALReadMultiTimestampsFiles(const LALStringVector *fnames)
backwards compatible wrapper to XLALReadMultiTimestampsFilesConstrained() without GPS-time constraint...
COORDINATESYSTEM_EQUATORIAL
LALStringVector * XLALCreateStringVector(const CHAR *str1,...)
int write_InjParams_to_fp(FILE *fp, const InjParams_t *par, const UINT4 dataStartGPS, const BOOLEAN outputMmunuX, const UINT4 numDetectors)
Write an injection-parameters structure to the given file-pointer, adding one line with the injection...
MultiFstatAtomVector * XLALSynthesizeTransientAtoms(InjParams_t *injParamsOut, SkyPosition skypos, AmplitudePrior_t AmpPrior, transientWindowRange_t transientInjectRange, const MultiDetectorStateSeries *multiDetStates, BOOLEAN SignalOnly, multiAMBuffer_t *multiAMBuffer, gsl_rng *rng, INT4 lineX, const MultiNoiseWeights *multiNoiseWeights)
Generates a multi-Fstat atoms vector for given parameters, drawing random parameters wherever require...
@ AMP_PRIOR_TYPE_CANONICAL
'canonical' priors: uniform in A^mu up to h_max
@ AMP_PRIOR_TYPE_PHYSICAL
'physical' priors: isotropic pdf{cosi,psi,phi0} AND flat pdf(h0)
transientFstatMap_t * XLALComputeTransientFstatMap(const MultiFstatAtomVector *multiFstatAtoms, transientWindowRange_t windowRange, BOOLEAN useFReg)
Function to compute transient-window "F-statistic map" over start-time and timescale {t0,...
void XLALDestroyTransientFstatMap(transientFstatMap_t *FstatMap)
Standard destructor for transientFstatMap_t Fully NULL-robust as usual.
int XLALParseTransientWindowName(const char *windowName)
Parse a transient window name string into the corresponding transientWindowType.
int write_transientCandidate_to_fp(LALFILE *fp, const transientCandidate_t *thisCand, const char timeUnit)
Write one line for given transient CW candidate into output file.
#define DAY24
standard 24h day = 86400 seconds ==> this is what's used in the definition of 'tauDays'
void XLALDestroyExpLUT(void)
Destructor function for expLUT_t lookup table.
int write_MultiFstatAtoms_to_fp(LALFILE *fp, const MultiFstatAtomVector *multiAtoms)
Write multi-IFO F-stat atoms 'multiAtoms' into output stream 'fstat'.
#define TRANSIENT_EXP_EFOLDING
e-folding parameter for exponential window, after which we truncate the window for efficiency.
pdf1D_t * XLALComputeTransientPosterior_t0(transientWindowRange_t windowRange, const transientFstatMap_t *FstatMap)
Compute transient-CW posterior (normalized) on start-time t0, using given type and parameters of tran...
pdf1D_t * XLALComputeTransientPosterior_tau(transientWindowRange_t windowRange, const transientFstatMap_t *FstatMap)
Compute transient-CW posterior (normalized) on timescale tau, using given type and parameters of tran...
REAL8 XLALComputeTransientBstat(transientWindowRange_t windowRange, const transientFstatMap_t *FstatMap)
Compute transient-CW Bayes-factor B_SG = P(x|HypS)/P(x|HypG) (where HypG = Gaussian noise hypothesis)...
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
void XLALDestroyUserVars(void)
#define XLALRegisterUvarMember(name, type, option, category,...)
void CHAR * XLALUserVarGetLog(UserVarLogFormat format)
int XLALUserVarWasSet(const void *cvar)
UVAR_LOGFMT_CMDLINE
#define xlalErrno
#define XLAL_ERROR(...)
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EIO
XLAL_EINVAL
list y
size_t numDetectors
Signal (amplitude) parameter ranges.
pdf1D_t * pdf_h0Nat
pdf for h0/sqrt{Sn}
REAL8 fixedSNR
alternative 1: adjust h0 to fix the optimal SNR of the signal
BOOLEAN fixRhohMax
alternative 2: draw h0 with fixed rhohMax = h0Max * (detM)^(1/8) <==> canonical Fstat prior
pdf1D_t * pdf_phi0
pdf(phi0)
pdf1D_t * pdf_cosi
pdf(cosi)
pdf1D_t * pdf_psi
pdf(psi)
Configuration settings required for and defining a coherent pulsar search.
gsl_rng * rng
gsl random-number generator
AmplitudePrior_t AmpPrior
amplitude-parameter priors to draw signals from
SkyPosition skypos
(Alpha,Delta,system).
CHAR * logString
logstring for file-output, containing cmdline-options + code VCS version info
BOOLEAN SignalOnly
dont generate noise-draws: will result in non-random 'signal only' values of F and B
MultiDetectorStateSeries * multiDetStates
multi-detector state series covering observation time
transientWindowRange_t transientSearchRange
transient-window range for the search (flat priors)
transientWindowRange_t transientInjectRange
transient-window range for injections (flat priors)
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
Hold all (generally) randomly drawn injection parameters: skypos, amplitude-params,...
Multi-IFO time-series of DetectorStates.
A multi-detector vector of FstatAtomVector.
Definition: ComputeFstat.h:187
array of detectors definitions 'LALDetector'
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
Definition: SFTfileIO.h:198
REAL8 longitude
REAL8 latitude
CoordinateSystem system
struct for buffering of AM-coeffs, if signal for same sky-position is injected
Struct holding a transient CW candidate.
Struct holding a transient-window "F-statistic map" over start-time and timescale {t0,...
REAL8 maxF
maximal F-value obtained over transientWindowRange
Struct defining a range of transient windows.
int main(int argc, char *argv[])
MAIN function Generates samples of B-stat and F-stat according to their pdfs for given signal-params.
int XLALInitCode(ConfigVariables *cfg, const UserInput_t *uvar)
Initialize Fstat-code: handle user-input and set everything up.
int XLALInitAmplitudePrior(AmplitudePrior_t *AmpPrior, const UserInput_t *uvar)
Initialize amplitude-prior pdfs from the user-input.
int vrbflg
defined in lal/lib/std/LALError.c
#define DEFAULT_TRANSIENT
#define CUBE(x)
#define SQ(x)
#define min(a, b)
#define DEFAULT_IFO
int XLALInitUserVars(UserInput_t *uvar)
Register all our "user-variables" that can be specified from cmd-line and/or config-file.
#define max(a, b)