update method for calculating expected number of photons from shower and delta rays

This commit introduces a new method for integrating over the particle track to
calculate the number of shower and delta ray photons expected at each PMT. The
reason for introducing a new method was that the previous method of just using
the trapezoidal rule was both inaccurate and not stable. By inaccurate I mean
that the trapezoidal rule was not producing a very good estimate of the true
integral and by not stable I mean that small changes in the fit parameters
(like theta and phi) could produce wildly different results. This meant that
the likelihood function was very noisy and was causing the minimizers to not be
able to find the global minimum.

The new integration method works *much* better than the trapezoidal rule for
the specific functions we are dealing with. The problem is essentially to
integrate the product of two functions over some interval, one of which is very
"peaky", i.e. we want to find:

    \int f(x) g(x) dx

where f(x) is peaked around some region and g(x) is relatively smooth. For our
case, f(x) represents the angular distribution of the Cerenkov light and g(x)
represents the factors like solid angle, absorption, etc.

The technique I discovered was that you can approximate this integral via a
discrete sum:

    constant \sum_i g(x_i)

where the x_i are chosen to have equal spacing along the range of the integral
of f(x), i.e.

    x_i = F^(-1)(i*constant)

This new method produces likelihood functions which are *much* more smooth and
accurate than previously.

In addition, there are a few other fixes in this commit:

- switch from specifying a step size for the shower integration to a number of
  points, i.e. dx_shower -> number of shower points

- only integrate to the PSUP

    I realized that previously we were integrating to the end of the track even
    if the particle left the PSUP, and that there was no code to deal with the
    fact that light emitted beyond the PSUP can't make it back to the PMTs.

- only integrate to the Cerenkov threshold

    When integrating over the particle track to calculate the expected number
    of direct Cerenkov photons, we now only integrate the track up to the point
    where the particle's velocity is 1/index. This should hopefully make the
    likelihood smoother because previously the estimate would depend on exactly
    whether the points we sampled the track were above or below this point.

- add a minimum theta0 value based on the angular width of the PMT

    When calculating the expected number of Cerenkov photons we assumed that
    the angular distribution was constant over the whole PMT. This is a bad
    assumption when the particle is very close to the PMT. Really we should
    average the function over all the angles of the PMT, but that would be too
    computationally expensive so instead we just calculate a minimum theta0
    value which depends on the distance and angle to the PMT. This seems to
    make the likelihood much smoother for particles near the PSUP.

- add a factor of sin(theta) when checking if we can skip calculating the
  charge in get_expected_charge()

- fix a nan in beta_root() when the momentum is negative

- update PSUP_RADIUS from 800 cm -> 840 cm

7 files changed, 443 insertions, 239 deletions

diff --git a/src/fit.c b/src/fit.c
index 88a4a1f..951f614 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -66,7 +66,7 @@ static size_t iter;
 typedef struct fitParams {
     event *ev;
     double dx;
-    double dx_shower;
+    double ns;
     int fast;
     int print;
     int id[MAX_VERTICES];
@@ -5036,7 +5036,7 @@ static double nlopt_nll2(unsigned int n, const double *x, double *grad, void *pa
     }
 
     gettimeofday(&tv_start, NULL);
-    fval = nll(fpars->ev, v, fpars->n, fpars->dx, fpars->dx_shower, fpars->fast, fpars->charge_only, fpars->hit_only);
+    fval = nll(fpars->ev, v, fpars->n, fpars->dx, fpars->ns, fpars->fast, fpars->charge_only, fpars->hit_only);
     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;
@@ -5093,7 +5093,7 @@ static double nlopt_nll(unsigned int n, const double *x, double *grad, void *par
     v.n = 1;
 
     gettimeofday(&tv_start, NULL);
-    fval = nll(fpars->ev, &v, 1, fpars->dx, fpars->dx_shower, fpars->fast, fpars->charge_only, fpars->hit_only);
+    fval = nll(fpars->ev, &v, 1, fpars->dx, fpars->ns, fpars->fast, fpars->charge_only, fpars->hit_only);
     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;
@@ -5334,7 +5334,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
     /* First we do a set of "quick" minimizations to try and start the
      * minimizer close to the minimum. */
     fpars.dx = 1.0;
-    fpars.dx_shower = 10.0;
+    fpars.ns = 100;
     fpars.fast = 1;
     fpars.hit_only = 0;
     fpars.print = 0;
@@ -5453,7 +5453,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
 
     /* Now, we do the "real" minimization. */
     fpars.dx = 1.0;
-    fpars.dx_shower = 10.0;
+    fpars.ns = 100;
     fpars.fast = 0;
     fpars.print = 1;
     fpars.charge_only = 0;
@@ -5617,7 +5617,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id, double maxtime)
      * maximum number of function evaluations to 100, and the relative
      * tolerance on the numerical integration to 10%. */
     fpars.dx = 1.0;
-    fpars.dx_shower = 10.0;
+    fpars.ns = 100;
     fpars.fast = 1;
     fpars.print = 0;
     fpars.charge_only = 0;
@@ -5713,7 +5713,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id, double maxtime)
 
     /* Now, we do the "real" minimization. */
     fpars.dx = 1.0;
-    fpars.dx_shower = 10.0;
+    fpars.ns = 100;
     fpars.fast = 0;
     fpars.print = 1;
     fpars.charge_only = 0;
diff --git a/src/likelihood.c b/src/likelihood.c
index 6f79325..34f379a 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -39,6 +39,8 @@
 #include "id_particles.h"
 #include <gsl/gsl_cdf.h>
 #include "find_peaks.h"
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_spline.h>
 
 particle *particle_init(int id, double T0, size_t n)
 {
@@ -58,7 +60,7 @@ particle *particle_init(int id, double T0, size_t n)
      * which we solve by dividing the track up into `n` segments and then
      * numerically computing the energy at each step. */
     size_t i;
-    double dEdx, rad;
+    double dEdx, rad, norm, pdf, pdf_last;
 
     particle *p = malloc(sizeof(particle));
     p->id = id;
@@ -67,9 +69,18 @@ particle *particle_init(int id, double T0, size_t n)
     p->n = n;
     p->a = 0.0;
     p->b = 0.0;
+    p->cos_theta = calloc(n, sizeof(double));
+    p->cdf_shower = calloc(n, sizeof(double));
+    p->spline_shower = gsl_spline_alloc(gsl_interp_linear, n);
+    p->acc_shower = gsl_interp_accel_alloc();
     p->shower_photons = 0.0;
     p->delta_ray_a = 0.0;
     p->delta_ray_b = 0.0;
+    p->cdf_delta = calloc(n, sizeof(double));
+    p->x_shower = calloc(n, sizeof(double));
+    p->gamma_pdf = calloc(n, sizeof(double));
+    p->spline_delta = gsl_spline_alloc(gsl_interp_linear, n);
+    p->acc_delta = gsl_interp_accel_alloc();
     p->delta_ray_photons = 0.0;
 
     p->x[0] = 0;
@@ -86,6 +97,8 @@ particle *particle_init(int id, double T0, size_t n)
         p->a = electron_get_angular_distribution_alpha(T0);
         p->b = electron_get_angular_distribution_beta(T0);
 
+        electron_get_position_distribution_parameters(T0, &p->pos_a, &p->pos_b);
+
         p->delta_ray_photons = 0.0;
 
         /* Now we loop over the points along the track and assume dE/dx is
@@ -109,6 +122,34 @@ particle *particle_init(int id, double T0, size_t n)
 
         p->shower_photons = electron_get_shower_photons(T0, rad);
 
+        norm = electron_get_angular_pdf_norm(p->a, p->b, 1/avg_index_d2o);
+
+        for (i = 0; i < n; i++) {
+            p->cos_theta[i] = -1.0 + 2.0*i/(n-1);
+            /* Note: We assume here that the peak of the angular distribution
+             * is at the Cerenkov angle for a particle with beta = 1. This
+             * seems to be an OK assumption for high energy showers, but is not
+             * exactly correct for showers from lower energy electrons or for
+             * delta rays from lower energy muons. It seems good enough, but in
+             * the future it would be nice to parameterize this. */
+            pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/avg_index_d2o)/norm;
+            if (i > 0)
+                p->cdf_shower[i] = p->cdf_shower[i-1] + (pdf_last + pdf)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0;
+            else
+                p->cdf_shower[i] = 0.0;
+            pdf_last = pdf;
+        }
+
+        gsl_spline_init(p->spline_shower, p->cdf_shower, p->cos_theta, n);
+
+        p->xlo_shower = gsl_cdf_gamma_Pinv(0.01,p->pos_a,p->pos_b);
+        p->xhi_shower = gsl_cdf_gamma_Pinv(0.99,p->pos_a,p->pos_b);
+
+        for (i = 0; i < n; i++) {
+            p->x_shower[i] = p->xlo_shower + (p->xhi_shower-p->xlo_shower)*i/(n-1);
+            p->gamma_pdf[i] = gamma_pdf(p->x_shower[i],p->pos_a,p->pos_b);
+        }
+
         break;
     case IDP_MU_MINUS:
     case IDP_MU_PLUS:
@@ -123,6 +164,8 @@ particle *particle_init(int id, double T0, size_t n)
         p->a = muon_get_angular_distribution_alpha(T0);
         p->b = muon_get_angular_distribution_beta(T0);
 
+        muon_get_position_distribution_parameters(T0, &p->pos_a, &p->pos_b);
+
         muon_get_delta_ray_distribution_parameters(T0, &p->delta_ray_a, &p->delta_ray_b);
         p->delta_ray_photons = muon_get_delta_ray_photons(T0);
 
@@ -145,7 +188,29 @@ particle *particle_init(int id, double T0, size_t n)
          * range here using the dE/dx table instead of reading in the range. */
         p->T[n-1] = 0;
 
-        p->shower_photons = rad*MUON_PHOTONS_PER_MEV;
+        /* FIXME: add shower photons for muons. */
+        p->shower_photons = 0.0;
+
+        /* Calculate the normalization constant for the angular distribution. */
+        norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o);
+
+        for (i = 0; i < n; i++) {
+            p->cos_theta[i] = -1.0 + 2.0*i/(n-1);
+            /* Note: We assume here that the peak of the angular distribution
+             * is at the Cerenkov angle for a particle with beta = 1. This
+             * seems to be an OK assumption for high energy showers, but is not
+             * exactly correct for showers from lower energy electrons or for
+             * delta rays from lower energy muons. It seems good enough, but in
+             * the future it would be nice to parameterize this. */
+            pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o)/norm;
+            if (i > 0)
+                p->cdf_delta[i] = p->cdf_delta[i-1] + (pdf + pdf_last)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0;
+            else
+                p->cdf_delta[i] = 0.0;
+            pdf_last = pdf;
+        }
+
+        gsl_spline_init(p->spline_delta, p->cdf_delta, p->cos_theta, n);
 
         break;
     case IDP_PROTON:
@@ -154,6 +219,7 @@ particle *particle_init(int id, double T0, size_t n)
          * tables for heavy water for protons. */
         p->range = proton_get_range(T0, WATER_DENSITY);
 
+        /* FIXME: add delta ray photons for muons. */
         p->delta_ray_photons = 0.0;
 
         /* Now we loop over the points along the track and assume dE/dx is
@@ -175,6 +241,7 @@ particle *particle_init(int id, double T0, size_t n)
          * range here using the dE/dx table instead of reading in the range. */
         p->T[n-1] = 0;
 
+        /* FIXME: Is this correct? */
         p->shower_photons = rad*PROTON_PHOTONS_PER_MEV;
 
         break;
@@ -186,77 +253,243 @@ particle *particle_init(int id, double T0, size_t n)
     return p;
 }
 
-double particle_get_energy(double x, particle *p)
+/* Returns the cosine of the angle between a particle travelling in a straight
+ * line and a PMT after a distance `x` in cm.
+ *
+ * `R` is the distance between the particle and the PMT at the start of the
+ * path, and `cos_gamma` is the cosine of the angle between the particle
+ * direction and the PMT at the start of the path.  */
+double x_to_cos_theta(double x, double R, double cos_gamma)
 {
-    /* Returns the approximate kinetic energy of a particle in water after
-     * travelling `x` cm with an initial kinetic energy `T`.
-     *
-     * Return value is in MeV. */
-    double T;
+    return (R*cos_gamma-x)/sqrt(R*R+x*x-2*x*R*cos_gamma);
+}
 
-    T = interp1d(x,p->x,p->T,p->n);
+/* Returns the x position along a straight particle path given the cosine of
+ * the angle between the particle direction and a PMT.
+ *
+ * `R` is the distance between the particle and the PMT at the start of the
+ * path, `cos_gamma` is the cosine of the angle between the particle
+ * direction and the PMT at the start of the path, and `sin_gamma` is the sine
+ * of the angle (this is passed simply to avoid an extra call to sin()). */
+double cos_theta_to_x(double cos_theta, double R, double cos_gamma, double sin_gamma)
+{
+    double sin_theta;
 
-    if (T < 0) return 0;
+    /* If cos(theta) == 1.0 we must be at the beginning of the track. */
+    if (cos_theta == 1.0) return R*cos_gamma;
 
-    return T;
-}
+    /* We don't care about the sign here because theta is in the range [0,pi]
+     * and so sin(theta) is always positive. */
+    sin_theta = sqrt(1-cos_theta*cos_theta);
 
-void particle_free(particle *p)
-{
-    free(p->x);
-    free(p->T);
-    free(p);
+    /* This result comes from using the law of sines. */
+    return R*(cos_gamma - cos_theta*sin_gamma/sin_theta);
 }
 
-static void get_expected_charge_delta_ray(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
+/* Returns the x points and weights for numerically integrating the delta ray
+ * photons.
+ *
+ * The integral we want to solve here is:
+ *
+ *     \int x n(x) f(x) g(x)
+ *
+ * where n(x) is the number of delta ray photons emitted per cm, f(x) is the
+ * angular distribution of the delta ray photons, and g(x) contains all the
+ * other factors like solid angle, absorption, quantum efficiency, etc.
+ *
+ * We approximate this integral as:
+ *
+ *      \sum_i w_i g(x_i)
+ *
+ * where the weights and points are returned by this function.
+ *
+ * We currently assume the number of delta ray photons is constant along the
+ * particle track, i.e.
+ *
+ *     n(x) = 1/range
+ *
+ * This is an OK approximation, but in the future it might be nice to
+ * parameterize the number of photons produced as a function of x.
+ *
+ * `distance_to_psup` is the distance to the PSUP. If the distance to the PSUP
+ * is shorter than the range then we only integrate to the PSUP.
+ *
+ * Other arguments:
+ *
+ *     R - distance to PMT from start of track
+ *     cos_gamma - cosine of angle between track and PMT at start of track
+ *     sin_gamma - sine of angle between track and PMT at start of track
+ *
+ * Returns 1 if the integration can be skipped, and 0 otherwise. */
+int get_delta_ray_weights(particle *p, double range, double distance_to_psup, double *x, double *w, size_t n, double R, double cos_gamma, double sin_gamma)
 {
-    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o;
+    size_t i;
+    double x1, x2, cdf1, cdf2, cdf, delta, cos_theta, cos_theta_last, xcdf, xcdf_last;
 
-    SUB(pmt_dir,pmts[pmt].pos,pos);
+    /* First we need to figure out the bounds of integration in x. */
+    x1 = 0.0;
 
-    normalize(pmt_dir);
+    if (range < distance_to_psup)
+        x2 = range;
+    else
+        x2 = distance_to_psup;
+
+    /* Map bounds of integration in x -> cos(theta) and then find the values of
+     * the CDF at the endpoints. */
+    cdf1 = interp1d(x_to_cos_theta(x1,R,cos_gamma),p->cos_theta,p->cdf_delta,p->n);
+    cdf2 = interp1d(x_to_cos_theta(x2,R,cos_gamma),p->cos_theta,p->cdf_delta,p->n);
+
+    delta = (cdf2-cdf1)/n;
+
+    /* Now, we loop over different values of the CDF spaced evenly by delta and
+     * compute cos(theta) for each of these values, map that value back to x,
+     * and then compute the weight associated with that point. */
+    for (i = 0; i < n + 1; i++) {
+        cdf = cdf1 + i*delta;
+        cos_theta = gsl_spline_eval(p->spline_delta,cdf,p->acc_delta);
+        xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma);
+        if (i > 0) {
+            /* For each interval we approximate the integral as:
+             *
+             *     \int f(x) g(x) ~ <f(x)> \int g(x)
+             *
+             * and so we are free to choose best how to approximate the second
+             * integral. Since g(x) is assumed to be slowly varying as a
+             * function of x, it shouldn't matter too much, however based on
+             * some testing I found that it was best to use the midpoint, i.e.
+             *
+             *     \int g(x) ~ g((a+b)/2)
+             *
+             * as opposed to evaluating g(x) at the beginning, end, or
+             * averaging the function at the two points. */
+            x[i-1] = (xcdf + xcdf_last)/2;
+            w[i-1] = p->delta_ray_photons*delta*(1/range)*(xcdf - xcdf_last)/(cos_theta - cos_theta_last);
+        }
+        cos_theta_last = cos_theta;
+        xcdf_last = xcdf;
+    }
 
-    cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);
+    return 0;
+}
 
-    *t = 0.0;
-    *q = 0.0;
-    *reflected = 0.0;
-    if (cos_theta_pmt > 0) return;
+/* Returns the x points and weights for numerically integrating the shower
+ * photons.
+ *
+ * The integral we want to solve here is:
+ *
+ *     \int x n(x) f(x) g(x)
+ *
+ * where n(x) is the number of shower photons emitted per cm, f(x) is the
+ * angular distribution of the shower photons, and g(x) contains all the other
+ * factors like solid angle, absorption, quantum efficiency, etc.
+ *
+ * We approximate this integral as:
+ *
+ *      \sum_i w_i g(x_i)
+ *
+ * where the weights and points are returned by this function.
+ *
+ * `distance_to_psup` is the distance to the PSUP. If the distance to the PSUP is
+ * shorter than the shower range then we only integrate to the PSUP.
+ *
+ * Other arguments:
+ *
+ *     R - distance to PMT from start of track
+ *     cos_gamma - cosine of angle between track and PMT at start of track
+ *     sin_gamma - sine of angle between track and PMT at start of track
+ *
+ * Returns 1 if the integration can be skipped, and 0 otherwise. */
+int get_shower_weights(particle *p, double distance_to_psup, double *x, double *w, size_t n, double R, double cos_gamma, double sin_gamma)
+{
+    size_t i;
+    double x1, x2, cdf1, cdf2, cdf, delta, cos_theta, cos_theta_last, xcdf, xcdf_last;
 
-    /* Calculate the cosine of the angle between the track direction and the
-     * vector to the PMT. */
-    cos_theta = DOT(dir,pmt_dir);
+    /* First we need to figure out the bounds of integration in x. */
 
-    omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);
+    /* If the start of the shower is past the PSUP, no light will reach the
+     * PMTs, so we just return 1. */
+    if (p->xlo_shower > distance_to_psup)
+        return 1;
 
-    f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
-    f = get_weighted_pmt_response(-cos_theta_pmt);
+    x1 = p->xlo_shower;
 
-    get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
+    if (p->xhi_shower > distance_to_psup)
+        x2 = distance_to_psup;
+    else
+        x2 = p->xhi_shower;
+
+    /* Map bounds of integration in x -> cos(theta) and then find the values of
+     * the CDF at the endpoints. */
+    cdf1 = interp1d(x_to_cos_theta(x1,R,cos_gamma),p->cos_theta,p->cdf_shower,p->n);
+    cdf2 = interp1d(x_to_cos_theta(x2,R,cos_gamma),p->cos_theta,p->cdf_shower,p->n);
+
+    delta = (cdf2-cdf1)/n;
+
+    /* Now, we loop over different values of the CDF spaced evenly by delta and
+     * compute cos(theta) for each of these values, map that value back to x,
+     * and then compute the weight associated with that point. */
+    for (i = 0; i < n + 1; i++) {
+        cdf = cdf1 + i*delta;
+        cos_theta = gsl_spline_eval(p->spline_shower,cdf,p->acc_shower);
+        xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma);
+        if (i > 0) {
+            /* For each interval we approximate the integral as:
+             *
+             *     \int f(x) g(x) ~ <f(x)> \int g(x)
+             *
+             * and so we are free to choose best how to approximate the second
+             * integral. Since g(x) is assumed to be slowly varying as a
+             * function of x, it shouldn't matter too much, however based on
+             * some testing I found that it was best to use the midpoint, i.e.
+             *
+             *     \int g(x) ~ g((a+b)/2)
+             *
+             * as opposed to evaluating g(x) at the beginning, end, or
+             * averaging the function at the two points. */
+            x[i-1] = (xcdf + xcdf_last)/2;
+            w[i-1] = p->shower_photons*delta*interp1d(x[i-1],p->x_shower,p->gamma_pdf,p->n)*(xcdf-xcdf_last)/(cos_theta - cos_theta_last);
+        }
+        cos_theta_last = cos_theta;
+        xcdf_last = xcdf;
+    }
 
-    prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS);
-    prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS);
+    return 0;
+}
 
-    constant = get_weighted_quantum_efficiency()*omega/(2*M_PI);
+double particle_get_energy(double x, particle *p)
+{
+    /* Returns the approximate kinetic energy of a particle in water after
+     * travelling `x` cm with an initial kinetic energy `T`.
+     *
+     * Return value is in MeV. */
+    double T;
 
-    /* Note: We assume here that the peak of the angular distribution is at the
-     * Cerenkov angle for a particle with beta = 1. This seems to be an OK
-     * assumption for high energy showers, but is not exactly correct for
-     * showers from lower energy electrons or for delta rays from lower energy
-     * muons. It seems good enough, but in the future it would be nice to
-     * parameterize this. */
-    charge = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/avg_index_d2o);
+    T = interp1d(x,p->x,p->T,p->n);
 
-    *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
+    if (T < 0) return 0;
 
-    *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+    return T;
+}
 
-    *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
+void particle_free(particle *p)
+{
+    free(p->x);
+    free(p->T);
+    free(p->cos_theta);
+    free(p->cdf_shower);
+    free(p->x_shower);
+    free(p->gamma_pdf);
+    free(p->cdf_delta);
+    gsl_spline_free(p->spline_shower);
+    gsl_spline_free(p->spline_delta);
+    gsl_interp_accel_free(p->acc_shower);
+    gsl_interp_accel_free(p->acc_delta);
+    free(p);
 }
 
 static void get_expected_charge_shower(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
 {
-    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o;
+    double pmt_dir[3], omega, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o;
 
     SUB(pmt_dir,pmts[pmt].pos,pos);
 
@@ -264,15 +497,6 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
 
     cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);
 
-    *t = 0.0;
-    *q = 0.0;
-    *reflected = 0.0;
-    if (cos_theta_pmt > 0) return;
-
-    /* Calculate the cosine of the angle between the track direction and the
-     * vector to the PMT. */
-    cos_theta = DOT(dir,pmt_dir);
-
     omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);
 
     f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
@@ -283,15 +507,7 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
     prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS);
     prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS);
 
-    constant = get_weighted_quantum_efficiency()*omega/(2*M_PI);
-
-    /* Note: We assume here that the peak of the angular distribution is at the
-     * Cerenkov angle for a particle with beta = 1. This seems to be an OK
-     * assumption for high energy showers, but is not exactly correct for
-     * showers from lower energy electrons or for delta rays from lower energy
-     * muons. It seems good enough, but in the future it would be nice to
-     * parameterize this. */
-    charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/avg_index_d2o);
+    charge = get_weighted_quantum_efficiency()*omega/(2*M_PI);
 
     *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
 
@@ -300,9 +516,14 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
     *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
 }
 
+double get_theta0_min(double R, double r, double sin_theta_pmt)
+{
+    return fast_acos((R-r*sin_theta_pmt)/sqrt(r*r + R*R - 2*r*R*sin_theta_pmt));
+}
+
 static void get_expected_charge(double beta, double theta0, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
 {
-    double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o;
+    double pmt_dir[3], cos_theta, sin_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o;
 
     z = 1.0;
 
@@ -322,10 +543,14 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
     /* Calculate the cosine of the angle between the track direction and the
      * vector to the PMT. */
     cos_theta = DOT(dir,pmt_dir);
-
-    if (fabs(cos_theta-1.0/(n*beta))/theta0 > 5) return;
+    sin_theta = sqrt(1-cos_theta*cos_theta);
 
     cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);
+    sin_theta_pmt = sqrt(1-cos_theta_pmt*cos_theta_pmt);
+
+    theta0 = fmax(theta0,get_theta0_min(R,PMT_RADIUS,sin_theta_pmt));
+
+    if (fabs(cos_theta-1.0/(n*beta))/(sin_theta*theta0) > 5) return;
 
     if (cos_theta_pmt > 0) return;
 
@@ -398,7 +623,7 @@ double time_pdf(double t, double mu_noise, double mu_indirect, double *mu, size_
      * function around the times `ts`, where each element of `ts` comes from a
      * different particle.  This assumption is probably valid for particles
      * like muons which don't scatter much, and the hope is that it is *good
-     * enough* for electrons too.  The probability distribution for indirect
+     * enough* for electrons too. The probability distribution for indirect
      * light is assumed to be a step function past some time `tmean`.
      *
      * The probability returned is calculated by taking the sum of these three
@@ -437,66 +662,10 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
         return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma));
 }
 
-static void integrate_path_delta_ray(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
-{
-    size_t i;
-    double pos[3], t, q, r, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
-
-    q_sum = 0.0;
-    r_sum = 0.0;
-    t_sum = 0.0;
-    t2_sum = 0.0;
-    for (i = 0; i < n; i++) {
-        pos[0] = pos0[0] + x[i]*dir0[0];
-        pos[1] = pos0[1] + x[i]*dir0[1];
-        pos[2] = pos0[2] + x[i]*dir0[2];
-
-        get_expected_charge_delta_ray(p, pos, dir0, pmt, &q, &r, &t);
-
-        q = q*pdf[i];
-        r = r*pdf[i];
-        t += x[i]/SPEED_OF_LIGHT;
-
-        if (i == 0 || i == (n - 1)) {
-            q_sum += q;
-            r_sum += r;
-            t_sum += t*q;
-            t2_sum += t*t*q;
-        } else {
-            q_sum += 2*q;
-            r_sum += 2*r;
-            t_sum += 2*t*q;
-            t2_sum += 2*t*t*q;
-        }
-    }
-
-    dx = x[1] - x[0];
-    *mu_direct = q_sum*dx*0.5;
-    *mu_indirect = r_sum*dx*0.5;
-
-    if (*mu_direct == 0.0) {
-        *time = 0.0;
-        *sigma = PMT_TTS;
-    } else {
-        *time = t_sum*dx*0.5/(*mu_direct);
-
-        /* Variance in the time = E(t^2) - E(t)^2. */
-        tmp = t2_sum*dx*0.5/(*mu_direct) - (*time)*(*time);
-
-        if (tmp >= 0) {
-            *sigma = sqrt(tmp + PMT_TTS*PMT_TTS);
-        } else {
-            /* This should never happen but does occasionally due to floating
-             * point rounding error. */
-            *sigma = PMT_TTS;
-        }
-    }
-}
-
-static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
+static void integrate_path_shower(particle *p, double *x, double *w, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
 {
     size_t i;
-    double pos[3], t, q, r, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
+    double pos[3], t, q, r, tmp, q_sum, r_sum, t_sum, t2_sum;
 
     q_sum = 0.0;
     r_sum = 0.0;
@@ -509,35 +678,25 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
 
         get_expected_charge_shower(p, pos, dir0, pmt, &q, &r, &t);
 
-        q = q*pdf[i];
-        r = r*pdf[i];
         t += x[i]/SPEED_OF_LIGHT;
 
-        if (i == 0 || i == (n - 1)) {
-            q_sum += q;
-            r_sum += r;
-            t_sum += t*q;
-            t2_sum += t*t*q;
-        } else {
-            q_sum += 2*q;
-            r_sum += 2*r;
-            t_sum += 2*t*q;
-            t2_sum += 2*t*t*q;
-        }
+        q_sum += q*w[i];
+        r_sum += r*w[i];
+        t_sum += t*q*w[i];
+        t2_sum += t*t*q*w[i];
     }
 
-    dx = x[1]-x[0];
-    *mu_direct = q_sum*dx*0.5;
-    *mu_indirect = r_sum*dx*0.5;
+    *mu_direct = q_sum;
+    *mu_indirect = r_sum;
 
     if (*mu_direct == 0.0) {
         *time = 0.0;
         *sigma = PMT_TTS;
     } else {
-        *time = t_sum*dx*0.5/(*mu_direct);
+        *time = t_sum/(*mu_direct);
 
         /* Variance in the time = E(t^2) - E(t)^2. */
-        tmp = t2_sum*dx*0.5/(*mu_direct) - (*time)*(*time);
+        tmp = t2_sum/(*mu_direct) - (*time)*(*time);
 
         if (tmp >= 0) {
             *sigma = sqrt(tmp + PMT_TTS*PMT_TTS);
@@ -757,7 +916,7 @@ 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, mom, beta;
+    double T, E, mom, beta, tmp;
     betaRootParams *pars;
 
     pars = (betaRootParams *) params;
@@ -766,8 +925,15 @@ static double beta_root(double x, void *params)
 
     /* Calculate total energy */
     E = T + pars->p->mass;
-    mom = sqrt(E*E - pars->p->mass*pars->p->mass);
-    beta = mom/E;
+
+    tmp = E*E - pars->p->mass*pars->p->mass;
+
+    if (tmp <= 0) {
+        beta = 0.0;
+    } else {
+        mom = sqrt(tmp);
+        beta = mom/E;
+    }
 
     return beta - pars->beta_min;
 }
@@ -801,7 +967,7 @@ static int get_smax(particle *p, double beta_min, double range, double *smax)
         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);
+        status = gsl_root_test_interval(x_lo, x_hi, 1e-10, 0);
 
         if (status == GSL_SUCCESS) break;
     } while (status == GSL_CONTINUE && iter < max_iter);
@@ -964,19 +1130,23 @@ double nll_best(event *ev)
     return kahan_sum(nll,nhit);
 }
 
-double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const int fast, int charge_only, int hit_only)
+/* Returns the negative log likelihood for event `ev` given a particle with id
+ * `id`, initial kinetic energy `T0`, position `pos`, direction `dir` and
+ * initial time `t0`.
+ *
+ * `dx` is the distance the track is sampled at to compute the likelihood.
+ *
+ * `ns` is the number of points to use when calculating the number of photons
+ * from shower and delta ray particles.
+ *
+ * If `fast` is nonzero, returns an approximate faster version of the
+ * likelihood.
+ *
+ * `z1` and `z2` should be arrays of length `n` and represent coefficients in
+ * the Karhunen Loeve expansion of the track path. These are only used when
+ * `fast` is zero. */
+double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, int charge_only, int hit_only)
 {
-    /* Returns the negative log likelihood for event `ev` given a particle with
-     * id `id`, initial kinetic energy `T0`, position `pos`, direction `dir` and
-     * initial time `t0`.
-     *
-     * `dx` is the distance the track is sampled at to compute the likelihood.
-     *
-     * If `fast` is nonzero, returns an approximate faster version of the likelihood.
-     *
-     * `z1` and `z2` should be arrays of length `n` and represent coefficients
-     * in the Karhunen Loeve expansion of the track path. These are only used
-     * when `fast` is zero. */
     size_t i, j, k, nhit;
     static double logp[MAX_PE], nll[MAX_PMTS];
     double range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
@@ -1007,26 +1177,38 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in
     static double mu_sum[MAX_PMTS];
     double mu_noise, mu_indirect[MAX_VERTICES], mu_indirect_total;
 
-    static double x[MAX_NPOINTS], pdf[MAX_NPOINTS];
-    static double x_delta[MAX_NPOINTS], pdf_delta[MAX_NPOINTS];
+    /* Points and weights to integrate the shower and delta ray photon
+     * distributions. */
+    static double x[MAX_NPOINTS], w[MAX_NPOINTS];
 
     double result;
 
     double min_ratio;
 
-    double pos_a, pos_b;
+    size_t npoints;
 
-    double xlo, xhi;
+    double log_pt_fast;
 
-    size_t npoints, npoints_shower;
+    double distance_to_psup;
 
-    double log_pt_fast;
+    double q_indirect;
+
+    double R, cos_gamma, sin_gamma;
+
+    double pmt_dir[3];
+
+    double cerenkov_range;
 
     if (n > MAX_VERTICES) {
         fprintf(stderr, "maximum number of vertices is %i\n", MAX_VERTICES);
         exit(1);
     }
 
+    if (ns > MAX_NPOINTS) {
+        fprintf(stderr, "maximum number of points is %i\n", MAX_NPOINTS);
+        exit(1);
+    }
+
     /* Initialize static arrays. */
     for (i = 0; i < MAX_PMTS; i++) {
         mu_sum[i] = 0.0;
@@ -1040,6 +1222,18 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in
     }
 
     for (j = 0; j < n; j++) {
+        /* Find the distance from the particle's starting position to the PSUP.
+         *
+         * We calculate this here to make sure that we don't integrate outside
+         * of the PSUP because otherwise there is no code to model the fact
+         * that the PSUP is optically separated from the water outside. */
+        if (!intersect_sphere(v[j].pos,v[j].dir,PSUP_RADIUS,&distance_to_psup)) {
+            /* If the particle is outside the PSUP we just assume it produces
+             * no light. */
+            fprintf(stderr, "particle doesn't intersect PSUP!\n");
+            continue;
+        }
+
         if (fast)
             p = particle_init(v[j].id, v[j].T0, 1000);
         else
@@ -1047,12 +1241,6 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in
 
         range = p->range;
 
-        /* Number of points to sample along the particle track. */
-        npoints = (size_t) (range/dx + 0.5);
-
-        if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS;
-        if (npoints > MAX_NPOINTS) npoints = MAX_NPOINTS;
-
         /* Calculate total energy */
         E0 = v[j].T0 + p->mass;
         p0 = sqrt(E0*E0 - p->mass*p->mass);
@@ -1061,63 +1249,23 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in
         /* FIXME: is this formula valid for muons? */
         theta0 = get_scattering_rms(range/2,p0,beta0,1.0)/sqrt(range/2);
 
-        if (!fast) {
-            path = path_init(v[j].pos, v[j].dir, v[j].T0, range, npoints, theta0, getKineticEnergy, p, v[j].z1, v[j].z2, v[j].n, p->mass);
-
-            switch (v[j].id) {
-            case IDP_E_MINUS:
-            case IDP_E_PLUS:
-                electron_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b);
-                break;
-            case IDP_MU_MINUS:
-            case IDP_MU_PLUS:
-                muon_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b);
-                break;
-            default:
-                fprintf(stderr, "unknown particle id %i\n", v[j].id);
-                exit(1);
-            }
-
-            /* Lower and upper limits to integrate for the electromagnetic
-             * shower. */
-            xlo = 0.0;
-            xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b);
-
-            /* Number of points to sample along the longitudinal direction for
-             * the electromagnetic shower. */
-            npoints_shower = (size_t) ((xhi-xlo)/dx_shower + 0.5);
-
-            if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS;
-            if (npoints_shower > MAX_NPOINTS) npoints_shower = MAX_NPOINTS;
+        /* Compute the distance at which the particle falls below the cerenkov
+         * threshold. */
+        if (beta0 > 1/avg_index_d2o)
+            get_smax(p, 1/avg_index_d2o, range, &cerenkov_range);
+        else
+            cerenkov_range = range;
 
-            for (i = 0; i < npoints_shower; i++) {
-                x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1);
-                pdf[i] = gamma_pdf(x[i],pos_a,pos_b);
-            }
+        if (distance_to_psup < cerenkov_range)
+            cerenkov_range = distance_to_psup;
 
-            /* Make sure the PDF is normalized. */
-            double norm = trapz(pdf,x[1]-x[0],npoints_shower);
-            for (i = 0; i < npoints_shower; i++) {
-                pdf[i] /= norm;
-            }
+        /* Number of points to sample along the particle track. */
+        npoints = (size_t) (cerenkov_range/dx + 0.5);
 
-            for (i = 0; i < npoints_shower; i++) {
-                x_delta[i] = i*p->range/(npoints_shower-1);
-                /* For now, we assume the delta rays are produced uniformly
-                 * along the path. This isn't quite correct, but the
-                 * distribution is close to flat.
-                 *
-                 * FIXME: Figure out a good parameterization for the
-                 * distribution of delta ray photons along the path. */
-                pdf_delta[i] = 1.0;
-            }
+        if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS;
 
-            /* Make sure the PDF is normalized. */
-            norm = trapz(pdf_delta,x_delta[1]-x_delta[0],npoints_shower);
-            for (i = 0; i < npoints_shower; i++) {
-                pdf_delta[i] /= norm;
-            }
-        }
+        if (!fast)
+            path = path_init(v[j].pos, v[j].dir, v[j].T0, cerenkov_range, npoints, theta0, getKineticEnergy, p, v[j].z1, v[j].z2, v[j].n, p->mass);
 
         if (beta0 > BETA_MIN)
             get_smax(p, BETA_MIN, range, &smax);
@@ -1133,6 +1281,9 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in
 
         mu_indirect[j] = 0.0;
         for (i = 0; i < MAX_PMTS; i++) {
+            /* Note: We loop over all normal PMTs here, even ones that are
+             * flagged since they will contribute to the reflected and
+             * scattered light. */
             if (pmts[i].pmt_type != PMT_NORMAL) continue;
 
             /* Skip PMTs which weren't hit when doing the "fast" likelihood
@@ -1146,7 +1297,6 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in
                 continue;
             }
 
-            double q_indirect;
             integrate_path(path, i, &mu[i][j][0], &q_indirect, &result);
             mu_indirect[j] += q_indirect;
 
@@ -1166,14 +1316,26 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in
                 ts[i][j][0] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,cos_theta_cerenkov,sin_theta_cerenkov);
             }
 
+            /* Calculate the distance and angle to the PMT from the starting
+             * position. These are used to convert between the particle's x
+             * position along the straight track and the cosine of the angle
+             * with the PMT. */
+            SUB(pmt_dir,pmts[i].pos,v[j].pos);
+
+            R = NORM(pmt_dir);
+
+            cos_gamma = DOT(v[j].dir,pmt_dir)/R;
+            sin_gamma = sqrt(1-cos_gamma*cos_gamma);
+
             if (p->shower_photons > 0) {
-                integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu[i][j][1],&q_indirect,&result,&ts_sigma[i][j][1]);
+                get_shower_weights(p,distance_to_psup,x,w,ns,R,cos_gamma,sin_gamma);
+                integrate_path_shower(p,x,w,v[j].T0,v[j].pos,v[j].dir,i,ns,&mu[i][j][1],&q_indirect,&result,&ts_sigma[i][j][1]);
                 mu_indirect[j] += q_indirect;
                 ts[i][j][1] = v[j].t0 + result;
             }
 
-            if (p->delta_ray_photons > 0) {
-                integrate_path_delta_ray(p,x_delta,pdf_delta,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu[i][j][2],&q_indirect,&result,&ts_sigma[i][j][2]);
+            if (p->delta_ray_photons > 0 && !get_delta_ray_weights(p,range,distance_to_psup,x,w,ns,R,cos_gamma,sin_gamma)) {
+                integrate_path_shower(p,x,w,v[j].T0,v[j].pos,v[j].dir,i,ns,&mu[i][j][2],&q_indirect,&result,&ts_sigma[i][j][2]);
                 mu_indirect[j] += q_indirect;
                 ts[i][j][2] = v[j].t0 + result;
             }
diff --git a/src/likelihood.h b/src/likelihood.h
index bda46f0..d739802 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -19,6 +19,8 @@
 
 #include "event.h"
 #include <stddef.h> /* for size_t */
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_spline.h>
 
 #define PSUP_REFLECTION_TIME 80.0
 
@@ -80,6 +82,20 @@ typedef struct particle {
     double range;
     double *x;
     double *T;
+    double *cos_theta;
+    double *cdf_shower;
+    /* Spline for looking up cdf -> cos(theta). */
+    gsl_spline *spline_shower;
+    gsl_interp_accel *acc_shower;
+    double *cdf_delta;
+    double *x_shower;
+    double *gamma_pdf;
+    double xlo_shower;
+    double xhi_shower;
+    double pos_a;
+    double pos_b;
+    gsl_spline *spline_delta;
+    gsl_interp_accel *acc_delta;
     size_t n;
     double a;
     double b;
@@ -95,6 +111,6 @@ void particle_free(particle *p);
 double time_pdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma);
 double time_cdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma);
 double nll_best(event *ev);
-double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast, int charge_only, int hit_only);
+double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, int charge_only, int hit_only);
 
 #endif
diff --git a/src/misc.c b/src/misc.c
index e6939f7..f705c74 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -862,3 +862,22 @@ void get_dir(double *dir, double theta, double phi)
     dir[1] = sin_theta*sin_phi;
     dir[2] = cos_theta;
 }
+
+/* Fast version of acos() which uses a lookup table computed on the first call. */
+double fast_acos(double x)
+{
+    size_t i;
+    static int initialized = 0;
+    static double xs[N_ACOS];
+    static double ys[N_ACOS];
+
+    if (!initialized) {
+        for (i = 0; i < LEN(xs); i++) {
+            xs[i] = -1.0 + 2.0*i/(LEN(xs)-1);
+            ys[i] = acos(xs[i]);
+        }
+        initialized = 1;
+    }
+
+    return interp1d(x,xs,ys,LEN(xs));
+}
diff --git a/src/misc.h b/src/misc.h
index e06e323..60e54ea 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -27,6 +27,9 @@
 #define LN_MAX 100
 #define LNFACT_MAX 100
 
+/* Number of points to create a lookup table for the fast_acos() function. */
+#define N_ACOS 10000
+
 double trapz(const double *y, double dx, size_t n);
 int intersect_sphere(double *pos, double *dir, double R, double *l);
 void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2);
@@ -51,5 +54,6 @@ void combinations_with_replacement(size_t n, size_t r, size_t *result, size_t *l
 size_t argmax(double *a, size_t n);
 size_t argmin(double *a, size_t n);
 void get_dir(double *dir, double theta, double phi);
+double fast_acos(double x);
 
 #endif
diff --git a/src/scattering.c b/src/scattering.c
index 92e5df0..5a7f7a9 100644
--- a/src/scattering.c
+++ b/src/scattering.c
@@ -236,4 +236,7 @@ void free_interpolation(void)
     gsl_spline2d_free(spline);
     gsl_interp_accel_free(xacc);
     gsl_interp_accel_free(yacc);
+
+    gsl_spline_free(spline2);
+    gsl_interp_accel_free(xacc2);
 }
diff --git a/src/sno.h b/src/sno.h
index bf0a1eb..99847c5 100644
--- a/src/sno.h
+++ b/src/sno.h
@@ -33,6 +33,6 @@
 #define AV_RADIUS 603.25
 #define AV_THICKNESS 6.5
 /* FIXME: is this right? */
-#define PSUP_RADIUS 800.0
+#define PSUP_RADIUS 840.0
 
 #endif