update find_peaks algorithm

This commit updates the find peaks algorithm with several improvements which
together drastically improve its ability to find Cerenkov rings:

- when computing the Hough transform, instead of charge we weight each PMT hit
  by the probability that it is a multi-photon PMT hit
- we don't subtract off previously found rings (this makes the code simpler and
  I don't think it previously had a huge effect)
- ignore PMT hits who are within approximately 5 degrees of any previously
  found ring (previously we ignored all hits within the center of previously
  found rings)
- ignore PMT hits which have a time residual of more than 10 nanoseconds to
  hopefully ignore more reflected and/or scattered light
- switch from weighting the Hough transform by exp(-fabs(cos(theta)-1/n)/0.1)
  -> exp(-pow(cos(theta)-1/n,2)/0.01). I'm still not sure if this has a huge
  effect, but the reason I switched is that the PDF for Cerenkov light looks
  closer to the second form.
- switch to calling quad with f = 1.0 in test-find-peaks (I still need to add
  this update to fit.c but will do that in a later commit).

1 files changed, 56 insertions, 50 deletions

diff --git a/src/find_peaks.c b/src/find_peaks.c
index 6109896..14c33e5 100644
--- a/src/find_peaks.c
+++ b/src/find_peaks.c
@@ -24,6 +24,9 @@
 #include <math.h> /* for exp() */
 #include "misc.h"
 #include "pmt.h"
+#include "sno_charge.h"
+#include <gsl/gsl_randist.h>
+#include "pdg.h"
 
 typedef struct peak {
     size_t i;
@@ -133,33 +136,7 @@ end:
     return;
 }
 
-double get_qhs_avg(event *ev, double *pos, double *dir)
-{
-    /* Returns the average QHS value for all PMTs within the Cerenkov cone of a
-     * particle at position `pos` travelling in direction `dir`. */
-    size_t i;
-    double qhs_sum, n_d2o;
-    double pmt_dir[3];
-    size_t n;
-
-    n_d2o = get_index_snoman_d2o(400.0);
-
-    n = 0;
-    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;
-        SUB(pmt_dir,pmts[i].pos,pos);
-        normalize(pmt_dir);
-        if (DOT(pmt_dir,dir) > 1/n_d2o) {
-            qhs_sum += ev->pmt_hits[i].q;
-            n += 1;
-        }
-    }
-
-    return qhs_sum/n;
-}
-
-void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, size_t m, double *result, double *last, size_t len)
+void get_hough_transform(event *ev, double *pos, double t0, double *x, double *y, size_t n, size_t m, double *result, double *last, size_t len)
 {
     /* Computes the "Hough transform" of the event `ev` and stores it in `result`. */
     size_t i, j, k;
@@ -167,7 +144,11 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
     int skip;
     double *sin_theta, *cos_theta, *sin_phi, *cos_phi;
     double wavelength0, n_d2o;
-    double qhs[MAX_PMTS], qhs_avg;
+    static double w[MAX_PMTS];
+    double expected_pe;
+    int nhit;
+    double pq;
+    double dt;
 
     sin_theta = malloc(sizeof(double)*n);
     cos_theta = malloc(sizeof(double)*n);
@@ -189,27 +170,46 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
         cos_phi[i] = cos(y[i]);
     }
 
+    expected_pe = 0.0;
+    nhit = 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[i] = ev->pmt_hits[i].q;
-    }
-
-    /* Subtract off previous rings. */
-    for (i = 0; i < len; i++) {
-        qhs_avg = get_qhs_avg(ev,pos,last+3*i);
-
-        for (j = 0; j < MAX_PMTS; j++) {
-            if (ev->pmt_hits[j].flags || pmts[j].pmt_type != PMT_NORMAL || !ev->pmt_hits[j].hit) continue;
-
-            SUB(pmt_dir,pmts[j].pos,pos);
-
-            normalize(pmt_dir);
-
-            cos_theta2 = DOT(pmt_dir,last+3*i);
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
 
-            qhs[j] -= qhs_avg*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1);
+        if (ev->pmt_hits[i].hit) {
+            expected_pe += ev->pmt_hits[i].q;
+            nhit += 1;
+        }
+    }
+    expected_pe /= nhit;
 
-            if (qhs[j] < 0.0) qhs[j] = 0.0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+
+        if (ev->pmt_hits[i].hit) {
+            /* Weights are equal to the probability that the hit is > 1 photon.
+             * This is supposed to be a rough proxy for the probability that
+             * it's not reflected and/or scattered light (which is typically
+             * single photons). */
+            if (ev->pmt_hits[i].q > FIND_PEAKS_MAX_PE*get_qmean()/2) {
+                /* If the charge is greater than FIND_PEAKS_MAX_PE, it's almost
+                 * certainly multiple photons so we don't waste time
+                 * calculating P(1 PE|q). */
+                w[i] = 1.0;
+            } else {
+                /* We want to calculate P(multiple photons|q) which we calculate as:
+                 *
+                 *     P(multiple photons|q) = 1 - P(1 PE|q) = 1 - P(q|1 PE)P(1 PE)/P(q)
+                 *
+                 * To calculate the second two quantities we assume the number
+                 * of PE is Poisson distributed with a mean equal to
+                 * expected_photons. */
+                pq = 0.0;
+                for (j = 1; j < FIND_PEAKS_MAX_PE; j++) {
+                    pq += get_pq(ev->pmt_hits[i].q,j)*gsl_ran_poisson_pdf(j,expected_pe);
+                }
+
+                w[i] = 1-get_pq(ev->pmt_hits[i].q,1)*gsl_ran_poisson_pdf(1,expected_pe)/pq;
+            }
         }
     }
 
@@ -221,12 +221,18 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
 
         SUB(pmt_dir,pmts[i].pos,pos);
 
+        dt = ev->pmt_hits[i].t - t0 - NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT;
+
+        /* Skip PMT hits with a time residual of more than 10 nanoseconds since
+         * they are most likely to be scattered and/or reflected light. */
+        if (dt > 10) continue;
+
         normalize(pmt_dir);
 
         /* Ignore PMT hits within the Cerenkov cone of previously found peaks. */
         skip = 0;
         for (j = 0; j < len; j++) {
-            if (DOT(pmt_dir,last+3*j) > 1/n_d2o) {
+            if (fabs(DOT(pmt_dir,last+3*j) - 1/n_d2o) < 0.1) {
                 skip = 1;
                 break;
             }
@@ -242,7 +248,7 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
 
                 cos_theta2 = DOT(pmt_dir,dir);
 
-                result[j*n + k] += qhs[i]*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1);
+                result[j*n + k] += w[i]*exp(-pow(cos_theta2-1.0/n_d2o,2)/0.01);
             }
         }
     }
@@ -268,7 +274,7 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
  * returned. So, for example if three peaks are found but the first two are
  * very close to each other, this function will only return the first and the
  * third peaks. */
-void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, double *peak_phi, size_t *npeaks, size_t max_peaks, double delta)
+void find_peaks(event *ev, double *pos, double t0, size_t n, size_t m, double *peak_theta, double *peak_phi, size_t *npeaks, size_t max_peaks, double delta)
 {
     size_t i, j;
     double *x, *y, *result, *dir;
@@ -296,7 +302,7 @@ void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta,
     *npeaks = 0;
 
     for (i = 0; i < max_peaks; i++) {
-        get_hough_transform(ev,pos,x,y,n,m,result,dir,i);
+        get_hough_transform(ev,pos,t0,x,y,n,m,result,dir,i);
         max = argmax(result,n*m);
         theta = x[max/m];
         phi = y[max % m];