diff options
| -rw-r--r-- | src/likelihood.c | 14 |
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); |