Loading [MathJax]/jax/output/HTML-CSS/config.js
LALPulsar 7.1.1.1-ea7c608
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
candidates.c
Go to the documentation of this file.
1/*
2* Copyright (C) 2010, 2011, 2013, 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#include <sys/stat.h>
21#include <lal/UserInput.h>
22#include "TwoSpect.h"
23#include "candidates.h"
24#include "falsealarm.h"
25#include "templates.h"
26
27/**
28 * Allocate a candidateVector
29 * \param [in] length Length of the candidateVector
30 * \return Pointer to the allocated candidateVector
31 */
32candidateVector *createcandidateVector( const UINT4 length )
33{
34
35 candidateVector *vector = NULL;
36 XLAL_CHECK_NULL( ( vector = XLALMalloc( sizeof( *vector ) ) ) != NULL, XLAL_ENOMEM );
37
38 vector->length = length;
39 vector->numofcandidates = 0;
40 if ( length == 0 ) {
41 vector->data = NULL;
42 } else {
43 XLAL_CHECK_NULL( ( vector->data = XLALMalloc( length * sizeof( *vector->data ) ) ) != NULL, XLAL_ENOMEM );
44 }
45
46 for ( UINT4 ii = 0; ii < vector->length; ii++ ) {
47 vector->data[ii].prob = 0.0;
48 }
49
50 return vector;
51
52} /* createcandidateVector() */
53
54
55/**
56 * Resize a candidateVector
57 * \param [in] vector Pointer of vector to resize
58 * \param [in] length New length of candidateVector
59 * \return Pointer to resized vector
60 */
61candidateVector *resizecandidateVector( candidateVector *vector, const UINT4 length )
62{
63
64 if ( vector == NULL ) {
65 return createcandidateVector( length );
66 }
67 if ( length == 0 ) {
68 destroycandidateVector( vector );
69 return NULL;
70 }
71
72 XLAL_CHECK_NULL( ( vector->data = XLALRealloc( vector->data, length * sizeof( *vector->data ) ) ) != NULL, XLAL_ENOMEM );
73
74 if ( length > vector->length ) for ( UINT4 ii = vector->length; ii < length; ii++ ) {
75 vector->data[ii].prob = 0.0;
76 }
77
78 vector->length = length;
79
80 return vector;
81
82} /* resizecandidateVector() */
83
84
85/**
86 * Free a candidateVector
87 * \param [in] vector Pointer of candidateVector to be freed
88 */
89void destroycandidateVector( candidateVector *vector )
90{
91 if ( vector == NULL ) {
92 return;
93 }
94 if ( ( !vector->length || !vector->data ) && ( vector->length || vector->data ) ) {
95 XLAL_ERROR_VOID( XLAL_EINVAL );
96 }
97 if ( vector->data ) {
98 XLALFree( ( candidate * )vector->data );
99 }
100 vector->data = NULL;
101 XLALFree( ( candidateVector * )vector );
102 return;
103} /* destroycandidateVector() */
104
105
106/**
107 * Load candidate data
108 * \param [out] output Pointer to candidate
109 * \param [in] fsig Frequency of candidate
110 * \param [in] period Orbital period of candidate
111 * \param [in] moddepth Modulation depth of candidate
112 * \param [in] ra Right ascension of candidate
113 * \param [in] dec Declination of candidate
114 * \param [in] statval Detection statistic
115 * \param [in] h0 Estimated strain amplitude
116 * \param [in] prob False alarm probability
117 * \param [in] proberrcode Davies' method error code
118 * \param [in] normalization Time-frequency normalization
119 * \param [in] templateVectorIndex Index value of the template in a templateVector (can be -1 if not from a vector)
120 * \param [in] lineContamination Boolean flag to indicate 0 = no contamination from lines or 1 = likely contaminated by one or more lines
121 */
122void 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 )
123{
124 XLAL_CHECK_VOID( output != NULL, XLAL_EINVAL );
125 output->fsig = fsig;
126 output->period = period;
127 output->moddepth = moddepth;
128 output->ra = ra;
129 output->dec = dec;
130 output->stat = statval;
131 output->h0 = h0;
132 output->prob = prob;
133 output->proberrcode = proberrcode;
134 output->normalization = normalization;
135 output->templateVectorIndex = templateVectorIndex;
136 output->lineContamination = lineContamination;
137} // loadCandidateData()
138
139
140/**
141 * Analyze a single template
142 * \param [out] output Pointer to candidate structure
143 * \param [in] input Pointer to candidate structure
144 * \param [in] ffdata Pointer to ffdataStruct
145 * \param [in] aveNoise Pointer to REAL4VectorAligned of expected 2nd FFT background powers
146 * \param [in] aveTFnoisePerFbinRatio Pointer to REAL4VectorAligned of normalized power across the frequency band
147 * \param [in] params Pointer to UserInput_t
148 * \param [in] plan Pointer to REAL4FFTPlan
149 * \param [in] rng Pointer to gsl_rng
150 * \param [in] exactflag Boolean value to indicate using exact templates
151 * \return Status value
152 */
153INT4 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 )
154{
155
156 XLAL_CHECK( output != NULL && input != NULL && ffdata != NULL && aveNoise != NULL && aveTFnoisePerFbinRatio != NULL && params != NULL && plan != NULL && rng != NULL, XLAL_EINVAL );
157
158 INT4 proberrcode = 0;
159
160 //Allocate and make the template
161 TwoSpectTemplate *template = NULL;
162 XLAL_CHECK( ( template = createTwoSpectTemplate( params->maxTemplateLength ) ) != NULL, XLAL_EFUNC );
163 resetTwoSpectTemplate( template );
164 if ( exactflag ) {
165 XLAL_CHECK( makeTemplate( template, *input, params, plan ) == XLAL_SUCCESS, XLAL_EFUNC );
166 } else {
167 XLAL_CHECK( makeTemplateGaussians( template, *input, params ) == XLAL_SUCCESS, XLAL_EFUNC );
168 }
169
170 //Calculate R from the template and the data
171 REAL8 R = calculateR( ffdata->ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
172 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
173
174 //Calculate FAP
175 REAL8 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
176 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
177
178 //Estimate the h0 if R>0.0
179 REAL8 h0 = 0.0;
180 if ( R > 0.0 ) {
181 h0 = 2.7426 * pow( R / ( params->Tsft * params->Tobs ), 0.25 );
182 }
183
184 loadCandidateData( output, input->fsig, input->period, input->moddepth, input->ra, input->dec, R, h0, prob, proberrcode, 1.0, -1, 0 );
185
186 destroyTwoSpectTemplate( template );
187
188 return XLAL_SUCCESS;
189}
190
191INT4 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 )
192{
193
194 XLAL_CHECK( output != NULL && input != NULL && vector != NULL && ffdata != NULL && aveNoise != NULL && aveTFnoisePerFbinRatio != NULL && params != NULL && rng != NULL, XLAL_EINVAL );
195
196 TwoSpectTemplate *template = NULL;
197 XLAL_CHECK( ( template = createTwoSpectTemplate( templateLen ) ) != NULL, XLAL_EFUNC );
198
199 for ( UINT4 ii = 0; ii < input->length; ii++ ) {
200 INT4 proberrcode = 0;
201
202 //First convert the template to the right pixels
203 resetTwoSpectTemplate( template );
204 XLAL_CHECK( convertTemplateForSpecificFbin( template, vector->data[input->data[ii].templateVectorIndex], input->data[ii].fsig, params ) == XLAL_SUCCESS, XLAL_EFUNC );
205
206 //Calculate R from the template and the data
207 REAL8 R = calculateR( ffdata->ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
208 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
209
210 //Calculate FAP
211 REAL8 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
212 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
213
214 //Estimate the h0 if R>0.0
215 REAL8 h0 = 0.0;
216 if ( R > 0.0 ) {
217 h0 = 2.7426 * pow( R / ( params->Tsft * params->Tobs ), 0.25 );
218 }
219
220 if ( prob < output->data[output->length - 1].prob ) {
221 UINT4 insertionPoint = output->length - 1;
222 while ( insertionPoint > 0 && prob < output->data[insertionPoint - 1].prob ) {
223 insertionPoint--;
224 }
225 for ( INT4 kk = ( INT4 )output->length - 2; kk >= ( INT4 )insertionPoint; kk-- ) {
226 loadCandidateData( &( output->data[kk + 1] ), output->data[kk].fsig, output->data[kk].period, output->data[kk].moddepth, output->data[kk].ra, output->data[kk].dec, output->data[kk].stat, output->data[kk].h0, output->data[kk].prob, output->data[kk].proberrcode, output->data[kk].normalization, output->data[kk].templateVectorIndex, 0 );
227 }
228 loadCandidateData( &( output->data[insertionPoint] ), template->f0, template->period, template->moddepth, input->data[ii].ra, input->data[ii].dec, R, h0, prob, proberrcode, ffdata->tfnormalization, input->data[ii].templateVectorIndex, 0 );
229 if ( output->numofcandidates < output->length ) {
230 output->numofcandidates++;
231 }
232 }
233 }
234
235 destroyTwoSpectTemplate( template );
236
237 return XLAL_SUCCESS;
238}
239
240/**
241 * A brute force template search to find the most significant template around a candidate
242 * \param [out] output Pointer to candidate structure
243 * \param [in] input Input candidate structure
244 * \param [in] paramspace Pointer to TwoSpectParamSpaceSearchVals containing the parameter space to be searched
245 * \param [in] params Pointer to UserInput_t
246 * \param [in] ffdata Pointer to ffdataStruct
247 * \param [in] aveNoise Pointer to REAL4VectorAligned of 2nd FFT background powers
248 * \param [in] aveTFnoisePerFbinRatio Pointer to REAL4VectorAligned of normalized power across the frequency band
249 * \param [in] secondFFTplan Pointer to REAL4FFTPlan
250 * \param [in] rng Pointer to gsl_rng
251 * \param [in] useExactTemplates Boolean of 0 (use Gaussian templates) or 1 (use exact templates)
252 * \return Status value
253 */
254INT4 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 )
255{
256
257 XLAL_CHECK( output != NULL && params != NULL && aveNoise != NULL && aveTFnoisePerFbinRatio != NULL && secondFFTplan != NULL, XLAL_EINVAL );
258
259 fprintf( stderr, "Performing brute force template search... " );
260
261 REAL8Vector *trialf, *trialb, *trialp;
262 REAL8 fstepsize, dfstepsize;
263 //REAL4 tcohfactor = 1.49e-3*params->Tsft + 1.76; //From in-text equation after Eq. 23 of E.G. and K.R. 2011
264 REAL4 alpha0 = 45.0 * ( params->Tsft / 1800.0 ) + 30.0;
265 REAL8 log10templatefar = log10( params->tmplfar );
266
267 TwoSpectParamSpaceSearchVals search = {paramspace->fminimum, paramspace->fmaximum, paramspace->numfsteps, paramspace->numperiodslonger, paramspace->numperiodsshorter, paramspace->periodSpacingFactor, paramspace->dfmin, paramspace->dfmax, paramspace->numdfsteps};
268
269 //Set up parameters of modulation depth search
270 if ( search.dfmin < params->dfmin ) {
271 search.dfmin = params->dfmin;
272 }
273 if ( search.dfmax > params->dfmax ) {
274 search.dfmax = params->dfmax;
275 }
276 XLAL_CHECK( ( trialb = XLALCreateREAL8Vector( search.numdfsteps ) ) != NULL, XLAL_EFUNC );
277 if ( search.numdfsteps > 1 ) {
278 dfstepsize = ( search.dfmax - search.dfmin ) / ( REAL8 )( search.numdfsteps - 1 );
279 for ( UINT4 ii = 0; ii < search.numdfsteps; ii++ ) {
280 trialb->data[ii] = search.dfmin + dfstepsize * ii;
281 }
282 } else {
283 trialb->data[0] = 0.5 * ( search.dfmin + search.dfmax );
284 }
285
286 //Set up parameters of signal frequency search
287 if ( search.fminimum < params->fmin ) {
288 search.fminimum = params->fmin;
289 }
290 if ( search.fmaximum > params->fmin + params->fspan ) {
291 search.fmaximum = params->fmin + params->fspan;
292 }
293 XLAL_CHECK( ( trialf = XLALCreateREAL8Vector( search.numfsteps ) ) != NULL, XLAL_EFUNC );
294 if ( search.numfsteps > 1 ) {
295 fstepsize = ( search.fmaximum - search.fminimum ) / ( REAL8 )( search.numfsteps - 1 );
296 for ( UINT4 ii = 0; ii < search.numfsteps; ii++ ) {
297 trialf->data[ii] = search.fminimum + fstepsize * ii;
298 }
299 } else {
300 trialf->data[0] = 0.5 * ( search.fminimum + search.fmaximum );
301 }
302
303 //Search over numperiods different periods
304 XLAL_CHECK( ( trialp = XLALCreateREAL8Vector( search.numperiodslonger + search.numperiodsshorter + 1 ) ) != NULL, XLAL_EFUNC );
305
306 //Now search over the parameter space. Frequency, then modulation depth, then period
307 //Initialze best values as the initial point we are searching around
308 INT4 bestproberrcode = 0;
309 REAL8 bestf = 0.0, bestp = 0.0, bestdf = 0.0, bestR = 0.0, besth0 = 0.0, bestProb = 0.0;
310 candidate cand;
311 TwoSpectTemplate *template = NULL;
312 XLAL_CHECK( ( template = createTwoSpectTemplate( params->maxTemplateLength ) ) != NULL, XLAL_EFUNC );
313 farStruct *farval = NULL;
314 if ( params->calcRthreshold ) {
315 XLAL_CHECK( ( farval = createfarStruct() ) != NULL, XLAL_EFUNC );
316 }
317
318 INT4 startposition = search.numperiodsshorter, proberrcode = 0;
319 //Search over frequency
320 for ( UINT4 ii = 0; ii < trialf->length; ii++ ) {
321 //Search over modulation depth
322 for ( UINT4 jj = 0; jj < trialb->length; jj++ ) {
323 //Start with period of the first guess, then determine nearest neighbor from the
324 //modulation depth amplitude to find the other period guesses. These parameters
325 //are determined from simulation to scale the N.N. distance w.r.t. mod. depth with
326 //20% mismatch parameter
327 trialp->data[startposition] = input.period;
328 for ( UINT4 kk = 0; kk < search.numperiodsshorter; kk++ ) {
329 REAL8 nnp = search.periodSpacingFactor * trialp->data[startposition - kk] * trialp->data[startposition - kk] * ( 1 + trialp->data[startposition - kk] / ( alpha0 * sqrt( trialb->data[jj] ) * params->Tobs ) ) / ( alpha0 * params->Tobs * sqrt( trialb->data[jj] ) );
330 trialp->data[startposition - ( kk + 1 )] = trialp->data[startposition - kk] - nnp;
331 }
332 for ( UINT4 kk = 0; kk < search.numperiodslonger; kk++ ) {
333 REAL8 nnp = search.periodSpacingFactor * trialp->data[startposition + kk] * trialp->data[startposition + kk] * ( 1 + trialp->data[startposition + kk] / ( alpha0 * sqrt( trialb->data[jj] ) * params->Tobs ) ) / ( alpha0 * params->Tobs * sqrt( trialb->data[jj] ) );
334 trialp->data[startposition + ( kk + 1 )] = trialp->data[startposition + kk] + nnp;
335 }
336
337 //Search over period
338 for ( UINT4 kk = 0; kk < trialp->length; kk++ ) {
339 //Within boundaries?
340 if ( trialf->data[ii] >= params->fmin &&
341 trialf->data[ii] < ( params->fmin + params->fspan ) &&
342 trialb->data[jj] < maxModDepth( trialp->data[kk], params->Tsft ) &&
343 trialp->data[kk] > minPeriod( trialb->data[jj], params->Tsft ) &&
344 trialp->data[kk] <= ( 0.2 * params->Tobs ) &&
345 trialp->data[kk] >= ( 4.0 * params->Tsft ) &&
346 trialb->data[jj] >= params->dfmin &&
347 trialb->data[jj] <= params->dfmax &&
348 trialp->data[kk] <= params->Pmax &&
349 trialp->data[kk] >= params->Pmin ) {
350
351 loadCandidateData( &cand, trialf->data[ii], trialp->data[kk], trialb->data[jj], input.ra, input.dec, 0, 0, 0.0, 0, 0.0, -1, 0 );
352
353 resetTwoSpectTemplate( template );
354
355 if ( useExactTemplates ) {
356 XLAL_CHECK( makeTemplate( template, cand, params, secondFFTplan ) == XLAL_SUCCESS, XLAL_EFUNC );
357 } else {
358 XLAL_CHECK( makeTemplateGaussians( template, cand, params ) == XLAL_SUCCESS, XLAL_EFUNC );
359 }
360
361 if ( params->calcRthreshold && bestProb == 0.0 ) {
362 XLAL_CHECK( numericFAR( farval, template, params->tmplfar, aveNoise, aveTFnoisePerFbinRatio, params, rng, params->BrentsMethod ) == XLAL_SUCCESS, XLAL_EFUNC );
363 }
364
365 REAL8 R = calculateR( ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
366 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
367 REAL8 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
368 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
369 REAL8 h0 = 0.0;
370 if ( R > 0.0 ) {
371 h0 = 2.7426 * pow( R / ( params->Tsft * params->Tobs ), 0.25 );
372 }
373
374 //fprintf(stderr, "%.8g %.9g %g %.14g\n", trialf->data[ii], trialp->data[kk], trialb->data[jj], R);
375
376 if ( ( bestProb != 0.0 && prob < bestProb ) || ( bestProb == 0.0 && !params->calcRthreshold && prob < log10templatefar ) || ( bestProb == 0.0 && params->calcRthreshold && R > farval->far ) ) {
377 bestf = trialf->data[ii];
378 bestp = trialp->data[kk];
379 bestdf = trialb->data[jj];
380 bestR = R;
381 besth0 = h0;
382 bestProb = prob;
383 bestproberrcode = proberrcode;
384 }
385
386 } /* if within boundaries */
387 } /* for kk < trialp */
388 } /* for jj < trialb */
389 } /* for ii < trialf */
390 destroyTwoSpectTemplate( template );
391 template = NULL;
392 if ( params->calcRthreshold ) {
393 destroyfarStruct( farval );
394 farval = NULL;
395 }
396 XLALDestroyREAL8Vector( trialf );
397 XLALDestroyREAL8Vector( trialb );
398 XLALDestroyREAL8Vector( trialp );
399 trialf = NULL;
400 trialb = NULL;
401 trialp = NULL;
402
403 if ( bestProb == 0.0 ) {
404 loadCandidateData( output, input.fsig, input.period, input.moddepth, input.ra, input.dec, input.stat, input.h0, input.prob, input.proberrcode, input.normalization, input.templateVectorIndex, input.lineContamination );
405 } else {
406 loadCandidateData( output, bestf, bestp, bestdf, input.ra, input.dec, bestR, besth0, bestProb, bestproberrcode, input.normalization, input.templateVectorIndex, input.lineContamination );
407 }
408
409 fprintf( stderr, "done\n" );
410
411 return XLAL_SUCCESS;
412
413}
414
415/**
416 * A brute force template search to test templates around a candidate
417 * \param [out] output Pointer to a pointer of a candidateVector
418 * \param [in] input Input candidate structure
419 * \param [in] paramspace Pointer to TwoSpectParamSpaceSearchVals containing the parameter space to be searched
420 * \param [in] params Pointer to UserInput_t
421 * \param [in] ffdata Pointer to ffdataStruct
422 * \param [in] aveNoise Pointer to REAL4VectorAligned of 2nd FFT background powers
423 * \param [in] aveTFnoisePerFbinRatio Pointer to REAL4VectorAligned of normalized power across the frequency band
424 * \param [in] secondFFTplan Pointer to REAL4FFTPlan
425 * \param [in] rng Pointer to gsl_rng
426 * \param [in] useExactTemplates Boolean of 0 (use Gaussian templates) or 1 (use exact templates)
427 * \return Status value
428 */
429INT4 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 )
430{
431
432 XLAL_CHECK( *output != NULL && paramspace != NULL && params != NULL && ffdata != NULL && aveNoise != NULL && aveTFnoisePerFbinRatio != NULL && secondFFTplan != NULL && rng != NULL, XLAL_EINVAL );
433
434 REAL8Vector *trialf, *trialb, *trialp;
435 REAL8 fstepsize, dfstepsize;
436 REAL4 tcohfactor = 1.49e-3 * params->Tsft + 1.76; //From in-text equation after Eq. 23 of E.G. and K.R. 2011
437
438 TwoSpectParamSpaceSearchVals search = {paramspace->fminimum, paramspace->fmaximum, paramspace->numfsteps, paramspace->numperiodslonger, paramspace->numperiodsshorter, paramspace->periodSpacingFactor, paramspace->dfmin, paramspace->dfmax, paramspace->numdfsteps};
439
440 //Set up parameters of modulation depth search
441 if ( search.dfmin < params->dfmin ) {
442 search.dfmin = params->dfmin;
443 }
444 if ( search.dfmax > params->dfmax ) {
445 search.dfmax = params->dfmax;
446 }
447 XLAL_CHECK( ( trialb = XLALCreateREAL8Vector( search.numdfsteps ) ) != NULL, XLAL_EFUNC );
448 if ( search.numdfsteps > 1 ) {
449 dfstepsize = ( search.dfmax - search.dfmin ) / ( REAL8 )( search.numdfsteps - 1 );
450 for ( UINT4 ii = 0; ii < search.numdfsteps; ii++ ) {
451 trialb->data[ii] = search.dfmin + dfstepsize * ii;
452 }
453 } else {
454 trialb->data[0] = 0.5 * ( search.dfmin + search.dfmax );
455 }
456
457 //Set up parameters of signal frequency search
458 if ( search.fminimum < params->fmin ) {
459 search.fminimum = params->fmin;
460 }
461 if ( search.fmaximum > params->fmin + params->fspan ) {
462 search.fmaximum = params->fmin + params->fspan;
463 }
464 XLAL_CHECK( ( trialf = XLALCreateREAL8Vector( search.numfsteps ) ) != NULL, XLAL_EFUNC );
465 if ( search.numfsteps > 1 ) {
466 fstepsize = ( search.fmaximum - search.fminimum ) / ( REAL8 )( search.numfsteps - 1 );
467 for ( UINT4 ii = 0; ii < search.numfsteps; ii++ ) {
468 trialf->data[ii] = search.fminimum + fstepsize * ii;
469 }
470 } else {
471 trialf->data[0] = 0.5 * ( search.fminimum + search.fmaximum );
472 }
473
474 //Search over numperiods different periods
475 XLAL_CHECK( ( trialp = XLALCreateREAL8Vector( search.numperiodslonger + search.numperiodsshorter + 1 ) ) != NULL, XLAL_EFUNC );
476
477 //Now search over the parameter space. Frequency, then modulation depth, then period
478 candidate cand;
479 TwoSpectTemplate *template = NULL;
480 XLAL_CHECK( ( template = createTwoSpectTemplate( params->maxTemplateLength ) ) != NULL, XLAL_EFUNC );
481
482 INT4 startposition = search.numperiodsshorter, proberrcode = 0;
483 //Search over frequency
484 for ( UINT4 ii = 0; ii < trialf->length; ii++ ) {
485 //Search over modulation depth
486 for ( UINT4 jj = 0; jj < trialb->length; jj++ ) {
487 //Start with period of the first guess, then determine nearest neighbor from the
488 //modulation depth amplitude to find the other period guesses. These parameters
489 //are determined from simulation to scale the N.N. distance w.r.t. mod. depth with
490 //20% mismatch parameter
491 trialp->data[startposition] = input.period;
492 for ( UINT4 kk = 0; kk < search.numperiodsshorter; kk++ ) {
493 REAL8 nnp = search.periodSpacingFactor * trialp->data[startposition - kk] * trialp->data[startposition - kk] * ( 1 + trialp->data[startposition - kk] / tcohfactor / params->Tobs ) / tcohfactor / params->Tobs * sqrt( 3.6e-3 / trialb->data[jj] );
494 trialp->data[startposition - ( kk + 1 )] = trialp->data[startposition - kk] - nnp;
495 }
496 for ( UINT4 kk = 0; kk < search.numperiodslonger; kk++ ) {
497 REAL8 nnp = search.periodSpacingFactor * trialp->data[startposition + kk] * trialp->data[startposition + kk] * ( 1 + trialp->data[startposition + kk] / tcohfactor / params->Tobs ) / tcohfactor / params->Tobs * sqrt( 3.6e-3 / trialb->data[jj] );
498 trialp->data[startposition + ( kk + 1 )] = trialp->data[startposition + kk] + nnp;
499 }
500
501 //Search over period
502 for ( UINT4 kk = 0; kk < trialp->length; kk++ ) {
503 //Within boundaries?
504 if ( trialf->data[ii] >= params->fmin &&
505 trialf->data[ii] < ( params->fmin + params->fspan ) &&
506 trialb->data[jj] < maxModDepth( trialp->data[kk], params->Tsft ) &&
507 trialp->data[kk] > minPeriod( trialb->data[jj], params->Tsft ) &&
508 trialp->data[kk] <= ( 0.2 * params->Tobs ) &&
509 trialp->data[kk] >= ( 4.0 * params->Tsft ) &&
510 trialb->data[jj] >= params->dfmin &&
511 trialb->data[jj] <= params->dfmax &&
512 trialp->data[kk] <= params->Pmax &&
513 trialp->data[kk] >= params->Pmin ) {
514
515 loadCandidateData( &cand, trialf->data[ii], trialp->data[kk], trialb->data[jj], input.ra, input.dec, 0, 0, 0.0, 0, 0.0, -1, 0 );
516
517 resetTwoSpectTemplate( template );
518
519 if ( useExactTemplates ) {
520 XLAL_CHECK( makeTemplate( template, cand, params, secondFFTplan ) == XLAL_SUCCESS, XLAL_EFUNC );
521 } else {
522 XLAL_CHECK( makeTemplateGaussians( template, cand, params ) == XLAL_SUCCESS, XLAL_EFUNC );
523 }
524
525 REAL8 R = calculateR( ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
526 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
527 REAL8 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
528 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
529 REAL8 h0 = 0.0;
530 if ( R > 0.0 ) {
531 h0 = 2.7426 * pow( R / ( params->Tsft * params->Tobs ), 0.25 );
532 }
533
534 //Resize the output candidate vector if necessary
535 if ( ( *output )->numofcandidates == ( *output )->length - 1 ) {
536 XLAL_CHECK( ( *output = resizecandidateVector( *output, 2 * ( *output )->length ) ) != NULL, XLAL_EFUNC );
537 }
538
539 loadCandidateData( &( ( *output )->data[( *output )->numofcandidates] ), trialf->data[ii], trialp->data[kk], trialb->data[jj], input.ra, input.dec, R, h0, prob, proberrcode, input.normalization, input.templateVectorIndex, input.lineContamination );
540 ( *output )->numofcandidates++;
541
542 } /* if within boundaries */
543 } /* for kk < trialp */
544 } /* for jj < trialb */
545 } /* for ii < trialf */
546 destroyTwoSpectTemplate( template );
547 XLALDestroyREAL8Vector( trialf );
548 XLALDestroyREAL8Vector( trialb );
549 XLALDestroyREAL8Vector( trialp );
550
551 return XLAL_SUCCESS;
552
553}
554
555
556/**
557 * A brute force template search to find the most significant template around a putative source whose parameters are somewhat constrained
558 * \param [out] output Pointer to a pointer of a candidateVector
559 * \param [in] fminimum Lower frequency bound to search (inclusive)
560 * \param [in] fspan Span of the frequency band (inclusive of endpoint)
561 * \param [in] period Specific orbital period (measured in seconds)
562 * \param [in] asini Specific projected semi-major axis (measured in light seconds)
563 * \param [in] asinisigma Uncertainty on the specific asini value (measured in light seconds)
564 * \param [in] skypos SkyPosition struct of the sky position (in RA and DEC) being searched
565 * \param [in] params Pointer to UserInput_t
566 * \param [in] ffdata Pointer to ffdataStruct
567 * \param [in] aveNoise Pointer to REAL4VectorAligned of 2nd FFT background powers
568 * \param [in] aveTFnoisePerFbinRatio Pointer to REAL4VectorAligned of normalized power across the frequency band
569 * \param [in] trackedlines Pointer to REAL4VectorSequence of lines (allowed to be NULL if no lines)
570 * \param [in] secondFFTplan Pointer to REAL4FFTPlan
571 * \param [in] rng Pointer to gsl_rng
572 * \param [in] useExactTemplates Boolean of 0 (use Gaussian templates) or 1 (use exact templates)
573 * \return Status value
574 */
575INT4 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 )
576{
577
578 XLAL_CHECK( *output != NULL && params != NULL && ffdata != NULL && aveNoise != NULL && aveTFnoisePerFbinRatio != NULL && secondFFTplan != NULL && rng != NULL, XLAL_EINVAL );
579
580 REAL8Vector *trialf;
581 REAL8Vector *trialdf;
582 REAL8 fstepsize;
583 REAL8 dfstepsize;
584
585 //Set up parameters of signal frequency search
586 INT4 numfsteps = ( INT4 )round( 2.0 * fspan * params->Tsft ) + 1;
587 XLAL_CHECK( ( trialf = XLALCreateREAL8Vector( numfsteps ) ) != NULL, XLAL_EFUNC );
588 fstepsize = fspan / ( REAL8 )( numfsteps - 1 );
589 for ( INT4 ii = 0; ii < numfsteps; ii++ ) {
590 trialf->data[ii] = fminimum + fstepsize * ii;
591 }
592
593 //Now search over the frequencies
594 INT4 proberrcode = 0;
595 candidate cand;
596 TwoSpectTemplate *template = NULL;
597 XLAL_CHECK( ( template = createTwoSpectTemplate( params->maxTemplateLength ) ) != NULL, XLAL_EFUNC );
598
599 //Search over frequency
600 for ( UINT4 ii = 0; ii < trialf->length; ii++ ) {
601
602 //Set up parameters of signal modulation depth search
603 /* Modulation depth is 2*pi*f*asini*period, or rearranged
604 0.8727*(f/1000.0)*(7200.0/period)*asini
605 Assuming sigma = 0.18 uncertainty in an asini of 1.44 for
606 Scorpius X-1 and giving +/- 3*sigma leeway, the conservative
607 number of df steps should cover
608 0.8727*(fmax/1000.0)*(7200.0/period)*6*0.18
609 with, as empirical testing has found necessary, a spacing of
610 4*Tsft */
611 /* While this initialization could be moved inside the loop
612 that searches over frequency, it is slightly faster not to have to
613 recalculate these variables every time,
614 and it gives us a bit of extra data*/
615 //REAL8 asinisigma = 0.18; typical for Scorpius X-1 with 2014 data
616 REAL8 moddepth = 0.8727 * ( trialf->data[ii] / 1000.0 ) * ( 7200.0 / period ) * asini;
617 //fprintf(stderr,"Making the first computation involving asinisigma, for moddepthmin\n");
618 //Note, 6*asinsigma for a span of plus/minus 3*asinisigma
619 REAL8 moddepthspan = 0.8727 * ( trialf->data[numfsteps - 1] / 1000.0 ) * ( 7200.0 / period ) * 6 * asinisigma;
620 //fprintf(stderr,"intended moddepthspan: %f \n", moddepthspan);
621 //fprintf(stderr,"Done with moddepthspan, making moddepthmin\n");
622 REAL8 moddepthmin = moddepth - 0.5 * moddepthspan;
623 //fprintf(stderr,"intended moddepthmin: %f \n", moddepthmin);
624 INT4 numdfsteps = ( INT4 )round( 4.0 * moddepthspan * params->Tsft ) + 1;
625 //fprintf(stderr,"intended numdfsteps: %d \n", numdfsteps);
626 trialdf = XLALCreateREAL8Vector( numdfsteps );
627 XLAL_CHECK( trialdf != NULL, XLAL_EFUNC );
628 if ( numdfsteps > 1 ) {
629 dfstepsize = moddepthspan / ( REAL8 )( numdfsteps - 1 );
630 } else {
631 dfstepsize = 0;
632 }
633 for ( INT4 jj = 0; jj < numdfsteps; jj++ ) {
634 trialdf->data[jj] = moddepthmin + dfstepsize * jj;
635 }
636
637 for ( UINT4 jj = 0; jj < trialdf->length; jj++ ) {
638 //Determine modulation depth
639 //REAL8 moddepth = 0.8727*(trialf->data[ii]/1000.0)*(7200.0/period)*asini;
640
641 //load candidate
642 //printf(stderr,"Loading candidate. Remember to get the RA and dec from outside in production run\n");
643 loadCandidateData( &cand, trialf->data[ii], period, trialdf->data[jj], skypos.longitude, skypos.latitude, 0, 0, 0.0, 0, 0.0, -1, 0 );
644
645 //Make the template
646 resetTwoSpectTemplate( template );
647 if ( useExactTemplates != 0 ) {
648 XLAL_CHECK( makeTemplate( template, cand, params, secondFFTplan ) == XLAL_SUCCESS, XLAL_EFUNC );
649 } else {
650 XLAL_CHECK( makeTemplateGaussians( template, cand, params ) == XLAL_SUCCESS, XLAL_EFUNC );
651 }
652
653 REAL8 R = calculateR( ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
654 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
655 REAL8 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
656 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
657 REAL8 h0 = 0.0;
658 if ( R > 0.0 ) {
659 h0 = 2.7426 * pow( R / ( params->Tsft * params->Tobs ), 0.25 );
660 }
661
662 //Line contamination?
663 BOOLEAN lineContamination = 0;
664 if ( trackedlines != NULL ) {
665 UINT4 kk = 0;
666 while ( kk < trackedlines->length && lineContamination == 0 ) {
667 if ( 2.0 * trialdf->data[jj] >= ( trackedlines->data[kk * 3 + 2] - trackedlines->data[kk * 3 + 1] ) ) {
668 if ( ( trackedlines->data[kk * 3 + 2] >= ( REAL4 )( trialf->data[ii] - trialdf->data[jj] ) && trackedlines->data[kk * 3 + 2] <= ( REAL4 )( trialf->data[ii] + trialdf->data[jj] ) ) ||
669 ( trackedlines->data[kk * 3 + 1] >= ( REAL4 )( trialf->data[ii] - trialdf->data[jj] ) && trackedlines->data[kk * 3 + 1] <= ( REAL4 )( trialf->data[ii] + trialdf->data[jj] ) ) ) {
670 lineContamination = 1;
671 }
672 } // if the band spanned by the line is smaller than the band spanned by the signal
673 else {
674 if ( ( ( REAL4 )( trialf->data[ii] + trialdf->data[jj] ) >= trackedlines->data[kk * 3 + 1] && ( REAL4 )( trialf->data[ii] + trialdf->data[jj] ) <= trackedlines->data[kk * 3 + 2] ) ||
675 ( ( REAL4 )( trialf->data[ii] - trialdf->data[jj] ) >= trackedlines->data[kk * 3 + 1] && ( REAL4 )( trialf->data[ii] - trialdf->data[jj] ) <= trackedlines->data[kk * 3 + 2] ) ) {
676 lineContamination = 1;
677 }
678 } // instead if the band spanned by the line is larger than the band spanned by the signal
679 kk++;
680 } // while kk < trackedlines->length && lineContamination==0
681 } // if trackedlines != NULL
682
683 //Resize the output candidate vector if necessary
684 if ( ( *output )->numofcandidates == ( *output )->length - 1 ) {
685 *output = resizecandidateVector( *output, 2 * ( ( *output )->length ) );
686 XLAL_CHECK( *output != NULL, XLAL_EFUNC );
687 }
688
689 loadCandidateData( &( ( *output )->data[( *output )->numofcandidates] ), trialf->data[ii], period, trialdf->data[jj], skypos.longitude, skypos.latitude, R, h0, prob, proberrcode, 0.0, -1, lineContamination );
690 ( *output )->numofcandidates++;
691 } /* for jj < trialdf */
692 XLALDestroyREAL8Vector( trialdf );
693 trialdf = NULL;
694 } /* for ii < trialf */
695 destroyTwoSpectTemplate( template );
696 template = NULL;
697 XLALDestroyREAL8Vector( trialf );
698 trialf = NULL;
699
700 return XLAL_SUCCESS;
701
702}
703
704/**
705 * A brute force template search to find the most significant template at a fixed modulation depth around a putative source whose parameters are somewhat constrained
706 * \param [out] output Pointer to a pointer of a candidateVector
707 * \param [in] dffixed Modulation depth (fixed by user)
708 * \param [in] fminimum Lower frequency bound to search (inclusive)
709 * \param [in] fspan Span of the frequency band (inclusive of endpoint)
710 * \param [in] period Specific orbital period (measured in seconds)
711 * \param [in] skypos SkyPosition struct of the sky position (in RA and DEC) being searched
712 * \param [in] params Pointer to UserInput_t
713 * \param [in] ffdata Pointer to ffdataStruct
714 * \param [in] aveNoise Pointer to REAL4VectorAligned of 2nd FFT background powers
715 * \param [in] aveTFnoisePerFbinRatio Pointer to REAL4VectorAligned of normalized power across the frequency band
716* \param [in] trackedlines Pointer to REAL4VectorSequence of lines (allowed to be NULL if no lines)
717 * \param [in] secondFFTplan Pointer to REAL4FFTPlan
718 * \param [in] rng Pointer to gsl_rng
719 * \param [in] useExactTemplates Boolean of 0 (use Gaussian templates) or 1 (use exact templates)
720 * \return Status value
721 */
722
723
724INT4 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 )
725{
726
727 XLAL_CHECK( *output != NULL && dffixed != NULL && params != NULL && ffdata != NULL && aveNoise != NULL && aveTFnoisePerFbinRatio != NULL && secondFFTplan != NULL && rng != NULL, XLAL_EINVAL );
728
729 REAL8Vector *trialf;
730 REAL8Vector *trialdf;
731 REAL8 fstepsize;
732
733 // Create the vector trialdf
734 XLAL_CHECK( ( trialdf = XLALCreateREAL8Vector( ( dffixed->length ) ) ) != NULL, XLAL_EFUNC );
735 for ( UINT4 ii = 0; ii < dffixed->length; ii++ ) {
736 XLAL_CHECK( XLALParseStringValueAsREAL8( &( trialdf->data[ii] ), dffixed->data[ii] ) == XLAL_SUCCESS, XLAL_EFUNC );
737
738 // Check that the specified df is ok; if not, return an error and end the program.
739 XLAL_CHECK( ( trialdf->data[ii] < maxModDepth( period, params->Tsft ) ) && ( trialdf->data[ii] > 0.5 / params->Tsft ), XLAL_EFAILED, "ERROR: Modulation depth must be between %.5f and %.5f.\n", 0.5 / params->Tsft, maxModDepth( period, params->Tsft ) );
740 }
741
742 //Set up parameters of signal frequency search
743 INT4 numfsteps = ( INT4 )round( 2.0 * fspan * params->Tsft ) + 1;
744 XLAL_CHECK( ( trialf = XLALCreateREAL8Vector( numfsteps ) ) != NULL, XLAL_EFUNC );
745 fstepsize = fspan / ( REAL8 )( numfsteps - 1 );
746 for ( INT4 ii = 0; ii < numfsteps; ii++ ) {
747 trialf->data[ii] = fminimum + fstepsize * ii;
748 }
749
750 //Now search over the frequencies
751 INT4 proberrcode = 0;
752 candidate cand;
753 TwoSpectTemplate *template = NULL;
754 XLAL_CHECK( ( template = createTwoSpectTemplate( params->maxTemplateLength ) ) != NULL, XLAL_EFUNC );
755
756 // loop over dfs
757 for ( UINT4 jj = 0; jj < trialdf->length; jj++ ) {
758
759 //Search over frequency
760 for ( UINT4 ii = 0; ii < trialf->length; ii++ ) {
761
762 loadCandidateData( &cand, trialf->data[ii], period, trialdf->data[jj], skypos.longitude, skypos.latitude, 0, 0, 0.0, 0, 0.0, -1, 0 );
763
764 //Make the template
765 resetTwoSpectTemplate( template );
766 if ( useExactTemplates != 0 ) {
767 XLAL_CHECK( makeTemplate( template, cand, params, secondFFTplan ) == XLAL_SUCCESS, XLAL_EFUNC );
768 } else {
769 XLAL_CHECK( makeTemplateGaussians( template, cand, params ) == XLAL_SUCCESS, XLAL_EFUNC );
770 }
771
772 REAL8 R = calculateR( ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
773 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
774 REAL8 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
775 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
776 REAL8 h0 = 0.0;
777 if ( R > 0.0 ) {
778 h0 = 2.7426 * pow( R / ( params->Tsft * params->Tobs ), 0.25 );
779 }
780
781 //Line contamination?
782 BOOLEAN lineContamination = 0;
783 if ( trackedlines != NULL ) {
784 UINT4 kk = 0;
785 while ( kk < trackedlines->length && lineContamination == 0 ) {
786 if ( 2.0 * trialdf->data[jj] >= ( trackedlines->data[kk * 3 + 2] - trackedlines->data[kk * 3 + 1] ) ) {
787 if ( ( trackedlines->data[kk * 3 + 2] >= ( REAL4 )( trialf->data[ii] - trialdf->data[jj] ) && trackedlines->data[kk * 3 + 2] <= ( REAL4 )( trialf->data[ii] + trialdf->data[jj] ) ) ||
788 ( trackedlines->data[kk * 3 + 1] >= ( REAL4 )( trialf->data[ii] - trialdf->data[jj] ) && trackedlines->data[kk * 3 + 1] <= ( REAL4 )( trialf->data[ii] + trialdf->data[jj] ) ) ) {
789 lineContamination = 1;
790 }
791 } // if the band spanned by the line is smaller than the band spanned by the signal
792 else {
793 if ( ( ( REAL4 )( trialf->data[ii] + trialdf->data[jj] ) >= trackedlines->data[kk * 3 + 1] && ( REAL4 )( trialf->data[ii] + trialdf->data[jj] ) <= trackedlines->data[kk * 3 + 2] ) ||
794 ( ( REAL4 )( trialf->data[ii] - trialdf->data[jj] ) >= trackedlines->data[kk * 3 + 1] && ( REAL4 )( trialf->data[ii] - trialdf->data[jj] ) <= trackedlines->data[kk * 3 + 2] ) ) {
795 lineContamination = 1;
796 }
797 } // instead if the band spanned by the line is larger than the band spanned by the signal
798 kk++;
799 } // while kk < trackedlines->length && lineContamination==0
800 } // if trackedlines != NULL
801
802 //Resize the output candidate vector if necessary
803 if ( ( *output )->numofcandidates == ( *output )->length - 1 ) {
804 *output = resizecandidateVector( *output, 2 * ( ( *output )->length ) );
805 XLAL_CHECK( *output != NULL, XLAL_EFUNC );
806 }
807
808 loadCandidateData( &( ( *output )->data[( *output )->numofcandidates] ), trialf->data[ii], period, trialdf->data[jj], skypos.longitude, skypos.latitude, R, h0, prob, proberrcode, 0.0, -1, lineContamination );
809 ( *output )->numofcandidates++;
810
811 } /* for ii < trialf */
812 } /* for jj < trialdf */
813
814 XLALDestroyREAL8Vector( trialdf );
815 trialdf = NULL;
816 destroyTwoSpectTemplate( template );
817 template = NULL;
818 XLALDestroyREAL8Vector( trialf );
819 trialf = NULL;
820
821 return XLAL_SUCCESS;
822
823}
824
825/**
826 * Cluster candidates by frequency, period, and modulation depth using templates
827 * \param [out] output Pointer to pointer of a candidateVector
828 * \param [in] input Pointer to a candidateVector
829 * \param [in] ffdata Pointer to ffdataStruct
830 * \param [in] params Pointer to UserInput_t
831 * \param [in] ffplanenoise Pointer to REAL4VectorAligned of 2nd FFT background powers
832 * \param [in] fbinaveratios Pointer to REAL4VectorAligned of normalized SFT background
833 * \param [in] rng Pointer to gsl_rng
834 * \param [in] exactflag Flag to use Gaussian templates (0) or exact templates (1)
835 * \return Status value
836 */
837INT4 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 )
838{
839
840 XLAL_CHECK( *output != NULL && input != NULL && ffdata != NULL && params != NULL && ffplanenoise != NULL && fbinaveratios != NULL && rng != NULL, XLAL_EINVAL );
841
842 UINT4 loc, loc2;
843 REAL8 avefsig, aveperiod, mindf, maxdf;
844
845 //Allocate int vectors for storage
846 INT4Vector *locs = NULL, *locs2 = NULL, *usedcandidate = NULL;
847 XLAL_CHECK( ( locs = XLALCreateINT4Vector( input->numofcandidates ) ) != NULL, XLAL_EFUNC );
848 XLAL_CHECK( ( locs2 = XLALCreateINT4Vector( input->numofcandidates ) ) != NULL, XLAL_EFUNC );
849 XLAL_CHECK( ( usedcandidate = XLALCreateINT4Vector( input->numofcandidates ) ) != NULL, XLAL_EFUNC );
850
851 //Initialize arrays
852 for ( UINT4 ii = 0; ii < input->numofcandidates; ii++ ) {
853 locs->data[ii] = -1;
854 locs2->data[ii] = -1;
855 usedcandidate->data[ii] = 0;
856 }
857
858 //Make FFT plan if exactflag is given
859 REAL4FFTPlan *plan = NULL;
860 if ( exactflag == 1 ) {
861 XLAL_CHECK( ( plan = XLALCreateForwardREAL4FFTPlan( ffdata->numffts, 1 ) ) != NULL, XLAL_EFUNC );
862 }
863
864 for ( UINT4 ii = 0; ii < input->numofcandidates; ii++ ) {
865
866 //Make note of first candidate available
867 locs->data[0] = ii;
868 loc = 1;
869
870 INT4 foundany = 0; //Switch to determine if any other candidates in the group. 1 if true
871 INT4 iter = 1;
872 //Find any in the list that are within +1/2 bin in first FFT frequency
873 for ( UINT4 jj = ii + 1; jj < input->numofcandidates; jj++ ) {
874 if ( usedcandidate->data[jj] == 0 && ( input->data[jj].fsig - input->data[locs->data[0]].fsig <= 0.5 * iter / params->Tsft + 1.0e-6 && input->data[jj].fsig - input->data[locs->data[0]].fsig >= -0.25 * iter / params->Tsft ) ) {
875 locs->data[loc] = jj;
876 loc++;
877 if ( foundany == 0 ) {
878 foundany = 1;
879 }
880 }
881 } /* for jj < input->numofcandidates */
882 //Keep checking as long as there are more connected frequencies going higher in frequency
883 while ( foundany == 1 ) {
884 foundany = 0;
885 iter++;
886 for ( UINT4 jj = ii + 1; jj < input->numofcandidates; jj++ ) {
887 if ( usedcandidate->data[jj] == 0 && ( input->data[jj].fsig - input->data[locs->data[0]].fsig - 0.25 / params->Tsft <= 0.5 * iter / params->Tsft && input->data[jj].fsig - input->data[locs->data[0]].fsig + 0.25 / params->Tsft >= 0.5 * iter / params->Tsft ) ) {
888 locs->data[loc] = jj;
889 loc++;
890 if ( foundany == 0 ) {
891 foundany = 1;
892 }
893 }
894 } /* for jj < input->numofcandidates */
895 } /* while foundany==1 */
896 //Now check frequencies 1/2 bin below and keep going as long as there are more connected frequencies
897 foundany = 1;
898 iter = 0;
899 while ( foundany == 1 ) {
900 foundany = 0;
901 iter++;
902 for ( UINT4 jj = ii + 1; jj < input->numofcandidates; jj++ ) {
903 if ( usedcandidate->data[jj] == 0 && ( input->data[locs->data[0]].fsig - input->data[jj].fsig - 0.25 / params->Tsft <= 0.5 * iter / params->Tsft && input->data[locs->data[0]].fsig - input->data[jj].fsig + 0.25 / params->Tsft >= 0.5 * iter / params->Tsft ) ) {
904 locs->data[loc] = jj;
905 loc++;
906 if ( foundany == 0 ) {
907 foundany = 1;
908 }
909 }
910 } /* for jj < input->numofcandidates */
911 } /* while foundany==1 */
912
913 //Using the list of locations, find the subset that have periods within 1 bin
914 //of the second FFT frequencies
915 UINT4 subsetloc = 0, nextsubsetloc = 0, subsetlocset = 0;
916 loc2 = 0;
917 for ( UINT4 jj = subsetloc; jj < loc; jj++ ) {
918 if ( usedcandidate->data[locs->data[jj]] == 0 && fabs( params->Tobs / input->data[locs->data[jj]].period - params->Tobs / input->data[locs->data[subsetloc]].period ) <= 1.0 ) {
919 locs2->data[loc2] = locs->data[jj];
920 loc2++;
921 } else if ( usedcandidate->data[locs->data[jj]] == 0 && subsetlocset == 0 ) {
922 subsetlocset = 1;
923 nextsubsetloc = jj;
924 }
925
926 if ( jj + 1 == loc ) {
927 if ( subsetlocset == 1 ) {
928 subsetloc = nextsubsetloc; //Reset subsetloc and jj to the next candidate period
929 jj = subsetloc - 1;
930 subsetlocset = 0; //Reset the logic of whether there are any more periods to go
931 }
932
933 //find best candidate moddepth
934 fprintf( stderr, "Finding best modulation depth with number to try %d\n", loc2 );
935 avefsig = 0.0, aveperiod = 0.0, mindf = 0.0, maxdf = 0.0;
936 REAL8 weight = 0.0, bestmoddepth = 0.0, bestR = 0.0, besth0 = 0.0, bestProb = 0.0;
937 INT4 bestproberrcode = 0;
938 for ( UINT4 kk = 0; kk < loc2; kk++ ) {
939 avefsig += input->data[locs2->data[kk]].fsig * ( input->data[locs2->data[kk]].prob * input->data[locs2->data[kk]].prob );
940 aveperiod += input->data[locs2->data[kk]].period * ( input->data[locs2->data[kk]].prob * input->data[locs2->data[kk]].prob );
941 weight += input->data[locs2->data[kk]].prob * input->data[locs2->data[kk]].prob;
942 if ( mindf > input->data[locs2->data[kk]].moddepth || mindf == 0.0 ) {
943 mindf = input->data[locs2->data[kk]].moddepth;
944 }
945 if ( maxdf < input->data[locs2->data[kk]].moddepth ) {
946 maxdf = input->data[locs2->data[kk]].moddepth;
947 }
948
949 if ( loc2 == 1 && input->data[locs2->data[kk]].fsig >= params->fmin && input->data[locs2->data[kk]].fsig < ( params->fmin + params->fspan ) && input->data[locs2->data[kk]].period >= params->Pmin && input->data[locs2->data[kk]].period <= params->Pmax ) {
950 besth0 = input->data[locs2->data[kk]].h0;
951 bestmoddepth = input->data[locs2->data[kk]].moddepth;
952 bestR = input->data[locs2->data[kk]].stat;
953 bestProb = input->data[locs2->data[kk]].prob;
954 bestproberrcode = input->data[locs2->data[kk]].proberrcode;
955 }
956
957 usedcandidate->data[locs2->data[kk]] = 1;
958 } /* for kk < loc2 */
959 avefsig = avefsig / weight;
960 aveperiod = aveperiod / weight;
961
962 INT4 proberrcode = 0;
963
964 if ( loc2 > 1 && aveperiod >= params->Pmin - 1.0 && aveperiod <= params->Pmax + 1.0 ) {
965 UINT4 numofmoddepths = ( UINT4 )floorf( 2 * ( maxdf - mindf ) * params->Tsft ) + 1;
966 candidate cand;
967 TwoSpectTemplate *template = NULL;
968 XLAL_CHECK( ( template = createTwoSpectTemplate( params->maxTemplateLength ) ) != NULL, XLAL_EFUNC );
969
970 for ( UINT4 kk = 0; kk < numofmoddepths; kk++ ) {
971 if ( ( mindf + kk * 0.5 / params->Tsft ) >= params->dfmin && ( mindf + kk * 0.5 / params->Tsft ) <= params->dfmax ) {
972
973 loadCandidateData( &cand, avefsig, aveperiod, mindf + kk * 0.5 / params->Tsft, input->data[0].ra, input->data[0].dec, 0, 0, 0.0, 0, 0.0, -1, 0 );
974
975 if ( exactflag == 1 ) {
976 XLAL_CHECK( makeTemplate( template, cand, params, plan ) == XLAL_SUCCESS, XLAL_EFUNC );
977 } else {
978 XLAL_CHECK( makeTemplateGaussians( template, cand, params ) == XLAL_SUCCESS, XLAL_EFUNC );
979 }
980
981 REAL8 R = calculateR( ffdata->ffdata, template, ffplanenoise, fbinaveratios );
982 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
983 REAL8 prob = probR( template, ffplanenoise, fbinaveratios, R, params, rng, &proberrcode );
984 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
985
986 if ( prob < bestProb ) {
987 bestmoddepth = mindf + kk * 0.5 / params->Tsft;
988 bestR = R;
989 bestProb = prob;
990 bestproberrcode = proberrcode;
991 }
992 } /* if test moddepth is within user specified range */
993 } /* for kk < numofmoddepths */
994
995 destroyTwoSpectTemplate( template );
996 template = NULL;
997 } /* if loc2 > 1 ... */
998
999 if ( bestProb != 0.0 ) {
1000 if ( bestR > 0.0 ) {
1001 besth0 = 2.7426 * pow( bestR / ( params->Tsft * params->Tobs ), 0.25 );
1002 } else {
1003 besth0 = 0.0;
1004 }
1005
1006 if ( ( *output )->numofcandidates == ( *output )->length - 1 ) {
1007 XLAL_CHECK( ( *output = resizecandidateVector( *output, 2 * ( *output )->length ) ) != NULL, XLAL_EFUNC );
1008 }
1009 loadCandidateData( &( ( *output )->data[( *output )->numofcandidates] ), avefsig, aveperiod, bestmoddepth, input->data[0].ra, input->data[0].dec, bestR, besth0, bestProb, bestproberrcode, input->data[0].normalization, -1, 0 );
1010 ( *output )->numofcandidates++;
1011 }
1012
1013 loc2 = 0;
1014 } /* if jj+1 == loc */
1015 } /* for jj < loc */
1016
1017 //Find location of first entry to be searched next time or finish the cluster search
1018 for ( UINT4 jj = ii; jj < input->numofcandidates; jj++ ) {
1019 if ( usedcandidate->data[jj] == 0 ) {
1020 ii = jj - 1;
1021 jj = input->numofcandidates - 1;
1022 } else if ( jj == input->numofcandidates - 1 ) {
1023 ii = input->numofcandidates - 1;
1024 }
1025 }
1026
1027 //Reinitialize values, just in case
1028 for ( UINT4 jj = 0; jj < locs->length; jj++ ) {
1029 locs->data[jj] = -1;
1030 locs2->data[jj] = -1;
1031 }
1032 } /* for ii < numofcandidates */
1033
1034 //Destroy stuff
1035 XLALDestroyINT4Vector( locs );
1036 XLALDestroyINT4Vector( locs2 );
1037 XLALDestroyINT4Vector( usedcandidate );
1038 if ( exactflag == 1 ) {
1039 XLALDestroyREAL4FFTPlan( plan );
1040 }
1041
1042 fprintf( stderr, "Clustering done with candidates = %d\n", ( *output )->numofcandidates );
1043 fprintf( LOG, "Clustering done with candidates = %d\n", ( *output )->numofcandidates );
1044
1045 return XLAL_SUCCESS;
1046
1047} /* clusterCandidates() */
1048
1049
1050/**
1051 * Function to test the IHS candidates against Gaussian templates
1052 * \param [out] output Pointer to pointer of a candidateVector
1053 * \param [in] ihsCandidates Pointer to candidateVector of IHS candidates
1054 * \param [in] ffdata Pointer to ffdataStruct
1055 * \param [in] aveNoise Pointer to REAL4VectorAligned of 2nd FFT background powers
1056 * \param [in] aveTFnoisePerFbinRatio Pointer to REAL4VectorAligned of normalized SFT background spectra
1057 * \param [in] pos The current sky position
1058 * \param [in] params Pointer to UserInput_t
1059 * \param [in] rng Pointer to gsl_rng
1060 * \return Status value
1061 */
1062INT4 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 )
1063{
1064
1065 XLAL_CHECK( *output != NULL && ihsCandidates != NULL && ffdata != NULL && aveNoise != NULL && aveTFnoisePerFbinRatio != NULL && params != NULL && rng != NULL, XLAL_EINVAL );
1066
1067 //R probability calculator errorcode
1068 INT4 proberrcode = 0;
1069
1070 //Allocate memory for FAR struct
1071 farStruct *farval = NULL;
1072 XLAL_CHECK( ( farval = createfarStruct() ) != NULL, XLAL_EFUNC );
1073
1074 //Allocate memory for template
1075 TwoSpectTemplate *template = NULL;
1076 XLAL_CHECK( ( template = createTwoSpectTemplate( params->maxTemplateLength ) ) != NULL, XLAL_EFUNC );
1077
1078 INT4 candidatesoutsideofmainULrange = 0;
1079 REAL8 log10templatefar = log10( params->tmplfar );
1080
1081 for ( UINT4 ii = 0; ii < ihsCandidates->numofcandidates; ii++ ) {
1082 //Assess the IHS candidate if the signal is away from the band edges, the modulation depth is greater or equal to minimum allowed and less than or equal to the maximum allowed, and if the period/modulation depth combo is within allowable limits for a template to be made. We will cut the period space in the next step.
1083 if ( ihsCandidates->data[ii].fsig >= params->fmin && ihsCandidates->data[ii].fsig < ( params->fmin + params->fspan ) ) {
1084 if ( params->followUpOutsideULrange || ( ihsCandidates->data[ii].fsig >= params->ULfmin && ihsCandidates->data[ii].fsig <= ( params->ULfmin + params->ULfspan ) && ihsCandidates->data[ii].moddepth >= params->ULminimumDeltaf && ihsCandidates->data[ii].moddepth <= params->ULmaximumDeltaf ) ) {
1085
1086 resetTwoSpectTemplate( template );
1087
1088 REAL8 R, prob, bestPeriod = 0.0, bestR = 0.0, bestProb = 0.0;
1089 INT4 bestproberrcode = 0;
1090
1091 if ( ihsCandidates->data[ii].period >= fmax( 4.0 * params->Tsft, minPeriod( ihsCandidates->data[ii].moddepth, params->Tsft ) ) && ihsCandidates->data[ii].period <= ( 0.2 * params->Tobs ) ) {
1092 //Make a Gaussian train template
1093 XLAL_CHECK( makeTemplateGaussians( template, ihsCandidates->data[ii], params ) == XLAL_SUCCESS, XLAL_EFUNC );
1094
1095 //Estimate the FAR for these bin weights if the option was given
1096 if ( params->calcRthreshold ) {
1097 XLAL_CHECK( numericFAR( farval, template, params->tmplfar, aveNoise, aveTFnoisePerFbinRatio, params, rng, params->BrentsMethod ) == XLAL_SUCCESS, XLAL_EFUNC );
1098 }
1099
1100 //Caclulate R and probability noise caused the candidate
1101 R = calculateR( ffdata->ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
1102 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
1103 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
1104 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
1105
1106 /* Note the candidate if R exceeds the FAR or check other possibilities of different
1107 periods */
1108 if ( ( !params->calcRthreshold && prob < log10templatefar ) || ( params->calcRthreshold && R > farval->far ) ) {
1109 bestR = R;
1110 bestProb = prob;
1111 bestPeriod = ihsCandidates->data[ii].period;
1112 } /* if prob<log10templatefar || R > farval->far */
1113 } // if within moddepth/period range
1114
1115 // longer or shorter
1116 REAL8 periodfact = 0.0;
1117 for ( UINT4 jj = 0; jj <= 1; jj++ ) {
1118 //Shift by harmonics
1119 for ( INT4 kk = 2; kk <= params->periodHarmToCheck; kk++ ) {
1120 if ( jj == 0 ) {
1121 periodfact = 1.0 / ( REAL8 )kk;
1122 } else {
1123 periodfact = ( REAL8 )kk;
1124 }
1125 if ( ihsCandidates->data[ii].period * periodfact >= fmax( params->Pmin, minPeriod( ihsCandidates->data[ii].moddepth, params->Tsft ) ) && ihsCandidates->data[ii].period * periodfact <= fmin( params->Pmax, params->Tobs * 0.2 ) ) {
1126 ihsCandidates->data[ii].period *= periodfact;
1127 XLAL_CHECK( makeTemplateGaussians( template, ihsCandidates->data[ii], params ) == XLAL_SUCCESS, XLAL_EFUNC );
1128 R = calculateR( ffdata->ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
1129 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
1130 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
1131 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
1132 if ( params->calcRthreshold && bestProb == 0.0 ) {
1133 XLAL_CHECK( numericFAR( farval, template, params->tmplfar, aveNoise, aveTFnoisePerFbinRatio, params, rng, params->BrentsMethod ) == XLAL_SUCCESS, XLAL_EFUNC );
1134 }
1135 if ( ( bestProb != 0.0 && prob < bestProb ) || ( bestProb == 0.0 && !params->calcRthreshold && prob < log10templatefar ) || ( bestProb == 0.0 && params->calcRthreshold && R > farval->far ) ) {
1136 bestPeriod = ihsCandidates->data[ii].period;
1137 bestR = R;
1138 bestProb = prob;
1139 bestproberrcode = proberrcode;
1140 }
1141 ihsCandidates->data[ii].period /= periodfact; //reset the period back to the original value
1142 } // in range?
1143 } // shift by harmonics for kk <= inputParams->periodHarmToCheck (harmonics)
1144
1145 //shift by fractions
1146 for ( INT4 kk = 1; kk <= params->periodFracToCheck; kk++ ) {
1147 if ( jj == 0 ) {
1148 periodfact = ( kk + 1.0 ) / ( kk + 2.0 );
1149 } else {
1150 periodfact = ( kk + 2.0 ) / ( kk + 1.0 );
1151 }
1152 if ( ihsCandidates->data[ii].period * periodfact >= fmax( params->Pmin, minPeriod( ihsCandidates->data[ii].moddepth, params->Tsft ) ) && ihsCandidates->data[ii].period * periodfact <= fmin( params->Pmax, params->Tobs * 0.2 ) ) {
1153 ihsCandidates->data[ii].period *= periodfact;
1154 XLAL_CHECK( makeTemplateGaussians( template, ihsCandidates->data[ii], params ) == XLAL_SUCCESS, XLAL_EFUNC );
1155 R = calculateR( ffdata->ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
1156 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
1157 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
1158 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
1159 if ( params->calcRthreshold && bestProb == 0.0 ) {
1160 XLAL_CHECK( numericFAR( farval, template, params->tmplfar, aveNoise, aveTFnoisePerFbinRatio, params, rng, params->BrentsMethod ) == XLAL_SUCCESS, XLAL_EFUNC );
1161 }
1162 if ( ( bestProb != 0.0 && prob < bestProb ) || ( bestProb == 0.0 && !params->calcRthreshold && prob < log10templatefar ) || ( bestProb == 0.0 && params->calcRthreshold && R > farval->far ) ) {
1163 bestPeriod = ihsCandidates->data[ii].period;
1164 bestR = R;
1165 bestProb = prob;
1166 bestproberrcode = proberrcode;
1167 }
1168 ihsCandidates->data[ii].period /= periodfact; //reset the period back to the original value
1169 } // in range?
1170 } // shift by fractions kk <= inputParams->periodFracToCheck
1171 } // longer or shorter
1172
1173 if ( bestProb != 0.0 ) {
1174 REAL8 h0 = 0.0;
1175 if ( bestR > 0.0 ) {
1176 h0 = 2.7426 * sqrt( sqrt( bestR / ( params->Tsft * params->Tobs ) ) ); //Now compute the h0 value
1177 }
1178
1179 if ( ( *output )->numofcandidates == ( *output )->length - 1 ) {
1180 XLAL_CHECK( ( *output = resizecandidateVector( *output, 2 * ( *output )->length ) ) != NULL, XLAL_EFUNC );
1181 }
1182 loadCandidateData( &( ( *output )->data[( *output )->numofcandidates] ), ihsCandidates->data[ii].fsig, bestPeriod, ihsCandidates->data[ii].moddepth, pos.longitude, pos.latitude, bestR, h0, bestProb, bestproberrcode, ihsCandidates->data[ii].normalization, -1, ihsCandidates->data[ii].lineContamination );
1183 ( *output )->numofcandidates++;
1184
1185 } /* if bestR != 0.0, add candidate or replace if something better is found */
1186 } /* if within UL boundaries */
1187 else {
1188 candidatesoutsideofmainULrange++;
1189 }
1190 } /* if within outer boundaries */
1191 } /* for ii < numofcandidates */
1192
1193 //Destroy allocated memory
1194 destroyTwoSpectTemplate( template );
1195 template = NULL;
1196 destroyfarStruct( farval );
1197 farval = NULL;
1198
1199 fprintf( stderr, "%d remaining candidate(s) inside UL range.\n", ihsCandidates->numofcandidates - candidatesoutsideofmainULrange );
1200 fprintf( stderr, "Initial stage done with candidates = %d\n", ( *output )->numofcandidates );
1201 fprintf( LOG, "Initial stage done with candidates = %d\n", ( *output )->numofcandidates );
1202
1203 return XLAL_SUCCESS;
1204
1205} /* testIHScandidates() */
1206
1207
1208/**
1209 * Test each of the templates in a TwoSpectTemplateVector and keep the top 10
1210 * This will not check the false alarm probability of any R value less than 0.
1211 * \param [out] output Pointer to pointer of a candidateVector storing a list of all candidates
1212 * \param [in] templateVec Pointer to a TwoSpectTemplateVector containing all the templates to be searched
1213 * \param [in] ffdata Pointer to ffdataStruct
1214 * \param [in] aveNoise Pointer to REAL4VectorAligned of 2nd FFT background powers
1215 * \param [in] aveTFnoisePerFbinRatio Pointer to REAL4VectorAligned of normalized SFT background spectra
1216 * \param [in] skypos The current sky position
1217 * \param [in] params Pointer to UserInput_t
1218 * \param [in] rng Pointer to gsl_rng
1219 * \param [in] templateLen Maximum length of a template
1220 * \return Status value
1221 */
1222INT4 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 )
1223{
1224
1225 XLAL_CHECK( output != NULL && templateVec != NULL && ffdata != NULL && aveNoise != NULL && aveTFnoisePerFbinRatio != NULL && params != NULL && rng != NULL, XLAL_EINVAL );
1226
1227 fprintf( stderr, "Testing TwoSpectTemplateVector... " );
1228
1229 TwoSpectTemplate *template = NULL;
1230 XLAL_CHECK( ( template = createTwoSpectTemplate( templateLen ) ) != NULL, XLAL_EFUNC );
1231
1232 INT4 proberrcode = 0;
1233
1234 FILE *RVALS = NULL;
1235 if ( XLALUserVarWasSet( &params->saveRvalues ) ) {
1236 XLAL_CHECK( ( RVALS = fopen( params->saveRvalues, "w" ) ) != NULL, XLAL_EIO, "Couldn't open %s for writing", params->saveRvalues );
1237 }
1238
1239 UINT4 numfbins = ( UINT4 )round( params->fspan * params->Tsft );
1240 for ( UINT4 ii = 0; ii < numfbins; ii++ ) {
1241 REAL8 freq = params->fmin + ii / params->Tsft;
1242 for ( UINT4 jj = 0; jj < templateVec->length; jj++ ) {
1243 if ( templateVec->data[jj]->templatedata->data[0] == 0.0 ) {
1244 break;
1245 }
1246
1247 XLAL_CHECK( convertTemplateForSpecificFbin( template, templateVec->data[jj], freq, params ) == XLAL_SUCCESS, XLAL_EFUNC );
1248
1249 REAL8 R = calculateR( ffdata->ffdata, template, aveNoise, aveTFnoisePerFbinRatio );
1250 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
1251 REAL8 prob = 0.0, h0 = 0.0;
1252 if ( R > 0.0 ) {
1253 prob = probR( template, aveNoise, aveTFnoisePerFbinRatio, R, params, rng, &proberrcode );
1254 XLAL_CHECK( xlalErrno == 0, XLAL_EFUNC );
1255 h0 = 2.7426 * pow( R / ( params->Tsft * params->Tobs ), 0.25 );
1256 }
1257
1258 if ( XLALUserVarWasSet( &params->saveRvalues ) ) {
1259 fprintf( RVALS, "%g\n", R );
1260 }
1261
1262 if ( prob < output->data[output->length - 1].prob ) {
1263 UINT4 insertionPoint = output->length - 1;
1264 while ( insertionPoint > 0 && prob < output->data[insertionPoint - 1].prob ) {
1265 insertionPoint--;
1266 }
1267 for ( INT4 kk = ( INT4 )output->length - 2; kk >= ( INT4 )insertionPoint; kk-- ) {
1268 loadCandidateData( &( output->data[kk + 1] ), output->data[kk].fsig, output->data[kk].period, output->data[kk].moddepth, output->data[kk].ra, output->data[kk].dec, output->data[kk].stat, output->data[kk].h0, output->data[kk].prob, output->data[kk].proberrcode, output->data[kk].normalization, output->data[kk].templateVectorIndex, output->data[kk].lineContamination );
1269 }
1270 loadCandidateData( &( output->data[insertionPoint] ), template->f0, template->period, template->moddepth, skypos.longitude, skypos.latitude, R, h0, prob, proberrcode, ffdata->tfnormalization, jj, 0 );
1271 if ( output->numofcandidates < output->length ) {
1272 output->numofcandidates++;
1273 }
1274 }
1275 }
1276 }
1277
1278 destroyTwoSpectTemplate( template );
1279
1280 if ( XLALUserVarWasSet( &params->saveRvalues ) ) {
1281 fclose( RVALS );
1282 }
1283
1284 fprintf( stderr, "done\n" );
1285
1286 return XLAL_SUCCESS;
1287
1288}
1289
1290
1291/**
1292 * Keep the most significant candidates, potentially reducing the number of candidates if there are more than allowed
1293 * \param [in] input Pointer to input candidateVector
1294 * \param [in] params Pointer to UserInput_t
1295 * \return Pointer to newly allocated candidateVector containing reduced number of candidates
1296 */
1297candidateVector *keepMostSignificantCandidates( const candidateVector *input, const UserInput_t *params )
1298{
1299
1300 XLAL_CHECK_NULL( input != NULL && params != NULL, XLAL_EINVAL );
1301
1302 fprintf( stderr, "Reducing total number of IHS candidates %d to user input %d\n", input->numofcandidates, params->keepOnlyTopNumIHS );
1303 fprintf( LOG, "Reducing total number of IHS candidates %d to user input %d\n", input->numofcandidates, params->keepOnlyTopNumIHS );
1304
1305 candidateVector *output = NULL;
1306
1307 //If the number to keep is > 0 and the number of candidates is less than the number to keep,
1308 //just move the input vector to the output vector
1309 if ( params->keepOnlyTopNumIHS > 0 && ( INT4 )input->numofcandidates <= params->keepOnlyTopNumIHS ) {
1310 XLAL_CHECK_NULL( ( output = createcandidateVector( input->numofcandidates ) ) != NULL, XLAL_EFUNC );
1311
1312 for ( UINT4 ii = 0; ii < input->numofcandidates; ii++ ) {
1313 loadCandidateData( &( output->data[ii] ), input->data[ii].fsig, input->data[ii].period, input->data[ii].moddepth, input->data[ii].ra, input->data[ii].dec, input->data[ii].stat, input->data[ii].h0, input->data[ii].prob, input->data[ii].proberrcode, input->data[ii].normalization, input->data[ii].templateVectorIndex, input->data[ii].lineContamination );
1314 }
1315 output->numofcandidates = input->numofcandidates;
1316
1317 } else if ( params->keepOnlyTopNumIHS > 0 && ( INT4 )input->numofcandidates > params->keepOnlyTopNumIHS ) {
1318 //If keep is > 0 and the number of candidates is > the number to keep,
1319 //we sort through the list and find the most significant candidates to keep
1320 XLAL_CHECK_NULL( ( output = createcandidateVector( params->keepOnlyTopNumIHS ) ) != NULL, XLAL_EFUNC );
1321
1322 for ( UINT4 ii = 0; ii < output->length; ii++ ) {
1323 REAL8 highestsignificance = 0.0;
1324 INT4 candidateWithHighestSignificance = 0;
1325 for ( UINT4 jj = 0; jj < input->numofcandidates; jj++ ) {
1326 if ( input->data[jj].prob > highestsignificance ) {
1327 highestsignificance = input->data[jj].prob;
1328 candidateWithHighestSignificance = jj;
1329 }
1330 }
1331
1332 loadCandidateData( &( output->data[ii] ), input->data[candidateWithHighestSignificance].fsig, input->data[candidateWithHighestSignificance].period, input->data[candidateWithHighestSignificance].moddepth, input->data[candidateWithHighestSignificance].ra, input->data[candidateWithHighestSignificance].dec, input->data[candidateWithHighestSignificance].stat, input->data[candidateWithHighestSignificance].h0, input->data[candidateWithHighestSignificance].prob, input->data[candidateWithHighestSignificance].proberrcode, input->data[candidateWithHighestSignificance].normalization, input->data[candidateWithHighestSignificance].templateVectorIndex, input->data[candidateWithHighestSignificance].lineContamination );
1333
1334 input->data[candidateWithHighestSignificance].prob = 0.0;
1335
1336 }
1337 output->numofcandidates = params->keepOnlyTopNumIHS;
1338
1339 } else {
1340 //Otherwise, we need to fail
1341 fprintf( stderr, "%s: keepOnlyTopNumIHS given (%d) is not greater than 0, but it should be to use this function.\n", __func__, params->keepOnlyTopNumIHS );
1342 XLAL_ERROR_NULL( XLAL_EINVAL );
1343 }
1344
1345 return output;
1346
1347} /* keepMostSignificantCandidates() */
1348
1349
1350/**
1351 * Calculate the R statistic from equation 13 of E. Goetz and K. Riles (2011)
1352 * \param [in] ffdata Pointer to REAL4VectorAligned of the 2nd FFT data
1353 * \param [in] template Pointer to the template
1354 * \param [in] noise Pointer to the REAL4VectorAligned containing the background 2nd FFT powers
1355 * \param [in] fbinaveratios Pointer to the REAL4VectorAligned of normalized SFT background powers
1356 * \return Value of the R statistic
1357 */
1358REAL8 calculateR( const REAL4VectorAligned *ffdata, const TwoSpectTemplate *template, const REAL4VectorAligned *noise, const REAL4VectorAligned *fbinaveratios )
1359{
1360
1361 XLAL_CHECK_REAL8( ffdata != NULL && template != NULL && noise != NULL && fbinaveratios != NULL, XLAL_EINVAL );
1362
1363 REAL8 sumofsqweights = 0.0;
1364 for ( UINT4 ii = 0; ii < template->templatedata->length; ii++ ) if ( template->templatedata->data[ii] != 0.0 ) {
1365 sumofsqweights += ( template->templatedata->data[ii] * template->templatedata->data[ii] );
1366 }
1367 XLAL_CHECK_REAL8( sumofsqweights != 0.0, XLAL_EFPDIV0 );
1368 REAL8 sumofsqweightsinv = 1.0 / sumofsqweights;
1369
1370 UINT4 numfprbins = noise->length;
1371
1372 REAL8 R = 0.0;
1373 for ( UINT4 ii = 0; ii < template->templatedata->length; ii++ ) {
1374 if ( template->templatedata->data[ii] != 0.0 ) {
1375 UINT4 firstfreqbin = template->pixellocations->data[ii] / numfprbins;
1376 UINT4 secfreqbin = template->pixellocations->data[ii] - firstfreqbin * numfprbins;
1377 R += ( ffdata->data[ template->pixellocations->data[ii] ] - noise->data[secfreqbin] * fbinaveratios->data[firstfreqbin] ) * template->templatedata->data[ii] * sumofsqweightsinv;
1378 }
1379 }
1380
1381 return R;
1382
1383} /* calculateR() */
1384
1385INT4 writeCandidateVector2File( const CHAR *outputfile, const candidateVector *input )
1386{
1387 XLAL_CHECK( outputfile != NULL && input != NULL, XLAL_EINVAL );
1388
1389 FILE *CANDFILE = NULL;
1390 struct stat XLAL_INIT_DECL( buf );
1391 if ( stat( outputfile, &buf ) == -1 && errno == ENOENT ) {
1392 XLAL_CHECK( ( CANDFILE = fopen( outputfile, "w" ) ) != NULL, XLAL_EIO, "Couldn't fopen file %s to output candidates\n", outputfile );
1393 fprintf( CANDFILE, "# TwoSpect candidates output file\n" );
1394 fprintf( CANDFILE, "# Freq Period Mod. depth RA DEC Statistic val. Est. h0 False alarm prob. TF normalization Template vec. num. Line contam.\n" );
1395 } else {
1396 XLAL_CHECK( ( CANDFILE = fopen( outputfile, "a" ) ) != NULL, XLAL_EIO, "Couldn't fopen file %s to output candidates\n", outputfile );
1397 }
1398
1399 for ( UINT4 ii = 0; ii < input->numofcandidates; ii++ ) {
1400 fprintf( CANDFILE, "%.6f %.6f %.7f %.4f %.4f %.4f %g %.4f %g %d %d\n", input->data[ii].fsig, input->data[ii].period, input->data[ii].moddepth, input->data[ii].ra, input->data[ii].dec, input->data[ii].stat, input->data[ii].h0, input->data[ii].prob, input->data[ii].normalization, input->data[ii].templateVectorIndex, input->data[ii].lineContamination );
1401 }
1402
1403 fclose( CANDFILE );
1404
1405 return XLAL_SUCCESS;
1406}
1407
1408/**
1409 * Calculates maximum modulation depth allowed, equation 6 of E. Goetz and K. Riles (2011)
1410 * \param [in] period Orbital period value
1411 * \param [in] cohtime SFT coherence length
1412 * \return Maximum modulation depth allowed
1413 */
1414REAL8 maxModDepth( const REAL8 period, const REAL8 cohtime )
1415{
1416 REAL8 maxB = 0.5 * period / ( cohtime * cohtime );
1417 return maxB;
1418} /* maxModDepth() */
1419
1420
1421/**
1422 * Calculates minimum period allowed, equation 6 of E. Goetz and K. Riles (2011)
1423 * \param [in] moddepth Modulation depth value
1424 * \param [in] cohtime SFT coherence length
1425 * \return Maximum modulation depth allowed
1426 */
1427REAL8 minPeriod( const REAL8 moddepth, const REAL8 cohtime )
1428{
1429 REAL8 minP = 2.0 * moddepth * cohtime * cohtime;
1430 return minP;
1431} /* minPeriod() */
__func__
#define __func__
log an I/O error, i.e.
Definition: CrossCorrToplist.c:865
fprintf
#define fprintf
LOG
FILE * LOG
Definition: TwoSpect.c:55
templateSearch_scox1Style
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
minPeriod
REAL8 minPeriod(const REAL8 moddepth, const REAL8 cohtime)
Calculates minimum period allowed, equation 6 of E.
Definition: candidates.c:1427
keepMostSignificantCandidates
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
createcandidateVector
candidateVector * createcandidateVector(const UINT4 length)
Allocate a candidateVector.
Definition: candidates.c:32
testIHScandidates
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
analyzeCandidatesTemplateFromVector
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
calculateR
REAL8 calculateR(const REAL4VectorAligned *ffdata, const TwoSpectTemplate *template, const REAL4VectorAligned *noise, const REAL4VectorAligned *fbinaveratios)
Calculate the R statistic from equation 13 of E.
Definition: candidates.c:1358
resizecandidateVector
candidateVector * resizecandidateVector(candidateVector *vector, const UINT4 length)
Resize a candidateVector.
Definition: candidates.c:61
writeCandidateVector2File
INT4 writeCandidateVector2File(const CHAR *outputfile, const candidateVector *input)
Definition: candidates.c:1385
bruteForceTemplateSearch
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
maxModDepth
REAL8 maxModDepth(const REAL8 period, const REAL8 cohtime)
Calculates maximum modulation depth allowed, equation 6 of E.
Definition: candidates.c:1414
bruteForceTemplateTest
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
testTwoSpectTemplateVector
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
clusterCandidates
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
analyzeOneTemplate
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
destroycandidateVector
void destroycandidateVector(candidateVector *vector)
Free a candidateVector.
Definition: candidates.c:89
loadCandidateData
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
templateSearch_fixedDf
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
createfarStruct
farStruct * createfarStruct(void)
Allocate memory for farStruct.
Definition: falsealarm.c:32
probR
REAL8 probR(const TwoSpectTemplate *template, const REAL4VectorAligned *ffplanenoise, const REAL4VectorAligned *fbinaveratios, const REAL8 R, const UserInput_t *params, const gsl_rng *rng, INT4 *errcode)
Analytically calculate the probability of a true signal using the Davies' method.
Definition: falsealarm.c:350
numericFAR
INT4 numericFAR(farStruct *output, const TwoSpectTemplate *template, const REAL8 thresh, const REAL4VectorAligned *ffplanenoise, const REAL4VectorAligned *fbinaveratios, const UserInput_t *inputParams, const gsl_rng *rng, const INT4 method)
Numerically solve for the FAR of the R statistic from the weights using the Davies algorithm and a ro...
Definition: falsealarm.c:140
destroyfarStruct
void destroyfarStruct(farStruct *farstruct)
Destroy an farStruct.
Definition: falsealarm.c:46
BOOLEAN
unsigned char BOOLEAN
REAL8
double REAL8
XLAL_INIT_DECL
#define XLAL_INIT_DECL(var,...)
CHAR
char CHAR
UINT4
uint32_t UINT4
INT4
int32_t INT4
REAL4
float REAL4
XLALMalloc
void * XLALMalloc(size_t n)
XLALRealloc
void * XLALRealloc(void *p, size_t n)
XLALFree
void XLALFree(void *p)
XLALCreateForwardREAL4FFTPlan
REAL4FFTPlan * XLALCreateForwardREAL4FFTPlan(UINT4 size, int measurelvl)
XLALDestroyREAL4FFTPlan
void XLALDestroyREAL4FFTPlan(REAL4FFTPlan *plan)
XLALUserVarWasSet
int XLALUserVarWasSet(const void *cvar)
XLALParseStringValueAsREAL8
int XLALParseStringValueAsREAL8(REAL8 *valREAL8, const char *valString)
XLALCreateINT4Vector
INT4Vector * XLALCreateINT4Vector(UINT4 length)
XLALDestroyINT4Vector
void XLALDestroyINT4Vector(INT4Vector *vector)
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLAL_ERROR_VOID
#define XLAL_ERROR_VOID(...)
XLAL_ERROR_NULL
#define XLAL_ERROR_NULL(...)
xlalErrno
#define xlalErrno
XLAL_CHECK
#define XLAL_CHECK(assertion,...)
XLAL_CHECK_VOID
#define XLAL_CHECK_VOID(assertion,...)
XLAL_CHECK_REAL8
#define XLAL_CHECK_REAL8(assertion,...)
XLAL_CHECK_NULL
#define XLAL_CHECK_NULL(assertion,...)
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLAL_EIO
XLAL_EIO
XLAL_EINVAL
XLAL_EINVAL
XLAL_EFAILED
XLAL_EFAILED
XLAL_EFPDIV0
XLAL_EFPDIV0
data
float data[BLOCKSIZE]
Definition: hwinject.c:360
pos
pos
lalpulsar_PiecewiseSearch.fmax
fmax
Definition: lalpulsar_PiecewiseSearch.py:496
lalpulsar_PiecewiseSearch.h0
h0
Definition: lalpulsar_PiecewiseSearch.py:601
lalpulsar_PiecewiseSearch.fmin
fmin
Definition: lalpulsar_PiecewiseSearch.py:495
make_frame_cache.output
output
Definition: make_frame_cache.py:53
test_flat_prior.stat
stat
Definition: test_flat_prior.py:108
test_simulateCW.freq
float freq
Definition: test_simulateCW.py:54
INT4Vector::data
INT4 * data
INT4Vector::length
UINT4 length
REAL8Vector::data
REAL8 * data
SkyPosition::longitude
REAL8 longitude
SkyPosition::latitude
REAL8 latitude
TwoSpectTemplate::templatedata
REAL4VectorAligned * templatedata
Definition: TwoSpectTypes.h:240
TwoSpectTemplateVector::data
TwoSpectTemplate ** data
Definition: TwoSpectTypes.h:248
UserInput_t
candidate::templateVectorIndex
INT4 templateVectorIndex
Definition: TwoSpectTypes.h:169
candidate::moddepth
REAL8 moddepth
Definition: TwoSpectTypes.h:161
candidate::lineContamination
BOOLEAN lineContamination
Definition: TwoSpectTypes.h:170
candidate::normalization
REAL8 normalization
Definition: TwoSpectTypes.h:168
candidate::proberrcode
INT4 proberrcode
Definition: TwoSpectTypes.h:167
candidate::period
REAL8 period
Definition: TwoSpectTypes.h:160
candidateVector::data
candidate * data
Definition: TwoSpectTypes.h:174
ffdataStruct::ffdata
REAL4VectorAligned * ffdata
Definition: TwoSpectTypes.h:150
ffdataStruct::tfnormalization
REAL8 tfnormalization
Definition: TwoSpectTypes.h:151
makeTemplate
INT4 makeTemplate(TwoSpectTemplate *output, const candidate input, const UserInput_t *params, const REAL4FFTPlan *plan)
Make an template based on FFT of sinc squared functions.
Definition: templates.c:727
resetTwoSpectTemplate
void resetTwoSpectTemplate(TwoSpectTemplate *template)
Reset the values in a TwoSpectTemplate.
Definition: templates.c:60
createTwoSpectTemplate
TwoSpectTemplate * createTwoSpectTemplate(const UINT4 length)
Allocate a new TwoSpectTemplate.
Definition: templates.c:35
destroyTwoSpectTemplate
void destroyTwoSpectTemplate(TwoSpectTemplate *template)
Free a TwoSpectTemplate.
Definition: templates.c:78
convertTemplateForSpecificFbin
INT4 convertTemplateForSpecificFbin(TwoSpectTemplate *output, const TwoSpectTemplate *input, const REAL8 freq, const UserInput_t *params)
Convert an arbitrary frequency bin template into a template for a specific frequency bin.
Definition: templates.c:330
makeTemplateGaussians
INT4 makeTemplateGaussians(TwoSpectTemplate *output, const candidate input, const UserInput_t *params)
Make an estimated template based on FFT of train of Gaussians.
Definition: templates.c:289