From d1ce12ce69604c6979f4ffa45e4908eb71b19c60 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 13 Apr 2020 11:40:13 -0500 Subject: update find_peaks algorithm This commit updates the find peaks algorithm with several improvements which together drastically improve its ability to find Cerenkov rings: - when computing the Hough transform, instead of charge we weight each PMT hit by the probability that it is a multi-photon PMT hit - we don't subtract off previously found rings (this makes the code simpler and I don't think it previously had a huge effect) - ignore PMT hits who are within approximately 5 degrees of any previously found ring (previously we ignored all hits within the center of previously found rings) - ignore PMT hits which have a time residual of more than 10 nanoseconds to hopefully ignore more reflected and/or scattered light - switch from weighting the Hough transform by exp(-fabs(cos(theta)-1/n)/0.1) -> exp(-pow(cos(theta)-1/n,2)/0.01). I'm still not sure if this has a huge effect, but the reason I switched is that the PDF for Cerenkov light looks closer to the second form. - switch to calling quad with f = 1.0 in test-find-peaks (I still need to add this update to fit.c but will do that in a later commit). --- src/find_peaks.c | 106 ++++++++++++++++++++++++++------------------------ src/find_peaks.h | 7 +++- src/fit.c | 2 +- src/test-find-peaks.c | 6 +-- 4 files changed, 65 insertions(+), 56 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 /* for exp() */ #include "misc.h" #include "pmt.h" +#include "sno_charge.h" +#include +#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]; diff --git a/src/find_peaks.h b/src/find_peaks.h index 0219561..584774a 100644 --- a/src/find_peaks.h +++ b/src/find_peaks.h @@ -20,8 +20,11 @@ #include "event.h" #include /* for size_t */ +/* Maximum number of PE to consider when calculating P(> 1 hit|q) */ +#define FIND_PEAKS_MAX_PE 10 + void find_peaks_array(double *x, size_t n, size_t m, size_t *imax, size_t *jmax, size_t *npeaks, size_t max_peaks, double threshold); -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 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 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); +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); #endif diff --git a/src/fit.c b/src/fit.c index 0fa8766..6c6b66d 100644 --- a/src/fit.c +++ b/src/fit.c @@ -5488,7 +5488,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double ub[3] = GTVALID; /* Find the peaks in the Hough transform of the event. */ - find_peaks(ev, pos, 100, 100, peak_theta, peak_phi, &npeaks, LEN(peak_theta), 0.1); + find_peaks(ev, pos, t0, 100, 100, peak_theta, peak_phi, &npeaks, LEN(peak_theta), 0.1); /* Don't fit more than 3 peaks for now. */ if (npeaks > 3) npeaks = 3; diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c index 2ccc90d..fd8d0cd 100644 --- a/src/test-find-peaks.c +++ b/src/test-find-peaks.c @@ -362,7 +362,7 @@ int main(int argc, char **argv) } /* Guess the position and t0 of the event using the QUAD fitter. */ - status = quad(&ev,pos,&t0,10000,0.1); + status = quad(&ev,pos,&t0,10000,1.0); if (status || t0 < 0 || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) { /* If the QUAD fitter fails or returns something outside the PSUP or @@ -374,9 +374,9 @@ int main(int argc, char **argv) pos[2] = 0.0; } - get_hough_transform(&ev,pos,x,y,n,m,result,0,0); + get_hough_transform(&ev,pos,t0,x,y,n,m,result,0,0); - find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,max_npeaks,0.1); + find_peaks(&ev,pos,t0,n,m,peak_theta,peak_phi,&npeaks,max_npeaks,0.1); printf("gtid %i\n", ev.gtid); for (i = 0; i < npeaks; i++) { -- cgit