102 ProcessParamsTable *
ppt = NULL, *ppt2 = NULL;
103 ProcessParamsTable *commandLine = runState->
commandLine;
106 CHAR *inputfile = NULL;
108 CHAR *filestr = NULL;
110 CHAR *efile = NULL, *sfile = NULL, *
tfile = NULL;
113 CHAR *tempdets = NULL;
114 CHAR *tempdet = NULL;
117 CHAR *fakestarts = NULL, *fakelengths = NULL, *fakedt = NULL;
118 REAL8 *fstarts = NULL, *flengths = NULL, *fdt = NULL;
121 INT4 numDets = 1,
i = 0, numPsds = 0, numLengths = 0, numStarts = 0;
138 CHAR *parFile = NULL;
141 runState->
data = NULL;
148 struct timeval time1, time2;
152 gettimeofday( &time1, NULL );
163 CHAR *tmpharms = NULL, *tmpharm = NULL, harmval[256];
169 for (
i = 0;
i < ml;
i++ ) {
173 modelFreqFactors->
data[
i] = atof( harmval );
180 fprintf( stderr,
"Must specify --par-file!\n" );
183 parFile =
ppt->value;
191 if (
ppt && !ppt2 ) {
198 fprintf( stderr,
"Error... too many detectors specified. Increase MAXDETS to be greater than %d if necessary.\n",
MAXDETS );
202 for (
i = 0;
i < numDets;
i++ ) {
210 else if ( ppt2 && !
ppt ) {
216 fprintf( stderr,
"Error... too many detectors specified. Increase MAXDETS to be greater than %d if necessary.\n",
MAXDETS );
224 CHAR *psds = NULL, *tmppsds = NULL, *tmppsd = NULL, psdval[256];
232 if ( ( numPsds =
count_csv( psds ) ) != ml * numDets ) {
233 fprintf( stderr,
"Error... for %d harmonics the number of PSDs for fake data must be %d times the number of \
234detectors specified (no. dets = %d)\n", ml, ml, numDets );
238 for (
i = 0;
i < ml * numDets;
i++ ) {
243 fpsds[
i] = atof( psdval );
249 if ( ( tmpstr = strstr( tempdet,
"A" ) ) != NULL ) {
270 for (
i = 0;
i < ml * numDets;
i++ ) {
280 if ( ( tmpstr = strstr( tempdet,
"A" ) ) != NULL ) {
285 }
else if ( !strcmp(
dets[
FACTOR(
i, ml )],
"V1" ) ) {
288 fprintf( stderr,
"Error... trying to use Advanced detector that is not available!\n" );
301 }
else if ( !strcmp(
dets[
FACTOR(
i, ml )],
"V1" ) ) {
303 }
else if ( !strcmp(
dets[
FACTOR(
i, ml )],
"G1" ) ) {
305 }
else if ( !strcmp(
dets[
FACTOR(
i, ml )],
"T1" ) ) {
308 fprintf( stderr,
"Error... trying to use detector that is not available!\n" );
322 CHAR *tmpstarts = NULL, *tmpstart = NULL, startval[256];
323 fakestarts =
ppt->value;
325 if ( ( numStarts =
count_csv( fakestarts ) ) != numDets * ml ) {
326 fprintf( stderr,
"Error... for %d harmonics the number of start times for fake data must be %d times the number \
327of detectors specified (no. dets = %d)\n", ml, ml, numDets );
332 for (
i = 0;
i < ml * numDets;
i++ ) {
336 fstarts[
i] = atof( startval );
339 for (
i = 0;
i < ml * numDets;
i++ ) {
340 fstarts[
i] = 900000000.;
347 CHAR *tmplengths = NULL, *tmplength = NULL, lengthval[256];
348 fakelengths =
ppt->value;
350 if ( ( numLengths =
count_csv( fakelengths ) ) != numDets * ml ) {
351 fprintf( stderr,
"Error... for %d harmonics the number of data lengths for fake data must be %d times the \
352number of detectors specified (no. dets = %d)\n", ml, ml, numDets );
357 for (
i = 0;
i < ml * numDets;
i++ ) {
360 flengths[
i] = atof( lengthval );
363 for (
i = 0;
i < ml * numDets;
i++ ) {
364 flengths[
i] = 86400.;
371 CHAR *tmpdts = NULL, *tmpdt = NULL, dtval[256];
374 if ( ( numDt =
count_csv( fakedt ) ) != ml * numDets ) {
375 fprintf( stderr,
"Error... for %d harmonics the number of sample time steps for fake data must be %d times the \
376number of detectors specified (no. dets =%d)\n", ml, ml, numDets );
382 for (
i = 0;
i < ml * numDets;
i++ ) {
385 fdt[
i] = atof( dtval );
388 for (
i = 0;
i < ml * numDets;
i++ ) {
396 fprintf( stderr,
"Error... --detectors OR --fake-data needs to be set.\n" );
406 UINT4 nstreams = ml * numDets;
411 fprintf( stderr,
"Error... --outfile needs to be set.\n" );
418 inputfile =
ppt->value;
421 if ( ( inputfile == NULL || strlen( inputfile ) == 0 ) && !ppt2 ) {
422 fprintf( stderr,
"Error... an input file or fake data needs to be set.\n" );
431 if (
count != ml * numDets ) {
432 fprintf( stderr,
"Error... for %d harmonics the number of input files given must be %d times the number of \
433detectors specified (no. dets =%d)\n", ml, ml, numDets );
452 INT4 inputsigma = 0, computesigma = 0, gaussianLike = 0;
462 fprintf( stderr,
"Error... --reheterodyne needs frequency as argument.\n" );
463 fprintf( stderr,
"Provide argument or remove flag\n." );
474 vetothresh = atof(
ppt->value );
486 REAL8 startTimeValue = 0., endTimeValue = INFINITY;
494 if ( endTimeValue <= startTimeValue ) {
495 fprintf( stderr,
"Error... start time is before end time.\n" );
499 if ( ( startTimeValue == 0 ) && ( endTimeValue == INFINITY ) ) {
504 char truncationFlags[][25] = {
"--truncate-time",
"--truncate-samples",
"--truncate-fraction" };
507 for ( trunci = 0; trunci <
sizeof( truncationFlags ) /
sizeof( truncationFlags[0] ); trunci++ ) {
513 fprintf( stderr,
"Error... truncation option requires a value.\n" );
520 if ( truncate > 1 ) {
521 fprintf( stderr,
"Error... can only take one truncation flag.\n" );
528 for (
i = 0, prev = NULL, prevmodel = NULL ;
i < ml * numDets ;
i++, prev = ifodata, prevmodel = ifomodel ) {
532 REAL8 dataValsRe = 0., dataValsIm = 0., sigmaVals = 0.;
534 UINT4 j = 0,
k = 0, datalength = 0;
535 ProcessParamsTable *ppte = NULL, *ppts = NULL, *pptt = NULL;
537 CHAR *filebuf = NULL;
551 ifodata->
next = NULL;
556 ifomodel->
next = NULL;
561 memcpy( freqFactorsCopy->
data, modelFreqFactors->
data,
sizeof(
REAL8 )*modelFreqFactors->
length );
570 CHAR *nonGRmodel = NULL;
577 runState->
data = ifodata;
581 prev->
next = ifodata;
582 prevmodel->
next = ifomodel;
624 fprintf( stderr,
"Error... could not convert data into separate lines.\n" );
632 if ( strchr( tlist->
tokens[
k],
'#' ) || strchr( tlist->
tokens[
k],
'%' ) ) {
643 if ( nvals == 3 && gaussianLike ) {
652 int rc = sscanf( tlist->
tokens[
k],
"%lf%lf%lf", ×, &dataValsRe, &dataValsIm );
656 }
else if ( nvals == 4 ) {
657 int rc = sscanf( tlist->
tokens[
k],
"%lf%lf%lf%lf", ×, &dataValsRe, &dataValsIm, &sigmaVals );
662 fprintf( stderr,
"Error... unrecognised number of values in first line of data file %s.\n",
datafile );
666 if ( fabs( dataValsRe ) > vetothresh || fabs( dataValsIm ) > vetothresh ) {
673 for (
UINT4 ucount = 0; ucount <
j; ucount++ ) {
674 if ( times == temptimes->
data[ucount] ) {
686 if ( !isinf( endTimeValue ) || ( startTimeValue != 0. ) ) {
687 if ( times < startTimeValue ) {
694 if ( times > endTimeValue ) {
702 if (
j >= (
UINT4 )endTimeValue ) {
719 temptimes->
data[
j - 1] = times;
722 if ( rehetfreq != 0 ) {
724 REAL8 deltaPhase = 2. *
LAL_PI * rehetfreq * times;
725 COMPLEX16 dataTemp = ( dataValsRe * cos( -deltaPhase ) - dataValsIm * sin( -deltaPhase ) ) +
726 I * ( dataValsRe * sin( -deltaPhase ) + dataValsIm * cos( -deltaPhase ) );
738 fprintf( stderr,
"Error... nothing read in from data file %s.\n",
datafile );
753 if ( endTimeValue <= 0 || endTimeValue > 1 ) {
754 fprintf( stderr,
"Error... truncation fraction must be between 0 and 1" );
758 int truncationIndex = floor( endTimeValue * datalength );
772 for (
k = 0;
k < datalength;
k++ ) {
777 gsl_sort( temptimes->
data, 1, datalength );
780 for (
k = 2;
k < datalength;
k++ ) {
781 if ( temptimes->
data[
i] - temptimes->
data[
i - 1] < sampledt ) {
782 sampledt = temptimes->
data[
i] - temptimes->
data[
i - 1];
795 gsl_rng_set( runState->
GSLrandom, randshufseed );
797 gsl_rng_set( runState->
GSLrandom, prevseed );
806 if ( gaussianLike ) {
810 datalength = flengths[
i] / fdt[
i];
838 psdscale = sqrt( ( fpsds[
i] / 2. ) / ( 2. * fdt[
i] ) );
841 for (
k = 0;
k < datalength;
k++ ) {
864 if ( ppte && ppts ) {
886 if ( fopen( sfile,
"r" ) == NULL || fopen( efile,
"r" ) == NULL ) {
887 fprintf( stderr,
"Error... ephemeris files not, or incorrectly, defined!\n" );
892 IFO_XTRA_DATA( ifomodel )->times->data[0].gpsSeconds,
IFO_XTRA_DATA( ifomodel )->times->data[datalength - 1].gpsSeconds ) ) ) {
893 fprintf( stderr,
"Error... not been able to set ephemeris files!\n" );
935 UINT4 chunkMin, chunkMax;
940 chunkMin = atoi(
ppt->value );
947 chunkMax = atoi(
ppt->value );
983 UINT4 outputchunks = 0;
987 chunkLength =
chop_n_merge( datatmp, chunkMin, chunkMax, outputchunks );
1001 if ( computesigma ) {
1012 if ( gaussianLike ) {
1016 datatmp = datatmp->
next;
1017 modeltmp = modeltmp->
next;
1029 gettimeofday( &time2, NULL );
1033 tottime = (
REAL8 )( ( time2.tv_sec + time2.tv_usec * 1.e-6 ) - ( time1.tv_sec + time1.tv_usec * 1.e-6 ) );
1091 ProcessParamsTable *
ppt = NULL;
1095 if (
ppt == NULL ) {
1096 fprintf( stderr,
"Must specify --par-file!\n" );
1136 CHAR *binarymodel = NULL;
1157 if (
dt > DeltaT ) {
1163 for (
j = 0;
j < freqFactors->
length;
j++ ) {
1164 REAL8Vector *dts = NULL, *bdts = NULL, *glitchphase = NULL;
1167 if ( freqFactors->
length == 2 ) {
1178 if ( binarymodel != NULL ) {
1189 if ( binarymodel != NULL ) {
1195 if ( binarymodel != NULL ) {
1206 if ( bdts != NULL ) {
1209 if ( glitchphase != NULL ) {
1214 ifo_model = ifo_model->
next;
1219 REAL8 df = 1. / ( 2.*DeltaT );
1252 ProcessParamsTable *
ppt = NULL;
1256 UINT4 i = 0,
k = 0, nsamps = 0, nnlive = 0,
n = 0;
1258 CHAR *nlivevals = NULL, *templives = NULL, *templive = NULL;
1260 CHAR *sampfile = NULL;
1261 CHAR *tempsamps = NULL, *tempsamp = NULL;
1271 if (
ppt != NULL ) {
1279 for (
i = 0;
i < nsamps;
i++ ) {
1288 if (
ppt != NULL ) {
1295 if ( nnlive != nsamps ) {
1299 for (
i = 0;
i < nnlive;
i++ ) {
1302 nlive->data[
i] = atoi( templive );
1308 fprintf( stderr,
"Must set the number of live points used in the input nested samples file.\n\n" );
1317 for (
n = 0;
n < nsamps;
n++ ) {
1330 if ( fopen( namefile,
"r" ) == NULL || fopen( sampfilenames->
data[
n],
"r" ) == NULL ) {
1335 fp = fopen( namefile,
"r" );
1343 fp = fopen( sampfilenames->
data[
n],
"r" );
1344 while ( !feof(
fp ) ) {
1348 for (
j = 0;
j < paramNames->
length;
j++ ) {
1364 for (
j = 0;
j < paramNames->
length;
j++ ) {
1388 for (
k = 1;
k <
i;
k++ ) {
1402 item1 = item1->
next;
1406 Nsamps->data[
n] =
i;
1416 if (
ppt != NULL ) {
1417 Ncell = atoi(
ppt->value );
#define __func__
log an I/O error, i.e.
int PulsarCheckParam(const PulsarParameters *pars, const CHAR *name)
Check for the existence of the parameter name in the PulsarParameters structure.
PulsarParameters * XLALReadTEMPOParFile(const CHAR *pulsarAndPath)
Read in the parameters from a TEMPO(2) parameter file into a PulsarParameters structure.
REAL8 PulsarGetREAL8Param(const PulsarParameters *pars, const CHAR *name)
Return a REAL8 parameter.
const CHAR * PulsarGetStringParam(const PulsarParameters *pars, const CHAR *name)
Return a string parameter.
void PulsarFreeParams(PulsarParameters *pars)
Function to free memory from pulsar parameters.
REAL8 PulsarGetREAL8VectorParamIndividual(const PulsarParameters *pars, const CHAR *name)
Return an individual REAL8 value from a REAL8Vector parameter.
char * XLALFileLoad(const char *path)
TimeCorrectionData * XLALInitTimeCorrections(const CHAR *timeCorrectionFile)
An XLAL interface for reading a time correction file containing a table of values for converting betw...
EphemerisData * XLALInitBarycenter(const CHAR *earthEphemerisFile, const CHAR *sunEphemerisFile)
XLAL interface to reading ephemeris files 'earth' and 'sun', and return ephemeris-data in old backwar...
TimeCorrectionType
Enumerated type denoting the time system type to be produced in the solar system barycentring routine...
@ TIMECORRECTION_ORIGINAL
void LALInferenceAddVariable(LALInferenceVariables *vars, const char *name, const void *value, LALInferenceVariableType type, LALInferenceParamVaryType vary)
void LALInferenceRemoveVariable(LALInferenceVariables *vars, const char *name)
LALInferenceParamVaryType
void * LALInferenceGetVariable(const LALInferenceVariables *vars, const char *name)
ProcessParamsTable * LALInferenceGetProcParamVal(ProcessParamsTable *procparams, const char *name)
LALInferenceVariableItem * LALInferenceGetItem(const LALInferenceVariables *vars, const char *name)
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
LALInferenceParamVaryType LALInferenceGetVariableVaryType(LALInferenceVariables *vars, const char *name)
void LALInferenceSetVariable(LALInferenceVariables *vars, const char *name, const void *value)
LALINFERENCE_PARAM_OUTPUT
LALINFERENCE_REAL8Vector_t
LALINFERENCE_UINT4Vector_t
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
void * XLALRealloc(void *p, size_t n)
double XLALSimNoisePSDGEO(double f)
double XLALSimNoisePSDTAMA(double f)
double XLALSimNoisePSDaLIGOZeroDetHighPower(double f)
double XLALSimNoisePSDAdvVirgo(double f)
double XLALSimNoisePSDiLIGOSRD(double f)
double XLALSimNoisePSDVirgo(double f)
char char * XLALStringDuplicate(const char *s)
size_t XLALStringCopy(char *dst, const char *src, size_t size)
int char * XLALStringAppend(char *s, const char *append)
char * XLALStringToken(char **s, const char *delim, int empty)
RandomParams * XLALCreateRandomParams(INT4 seed)
void XLALDestroyRandomParams(RandomParams *params)
int XLALNormalDeviates(REAL4Vector *deviates, RandomParams *params)
const LALDetector * XLALGetSiteInfo(const CHAR *name)
Find the site geometry-information 'LALDetector' for given a detector name (or prefix).
LIGOTimeGPSVector * XLALCreateTimestampVector(UINT4 len)
Allocate a LIGOTimeGPSVector.
COORDINATESYSTEM_EQUATORIAL
LALStringVector * XLALAppendString2Vector(LALStringVector *vect, const CHAR *string)
COMPLEX16TimeSeries * XLALCreateCOMPLEX16TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
REAL8TimeSeries * XLALResizeREAL8TimeSeries(REAL8TimeSeries *series, int first, size_t length)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
COMPLEX16TimeSeries * XLALResizeCOMPLEX16TimeSeries(COMPLEX16TimeSeries *series, int first, size_t length)
const LALUnit lalSecondUnit
REAL8Vector * XLALResizeREAL8Vector(REAL8Vector *vector, UINT4 length)
REAL4Vector * XLALCreateREAL4Vector(UINT4 length)
UINT4Vector * XLALResizeUINT4Vector(UINT4Vector *vector, UINT4 length)
void XLALDestroyREAL4Vector(REAL4Vector *vector)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
UINT4Vector * XLALCreateUINT4Vector(UINT4 length)
#define XLAL_ERROR_VOID(...)
#define XLAL_CHECK_VOID(assertion,...)
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
void add_initial_variables(LALInferenceVariables *ini, PulsarParameters *pars)
Set up all the allowed variables for a known pulsar search This functions sets up all possible variab...
void setup_lookup_tables(LALInferenceRunState *runState, LALSource *source)
Sets the time angle antenna response lookup table.
void create_kdtree_prior(LALInferenceRunState *runState)
Create a k-d tree from prior samples.
void ns_to_posterior(LALInferenceRunState *runState)
Convert an array of nested samples to posterior samples.
REAL8Vector * get_glitch_phase(PulsarParameters *pars, LIGOTimeGPSVector *datatimes, REAL8Vector *dts, REAL8Vector *bdts)
Computes the phase from the glitch model.
REAL8Vector * get_bsb_delay(PulsarParameters *pars, LIGOTimeGPSVector *datatimes, REAL8Vector *dts, EphemerisData *edat)
Computes the delay between a pulsar in a binary system and the barycentre of the system.
REAL8Vector * get_ssb_delay(PulsarParameters *pars, LIGOTimeGPSVector *datatimes, EphemerisData *ephem, TimeCorrectionData *tdat, TimeCorrectionType ttype, LALDetector *detector)
Computes the delay between a GPS time at Earth and the solar system barycentre.
#define IFO_XTRA_DATA(ifo)
void samples_prior(LALInferenceRunState *runState)
Read in an ascii text file of nested samples, convert to posterior samples and create k-d tree.
void setup_from_par_file(LALInferenceRunState *runState)
Reads in the parameters of the pulsar being searched for.
void read_pulsar_data(LALInferenceRunState *runState)
Reads in the input gravitational wave data files, or creates fake data sets.
Header file for the data reading functions for the parameter estimation code for known pulsar searche...
void check_and_add_fixed_variable(LALInferenceVariables *vars, const char *name, void *value, LALInferenceVariableType type)
Add a variable, checking first if it already exists and is of type LALINFERENCE_PARAM_FIXED and if so...
void compute_variance(LALInferenceIFOData *data, LALInferenceIFOModel *model)
Compute the noise variance for each data segment.
UINT4Vector * get_chunk_lengths(LALInferenceIFOModel *ifo, UINT4 chunkMax)
Split the data into segments.
TimeCorrectionType XLALAutoSetEphemerisFiles(CHAR **efile, CHAR **sfile, CHAR **tfile, PulsarParameters *pulsar, INT4 gpsstart, INT4 gpsend)
Automatically set the solar system ephemeris file based on environment variables and data time span.
UINT4Vector * chop_n_merge(LALInferenceIFOData *data, UINT4 chunkMin, UINT4 chunkMax, UINT4 outputchunks)
Chops and remerges data into stationary segments.
INT4 count_csv(CHAR *csvline)
Counts the number of comma separated values in a string.
#define MAXDETS
The maximum number of different detectors allowable.
#define FACTOR(x, y)
Macro that gives the integer number of times that x goes in to y.
#define CHUNKMAX
Default value of the maximum length into which the data can be split.
#define CHUNKMIN
Default value of the minimum length into which the data can be split.
LALDetector detectors[LAL_NUM_DETECTORS]
This structure contains all information about the center-of-mass positions of the Earth and Sun,...
REAL8TimeSeries * varTimeData
struct tagLALInferenceIFOData * next
COMPLEX16TimeSeries * compTimeData
struct tagLALInferenceIFOModel * next
LALInferenceVariables * params
COMPLEX16TimeSeries * compTimeSignal
REAL8 * ifo_loglikelihoods
LALInferenceVariables * params
LALInferenceIFOModel * ifo
LALSimulationDomain domain
ProcessParamsTable * commandLine
LALInferenceVariables * algorithmParams
struct tagLALInferenceIFOData * data
LALInferenceVariables * priorArgs
LALInferenceThreadState * threads
LALInferenceVariables * currentParams
LALInferenceModel * model
struct tagVariableItem * next
LALInferenceParamVaryType vary
SkyPosition equatorialCoords
The PulsarParameters structure to contain a set of pulsar parameters.