update test-find-peaks to plot cerenkov rings

This commit updates the test-find-peaks script to plot Cerenkov rings for each
of the peaks. It also updates the script to use quad to find the position
instead of using the MC information. Finally, I added a -n argument to the
script to specify how many peaks to draw.

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);