diff options
Diffstat (limited to 'src/quad.c')
| -rw-r--r-- | src/quad.c | 126 |
1 files changed, 95 insertions, 31 deletions
@@ -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); |