diff options
Diffstat (limited to 'src/find_peaks.c')
| -rw-r--r-- | src/find_peaks.c | 106 |
1 files changed, 56 insertions, 50 deletions
diff --git a/src/find_peaks.c b/src/find_peaks.c index 6109896..14c33e5 100644 --- a/src/find_peaks.c +++ b/src/find_peaks.c @@ -24,6 +24,9 @@ #include <math.h> /* for exp() */ #include "misc.h" #include "pmt.h" +#include "sno_charge.h" +#include <gsl/gsl_randist.h> +#include "pdg.h" typedef struct peak { size_t i; @@ -133,33 +136,7 @@ end: return; } -double get_qhs_avg(event *ev, double *pos, double *dir) -{ - /* Returns the average QHS value for all PMTs within the Cerenkov cone of a - * particle at position `pos` travelling in direction `dir`. */ - size_t i; - double qhs_sum, n_d2o; - double pmt_dir[3]; - size_t n; - - n_d2o = get_index_snoman_d2o(400.0); - - n = 0; - qhs_sum = 0.0; - for (i = 0; i < MAX_PMTS; i++) { - if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue; - SUB(pmt_dir,pmts[i].pos,pos); - normalize(pmt_dir); - if (DOT(pmt_dir,dir) > 1/n_d2o) { - qhs_sum += ev->pmt_hits[i].q; - n += 1; - } - } - - return qhs_sum/n; -} - -void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, size_t m, double *result, double *last, size_t len) +void get_hough_transform(event *ev, double *pos, double t0, double *x, double *y, size_t n, size_t m, double *result, double *last, size_t len) { /* Computes the "Hough transform" of the event `ev` and stores it in `result`. */ size_t i, j, k; @@ -167,7 +144,11 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, int skip; double *sin_theta, *cos_theta, *sin_phi, *cos_phi; double wavelength0, n_d2o; - double qhs[MAX_PMTS], qhs_avg; + static double w[MAX_PMTS]; + double expected_pe; + int nhit; + double pq; + double dt; sin_theta = malloc(sizeof(double)*n); cos_theta = malloc(sizeof(double)*n); @@ -189,27 +170,46 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, cos_phi[i] = cos(y[i]); } + 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 || !ev->pmt_hits[i].hit) continue; - qhs[i] = ev->pmt_hits[i].q; - } - - /* 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[j].flags || pmts[j].pmt_type != PMT_NORMAL || !ev->pmt_hits[j].hit) continue; - - SUB(pmt_dir,pmts[j].pos,pos); - - normalize(pmt_dir); - - cos_theta2 = DOT(pmt_dir,last+3*i); + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; - qhs[j] -= qhs_avg*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1); + if (ev->pmt_hits[i].hit) { + expected_pe += ev->pmt_hits[i].q; + nhit += 1; + } + } + expected_pe /= nhit; - if (qhs[j] < 0.0) qhs[j] = 0.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) { + /* 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 > FIND_PEAKS_MAX_PE*get_qmean()/2) { + /* If the charge is greater than FIND_PEAKS_MAX_PE, it's almost + * certainly multiple photons so we don't waste time + * calculating P(1 PE|q). */ + w[i] = 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 < FIND_PEAKS_MAX_PE; j++) { + pq += get_pq(ev->pmt_hits[i].q,j)*gsl_ran_poisson_pdf(j,expected_pe); + } + + w[i] = 1-get_pq(ev->pmt_hits[i].q,1)*gsl_ran_poisson_pdf(1,expected_pe)/pq; + } } } @@ -221,12 +221,18 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, SUB(pmt_dir,pmts[i].pos,pos); + dt = ev->pmt_hits[i].t - t0 - NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT; + + /* Skip PMT hits with a time residual of more than 10 nanoseconds since + * they are most likely to be scattered and/or reflected light. */ + if (dt > 10) continue; + normalize(pmt_dir); /* Ignore PMT hits within the Cerenkov cone of previously found peaks. */ skip = 0; for (j = 0; j < len; j++) { - if (DOT(pmt_dir,last+3*j) > 1/n_d2o) { + if (fabs(DOT(pmt_dir,last+3*j) - 1/n_d2o) < 0.1) { skip = 1; break; } @@ -242,7 +248,7 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, cos_theta2 = DOT(pmt_dir,dir); - result[j*n + k] += qhs[i]*exp(-fabs(cos_theta2-1.0/n_d2o)/0.1); + result[j*n + k] += w[i]*exp(-pow(cos_theta2-1.0/n_d2o,2)/0.01); } } } @@ -268,7 +274,7 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, * returned. So, for example if three peaks are found but the first two are * very close to each other, this function will only return the first and the * third peaks. */ -void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, double *peak_phi, size_t *npeaks, size_t max_peaks, double delta) +void find_peaks(event *ev, double *pos, double t0, size_t n, size_t m, double *peak_theta, double *peak_phi, size_t *npeaks, size_t max_peaks, double delta) { size_t i, j; double *x, *y, *result, *dir; @@ -296,7 +302,7 @@ void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, *npeaks = 0; for (i = 0; i < max_peaks; i++) { - get_hough_transform(ev,pos,x,y,n,m,result,dir,i); + get_hough_transform(ev,pos,t0,x,y,n,m,result,dir,i); max = argmax(result,n*m); theta = x[max/m]; phi = y[max % m]; |