#!/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()