From a3fab763491be8d8d6e4040895615431827c4e58 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Sun, 24 Nov 2019 13:48:41 -0600 Subject: update get_expected_charge to not skip calculating the charge if the angle is large Previously to achieve a large speedup in the likelihood calculation I added a line to skip calculating the charge if: abs((cos(theta)-cos_theta_cerenkov)/(sin_theta*theta0)) > 5 However I noticed that this was causing discontinuities in the likelihood function when fitting low energy muons so I'm putting it behind a compile time flag for now. --- src/likelihood.c | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/likelihood.c b/src/likelihood.c index 30a745c..f77cab6 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -607,7 +607,34 @@ static void get_expected_charge(double beta, double theta0, double *pos, double theta0 = fmax(theta0,get_theta0_min(distance_to_pmt,PMT_RADIUS,sin_theta_pmt)); +#ifdef FAST_GET_EXPECTED_CHARGE + /* This next line is used to skip out of calculating the expected charge if + * abs((cos(theta)-cos_theta_cherenkov)/(sin(theta)*theta0)) > 5. The idea + * here is that later in the likelihood calculation the PDF for Cerenkov + * light has a term like + * + * exp(((cos(theta)-cos_theta_cherenkov)/(sin(theta)*theta0))**2) + * + * which will be less than 10^-25 if the following inequality holds and + * therefore the charge should be negligible. + * + * However! I noticed that this line was causing discontinuities in the + * likelihood when fitting low energy muons. I realized there are two + * potential issues with this. One is that the PDF is multiplied by two + * other terms: 1/(sin_theta*theta0) and the solid angle of the PMT. + * Therefore if these are large it may have an impact. Second, the PDF + * actually integrates over the Cerenkov spectrum and correctly takes into + * account the change in the index as a function of wavelength. + * + * For now we don't use this to prevent the discontinuities. However, it + * would be nice to use it in the future since it provides a massive ~2.5x + * speedup. One idea for trying to eventually include it again is to update + * the PDF for Cerenkov light to not include the wavelength dependence + * (which would make it slightly less accurate, but it may be worth the + * 2.5x speedup). */ + if (fabs(cos_theta-cos_theta_cerenkov) > 5*sin_theta*theta0) return; +#endif omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS); -- cgit