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/test-find-peaks.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src/test-find-peaks.c') 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