Commit c3b731c6 authored by David Keitel's avatar David Keitel
Browse files

HierarchSearchGCT: fix bug with --outputFX when not all segments have data for all detectors

Combination of 7 commits, original commit messages as follows:
Made error messages in XLALComputeFstatFromAtoms more explicit
HierarchSearchGCT: for outputFX option, fix case of segments where not all detectors have data
-changes in XLALComputeExtraStatsSemiCoherent to skip these segments for concerned detectors
-new XLALGetNumDetectors to make sure maximum detector number is found
-somewhat dirty fix right now, L1 and H1 hardcoded as detector IDs
-test script now has some L1 data missing to check these changes
Fix test script, simplify code for segment/detector fix
-testGCTLV.sh: use global ephemeris files
-LineVeto.c, XLALComputeExtraStatsSemiCoherent: removed unnecessary check, temp variable instead of array for numDetectors per segment, simplified a loop
-LineVeto.c, XLALGetNumDetectors: removed double const qualifier
HierarchSearchGCT with --outputFX: more flexible detector name handling
-new data type NameVector with XLALCreateNameVector and XLALDestroyNameVector
-new function XLALGetDetectorIDs that reads in detector names from multiSFTVector
-changes in XLALComputeExtraStatsSemiCoherent to use this
HierarchSearchGCT with --outputFX: allow for arbitrary detector names and numbers
-get all detectors from SFTs
-changed data type NameVector to NameList with .elems <= .length
-XLALGetDetectorIDs now also gets total numDetectors
-corresponding changes in XLALComputeExtraStatsSemiCoherent
-test script now has L1 in first segment, H1 in last, both in others
further changes in detector name handling
-now get IDs only once outside toplist candidate loop
-therefore moved from XLALComputeExtraStatsSemiCoherent to XLALComputeExtraStatsForToplist
-some smaller simplifications
GCT Line Veto - make ready for patch
-LineVeto.c: removed debug printf statements
-LineVeto.h: removed obsolete prototype
-testGCTLV.sh: add timing, simplify H1/L1 treatment, reset nCand etc.
Original: 64f228736b657b4eebd11a5eca40919b4c338d9e
parent f25b0f0a
......@@ -97,14 +97,29 @@ int XLALComputeExtraStatsForToplist ( toplist_t *list,
PulsarDopplerParams candidateDopplerParams = empty_PulsarDopplerParams; /* struct containing sky position, frequency and fdot for the current candidate */
PulsarSpins fkdotTMP; /* temporary spin parameters for XLALExtrapolatePulsarSpins */
REAL8 deltaTau; /* temporary variable to convert LIGOTimeGPS into real number difference for XLALExtrapolatePulsarSpins */
UINT4 X;
/* initialize doppler parameters */
candidateDopplerParams.refTime = tMidGPS; /* spin parameters will be given to ComputeFStat at this refTime */
INIT_MEM( fkdotTMP );
deltaTau = XLALGPSDiff( &candidateDopplerParams.refTime, &refTimeGPS );
/* initialise detector name vector for later identification */
NameList *detectorIDs;
if ( ( detectorIDs = XLALCreateNameList ( 5 ) ) == NULL ) {
XLALPrintError ("%s: failed to XLALCreateNameList( %d )\n", fn, 5 );
XLAL_ERROR ( fn, XLAL_EFUNC );
}
for ( X = 0; X < detectorIDs->length; X++)
{
CHAR name[LALNameLength] = "INVALIDNAME";
detectorIDs->data[X] = name;
}
XLALGetDetectorIDs ( detectorIDs, multiSFTsV ); /* fill detector name vector with all detectors present in any data sements */
UINT4 numDetectors = detectorIDs->elems;
/* initialise LVcomponents structure and allocate memory */
UINT4 numDetectors = multiSFTsV->data[0]->length; /* number of different single detectors */
LVcomponents lineVeto; /* struct containing multi-detector Fstat, single-detector Fstats, Line Veto stat */
if ( (lineVeto.TwoFX = XLALCreateREAL8Vector ( numDetectors )) == NULL ) {
XLALPrintError ("%s: failed to XLALCreateREAL8Vector( %d )\n", fn, numDetectors );
......@@ -136,7 +151,7 @@ int XLALComputeExtraStatsForToplist ( toplist_t *list,
}
/* recalculate multi- and single-IFO Fstats for all segments for this candidate */
XLALComputeExtraStatsSemiCoherent( &lineVeto, &candidateDopplerParams, multiSFTsV, multiNoiseWeightsV, multiDetStatesV, CFparams, SignalOnly );
XLALComputeExtraStatsSemiCoherent( &lineVeto, &candidateDopplerParams, multiSFTsV, multiNoiseWeightsV, multiDetStatesV, detectorIDs, CFparams, SignalOnly );
if ( xlalErrno != 0 ) {
XLALPrintError ("\nError in function %s, line %d : Failed call to XLALComputeLineVetoSemiCoherent().\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EFUNC );
......@@ -144,7 +159,6 @@ int XLALComputeExtraStatsForToplist ( toplist_t *list,
/* save values in toplist */
elem->sumTwoFnew = lineVeto.TwoF;
UINT4 X;
for ( X = 0; X < numDetectors; X ++ )
elem->sumTwoFX->data[X] = lineVeto.TwoFX->data[X];
......@@ -152,6 +166,7 @@ int XLALComputeExtraStatsForToplist ( toplist_t *list,
/* free temporary structures */
XLALDestroyREAL8Vector ( lineVeto.TwoFX );
XLALDestroyNameList ( detectorIDs );
return (XLAL_SUCCESS);
......@@ -169,6 +184,7 @@ int XLALComputeExtraStatsSemiCoherent ( LVcomponents *lineVeto,
const MultiSFTVectorSequence *multiSFTsV, /**< data files (SFTs) for all detectors and segments */
const MultiNoiseWeightsSequence *multiNoiseWeightsV, /**< noise weights for all detectors and segments */
const MultiDetectorStateSeriesSequence *multiDetStatesV,/**< some state info for all detectors */
const NameList *detectorIDs, /**< name strings of all detectors present in multiSFTsV */
const ComputeFParams *CFparams, /**< additional parameters needed for ComputeFStat */
const BOOLEAN SignalOnly /**< flag for case with no noise, makes some extra checks necessary */
)
......@@ -180,27 +196,24 @@ int XLALComputeExtraStatsSemiCoherent ( LVcomponents *lineVeto,
XLALPrintError ("\nError in function %s, line %d : Empty pointer as input parameter!\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EFAULT);
}
if ( !multiSFTsV->data[0] ) {
XLALPrintError ("\nError in function %s, line %d : Input multiSFT vector has no elements!\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EFAULT);
}
if ( multiSFTsV->length == 0 ) {
XLALPrintError ("\nError in function %s, line %d : Input multiSFT vector over segments has zero length!\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EBADLEN);
}
/* fake LAL status structure, needed as long as ComputeFStat is LAL function and not XLAL */
LALStatus fakeStatus = blank_status;
Fcomponents Fstat; /* temporary struct for ComputeFStat */
UINT4 X, Y;
UINT4 numSegments = multiSFTsV->length;
UINT4 numDetectors = multiSFTsV->data[0]->length;
UINT4 numDetectors = detectorIDs->elems;
if ( lineVeto->TwoFX->length != numDetectors ) {
XLALPrintError ("\%s, line %d : Inconsistent number of detector in TwoFX: %d, while multiSFTsV: %d!\n\n", fn, __LINE__, lineVeto->TwoFX->length, numDetectors );
XLALPrintError ("\%s, line %d : Inconsistent number of detectors: TwoFX vector has length %d, while detectorID list contains %d elements!\n\n", fn, __LINE__, lineVeto->TwoFX->length, numDetectors );
XLAL_ERROR ( fn, XLAL_EBADLEN );
}
/* fake LAL status structure, needed as long as ComputeFStat is LAL function and not XLAL */
LALStatus fakeStatus = blank_status;
Fcomponents Fstat; /* temporary struct for ComputeFStat */
/* temporary copy of Fstatistic parameters structure, needed to change returnAtoms for function scope only */
ComputeFParams CFparams_internal = (*CFparams);
CFparams_internal.returnAtoms = TRUE;
......@@ -208,13 +221,21 @@ int XLALComputeExtraStatsSemiCoherent ( LVcomponents *lineVeto,
/* initialiase LVcomponents structure */
lineVeto->TwoF = 0.0;
lineVeto->LV = 0.0;
UINT4 X;
for (X = 0; X < numDetectors; X++) {
lineVeto->TwoFX->data[X] = 0.0;
}
REAL8 Tsft = 1.0 / multiSFTsV->data[0]->data[0]->data[0].deltaF; /* get SFT duration */
/* variables necessary to catch segments where not all detectors have data */
INT4 detid = -1; /* count through detector IDs for matching with name strings */
UINT4 numDetectorsSeg = 0; /* number of detectors with data might be different for each segment */
UINT4 numSegmentsX[numDetectors]; /* number of segments with data might be different for each detector */
for (X = 0; X < numDetectors; X++)
{
numSegmentsX[X] = 0;
}
/* compute single- and multi-detector Fstats for each data segment and sum up */
UINT4 k;
for (k = 0; k < numSegments; k++)
......@@ -234,23 +255,38 @@ int XLALComputeExtraStatsSemiCoherent ( LVcomponents *lineVeto,
lineVeto->TwoF += 2.0 * Fstat.F; /* sum up multi-detector Fstat for this segment*/
numDetectorsSeg = multiSFTsV->data[k]->length; /* for each segment, could be smaller than overall number */
/* recompute single-detector Fstats from atoms */
for (X = 0; X < numDetectors; X++)
for (X = 0; X < numDetectorsSeg; X++)
{
REAL8 twoFX = 2.0 * XLALComputeFstatFromAtoms ( &Fstat, X );
if ( xlalErrno != 0 ) {
XLALPrintError ("\nError in function %s, line %d : Failed call to XLALComputeFstatFromAtoms().\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EFUNC );
}
XLALPrintError ("\nError in function %s, line %d : Failed call to XLALComputeFstatFromAtoms().\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EFUNC );
}
if ( SignalOnly ) { /* normalization factor correction (TwoF=2.0*F has been done before, this time!) */
twoFX *= 4.0 / Tsft;
twoFX += 4;
}
lineVeto->TwoFX->data[X] += twoFX; /* sum up single-detector Fstat for this segment*/
/* match detector ID in this segment to one from detectorIDs list, sum up the corresponding twoFX */
detid = -1;
for (Y = 0; Y < numDetectors; Y++)
{
if ( strcmp( multiSFTsV->data[k]->data[X]->data[0].name, detectorIDs->data[Y] ) == 0 )
detid = Y;
}
if ( detid == -1 )
{
XLALPrintError ("\nError in function %s, line %d : For segment k=%d, detector X=%d, could not match detector ID %s.\n\n", fn, __LINE__, k, X, multiSFTsV->data[k]->data[X]->data[0].name);
XLAL_ERROR ( fn, XLAL_EFAILED );
}
numSegmentsX[detid] += 1; /* have to keep this for correct averaging */
lineVeto->TwoFX->data[detid] += twoFX; /* sum up single-detector Fstat for this segment*/
} /* for X < numDetectors X */
} /* for X < numDetectorsSeg */
/* free memory for atoms that was allocated within ComputeFStat */
XLALDestroyMultiFstatAtomVector ( Fstat.multiFstatAtoms );
......@@ -260,7 +296,7 @@ int XLALComputeExtraStatsSemiCoherent ( LVcomponents *lineVeto,
/* get average stats over all segments */
lineVeto->TwoF /= numSegments;
for (X = 0; X < numDetectors; X++) {
lineVeto->TwoFX->data[X] /= numSegments;
lineVeto->TwoFX->data[X] /= numSegmentsX[X];
}
return(XLAL_SUCCESS);
......@@ -282,17 +318,17 @@ REAL8 XLALComputeFstatFromAtoms ( const Fcomponents *Fstat, /**< multi-detecto
}
if ( Fstat->multiFstatAtoms->length == 0 ) {
XLALPrintError ("\nError in function %s, line %d : Input MultiFstatAtomVector has zero length!\n\n", fn, __LINE__);
XLALPrintError ("\nError in function %s, line %d : Input MultiFstatAtomVector has zero length! (no detectors)\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EBADLEN );
}
if ( X > Fstat->multiFstatAtoms->length-1 ) {
XLALPrintError ("\nError in function %s, line %d : Invalid detector number!\n\n", fn, __LINE__);
XLALPrintError ("\nError in function %s, line %d : Invalid detector number!\nRequested X=%d, but FstatAtoms only have length %d.\n\n", fn, __LINE__, X, Fstat->multiFstatAtoms->length);
XLAL_ERROR ( fn, XLAL_EDOM );
}
if ( Fstat->multiFstatAtoms->data[X]->length == 0 ) {
XLALPrintError ("\nError in function %s, line %d : Input FstatAtomVector has zero length!\n\n", fn, __LINE__);
XLALPrintError ("\nError in function %s, line %d : Input FstatAtomVector has zero length! (no timestamps for detector X=%d)\n\n", fn, __LINE__, X);
XLAL_ERROR ( fn, XLAL_EDOM );
}
......@@ -399,3 +435,113 @@ REAL8 XLALComputeLineVeto ( const REAL8 TwoF, /**< multi-detector Fsta
} /* XLALComputeLineVeto() */
/** XLAL function to get a list of detector IDs from multi-segment multiSFT vectors
* returns all unique detector IDs for cases with some detectors switching on and off
* detectorIDs structure needs to be properly initialised with length >= total number of detectors present in data
*/
int XLALGetDetectorIDs ( NameList *detectorIDs, /** [out] vector of IFO name strings */
const MultiSFTVectorSequence *multiSFTsV /**< data files (SFTs) for all detectors and segments */
)
{
const char *fn = __func__;
/* check input parameters and report errors */
if ( !multiSFTsV || !multiSFTsV->data ) {
XLALPrintError ("\nError in function %s, line %d : Empty pointer as input parameter!\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EFAULT);
}
if ( multiSFTsV->length == 0 ) {
XLALPrintError ("\nError in function %s, line %d : Input multiSFT vector contains no segments (zero length)!\n\n", fn, __LINE__);
XLAL_ERROR ( fn, XLAL_EBADLEN);
}
/* set up variables */
UINT4 k, X, Y;
UINT4 numSegments = multiSFTsV->length; /* number of segments with SFTs from any detector */
UINT4 numDetectorsSeg = 0; /* number of detectors with data might be different for each segment */
BOOLEAN detIDfound = FALSE; /* for each k and X, check if the detector ID was already in the list */
for (k = 0; k < numSegments; k++)
{
if ( !multiSFTsV->data[k] || (multiSFTsV->data[k]->length == 0) ) {
XLALPrintError ("\nError in function %s, line %d : Input multiSFT vector, segment k=%d has no elements (0 detectors)!\n\n", fn, __LINE__, k);
XLAL_ERROR ( fn, XLAL_EFAULT);
}
numDetectorsSeg = multiSFTsV->data[k]->length; /* for each segment, could be smaller than overall number */
for (X = 0; X < numDetectorsSeg; X++)
{
/* check if this k, X combination yields a new detector ID */
Y = 0;
detIDfound = FALSE;
while ( ( Y < detectorIDs->length ) && ( detIDfound == FALSE ) )
{
if ( strcmp( multiSFTsV->data[k]->data[X]->data[0].name, detectorIDs->data[Y] ) == 0 )
detIDfound = TRUE; /* ID already in list, do nothing */
Y++;
}
if ( detIDfound == FALSE ) /* ID not in list yet, add it */
{
detectorIDs->data[detectorIDs->elems] = multiSFTsV->data[k]->data[X]->data[0].name;
detectorIDs->elems++;
if ( detectorIDs->elems > detectorIDs->length )
{
XLALPrintError ("\nError in function %s, line %d : already found %d detector IDs, exceeded ID list length %d. This happened for segment k=%d, detector X=%d!\n\n", fn, __LINE__, detectorIDs->elems, detectorIDs->length, k, X);
XLAL_ERROR ( fn, XLAL_EBADLEN);
}
}
} /* for X < numDetectorsSeg */
} /* for k < numSegments */
return(XLAL_SUCCESS);
} /* XLALGetDetectorIDs() */
NameList * XLALCreateNameList ( UINT4 length )
{
NameList *list;
list = LALMalloc( sizeof( *list ) );
if ( ! list )
XLAL_ERROR_NULL( "XLALCreateNameList", XLAL_ENOMEM );
list->length = length;
list->elems = 0;
if ( ! length ) /* zero length: set data pointer to be NULL */
list->data = NULL;
else /* non-zero length: allocate memory for data */
{
list->data = LALMalloc( length * sizeof( *list->data ) );
if ( ! list->data )
{
LALFree( list );
XLAL_ERROR_NULL( "XLALCreateNameList", XLAL_ENOMEM );
}
}
return list;
}
void XLALDestroyNameList ( NameList *list )
{
if ( ! list )
return;
if ( ( ! list->length || ! list->data ) && ( list->length || list->data ) )
XLAL_ERROR_VOID( "XLALDestroyNameList", XLAL_EINVAL );
if ( list->data )
LALFree( list->data );
list->data = NULL; /* leave length non-zero to detect repeated frees */
LALFree( list );
return;
}
......@@ -79,6 +79,12 @@ extern "C" {
REAL8Vector *TwoFX; /**< vector of single-detector F-statistic values */
} LVcomponents;
/** Type containing a vector of IFO name strings */
typedef struct tagNameList {
UINT4 length; /**< length (maximal number of entries) of the list (how many detectors expected in theory?) */
UINT4 elems; /**< number of elements really in the list (how many detectors where found in segment SFTs?) */
CHAR **data; /**< name strings of the ifos */
} NameList;
/*---------- exported Global variables ----------*/
/* empty init-structs for the types defined in here */
......@@ -101,6 +107,7 @@ XLALComputeExtraStatsSemiCoherent ( LVcomponents *lineVeto,
const MultiSFTVectorSequence *multiSFTs,
const MultiNoiseWeightsSequence *multiNoiseWeights,
const MultiDetectorStateSeriesSequence *multiDetStates,
const NameList *detectorIDs,
const ComputeFParams *CFparams,
const BOOLEAN SignalOnly);
......@@ -114,6 +121,16 @@ XLALComputeLineVeto ( const REAL8 TwoF,
const REAL8 rhomax,
const REAL8Vector *priorX );
int
XLALGetDetectorIDs ( NameList *detectorIDs,
const MultiSFTVectorSequence *multiSFTsV );
NameList
*XLALCreateNameList ( UINT4 length );
void
XLALDestroyNameList ( NameList *list );
#ifdef __cplusplus
}
#endif
......
......@@ -83,8 +83,8 @@ gct_dFreq="0.000002" #"2.0e-6"
gct_dF1dot="1.0e-10"
gct_nCands="1000" # number of candidates in output
edat="earth05-09.dat"
sdat="sun05-09.dat"
edat="earth09-11.dat"
sdat="sun09-11.dat"
RngMedWindow=101 # running median window, needs to be equal for HSGCT and CFS
## --------- Generate fake data set time stamps -------------
......@@ -94,16 +94,17 @@ echo "----------------------------------------------------------------------"
echo
Tsft="1800"
startTime="852443819"
refTime="862999869"
startTime="952443819"
refTime="962999869"
Tsegment="90000"
Nsegments="14"
## insert gaps between segments
seggap=$(echo $Tsegment | awk '{printf "%f", $1 * 1.12345}');
# midTime="853731034.5"
tsfile="./timestamps.txt" # for makefakedata
segfile="./segments.txt" # for HSGCT
rm -rf $tsfile $segfile
tsfileH1="./timestampsH1.dat" # for makefakedata
tsfileL1="./timestampsL1.dat" # for makefakedata
segfile="./segments.dat" # for HSGCT
rm -rf $tsfileH1 $tsfileL1 $segfile
tmpTime=$startTime
ic1="1"
while [ "$ic1" -le "$Nsegments" ];
......@@ -112,6 +113,12 @@ do
t1=`echo $t0 $Tsegment | awk '{print $1 + $2}'`
TspanHours=`echo $Tsegment | awk '{printf "%.7f", $1 / 3600.0 }'`
NSFT=`echo $Tsegment $Tsft | awk '{print int(2.0 * $1 / $2 + 0.5) }'`
if [ `echo $ic1 | awk '{if($1==1) {print "1"}}'` ];then
NSFT=`echo $Tsegment $Tsft | awk '{print int(1.0 * $1 / $2 + 0.5) }'`
fi
if [ `echo $ic1" "$Nsegments | awk '{if($1==$2) {print "1"}}'` ];then
NSFT=`echo $Tsegment $Tsft | awk '{print int(1.0 * $1 / $2 + 0.5) }'`
fi
echo "$t0 $t1 $TspanHours $NSFT" >> $segfile
segs[${ic1}]=$tmpTime # save seg's beginning for later use
echo "Segment: "$ic1" of "$Nsegments" GPS start time: "${segs[${ic1}]}
......@@ -119,7 +126,12 @@ do
ic2=$Tsft
while [ "$ic2" -le "$Tsegment" ];
do
echo ${tmpTime}" 0" >> $tsfile
if [ `echo $ic1 | awk '{if($1>1) {print "1"}}'` ];then
echo ${tmpTime}" 0" >> $tsfileH1
fi
if [ `echo $ic1" "$Nsegments | awk '{if($1<$2) {print "1"}}'` ];then
echo ${tmpTime}" 0" >> $tsfileL1
fi
tmpTime=$(echo $tmpTime $Tsft | awk '{printf "%.0f", $1 + $2}');
ic2=$(echo $ic2 $Tsft | awk '{printf "%.0f", $1 + $2}');
done
......@@ -141,10 +153,10 @@ else
fi
## construct MFD cmdline
mfd_CL=" --fmin=$mfd_fmin --Band=$mfd_FreqBand --Freq=$Freq --outSFTbname=$SFTdir --f1dot=$f1dot --Alpha=$Alpha --Delta=$Delta --psi=$psi --phi0=$phi0 --cosi=$cosi --ephemYear=05-09 --generationMode=1 --timestampsFile=$tsfile --refTime=$refTime --Tsft=$Tsft --noiseSqrtSh=$noiseSqrtSh"
mfd_CL=" --fmin=$mfd_fmin --Band=$mfd_FreqBand --Freq=$Freq --outSFTbname=$SFTdir --f1dot=$f1dot --Alpha=$Alpha --Delta=$Delta --psi=$psi --phi0=$phi0 --cosi=$cosi --ephemYear=09-11 --generationMode=1 --refTime=$refTime --Tsft=$Tsft --noiseSqrtSh=$noiseSqrtSh"
## detector H1
cmdline="$mfd_code $mfd_CL --IFO=H1 --h0=$h0H1 --randSeed=1000";
cmdline="$mfd_code $mfd_CL --IFO=H1 --h0=$h0H1 --timestampsFile=$tsfileH1 --randSeed=1000";
echo $cmdline;
if ! eval $cmdline; then
echo "Error.. something failed when running '$mfd_code' ..."
......@@ -152,7 +164,7 @@ if ! eval $cmdline; then
fi
## detector L1:
cmdline="$mfd_code $mfd_CL --IFO=L1 --h0=$h0L1 --randSeed=1001";
cmdline="$mfd_code $mfd_CL --IFO=L1 --h0=$h0L1 --timestampsFile=$tsfileL1 --randSeed=1001";
echo $cmdline;
if ! eval $cmdline; then
echo "Error.. something failed when running '$mfd_code' ..."
......@@ -172,7 +184,7 @@ fi
outfile_gct="HS_GCT_LV.dat"
gct_CL=" -d1 --fnameout=$outfile_gct --gridType1=3 --nCand1=$gct_nCands --skyRegion='allsky' --Freq=$Freq --DataFiles='$SFTfiles' --ephemE=$edat --ephemS=$sdat --skyGridFile='./$skygridfile' --printCand1 --semiCohToplist --df1dot=$gct_dF1dot --f1dot=$f1dot --f1dotBand=$gct_F1dotBand --dFreq=$gct_dFreq --FreqBand=$gct_FreqBand --refTime=$refTime --segmentList=$segfile --outputFX --blocksRngMed=$RngMedWindow"
gct_CL=" -d1 --fnameout=$outfile_gct --gridType1=3 --nCand1=$gct_nCands --skyRegion='allsky' --Freq=$Freq --DataFiles='$SFTfiles' --ephemE=$edat --ephemS=$sdat --skyGridFile='./$skygridfile' --printCand1 --semiCohToplist --df1dot=$gct_dF1dot --f1dot=$f1dot --f1dotBand=$gct_F1dotBand --dFreq=$gct_dFreq --FreqBand=$gct_FreqBand --refTime=$refTime --segmentList=$segfile --outputFX --blocksRngMed=$RngMedWindow --outputTiming='./timing.dat'"
cmdline="$gct_code $gct_CL > >(tee stdout.log) 2> >(tee stderr.log >&2)"
echo $cmdline
......@@ -194,8 +206,8 @@ echo
outfile_cfs="fstat_loudest.dat"
## initialise summation variables
twoFcfs="0"
twoFcfsX1="0"
twoFcfsX2="0"
twoFcfsH1="0"
twoFcfsL1="0"
for ((x=1; x <= $Nsegments; x++))
do
......@@ -206,7 +218,7 @@ for ((x=1; x <= $Nsegments; x++))
echo "Segment: "$x" "$startGPS" "$endGPS
## construct ComputeFStatistic command lines
cfs_CL="--DataFiles='$SFTfiles' --outputLoudest='$outfile_cfs' --TwoFthreshold=0.0 --ephemYear=05-09 --refTime=$refTime --Alpha=$AlphaSearch --Delta=$DeltaSearch --Freq=$freqGCT --f1dot=$f1dotGCT --minStartTime=$startGPS --maxEndTime=$endGPS --RngMedWindow=$RngMedWindow"
cfs_CL="--DataFiles='$SFTfiles' --outputLoudest='$outfile_cfs' --TwoFthreshold=0.0 --ephemYear=09-11 --refTime=$refTime --Alpha=$AlphaSearch --Delta=$DeltaSearch --Freq=$freqGCT --f1dot=$f1dotGCT --minStartTime=$startGPS --maxEndTime=$endGPS --RngMedWindow=$RngMedWindow"
## multi-IFO
cmdline="$cfs_code -v0 $cfs_CL"
......@@ -219,32 +231,35 @@ for ((x=1; x <= $Nsegments; x++))
twoFcfs=$(echo $twoFcfs $twoFcfs_seg | awk '{printf "%.6f", $1 + $2}');
## detector H1
cmdline="$cfs_code -v0 $cfs_CL --IFO='H1'"
echo $cmdline
if ! eval $cmdline; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1
fi
twoFcfs_seg=$(sed 's/\;//' $outfile_cfs | awk '{if($1=="twoF"){printf "%.6f",$3}}')
twoFcfsX1=$(echo $twoFcfsX1 $twoFcfs_seg | awk '{printf "%.6f", $1 + $2}');
if [ `echo $x | awk '{if($1>1) {print "1"}}'` ];then
cmdline="$cfs_code -v0 $cfs_CL --IFO='H1'"
echo $cmdline
if ! eval $cmdline; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1
fi
twoFcfs_seg=$(sed 's/\;//' $outfile_cfs | awk '{if($1=="twoF"){printf "%.6f",$3}}')
twoFcfsH1=$(echo $twoFcfsH1 $twoFcfs_seg | awk '{printf "%.6f", $1 + $2}');
fi
## detector L1
cmdline="$cfs_code -v0 $cfs_CL --IFO='L1'"
echo $cmdline
if ! eval $cmdline; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1
if [ `echo $x" "$Nsegments | awk '{if($1<$2) {print "1"}}'` ];then
cmdline="$cfs_code -v0 $cfs_CL --IFO='L1'"
echo $cmdline
if ! eval $cmdline; then
echo "Error.. something failed when running '$cfs_code' ..."
exit 1
fi
twoFcfs_seg=$(sed 's/\;//' $outfile_cfs | awk '{if($1=="twoF"){printf "%.6f",$3}}')
twoFcfsL1=$(echo $twoFcfsL1 $twoFcfs_seg | awk '{printf "%.6f", $1 + $2}');
fi
twoFcfs_seg=$(sed 's/\;//' $outfile_cfs | awk '{if($1=="twoF"){printf "%.6f",$3}}')
twoFcfsX2=$(echo $twoFcfsX2 $twoFcfs_seg | awk '{printf "%.6f", $1 + $2}');
done
## get averages
twoFcfs=$(echo $twoFcfs $Nsegments | awk '{printf "%.6f", $1/$2}');
twoFcfsX1=$(echo $twoFcfsX1 $Nsegments | awk '{printf "%.6f", $1/$2}');
twoFcfsX2=$(echo $twoFcfsX2 $Nsegments | awk '{printf "%.6f", $1/$2}');
echo "==> CFS at f="$freqGCT": F ="$twoFcfs" F1="$twoFcfsX1" F2="$twoFcfsX2
twoFcfsH1=$(echo $twoFcfsH1 $Nsegments | awk '{printf "%.6f", $1/($2-1)}');
twoFcfsL1=$(echo $twoFcfsL1 $Nsegments | awk '{printf "%.6f", $1/($2-1)}');
echo
echo "----------------------------------------------------------------------"
......@@ -274,66 +289,66 @@ echo
## get Fstats values from GCT output file
twoFGCT=$(sed -e '/%/d;' $outfile_gct | sort -nr -k6,6 | head -1 | awk '{print $6}')
twoFGCTLV=$(sed -e '/%/d;' $outfile_gct | sort -nr -k6,6 | head -1 | awk '{print $7}')
twoFGCTLVX1=$(sed -e '/%/d;' $outfile_gct | sort -nr -k6,6 | head -1 | awk '{print $8}')
twoFGCTLVX2=$(sed -e '/%/d;' $outfile_gct | sort -nr -k6,6 | head -1 | awk '{print $9}')
twoFGCTLVH1=$(sed -e '/%/d;' $outfile_gct | sort -nr -k6,6 | head -1 | awk '{print $9}')
twoFGCTLVL1=$(sed -e '/%/d;' $outfile_gct | sort -nr -k6,6 | head -1 | awk '{print $8}')
## compute relative errors of Fstats from CFS, GTC plain and GTCLV
reldev_cfs_gct=$(echo $twoFcfs $twoFGCT | awk '{ if(($1-$2)>=0) {printf "%.12f", ($1-$2)/(0.5*($1+$2))} else {printf "%.12f", (-1)*($1-$2)/(0.5*($1+$2))}}');
reldev_gct_LV=$(echo $twoFGCT $twoFGCTLV | awk '{ if(($1-$2)>=0) {printf "%.12f", ($1-$2)/(0.5*($1+$2))} else {printf "%.12f", (-1)*($1-$2)/(0.5*($1+$2))}}');
reldev_cfs_LV=$(echo $twoFcfs $twoFGCTLV | awk '{ if(($1-$2)>=0) {printf "%.12f", ($1-$2)/(0.5*($1+$2))} else {printf "%.12f", (-1)*($1-$2)/(0.5*($1+$2))}}');
reldev_cfs_LVX1=$(echo $twoFcfsX1 $twoFGCTLVX1 | awk '{ if(($1-$2)>=0) {printf "%.12f", ($1-$2)/(0.5*($1+$2))} else {printf "%.12f", (-1)*($1-$2)/(0.5*($1+$2))}}');
reldev_cfs_LVX2=$(echo $twoFcfsX2 $twoFGCTLVX2 | awk '{ if(($1-$2)>=0) {printf "%.12f", ($1-$2)/(0.5*($1+$2))} else {printf "%.12f", (-1)*($1-$2)/(0.5*($1+$2))}}');
reldev_cfs_LVH1=$(echo $twoFcfsH1 $twoFGCTLVH1 | awk '{ if(($1-$2)>=0) {printf "%.12f", ($1-$2)/(0.5*($1+$2))} else {printf "%.12f", (-1)*($1-$2)/(0.5*($1+$2))}}');
reldev_cfs_LVL1=$(echo $twoFcfsL1 $twoFGCTLVL1 | awk '{ if(($1-$2)>=0) {printf "%.12f", ($1-$2)/(0.5*($1+$2))} else {printf "%.12f", (-1)*($1-$2)/(0.5*($1+$2))}}');
## get maximum deviations (over all candidates) of freq, Fstat between plain GCT and LV
maxdev_gct_LV=$(sed '/^%.*/d' $outfile_gct | awk 'BEGIN { maxDiff = 0; } { relDiff=($6-$7)/($6+$7); if( relDiff^2 > maxDiff^2 ) {maxDiff=relDiff}; } END {printf "%.12f", maxDiff}' )
maxdevfreq_gct_LV=$(sed '/^%.*/d' $outfile_gct | awk 'BEGIN { maxDiff = 0; maxFreq = 0;} { relDiff=($6-$7)/($6+$7); if( relDiff^2 > maxDiff^2 ) {maxDiff=relDiff; maxFreq=$1}; } END {printf "%.12f", maxFreq}' )
echo "==> CFS at f="$freqGCT" : F="$twoFcfs" F1="$twoFcfsX1" F2="$twoFcfsX2
echo "==> CFS at f="$freqGCT" : F="$twoFcfs" FH1="$twoFcfsH1" FL1="$twoFcfsL1
## check plain GCT search code output 2F against externally computed 2F
if [ `echo $reldev_cfs_gct" "$Tolerance | awk '{if($1>$2) {print "1"}}'` ];then
echo "==> GCT plain : F ="$twoFGCT" (diff vs cfs: "$reldev_cfs_gct")"
echo "==> GCT plain : F ="$twoFGCT" (diff vs cfs: "$reldev_cfs_gct")"
echo "OUCH... results differ by more than tolerance limit. Something might be wrong..."
exit 2
else
echo "==> GCT plain : F ="$twoFGCT" (diff vs cfs: "$reldev_cfs_gct") OK."
echo "==> GCT plain : F ="$twoFGCT" (diff vs cfs: "$reldev_cfs_gct") OK."
fi
## check LV-recomputed 2F against plain GCT value
if [ `echo $reldev_gct_LV" "$Tolerance | awk '{if($1>$2) {print "1"}}'` ];then
echo "==> GCT LV : F ="$twoFGCTLV" (diff vs plain: "$reldev_gct_LV")"
echo "==> GCT LV : F ="$twoFGCTLV" (diff vs plain: "$reldev_gct_LV")"
echo " (maximum diff: "$maxdev_gct_LV" at freq="$maxdevfreq_gct_LV")"
echo "OUCH... results differ by more than tolerance limit. Something might be wrong..."
exit 2
else
echo "==> GCT LV : F ="$twoFGCTLV" (diff vs plain: "$reldev_gct_LV") OK."
echo " (maximum diff: "$maxdev_gct_LV" at freq="$maxdevfreq_gct_LV")"
echo "==> GCT LV : F ="$twoFGCTLV" (diff vs plain: "$reldev_gct_LV") OK."
echo " (maximum diff: "$maxdev_gct_LV" at freq="$maxdevfreq_gct_LV")"
fi
## check recomputed 2F against external CFS
if [ `echo $reldev_cfs_LV" "$Tolerance | awk '{if($1>$2) {print "1"}}'` ];then
echo " (diff vs cfs: "$reldev_cfs_LV")"
echo " (diff vs cfs: "$reldev_cfs_LV")"
echo "OUCH... results differ by more than tolerance limit. Something might be wrong..."
exit 2
else
echo " (diff vs cfs: "$reldev_cfs_LV") OK."
echo " (diff vs cfs: "$reldev_cfs_LV") OK."
fi
## check single-detector F1 against external CFS (at GCT freq)
if [ `echo $reldev_cfs_LVX1" "$Tolerance | awk '{if($1>$2) {print "1"}}'` ];then
echo "==> GCT LV : F1="$twoFGCTLVX1" (diff vs cfs: "$reldev_cfs_LVX1")"
if [ `echo $reldev_cfs_LVH1" "$Tolerance | awk '{if($1>$2) {print "1"}}'` ];then
echo "==> GCT LV : FH1="$twoFGCTLVH1" (diff vs cfs: "$reldev_cfs_LVH1")"
echo "OUCH... results differ by more than tolerance limit. Something might be wrong..."
exit 2
else
echo "==> GCT LV : F1="$twoFGCTLVX1" (diff vs cfs: "$reldev_cfs_LVX1") OK."
echo "==> GCT LV : FH1="$twoFGCTLVH1" (diff vs cfs: "$reldev_cfs_LVH1") OK."
fi
## check single-detector F2 against external CFS (at GCT freq)
if [ `echo $reldev_cfs_LVX2" "$Tolerance | awk '{if($1>$2) {print "1"}}'` ];then
echo "==> GCT LV : F2="$twoFGCTLVX2" (diff vs cfs: "$reldev_cfs_LVX2")"
if [ `echo $reldev_cfs_LVL1" "$Tolerance | awk '{if($1>$2) {print "1"}}'` ];then
echo "==> GCT LV : FL1="$twoFGCTLVL1" (diff vs cfs: "$reldev_cfs_LVL1")"
echo "OUCH... results differ by more than tolerance limit. Something might be wrong..."
exit 2
else
echo "==> GCT LV : F2="$twoFGCTLVX2" (diff vs cfs: "$reldev_cfs_LVX2") OK."
echo "==> GCT LV : FL1="$twoFGCTLVL1" (diff vs cfs: "$reldev_cfs_LVL1") OK."
fi
......@@ -341,6 +356,6 @@ echo "----------------------------------------------------------------------"
## clean up files
if [ -z "$NOCLEANUP" ]; then
rm -rf $SFTdir $skygridfile $tsfile $segfile $outfile_gct $outfile_cfs checkpoint.cpt stderr.log stdout.log
rm -rf $SFTdir $skygridfile $tsfileH1 $tsfileL1 $segfile $outfile_gct $outfile_cfs checkpoint.cpt stderr.log stdout.log
echo "Cleaned up."
fi