From e8921bc1865c6a03900408c0ef823fb2f5c9138c Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 18 Nov 2019 20:11:35 -0600 Subject: 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. --- src/likelihood.c | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) (limited to 'src/likelihood.c') 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); -- cgit