update likelihood to calculate absorption length correctly

1 files changed, 6 insertions, 3 deletions

diff --git a/src/muon.c b/src/muon.c
index 5eb1f54..0ab5cd8 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -258,7 +258,7 @@ double get_dEdx(double T, double rho)
 
 double get_expected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
 {
-    double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length, distance;
+    double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, l_h2o, l_d2o, distance;
 
     z = 1.0;
 
@@ -298,8 +298,11 @@ double get_expected_charge(double x, double T, double theta0, double *pos, doubl
 
     f = get_weighted_pmt_response(acos(-cos_theta_pmt));
 
-    absorption_length = get_absorption_length_snoman_d2o(wavelength0);
+    absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0);
+    absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0);
+
+    get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
 
     /* FIXME: add angular response and scattering/absorption. */
-    return f*exp(-distance/absorption_length)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
+    return f*exp(-l_d2o/absorption_length_d2o)*exp(-l_h2o/absorption_length_h2o)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
 }