#!/usr/bin/env python
"""
Script for plotting solar fluxes from
http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0403i/index.html.
"""
from __future__ import print_function, division
import numpy as np
import os
from sddm.plot import despine
from sddm import splitext
# Data is in the form: Fluxes in bins of neutrino energy (equally spaced bins
# in logE with 10 bins per decade with the low edge of the first bin at 100
# MeV) and zenith angle (20 bins equally spaced in cos(zenith) with bin width
# 0.1), integrated over azimuth. Note logE means log to base e. Fluxes below 10
# GeV are from the 3D calculation.
# The files all have the same format. After the initial comment lines (starting
# with a # character), the files contain one line per bin. No smoothoing
# between bins has been done. The columns are:
#
# 1. Energy at centre of bin in GeV
# 2. Zenith files: Cos(zenith) at centre of bin
# Azimuth files: Azimuth at centre of bin (degrees)
# 3. Flux in dN/dlogE in /m**2/steradian/sec
# 4. MC Statistical error on the flux
# 5. Number of unweighted events entering the bin (not too useful)
if __name__ == '__main__':
import argparse
import matplotlib
import glob
from sddm import setup_matplotlib
parser = argparse.ArgumentParser("plot solar fluxes")
parser.add_argument("filenames", nargs='+', help="filenames of flux files")
parser.add_argument("--save", action="store_true", default=False, help="save plots")
args = parser.parse_args()
setup_matplotlib(args.save)
pre { line-height: 125%; }
td.linenos .normal { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; }
span.linenos { color: inherit; background-color: transparent; padding-left: 5px; padding-right: 5px; }
td.linenos .special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }
span.linenos.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }
.highlight .hll { background-color: #ffffcc }
.highlight .c { color: #888888 } /* Comment */
.highlight .err { color: #a61717; background-color: #e3d2d2 } /* Error */
.highlight .k { color: #008800; font-weight: bold } /* Keyword */
.highlight .ch { color: #888888 } /* Comment.Hashbang */
.highlight .cm { color: #888888 } /* Comment.Multiline */
.highlight .cp { color: #cc0000; font-weight: bold } /* Comment.Preproc */
.highlight .cpf { color: #888888 } /* Comment.PreprocFile */
.highlight .c1 { color: #888888 } /* Comment.Single */
.highlight .cs { color: #cc0000; font-weight: bold; background-color: #fff0f0 } /* Comment.Special */
.highlight .gd { color: #000000; background-color: #ffdddd } /* Generic.Deleted */
.highlight .ge { font-style: italic } /* Generic.Emph */
.highlight .ges { font-weight: bold; font-style: italic } /* Generic.EmphStrong */
.highlight .gr { color: #aa0000 } /* Generic.Error */
.highlight .gh { color: #333333 } /* Generic.Heading */
.highlight .gi { color: #000000; background-color: #ddffdd } /* Generic.Inserted */
.highlight .go { color: #888888 } /* Generic.Output */
.highlight .gp { color: #555555 } /* Generic.Prompt */
.highlight .gs { font-weight: bold } /* Generic.Strong */
.highlight .gu { color: #666666 } /* Generic.Subheading */
.highlight .gt { color: #aa0000 } /* Generic.Traceback */
.highlight .kc { color: #008800; font-weight: bold } /* Keyword.Constant */
.highlight .kd { color: #008800; font-weight: bold } /* Keyword.Declaration */
.highlight .kn { color: #008800; font-weight: bold } /* Keyword.Namespace */
.highlight .kp { color: #008800 } /* Keyword.Pseudo */
.highlight .kr { color: #008800; font-weight: bold } /* Keyword.Reserved */
.highlight .kt { color: #888888; font-weight: bold } /* Keyword.Type */
.highlight .m { color: #0000DD; font-weight: bold } /* Literal.Number */
.highlight .s { color: #dd2200; background-color: #fff0f0 } /* Literal.String */
.highlight .na { color: #336699 } /* Name.Attribute */
.highlight .nb { color: #003388 } /* Name.Builtin */
.highlight .nc { color: #bb0066; font-weight: bold } /* Name.Class */
.highlight .no { color: #003366; font-weight: bold } /* Name.Constant */
.highlight .nd { color: #555555 } /* Name.Decorator */
.highlight .ne { color: #bb0066; font-weight: bold } /* Name.Exception */
.highlight .nf { color: #0066bb; font-weight: bold } /* Name.Function */
.highlight .nl { color: #336699; font-style: italic } /* Name.Label */
.highlight .nn { color: #bb0066; font-weight: bold } /* Name.Namespace */
.highlight .py { color: #336699; font-weight: bold } /* Name.Property */
.highlight .nt { color: #bb0066; font-weight: bold } /* Name.Tag */
.highlight .nv { color: #336699 } /* Name.Variable */
.highlight .ow { color: #008800 } /* Operator.Word */
.highlight .w { color: #bbbbbb } /* Text.Whitespace */
.highlight .mb { color: #0000DD; font-weight: bold } /* Literal.Number.Bin */
.highlight .mf { color: #0000DD; font-weight: bold } /* Literal.Number.Float */
.highlight .mh { color: #0000DD; font-weight: bold } /* Literal.Number.Hex */
.highlight .mi { color: #0000DD; font-weight: bold } /* Literal.Number.Integer */
.highlight .mo { color: #0000DD; font-weight: bold } /* Literal.Number.Oct */
.highlight .sa { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Affix */
.highlight .sb { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Backtick */
.highlight .sc { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Char */
.highlight .dl { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Delimiter */
.highlight .sd { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Doc */
.highlight .s2 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Double */
.highlight .se { color: #0044dd; background-color: #fff0f0 } /* Literal.String.Escape */
.highlight .sh { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Heredoc */
.highlight .si { color: #3333bb; background-color: #fff0f0 } /* Literal.String.Interpol */
.highlight .sx { color: #22bb22; background-color: #f0fff0 } /* Literal.String.Other */
.highlight .sr { color: #008800; background-color: #fff0ff } /* Literal.String.Regex */
.highlight .s1 { color: #dd2200; background-color: #fff0f0 } /* Literal.String.Single */
.highlight .ss { color: #aa6600; background-color: #fff0f0 } /* Literal.String.Symbol */
.highlight .bp { color: #003388 } /* Name.Builtin.Pseudo */
.highlight .fm { color: #0066bb; font-weight: bold } /* Name.Function.Magic */
.highlight .vc { color: #336699 } /* Name.Variable.Class */
.highlight .vg { color: #dd7700 } /* Name.Variable.Global */
.highlight .vi { color: #3333bb } /* Name.Variable.Instance */
.highlight .vm { color: #336699 } /* Name.Variable.Magic */
.highlight .il { color: #0000DD; font-weight: bold } /* Literal.Number.Integer.Long */#!/usr/bin/env python
# Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option)
# any later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
# more details.
#
# You should have received a copy of the GNU General Public License along with
# this program. If not, see <https://www.gnu.org/licenses/>.
"""
Very simple script to analyze GENIE MC to estimate the expected event rate from
atmospheric neutrino events before and after different cuts. To run it simply
pass the filename of a GENIE ntuple file:
$ ./analyze-genie-mc mc_atm_nu_genie_010000_0.gst.root
Note: The input files should be in the GENIE "ntuple" format which you can
create by using the gntpc utility. For example:
$ gntpc -i /sno/output/d2o-leta/mc_atmospherics/atmospheric_neutrinos_genie/root/mc_atm_nu_genie_010000_0.root -f gst
will create a new file mc_atm_nu_genie_010000_0.gst.root in the current working
directory in the ntuple file format.
"""
from __future__ import print_function, division
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()
