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, 134 insertions, 28 deletions

diff --git a/test.c b/test.c
index 7c5dd33..747e664 100644
--- a/test.c
+++ b/test.c
@@ -3,6 +3,11 @@
 #include <stdio.h>
 #include "optics.h"
 #include "muon.h"
+#include "mt19937ar.h"
+#include <gsl/gsl_sf_gamma.h>
+#include "misc.h"
+#include "sno_charge.h"
+#include <gsl/gsl_integration.h>
 
 typedef int testFunction(char *err);
 
@@ -19,30 +24,30 @@ struct refractive_index_results {
     double wavelength;
     double n;
 } refractive_index_results[] = {
-    {1000.0, 0 ,  226.50, 1.39468},
-    {1000.0, 10,  226.50, 1.39439},
-    {1000.0, 20,  226.50, 1.39353},
-    {1000.0, 30,  226.50, 1.39224},
-    {1000.0, 0 ,  404.41, 1.34431},
-    {1000.0, 10,  404.41, 1.34404},
-    {1000.0, 20,  404.41, 1.34329},
-    {1000.0, 30,  404.41, 1.34218},
-    {1000.0, 0 ,  589.00, 1.33447},
-    {1000.0, 10,  589.00, 1.33422},
-    {1000.0, 20,  589.00, 1.33350},
-    {1000.0, 30,  589.00, 1.33243},
-    {1000.0, 0 ,  632.80, 1.33321},
-    {1000.0, 10,  632.80, 1.33296},
-    {1000.0, 20,  632.80, 1.33224},
-    {1000.0, 30,  632.80, 1.33118},
-    {1000.0, 0 , 1013.98, 1.32626},
-    {1000.0, 10, 1013.98, 1.32604},
-    {1000.0, 20, 1013.98, 1.32537},
-    {1000.0, 30, 1013.98, 1.32437},
-    {1000.0, 0 , 2325.42, 1.27663},
-    {1000.0, 10, 2325.42, 1.27663},
-    {1000.0, 20, 2325.42, 1.27627},
-    {1000.0, 30, 2325.42, 1.27563},
+    {1.0, 0 ,  226.50, 1.39468},
+    {1.0, 10,  226.50, 1.39439},
+    {1.0, 20,  226.50, 1.39353},
+    {1.0, 30,  226.50, 1.39224},
+    {1.0, 0 ,  404.41, 1.34431},
+    {1.0, 10,  404.41, 1.34404},
+    {1.0, 20,  404.41, 1.34329},
+    {1.0, 30,  404.41, 1.34218},
+    {1.0, 0 ,  589.00, 1.33447},
+    {1.0, 10,  589.00, 1.33422},
+    {1.0, 20,  589.00, 1.33350},
+    {1.0, 30,  589.00, 1.33243},
+    {1.0, 0 ,  632.80, 1.33321},
+    {1.0, 10,  632.80, 1.33296},
+    {1.0, 20,  632.80, 1.33224},
+    {1.0, 30,  632.80, 1.33118},
+    {1.0, 0 , 1013.98, 1.32626},
+    {1.0, 10, 1013.98, 1.32604},
+    {1.0, 20, 1013.98, 1.32537},
+    {1.0, 30, 1013.98, 1.32437},
+    {1.0, 0 , 2325.42, 1.27663},
+    {1.0, 10, 2325.42, 1.27663},
+    {1.0, 20, 2325.42, 1.27627},
+    {1.0, 30, 2325.42, 1.27563},
 };
 
 /* Table of the values of solid angle for various values of r0/r and L/r.
@@ -109,7 +114,7 @@ int isclose(double a, double b, double rel_tol, double abs_tol)
     return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol);
 }
 
-int test_muon_get_E(char *err)
+int test_muon_get_T(char *err)
 {
     /* A very simple test to make sure the energy as a function of distance
      * along the track makes sense. Should include more detailed tests later. */
@@ -117,7 +122,7 @@ int test_muon_get_E(char *err)
 
     /* Assume initial kinetic energy is 1 GeV. */
     T = 1000.0;
-    E = get_E(T,1e-9,1.0);
+    E = get_T(T,1e-9,1.0);
 
     /* At the beginning of the track we should have roughly the same energy. */
     if (!isclose(E, T, 1e-5, 0)) {
@@ -128,7 +133,7 @@ int test_muon_get_E(char *err)
     /* Assume initial kinetic energy is 1 GeV. */
     T = 1000.0;
     range = get_range(T,1.0);
-    E = get_E(T,range,1.0);
+    E = get_T(T,range,1.0);
 
     /* At the end of the track we should have roughly the same energy. */
     if (!isclose(E, 0, 1e-5, 1e-5)) {
@@ -189,6 +194,38 @@ int test_refractive_index(char *err)
     return 0;
 }
 
+int test_solid_angle_approx(char *err)
+{
+    /* Tests the get_solid_angle_approx() function. */
+    int i;
+    double pmt[3] = {0,0,0};
+    double pos[3] = {0,0,1};
+    double n[3] = {0,0,1};
+    double r = 1.0;
+    double solid_angle;
+
+    solid_angle = get_solid_angle_approx(pos,pmt,n,r);
+
+    if (!isclose(solid_angle, 2*M_PI*(1-1/sqrt(2)), 1e-2, 0)) {
+        sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, 2*M_PI*(1-1/sqrt(2)));
+        return 1;
+    }
+
+    for (i = 0; i < sizeof(solid_angle_results)/sizeof(struct solid_angle_results); i++) {
+        pos[0] = solid_angle_results[i].r0*r;
+        pos[2] = solid_angle_results[i].L*r;
+
+        solid_angle = get_solid_angle_approx(pos,pmt,n,r);
+
+        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-2, 0)) {
+            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, solid_angle_results[i].omega);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
 int test_solid_angle(char *err)
 {
     /* Tests the get_solid_angle() function. */
@@ -221,15 +258,84 @@ int test_solid_angle(char *err)
     return 0;
 }
 
+static double sno_charge(double q, void *params)
+{
+    int n = ((int *) params)[0];
+    return pq(q,n);
+}
+
+int test_sno_charge_integral(char *err)
+{
+    /* Test that the single PE charge distribution integrates to one. */
+    double result, error;
+    gsl_function F;
+    size_t nevals;
+    gsl_integration_cquad_workspace *w;
+    int params[1];
+
+    w = gsl_integration_cquad_workspace_alloc(100);
+
+    F.function = &sno_charge;
+    params[0] = 1;
+    F.params = params;
+
+    init_charge();
+
+    gsl_integration_cquad(&F, -10.0, 1000.0, 0, 1e-9, w, &result, &error, &nevals);
+
+    if (!isclose(result, 1.0, 1e-9, 1e-9)) {
+        sprintf(err, "integral of single PE charge distribution is %.2f", result);
+        goto err;
+    }
+
+    gsl_integration_cquad_workspace_free(w);
+
+    return 0;
+
+err:
+    gsl_integration_cquad_workspace_free(w);
+    return 1;
+}
+
+int test_logsumexp(char *err)
+{
+    /* Tests the logsumexp() function. */
+    int i, j;
+    double logp[100], mu, result, expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        mu = genrand_real2();
+
+        for (j = 0; j < sizeof(logp)/sizeof(double); j++) {
+            logp[j] = -mu + j*log(mu) - gsl_sf_lnfact(j);
+        }
+        result = logsumexp(logp, sizeof(logp)/sizeof(double));
+
+        expected = 0.0;
+
+        if (!isclose(result, expected, 1e-9, 1e-9)) {
+            sprintf(err, "logsumexp(logp) = %.5g, but expected %.5g", result, expected);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
 struct tests {
     testFunction *test;
     char *name;
 } tests[] = {
     {test_solid_angle, "test_solid_angle"},
+    {test_solid_angle_approx, "test_solid_angle_approx"},
     {test_refractive_index, "test_refractive_index"},
     {test_muon_get_dEdx, "test_muon_get_dEdx"},
     {test_muon_get_range, "test_muon_get_range"},
-    {test_muon_get_E, "test_muon_get_E"}
+    {test_muon_get_T, "test_muon_get_T"},
+    {test_logsumexp, "test_logsumexp"},
+    {test_sno_charge_integral, "test_sno_charge_integral"}
 };
 
 int main(int argc, char **argv)