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, 40 insertions, 0 deletions

diff --git a/src/random.c b/src/random.c
index 0112ebd..482e102 100644
--- a/src/random.c
+++ b/src/random.c
@@ -16,6 +16,10 @@
 
 #include "random.h"
 #include <math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+static gsl_rng *r;
 
 double randn(void)
 {
@@ -44,3 +48,39 @@ void rand_sphere(double *dir)
     dir[1] = sin(theta)*sin(phi);
     dir[2] = cos(theta);
 }
+
+/* Randomly choose `k` elements from array `src` with `n` elements with
+ * probabilities proportional to `w`. The elements are sampled without
+ * replacement and placed into the `dest` array. */
+void ran_choose_weighted(int *dest, double *w, size_t k, int *src, size_t n)
+{
+    int i, j;
+    gsl_ran_discrete_t *g;
+    size_t index, result;
+    int duplicate;
+
+    if (!r)
+        r = gsl_rng_alloc(gsl_rng_default);
+
+    g = gsl_ran_discrete_preproc(n, w);
+
+    i = 0;
+    while (i < k) {
+        index = gsl_ran_discrete(r, g);
+
+        result = src[index];
+
+        duplicate = 0;
+        for (j = 0; j < i; j++) {
+            if (result == dest[j]) {
+                duplicate = 1;
+                break;
+            }
+        }
+
+        if (!duplicate) dest[i++] = result;
+    }
+
+    gsl_ran_discrete_free(g);
+}
+