1 files changed, 40 insertions, 16 deletions

diff --git a/src/find_peaks.c b/src/find_peaks.c
index 04bf40a..135e4e3 100644
--- a/src/find_peaks.c
+++ b/src/find_peaks.c
@@ -252,21 +252,29 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
     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)
+/* 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 after subtracting off the previous ring and ignoring any PMTs
+ * within the Cerenkov cone of the first peak. This process is repeated until
+ * `max_peaks` peaks are found.
+ *
+ * In addition, only peaks which are `delta` radians away from other peaks are
+ * 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)
 {
-    /* 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;
+    size_t i, j;
     double *x, *y, *result, *dir;
     size_t *imax, *jmax;
     size_t max;
+    double theta, phi;
+    int unique;
 
     x = calloc(n,sizeof(double));
     y = calloc(m,sizeof(double));
@@ -284,15 +292,31 @@ void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta,
     imax = calloc(max_peaks,sizeof(size_t));
     jmax = calloc(max_peaks,sizeof(size_t));
 
+    *npeaks = 0;
+
     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]);
-    }
+        theta = x[max/m];
+        phi = y[max % m];
+        get_dir(dir+i*3,theta,phi);
+
+        /* Only add the latest peak to the results if it's more than `delta`
+         * radians away from all other results so far. */
+        unique = 1;
+        for (j = 0; j < i; j++) {
+            /* If the angle between the latest peak and a previous peak is
+             * within `delta` radians, it's not unique. */
+            if (acos(DOT(dir+i*3,dir+j*3)) < delta) unique = 0;
+        }
 
-    *npeaks = max_peaks;
+        /* Add it to the results if it's unique. */
+        if (unique) {
+            peak_theta[*npeaks] = theta;
+            peak_phi[*npeaks] = phi;
+            *npeaks += 1;
+        }
+    }
 
     free(imax);
     free(jmax);