add sub_run variable to the events array in the HDF5 file

This commit adds the sub_run variable to the ev array in the HDF5 output file
and updates plot-energy to order the events using the run and sub_run
variables. This fixes a potential issue where I was sorting by GTID before, but
the GTID can wrap around and so isn't guaranteed to put the events in the right
order.

8 files changed, 41 insertions, 13 deletions

diff --git a/src/event.h b/src/event.h
index 88fcbe1..fbb1658 100644
--- a/src/event.h
+++ b/src/event.h
@@ -73,6 +73,7 @@ typedef struct pmt_hit {
 
 typedef struct event {
     int run;
+    int sub_run;
     /* Global trigger time in ns. */
     double trigger_time;
     /* Number of normal hit PMTs. */
diff --git a/src/fit.c b/src/fit.c
index 056f6b0..7fdf994 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -6299,6 +6299,7 @@ skip_mc:
 
             if (output) {
                 hdf5_events[nevents].run = ev.run;
+                hdf5_events[nevents].sub_run = ev.sub_run;
                 hdf5_events[nevents].evn = bmc.evn;
                 hdf5_events[nevents].gtr = ev.trigger_time;
                 hdf5_events[nevents].nhit = ev.nhit;
diff --git a/src/hdf5_utils.c b/src/hdf5_utils.c
index 05d1002..2a8f396 100644
--- a/src/hdf5_utils.c
+++ b/src/hdf5_utils.c
@@ -102,6 +102,7 @@ int save_output(const char *output, HDF5Event *hdf5_events, size_t nevents, HDF5
      */
     hdf5_events_tid = H5Tcreate(H5T_COMPOUND, sizeof(HDF5Event));
     H5Tinsert(hdf5_events_tid, "run",        HOFFSET(HDF5Event, run),        H5T_NATIVE_INT);
+    H5Tinsert(hdf5_events_tid, "sub_run",    HOFFSET(HDF5Event, sub_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);
diff --git a/src/hdf5_utils.h b/src/hdf5_utils.h
index 4f6ad32..850355f 100644
--- a/src/hdf5_utils.h
+++ b/src/hdf5_utils.h
@@ -13,6 +13,7 @@
  * will be set to nan. */
 typedef struct HDF5Event {
     int run;
+    int sub_run;
     int evn;
     double gtr;
     int nhit;
diff --git a/src/zdab-cat.c b/src/zdab-cat.c
index dfb295a..7b0b714 100644
--- a/src/zdab-cat.c
+++ b/src/zdab-cat.c
@@ -170,6 +170,7 @@ int main(int argc, char **argv)
              * instead of a MAST record. */
             continue;
         } else {
+            continue;
             fprintf(stderr, "logical record starts with unknown bank '%s'\n", bmast.name);
             goto err;
         }
@@ -323,6 +324,7 @@ skip_mc:
 
             if (output) {
                 hdf5_events[nevents].run = ev.run;
+                hdf5_events[nevents].sub_run = ev.sub_run;
                 hdf5_events[nevents].evn = bmc.evn;
                 hdf5_events[nevents].gtr = ev.trigger_time;
                 hdf5_events[nevents].nhit = ev.nhit;
diff --git a/src/zdab_utils.c b/src/zdab_utils.c
index cf1a01c..8b6cbc8 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -103,6 +103,7 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
     unpack_ev(bev->data, &ev_bank);
 
     ev->run = ev_bank.run;
+    ev->sub_run = ev_bank.sub_run;
     ev->gtid = ev_bank.gtr_id;
     ev->trigger_type = ev_bank.trg_type;
     ev->trigger_time = ev_bank.gtr;
diff --git a/src/zdab_utils.h b/src/zdab_utils.h
index ce3661c..91f8898 100644
--- a/src/zdab_utils.h
+++ b/src/zdab_utils.h
@@ -526,7 +526,7 @@ typedef struct EVBank {
     /* Integrated charge in the event. */
     float nph;
     /* Sub-run number, or -1 if unknown. Only defined for run number >= 10614. */
-    uint32_t sub_run;
+    int32_t sub_run;
     /* Packed word. */
     uint32_t mc_pck;
     /* ZDAB input/output record type, e.g. PMT, NCD, CMA, RUN, etc. */
diff --git a/utils/plot-energy b/utils/plot-energy
index 6d2da9d..896a6e5 100755
--- a/utils/plot-energy
+++ b/utils/plot-energy
@@ -183,12 +183,40 @@ if __name__ == '__main__':
     # could screw up the 50 MHz clock unwrapping
     ev = ev[ev.dc & DC_JUNK == 0]
 
-    # make sure events are sorted by GTID
-    ev = ev.sort_values(by=['run','gtid'])
+    # We need the events to be in time order here in order to calculate the
+    # delta t between events. It's not obvious exactly how to do this. You
+    # could sort by GTID, but that wraps around. Similarly we can't sort by the
+    # 50 MHz clock because it also wraps around. Finally, I'm hesitant to sort
+    # by the 10 MHz clock since it can be unreliable.
+    #
+    # Therefore, we just assume that the events are already in order in the
+    # zdab file. Later, we check that the GTIDs make sense.
+    #
+    # Note: There are no sub run IDs in the runs before run 10614, so it's
+    # necessary to specify those on the command line in the order of the
+    # subrun.
+    #
+    # Note: We use mergesort here since it's the only stable sort and we want
+    # to preserve the order of the events.
+    ev = ev.sort_values(by=['run','sub_run'],kind='mergesort')
+
+    for run, ev_run in ev.groupby('run'):
+        # Warn about 50 MHz clock jumps since they could indicate that the
+        # events aren't in order.
+        dt = np.diff(ev_run.gtr)
+        if np.count_nonzero(np.abs(dt) > 1e9):
+            print_warning("Warning: %i 50 MHz clock jumps in run %i. Are the events in order?" % (np.count_nonzero(np.abs(dt) > 1e9),run))
 
     # unwrap the 50 MHz clock within each run
     ev.gtr = ev.groupby(['run'],as_index=False)['gtr'].transform(unwrap_50_mhz_clock)
 
+    for run, ev_run in ev.groupby('run'):
+        # Warn about GTID jumps since we could be missing a potential flasher
+        # and/or breakdown, and we need all the events in order to do a
+        # retrigger cut
+        if np.count_nonzero(np.diff(ev_run.gtid) != 1):
+            print_warning("Warning: %i GTID jumps in run %i" % (np.count_nonzero(np.diff(ev_run.gtid) != 1),run))
+
     # calculate the time difference between each event and the previous event
     # so we can tag retrigger events
     ev['dt'] = ev.groupby(['run'],as_index=False)['gtr'].transform(lambda x: np.concatenate(([1e9],np.diff(x.values))))
@@ -231,23 +259,16 @@ if __name__ == '__main__':
     fits.loc[fits['n'] == 3, 'id'] = fits['id1']*10000 + fits['id2']*100 + fits['id3']
     fits['theta'] = fits['theta1']
 
-    # Make sure events are in order. We use run number and GTID here which
-    # should be in time order as well except for bit flips in the GTID
-    # This is important for later when we look at time differences in the 50
-    # MHz clock. We should perhaps do a check here based on the 10 MHz clock to
-    # make sure things are in order
-    ev = ev.sort_values(by=['run','gtid'])
-
     print("number of events = %i" % len(ev))
 
     # flasher follower cut
-    ev = ev.groupby(['run'],as_index=False).apply(flasher_follower_cut)
+    ev = ev.groupby('run',as_index=False).apply(flasher_follower_cut).reset_index(level=0,drop=True)
 
     # breakdown follower cut
-    ev = ev.groupby(['run'],as_index=False).apply(breakdown_follower_cut)
+    ev = ev.groupby('run',as_index=False).apply(breakdown_follower_cut).reset_index(level=0,drop=True)
 
     # retrigger cut
-    ev = ev.groupby(['run'],as_index=False).apply(retrigger_cut)
+    ev = ev.groupby('run',as_index=False).apply(retrigger_cut).reset_index(level=0,drop=True)
 
     # First, do basic data cleaning which is done for all events.
     ev = ev[ev.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_FLASHER | DC_NECK | DC_ITC | DC_BREAKDOWN) == 0]