Commit 62cf1fc6 authored by Reinhard Prix's avatar Reinhard Prix

Simplify types carrying multi-IFO detector and noise-floor info

   - split MultiDetectorInfo type into MultiLALDetector and MultiNoiseFloor types
     (current usage indicates that it's better to separate these data structures)
   - use fixed-size MultiXXX arrays to simplify memory handling
   - adapted all affected lalsuite codes
   - adjusted octapps to this API change in issue refs #1090
   - refs #1084
Original: 44644a4ff70e5c89ffa6deb94eb33a144546bd5b
parent 35181c58
......@@ -561,20 +561,8 @@ XLALInitCode ( ConfigVariables *cfg, const UserInput_t *uvar )
}
UINT4 numDetectors = uvar->IFOs->length;
MultiLALDetector *multiDet;
if ( (multiDet = XLALCreateMultiLALDetector ( numDetectors )) == NULL ) {
XLALPrintError ("%s: XLALCreateMultiLALDetector(1) failed with errno=%d\n", __func__, xlalErrno );
XLAL_ERROR ( XLAL_EFUNC );
}
LALDetector *site = NULL;
for ( UINT4 X=0; X < numDetectors; X++ ) {
if ( (site = XLALGetSiteInfo ( uvar->IFOs->data[X] )) == NULL ) {
XLALPrintError ("%s: Failed to get site-info for detector '%s'\n", __func__, uvar->IFOs->data[X] );
XLAL_ERROR ( XLAL_EFUNC );
}
multiDet->data[X] = (*site); /* copy! */
XLALFree ( site );
}
MultiLALDetector multiDet;
XLAL_CHECK ( XLALParseMultiLALDetector ( &multiDet, uvar->IFOs ) == XLAL_SUCCESS, XLAL_EFUNC );
/* init timestamps vector covering observation time */
UINT4 numSteps = (UINT4) ceil ( uvar->dataDuration / uvar->TAtom );
......@@ -598,7 +586,7 @@ XLALInitCode ( ConfigVariables *cfg, const UserInput_t *uvar )
}
/* get detector states */
if ( (cfg->multiDetStates = XLALGetMultiDetectorStates ( multiTS, multiDet, edat, 0.5 * uvar->TAtom )) == NULL ) {
if ( (cfg->multiDetStates = XLALGetMultiDetectorStates ( multiTS, &multiDet, edat, 0.5 * uvar->TAtom )) == NULL ) {
XLALPrintError ( "%s: XLALGetMultiDetectorStates() failed.\n", __func__ );
XLAL_ERROR ( XLAL_EFUNC );
}
......@@ -653,7 +641,6 @@ XLALInitCode ( ConfigVariables *cfg, const UserInput_t *uvar )
} /* if ( uvar->computeLV && uvar->LVlX ) */
/* get rid of all temporary memory allocated for this step */
XLALDestroyMultiLALDetector ( multiDet );
XLALDestroyEphemerisData ( edat );
XLALDestroyMultiTimestamps ( multiTS );
multiTS = NULL;
......
......@@ -697,12 +697,9 @@ XLALInitCode ( ConfigVariables *cfg, const UserInput_t *uvar )
XLALPrintError ("%s: Failed to get site-info for detector '%s'\n", __func__, uvar->IFO );
XLAL_ERROR ( XLAL_EFUNC );
}
MultiLALDetector *multiDet;
if ( (multiDet = XLALCreateMultiLALDetector ( 1 )) == NULL ) { /* currently only implemented for single-IFO case */
XLALPrintError ("%s: XLALCreateMultiLALDetector(1) failed with errno=%d\n", __func__, xlalErrno );
XLAL_ERROR ( XLAL_EFUNC );
}
multiDet->data[0] = (*site); /* copy! */
MultiLALDetector multiDet;
multiDet.length = 1;
multiDet.sites[0] = (*site); /* copy! */
XLALFree ( site );
/* init timestamps vector covering observation time */
......@@ -728,13 +725,12 @@ XLALInitCode ( ConfigVariables *cfg, const UserInput_t *uvar )
}
/* get detector states */
if ( (cfg->multiDetStates = XLALGetMultiDetectorStates ( multiTS, multiDet, edat, 0.5 * uvar->TAtom )) == NULL ) {
if ( (cfg->multiDetStates = XLALGetMultiDetectorStates ( multiTS, &multiDet, edat, 0.5 * uvar->TAtom )) == NULL ) {
XLALPrintError ( "%s: XLALGetMultiDetectorStates() failed.\n", __func__ );
XLAL_ERROR ( XLAL_EFUNC );
}
/* get rid of all temporary memory allocated for this step */
XLALDestroyMultiLALDetector ( multiDet );
XLALDestroyEphemerisData ( edat );
XLALDestroyMultiTimestamps ( multiTS );
multiTS = NULL;
......
......@@ -76,9 +76,10 @@
/** configuration-variables derived from user-variables */
typedef struct
{
EphemerisData *edat; /**< ephemeris-data */
MultiLIGOTimeGPSVector *multiTimestamps;/**< a vector of timestamps to generate time-series/SFTs for */
MultiDetectorInfo detInfo; //!< detectors and noise-floors (for Gaussian noise) to generate data for
EphemerisData *edat; /**< ephemeris-data */
MultiLIGOTimeGPSVector *multiTimestamps; /**< a vector of timestamps to generate time-series/SFTs for */
MultiLALDetector multiIFO; //!< detectors to generate data for
MultiNoiseFloor multiNoiseFloor; //!< ... and corresponding noise-floors to generate Gaussian white noise for
SFTCatalog *noiseCatalog; /**< catalog of noise-SFTs */
MultiSFTCatalogView *multiNoiseCatalogView; /**< multi-IFO 'view' of noise-SFT catalogs */
......@@ -186,7 +187,8 @@ main(int argc, char *argv[])
CWMFDataParams DataParams = empty_CWMFDataParams;
DataParams.fMin = uvar.fmin;
DataParams.Band = uvar.Band;
DataParams.detInfo = GV.detInfo;
DataParams.multiIFO = GV.multiIFO;
DataParams.multiNoiseFloor = GV.multiNoiseFloor;
DataParams.multiTimestamps = (*GV.multiTimestamps);
DataParams.randSeed = uvar.randSeed;
DataParams.SFTWindowType = uvar.SFTWindowType;
......@@ -365,7 +367,7 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
XLAL_CHECK ( ! ( have_IFOs && have_noiseSFTs ), XLAL_EINVAL, "Allow only *one* of --IFOs or --noiseSFTs to determine detectors\n");
if ( have_IFOs ) {
XLAL_CHECK ( XLALParseMultiDetectorInfo ( &(cfg->detInfo), uvar->IFOs, uvar->sqrtSX ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALParseMultiLALDetector ( &(cfg->multiIFO), uvar->IFOs ) == XLAL_SUCCESS, XLAL_EFUNC );
}
// TIMESTAMPS: either from --timestampsFiles, --startTime+duration, or --noiseSFTs
......@@ -406,7 +408,7 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
// extract multi-timestamps from the multi-SFT-catalog view
XLAL_CHECK ( (cfg->multiTimestamps = XLALTimestampsFromMultiSFTCatalogView ( cfg->multiNoiseCatalogView )) != NULL, XLAL_EFUNC );
// extract IFOs from multi-SFT catalog
XLAL_CHECK ( XLALMultiDetectorInfoFromMultiSFTCatalogView ( &(cfg->detInfo), cfg->multiNoiseCatalogView ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALMultiLALDetectorFromMultiSFTCatalogView ( &(cfg->multiIFO), cfg->multiNoiseCatalogView ) == XLAL_SUCCESS, XLAL_EFUNC );
} // endif have_noiseSFTs
else if ( have_timestampsFiles )
......@@ -425,12 +427,20 @@ XLALInitMakefakedata ( ConfigVars_t *cfg, UserVariables_t *uvar )
{
LIGOTimeGPS tStart;
XLALGPSSetREAL8 ( &tStart, uvar->startTime );
XLAL_CHECK ( ( cfg->multiTimestamps = XLALMakeMultiTimestamps ( tStart, uvar->duration, uvar->Tsft, uvar->SFToverlap, cfg->detInfo.length )) != NULL, XLAL_EFUNC );
XLAL_CHECK ( ( cfg->multiTimestamps = XLALMakeMultiTimestamps ( tStart, uvar->duration, uvar->Tsft, uvar->SFToverlap, cfg->multiIFO.length )) != NULL, XLAL_EFUNC );
} // endif have_startTime
else {
XLAL_ERROR (XLAL_EFAILED, "Something went wrong with my internal logic ..\n");
}
// check if the user asked for Gaussian white noise to be produced (sqrtSn[X]!=0), otherwise leave noise-floors at 0
if ( uvar->sqrtSX != NULL ) {
XLAL_CHECK ( XLALParseMultiNoiseFloor ( &(cfg->multiNoiseFloor), uvar->sqrtSX ) == XLAL_SUCCESS, XLAL_EFUNC );
} else {
cfg->multiNoiseFloor.length = cfg->multiIFO.length;
// values remain at their default sqrtSn[X] = 0;
}
/* Init ephemerides */
XLAL_CHECK ( (cfg->edat = XLALInitBarycenter ( uvar->ephemEarth, uvar->ephemSun )) != NULL, XLAL_EFUNC );
......
......@@ -317,19 +317,11 @@ XLALInitCode ( ConfigVariables *cfg, const UserVariables_t *uvar, const char *ap
cfg->averageABCD = uvar->averageABCD;
/* convert detector names into site-info */
MultiLALDetector *multiDet;
XLAL_CHECK ( (multiDet = XLALCreateMultiLALDetector ( numDetectors )) != NULL, XLAL_EFUNC, "XLALCreateMultiLALDetector(1) failed." );
LALDetector *site = NULL;
for ( UINT4 X=0; X < numDetectors; X++ ) {
XLAL_CHECK ( (site = XLALGetSiteInfo ( uvar->IFOs->data[X] )) != NULL, XLAL_EFUNC, "Failed to get site-info for detector '%s'.", uvar->IFOs->data[X] );
multiDet->data[X] = (*site); /* copy! */
XLALFree ( site );
}
MultiLALDetector multiDet;
XLAL_CHECK ( XLALParseMultiLALDetector ( &multiDet, uvar->IFOs ) == XLAL_SUCCESS, XLAL_EFUNC );
/* get detector states */
XLAL_CHECK ( (cfg->multiDetStates = XLALGetMultiDetectorStates ( cfg->multiTimestamps, multiDet, cfg->edat, 0.5 * uvar->Tsft )) != NULL, XLAL_EFUNC, "XLALGetDetectorStates() failed." );
XLALDestroyMultiLALDetector ( multiDet );
XLAL_CHECK ( (cfg->multiDetStates = XLALGetMultiDetectorStates ( cfg->multiTimestamps, &multiDet, cfg->edat, 0.5 * uvar->Tsft )) != NULL, XLAL_EFUNC, "XLALGetDetectorStates() failed." );
cfg->multiWeights = NULL; // FIXME: noise weights not supported (yet)
......
......@@ -251,9 +251,10 @@ int main(int argc, char *argv[])
CWMFDataParams DataParams;
DataParams.fMin = round(inputParams->fmin*inputParams->Tcoh - inputParams->dfmax*inputParams->Tcoh - 0.5*(inputParams->blksize-1) - (REAL8)(inputParams->maxbinshift) - 6.0)/inputParams->Tcoh;
DataParams.Band = round(inputParams->fspan*inputParams->Tcoh + 2.0*inputParams->dfmax*inputParams->Tcoh + (inputParams->blksize-1) + (REAL8)(2.0*inputParams->maxbinshift) + 12.0)/inputParams->Tcoh;
DataParams.detInfo.length = 1;
DataParams.detInfo.sites[0] = *(inputParams->det);
DataParams.detInfo.sqrtSn[0] = 0.0;
DataParams.multiIFO.length = 1;
DataParams.multiIFO.sites[0] = *(inputParams->det);
DataParams.multiNoiseFloor.length = 1;
DataParams.multiNoiseFloor.sqrtSn[0] = args_info.avesqrtSh_arg;
DataParams.multiTimestamps = *multiTimestamps;
DataParams.randSeed = args_info.injRandSeed_arg;
DataParams.SFTWindowType = "Hann";
......@@ -272,7 +273,7 @@ int main(int argc, char *argv[])
//If not signal only, create sfts that include noise or extract a band from real data
if (!inputParams->signalOnly) {
if (args_info.gaussNoiseWithSFTgaps_given || args_info.timestampsFile_given || args_info.segmentFile_given || !(args_info.sftDir_given || args_info.sftFile_given)) {
DataParams.detInfo.sqrtSn[0] = args_info.avesqrtSh_arg;
DataParams.multiNoiseFloor.sqrtSn[0] = args_info.avesqrtSh_arg;
XLAL_CHECK( XLALCWMakeFakeMultiData(&sftvector, NULL, NULL, &DataParams, edat) == XLAL_SUCCESS, XLAL_EFUNC );
} else {
SFTConstraints constraints = empty_constraints;
......@@ -296,7 +297,7 @@ int main(int argc, char *argv[])
//If printing the data outputs, then do that here
if ((args_info.printSignalData_given || args_info.printMarginalizedSignalData_given) && args_info.injectionSources_given) {
DataParams.detInfo.sqrtSn[0] = 0.0;
DataParams.multiNoiseFloor.sqrtSn[0] = 0.0;
PulsarParamsVector *oneSignal = NULL;
XLAL_CHECK( (oneSignal = XLALCreatePulsarParamsVector(1)) != NULL, XLAL_EFUNC );
......
......@@ -35,6 +35,7 @@
#include <lal/SFTutils.h>
#include <lal/LogPrintf.h>
#include <lal/DopplerScan.h>
#include <lal/DetectorStates.h>
#include <lal/ExtrapolatePulsarSpins.h>
#include <lal/LALInitBarycenter.h>
#include <lal/NormalizeSFTRngMed.h>
......@@ -106,7 +107,7 @@ int main(int argc, char *argv[]){
MultiPSDVector *multiPSDs = NULL;
MultiNoiseWeights *multiWeights = NULL;
MultiLIGOTimeGPSVector *multiTimes = NULL;
MultiLALDetector *multiDetectors = NULL;
MultiLALDetector multiDetectors;
MultiDetectorStateSeries *multiStates = NULL;
MultiAMCoeffs *multiCoeffs = NULL;
SFTIndexList *sftIndices = NULL;
......@@ -259,13 +260,13 @@ int main(int argc, char *argv[]){
}
/* read the detector information from the SFTs */
if ((multiDetectors = XLALExtractMultiLALDetectorFromSFTs ( inputSFTs )) == NULL){
LogPrintf ( LOG_CRITICAL, "%s: XLALExtractMultiLALDetectorFromSFTs() failed with errno=%d\n", __func__, xlalErrno );
if ( XLALMultiLALDetectorFromMultiSFTs ( &multiDetectors, inputSFTs ) != XLAL_SUCCESS){
LogPrintf ( LOG_CRITICAL, "%s: XLALMultiLALDetectorFromMultiSFTs() failed with errno=%d\n", __func__, xlalErrno );
XLAL_ERROR( XLAL_EFUNC );
}
/* Find the detector state for each SFT */
if ((multiStates = XLALGetMultiDetectorStates ( multiTimes, multiDetectors, config.edat, 0.0 )) == NULL){
if ((multiStates = XLALGetMultiDetectorStates ( multiTimes, &multiDetectors, config.edat, 0.0 )) == NULL){
LogPrintf ( LOG_CRITICAL, "%s: XLALGetMultiDetectorStates() failed with errno=%d\n", __func__, xlalErrno );
XLAL_ERROR( XLAL_EFUNC );
}
......@@ -459,7 +460,6 @@ int main(int argc, char *argv[]){
XLALDestroySFTIndexList( sftIndices );
XLALDestroyMultiAMCoeffs ( multiCoeffs );
XLALDestroyMultiDetectorStateSeries ( multiStates );
XLALDestroyMultiLALDetector( multiDetectors );
XLALDestroyMultiTimestamps ( multiTimes );
XLALDestroyMultiNoiseWeights ( multiWeights );
XLALDestroyMultiPSDVector ( multiPSDs );
......
......@@ -121,7 +121,8 @@ typedef struct
EphemerisData *edat; /**< ephemeris data (from XLALInitBarycenter()) */
LALSegList segmentList; /**< segment list contains start- and end-times of all segments */
PulsarParams signalParams; /**< GW signal parameters: Amplitudes + doppler */
MultiDetectorInfo detInfo; /**< (multi-)detector info */
MultiLALDetector multiIFO; /**< (multi-)detector info */
MultiNoiseFloor multiNoiseFloor; /**< ... and corresponding noise-floors to be used as weights */
DopplerCoordinateSystem coordSys; /**< array of enums describing Doppler-coordinates to compute metric in */
ResultHistory_t *history; /**< history trail leading up to and including this application */
} ConfigVariables;
......@@ -262,7 +263,8 @@ main(int argc, char *argv[])
metricParams.segmentList = config.segmentList;
metricParams.coordSys = config.coordSys;
metricParams.metricType = uvar.metricType;
metricParams.detInfo = config.detInfo;
metricParams.multiIFO = config.multiIFO;
metricParams.multiNoiseFloor = config.multiNoiseFloor;
metricParams.signalParams = config.signalParams;
metricParams.projectCoord = uvar.projection - 1; /* user-input counts from 1, but interally we count 0=1st coord. (-1==no projection) */
metricParams.approxPhase = uvar.approxPhase;
......@@ -472,7 +474,14 @@ XLALInitCode ( ConfigVariables *cfg, const UserVariables_t *uvar, const char *ap
}
/* ----- initialize IFOs and (Multi-)DetectorStateSeries ----- */
XLAL_CHECK ( XLALParseMultiDetectorInfo ( &cfg->detInfo, uvar->IFOs, uvar->sqrtSX ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALParseMultiLALDetector ( &cfg->multiIFO, uvar->IFOs ) == XLAL_SUCCESS, XLAL_EFUNC );
UINT4 numDet = cfg->multiIFO.length;
XLAL_CHECK ( numDet >= 1, XLAL_EINVAL );
if ( uvar->IFOs ) {
XLAL_CHECK ( XLALParseMultiNoiseFloor ( &cfg->multiNoiseFloor, uvar->IFOs ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( numDet == cfg->multiNoiseFloor.length, XLAL_EINVAL );
}
/* ---------- translate coordinate system into internal representation ---------- */
if ( XLALDopplerCoordinateNames2System ( &cfg->coordSys, uvar->coords ) ) {
......@@ -607,17 +616,17 @@ XLALOutputDopplerMetric ( FILE *fp, const DopplerMetric *metric, const ResultHis
fprintf ( fp, "Tspan = %.1f;\n", Tspan );
fprintf ( fp, "Nseg = %d;\n", Nseg );
fprintf ( fp, "detectors = {");
for ( i=0; i < meta->detInfo.length; i ++ )
for ( i=0; i < meta->multiIFO.length; i ++ )
{
if ( i > 0 ) fprintf ( fp, ", ");
fprintf ( fp, "\"%s\"", meta->detInfo.sites[i].frDetector.name );
fprintf ( fp, "\"%s\"", meta->multiIFO.sites[i].frDetector.name );
}
fprintf ( fp, "};\n");
fprintf ( fp, "detectorWeights = [");
for ( i=0; i < meta->detInfo.length; i ++ )
for ( i=0; i < meta->multiNoiseFloor.length; i ++ )
{
if ( i > 0 ) fprintf ( fp, ", ");
fprintf ( fp, "%f", meta->detInfo.detWeights[i] );
fprintf ( fp, "%f", meta->multiNoiseFloor.sqrtSn[i] );
}
fprintf ( fp, "];\n");
......
......@@ -80,9 +80,10 @@ XLALCWMakeFakeMultiData ( MultiSFTVector **multiSFTs, ///< [out] pointer to op
const MultiLIGOTimeGPSVector *multiTimestamps = &(dataParams->multiTimestamps);
// check multi-detector input
XLAL_CHECK ( dataParams->detInfo.length >= 1, XLAL_EINVAL );
UINT4 numDet = dataParams->detInfo.length;
XLAL_CHECK ( dataParams->multiIFO.length >= 1, XLAL_EINVAL );
UINT4 numDet = dataParams->multiIFO.length;
XLAL_CHECK ( multiTimestamps->length == numDet, XLAL_EINVAL, "Inconsistent number of IFOs: detInfo says '%d', multiTimestamps says '%d'\n", numDet, multiTimestamps->length );
XLAL_CHECK ( dataParams->multiNoiseFloor.length == numDet, XLAL_EINVAL );
// check Tsft, consistent over detectors
REAL8 Tsft = multiTimestamps->data[0]->deltaT;
......@@ -111,9 +112,10 @@ XLALCWMakeFakeMultiData ( MultiSFTVector **multiSFTs, ///< [out] pointer to op
{
/* detector params */
CWMFDataParams dataParamsX = (*dataParams); // struct-copy for general settings
dataParamsX.detInfo.length = 1;
dataParamsX.detInfo.sites[0] = dataParams->detInfo.sites[X];
dataParamsX.detInfo.sqrtSn[0] = dataParams->detInfo.sqrtSn[X];
dataParamsX.multiIFO.length = 1;
dataParamsX.multiIFO.sites[0] = dataParams->multiIFO.sites[X];
dataParamsX.multiNoiseFloor.length = 1;
dataParamsX.multiNoiseFloor.sqrtSn[0] = dataParams->multiNoiseFloor.sqrtSn[X];
MultiLIGOTimeGPSVector mTimestamps = empty_MultiLIGOTimeGPSVector;
mTimestamps.length = 1;
mTimestamps.data = &(multiTimestamps->data[X]); // such that pointer mTimestamps.data[0] = multiTimestamps->data[X]
......@@ -164,7 +166,8 @@ XLALCWMakeFakeData ( SFTVector **SFTvect,
XLAL_CHECK ( dataParams != NULL, XLAL_EINVAL );
XLAL_CHECK ( edat != NULL, XLAL_EINVAL );
XLAL_CHECK ( dataParams->detInfo.length ==1, XLAL_EINVAL );
XLAL_CHECK ( dataParams->multiIFO.length == 1, XLAL_EINVAL );
XLAL_CHECK ( dataParams->multiNoiseFloor.length == 1, XLAL_EINVAL );
XLAL_CHECK ( dataParams->multiTimestamps.length == 1, XLAL_EINVAL );
XLAL_CHECK ( dataParams->multiTimestamps.data[0] != NULL, XLAL_EINVAL );
......@@ -174,7 +177,7 @@ XLALCWMakeFakeData ( SFTVector **SFTvect,
REAL8 fSamp = 2.0 * fBand;
const LIGOTimeGPSVector *timestamps = dataParams->multiTimestamps.data[0];
const LALDetector *site = &dataParams->detInfo.sites[0];
const LALDetector *site = &dataParams->multiIFO.sites[0];
REAL8 Tsft = timestamps->deltaT;
// if SFT output requested: need *effective* fMin and Band consistent with SFT bins
......@@ -276,10 +279,10 @@ XLALCWMakeFakeData ( SFTVector **SFTvect,
XLAL_CHECK ( (Tseries_sum = XLALAddREAL4TimeSeries ( Tseries_sum, Tseries_i )) != NULL, XLAL_EFUNC );
XLALDestroyREAL4TimeSeries ( Tseries_i );
}
} // for iInj = 1 ... (numPulsars-1)
} // for iInj < numSources
// add Gaussian noise if requested
REAL8 sqrtSn = dataParams->detInfo.sqrtSn[0];
/* add Gaussian noise if requested */
REAL8 sqrtSn = dataParams->multiNoiseFloor.sqrtSn[0];
if ( sqrtSn > 0)
{
REAL8 noiseSigma = sqrtSn * sqrt ( fBand );
......
......@@ -72,7 +72,8 @@ typedef struct tagCWMFDataParams
{
REAL8 fMin; //!< smallest frequency guaranteed to be generated [returned fMin can be smaller]
REAL8 Band; //!< smallest frequency band guaranteed to be generated [returned Band can be larger]
MultiDetectorInfo detInfo; //!< detectors and noise-floors (for Gaussian noise) to generate data for
MultiLALDetector multiIFO; //!< detectors to generate data for
MultiNoiseFloor multiNoiseFloor; //!< ... and corresponding noise-floors to generate Gaussian white noise for
MultiLIGOTimeGPSVector multiTimestamps; //!< timestamps to generate SFTs for
const char *SFTWindowType; //!< window to apply to the SFT timeseries
REAL8 SFTWindowBeta; //!< 'beta' parameter required for *some* windows [otherwise must be 0]
......
This diff is collapsed.
......@@ -99,28 +99,30 @@ typedef struct tagDetectorArm
typedef DetectorArm Detector3Arms[3]; /**< used to allow functions some type/size checking */
/**
* simple multi-IFO array of detector-information, standard LAL-vector
* array of detectors definitions 'LALDetector'
*
*/
typedef struct tagMultiLALDetector
{
UINT4 length; /**< number of IFOs */
LALDetector *data; /**< array of LALDetector structs */
UINT4 length; //!< number of detectors \f$N\f$
LALDetector sites[PULSAR_MAX_DETECTORS]; //!< array of site information
} MultiLALDetector;
/**
* Struct describing a set of detectors with their PSDs and derived noise-weights
*
*/
typedef struct tagMultiDetectorInfo
//! array of detector-specific 'noise floors' (ie PSD values), assumed constant
//! over the frequency-band of interest
typedef struct tagMultiNoiseFloor
{
UINT4 length; //!< number of detectors \f$N\f$
LALDetector sites[PULSAR_MAX_DETECTORS]; //!< array of site information
REAL8 sqrtSn[PULSAR_MAX_DETECTORS]; //!< per-IFO sqrt{Sn} values, \f$\sqrt{S_X}\f$
REAL8 detWeights[PULSAR_MAX_DETECTORS]; //!< (derived) noise-weights, defined as \f$w_X = \frac{S_X^{-1}}{\mathcal{S}^{-1}}\f$
REAL8 calS; //!< noise normalization constant \f$\mathcal{S}^{-1}= \frac{1}{N}\sum_{X=1}^{N} S_X^{-1}\f$
//!< such that \f$\sum_{X=1}^N w_X = N\f$
} MultiDetectorInfo;
UINT4 length; //!< number of detectors \f$N_{\mathrm{det}}\f$
REAL8 sqrtSn[PULSAR_MAX_DETECTORS]; //!< per-IFO sqrt(PSD) values \f$\sqrt{S_X}\f$, where
//!< \f$S_X^{-1}\equiv\frac{1}{N_{\mathrm{sft}}^X} \sum_{\alpha=0}^{N_{\mathrm{sft}}^X-1} S_{X\alpha}^{-1}\f$
//!< with \f$N_{\mathrm{sft}}^X\f$ the number of SFTs (labeled by \f$\alpha\f$) from detector \f$X\f$
REAL8 sqrtSnTotal; //!< overall 'noise-floor': harmonic mean
//!< \f$\mathcal{S}^{-1}\equiv\frac{1}{N_{\mathrm{sft}}}\sum_{X\alpha}^{N_{\mathrm{sft}}-1} S_{X\alpha}^{-1}\f$
//!< with \f$N_{\mathrm{sft}} = \sum_{X=0}^{N_{\mathrm{det}}-1} N_{\mathrm{sft}}^X\f$ the total number of SFTs.
//!< Note: only equals the harmonic mean over the \f$S_X\f$ if \f$N_{\mathrm{sft}}^X=N_{\mathrm{sft}}/N_{\mathrm{det}}\f$ for all \f$X\f$,
//!< while generally: \f$S^{-1} = \frac{1}{N_{\mathrm{det}}}\sum_X \frac{N_{\mathrm{det}} N_{\mathrm{sft}}^X}{N_{\mathrm{sft}}} S_X^{-1}\f$
} MultiNoiseFloor;
/* ----- Output types for LALGetDetectorStates() */
/**
......@@ -180,13 +182,13 @@ LALGetMultiDetectorStates( LALStatus *,
void LALCreateDetectorStateSeries (LALStatus *, DetectorStateSeries **vect, UINT4 length );
DetectorStateSeries*
XLALGetDetectorStates ( const LIGOTimeGPSVector *timestamps, const LALDetector *detector, const EphemerisData *edat, REAL8 tOffset );
MultiDetectorStateSeries*
XLALGetMultiDetectorStates( const MultiLIGOTimeGPSVector *multiTS, const MultiLALDetector *multiIFO, const EphemerisData *edat, REAL8 tOffset );
DetectorStateSeries* XLALGetDetectorStates ( const LIGOTimeGPSVector *timestamps, const LALDetector *detector, const EphemerisData *edat, REAL8 tOffset );
MultiDetectorStateSeries* XLALGetMultiDetectorStates( const MultiLIGOTimeGPSVector *multiTS, const MultiLALDetector *multiIFO, const EphemerisData *edat, REAL8 tOffset );
int XLALParseMultiDetectorInfo ( MultiDetectorInfo *detInfo, const LALStringVector *detNames, const LALStringVector *sqrtSX );
int XLALMultiDetectorInfoFromMultiSFTCatalogView ( MultiDetectorInfo *multiDetInfo, const MultiSFTCatalogView *multiView );
int XLALParseMultiLALDetector ( MultiLALDetector *multiIFO, const LALStringVector *detNames );
int XLALParseMultiNoiseFloor ( MultiNoiseFloor *multiNoiseFloor, const LALStringVector *sqrtSX );
int XLALMultiLALDetectorFromMultiSFTCatalogView ( MultiLALDetector *multiIFO, const MultiSFTCatalogView *multiView );
int XLALMultiLALDetectorFromMultiSFTs ( MultiLALDetector *multiIFO, const MultiSFTVector *multiSFTs );
int XLALAddSymmTensor3s ( SymmTensor3 *sum, const SymmTensor3 *aT, const SymmTensor3 *bT );
int XLALSubtractSymmTensor3s ( SymmTensor3 *diff, const SymmTensor3 *aT, const SymmTensor3 *bT );
......@@ -196,16 +198,12 @@ int XLALSymmetricTensorProduct3 ( SymmTensor3 *vxw, REAL4 v[3], REAL4 w[3] );
REAL4 XLALContractSymmTensor3s ( const SymmTensor3 *T1, const SymmTensor3 *T2 );
/* creators */
MultiLALDetector *XLALCreateMultiLALDetector ( UINT4 numDetectors );
MultiLALDetector *XLALExtractMultiLALDetectorFromSFTs ( const MultiSFTVector *multiSFTs );
DetectorStateSeries *XLALCreateDetectorStateSeries ( UINT4 length );
/* destructors */
void XLALDestroyDetectorStateSeries ( DetectorStateSeries *detStates );
void LALDestroyDetectorStateSeries(LALStatus *, DetectorStateSeries **vect );
void XLALDestroyMultiDetectorStateSeries ( MultiDetectorStateSeries *mdetStates );
void XLALDestroyMultiLALDetector ( MultiLALDetector *multiIFO );
/* helpers */
......
......@@ -48,6 +48,7 @@
*/
/*---------- SCALING FACTORS ----------*/
#define INIT_MEM(x) memset(&(x), 0, sizeof((x)))
/** shortcuts for integer powers */
#define POW2(a) ( (a) * (a) )
......@@ -235,7 +236,7 @@ double CW_am1_am2_Phi_i_Phi_j ( double tt, void *params );
double CWPhaseDeriv_i ( double tt, void *params );
double XLALAverage_am1_am2_Phi_i_Phi_j ( const intparams_t *params, double *relerr_max );
double CWPhase_cov_Phi_ij ( const MultiDetectorInfo *detInfo, const intparams_t *params, double *relerr_max );
double CWPhase_cov_Phi_ij ( const MultiLALDetector *multiIFO, const MultiNoiseFloor *multiNoiseFloor, const intparams_t *params, double *relerr_max );
UINT4 findHighestGCSpinOrder ( const DopplerCoordinateSystem *coordSys );
......@@ -785,11 +786,26 @@ XLALPtolemaicPosVel ( PosVel3D_t *posvel, /**< [out] instantaneous position and
/**
* Compute a pure phase-deriv covariance \f$[\phi_i, \phi_j] = \langle phi_i phi_j\rangle - \langle phi_i\rangle\langle phi_j\rangle\f$
* which gives a component of the "phase metric"
* which gives a component of the "phase metric".
*
* NOTE: for passing unit noise-weights, set MultiNoiseFloor->length=0 (but multiNoiseFloor==NULL is invalid)
*/
double
CWPhase_cov_Phi_ij ( const MultiDetectorInfo *detInfo, const intparams_t *params, double* relerr_max )
CWPhase_cov_Phi_ij ( const MultiLALDetector *multiIFO, //!< [in] detectors to use
const MultiNoiseFloor *multiNoiseFloor, //!< [in] corresponding noise floors for weights, NULL means unit-weights
const intparams_t *params, //!< [in] integration parameters
double* relerr_max //!< [in] maximal error for integration
)
{
XLAL_CHECK_REAL8 ( multiIFO != NULL, XLAL_EINVAL );
UINT4 numDet = multiIFO->length;
XLAL_CHECK_REAL8 ( numDet > 0, XLAL_EINVAL );
// either no noise-weights given (multiNoiseFloor->length=0) or same number of detectors
XLAL_CHECK_REAL8 ( multiNoiseFloor != NULL, XLAL_EINVAL );
BOOLEAN haveNoiseWeights = (multiNoiseFloor->length > 0);
XLAL_CHECK_REAL8 ( !haveNoiseWeights || (multiNoiseFloor->length == numDet), XLAL_EINVAL );
gsl_function integrand;
intparams_t par = (*params); /* struct-copy, as the 'deriv' field has to be changeable */
......@@ -806,7 +822,6 @@ CWPhase_cov_Phi_ij ( const MultiDetectorInfo *detInfo, const intparams_t *params
REAL8 Tseg = params->Tseg;
UINT4 Nseg = (UINT4) ceil ( params->Tspan / Tseg );
UINT4 n;
REAL8 dT = 1.0 / Nseg;
double epsrel = params->epsrel;
......@@ -827,16 +842,16 @@ CWPhase_cov_Phi_ij ( const MultiDetectorInfo *detInfo, const intparams_t *params
// loop over detectors
REAL8 total_weight = 0;
for (UINT4 det = 0; det < detInfo->length; ++det) {
for (UINT4 X = 0; X < numDet; X ++) {
// set detector for phase integrals
par.site = &detInfo->sites[det];
par.site = &multiIFO->sites[X];
// accumulate detector weights
const REAL8 weight = detInfo->detWeights[det];
const REAL8 weight = haveNoiseWeights ? multiNoiseFloor->sqrtSn[X] : 1.0;
total_weight += weight;
for (n=0; n < Nseg; n ++ ) {
for (UINT4 n=0; n < Nseg; n ++ ) {
REAL8 ti = 1.0 * n * dT;
REAL8 tf = MYMIN( (n+1.0) * dT, 1.0 );
......@@ -885,12 +900,10 @@ CWPhase_cov_Phi_ij ( const MultiDetectorInfo *detInfo, const intparams_t *params
} /* for i < Nseg */
} // for det < detInfo->length
} // for X < numDet
// raise error if no detector weights were given
if (total_weight == 0) {
XLAL_ERROR( XLAL_EDOM, "Detectors weights are all zero!" );
}
XLAL_CHECK_REAL8 (total_weight > 0, XLAL_EDOM, "Detectors noise-floors given but all zero!" );
// normalise by total weight
const REAL8 inv_total_weight = 1.0 / total_weight;
......@@ -1041,7 +1054,7 @@ XLALDopplerPhaseMetric ( const DopplerMetricParams *metricParams, /**< input p
/* g_ij */
intparams.deriv1 = coordSys->coordIDs[i];
intparams.deriv2 = coordSys->coordIDs[j];
REAL8 gg = CWPhase_cov_Phi_ij ( &metricParams->detInfo, &intparams, &err ); /* [Phi_i, Phi_j] */
REAL8 gg = CWPhase_cov_Phi_ij ( &metricParams->multiIFO, &metricParams->multiNoiseFloor, &intparams, &err ); /* [Phi_i, Phi_j] */
maxrelerr = MYMAX ( maxrelerr, err );
if ( xlalErrno ) {
XLALPrintError ("\n%s: Integration of g_ij (i=%d, j=%d) failed. errno = %d\n", __func__, i, j, xlalErrno );
......@@ -1328,6 +1341,8 @@ XLALDopplerFstatMetricCoh ( const DopplerMetricParams *metricParams, /**< inpu
/**
* Function to the compute the FmetricAtoms_t, from which the F-metric and Fisher-matrix can be computed.
*
* NOTE: if MultiNoiseFloor.length=0, unit noise weights are assumed.
*/
FmetricAtoms_t*
XLALComputeAtomsForFmetric ( const DopplerMetricParams *metricParams, /**< input parameters determining the metric calculation */
......@@ -1337,7 +1352,7 @@ XLALComputeAtomsForFmetric ( const DopplerMetricParams *metricParams, /**< inp
FmetricAtoms_t *ret; /* return struct */
intparams_t intparams = empty_intparams;
UINT4 dim, numDet, i=-1, j=-1, X; /* index counters */
UINT4 dim, i=-1, j=-1, X; /* index counters */
REAL8 A, B, C; /* antenna-pattern coefficients (gsl-integrated) */
const DopplerCoordinateSystem *coordSys;
......@@ -1350,10 +1365,15 @@ XLALComputeAtomsForFmetric ( const DopplerMetricParams *metricParams, /**< inp
XLALPrintError ("\n%s: Illegal NULL pointer passed!\n\n", __func__);
XLAL_ERROR_NULL( XLAL_EINVAL );
}
UINT4 numDet = metricParams->multiIFO.length;
XLAL_CHECK_NULL ( numDet >= 1, XLAL_EINVAL );
XLAL_CHECK_NULL ( XLALSegListIsInitialized ( &(metricParams->segmentList) ), XLAL_EINVAL, "Passed un-initialzied segment list 'metricParams->segmentList'\n");
UINT4 Nseg = metricParams->segmentList.length;
XLAL_CHECK_NULL ( Nseg == 1, XLAL_EINVAL, "Segment list must only contain Nseg=1 segments, got Nseg=%d", Nseg );
BOOLEAN haveNoiseWeights = (metricParams->multiNoiseFloor.length > 0);
XLAL_CHECK_NULL ( !haveNoiseWeights || metricParams->multiNoiseFloor.length == numDet, XLAL_EINVAL );
LIGOTimeGPS *startTime = &(metricParams->segmentList.segs[0].start);
LIGOTimeGPS *endTime = &(metricParams->segmentList.segs[0].end);
REAL8 Tspan = XLALGPSDiff( endTime, startTime );
......@@ -1361,7 +1381,6 @@ XLALComputeAtomsForFmetric ( const DopplerMetricParams *metricParams, /**< inp
const LIGOTimeGPS *refTime = &(metricParams->signalParams.Doppler.refTime);
dim = metricParams->coordSys.dim; /* shorthand: number of Doppler dimensions */
numDet = metricParams->detInfo.length;
coordSys = &(metricParams->coordSys);
/* ----- create output structure ---------- */
......@@ -1406,10 +1425,10 @@ XLALComputeAtomsForFmetric ( const DopplerMetricParams *metricParams, /**< inp
A = B = C = 0;
for ( X = 0; X < numDet; X ++ )
{
REAL8 weight = metricParams->detInfo.detWeights[X];
REAL8 weight = haveNoiseWeights ? metricParams->multiNoiseFloor.sqrtSn[X] : 1.0;
sum_weights += weight;
REAL8 av, relerr;
intparams.site = &(metricParams->detInfo.sites[X]);
intparams.site = &(metricParams->multiIFO.sites[X]);
intparams.deriv1 = DOPPLERCOORD_NONE;
intparams.deriv2 = DOPPLERCOORD_NONE;
......@@ -1465,9 +1484,9 @@ XLALComputeAtomsForFmetric ( const DopplerMetricParams *metricParams, /**< inp
for ( X = 0; X < numDet; X ++ )
{
REAL8 weight = metricParams->detInfo.detWeights[X];
REAL8 weight = haveNoiseWeights ? metricParams->multiNoiseFloor.sqrtSn[X] : 1.0;
REAL8 av, relerr;
intparams.site = &(metricParams->detInfo.sites[X]);
intparams.site = &(metricParams->multiIFO.sites[X]);
/* ------------------------------ */