Commit d3358b80 authored by Salvatore Vitale's avatar Salvatore Vitale Committed by John Douglas Veitch

Fix LALInferencePrintInjectionSample

parent 66906277
......@@ -156,6 +156,9 @@ int main(int argc, char *argv[]){
/* write injection with noise evidence information from algorithm */
LALInferencePrintInjectionSample(state);
if (LALInferenceGetProcParamVal(state->commandLine, "--dont-sample")){
return(0);
}
/* Call nested sampling algorithm */
state->algorithm(state);
......
......@@ -1805,6 +1805,7 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
char SNRpath[FILENAME_MAX+50];
ProcessParamsTable *ppt=NULL;
int flipped_masses=0;
LALInferenceIFOData *dataPtr;
ppt = LALInferenceGetProcParamVal(commandLine,"--outfile");
if (ppt)
......@@ -1895,6 +1896,8 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
*/
}
/* FIXME: One also need to code the same manipulations done to f_max in LALInferenceTemplateXLALSimInspiralChooseWaveform line 856 circa*/
/* Set up LAL dictionary */
LALDict* LALpars=XLALCreateDict();
......@@ -1931,7 +1934,17 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
if(LALInferenceGetProcParamVal(commandLine,"--inj-fref")) {
fref = atoi(LALInferenceGetProcParamVal(commandLine,"--inj-fref")->value);
}
if (approximant == TaylorF2)
f_max = 0.0; /* this will stop at ISCO */
else{
/* get the max f_max as done in Template. Has to do this since I cannot access LALInferenceModel here*/
dataPtr=IFOdata;
while(dataPtr){
if(dataPtr->fHigh>f_max) f_max=dataPtr->fHigh;
dataPtr=dataPtr->next;
}
}
/* Print a line with information about approximant, amp_order, phaseorder, tide order and spin order */
fprintf(stdout,"\n\n---\t\t ---\n");
fprintf(stdout,"Injection will run using Approximant %i (%s), phase order %i, amp order %i, spin order %i, tidal order %i, in the frequency domain.\n",approximant,XLALSimInspiralGetStringFromApproximant(approximant),phase_order,amp_order,(int) spinO,(int) tideO);
......@@ -1959,7 +1972,6 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
exit(1);
}
LALInferenceIFOData *dataPtr;
REAL8 Fplus, Fcross;
REAL8 plainTemplateReal, plainTemplateImag;
REAL8 templateReal, templateImag;
......@@ -2099,7 +2111,7 @@ static void PrintSNRsToFile(LALInferenceIFOData *IFOdata , char SNRpath[] ){
*/
void LALInferenceInjectionToVariables(SimInspiralTable *theEventTable, LALInferenceVariables *vars)
{
UINT4 spinCheck=0;
//UINT4 spinCheck=0;
if(!vars) {
XLALPrintError("Encountered NULL variables pointer");
XLAL_ERROR_VOID(XLAL_EINVAL);
......@@ -2108,46 +2120,13 @@ void LALInferenceInjectionToVariables(SimInspiralTable *theEventTable, LALInfere
REAL8 q = theEventTable->mass2 / theEventTable->mass1;
if (q > 1.0) q = 1.0/q;
REAL8 sx = theEventTable->spin1x;
REAL8 sy = theEventTable->spin1y;
REAL8 s1z = theEventTable->spin1z;
REAL8 a_spin1 = sqrt(sx*sx + sy*sy + s1z*s1z);
REAL8 theta_spin1, phi_spin1;
if (a_spin1 == 0.0) {
theta_spin1 = 0.0;
phi_spin1 = 0.0;
} else {
theta_spin1 = acos(s1z / a_spin1);
phi_spin1 = atan2(sy, sx);
if (phi_spin1 < 0.0) phi_spin1 += 2.0*M_PI;
}
sx = theEventTable->spin2x;
sy = theEventTable->spin2y;
REAL8 s2z = theEventTable->spin2z;
REAL8 a_spin2 = sqrt(sx*sx + sy*sy + s2z*s2z), theta_spin2, phi_spin2;
if (a_spin2 == 0.0) {
theta_spin2 = 0.0;
phi_spin2 = 0.0;
} else {
theta_spin2 = acos(s2z / a_spin2);
phi_spin2 = atan2(sy, sx);
if (phi_spin2 < 0.0) phi_spin2 += 2.0*M_PI;
}
/* Check for presence of spin in the injection */
if(a_spin1!=0.0 || a_spin2!=0.0) spinCheck=1;
REAL8 psi = theEventTable->polarization;
if (psi>=M_PI) psi -= M_PI;
REAL8 injGPSTime = XLALGPSGetREAL8(&(theEventTable->geocent_end_time));
REAL8 dist = theEventTable->distance;
REAL8 cosinclination = cos(theEventTable->inclination);
//REAL8 cosinclination = cos(theEventTable->inclination);
REAL8 phase = theEventTable->coa_phase;
REAL8 dec = theEventTable->latitude;
REAL8 ra = theEventTable->longitude;
......@@ -2158,13 +2137,18 @@ void LALInferenceInjectionToVariables(SimInspiralTable *theEventTable, LALInfere
REAL8 m1=theEventTable->mass1;
REAL8 m2=theEventTable->mass2;
REAL8 chirpmass = theEventTable->mchirp;
LALInferenceAddVariable(vars, "mass1", &m1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "mass2", &m2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "chirpmass", &chirpmass, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "q", &q, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "time", &injGPSTime, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "distance", &dist, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "costheta_jn", &cosinclination, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"distance"))
LALInferenceAddVariable(vars, "distance", &dist, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
else if (LALInferenceCheckVariable(vars,"logdistance")){
REAL8 logdistance=log(dist);
LALInferenceAddVariable(vars, "logdistance", &logdistance, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
//LALInferenceAddVariable(vars, "costheta_jn", &cosinclination, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "polarisation", &(psi), LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "phase", &phase, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "declination", &dec, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
......@@ -2172,20 +2156,34 @@ void LALInferenceInjectionToVariables(SimInspiralTable *theEventTable, LALInfere
LALInferenceAddVariable(vars, "LAL_APPROXIMANT", &injapprox, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "LAL_PNORDER",&order,LALINFERENCE_INT4_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "LAL_AMPORDER", &(theEventTable->amp_order), LALINFERENCE_INT4_t, LALINFERENCE_PARAM_FIXED);
if(spinCheck){
if (theEventTable->spin1x==0 && theEventTable->spin1y==0.0 && theEventTable->spin2x==0.0 && theEventTable->spin2y==0.0){
LALInferenceAddVariable(vars, "a_spin1", &s1z, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "a_spin2", &s2z, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
else{
LALInferenceAddVariable(vars, "a_spin1", &a_spin1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "a_spin2", &a_spin2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "theta_spin1", &theta_spin1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "theta_spin2", &theta_spin2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "phi_spin1", &phi_spin1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
LALInferenceAddVariable(vars, "phi_spin2", &phi_spin2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
}
REAL8 thetaJN,phiJL,theta1,theta2,phi12,chi1,chi2;
/* Convert cartesian spin coordinates to system-frame variables*/
/* This fref is --inj-fref, which has been set previously by LALInferencePrintInjectionSample
I don't call it inj_fref since LALInferenceTemplate looks for fref, and that is what will be called to calculate
the logL at injval
*/
REAL8 fref=100.0;
if (LALInferenceCheckVariable(vars,"f_ref"))
fref= *(REAL8*) LALInferenceGetVariable(vars,"f_ref");
XLALSimInspiralTransformPrecessingWvf2PE(&thetaJN,&phiJL,&theta1,&theta2,&phi12,&chi1,&chi2,theEventTable->inclination,theEventTable->spin1x,theEventTable->spin1y,theEventTable->spin1z, theEventTable->spin2x, theEventTable->spin2y, theEventTable->spin2z,m1,m2,fref,phase);
if (LALInferenceCheckVariable(vars,"a_spin1"))
LALInferenceAddVariable(vars,"a_spin1", &chi1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"a_spin2"))
LALInferenceAddVariable(vars,"a_spin2", &chi2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"tilt_spin1"))
LALInferenceAddVariable(vars,"tilt_spin1", &theta1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"tilt_spin2"))
LALInferenceAddVariable(vars,"tilt_spin2", &theta2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"phi_jl"))
LALInferenceAddVariable(vars,"phi_jl", &phiJL, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
if (LALInferenceCheckVariable(vars,"phi12"))
LALInferenceAddVariable(vars,"phi12", &phi12, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
REAL8 costhetajn=cos(thetaJN);
LALInferenceAddVariable(vars, "costheta_jn", &costhetajn, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}
......@@ -2254,10 +2252,10 @@ void LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
}
}
LALInferenceAddVariable(injparams, "logL", (void *)&injL,LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
if (LALInferenceCheckVariable(runState->algorithmParams, "logZnoise")){
REAL8 tmp=injL-*(REAL8 *)LALInferenceGetVariable(runState->algorithmParams,"logZnoise");
LALInferenceAddVariable(injparams,"deltalogL",(void *)&tmp,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
}
REAL8 logZnoise=LALInferenceNullLogLikelihood(runState->data);
REAL8 tmp2=injL-logZnoise;
LALInferenceAddVariable(injparams,"deltalogL",(void *)&tmp2,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
LALInferenceIFOData *data=runState->data;
while(data) {
char tmpName[320];
......@@ -2266,6 +2264,13 @@ void LALInferencePrintInjectionSample(LALInferenceRunState *runState) {
LALInferenceAddVariable(injparams, tmpName, &tmp, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
data=data->next;
}
REAL8 fref = 100.;
if(LALInferenceGetProcParamVal(runState->commandLine,"--inj-fref")) {
fref = atoi(LALInferenceGetProcParamVal(runState->commandLine,"--inj-fref")->value);
}
LALInferenceAddVariable(injparams,"f_ref",(void *)&fref,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
/* Save to file */
outfile=fopen(fname,"w");
if(!outfile) {fprintf(stderr,"ERROR: Unable to open file %s for injection saving\n",fname); exit(1);}
......
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