fix delta ray charge calculation

Previously I was calculating the expected number of delta ray photons when
integrating over the shower path, but since the delta rays are produced along
the particle path and not further out like the shower photons, this wasn't
correct. The normalization of the probability distribution for the photons
produced along the path was also not handled correctly.

This commit adds a new function called integrate_path_delta_ray() to compute
the expected number of photons from delta rays hitting each PMT. Currently this
means that the likelihood function for muons will be significantly slower than
previously, but hopefully I can speed it up again in the future (for example by
skipping the shower calculation which is negligible for lower energy muons).

1 files changed, 193 insertions, 76 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 0228bd8..1919a53 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -208,7 +208,7 @@ void particle_free(particle *p)
     free(p);
 }
 
-static void get_expected_charge_shower(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *q_delta_ray, double *q_indirect_delta_ray, double *t)
+static void get_expected_charge_delta_ray(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
 {
     double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o;
 
@@ -221,8 +221,6 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
     *t = 0.0;
     *q = 0.0;
     *reflected = 0.0;
-    *q_delta_ray = 0.0;
-    *q_indirect_delta_ray = 0.0;
     if (cos_theta_pmt > 0) return;
 
     /* Calculate the cosine of the angle between the track direction and the
@@ -247,17 +245,55 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
      * showers from lower energy electrons or for delta rays from lower energy
      * muons. It seems good enough, but in the future it would be nice to
      * parameterize this. */
-    charge = 0.0;
-    if (p->shower_photons > 0)
-        charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/avg_index_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/avg_index_d2o);
+    charge = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/avg_index_d2o);
 
     *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
-    *q_indirect_delta_ray = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*(*q_delta_ray) + prob_sct*(*q_delta_ray);
 
-    *q_delta_ray *= (1.0-prob_abs)*(1.0-prob_sct)*f;
+    *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+
+    *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
+}
+
+static void get_expected_charge_shower(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
+{
+    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o;
+
+    SUB(pmt_dir,pmts[pmt].pos,pos);
+
+    normalize(pmt_dir);
+
+    cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);
+
+    *t = 0.0;
+    *q = 0.0;
+    *reflected = 0.0;
+    if (cos_theta_pmt > 0) return;
+
+    /* Calculate the cosine of the angle between the track direction and the
+     * vector to the PMT. */
+    cos_theta = DOT(dir,pmt_dir);
+
+    omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);
+
+    f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
+    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);
+
+    constant = get_weighted_quantum_efficiency()*omega/(2*M_PI);
+
+    /* Note: We assume here that the peak of the angular distribution is at the
+     * Cerenkov angle for a particle with beta = 1. This seems to be an OK
+     * assumption for high energy showers, but is not exactly correct for
+     * showers from lower energy electrons or for delta rays from lower energy
+     * muons. It seems good enough, but in the future it would be nice to
+     * parameterize this. */
+    charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/avg_index_d2o);
+
+    *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
 
     *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
 
@@ -327,7 +363,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
     *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
 }
 
-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)
+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;
@@ -344,20 +380,15 @@ double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct
 
     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];
+        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;
 }
 
-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)
+double time_pdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma)
 {
     /* Returns the probability that a photon is detected at time `t`.
      *
@@ -376,24 +407,19 @@ double time_pdf(double t, double mu_noise, double mu_indirect, double *mu_direct
     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);
+    p = mu_noise/GTVALID + mu_indirect*(erf((t-tmean)/(sqrt(2)*PMT_TTS))-erf((t-tmean-PSUP_REFLECTION_TIME)/(sqrt(2)*PMT_TTS)))/(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];
+        if (mu[i] == 0.0) continue;
+        p += mu[i]*norm(t,ts[i],ts_sigma[i]);
+        mu_total += mu[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)
+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)
 {
     /* Returns the first order statistic for observing a PMT hit at time `t`
      * given `n` hits.
@@ -406,15 +432,71 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
      * 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));
+        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_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));
+        return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma));
+}
+
+static void integrate_path_delta_ray(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
+{
+    size_t i;
+    double pos[3], t, q, r, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
+
+    q_sum = 0.0;
+    r_sum = 0.0;
+    t_sum = 0.0;
+    t2_sum = 0.0;
+    for (i = 0; i < n; i++) {
+        pos[0] = pos0[0] + x[i]*dir0[0];
+        pos[1] = pos0[1] + x[i]*dir0[1];
+        pos[2] = pos0[2] + x[i]*dir0[2];
+
+        get_expected_charge_delta_ray(p, pos, dir0, pmt, &q, &r, &t);
+
+        q = q*pdf[i];
+        r = r*pdf[i];
+        t += x[i]/SPEED_OF_LIGHT;
+
+        if (i == 0 || i == (n - 1)) {
+            q_sum += q;
+            r_sum += r;
+            t_sum += t*q;
+            t2_sum += t*t*q;
+        } else {
+            q_sum += 2*q;
+            r_sum += 2*r;
+            t_sum += 2*t*q;
+            t2_sum += 2*t*t*q;
+        }
+    }
+
+    dx = x[1] - x[0];
+    *mu_direct = q_sum*dx*0.5;
+    *mu_indirect = r_sum*dx*0.5;
+
+    if (*mu_direct == 0.0) {
+        *time = 0.0;
+        *sigma = PMT_TTS;
+    } else {
+        *time = t_sum*dx*0.5/(*mu_direct);
+
+        /* Variance in the time = E(t^2) - E(t)^2. */
+        tmp = t2_sum*dx*0.5/(*mu_direct) - (*time)*(*time);
+
+        if (tmp >= 0) {
+            *sigma = sqrt(tmp + PMT_TTS*PMT_TTS);
+        } else {
+            /* This should never happen but does occasionally due to floating
+             * point rounding error. */
+            *sigma = PMT_TTS;
+        }
+    }
 }
 
 static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
 {
     size_t i;
-    double pos[3], t, q, r, qd, rd, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
+    double pos[3], t, q, r, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
 
     q_sum = 0.0;
     r_sum = 0.0;
@@ -425,10 +507,10 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
         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, &qd, &rd, &t);
+        get_expected_charge_shower(p, pos, dir0, pmt, &q, &r, &t);
 
-        q = q*pdf[i] + qd/p->range;
-        r = r*pdf[i] + rd/p->range;
+        q = q*pdf[i];
+        r = r*pdf[i];
         t += x[i]/SPEED_OF_LIGHT;
 
         if (i == 0 || i == (n - 1)) {
@@ -803,7 +885,7 @@ double nll_best(event *ev)
      * 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;
+    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;
@@ -840,8 +922,6 @@ double nll_best(event *ev)
         ts[i] = ev->pmt_hits[i].t;
     }
 
-    mu_shower = 0;
-    ts_shower = 0;
     ts_sigma = PMT_TTS;
 
     mu_indirect_total = 0.0;
@@ -861,7 +941,7 @@ double nll_best(event *ev)
 
         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);
+                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], 1, &ts[i], ts[i], &ts_sigma);
 
                 if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
 
@@ -897,7 +977,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
      * `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;
+    size_t i, j, k, nhit;
     static double logp[MAX_PE], nll[MAX_PMTS];
     double range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
     particle *p;
@@ -906,15 +986,29 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
     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];
+    /* 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;
 
-    double *x, *pdf;
+    static double x[MAX_NPOINTS], pdf[MAX_NPOINTS];
+    static double x_delta[MAX_NPOINTS], pdf_delta[MAX_NPOINTS];
 
     double result;
 
@@ -935,13 +1029,13 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
 
     /* Initialize static arrays. */
     for (i = 0; i < MAX_PMTS; i++) {
-        mu[i] = 0.0;
+        mu_sum[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 (k = 0; k < 3; k++) {
+                mu[i][j][k] = 0.0;
+                ts[i][j][k] = 0.0;
+                ts_sigma[i][j][k] = PMT_TTS;
+            }
         }
     }
 
@@ -996,8 +1090,6 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
             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);
@@ -1008,6 +1100,23 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
             for (i = 0; i < npoints_shower; i++) {
                 pdf[i] /= norm;
             }
+
+            for (i = 0; i < npoints_shower; i++) {
+                x_delta[i] = i*p->range/(npoints_shower-1);
+                /* For now, we assume the delta rays are produced uniformly
+                 * along the path. This isn't quite correct, but the
+                 * distribution is close to flat.
+                 *
+                 * FIXME: Figure out a good parameterization for the
+                 * distribution of delta ray photons along the path. */
+                pdf_delta[i] = 1.0;
+            }
+
+            /* Make sure the PDF is normalized. */
+            norm = trapz(pdf_delta,x_delta[1]-x_delta[0],npoints_shower);
+            for (i = 0; i < npoints_shower; i++) {
+                pdf_delta[i] /= norm;
+            }
         }
 
         if (beta0 > BETA_MIN)
@@ -1031,18 +1140,18 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
             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, cos_theta_cerenkov, sin_theta_cerenkov);
-                ts[i][j] += v[j].t0;
+                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;
             }
 
             double q_indirect;
-            integrate_path(path, i, &mu_direct[i][j], &q_indirect, &result);
+            integrate_path(path, i, &mu[i][j][0], &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];
+            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
@@ -1054,19 +1163,24 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
                  * 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,cos_theta_cerenkov,sin_theta_cerenkov);
+                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);
             }
 
-            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 (p->shower_photons > 0) {
+                integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu[i][j][1],&q_indirect,&result,&ts_sigma[i][j][1]);
+                mu_indirect[j] += q_indirect;
+                ts[i][j][1] = v[j].t0 + result;
+            }
+
+            if (p->delta_ray_photons > 0) {
+                integrate_path_delta_ray(p,x_delta,pdf_delta,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu[i][j][2],&q_indirect,&result,&ts_sigma[i][j][2]);
+                mu_indirect[j] += q_indirect;
+                ts[i][j][2] = v[j].t0 + result;
+            }
         }
 
         if (!fast) {
             path_free(path);
-
-            free(x);
-            free(pdf);
         }
 
         particle_free(p);
@@ -1089,8 +1203,11 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
          * calculation. */
         if (hit_only && !ev->pmt_hits[i].hit) continue;
 
+        mu_sum[i] = mu_indirect_total + mu_noise;
         for (j = 0; j < n; j++) {
-            mu[i] += mu_direct[i][j] + mu_shower[i][j] + mu_indirect_total + mu_noise;
+            for (k = 0; k < 3; k++) {
+                mu_sum[i] += mu[i][j][k];
+            }
         }
     }
 
@@ -1108,20 +1225,20 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
          * calculation. */
         if (hit_only && !ev->pmt_hits[i].hit) continue;
 
-        log_mu = log(mu[i]);
+        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_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]);
+            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].qhs,j) - mu[i] + j*log_mu - lnfact(j);
+                logp[j] = get_log_pq(ev->pmt_hits[i].qhs,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_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]);
+                        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];
@@ -1134,13 +1251,13 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
 
             nll[nhit++] = -logsumexp(logp+1, j-1);
         } else {
-            logp[0] = -mu[i];
+            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[i] + j*log_mu - lnfact(j);
+                logp[j] = get_log_pmiss(j) - mu_sum[i] + j*log_mu - lnfact(j);
             }
             nll[nhit++] = -logsumexp(logp, MAX_PE_NO_HIT);
         }