aboutsummaryrefslogtreecommitdiff
path: root/utils/plot-fit-results
diff options
context:
space:
mode:
Diffstat (limited to 'utils/plot-fit-results')
-rwxr-xr-xutils/plot-fit-results321
1 files changed, 268 insertions, 53 deletions
diff --git a/utils/plot-fit-results b/utils/plot-fit-results
index ab08586..450de72 100755
--- a/utils/plot-fit-results
+++ b/utils/plot-fit-results
@@ -16,8 +16,9 @@
from __future__ import print_function, division
import numpy as np
-from scipy.stats import iqr
+from scipy.stats import iqr, norm
from matplotlib.lines import Line2D
+import itertools
# on retina screens, the default plots are way too small
# by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1
@@ -27,6 +28,12 @@ matplotlib.use("Qt5Agg")
matplotlib.rc('font', size=22)
+font = {'family':'serif', 'serif': ['computer modern roman']}
+matplotlib.rc('font',**font)
+
+matplotlib.rc('text', usetex=True)
+
+
IDP_E_MINUS = 20
IDP_MU_MINUS = 22
@@ -66,6 +73,59 @@ def get_stats(x):
error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std)
return mean, std/np.sqrt(n), std, error
+def iqr_std_err(x):
+ """
+ Returns the approximate standard deviation assuming the central part of the
+ distribution is gaussian.
+ """
+ x = x.dropna()
+ n = len(x)
+ if n == 0:
+ return np.nan
+ # see https://stats.stackexchange.com/questions/110902/error-on-interquartile-range
+ std = iqr(x.values)/1.3489795
+ return 1.573*std/np.sqrt(n)
+
+def iqr_std(x):
+ """
+ Returns the approximate standard deviation assuming the central part of the
+ distribution is gaussian.
+ """
+ x = x.dropna()
+ n = len(x)
+ if n == 0:
+ return np.nan
+ return iqr(x.values)/1.3489795
+
+def median(x):
+ """
+ Returns the median for all the non NaN values in a Series `x`.
+ """
+ x = x.dropna()
+ n = len(x)
+ if n == 0:
+ return np.nan
+ return np.median(x.values)
+
+def median_err(x):
+ """
+ Returns the approximate error on the median for all the non NaN values in a
+ Series `x`. The error on the median is approximated assuming the central
+ part of the distribution is gaussian.
+ """
+ x = x.dropna()
+ n = len(x)
+ if n == 0:
+ return np.nan
+ # First we estimate the standard deviation using the interquartile range.
+ # Here we are essentially assuming the central part of the distribution is
+ # gaussian.
+ std = iqr(x.values)/1.3489795
+ median = np.median(x.values)
+ # Now we estimate the error on the median for a gaussian
+ # See https://stats.stackexchange.com/questions/45124/central-limit-theorem-for-sample-medians.
+ return 1/(2*np.sqrt(n)*norm.pdf(median,median,std))
+
def std_err(x):
x = x.dropna()
mean = np.mean(x)
@@ -77,6 +137,108 @@ def std_err(x):
error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std)
return error
+# Taken from https://raw.githubusercontent.com/mwaskom/seaborn/c73055b2a9d9830c6fbbace07127c370389d04dd/seaborn/utils.py
+def despine(fig=None, ax=None, top=True, right=True, left=False,
+ bottom=False, offset=None, trim=False):
+ """Remove the top and right spines from plot(s).
+
+ fig : matplotlib figure, optional
+ Figure to despine all axes of, default uses current figure.
+ ax : matplotlib axes, optional
+ Specific axes object to despine.
+ top, right, left, bottom : boolean, optional
+ If True, remove that spine.
+ offset : int or dict, optional
+ Absolute distance, in points, spines should be moved away
+ from the axes (negative values move spines inward). A single value
+ applies to all spines; a dict can be used to set offset values per
+ side.
+ trim : bool, optional
+ If True, limit spines to the smallest and largest major tick
+ on each non-despined axis.
+
+ Returns
+ -------
+ None
+
+ """
+ # Get references to the axes we want
+ if fig is None and ax is None:
+ axes = plt.gcf().axes
+ elif fig is not None:
+ axes = fig.axes
+ elif ax is not None:
+ axes = [ax]
+
+ for ax_i in axes:
+ for side in ["top", "right", "left", "bottom"]:
+ # Toggle the spine objects
+ is_visible = not locals()[side]
+ ax_i.spines[side].set_visible(is_visible)
+ if offset is not None and is_visible:
+ try:
+ val = offset.get(side, 0)
+ except AttributeError:
+ val = offset
+ _set_spine_position(ax_i.spines[side], ('outward', val))
+
+ # Potentially move the ticks
+ if left and not right:
+ maj_on = any(
+ t.tick1line.get_visible()
+ for t in ax_i.yaxis.majorTicks
+ )
+ min_on = any(
+ t.tick1line.get_visible()
+ for t in ax_i.yaxis.minorTicks
+ )
+ ax_i.yaxis.set_ticks_position("right")
+ for t in ax_i.yaxis.majorTicks:
+ t.tick2line.set_visible(maj_on)
+ for t in ax_i.yaxis.minorTicks:
+ t.tick2line.set_visible(min_on)
+
+ if bottom and not top:
+ maj_on = any(
+ t.tick1line.get_visible()
+ for t in ax_i.xaxis.majorTicks
+ )
+ min_on = any(
+ t.tick1line.get_visible()
+ for t in ax_i.xaxis.minorTicks
+ )
+ ax_i.xaxis.set_ticks_position("top")
+ for t in ax_i.xaxis.majorTicks:
+ t.tick2line.set_visible(maj_on)
+ for t in ax_i.xaxis.minorTicks:
+ t.tick2line.set_visible(min_on)
+
+ if trim:
+ # clip off the parts of the spines that extend past major ticks
+ xticks = ax_i.get_xticks()
+ if xticks.size:
+ firsttick = np.compress(xticks >= min(ax_i.get_xlim()),
+ xticks)[0]
+ lasttick = np.compress(xticks <= max(ax_i.get_xlim()),
+ xticks)[-1]
+ ax_i.spines['bottom'].set_bounds(firsttick, lasttick)
+ ax_i.spines['top'].set_bounds(firsttick, lasttick)
+ newticks = xticks.compress(xticks <= lasttick)
+ newticks = newticks.compress(newticks >= firsttick)
+ ax_i.set_xticks(newticks)
+
+ yticks = ax_i.get_yticks()
+ if yticks.size:
+ firsttick = np.compress(yticks >= min(ax_i.get_ylim()),
+ yticks)[0]
+ lasttick = np.compress(yticks <= max(ax_i.get_ylim()),
+ yticks)[-1]
+ ax_i.spines['left'].set_bounds(firsttick, lasttick)
+ ax_i.spines['right'].set_bounds(firsttick, lasttick)
+ newticks = yticks.compress(yticks <= lasttick)
+ newticks = newticks.compress(newticks >= firsttick)
+ ax_i.set_yticks(newticks)
+
if __name__ == '__main__':
import argparse
import matplotlib.pyplot as plt
@@ -137,6 +299,7 @@ if __name__ == '__main__':
data['theta'] = np.degrees(np.arccos((true_dir*dir).sum(axis=-1)))
+
# only select fits which have at least 2 fits
data = data.groupby(['filename','run','gtid']).filter(lambda x: len(x) > 1)
data_true = data[data['fit_id1'] == data['mcgn_id']]
@@ -155,67 +318,119 @@ if __name__ == '__main__':
# 100 bins between 50 MeV and 1 GeV
bins = np.arange(50,1000,100)
+ markers = itertools.cycle(('o', 'v'))
+
+ # Default figure size. Currently set to my monitor width and height so that
+ # things are properly formatted
+ FIGSIZE = (13.78,7.48)
+
+ fig3, ax3 = plt.subplots(3,1,figsize=FIGSIZE,num=3,sharex=True)
+ fig4, ax4 = plt.subplots(3,1,figsize=FIGSIZE,num=4,sharex=True)
+
for id in [IDP_E_MINUS, IDP_MU_MINUS]:
events = data_true[data_true['mcgn_id'] == id]
pd_bins = pd.cut(events['T'],bins)
- dT = events.groupby(pd_bins)['dT'].agg(['mean','sem','std',std_err])
- dx = events.groupby(pd_bins)['dx'].agg(['mean','sem','std',std_err])
- dy = events.groupby(pd_bins)['dy'].agg(['mean','sem','std',std_err])
- dz = events.groupby(pd_bins)['dz'].agg(['mean','sem','std',std_err])
- theta = events.groupby(pd_bins)['theta'].agg(['mean','sem','std',std_err])
+ dT = events.groupby(pd_bins)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
+ dx = events.groupby(pd_bins)['dx'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
+ dy = events.groupby(pd_bins)['dy'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
+ dz = events.groupby(pd_bins)['dz'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
+ theta = events.groupby(pd_bins)['theta'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
label = 'Muon' if id == IDP_MU_MINUS else 'Electron'
T = (bins[1:] + bins[:-1])/2
- plt.figure(1)
- plt.errorbar(T,dT['mean']*100/T,yerr=dT['sem']*100/T,fmt='o',label=label)
- plt.xlabel("Kinetic Energy (MeV)")
- plt.ylabel("Energy bias (%)")
- plt.title("Energy Bias")
- plt.legend()
-
- plt.figure(2)
- plt.errorbar(T,dT['std']*100/T,yerr=dT['std_err']*100/T,fmt='o',label=label)
- plt.xlabel("Kinetic Energy (MeV)")
- plt.ylabel("Energy resolution (%)")
- plt.title("Energy Resolution")
- plt.legend()
-
- plt.figure(3)
- plt.errorbar(T,dx['mean'],yerr=dx['sem'],fmt='o',label='%s (x)' % label)
- plt.errorbar(T,dy['mean'],yerr=dy['sem'],fmt='o',label='%s (y)' % label)
- plt.errorbar(T,dz['mean'],yerr=dz['sem'],fmt='o',label='%s (z)' % label)
- plt.xlabel("Kinetic Energy (MeV)")
- plt.ylabel("Position bias (cm)")
- plt.title("Position Bias")
- plt.legend()
-
- plt.figure(4)
- plt.errorbar(T,dx['std'],yerr=dx['std_err'],fmt='o',label='%s (x)' % label)
- plt.errorbar(T,dy['std'],yerr=dy['std_err'],fmt='o',label='%s (y)' % label)
- plt.errorbar(T,dz['std'],yerr=dz['std_err'],fmt='o',label='%s (z)' % label)
- plt.xlabel("Kinetic Energy (MeV)")
- plt.ylabel("Position resolution (cm)")
- plt.title("Position Resolution")
- plt.ylim((0,plt.gca().get_ylim()[1]))
- plt.legend()
-
- plt.figure(5)
- plt.errorbar(T,theta['std'],yerr=theta['std_err'],fmt='o',label=label)
- plt.xlabel("Kinetic Energy (MeV)")
- plt.ylabel("Angular resolution (deg)")
- plt.title("Angular Resolution")
- plt.ylim((0,plt.gca().get_ylim()[1]))
- plt.legend()
-
- plt.figure(6)
- plt.scatter(events['Te'],events['ratio'],label=label)
- plt.xlabel("Reconstructed Electron Energy (MeV)")
- plt.ylabel(r"Log Likelihood Ratio (e/$\mu$)")
- plt.title("Log Likelihood Ratio vs Reconstructed Electron Energy")
- plt.legend()
+ marker = markers.next()
+
+ plt.figure(1,figsize=FIGSIZE)
+ plt.errorbar(T,dT['median']*100/T,yerr=dT['median_err']*100/T,fmt=marker,label=label)
+
+ plt.figure(2,figsize=FIGSIZE)
+ plt.errorbar(T,dT['iqr_std']*100/T,yerr=dT['iqr_std_err']*100/T,fmt=marker,label=label)
+
+ ax3[0].errorbar(T,dx['median'],yerr=dx['median_err'],fmt=marker,label=label)
+ ax3[1].errorbar(T,dy['median'],yerr=dy['median_err'],fmt=marker,label=label)
+ ax3[2].errorbar(T,dz['median'],yerr=dz['median_err'],fmt=marker,label=label)
+
+ ax4[0].errorbar(T,dx['iqr_std'],yerr=dx['iqr_std_err'],fmt=marker,label=label)
+ ax4[1].errorbar(T,dy['iqr_std'],yerr=dy['iqr_std_err'],fmt=marker,label=label)
+ ax4[2].errorbar(T,dz['iqr_std'],yerr=dz['iqr_std_err'],fmt=marker,label=label)
+
+ plt.figure(5,figsize=FIGSIZE)
+ plt.errorbar(T,theta['iqr_std'],yerr=theta['iqr_std_err'],fmt=marker,label=label)
+
+ plt.figure(6,figsize=FIGSIZE)
+ plt.scatter(events['Te'],events['ratio'],marker=marker,label=label)
+
+ fig = plt.figure(1)
+ despine(fig,trim=True)
+ plt.xlabel("Kinetic Energy (MeV)")
+ plt.ylabel(r"Energy bias (\%)")
+ plt.title("Energy Bias")
+ plt.legend()
+ plt.tight_layout()
+
+ fig = plt.figure(2)
+ despine(fig,trim=True)
+ plt.xlabel("Kinetic Energy (MeV)")
+ plt.ylabel(r"Energy resolution (\%)")
+ plt.title("Energy Resolution")
+ plt.legend()
+ plt.tight_layout()
+
+ ax3[0].set_ylabel("X")
+ ax3[0].set_ylim((-5,5))
+ ax3[1].set_ylabel("Y")
+ ax3[1].set_ylim((-5,5))
+ ax3[2].set_xlabel("Kinetic Energy (MeV)")
+ ax3[2].set_ylabel("Z")
+ ax3[2].set_ylim((-5,5))
+ despine(ax=ax3[0],trim=True)
+ despine(ax=ax3[1],trim=True)
+ despine(ax=ax3[2],trim=True)
+ fig3.suptitle("Position Bias (cm)")
+ h,l = ax3[0].get_legend_handles_labels()
+ fig3.legend(h,l,loc='upper right')
+ fig3.subplots_adjust(right=0.75)
+ fig3.tight_layout()
+ fig3.subplots_adjust(top=0.9)
+
+ ax4[0].set_ylabel("X")
+ ax4[0].set_ylim((0,ax4[0].get_ylim()[1]))
+ ax4[1].set_ylabel("Y")
+ ax4[1].set_ylim((0,ax4[1].get_ylim()[1]))
+ ax4[2].set_xlabel("Kinetic Energy (MeV)")
+ ax4[2].set_ylabel("Z")
+ ax4[2].set_ylim((0,ax4[2].get_ylim()[1]))
+ despine(ax=ax4[0],trim=True)
+ despine(ax=ax4[1],trim=True)
+ despine(ax=ax4[2],trim=True)
+ fig4.suptitle("Position Resolution (cm)")
+ h,l = ax4[0].get_legend_handles_labels()
+ fig4.legend(h,l,loc='upper right')
+ fig4.subplots_adjust(right=0.75)
+ fig4.tight_layout()
+ fig4.subplots_adjust(top=0.9)
+
+ fig = plt.figure(5)
+ despine(fig,trim=True)
+ plt.xlabel("Kinetic Energy (MeV)")
+ plt.ylabel("Angular resolution (deg)")
+ plt.title("Angular Resolution")
+ plt.ylim((0,plt.gca().get_ylim()[1]))
+ plt.legend()
+ plt.tight_layout()
+
+ fig = plt.figure(6)
+ plt.xticks(range(0,1250,100))
+ plt.hlines(0,0,1200,color='k',linestyles='--',alpha=0.5)
+ despine(fig,trim=True)
+ plt.xlabel("Reconstructed Electron Energy (MeV)")
+ plt.ylabel(r"Log Likelihood Ratio (e/$\mu$)")
+ plt.title("Log Likelihood Ratio vs Reconstructed Electron Energy")
+ plt.legend()
+ plt.tight_layout()
plt.show()