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>
40#include <gsl/gsl_math.h>
60int main(
int argc,
char *argv[] )
64 time_t programstarttime, programendtime;
67 time( &programstarttime );
68 ptm = localtime( &programstarttime );
71 gsl_set_error_handler_off();
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 ) ;
82 CHARVector *LOGFILENAME = NULL, *ULFILENAME = NULL, *CANDFILENAME = NULL;
84 sprintf( LOGFILENAME->
data,
"%s/%s", uvar.outdirectory, uvar.outfilename );
86 sprintf( ULFILENAME->data,
"%s/%s", uvar.outdirectory, uvar.ULfilename );
88 sprintf( CANDFILENAME->data,
"%s/%s", uvar.outdirectory, uvar.candidatesFilename );
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;
98 fprintf( INPUTVALS,
"%s", cfgFileStringCopy );
110 fprintf( stderr,
"%s\n", VCSInfoString );
114 fprintf( stderr,
"Program executed on %s", asctime( ptm ) );
115 fprintf(
LOG,
"Program executed on %s", asctime( ptm ) );
118 fprintf( stderr,
"Output directory: %s\n", uvar.outdirectory );
119 fprintf(
LOG,
"Output directory: %s\n", uvar.outdirectory );
131 if ( uvar.chooseSeed ) {
132 UINT8 IFOmultiplier = 0;
133 if ( strcmp( uvar.IFO->data[0],
"H1" ) == 0 ) {
135 }
else if ( strcmp( uvar.IFO->data[0],
"L1" ) == 0 ) {
140 uvar.randSeed = IFOmultiplier * (
UINT8 )fabs( round( uvar.fmin + uvar.fspan + uvar.Pmin + uvar.Pmax + uvar.dfmin + uvar.dfmax ) );
142 gsl_rng_set( rng, uvar.randSeed );
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 );
179 if ( detectorVmax > Vmax ) {
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;
199 scanInit.skyRegionString = uvar.skyRegion;
200 scanInit.numSkyPartitions = 1;
201 scanInit.Freq = uvar.fmin + 0.5 * uvar.fspan;
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 );
208 scanInit.skyGridFile = uvar.skyRegionFile;
209 scanInit.numSkyPartitions = 1;
210 scanInit.Freq = uvar.fmin + 0.5 * uvar.fspan;
211 fprintf(
LOG,
"Sky file = %s\n", uvar.skyRegionFile );
212 fprintf( stderr,
"Sky file = %s\n", uvar.skyRegionFile );
223 REAL8 tempfspan = uvar.fspan + 2.0 * uvar.dfmax + ( uvar.blksize - 1 + 12 ) / uvar.Tsft;
224 INT4 tempnumfbins = (
INT4 )round( tempfspan * uvar.Tsft ) + 1;
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;
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 );
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__ );
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 ) );
275 if ( uvar.printSFTtimes ) {
290 memset( backgroundRatioMeans->data[ii]->data, 0,
sizeof(
REAL4 )*backgroundRatioMeans->data[ii]->length );
296 if ( !uvar.signalOnly ) {
314 if ( !uvar.signalOnly ) {
319 memset( background0->data, 0,
sizeof(
REAL4 )*background->
length );
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];
329 backgroundRatioMeans->data[ii]->data[jj] =
calcMean( oneSFTbackgroundRatio );
332 backgroundRatioMeans->data[ii]->data[jj] = 0.0;
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 ) ) );
346 memset( &( background->
data[( ffdata->
numfbins + 2 * ( maxbinshift - 10 ) )*ii] ), 0,
sizeof(
REAL4 ) * ( ffdata->
numfbins + 2 * ( maxbinshift - 10 ) ) );
351 memcpy( background0->data, background->
data,
sizeof(
REAL4 )*background->
length );
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;
368 INT4 heavilyContaminatedBand = 0;
372 if ( lines != NULL ) {
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" );
381 if ( heavilyContaminatedBand ) {
397 REAL4 frac_tobs_complete = 1.0;
398 if ( !uvar.signalOnly &&
detectors->length == 1 ) {
400 INT4 totalincludedsftnumber = 0;
401 for ( ii = 0; ii < (
INT4 )sftexist->
length; ii++ )
if ( sftexist->
data[ii] == 1 ) {
402 totalincludedsftnumber++;
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__ );
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 ) );
422 for ( ii = 0; ii < (
INT4 )sftexist->
length; ii++ ) {
423 if ( sftexist->
data[ii] == 1 ) {
424 indexValuesOfExistingSFTs->
data[jj] = ii;
431 REAL8 backgroundmeannormfactor = 1.0;
432 if ( !uvar.signalOnly &&
detectors->length == 1 ) {
433 REAL4 harmonicMean = 0.0;
435 backgroundmeannormfactor = 1.0 / harmonicMean;
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 ) );
454 if ( uvar.printData &&
detectors->length == 1 ) {
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 );
477 candidateVector *gaussCandidates1 = NULL, *gaussCandidates2 = NULL, *gaussCandidates3 = NULL, *gaussCandidates4 = NULL, *exactCandidates1 = NULL, *exactCandidates2 = NULL, *ihsCandidates = NULL;
489 REAL4FFTPlan *secondFFTplan = NULL;
508 REAL4 antweightsrms = 0.0;
532 memcpy( backgroundForihs2h0->
data, background0->data,
sizeof(
REAL4 )*background->
length );
533 memcpy( backgroundScalingForihs2h0->data, backgroundScaling->
data,
sizeof(
REAL4 )*backgroundScaling->
length );
542 NSparams.assumeNSpos = skypos0;
544 NSparams.assumeNScosi = &uvar.assumeNScosi;
547 NSparams.assumeNSpsi = &uvar.assumeNSpsi;
550 NSparams.assumeNSGWfreq = &uvar.assumeNSGWfreq;
551 NSparams.assumeNSorbitP = &uvar.assumeNSorbitP;
552 NSparams.assumeNSasini = &uvar.assumeNSasini;
553 NSparams.assumeNSorbitTp = &uvar.assumeNSorbitTp;
555 NSparams.assumeNSrefTime = &uvar.assumeNSrefTime;
558 XLAL_CHECK( ( tmpTFdata =
coherentlyAddSFTs( multiSFTvector, multissb, multiAMcoefficients, jointSFTtimestamps, backgroundRatioMeans, uvar.cosiSignCoherent, &NSparams, &uvar, backgroundScalingForihs2h0 ) ) != NULL,
XLAL_EFUNC );
560 INT4 totalincludedsftnumber = 0;
561 for ( ii = 0; ii < (
INT4 )sftexistForihs2h0->
length; ii++ )
if ( sftexistForihs2h0->
data[ii] == 1 ) {
562 totalincludedsftnumber++;
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__ );
587 REAL8 skypointffnormalization = 1.0;
588 XLAL_CHECK(
ffPlaneNoise( aveNoise, &uvar, sftexistForihs2h0, aveNoiseInTimeForihs2h0, antweightsforihs2h0, backgroundScalingForihs2h0, secondFFTplan, expRandVals, rng, &( skypointffnormalization ) ) ==
XLAL_SUCCESS,
XLAL_EFUNC );
591 fprintf(
LOG,
"Starting TwoSpect analysis...\n" );
592 fprintf( stderr,
"Starting TwoSpect analysis...\n" );
597 fprintf( stderr,
"Sky location: RA = %g, DEC = %g\n", dopplerpos.
Alpha, dopplerpos.
Delta );
605 if ( uvar.antennaOff )
for ( ii = 0; ii < (
INT4 )antweights->
length; ii++ ) {
606 antweights->
data[ii] = 1.0;
636 NSparams.assumeNSpos = skypos;
638 NSparams.assumeNScosi = &uvar.assumeNScosi;
641 NSparams.assumeNSpsi = &uvar.assumeNSpsi;
644 NSparams.assumeNSGWfreq = &uvar.assumeNSGWfreq;
645 NSparams.assumeNSorbitP = &uvar.assumeNSorbitP;
646 NSparams.assumeNSasini = &uvar.assumeNSasini;
647 NSparams.assumeNSorbitTp = &uvar.assumeNSorbitTp;
649 NSparams.assumeNSrefTime = &uvar.assumeNSrefTime;
652 XLAL_CHECK( ( tfdata =
coherentlyAddSFTs( multiSFTvector, multissb, multiAMcoefficients, jointSFTtimestamps, backgroundRatioMeans, uvar.cosiSignCoherent, &NSparams, &uvar, backgroundScaling ) ) != NULL,
XLAL_EFUNC );
662 if ( lines != NULL ) {
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" );
670 if ( heavilyContaminatedBand ) {
675 memcpy( background->
data, background0->data,
sizeof(
REAL4 )*background->
length );
677 if ( !uvar.signalOnly ) {
680 INT4 totalincludedsftnumber = 0;
681 for ( ii = 0; ii < (
INT4 )sftexist->
length; ii++ )
if ( sftexist->
data[ii] == 1 ) {
682 totalincludedsftnumber++;
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__ );
696 for ( ii = 0; ii < (
INT4 )sftexist->
length; ii++ ) {
697 if ( sftexist->
data[ii] == 1 ) {
698 indexValuesOfExistingSFTs->
data[jj] = ii;
703 REAL4 harmonicMean = 0.0;
705 backgroundmeannormfactor = 1.0 / harmonicMean;
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 ) );
718 if ( uvar.printData ) {
732 REAL4 fbin0 = (
REAL4 )( round( uvar.fmin * uvar.Tsft - uvar.dfmax * uvar.Tsft - 6.0 - 0.5 * ( uvar.blksize - 1 ) - (
REAL8 )( maxbinshift ) ) / uvar.Tsft );
735 if ( lines != NULL ) {
741 REAL4 currentAntWeightsRMS = 0.0;
745 REAL4VectorAligned *TFdata_slided = NULL, *background_slided = NULL, *backgroundScaling_slided = NULL;
758 if ( uvar.printData ) {
764 if ( antweightsrms == 0.0 ) {
765 antweightsrms = currentAntWeightsRMS;
767 if ( fabs( currentAntWeightsRMS - antweightsrms ) / antweightsrms >= 0.01 ) {
769 antweightsrms = currentAntWeightsRMS;
777 XLAL_CHECK(
ffPlaneNoise( aveNoise, &uvar, sftexist, aveNoiseInTime, antweights, backgroundScaling_slided, secondFFTplan, expRandVals, rng, &( ffdata->
ffnormalization ) ) ==
XLAL_SUCCESS,
XLAL_EFUNC );
778 if ( uvar.printData ) {
790 if ( uvar.printData ) {
811 REAL4 secFFTmean = 0.0, secFFTsigma = 0.0;
818 if ( uvar.printData ) {
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 );
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" );
829 if ( uvar.templateTest ) {
830 if ( uvar.printData ) {
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 );
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 );
840 if ( exactCandidates2->numofcandidates == exactCandidates2->length - 1 ) {
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 )++;
849 exactCandidates2->data[exactCandidates2->numofcandidates - 1].h0 /= sqrt( ffdata->
tfnormalization ) * pow( frac_tobs_complete * ffdata->
ffnormalization / skypointffnormalization, 0.25 );
851 }
else if ( uvar.bruteForceTemplateTest ) {
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
857 XLAL_CHECK(
bruteForceTemplateTest( &( exactCandidates2 ), cand, ¶mspace, &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 );
864 if ( uvar.templateSearch ) {
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 );
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 );
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 );
888 for ( ii = 0; ii < (
INT4 )subsetVec->length; ii++ ) {
889 if ( exactCandidates2->numofcandidates == exactCandidates2->length ) {
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 );
894 ( exactCandidates2->numofcandidates )++;
901 if ( uvar.signalOnly ) {
908 if ( !
XLALUserVarWasSet( &uvar.templatebankfile ) && !uvar.templateTest && !uvar.templateSearch && !uvar.bruteForceTemplateTest && !uvar.templateSearchFixedDf ) {
921 if (
XLALUserVarWasSet( &uvar.keepOnlyTopNumIHS ) && (
INT4 )ihsCandidates->numofcandidates > uvar.keepOnlyTopNumIHS ) {
924 ihsCandidates->numofcandidates = 0;
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 )++;
936 if ( uvar.IHSonly && !uvar.templateTest && !uvar.templateSearch && !uvar.templateSearchFixedDf && !
XLALUserVarWasSet( &uvar.templatebankfile ) ) {
938 if ( exactCandidates2->length < exactCandidates2->numofcandidates + ihsCandidates->numofcandidates ) {
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 );
947 ( exactCandidates2->numofcandidates )++;
950 }
else if ( !uvar.templateTest && !uvar.templateSearch && !uvar.templateSearchFixedDf && !
XLALUserVarWasSet( &uvar.templatebankfile ) ) {
962 ihsCandidates->numofcandidates = 0;
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 );
977 for ( ii = 0; ii < (
INT4 )gaussCandidates2->numofcandidates; ii++ ) {
979 if ( gaussCandidates3->numofcandidates == gaussCandidates3->length - 1 ) {
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
988 gaussCandidates3->numofcandidates++;
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 );
997 gaussCandidates2->numofcandidates = 0;
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 );
1007 gaussCandidates3->numofcandidates = 0;
1010 for ( ii = 0; ii < (
INT4 )gaussCandidates4->numofcandidates; ii++ ) {
1015 if ( cand.
prob < log10( uvar.tmplfar ) ) {
1016 if ( exactCandidates1->numofcandidates == exactCandidates1->length - 1 ) {
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++;
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 );
1031 gaussCandidates4->numofcandidates = 0;
1034 for ( ii = 0; ii < (
INT4 )exactCandidates1->numofcandidates; ii++ ) {
1036 if ( exactCandidates2->numofcandidates == exactCandidates2->length - 1 ) {
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
1043 if ( !uvar.gaussTemplatesOnly ) {
1044 XLAL_CHECK(
bruteForceTemplateSearch( &( exactCandidates2->data[exactCandidates2->numofcandidates] ), exactCandidates1->data[ii], ¶mspace, &uvar, ffdata->
ffdata, aveNoise, aveTFnoisePerFbinRatio, secondFFTplan, rng, 1 ) ==
XLAL_SUCCESS,
XLAL_EFUNC );
1046 XLAL_CHECK(
bruteForceTemplateSearch( &( exactCandidates2->data[exactCandidates2->numofcandidates] ), exactCandidates1->data[ii], ¶mspace, &uvar, ffdata->
ffdata, aveNoise, aveTFnoisePerFbinRatio, secondFFTplan, rng, 0 ) ==
XLAL_SUCCESS,
XLAL_EFUNC );
1048 exactCandidates2->data[exactCandidates2->numofcandidates].h0 /= sqrt( ffdata->
tfnormalization ) * pow( frac_tobs_complete * ffdata->
ffnormalization / skypointffnormalization, 0.25 );
1049 exactCandidates2->numofcandidates++;
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 );
1056 exactCandidates1->numofcandidates = 0;
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 );
1064 if ( !uvar.ULoff && !uvar.templateTest && !uvar.templateSearch && !uvar.templateSearchFixedDf && !
XLALUserVarWasSet( &uvar.templatebankfile ) ) {
1080 if ( trackedlines != NULL ) {
1083 if (
detectors->length > 1 && !uvar.signalOnly ) {
1093 if ( exactCandidates2->numofcandidates != 0 ) {
1094 fprintf(
LOG,
"\n**Report of candidates:**\n" );
1095 fprintf( stderr,
"\n**Report of candidates:**\n" );
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 );
1104 if ( !uvar.ULoff ) {
1105 for ( ii = 0; ii < (
INT4 )upperlimits->
length - 1; ii++ ) {
1118 if ( !uvar.signalOnly ) {
1127 gsl_rng_free( rng );
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 ) );
1236 fprintf( stderr,
"Running line detection algorithm... " );
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 ) );
1260 weights->data[ii] = 1.0 / ( stddev * stddev );
1261 sumweights +=
weights->data[ii];
1264 REAL4 invsumweights = 1.0 / sumweights;
1268 REAL4VectorAligned *testRMSvals = NULL, *testRngMedian = NULL, *testTSofPowers = NULL;
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;
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 ) {
1289 lines->data[numlines] = ii + ( blksize - 1 ) / 2;
1302 if ( lines != NULL ) {
1303 fprintf( stderr,
"WARNING: %d line(s) found.\n", lines->length );
1323 fprintf( stderr,
"Tracking lines... " );
1331 UINT4 maxshiftindex = 0, minshiftindex = 0;
1333 INT4 maxshift = binshifts->
data[maxshiftindex], minshift = binshifts->
data[minshiftindex];
1339 output->data[ii * 3 + 1] = ( lines->
data[ii] + ( minshift - 1 ) ) *
df + minfbin;
1340 output->data[ii * 3 + 2] = ( lines->
data[ii] + ( maxshift + 1 ) ) *
df + minfbin;
1362 if ( lines == NULL ) {
1369 REAL8 prevnoiseval = 0.0;
1370 for (
INT4 ii = 0; ii < numffts; ii++ ) {
1371 if ( TFdata->
data[ii * numfbins] != 0.0 ) {
1373 if ( lines->
data[jj] - (
params->blksize - 1 ) / 2 >= 0 ) {
1376 if ( ii > 0 && TFdata->
data[( ii - 1 ) * numfbins] != 0.0 ) {
1377 noiseval *= ( 1.0 - 0.167 * 0.167 );
1378 noiseval += 0.167 * prevnoiseval;
1380 TFdata->
data[ii * numfbins + lines->
data[jj]] = noiseval;
1381 prevnoiseval = noiseval;
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] );
1414 aveTFnoisePerFbinRatio->
data[ii] = 0.0;
1415 for (
UINT4 jj = 0; jj < numffts; jj++ ) {
1416 aveTFnoisePerFbinRatio->
data[ii] += TSofPowers->data[jj];
1418 aveTFnoisePerFbinRatio->
data[ii] /= totalweightval;
1420 REAL4 aveTFaveinv = 1.0 /
calcMean( aveTFnoisePerFbinRatio );
1421 for (
UINT4 ii = 0; ii < numfbins; ii++ ) {
1422 aveTFnoisePerFbinRatio->
data[ii] *= aveTFaveinv;
1425 return aveTFnoisePerFbinRatio;
1440 fprintf( stderr,
"Computing second FFT over SFTs... " );
1442 REAL8 winFactor = 8.0 / 3.0;
1451 memcpy( windowData->data,
win->data->data,
sizeof(
REAL4 )*windowData->length );
1453 for (
INT4 ii = 0; ii <
output->numfbins; ii++ ) {
1458 for (
UINT4 jj = 0; jj <
x->length; jj++ ) {
1459 x->data[jj] = tfdata->
data[ii + jj *
output->numfbins];
1467 if ( GSL_IS_EVEN(
x->length ) == 1 ) {
1468 psd->data[0] *= 2.0;
1469 psd->data[
psd->length - 1] *= 2.0;
1471 psd->data[0] *= 2.0;
1479 for (
UINT4 jj = 0; jj <
psd->length; jj++ ) {
1480 output->ffdata->data[
psd->length * ii + jj] = (
REAL4 )(
psd->data[jj] * winFactor *
output->ffnormalization );
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 );
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 );
1555 REAL4 rmsTFdata = 0.0;
1569 memset( aveNoiseInTime->
data, 0,
sizeof(
REAL4 )*aveNoiseInTime->
length );
1572 for (
UINT4 ii = 0; ii < sftexist->
length; ii++ ) {
1573 if ( sftexist->
data[ii] != 0 ) {
1599 XLAL_CHECK( aveNoise != NULL &&
params != NULL && sftexist != NULL && aveNoiseInTime != NULL && antweights != NULL && backgroundScaling != NULL && plan != NULL && normalization != NULL,
XLAL_EINVAL );
1601 fprintf( stderr,
"Computing noise background estimate... " );
1603 REAL8 invsumofweights = 0.0, sumofweights = 0.0;
1606 UINT4 numfbins = backgroundScaling->
length / numffts;
1613 if ( !
params->signalOnly ) {
1619 REAL4 winFactor = 8.0 / 3.0;
1620 REAL8 dutyfactor = 0.0, dutyfactorincrement = 1.0 / (
REAL8 )numffts;
1623 REAL4VectorAligned *backgroundScalingOverBand = NULL, *aveBackgroundScalingInTime = NULL, *scaledAveNoiseInTime = NULL;
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 );
1636 REAL4 bandMeanBackgroundScalingSq =
calcMean( backgroundScalingOverBand );
1639 if ( !
params->noiseWeightOff ) {
1640 sumofweights += ( antweights->
data[ii] * antweights->
data[ii] * bandMeanBackgroundScalingSq ) / ( aveNoiseInTime->
data[ii] * aveNoiseInTime->
data[ii] );
1642 sumofweights += ( antweights->
data[ii] * antweights->
data[ii] * bandMeanBackgroundScalingSq );
1645 dutyfactor += dutyfactorincrement;
1648 invsumofweights = 1.0 / sumofweights;
1655 memset( multiplicativeFactor->data, 0,
sizeof(
REAL4 )*multiplicativeFactor->length );
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;
1663 multiplicativeFactor->data[ii] =
win->data->data[ii] * antweights->
data[ii] * invsumofweights;
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 )];
1676 correlationfactor /=
win->sumofsquares;
1677 REAL8 corrfactorsquared = correlationfactor * correlationfactor;
1679 for (
UINT4 ii = 0; ii < 1000; ii++ ) {
1680 memset(
x->data, 0,
sizeof(
REAL4 )*
x->length );
1685 for (
UINT4 jj = 0; jj <
x->length; jj++ )
if ( sftexist->
data[jj] != 0 ) {
1686 x->data[jj] =
x->data[jj] - scaledAveNoiseInTime->data[jj];
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];
1711 REAL4 averageRescaleFactor = 1.0e-3 * psdfactor;
1715 if ( GSL_IS_EVEN(
x->length ) == 1 ) {
1716 aveNoise->
data[0] *= 2.0;
1719 aveNoise->
data[0] *= 2.0;
1723 *( normalization ) = 1.0 / (
calcMean( aveNoise ) );
1724 for (
UINT4 ii = 0; ii < aveNoise->
length; ii++ ) {
1725 aveNoise->
data[ii] *= *( normalization );
1735 *normalization /= ( 1.0 + pow( dutyfactor, 0.25 ) / sqrt( (
REAL8 )
params->blksize ) ) * sqrt( 1.0 + 2.0 * corrfactorsquared );
1746 *( normalization ) = 1.0;
1764 for (
UINT4 ii = 0; ii <
IFO->length; ii++ ) {
1766 fprintf( stderr,
"IFO = %s",
IFO->data[ii] );
1770 if ( ii < IFO->length - 1 ) {
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 ) {
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] );
1811 uvar->blksize = 101;
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;
1826 uvar->dopplerMultiplier = 1.0;
1828 uvar->tmplfar = 1.0;
1829 uvar->ihsfomfar = 1.0;
1831 uvar->keepOnlyTopNumIHS = -1;
1832 uvar->lineDetection = -1.0;
1833 uvar->cosiSignCoherent = 0;
1843 XLALRegisterUvarMember( avesqrtSh, STRINGVector, 0, REQUIRED,
"CSV list of expected average of square root of Sh for each detector" );
1848 XLALRegisterUvarMember( blksize,
INT4, 0, OPTIONAL,
"Blocksize for running median to determine expected noise of input SFTs" );
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>'" );
1857 XLALRegisterUvarMember( skyRegion,
STRING, 0, OPTIONAL,
"Region of the sky to search (e.g. (ra1,dec1),(ra2,dec2),(ra3,dec3)...) or allsky" );
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)" );
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)" );
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)" );
1883 XLALRegisterUvarMember( keepOnlyTopNumIHS,
INT4, 0, OPTIONAL,
"Keep the top <number> of IHS candidates based on significance" );
1890 XLALRegisterUvarMember( allULvalsPerSkyLoc,
BOOLEAN, 0, OPTIONAL,
"Print all UL values in the band specified by ULminimumDeltaf and ULmaximumDeltaf (default prints only the maximum UL)" );
1892 XLALRegisterUvarMember( lineDetection,
REAL8, 0, OPTIONAL,
"Detect stationary lines above threshold, and, if any present, set upper limit only, no template follow-up" );
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)" );
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" );
1903 XLALRegisterUvarMember( injRandSeed,
INT4, 0, OPTIONAL,
"Random seed value for reproducable noise, conflicts with inputSFTs" );
1909 XLALRegisterUvarMember( templateTestDf,
REAL8, 0, DEVELOPER,
"The template test modulation depth; templateTest flag is required" );
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" );
1918 XLALRegisterUvarMember( gaussTemplatesOnly,
BOOLEAN, 0, DEVELOPER,
"Gaussian templates only throughout the pipeline if this flag is used" );
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" );
1933 if ( should_exit ) {
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 );
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" );
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" );
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" );
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" );
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;
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 );
1971 if ( newfspan != uvar->fspan ) {
1972 uvar->fspan = newfspan;
1973 fprintf( stderr,
"WARNING! Adjusting fspan to %g\n", newfspan );
1980 if ( uvar->blksize % 2 != 1 ) {
1999 if ( uvar->templateSearchFixedDf ||
XLALUserVarWasSet( &uvar->templateSearchDf ) ) {
2001 XLAL_ERROR(
XLAL_FAILURE,
"Must specify all or none of templateSearchFixedDf, templateSearchDf, templateSearchP\n" );
2008 XLAL_ERROR(
XLAL_FAILURE,
"Must specify all or none of templateSearch, templateSearchP, templateSearchAsini, and templateSearchAsiniSigma\n" );
2013 if ( uvar->templateTest ) {
2014 XLAL_CHECK( uvar->dfmax >= uvar->templateTestDf,
XLAL_EINVAL,
"templateTestDf is larger than dfmax\n" );
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" );
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" );
2029 if ( uvar->ihsfar < 0.0 || uvar->ihsfar > 1.0 ) {
2032 if ( uvar->ihsfomfar < 0.0 || uvar->ihsfomfar > 1.0 ) {
2035 if ( uvar->tmplfar < 0.0 || uvar->tmplfar > 1.0 ) {
2043 uvar->ULfmin = uvar->fmin;
2046 uvar->ULfspan = uvar->fspan;
2049 uvar->ULminimumDeltaf = uvar->dfmin;
2052 uvar->ULmaximumDeltaf = uvar->dfmax;
2056 if ( uvar->vectorMath > 2 || uvar->vectorMath < 0 ) {
2061 if ( uvar->templateTest && uvar->bruteForceTemplateTest ) {
2066 XLAL_ERROR(
XLAL_FAILURE,
"Must specify all or none of templateTest/bruteForceTemplateTest, templateTestF, templateTestP, and templateTestDf\n" );
2072 XLAL_ERROR(
XLAL_EINVAL,
"When specifying injRandSeed, inputSFTs cannot be used unless specifying gaussNoiseWithSFTgaps\n" );
2093 FILE *OUTPUT = NULL;
#define __func__
log an I/O error, i.e.
void InitDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan, const DopplerSkyScanInit *init)
void FreeDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan)
int XLALNextDopplerSkyPos(PulsarDopplerParams *pos, DopplerSkyScanState *skyScan)
NextDopplerSkyPos(): step through sky-grid return 0 = OK, -1 = ERROR.
@ STATE_FINISHED
all templates have been read
@ GRID_FILE_SKYGRID
read skygrid from a file
@ GRID_ISOTROPIC
approximately isotropic sky-grid
void destroyihsfarStruct(ihsfarStruct *ihsfarstruct)
Destroy ihsfarStruct struct.
void destroyihsMaxima(ihsMaximaStruct *data)
Destroy vectors and the IHS maxima struct.
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.
ihsfarStruct * createihsfarStruct(const UINT4 rows, const UserInput_t *params)
Allocate memory for ihsfarStruct struct.
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.
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.
ihsMaximaStruct * createihsMaxima(const UINT4 fbins, const UINT4 rows)
Create vectors for IHS maxima struct.
void LALCheckMemoryLeaks(void)
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
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)
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.
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.
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.
int main(int argc, char *argv[])
Main program.
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.
INT4 makeSecondFFT(ffdataStruct *output, REAL4VectorAligned *tfdata, const REAL4FFTPlan *plan)
Compute the second Fourier transform for TwoSpect.
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.
void destroyffdata(ffdataStruct *data)
Free the frequency-frequency data structure.
INT4 medianBackgroundBandInTime(REAL4VectorAligned *aveNoiseInTime, const REAL4VectorAligned *backgrnd, const INT4Vector *sftexist)
REAL4VectorSequence * trackLines(const INT4Vector *lines, const INT4Vector *binshifts, const REAL4 minfbin, const REAL4 df)
Track the lines for the sky position.
MultiLALDetector * setupMultiLALDetector(LALStringVector *IFO)
INT4 cleanLines(REAL4VectorAligned *TFdata, const REAL4VectorAligned *background, const INT4Vector *lines, const UserInput_t *params, const gsl_rng *rng)
Test algorithm to clean lines.
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.
INT4Vector * detectLines_simple(const REAL4VectorAligned *TFdata, const ffdataStruct *ffdata, const UserInput_t *params)
Line detection algorithm.
INT4 printREAL4Vector2File(const REAL4Vector *vector, const CHAR *directory, const CHAR *filename)
Print REAL4Vector values to an ASCII file.
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.
ffdataStruct * createffdata(const UserInput_t *params)
Create a new frequency-frequency data structure for the TwoSpect analysis.
INT4 readTwoSpectInputParams(UserInput_t *uvar, int argc, char *argv[])
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.
INT4 CompAntennaPatternWeights2(REAL4VectorAligned *output, const SkyPosition skypos, const LIGOTimeGPSVector *timestamps, const LALDetector det, const REAL8 *cosi, const REAL8 *psi)
REAL4 CompDetectorVmax2(const LIGOTimeGPSVector *timestamps, const LALDetector det, EphemerisData *edat)
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...
candidateVector * keepMostSignificantCandidates(const candidateVector *input, const UserInput_t *params)
Keep the most significant candidates, potentially reducing the number of candidates if there are more...
candidateVector * createcandidateVector(const UINT4 length)
Allocate a candidateVector.
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.
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)
candidateVector * resizecandidateVector(candidateVector *vector, const UINT4 length)
Resize a candidateVector.
INT4 writeCandidateVector2File(const CHAR *outputfile, const candidateVector *input)
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.
REAL8 maxModDepth(const REAL8 period, const REAL8 cohtime)
Calculates maximum modulation depth allowed, equation 6 of E.
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.
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...
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.
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.
void destroycandidateVector(candidateVector *vector)
Free a candidateVector.
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.
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...
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.
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
#define XLAL_INIT_DECL(var,...)
void * XLALMalloc(size_t n)
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.
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,...
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.
MultiSSBtimes * XLALGetMultiSSBtimes(const MultiDetectorStateSeries *multiDetStates, SkyPosition skypos, LIGOTimeGPS refTime, SSBprecision precision)
Multi-IFO version of XLALGetSSBtimes().
void XLALDestroyMultiSSBtimes(MultiSSBtimes *multiSSB)
Destroy a MultiSSBtimes structure.
@ SSBPREC_RELATIVISTICOPT
optimized relativistic, numerically equivalent to SSBPREC_RELATIVISTIC, but faster
COORDINATESYSTEM_EQUATORIAL
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 XLAL_CHECK(assertion,...)
#define XLAL_CHECK_REAL4(assertion,...)
#define XLAL_CHECK_NULL(assertion,...)
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
LALDetector detectors[LAL_NUM_DETECTORS]
initialization-structure passed to InitDopplerSkyScan()
this structure reflects the current state of a DopplerSkyScan
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
A vector of 'timestamps' of type LIGOTimeGPS.
UINT4 length
number of timestamps
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
Multi-IFO time-series of DetectorStates.
array of detectors definitions 'LALDetector'
A collection of (multi-IFO) LIGOTimeGPSVector time-stamps vectors.
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.
SFTVector ** data
sftvector for each ifo
Multi-IFO container for SSB timings.
SSBtimes ** data
array of SSBtimes (pointers)
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
BOOLEAN lineContamination
REAL4VectorAligned * ffdata
REAL4VectorAligned * ihsfar
void destroyTwoSpectTemplateVector(TwoSpectTemplateVector *vector)
Free a TwoSpectTemplateVector.
TwoSpectTemplateVector * readTwoSpectTemplateVector(const CHAR *filename)
Read a TwoSpectTemplateVector from a binary file.
void destroyUpperLimitVector(UpperLimitVector *vector)
Free an UpperLimitVector.
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...
UpperLimitVector * resizeUpperLimitVector(UpperLimitVector *vector, const UINT4 length)
Resize an UpperLimitVector.
UpperLimitVector * createUpperLimitVector(const UINT4 length)
Allocate a new UpperLimitVector.
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,...
REAL4VectorAlignedArray * createREAL4VectorAlignedArray(const UINT4 length, const UINT4 vectorLength, const size_t align)
void destroyREAL4VectorAlignedArray(REAL4VectorAlignedArray *array)