Commit bb23f7e4 authored by Karl Wette's avatar Karl Wette

Merge branch 'fstat-check-max-sft-length' into 'master'

ComputeFstat: check that SFT length is within allowed maximum

Closes CW/lalsuite#41

See merge request lscsoft/lalsuite!492
parents d59255b6 50c73dc5
Pipeline #36177 passed with stages
in 155 minutes and 50 seconds
......@@ -68,9 +68,9 @@ pcc_CL="--startTime=$startTime --endTime=$endTime --sftLocation='$SFTdir/*.sft'
## ---------- Run MFDv4 ----------
cmdline="$mfd_path $mfd_CL1";
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${mfd_code} ... "
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "FAILED:"
echo $cmdline
exit 1
......@@ -79,9 +79,9 @@ else
fi
cmdline="$mfd_path $mfd_CL2";
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${mfd_code} ... "
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "FAILED:"
echo $cmdline
exit 1
......@@ -91,7 +91,7 @@ fi
## ---------- Run PulsarCrossCorr_v2 ----------
cmdline="$pcc_path $pcc_CL"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${pcc_code} ... "
if ! tmp=`eval $cmdline 2> /dev/null`; then
echo "FAILED:"
......
......@@ -1422,6 +1422,10 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
XLALFree ( detector );
}
/* maximum ranges of binary orbit parameters */
const REAL8 binaryMaxAsini = MYMAX( uvar->orbitasini, uvar->orbitasini + uvar->orbitasiniBand );
const REAL8 binaryMinPeriod = MYMIN( uvar->orbitPeriod, uvar->orbitPeriod + uvar->orbitPeriodBand );
const REAL8 binaryMaxEcc = MYMAX( uvar->orbitEcc, uvar->orbitEcc + uvar->orbitEccBand );
{ /* ----- What frequency-band do we need to read from the SFTs?
* propagate spin-range from refTime to startTime and endTime of observation
......@@ -1443,13 +1447,18 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
memcpy ( cfg->searchRegion.fkdotBand, spinRangeRef.fkdotBand, sizeof(cfg->searchRegion.fkdotBand) );
// Compute covering frequency range, accounting for Doppler modulation due to the Earth and any binary companion */
const REAL8 binaryMaxAsini = MYMAX( uvar->orbitasini, uvar->orbitasini + uvar->orbitasiniBand );
const REAL8 binaryMinPeriod = MYMIN( uvar->orbitPeriod, uvar->orbitPeriod + uvar->orbitPeriodBand );
const REAL8 binaryMaxEcc = MYMAX( uvar->orbitEcc, uvar->orbitEcc + uvar->orbitEccBand );
XLALCWSignalCoveringBand( &fCoverMin, &fCoverMax, &cfg->startTime, &endTime, &spinRangeRef, binaryMaxAsini, binaryMinPeriod, binaryMaxEcc );
} /* extrapolate spin-range */
/* check that SFT length is within allowed maximum */
{
/* use fCoverMax here to work with loading grid from file (--gridType=6) */
const REAL8 Tsft_max = XLALFstatMaximumSFTLength( fCoverMax, binaryMaxAsini, binaryMinPeriod );
XLAL_CHECK ( !XLALIsREAL8FailNaN( Tsft_max ), XLAL_EINVAL );
XLAL_CHECK ( cfg->Tsft < Tsft_max, XLAL_EINVAL, "Length of input SFTs (%g s) must be less than %g s for CW signal with frequency = %g, binary asini = %g, period = %g", cfg->Tsft, Tsft_max, fCoverMax, binaryMaxAsini, binaryMinPeriod );
}
/* if single-only flag is given, assume a PSD with sqrt(S) = 1.0 */
MultiNoiseFloor s_assumeSqrtSX, *assumeSqrtSX;
if ( uvar->SignalOnly ) {
......
......@@ -16,13 +16,6 @@ export LC_ALL=C
builddir="./";
injectdir="../Injections/"
## ----- user-controlled level of debug-output detail
if [ -n "$DEBUG" ]; then
debug=${DEBUG}
else
debug=0 ## default=quiet
fi
## ----- allow user-control of hotloop variant to use
if [ -n "$FSTAT_METHOD" ]; then
FstatMethod="--FstatMethod=${FSTAT_METHOD}"
......@@ -116,7 +109,7 @@ fi
cmdline="$mfd_code $mfd_CL --randSeed=1"
echo $cmdline;
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$mfd_code' ..."
exit 1
fi
......@@ -145,7 +138,7 @@ fi
cmdline="$cfsv2_code $cfs_CL --outputFstat=${outfile_Fstat} --TwoFthreshold=0 --FreqBand=$cfs_FreqBand"
echo $cmdline;
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1;
fi
......
......@@ -36,9 +36,7 @@ f1dot=1e-09
IFO=H1
## fixed binary parameters
##orbitasini=2.94
orbitasini=4.94036957341473
##orbitArgp=5.2
orbitasini=2.94
orbitArgp=2.88
orbitTp=22665.001234567GPS
......@@ -73,9 +71,9 @@ outfileCFS="CFS.dat";
mfd_CL1="--refTime=${refTime} --Alpha=$Alpha --Delta=$Delta --Tsft=$Tsft --startTime=$startTime --duration=$duration --h0=$h0 --cosi=$cosi --psi=$psi --phi0=$phi0 --f1dot=$f1dot --orbitasini=$orbitasini --orbitPeriod=$orbitPeriod --orbitEcc=$orbitEcc --orbitArgp=$orbitArgp --orbitTp=$orbitTp"
mfd_CL="${mfd_CL1} --fmin=$mfd_fmin --Band=$mfd_FreqBand --Freq=$Freq --outSFTbname=$SFTdir/catsft.sft --outSingleSFT --IFO=$IFO"
cmdline="$mfd_code $mfd_CL"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running $cmdline ... "
if ! eval "$cmdline &> /dev/null"; then
if ! eval "$cmdline"; then
echo "lalapps_Makefakedata_v4 failed"
exit 1;
fi
......@@ -83,7 +81,7 @@ echo "done."
## ---------- predict F-stat value 'PFS'
cmdline="${pfs_code} --Alpha=$Alpha --Delta=$Delta --Freq=$Freq --h0=$h0 --cosi=$cosi --psi=$psi --phi0=$phi0 --DataFiles='$SFTdir/*.sft' --PureSignal --assumeSqrtSX=1"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running $cmdline ... "
resPFS0=`$cmdline`
resPFS=`echo $resPFS0 | awk '{printf "%g", $1}'`
......@@ -92,9 +90,9 @@ echo "done: 2F_PFS = $resPFS"
## ---------- targeted CFSv2 search
cfs_CL=" --refTime=${refTime} --Alpha=$Alpha --Delta=$Delta --Freq=$Freq --f1dot=$f1dot --orbitasini=$orbitasini --orbitPeriod=$orbitPeriod --orbitEcc=$orbitEcc --orbitArgp=$orbitArgp --orbitTp=$orbitTp --DataFiles='$SFTdir/*.sft' --assumeSqrtSX=1"
cmdline="$cfs_code $cfs_CL --outputLoudest=$outfileCFS"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running $cmdline ... "
if ! eval "$cmdline &> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1
fi
......
......@@ -16,12 +16,6 @@ export LC_ALL=C
builddir="./";
injectdir="../Injections/"
## ----- user-controlled level of debug-output detail
debug=0 ## default=quiet
if [ -n "$DEBUG" ]; then
debug=${DEBUG}
fi
##---------- names of codes and input/output files
mfd_code="${injectdir}lalapps_Makefakedata_v5"
cmp_code="${builddir}lalapps_compareFstats"
......@@ -117,7 +111,7 @@ dataSpec="--Tsft=$Tsft --startTime=$startTime --duration=$duration --sqrtSX=${sq
mfd_CL_comp="${dataSpec} --injectionSources='${injectionSources}' --outSingleSFT --outSFTdir=${SFTdir2} --randSeed=1 --IFOs=H1,L1"
cmdline="$mfd_code $mfd_CL_comp "
echo $cmdline;
if ! eval "$cmdline &> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$mfd_code' ..."
exit 1
fi
......@@ -126,7 +120,7 @@ echo
mfd_CL="${dataSpec} --outSingleSFT --outSFTdir=${SFTdir} --randSeed=1 --IFOs=H1,L1"
cmdline="$mfd_code $mfd_CL "
echo $cmdline;
if ! eval "$cmdline &> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$mfd_code' ..."
exit 1
fi
......@@ -142,7 +136,7 @@ echo "----------------------------------------------------------------------"
echo
cmdline="$cfs_code $cfs_CL --FstatMethod=DemodBest --DataFiles='$SFTdir2/*.sft' --outputFstat=$outfile_Comp --outputTiming=$timefile_Comp --outputLoudest=${loudest_Comp} "
echo $cmdline;
if ! eval "$cmdline &> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1
fi
......@@ -154,7 +148,7 @@ echo "----------------------------------------------------------------------"
echo
cmdline="$cfs_code $cfs_CL --injectionSources='${injectionSources}' --FstatMethod=DemodBest --DataFiles='$SFTdir/*.sft' --outputFstat=$outfile_Demod --outputTiming=$timefile_Demod --outputLoudest=${loudest_Demod} "
echo $cmdline;
if ! eval "$cmdline &> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1
fi
......@@ -166,7 +160,7 @@ echo "----------------------------------------------------------------------"
echo
cmdline="$cfs_code $cfs_CL --injectionSources='${injectionSources}' --FstatMethod=ResampBest --DataFiles='$SFTdir/*.sft' --outputFstat=$outfile_Resamp --outputTiming=$timefile_Resamp --outputLoudest=${loudest_Resamp}"
echo $cmdline;
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1;
fi
......
......@@ -16,13 +16,6 @@ export LC_ALL=C
builddir="./";
injectdir="../Injections/"
## ----- user-controlled level of debug-output detail
if [ -n "$DEBUG" ]; then
debug=${DEBUG}
else
debug=0 ## default=quiet
fi
## ----- allow user-control of hotloop variant to use
if [ -n "$FSTAT_METHOD" ]; then
FstatMethod="--FstatMethod=${FSTAT_METHOD}"
......@@ -139,7 +132,7 @@ fi
cmdline="$mfd_code $mfd_CL --randSeed=1"
echo $cmdline;
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$mfd_code' ..."
exit 1
fi
......@@ -162,7 +155,7 @@ cfs_CL_run1="${cfs_CL_base} ${cfs_searchBand} --outputFstat=${outfile_Fstat1} --
cmdline="$cfsv2_code $cfs_CL_run1"
echo $cmdline;
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfsv2_code' ..."
exit 1;
fi
......@@ -260,7 +253,7 @@ t0Band=$(echo $Tdata $taumin | awk '{printf "%d", $1-$2 }');
cfs_CL_run2="${cfs_CL_base} ${cfs_searchBand} --outputFstat=${outfile_Fstat2} --outputLoudest=${outfile_Loudest2} --outputTransientStats=${outfile_transient2} --transient-WindowType=rect --transient-t0Epoch=$t0min --transient-t0Band=$t0Band --transient-dt0=$Tsft --transient-tau=$taumin --transient-tauBand=$tauBand --transient-dtau=$Tsft"
cmdline="$cfsv2_code $cfs_CL_run2"
echo $cmdline;
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfsv2_code' ..."
exit 1;
fi
......@@ -398,7 +391,7 @@ echo
cfs_CL_run3="${cfs_CL_base} --Freq=$Freq --FreqBand=0 --dFreq=0 --f1dot=$f1dot --f1dotBand=0 --df1dot=0 --outputTransientStatsAll=${outfile_transientMap3} --transient-WindowType=rect --transient-t0Epoch=$t0min --transient-t0Band=$t0Band --transient-dt0=$Tsft --transient-tau=$taumin --transient-tauBand=$tauBand --transient-dtau=$Tsft"
cmdline="$cfsv2_code $cfs_CL_run3"
echo $cmdline;
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfsv2_code' ..."
exit 1;
fi
......
......@@ -82,9 +82,9 @@ mfd_CL="${saf_CL} --fmin=$mfd_fmin --Band=$mfd_FreqBand --Freq=$Freq --outSFTbna
## ---------- Run MFDv4 ----------
cmdline="$mfd_path $mfd_CL";
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${mfd_code} ... "
if ! eval "$cmdline 2> /dev/null"; then
if ! eval "$cmdline"; then
echo "FAILED:"
echo $cmdline
exit 1
......@@ -94,7 +94,7 @@ fi
## ---------- Run SemiAnalyticF ----------
cmdline="$saf_path $saf_CL --sqrtSh=$noiseSqrtSh"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${saf_code} ... "
if ! tmp=`eval $cmdline 2> /dev/null`; then
echo "FAILED:"
......@@ -110,7 +110,7 @@ pfs_CL_common=" --Alpha=$Alpha --Delta=$Delta --cosi=$cosi --h0=$h0 --psi=$psi"
outfile_pfs1="__tmp_PFS1.dat";
pfs_CL="${pfs_CL_common} --DataFiles=\"${SFTdir}/*.sft\" --Freq=$Freq --outputFstat=$outfile_pfs1"
cmdline="$pfs_path $pfs_CL"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${pfs_code}{NoiseWeights} ... "
if ! tmp=`eval ${cmdline}`; then
echo "FAILED:"
......@@ -125,7 +125,7 @@ resPFS1=`echo $tmp | awk '{printf "%g", $1}'`
outfile_pfs0="__tmp_PFS0.dat";
pfs_CL="${pfs_CL_common} --DataFiles=\"${SFTdir}/*.sft\" --outputFstat=$outfile_pfs0 --assumeSqrtSX=${noiseSqrtSh}"
cmdline="$pfs_path $pfs_CL"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${pfs_code}{assumeSqrtSX} ... "
if ! tmp=`eval $cmdline`; then
echo "FAILED:"
......@@ -140,7 +140,7 @@ resPFS0=`echo $tmp | awk '{printf "%g", $1}'`
outfile_pfs2="__tmp_PFS2.dat";
pfs_CL="${pfs_CL_common} --timestampsFiles=${timestamps} --Tsft=$Tsft --outputFstat=$outfile_pfs2 --assumeSqrtSX=${noiseSqrtSh} --IFOs=$IFO"
cmdline="$pfs_path $pfs_CL"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${pfs_code}{timestamps,assumeSqrtSX} ... "
if ! tmp=`eval $cmdline`; then
echo "FAILED:"
......@@ -155,7 +155,7 @@ resPFS2=`echo $tmp | awk '{printf "%g", $1}'`
outfile_pfs3="__tmp_PFS3.dat";
pfs_CL="${pfs_CL_common} --minStartTime=$startTime --maxStartTime=$endTime --Tsft=$Tsft --outputFstat=$outfile_pfs3 --assumeSqrtSX=${noiseSqrtSh} --IFOs=$IFO"
cmdline="$pfs_path $pfs_CL"
if [ "$DEBUG" ]; then echo $cmdline; fi
echo $cmdline
echo -n "Running ${pfs_code}{minStartTime+maxStartTime,assumeSqrtSX} ... "
if ! tmp=`eval $cmdline`; then
echo "FAILED:"
......
......@@ -170,11 +170,6 @@ while [ $iFreq -le $numFreqBands ]; do
if [ $iFreq -eq 2 ]; then
cmdline="$cmdline --injectionSources=\"{Freq=$Freq; f1dot=$f1dot; f2dot=$f2dot; Alpha=$Alpha; Delta=$Delta; psi=$psi; phi0=$phi0; h0=$h0; cosi=$cosi; refTime=$refTime}\""
fi
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo "$cmdline";
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$mfd_code' ..."
......@@ -212,11 +207,6 @@ if [ ! -r "$outfile_cfs" ]; then
# ----- get multi-IFO + single-IFO F-stat values
cmdline="$cfs_CL --DataFiles='$SFTfiles'"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo "$cmdline"
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$cfs_code' ..."
......@@ -277,11 +267,6 @@ outfile_GCT_RS="${testDir}/GCT_RS.dat"
timingsfile_RS="${testDir}/timing_RS.dat"
cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS' --outputTiming='$timingsfile_RS' ${BSGL_flags}"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo "$cmdline"
if ! eval "$cmdline"; then
echo "Error.. something failed when running '$gct_code' ..."
......@@ -307,11 +292,6 @@ outfile_GCT_DM="${testDir}/GCT_DM.dat"
timingsfile_DM="${testDir}/timing_DM.dat"
cmdline="$gct_code $gct_CL_common --FstatMethod=DemodBest --fnameout='$outfile_GCT_DM' --outputTiming='$timingsfile_DM' ${BSGL_flags}"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......@@ -339,11 +319,6 @@ outfile_GCT_DM_BSGL="${testDir}/GCT_DM_BSGL.dat"
timingsfile_DM_BSGL="${testDir}/timing_DM_BSGL.dat"
cmdline="$gct_code $gct_CL_common --FstatMethod=DemodBest ${BSGL_flags} --SortToplist=2 --fnameout='$outfile_GCT_DM_BSGL' --outputTiming='$timingsfile_DM_BSGL'"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......@@ -368,11 +343,6 @@ outfile_GCT_DM_DUAL="${testDir}/GCT_DM_DUAL.dat"
timingsfile_DM_DUAL="${testDir}/timing_DM_DUAL.dat"
cmdline="$gct_code $gct_CL_common --FstatMethod=DemodBest --SortToplist=3 ${BSGL_flags} --fnameout='$outfile_GCT_DM_DUAL' --outputTiming='$timingsfile_DM_DUAL'"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......@@ -408,11 +378,6 @@ export LAL_FSTAT_FFT_PLAN_TIMEOUT=30
cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS_triple' --outputTiming='$timingsfile_RS_triple' ${BSGL_flags} --getMaxFperSeg --loudestSegOutput --SortToplist=6"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......@@ -424,11 +389,6 @@ fi
# second variant of triple toplists, with 2F sorted first toplist
cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS_triple2' --outputTiming='$timingsfile_RS_triple2' ${BSGL_flags} --getMaxFperSeg --loudestSegOutput --SortToplist=7"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......@@ -444,11 +404,6 @@ fi
outfile_GCT_RS_triple="${testDir}/GCT_RS_triple_0.dat"
cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS_triple' ${BSGL_flags} --getMaxFperSeg --loudestSegOutput --SortToplist=0"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......@@ -461,11 +416,6 @@ fi
outfile_GCT_RS_triple="${testDir}/GCT_RS_triple_1.dat"
cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS_triple' ${BSGL_flags} --getMaxFperSeg --loudestSegOutput --SortToplist=2"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......@@ -476,11 +426,6 @@ fi
outfile_GCT_RS_triple="${testDir}/GCT_RS_triple_2.dat"
cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS_triple' ${BSGL_flags} --getMaxFperSeg --loudestSegOutput --SortToplist=4"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......@@ -491,11 +436,6 @@ fi
outfile_GCT_RS_triple="${testDir}/GCT_RS_triple_3.dat"
cmdline="$gct_code $gct_CL_common --FstatMethod=ResampGeneric --fnameout='$outfile_GCT_RS_triple' ${BSGL_flags} --getMaxFperSeg --loudestSegOutput --SortToplist=5"
if [ -n "$DEBUG" ]; then
cmdline="$cmdline"
else
cmdline="$cmdline &> /dev/null"
fi
echo $cmdline
if ! eval "$cmdline"; then
......
......@@ -328,3 +328,18 @@ archivePrefix = "arXiv",
editor = {Meshkov, S.},
url = {http://arxiv.org/abs/gr-qc/9912029}
}
@article{LeaciPrix2015,
title = {{Directed searches for continuous gravitational waves from binary systems: Parameter-space metrics and optimal Scorpius X-1 sensitivity}},
author = {Leaci, Paola and Prix, Reinhard},
journal = {Phys. Rev. D},
volume = {91},
issue = {10},
pages = {102003},
numpages = {26},
year = {2015},
month = {May},
publisher = {American Physical Society},
doi = {10.1103/PhysRevD.91.102003},
url = {https://link.aps.org/doi/10.1103/PhysRevD.91.102003}
}
......@@ -23,6 +23,7 @@
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#define _COMPUTE_FSTAT_C
#include "ComputeFstat_internal.h"
......@@ -37,6 +38,7 @@
// Internal definition of input data structure
struct tagFstatInput {
REAL8 Tsft; // Length of input SFTs (for maximum length checking)
REAL8 minFreqFull; // Minimum frequency loaded from input SFTs
REAL8 maxFreqFull; // Maximum frequency loaded from input SFTs
int singleFreqBin; // True if XLALComputeFstat() can only compute a single frequency bin, due to zero dFreq being passed to XLALCreateFstatInput()
......@@ -84,6 +86,73 @@ const FstatOptionalArgs FstatOptionalArgsDefaults = {
// ==================== Function definitions =================== //
///
/// Compute the maximum SFT length that can safely be used as input to XLALComputeFstat(), given
/// the desired range limits in search parameters.
///
/// The \f$\mathcal{F}\f$-statistic algorithms implemented in this module assume that the input
/// SFTs supplied to XLALCreateFstatInput() are of such a length that, within the time span of
/// each any SFT, the spectrum of any CW signal being search for will be mostly contained within
/// one SFT bin:
/// * The \a Demod algorithm make explicit use of this assumption in order to sum up only a few
/// bins in the Dirichlet kernel, thereby speeding up the computation.
/// * While the \a Resamp algorithm is in principle independent of the SFT length (as the data
/// are converted from SFTs back into a time series) it practise it makes use of this assumption
/// for convenience to linearly approximate the phase over the time span of a single SFT.
///
/// In order for this assumption to be valid, the SFT bins must be of a certain minimum size, and
/// hence the time spans of the SFTs must be of a certain maximum length. This maximum length is
/// dependent upon the parameter space being searched: searches for isolated CW sources rarely
/// encounter this restriction by using standard 1800-second SFTs, but searches for CW sources in
/// binary systems typically must use shorter SFTs.
///
/// An expression giving the maximum allowed SFT length is given by \cite LeaciPrix2015 , eq. (C2):
/// \f[
/// T_{\textrm{SFT-max}} = \sqrt{ \frac{
/// 6 \sqrt{ 5 \mu_{\textrm{SFT}} }
/// }{
/// \pi (a \sin\iota / c)_{\textrm{max}} f_{\textrm{max}} \Omega_{\textrm{max}}^2
/// } } \,,
/// \f]
/// This function computes this expression with \f$\mu_{\textrm{SFT}} = 0.01\f$, and
/// \f$(a \sin\iota / c)_{\textrm{max}}\f$ and \f$\Omega_{\textrm{max}}\f$ set to either the binary
/// motion of the source, as specified by \p binaryMaxAsini and \p binaryMinPeriod, or the sidereal
/// motion of the Earth, i.e. with \f$(a \sin\iota / c) \sim 0.02~\textrm{s}\f$ and
/// \f$\Omega \sim 2\pi / 1~\textrm{day}\f$.
///
REAL8 XLALFstatMaximumSFTLength ( const REAL8 maxFreq, /**< [in] Maximum signal frequency */
const REAL8 binaryMaxAsini, /**< [in] Maximum projected semi-major axis a*sini/c (= 0 for isolated sources) */
const REAL8 binaryMinPeriod /**< [in] Minimum orbital period (s) */
)
{
// Check input
XLAL_CHECK_REAL8( maxFreq > 0, XLAL_EINVAL );
XLAL_CHECK_REAL8( binaryMaxAsini >= 0, XLAL_EINVAL );
XLAL_CHECK_REAL8( (binaryMaxAsini == 0) || (binaryMinPeriod > 0), XLAL_EINVAL ); // if binary: P>0
// Use 1% mismatch in maximum allowed SFT length expression
const REAL8 mu_SFT = 0.01;
// Compute maximum allowed SFT length due to sidereal motion of the Earth
const REAL8 earthAsini = LAL_REARTH_SI / LAL_C_SI;
const REAL8 earthOmega = LAL_TWOPI / LAL_DAYSID_SI;
const REAL8 Tsft_max_earth = sqrt( ( 6.0 * sqrt( 5.0 * mu_SFT ) ) / ( LAL_PI * earthAsini * maxFreq * earthOmega * earthOmega ) );
if ( binaryMaxAsini == 0 ) {
return Tsft_max_earth;
}
// Compute maximum allowed SFT length due to binary motion of the source
const REAL8 binaryMaxOmega = LAL_TWOPI / binaryMinPeriod;
const REAL8 Tsft_max_binary = sqrt( ( 6.0 * sqrt( 5.0 * mu_SFT ) ) / ( LAL_PI * binaryMaxAsini * maxFreq * binaryMaxOmega * binaryMaxOmega ) );
// Compute maximum allowed SFT length
const REAL8 Tsft_max = GSL_MIN( Tsft_max_earth, Tsft_max_binary );
return Tsft_max;
} // XLALFstatMaximumSFTLength()
///
/// Create a #FstatInputVector of the given length, for example for setting up
/// F-stat searches over several segments.
......@@ -344,8 +413,8 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Cata
}
// Determine the time baseline of an SFT
const REAL8 Tsft = 1.0 / SFTcatalog->data[0].header.deltaF;
// Determine the length of an SFT
input->Tsft = 1.0 / SFTcatalog->data[0].header.deltaF;
// Compute the mid-time and time-span of the SFTs
double Tspan = 0;
......@@ -353,7 +422,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Cata
const LIGOTimeGPS startTime = SFTcatalog->data[0].header.epoch;
const LIGOTimeGPS endTime = SFTcatalog->data[SFTcatalog->length - 1].header.epoch;
common->midTime = startTime;
Tspan = Tsft + XLALGPSDiff( &endTime, &startTime );
Tspan = input->Tsft + XLALGPSDiff( &endTime, &startTime );
XLALGPSAdd ( &common->midTime, 0.5 * Tspan );
}
......@@ -371,11 +440,11 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Cata
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;
const REAL8 extraFreqMethod = extraBinsMethod / input->Tsft;
minFreqMethod = minCoverFreq - extraFreqMethod;
maxFreqMethod = maxCoverFreq + extraFreqMethod;
const REAL8 extraFreqFull = extraBinsFull / Tsft;
const REAL8 extraFreqFull = extraBinsFull / input->Tsft;
input->minFreqFull = minCoverFreq - extraFreqFull;
input->maxFreqFull = maxCoverFreq + extraFreqFull;
......@@ -470,7 +539,7 @@ XLALCreateFstatInput ( const SFTCatalog *SFTcatalog, ///< [in] Cata
XLAL_CHECK_NULL ( XLALMultiSFTVectorResizeBand ( multiSFTs, minFreqMethod, maxFreqMethod - minFreqMethod ) == XLAL_SUCCESS, XLAL_EFUNC );
// Get detector states, with a timestamp shift of Tsft/2
const REAL8 tOffset = 0.5 * Tsft;
const REAL8 tOffset = 0.5 * input->Tsft;
XLAL_CHECK_NULL ( (common->multiDetectorStates = XLALGetMultiDetectorStates ( common->multiTimestamps, &common->detectors, ephemerides, tOffset )) != NULL, XLAL_EFUNC );
// Save ephemerides and SSB precision
......@@ -623,6 +692,14 @@ XLALComputeFstat ( FstatResults **Fstats, ///< [in/out] Address of
XLAL_CHECK ( !input->singleFreqBin || numFreqBins == 1, XLAL_EINVAL, "numFreqBins must be 1 if XLALCreateFstatInput() was passed zero dFreq" );
XLAL_CHECK ( whatToCompute < FSTATQ_LAST, XLAL_EINVAL);
// Check that SFT length is within allowed maximum
{
const REAL8 maxFreq = doppler->fkdot[0] + input->common.dFreq * numFreqBins;
const REAL8 Tsft_max = XLALFstatMaximumSFTLength( maxFreq, doppler->asini, doppler->period );
XLAL_CHECK ( !XLALIsREAL8FailNaN( Tsft_max ), XLAL_EINVAL );
XLAL_CHECK ( input->Tsft < Tsft_max, XLAL_EINVAL, "Length of input SFTs (%g s) must be less than %g s for CW signal with frequency = %g, binary asini = %g, period = %g", input->Tsft, Tsft_max, maxFreq, doppler->asini, doppler->period );
}
// Allocate results struct, if needed
if ( (*Fstats) == NULL ) {
XLAL_CHECK ( ((*Fstats) = XLALCalloc ( 1, sizeof(**Fstats) )) != NULL, XLAL_ENOMEM );
......
......@@ -329,6 +329,8 @@ typedef struct tagFstatTimingModel
} FstatTimingModel;
// ---------- API function prototypes ----------
REAL8 XLALFstatMaximumSFTLength ( const REAL8 maxFreq, const REAL8 binaryMaxAsini, const REAL8 binaryMinPeriod );
int XLALFstatMethodIsAvailable ( FstatMethodType method );
const CHAR *XLALFstatMethodName ( FstatMethodType method );
const UserChoices *XLALFstatMethodChoices ( void );
......
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