diff options
Diffstat (limited to 'src/test-find-peaks.c')
| -rw-r--r-- | src/test-find-peaks.c | 102 |
1 files changed, 92 insertions, 10 deletions
diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c index e14a794..bf623b2 100644 --- a/src/test-find-peaks.c +++ b/src/test-find-peaks.c @@ -27,11 +27,17 @@ #include <string.h> /* for strerror() */ #include "vector.h" #include "util.h" +#include "sno.h" +#include "quad.h" +#include "likelihood.h" +#include "optics.h" #define EV_RECORD 0x45562020 #define MCTK_RECORD 0x4d43544b #define MCVX_RECORD 0x4d435658 +#define NUM_ANGLES 1000 + void plot_hough_transform(double *x, double *y, double *z, size_t n, size_t m) { size_t i, j; @@ -67,9 +73,16 @@ void plot_hough_transform(double *x, double *y, double *z, size_t n, size_t m) } } -void plot_find_peaks(event *ev, double *peak_theta, double *peak_phi, size_t n) +void plot_find_peaks(event *ev, double *pos, double *peak_theta, double *peak_phi, size_t n) { - size_t i; + size_t i, j; + double r, theta, phi, last_phi; + double dir[3], k[3], tmp[3], tmp2[3], hit[3]; + double l; + double wavelength0, n_d2o; + + wavelength0 = 400.0; + n_d2o = get_index_snoman_d2o(wavelength0); FILE *pipe = popen("gnuplot -p", "w"); @@ -85,21 +98,67 @@ void plot_find_peaks(event *ev, double *peak_theta, double *peak_phi, size_t n) * http://lowrank.net/gnuplot/plot3d2-e.html. */ fprintf(pipe, "set xrange [0:6.28]\n"); fprintf(pipe, "set yrange [3.14:0]\n"); - fprintf(pipe, "plot '-' u 1:2:3:4 with circles palette fillstyle solid, '-' u 1:2 with circles lc rgb \"red\" lw 2\n"); + fprintf(pipe, "plot '-' u 1:2:3:4 with circles palette fillstyle solid, '-' u 1:2 with lines lc rgb \"red\" lw 2\n"); for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].hit) { - double r = NORM(pmts[i].pos); - double theta = acos(pmts[i].pos[2]/r); - double phi = atan2(pmts[i].pos[1],pmts[i].pos[0]); + r = NORM(pmts[i].pos); + theta = acos(pmts[i].pos[2]/r); + phi = atan2(pmts[i].pos[1],pmts[i].pos[0]); phi = (phi < 0) ? phi + 2*M_PI: phi; fprintf(pipe, "%.10g %.10g %.10g %.10g\n", phi, theta, 0.01, ev->pmt_hits[i].qhs); } } fprintf(pipe,"e\n"); + /* Now we draw the Cerenkov rings for each peak. We do this by computing + * rays which lie along the Cerenkov cone, and intersecting each of these + * with the PSUP. */ for (i = 0; i < n; i++) { - fprintf(pipe, "%.10g %.10g\n", peak_phi[i], peak_theta[i]); + /* Compute the direction of the peak. */ + dir[0] = sin(peak_theta[i])*cos(peak_phi[i]); + dir[1] = sin(peak_theta[i])*sin(peak_phi[i]); + dir[2] = cos(peak_theta[i]); + + k[0] = 0; + k[1] = 0; + k[2] = 1; + + CROSS(tmp,dir,k); + + normalize(tmp); + + /* Rotate the direction by the Cerenkov angle. */ + rotate(tmp2,dir,tmp,acos(1/n_d2o)); + + for (j = 0; j < NUM_ANGLES; j++) { + /* Rotate the new direction around the peak by 2 pi. */ + rotate(tmp,tmp2,dir,j*2*M_PI/(NUM_ANGLES-1)); + + /* Calculate the intersection of the ray with the sphere. */ + if (!intersect_sphere(pos,tmp,PSUP_RADIUS,&l)) continue; + + hit[0] = pos[0] + l*tmp[0]; + hit[1] = pos[1] + l*tmp[1]; + hit[2] = pos[2] + l*tmp[2]; + + /* Calculate the theta and phi angles for the intersection. */ + theta = acos(hit[2]/NORM(hit)); + phi = atan2(hit[1],hit[0]); + + /* Since the phi angle is periodic, avoid drawing a horizontal line + * when we switch from the left side of the plot to the right side + * or vice versa. */ + if ((j > 0) && ((last_phi < 0 && phi > 0) || (last_phi > 0 && phi < 0))) + fprintf(pipe,"\n"); + + last_phi = phi; + + phi = (phi < 0) ? phi + 2*M_PI: phi; + + fprintf(pipe, "%.10g %.10g\n", phi, theta); + } + fprintf(pipe,"\n"); } fprintf(pipe,"e\n"); @@ -162,6 +221,7 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev) void usage(void) { fprintf(stderr,"Usage: ./test-find-peaks [options] FILENAME\n"); + fprintf(stderr," -n number of peaks to plot\n"); fprintf(stderr," --plot plot hough transform\n"); fprintf(stderr," -h display this help message\n"); exit(1); @@ -180,8 +240,10 @@ int main(int argc, char **argv) double pos[3]; double peak_theta[10]; double peak_phi[10]; - size_t npeaks; + size_t npeaks, max_npeaks = 10; int plot = 0; + int status; + double t0; for (i = 1; i < argc; i++) { if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) { @@ -191,6 +253,13 @@ int main(int argc, char **argv) } } else if (argv[i][0] == '-') { switch (argv[i][1]) { + case 'n': + max_npeaks = atoi(argv[++i]); + if (max_npeaks > LEN(peak_theta)) { + fprintf(stderr, "number of peaks must be less than %lu\n", LEN(peak_theta)); + exit(1); + } + break; case 'h': usage(); default: @@ -326,9 +395,22 @@ int main(int argc, char **argv) y[i] = i*2*M_PI/(m-1); } + /* Guess the position and t0 of the event using the QUAD fitter. */ + status = quad(&ev,pos,&t0,10000); + + if (status || t0 < 0 || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) { + /* If the QUAD fitter fails or returns something outside the PSUP or + * with an invalid time we just assume it's at the center. */ + fprintf(stderr, "quad returned pos = %.2f, %.2f, %.2f t0 = %.2f. Assuming vertex is at the center!\n", + pos[0], pos[1], pos[2], t0); + pos[0] = 0.0; + pos[1] = 0.0; + pos[2] = 0.0; + } + get_hough_transform(&ev,pos,x,y,n,m,result,0,0); - find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,LEN(peak_theta)); + find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,max_npeaks); printf("gtid %i\n", ev.gtid); for (i = 0; i < npeaks; i++) { @@ -336,7 +418,7 @@ int main(int argc, char **argv) } if (plot) - plot_find_peaks(&ev,peak_theta,peak_phi,npeaks); + plot_find_peaks(&ev,pos,peak_theta,peak_phi,npeaks); } zebra_close(f); |