Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-b246709
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
TwoSpect.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2010 -- 2014 Evan Goetz
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 * \defgroup lalpulsar_bin_TwoSpect TwoSpect Search Application
22 * \ingroup lalpulsar_bin_Apps
23 */
24
25/**
26 * \file
27 * \ingroup lalpulsar_bin_TwoSpect
28 * \author Evan Goetz
29 */
30
31#include "config.h"
32#include <sys/stat.h>
33
34#include <lal/UserInput.h>
35#include <lal/LALString.h>
36#include <lal/Window.h>
37#include <lal/DopplerScan.h>
38#include <lal/LALPulsarVCSInfo.h>
39
40#include <gsl/gsl_math.h>
41
42#include "SFTfunctions.h"
43#include "IHS.h"
44#include "candidates.h"
45#include "antenna.h"
46#include "templates.h"
47#include "TwoSpect.h"
48#include "statistics.h"
49#include "upperlimits.h"
50#include "vectormath.h"
51#include "falsealarm.h"
52
53
54//Global variables
55FILE *LOG = NULL;
56
57/**
58 * Main program
59 */
60int main( int argc, char *argv[] )
61{
62 INT4 ii, jj; //counter variables
64 time_t programstarttime, programendtime;
65 struct tm *ptm;
66
67 time( &programstarttime );
68 ptm = localtime( &programstarttime );
69
70 //Turn off gsl error handler
71 gsl_set_error_handler_off();
72
73 //Set up shop and read in the user input values
76
77 //Create directory
78 INT4 dirstatus = mkdir( uvar.outdirectory, 0777 );
79 XLAL_CHECK( dirstatus == 0 || ( dirstatus == -1 && errno == EEXIST ), XLAL_EIO, "Couldn't create directory %s\n", uvar.outdirectory ) ;
80
81 //Filenames for logfile and ULfile
82 CHARVector *LOGFILENAME = NULL, *ULFILENAME = NULL, *CANDFILENAME = NULL;
83 XLAL_CHECK( ( LOGFILENAME = XLALCreateCHARVector( strlen( uvar.outdirectory ) + strlen( uvar.outfilename ) + 3 ) ) != NULL, XLAL_EFUNC );
84 sprintf( LOGFILENAME->data, "%s/%s", uvar.outdirectory, uvar.outfilename );
85 XLAL_CHECK( ( ULFILENAME = XLALCreateCHARVector( strlen( uvar.outdirectory ) + strlen( uvar.ULfilename ) + 3 ) ) != NULL, XLAL_EFUNC );
86 sprintf( ULFILENAME->data, "%s/%s", uvar.outdirectory, uvar.ULfilename );
87 XLAL_CHECK( ( CANDFILENAME = XLALCreateCHARVector( strlen( uvar.outdirectory ) + strlen( uvar.candidatesFilename ) + 3 ) ) != NULL, XLAL_EFUNC );
88 sprintf( CANDFILENAME->data, "%s/%s", uvar.outdirectory, uvar.candidatesFilename );
89
90 //Save args_info
91 CHARVector *configfilecopy = NULL;
92 XLAL_CHECK( ( configfilecopy = XLALCreateCHARVector( strlen( uvar.outdirectory ) + strlen( uvar.configCopy ) + 3 ) ) != NULL, XLAL_EFUNC );
93 sprintf( configfilecopy->data, "%s/%s", uvar.outdirectory, uvar.configCopy );
94 FILE *INPUTVALS = NULL;
95 XLAL_CHECK( ( INPUTVALS = fopen( configfilecopy->data, "w" ) ) != NULL, XLAL_EIO, "Failed to fopen %s for writing input parameter values\n", configfilecopy->data );
96 CHAR *cfgFileStringCopy = NULL;
97 XLAL_CHECK( ( cfgFileStringCopy = XLALUserVarGetLog( UVAR_LOGFMT_CFGFILE ) ) != NULL, XLAL_EFUNC );
98 fprintf( INPUTVALS, "%s", cfgFileStringCopy );
99 fclose( INPUTVALS );
100 XLALFree( cfgFileStringCopy );
101 XLALDestroyCHARVector( configfilecopy );
102
103 //Open log file
104 XLAL_CHECK( ( LOG = fopen( LOGFILENAME->data, "w" ) ) != NULL, XLAL_EIO, "Failed to fopen %s for writing log file\n", LOGFILENAME->data );
105
106 //print VCS info
107 CHAR *VCSInfoString;
108 XLAL_CHECK( ( VCSInfoString = XLALVCSInfoString( lalPulsarVCSInfoList, 0, "%% " ) ) != NULL, XLAL_EFUNC );
109 fprintf( LOG, "%s\n", VCSInfoString );
110 fprintf( stderr, "%s\n", VCSInfoString );
111 XLALFree( VCSInfoString );
112
113 //print start time
114 fprintf( stderr, "Program executed on %s", asctime( ptm ) );
115 fprintf( LOG, "Program executed on %s", asctime( ptm ) );
116
117 //Print the output directory
118 fprintf( stderr, "Output directory: %s\n", uvar.outdirectory );
119 fprintf( LOG, "Output directory: %s\n", uvar.outdirectory );
120
121 //Set up the MultiLALDetector structure
123 XLAL_CHECK( ( detectors = setupMultiLALDetector( uvar.IFO ) ) != NULL, XLAL_EFUNC );
124
125 //Allocate GSL RNG
126 gsl_rng *rng = NULL;
127 XLAL_CHECK( ( rng = gsl_rng_alloc( gsl_rng_mt19937 ) ) != NULL, XLAL_EFUNC );
128 //Random seed value settings for random number generator
129 //If the chooseSeed option was given, then:
130 //seed = (IFO multiplier)*fabs(round(fmin + fspan + Pmin + Pmax + dfmin + dfmax + alpha + delta))
131 if ( uvar.chooseSeed ) {
132 UINT8 IFOmultiplier = 0;
133 if ( strcmp( uvar.IFO->data[0], "H1" ) == 0 ) {
134 IFOmultiplier = 1; //H1 gets multiplier 1
135 } else if ( strcmp( uvar.IFO->data[0], "L1" ) == 0 ) {
136 IFOmultiplier = 2; //L1 gets multiplier 2
137 } else {
138 IFOmultiplier = 3; //V1 gets multiplier 3
139 }
140 uvar.randSeed = IFOmultiplier * ( UINT8 )fabs( round( uvar.fmin + uvar.fspan + uvar.Pmin + uvar.Pmax + uvar.dfmin + uvar.dfmax ) );
141 }
142 gsl_rng_set( rng, uvar.randSeed ); //Set the random number generator with the given seed
143
144 //Print to stderr the parameters of the search
145 fprintf( stderr, "Tobs = %f sec\n", uvar.Tobs );
146 fprintf( stderr, "Tsft = %f sec\n", uvar.Tsft );
147 fprintf( stderr, "SFToverlap = %f sec\n", uvar.SFToverlap );
148 fprintf( stderr, "fmin = %f Hz\n", uvar.fmin );
149 fprintf( stderr, "fspan = %f Hz\n", uvar.fspan );
150 fprintf( stderr, "Pmin = %f s\n", uvar.Pmin );
151 fprintf( stderr, "Pmax = %f s\n", uvar.Pmax );
152 fprintf( stderr, "dfmin = %f Hz\n", uvar.dfmin );
153 fprintf( stderr, "dfmax = %f Hz\n", uvar.dfmax );
154 fprintf( stderr, "Running median blocksize = %d\n", uvar.blksize );
155 fprintf( stderr, "FFT plan flag = %d\n", uvar.FFTplanFlag );
156 fprintf( stderr, "RNG seed: %d\n", uvar.randSeed );
157 fprintf( LOG, "FAR for templates = %g\n", uvar.tmplfar );
158 fprintf( stderr, "FAR for templates = %g\n", uvar.tmplfar );
159
160 //Initialize ephemeris data structure
161 EphemerisData *edat = NULL;
162 XLAL_CHECK( ( edat = XLALInitBarycenter( uvar.ephemEarth, uvar.ephemSun ) ) != NULL, XLAL_EFUNC );
163
164 //Initialize timestamps and state series and set refTime
165 LIGOTimeGPS tStart;
167 MultiDetectorStateSeries *multiStateSeries = NULL;
168 XLALGPSSetREAL8( &tStart, uvar.t0 );
169 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC, "XLALGPSSetREAL8 failed\n" );
170 XLAL_CHECK( ( multiTimestamps = XLALMakeMultiTimestamps( tStart, uvar.Tobs, uvar.Tsft, uvar.SFToverlap, detectors->length ) ) != NULL, XLAL_EFUNC );
171 XLAL_CHECK( ( multiStateSeries = XLALGetMultiDetectorStates( multiTimestamps, detectors, edat, uvar.SFToverlap ) ) != NULL, XLAL_EFUNC );
172 LIGOTimeGPS refTime = multiTimestamps->data[0]->data[0];
173
174 //Maximum detector velocity in units of c from start of observation time - Tsft to end of observation + Tsft
175 REAL4 Vmax = 0.0;
176 for ( ii = 0; ii < ( INT4 )detectors->length; ii++ ) {
177 REAL4 detectorVmax = CompDetectorVmax2( multiTimestamps->data[ii], detectors->sites[ii], edat );
178 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC, "CompDetectorVmax2() failed\n" );
179 if ( detectorVmax > Vmax ) {
180 Vmax = detectorVmax;
181 }
182 }
183
184 //Assume maximum possible bin shift
185 INT4 maxbinshift = ( INT4 )round( Vmax * ( uvar.fmin + uvar.fspan ) * uvar.Tsft ) + 1;
186 if ( uvar.dopplerMultiplier > 1.0 ) {
187 maxbinshift = ( INT4 )round( Vmax * ( uvar.fmin + uvar.fspan ) * uvar.Tsft * uvar.dopplerMultiplier ) + 1;
188 }
189 if ( detectors->length > 1 ) {
190 maxbinshift += 10; //Add 10 bins for when there is more than 1 detector
191 }
192
193 //Parameters for the sky-grid from a point/polygon or a sky-grid file
196 PulsarDopplerParams dopplerpos;
197 if ( XLALUserVarWasSet( &uvar.skyRegion ) ) {
198 scanInit.gridType = GRID_ISOTROPIC; //Default value for an approximate-isotropic grid
199 scanInit.skyRegionString = uvar.skyRegion; //"allsky" = Default value for all-sky search
200 scanInit.numSkyPartitions = 1; //Default value so sky is not broken into chunks
201 scanInit.Freq = uvar.fmin + 0.5 * uvar.fspan; //Midpoint of the frequency band
202 scanInit.dAlpha = 0.5 / ( ( uvar.fmin + uvar.fspan ) * uvar.Tsft * Vmax );
203 scanInit.dDelta = scanInit.dAlpha;
204 fprintf( LOG, "Sky region = %s\n", uvar.skyRegion );
205 fprintf( stderr, "Sky region = %s\n", uvar.skyRegion );
206 } else {
207 scanInit.gridType = GRID_FILE_SKYGRID;
208 scanInit.skyGridFile = uvar.skyRegionFile;
209 scanInit.numSkyPartitions = 1; //Default value so sky is not broken into chunks
210 scanInit.Freq = uvar.fmin + 0.5 * uvar.fspan; //Midpoint of the frequency band
211 fprintf( LOG, "Sky file = %s\n", uvar.skyRegionFile );
212 fprintf( stderr, "Sky file = %s\n", uvar.skyRegionFile );
213 }
214
215 //Initialize the sky-grid
216 InitDopplerSkyScan( &status, &scan, &scanInit );
217 XLAL_CHECK( status.statusCode == 0, XLAL_EFUNC );
218
219 //Start at first location
220 XLAL_CHECK( XLALNextDopplerSkyPos( &dopplerpos, &scan ) == XLAL_SUCCESS, XLAL_EFUNC );
221
222 //Basic units
223 REAL8 tempfspan = uvar.fspan + 2.0 * uvar.dfmax + ( uvar.blksize - 1 + 12 ) / uvar.Tsft; //= fspan + 2*dfmax + running median blocksize-1 + extra bins (Hz)
224 INT4 tempnumfbins = ( INT4 )round( tempfspan * uvar.Tsft ) + 1; //= number of bins in tempfspan
225
226 //Determine band size to get the SFT data (remember to get extra bins because of the running median and the bin shifts due to detector velocity) with nudge of 0.1/Tsft for rounding issues
227 REAL8 minfbin = round( uvar.fmin * uvar.Tsft - uvar.dfmax * uvar.Tsft - 0.5 * ( uvar.blksize - 1 ) - ( REAL8 )( maxbinshift ) - 6.0 ) / uvar.Tsft + 0.1 / uvar.Tsft;
228 REAL8 maxfbin = round( ( uvar.fmin + uvar.fspan ) * uvar.Tsft + uvar.dfmax * uvar.Tsft + 0.5 * ( uvar.blksize - 1 ) + ( REAL8 )( maxbinshift ) + 6.0 ) / uvar.Tsft - 0.1 / uvar.Tsft;
229
230 //Allocate memory for ffdata structure
231 ffdataStruct *ffdata = NULL;
232 XLAL_CHECK( ( ffdata = createffdata( &uvar ) ) != NULL, XLAL_EFUNC );
233
234 //Maximum number of IHS values to sum = twice the maximum modulation depth
235 //Minimum number of IHS values to sum = twice the minimum modulation depth
236 INT4 maxrows = ( INT4 )round( 2.0 * uvar.dfmax * uvar.Tsft ) + 1;
237 fprintf( LOG, "Maximum row width to be searched = %d\n", maxrows );
238 fprintf( stderr, "Maximum row width to be searched = %d\n", maxrows );
239
240 //parse noise list
241 MultiNoiseFloor multiNoiseFloor;
242 XLAL_CHECK( XLALParseMultiNoiseFloor( &multiNoiseFloor, uvar.avesqrtSh, uvar.IFO->length ) == XLAL_SUCCESS, XLAL_EFUNC );
243
244 //TF normalization
245 ffdata->tfnormalization = 2.0 / uvar.Tsft / ( multiNoiseFloor.sqrtSn[0] * multiNoiseFloor.sqrtSn[0] );
246
247 //Read in the data from SFTs or generate SFTs
248 MultiSFTVector *multiSFTvector = NULL;
249 if ( !XLALUserVarWasSet( &uvar.injectionSources ) && XLALUserVarWasSet( &uvar.inputSFTs ) && !uvar.gaussNoiseWithSFTgaps && !XLALUserVarWasSet( &uvar.timestampsFile ) && !XLALUserVarWasSet( &uvar.segmentFile ) ) {
250 XLAL_CHECK( ( multiSFTvector = getMultiSFTVector( &uvar, minfbin, maxfbin ) ) != NULL, XLAL_EFUNC );
251 } else {
252 XLAL_CHECK( ( multiSFTvector = generateSFTdata( &uvar, detectors, edat, maxbinshift, rng ) ) != NULL, XLAL_EFUNC );
253 }
254
255 //Get timestamps of SFTs and determine joint timestamps
256 MultiLIGOTimeGPSVector *multiSFTtimestamps = NULL;
257 XLAL_CHECK( ( multiSFTtimestamps = XLALExtractMultiTimestampsFromSFTs( multiSFTvector ) ) != NULL, XLAL_EFUNC );
258 LIGOTimeGPSVector *jointSFTtimestamps = NULL;
259 XLAL_CHECK( ( jointSFTtimestamps = jointTimestampsFromMultiTimestamps( multiSFTtimestamps ) ) != NULL, XLAL_EFUNC );
260 XLALDestroyMultiTimestamps( multiSFTtimestamps );
261
262 if ( ( REAL8 )jointSFTtimestamps->length / ( REAL8 )ffdata->numffts < 0.1 ) {
263 fprintf( stderr, "%s: The useable SFTs cover less than 10 percent of the total observation time\n", __func__ );
264 fprintf( LOG, "%s: The useable SFTs cover less than 10 percent of the total observation time\n", __func__ );
265 //print end time
266 time( &programendtime );
267 ptm = localtime( &programendtime );
268 fprintf( stderr, "Program finished on %s", asctime( ptm ) );
269 fprintf( LOG, "Program finished on %s", asctime( ptm ) );
270 fclose( LOG );
271 return 0;
272 }
273
274 //Print SFT times, if requested by user
275 if ( uvar.printSFTtimes ) {
276 XLAL_CHECK( printSFTtimestamps2File( multiSFTvector, &uvar ) == XLAL_SUCCESS, XLAL_EFUNC );
277 }
278
279 //Next the running mean is calculated and we need each detector SFT data converted to powers for this step
280 //Calculate the running mean values of the SFTs (output here is smaller than the data read in for the SFTs). Here,
281 //numfbins needs to be the bins you expect to come out of the running means -- the band you are going
282 //to search!
283 REAL4VectorAligned *tfdata = NULL, *oneSFTbackgroundRatio = NULL;
284 REAL4VectorAlignedArray *backgrounds = NULL, *backgroundRatioMeans = NULL;
285 XLAL_CHECK( ( oneSFTbackgroundRatio = XLALCreateREAL4VectorAligned( ffdata->numfbins + 2 * maxbinshift, 32 ) ) != NULL, XLAL_EFUNC );
286 XLAL_CHECK( ( backgrounds = createREAL4VectorAlignedArray( detectors->length, ffdata->numffts * ( ffdata->numfbins + 2 * maxbinshift ), 32 ) ) != NULL, XLAL_EFUNC );
287 XLAL_CHECK( ( backgroundRatioMeans = createREAL4VectorAlignedArray( detectors->length, ffdata->numffts, 32 ) ) != NULL, XLAL_EFUNC ); //mean value of background ratio, 1 per SFT
288 for ( ii = 0; ii < ( INT4 )detectors->length; ii++ ) {
289 memset( backgrounds->data[ii]->data, 0, sizeof( REAL4 )*backgrounds->data[ii]->length );
290 memset( backgroundRatioMeans->data[ii]->data, 0, sizeof( REAL4 )*backgroundRatioMeans->data[ii]->length );
291
292 //Convert SFT data to powers
293 REAL4VectorAligned *tmpTFdata = NULL;
294 XLAL_CHECK( ( tmpTFdata = convertSFTdataToPowers( multiSFTvector->data[ii], &uvar, 2.0 / ( multiNoiseFloor.sqrtSn[0] * multiNoiseFloor.sqrtSn[0] ) / uvar.Tsft ) ) != NULL, XLAL_EFUNC );
295 //Determine running means
296 if ( !uvar.signalOnly ) {
297 XLAL_CHECK( tfRngMeans( backgrounds->data[ii], tmpTFdata, ffdata->numffts, ffdata->numfbins + 2 * maxbinshift, uvar.blksize ) == XLAL_SUCCESS, XLAL_EFUNC );
298 } else {
299 memset( backgrounds->data[ii]->data, 0, sizeof( REAL4 )*backgrounds->data[ii]->length );
300 }
301 if ( detectors->length == 1 ) {
302 XLAL_CHECK( ( tfdata = XLALCreateREAL4VectorAligned( tmpTFdata->length, 32 ) ) != NULL, XLAL_EFUNC );
303 memcpy( tfdata->data, tmpTFdata->data, sizeof( REAL4 )*tmpTFdata->length );
304 XLALDestroyMultiSFTVector( multiSFTvector );
305 }
307 }
308 //Replace any zeros in detector 0 with data from first detector available
310 //copy background data from detector 0
311 REAL4VectorAligned *background = NULL, *background0 = NULL;
312 XLAL_CHECK( ( background = XLALCreateREAL4VectorAligned( ffdata->numffts * ( ffdata->numfbins + 2 * maxbinshift ), 32 ) ) != NULL, XLAL_EFUNC );
313 XLAL_CHECK( ( background0 = XLALCreateREAL4VectorAligned( ffdata->numffts * ( ffdata->numfbins + 2 * maxbinshift ), 32 ) ) != NULL, XLAL_EFUNC );
314 if ( !uvar.signalOnly ) {
315 memcpy( background->data, backgrounds->data[0]->data, sizeof( REAL4 )*backgrounds->data[0]->length );
316 memcpy( background0->data, backgrounds->data[0]->data, sizeof( REAL4 )*backgrounds->data[0]->length );
317 } else {
318 memset( background->data, 0, sizeof( REAL4 )*background->length );
319 memset( background0->data, 0, sizeof( REAL4 )*background->length );
320 }
321 //compute ratios, this only gets executed if detectors->length>1
322 for ( ii = 1; ii < ( INT4 )detectors->length; ii++ ) {
323 for ( jj = 0; jj < ffdata->numffts; jj++ ) {
324 if ( backgrounds->data[ii]->data[jj * ( ffdata->numfbins + 2 * maxbinshift )] != 0.0 ) {
325 memcpy( oneSFTbackgroundRatio->data, &( backgrounds->data[0]->data[jj * ( ffdata->numfbins + 2 * maxbinshift )] ), sizeof( REAL4 )*oneSFTbackgroundRatio->length );
326 for ( UINT4 kk = 0; kk < oneSFTbackgroundRatio->length; kk++ ) {
327 oneSFTbackgroundRatio->data[kk] /= backgrounds->data[ii]->data[jj * ( ffdata->numfbins + 2 * maxbinshift ) + kk];
328 }
329 backgroundRatioMeans->data[ii]->data[jj] = calcMean( oneSFTbackgroundRatio );
331 } else {
332 backgroundRatioMeans->data[ii]->data[jj] = 0.0;
333 }
334 }
335 }
336 XLALDestroyREAL4VectorAligned( oneSFTbackgroundRatio );
337 destroyREAL4VectorAlignedArray( backgrounds );
338 //Need to do a little data manipulation here for the extra background data
339 if ( detectors->length > 1 ) {
340 XLALDestroyREAL4VectorAligned( background );
341 XLAL_CHECK( ( background = XLALCreateREAL4VectorAligned( ffdata->numffts * ( ffdata->numfbins + 2 * ( maxbinshift - 10 ) ), 32 ) ) != NULL, XLAL_EFUNC );
342 for ( ii = 0; ii < ffdata->numffts; ii++ ) {
343 if ( background0->data[( ffdata->numfbins + 2 * maxbinshift ) * ii] != 0.0 ) {
344 memcpy( &( background->data[( ffdata->numfbins + 2 * ( maxbinshift - 10 ) )*ii] ), &( background0->data[( ffdata->numfbins + 2 * maxbinshift )*ii + 10] ), sizeof( REAL4 ) * ( ffdata->numfbins + 2 * ( maxbinshift - 10 ) ) );
345 } else {
346 memset( &( background->data[( ffdata->numfbins + 2 * ( maxbinshift - 10 ) )*ii] ), 0, sizeof( REAL4 ) * ( ffdata->numfbins + 2 * ( maxbinshift - 10 ) ) );
347 }
348 }
349 XLALDestroyREAL4VectorAligned( background0 );
350 XLAL_CHECK( ( background0 = XLALCreateREAL4VectorAligned( background->length, 32 ) ) != NULL, XLAL_EFUNC );
351 memcpy( background0->data, background->data, sizeof( REAL4 )*background->length );
352 }
353
354 if ( detectors->length > 1 ) {
355 maxbinshift -= 10;
356 }
357
358 //Background scaling and normalization
359 REAL4VectorAligned *backgroundScaling = NULL;
360 XLAL_CHECK( ( backgroundScaling = XLALCreateREAL4VectorAligned( background->length, 32 ) ) != NULL, XLAL_EFUNC );
361 memset( backgroundScaling->data, 0, sizeof( REAL4 )*backgroundScaling->length );
362 for ( ii = 0; ii < ( INT4 )background->length; ii++ ) if ( background->data[ii] != 0.0 ) {
363 backgroundScaling->data[ii] = 1.0;
364 }
365
366 //Line detection only if there is valid noise to be looking at
367 INT4Vector *lines = NULL;
368 INT4 heavilyContaminatedBand = 0;
369 if ( XLALUserVarWasSet( &uvar.lineDetection ) && !uvar.signalOnly && detectors->length == 1 ) {
370 lines = detectLines_simple( tfdata, ffdata, &uvar );
372 if ( lines != NULL ) {
373 if ( ( REAL4 )lines->length / ( ffdata->numfbins + 2 * maxbinshift ) >= 0.1 ) {
374 heavilyContaminatedBand = 1;
375 fprintf( LOG, "WARNING: Band is heavily contaminated by artifacts.\n" );
376 fprintf( stderr, "WARNING: Band is heavily contaminated by artifacts.\n" );
377 }
378 }
379 }
380 //If the band is heavily contaminated by lines, don't do any follow up.
381 if ( heavilyContaminatedBand ) {
382 uvar.IHSonly = 1;
383 }
384
385 //TEST: Try cleaning lines
386 /* XLAL_CHECK( cleanLines(tfdata, background, lines, &uvar) == XLAL_SUCCESS, XLAL_EFUNC ); */
387 /* if (lines!=NULL) { */
388 /* if (!uvar.signalOnly) XLAL_CHECK( tfRngMeans(background, tfdata, ffdata->numffts, ffdata->numfbins + 2*maxbinshift, uvar.blksize) == XLAL_SUCCESS, XLAL_EFUNC ); */
389 /* else memset(background->data, 0, sizeof(REAL4)*background->length); */
390 /* } */
391 /* XLALDestroyINT4Vector(lines); */
392 /* lines = NULL; */
393
394 //Existing SFTs listed in this vector
395 INT4Vector *sftexist = NULL;
396 INT4Vector *indexValuesOfExistingSFTs = NULL;
397 REAL4 frac_tobs_complete = 1.0;
398 if ( !uvar.signalOnly && detectors->length == 1 ) {
399 XLAL_CHECK( ( sftexist = existingSFTs( tfdata, ( UINT4 )ffdata->numffts ) ) != NULL, XLAL_EFUNC );
400 INT4 totalincludedsftnumber = 0;
401 for ( ii = 0; ii < ( INT4 )sftexist->length; ii++ ) if ( sftexist->data[ii] == 1 ) {
402 totalincludedsftnumber++;
403 }
404 frac_tobs_complete = ( REAL4 )totalincludedsftnumber / ( REAL4 )sftexist->length;
405 fprintf( LOG, "Duty factor of usable SFTs = %f\n", frac_tobs_complete );
406 fprintf( stderr, "Duty factor of usable SFTs = %f\n", frac_tobs_complete );
407 if ( frac_tobs_complete < 0.1 ) {
408 fprintf( stderr, "%s: The useable SFTs cover less than 10 percent of the total observation time\n", __func__ );
409 fprintf( LOG, "%s: The useable SFTs cover less than 10 percent of the total observation time\n", __func__ );
410 //print end time
411 time( &programendtime );
412 ptm = localtime( &programendtime );
413 fprintf( stderr, "Program finished on %s", asctime( ptm ) );
414 fprintf( LOG, "Program finished on %s", asctime( ptm ) );
415 fclose( LOG );
416 return 0;
417 }
418
419 //Index values of existing SFTs
420 XLAL_CHECK( ( indexValuesOfExistingSFTs = XLALCreateINT4Vector( totalincludedsftnumber ) ) != NULL, XLAL_EFUNC );
421 jj = 0;
422 for ( ii = 0; ii < ( INT4 )sftexist->length; ii++ ) {
423 if ( sftexist->data[ii] == 1 ) {
424 indexValuesOfExistingSFTs->data[jj] = ii;
425 jj++;
426 }
427 }
428 }
429
430 //I wrote this to compensate for a bad input of the expected noise floor
431 REAL8 backgroundmeannormfactor = 1.0;
432 if ( !uvar.signalOnly && detectors->length == 1 ) {
433 REAL4 harmonicMean = 0.0;
434 XLAL_CHECK( calcHarmonicMean( &harmonicMean, background, ffdata->numfbins + 2 * maxbinshift, ffdata->numffts ) == XLAL_SUCCESS, XLAL_EFUNC );
435 backgroundmeannormfactor = 1.0 / harmonicMean;
436 ffdata->tfnormalization *= backgroundmeannormfactor;
437 }
438
439 //Need to reduce the original TF data to remove the excess bins used for running median calculation. Normalize the TF as the same as the background was normalized
440 REAL4VectorAligned *usableTFdata = NULL;
441 if ( detectors->length == 1 ) {
442 XLAL_CHECK( ( usableTFdata = XLALCreateREAL4VectorAligned( background->length, 32 ) ) != NULL, XLAL_EFUNC );
443 for ( ii = 0; ii < ffdata->numffts; ii++ ) {
444 memcpy( &( usableTFdata->data[ii * ( ffdata->numfbins + 2 * maxbinshift )] ), &( tfdata->data[ii * ( tempnumfbins + 2 * maxbinshift ) + ( INT4 )round( 0.5 * ( uvar.blksize - 1 ) )] ), sizeof( REAL4 ) * ( ffdata->numfbins + 2 * maxbinshift ) );
445 }
446 XLAL_CHECK( XLALVectorScaleREAL4( usableTFdata->data, backgroundmeannormfactor, usableTFdata->data, usableTFdata->length ) == XLAL_SUCCESS, XLAL_EFUNC );
447 XLAL_CHECK( XLALVectorScaleREAL4( background->data, backgroundmeannormfactor, background->data, background->length ) == XLAL_SUCCESS, XLAL_EFUNC );
448 //At this point the TF plane and the running median calculation are the same size=numffts*(numfbins + 2*maxbinshift)
449 //We can delete the originally loaded SFTs since we have the usableTFdata saved
451 }
452
453 //Print out data product if requested
454 if ( uvar.printData && detectors->length == 1 ) {
455 XLAL_CHECK( printREAL4Vector2File( ( REAL4Vector * )usableTFdata, uvar.outdirectory, "tfdata.dat" ) == XLAL_SUCCESS, XLAL_EFUNC );
456 }
457
458 //Do mean subtraction of TFdata here--modifies the usableTFdata vector!!!
459 if ( detectors->length == 1 ) {
460 XLAL_CHECK( tfMeanSubtract( usableTFdata, background, backgroundScaling, ffdata->numffts, ffdata->numfbins + 2 * maxbinshift ) == XLAL_SUCCESS, XLAL_EFUNC );
461 }
462
463 //Initialize IHS structures
464 ihsMaximaStruct *ihsmaxima = NULL;
465 ihsfarStruct *ihsfarstruct = NULL;
466 XLAL_CHECK( ( ihsmaxima = createihsMaxima( ffdata->numfbins, maxrows ) ) != NULL, XLAL_EFUNC );
467 XLAL_CHECK( ( ihsfarstruct = createihsfarStruct( maxrows, &uvar ) ) != NULL, XLAL_EFUNC );
468
469 //If necessary, read in template bank file
470 TwoSpectTemplateVector *templateVec = NULL;
471 if ( XLALUserVarWasSet( &( uvar.templatebankfile ) ) ) {
472 XLAL_CHECK( ( templateVec = readTwoSpectTemplateVector( uvar.templatebankfile ) ) != NULL, XLAL_EFUNC );
473 XLAL_CHECK( templateVec->Tsft == uvar.Tsft && templateVec->SFToverlap == uvar.SFToverlap && templateVec->Tobs == uvar.Tobs, XLAL_EINVAL, "Template bank %s doesn't match input parameters\n", uvar.templatebankfile );
474 }
475
476 //Initialize candidate vectors and upper limits
477 candidateVector *gaussCandidates1 = NULL, *gaussCandidates2 = NULL, *gaussCandidates3 = NULL, *gaussCandidates4 = NULL, *exactCandidates1 = NULL, *exactCandidates2 = NULL, *ihsCandidates = NULL;
478 UpperLimitVector *upperlimits = NULL;
479 XLAL_CHECK( ( gaussCandidates1 = createcandidateVector( 100 ) ) != NULL, XLAL_EFUNC, "createCandidateVector(%d) failed\n", 100 );
480 XLAL_CHECK( ( gaussCandidates2 = createcandidateVector( 100 ) ) != NULL, XLAL_EFUNC, "createCandidateVector(%d) failed\n", 100 );
481 XLAL_CHECK( ( gaussCandidates3 = createcandidateVector( 100 ) ) != NULL, XLAL_EFUNC, "createCandidateVector(%d) failed\n", 100 );
482 XLAL_CHECK( ( gaussCandidates4 = createcandidateVector( 100 ) ) != NULL, XLAL_EFUNC, "createCandidateVector(%d) failed\n", 100 );
483 XLAL_CHECK( ( exactCandidates1 = createcandidateVector( 100 ) ) != NULL, XLAL_EFUNC, "createCandidateVector(%d) failed\n", 100 );
484 XLAL_CHECK( ( exactCandidates2 = createcandidateVector( 100 ) ) != NULL, XLAL_EFUNC, "createCandidateVector(%d) failed\n", 100 );
485 XLAL_CHECK( ( ihsCandidates = createcandidateVector( 100 ) ) != NULL, XLAL_EFUNC, "createCandidateVector(%d) failed\n", 100 );
486 XLAL_CHECK( ( upperlimits = createUpperLimitVector( 1 ) ) != NULL, XLAL_EFUNC, "createUpperLimitVector(%d) failed.\n", 1 );
487
488 //Initialize second FFT plan
489 REAL4FFTPlan *secondFFTplan = NULL;
490 XLAL_CHECK( ( secondFFTplan = XLALCreateForwardREAL4FFTPlan( ffdata->numffts, uvar.FFTplanFlag ) ) != NULL, XLAL_EFUNC );
491
492 //Initialize aveNoise and binshifts vectors
493 REAL4VectorAligned *aveNoise = NULL;
494 INT4Vector *binshifts = NULL;
495 XLAL_CHECK( ( aveNoise = XLALCreateREAL4VectorAligned( ffdata->numfprbins, 32 ) ) != NULL, XLAL_EFUNC );
496 XLAL_CHECK( ( binshifts = XLALCreateINT4Vector( ffdata->numffts ) ) != NULL, XLAL_EFUNC );
497
498 //Create exponentially distributed random values for background noise estimate to sample from
499 REAL4VectorAligned *expRandVals = NULL;
500 XLAL_CHECK( ( expRandVals = expRandNumVector( 100 * ffdata->numffts, 1.0, rng ) ) != NULL, XLAL_EFUNC );
501
502 //Allocate aveNoiseInTime vector
503 REAL4VectorAligned *aveNoiseInTime = NULL;
504 XLAL_CHECK( ( aveNoiseInTime = XLALCreateREAL4VectorAligned( ffdata->numffts, 32 ) ) != NULL, XLAL_EFUNC );
505
506 //Initialize to -1.0 for far just at the start
507 ihsfarstruct->ihsfar->data[0] = -1.0;
508 REAL4 antweightsrms = 0.0;
509
510 //Set TF normalization
511 if ( detectors->length == 1 ) {
512 ffdata->tfnormalization *= 0.5 * uvar.Tsft;
513 }
514
515 //Antenna normalization (determined from injections on H1 at ra=0, dec=0, with circular polarization)
516 //When doing linear polarizations, the IHS factor needs to be 25.2*1.082 and this antenna weights
517 //function needs to be set to use linear polarization.
518 SkyPosition skypos0 = {0.0, 0.0, COORDINATESYSTEM_EQUATORIAL};
519 REAL4VectorAligned *antweightsforihs2h0 = NULL;
520 XLAL_CHECK( ( antweightsforihs2h0 = XLALCreateREAL4VectorAligned( ffdata->numffts, 32 ) ) != NULL, XLAL_EFUNC );
521 if ( XLALUserVarWasSet( &uvar.assumeNScosi ) && XLALUserVarWasSet( &uvar.assumeNSpsi ) ) {
522 XLAL_CHECK( CompAntennaPatternWeights2( antweightsforihs2h0, skypos0, multiTimestamps->data[0], lalCachedDetectors[LAL_LHO_4K_DETECTOR], &uvar.assumeNScosi, &uvar.assumeNSpsi ) == XLAL_SUCCESS, XLAL_EFUNC );
523 } else if ( XLALUserVarWasSet( &uvar.assumeNScosi ) ) {
524 XLAL_CHECK( CompAntennaPatternWeights2( antweightsforihs2h0, skypos0, multiTimestamps->data[0], lalCachedDetectors[LAL_LHO_4K_DETECTOR], &uvar.assumeNScosi, NULL ) == XLAL_SUCCESS, XLAL_EFUNC );
525 } else {
526 REAL8 cosi = 1.0, psi = 0.0;
528 }
529 REAL4VectorAligned *backgroundForihs2h0 = NULL, *backgroundScalingForihs2h0 = NULL;
530 XLAL_CHECK( ( backgroundForihs2h0 = XLALCreateREAL4VectorAligned( background->length, 32 ) ) != NULL, XLAL_EFUNC );
531 XLAL_CHECK( ( backgroundScalingForihs2h0 = XLALCreateREAL4VectorAligned( backgroundScaling->length, 32 ) ) != NULL, XLAL_EFUNC );
532 memcpy( backgroundForihs2h0->data, background0->data, sizeof( REAL4 )*background->length );
533 memcpy( backgroundScalingForihs2h0->data, backgroundScaling->data, sizeof( REAL4 )*backgroundScaling->length );
534 INT4Vector *sftexistForihs2h0 = NULL;
535 if ( detectors->length > 1 ) {
536 MultiSSBtimes *multissb = NULL;
537 XLAL_CHECK( ( multissb = XLALGetMultiSSBtimes( multiStateSeries, skypos0, refTime, SSBPREC_RELATIVISTICOPT ) ) != NULL, XLAL_EFUNC );
538 MultiAMCoeffs *multiAMcoefficients = NULL;
539 XLAL_CHECK( ( multiAMcoefficients = XLALComputeMultiAMCoeffs( multiStateSeries, NULL, skypos0 ) ) != NULL, XLAL_EFUNC );
540 REAL4VectorAligned *tmpTFdata = NULL;
541 assumeNSparams XLAL_INIT_DECL( NSparams );
542 NSparams.assumeNSpos = skypos0;
543 if ( XLALUserVarWasSet( &uvar.assumeNScosi ) ) {
544 NSparams.assumeNScosi = &uvar.assumeNScosi;
545 }
546 if ( XLALUserVarWasSet( &uvar.assumeNSpsi ) ) {
547 NSparams.assumeNSpsi = &uvar.assumeNSpsi;
548 }
549 if ( XLALUserVarWasSet( &uvar.assumeNSGWfreq ) ) {
550 NSparams.assumeNSGWfreq = &uvar.assumeNSGWfreq;
551 NSparams.assumeNSorbitP = &uvar.assumeNSorbitP;
552 NSparams.assumeNSasini = &uvar.assumeNSasini;
553 NSparams.assumeNSorbitTp = &uvar.assumeNSorbitTp;
554 if ( XLALUserVarWasSet( &uvar.assumeNSrefTime ) ) {
555 NSparams.assumeNSrefTime = &uvar.assumeNSrefTime;
556 }
557 }
558 XLAL_CHECK( ( tmpTFdata = coherentlyAddSFTs( multiSFTvector, multissb, multiAMcoefficients, jointSFTtimestamps, backgroundRatioMeans, uvar.cosiSignCoherent, &NSparams, &uvar, backgroundScalingForihs2h0 ) ) != NULL, XLAL_EFUNC );
559 XLAL_CHECK( ( sftexistForihs2h0 = existingSFTs( tmpTFdata, ( UINT4 )ffdata->numffts ) ) != NULL, XLAL_EFUNC );
560 INT4 totalincludedsftnumber = 0;
561 for ( ii = 0; ii < ( INT4 )sftexistForihs2h0->length; ii++ ) if ( sftexistForihs2h0->data[ii] == 1 ) {
562 totalincludedsftnumber++;
563 }
564 frac_tobs_complete = ( REAL4 )totalincludedsftnumber / ( REAL4 )sftexistForihs2h0->length;
565 fprintf( LOG, "Duty factor of usable SFTs = %f\n", frac_tobs_complete );
566 fprintf( stderr, "Duty factor of usable SFTs = %f\n", frac_tobs_complete );
567 if ( frac_tobs_complete < 0.1 ) {
568 fprintf( stderr, "%s: The useable SFTs cover less than 10 percent of the total observation time\n", __func__ );
569 fprintf( LOG, "%s: The useable SFTs cover less than 10 percent of the total observation time\n", __func__ );
570 return 0;
571 }
572 XLAL_CHECK( checkBackgroundScaling( backgroundForihs2h0, backgroundScalingForihs2h0, sftexistForihs2h0 ) == XLAL_SUCCESS, XLAL_EFUNC );
574 XLALDestroyMultiAMCoeffs( multiAMcoefficients );
575 XLALDestroyMultiSSBtimes( multissb );
576 } else {
577 XLAL_CHECK( ( sftexistForihs2h0 = XLALCreateINT4Vector( sftexist->length ) ) != NULL, XLAL_EFUNC );
578 memcpy( sftexistForihs2h0->data, sftexist->data, sizeof( INT4 )*sftexist->length );
579 }
580
581 //Compute median values in time of the background
582 REAL4VectorAligned *aveNoiseInTimeForihs2h0 = NULL;
583 XLAL_CHECK( ( aveNoiseInTimeForihs2h0 = XLALCreateREAL4VectorAligned( aveNoiseInTime->length, 32 ) ) != NULL, XLAL_EFUNC );
584 XLAL_CHECK( medianBackgroundBandInTime( aveNoiseInTimeForihs2h0, backgroundForihs2h0, sftexistForihs2h0 ) == XLAL_SUCCESS, XLAL_EFUNC );
585
586 //Antenna normalization for different sky locations
587 REAL8 skypointffnormalization = 1.0;
588 XLAL_CHECK( ffPlaneNoise( aveNoise, &uvar, sftexistForihs2h0, aveNoiseInTimeForihs2h0, antweightsforihs2h0, backgroundScalingForihs2h0, secondFFTplan, expRandVals, rng, &( skypointffnormalization ) ) == XLAL_SUCCESS, XLAL_EFUNC );
589
590 //Print message that we start the analysis
591 fprintf( LOG, "Starting TwoSpect analysis...\n" );
592 fprintf( stderr, "Starting TwoSpect analysis...\n" );
593
594 //Search over the sky region (outer loop of single TwoSpect program instance)
595 while ( scan.state != STATE_FINISHED ) {
596 fprintf( LOG, "Sky location: RA = %g, DEC = %g\n", dopplerpos.Alpha, dopplerpos.Delta );
597 fprintf( stderr, "Sky location: RA = %g, DEC = %g\n", dopplerpos.Alpha, dopplerpos.Delta );
598
599 SkyPosition skypos = {dopplerpos.Alpha, dopplerpos.Delta, COORDINATESYSTEM_EQUATORIAL};
600
601 //Compute antenna pattern weights. If antennaOff input flag is given, then set all values equal to 1.0
602 REAL4VectorAligned *antweights = NULL;
603 XLAL_CHECK( ( antweights = XLALCreateREAL4VectorAligned( ffdata->numffts, 32 ) ) != NULL, XLAL_EFUNC );
604 memset( antweights->data, 0, antweights->length * sizeof( REAL4 ) );
605 if ( uvar.antennaOff ) for ( ii = 0; ii < ( INT4 )antweights->length; ii++ ) {
606 antweights->data[ii] = 1.0;
607 } else {
608 if ( XLALUserVarWasSet( &uvar.assumeNScosi ) && XLALUserVarWasSet( &uvar.assumeNSpsi ) ) {
609 XLAL_CHECK( CompAntennaPatternWeights2( antweights, skypos, multiTimestamps->data[0], detectors->sites[0], &uvar.assumeNScosi, &uvar.assumeNSpsi ) == XLAL_SUCCESS, XLAL_EFUNC );
610 } else if ( XLALUserVarWasSet( &uvar.assumeNScosi ) ) {
611 XLAL_CHECK( CompAntennaPatternWeights2( antweights, skypos, multiTimestamps->data[0], detectors->sites[0], &uvar.assumeNScosi, NULL ) == XLAL_SUCCESS, XLAL_EFUNC );
612 } else {
613 REAL8 cosi = 1.0, psi = 0.0;
614 XLAL_CHECK( CompAntennaPatternWeights2( antweights, skypos, multiTimestamps->data[0], detectors->sites[0], &cosi, &psi ) == XLAL_SUCCESS, XLAL_EFUNC );
615 }
616 }
617
618 //Get SSB times
619 MultiSSBtimes *multissb = NULL;
620 XLAL_CHECK( ( multissb = XLALGetMultiSSBtimes( multiStateSeries, skypos, refTime, SSBPREC_RELATIVISTICOPT ) ) != NULL, XLAL_EFUNC );
621
622 //Compute the bin shifts for each SFT
623 XLAL_CHECK( CompBinShifts( binshifts, multissb->data[0], uvar.fmin + 0.5 * uvar.fspan, uvar.Tsft, uvar.dopplerMultiplier ) == XLAL_SUCCESS, XLAL_EFUNC );
624
625 //Coherent combination
626 if ( detectors->length >= 2 ) {
627 //TF normalization must be reset every time
628 ffdata->tfnormalization = 2.0 / uvar.Tsft / ( multiNoiseFloor.sqrtSn[0] * multiNoiseFloor.sqrtSn[0] );
629
630 //Get multiAMcoefficients
631 MultiAMCoeffs *multiAMcoefficients = NULL;
632 XLAL_CHECK( ( multiAMcoefficients = XLALComputeMultiAMCoeffs( multiStateSeries, NULL, skypos ) ) != NULL, XLAL_EFUNC );
633
634 //Coherenly combine the SFT data
635 assumeNSparams XLAL_INIT_DECL( NSparams );
636 NSparams.assumeNSpos = skypos;
637 if ( XLALUserVarWasSet( &uvar.assumeNScosi ) ) {
638 NSparams.assumeNScosi = &uvar.assumeNScosi;
639 }
640 if ( XLALUserVarWasSet( &uvar.assumeNSpsi ) ) {
641 NSparams.assumeNSpsi = &uvar.assumeNSpsi;
642 }
643 if ( XLALUserVarWasSet( &uvar.assumeNSGWfreq ) ) {
644 NSparams.assumeNSGWfreq = &uvar.assumeNSGWfreq;
645 NSparams.assumeNSorbitP = &uvar.assumeNSorbitP;
646 NSparams.assumeNSasini = &uvar.assumeNSasini;
647 NSparams.assumeNSorbitTp = &uvar.assumeNSorbitTp;
648 if ( XLALUserVarWasSet( &uvar.assumeNSrefTime ) ) {
649 NSparams.assumeNSrefTime = &uvar.assumeNSrefTime;
650 }
651 }
652 XLAL_CHECK( ( tfdata = coherentlyAddSFTs( multiSFTvector, multissb, multiAMcoefficients, jointSFTtimestamps, backgroundRatioMeans, uvar.cosiSignCoherent, &NSparams, &uvar, backgroundScaling ) ) != NULL, XLAL_EFUNC );
653
654 XLALDestroyMultiAMCoeffs( multiAMcoefficients );
655
656 //Continue like there is just one detector
657
658 //Check for lines
659 if ( XLALUserVarWasSet( &uvar.lineDetection ) && !uvar.signalOnly ) {
660 lines = detectLines_simple( tfdata, ffdata, &uvar );
662 if ( lines != NULL ) {
663 if ( ( REAL4 )lines->length / ( ffdata->numfbins + 2 * maxbinshift ) >= 0.1 ) {
664 heavilyContaminatedBand = 1;
665 fprintf( LOG, "WARNING: Band is heavily contaminated by artifacts.\n" );
666 fprintf( stderr, "WARNING: Band is heavily contaminated by artifacts.\n" );
667 }
668 }
669 }
670 if ( heavilyContaminatedBand ) {
671 uvar.IHSonly = 1;
672 }
673
674 //Copy background from storage
675 memcpy( background->data, background0->data, sizeof( REAL4 )*background->length );
676
677 if ( !uvar.signalOnly ) {
678 XLAL_CHECK( ( sftexist = existingSFTs( tfdata, ( UINT4 )ffdata->numffts ) ) != NULL, XLAL_EFUNC );
679 XLAL_CHECK( checkBackgroundScaling( background, backgroundScaling, sftexist ) == XLAL_SUCCESS, XLAL_EFUNC );
680 INT4 totalincludedsftnumber = 0;
681 for ( ii = 0; ii < ( INT4 )sftexist->length; ii++ ) if ( sftexist->data[ii] == 1 ) {
682 totalincludedsftnumber++;
683 }
684 frac_tobs_complete = ( REAL4 )totalincludedsftnumber / ( REAL4 )sftexist->length;
685 fprintf( LOG, "Duty factor of usable SFTs = %f\n", frac_tobs_complete );
686 fprintf( stderr, "Duty factor of usable SFTs = %f\n", frac_tobs_complete );
687 if ( frac_tobs_complete < 0.1 ) {
688 fprintf( stderr, "%s: The useable SFTs cover less than 10 percent of the total observation time\n", __func__ );
689 fprintf( LOG, "%s: The useable SFTs cover less than 10 percent of the total observation time\n", __func__ );
690 return 0;
691 }
692
693 //Index values of existing SFTs
694 XLAL_CHECK( ( indexValuesOfExistingSFTs = XLALCreateINT4Vector( totalincludedsftnumber ) ) != NULL, XLAL_EFUNC );
695 jj = 0;
696 for ( ii = 0; ii < ( INT4 )sftexist->length; ii++ ) {
697 if ( sftexist->data[ii] == 1 ) {
698 indexValuesOfExistingSFTs->data[jj] = ii;
699 jj++;
700 }
701 }
702
703 REAL4 harmonicMean = 0.0;
704 XLAL_CHECK( calcHarmonicMean( &harmonicMean, background, ffdata->numfbins + 2 * maxbinshift, ffdata->numffts ) == XLAL_SUCCESS, XLAL_EFUNC );
705 backgroundmeannormfactor = 1.0 / harmonicMean;
706 ffdata->tfnormalization *= backgroundmeannormfactor;
707 }
708
709 XLAL_CHECK( ( usableTFdata = XLALCreateREAL4VectorAligned( background->length, 32 ) ) != NULL, XLAL_EFUNC );
710 for ( ii = 0; ii < ffdata->numffts; ii++ ) {
711 memcpy( &( usableTFdata->data[ii * ( ffdata->numfbins + 2 * maxbinshift )] ), &( tfdata->data[ii * ( tempnumfbins + 2 * maxbinshift ) + ( INT4 )round( 0.5 * ( uvar.blksize - 1 ) )] ), sizeof( REAL4 ) * ( ffdata->numfbins + 2 * maxbinshift ) );
712 }
713 XLAL_CHECK( XLALVectorScaleREAL4( usableTFdata->data, backgroundmeannormfactor, usableTFdata->data, usableTFdata->length ) == XLAL_SUCCESS, XLAL_EFUNC );
714 XLAL_CHECK( XLALVectorScaleREAL4( background->data, backgroundmeannormfactor, background->data, background->length ) == XLAL_SUCCESS, XLAL_EFUNC );
716
717 //Print out data product if requested
718 if ( uvar.printData ) {
719 XLAL_CHECK( printREAL4Vector2File( ( REAL4Vector * )usableTFdata, uvar.outdirectory, "tfdata.dat" ) == XLAL_SUCCESS, XLAL_EFUNC );
720 }
721
722 //Mean subtraction
723 XLAL_CHECK( tfMeanSubtract( usableTFdata, background, backgroundScaling, ffdata->numffts, ffdata->numfbins + 2 * maxbinshift ) == XLAL_SUCCESS, XLAL_EFUNC );
724
725 ffdata->tfnormalization *= 0.5 * uvar.Tsft;
726 }
727 /////
728
729 XLALDestroyMultiSSBtimes( multissb );
730
731 //Track identified lines
732 REAL4 fbin0 = ( REAL4 )( round( uvar.fmin * uvar.Tsft - uvar.dfmax * uvar.Tsft - 6.0 - 0.5 * ( uvar.blksize - 1 ) - ( REAL8 )( maxbinshift ) ) / uvar.Tsft );
733 REAL4 df = 1.0 / uvar.Tsft;
734 REAL4VectorSequence *trackedlines = NULL;
735 if ( lines != NULL ) {
736 XLAL_CHECK( ( trackedlines = trackLines( lines, binshifts, fbin0, df ) ) != NULL, XLAL_EFUNC );
737 XLALDestroyINT4Vector( lines );
738 }
739
740 //Calculate antenna RMS value
741 REAL4 currentAntWeightsRMS = 0.0;
742 XLAL_CHECK( calcRms( &currentAntWeightsRMS, antweights ) == XLAL_SUCCESS, XLAL_EFUNC );
743
744 //Slide SFTs here -- need to slide the data and the estimated background
745 REAL4VectorAligned *TFdata_slided = NULL, *background_slided = NULL, *backgroundScaling_slided = NULL;
746 XLAL_CHECK( ( TFdata_slided = XLALCreateREAL4VectorAligned( ffdata->numffts * ffdata->numfbins, 32 ) ) != NULL, XLAL_EFUNC );
747 XLAL_CHECK( ( background_slided = XLALCreateREAL4VectorAligned( TFdata_slided->length, 32 ) ) != NULL, XLAL_EFUNC );
748 XLAL_CHECK( ( backgroundScaling_slided = XLALCreateREAL4VectorAligned( TFdata_slided->length, 32 ) ) != NULL, XLAL_EFUNC );
749 XLAL_CHECK( slideTFdata( TFdata_slided, &uvar, usableTFdata, binshifts ) == XLAL_SUCCESS, XLAL_EFUNC );
750 XLAL_CHECK( slideTFdata( background_slided, &uvar, background, binshifts ) == XLAL_SUCCESS, XLAL_EFUNC );
751 XLAL_CHECK( slideTFdata( backgroundScaling_slided, &uvar, backgroundScaling, binshifts ) == XLAL_SUCCESS, XLAL_EFUNC );
752
753 if ( detectors->length > 1 ) {
754 XLALDestroyREAL4VectorAligned( usableTFdata );
755 }
756
757 //Print out data product if requested
758 if ( uvar.printData ) {
759 XLAL_CHECK( printREAL4Vector2File( ( REAL4Vector * )background_slided, uvar.outdirectory, "tfbackground.dat" ) == XLAL_SUCCESS, XLAL_EFUNC );
760 XLAL_CHECK( printREAL4Vector2File( ( REAL4Vector * )TFdata_slided, uvar.outdirectory, "tfdata_slided.dat" ) == XLAL_SUCCESS, XLAL_EFUNC );
761 }
762
763 //Check the RMS of the antenna weights; if a deviation of more than 1%, then reset the IHS FAR
764 if ( antweightsrms == 0.0 ) {
765 antweightsrms = currentAntWeightsRMS;
766 }
767 if ( fabs( currentAntWeightsRMS - antweightsrms ) / antweightsrms >= 0.01 ) {
768 ihsfarstruct->ihsfar->data[0] = -1.0;
769 antweightsrms = currentAntWeightsRMS;
770 }
771
772 //Compute median values
773 XLAL_CHECK( medianBackgroundBandInTime( aveNoiseInTime, background_slided, sftexist ) == XLAL_SUCCESS, XLAL_EFUNC );
774
775 //Average noise floor of FF plane for each 1st FFT frequency bin
776 ffdata->ffnormalization = 1.0;
777 XLAL_CHECK( ffPlaneNoise( aveNoise, &uvar, sftexist, aveNoiseInTime, antweights, backgroundScaling_slided, secondFFTplan, expRandVals, rng, &( ffdata->ffnormalization ) ) == XLAL_SUCCESS, XLAL_EFUNC );
778 if ( uvar.printData ) {
779 XLAL_CHECK( printREAL4Vector2File( ( REAL4Vector * )aveNoise, uvar.outdirectory, "ffbackground.dat" ) == XLAL_SUCCESS, XLAL_EFUNC );
780 }
781
782 //Compute the weighted TF data
783 REAL4VectorAligned *TFdata_weighted = NULL;
784 XLAL_CHECK( ( TFdata_weighted = XLALCreateREAL4VectorAligned( ffdata->numffts * ffdata->numfbins, 32 ) ) != NULL, XLAL_EFUNC );
785 XLAL_CHECK( tfWeight( TFdata_weighted, TFdata_slided, background_slided, antweights, backgroundScaling_slided, indexValuesOfExistingSFTs, &uvar ) == XLAL_SUCCESS, XLAL_EFUNC );
786 XLALDestroyREAL4VectorAligned( TFdata_slided );
787 XLALDestroyREAL4VectorAligned( antweights );
788
789 //Print out data product if requested
790 if ( uvar.printData ) {
791 XLAL_CHECK( printREAL4Vector2File( ( REAL4Vector * )TFdata_weighted, uvar.outdirectory, "procTFdata.dat" ) == XLAL_SUCCESS, XLAL_EFUNC );
792 }
793
794 //Calculation of average TF noise per frequency bin ratio to total mean
795 //this block of code does not avoid lines when computing the average F-bin ratio. Upper limits remain virtually unchanged
796 //when comaring runs that have line finding enabled or disabled
797 REAL4VectorAligned *aveTFnoisePerFbinRatio = NULL;
798 XLAL_CHECK( ( aveTFnoisePerFbinRatio = calcAveTFnoisePerFbinRatio( background_slided, backgroundScaling_slided, ffdata->numffts ) ) != NULL, XLAL_EFUNC );
799 XLALDestroyREAL4VectorAligned( background_slided );
800 XLALDestroyREAL4VectorAligned( backgroundScaling_slided );
801
802 //Do the second FFT
803 XLAL_CHECK( makeSecondFFT( ffdata, TFdata_weighted, secondFFTplan ) == XLAL_SUCCESS, XLAL_EFUNC );
804 //Normalize according to LAL PSD spec (also done in ffPlaneNoise() so this doesn't change anything)
805 //There is a secret divide by numffts in the weighting of the TF data (sumofweights), so we don't need to do it here;
806 //the numffts divisor gets squared when taking the PSD, so it is not applied here
807 for ( ii = 0; ii < ( INT4 )ffdata->ffdata->length; ii++ ) {
808 ffdata->ffdata->data[ii] *= uvar.Tobs;
809 }
810
811 REAL4 secFFTmean = 0.0, secFFTsigma = 0.0;
812 secFFTmean = calcMean( ffdata->ffdata );
813 XLAL_CHECK( calcStddev( &secFFTsigma, ffdata->ffdata ) == XLAL_SUCCESS, XLAL_EFUNC );
814
815 XLALDestroyREAL4VectorAligned( TFdata_weighted );
816
817 //Print out data product if requested
818 if ( uvar.printData ) {
819 XLAL_CHECK( printREAL4Vector2File( ( REAL4Vector * )( ffdata->ffdata ), uvar.outdirectory, "ffdata.dat" ) == XLAL_SUCCESS, XLAL_EFUNC );
820 }
821
822 fprintf( LOG, "2nd FFT ave = %g, 2nd FFT stddev = %g, expected ave = 1.0\n", secFFTmean, secFFTsigma );
823 fprintf( stderr, "2nd FFT ave = %g, 2nd FFT stddev = %g, expected ave = 1.0\n", secFFTmean, secFFTsigma );
824
825 //Exit with failure if there are no SFTs (probably this doesn't get hit)
826 XLAL_CHECK( secFFTmean != 0.0, XLAL_FAILURE, "Average second FFT power is 0.0. Perhaps no SFTs are remaining? Program exiting with failure.\n" );
827
828 //If the user wants to test a single, exact template, then we do that here
829 if ( uvar.templateTest ) {
830 if ( uvar.printData ) {
831 fprintf( stderr, "numfbins=%d, maxbinshift=%d, numffts=%d, numfprbins=%d\n", ffdata->numfbins, maxbinshift, ffdata->numffts, ffdata->numfprbins );
832 }
833 fprintf( stderr, "Testing template f=%f, P=%f, Df=%f\n", uvar.templateTestF, uvar.templateTestP, uvar.templateTestDf );
834 fprintf( LOG, "Testing template f=%f, P=%f, Df=%f\n", uvar.templateTestF, uvar.templateTestP, uvar.templateTestDf );
835
836 //Load template quantities into a test candidate
837 loadCandidateData( &( exactCandidates1->data[0] ), uvar.templateTestF, uvar.templateTestP, uvar.templateTestDf, dopplerpos.Alpha, dopplerpos.Delta, 0.0, 0.0, 0.0, 0, 0.0, -1, 0 );
838
839 //Resize the output candidate vector if necessary
840 if ( exactCandidates2->numofcandidates == exactCandidates2->length - 1 ) {
841 XLAL_CHECK( ( exactCandidates2 = resizecandidateVector( exactCandidates2, 2 * exactCandidates2->length ) ) != NULL, XLAL_EFUNC );
842 }
843
844 //Analyze the template stored in the test candidate
845 XLAL_CHECK( analyzeOneTemplate( &( exactCandidates2->data[exactCandidates2->numofcandidates] ), &( exactCandidates1->data[0] ), ffdata, aveNoise, aveTFnoisePerFbinRatio, &uvar, secondFFTplan, rng, !uvar.gaussTemplatesOnly ) == XLAL_SUCCESS, XLAL_EFUNC );
846 ( exactCandidates2->numofcandidates )++;
847
848 //Rescale the h0 output from the normaliztions and amount of observation time present
849 exactCandidates2->data[exactCandidates2->numofcandidates - 1].h0 /= sqrt( ffdata->tfnormalization ) * pow( frac_tobs_complete * ffdata->ffnormalization / skypointffnormalization, 0.25 );
850
851 } else if ( uvar.bruteForceTemplateTest ) {
852 candidate cand;
853 loadCandidateData( &cand, uvar.templateTestF, uvar.templateTestP, uvar.templateTestDf, dopplerpos.Alpha, dopplerpos.Delta, 0.0, 0.0, 0.0, 0, 0.0, -1, 0 );
854 TwoSpectParamSpaceSearchVals paramspace = {uvar.templateTestF - 2.0 / uvar.Tsft, uvar.templateTestF + 2.0 / uvar.Tsft, 21, 5, 5, 0.5, uvar.templateTestDf - 2.0 / uvar.Tsft,
855 uvar.templateTestDf + 2.0 / uvar.Tsft, 21
856 };
857 XLAL_CHECK( bruteForceTemplateTest( &( exactCandidates2 ), cand, &paramspace, &uvar, ffdata->ffdata, aveNoise, aveTFnoisePerFbinRatio, secondFFTplan, rng, 1 ) == XLAL_SUCCESS, XLAL_EFUNC );
858 for ( ii = 0; ii < ( INT4 )exactCandidates2->numofcandidates; ii++ ) {
859 exactCandidates2->data[ii].h0 /= sqrt( ffdata->tfnormalization ) * pow( frac_tobs_complete * ffdata->ffnormalization / skypointffnormalization, 0.25 );
860 }
861 }
862
863 //If the user wants to do a template search, that is done here
864 if ( uvar.templateSearch ) {
865 //fprintf(stderr, "Calling templateSearch\n (in development, last edited 2014-06-09)\n");
866 XLAL_CHECK( templateSearch_scox1Style( &exactCandidates2, uvar.fmin, uvar.fspan, uvar.templateSearchP, uvar.templateSearchAsini, uvar.templateSearchAsiniSigma, skypos, &uvar, ffdata->ffdata, aveNoise, aveTFnoisePerFbinRatio, trackedlines, secondFFTplan, rng, 1 ) == XLAL_SUCCESS, XLAL_EFUNC );
867 //fprintf(stderr, "Done calling templateSearch\n");
868 for ( ii = 0; ii < ( INT4 )exactCandidates2->numofcandidates; ii++ ) {
869 exactCandidates2->data[ii].h0 /= sqrt( ffdata->tfnormalization ) * pow( frac_tobs_complete * ffdata->ffnormalization / skypointffnormalization, 0.25 );
870 }
871 }
872
873 //If the user wants to do a templated search with fixed Df, that is done here
874 if ( uvar.templateSearchFixedDf ) {
875 XLAL_CHECK( templateSearch_fixedDf( &exactCandidates2, uvar.templateSearchDf, uvar.fmin, uvar.fspan, uvar.templateSearchP, skypos, &uvar, ffdata->ffdata, aveNoise, aveTFnoisePerFbinRatio, trackedlines, secondFFTplan, rng, 1 ) == XLAL_SUCCESS, XLAL_EFUNC );
876 for ( ii = 0; ii < ( INT4 )exactCandidates2->numofcandidates; ii++ ) {
877 exactCandidates2->data[ii].h0 /= sqrt( ffdata->tfnormalization ) * pow( frac_tobs_complete * ffdata->ffnormalization / skypointffnormalization, 0.25 );
878 }
879 }
880
881 //Template bank analysis
882 if ( XLALUserVarWasSet( &( uvar.templatebankfile ) ) ) {
883 candidateVector *candVec = NULL, *subsetVec = NULL;
884 XLAL_CHECK( ( candVec = createcandidateVector( 50 ) ) != NULL, XLAL_EFUNC );
885 XLAL_CHECK( ( subsetVec = createcandidateVector( 10 ) ) != NULL, XLAL_EFUNC );
886 XLAL_CHECK( testTwoSpectTemplateVector( candVec, templateVec, ffdata, aveNoise, aveTFnoisePerFbinRatio, skypos, &uvar, rng, 10 ) == XLAL_SUCCESS, XLAL_EFUNC );
887 XLAL_CHECK( analyzeCandidatesTemplateFromVector( subsetVec, candVec, templateVec, ffdata, aveNoise, aveTFnoisePerFbinRatio, &uvar, rng, 500 ) == XLAL_SUCCESS, XLAL_EFUNC );
888 for ( ii = 0; ii < ( INT4 )subsetVec->length; ii++ ) {
889 if ( exactCandidates2->numofcandidates == exactCandidates2->length ) {
890 XLAL_CHECK( ( exactCandidates2 = resizecandidateVector( exactCandidates2, 2 * exactCandidates2->length ) ) != NULL, XLAL_EFUNC );
891 }
892 loadCandidateData( &( exactCandidates2->data[exactCandidates2->numofcandidates] ), subsetVec->data[ii].fsig, subsetVec->data[ii].period, subsetVec->data[ii].moddepth, subsetVec->data[ii].ra, subsetVec->data[ii].dec, subsetVec->data[ii].stat, subsetVec->data[ii].h0, subsetVec->data[ii].prob, subsetVec->data[ii].proberrcode, subsetVec->data[ii].normalization, subsetVec->data[ii].templateVectorIndex, subsetVec->data[ii].lineContamination );
893 exactCandidates2->data[exactCandidates2->numofcandidates + ii].h0 /= sqrt( ffdata->tfnormalization ) * pow( frac_tobs_complete * ffdata->ffnormalization / skypointffnormalization, 0.25 ); //Scaling here
894 ( exactCandidates2->numofcandidates )++;
895 }
896
897 destroycandidateVector( candVec );
898 destroycandidateVector( subsetVec );
899 }
900
901 if ( uvar.signalOnly ) {
902 return 0;
903 }
904
905 //Start of the IHS step!
906 //Find the FAR of IHS sum -- only if the templateTest has not been given
907 candidateVector *ihsCandidates_reduced = NULL;
908 if ( !XLALUserVarWasSet( &uvar.templatebankfile ) && !uvar.templateTest && !uvar.templateSearch && !uvar.bruteForceTemplateTest && !uvar.templateSearchFixedDf ) {
909 //If the false alarm thresholds need to be computed
910 if ( ihsfarstruct->ihsfar->data[0] < 0.0 ) {
911 XLAL_CHECK( genIhsFar( ihsfarstruct, &uvar, maxrows, aveNoise, rng ) == XLAL_SUCCESS, XLAL_EFUNC );
912 }
913
914 //Run the IHS algorithm on the data
915 XLAL_CHECK( runIHS( ihsmaxima, ffdata, ihsfarstruct, &uvar, maxrows, aveNoise, aveTFnoisePerFbinRatio ) == XLAL_SUCCESS, XLAL_EFUNC );
916
917 //Find any IHS candidates
918 XLAL_CHECK( findIHScandidates( &ihsCandidates, ihsfarstruct, &uvar, ffdata, ihsmaxima, aveTFnoisePerFbinRatio, trackedlines ) == XLAL_SUCCESS, XLAL_EFUNC );
919
920 //If requested, keep only the most significant IHS candidates
921 if ( XLALUserVarWasSet( &uvar.keepOnlyTopNumIHS ) && ( INT4 )ihsCandidates->numofcandidates > uvar.keepOnlyTopNumIHS ) {
922 XLAL_CHECK( ( ihsCandidates_reduced = keepMostSignificantCandidates( ihsCandidates, &uvar ) ) != NULL, XLAL_EFUNC );
923 //Put ihsCandidates_reduced back into a reset ihsCandidates
924 ihsCandidates->numofcandidates = 0;
925 for ( ii = 0; ii < ( INT4 )ihsCandidates_reduced->numofcandidates; ii++ ) {
926 loadCandidateData( &( ihsCandidates->data[ii] ), ihsCandidates_reduced->data[ii].fsig, ihsCandidates_reduced->data[ii].period, ihsCandidates_reduced->data[ii].moddepth, dopplerpos.Alpha, dopplerpos.Delta, ihsCandidates_reduced->data[ii].stat, ihsCandidates_reduced->data[ii].h0, ihsCandidates_reduced->data[ii].prob, 0, ihsCandidates_reduced->data[ii].normalization, ihsCandidates_reduced->data[ii].templateVectorIndex, ihsCandidates_reduced->data[ii].lineContamination );
927 ( ihsCandidates->numofcandidates )++;
928 }
929 destroycandidateVector( ihsCandidates_reduced );
930 }
931 }
932 //End of the IHS step
933
934 //Start of the Gaussian template search!
935 //First check to see if the IHSonly or templateTest or templateSearch was given
936 if ( uvar.IHSonly && !uvar.templateTest && !uvar.templateSearch && !uvar.templateSearchFixedDf && !XLALUserVarWasSet( &uvar.templatebankfile ) ) {
937 //Check the length of the exactCandidates2 vector is large enough and resize if necessary
938 if ( exactCandidates2->length < exactCandidates2->numofcandidates + ihsCandidates->numofcandidates ) {
939 XLAL_CHECK( ( exactCandidates2 = resizecandidateVector( exactCandidates2, exactCandidates2->numofcandidates + ihsCandidates->numofcandidates ) ) != NULL, XLAL_EFUNC );
940 }
941
942 //Use the typical list
943 INT4 numofcandidatesalready = exactCandidates2->numofcandidates;
944 for ( ii = 0; ii < ( INT4 )ihsCandidates->numofcandidates; ii++ ) {
945 loadCandidateData( &( exactCandidates2->data[ii + numofcandidatesalready] ), ihsCandidates->data[ii].fsig, ihsCandidates->data[ii].period, ihsCandidates->data[ii].moddepth, dopplerpos.Alpha, dopplerpos.Delta, ihsCandidates->data[ii].stat, ihsCandidates->data[ii].h0, ihsCandidates->data[ii].prob, 0, ihsCandidates->data[ii].normalization, ihsCandidates->data[ii].templateVectorIndex, ihsCandidates->data[ii].lineContamination );
946 exactCandidates2->data[ii + numofcandidatesalready].h0 /= sqrt( ffdata->tfnormalization ) * pow( frac_tobs_complete * ffdata->ffnormalization / skypointffnormalization, 0.25 ); //Scaling here
947 ( exactCandidates2->numofcandidates )++;
948 }
949
950 } else if ( !uvar.templateTest && !uvar.templateSearch && !uvar.templateSearchFixedDf && !XLALUserVarWasSet( &uvar.templatebankfile ) ) {
951
952 //Test the IHS candidates against Gaussian templates in this function
953 XLAL_CHECK( testIHScandidates( &gaussCandidates1, ihsCandidates, ffdata, aveNoise, aveTFnoisePerFbinRatio, skypos, &uvar, rng ) == XLAL_SUCCESS, XLAL_EFUNC );
954
955 for ( ii = 0; ii < ( INT4 )gaussCandidates1->numofcandidates; ii++ ) {
956 fprintf( stderr, "Candidate %d: f0=%g, P=%g, df=%g\n", ii, gaussCandidates1->data[ii].fsig, gaussCandidates1->data[ii].period, gaussCandidates1->data[ii].moddepth );
957 }
958 } /* if IHSonly is not given && templateTest not given and templateSearch not given */
959 //End of the Gaussian template search
960
961 //Reset IHS candidates, but keep length the same (doesn't reset actual values in the vector)
962 ihsCandidates->numofcandidates = 0;
963
964 //Search the candidates further if the number of candidates passing the first Gaussian template test is greater than 0
965 if ( gaussCandidates1->numofcandidates > 0 ) {
966 //Start clustering! Note that the clustering algorithm takes care of the period range of parameter space
967 XLAL_CHECK( clusterCandidates( &gaussCandidates2, gaussCandidates1, ffdata, &uvar, aveNoise, aveTFnoisePerFbinRatio, rng, 0 ) == XLAL_SUCCESS, XLAL_EFUNC );
968
969 for ( ii = 0; ii < ( INT4 )gaussCandidates2->numofcandidates; ii++ ) {
970 fprintf( stderr, "Candidate %d: f0=%g, P=%g, df=%g\n", ii, gaussCandidates2->data[ii].fsig, gaussCandidates2->data[ii].period, gaussCandidates2->data[ii].moddepth );
971 }
972
973 //Reset first set of Gaussian template candidates, but keep length the same (doesn't reset actual values in the vector)
974 gaussCandidates1->numofcandidates = 0;
975
976 //Start detailed Gaussian template search!
977 for ( ii = 0; ii < ( INT4 )gaussCandidates2->numofcandidates; ii++ ) {
978
979 if ( gaussCandidates3->numofcandidates == gaussCandidates3->length - 1 ) {
980 XLAL_CHECK( ( gaussCandidates3 = resizecandidateVector( gaussCandidates3, 2 * gaussCandidates3->length ) ) != NULL, XLAL_EFUNC );
981 }
982
983 TwoSpectParamSpaceSearchVals paramspace = {gaussCandidates2->data[ii].fsig - 1.0 / uvar.Tsft, gaussCandidates2->data[ii].fsig + 1.0 / uvar.Tsft, 5, 2, 2, 1.0,
984 gaussCandidates2->data[ii].moddepth - 1.0 / uvar.Tsft, gaussCandidates2->data[ii].moddepth + 1.0 / uvar.Tsft, 5
985 };
986 XLAL_CHECK( bruteForceTemplateSearch( &( gaussCandidates3->data[gaussCandidates3->numofcandidates] ), gaussCandidates2->data[ii], &paramspace, &uvar, ffdata->ffdata,
987 aveNoise, aveTFnoisePerFbinRatio, secondFFTplan, rng, 0 ) == XLAL_SUCCESS, XLAL_EFUNC );
988 gaussCandidates3->numofcandidates++;
989
990 } /* for ii < numofcandidates */
991
992 for ( ii = 0; ii < ( INT4 )gaussCandidates3->numofcandidates; ii++ ) {
993 fprintf( stderr, "Candidate %d: f0=%g, P=%g, df=%g\n", ii, gaussCandidates3->data[ii].fsig, gaussCandidates3->data[ii].period, gaussCandidates3->data[ii].moddepth );
994 }
995
996 //Reset 2nd round of Gaussian template candidates, but keep length the same (doesn't reset actual values in the vector)
997 gaussCandidates2->numofcandidates = 0;
998
999 //Start clustering!
1000 XLAL_CHECK( clusterCandidates( &gaussCandidates4, gaussCandidates3, ffdata, &uvar, aveNoise, aveTFnoisePerFbinRatio, rng, 0 ) == XLAL_SUCCESS, XLAL_EFUNC );
1001
1002 for ( ii = 0; ii < ( INT4 )gaussCandidates4->numofcandidates; ii++ ) {
1003 fprintf( stderr, "Candidate %d: f0=%g, P=%g, df=%g\n", ii, gaussCandidates4->data[ii].fsig, gaussCandidates4->data[ii].period, gaussCandidates4->data[ii].moddepth );
1004 }
1005
1006 //Reset 3rd set of Gaussian template candidates
1007 gaussCandidates3->numofcandidates = 0;
1008
1009 //Initial check using "exact" template
1010 for ( ii = 0; ii < ( INT4 )gaussCandidates4->numofcandidates; ii++ ) {
1011 candidate cand;
1012 XLAL_CHECK( analyzeOneTemplate( &cand, &( gaussCandidates4->data[ii] ), ffdata, aveNoise, aveTFnoisePerFbinRatio, &uvar, secondFFTplan, rng, 1 ) == XLAL_SUCCESS, XLAL_EFUNC );
1013
1014 //Check that we are above threshold before storing the candidate
1015 if ( cand.prob < log10( uvar.tmplfar ) ) {
1016 if ( exactCandidates1->numofcandidates == exactCandidates1->length - 1 ) {
1017 XLAL_CHECK( ( exactCandidates1 = resizecandidateVector( exactCandidates1, 2 * exactCandidates1->length ) ) != NULL, XLAL_EFUNC );
1018 }
1019 loadCandidateData( &exactCandidates1->data[exactCandidates1->numofcandidates], cand.fsig, cand.period, cand.moddepth, ( REAL4 )dopplerpos.Alpha, ( REAL4 )dopplerpos.Delta, cand.stat, cand.h0, cand.prob, cand.proberrcode, cand.normalization, cand.templateVectorIndex, cand.lineContamination );
1020 exactCandidates1->numofcandidates++;
1021 }
1022 } /* for ii < numofcandidates */
1023 fprintf( LOG, "Number of candidates confirmed with exact templates = %d\n", exactCandidates1->numofcandidates );
1024 fprintf( stderr, "Number of candidates confirmed with exact templates = %d\n", exactCandidates1->numofcandidates );
1025 for ( ii = 0; ii < ( INT4 )exactCandidates1->numofcandidates; ii++ ) {
1026 fprintf( stderr, "Candidate %d: f0=%g, P=%g, df=%g\n", ii, exactCandidates1->data[ii].fsig, exactCandidates1->data[ii].period, exactCandidates1->data[ii].moddepth );
1027 }
1028 //Done with initial check
1029
1030 //Reset 4th set of Gaussian template candidates, but keep length the same (doesn't reset actual values in the vector)
1031 gaussCandidates4->numofcandidates = 0;
1032
1033 //Start detailed "exact" template search!
1034 for ( ii = 0; ii < ( INT4 )exactCandidates1->numofcandidates; ii++ ) {
1035
1036 if ( exactCandidates2->numofcandidates == exactCandidates2->length - 1 ) {
1037 XLAL_CHECK( ( exactCandidates2 = resizecandidateVector( exactCandidates2, 2 * exactCandidates2->length ) ) != NULL, XLAL_EFUNC );
1038 }
1039
1040 TwoSpectParamSpaceSearchVals paramspace = {exactCandidates1->data[ii].fsig - 0.5 / uvar.Tsft, exactCandidates1->data[ii].fsig + 0.5 / uvar.Tsft, 3, 1, 1, 1.0,
1041 exactCandidates1->data[ii].moddepth - 0.5 / uvar.Tsft, exactCandidates1->data[ii].moddepth + 0.5 / uvar.Tsft, 3
1042 };
1043 if ( !uvar.gaussTemplatesOnly ) {
1044 XLAL_CHECK( bruteForceTemplateSearch( &( exactCandidates2->data[exactCandidates2->numofcandidates] ), exactCandidates1->data[ii], &paramspace, &uvar, ffdata->ffdata, aveNoise, aveTFnoisePerFbinRatio, secondFFTplan, rng, 1 ) == XLAL_SUCCESS, XLAL_EFUNC );
1045 } else {
1046 XLAL_CHECK( bruteForceTemplateSearch( &( exactCandidates2->data[exactCandidates2->numofcandidates] ), exactCandidates1->data[ii], &paramspace, &uvar, ffdata->ffdata, aveNoise, aveTFnoisePerFbinRatio, secondFFTplan, rng, 0 ) == XLAL_SUCCESS, XLAL_EFUNC );
1047 }
1048 exactCandidates2->data[exactCandidates2->numofcandidates].h0 /= sqrt( ffdata->tfnormalization ) * pow( frac_tobs_complete * ffdata->ffnormalization / skypointffnormalization, 0.25 ); //Scaling here
1049 exactCandidates2->numofcandidates++;
1050
1051 fprintf( stderr, "Candidate %d: f0=%g, P=%g, df=%g\n", ii, exactCandidates2->data[ii].fsig, exactCandidates2->data[ii].period, exactCandidates2->data[ii].moddepth );
1052 } /* for ii < numofcandidates */
1053 //End of detailed search
1054
1055 //Reset first round of exact template candidates, but keep length the same (doesn't reset actual values in the vector)
1056 exactCandidates1->numofcandidates = 0;
1057
1058 fprintf( LOG, "Exact step is done with the total number of candidates = %d\n", exactCandidates2->numofcandidates );
1059 fprintf( stderr, "Exact step is done with the total number of candidates = %d\n", exactCandidates2->numofcandidates );
1060
1061 } /* if gaussCandidates1->numofcandidates > 0 */
1062
1063 //Determine upper limits, if the ULoff has not been set
1064 if ( !uvar.ULoff && !uvar.templateTest && !uvar.templateSearch && !uvar.templateSearchFixedDf && !XLALUserVarWasSet( &uvar.templatebankfile ) ) {
1065 upperlimits->data[upperlimits->length - 1].alpha = ( REAL4 )dopplerpos.Alpha;
1066 upperlimits->data[upperlimits->length - 1].delta = ( REAL4 )dopplerpos.Delta;
1067 upperlimits->data[upperlimits->length - 1].normalization = ffdata->tfnormalization;
1068
1069 XLAL_CHECK( skypoint95UL( &( upperlimits->data[upperlimits->length - 1] ), &uvar, ffdata, ihsmaxima, ihsfarstruct, aveTFnoisePerFbinRatio ) == XLAL_SUCCESS, XLAL_EFUNC );
1070
1071 for ( ii = 0; ii < ( INT4 )upperlimits->data[upperlimits->length - 1].ULval->length; ii++ ) {
1072 upperlimits->data[upperlimits->length - 1].ULval->data[ii] /= sqrt( ffdata->tfnormalization ) * pow( frac_tobs_complete * ffdata->ffnormalization / skypointffnormalization, 0.25 ); //Compensation for different duty cycle and antenna pattern weights
1073 }
1074
1075 XLAL_CHECK( ( upperlimits = resizeUpperLimitVector( upperlimits, upperlimits->length + 1 ) ) != NULL, XLAL_EFUNC );
1076 } /* if producing UL */
1077
1078 //Destroy stuff
1079 XLALDestroyREAL4VectorAligned( aveTFnoisePerFbinRatio );
1080 if ( trackedlines != NULL ) {
1081 XLALDestroyREAL4VectorSequence( trackedlines );
1082 }
1083 if ( detectors->length > 1 && !uvar.signalOnly ) {
1084 XLALDestroyINT4Vector( sftexist );
1085 XLALDestroyINT4Vector( indexValuesOfExistingSFTs );
1086 }
1087
1088 //Iterate to next sky location
1089 XLAL_CHECK( XLALNextDopplerSkyPos( &dopplerpos, &scan ) == XLAL_SUCCESS, XLAL_EFUNC );
1090
1091 } /* while sky scan is not finished */
1092
1093 if ( exactCandidates2->numofcandidates != 0 ) {
1094 fprintf( LOG, "\n**Report of candidates:**\n" );
1095 fprintf( stderr, "\n**Report of candidates:**\n" );
1096
1097 for ( ii = 0; ii < ( INT4 )exactCandidates2->numofcandidates; ii++ ) {
1098 fprintf( LOG, "fsig = %.6f, period = %.6f, df = %.7f, RA = %.4f, DEC = %.4f, R = %.4f, h0 = %g, Prob = %.4f, TF norm = %g, template vec num = %d, line contamination = %d\n", exactCandidates2->data[ii].fsig, exactCandidates2->data[ii].period, exactCandidates2->data[ii].moddepth, exactCandidates2->data[ii].ra, exactCandidates2->data[ii].dec, exactCandidates2->data[ii].stat, exactCandidates2->data[ii].h0, exactCandidates2->data[ii].prob, ffdata->tfnormalization, exactCandidates2->data[ii].templateVectorIndex, exactCandidates2->data[ii].lineContamination );
1099 fprintf( stderr, "fsig = %.6f, period = %.6f, df = %.7f, RA = %.4f, DEC = %.4f, R = %.4f, h0 = %g, Prob = %.4f, TF norm = %g, template vec num = %d, line contamination = %d\n", exactCandidates2->data[ii].fsig, exactCandidates2->data[ii].period, exactCandidates2->data[ii].moddepth, exactCandidates2->data[ii].ra, exactCandidates2->data[ii].dec, exactCandidates2->data[ii].stat, exactCandidates2->data[ii].h0, exactCandidates2->data[ii].prob, ffdata->tfnormalization, exactCandidates2->data[ii].templateVectorIndex, exactCandidates2->data[ii].lineContamination );
1100 } /* for ii < exactCandidates2->numofcandidates */
1101 } /* if exactCandidates2->numofcandidates != 0 */
1102
1103 //Output upper limits to a file, if ULoff is not given
1104 if ( !uvar.ULoff ) {
1105 for ( ii = 0; ii < ( INT4 )upperlimits->length - 1; ii++ ) {
1106 XLAL_CHECK( outputUpperLimitToFile( ULFILENAME->data, upperlimits->data[ii], uvar.allULvalsPerSkyLoc ) == XLAL_SUCCESS, XLAL_EFUNC );
1107 }
1108 }
1109
1110 //Output candidates to a file
1111 XLAL_CHECK( writeCandidateVector2File( CANDFILENAME->data, exactCandidates2 ) == XLAL_SUCCESS, XLAL_EFUNC );
1112
1113 //Destroy varaibles
1114 if ( detectors->length > 1 ) {
1115 XLALDestroyMultiSFTVector( multiSFTvector );
1116 } else {
1117 XLALDestroyREAL4VectorAligned( usableTFdata );
1118 if ( !uvar.signalOnly ) {
1119 XLALDestroyINT4Vector( sftexist );
1120 XLALDestroyINT4Vector( indexValuesOfExistingSFTs );
1121 }
1122 }
1123 if ( XLALUserVarWasSet( &( uvar.templatebankfile ) ) ) {
1124 destroyTwoSpectTemplateVector( templateVec );
1125 }
1127 gsl_rng_free( rng );
1128 XLALDestroyREAL4VectorAligned( background );
1129 XLALDestroyREAL4VectorAligned( background0 );
1130 XLALDestroyREAL4VectorAligned( backgroundForihs2h0 );
1131 XLALDestroyREAL4VectorAligned( antweightsforihs2h0 );
1132 XLALDestroyINT4Vector( sftexistForihs2h0 );
1134 XLALDestroyREAL4VectorAligned( aveNoiseInTime );
1135 XLALDestroyREAL4VectorAligned( aveNoiseInTimeForihs2h0 );
1136 XLALDestroyREAL4VectorAligned( expRandVals );
1137 XLALDestroyREAL4VectorAligned( backgroundScaling );
1138 XLALDestroyREAL4VectorAligned( backgroundScalingForihs2h0 );
1139 destroyREAL4VectorAlignedArray( backgroundRatioMeans );
1140 XLALDestroyINT4Vector( binshifts );
1142 XLALDestroyTimestampVector( jointSFTtimestamps );
1143 XLALDestroyMultiDetectorStateSeries( multiStateSeries );
1144 destroyffdata( ffdata );
1145 destroyihsfarStruct( ihsfarstruct );
1146 destroyihsMaxima( ihsmaxima );
1147 XLALDestroyREAL4FFTPlan( secondFFTplan );
1150 XLALDestroyCHARVector( ULFILENAME );
1151 XLALDestroyCHARVector( LOGFILENAME );
1152 XLALDestroyCHARVector( CANDFILENAME );
1153 destroyUpperLimitVector( upperlimits );
1154 destroycandidateVector( ihsCandidates );
1155 destroycandidateVector( gaussCandidates1 );
1156 destroycandidateVector( gaussCandidates2 );
1157 destroycandidateVector( gaussCandidates3 );
1158 destroycandidateVector( gaussCandidates4 );
1159 destroycandidateVector( exactCandidates1 );
1160 destroycandidateVector( exactCandidates2 );
1161 FreeDopplerSkyScan( &status, &scan );
1162
1163 //print end time
1164 time( &programendtime );
1165 ptm = localtime( &programendtime );
1166 fprintf( stderr, "Program finished on %s", asctime( ptm ) );
1167 fprintf( LOG, "Program finished on %s", asctime( ptm ) );
1168
1169 fclose( LOG );
1170
1171 //Check for leaks
1173
1174 //The end!
1175 return 0;
1176
1177} /* main() */
1178
1179
1180/**
1181 * Create a new frequency-frequency data structure for the TwoSpect analysis
1182 * \param [in] params Pointer to the UserInput_t
1183 * \return Pointer to new ffdataStruct
1184 */
1186{
1187
1189
1190 ffdataStruct *ffdata = NULL;
1191 XLAL_CHECK_NULL( ( ffdata = XLALMalloc( sizeof( *ffdata ) ) ) != NULL, XLAL_ENOMEM );
1192
1193 ffdata->numfbins = ( INT4 )( round( params->fspan * params->Tsft + 2.0 * params->dfmax * params->Tsft ) + 12 + 1 );
1194 ffdata->numffts = ( INT4 )floor( params->Tobs / ( params->Tsft - params->SFToverlap ) - 1 );
1195 ffdata->numfprbins = ( INT4 )floorf( ffdata->numffts * 0.5 ) + 1;
1196 ffdata->tfnormalization = ffdata->ffnormalization = 0.0;
1197
1198 XLAL_CHECK_NULL( ( ffdata->ffdata = XLALCreateREAL4VectorAligned( ffdata->numfbins * ffdata->numfprbins, 32 ) ) != NULL, XLAL_EFUNC );
1199
1200 return ffdata;
1201
1202} // createffdata()
1203
1204
1205/**
1206 * Free the frequency-frequency data structure
1207 * \param [in] data Pointer to the ffdataStruct
1208 */
1210{
1211 if ( data ) {
1213 XLALFree( ( ffdataStruct * )data );
1214 }
1215} // destroyffdata()
1216
1217
1218/**
1219 * Line detection algorithm
1220 *
1221 * 1. Calculate RMS for each fbin as function of time
1222 * 2. Calculate running median of RMS values
1223 * 3. Divide RMS values by running median. Lines stick out by being >>1
1224 * 4. Add up number of lines above threshold and report
1225 *
1226 * \param [in] TFdata Pointer to REAL4VectorAligned of SFT powers
1227 * \param [in] ffdata Pointer to ffdataStruct
1228 * \param [in] params Pointer to UserInput_t
1229 * \return Pointer to INT4Vector of bin number of the lines
1230 */
1232{
1233
1234 XLAL_CHECK_NULL( TFdata != NULL && ffdata != NULL && params != NULL, XLAL_EINVAL );
1235
1236 fprintf( stderr, "Running line detection algorithm... " );
1237
1239 INT4 blksize = 11;
1240 INT4 numlines = 0;
1241 INT4Vector *lines = NULL;
1242 //INT4 totalnumfbins = ffdata->numfbins+(params->blksize-1)+2*params->maxbinshift;
1243 INT4 numffts = ( INT4 )floor( params->Tobs / ( params->Tsft - params->SFToverlap ) - 1 ); //Number of FFTs
1244 INT4 totalnumfbins = ( INT4 )TFdata->length / numffts;
1245
1246 //Blocksize of running median
1247 LALRunningMedianPar block = {blksize};
1248
1249 //Compute weights
1250 REAL4VectorAligned *sftdata = NULL, *weights = NULL;
1251 XLAL_CHECK_NULL( ( sftdata = XLALCreateREAL4VectorAligned( totalnumfbins, 32 ) ) != NULL, XLAL_EFUNC );
1252 XLAL_CHECK_NULL( ( weights = XLALCreateREAL4VectorAligned( ffdata->numffts, 32 ) ) != NULL, XLAL_EFUNC );
1253 memset( weights->data, 0, ffdata->numffts * sizeof( REAL4 ) );
1254 REAL4 sumweights = 0.0;
1255 for ( INT4 ii = 0; ii < ffdata->numffts; ii++ ) {
1256 if ( TFdata->data[ii * totalnumfbins] != 0.0 ) {
1257 memcpy( sftdata->data, &( TFdata->data[ii * totalnumfbins] ), totalnumfbins * sizeof( REAL4 ) );
1258 REAL4 stddev = 0.0;
1259 XLAL_CHECK_NULL( calcStddev( &stddev, sftdata ) == XLAL_SUCCESS, XLAL_EFUNC );
1260 weights->data[ii] = 1.0 / ( stddev * stddev );
1261 sumweights += weights->data[ii];
1262 }
1263 }
1264 REAL4 invsumweights = 1.0 / sumweights;
1266
1267 //Compute RMS for each frequency bin as a function of time
1268 REAL4VectorAligned *testRMSvals = NULL, *testRngMedian = NULL, *testTSofPowers = NULL;
1269 XLAL_CHECK_NULL( ( testRMSvals = XLALCreateREAL4VectorAligned( totalnumfbins, 32 ) ) != NULL, XLAL_EFUNC );
1270 XLAL_CHECK_NULL( ( testRngMedian = XLALCreateREAL4VectorAligned( totalnumfbins - ( blksize - 1 ), 32 ) ) != NULL, XLAL_EFUNC );
1271 XLAL_CHECK_NULL( ( testTSofPowers = XLALCreateREAL4VectorAligned( ffdata->numffts, 32 ) ) != NULL, XLAL_EFUNC );
1272
1273 for ( UINT4 ii = 0; ii < testRMSvals->length; ii++ ) {
1274 for ( INT4 jj = 0; jj < ffdata->numffts; jj++ ) {
1275 testTSofPowers->data[jj] = TFdata->data[jj * testRMSvals->length + ii] * weights->data[jj] * invsumweights;
1276 }
1277 XLAL_CHECK_NULL( calcRms( &( testRMSvals->data[ii] ), testTSofPowers ) == XLAL_SUCCESS, XLAL_EFUNC ); //This approaches calcMean(TSofPowers) for stationary noise
1278 }
1279
1280 //Running median of RMS values
1281 LALSRunningMedian2( &status, ( REAL4Vector * )testRngMedian, ( REAL4Vector * )testRMSvals, block );
1282 XLAL_CHECK_NULL( status.statusCode == 0, XLAL_EFUNC );
1283
1284 //Determine which bins are above the threshold and store the bin number of the line
1285 for ( UINT4 ii = 0; ii < testRngMedian->length; ii++ ) {
1286 REAL4 normrmsval = testRMSvals->data[ii + ( blksize - 1 ) / 2] / testRngMedian->data[ii];
1287 if ( ( INT4 )( ii + ( blksize - 1 ) / 2 ) > ( ( params->blksize - 1 ) / 2 ) && normrmsval > params->lineDetection ) {
1288 XLAL_CHECK_NULL( ( lines = XLALResizeINT4Vector( lines, numlines + 1 ) ) != NULL, XLAL_EFUNC );
1289 lines->data[numlines] = ii + ( blksize - 1 ) / 2;
1290 numlines++;
1291 }
1292 }
1293
1294 //Destroy stuff
1295 XLALDestroyREAL4VectorAligned( testTSofPowers );
1296 XLALDestroyREAL4VectorAligned( testRngMedian );
1297 XLALDestroyREAL4VectorAligned( testRMSvals );
1299
1300 fprintf( stderr, "done\n" );
1301
1302 if ( lines != NULL ) {
1303 fprintf( stderr, "WARNING: %d line(s) found.\n", lines->length );
1304 }
1305
1306 return lines;
1307
1308} // detectLines_simple()
1309
1310/**
1311 * Track the lines for the sky position
1312 * \param [in] lines Pointer to INT4Vector with list of lines
1313 * \param [in] binshifts Pointer to INT4Vector of SFT bin shifts
1314 * \param [in] minfbin Frequency value of the lowest frequency bin
1315 * \param [in] df Spacing of frequency bins (typically 1/Tsft)
1316 * \return Pointer to REAL4VectorSequence containing lower/mid/upper frequency for each line
1317 */
1318REAL4VectorSequence *trackLines( const INT4Vector *lines, const INT4Vector *binshifts, const REAL4 minfbin, const REAL4 df )
1319{
1320
1321 XLAL_CHECK_NULL( lines != NULL && binshifts != NULL, XLAL_EINVAL );
1322
1323 fprintf( stderr, "Tracking lines... " );
1324
1327
1328 //REAL4 df = 1.0/params->Tsft;
1329 //REAL4 minfbin = (REAL4)(round(params->fmin*params->Tsft - params->dfmax*params->Tsft - 6.0 - 0.5*(params->blksize-1) - (REAL8)(maxbinshift))/params->Tsft);
1330
1331 UINT4 maxshiftindex = 0, minshiftindex = 0;
1332 min_max_index_INT4Vector( binshifts, &minshiftindex, &maxshiftindex );
1333 INT4 maxshift = binshifts->data[maxshiftindex], minshift = binshifts->data[minshiftindex];
1334
1335 for ( UINT4 ii = 0; ii < lines->length; ii++ ) {
1336 output->data[ii * 3] = lines->data[ii] * df + minfbin;
1337 //output->data[ii*3 + 1] = (lines->data[ii] + minshift)*df + minfbin;
1338 //output->data[ii*3 + 2] = (lines->data[ii] + maxshift)*df + minfbin;
1339 output->data[ii * 3 + 1] = ( lines->data[ii] + ( minshift - 1 ) ) * df + minfbin; //Add one extra bin for buffer
1340 output->data[ii * 3 + 2] = ( lines->data[ii] + ( maxshift + 1 ) ) * df + minfbin; //Add one extra bin for buffer
1341 }
1342
1343 fprintf( stderr, "done\n" );
1344
1345 return output;
1346
1347} // trackLines()
1348
1349
1350/**
1351 * Test algorithm to clean lines. NOT FULLY TESTED!
1352 * \param [in,out] TFdata Pointer to time-frequency data in a REAL4Vector to be cleaned
1353 * \param [in] background Pointer to REAL4Vector of running mean data
1354 * \param [in] lines Pointer to INT4Vector of lines
1355 * \param [in] params Pointer to UserInput_t
1356 * \param [in] rng Pointer to gsl_rng
1357 * \return Status value
1358 */
1359INT4 cleanLines( REAL4VectorAligned *TFdata, const REAL4VectorAligned *background, const INT4Vector *lines, const UserInput_t *params, const gsl_rng *rng )
1360{
1361
1362 if ( lines == NULL ) {
1363 return XLAL_SUCCESS;
1364 }
1365
1366 INT4 numffts = ( INT4 )floor( params->Tobs / ( params->Tsft - params->SFToverlap ) - 1 ); //Number of FFTs
1367 INT4 numfbins = ( INT4 )TFdata->length / numffts; //Number of frequency bins
1368
1369 REAL8 prevnoiseval = 0.0;
1370 for ( INT4 ii = 0; ii < numffts; ii++ ) {
1371 if ( TFdata->data[ii * numfbins] != 0.0 ) {
1372 for ( UINT4 jj = 0; jj < lines->length; jj++ ) {
1373 if ( lines->data[jj] - ( params->blksize - 1 ) / 2 >= 0 ) {
1374 REAL8 noiseval = expRandNum( background->data[ii * ( numfbins - ( params->blksize - 1 ) ) + lines->data[jj]], rng );
1375 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC, "ii=%d, jj=%d, background=%g", ii, jj, background->data[ii * ( numfbins - ( params->blksize - 1 ) ) + lines->data[jj]] );
1376 if ( ii > 0 && TFdata->data[( ii - 1 ) * numfbins] != 0.0 ) {
1377 noiseval *= ( 1.0 - 0.167 * 0.167 );
1378 noiseval += 0.167 * prevnoiseval;
1379 }
1380 TFdata->data[ii * numfbins + lines->data[jj]] = noiseval;
1381 prevnoiseval = noiseval;
1382 }
1383 } //end loop over vector of lines
1384 } //if SFT exists
1385 } //end loop over SFTs
1386
1387 return XLAL_SUCCESS;
1388
1389} // cleanLines()
1390
1391/**
1392 * Calculate the ratio of the average SFT noise to the mean of the average SFT noise
1393 * \param [in] background Pointer to REAL4VectorAligned of the running means
1394 * \param [in] backgroundScaling Pointer to REAL4VectorAligned of background scaling values
1395 * \param [in] numffts Number of SFTs in the total observation time
1396 * \return Pointer to REAL4VectorAligned of the ratio of values to the band mean
1397 */
1398REAL4VectorAligned *calcAveTFnoisePerFbinRatio( const REAL4VectorAligned *background, const REAL4VectorAligned *backgroundScaling, const UINT4 numffts )
1399{
1400 XLAL_CHECK_NULL( background != NULL, XLAL_EINVAL );
1401 UINT4 numfbins = background->length / numffts;
1402 REAL4VectorAligned *aveTFnoisePerFbinRatio = NULL, *TSofPowers = NULL;
1403 XLAL_CHECK_NULL( ( aveTFnoisePerFbinRatio = XLALCreateREAL4VectorAligned( numfbins, 32 ) ) != NULL, XLAL_EFUNC );
1404 XLAL_CHECK_NULL( ( TSofPowers = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1405 memset( TSofPowers->data, 0, sizeof( REAL4 )*TSofPowers->length );
1406 for ( UINT4 ii = 0; ii < numfbins; ii++ ) {
1407 REAL4 totalweightval = 0.0;
1408 for ( UINT4 jj = 0; jj < numffts; jj++ ) {
1409 if ( background->data[jj * numfbins + ii] != 0.0 && backgroundScaling->data[jj * numfbins + ii] != 0.0 ) {
1410 TSofPowers->data[jj] = 1.0 / ( background->data[jj * numfbins + ii] * backgroundScaling->data[jj * numfbins + ii] );
1411 totalweightval += 1.0 / ( background->data[jj * numfbins + ii] * background->data[jj * numfbins + ii] * backgroundScaling->data[jj * numfbins + ii] * backgroundScaling->data[jj * numfbins + ii] );
1412 }
1413 }
1414 aveTFnoisePerFbinRatio->data[ii] = 0.0;
1415 for ( UINT4 jj = 0; jj < numffts; jj++ ) {
1416 aveTFnoisePerFbinRatio->data[ii] += TSofPowers->data[jj];
1417 }
1418 aveTFnoisePerFbinRatio->data[ii] /= totalweightval;
1419 }
1420 REAL4 aveTFaveinv = 1.0 / calcMean( aveTFnoisePerFbinRatio );
1421 for ( UINT4 ii = 0; ii < numfbins; ii++ ) {
1422 aveTFnoisePerFbinRatio->data[ii] *= aveTFaveinv;
1423 }
1424 XLALDestroyREAL4VectorAligned( TSofPowers );
1425 return aveTFnoisePerFbinRatio;
1426}
1427
1428/**
1429 * Compute the second Fourier transform for TwoSpect
1430 * \param [out] output Pointer to the ffdataStruct to the containers for the second FFT
1431 * \param [in] tfdata Pointer REAL4VectorAligned of mean subtracted and weighted data
1432 * \param [in] plan Pointer to REAL4FFTPlan
1433 * \return Status value
1434 */
1435INT4 makeSecondFFT( ffdataStruct *output, REAL4VectorAligned *tfdata, const REAL4FFTPlan *plan )
1436{
1437
1438 XLAL_CHECK( output != NULL && tfdata != NULL && plan != NULL, XLAL_EINVAL );
1439
1440 fprintf( stderr, "Computing second FFT over SFTs... " );
1441
1442 REAL8 winFactor = 8.0 / 3.0;
1443
1444 //Do the second FFT
1445 REAL4VectorAligned *x = NULL, *psd = NULL, *windowData = NULL;
1446 REAL4Window *win = NULL;
1447 XLAL_CHECK( ( x = XLALCreateREAL4VectorAligned( output->numffts, 32 ) ) != NULL, XLAL_EFUNC );
1448 XLAL_CHECK( ( psd = XLALCreateREAL4VectorAligned( ( UINT4 )floor( x->length * 0.5 ) + 1, 32 ) ) != NULL, XLAL_EFUNC );
1449 XLAL_CHECK( ( win = XLALCreateHannREAL4Window( x->length ) ) != NULL, XLAL_EFUNC );
1450 XLAL_CHECK( ( windowData = XLALCreateREAL4VectorAligned( win->data->length, 32 ) ) != NULL, XLAL_EFUNC );
1451 memcpy( windowData->data, win->data->data, sizeof( REAL4 )*windowData->length );
1452
1453 for ( INT4 ii = 0; ii < output->numfbins; ii++ ) {
1454
1455 //Next, loop over times and pick the right frequency bin for each FFT and window
1456 //for (UINT4 jj=0; jj<x->length; jj++) x->data[jj] = tfdata->data[ii + jj*output->numfbins] * win->data->data[jj];
1457 //fastSSVectorMultiply_with_stride_and_offset(x, tfdata, windowData, output->numfbins, 1, ii, 0);
1458 for ( UINT4 jj = 0; jj < x->length; jj++ ) {
1459 x->data[jj] = tfdata->data[ii + jj * output->numfbins];
1460 }
1461 XLAL_CHECK( XLALVectorMultiplyREAL4( x->data, x->data, windowData->data, x->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1462
1463 //Make the FFT
1465
1466 //Fix beginning and end values if even, otherwise just the beginning if odd
1467 if ( GSL_IS_EVEN( x->length ) == 1 ) {
1468 psd->data[0] *= 2.0;
1469 psd->data[psd->length - 1] *= 2.0;
1470 } else {
1471 psd->data[0] *= 2.0;
1472 }
1473
1474 //Scale the data points by 1/N and window factor and (1/fs)
1475 //Order of vector is by second frequency then first frequency
1476 //It is possible that when dealing with very loud signals, lines, injections, etc. (e.g., far above the background)
1477 //then the output power here can be "rounded" because of the cast to nearby integer values.
1478 //For high (but not too high) power values, this may not be noticed because the cast can round to nearby decimal values.
1479 for ( UINT4 jj = 0; jj < psd->length; jj++ ) {
1480 output->ffdata->data[psd->length * ii + jj] = ( REAL4 )( psd->data[jj] * winFactor * output->ffnormalization );
1481 }
1482
1483 } /* for ii < numfbins */
1484
1485 //Destroy stuff
1488 XLALDestroyREAL4VectorAligned( windowData );
1490
1491 fprintf( stderr, "done\n" );
1492
1493 return XLAL_SUCCESS;
1494
1495} /* makeSecondFFT() */
1496
1497
1498/**
1499 * Determine the average of the noise power in each frequency bin across the band
1500 * \param [in] backgrnd Pointer to REAL4VectorAligned of the running mean values
1501 * \param [in] numfbins Number of frequency bins in the SFTs
1502 * \param [in] numffts Number of SFTs in the observation time
1503 * \param [in] binmin Minimum SFT bin to look at with this algorithm
1504 * \param [in] binmax Maximum SFT bin to look at with this algorithm
1505 * \return The average of the noise power across the frequency band
1506 */
1507REAL4 avgTFdataBand( const REAL4VectorAligned *backgrnd, UINT4 numfbins, UINT4 numffts, UINT4 binmin, UINT4 binmax )
1508{
1509
1510 XLAL_CHECK_REAL4( backgrnd != NULL && binmax > binmin, XLAL_EINVAL );
1511
1512 REAL4VectorAligned *aveNoiseInTime = NULL, *rngMeansOverBand = NULL;
1513 XLAL_CHECK_REAL4( ( aveNoiseInTime = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1514 XLAL_CHECK_REAL4( ( rngMeansOverBand = XLALCreateREAL4VectorAligned( binmax - binmin, 32 ) ) != NULL, XLAL_EFUNC );
1515
1516 for ( UINT4 ii = 0; ii < aveNoiseInTime->length; ii++ ) {
1517 memcpy( rngMeansOverBand->data, &( backgrnd->data[ii * numfbins + binmin] ), sizeof( *rngMeansOverBand->data )*rngMeansOverBand->length );
1518 aveNoiseInTime->data[ii] = calcMean( rngMeansOverBand );
1519 } /* for ii < aveNoiseInTime->length */
1520
1521 REAL4 avgTFdata = calcMean( aveNoiseInTime );
1522
1523 //Destroy stuff
1524 XLALDestroyREAL4VectorAligned( aveNoiseInTime );
1525 XLALDestroyREAL4VectorAligned( rngMeansOverBand );
1526
1527 return avgTFdata;
1528
1529} /* avgTFdataBand() */
1530
1531
1532/**
1533 * Determine the rms of the noise power in each frequency bin across the band
1534 * \param [in] backgrnd Pointer to REAL4Vector of the running mean values
1535 * \param [in] numfbins Number of frequency bins in the SFTs
1536 * \param [in] numffts Number of SFTs in the observation time
1537 * \param [in] binmin Minimum SFT bin to look at with this algorithm
1538 * \param [in] binmax Maximum SFT bin to look at with this algorithm
1539 * \return The RMS of the noise power across the frequency band
1540 */
1541REAL4 rmsTFdataBand( const REAL4VectorAligned *backgrnd, UINT4 numfbins, UINT4 numffts, UINT4 binmin, UINT4 binmax )
1542{
1543
1544 XLAL_CHECK_REAL4( backgrnd != NULL && binmax > binmin, XLAL_EINVAL );
1545
1546 REAL4VectorAligned *aveNoiseInTime = NULL, *rngMeansOverBand = NULL;
1547 XLAL_CHECK_REAL4( ( aveNoiseInTime = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1548 XLAL_CHECK_REAL4( ( rngMeansOverBand = XLALCreateREAL4VectorAligned( binmax - binmin, 32 ) ) != NULL, XLAL_EFUNC );
1549
1550 for ( UINT4 ii = 0; ii < aveNoiseInTime->length; ii++ ) {
1551 memcpy( rngMeansOverBand->data, &( backgrnd->data[ii * numfbins + binmin] ), sizeof( *rngMeansOverBand->data )*rngMeansOverBand->length );
1552 aveNoiseInTime->data[ii] = calcMean( rngMeansOverBand );
1553 } /* for ii < aveNoiseInTime->length */
1554
1555 REAL4 rmsTFdata = 0.0;
1556 XLAL_CHECK_REAL4( calcRms( &rmsTFdata, aveNoiseInTime ) == XLAL_SUCCESS, XLAL_EFUNC );
1557
1558 //Destroy stuff
1559 XLALDestroyREAL4VectorAligned( aveNoiseInTime );
1560 XLALDestroyREAL4VectorAligned( rngMeansOverBand );
1561
1562 return rmsTFdata;
1563
1564} /* rmsTFdataBand() */
1565
1566INT4 medianBackgroundBandInTime( REAL4VectorAligned *aveNoiseInTime, const REAL4VectorAligned *backgrnd, const INT4Vector *sftexist )
1567{
1568 XLAL_CHECK( aveNoiseInTime != NULL && backgrnd != NULL && sftexist != NULL, XLAL_EINVAL );
1569 memset( aveNoiseInTime->data, 0, sizeof( REAL4 )*aveNoiseInTime->length );
1570 REAL4VectorAligned *band = NULL;
1571 XLAL_CHECK( ( band = XLALCreateREAL4VectorAligned( backgrnd->length / sftexist->length, 32 ) ) != NULL, XLAL_EFUNC );
1572 for ( UINT4 ii = 0; ii < sftexist->length; ii++ ) {
1573 if ( sftexist->data[ii] != 0 ) {
1574 memcpy( band->data, &( backgrnd->data[ii * band->length] ), sizeof( REAL4 )*band->length );
1575 XLAL_CHECK( calcMedian( &( aveNoiseInTime->data[ii] ), band ) == XLAL_SUCCESS, XLAL_EFUNC );
1576 }
1577 }
1579 return XLAL_SUCCESS;
1580}
1581
1582/**
1583 * Measure of the average noise power in each 2nd FFT frequency bin
1584 * \param [out] aveNoise Pointer to REAL4VectorAligned of the expected 2nd FFT powers
1585 * \param [in] params Pointer to UserInput_t
1586 * \param [in] sftexist Pointer to INT4Vector of SFTs existing or not
1587 * \param [in] aveNoiseInTime Pointer to REAL4VectorAligned of running means
1588 * \param [in] antweights Pointer to REAL4VectorAligned of antenna pattern weights
1589 * \param [in] backgroundScaling Pointer to REAL4VectorAligned of background scaling values
1590 * \param [in] plan Pointer to REAL4FFTPlan
1591 * \param [in] expDistVals Pointer to REAL4VectorAligned of precomputed exponentially distributed random values to sample from
1592 * \param [in] rng Pointer to gsl_rng
1593 * \param [in,out] normalization Pointer to REAL8 value of the normalization for the 2nd FFT
1594 * \return Status value
1595 */
1596INT4 ffPlaneNoise( REAL4VectorAligned *aveNoise, const UserInput_t *params, const INT4Vector *sftexist, const REAL4VectorAligned *aveNoiseInTime, const REAL4VectorAligned *antweights, const REAL4VectorAligned *backgroundScaling, const REAL4FFTPlan *plan, const REAL4VectorAligned *expDistVals, const gsl_rng *rng, REAL8 *normalization )
1597{
1598
1599 XLAL_CHECK( aveNoise != NULL && params != NULL && sftexist != NULL && aveNoiseInTime != NULL && antweights != NULL && backgroundScaling != NULL && plan != NULL && normalization != NULL, XLAL_EINVAL );
1600
1601 fprintf( stderr, "Computing noise background estimate... " );
1602
1603 REAL8 invsumofweights = 0.0, sumofweights = 0.0;
1604
1605 UINT4 numffts = antweights->length;
1606 UINT4 numfbins = backgroundScaling->length / numffts;
1607 UINT4 numfprbins = aveNoise->length;
1608
1609 //Set up for making the PSD
1610 memset( aveNoise->data, 0, sizeof( REAL4 )*aveNoise->length );
1611
1612 //If the user has said there is signal only and no noise in the SFTs, then the noise background of the FF plane will be filled with zeros
1613 if ( !params->signalOnly ) {
1614 //Window and psd allocation
1615 REAL4Window *win = NULL;
1616 REAL4VectorAligned *psd = NULL;
1617 XLAL_CHECK( ( win = XLALCreateHannREAL4Window( numffts ) ) != NULL, XLAL_EFUNC ); //Window function
1618 XLAL_CHECK( ( psd = XLALCreateREAL4VectorAligned( numfprbins, 32 ) ) != NULL, XLAL_EFUNC ); //Current PSD calculation
1619 REAL4 winFactor = 8.0 / 3.0;
1620 REAL8 dutyfactor = 0.0, dutyfactorincrement = 1.0 / ( REAL8 )numffts;
1621
1622 //Average each SFT across the frequency band, also compute normalization factor
1623 REAL4VectorAligned *backgroundScalingOverBand = NULL, *aveBackgroundScalingInTime = NULL, *scaledAveNoiseInTime = NULL;
1624 XLAL_CHECK( ( backgroundScalingOverBand = XLALCreateREAL4VectorAligned( numfbins, 32 ) ) != NULL, XLAL_EFUNC );
1625 XLAL_CHECK( ( aveBackgroundScalingInTime = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1626 XLAL_CHECK( ( scaledAveNoiseInTime = XLALCreateREAL4VectorAligned( numffts, 32 ) ) != NULL, XLAL_EFUNC );
1627
1628 memset( aveBackgroundScalingInTime->data, 0, sizeof( REAL4 )*numffts );
1629 for ( UINT4 ii = 0; ii < aveNoiseInTime->length; ii++ ) {
1630 if ( sftexist->data[ii] != 0 ) {
1631 memcpy( backgroundScalingOverBand->data, &( backgroundScaling->data[ii * numfbins] ), sizeof( REAL4 )*backgroundScalingOverBand->length );
1632 aveBackgroundScalingInTime->data[ii] = calcMean( backgroundScalingOverBand );
1634
1635 XLAL_CHECK( XLALVectorMultiplyREAL4( backgroundScalingOverBand->data, backgroundScalingOverBand->data, backgroundScalingOverBand->data, backgroundScalingOverBand->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1636 REAL4 bandMeanBackgroundScalingSq = calcMean( backgroundScalingOverBand );
1638
1639 if ( !params->noiseWeightOff ) {
1640 sumofweights += ( antweights->data[ii] * antweights->data[ii] * bandMeanBackgroundScalingSq ) / ( aveNoiseInTime->data[ii] * aveNoiseInTime->data[ii] );
1641 } else {
1642 sumofweights += ( antweights->data[ii] * antweights->data[ii] * bandMeanBackgroundScalingSq );
1643 }
1644
1645 dutyfactor += dutyfactorincrement;
1646 }
1647 } /* for ii < aveNoiseInTime->length */
1648 invsumofweights = 1.0 / sumofweights;
1649
1650 //Load time series of powers, normalize, mean subtract and Hann window
1651 REAL4VectorAligned *x = NULL, *multiplicativeFactor = NULL;
1652 XLAL_CHECK( ( x = XLALCreateREAL4VectorAligned( aveNoiseInTime->length, 32 ) ) != NULL, XLAL_EFUNC );
1653 XLAL_CHECK( ( multiplicativeFactor = XLALCreateREAL4VectorAligned( aveNoiseInTime->length, 32 ) ) != NULL, XLAL_EFUNC );
1654
1655 memset( multiplicativeFactor->data, 0, sizeof( REAL4 )*multiplicativeFactor->length );
1656 REAL4 psdfactor = winFactor * params->Tobs; //Only multiply by Tobs instead of Tobs/numffts^2 because of the weighting normalization
1657
1658 for ( UINT4 ii = 0; ii < x->length; ii++ ) {
1659 if ( aveNoiseInTime->data[ii] != 0.0 ) {
1660 if ( !params->noiseWeightOff ) {
1661 multiplicativeFactor->data[ii] = win->data->data[ii] * antweights->data[ii] / ( aveNoiseInTime->data[ii] * aveNoiseInTime->data[ii] ) * invsumofweights;
1662 } else {
1663 multiplicativeFactor->data[ii] = win->data->data[ii] * antweights->data[ii] * invsumofweights;
1664 }
1665 }
1666 }
1667
1668 //Multiply the aveNoiseInTime by aveBackgroundScalingInTime
1669 XLAL_CHECK( XLALVectorMultiplyREAL4( scaledAveNoiseInTime->data, aveNoiseInTime->data, aveBackgroundScalingInTime->data, aveNoiseInTime->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1670
1671 //New version. Computes expected background with correlation estimate from Hann windowed and overlapped (must be 50% or 0%) SFTs
1672 REAL8 correlationfactor = 0.0;
1673 for ( INT4 ii = 0; ii < ( INT4 )floor( win->data->length * ( params->SFToverlap / params->Tsft ) - 1 ); ii++ ) {
1674 correlationfactor += win->data->data[ii] * win->data->data[ii + ( INT4 )( ( 1.0 - ( params->SFToverlap / params->Tsft ) ) * win->data->length )];
1675 }
1676 correlationfactor /= win->sumofsquares;
1677 REAL8 corrfactorsquared = correlationfactor * correlationfactor;
1678 //REAL8 prevnoiseval = 0.0, noiseval = 0.0;
1679 for ( UINT4 ii = 0; ii < 1000; ii++ ) {
1680 memset( x->data, 0, sizeof( REAL4 )*x->length );
1681
1682 XLAL_CHECK( sampleREAL4VectorAligned( x, expDistVals, rng ) == XLAL_SUCCESS, XLAL_EFUNC );
1683 XLAL_CHECK( XLALVectorMultiplyREAL4( x->data, x->data, scaledAveNoiseInTime->data, x->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1684
1685 for ( UINT4 jj = 0; jj < x->length; jj++ ) if ( sftexist->data[jj] != 0 ) {
1686 x->data[jj] = x->data[jj] - scaledAveNoiseInTime->data[jj];
1687 }
1688
1689 //Window and rescale because of antenna and noise weights
1690 XLAL_CHECK( XLALVectorMultiplyREAL4( x->data, x->data, multiplicativeFactor->data, x->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1691
1692 //Correlations (forward)
1693 for ( UINT4 jj = 0; jj < x->length; jj++ ) {
1694 if ( sftexist->data[jj] != 0 ) {
1695 if ( jj > 0 && sftexist->data[jj - 1] != 0 ) {
1696 x->data[jj] *= ( 1.0 - corrfactorsquared );
1697 x->data[jj] += corrfactorsquared * x->data[jj - 1];
1698 }
1699 }
1700 }
1701
1702 //Do the FFT
1704
1705 //Sum into the bins
1706 XLAL_CHECK( XLALVectorAddREAL4( aveNoise->data, aveNoise->data, psd->data, psd->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1707 } /* for ii < 1000 */
1708
1709 //Average and rescale
1710 //REAL4 averageRescaleFactor = 2.5e-4*psdfactor*(1.0+2.0*corrfactorsquared);
1711 REAL4 averageRescaleFactor = 1.0e-3 * psdfactor; //*(1.0+2.0*corrfactorsquared);
1712 XLAL_CHECK( XLALVectorScaleREAL4( aveNoise->data, averageRescaleFactor, aveNoise->data, aveNoise->length ) == XLAL_SUCCESS, XLAL_EFUNC );
1713
1714 //Fix 0th and end bins (0 only for odd x->length, 0 and end for even x->length)
1715 if ( GSL_IS_EVEN( x->length ) == 1 ) {
1716 aveNoise->data[0] *= 2.0;
1717 aveNoise->data[aveNoise->length - 1] *= 2.0;
1718 } else {
1719 aveNoise->data[0] *= 2.0;
1720 }
1721
1722 //Compute normalization
1723 *( normalization ) = 1.0 / ( calcMean( aveNoise ) );
1724 for ( UINT4 ii = 0; ii < aveNoise->length; ii++ ) {
1725 aveNoise->data[ii] *= *( normalization );
1726 }
1727
1728 //Extra factor for normalization to 1.0 (empirically determined) (Old method)
1729 //*normalization /= 1.08;
1730
1731 //Extra factor for normalization to 1.0
1732 //duty factor because there are zero-valued SFTs (fourth root term)
1733 //square root term because the uncertainty on the background improves with the size of the running median block
1734 //correlation term because neighboring sfts are correlated
1735 *normalization /= ( 1.0 + pow( dutyfactor, 0.25 ) / sqrt( ( REAL8 )params->blksize ) ) * sqrt( 1.0 + 2.0 * corrfactorsquared );
1736
1737 //Destroy stuff
1741 XLALDestroyREAL4VectorAligned( aveBackgroundScalingInTime );
1742 XLALDestroyREAL4VectorAligned( scaledAveNoiseInTime );
1743 XLALDestroyREAL4VectorAligned( backgroundScalingOverBand );
1744 XLALDestroyREAL4VectorAligned( multiplicativeFactor );
1745 } else {
1746 *( normalization ) = 1.0;
1747 }
1748
1749 //fprintf(stderr, "Noise average = %g ", calcMean((REAL4Vector*)aveNoise));
1750
1751 fprintf( stderr, "done\n" );
1752
1753 return XLAL_SUCCESS;
1754
1755} /* ffPlaneNoise() */
1756
1757/**
1758 * \param [in] IFO Pointer to LALStringVector of the IFO names like H1/L1/V1/etc.
1759 * \return Pointer to MultiLALDetector
1760 */
1762{
1763 XLAL_CHECK_NULL( IFO != NULL, XLAL_EINVAL );
1764 for ( UINT4 ii = 0; ii < IFO->length; ii++ ) {
1765 if ( ii == 0 ) {
1766 fprintf( stderr, "IFO = %s", IFO->data[ii] );
1767 } else {
1768 fprintf( stderr, "%s", IFO->data[ii] );
1769 }
1770 if ( ii < IFO->length - 1 ) {
1771 fprintf( stderr, "," );
1772 } else {
1773 fprintf( stderr, "\n" );
1774 }
1775 }
1777 XLAL_CHECK_NULL( ( detectors = XLALMalloc( sizeof( MultiLALDetector ) ) ) != NULL, XLAL_ENOMEM );
1778 //XLAL_CHECK_NULL( XLALParseMultiLALDetector(detectors, IFO) == XLAL_SUCCESS, XLAL_EFUNC ); //Normally use this but we need the special use case of H2r
1779 detectors->length = IFO->length;
1780 for ( UINT4 ii = 0; ii < detectors->length; ii++ ) {
1781 if ( strcmp( "H1", IFO->data[ii] ) == 0 ) {
1783 } else if ( strcmp( "L1", IFO->data[ii] ) == 0 ) {
1785 } else if ( strcmp( "V1", IFO->data[ii] ) == 0 ) {
1787 } else if ( strcmp( "H2", IFO->data[ii] ) == 0 ) {
1789 } else if ( strcmp( "H2r", IFO->data[ii] ) == 0 ) {
1793 memset( &( H2.frDetector.name ), 0, sizeof( CHAR )*LALNameLength );
1794 snprintf( H2.frDetector.name, LALNameLength, "%s", "LHO_2k_rotatedPiOver4" );
1795 XLAL_CHECK_NULL( ( XLALCreateDetector( &( detectors->sites[ii] ), &( H2.frDetector ), LALDETECTORTYPE_IFODIFF ) ) != NULL, XLAL_EFUNC );
1796 } else {
1797 XLAL_ERROR_NULL( XLAL_EINVAL, "Not using valid interferometer! Expected 'H1', 'H2', 'H2r' (rotated H2), 'L1', or 'V1' not %s.\n", IFO->data[ii] );
1798 }
1799 }
1800 return detectors;
1801}
1802
1803INT4 readTwoSpectInputParams( UserInput_t *uvar, int argc, char *argv[] )
1804{
1805
1806 XLAL_CHECK( uvar != NULL, XLAL_EINVAL, "Invalid NULL input 'uvar'\n" );
1807 XLAL_CHECK( argv != NULL, XLAL_EINVAL, "Invalid NULL input 'argv'\n" );
1808
1809 uvar->ephemEarth = XLALStringDuplicate( "earth00-40-DE405.dat.gz" );
1810 uvar->ephemSun = XLALStringDuplicate( "sun00-40-DE405.dat.gz" );
1811 uvar->blksize = 101;
1812 uvar->outfilename = XLALStringDuplicate( "logfile.txt" );
1813 uvar->configCopy = XLALStringDuplicate( "input_args.conf" );
1814 uvar->ULfilename = XLALStringDuplicate( "uls.dat" );
1815 uvar->candidatesFilename = XLALStringDuplicate( "candidates.dat" );
1816 uvar->harmonicNumToSearch = 1;
1817 uvar->periodHarmToCheck = 5;
1818 uvar->periodFracToCheck = 3;
1819 uvar->ihsfactor = 5;
1820 uvar->minTemplateLength = 1;
1821 uvar->maxTemplateLength = 500;
1822 uvar->FFTplanFlag = 1;
1823 uvar->vectorMath = 0;
1824 uvar->injRandSeed = 0;
1825 uvar->ULsolver = 0;
1826 uvar->dopplerMultiplier = 1.0;
1827 uvar->ihsfar = 1.0;
1828 uvar->tmplfar = 1.0;
1829 uvar->ihsfomfar = 1.0;
1830 uvar->ihsfom = 0.0;
1831 uvar->keepOnlyTopNumIHS = -1;
1832 uvar->lineDetection = -1.0;
1833 uvar->cosiSignCoherent = 0;
1834
1835 XLALRegisterUvarMember( outdirectory, STRING, 0, REQUIRED, "Output directory" );
1836 XLALRegisterUvarMember( IFO, STRINGVector, 0, REQUIRED, "CSV list of detectors, eg. \"H1,H2,L1,G1, ...\" " );
1837 XLALRegisterUvarMember( Tobs, REAL8, 0, REQUIRED, "Total observation time (in seconds)" );
1838 XLALRegisterUvarMember( Tsft, REAL8, 0, REQUIRED, "SFT coherence time (in seconds)" );
1839 XLALRegisterUvarMember( SFToverlap, REAL8, 0, REQUIRED, "SFT overlap (in seconds), usually Tsft/2" );
1840 XLALRegisterUvarMember( t0, REAL8, 0, REQUIRED, "Start time of the search (in GPS seconds)" );
1841 XLALRegisterUvarMember( fmin, REAL8, 0, REQUIRED, "Minimum frequency of band (Hz)" );
1842 XLALRegisterUvarMember( fspan, REAL8, 0, REQUIRED, "Frequency span of band (Hz)" );
1843 XLALRegisterUvarMember( avesqrtSh, STRINGVector, 0, REQUIRED, "CSV list of expected average of square root of Sh for each detector" );
1844 XLALRegisterUvarMember( Pmin, REAL8, 0, REQUIRED, "Minimum period to be searched (in seconds)" );
1845 XLALRegisterUvarMember( Pmax, REAL8, 0, REQUIRED, "Maximum period to be searched (in seconds)" );
1846 XLALRegisterUvarMember( dfmin, REAL8, 0, REQUIRED, "Minimum modulation depth to search (Hz)" );
1847 XLALRegisterUvarMember( dfmax, REAL8, 0, REQUIRED, "Maximum modulation depth to search (Hz)" );
1848 XLALRegisterUvarMember( blksize, INT4, 0, OPTIONAL, "Blocksize for running median to determine expected noise of input SFTs" );
1849 XLALRegisterUvarMember( outfilename, STRING, 0, OPTIONAL, "Output file name" );
1850 XLALRegisterUvarMember( configCopy, STRING, 0, OPTIONAL, "Copy of the input values" );
1851 XLALRegisterUvarMember( ULfilename, STRING, 0, OPTIONAL, "Upper limit file name" );
1852 XLALRegisterUvarMember( candidatesFilename, STRING, 0, OPTIONAL, "Candidates file name" );
1853 XLALRegisterUvarMember( inputSFTs, STRING, 0, OPTIONAL, "Path and filename of SFTs, conflicts with timestampsFile and segmentFile. Possibilities are:\n"
1854 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" );
1855 XLALRegisterUvarMember( ephemEarth, STRING, 0, OPTIONAL, "Earth ephemeris file to use" );
1856 XLALRegisterUvarMember( ephemSun, STRING, 0, OPTIONAL, "Sun ephemeris file to use" );
1857 XLALRegisterUvarMember( skyRegion, STRING, 0, OPTIONAL, "Region of the sky to search (e.g. (ra1,dec1),(ra2,dec2),(ra3,dec3)...) or allsky" );
1858 XLALRegisterUvarMember( skyRegionFile, STRING, 0, OPTIONAL, "File with the grid points" );
1859 XLALRegisterUvarMember( linPolAngle, REAL8, 0, OPTIONAL, "Polarization angle to search using linear polarization of Fplus (when unspecified default is circular polarization)" );
1860 XLALRegisterUvarMember( templatebankfile, STRING, 0, OPTIONAL, "File containing the template data (generated by lalpulsar_TwoSpectTemplateBank)" );
1861 XLALRegisterUvarMember( harmonicNumToSearch, INT4, 0, OPTIONAL, "Number of harmonics of the Pmin to Pmax range to search" );
1862 XLALRegisterUvarMember( periodHarmToCheck, INT4, 0, OPTIONAL, "Number of harmonics/sub-harmonics of the IHS candidates to test" );
1863 XLALRegisterUvarMember( periodFracToCheck, INT4, 0, OPTIONAL, "Number of fractional periods to check in the sense of [(1...N)+1]/[(1...N)+2]" );
1864 XLALRegisterUvarMember( templateSearch, BOOLEAN, 0, OPTIONAL, "Flag for doing a pure template-based search on search region specified by (sky,f,fspan,P, Asini +- 3 AsiniSigma)" );
1865 XLALRegisterUvarMember( templateSearchP, REAL8, 0, OPTIONAL, "The template search period; templateSearch flag is required" );
1866 XLALRegisterUvarMember( templateSearchAsini, REAL8, 0, OPTIONAL, "The template search Asini; templateSearch flag is required" );
1867 XLALRegisterUvarMember( templateSearchAsiniSigma, REAL8, 0, OPTIONAL, "The template search uncertainty in Asini; templateSearch flag is required" );
1868 XLALRegisterUvarMember( templateSearchFixedDf, BOOLEAN, 0, OPTIONAL, "Flag for doing a template-based search on search region specified by (df,sky,f,fspan,P)" );
1869 XLALRegisterUvarMember( templateSearchDf, STRINGVector, 0, OPTIONAL, "The (list of) template search df(s); templateSearchFixedDf flag is required" );
1870 XLALRegisterUvarMember( assumeNScosi, REAL8, 0, OPTIONAL, "Assume cosi orientation of the source (only used when specifying more than 1 detector)" );
1871 XLALRegisterUvarMember( assumeNSpsi, REAL8, 0, OPTIONAL, "Assume psi polarization angle of the source (only used when specifying more than 1 detector)" );
1872 XLALRegisterUvarMember( assumeNSGWfreq, REAL8, 0, OPTIONAL, "Assume GW frequency of the source (only used when specifying more than 1 detector)" );
1873 XLALRegisterUvarMember( assumeNSorbitP, REAL8, 0, OPTIONAL, "Assume orbital period of the source (only used when specifying more than 1 detector)" );
1874 XLALRegisterUvarMember( assumeNSasini, REAL8, 0, OPTIONAL, "Assume projected semi-major axis (units of lt-sec) of the source (only used when specifying more than 1 detector)" );
1875 XLALRegisterUvarMember( assumeNSorbitTp, EPOCH, 0, OPTIONAL, "Assume NS time of periapsis passage of the source (only used when specifying more than 1 detector)" );
1876 XLALRegisterUvarMember( assumeNSrefTime, EPOCH, 0, OPTIONAL, "Assume NS spin reference time (if not given, assume start time of search) (only used when specifying more than 1 detector)" );
1877 XLALRegisterUvarMember( cosiSignCoherent, INT4, 0, OPTIONAL, "For coherent analysis assume [-1,1] values (0), [0,1] values (1), or [-1,0] values (-1) for cosi (Note: unused when assumeNScosi is specified)" );
1878 XLALRegisterUvarMember( ihsfactor, INT4, 0, OPTIONAL, "Number of harmonics to sum in IHS algorithm" );
1879 XLALRegisterUvarMember( ihsfar, REAL8, 0, OPTIONAL, "IHS FAR threshold" );
1880 XLALRegisterUvarMember( ihsfom, REAL8, 0, OPTIONAL, "IHS FOM = 12*(L_IHS_loc - U_IHS_loc)^2" );
1881 XLALRegisterUvarMember( ihsfomfar, REAL8, 0, OPTIONAL, "IHS FOM FAR threshold" );
1882 XLALRegisterUvarMember( tmplfar, REAL8, 0, OPTIONAL, "Template FAR threshold" );
1883 XLALRegisterUvarMember( keepOnlyTopNumIHS, INT4, 0, OPTIONAL, "Keep the top <number> of IHS candidates based on significance" );
1884 XLALRegisterUvarMember( minTemplateLength, INT4, 0, OPTIONAL, "Minimum number of pixels to use in the template" );
1885 XLALRegisterUvarMember( maxTemplateLength, INT4, 0, OPTIONAL, "Maximum number of pixels to use in the template" );
1886 XLALRegisterUvarMember( ULfmin, REAL8, 0, OPTIONAL, "Minimum signal frequency considered for the upper limit value (Hz)" );
1887 XLALRegisterUvarMember( ULfspan, REAL8, 0, OPTIONAL, "Span of signal frequencies considered for the upper limit value (Hz)" );
1888 XLALRegisterUvarMember( ULminimumDeltaf, REAL8, 0, OPTIONAL, "Minimum modulation depth counted in the upper limit value (Hz)" );
1889 XLALRegisterUvarMember( ULmaximumDeltaf, REAL8, 0, OPTIONAL, "Maximum modulation depth counted in the upper limit value (Hz)" );
1890 XLALRegisterUvarMember( allULvalsPerSkyLoc, BOOLEAN, 0, OPTIONAL, "Print all UL values in the band specified by ULminimumDeltaf and ULmaximumDeltaf (default prints only the maximum UL)" );
1891 XLALRegisterUvarMember( markBadSFTs, BOOLEAN, 0, OPTIONAL, "Mark bad SFTs" );
1892 XLALRegisterUvarMember( lineDetection, REAL8, 0, OPTIONAL, "Detect stationary lines above threshold, and, if any present, set upper limit only, no template follow-up" );
1893 XLALRegisterUvarMember( FFTplanFlag, INT4, 0, OPTIONAL, "0=Estimate, 1=Measure, 2=Patient, 3=Exhaustive" );
1894 XLALRegisterUvarMember( fastchisqinv, BOOLEAN, 0, OPTIONAL, "Use a faster central chi-sq inversion function (roughly float precision instead of double)" );
1895 XLALRegisterUvarMember( vectorMath, INT4, 0, OPTIONAL, "Vector math functions: 0=None, 1=SSE, 2=AVX/SSE (Note that user needs to have compiled for SSE or AVX/SSE or program fails)" );
1896 XLALRegisterUvarMember( followUpOutsideULrange, BOOLEAN, 0, OPTIONAL, "Follow up outliers outside the range of the UL values" );
1897 XLALRegisterUvarMember( timestampsFile, STRINGVector, 0, OPTIONAL, "CSV list of files with timestamps, file-format: lines of <GPSsec> <GPSnsec>, conflicts with inputSFTs and segmentFile" );
1898 XLALRegisterUvarMember( segmentFile, STRINGVector, 0, OPTIONAL, "CSV list of files with segments, file-format: lines with <GPSstart> <GPSend>, conflicts with inputSFTs and timestampsFile" );
1899 XLALRegisterUvarMember( gaussNoiseWithSFTgaps, BOOLEAN, 0, OPTIONAL, "Gaussian noise using avesqrtSh with same gaps as inputSFTs, conflicts with timstampsFile and segmentFile" );
1900 XLALRegisterUvarMember( injectionSources, STRINGVector, 0, OPTIONAL, "CSV file list containing sources to inject or '{Alpha=0;Delta=0;...}'" );
1901 XLALRegisterUvarMember( injFmin, REAL8, 0, OPTIONAL, "Minimum frequency of band to create in TwoSpect" );
1902 XLALRegisterUvarMember( injBand, REAL8, 0, OPTIONAL, "Width of band to create in TwoSpect" );
1903 XLALRegisterUvarMember( injRandSeed, INT4, 0, OPTIONAL, "Random seed value for reproducable noise, conflicts with inputSFTs" );
1904 XLALRegisterUvarMember( weightedIHS, BOOLEAN, 0, DEVELOPER, "Use the noise-weighted IHS scheme" );
1905 XLALRegisterUvarMember( signalOnly, BOOLEAN, 0, DEVELOPER, "SFTs contain only signal, no noise" );
1906 XLALRegisterUvarMember( templateTest, BOOLEAN, 0, DEVELOPER, "Test the doubly-Fourier transformed data against a single, exact template" );
1907 XLALRegisterUvarMember( templateTestF, REAL8, 0, DEVELOPER, "The template test frequency; templateTest flag is required" );
1908 XLALRegisterUvarMember( templateTestP, REAL8, 0, DEVELOPER, "The template test period; templateTest flag is required" );
1909 XLALRegisterUvarMember( templateTestDf, REAL8, 0, DEVELOPER, "The template test modulation depth; templateTest flag is required" );
1910 XLALRegisterUvarMember( bruteForceTemplateTest, BOOLEAN, 0, DEVELOPER, "Test a number of different templates using templateTest parameters" );
1911 XLALRegisterUvarMember( ULsolver, INT4, 0, DEVELOPER, "Solver function for the upper limit calculation: 0=gsl_ncx2cdf_float_withouttinyprob_solver, 1=gsl_ncx2cdf_withouttinyprob_solver, 2=gsl_ncx2cdf_float_solver, 3=gsl_ncx2cdf_solver, 4=ncx2cdf_float_withouttinyprob_withmatlabchi2cdf_solver, 5=ncx2cdf_withouttinyprob_withmatlabchi2cdf_solver" );
1912 XLALRegisterUvarMember( dopplerMultiplier, REAL8, 0, DEVELOPER, "Multiplier for the Doppler velocity" );
1913 XLALRegisterUvarMember( IHSonly, BOOLEAN, 0, DEVELOPER, "IHS stage only is run. Output statistic is the IHS statistic" );
1914 XLALRegisterUvarMember( noNotchHarmonics, BOOLEAN, 0, DEVELOPER, "Do not notch the daily/sidereal harmonics in the IHS step" );
1915 XLALRegisterUvarMember( calcRthreshold, BOOLEAN, 0, DEVELOPER, "Calculate the threshold value for R given the template false alarm rate" );
1916 XLALRegisterUvarMember( antennaOff, BOOLEAN, 0, DEVELOPER, "Antenna pattern weights are /NOT/ used if this flag is used" );
1917 XLALRegisterUvarMember( noiseWeightOff, BOOLEAN, 0, DEVELOPER, "Turn off noise weighting if this flag is used" );
1918 XLALRegisterUvarMember( gaussTemplatesOnly, BOOLEAN, 0, DEVELOPER, "Gaussian templates only throughout the pipeline if this flag is used" );
1919 XLALRegisterUvarMember( ULoff, BOOLEAN, 0, DEVELOPER, "Turn off upper limits computation" );
1920 XLALRegisterUvarMember( BrentsMethod, BOOLEAN, 0, DEVELOPER, "Use Brent's method for root finding of the R threshold value" );
1921 XLALRegisterUvarMember( printSFTtimes, BOOLEAN, 0, DEVELOPER, "Output a list <GPS sec> <GPS nanosec> of SFT start times of SFTs" );
1922 XLALRegisterUvarMember( printData, BOOLEAN, 0, DEVELOPER, "Print to ASCII files the data values" );
1923 XLALRegisterUvarMember( saveRvalues, STRING, 0, DEVELOPER, "Print all R values from template bank file" );
1924 XLALRegisterUvarMember( printSignalData, STRING, 0, DEVELOPER, "Print f0 and h0 per SFT of the signal, used only with injectionSources" );
1925 XLALRegisterUvarMember( printMarginalizedSignalData, STRING, 0, DEVELOPER, "Print f0 and h0 per SFT of the signal, used only with injectionSources" );
1926 XLALRegisterUvarMember( randSeed, INT4, 0, DEVELOPER, "Random seed value" );
1927 XLALRegisterUvarMember( chooseSeed, BOOLEAN, 0, DEVELOPER, "The random seed value is chosen based on the input search parameters" );
1928
1929 //Read all the input from config file and command line (command line has priority)
1930 //Also checks required variables unless help is requested
1931 BOOLEAN should_exit = 0;
1933 if ( should_exit ) {
1934 exit( 1 );
1935 }
1936
1937 //Check analysis parameters
1938 if ( ceil( uvar->t0 / uvar->SFToverlap ) * uvar->SFToverlap - uvar->t0 != 0.0 ) {
1939 REAL8 oldstart = uvar->t0;
1940 uvar->t0 = ceil( uvar->t0 / uvar->SFToverlap ) * uvar->SFToverlap;
1941 fprintf( stderr, "WARNING! Adjusting start time from %f to %f\n", oldstart, uvar->t0 );
1942 }
1943 XLAL_CHECK( uvar->Pmax >= uvar->Pmin, XLAL_EINVAL, "Pmax is smaller than Pmin\n" );
1944 XLAL_CHECK( uvar->dfmax >= uvar->dfmin, XLAL_EINVAL, "dfmax is smaller than dfmin\n" );
1945 if ( uvar->Pmax > 0.2 * uvar->Tobs ) {
1946 uvar->Pmax = 0.2 * uvar->Tobs;
1947 fprintf( stderr, "WARNING! Adjusting input maximum period to 1/5 the observation time!\n" );
1948 }
1949 if ( uvar->Pmin < 4.0 * uvar->Tsft ) {
1950 uvar->Pmin = 4.0 * uvar->Tsft;
1951 fprintf( stderr, "WARNING! Adjusting input minimum period to 4 coherence times!\n" );
1952 }
1953 if ( uvar->dfmax > maxModDepth( uvar->Pmax, uvar->Tsft ) ) {
1954 uvar->dfmax = floor( maxModDepth( uvar->Pmax, uvar->Tsft ) * ( uvar->Tsft ) ) / ( uvar->Tsft );
1955 fprintf( stderr, "WARNING! Adjusting input maximum modulation depth due to maximum modulation depth allowed\n" );
1956 }
1957 if ( uvar->dfmin < 0.5 / uvar->Tsft ) {
1958 uvar->dfmin = 0.5 / uvar->Tsft;
1959 fprintf( stderr, "WARNING! Adjusting input minimum modulation depth to 1/2 a frequency bin!\n" );
1960 }
1961 uvar->dfmin = 0.5 * round( 2.0 * uvar->dfmin * uvar->Tsft ) / uvar->Tsft;
1962 uvar->dfmax = 0.5 * round( 2.0 * uvar->dfmax * uvar->Tsft ) / uvar->Tsft;
1963 UINT4 firstbin, numbins;
1964 XLAL_CHECK( XLALFindCoveringSFTBins( &firstbin, &numbins, uvar->fmin - uvar->dfmax - 6.0 / uvar->Tsft, uvar->fspan + 2.0 * uvar->dfmax + 12.0 / uvar->Tsft, uvar->Tsft ) == XLAL_SUCCESS, XLAL_EFUNC );
1965 REAL8 newfmin = ( firstbin + 6 + uvar->dfmax * uvar->Tsft ) / uvar->Tsft;
1966 REAL8 newfspan = ( numbins - 2 * uvar->dfmax * uvar->Tsft - 13 ) / uvar->Tsft;
1967 if ( newfmin != uvar->fmin ) {
1968 uvar->fmin = newfmin;
1969 fprintf( stderr, "WARNING! Adjusting fmin to %g\n", newfmin );
1970 }
1971 if ( newfspan != uvar->fspan ) {
1972 uvar->fspan = newfspan;
1973 fprintf( stderr, "WARNING! Adjusting fspan to %g\n", newfspan );
1974 }
1975
1976 //If specifying coherent addition assumed frequency parameters, then check all are present
1977 XLAL_CHECK( !( XLALUserVarWasSet( &uvar->assumeNSGWfreq ) || XLALUserVarWasSet( &uvar->assumeNSorbitP ) || XLALUserVarWasSet( &uvar->assumeNSasini ) || XLALUserVarWasSet( &uvar->assumeNSorbitTp ) ) || ( XLALUserVarWasSet( &uvar->assumeNSGWfreq ) && XLALUserVarWasSet( &uvar->assumeNSorbitP ) && XLALUserVarWasSet( &uvar->assumeNSasini ) && XLALUserVarWasSet( &uvar->assumeNSorbitTp ) ), XLAL_EINVAL, "Need to specify all parameters for frequency evolution (assumeNSrefTime is optional)" );
1978
1979 //Check blocksize settings for running mean
1980 if ( uvar->blksize % 2 != 1 ) {
1981 uvar->blksize += 1;
1982 }
1983
1984 //Check timestampsFile/inputSFTs/segmentFile conflicts
1985 if ( XLALUserVarWasSet( &uvar->timestampsFile ) && XLALUserVarWasSet( &uvar->inputSFTs ) ) {
1986 XLAL_ERROR( XLAL_FAILURE, "Cannot specify timestampsFile and inputSFTs\n" );
1987 } else if ( XLALUserVarWasSet( &uvar->segmentFile ) && XLALUserVarWasSet( &uvar->inputSFTs ) ) {
1988 XLAL_ERROR( XLAL_FAILURE, "Cannot specify segmentFile and inputSFTs\n" );
1989 } else if ( XLALUserVarWasSet( &uvar->segmentFile ) && XLALUserVarWasSet( &uvar->timestampsFile ) ) {
1990 XLAL_ERROR( XLAL_FAILURE, "Cannot specify segmentFile and timestampsFile\n" );
1991 }
1992
1993 //Check skyRegion and skyRegionFile options
1994 if ( ( XLALUserVarWasSet( &uvar->skyRegion ) && XLALUserVarWasSet( &uvar->skyRegionFile ) ) || ( !XLALUserVarWasSet( &uvar->skyRegion ) && !XLALUserVarWasSet( &uvar->skyRegionFile ) ) ) {
1995 XLAL_ERROR( XLAL_EINVAL, "Specify only one of skyRegion or skyRegionFile\n" );
1996 }
1997
1998 //Check required options if specifying templateSearchFixedDf
1999 if ( uvar->templateSearchFixedDf || XLALUserVarWasSet( &uvar->templateSearchDf ) ) {
2000 if ( !( uvar->templateSearchFixedDf && XLALUserVarWasSet( &uvar->templateSearchDf ) && XLALUserVarWasSet( &uvar->templateSearchP ) ) ) {
2001 XLAL_ERROR( XLAL_FAILURE, "Must specify all or none of templateSearchFixedDf, templateSearchDf, templateSearchP\n" );
2002 }
2003 }
2004
2005 //Check required options if specifying templateSearch
2006 if ( uvar->templateSearch || XLALUserVarWasSet( &uvar->templateSearchP ) || XLALUserVarWasSet( &uvar->templateSearchAsini ) || XLALUserVarWasSet( &uvar->templateSearchAsiniSigma ) ) {
2007 if ( !( uvar->templateSearchFixedDf ) && !( uvar->templateSearch && XLALUserVarWasSet( &uvar->templateSearchP ) && XLALUserVarWasSet( &uvar->templateSearchAsini ) && XLALUserVarWasSet( &uvar->templateSearchAsiniSigma ) ) ) {
2008 XLAL_ERROR( XLAL_FAILURE, "Must specify all or none of templateSearch, templateSearchP, templateSearchAsini, and templateSearchAsiniSigma\n" );
2009 }
2010 }
2011
2012 //Error if dfmax is smaller than templateTestDf or if dfmax is smaller than the templateSearch or bruteForceTemplateTest largest modulation depth
2013 if ( uvar->templateTest ) {
2014 XLAL_CHECK( uvar->dfmax >= uvar->templateTestDf, XLAL_EINVAL, "templateTestDf is larger than dfmax\n" );
2015 }
2016 if ( uvar->templateSearch ) {
2017 XLAL_CHECK( uvar->dfmax >= LAL_TWOPI * ( uvar->fmin + uvar->fspan ) * ( uvar->templateSearchAsini + 3.0 * uvar->templateSearchAsiniSigma ) / uvar->templateSearchP, XLAL_EINVAL, "templateSearch parameters would make the largest modulation depth larger than dfmax\n" );
2018 }
2019 if ( uvar->bruteForceTemplateTest ) {
2020 XLAL_CHECK( uvar->dfmax >= uvar->templateTestDf + 2.0 / uvar->Tsft, XLAL_EINVAL, "templateTestDf+2/Tsft is larger than dfmax\n" );
2021 }
2022
2023 //Check IHS FOM threshold conflicts
2024 if ( XLALUserVarWasSet( &uvar->ihsfomfar ) && XLALUserVarWasSet( &uvar->ihsfom ) ) {
2025 XLAL_ERROR( XLAL_FAILURE, "Only choose only one of ihsfomfar or ihsfom\n" );
2026 }
2027
2028 //Check ranges of settings for FAR values
2029 if ( uvar->ihsfar < 0.0 || uvar->ihsfar > 1.0 ) {
2030 XLAL_ERROR( XLAL_EINVAL, "ihsfar must lie between 0 and 1 (inclusive)\n" );
2031 }
2032 if ( uvar->ihsfomfar < 0.0 || uvar->ihsfomfar > 1.0 ) {
2033 XLAL_ERROR( XLAL_EINVAL, "ihsfomfar must lie between 0 and 1 (inclusive)\n" );
2034 }
2035 if ( uvar->tmplfar < 0.0 || uvar->tmplfar > 1.0 ) {
2036 XLAL_ERROR( XLAL_EINVAL, "tmplfar must lie between 0 and 1 (inclusive)\n" );
2037 }
2038
2039 //params->log10templatefar = log10(params->templatefar);
2040
2041 //Upper limit settings take span of search values unless specified
2042 if ( !XLALUserVarWasSet( &uvar->ULfmin ) ) {
2043 uvar->ULfmin = uvar->fmin; //Upper limit minimum frequency (Hz)
2044 }
2045 if ( !XLALUserVarWasSet( &uvar->ULfspan ) ) {
2046 uvar->ULfspan = uvar->fspan; //Upper limit frequency span (Hz)
2047 }
2048 if ( !XLALUserVarWasSet( &uvar->ULminimumDeltaf ) ) {
2049 uvar->ULminimumDeltaf = uvar->dfmin; //Upper limit minimum modulation depth (Hz)
2050 }
2051 if ( !XLALUserVarWasSet( &uvar->ULmaximumDeltaf ) ) {
2052 uvar->ULmaximumDeltaf = uvar->dfmax; //Upper limit maximum modulation depth (Hz)
2053 }
2054
2055 //Check SSE/AVX settings
2056 if ( uvar->vectorMath > 2 || uvar->vectorMath < 0 ) {
2057 XLAL_ERROR( XLAL_FAILURE, "Must specify vectorMath to be 0, 1, or 2" );
2058 }
2059
2060 //Developer options
2061 if ( uvar->templateTest && uvar->bruteForceTemplateTest ) {
2062 XLAL_ERROR( XLAL_FAILURE, "Specify one of templateTest or bruteForceTemplateTest\n" );
2063 }
2064 if ( ( uvar->templateTest || uvar->bruteForceTemplateTest ) || XLALUserVarWasSet( &uvar->templateTestF ) || XLALUserVarWasSet( &uvar->templateTestP ) || XLALUserVarWasSet( &uvar->templateTestDf ) ) {
2065 if ( !( ( uvar->templateTest || uvar->bruteForceTemplateTest ) && XLALUserVarWasSet( &uvar->templateTestF ) && XLALUserVarWasSet( &uvar->templateTestP ) && XLALUserVarWasSet( &uvar->templateTestDf ) ) ) {
2066 XLAL_ERROR( XLAL_FAILURE, "Must specify all or none of templateTest/bruteForceTemplateTest, templateTestF, templateTestP, and templateTestDf\n" );
2067 }
2068 }
2069
2070 //SFT input conflicts with injRandSeed option
2071 if ( XLALUserVarWasSet( &uvar->injRandSeed ) && XLALUserVarWasSet( &uvar->inputSFTs ) && !uvar->gaussNoiseWithSFTgaps ) {
2072 XLAL_ERROR( XLAL_EINVAL, "When specifying injRandSeed, inputSFTs cannot be used unless specifying gaussNoiseWithSFTgaps\n" );
2073 }
2074
2075 return XLAL_SUCCESS;
2076
2077}
2078
2079
2080/**
2081 * Print REAL4Vector values to an ASCII file
2082 * \param [in] vector Pointer to the REAL4Vector to be printed
2083 * \param [in] directory CHAR array of the file path
2084 * \param [in] filename CHAR array of the file name
2085 * \return Status value
2086 */
2088{
2089
2090 CHARVector *outputfile = NULL;
2091 XLAL_CHECK( ( outputfile = XLALCreateCHARVector( strlen( directory ) + strlen( filename ) + 3 ) ) != NULL, XLAL_EFUNC );
2092 sprintf( outputfile->data, "%s/%s", directory, filename );
2093 FILE *OUTPUT = NULL;
2094 XLAL_CHECK( ( OUTPUT = fopen( outputfile->data, "w" ) ) != NULL, XLAL_EIO, "Couldn't fopen %s/%s for writing.\n", directory, filename );
2095 for ( UINT4 ii = 0; ii < vector->length; ii++ ) {
2096 fprintf( OUTPUT, "%g\n", vector->data[ii] );
2097 }
2098 fclose( OUTPUT );
2099 XLALDestroyCHARVector( outputfile );
2100
2101 return XLAL_SUCCESS;
2102
2103}
#define __func__
log an I/O error, i.e.
void InitDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan, const DopplerSkyScanInit *init)
Definition: DopplerScan.c:273
void FreeDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan)
Definition: DopplerScan.c:324
int XLALNextDopplerSkyPos(PulsarDopplerParams *pos, DopplerSkyScanState *skyScan)
NextDopplerSkyPos(): step through sky-grid return 0 = OK, -1 = ERROR.
Definition: DopplerScan.c:127
@ STATE_FINISHED
all templates have been read
Definition: DopplerScan.h:84
@ GRID_FILE_SKYGRID
read skygrid from a file
Definition: DopplerScan.h:94
@ GRID_ISOTROPIC
approximately isotropic sky-grid
Definition: DopplerScan.h:92
#define IFO
void destroyihsfarStruct(ihsfarStruct *ihsfarstruct)
Destroy ihsfarStruct struct.
Definition: IHS.c:316
void destroyihsMaxima(ihsMaximaStruct *data)
Destroy vectors and the IHS maxima struct.
Definition: IHS.c:65
INT4 runIHS(ihsMaximaStruct *output, const ffdataStruct *input, const ihsfarStruct *ihsfarinput, const UserInput_t *params, const UINT4 rows, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *FbinMean)
Run the IHS algorithm.
Definition: IHS.c:89
ihsfarStruct * createihsfarStruct(const UINT4 rows, const UserInput_t *params)
Allocate memory for ihsfarStruct struct.
Definition: IHS.c:289
INT4 genIhsFar(ihsfarStruct *output, const UserInput_t *params, const UINT4 rows, const REAL4VectorAligned *aveNoise, const gsl_rng *rng)
Compute the IHS FAR for a sum of a number of rows.
Definition: IHS.c:340
INT4 findIHScandidates(candidateVector **candlist, const ihsfarStruct *ihsfarstruct, const UserInput_t *params, const ffdataStruct *ffdata, const ihsMaximaStruct *ihsmaxima, const REAL4VectorAligned *fbinavgs, const REAL4VectorSequence *trackedlines)
Find IHS candidates above thresholds.
Definition: IHS.c:809
ihsMaximaStruct * createihsMaxima(const UINT4 fbins, const UINT4 rows)
Create vectors for IHS maxima struct.
Definition: IHS.c:37
void LALCheckMemoryLeaks(void)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
#define STRING(a)
INT4 tfRngMeans(REAL4VectorAligned *output, const REAL4VectorAligned *tfdata, const UINT4 numffts, const UINT4 numfbins, const UINT4 blksize)
Determine the running mean of each SFT.
LIGOTimeGPSVector * jointTimestampsFromMultiTimestamps(const MultiLIGOTimeGPSVector *multiTimestamps)
Definition: SFTfunctions.c:981
INT4 replaceTFdataWithSubsequentTFdata(REAL4VectorAlignedArray *tfdataarray, const UINT4 numffts)
INT4Vector * markBadSFTs(const REAL4VectorAligned *tfdata, const UserInput_t *params)
Mark the non-Gaussian SFTs using K-S and Kuiper's tests.
INT4 tfMeanSubtract(REAL4VectorAligned *tfdata, const REAL4VectorAligned *rngMeans, const REAL4VectorAligned *backgroundScaling, const UINT4 numffts, const UINT4 numfbins)
Subtract running mean values from the SFT data, modifying input time-frequency data.
INT4 printSFTtimestamps2File(const MultiSFTVector *multiSFTvector, const UserInput_t *params)
REAL4VectorAligned * convertSFTdataToPowers(const SFTVector *sfts, const UserInput_t *params, const REAL8 normalization)
Convert a SFTVector of sfts into powers.
Definition: SFTfunctions.c:753
INT4 tfWeight(REAL4VectorAligned *output, const REAL4VectorAligned *tfdata, REAL4VectorAligned *rngMeans, REAL4VectorAligned *antPatternWeights, const REAL4VectorAligned *backgroundScaling, const INT4Vector *indexValuesOfExistingSFTs, const UserInput_t *params)
Weight the SFTs based on antenna pattern and noise variance (Equation 11, assuming the input time-fre...
MultiSFTVector * generateSFTdata(UserInput_t *uvar, const MultiLALDetector *detectors, const EphemerisData *edat, const INT4 maxbinshift, const gsl_rng *rng)
INT4 slideTFdata(REAL4VectorAligned *output, const UserInput_t *params, const REAL4VectorAligned *tfdata, const INT4Vector *binshifts)
Slide the time-frequency data to account for detector motion.
MultiSFTVector * getMultiSFTVector(UserInput_t *params, const REAL8 minfbin, const REAL8 maxfbin)
Get a MultiSFTVector from the user-input values.
Definition: SFTfunctions.c:102
INT4 checkBackgroundScaling(const REAL4VectorAligned *background, const REAL4VectorAligned *backgroundScaling, const INT4Vector *SFTexistVector)
Go through the backgroundScaling vector and zero out if the SFTexistVector has a 0 in that SFT locati...
INT4Vector * existingSFTs(const REAL4VectorAligned *tfdata, const UINT4 numffts)
Determine if the SFTs are existing based on whether the first data point of power in each SFT is 0.
REAL4VectorAligned * coherentlyAddSFTs(const MultiSFTVector *multiSFTvector, const MultiSSBtimes *multissb, const MultiAMCoeffs *multiAMcoeffs, const LIGOTimeGPSVector *jointTimestamps, const REAL4VectorAlignedArray *backgroundRatio, const INT4 cosiSign, const assumeNSparams *NSparams, const UserInput_t *params, REAL4VectorAligned *backgroundScaling)
Add SFTs together from a MultiSFTVector.
Definition: SFTfunctions.c:147
#define fprintf
int main(int argc, char *argv[])
Main program.
Definition: TwoSpect.c:60
REAL4VectorAligned * calcAveTFnoisePerFbinRatio(const REAL4VectorAligned *background, const REAL4VectorAligned *backgroundScaling, const UINT4 numffts)
Calculate the ratio of the average SFT noise to the mean of the average SFT noise.
Definition: TwoSpect.c:1398
INT4 makeSecondFFT(ffdataStruct *output, REAL4VectorAligned *tfdata, const REAL4FFTPlan *plan)
Compute the second Fourier transform for TwoSpect.
Definition: TwoSpect.c:1435
REAL4 rmsTFdataBand(const REAL4VectorAligned *backgrnd, UINT4 numfbins, UINT4 numffts, UINT4 binmin, UINT4 binmax)
Determine the rms of the noise power in each frequency bin across the band.
Definition: TwoSpect.c:1541
void destroyffdata(ffdataStruct *data)
Free the frequency-frequency data structure.
Definition: TwoSpect.c:1209
INT4 medianBackgroundBandInTime(REAL4VectorAligned *aveNoiseInTime, const REAL4VectorAligned *backgrnd, const INT4Vector *sftexist)
Definition: TwoSpect.c:1566
REAL4VectorSequence * trackLines(const INT4Vector *lines, const INT4Vector *binshifts, const REAL4 minfbin, const REAL4 df)
Track the lines for the sky position.
Definition: TwoSpect.c:1318
MultiLALDetector * setupMultiLALDetector(LALStringVector *IFO)
Definition: TwoSpect.c:1761
INT4 cleanLines(REAL4VectorAligned *TFdata, const REAL4VectorAligned *background, const INT4Vector *lines, const UserInput_t *params, const gsl_rng *rng)
Test algorithm to clean lines.
Definition: TwoSpect.c:1359
INT4 ffPlaneNoise(REAL4VectorAligned *aveNoise, const UserInput_t *params, const INT4Vector *sftexist, const REAL4VectorAligned *aveNoiseInTime, const REAL4VectorAligned *antweights, const REAL4VectorAligned *backgroundScaling, const REAL4FFTPlan *plan, const REAL4VectorAligned *expDistVals, const gsl_rng *rng, REAL8 *normalization)
Measure of the average noise power in each 2nd FFT frequency bin.
Definition: TwoSpect.c:1596
INT4Vector * detectLines_simple(const REAL4VectorAligned *TFdata, const ffdataStruct *ffdata, const UserInput_t *params)
Line detection algorithm.
Definition: TwoSpect.c:1231
INT4 printREAL4Vector2File(const REAL4Vector *vector, const CHAR *directory, const CHAR *filename)
Print REAL4Vector values to an ASCII file.
Definition: TwoSpect.c:2087
REAL4 avgTFdataBand(const REAL4VectorAligned *backgrnd, UINT4 numfbins, UINT4 numffts, UINT4 binmin, UINT4 binmax)
Determine the average of the noise power in each frequency bin across the band.
Definition: TwoSpect.c:1507
ffdataStruct * createffdata(const UserInput_t *params)
Create a new frequency-frequency data structure for the TwoSpect analysis.
Definition: TwoSpect.c:1185
FILE * LOG
Definition: TwoSpect.c:55
INT4 readTwoSpectInputParams(UserInput_t *uvar, int argc, char *argv[])
Definition: TwoSpect.c:1803
INT4 CompBinShifts(INT4Vector *output, const SSBtimes *ssbTimes, const REAL8 freq, const REAL8 Tsft, const REAL4 dopplerMultiplier)
Compute the number of integer bin shifts per SFT.
Definition: antenna.c:38
INT4 CompAntennaPatternWeights2(REAL4VectorAligned *output, const SkyPosition skypos, const LIGOTimeGPSVector *timestamps, const LALDetector det, const REAL8 *cosi, const REAL8 *psi)
Definition: antenna.c:93
REAL4 CompDetectorVmax2(const LIGOTimeGPSVector *timestamps, const LALDetector det, EphemerisData *edat)
Definition: antenna.c:256
INT4 calcHarmonicMean(REAL4 *harmonicMean, const REAL4VectorAligned *vector, const UINT4 numfbins, const UINT4 numffts)
Compute the harmonic mean value of a REAL4VectorAligned of SFT values.
INT4 calcMedian(REAL4 *median, const REAL4VectorAligned *vector)
Calculate the median value from a REAL4VectorAligned.
INT4 min_max_index_INT4Vector(const INT4Vector *inputvector, UINT4 *min_index_out, UINT4 *max_index_out)
Determine the index value of the maximum and minimum values in an INT4Vector.
INT4 sampleREAL4VectorAligned(REAL4VectorAligned *output, const REAL4VectorAligned *input, const gsl_rng *rng)
Sample a number (sampleSize) of values from a REAL4VectorAligned (input) randomly.
INT4 calcRms(REAL4 *rms, const REAL4VectorAligned *vector)
Compute the RMS value of a REAL4VectorAligned.
REAL8 expRandNum(const REAL8 mu, const gsl_rng *ptrToGenerator)
Create a exponentially distributed noise value.
INT4 calcStddev(REAL4 *sigma, const REAL4VectorAligned *vector)
Compute the standard deviation of a REAL4VectorAligned.
REAL4VectorAligned * expRandNumVector(const UINT4 length, const REAL8 mu, const gsl_rng *ptrToGenerator)
REAL4 calcMean(const REAL4VectorAligned *vector)
Compute the mean value of a REAL4VectorAligned, computed via recursion like in GSL.
INT4 templateSearch_scox1Style(candidateVector **output, const REAL8 fminimum, const REAL8 fspan, const REAL8 period, const REAL8 asini, const REAL8 asinisigma, const SkyPosition skypos, const UserInput_t *params, const REAL4VectorAligned *ffdata, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *aveTFnoisePerFbinRatio, const REAL4VectorSequence *trackedlines, const REAL4FFTPlan *secondFFTplan, const gsl_rng *rng, BOOLEAN useExactTemplates)
A brute force template search to find the most significant template around a putative source whose pa...
Definition: candidates.c:575
candidateVector * keepMostSignificantCandidates(const candidateVector *input, const UserInput_t *params)
Keep the most significant candidates, potentially reducing the number of candidates if there are more...
Definition: candidates.c:1297
candidateVector * createcandidateVector(const UINT4 length)
Allocate a candidateVector.
Definition: candidates.c:32
INT4 testIHScandidates(candidateVector **output, const candidateVector *ihsCandidates, const ffdataStruct *ffdata, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *aveTFnoisePerFbinRatio, const SkyPosition pos, const UserInput_t *params, const gsl_rng *rng)
Function to test the IHS candidates against Gaussian templates.
Definition: candidates.c:1062
INT4 analyzeCandidatesTemplateFromVector(candidateVector *output, const candidateVector *input, const TwoSpectTemplateVector *vector, const ffdataStruct *ffdata, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *aveTFnoisePerFbinRatio, const UserInput_t *params, const gsl_rng *rng, const UINT4 templateLen)
Definition: candidates.c:191
candidateVector * resizecandidateVector(candidateVector *vector, const UINT4 length)
Resize a candidateVector.
Definition: candidates.c:61
INT4 writeCandidateVector2File(const CHAR *outputfile, const candidateVector *input)
Definition: candidates.c:1385
INT4 bruteForceTemplateSearch(candidate *output, const candidate input, const TwoSpectParamSpaceSearchVals *paramspace, const UserInput_t *params, const REAL4VectorAligned *ffdata, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *aveTFnoisePerFbinRatio, const REAL4FFTPlan *secondFFTplan, const gsl_rng *rng, const BOOLEAN useExactTemplates)
A brute force template search to find the most significant template around a candidate.
Definition: candidates.c:254
REAL8 maxModDepth(const REAL8 period, const REAL8 cohtime)
Calculates maximum modulation depth allowed, equation 6 of E.
Definition: candidates.c:1414
INT4 bruteForceTemplateTest(candidateVector **output, const candidate input, const TwoSpectParamSpaceSearchVals *paramspace, const UserInput_t *params, const REAL4VectorAligned *ffdata, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *aveTFnoisePerFbinRatio, const REAL4FFTPlan *secondFFTplan, const gsl_rng *rng, const BOOLEAN useExactTemplates)
A brute force template search to test templates around a candidate.
Definition: candidates.c:429
INT4 testTwoSpectTemplateVector(candidateVector *output, const TwoSpectTemplateVector *templateVec, const ffdataStruct *ffdata, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *aveTFnoisePerFbinRatio, const SkyPosition skypos, const UserInput_t *params, const gsl_rng *rng, const UINT4 templateLen)
Test each of the templates in a TwoSpectTemplateVector and keep the top 10 This will not check the fa...
Definition: candidates.c:1222
INT4 clusterCandidates(candidateVector **output, const candidateVector *input, const ffdataStruct *ffdata, const UserInput_t *params, const REAL4VectorAligned *ffplanenoise, const REAL4VectorAligned *fbinaveratios, const gsl_rng *rng, const BOOLEAN exactflag)
Cluster candidates by frequency, period, and modulation depth using templates.
Definition: candidates.c:837
INT4 analyzeOneTemplate(candidate *output, const candidate *input, const ffdataStruct *ffdata, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *aveTFnoisePerFbinRatio, const UserInput_t *params, const REAL4FFTPlan *plan, const gsl_rng *rng, const BOOLEAN exactflag)
Analyze a single template.
Definition: candidates.c:153
void destroycandidateVector(candidateVector *vector)
Free a candidateVector.
Definition: candidates.c:89
void loadCandidateData(candidate *output, const REAL8 fsig, const REAL8 period, const REAL8 moddepth, const REAL4 ra, const REAL4 dec, const REAL8 statval, const REAL8 h0, const REAL8 prob, const INT4 proberrcode, const REAL8 normalization, const INT4 templateVectorIndex, const BOOLEAN lineContamination)
Load candidate data.
Definition: candidates.c:122
INT4 templateSearch_fixedDf(candidateVector **output, const LALStringVector *dffixed, const REAL8 fminimum, const REAL8 fspan, const REAL8 period, const SkyPosition skypos, const UserInput_t *params, const REAL4VectorAligned *ffdata, const REAL4VectorAligned *aveNoise, const REAL4VectorAligned *aveTFnoisePerFbinRatio, const REAL4VectorSequence *trackedlines, const REAL4FFTPlan *secondFFTplan, const gsl_rng *rng, BOOLEAN useExactTemplates)
A brute force template search to find the most significant template at a fixed modulation depth aroun...
Definition: candidates.c:724
void XLALDestroyREAL4VectorSequence(REAL4VectorSequence *vecseq)
REAL4VectorSequence * XLALCreateREAL4VectorSequence(UINT4 length, UINT4 veclen)
LALDetector * XLALCreateDetector(LALDetector *detector, const LALFrDetector *frDetector, LALDetectorType type)
const LALDetector lalCachedDetectors[LAL_NUM_DETECTORS]
void XLALDestroyMultiDetectorStateSeries(MultiDetectorStateSeries *mdetStates)
Helper function to get rid of a multi-IFO DetectorStateSeries Note, this is "NULL-robust" in the sens...
MultiDetectorStateSeries * XLALGetMultiDetectorStates(const MultiLIGOTimeGPSVector *multiTS, const MultiLALDetector *multiIFO, const EphemerisData *edat, REAL8 tOffset)
Get the detector-time series for the given MultiLIGOTimeGPSVector.
int XLALParseMultiNoiseFloor(MultiNoiseFloor *multiNoiseFloor, const LALStringVector *sqrtSX, UINT4 numDetectors)
Parse string-vectors (typically input by user) of N detector noise-floors for detectors ,...
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
void XLALDestroyMultiAMCoeffs(MultiAMCoeffs *multiAMcoef)
Destroy a MultiAMCoeffs structure.
Definition: LALComputeAM.c:469
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
Definition: LALComputeAM.c:379
#define LAL_PI
#define LAL_TWOPI
unsigned char BOOLEAN
uint64_t UINT8
double REAL8
#define XLAL_INIT_DECL(var,...)
char CHAR
uint32_t UINT4
int32_t INT4
float REAL4
LALNameLength
LALDETECTORTYPE_IFODIFF
LAL_LLO_4K_DETECTOR
LAL_VIRGO_DETECTOR
LAL_LHO_2K_DETECTOR
LAL_LHO_4K_DETECTOR
void * XLALMalloc(size_t n)
void XLALFree(void *p)
void LALSRunningMedian2(LALStatus *status, REAL4Sequence *medians, const REAL4Sequence *input, LALRunningMedianPar param)
char char * XLALStringDuplicate(const char *s)
char * XLALVCSInfoString(const LALVCSInfoList vcs_list, const int verbose, const char *prefix)
REAL4FFTPlan * XLALCreateForwardREAL4FFTPlan(UINT4 size, int measurelvl)
void XLALDestroyREAL4FFTPlan(REAL4FFTPlan *plan)
int XLALREAL4PowerSpectrum(REAL4Vector *_LAL_RESTRICT_ spec, const REAL4Vector *_LAL_RESTRICT_ data, const REAL4FFTPlan *plan)
void XLALDestroyMultiSFTVector(MultiSFTVector *multvect)
Destroy a multi SFT-vector.
Definition: SFTtypes.c:424
MultiLIGOTimeGPSVector * XLALMakeMultiTimestamps(LIGOTimeGPS tStart, REAL8 Tspan, REAL8 Tsft, REAL8 Toverlap, UINT4 numDet)
Same as XLALMakeTimestamps() just for several detectors, additionally specify the number of detectors...
int XLALFindCoveringSFTBins(UINT4 *firstBin, UINT4 *numBins, REAL8 fMinIn, REAL8 BandIn, REAL8 Tsft)
Return the 'effective' frequency-band [fMinEff, fMaxEff] = [firstBin, lastBin] * 1/Tsft,...
Definition: SFTtypes.c:106
MultiLIGOTimeGPSVector * XLALExtractMultiTimestampsFromSFTs(const MultiSFTVector *multiSFTs)
Given a multi-SFT vector, return a MultiLIGOTimeGPSVector holding the SFT timestamps.
void XLALDestroyMultiTimestamps(MultiLIGOTimeGPSVector *multiTS)
Destroy a MultiLIGOTimeGPSVector timestamps vector.
void XLALDestroyTimestampVector(LIGOTimeGPSVector *vect)
De-allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:69
MultiSSBtimes * XLALGetMultiSSBtimes(const MultiDetectorStateSeries *multiDetStates, SkyPosition skypos, LIGOTimeGPS refTime, SSBprecision precision)
Multi-IFO version of XLALGetSSBtimes().
Definition: SSBtimes.c:657
void XLALDestroyMultiSSBtimes(MultiSSBtimes *multiSSB)
Destroy a MultiSSBtimes structure.
Definition: SSBtimes.c:845
@ SSBPREC_RELATIVISTICOPT
optimized relativistic, numerically equivalent to SSBPREC_RELATIVISTIC, but faster
Definition: SSBtimes.h:48
COORDINATESYSTEM_EQUATORIAL
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
void XLALDestroyUserVars(void)
#define XLALRegisterUvarMember(name, type, option, category,...)
void CHAR * XLALUserVarGetLog(UserVarLogFormat format)
int XLALUserVarWasSet(const void *cvar)
UVAR_LOGFMT_CFGFILE
INT4Vector * XLALResizeINT4Vector(INT4Vector *vector, UINT4 length)
INT4Vector * XLALCreateINT4Vector(UINT4 length)
void XLALDestroyINT4Vector(INT4Vector *vector)
void XLALDestroyCHARVector(CHARVector *vector)
CHARVector * XLALCreateCHARVector(UINT4 length)
int XLALVectorScaleREAL4(REAL4 *out, REAL4 scalar, const REAL4 *in, const UINT4 len)
void XLALDestroyREAL4VectorAligned(REAL4VectorAligned *in)
REAL4VectorAligned * XLALCreateREAL4VectorAligned(const UINT4 length, const UINT4 align)
int XLALVectorAddREAL4(REAL4 *out, const REAL4 *in1, const REAL4 *in2, const UINT4 len)
int XLALVectorMultiplyREAL4(REAL4 *out, const REAL4 *in1, const REAL4 *in2, const UINT4 len)
void XLALDestroyREAL4Window(REAL4Window *window)
REAL4Window * XLALCreateHannREAL4Window(UINT4 length)
#define XLAL_ERROR_NULL(...)
#define xlalErrno
#define XLAL_ERROR(...)
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_REAL4(assertion,...)
#define XLAL_CHECK_NULL(assertion,...)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EIO
XLAL_EINVAL
XLAL_FAILURE
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
float data[BLOCKSIZE]
Definition: hwinject.c:360
char * directory
Definition: hwinject.c:363
psd
list weights
LALDetector detectors[LAL_NUM_DETECTORS]
double t0
CHAR * data
initialization-structure passed to InitDopplerSkyScan()
Definition: DopplerScan.h:134
this structure reflects the current state of a DopplerSkyScan
Definition: DopplerScan.h:152
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
INT4 * data
UINT4 length
LALFrDetector frDetector
REAL4 xArmAzimuthRadians
REAL4 yArmAzimuthRadians
CHAR name[LALNameLength]
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
UINT4 length
number of timestamps
Definition: SFTfileIO.h:192
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
Definition: LALComputeAM.h:137
Multi-IFO time-series of DetectorStates.
array of detectors definitions 'LALDetector'
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 sqrtSn[PULSAR_MAX_DETECTORS]
per-IFO sqrt(PSD) values , where
UINT4 length
number of detectors
A collection of SFT vectors – one for each IFO in a multi-IFO search.
Definition: SFTfileIO.h:179
SFTVector ** data
sftvector for each ifo
Definition: SFTfileIO.h:184
Multi-IFO container for SSB timings.
Definition: SSBtimes.h:67
SSBtimes ** data
array of SSBtimes (pointers)
Definition: SSBtimes.h:72
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
REAL8 Delta
Sky position: DEC (latitude) in equatorial coords and radians.
REAL8 Alpha
Sky position: RA (longitude) in equatorial coords and radians.
REAL4VectorAligned ** data
REAL4 * data
REAL8 * data
REAL8Vector * ULval
REAL8 normalization
UpperLimit * data
INT4 templateVectorIndex
REAL8 moddepth
BOOLEAN lineContamination
REAL8 normalization
INT4 proberrcode
REAL8 period
candidate * data
REAL8 ffnormalization
REAL4VectorAligned * ffdata
REAL8 tfnormalization
REAL4VectorAligned * ihsfar
void destroyTwoSpectTemplateVector(TwoSpectTemplateVector *vector)
Free a TwoSpectTemplateVector.
Definition: templates.c:120
TwoSpectTemplateVector * readTwoSpectTemplateVector(const CHAR *filename)
Read a TwoSpectTemplateVector from a binary file.
Definition: templates.c:250
void destroyUpperLimitVector(UpperLimitVector *vector)
Free an UpperLimitVector.
Definition: upperlimits.c:87
INT4 skypoint95UL(UpperLimit *ul, const UserInput_t *params, const ffdataStruct *ffdata, const ihsMaximaStruct *ihsmaxima, const ihsfarStruct *ihsfar, const REAL4VectorAligned *fbinavgs)
Determine the 95% confidence level upper limit at a particular sky location from the loudest IHS valu...
Definition: upperlimits.c:155
UpperLimitVector * resizeUpperLimitVector(UpperLimitVector *vector, const UINT4 length)
Resize an UpperLimitVector.
Definition: upperlimits.c:59
UpperLimitVector * createUpperLimitVector(const UINT4 length)
Allocate a new UpperLimitVector.
Definition: upperlimits.c:32
INT4 outputUpperLimitToFile(const CHAR *outputfile, const UpperLimit ul, const BOOLEAN printAllULvalues)
Output the highest upper limit to a file unless printAllULvalues==1 in which case,...
Definition: upperlimits.c:364
REAL4VectorAlignedArray * createREAL4VectorAlignedArray(const UINT4 length, const UINT4 vectorLength, const size_t align)
Definition: vectormath.c:99
void destroyREAL4VectorAlignedArray(REAL4VectorAlignedArray *array)
Definition: vectormath.c:110
double df