initial commit of likelihood fit for muons

This commit contains code to fit for the energy, position, and direction of
muons in the SNO detector. Currently, we read events from SNOMAN zebra files
and fill an event struct containing the PMT hits and fit it with the Nelder
Mead simplex algorithm from GSL.

I've also added code to read in ZEBRA title bank files to read in the DQXX
files for a specific run. Any problems with channels in the DQCH and DQCR banks
are flagged in the event struct by masking in a bit in the flags variable and
these PMT hits are not included in the likelihood calculation.

The likelihood for an event is calculated by integrating along the particle
track for each PMT and computing the expected number of PE. The charge
likelihood is then calculated by looping over all possible number of PE and
computing:

    P(q|n)*P(n|mu)

where q is the calibrated QHS charge, n is the number of PE, and mu is the
expected number of photoelectrons. The latter is calculated assuming the
distribution of PE at a given PMT follows a Poisson distribution (which I think
should be correct given the track, but is probably not perfect for tracks which
scatter a lot).

The time part of the likelihood is calculated by integrating over the track for
each PMT and calculating the average time at which the PMT is hit. We then
assume the PDF for the photons to arrive is approximately a delta function and
compute the first order statistic for a given time to compute the probability
that the first photon arrived at a given time. So far I've only tested this
with single tracks but the method was designed to be easy to use when you are
fitting for multiple particles.

1 files changed, 67 insertions, 16 deletions

diff --git a/muon.c b/muon.c
index 12a68b2..4b46769 100644
--- a/muon.c
+++ b/muon.c
@@ -5,6 +5,14 @@
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_spline.h>
 #include <math.h>
+#include "optics.h"
+#include "quantum_efficiency.h"
+#include "solid_angle.h"
+#include "pdg.h"
+#include "vector.h"
+#include "muon.h"
+#include "sno.h"
+#include "scattering.h"
 
 static int initialized = 0;
 
@@ -95,8 +103,6 @@ static int init()
         j = 0;
         while (str) {
             value = strtod(str, NULL);
-            /* According to the file, the values are stored for wavelengths
-             * between 230 and 700 in 1 nm increments. */
             switch (j) {
             case 0:
                 x[n] = value;
@@ -156,17 +162,17 @@ double get_range(double T, double rho)
     return gsl_spline_eval(spline_range, T, acc_range)/rho;
 }
 
-double get_E(double T, double x, double rho)
+double get_T(double T0, double x, double rho)
 {
-    /* Returns the approximate energy of a muon in water after travelling `x`
-     * cm with an initial kinetic energy `T`.
+    /* Returns the approximate kinetic energy of a muon in water after
+     * travelling `x` cm with an initial kinetic energy `T`.
      *
      * `T` should be in MeV, `x` in cm, and `rho` in g/cm^3.
      *
      * Return value is in MeV.
      *
      * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
-    double a, b, range, E;
+    double a, b, range, T;
 
     if (!initialized) {
         if (init()) {
@@ -174,19 +180,21 @@ double get_E(double T, double x, double rho)
         }
     }
 
-    range = gsl_spline_eval(spline_range, T, acc_range)/rho;
-    /* FIXME: Need to double check if it's kosher to be using kinetic energies
-     * here instead of the total energy. Equation 1 in the document uses the
-     * total energy, but here I'm using the critical energy in kinetic energy,
-     * so I should check to see if I need to convert both. */
-    b = log(1 + T/MUON_CRITICAL_ENERGY)/range;
-    a = MUON_CRITICAL_ENERGY*b;
+    range = get_range(T0, rho);
 
-    E = T + a*(1-exp(b*x))/b;
+    /* This comes from Equation 33.42 in the PDG Passage of Particles Through
+     * Matter article. */
+    b = log(1 + T0/MUON_CRITICAL_ENERGY)/range;
+    /* Now we compute the ionization energy loss from the known range and b. */
+    a = b*T0/(exp(b*range)-1.0);
 
-    if (E < 0) return 0;
+    /* Compute the kinetic energy after travelling a distance `x` in the
+     * continuous slowing down approximation. */
+    T = -a/b + (T0+a/b)*exp(-b*x);
 
-    return E;
+    if (T < 0) return 0;
+
+    return T;
 }
 
 double get_dEdx(double T, double rho)
@@ -207,3 +215,46 @@ double get_dEdx(double T, double rho)
 
     return gsl_spline_eval(spline_dEdx, T, acc_dEdx)/rho;
 }
+
+double get_expected_charge(double x, double T, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
+{
+    double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R;
+
+    z = 1.0;
+
+    SUB(pmt_dir,pmt_pos,pos);
+    normalize(pmt_dir);
+
+    if (DOT(pmt_dir,pmt_normal) > 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);
+
+    if (R <= AV_RADIUS) {
+        rho = HEAVY_WATER_DENSITY;
+    } else {
+        rho = WATER_DENSITY;
+    }
+
+    /* FIXME: I just calculate delta assuming 400 nm light. */
+    wavelength0 = 400.0;
+    n = get_index(rho, wavelength0, 10.0);
+
+    if (beta < 1/n) return 0;
+
+    /* FIXME: is this formula valid for muons? */
+    theta0 = get_scattering_rms(x,p,beta,z,rho);
+
+    /* FIXME: add angular response and scattering/absorption. */
+    return 2*omega*2*M_PI*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0)/(sqrt(2*M_PI)*theta0);
+}