Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org starting 2 March 2020 at approximately 8am MST. It is expected to take around 10 minutes and will include a short period of downtime towards the end of the maintenance window. Please direct any comments, concerns, or questions to computing-help@igwn.org.

Commit 1f4733d8 authored by John Douglas Veitch's avatar John Douglas Veitch

Hacky version compiling

parent b4c37584
......@@ -7,12 +7,34 @@
#include "LALInferenceDistanceMarg.h"
#include <math.h>
double dist_integral(double rho_opt, double rho_match);
#include <gsl/gsl_integration.h>
double dist_snr_pdf(double rho_opt, double rho_match, double dL)
struct integrand_args
{
double A,B;
};
double dist_integral(double rho_opt, double rho_match, double dist_min, double dist_max)
{
size_t limit=100;
double result,abserr;
gsl_function F;
gsl_integration_workspace *workspace = gsl_integration_workspace_alloc (limit);
F.function = &dist_snr_pdf;
struct integrand_args args = {rho_opt, rho_match};
F.params = &args;
gsl_integration_qag (&F, dist_min, dist_max, 0, 1e-6, limit, GSL_INTEG_GAUSS61, workspace, &result, &abserr);
gsl_integration_workspace_free(workspace);
return(result / (dist_max-dist_min));
}
double dist_snr_pdf(double dL, void *args)
{
double A = (rho_opt / dL)*(rho_opt / dL);
double B = rho_match/dL;
return exp(-A + B);
struct integrand_args *a = (struct integrand_args *)args;
return exp(-a->A/dL/dL + a->B/dL)*dL*dL;
}
......@@ -17,8 +17,8 @@
* rho_opt is optimal SNR at 1 Mpc
* rho_mf is <d|h> at 1 Mpc
*/
double dist_integral(double rho_opt, double rho_match);
double dist_integral(double rho_opt, double rho_match, double dist_min, double dist_max);
double dist_snr_pdf(double rho_opt, double rho_match, double dL);
double dist_snr_pdf(double dL, void *args);
#endif /* LALInferenceDistanceMarg_h */
......@@ -1288,8 +1288,15 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
Dinitial=atof(ppt->value);
distanceVary = LALINFERENCE_PARAM_FIXED;
}
LALInferenceRegisterUniformVariableREAL8(state, model->params, "logdistance", log(Dinitial), log(Dmin), log(Dmax), distanceVary);
if(LALInferenceGetProcParamVal(commandLine,"--margdist"))
{
REAL8 a = log(Dmin), b=log(Dmax);
LALInferenceAddMinMaxPrior(model->params, "logdistance", &a, &b, LALINFERENCE_REAL8_t);
}
else
LALInferenceRegisterUniformVariableREAL8(state, model->params, "logdistance", log(Dinitial), log(Dmin), log(Dmax), distanceVary) ;
LALInferenceRegisterUniformVariableREAL8(state, model->params, "polarisation", zero, psiMin, psiMax, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "costheta_jn", zero, costhetaJNmin, costhetaJNmax,LALINFERENCE_PARAM_LINEAR);
......
......@@ -33,6 +33,7 @@
#include <lal/Sequence.h>
#include <lal/FrequencySeries.h>
#include <lal/TimeFreqFFT.h>
#include <lal/LALInferenceDistanceMarg.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_dawson.h>
......@@ -1200,9 +1201,19 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
REAL8 OptimalSNR=sqrt(2.0*S);
REAL8 MatchedFilterSNR = 0.;
/* Avoid nan's, since noise-only model has OptimalSNR == 0. */
if (OptimalSNR > 0.)
MatchedFilterSNR = 2.0*d_inner_h/OptimalSNR;
if(1){
double dist_min, dist_max;
LALInferenceGetMinMaxPrior(model->params, "logdistance", &dist_min, &dist_max);
double marg_l = dist_integral(OptimalSNR, 2.0*d_inner_h, exp(dist_min), exp(dist_max));
loglikelihood = -D + log(marg_l);
}
LALInferenceAddVariable(currentParams,"optimal_snr",&OptimalSNR,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
LALInferenceAddVariable(currentParams,"matched_filter_snr",&MatchedFilterSNR,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
......
......@@ -709,8 +709,11 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
exit(0);
}
distance = exp(LALInferenceGetREAL8Variable(model->params, "logdistance"))* LAL_PC_SI * 1.0e6; /* distance (1 Mpc) in units of metres */
if(!LALInferenceCheckVariable(model->params,"logdistance")) distance=LAL_PC_SI * 1e6; /* If distance not given, 1Mpc used */
else
{
distance = exp(LALInferenceGetREAL8Variable(model->params, "logdistance"))* LAL_PC_SI * 1.0e6; /* distance (1 Mpc) in units of metres */
}
phi0 = LALInferenceGetREAL8Variable(model->params, "phase"); /* START phase as per lalsimulation convention, radians*/
/* Zenith angle between J and N in radians. Also known as inclination angle when spins are aligned */
......
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