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);
             }
         }
     }