add guess_time() function to approximate the PMT hit time

This function is only used when the expected number of photons reaching a PMT
is *very* small. In this case, we still need to estimate the PMT hit time PDF
for indirect light which is modelled as a flat distribution starting at the
time where the PMT is most likely to be hit from direct light. Since we compute
the most likely time for a PMT to be hit from direct light by computing the
integral of the expected charge times the time and then dividing by the total
charge, when the total charge is very small this can introduce large errors.

Note that this code already existed but it was computed in the likelihood
function. This commit just moves it to its own function to make things look
nicer.

1 files changed, 50 insertions, 63 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 8116261..4e44b00 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -642,6 +642,54 @@ static int get_smax(particle *p, double beta_min, double range, double *smax)
     return status;
 }
 
+static double guess_time(double *pos, double *dir, double *pmt_pos, double smax, double n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov)
+{
+    /* Returns an approximate time at which a PMT is most likely to get hit
+     * from Cerenkov light.
+     *
+     * To do this, we try to compute the distance along the track where the PMT
+     * is at the Cerenkov angle. */
+    double pmt_dir[3], tmp[3];
+    double R, cos_theta, sin_theta, s, l_d2o, l_h2o;
+
+    /* First, we find the point along the track where the PMT is at the
+     * Cerenkov angle. */
+    SUB(pmt_dir,pmt_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);
+    sin_theta = sqrt(1-pow(cos_theta,2));
+
+    /* 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*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov;
+
+    /* Make sure that the point is somewhere along the track between 0 and
+     * `smax`. */
+    if (s < 0) s = 0.0;
+    else if (s > smax) s = smax;
+
+    /* Compute the position of the particle at a distance `s` along the track. */
+    tmp[0] = pos[0] + s*dir[0];
+    tmp[1] = pos[1] + s*dir[1];
+    tmp[2] = pos[2] + s*dir[2];
+
+    SUB(pmt_dir,pmt_pos,tmp);
+
+    get_path_length(tmp,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
+
+    /* Assume the particle is travelling at the speed of light. */
+    return s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+}
+
 static double getKineticEnergy(double x, void *p)
 {
     return particle_get_energy(x, (particle *) p);
@@ -664,10 +712,8 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
     static double logp[MAX_PE], nll[MAX_PMTS];
     double range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
     particle *p;
-    double pmt_dir[3], R, cos_theta, sin_theta, wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov, s, l_h2o, l_d2o;
+    double wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov;
     double logp_path;
-    double a, b;
-    double tmp[3];
     double mu_reflected;
     path *path;
 
@@ -764,53 +810,6 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
             continue;
         }
 
-        /* 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 only integrate 1 meter on both sides of the
-         * point along the track where the PMT is at the Cerenkov angle. */
-
-        /* 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);
-        sin_theta = sqrt(1-pow(cos_theta,2));
-
-        /* 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*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov;
-
-        /* Make sure that the point is somewhere along the track between 0 and
-         * `smax`. */
-        if (s < 0) s = 0.0;
-        else if (s > smax) s = smax;
-
-        /* Compute the endpoints of the integration. */
-        a = s - 100.0;
-        b = s + 100.0;
-
-        if (a < 0.0) a = 0.0;
-        if (b > range) b = range;
-
         double q_indirect;
         integrate_path(path, i, mu_direct+i, &q_indirect, &result);
         mu_indirect += q_indirect;
@@ -823,19 +822,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
              * in the denominator. Therefore, we estimate the time using
              * the same technique as in get_total_charge_approx(). See that
              * function for more details. */
-            n_h2o = get_index_snoman_h2o(wavelength0);
-
-            /* Compute the position of the particle at a distance `s` along the track. */
-            tmp[0] = pos[0] + s*dir[0];
-            tmp[1] = pos[1] + s*dir[1];
-            tmp[2] = pos[2] + s*dir[2];
-
-            SUB(pmt_dir,pmts[i].pos,tmp);
-
-            get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);
-
-            /* Assume the particle is travelling at the speed of light. */
-            ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+            ts[i] = t0 + guess_time(pos,dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov);
         }
 
         if (id == IDP_E_MINUS || id == IDP_E_PLUS) {