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