aboutsummaryrefslogtreecommitdiff
path: root/utils/analyze-genie-mc
diff options
context:
space:
mode:
authortlatorre <tlatorre@uchicago.edu>2018-10-03 09:52:19 -0500
committertlatorre <tlatorre@uchicago.edu>2018-10-03 09:52:19 -0500
commit86e89541032a0a345ae70d4b31aec474e4899890 (patch)
tree8257c8f1a7ca847a5f8d4cfcdee055ed02d1049a /utils/analyze-genie-mc
parent405f0d448eb15b4f686c1f7ee47794d82eaaf62c (diff)
downloadsddm-86e89541032a0a345ae70d4b31aec474e4899890.tar.gz
sddm-86e89541032a0a345ae70d4b31aec474e4899890.tar.bz2
sddm-86e89541032a0a345ae70d4b31aec474e4899890.zip
move python scripts into utils/ directory
Diffstat (limited to 'utils/analyze-genie-mc')
-rwxr-xr-xutils/analyze-genie-mc373
1 files changed, 373 insertions, 0 deletions
diff --git a/utils/analyze-genie-mc b/utils/analyze-genie-mc
new file mode 100755
index 0000000..9739c65
--- /dev/null
+++ b/utils/analyze-genie-mc
@@ -0,0 +1,373 @@
+#!/usr/bin/env python3
+import ROOT
+import numpy as np
+
+# on retina screens, the default plots are way too small
+# by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1
+# Qt5 will scale everything using the dpi in ~/.Xresources
+import matplotlib
+matplotlib.use("Qt5Agg")
+
+def tick_formatter(x, pos):
+ if x < 1:
+ return '%.1f' % x
+ else:
+ return '%.0f' % x
+
+# FIXME: What is this for the salt phase?
+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])
+ if Z == 1:
+ atom = 'H'
+ elif Z == 8:
+ atom = 'O'
+ elif Z == 6:
+ atom = 'C'
+ else:
+ raise NotImplementedError("unknown atom %i" % Z)
+ return '%i%s' % (A,atom)
+
+def get_reaction(event):
+ reactants = []
+ products = []
+ if event.neu == 12:
+ reactants.append('ve')
+ elif event.neu == -12:
+ reactants.append('vebar')
+ elif event.neu == 14:
+ reactants.append('vu')
+ elif event.neu == -14:
+ reactants.append('vubar')
+ elif event.neu == 16:
+ reactants.append('vt')
+ elif event.neu == -16:
+ reactants.append('vtbar')
+
+ if event.hitnuc == 2212:
+ reactants.append('p')
+ elif event.hitnuc == 2112:
+ reactants.append('n')
+ elif event.hitnuc == 0:
+ reactants.append(pdg_code_to_string(event.tgt))
+ else:
+ print("unknown nucleon %i" % event.hitnuc)
+
+ if event.cc:
+ if event.neu == 12:
+ products.append('e-')
+ elif event.neu == -12:
+ products.append('e+')
+ elif event.neu == 14:
+ products.append('u-')
+ elif event.neu == -14:
+ products.append('u+')
+ elif event.neu == 16:
+ products.append('t-')
+ elif event.neu == -16:
+ products.append('t+')
+ elif event.nc:
+ if event.neu == 12:
+ products.append('ve')
+ elif event.neu == -12:
+ products.append('vebar')
+ elif event.neu == 14:
+ products.append('vu')
+ elif event.neu == -14:
+ products.append('vubar')
+ elif event.neu == 16:
+ products.append('vt')
+ elif event.neu == -16:
+ products.append('vtbar')
+ else:
+ products.append("???")
+
+ for pdg in event.pdgf:
+ if pdg == 2112:
+ products.append('n')
+ elif abs(pdg) == 11:
+ # e- or e+
+ if pdg == 11:
+ products.append('e-')
+ else:
+ products.append('e+')
+ elif pdg == 22:
+ # gamma
+ products.append('gamma')
+ elif pdg == 111:
+ # pi0
+ products.append('pi0')
+ elif abs(pdg) == 211:
+ # pi+/-
+ if pdg == 211:
+ products.append('pi+')
+ else:
+ products.append('pi-')
+ elif abs(pdg) == 311:
+ if pdg == 311:
+ products.append('K0')
+ else:
+ products.append('K0bar')
+ elif abs(pdg) == 321:
+ # K+/-
+ if pdg == 321:
+ products.append('K+')
+ else:
+ products.append('K-')
+ elif abs(pdg) == 3222:
+ products.append('Sigma+')
+ elif abs(pdg) == 3112:
+ products.append('Sigma-')
+ elif abs(pdg) == 3122:
+ products.append('Delta')
+ elif pdg == 2212:
+ products.append('p')
+ elif int(("%010i" % abs(pdg))[0]) == 1:
+ products.append(pdg_code_to_string(pdg))
+ else:
+ print("unknown pdg code %i" % pdg)
+
+ return ' + '.join(reactants) + ' -> ' + ' + '.join(products)
+
+if __name__ == '__main__':
+ import argparse
+ import matplotlib.pyplot as plt
+ from collections import Counter
+
+ parser = argparse.ArgumentParser("script to analyze GENIE 'ntuple' ROOT files")
+ parser.add_argument("filenames", nargs='+', help="GENIE ROOT files")
+ args = parser.parse_args()
+
+ bins = np.logspace(-1,2,100)
+
+ El = []
+ total_neutrons = []
+ total_neutrons_detected = []
+ E = []
+ KE = []
+ r = []
+ total_nrings = []
+ total_e_like_rings = []
+ total_u_like_rings = []
+ reactions = Counter()
+ for filename in args.filenames:
+ print("analyzing %s" % filename)
+ f = ROOT.TFile(filename)
+ T = f.Get("gst")
+
+ for event in T:
+ neutrons = 0
+ nrings = 0
+ e_like_rings = 0
+ u_like_rings = 0
+ ke = 0
+
+ if event.cc:
+ if abs(event.neu) == 12:
+ e_like_rings = 1
+ else:
+ u_like_rings = 1
+ nrings = 1
+ ke += event.El
+ elif event.nc:
+ pass
+ else:
+ print("event is not cc or nc!")
+ continue
+
+ for i, pdg in enumerate(event.pdgf):
+ if pdg == 2112:
+ neutrons += 1
+ elif abs(pdg) == 11:
+ # e- or e+
+ if event.Ef[i] > 0.1:
+ # for now assume we only count rings from electrons
+ # with > 100 MeV
+ nrings += 1
+ e_like_rings += 1
+ ke += event.Ef[i]
+ elif pdg == 22:
+ # gamma
+ if event.Ef[i] > 0.1:
+ # for now assume we only count rings from gammas with >
+ # 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
+ 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
+ # 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)
+ # 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:
+ # 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
+ pass
+ else:
+ print("unknown pdg code %i" % pdg)
+ total_neutrons.append(neutrons)
+ total_neutrons_detected.append(np.random.binomial(neutrons,NEUTRON_DETECTION_EFFICIENCY))
+ total_nrings.append(nrings)
+ total_e_like_rings.append(e_like_rings)
+ total_u_like_rings.append(u_like_rings)
+ El.append(event.El)
+ E.append(event.calresp0)
+ 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)):
+ reactions.update([get_reaction(event)])
+
+ total = sum(reactions.values())
+ for reaction, count in reactions.most_common(10):
+ print("%.0f%% %s" % (count*100/total, reaction))
+
+ El = np.array(El)
+ 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)
+ total_u_like_rings = np.array(total_u_like_rings)
+
+ 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]
+ KE1 = KE[cut1]
+ KE2 = KE[cut2]
+ KE3 = KE[cut3]
+
+ plt.figure()
+ bincenters = (bins[1:] + bins[:-1])/2
+ x = np.repeat(bins,2)
+ El_hist, _ = np.histogram(El, bins=bins)
+ total_events = El_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.
+ 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*230/total_events
+ y1 = np.concatenate([[0],np.repeat(El1_hist,2),[0]])
+ El2_hist, _ = np.histogram(El2, bins=bins)
+ 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*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")
+ 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")
+ 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("Primary Lepton Energy")
+
+ 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
+ # the D2O phase we expect a total of 192.4 events within the acrylic vessel
+ # and 504.5 events within the PSUP.
+ 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]])
+ 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")
+ 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 Visible Energy")
+
+ plt.figure()
+ plt.hist(r, bins=np.linspace(0,8,100), histtype='step')
+ plt.xlabel("R (m)")
+ plt.title("Radius of Events")
+
+ plt.figure()
+ plt.hist(total_neutrons, bins=np.arange(11)-0.5, histtype='step')
+ plt.xlabel("Number of Neutrons")
+ plt.title("Number of Neutrons")
+
+ plt.figure()
+ plt.hist(total_nrings, bins=np.arange(11)-0.5, histtype='step')
+ plt.xlabel("Number of Rings")
+ plt.title("Number of Rings (approximate)")
+ plt.show()