#!/usr/bin/env python
from __future__ import print_function, division
import yaml
import numpy as np
from scipy.stats import iqr
from matplotlib.lines import Line2D
# 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")
IDP_E_MINUS = 20
IDP_MU_MINUS = 22
SNOMAN_MASS = {
20: 0.511,
21: 0.511,
22: 105.658,
23: 105.658
}
def plot_hist(x, label=None):
# determine the bin width using the Freedman Diaconis rule
# see https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule
h = 2*iqr(x)/len(x)**(1/3)
n = max(int((np.max(x)-np.min(x))/h),10)
bins = np.linspace(np.min(x),np.max(x),n)
plt.hist(x, bins=bins, histtype='step', label=label)
def plot_legend(n):
plt.figure(n)
ax = plt.gca()
handles, labels = ax.get_legend_handles_labels()
new_handles = [Line2D([],[],c=h.get_edgecolor()) for h in handles]
plt.legend(handles=new_handles,labels=labels)
def get_stats(x):
"""
Returns a tuple (mean, error mean, std, error std) for the values in x.
The formula for the standard error on the standard deviation comes from
https://stats.stackexchange.com/questions/156518.
"""
mean = np.mean(x)
std = np.std(x)
n = len(x)
u4 = np.mean((x-mean)**4)
error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std)
return mean, std/np.sqrt(n), std, error
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 = []
dT = []
thetas = []
likelihood_ratio = []
t_electron = []
t_muon = []
for event in data['data']:
# get the particle ID
id = event['mctk'][-1]['id']
if id not in event['ev'][0]['fit']:
continue
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'][id]['energy']
dT.append(energy-ke)
true_posx = event['mcvx'][0]['posx']
posx = event['ev'][0]['fit'][id]['posx']
dx.append(posx-true_posx)
true_posy = event['mcvx'][0]['posy']
posy = event['ev'][0]['fit'][id]['posy']
dy.append(posy-true_posy)
true_posz = event['mcvx'][0]['posz']
posz = event['ev'][0]['fit'][id]['posz']
dz.append(posz-true_posz)
dirx = event['mctk'][-1]['dirx']
diry = event['mctk'][-1]['diry']
dirz = event['mctk'][-1]['dirz']
true_dir = [dirx,diry,dirz]
true_dir = np.array(true_dir)/np.linalg.norm(true_dir)
theta = event['ev'][0]['fit'][id]['theta']
phi = event['ev'][0]['fit'][id]['phi']
dir = [np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)]
dir = np.array(dir)/np.linalg.norm(dir)
thetas.append(np.degrees(np.arccos(np.dot(true_dir,dir))))
if IDP_E_MINUS in event['ev'][0]['fit'] and IDP_MU_MINUS in event['ev'][0]['fit']:
fmin_electron = event['ev'][0]['fit'][IDP_E_MINUS]['fmin']
fmin_muon = event['ev'][0]['fit'][IDP_MU_MINUS]['fmin']
likelihood_ratio.append(fmin_muon-fmin_electron)
if IDP_E_MINUS in event['ev'][0]['fit']:
t_electron.append(event['ev'][0]['fit'][IDP_E_MINUS]['time'])
if IDP_MU_MINUS in event['ev'][0]['fit']:
t_muon.append(event['ev'][0]['fit'][IDP_MU_MINUS]['time'])
mean, mean_error, std, std_error = get_stats(dT)
print("dT = %.2g +/- %.2g" % (mean, mean_error))
print("std(dT) = %.2g +/- %.2g" % (std, std_error))
mean, mean_error, std, std_error = get_stats(dx)
print("dx = %4.2g +/- %.2g" % (mean, mean_error))
print("std(dx) = %4.2g +/- %.2g" % (std, std_error))
mean, mean_error, std, std_error = get_stats(dy)
print("dy = %4.2g +/- %.2g" % (mean, mean_error))
print("std(dy) = %4.2g +/- %.2g" % (std, std_error))
mean, mean_error, std, std_error = get_stats(dz)
print("dz = %4.2g +/- %.2g" % (mean, mean_error))
print("std(dz) = %4.2g +/- %.2g" % (std, std_error))
mean, mean_error, std, std_error = get_stats(thetas)
print("std(theta) = %4.2g +/- %.2g" % (std, std_error))
plt.figure(1)
plot_hist(dT, label=filename)
plt.xlabel("Kinetic Energy difference (MeV)")
plt.figure(2)
plot_hist(dx, label=filename)
plt.xlabel("X Position difference (cm)")
plt.figure(3)
plot_hist(dy, label=filename)
plt.xlabel("Y Position difference (cm)")
plt.figure(4)
plot_hist(dz, label=filename)
plt.xlabel("Z Position difference (cm)")
plt.figure(5)
plot_hist(thetas, label=filename)
plt.xlabel(r"$\theta$ (deg)")
plt.figure(6)
plot_hist(likelihood_ratio, label=filename)
plt.xlabel(r"Log Likelihood Ratio ($e/\mu$)")
plt.figure(7)
plot_hist(np.array(t_electron)/1e3/60.0, label=filename)
plt.xlabel(r"Electron Fit time (minutes)")
plt.figure(8)
plot_hist(np.array(t_muon)/1e3/60.0, label=filename)
plt.xlabel(r"Muon Fit time (minutes)")
plot_legend(1)
plot_legend(2)
plot_legend(3)
plot_legend(4)
plot_legend(5)
plot_legend(6)
plot_legend(7)
plot_legend(8)
plt.show()
' href='#n114'>114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#!/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/>.
"""
Script to plot the energy and time difference distribution for neutrons. To run
it just run:
$ ./plot-neutrons [list of fit results]
"""
from __future__ import print_function, division
import numpy as np
from scipy.stats import iqr, poisson
from scipy.stats import iqr, norm, beta
from sddm.stats import *
import emcee
from sddm.dc import estimate_errors, EPSILON
import nlopt
from sddm import printoptions
particle_id = {20: 'e', 22: 'u'}
if __name__ == '__main__':
import argparse
import numpy as np
import pandas as pd
import sys
import h5py
from sddm.plot_energy import *
from sddm.plot import *
from sddm import setup_matplotlib
from sddm.utils import correct_energy_bias
parser = argparse.ArgumentParser("plot fit results")
parser.add_argument("filenames", nargs='+', help="input files")
parser.add_argument("--save", action='store_true', default=False, help="save corner plots for backgrounds")
parser.add_argument("--mc", nargs='+', required=True, help="atmospheric MC files")
args = parser.parse_args()
setup_matplotlib(args.save)
import matplotlib.pyplot as plt
# Loop over runs to prevent using too much memory
evs = []
rhdr = pd.concat([read_hdf(filename, "rhdr").assign(filename=filename) for filename in args.filenames],ignore_index=True)
for run, df in rhdr.groupby('run'):
evs.append(get_events(df.filename.values, merge_fits=True))
ev = pd.concat(evs)
ev = correct_energy_bias(ev)
# Note: We loop over the MC filenames here instead of just passing the
# whole list to get_events() because I had to rerun some of the MC events
# using SNOMAN and so most of the runs actually have two different files
# and otherwise the GTIDs will clash
ev_mcs = []
for filename in args.mc:
ev_mcs.append(get_events([filename], merge_fits=True, mc=True))
ev_mc = pd.concat(ev_mcs)
ev_mc = correct_energy_bias(ev_mc)
ev = ev.reset_index()
ev_mc = ev_mc.reset_index()
# remove events 200 microseconds after a muon
ev = ev.groupby('run',group_keys=False).apply(muon_follower_cut)
# 00-orphan cut
ev = ev[(ev.gtid & 0xff) != 0]
ev_mc = ev_mc[(ev_mc.gtid & 0xff) != 0]
neutrons = ev[ev.neutron]
neutrons_mc = ev_mc[ev_mc.neutron]
atm = ev[ev.signal & ev.prompt & ev.atm]
atm_mc = ev_mc[ev_mc.signal & ev_mc.prompt & ev_mc.atm]
# Drop events without fits
atm = atm[~np.isnan(atm.fmin)]
atm_mc = atm_mc[~np.isnan(atm_mc.fmin)]
atm = atm[atm.psi < 6]
atm_mc = atm_mc[atm_mc.psi < 6]
atm = pd.merge(atm,neutrons,left_on=['run','gtid'],right_on=['run','atm_gtid'],suffixes=('','_neutron'))
atm_mc = pd.merge(atm_mc,neutrons_mc,left_on=['run','gtid'],right_on=['run','atm_gtid'],suffixes=('','_neutron'))
print("neutrons with nhit > 100")
print(atm[atm.nhit_neutron >= 100][['run','gtid_neutron','nhit_neutron','nhit_cal_neutron']])
fig = plt.figure(1)
plt.hist(atm.nhit_cal_neutron.values,bins=np.linspace(0,100,101),histtype='step',color='C0',label="Data")
weights = np.tile(len(atm)/len(atm_mc),len(atm_mc))
plt.hist(atm_mc.nhit_cal_neutron.values,bins=np.linspace(0,100,101),weights=weights,histtype='step',color='C1',label="Monte Carlo")
plt.xlabel("Nhit")
despine(fig,trim=True)
plt.tight_layout()
plot_legend(1)
if args.save:
plt.savefig("neutron_nhit_cal.pdf")
plt.savefig("neutron_nhit_cal.eps")
else:
plt.title("Neutron Nhit Distribution")
fig = plt.figure(2)
dt = (atm.gtr_neutron - atm.gtr)/1e6;
bins = np.linspace(20e-3,250,101)
plt.hist(dt,bins=bins,histtype='step',color='C0',label="Data")
dt = (atm_mc.gtr_neutron - atm_mc.gtr)/1e6;
plt.hist(dt,bins=bins,weights=weights,histtype='step',color='C1',label="Monte Carlo")
plt.xlabel(r"$\Delta$ t (ms)")
despine(fig,trim=True)
plt.tight_layout()
plot_legend(1)
if args.save:
plt.savefig("neutron_delta_t.pdf")
plt.savefig("neutron_delta_t.eps")
else:
plt.title(r"Neutron $\Delta t$ Distribution")
plt.show()
|