update fast likelihood function to include the pmt response and absorption

1 files changed, 11 insertions, 3 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index fc1cf54..92c8306 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -17,6 +17,7 @@
 #include "solid_angle.h"
 #include <gsl/gsl_roots.h>
 #include <gsl/gsl_errno.h>
+#include "pmt_response.h"
 
 typedef struct intParams {
     path *p;
@@ -152,7 +153,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl
      *
      * `smax` is currently calculated as the point where the particle velocity
      * drops to 0.8 times the speed of light. */
-    double pmt_dir[3], tmp[3], R, cos_theta, theta, x, z, s, a, b, beta, E, p, T, omega, theta_cerenkov, n, sin_theta, E0, p0, beta0;
+    double pmt_dir[3], tmp[3], R, cos_theta, theta, x, z, s, a, b, beta, E, p, T, omega, theta_cerenkov, n, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length, wavelength0;
 
     /* Calculate beta at the start of the track. */
     E0 = T0 + MUON_MASS;
@@ -173,7 +174,8 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl
     theta = acos(cos_theta);
 
     /* Compute the Cerenkov angle at the start of the track. */
-    n = get_index_snoman_d2o(400.0);
+    wavelength0 = 400.0;
+    n = get_index_snoman_d2o(wavelength0);
     theta_cerenkov = acos(1/(n*beta0));
 
     /* Now, we compute the distance along the track where the PMT is at the
@@ -236,6 +238,8 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl
 
     SUB(pmt_dir,pmts[i].pos,tmp);
 
+    cos_theta_pmt = DOT(pmt_dir,pmts[i].normal)/NORM(pmt_dir);
+
     /* Calculate the sine of the angle between the track direction and the PMT
      * at the position `s` along the track. */
     sin_theta = fabs(sin(acos(DOT(dir,pmt_dir)/NORM(pmt_dir))));
@@ -247,7 +251,11 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl
 
     double frac = sqrt(2)*n*x*beta0*theta0;
 
-    return n*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n*(smax-z)*beta0)/frac) + erf((-a+b*s+n*z*beta0)/frac))/(b+n*beta0)/(4*M_PI);
+    f = get_weighted_pmt_response(acos(-cos_theta_pmt));
+
+    absorption_length = get_absorption_length_snoman_d2o(wavelength0);
+
+    return f*exp(-a/absorption_length)*n*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n*(smax-z)*beta0)/frac) + erf((-a+b*s+n*z*beta0)/frac))/(b+n*beta0)/(4*M_PI);
 }
 
 double getKineticEnergy(double x, double T0)