unwrap the 50 MHz clock in plot-energy

1 files changed, 53 insertions, 6 deletions

diff --git a/utils/plot-energy b/utils/plot-energy
index 3c645d5..917b5d6 100755
--- a/utils/plot-energy
+++ b/utils/plot-energy
@@ -101,6 +101,45 @@ def chunks(l, n):
 def print_warning(msg):
     print(bcolors.FAIL + msg + bcolors.ENDC,file=sys.stderr)
 
+def unwrap(p, delta, axis=-1):
+    """
+    A modified version of np.unwrap() useful for unwrapping the 50 MHz clock.
+    It unwraps discontinuities bigger than delta/2 by delta.
+
+    Example:
+
+        >>> a = np.arange(10) % 5
+        >>> a
+        array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4])
+        >>> unwrap(a,5)
+        array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.])
+
+    In the case of the 50 MHz clock delta should be 0x7ffffffffff*20.0.
+    """
+    p = np.asarray(p)
+    nd = p.ndim
+    dd = np.diff(p, axis=axis)
+    slice1 = [slice(None, None)]*nd     # full slices
+    slice1[axis] = slice(1, None)
+    slice1 = tuple(slice1)
+    ddmod = np.mod(dd + delta/2, delta) - delta/2
+    np.copyto(ddmod, delta/2, where=(ddmod == -delta/2) & (dd > 0))
+    ph_correct = ddmod - dd
+    np.copyto(ph_correct, 0, where=abs(dd) < delta/2)
+    up = np.array(p, copy=True, dtype='d')
+    up[slice1] = p[slice1] + ph_correct.cumsum(axis)
+    return up
+
+def unwrap_50_mhz_clock(gtr):
+    """
+    Unwrap an array with 50 MHz clock times. These times should all be in
+    nanoseconds and come from the KEV_GTR variable in the EV bank.
+
+    Note: We assume here that the events are already ordered contiguously by
+    GTID, so you shouldn't pass an array with multiple runs!
+    """
+    return unwrap(gtr,0x7ffffffffff*20.0)
+
 if __name__ == '__main__':
     import argparse
     import matplotlib.pyplot as plt
@@ -116,6 +155,16 @@ if __name__ == '__main__':
     ev = pd.concat([pd.read_hdf(filename, "ev") for filename in args.filenames],ignore_index=True)
     fits = pd.concat([pd.read_hdf(filename, "fits") for filename in args.filenames],ignore_index=True)
 
+    # First, remove junk events since orphans won't have a 50 MHz clock and so
+    # 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'])
+
+    # unwrap the 50 MHz clock within each run
+    ev.gtr = ev.groupby(['run'],as_index=False)['gtr'].transform(unwrap_50_mhz_clock)
+
     # This is a bit of a hack. It appears that many times the fit will
     # actually do much better by including a very low energy electron or
     # muon. I believe the reason for this is that of course my likelihood
@@ -179,9 +228,6 @@ if __name__ == '__main__':
     #
     # We define a prompt event here as any event with an NHIT > 100 and whose
     # previous > 100 nhit event was more than 250 ms ago
-    #
-    # Should I use 10 MHz clock here since the time isn't critical and then I
-    # don't have to worry about 50 MHz clock rollover?
     prompt = ev[ev.nhit >= 100]
     prompt_mask = np.concatenate(([True],np.diff(prompt.gtr.values) > 250e6))
 
@@ -236,7 +282,10 @@ if __name__ == '__main__':
             else:
                 michel_sum = michel_sum.append(michel_run)
 
-        michel = michel_sum
+        if michel_sum:
+            michel = michel_sum
+        else:
+            michel = pd.DataFrame(columns=follower.columns)
     else:
         michel = pd.DataFrame(columns=follower.columns)
 
@@ -249,7 +298,6 @@ if __name__ == '__main__':
         neutron = neutron[neutron.rsp_energy > 4.0]
 
         # neutron events accepted after 20 microseconds and before 250 ms (50 ms during salt)
-        # FIXME: need to deal with 50 MHz clock rollover
         prompt_sum = None
         atm_sum = None
         for run, prompt_run in prompt.groupby(['run']):
@@ -282,7 +330,6 @@ if __name__ == '__main__':
     print("number of events after muon cut = %i" % len(prompt))
 
     if muons.size:
-        # FIXME: need to deal with 50 MHz clock rollover
         # remove events 200 microseconds after a muon
         prompt = prompt[~np.any((prompt.run.values == muons.run.values[:,np.newaxis]) & \
                                 (prompt.gtr.values > muons.gtr.values[:,np.newaxis]) & \