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, 13 insertions, 7 deletions

diff --git a/src/fit.c b/src/fit.c
index 9cf1809..5f3db19 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5298,11 +5298,6 @@ int bisect_energy(double qhs_sum, double distance, int id, double cos_theta, dou
     pars.qhs_sum = qhs_sum;
     pars.cos_theta = cos_theta;
 
-    s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
-
-    F.function = &get_expected_photons_root;
-    F.params = &pars;
-
     switch (id) {
     case IDP_E_MINUS:
     case IDP_E_PLUS:
@@ -5326,6 +5321,17 @@ int bisect_energy(double qhs_sum, double distance, int id, double cos_theta, dou
     /* Calculate the Cerenkov threshold for a muon. */
     Tmin = mass/sqrt(1.0-1.0/(n_d2o*n_d2o)) - mass;
 
+    if (get_expected_photons(id,Tmin,distance,cos_theta) > qhs_sum)
+        return Tmin;
+
+    if (get_expected_photons(id,Tmax,distance,cos_theta) < qhs_sum)
+        return Tmax;
+
+    s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
+
+    F.function = &get_expected_photons_root;
+    F.params = &pars;
+
     gsl_root_fsolver_set(s, &F, Tmin, Tmax);
 
     do {
@@ -5432,7 +5438,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
     nlopt_set_min_objective(opt,nlopt_nll2,&fpars);
 
     /* Guess the position and t0 of the event using the QUAD fitter. */
-    status = quad(ev,pos,&t0,10000);
+    status = quad(ev,pos,&t0,10000,0.1);
 
     if (status || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) {
         /* If the QUAD fitter fails or returns something outside the PSUP or
@@ -5474,7 +5480,7 @@ 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),0.1);
+    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;