aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/random.c21
-rw-r--r--src/random.h1
-rw-r--r--src/test.c143
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)
generated by cgit (git 2.34.1) at 2025-05-06 16:15:31 +0000