Commit bf0f9bb6 authored by Reinhard Prix's avatar Reinhard Prix
Browse files

ComputeFstat setup API: packed all optional function arguments into FstatOptionalArgs

- less cluttered function interface, easier to pass 'sane defaults'
- provide global initializer 'FstatOptionalArgsDefaults'
- refs #1936
Original: 9abb52cdd4f73e43778dcfc2e3a740f4da763241
parent a99131d8
......@@ -1374,13 +1374,16 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
PulsarParamsVector *injectSources = NULL;
MultiNoiseFloor *injectSqrtSX = NULL;
FstatExtraParams XLAL_INIT_DECL(extraParams);
extraParams.Dterms = uvar->Dterms;
extraParams.SSBprec = uvar->SSBprecision;
cfg->Fstat_in = XLALCreateFstatInput( catalog, fCoverMin, fCoverMax,
injectSources, injectSqrtSX, assumeSqrtSX, uvar->RngMedWindow,
cfg->ephemeris, cfg->FstatMethod, &extraParams );
XLAL_CHECK ( cfg->Fstat_in != NULL, XLAL_EFUNC, "XLALCreateFstatInput() failed with errno=%d", xlalErrno);
FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
optionalArgs.Dterms = uvar->Dterms;
optionalArgs.SSBprec = uvar->SSBprecision;
optionalArgs.runningMedianWindow = uvar->RngMedWindow;
optionalArgs.injectSources = injectSources;
optionalArgs.injectSqrtSX = injectSqrtSX;
optionalArgs.assumeSqrtSX = assumeSqrtSX;
optionalArgs.FstatMethod = cfg->FstatMethod;
XLAL_CHECK ( (cfg->Fstat_in = XLALCreateFstatInput( catalog, fCoverMin, fCoverMax, cfg->ephemeris, &optionalArgs )) != NULL, XLAL_EFUNC );
XLALDestroySFTCatalog(catalog);
cfg->Fstat_what = FSTATQ_2F; // always calculate multi-detector 2F
......
......@@ -2025,38 +2025,35 @@ void SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */
}
/* loop over segments and read sfts */
in->nSFTs = 0;
for (UINT4 X = 0; X < numDetectors; X++) {
in->NSegmentsInvX[X] = 0;
}
FstatExtraParams XLAL_INIT_DECL(extraParams);
extraParams.SSBprec = in->SSBprec;
extraParams.Dterms = in->Dterms;
FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
optionalArgs.SSBprec = in->SSBprec;
optionalArgs.Dterms = in->Dterms;
optionalArgs.runningMedianWindow = in->blocksRngMed;
optionalArgs.FstatMethod = in->Fmethod;
/* loop over segments and read sfts */
for (k = 0; k < in->nStacks; k++) {
/* if single-only flag is given, assume a PSD with sqrt(S) = 1.0 */
MultiNoiseFloor s_assumeSqrtSX, *assumeSqrtSX;
MultiNoiseFloor s_assumeSqrtSX;
if ( in->SignalOnly ) {
const SFTCatalog *catalog_k = &(catalogSeq.data[k]);
s_assumeSqrtSX.length = XLALCountIFOsInCatalog ( catalog_k );
for (UINT4 X = 0; X < s_assumeSqrtSX.length; ++X) {
s_assumeSqrtSX.sqrtSn[X] = 1.0;
}
assumeSqrtSX = &s_assumeSqrtSX;
optionalArgs.assumeSqrtSX = &s_assumeSqrtSX;
} else {
assumeSqrtSX = NULL;
optionalArgs.assumeSqrtSX = NULL;
}
PulsarParamsVector *injectSources = NULL;
MultiNoiseFloor *injectSqrtSX = NULL;
/* ----- create Fstat input data struct ----- */
(*p_Fstat_in_vec)->data[k] = XLALCreateFstatInput ( &catalogSeq.data[k], freqmin, freqmax,
injectSources, injectSqrtSX, assumeSqrtSX, in->blocksRngMed,
in->edat, in->Fmethod, &extraParams );
(*p_Fstat_in_vec)->data[k] = XLALCreateFstatInput ( &catalogSeq.data[k], freqmin, freqmax, in->edat, &optionalArgs );
if ( (*p_Fstat_in_vec)->data[k] == NULL ) {
XLALPrintError("%s: XLALCreateFstatInput() failed with errno=%d", __func__, xlalErrno);
ABORT ( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
......
......@@ -1426,17 +1426,17 @@ void SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */
nSFTs = 0;
#endif
FstatExtraParams XLAL_INIT_DECL(extraParams);
extraParams.SSBprec = in->SSBprec;
extraParams.Dterms = in->Dterms;
FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
optionalArgs.SSBprec = in->SSBprec;
optionalArgs.Dterms = in->Dterms;
optionalArgs.runningMedianWindow = in->blocksRngMed;
optionalArgs.FstatMethod = FMETHOD_DEMOD_BEST;
/* loop over stacks and read sfts */
for (k = 0; k < in->nStacks; k++) {
/* create Fstat input data struct for Fstat-computation */
(*p_Fstat_in_vec)->data[k] = XLALCreateFstatInput( &catalogSeq.data[k], fMin, fMax,
NULL, NULL, NULL, in->blocksRngMed,
in->edat, FMETHOD_DEMOD_BEST, &extraParams );
(*p_Fstat_in_vec)->data[k] = XLALCreateFstatInput( &catalogSeq.data[k], fMin, fMax, in->edat, &optionalArgs );
if ( (*p_Fstat_in_vec)->data[k] == NULL ) {
XLALPrintError("%s: XLALCreateFstatInput() failed with errno=%d", __func__, xlalErrno);
ABORT ( status, HIERARCHICALSEARCH_EXLAL, HIERARCHICALSEARCH_MSGEXLAL );
......
......@@ -79,14 +79,17 @@ struct tagFstatInput {
// ---------- Include F-statistic method implementations ---------- //
#if CFS_HAVE_SSE
const int FMETHOD_DEMOD_BEST = FMETHOD_DEMOD_SSE;
#define DEF_FMETHOD_DEMOD_BEST FMETHOD_DEMOD_SSE
#elif CFS_HAVE_ALTIVEC
const int FMETHOD_DEMOD_BEST = FMETHOD_DEMOD_ALTIVEC;
#define DEF_FMETHOD_DEMOD_BEST FMETHOD_DEMOD_ALTIVEC
#else
const int FMETHOD_DEMOD_BEST = FMETHOD_DEMOD_OPTC;
#define DEF_FMETHOD_DEMOD_BEST FMETHOD_DEMOD_OPTC
#endif
#define DEF_FMETHOD_RESAMP_BEST FMETHOD_RESAMP_GENERIC
const int FMETHOD_RESAMP_BEST = FMETHOD_RESAMP_GENERIC;
// these are for exporting
const int FMETHOD_DEMOD_BEST = DEF_FMETHOD_DEMOD_BEST;
const int FMETHOD_RESAMP_BEST = DEF_FMETHOD_RESAMP_BEST;
static const struct {
const char *const name;
......@@ -100,6 +103,22 @@ static const struct {
[FMETHOD_RESAMP_GENERIC] = {"ResampGeneric", 1 }
} ;
// ----- global variables ----------
/// global initializer for setting FstatOptionalArgs to default values
const FstatOptionalArgs FstatOptionalArgsDefaults = {
.randSeed = 0,
.SSBprec = SSBPREC_RELATIVISTICOPT,
.Dterms = 8,
.runningMedianWindow = 50,
.FstatMethod = DEF_FMETHOD_DEMOD_BEST,
.injectSources = NULL,
.injectSqrtSX = NULL,
.assumeSqrtSX = NULL
};
#include "ComputeFstat_Demod.c"
#include "ComputeFstat_Resamp.c"
......@@ -240,26 +259,11 @@ XLALDestroyMultiFstatAtomVector ( MultiFstatAtomVector *multiAtoms ///< [in] #M
///
FstatInput *
XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFTs to either load from files, or generate in memory.
///< The \c locator field of each ::SFTDescriptor must be \c !NULL for SFT loading,
///< and \c NULL for SFT generation.
///< The \c locator field of each ::SFTDescriptor must be \c !=NULL for SFT loading, and \c ==NULL for SFT generation.
const REAL8 minCoverFreq, ///< [in] Minimum instantaneous frequency which will be covered over the SFT time span.
const REAL8 maxCoverFreq, ///< [in] Maximum instantaneous frequency which will be covered over the SFT time span.
const PulsarParamsVector *injectSources, ///< [in] Optional vector of parameters of CW signals to simulate and inject.
const MultiNoiseFloor *injectSqrtSX, ///< [in] Optional array of single-sided PSD values governing fake Gaussian noise generation.
///< If supplied, then fake Gaussian noise with the given PSD values will be added to the SFTs.
const MultiNoiseFloor *assumeSqrtSX, ///< [in] Optional array of single-sided PSD values governing the calculation of SFT noise weights.
///< If supplied, then SFT noise weights are calculated from constant spectra with the given PSD
///< values; otherwise, SFT noise weights are calculated from PSDs computed from a running median
///< of the SFTs themselves.
const UINT4 runningMedianWindow, ///< [in] If SFT noise weights are calculated from the SFTs, the running median window length to use.
const EphemerisData *ephemerides, ///< [in] Ephemerides for the time-span of the SFTs.
const FstatMethodType FstatMethod, ///< [in] Method to use for computing the \f$\mathcal{F}\f$-statistic.
const FstatExtraParams *extraParams ///< [in] Minor tuning or method-specific parameters.
const FstatOptionalArgs *optionalArgs ///< [in] Optional 'advanced-level' and method-specific extra arguments; NULL: use defaults from FstatOptionalArgsDefaults.
)
{
// Check catalog
......@@ -271,19 +275,28 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
"All 'locator' fields of SFTDescriptors in 'SFTcatalog' must be either NULL or !NULL." );
}
// Check remaining parameters
// Check remaining required parameters
XLAL_CHECK_NULL ( isfinite(minCoverFreq) && minCoverFreq > 0, XLAL_EINVAL );
XLAL_CHECK_NULL ( isfinite(maxCoverFreq) && maxCoverFreq > 0, XLAL_EINVAL );
XLAL_CHECK_NULL ( maxCoverFreq > minCoverFreq, XLAL_EINVAL );
XLAL_CHECK_NULL ( injectSources == NULL || injectSources->length > 0, XLAL_EINVAL );
XLAL_CHECK_NULL ( injectSources == NULL || injectSources->data != NULL, XLAL_EINVAL );
XLAL_CHECK_NULL ( ephemerides != NULL, XLAL_EINVAL );
XLAL_CHECK_NULL ( extraParams != NULL, XLAL_EINVAL );
XLAL_CHECK_NULL ( extraParams->SSBprec < SSBPREC_LAST, XLAL_EINVAL );
// handle optional arguments, if given
const FstatOptionalArgs *optArgs;
if ( optionalArgs != NULL ) {
optArgs = optionalArgs;
} else {
optArgs = &FstatOptionalArgsDefaults;
}
// check optional arguments sanity
XLAL_CHECK_NULL ( (optArgs->injectSources == NULL) || ((optArgs->injectSources->length > 0) && (optArgs->injectSources->data != NULL)), XLAL_EINVAL );
XLAL_CHECK_NULL ( (optArgs->injectSqrtSX == NULL) || (optArgs->injectSqrtSX->length > 0), XLAL_EINVAL );
XLAL_CHECK_NULL ( (optArgs->assumeSqrtSX == NULL) || (optArgs->assumeSqrtSX->length > 0), XLAL_EINVAL );
XLAL_CHECK_NULL ( optArgs->SSBprec < SSBPREC_LAST, XLAL_EINVAL );
// Determine whether to load and/or generate SFTs
const BOOLEAN loadSFTs = (SFTcatalog->data[0].locator != NULL);
const BOOLEAN generateSFTs = (injectSources != NULL) || (injectSqrtSX != NULL);
const BOOLEAN generateSFTs = (optArgs->injectSources != NULL) || (optArgs->injectSqrtSX != NULL);
XLAL_CHECK_NULL ( loadSFTs || generateSFTs, XLAL_EINVAL, "Can neither load nor generate SFTs with given parameters" );
// Create top-level input data struct
......@@ -295,20 +308,20 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
FstatInput_Common *const common = input->common; // handy shortcut
// create method-specific input data
if ( XLALFstatMethodClassIsDemod ( FstatMethod ) )
if ( XLALFstatMethodClassIsDemod ( optArgs->FstatMethod ) )
{
XLAL_CHECK_NULL ( (input->demod = XLALCalloc ( 1, sizeof(FstatInput_Demod) )) != NULL, XLAL_ENOMEM );
input->demod->Dterms = extraParams->Dterms;
input->demod->Dterms = optArgs->Dterms;
}
else if ( XLALFstatMethodClassIsResamp ( FstatMethod ) )
else if ( XLALFstatMethodClassIsResamp ( optArgs->FstatMethod ) )
{
XLAL_CHECK_NULL ( (input->resamp = XLALCalloc(1, sizeof(FstatInput_Resamp))) != NULL, XLAL_ENOMEM );
}
else
{
XLAL_ERROR_NULL ( XLAL_EINVAL, "Received invalid Fstat method enum '%d'\n", FstatMethod );
XLAL_ERROR_NULL ( XLAL_EINVAL, "Received invalid Fstat method enum '%d'\n", optArgs->FstatMethod );
}
common->FstatMethod = FstatMethod;
common->FstatMethod = optArgs->FstatMethod;
// Determine the time baseline of an SFT
const REAL8 Tsft = 1.0 / SFTcatalog->data[0].header.deltaF;
......@@ -337,7 +350,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
XLAL_CHECK_NULL ( extraBinsMethod >= 0, XLAL_EFAILED );
// Add number of extra frequency bins required by running median
int extraBinsFull = extraBinsMethod + runningMedianWindow/2 + 1; // NOTE: running-median window needed irrespective of assumeSqrtSX!
int extraBinsFull = extraBinsMethod + optArgs->runningMedianWindow/2 + 1; // NOTE: running-median window needed irrespective of assumeSqrtSX!
// Extend frequency range by number of extra bins times SFT bin width
const REAL8 extraFreqMethod = extraBinsMethod / Tsft;
......@@ -377,8 +390,8 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
} // end: if !loadSFTs
// Check length of multi-noise floor arrays
XLAL_CHECK_NULL ( injectSqrtSX == NULL || injectSqrtSX->length == common->detectors.length, XLAL_EINVAL );
XLAL_CHECK_NULL ( assumeSqrtSX == NULL || assumeSqrtSX->length == common->detectors.length, XLAL_EINVAL );
XLAL_CHECK_NULL ( (optArgs->injectSqrtSX == NULL) || (optArgs->injectSqrtSX->length == common->detectors.length), XLAL_EINVAL );
XLAL_CHECK_NULL ( (optArgs->assumeSqrtSX == NULL) || (optArgs->assumeSqrtSX->length == common->detectors.length), XLAL_EINVAL );
// Generate SFTs with injections and noise, if required
if (generateSFTs)
......@@ -389,18 +402,18 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
MFDparams.Band = maxFreqFull - minFreqFull;
MFDparams.multiIFO = common->detectors;
MFDparams.multiTimestamps = *(common->multiTimestamps);
MFDparams.randSeed = extraParams->randSeed;
MFDparams.randSeed = optArgs->randSeed;
// Set noise floors if sqrtSX is given; otherwise noise floors are zero
if ( injectSqrtSX != NULL ) {
MFDparams.multiNoiseFloor = (*injectSqrtSX);
if ( optArgs->injectSqrtSX != NULL ) {
MFDparams.multiNoiseFloor = *(optArgs->injectSqrtSX);
} else {
MFDparams.multiNoiseFloor.length = common->detectors.length;
}
// Generate SFTs with injections
MultiSFTVector *fakeMultiSFTs = NULL;
XLAL_CHECK_NULL ( XLALCWMakeFakeMultiData ( &fakeMultiSFTs, NULL, injectSources, &MFDparams, ephemerides ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK_NULL ( XLALCWMakeFakeMultiData ( &fakeMultiSFTs, NULL, optArgs->injectSources, &MFDparams, ephemerides ) == XLAL_SUCCESS, XLAL_EFUNC );
// If SFTs were loaded, add generated SFTs to then, otherwise just used generated SFTs
if (multiSFTs != NULL) {
......@@ -419,10 +432,10 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
// Normalise SFTs using either running median or assumed PSDs
MultiPSDVector *runningMedian;
XLAL_CHECK_NULL ( (runningMedian = XLALNormalizeMultiSFTVect ( multiSFTs, runningMedianWindow, assumeSqrtSX )) != NULL, XLAL_EFUNC );
XLAL_CHECK_NULL ( (runningMedian = XLALNormalizeMultiSFTVect ( multiSFTs, optArgs->runningMedianWindow, optArgs->assumeSqrtSX )) != NULL, XLAL_EFUNC );
// Calculate SFT noise weights from PSD
XLAL_CHECK_NULL ( (common->multiNoiseWeights = XLALComputeMultiNoiseWeights ( runningMedian, runningMedianWindow, 0 )) != NULL, XLAL_EFUNC );
XLAL_CHECK_NULL ( (common->multiNoiseWeights = XLALComputeMultiNoiseWeights ( runningMedian, optArgs->runningMedianWindow, 0 )) != NULL, XLAL_EFUNC );
// at this point we're done with running-median noise estimation and can 'trim' the SFTs back to
// the width actually required by the Fstat-methods *methods*.
......@@ -436,7 +449,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Catalog of SFT
// Save ephemerides and SSB precision
common->ephemerides = ephemerides;
common->SSBprec = extraParams->SSBprec;
common->SSBprec = optArgs->SSBprec;
// Call the appropriate method function to setup their input data structures
// - The method input data structures are expected to take ownership of the
......
......@@ -145,14 +145,22 @@ extern const int FMETHOD_DEMOD_BEST;
extern const int FMETHOD_RESAMP_BEST;
///
/// Struct for collecting 'lower-level' tuning and method-specific parameters to be passed to the
/// \f$\mathcal{F}\f$-statistic setup function.
/// Struct of optional 'advanced level' and (potentially method-specific) arguments to be passed to the
/// \f$\mathcal{F}\f$-statistic setup function XLALCreateFstatInput().
///
typedef struct tagFstatExtraParams {
UINT4 randSeed; ///< Random-number seed value used for fake Gaussian noise generation.
SSBprecision SSBprec; ///< Barycentric transformation precision.
UINT4 Dterms; ///< Number of Dirichlet kernel terms, used by some \a Demod methods; see #FstatMethodType
} FstatExtraParams;
typedef struct tagFstatOptionalArgs {
UINT4 randSeed; ///< Random-number seed value used in case of fake Gaussian noise generation (\c injectSqrtSX)
SSBprecision SSBprec; ///< Barycentric transformation precision.
UINT4 Dterms; ///< Number of Dirichlet kernel terms, used by some \a Demod methods; see #FstatMethodType.
UINT4 runningMedianWindow; ///< If SFT noise weights are calculated from the SFTs, the running median window length to use.
FstatMethodType FstatMethod; ///< Method to use for computing the \f$\mathcal{F}\f$-statistic.
PulsarParamsVector *injectSources; ///< Vector of parameters of CW signals to simulate and inject.
MultiNoiseFloor *injectSqrtSX; ///< Single-sided PSD values for fake Gaussian noise to be added to SFT data.
MultiNoiseFloor *assumeSqrtSX; ///< Single-sided PSD values to be used for computing SFT noise weights instead of from a running median of the SFTs themselves.
} FstatOptionalArgs;
/// global initializer for setting FstatOptionalArgs to default values
extern const FstatOptionalArgs FstatOptionalArgsDefaults;
///
/// An \f$\mathcal{F}\f$-statistic 'atom', i.e. the elementary per-SFT quantities required to compute the
......@@ -289,11 +297,7 @@ MultiFstatAtomVector* XLALCreateMultiFstatAtomVector ( const UINT4 length );
void XLALDestroyMultiFstatAtomVector ( MultiFstatAtomVector *atoms );
FstatInput *
XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, const REAL8 minCoverFreq, const REAL8 maxCoverFreq,
const PulsarParamsVector *injectSources, const MultiNoiseFloor *injectSqrtSX,
const MultiNoiseFloor *assumeSqrtSX, const UINT4 runningMedianWindow,
const EphemerisData *ephemerides,
const FstatMethodType FstatMethod, const FstatExtraParams *extraParams );
XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, const REAL8 minCoverFreq, const REAL8 maxCoverFreq, const EphemerisData *ephemerides, const FstatOptionalArgs *optionalArgs );
const MultiLALDetector* XLALGetFstatInputDetectors ( const FstatInput* input );
const MultiLIGOTimeGPSVector* XLALGetFstatInputTimestamps ( const FstatInput* input );
......
......@@ -120,9 +120,6 @@ main ( int argc, char *argv[] )
} // for l < numSegments
LIGOTimeGPS endTime = endTime_l;
UINT4 rngMed = 50;
UINT4 Dterms = 8;
PulsarSpinRange XLAL_INIT_DECL(spinRange);
LIGOTimeGPS refTime = { startTime.gpsSeconds - 2.3 * uvar->Tseg, 0 };
spinRange.refTime = refTime;
......@@ -134,7 +131,7 @@ main ( int argc, char *argv[] )
REAL8 minCoverFreq, maxCoverFreq;
XLAL_CHECK ( XLALCWSignalCoveringBand ( &minCoverFreq, &maxCoverFreq, &startTime, &endTime, &spinRange, asini, Period, ecc ) == XLAL_SUCCESS, XLAL_EFUNC );
UINT4 numBinsSFT = ceil ( (maxCoverFreq - minCoverFreq) * uvar->Tsft + 2 * Dterms );
UINT4 numBinsSFT = ceil ( (maxCoverFreq - minCoverFreq) * uvar->Tsft + 2 * 8 );
UINT4 numSFTsPerSeg = catalogs[0]->length;
REAL8 memSFTs = uvar->numSegments * numSFTsPerSeg * ( sizeof(SFTtype) + numBinsSFT * sizeof(COMPLEX8)) / 1e6;
......@@ -147,23 +144,22 @@ main ( int argc, char *argv[] )
Doppler.ecc = ecc;
Doppler.asini = asini;
// ----- setup extra Fstat method params
FstatExtraParams XLAL_INIT_DECL(extraParams);
extraParams.randSeed = 1;
extraParams.SSBprec = SSBPREC_RELATIVISTICOPT;
extraParams.Dterms = Dterms; // constant value that works for all Demod methods
// ----- setup optional Fstat arguments
FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
MultiNoiseFloor XLAL_INIT_DECL(injectSqrtSX);
injectSqrtSX.length = numDetectors;
for ( UINT4 X=0; X < numDetectors; X ++ ) {
injectSqrtSX.sqrtSn[X] = 1;
}
optionalArgs.injectSqrtSX = &injectSqrtSX;
optionalArgs.FstatMethod = FstatMethod;
FstatInputVector *inputs;
XLAL_CHECK ( (inputs = XLALCreateFstatInputVector ( uvar->numSegments )) != NULL, XLAL_EFUNC );
for ( INT4 l = 0; l < uvar->numSegments; l ++ )
{
XLAL_CHECK ( (inputs->data[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, NULL, &injectSqrtSX, NULL, rngMed, ephem, FstatMethod, &extraParams )) != NULL, XLAL_EFUNC );
XLAL_CHECK ( (inputs->data[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, ephem, &optionalArgs )) != NULL, XLAL_EFUNC );
}
// ----- compute Fstatistics over segments
......
......@@ -144,11 +144,11 @@ main ( int argc, char *argv[] )
REAL8 minCoverFreq, maxCoverFreq;
XLAL_CHECK ( XLALCWSignalCoveringBand ( &minCoverFreq, &maxCoverFreq, &startTime, &endTime, &spinRange, asini, Period, ecc ) == XLAL_SUCCESS, XLAL_EFUNC );
// ----- setup extra Fstat method params
FstatExtraParams XLAL_INIT_DECL(extraParams);
extraParams.randSeed = 1;
extraParams.SSBprec = SSBPREC_RELATIVISTICOPT;
extraParams.Dterms = 8; // constant value that works for all Demod methods
// ----- setup optional Fstat arguments
FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
optionalArgs.injectSources = injectSources;
optionalArgs.injectSqrtSX = &injectSqrtSX;
optionalArgs.assumeSqrtSX = &assumeSqrtSX;
// ----- prepare input data with injection for all available methods
FstatInput *input[FMETHOD_END];
......@@ -159,8 +159,8 @@ main ( int argc, char *argv[] )
if ( !XLALFstatMethodIsAvailable(iMethod) ) {
continue;
}
input[iMethod] = XLALCreateFstatInput ( catalog, minCoverFreq, maxCoverFreq, injectSources, &injectSqrtSX, &assumeSqrtSX, 50, ephem, iMethod, &extraParams );
XLAL_CHECK ( input[iMethod] != NULL, XLAL_EFUNC );
optionalArgs.FstatMethod = iMethod;
XLAL_CHECK ( (input[iMethod] = XLALCreateFstatInput ( catalog, minCoverFreq, maxCoverFreq, ephem, &optionalArgs )) != NULL, XLAL_EFUNC );
}
FstatQuantities whatToCompute = (FSTATQ_2F | FSTATQ_FAFB);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment