Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-5e288d3
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
compareCandidates.c
Go to the documentation of this file.
1/*
2* Copyright (C) 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 <stdio.h>
21#include <stdlib.h>
22#include <math.h>
23
24#include <gsl/gsl_sort.h>
25#include <gsl/gsl_roots.h>
26
27#include <lal/UserInput.h>
28#include <lal/LALStdlib.h>
29#include <lal/SeqFactories.h>
30
31typedef struct {
32 CHAR *infile1;
33 CHAR *infile2;
34 CHAR *outfile1;
35 CHAR *outfile2;
36 CHAR *outfile3;
37 CHAR *finalOutfile;
38 REAL8 Tobs;
39 REAL8 Tcoh;
40 REAL8 fdiff_allowed;
41 REAL8 dfdiff_allowed;
42 REAL8 skydiff_allowed;
43} UserVariables_t;
44
45struct solver_params {
46 double fvalue;
47 double *data_array;
48};
49
50INT4 InitUserVars( UserVariables_t *uvar, int argc, char *argv[] );
51REAL8 f_diff( double indexval, void *params );
52INT4 readInCoincidentOutliers( double **output, int **output_job, int *output_numOutliers, const char *infilename );
53INT4 readInIFOoutliersAndSort( double **output, int **output_job, int *output_numCoincident, const char *infilename );
54
55INT4 InitUserVars( UserVariables_t *uvar, int argc, char *argv[] )
56{
57 XLAL_CHECK( uvar != NULL, XLAL_EINVAL, "Invalid NULL input 'uvar'\n" );
58 XLAL_CHECK( argv != NULL, XLAL_EINVAL, "Invalid NULL input 'argv'\n" );
59
60 uvar->Tcoh = 1800;
61
62 XLALRegisterUvarMember( infile1, STRING, 0, REQUIRED, "Input file 1" );
63 XLALRegisterUvarMember( infile2, STRING, 0, REQUIRED, "Input file 2" );
64 XLALRegisterUvarMember( outfile1, STRING, 0, REQUIRED, "Temporary output file with all coincident candidates" );
65 XLALRegisterUvarMember( outfile2, STRING, 0, REQUIRED, "Temporary output file with subset of coincient outliers" );
66 XLALRegisterUvarMember( outfile3, STRING, 0, REQUIRED, "Temporary output file with alternate subset of coincident outliers" );
67 XLALRegisterUvarMember( finalOutfile, STRING, 0, REQUIRED, "Final output file of coincident outliers" );
68 XLALRegisterUvarMember( Tobs, REAL8, 0, REQUIRED, "Total observation time (in seconds)" );
69 XLALRegisterUvarMember( Tcoh, REAL8, 0, OPTIONAL, "SFT coherence time (in seconds)" );
70 XLALRegisterUvarMember( fdiff_allowed, REAL8, 0, REQUIRED, "Difference in frequencies allowed (in Hz)" );
71 XLALRegisterUvarMember( dfdiff_allowed, REAL8, 0, REQUIRED, "Difference in modulation depth allowed (in Hz)" );
72 XLALRegisterUvarMember( skydiff_allowed, REAL8, 0, REQUIRED, "Difference in sky location allowed (in radians) at fiducial frequency 200 Hz" );
73
74 BOOLEAN should_exit = 0;
75 XLAL_CHECK( XLALUserVarReadAllInput( &should_exit, argc, argv, lalPulsarVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
76 if ( should_exit ) {
77 exit( 1 );
78 }
79
80 return XLAL_SUCCESS;
81}
82double f_diff( double indexval, void *params )
83{
84 struct solver_params *p = ( struct solver_params * )params;
85 return p->fvalue - p->data_array[( int )round( indexval ) * 9];
86}
87INT4 readInCoincidentOutliers( double **output, int **output_job, int *output_numCoincident, const char *infilename )
88{
89 FILE *CANDS = NULL;
90 XLAL_CHECK( ( CANDS = fopen( infilename, "r" ) ) != NULL, XLAL_EIO, "Can't fopen %s", infilename );
91 //Determines number of candidates in the file
92 INT4 ch, count = 0;
93 do {
94 ch = fgetc( CANDS );
95 if ( ch == '\n' ) {
96 count++;
97 }
98 } while ( ch != EOF );
99 double *allcands = NULL;
100 XLAL_CHECK( ( allcands = ( double * )XLALMalloc( sizeof( double ) * count * 18 ) ) != NULL, XLAL_ENOMEM );
101 int *allcands_job = NULL;
102 XLAL_CHECK( ( allcands_job = ( int * )XLALMalloc( sizeof( int ) * count * 2 ) ) != NULL, XLAL_ENOMEM );
103 rewind( CANDS );
104 //Put the data into the array
105 for ( INT4 ii = 0; ii < count; ii++ ) {
106 fscanf( CANDS, "%la %la %la %la %la %la %la %la %la %d %la %la %la %la %la %la %la %la %la %d", &( allcands[ii * 18] ), &( allcands[ii * 18 + 1] ), &( allcands[ii * 18 + 2] ), &( allcands[ii * 18 + 3] ), &( allcands[ii * 18 + 4] ), &( allcands[ii * 18 + 5] ), &( allcands[ii * 18 + 6] ), &( allcands[ii * 18 + 7] ), &( allcands[ii * 18 + 8] ), &( allcands_job[2 * ii] ), &( allcands[ii * 18 + 9] ), &( allcands[ii * 18 + 10] ), &( allcands[ii * 18 + 11] ), &( allcands[ii * 18 + 12] ), &( allcands[ii * 18 + 13] ), &( allcands[ii * 18 + 14] ), &( allcands[ii * 18 + 15] ), &( allcands[ii * 18 + 16] ), &( allcands[ii * 18 + 17] ), &( allcands_job[ii * 2 + 1] ) );
107 }
108 fclose( CANDS );
109 ( *output ) = allcands;
110 ( *output_job ) = allcands_job;
111 ( *output_numCoincident ) = count;
112 return XLAL_SUCCESS;
113}
114INT4 readInIFOoutliersAndSort( double **output, int **output_job, int *output_numOutliers, const char *infilename )
115{
116 FILE *IFOCANDS = NULL;
117 XLAL_CHECK( ( IFOCANDS = fopen( infilename, "r" ) ) != NULL, XLAL_EIO, "Can't fopen %s", infilename );
118 //Determines number of candidates in the file
119 INT4 ch, ifocount = 0;
120 do {
121 ch = fgetc( IFOCANDS );
122 if ( ch == '\n' ) {
123 ifocount++;
124 }
125 } while ( ch != EOF );
126 double *allIFOcands = NULL;
127 XLAL_CHECK( ( allIFOcands = ( double * )XLALMalloc( sizeof( double ) * ifocount * 9 ) ) != NULL, XLAL_ENOMEM );
128 int *allIFOcands_job = NULL;
129 XLAL_CHECK( ( allIFOcands_job = ( int * )XLALMalloc( sizeof( int ) * ifocount ) ) != NULL, XLAL_ENOMEM );
130 rewind( IFOCANDS );
131 //Put data in array
132 for ( INT4 ii = 0; ii < ifocount; ii++ ) {
133 fscanf( IFOCANDS, "%la %la %la %la %la %la %la %la %la %d", &( allIFOcands[ii * 9] ), &( allIFOcands[ii * 9 + 1] ), &( allIFOcands[ii * 9 + 2] ), &( allIFOcands[ii * 9 + 3] ), &( allIFOcands[ii * 9 + 4] ), &( allIFOcands[ii * 9 + 5] ), &( allIFOcands[ii * 9 + 6] ), &( allIFOcands[ii * 9 + 7] ), &( allIFOcands[ii * 9 + 8] ), &( allIFOcands_job[ii] ) );
134 }
135 fclose( IFOCANDS );
136 //Sort the array based on the frequency
137 size_t *sorted_index = NULL;
138 XLAL_CHECK( ( sorted_index = ( size_t * )XLALMalloc( sizeof( size_t ) * ifocount ) ) != NULL, XLAL_ENOMEM );
139 double *allIFOcands_sorted = NULL;
140 XLAL_CHECK( ( allIFOcands_sorted = ( double * )XLALMalloc( sizeof( double ) * ifocount * 9 ) ) != NULL, XLAL_ENOMEM );
141 int *allIFOcands_job_sorted = NULL;
142 XLAL_CHECK( ( allIFOcands_job_sorted = ( int * )XLALMalloc( sizeof( int ) * ifocount ) ) != NULL, XLAL_ENOMEM );
143 gsl_sort_index( sorted_index, allIFOcands, 9, ifocount );
144 for ( INT4 ii = 0; ii < ifocount; ii++ ) {
145 memcpy( &( allIFOcands_sorted[ii * 9] ), &( allIFOcands[sorted_index[ii] * 9] ), sizeof( double ) * 9 );
146 allIFOcands_job_sorted[ii] = allIFOcands_job[sorted_index[ii]];
147 }
148 XLALFree( allIFOcands );
149 XLALFree( sorted_index );
150 XLALFree( allIFOcands_job );
151 ( *output ) = allIFOcands_sorted;
152 ( *output_job ) = allIFOcands_job_sorted;
153 ( *output_numOutliers ) = ifocount;
154 return XLAL_SUCCESS;
155}
156
157INT4 main( int argc, char *argv[] )
158{
159
160 //Turn off gsl error handler
161 gsl_set_error_handler_off();
162
163 UserVariables_t XLAL_INIT_DECL( uvar );
164 XLAL_CHECK( InitUserVars( &uvar, argc, argv ) == XLAL_SUCCESS, XLAL_EFUNC );
165
166 double *allIFO1cands_sorted = NULL, *allIFO2cands_sorted = NULL;
167 int *allIFO1cands_job_sorted = NULL, *allIFO2cands_job_sorted = NULL;
168 int ifo1count = 0, ifo2count = 0;
169 XLAL_CHECK( readInIFOoutliersAndSort( &allIFO1cands_sorted, &allIFO1cands_job_sorted, &ifo1count, uvar.infile1 ) == XLAL_SUCCESS, XLAL_EFUNC );
170 XLAL_CHECK( readInIFOoutliersAndSort( &allIFO2cands_sorted, &allIFO2cands_job_sorted, &ifo2count, uvar.infile2 ) == XLAL_SUCCESS, XLAL_EFUNC );
171
172 //Open a file to save the output data
173 FILE *CANDS = NULL;
174 XLAL_CHECK( ( CANDS = fopen( uvar.outfile1, "w" ) ) != NULL, XLAL_EIO, "Can't fopen %s", uvar.outfile1 );
175
176 //Setup and allocate the solver
177 int status = GSL_CONTINUE;
178 int max_iter = 100;
179 double x_lo = 0.0, x_hi = ( double )( ifo2count - 1 );
180 double foundIndex = -1.0;
181 gsl_function F;
182 F.function = &f_diff;
183 const gsl_root_fsolver_type *T = gsl_root_fsolver_brent;
184 gsl_root_fsolver *s = NULL;
185 XLAL_CHECK( ( s = gsl_root_fsolver_alloc( T ) ) != NULL, XLAL_EFUNC );
186
187 double tobs = uvar.Tobs;
188 double fdiff_allowed = uvar.fdiff_allowed;
189 double dfdiff_allowed = uvar.dfdiff_allowed;
190 double skydiff_allowed = uvar.skydiff_allowed * 200.0;
191 double periodTcohFactor = 2.7 * ( uvar.Tcoh / 1800.0 ) + 1.8;
192 int numpassingf = 0, numpassingdf = 0, numpassingP = 0, numpassingskyloc = 0;
193 INT4 ii, jj;
194 for ( ii = 0; ii < ifo1count; ii++ ) {
195
196 //Check that the frequency of the H1 candidate we look at is within the frequency span of the L1 candidates
197 if ( allIFO1cands_sorted[ii * 9] >= allIFO2cands_sorted[0] && allIFO1cands_sorted[ii * 9] <= allIFO2cands_sorted[( ifo2count - 1 ) * 9] ) {
198 if ( foundIndex < 0.0 || status != GSL_SUCCESS || ( status == GSL_SUCCESS && fabs( allIFO1cands_sorted[ii * 9] - allIFO2cands_sorted[( int )round( gsl_root_fsolver_root( s ) ) * 9] ) > fdiff_allowed ) ) {
199 //Do a root finding search for the closest L1 candidate in frequency to the H1 candidate frequency
200 int iter = 0;
201 struct solver_params params = {allIFO1cands_sorted[ii * 9], allIFO2cands_sorted};
202 F.params = &params;
203 XLAL_CHECK( gsl_root_fsolver_set( s, &F, x_lo, x_hi ) == GSL_SUCCESS, XLAL_EFUNC );
204 do {
205 iter++;
206 status = gsl_root_fsolver_iterate( s );
207 XLAL_CHECK( status == GSL_SUCCESS, XLAL_EFUNC );
208 foundIndex = gsl_root_fsolver_root( s );
209 status = gsl_root_test_residual( f_diff( foundIndex, &params ), 1.05 * fdiff_allowed );
210 XLAL_CHECK( status == GSL_SUCCESS || status == GSL_CONTINUE, XLAL_EFUNC );
211 } while ( status == GSL_CONTINUE && iter < max_iter );
212 }
213
214 //If the search was successful, then we step through the L1 candidates to find matching candidates
215 if ( status == GSL_SUCCESS ) {
216 jj = ( int )round( foundIndex ); //start at the index of the L1 candidate found in the root finding
217
218 //Step backwards in L1 candidates until we are definitely below the H1 candidate in frequency (or at the start of the L1 list)
219 while ( jj > 0 && ( allIFO1cands_sorted[ii * 9] - allIFO2cands_sorted[jj * 9] ) < 1.05 * fdiff_allowed ) {
220 jj--;
221 }
222
223 //Set the foundIndex value to the start of where we should search from
224 foundIndex = jj;
225
226 //Starting from the L1 candidate below the H1 candidate frequency
227 for ( ; jj < ifo2count; jj++ ) {
228 //Check that if the frequency of L1 candidate is above the H1 value by greater than the allowed value, break the loop
229 if ( allIFO2cands_sorted[jj * 9] - allIFO1cands_sorted[ii * 9] > 1.05 * fdiff_allowed ) {
230 break;
231 }
232
233 //If the H1 and L1 frequency values are near enough, proceed with checking more parameter values
234 if ( fabs( allIFO1cands_sorted[ii * 9] - allIFO2cands_sorted[jj * 9] ) <= fdiff_allowed ) {
235 numpassingf++;
236 //Check the modulation depth
237 if ( fabs( allIFO1cands_sorted[ii * 9 + 2] - allIFO2cands_sorted[jj * 9 + 2] ) <= dfdiff_allowed ) {
238 numpassingdf++;
239 //Check the period and harmonic values
240 double Pdiff_allowed = allIFO1cands_sorted[ii * 9 + 1] * allIFO1cands_sorted[ii * 9 + 1] * sqrt( 3.6e-3 / allIFO1cands_sorted[ii * 9 + 2] ) / ( periodTcohFactor * tobs );
241 double Pdiff_allowed_2 = allIFO2cands_sorted[jj * 9 + 1] * allIFO2cands_sorted[jj * 9 + 1] * sqrt( 3.6e-3 / allIFO2cands_sorted[jj * 9 + 2] ) / ( periodTcohFactor * tobs );
242 int foundmatch = 0, passedtestP = 0;
243 for ( int kk = 1; kk <= 7; kk++ ) {
244 double P1factor = 0.0;
245 if ( kk == 1 ) {
246 P1factor = 1.0;
247 } else if ( kk < 5 ) {
248 P1factor = 1.0 / kk;
249 } else {
250 P1factor = ( double )( kk - 3 );
251 }
252
253 for ( int ll = 1; ll <= 7; ll++ ) {
254 double P2factor = 0.0;
255 if ( ll == 1 ) {
256 P2factor = 1.0;
257 } else if ( ll < 5 ) {
258 P2factor = 1.0 / ll;
259 } else {
260 P2factor = ( double )( ll - 3 );
261 }
262
263 if ( fabs( P1factor * allIFO1cands_sorted[ii * 9 + 1] - P2factor * allIFO2cands_sorted[jj * 9 + 1] ) <= Pdiff_allowed * ( P1factor * P1factor ) || fabs( P1factor * allIFO1cands_sorted[ii * 9 + 1] - P2factor * allIFO2cands_sorted[jj * 9 + 1] ) <= Pdiff_allowed_2 * ( P2factor * P2factor ) ) {
264 if ( !passedtestP ) {
265 numpassingP++;
266 passedtestP = 1;
267 }
268 //Check the sky location
269 double absd1mPo2 = fabs( allIFO1cands_sorted[ii * 9 + 4] - M_PI_2 );
270 double absd2mPo2 = fabs( allIFO2cands_sorted[jj * 9 + 4] - M_PI_2 );
271 double dist = acos( sin( absd1mPo2 ) * sin( absd2mPo2 ) * cos( allIFO1cands_sorted[ii * 9 + 3] - allIFO2cands_sorted[jj * 9 + 3] ) + cos( absd1mPo2 ) * cos( absd2mPo2 ) );
272
273 if ( dist <= 2.0 * skydiff_allowed / ( allIFO1cands_sorted[ii * 9] + allIFO2cands_sorted[jj * 9] ) ) {
274 foundmatch = 1;
275 numpassingskyloc++;
276 fprintf( CANDS, "%f %f %f %f %f %f %g %f %g %d %f %f %f %f %f %f %g %f %g %d\n", allIFO1cands_sorted[ii * 9], allIFO1cands_sorted[ii * 9 + 1], allIFO1cands_sorted[ii * 9 + 2], allIFO1cands_sorted[ii * 9 + 3], allIFO1cands_sorted[ii * 9 + 4], allIFO1cands_sorted[ii * 9 + 5], allIFO1cands_sorted[ii * 9 + 6], allIFO1cands_sorted[ii * 9 + 7], allIFO1cands_sorted[ii * 9 + 8], allIFO1cands_job_sorted[ii], allIFO2cands_sorted[jj * 9], allIFO2cands_sorted[jj * 9 + 1], allIFO2cands_sorted[jj * 9 + 2], allIFO2cands_sorted[jj * 9 + 3], allIFO2cands_sorted[jj * 9 + 4], allIFO2cands_sorted[jj * 9 + 5], allIFO2cands_sorted[jj * 9 + 6], allIFO2cands_sorted[jj * 9 + 7], allIFO2cands_sorted[jj * 9 + 8], allIFO2cands_job_sorted[jj] );
277 } //end sky check
278 } //end period check
279 if ( foundmatch ) {
280 break;
281 }
282 } //end test different P2factors
283 if ( foundmatch ) {
284 break;
285 }
286 } //end test different P1factors
287 } //end modulation depth check
288 } //end frequency check
289 } //end test against L1 values
290 } //end successful search
291 } else if ( ( allIFO2cands_sorted[0] - allIFO1cands_sorted[ii * 9] ) <= fdiff_allowed && allIFO1cands_sorted[ii * 9] <= allIFO2cands_sorted[0] ) {
292
293 for ( jj = 0; jj < ifo2count; jj++ ) {
294 if ( allIFO2cands_sorted[jj * 9] - allIFO1cands_sorted[ii * 9] > 1.05 * fdiff_allowed ) {
295 break;
296 }
297
298 if ( fabs( allIFO1cands_sorted[ii * 9] - allIFO2cands_sorted[jj * 9] ) <= fdiff_allowed ) {
299 numpassingf++;
300 if ( fabs( allIFO1cands_sorted[ii * 9 + 2] - allIFO2cands_sorted[jj * 9 + 2] ) <= dfdiff_allowed ) {
301 numpassingdf++;
302 double Pdiff_allowed = allIFO1cands_sorted[ii * 9 + 1] * allIFO1cands_sorted[ii * 9 + 1] * sqrt( 3.6e-3 / allIFO1cands_sorted[ii * 9 + 2] ) / ( periodTcohFactor * tobs );
303 double Pdiff_allowed_2 = allIFO2cands_sorted[jj * 9 + 1] * allIFO2cands_sorted[jj * 9 + 1] * sqrt( 3.6e-3 / allIFO2cands_sorted[jj * 9 + 2] ) / ( periodTcohFactor * tobs );
304 int foundmatch = 0, passedtestP = 0;
305 for ( int kk = 1; kk <= 7; kk++ ) {
306 double P1factor = 0.0;
307 if ( kk == 1 ) {
308 P1factor = 1.0;
309 } else if ( kk < 5 ) {
310 P1factor = 1.0 / kk;
311 } else {
312 P1factor = ( double )( kk - 3 );
313 }
314
315 for ( int ll = 1; ll <= 7; ll++ ) {
316 double P2factor = 0.0;
317 if ( ll == 1 ) {
318 P2factor = 1.0;
319 } else if ( ll < 5 ) {
320 P2factor = 1.0 / ll;
321 } else {
322 P2factor = ( double )( ll - 3 );
323 }
324
325 if ( fabs( P1factor * allIFO1cands_sorted[ii * 9 + 1] - P2factor * allIFO2cands_sorted[jj * 9 + 1] ) <= Pdiff_allowed * ( P1factor * P1factor ) || fabs( P1factor * allIFO1cands_sorted[ii * 9 + 1] - P2factor * allIFO2cands_sorted[jj * 9 + 1] ) <= Pdiff_allowed_2 * ( P2factor * P2factor ) ) {
326 if ( !passedtestP ) {
327 numpassingP++;
328 passedtestP = 1;
329 }
330 //Check the sky location
331 double absd1mPo2 = fabs( allIFO1cands_sorted[ii * 9 + 4] - M_PI_2 );
332 double absd2mPo2 = fabs( allIFO2cands_sorted[jj * 9 + 4] - M_PI_2 );
333 double dist = acos( sin( absd1mPo2 ) * sin( absd2mPo2 ) * cos( allIFO1cands_sorted[ii * 9 + 3] - allIFO2cands_sorted[jj * 9 + 3] ) + cos( absd1mPo2 ) * cos( absd2mPo2 ) );
334
335 if ( dist <= 2.0 * skydiff_allowed / ( allIFO1cands_sorted[ii * 9] + allIFO2cands_sorted[jj * 9] ) ) {
336 foundmatch = 1;
337 numpassingskyloc++;
338 fprintf( CANDS, "%f %f %f %f %f %f %g %f %g %d %f %f %f %f %f %f %g %f %g %d\n", allIFO1cands_sorted[ii * 9], allIFO1cands_sorted[ii * 9 + 1], allIFO1cands_sorted[ii * 9 + 2], allIFO1cands_sorted[ii * 9 + 3], allIFO1cands_sorted[ii * 9 + 4], allIFO1cands_sorted[ii * 9 + 5], allIFO1cands_sorted[ii * 9 + 6], allIFO1cands_sorted[ii * 9 + 7], allIFO1cands_sorted[ii * 9 + 8], allIFO1cands_job_sorted[ii], allIFO2cands_sorted[jj * 9], allIFO2cands_sorted[jj * 9 + 1], allIFO2cands_sorted[jj * 9 + 2], allIFO2cands_sorted[jj * 9 + 3], allIFO2cands_sorted[jj * 9 + 4], allIFO2cands_sorted[jj * 9 + 5], allIFO2cands_sorted[jj * 9 + 6], allIFO2cands_sorted[jj * 9 + 7], allIFO2cands_sorted[jj * 9 + 8], allIFO2cands_job_sorted[jj] );
339 } //end sky check
340 } //end period check
341 if ( foundmatch ) {
342 break;
343 }
344 } //end test different P2factors
345 if ( foundmatch ) {
346 break;
347 }
348 } //end test different P1factors
349 } //end modulation depth check
350 } //end frequency check
351 } //end test against L1 values
352 } else if ( ( allIFO1cands_sorted[ii * 9] - allIFO2cands_sorted[( ifo2count - 1 ) * 9] ) <= fdiff_allowed && allIFO1cands_sorted[ii * 9] >= allIFO2cands_sorted[( ifo2count - 1 ) * 9] ) {
353 jj = ifo2count - 1;
354 while ( jj > 0 && ( allIFO1cands_sorted[ii * 9] - allIFO2cands_sorted[jj * 9] ) < 1.05 * fdiff_allowed ) {
355 jj--;
356 }
357
358 for ( ; jj < ifo2count; jj++ ) {
359 if ( fabs( allIFO1cands_sorted[ii * 9] - allIFO2cands_sorted[jj * 9] ) <= fdiff_allowed ) {
360 numpassingf++;
361 if ( fabs( allIFO1cands_sorted[ii * 9 + 2] - allIFO2cands_sorted[jj * 9 + 2] ) <= dfdiff_allowed ) {
362 numpassingdf++;
363 double Pdiff_allowed = allIFO1cands_sorted[ii * 9 + 1] * allIFO1cands_sorted[ii * 9 + 1] * sqrt( 3.6e-3 / allIFO1cands_sorted[ii * 9 + 2] ) / ( periodTcohFactor * tobs );
364 double Pdiff_allowed_2 = allIFO2cands_sorted[jj * 9 + 1] * allIFO2cands_sorted[jj * 9 + 1] * sqrt( 3.6e-3 / allIFO2cands_sorted[jj * 9 + 2] ) / ( periodTcohFactor * tobs );
365 int foundmatch = 0, passedtestP = 0;
366 for ( int kk = 1; kk <= 7; kk++ ) {
367 double P1factor = 0.0;
368 if ( kk == 1 ) {
369 P1factor = 1.0;
370 } else if ( kk < 5 ) {
371 P1factor = 1.0 / kk;
372 } else {
373 P1factor = ( double )( kk - 3 );
374 }
375
376 for ( int ll = 1; ll <= 7; ll++ ) {
377 double P2factor = 0.0;
378 if ( ll == 1 ) {
379 P2factor = 1.0;
380 } else if ( ll < 5 ) {
381 P2factor = 1.0 / ll;
382 } else {
383 P2factor = ( double )( ll - 3 );
384 }
385
386 if ( fabs( P1factor * allIFO1cands_sorted[ii * 9 + 1] - P2factor * allIFO2cands_sorted[jj * 9 + 1] ) <= Pdiff_allowed * ( P1factor * P1factor ) || fabs( P1factor * allIFO1cands_sorted[ii * 9 + 1] - P2factor * allIFO2cands_sorted[jj * 9 + 1] ) <= Pdiff_allowed_2 * ( P2factor * P2factor ) ) {
387 if ( !passedtestP ) {
388 numpassingP++;
389 passedtestP = 1;
390 }
391 //Check the sky location
392 double absd1mPo2 = fabs( allIFO1cands_sorted[ii * 9 + 4] - M_PI_2 );
393 double absd2mPo2 = fabs( allIFO2cands_sorted[jj * 9 + 4] - M_PI_2 );
394 double dist = acos( sin( absd1mPo2 ) * sin( absd2mPo2 ) * cos( allIFO1cands_sorted[ii * 9 + 3] - allIFO2cands_sorted[jj * 9 + 3] ) + cos( absd1mPo2 ) * cos( absd2mPo2 ) );
395
396 if ( dist <= 2.0 * skydiff_allowed / ( allIFO1cands_sorted[ii * 9] + allIFO2cands_sorted[jj * 9] ) ) {
397 foundmatch = 1;
398 numpassingskyloc++;
399 fprintf( CANDS, "%f %f %f %f %f %f %g %f %g %d %f %f %f %f %f %f %g %f %g %d\n", allIFO1cands_sorted[ii * 9], allIFO1cands_sorted[ii * 9 + 1], allIFO1cands_sorted[ii * 9 + 2], allIFO1cands_sorted[ii * 9 + 3], allIFO1cands_sorted[ii * 9 + 4], allIFO1cands_sorted[ii * 9 + 5], allIFO1cands_sorted[ii * 9 + 6], allIFO1cands_sorted[ii * 9 + 7], allIFO1cands_sorted[ii * 9 + 8], allIFO1cands_job_sorted[ii], allIFO2cands_sorted[jj * 9], allIFO2cands_sorted[jj * 9 + 1], allIFO2cands_sorted[jj * 9 + 2], allIFO2cands_sorted[jj * 9 + 3], allIFO2cands_sorted[jj * 9 + 4], allIFO2cands_sorted[jj * 9 + 5], allIFO2cands_sorted[jj * 9 + 6], allIFO2cands_sorted[jj * 9 + 7], allIFO2cands_sorted[jj * 9 + 8], allIFO2cands_job_sorted[jj] );
400 } //end sky check
401 } //end period check
402 if ( foundmatch ) {
403 break;
404 }
405 } //end test different P2factors
406 if ( foundmatch ) {
407 break;
408 }
409 } //end test different P1factors
410 } //end modulation depth check
411 } //end frequency check
412 } //end check against L1 values
413 } //end if H1 candidate is barely outside L1 frequency range
414 } //end loop over H1 values
415
416 //Close combined list file
417 fclose( CANDS );
418 CANDS = NULL;
419
420 gsl_root_fsolver_free( s );
421 XLALFree( allIFO1cands_sorted );
422 XLALFree( allIFO2cands_sorted );
423 XLALFree( allIFO1cands_job_sorted );
424 XLALFree( allIFO2cands_job_sorted );
425
426 fprintf( stderr, "Passed f = %d, passed df = %d, passed P = %d, passed sky loc %d\n", numpassingf, numpassingdf, numpassingP, numpassingskyloc );
427
428
429 /// PART TWO: ///
430
431 double *allcands = NULL;
432 int *allcands_job = NULL;
433 int count = 0;
434 XLAL_CHECK( readInCoincidentOutliers( &allcands, &allcands_job, &count, uvar.outfile1 ) == XLAL_SUCCESS, XLAL_EFUNC );
435
436 //Allocate usedvalue array
437 INT4Vector *usedvalue = NULL;
438 XLAL_CHECK( ( usedvalue = XLALCreateINT4Vector( count ) ) != NULL, XLAL_EFUNC );
439 memset( usedvalue->data, 0, sizeof( INT4 )*count );
440
441 //Open a file to save the output data
442 FILE *NEWCANDS = NULL;
443 XLAL_CHECK( ( NEWCANDS = fopen( uvar.outfile2, "w" ) ) != NULL, XLAL_EIO, "Couldn't fopen %s\n", uvar.outfile2 );
444
445 for ( ii = 0; ii < count; ii++ ) {
446 if ( !usedvalue->data[ii] ) {
447 int bestcand = ii;
448 double bestcandprob = allcands[ii * 18 + 16];
449 for ( jj = 0; jj < count; jj++ ) {
450 if ( usedvalue->data[jj] || jj == ii ) {
451 continue;
452 }
453 if ( allcands[jj * 18] == allcands[bestcand * 18] && allcands[jj * 18 + 1] == allcands[bestcand * 18 + 1] && allcands[jj * 18 + 2] == allcands[bestcand * 18 + 2] && allcands[jj * 18 + 3] == allcands[bestcand * 18 + 3] && allcands[jj * 18 + 4] == allcands[bestcand * 18 + 4] ) {
454 if ( allcands[jj * 18 + 16] < bestcandprob ) {
455 usedvalue->data[bestcand] = 1;
456 bestcandprob = allcands[jj * 18 + 16];
457 bestcand = jj;
458 } else {
459 usedvalue->data[jj] = 1;
460 }
461 }
462 }
463 fprintf( NEWCANDS, "%f %f %f %f %f %f %g %f %g %d %f %f %f %f %f %f %g %f %g %d\n", allcands[bestcand * 18], allcands[bestcand * 18 + 1], allcands[bestcand * 18 + 2], allcands[bestcand * 18 + 3], allcands[bestcand * 18 + 4], allcands[bestcand * 18 + 5], allcands[bestcand * 18 + 6], allcands[bestcand * 18 + 7], allcands[bestcand * 18 + 8], allcands_job[bestcand * 2], allcands[bestcand * 18 + 9], allcands[bestcand * 18 + 10], allcands[bestcand * 18 + 11], allcands[bestcand * 18 + 12], allcands[bestcand * 18 + 13], allcands[bestcand * 18 + 14], allcands[bestcand * 18 + 15], allcands[bestcand * 18 + 16], allcands[bestcand * 18 + 17], allcands_job[bestcand * 2 + 1] );
464 usedvalue->data[bestcand] = 1;
465 }
466 }
467
468 fclose( NEWCANDS );
469 NEWCANDS = NULL;
470
471 XLALFree( allcands );
472 XLALFree( allcands_job );
473 XLALDestroyINT4Vector( usedvalue );
474
475 /// PART THREE: ///
476
477 XLAL_CHECK( readInCoincidentOutliers( &allcands, &allcands_job, &count, uvar.outfile2 ) == XLAL_SUCCESS, XLAL_EFUNC );
478
479 XLAL_CHECK( ( usedvalue = XLALCreateINT4Vector( count ) ) != NULL, XLAL_EFUNC );
480 memset( usedvalue->data, 0, sizeof( INT4 )*count );
481
482 //Open a file to save the output data
483 XLAL_CHECK( ( NEWCANDS = fopen( uvar.finalOutfile, "w" ) ) != NULL, XLAL_EIO, "Couldn't fopen %s\n", uvar.finalOutfile );
484
485 for ( ii = 0; ii < count; ii++ ) {
486 if ( !usedvalue->data[ii] ) {
487 int bestcand = ii;
488 double bestcandprob = allcands[ii * 18 + 7];
489 for ( jj = 0; jj < count; jj++ ) {
490 if ( usedvalue->data[jj] || jj == ii ) {
491 continue;
492 }
493 if ( allcands[jj * 18 + 9] == allcands[bestcand * 18 + 9] && allcands[jj * 18 + 10] == allcands[bestcand * 18 + 10] && allcands[jj * 18 + 11] == allcands[bestcand * 18 + 11] && allcands[jj * 18 + 12] == allcands[bestcand * 18 + 12] && allcands[jj * 18 + 13] == allcands[bestcand * 18 + 13] ) {
494 if ( allcands[jj * 18 + 7] < bestcandprob ) {
495 usedvalue->data[bestcand] = 1;
496 bestcandprob = allcands[jj * 18 + 7];
497 bestcand = jj;
498 } else {
499 usedvalue->data[jj] = 1;
500 }
501 }
502 }
503 fprintf( NEWCANDS, "%f %f %f %f %f %f %g %f %g %d %f %f %f %f %f %f %g %f %g %d\n", allcands[bestcand * 18], allcands[bestcand * 18 + 1], allcands[bestcand * 18 + 2], allcands[bestcand * 18 + 3], allcands[bestcand * 18 + 4], allcands[bestcand * 18 + 5], allcands[bestcand * 18 + 6], allcands[bestcand * 18 + 7], allcands[bestcand * 18 + 8], allcands_job[bestcand * 2], allcands[bestcand * 18 + 9], allcands[bestcand * 18 + 10], allcands[bestcand * 18 + 11], allcands[bestcand * 18 + 12], allcands[bestcand * 18 + 13], allcands[bestcand * 18 + 14], allcands[bestcand * 18 + 15], allcands[bestcand * 18 + 16], allcands[bestcand * 18 + 17], allcands_job[bestcand * 2 + 1] );
504 usedvalue->data[bestcand] = 1;
505 }
506 }
507
508 fclose( NEWCANDS );
509
510 XLALFree( allcands );
511 XLALFree( allcands_job );
512 XLALDestroyINT4Vector( usedvalue );
513
514
515 //PART FOUR: Swap input values
516 XLAL_CHECK( readInCoincidentOutliers( &allcands, &allcands_job, &count, uvar.outfile1 ) == XLAL_SUCCESS, XLAL_EFUNC );
517
518 XLAL_CHECK( ( usedvalue = XLALCreateINT4Vector( count ) ) != NULL, XLAL_EFUNC );
519 memset( usedvalue->data, 0, sizeof( INT4 )*count );
520
521 //Open a file to save the output data
522 XLAL_CHECK( ( NEWCANDS = fopen( uvar.outfile3, "w" ) ) != NULL, XLAL_EIO, "Couldn't fopen %s\n", uvar.outfile3 );
523
524 for ( ii = 0; ii < count; ii++ ) {
525 if ( !usedvalue->data[ii] ) {
526 int bestcand = ii;
527 double bestcandprob = allcands[ii * 18 + 7];
528 for ( jj = 0; jj < count; jj++ ) {
529 if ( usedvalue->data[jj] || jj == ii ) {
530 continue;
531 }
532 if ( allcands[jj * 18 + 9] == allcands[bestcand * 18 + 9] && allcands[jj * 18 + 10] == allcands[bestcand * 18 + 10] && allcands[jj * 18 + 11] == allcands[bestcand * 18 + 11] && allcands[jj * 18 + 12] == allcands[bestcand * 18 + 12] && allcands[jj * 18 + 13] == allcands[bestcand * 18 + 13] ) {
533 if ( allcands[jj * 18 + 7] < bestcandprob ) {
534 usedvalue->data[bestcand] = 1;
535 bestcandprob = allcands[jj * 18 + 7];
536 bestcand = jj;
537 } else {
538 usedvalue->data[jj] = 1;
539 }
540 }
541 }
542 fprintf( NEWCANDS, "%f %f %f %f %f %f %g %f %g %d %f %f %f %f %f %f %g %f %g %d\n", allcands[bestcand * 18], allcands[bestcand * 18 + 1], allcands[bestcand * 18 + 2], allcands[bestcand * 18 + 3], allcands[bestcand * 18 + 4], allcands[bestcand * 18 + 5], allcands[bestcand * 18 + 6], allcands[bestcand * 18 + 7], allcands[bestcand * 18 + 8], allcands_job[bestcand * 2], allcands[bestcand * 18 + 9], allcands[bestcand * 18 + 10], allcands[bestcand * 18 + 11], allcands[bestcand * 18 + 12], allcands[bestcand * 18 + 13], allcands[bestcand * 18 + 14], allcands[bestcand * 18 + 15], allcands[bestcand * 18 + 16], allcands[bestcand * 18 + 17], allcands_job[bestcand * 2 + 1] );
543 usedvalue->data[bestcand] = 1;
544 }
545 }
546
547 fclose( NEWCANDS );
548 NEWCANDS = NULL;
549
550 XLALFree( allcands );
551 XLALFree( allcands_job );
552 XLALDestroyINT4Vector( usedvalue );
553
554
555 XLAL_CHECK( readInCoincidentOutliers( &allcands, &allcands_job, &count, uvar.outfile3 ) == XLAL_SUCCESS, XLAL_EFUNC );
556
557 XLAL_CHECK( ( usedvalue = XLALCreateINT4Vector( count ) ) != NULL, XLAL_EFUNC );
558 memset( usedvalue->data, 0, sizeof( INT4 )*count );
559
560 //Open a file to save the output data
561 XLAL_CHECK( ( NEWCANDS = fopen( uvar.finalOutfile, "a" ) ) != NULL, XLAL_EIO, "Couldn't fopen %s\n", uvar.finalOutfile );
562
563 for ( ii = 0; ii < count; ii++ ) {
564 if ( !usedvalue->data[ii] ) {
565 int bestcand = ii;
566 double bestcandprob = allcands[ii * 18 + 16];
567 for ( jj = 0; jj < count; jj++ ) {
568 if ( usedvalue->data[jj] || jj == ii ) {
569 continue;
570 }
571 if ( allcands[jj * 18] == allcands[bestcand * 18] && allcands[jj * 18 + 1] == allcands[bestcand * 18 + 1] && allcands[jj * 18 + 2] == allcands[bestcand * 18 + 2] && allcands[jj * 18 + 3] == allcands[bestcand * 18 + 3] && allcands[jj * 18 + 4] == allcands[bestcand * 18 + 4] ) {
572 if ( allcands[jj * 18 + 16] < bestcandprob ) {
573 usedvalue->data[bestcand] = 1;
574 bestcandprob = allcands[jj * 18 + 16];
575 bestcand = jj;
576 } else {
577 usedvalue->data[jj] = 1;
578 }
579 }
580 }
581 fprintf( NEWCANDS, "%f %f %f %f %f %f %g %f %g %d %f %f %f %f %f %f %g %f %g %d\n", allcands[bestcand * 18], allcands[bestcand * 18 + 1], allcands[bestcand * 18 + 2], allcands[bestcand * 18 + 3], allcands[bestcand * 18 + 4], allcands[bestcand * 18 + 5], allcands[bestcand * 18 + 6], allcands[bestcand * 18 + 7], allcands[bestcand * 18 + 8], allcands_job[bestcand * 2], allcands[bestcand * 18 + 9], allcands[bestcand * 18 + 10], allcands[bestcand * 18 + 11], allcands[bestcand * 18 + 12], allcands[bestcand * 18 + 13], allcands[bestcand * 18 + 14], allcands[bestcand * 18 + 15], allcands[bestcand * 18 + 16], allcands[bestcand * 18 + 17], allcands_job[bestcand * 2 + 1] );
582 usedvalue->data[bestcand] = 1;
583 }
584 }
585
586 fclose( NEWCANDS );
587 NEWCANDS = NULL;
588
589 XLALFree( allcands );
590 XLALFree( allcands_job );
591 XLALDestroyINT4Vector( usedvalue );
592
593 XLALDestroyUserVars();
594
595 return 0;
596}
lalPulsarVCSInfoList
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
STRING
#define STRING(a)
fscanf
#define fscanf
fprintf
#define fprintf
f_diff
REAL8 f_diff(double indexval, void *params)
Definition: compareCandidates.c:82
readInCoincidentOutliers
INT4 readInCoincidentOutliers(double **output, int **output_job, int *output_numOutliers, const char *infilename)
Definition: compareCandidates.c:87
readInIFOoutliersAndSort
INT4 readInIFOoutliersAndSort(double **output, int **output_job, int *output_numCoincident, const char *infilename)
Definition: compareCandidates.c:114
InitUserVars
INT4 InitUserVars(UserVariables_t *uvar, int argc, char *argv[])
Definition: compareCandidates.c:55
main
INT4 main(int argc, char *argv[])
Definition: compareCandidates.c:157
BOOLEAN
unsigned char BOOLEAN
REAL8
double REAL8
XLAL_INIT_DECL
#define XLAL_INIT_DECL(var,...)
CHAR
char CHAR
INT4
int32_t INT4
XLALMalloc
void * XLALMalloc(size_t n)
XLALFree
void XLALFree(void *p)
XLALUserVarReadAllInput
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
XLALDestroyUserVars
void XLALDestroyUserVars(void)
XLALRegisterUvarMember
#define XLALRegisterUvarMember(name, type, option, category,...)
XLALCreateINT4Vector
INT4Vector * XLALCreateINT4Vector(UINT4 length)
XLALDestroyINT4Vector
void XLALDestroyINT4Vector(INT4Vector *vector)
XLAL_CHECK
#define XLAL_CHECK(assertion,...)
XLAL_ENOMEM
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLAL_EIO
XLAL_EIO
XLAL_EINVAL
XLAL_EINVAL
count
long long count
Definition: hwinject.c:371
T
T
infilename
def infilename
dist
dist
lalpulsar.git_version.status
string status
Definition: git_version.py:32
lalpulsar_MakeSFTDAG.int
int
Definition: lalpulsar_MakeSFTDAG.py:413
lalpulsar_PiecewiseSearch.s
s
Definition: lalpulsar_PiecewiseSearch.py:957
lalpulsar_knope_automation_script.p
p
RUN ANALYSIS SCRIPT ###.
Definition: lalpulsar_knope_automation_script.py:709
make_frame_cache.output
output
Definition: make_frame_cache.py:53
INT4Vector::data
INT4 * data
UserVariables_t
user input variables
Definition: compareFstats.c:51
UserVariables_t::Tobs
REAL8 Tobs
Definition: compareCandidates.c:38
UserVariables_t::fdiff_allowed
REAL8 fdiff_allowed
Definition: compareCandidates.c:40
UserVariables_t::infile1
CHAR * infile1
Definition: compareCandidates.c:32
UserVariables_t::outfile3
CHAR * outfile3
Definition: compareCandidates.c:36
UserVariables_t::skydiff_allowed
REAL8 skydiff_allowed
Definition: compareCandidates.c:42
UserVariables_t::dfdiff_allowed
REAL8 dfdiff_allowed
Definition: compareCandidates.c:41
UserVariables_t::outfile2
CHAR * outfile2
Definition: compareCandidates.c:35
UserVariables_t::Tcoh
REAL8 Tcoh
Definition: compareCandidates.c:39
UserVariables_t::outfile1
CHAR * outfile1
Definition: compareCandidates.c:34
UserVariables_t::finalOutfile
CHAR * finalOutfile
Definition: compareCandidates.c:37
UserVariables_t::infile2
CHAR * infile2
Definition: compareCandidates.c:33
params::params
LALDict * params
solver_params
Definition: compareCandidates.c:45
solver_params::fvalue
double fvalue
Definition: compareCandidates.c:46
solver_params::data_array
double * data_array
Definition: compareCandidates.c:47