1 files changed, 97 insertions, 52 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 786b221..8e12f57 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -42,16 +42,19 @@ particle *particle_init(int id, double T0, size_t n)
      * 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;
+    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->rad = 0.0;
     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;
@@ -67,14 +70,17 @@ particle *particle_init(int id, double T0, size_t n)
         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);
-            p->rad += dEdx*(p->x[i]-p->x[i-1]);
+            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
@@ -85,6 +91,8 @@ particle *particle_init(int id, double T0, size_t n)
          * 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:
@@ -96,14 +104,21 @@ particle *particle_init(int id, double T0, size_t n)
          * 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);
-            p->rad += dEdx*(p->x[i]-p->x[i-1]);
+            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 +129,8 @@ particle *particle_init(int id, double T0, size_t n)
          * 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;
@@ -121,14 +138,17 @@ particle *particle_init(int id, double T0, size_t n)
          * 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);
-            p->rad += dEdx*(p->x[i]-p->x[i-1]);
+            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
@@ -139,6 +159,8 @@ particle *particle_init(int id, double T0, size_t n)
          * 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);
@@ -170,9 +192,9 @@ 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)
+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;
+    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);
 
@@ -180,6 +202,9 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p
 
     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;
 
@@ -201,11 +226,26 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p
 
     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);
+    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;
 }
@@ -345,14 +385,14 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
         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, double rad)
+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;
+    double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, q_delta_ray, q_indirect_delta_ray;
 
     if (n > MAX_NPOINTS) {
         fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
@@ -364,6 +404,8 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
     n_d2o = get_index_snoman_d2o(wavelength0);
     n_h2o = get_index_snoman_h2o(wavelength0);
 
+    double range = x[n-1] - x[0];
+
     for (i = 0; i < n; i++) {
         pos[0] = pos0[0] + x[i]*dir0[0];
         pos[1] = pos0[1] + x[i]*dir0[1];
@@ -373,8 +415,9 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
 
         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, p->a, p->b, rad)*pdf[i];
-        q_indirects[i] *= pdf[i];
+        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/range;
+        q_indirects[i] = q_indirects[i]*pdf[i] + q_indirect_delta_ray/range;
         ts[i] = t*qs[i];
         ts2[i] = t*t*qs[i];
     }
@@ -593,10 +636,11 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
 
     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 for electrons
-     * and positrons. */
-    if (p->id == IDP_E_MINUS || p->id == IDP_E_PLUS)
-        charge += abs_prob*get_weighted_quantum_efficiency()*p->rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*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;
 
@@ -891,6 +935,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
         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;
@@ -903,32 +948,43 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
         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);
 
-            if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) {
+            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);
+            /* 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);
+            /* 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 < 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);
-                }
+            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;
-                }
+            /* 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;
             }
         }
 
@@ -982,27 +1038,16 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
                 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);
             }
 
-            if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) {
-                /* If we are fitting for an electron or positron, we need to
-                 * estimate the expected number of photons produced by the
-                 * electromagnetic shower.
-                 *
-                 * Technically we should also have a term for the electromagnetic
-                 * shower produced by muons, but for the energy range I'm
-                 * considering right now, it is very small. */
-                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],p->rad);
-                mu_indirect[j] += q_indirect;
-                ts_shower[i][j] = v[j].t0 + result;
-            }
+            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);
 
-            if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) {
-                free(x);
-                free(pdf);
-            }
+            free(x);
+            free(pdf);
         }
 
         particle_free(p);