plot kinetic energy instead of calresp0 variable. Also update the expected number of events to 230 per year since that's what's in Richie's thesis

1 files changed, 36 insertions, 19 deletions

diff --git a/analyze_genie_mc.py b/analyze_genie_mc.py
index 0a83fbb..8c006a3 100755
--- a/analyze_genie_mc.py
+++ b/analyze_genie_mc.py
@@ -149,6 +149,7 @@ if __name__ == '__main__':
     total_neutrons = []
     total_neutrons_detected = []
     E = []
+    KE = []
     r = []
     total_nrings = []
     total_e_like_rings = []
@@ -164,6 +165,7 @@ if __name__ == '__main__':
             nrings = 0
             e_like_rings = 0
             u_like_rings = 0
+            ke = 0
 
             if event.cc:
                 if abs(event.neu) == 12:
@@ -171,6 +173,7 @@ if __name__ == '__main__':
                 else:
                     u_like_rings = 1
                 nrings = 1
+                ke += event.El
             elif event.nc:
                 pass
             else:
@@ -187,6 +190,7 @@ if __name__ == '__main__':
                         # with > 100 MeV
                         nrings += 1
                         e_like_rings += 1
+                    ke += event.Ef[i]
                 elif pdg == 22:
                     # gamma
                     if event.Ef[i] > 0.1:
@@ -194,10 +198,12 @@ if __name__ == '__main__':
                         # 100 MeV
                         nrings += 1
                         e_like_rings += 1
+                    ke += event.Ef[i]
                 elif pdg == 111:
                     # pi0
                     nrings += 1
                     e_like_rings += 1
+                    ke += event.Ef[i]
                 elif abs(pdg) == 211:
                     # pi+/-
                     # momentum of ith particle in hadronic system
@@ -215,6 +221,10 @@ if __name__ == '__main__':
                         # produces 1 muon like ring
                         nrings += 1
                         u_like_rings += 1
+                    # FIXME: should actually be a beta distribution
+                    p = np.sqrt(event.pxf[i]**2 + event.pyf[i]**2 + event.pzf[i]**2)
+                    m = np.sqrt(event.Ef[i]**2 - p**2)
+                    ke += event.Ef[i] - m
                 elif abs(pdg) in [2212,3222,311,321,3122,3112]:
                     # momentum of ith particle in hadronic system
                     p = np.sqrt(event.pxf[i]**2 + event.pyf[i]**2 + event.pzf[i]**2)
@@ -225,6 +235,8 @@ if __name__ == '__main__':
                         # above cerenkov threshold
                         nrings += 1
                         u_like_rings += 1
+                        m = np.sqrt(event.Ef[i]**2 - p**2)
+                        ke += event.Ef[i] - m
                 elif int(("%010i" % abs(pdg))[0]) == 1:
                     # usually just excited 16O atom which won't produce a lot
                     # of light
@@ -238,6 +250,7 @@ if __name__ == '__main__':
             total_u_like_rings.append(u_like_rings)
             El.append(event.El)
             E.append(event.calresp0)
+            KE.append(ke)
             r.append(np.sqrt(event.vtxx**2 + event.vtxy**2 + event.vtxz**2))
 
             if total_neutrons_detected[-1] == 0 and nrings >= 2 and ((e_like_rings == 0) or (u_like_rings == 0)):
@@ -251,6 +264,7 @@ if __name__ == '__main__':
     total_neutrons = np.array(total_neutrons)
     total_neutrons_detected = np.array(total_neutrons_detected)
     E = np.array(E)
+    KE = np.array(KE)
     r = np.array(r)
     total_nrings = np.array(total_nrings)
     total_e_like_rings = np.array(total_e_like_rings)
@@ -266,6 +280,9 @@ if __name__ == '__main__':
     E1 = E[cut1]
     E2 = E[cut2]
     E3 = E[cut3]
+    KE1 = KE[cut1]
+    KE2 = KE[cut2]
+    KE3 = KE[cut3]
 
     plt.figure()
     bincenters = (bins[1:] + bins[:-1])/2
@@ -276,16 +293,16 @@ if __name__ == '__main__':
     # years based on Richie's thesis which says "Over the 306.4 live days of
     # the D2O phase we expect a total of 192.4 events within the acrylic vessel
     # and 504.5 events within the PSUP.
-    El_hist = El_hist*1000/total_events
+    El_hist = El_hist*230/total_events
     y = np.concatenate([[0],np.repeat(El_hist,2),[0]])
     El1_hist, _ = np.histogram(El1, bins=bins)
-    El1_hist = El1_hist*1000/total_events
+    El1_hist = El1_hist*230/total_events
     y1 = np.concatenate([[0],np.repeat(El1_hist,2),[0]])
     El2_hist, _ = np.histogram(El2, bins=bins)
-    El2_hist = El2_hist*1000/total_events
+    El2_hist = El2_hist*230/total_events
     y2 = np.concatenate([[0],np.repeat(El2_hist,2),[0]])
     El3_hist, _ = np.histogram(El3, bins=bins)
-    El3_hist = El3_hist*1000/total_events
+    El3_hist = El3_hist*230/total_events
     y3 = np.concatenate([[0],np.repeat(El3_hist,2),[0]])
     plt.plot(x, y, label="All events")
     plt.step(x, y1, where='mid', label="n")
@@ -301,35 +318,35 @@ if __name__ == '__main__':
     plt.title("Primary Lepton Energy")
 
     plt.figure()
-    E_hist, _ = np.histogram(E, bins=bins)
-    total_events = E_hist.sum()
+    KE_hist, _ = np.histogram(KE, bins=bins)
+    total_events = KE_hist.sum()
     # FIXME: this is just a rough estimate of how many events we expect in 3
     # years based on Richie's thesis which says "Over the 306.4 live days of
     # the D2O phase we expect a total of 192.4 events within the acrylic vessel
     # and 504.5 events within the PSUP.
-    E_hist = E_hist*1000/total_events
-    y = np.concatenate([[0],np.repeat(E_hist,2),[0]])
-    E1_hist, _ = np.histogram(E1, bins=bins)
-    E1_hist = E1_hist*1000/total_events
-    y1 = np.concatenate([[0],np.repeat(E1_hist,2),[0]])
-    E2_hist, _ = np.histogram(E2, bins=bins)
-    E2_hist = E2_hist*1000/total_events
-    y2 = np.concatenate([[0],np.repeat(E2_hist,2),[0]])
-    E3_hist, _ = np.histogram(E3, bins=bins)
-    E3_hist = E3_hist*1000/total_events
-    y3 = np.concatenate([[0],np.repeat(E3_hist,2),[0]])
+    KE_hist = KE_hist*230/total_events
+    y = np.concatenate([[0],np.repeat(KE_hist,2),[0]])
+    KE1_hist, _ = np.histogram(KE1, bins=bins)
+    KE1_hist = KE1_hist*230/total_events
+    y1 = np.concatenate([[0],np.repeat(KE1_hist,2),[0]])
+    KE2_hist, _ = np.histogram(KE2, bins=bins)
+    KE2_hist = KE2_hist*230/total_events
+    y2 = np.concatenate([[0],np.repeat(KE2_hist,2),[0]])
+    KE3_hist, _ = np.histogram(KE3, bins=bins)
+    KE3_hist = KE3_hist*230/total_events
+    y3 = np.concatenate([[0],np.repeat(KE3_hist,2),[0]])
     plt.plot(x, y, label="All events")
     plt.plot(x, y1, label="n")
     plt.plot(x, y2, label="n + nrings")
     plt.plot(x, y3, label="n + nrings + same flavor")
     plt.xlabel("Energy (GeV)")
-    plt.ylabel("Events/year")
+    plt.ylabel(r"Expected Event Rate (year$^{-1}$)")
     plt.gca().set_xscale("log")
     plt.gca().xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(tick_formatter))
     plt.xlim((0.02,bins[-1]))
     plt.ylim((0,None))
     plt.legend()
-    plt.title("Approximate Total Energy")
+    plt.title("Approximate Visible Energy")
 
     plt.figure()
     plt.hist(r, bins=np.linspace(0,8,100), histtype='step')