Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-8a6b96f
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
MCInjectHoughMulti.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2005 Alicia Sintes, Badri Krishnan,
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 * \file
22 * \ingroup lalpulsar_bin_Hough
23 * \author Alicia Sintes, Badri Krishnan
24 * \brief
25 * Monte Carlo signal injections for several h_0 values and
26 * compute the Hough transform for a small number of point in parameter space each time
27 */
28
29/*
30 The idea is that we would like to analize a 300 Hz band on a cluster of
31 machines. Each process should analyze 1 Hz band (or whatever).
32
33 - Read the band to be analized and the wings needed to read the originals SFTs.
34 -Read the h_0 values to be analyzed in one go
35 -loop over the MC injections:
36 + Generate random parameters (f, f', alpha, delata, i...)
37 + generate h(t), produce its FFT
38 + Add h(f) to SFT for a given h_o value (and all of them)
39 + get number count
40
41 Input shoud be from
42 SFT files
43 band, nh_0, h_01, h_02....
44 ephemeris info
45
46 This code will output files containing the MC results and info about injected
47 signals.
48*/
49
50#include "MCInjectHoughMulti.h"
51
52
53#define EARTHEPHEMERIS "./earth05-09.dat"
54#define SUNEPHEMERIS "./sun05-09.dat"
55
56#define MAXFILENAMELENGTH 512 /* maximum # of characters of a filename */
57
58#define ACCURACY 0.00000001 /* of the velocity calculation */
59#define MAXFILES 3000 /* maximum number of files to read in a directory */
60
61#define THRESHOLD 1.6 /* thresold for peak selection, with respect to the
62 the averaged power in the search band */
63#define F0 250.0 /* frequency to build the LUT and start search */
64#define FBAND 2.0 /* search frequency band (in Hz) */
65#define ALPHA 0.0 /* center of the sky patch (in radians) */
66#define DELTA (-LAL_PI_2)
67#define PATCHSIZEX (LAL_PI*0.99) /* patch size */
68#define PATCHSIZEY (LAL_PI*0.99)
69#define NFSIZE 21 /* n-freq. span of the cylinder, to account for spin-down
70 search */
71#define BLOCKSRNGMED 101 /* Running median window size */
72#define NH0 2 /* number of h0 values to be analyzed */
73#define H0MIN 1.0e-23
74#define H0MAX 1.0e-22
75#define NMCLOOP 2 /* number of Monte-Carlos */
76#define NTEMPLATES 16 /* number templates for each Monte-Carlo */
77#define DTERMS 5
78
79#define SFTDIRECTORY "/local_data/sintes/SFT-S5-120-130/*SFT*.*"
80/* */
81#define DIROUT "./" /* output directory */
82#define FILEOUT "/HoughMC" /* prefix file output */
83
84#define TRUE (1==1)
85#define FALSE (1==0)
86
87/******************************************/
88void ComputeFoft_NM( LALStatus *status,
89 REAL8Vector *foft,
90 HoughTemplate *pulsarTemplate,
91 REAL8Vector *timeDiffV,
92 REAL8Cart3CoorVector *velV );
93
94
95void PrintLogFile2( LALStatus *status,
96 CHAR *dir,
97 CHAR *basename,
98 LALStringVector *linefiles,
99 CHAR *executable );
100
101/******************************************/
102
103
104/* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv------------------------------------ */
105int main( int argc, char *argv[] )
106{
107
108 static LALStatus status;
109
110 EphemerisData *edat = NULL;
111 static REAL8Cart3CoorVector velV;
112 static LIGOTimeGPSVector timeV;
113 MultiLIGOTimeGPSVector *multiIniTimeV = NULL;
114 static REAL8Vector timeDiffV;
115 LIGOTimeGPS firstTimeStamp, lastTimeStamp;
116 REAL8 tObs;
117
118 static REAL8Vector foft;
119 static REAL8Vector foftV[NTEMPLATES];
120 static REAL8Vector h0V;
121
122 static HoughInjectParams injectPar;
123 static PulsarData pulsarInject;
124 static HoughTemplate pulsarTemplate;
125 static HoughNearTemplates closeTemplates;
126 SkyConstAndZeroPsiAMResponse *pSkyConstAndZeroPsiAMResponse = NULL;
127 SFTandSignalParams *pSFTandSignalParams = NULL;
128
129
130 /* standard pulsar sft types */
131 MultiSFTVector *inputSFTs = NULL;
132 MultiSFTVector *sumSFTs = NULL;
133 MultiSFTVector *signalSFTs = NULL;
134
135 /* information about all the ifos */
136 MultiDetectorStateSeries *mdetStates = NULL;
137 UINT4 numifo;
138 UINT4 binsSFT;
139
140 /* vector of weights */
141 REAL8Vector weightsV, weightsNoise;
142 REAL8Vector XLAL_INIT_DECL( weightsAM );
143 REAL8 alphaPeak, meanN, sigmaN, significance;
144
145 REAL4TimeSeries *signalTseries = NULL;
146 static PulsarSignalParams params;
147 static SFTParams sftParams;
148 static UCHARPeakGram pg1;
149
150 UINT4 msp = 1; /*number of spin-down parameters */
151 REAL8 numberCount, maxNumberCount;
152 REAL8 numberCountV[NTEMPLATES];
153 UINT4 nTemplates;
154
155 UINT4 mObsCoh;
156 REAL8 normalizeThr;
157 REAL8 timeBase, deltaF;
158 REAL8 h0scale;
159
160 INT4 sftFminBin;
161 REAL8 fHeterodyne;
162 REAL8 tSamplingRate;
163
164 INT4 MCloopId;
165 INT4 h0loop;
166
167 FILE *fpPar = NULL;
168 FILE *fpH0 = NULL;
169 FILE *fpNc = NULL;
170
171 /* sft constraint variables */
172 LIGOTimeGPS startTimeGPS, endTimeGPS;
173 LIGOTimeGPSVector *inputTimeStampsVector = NULL;
174
175 /******************************************************************/
176 /* user input variables */
177 /******************************************************************/
178 BOOLEAN uvar_weighAM, uvar_weighNoise, uvar_printLog, uvar_fast;
179 INT4 uvar_blocksRngMed, uvar_nh0, uvar_nMCloop, uvar_AllSkyFlag;
180 INT4 uvar_nfSizeCylinder, uvar_maxBinsClean, uvar_Dterms;
181 REAL8 uvar_f0, uvar_fSearchBand, uvar_peakThreshold, uvar_h0Min, uvar_h0Max, uvar_maxSpin, uvar_minSpin;
182 REAL8 uvar_alpha, uvar_delta, uvar_patchSizeAlpha, uvar_patchSizeDelta, uvar_pixelFactor;
183 CHAR *uvar_earthEphemeris = NULL;
184 CHAR *uvar_sunEphemeris = NULL;
185 CHAR *uvar_sftDir = NULL;
186 CHAR *uvar_dirnameOut = NULL;
187 CHAR *uvar_fnameOut = NULL;
188 LALStringVector *uvar_linefiles = NULL;
189 REAL8 uvar_startTime, uvar_endTime;
190 CHAR *uvar_timeStampsFile = NULL;
191 /******************************************************************/
192 /* set up the default parameters */
193 /******************************************************************/
194 /* LAL error-handler */
195 lal_errhandler = LAL_ERR_EXIT;
196
197 uvar_AllSkyFlag = 1;
198
199 uvar_weighAM = TRUE;
200 uvar_weighNoise = TRUE;
201 uvar_printLog = FALSE;
202 uvar_fast = TRUE;
203
204 uvar_nh0 = NH0;
205 uvar_h0Min = H0MIN;
206 uvar_h0Max = H0MAX;
207
208 uvar_nMCloop = NMCLOOP;
209 nTemplates = NTEMPLATES;
210 uvar_alpha = ALPHA;
211 uvar_delta = DELTA;
212 uvar_f0 = F0;
213 uvar_fSearchBand = FBAND;
214 uvar_peakThreshold = THRESHOLD;
215 uvar_nfSizeCylinder = NFSIZE;
216 uvar_blocksRngMed = BLOCKSRNGMED;
217 uvar_maxBinsClean = 100;
218 uvar_Dterms = DTERMS;
219
220 uvar_pixelFactor = PIXELFACTOR;
221 uvar_maxSpin = 1e-9;
222 uvar_minSpin = -1e-8;
223
224 uvar_patchSizeAlpha = PATCHSIZEX;
225 uvar_patchSizeDelta = PATCHSIZEY;
226
227 uvar_earthEphemeris = ( CHAR * )LALCalloc( MAXFILENAMELENGTH, sizeof( CHAR ) );
228 strcpy( uvar_earthEphemeris, EARTHEPHEMERIS );
229
230 uvar_sunEphemeris = ( CHAR * )LALCalloc( MAXFILENAMELENGTH, sizeof( CHAR ) );
231 strcpy( uvar_sunEphemeris, SUNEPHEMERIS );
232
233 uvar_sftDir = ( CHAR * )LALCalloc( MAXFILENAMELENGTH, sizeof( CHAR ) );
234 strcpy( uvar_sftDir, SFTDIRECTORY );
235
236 uvar_dirnameOut = ( CHAR * )LALCalloc( MAXFILENAMELENGTH, sizeof( CHAR ) );
237 strcpy( uvar_dirnameOut, DIROUT );
238
239 uvar_fnameOut = ( CHAR * )LALCalloc( MAXFILENAMELENGTH, sizeof( CHAR ) );
240 strcpy( uvar_fnameOut, FILEOUT );
241
242
243 /******************************************************************/
244 /* register user input variables */
245 /******************************************************************/
246 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_blocksRngMed, "blocksRngMed", INT4, 'w', OPTIONAL, "RngMed block size" ) == XLAL_SUCCESS, XLAL_EFUNC );
247 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f0, "f0", REAL8, 'f', OPTIONAL, "Start search frequency" ) == XLAL_SUCCESS, XLAL_EFUNC );
248 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_fSearchBand, "fSearchBand", REAL8, 'b', OPTIONAL, "Search frequency band" ) == XLAL_SUCCESS, XLAL_EFUNC );
249 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_peakThreshold, "peakThreshold", REAL8, 't', OPTIONAL, "Peak selection threshold" ) == XLAL_SUCCESS, XLAL_EFUNC );
250 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_earthEphemeris, "earthEphemeris", STRING, 'E', OPTIONAL, "Earth Ephemeris file" ) == XLAL_SUCCESS, XLAL_EFUNC );
251 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_sunEphemeris, "sunEphemeris", STRING, 'S', OPTIONAL, "Sun Ephemeris file" ) == XLAL_SUCCESS, XLAL_EFUNC );
252 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_sftDir, "sftDir", STRING, 'D', OPTIONAL, "SFT Directory. Possibilities are:\n"
253 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" ) == XLAL_SUCCESS, XLAL_EFUNC );
254 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dirnameOut, "dirnameOut", STRING, 'o', OPTIONAL, "Output directory" ) == XLAL_SUCCESS, XLAL_EFUNC );
255 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_fnameOut, "fnameout", STRING, '0', OPTIONAL, "Output file prefix" ) == XLAL_SUCCESS, XLAL_EFUNC );
256 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_alpha, "alpha", REAL8, 'r', OPTIONAL, "Right ascension" ) == XLAL_SUCCESS, XLAL_EFUNC );
257 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_delta, "delta", REAL8, 'l', OPTIONAL, "Declination" ) == XLAL_SUCCESS, XLAL_EFUNC );
258 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_patchSizeAlpha, "patchSizeAlpha", REAL8, 'R', OPTIONAL, "Patch size in right ascension" ) == XLAL_SUCCESS, XLAL_EFUNC );
259 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_patchSizeDelta, "patchSizeDelta", REAL8, 'L', OPTIONAL, "Patch size in declination" ) == XLAL_SUCCESS, XLAL_EFUNC );
260 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_AllSkyFlag, "patch", INT4, 'P', OPTIONAL, "Inject in patch if 0" ) == XLAL_SUCCESS, XLAL_EFUNC );
261 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nMCloop, "nMCloop", INT4, 'N', OPTIONAL, "Number of MC injections" ) == XLAL_SUCCESS, XLAL_EFUNC );
262 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_h0Min, "h0Min", REAL8, 'm', OPTIONAL, "Smallest h0 to inject" ) == XLAL_SUCCESS, XLAL_EFUNC );
263 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_h0Max, "h0Max", REAL8, 'M', OPTIONAL, "Largest h0 to inject" ) == XLAL_SUCCESS, XLAL_EFUNC );
264 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nh0, "nh0", INT4, 'n', OPTIONAL, "Number of h0 values to inject" ) == XLAL_SUCCESS, XLAL_EFUNC );
265 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_linefiles, "linefiles", STRINGVector, 0, OPTIONAL, "list of linefiles separated by commas" ) == XLAL_SUCCESS, XLAL_EFUNC );
266 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_weighAM, "weighAM", BOOLEAN, 0, OPTIONAL, "Use amplitude modulation weights" ) == XLAL_SUCCESS, XLAL_EFUNC );
267 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_weighNoise, "weighNoise", BOOLEAN, 0, OPTIONAL, "Use SFT noise weights" ) == XLAL_SUCCESS, XLAL_EFUNC );
268 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_printLog, "printLog", BOOLEAN, 0, OPTIONAL, "Print Log file" ) == XLAL_SUCCESS, XLAL_EFUNC );
269 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_fast, "fast", BOOLEAN, 0, OPTIONAL, "Use fast frequency domain SFT injections" ) == XLAL_SUCCESS, XLAL_EFUNC );
270 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_startTime, "startTime", REAL8, 0, OPTIONAL, "GPS start time of observation (SFT timestamps must be >= this)" ) == XLAL_SUCCESS, XLAL_EFUNC );
271 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_endTime, "endTime", REAL8, 0, OPTIONAL, "GPS end time of observation (SFT timestamps must be < this)" ) == XLAL_SUCCESS, XLAL_EFUNC );
272 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_timeStampsFile, "timeStampsFile", STRING, 0, OPTIONAL, "Input time-stamps file" ) == XLAL_SUCCESS, XLAL_EFUNC );
273 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nfSizeCylinder, "nfSizeCylinder", INT4, 0, OPTIONAL, "Size of cylinder of PHMDs" ) == XLAL_SUCCESS, XLAL_EFUNC );
274 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_maxBinsClean, "maxBinsClean", INT4, 0, OPTIONAL, "Maximum number of bins in cleaning" ) == XLAL_SUCCESS, XLAL_EFUNC );
275 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_Dterms, "Dterms", INT4, 0, OPTIONAL, "Number of f-bins in MC injection" ) == XLAL_SUCCESS, XLAL_EFUNC );
276 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_pixelFactor, "pixelfactor", REAL8, 0, OPTIONAL, "Pixelfactor used for sky resolution" ) == XLAL_SUCCESS, XLAL_EFUNC );
277 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_maxSpin, "MaxSpin", REAL8, 0, OPTIONAL, "Maximum Spin in PHMDs" ) == XLAL_SUCCESS, XLAL_EFUNC );
278 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_minSpin, "MinSpin", REAL8, 0, OPTIONAL, "Minimum Spin in PHMDs" ) == XLAL_SUCCESS, XLAL_EFUNC );
279
280 /******************************************************************/
281 /* read all command line variables */
282 /******************************************************************/
283 BOOLEAN should_exit = 0;
284 XLAL_CHECK_MAIN( XLALUserVarReadAllInput( &should_exit, argc, argv, lalPulsarVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
285 if ( should_exit ) {
286 exit( 1 );
287 }
288
289
290 /* very basic consistency checks on user input */
291 if ( uvar_f0 < 0 ) {
292 fprintf( stderr, "start frequency must be positive\n" );
293 exit( 1 );
294 }
295
296 if ( uvar_fSearchBand < 0 ) {
297 fprintf( stderr, "search frequency band must be positive\n" );
298 exit( 1 );
299 }
300
301 if ( uvar_peakThreshold < 0 ) {
302 fprintf( stderr, "peak selection threshold must be positive\n" );
303 exit( 1 );
304 }
305
306 if ( uvar_maxSpin < uvar_minSpin ) {
307 fprintf( stderr, "Maximum frequency derivative must be greater than minimum frequency derivative\n" );
308 exit( 1 );
309 }
310
311 LAL_CALL( LALRngMedBias( &status, &normalizeThr, uvar_blocksRngMed ), &status );
312
313 /******************************************************************/
314 /* write log file with command line arguments, cvs tags, and contents of skypatch file */
315 /******************************************************************/
316 if ( uvar_printLog ) {
317 LAL_CALL( PrintLogFile2( &status, uvar_dirnameOut, uvar_fnameOut, uvar_linefiles, argv[0] ), &status );
318 }
319
320 /******************************************************************/
321 /* set fullsky flag */
322 /******************************************************************/
323
324 injectPar.fullSky = 1;
325 if ( uvar_AllSkyFlag == 0 ) {
326 injectPar.fullSky = 0; /* patch case */
327 }
328
329 /******************************************************************/
330 /* computing h0 values */
331 /******************************************************************/
332 h0V.length = uvar_nh0;
333 h0V.data = NULL;
334 h0V.data = ( REAL8 * )LALMalloc( uvar_nh0 * sizeof( REAL8 ) );
335 h0V.data[0] = uvar_h0Min;
336
337 if ( uvar_nh0 > 1 ) {
338 INT4 k;
339 REAL8 steph0;
340 steph0 = ( uvar_h0Max - uvar_h0Min ) / ( uvar_nh0 - 1. );
341 for ( k = 1; k < uvar_nh0; ++k ) {
342 h0V.data[k] = h0V.data[k - 1] + steph0;
343 }
344 }
345
346 /******************************************************************/
347 /* preparing output files */
348 /******************************************************************/
349 {
350 INT4 k;
351 CHAR filename[MAXFILENAMELENGTH];
352
353 /* the paramerter file */
354 strcpy( filename, uvar_dirnameOut );
355
356 strcat( filename, uvar_fnameOut );
357 strcat( filename, "_par" );
358 fpPar = fopen( filename, "w" ); /* where to write the parameters */
359 /*setlinebuf(fpPar);*/ /* line buffered on */
360 setvbuf( fpPar, ( char * )NULL, _IOLBF, 0 );
361
362 /* the file with the h0 values */
363 strcpy( filename, uvar_dirnameOut );
364 strcat( filename, uvar_fnameOut );
365 strcat( filename, "_h0" );
366 fpH0 = fopen( filename, "w" ); /* where to write the parameters */
367 /*setlinebuf(fpH0); */ /* line buffered on */
368 setvbuf( fpH0, ( char * )NULL, _IOLBF, 0 );
369
370 /* the file with the the number-counts for different h0 values */
371 strcpy( filename, uvar_dirnameOut );
372 strcat( filename, uvar_fnameOut );
373 strcat( filename, "_nc" );
374 fpNc = fopen( filename, "w" ); /* where to write the parameters */
375 /*setlinebuf(fpNc);*/ /* line buffered on */
376 setvbuf( fpNc, ( char * )NULL, _IOLBF, 0 );
377
378 for ( k = 0; k < uvar_nh0; ++k ) {
379 fprintf( fpH0, "%g \n", h0V.data[k] );
380 }
381 fclose( fpH0 );
382
383 }
384
385 /******************************************************************/
386 /* sft reading */
387 /******************************************************************/
388
389 {
390 /* new SFT I/O data types */
391 SFTCatalog *catalog = NULL;
392 static SFTConstraints constraints;
393
394 REAL8 doppWings, f_min, f_max;
395
396 /* set detector constraint */
397 constraints.detector = NULL;
398
399 if ( XLALUserVarWasSet( &uvar_startTime ) ) {
400 XLALGPSSetREAL8( &startTimeGPS, uvar_startTime );
401 constraints.minStartTime = &startTimeGPS;
402 }
403
404 if ( XLALUserVarWasSet( &uvar_endTime ) ) {
405 XLALGPSSetREAL8( &endTimeGPS, uvar_endTime );
406 constraints.maxStartTime = &endTimeGPS;
407 }
408
409 if ( XLALUserVarWasSet( &uvar_timeStampsFile ) ) {
410 XLAL_CHECK_MAIN( ( inputTimeStampsVector = XLALReadTimestampsFile( uvar_timeStampsFile ) ) != NULL, XLAL_EFUNC );
411 constraints.timestamps = inputTimeStampsVector;
412 }
413
414 /* get sft catalog */
415 XLAL_CHECK_MAIN( ( catalog = XLALSFTdataFind( uvar_sftDir, &constraints ) ) != NULL, XLAL_EFUNC );
416 if ( ( catalog == NULL ) || ( catalog->length == 0 ) ) {
417 fprintf( stderr, "Unable to match any SFTs with pattern '%s'\n", uvar_sftDir );
418 exit( 1 );
419 }
420
421 /* now we can free the inputTimeStampsVector */
422 if ( XLALUserVarWasSet( &uvar_timeStampsFile ) ) {
423 XLALDestroyTimestampVector( inputTimeStampsVector );
424 }
425
426
427 /* get some sft parameters */
428 mObsCoh = catalog->length; /* number of sfts */
429 deltaF = catalog->data->header.deltaF; /* frequency resolution */
430 timeBase = 1.0 / deltaF; /* coherent integration time */
431
432 /* catalog is ordered in time so we can get start, end time and tObs*/
433 firstTimeStamp = catalog->data[0].header.epoch;
434 lastTimeStamp = catalog->data[mObsCoh - 1].header.epoch;
435 tObs = XLALGPSDiff( &lastTimeStamp, &firstTimeStamp ) + timeBase;
436
437 /* add wings for Doppler modulation and running median block size*/
438 doppWings = ( uvar_f0 + uvar_fSearchBand ) * VTOT;
439 f_min = uvar_f0 - doppWings - ( uvar_blocksRngMed + uvar_nfSizeCylinder ) * deltaF;
440 f_max = uvar_f0 + uvar_fSearchBand + doppWings + ( uvar_blocksRngMed + uvar_nfSizeCylinder ) * deltaF;
441
442 /* read sft files making sure to add extra bins for running median */
443 XLAL_CHECK_MAIN( ( inputSFTs = XLALLoadMultiSFTs( catalog, f_min, f_max ) ) != NULL, XLAL_EFUNC );
444
445 /* SFT info -- assume all SFTs have same length */
446 numifo = inputSFTs->length;
447 binsSFT = inputSFTs->data[0]->data->data->length;
448
449 /* some more sft parameetrs */
450 fHeterodyne = inputSFTs->data[0]->data[0].f0;
451 sftFminBin = ( INT4 ) floor( fHeterodyne * timeBase + 0.5 );
452 tSamplingRate = 2.0 * deltaF * ( binsSFT - 1. );
453
454 /* free memory */
455 XLALDestroySFTCatalog( catalog );
456 }
457
458 /******************************************************************/
459 /* allocate memory for sumSFTs of the same size of inputSFTs and place for
460 signal only sfts. */
461 /******************************************************************/
462
463 {
464 UINT4Vector numsft;
465 UINT4 iIFO;
466
467 numsft.length = numifo;
468 numsft.data = NULL;
469 numsft.data = ( UINT4 * )LALCalloc( numifo, sizeof( UINT4 ) );
470
471 for ( iIFO = 0; iIFO < numifo; iIFO++ ) {
472 numsft.data[iIFO] = inputSFTs->data[iIFO]->length;
473 }
474
475 XLAL_CHECK_MAIN( ( sumSFTs = XLALCreateMultiSFTVector( binsSFT, &numsft ) ) != NULL, XLAL_EFUNC );
476 XLAL_CHECK_MAIN( ( signalSFTs = XLALCreateMultiSFTVector( binsSFT, &numsft ) ) != NULL, XLAL_EFUNC );
477 LALFree( numsft.data );
478
479 }
480
481 /******************************************************************/
482 /* allocate memory for velocity vector and timestamps */
483 /******************************************************************/
484
485 velV.length = mObsCoh;
486 velV.data = NULL;
487 velV.data = ( REAL8Cart3Coor * )LALCalloc( mObsCoh, sizeof( REAL8Cart3Coor ) );
488
489 /* allocate memory for timestamps vector */
490 timeV.length = mObsCoh;
491 timeV.data = NULL;
492 timeV.data = ( LIGOTimeGPS * )LALCalloc( mObsCoh, sizeof( LIGOTimeGPS ) );
493
494 /* allocate memory for vector of time differences from start */
495 timeDiffV.length = mObsCoh;
496 timeDiffV.data = NULL;
497 timeDiffV.data = ( REAL8 * )LALCalloc( mObsCoh, sizeof( REAL8 ) );
498
499 /* start time of the sfts sorted for the different IFOs */
500 multiIniTimeV = ( MultiLIGOTimeGPSVector * )LALMalloc( sizeof( MultiLIGOTimeGPSVector ) );
501 multiIniTimeV->length = numifo;
502 multiIniTimeV->data = ( LIGOTimeGPSVector ** )LALCalloc( numifo, sizeof( LIGOTimeGPSVector * ) );
503
504 /******************************************************************/
505 /* get detector velocities and timestamps */
506 /******************************************************************/
507
508 {
509 UINT4 iIFO, iSFT, numsft, j;
510
511 XLAL_CHECK_MAIN( ( edat = XLALInitBarycenter( uvar_earthEphemeris, uvar_sunEphemeris ) ) != NULL, XLAL_EFUNC );
512
513 /* get information about all detectors including velocity and timestamps */
514 /* note that this function returns the velocity at the
515 mid-time of the SFTs --CAREFULL later on with the time stamps!!! velocity
516 is ok */
517
518 const REAL8 tOffset = 0.5 / inputSFTs->data[0]->data[0].deltaF;
519 XLAL_CHECK_MAIN( ( mdetStates = XLALGetMultiDetectorStatesFromMultiSFTs( inputSFTs, edat, tOffset ) ) != NULL, XLAL_EFUNC );
520
521 /* copy the timestamps and velocity vector */
522 for ( j = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
523 numsft = mdetStates->data[iIFO]->length;
524 for ( iSFT = 0; iSFT < numsft; iSFT++, j++ ) {
525 velV.data[j].x = mdetStates->data[iIFO]->data[iSFT].vDetector[0];
526 velV.data[j].y = mdetStates->data[iIFO]->data[iSFT].vDetector[1];
527 velV.data[j].z = mdetStates->data[iIFO]->data[iSFT].vDetector[2];
528 /* mid time of sfts */
529 timeV.data[j] = mdetStates->data[iIFO]->data[iSFT].tGPS;
530 } /* loop over SFTs */
531
532 multiIniTimeV->data[iIFO] = NULL;
533 XLAL_CHECK_MAIN( ( multiIniTimeV->data[iIFO] = XLALExtractTimestampsFromSFTs(
534 inputSFTs->data[iIFO] ) ) != NULL, XLAL_EFUNC );
535
536 } /* loop over IFOs */
537
538 /* compute the time difference relative to startTime for all SFT */
539 for ( j = 0; j < mObsCoh; j++ ) {
540 timeDiffV.data[j] = XLALGPSDiff( timeV.data + j, &firstTimeStamp );
541 }
542
543 /* removing mid time-stamps, no longer needed now */
544 LALFree( timeV.data );
545
546 }
547
548
549 /******************************************************************/
550 /* probability of selecting a peak and expected mean for noise only */
551 /******************************************************************/
552 alphaPeak = exp( - uvar_peakThreshold );
553 meanN = mObsCoh * alphaPeak;
554
555 /******************************************************************/
556 /* initialization of fast injections parameters TO BE FIXED for multiIFO
557 and memory should be dealocated at the end */
558 /******************************************************************/
559
560 pSkyConstAndZeroPsiAMResponse = ( SkyConstAndZeroPsiAMResponse * )
561 LALMalloc( sizeof( SkyConstAndZeroPsiAMResponse ) * numifo );
562 {
563 UINT4 numsft, iIFO;
564
565 for ( iIFO = 0; iIFO < numifo; iIFO++ ) {
566
567 numsft = mdetStates->data[iIFO]->length;
568
569 pSkyConstAndZeroPsiAMResponse[iIFO].skyConst = ( REAL8 * )LALMalloc( ( 2 * msp * ( numsft + 1 ) + 2 * numsft + 3 ) * sizeof( REAL8 ) );
570 pSkyConstAndZeroPsiAMResponse[iIFO].fPlusZeroPsi = ( REAL4 * )LALMalloc( numsft * sizeof( REAL4 ) );
571 pSkyConstAndZeroPsiAMResponse[iIFO].fCrossZeroPsi = ( REAL4 * )LALMalloc( numsft * sizeof( REAL4 ) );
572 }
573 }
574
575 {
576 INT4 k;
577
578 pSFTandSignalParams = ( SFTandSignalParams * )LALMalloc( sizeof( SFTandSignalParams ) );
579 /* create lookup table (LUT) values for doing trig */
580 /* pSFTandSignalParams->resTrig = 64; */ /* length sinVal and cosVal; resolution of trig functions = 2pi/resTrig */
581 pSFTandSignalParams->resTrig = 0; /* 08/02/04 gam; avoid serious bug when using LUT for trig calls */
582 pSFTandSignalParams->trigArg = ( REAL8 * )LALMalloc( ( pSFTandSignalParams->resTrig + 1 ) * sizeof( REAL8 ) );
583 pSFTandSignalParams->sinVal = ( REAL8 * )LALMalloc( ( pSFTandSignalParams->resTrig + 1 ) * sizeof( REAL8 ) );
584 pSFTandSignalParams->cosVal = ( REAL8 * )LALMalloc( ( pSFTandSignalParams->resTrig + 1 ) * sizeof( REAL8 ) );
585
586 for ( k = 0; k <= pSFTandSignalParams->resTrig; k++ ) {
587 pSFTandSignalParams->trigArg[k] = ( ( REAL8 )LAL_TWOPI ) * ( ( REAL8 )k ) / ( ( REAL8 )pSFTandSignalParams->resTrig );
588 pSFTandSignalParams->sinVal[k] = sin( pSFTandSignalParams->trigArg[k] );
589 pSFTandSignalParams->cosVal[k] = cos( pSFTandSignalParams->trigArg[k] );
590 }
591
592 pSFTandSignalParams->pSigParams = &params; /* as defined in Hough*/
593 pSFTandSignalParams->pSFTParams = &sftParams; /* as defined in Hough*/
594 pSFTandSignalParams->nSamples = ( INT4 )( 0.5 * timeBase ); /* nsample to get version 2 sfts */
595 pSFTandSignalParams->Dterms = uvar_Dterms;
596
597 }
598
599 /******************************************************************/
600 /* setting of parameters */
601 /******************************************************************/
602 injectPar.h0 = uvar_h0Min;
603 injectPar.fmin = uvar_f0;
604 injectPar.fSearchBand = uvar_fSearchBand;
605 injectPar.deltaF = deltaF;
606 injectPar.alpha = uvar_alpha; /* patch center if not full sky */
607 injectPar.delta = uvar_delta;
608 injectPar.patchSizeAlpha = uvar_patchSizeAlpha; /* patch size if not full sky */
609 injectPar.patchSizeDelta = uvar_patchSizeDelta;
610 injectPar.pixelFactor = uvar_pixelFactor;
611 injectPar.vTotC = VTOT;
612 injectPar.timeObs = tObs;
613 injectPar.spnFmax.data = NULL;
614 injectPar.spnFmax.length = msp; /*only 1 spin */
615 injectPar.spnFmax.data = ( REAL8 * )LALMalloc( msp * sizeof( REAL8 ) );
616 injectPar.spnFmax.data[0] = uvar_maxSpin;
617 injectPar.spnFmin.data = NULL;
618 injectPar.spnFmin.length = msp; /*only 1 spin */
619 injectPar.spnFmin.data = ( REAL8 * )LALMalloc( msp * sizeof( REAL8 ) );
620 injectPar.spnFmin.data[0] = uvar_minSpin;
621
622 pulsarInject.spindown.length = msp;
623 pulsarTemplate.spindown.length = msp;
624 pulsarInject.spindown.data = NULL;
625 pulsarTemplate.spindown.data = NULL;
626 pulsarInject.spindown.data = ( REAL8 * )LALMalloc( msp * sizeof( REAL8 ) );
627 pulsarTemplate.spindown.data = ( REAL8 * )LALMalloc( msp * sizeof( REAL8 ) );
628
629 sftParams.Tsft = timeBase;
630 sftParams.noiseSFTs = NULL;
631
632 params.orbit.asini = 0 /* isolated pulsar */;
633 /* params.transferFunction = NULL; */
634 params.ephemerides = edat;
635 params.startTimeGPS.gpsSeconds = firstTimeStamp.gpsSeconds; /* start time of output time series */
636 params.startTimeGPS.gpsNanoSeconds = firstTimeStamp.gpsNanoSeconds; /* start time of output time series */
637 params.duration = injectPar.timeObs; /* length of time series in seconds */
638 params.samplingRate = tSamplingRate;
639 params.fHeterodyne = fHeterodyne;
640 params.pulsar.refTime.gpsSeconds = firstTimeStamp.gpsSeconds;
641 params.pulsar.refTime.gpsNanoSeconds = firstTimeStamp.gpsNanoSeconds;
642 /* ****************************************************************/
643
644 /* WE SHOULD LOOP OVER MC SIGNAL INJECTION HERE
645 BEFORE THAT :
646 -for each different h0 value create a file containing the h0 value
647 LOOP over xxx Monte-Carlo signal Injections:
648 - Generate signal injections parameters and the corresponding template
649 parameters allowing some mismatch
650 - Compute the frequency path for the template parameters
651 -Generate the time series for injected signals and the
652 corresponding SFTs with no added noise (for all times).
653
654 LOOP over the different h0 values
655 -Add SFT with the signal normalized to the SFT original noise
656 -clean lines, compute weights, select peaks
657 -get the significance
658 */
659
660 /* ****************************************************************/
661
662 pg1.length = binsSFT; /*equal to binsSFT */
663 pg1.data = NULL;
664 pg1.data = ( UCHAR * )LALMalloc( binsSFT * sizeof( UCHAR ) );
665
666 /* ****************************************************************/
667 foft.length = mObsCoh;
668 foft.data = NULL;
669 foft.data = ( REAL8 * )LALMalloc( mObsCoh * sizeof( REAL8 ) );
670 {
671 UINT4 j;
672 for ( j = 0; j < nTemplates; ++j ) {
673 foftV[j].length = mObsCoh;
674 foftV[j].data = NULL;
675 foftV[j].data = ( REAL8 * )LALMalloc( mObsCoh * sizeof( REAL8 ) );
676 }
677 }
678
679 /* ****************************************************************/
680 /* HERE SHOULD START THE MONTE-CARLO */
681
682 for ( MCloopId = 0; MCloopId < uvar_nMCloop; ++MCloopId ) {
683
684 LAL_CALL( GenerateInjectParamsNoVeto( &status, &pulsarInject, &pulsarTemplate,
685 &closeTemplates, &injectPar ), &status );
686
687 /* writing the parameters into fpPar, following the format
688 MCloopId I.f0 H.f0 I.f1 H.f1 I.alpha H.alpha I.delta H.delta I.phi0 I.psi
689 (not cos iota) */
690
691 fprintf( fpPar, "%d %f %f %g %g %f %f %f %f %f %f %f\n",
692 MCloopId, pulsarInject.f0, pulsarTemplate.f0,
693 pulsarInject.spindown.data[0], pulsarTemplate.spindown.data[0],
694 pulsarInject.longitude, pulsarTemplate.longitude,
695 pulsarInject.latitude, pulsarTemplate.latitude,
696 pulsarInject.phi0, pulsarInject.psi, ( pulsarInject.aCross ) / injectPar.h0
697 );
698
699
700 /* ****************************************************************/
701 /* Computing the frequency path f(t) = f0(t)* (1+v/c.n) for all the different templates */
702
703 /* the geometrically nearest template */
704 LAL_CALL( ComputeFoft_NM( &status, &foft, &pulsarTemplate, &timeDiffV, &velV ), &status );
705
706 /* for all the 16 near templates */
707 {
708 UINT4 j, i, k, itemplate;
709
710 itemplate = 0;
711 for ( j = 0; j < 2; ++j ) {
712 pulsarTemplate.f0 = closeTemplates.f0[j];
713 for ( i = 0; i < 2; ++i ) {
714 pulsarTemplate.spindown.data[0] = closeTemplates.f1[i];
715 for ( k = 0; k < 4; ++k ) {
716 pulsarTemplate.latitude = closeTemplates.skytemp[k].delta;
717 pulsarTemplate.longitude = closeTemplates.skytemp[k].alpha;
718 LAL_CALL( ComputeFoft_NM( &status, &( foftV[itemplate] ),
719 &pulsarTemplate, &timeDiffV, &velV ), &status );
720 ++itemplate;
721 }
722 }
723 }
724 }
725
726 /* ****************************************************************/
727
728 /* params.pulsar.TRefSSB= ? ; */
729 params.pulsar.position.longitude = pulsarInject.longitude;
730 params.pulsar.position.latitude = pulsarInject.latitude ;
731 params.pulsar.position.system = COORDINATESYSTEM_EQUATORIAL;
732 params.pulsar.psi = pulsarInject.psi;
733 params.pulsar.aPlus = pulsarInject.aPlus;
734 params.pulsar.aCross = pulsarInject.aCross;
735 params.pulsar.phi0 = pulsarInject.phi0;
736 params.pulsar.f0 = pulsarInject.f0;
737 params.pulsar.spindown = &pulsarInject.spindown ;
738
739 {
740 UINT4 iIFO, numsft, iSFT, j;
741
742 if ( uvar_fast ) {
743
744 for ( iIFO = 0; iIFO < numifo; iIFO++ ) {
745 params.site = &( mdetStates->data[iIFO]->detector );
746 sftParams.timestamps = multiIniTimeV->data[iIFO];
747 numsft = mdetStates->data[iIFO]->length;
748
749 /* initialize data to zero */
750 for ( iSFT = 0; iSFT < numsft; iSFT++ ) {
751 for ( j = 0; j < binsSFT; j++ ) {
752 signalSFTs->data[iIFO]->data[iSFT].data->data[j] = 0.0;
753 }
754 }
755
756 LAL_CALL( LALComputeSkyAndZeroPsiAMResponse( &status,
757 &pSkyConstAndZeroPsiAMResponse[iIFO], pSFTandSignalParams ), &status );
758 LAL_CALL( LALFastGeneratePulsarSFTs( &status, &signalSFTs->data[iIFO],
759 &pSkyConstAndZeroPsiAMResponse[iIFO], pSFTandSignalParams ), &status );
760 }
761 } else {
762
763 for ( iIFO = 0; iIFO < numifo; iIFO++ ) {
764 params.site = &( mdetStates->data[iIFO]->detector );
765 sftParams.timestamps = multiIniTimeV->data[iIFO];
766
767 XLALDestroySFTVector( signalSFTs->data[iIFO] );
768 signalSFTs->data[iIFO] = NULL;
769
770 LAL_CALL( LALGeneratePulsarSignal( &status, &signalTseries, &params ), &status );
771 LAL_CALL( LALSignalToSFTs( &status, &signalSFTs->data[iIFO], signalTseries, &sftParams ), &status );
772
773 LALFree( signalTseries->data->data );
774 LALFree( signalTseries->data );
775 LALFree( signalTseries );
776 signalTseries = NULL;
777 }
778 }
779 }
780
781
782 /******************************************************************/
783 /* initialize all weights to unity */
784 /******************************************************************/
785
786 /* set up weights -- this should be done before normalizing the sfts */
787 weightsV.length = mObsCoh;
788 weightsV.data = ( REAL8 * )LALCalloc( mObsCoh, sizeof( REAL8 ) );
789
790 weightsNoise.length = mObsCoh;
791 weightsNoise.data = ( REAL8 * )LALCalloc( mObsCoh, sizeof( REAL8 ) );
792
793 /* initialize all weights to unity */
794 LAL_CALL( LALHOUGHInitializeWeights( &status, &weightsNoise ), &status );
795 LAL_CALL( LALHOUGHInitializeWeights( &status, &weightsV ), &status );
796
797 /******************************************************************/
798 /* setting the weights considering only the AM coefficients to be only
799 computed just where we need*/
800 /******************************************************************/
801 if ( uvar_weighAM ) {
802 SkyPosition skypos;
803 UINT4 iIFO, iSFT;
804 UINT4 k, numsft;
805 MultiAMCoeffs *multiAMcoef = NULL;
806
807 weightsAM.length = mObsCoh;
808 weightsAM.data = ( REAL8 * )LALCalloc( mObsCoh, sizeof( REAL8 ) );
809 skypos.system = COORDINATESYSTEM_EQUATORIAL;
810
811 skypos.longitude = pulsarInject.longitude;
812 skypos.latitude = pulsarInject.latitude;
813 XLAL_CHECK_MAIN( ( multiAMcoef = XLALComputeMultiAMCoeffs( mdetStates, NULL, skypos ) ) != NULL, XLAL_EFUNC );
814
815 /* loop over the weights and set them by the appropriate AM coefficients */
816 for ( k = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
817 numsft = mdetStates->data[iIFO]->length;
818 for ( iSFT = 0; iSFT < numsft; iSFT++, k++ ) {
819 REAL8 a, b;
820
821 a = multiAMcoef->data[iIFO]->a->data[iSFT];
822 b = multiAMcoef->data[iIFO]->b->data[iSFT];
823 weightsAM.data[k] = ( a * a + b * b );
824 } /* loop over SFTs */
825 } /* loop over IFOs */
826
827 XLALDestroyMultiAMCoeffs( multiAMcoef );
828
829 }
830
831
832 /* ****************************************************************/
833 /* HERE THE LOOP FOR DIFFERENT h0 VALUES */
834
835 fprintf( fpNc, " %d ", MCloopId );
836
837 for ( h0loop = 0; h0loop < uvar_nh0; ++h0loop ) {
838
839 UINT4 j, ind, itemplate;
840 UINT4 iIFO, numsft, iSFT;
841 COMPLEX8 *noiseSFT;
842 COMPLEX8 *signalSFT;
843 COMPLEX8 *sumSFT;
844
845 numberCount = 0.0;
846 for ( itemplate = 0; itemplate < nTemplates; ++itemplate ) {
847 numberCountV[itemplate] = 0.0;
848 }
849
850 h0scale = h0V.data[h0loop] / h0V.data[0]; /* different for different h0 values */
851
852 /* ****************************************************************/
853 /* adding signal+ noise SFT, TO BE CHECKED */
854
855 for ( iIFO = 0; iIFO < numifo; iIFO++ ) {
856 numsft = inputSFTs->data[iIFO]->length;
857 for ( iSFT = 0; iSFT < numsft; iSFT++ ) {
858
859 noiseSFT = inputSFTs->data[iIFO]->data[iSFT].data->data;
860 signalSFT = signalSFTs->data[iIFO]->data[iSFT].data->data;
861 sumSFT = sumSFTs->data[iIFO]->data[iSFT].data->data;
862
863 for ( j = 0; j < binsSFT; j++ ) {
864 *( sumSFT ) = crectf( crealf( *noiseSFT ) + h0scale * crealf( *signalSFT ), cimagf( *noiseSFT ) + h0scale * cimagf( *signalSFT ) );
865 ++noiseSFT;
866 ++signalSFT;
867 ++sumSFT;
868 }
869 }
870 }
871
872 /* ****************************************************************/
873 /* clean sfts if required */
874 if ( XLALUserVarWasSet( &uvar_linefiles ) ) {
875 RandomParams *randPar = NULL;
876 FILE *fpRand = NULL;
877 INT4 seed, ranCount;
878
879 if ( ( fpRand = fopen( "/dev/urandom", "r" ) ) == NULL ) {
880 fprintf( stderr, "Error in opening /dev/urandom" );
881 exit( 1 );
882 }
883
884 if ( ( ranCount = fread( &seed, sizeof( seed ), 1, fpRand ) ) != 1 ) {
885 fprintf( stderr, "Error in getting random seed" );
886 exit( 1 );
887 }
888
889 LAL_CALL( LALCreateRandomParams( &status, &randPar, seed ), &status );
890
891 LAL_CALL( LALRemoveKnownLinesInMultiSFTVector( &status, sumSFTs, uvar_maxBinsClean, uvar_blocksRngMed, uvar_linefiles, randPar ), &status );
892
893 LAL_CALL( LALDestroyRandomParams( &status, &randPar ), &status );
894 fclose( fpRand );
895 } /* end cleaning */
896
897
898 /* ****************************************************************/
899 /* normalize sfts compute weights */
900 {
901 MultiNoiseWeights *multweight = NULL;
902 MultiPSDVector *multPSD = NULL;
903 REAL8 sumWeightSquare;
904
905
906 /* initialize all weights to unity each time*/
907 LAL_CALL( LALHOUGHInitializeWeights( &status, &weightsNoise ), &status );
908 LAL_CALL( LALHOUGHInitializeWeights( &status, &weightsV ), &status );
909
910
911 /* normalize sfts */
912 XLAL_CHECK_MAIN( ( multPSD = XLALNormalizeMultiSFTVect( sumSFTs, uvar_blocksRngMed, NULL ) ) != NULL, XLAL_EFUNC );
913
914 /* compute multi noise weights */
915 if ( uvar_weighNoise ) {
916 XLAL_CHECK_MAIN( ( multweight = XLALComputeMultiNoiseWeights( multPSD, uvar_blocksRngMed, 0 ) ) != NULL, XLAL_EFUNC );
917 }
918
919 /* we are now done with the psd */
920 XLALDestroyMultiPSDVector( multPSD );
921
922 /* copy weights */
923 if ( uvar_weighNoise ) {
924 for ( j = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
925 numsft = mdetStates->data[iIFO]->length;
926 for ( iSFT = 0; iSFT < numsft; iSFT++, j++ ) {
927 weightsNoise.data[j] = multweight->data[iIFO]->data[iSFT];
928 } /* loop over SFTs */
929 } /* loop over IFOs */
930
931 XLALDestroyMultiNoiseWeights( multweight );
932 memcpy( weightsV.data, weightsNoise.data, mObsCoh * sizeof( REAL8 ) );
933 }
934
935 if ( uvar_weighAM && weightsAM.data ) {
936 for ( j = 0; j < mObsCoh; j++ ) {
937 weightsV.data[j] = weightsV.data[j] * weightsAM.data[j];
938 }
939 }
940
941 LAL_CALL( LALHOUGHNormalizeWeights( &status, &weightsV ), &status );
942
943 /* calculate the sum of the weights squared */
944 sumWeightSquare = 0.0;
945 for ( j = 0; j < mObsCoh; j++ ) {
946 sumWeightSquare += weightsV.data[j] * weightsV.data[j];
947 }
948
949 /* standard deviation for noise only */
950 sigmaN = sqrt( sumWeightSquare * alphaPeak * ( 1.0 - alphaPeak ) );
951 }
952
953 /* ****************************************************************/
954 /* loop over SFT, generate peakgram and get number count */
955 {
956 SFTtype *sft;
957
958 for ( j = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
959 numsft = mdetStates->data[iIFO]->length;
960
961 for ( iSFT = 0; iSFT < numsft; iSFT++, j++ ) {
962
963 sft = sumSFTs->data[iIFO]->data + iSFT;
964 LAL_CALL( SFTtoUCHARPeakGram( &status, &pg1, sft, uvar_peakThreshold ), &status );
965 ind = floor( foft.data[j] * timeBase - sftFminBin + 0.5 );
966 numberCount += pg1.data[ind] * weightsV.data[j]; /* adds 0 or 1 to the counter*/
967
968 for ( itemplate = 0; itemplate < nTemplates; ++itemplate ) {
969 ind = floor( foftV[itemplate].data[j] * timeBase - sftFminBin + 0.5 );
970 numberCountV[itemplate] += pg1.data[ind] * weightsV.data[j];
971 }
972 } /* loop over SFTs */
973 } /* loop over IFOs */
974
975 }
976
977 /* ****************************************************************/
978 /*check the max number count */
979 maxNumberCount = numberCount;
980 for ( itemplate = 0; itemplate < nTemplates; ++itemplate ) {
981 if ( numberCountV[itemplate] > maxNumberCount ) {
982 maxNumberCount = numberCountV[itemplate];
983 }
984 }
985
986 /******************************************************************/
987 /* printing the significance in the proper file */
988 /******************************************************************/
989 significance = ( maxNumberCount - meanN ) / sigmaN;
990 fprintf( fpNc, " %f ", significance );
991
992 } /* closing loop for different h0 values */
993 fprintf( fpNc, " \n" );
994
995
996 LALFree( weightsV.data );
997 LALFree( weightsNoise.data );
998 if ( uvar_weighAM && weightsAM.data ) {
999 LALFree( weightsAM.data );
1000 }
1001
1002
1003 } /* Closing MC loop */
1004
1005 /******************************************************************/
1006 /* Closing files */
1007 /******************************************************************/
1008 fclose( fpPar );
1009 fclose( fpNc );
1010
1011
1012 /******************************************************************/
1013 /* Free memory and exit */
1014 /******************************************************************/
1015
1016 /* LALFree(fp); */
1017 LALFree( pg1.data );
1018
1019 LALFree( velV.data );
1020 LALFree( timeDiffV.data );
1021
1022 LALFree( pSFTandSignalParams->cosVal );
1023 LALFree( pSFTandSignalParams->sinVal );
1024 LALFree( pSFTandSignalParams->trigArg );
1025 LALFree( pSFTandSignalParams );
1026
1027
1028 {
1029 UINT4 iIFO;
1030
1031 for ( iIFO = 0; iIFO < numifo; iIFO++ ) {
1032 XLALDestroyTimestampVector( multiIniTimeV->data[iIFO] );
1033 LALFree( pSkyConstAndZeroPsiAMResponse[iIFO].skyConst );
1034 LALFree( pSkyConstAndZeroPsiAMResponse[iIFO].fPlusZeroPsi );
1035 LALFree( pSkyConstAndZeroPsiAMResponse[iIFO].fCrossZeroPsi );
1036 }
1037 }
1038
1039 LALFree( multiIniTimeV->data );
1040 LALFree( multiIniTimeV );
1041 LALFree( pSkyConstAndZeroPsiAMResponse );
1042
1043 XLALDestroyMultiDetectorStateSeries( mdetStates );
1044
1045 LALFree( foft.data );
1046 LALFree( h0V.data );
1047
1048 {
1049 UINT4 j;
1050 for ( j = 0; j < nTemplates; ++j ) {
1051 LALFree( foftV[j].data );
1052 }
1053 }
1054
1055
1056 LALFree( injectPar.spnFmax.data );
1057 LALFree( pulsarInject.spindown.data );
1058 LALFree( pulsarTemplate.spindown.data );
1059
1060 XLALDestroyEphemerisData( edat );
1061
1062 XLALDestroyMultiSFTVector( inputSFTs );
1063 XLALDestroyMultiSFTVector( sumSFTs );
1064 XLALDestroyMultiSFTVector( signalSFTs );
1065
1066 XLALDestroyUserVars();
1067
1068 LALCheckMemoryLeaks();
1069
1070 if ( lalDebugLevel ) {
1071 REPORTSTATUS( &status );
1072 }
1073
1074 return status.statusCode;
1075}
1076
1077
1078
1079
1080/***************************************************************************/
1081void GenerateInjectParams( LALStatus *status,
1082 PulsarData *injectPulsar,
1083 HoughTemplate *templatePulsar,
1084 HoughNearTemplates *closeTemplates,
1085 HoughInjectParams *params,
1086 LineNoiseInfo *lines )
1087{
1088
1089 INT4 seed = 0; /* seed generated using current time */
1090 REAL4 randval;
1091 RandomParams *randPar = NULL;
1092 FILE *fpRandom;
1093 INT4 count;
1094
1095 REAL4 cosiota, h0;
1096 REAL8 f0, deltaF, deltaX;
1097 REAL8 latitude, longitude; /* of the source in radians */
1098 INT8 f0bin;
1099 UINT4 msp;
1100
1101 /* --------------------------------------------- */
1102 INITSTATUS( status );
1103 ATTATCHSTATUSPTR( status );
1104
1105 /* Make sure the arguments are not NULL: */
1106 ASSERT( injectPulsar, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1107 ASSERT( templatePulsar, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1108 ASSERT( params, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1109 ASSERT( lines, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1110
1111 /* ++++++++++++++++++from makefakedata
1112 * Modified so as to not create random number parameters with seed
1113 * drawn from clock. Seconds don't change fast enough and sft's
1114 * look alike. We open /dev/urandom and read a 4 byte integer from
1115 * it and use that as our seed. Note: /dev/random is slow after the
1116 * first, few accesses.
1117 */
1118
1119 fpRandom = fopen( "/dev/urandom", "r" );
1120 ASSERT( fpRandom, status, DRIVEHOUGHCOLOR_EFILE, DRIVEHOUGHCOLOR_MSGEFILE );
1121
1122 count = fread( &seed, sizeof( INT4 ), 1, fpRandom );
1123 if ( count == 0 ) {
1124 ABORT( status, DRIVEHOUGHCOLOR_EARG, DRIVEHOUGHCOLOR_MSGEARG );
1125 }
1126
1127 fclose( fpRandom );
1128
1129 TRY( LALCreateRandomParams( status->statusPtr, &randPar, seed ), status );
1130
1131 /*
1132 * to create a single random deviate distributed uniforly between zero and unity
1133 * TRY( LALUniformDeviate(status->statusPtr, &randval, randPar), status);
1134 */
1135
1136
1137 /* get random value phi0 [0, 2 pi] */
1138 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1139 injectPulsar->phi0 = randval * LAL_TWOPI;
1140
1141 /* get random value cos iota [-1,1] */
1142 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1143 cosiota = 2.0 * randval - 1.0;
1144
1145 h0 = params->h0;
1146 injectPulsar->aCross = h0 * cosiota;
1147 injectPulsar->aPlus = 0.5 * h0 * ( 1.0 + cosiota * cosiota );
1148
1149 /* get random value psi [0, 2 pi] */
1150 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1151 injectPulsar->psi = randval * LAL_TWOPI;
1152
1153 /* getting random number for the frequency (and mismatch)*/
1154 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1155 f0 = params->fmin + ( params->fSearchBand ) * randval;
1156
1157 /* veto the frequency if it is affected by a line */
1158 {
1159 INT4 veto = 1;
1160 while ( veto > 0 ) {
1161
1162 TRY( LALCheckLines( status->statusPtr, &veto, lines, f0 ), status );
1163 if ( veto > 0 ) {
1164 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1165 f0 = params->fmin + ( params->fSearchBand ) * randval;
1166 }
1167 } /* end of while loop */
1168 }
1169
1170 injectPulsar->f0 = f0;
1171 deltaF = params->deltaF;
1172 f0bin = floor( f0 / deltaF + 0.5 );
1173 templatePulsar->f0 = f0bin * deltaF;
1174 closeTemplates->f0[0] = floor( f0 / deltaF ) * deltaF;
1175 closeTemplates->f0[1] = ceil( f0 / deltaF ) * deltaF;
1176
1177 /* sky location, depending if full sky or small patch is analyzed */
1178 /*
1179 * deltaX = deltaF/(params->vTotC * params->pixelFactor *
1180 * (params->fmin + params->fSearchBand) );
1181 */
1182 deltaX = deltaF / ( params->vTotC * params->pixelFactor * f0 );
1183
1184
1185 if ( params->fullSky ) { /*full sky*/
1186 REAL8 kkcos;
1187
1188 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1189 longitude = randval * LAL_TWOPI;
1190 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1191 kkcos = 2.0 * randval - 1.0;
1192 latitude = acos( kkcos ) - LAL_PI_2;
1193 } else { /*small patch */
1194 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1195 longitude = params->alpha + ( params->patchSizeAlpha ) * ( randval - 0.5 );
1196 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1197 latitude = params->delta + ( params->patchSizeDelta ) * ( randval - 0.5 );
1198 }
1199
1200 injectPulsar->longitude = longitude;
1201 injectPulsar->latitude = latitude;
1202
1203 {
1204 REAL8UnitPolarCoor template, par;
1205 REAL8UnitPolarCoor templRotated;
1206 REAL8Cart2Coor templProjected;
1207 REAL8 dX1[2], dX2[2];
1208 INT4 ii, jj, kk;
1209
1210 par.alpha = injectPulsar->longitude;
1211 par.delta = injectPulsar->latitude;
1212
1213 /* mismatch with the template in stereographic plane */
1214 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1215 templProjected.x = dX1[0] = deltaX * ( randval - 0.5 );
1216 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1217 templProjected.y = dX2[0] = deltaX * ( randval - 0.5 );
1218
1219 if ( dX1[0] < 0.0 ) {
1220 dX1[1] = dX1[0] + deltaX;
1221 } else {
1222 dX1[1] = dX1[0] - deltaX;
1223 }
1224
1225 if ( dX2[0] < 0.0 ) {
1226 dX2[1] = dX2[0] + deltaX;
1227 } else {
1228 dX2[1] = dX2[0] - deltaX;
1229 }
1230
1231 /* invert the stereographic projection for a point on the projected plane */
1232 TRY( LALStereoInvProjectCart( status->statusPtr,
1233 &templRotated, &templProjected ), status );
1234 /* inverse rotate the mismatch from the south pole to desired location */
1235 TRY( LALInvRotatePolarU( status->statusPtr, &template, &templRotated, &par ), status );
1236 templatePulsar->longitude = template.alpha;
1237 templatePulsar->latitude = template.delta;
1238
1239 kk = 0;
1240 for ( ii = 0; ii < 2; ii++ ) {
1241 for ( jj = 0; jj < 2; jj++ ) {
1242 templProjected.x = dX1[ii];
1243 templProjected.y = dX2[jj];
1244 TRY( LALStereoInvProjectCart( status->statusPtr,
1245 &templRotated, &templProjected ), status );
1246 TRY( LALInvRotatePolarU( status->statusPtr, &( closeTemplates->skytemp[kk] ), &templRotated,
1247 &par ), status );
1248 ++kk;
1249 }
1250 }
1251
1252 }
1253
1254 /* now the spindown if any */
1255 msp = params->spnFmax.length ;
1256 closeTemplates->f1[0] = 0.0;
1257 closeTemplates->f1[1] = 0.0;
1258
1259 ASSERT( templatePulsar->spindown.length == msp, status, DRIVEHOUGHCOLOR_EBAD,
1260 DRIVEHOUGHCOLOR_MSGEBAD );
1261 ASSERT( injectPulsar->spindown.length == msp, status, DRIVEHOUGHCOLOR_EBAD,
1262 DRIVEHOUGHCOLOR_MSGEBAD );
1263
1264 if ( msp ) { /*if there are spin-down values */
1265 REAL8 deltaFk, spink;
1266 REAL8 timeObsInv;
1267 UINT4 i;
1268 ASSERT( injectPulsar->spindown.data, status, DRIVEHOUGHCOLOR_ENULL,
1269 DRIVEHOUGHCOLOR_MSGENULL );
1270 ASSERT( templatePulsar->spindown.data, status, DRIVEHOUGHCOLOR_ENULL,
1271 DRIVEHOUGHCOLOR_MSGENULL );
1272 ASSERT( params->spnFmax.data, status, DRIVEHOUGHCOLOR_ENULL,
1273 DRIVEHOUGHCOLOR_MSGENULL );
1274
1275 /* delta f_k = k! deltaF/ [T_Obs}^k spd grid resolution*/
1276 timeObsInv = 1.0 / params->timeObs;
1277 deltaFk = deltaF * timeObsInv;
1278
1279 /* first spin-down parameter, (only spin-down) */
1280 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1281 spink = randval * ( params->spnFmax.data[0] - params->spnFmin.data[0] ) + params->spnFmin.data[0];
1282
1283 injectPulsar->spindown.data[0] = spink;
1284 templatePulsar->spindown.data[0] = floor( spink / deltaFk + 0.5 ) * deltaFk;
1285
1286 closeTemplates->f1[0] = floor( spink / deltaFk ) * deltaFk;
1287 closeTemplates->f1[1] = ceil( spink / deltaFk ) * deltaFk;
1288
1289 /* the rest of the spin orders */
1290 for ( i = 1; i < msp; ++i ) {
1291 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1292 spink = params->spnFmax.data[i] * ( 2.0 * randval - 1.0 );
1293 injectPulsar->spindown.data[i] = spink;
1294 deltaFk = deltaFk * timeObsInv * ( i + 1.0 );
1295 templatePulsar->spindown.data[i] = floor( spink / deltaFk + 0.5 ) * deltaFk;
1296 }
1297 }
1298 /* free memory */
1299 TRY( LALDestroyRandomParams( status->statusPtr, &randPar ), status );
1300
1301 DETATCHSTATUSPTR( status );
1302 /* normal exit */
1303 RETURN( status );
1304}
1305
1306
1307/***************************************************************************/
1308void GenerateInjectParamsNoVeto( LALStatus *status,
1309 PulsarData *injectPulsar,
1310 HoughTemplate *templatePulsar,
1311 HoughNearTemplates *closeTemplates,
1312 HoughInjectParams *params )
1313{
1314
1315 INT4 seed = 0; /* seed generated using current time */
1316 REAL4 randval;
1317 RandomParams *randPar = NULL;
1318 FILE *fpRandom;
1319 INT4 count;
1320
1321 REAL4 cosiota, h0;
1322 REAL8 f0, deltaF, deltaX;
1323 REAL8 latitude, longitude; /* of the source in radians */
1324 INT8 f0bin;
1325 UINT4 msp;
1326
1327 /* --------------------------------------------- */
1328 INITSTATUS( status );
1329 ATTATCHSTATUSPTR( status );
1330
1331 /* Make sure the arguments are not NULL: */
1332 ASSERT( injectPulsar, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1333 ASSERT( templatePulsar, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1334 ASSERT( params, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1335
1336 /* ++++++++++++++++++from makefakedata
1337 * Modified so as to not create random number parameters with seed
1338 * drawn from clock. Seconds don't change fast enough and sft's
1339 * look alike. We open /dev/urandom and read a 4 byte integer from
1340 * it and use that as our seed. Note: /dev/random is slow after the
1341 * first, few accesses.
1342 */
1343
1344 fpRandom = fopen( "/dev/urandom", "r" );
1345 ASSERT( fpRandom, status, DRIVEHOUGHCOLOR_EFILE, DRIVEHOUGHCOLOR_MSGEFILE );
1346
1347 count = fread( &seed, sizeof( INT4 ), 1, fpRandom );
1348 if ( count == 0 ) {
1349 ABORT( status, DRIVEHOUGHCOLOR_EARG, DRIVEHOUGHCOLOR_MSGEARG );
1350 }
1351
1352 fclose( fpRandom );
1353
1354 TRY( LALCreateRandomParams( status->statusPtr, &randPar, seed ), status );
1355
1356 /*
1357 * to create a single random deviate distributed uniforly between zero and unity
1358 * TRY( LALUniformDeviate(status->statusPtr, &randval, randPar), status);
1359 */
1360
1361
1362 /* get random value phi0 [0, 2 pi] */
1363 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1364 injectPulsar->phi0 = randval * LAL_TWOPI;
1365
1366 /* get random value cos iota [-1,1] */
1367 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1368 cosiota = 2.0 * randval - 1.0;
1369
1370 h0 = params->h0;
1371 injectPulsar->aCross = h0 * cosiota;
1372 injectPulsar->aPlus = 0.5 * h0 * ( 1.0 + cosiota * cosiota );
1373
1374 /* get random value psi [0, 2 pi] */
1375 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1376 injectPulsar->psi = randval * LAL_TWOPI;
1377
1378 /* getting random number for the frequency (and mismatch)*/
1379 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1380 f0 = params->fmin + ( params->fSearchBand ) * randval;
1381
1382 injectPulsar->f0 = f0;
1383 deltaF = params->deltaF;
1384 f0bin = floor( f0 / deltaF + 0.5 );
1385 templatePulsar->f0 = f0bin * deltaF;
1386 closeTemplates->f0[0] = floor( f0 / deltaF ) * deltaF;
1387 closeTemplates->f0[1] = ceil( f0 / deltaF ) * deltaF;
1388
1389 /* sky location, depending if full sky or small patch is analyzed */
1390 deltaX = deltaF / ( params->vTotC * params->pixelFactor *
1391 ( params->fmin + params->fSearchBand ) );
1392
1393
1394 if ( params->fullSky ) { /*full sky*/
1395 REAL8 kkcos;
1396
1397 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1398 longitude = randval * LAL_TWOPI;
1399 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1400 kkcos = 2.0 * randval - 1.0;
1401 latitude = acos( kkcos ) - LAL_PI_2;
1402 } else { /*small patch */
1403 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1404 longitude = params->alpha + ( params->patchSizeAlpha ) * ( randval - 0.5 );
1405 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1406 latitude = params->delta + ( params->patchSizeDelta ) * ( randval - 0.5 );
1407 }
1408
1409 injectPulsar->longitude = longitude;
1410 injectPulsar->latitude = latitude;
1411
1412 {
1413 REAL8UnitPolarCoor template, par;
1414 REAL8UnitPolarCoor templRotated;
1415 REAL8Cart2Coor templProjected;
1416 REAL8 dX1[2], dX2[2];
1417 INT4 ii, jj, kk;
1418
1419 par.alpha = injectPulsar->longitude;
1420 par.delta = injectPulsar->latitude;
1421
1422 /* mismatch with the template in stereographic plane */
1423 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1424 templProjected.x = dX1[0] = deltaX * ( randval - 0.5 );
1425 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1426 templProjected.y = dX2[0] = deltaX * ( randval - 0.5 );
1427
1428 if ( dX1[0] < 0.0 ) {
1429 dX1[1] = dX1[0] + deltaX;
1430 } else {
1431 dX1[1] = dX1[0] - deltaX;
1432 }
1433
1434 if ( dX2[0] < 0.0 ) {
1435 dX2[1] = dX2[0] + deltaX;
1436 } else {
1437 dX2[1] = dX2[0] - deltaX;
1438 }
1439
1440 /* invert the stereographic projection for a point on the projected plane */
1441 TRY( LALStereoInvProjectCart( status->statusPtr,
1442 &templRotated, &templProjected ), status );
1443 /* inverse rotate the mismatch from the south pole to desired location */
1444 TRY( LALInvRotatePolarU( status->statusPtr, &template, &templRotated, &par ), status );
1445 templatePulsar->longitude = template.alpha;
1446 templatePulsar->latitude = template.delta;
1447
1448 kk = 0;
1449 for ( ii = 0; ii < 2; ii++ ) {
1450 for ( jj = 0; jj < 2; jj++ ) {
1451 templProjected.x = dX1[ii];
1452 templProjected.y = dX2[jj];
1453 TRY( LALStereoInvProjectCart( status->statusPtr,
1454 &templRotated, &templProjected ), status );
1455 TRY( LALInvRotatePolarU( status->statusPtr, &( closeTemplates->skytemp[kk] ), &templRotated,
1456 &par ), status );
1457 ++kk;
1458 }
1459 }
1460
1461 }
1462
1463 /* now the spindown if any */
1464 msp = params->spnFmax.length ;
1465 closeTemplates->f1[0] = 0.0;
1466 closeTemplates->f1[1] = 0.0;
1467
1468 ASSERT( templatePulsar->spindown.length == msp, status, DRIVEHOUGHCOLOR_EBAD,
1469 DRIVEHOUGHCOLOR_MSGEBAD );
1470 ASSERT( injectPulsar->spindown.length == msp, status, DRIVEHOUGHCOLOR_EBAD,
1471 DRIVEHOUGHCOLOR_MSGEBAD );
1472
1473 if ( msp ) { /*if there are spin-down values */
1474 REAL8 deltaFk, spink;
1475 REAL8 timeObsInv;
1476 UINT4 i;
1477 ASSERT( injectPulsar->spindown.data, status, DRIVEHOUGHCOLOR_ENULL,
1478 DRIVEHOUGHCOLOR_MSGENULL );
1479 ASSERT( templatePulsar->spindown.data, status, DRIVEHOUGHCOLOR_ENULL,
1480 DRIVEHOUGHCOLOR_MSGENULL );
1481 ASSERT( params->spnFmax.data, status, DRIVEHOUGHCOLOR_ENULL,
1482 DRIVEHOUGHCOLOR_MSGENULL );
1483
1484 /* delta f_k = k! deltaF/ [T_Obs}^k spd grid resolution*/
1485 timeObsInv = 1.0 / params->timeObs;
1486 deltaFk = deltaF * timeObsInv;
1487
1488 /* first spin-down parameter, (only spin-down) */
1489 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1490 spink = randval * ( params->spnFmax.data[0] - params->spnFmin.data[0] ) + params->spnFmin.data[0];
1491
1492 injectPulsar->spindown.data[0] = spink;
1493 templatePulsar->spindown.data[0] = floor( spink / deltaFk + 0.5 ) * deltaFk;
1494
1495 closeTemplates->f1[0] = floor( spink / deltaFk ) * deltaFk;
1496 closeTemplates->f1[1] = ceil( spink / deltaFk ) * deltaFk;
1497
1498 /* the rest of the spin orders */
1499 for ( i = 1; i < msp; ++i ) {
1500 TRY( LALUniformDeviate( status->statusPtr, &randval, randPar ), status );
1501 spink = params->spnFmax.data[i] * ( 2.0 * randval - 1.0 );
1502 injectPulsar->spindown.data[i] = spink;
1503 deltaFk = deltaFk * timeObsInv * ( i + 1.0 );
1504 templatePulsar->spindown.data[i] = floor( spink / deltaFk + 0.5 ) * deltaFk;
1505 }
1506 }
1507 /* free memory */
1508 TRY( LALDestroyRandomParams( status->statusPtr, &randPar ), status );
1509
1510 DETATCHSTATUSPTR( status );
1511 /* normal exit */
1512 RETURN( status );
1513}
1514
1515
1516/* ****************************************************************/
1517/* Computing the frequency path f(t) = f0(t)* (1+v/c.n) */
1518/* without mismatch */
1519/* ****************************************************************/
1520/******************************************************************/
1521void ComputeFoft_NM( LALStatus *status,
1522 REAL8Vector *foft,
1523 HoughTemplate *pulsarTemplate,
1524 REAL8Vector *timeDiffV,
1525 REAL8Cart3CoorVector *velV )
1526{
1527
1528 INT4 mObsCoh;
1529 REAL8 f0new, vcProdn, timeDiffN;
1530 REAL8 sourceDelta, sourceAlpha, cosDelta;
1531 INT4 j, i, nspin, factorialN;
1532 REAL8Cart3Coor sourceLocation;
1533
1534 /* --------------------------------------------- */
1535 INITSTATUS( status );
1536 ATTATCHSTATUSPTR( status );
1537
1538 /* Make sure the arguments are not NULL: */
1539 ASSERT( foft, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1540 ASSERT( pulsarTemplate, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1541 ASSERT( timeDiffV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1542 ASSERT( velV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1543
1544 ASSERT( foft->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1545 ASSERT( timeDiffV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1546 ASSERT( velV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1547
1548 sourceDelta = pulsarTemplate->latitude;
1549 sourceAlpha = pulsarTemplate->longitude;
1550 cosDelta = cos( sourceDelta );
1551
1552 sourceLocation.x = cosDelta * cos( sourceAlpha );
1553 sourceLocation.y = cosDelta * sin( sourceAlpha );
1554 sourceLocation.z = sin( sourceDelta );
1555
1556 mObsCoh = foft->length;
1557 nspin = pulsarTemplate->spindown.length;
1558
1559 for ( j = 0; j < mObsCoh; ++j ) { /* loop for all different time stamps */
1560 vcProdn = velV->data[j].x * sourceLocation.x
1561 + velV->data[j].y * sourceLocation.y
1562 + velV->data[j].z * sourceLocation.z;
1563 f0new = pulsarTemplate->f0;
1564 factorialN = 1;
1565 timeDiffN = timeDiffV->data[j];
1566
1567 for ( i = 0; i < nspin; ++i ) { /* loop for spin-down values */
1568 factorialN *= ( i + 1 );
1569 f0new += pulsarTemplate->spindown.data[i] * timeDiffN / factorialN;
1570 timeDiffN *= timeDiffN;
1571 }
1572 foft->data[j] = f0new * ( 1.0 + vcProdn );
1573 }
1574
1575 DETATCHSTATUSPTR( status );
1576 /* normal exit */
1577 RETURN( status );
1578}
1579
1580/******************************************************************/
1581
1582
1583/* ****************************************************************/
1584/* PrintLogFile2 (like in the Driver, but this one doesn't */
1585/* copy the contents of skypatch file) */
1586/* ****************************************************************/
1587/******************************************************************/
1588
1589
1590void PrintLogFile2( LALStatus *status,
1591 CHAR *dir,
1592 CHAR *basename,
1593 LALStringVector *linefiles,
1594 CHAR *executable )
1595{
1596 CHAR *fnameLog = NULL;
1597 FILE *fpLog = NULL;
1598 CHAR *logstr = NULL;
1599 UINT4 k;
1600
1601 INITSTATUS( status );
1602 ATTATCHSTATUSPTR( status );
1603
1604 /* open log file for writing */
1605 fnameLog = ( CHAR * )LALCalloc( MAXFILENAMELENGTH, sizeof( CHAR ) );
1606 strcpy( fnameLog, dir );
1607 strcat( fnameLog, "/logfiles/" );
1608 /* now create directory fdirOut/logfiles using mkdir */
1609 errno = 0;
1610 {
1611 /* check whether file can be created or if it exists already
1612 if not then exit */
1613 INT4 mkdir_result;
1614 mkdir_result = mkdir( fnameLog, S_IRWXU | S_IRWXG | S_IRWXO );
1615 if ( ( mkdir_result == -1 ) && ( errno != EEXIST ) ) {
1616 fprintf( stderr, "unable to create logfiles directory %s\n", fnameLog );
1617 LALFree( fnameLog );
1618 exit( 1 ); /* stop the program */
1619 }
1620 }
1621
1622 /* create the logfilename in the logdirectory */
1623 strcat( fnameLog, basename );
1624 strcat( fnameLog, ".log" );
1625 /* open the log file for writing */
1626 if ( ( fpLog = fopen( fnameLog, "w" ) ) == NULL ) {
1627 fprintf( stderr, "Unable to open file %s for writing\n", fnameLog );
1628 LALFree( fnameLog );
1629 exit( 1 );
1630 }
1631
1632 /* get the log string */
1633 XLAL_CHECK_LAL( status, ( logstr = XLALUserVarGetLog( UVAR_LOGFMT_CFGFILE ) ) != NULL, XLAL_EFUNC );
1634
1635 fprintf( fpLog, "## LOG FILE FOR MC Inject Hough\n\n" );
1636 fprintf( fpLog, "# User Input:\n" );
1637 fprintf( fpLog, "#-------------------------------------------\n" );
1638 fprintf( fpLog, "%s", logstr );
1639 LALFree( logstr );
1640
1641 /* copy contents of linefile if necessary */
1642 if ( linefiles ) {
1643
1644 for ( k = 0; k < linefiles->length; k++ ) {
1645
1646 if ( ( fpLog = fopen( fnameLog, "a" ) ) != NULL ) {
1647 CHAR command[1024] = "";
1648 fprintf( fpLog, "\n\n# Contents of linefile %s :\n", linefiles->data[k] );
1649 fprintf( fpLog, "# -----------------------------------------\n" );
1650 fclose( fpLog );
1651 sprintf( command, "cat %s >> %s", linefiles->data[k], fnameLog );
1652 if ( system( command ) ) {
1653 fprintf( stderr, "\nsystem('%s') returned non-zero status!\n\n", command );
1654 }
1655 }
1656 }
1657 }
1658
1659 /* append an ident-string defining the exact CVS-version of the code used */
1660 if ( ( fpLog = fopen( fnameLog, "a" ) ) != NULL ) {
1661 CHAR command[1024] = "";
1662 fprintf( fpLog, "\n\n# CVS-versions of executable:\n" );
1663 fprintf( fpLog, "# -----------------------------------------\n" );
1664 fclose( fpLog );
1665
1666 sprintf( command, "ident %s | sort -u >> %s", executable, fnameLog );
1667 /* we don't check this. If it fails, we assume that */
1668 /* one of the system-commands was not available, and */
1669 /* therefore the CVS-versions will not be logged */
1670 if ( system( command ) ) {
1671 fprintf( stderr, "\nsystem('%s') returned non-zero status!\n\n", command );
1672 }
1673 }
1674
1675 LALFree( fnameLog );
1676
1677 DETATCHSTATUSPTR( status );
1678 /* normal exit */
1679 RETURN( status );
1680}
DRIVEHOUGHCOLOR_MSGEBAD
#define DRIVEHOUGHCOLOR_MSGEBAD
Definition: DriveHoughColor.h:102
DRIVEHOUGHCOLOR_MSGENULL
#define DRIVEHOUGHCOLOR_MSGENULL
Definition: DriveHoughColor.h:105
DRIVEHOUGHCOLOR_MSGEARG
#define DRIVEHOUGHCOLOR_MSGEARG
Definition: DriveHoughColor.h:101
DRIVEHOUGHCOLOR_EFILE
#define DRIVEHOUGHCOLOR_EFILE
Definition: DriveHoughColor.h:94
DRIVEHOUGHCOLOR_ENULL
#define DRIVEHOUGHCOLOR_ENULL
Definition: DriveHoughColor.h:96
DRIVEHOUGHCOLOR_EBAD
#define DRIVEHOUGHCOLOR_EBAD
Definition: DriveHoughColor.h:93
DRIVEHOUGHCOLOR_MSGEFILE
#define DRIVEHOUGHCOLOR_MSGEFILE
Definition: DriveHoughColor.h:103
DRIVEHOUGHCOLOR_EARG
#define DRIVEHOUGHCOLOR_EARG
Definition: DriveHoughColor.h:92
REPORTSTATUS
void REPORTSTATUS(LALStatus *status)
lal_errhandler
lal_errhandler_t lal_errhandler
LAL_ERR_EXIT
int LAL_ERR_EXIT(LALStatus *stat, const char *func, const char *file, const int line, volatile const char *id)
LAL_CALL
#define LAL_CALL(function, statusptr)
j
int j
k
int k
LALCheckMemoryLeaks
void LALCheckMemoryLeaks(void)
LALCalloc
#define LALCalloc(m, n)
LALMalloc
#define LALMalloc(n)
LALFree
#define LALFree(p)
lalPulsarVCSInfoList
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
ABORT
#define ABORT(statusptr, code, mesg)
XLAL_CHECK_LAL
#define XLAL_CHECK_LAL(sp, assertion,...)
TRY
#define TRY(func, statusptr)
ATTATCHSTATUSPTR
#define ATTATCHSTATUSPTR(statusptr)
ASSERT
#define ASSERT(assertion, statusptr, code, mesg)
DETATCHSTATUSPTR
#define DETATCHSTATUSPTR(statusptr)
INITSTATUS
#define INITSTATUS(statusptr)
RETURN
#define RETURN(statusptr)
H0MIN
#define H0MIN
Definition: MCInjectHoughMulti.c:71
main
int main(int argc, char *argv[])
Definition: MCInjectHoughMulti.c:103
DTERMS
#define DTERMS
Definition: MCInjectHoughMulti.c:75
F0
#define F0
Definition: MCInjectHoughMulti.c:62
NFSIZE
#define NFSIZE
Definition: MCInjectHoughMulti.c:68
PATCHSIZEX
#define PATCHSIZEX
Definition: MCInjectHoughMulti.c:66
DELTA
#define DELTA
Definition: MCInjectHoughMulti.c:65
THRESHOLD
#define THRESHOLD
Definition: MCInjectHoughMulti.c:61
GenerateInjectParams
void GenerateInjectParams(LALStatus *status, PulsarData *injectPulsar, HoughTemplate *templatePulsar, HoughNearTemplates *closeTemplates, HoughInjectParams *params, LineNoiseInfo *lines)
Definition: MCInjectHoughMulti.c:1079
EARTHEPHEMERIS
#define EARTHEPHEMERIS
Definition: MCInjectHoughMulti.c:53
DIROUT
#define DIROUT
Definition: MCInjectHoughMulti.c:79
NMCLOOP
#define NMCLOOP
Definition: MCInjectHoughMulti.c:73
NTEMPLATES
#define NTEMPLATES
Definition: MCInjectHoughMulti.c:74
TRUE
#define TRUE
Definition: MCInjectHoughMulti.c:82
FALSE
#define FALSE
Definition: MCInjectHoughMulti.c:83
NH0
#define NH0
Definition: MCInjectHoughMulti.c:70
PrintLogFile2
void PrintLogFile2(LALStatus *status, CHAR *dir, CHAR *basename, LALStringVector *linefiles, CHAR *executable)
Definition: MCInjectHoughMulti.c:1588
GenerateInjectParamsNoVeto
void GenerateInjectParamsNoVeto(LALStatus *status, PulsarData *injectPulsar, HoughTemplate *templatePulsar, HoughNearTemplates *closeTemplates, HoughInjectParams *params)
Definition: MCInjectHoughMulti.c:1306
FILEOUT
#define FILEOUT
Definition: MCInjectHoughMulti.c:80
PATCHSIZEY
#define PATCHSIZEY
Definition: MCInjectHoughMulti.c:67
BLOCKSRNGMED
#define BLOCKSRNGMED
Definition: MCInjectHoughMulti.c:69
ComputeFoft_NM
void ComputeFoft_NM(LALStatus *status, REAL8Vector *foft, HoughTemplate *pulsarTemplate, REAL8Vector *timeDiffV, REAL8Cart3CoorVector *velV)
Definition: MCInjectHoughMulti.c:1519
FBAND
#define FBAND
Definition: MCInjectHoughMulti.c:63
SUNEPHEMERIS
#define SUNEPHEMERIS
Definition: MCInjectHoughMulti.c:54
MAXFILENAMELENGTH
#define MAXFILENAMELENGTH
Definition: MCInjectHoughMulti.c:56
SFTDIRECTORY
#define SFTDIRECTORY
Definition: MCInjectHoughMulti.c:77
H0MAX
#define H0MAX
Definition: MCInjectHoughMulti.c:72
ALPHA
#define ALPHA
Definition: MCInjectHoughMulti.c:64
SFTtoUCHARPeakGram
void SFTtoUCHARPeakGram(LALStatus *status, UCHARPeakGram *pg, const SFTtype *sft, REAL8 thr)
Convert a normalized sft into a peakgram by selecting bins in which power exceeds a threshold.
Definition: PeakSelect.c:371
STRING
#define STRING(a)
fprintf
#define fprintf
e
double e
XLALDestroyMultiDetectorStateSeries
void XLALDestroyMultiDetectorStateSeries(MultiDetectorStateSeries *mdetStates)
Helper function to get rid of a multi-IFO DetectorStateSeries Note, this is "NULL-robust" in the sens...
Definition: DetectorStates.c:301
XLALGetMultiDetectorStatesFromMultiSFTs
MultiDetectorStateSeries * XLALGetMultiDetectorStatesFromMultiSFTs(const MultiSFTVector *multiSFTs, const EphemerisData *edat, REAL8 tOffset)
Get the 'detector state' (ie detector-tensor, position, velocity, etc) for the given multi-vector of ...
Definition: DetectorStates.c:537
LALGeneratePulsarSignal
void LALGeneratePulsarSignal(LALStatus *status, REAL4TimeSeries **signalvec, const PulsarSignalParams *params)
Definition: GeneratePulsarSignal.c:228
LALFastGeneratePulsarSFTs
void LALFastGeneratePulsarSFTs(LALStatus *status, SFTVector **outputSFTs, const SkyConstAndZeroPsiAMResponse *input, const SFTandSignalParams *params)
Fast generation of Fake SFTs for a pure pulsar signal.
Definition: GeneratePulsarSignal.c:573
LALSignalToSFTs
void LALSignalToSFTs(LALStatus *status, SFTVector **outputSFTs, const REAL4TimeSeries *signalvec, const SFTParams *params)
Definition: GeneratePulsarSignal.c:436
LALComputeSkyAndZeroPsiAMResponse
void LALComputeSkyAndZeroPsiAMResponse(LALStatus *status, SkyConstAndZeroPsiAMResponse *output, const SFTandSignalParams *params)
/deprecated Move to attic? Wrapper for LALComputeSky() and LALComputeDetAMResponse() that finds the s...
Definition: GeneratePulsarSignal.c:475
XLALInitBarycenter
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
Definition: LALInitBarycenter.c:206
XLALDestroyEphemerisData
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
Definition: LALInitBarycenter.c:322
XLALDestroyMultiAMCoeffs
void XLALDestroyMultiAMCoeffs(MultiAMCoeffs *multiAMcoef)
Destroy a MultiAMCoeffs structure.
Definition: LALComputeAM.c:469
XLALComputeMultiAMCoeffs
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
Definition: LALComputeAM.c:379
LAL_PI_2
#define LAL_PI_2
LAL_TWOPI
#define LAL_TWOPI
BOOLEAN
unsigned char BOOLEAN
UCHAR
unsigned char UCHAR
REAL8
double REAL8
XLAL_INIT_DECL
#define XLAL_INIT_DECL(var,...)
INT8
int64_t INT8
crectf
#define crectf(re, im)
CHAR
char CHAR
UINT4
uint32_t UINT4
COMPLEX8
float complex COMPLEX8
INT4
int32_t INT4
REAL4
float REAL4
lalDebugLevel
#define lalDebugLevel
LALHOUGHNormalizeWeights
void LALHOUGHNormalizeWeights(LALStatus *status, REAL8Vector *weightV)
Normalizes weight factors so that their sum is N.
Definition: DriveHough.c:668
LALHOUGHInitializeWeights
void LALHOUGHInitializeWeights(LALStatus *status, REAL8Vector *weightV)
Initializes weight factors to unity.
Definition: DriveHough.c:633
VTOT
#define VTOT
Total detector velocity/c TO BE CHANGED DEPENDING ON DETECTOR.
Definition: LUT.h:207
LALInvRotatePolarU
void LALInvRotatePolarU(LALStatus *status, REAL8UnitPolarCoor *out, REAL8UnitPolarCoor *in, REAL8UnitPolarCoor *par)
Definition: Stereographic.c:177
LALStereoInvProjectCart
void LALStereoInvProjectCart(LALStatus *status, REAL8UnitPolarCoor *out, REAL8Cart2Coor *in)
Definition: Stereographic.c:357
XLALNormalizeMultiSFTVect
MultiPSDVector * XLALNormalizeMultiSFTVect(MultiSFTVector *multsft, UINT4 blockSize, const MultiNoiseFloor *assumeSqrtSX)
Function for normalizing a multi vector of SFTs in a multi IFO search and returns the running-median ...
Definition: NormalizeSFTRngMed.c:159
XLALDestroyMultiPSDVector
void XLALDestroyMultiPSDVector(MultiPSDVector *multvect)
Destroy a multi PSD-vector.
Definition: PSDutils.c:102
XLALComputeMultiNoiseWeights
MultiNoiseWeights * XLALComputeMultiNoiseWeights(const MultiPSDVector *rngmed, UINT4 blocksRngMed, UINT4 excludePercentile)
Computes weight factors arising from MultiSFTs with different noise floors.
Definition: PSDutils.c:285
XLALDestroyMultiNoiseWeights
void XLALDestroyMultiNoiseWeights(MultiNoiseWeights *weights)
Destroy a MultiNoiseWeights object.
Definition: PSDutils.c:172
LALCreateRandomParams
void LALCreateRandomParams(LALStatus *status, RandomParams **params, INT4 seed)
LALUniformDeviate
void LALUniformDeviate(LALStatus *status, REAL4 *deviate, RandomParams *params)
a
static const INT4 a
LALDestroyRandomParams
void LALDestroyRandomParams(LALStatus *status, RandomParams **params)
LALRngMedBias
void LALRngMedBias(LALStatus *status, REAL8 *biasFactor, INT4 blkSize)
LALRemoveKnownLinesInMultiSFTVector
void LALRemoveKnownLinesInMultiSFTVector(LALStatus *status, MultiSFTVector *MultiSFTVect, INT4 width, INT4 window, LALStringVector *linefiles, RandomParams *randPar)
top level function to remove lines from a multi sft vector given a list of files containing list of k...
Definition: lib/SFTClean.c:953
LALCheckLines
void LALCheckLines(LALStatus *status, INT4 *countL, LineNoiseInfo *lines, REAL8 freq)
Function to count how many lines affect a given frequency.
Definition: lib/SFTClean.c:464
XLALDestroySFTVector
void XLALDestroySFTVector(SFTVector *vect)
XLAL interface to destroy an SFTVector.
Definition: SFTtypes.c:300
XLALDestroySFTCatalog
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
XLALLoadMultiSFTs
MultiSFTVector * XLALLoadMultiSFTs(const SFTCatalog *inputCatalog, REAL8 fMin, REAL8 fMax)
Function to load a catalog of SFTs from possibly different detectors.
Definition: SFTfileIO.c:416
XLALDestroyMultiSFTVector
void XLALDestroyMultiSFTVector(MultiSFTVector *multvect)
Destroy a multi SFT-vector.
Definition: SFTtypes.c:424
XLALExtractTimestampsFromSFTs
LIGOTimeGPSVector * XLALExtractTimestampsFromSFTs(const SFTVector *sfts)
Extract timstamps-vector from the given SFTVector.
Definition: SFTtimestamps.c:416
XLALCreateMultiSFTVector
MultiSFTVector * XLALCreateMultiSFTVector(UINT4 length, UINT4Vector *numsft)
Create a multi-IFO SFT vector with a given number of bins per SFT and number of SFTs per IFO (which w...
Definition: SFTtypes.c:362
XLALSFTdataFind
SFTCatalog * XLALSFTdataFind(const CHAR *file_pattern, const SFTConstraints *constraints)
Find the list of SFTs matching the file_pattern and satisfying the given constraints,...
Definition: SFTcatalog.c:71
XLALReadTimestampsFile
LIGOTimeGPSVector * XLALReadTimestampsFile(const CHAR *fname)
backwards compatible wrapper to XLALReadTimestampsFileConstrained() without GPS-time constraints
Definition: SFTtimestamps.c:297
XLALDestroyTimestampVector
void XLALDestroyTimestampVector(LIGOTimeGPSVector *vect)
De-allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:69
COORDINATESYSTEM_EQUATORIAL
COORDINATESYSTEM_EQUATORIAL
XLALUserVarReadAllInput
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
XLALDestroyUserVars
void XLALDestroyUserVars(void)
XLALUserVarGetLog
void CHAR * XLALUserVarGetLog(UserVarLogFormat format)
XLALRegisterNamedUvar
#define XLALRegisterNamedUvar(cvar, name, type, option, category,...)
XLALUserVarWasSet
int XLALUserVarWasSet(const void *cvar)
UVAR_LOGFMT_CFGFILE
UVAR_LOGFMT_CFGFILE
XLAL_CHECK_MAIN
#define XLAL_CHECK_MAIN(assertion,...)
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLALGPSSetREAL8
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
XLALGPSDiff
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)
count
long long count
Definition: hwinject.c:371
data
float data[BLOCKSIZE]
Definition: hwinject.c:360
deltaF
int deltaF
lalpulsar.git_version.status
string status
Definition: git_version.py:32
lalpulsar_PiecewiseSearch.filename
filename
Definition: lalpulsar_PiecewiseSearch.py:406
lalpulsar_PiecewiseSearch.h0
h0
Definition: lalpulsar_PiecewiseSearch.py:601
lalpulsar_PiecewiseSearch.seed
seed
Definition: lalpulsar_PiecewiseSearch.py:684
lalpulsar_knope_result_page.par
par
Definition: lalpulsar_knope_result_page.py:2446
lalpulsar_knope_result_page.f0
f0
Definition: lalpulsar_knope_result_page.py:2463
make_frame_cache.i
int i
Definition: make_frame_cache.py:77
test_computePSD.edat
edat
Definition: test_computePSD.py:65
uvar_sftDir
CHAR * uvar_sftDir
Definition: pulsar_crosscorr.c:82
uvar_startTime
REAL8 uvar_startTime
Definition: pulsar_crosscorr.c:57
uvar_f0
REAL8 uvar_f0
Definition: pulsar_crosscorr.c:58
uvar_dirnameOut
CHAR * uvar_dirnameOut
Definition: pulsar_crosscorr.c:83
uvar_blocksRngMed
INT4 uvar_blocksRngMed
Definition: pulsar_crosscorr.c:55
uvar_endTime
REAL8 uvar_endTime
Definition: pulsar_crosscorr.c:57
PIXELFACTOR
#define PIXELFACTOR
Definition: pulsar_crosscorr.h:113
AMCoeffs::b
REAL4Vector * b
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:65
AMCoeffs::a
REAL4Vector * a
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:64
COMPLEX8FrequencySeries::data
COMPLEX8Sequence * data
COMPLEX8FrequencySeriesVector::data
COMPLEX8FrequencySeries * data
Pointer to the data array.
Definition: PulsarDataTypes.h:75
COMPLEX8FrequencySeriesVector::length
UINT4 length
Number of elements in array.
Definition: PulsarDataTypes.h:74
COMPLEX8Vector::data
COMPLEX8 * data
DetectorState::vDetector
REAL8 vDetector[3]
Cart.
Definition: DetectorStates.h:130
DetectorState::tGPS
LIGOTimeGPS tGPS
GPS timestamps corresponding to this entry.
Definition: DetectorStates.h:128
DetectorStateSeries::data
DetectorState * data
array of DetectorState entries
Definition: DetectorStates.h:147
DetectorStateSeries::length
UINT4 length
total number of entries
Definition: DetectorStates.h:146
DetectorStateSeries::detector
LALDetector detector
detector-info corresponding to this timeseries
Definition: DetectorStates.h:148
EphemerisData
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
Definition: LALBarycenter.h:127
HoughNearTemplates::skytemp
REAL8UnitPolarCoor skytemp[4]
Definition: MCInjectHoughMulti.h:103
HoughTemplate::spindown
REAL8Vector spindown
Definition: MCInjectHoughMulti.h:97
LIGOTimeGPS::gpsNanoSeconds
INT4 gpsNanoSeconds
LIGOTimeGPSVector
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
LIGOTimeGPSVector::data
LIGOTimeGPS * data
array of timestamps
Definition: SFTfileIO.h:193
LIGOTimeGPSVector::length
UINT4 length
number of timestamps
Definition: SFTfileIO.h:192
LineNoiseInfo
structure for storing list of spectral lines – constructed by expanding list of harmonics
Definition: SFTClean.h:132
MultiAMCoeffs
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
Definition: LALComputeAM.h:137
MultiAMCoeffs::length
UINT4 length
number of IFOs
Definition: LALComputeAM.h:141
MultiAMCoeffs::data
AMCoeffs ** data
noise-weighted AM-coeffs , and
Definition: LALComputeAM.h:142
MultiDetectorStateSeries
Multi-IFO time-series of DetectorStates.
Definition: DetectorStates.h:154
MultiDetectorStateSeries::data
DetectorStateSeries ** data
vector of pointers to DetectorStateSeries
Definition: DetectorStates.h:159
MultiLIGOTimeGPSVector
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
Definition: SFTfileIO.h:198
MultiLIGOTimeGPSVector::length
UINT4 length
number of timestamps vectors or ifos
Definition: SFTfileIO.h:202
MultiLIGOTimeGPSVector::data
LIGOTimeGPSVector ** data
timestamps vector for each ifo
Definition: SFTfileIO.h:203
MultiNoiseWeights
One noise-weight (number) per SFT (therefore indexed over IFOs and SFTs.
Definition: PSDutils.h:71
MultiNoiseWeights::data
REAL8Vector ** data
weights-vector for each detector
Definition: PSDutils.h:76
MultiPSDVector
A collection of PSD vectors – one for each IFO in a multi-IFO search.
Definition: PSDutils.h:62
MultiSFTVector
A collection of SFT vectors – one for each IFO in a multi-IFO search.
Definition: SFTfileIO.h:179
MultiSFTVector::length
UINT4 length
number of ifos
Definition: SFTfileIO.h:183
MultiSFTVector::data
SFTVector ** data
sftvector for each ifo
Definition: SFTfileIO.h:184
PulsarData::spindown
REAL8Vector spindown
Definition: MCInjectComputeHough.h:143
PulsarSignalParams
Input parameters to GeneratePulsarSignal(), defining the source and the time-series.
Definition: GeneratePulsarSignal.h:279
REAL4TimeSeries::data
REAL4Sequence * data
REAL4Vector::data
REAL4 * data
REAL8Cart2Coor
Two dimensional Cartessian coordinates.
Definition: LUT.h:315
REAL8Cart2Coor::y
REAL8 y
Definition: LUT.h:317
REAL8Cart2Coor::x
REAL8 x
Definition: LUT.h:316
REAL8Cart3Coor
Three dimensional Cartessian coordinates.
Definition: LUT.h:308
REAL8Cart3Coor::y
REAL8 y
Definition: LUT.h:310
REAL8Cart3Coor::x
REAL8 x
Definition: LUT.h:309
REAL8Cart3Coor::z
REAL8 z
Definition: LUT.h:311
REAL8Cart3CoorVector::data
REAL8Cart3Coor * data
x.y.z
Definition: DriveHoughColor.h:148
REAL8Cart3CoorVector::length
UINT4 length
number of elements
Definition: DriveHoughColor.h:147
REAL8UnitPolarCoor
Polar coordinates of a unitary vector on the sphere.
Definition: LUT.h:327
REAL8UnitPolarCoor::alpha
REAL8 alpha
any value
Definition: LUT.h:328
REAL8UnitPolarCoor::delta
REAL8 delta
In the interval [ ].
Definition: LUT.h:329
REAL8Vector::data
REAL8 * data
SFTCatalog
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238
SFTCatalog::data
SFTDescriptor * data
array of data-entries describing matched SFTs
Definition: SFTfileIO.h:243
SFTCatalog::length
UINT4 length
number of SFTs in catalog
Definition: SFTfileIO.h:242
SFTConstraints
'Constraints' for SFT-matching: which detector, within which time-stretch and which timestamps exactl...
Definition: SFTfileIO.h:212
SFTConstraints::detector
CHAR * detector
2-char channel-prefix describing the detector (eg 'H1', 'H2', 'L1', 'G1' etc)
Definition: SFTfileIO.h:213
SFTConstraints::timestamps
LIGOTimeGPSVector * timestamps
list of timestamps
Definition: SFTfileIO.h:216
SFTConstraints::maxStartTime
LIGOTimeGPS * maxStartTime
only include SFTs whose epoch is < maxStartTime
Definition: SFTfileIO.h:215
SFTConstraints::minStartTime
LIGOTimeGPS * minStartTime
only include SFTs whose epoch is >= minStartTime
Definition: SFTfileIO.h:214
SFTDescriptor::header
SFTtype header
SFT-header info.
Definition: SFTfileIO.h:228
SFTParams
Parameters defining the SFTs to be returned from LALSignalToSFTs().
Definition: GeneratePulsarSignal.h:321
SFTParams::Tsft
REAL8 Tsft
length of each SFT in seconds
Definition: GeneratePulsarSignal.h:322
SFTParams::timestamps
const LIGOTimeGPSVector * timestamps
timestamps to produce SFTs for (can be NULL)
Definition: GeneratePulsarSignal.h:323
SFTParams::noiseSFTs
const SFTVector * noiseSFTs
noise SFTs to be added (can be NULL)
Definition: GeneratePulsarSignal.h:324
SFTandSignalParams
Parameters defining the pulsar signal and SFTs used by LALFastGeneratePulsarSFTs().
Definition: GeneratePulsarSignal.h:333
SFTandSignalParams::sinVal
REAL8 * sinVal
sinVal holds lookup table (LUT) values for doing trig sin calls
Definition: GeneratePulsarSignal.h:340
SFTandSignalParams::nSamples
INT4 nSamples
nsample from noise SFT header; 2x this equals effective number of time samples
Definition: GeneratePulsarSignal.h:336
SFTandSignalParams::pSigParams
PulsarSignalParams * pSigParams
Definition: GeneratePulsarSignal.h:334
SFTandSignalParams::cosVal
REAL8 * cosVal
cosVal holds lookup table (LUT) values for doing trig cos calls
Definition: GeneratePulsarSignal.h:341
SFTandSignalParams::trigArg
REAL8 * trigArg
array of arguments to hold lookup table (LUT) values for doing trig calls
Definition: GeneratePulsarSignal.h:339
SFTandSignalParams::Dterms
INT4 Dterms
use this to fill in SFT bins with fake data as per LALDemod else fill in bin with zero
Definition: GeneratePulsarSignal.h:337
SFTandSignalParams::resTrig
INT4 resTrig
length sinVal, cosVal; domain: -2pi to 2pi; resolution = 4pi/resTrig
Definition: GeneratePulsarSignal.h:338
SkyConstAndZeroPsiAMResponse
Sky Constants and beam pattern response functions used by LALFastGeneratePulsarSFTs().
Definition: GeneratePulsarSignal.h:349
SkyConstAndZeroPsiAMResponse::skyConst
REAL8 * skyConst
vector of A and B sky constants
Definition: GeneratePulsarSignal.h:350
SkyConstAndZeroPsiAMResponse::fPlusZeroPsi
REAL4 * fPlusZeroPsi
vector of Fplus values for psi = 0 at midpoint of each SFT
Definition: GeneratePulsarSignal.h:351
SkyConstAndZeroPsiAMResponse::fCrossZeroPsi
REAL4 * fCrossZeroPsi
vector of Fcross values for psi = 0 at midpoint of each SFT
Definition: GeneratePulsarSignal.h:352
SkyPosition::longitude
REAL8 longitude
SkyPosition::latitude
REAL8 latitude
SkyPosition::system
CoordinateSystem system
UCHARPeakGram
Explicit peakgram structure – 1 if power in bin is above threshold and 0 if below.
Definition: PeakSelect.h:114
UCHARPeakGram::data
UCHAR * data
pointer to the data {0,1}
Definition: PeakSelect.h:120
UCHARPeakGram::length
INT4 length
number of elements in data
Definition: PeakSelect.h:118
UINT4Vector::data
UINT4 * data
params::duration
double duration
params::psi
double psi
f_min
double f_min
f_max
double f_max