Commit 334d43b8 authored by David Keitel's avatar David Keitel Committed by Reinhard Prix

CFSv2: additional output option for full transient-CW grid in (t0,tau)

 -signal starttime t0 as GPS epoch [seconds]
 -signal duration tau also in seconds
Original: 12dac9f9c315c00b55d13bf1b50c7e9adbe44f6c
parent 05af5dd2
......@@ -293,6 +293,7 @@ typedef struct {
REAL8 BSGLthreshold; /**< output threshold on BSGL */
CHAR *outputTransientStats; /**< output file for transient B-stat values */
CHAR *outputTransientStatsAll; /**< output file for transient B-stat values */
CHAR *transient_WindowType; /**< name of transient window ('none', 'rect', 'exp',...) */
REAL8 transient_t0Days; /**< earliest GPS start-time for transient window search, as offset in days from dataStartGPS */
REAL8 transient_t0DaysBand; /**< Range of GPS start-times to search in transient search, in days */
......@@ -360,7 +361,7 @@ BOOLEAN XLALCenterIsLocalMax ( const scanlineWindow_t *scanWindow, const UINT4 s
*/
int main(int argc,char *argv[])
{
FILE *fpFstat = NULL, *fpTransientStats = NULL;
FILE *fpFstat = NULL, *fpTransientStats = NULL, *fpTransientStatsAll = NULL;
REAL8 numTemplates, templateCounter;
time_t clock0;
PulsarDopplerParams XLAL_INIT_DECL(dopplerpos);
......@@ -446,6 +447,13 @@ int main(int argc,char *argv[])
XLAL_CHECK_MAIN ( write_transientCandidate_to_fp ( fpTransientStats, NULL ) == XLAL_SUCCESS, XLAL_EFUNC ); /* write header-line comment */
}
if ( uvar.outputTransientStatsAll )
{
XLAL_CHECK_MAIN ( (fpTransientStatsAll = fopen (uvar.outputTransientStatsAll, "wb")) != NULL, XLAL_ESYS, "\nError opening file '%s' for writing..\n\n", uvar.outputTransientStatsAll );
fprintf (fpTransientStatsAll, "%s", GV.logstring ); /* write search log comment */
XLAL_CHECK_MAIN ( write_transientCandidateAll_to_fp ( fpTransientStatsAll, NULL ) == XLAL_SUCCESS, XLAL_EFUNC ); /* write header-line comment */
}
/* start Fstatistic histogram with a single empty bin */
if (uvar.outputFstatHist) {
XLAL_CHECK_MAIN ((Fstat_histogram = gsl_vector_int_alloc(1)) != NULL, XLAL_ENOMEM );
......@@ -688,7 +696,7 @@ int main(int argc,char *argv[])
} /* if outputFstatAtoms */
/* ----- compute transient-CW statistics if their output was requested ----- */
if ( fpTransientStats )
if ( fpTransientStats || fpTransientStatsAll )
{
transientCandidate_t XLAL_INIT_DECL(transientCand);
......@@ -712,11 +720,14 @@ int main(int argc,char *argv[])
XLAL_CHECK_MAIN ( (pdf_t0 = XLALComputeTransientPosterior_t0 ( GV.transientWindowRange, transientCand.FstatMap )) != NULL, XLAL_EFUNC );
XLAL_CHECK_MAIN ( (pdf_tau = XLALComputeTransientPosterior_tau ( GV.transientWindowRange, transientCand.FstatMap )) != NULL, XLAL_EFUNC );
/* get maximum-posterior estimate (MP) from the modes of these pdfs */
transientCand.t0_MP = XLALFindModeOfPDF1D ( pdf_t0 );
XLAL_CHECK_MAIN ( xlalErrno == 0, XLAL_EFUNC, "mode-estimation failed for pdf_t0. xlalErrno = %d\n", xlalErrno );
transientCand.tau_MP = XLALFindModeOfPDF1D ( pdf_tau );
XLAL_CHECK_MAIN ( xlalErrno == 0, XLAL_EFUNC, "mode-estimation failed for pdf_tau. xlalErrno = %d\n", xlalErrno );
if ( fpTransientStats )
{
/* get maximum-posterior estimate (MP) from the modes of these pdfs */
transientCand.t0_MP = XLALFindModeOfPDF1D ( pdf_t0 );
XLAL_CHECK_MAIN ( xlalErrno == 0, XLAL_EFUNC, "mode-estimation failed for pdf_t0. xlalErrno = %d\n", xlalErrno );
transientCand.tau_MP = XLALFindModeOfPDF1D ( pdf_tau );
XLAL_CHECK_MAIN ( xlalErrno == 0, XLAL_EFUNC, "mode-estimation failed for pdf_tau. xlalErrno = %d\n", xlalErrno );
}
/* record timing-relevant transient search params */
timing.tauMin = GV.transientWindowRange.tau;
......@@ -732,15 +743,24 @@ int main(int argc,char *argv[])
if ( uvar.SignalOnly )
transientCand.FstatMap->maxF += 2;
/* output everything into stats-file (one line per candidate) */
XLAL_CHECK_MAIN ( write_transientCandidate_to_fp ( fpTransientStats, &transientCand ) == XLAL_SUCCESS, XLAL_EFUNC );
if ( fpTransientStats )
{
/* output everything into stats-file (one line per candidate) */
XLAL_CHECK_MAIN ( write_transientCandidate_to_fp ( fpTransientStats, &transientCand ) == XLAL_SUCCESS, XLAL_EFUNC );
}
if ( fpTransientStatsAll )
{
/* output everything into stats-file (one block for the whole (t0,tau) grid per candidate) */
XLAL_CHECK_MAIN ( write_transientCandidateAll_to_fp ( fpTransientStatsAll, &transientCand ) == XLAL_SUCCESS, XLAL_EFUNC );
}
/* free dynamically allocated F-stat map */
XLALDestroyTransientFstatMap ( transientCand.FstatMap );
XLALDestroyPDF1D ( pdf_t0 );
XLALDestroyPDF1D ( pdf_tau );
} /* if fpTransientStats */
} /* if ( fpTransientStats || fpTransientStatsAll ) */
} // for ( iFreq < numFreqBins_FBand )
......@@ -812,6 +832,7 @@ int main(int argc,char *argv[])
fpFstat = NULL;
}
if ( fpTransientStats ) fclose (fpTransientStats);
if ( fpTransientStatsAll ) fclose (fpTransientStatsAll);
/* ----- estimate amplitude-parameters for the loudest canidate and output into separate file ----- */
if ( uvar.outputLoudest )
......@@ -1052,6 +1073,7 @@ initUserVars ( UserInput_t *uvar )
// --------------------------------------------
XLALRegisterUvarMember(outputTransientStats,STRING, 0, OPTIONAL, "TransientCW: Output filename for transient-CW statistics.");
XLALRegisterUvarMember(outputTransientStatsAll,STRING, 0, DEVELOPER, "TransientCW: Output filename for transient-CW statistics -- including the whole (t0,tau) grid for each candidate -- WARNING: CAN BE HUGE!.");
XLALRegisterUvarMember( transient_WindowType,STRING, 0,OPTIONAL, "TransientCW: Type of transient signal window to use. ('none', 'rect', 'exp').");
XLALRegisterUvarMember( transient_t0Days, REAL8, 0, OPTIONAL, "TransientCW: Earliest GPS start-time for transient window search, as offset in days from dataStartGPS");
XLALRegisterUvarMember( transient_t0DaysBand,REAL8, 0,OPTIONAL, "TransientCW: Range of GPS start-times to search in transient search, in days");
......@@ -1481,7 +1503,7 @@ InitFstat ( ConfigVariables *cfg, const UserInput_t *uvar )
}
/* get atoms back from Fstat-computing, either if atoms-output or transient-Bstat output was requested */
if ( ( uvar->outputFstatAtoms != NULL ) || ( uvar->outputTransientStats != NULL ) ) {
if ( ( uvar->outputFstatAtoms != NULL ) || ( uvar->outputTransientStats != NULL ) || ( uvar->outputTransientStatsAll != NULL ) ) {
cfg->Fstat_what |= FSTATQ_ATOMS_PER_DET;
}
......
......@@ -999,7 +999,60 @@ write_transientCandidate_to_fp ( FILE *fp, const transientCandidate_t *thisCand
return XLAL_SUCCESS;
} /* write_TransCandidate_to_fp() */
} /* write_transientCandidate_to_fp() */
/**
* Write full set of t0 and tau grid points for given transient CW candidate into output file.
*
* NOTE: if input thisCand == NULL, we write a header comment-line explaining the fields
*
*/
int
write_transientCandidateAll_to_fp ( FILE *fp, const transientCandidate_t *thisCand )
{
/* sanity checks */
if ( !fp ) {
XLALPrintError ( "%s: invalid NULL filepointer input.\n", __func__ );
XLAL_ERROR ( XLAL_EINVAL );
}
if ( thisCand == NULL ) /* write header-line comment */
{
fprintf (fp, "%%%% Freq[Hz] Alpha[rad] Delta[rad] fkdot[1] fkdot[2] fkdot[3] t0[s] tau[s] twoF\n");
}
else
{
if ( !thisCand->FstatMap ) {
XLALPrintError ("%s: incomplete: transientCand->FstatMap == NULL!\n", __func__ );
XLAL_ERROR ( XLAL_EINVAL );
}
UINT4 t0 = thisCand->windowRange.t0;
UINT4 N_t0Range = (UINT4) floor ( thisCand->windowRange.t0Band / thisCand->windowRange.dt0 ) + 1;
UINT4 N_tauRange = (UINT4) floor ( thisCand->windowRange.tauBand / thisCand->windowRange.dtau ) + 1;
LogPrintf ( LOG_DETAIL, "N_t0Range=%d, N_tauRange=%d\n", N_t0Range, N_tauRange);
/* ----- OUTER loop over start-times [t0,t0+t0Band] ---------- */
for ( UINT4 m = 0; m < N_t0Range; m ++ ) /* m enumerates 'binned' t0 start-time indices */
{
UINT4 this_t0 = t0 + m * thisCand->windowRange.dt0;
/* ----- INNER loop over timescale-parameter tau ---------- */
for ( UINT4 n = 0; n < N_tauRange; n ++ )
{
UINT4 this_tau = thisCand->windowRange.tau + n * thisCand->windowRange.dtau;
fprintf (fp, " %- 18.16f %- 19.16f %- 19.16f %- 9.6g %- 9.5g %- 9.5g %10d %10d %- 11.8g\n",
thisCand->doppler.fkdot[0], thisCand->doppler.Alpha, thisCand->doppler.Delta,
thisCand->doppler.fkdot[1], thisCand->doppler.fkdot[2], thisCand->doppler.fkdot[3],
this_t0, this_tau, 2.0 * gsl_matrix_get ( thisCand->FstatMap->F_mn, m, n )
);
}
}
}
return XLAL_SUCCESS;
} /* write_transientCandidateAll_to_fp() */
/**
......
......@@ -109,6 +109,8 @@ int XLALApplyTransientWindow2NoiseWeights ( MultiNoiseWeights *multiNoiseWeights
int write_transientCandidate_to_fp ( FILE *fp, const transientCandidate_t *thisTransCand );
int write_transientCandidateAll_to_fp ( FILE *fp, const transientCandidate_t *thisTransCand );
transientFstatMap_t *XLALComputeTransientFstatMap ( const MultiFstatAtomVector *multiFstatAtoms,
transientWindowRange_t windowRange,
......
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