update likelihood function to fit electrons!

To characterize the angular distribution of photons from an electromagnetic
shower I came up with the following functional form:

    f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)

and fit this to data simulated using RAT-PAC at several different energies. I
then fit the alpha and beta coefficients as a function of energy to the
functional form:

    alpha = c0 + c1/log(c2*T0 + c3)
    beta  = c0 + c1/log(c2*T0 + c3).

where T0 is the initial energy of the electron in MeV and c0, c1, c2, and c3
are parameters which I fit.

The longitudinal distribution of the photons generated from an electromagnetic
shower is described by a gamma distribution:

    f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a).

This parameterization comes from the PDG "Passage of particles through matter"
section 32.5. I also fit the data from my RAT-PAC simulation, but currently I
am not using it, and instead using a simpler form to calculate the coefficients
from the PDG (although I estimated the b parameter from the RAT-PAC data).

I also sped up the calculation of the solid angle by making a lookup table
since it was taking a significant fraction of the time to compute the
likelihood function.

1 files changed, 230 insertions, 88 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 2340fb1..f769216 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -21,6 +21,7 @@
 #include "electron.h"
 #include "proton.h"
 #include "id_particles.h"
+#include <gsl/gsl_cdf.h>
 
 particle *particle_init(int id, double T0, size_t n)
 {
@@ -46,6 +47,7 @@ particle *particle_init(int id, double T0, size_t n)
     p->x = malloc(sizeof(double)*n);
     p->T = malloc(sizeof(double)*n);
     p->n = n;
+    p->rad = 0.0;
 
     p->x[0] = 0;
     p->T[0] = T0;
@@ -64,6 +66,8 @@ particle *particle_init(int id, double T0, size_t n)
             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);
+            p->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
@@ -91,6 +95,8 @@ particle *particle_init(int id, double T0, size_t n)
             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);
+            p->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
@@ -114,6 +120,8 @@ particle *particle_init(int id, double T0, size_t n)
             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);
+            p->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
@@ -155,6 +163,50 @@ void particle_free(particle *p)
     free(p);
 }
 
+static double get_expected_charge_shower(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 alpha, double beta, double rad)
+{
+    double pmt_dir[3], cos_theta, n, omega, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge;
+
+    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_fast(pos,pmt_pos,pmt_normal,r);
+
+    R = NORM(pos);
+
+    if (R <= AV_RADIUS) {
+        n = n_d2o;
+    } else {
+        n = n_h2o;
+    }
+
+    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)*get_weighted_quantum_efficiency()*rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,alpha,beta,1/n_d2o)*omega/(2*M_PI);
+
+    *reflected = f_reflec*charge;
+
+    return f*charge;
+}
+
 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;
@@ -174,7 +226,7 @@ static double get_expected_charge(double x, double beta, double theta0, double *
      * vector to the PMT. */
     cos_theta = DOT(dir,pmt_dir);
 
-    omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r);
+    omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r);
 
     R = NORM(pos);
 
@@ -204,7 +256,7 @@ static double get_expected_charge(double x, double beta, double theta0, double *
     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)
+double F(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;
@@ -217,11 +269,15 @@ double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_
         p += mu_direct[i]*norm_cdf(t,ts[i],sigma);
         mu_total += mu_direct[i];
     }
+    for (i = 0; i < n; i++) {
+        p += mu_shower[i]*norm_cdf(t,ts_shower[i],ts_sigma[i]);
+        mu_total += mu_shower[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)
+double f(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`.
      *
@@ -247,11 +303,15 @@ double f(double t, double mu_noise, double mu_indirect, double *mu_direct, size_
         p += mu_direct[i]*norm(t,ts[i],sigma);
         mu_total += mu_direct[i];
     }
+    for (i = 0; i < n; i++) {
+        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, size_t n2, double *ts, double tmean, double sigma)
+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.
@@ -263,18 +323,73 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
      * 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));
+    return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
+}
+
+static void integrate_path_shower(double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad, double pos_a, double pos_b)
+{
+    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, alpha, beta;
+
+    if (n > MAX_NPOINTS) {
+        fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
+        exit(1);
+    }
+
+    alpha = electron_get_angular_distribution_alpha(T0);
+    beta = electron_get_angular_distribution_beta(T0);
+
+    /* 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(pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, alpha, beta, rad)*pdf[i];
+        q_indirects[i] *= pdf[i];
+        ts[i] = t*qs[i];
+        ts2[i] = t*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)/(*mu_direct);
+    double t2 = trapz(ts2,(b-a)/(n-1),n)/(*mu_direct);
+
+    *sigma = sqrt(t2-(*time)*(*time));
+
+    if (*mu_direct == 0.0) {
+        *time = 0.0;
+        *sigma = PMT_TTS;
+    } else {
+        *sigma = sqrt((*sigma)*(*sigma) + PMT_TTS*PMT_TTS);
+    }
 }
 
 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;
+    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;
 
-    qs = malloc(sizeof(double)*n);
-    q_indirects = malloc(sizeof(double)*n);
-    ts = malloc(sizeof(double)*n);
+    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;
@@ -297,10 +412,6 @@ static void integrate_path(path *p, int pmt, double a, double b, size_t n, doubl
     *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)
@@ -424,7 +535,7 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti
     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_approx(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS);
+    omega = get_solid_angle_fast(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS);
 
     theta0 = fmax(theta0*sqrt(s),MIN_THETA0);
 
@@ -526,7 +637,8 @@ static double getKineticEnergy(double x, void *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;
+    static double logp[MAX_PE], nll[MAX_PMTS];
+    double 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;
@@ -535,11 +647,31 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
     double mu_reflected;
     path *path;
 
-    double mu_direct[MAX_PMTS];
-    double ts[MAX_PMTS];
-    double mu[MAX_PMTS];
+    static double mu_direct[MAX_PMTS];
+    static double mu_shower[MAX_PMTS] = {0};
+    static double ts[MAX_PMTS];
+    static double ts_shower[MAX_PMTS] = {0};
+    static double ts_sigma[MAX_PMTS] = {0};
+    static double mu[MAX_PMTS];
     double mu_noise, mu_indirect;
 
+    for (i = 0; i < MAX_PMTS; i++) ts_sigma[i] = PMT_TTS;
+
+    double pos_a, pos_b;
+
+    electron_get_position_distribution_parameters(T0, &pos_a, &pos_b);
+
+    /* Upper limit to integrate for the electromagnetic shower. */
+    double xlo = 0.0;
+    double xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b);
+
+    double *x = malloc(sizeof(double)*npoints);
+    double *pdf = malloc(sizeof(double)*npoints);
+    for (i = 0; i < npoints; i++) {
+        x[i] = xlo + i*(xhi-xlo)/(npoints-1);
+        pdf[i] = gamma_pdf(x[i],pos_a,pos_b);
+    }
+
     double result;
 
     double min_ratio;
@@ -582,80 +714,87 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
             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;
+            continue;
+        }
+
+        /* 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 {
-            /* 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 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);
 
-            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];
 
-                /* 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);
 
-                SUB(pmt_dir,pmts[i].pos,tmp);
+            get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);
 
-                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;
+        }
 
-                /* 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 (id == IDP_E_MINUS || id == IDP_E_PLUS) {
+            integrate_path_shower(x,pdf,T0,pos,dir,i,0,xhi,npoints,mu_shower+i,&q_indirect,&result,&ts_sigma[i],p->rad,pos_a,pos_b);
+            mu_indirect += q_indirect;
+            ts_shower[i] = t0 + result;
         }
     }
 
@@ -677,7 +816,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
         if (fast)
             mu[i] = mu_direct[i] + mu_noise;
         else
-            mu[i] = mu_direct[i] + mu_indirect + mu_noise;
+            mu[i] = mu_direct[i] + mu_shower[i] + mu_indirect + mu_noise;
     }
 
     if (fast)
@@ -697,7 +836,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
 
         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);
+                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], &mu_shower[i], 1, &ts[i], &ts_shower[i], ts[i], PMT_TTS, &ts_sigma[i]);
                 if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
 
                 if (logp[j] - max_logp < min_ratio*ln(10)) {
@@ -724,5 +863,8 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
     for (i = 0; i < n; i++)
         logp_path += log_norm(z1[i],0,1) + log_norm(z2[i],0,1);
 
+    free(x);
+    free(pdf);
+
     return kahan_sum(nll,nhit) - logp_path;
 }