diff options
Diffstat (limited to 'src/test.c')
| -rw-r--r-- | src/test.c | 210 |
1 files changed, 210 insertions, 0 deletions
@@ -22,6 +22,8 @@ #include <gsl/gsl_errno.h> /* for gsl_strerror() */ #include <gsl/gsl_randist.h> #include <gsl/gsl_cdf.h> +#include "sno.h" +#include "quad.h" typedef int testFunction(char *err); @@ -1358,6 +1360,212 @@ err: return 1; } +int test_quad(char *err) +{ + /* Tests the quad fitter without noise. We draw 1000 random hits, compute + * the time each PMT is hit (without noise) and then make sure that the + * positions returned by the quad fitter are all within 1 cm and the + * initial time is within 1 ns. */ + size_t i, j; + double x0[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 wavelength0; + 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; + } + + wavelength0 = 400.0; + n_d2o = get_index_snoman_d2o(wavelength0); + + 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; + + /* 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]].qhs = genrand_real2(); + ev.pmt_hits[hits[j]].t = t0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT; + } + + if (quad(&ev, fit_pos, &fit_t0, 10000)) { + sprintf(err, 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; +} + +int test_quad_noise(char *err) +{ + /* Tests the quad fitter with noise. We draw 1000 random hits, compute + * the time each PMT is hit, add a gaussian sample with standard deviation + * PMT_TTS and then make sure that the positions returned by the quad + * fitter are all within 1 m and the initial time is within 10 ns. */ + size_t i, j; + double x0[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 wavelength0; + 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; + } + + wavelength0 = 400.0; + n_d2o = get_index_snoman_d2o(wavelength0); + + 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; + + /* 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]].qhs = genrand_real2(); + 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)) { + sprintf(err, quad_err); + goto err; + } + + if (!isclose(fit_pos[0], x0[0], 0, 100.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, 100.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, 100.0)) { + sprintf(err, "quad returned z = %.5f, but expected %.5f", fit_pos[2], x0[2]); + goto err; + } + + if (!isclose(fit_t0, t0, 0, 10.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; @@ -1396,6 +1604,8 @@ struct tests { {test_norm_cdf, "test_norm_cdf"}, {test_time_pdf_norm, "test_time_pdf_norm"}, {test_time_cdf, "test_time_cdf"}, + {test_quad, "test_quad"}, + {test_quad_noise, "test_quad_noise"}, }; int main(int argc, char **argv) |