diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | src/Makefile | 8 | ||||
| -rw-r--r-- | src/find_peaks.c | 141 | ||||
| -rw-r--r-- | src/find_peaks.h | 11 | ||||
| -rw-r--r-- | src/likelihood.c | 1 | ||||
| -rw-r--r-- | src/test-find-peaks.c | 229 | ||||
| -rw-r--r-- | src/test.c | 69 |
7 files changed, 457 insertions, 3 deletions
@@ -12,3 +12,4 @@ *.nav *.snm *.toc +test-find-peaks diff --git a/src/Makefile b/src/Makefile index 9739428..d8dfdca 100644 --- a/src/Makefile +++ b/src/Makefile @@ -3,7 +3,7 @@ release_hdr := $(shell sh -c './mkreleasehdr.sh') CFLAGS=-O2 -Wall -g -DSWAP_BYTES LDLIBS=-lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++ -all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra +all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra test-find-peaks Makefile.dep: -$(CC) -MM *.c > Makefile.dep @@ -14,7 +14,7 @@ calculate_limits: calculate_limits.c solid_angle.o: solid_angle.c -test: test.o solid_angle.o optics.o muon.o vector.o quantum_efficiency.o pdg.o scattering.o misc.o mt19937ar.o sno_charge.o path.o random.o pmt_response.o db.o dict.o siphash.o electron.o proton.o likelihood.o pmt.o quad.o +test: test.o solid_angle.o optics.o muon.o vector.o quantum_efficiency.o pdg.o scattering.o misc.o mt19937ar.o sno_charge.o path.o random.o pmt_response.o db.o dict.o siphash.o electron.o proton.o likelihood.o pmt.o quad.o find_peaks.o test-vector: test-vector.o vector.o mt19937ar.o @@ -28,7 +28,9 @@ test-time-pdf: test-time-pdf.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o test-zebra: test-zebra.o zebra.o pack2b.o -fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o +fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o + +test-find-peaks: test-find-peaks.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o calculate-csda-range: calculate-csda-range.o diff --git a/src/find_peaks.c b/src/find_peaks.c new file mode 100644 index 0000000..e3a2fc6 --- /dev/null +++ b/src/find_peaks.c @@ -0,0 +1,141 @@ +#include "find_peaks.h" +#include "vector.h" +#include "event.h" +#include "pmt.h" +#include <stddef.h> /* for size_t */ +#include <stdlib.h> /* for malloc() */ +#include "optics.h" +#include <math.h> /* for exp() */ + +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) +{ + /* Find a maximum of `max_peaks` in the 2D array `x` indexed as: + * + * x[i,j] = x[i*m + j] + * + * and store the indices in the two arrays `imax` and `jmax`. + * + * Peaks are defined as any array element which is greater than all 8 of + * it's neighbors. Peaks are also required to be at least max*threshold + * times as high as the highest peak. */ + size_t i, j; + double max; + + if (n*m == 0) { + *npeaks = 0; + return; + } + + /* First, find the highest value in the array (which is guaranteed to be a peak). */ + max = x[0]; + imax[0] = 0; + jmax[0] = 0; + for (i = 0; i < n; i++) { + for (j = 0; j < m; j++) { + if (x[i*m + j] > max) { + max = x[i*m + j]; + imax[0] = i; + jmax[0] = j; + } + } + } + + *npeaks = 1; + + /* Now we look for other peaks which are at least max*threshold high. */ + for (i = 0; i < n; i++) { + for (j = 0; j < m; j++) { + if (x[i*m + j] <= max*threshold) continue; + + if (i == imax[0] && j == jmax[0]) continue; + + /* Check to see if it is actually a peak. */ + if (x[i*m + (j+1) % m] < x[i*m + j] && /* 0 +1 */ + x[i*m + (j+m-1) % m] < x[i*m + j] && /* 0 -1 */ + x[((i+1) % n)*m + j] < x[i*m + j] && /* +1 0 */ + x[((i+1) % n)*m + (j+1) % m] < x[i*m + j] && /* +1 +1 */ + x[((i+1) % n)*m + (j+m-1) % m] < x[i*m + j] && /* +1 -1 */ + x[((i+n-1) % n)*m + j] < x[i*m + j] && /* -1 0 */ + x[((i+n-1) % n)*m + (j+1) % m] < x[i*m + j] && /* -1 +1 */ + x[((i+n-1) % n)*m + (j+m-1) % m] < x[i*m + j]) { /* -1 -1 */ + imax[*npeaks] = i; + jmax[*npeaks] = j; + *npeaks += 1; + + if (*npeaks >= max_peaks) goto end; + } + } + } + +end: + return; +} + +void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n, size_t m, double *result) +{ + /* Computes the "Hough transform" of the event `ev` and stores it in `result`. */ + size_t i, j, k; + double dir[3], pmt_dir[3], cos_theta; + + double wavelength0 = 400.0; + double n_d2o = get_index_snoman_d2o(wavelength0); + + for (i = 0; i < MAX_PMTS; i++) { + if (pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue; + + SUB(pmt_dir,pmts[i].pos,pos); + + normalize(pmt_dir); + + for (j = 0; j < n; j++) { + for (k = 0; k < m; k++) { + double r = 1+x[j]*x[j]+y[k]*y[k]; + dir[0] = 2*x[j]/r; + dir[1] = 2*y[k]/r; + dir[2] = (-1+x[j]*x[j]+y[k]*y[k])/r; + + cos_theta = DOT(pmt_dir,dir); + + result[j*m + k] += ev->pmt_hits[i].qhs*exp(-fabs(cos_theta-1.0/n_d2o)/0.1); + } + } + } +} + +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 threshold) +{ + size_t i; + + double *x = calloc(n,sizeof(double)); + double *y = calloc(m,sizeof(double)); + double *result = calloc(n*m,sizeof(double)); + + for (i = 0; i < n; i++) { + x[i] = -10 + 20.0*i/(n-1); + } + + for (i = 0; i < m; i++) { + y[i] = -10 + 20.0*i/(m-1); + } + + get_hough_transform(ev,pos,x,y,n,m,result); + + size_t *imax = calloc(max_peaks,sizeof(size_t)); + size_t *jmax = calloc(max_peaks,sizeof(size_t)); + + find_peaks_array(result,n,m,imax,jmax,npeaks,max_peaks,threshold); + + for (i = 0; i < *npeaks; i++) { + double R, theta; + R = sqrt(x[imax[i]]*x[imax[i]]+y[jmax[i]]*y[jmax[i]]); + theta = atan2(y[jmax[i]],x[imax[i]]); + peak_theta[i] = 2*atan(1/R); + peak_phi[i] = theta; + } + + free(imax); + free(jmax); + free(x); + free(y); + free(result); +} diff --git a/src/find_peaks.h b/src/find_peaks.h new file mode 100644 index 0000000..2abf6fe --- /dev/null +++ b/src/find_peaks.h @@ -0,0 +1,11 @@ +#ifndef FIND_PEAKS +#define FIND_PEAKS + +#include "event.h" +#include <stddef.h> /* for size_t */ + +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); +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 threshold); + +#endif diff --git a/src/likelihood.c b/src/likelihood.c index 64f28d4..da0614d 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -22,6 +22,7 @@ #include "proton.h" #include "id_particles.h" #include <gsl/gsl_cdf.h> +#include "find_peaks.h" particle *particle_init(int id, double T0, size_t n) { diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c new file mode 100644 index 0000000..c7c2c1d --- /dev/null +++ b/src/test-find-peaks.c @@ -0,0 +1,229 @@ +#include <stdio.h> +#include "find_peaks.h" +#include "misc.h" +#include "zebra.h" +#include "zdab_utils.h" +#include <stddef.h> /* for size_t */ +#include <unistd.h> /* for exit() */ +#include "pmt.h" +#include <math.h> /* for M_PI */ +#include <errno.h> /* for errno */ +#include <string.h> /* for strerror() */ + +#define EV_RECORD 0x45562020 +#define MCTK_RECORD 0x4d43544b +#define MCVX_RECORD 0x4d435658 + +void plot_3d(double *x, double *y, double *z, size_t n, size_t m) +{ + size_t i, j; + + FILE *pipe = popen("gnuplot -p", "w"); + + if (!pipe) { + fprintf(stderr, "error running gnuplot command: %s\n", strerror(errno)); + exit(1); + } + + fprintf(pipe, "set title 'Hough Transform'\n"); + /* Not entirely sure what these do, but following the instructions from + * http://lowrank.net/gnuplot/plot3d2-e.html. */ + fprintf(pipe, "set dgrid3d %zu,%zu\n",n,m); + fprintf(pipe, "set hidden3d\n"); + fprintf(pipe, "set contour\n"); + fprintf(pipe, "splot '-' u 1:2:3 with lines\n"); + + for (i = 0; i < n; i++) { + for (j = 0; j < m; j++) { + fprintf(pipe, "%.10g %.10g %.10g\n", x[i], y[j], z[i*m+j]); + } + } + fprintf(pipe,"e\n"); + + /* Pause so that you can rotate the 3D graph. */ + fprintf(pipe,"pause mouse keypress\n"); + + if (pclose(pipe)) { + fprintf(stderr, "error closing gnuplot command: %s\n", strerror(errno)); + exit(1); + } +} + +int get_event(zebraFile *f, event *ev, bank *b) +{ + /* Read all the PMT banks from the zebra file and update `ev`. + * + * Returns the last return value from get_bank(). */ + int i, rv; + PMTBank bpmt; + int id, crate, card, channel; + + for (i = 0; i < MAX_PMTS; i++) { + ev->pmt_hits[i].hit = 0; + } + + while (1) { + rv = next_bank(f, b); + + if (rv == -1) { + break; + } else if (rv == 1) { + /* EOF */ + break; + } + + switch (b->name) { + case PMT_RECORD: + unpack_pmt(b->data, &bpmt); + card = bpmt.pin/1024; + crate = (bpmt.pin % 1024)/32; + channel = bpmt.pin % 32; + id = crate*512 + card*32 + channel; + ev->pmt_hits[id].hit = 1; + ev->pmt_hits[id].t = bpmt.pt; + ev->pmt_hits[id].qhl = bpmt.phl; + ev->pmt_hits[id].qhs = bpmt.phs; + ev->pmt_hits[id].qlx = bpmt.plx; + ev->pmt_hits[id].flags |= bpmt.pf & KPF_DIS; + break; + case EV_RECORD: + case MAST_RECORD: + goto end; + } + } + +end: + + return rv; +} + +void usage(void) +{ + fprintf(stderr,"Usage: ./test-find-peaks [options] FILENAME\n"); + fprintf(stderr," --plot plot hough transform\n"); + fprintf(stderr," -h display this help message\n"); + exit(1); +} + +int main(int argc, char **argv) +{ + int i; + char *filename = NULL; + zebraFile *f; + bank b; + EVBank bev; + event ev = {0}; + int rv; + MCVXBank bmcvx; + double pos[3]; + double peak_theta[10]; + double peak_phi[10]; + size_t npeaks; + int plot = 0; + + for (i = 1; i < argc; i++) { + if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) { + if (!strcmp(argv[i]+2,"plot")) { + plot = 1; + continue; + } + } else if (argv[i][0] == '-') { + switch (argv[i][1]) { + case 'h': + usage(); + default: + fprintf(stderr, "unrecognized option '%s'\n", argv[i]); + exit(1); + } + } else { + filename = argv[i]; + } + } + + if (!filename) + usage(); + + load_pmt_info(); + + f = zebra_open(filename); + + if (!f) { + fprintf(stderr, "%s\n", zebra_err); + return 1; + } + + int skip = 0; + while (1) { + if (!skip) + rv = next_bank(f, &b); + + skip = 0; + + if (rv == -1) { + fprintf(stderr, "error getting bank: %s\n", zebra_err); + goto err; + } else if (rv == 1) { + /* EOF */ + break; + } + + switch (b.name) { + case MCVX_RECORD: + /* New MC vertex. */ + if (b.status & KMCVX_CRE) { + unpack_mcvx(b.data, &bmcvx); + } + break; + case EV_RECORD: + unpack_ev(b.data, &bev); + ev.run = bev.run; + ev.gtid = bev.gtr_id; + + rv = get_event(f,&ev,&b); + + pos[0] = bmcvx.x; + pos[1] = bmcvx.y; + pos[2] = bmcvx.z; + + size_t n = 100; + size_t m = 100; + + double *x = calloc(n,sizeof(double)); + double *y = calloc(m,sizeof(double)); + double *result = calloc(n*m,sizeof(double)); + + for (i = 0; i < n; i++) { + x[i] = -10 + 20.0*i/(n-1); + } + + for (i = 0; i < m; i++) { + y[i] = -10 + 20.0*i/(m-1); + } + + get_hough_transform(&ev,pos,x,y,n,m,result); + + find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,LEN(peak_theta),0.5); + + printf("gtid %i\n", ev.gtid); + for (i = 0; i < npeaks; i++) { + printf("%.2f %.2f\n", peak_theta[i], peak_phi[i]); + } + + if (plot) + plot_3d(x,y,result,n,m); + + /* Skip reading in the next bank since get_event() already read in + * the next bank. */ + skip = 1; + } + } + + zebra_close(f); + + return 0; + +err: + zebra_close(f); + + return 1; +} @@ -24,6 +24,7 @@ #include <gsl/gsl_cdf.h> #include "sno.h" #include "quad.h" +#include "find_peaks.h" typedef int testFunction(char *err); @@ -1566,6 +1567,73 @@ err: return 1; } +int test_find_peaks_array(char *err) +{ + /* Tests the find_peaks_array() function. */ + size_t i; + double x[4] = {0,1,0,0}; + size_t imax[10], jmax[10], npeaks; + + find_peaks_array(x,2,2,imax,jmax,&npeaks,LEN(imax),0.1); + + if (npeaks != 1) { + sprintf(err, "number of peaks = %zu, but expected %i", npeaks, 1); + return 1; + } + + if (imax[0] != 0) { + sprintf(err, "imax[0] = %zu, but expected %i", imax[0], 0); + return 1; + } + + if (jmax[0] != 1) { + sprintf(err, "jmax[0] = %zu, but expected %i", jmax[0], 1); + return 1; + } + + double y[10][10] = { + {0,0,2,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0}, + {4,0,0,0,0,0,0,0,0,5}, + {0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0}, + {0,0,3,0,0,0,0,0,0,0}, + }; + + find_peaks_array(y,10,10,imax,jmax,&npeaks,LEN(imax),0.1); + + if (npeaks != 2) { + sprintf(err, "number of peaks = %zu, but expected %i", npeaks, 2); + return 1; + } + + if (imax[0] != 2) { + sprintf(err, "imax[0] = %zu, but expected %i", imax[0], 2); + return 1; + } + + if (jmax[0] != 9) { + sprintf(err, "jmax[0] = %zu, but expected %i", jmax[0], 9); + return 1; + } + + if (imax[1] != 9) { + sprintf(err, "imax[1] = %zu, but expected %i", imax[1], 9); + return 1; + } + + if (jmax[1] != 2) { + sprintf(err, "jmax[1] = %zu, but expected %i", jmax[1], 2); + return 1; + } + + return 0; +} + struct tests { testFunction *test; char *name; @@ -1606,6 +1674,7 @@ struct tests { {test_time_cdf, "test_time_cdf"}, {test_quad, "test_quad"}, {test_quad_noise, "test_quad_noise"}, + {test_find_peaks_array, "test_find_peaks_array"}, }; int main(int argc, char **argv) |