diff options
Diffstat (limited to 'src/test-find-peaks.c')
| -rw-r--r-- | src/test-find-peaks.c | 66 |
1 files changed, 61 insertions, 5 deletions
diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c index 11fa4b4..2ccc90d 100644 --- a/src/test-find-peaks.c +++ b/src/test-find-peaks.c @@ -31,12 +31,16 @@ #include "quad.h" #include "likelihood.h" #include "optics.h" +#include "sno_charge.h" +#include "db.h" +#include "pmt_response.h" +#include "dqxx.h" #define EV_RECORD 0x45562020 #define MCTK_RECORD 0x4d43544b #define MCVX_RECORD 0x4d435658 -#define NUM_ANGLES 1000 +#define TEST_FIND_PEAKS_NUM_ANGLES 1000 void plot_hough_transform(double *x, double *y, double *z, size_t n, size_t m) { @@ -101,6 +105,8 @@ void plot_find_peaks(event *ev, double *pos, double *peak_theta, double *peak_ph 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].flags || pmts[i].pmt_type != PMT_NORMAL) continue; + if (ev->pmt_hits[i].hit) { r = NORM(pmts[i].pos); theta = acos(pmts[i].pos[2]/r); @@ -131,9 +137,9 @@ void plot_find_peaks(event *ev, double *pos, double *peak_theta, double *peak_ph /* Rotate the direction by the Cerenkov angle. */ rotate(tmp2,dir,tmp,acos(1/n_d2o)); - for (j = 0; j < NUM_ANGLES; j++) { + for (j = 0; j < TEST_FIND_PEAKS_NUM_ANGLES; j++) { /* Rotate the new direction around the peak by 2 pi. */ - rotate(tmp,tmp2,dir,j*2*M_PI/(NUM_ANGLES-1)); + rotate(tmp,tmp2,dir,j*2*M_PI/(TEST_FIND_PEAKS_NUM_ANGLES-1)); /* Calculate the intersection of the ray with the sphere. */ if (!intersect_sphere(pos,tmp,PSUP_RADIUS,&l)) continue; @@ -197,6 +203,8 @@ int main(int argc, char **argv) int plot = 0; int status; double t0; + int last_run; + char dqxx_file[256]; int32_t gtid = -1; for (i = 1; i < argc; i++) { @@ -233,6 +241,8 @@ int main(int argc, char **argv) load_pmt_info(); + init_charge(); + f = zebra_open(filename); if (!f) { @@ -240,6 +250,30 @@ int main(int argc, char **argv) return 1; } + dict *db = db_init(); + + last_run = -1; + + if (load_file(db, "pmt_response_qoca_d2o_20060216.dat", 0)) { + fprintf(stderr, "failed to load pmt_response_qoca_d2o_20060216.dat: %s\n", db_err); + goto err; + } + + if (load_file(db, "rsp_rayleigh.dat", 0)) { + fprintf(stderr, "failed to load rsp_rayleigh.dat: %s\n", db_err); + goto err; + } + + if (pmt_response_init(db)) { + fprintf(stderr, "failed to initialize PMTR bank: %s\n", pmtr_err); + goto err; + } + + if (optics_init()) { + fprintf(stderr, "failed to initialize optics: %s\n", optics_err); + goto err; + } + while (1) { rv = zebra_read_next_logical_record(f); @@ -251,7 +285,12 @@ int main(int argc, char **argv) break; } - rv = zebra_get_bank(f, &bmast, f->first_bank); + if (f->mast_bank == -1) { + fprintf(stderr, "no MAST bank in logical record! Skipping...\n"); + continue; + } + + rv = zebra_get_bank(f, &bmast, f->mast_bank); if (rv) { fprintf(stderr, "error getting MAST bank: %s\n", zebra_err); @@ -288,6 +327,23 @@ int main(int argc, char **argv) if (gtid > 0 && ev.gtid != gtid) continue; + if (bev.run != last_run) { + printf("loading DQXX file for run %010i\n", bev.run); + + sprintf(dqxx_file, "DQXX_%010i.dat", bev.run); + if (load_file(db, dqxx_file, 1)) { + fprintf(stderr, "failed to load %s: %s\n", dqxx_file, db_err); + goto err; + } + + if (dqxx_init(db, &ev)) { + fprintf(stderr, "failed to initialize DQXX bank: %s\n", dqxx_err); + goto err; + } + + last_run = bev.run; + } + rv = get_event(f,&ev,&b); size_t n = 100; @@ -306,7 +362,7 @@ int main(int argc, char **argv) } /* Guess the position and t0 of the event using the QUAD fitter. */ - status = quad(&ev,pos,&t0,10000); + status = quad(&ev,pos,&t0,10000,0.1); if (status || t0 < 0 || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) { /* If the QUAD fitter fails or returns something outside the PSUP or |