LALPulsar 7.1.1.1-eeff03c
DriveHoughMulti.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2005 Badri Krishnan, Alicia Sintes
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 USA2
18 */
19
20/**
21 * \file
22 * \ingroup lalpulsar_bin_Hough
23 * \author Badri Krishnan, Alicia Sintes, Llucia Sancho
24 * \brief Driver code for performing Hough transform search on non-demodulated
25 * data using SFTs from possible multiple IFOs
26 *
27 * History: Created by Sintes and Krishnan July 04, 2003
28 * Modifications for S4 January 2006
29 * Modifications for S5 November 2007
30 *
31 * \par Description
32 *
33 * This is the main driver for the Hough transform routines. It takes as input
34 * a set of SFTs from possibly more than one IFO and outputs the number counts
35 * using the Hough transform. For a single IFO, this should be essentially equivalent
36 * to DriveHough_v3. validatehoughmultichi2
37 *
38 * This code just does spin-downs values.
39 *
40 * \par User input
41 *
42 * The user inputs the following parameters:
43 *
44 * - Search frequency range
45 *
46 * - A file containing list of skypatches to search over. For each skypatch,
47 * the information is:
48 * - RA and dec of skypatch center.
49 * - Size in RA and dec.
50 *
51 * - Location of Directory containing the SFTs (must be v2 SFTs).
52 *
53 * - Interferometer (optional)
54 *
55 * - List of linefiles containing information about known spectral disturbanclves
56 *
57 * - Location of output directory and basename of output files.
58 *
59 * - Block size of running median for estimating psd.
60 *
61 * - The parameter nfSizeCylinder which determines the range of spindown parameters
62 * to be searched over.
63 *
64 * - Boolean variable for deciding if the SFTs should be inverse noise weighed.
65 *
66 * - Boolean variable for deciding whether amplitude modulation weights should be used.
67 *
68 * - Boolean variables for deciding whether the Hough maps, the statistics, list of
69 * events above a threshold, and logfile should be written
70 *
71 * /par Output
72 *
73 * The output is written in several sub-directories of the specified output directory. The
74 * first two items are default while the rest are written according to the user input:
75 *
76 * - A directory called logfiles records the user input, contents of the skypatch file
77 * and cvs tags contained in the executable (if the user has required logging)
78 *
79 * - A directory called nstarfiles containing the loudest event for each search frequency
80 * bin maximised over all sky-locations and spindown parameters. An event is said to be
81 * the loudest if it has the maximum significance defined as: (number count - mean)/sigma.
82 *
83 * - A directory for each skypatch containing the number count statistics, the histogram,
84 * the list of events, and the Hough maps
85 */
86
87#include "config.h"
88
89#include <gsl/gsl_permutation.h>
90
91#include <lal/LogPrintf.h>
92#include <lal/DopplerScan.h>
93#include <lal/LALPulsarVCSInfo.h>
94
95#include "DriveHoughColor.h"
96#include "MCInjectHoughMulti.h"
97#include "FstatToplist.h"
98
99/* globals, constants and defaults */
100
101/* boolean global variables for controlling output */
102BOOLEAN uvar_EnableExtraInfo, uvar_EnableChi2, uvar_EnableToplistPatch;
103
104/* #define EARTHEPHEMERIS "./earth05-09.dat" */
105/* #define SUNEPHEMERIS "./sun05-09.dat" */
106
107#define EARTHEPHEMERIS "/home/badkri/lscsoft/share/lal/earth05-09.dat"
108#define SUNEPHEMERIS "/home/badkri/lscsoft/share/lal/sun05-09.dat"
109
110#define HOUGHMAXFILENAMELENGTH 512 /* maximum # of characters of a filename */
111
112#define DIROUT "./outMulti" /* output directory */
113#define BASENAMEOUT "HM" /* prefix file output */
114
115#define THRESHOLD 1.6 /* thresold for peak selection, with respect to the
116the averaged power in the search band */
117#define SKYFILE "./skypatchfile"
118#define F0 310.0 /* frequency to build the LUT and start search */
119#define FBAND 0.05 /* search frequency band */
120#define NFSIZE 21 /* n-freq. span of the cylinder, to account for spin-down search */
121#define NSPINUP 0 /* n-bins for spin-up search */
122#define BLOCKSRNGMED 101 /* Running median window size */
123
124#define SKYREGION "allsky" /**< default sky region to search over -- just a single point*/
125
126#define TRUE (1==1)
127#define FALSE (1==0)
128
129#define NBLOCKSTEST 8 /* number of data blocks to do Chi2 test */
130
131#define SPINDOWNJUMP 1 /* "Jump" to the next spin-down being analyzed (to avoid doing them all) */
132
133/* ****************************************
134 * Structure, HoughParamsTest, typedef
135 */
136
137typedef struct tagHoughParamsTest {
138 UINT4 length; /* number p of blocks to split the data into */
139 UINT4 *numberSFTp; /* Ni SFTs in each data block */
140 REAL8 *sumWeight; /* Pointer to the sumWeight of each block of data */
141 REAL8 *sumWeightSquare; /* Pointer to the sumWeightSquare of each block of data */
142} HoughParamsTest;
143
144
145typedef struct tagUCHARPeakGramVector {
146 UINT4 length; /**< number of elements */
147 UCHARPeakGram *upg; /**< expanded Peakgrams */
148} UCHARPeakGramVector;
149
150/******************************************/
151
152/* local function prototype */
153
154void SplitSFTs( LALStatus *status, REAL8Vector *weightsV, HoughParamsTest *chi2Params );
155
156void ComputeFoft_NM( LALStatus *status, REAL8Vector *foft, HoughTemplate *pulsarTemplate, REAL8Vector *timeDiffV, REAL8Cart3CoorVector *velV );
157
158void ComputeandPrintChi2( LALStatus *status, toplist_t *tl, REAL8Vector *timeDiffV, REAL8Cart3CoorVector *velV, INT4 skyCounter, INT4 nSkyPatches, INT4 p, REAL8 alphaPeak, MultiDetectorStateSeries *mdetStates, REAL8Vector *weightsNoise, UCHARPeakGramVector *upgV );
159
160void GetPeakGramFromMultSFTVector_NondestroyPg1( LALStatus *status, HOUGHPeakGramVector *out, UCHARPeakGramVector *upgV, MultiSFTVector *in, REAL8 thr ) ;
161
162void PrintLogFile( LALStatus *status, CHAR *dir, CHAR *basename, CHAR *skyfile, LALStringVector *linefiles, CHAR *executable );
163
164int CreateSkypatchDirs( CHAR *filestats, CHAR *base, INT4 ind );
165
166int PrintHistogram( UINT8Vector *hist, CHAR *fnameOut, REAL8 minSignificance, REAL8 maxSignificance );
167
168int PrintHmap2m_file( HOUGHMapTotal *ht, CHAR *fnameOut, INT4 iHmap );
169
170int PrintHmap2file( HOUGHMapTotal *ht, CHAR *fnameOut, INT4 iHmap );
171
172int OpenExtraInfoFiles( CHAR *fileMaps, FILE **fp1_ptr, CHAR *filehisto, CHAR *dirname, CHAR *basename, INT4 ind );
173
174int PrintExtraInfo( CHAR *fileMaps, FILE **fp1_ptr, INT4 iHmap, HOUGHMapTotal *ht, REAL8UnitPolarCoor *sourceLocation, HoughStats *stats, INT8 fBinSearch, REAL8 deltaF );
175
176void ReadTimeStampsFile( LALStatus *status, LIGOTimeGPSVector *ts, CHAR *filename );
177
178void LALHoughHistogramSignificance( LALStatus *status, UINT8Vector *out, HOUGHMapTotal *in,
179 REAL8 mean, REAL8 sigma, REAL8 minSignificance, REAL8 maxSignificance );
180
181void GetSFTVelTime( LALStatus *status, REAL8Cart3CoorVector *velV, LIGOTimeGPSVector *timeV, MultiDetectorStateSeries *in );
182
183void GetSFTNoiseWeights( LALStatus *status, REAL8Vector *out, MultiNoiseWeights *in );
184
185void GetPeakGramFromMultSFTVector( LALStatus *status, HOUGHPeakGramVector *out, MultiSFTVector *in, REAL8 thr );
186
187void SetUpSkyPatches( LALStatus *status, HoughSkyPatchesInfo *out, CHAR *skyFileName, CHAR *skyRegion, REAL8 dAlpha, REAL8 dDelta, INT4 numSkyPartitions, INT4 partitionIndex );
188
189void GetAMWeights( LALStatus *status, REAL8Vector *out, MultiDetectorStateSeries *mdetStates, REAL8 alpha, REAL8 delta );
190
191void SelectBestStuff( LALStatus *status, BestVariables *out, BestVariables *in, UINT4 mObsCohBest );
192
193void DuplicateBestStuff( LALStatus *status, BestVariables *out, BestVariables *in );
194
195void GetToplistFromHoughmap( LALStatus *status, toplist_t *list, HOUGHMapTotal *ht, HOUGHPatchGrid *patch, HOUGHDemodPar *parDem, REAL8 mean, REAL8 sigma );
196
197
198void LALHOUGHCreateLUTVector( LALStatus *status,
199 HOUGHptfLUTVector *lutV,
200 UINT4 length );
201
202void LALHOUGHCreateLUTs( LALStatus *status,
203 HOUGHptfLUTVector *lutV,
204 UINT2 maxNBins,
205 UINT2 maxNBorders,
206 UINT2 ySide );
207
208void LALHOUGHCreatePHMDVS( LALStatus *status,
209 PHMDVectorSequence *phmdVS,
210 UINT4 length,
211 UINT4 nfSize );
212
213void LALHOUGHCreatePHMDs( LALStatus *status,
214 PHMDVectorSequence *phmdVS,
215 UINT2 maxNBins,
216 UINT2 maxNBorders,
217 UINT2 ySide );
218
219void LALHOUGHDestroyPHMDs( LALStatus *status,
220 PHMDVectorSequence *phmdVS );
221
222void LALHOUGHCreateHT( LALStatus *status,
223 HOUGHMapTotal *ht,
224 UINT2 xSide,
225 UINT2 ySide );
226
227void LALHOUGHCreateFreqIndVector( LALStatus *status,
228 UINT8FrequencyIndexVector *freqInd,
229 UINT4 length,
230 REAL8 deltaF );
231
232
233void LALHOUGHDestroyLUTs( LALStatus *status,
234 HOUGHptfLUTVector *lut );
235
236/******************************************/
237
238int main( int argc, char *argv[] )
239{
240
241 /* LALStatus pointer */
242 static LALStatus status;
243
244 /* time and velocity */
245 LIGOTimeGPSVector *timeV = NULL;
246 REAL8Vector *timeDiffV = NULL;
247 static REAL8Cart3CoorVector velV;
248
249
250 LIGOTimeGPS firstTimeStamp, lastTimeStamp;
251 REAL8 tObs;
252
253 /* standard pulsar sft types */
254 MultiSFTVector *inputSFTs = NULL;
255 UINT4 numSearchBins;
256
257 /* information about all the ifos */
258 MultiDetectorStateSeries *mdetStates = NULL;
259 UINT4 numifo;
260
261 /* vector of weights */
262 REAL8Vector *weightsV = NULL, *weightsNoise = NULL;
263
264 /* ephemeris */
265 EphemerisData *edat = NULL;
266
267 /* struct containing subset of weights, timediff, velocity and peakgrams */
268 static BestVariables best;
269
270 /* hough structures */
271 static HOUGHptfLUTVector lutV; /* the Look Up Table vector*/
272 static HOUGHPeakGramVector pgV; /* vector of peakgrams */
273 static UCHARPeakGramVector upgV; /* vector of expanded peakgrams */
274 static PHMDVectorSequence phmdVS; /* the partial Hough map derivatives */
275 static UINT8FrequencyIndexVector freqInd; /* for trajectory in time-freq plane */
276 static HOUGHResolutionPar parRes; /* patch grid information */
277 static HOUGHPatchGrid patch; /* Patch description */
278 static HOUGHParamPLUT parLut; /* parameters needed to build lut */
279 static HOUGHDemodPar parDem; /* demodulation parameters or */
280 static HOUGHSizePar parSize;
281 static HOUGHMapTotal ht; /* the total Hough map */
282 static UINT8Vector *hist; /* histogram of number counts for a single map */
283 static UINT8Vector *histTotal; /* number count histogram for all maps */
284 static HoughStats stats; /* statistical information about a Hough map */
285
286 /* skypatch info */
287 REAL8 *skyAlpha = NULL, *skyDelta = NULL, *skySizeAlpha = NULL, *skySizeDelta = NULL;
288 INT4 nSkyPatches, skyCounter = 0;
289 static HoughSkyPatchesInfo skyInfo;
290
291 /* output filenames and filepointers */
292 CHAR filehisto[ HOUGHMAXFILENAMELENGTH ];
293 /* CHAR filestats[ HOUGHMAXFILENAMELENGTH ];*/
294 CHAR fileMaps[ HOUGHMAXFILENAMELENGTH ];
295 CHAR fileSigma[ HOUGHMAXFILENAMELENGTH ];
296 FILE *fp1 = NULL;
297 FILE *fpSigma = NULL;
298
299 /* the maximum number count */
300 static HoughSignificantEventVector nStarEventVec;
301
302 /* miscellaneous */
303 INT4 iHmap, nSpin1Max ;
304 UINT4 mObsCoh, mObsCohBest;
305 INT8 f0Bin, fLastBin, fBin;
306 REAL8 alpha, delta, timeBase, deltaF, f1jump;
307 REAL8 patchSizeX, patchSizeY;
308 REAL8 alphaPeak, minSignificance, maxSignificance;
309 REAL8 meanN, sigmaN;
310 UINT2 xSide, ySide;
311 UINT2 maxNBins, maxNBorders;
312
313 /* output toplist candidate structure */
314 toplist_t *toplist = NULL;
315
316 /* sft constraint variables */
317 LIGOTimeGPS startTimeGPS, endTimeGPS;
318 LIGOTimeGPSVector inputTimeStampsVector;
319
320 /* user input variables */
321 BOOLEAN uvar_weighAM, uvar_weighNoise, uvar_printLog;
322 INT4 uvar_blocksRngMed, uvar_nfSizeCylinder, uvar_nSpinUp, uvar_maxBinsClean, uvar_binsHisto;
323 INT4 uvar_keepBestSFTs = 1;
324 INT4 uvar_numCand;
325 REAL8 uvar_startTime, uvar_endTime;
326 REAL8 uvar_pixelFactor;
327 REAL8 uvar_f0, uvar_peakThreshold, uvar_freqBand;
328 REAL8 uvar_dAlpha, uvar_dDelta; /* resolution for isotropic sky-grid */
329 CHAR *uvar_earthEphemeris = NULL;
330 CHAR *uvar_sunEphemeris = NULL;
331 CHAR *uvar_sftData = NULL;
332 CHAR *uvar_dirnameOut = NULL;
333 CHAR *uvar_fbasenameOut = NULL;
334 CHAR *uvar_skyfile = NULL;
335 CHAR *uvar_timeStampsFile = NULL;
336 CHAR *uvar_skyRegion = NULL;
337 LALStringVector *uvar_linefiles = NULL;
338
339 INT4 uvar_chiSqBins;
340
341 INT4 uvar_spindownJump;
342
343 INT4 uvar_nfLUTvalidity = 0;
344 INT4 uvar_numSkyPartitions = 0;
345 INT4 uvar_partitionIndex = 0;
346
347 LIGOTimeGPS refTimeGPS; /* reference time */
348 REAL8 uvar_refTime;
349
350 REAL8 uvar_deltaF1dot;
351
352 /* Set up the default parameters */
353
354 /* LAL error-handler */
355 lal_errhandler = LAL_ERR_EXIT;
356
357 uvar_weighAM = TRUE;
358 uvar_weighNoise = TRUE;
359 uvar_printLog = FALSE;
360 uvar_blocksRngMed = BLOCKSRNGMED;
361 uvar_nfSizeCylinder = NFSIZE;
362 uvar_nSpinUp = NSPINUP;
363 uvar_f0 = F0;
364 uvar_freqBand = FBAND;
365 uvar_peakThreshold = THRESHOLD;
366 uvar_maxBinsClean = 100;
367 uvar_binsHisto = 1000;
368 uvar_pixelFactor = PIXELFACTOR;
369 uvar_dAlpha = 0.2;
370 uvar_dDelta = 0.2;
371 uvar_numCand = 1;
372 uvar_EnableExtraInfo = FALSE;
373 uvar_EnableChi2 = FALSE;
374 uvar_chiSqBins = NBLOCKSTEST;
375 uvar_spindownJump = SPINDOWNJUMP;
376
377 uvar_EnableToplistPatch = FALSE;
378
379
380
381 uvar_earthEphemeris = ( CHAR * )LALCalloc( HOUGHMAXFILENAMELENGTH, sizeof( CHAR ) );
382 strcpy( uvar_earthEphemeris, EARTHEPHEMERIS );
383
384 uvar_sunEphemeris = ( CHAR * )LALCalloc( HOUGHMAXFILENAMELENGTH, sizeof( CHAR ) );
385 strcpy( uvar_sunEphemeris, SUNEPHEMERIS );
386
387 uvar_dirnameOut = ( CHAR * )LALCalloc( HOUGHMAXFILENAMELENGTH, sizeof( CHAR ) );
388 strcpy( uvar_dirnameOut, DIROUT );
389
390 uvar_fbasenameOut = ( CHAR * )LALCalloc( HOUGHMAXFILENAMELENGTH, sizeof( CHAR ) );
391 strcpy( uvar_fbasenameOut, BASENAMEOUT );
392
393 uvar_skyfile = ( CHAR * )LALCalloc( HOUGHMAXFILENAMELENGTH, sizeof( CHAR ) );
394 strcpy( uvar_skyfile, SKYFILE );
395
396 /* register user input variables */
397 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_f0, "f0", REAL8, 'f', OPTIONAL, "Start search frequency" ) == XLAL_SUCCESS, XLAL_EFUNC );
398 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_freqBand, "freqBand", REAL8, 'b', OPTIONAL, "Search frequency band" ) == XLAL_SUCCESS, XLAL_EFUNC );
399 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_startTime, "startTime", REAL8, 0, OPTIONAL, "GPS start time of observation (SFT timestamps must be >= this)" ) == XLAL_SUCCESS, XLAL_EFUNC );
400 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_endTime, "endTime", REAL8, 0, OPTIONAL, "GPS end time of observation (SFT timestamps must be < this)" ) == XLAL_SUCCESS, XLAL_EFUNC );
401 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_timeStampsFile, "timeStampsFile", STRING, 0, OPTIONAL, "Input time-stamps file" ) == XLAL_SUCCESS, XLAL_EFUNC );
402 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_skyRegion, "skyRegion", STRING, 0, OPTIONAL, "sky-region polygon (or 'allsky')" ) == XLAL_SUCCESS, XLAL_EFUNC );
403 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dAlpha, "dAlpha", REAL8, 0, OPTIONAL, "Resolution for flat or isotropic coarse grid (rad)" ) == XLAL_SUCCESS, XLAL_EFUNC );
404 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dDelta, "dDelta", REAL8, 0, OPTIONAL, "Resolution for flat or isotropic coarse grid (rad)" ) == XLAL_SUCCESS, XLAL_EFUNC );
405 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_skyfile, "skyfile", STRING, 0, OPTIONAL, "Alternative: input skypatch file" ) == XLAL_SUCCESS, XLAL_EFUNC );
406 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_peakThreshold, "peakThreshold", REAL8, 0, OPTIONAL, "Peak selection threshold" ) == XLAL_SUCCESS, XLAL_EFUNC );
407 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_weighAM, "weighAM", BOOLEAN, 0, OPTIONAL, "Use amplitude modulation weights" ) == XLAL_SUCCESS, XLAL_EFUNC );
408 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_weighNoise, "weighNoise", BOOLEAN, 0, OPTIONAL, "Use SFT noise weights" ) == XLAL_SUCCESS, XLAL_EFUNC );
409 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_keepBestSFTs, "keepBestSFTs", INT4, 0, OPTIONAL, "Number of best SFTs to use (default--keep all)" ) == XLAL_SUCCESS, XLAL_EFUNC );
410 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_printLog, "printLog", BOOLEAN, 0, OPTIONAL, "Print Log file" ) == XLAL_SUCCESS, XLAL_EFUNC );
411 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_earthEphemeris, "earthEphemeris", STRING, 'E', OPTIONAL, "Earth Ephemeris file" ) == XLAL_SUCCESS, XLAL_EFUNC );
412 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_sunEphemeris, "sunEphemeris", STRING, 'S', OPTIONAL, "Sun Ephemeris file" ) == XLAL_SUCCESS, XLAL_EFUNC );
413 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_sftData, "sftData", STRING, 'D', REQUIRED, "SFT filename pattern. Possibilities are:\n"
414 " - '<SFT file>;<SFT file>;...', where <SFT file> may contain wildcards\n - 'list:<file containing list of SFT files>'" ) == XLAL_SUCCESS, XLAL_EFUNC );
415 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_dirnameOut, "dirnameOut", STRING, 'o', OPTIONAL, "Output directory" ) == XLAL_SUCCESS, XLAL_EFUNC );
416 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_fbasenameOut, "fbasenameOut", STRING, 0, OPTIONAL, "Output file basename" ) == XLAL_SUCCESS, XLAL_EFUNC );
417 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_binsHisto, "binsHisto", INT4, 0, OPTIONAL, "No. of bins for histogram" ) == XLAL_SUCCESS, XLAL_EFUNC );
418 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_linefiles, "linefiles", STRINGVector, 0, OPTIONAL, "Comma separated List of linefiles (filenames must contain IFO name)" ) == XLAL_SUCCESS, XLAL_EFUNC );
419 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nfSizeCylinder, "nfSizeCylinder", INT4, 0, OPTIONAL, "Size of cylinder of PHMDs" ) == XLAL_SUCCESS, XLAL_EFUNC );
420 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nfLUTvalidity, "nfLUTvalidity", INT4, 0, OPTIONAL, "Frequency bins validity of LUT (default most restrictive value)" ) == XLAL_SUCCESS, XLAL_EFUNC );
421 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_nSpinUp, "nSpinUp", INT4, 0, OPTIONAL, "Num of bins for Spin-up in PHMDs" ) == XLAL_SUCCESS, XLAL_EFUNC );
422 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_pixelFactor, "pixelFactor", REAL8, 'p', OPTIONAL, "sky resolution=1/v*pixelFactor*f*Tcoh" ) == XLAL_SUCCESS, XLAL_EFUNC );
423 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_numCand, "numCand", INT4, 0, OPTIONAL, "No. of toplist candidates" ) == XLAL_SUCCESS, XLAL_EFUNC );
424 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_EnableExtraInfo, "printExtraInfo", BOOLEAN, 0, OPTIONAL, "Print HoughMaps, HoughStatistics, expected number count stdev" ) == XLAL_SUCCESS, XLAL_EFUNC );
425 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_chiSqBins, "chiSqBins", INT4, 0, OPTIONAL, "Number of chi-square bins for veto tests" ) == XLAL_SUCCESS, XLAL_EFUNC );
426 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_EnableChi2, "enableChi2", BOOLEAN, 0, OPTIONAL, "Print Chi2 value for each element in the Toplist" ) == XLAL_SUCCESS, XLAL_EFUNC );
427 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_spindownJump, "spindownJump", INT4, 0, OPTIONAL, "Jump to the next spin-down being analyzed (to avoid doing them all)" ) == XLAL_SUCCESS, XLAL_EFUNC );
428 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_numSkyPartitions, "numSkyPartitions", INT4, 0, OPTIONAL, "Number of (equi-)partitions to split skygrid into" ) == XLAL_SUCCESS, XLAL_EFUNC );
429 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_partitionIndex, "partitionIndex", INT4, 0, OPTIONAL, "Index [0,xnumSkyPartitions-1] of sky-partition to generate" ) == XLAL_SUCCESS, XLAL_EFUNC );
430 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_refTime, "refTime", REAL8, 0, OPTIONAL, "GPS reference time of observation" ) == XLAL_SUCCESS, XLAL_EFUNC );
431 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_deltaF1dot, "deltaF1dot", REAL8, 0, OPTIONAL, "(Step size for f1dot)*Tcoh [Default: 1/Tobs]" ) == XLAL_SUCCESS, XLAL_EFUNC );
432
433 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_EnableToplistPatch, "EnableToplistPatch", BOOLEAN, 0, OPTIONAL, "Enables a toplist per Patch, requires to enableChi2" ) == XLAL_SUCCESS, XLAL_EFUNC );
434
435
436 /* developer input variables */
437 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_blocksRngMed, "blocksRngMed", INT4, 0, DEVELOPER, "Running Median block size" ) == XLAL_SUCCESS, XLAL_EFUNC );
438 XLAL_CHECK_MAIN( XLALRegisterNamedUvar( &uvar_maxBinsClean, "maxBinsClean", INT4, 0, DEVELOPER, "Maximum number of bins in cleaning" ) == XLAL_SUCCESS, XLAL_EFUNC );
439
440
441 /* read all command line variables */
442 BOOLEAN should_exit = 0;
443 XLAL_CHECK_MAIN( XLALUserVarReadAllInput( &should_exit, argc, argv, lalPulsarVCSInfoList ) == XLAL_SUCCESS, XLAL_EFUNC );
444 if ( should_exit ) {
445 exit( 1 );
446 }
447
448 /* very basic consistency checks on user input */
449 if ( uvar_f0 < 0 ) {
450 LogPrintf( LOG_CRITICAL, "start frequency must be positive\n" );
451 exit( 1 );
452 }
453
454 if ( uvar_freqBand < 0 ) {
455 LogPrintf( LOG_CRITICAL, "search frequency band must be positive\n" );
456 exit( 1 );
457 }
458
459 if ( uvar_peakThreshold < 0 ) {
460 LogPrintf( LOG_CRITICAL, "peak selection threshold must be positive\n" );
461 exit( 1 );
462 }
463
464 /* probability of peak selection */
465 alphaPeak = exp( -uvar_peakThreshold );
466
467
468 if ( uvar_binsHisto < 1 ) {
469 LogPrintf( LOG_CRITICAL, "binsHisto must be at least 1\n" );
470 exit( 1 );
471 }
472
473 if ( uvar_nfSizeCylinder < uvar_nSpinUp + 1 ) {
474 LogPrintf( LOG_CRITICAL, "Number of search spin must bigger than the positive ones\n" );
475 exit( 1 );
476 }
477
478 if ( uvar_keepBestSFTs < 1 ) {
479 LogPrintf( LOG_CRITICAL, "must keep at least 1 SFT\n" );
480 exit( 1 );
481 }
482
483 /* write log file with command line arguments, cvs tags, and contents of skypatch file */
484 if ( uvar_printLog ) {
485 LAL_CALL( PrintLogFile( &status, uvar_dirnameOut, uvar_fbasenameOut, uvar_skyfile, uvar_linefiles, argv[0] ), &status );
486 }
487
488
489 /***** start main calculations *****/
490
491 LogPrintf( LOG_NORMAL, "Setting up sky-patches..." );
492 /* set up skypatches */
493 LAL_CALL( SetUpSkyPatches( &status, &skyInfo, uvar_skyfile, uvar_skyRegion, uvar_dAlpha, uvar_dDelta, uvar_numSkyPartitions, uvar_partitionIndex ), &status );
494 nSkyPatches = skyInfo.numSkyPatches;
495 skyAlpha = skyInfo.alpha;
496 skyDelta = skyInfo.delta;
497 skySizeAlpha = skyInfo.alphaSize;
498 skySizeDelta = skyInfo.deltaSize;
499 LogPrintfVerbatim( LOG_NORMAL, "done\n" );
500
501 /* set up toplist */
502 /* create toplist -- semiCohToplist has the same structure
503 as a fstat candidate, so treat it as a fstat candidate */
504 if ( create_fstat_toplist( &toplist, uvar_numCand ) != 0 ) {
505 LogPrintf( LOG_CRITICAL, "Unable to create toplist\n" );
506 }
507
508
509 LogPrintf( LOG_NORMAL, "Reading SFTs..." );
510 /* read sft Files and set up weights */
511 {
512 /* new SFT I/O data types */
513 SFTCatalog *catalog = NULL;
514 static SFTConstraints constraints;
515
516 REAL8 doppWings, f_min, f_max;
517 INT4 k;
518
519 /* set detector constraint */
520 constraints.detector = NULL;
521
522 if ( XLALUserVarWasSet( &uvar_startTime ) ) {
523 XLALGPSSetREAL8( &startTimeGPS, uvar_startTime );
524 constraints.minStartTime = &startTimeGPS;
525 }
526
527 if ( XLALUserVarWasSet( &uvar_endTime ) ) {
528 XLALGPSSetREAL8( &endTimeGPS, uvar_endTime );
529 constraints.maxStartTime = &endTimeGPS;
530 }
531
532 if ( XLALUserVarWasSet( &uvar_timeStampsFile ) ) {
533 LAL_CALL( ReadTimeStampsFile( &status, &inputTimeStampsVector, uvar_timeStampsFile ), &status );
534 constraints.timestamps = &inputTimeStampsVector;
535 }
536
537 /* get sft catalog */
538 XLAL_CHECK_MAIN( ( catalog = XLALSFTdataFind( uvar_sftData, &constraints ) ) != NULL, XLAL_EFUNC );
539 if ( ( catalog == NULL ) || ( catalog->length == 0 ) ) {
540 LogPrintf( LOG_CRITICAL, "Unable to match any SFTs with pattern '%s'\n", uvar_sftData );
541 exit( 1 );
542 }
543
544 /* now we can free the inputTimeStampsVector */
545 if ( XLALUserVarWasSet( &uvar_timeStampsFile ) ) {
546 LALFree( inputTimeStampsVector.data );
547 }
548
549
550
551 /* first some sft parameters */
552 deltaF = catalog->data[0].header.deltaF; /* frequency resolution */
553 timeBase = 1.0 / deltaF; /* coherent integration time */
554 f0Bin = floor( uvar_f0 * timeBase + 0.5 ); /* initial search frequency */
555 numSearchBins = uvar_freqBand * timeBase; /* total number of search bins - 1 */
556 fLastBin = f0Bin + numSearchBins; /* final frequency bin to be analyzed */
557
558
559 /* read sft files making sure to add extra bins for running median */
560 /* add wings for Doppler modulation and running median block size*/
561 doppWings = ( uvar_f0 + uvar_freqBand ) * VTOT;
562 f_min = uvar_f0 - doppWings - ( uvar_blocksRngMed + uvar_nfSizeCylinder ) * deltaF;
563 f_max = uvar_f0 + uvar_freqBand + doppWings + ( uvar_blocksRngMed + uvar_nfSizeCylinder ) * deltaF;
564
565 /* read the sfts */
566 XLAL_CHECK_MAIN( ( inputSFTs = XLALLoadMultiSFTs( catalog, f_min, f_max ) ) != NULL, XLAL_EFUNC );
567 numifo = inputSFTs->length;
568
569 /* find number of sfts */
570 /* loop over ifos and calculate number of sfts */
571 /* note that we can't use the catalog to determine the number of SFTs
572 because SFTs might be segmented in frequency */
573 mObsCoh = 0; /* initialization */
574 for ( k = 0; k < ( INT4 )numifo; k++ ) {
575 mObsCoh += inputSFTs->data[k]->length;
576 }
577
578 /* set number of SFTs to be kept */
579 /* currently mobscohbest is set equal to mobscoh if no weights are used
580 -- this probably will be changed in the future */
581 if ( XLALUserVarWasSet( &uvar_keepBestSFTs ) && ( uvar_weighNoise || uvar_weighAM ) ) {
582 mObsCohBest = uvar_keepBestSFTs;
583
584 /* set to mobscoh if it is more than number of sfts */
585 if ( mObsCohBest > mObsCoh ) {
586 mObsCohBest = mObsCoh;
587 }
588 } else {
589 mObsCohBest = mObsCoh;
590 }
591
592
593 /* catalog is ordered in time so we can get start, end time and tObs*/
594 firstTimeStamp = catalog->data[0].header.epoch;
595 lastTimeStamp = catalog->data[catalog->length - 1].header.epoch;
596
597 if ( XLALUserVarWasSet( &uvar_refTime ) ) {
598 XLALGPSSetREAL8( &refTimeGPS, uvar_refTime );
599 tObs = XLALGPSDiff( &lastTimeStamp, &refTimeGPS ) + timeBase;
600 } else {
601 tObs = XLALGPSDiff( &lastTimeStamp, &firstTimeStamp ) + timeBase;
602 }
603
604
605
606
607 /* clean sfts if required */
608 if ( XLALUserVarWasSet( &uvar_linefiles ) ) {
609
610 RandomParams *randPar = NULL;
611 FILE *fpRand = NULL;
612 INT4 seed, ranCount;
613
614 LogPrintfVerbatim( LOG_NORMAL, "...cleaning SFTs..." );
615 if ( ( fpRand = fopen( "/dev/urandom", "r" ) ) == NULL ) {
616 LogPrintf( LOG_CRITICAL, "error in opening /dev/urandom" );
617 exit( 1 );
618 }
619
620 if ( ( ranCount = fread( &seed, sizeof( seed ), 1, fpRand ) ) != 1 ) {
621 LogPrintf( LOG_CRITICAL, "error in getting random seed" );
622 exit( 1 );
623 }
624
625 LAL_CALL( LALCreateRandomParams( &status, &randPar, seed ), &status );
626
627 LAL_CALL( LALRemoveKnownLinesInMultiSFTVector( &status, inputSFTs, uvar_maxBinsClean,
628 uvar_blocksRngMed, uvar_linefiles, randPar ), &status );
629
630 LAL_CALL( LALDestroyRandomParams( &status, &randPar ), &status );
631 fclose( fpRand );
632 } /* end cleaning */
633
634
635 XLALDestroySFTCatalog( catalog );
636
637 } /* end of sft reading block */
638 LogPrintfVerbatim( LOG_NORMAL, "done\n" );
639
640
641 /** some memory allocations */
642
643 /* allocate memory for velocity vector */
644 velV.length = mObsCoh;
645 velV.data = NULL;
646 velV.data = ( REAL8Cart3Coor * )LALCalloc( 1, mObsCoh * sizeof( REAL8Cart3Coor ) );
647
648 /* allocate memory for timestamps and timediff vectors */
649 timeV = XLALCreateTimestampVector( mObsCoh );
650 timeDiffV = XLALCreateREAL8Vector( mObsCoh );
651
652 /* allocate and initialize noise and AMweights vectors */
653 weightsV = XLALCreateREAL8Vector( mObsCoh );
654 weightsNoise = XLALCreateREAL8Vector( mObsCoh );
655 LAL_CALL( LALHOUGHInitializeWeights( &status, weightsNoise ), &status );
656 LAL_CALL( LALHOUGHInitializeWeights( &status, weightsV ), &status );
657
658
659 LogPrintf( LOG_NORMAL, "Setting up weights..." );
660 /* get detector velocities weights vector, and timestamps */
661 {
662 MultiNoiseWeights *multweight = NULL;
663 MultiPSDVector *multPSD = NULL;
664 UINT4 j;
665
666 /* get ephemeris */
667 XLAL_CHECK_MAIN( ( edat = XLALInitBarycenter( uvar_earthEphemeris, uvar_sunEphemeris ) ) != NULL, XLAL_EFUNC );
668
669
670 /* normalize sfts */
671 XLAL_CHECK_MAIN( ( multPSD = XLALNormalizeMultiSFTVect( inputSFTs, uvar_blocksRngMed, NULL ) ) != NULL, XLAL_EFUNC );
672
673 /* compute multi noise weights */
674 if ( uvar_weighNoise ) {
675 XLAL_CHECK_MAIN( ( multweight = XLALComputeMultiNoiseWeights( multPSD, uvar_blocksRngMed, 0 ) ) != NULL, XLAL_EFUNC );
676 }
677
678 /* we are now done with the psd */
679 XLALDestroyMultiPSDVector( multPSD );
680
681 /* get information about all detectors including velocity and timestamps */
682 /* note that this function returns the velocity at the
683 mid-time of the SFTs -- should not make any difference */
684 const REAL8 tOffset = 0.5 / inputSFTs->data[0]->data[0].deltaF;
685 XLAL_CHECK_MAIN( ( mdetStates = XLALGetMultiDetectorStatesFromMultiSFTs( inputSFTs, edat, tOffset ) ) != NULL, XLAL_EFUNC );
686
687 LAL_CALL( GetSFTVelTime( &status, &velV, timeV, mdetStates ), &status );
688
689 /* copy the noise-weights vector if required*/
690 if ( uvar_weighNoise ) {
691
692 LAL_CALL( GetSFTNoiseWeights( &status, weightsNoise, multweight ), &status );
693
694 XLALDestroyMultiNoiseWeights( multweight );
695 }
696
697 /* compute the time difference relative to startTime for all SFTs */
698 if ( XLALUserVarWasSet( &uvar_refTime ) ) {
699 for ( j = 0; j < mObsCoh; j++ ) {
700 timeDiffV->data[j] = XLALGPSDiff( timeV->data + j, &refTimeGPS );
701 }
702 } else {
703 for ( j = 0; j < mObsCoh; j++ ) {
704 timeDiffV->data[j] = XLALGPSDiff( timeV->data + j, &firstTimeStamp );
705 }
706 }
707
708 } /* end block for weights, velocity and time */
709 LogPrintfVerbatim( LOG_NORMAL, "done\n" );
710
711
712 LogPrintf( LOG_NORMAL, "Generating peakgrams..." );
713 /* generating peakgrams */
714 pgV.length = mObsCoh;
715 pgV.pg = NULL;
716 pgV.pg = ( HOUGHPeakGram * )LALCalloc( 1, mObsCoh * sizeof( HOUGHPeakGram ) );
717
718 if ( uvar_EnableChi2 ) {
719 upgV.length = mObsCoh;
720 upgV.upg = NULL;
721 upgV.upg = ( UCHARPeakGram * )LALCalloc( 1, mObsCoh * sizeof( UCHARPeakGram ) );
722
723 LAL_CALL( GetPeakGramFromMultSFTVector_NondestroyPg1( &status, &pgV, &upgV, inputSFTs, uvar_peakThreshold ), &status );
724 } else {
725 LAL_CALL( GetPeakGramFromMultSFTVector( &status, &pgV, inputSFTs, uvar_peakThreshold ), &status );
726 }
727
728
729 /* we are done with the sfts and ucharpeakgram now */
730 XLALDestroyMultiSFTVector( inputSFTs );
731 LogPrintfVerbatim( LOG_NORMAL, "done\n" );
732
733 /* if we want to print expected sigma for each skypatch */
734 if ( uvar_EnableExtraInfo ) {
735 strcpy( fileSigma, uvar_dirnameOut );
736 strcat( fileSigma, "/" );
737 strcat( fileSigma, uvar_fbasenameOut );
738 strcat( fileSigma, "sigma" );
739
740 if ( ( fpSigma = fopen( fileSigma, "w" ) ) == NULL ) {
741 LogPrintf( LOG_CRITICAL, "Unable to find file %s for writing\n", fileSigma );
742 return DRIVEHOUGHCOLOR_EFILE;
743 }
744 } /* end if( uvar_EnableExtraInfo) */
745
746
747 /* min and max values significance that are possible */
748 minSignificance = -sqrt( mObsCohBest * alphaPeak / ( 1 - alphaPeak ) );
749 maxSignificance = sqrt( mObsCohBest * ( 1 - alphaPeak ) / alphaPeak );
750
751
752
753
754 if ( uvar_EnableToplistPatch ) {
755 LogPrintf( LOG_NORMAL, "Starting loop over skypatches and chi-square follow-up of top candidates per patch: " );
756 } else {
757 LogPrintf( LOG_NORMAL, "Starting loop over skypatches..." );
758 }
759
760 /* loop over sky patches -- main Hough calculations */
761 for ( skyCounter = 0; skyCounter < nSkyPatches; skyCounter++ ) {
762 UINT4 k;
763 REAL8 sumWeightSquare;
764 /* REAL8 meanN, sigmaN;*/
765 BestVariables temp;
766
767 LogPrintfVerbatim( LOG_NORMAL, "%d/%d,", skyCounter, nSkyPatches );
768
769 /* set sky positions and skypatch sizes */
770 alpha = skyAlpha[skyCounter];
771 delta = skyDelta[skyCounter];
772 patchSizeX = skySizeDelta[skyCounter];
773 patchSizeY = skySizeAlpha[skyCounter];
774
775 /* copy noise weights if required */
776 if ( uvar_weighNoise ) {
777 memcpy( weightsV->data, weightsNoise->data, mObsCoh * sizeof( REAL8 ) );
778 }
779
780 /* calculate amplitude modulation weights if required */
781 if ( uvar_weighAM ) {
782 LAL_CALL( GetAMWeights( &status, weightsV, mdetStates, alpha, delta ), &status );
783 }
784
785 /* sort weights vector to get the best sfts */
786 temp.length = mObsCoh;
787 temp.weightsV = weightsV;
788 temp.timeDiffV = timeDiffV;
789 temp.velV = &velV;
790 temp.pgV = &pgV;
791
792 if ( uvar_weighAM || uvar_weighNoise ) {
793 LAL_CALL( SelectBestStuff( &status, &best, &temp, mObsCohBest ), &status );
794 } else {
795 LAL_CALL( DuplicateBestStuff( &status, &best, &temp ), &status );
796 }
797
798 /* Normalize the Best SFTs weights */
799 LAL_CALL( LALHOUGHNormalizeWeights( &status, best.weightsV ), &status );
800
801 /* calculate the sum of the weights squared */
802 sumWeightSquare = 0.0;
803 for ( k = 0; k < mObsCohBest; k++ ) {
804 sumWeightSquare += best.weightsV->data[k] * best.weightsV->data[k];
805 }
806
807 /* probability of selecting a peak expected mean and standard deviation for noise only */
808 meanN = mObsCohBest * alphaPeak;
809 sigmaN = sqrt( sumWeightSquare * alphaPeak * ( 1.0 - alphaPeak ) );
810
811
812 if ( uvar_EnableExtraInfo ) {
813 fprintf( fpSigma, "%f\n", sigmaN );
814 if ( OpenExtraInfoFiles( fileMaps, &fp1, filehisto, uvar_dirnameOut, uvar_fbasenameOut, skyCounter ) ) {
815 return DRIVEHOUGHCOLOR_EFILE;
816 }
817 }
818
819 /**** general parameter settings and 1st memory allocation ****/
820
821 LAL_CALL( LALHOUGHCreateLUTVector( &status, &lutV, mObsCohBest ), &status );
822
823 LAL_CALL( LALHOUGHCreatePHMDVS( &status, &phmdVS, mObsCohBest, uvar_nfSizeCylinder ), &status );
824 phmdVS.deltaF = deltaF;
825
826 LAL_CALL( LALHOUGHCreateFreqIndVector( &status, &freqInd, mObsCohBest, deltaF ), &status );
827
828 /* allocating histogram of the number-counts in the Hough maps */
829 if ( uvar_EnableExtraInfo ) {
830 UINT4 k0;
831 hist = XLALCreateUINT8Vector( uvar_binsHisto );
832 histTotal = XLALCreateUINT8Vector( uvar_binsHisto );
833
834 /* initialize to 0 */
835 /* memset(histTotal->data, histTotal->length*sizeof(histTotal->data[0]), 0); */
836 for ( k0 = 0; k0 < histTotal->length; k0++ ) {
837 histTotal->data[k0] = 0;
838 hist->data[k0] = 0;
839 }
840 }
841
842 /* set demodulation pars for non-demodulated data (SFT input)*/
843 parDem.deltaF = deltaF;
844 parDem.skyPatch.alpha = alpha;
845 parDem.skyPatch.delta = delta;
846 parDem.timeDiff = 0.0;
847 parDem.spin.length = 0;
848 parDem.spin.data = NULL;
849 parDem.positC.x = 0.0;
850 parDem.positC.y = 0.0;
851 parDem.positC.z = 0.0;
852
853 /* sky-resolution parameters **/
854 parRes.deltaF = deltaF;
855 parRes.patchSkySizeX = patchSizeX;
856 parRes.patchSkySizeY = patchSizeY;
857 parRes.pixelFactor = uvar_pixelFactor;
858 parRes.pixErr = PIXERR;
859 parRes.linErr = LINERR;
860 parRes.vTotC = VTOT;
861
862
863
864 fBin = f0Bin;
865 iHmap = 0;
866
867 /* ***** for spin-down case ****/
868 nSpin1Max = uvar_nfSizeCylinder - 1 - uvar_nSpinUp;
869 /* nSpin1Max = floor(uvar_nfSizeCylinder/2.0) ;*/
870
871
872 if ( XLALUserVarWasSet( &uvar_deltaF1dot ) ) {
873 f1jump = uvar_deltaF1dot * uvar_spindownJump;
874 } else {
875 f1jump = 1. / tObs * uvar_spindownJump;
876 }
877
878
879 /* start of main loop over search frequency bins */
880 /********** starting the search from f0Bin to fLastBin.
881 Note one set LUT might not cover all the interval.
882 This is taken into account *******************/
883
884 while ( fBin <= fLastBin ) {
885 INT8 fBinSearch, fBinSearchMax;
886 UINT4 j;
887 REAL8UnitPolarCoor sourceLocation;
888
889
890 parRes.f0Bin = fBin;
891 LAL_CALL( LALHOUGHComputeNDSizePar( &status, &parSize, &parRes ), &status );
892 xSide = parSize.xSide;
893 ySide = parSize.ySide;
894 maxNBins = parSize.maxNBins;
895 maxNBorders = parSize.maxNBorders;
896 if ( uvar_nfLUTvalidity ) {
897 parSize.nFreqValid = uvar_nfLUTvalidity;
898 }
899
900 /* *******************create patch grid at fBin **************** */
901 patch.xSide = xSide;
902 patch.ySide = ySide;
903 patch.xCoor = NULL;
904 patch.yCoor = NULL;
905 patch.xCoor = ( REAL8 * )LALCalloc( 1, xSide * sizeof( REAL8 ) );
906 patch.yCoor = ( REAL8 * )LALCalloc( 1, ySide * sizeof( REAL8 ) );
907 LAL_CALL( LALHOUGHFillPatchGrid( &status, &patch, &parSize ), &status );
908
909 /*************** other memory allocation and settings************ */
910
911 LAL_CALL( LALHOUGHCreateLUTs( &status, &lutV, maxNBins, maxNBorders, ySide ), &status );
912
913 LAL_CALL( LALHOUGHCreatePHMDs( &status, &phmdVS, maxNBins, maxNBorders, ySide ), &status );
914
915
916 /* ************* create all the LUTs at fBin ******************** */
917 for ( j = 0; j < mObsCohBest; ++j ) { /* create all the LUTs */
918 parDem.veloC.x = best.velV->data[j].x;
919 parDem.veloC.y = best.velV->data[j].y;
920 parDem.veloC.z = best.velV->data[j].z;
921 /* calculate parameters needed for buiding the LUT */
922 LAL_CALL( LALNDHOUGHParamPLUT( &status, &parLut, &parSize, &parDem ), &status );
923 /* build the LUT */
924 LAL_CALL( LALHOUGHConstructPLUT( &status, &( lutV.lut[j] ), &patch, &parLut ),
925 &status );
926 }
927
928 /************* build the set of PHMD centered around fBin***********/
929 phmdVS.fBinMin = fBin - uvar_nfSizeCylinder + 1 + uvar_nSpinUp;
930 /*phmdVS.fBinMin = fBin - floor( uvar_nfSizeCylinder/2.) ;*/
931
932 LAL_CALL( LALHOUGHConstructSpacePHMD( &status, &phmdVS, best.pgV, &lutV ), &status );
933 if ( uvar_weighAM || uvar_weighNoise ) {
934 LAL_CALL( LALHOUGHWeighSpacePHMD( &status, &phmdVS, best.weightsV ), &status );
935 }
936
937 /* ************ initializing the Total Hough map space *********** */
938
939 LAL_CALL( LALHOUGHCreateHT( &status, &ht, xSide, ySide ), &status );
940 ht.mObsCoh = mObsCohBest;
941 ht.deltaF = deltaF;
942
943
944 /* Search frequency interval possible using the same LUTs */
945 fBinSearch = fBin;
946 fBinSearchMax = fBin + parSize.nFreqValid - 1 - uvar_nSpinUp;
947 /*fBinSearchMax = fBin + parSize.nFreqValid - 1 - floor((uvar_nfSizeCylinder/2. ) ;*/
948
949
950 /* Study all possible frequencies with one set of LUT */
951
952 while ( ( fBinSearch <= fLastBin ) && ( fBinSearch < fBinSearchMax ) ) {
953
954 /**** study 1 spin-down. at fBinSearch ****/
955
956 INT4 n;
957 REAL8 f1dis;
958
959 ht.f0Bin = fBinSearch;
960 ht.spinRes.length = 1;
961 ht.spinRes.data = NULL;
962 ht.spinRes.data = ( REAL8 * )LALCalloc( ht.spinRes.length, sizeof( REAL8 ) );
963 for ( n = floor( uvar_nSpinUp / uvar_spindownJump ); n >= - floor( nSpin1Max / uvar_spindownJump ); --n ) {
964 /*for ( n = 0; n <= floor(nSpin1Max/uvar_spindownJump); ++n) {*/
965 /* f1dis = - n * f1jump; */
966 /*loop over all spindown values */
967
968 f1dis = + n * f1jump;
969 ht.spinRes.data[0] = f1dis * deltaF;
970
971 /* construct path in time-freq plane */
972 for ( j = 0 ; j < mObsCohBest; ++j ) {
973 freqInd.data[j] = fBinSearch + floor( best.timeDiffV->data[j] * f1dis + 0.5 );
974 }
975
976 if ( uvar_weighAM || uvar_weighNoise ) {
977 LAL_CALL( LALHOUGHConstructHMT_W( &status, &ht, &freqInd, &phmdVS ), &status );
978 } else {
979 LAL_CALL( LALHOUGHConstructHMT( &status, &ht, &freqInd, &phmdVS ), &status );
980 }
981
982
983 /* ********************* perfom stat. analysis on the maps ****************** */
984
985 if ( uvar_EnableExtraInfo ) {
986
987 LAL_CALL( LALHoughStatistics( &status, &stats, &ht ), &status );
988 LAL_CALL( LALStereo2SkyLocation( &status, &sourceLocation,
989 stats.maxIndex[0], stats.maxIndex[1], &patch, &parDem ), &status );
990
991 /*LAL_CALL( LALHoughHistogram ( &status, &hist, &ht), &status);*/
992 LAL_CALL( LALHoughHistogramSignificance( &status, hist, &ht, meanN, sigmaN,
993 minSignificance, maxSignificance ), &status );
994
995 for ( j = 0; j < histTotal->length; j++ ) {
996 histTotal->data[j] += hist->data[j];
997 }
998 }
999
1000 /* select candidates from hough maps */
1001 LAL_CALL( GetToplistFromHoughmap( &status, toplist, &ht, &patch, &parDem, meanN, sigmaN ), &status );
1002
1003
1004 /* ***** print results *********************** */
1005
1006 if ( uvar_EnableExtraInfo ) {
1007 if ( PrintExtraInfo( fileMaps, &fp1, iHmap, &ht, &sourceLocation, &stats, fBinSearch, deltaF ) ) {
1008 return DRIVEHOUGHCOLOR_EFILE;
1009 }
1010 }
1011
1012 ++iHmap;
1013 } /* end loop over spindown values */
1014
1015 LALFree( ht.spinRes.data );
1016
1017
1018 /***** shift the search freq. & PHMD structure 1 freq.bin ****** */
1019 ++fBinSearch;
1020
1021 LAL_CALL( LALHOUGHupdateSpacePHMDup( &status, &phmdVS, best.pgV, &lutV ), &status );
1022
1023 if ( uvar_weighAM || uvar_weighNoise ) {
1024 LAL_CALL( LALHOUGHWeighSpacePHMD( &status, &phmdVS, best.weightsV ), &status );
1025 }
1026
1027 } /*closing second while */
1028
1029 fBin = fBinSearch;
1030
1031 /* ******************** Free partial memory ******************* */
1032 LALFree( patch.xCoor );
1033 LALFree( patch.yCoor );
1034 LALFree( ht.map );
1035
1036 LALHOUGHDestroyLUTs( &status, &lutV );
1037
1038 LALHOUGHDestroyPHMDs( &status, &phmdVS );
1039
1040
1041 } /* closing while */
1042
1043 /* printing toplist per patch and free toplist memory */
1044 if ( uvar_EnableToplistPatch ) {
1045 if ( uvar_EnableChi2 ) {
1046 LAL_CALL( ComputeandPrintChi2( &status, toplist, timeDiffV, &velV, skyCounter, nSkyPatches, uvar_chiSqBins, alphaPeak, mdetStates, weightsNoise, &upgV ), &status );
1047 free_fstat_toplist( &toplist );
1048 if ( create_fstat_toplist( &toplist, uvar_numCand ) != 0 ) {
1049 LogPrintf( LOG_CRITICAL, "Unable to create toplist\n" );
1050 }
1051 }
1052 }
1053
1054 /* printing total histogram */
1055 if ( uvar_EnableExtraInfo ) {
1056 if ( PrintHistogram( histTotal, filehisto, minSignificance, maxSignificance ) ) {
1057 return 7;
1058 }
1059 }
1060
1061 /* --------------------------------------------------*/
1062 /* Closing files with statistics results and events*/
1063 if ( uvar_EnableExtraInfo ) {
1064 fclose( fp1 );
1065 }
1066
1067 /* Free memory allocated inside skypatches loop */
1068 LALFree( lutV.lut );
1069 lutV.lut = NULL;
1070
1071 LALFree( phmdVS.phmd );
1072 phmdVS.phmd = NULL;
1073
1074 LALFree( freqInd.data );
1075 freqInd.data = NULL;
1076
1077 if ( uvar_EnableExtraInfo ) {
1078 XLALDestroyUINT8Vector( hist );
1079 XLALDestroyUINT8Vector( histTotal );
1080 }
1081
1082 } /* finish loop over skypatches */
1083 LogPrintfVerbatim( LOG_NORMAL, "...done\n" );
1084
1085
1086 /* close sigma file */
1087 if ( uvar_EnableExtraInfo ) {
1088 fclose( fpSigma );
1089 }
1090
1091 /*********************************************************/
1092 /* print toplist */
1093 /********************************************************/
1094
1095 /* If we want to print Chi2 value */
1096
1097 if ( uvar_EnableChi2 ) {
1098 if ( !uvar_EnableToplistPatch ) {
1099 LogPrintf( LOG_NORMAL, "Starting chi-square follow-up of top candidates..." );
1100 LAL_CALL( ComputeandPrintChi2( &status, toplist, timeDiffV, &velV, -1, 0, uvar_chiSqBins, alphaPeak, mdetStates, weightsNoise, &upgV ), &status );
1101 LogPrintfVerbatim( LOG_NORMAL, "done\n" );
1102 }
1103 } else {
1104
1105 FILE *fpToplist = NULL;
1106
1107 LogPrintf( LOG_NORMAL, "Sort and print toplist..." );
1108 fpToplist = fopen( "hough_top.dat", "w" );
1109
1110 sort_fstat_toplist( toplist );
1111
1112 if ( write_fstat_toplist_to_fp( toplist, fpToplist, NULL ) < 0 ) {
1113 LogPrintf( LOG_CRITICAL, "error in writing toplist to file...\n" );
1114 }
1115
1116 if ( fprintf( fpToplist, "%%DONE\n" ) < 0 ) {
1117 LogPrintf( LOG_CRITICAL, "error writing end marker...\n" );
1118 }
1119
1120 fclose( fpToplist );
1121 LogPrintfVerbatim( LOG_NORMAL, "done\n" );
1122 }
1123 /********************************************************/
1124
1125 LogPrintf( LOG_NORMAL, "Free memory and exit..." );
1126 {
1127 UINT4 j;
1128 for ( j = 0; j < mObsCoh; ++j ) {
1129 LALFree( pgV.pg[j].peak );
1130 }
1131 }
1132
1133 LALFree( pgV.pg );
1134
1135 if ( uvar_EnableChi2 ) {
1136 {
1137 UINT4 j;
1138 for ( j = 0; j < mObsCoh; ++j ) {
1139 LALFree( upgV.upg[j].data );
1140 }
1141 }
1142 LALFree( upgV.upg );
1143
1144 }
1145
1146 XLALDestroyTimestampVector( timeV );
1147 XLALDestroyREAL8Vector( timeDiffV );
1148
1149
1150 LALFree( velV.data );
1151
1152 XLALDestroyREAL8Vector( weightsV );
1153 XLALDestroyREAL8Vector( weightsNoise );
1154
1155 XLALDestroyMultiDetectorStateSeries( mdetStates );
1156
1157 XLALDestroyEphemerisData( edat );
1158
1159 LALFree( skyAlpha );
1160 LALFree( skyDelta );
1161 LALFree( skySizeAlpha );
1162 LALFree( skySizeDelta );
1163
1164 if ( nStarEventVec.event ) {
1165 LALFree( nStarEventVec.event );
1166 }
1167
1168 LALFree( best.weightsV->data );
1169 LALFree( best.weightsV );
1170 LALFree( best.timeDiffV->data );
1171 LALFree( best.timeDiffV );
1172 LALFree( best.velV->data );
1173 LALFree( best.velV );
1174 LALFree( best.pgV->pg );
1175 LALFree( best.pgV );
1176
1177 free_fstat_toplist( &toplist );
1178
1179 XLALDestroyUserVars();
1180
1181 LALCheckMemoryLeaks();
1182
1183 LogPrintfVerbatim( LOG_NORMAL, "bye\n" );
1184
1185 /* if ( lalDebugLevel ) */
1186 /* REPORTSTATUS ( &status); */
1187
1188 return status.statusCode;
1189}
1190
1191
1192
1193
1194/******************************************************************/
1195/* printing the Histogram of all maps into a file */
1196/******************************************************************/
1197
1198int CreateSkypatchDirs( CHAR *filestats,
1199 CHAR *base,
1200 INT4 ind )
1201{
1202 CHAR tempstr[16];
1203 INT4 mkdir_result;
1204
1205 /* create the directory name uvar_dirnameOut/skypatch_$j */
1206 strcpy( filestats, base );
1207
1208 strcat( filestats, "/skypatch_" );
1209 sprintf( tempstr, "%d", ind );
1210 strcat( filestats, tempstr );
1211 strcat( filestats, "/" );
1212
1213 /* check whether file can be created or if it exists already
1214 if not then exit */
1215 errno = 0;
1216 mkdir_result = mkdir( filestats, S_IRWXU | S_IRWXG | S_IRWXO );
1217
1218 if ( ( mkdir_result == -1 ) && ( errno != EEXIST ) ) {
1219 fprintf( stderr, "unable to create skypatch directory %d\n", ind );
1220 return 1;
1221 }
1222
1223 return 0;
1224
1225}
1226
1227
1228int PrintHistogram( UINT8Vector *hist, CHAR *fnameOut, REAL8 minSignificance, REAL8 maxSignificance )
1229{
1230
1231 INT4 binsHisto;
1232 REAL8 dSig;
1233 FILE *fp = NULL; /* Output file */
1234 char filename[ HOUGHMAXFILENAMELENGTH ];
1235 INT4 i ;
1236
1237 strcpy( filename, fnameOut );
1238 strcat( filename, "histo" );
1239
1240 binsHisto = hist->length;
1241 dSig = ( maxSignificance - minSignificance ) / binsHisto;
1242
1243 if ( ( fp = fopen( filename, "w" ) ) == NULL ) {
1244 fprintf( stderr, "Unable to find file %s\n", filename );
1245 return DRIVEHOUGHCOLOR_EFILE;
1246 }
1247
1248 for ( i = 0; i < binsHisto; i++ ) {
1249 fprintf( fp, "%g %" LAL_UINT8_FORMAT "\n", minSignificance + i * dSig, hist->data[i] );
1250 }
1251
1252 fclose( fp );
1253 return 0;
1254
1255}
1256
1257
1258
1259/******************************************************************/
1260/* printing the HM into a file */
1261/******************************************************************/
1262
1263int PrintHmap2file( HOUGHMapTotal *ht, CHAR *fnameOut, INT4 iHmap )
1264{
1265
1266 FILE *fp = NULL; /* Output file */
1267 char filename[ HOUGHMAXFILENAMELENGTH ], filenumber[16];
1268 INT4 k, i ;
1269 UINT2 xSide, ySide;
1270
1271 strcpy( filename, fnameOut );
1272 sprintf( filenumber, ".%06d", iHmap );
1273 strcat( filename, filenumber );
1274
1275 if ( ( fp = fopen( filename, "w" ) ) == NULL ) {
1276 fprintf( stderr, "Unable to find file %s\n", filename );
1277 return DRIVEHOUGHCOLOR_EFILE;
1278 }
1279
1280 ySide = ht->ySide;
1281 xSide = ht->xSide;
1282
1283 for ( k = ySide - 1; k >= 0; --k ) {
1284 for ( i = 0; i < xSide; ++i ) {
1285 fprintf( fp, " %f", ( REAL4 )ht->map[k * xSide + i] );
1286 fflush( fp );
1287 }
1288 fprintf( fp, " \n" );
1289 fflush( fp );
1290 }
1291
1292 fclose( fp );
1293 return 0;
1294}
1295
1296
1297/******************************************************************/
1298/* printing the HM into a m_file */
1299/******************************************************************/
1300
1301int PrintHmap2m_file( HOUGHMapTotal *ht, CHAR *fnameOut, INT4 iHmap )
1302{
1303
1304 FILE *fp = NULL; /* Output file */
1305 char filename[ HOUGHMAXFILENAMELENGTH ], filenumber[16];
1306 INT4 k, i ;
1307 UINT2 xSide, ySide;
1308 REAL8 f0, f1;
1309
1310 strcpy( filename, fnameOut );
1311 sprintf( filenumber, "%06d.m", iHmap );
1312 strcat( filename, filenumber );
1313 fp = fopen( filename, "w" );
1314
1315 if ( !fp ) {
1316 fprintf( stderr, "Unable to find file %s\n", filename );
1317 return DRIVEHOUGHCOLOR_EFILE;
1318 }
1319
1320 ySide = ht->ySide;
1321 xSide = ht->xSide;
1322 f0 = ht->f0Bin * ht->deltaF;
1323 f1 = 0.0;
1324 if ( ht->spinRes.length ) {
1325 f1 = ht->spinRes.data[0];
1326 }
1327
1328 /* printing into matlab format */
1329
1330 fprintf( fp, "f0= %f ; \n", f0 );
1331 fprintf( fp, "f1= %g ; \n", f1 );
1332 fprintf( fp, "map= [ \n" );
1333
1334 for ( k = ySide - 1; k >= 0; --k ) {
1335 for ( i = 0; i < xSide; ++i ) {
1336 fprintf( fp, " %f", ( REAL4 )ht->map[k * xSide + i] );
1337 fflush( fp );
1338 }
1339 fprintf( fp, " \n" );
1340 fflush( fp );
1341 }
1342
1343 fprintf( fp, " ]; \n" );
1344 fclose( fp );
1345 return 0;
1346}
1347
1348/******************************************************************************************/
1349
1350void PrintLogFile( LALStatus *status,
1351 CHAR *dir,
1352 CHAR *basename,
1353 CHAR *skyfile,
1354 LALStringVector *linefiles,
1355 CHAR *executable )
1356{
1357 CHAR *fnameLog = NULL;
1358 FILE *fpLog = NULL;
1359 CHAR *logstr = NULL;
1360 UINT4 k;
1361
1362 INITSTATUS( status );
1363 ATTATCHSTATUSPTR( status );
1364
1365 /* open log file for writing */
1366 fnameLog = ( CHAR * )LALCalloc( HOUGHMAXFILENAMELENGTH, sizeof( CHAR ) );
1367 strcpy( fnameLog, dir );
1368 strcat( fnameLog, "/logfiles/" );
1369 /* now create directory fdirOut/logfiles using mkdir */
1370 errno = 0;
1371 {
1372 /* check whether file can be created or if it exists already
1373 if not then exit */
1374 INT4 mkdir_result;
1375 mkdir_result = mkdir( fnameLog, S_IRWXU | S_IRWXG | S_IRWXO );
1376 if ( ( mkdir_result == -1 ) && ( errno != EEXIST ) ) {
1377 fprintf( stderr, "unable to create logfiles directory %s\n", fnameLog );
1378 LALFree( fnameLog );
1379 exit( 1 ); /* stop the program */
1380 }
1381 }
1382
1383 /* create the logfilename in the logdirectory */
1384 strcat( fnameLog, basename );
1385 strcat( fnameLog, ".log" );
1386 /* open the log file for writing */
1387 if ( ( fpLog = fopen( fnameLog, "w" ) ) == NULL ) {
1388 fprintf( stderr, "Unable to open file %s for writing\n", fnameLog );
1389 LALFree( fnameLog );
1390 exit( 1 );
1391 }
1392
1393 /* get the log string */
1394 XLAL_CHECK_LAL( status, ( logstr = XLALUserVarGetLog( UVAR_LOGFMT_CFGFILE ) ) != NULL, XLAL_EFUNC );
1395
1396 fprintf( fpLog, "## LOG FILE FOR Hough Driver\n\n" );
1397 fprintf( fpLog, "# User Input:\n" );
1398 fprintf( fpLog, "#-------------------------------------------\n" );
1399 fprintf( fpLog, "%s", logstr );
1400 LALFree( logstr );
1401
1402 /* copy contents of skypatch file into logfile */
1403 fprintf( fpLog, "\n\n# Contents of skypatch file:\n" );
1404 fclose( fpLog );
1405 {
1406 CHAR command[1024] = "";
1407 sprintf( command, "cat %s >> %s", skyfile, fnameLog );
1408 if ( system( command ) ) {
1409 fprintf( stderr, "\nsystem('%s') returned non-zero status!\n\n", command );
1410 }
1411 }
1412
1413
1414 /* copy contents of linefile if necessary */
1415 if ( linefiles ) {
1416
1417 for ( k = 0; k < linefiles->length; k++ ) {
1418
1419 if ( ( fpLog = fopen( fnameLog, "a" ) ) != NULL ) {
1420 CHAR command[1024] = "";
1421 fprintf( fpLog, "\n\n# Contents of linefile %s :\n", linefiles->data[k] );
1422 fprintf( fpLog, "# -----------------------------------------\n" );
1423 fclose( fpLog );
1424 sprintf( command, "cat %s >> %s", linefiles->data[k], fnameLog );
1425 if ( system( command ) ) {
1426 fprintf( stderr, "\nsystem('%s') returned non-zero status!\n\n", command );
1427 }
1428 }
1429 }
1430 }
1431
1432 /* append an ident-string defining the exact CVS-version of the code used */
1433 if ( ( fpLog = fopen( fnameLog, "a" ) ) != NULL ) {
1434 CHAR command[1024] = "";
1435 fprintf( fpLog, "\n\n# CVS-versions of executable:\n" );
1436 fprintf( fpLog, "# -----------------------------------------\n" );
1437 fclose( fpLog );
1438
1439 sprintf( command, "ident %s | sort -u >> %s", executable, fnameLog );
1440 /* we don't check this. If it fails, we assume that */
1441 /* one of the system-commands was not available, and */
1442 /* therefore the CVS-versions will not be logged */
1443 if ( system( command ) ) {
1444 fprintf( stderr, "\nsystem('%s') returned non-zero status!\n\n", command );
1445 }
1446 }
1447
1448 LALFree( fnameLog );
1449
1450 DETATCHSTATUSPTR( status );
1451 /* normal exit */
1452 RETURN( status );
1453}
1454
1455
1456/* read timestamps file */
1457void ReadTimeStampsFile( LALStatus *status,
1458 LIGOTimeGPSVector *ts,
1459 CHAR *filename )
1460{
1461
1462 FILE *fp = NULL;
1463 INT4 numTimeStamps, r;
1464 UINT4 j;
1465 REAL8 temp1, temp2;
1466
1467 INITSTATUS( status );
1468 ATTATCHSTATUSPTR( status );
1469
1470 ASSERT( ts, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1471 ASSERT( ts->data == NULL, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1472 ASSERT( ts->length == 0, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1473 ASSERT( filename, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1474
1475 if ( ( fp = fopen( filename, "r" ) ) == NULL ) {
1476 ABORT( status, DRIVEHOUGHCOLOR_EFILE, DRIVEHOUGHCOLOR_MSGEFILE );
1477 }
1478
1479 /* count number of timestamps */
1480 numTimeStamps = 0;
1481
1482 do {
1483 r = fscanf( fp, "%lf%lf\n", &temp1, &temp2 );
1484 /* make sure the line has the right number of entries or is EOF */
1485 if ( r == 2 ) {
1486 numTimeStamps++;
1487 }
1488 } while ( r != EOF );
1489 rewind( fp );
1490
1491 ts->length = numTimeStamps;
1492 ts->data = LALCalloc( 1, numTimeStamps * sizeof( LIGOTimeGPS ) );;
1493 if ( ts->data == NULL ) {
1494 fclose( fp );
1495 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1496 }
1497
1498 for ( j = 0; j < ts->length; j++ ) {
1499 r = fscanf( fp, "%lf%lf\n", &temp1, &temp2 );
1500 ts->data[j].gpsSeconds = ( INT4 )temp1;
1501 ts->data[j].gpsNanoSeconds = ( INT4 )temp2;
1502 }
1503
1504 fclose( fp );
1505
1506 DETATCHSTATUSPTR( status );
1507 /* normal exit */
1508 RETURN( status );
1509}
1510
1511
1512
1513/**
1514 * given a total hough map, this function produces a histogram of
1515 * the number count significance
1516 */
1517void LALHoughHistogramSignificance( LALStatus *status,
1518 UINT8Vector *out,
1519 HOUGHMapTotal *in,
1520 REAL8 mean,
1521 REAL8 sigma,
1522 REAL8 minSignificance,
1523 REAL8 maxSignificance )
1524{
1525
1526 INT4 i, j, binsHisto, xSide, ySide, binIndex;
1527 REAL8 temp;
1528
1529 INITSTATUS( status );
1530 ATTATCHSTATUSPTR( status );
1531
1532 /* make sure arguments are not null */
1533 ASSERT( in, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1534 ASSERT( in->map, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1535 ASSERT( out, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1536
1537 /* make sure input hough map is ok*/
1538 ASSERT( in->xSide > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1539 ASSERT( in->ySide > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1540 ASSERT( in->mObsCoh > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1541
1542 ASSERT( out->length > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1543
1544 binsHisto = out->length;
1545 xSide = in->xSide;
1546 ySide = in->ySide;
1547
1548 /* initialize histogram vector*/
1549 for ( i = 0; i < binsHisto; i++ ) {
1550 out->data[i] = 0;
1551 }
1552
1553 /* loop over hough map and find histogram */
1554 for ( i = 0; i < ySide; i++ ) {
1555 for ( j = 0; j < xSide; j++ ) {
1556
1557 /* calculate significance of number count */
1558 temp = ( in->map[i * xSide + j] - mean ) / sigma;
1559
1560 /* make sure temp is in proper range */
1561 ASSERT( temp > minSignificance, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1562 ASSERT( temp < maxSignificance, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1563
1564 binIndex = ( INT4 ) binsHisto * ( temp - minSignificance ) / ( maxSignificance - minSignificance );
1565
1566 /* make sure binIndex is in proper range */
1567 ASSERT( binIndex >= 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1568 ASSERT( binIndex <= binsHisto - 1, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1569
1570 /* add to relevant entry in histogram */
1571 out->data[binIndex] += 1;
1572 }
1573 }
1574
1575 DETATCHSTATUSPTR( status );
1576
1577 /* normal exit */
1578 RETURN( status );
1579}
1580
1581
1582
1583void GetSFTVelTime( LALStatus *status,
1584 REAL8Cart3CoorVector *velV,
1585 LIGOTimeGPSVector *timeV,
1586 MultiDetectorStateSeries *in )
1587{
1588
1589 UINT4 numifo, numsft, iIFO, iSFT, j;
1590
1591 INITSTATUS( status );
1592 ATTATCHSTATUSPTR( status );
1593
1594 ASSERT( in, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1595 ASSERT( in->length > 0, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1596
1597 ASSERT( velV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1598 ASSERT( velV->length > 0, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1599 ASSERT( velV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1600
1601 ASSERT( timeV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1602 ASSERT( timeV->length > 0, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1603 ASSERT( timeV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1604
1605 ASSERT( velV->length == timeV->length, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1606
1607 numifo = in->length;
1608
1609 /* copy the timestamps, weights, and velocity vector */
1610 for ( j = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
1611
1612 ASSERT( in->data[iIFO], status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1613
1614 numsft = in->data[iIFO]->length;
1615 ASSERT( numsft > 0, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1616
1617 for ( iSFT = 0; iSFT < numsft; iSFT++, j++ ) {
1618
1619 velV->data[j].x = in->data[iIFO]->data[iSFT].vDetector[0];
1620 velV->data[j].y = in->data[iIFO]->data[iSFT].vDetector[1];
1621 velV->data[j].z = in->data[iIFO]->data[iSFT].vDetector[2];
1622
1623 /* mid time of sfts */
1624 timeV->data[j] = in->data[iIFO]->data[iSFT].tGPS;
1625
1626 } /* loop over SFTs */
1627
1628 } /* loop over IFOs */
1629
1630 DETATCHSTATUSPTR( status );
1631
1632 /* normal exit */
1633 RETURN( status );
1634}
1635
1636
1637
1638
1639void GetSFTNoiseWeights( LALStatus *status,
1640 REAL8Vector *out,
1641 MultiNoiseWeights *in )
1642{
1643
1644 UINT4 numifo, numsft, iIFO, iSFT, j;
1645
1646 INITSTATUS( status );
1647 ATTATCHSTATUSPTR( status );
1648
1649 ASSERT( in, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1650 ASSERT( in->length > 0, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1651
1652 ASSERT( out, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1653 ASSERT( out->length > 0, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1654 ASSERT( out->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1655
1656 numifo = in->length;
1657
1658 /* copy the timestamps, weights, and velocity vector */
1659 for ( j = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
1660
1661 ASSERT( in->data[iIFO], status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1662
1663 numsft = in->data[iIFO]->length;
1664 ASSERT( numsft > 0, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1665
1666 for ( iSFT = 0; iSFT < numsft; iSFT++, j++ ) {
1667
1668 out->data[j] = in->data[iIFO]->data[iSFT];
1669
1670 } /* loop over SFTs */
1671
1672 } /* loop over IFOs */
1673
1674 TRY( LALHOUGHNormalizeWeights( status->statusPtr, out ), status );
1675
1676 DETATCHSTATUSPTR( status );
1677
1678 /* normal exit */
1679 RETURN( status );
1680}
1681
1682
1683/** Loop over SFTs and set a threshold to get peakgrams. SFTs must be normalized. */
1684void GetPeakGramFromMultSFTVector( LALStatus *status, /**< pointer to LALStatus structure */
1685 HOUGHPeakGramVector *out, /**< Output peakgrams */
1686 MultiSFTVector *in, /**< Input SFTs */
1687 REAL8 thr /**< Threshold on SFT power */ )
1688{
1689
1690 SFTtype *sft;
1691 UCHARPeakGram pg1;
1692 INT4 nPeaks;
1693 UINT4 iIFO, iSFT, numsft, numifo, j, binsSFT;
1694
1695 INITSTATUS( status );
1696 ATTATCHSTATUSPTR( status );
1697
1698 numifo = in->length;
1699 binsSFT = in->data[0]->data->data->length;
1700
1701 pg1.length = binsSFT;
1702 pg1.data = NULL;
1703 pg1.data = ( UCHAR * )LALCalloc( 1, binsSFT * sizeof( UCHAR ) );
1704
1705 /* loop over sfts and select peaks */
1706 for ( j = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
1707
1708 numsft = in->data[iIFO]->length;
1709
1710 for ( iSFT = 0; iSFT < numsft; iSFT++, j++ ) {
1711
1712 sft = in->data[iIFO]->data + iSFT;
1713
1714 TRY( SFTtoUCHARPeakGram( status->statusPtr, &pg1, sft, thr ), status );
1715
1716 nPeaks = pg1.nPeaks;
1717
1718 /* compress peakgram */
1719 out->pg[j].length = nPeaks;
1720 out->pg[j].peak = NULL;
1721 out->pg[j].peak = ( INT4 * )LALCalloc( 1, nPeaks * sizeof( INT4 ) );
1722
1723 TRY( LALUCHAR2HOUGHPeak( status->statusPtr, &( out->pg[j] ), &pg1 ), status );
1724
1725 } /* loop over SFTs */
1726
1727 } /* loop over IFOs */
1728
1729 LALFree( pg1.data );
1730
1731 DETATCHSTATUSPTR( status );
1732
1733 /* normal exit */
1734 RETURN( status );
1735
1736}
1737
1738
1739/**
1740 * Set up location of skypatch centers and sizes
1741 * If user specified skyRegion then use DopplerScan function
1742 * to construct an isotropic grid. Otherwise use skypatch file.
1743 */
1744void SetUpSkyPatches( LALStatus *status, /**< pointer to LALStatus structure */
1745 HoughSkyPatchesInfo *out, /**< output skypatches info */
1746 CHAR *skyFileName, /**< name of skypatch file */
1747 CHAR *skyRegion, /**< skyregion (if isotropic grid is to be constructed) */
1748 REAL8 dAlpha, /**< alpha resolution (if isotropic grid is to be constructed) */
1749 REAL8 dDelta, /**< delta resolution (if isotropic grid is to be constructed) */
1750 INT4 numSkyPartitions,/**<Number of (equi-)partitions to split skygrid into */
1751 INT4 partitionIndex ) /**< Index [0,numSkyPartitions-1] of sky-partition to generate */
1752{
1753
1754 DopplerSkyScanInit XLAL_INIT_DECL( scanInit ); /* init-structure for DopperScanner */
1755 DopplerSkyScanState XLAL_INIT_DECL( thisScan ); /* current state of the Doppler-scan */
1756 UINT4 nSkyPatches, skyCounter;
1757 PulsarDopplerParams dopplerpos;
1758
1759 INITSTATUS( status );
1760 ATTATCHSTATUSPTR( status );
1761
1762 ASSERT( out, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1763 ASSERT( dAlpha > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1764 ASSERT( dDelta > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1765
1766 if ( skyRegion ) {
1767
1768 scanInit.dAlpha = dAlpha;
1769 scanInit.dDelta = dDelta;
1770 scanInit.gridType = GRID_ISOTROPIC;
1771 scanInit.metricType = LAL_PMETRIC_NONE;
1772 /* scanInit.metricMismatch = 0; */
1773 /* scanInit.projectMetric = TRUE; */
1774 /* scanInit.obsDuration = tStack; */
1775 /* scanInit.obsBegin = tMidGPS; */
1776 /* scanInit.Detector = &(stackMultiDetStates.data[0]->data[0]->detector); */
1777 /* scanInit.ephemeris = edat; */
1778 /* scanInit.skyGridFile = uvar_skyGridFile; */
1779 scanInit.skyRegionString = ( CHAR * )LALCalloc( 1, strlen( skyRegion ) + 1 );
1780 strcpy( scanInit.skyRegionString, skyRegion );
1781 /* scanInit.Freq = usefulParams.spinRange_midTime.fkdot[0] + usefulParams.spinRange_midTime.fkdotBand[0]; */
1782
1783 scanInit.numSkyPartitions = numSkyPartitions;
1784 scanInit.partitionIndex = partitionIndex;
1785
1786 /* set up the grid */
1787 TRY( InitDopplerSkyScan( status->statusPtr, &thisScan, &scanInit ), status );
1788
1789 nSkyPatches = out->numSkyPatches = thisScan.numSkyGridPoints;
1790
1791 out->alpha = ( REAL8 * )LALCalloc( 1, nSkyPatches * sizeof( REAL8 ) );
1792 out->delta = ( REAL8 * )LALCalloc( 1, nSkyPatches * sizeof( REAL8 ) );
1793 out->alphaSize = ( REAL8 * )LALCalloc( 1, nSkyPatches * sizeof( REAL8 ) );
1794 out->deltaSize = ( REAL8 * )LALCalloc( 1, nSkyPatches * sizeof( REAL8 ) );
1795
1796 /* loop over skygrid points */
1797 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1798
1799 skyCounter = 0;
1800 while ( thisScan.state != STATE_FINISHED ) {
1801
1802 out->alpha[skyCounter] = dopplerpos.Alpha;
1803 out->delta[skyCounter] = dopplerpos.Delta;
1804 out->alphaSize[skyCounter] = dAlpha;
1805 out->deltaSize[skyCounter] = dDelta;
1806
1807 /* if ((dopplerpos.Delta>0) && (dopplerpos.Delta < atan(4*LAL_PI/dAlpha/dDelta) ))
1808 out->alphaSize[skyCounter] = dAlpha*cos(dopplerpos.Delta -0.5*dDelta)/cos(dopplerpos.Delta);
1809
1810 if ((dopplerpos.Delta<0) && (dopplerpos.Delta > -atan(4*LAL_PI/dAlpha/dDelta) ))
1811 out->alphaSize[skyCounter] = dAlpha*cos(dopplerpos.Delta +0.5*dDelta)/cos(dopplerpos.Delta);
1812 */
1813
1814 XLALNextDopplerSkyPos( &dopplerpos, &thisScan );
1815 skyCounter++;
1816
1817 } /* end while loop over skygrid */
1818
1819 /* free dopplerscan stuff */
1820 TRY( FreeDopplerSkyScan( status->statusPtr, &thisScan ), status );
1821 if ( scanInit.skyRegionString ) {
1822 LALFree( scanInit.skyRegionString );
1823 }
1824
1825 } else {
1826
1827 /* read skypatch info */
1828 {
1829 FILE *fpsky = NULL;
1830 INT4 r;
1831 REAL8 temp1, temp2, temp3, temp4;
1832
1833 ASSERT( skyFileName, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1834
1835 if ( ( fpsky = fopen( skyFileName, "r" ) ) == NULL ) {
1836 ABORT( status, DRIVEHOUGHCOLOR_EFILE, DRIVEHOUGHCOLOR_MSGEFILE );
1837 }
1838
1839 nSkyPatches = 0;
1840 do {
1841 r = fscanf( fpsky, "%lf%lf%lf%lf\n", &temp1, &temp2, &temp3, &temp4 );
1842 /* make sure the line has the right number of entries or is EOF */
1843 if ( r == 4 ) {
1844 nSkyPatches++;
1845 }
1846 } while ( r != EOF );
1847 rewind( fpsky );
1848
1849 out->numSkyPatches = nSkyPatches;
1850 out->alpha = ( REAL8 * )LALCalloc( nSkyPatches, sizeof( REAL8 ) );
1851 out->delta = ( REAL8 * )LALCalloc( nSkyPatches, sizeof( REAL8 ) );
1852 out->alphaSize = ( REAL8 * )LALCalloc( nSkyPatches, sizeof( REAL8 ) );
1853 out->deltaSize = ( REAL8 * )LALCalloc( nSkyPatches, sizeof( REAL8 ) );
1854
1855 for ( skyCounter = 0; skyCounter < nSkyPatches; skyCounter++ ) {
1856 r = fscanf( fpsky, "%lf%lf%lf%lf\n", out->alpha + skyCounter, out->delta + skyCounter,
1857 out->alphaSize + skyCounter, out->deltaSize + skyCounter );
1858 }
1859
1860 fclose( fpsky );
1861 } /* end skypatchfile reading block */
1862 } /* end setting up of skypatches */
1863
1864
1865 DETATCHSTATUSPTR( status );
1866
1867 /* normal exit */
1868 RETURN( status );
1869
1870}
1871
1872
1873void GetAMWeights( LALStatus *status,
1874 REAL8Vector *out,
1875 MultiDetectorStateSeries *mdetStates,
1876 REAL8 alpha,
1877 REAL8 delta )
1878{
1879
1880 MultiAMCoeffs *multiAMcoef = NULL;
1881 UINT4 iIFO, iSFT, k, numsft, numifo;
1882 REAL8 a, b;
1883 SkyPosition skypos;
1884
1885 INITSTATUS( status );
1886 ATTATCHSTATUSPTR( status );
1887
1888 /* get the amplitude modulation coefficients */
1889 skypos.longitude = alpha;
1890 skypos.latitude = delta;
1891 skypos.system = COORDINATESYSTEM_EQUATORIAL;
1892 XLAL_CHECK_LAL( status, ( multiAMcoef = XLALComputeMultiAMCoeffs( mdetStates, NULL, skypos ) ) != NULL, XLAL_EFUNC );
1893
1894 numifo = mdetStates->length;
1895
1896 /* loop over the weights and multiply them by the appropriate
1897 AM coefficients */
1898 for ( k = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
1899
1900 numsft = mdetStates->data[iIFO]->length;
1901
1902 for ( iSFT = 0; iSFT < numsft; iSFT++, k++ ) {
1903
1904 a = multiAMcoef->data[iIFO]->a->data[iSFT];
1905 b = multiAMcoef->data[iIFO]->b->data[iSFT];
1906 out->data[k] *= ( a * a + b * b );
1907
1908 } /* loop over SFTs */
1909
1910 } /* loop over IFOs */
1911
1912 TRY( LALHOUGHNormalizeWeights( status->statusPtr, out ), status );
1913
1914 XLALDestroyMultiAMCoeffs( multiAMcoef );
1915
1916
1917 DETATCHSTATUSPTR( status );
1918
1919 /* normal exit */
1920 RETURN( status );
1921
1922}
1923
1924
1925
1926void SelectBestStuff( LALStatus *status,
1927 BestVariables *out,
1928 BestVariables *in,
1929 UINT4 mObsCohBest )
1930{
1931
1932 size_t *ind = NULL;
1933 UINT4 k, mObsCoh;
1934
1935 INITSTATUS( status );
1936 ATTATCHSTATUSPTR( status );
1937
1938 /* check consistency of input */
1939 ASSERT( out, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1940
1941 ASSERT( in, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1942 ASSERT( in->length > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1943 mObsCoh = in->length;
1944
1945 ASSERT( in->weightsV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1946 ASSERT( in->weightsV->length == mObsCoh, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1947
1948 ASSERT( in->timeDiffV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1949 ASSERT( in->timeDiffV->length == mObsCoh, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1950
1951 ASSERT( in->velV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
1952 ASSERT( in->velV->length == mObsCoh, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1953
1954 ASSERT( mObsCohBest > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1955 ASSERT( mObsCohBest <= in->length, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
1956
1957 /* memory allocation for output -- use reallocs because this
1958 function may be called within the loop over sky positions, so memory might
1959 already have been allocated previously */
1960 out->length = mObsCohBest;
1961
1962
1963 if ( out->weightsV == NULL ) {
1964 out->weightsV = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
1965 out->weightsV->length = mObsCohBest;
1966 out->weightsV->data = ( REAL8 * )LALCalloc( 1, mObsCohBest * sizeof( REAL8 ) );
1967 }
1968
1969 if ( out->timeDiffV == NULL ) {
1970 out->timeDiffV = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
1971 out->timeDiffV->length = mObsCohBest;
1972 out->timeDiffV->data = ( REAL8 * )LALCalloc( 1, mObsCohBest * sizeof( REAL8 ) );
1973 }
1974
1975 if ( out->velV == NULL ) {
1976 out->velV = ( REAL8Cart3CoorVector * )LALCalloc( 1, sizeof( REAL8Cart3CoorVector ) );
1977 out->velV->length = mObsCohBest;
1978 out->velV->data = ( REAL8Cart3Coor * )LALCalloc( 1, mObsCohBest * sizeof( REAL8Cart3Coor ) );
1979 }
1980
1981 if ( out->pgV == NULL ) {
1982 out->pgV = ( HOUGHPeakGramVector * )LALCalloc( 1, sizeof( HOUGHPeakGramVector ) );
1983 out->pgV->length = mObsCohBest;
1984 out->pgV->pg = ( HOUGHPeakGram * )LALCalloc( 1, mObsCohBest * sizeof( HOUGHPeakGram ) );
1985 }
1986
1987 ind = LALCalloc( 1, mObsCohBest * sizeof( size_t ) );
1988 gsl_sort_largest_index( ind, mObsCohBest, in->weightsV->data, 1, mObsCoh );
1989
1990 for ( k = 0; k < mObsCohBest; k++ ) {
1991
1992 out->weightsV->data[k] = in->weightsV->data[ind[k]];
1993 out->timeDiffV->data[k] = in->timeDiffV->data[ind[k]];
1994
1995 out->velV->data[k].x = in->velV->data[ind[k]].x;
1996 out->velV->data[k].y = in->velV->data[ind[k]].y;
1997 out->velV->data[k].z = in->velV->data[ind[k]].z;
1998
1999 /* this copies the pointers to the peakgram data from the input */
2000 out->pgV->pg[k] = in->pgV->pg[ind[k]];
2001
2002 }
2003
2004 /* gsl_sort_index( index, timeDiffV.data, 1, mObsCohBest); */
2005
2006 LALFree( ind );
2007
2008 DETATCHSTATUSPTR( status );
2009
2010 /* normal exit */
2011 RETURN( status );
2012
2013}
2014
2015
2016/* copy data in BestVariables struct */
2017void DuplicateBestStuff( LALStatus *status,
2018 BestVariables *out,
2019 BestVariables *in )
2020{
2021
2022 UINT4 mObsCoh;
2023
2024 INITSTATUS( status );
2025 ATTATCHSTATUSPTR( status );
2026
2027 /* check consistency of input */
2028 ASSERT( out, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2029
2030 ASSERT( in, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2031 ASSERT( in->length > 0, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2032 mObsCoh = in->length;
2033
2034 ASSERT( in->weightsV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2035 ASSERT( in->weightsV->length == mObsCoh, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2036
2037 ASSERT( in->timeDiffV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2038 ASSERT( in->timeDiffV->length == mObsCoh, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2039
2040 ASSERT( in->velV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2041 ASSERT( in->velV->length == mObsCoh, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2042
2043 /* memory allocation for output -- check output is null because
2044 function may be called within the loop over sky positions, so memory might
2045 already have been allocated previously */
2046 out->length = mObsCoh;
2047
2048 if ( out->weightsV == NULL ) {
2049 out->weightsV = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
2050 out->weightsV->length = mObsCoh;
2051 out->weightsV->data = ( REAL8 * )LALCalloc( 1, mObsCoh * sizeof( REAL8 ) );
2052 }
2053
2054 if ( out->timeDiffV == NULL ) {
2055 out->timeDiffV = ( REAL8Vector * )LALCalloc( 1, sizeof( REAL8Vector ) );
2056 out->timeDiffV->length = mObsCoh;
2057 out->timeDiffV->data = ( REAL8 * )LALCalloc( 1, mObsCoh * sizeof( REAL8 ) );
2058 }
2059
2060 if ( out->velV == NULL ) {
2061 out->velV = ( REAL8Cart3CoorVector * )LALCalloc( 1, sizeof( REAL8Cart3CoorVector ) );
2062 out->velV->length = mObsCoh;
2063 out->velV->data = ( REAL8Cart3Coor * )LALCalloc( 1, mObsCoh * sizeof( REAL8Cart3Coor ) );
2064 }
2065
2066 if ( out->pgV == NULL ) {
2067 out->pgV = ( HOUGHPeakGramVector * )LALCalloc( 1, sizeof( HOUGHPeakGramVector ) );
2068 out->pgV->length = mObsCoh;
2069 out->pgV->pg = ( HOUGHPeakGram * )LALCalloc( 1, mObsCoh * sizeof( HOUGHPeakGram ) );
2070 }
2071
2072 memcpy( out->weightsV->data, in->weightsV->data, mObsCoh * sizeof( REAL8 ) );
2073 memcpy( out->timeDiffV->data, in->timeDiffV->data, mObsCoh * sizeof( REAL8 ) );
2074 memcpy( out->velV->data, in->velV->data, mObsCoh * sizeof( REAL8Cart3Coor ) );
2075 memcpy( out->pgV->pg, in->pgV->pg, mObsCoh * sizeof( HOUGHPeakGram ) );
2076
2077 DETATCHSTATUSPTR( status );
2078
2079 /* normal exit */
2080 RETURN( status );
2081
2082}
2083
2084
2085
2086
2087
2088/** helper function for creating LUT vector */
2089void LALHOUGHCreateLUTVector( LALStatus *status,
2090 HOUGHptfLUTVector *lutV,
2091 UINT4 length )
2092{
2093
2094 INITSTATUS( status );
2095 ATTATCHSTATUSPTR( status );
2096
2097 /* check input pars are ok */
2098 if ( lutV == NULL ) {
2099 ABORT( status, DRIVEHOUGHCOLOR_ENONULL, DRIVEHOUGHCOLOR_MSGENONULL );
2100 }
2101
2102 if ( lutV->lut != NULL ) {
2103 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2104 }
2105
2106 if ( length <= 0 ) {
2107 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2108 }
2109
2110 /* allocate memory */
2111 lutV->length = length;
2112 lutV->lut = ( HOUGHptfLUT * )LALCalloc( lutV->length, sizeof( HOUGHptfLUT ) );
2113
2114 DETATCHSTATUSPTR( status );
2115
2116 /* normal exit */
2117 RETURN( status );
2118
2119}
2120
2121
2122void LALHOUGHCreateLUTs( LALStatus *status,
2123 HOUGHptfLUTVector *lutV,
2124 UINT2 maxNBins,
2125 UINT2 maxNBorders,
2126 UINT2 ySide )
2127{
2128
2129 UINT4 j, i;
2130
2131 INITSTATUS( status );
2132 ATTATCHSTATUSPTR( status );
2133
2134 /* check input pars are ok */
2135 if ( lutV == NULL ) {
2136 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2137 }
2138
2139 if ( lutV->lut == NULL ) {
2140 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2141 }
2142
2143 if ( maxNBins <= 0 ) {
2144 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2145 }
2146
2147 if ( maxNBorders <= 0 ) {
2148 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2149 }
2150
2151 if ( ySide <= 0 ) {
2152 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2153 }
2154
2155
2156 /* loop over luts and allocate memory */
2157 for ( j = 0; j < lutV->length; j++ ) {
2158 lutV->lut[j].maxNBins = maxNBins;
2159 lutV->lut[j].maxNBorders = maxNBorders;
2160 lutV->lut[j].border = ( HOUGHBorder * )LALCalloc( maxNBorders, sizeof( HOUGHBorder ) );
2161 lutV->lut[j].bin = ( HOUGHBin2Border * )LALCalloc( maxNBins, sizeof( HOUGHBin2Border ) );
2162
2163 for ( i = 0; i < maxNBorders; i++ ) {
2164 lutV->lut[j].border[i].ySide = ySide;
2165 lutV->lut[j].border[i].xPixel = ( COORType * )LALCalloc( ySide, sizeof( COORType ) );
2166 }
2167
2168 }
2169
2170 DETATCHSTATUSPTR( status );
2171
2172 /* normal exit */
2173 RETURN( status );
2174
2175}
2176
2177
2178void LALHOUGHDestroyLUTs( LALStatus *status,
2179 HOUGHptfLUTVector *lutV )
2180{
2181
2182 UINT4 j, i;
2183
2184 INITSTATUS( status );
2185 ATTATCHSTATUSPTR( status );
2186
2187 for ( j = 0; j < lutV->length ; ++j ) {
2188 for ( i = 0; i < lutV->lut[j].maxNBorders; ++i ) {
2189
2190 if ( lutV->lut[j].border[i].xPixel ) {
2191 LALFree( lutV->lut[j].border[i].xPixel );
2192 }
2193 }
2194
2195 if ( lutV->lut[j].border ) {
2196 LALFree( lutV->lut[j].border );
2197 }
2198
2199 if ( lutV->lut[j].bin ) {
2200 LALFree( lutV->lut[j].bin );
2201 }
2202 }
2203
2204 DETATCHSTATUSPTR( status );
2205
2206 /* normal exit */
2207 RETURN( status );
2208
2209}
2210
2211
2212void LALHOUGHCreatePHMDVS( LALStatus *status,
2213 PHMDVectorSequence *phmdVS,
2214 UINT4 length,
2215 UINT4 nfSize )
2216{
2217
2218 INITSTATUS( status );
2219 ATTATCHSTATUSPTR( status );
2220
2221 /* check input pars are ok */
2222 if ( phmdVS == NULL ) {
2223 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2224 }
2225
2226 if ( phmdVS->phmd != NULL ) {
2227 ABORT( status, DRIVEHOUGHCOLOR_ENONULL, DRIVEHOUGHCOLOR_MSGENONULL );
2228 }
2229
2230 if ( length <= 0 ) {
2231 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2232 }
2233
2234 if ( nfSize <= 0 ) {
2235 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2236 }
2237
2238 phmdVS->length = length;
2239 phmdVS->nfSize = nfSize;
2240 phmdVS->deltaF = 0; /* initialization */
2241
2242 phmdVS->phmd = ( HOUGHphmd * )LALCalloc( 1, length * nfSize * sizeof( HOUGHphmd ) );
2243
2244 DETATCHSTATUSPTR( status );
2245
2246 /* normal exit */
2247 RETURN( status );
2248
2249}
2250
2251
2252
2253void LALHOUGHCreatePHMDs( LALStatus *status,
2254 PHMDVectorSequence *phmdVS,
2255 UINT2 maxNBins,
2256 UINT2 maxNBorders,
2257 UINT2 ySide )
2258{
2259
2260 UINT4 j;
2261
2262 INITSTATUS( status );
2263 ATTATCHSTATUSPTR( status );
2264
2265 /* check input pars are ok */
2266 if ( phmdVS == NULL ) {
2267 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2268 }
2269
2270 if ( phmdVS->phmd == NULL ) {
2271 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2272 }
2273
2274 if ( maxNBins <= 0 ) {
2275 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2276 }
2277
2278 if ( maxNBorders <= 0 ) {
2279 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2280 }
2281
2282 if ( ySide <= 0 ) {
2283 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2284 }
2285
2286 for ( j = 0; j < phmdVS->length * phmdVS->nfSize; j++ ) {
2287 phmdVS->phmd[j].maxNBorders = maxNBorders;
2288 phmdVS->phmd[j].leftBorderP = ( HOUGHBorder ** )LALCalloc( maxNBorders, sizeof( HOUGHBorder * ) );
2289 phmdVS->phmd[j].rightBorderP = ( HOUGHBorder ** )LALCalloc( maxNBorders, sizeof( HOUGHBorder * ) );
2290
2291 phmdVS->phmd[j].ySide = ySide;
2292 phmdVS->phmd[j].firstColumn = ( UCHAR * )LALCalloc( ySide, sizeof( UCHAR ) );
2293 }
2294
2295
2296 DETATCHSTATUSPTR( status );
2297
2298 /* normal exit */
2299 RETURN( status );
2300
2301}
2302
2303
2304void LALHOUGHDestroyPHMDs( LALStatus *status,
2305 PHMDVectorSequence *phmdVS )
2306{
2307 UINT4 j;
2308
2309 INITSTATUS( status );
2310 ATTATCHSTATUSPTR( status );
2311
2312 for ( j = 0; j < phmdVS->length * phmdVS->nfSize; j++ ) {
2313 LALFree( phmdVS->phmd[j].leftBorderP );
2314 LALFree( phmdVS->phmd[j].rightBorderP );
2315 LALFree( phmdVS->phmd[j].firstColumn );
2316 }
2317
2318 DETATCHSTATUSPTR( status );
2319
2320 /* normal exit */
2321 RETURN( status );
2322
2323}
2324
2325
2326void LALHOUGHCreateHT( LALStatus *status,
2327 HOUGHMapTotal *ht,
2328 UINT2 xSide,
2329 UINT2 ySide )
2330{
2331
2332 INITSTATUS( status );
2333 ATTATCHSTATUSPTR( status );
2334
2335 /* check input pars are ok */
2336 if ( ht == NULL ) {
2337 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2338 }
2339
2340 if ( xSide <= 0 ) {
2341 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2342 }
2343
2344 if ( ySide <= 0 ) {
2345 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2346 }
2347
2348
2349 ht->xSide = xSide;
2350 ht->ySide = ySide;
2351
2352 ht->map = NULL;
2353 ht->map = ( HoughTT * )LALCalloc( xSide * ySide, sizeof( HoughTT ) );
2354
2355 ht->deltaF = 0; /* initialization */
2356 ht->mObsCoh = 0; /* initialization */
2357
2358 DETATCHSTATUSPTR( status );
2359
2360 /* normal exit */
2361 RETURN( status );
2362
2363}
2364
2365
2366
2367void LALHOUGHCreateFreqIndVector( LALStatus *status,
2368 UINT8FrequencyIndexVector *freqInd,
2369 UINT4 length,
2370 REAL8 deltaF )
2371{
2372
2373 INITSTATUS( status );
2374 ATTATCHSTATUSPTR( status );
2375
2376 /* check input pars are ok */
2377 if ( freqInd == NULL ) {
2378 ABORT( status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2379 }
2380
2381 if ( length <= 0 ) {
2382 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2383 }
2384
2385 if ( deltaF <= 0 ) {
2386 ABORT( status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2387 }
2388
2389 freqInd->deltaF = deltaF;
2390 freqInd->length = length;
2391 freqInd->data = NULL;
2392 freqInd->data = ( UINT8 * )LALCalloc( 1, freqInd->length * sizeof( UINT8 ) );
2393
2394 DETATCHSTATUSPTR( status );
2395
2396 /* normal exit */
2397 RETURN( status );
2398
2399}
2400
2401
2402
2403/** Get Hough candidates as a toplist */
2404void GetToplistFromHoughmap( LALStatus *status,
2405 toplist_t *list,
2406 HOUGHMapTotal *ht,
2407 HOUGHPatchGrid *patch,
2408 HOUGHDemodPar *parDem,
2409 REAL8 mean,
2410 REAL8 sigma )
2411{
2412 REAL8UnitPolarCoor sourceLocation;
2413 REAL8 deltaF, f0, fdot;
2414 INT8 f0Bin;
2415 INT4 i, j, xSide, ySide;
2416 FstatOutputEntry candidate;
2417
2418 INITSTATUS( status );
2419 ATTATCHSTATUSPTR( status );
2420
2421 deltaF = ht->deltaF;
2422 f0Bin = ht->f0Bin;
2423 f0 = f0Bin * deltaF;
2424
2425 fdot = ht->spinRes.data[0];
2426
2427 xSide = ht->xSide;
2428 ySide = ht->ySide;
2429
2430 candidate.Freq = f0;
2431 candidate.f1dot = fdot;
2432
2433 for ( i = 0; i < ySide; i++ ) {
2434 for ( j = 0; j < xSide; j++ ) {
2435 /* get sky location of pixel */
2436 TRY( LALStereo2SkyLocation( status->statusPtr, &sourceLocation,
2437 j, i, patch, parDem ), status );
2438
2439 candidate.Alpha = sourceLocation.alpha;
2440 candidate.Delta = sourceLocation.delta;
2441 candidate.Fstat = ( ht->map[i * xSide + j] - mean ) / sigma;
2442
2443 insert_into_fstat_toplist( list, candidate );
2444
2445 } /* end loop over xSide */
2446 } /* end loop over ySide */
2447
2448 DETATCHSTATUSPTR( status );
2449 RETURN( status );
2450
2451} /* end hough toplist selection */
2452
2453
2454/* *********************************************************************************************/
2455/* Print some extra information: HoughMaps, HoughStatistics, sigma */
2456int OpenExtraInfoFiles( CHAR *fileMaps,
2457 FILE **fp1_ptr,
2458 CHAR *filehisto,
2459 CHAR *dirname,
2460 CHAR *basename,
2461 INT4 ind )
2462{
2463 CHAR filestats[ MAXFILENAMELENGTH ];
2464 /*CHAR fileMaps[ MAXFILENAMELENGTH ];*/
2465 FILE *fp1 = NULL;
2466
2467
2468 /* --------------------------------------------------*/
2469 /* Create directory fnameout/skypatch_$j using mkdir if required */
2470 if ( CreateSkypatchDirs( filestats, dirname, ind ) ) {
2471 return DRIVEHOUGHCOLOR_EFILE;
2472 }
2473
2474 /* --------------------------------------------------*/
2475 /* Create the base filenames for the stats, histo and template files.*/
2476 strcat( filestats, basename );
2477 strcpy( filehisto, filestats );
2478 strcpy( fileMaps, filestats );
2479
2480 /* --------------------------------------------------*/
2481 /* Create and open the stats file for writing */
2482 strcat( filestats, "stats" );
2483 if ( ( fp1 = fopen( filestats, "w" ) ) == NULL ) {
2484 fprintf( stderr, "Unable to find file %s for writing\n", filestats );
2485 return DRIVEHOUGHCOLOR_EFILE;
2486 }
2487 setvbuf( fp1, ( char * )NULL, _IOLBF, 0 );
2488 *fp1_ptr = fp1;
2489 return ( 0 );
2490
2491} /* end OpenExtraInfoFiles */
2492
2493
2494/* *********************************************************************************************/
2495/** Print some extra information: HoughMaps, HoughStatistics, sigma */
2496int PrintExtraInfo( CHAR *fileMaps,
2497 FILE **fp1_ptr,
2498 INT4 iHmap,
2499 HOUGHMapTotal *ht,
2500 REAL8UnitPolarCoor *sourceLocation,
2501 HoughStats *stats,
2502 INT8 fBinSearch,
2503 REAL8 deltaF )
2504{
2505 FILE *fp1 = *fp1_ptr;
2506
2507 /* --------------------------------------------------*/
2508 /* Print results */
2509
2510 /* Printing Hough Maps */
2511 if ( PrintHmap2m_file( ht, fileMaps, iHmap ) ) {
2512 return 5;
2513 }
2514
2515 /* Printing Statistics */
2516 fprintf( fp1, "%d %f %f %f %f %f %f %f %g\n",
2517 iHmap, sourceLocation->alpha, sourceLocation->delta,
2518 ( REAL4 )stats->maxCount, ( REAL4 )stats->minCount, stats->avgCount, stats->stdDev,
2519 ( fBinSearch * deltaF ), ht->spinRes.data[0] );
2520
2521 return ( 0 );
2522
2523} /* end pirntExtraInfo */
2524
2525
2526
2527
2528/********************************************************************************/
2529/* Computing the frequency path with no mismatch */
2530/********************************************************************************/
2531void ComputeFoft_NM( LALStatus *status,
2532 REAL8Vector *foft,
2533 HoughTemplate *pulsarTemplate,
2534 REAL8Vector *timeDiffV,
2535 REAL8Cart3CoorVector *velV )
2536{
2537
2538 INT4 mObsCoh;
2539 REAL8 f0new, vcProdn, timeDiffN;
2540 REAL8 sourceDelta, sourceAlpha, cosDelta;
2541 INT4 j, i, nspin, factorialN;
2542 REAL8Cart3Coor sourceLocation;
2543
2544 /* --------------------------------------------- */
2545 INITSTATUS( status );
2546 ATTATCHSTATUSPTR( status );
2547
2548 /* Make sure the arguments are not NULL: */
2549 ASSERT( foft, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2550 ASSERT( pulsarTemplate, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2551 ASSERT( timeDiffV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2552 ASSERT( velV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2553
2554 ASSERT( foft->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2555 ASSERT( timeDiffV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2556 ASSERT( velV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2557
2558 sourceDelta = pulsarTemplate->latitude;
2559 sourceAlpha = pulsarTemplate->longitude;
2560 cosDelta = cos( sourceDelta );
2561
2562 sourceLocation.x = cosDelta * cos( sourceAlpha );
2563 sourceLocation.y = cosDelta * sin( sourceAlpha );
2564 sourceLocation.z = sin( sourceDelta );
2565
2566 mObsCoh = foft->length;
2567 nspin = pulsarTemplate->spindown.length;
2568
2569 for ( j = 0; j < mObsCoh; ++j ) { /* loop for all different time stamps */
2570 vcProdn = velV->data[j].x * sourceLocation.x
2571 + velV->data[j].y * sourceLocation.y
2572 + velV->data[j].z * sourceLocation.z;
2573 f0new = pulsarTemplate->f0;
2574 factorialN = 1;
2575 timeDiffN = timeDiffV->data[j];
2576
2577 for ( i = 0; i < nspin; ++i ) { /* loop for spin-down values */
2578 factorialN *= ( i + 1 );
2579 f0new += pulsarTemplate->spindown.data[i] * timeDiffN / factorialN;
2580 timeDiffN *= timeDiffN;
2581 }
2582 foft->data[j] = f0new * ( 1.0 + vcProdn );
2583 /*foft->data[j] = floor(f0new*1800+0.5) * (1.0 +vcProdn)/1800;*/
2584 }
2585
2586
2587 DETATCHSTATUSPTR( status );
2588 /* normal exit */
2589 RETURN( status );
2590
2591} /* end Computefoft */
2592
2593
2594
2595
2596/* *********************************************************************/
2597
2598void SplitSFTs( LALStatus *status,
2599 REAL8Vector *weightsV,
2600 HoughParamsTest *chi2Params )
2601{
2602
2603 UINT4 j = 0; /* index of each block. It runs betwen 0 and p */
2604 UINT4 iSFT;
2605 REAL8 *weights_ptr; /* pointer to weightsV.data */
2606 REAL8 sumWeightpMax; /* Value of sumWeight we want to fix in each set of SFTs */
2607 UINT4 numberSFT; /* Counter with the # of SFTs in each block */
2608 UINT4 mObsCoh, p;
2609 REAL8 partialsumWeightp, partialsumWeightSquarep;
2610
2611 /* --------------------------------------------- */
2612 INITSTATUS( status );
2613 ATTATCHSTATUSPTR( status );
2614
2615 /* Make sure the arguments are not NULL: */
2616 ASSERT( weightsV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2617 ASSERT( chi2Params, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2618
2619 ASSERT( weightsV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2620 ASSERT( chi2Params->length, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2621 ASSERT( chi2Params->numberSFTp, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2622 ASSERT( chi2Params->sumWeight, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2623 ASSERT( chi2Params->sumWeightSquare, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2624 ASSERT( chi2Params->length < weightsV->length, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2625
2626 mObsCoh = weightsV->length;
2627 p = chi2Params->length;
2628
2629 sumWeightpMax = ( REAL8 )( mObsCoh ) / p; /* Compute the value of the sumWeight we want to fix in each set of SFT's */
2630 weights_ptr = weightsV->data; /* Make the pointer to point to the first position of the vector weightsV.data */
2631
2632 iSFT = 0;
2633 for ( j = 0; j < p; j++ ) {
2634
2635 partialsumWeightSquarep = 0;
2636 partialsumWeightp = 0;
2637
2638 for ( numberSFT = 0; ( partialsumWeightp < sumWeightpMax ) && ( iSFT < mObsCoh ); numberSFT++, iSFT++ ) {
2639
2640 partialsumWeightp += *weights_ptr;
2641 partialsumWeightSquarep += ( *weights_ptr ) * ( *weights_ptr );
2642 weights_ptr++;
2643
2644 } /* loop over SFTs */
2645
2646 ASSERT( ( UINT4 )j < p, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2647
2648 chi2Params->numberSFTp[j] = numberSFT;
2649 chi2Params->sumWeight[j] = partialsumWeightp;
2650 chi2Params->sumWeightSquare[j] = partialsumWeightSquarep;
2651
2652 } /* loop over the p blocks of data */
2653
2654 ASSERT( iSFT == mObsCoh, status, DRIVEHOUGHCOLOR_EBAD, DRIVEHOUGHCOLOR_MSGEBAD );
2655
2656 DETATCHSTATUSPTR( status );
2657 /* normal exit */
2658 RETURN( status );
2659
2660}/* SplitSFTs */
2661
2662
2663
2664
2665/**********************************************************************/
2666
2667void ComputeandPrintChi2( LALStatus *status,
2668 toplist_t *tl,
2669 REAL8Vector *timeDiffV,
2670 REAL8Cart3CoorVector *velV,
2671 INT4 skyCounter,
2672 INT4 nSkyPatches,
2673 INT4 p,
2674 REAL8 alphaPeak,
2675 MultiDetectorStateSeries *mdetStates,
2676 REAL8Vector *weightsNoise,
2677 UCHARPeakGramVector *upgV /* Expanded (UCHAR) peakgrams */ )
2678{
2679 HoughTemplate pulsarTemplate;
2680 UINT4 mObsCoh;
2681 UINT4 k, i, j, ii, numberSFTp ;
2682 INT4 ind;
2683 REAL8 sumWeightSquare, meanN, sigmaN;
2684 REAL8 eta, numberCount;
2685 REAL8 nj, sumWeightj, sumWeightSquarej;
2686 FstatOutputEntry readTopList;
2687 REAL8Vector foft;
2688 REAL8Vector weightsV;
2689 REAL8 timeBase;
2690 UINT4 sftFminBin ;
2691 FILE *fpChi2 = NULL;
2692 REAL8 oldSig;
2693
2694 /* Chi2Test parameters */
2695 HoughParamsTest chi2Params;
2696 REAL8 numberCountTotal; /* Sum over all the numberCounts */
2697 REAL8 chi2;
2698 REAL8Vector numberCountV; /* Vector with the number count of each block inside */
2699
2700 /* --------------------------------------------- */
2701 INITSTATUS( status );
2702 ATTATCHSTATUSPTR( status );
2703
2704 /* Make sure the arguments are not NULL: */
2705 ASSERT( timeDiffV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2706 ASSERT( velV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2707 ASSERT( mdetStates, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2708 ASSERT( weightsNoise, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2709 ASSERT( upgV, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2710
2711 ASSERT( weightsNoise->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2712 ASSERT( timeDiffV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2713 ASSERT( velV->data, status, DRIVEHOUGHCOLOR_ENULL, DRIVEHOUGHCOLOR_MSGENULL );
2714
2715 mObsCoh = velV->length;
2716 timeBase = upgV->upg[0].timeBase;
2717 sftFminBin = upgV->upg[0].fminBinIndex;
2718
2719 /* memory for weightsV */
2720 weightsV.length = mObsCoh;
2721 weightsV.data = ( REAL8 * )LALCalloc( 1, mObsCoh * sizeof( REAL8 ) );
2722
2723 /* memory for chi2Params */
2724 chi2Params.length = p;
2725 chi2Params.numberSFTp = NULL;
2726 chi2Params.sumWeight = NULL;
2727 chi2Params.sumWeightSquare = NULL;
2728 chi2Params.numberSFTp = ( UINT4 * )LALMalloc( p * sizeof( UINT4 ) );
2729 chi2Params.sumWeight = ( REAL8 * )LALMalloc( p * sizeof( REAL8 ) );
2730 chi2Params.sumWeightSquare = ( REAL8 * )LALMalloc( p * sizeof( REAL8 ) );
2731
2732 /* memory for number Count Vector */
2733 numberCountV.length = p;
2734 numberCountV.data = NULL;
2735 numberCountV.data = ( REAL8 * )LALMalloc( p * sizeof( REAL8 ) );
2736
2737 /* Memory for one spindown */
2738 pulsarTemplate.spindown.length = 1 ;
2739 pulsarTemplate.spindown.data = NULL;
2740 pulsarTemplate.spindown.data = ( REAL8 * )LALMalloc( sizeof( REAL8 ) );
2741
2742 /* Memory for f(t) vector */
2743 foft.length = mObsCoh;
2744 foft.data = NULL;
2745 foft.data = ( REAL8 * )LALMalloc( mObsCoh * sizeof( REAL8 ) );
2746
2747 /* Open file to write the toplist with 2 new columns: significance and chi2 */
2748
2749 fpChi2 = fopen( "hough_top.dat", "w" );
2750 char path[512];
2751 snprintf( path, sizeof( path ), "houghtop_chi2_%d_%d.dat", skyCounter + 1, nSkyPatches );
2752 fpChi2 = fopen( path, "w" );
2753
2754
2755 /* ----------------------------------------------------------------------------------*/
2756 /* Loop over all the elements in the TopList */
2757
2758 for ( i = 0; i < tl->elems; i++ ) {
2759
2760 readTopList = *( ( FstatOutputEntry * )( void * )( tl->heap[i] ) );
2761
2762 /* Copy template parameters from the TopList */
2763 pulsarTemplate.f0 = readTopList.Freq ;
2764 pulsarTemplate.spindown.data[0] = readTopList.f1dot ;
2765 pulsarTemplate.latitude = readTopList.Delta ;
2766 pulsarTemplate.longitude = readTopList.Alpha ;
2767 oldSig = readTopList.Fstat;
2768
2769 /* copy noise weights */
2770 memcpy( weightsV.data, weightsNoise->data, mObsCoh * sizeof( REAL8 ) );
2771
2772 /* calculate amplitude modulation weights */
2773 TRY( GetAMWeights( status->statusPtr, &weightsV, mdetStates, pulsarTemplate.longitude, pulsarTemplate.latitude ), status );
2774
2775
2776 /**********************************************************************************/
2777 /* Split the SFTs into p blocks and calculate the number count in each block */
2778 /**********************************************************************************/
2779
2780 /* compute mean and sigma for noise only */
2781 /* first calculate the sum of the weights squared */
2782 sumWeightSquare = 0.0;
2783
2784 for ( k = 0; k < ( UINT4 )mObsCoh; k++ ) {
2785 sumWeightSquare += weightsV.data[k] * weightsV.data[k];
2786 }
2787
2788 meanN = mObsCoh * alphaPeak;
2789 sigmaN = sqrt( sumWeightSquare * alphaPeak * ( 1.0 - alphaPeak ) );
2790
2791 /* the received frequency as a function of time */
2792 TRY( ComputeFoft_NM( status->statusPtr, &foft, &pulsarTemplate, timeDiffV, velV ), status );
2793
2794
2795 /* Split the SFTs into p blocs */
2796 TRY( SplitSFTs( status->statusPtr, &weightsV, &chi2Params ), status );
2797
2798
2799 /* loop over SFT, generate peakgram and get number count */
2800
2801 j = 0;
2802
2803 for ( k = 0 ; k < ( UINT4 )p ; k++ ) {
2804
2805 numberSFTp = chi2Params.numberSFTp[k];
2806 numberCount = 0;
2807
2808 for ( ii = 0 ; ( ii < numberSFTp ) ; ii++ ) {
2809
2810 ind = floor( foft.data[j] * timeBase - sftFminBin + 0.5 );
2811
2812 numberCount += upgV->upg[j].data[ind] * weightsV.data[j];
2813
2814 j++;
2815
2816 } /* loop over SFTs */
2817
2818 numberCountV.data[k] = numberCount;
2819
2820 }/* loop over blocks */
2821
2822
2823 /*-----------------------------*/
2824 /* Chi2 Test */
2825
2826 numberCountTotal = 0;
2827 chi2 = 0;
2828
2829 for ( k = 0; k < ( UINT4 )p ; k++ ) {
2830 numberCountTotal += numberCountV.data[k];
2831 }
2832
2833 eta = numberCountTotal / mObsCoh;
2834
2835 for ( j = 0 ; j < ( ( UINT4 )p ) ; j++ ) {
2836
2837 nj = numberCountV.data[j];
2838 sumWeightj = chi2Params.sumWeight[j];
2839 sumWeightSquarej = chi2Params.sumWeightSquare[j];
2840
2841 chi2 += ( nj - sumWeightj * eta ) * ( nj - sumWeightj * eta ) / ( sumWeightSquarej * eta * ( 1 - eta ) );
2842 }
2843
2844 setvbuf( fpChi2, ( char * )NULL, _IOLBF, 0 );
2845 fprintf( fpChi2, "%.13g %.7g %.7g %.5g %.6g %.6g %.7g \n", pulsarTemplate.f0, pulsarTemplate.longitude, pulsarTemplate.latitude, pulsarTemplate.spindown.data[0], ( numberCountTotal - meanN ) / sigmaN, oldSig, chi2 );
2846
2847 /*-----------------------------*/
2848
2849 } /* End of loop over top list elements */
2850
2851 LALFree( pulsarTemplate.spindown.data );
2852 LALFree( foft.data );
2853 LALFree( numberCountV.data );
2854 LALFree( chi2Params.numberSFTp );
2855 LALFree( chi2Params.sumWeight );
2856 LALFree( chi2Params.sumWeightSquare );
2857 LALFree( weightsV.data );
2858 weightsV.data = NULL;
2859
2860 fclose( fpChi2 );
2861
2862
2863 DETATCHSTATUSPTR( status );
2864
2865 /* normal exit */
2866 RETURN( status );
2867
2868} /* End of ComputeChi2 */
2869
2870
2871
2872/************************************************************************************/
2873/** Loop over SFTs and set a threshold to get peakgrams. SFTs must be normalized. */
2874/** This function will create a vector with the uncompressed PeakGrams in it */
2875/** (this is necesary if we want to compute the chi2 later ) */
2876void GetPeakGramFromMultSFTVector_NondestroyPg1( LALStatus *status, /**< pointer to LALStatus structure */
2877 HOUGHPeakGramVector *out, /**< Output compressed peakgrams */
2878 UCHARPeakGramVector *upgV, /**< Output uncompressed peakgrams */
2879 MultiSFTVector *in, /**< Input SFTs */
2880 REAL8 thr /**< Threshold on SFT power */ )
2881{
2882 SFTtype *sft;
2883 INT4 nPeaks;
2884 UINT4 iIFO, iSFT, numsft, numifo, j, binsSFT;
2885
2886 INITSTATUS( status );
2887 ATTATCHSTATUSPTR( status );
2888
2889 numifo = in->length;
2890 binsSFT = in->data[0]->data->data->length;
2891
2892 /* loop over sfts and select peaks */
2893 for ( j = 0, iIFO = 0; iIFO < numifo; iIFO++ ) {
2894
2895 numsft = in->data[iIFO]->length;
2896
2897 for ( iSFT = 0; iSFT < numsft; iSFT++, j++ ) {
2898
2899 sft = in->data[iIFO]->data + iSFT;
2900
2901 /* Store the expanded PeakGrams */
2902 upgV->upg[j].length = binsSFT;
2903 upgV->upg[j].data = NULL;
2904 upgV->upg[j].data = ( UCHAR * )LALCalloc( 1, binsSFT * sizeof( UCHAR ) );
2905
2906 TRY( SFTtoUCHARPeakGram( status->statusPtr, &( upgV->upg[j] ), sft, thr ), status );
2907
2908 nPeaks = upgV->upg[j].nPeaks;
2909
2910 /* compress peakgram */
2911 out->pg[j].length = nPeaks;
2912 out->pg[j].peak = NULL;
2913 out->pg[j].peak = ( INT4 * )LALCalloc( 1, nPeaks * sizeof( INT4 ) );
2914
2915 TRY( LALUCHAR2HOUGHPeak( status->statusPtr, &( out->pg[j] ), &( upgV->upg[j] ) ), status );
2916
2917 } /* loop over SFTs */
2918 } /* loop over IFOs */
2919
2920
2921
2922 DETATCHSTATUSPTR( status );
2923
2924 /* normal exit */
2925 RETURN( status );
2926
2927}
InitDopplerSkyScan
void InitDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan, const DopplerSkyScanInit *init)
Definition: DopplerScan.c:273
FreeDopplerSkyScan
void FreeDopplerSkyScan(LALStatus *status, DopplerSkyScanState *skyScan)
Definition: DopplerScan.c:324
XLALNextDopplerSkyPos
int XLALNextDopplerSkyPos(PulsarDopplerParams *pos, DopplerSkyScanState *skyScan)
NextDopplerSkyPos(): step through sky-grid return 0 = OK, -1 = ERROR.
Definition: DopplerScan.c:127
STATE_FINISHED
@ STATE_FINISHED
all templates have been read
Definition: DopplerScan.h:84
GRID_ISOTROPIC
@ GRID_ISOTROPIC
approximately isotropic sky-grid
Definition: DopplerScan.h:92
DriveHoughColor.h
Header file for non-demodulated Hough search.
DRIVEHOUGHCOLOR_ENONULL
#define DRIVEHOUGHCOLOR_ENONULL
Definition: DriveHoughColor.h:97
DRIVEHOUGHCOLOR_MSGEBAD
#define DRIVEHOUGHCOLOR_MSGEBAD
Definition: DriveHoughColor.h:102
DRIVEHOUGHCOLOR_MSGENULL
#define DRIVEHOUGHCOLOR_MSGENULL
Definition: DriveHoughColor.h:105
DRIVEHOUGHCOLOR_EFILE
#define DRIVEHOUGHCOLOR_EFILE
Definition: DriveHoughColor.h:94
DRIVEHOUGHCOLOR_ENULL
#define DRIVEHOUGHCOLOR_ENULL
Definition: DriveHoughColor.h:96
DRIVEHOUGHCOLOR_EBAD
#define DRIVEHOUGHCOLOR_EBAD
Definition: DriveHoughColor.h:93
DRIVEHOUGHCOLOR_MSGEFILE
#define DRIVEHOUGHCOLOR_MSGEFILE
Definition: DriveHoughColor.h:103
DRIVEHOUGHCOLOR_MSGENONULL
#define DRIVEHOUGHCOLOR_MSGENONULL
Definition: DriveHoughColor.h:106
ComputeandPrintChi2
void ComputeandPrintChi2(LALStatus *status, toplist_t *tl, REAL8Vector *timeDiffV, REAL8Cart3CoorVector *velV, INT4 skyCounter, INT4 nSkyPatches, INT4 p, REAL8 alphaPeak, MultiDetectorStateSeries *mdetStates, REAL8Vector *weightsNoise, UCHARPeakGramVector *upgV)
Definition: DriveHoughMulti.c:2666
SKYFILE
#define SKYFILE
Definition: DriveHoughMulti.c:116
main
int main(int argc, char *argv[])
Definition: DriveHoughMulti.c:237
PrintExtraInfo
int PrintExtraInfo(CHAR *fileMaps, FILE **fp1_ptr, INT4 iHmap, HOUGHMapTotal *ht, REAL8UnitPolarCoor *sourceLocation, HoughStats *stats, INT8 fBinSearch, REAL8 deltaF)
Print some extra information: HoughMaps, HoughStatistics, sigma.
Definition: DriveHoughMulti.c:2495
SplitSFTs
void SplitSFTs(LALStatus *status, REAL8Vector *weightsV, HoughParamsTest *chi2Params)
Definition: DriveHoughMulti.c:2597
NSPINUP
#define NSPINUP
Definition: DriveHoughMulti.c:120
uvar_EnableToplistPatch
BOOLEAN uvar_EnableToplistPatch
Definition: DriveHoughMulti.c:102
GetPeakGramFromMultSFTVector_NondestroyPg1
void GetPeakGramFromMultSFTVector_NondestroyPg1(LALStatus *status, HOUGHPeakGramVector *out, UCHARPeakGramVector *upgV, MultiSFTVector *in, REAL8 thr)
Loop over SFTs and set a threshold to get peakgrams.
Definition: DriveHoughMulti.c:2875
SelectBestStuff
void SelectBestStuff(LALStatus *status, BestVariables *out, BestVariables *in, UINT4 mObsCohBest)
Definition: DriveHoughMulti.c:1925
F0
#define F0
Definition: DriveHoughMulti.c:117
NFSIZE
#define NFSIZE
Definition: DriveHoughMulti.c:119
GetToplistFromHoughmap
void GetToplistFromHoughmap(LALStatus *status, toplist_t *list, HOUGHMapTotal *ht, HOUGHPatchGrid *patch, HOUGHDemodPar *parDem, REAL8 mean, REAL8 sigma)
Get Hough candidates as a toplist.
Definition: DriveHoughMulti.c:2403
GetAMWeights
void GetAMWeights(LALStatus *status, REAL8Vector *out, MultiDetectorStateSeries *mdetStates, REAL8 alpha, REAL8 delta)
Definition: DriveHoughMulti.c:1872
THRESHOLD
#define THRESHOLD
Definition: DriveHoughMulti.c:115
SetUpSkyPatches
void SetUpSkyPatches(LALStatus *status, HoughSkyPatchesInfo *out, CHAR *skyFileName, CHAR *skyRegion, REAL8 dAlpha, REAL8 dDelta, INT4 numSkyPartitions, INT4 partitionIndex)
Set up location of skypatch centers and sizes If user specified skyRegion then use DopplerScan functi...
Definition: DriveHoughMulti.c:1743
NBLOCKSTEST
#define NBLOCKSTEST
Definition: DriveHoughMulti.c:128
uvar_EnableChi2
BOOLEAN uvar_EnableChi2
Definition: DriveHoughMulti.c:102
OpenExtraInfoFiles
int OpenExtraInfoFiles(CHAR *fileMaps, FILE **fp1_ptr, CHAR *filehisto, CHAR *dirname, CHAR *basename, INT4 ind)
Definition: DriveHoughMulti.c:2455
DuplicateBestStuff
void DuplicateBestStuff(LALStatus *status, BestVariables *out, BestVariables *in)
Definition: DriveHoughMulti.c:2016
GetPeakGramFromMultSFTVector
void GetPeakGramFromMultSFTVector(LALStatus *status, HOUGHPeakGramVector *out, MultiSFTVector *in, REAL8 thr)
Loop over SFTs and set a threshold to get peakgrams.
Definition: DriveHoughMulti.c:1683
LALHOUGHDestroyLUTs
void LALHOUGHDestroyLUTs(LALStatus *status, HOUGHptfLUTVector *lut)
Definition: DriveHoughMulti.c:2177
LALHOUGHCreatePHMDVS
void LALHOUGHCreatePHMDVS(LALStatus *status, PHMDVectorSequence *phmdVS, UINT4 length, UINT4 nfSize)
Definition: DriveHoughMulti.c:2211
EARTHEPHEMERIS
#define EARTHEPHEMERIS
Definition: DriveHoughMulti.c:107
LALHOUGHCreateLUTs
void LALHOUGHCreateLUTs(LALStatus *status, HOUGHptfLUTVector *lutV, UINT2 maxNBins, UINT2 maxNBorders, UINT2 ySide)
Definition: DriveHoughMulti.c:2121
LALHOUGHCreateHT
void LALHOUGHCreateHT(LALStatus *status, HOUGHMapTotal *ht, UINT2 xSide, UINT2 ySide)
Definition: DriveHoughMulti.c:2325
DIROUT
#define DIROUT
Definition: DriveHoughMulti.c:112
LALHoughHistogramSignificance
void LALHoughHistogramSignificance(LALStatus *status, UINT8Vector *out, HOUGHMapTotal *in, REAL8 mean, REAL8 sigma, REAL8 minSignificance, REAL8 maxSignificance)
given a total hough map, this function produces a histogram of the number count significance
Definition: DriveHoughMulti.c:1516
LALHOUGHDestroyPHMDs
void LALHOUGHDestroyPHMDs(LALStatus *status, PHMDVectorSequence *phmdVS)
Definition: DriveHoughMulti.c:2303
CreateSkypatchDirs
int CreateSkypatchDirs(CHAR *filestats, CHAR *base, INT4 ind)
Definition: DriveHoughMulti.c:1197
ReadTimeStampsFile
void ReadTimeStampsFile(LALStatus *status, LIGOTimeGPSVector *ts, CHAR *filename)
Definition: DriveHoughMulti.c:1456
PrintHmap2m_file
int PrintHmap2m_file(HOUGHMapTotal *ht, CHAR *fnameOut, INT4 iHmap)
Definition: DriveHoughMulti.c:1300
PrintLogFile
void PrintLogFile(LALStatus *status, CHAR *dir, CHAR *basename, CHAR *skyfile, LALStringVector *linefiles, CHAR *executable)
Definition: DriveHoughMulti.c:1349
TRUE
#define TRUE
Definition: DriveHoughMulti.c:125
FALSE
#define FALSE
Definition: DriveHoughMulti.c:126
HOUGHMAXFILENAMELENGTH
#define HOUGHMAXFILENAMELENGTH
Definition: DriveHoughMulti.c:110
BLOCKSRNGMED
#define BLOCKSRNGMED
Definition: DriveHoughMulti.c:121
uvar_EnableExtraInfo
BOOLEAN uvar_EnableExtraInfo
Definition: DriveHoughMulti.c:102
PrintHistogram
int PrintHistogram(UINT8Vector *hist, CHAR *fnameOut, REAL8 minSignificance, REAL8 maxSignificance)
Definition: DriveHoughMulti.c:1227
ComputeFoft_NM
void ComputeFoft_NM(LALStatus *status, REAL8Vector *foft, HoughTemplate *pulsarTemplate, REAL8Vector *timeDiffV, REAL8Cart3CoorVector *velV)
Definition: DriveHoughMulti.c:2530
LALHOUGHCreatePHMDs
void LALHOUGHCreatePHMDs(LALStatus *status, PHMDVectorSequence *phmdVS, UINT2 maxNBins, UINT2 maxNBorders, UINT2 ySide)
Definition: DriveHoughMulti.c:2252
FBAND
#define FBAND
Definition: DriveHoughMulti.c:118
BASENAMEOUT
#define BASENAMEOUT
Definition: DriveHoughMulti.c:113
SUNEPHEMERIS
#define SUNEPHEMERIS
Definition: DriveHoughMulti.c:108
GetSFTNoiseWeights
void GetSFTNoiseWeights(LALStatus *status, REAL8Vector *out, MultiNoiseWeights *in)
Definition: DriveHoughMulti.c:1638
PrintHmap2file
int PrintHmap2file(HOUGHMapTotal *ht, CHAR *fnameOut, INT4 iHmap)
Definition: DriveHoughMulti.c:1262
LALHOUGHCreateLUTVector
void LALHOUGHCreateLUTVector(LALStatus *status, HOUGHptfLUTVector *lutV, UINT4 length)
helper function for creating LUT vector
Definition: DriveHoughMulti.c:2088
GetSFTVelTime
void GetSFTVelTime(LALStatus *status, REAL8Cart3CoorVector *velV, LIGOTimeGPSVector *timeV, MultiDetectorStateSeries *in)
Definition: DriveHoughMulti.c:1582
SPINDOWNJUMP
#define SPINDOWNJUMP
Definition: DriveHoughMulti.c:130
LALHOUGHCreateFreqIndVector
void LALHOUGHCreateFreqIndVector(LALStatus *status, UINT8FrequencyIndexVector *freqInd, UINT4 length, REAL8 deltaF)
Definition: DriveHoughMulti.c:2366
sort_fstat_toplist
void sort_fstat_toplist(toplist_t *l)
sorts the toplist with an internal sorting function, used before finally writing it
Definition: FstatToplist.c:140
free_fstat_toplist
void free_fstat_toplist(toplist_t **l)
frees the space occupied by the toplist
Definition: FstatToplist.c:118
create_fstat_toplist
int create_fstat_toplist(toplist_t **tl, UINT8 length)
creates a toplist with length elements, returns -1 on error (usually out of memory),...
Definition: FstatToplist.c:111
write_fstat_toplist_to_fp
int write_fstat_toplist_to_fp(toplist_t *tl, FILE *fp, UINT4 *checksum)
Writes the toplist to an (already open) filepointer Returns the number of written charactes sets the ...
Definition: FstatToplist.c:352
insert_into_fstat_toplist
int insert_into_fstat_toplist(toplist_t *tl, FstatOutputEntry elem)
Inserts an element in to the toplist either if there is space left or the element is larger than the ...
Definition: FstatToplist.c:129
lal_errhandler
lal_errhandler_t lal_errhandler
LAL_ERR_EXIT
int LAL_ERR_EXIT(LALStatus *stat, const char *func, const char *file, const int line, volatile const char *id)
LAL_CALL
#define LAL_CALL(function, statusptr)
j
int j
k
int k
LALCheckMemoryLeaks
void LALCheckMemoryLeaks(void)
LALCalloc
#define LALCalloc(m, n)
LALMalloc
#define LALMalloc(n)
LALFree
#define LALFree(p)
lalPulsarVCSInfoList
const LALVCSInfoList lalPulsarVCSInfoList
NULL-terminated list of VCS and build information for LALPulsar and its dependencies
delta
static double double delta
ABORT
#define ABORT(statusptr, code, mesg)
XLAL_CHECK_LAL
#define XLAL_CHECK_LAL(sp, assertion,...)
TRY
#define TRY(func, statusptr)
ATTATCHSTATUSPTR
#define ATTATCHSTATUSPTR(statusptr)
ASSERT
#define ASSERT(assertion, statusptr, code, mesg)
DETATCHSTATUSPTR
#define DETATCHSTATUSPTR(statusptr)
INITSTATUS
#define INITSTATUS(statusptr)
RETURN
#define RETURN(statusptr)
LALUCHAR2HOUGHPeak
void LALUCHAR2HOUGHPeak(LALStatus *status, HOUGHPeakGram *pgOut, UCHARPeakGram *pgIn)
Compress explicit peak gram.
Definition: PeakSelect.c:160
SFTtoUCHARPeakGram
void SFTtoUCHARPeakGram(LALStatus *status, UCHARPeakGram *pg, const SFTtype *sft, REAL8 thr)
Convert a normalized sft into a peakgram by selecting bins in which power exceeds a threshold.
Definition: PeakSelect.c:371
STRING
#define STRING(a)
fscanf
#define fscanf
fprintf
#define fprintf
k0
coeffs k0
XLALDestroyMultiDetectorStateSeries
void XLALDestroyMultiDetectorStateSeries(MultiDetectorStateSeries *mdetStates)
Helper function to get rid of a multi-IFO DetectorStateSeries Note, this is "NULL-robust" in the sens...
Definition: DetectorStates.c:301
XLALGetMultiDetectorStatesFromMultiSFTs
MultiDetectorStateSeries * XLALGetMultiDetectorStatesFromMultiSFTs(const MultiSFTVector *multiSFTs, const EphemerisData *edat, REAL8 tOffset)
Get the 'detector state' (ie detector-tensor, position, velocity, etc) for the given multi-vector of ...
Definition: DetectorStates.c:537
LALStereo2SkyLocation
void LALStereo2SkyLocation(LALStatus *status, REAL8UnitPolarCoor *sourceLocation, UINT2 xPos, UINT2 yPos, HOUGHPatchGrid *patch, HOUGHDemodPar *parDem)
Find source sky location given stereographic coordinates indexes.
Definition: HoughMap.c:405
HoughTT
REAL8 HoughTT
Total Hough Map pixel type.
Definition: HoughMap.h:113
XLALInitBarycenter
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
Definition: LALInitBarycenter.c:206
XLALDestroyEphemerisData
void XLALDestroyEphemerisData(EphemerisData *edat)
Destructor for EphemerisData struct, NULL robust.
Definition: LALInitBarycenter.c:322
XLALDestroyMultiAMCoeffs
void XLALDestroyMultiAMCoeffs(MultiAMCoeffs *multiAMcoef)
Destroy a MultiAMCoeffs structure.
Definition: LALComputeAM.c:469
XLALComputeMultiAMCoeffs
MultiAMCoeffs * XLALComputeMultiAMCoeffs(const MultiDetectorStateSeries *multiDetStates, const MultiNoiseWeights *multiWeights, SkyPosition skypos)
Multi-IFO version of XLALComputeAMCoeffs().
Definition: LALComputeAM.c:379
BOOLEAN
unsigned char BOOLEAN
UCHAR
unsigned char UCHAR
UINT8
uint64_t UINT8
REAL8
double REAL8
XLAL_INIT_DECL
#define XLAL_INIT_DECL(var,...)
INT8
int64_t INT8
UINT2
uint16_t UINT2
CHAR
char CHAR
UINT4
uint32_t UINT4
INT4
int32_t INT4
REAL4
float REAL4
LALHOUGHConstructHMT_W
void LALHOUGHConstructHMT_W(LALStatus *status, HOUGHMapTotal *ht, UINT8FrequencyIndexVector *freqInd, PHMDVectorSequence *phmdVS)
Calculates the total hough map for a given trajectory in the time-frequency plane and a set of partia...
Definition: DriveHough.c:474
LALHOUGHConstructSpacePHMD
void LALHOUGHConstructSpacePHMD(LALStatus *status, PHMDVectorSequence *phmdVS, HOUGHPeakGramVector *pgV, HOUGHptfLUTVector *lutV)
constructs the space of phmd PHMDVectorSequence *phmdVS, given a HOUGHPeakGramVector *pgV and HOUGHpt...
Definition: DriveHough.c:54
LALHOUGHNormalizeWeights
void LALHOUGHNormalizeWeights(LALStatus *status, REAL8Vector *weightV)
Normalizes weight factors so that their sum is N.
Definition: DriveHough.c:668
LALHOUGHWeighSpacePHMD
void LALHOUGHWeighSpacePHMD(LALStatus *status, PHMDVectorSequence *phmdVS, REAL8Vector *weightV)
Adds weight factors for set of partial hough map derivatives – the weights must be calculated outside...
Definition: DriveHough.c:580
LALHOUGHInitializeWeights
void LALHOUGHInitializeWeights(LALStatus *status, REAL8Vector *weightV)
Initializes weight factors to unity.
Definition: DriveHough.c:633
LALHOUGHupdateSpacePHMDup
void LALHOUGHupdateSpacePHMDup(LALStatus *status, PHMDVectorSequence *phmdVS, HOUGHPeakGramVector *pgV, HOUGHptfLUTVector *lutV)
This function updates the space of phmd increasing the frequency phmdVS->fBinMin by one.
Definition: DriveHough.c:130
LALHOUGHConstructHMT
void LALHOUGHConstructHMT(LALStatus *status, HOUGHMapTotal *ht, UINT8FrequencyIndexVector *freqInd, PHMDVectorSequence *phmdVS)
Given PHMDVectorSequence *phmdVS, the space of phmd, and UINT8FrequencyIndexVector *freqInd,...
Definition: DriveHough.c:301
LAL_UINT8_FORMAT
#define LAL_UINT8_FORMAT
LALNDHOUGHParamPLUT
void LALNDHOUGHParamPLUT(LALStatus *status, HOUGHParamPLUT *out, HOUGHSizePar *sizePar, HOUGHDemodPar *par)
Definition: NDParamPLUT.c:72
VTOT
#define VTOT
Total detector velocity/c TO BE CHANGED DEPENDING ON DETECTOR.
Definition: LUT.h:207
PIXERR
#define PIXERR
Maximum `‘error’' (as a fraction of the width of the thinnest annulus) which allows to consider two b...
Definition: LUT.h:195
LINERR
#define LINERR
Maximum `‘error’' (as a fraction of the width of the thinnest annulus) which allows to represent a ci...
Definition: LUT.h:188
LALHOUGHComputeNDSizePar
void LALHOUGHComputeNDSizePar(LALStatus *status, HOUGHSizePar *out, HOUGHResolutionPar *in)
Definition: PatchGrid.c:204
LALHOUGHFillPatchGrid
void LALHOUGHFillPatchGrid(LALStatus *status, HOUGHPatchGrid *out, HOUGHSizePar *par)
Definition: PatchGrid.c:331
LALHOUGHConstructPLUT
void LALHOUGHConstructPLUT(LALStatus *status, HOUGHptfLUT *lut, HOUGHPatchGrid *patch, HOUGHParamPLUT *par)
Definition: ConstructPLUT.c:149
COORType
INT2 COORType
To be changed to {INT2 COORType} if the number of pixels in the x-direction exceeds 255.
Definition: LUT.h:218
LogPrintfVerbatim
void void LogPrintfVerbatim(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
LogPrintf
void LogPrintf(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
LOG_CRITICAL
LOG_CRITICAL
LOG_NORMAL
LOG_NORMAL
XLALNormalizeMultiSFTVect
MultiPSDVector * XLALNormalizeMultiSFTVect(MultiSFTVector *multsft, UINT4 blockSize, const MultiNoiseFloor *assumeSqrtSX)
Function for normalizing a multi vector of SFTs in a multi IFO search and returns the running-median ...
Definition: NormalizeSFTRngMed.c:159
XLALDestroyMultiPSDVector
void XLALDestroyMultiPSDVector(MultiPSDVector *multvect)
Destroy a multi PSD-vector.
Definition: PSDutils.c:102
XLALComputeMultiNoiseWeights
MultiNoiseWeights * XLALComputeMultiNoiseWeights(const MultiPSDVector *rngmed, UINT4 blocksRngMed, UINT4 excludePercentile)
Computes weight factors arising from MultiSFTs with different noise floors.
Definition: PSDutils.c:285
XLALDestroyMultiNoiseWeights
void XLALDestroyMultiNoiseWeights(MultiNoiseWeights *weights)
Destroy a MultiNoiseWeights object.
Definition: PSDutils.c:172
LAL_PMETRIC_NONE
@ LAL_PMETRIC_NONE
Definition: PtoleMetric.h:96
r
static const INT4 r
LALCreateRandomParams
void LALCreateRandomParams(LALStatus *status, RandomParams **params, INT4 seed)
a
static const INT4 a
LALDestroyRandomParams
void LALDestroyRandomParams(LALStatus *status, RandomParams **params)
LALRemoveKnownLinesInMultiSFTVector
void LALRemoveKnownLinesInMultiSFTVector(LALStatus *status, MultiSFTVector *MultiSFTVect, INT4 width, INT4 window, LALStringVector *linefiles, RandomParams *randPar)
top level function to remove lines from a multi sft vector given a list of files containing list of k...
Definition: lib/SFTClean.c:953
XLALDestroySFTCatalog
void XLALDestroySFTCatalog(SFTCatalog *catalog)
Free an 'SFT-catalogue'.
Definition: SFTcatalog.c:329
XLALLoadMultiSFTs
MultiSFTVector * XLALLoadMultiSFTs(const SFTCatalog *inputCatalog, REAL8 fMin, REAL8 fMax)
Function to load a catalog of SFTs from possibly different detectors.
Definition: SFTfileIO.c:416
XLALDestroyMultiSFTVector
void XLALDestroyMultiSFTVector(MultiSFTVector *multvect)
Destroy a multi SFT-vector.
Definition: SFTtypes.c:424
XLALCreateTimestampVector
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:47
XLALSFTdataFind
SFTCatalog * XLALSFTdataFind(const CHAR *file_pattern, const SFTConstraints *constraints)
Find the list of SFTs matching the file_pattern and satisfying the given constraints,...
Definition: SFTcatalog.c:71
XLALDestroyTimestampVector
void XLALDestroyTimestampVector(LIGOTimeGPSVector *vect)
De-allocate a LIGOTimeGPSVector.
Definition: SFTtimestamps.c:69
COORDINATESYSTEM_EQUATORIAL
COORDINATESYSTEM_EQUATORIAL
LALHoughStatistics
void LALHoughStatistics(LALStatus *status, HoughStats *out, HOUGHMapTotal *in)
This function calculates the maximum number count, minimum number count, average and standard deviati...
Definition: lib/Statistics.c:33
XLALUserVarReadAllInput
int XLALUserVarReadAllInput(BOOLEAN *should_exit, int argc, char *argv[], const LALVCSInfoList vcs_list)
XLALDestroyUserVars
void XLALDestroyUserVars(void)
XLALUserVarGetLog
void CHAR * XLALUserVarGetLog(UserVarLogFormat format)
XLALRegisterNamedUvar
#define XLALRegisterNamedUvar(cvar, name, type, option, category,...)
XLALUserVarWasSet
int XLALUserVarWasSet(const void *cvar)
UVAR_LOGFMT_CFGFILE
UVAR_LOGFMT_CFGFILE
XLALDestroyUINT8Vector
void XLALDestroyUINT8Vector(UINT8Vector *vector)
XLALCreateREAL8Vector
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
XLALDestroyREAL8Vector
void XLALDestroyREAL8Vector(REAL8Vector *vector)
XLALCreateUINT8Vector
UINT8Vector * XLALCreateUINT8Vector(UINT4 length)
XLAL_CHECK_MAIN
#define XLAL_CHECK_MAIN(assertion,...)
XLAL_SUCCESS
XLAL_SUCCESS
XLAL_EFUNC
XLAL_EFUNC
XLALGPSSetREAL8
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
XLALGPSDiff
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)
temp
temp
eta
eta
path
path
n
n
out
out
deltaF
int deltaF
chi2
chi2
lalpulsar.git_version.status
string status
Definition: git_version.py:32
lalpulsar_PiecewiseSearch.filename
filename
Definition: lalpulsar_PiecewiseSearch.py:406
lalpulsar_PiecewiseSearch.seed
seed
Definition: lalpulsar_PiecewiseSearch.py:684
lalpulsar_PiecewiseSearchTemplateBank.stats
stats
Definition: lalpulsar_PiecewiseSearchTemplateBank.py:118
lalpulsar_knope_automation_script.p
p
RUN ANALYSIS SCRIPT ###.
Definition: lalpulsar_knope_automation_script.py:709
lalpulsar_knope_result_page.f0
f0
Definition: lalpulsar_knope_result_page.py:2463
lalpulsar_knope_result_page.f1
f1
Definition: lalpulsar_knope_result_page.py:2471
lalpulsar_knope.fp
fp
Definition: lalpulsar_knope.py:137
make_frame_cache.i
int i
Definition: make_frame_cache.py:77
test_computePSD.edat
edat
Definition: test_computePSD.py:65
ts
ts
uvar_skyRegion
CHAR * uvar_skyRegion
Definition: pulsar_crosscorr.c:85
uvar_dAlpha
REAL8 uvar_dAlpha
Definition: pulsar_crosscorr.c:61
uvar_refTime
REAL8 uvar_refTime
Definition: pulsar_crosscorr.c:64
uvar_dDelta
REAL8 uvar_dDelta
Definition: pulsar_crosscorr.c:61
uvar_startTime
REAL8 uvar_startTime
Definition: pulsar_crosscorr.c:57
uvar_f0
REAL8 uvar_f0
Definition: pulsar_crosscorr.c:58
uvar_skyfile
CHAR * uvar_skyfile
Definition: pulsar_crosscorr.c:84
uvar_dirnameOut
CHAR * uvar_dirnameOut
Definition: pulsar_crosscorr.c:83
uvar_blocksRngMed
INT4 uvar_blocksRngMed
Definition: pulsar_crosscorr.c:55
MAXFILENAMELENGTH
#define MAXFILENAMELENGTH
Definition: pulsar_crosscorr.c:93
uvar_endTime
REAL8 uvar_endTime
Definition: pulsar_crosscorr.c:57
PIXELFACTOR
#define PIXELFACTOR
Definition: pulsar_crosscorr.h:113
alpha
double alpha
AMCoeffs::b
REAL4Vector * b
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:65
AMCoeffs::a
REAL4Vector * a
(weighted) per-SFT antenna-pattern function
Definition: LALComputeAM.h:64
BestVariables
struct fo storing all the variables affected by the selection of a subset of SFTs
Definition: DriveHoughColor.h:178
BestVariables::timeDiffV
REAL8Vector * timeDiffV
the vector of time diffs
Definition: DriveHoughColor.h:181
BestVariables::pgV
HOUGHPeakGramVector * pgV
the vector of peakgrams
Definition: DriveHoughColor.h:183
BestVariables::weightsV
REAL8Vector * weightsV
noise and AM weights
Definition: DriveHoughColor.h:180
BestVariables::velV
REAL8Cart3CoorVector * velV
vector of detector velocities
Definition: DriveHoughColor.h:182
BestVariables::length
UINT4 length
the number of SFTs to be selected
Definition: DriveHoughColor.h:179
COMPLEX8FrequencySeries::data
COMPLEX8Sequence * data
COMPLEX8FrequencySeriesVector::data
COMPLEX8FrequencySeries * data
Pointer to the data array.
Definition: PulsarDataTypes.h:75
COMPLEX8FrequencySeriesVector::length
UINT4 length
Number of elements in array.
Definition: PulsarDataTypes.h:74
DetectorState::vDetector
REAL8 vDetector[3]
Cart.
Definition: DetectorStates.h:130
DetectorState::tGPS
LIGOTimeGPS tGPS
GPS timestamps corresponding to this entry.
Definition: DetectorStates.h:128
DetectorStateSeries::data
DetectorState * data
array of DetectorState entries
Definition: DetectorStates.h:147
DetectorStateSeries::length
UINT4 length
total number of entries
Definition: DetectorStates.h:146
DopplerSkyScanInit
initialization-structure passed to InitDopplerSkyScan()
Definition: DopplerScan.h:134
DopplerSkyScanState
this structure reflects the current state of a DopplerSkyScan
Definition: DopplerScan.h:152
EphemerisData
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
Definition: LALBarycenter.h:127
FstatOutputEntry
Type to hold the fields that will be output in unclustered output file
Definition: FstatToplist.h:32
FstatOutputEntry::Alpha
REAL8 Alpha
Skyposition: longitude in equatorial coords, radians.
Definition: FstatToplist.h:35
FstatOutputEntry::Fstat
REAL8 Fstat
value of 2F
Definition: FstatToplist.h:37
FstatOutputEntry::Delta
REAL8 Delta
skyposition: latitude
Definition: FstatToplist.h:36
FstatOutputEntry::Freq
REAL8 Freq
Frequency at maximum (?) of the cluster.
Definition: FstatToplist.h:33
FstatOutputEntry::f1dot
REAL8 f1dot
spindown value f1dot = df/dt
Definition: FstatToplist.h:34
HOUGHBin2Border
This structure stores the border indexes corresponding to one frequency bin plus the corrections to b...
Definition: LUT.h:234
HOUGHBorder
This structure stores the border of a circle clipped on the projected plane.
Definition: LUT.h:221
HOUGHBorder::xPixel
COORType * xPixel
x pixel index to be marked
Definition: LUT.h:226
HOUGHBorder::ySide
UINT2 ySide
length of xPixel
Definition: LUT.h:225
HOUGHDemodPar
Demodulation parameters needed for the Hough transform; all coordinates are assumed to be with respec...
Definition: LUT.h:353
HOUGHDemodPar::positC
REAL8Cart3Coor positC
(x,y,z): Position of the detector
Definition: LUT.h:357
HOUGHDemodPar::deltaF
REAL8 deltaF
Frequency resolution: df=1/TCOH
Definition: LUT.h:354
HOUGHDemodPar::timeDiff
REAL8 timeDiff
: time difference
Definition: LUT.h:358
HOUGHDemodPar::veloC
REAL8Cart3Coor veloC
(x,y,z): Relative detector velocity
Definition: LUT.h:356
HOUGHDemodPar::spin
REAL8Vector spin
Spin down information.
Definition: LUT.h:359
HOUGHDemodPar::skyPatch
REAL8UnitPolarCoor skyPatch
(alpha, delta): position of the center of the patch
Definition: LUT.h:355
HOUGHMapTotal
This structure stores the Hough map.
Definition: HoughMap.h:130
HOUGHMapTotal::ySide
UINT2 ySide
number of physical pixels in the y direction
Definition: HoughMap.h:142
HOUGHMapTotal::xSide
UINT2 xSide
number of physical pixels in the x direction
Definition: HoughMap.h:141
HOUGHMapTotal::map
HoughTT * map
the pixel counts; the number of elements to allocate is ySide*xSide
Definition: HoughMap.h:143
HOUGHMapTotal::spinRes
REAL8Vector spinRes
Refined spin parameters used in the Hough transform.
Definition: HoughMap.h:139
HOUGHMapTotal::f0Bin
INT8 f0Bin
frequency bin for which it has been constructed
Definition: HoughMap.h:131
HOUGHMapTotal::mObsCoh
UINT4 mObsCoh
ratio of observation time and coherent timescale
Definition: HoughMap.h:133
HOUGHMapTotal::deltaF
REAL8 deltaF
frequency resolution
Definition: HoughMap.h:132
HOUGHParamPLUT
Parameters needed to construct the partial look up table.
Definition: LUT.h:333
HOUGHPatchGrid
This structure stores patch-frequency grid information.
Definition: LUT.h:264
HOUGHPatchGrid::ySide
UINT2 ySide
Real number of pixels in the y-direction (in the projected plane).
Definition: LUT.h:275
HOUGHPatchGrid::xCoor
REAL8 * xCoor
Coordinates of the pixel centers.
Definition: LUT.h:271
HOUGHPatchGrid::xSide
UINT2 xSide
Real number of pixels in the x direction (in the projected plane); it should be less than or equal to...
Definition: LUT.h:270
HOUGHPatchGrid::yCoor
REAL8 * yCoor
Coordinates of the pixel centers.
Definition: LUT.h:278
HOUGHPeakGram
This structure stores the `‘peak-gram’'.
Definition: PHMD.h:129
HOUGHPeakGram::peak
INT4 * peak
The peak indices relative to fBinIni, i.e., the zero peak corresponds to fBinIni.
Definition: PHMD.h:135
HOUGHPeakGramVector
This structure contains a vector of peak-grams (for the different time stamps)
Definition: LALHough.h:167
HOUGHPeakGramVector::length
UINT4 length
number of elements
Definition: LALHough.h:168
HOUGHPeakGramVector::pg
HOUGHPeakGram * pg
the Peakgrams
Definition: LALHough.h:169
HOUGHResolutionPar
parameters needed for gridding the patch
Definition: LUT.h:282
HOUGHResolutionPar::deltaF
REAL8 deltaF
Frequency resolution: df=1/TCOH
Definition: LUT.h:284
HOUGHResolutionPar::pixErr
REAL8 pixErr
for validity of LUT as PIXERR
Definition: LUT.h:288
HOUGHResolutionPar::f0Bin
INT8 f0Bin
Frequency bin at which construct the grid.
Definition: LUT.h:283
HOUGHResolutionPar::patchSkySizeY
REAL8 patchSkySizeY
Patch size in radians along y-axis.
Definition: LUT.h:286
HOUGHResolutionPar::patchSkySizeX
REAL8 patchSkySizeX
Patch size in radians along x-axis.
Definition: LUT.h:285
HOUGHResolutionPar::pixelFactor
REAL8 pixelFactor
number of pixel that fit in the thinnest annulus
Definition: LUT.h:287
HOUGHResolutionPar::vTotC
REAL8 vTotC
estimate value of v-total/C as VTOT
Definition: LUT.h:290
HOUGHResolutionPar::linErr
REAL8 linErr
as LINERR circle ->line
Definition: LUT.h:289
HOUGHSizePar
required for constructing patch
Definition: LUT.h:294
HOUGHSizePar::maxNBorders
UINT2 maxNBorders
maximum number of borders affecting the patch; for memory allocation
Definition: LUT.h:302
HOUGHSizePar::maxNBins
UINT2 maxNBins
maximum number of bins affecting the patch; for memory allocation
Definition: LUT.h:301
HOUGHSizePar::ySide
UINT2 ySide
number of pixels in the y direction
Definition: LUT.h:300
HOUGHSizePar::nFreqValid
INT8 nFreqValid
number of frequencies where the LUT is valid
Definition: LUT.h:303
HOUGHSizePar::xSide
UINT2 xSide
number of pixels in the x direction (projected plane)
Definition: LUT.h:299
HOUGHphmd
This structure stores a partial Hough map derivative.
Definition: PHMD.h:141
HOUGHphmd::firstColumn
UCHAR * firstColumn
Number of elements of firstColumn.
Definition: PHMD.h:151
HOUGHphmd::ySide
UINT2 ySide
number of elements of firstColumn
Definition: PHMD.h:150
HOUGHphmd::rightBorderP
HOUGHBorder ** rightBorderP
Pointers to borders.
Definition: PHMD.h:149
HOUGHphmd::maxNBorders
UINT2 maxNBorders
Maximun number of borders of each type (for memory allocation purposes), i.e. length of *leftBorderP ...
Definition: PHMD.h:145
HOUGHphmd::leftBorderP
HOUGHBorder ** leftBorderP
Pointers to borders.
Definition: PHMD.h:148
HOUGHptfLUT
This structure stores the patch-time-frequency look up table.
Definition: LUT.h:246
HOUGHptfLUT::maxNBins
UINT2 maxNBins
Maximum number of bins affecting the patch (for memory allocation purposes)
Definition: LUT.h:256
HOUGHptfLUT::border
HOUGHBorder * border
The annulus borders.
Definition: LUT.h:258
HOUGHptfLUT::bin
HOUGHBin2Border * bin
Bin to border correspondence.
Definition: LUT.h:259
HOUGHptfLUT::maxNBorders
UINT2 maxNBorders
Maximum number of borders affecting the patch (for memory allocation purposes)
Definition: LUT.h:257
HOUGHptfLUTVector
This structure contains a vector of partial look up tables (for the different time stamps)
Definition: LALHough.h:173
HOUGHptfLUTVector::lut
HOUGHptfLUT * lut
the partial Look Up Tables
Definition: LALHough.h:175
HOUGHptfLUTVector::length
UINT4 length
number of elements
Definition: LALHough.h:174
HoughSignificantEventVector::event
HoughSignificantEvent * event
Definition: DriveHoughColor.h:162
HoughStats
Structure for storing statistical information about a Hough map.
Definition: lib/Statistics.h:76
HoughTemplate::spindown
REAL8Vector spindown
Definition: MCInjectHoughMulti.h:97
LIGOTimeGPSVector
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
LIGOTimeGPSVector::data
LIGOTimeGPS * data
array of timestamps
Definition: SFTfileIO.h:193
LIGOTimeGPSVector::length
UINT4 length
number of timestamps
Definition: SFTfileIO.h:192
MultiAMCoeffs
Multi-IFO container for antenna-pattern coefficients and atenna-pattern matrix .
Definition: LALComputeAM.h:137
MultiAMCoeffs::data
AMCoeffs ** data
noise-weighted AM-coeffs , and
Definition: LALComputeAM.h:142
MultiDetectorStateSeries
Multi-IFO time-series of DetectorStates.
Definition: DetectorStates.h:154
MultiDetectorStateSeries::length
UINT4 length
number of detectors
Definition: DetectorStates.h:158
MultiDetectorStateSeries::data
DetectorStateSeries ** data
vector of pointers to DetectorStateSeries
Definition: DetectorStates.h:159
MultiNoiseWeights
One noise-weight (number) per SFT (therefore indexed over IFOs and SFTs.
Definition: PSDutils.h:71
MultiNoiseWeights::length
UINT4 length
number of detectors
Definition: PSDutils.h:75
MultiNoiseWeights::data
REAL8Vector ** data
weights-vector for each detector
Definition: PSDutils.h:76
MultiPSDVector
A collection of PSD vectors – one for each IFO in a multi-IFO search.
Definition: PSDutils.h:62
MultiSFTVector
A collection of SFT vectors – one for each IFO in a multi-IFO search.
Definition: SFTfileIO.h:179
MultiSFTVector::length
UINT4 length
number of ifos
Definition: SFTfileIO.h:183
MultiSFTVector::data
SFTVector ** data
sftvector for each ifo
Definition: SFTfileIO.h:184
PHMDVectorSequence
This structure contains a vector sequence of partial-Hough maps derivatives (for different time stamp...
Definition: LALHough.h:189
PHMDVectorSequence::fBinMin
UINT8 fBinMin
frequency index of smallest intrinsic frequency in circular buffer
Definition: LALHough.h:192
PHMDVectorSequence::nfSize
UINT4 nfSize
number of different frequencies
Definition: LALHough.h:190
PHMDVectorSequence::deltaF
REAL8 deltaF
frequency resolution
Definition: LALHough.h:193
PHMDVectorSequence::phmd
HOUGHphmd * phmd
the partial Hough map derivatives
Definition: LALHough.h:195
PHMDVectorSequence::length
UINT4 length
number of elements for each frequency
Definition: LALHough.h:191
PulsarDopplerParams
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
Definition: PulsarDataTypes.h:141
PulsarDopplerParams::Delta
REAL8 Delta
Sky position: DEC (latitude) in equatorial coords and radians.
Definition: PulsarDataTypes.h:144
PulsarDopplerParams::Alpha
REAL8 Alpha
Sky position: RA (longitude) in equatorial coords and radians.
Definition: PulsarDataTypes.h:143
REAL4Vector::data
REAL4 * data
REAL8Cart3Coor
Three dimensional Cartessian coordinates.
Definition: LUT.h:308
REAL8Cart3Coor::y
REAL8 y
Definition: LUT.h:310
REAL8Cart3Coor::x
REAL8 x
Definition: LUT.h:309
REAL8Cart3Coor::z
REAL8 z
Definition: LUT.h:311
REAL8Cart3CoorVector::data
REAL8Cart3Coor * data
x.y.z
Definition: DriveHoughColor.h:148
REAL8Cart3CoorVector::length
UINT4 length
number of elements
Definition: DriveHoughColor.h:147
REAL8UnitPolarCoor
Polar coordinates of a unitary vector on the sphere.
Definition: LUT.h:327
REAL8UnitPolarCoor::alpha
REAL8 alpha
any value
Definition: LUT.h:328
REAL8UnitPolarCoor::delta
REAL8 delta
In the interval [ ].
Definition: LUT.h:329
REAL8Vector::data
REAL8 * data
SFTCatalog
An "SFT-catalogue": a vector of SFTdescriptors, as returned by XLALSFTdataFind()
Definition: SFTfileIO.h:238
SFTCatalog::data
SFTDescriptor * data
array of data-entries describing matched SFTs
Definition: SFTfileIO.h:243
SFTCatalog::length
UINT4 length
number of SFTs in catalog
Definition: SFTfileIO.h:242
SFTConstraints
'Constraints' for SFT-matching: which detector, within which time-stretch and which timestamps exactl...
Definition: SFTfileIO.h:212
SFTConstraints::detector
CHAR * detector
2-char channel-prefix describing the detector (eg 'H1', 'H2', 'L1', 'G1' etc)
Definition: SFTfileIO.h:213
SFTConstraints::timestamps
LIGOTimeGPSVector * timestamps
list of timestamps
Definition: SFTfileIO.h:216
SFTConstraints::maxStartTime
LIGOTimeGPS * maxStartTime
only include SFTs whose epoch is < maxStartTime
Definition: SFTfileIO.h:215
SFTConstraints::minStartTime
LIGOTimeGPS * minStartTime
only include SFTs whose epoch is >= minStartTime
Definition: SFTfileIO.h:214
SFTDescriptor::header
SFTtype header
SFT-header info.
Definition: SFTfileIO.h:228
SkyPosition::longitude
REAL8 longitude
SkyPosition::latitude
REAL8 latitude
SkyPosition::system
CoordinateSystem system
UCHARPeakGram
Explicit peakgram structure – 1 if power in bin is above threshold and 0 if below.
Definition: PeakSelect.h:114
UCHARPeakGram::data
UCHAR * data
pointer to the data {0,1}
Definition: PeakSelect.h:120
UCHARPeakGram::timeBase
REAL8 timeBase
coherent time baseline used to construct peakgram
Definition: PeakSelect.h:116
UCHARPeakGram::nPeaks
INT4 nPeaks
number of peaks selected in data
Definition: PeakSelect.h:119
UCHARPeakGram::length
INT4 length
number of elements in data
Definition: PeakSelect.h:118
UCHARPeakGram::fminBinIndex
INT4 fminBinIndex
first frequency bin of peakgram
Definition: PeakSelect.h:117
UCHARPeakGramVector::upg
UCHARPeakGram * upg
expanded Peakgrams
Definition: DriveHoughMulti.c:146
UCHARPeakGramVector::length
UINT4 length
number of elements
Definition: DriveHoughMulti.c:145
UINT8FrequencyIndexVector
This structure stores the frequency indexes of the partial-Hough map derivatives at different time st...
Definition: LALHough.h:149
UINT8FrequencyIndexVector::data
UINT8 * data
the frequency indexes
Definition: LALHough.h:152
UINT8FrequencyIndexVector::deltaF
REAL8 deltaF
frequency resolution
Definition: LALHough.h:151
UINT8FrequencyIndexVector::length
UINT4 length
number of elements
Definition: LALHough.h:150
UINT8Vector::data
UINT8 * data
toplist_t
Definition: HeapToplist.h:35
toplist_t::heap
char ** heap
Definition: HeapToplist.h:40
toplist_t::elems
size_t elems
Definition: HeapToplist.h:37
f_min
double f_min
f_max
double f_max