update likelihood function to calculate time of flight of photons using get_path_length()

1 files changed, 15 insertions, 21 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index fca14cc..7e1bafe 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -89,25 +89,18 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
 static double gsl_muon_time(double x, void *params)
 {
     intParams *pars = (intParams *) params;
-    double dir[3], pos[3], pmt_dir[3], R, n, wavelength0, T, t, theta0, distance;
+    double dir[3], pos[3], n_d2o, n_h2o, wavelength0, T, t, theta0, l_d2o, l_h2o;
 
     path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
 
-    R = NORM(pos);
-
-    SUB(pmt_dir,pmts[pars->i].pos,pos);
-
-    distance = NORM(pmt_dir);
+    get_path_length(pos,pmts[pars->i].pos,AV_RADIUS,&l_d2o,&l_h2o);
 
     /* FIXME: I just calculate delta assuming 400 nm light. */
     wavelength0 = 400.0;
-    if (R <= AV_RADIUS) {
-        n = get_index_snoman_d2o(wavelength0);
-    } else {
-        n = get_index_snoman_h2o(wavelength0);
-    }
+    n_d2o = get_index_snoman_d2o(wavelength0);
+    n_h2o = get_index_snoman_h2o(wavelength0);
 
-    t += distance*n/SPEED_OF_LIGHT;
+    t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
 
     return t*get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
 }
@@ -153,7 +146,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
      *
      * `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, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, l_h2o, l_d2o, wavelength0;
+    double pmt_dir[3], tmp[3], R, cos_theta, theta, x, z, s, a, b, beta, E, p, T, omega, theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, l_h2o, l_d2o, wavelength0;
 
     /* Calculate beta at the start of the track. */
     E0 = T0 + MUON_MASS;
@@ -175,8 +168,9 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
 
     /* Compute the Cerenkov angle at the start of the track. */
     wavelength0 = 400.0;
-    n = get_index_snoman_d2o(wavelength0);
-    theta_cerenkov = acos(1/(n*beta0));
+    n_d2o = get_index_snoman_d2o(wavelength0);
+    n_h2o = get_index_snoman_h2o(wavelength0);
+    theta_cerenkov = acos(1/(n_d2o*beta0));
 
     /* Now, we compute the distance along the track where the PMT is at the
      * Cerenkov angle.
@@ -205,8 +199,10 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
      * `s`. */
     a = NORM(tmp);
 
+    get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);
+
     /* Assume the particle is travelling at the speed of light. */
-    *t = s/SPEED_OF_LIGHT + a*n/SPEED_OF_LIGHT;
+    *t = s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
 
     /* `z` is the distance to the PMT projected onto the track direction. */
     z = R*cos_theta;
@@ -225,7 +221,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
     p = sqrt(E*E - MUON_MASS*MUON_MASS);
     beta = p/E;
 
-    if (beta < 1/n) return 0.0;
+    if (beta < 1/n_d2o) return 0.0;
 
     /* `prob` is the number of photons emitted per cm by the particle at a
      * distance `s` along the track. */
@@ -249,16 +245,14 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
 
     theta0 = fmax(theta0*sqrt(s),MIN_THETA0);
 
-    double frac = sqrt(2)*n*x*beta0*theta0;
+    double frac = sqrt(2)*n_d2o*x*beta0*theta0;
 
     f = get_weighted_pmt_response(acos(-cos_theta_pmt));
 
     absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0);
     absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0);
 
-    get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);
-
-    return f*exp(-l_d2o/absorption_length_d2o)*exp(-l_h2o/absorption_length_h2o)*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);
+    return f*exp(-l_d2o/absorption_length_d2o)*exp(-l_h2o/absorption_length_h2o)*n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI);
 }
 
 typedef struct betaRootParams {