From 0a4e568aa15fbf55dd46a2f9e232f61d93d5ca2f Mon Sep 17 00:00:00 2001
From: Tooba Ali <alit1@myumanitoba.ca>
Date: Mon, 24 Apr 2023 21:11:59 +0000
Subject: [PATCH] Update benchmarks/dis/analysis/truth_reconstruction.py
---
.../dis/analysis/truth_reconstruction.py | 407 ++++++++----------
1 file changed, 170 insertions(+), 237 deletions(-)
diff --git a/benchmarks/dis/analysis/truth_reconstruction.py b/benchmarks/dis/analysis/truth_reconstruction.py
index 7d24cf7d..ba79e47a 100644
--- a/benchmarks/dis/analysis/truth_reconstruction.py
+++ b/benchmarks/dis/analysis/truth_reconstruction.py
@@ -86,22 +86,17 @@ if Nevents == 100:
ssize = 1
else:
ssize = 0.01
+text_size = 8
#Particle type for Single events
particle = config.split('-')[0].strip()
particle_dict = {'e':[boolean_electron,'Electrons'],'pi':[boolean_pion,'Pions']}
-
-def error_bars(plot_x, plot_y, x_range):
- fig = plt.figure()
- fig.set_figwidth(15)
- fig.set_figheight(6)
- x_axis_range = x_range
+def error_bars(plot_x, plot_y, x_axis_range):
if i == 0 or title == 'difference vs momentum':
xbins = np.geomspace(x_axis_range[0],x_axis_range[-1],12)
else:
xbins = 11
- plt.xlim(x_axis_range)
if np.any(plot_x):
plot_x, plot_y = zip(*sorted(zip(plot_x, plot_y)))
n, xedges = np.histogram(plot_x, bins=xbins, range = x_axis_range)
@@ -114,19 +109,13 @@ def error_bars(plot_x, plot_y, x_range):
start = start+n[index]
sum_yy, xedges = np.histogram(plot_x, bins=xbins, range = x_axis_range, weights=(plot_y-mean_list)**2)
std = np.sqrt(sum_yy/(n-1))
- no_nan_std = std[np.invert(np.isnan(std))]
+ no_nan_std = std[np.invert(np.logical_or(np.isnan(std),std == 0))]
if np.any(no_nan_std):
min_std = no_nan_std.min()
- if i == 0 and title == 'ratio':
- plt.ylim(1-(min_std*10),1+(min_std*10))
- else:
- plt.ylim(0-(min_std*10),0+(min_std*10))
else:
min_std = np.nan
- plt.scatter(plot_x, plot_y, s = ssize)
bin_Midpoint = (xedges[1:] + xedges[:-1])/2
bin_HalfWidth = (xedges[1:] - xedges[:-1])/2
- plt.errorbar(bin_Midpoint, mean, yerr=std, xerr=bin_HalfWidth ,fmt='None', ecolor = 'orange', elinewidth = 1)
return bin_Midpoint, mean, bin_HalfWidth, std, min_std
def same_length_lists(plot_x, plot_y):
@@ -137,7 +126,7 @@ def same_length_lists(plot_x, plot_y):
else:
F_boolean = np.ones_like(plot_x) == 1
return F_boolean
-
+
for i in range(len(MC_list)): #Repeat the following steps for each variable (momentum,theta,phi,eta)
MCparts = MC_list[i] #MCParticles events to be plotted on x-axis
RCparts = RC_list[i] #ReconstructedParticles events
@@ -157,6 +146,7 @@ for i in range(len(MC_list)): #Repeat the following steps for each variable (mom
Y_plot = list(np.zeros(len(X_list)))
Y_error = list(np.zeros(len(X_list)))
+
####################################################################################################
#Ratio
####################################################################################################
@@ -171,125 +161,84 @@ for i in range(len(MC_list)): #Repeat the following steps for each variable (mom
ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
X_plot[j] = X_s
Y_plot[j] = ratio
-
- fig = plt.figure()
- gs = fig.add_gridspec(3, 2, wspace=0)
- (ax1, ax2), (ax3, ax4),(ax5, ax6) = gs.subplots(sharex=True, sharey=True)
- # fig.suptitle('')
if i == 1: # for theta
X_plot[0],X_plot[1],X_plot[2],X_plot[3],X_plot[4],X_plot[5] = -X_plot[0],-X_plot[1],-X_plot[2],-X_plot[3],-X_plot[4],-X_plot[5]
- ax1.scatter(X_plot[0], Y_plot[0], s = ssize)
- ax2.scatter(X_plot[1], Y_plot[1], s = ssize)
- ax3.scatter(X_plot[2], Y_plot[2], s = ssize)
- ax4.scatter(X_plot[3], Y_plot[3], s = ssize)
- ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
- ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
-
- if i == 0: # for momentum
- ax1.set_ylabel('rc/mc') #ratio
- ax3.set_ylabel('rc/mc')
- ax5.set_ylabel('rc/mc')
- title ='ratio'
- ax1.set_yscale('log')
- ax1.set_xscale('log')
- else: # for angles
- ax1.set_ylabel('rc-mc') #difference
- ax3.set_ylabel('rc-mc')
- ax5.set_ylabel('rc-mc')
- title ='difference'
- ax2.set_title('Pions')
- ax3.set_title('Protons')
- ax4.set_title('Electrons')
- ax5.set_title('Neutrons')
- ax6.set_title('Photons')
-
- if i == 3: #Eta
- ax1.set_xlim(-5.5,5.5)
- if i == 1: #Theta
- ax1.set_xlim(-np.pi,0)
- ax5.set_xlabel('- %s mc'%(title_list[i]))
- ax6.set_xlabel('- %s mc'%(title_list[i]))
- x_range = [0,np.pi]
- else:
- ax5.set_xlabel('%s mc'%(title_list[i]))
- ax6.set_xlabel('%s mc'%(title_list[i]))
- x_range = list(ax1.get_xlim())
- fig.set_figwidth(20)
- fig.set_figheight(10)
- ax1.set_title('%s %s %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
- plt.savefig(os.path.join(r_path, '%s_%s_%s.png' % (title_list_n[i],title,config)))
- plt.close()
+
+ particle_nlist = ['All','Pions','Protons','Electrons','Neutrons','Photons']
+ for iterate in [0,1]:
+ fig = plt.figure()
+ gs = fig.add_gridspec(3, 2, wspace=0)
+ (ax1, ax2), (ax3, ax4),(ax5, ax6) = gs.subplots(sharex=True, sharey=True)
+ ax1.scatter(X_plot[0], Y_plot[0], s = ssize)
+ ax2.scatter(X_plot[1], Y_plot[1], s = ssize)
+ ax3.scatter(X_plot[2], Y_plot[2], s = ssize)
+ ax4.scatter(X_plot[3], Y_plot[3], s = ssize)
+ ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
+ ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
+
+ if i == 0: # for momentum
+ ax1.set_ylabel('rc/mc') #ratio
+ ax3.set_ylabel('rc/mc')
+ ax5.set_ylabel('rc/mc')
+ title ='ratio'
+ ax1.set_yscale('log')
+ ax1.set_xscale('log')
+ else: # for angles
+ ax1.set_ylabel('rc-mc') #difference
+ ax3.set_ylabel('rc-mc')
+ ax5.set_ylabel('rc-mc')
+ title ='difference'
+ ax2.set_title('Pions')
+ ax3.set_title('Protons')
+ ax4.set_title('Electrons')
+ ax5.set_title('Neutrons')
+ ax6.set_title('Photons')
+ if i == 3: #Eta
+ ax1.set_xlim(-5.5,5.5)
+ if i == 1: #Theta
+ ax1.set_xlim(-np.pi,0)
+ ax5.set_xlabel('- %s mc'%(title_list[i]))
+ ax6.set_xlabel('- %s mc'%(title_list[i]))
+ x_range = [0,np.pi]
+ else:
+ ax5.set_xlabel('%s mc'%(title_list[i]))
+ ax6.set_xlabel('%s mc'%(title_list[i]))
+ x_range = list(ax1.get_xlim())
+ if i == 0:
+ momentum_range = x_range
+ fig.set_figwidth(20)
+ fig.set_figheight(10)
- particle_nlist = ['All','Pions','Protons','Electrons','Neutrons','Photons']
- for j in range(len(X_list)):#Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
- if i == 1: #theta
- Y_error[j] = error_bars(X_plot[j], Y_plot[j], [-np.pi,0])
+ if iterate == 0:
+ for j in range(len(X_list)):#Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
+ if i == 1: #theta
+ Y_error[j] = error_bars(X_plot[j], Y_plot[j], [-np.pi,0])
+ else:
+ Y_error[j] = error_bars(X_plot[j], Y_plot[j], x_range)
+ ax1.errorbar(Y_error[0][0], Y_error[0][1], yerr=Y_error[0][3], xerr=Y_error[0][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax2.errorbar(Y_error[1][0], Y_error[1][1], yerr=Y_error[1][3], xerr=Y_error[1][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax3.errorbar(Y_error[2][0], Y_error[2][1], yerr=Y_error[2][3], xerr=Y_error[2][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax4.errorbar(Y_error[3][0], Y_error[3][1], yerr=Y_error[3][3], xerr=Y_error[3][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax5.errorbar(Y_error[4][0], Y_error[4][1], yerr=Y_error[4][3], xerr=Y_error[4][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax6.errorbar(Y_error[5][0], Y_error[5][1], yerr=Y_error[5][3], xerr=Y_error[5][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+
+ if i == 0: # for momentum
+ ax1.set_ylim(1-(Y_error[0][4]*10),1+(Y_error[0][4]*10))
+ center = 1
+ else: # for angles
+ ax1.set_ylim(0-(Y_error[0][4]*10),0+(Y_error[0][4]*10))
+ center = 0
+ for each_bin in range(len(Y_error[0][0])):
+ ax1.text(x=Y_error[0][0][each_bin],y=center + Y_error[0][4]*7, s= '\u03BC = %.3f\n\u03C3 = %.3f' % (Y_error[0][1][each_bin],Y_error[0][3][each_bin]),size=text_size,horizontalalignment='center',verticalalignment='top')
+
+ ax1.set_title('%s %s with error bars %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
+ plt.savefig(os.path.join(r_path, '%s_%s_2_error_%s.png' % (title_list_n[i],title,config)))
else:
- Y_error[j] = error_bars(X_plot[j], Y_plot[j], x_range)
- plt.title('%s %s %s\n %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,particle_nlist[j],config,Nevents,Dconfig))
- # plt.savefig(os.path.join(r_path, 'error_bars_%s_%s_%s_%s.png' % (title_list[i],title,particle_nlist[j],config)))
- if i == 0:
- if np.any(Y_plot[j]):
- plt.yscale('log')
- plt.xscale('log')
- if i == 1:
- plt.xlim(-np.pi,0)
- plt.close()
- fig = plt.figure()
- gs = fig.add_gridspec(3, 2, wspace=0)
- (ax1, ax2), (ax3, ax4),(ax5, ax6) = gs.subplots(sharex=True, sharey=True)
- # if i == 1: # for theta
- # X_plot[0],X_plot[1],X_plot[2],X_plot[3],X_plot[4],X_plot[5] = -X_plot[0],-X_plot[1],-X_plot[2],-X_plot[3],-X_plot[4],-X_plot[5]
- ax1.scatter(X_plot[0], Y_plot[0], s = ssize)
- ax2.scatter(X_plot[1], Y_plot[1], s = ssize)
- ax3.scatter(X_plot[2], Y_plot[2], s = ssize)
- ax4.scatter(X_plot[3], Y_plot[3], s = ssize)
- ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
- ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
- ax1.errorbar(Y_error[0][0], Y_error[0][1], yerr=Y_error[0][3], xerr=Y_error[0][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax2.errorbar(Y_error[1][0], Y_error[1][1], yerr=Y_error[1][3], xerr=Y_error[1][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax3.errorbar(Y_error[2][0], Y_error[2][1], yerr=Y_error[2][3], xerr=Y_error[2][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax4.errorbar(Y_error[3][0], Y_error[3][1], yerr=Y_error[3][3], xerr=Y_error[3][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax5.errorbar(Y_error[4][0], Y_error[4][1], yerr=Y_error[4][3], xerr=Y_error[4][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax6.errorbar(Y_error[5][0], Y_error[5][1], yerr=Y_error[5][3], xerr=Y_error[5][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- if i == 0: # for momentum
- ax1.set_ylabel('rc/mc') #ratio
- ax3.set_ylabel('rc/mc')
- ax5.set_ylabel('rc/mc')
- title ='ratio'
- ax1.set_yscale('log')
- ax1.set_xscale('log')
- ax1.set_ylim(1-(Y_error[0][4]*10),1+(Y_error[0][4]*10))
- else: # for angles
- ax1.set_ylabel('rc-mc') #difference
- ax3.set_ylabel('rc-mc')
- ax5.set_ylabel('rc-mc')
- title ='difference'
- ax1.set_ylim(0-(Y_error[0][4]*10),0+(Y_error[0][4]*10))
- ax2.set_title('Pions')
- ax3.set_title('Protons')
- ax4.set_title('Electrons')
- ax5.set_title('Neutrons')
- ax6.set_title('Photons')
- if i == 3: #Eta
- ax1.set_xlim(-5.5,5.5)
- if i == 1: #Theta
- ax1.set_xlim(-np.pi,0)
- ax5.set_xlabel('- %s mc'%(title_list[i]))
- ax6.set_xlabel('- %s mc'%(title_list[i]))
- x_range = [0,np.pi]
- else:
- ax5.set_xlabel('%s mc'%(title_list[i]))
- ax6.set_xlabel('%s mc'%(title_list[i]))
- x_range = list(ax1.get_xlim())
- fig.set_figwidth(20)
- fig.set_figheight(10)
- ax1.set_title('%s %s with error bars %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
- plt.savefig(os.path.join(r_path, '%s_error_%s_%s.png' % (title_list_n[i],title,config)))
- plt.close()
-
-
+ ax1.set_title('%s %s %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
+ plt.savefig(os.path.join(r_path, '%s_%s_1_%s.png' % (title_list_n[i],title,config)))
+
+
###################################################################################################
#Ratio vs momentum
###################################################################################################
@@ -307,98 +256,65 @@ for i in range(len(MC_list)): #Repeat the following steps for each variable (mom
X_plot[j] = M_s
Y_plot[j] = ratio
- fig = plt.figure()
- gs = fig.add_gridspec(3, 2, wspace=0)
- (ax1, ax2), (ax3, ax4),(ax5, ax6) = gs.subplots(sharex=True, sharey=True)
- ax1.scatter(X_plot[0], Y_plot[0], s = ssize)
- ax2.scatter(X_plot[1], Y_plot[1], s = ssize)
- ax3.scatter(X_plot[2], Y_plot[2], s = ssize)
- ax4.scatter(X_plot[3], Y_plot[3], s = ssize)
- ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
- ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
- ax1.set_xscale('log')
- ax1.set_ylabel('rc-mc')
- ax3.set_ylabel('rc-mc')
- ax5.set_ylabel('rc-mc')
- ax5.set_xlabel('Momentum mc')
- ax6.set_xlabel('Momentum mc')
- ax2.set_title('Pions')
- ax3.set_title('Protons')
- ax4.set_title('Electrons')
- ax5.set_title('Neutrons')
- ax6.set_title('Photons')
- fig.set_figwidth(20)
- fig.set_figheight(10)
- x_range = list(ax1.get_xlim())
- ax1.set_title('%s Difference Vs Momentum %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],config,Nevents,Dconfig))
- plt.savefig(os.path.join(r_path, '%s_difference_vs_momentum_%s.png' % (title_list_n[i],config)))
- plt.close()
-
title ='difference vs momentum'
particle_nlist = ['All','Pions','Protons','Electrons','Neutrons','Photons']
- for j in range(len(X_list)):#Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
- Y_error[j] = error_bars(X_plot[j], Y_plot[j], x_range)
- plt.title('%s %s %s\n %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,particle_nlist[j],config,Nevents,Dconfig))
- # plt.savefig(os.path.join(r_path, 'error_bars_%s_%s_%s_%s.png' % (title_list[i],title,particle_nlist[j],config)))
- plt.xscale('log')
- plt.show()
- plt.close()
- fig = plt.figure()
- gs = fig.add_gridspec(3, 2, wspace=0)
- (ax1, ax2), (ax3, ax4),(ax5, ax6) = gs.subplots(sharex=True, sharey=True)
- ax1.scatter(X_plot[0], Y_plot[0], s = ssize)
- ax2.scatter(X_plot[1], Y_plot[1], s = ssize)
- ax3.scatter(X_plot[2], Y_plot[2], s = ssize)
- ax4.scatter(X_plot[3], Y_plot[3], s = ssize)
- ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
- ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
- ax1.errorbar(Y_error[0][0], Y_error[0][1], yerr=Y_error[0][3], xerr=Y_error[0][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax2.errorbar(Y_error[1][0], Y_error[1][1], yerr=Y_error[1][3], xerr=Y_error[1][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax3.errorbar(Y_error[2][0], Y_error[2][1], yerr=Y_error[2][3], xerr=Y_error[2][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax4.errorbar(Y_error[3][0], Y_error[3][1], yerr=Y_error[3][3], xerr=Y_error[3][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax5.errorbar(Y_error[4][0], Y_error[4][1], yerr=Y_error[4][3], xerr=Y_error[4][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- ax6.errorbar(Y_error[5][0], Y_error[5][1], yerr=Y_error[5][3], xerr=Y_error[5][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
- # for angles
- ax1.set_ylabel('rc-mc') #difference
- ax3.set_ylabel('rc-mc')
- ax5.set_ylabel('rc-mc')
- ax5.set_xlabel('Momentum mc')
- ax6.set_xlabel('Momentum mc')
- ax1.set_xscale('log')
- ax1.set_ylim(0-(Y_error[0][4]*10),0+(Y_error[0][4]*10))
- ax2.set_title('Pions')
- ax3.set_title('Protons')
- ax4.set_title('Electrons')
- ax5.set_title('Neutrons')
- ax6.set_title('Photons')
-
- fig.set_figwidth(20)
- fig.set_figheight(10)
- ax1.set_title('%s %s with error bars %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
- plt.savefig(os.path.join(r_path, '%s_difference_vs_momentum_error_%s.png' % (title_list_n[i],config)))
- plt.close()
-
+ for iterate in [0,1]:
+ fig = plt.figure()
+ gs = fig.add_gridspec(3, 2, wspace=0)
+ (ax1, ax2), (ax3, ax4),(ax5, ax6) = gs.subplots(sharex=True, sharey=True)
+ ax1.scatter(X_plot[0], Y_plot[0], s = ssize)
+ ax2.scatter(X_plot[1], Y_plot[1], s = ssize)
+ ax3.scatter(X_plot[2], Y_plot[2], s = ssize)
+ ax4.scatter(X_plot[3], Y_plot[3], s = ssize)
+ ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
+ ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
+ ax1.set_xscale('log')
+ ax1.set_ylabel('rc-mc')
+ ax3.set_ylabel('rc-mc')
+ ax5.set_ylabel('rc-mc')
+ ax5.set_xlabel('Momentum mc')
+ ax6.set_xlabel('Momentum mc')
+ ax2.set_title('Pions')
+ ax3.set_title('Protons')
+ ax4.set_title('Electrons')
+ ax5.set_title('Neutrons')
+ ax6.set_title('Photons')
+ fig.set_figwidth(20)
+ fig.set_figheight(10)
+ if iterate == 0:
+ for j in range(len(X_list)):#Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
+ Y_error[j] = error_bars(X_plot[j], Y_plot[j], momentum_range)
+ ax1.errorbar(Y_error[0][0], Y_error[0][1], yerr=Y_error[0][3], xerr=Y_error[0][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax2.errorbar(Y_error[1][0], Y_error[1][1], yerr=Y_error[1][3], xerr=Y_error[1][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax3.errorbar(Y_error[2][0], Y_error[2][1], yerr=Y_error[2][3], xerr=Y_error[2][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax4.errorbar(Y_error[3][0], Y_error[3][1], yerr=Y_error[3][3], xerr=Y_error[3][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax5.errorbar(Y_error[4][0], Y_error[4][1], yerr=Y_error[4][3], xerr=Y_error[4][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax6.errorbar(Y_error[5][0], Y_error[5][1], yerr=Y_error[5][3], xerr=Y_error[5][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax1.set_ylim(0-(Y_error[0][4]*10),0+(Y_error[0][4]*10))
+ for each_bin in range(len(Y_error[0][0])):
+ ax1.text(x=Y_error[0][0][each_bin],y=0 + Y_error[0][4]*7,
+ s= '\u03BC = %.3f\n\u03C3 = %.3f' % (Y_error[0][1][each_bin],Y_error[0][3][each_bin]),size=text_size,horizontalalignment='center',verticalalignment='top')
+ ax1.set_title('%s %s with error bars %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
+ plt.savefig(os.path.join(r_path, '%s_difference_vs_momentum_2_error_%s.png' % (title_list_n[i],config)))
+ else:
+ ax1.set_title('%s Difference Vs Momentum %s %s events\n DETECTOR_CONFIG: %s'%(title_list[i],config,Nevents,Dconfig))
+ plt.savefig(os.path.join(r_path, '%s_difference_vs_momentum_1_%s.png' % (title_list_n[i],config)))
+
###################################################################################################
#Correlation
###################################################################################################
#Repeat the following steps for each variable (momentum,theta,phi,eta)
- X_len = ak.count(MCparts,axis=None)
- Y_len = ak.count(RCparts,axis=None)
- if X_len > Y_len:
- F_boolean = np.ones_like(RCparts) == 1
- else:
- F_boolean = np.ones_like(MCparts) == 1
- X_s = np.array(ak.flatten(MCparts[F_boolean]))
+ F_boolean = same_length_lists(MCparts,RCparts)
+ X_s = np.array(ak.flatten(MCparts[F_boolean])) #Filtered lists
Y_s = np.array(ak.flatten(RCparts[F_boolean]))
-
+
#Histogram
if i == 0 and particle in particle_dict.keys(): #Momentum in Single events
- h, xedges, yedges, image = plt.hist2d(x=X_s,y= Y_s, bins = 11)
+ h, xedges, yedges = np.histogram2d(x=X_s,y= Y_s, bins = 11)
else:
- h, xedges, yedges, image = plt.hist2d(x=X_s,y= Y_s, bins = 11, range = [x_range,x_range])
- plt.close()
+ h, xedges, yedges = np.histogram2d(x=X_s,y= Y_s, bins = 11, range = [x_range,x_range])
col_sum = ak.sum(h,axis=-1) #number of events in each (verticle) column
norm_h = [] #norm_h is the normalized matrix
@@ -418,15 +334,14 @@ for i in range(len(MC_list)): #Repeat the following steps for each variable (mom
else:
axs[0].hist2d(x=X_s,y=Y_s, bins = 11,range = [x_range,x_range])
mplhep.hist2dplot(H=norm_h,norm=mpl.colors.LogNorm(vmin= 1e-4, vmax= 1),labels=norm_h_text, xbins = xedges, ybins = yedges, ax=axs[1])
-
axs[0].set_title('%s Histogram'%(title_list[i]))
- axs[0].set_xlabel('%s_mc'%(title_list[i]))
axs[0].set_ylabel('%s_rc'%(title_list[i]))
- axs[1].set_xlabel('%s_mc'%(title_list[i]))
axs[1].set_ylabel('%s_rc'%(title_list[i]))
+ axs[0].set_xlabel('%s_mc'%(title_list[i]))
+ axs[1].set_xlabel('%s_mc'%(title_list[i]))
axs[1].set_title('%s Correlation'%(title_list[i]))
fig.suptitle('%s %s events\n DETECTOR_CONFIG: %s'%(config,Nevents,Dconfig))
- # plt.savefig(os.path.join(r_path, '%s%s_%s.png' % (title_list_n[i],args.config.split('_epic_')[1].strip(),config)))
+ plt.savefig(os.path.join(r_path, '%s%s_%s.png' % (title_list_n[i],args.config.split('_epic_')[1].strip(),config)))
###################################################################################################
@@ -439,45 +354,68 @@ def particle_plots(boolean_particle):
theta_rc_fil = ak.Array(theta_rc[recID][booll])[boolean_particle]
phi_mc_fil = ak.Array(phi_mc[simID][booll])[boolean_particle]
phi_rc_fil = ak.Array(phi_rc[recID][booll])[boolean_particle]
-
- theta_mc_fil_len = ak.count(theta_mc_fil,axis=None)
- theta_rc_fil_len = ak.count(theta_rc_fil,axis=None)
- if theta_mc_fil_len > theta_rc_fil_len:
- F_boolean = np.ones_like(theta_rc_fil) == 1
- else:
- F_boolean = np.ones_like(theta_mc_fil) == 1
+
+ F_boolean = same_length_lists(theta_mc_fil, theta_rc_fil)
#filtered lists w.r.t length
theta_mc_F = np.array(ak.flatten(theta_mc_fil[F_boolean]))
theta_rc_F = np.array(ak.flatten(theta_rc_fil[F_boolean]))
phi_mc_F = np.array(ak.flatten(phi_mc_fil[F_boolean]))
phi_rc_F = np.array(ak.flatten(phi_rc_fil[F_boolean]))
ratio = np.array((ak.Array(theta_rc_F)-(ak.Array(theta_mc_F))))
-
+ x_range = [0,np.pi]
+ Y_error = [error_bars(theta_mc_F, ratio, x_range),error_bars(theta_rc_F, ratio, x_range)]
fig = plt.figure()
- gs = fig.add_gridspec(2, 2, wspace=0.01)
- (ax1, ax2), (ax3, ax4) = gs.subplots(sharex=True, sharey=True)
+ gs = fig.add_gridspec(3, 2, wspace=0, hspace = 0.3)
+ (ax1, ax2), (ax3, ax4), (ax5, ax6) = gs.subplots(sharex=True, sharey='row')
ax1.scatter(-theta_mc_F, ratio, s = ssize)
ax2.scatter(-theta_rc_F, ratio, s = ssize)
- ax3.scatter(-theta_mc_F, phi_mc_F, s = ssize)
- ax4.scatter(-theta_rc_F, phi_rc_F, s = ssize)
- ax1.set_ylabel('rc-mc')
- ax2.set_ylabel('rc-mc')
- ax3.set_ylabel('Phi mc')
- ax4.set_ylabel('Phi rc')
+ ax3.scatter(-theta_mc_F, ratio, s = ssize)
+ ax4.scatter(-theta_rc_F, ratio, s = ssize)
+ ax5.scatter(-theta_mc_F, phi_mc_F, s = ssize)
+ ax6.scatter(-theta_rc_F, phi_rc_F, s = ssize)
+ ax1.set_ylabel('Theta rc-mc')
+ ax2.set_ylabel('Theta rc-mc')
+ ax3.set_ylabel('Theta rc-mc')
+ ax4.set_ylabel('Theta rc-mc')
+ ax5.set_ylabel('Phi mc')
+ ax6.set_ylabel('Phi rc')
ax1.set_xlabel('- Theta mc')
ax2.set_xlabel('- Theta rc')
ax3.set_xlabel('- Theta mc')
ax4.set_xlabel('- Theta rc')
+ ax5.set_xlabel('- Theta mc')
+ ax6.set_xlabel('- Theta rc')
+ ax1.set_title('Zoom-in')
+ ax2.set_title('Zoom-in')
+ ax3.set_title('Zoom-out')
+ ax4.set_title('Zoom-out')
+ ax5.set_title('Phi vs Theta mc')
+ ax6.set_title('Phi vs Theta rc')
+ ax1.errorbar(-Y_error[0][0], Y_error[0][1], yerr=Y_error[0][3], xerr=Y_error[0][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax2.errorbar(-Y_error[1][0], Y_error[1][1], yerr=Y_error[1][3], xerr=Y_error[1][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax3.errorbar(-Y_error[0][0], Y_error[0][1], yerr=Y_error[0][3], xerr=Y_error[0][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ ax4.errorbar(-Y_error[1][0], Y_error[1][1], yerr=Y_error[1][3], xerr=Y_error[1][2] ,fmt='None', ecolor = 'orange', elinewidth = 1)
+ y_limits = ax3.get_ylim()
+ for each_bin in range(len(Y_error[0][0])):
+ if not np.isnan(Y_error[0][1][each_bin]):
+ ax3.text(x=-Y_error[0][0][each_bin],y=y_limits[1],
+ s= '\u03BC = %.3f\n\u03C3 = %.3f' % (Y_error[0][1][each_bin],Y_error[0][3][each_bin]),size=text_size,horizontalalignment='center',verticalalignment='top')
+ if not np.isnan(Y_error[1][1][each_bin]):
+ ax4.text(x=-Y_error[1][0][each_bin],y=y_limits[1],
+ s= '\u03BC = %.3f\n\u03C3 = %.3f' % (Y_error[1][1][each_bin],Y_error[1][3][each_bin]),size=text_size,horizontalalignment='center',verticalalignment='top')
+ if not np.isnan(Y_error[0][4]):
+ ax1.set_ylim(0-(Y_error[1][4]*10),0+(Y_error[1][4]*10))
+ ax2.set_ylim(0-(Y_error[1][4]*10),0+(Y_error[1][4]*10))
fig.set_figwidth(20)
fig.set_figheight(10)
-
+title ='difference'
if particle in particle_dict.keys():
boolean_particle = particle_dict[particle][0]
particle_name = particle_dict[particle][1]
particle_plots(boolean_particle)
plt.suptitle('%s in %s %s events\n DETECTOR_CONFIG: %s'%(particle_name,config,Nevents,Dconfig))
- # plt.savefig(os.path.join(r_path, '%s_%s.png' % (particle_name_n[particle_name],config)))
+ plt.savefig(os.path.join(r_path, '%s_%s.png' % (particle_name_n[particle_name],config)))
else:
for i in [[boolean_photon,'Photons'],[boolean_electron,'Electrons'],[boolean_pion,'Pions']]:
boolean_particle = i[0]
@@ -485,9 +423,4 @@ else:
particle_plots(boolean_particle)
plt.suptitle('%s in %s %s events\n DETECTOR_CONFIG: %s'%(particle_name,config,Nevents,Dconfig))
- # plt.savefig(os.path.join(r_path, '%s_%s.png' % (particle_name_n[particle_name],config)))
-
-
-
-
-
+ plt.savefig(os.path.join(r_path, '%s_%s.png' % (particle_name_n[particle_name],config)))
--
GitLab