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
ComputeFstatTest.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2014 Reinhard Prix
3*
4* This program is free software; you can redistribute it and/or modify
5* it under the terms of the GNU General Public License as published by
6* the Free Software Foundation; either version 2 of the License, or
7* (at your option) any later version.
8*
9* This program is distributed in the hope that it will be useful,
10* but WITHOUT ANY WARRANTY; without even the implied warranty of
11* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12* GNU General Public License for more details.
13*
14* You should have received a copy of the GNU General Public License
15* along with with program; see the file COPYING. If not, write to the
16* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17* MA 02110-1301 USA
18*/
19
20#include <lal/XLALError.h>
21#include <lal/LALBarycenter.h>
22#include <lal/LALInitBarycenter.h>
23#include <lal/UniversalDopplerMetric.h>
24#include <lal/LogPrintf.h>
25#include <lal/CWMakeFakeData.h>
26#include <lal/LALConstants.h>
27#include <lal/ExtrapolatePulsarSpins.h>
28#include <lal/ComputeFstat.h>
29#include <lal/DetectorStates.h>
30#include <lal/LFTandTSutils.h>
31#include <lal/LALString.h>
32
33// basic consistency checks of the ComputeFstat module: compare F-stat results for all
34// *available* Fstat methods against each other
35
36static int compareFstatResults( const FstatResults *result1, const FstatResults *result2 );
37static int printFstatResults2File( const FstatResults *results, const char *MethodName, const INT4 iSky, const INT4 if1dot, const INT4 iPeriod, FstatQuantities whatToCompute );
38// ---------- main ----------
39int
40main( int argc, char *argv[] )
41{
42 XLAL_CHECK( argc == 1, XLAL_EINVAL, "No input arguments allowed.\n" );
43 XLAL_CHECK( argv != NULL, XLAL_EINVAL );
44
45 // ----- load ephemeris
46 EphemerisData *ephem;
47 XLAL_CHECK( ( ephem = XLALInitBarycenter( TEST_PKG_DATA_DIR "earth00-40-DE405.dat.gz", TEST_PKG_DATA_DIR "sun00-40-DE405.dat.gz" ) ) != NULL, XLAL_EFUNC );
48
49 // ----- setup injection and data parameters
50 LALStringVector *detNames = NULL;
51 XLAL_CHECK( ( detNames = XLALCreateStringVector( "H1", "L1", NULL ) ) != NULL, XLAL_EFUNC );
52 UINT4 numDetectors = detNames->length;
53
54 // generate and assume some gaussian noise floors
55 MultiNoiseFloor XLAL_INIT_DECL( injectSqrtSX );
56 MultiNoiseFloor XLAL_INIT_DECL( assumeSqrtSX );
57 injectSqrtSX.length = numDetectors;
58 assumeSqrtSX.length = numDetectors;
59 for ( UINT4 X = 0; X < numDetectors; X ++ ) {
60 injectSqrtSX.sqrtSn[X] = 0; // don't inject random noise to keep errors deterministic and informative (resampling differs much more on noise)
61 assumeSqrtSX.sqrtSn[X] = 1.0 + 2.0 * X;
62 }
63
64 LIGOTimeGPS startTime = {711595934, 0};
65 REAL8 Tspan = 20 * 3600;
66 LIGOTimeGPS endTime = startTime;
67 XLALGPSAdd( &endTime, Tspan );
68 REAL8 Tsft = 1800;
69
70 LIGOTimeGPS refTime = { startTime.gpsSeconds - 2.3 * Tspan, 0 };
71
73 XLAL_CHECK( ( multiTimestamps = XLALCalloc( 1, sizeof( *multiTimestamps ) ) ) != NULL, XLAL_ENOMEM );
74 XLAL_CHECK( ( multiTimestamps->data = XLALCalloc( numDetectors, sizeof( multiTimestamps->data[0] ) ) ) != NULL, XLAL_ENOMEM );
76 LIGOTimeGPS startTimeX = startTime;
77 for ( UINT4 X = 0; X < numDetectors; X ++ ) {
78 XLAL_CHECK( ( multiTimestamps->data[X] = XLALMakeTimestamps( startTimeX, Tspan, Tsft, 0 ) ) != NULL, XLAL_EFUNC );
79 XLALGPSAdd( &startTimeX, 0.5 * Tspan ); // shift start-times by 1/2 Tspan for each detector
80 Tspan *= 2.0;
81 } // for X < numDetectors
82
83 // shift a few timestamps around to create gaps
84 UINT4 numSFTsPerDet = multiTimestamps->data[0]->length;
85 multiTimestamps->data[0]->data[numSFTsPerDet - 1].gpsSeconds += 10000;
86 multiTimestamps->data[0]->data[numSFTsPerDet - 2].gpsSeconds += 5000;
87 multiTimestamps->data[1]->data[0].gpsSeconds -= 10000;
88 multiTimestamps->data[1]->data[1].gpsSeconds -= 5000;
89
90 SFTCatalog *catalog;
91 XLAL_CHECK( ( catalog = XLALMultiAddToFakeSFTCatalog( NULL, detNames, multiTimestamps ) ) != NULL, XLAL_EFUNC );
92
93 // ----- CW sources to injet ----------
94 REAL8 Freq = 100.0;
95 REAL8 h0 = 1.0;
96 REAL8 cosi = 0.5;
97
98 PulsarParamsVector *injectSources;
99 XLAL_CHECK( ( injectSources = XLALCreatePulsarParamsVector( 1 ) ) != NULL, XLAL_EFUNC );
100
101 injectSources->data[0].Amp.aPlus = 0.5 * h0 * ( 1.0 + cosi * cosi );
102 injectSources->data[0].Amp.aCross = h0 * cosi;
103 injectSources->data[0].Amp.psi = 0.1;
104 injectSources->data[0].Amp.phi0 = 1.2;
105
106 REAL8 asini = 0; // 1.4; // sco-X1 like
107 REAL8 Period = 0; // 19 * 3600;// sco-X1 like
108 REAL8 ecc = 0; // 0.1; // much larger than ScoX1
110 Doppler.Alpha = 0.5;
111 Doppler.Delta = -0.5;
112 Doppler.fkdot[0] = Freq;
113 Doppler.fkdot[1] = -1e-9;
114 Doppler.refTime = refTime;
115
116 Doppler.asini = asini;
117 Doppler.ecc = ecc;
118 Doppler.tp = startTime;
119 Doppler.period = Period;
120 Doppler.argp = 0.5;
121
122 injectSources->data[0].Doppler = Doppler;
123
124 REAL8 dFreq = 0.1 / Tspan; // 10x finer than native FFT resolution
125 REAL8 mis = 0.5;
126 REAL8 df1dot = sqrt( 720.0 * mis ) / ( LAL_PI * Tspan * Tspan ); // metric (f-projected) stepsize for given mismatch mis
127 REAL8 dSky = 1e4 / ( Freq * Tspan ); // rough estimate of a 'metric' sky step, eg. Eq.(118) in \cite Prix07
128
129 REAL8 dPeriod = 3600;
130 UINT4 numFreqBins = 1000;
131
132 UINT4 numf1dotPoints = 2;
133 UINT4 numSkyPoints = 2;
134 UINT4 numPeriodPoints = 2;
135
136 PulsarSpinRange XLAL_INIT_DECL( spinRange );
137 spinRange.refTime = refTime;
138 memcpy( &spinRange.fkdot, &injectSources->data[0].Doppler.fkdot, sizeof( spinRange.fkdot ) );
139 spinRange.fkdotBand[0] = ( numFreqBins - 1 ) * dFreq - 10 * LAL_REAL8_EPS;
140 spinRange.fkdotBand[1] = ( numf1dotPoints - 1 ) * df1dot - 10 * LAL_REAL8_EPS;
141
142 Doppler.fkdot[0] -= 0.4 * spinRange.fkdotBand[0];
143
144 REAL8 minCoverFreq, maxCoverFreq;
145 XLAL_CHECK( XLALCWSignalCoveringBand( &minCoverFreq, &maxCoverFreq, &startTime, &endTime, &spinRange, asini, Period, ecc ) == XLAL_SUCCESS, XLAL_EFUNC );
146
147 // ----- setup optional Fstat arguments
149 optionalArgs.injectSources = injectSources;
150 optionalArgs.injectSqrtSX = &injectSqrtSX;
151 optionalArgs.assumeSqrtSX = &assumeSqrtSX;
152
153 // ----- prepare input data with injection for all available methods
154 FstatInput *input_seg1[FMETHOD_END], *input_seg2[FMETHOD_END];
155 FstatResults *results_seg1[FMETHOD_END], *results_seg2[FMETHOD_END];
156 for ( UINT4 iMethod = FMETHOD_START; iMethod < FMETHOD_END; iMethod ++ ) {
157 results_seg1[iMethod] = results_seg2[iMethod] = NULL;
158 if ( !XLALFstatMethodIsAvailable( iMethod ) ) {
159 continue;
160 }
161 optionalArgs.FstatMethod = iMethod;
162 optionalArgs.prevInput = NULL;
163 optionalArgs.resampFFTPowerOf2 = ( 1 == 1 );
164 XLAL_CHECK( ( input_seg1[iMethod] = XLALCreateFstatInput( catalog, minCoverFreq, maxCoverFreq, dFreq, ephem, &optionalArgs ) ) != NULL, XLAL_EFUNC );
165 optionalArgs.prevInput = input_seg1[iMethod];
166 optionalArgs.resampFFTPowerOf2 = ( 1 == 0 );
167 XLAL_CHECK( ( input_seg2[iMethod] = XLALCreateFstatInput( catalog, minCoverFreq - 0.01, maxCoverFreq + 0.01, dFreq, ephem, &optionalArgs ) ) != NULL, XLAL_EFUNC );
168 }
169
170
171
172 FstatQuantities whatToCompute = ( FSTATQ_2F | FSTATQ_FAFB );
173 // ----- loop over all templates {sky, f1dot, period}
174 for ( UINT4 iSky = 0; iSky < numSkyPoints; iSky ++ ) {
175 for ( UINT4 if1dot = 0; if1dot < numf1dotPoints; if1dot ++ ) {
176 for ( UINT4 iPeriod = 0; iPeriod < numPeriodPoints; iPeriod ++ ) {
177 MultiCOMPLEX8TimeSeries *first_SRC_a = NULL;
178 MultiCOMPLEX8TimeSeries *first_SRC_b = NULL;
179
180 // ----- loop over all available methods and compare Fstat results
181 FstatMethodType firstMethod = FMETHOD_START;
182 for ( UINT4 iMethod = FMETHOD_START; iMethod < FMETHOD_END; iMethod ++ ) {
183 if ( !XLALFstatMethodIsAvailable( iMethod ) ) {
184 continue;
185 }
186 if ( firstMethod == FMETHOD_START ) { // keep track of first available method found
187 firstMethod = iMethod;
188 }
189 if ( ( iMethod == FMETHOD_DEMOD_BEST ) || ( iMethod == FMETHOD_RESAMP_BEST ) ) {
190 continue; // avoid re-running comparisons for same method because labelled 'best'
191 }
192
193 XLAL_CHECK( XLALComputeFstat( &results_seg1[iMethod], input_seg1[iMethod], &Doppler, numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
194 XLAL_CHECK( XLALComputeFstat( &results_seg2[iMethod], input_seg2[iMethod], &Doppler, numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
195
196 if ( lalDebugLevel & LALINFOBIT ) {
197 printFstatResults2File( results_seg1[iMethod], XLALGetFstatInputMethodName( input_seg1[iMethod] ), iSky, if1dot, iPeriod, whatToCompute );
198 } // if info
199
200 // compare to first result
201 if ( iMethod != firstMethod ) {
202 XLALPrintInfo( "Comparing results between method '%s' and '%s'\n", XLALGetFstatInputMethodName( input_seg1[firstMethod] ), XLALGetFstatInputMethodName( input_seg1[iMethod] ) );
203 if ( compareFstatResults( results_seg1[firstMethod], results_seg1[iMethod] ) != XLAL_SUCCESS ) {
204 XLALPrintError( "Comparison between method '%s' and '%s' failed on 'seg1'\n", XLALGetFstatInputMethodName( input_seg1[firstMethod] ), XLALGetFstatInputMethodName( input_seg1[iMethod] ) );
206 }
207 if ( compareFstatResults( results_seg2[firstMethod], results_seg2[iMethod] ) != XLAL_SUCCESS ) {
208 XLALPrintError( "Comparison between method '%s' and '%s' failed on 'seg2'\n", XLALGetFstatInputMethodName( input_seg1[firstMethod] ), XLALGetFstatInputMethodName( input_seg1[iMethod] ) );
210 }
211
212 // for resampling methods, check time series extraction and consistency
213 if ( iMethod >= FMETHOD_RESAMP_GENERIC ) {
214 if ( first_SRC_a == NULL ) {
215 XLAL_CHECK( XLALExtractResampledTimeseries( &first_SRC_a, &first_SRC_b, input_seg2[iMethod] ) == XLAL_SUCCESS, XLAL_EFUNC );
216 XLAL_CHECK( first_SRC_a != NULL, XLAL_EFAULT );
217 XLAL_CHECK( first_SRC_b != NULL, XLAL_EFAULT );
218 XLAL_CHECK( first_SRC_a->length == first_SRC_b->length, XLAL_EFAILED );
219 for ( UINT4 X = 0; X < first_SRC_a->length; ++X ) {
220 XLAL_CHECK( first_SRC_a->data[X] != NULL, XLAL_EFAULT );
221 XLAL_CHECK( first_SRC_b->data[X] != NULL, XLAL_EFAULT );
222 XLAL_CHECK( first_SRC_a->data[X]->data != NULL, XLAL_EFAULT );
223 XLAL_CHECK( first_SRC_b->data[X]->data != NULL, XLAL_EFAULT );
224 XLAL_CHECK( XLALGPSCmp( &first_SRC_a->data[X]->epoch, &first_SRC_b->data[X]->epoch ) == 0, XLAL_EFAILED );
225 XLAL_CHECK( fabs( first_SRC_a->data[X]->f0 - first_SRC_b->data[X]->f0 ) < 10 * LAL_REAL8_EPS, XLAL_EFAILED );
226 XLAL_CHECK( fabs( first_SRC_a->data[X]->deltaT - first_SRC_b->data[X]->deltaT ) < 10 * LAL_REAL8_EPS, XLAL_EFAILED );
227 XLAL_CHECK( first_SRC_a->data[X]->data->data != NULL, XLAL_EFAULT );
228 XLAL_CHECK( first_SRC_b->data[X]->data->data != NULL, XLAL_EFAULT );
229 XLAL_CHECK( first_SRC_a->data[X]->data->length == first_SRC_b->data[X]->data->length, XLAL_EFAILED );
230 }
231 } else {
232 MultiCOMPLEX8TimeSeries *SRC_a = NULL;
233 MultiCOMPLEX8TimeSeries *SRC_b = NULL;
234 XLAL_CHECK( XLALExtractResampledTimeseries( &SRC_a, &SRC_b, input_seg2[iMethod] ) == XLAL_SUCCESS, XLAL_EFUNC );
235 XLAL_CHECK( SRC_a != NULL, XLAL_EFAULT );
236 XLAL_CHECK( SRC_b != NULL, XLAL_EFAULT );
237 XLAL_CHECK( SRC_a->length == first_SRC_a->length, XLAL_EFAILED );
238 XLAL_CHECK( SRC_b->length == first_SRC_b->length, XLAL_EFAILED );
239 for ( UINT4 X = 0; X < first_SRC_a->length; ++X ) {
240 XLAL_CHECK( SRC_a->data[X] != NULL, XLAL_EFAULT );
241 XLAL_CHECK( SRC_b->data[X] != NULL, XLAL_EFAULT );
242 XLAL_CHECK( SRC_a->data[X]->data != NULL, XLAL_EFAULT );
243 XLAL_CHECK( SRC_b->data[X]->data != NULL, XLAL_EFAULT );
244 XLAL_CHECK( XLALGPSCmp( &SRC_a->data[X]->epoch, &first_SRC_a->data[X]->epoch ) == 0, XLAL_EFAILED );
245 XLAL_CHECK( fabs( SRC_a->data[X]->f0 - first_SRC_a->data[X]->f0 ) < 10 * LAL_REAL8_EPS, XLAL_EFAILED );
246 XLAL_CHECK( fabs( SRC_a->data[X]->deltaT - first_SRC_a->data[X]->deltaT ) < 10 * LAL_REAL8_EPS, XLAL_EFAILED );
247 XLAL_CHECK( SRC_a->data[X]->data->data != NULL, XLAL_EFAULT );
248 XLAL_CHECK( SRC_b->data[X]->data->data != NULL, XLAL_EFAULT );
249 XLAL_CHECK( SRC_a->data[X]->data->length == first_SRC_a->data[X]->data->length, XLAL_EFAILED );
250 VectorComparison tol = { .relErr_L1 = 5e-5, .relErr_L2 = 5e-5, .angleV = 5e-5, .relErr_atMaxAbsx = 5e-5, .relErr_atMaxAbsy = 5e-5 };
252 XLAL_CHECK( XLALCompareCOMPLEX8Vectors( &cmp, SRC_a->data[X]->data, first_SRC_a->data[X]->data, &tol ) == XLAL_SUCCESS, XLAL_EFUNC );
253 }
254 }
255 }
256
257 }
258
259 } // for i < FMETHOD_END
260
261 Doppler.period += dPeriod;
262
263 } // for iPeriod < numPeriodPoints
264
265 Doppler.fkdot[1] += df1dot;
266
267 } // for if1dot < numf1dotPoints
268
269 Doppler.Alpha += dSky;
270
271 } // for iSky < numSkyPoints
272
273 // ----- test XLALFstatInputTimeslice()
274 // setup optional Fstat arguments
275 optionalArgs.FstatMethod = FMETHOD_DEMOD_BEST; // only use demod best
276 optionalArgs.prevInput = NULL;
277 // find overlap of all detectors
278 LIGOTimeGPS startTime_slice = multiTimestamps->data[0]->data[0];
279 LIGOTimeGPS endTime_slice = multiTimestamps->data[0]->data[multiTimestamps->data[0]->length - 1];
280 for ( UINT4 X = 1; X < numDetectors; X++ ) {
281 if ( XLALGPSCmp( &startTime_slice, &multiTimestamps->data[X]->data[0] ) == -1 ) {
282 startTime_slice = multiTimestamps->data[X]->data[0];
283 }
284 if ( XLALGPSCmp( &endTime_slice, &multiTimestamps->data[X]->data[multiTimestamps->data[X]->length - 1] ) == 1 ) {
285 endTime_slice = multiTimestamps->data[X]->data[multiTimestamps->data[X]->length - 1];
286 }
287 } // for X < numDetectors
288
289 // ----- prepare input data with XLALFstatInputTimeslice and XLALSFTCatalogTimeslice
290 FstatResults *results_input_slice = NULL, *results_sft_slice = NULL;
291 FstatInput *input_sft_slice = NULL, *input_slice = NULL;
292 // slice catalog and create FstatInput structure
293 SFTCatalog *catalog_slice = NULL;
294 XLAL_CHECK( ( catalog_slice = XLALCalloc( 1, sizeof( *catalog ) ) ) != NULL, XLAL_ENOMEM );
295 XLAL_CHECK( ( XLALSFTCatalogTimeslice( catalog_slice, catalog, &startTime_slice, &endTime_slice ) ) == XLAL_SUCCESS, XLAL_EFUNC );
296 XLAL_CHECK( ( input_sft_slice = XLALCreateFstatInput( catalog_slice, minCoverFreq, maxCoverFreq, dFreq, ephem, &optionalArgs ) ) != NULL, XLAL_EFUNC );
297 // generate a timeslice of the FstatInput structure from previous test
298 XLAL_CHECK( ( XLALFstatInputTimeslice( &input_slice, input_seg1[FMETHOD_DEMOD_BEST], &startTime_slice, &endTime_slice ) ) == XLAL_SUCCESS, XLAL_EFUNC );
299
300 // Compute Fstatistic for both
301 XLAL_CHECK( XLALComputeFstat( &results_sft_slice, input_sft_slice, &Doppler, numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
302 XLAL_CHECK( XLALComputeFstat( &results_input_slice, input_slice, &Doppler, numFreqBins, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
303
304 if ( lalDebugLevel & LALINFOBIT ) {
305 printFstatResults2File( results_input_slice, "FstatInputTimeslice", numSkyPoints, numf1dotPoints, numPeriodPoints, whatToCompute );
306 printFstatResults2File( results_sft_slice, "SFTCatalogTimeslice", numSkyPoints, numf1dotPoints, numPeriodPoints, whatToCompute );
307 } // if info
308
309 XLALPrintInfo( "Comparing results between SFTCatalogTimeslice and FstatInputTimeslice\n" );
310 if ( compareFstatResults( results_sft_slice, results_input_slice ) != XLAL_SUCCESS ) {
311 XLALPrintError( "Comparison between SFTCatalogTimeslice and FstatInputTimeslice failed\n" );
313 }
314
315 // free remaining memory
316 for ( UINT4 iMethod = FMETHOD_START; iMethod < FMETHOD_END; iMethod ++ ) {
317 if ( !XLALFstatMethodIsAvailable( iMethod ) ) {
318 continue;
319 }
320 XLALDestroyFstatInput( input_seg1[iMethod] );
321 XLALDestroyFstatInput( input_seg2[iMethod] );
322 XLALDestroyFstatResults( results_seg1[iMethod] );
323 XLALDestroyFstatResults( results_seg2[iMethod] );
324 } // for i < FMETHOD_END
325
326 XLALDestroyFstatInput( input_slice );
327 XLALDestroyFstatInput( input_sft_slice );
328 XLALDestroyFstatResults( results_input_slice );
329 XLALDestroyFstatResults( results_sft_slice );
330 XLALFree( catalog_slice );
331
332 XLALDestroyPulsarParamsVector( injectSources );
333 XLALDestroySFTCatalog( catalog );
335 XLALDestroyStringVector( detNames );
337
339
340 return XLAL_SUCCESS;
341
342} // main()
343
344static int
345compareFstatResults( const FstatResults *result1, const FstatResults *result2 )
346{
347 XLAL_CHECK( ( result1 != NULL ) && ( result2 != NULL ), XLAL_EINVAL );
348
349 XLAL_CHECK( result1->whatWasComputed == result2->whatWasComputed, XLAL_EINVAL );
350 XLAL_CHECK( result1->dFreq == result2->dFreq, XLAL_EINVAL );
351 XLAL_CHECK( result1->numFreqBins == result2->numFreqBins, XLAL_EINVAL );
352 XLAL_CHECK( result1->numDetectors == result2->numDetectors, XLAL_EINVAL );
353
354 // ----- set tolerance levels for comparisons ----------
356 tol.relErr_L1 = 2.5e-2;
357 tol.relErr_L2 = 2.2e-2;
358 tol.angleV = 0.02; // rad
359 tol.relErr_atMaxAbsx = 2.1e-2;
360 tol.relErr_atMaxAbsy = 2.1e-2;
361
362 UINT4 numFreqBins = result1->numFreqBins;
364
365 if ( result1->whatWasComputed & FSTATQ_2F ) {
366 // ----- package twoF arrays into REAL4Vectors
369 v1.length = numFreqBins;
370 v2.length = numFreqBins;
371 v1.data = result1->twoF;
372 v2.data = result2->twoF;
373
374 XLALPrintInfo( "Comparing 2F values:\n" );
376
377 // test comparison sanity with identical vectors should yield 0 differences
379 // we saw 6.0e-09 on arm64
380 tol0.angleV = 5e-8;
381 XLAL_CHECK( XLALCompareREAL4Vectors( &cmp, &v1, &v1, &tol0 ) == XLAL_SUCCESS, XLAL_EFUNC );
382 XLAL_CHECK( XLALCompareREAL4Vectors( &cmp, &v2, &v2, &tol0 ) == XLAL_SUCCESS, XLAL_EFUNC );
383 }
384
385 if ( result1->whatWasComputed & FSTATQ_FAFB ) {
386 // ----- package Fa,Fb arrays int COMPLEX8Vectors
389 c1.length = numFreqBins;
390 c2.length = numFreqBins;
391 // Fa
392 c1.data = result1->Fa;
393 c2.data = result2->Fa;
394 XLALPrintInfo( "Comparing Fa values:\n" );
395 XLAL_CHECK( XLALCompareCOMPLEX8Vectors( &cmp, &c1, &c2, &tol ) == XLAL_SUCCESS, XLAL_EFUNC ); // FIXME: deactivated test
396 // Fb
397 c1.data = result1->Fb;
398 c2.data = result2->Fb;
399 XLALPrintInfo( "Comparing Fb values:\n" );
400 XLAL_CHECK( XLALCompareCOMPLEX8Vectors( &cmp, &c1, &c2, &tol ) == XLAL_SUCCESS, XLAL_EFUNC ); // FIXME: deactivated test
401 }
402
403 return XLAL_SUCCESS;
404
405} // compareFstatResults()
406
407static int
408printFstatResults2File( const FstatResults *results, const char *MethodName, const INT4 iSky, const INT4 if1dot, const INT4 iPeriod, FstatQuantities whatToCompute )
409{
410 FILE *fp;
411 char fname[1024];
412 XLAL_INIT_MEM( fname );
413 snprintf( fname, sizeof( fname ) - 1, "twoF%s-iSky%02d-if1dot%02d-iPeriod%02d.dat", MethodName, iSky, if1dot, iPeriod );
414 XLAL_CHECK( ( fp = fopen( fname, "wb" ) ) != NULL, XLAL_EFUNC );
415 for ( UINT4 k = 0; k < results->numFreqBins; k ++ ) {
416 REAL8 Freq0 = results->doppler.fkdot[0];
417 REAL8 Freq_k = Freq0 + k * results->dFreq;
418 if ( whatToCompute & FSTATQ_FAFB ) {
419 fprintf( fp, "%20.16g %10.4g %10.4g %10.4g %10.4g %10.4g\n",
420 Freq_k, results->twoF[k],
421 crealf( results->Fa[k] ), cimagf( results->Fa[k] ),
422 crealf( results->Fb[k] ), cimagf( results->Fb[k] )
423 );
424 } else {
425 fprintf( fp, "%20.16g %10.4g\n",
426 Freq_k, results->twoF[k] );
427 }
428 } // for k < numFreqBins
429 fclose( fp );
430
431 return XLAL_SUCCESS;
432} // printFstatResults2File()
int main(int argc, char *argv[])
static int compareFstatResults(const FstatResults *result1, const FstatResults *result2)
static int printFstatResults2File(const FstatResults *results, const char *MethodName, const INT4 iSky, const INT4 if1dot, const INT4 iPeriod, FstatQuantities whatToCompute)
int k
void LALCheckMemoryLeaks(void)
const double c1
const double c2
REAL8 tol
#define fprintf
double e
PulsarParamsVector * XLALCreatePulsarParamsVector(UINT4 numPulsars)
Create zero-initialized PulsarParamsVector for numPulsars.
void XLALDestroyPulsarParamsVector(PulsarParamsVector *ppvect)
Destructor for PulsarParamsVector type.
FstatMethodType
Different methods available to compute the -statistic, falling into two broad classes:
Definition: ComputeFstat.h:111
const CHAR * XLALGetFstatInputMethodName(const FstatInput *input)
Returns the human-readable name of the -statistic method being used by a FstatInput structure.
Definition: ComputeFstat.c:672
const FstatOptionalArgs FstatOptionalArgsDefaults
Global initializer for setting FstatOptionalArgs to default values.
Definition: ComputeFstat.c:90
int XLALFstatInputTimeslice(FstatInput **slice, const FstatInput *input, const LIGOTimeGPS *minStartGPS, const LIGOTimeGPS *maxStartGPS)
Create and return an FstatInput 'timeslice' for given input FstatInput object [must be using LALDemod...
int XLALExtractResampledTimeseries(MultiCOMPLEX8TimeSeries **multiTimeSeries_SRC_a, MultiCOMPLEX8TimeSeries **multiTimeSeries_SRC_b, FstatInput *input)
Extracts the resampled timeseries from a given FstatInput structure (must have been initialized previ...
void XLALDestroyFstatInput(FstatInput *input)
Free all memory associated with a FstatInput structure.
Definition: ComputeFstat.c:885
FstatQuantities
Bit-field of -statistic quantities which can be computed by XLALComputeFstat().
Definition: ComputeFstat.h:94
FstatInput * XLALCreateFstatInput(const SFTCatalog *SFTcatalog, const REAL8 minCoverFreq, const REAL8 maxCoverFreq, const REAL8 dFreq, const EphemerisData *ephemerides, const FstatOptionalArgs *optionalArgs)
Create a fully-setup FstatInput structure for computing the -statistic using XLALComputeFstat().
Definition: ComputeFstat.c:359
int XLALComputeFstat(FstatResults **Fstats, FstatInput *input, const PulsarDopplerParams *doppler, const UINT4 numFreqBins, const FstatQuantities whatToCompute)
Compute the -statistic over a band of frequencies.
Definition: ComputeFstat.c:760
int XLALFstatMethodIsAvailable(FstatMethodType method)
Return true if given FstatMethodType corresponds to a valid and available Fstat method,...
void XLALDestroyFstatResults(FstatResults *Fstats)
Free all memory associated with a FstatResults structure.
Definition: ComputeFstat.c:934
@ FMETHOD_RESAMP_GENERIC
Resamp: generic implementation
Definition: ComputeFstat.h:123
@ FMETHOD_DEMOD_BEST
Demod: best guess of the fastest available hotloop
Definition: ComputeFstat.h:121
@ FMETHOD_RESAMP_BEST
Resamp: best guess of the fastest available implementation
Definition: ComputeFstat.h:125
@ FSTATQ_2F
Compute multi-detector .
Definition: ComputeFstat.h:96
@ FSTATQ_FAFB
Compute multi-detector and .
Definition: ComputeFstat.h:97
int XLALCWSignalCoveringBand(REAL8 *minCoverFreq, REAL8 *maxCoverFreq, const LIGOTimeGPS *time1, const LIGOTimeGPS *time2, const PulsarSpinRange *spinRange, const REAL8 binaryMaxAsini, const REAL8 binaryMinPeriod, const REAL8 binaryMaxEcc)
Determines a frequency band which covers the frequency evolution of a band of CW signals between two ...
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
#define LAL_REAL8_EPS
#define LAL_PI
#define XLAL_INIT_MEM(x)
double REAL8
#define XLAL_INIT_DECL(var,...)
uint32_t UINT4
int32_t INT4
#define lalDebugLevel
LALINFOBIT
void * XLALCalloc(size_t m, size_t n)
void XLALFree(void *p)
int XLALCompareREAL4Vectors(VectorComparison *result, const REAL4Vector *x, const REAL4Vector *y, const VectorComparison *tol)
Compare two REAL4 vectors using various different comparison metrics.
int XLALCompareCOMPLEX8Vectors(VectorComparison *result, const COMPLEX8Vector *x, const COMPLEX8Vector *y, const VectorComparison *tol)
Compare two COMPLEX8 vectors using various different comparison metrics.
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
LIGOTimeGPSVector * XLALMakeTimestamps(LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap)
Given a start-time, Tspan, Tsft and Toverlap, returns a list of timestamps covering this time-stretch...
SFTCatalog * XLALMultiAddToFakeSFTCatalog(SFTCatalog *catalog, const LALStringVector *detectors, const MultiLIGOTimeGPSVector *timestamps)
Multi-detector and multi-timestamp wrapper of XLALAddToFakeSFTCatalog().
Definition: SFTcatalog.c:749
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
int XLALSFTCatalogTimeslice(SFTCatalog *slice, const SFTCatalog *catalog, const LIGOTimeGPS *minStartGPS, const LIGOTimeGPS *maxStartGPS)
Set a SFT catalog 'slice' to a timeslice of a larger SFT catalog 'catalog', with entries restricted t...
Definition: SFTcatalog.c:624
void XLALDestroyStringVector(LALStringVector *vect)
LALStringVector * XLALCreateStringVector(const CHAR *str1,...)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_ERROR(...)
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFAULT
XLAL_EFUNC
XLAL_EINVAL
XLAL_EFAILED
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
size_t numDetectors
LIGOTimeGPS epoch
COMPLEX8Sequence * data
COMPLEX8 * data
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
Struct of optional 'advanced level' and (potentially method-specific) arguments to be passed to the ...
Definition: ComputeFstat.h:137
MultiNoiseFloor * injectSqrtSX
Single-sided PSD values for fake Gaussian noise to be added to SFT data.
Definition: ComputeFstat.h:144
PulsarParamsVector * injectSources
Vector of parameters of CW signals to simulate and inject.
Definition: ComputeFstat.h:143
MultiNoiseFloor * assumeSqrtSX
Single-sided PSD values to be used for computing SFT noise weights instead of from a running median o...
Definition: ComputeFstat.h:145
BOOLEAN resampFFTPowerOf2
Resamp: round up FFT lengths to next power of 2; see FstatMethodType.
Definition: ComputeFstat.h:148
FstatInput * prevInput
An FstatInput structure from a previous call to XLALCreateFstatInput(); may contain common workspace ...
Definition: ComputeFstat.h:146
FstatMethodType FstatMethod
Method to use for computing the -statistic.
Definition: ComputeFstat.h:142
XLALComputeFstat() computed results structure.
Definition: ComputeFstat.h:202
COMPLEX8 * Fb
Definition: ComputeFstat.h:260
UINT4 numDetectors
Number of detectors over which the were computed.
Definition: ComputeFstat.h:221
FstatQuantities whatWasComputed
Bit-field of which -statistic quantities were computed.
Definition: ComputeFstat.h:234
REAL8 dFreq
Spacing in frequency between each computed -statistic.
Definition: ComputeFstat.h:214
UINT4 numFreqBins
Number of frequencies at which the were computed.
Definition: ComputeFstat.h:217
REAL4 * twoF
If whatWasComputed & FSTATQ_2F is true, the multi-detector values computed at numFreqBins frequencie...
Definition: ComputeFstat.h:242
COMPLEX8 * Fa
If whatWasComputed & FSTATQ_PARTS is true, the multi-detector and computed at numFreqBins frequenci...
Definition: ComputeFstat.h:259
PulsarDopplerParams doppler
Doppler parameters, including the starting frequency, at which the were computed.
Definition: ComputeFstat.h:206
Multi-IFO container for COMPLEX8 resampled timeseries.
Definition: LFTandTSutils.h:52
COMPLEX8TimeSeries ** data
array of COMPLEX8TimeSeries (pointers)
Definition: LFTandTSutils.h:57
UINT4 length
number of IFOs
Definition: LFTandTSutils.h:56
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
Definition: SFTfileIO.h:198
array of detector-specific 'noise floors' (ie PSD values), assumed constant over the frequency-band o...
REAL8 aCross
Signal amplitude (cross)
REAL8 psi
polarization angle psi
REAL8 aPlus
Signal amplitude (plus)
REAL8 phi0
initial signal-phase (at some reference time)
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ... ] where fkdot = d^kFreq/dt^k.
PulsarAmplitudeParams Amp
'Amplitude-parameters': h0, cosi, phi0, psi
PulsarDopplerParams Doppler
'Phase-evolution parameters': {skypos, fkdot, orbital params }
Straightforward vector type of N PulsarParams structs.
PulsarParams * data
array of pulsar-signal parameters
Contains a "spin-range", ie spins and corresponding bands at a given (SSB) reference GPS-time .
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238
Struct holding the results of comparing two floating-point vectors (real-valued or complex),...
Definition: LFTandTSutils.h:64