diff --git a/lalapps/src/pulsar/Xray/SemiCoherentBinary_v2.c b/lalapps/src/pulsar/Xray/SemiCoherentBinary_v2.c
index 9e7c498b5c65a7f2487a31dce64e758b12bd3054..49fc9a218e804c0a31dfd1f1782748541c00dd59 100644
--- a/lalapps/src/pulsar/Xray/SemiCoherentBinary_v2.c
+++ b/lalapps/src/pulsar/Xray/SemiCoherentBinary_v2.c
@@ -95,6 +95,7 @@ typedef struct {
REAL8 tsft; /**< the length of the input sfts */
CHAR *comment;
CHAR *tempdir; /**< a temporary directory for keeping the results */
+ BOOLEAN with_chi2_renorm; /**< switch on chi^2 renormalisation */
BOOLEAN with_xbins; /**< enable fast summing of extra bins */
} UserInput_t;
@@ -113,7 +114,7 @@ extern int vrbflg; /**< defined in lalapps.c */
/* define functions */
int main(int argc,char *argv[]);
int XLALReadUserVars(int argc,char *argv[],UserInput_t *uvar, CHAR **clargs);
-int XLALComputeSemiCoherentStat(FILE *fp,REAL4DemodulatedPowerVector *power,ParameterSpace *pspace,GridParametersVector *fgrid,GridParameters *bingrid,INT4 ntoplist,BOOLEAN with_xbins);
+int XLALComputeSemiCoherentStat(FILE *fp,REAL4DemodulatedPowerVector *power,ParameterSpace *pspace,GridParametersVector *fgrid,GridParameters *bingrid,INT4 ntoplist,BOOLEAN with_chi2_renorm,BOOLEAN with_xbins);
int XLALDefineBinaryParameterSpace(REAL8Space **, LIGOTimeGPS, REAL8, UserInput_t *);
int XLALOpenSemiCoherentResultsFile(FILE **,CHAR *,ParameterSpace *,CHAR *,UserInput_t *);
int XLALtoplist(REAL4 x,Template *params,toplist *TL);
@@ -369,7 +370,7 @@ int main( int argc, char *argv[] ) {
/**********************************************************************************/
/* compute the semi-coherent detection statistic on the fine grid */
- if (XLALComputeSemiCoherentStat(sfp,dmpower,&pspace,freqgridparams,bingridparams,uvar.ntoplist,uvar.with_xbins)) {
+ if (XLALComputeSemiCoherentStat(sfp,dmpower,&pspace,freqgridparams,bingridparams,uvar.ntoplist,uvar.with_chi2_renorm,uvar.with_xbins)) {
LogPrintf(LOG_CRITICAL,"%s : XLALComputeSemiCoherentStat() failed with error = %d\n",__func__,xlalErrno);
return 1;
}
@@ -464,6 +465,7 @@ int XLALReadUserVars(int argc, /**< [in] the command line argument co
uvar->tsft = 256;
uvar->seed = 1;
uvar->tempdir = NULL;
+ uvar->with_chi2_renorm = 1;
uvar->with_xbins = 1;
/* initialise all parameter space ranges to zero */
@@ -498,6 +500,7 @@ int XLALReadUserVars(int argc, /**< [in] the command line argument co
XLALRegisterUvarMember(seed, INT4, 'X', OPTIONAL, "The random number seed (0 = clock)");
XLALRegisterUvarMember(gpsstart, INT4, 's', OPTIONAL, "The minimum start time (GPS sec)");
XLALRegisterUvarMember(gpsend, INT4, 'e', OPTIONAL, "The maximum end time (GPS sec)");
+ XLALRegisterUvarMember(with_chi2_renorm, BOOLEAN, 0, OPTIONAL, "Switch on chi^2 renormalisation");
XLALRegisterUvarMember(with_xbins, BOOLEAN, 0, DEVELOPER, "Enable fast summing of extra bins");
/* do ALL cmdline and cfgfile handling */
@@ -537,6 +540,7 @@ int XLALComputeSemiCoherentStat(FILE *fp, /**< [i
GridParametersVector *fgrid, /**< UNDOCUMENTED */
GridParameters *bingrid, /**< [in] the grid parameters */
INT4 ntoplist, /**< UNDOCUMENTED */
+ BOOLEAN with_chi2_renorm, /**< switch on chi^2 renormalisation */
BOOLEAN with_xbins /**< enable fast summing of extra bins */
)
{
@@ -777,8 +781,10 @@ int XLALComputeSemiCoherentStat(FILE *fp, /**< [i
var = var/(REAL8)(Ntemp-1);
var = (var - mean*mean);
- /* modify toplist to have correct mean and variance */
- for (i=0;i<(UINT4)TL.n;i++) TL.data[i] = (TL.data[i] - mean)*sqrt(4.0*power->length)/sqrt(var) + 2.0*power->length;
+ if ( with_chi2_renorm ) {
+ /* modify toplist to have correct mean and variance */
+ for (i=0;i<(UINT4)TL.n;i++) TL.data[i] = (TL.data[i] - mean)*sqrt(4.0*power->length)/sqrt(var) + 2.0*power->length;
+ }
/* output toplist to file */
for (i=0;i<(UINT4)TL.n;i++) {
diff --git a/lalapps/src/pulsar/Xray/SemiCoherentUtils.c b/lalapps/src/pulsar/Xray/SemiCoherentUtils.c
index 3882dfc2df1e2a01334e2ace50a6003b465a5e85..6ebc894f0467561a1c499cd718aa9cec8045ac38 100644
--- a/lalapps/src/pulsar/Xray/SemiCoherentUtils.c
+++ b/lalapps/src/pulsar/Xray/SemiCoherentUtils.c
@@ -250,17 +250,17 @@ int XLALGetNextRandomBinaryTemplate(Template **temp, /**<
while (!flag) {
REAL8 temp1 = gsl_ran_flat((gsl_rng*)r,0,1);
- nu = n2*pow((pow(n1/n2,4.0) + (1.0 - pow(n1/n2,4.0))*temp1),1.0/4.0);
+ nu = (n1 < n2) ? n2*pow((pow(n1/n2,4.0) + (1.0 - pow(n1/n2,4.0))*temp1),1.0/4.0) : n1;
REAL8 temp2 = gsl_ran_flat((gsl_rng*)r,0,1);
- a = a2*pow((pow(a1/a2,3.0) + (1.0 - pow(a1/a2,3.0))*temp2),1.0/3.0);
+ a = (a1 < a2) ? a2*pow((pow(a1/a2,3.0) + (1.0 - pow(a1/a2,3.0))*temp2),1.0/3.0) : a1;
REAL8 temp3 = gsl_ran_flat((gsl_rng*)r,0,1);
- tasc = t1 + (t2-t1)*temp3;
+ tasc = (t1 < t2) ? t1 + (t2-t1)*temp3 : t1;
REAL8 temp4 = gsl_ran_flat((gsl_rng*)r,0,1);
- Om = O2*pow((pow(O1/O2,5.0) + (1.0 - pow(O1/O2,5.0))*temp4),1.0/5.0);
+ Om = (O1 < O2) ? O2*pow((pow(O1/O2,5.0) + (1.0 - pow(O1/O2,5.0))*temp4),1.0/5.0) : O1;
/* LogPrintf(LOG_DEBUG,"%f (%f %f) %f (%f %f) %f (%f %f) %f (%f %f)\n",nu,n1,n2,a,a1,a2,tasc,t1,t2,Om,O1,O2); */
- if ((nu>n1)&&(nu<n2)&&(a>a1)&&(a<a2)&&(tasc>t1)&&(tasc<t2)&&(Om>O1)&&(Om<O2)) flag = 1;
+ if ((nu>=n1)&&(nu<=n2)&&(a>=a1)&&(a<=a2)&&(tasc>=t1)&&(tasc<=t2)&&(Om>=O1)&&(Om<=O2)) flag = 1;
}
@@ -1313,19 +1313,46 @@ int XLALComputeBinaryGridParams(GridParameters **binarygridparams, /**< [out] t
(*binarygridparams)->coverage = coverage;
if ((*binarygridparams)->coverage>0) {
+ ndim = 4;
REAL8 Vn = pow(LAL_PI,ndim/2.0)/gsl_sf_gamma(1.0+ndim/2.0);
+ REAL8 Vsr = 1;
REAL8 G11 = pow(LAL_PI,2.0)*DT*DT/3;
+ REAL8 fmin = space->data[0].min;
+ REAL8 fmax = space->data[0].max;
+ REAL8 fmid = 0.5*(fmin + fmax);
+ fmin = MYMIN(fmin, fmid - sqrt(mu/G11));
+ fmax = MYMAX(fmax, fmid + sqrt(mu/G11));
+ Vsr *= sqrt(G11/mu) * (pow(fmax,4.0) - pow(fmin,4.0)) / 4.0;
+
REAL8 G22 = pow(LAL_PI,2.0)*DT*DT/6;
+ REAL8 amin = space->data[1].min;
+ REAL8 amax = space->data[1].max;
+ REAL8 amid = 0.5*(amin + amax);
+ amin = MYMIN(amin, amid - sqrt(mu/G22));
+ amax = MYMAX(amax, amid + sqrt(mu/G22));
+ Vsr *= sqrt(G22/mu) * (pow(amax,3.0) - pow(amin,3.0)) / 3.0;
+
REAL8 G33 = pow(LAL_PI,2.0)*DT*DT/6;
+ REAL8 tascmin = space->data[2].min;
+ REAL8 tascmax = space->data[2].max;
+ REAL8 tascmid = 0.5*(tascmin + tascmax);
+ tascmin = MYMIN(tascmin, tascmid - sqrt(mu/G33));
+ tascmax = MYMAX(tascmax, tascmid + sqrt(mu/G33));
+ Vsr *= sqrt(G33/mu) * (tascmax - tascmin);
+
REAL8 G44 = pow(LAL_PI,2.0)*DT*DT*T*T/72;
+ REAL8 omegamin = space->data[3].min;
+ REAL8 omegamax = space->data[3].max;
+ REAL8 omegamid = 0.5*(omegamin + omegamax);
+ omegamin = MYMIN(omegamin, omegamid - sqrt(mu/G44));
+ omegamax = MYMAX(omegamax, omegamid + sqrt(mu/G44));
+ Vsr *= sqrt(G44/mu) * (pow(omegamax,5.0) - pow(omegamin,5.0)) / 5.0;
- REAL8 dVr = sqrt(G11*G22*G33*G44);
- REAL8 Vsr = dVr*space->data[2].span*(1.0/60.0)*(pow(space->data[3].max,5.0)-pow(space->data[3].min,5.0))*(pow(space->data[0].max,4.0)-pow(space->data[0].min,4.0))*(pow(space->data[1].max,3.0)-pow(space->data[1].min,3.0));
+ (*binarygridparams)->Nr = (UINT8)ceil( (1.0/Vn) * log(1.0/(1.0-coverage)) * Vsr );
- (*binarygridparams)->Nr = (UINT8)ceil((1.0/Vn)*log(1.0/(1.0-coverage))*(pow(mu,-ndim/2.0))*Vsr);
- LogPrintf(LOG_DEBUG,"%s : computed the number of random binary templates to be %"LAL_UINT8_FORMAT".\n",__func__,(*binarygridparams)->Nr);
- LogPrintf(LOG_DEBUG,"%s : to be compared to the total number of cubic templates %"LAL_UINT8_FORMAT" (%.6f).\n", __func__, (*binarygridparams)->max, (REAL8)(*binarygridparams)->max/(REAL8)(*binarygridparams)->Nr);
+ LogPrintf(LOG_NORMAL,"%s : computed the number of random binary templates to be %"LAL_UINT8_FORMAT".\n",__func__,(*binarygridparams)->Nr);
+ LogPrintf(LOG_NORMAL,"%s : to be compared to the total number of cubic templates %"LAL_UINT8_FORMAT" (%.6f).\n", __func__, (*binarygridparams)->max, (REAL8)(*binarygridparams)->max/(REAL8)(*binarygridparams)->Nr);
}
@@ -1525,8 +1552,22 @@ int XLALBinaryToSFTVector(SFTVector **SFTvect, /**< [out] copied SFT (needs
/**********************************************************************************/
/* extract effective band from this, if neccessary (ie if faster-sampled output SFTs) */
- SFTVector *tempSFTvect = XLALExtractBandFromSFTVector ( sftvect, par->freq, par->freqband-1.0/par->tsft ); /* the last bin has to be avoided */
+ REAL8 df = sftvect->data[0].deltaF;
+ long bin0 = lround( ( par->freq - sftvect->data[0].f0 ) / df);
+ long bin1 = lround( ( par->freq + par->freqband - sftvect->data[0].f0 ) / df);
+ SFTVector *tempSFTvect = NULL;
+ XLAL_CHECK ( (tempSFTvect = XLALCreateSFTVector ( sftvect->length, 0 )) != NULL, XLAL_EFUNC );
+ for ( UINT4 alpha = 0; alpha < sftvect->length; ++alpha ) {
+ SFTtype *dest = &(tempSFTvect->data[alpha]);
+ SFTtype *src = &(sftvect->data[alpha]);
+ *dest = *src;
+ dest->f0 += bin0 * df;
+ XLAL_CHECK ( (dest->data = XLALCreateCOMPLEX8Vector ( bin1 - bin0 )) != NULL, XLAL_EFUNC );
+ memcpy( dest->data->data, src->data->data + bin0, ( bin1 - bin0 ) * sizeof(dest->data->data[0]) );
+ }
XLALDestroySFTVector(sftvect);
+ LogPrintf(LOG_DEBUG, "extracted frequency band %0.16g to %0.16g, bin %li to %li\n",
+ tempSFTvect->data[0].f0, tempSFTvect->data[0].f0 + tempSFTvect->data[0].deltaF * tempSFTvect->data[0].data->length, bin0, bin1);
/**********************************************************************************/
/* append these SFTs to the full list of SFTs */