use the PMT response table to calculate the amount of reflected light

To calculate the expected number of photons from reflected light we now
integrate over the track and use the PMT response table to calculate what
fraction of the light is reflected. Previously we were just using a constant
fraction of the total detected light which was faster since we only had to
integrate over the track once, but this should be more accurate.

4 files changed, 82 insertions, 22 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 000389b..daad9cf 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -105,6 +105,16 @@ static double gsl_muon_time(double x, void *params)
     return t*get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
 }
 
+static double gsl_muon_reflected_charge(double x, void *params)
+{
+    intParams *pars = (intParams *) params;
+    double dir[3], pos[3], T, t, theta0;
+
+    path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+
+    return get_expected_reflected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
+}
+
 static double gsl_muon_charge(double x, void *params)
 {
     intParams *pars = (intParams *) params;
@@ -115,7 +125,7 @@ static double gsl_muon_charge(double x, void *params)
     return get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
 }
 
-double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy *m, int i, double smax, double theta0, double *t)
+double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy *m, int i, double smax, double theta0, double *t, double *mu_reflected)
 {
     /* Returns the approximate expected number of photons seen by PMT `i` using
      * an analytic formula.
@@ -146,7 +156,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
      *
      * `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, p, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0;
+    double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, p, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0, f_reflected;
 
     /* Calculate beta at the start of the track. */
     E0 = T0 + MUON_MASS;
@@ -244,6 +254,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
     double frac = sqrt(2)*n_d2o*x*beta0*theta0;
 
     f = get_weighted_pmt_response(acos(-cos_theta_pmt));
+    f_reflected = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt));
 
     absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0);
     absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0);
@@ -256,7 +267,11 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
     /* 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;
 
-    return f*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 charge = f*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 {
@@ -333,14 +348,13 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
 {
     size_t i, j, nhit;
     intParams params;
-    double total_charge;
     double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
-    double tmean = 0.0;
     muon_energy *m;
     double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, n_h2o, theta_cerenkov, s, l_h2o, l_d2o;
     double logp_path;
     double a, b;
     double tmp[3];
+    double mu_reflected;
 
     double mu_direct[MAX_PMTS];
     double ts[MAX_PMTS];
@@ -374,15 +388,16 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
     else
         smax = 0.0;
 
-    total_charge = 0.0;
+    mu_indirect = 0.0;
     for (i = 0; i < MAX_PMTS; i++) {
         if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
 
         params.i = i;
 
         if (fast) {
-            mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i]);
+            mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i], &mu_reflected);
             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
@@ -442,6 +457,11 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
             gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
             mu_direct[i] = result;
 
+            F.function = &gsl_muon_reflected_charge;
+
+            gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
+            mu_indirect += result;
+
             F.function = &gsl_muon_time;
 
             if (mu_direct[i] > 1e-9) {
@@ -468,23 +488,16 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
                 ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
             }
         }
-
-        tmean += ts[i]*mu_direct[i];
-
-        total_charge += mu_direct[i];
     }
 
     path_free(params.p);
 
     muon_free_energy(m);
 
-    if (total_charge > 0)
-        tmean /= total_charge;
-
     gsl_integration_cquad_workspace_free(w);
 
     mu_noise = DARK_RATE*GTVALID*1e-9;
-    mu_indirect = total_charge*CHARGE_FRACTION/10000.0;
+    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;
diff --git a/src/likelihood.h b/src/likelihood.h
index d53b777..b4b3be9 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -18,13 +18,8 @@
  * value for n. */
 #define MIN_RATIO -10
 
-/* Roughly equal to the fraction of light which is reflected at each PMT.
- *
- * Note: The proper thing to do here would be to separately integrate over the
- * track to calculate the amount of reflected light using the data from the
- * PMTR bank. However, this would make the likelihood function almost twice as
- * long, so we just assume a constant fraction of the light is reflected now. */
-#define CHARGE_FRACTION 0.3
+/* The fraction of reflected light which is detected. */
+#define CHARGE_FRACTION 1.0
 
 /* Dark rate of the PMTs (Hz).
  *
diff --git a/src/muon.c b/src/muon.c
index 575158b..1cccad6 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -256,6 +256,57 @@ double get_dEdx(double T, double rho)
     return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
 }
 
+double get_expected_reflected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
+{
+    double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
+
+    z = 1.0;
+
+    SUB(pmt_dir,pmt_pos,pos);
+
+    normalize(pmt_dir);
+
+    cos_theta_pmt = DOT(pmt_dir,pmt_normal);
+
+    if (cos_theta_pmt > 0) return 0;
+
+    /* Calculate the cosine of the angle between the track direction and the
+     * vector to the PMT. */
+    cos_theta = DOT(dir,pmt_dir);
+
+    /* Calculate total energy */
+    E = T + MUON_MASS;
+    p = sqrt(E*E - MUON_MASS*MUON_MASS);
+    beta = p/E;
+
+    omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r);
+
+    R = NORM(pos);
+
+    /* FIXME: I just calculate delta assuming 400 nm light. */
+    wavelength0 = 400.0;
+
+    if (R <= AV_RADIUS) {
+        n = get_index_snoman_d2o(wavelength0);
+    } else {
+        n = get_index_snoman_h2o(wavelength0);
+    }
+
+    if (beta < 1/n) return 0;
+
+    f = get_weighted_pmt_reflectivity(acos(-cos_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();
+
+    get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
+
+    l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
+
+    return f*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);
+}
+
 double get_expected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
 {
     double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
diff --git a/src/muon.h b/src/muon.h
index 253a9df..ce3db29 100644
--- a/src/muon.h
+++ b/src/muon.h
@@ -16,6 +16,7 @@ double muon_get_energy(double x, muon_energy *m);
 void muon_free_energy(muon_energy *m);
 double get_range(double T, double rho);
 double get_dEdx(double T, double rho);
+double get_expected_reflected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r);
 double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r);
 
 #endif