add shower photons to fast likelihood calculation

1 files changed, 18 insertions, 8 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index a24fb41..724cc94 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -44,10 +44,13 @@ particle *particle_init(int id, double T0, size_t n)
     double dEdx;
 
     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->x[0] = 0;
     p->T[0] = T0;
@@ -60,6 +63,9 @@ particle *particle_init(int id, double T0, size_t n)
          * 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);
+
         /* Now we loop over the points along the track and assume dE/dx is
          * constant between points. */
         for (i = 1; i < n; i++) {
@@ -320,23 +326,20 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
     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, 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, double rad)
 {
     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;
+    double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o;
 
     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);
@@ -351,7 +354,7 @@ static void integrate_path_shower(double *x, double *pdf, double T0, double *pos
 
         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];
+        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];
         ts[i] = t*qs[i];
         ts2[i] = t*t*qs[i];
@@ -570,7 +573,14 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti
     /* 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);
+    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 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);
 
     *mu_reflected = f_reflected*charge;
 
@@ -841,7 +851,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
              * 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(x,pdf,T0,pos,dir,i,npoints_shower,mu_shower+i,&q_indirect,&result,&ts_sigma[i],p->rad);
+            integrate_path_shower(p,x,pdf,T0,pos,dir,i,npoints_shower,mu_shower+i,&q_indirect,&result,&ts_sigma[i],p->rad);
             mu_indirect += q_indirect;
             ts_shower[i] = t0 + result;
         }