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.

11 files changed, 317 insertions, 56 deletions

diff --git a/src/Makefile b/src/Makefile
index 7787558..cb701cb 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -32,9 +32,9 @@ test-time-pdf: test-time-pdf.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o
 
 test-zebra: test-zebra.o zebra.o pack2b.o
 
-fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o quad.o dc.o sort.o util.o hdf5_utils.o
+fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o quad.o dc.o sort.o util.o hdf5_utils.o random.o mt19937ar.o
 
-test-find-peaks: test-find-peaks.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o util.o quad.o
+test-find-peaks: test-find-peaks.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o util.o quad.o random.o mt19937ar.o
 
 calculate-csda-range: calculate-csda-range.o
 
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;
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;
diff --git a/src/quad.c b/src/quad.c
index 0c1f349..ca5483c 100644
--- a/src/quad.c
+++ b/src/quad.c
@@ -31,6 +31,10 @@
 #include <gsl/gsl_cblas.h>
 #include "misc.h"
 #include <gsl/gsl_errno.h>
+#include "sno_charge.h"
+#include "sno.h"
+#include <gsl/gsl_permute.h>
+#include "random.h"
 
 char quad_err[256];
 
@@ -54,6 +58,11 @@ void my_handler(const char *reason, const char *file, int line, int gsl_errno)
  * density of points, but since we are dealing with mostly high nhit events we
  * just take the median because it is easier and quicker.
  *
+ * Update: This version of quad now also weights the PMT hits it selects by the
+ * probability that they are more than 1 PE. The idea here is that we'd like to
+ * ignore scattered and reflected light when selecting these points and most
+ * scattered and reflected light is single photons.
+ *
  * `ev` should be a pointer to the event.
  *
  * `pos` is a pointer to a double array of length 3.
@@ -64,6 +73,15 @@ void my_handler(const char *reason, const char *file, int line, int gsl_errno)
  * try to compute this many points, but it will stop trying after 10*npoints
  * times to avoid an infinite loop.
  *
+ * `f` is the quantile of t0 to cut on. It should be between 0 and 1. The
+ * reason to cut on a quantile of t0 is that for particles with tracks much
+ * longer than a centimeter, the quad cloud of points generally follows the
+ * whole track. Since we are interested in finding the position and time of the
+ * start of the track, we'd like to only select the quad points at the start of
+ * the track. To do this we only include quad cloud points in the first `f`
+ * quantile when computing the position and time. For the default quad behavior
+ * without cutting on t0, you should set `f` to 1.0.
+ *
  * Returns 0 on success and `pos` and `t0` will be set to the approximate
  * vertex position and time respectively. On error, returns -1 and sets the
  * `quad_err` string.
@@ -80,11 +98,10 @@ void my_handler(const char *reason, const char *file, int line, int gsl_errno)
  * is Stephen Brice's PHD thesis which cites three SNO technical reports, but I
  * haven't been able to find these. The first two of these reports are by Bill
  * Frati who I assume wrote the initial implementation. */
-int quad(event *ev, double *pos, double *t0, size_t npoints)
+int quad(event *ev, double *pos, double *t0, size_t npoints, double f)
 {
     size_t i, j, k;
     static int index[MAX_PMTS];
-    static int index2[MAX_PMTS];
     int nhit;
     const gsl_rng_type *T;
     gsl_rng *r;
@@ -96,21 +113,20 @@ int quad(event *ev, double *pos, double *t0, size_t npoints)
     double pos1[3], pos2[3];
     static double results[QUAD_MAX][4];
     double c2;
-    double wavelength0;
     double n_d2o;
     double t1, t2;
     int s;
     double pmt_dir[3];
     double expected;
-    static double t[MAX_PMTS];
-    static size_t p2[MAX_PMTS];
+    static double w[MAX_PMTS];
     double tmin;
     size_t max_tries;
     gsl_error_handler_t *old_handler;
     int status;
+    size_t reorder[QUAD_MAX];
+    double pq = 0.0;
 
-    wavelength0 = 400.0;
-    n_d2o = get_index_snoman_d2o(wavelength0);
+    n_d2o = get_avg_index_d2o();
 
     c2 = SPEED_OF_LIGHT*SPEED_OF_LIGHT/pow(n_d2o,2);
 
@@ -119,31 +135,54 @@ int quad(event *ev, double *pos, double *t0, size_t npoints)
         return 1;
     }
 
+    double expected_pe = 0.0;
     nhit = 0;
     for (i = 0; i < MAX_PMTS; i++) {
         if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
 
         if (ev->pmt_hits[i].hit) {
-            index[nhit] = i;
-            t[nhit] = ev->pmt_hits[i].t;
-            nhit++;
+            expected_pe += ev->pmt_hits[i].q;
+            nhit += 1;
         }
     }
+    expected_pe /= nhit;
 
-    /* Only use the first 50% of the hits to avoid counting late/reflected
-     * light. */
-    k = nhit > 10 ? ceil(nhit/2.0) : nhit;
-    gsl_sort_smallest_index(p2,k,t,1,nhit);
+    nhit = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
 
-    /* Reorder the hits. */
-    for (i = 0; i < k; i++) {
-        index2[i] = index[p2[i]];
-    }
+        if (ev->pmt_hits[i].hit) {
+            index[nhit] = i;
 
-    nhit = k;
+            /* Weights are equal to the probability that the hit is > 1 photon.
+             * This is supposed to be a rough proxy for the probability that
+             * it's not reflected and/or scattered light (which is typically
+             * single photons). */
+            if (ev->pmt_hits[i].q > QUAD_MAX_PE*get_qmean()/2) {
+                /* If the charge is greater than QUAD_MAX_PE, it's almost
+                 * certainly multiple photons so we don't waste time calculating P(1 PE|q). */
+                w[nhit] = 1.0;
+            } else {
+                /* We want to calculate P(multiple photons|q) which we calculate as:
+                 *
+                 *     P(multiple photons|q) = 1 - P(1 PE|q) = 1 - P(q|1 PE)P(1 PE)/P(q)
+                 *
+                 * To calculate the second two quantities we assume the number
+                 * of PE is Poisson distributed with a mean equal to
+                 * expected_photons. */
+                pq = 0.0;
+                for (j = 1; j < QUAD_MAX_PE; j++) {
+                    pq += get_pq(ev->pmt_hits[i].q,j)*gsl_ran_poisson_pdf(j,expected_pe);
+                }
+
+                w[nhit] = 1-get_pq(ev->pmt_hits[i].q,1)*gsl_ran_poisson_pdf(1,expected_pe)/pq;
+            }
+            nhit++;
+        }
+    }
 
     if (nhit < 5) {
-        sprintf(quad_err, "only %i pmt hits. quad needs at least 5 points!", nhit);
+        sprintf(quad_err, "only %i pmt hit(s). quad needs at least 5 points!", nhit);
         return 1;
     }
 
@@ -157,7 +196,7 @@ int quad(event *ev, double *pos, double *t0, size_t npoints)
     nresults = 0;
     while (nresults < npoints && i++ < max_tries) {
         /* Choose 4 random hits. */
-        gsl_ran_choose(r,xs,4,index2,nhit,sizeof(int));
+        ran_choose_weighted(xs,w,4,index,nhit);
 
         /* Shuffle them since GSL always returns the random choices in order.
          *
@@ -224,11 +263,13 @@ int quad(event *ev, double *pos, double *t0, size_t npoints)
         }
 
         if (t1 < tmin && isclose(ev->pmt_hits[xs[0]].t,expected,0,1e-2)) {
-            results[nresults][0] = pos1[0];
-            results[nresults][1] = pos1[1];
-            results[nresults][2] = pos1[2];
-            results[nresults][3] = t1;
-            nresults++;
+            if (NORM(pos1) < PSUP_RADIUS) {
+                results[nresults][0] = pos1[0];
+                results[nresults][1] = pos1[1];
+                results[nresults][2] = pos1[2];
+                results[nresults][3] = t1;
+                nresults++;
+            }
         } else {
             /* Check the second solution. */
             t2 = (-b - sqrt(b*b - 4*a*c))/(2*a);
@@ -240,11 +281,13 @@ int quad(event *ev, double *pos, double *t0, size_t npoints)
             expected = t2 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT;
 
             if (t2 < tmin && isclose(ev->pmt_hits[xs[0]].t,expected,0,1e-2)) {
-                results[nresults][0] = pos2[0];
-                results[nresults][1] = pos2[1];
-                results[nresults][2] = pos2[2];
-                results[nresults][3] = t2;
-                nresults++;
+                if (NORM(pos2) < PSUP_RADIUS) {
+                    results[nresults][0] = pos2[0];
+                    results[nresults][1] = pos2[1];
+                    results[nresults][2] = pos2[2];
+                    results[nresults][3] = t2;
+                    nresults++;
+                }
             }
         }
     }
@@ -254,6 +297,27 @@ int quad(event *ev, double *pos, double *t0, size_t npoints)
         goto err;
     }
 
+    /* Compute the permutation required to sort the results by t0. */
+    gsl_sort_index(reorder,&results[0][3],4,nresults);
+
+    /* Sort the results by t0. */
+    gsl_permute(reorder,&results[0][0],4,nresults);
+    gsl_permute(reorder,&results[0][1],4,nresults);
+    gsl_permute(reorder,&results[0][2],4,nresults);
+    gsl_permute(reorder,&results[0][3],4,nresults);
+
+    if (f > 0.0 && f < 1.0) {
+        /* Now, we filter only the results with t0 less than the quantile given
+         * by `f`. The idea here is that for high energy particles which travel
+         * macroscopic distances in the detector we want to only sample the
+         * quad points near the start of the track, i.e. the points with the
+         * earliest times. */
+        nresults = (int) f*nresults;
+    }
+
+    /* Sort the x, y, z, and t0 columns so we can calculate the median. Note:
+     * The rows of the results array don't represent the quad cloud anymore
+     * since x, y, z, and t0 are mixed up! */
     gsl_sort(&results[0][0],4,nresults);
     gsl_sort(&results[0][1],4,nresults);
     gsl_sort(&results[0][2],4,nresults);
diff --git a/src/quad.h b/src/quad.h
index 51ff6ee..e95ae37 100644
--- a/src/quad.h
+++ b/src/quad.h
@@ -19,12 +19,16 @@
 
 #include "event.h"
 #include <stdlib.h> /* for size_t */
+#include <gsl/gsl_rng.h>
 
 /* Maximum number of quad cloud points. */
 #define QUAD_MAX 10000
 
+/* Maximum number of PE to consider when calculating P(> 1 hit|q) */
+#define QUAD_MAX_PE 10
+
 extern char quad_err[256];
 
-int quad(event *ev, double *pos, double *t0, size_t npoints);
+int quad(event *ev, double *pos, double *t0, size_t npoints, double f);
 
 #endif
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);
+}
+
diff --git a/src/random.h b/src/random.h
index dac48f8..747cfeb 100644
--- a/src/random.h
+++ b/src/random.h
@@ -18,8 +18,10 @@
 #define RANDOM_H
 
 #include "mt19937ar.h"
+#include <stdlib.h> /* for size_t */
 
 double randn(void);
 void rand_sphere(double *dir);
+void ran_choose_weighted(int *dest, double *w, size_t k, int *src, size_t n);
 
 #endif
diff --git a/src/sno_charge.c b/src/sno_charge.c
index f0d38ec..130df9f 100644
--- a/src/sno_charge.c
+++ b/src/sno_charge.c
@@ -116,6 +116,26 @@ double get_qhi(void)
     return qhi;
 }
 
+double get_qmean(void)
+{
+    if (!initialized) {
+        fprintf(stderr, "charge interpolation hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return qmean;
+}
+
+double get_qstd(void)
+{
+    if (!initialized) {
+        fprintf(stderr, "charge interpolation hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return qstd;
+}
+
 static double nlopt_log_pq(unsigned int n, const double *x, double *grad, void *params)
 {
     int i;
diff --git a/src/sno_charge.h b/src/sno_charge.h
index 7c6513d..b54f844 100644
--- a/src/sno_charge.h
+++ b/src/sno_charge.h
@@ -19,6 +19,8 @@
 
 double get_qlo(void);
 double get_qhi(void);
+double get_qmean(void);
+double get_qstd(void);
 double get_most_likely_mean_pe(double q);
 void init_charge(void);
 double get_log_pq(double q, int n);
diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c
index 11fa4b4..2ccc90d 100644
--- a/src/test-find-peaks.c
+++ b/src/test-find-peaks.c
@@ -31,12 +31,16 @@
 #include "quad.h"
 #include "likelihood.h"
 #include "optics.h"
+#include "sno_charge.h"
+#include "db.h"
+#include "pmt_response.h"
+#include "dqxx.h"
 
 #define EV_RECORD      0x45562020
 #define MCTK_RECORD    0x4d43544b
 #define MCVX_RECORD    0x4d435658
 
-#define NUM_ANGLES 1000
+#define TEST_FIND_PEAKS_NUM_ANGLES 1000
 
 void plot_hough_transform(double *x, double *y, double *z, size_t n, size_t m)
 {
@@ -101,6 +105,8 @@ void plot_find_peaks(event *ev, double *pos, double *peak_theta, double *peak_ph
     fprintf(pipe, "plot '-' u 1:2:3:4 with circles palette fillstyle solid, '-' u 1:2 with lines lc rgb \"red\" lw 2\n");
 
     for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+
         if (ev->pmt_hits[i].hit) {
             r = NORM(pmts[i].pos);
             theta = acos(pmts[i].pos[2]/r);
@@ -131,9 +137,9 @@ void plot_find_peaks(event *ev, double *pos, double *peak_theta, double *peak_ph
         /* Rotate the direction by the Cerenkov angle. */
         rotate(tmp2,dir,tmp,acos(1/n_d2o));
 
-        for (j = 0; j < NUM_ANGLES; j++) {
+        for (j = 0; j < TEST_FIND_PEAKS_NUM_ANGLES; j++) {
             /* Rotate the new direction around the peak by 2 pi. */
-            rotate(tmp,tmp2,dir,j*2*M_PI/(NUM_ANGLES-1));
+            rotate(tmp,tmp2,dir,j*2*M_PI/(TEST_FIND_PEAKS_NUM_ANGLES-1));
 
             /* Calculate the intersection of the ray with the sphere. */
             if (!intersect_sphere(pos,tmp,PSUP_RADIUS,&l)) continue;
@@ -197,6 +203,8 @@ int main(int argc, char **argv)
     int plot = 0;
     int status;
     double t0;
+    int last_run;
+    char dqxx_file[256];
     int32_t gtid = -1;
 
     for (i = 1; i < argc; i++) {
@@ -233,6 +241,8 @@ int main(int argc, char **argv)
 
     load_pmt_info();
 
+    init_charge();
+
     f = zebra_open(filename);
 
     if (!f) {
@@ -240,6 +250,30 @@ int main(int argc, char **argv)
         return 1;
     }
 
+    dict *db = db_init();
+
+    last_run = -1;
+
+    if (load_file(db, "pmt_response_qoca_d2o_20060216.dat", 0)) {
+        fprintf(stderr, "failed to load pmt_response_qoca_d2o_20060216.dat: %s\n", db_err);
+        goto err;
+    }
+
+    if (load_file(db, "rsp_rayleigh.dat", 0)) {
+        fprintf(stderr, "failed to load rsp_rayleigh.dat: %s\n", db_err);
+        goto err;
+    }
+
+    if (pmt_response_init(db)) {
+        fprintf(stderr, "failed to initialize PMTR bank: %s\n", pmtr_err);
+        goto err;
+    }
+
+    if (optics_init()) {
+        fprintf(stderr, "failed to initialize optics: %s\n", optics_err);
+        goto err;
+    }
+
     while (1) {
         rv = zebra_read_next_logical_record(f);
 
@@ -251,7 +285,12 @@ int main(int argc, char **argv)
             break;
         }
 
-        rv = zebra_get_bank(f, &bmast, f->first_bank);
+        if (f->mast_bank == -1) {
+            fprintf(stderr, "no MAST bank in logical record! Skipping...\n");
+            continue;
+        }
+
+        rv = zebra_get_bank(f, &bmast, f->mast_bank);
 
         if (rv) {
             fprintf(stderr, "error getting MAST bank: %s\n", zebra_err);
@@ -288,6 +327,23 @@ int main(int argc, char **argv)
 
         if (gtid > 0 && ev.gtid != gtid) continue;
 
+        if (bev.run != last_run) {
+            printf("loading DQXX file for run %010i\n", bev.run);
+
+            sprintf(dqxx_file, "DQXX_%010i.dat", bev.run);
+            if (load_file(db, dqxx_file, 1)) {
+                fprintf(stderr, "failed to load %s: %s\n", dqxx_file, db_err);
+                goto err;
+            }
+
+            if (dqxx_init(db, &ev)) {
+                fprintf(stderr, "failed to initialize DQXX bank: %s\n", dqxx_err);
+                goto err;
+            }
+
+            last_run = bev.run;
+        }
+
         rv = get_event(f,&ev,&b);
 
         size_t n = 100;
@@ -306,7 +362,7 @@ int main(int argc, char **argv)
         }
 
         /* 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 < 0 || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) {
             /* If the QUAD fitter fails or returns something outside the PSUP or
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)