path: root/src/likelihood.c

/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
 *
 * This program is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option)
 * any later version.

 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details.

 * You should have received a copy of the GNU General Public License along with
 * this program. If not, see <https://www.gnu.org/licenses/>.
 */

#include "likelihood.h"
#include <stdlib.h> /* for size_t */
#include "pmt.h"
#include <gsl/gsl_integration.h>
#include "muon.h"
#include "misc.h"
#include <gsl/gsl_sf_gamma.h>
#include "sno.h"
#include "vector.h"
#include "event.h"
#include "optics.h"
#include "sno_charge.h"
#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>
#include "pmt_response.h"
#include "electron.h"
#include "proton.h"
#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)
{
    /* Returns a struct with arrays of the particle position and kinetic
     * energy. This struct can then be passed to particle_get_energy() to
     * interpolate the particle's kinetic energy at any point along the track.
     * For example:
     *
     *     particle *p = particle_init(IDP_MU_MINUS, 1000.0, 100);
     *     double T = particle_get_energy(x, p);
     *
     * To compute the kinetic energy as a function of distance we need to solve
     * the following integral equation:
     *
     *     T(x) = T0 - \int_0^x dT(T(x))/dx
     *
     * 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, norm, pdf, pdf_last;

    particle *p = malloc(sizeof(particle));
    p->id = id;
    p->x = malloc(sizeof(double)*n);
    p->T = malloc(sizeof(double)*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;
    p->T[0] = T0;

    switch (id) {
    case IDP_E_MINUS:
    case IDP_E_PLUS:
        p->mass = ELECTRON_MASS;
        /* We use the density of light water here since we don't have the
         * tables for heavy water for electrons. */
        p->range = electron_get_range(T0, WATER_DENSITY);

        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
         * constant between points. */
        rad = 0.0;
        for (i = 1; i < n; i++) {
            p->x[i] = p->range*i/(n-1);
            dEdx = electron_get_dEdx(p->T[i-1], WATER_DENSITY);
            p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
            dEdx = electron_get_dEdx_rad(p->T[i-1], WATER_DENSITY);
            rad += dEdx*(p->x[i]-p->x[i-1]);
        }
        /* Make sure that the energy is zero at the last step. This is so that
         * when we try to bisect the point along the track where the speed of
         * the particle is equal to BETA_MIN, we guarantee that there is a
         * point along the track where the speed drops below BETA_MIN.
         *
         * A possible future improvement would be to dynamically compute the
         * range here using the dE/dx table instead of reading in the range. */
        p->T[n-1] = 0;

        p->shower_photons = electron_get_shower_photons(T0, rad);

        break;
    case IDP_MU_MINUS:
    case IDP_MU_PLUS:
        p->mass = MUON_MASS;
        /* We use the density of heavy water here since we only load the heavy
         * water table for muons.
         *
         * FIXME: It would be nice in the future to load both tables and then
         * load the correct table based on the media. */
        p->range = muon_get_range(T0, HEAVY_WATER_DENSITY);

        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);

        /* Now we loop over the points along the track and assume dE/dx is
         * constant between points. */
        rad = 0.0;
        for (i = 1; i < n; i++) {
            p->x[i] = p->range*i/(n-1);
            dEdx = muon_get_dEdx(p->T[i-1], HEAVY_WATER_DENSITY);
            p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
            dEdx = muon_get_dEdx_rad(p->T[i-1], HEAVY_WATER_DENSITY);
            rad += dEdx*(p->x[i]-p->x[i-1]);
        }
        /* Make sure that the energy is zero at the last step. This is so that
         * when we try to bisect the point along the track where the speed of
         * the particle is equal to BETA_MIN, we guarantee that there is a
         * point along the track where the speed drops below BETA_MIN.
         *
         * A possible future improvement would be to dynamically compute the
         * range here using the dE/dx table instead of reading in the range. */
        p->T[n-1] = 0;

        p->shower_photons = muon_get_shower_photons(T0, rad);

        break;
    case IDP_PROTON:
        p->mass = PROTON_MASS;
        /* We use the density of light water here since we don't have the
         * 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
         * constant between points. */
        rad = 0.0;
        for (i = 1; i < n; i++) {
            p->x[i] = p->range*i/(n-1);
            dEdx = proton_get_dEdx(p->T[i-1], WATER_DENSITY);
            p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
            dEdx = proton_get_dEdx_rad(p->T[i-1], WATER_DENSITY);
            rad += dEdx*(p->x[i]-p->x[i-1]);
        }
        /* Make sure that the energy is zero at the last step. This is so that
         * when we try to bisect the point along the track where the speed of
         * the particle is equal to BETA_MIN, we guarantee that there is a
         * point along the track where the speed drops below BETA_MIN.
         *
         * A possible future improvement would be to dynamically compute the
         * range here using the dE/dx table instead of reading in the range. */
        p->T[n-1] = 0;

        /* FIXME: Add shower photons for protons. */
        p->shower_photons = 0.0;

        break;
    default:
        fprintf(stderr, "unknown particle id %i\n", id);
        exit(1);
    }

    if (p->shower_photons) {
        norm = electron_get_angular_pdf_norm(p->a, p->b, 1/get_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.
             *
             * Note: We add EPSILON to the pdf since the cdf values are
             * required to be strictly increasing in order to use gsl splines. */
            pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/get_avg_index_d2o())/norm + EPSILON;
            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);
        }
    }

    if (p->delta_ray_photons) {
        /* Calculate the normalization constant for the angular distribution. */
        norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/get_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.
             *
             * Note: We add EPSILON to the pdf since the cdf values are
             * required to be strictly increasing in order to use gsl splines. */
            pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/get_avg_index_d2o())/norm + EPSILON;
            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);
    }

    return 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)
{
    return (R*cos_gamma-x)/sqrt(R*R+x*x-2*x*R*cos_gamma);
}

/* 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 cos(theta) == 1.0 we must be at the beginning of the track. */
    if (cos_theta == 1.0) return R*cos_gamma;

    /* 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);

    /* This result comes from using the law of sines. */
    return R*(cos_gamma - cos_theta*sin_gamma/sin_theta);
}

/* 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)
{
    size_t i;
    double x1, x2, cdf1, cdf2, cdf, delta, cos_theta, cos_theta_last, xcdf, xcdf_last;

    /* First we need to figure out the bounds of integration in x. */
    x1 = 0.0;

    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++) {
        if (i < n)
            cdf = cdf1 + i*delta;
        else
            cdf = cdf2;

        cos_theta = gsl_spline_eval(p->spline_delta,cdf,p->acc_delta);
        xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma);

        /* Occasionally due to floating point rounding error, it's possible
         * that cdf > cdf2 and so we end up with xcdf = inf. To prevent this,
         * we just make sure that xcdf is bounded by x1 and x2. */
        if (xcdf < x1) xcdf = x1;
        if (xcdf > x2) xcdf = x2;

        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;
    }

    return 0;
}

/* 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;

    /* First we need to figure out the bounds of integration in x. */

    /* 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;

    x1 = p->xlo_shower;

    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++) {
        if (i < n)
            cdf = cdf1 + i*delta;
        else
            cdf = cdf2;

        cos_theta = gsl_spline_eval(p->spline_shower,cdf,p->acc_shower);
        xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma);

        /* Occasionally due to floating point rounding error, it's possible
         * that cdf > cdf2 and so we end up with xcdf = inf. To prevent this,
         * we just make sure that xcdf is bounded by x1 and x2. */
        if (xcdf < x1) xcdf = x1;
        if (xcdf > x2) xcdf = x2;

        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;
    }

    return 0;
}

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;

    T = interp1d(x,p->x,p->T,p->n);

    if (T < 0) return 0;

    return T;
}

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], omega, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o;

    SUB(pmt_dir,pmts[pmt].pos,pos);

    normalize(pmt_dir);

    cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);

    omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);

    f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
    f = get_weighted_pmt_response(-cos_theta_pmt);

    get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);

    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);

    charge = get_weighted_quantum_efficiency()*omega/(2*M_PI);

    /* Note: We multiply the amount of Rayleigh scattered light here by 2 since
     * we are only calculating the fraction of Rayleigh scattered light that
     * would hit the PMT and the coverage in SNO+ is about 50%. */
    *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + 2*prob_sct*charge;

    *t = (l_d2o*get_avg_index_d2o() + l_h2o*get_avg_index_h2o())/SPEED_OF_LIGHT;

    *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
}

/* Returns the angular width of the PMT bucket from the center of the PMT to
 * the edge of the bucket in the plane formed by the particle direction and the
 * direction to the PMT.
 *
 * `R` is the distance to the PMT, `r` is the PMT radius, and `sin_theta_pmt`
 * is the sine of the angle between the vector from the particle position to
 * the PMT and the PMT normal vector.
 *
 * This function is called get_theta0_min because we use this angular width to
 * set a minimum value for the RMS scattering angle of the particle as a kind
 * of hack to deal with the fact that we assume that the angular distribution
 * is constant across the face of the PMT. By introducing a minimum value for
 * the scattering RMS we broaden the angular distribution such that it
 * effectively averages across the face of a PMT. */
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, sin_theta, n, omega, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o, cos_theta_cerenkov, distance_to_pmt;

    /* Previously the index of refraction was calculated based on the position,
     * i.e.:
     *
     *     n = (NORM(pos) <= AV_RADIUS) ? get_avg_index_d2o() : get_avg_index_h2o();
     *
     * but this was causing a discontinuity in the likelihood function when the
     * particle track crossed the AV since we're numerically integrating over a
     * finite set of points. */
    n = get_avg_index_d2o();

    cos_theta_cerenkov = 1/(beta*n);

    *q = 0.0;
    *t = 0.0;
    *reflected = 0.0;
    if (cos_theta_cerenkov > 1) return;

    SUB(pmt_dir,pmts[pmt].pos,pos);

    distance_to_pmt = 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);
    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(distance_to_pmt,PMT_RADIUS,sin_theta_pmt));

#ifdef FAST_GET_EXPECTED_CHARGE
    /* This next line is used to skip out of calculating the expected charge if
     * abs((cos(theta)-cos_theta_cherenkov)/(sin(theta)*theta0)) > 5. The idea
     * here is that later in the likelihood calculation the PDF for Cerenkov
     * light has a term like
     *
     *     exp(((cos(theta)-cos_theta_cherenkov)/(sin(theta)*theta0))**2)
     *
     * which will be less than 10^-25 if the following inequality holds and
     * therefore the charge should be negligible.
     *
     * However! I noticed that this line was causing discontinuities in the
     * likelihood when fitting low energy muons. I realized there are two
     * potential issues with this. One is that the PDF is multiplied by two
     * other terms: 1/(sin_theta*theta0) and the solid angle of the PMT.
     * Therefore if these are large it may have an impact. Second, the PDF
     * actually integrates over the Cerenkov spectrum and correctly takes into
     * account the change in the index as a function of wavelength.
     *
     * For now we don't use this to prevent the discontinuities. However, it
     * would be nice to use it in the future since it provides a massive ~2.5x
     * speedup. One idea for trying to eventually include it again is to update
     * the PDF for Cerenkov light to not include the wavelength dependence
     * (which would make it slightly less accurate, but it may be worth the
     * 2.5x speedup). */

    if (fabs(cos_theta-cos_theta_cerenkov) > 5*sin_theta*theta0) return;
#endif

    omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);

    f_reflec = get_weighted_pmt_reflectivity(cos_theta_pmt);
    f = get_weighted_pmt_response(cos_theta_pmt);

    get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);

    *t = (l_d2o*get_avg_index_d2o() + l_h2o*get_avg_index_h2o())/SPEED_OF_LIGHT;

    /* Probability that a photon is absorbed. We calculate this by computing:
     *
     *     1.0 - P(not absorbed in D2O)*P(not absorbed in H2O)*P(not absorbed in acrylic)
     *
     * since if we worked with the absorption probabilities directly it would
     * be more complicated, i.e.
     *
     *     P(absorbed in D2O) + P(absorbed in acrylic|not absorbed in D2O)*P(not absorbed in D2O) + ...
     *
     */
    prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS);
    /* Similiar calculation for the probability that a photon is scattered.
     *
     * Technically we should compute this conditionally on the probability that
     * a photon is not absorbed, but since the probability of scattering is
     * pretty low, this is expected to be a very small effect. */
    prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS);

    charge = omega*(1-cos_theta_cerenkov*cos_theta_cerenkov)*get_probability(beta, cos_theta, sin_theta, theta0);

    /* Note: We multiply the amount of Rayleigh scattered light here by 2 since
     * we are only calculating the fraction of Rayleigh scattered light that
     * would hit the PMT and the coverage in SNO+ is about 50%. */
    *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + 2*prob_sct*charge;

    *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
}

double time_cdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma)
{
    /* Returns the CDF for the time distribution of photons at time `t`. */
    size_t i;
    double p, mu_total;

    p = mu_noise*t/GTVALID;

    if (t < tmean)
        ;
    else if (t < tmean + 2*PSUP_REFLECTION_TIME)
        p +=  mu_indirect*(t-tmean)/(2*PSUP_REFLECTION_TIME);
    else
        p +=  mu_indirect;

    mu_total = mu_noise + mu_indirect;
    for (i = 0; i < n; i++) {
        if (mu[i] == 0.0) continue;
        p += mu[i]*norm_cdf(t,ts[i],ts_sigma[i]);
        mu_total += mu[i];
    }

    return p/mu_total;
}

/* Returns the probability that a photon is detected at time `t`.
 *
 * The probability distribution is the sum of three different components: dark
 * noise, indirect light, and direct light. The dark noise is assumed to be
 * constant in time. The direct light is assumed to be a delta 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 light is assumed to be a step
 * function lasting 2*PSUP_REFLECTION_TIME past some time `tmean`.
 *
 * Note: Initially I had the indirect light be a step function which only
 * lasted PSUP_REFLECTION_TIME (which is set to the time it takes light to
 * cross from one side of the PSUP to the other) because I thought that this
 * represented the maximum amount of time it took for light to traverse the
 * detector (and the PMT hit time plot for a laserball run has a bump after the
 * main peak which lasts for about 80 ns). However, I later realized that for
 * events near the edge of the detector, like say a flasher or muon, the
 * reflected light from the other side of the detector will take twice as long.
 * So, now we use twice the PSUP reflection time.
 *
 * In the future it might be cool to experiment with new ways of defining the
 * time PDF. For example, we could get a much more accurate approximation for
 * the reflected light time distribution by using something like a kernel
 * density estimator and for each PMT we could loop over the other PMTs and add
 * a point to our PDF at the light travel time between the PMTs weighted by the
 * second PMT's charge. I think currently this would be way too slow, and it's
 * not obvious how much it would help.
 *
 * The probability returned is calculated by taking the sum of these three
 * components and convolving it with a gaussian with standard deviation `sigma`
 * which should typically be the PMT transit time spread. */
double time_pdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma)
{
    size_t i;
    double p, mu_total;

    p = mu_noise/GTVALID + mu_indirect*(erf((t-tmean)/(sqrt(2)*PMT_TTS))-erf((t-tmean-2*PSUP_REFLECTION_TIME)/(sqrt(2)*PMT_TTS)))/(4*PSUP_REFLECTION_TIME);

    mu_total = mu_noise + mu_indirect;
    for (i = 0; i < n; i++) {
        if (mu[i] == 0.0) continue;
        p += mu[i]*norm(t,ts[i],ts_sigma[i]);
        mu_total += mu[i];
    }

    return p/mu_total;
}

/* Returns the first order statistic for observing a PMT hit at time `t` given
 * `n` hits.
 *
 * The first order statistic is computed from the probability distribution
 * above. It's not obvious whether one should take the first order statistic
 * before or after convolving with the PMT transit time spread. Since at least
 * some of the transit time spread in SNO comes from the different transit
 * times across the face of the PMT, it seems better to convolve first which is
 * what we do here. In addition, the problem is not analytically tractable if
 * you do things the other way around.
 *
 * Note: Technically we should include a final convolution representing both
 * the ECA+PCA uncertainty and the digitization noise. I was initially worried
 * about this causing an issue because as you take higher and higher order
 * statistics the width of the time PDF becomes smaller and smaller and these
 * other uncertainties may become relevant. However, the maximum number of PE
 * we should be calculating the time PDF for is when QLX rails which is at
 * approximately:
 *
 *     (4096-600)*12/30 ~= 1000
 *
 * and the width at that point is ~0.6 ns. The digitization noise is
 * approximately 4096/400.0 ~= 0.1 ns and according to Javi:
 *
 * > The ECA+PCA uncertainties are negligible with respect to the TTS.
 * > Basically, the width of the prompt peak for N16 calibration events (mainly
 * > SPEs) is 1.6ns, compatible with the PMT TTS.
 * >
 * > That said, long ago I estimated the precision of the ECA calibration to be
 * > about 0.4ns. The PCA delays are pretty stable and precisely measured, so
 * > probably negligible compared to the ECA uncertainties.
 *
 * So therefore, even at the highest possible charge we should still be
 * dominated by the width of the PDF.
 *
 * It is probably too expensive to try and actually calculate another
 * convolution at this point since we'd have to do it numerically. I think a
 * better option would just be to cut off the first order statistic at some
 * point, i.e. to calculate:
 *
 *    logp[j] += log_pt(ev->pmt_hits[i].t, j > MAX_PE_TIME_PDF ? MAX_PE_TIME_PDF : j, mu_noise, mu_indirect_total, &mu[i][0][0], n*3, &ts[i][0][0], ts[i][0][0], &ts_sigma[i][0][0]);
 *
 */
double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu, size_t n2, double *ts, double tmean, double *ts_sigma)
{
    if (n == 1)
        return ln(n) + log(time_pdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma));
    else
        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_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, 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_shower(p, pos, dir0, pmt, &q, &r, &t);

        t += x[i]/SPEED_OF_LIGHT;

        q_sum += q*w[i];
        r_sum += r*w[i];
        t_sum += t*q*w[i];
        t2_sum += t*t*q*w[i];
    }

    *mu_direct = q_sum;
    *mu_indirect = r_sum;

    if (*mu_direct == 0.0) {
        *time = 0.0;
        *sigma = PMT_TTS;
    } else {
        *time = t_sum/(*mu_direct);

        /* Variance in the time = E(t^2) - E(t)^2. */
        tmp = t2_sum/(*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(path *p, int pmt, double *mu_direct, double *mu_indirect, double *time)
{
    size_t i;
    double dir[3], pos[3], t, theta0, beta, q, r, q_sum, r_sum, t_sum, dx;

    q_sum = 0.0;
    r_sum = 0.0;
    t_sum = 0.0;
    for (i = 0; i < p->len; i++) {
        pos[0] = p->x[i];
        pos[1] = p->y[i];
        pos[2] = p->z[i];

        dir[0] = p->dx[i];
        dir[1] = p->dy[i];
        dir[2] = p->dz[i];

        beta = p->beta[i];

        if (p->n > 0) {
            theta0 = p->theta0;
        } else {
            theta0 = fmax(MIN_THETA0,p->theta0*sqrt(p->s[i]));
            theta0 = fmin(MAX_THETA0,theta0);
        }

        get_expected_charge(beta, theta0, pos, dir, pmt, &q, &r, &t);

        t += p->t[i];

        if (i == 0 || i == p->len - 1) {
            q_sum += q;
            r_sum += r;
            t_sum += t*q;
        } else {
            q_sum += 2*q;
            r_sum += 2*r;
            t_sum += 2*t*q;
        }
    }

    dx = p->s[1] - p->s[0];
    *mu_direct = q_sum*dx*0.5;
    *mu_indirect = r_sum*dx*0.5;
    *time = t_sum*dx*0.5;
}

static double get_total_charge_approx(double beta0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected, double cos_theta_cerenkov, double sin_theta_cerenkov)
{
    /* 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, x, z, s, a, b, beta, E, mom, T, omega, sin_theta, f, cos_theta_pmt, l_h2o, l_d2o, f_reflec, charge, prob, frac, prob_abs, prob_sct;

    /* 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. */
    if (sin_theta_cerenkov != 0)
        s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov;
    else
        s = R*cos_theta;

    /* 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);

    /* `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*sqrt(1-pow(cos_theta,2));

    /* `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 = particle_get_energy(s,p);

    /* Calculate the particle velocity at the point `s`. */
    E = T + p->mass;
    mom = sqrt(E*E - p->mass*p->mass);
    beta = mom/E;

    *t = 0.0;
    *mu_reflected = 0.0;
    if (beta < 1/get_avg_index_d2o()) return 0.0;

    /* `prob` is the number of photons emitted per cm by the particle at a
     * distance `s` along the track. */
    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);

    normalize(pmt_dir);

    cos_theta_pmt = DOT(pmt_dir,pmts[i].normal);

    *t = 0.0;
    *mu_reflected = 0.0;
    if (cos_theta_pmt > 0) return 0.0;

    /* Calculate the sine of the angle between the track direction and the PMT
     * at the position `s` along the track. */
    cos_theta = DOT(dir,pmt_dir);
    sin_theta = sqrt(1-pow(cos_theta,2));

    /* Get the solid angle of the PMT at the position `s` along the track. */
    omega = get_solid_angle_fast(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS);

    theta0 = fmax(theta0*sqrt(s),MIN_THETA0);

    frac = sqrt(2)*get_avg_index_d2o()*x*beta0*theta0;

    f = get_weighted_pmt_response(-cos_theta_pmt);
    f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);

    get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);

    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);

    /* Assume the particle is travelling at the speed of light. */
    *t = s/SPEED_OF_LIGHT + l_d2o*get_avg_index_d2o()/SPEED_OF_LIGHT + l_h2o*get_avg_index_h2o()/SPEED_OF_LIGHT;

    charge = get_avg_index_d2o()*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+get_avg_index_d2o()*(smax-z)*beta0)/frac) + erf((-a+b*s+get_avg_index_d2o()*z*beta0)/frac))/(b+get_avg_index_d2o()*beta0)/(4*M_PI);

    /* Add expected number of photons from electromagnetic shower. */
    if (p->shower_photons > 0)
        charge += get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/get_avg_index_d2o())*omega/(2*M_PI);
    if (p->delta_ray_photons > 0)
        charge += get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/get_avg_index_d2o())*omega/(2*M_PI);

    /* Note: We multiply the amount of Rayleigh scattered light here by 2 since
     * we are only calculating the fraction of Rayleigh scattered light that
     * would hit the PMT and the coverage in SNO+ is about 50%. */
    *mu_reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + 2*prob_sct*charge;

    return (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
}

typedef struct betaRootParams {
    particle *p;
    double beta_min;
} betaRootParams;

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, tmp;
    betaRootParams *pars;

    pars = (betaRootParams *) params;

    T = particle_get_energy(x, pars->p);

    /* Calculate total energy */
    E = T + pars->p->mass;

    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;
}

static int get_smax(particle *p, 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;
    betaRootParams pars;
    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);

    pars.p = p;
    pars.beta_min = beta_min;

    F.function = &beta_root;
    F.params = &pars;

    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-10, 0);

        if (status == GSL_SUCCESS) break;
    } while (status == GSL_CONTINUE && iter < max_iter);

    gsl_root_fsolver_free(s);

    *smax = r;

    return status;
}

static double guess_time(double *pos, double *dir, double *pmt_pos, double smax, 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. */
    if (sin_theta_cerenkov != 0)
        s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov;
    else
        s = R*cos_theta;

    /* 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*get_avg_index_d2o()/SPEED_OF_LIGHT + l_h2o*get_avg_index_h2o()/SPEED_OF_LIGHT;
}

static double getKineticEnergy(double x, void *p)
{
    return particle_get_energy(x, (particle *) p);
}

double nll_best(event *ev)
{
    /* Returns the negative log likelihood of the "best hypothesis" for the event `ev`.
     *
     * By "best hypothesis" I mean we restrict the model to assign a single
     * mean number of PE for each PMT in the event assuming that the number of
     * PE hitting each PMT is poisson distributed. In addition, the model picks
     * a single time for the light to arrive with the PDF being given by a sum
     * of dark noise and a Gaussian around this time with a width equal to the
     * single PE transit time spread.
     *
     * This calculation is intended to be used as a goodness of fit test for
     * the fitter by computing a likelihood ratio between the best fit
     * hypothesis and this ideal case. See Chapter 9 in Jayne's "Probability
     * Theory: The Logic of Science" for more details. */
    size_t i, j, nhit;
    static double logp[MAX_PE], nll[MAX_PMTS], mu[MAX_PMTS], ts[MAX_PMTS], ts_sigma;
    double log_mu, max_logp, mu_noise, mu_indirect_total, min_ratio;

    mu_noise = DARK_RATE*GTVALID*1e-9;

    /* Compute the "best" number of expected PE for each PMT. */
    for (i = 0; i < MAX_PMTS; i++) {
        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;

        if (!ev->pmt_hits[i].hit) {
            /* If the PMT wasn't hit we assume the expected number of PE just
             * comes from noise hits. */
            mu[i] = mu_noise;
            continue;
        }

        /* Set the expected number of PE to whatever maximizes P(q|mu), i.e.
         *
         *     mu[i] = max(p(q|mu))
         *
         */
        mu[i] = get_most_likely_mean_pe(ev->pmt_hits[i].q);
        /* We want to estimate the mean time which is most likely to produce a
         * first order statistic of the actual hit time given there were
         * approximately mu[i] PE. As far as I know there are no closed form
         * solutions for the first order statistic of a Gaussian, but there is
         * an approximate form in the paper "Expected Normal Order Statistics
         * (Exact and Approximate)" by J. P. Royston which is what I use here.
         *
         * I found this via the stack exchange question here:
         * https://stats.stackexchange.com/questions/9001/approximate-order-statistics-for-normal-random-variables. */
        double alpha = 0.375;
        ts[i] = ev->pmt_hits[i].t - gsl_cdf_ugaussian_Pinv((1-alpha)/(mu[i]-2*alpha+1))*PMT_TTS;
    }

    ts_sigma = PMT_TTS;

    mu_indirect_total = 0.0;

    min_ratio = MIN_RATIO;

    /* Now, we actually compute the negative log likelihood for the best
     * hypothesis.
     *
     * Currently this calculation is identical to the one in nll() so it should
     * probably be made into a function. */
    nhit = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;

        log_mu = log(mu[i]);

        if (ev->pmt_hits[i].hit) {
            for (j = 1; j < MAX_PE; j++) {
                logp[j] = get_log_pq(ev->pmt_hits[i].q,j) + get_log_phit(j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], 1, &ts[i], ts[i], &ts_sigma);

                if (j == 1 || logp[j] > max_logp) max_logp = logp[j];

                if (logp[j] - max_logp < min_ratio*ln(10)) {
                    j++;
                    break;
                }
            }

            nll[nhit++] = -logsumexp(logp+1, j-1);
        } else {
            logp[0] = -mu[i];
            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, MAX_PE_NO_HIT);
        }
    }

    return kahan_sum(nll,nhit);
}

/* 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)
{
    size_t i, j, k, nhit;
    static double logp[MAX_PE], nll[MAX_PMTS], tmp[2];
    double range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
    particle *p;
    double cos_theta_cerenkov, sin_theta_cerenkov;
    double logp_path;
    double mu_reflected;
    path *path;

    /* Array representing the expected number of photons at each PMT for each
     * vertex. The last axis represents the contributions from direct, shower,
     * and delta ray photons respectively, i.e.
     *
     *     mu[0][2][0] - expected number of direct photons from vertex 2 at PMT 0
     *     mu[0][2][1] - expected number of shower photons from vertex 2 at PMT 0
     *     mu[0][2][2] - expected number of delta ray photons from vertex 2 at PMT 0
     *     ...
     *     etc.
     *
     */
    static double mu[MAX_PMTS][MAX_VERTICES][3];
    /* Array representing the average time of arrival at each PMT for each
     * vertex. See above for the how the array is indexed. */
    static double ts[MAX_PMTS][MAX_VERTICES][3];
    /* Array representing the standard deviation of the time at each PMT for
     * each vertex. See above for the how the array is indexed. */
    static double ts_sigma[MAX_PMTS][MAX_VERTICES][3];
    static double mu_sum[MAX_PMTS];
    double mu_noise, mu_indirect[MAX_VERTICES], mu_indirect_total;

    /* 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;

    size_t npoints;

    double log_pt_fast;

    double distance_to_psup;

    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);
    }

    for (i = 0; i < LEN(mu_indirect); i++) {
        mu_indirect[i] = 0.0;
    }

    /* Initialize static arrays. */
    for (i = 0; i < MAX_PMTS; i++) {
        mu_sum[i] = 0.0;
        for (j = 0; j < n; j++) {
            for (k = 0; k < 3; k++) {
                mu[i][j][k] = 0.0;
                ts[i][j][k] = 0.0;
                ts_sigma[i][j][k] = PMT_TTS;
            }
        }
    }

    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 (NORM(v[j].pos) > PSUP_RADIUS || !intersect_sphere(v[j].pos,v[j].dir,PSUP_RADIUS,&distance_to_psup)) {
            /* If the particle is outside the PSUP or it doesn't intersect 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
            p = particle_init(v[j].id, v[j].T0, 10000);

        range = p->range;

        /* Calculate total energy */
        E0 = v[j].T0 + p->mass;
        p0 = sqrt(E0*E0 - p->mass*p->mass);
        beta0 = p0/E0;

        /* FIXME: is this formula valid for muons? */
        theta0 = get_scattering_rms(range/2,p0,beta0,1.0)/sqrt(range/2);

        /* Compute the distance at which the particle falls below the cerenkov
         * threshold. */
        if (beta0 > 1/get_avg_index_d2o())
            get_smax(p, 1/get_avg_index_d2o(), range, &cerenkov_range);
        else
            cerenkov_range = range;

        if (distance_to_psup < cerenkov_range)
            cerenkov_range = distance_to_psup;

        /* Number of points to sample along the particle track. */
        npoints = (size_t) (cerenkov_range/dx + 0.5);

        if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS;

        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);
        else
            smax = 0.0;

        /* Compute the Cerenkov angle at the start of the track.
         *
         * These are computed at the beginning of this function and then passed to
         * the different functions to avoid recomputing them on the fly. */
        cos_theta_cerenkov = 1/(get_avg_index_d2o()*beta0);
        sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2));

        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
             * calculation. */
            if (hit_only && !ev->pmt_hits[i].hit) continue;

            if (fast) {
                mu[i][j][0] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j][0], &mu_reflected, cos_theta_cerenkov, sin_theta_cerenkov);
                ts[i][j][0] += v[j].t0;
                mu_indirect[j] += mu_reflected;
                continue;
            }

            integrate_path(path, i, &mu[i][j][0], &q_indirect, &result);
            mu_indirect[j] += q_indirect;

            if (mu[i][j][0] > 1e-9) {
                ts[i][j][0] = v[j].t0 + result/mu[i][j][0];
            } else {
                /* If the expected number of PE is very small, our estimate of
                 * the time can be crazy since the error in the total charge is
                 * in the denominator. Therefore, we estimate the time using
                 * the same technique as in get_total_charge_approx(). See that
                 * function for more details.
                 *
                 * Note: I don't really know how much this affects the likelihood.
                 * I should test it to see if we can get away with something even
                 * simpler (like just computing the distance from the start of the
                 * track to the PMT). */
                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) {
                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 && !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;
            }
        }

        if (!fast) {
            path_free(path);
        }

        particle_free(p);
    }

    mu_noise = DARK_RATE*GTVALID*1e-9;

    mu_indirect_total = 0.0;
    for (j = 0; j < n; j++) {
        mu_indirect_total += mu_indirect[j]*CHARGE_FRACTION/10000.0;
    }

    /* Compute the expected number of photons reaching each PMT by adding up
     * the contributions from the noise hits and the direct, indirect, and
     * shower light. */
    for (i = 0; i < MAX_PMTS; i++) {
        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;

        /* Skip PMTs which weren't hit when doing the "fast" likelihood
         * calculation. */
        if (hit_only && !ev->pmt_hits[i].hit) continue;

        mu_sum[i] = mu_indirect_total + mu_noise;
        for (j = 0; j < n; j++) {
            for (k = 0; k < 3; k++) {
                mu_sum[i] += mu[i][j][k];
            }
        }
    }

    if (fast)
        min_ratio = MIN_RATIO_FAST;
    else
        min_ratio = MIN_RATIO;

    /* Now, we actually compute the negative log likelihood. */
    nhit = 0;
    for (i = 0; i < MAX_PMTS; i++) {
        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;

        /* Skip PMTs which weren't hit when doing the "fast" likelihood
         * calculation. */
        if (hit_only && !ev->pmt_hits[i].hit) continue;

        log_mu = log(mu_sum[i]);

        if (fast)
            log_pt_fast = log_pt(ev->pmt_hits[i].t, 1, mu_noise, mu_indirect_total, &mu[i][0][0], n*3, &ts[i][0][0], ts[i][0][0], &ts_sigma[i][0][0]);

        if (ev->pmt_hits[i].hit) {
            for (j = 1; j < MAX_PE; j++) {
                logp[j] = get_log_pq(ev->pmt_hits[i].q,j) + get_log_phit(j) - mu_sum[i] + j*log_mu - lnfact(j);

                if (!charge_only) {
                    if (fast)
                        logp[j] += log_pt_fast;
                    else
                        logp[j] += log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i][0][0], n*3, &ts[i][0][0], ts[i][0][0], &ts_sigma[i][0][0]);
                }

                if (j == 1 || logp[j] > max_logp) max_logp = logp[j];

                if (logp[j] - max_logp < min_ratio*ln(10)) {
                    j++;
                    break;
                }
            }

            /* Now, we include the probability that the channel is miscalibrated as:
             *
             *     P(q,t) = P(q,t|calibrated)P(calibrated) + P(q,t|miscalibrated)P(miscalibrated)
             *
             * For P(q,t|miscalibrated) we just assume a uniform distribution
             * over all 4096 possible values. */
            tmp[0] = logsumexp(logp+1, j-1) + log(1-P_MISCALIBRATION);
            tmp[1] = log(1/4096.0) + log(1/4096.0) + log(P_MISCALIBRATION);
            nll[nhit++] = -logsumexp(tmp, 2);
        } else {
            logp[0] = -mu_sum[i];
            if (fast) {
                nll[nhit++] = -logp[0];
                continue;
            }
            for (j = 1; j < MAX_PE_NO_HIT; j++) {
                logp[j] = get_log_pmiss(j) - mu_sum[i] + j*log_mu - lnfact(j);
            }

            /* Now, we include the probability that the channel is miscalibrated as:
             *
             *     P(not hit) = P(not hit|calibrated)P(calibrated) + P(not hit|miscalibrated)P(miscalibrated)
             *
             * where by miscalibrated here we really mean a PMT that wasn't
             * read out or included in the event for whatever reason. */
            tmp[0] = logsumexp(logp, MAX_PE_NO_HIT) + log(1-P_MISCALIBRATION);
            tmp[1] = log(P_MISCALIBRATION);
            nll[nhit++] = -logsumexp(tmp, 2);
        }
    }

    logp_path = 0.0;
    for (j = 0; j < n; j++)
        for (i = 0; i < v[j].n; i++)
            logp_path += log_norm(v[j].z1[i],0,1) + log_norm(v[j].z2[i],0,1);

    /* Add up the negative log likelihood terms from each PMT. I use a kahan
     * sum here instead of just summing all the values to help prevent round
     * off error. I actually don't think this makes any difference here since
     * the numbers are all of the same order of magnitude, so I should actually
     * test to see if it makes any difference. */
    return kahan_sum(nll,nhit) - logp_path;
}