diff options
| -rw-r--r-- | src/Makefile | 2 | ||||
| -rw-r--r-- | src/misc.c | 48 | ||||
| -rw-r--r-- | src/misc.h | 1 | ||||
| -rw-r--r-- | src/test-find-peaks.c | 102 |
4 files changed, 142 insertions, 11 deletions
diff --git a/src/Makefile b/src/Makefile index ca4c39e..f1b4f0d 100644 --- a/src/Makefile +++ b/src/Makefile @@ -30,7 +30,7 @@ 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 find_peaks.o quad.o dc.o sort.o util.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 util.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 util.o quad.o calculate-csda-range: calculate-csda-range.o @@ -256,6 +256,54 @@ double trapz(const double *y, double dx, size_t n) return sum*dx/2.0; } +/* Compute the first intersection of a ray starting at `pos` with direction + * `dir` and a sphere centered at the origin with radius `R`. The distance to + * the intersection is stored in `l`. + * + * Returns 1 if the ray intersects the sphere, and 0 if it doesn't. + * + * Example: + * + * double l; + * double pos[0] = {0,0,0}; + * double dir[3] = {1,0,0}; + * + * if (intersect_sphere(pos,dir,PSUP_RADIUS,&l)) { + * hit[0] = pos[0] + l*dir[0]; + * hit[1] = pos[1] + l*dir[1]; + * hit[2] = pos[2] + l*dir[2]; + * printf("ray intersects sphere at %.2f %.2f %.2f\n", hit[0], hit[1], hit[2]); + * } else { + * printf("ray didn't intersect sphere\n"); + * } + * + */ +int intersect_sphere(double *pos, double *dir, double R, double *l) +{ + double b, c; + + b = 2*DOT(dir,pos); + c = DOT(pos,pos) - R*R; + + if (b*b - 4*c <= 0) { + /* Ray doesn't intersect the sphere. */ + *l = 0.0; + return 0; + } + + /* First, check the shorter solution. */ + *l = (-b - sqrt(b*b - 4*c))/2; + + /* If the shorter solution is less than 0, check the second solution. */ + if (*l < 0) + *l = (-b + sqrt(b*b - 4*c))/2; + + /* If the distance is still negative, we didn't intersect the sphere. */ + if (*l < 0) return 0; + + return 1; +} + void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2) { /* Returns the path length inside and outside a circle of radius `R` for a @@ -28,6 +28,7 @@ #define LNFACT_MAX 100 double trapz(const double *y, double dx, size_t n); +int intersect_sphere(double *pos, double *dir, double R, double *l); void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2); double ln(unsigned int n); double lnfact(unsigned int n); 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); |