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"

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;

    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->shower_photons = 0.0;
    p->delta_ray_a = 0.0;
    p->delta_ray_b = 0.0;
    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);

        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 = rad*ELECTRON_PHOTONS_PER_MEV;

        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_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], 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 = rad*MUON_PHOTONS_PER_MEV;

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

        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;

        p->shower_photons = rad*PROTON_PHOTONS_PER_MEV;

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

    return p;
}

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

static double get_expected_charge_shower(particle *p, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double *q_delta_ray, double *q_indirect_delta_ray)
{
    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter, constant;

    SUB(pmt_dir,pmt_pos,pos);

    normalize(pmt_dir);

    cos_theta_pmt = DOT(pmt_dir,pmt_normal);

    charge = 0.0;
    *q_delta_ray = 0.0;
    *q_indirect_delta_ray = 0.0;
    *reflected = 0.0;
    if (cos_theta_pmt > 0) return 0.0;

    /* 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,pmt_pos,pmt_normal,r);

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

    absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
    absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
    absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
    scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o();
    scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o();

    l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;

    constant = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*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. */
    if (p->shower_photons > 0)
        charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/n_d2o);

    if (p->delta_ray_photons > 0)
        *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/n_d2o);

    scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o);

    *reflected = (f_reflec + 1.0 - scatter)*charge;
    *q_indirect_delta_ray = (f_reflec + 1.0 - scatter)*(*q_delta_ray);

    *q_delta_ray *= f*scatter;

    return f*charge*scatter;
}

static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o)
{
    double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter;

    z = 1.0;

    R = NORM(pos);

    n = (R <= AV_RADIUS) ?  n_d2o : n_h2o;

    *reflected = 0.0;
    if (beta < 1/n) return 0.0;

    SUB(pmt_dir,pmt_pos,pos);

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

    *reflected = 0.0;
    if (fabs(cos_theta-1.0/(n_d2o*beta))/theta0 > 5) return 0.0;

    cos_theta_pmt = DOT(pmt_dir,pmt_normal);

    if (cos_theta_pmt > 0) return 0.0;

    omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r);

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

    absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
    absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
    absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
    scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o();
    scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o();

    l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;

    charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);

    scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o);

    *reflected = (f_reflec + 1.0 - scatter)*charge;

    return f*charge*scatter;
}

double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, 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 + PSUP_REFLECTION_TIME)
        p +=  mu_indirect*(t-tmean)/PSUP_REFLECTION_TIME;
    else
        p +=  mu_indirect;

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

    return p/mu_total;
}

double time_pdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
{
    /* 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 past some time `tmean`.
     *
     * 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. */
    size_t i;
    double p, mu_total;

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

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

    return p/mu_total;
}

double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n2, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
{
    /* 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. */
    if (n == 1)
        return ln(n) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
    else
        return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
}

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)
{
    size_t i;
    static double qs[MAX_NPOINTS];
    static double q_indirects[MAX_NPOINTS];
    static double ts[MAX_NPOINTS];
    static double ts2[MAX_NPOINTS];
    double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, q_delta_ray, q_indirect_delta_ray, dx, tmp;

    if (n > MAX_NPOINTS) {
        fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
        exit(1);
    }

    /* FIXME: I just calculate delta assuming 400 nm light. */
    wavelength0 = 400.0;
    n_d2o = get_index_snoman_d2o(wavelength0);
    n_h2o = get_index_snoman_h2o(wavelength0);

    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_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);

        t = x[i]/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;

        qs[i] = get_expected_charge_shower(p, pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, &q_delta_ray, &q_indirect_delta_ray);
        qs[i] = qs[i]*pdf[i] + q_delta_ray/p->range;
        q_indirects[i] = q_indirects[i]*pdf[i] + q_indirect_delta_ray/p->range;
        ts[i] = t*qs[i];
        ts2[i] = t*t*qs[i];
    }

    dx = x[1]-x[0];
    *mu_direct = trapz(qs,dx,n);
    *mu_indirect = trapz(q_indirects,dx,n);

    if (*mu_direct == 0.0) {
        *time = 0.0;
        *sigma = PMT_TTS;
    } else {
        *time = trapz(ts,dx,n)/(*mu_direct);

        /* Variance in the time = E(t^2) - E(t)^2. */
        tmp = trapz(ts2,dx,n)/(*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;
    static double qs[MAX_NPOINTS];
    static double q_indirects[MAX_NPOINTS];
    static double ts[MAX_NPOINTS];
    double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x;

    if (p->len > MAX_NPOINTS) {
        fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
        exit(1);
    }

    /* FIXME: I just calculate delta assuming 400 nm light. */
    wavelength0 = 400.0;
    n_d2o = get_index_snoman_d2o(wavelength0);
    n_h2o = get_index_snoman_h2o(wavelength0);

    for (i = 0; i < p->len; i++) {
        x = p->s[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];

        t = p->t[i];

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

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

        t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;

        qs[i] = get_expected_charge(x, beta, theta0, pos, dir, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o);
        ts[i] = t*qs[i];
    }

    double dx = p->s[1]-p->s[0];
    *mu_direct = trapz(qs,dx,p->len);
    *mu_indirect = trapz(q_indirects,dx,p->len);
    *time = trapz(ts,dx,p->len);
}

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 n_d2o, double n_h2o, 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, theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, f_reflected, charge, prob, frac;

    /* 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/n_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)*n_d2o*x*beta0*theta0;

    theta_pmt = acos(-cos_theta_pmt);
    f = get_weighted_pmt_response(theta_pmt);
    f_reflected = get_weighted_pmt_reflectivity(theta_pmt);

    absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
    absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
    absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();

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

    l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;

    /* Assume the particle is travelling at the speed of light. */
    *t = s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;

    double abs_prob = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic);

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

    /* Add expected number of photons from electromagnetic shower. */
    if (p->shower_photons > 0)
        charge += abs_prob*get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI);
    if (p->delta_ray_photons > 0)
        charge += abs_prob*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/n_d2o)*omega/(2*M_PI);

    *mu_reflected = f_reflected*charge;

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

    pars = (betaRootParams *) params;

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

    /* Calculate total energy */
    E = T + pars->p->mass;
    mom = sqrt(E*E - pars->p->mass*pars->p->mass);
    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-1, 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 n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov)
{
    /* Returns an approximate time at which a PMT is most likely to get hit
     * from Cerenkov light.
     *
     * To do this, we try to compute the distance along the track where the PMT
     * is at the Cerenkov angle. */
    double pmt_dir[3], tmp[3];
    double R, cos_theta, sin_theta, s, l_d2o, l_h2o;

    /* First, we find the point along the track where the PMT is at the
     * Cerenkov angle. */
    SUB(pmt_dir,pmt_pos,pos);
    /* Compute the distance to the PMT. */
    R = NORM(pmt_dir);
    normalize(pmt_dir);

    /* Calculate the cosine of the angle between the track direction and the
     * vector to the PMT at the start of the track. */
    cos_theta = DOT(dir,pmt_dir);
    sin_theta = sqrt(1-pow(cos_theta,2));

    /* Now, we compute the distance along the track where the PMT is at the
     * Cerenkov angle.
     *
     * Note: This formula comes from using the "Law of sines" where the three
     * vertices of the triangle are the starting position of the track, the
     * point along the track that we want to find, and the PMT position. */
    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*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
}

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

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, maxj;
    static double logp[MAX_PE], nll[MAX_PMTS], mu[MAX_PMTS], ts[MAX_PMTS], ts_shower, ts_sigma, mu_shower;
    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
             * come from noise hits. */
            mu[i] = mu_noise;
            continue;
        }

        /* Technically we should be computing:
         *
         *     mu[i] = max(p(q|mu))
         *
         * where by max I mean we find the expected number of PE which
         * maximizes p(q|mu). However, I don't know of any easy analytical
         * solution to this. So instead we just compute the number of PE which
         * maximizes p(q|n). */
        for (j = 1; j < MAX_PE; j++) {
            logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j);

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

        mu[i] = maxj;
        ts[i] = ev->pmt_hits[i].t;
    }

    mu_shower = 0;
    ts_shower = 0;
    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].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], &mu_shower, 1, &ts[i], &ts_shower, ts[i], PMT_TTS, &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);
}

double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, 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, nhit;
    static double logp[MAX_PE], nll[MAX_PMTS];
    double range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
    particle *p;
    double wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov;
    double logp_path;
    double mu_reflected;
    path *path;

    static double mu_direct[MAX_PMTS][MAX_VERTICES];
    static double mu_shower[MAX_PMTS][MAX_VERTICES];
    static double ts[MAX_PMTS][MAX_VERTICES];
    static double ts_shower[MAX_PMTS][MAX_VERTICES];
    static double ts_sigma[MAX_PMTS][MAX_VERTICES];
    static double mu[MAX_PMTS];
    double mu_noise, mu_indirect[MAX_VERTICES], mu_indirect_total;

    double *x, *pdf;

    double result;

    double min_ratio;

    double pos_a, pos_b;

    double xlo, xhi;

    size_t npoints, npoints_shower;

    double log_pt_fast;

    if (n > MAX_VERTICES) {
        fprintf(stderr, "maximum number of vertices is %i\n", MAX_VERTICES);
        exit(1);
    }

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

    for (j = 0; j < n; j++) {
        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;

        /* 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);
        beta0 = p0/E0;

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

            x = malloc(sizeof(double)*npoints_shower);
            pdf = malloc(sizeof(double)*npoints_shower);
            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);
            }

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

        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. */
        wavelength0 = 400.0;
        n_d2o = get_index_snoman_d2o(wavelength0);
        n_h2o = get_index_snoman_h2o(wavelength0);
        cos_theta_cerenkov = 1/(n_d2o*beta0);
        sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2));

        mu_indirect[j] = 0.0;
        for (i = 0; i < MAX_PMTS; i++) {
            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_direct[i][j] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov);
                ts[i][j] += v[j].t0;
                mu_indirect[j] += mu_reflected;
                continue;
            }

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

            if (mu_direct[i][j] > 1e-9) {
                ts[i][j] = v[j].t0 + result/mu_direct[i][j];
            } 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] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov);
            }

            integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j]);
            mu_indirect[j] += q_indirect;
            ts_shower[i][j] = v[j].t0 + result;
        }

        if (!fast) {
            path_free(path);

            free(x);
            free(pdf);
        }

        particle_free(p);
    }

    mu_noise = DARK_RATE*GTVALID*1e-9;

    mu_indirect_total = 0.0;
    for (j = 0; j < n; j++) {
        if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS)
            mu_indirect_total += mu_indirect[j]*CHARGE_FRACTION_ELECTRON/10000.0;
        else
            mu_indirect_total += mu_indirect[j]*CHARGE_FRACTION_MUON/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;

        for (j = 0; j < n; j++) {
            if (fast)
                mu[i] += mu_direct[i][j] + mu_noise;
            else
                mu[i] += mu_direct[i][j] + mu_shower[i][j] + mu_indirect_total + mu_noise;
        }
    }

    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[i]);

        if (fast)
            log_pt_fast = log_pt(ev->pmt_hits[i].t, 1, mu_noise, mu_indirect_total, &mu_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]);

        if (ev->pmt_hits[i].hit) {
            for (j = 1; j < MAX_PE; j++) {
                logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[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_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]);
                }

                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];
            if (fast) {
                nll[nhit++] = -logp[0];
                continue;
            }
            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);
        }
    }

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