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.

1 files changed, 27 insertions, 0 deletions

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);