update find_peaks algorithm to subtract previous rings

Previously the find peaks algorithm would ignore any PMT hits within the
Cerenkov ring of previously found rings. This had the problem that occasionally
the algorithm would repeatedly find the same direction due to hits outside of
the Cerenkov cone. The new update was inspired by how SuperK does this and
instead we "subtract" off the previous rings by subtracting the average qhs
times e^(-cos(theta-1/n)/0.1) from each PMT for each previous ring.

Based on some quick testing this seems a lot better than the previous
algorithm, but still probably needs some tweaking.

1 files changed, 52 insertions, 1 deletions

diff --git a/src/find_peaks.c b/src/find_peaks.c
index e9030a0..008e549 100644
--- a/src/find_peaks.c
+++ b/src/find_peaks.c
@@ -7,6 +7,7 @@
 #include "optics.h"
 #include <math.h> /* for exp() */
 #include "misc.h"
+#include "pmt.h"
 
 typedef struct peak {
     size_t i;
@@ -116,6 +117,32 @@ 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].qhs;
+            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)
 {
     /* Computes the "Hough transform" of the event `ev` and stores it in `result`. */
@@ -124,6 +151,7 @@ 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;
 
     sin_theta = malloc(sizeof(double)*n);
     cos_theta = malloc(sizeof(double)*n);
@@ -145,6 +173,29 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
         cos_phi[i] = cos(y[i]);
     }
 
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
+        qhs[i] = ev->pmt_hits[i].qhs;
+    }
+
+    /* 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 (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);
+
+            qhs[j] -= qhs_avg*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1);
+
+            if (qhs[j] < 0.0) qhs[j] = 0.0;
+        }
+    }
+
     /* Zero out the result array. */
     for (i = 0; i < n*m; i++) result[i] = 0.0;
 
@@ -174,7 +225,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] += ev->pmt_hits[i].qhs*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1);
+                result[j*n + k] += qhs[i]*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1);
             }
         }
     }