update fit to guess energy, direction, and t0

This commit updates the initial guess for the energy using a simple heuristic
of ~6 hits/MeV. I also updated the initial phase where we do a bunch of "quick"
minimizations to loop over a series of starting positions and automatically
calculate the approximate direction and t0 for the event.

2 files changed, 147 insertions, 28 deletions

diff --git a/src/fit.c b/src/fit.c
index 045141d..c26b27e 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -17,6 +17,7 @@
 #include <errno.h> /* for errno */
 #include "pdg.h"
 #include "optics.h"
+#include "vector.h"
 
 #define EV_RECORD    0x45562020 // 'EV  '  (as written to ZDAB file)
 
@@ -36,15 +37,34 @@ static struct startingParameters {
     double x;
     double y;
     double z;
-    double theta;
-    double phi;
 } startingParameters[] = {
-    {0.0, 0.0, 0.0,    0.0, 0.0},
-    {0.0, 0.0, 0.0, M_PI/2, 0.0},
-    {0.0, 0.0, 0.0, M_PI/2, M_PI/2},
-    {0.0, 0.0, 0.0, M_PI/2, M_PI},
-    {0.0, 0.0, 0.0, M_PI/2, 3*M_PI/2},
-    {0.0, 0.0, 0.0,   M_PI, 0.0},
+    {   0.0,   0.0,    0.0},
+    { 100.0,   0.0,    0.0},
+    { 100.0, 100.0,    0.0},
+    {   0.0, 100.0,    0.0},
+    {-100.0, 100.0,    0.0},
+    {-100.0,   0.0,    0.0},
+    {-100.0,-100.0,    0.0},
+    {   0.0,-100.0,    0.0},
+    { 100.0,-100.0,    0.0},
+    {   0.0,   0.0,  100.0},
+    { 100.0,   0.0,  100.0},
+    { 100.0, 100.0,  100.0},
+    {   0.0, 100.0,  100.0},
+    {-100.0, 100.0,  100.0},
+    {-100.0,   0.0,  100.0},
+    {-100.0,-100.0,  100.0},
+    {   0.0,-100.0,  100.0},
+    { 100.0,-100.0,  100.0},
+    {   0.0,   0.0, -100.0},
+    { 100.0,   0.0, -100.0},
+    { 100.0, 100.0, -100.0},
+    {   0.0, 100.0, -100.0},
+    {-100.0, 100.0, -100.0},
+    {-100.0,   0.0, -100.0},
+    {-100.0,-100.0, -100.0},
+    {   0.0,-100.0, -100.0},
+    { 100.0,-100.0, -100.0},
 };
 
 double nll(unsigned int n, const double *x, double *grad, void *params)
@@ -96,25 +116,127 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
     return fval;
 }
 
+double guess_t0(event *ev, double *pos)
+{
+    /* Returns a guess for the t0 of event `ev` given a position `pos`. The t0
+     * is calculated by computing the charge weighted average of the time
+     * difference between each PMT's hit time and the time it would take a
+     * photon to travel from `pos` to the PMT. */
+    size_t i;
+    double pmt_dir[3], distance, t0, n, qhs_sum;
+
+    /* Compute the index of refraction for heavy water at 400 nm. */
+    n = get_index_snoman_d2o(400.0);
+
+    t0 = 0.0;
+    qhs_sum = 0.0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        /* Only look at normal PMTs that were hit and aren't flagged. */
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
+
+        /* Compute the vector between `pos` and the PMT. */
+        SUB(pmt_dir,pmts[i].pos,pos);
+
+        /* Compute the distance to the PMT. */
+        distance = NORM(pmt_dir);
+
+        /* Add the charge weighted time difference between the PMT hit time and
+         * the time it would take a photon to hit the PMT from `pos`. */
+        t0 += ev->pmt_hits[i].qhs*(ev->pmt_hits[i].t - distance*n/SPEED_OF_LIGHT);
+
+        qhs_sum += ev->pmt_hits[i].qhs;
+    }
+
+    /* Divide by the total QHS sum. */
+    t0 /= qhs_sum;
+
+    return t0;
+}
+
+void guess_direction(event *ev, double *pos, double *theta, double *phi)
+{
+    /* Returns the approximate direction of the event from a given position.
+     * The direction is computed by taking the charge weighted average of the
+     * vectors from `pos` to each hit PMT.
+     *
+     * A possible improvement here is to do something like a Hough transform
+     * where we map PMT hits to cones with a 42 degree angle and then search
+     * for the highest point in the transformed space. This method might also
+     * generalize better when searching for a second ring since we can look for
+     * secondary peaks not near the highest peak. */
+    size_t i;
+    double pmt_dir[3], dir[3];
+
+    dir[0] = 0.0;
+    dir[1] = 0.0;
+    dir[2] = 0.0;
+
+    for (i = 0; i < MAX_PMTS; i++) {
+        /* Only look at normal PMTs that were hit and aren't flagged. */
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
+
+        /* Compute the vector between `pos` and the PMT. */
+        SUB(pmt_dir,pmts[i].pos,pos);
+
+        /* Normalize it. */
+        normalize(pmt_dir);
+
+        /* Multiply the vector by the QHS charge in the PMT. */
+        MUL(pmt_dir,ev->pmt_hits[i].qhs);
+
+        /* Add this to the estimated direction. */
+        ADD(dir,dir,pmt_dir);
+    }
+
+    /* Normalize the charge weighted sum. */
+    normalize(dir);
+
+    /* Compute the spherical coordinates for `dir`. */
+    *theta = acos(dir[2]);
+    *phi = atan2(dir[1],dir[0]);
+}
+
 int fit_event(event *ev, double *xopt, double *fmin)
 {
     size_t i;
     fitParams fpars;
-    double x[9], ss[9], lb[9], ub[9], fval, n;
+    double x[9], ss[9], lb[9], ub[9], fval, n, qhs_sum, x0[9], T0, Tmin;
     int rv;
 
     nlopt_opt opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
     nlopt_set_min_objective(opt,nll,&fpars);
 
-    x[0] = 1200.0;
-    x[1] = 0.0;
-    x[2] = 0.0;
-    x[3] = 0.0;
-    x[4] = 1.57;
-    x[5] = 0.0;
-    x[6] = 130.0;
-    x[7] = 0.0;
-    x[8] = 0.0;
+    /* Make a guess as to the energy. Right now we just use a simple
+     * approximation that the muon produces approximately 6 photons/MeV.
+     *
+     * FIXME: Should update this to something better. */
+
+    qhs_sum = 0.0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
+        qhs_sum += ev->pmt_hits[i].qhs;
+    }
+
+    T0 = qhs_sum/6.0;
+
+    n = get_index_snoman_d2o(400.0);
+
+    /* Calculate the Cerenkov threshold for a muon. */
+    Tmin = MUON_MASS/sqrt(1.0-1.0/(n*n)) - MUON_MASS;
+
+    /* If our guess is below the Cerenkov threshold, start at the Cerenkov
+     * threshold. */
+    if (T0 < Tmin) T0 = Tmin;
+
+    x0[0] = T0;
+    x0[1] = 0.0;
+    x0[2] = 0.0;
+    x0[3] = 0.0;
+    x0[4] = 1.57;
+    x0[5] = 0.0;
+    x0[6] = 130.0;
+    x0[7] = 0.0;
+    x0[8] = 0.0;
 
     ss[0] = 10.0;
     ss[1] = 10.0;
@@ -126,9 +248,7 @@ int fit_event(event *ev, double *xopt, double *fmin)
     ss[7] = 0.1;
     ss[8] = 0.1;
 
-    n = get_index_snoman_d2o(400.0);
-
-    lb[0] = MUON_MASS/sqrt(1.0-1.0/(n*n));
+    lb[0] = Tmin;
     lb[1] = -1000.0;
     lb[2] = -1000.0;
     lb[3] = -1000.0;
@@ -167,15 +287,14 @@ int fit_event(event *ev, double *xopt, double *fmin)
     nlopt_set_maxeval(opt, 100);
 
     for (i = 0; i < sizeof(startingParameters)/sizeof(startingParameters[0]); i++) {
-        x[0] = 1200.0;
+        memcpy(x,x0,sizeof(x));
         x[1] = startingParameters[i].x;
         x[2] = startingParameters[i].y;
         x[3] = startingParameters[i].z;
-        x[4] = startingParameters[i].theta;
-        x[5] = startingParameters[i].phi;
-        x[6] = 130.0;
-        x[7] = 0.0;
-        x[8] = 0.0;
+
+        guess_direction(ev,x+1,&x[4],&x[5]);
+
+        x[6] = guess_t0(ev,x+1);
 
         rv = nlopt_optimize(opt,x,&fval);
 
diff --git a/src/path.c b/src/path.c
index a15cba6..80bba50 100644
--- a/src/path.c
+++ b/src/path.c
@@ -8,7 +8,7 @@
 #include <stdlib.h> /* for malloc(), calloc(), etc. */
 #include "misc.h"
 
-#define N 100
+#define N 1000
 
 static double foo(double s, double range, double theta0, int k)
 {