path: root/src/likelihood.c

#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"

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;

    particle *p = malloc(sizeof(particle));
    p->x = malloc(sizeof(double)*n);
    p->T = malloc(sizeof(double)*n);
    p->n = n;

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

        /* Now we loop over the points along the track and assume dE/dx is
         * constant between points. */
        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]);
        }
        /* 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;

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

        /* Now we loop over the points along the track and assume dE/dx is
         * constant between points. */
        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]);
        }
        /* 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;

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

        /* Now we loop over the points along the track and assume dE/dx is
         * constant between points. */
        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]);
        }
        /* 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;

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

    z = 1.0;

    SUB(pmt_dir,pmt_pos,pos);

    normalize(pmt_dir);

    cos_theta_pmt = DOT(pmt_dir,pmt_normal);

    *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_approx(pos,pmt_pos,pmt_normal,r);

    R = NORM(pos);

    if (R <= AV_RADIUS) {
        n = n_d2o;
    } else {
        n = n_h2o;
    }

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

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

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

    *reflected = f_reflec*charge;

    return f*charge;
}

double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double 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 + mu_indirect*(pow(sigma,2)*norm(tmean,t,sigma) + (t-tmean)*norm_cdf(t,tmean,sigma))/(GTVALID-tmean);

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

    return p/mu_total;
}

double f(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double 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*norm_cdf(t,tmean,sigma)/(GTVALID-tmean);

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

    return p/mu_total;
}

double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu_direct, size_t n2, double *ts, double tmean, double 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. */
    return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma));
}

static void integrate_path(path *p, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time)
{
    size_t i;
    double *qs, *q_indirects, *ts;
    double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x;

    qs = malloc(sizeof(double)*n);
    q_indirects = malloc(sizeof(double)*n);
    ts = malloc(sizeof(double)*n);

    /* 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++) {
        x = a + i*(b-a)/(n-1);

        path_eval(p, x, pos, dir, &beta, &t, &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];
    }

    *mu_direct = trapz(qs,(b-a)/(n-1),n);
    *mu_indirect = trapz(q_indirects,(b-a)/(n-1),n);
    *time = trapz(ts,(b-a)/(n-1),n);

    free(qs);
    free(q_indirects);
    free(ts);
}

static double get_total_charge_approx(double T0, 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, E0, p0, beta0, f, cos_theta_pmt, theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0, f_reflected, charge;

    /* Calculate beta at the start of the track. */
    E0 = T0 + p->mass;
    p0 = sqrt(E0*E0 - p->mass*p->mass);
    beta0 = p0/E0;

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

    /* Calculate the cosine of the angle between the track direction and the
     * vector to the PMT at the start of the track. */
    cos_theta = DOT(dir,pmt_dir);
    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. */
    s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov;

    /* Make sure that the point is somewhere along the track between 0 and
     * `smax`. */
    if (s < 0) s = 0.0;
    else if (s > smax) s = smax;

    /* Compute the vector from the point `s` along the track to the PMT. */
    tmp[0] = pmts[i].pos[0] - (pos[0] + s*dir[0]);
    tmp[1] = pmts[i].pos[1] - (pos[1] + s*dir[1]);
    tmp[2] = pmts[i].pos[2] - (pos[2] + s*dir[2]);

    /* To do the integral analytically, we expand the distance to the PMT along
     * the track in a Taylor series around `s0`, i.e.
     *
     *     r(s) = a + b*(s-s0)
     *
     * Here, we calculate `a` which is the distance to the PMT at the point
     * `s`. */
    a = NORM(tmp);

    /* `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. */
    double prob = get_probability2(beta);

    /* Compute the position of the particle at a distance `s` along the track. */
    tmp[0] = pos[0] + s*dir[0];
    tmp[1] = pos[1] + s*dir[1];
    tmp[2] = pos[2] + s*dir[2];

    SUB(pmt_dir,pmts[i].pos,tmp);

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

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

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

    wavelength0 = 400.0;
    absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0);
    absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0);
    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;

    charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*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);

    *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 getKineticEnergy(double x, void *p)
{
    return particle_get_energy(x, (particle *) p);
}

double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, size_t npoints, int fast)
{
    size_t i, j, nhit;
    double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
    particle *p;
    double pmt_dir[3], R, cos_theta, sin_theta, wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov, s, l_h2o, l_d2o;
    double logp_path;
    double a, b;
    double tmp[3];
    double mu_reflected;
    path *path;

    double mu_direct[MAX_PMTS];
    double ts[MAX_PMTS];
    double mu[MAX_PMTS];
    double mu_noise, mu_indirect;

    double result;

    double min_ratio;

    p = particle_init(id, T0, 10000);
    range = p->range;

    /* Calculate total energy */
    E0 = 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(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, 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. */
    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 = 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 (fast && !ev->pmt_hits[i].hit) continue;

        if (fast) {
            mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov);
            ts[i] += t0;
            mu_indirect += mu_reflected;
        } else {
            /* First, we try to compute the distance along the track where the
             * PMT is at the Cerenkov angle. The reason for this is because for
             * heavy particles like muons which don't scatter much, the
             * probability distribution for getting a photon hit along the
             * track looks kind of like a delta function, i.e. the PMT is only
             * hit over a very narrow window when the angle between the track
             * direction and the PMT is *very* close to the Cerenkov angle
             * (it's not a perfect delta function since there is some width due
             * to dispersion). In this case, it's possible that the numerical
             * integration completely skips over the delta function and so
             * predicts an expected charge of 0.
             *
             * To fix this, we only integrate 1 meter on both sides of the
             * point along the track where the PMT is at the Cerenkov angle. */

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

            /* Make sure that the point is somewhere along the track between 0 and
             * `smax`. */
            if (s < 0) s = 0.0;
            else if (s > smax) s = smax;

            /* Compute the endpoints of the integration. */
            a = s - 100.0;
            b = s + 100.0;

            if (a < 0.0) a = 0.0;
            if (b > range) b = range;

            double q_indirect;
            integrate_path(path, i, a, b, npoints, mu_direct+i, &q_indirect, &result);
            mu_indirect += q_indirect;

            if (mu_direct[i] > 1e-9) {
                ts[i] = t0 + result/mu_direct[i];
            } 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. */
                n_h2o = get_index_snoman_h2o(wavelength0);

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

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

                /* Assume the particle is travelling at the speed of light. */
                ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
            }
        }
    }

    if (!fast)
        path_free(path);

    particle_free(p);

    mu_noise = DARK_RATE*GTVALID*1e-9;
    mu_indirect *= CHARGE_FRACTION/10000.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 (fast && !ev->pmt_hits[i].hit) continue;

        if (fast)
            mu[i] = mu_direct[i] + mu_noise;
        else
            mu[i] = mu_direct[i] + mu_indirect + mu_noise;
    }

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

    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 (fast && !ev->pmt_hits[i].hit) continue;

        log_mu = log(mu[i]);

        if (ev->pmt_hits[i].hit) {
            for (j = 1; j < MAX_PE; j++) {
                logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], ts[i], PMT_TTS);
                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 (i = 0; i < n; i++)
        logp_path += log_norm(z1[i],0,1) + log_norm(z2[i],0,1);

    return kahan_sum(nll,nhit) - logp_path;
}