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.

1 files changed, 33 insertions, 12 deletions

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]