diff options
Diffstat (limited to 'utils')
| -rwxr-xr-x | utils/analyze-genie-mc | 373 | ||||
| -rwxr-xr-x | utils/calculate_limits.py | 419 | ||||
| -rwxr-xr-x | utils/convert-genie-to-gst | 21 | ||||
| -rwxr-xr-x | utils/convert-ratdb | 43 | ||||
| -rw-r--r-- | utils/index.py | 85 | ||||
| -rwxr-xr-x | utils/plot.py | 78 |
6 files changed, 1019 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() 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() |