add a new test for the quad fitter

This commit adds a new test to test the quad fitter when the t0 quantile
argument is less than 1.

3 files changed, 136 insertions, 29 deletions

diff --git a/src/random.c b/src/random.c
index 482e102..8f005eb 100644
--- a/src/random.c
+++ b/src/random.c
@@ -18,13 +18,14 @@
 #include <math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include "vector.h"
 
 static gsl_rng *r;
 
+/* Generates a random number from a normal distribution using the Box-Muller
+ * transform. */
 double randn(void)
 {
-    /* Generates a random number from a normal distribution using the
-     * Box-Muller transform. */
     double u1, u2;
 
     u1 = genrand_real1();
@@ -33,9 +34,23 @@ double randn(void)
     return sqrt(-2*log(u1))*cos(2*M_PI*u2);
 }
 
+/* Generates a random point in the unit sphere. */
+void rand_ball(double *pos, double radius)
+{
+    pos[0] = (genrand_real2()*2 - 1)*radius;
+    pos[1] = (genrand_real2()*2 - 1)*radius;
+    pos[2] = (genrand_real2()*2 - 1)*radius;
+
+    while (NORM(pos) >= radius) {
+        pos[0] = (genrand_real2()*2 - 1)*radius;
+        pos[1] = (genrand_real2()*2 - 1)*radius;
+        pos[2] = (genrand_real2()*2 - 1)*radius;
+    }
+}
+
+/* Generates a random point on the unit sphere. */
 void rand_sphere(double *dir)
 {
-    /* Generates a random point on the unit sphere. */
     double u, v, theta, phi;
 
     u = genrand_real1();
diff --git a/src/random.h b/src/random.h
index 747cfeb..e2cba7d 100644
--- a/src/random.h
+++ b/src/random.h
@@ -21,6 +21,7 @@
 #include <stdlib.h> /* for size_t */
 
 double randn(void);
+void rand_ball(double *pos, double radius);
 void rand_sphere(double *dir);
 void ran_choose_weighted(int *dest, double *w, size_t k, int *src, size_t n);
 
diff --git a/src/test.c b/src/test.c
index a0b22ae..67b15bb 100644
--- a/src/test.c
+++ b/src/test.c
@@ -1440,19 +1440,8 @@ int test_quad(char *err)
     n_d2o = get_avg_index_d2o();
 
     for (i = 0; i < 100; i++) {
-        /* Generate a random position within a cube the size of the AV.
-         *
-         * Note: This does produce some points which are outside of the PSUP,
-         * but I think the quad fitter should still be able to fit these. */
-        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;
-        }
+        /* Generate a random position within a cube the size of the AV. */
+        rand_ball(x0,AV_RADIUS);
 
         t0 = genrand_real2()*GTVALID;
 
@@ -1550,20 +1539,10 @@ int test_quad_noise(char *err)
     n_d2o = get_avg_index_d2o();
 
     for (i = 0; i < 100; i++) {
-        /* Generate a random position within a cube the size of the AV.
-         *
-         * Note: This does produce some points which are outside of the PSUP,
-         * but I think the quad fitter should still be able to fit these. */
-        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;
+        /* Generate a random position within a cube the size of the AV. */
+        rand_ball(x0,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. */
         for (j = 0; j < LEN(ev.pmt_hits); j++) {
@@ -2194,6 +2173,117 @@ err:
     return 1;
 }
 
+/* Tests that the quad fitter works when we have multiple vertices. */
+int test_quad_muon(char *err)
+{
+    size_t i, j;
+    double x0[3], pos2[3];
+    double t0;
+    const gsl_rng_type *T;
+    gsl_rng *r;
+    event ev;
+    int index[MAX_PMTS];
+    int hits[1000];
+    size_t valid_pmts;
+    double fit_pos[3];
+    double pmt_dir[3];
+    double fit_t0;
+    double n_d2o;
+
+    init_genrand(0);
+
+    T = gsl_rng_default;
+    r = gsl_rng_alloc(T);
+
+    load_pmt_info();
+
+    valid_pmts = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (pmts[i].pmt_type != PMT_NORMAL) continue;
+
+        index[valid_pmts++] = i;
+    }
+
+    n_d2o = get_avg_index_d2o();
+
+    for (i = 0; i < 100; i++) {
+        /* Generate a random position within the AV. */
+        rand_ball(x0,AV_RADIUS);
+
+        t0 = genrand_real2()*GTVALID;
+
+        /* Zero out all PMTs. */
+        for (j = 0; j < LEN(ev.pmt_hits); j++) {
+            ev.pmt_hits[j].hit = 0;
+            ev.pmt_hits[j].flags = 0;
+        }
+
+        /* Choose LEN(hits) random PMTs which are hit. */
+        gsl_ran_choose(r,hits,LEN(hits),index,valid_pmts,sizeof(int));
+
+        /* Calculate the time the PMT got hit. */
+        for (j = 0; j < LEN(hits); j++) {
+            SUB(pmt_dir,pmts[hits[j]].pos,x0);
+            ev.pmt_hits[hits[j]].hit = 1;
+            ev.pmt_hits[hits[j]].q = genrand_real2();
+            ev.pmt_hits[hits[j]].t = t0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT;
+        }
+
+        /* Choose another position. */
+        rand_ball(pos2,AV_RADIUS);
+
+        /* Choose LEN(hits) random PMTs which are hit. */
+        gsl_ran_choose(r,hits,LEN(hits),index,valid_pmts,sizeof(int));
+
+        /* Calculate the time the PMT got hit. */
+        for (j = 0; j < LEN(hits); j++) {
+            SUB(pmt_dir,pmts[hits[j]].pos,pos2);
+            ev.pmt_hits[hits[j]].hit = 1;
+            ev.pmt_hits[hits[j]].q = genrand_real2();
+            /* Assume this vertex occurs 100 ns later. */
+            ev.pmt_hits[hits[j]].t = t0 + 100.0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT;
+        }
+
+        if (quad(&ev, fit_pos, &fit_t0, 10000, 0.1)) {
+            sprintf(err, "%s", quad_err);
+            goto err;
+        }
+
+        /* Since there's no noise, these results should be exact. We test to
+         * see that all of the positions are within 1 cm and the time is within
+         * 1 ns. */
+
+        if (!isclose(fit_pos[0], x0[0], 0, 1.0)) {
+            sprintf(err, "quad returned x = %.5f, but expected %.5f", fit_pos[0], x0[0]);
+            goto err;
+        }
+
+        if (!isclose(fit_pos[1], x0[1], 0, 1.0)) {
+            sprintf(err, "quad returned y = %.5f, but expected %.5f", fit_pos[1], x0[1]);
+            goto err;
+        }
+
+        if (!isclose(fit_pos[2], x0[2], 0, 1.0)) {
+            sprintf(err, "quad returned z = %.5f, but expected %.5f", fit_pos[2], x0[2]);
+            goto err;
+        }
+
+        if (!isclose(fit_t0, t0, 0, 1.0)) {
+            sprintf(err, "quad returned t0 = %.5f, but expected %.5f", fit_t0, t0);
+            goto err;
+        }
+    }
+
+    gsl_rng_free(r);
+
+    return 0;
+
+err:
+    gsl_rng_free(r);
+
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -2248,6 +2338,7 @@ struct tests {
     {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"},
+    {test_quad_muon, "test_quad_muon"},
 };
 
 int main(int argc, char **argv)