diff options
| -rwxr-xr-x | analyze_genie_mc.py | 10 |
1 files changed, 9 insertions, 1 deletions
diff --git a/analyze_genie_mc.py b/analyze_genie_mc.py index 8c006a3..9739c65 100755 --- a/analyze_genie_mc.py +++ b/analyze_genie_mc.py @@ -20,6 +20,10 @@ NEUTRON_DETECTION_EFFICIENCY = 0.5 # FIXME: What is the index for D2O? INDEX_HEAVY_WATER = 1.3 +# fractional energy resolution +# from Richie's thesis page 134 +ENERGY_RESOLUTION = 0.05 + def pdg_code_to_string(pdg): A = int(("%010i" % event.tgt)[-4:-1]) Z = int(("%010i" % event.tgt)[-7:-4]) @@ -250,7 +254,7 @@ if __name__ == '__main__': total_u_like_rings.append(u_like_rings) El.append(event.El) E.append(event.calresp0) - KE.append(ke) + KE.append(ke + np.random.randn()*ke*ENERGY_RESOLUTION) 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)): @@ -319,6 +323,7 @@ if __name__ == '__main__': plt.figure() KE_hist, _ = np.histogram(KE, bins=bins) + KE_signal, _ = np.histogram(np.random.randn(1000)*1.0*ENERGY_RESOLUTION + 1.0, 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 @@ -335,10 +340,13 @@ if __name__ == '__main__': KE3_hist, _ = np.histogram(KE3, bins=bins) KE3_hist = KE3_hist*230/total_events y3 = np.concatenate([[0],np.repeat(KE3_hist,2),[0]]) + KE_signal = KE_signal*10/np.sum(KE_signal) + y4 = np.concatenate([[0],np.repeat(KE_signal,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.plot(x, y4, label="1 GeV signal") plt.xlabel("Energy (GeV)") plt.ylabel(r"Expected Event Rate (year$^{-1}$)") plt.gca().set_xscale("log") |