diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/.gitignore | 1 | ||||
| -rw-r--r-- | src/Makefile | 6 | ||||
| -rw-r--r-- | src/fit.c | 85 | ||||
| -rw-r--r-- | src/test-find-peaks.c | 47 | ||||
| -rw-r--r-- | src/zdab-cat.c | 330 | ||||
| -rw-r--r-- | src/zdab_utils.c | 81 | ||||
| -rw-r--r-- | src/zdab_utils.h | 4 |
7 files changed, 423 insertions, 131 deletions
diff --git a/src/.gitignore b/src/.gitignore index f613ee3..47de265 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -9,3 +9,4 @@ test-path release.h calculate-csda-range test-time-pdf +zdab-cat diff --git a/src/Makefile b/src/Makefile index f1b4f0d..bd77dfb 100644 --- a/src/Makefile +++ b/src/Makefile @@ -3,7 +3,7 @@ release_hdr := $(shell sh -c './mkreleasehdr.sh') CFLAGS=-fdiagnostics-color -O2 -Wall -g -DSWAP_BYTES LDLIBS=-fdiagnostics-color -lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++ -all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra test-find-peaks +all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra test-find-peaks zdab-cat Makefile.dep: -$(CC) -MM *.c > Makefile.dep @@ -34,7 +34,9 @@ test-find-peaks: test-find-peaks.o zebra.o likelihood.o pmt.o vector.o misc.o mu calculate-csda-range: calculate-csda-range.o +zdab-cat: zdab-cat.o zebra.o pmt.o vector.o misc.o zdab_utils.o pack2b.o db.o dqxx.o dict.o siphash.o release.o dc.o sort.o util.o + clean: - rm -f *.o calculate_limits test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra test-find-peaks Makefile.dep + rm -f *.o calculate_limits test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra test-find-peaks Makefile.dep zdab-cat .PHONY: all clean @@ -57,9 +57,9 @@ char *GitDirty(void); static int stop = 0; static nlopt_opt opt; -#define EV_RECORD 0x45562020 -#define MCTK_RECORD 0x4d43544b -#define MCVX_RECORD 0x4d435658 +#define EV_RECORD 0x45562020 // 'EV ' +#define MCTK_RECORD 0x4d43544b // 'MCTK' +#define MCVX_RECORD 0x4d435658 // 'MCVX' char *flikelihood; @@ -5812,85 +5812,6 @@ void sigint_handler(int dummy) stop = 1; } -size_t get_nhit(event *ev) -{ - /* Returns the number of PMT hits in event `ev`. - * - * Note: Only hits on normal PMTs which aren't flagged are counted. */ - size_t i, nhit; - - nhit = 0; - 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) continue; - - nhit++; - } - - return nhit; -} - -int get_event(zebraFile *f, event *ev, zebraBank *bev) -{ - /* Read all the PMT banks from the zebra file and update `ev`. - * - * Returns 0 on success, -1 on error. */ - int i, rv; - PMTBank bpmt; - zebraBank b; - int id, crate, card, channel; - - for (i = 0; i < MAX_PMTS; i++) { - ev->pmt_hits[i].hit = 0; - } - - if (bev->links[KEV_PMT-1] == 0) { - /* If the PMT link is zero, we assume it's just a 0 nhit event. */ - return 0; - } - - rv = zebra_get_bank(f,&b,bev->links[KEV_PMT-1]); - - if (rv) { - fprintf(stderr, "error getting PMT bank: %s\n", zebra_err); - return -1; - } - - while (1) { - unpack_pmt(b.data, &bpmt); - card = bpmt.pin/1024; - crate = (bpmt.pin % 1024)/32; - channel = bpmt.pin % 32; - id = crate*512 + card*32 + channel; - ev->pmt_hits[id].hit = 1; - ev->pmt_hits[id].t = bpmt.pt; - ev->pmt_hits[id].qihl = bpmt.pihl; - ev->pmt_hits[id].qihs = bpmt.pihs; - ev->pmt_hits[id].qilx = bpmt.pilx; - ev->pmt_hits[id].qhl = bpmt.phl; - ev->pmt_hits[id].qhs = bpmt.phs; - ev->pmt_hits[id].qlx = bpmt.plx; - /* Clear the PMT_FLAG_DIS bit. */ - ev->pmt_hits[id].flags &= ~PMT_FLAG_DIS; - if (bpmt.pf & KPF_DIS) - ev->pmt_hits[id].flags |= PMT_FLAG_DIS; - - if (!b.next) break; - - rv = zebra_get_bank(f,&b,b.next); - - if (rv) { - fprintf(stderr, "error getting PMT bank: %s\n", zebra_err); - return -1; - } - } - - ev->nhit = get_nhit(ev); - - return 0; -} - void sprintf_particle_string(int *id, size_t n, char *str) { /* Convert a list of particle id codes to a string. diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c index bf623b2..a825f4b 100644 --- a/src/test-find-peaks.c +++ b/src/test-find-peaks.c @@ -171,53 +171,6 @@ void plot_find_peaks(event *ev, double *pos, double *peak_theta, double *peak_ph } } -int get_event(zebraFile *f, event *ev, zebraBank *bev) -{ - /* Read all the PMT banks from the zebra file and update `ev`. - * - * Returns 0 on success, -1 on error. */ - int i, rv; - PMTBank bpmt; - zebraBank b; - int id, crate, card, channel; - - for (i = 0; i < MAX_PMTS; i++) { - ev->pmt_hits[i].hit = 0; - } - - rv = zebra_get_bank(f,&b,bev->links[KEV_PMT-1]); - - if (rv) { - fprintf(stderr, "error getting PMT bank: %s\n", zebra_err); - return -1; - } - - while (1) { - unpack_pmt(b.data, &bpmt); - card = bpmt.pin/1024; - crate = (bpmt.pin % 1024)/32; - channel = bpmt.pin % 32; - id = crate*512 + card*32 + channel; - ev->pmt_hits[id].hit = 1; - ev->pmt_hits[id].t = bpmt.pt; - ev->pmt_hits[id].qhl = bpmt.phl; - ev->pmt_hits[id].qhs = bpmt.phs; - ev->pmt_hits[id].qlx = bpmt.plx; - ev->pmt_hits[id].flags |= bpmt.pf & KPF_DIS; - - if (!b.next) break; - - rv = zebra_get_bank(f,&b,b.next); - - if (rv) { - fprintf(stderr, "error getting PMT bank: %s\n", zebra_err); - return -1; - } - } - - return 0; -} - void usage(void) { fprintf(stderr,"Usage: ./test-find-peaks [options] FILENAME\n"); diff --git a/src/zdab-cat.c b/src/zdab-cat.c new file mode 100644 index 0000000..294e014 --- /dev/null +++ b/src/zdab-cat.c @@ -0,0 +1,330 @@ +/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago> + * + * 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 <https://www.gnu.org/licenses/>. + */ + +#include <stdio.h> +#include "zebra.h" +#include "event.h" +#include "zdab_utils.h" +#include "pmt.h" +#include "db.h" +#include "dqxx.h" +#include <inttypes.h> /* for PRIu32 macro */ +#include <string.h> /* for memcpy() */ +#include <errno.h> /* for errno */ +#include "release.h" + +#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," -h display this help message\n"); + exit(1); +} + +int main(int argc, char **argv) +{ + int i; + zebraFile *f; + zebraBank bmast, bmc, bmcgn, mctk, b; + int rv; + EVBank bev; + MCTKBank bmctk; + MCVXBank bmcvx; + event ev = {0}; + char *filename = NULL; + char *output = NULL; + FILE *fout = NULL; + size_t nhit; + int last_run; + char dqxx_file[256]; + + for (i = 1; i < argc; i++) { + if (argv[i][0] == '-') { + switch (argv[i][1]) { + case 'o': + output = argv[++i]; + 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 (output) { + fout = fopen(output, "w"); + + if (!fout) { + fprintf(stderr, "failed to open '%s': %s\n", output, strerror(errno)); + return 1; + } + + fprintf(fout, "git_sha1: %s\n", GitSHA1()); + fprintf(fout, "git_dirty: %s\n", GitDirty()); + fprintf(fout, "data:\n"); + } + + if (load_pmt_info()) { + zebra_close(f); + if (output) fclose(fout); + 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 (fout) fprintf(fout, " -\n"); + + if (bmast.links[KMAST_MC-1] == 0) goto skip_mc; + + rv = zebra_get_bank(f,&bmc,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; + } + + rv = zebra_get_bank(f,&bmcgn,bmc.links[KMC_MCGN-1]); + + if (rv) { + fprintf(stderr, "error getting MCGN bank: %s\n", zebra_err); + goto err; + } + + if (fout) fprintf(fout, " mcgn:\n"); + 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 (fout) { + fprintf(fout, " -\n"); + fprintf(fout, " id: %" PRIu32 "\n", bmctk.idp); + fprintf(fout, " energy: %.2f\n", bmctk.ene); + fprintf(fout, " posx: %.2f\n", bmcvx.x); + fprintf(fout, " posy: %.2f\n", bmcvx.y); + fprintf(fout, " posz: %.2f\n", bmcvx.z); + fprintf(fout, " dirx: %.4f\n", bmctk.drx); + fprintf(fout, " diry: %.4f\n", bmctk.dry); + fprintf(fout, " dirz: %.4f\n", bmctk.drz); + } + + 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: + + if (bmast.links[KMAST_EV-1] == 0) { + /* First logical record in SNOCR files doesn't have an EV bank. */ + continue; + } + + 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; + } + } + + if (fout) fprintf(fout, " ev:\n"); + while (1) { + unpack_ev(b.data, &bev); + ev.run = bev.run; + ev.gtid = bev.gtr_id; + + if (ev.run != last_run) { + printf("loading DQXX file for run %010i\n", ev.run); + + sprintf(dqxx_file, "DQXX_%010i.dat", ev.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 = ev.run; + } + + rv = get_event(f,&ev,&b); + + nhit = get_nhit(&ev); + + if (fout) { + fprintf(fout, " - gtid: %i\n", ev.gtid); + fprintf(fout, " nhit: %zu\n", nhit); + } + + /* 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) break; + + rv = zebra_get_bank(f,&b,b.orig); + + if (rv) { + fprintf(stderr, "error getting EV bank: %s\n", zebra_err); + goto err; + } + } + } + + db_free(db); + + if (fout) fclose(fout); + + zebra_close(f); + + return 0; + +err: + db_free(db); + + if (fout) fclose(fout); + + zebra_close(f); + + return 1; +} diff --git a/src/zdab_utils.c b/src/zdab_utils.c index dd377b9..37d444c 100644 --- a/src/zdab_utils.c +++ b/src/zdab_utils.c @@ -20,6 +20,87 @@ #include "pack2b.h" #include <stdlib.h> /* for size_t */ #include <stdio.h> /* for fprintf() */ +#include "event.h" +#include "zebra.h" + +size_t get_nhit(event *ev) +{ + /* Returns the number of PMT hits in event `ev`. + * + * Note: Only hits on normal PMTs which aren't flagged are counted. */ + size_t i, nhit; + + nhit = 0; + 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) continue; + + nhit++; + } + + return nhit; +} + +int get_event(zebraFile *f, event *ev, zebraBank *bev) +{ + /* Read all the PMT banks from the zebra file and update `ev`. + * + * Returns 0 on success, -1 on error. */ + int i, rv; + PMTBank bpmt; + zebraBank b; + int id, crate, card, channel; + + for (i = 0; i < MAX_PMTS; i++) { + ev->pmt_hits[i].hit = 0; + } + + if (bev->links[KEV_PMT-1] == 0) { + /* If the PMT link is zero, we assume it's just a 0 nhit event. */ + return 0; + } + + rv = zebra_get_bank(f,&b,bev->links[KEV_PMT-1]); + + if (rv) { + fprintf(stderr, "error getting PMT bank: %s\n", zebra_err); + return -1; + } + + while (1) { + unpack_pmt(b.data, &bpmt); + card = bpmt.pin/1024; + crate = (bpmt.pin % 1024)/32; + channel = bpmt.pin % 32; + id = crate*512 + card*32 + channel; + ev->pmt_hits[id].hit = 1; + ev->pmt_hits[id].t = bpmt.pt; + ev->pmt_hits[id].qihl = bpmt.pihl; + ev->pmt_hits[id].qihs = bpmt.pihs; + ev->pmt_hits[id].qilx = bpmt.pilx; + ev->pmt_hits[id].qhl = bpmt.phl; + ev->pmt_hits[id].qhs = bpmt.phs; + ev->pmt_hits[id].qlx = bpmt.plx; + /* Clear the PMT_FLAG_DIS bit. */ + ev->pmt_hits[id].flags &= ~PMT_FLAG_DIS; + if (bpmt.pf & KPF_DIS) + ev->pmt_hits[id].flags |= PMT_FLAG_DIS; + + if (!b.next) break; + + rv = zebra_get_bank(f,&b,b.next); + + if (rv) { + fprintf(stderr, "error getting PMT bank: %s\n", zebra_err); + return -1; + } + } + + ev->nhit = get_nhit(ev); + + return 0; +} int isOrphan(aPmtEventRecord *pmtRecord) { diff --git a/src/zdab_utils.h b/src/zdab_utils.h index ee4df6e..a2804e1 100644 --- a/src/zdab_utils.h +++ b/src/zdab_utils.h @@ -20,6 +20,8 @@ #include "Record_Info.h" #include <stdint.h> #include <stdlib.h> /* for size_t */ +#include "event.h" +#include "zebra.h" // the builder won't put out events with NHIT > 10000 // (note that these are possible due to hardware problems) @@ -412,6 +414,8 @@ void unpack_mctk(uint32_t *data, MCTKBank *b); void unpack_ev(uint32_t *data, EVBank *b); void unpack_pmt(uint32_t *data, PMTBank *b); +size_t get_nhit(event *ev); +int get_event(zebraFile *f, event *ev, zebraBank *bev); int isOrphan(aPmtEventRecord *pmtRecord); void swap_int32(int32_t *val_pt, int count); void swap_int16(int16_t *val_pt, int count); |