diff options
Diffstat (limited to 'analyze_genie_mc.py')
| -rwxr-xr-x | analyze_genie_mc.py | 41 |
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") |