dis/analysis/truth_reconstruction.py

0001 import os
0002 import numpy as np
0003 import uproot as ur
0004 import awkward as ak
0005 import matplotlib.pyplot as plt
0006 import matplotlib as mpl
0007 import mplhep
0008 import argparse
0009
0010 parser = argparse.ArgumentParser()
0011 parser.add_argument('--rec_file', type=str, help='Reconstructed track file.')
0012 parser.add_argument('--config', type=str, help='Momentum configuration.')
0013 parser.add_argument('--nevents', type=float, help='Number of events to process.')
0014 parser.add_argument('--results_path', type=str, help='Output directory.')
0015 args = parser.parse_args()
0016 kwargs = vars(args)
0017
0018 rec_file = args.rec_file
0019 Nevents = int(args.nevents)
0020 r_path = args.results_path + '/truth_reconstruction/' #Path for output figures and file.
0021 Dconfig = 'epic' + args.config.split('_epic')[1].strip() #Detector config
0022 config = args.config.split('_epic')[0].strip()
0023
0024 for array in ur.iterate(rec_file + ':events',['MCParticles/MCParticles.generatorStatus',
0025                                         'MCParticles/MCParticles.PDG',
0026                                         'MCParticles/MCParticles.momentum.x',
0027                                         'MCParticles/MCParticles.momentum.y',
0028                                         'MCParticles/MCParticles.momentum.z',
0029                                         'ReconstructedChargedParticles/ReconstructedChargedParticles.PDG',
0030                                         'ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.x',
0031                                         'ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.y',
0032                                         'ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.z',
0033                                         'ReconstructedChargedParticleAssociations/ReconstructedChargedParticleAssociations.simID',
0034                                         'ReconstructedChargedParticleAssociations/ReconstructedChargedParticleAssociations.recID'],step_size=Nevents):
0035     PDG_mc = array['MCParticles/MCParticles.PDG']  #Monte Carlo (MC) particle numbering scheme.
0036     px_mc = array['MCParticles/MCParticles.momentum.x']
0037     py_mc = array['MCParticles/MCParticles.momentum.y']
0038     pz_mc = array['MCParticles/MCParticles.momentum.z']
0039     PDG_rc = array['ReconstructedChargedParticles/ReconstructedChargedParticles.PDG']
0040     px_rc = array['ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.x']
0041     py_rc = array['ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.y']
0042     pz_rc = array['ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.z']
0043     simID = array['ReconstructedChargedParticleAssociations/ReconstructedChargedParticleAssociations.simID']
0044     recID = array['ReconstructedChargedParticleAssociations/ReconstructedChargedParticleAssociations.recID']
0045     #SimID and recID contain the indices of the MCParticles and ReconstructedParticles entry for that event.
0046
0047 ### MCParticles Variables
0048 momentum_mc = np.sqrt(((px_mc**2)+(py_mc**2)+(pz_mc**2)))
0049 theta_mc = np.arctan2(np.sqrt(px_mc**2+py_mc**2), pz_mc)
0050 phi_mc = np.arctan2(py_mc, px_mc)
0051 ### ReconstructedParticles Variables
0052 momentum_rc = np.sqrt(((px_rc**2)+(py_rc**2)+(pz_rc**2)))
0053 theta_rc = np.arctan2(np.sqrt(px_rc**2+py_rc**2), pz_rc)
0054 phi_rc = np.arctan2(py_rc, px_rc)
0055
0056 booll = (PDG_mc[simID])==(PDG_rc[recID]) #boolean that allows events where the same particle is reconstructed
0057 boolean_pion = np.logical_or(ak.Array(PDG_mc[simID][booll])==-211, ak.Array(PDG_mc[simID][booll])==+211) #boolean that allows events involving pions
0058 boolean_proton = np.logical_or(ak.Array(PDG_mc[simID][booll])==-2212, ak.Array(PDG_mc[simID][booll])==+2212) #boolean that allows events involving protons
0059 boolean_electron = ak.Array(PDG_mc[simID][booll])==11 #boolean that allows events involving electrons
0060 boolean_neutron = ak.Array(PDG_mc[simID][booll])==2112 #boolean that allows events involving neutrons
0061 boolean_photon = ak.Array(PDG_mc[simID][booll])==22 #boolean that allows events involving photons
0062
0063 ### MCParticles variables list
0064 MC_list = [ak.Array(momentum_mc[simID][booll]),                     #Momentum
0065            ak.Array(theta_mc[simID][booll]),                        #Theta
0066            ak.Array(phi_mc[simID][booll]),                          #Phi
0067            -np.log(np.tan((ak.Array(theta_mc[simID][booll]))/2))]   #Eta
0068 ### ReconstructedParticles variables list
0069 RC_list = [ak.Array(momentum_rc[recID][booll]),                     #Momentum
0070            ak.Array(theta_rc[recID][booll]),                        #Theta
0071            ak.Array(phi_rc[recID][booll]),                          #Phi
0072            -np.log(np.tan((ak.Array(theta_rc[recID][booll]))/2))]   #Eta
0073 title_list = ['Momentum','Theta','Phi','Eta']
0074 ### MC Momentum for different particles list
0075 M_list = [ak.Array(momentum_mc[simID][booll]),
0076           ak.Array(momentum_mc[simID][booll][boolean_pion]),
0077           ak.Array(momentum_mc[simID][booll][boolean_proton]),
0078           ak.Array(momentum_mc[simID][booll][boolean_electron]),
0079           ak.Array(momentum_mc[simID][booll][boolean_neutron]),
0080           ak.Array(momentum_mc[simID][booll][boolean_photon])]
0081
0082 #Marker Size in plots
0083 if Nevents == 100:
0084     ssize = 1
0085 else:
0086     ssize = 0.01
0087 text_size = 8
0088
0089 #Particle type for Single events
0090 particle = config.split('-')[0].strip()
0091 particle_dict = {'e':[boolean_electron,'Electrons'],'pi':[boolean_pion,'Pions']}
0092
0093 def error_bars(plot_x, plot_y, x_axis_range):
0094     if i == 0 or title == 'difference vs momentum':
0095         xbins = np.geomspace(x_axis_range[0],x_axis_range[-1],12)
0096     else:
0097         xbins = 11
0098     if np.any(plot_x):
0099         plot_x, plot_y = zip(*sorted(zip(plot_x, plot_y)))
0100     n, xedges = np.histogram(plot_x, bins=xbins, range = x_axis_range)
0101     sum_y, xedges = np.histogram(plot_x, bins=xbins, range = x_axis_range, weights=plot_y)
0102     mean = sum_y / n
0103     mean_list = np.zeros(len(plot_y))
0104     start = 0
0105     for index in range(len(n)):
0106         mean_list[start:start+n[index]] = mean[index]
0107         start = start+n[index]
0108     sum_yy, xedges = np.histogram(plot_x, bins=xbins, range = x_axis_range, weights=(plot_y-mean_list)**2)
0109     std = np.sqrt(sum_yy/(n-1))
0110     no_nan_std = std[np.invert(np.logical_or(np.isnan(std),std == 0))]
0111     if np.any(no_nan_std):
0112         min_std = no_nan_std.min()
0113     else:
0114         min_std = np.nan
0115     bin_Midpoint = (xedges[1:] + xedges[:-1])/2
0116     bin_HalfWidth = (xedges[1:] - xedges[:-1])/2
0117     return bin_Midpoint, mean, bin_HalfWidth, std, min_std
0118
0119 def same_length_lists(plot_x, plot_y):
0120     X_length = ak.count(plot_x,axis=None)
0121     Y_length = ak.count(plot_y,axis=None)
0122     if X_length > Y_length:
0123         F_boolean = np.ones_like(plot_y) == 1
0124     else:
0125         F_boolean = np.ones_like(plot_x) == 1
0126     return F_boolean
0127
0128 for i in range(len(MC_list)): #Repeat the following steps for each variable (momentum,theta,phi,eta)
0129     MCparts = MC_list[i] #MCParticles events to be plotted on x-axis
0130     RCparts = RC_list[i] #ReconstructedParticles events
0131     X_list = [ak.Array(MCparts),
0132           ak.Array(MCparts[boolean_pion]),
0133           ak.Array(MCparts[boolean_proton]),
0134           ak.Array(MCparts[boolean_electron]),
0135           ak.Array(MCparts[boolean_neutron]),
0136           ak.Array(MCparts[boolean_photon])]
0137     Y_list = [ak.Array(RCparts),
0138           ak.Array(RCparts[boolean_pion]),
0139           ak.Array(RCparts[boolean_proton]),
0140           ak.Array(RCparts[boolean_electron]),
0141           ak.Array(RCparts[boolean_neutron]),
0142           ak.Array(RCparts[boolean_photon])]
0143     X_plot = list(np.zeros(len(X_list)))
0144     Y_plot = list(np.zeros(len(X_list)))
0145     Y_error = list(np.zeros(len(X_list)))
0146
0147
0148 ####################################################################################################
0149     #Ratio
0150 ####################################################################################################
0151
0152     for j in range(len(X_list)): #Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
0153         F_boolean = same_length_lists(X_list[j],Y_list[j])
0154         X_s = np.array(ak.flatten(X_list[j][F_boolean])) #Filtered lists
0155         Y_s = np.array(ak.flatten(Y_list[j][F_boolean]))
0156         if i == 0:   #Momentum
0157             ratio = np.array((ak.Array(Y_s)/ak.Array(X_s)))
0158         else: #Angle difference
0159             ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
0160         X_plot[j] = X_s
0161         Y_plot[j] = ratio
0162     if i == 1:   # for theta
0163         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]
0164
0165     particle_nlist = ['All','Pions','Protons','Electrons','Neutrons','Photons']
0166     for iterate in [0,1]:
0167         fig = plt.figure()
0168         gs = fig.add_gridspec(3, 2, wspace=0)
0169         (ax1, ax2), (ax3, ax4),(ax5, ax6) = gs.subplots(sharex=True, sharey=True)
0170         ax1.scatter(X_plot[0], Y_plot[0], s = ssize)
0171         ax2.scatter(X_plot[1], Y_plot[1], s = ssize)
0172         ax3.scatter(X_plot[2], Y_plot[2], s = ssize)
0173         ax4.scatter(X_plot[3], Y_plot[3], s = ssize)
0174         ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
0175         ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
0176
0177         if i == 0:  # for momentum
0178             ax1.set_ylabel('rc/mc') #ratio
0179             ax3.set_ylabel('rc/mc')
0180             ax5.set_ylabel('rc/mc')
0181             title ='ratio'
0182             ax1.set_yscale('log')
0183             ax1.set_xscale('log')
0184         else:       # for angles
0185             ax1.set_ylabel('rc-mc') #difference
0186             ax3.set_ylabel('rc-mc')
0187             ax5.set_ylabel('rc-mc')
0188             title ='difference'
0189         ax2.set_title('Pions')
0190         ax3.set_title('Protons')
0191         ax4.set_title('Electrons')
0192         ax5.set_title('Neutrons')
0193         ax6.set_title('Photons')
0194
0195         if i == 3: #Eta
0196             ax1.set_xlim(-5.5,5.5)
0197         if i == 1: #Theta
0198             ax1.set_xlim(-np.pi,0)
0199             ax5.set_xlabel('- %s mc'%(title_list[i]))
0200             ax6.set_xlabel('- %s mc'%(title_list[i]))
0201             x_range = [0,np.pi]
0202         else:
0203             ax5.set_xlabel('%s mc'%(title_list[i]))
0204             ax6.set_xlabel('%s mc'%(title_list[i]))
0205             x_range = list(ax1.get_xlim())
0206             if i == 0:
0207                 momentum_range = x_range
0208         fig.set_figwidth(20)
0209         fig.set_figheight(10)
0210
0211         if iterate == 0:
0212             for j in range(len(X_list)):#Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
0213                 if i == 1: #theta
0214                     Y_error[j] = error_bars(X_plot[j], Y_plot[j], [-np.pi,0])
0215                 else:
0216                     Y_error[j] = error_bars(X_plot[j], Y_plot[j], x_range)
0217             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)
0218             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)
0219             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)
0220             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)
0221             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)
0222             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)
0223
0224             if i == 0:  # for momentum
0225                 ax1.set_ylim(1-(Y_error[0][4]*10),1+(Y_error[0][4]*10))
0226                 center = 1
0227             else:       # for angles
0228                 ax1.set_ylim(0-(Y_error[0][4]*10),0+(Y_error[0][4]*10))
0229                 center = 0
0230             for each_bin in range(len(Y_error[0][0])):
0231                 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')
0232
0233             ax1.set_title('%s %s with error bars %s  %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
0234             plt.savefig(os.path.join(r_path, '%s_%s_error_%s.png' %  (title_list[i],title,config)))
0235         else:
0236             ax1.set_title('%s %s  %s  %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
0237             plt.savefig(os.path.join(r_path, '%s_%s_%s.png' %  (title_list[i],title,config)))
0238         plt.close()
0239
0240
0241 ###################################################################################################
0242      #Ratio vs momentum
0243 ###################################################################################################
0244
0245     if i > 0: #for each variable theta, phi, and eta
0246         for j in range(len(M_list)): #Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
0247             X = X_list[j]
0248             Y = Y_list[j]
0249             M_mc = M_list[j]
0250             boolean_M = np.ones_like(M_mc) == 1
0251             X_s = np.array(ak.flatten(X[boolean_M]))
0252             Y_s = np.array(ak.flatten(Y[boolean_M]))
0253             M_s = np.array(ak.flatten(M_mc))
0254             ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
0255             X_plot[j] = M_s
0256             Y_plot[j] = ratio
0257
0258         title ='difference vs momentum'
0259         particle_nlist = ['All','Pions','Protons','Electrons','Neutrons','Photons']
0260         for iterate in [0,1]:
0261             fig = plt.figure()
0262             gs = fig.add_gridspec(3, 2, wspace=0)
0263             (ax1, ax2), (ax3, ax4),(ax5, ax6) = gs.subplots(sharex=True, sharey=True)
0264             ax1.scatter(X_plot[0], Y_plot[0], s = ssize)
0265             ax2.scatter(X_plot[1], Y_plot[1], s = ssize)
0266             ax3.scatter(X_plot[2], Y_plot[2], s = ssize)
0267             ax4.scatter(X_plot[3], Y_plot[3], s = ssize)
0268             ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
0269             ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
0270             ax1.set_xscale('log')
0271             ax1.set_ylabel('rc-mc')
0272             ax3.set_ylabel('rc-mc')
0273             ax5.set_ylabel('rc-mc')
0274             ax5.set_xlabel('Momentum mc')
0275             ax6.set_xlabel('Momentum mc')
0276             ax2.set_title('Pions')
0277             ax3.set_title('Protons')
0278             ax4.set_title('Electrons')
0279             ax5.set_title('Neutrons')
0280             ax6.set_title('Photons')
0281             fig.set_figwidth(20)
0282             fig.set_figheight(10)
0283             if iterate == 0:
0284                 for j in range(len(X_list)):#Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
0285                     Y_error[j] = error_bars(X_plot[j], Y_plot[j], momentum_range)
0286                 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)
0287                 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)
0288                 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)
0289                 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)
0290                 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)
0291                 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)
0292                 ax1.set_ylim(0-(Y_error[0][4]*10),0+(Y_error[0][4]*10))
0293                 for each_bin in range(len(Y_error[0][0])):
0294                     ax1.text(x=Y_error[0][0][each_bin],y=0 + Y_error[0][4]*7,
0295                         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')
0296                 ax1.set_title('%s %s with error bars %s  %s events\n DETECTOR_CONFIG: %s'%(title_list[i],title,config,Nevents,Dconfig))
0297                 plt.savefig(os.path.join(r_path, '%s_difference_vs_momentum_error_%s.png' %  (title_list[i],config)))
0298             else:
0299                 ax1.set_title('%s Difference Vs Momentum  %s  %s events\n DETECTOR_CONFIG: %s'%(title_list[i],config,Nevents,Dconfig))
0300                 plt.savefig(os.path.join(r_path, '%s_difference_vs_momentum_%s.png' %  (title_list[i],config)))
0301             plt.close()
0302
0303
0304 ###################################################################################################
0305      #Correlation
0306 ###################################################################################################
0307
0308     #Repeat the following steps for each variable (momentum,theta,phi,eta)
0309     F_boolean = same_length_lists(MCparts,RCparts)
0310     X_s = np.array(ak.flatten(MCparts[F_boolean])) #Filtered lists
0311     Y_s = np.array(ak.flatten(RCparts[F_boolean]))
0312
0313     #Histogram
0314     if i == 0 and particle in particle_dict.keys(): #Momentum in Single events
0315         h, xedges, yedges = np.histogram2d(x=X_s,y= Y_s,  bins = 11)
0316     else:
0317         h, xedges, yedges = np.histogram2d(x=X_s,y= Y_s,  bins = 11, range = [x_range,x_range])
0318
0319     col_sum = ak.sum(h,axis=-1) #number of events in each (verticle) column
0320     norm_h = [] #norm_h is the normalized matrix
0321     norm_h_text = []
0322     for j in range(len(col_sum)):
0323         if col_sum[j] != 0:
0324             norm_c = h[j]/col_sum[j] #normalized column = column values divide by sum of the column
0325         else:
0326             norm_c = h[j]
0327         norm_h.append(norm_c)
0328         norm_c_text = [ '%.3f' % elem for elem in norm_c ] #display value to 3 dp
0329         norm_h_text.append(norm_c_text)
0330
0331     fig, axs = plt.subplots(1, 2, figsize=(20, 10))
0332     if i == 0 and particle in particle_dict.keys(): #Momentum in Single events
0333         axs[0].hist2d(x=X_s,y=Y_s, bins = 11)
0334     else:
0335         axs[0].hist2d(x=X_s,y=Y_s, bins = 11,range = [x_range,x_range])
0336     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])
0337     axs[0].set_title('%s Histogram'%(title_list[i]))
0338     axs[0].set_ylabel('%s_rc'%(title_list[i]))
0339     axs[1].set_ylabel('%s_rc'%(title_list[i]))
0340     axs[0].set_xlabel('%s_mc'%(title_list[i]))
0341     axs[1].set_xlabel('%s_mc'%(title_list[i]))
0342     axs[1].set_title('%s Correlation'%(title_list[i]))
0343     fig.suptitle('%s  %s events\n DETECTOR_CONFIG: %s'%(config,Nevents,Dconfig))
0344     plt.savefig(os.path.join(r_path, '%s_correlation_%s.png' %  (title_list[i],config)))
0345     plt.close()
0346
0347
0348 ###################################################################################################
0349      #Phi vs Theta plots
0350 ###################################################################################################
0351
0352 def particle_plots(boolean_particle):
0353     #filtered lists w.r.t the particle
0354     theta_mc_fil = ak.Array(theta_mc[simID][booll])[boolean_particle]
0355     theta_rc_fil = ak.Array(theta_rc[recID][booll])[boolean_particle]
0356     phi_mc_fil = ak.Array(phi_mc[simID][booll])[boolean_particle]
0357     phi_rc_fil = ak.Array(phi_rc[recID][booll])[boolean_particle]
0358
0359     F_boolean = same_length_lists(theta_mc_fil, theta_rc_fil)
0360     #filtered lists w.r.t length
0361     theta_mc_F = np.array(ak.flatten(theta_mc_fil[F_boolean]))
0362     theta_rc_F = np.array(ak.flatten(theta_rc_fil[F_boolean]))
0363     phi_mc_F = np.array(ak.flatten(phi_mc_fil[F_boolean]))
0364     phi_rc_F = np.array(ak.flatten(phi_rc_fil[F_boolean]))
0365     ratio = np.array((ak.Array(theta_rc_F)-(ak.Array(theta_mc_F))))
0366     x_range = [0,np.pi]
0367     Y_error = [error_bars(theta_mc_F, ratio, x_range),error_bars(theta_rc_F, ratio, x_range)]
0368     fig = plt.figure()
0369     gs = fig.add_gridspec(3, 2, wspace=0, hspace = 0.3)
0370     (ax1, ax2), (ax3, ax4), (ax5, ax6) = gs.subplots(sharex=True, sharey='row')
0371     ax1.scatter(-theta_mc_F, ratio, s = ssize)
0372     ax2.scatter(-theta_rc_F, ratio, s = ssize)
0373     ax3.scatter(-theta_mc_F, ratio, s = ssize)
0374     ax4.scatter(-theta_rc_F, ratio, s = ssize)
0375     ax5.scatter(-theta_mc_F, phi_mc_F, s = ssize)
0376     ax6.scatter(-theta_rc_F, phi_rc_F, s = ssize)
0377     ax1.set_ylabel('Theta rc-mc')
0378     ax2.set_ylabel('Theta rc-mc')
0379     ax3.set_ylabel('Theta rc-mc')
0380     ax4.set_ylabel('Theta rc-mc')
0381     ax5.set_ylabel('Phi mc')
0382     ax6.set_ylabel('Phi rc')
0383     ax1.set_xlabel('- Theta mc')
0384     ax2.set_xlabel('- Theta rc')
0385     ax3.set_xlabel('- Theta mc')
0386     ax4.set_xlabel('- Theta rc')
0387     ax5.set_xlabel('- Theta mc')
0388     ax6.set_xlabel('- Theta rc')
0389     ax1.set_title('Zoom-in')
0390     ax2.set_title('Zoom-in')
0391     ax3.set_title('Zoom-out')
0392     ax4.set_title('Zoom-out')
0393     ax5.set_title('Phi vs Theta mc')
0394     ax6.set_title('Phi vs Theta rc')
0395     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)
0396     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)
0397     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)
0398     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)
0399     y_limits = ax3.get_ylim()
0400     for each_bin in range(len(Y_error[0][0])):
0401         if not np.isnan(Y_error[0][1][each_bin]):
0402             ax3.text(x=-Y_error[0][0][each_bin],y=y_limits[1],
0403                     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')
0404         if not np.isnan(Y_error[1][1][each_bin]):
0405             ax4.text(x=-Y_error[1][0][each_bin],y=y_limits[1],
0406                     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')
0407     if not np.isnan(Y_error[0][4]):
0408         ax1.set_ylim(0-(Y_error[1][4]*10),0+(Y_error[1][4]*10))
0409         ax2.set_ylim(0-(Y_error[1][4]*10),0+(Y_error[1][4]*10))
0410     fig.set_figwidth(20)
0411     fig.set_figheight(10)
0412 title ='difference'
0413 if particle in particle_dict.keys():
0414     boolean_particle = particle_dict[particle][0]
0415     particle_name = particle_dict[particle][1]
0416     particle_plots(boolean_particle)
0417
0418     plt.suptitle('%s in %s  %s events\n DETECTOR_CONFIG: %s'%(particle_name,config,Nevents,Dconfig))
0419     plt.savefig(os.path.join(r_path, '%s_%s.png' %  (particle_name,config)))
0420 else:
0421     for i in [[boolean_photon,'Photons'],[boolean_electron,'Electrons'],[boolean_pion,'Pions']]:
0422         boolean_particle = i[0]
0423         particle_name = i[1]
0424         particle_plots(boolean_particle)
0425
0426         plt.suptitle('%s in %s  %s events\n DETECTOR_CONFIG: %s'%(particle_name,config,Nevents,Dconfig))
0427         plt.savefig(os.path.join(r_path, '%s_%s.png' %  (particle_name,config)))