update find_peaks() to only return unique peaks

This commit updates find_peaks() to only return peaks which are at least a
certain number of degrees apart from each other. This is because I found that
for many events the first few peaks would all be essentially the same direction
and so the fit was taking a lot of time fitting essentially the same seed
points. Since I now have to only try 3 peaks in order to get my grid jobs to
run for less than a few hours it's necessary to make sure we aren't just
fitting the same three directions for the "quick" minimization.

I also updated the fit to only use a maximum of 3 seed directions.

4 files changed, 47 insertions, 20 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);
diff --git a/src/find_peaks.h b/src/find_peaks.h
index 8511af6..0219561 100644
--- a/src/find_peaks.h
+++ b/src/find_peaks.h
@@ -22,6 +22,6 @@
 
 void find_peaks_array(double *x, size_t n, size_t m, size_t *imax, size_t *jmax, size_t *npeaks, size_t max_peaks, double threshold);
 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);
-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);
+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);
 
 #endif
diff --git a/src/fit.c b/src/fit.c
index b1882d8..ac385e3 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -47,7 +47,7 @@
 /* Maximum number of fit parameters. Should be at least 4 + 3*MAX_VERTICES. */
 #define MAX_PARS 100
 /* Maximum number of peaks to search for in Hough transform. */
-#define MAX_NPEAKS 5
+#define MAX_NPEAKS 10
 /* Maximum kinetic energy for any particle. */
 #define MAX_ENERGY 10000
 
@@ -5316,7 +5316,10 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
     ub[3] = GTVALID;
 
     /* Find the peaks in the Hough transform of the event. */
-    find_peaks(ev, pos, 100, 100, peak_theta, peak_phi, &npeaks, LEN(peak_theta));
+    find_peaks(ev, pos, 100, 100, peak_theta, peak_phi, &npeaks, LEN(peak_theta),0.1);
+
+    /* Don't fit more than 3 peaks for now. */
+    if (npeaks > 3) npeaks = 3;
 
     result = malloc(sizeof(size_t)*n*ipow(npeaks,n));
 
diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c
index a825f4b..dd3a9d8 100644
--- a/src/test-find-peaks.c
+++ b/src/test-find-peaks.c
@@ -363,7 +363,7 @@ int main(int argc, char **argv)
 
         get_hough_transform(&ev,pos,x,y,n,m,result,0,0);
 
-        find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,max_npeaks);
+        find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,max_npeaks,0.1);
 
         printf("gtid %i\n", ev.gtid);
         for (i = 0; i < npeaks; i++) {