split up the track integral into two intervals

This commit updates the likelihood calculation to split up the track integral
into two intervals in some cases. I noticed when fitting some events that the
likelihood value would change drastically for a very small change in the fit
parameters. I eventually tracked it down to the fact that the track integral
was occasionally returning a very small charge for a PMT which should have a
very high charge. This was happening because the region of the track which was
hitting the PMT was very small and the cquad integration routine was completely
skipping it.

The solution to this problem is a bit of a hack, but it seems to work. I first
calculate where along the track (for a straight track) the PMT would be at the
Cerenkov angle from the track. If this point is somewhere along the track then
we split up the integral into two intervals: one going from the start of the
track to this point and the other from the point to the end of the track. Since
the cquad routine always samples points near the end of the intervals this
should prevent it from completely skipping over the point in the track where
the integrand is non-zero.

1 files changed, 66 insertions, 5 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index db59040..e045512 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -333,6 +333,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
     double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
     double tmean = 0.0;
     muon_energy *m;
+    double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, theta_cerenkov, s;
 
     double mu_direct[MAX_PMTS];
     double ts[MAX_PMTS];
@@ -376,15 +377,75 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
             mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i]);
             ts[i] += t0;
         } else {
+            /* First, we try to compute the distance along the track where the
+             * PMT is at the Cerenkov angle. The reason for this is because for
+             * heavy particles like muons which don't scatter much, the
+             * probability distribution for getting a photon hit along the
+             * track looks kind of like a delta function, i.e. the PMT is only
+             * hit over a very narrow window when the angle between the track
+             * direction and the PMT is *very* close to the Cerenkov angle
+             * (it's not a perfect delta function since there is some width due
+             * to dispersion). In this case, it's possible that the numerical
+             * integration completely skips over the delta function and so
+             * predicts an expected charge of 0. To fix this, we compute the
+             * integral in two steps, one up to the point along the track where
+             * the PMT is at the Cerenkov angle and another from that point to
+             * the end of the track. Since the integration routine always
+             * samples points near the beginning and end of the integral, this
+             * allows the routine to correctly compute that the integral is non
+             * zero. */
+
+            /* First, we find the point along the track where the PMT is at the
+             * Cerenkov angle. */
+            SUB(pmt_dir,pmts[i].pos,pos);
+            /* Compute the distance to the PMT. */
+            R = NORM(pmt_dir);
+            normalize(pmt_dir);
+
+            /* Calculate the cosine of the angle between the track direction and the
+             * vector to the PMT at the start of the track. */
+            cos_theta = DOT(dir,pmt_dir);
+            /* Compute the angle between the track direction and the PMT. */
+            theta = acos(cos_theta);
+
+            /* Compute the Cerenkov angle at the start of the track. */
+            wavelength0 = 400.0;
+            n_d2o = get_index_snoman_d2o(wavelength0);
+            theta_cerenkov = acos(1/(n_d2o*beta0));
+
+            /* Now, we compute the distance along the track where the PMT is at the
+             * Cerenkov angle.
+             *
+             * Note: This formula comes from using the "Law of sines" where the three
+             * vertices of the triangle are the starting position of the track, the
+             * point along the track that we want to find, and the PMT position. */
+            s = R*sin(theta_cerenkov-theta)/sin(theta_cerenkov);
+
             F.function = &gsl_muon_charge;
-            gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals);
-            mu_direct[i] = result;
+
+            if (s > 0 && s < range) {
+                gsl_integration_cquad(&F, 0, s, 0, epsrel, w, &result, &error, &nevals);
+                mu_direct[i] = result;
+                gsl_integration_cquad(&F, s, range, 0, epsrel, w, &result, &error, &nevals);
+                mu_direct[i] += result;
+            } else {
+                gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals);
+                mu_direct[i] = result;
+            }
+
+            F.function = &gsl_muon_time;
 
             ts[i] = t0;
             if (mu_direct[i] > 1e-9) {
-                F.function = &gsl_muon_time;
-                gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals);
-                ts[i] += result/mu_direct[i];
+                if (s > 0 && s < range) {
+                    gsl_integration_cquad(&F, 0, s, 0, epsrel, w, &result, &error, &nevals);
+                    ts[i] += result/mu_direct[i];
+                    gsl_integration_cquad(&F, s, range, 0, epsrel, w, &result, &error, &nevals);
+                    ts[i] += result/mu_direct[i];
+                } else {
+                    gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals);
+                    ts[i] += result/mu_direct[i];
+                }
             }
         }