update script to analyze GENIE Monte Carlo

- plot energy for just the neutron follower cut
- update pions to produce 2 muon like rings
- add a check to see if the pion is above the cerenkov threshold

commit 924d4277be5ff67938247cd1babc83eaf85e7852
Author: tlatorre <tlatorre@uchicago.edu>
Date:   Wed May 23 12:15:32 2018 -0400

    update

1 files changed, 32 insertions, 9 deletions

diff --git a/analyze_genie_mc.py b/analyze_genie_mc.py
index 86691c9..cdc839c 100755
--- a/analyze_genie_mc.py
+++ b/analyze_genie_mc.py
@@ -177,9 +177,21 @@ if __name__ == '__main__':
                     e_like_rings += 1
                 elif abs(pdg) == 211:
                     # pi+/-
-                    nrings += 2
-                    u_like_rings += 1
-                    e_like_rings += 1
+                    # momentum of ith particle in hadronic system
+                    p = np.sqrt(event.pxf[i]**2 + event.pyf[i]**2 + event.pzf[i]**2)
+                    # velocity of ith particle (in units of c)
+                    # FIXME: is energy total energy or kinetic energy?
+                    v = p/event.Ef[i]
+                    if v > 1/INDEX_HEAVY_WATER:
+                        # if the pion is above threshold, we assume that it
+                        # produces 2 muon like rings
+                        nrings += 2
+                        u_like_rings += 2
+                    else:
+                        # if the pion is not above threshold, we assume that it
+                        # produces 1 muon like ring
+                        nrings += 1
+                        u_like_rings += 1
                 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)
@@ -221,13 +233,16 @@ if __name__ == '__main__':
     total_e_like_rings = np.array(total_e_like_rings)
     total_u_like_rings = np.array(total_u_like_rings)
 
-    cut1 = (total_neutrons_detected == 0) & (total_nrings >= 2)
-    cut2 = (total_neutrons_detected == 0) & (total_nrings >= 2) & ((total_e_like_rings == 0) | (total_u_like_rings == 0))
+    cut1 = (total_neutrons_detected == 0)
+    cut2 = (total_neutrons_detected == 0) & (total_nrings >= 2)
+    cut3 = (total_neutrons_detected == 0) & (total_nrings >= 2) & ((total_e_like_rings == 0) | (total_u_like_rings == 0))
 
     El1 = El[cut1]
     El2 = El[cut2]
+    El3 = El[cut3]
     E1 = E[cut1]
     E2 = E[cut2]
+    E3 = E[cut3]
 
     plt.figure()
     bincenters = (bins[1:] + bins[:-1])/2
@@ -246,9 +261,13 @@ if __name__ == '__main__':
     El2_hist, _ = np.histogram(El2, bins=bins)
     El2_hist = El2_hist*1000/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
+    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 + nrings")
-    plt.step(x, y2, where='mid', label="n + nrings + same flavor")
+    plt.step(x, y1, where='mid', label="n")
+    plt.step(x, y2, where='mid', label="n + nrings")
+    plt.step(x, y3, where='mid', label="n + nrings + same flavor")
     plt.xlabel("Energy (GeV)")
     plt.ylabel("Events/year")
     plt.gca().set_xscale("log")
@@ -273,9 +292,13 @@ if __name__ == '__main__':
     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]])
     plt.plot(x, y, label="All events")
-    plt.plot(x, y1, label="n + nrings")
-    plt.plot(x, y2, label="n + nrings + same flavor")
+    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.gca().set_xscale("log")