diff options
Diffstat (limited to 'analyze_genie_mc.py')
| -rwxr-xr-x | analyze_genie_mc.py | 373 |
1 files changed, 0 insertions, 373 deletions
diff --git a/analyze_genie_mc.py b/analyze_genie_mc.py deleted file mode 100755 index 9739c65..0000000 --- a/analyze_genie_mc.py +++ /dev/null @@ -1,373 +0,0 @@ -#!/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() |