update get_expected_charge() to always use the index of refraction for d2o

This commit updates get_expected_charge() to always use the index of refraction
for d2o instead of choosing the index of d2o or h2o based on the position of
the particle. The reason for this is that selecting the index based on the
position was causing discontinuities in the likelihood function for muon tracks
which crossed the AV.

1 files changed, 10 insertions, 4 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index bf4148c..30a745c 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -572,11 +572,17 @@ double get_theta0_min(double R, double r, double sin_theta_pmt)
 
 static void get_expected_charge(double beta, double theta0, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
 {
-    double pmt_dir[3], cos_theta, sin_theta, n, omega, R, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o, cos_theta_cerenkov, distance_to_pmt;
+    double pmt_dir[3], cos_theta, sin_theta, n, omega, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o, cos_theta_cerenkov, distance_to_pmt;
 
-    R = NORM(pos);
-
-    n = (R <= AV_RADIUS) ? get_avg_index_d2o() : get_avg_index_h2o();
+    /* Previously the index of refraction was calculated based on the position,
+     * i.e.:
+     *
+     *     n = (NORM(pos) <= AV_RADIUS) ? get_avg_index_d2o() : get_avg_index_h2o();
+     *
+     * but this was causing a discontinuity in the likelihood function when the
+     * particle track crossed the AV since we're numerically integrating over a
+     * finite set of points. */
+    n = get_avg_index_d2o();
 
     cos_theta_cerenkov = 1/(beta*n);