add zdab-cat

This commit adds a new program called zdab-cat which is kind of like fit, but
just produces the YAML output without actually fitting anything.

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
diff --git a/src/fit.c b/src/fit.c
index 2c8c544..2e2fee3 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -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);