calculate livetime based on pulse gt and 10 MHz clock

This commit updates get_events() to calculate the livetime based on both
the number of pulse gt events and the 10 MHz clock and to return it in a
dictionary stored with the dataframe.

I also update dm-search so that the results are now reported as a
function of events/cm^3/s.

Also updated radius cut to be the AV radius.

3 files changed, 47 insertions, 11 deletions

diff --git a/utils/chi2 b/utils/chi2
index 0326e1a..39db081 100755
--- a/utils/chi2
+++ b/utils/chi2
@@ -480,10 +480,10 @@ if __name__ == '__main__':
     ev_mc = ev_mc[~np.isnan(ev_mc.fmin)]
     muon_mc = muon_mc[~np.isnan(muon_mc.fmin)]
 
-    # require (r/r_psup)^3 < 0.9
-    ev = ev[ev.r_psup < 0.9]
-    ev_mc = ev_mc[ev_mc.r_psup < 0.9]
-    muon_mc = muon_mc[muon_mc.r_psup < 0.9]
+    # require (r < av radius)
+    ev = ev[ev.r < AV_RADIUS]
+    ev_mc = ev_mc[ev_mc.r < AV_RADIUS]
+    muon_mc = muon_mc[muon_mc.r < AV_RADIUS]
 
     # require psi < 6
     ev = ev[ev.psi < 6]
diff --git a/utils/dm-search b/utils/dm-search
index a486593..25fa185 100755
--- a/utils/dm-search
+++ b/utils/dm-search
@@ -35,7 +35,7 @@ from sddm import printoptions
 from sddm.utils import fast_cdf, correct_energy_bias
 from scipy.integrate import quad
 from sddm.dm import *
-from sddm import SNOMAN_MASS
+from sddm import SNOMAN_MASS, AV_RADIUS
 import nlopt
 
 # Likelihood Fit Parameters
@@ -494,6 +494,18 @@ if __name__ == '__main__':
         evs.append(get_events(df.filename.values, merge_fits=True, nhit_thresh=args.nhit_thresh))
     ev = pd.concat(evs)
 
+    livetime = 0.0
+    livetime_pulse_gt = 0.0
+    for ev in evs:
+        if not np.isnan(ev.attrs['time_10_mhz']):
+            livetime += ev.attrs['time_10_mhz']
+        else:
+            livetime += ev.attrs['time_pulse_gt']
+        livetime_pulse_gt += ev.attrs['time_pulse_gt']
+
+    print("livetime            = %.2f" % livetime)
+    print("livetime (pulse gt) = %.2f" % livetime_pulse_gt)
+
     ev = correct_energy_bias(ev)
 
     # Note: We loop over the MC filenames here instead of just passing the
@@ -530,10 +542,12 @@ if __name__ == '__main__':
     ev_mc = ev_mc[~np.isnan(ev_mc.fmin)]
     muon_mc = muon_mc[~np.isnan(muon_mc.fmin)]
 
-    # require (r/r_psup)^3 < 0.9
-    ev = ev[ev.r_psup < 0.9]
-    ev_mc = ev_mc[ev_mc.r_psup < 0.9]
-    muon_mc = muon_mc[muon_mc.r_psup < 0.9]
+    # require (r < av radius)
+    ev = ev[ev.r < AV_RADIUS]
+    ev_mc = ev_mc[ev_mc.r < AV_RADIUS]
+    muon_mc = muon_mc[muon_mc.r < AV_RADIUS]
+
+    fiducial_volume = (4/3)*np.pi*(AV_RADIUS)**3
 
     # require psi < 6
     ev = ev[ev.psi < 6]
@@ -702,11 +716,11 @@ if __name__ == '__main__':
 
     fig = plt.figure()
     for color, dm_particle_id in zip(('C0','C1'),(2020,2222)):
-        plt.plot(DM_MASSES[dm_particle_id],limits[dm_particle_id],color=color,label='$' + ''.join([particle_id[int(''.join(x))] for x in grouper(str(dm_particle_id),2)]) + '$')
+        plt.plot(DM_MASSES[dm_particle_id],limits[dm_particle_id]/fiducial_volume/livetime,color=color,label='$' + ''.join([particle_id[int(''.join(x))] for x in grouper(str(dm_particle_id),2)]) + '$')
     plt.gca().set_xscale("log")
     despine(fig,trim=True)
     plt.xlabel("Energy (MeV)")
-    plt.ylabel("Event Rate Limit (events)")
+    plt.ylabel("Event Rate Limit (events/cm^3/s)")
     plt.legend()
     plt.tight_layout()
 
diff --git a/utils/sddm/plot_energy.py b/utils/sddm/plot_energy.py
index be13d78..b16e2cd 100644
--- a/utils/sddm/plot_energy.py
+++ b/utils/sddm/plot_energy.py
@@ -42,6 +42,7 @@ TRIG_NHIT_100_MED =      0x00000002
 TRIG_NHIT_100_HI  =      0x00000004
 TRIG_NHIT_20      =      0x00000008
 TRIG_NHIT_20_LB   =      0x00000010
+TRIG_PULSE_GT     =      0x00000400
 
 def unwrap(p, delta, axis=-1):
     """
@@ -591,6 +592,19 @@ def get_events(filenames, merge_fits=False, nhit_thresh=None, mc=False):
     ev['signal'] = ~(ev.noise | ev.neck | ev.flasher | ev.breakdown | ev.muon)
     ev['instrumental'] = ~ev.signal
 
+    n_pulse_gt = np.count_nonzero((ev.trg_type & TRIG_PULSE_GT) != 0)
+    # Pulse GT went at approximately 5 Hz for the duration of SNO
+    time_pulse_gt = n_pulse_gt/5.0
+
+    if 'jdy' in ev.columns:
+        # Calculate total time based on 10 MHz clock
+        t10 = ev['jdy']*24*3600 + ev['ut1'] + ev['ut2']/1e9
+        # Note: We don't take the max and min here because there might be bit flips
+        # in the 10 MHz clock
+        time_10_mhz = np.percentile(t10,99.9) - np.percentile(t10,0.1)
+    else:
+        time_10_mhz = np.nan
+
     # Try to identify Michel electrons. Currently, the event selection is based
     # on Richie's thesis. Here, we do the following:
     #
@@ -634,6 +648,14 @@ def get_events(filenames, merge_fits=False, nhit_thresh=None, mc=False):
                       np.cos(ev_single_particle.theta1)*ev_single_particle.z
         ev['udotr'] /= ev_single_particle.r
 
+        # Yes, this is super hacky, but I don't want to change what's returned from
+        # this function.
+        ev.attrs = {'time_pulse_gt': time_pulse_gt, 'time_10_mhz': time_10_mhz}
+
         return ev
 
+    # Yes, this is super hacky, but I don't want to change what's returned from
+    # this function.
+    ev.attrs = {'time_pulse_gt': time_pulse_gt, 'time_10_mhz': time_10_mhz}
+
     return ev, fits