From 86e89541032a0a345ae70d4b31aec474e4899890 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Wed, 3 Oct 2018 09:52:19 -0500 Subject: move python scripts into utils/ directory --- analyze-genie-mc | 373 ---------------------------------------- calculate_limits.py | 419 --------------------------------------------- convert-genie-to-gst | 21 --- convert-ratdb | 43 ----- index.py | 85 --------- src/plot.py | 78 --------- utils/analyze-genie-mc | 373 ++++++++++++++++++++++++++++++++++++++++ utils/calculate_limits.py | 419 +++++++++++++++++++++++++++++++++++++++++++++ utils/convert-genie-to-gst | 21 +++ utils/convert-ratdb | 43 +++++ utils/index.py | 85 +++++++++ utils/plot.py | 78 +++++++++ 12 files changed, 1019 insertions(+), 1019 deletions(-) delete mode 100755 analyze-genie-mc delete mode 100755 calculate_limits.py delete mode 100755 convert-genie-to-gst delete mode 100755 convert-ratdb delete mode 100644 index.py delete mode 100755 src/plot.py create mode 100755 utils/analyze-genie-mc create mode 100755 utils/calculate_limits.py create mode 100755 utils/convert-genie-to-gst create mode 100755 utils/convert-ratdb create mode 100644 utils/index.py create mode 100755 utils/plot.py diff --git a/analyze-genie-mc b/analyze-genie-mc deleted file mode 100755 index 9739c65..0000000 --- a/analyze-genie-mc +++ /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() diff --git a/calculate_limits.py b/calculate_limits.py deleted file mode 100755 index f297a69..0000000 --- a/calculate_limits.py +++ /dev/null @@ -1,419 +0,0 @@ -#!/usr/bin/env python3 -import numpy as np -from scipy.integrate import quad -from scipy.stats import expon -from scipy.special import spherical_jn, erf -from numpy import pi -from functools import lru_cache - -# 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") - -# speed of light (exact) -SPEED_OF_LIGHT = 299792458 # m/s - -# depth of the SNO detector (m) -# currently just converted 6800 feet -> meters -h = 2072.0 - -# radius of earth (m) -# from the "Earth Fact Sheet" at https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html -# we use the volumetric mean radius here -R = 6.371e6 - -# Approximate dark matter velocity in m/s. The true distribution is expected to -# be a Maxwell Boltzmann distribution which is modulated annually by the -# earth's rotation around the sun, but we just assume a single constant -# velocity here. From Lewin and Smith Appendix B -DM_VELOCITY = 230e3 - -# Approximate earth velocity in the galactic rest frame (?) -# from Lewin and Smith Equation 3.6 -EARTH_VELOCITY = 244e3 - -# Approximate dark matter density in GeV/m^3. From Tom Caldwell's thesis. -# from Lewin and Smith page 91 -DM_DENSITY = 0.4e6 - -# Number density of scatterers in the Earth. -# -# FIXME: Currently just set to the number density of atoms in water. Need to -# update this for rock, and in fact this will change near the detector since -# there is water outside the AV. -DENSITY_WATER = 1e3 # In kg/m^3 - -# mean density of norite rock which is the rock surrounding SNO -# probably conservative -NORITE_DENSITY = 3e3 # In kg/m^3 - -# atomic masses for various elements -# from https://www.angelo.edu/faculty/kboudrea/periodic/structure_mass.htm -element_mass = { - 'H':1.0, - 'C':12.011, - 'O':15.9994, - 'Na':22.98977, - 'Mg':24.305, - 'Al':26.98154, - 'Si':28.0855, - 'K':39.0983, - 'Ca':40.08, - 'Mn':54.9380, - 'Fe':55.847, - 'Ti':47.90 -} - -# composition and mean density of norite rock which is the rock surrounding SNO -# the composition is from Table 3.2 in the SNOLAB User's handbook -water = {'composition': - {'H':20, - 'O':80}, - 'density':1e3} - -# composition and mean density of the mantle -# from www.knowledgedoor.com/2/elements_handbook/element_abundances_in_the_earth_s_mantle.html -# density from hyperphysics.phy-astr.gsu.edu/hbase/Geophys/earthstruct.html -mantle = {'composition': - {'O':44.33, - 'Mg':22.17, - 'Si':21.22, - 'Fe':6.3}, - 'density':4.4e3} - -# composition and mean density of norite rock which is the rock surrounding SNO -# the composition is from Table 3.2 in the SNOLAB User's handbook -norite = {'composition': - {'H':0.15, - 'C':0.04, - 'O':46.0, - 'Na':2.2, - 'Mg':3.3, - 'Al':9.0, - 'Si':26.2, - 'K': 1.2, - 'Ca':5.2, - 'Mn':0.1, - 'Fe':6.2, - 'Ti':0.5}, - 'density':3e3} - -# Fiducial volume (m) -FIDUCIAL_RADIUS = 5 - -# Fiducial volume (m^3) -FIDUCIAL_VOLUME = 4*pi*FIDUCIAL_RADIUS**3/3 - -# proton mass from the PDG (2018) -PROTON_MASS = 0.938 # GeV - -# proton mass from the PDG (2018) -ATOMIC_MASS_UNIT = 0.931 # GeV - -# mass of Xenon in atomic units -XENON_MASS = 131.293 - -# mass of Neon in atomic units -NEON_MASS = 20.18 - -# mass of argon in atomic units -ARGON_MASS = 39.948 - -# mass of germanium in atomic units -GERMANIUM_MASS = 72.64 - -# mass of tungsten in atomic units -TUNGSTEN_MASS = 183.84 - -# mass of oxygen in atomic units -OXYGEN_MASS = 15.999 - -# mass of silicon in atomic units -SILICON_MASS = 28.0855 - -# mass of iron in atomic units -IRON_MASS = 55.845 - -# mass of magnesium in atomic units -MAGNESIUM_MASS = 24.305 - -# galactic escape velocity (m/s) -# from Tom Caldwell's thesis page 25 -ESCAPE_VELOCITY = 244e3 - -# conversion constant from PDG -HBARC = 197.326978812e-3 # GeV fm - -# Avogadros number (kg^-1) -N0 = 6.02214085774e26 - -def get_probability(r, l): - """ - Returns the probability of a dark photon decaying in the SNO detector from - a dark matter decay distributed uniformly in the Earth. Assumes that the - depth of SNO is much larger than the dimensions of the SNO detector. - """ - if r <= h: - theta_min = 0 - elif r <= 2*R - h: - theta_min = pi - np.arccos((h**2 + r**2 - 2*R*h)/(2*r*(R-h))) - else: - return 0 - return (1 + np.cos(theta_min))*np.exp(-r/l)/(2*l) - -@lru_cache() -def get_probability2(l): - p, err = quad(get_probability, 0, min(2*R-h,10*l), args=(l,), epsabs=0, epsrel=1e-7, limit=1000) - - return p - -def get_event_rate(m, cs0, l, A): - """ - Returns the event rate of leptons produced from self-destructing dark - matter in the SNO detector for a given dark matter mass m, a cross section - cs, and a mediator decay length l. - """ - # For now we assume the event rate is constant throughout the earth, so we - # are implicitly assuming that the cross section is pretty small. - flux = DM_VELOCITY*DM_DENSITY/m - - #p, err = quad(get_probability, 0, min(2*R-h,10*l), args=(l,), epsabs=0, epsrel=1e-7, limit=1000) - p = get_probability2(l) - - cs = get_cross_section(cs0, m, A) - - # FIXME: factor of 2 because the DM particle decays into two mediators? - return p*cs*(N0/A)*flux*FIDUCIAL_VOLUME - -def get_event_rate_sno(m, cs0, l, composition): - """ - Returns the event rate of leptons produced from self-destructing dark - matter in the SNO detector for a given dark matter mass m, a cross section - cs, and a mediator decay length l. - """ - rate = 0.0 - - for element, mass_fraction in composition['composition'].items(): - rate += mass_fraction/100*get_event_rate(m,cs0,l,element_mass[element]) - - return rate*composition['density'] - -def get_nuclear_form_factor(A, e): - """ - Returns the nuclear form factor for a WIMP-nucleus interaction. - - From Tom Caldwell's thesis page 24. - - Also used Mark Pepin's thesis page 50 - """ - # mass of xenon nucleus - mn = A*ATOMIC_MASS_UNIT - - # calculate approximate size of radius - s = 0.9 # fm - a = 0.52 # fm - c = 1.23*A**(1/3) - 0.60 # fm - r1 = np.sqrt(c**2 + (7/3)*pi**2*a**2 - 5*s**2) - - q = np.sqrt(2*mn*e) - - if q*r1/HBARC < 1e-10: - return 1.0 - - # Helm form factor - # from Mark Pepin's thesis page 50 - f = 3*spherical_jn(1,q*r1/HBARC)*np.exp(-(q*s/HBARC)**2/2)/(q*r1/HBARC) - - return f - -def get_cross_section(cs0, m, A): - """ - Returns the WIMP cross section from the target-independent WIMP-nucleon - cross section cs0 at zero momentum transfer. - - From Tom Caldwell's thesis page 21. - """ - # mass of xenon nucleus - mn = A*ATOMIC_MASS_UNIT - - # reduced mass of the nucleus and the WIMP - mr = (m*mn)/(m + mn) - - # reduced mass of the proton and the WIMP - mp = (m*PROTON_MASS)/(m + PROTON_MASS) - - return cs0*A**2*mr**2/mp**2 - -def get_differential_event_rate_xenon(e, m, cs0, A): - """ - Returns the event rate of WIMP scattering in the Xenon 100T detector. - """ - # mass of nucleus - mn = A*ATOMIC_MASS_UNIT - - # reduced mass of the nucleus and the WIMP - mr = (m*mn)/(m + mn) - - cs = get_cross_section(cs0, m, A) - - v0 = DM_VELOCITY - - vesc = ESCAPE_VELOCITY - - # earth's velocity through the galaxy - ve = EARTH_VELOCITY - - # minimum wimp velocity needed to produce a recoil of energy e - vmin = np.sqrt(mn*e/2)*(mn+m)*SPEED_OF_LIGHT/(mn*m) - - f = get_nuclear_form_factor(A, e) - - x = vmin/v0 - y = ve/v0 - z = vesc/v0 - - # Equation 3.49 in Mark Pepin's thesis - k0 = (pi*v0**2)**(3/2) - - # Equation 3.49 in Mark Pepin's thesis - k1 = k0*(erf(z)-(2/np.sqrt(pi))*z*np.exp(-z**2)) - - # From Mark Pepin's CDMS thesis page 59 - if x <= z - y: - I = (k0/k1)*(DM_DENSITY/(2*ve))*(erf(x+y) - erf(x-y) - (4/np.sqrt(pi))*y*np.exp(-z**2)) - elif x <= y + z: - I = (k0/k1)*(DM_DENSITY/(2*ve))*(erf(z) - erf(x-y) - (2/np.sqrt(pi))*(y+z-x)*np.exp(-z**2)) - else: - return 0 - - return (N0/A)*(mn/(2*m*mr**2))*cs*f**2*I*SPEED_OF_LIGHT**2 - -def get_event_rate_xenon(m, cs, A, threshold): - """ - Returns the event rate of WIMP scattering in a dark matter detector using - an element with atomic mass A. Rate is in Hz kg^-1. - """ - # mass of nucleus - mn = A*ATOMIC_MASS_UNIT - - # reduced mass of the nucleus and the WIMP - mr = (m*mn)/(m + mn) - - vesc = ESCAPE_VELOCITY - - # earth's velocity through the galaxy - ve = EARTH_VELOCITY - - # max recoil they looked for was 40 keV - emax = 2*mr**2*((ve+vesc)/SPEED_OF_LIGHT)**2/mn - - rate, err = quad(get_differential_event_rate_xenon, threshold, emax, args=(m,cs,A), epsabs=0, epsrel=1e-2, limit=1000) - - return rate - -if __name__ == '__main__': - import matplotlib.pyplot as plt - - ls = np.logspace(-1,8,1000) - - cs = 1e-50 # cm^2 - - # FIXME: should use water density for L < 1 m, silicon density for L ~ 1 - # km, and iron for L >> 1 km - rate = np.array([get_event_rate_sno(1, cs*1e-4, l, water) for l in ls]) - - colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] - - plt.figure() - plt.subplot(111) - plt.plot(ls, rate/np.max(rate),color=colors[0]) - plt.xlabel("Mediator Decay Length (m)") - plt.ylabel("Event Rate (arbitrary units)") - plt.axvline(x=FIDUCIAL_RADIUS, color=colors[1], ls='--', label="SNO radius") - plt.axvline(x=h, ls='--', color=colors[2], label="Depth of SNO") - plt.axvline(x=R, ls='--', color=colors[3], label="Earth radius") - plt.gca().set_xscale('log') - plt.gca().set_xlim((ls[0], ls[-1])) - plt.title("Event Rate of Self-Destructing Dark Matter in SNO") - plt.legend() - - # threshold is ~5 keVr - xenon_100t_threshold = 5e-6 - - # threshold is ~1 keVr - cdms_threshold = 1e-6 - - # threshold is ~100 eV - # FIXME: is this correct? - cresst_threshold = 1e-7 - - ms = np.logspace(-2,3,200) - cs0s = np.logspace(-50,-40,200) - - mm, cs0cs0 = np.meshgrid(ms, cs0s) - - rate1 = np.empty(mm.shape) - rate2 = np.empty(mm.shape) - rate3 = np.empty(mm.shape) - rate4 = np.empty(mm.shape) - rate5 = np.empty(mm.shape) - rate6 = np.empty(mm.shape) - for i in range(mm.shape[0]): - print("\r%i/%i" % (i+1,mm.shape[0]),end='') - for j in range(mm.shape[1]): - rate1[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, XENON_MASS, xenon_100t_threshold) - rate2[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1.0, water) - rate3[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, GERMANIUM_MASS, cdms_threshold) - rate4[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, TUNGSTEN_MASS, cresst_threshold) - rate5[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1e3, norite) - rate6[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1e6, mantle) - print() - - # Fiducial volume of the Xenon1T detector is 1042 +/- 12 kg - # from arxiv:1705.06655 - xenon_100t_fiducial_volume = 1042 # kg - - # Livetime of the Xenon1T results - # from arxiv:1705.06655 - xenon_100t_livetime = 34.2 # days - - plt.figure() - plt.subplot(111) - plt.gca().set_xscale('log') - plt.gca().set_yscale('log') - CS1 = plt.contour(mm,cs0cs0,rate1*3600*24*xenon_100t_fiducial_volume*xenon_100t_livetime,[10.0], colors=[colors[0]]) - CS2 = plt.contour(mm,cs0cs0,rate2*3600*24*668.8,[10.0],colors=[colors[1]]) - CS3 = plt.contour(mm,cs0cs0,rate3*3600*24*70.10,[10.0],colors=[colors[2]]) - # FIXME: I used 2.39 kg day because that's what CRESST-3 reports in their paper - # but! I only use Tungsten here so do I need to multiply by the mass fraction of tungsten? - CS4 = plt.contour(mm,cs0cs0,rate4*3600*24*2.39,[10.0],colors=[colors[3]]) - CS5 = plt.contour(mm,cs0cs0,rate5*3600*24*668.8,[10.0],colors=[colors[1]], linestyles=['dashed']) - CS6 = plt.contour(mm,cs0cs0,rate6*3600*24*668.8,[10.0],colors=[colors[1]], linestyles=['dotted']) - plt.clabel(CS1, inline=1, fmt="XENON1T", fontsize=10, use_clabeltext=True) - plt.clabel(CS2, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1 m)", fontsize=10) - plt.clabel(CS3, inline=1, fmt="CDMSLite", fontsize=10) - plt.clabel(CS4, inline=1, fmt="CRESST-3", fontsize=10) - plt.clabel(CS5, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1 km)", fontsize=10) - plt.clabel(CS6, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1000 km)", fontsize=10) - plt.xlabel(r"$m_\chi$ (GeV)") - plt.ylabel(r"WIMP-nucleon scattering cross section ($\mathrm{cm}^2$)") - plt.title("Self-Destructing Dark Matter Limits") - - x = np.linspace(0,300e-6,1000) - - # reproducing Figure 2.1 in Tom Caldwell's thesis - plt.figure() - plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(XENON_MASS, x)**2,x)), label="Xe") - plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(NEON_MASS, x)**2,x)), label="Ne") - plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(ARGON_MASS, x)**2,x)), label="Ar") - plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(GERMANIUM_MASS, x)**2,x)), label="Ge") - plt.xlabel("Recoil Energy (keV)") - plt.ylabel(r"$F^2(E)$") - plt.gca().set_yscale('log') - plt.legend() - plt.gca().set_ylim((1e-5,1)) - plt.gca().set_xlim((0,300)) - plt.title("Nuclear Form Factors for various elements") - plt.show() diff --git a/convert-genie-to-gst b/convert-genie-to-gst deleted file mode 100755 index b2861e5..0000000 --- a/convert-genie-to-gst +++ /dev/null @@ -1,21 +0,0 @@ -#!/usr/bin/env python -from subprocess import check_call -from os.path import split, splitext, join -import os - -if __name__ == '__main__': - import argparse - - parser = argparse.ArgumentParser("script to convert full GENIE root files to the reduced gst ROOT format") - parser.add_argument("filenames", nargs="+", help="GENIE root files") - parser.add_argument("--dest", required=True, help="destination directory") - args = parser.parse_args() - - for filename in args.filenames: - head, tail = split(filename) - root, ext = splitext(tail) - output = join(args.dest, root) + ".ntuple.root" - cmd = ["gntpc","-f","gst","-i",filename,"-o",output] - print(" ".join(cmd)) - with open(os.devnull,"w") as devnull: - check_call(cmd, stdout=devnull, stderr=devnull) diff --git a/convert-ratdb b/convert-ratdb deleted file mode 100755 index 5dd6ee4..0000000 --- a/convert-ratdb +++ /dev/null @@ -1,43 +0,0 @@ -#!/usr/bin/env python -""" -Script to convert the PMT database in RAT to a CSV file (with spaces instead of commas). -""" -from __future__ import print_function, division -import yaml - -if __name__ == '__main__': - import argparse - import numpy as np - - parser = argparse.ArgumentParser("convert a RATDB file -> CSV file") - parser.add_argument("filename", help="RATDB filename") - parser.add_argument("-o", "--output", help="output filename", required=True) - args = parser.parse_args() - - s = '' - with open(args.filename) as f: - for line in f: - if line.startswith('//'): - continue - s += line - - data = yaml.load(s) - - csv = np.zeros((len(data['x']),7),dtype=np.float64) - csv[:,0] = data['x'] - csv[:,1] = data['y'] - csv[:,2] = data['z'] - csv[:,3] = data['u'] - csv[:,4] = data['v'] - csv[:,5] = data['w'] - csv[:,6] = data['pmt_type'] - - # reverse PMT normals - csv[:,3] = -csv[:,3] - csv[:,4] = -csv[:,4] - csv[:,5] = -csv[:,5] - - header = "Format: x, y, z, u, v, w, pmt_type" - - np.savetxt(args.output, csv, fmt=['%10.2f']*6 + ['%i'], - header=header) diff --git a/index.py b/index.py deleted file mode 100644 index 4edcd59..0000000 --- a/index.py +++ /dev/null @@ -1,85 +0,0 @@ -#!/usr/bin/env python -""" -This is a script to calculate the index of refraction of water as a function of -wavelength, temperature, and pressure and fit it to a linear approximation in a -wavelength region where the PMTs are sensitive to. The data comes from [1]. - -[1] "Refractive index of Water and Steam as Function of Wavelength, Temperature, and Density". Schiebener et al. 1989. -""" -from __future__ import print_function, division -import numpy as np - -A0 = 0.243905091 -A1 = 9.53518094e-3 -A2 = -3.64358110e-3 -A3 = 2.65666426e-4 -A4 = 1.59189325e-3 -A5 = 2.45733798e-3 -A6 = 0.897478251 -A7 = -1.63066183e-2 -UV = 0.2292020 -IR = 5.432937 - -def get_index(p, wavelength, T): - """ - Returns the index of pure water for a given density, wavelength, and - temperature. The density should be in units of kg/m^3, the wavelength in nm, - and the temperature in Celsius. - """ - # normalize the density, temperature, and pressure - p = p/1000.0 - wavelength = wavelength/589.0 - T = (T+273.15)/273.15 - - # first we compute the right hand side of Equation 7 - c = A0 + A1*p + A2*T + A3*wavelength**2*T + A4/wavelength**2 + A5/(wavelength**2-UV**2) + A6/(wavelength**2-IR**2) + A7*p**2 - c *= p - - return np.sqrt((2*c+1)/(1-c)) - -def approximate_index(x, a, b): - """ - Returns the index of refraction using the approximation: - - n = a + b/x**2 - """ - return a + b/x**2 - -def chi2(pars, x, data): - """ - Returns the squares of the difference between the data and the model for - the inverse of the index of refraction, i.e. - - 1/n = lambda**2/(a*lambda**2+b) - """ - a, b = pars - return np.sum((approximate_index(x,a,b)-data)**2) - -if __name__ == '__main__': - import matplotlib.pyplot as plt - import argparse - from scipy.optimize import fmin - - parser = argparse.ArgumentParser("fit for an approximation to the inverse of the index of refraction of water") - parser.add_argument("--temperature", type=float, default=10.0) - parser.add_argument("--density", type=float, default=1000.0) - parser.add_argument("--low", type=float, default=300.0) - parser.add_argument("--high", type=float, default=600.0) - args = parser.parse_args() - - x = np.linspace(args.low,args.high,1000) - n = get_index(args.density,x,args.temperature) - xmin = fmin(chi2,[1,1],args=(x,1/n)) - - print(xmin) - - a, b = xmin - y = approximate_index(x,a,b) - - plt.plot(x,1/n,label="data") - plt.plot(x,y,label="fit") - plt.xlabel("Wavelength (nm)") - plt.ylabel("1/n") - plt.title("Index of Refraction of Water at 273.15 K") - plt.legend() - plt.show() diff --git a/src/plot.py b/src/plot.py deleted file mode 100755 index 6a2fe73..0000000 --- a/src/plot.py +++ /dev/null @@ -1,78 +0,0 @@ -#!/usr/bin/env python -from __future__ import print_function, division -import yaml - -SNOMAN_MASS = { - 20: 0.511, - 21: 0.511, - 22: 105.658, - 23: 105.658 -} - -if __name__ == '__main__': - import argparse - import matplotlib.pyplot as plt - import numpy as np - - parser = argparse.ArgumentParser("plot fit results") - parser.add_argument("filenames", nargs='+', help="input files") - args = parser.parse_args() - - for filename in args.filenames: - print(filename) - with open(filename) as f: - data = yaml.load(f.read()) - - dx = [] - dy = [] - dz = [] - e = [] - for event in data['data']: - # get the particle ID - id = event['mctk'][-1]['id'] - mass = SNOMAN_MASS[id] - # for some reason it's the *second* track which seems to contain the - # initial energy - true_energy = event['mctk'][-1]['energy'] - # The MCTK bank has the particle's total energy (except for neutrons) - # so we need to convert it into kinetic energy - ke = true_energy - mass - energy = event['ev'][0]['fit'][0]['energy'] - e.append(energy-ke) - true_posx = event['mcvx'][0]['posx'] - posx = event['ev'][0]['fit'][0]['posx'] - dx.append(posx-true_posx) - true_posy = event['mcvx'][0]['posy'] - posy = event['ev'][0]['fit'][0]['posy'] - dy.append(posy-true_posy) - true_posz = event['mcvx'][0]['posz'] - posz = event['ev'][0]['fit'][0]['posz'] - dz.append(posz-true_posz) - - print("dT = %.2g +/- %.2g" % (np.mean(e), np.std(e))) - print("dx = %4.2g +/- %.2g" % (np.mean(dx), np.std(dx))) - print("dy = %4.2g +/- %.2g" % (np.mean(dy), np.std(dy))) - print("dz = %4.2g +/- %.2g" % (np.mean(dz), np.std(dz))) - - plt.figure(1) - plt.hist(e, histtype='step', label=filename) - plt.xlabel("Kinetic Energy difference (MeV)") - plt.figure(2) - plt.hist(dx, histtype='step', label=filename) - plt.xlabel("X Position difference (cm)") - plt.figure(3) - plt.hist(dy, histtype='step', label=filename) - plt.xlabel("Y Position difference (cm)") - plt.figure(4) - plt.hist(dz, histtype='step', label=filename) - plt.xlabel("Z Position difference (cm)") - - plt.figure(1) - plt.legend() - plt.figure(2) - plt.legend() - plt.figure(3) - plt.legend() - plt.figure(4) - plt.legend() - plt.show() 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() diff --git a/utils/calculate_limits.py b/utils/calculate_limits.py new file mode 100755 index 0000000..f297a69 --- /dev/null +++ b/utils/calculate_limits.py @@ -0,0 +1,419 @@ +#!/usr/bin/env python3 +import numpy as np +from scipy.integrate import quad +from scipy.stats import expon +from scipy.special import spherical_jn, erf +from numpy import pi +from functools import lru_cache + +# 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") + +# speed of light (exact) +SPEED_OF_LIGHT = 299792458 # m/s + +# depth of the SNO detector (m) +# currently just converted 6800 feet -> meters +h = 2072.0 + +# radius of earth (m) +# from the "Earth Fact Sheet" at https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html +# we use the volumetric mean radius here +R = 6.371e6 + +# Approximate dark matter velocity in m/s. The true distribution is expected to +# be a Maxwell Boltzmann distribution which is modulated annually by the +# earth's rotation around the sun, but we just assume a single constant +# velocity here. From Lewin and Smith Appendix B +DM_VELOCITY = 230e3 + +# Approximate earth velocity in the galactic rest frame (?) +# from Lewin and Smith Equation 3.6 +EARTH_VELOCITY = 244e3 + +# Approximate dark matter density in GeV/m^3. From Tom Caldwell's thesis. +# from Lewin and Smith page 91 +DM_DENSITY = 0.4e6 + +# Number density of scatterers in the Earth. +# +# FIXME: Currently just set to the number density of atoms in water. Need to +# update this for rock, and in fact this will change near the detector since +# there is water outside the AV. +DENSITY_WATER = 1e3 # In kg/m^3 + +# mean density of norite rock which is the rock surrounding SNO +# probably conservative +NORITE_DENSITY = 3e3 # In kg/m^3 + +# atomic masses for various elements +# from https://www.angelo.edu/faculty/kboudrea/periodic/structure_mass.htm +element_mass = { + 'H':1.0, + 'C':12.011, + 'O':15.9994, + 'Na':22.98977, + 'Mg':24.305, + 'Al':26.98154, + 'Si':28.0855, + 'K':39.0983, + 'Ca':40.08, + 'Mn':54.9380, + 'Fe':55.847, + 'Ti':47.90 +} + +# composition and mean density of norite rock which is the rock surrounding SNO +# the composition is from Table 3.2 in the SNOLAB User's handbook +water = {'composition': + {'H':20, + 'O':80}, + 'density':1e3} + +# composition and mean density of the mantle +# from www.knowledgedoor.com/2/elements_handbook/element_abundances_in_the_earth_s_mantle.html +# density from hyperphysics.phy-astr.gsu.edu/hbase/Geophys/earthstruct.html +mantle = {'composition': + {'O':44.33, + 'Mg':22.17, + 'Si':21.22, + 'Fe':6.3}, + 'density':4.4e3} + +# composition and mean density of norite rock which is the rock surrounding SNO +# the composition is from Table 3.2 in the SNOLAB User's handbook +norite = {'composition': + {'H':0.15, + 'C':0.04, + 'O':46.0, + 'Na':2.2, + 'Mg':3.3, + 'Al':9.0, + 'Si':26.2, + 'K': 1.2, + 'Ca':5.2, + 'Mn':0.1, + 'Fe':6.2, + 'Ti':0.5}, + 'density':3e3} + +# Fiducial volume (m) +FIDUCIAL_RADIUS = 5 + +# Fiducial volume (m^3) +FIDUCIAL_VOLUME = 4*pi*FIDUCIAL_RADIUS**3/3 + +# proton mass from the PDG (2018) +PROTON_MASS = 0.938 # GeV + +# proton mass from the PDG (2018) +ATOMIC_MASS_UNIT = 0.931 # GeV + +# mass of Xenon in atomic units +XENON_MASS = 131.293 + +# mass of Neon in atomic units +NEON_MASS = 20.18 + +# mass of argon in atomic units +ARGON_MASS = 39.948 + +# mass of germanium in atomic units +GERMANIUM_MASS = 72.64 + +# mass of tungsten in atomic units +TUNGSTEN_MASS = 183.84 + +# mass of oxygen in atomic units +OXYGEN_MASS = 15.999 + +# mass of silicon in atomic units +SILICON_MASS = 28.0855 + +# mass of iron in atomic units +IRON_MASS = 55.845 + +# mass of magnesium in atomic units +MAGNESIUM_MASS = 24.305 + +# galactic escape velocity (m/s) +# from Tom Caldwell's thesis page 25 +ESCAPE_VELOCITY = 244e3 + +# conversion constant from PDG +HBARC = 197.326978812e-3 # GeV fm + +# Avogadros number (kg^-1) +N0 = 6.02214085774e26 + +def get_probability(r, l): + """ + Returns the probability of a dark photon decaying in the SNO detector from + a dark matter decay distributed uniformly in the Earth. Assumes that the + depth of SNO is much larger than the dimensions of the SNO detector. + """ + if r <= h: + theta_min = 0 + elif r <= 2*R - h: + theta_min = pi - np.arccos((h**2 + r**2 - 2*R*h)/(2*r*(R-h))) + else: + return 0 + return (1 + np.cos(theta_min))*np.exp(-r/l)/(2*l) + +@lru_cache() +def get_probability2(l): + p, err = quad(get_probability, 0, min(2*R-h,10*l), args=(l,), epsabs=0, epsrel=1e-7, limit=1000) + + return p + +def get_event_rate(m, cs0, l, A): + """ + Returns the event rate of leptons produced from self-destructing dark + matter in the SNO detector for a given dark matter mass m, a cross section + cs, and a mediator decay length l. + """ + # For now we assume the event rate is constant throughout the earth, so we + # are implicitly assuming that the cross section is pretty small. + flux = DM_VELOCITY*DM_DENSITY/m + + #p, err = quad(get_probability, 0, min(2*R-h,10*l), args=(l,), epsabs=0, epsrel=1e-7, limit=1000) + p = get_probability2(l) + + cs = get_cross_section(cs0, m, A) + + # FIXME: factor of 2 because the DM particle decays into two mediators? + return p*cs*(N0/A)*flux*FIDUCIAL_VOLUME + +def get_event_rate_sno(m, cs0, l, composition): + """ + Returns the event rate of leptons produced from self-destructing dark + matter in the SNO detector for a given dark matter mass m, a cross section + cs, and a mediator decay length l. + """ + rate = 0.0 + + for element, mass_fraction in composition['composition'].items(): + rate += mass_fraction/100*get_event_rate(m,cs0,l,element_mass[element]) + + return rate*composition['density'] + +def get_nuclear_form_factor(A, e): + """ + Returns the nuclear form factor for a WIMP-nucleus interaction. + + From Tom Caldwell's thesis page 24. + + Also used Mark Pepin's thesis page 50 + """ + # mass of xenon nucleus + mn = A*ATOMIC_MASS_UNIT + + # calculate approximate size of radius + s = 0.9 # fm + a = 0.52 # fm + c = 1.23*A**(1/3) - 0.60 # fm + r1 = np.sqrt(c**2 + (7/3)*pi**2*a**2 - 5*s**2) + + q = np.sqrt(2*mn*e) + + if q*r1/HBARC < 1e-10: + return 1.0 + + # Helm form factor + # from Mark Pepin's thesis page 50 + f = 3*spherical_jn(1,q*r1/HBARC)*np.exp(-(q*s/HBARC)**2/2)/(q*r1/HBARC) + + return f + +def get_cross_section(cs0, m, A): + """ + Returns the WIMP cross section from the target-independent WIMP-nucleon + cross section cs0 at zero momentum transfer. + + From Tom Caldwell's thesis page 21. + """ + # mass of xenon nucleus + mn = A*ATOMIC_MASS_UNIT + + # reduced mass of the nucleus and the WIMP + mr = (m*mn)/(m + mn) + + # reduced mass of the proton and the WIMP + mp = (m*PROTON_MASS)/(m + PROTON_MASS) + + return cs0*A**2*mr**2/mp**2 + +def get_differential_event_rate_xenon(e, m, cs0, A): + """ + Returns the event rate of WIMP scattering in the Xenon 100T detector. + """ + # mass of nucleus + mn = A*ATOMIC_MASS_UNIT + + # reduced mass of the nucleus and the WIMP + mr = (m*mn)/(m + mn) + + cs = get_cross_section(cs0, m, A) + + v0 = DM_VELOCITY + + vesc = ESCAPE_VELOCITY + + # earth's velocity through the galaxy + ve = EARTH_VELOCITY + + # minimum wimp velocity needed to produce a recoil of energy e + vmin = np.sqrt(mn*e/2)*(mn+m)*SPEED_OF_LIGHT/(mn*m) + + f = get_nuclear_form_factor(A, e) + + x = vmin/v0 + y = ve/v0 + z = vesc/v0 + + # Equation 3.49 in Mark Pepin's thesis + k0 = (pi*v0**2)**(3/2) + + # Equation 3.49 in Mark Pepin's thesis + k1 = k0*(erf(z)-(2/np.sqrt(pi))*z*np.exp(-z**2)) + + # From Mark Pepin's CDMS thesis page 59 + if x <= z - y: + I = (k0/k1)*(DM_DENSITY/(2*ve))*(erf(x+y) - erf(x-y) - (4/np.sqrt(pi))*y*np.exp(-z**2)) + elif x <= y + z: + I = (k0/k1)*(DM_DENSITY/(2*ve))*(erf(z) - erf(x-y) - (2/np.sqrt(pi))*(y+z-x)*np.exp(-z**2)) + else: + return 0 + + return (N0/A)*(mn/(2*m*mr**2))*cs*f**2*I*SPEED_OF_LIGHT**2 + +def get_event_rate_xenon(m, cs, A, threshold): + """ + Returns the event rate of WIMP scattering in a dark matter detector using + an element with atomic mass A. Rate is in Hz kg^-1. + """ + # mass of nucleus + mn = A*ATOMIC_MASS_UNIT + + # reduced mass of the nucleus and the WIMP + mr = (m*mn)/(m + mn) + + vesc = ESCAPE_VELOCITY + + # earth's velocity through the galaxy + ve = EARTH_VELOCITY + + # max recoil they looked for was 40 keV + emax = 2*mr**2*((ve+vesc)/SPEED_OF_LIGHT)**2/mn + + rate, err = quad(get_differential_event_rate_xenon, threshold, emax, args=(m,cs,A), epsabs=0, epsrel=1e-2, limit=1000) + + return rate + +if __name__ == '__main__': + import matplotlib.pyplot as plt + + ls = np.logspace(-1,8,1000) + + cs = 1e-50 # cm^2 + + # FIXME: should use water density for L < 1 m, silicon density for L ~ 1 + # km, and iron for L >> 1 km + rate = np.array([get_event_rate_sno(1, cs*1e-4, l, water) for l in ls]) + + colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] + + plt.figure() + plt.subplot(111) + plt.plot(ls, rate/np.max(rate),color=colors[0]) + plt.xlabel("Mediator Decay Length (m)") + plt.ylabel("Event Rate (arbitrary units)") + plt.axvline(x=FIDUCIAL_RADIUS, color=colors[1], ls='--', label="SNO radius") + plt.axvline(x=h, ls='--', color=colors[2], label="Depth of SNO") + plt.axvline(x=R, ls='--', color=colors[3], label="Earth radius") + plt.gca().set_xscale('log') + plt.gca().set_xlim((ls[0], ls[-1])) + plt.title("Event Rate of Self-Destructing Dark Matter in SNO") + plt.legend() + + # threshold is ~5 keVr + xenon_100t_threshold = 5e-6 + + # threshold is ~1 keVr + cdms_threshold = 1e-6 + + # threshold is ~100 eV + # FIXME: is this correct? + cresst_threshold = 1e-7 + + ms = np.logspace(-2,3,200) + cs0s = np.logspace(-50,-40,200) + + mm, cs0cs0 = np.meshgrid(ms, cs0s) + + rate1 = np.empty(mm.shape) + rate2 = np.empty(mm.shape) + rate3 = np.empty(mm.shape) + rate4 = np.empty(mm.shape) + rate5 = np.empty(mm.shape) + rate6 = np.empty(mm.shape) + for i in range(mm.shape[0]): + print("\r%i/%i" % (i+1,mm.shape[0]),end='') + for j in range(mm.shape[1]): + rate1[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, XENON_MASS, xenon_100t_threshold) + rate2[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1.0, water) + rate3[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, GERMANIUM_MASS, cdms_threshold) + rate4[i,j] = get_event_rate_xenon(mm[i,j], cs0cs0[i,j]*1e-4, TUNGSTEN_MASS, cresst_threshold) + rate5[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1e3, norite) + rate6[i,j] = get_event_rate_sno(mm[i,j], cs0cs0[i,j]*1e-4, 1e6, mantle) + print() + + # Fiducial volume of the Xenon1T detector is 1042 +/- 12 kg + # from arxiv:1705.06655 + xenon_100t_fiducial_volume = 1042 # kg + + # Livetime of the Xenon1T results + # from arxiv:1705.06655 + xenon_100t_livetime = 34.2 # days + + plt.figure() + plt.subplot(111) + plt.gca().set_xscale('log') + plt.gca().set_yscale('log') + CS1 = plt.contour(mm,cs0cs0,rate1*3600*24*xenon_100t_fiducial_volume*xenon_100t_livetime,[10.0], colors=[colors[0]]) + CS2 = plt.contour(mm,cs0cs0,rate2*3600*24*668.8,[10.0],colors=[colors[1]]) + CS3 = plt.contour(mm,cs0cs0,rate3*3600*24*70.10,[10.0],colors=[colors[2]]) + # FIXME: I used 2.39 kg day because that's what CRESST-3 reports in their paper + # but! I only use Tungsten here so do I need to multiply by the mass fraction of tungsten? + CS4 = plt.contour(mm,cs0cs0,rate4*3600*24*2.39,[10.0],colors=[colors[3]]) + CS5 = plt.contour(mm,cs0cs0,rate5*3600*24*668.8,[10.0],colors=[colors[1]], linestyles=['dashed']) + CS6 = plt.contour(mm,cs0cs0,rate6*3600*24*668.8,[10.0],colors=[colors[1]], linestyles=['dotted']) + plt.clabel(CS1, inline=1, fmt="XENON1T", fontsize=10, use_clabeltext=True) + plt.clabel(CS2, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1 m)", fontsize=10) + plt.clabel(CS3, inline=1, fmt="CDMSLite", fontsize=10) + plt.clabel(CS4, inline=1, fmt="CRESST-3", fontsize=10) + plt.clabel(CS5, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1 km)", fontsize=10) + plt.clabel(CS6, inline=1, fmt=r"SNO ($\mathrm{L}_V$ = 1000 km)", fontsize=10) + plt.xlabel(r"$m_\chi$ (GeV)") + plt.ylabel(r"WIMP-nucleon scattering cross section ($\mathrm{cm}^2$)") + plt.title("Self-Destructing Dark Matter Limits") + + x = np.linspace(0,300e-6,1000) + + # reproducing Figure 2.1 in Tom Caldwell's thesis + plt.figure() + plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(XENON_MASS, x)**2,x)), label="Xe") + plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(NEON_MASS, x)**2,x)), label="Ne") + plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(ARGON_MASS, x)**2,x)), label="Ar") + plt.plot(x*1e6, list(map(lambda x: get_nuclear_form_factor(GERMANIUM_MASS, x)**2,x)), label="Ge") + plt.xlabel("Recoil Energy (keV)") + plt.ylabel(r"$F^2(E)$") + plt.gca().set_yscale('log') + plt.legend() + plt.gca().set_ylim((1e-5,1)) + plt.gca().set_xlim((0,300)) + plt.title("Nuclear Form Factors for various elements") + plt.show() diff --git a/utils/convert-genie-to-gst b/utils/convert-genie-to-gst new file mode 100755 index 0000000..b2861e5 --- /dev/null +++ b/utils/convert-genie-to-gst @@ -0,0 +1,21 @@ +#!/usr/bin/env python +from subprocess import check_call +from os.path import split, splitext, join +import os + +if __name__ == '__main__': + import argparse + + parser = argparse.ArgumentParser("script to convert full GENIE root files to the reduced gst ROOT format") + parser.add_argument("filenames", nargs="+", help="GENIE root files") + parser.add_argument("--dest", required=True, help="destination directory") + args = parser.parse_args() + + for filename in args.filenames: + head, tail = split(filename) + root, ext = splitext(tail) + output = join(args.dest, root) + ".ntuple.root" + cmd = ["gntpc","-f","gst","-i",filename,"-o",output] + print(" ".join(cmd)) + with open(os.devnull,"w") as devnull: + check_call(cmd, stdout=devnull, stderr=devnull) diff --git a/utils/convert-ratdb b/utils/convert-ratdb new file mode 100755 index 0000000..5dd6ee4 --- /dev/null +++ b/utils/convert-ratdb @@ -0,0 +1,43 @@ +#!/usr/bin/env python +""" +Script to convert the PMT database in RAT to a CSV file (with spaces instead of commas). +""" +from __future__ import print_function, division +import yaml + +if __name__ == '__main__': + import argparse + import numpy as np + + parser = argparse.ArgumentParser("convert a RATDB file -> CSV file") + parser.add_argument("filename", help="RATDB filename") + parser.add_argument("-o", "--output", help="output filename", required=True) + args = parser.parse_args() + + s = '' + with open(args.filename) as f: + for line in f: + if line.startswith('//'): + continue + s += line + + data = yaml.load(s) + + csv = np.zeros((len(data['x']),7),dtype=np.float64) + csv[:,0] = data['x'] + csv[:,1] = data['y'] + csv[:,2] = data['z'] + csv[:,3] = data['u'] + csv[:,4] = data['v'] + csv[:,5] = data['w'] + csv[:,6] = data['pmt_type'] + + # reverse PMT normals + csv[:,3] = -csv[:,3] + csv[:,4] = -csv[:,4] + csv[:,5] = -csv[:,5] + + header = "Format: x, y, z, u, v, w, pmt_type" + + np.savetxt(args.output, csv, fmt=['%10.2f']*6 + ['%i'], + header=header) diff --git a/utils/index.py b/utils/index.py new file mode 100644 index 0000000..4edcd59 --- /dev/null +++ b/utils/index.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python +""" +This is a script to calculate the index of refraction of water as a function of +wavelength, temperature, and pressure and fit it to a linear approximation in a +wavelength region where the PMTs are sensitive to. The data comes from [1]. + +[1] "Refractive index of Water and Steam as Function of Wavelength, Temperature, and Density". Schiebener et al. 1989. +""" +from __future__ import print_function, division +import numpy as np + +A0 = 0.243905091 +A1 = 9.53518094e-3 +A2 = -3.64358110e-3 +A3 = 2.65666426e-4 +A4 = 1.59189325e-3 +A5 = 2.45733798e-3 +A6 = 0.897478251 +A7 = -1.63066183e-2 +UV = 0.2292020 +IR = 5.432937 + +def get_index(p, wavelength, T): + """ + Returns the index of pure water for a given density, wavelength, and + temperature. The density should be in units of kg/m^3, the wavelength in nm, + and the temperature in Celsius. + """ + # normalize the density, temperature, and pressure + p = p/1000.0 + wavelength = wavelength/589.0 + T = (T+273.15)/273.15 + + # first we compute the right hand side of Equation 7 + c = A0 + A1*p + A2*T + A3*wavelength**2*T + A4/wavelength**2 + A5/(wavelength**2-UV**2) + A6/(wavelength**2-IR**2) + A7*p**2 + c *= p + + return np.sqrt((2*c+1)/(1-c)) + +def approximate_index(x, a, b): + """ + Returns the index of refraction using the approximation: + + n = a + b/x**2 + """ + return a + b/x**2 + +def chi2(pars, x, data): + """ + Returns the squares of the difference between the data and the model for + the inverse of the index of refraction, i.e. + + 1/n = lambda**2/(a*lambda**2+b) + """ + a, b = pars + return np.sum((approximate_index(x,a,b)-data)**2) + +if __name__ == '__main__': + import matplotlib.pyplot as plt + import argparse + from scipy.optimize import fmin + + parser = argparse.ArgumentParser("fit for an approximation to the inverse of the index of refraction of water") + parser.add_argument("--temperature", type=float, default=10.0) + parser.add_argument("--density", type=float, default=1000.0) + parser.add_argument("--low", type=float, default=300.0) + parser.add_argument("--high", type=float, default=600.0) + args = parser.parse_args() + + x = np.linspace(args.low,args.high,1000) + n = get_index(args.density,x,args.temperature) + xmin = fmin(chi2,[1,1],args=(x,1/n)) + + print(xmin) + + a, b = xmin + y = approximate_index(x,a,b) + + plt.plot(x,1/n,label="data") + plt.plot(x,y,label="fit") + plt.xlabel("Wavelength (nm)") + plt.ylabel("1/n") + plt.title("Index of Refraction of Water at 273.15 K") + plt.legend() + plt.show() diff --git a/utils/plot.py b/utils/plot.py new file mode 100755 index 0000000..6a2fe73 --- /dev/null +++ b/utils/plot.py @@ -0,0 +1,78 @@ +#!/usr/bin/env python +from __future__ import print_function, division +import yaml + +SNOMAN_MASS = { + 20: 0.511, + 21: 0.511, + 22: 105.658, + 23: 105.658 +} + +if __name__ == '__main__': + import argparse + import matplotlib.pyplot as plt + import numpy as np + + parser = argparse.ArgumentParser("plot fit results") + parser.add_argument("filenames", nargs='+', help="input files") + args = parser.parse_args() + + for filename in args.filenames: + print(filename) + with open(filename) as f: + data = yaml.load(f.read()) + + dx = [] + dy = [] + dz = [] + e = [] + for event in data['data']: + # get the particle ID + id = event['mctk'][-1]['id'] + mass = SNOMAN_MASS[id] + # for some reason it's the *second* track which seems to contain the + # initial energy + true_energy = event['mctk'][-1]['energy'] + # The MCTK bank has the particle's total energy (except for neutrons) + # so we need to convert it into kinetic energy + ke = true_energy - mass + energy = event['ev'][0]['fit'][0]['energy'] + e.append(energy-ke) + true_posx = event['mcvx'][0]['posx'] + posx = event['ev'][0]['fit'][0]['posx'] + dx.append(posx-true_posx) + true_posy = event['mcvx'][0]['posy'] + posy = event['ev'][0]['fit'][0]['posy'] + dy.append(posy-true_posy) + true_posz = event['mcvx'][0]['posz'] + posz = event['ev'][0]['fit'][0]['posz'] + dz.append(posz-true_posz) + + print("dT = %.2g +/- %.2g" % (np.mean(e), np.std(e))) + print("dx = %4.2g +/- %.2g" % (np.mean(dx), np.std(dx))) + print("dy = %4.2g +/- %.2g" % (np.mean(dy), np.std(dy))) + print("dz = %4.2g +/- %.2g" % (np.mean(dz), np.std(dz))) + + plt.figure(1) + plt.hist(e, histtype='step', label=filename) + plt.xlabel("Kinetic Energy difference (MeV)") + plt.figure(2) + plt.hist(dx, histtype='step', label=filename) + plt.xlabel("X Position difference (cm)") + plt.figure(3) + plt.hist(dy, histtype='step', label=filename) + plt.xlabel("Y Position difference (cm)") + plt.figure(4) + plt.hist(dz, histtype='step', label=filename) + plt.xlabel("Z Position difference (cm)") + + plt.figure(1) + plt.legend() + plt.figure(2) + plt.legend() + plt.figure(3) + plt.legend() + plt.figure(4) + plt.legend() + plt.show() -- cgit