add a couple of improvements to the quad fitter and fix a bug in get_hough_transform()

This commit adds two improvements to the quad fitter:

1. I updated quad to weight the random PMT hit selection by the probability
that the PMT hit is a multiphoton hit. The idea here is that we really only
want to sample direct light and for high energy events the reflected and
scattered light is usually single photon.

2. I added an option to quad to only use points in the quad cloud which are
below a given quantile of t0. The idea here is that for particles like muons
which travel more than a few centimeters in the detector the quad cloud usually
looks like the whole track. Since we want the QUAD fitter to find the position
of the *start* of the track we select only those quad cloud points with an
early time so the position is closer to the position of the start of the track.

Also, I fixed a major bug in get_hough_transform() in which I was using the
wrong index variable when checking if a PMT was not flagged, a normal PMT, and
was hit. This was causing the algorithm to completely miss finding more than
one ring while I was testing it.

1 files changed, 3 insertions, 2 deletions

diff --git a/src/find_peaks.c b/src/find_peaks.c
index 7d350d5..6109896 100644
--- a/src/find_peaks.c
+++ b/src/find_peaks.c
@@ -197,8 +197,9 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
     /* 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 (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
+            if (ev->pmt_hits[j].flags || pmts[j].pmt_type != PMT_NORMAL || !ev->pmt_hits[j].hit) continue;
 
             SUB(pmt_dir,pmts[j].pos,pos);
 
@@ -304,7 +305,7 @@ void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta,
         /* 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++) {
+        for (j = 0; j < *npeaks; 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;