1 files changed, 231 insertions, 66 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 012243c..fc1cf54 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -13,6 +13,10 @@
 #include "pdg.h"
 #include "path.h"
 #include <stddef.h> /* for size_t */
+#include "scattering.h"
+#include "solid_angle.h"
+#include <gsl/gsl_roots.h>
+#include <gsl/gsl_errno.h>
 
 typedef struct intParams {
     path *p;
@@ -117,19 +121,208 @@ static double gsl_muon_charge(double x, void *params)
     return get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
 }
 
+double get_total_charge_approx(double T0, double *pos, double *dir, int i, double smax, double theta0, double *t)
+{
+    /* Returns the approximate expected number of photons seen by PMT `i` using
+     * an analytic formula.
+     *
+     * To come up with an analytic formula for the expected number of photons,
+     * it was necessary to make the following approximations:
+     *
+     * - the index of refraction is constant
+     * - the particle track is a straight line
+     * - the integral along the particle track is dominated by the gaussian
+     *   term describing the angular distribution of the light
+     *
+     * With these approximations and a few other ones (like using a Taylor
+     * expansion for the distance to the PMT), it is possible to pull
+     * everything out of the track integral and assume it's equal to it's value
+     * along the track where the exponent of the gaussian dominates.
+     *
+     * The point along the track where the exponent dominates is calculated by
+     * finding the point along the track where the angle between the track
+     * direction and the PMT is equal to the Cerenkov angle. If this point is
+     * before the start of the track, we use the start of the track and if it's
+     * past the end of `smax` we use `smax`.
+     *
+     * Since the integral over the track also contains a term like
+     * (1-1/(beta**2*n**2)) which is not constant near the end of the track, it
+     * is necessary to define `smax` as the point along the track where the
+     * particle velocity drops below some threshold.
+     *
+     * `smax` is currently calculated as the point where the particle velocity
+     * drops to 0.8 times the speed of light. */
+    double pmt_dir[3], tmp[3], R, cos_theta, theta, x, z, s, a, b, beta, E, p, T, omega, theta_cerenkov, n, sin_theta, E0, p0, beta0;
+
+    /* Calculate beta at the start of the track. */
+    E0 = T0 + MUON_MASS;
+    p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
+    beta0 = p0/E0;
+
+    /* 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. */
+    n = get_index_snoman_d2o(400.0);
+    theta_cerenkov = acos(1/(n*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);
+
+    /* 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 vector from the point `s` along the track to the PMT. */
+    tmp[0] = pmts[i].pos[0] - (pos[0] + s*dir[0]);
+    tmp[1] = pmts[i].pos[1] - (pos[1] + s*dir[1]);
+    tmp[2] = pmts[i].pos[2] - (pos[2] + s*dir[2]);
+
+    /* To do the integral analytically, we expand the distance to the PMT along
+     * the track in a Taylor series around `s0`, i.e.
+     *
+     *     r(s) = a + b*(s-s0)
+     *
+     * Here, we calculate `a` which is the distance to the PMT at the point
+     * `s`. */
+    a = NORM(tmp);
+
+    /* Assume the particle is travelling at the speed of light. */
+    *t = s/SPEED_OF_LIGHT + a*n/SPEED_OF_LIGHT;
+
+    /* `z` is the distance to the PMT projected onto the track direction. */
+    z = R*cos_theta;
+
+    /* `x` is the perpendicular distance from the PMT position to the track. */
+    x = R*fabs(sin(acos(cos_theta)));
+
+    /* `b` is the second coefficient in the Taylor expansion. */
+    b = (s-z)/a;
+
+    /* Compute the kinetic energy at the point `s` along the track. */
+    T = get_T(T0,s,HEAVY_WATER_DENSITY);
+
+    /* Calculate the particle velocity at the point `s`. */
+    E = T + MUON_MASS;
+    p = sqrt(E*E - MUON_MASS*MUON_MASS);
+    beta = p/E;
+
+    if (beta < 1/n) return 0.0;
+
+    /* `prob` is the number of photons emitted per cm by the particle at a
+     * distance `s` along the track. */
+    double prob = get_probability2(beta);
+
+    /* 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);
+
+    /* Calculate the sine of the angle between the track direction and the PMT
+     * at the position `s` along the track. */
+    sin_theta = fabs(sin(acos(DOT(dir,pmt_dir)/NORM(pmt_dir))));
+
+    /* Get the solid angle of the PMT at the position `s` along the track. */
+    omega = get_solid_angle_approx(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS);
+
+    theta0 = fmax(theta0*sqrt(s),MIN_THETA0);
+
+    double frac = sqrt(2)*n*x*beta0*theta0;
+
+    return n*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n*(smax-z)*beta0)/frac) + erf((-a+b*s+n*z*beta0)/frac))/(b+n*beta0)/(4*M_PI);
+}
+
 double getKineticEnergy(double x, double T0)
 {
     return get_T(T0, x, HEAVY_WATER_DENSITY);
 }
 
-double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel)
+static double beta_root(double x, void *params)
+{
+    /* Function used to find at what point along a track a particle crosses
+     * some threshold in beta. */
+    double T, E, p, beta, T0, beta_min;
+
+    T0 = ((double *) params)[0];
+    beta_min = ((double *) params)[1];
+
+    T = get_T(T0, x, HEAVY_WATER_DENSITY);
+
+    /* Calculate total energy */
+    E = T + MUON_MASS;
+    p = sqrt(E*E - MUON_MASS*MUON_MASS);
+    beta = p/E;
+
+    return beta - beta_min;
+}
+
+static int get_smax(double T0, double beta_min, double range, double *smax)
+{
+    /* Find the point along the track at which the particle's velocity drops to
+     * `beta_min`. */
+    int status;
+    double params[2];
+    gsl_root_fsolver *s;
+    gsl_function F;
+    int iter = 0, max_iter = 100;
+    double r, x_lo, x_hi;
+
+    s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
+
+    params[0] = T0;
+    params[1] = beta_min;
+
+    F.function = &beta_root;
+    F.params = params;
+
+    gsl_root_fsolver_set(s, &F, 0.0, range);
+
+    do {
+        iter++;
+        status = gsl_root_fsolver_iterate(s);
+        r = gsl_root_fsolver_root(s);
+        x_lo = gsl_root_fsolver_x_lower(s);
+        x_hi = gsl_root_fsolver_x_upper(s);
+
+        /* Find the root to within 1 mm. */
+        status = gsl_root_test_interval(x_lo, x_hi, 1e-1, 0);
+
+        if (status == GSL_SUCCESS) break;
+    } while (status == GSL_CONTINUE && iter < max_iter);
+
+    gsl_root_fsolver_free(s);
+
+    *smax = r;
+
+    return status;
+}
+
+double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast)
 {
     size_t i, j, nhit;
     intParams params;
     double total_charge;
-    double logp[MAX_PE], nll[MAX_PMTS], range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov, theta0, E0, p0, beta0;
+    double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu;
     double tmean = 0.0;
-    int npmt = 0;
+    int jmax;
 
     double mu_direct[MAX_PMTS];
     double ts[MAX_PMTS];
@@ -156,75 +349,42 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
 
     params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, z1, z2, n, MUON_MASS);
 
+    if (beta0 > BETA_MIN)
+        get_smax(T0, BETA_MIN, range, &smax);
+    else
+        smax = 0.0;
+
     total_charge = 0.0;
-    npmt = 0;
     for (i = 0; i < MAX_PMTS; i++) {
         if (ev->pmt_hits[i].flags || (pmts[i].pmt_type != PMT_NORMAL && pmts[i].pmt_type != PMT_OWL)) continue;
 
         params.i = i;
 
-        /* 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. */
-
-        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. */
-        cos_theta = DOT(dir,pmt_dir);
-        /* Compute the angle between the track direction and the PMT. */
-        theta = acos(cos_theta);
-        /* Compute the Cerenkov angle. Note that this isn't entirely correct
-         * since we aren't including the factor of beta, but since the point is
-         * just to split up the integral, we only need to find a point along
-         * the track close enough such that the integral isn't completely zero.
-         */
-        theta_cerenkov = acos(1/get_index_snoman_d2o(400.0));
-
-        /* Now, we compute the distance along the track where the PMT is at the
-         * Cerenkov angle. */
-        x = 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;
-
-        total_charge += mu_direct[i];
-
-        if (mu_direct[i] > 0.001) {
-            F.function = &gsl_muon_time;
-            gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals);
-            ts[i] = result;
-
-            ts[i] /= mu_direct[i];
+        if (fast) {
+            mu_direct[i] = get_total_charge_approx(T0, pos, dir, i, smax, theta0, &ts[i]);
             ts[i] += t0;
-            tmean += ts[i];
-            npmt += 1;
         } else {
-            ts[i] = 0.0;
+            F.function = &gsl_muon_charge;
+            gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals);
+            mu_direct[i] = result;
+
+            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];
+            }
         }
+
+        tmean += ts[i]*mu_direct[i];
+
+        total_charge += mu_direct[i];
     }
 
     path_free(params.p);
 
-    if (npmt)
-        tmean /= npmt;
+    if (total_charge > 0)
+        tmean /= total_charge;
 
     gsl_integration_cquad_workspace_free(w);
 
@@ -240,18 +400,23 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
     for (i = 0; i < MAX_PMTS; i++) {
         if (ev->pmt_hits[i].flags || (pmts[i].pmt_type != PMT_NORMAL && pmts[i].pmt_type != PMT_OWL)) continue;
 
+        log_mu = log(mu[i]);
+
         if (ev->pmt_hits[i].hit) {
-            logp[0] = -INFINITY;
-            for (j = 1; j < MAX_PE; j++) {
-                logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log(mu[i]) - gsl_sf_lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], tmean, 1.5);
+            jmax = (int) ceil(fmax(mu[i],1)*STD_MAX + fmax(mu[i],1));
+
+            if (jmax > MAX_PE) jmax = MAX_PE;
+
+            for (j = 1; j < jmax; j++) {
+                logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], tmean, 1.5);
             }
-            nll[nhit++] = -logsumexp(logp, sizeof(logp)/sizeof(double));
+            nll[nhit++] = -logsumexp(logp+1, jmax-1);
         } else {
             logp[0] = -mu[i];
-            for (j = 1; j < MAX_PE; j++) {
-                logp[j] = log(get_pmiss(j)) - mu[i] + j*log(mu[i]) - gsl_sf_lnfact(j);
+            for (j = 1; j < MAX_PE_NO_HIT; j++) {
+                logp[j] = get_log_pmiss(j) - mu[i] + j*log_mu - lnfact(j);
             }
-            nll[nhit++] = -logsumexp(logp, sizeof(logp)/sizeof(double));
+            nll[nhit++] = -logsumexp(logp, MAX_PE_NO_HIT);
         }
     }