switch from YAML output to HDF5 to speed things up

7 files changed, 497 insertions, 149 deletions

diff --git a/src/Makefile b/src/Makefile
index 99cff5e..e5ca1cd 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,7 +1,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++ -lz
+LDLIBS=-fdiagnostics-color -lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++ -lz -lhdf5
 
 PREFIX?=$(HOME)/local
 INSTALL_BIN=$(PREFIX)/bin
@@ -32,17 +32,18 @@ test-time-pdf: test-time-pdf.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o
 
 test-zebra: test-zebra.o zebra.o pack2b.o
 
-fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o quad.o dc.o sort.o util.o
+fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o quad.o dc.o sort.o util.o hdf5_utils.o
 
 test-find-peaks: test-find-peaks.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o util.o quad.o
 
 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
+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 hdf5_utils.o
 
 install:
 	@mkdir -p $(INSTALL_BIN)
 	$(INSTALL) fit $(INSTALL_BIN)
+	$(INSTALL) zdab-cat $(INSTALL_BIN)
 
 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 zdab-cat
diff --git a/src/fit.c b/src/fit.c
index adc6ebc..0360a8f 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -47,12 +47,21 @@
 #include "electron.h"
 #include "muon.h"
 #include "proton.h"
+#include "hdf5.h"
+#include "hdf5_utils.h"
 
 /* Maximum number of fit parameters. Should be at least 4 + 3*MAX_VERTICES. */
 #define MAX_PARS 100
+
 /* Maximum number of peaks to search for in Hough transform. */
 #define MAX_NPEAKS 10
 
+/* 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
+
 char *GitSHA1(void);
 char *GitDirty(void);
 
@@ -5869,12 +5878,13 @@ int main(int argc, char **argv)
 {
     int i, j, k;
     zebraFile *f;
-    zebraBank bmast, bmc, bmcgn, mctk, b;
+    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};
@@ -5882,7 +5892,6 @@ int main(int argc, char **argv)
     double xopt[MAX_PARS];
     char *filename = NULL;
     char *output = NULL;
-    FILE *fout = NULL;
     int skip_second_event = 0;
     struct timeval tv_start, tv_stop;
     long long elapsed;
@@ -5892,12 +5901,29 @@ int main(int argc, char **argv)
     int id[MAX_VERTICES];
     size_t combos[100];
     size_t len;
-    char tmp[256];
-    size_t nhit, min_nhit = 100;
+    size_t min_nhit = 100;
     int last_run;
     char dqxx_file[256];
     int32_t gtid = -1;
 
+    int nevents = 0;
+    /* Array of events to write out to HDF5 file.
+     *
+     * Note: Declared static since otherwise it will cause a stack overflow. */
+    static HDF5Event hdf5_events[MAX_NEVENTS];
+
+    int nmcgn = 0;
+    /* Array of primary seed tracks to write out to HDF5 file.
+     *
+     * Note: Declared static since otherwise it will cause a stack overflow. */
+    static HDF5MCGN hdf5_mcgn[MAX_NEVENTS];
+
+    int nfits = 0;
+    /* Array of fit results to write out to HDF5 file.
+     *
+     * Note: Declared static since otherwise it will cause a stack overflow. */
+    static HDF5Fit hdf5_fits[MAX_NEVENTS];
+
     for (i = 1; i < argc; i++) {
         if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) {
             if (!strcmp(argv[i]+2,"skip-second-event")) {
@@ -5953,22 +5979,8 @@ int main(int argc, char **argv)
         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;
     }
 
@@ -6038,30 +6050,29 @@ int main(int argc, char **argv)
             continue;
         }
 
-        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]);
+        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;
         }
 
+        unpack_mc(mc.data, &bmc);
+
         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]);
+        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;
         }
 
-        if (fout) fprintf(fout, "    mcgn:\n");
         while (1) {
             if (bmcgn.links[KMCGN_MCTK-1] == 0) {
                 fprintf(stderr, "MCTK link is zero!\n");
@@ -6110,17 +6121,18 @@ int main(int argc, char **argv)
 
             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);
-                fprintf(fout, "        time: %.4f\n", bmcvx.tim);
+            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) {
@@ -6156,7 +6168,6 @@ skip_mc:
             }
         }
 
-        if (fout) fprintf(fout, "    ev:\n");
         while (1) {
             unpack_ev(b.data, &bev);
             ev.run = bev.run;
@@ -6185,49 +6196,44 @@ skip_mc:
 
             if (get_event(f,&ev,&b)) goto skip_event;
 
-            nhit = get_nhit(&ev);
+            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 (fout) {
-                fprintf(fout, "      - run: %i\n", ev.run);
-                fprintf(fout, "        gtr: %.0f\n", ev.trigger_time);
-                fprintf(fout, "        nhit: %zu\n", nhit);
-                fprintf(fout, "        gtid: %i\n", ev.gtid);
-                fprintf(fout, "        trg_type: 0x%08x\n", ev.trigger_type);
-                fprintf(fout, "        dc: 0x%08x\n", get_dc_word(&ev, f, &bmast, &b));
-            }
-
-            if (fout) {
                 if (get_ftpv(f,&b,&bftpv)) {
                     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 {
-                    fprintf(fout, "        ftp:\n");
-                    fprintf(fout, "            x: %.2f\n", bftpv.x);
-                    fprintf(fout, "            y: %.2f\n", bftpv.y);
-                    fprintf(fout, "            z: %.2f\n", bftpv.z);
+                    hdf5_events[nevents].ftp_x = bftpv.x;
+                    hdf5_events[nevents].ftp_y = bftpv.y;
+                    hdf5_events[nevents].ftp_z = bftpv.z;
                 }
-            }
 
-            if (fout) {
                 if (get_ftxk(f,&b,&bftxk)) {
                     fprintf(stderr, "%s\n", zdab_err);
+                    hdf5_events[nevents].ftk_energy = NAN;
                 } else {
-                    fprintf(fout, "        ftk:\n");
-                    fprintf(fout, "            energy: %.2f\n", bftxk.energy);
+                    hdf5_events[nevents].ftk_energy = bftxk.energy;
                 }
-            }
 
-            if (fout) {
                 if (get_rsp(f,&b,&bftxr)) {
                     fprintf(stderr, "%s\n", zdab_err);
+                    hdf5_events[nevents].rsp_energy = NAN;
                 } else {
-                    fprintf(fout, "        rsp:\n");
-                    fprintf(fout, "            energy: %.2f\n", bftxr.ene);
+                    hdf5_events[nevents].rsp_energy = bftxr.ene;
                 }
-            }
 
-            if (nhit < min_nhit) goto skip_event;
+                nevents++;
+            }
 
-            if (fout) fprintf(fout, "        fit:\n");
+            if (ev.nhit < min_nhit) goto skip_event;
 
             fmin_best = nll_best(&ev);
 
@@ -6251,25 +6257,47 @@ skip_mc:
 
                     elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
 
-                    if (fout) {
-                        sprintf_particle_string(id,i,tmp);
-                        fprintf(fout, "          %s:\n", tmp);
-                        fprintf(fout, "            posx: %.2f\n", xopt[0]);
-                        fprintf(fout, "            posy: %.2f\n", xopt[1]);
-                        fprintf(fout, "            posz: %.2f\n", xopt[2]);
-                        fprintf(fout, "            t0: %.2f\n", xopt[3]);
-                        sprintf_yaml_list(xopt+4,i,3,tmp);
-                        fprintf(fout, "            energy: %s\n", tmp);
-                        sprintf_yaml_list(xopt+5,i,3,tmp);
-                        fprintf(fout, "            theta: %s\n", tmp);
-                        sprintf_yaml_list(xopt+6,i,3,tmp);
-                        fprintf(fout, "            phi: %s\n", tmp);
-                        //fprintf(fout, "            z1: %.2f\n", xopt[7]);
-                        //fprintf(fout, "            z2: %.2f\n", xopt[8]);
-                        fprintf(fout, "            fmin: %.2f\n", fmin);
-                        fprintf(fout, "            time: %lld\n", elapsed);
-                        fprintf(fout, "            psi: %.2f\n", fmin-fmin_best);
-                        fflush(fout);
+                    if (output) {
+                        hdf5_fits[nfits].run = ev.run;
+                        hdf5_fits[nfits].gtid = ev.gtid;
+                        hdf5_fits[nfits].n = i;
+                        hdf5_fits[nfits].x = xopt[0];
+                        hdf5_fits[nfits].y = xopt[1];
+                        hdf5_fits[nfits].z = xopt[2];
+                        hdf5_fits[nfits].t0 = xopt[3];
+                        hdf5_fits[nfits].id1 = id[0];
+                        hdf5_fits[nfits].energy1 = xopt[4];
+                        hdf5_fits[nfits].theta1 = xopt[5];
+                        hdf5_fits[nfits].phi1 = xopt[6];
+
+                        if (i > 1) {
+                            hdf5_fits[nfits].id2 = id[1];
+                            hdf5_fits[nfits].energy2 = xopt[7];
+                            hdf5_fits[nfits].theta2 = xopt[8];
+                            hdf5_fits[nfits].phi2 = xopt[9];
+                        } else {
+                            hdf5_fits[nfits].id2 = 0;
+                            hdf5_fits[nfits].energy2 = NAN;
+                            hdf5_fits[nfits].theta2 = NAN;
+                            hdf5_fits[nfits].phi2 = NAN;
+                        }
+
+                        if (i > 2) {
+                            hdf5_fits[nfits].id3 = id[2];
+                            hdf5_fits[nfits].energy3 = xopt[10];
+                            hdf5_fits[nfits].theta3 = xopt[11];
+                            hdf5_fits[nfits].phi3 = xopt[12];
+                        } else {
+                            hdf5_fits[nfits].id3 = 0;
+                            hdf5_fits[nfits].energy3 = NAN;
+                            hdf5_fits[nfits].theta3 = NAN;
+                            hdf5_fits[nfits].phi3 = NAN;
+                        }
+
+                        hdf5_fits[nfits].fmin = fmin;
+                        hdf5_fits[nfits].time = elapsed;
+                        hdf5_fits[nfits].psi = fmin-fmin_best;
+                        nfits++;
                     }
                 }
             }
@@ -6292,13 +6320,13 @@ skip_event:
 
 end:
 
+    if (output) save_output(output, hdf5_events, nevents, hdf5_mcgn, nmcgn, hdf5_fits, nfits);
+
     free_interpolation();
     pmt_response_free();
 
     db_free(db);
 
-    if (fout) fclose(fout);
-
     zebra_close(f);
 
     return 0;
@@ -6309,8 +6337,6 @@ err:
 
     db_free(db);
 
-    if (fout) fclose(fout);
-
     zebra_close(f);
 
     return 1;
diff --git a/src/hdf5_utils.c b/src/hdf5_utils.c
new file mode 100644
index 0000000..05d1002
--- /dev/null
+++ b/src/hdf5_utils.c
@@ -0,0 +1,193 @@
+#include "hdf5_utils.h"
+#include "hdf5.h"
+#include <string.h> /* for strlen() */
+
+char *GitSHA1(void);
+char *GitDirty(void);
+
+/* Write out data to an HDF5 file.
+ *
+ * This function writes out three datasets to an HDF5 file:
+ *
+ * 1. ev: the event dataset
+ *
+ *    This dataset contains event level information similar to the EV bank in
+ *    SNOMAN. It also contains a data cleaning word and the FTP, FTK, and RSP
+ *    fit values for the position and energy.
+ *
+ * 2. mcgn: the primary particle tracks
+ *
+ *     This dataset contains vertex and track information for all the primary
+ *     particle tracks used to simulate the event.
+ *
+ * 3. fits: the fit information
+ *
+ *     This dataset contains the reconstructed fit information for each event.
+ *     Note that this dataset will in general contain multiple rows for a
+ *     single event since we fit for different particle hypotheses.
+ *
+ * Currently the ev and mcgn datasets are fixed size, and the fits dataset is
+ * chunked and extensible to allow the cat-grid-jobs script to work.
+ *
+ */
+int save_output(const char *output, HDF5Event *hdf5_events, size_t nevents, HDF5MCGN *hdf5_mcgn, size_t nmcgn, HDF5Fit *hdf5_fits, size_t nfits)
+{
+    hid_t      hdf5_events_tid;     /* File datatype identifier */
+    hid_t      file, dataset, space; /* Handles */
+    herr_t     status;
+    hsize_t    dim[1];   /* Dataspace dimensions */
+    hid_t      hdf5_mcgn_tid;     /* File datatype identifier */
+    hid_t      dataset_mcgn, space_mcgn; /* Handles */
+    hsize_t    dim_mcgn[1];   /* Dataspace dimensions */
+    hid_t      hdf5_fits_tid;     /* File datatype identifier */
+    hid_t      dataset_fits, space_fits; /* Handles */
+    hsize_t    dim_fits[1];   /* Dataspace dimensions */
+    hsize_t    chunk_dims[1];   /* Dataspace dimensions */
+    hid_t      dcpl_id; /* Handles */
+    hsize_t    maxdims[1] = {H5S_UNLIMITED};
+
+    /*
+     * Create the data space.
+     */
+    dim[0] = nevents;
+    space = H5Screate_simple(1, dim, NULL);
+    dim_mcgn[0] = nmcgn;
+    space_mcgn = H5Screate_simple(1, dim_mcgn, NULL);
+    dim_fits[0] = nfits;
+    space_fits = H5Screate_simple(1, dim_fits, maxdims);
+    chunk_dims[0] = 10;
+    dcpl_id = H5Pcreate(H5P_DATASET_CREATE);
+    H5Pset_chunk(dcpl_id, 1, chunk_dims);
+
+    /*
+     * Create the file.
+     */
+    file = H5Fcreate(output, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+
+    hid_t   attr1, attr2; /* Attribute identifiers */
+    hid_t   aid1, aid2;    /* Attribute dataspace identifiers */
+    hid_t   atype1, atype2;    /* Attribute type */
+
+    char *git_sha1 = GitSHA1();
+    char *git_dirty = GitDirty();
+
+    /*
+     * Create string attribute.
+     */
+    aid1  = H5Screate(H5S_SCALAR);
+    atype1 = H5Tcopy(H5T_C_S1);
+             H5Tset_size(atype1, strlen(git_sha1));
+             H5Tset_strpad(atype1,H5T_STR_NULLTERM);
+    attr1 = H5Acreate2(file, "git_sha1", atype1, aid1, H5P_DEFAULT, H5P_DEFAULT);
+
+    /*
+     * Create string attribute.
+     */
+    aid2  = H5Screate(H5S_SCALAR);
+    atype2 = H5Tcopy(H5T_C_S1);
+             H5Tset_size(atype2, strlen(git_dirty));
+             H5Tset_strpad(atype2,H5T_STR_NULLTERM);
+    attr2 = H5Acreate2(file, "git_dirty", atype2, aid2, H5P_DEFAULT, H5P_DEFAULT);
+
+    herr_t  ret;                /* Return value */
+
+    /*
+     * Write string attribute.
+     */
+    ret = H5Awrite(attr1, atype1, git_sha1);
+    ret = H5Awrite(attr2, atype2, git_dirty);
+
+    /*
+     * Create the memory data type.
+     */
+    hdf5_events_tid = H5Tcreate(H5T_COMPOUND, sizeof(HDF5Event));
+    H5Tinsert(hdf5_events_tid, "run",        HOFFSET(HDF5Event, run),        H5T_NATIVE_INT);
+    H5Tinsert(hdf5_events_tid, "evn",        HOFFSET(HDF5Event, evn),        H5T_NATIVE_INT);
+    H5Tinsert(hdf5_events_tid, "gtr",        HOFFSET(HDF5Event, gtr),        H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_events_tid, "nhit",       HOFFSET(HDF5Event, nhit),       H5T_NATIVE_INT);
+    H5Tinsert(hdf5_events_tid, "gtid",       HOFFSET(HDF5Event, gtid),       H5T_NATIVE_UINT32);
+    H5Tinsert(hdf5_events_tid, "trg_type",   HOFFSET(HDF5Event, trg_type),   H5T_NATIVE_UINT32);
+    H5Tinsert(hdf5_events_tid, "dc",         HOFFSET(HDF5Event, dc),         H5T_NATIVE_UINT32);
+    H5Tinsert(hdf5_events_tid, "ftp_x",      HOFFSET(HDF5Event, ftp_x),      H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_events_tid, "ftp_y",      HOFFSET(HDF5Event, ftp_y),      H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_events_tid, "ftp_z",      HOFFSET(HDF5Event, ftp_z),      H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_events_tid, "ftk_energy", HOFFSET(HDF5Event, ftk_energy), H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_events_tid, "rsp_energy", HOFFSET(HDF5Event, rsp_energy), H5T_NATIVE_DOUBLE);
+
+    hdf5_mcgn_tid = H5Tcreate(H5T_COMPOUND, sizeof(HDF5MCGN));
+    H5Tinsert(hdf5_mcgn_tid, "run",    HOFFSET(HDF5MCGN, run),    H5T_NATIVE_INT);
+    H5Tinsert(hdf5_mcgn_tid, "evn",    HOFFSET(HDF5MCGN, evn),    H5T_NATIVE_INT);
+    H5Tinsert(hdf5_mcgn_tid, "id",     HOFFSET(HDF5MCGN, id),     H5T_NATIVE_INT);
+    H5Tinsert(hdf5_mcgn_tid, "energy", HOFFSET(HDF5MCGN, energy), H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_mcgn_tid, "x",      HOFFSET(HDF5MCGN, x),      H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_mcgn_tid, "y",      HOFFSET(HDF5MCGN, y),      H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_mcgn_tid, "z",      HOFFSET(HDF5MCGN, z),      H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_mcgn_tid, "dirx",   HOFFSET(HDF5MCGN, dirx),   H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_mcgn_tid, "diry",   HOFFSET(HDF5MCGN, diry),   H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_mcgn_tid, "dirz",   HOFFSET(HDF5MCGN, dirz),   H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_mcgn_tid, "time",   HOFFSET(HDF5MCGN, time),   H5T_NATIVE_DOUBLE);
+
+    hdf5_fits_tid = H5Tcreate(H5T_COMPOUND, sizeof(HDF5Fit));
+    H5Tinsert(hdf5_fits_tid, "run",     HOFFSET(HDF5Fit, run),     H5T_NATIVE_INT);
+    H5Tinsert(hdf5_fits_tid, "gtid",    HOFFSET(HDF5Fit, gtid),    H5T_NATIVE_UINT32);
+    H5Tinsert(hdf5_fits_tid, "n",       HOFFSET(HDF5Fit, n),       H5T_NATIVE_INT);
+    H5Tinsert(hdf5_fits_tid, "x",       HOFFSET(HDF5Fit, x),       H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "y",       HOFFSET(HDF5Fit, y),       H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "z",       HOFFSET(HDF5Fit, z),       H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "t0",      HOFFSET(HDF5Fit, t0),      H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "id1",     HOFFSET(HDF5Fit, id1),     H5T_NATIVE_INT);
+    H5Tinsert(hdf5_fits_tid, "energy1", HOFFSET(HDF5Fit, energy1), H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "theta1",  HOFFSET(HDF5Fit, theta1),  H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "phi1",    HOFFSET(HDF5Fit, phi1),    H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "id2",     HOFFSET(HDF5Fit, id2),     H5T_NATIVE_INT);
+    H5Tinsert(hdf5_fits_tid, "energy2", HOFFSET(HDF5Fit, energy2), H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "theta2",  HOFFSET(HDF5Fit, theta2),  H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "phi2",    HOFFSET(HDF5Fit, phi2),    H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "id3",     HOFFSET(HDF5Fit, id3),     H5T_NATIVE_INT);
+    H5Tinsert(hdf5_fits_tid, "energy3", HOFFSET(HDF5Fit, energy3), H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "theta3",  HOFFSET(HDF5Fit, theta3),  H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "phi3",    HOFFSET(HDF5Fit, phi3),    H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "fmin",    HOFFSET(HDF5Fit, fmin),    H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "time",    HOFFSET(HDF5Fit, time),    H5T_NATIVE_DOUBLE);
+    H5Tinsert(hdf5_fits_tid, "psi",     HOFFSET(HDF5Fit, psi),     H5T_NATIVE_DOUBLE);
+
+    /*
+     * Create the dataset.
+     */
+    dataset = H5Dcreate2(file, "ev", hdf5_events_tid, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    dataset_mcgn = H5Dcreate2(file, "mcgn", hdf5_mcgn_tid, space_mcgn, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    dataset_fits = H5Dcreate2(file, "fits", hdf5_fits_tid, space_fits, H5P_DEFAULT, dcpl_id, H5P_DEFAULT);
+
+    /*
+     * Wtite data to the dataset;
+     */
+    if (nevents)
+        status = H5Dwrite(dataset, hdf5_events_tid, H5S_ALL, H5S_ALL, H5P_DEFAULT, hdf5_events);
+    if (nmcgn)
+        status = H5Dwrite(dataset_mcgn, hdf5_mcgn_tid, H5S_ALL, H5S_ALL, H5P_DEFAULT, hdf5_mcgn);
+    if (nfits)
+        status = H5Dwrite(dataset_fits, hdf5_fits_tid, H5S_ALL, H5S_ALL, H5P_DEFAULT, hdf5_fits);
+
+    /*
+     * Release resources
+     */
+    H5Sclose(aid1);
+    H5Sclose(aid2);
+    H5Tclose(atype1);
+    H5Tclose(atype2);
+    H5Aclose(attr1);
+    H5Aclose(attr2);
+    H5Tclose(hdf5_events_tid);
+    H5Tclose(hdf5_mcgn_tid);
+    H5Tclose(hdf5_fits_tid);
+    H5Sclose(space);
+    H5Sclose(space_mcgn);
+    H5Sclose(space_fits);
+    H5Pclose(dcpl_id);
+    H5Dclose(dataset);
+    H5Dclose(dataset_mcgn);
+    H5Dclose(dataset_fits);
+    H5Fclose(file);
+
+    return 0;
+}
diff --git a/src/hdf5_utils.h b/src/hdf5_utils.h
new file mode 100644
index 0000000..4f6ad32
--- /dev/null
+++ b/src/hdf5_utils.h
@@ -0,0 +1,85 @@
+#ifndef HDF5_UTILS
+#define HDF5_UTILS
+
+#include <stdint.h>
+#include <stdlib.h> /* for size_t */
+
+/* Struct to hold information for each triggered event (similar to the EV bank,
+ * but with some reconstructed information).
+ *
+ * This is the struct that is written out to the 'ev' dataset in the HDF5 file.
+ *
+ * Note: If there is no FTP, FTK, or RSP info in the zdab file these values
+ * will be set to nan. */
+typedef struct HDF5Event {
+    int run;
+    int evn;
+    double gtr;
+    int nhit;
+    uint32_t gtid;
+    uint32_t trg_type;
+    uint32_t dc;
+    double ftp_x;
+    double ftp_y;
+    double ftp_z;
+    double ftk_energy;
+    double rsp_energy;
+} HDF5Event;
+
+/* Struct to hold information for each primary seed tracks that is written out
+ * as the `mcgn` dataset in the HDF5 file. To figure out which track
+ * corresponds to a given GTID you can match them up based on the evn column,
+ * which is filled with the MC event number in both datasets. */
+typedef struct HDF5MCGN {
+    int run;
+    int evn;
+    int id;
+    double energy;
+    double x;
+    double y;
+    double z;
+    double dirx;
+    double diry;
+    double dirz;
+    double time;
+} HDF5MCGN;
+
+/* Struct to hold information for each fit that is written out
+ * as the `fits` dataset in the HDF5 file.
+ *
+ * Note that the energy2, theta2, phi2 fields will only be filled in for fits
+ * with 2 tracks and similarly for energy3, theta3, and phi3 for 3 track fits. */
+typedef struct HDF5Fit {
+    int run;
+    uint32_t gtid;
+    double x;
+    double y;
+    double z;
+    double t0;
+
+    /* Number of tracks. */
+    int n;
+
+    int id1;
+    double energy1;
+    double theta1;
+    double phi1;
+
+    int id2;
+    double energy2;
+    double theta2;
+    double phi2;
+
+    int id3;
+    double energy3;
+    double theta3;
+    double phi3;
+
+    double fmin;
+    double time;
+    double psi;
+} HDF5Fit;
+
+int save_output(const char *output, HDF5Event *hdf5_events, size_t nevents, HDF5MCGN *hdf5_mcgn, size_t nmcgn, HDF5Fit *hdf5_fits, size_t nfits);
+
+#endif
diff --git a/src/zdab-cat.c b/src/zdab-cat.c
index 5b9bb3d..278236c 100644
--- a/src/zdab-cat.c
+++ b/src/zdab-cat.c
@@ -26,6 +26,15 @@
 #include <errno.h> /* for errno */
 #include "release.h"
 #include "dc.h"
+#include "hdf5.h"
+#include <math.h>
+#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'
@@ -47,22 +56,25 @@ int main(int argc, char **argv)
 {
     int i;
     zebraFile *f;
-    zebraBank bmast, bmc, bmcgn, mctk, b;
+    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;
-    FILE *fout = stdout;
     int skip_second_event = 0;
-    size_t nhit;
     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];
 
     for (i = 1; i < argc; i++) {
         if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) {
@@ -96,23 +108,8 @@ int main(int argc, char **argv)
         return 1;
     }
 
-    if (output) {
-        fout = fopen(output, "w");
-
-        if (!fout) {
-            fprintf(stderr, "failed to open '%s': %s\n", output, strerror(errno));
-            return 1;
-        }
-    }
-
-    if (fout) {
-        fprintf(fout, "git_sha1: %s\n", GitSHA1());
-        fprintf(fout, "git_dirty: %s\n", GitDirty());
-    }
-
     if (load_pmt_info()) {
         zebra_close(f);
-        if (output) fclose(fout);
         return 1;
     }
 
@@ -160,11 +157,9 @@ int main(int argc, char **argv)
             continue;
         }
 
-        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]);
+        rv = zebra_get_bank(f,&mc,bmast.links[KMAST_MC-1]);
 
         if (rv) {
             fprintf(stderr, "error getting MC bank: %s\n", zebra_err);
@@ -176,14 +171,15 @@ int main(int argc, char **argv)
             goto err;
         }
 
-        rv = zebra_get_bank(f,&bmcgn,bmc.links[KMC_MCGN-1]);
+        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;
         }
 
-        if (fout) fprintf(fout, "    mcgn:\n");
         while (1) {
             if (bmcgn.links[KMCGN_MCTK-1] == 0) {
                 fprintf(stderr, "MCTK link is zero!\n");
@@ -232,17 +228,18 @@ int main(int argc, char **argv)
 
             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);
-                fprintf(fout, "        time: %.4f\n", bmcvx.tim);
+            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) {
@@ -278,7 +275,6 @@ skip_mc:
             }
         }
 
-        if (fout) fprintf(fout, "    ev:\n");
         while (1) {
             unpack_ev(b.data, &bev);
             ev.run = bev.run;
@@ -305,44 +301,41 @@ skip_mc:
 
             rv = get_event(f,&ev,&b);
 
-            nhit = get_nhit(&ev);
+            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 (fout) {
-                fprintf(fout, "      - run: %i\n", ev.run);
-                fprintf(fout, "        gtr: %.0f\n", ev.trigger_time);
-                fprintf(fout, "        nhit: %zu\n", nhit);
-                fprintf(fout, "        gtid: %i\n", ev.gtid);
-                fprintf(fout, "        trg_type: 0x%08x\n", ev.trigger_type);
-                fprintf(fout, "        dc: 0x%08x\n", get_dc_word(&ev, f, &bmast, &b));
-            }
-
-            if (fout) {
                 if (get_ftpv(f,&b,&bftpv)) {
                     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 {
-                    fprintf(fout, "        ftp:\n");
-                    fprintf(fout, "            x: %.2f\n", bftpv.x);
-                    fprintf(fout, "            y: %.2f\n", bftpv.y);
-                    fprintf(fout, "            z: %.2f\n", bftpv.z);
+                    hdf5_events[nevents].ftp_x = bftpv.x;
+                    hdf5_events[nevents].ftp_y = bftpv.y;
+                    hdf5_events[nevents].ftp_z = bftpv.z;
                 }
-            }
 
-            if (fout) {
                 if (get_ftxk(f,&b,&bftxk)) {
                     fprintf(stderr, "%s\n", zdab_err);
+                    hdf5_events[nevents].ftk_energy = NAN;
                 } else {
-                    fprintf(fout, "        ftk:\n");
-                    fprintf(fout, "            energy: %.2f\n", bftxk.energy);
+                    hdf5_events[nevents].ftk_energy = bftxk.energy;
                 }
-            }
 
-            if (fout) {
                 if (get_rsp(f,&b,&bftxr)) {
                     fprintf(stderr, "%s\n", zdab_err);
+                    hdf5_events[nevents].rsp_energy = NAN;
                 } else {
-                    fprintf(fout, "        rsp:\n");
-                    fprintf(fout, "            energy: %.2f\n", bftxr.ene);
+                    hdf5_events[nevents].rsp_energy = bftxr.ene;
                 }
+
+                nevents++;
             }
 
             /* Note the origin link for the first EV bank points back to the
@@ -359,9 +352,11 @@ skip_mc:
         }
     }
 
-    db_free(db);
+    if (output) {
+        save_output(output, hdf5_events, nevents, hdf5_mcgn, nmcgn, NULL, 0);
+    }
 
-    if (fout) fclose(fout);
+    db_free(db);
 
     zebra_close(f);
 
@@ -370,8 +365,6 @@ skip_mc:
 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 716503c..416bd79 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -580,6 +580,24 @@ void unpack_ftk(uint32_t *data, FTKBank *b)
     unpack((uint8_t *) (data+25), "f",&b->spare5);
 }
 
+void unpack_mc(uint32_t *data, MCBank *b)
+{
+    unpack((uint8_t *) data,      "l",&b->jdy);
+    unpack((uint8_t *) (data+1),  "l",&b->ut1);
+    unpack((uint8_t *) (data+2),  "l",&b->ut2);
+    unpack((uint8_t *) (data+3),  "l",&b->dte);
+    unpack((uint8_t *) (data+4),  "l",&b->hmsc);
+    unpack((uint8_t *) (data+5),  "l",&b->seed1);
+    unpack((uint8_t *) (data+6),  "l",&b->seed2);
+    unpack((uint8_t *) (data+7),  "l",&b->seed_num);
+    unpack((uint8_t *) (data+8),  "l",&b->dtp);
+    unpack((uint8_t *) (data+9),  "f",&b->mcver);
+    unpack((uint8_t *) (data+10), "l",&b->evn);
+    unpack((uint8_t *) (data+11), "F",&b->gtr);
+    unpack((uint8_t *) (data+13), "l",&b->num_ge_err);
+    unpack((uint8_t *) (data+14), "f",&b->pmt_eff);
+}
+
 void unpack_mcgn(uint32_t *data, MCGNBank *b)
 {
     unpack((uint8_t *) data,      "l",&b->id);
diff --git a/src/zdab_utils.h b/src/zdab_utils.h
index d4411d9..043e185 100644
--- a/src/zdab_utils.h
+++ b/src/zdab_utils.h
@@ -378,6 +378,37 @@ typedef struct FTKBank {
     float spare5;
 } FTKBank;
 
+typedef struct MCBank {
+    /* Julian date. */
+    uint32_t jdy;
+    /* Universal time seconds. */
+    uint32_t ut1;
+    /* Universal time nanoseconds. */
+    uint32_t ut2;
+    /* Date (format: yyyymmdd). */
+    uint32_t dte;
+    /* Time (format: hhmmsscc  - cc is centisec). */
+    uint32_t hmsc;
+    /* First random number seed. */
+    uint32_t seed1;
+    /* Second random number seed. */
+    uint32_t seed2;
+    /* Random number number. */
+    uint32_t seed_num;
+    /* Data type. */
+    uint32_t dtp;
+    /* SNOMAN MC version number. */
+    float mcver;
+    /* Monte Carlo event number. */
+    uint32_t evn;
+    /* Generation time in nsec (first word). */
+    double gtr;
+    /* Number of suppressed Cerenkov photon tracking errors. */
+    uint32_t num_ge_err;
+    /* The current PMT collection efficiency. */
+    float pmt_eff;
+} MCBank;
+
 typedef struct MCGNBank {
     /* Particle id. */
     uint32_t id;
@@ -638,6 +669,7 @@ void unpack_ftpt(uint32_t *data, FTPTBank *b);
 void unpack_ftpv(uint32_t *data, FTPVBank *b);
 void unpack_ftxk(uint32_t *data, FTXKBank *b);
 void unpack_ftk(uint32_t *data, FTKBank *b);
+void unpack_mc(uint32_t *data, MCBank *b);
 void unpack_mcgn(uint32_t *data, MCGNBank *b);
 void unpack_mcvx(uint32_t *data, MCVXBank *b);
 void unpack_mctk(uint32_t *data, MCTKBank *b);