update find_peaks algorithm

Previously, the algorithm used to find peaks was to search for all peaks in the
Hough transform above some constant fraction of the highest peak. This
algorithm could have issues finding smaller peaks away from the highest peak.
The new algorithm instead finds the highest peak in the Hough transform and
then recomputes the Hough transform ignoring all PMT hits within the Cerenkov
cone of the first peak. The next peak is found from this transform and the
process is iteratively repeated until a certain number of peaks are found.

One disadvantage of this new system is that it will *always* find the same
number of peaks and this will usually be greater than the actual number of
rings in the event. This is not a problem though since when fitting the event
we loop over all possible peaks and do a quick fit to determine the starting
point and so false positives are OK because the real peaks will fit better
during this quick fit.

Another potential issue with this new method is that by rejecting all PMT hits
within the Cerenkov cone of the first peak we could miss a second peak very
close to the first peak. This is partially mitigated by the fact that when we
loop over all possible combinations of the particle ids and directions we allow
each peak to be used more than once. For example, when fitting for the
hypothesis that an event is caused by two electrons and one muon and given two
possible directions 1 and 2, we will fit for the following possible direction
combinations:

    1 1 1
    1 1 2
    1 2 1
    1 2 2
    2 2 1
    2 2 2

Therefore if there is a second ring close to the first it is possible to fit it
correctly since we will seed the quick fit with two particles pointing in the
same direction.

This commit also adds a few tests for new functions and changes the energy step
size during the quick fit to 10% of the starting energy value.

1 files changed, 79 insertions, 28 deletions

diff --git a/src/find_peaks.c b/src/find_peaks.c
index 05d8c15..e9030a0 100644
--- a/src/find_peaks.c
+++ b/src/find_peaks.c
@@ -6,6 +6,7 @@
 #include <stdlib.h> /* for malloc() */
 #include "optics.h"
 #include <math.h> /* for exp() */
+#include "misc.h"
 
 typedef struct peak {
     size_t i;
@@ -115,14 +116,37 @@ end:
     return;
 }
 
-void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, size_t m, double *result)
+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`. */
     size_t i, j, k;
-    double dir[3], pmt_dir[3], cos_theta;
+    double dir[3], pmt_dir[3], cos_theta2;
+    int skip;
+    double *sin_theta, *cos_theta, *sin_phi, *cos_phi;
+    double wavelength0, n_d2o;
 
-    double wavelength0 = 400.0;
-    double n_d2o = get_index_snoman_d2o(wavelength0);
+    sin_theta = malloc(sizeof(double)*n);
+    cos_theta = malloc(sizeof(double)*n);
+    sin_phi = malloc(sizeof(double)*n);
+    cos_phi = malloc(sizeof(double)*n);
+
+    wavelength0 = 400.0;
+    n_d2o = get_index_snoman_d2o(wavelength0);
+
+    /* Precompute sin(theta), cos(theta), sin(phi), and cos(phi) to speed
+     * things up. */
+    for (i = 0; i < n; i++) {
+        sin_theta[i] = sin(x[i]);
+        cos_theta[i] = cos(x[i]);
+    }
+
+    for (i = 0; i < m; i++) {
+        sin_phi[i] = sin(y[i]);
+        cos_phi[i] = cos(y[i]);
+    }
+
+    /* Zero out the result array. */
+    for (i = 0; i < n*m; i++) result[i] = 0.0;
 
     for (i = 0; i < MAX_PMTS; i++) {
         if (pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
@@ -131,55 +155,82 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
 
         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) {
+                skip = 1;
+                break;
+            }
+        }
+
+        if (skip) continue;
+
         for (j = 0; j < n; j++) {
             for (k = 0; k < m; k++) {
-                double r = 1+x[j]*x[j]+y[k]*y[k];
-                dir[0] = 2*x[j]/r;
-                dir[1] = 2*y[k]/r;
-                dir[2] = (-1+x[j]*x[j]+y[k]*y[k])/r;
+                dir[0] = sin_theta[j]*cos_phi[k];
+                dir[1] = sin_theta[j]*sin_phi[k];
+                dir[2] = cos_theta[j];
 
-                cos_theta = DOT(pmt_dir,dir);
+                cos_theta2 = DOT(pmt_dir,dir);
 
-                result[j*m + k] += ev->pmt_hits[i].qhs*exp(-fabs(cos_theta-1.0/n_d2o)/0.1);
+                result[j*n + k] += ev->pmt_hits[i].qhs*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1);
             }
         }
     }
+
+    free(sin_theta);
+    free(cos_theta);
+    free(sin_phi);
+    free(cos_phi);
 }
 
-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 threshold)
+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)
 {
+    /* Finds rings in the event `ev` by looking for peaks in the Hough
+     * transform. The directions of potential particle rings are stored in the
+     * arrays `peak_theta` and `peak_phi`. The number of rings found are stored
+     * in the variable `npeaks`.
+     *
+     * This function first computes the Hough transform of the event `ev` and
+     * looks for the highest peak. The next peak is then found by computing the
+     * Hough transform ignoring any PMTs within the Cerenkov cone of the first
+     * peak. This process is repeated until `max_peaks` peaks are found. */
     size_t i;
+    double *x, *y, *result, *dir;
+    size_t *imax, *jmax;
+    size_t max;
 
-    double *x = calloc(n,sizeof(double));
-    double *y = calloc(m,sizeof(double));
-    double *result = calloc(n*m,sizeof(double));
+    x = calloc(n,sizeof(double));
+    y = calloc(m,sizeof(double));
+    result = malloc(n*m*sizeof(double));
+    dir = calloc(max_peaks*3,sizeof(double));
 
     for (i = 0; i < n; i++) {
-        x[i] = -10 + 20.0*i/(n-1);
+        x[i] = i*M_PI/(n-1);
     }
 
     for (i = 0; i < m; i++) {
-        y[i] = -10 + 20.0*i/(m-1);
+        y[i] = i*2*M_PI/(m-1);
     }
 
-    get_hough_transform(ev,pos,x,y,n,m,result);
+    imax = calloc(max_peaks,sizeof(size_t));
+    jmax = calloc(max_peaks,sizeof(size_t));
 
-    size_t *imax = calloc(max_peaks,sizeof(size_t));
-    size_t *jmax = calloc(max_peaks,sizeof(size_t));
-
-    find_peaks_array(result,n,m,imax,jmax,npeaks,max_peaks,threshold);
-
-    for (i = 0; i < *npeaks; i++) {
-        double R, theta;
-        R = sqrt(x[imax[i]]*x[imax[i]]+y[jmax[i]]*y[jmax[i]]);
-        theta = atan2(y[jmax[i]],x[imax[i]]);
-        peak_theta[i] = 2*atan(1/R);
-        peak_phi[i] = theta;
+    for (i = 0; i < max_peaks; i++) {
+        get_hough_transform(ev,pos,x,y,n,m,result,dir,i);
+        max = argmax(result,n*m);
+        peak_theta[i] = x[max/m];
+        peak_phi[i] = y[max % m];
+        get_dir(dir+i*3,peak_theta[i],peak_phi[i]);
     }
 
+    *npeaks = max_peaks;
+
     free(imax);
     free(jmax);
     free(x);
     free(y);
     free(result);
+    free(dir);
 }