update retrigger cut

This commit updates the retrigger cut to cut events where the previous
event is missing, so that even if I forget to run the analysis with all
the orphan events included, we will cut events potentially coming after
an instrumental or muon that got cut by the junk cut.

1 files changed, 16 insertions, 4 deletions

diff --git a/utils/sddm/plot_energy.py b/utils/sddm/plot_energy.py
index 7852412..efbaa36 100644
--- a/utils/sddm/plot_energy.py
+++ b/utils/sddm/plot_energy.py
@@ -86,7 +86,7 @@ def retrigger_cut(ev):
     """
     Cuts all retrigger events.
     """
-    return ev[ev.dt > 500]
+    return ev[(ev.dt > 500) & (ev.gtid_diff == 0)]
 
 def breakdown_follower_cut(ev):
     """
@@ -459,7 +459,13 @@ def get_events(filenames, merge_fits=False, nhit_thresh=None, apply_nhit_trigger
 
     # First, remove orphans since they won't have a 50 MHz clock and so could
     # screw up the 50 MHz clock unwrapping
-    ev = ev[(ev.gtid == 0) & (ev.gtr == 0)]
+    ev = ev[~((ev.gtid == 0) & (ev.gtr == 0) & (ev.trg_type == 0))]
+
+    # Calculate the expected next GTID for all GTIDs. The algorithm here is
+    # taken from the builder's eb_make_entry() function in eb_queue_ph.c.
+    next_gtid = np.roll(np.arange(0x1000000),-1)
+    next_gtid[(next_gtid & 0xffff) == 0xffff] += 1
+    next_gtid[(next_gtid & 0xff0000) == 0xff0000] = 0
 
     # 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
@@ -498,13 +504,19 @@ def get_events(filenames, merge_fits=False, nhit_thresh=None, apply_nhit_trigger
         # 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))
+        jumps = ev_run.gtid.values[1:] - next_gtid[ev_run.gtid[:-1].values]
+        if np.count_nonzero(jumps != 0):
+            print(ev_run[jumps != 0][['run','gtid','gtid_diff']])
+            print_warning("Warning: %i GTID jumps in run %i" % (np.count_nonzero(jumps != 0),run))
 
     # calculate the time difference between each event and the previous event
     # so we can tag retrigger events
     ev['dt'] = ev.groupby(['run'],group_keys=False)['gtr'].transform(lambda x: np.concatenate(([1e9],np.diff(x.values))))
 
+    # calculate the time difference between each event and the previous event
+    # so we can tag retrigger events
+    ev['gtid_diff'] = ev.groupby(['run'],group_keys=False)['gtid'].transform(lambda x: np.concatenate(([0],x.values[1:]-next_gtid[x[:-1].values])))
+
     # Calculate the approximate Ockham factor.
     # See Chapter 20 in "Probability Theory: The Logic of Science" by Jaynes
     #