/* Copyright (c) 2019, Anthony Latorre * * This program is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) * any later version. * This program is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for * more details. * You should have received a copy of the GNU General Public License along with * this program. If not, see . */ #include #include "zebra.h" #include "event.h" #include "zdab_utils.h" #include "pmt.h" #include "db.h" #include "dqxx.h" #include /* for PRIu32 macro */ #include /* for memcpy() */ #include /* for errno */ #include "release.h" #include "dc.h" #include "hdf5.h" #include #include "hdf5_utils.h" /* Maximum number of events, primary seed particle tracks, and fits. * * FIXME: I should write out the events, tracks, and fits as they come up * instead of storing them until the end. */ #define MAX_NEVENTS 1000000 #define EV_RECORD 0x45562020 // 'EV ' #define MCTK_RECORD 0x4d43544b // 'MCTK' #define MCVX_RECORD 0x4d435658 // 'MCVX' char *GitSHA1(void); char *GitDirty(void); void usage(void) { fprintf(stderr,"Usage: ./zdab-cat [options] FILENAME\n"); fprintf(stderr," -o output file (default: stdout)\n"); fprintf(stderr," --skip-second-event only fit the first event after a MAST bank\n"); fprintf(stderr," -v verbose\n"); fprintf(stderr," -h display this help message\n"); exit(1); } int main(int argc, char **argv) { int i; zebraFile *f; zebraBank bmast, mc, bmcgn, mctk, b; int rv; EVBank bev; FTPVBank bftpv; FTXKBank bftxk; RSPBank bftxr; MCBank bmc; MCTKBank bmctk; MCVXBank bmcvx; event ev = {0}; char *filename = NULL; char *output = NULL; int skip_second_event = 0; int last_run; char dqxx_file[256]; int nevents = 0; static HDF5Event hdf5_events[MAX_NEVENTS]; int nmcgn = 0; static HDF5MCGN hdf5_mcgn[MAX_NEVENTS]; int verbose = 0; for (i = 1; i < argc; i++) { if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) { if (!strcmp(argv[i]+2,"skip-second-event")) { skip_second_event = 1; continue; } } else if (argv[i][0] == '-') { switch (argv[i][1]) { case 'o': output = argv[++i]; break; case 'v': verbose = 1; break; case 'h': usage(); default: fprintf(stderr, "unrecognized option '%s'\n", argv[i]); exit(1); } } else { filename = argv[i]; } } if (!filename) usage(); f = zebra_open(filename); if (!f) { fprintf(stderr, "%s\n", zebra_err); return 1; } if (load_pmt_info()) { zebra_close(f); return 1; } for (i = 0; i < MAX_PMTS; i++) { ev.pmt_hits[i].hit = 0; ev.pmt_hits[i].flags = 0; } ev.run = -1; dict *db = db_init(); last_run = 10000; if (load_file(db, "DQXX_0000010000.dat", 0)) { fprintf(stderr, "failed to load DQXX_0000010000.dat: %s\n", db_err); goto err; } if (dqxx_init(db, &ev)) { fprintf(stderr, "failed to initialize DQXX bank: %s\n", dqxx_err); goto err; } while (1) { rv = zebra_read_next_logical_record(f); if (rv == -1) { fprintf(stderr, "error getting logical record: %s\n", zebra_err); goto err; } else if (rv == 1) { /* EOF */ break; } rv = zebra_get_bank(f, &bmast, f->first_bank); if (rv) { fprintf(stderr, "error getting MAST bank: %s\n", zebra_err); goto err; } if (bmast.links[KMAST_EV-1] == 0) { /* First logical record in SNOCR files doesn't have an EV bank. */ continue; } if (bmast.links[KMAST_MC-1] == 0) goto skip_mc; rv = zebra_get_bank(f,&mc,bmast.links[KMAST_MC-1]); if (rv) { fprintf(stderr, "error getting MC bank: %s\n", zebra_err); goto err; } if (bmast.links[KMC_MCGN-1] == 0) { fprintf(stderr, "MCGN link is zero!\n"); goto err; } unpack_mc(mc.data, &bmc); rv = zebra_get_bank(f,&bmcgn,mc.links[KMC_MCGN-1]); if (rv) { fprintf(stderr, "error getting MCGN bank: %s\n", zebra_err); goto err; } while (1) { if (bmcgn.links[KMCGN_MCTK-1] == 0) { fprintf(stderr, "MCTK link is zero!\n"); goto err; } rv = zebra_get_bank(f,&mctk,bmcgn.links[KMCGN_MCTK-1]); if (rv) { fprintf(stderr, "error getting MCTK bank: %s\n", zebra_err); goto err; } if (mctk.orig == mctk.up - KMCVX_MCTK) { /* This is the first MCTK bank. */ unpack_mctk(mctk.data, &bmctk); } else { /* For some reason SNOMAN sometimes links to the second MCTK * from the MCGN bank. */ rv = zebra_get_bank(f,&b,mctk.orig); if (b.idh != MCTK_RECORD) { fprintf(stderr, "error following origin link from MCTK bank!\n"); goto err; } if (rv) { fprintf(stderr, "error getting MCTK bank: %s\n", zebra_err); goto err; } unpack_mctk(b.data, &bmctk); } if (mctk.up == 0) { fprintf(stderr, "MCVX link is zero!\n"); goto err; } rv = zebra_get_bank(f,&b,mctk.up); if (rv) { fprintf(stderr, "error getting MCVX bank: %s\n", zebra_err); goto err; } unpack_mcvx(b.data, &bmcvx); if (output) { hdf5_mcgn[nmcgn].evn = bmc.evn; hdf5_mcgn[nmcgn].id = bmctk.idp; hdf5_mcgn[nmcgn].energy = bmctk.ene; hdf5_mcgn[nmcgn].x = bmcvx.x; hdf5_mcgn[nmcgn].y = bmcvx.y; hdf5_mcgn[nmcgn].z = bmcvx.z; hdf5_mcgn[nmcgn].dirx = bmctk.drx; hdf5_mcgn[nmcgn].diry = bmctk.dry; hdf5_mcgn[nmcgn].dirz = bmctk.drz; hdf5_mcgn[nmcgn].time = bmcvx.tim; nmcgn++; } if (bmcgn.next) { rv = zebra_get_bank(f,&bmcgn,bmcgn.next); if (rv) { fprintf(stderr, "error getting MCGN bank: %s\n", zebra_err); goto err; } } else { break; } } skip_mc: rv = zebra_get_bank(f,&b,bmast.links[KMAST_EV-1]); if (rv) { fprintf(stderr, "error getting EV bank: %s\n", zebra_err); goto err; } /* Skip to the last event so we can traverse them in reverse order. The * reason for this is that for some reason SNOMAN puts the events in * reverse order within each logical record. */ while (b.next) { rv = zebra_get_bank(f,&b,b.next); if (rv) { fprintf(stderr, "error getting EV bank: %s\n", zebra_err); goto err; } } while (1) { unpack_ev(b.data, &bev); if (bev.run != last_run) { fprintf(stderr, "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); if (output) { hdf5_events[nevents].run = ev.run; hdf5_events[nevents].evn = bmc.evn; hdf5_events[nevents].gtr = ev.trigger_time; hdf5_events[nevents].nhit = ev.nhit; hdf5_events[nevents].gtid = ev.gtid; hdf5_events[nevents].trg_type = ev.trigger_type; hdf5_events[nevents].dc = get_dc_word(&ev, f, &bmast, &b); if (get_ftpv(f,&b,&bftpv)) { if (verbose) fprintf(stderr, "%s\n", zdab_err); hdf5_events[nevents].ftp_x = NAN; hdf5_events[nevents].ftp_y = NAN; hdf5_events[nevents].ftp_z = NAN; } else { hdf5_events[nevents].ftp_x = bftpv.x; hdf5_events[nevents].ftp_y = bftpv.y; hdf5_events[nevents].ftp_z = bftpv.z; } if (get_ftxk(f,&b,&bftxk)) { if (verbose) fprintf(stderr, "%s\n", zdab_err); hdf5_events[nevents].ftk_energy = NAN; } else { hdf5_events[nevents].ftk_energy = bftxk.energy; } if (get_rsp(f,&b,&bftxr)) { if (verbose) fprintf(stderr, "%s\n", zdab_err); hdf5_events[nevents].rsp_energy = NAN; } else { hdf5_events[nevents].rsp_energy = bftxr.ene; } nevents++; } /* Note the origin link for the first EV bank points back to the * structural link location in the MAST bank. These links are super * confusing! */ if ((b.orig == f->first_bank - KMAST_EV) || skip_second_event) break; rv = zebra_get_bank(f,&b,b.orig); if (rv) { fprintf(stderr, "error getting EV bank: %s\n", zebra_err); goto err; } } } if (output) { save_output(output, hdf5_events, nevents, hdf5_mcgn, nmcgn, NULL, 0); } db_free(db); zebra_close(f); return 0; err: db_free(db); zebra_close(f); return 1; }