From b78b16d32f4ed170f4a24d290348765da198d5f0 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Wed, 6 Nov 2019 11:28:27 -0600 Subject: add a couple of improvements to the quad fitter and fix a bug in get_hough_transform() This commit adds two improvements to the quad fitter: 1. I updated quad to weight the random PMT hit selection by the probability that the PMT hit is a multiphoton hit. The idea here is that we really only want to sample direct light and for high energy events the reflected and scattered light is usually single photon. 2. I added an option to quad to only use points in the quad cloud which are below a given quantile of t0. The idea here is that for particles like muons which travel more than a few centimeters in the detector the quad cloud usually looks like the whole track. Since we want the QUAD fitter to find the position of the *start* of the track we select only those quad cloud points with an early time so the position is closer to the position of the start of the track. Also, I fixed a major bug in get_hough_transform() in which I was using the wrong index variable when checking if a PMT was not flagged, a normal PMT, and was hit. This was causing the algorithm to completely miss finding more than one ring while I was testing it. --- src/find_peaks.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'src/find_peaks.c') diff --git a/src/find_peaks.c b/src/find_peaks.c index 7d350d5..6109896 100644 --- a/src/find_peaks.c +++ b/src/find_peaks.c @@ -197,8 +197,9 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, /* 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[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue; + 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); @@ -304,7 +305,7 @@ void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, /* Only add the latest peak to the results if it's more than `delta` * radians away from all other results so far. */ unique = 1; - for (j = 0; j < i; j++) { + for (j = 0; j < *npeaks; j++) { /* If the angle between the latest peak and a previous peak is * within `delta` radians, it's not unique. */ if (acos(DOT(dir+i*3,dir+j*3)) < delta) unique = 0; -- cgit