update dc scripts

- delete instrumental tags in both scripts since we already tag them in
  get_events()
- apply the same cuts to the MC as to the data
- add a comment about the fitted fraction
- quit if we don't have at least 1 type of each instrumental

1 files changed, 10 insertions, 19 deletions

diff --git a/utils/dc b/utils/dc
index 0cc4ba3..a267090 100755
--- a/utils/dc
+++ b/utils/dc
@@ -290,17 +290,19 @@ if __name__ == '__main__':
     ev = ev[ev.prompt]
     ev = ev[ev.nhit_cal > 100]
 
+    # Note: Technically we want to know the fitted fraction only for events
+    # which *would* reconstruct above 20 MeV. However, there is no way to know
+    # if the energy is above 20 MeV without fitting it. However, since I only
+    # skip fitting events based on the gtid, there shouldn't be any correlation
+    # with energy and so the fitted fraction here should be correct.
     fitted_fraction = {}
     for bg in ['signal','muon','noise','neck','flasher','breakdown']:
         if np.count_nonzero(ev[bg]):
             fitted_fraction[bg] = np.count_nonzero(ev[bg] & ~np.isnan(ev.fmin))/np.count_nonzero(ev[bg])
             print("Fitted fraction for %s: %.0f %%" % (bg,fitted_fraction[bg]*100))
-        elif bg != 'signal':
-            print_warning("Warning: No %s events in sample, assuming 10%% fitted fraction" % bg)
-            fitted_fraction[bg] = 0.1
         else:
-            print_warning("Warning: No signal events in sample, assuming 100%% fitted fraction")
-            fitted_fraction[bg] = 1.0
+            print_warning("Warning: No %s events in sample!" % bg)
+            sys.exit(1)
 
     ev = ev[~np.isnan(ev.fmin)]
     ev = ev[ev.ke > 20]
@@ -312,13 +314,6 @@ if __name__ == '__main__':
     ev = z_cut(ev)
     ev = udotr_cut(ev)
 
-    ev['noise'] = ev.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_ITC | DC_ESUM) != 0
-    ev['neck'] = ((ev.dc & DC_NECK) != 0) & ~ev.noise
-    ev['flasher'] = ((ev.dc & DC_FLASHER) != 0) & ~(ev.noise | ev.neck) & (ev.nhit < 1000)
-    ev['breakdown'] = ((ev.dc & (DC_FLASHER | DC_BREAKDOWN)) != 0) & ~(ev.noise | ev.neck) & (ev.nhit >= 1000)
-    ev['muon'] = ((ev.dc & DC_MUON) != 0) & ~(ev.noise | ev.neck | ev.flasher | ev.breakdown)
-    ev['signal'] = ~(ev.noise | ev.neck | ev.flasher | ev.breakdown | ev.muon)
-
     data = {}
     for bg in ['signal','muon','noise','neck','flasher','breakdown']:
         data[bg] = np.zeros((2,2,2,2),dtype=int)
@@ -328,6 +323,9 @@ if __name__ == '__main__':
     ev_mc = get_events(args.mc, merge_fits=True)
 
     ev_mc = ev_mc[ev_mc.prompt]
+    ev_mc = ev_mc[ev_mc.nhit_cal > 100]
+    ev_mc = ev_mc[~np.isnan(ev_mc.fmin)]
+    ev_mc = ev_mc[ev_mc.ke > 20]
 
     # figure out bins for high level variables
     ev_mc = radius_cut(ev_mc)
@@ -336,13 +334,6 @@ if __name__ == '__main__':
     ev_mc = z_cut(ev_mc)
     ev_mc = udotr_cut(ev_mc)
 
-    ev_mc['noise'] = ev_mc.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_ITC | DC_ESUM) != 0
-    ev_mc['neck'] = ((ev_mc.dc & DC_NECK) != 0) & ~ev_mc.noise
-    ev_mc['flasher'] = ((ev_mc.dc & DC_FLASHER) != 0) & ~(ev_mc.noise | ev_mc.neck) & (ev_mc.nhit < 1000)
-    ev_mc['breakdown'] = ((ev_mc.dc & (DC_FLASHER | DC_BREAKDOWN)) != 0) & ~(ev_mc.noise | ev_mc.neck) & (ev_mc.nhit >= 1000)
-    ev_mc['muon'] = ((ev_mc.dc & DC_MUON) != 0) & ~(ev_mc.noise | ev_mc.neck | ev_mc.flasher | ev_mc.breakdown)
-    ev_mc['signal'] = ~(ev_mc.noise | ev_mc.neck | ev_mc.flasher | ev_mc.breakdown | ev_mc.muon)
-
     # FIXME: Double check that what I'm calculating here matches with what I
     # expect
     sacrifice = {}