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, 74 insertions, 8 deletions

diff --git a/src/test.c b/src/test.c
index d07ee11..a0b22ae 100644
--- a/src/test.c
+++ b/src/test.c
@@ -1421,7 +1421,6 @@ int test_quad(char *err)
     double fit_pos[3];
     double pmt_dir[3];
     double fit_t0;
-    double wavelength0;
     double n_d2o;
 
     init_genrand(0);
@@ -1438,8 +1437,7 @@ int test_quad(char *err)
         index[valid_pmts++] = i;
     }
 
-    wavelength0 = 400.0;
-    n_d2o = get_index_snoman_d2o(wavelength0);
+    n_d2o = get_avg_index_d2o();
 
     for (i = 0; i < 100; i++) {
         /* Generate a random position within a cube the size of the AV.
@@ -1449,6 +1447,13 @@ int test_quad(char *err)
         x0[0] = (genrand_real2()*2 - 1)*AV_RADIUS;
         x0[1] = (genrand_real2()*2 - 1)*AV_RADIUS;
         x0[2] = (genrand_real2()*2 - 1)*AV_RADIUS;
+
+        while (NORM(x0) >= PSUP_RADIUS) {
+            x0[0] = (genrand_real2()*2 - 1)*AV_RADIUS;
+            x0[1] = (genrand_real2()*2 - 1)*AV_RADIUS;
+            x0[2] = (genrand_real2()*2 - 1)*AV_RADIUS;
+        }
+
         t0 = genrand_real2()*GTVALID;
 
         /* Zero out all PMTs. */
@@ -1468,7 +1473,7 @@ int test_quad(char *err)
             ev.pmt_hits[hits[j]].t = t0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT;
         }
 
-        if (quad(&ev, fit_pos, &fit_t0, 10000)) {
+        if (quad(&ev, fit_pos, &fit_t0, 10000, 1.0)) {
             sprintf(err, "%s", quad_err);
             goto err;
         }
@@ -1526,7 +1531,6 @@ int test_quad_noise(char *err)
     double fit_pos[3];
     double pmt_dir[3];
     double fit_t0;
-    double wavelength0;
     double n_d2o;
 
     init_genrand(0);
@@ -1543,8 +1547,7 @@ int test_quad_noise(char *err)
         index[valid_pmts++] = i;
     }
 
-    wavelength0 = 400.0;
-    n_d2o = get_index_snoman_d2o(wavelength0);
+    n_d2o = get_avg_index_d2o();
 
     for (i = 0; i < 100; i++) {
         /* Generate a random position within a cube the size of the AV.
@@ -1556,6 +1559,12 @@ int test_quad_noise(char *err)
         x0[2] = (genrand_real2()*2 - 1)*AV_RADIUS;
         t0 = genrand_real2()*GTVALID;
 
+        while (NORM(x0) >= PSUP_RADIUS) {
+            x0[0] = (genrand_real2()*2 - 1)*AV_RADIUS;
+            x0[1] = (genrand_real2()*2 - 1)*AV_RADIUS;
+            x0[2] = (genrand_real2()*2 - 1)*AV_RADIUS;
+        }
+
         /* Zero out all PMTs. */
         for (j = 0; j < LEN(ev.pmt_hits); j++) {
             ev.pmt_hits[j].hit = 0;
@@ -1573,7 +1582,7 @@ int test_quad_noise(char *err)
             ev.pmt_hits[hits[j]].t = t0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT + randn()*PMT_TTS;
         }
 
-        if (quad(&ev, fit_pos, &fit_t0, 10000)) {
+        if (quad(&ev, fit_pos, &fit_t0, 10000, 1.0)) {
             sprintf(err, "%s", quad_err);
             goto err;
         }
@@ -2129,6 +2138,62 @@ err:
     return 1;
 }
 
+/* Tests that the ran_choose_weighted() function returns elements proportional
+ * to the weights. */
+int test_ran_choose_weighted(char *err)
+{
+    size_t i, j, k;
+
+    int choices[10] = {0,1,2,3,4,5,6,7,8,9};
+    int results[1000000];
+    double w[10] = {0.0,0.0,0.0,0.0,0.0,0.2,0.2,0.2,0.2,0.2};
+    double fraction[10];
+    double sum;
+
+    for (j = 0; j < 100; j++) {
+        sum = 0.0;
+        for (i = 0; i < LEN(choices); i++) {
+            w[i] = genrand_real2();
+            sum += w[i];
+        }
+
+        for (i = 0; i < LEN(choices); i++) {
+            w[i] /= sum;
+        }
+
+        for (i = 0; i < LEN(results); i++) {
+            ran_choose_weighted(results+i,w,1,choices,LEN(choices));
+        }
+
+        sum = 0.0;
+        for (i = 0; i < LEN(choices); i++) {
+            fraction[i] = 0.0;
+            for (k = 0; k < LEN(results); k++) {
+                if (results[k] == choices[i]) {
+                    fraction[i] += 1.0;
+                }
+            }
+            sum += fraction[i];
+        }
+
+        for (i = 0; i < LEN(choices); i++) {
+            fraction[i] /= sum;
+        }
+
+        for (i = 0; i < LEN(choices); i++) {
+            if (!isclose(fraction[i],w[i],0,1e-2)) {
+                sprintf(err, "ran_choose_weighted() returned %zu %.5g%% of the time but expected %.5g", i, fraction[i], w[i]);
+                goto err;
+            }
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -2182,6 +2247,7 @@ struct tests {
     {test_electron_get_angular_pdf_norm, "test_electron_get_angular_pdf_norm"},
     {test_fast_acos, "test_fast_acos"},
     {test_get_most_likely_mean_pe, "test_get_most_likely_mean_pe"},
+    {test_ran_choose_weighted, "test_ran_choose_weighted"},
 };
 
 int main(int argc, char **argv)