diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/find_peaks.c | 56 | ||||
| -rw-r--r-- | src/find_peaks.h | 2 | ||||
| -rw-r--r-- | src/fit.c | 7 | ||||
| -rw-r--r-- | src/test-find-peaks.c | 2 |
4 files changed, 47 insertions, 20 deletions
diff --git a/src/find_peaks.c b/src/find_peaks.c index 04bf40a..135e4e3 100644 --- a/src/find_peaks.c +++ b/src/find_peaks.c @@ -252,21 +252,29 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, free(cos_phi); } -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) +/* Finds rings in the event `ev` by looking for peaks in the Hough transform. + * The directions of potential particle rings are stored in the arrays + * `peak_theta` and `peak_phi`. The number of rings found are stored in the + * variable `npeaks`. + * + * This function first computes the Hough transform of the event `ev` and looks + * for the highest peak. The next peak is then found by computing the Hough + * transform after subtracting off the previous ring and ignoring any PMTs + * within the Cerenkov cone of the first peak. This process is repeated until + * `max_peaks` peaks are found. + * + * In addition, only peaks which are `delta` radians away from other peaks are + * 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) { - /* Finds rings in the event `ev` by looking for peaks in the Hough - * transform. The directions of potential particle rings are stored in the - * arrays `peak_theta` and `peak_phi`. The number of rings found are stored - * in the variable `npeaks`. - * - * This function first computes the Hough transform of the event `ev` and - * looks for the highest peak. The next peak is then found by computing the - * Hough transform ignoring any PMTs within the Cerenkov cone of the first - * peak. This process is repeated until `max_peaks` peaks are found. */ - size_t i; + size_t i, j; double *x, *y, *result, *dir; size_t *imax, *jmax; size_t max; + double theta, phi; + int unique; x = calloc(n,sizeof(double)); y = calloc(m,sizeof(double)); @@ -284,15 +292,31 @@ void find_peaks(event *ev, double *pos, size_t n, size_t m, double *peak_theta, imax = calloc(max_peaks,sizeof(size_t)); jmax = calloc(max_peaks,sizeof(size_t)); + *npeaks = 0; + for (i = 0; i < max_peaks; i++) { get_hough_transform(ev,pos,x,y,n,m,result,dir,i); max = argmax(result,n*m); - peak_theta[i] = x[max/m]; - peak_phi[i] = y[max % m]; - get_dir(dir+i*3,peak_theta[i],peak_phi[i]); - } + theta = x[max/m]; + phi = y[max % m]; + get_dir(dir+i*3,theta,phi); + + /* 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++) { + /* 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; + } - *npeaks = max_peaks; + /* Add it to the results if it's unique. */ + if (unique) { + peak_theta[*npeaks] = theta; + peak_phi[*npeaks] = phi; + *npeaks += 1; + } + } free(imax); free(jmax); diff --git a/src/find_peaks.h b/src/find_peaks.h index 8511af6..0219561 100644 --- a/src/find_peaks.h +++ b/src/find_peaks.h @@ -22,6 +22,6 @@ 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); +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); #endif @@ -47,7 +47,7 @@ /* Maximum number of fit parameters. Should be at least 4 + 3*MAX_VERTICES. */ #define MAX_PARS 100 /* Maximum number of peaks to search for in Hough transform. */ -#define MAX_NPEAKS 5 +#define MAX_NPEAKS 10 /* Maximum kinetic energy for any particle. */ #define MAX_ENERGY 10000 @@ -5316,7 +5316,10 @@ 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)); + find_peaks(ev, pos, 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; result = malloc(sizeof(size_t)*n*ipow(npeaks,n)); diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c index a825f4b..dd3a9d8 100644 --- a/src/test-find-peaks.c +++ b/src/test-find-peaks.c @@ -363,7 +363,7 @@ int main(int argc, char **argv) get_hough_transform(&ev,pos,x,y,n,m,result,0,0); - find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,max_npeaks); + find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,max_npeaks,0.1); printf("gtid %i\n", ev.gtid); for (i = 0; i < npeaks; i++) { |