1 files changed, 52 insertions, 6 deletions

diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c
index c7c2c1d..c82f9ec 100644
--- a/src/test-find-peaks.c
+++ b/src/test-find-peaks.c
@@ -9,12 +9,13 @@
 #include <math.h> /* for M_PI */
 #include <errno.h> /* for errno */
 #include <string.h> /* for strerror() */
+#include "vector.h"
 
 #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)
+void plot_hough_transform(double *x, double *y, double *z, size_t n, size_t m)
 {
     size_t i, j;
 
@@ -49,6 +50,51 @@ void plot_3d(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)
+{
+    size_t i;
+
+    FILE *pipe = popen("gnuplot -p", "w");
+
+    if (!pipe) {
+        fprintf(stderr, "error running gnuplot command: %s\n", strerror(errno));
+        exit(1);
+    }
+
+    fprintf(pipe, "set macros\n");
+    fprintf(pipe, "load 'viridis.pal'\n");
+    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 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");
+
+    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]);
+            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");
+
+    for (i = 0; i < n; i++) {
+        fprintf(pipe, "%.10g %.10g\n", peak_phi[i], peak_theta[i]);
+    }
+    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`.
@@ -193,16 +239,16 @@ int main(int argc, char **argv)
             double *result = calloc(n*m,sizeof(double));
 
             for (i = 0; i < n; i++) {
-                x[i] = -10 + 20.0*i/(n-1);
+                x[i] = i*M_PI/(n-1);
             }
 
             for (i = 0; i < m; i++) {
-                y[i] = -10 + 20.0*i/(m-1);
+                y[i] = i*2*M_PI/(m-1);
             }
 
-            get_hough_transform(&ev,pos,x,y,n,m,result);
+            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),0.5);
+            find_peaks(&ev,pos,n,m,peak_theta,peak_phi,&npeaks,LEN(peak_theta));
 
             printf("gtid %i\n", ev.gtid);
             for (i = 0; i < npeaks; i++) {
@@ -210,7 +256,7 @@ int main(int argc, char **argv)
             }
 
             if (plot)
-                plot_3d(x,y,result,n,m);
+                plot_find_peaks(&ev,peak_theta,peak_phi,npeaks);
 
             /* Skip reading in the next bank since get_event() already read in
              * the next bank. */