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