add a fast likelihood function

This commit adds a fast function to calculate the expected number of PE at a
PMT without numerically integrating over the track. This calculation is *much*
faster than integrating over the track (~30 ms compared to several seconds) and
so we use it during the "quick" minimization phase of the fit to quickly find
the best position.

10 files changed, 466 insertions, 78 deletions

diff --git a/src/fit.c b/src/fit.c
index 50501b2..4635ccf 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -26,6 +26,7 @@ static size_t iter;
 typedef struct fitParams {
     event *ev;
     double epsrel;
+    int fast;
 } fitParams;
 
 /* In order to start the fitter close to the minimum, we first do a series of
@@ -94,7 +95,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
     z2[0] = x[8];
 
     gettimeofday(&tv_start, NULL);
-    fval = nll_muon(fpars->ev, T, pos, dir, t0, z1, z2, 1, fpars->epsrel);
+    fval = nll_muon(fpars->ev, T, pos, dir, t0, z1, z2, 1, fpars->epsrel, fpars->fast);
     gettimeofday(&tv_stop, NULL);
 
     long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
@@ -226,7 +227,8 @@ int fit_event(event *ev, double *xopt, double *fmin)
 
     /* If our guess is below the Cerenkov threshold, start at the Cerenkov
      * threshold. */
-    if (T0 < Tmin) T0 = Tmin;
+    double Tmin2 = sqrt(MUON_MASS*MUON_MASS/(1-pow(0.9,2)))-MUON_MASS;
+    if (T0 < Tmin2) T0 = Tmin2;
 
     x0[0] = T0;
     x0[1] = 0.0;
@@ -283,6 +285,7 @@ int fit_event(event *ev, double *xopt, double *fmin)
      * maximum number of function evaluations to 100, and the relative
      * tolerance on the numerical integration to 10%. */
     fpars.epsrel = 1e-1;
+    fpars.fast = 1;
     nlopt_set_ftol_abs(opt, 1.0);
     nlopt_set_maxeval(opt, 100);
 
@@ -306,12 +309,16 @@ int fit_event(event *ev, double *xopt, double *fmin)
 
     memcpy(x,xopt,sizeof(x));
 
+    /* Reset path coefficients. */
+    x[7] = 0.0;
+    x[8] = 0.0;
+
     /* Now, we do the "real" minimization. */
     fpars.epsrel = 1e-4;
+    fpars.fast = 0;
     nlopt_set_ftol_abs(opt, 1e-5);
     nlopt_set_maxeval(opt, 1000);
 
-
     do {
         *fmin = fval;
         memcpy(xopt,x,sizeof(x));
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);
         }
     }
 
diff --git a/src/likelihood.h b/src/likelihood.h
index 4871f42..8ddaf2a 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -4,14 +4,22 @@
 #include "event.h"
 #include <stddef.h> /* for size_t */
 
+/* Maximum number of PE to consider if the PMT was hit. */
 #define MAX_PE 100
+/* Maximum number of PE to consider if the PMT wasn't hit.
+ *
+ * Note: This must be less than MAX_PE. */
+#define MAX_PE_NO_HIT 10
+#define STD_MAX 10
+
 #define CHARGE_FRACTION 60000.0
 #define DARK_RATE 500.0
 /* Single PE transit time spread (ns). */
 #define PMT_TTS 1.5
 /* Event window (ns) */
 #define GTVALID 400.0
+#define BETA_MIN 0.8
 
-double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel);
+double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast);
 
 #endif
diff --git a/src/misc.c b/src/misc.c
index bff0ba9..82d114b 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -1,6 +1,126 @@
 #include "misc.h"
 #include <math.h>
 #include <stdlib.h> /* for size_t */
+#include <gsl/gsl_sf_gamma.h>
+
+static struct {
+    int n;
+    double f;
+} ln_fact_table[LNFACT_MAX + 1] = {
+    {0,0},
+    {1,0},
+    {2,0.69314718055994529},
+    {3,1.791759469228055},
+    {4,3.1780538303479458},
+    {5,4.7874917427820458},
+    {6,6.5792512120101012},
+    {7,8.5251613610654147},
+    {8,10.604602902745251},
+    {9,12.801827480081469},
+    {10,15.104412573075516},
+    {11,17.502307845873887},
+    {12,19.987214495661885},
+    {13,22.552163853123425},
+    {14,25.19122118273868},
+    {15,27.89927138384089},
+    {16,30.671860106080672},
+    {17,33.505073450136891},
+    {18,36.395445208033053},
+    {19,39.339884187199495},
+    {20,42.335616460753485},
+    {21,45.380138898476908},
+    {22,48.471181351835227},
+    {23,51.606675567764377},
+    {24,54.784729398112319},
+    {25,58.003605222980518},
+    {26,61.261701761002001},
+    {27,64.557538627006338},
+    {28,67.88974313718154},
+    {29,71.257038967168015},
+    {30,74.658236348830158},
+    {31,78.092223553315307},
+    {32,81.557959456115043},
+    {33,85.054467017581516},
+    {34,88.580827542197682},
+    {35,92.136175603687093},
+    {36,95.719694542143202},
+    {37,99.330612454787428},
+    {38,102.96819861451381},
+    {39,106.63176026064346},
+    {40,110.32063971475739},
+    {41,114.03421178146171},
+    {42,117.77188139974507},
+    {43,121.53308151543864},
+    {44,125.3172711493569},
+    {45,129.12393363912722},
+    {46,132.95257503561632},
+    {47,136.80272263732635},
+    {48,140.67392364823425},
+    {49,144.5657439463449},
+    {50,148.47776695177302},
+    {51,152.40959258449735},
+    {52,156.3608363030788},
+    {53,160.3311282166309},
+    {54,164.32011226319517},
+    {55,168.32744544842765},
+    {56,172.35279713916279},
+    {57,176.39584840699735},
+    {58,180.45629141754378},
+    {59,184.53382886144948},
+    {60,188.6281734236716},
+    {61,192.7390472878449},
+    {62,196.86618167289001},
+    {63,201.00931639928152},
+    {64,205.1681994826412},
+    {65,209.34258675253685},
+    {66,213.53224149456327},
+    {67,217.73693411395422},
+    {68,221.95644181913033},
+    {69,226.1905483237276},
+    {70,230.43904356577696},
+    {71,234.70172344281826},
+    {72,238.97838956183432},
+    {73,243.26884900298271},
+    {74,247.57291409618688},
+    {75,251.89040220972319},
+    {76,256.22113555000954},
+    {77,260.56494097186322},
+    {78,264.92164979855278},
+    {79,269.29109765101981},
+    {80,273.67312428569369},
+    {81,278.06757344036612},
+    {82,282.4742926876304},
+    {83,286.89313329542699},
+    {84,291.32395009427029},
+    {85,295.76660135076065},
+    {86,300.22094864701415},
+    {87,304.68685676566872},
+    {88,309.1641935801469},
+    {89,313.65282994987905},
+    {90,318.1526396202093},
+    {91,322.66349912672615},
+    {92,327.1852877037752},
+    {93,331.71788719692847},
+    {94,336.26118197919845},
+    {95,340.81505887079902},
+    {96,345.37940706226686},
+    {97,349.95411804077025},
+    {98,354.53908551944079},
+    {99,359.1342053695754},
+    {100,363.73937555556347},
+};
+
+double lnfact(unsigned int n)
+{
+    /* Returns the logarithm of n!.
+     *
+     * Uses a lookup table to return results for n < 100. */
+
+    if (n <= LNFACT_MAX)
+        return ln_fact_table[n].f;
+
+    return gsl_sf_lnfact(n);
+}
 
 double kahan_sum(double *x, size_t n)
 {
diff --git a/src/misc.h b/src/misc.h
index c93afba..eb62e46 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -3,6 +3,9 @@
 
 #include <stdlib.h> /* for size_t */
 
+#define LNFACT_MAX 100
+
+double lnfact(unsigned int n);
 double kahan_sum(double *x, size_t n);
 double interp1d(double x, double *xp, double *yp, size_t n);
 int isclose(double a, double b, double rel_tol, double abs_tol);
diff --git a/src/scattering.c b/src/scattering.c
index 66e8398..433387d 100644
--- a/src/scattering.c
+++ b/src/scattering.c
@@ -10,6 +10,9 @@
 #include <stdlib.h> /* for exit() */
 #include <stdio.h> /* for fprintf() */
 #include <gsl/gsl_errno.h> /* for gsl_strerror() */
+#include "pdg.h"
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_spline.h>
 
 static double xlo = -1.0;
 static double xhi = 1.0;
@@ -26,6 +29,16 @@ static gsl_spline2d *spline;
 static gsl_interp_accel *xacc;
 static gsl_interp_accel *yacc;
 
+static double x2lo = 0.0;
+static double x2hi = 1.0;
+static size_t nx2 = 1000;
+
+static double *x2;
+static double *y2;
+
+static gsl_spline *spline2;
+static gsl_interp_accel *xacc2;
+
 static double prob_scatter(double wavelength, void *params)
 {
     /* Calculate the number of photons emitted per unit solid angle per cm at
@@ -45,6 +58,21 @@ static double prob_scatter(double wavelength, void *params)
     return qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI);
 }
 
+static double prob_scatter2(double wavelength, void *params)
+{
+    /* Calculate the number of photons emitted per per cm for a particle
+     * travelling with velocity `beta`. */
+    double qe, index;
+
+    double beta = ((double *) params)[0];
+
+    qe = get_quantum_efficiency(wavelength);
+
+    index = get_index_snoman_d2o(wavelength);
+
+    return 2*M_PI*FINE_STRUCTURE_CONSTANT*(1-(1/(beta*beta*index*index)))*qe/pow(wavelength,2)*1e7;
+}
+
 void init_interpolation(void)
 {
     size_t i, j;
@@ -52,6 +80,8 @@ void init_interpolation(void)
     double result, error;
     size_t nevals;
     int status;
+    gsl_integration_cquad_workspace *w;
+    gsl_function F;
 
     x = malloc(nx*sizeof(double));
     y = malloc(ny*sizeof(double));
@@ -69,9 +99,8 @@ void init_interpolation(void)
         y[i] = ylo + (yhi-ylo)*i/(ny-1);
     }
 
-    gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(100);
+    w = gsl_integration_cquad_workspace_alloc(100);
 
-    gsl_function F;
     F.function = &prob_scatter;
     F.params = params;
 
@@ -92,6 +121,34 @@ void init_interpolation(void)
 
     gsl_spline2d_init(spline, x, y, z, nx, ny);
 
+    x2 = malloc(nx2*sizeof(double));
+    y2 = malloc(nx2*sizeof(double));
+
+    spline2 = gsl_spline_alloc(gsl_interp_linear, nx2);
+    xacc2 = gsl_interp_accel_alloc();
+
+    for (i = 0; i < nx2; i++) {
+        x2[i] = x2lo + (x2hi-x2lo)*i/(nx2-1);
+    }
+
+    F.function = &prob_scatter2;
+    F.params = params;
+
+    for (i = 0; i < nx2; i++) {
+        params[0] = x2[i];
+
+        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+        if (status) {
+            fprintf(stderr, "error integrating photon angular distribution: %s\n", gsl_strerror(status));
+            exit(1);
+        }
+
+        y2[i] = result;
+    }
+
+    gsl_spline_init(spline2, x2, y2, nx2);
+
     gsl_integration_cquad_workspace_free(w);
 }
 
@@ -106,6 +163,11 @@ double get_probability(double beta, double cos_theta, double theta0)
     return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/(theta0*sin_theta);
 }
 
+double get_probability2(double beta)
+{
+    return gsl_spline_eval(spline2, beta, xacc2);
+}
+
 void free_interpolation(void)
 {
     free(x);
diff --git a/src/scattering.h b/src/scattering.h
index f3a38d8..308f8ef 100644
--- a/src/scattering.h
+++ b/src/scattering.h
@@ -3,6 +3,7 @@
 
 void init_interpolation(void);
 double get_probability(double beta, double cos_theta, double theta0);
+double get_probability2(double beta);
 void free_interpolation(void);
 
 #endif
diff --git a/src/sno_charge.c b/src/sno_charge.c
index fe1e553..31de608 100644
--- a/src/sno_charge.c
+++ b/src/sno_charge.c
@@ -101,7 +101,7 @@ double pq(double q, int n)
     return gsl_spline_eval(splines[n-1], q, acc);
 }
 
-double get_pmiss(int n)
+double get_log_pmiss(int n)
 {
     if (!initialized) {
         fprintf(stderr, "charge interpolation hasn't been initialized!\n");
@@ -109,10 +109,9 @@ double get_pmiss(int n)
     }
 
     if (n == 0) {
-        return 1.0;
-    } else if (n > MAX_PE) {
-        /* Assume the distribution is gaussian by the central limit theorem. */
         return 0.0;
+    } else if (n > MAX_PE) {
+        return -INFINITY;
     }
 
     return pmiss[n-1];
@@ -223,7 +222,7 @@ void init_charge(void)
     for (i = 1; i <= MAX_PE; i++) {
         params[0] = i;
         gsl_integration_cquad(&F, 0, MEAN_THRESH/MEAN_HIPT, 0, 1e-9, w, &result, &error, &nevals);
-        pmiss[i-1] = result;
+        pmiss[i-1] = log(result);
     }
 
     F.function = &gsl_smear;
diff --git a/src/sno_charge.h b/src/sno_charge.h
index 545fdd2..bd95cc4 100644
--- a/src/sno_charge.h
+++ b/src/sno_charge.h
@@ -3,7 +3,7 @@
 
 void init_charge(void);
 double pq(double q, int n);
-double get_pmiss(int n);
+double get_log_pmiss(int n);
 double spe_pol2exp(double q);
 
 #endif
diff --git a/src/test.c b/src/test.c
index fb788ba..fa4256f 100644
--- a/src/test.c
+++ b/src/test.c
@@ -467,6 +467,28 @@ err:
     return 1;
 }
 
+int test_lnfact(char *err)
+{
+    /* Tests the lnfact() function. */
+    size_t i;
+    double x, expected;
+
+    for (i = 0; i < 1000; i++) {
+        x = lnfact(i);
+        expected = gsl_sf_lnfact(i);
+
+        if (!isclose(x, expected, 1e-9, 1e-9)) {
+            sprintf(err, "lnfact(%zu) returned %.5g, but expected %.5g", i, x, expected);
+            goto err;
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -482,6 +504,7 @@ struct tests {
     {test_path, "test_path"},
     {test_interp1d, "test_interp1d"},
     {test_kahan_sum, "test_kahan_sum"},
+    {test_lnfact, "test_lnfact"},
 };
 
 int main(int argc, char **argv)