diff --git a/benchmarks/dis/analysis/truth_reconstruction.py b/benchmarks/dis/analysis/truth_reconstruction.py
new file mode 100644
index 0000000000000000000000000000000000000000..9118ec8a35121356f8e1acd84160156a74dfcc32
--- /dev/null
+++ b/benchmarks/dis/analysis/truth_reconstruction.py
@@ -0,0 +1,750 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[1]:
+
+
+import os
+import numpy as np
+import uproot as ur
+import awkward as ak
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import mplhep 
+import argparse
+
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument('--rec_file', type=str, help='Reconstructed track file.')
+parser.add_argument('--ebeam', type=float, help='Electron beam energy.')
+parser.add_argument('--pbeam', type=float, help='Proton (or ion) beam energy.')
+parser.add_argument('--minq2', type=float, help='Minimum four-momentum transfer squared Q2.')
+parser.add_argument('--nevents', type=float, help='number of events to process.')
+parser.add_argument('-o', dest='outdir', default='results/dis/', help='Output directory.')
+args = parser.parse_args()
+kwargs = vars(args)
+
+rec_file = args.rec_file
+minq2 = int(args.minq2)
+k = int(args.ebeam)
+p = int(args.pbeam)
+Nevents = int(args.nevents)
+
+
+# In[16]:
+
+
+for array in ur.iterate(rec_file + ':events',['MCParticles/MCParticles.generatorStatus',
+                                          'MCParticles/MCParticles.mass',
+                                          'MCParticles/MCParticles.PDG',
+                                          'MCParticles/MCParticles.momentum.x',
+                                          'MCParticles/MCParticles.momentum.y',
+                                          'MCParticles/MCParticles.momentum.z',
+                                          'ReconstructedParticles/ReconstructedParticles.energy',
+                                          'ReconstructedParticles/ReconstructedParticles.PDG',
+                                          'ReconstructedParticles/ReconstructedParticles.momentum.x',
+                                          'ReconstructedParticles/ReconstructedParticles.momentum.y',
+                                          'ReconstructedParticles/ReconstructedParticles.momentum.z',
+                                          'ReconstructedParticlesAssoc/ReconstructedParticlesAssoc.simID',
+                                          'ReconstructedParticlesAssoc/ReconstructedParticlesAssoc.recID',
+                                         ],step_size=Nevents):
+    # generatorStatus = array['MCParticles/MCParticles.generatorStatus']
+    mass = array['MCParticles/MCParticles.mass']
+    PDG_mc = array['MCParticles/MCParticles.PDG']
+    px_mc = array['MCParticles/MCParticles.momentum.x']
+    py_mc = array['MCParticles/MCParticles.momentum.y']
+    pz_mc = array['MCParticles/MCParticles.momentum.z']
+    PDG_rc = array['ReconstructedParticles/ReconstructedParticles.PDG']
+    px_rc = array['ReconstructedParticles/ReconstructedParticles.momentum.x']
+    py_rc = array['ReconstructedParticles/ReconstructedParticles.momentum.y']
+    pz_rc = array['ReconstructedParticles/ReconstructedParticles.momentum.z']
+    simID = array['ReconstructedParticlesAssoc/ReconstructedParticlesAssoc.simID']
+    recID = array['ReconstructedParticlesAssoc/ReconstructedParticlesAssoc.recID']
+
+momentum_mc = np.sqrt(((px_mc**2)+(py_mc**2)+(pz_mc**2)))
+momentum_rc = np.sqrt(((px_rc**2)+(py_rc**2)+(pz_rc**2)))
+
+theta_mc = np.arctan2(np.sqrt(px_mc**2+py_mc**2), pz_mc)
+phi_mc = np.arctan2(py_mc, px_mc)
+theta_rc = np.arctan2(np.sqrt(px_rc**2+py_rc**2), pz_rc)
+phi_rc = np.arctan2(py_rc, px_rc)
+
+
+# In[ ]:
+
+
+booll = (PDG_mc[simID])==(PDG_rc[recID])
+MC_list = [ak.Array(momentum_mc[simID][booll]),
+           ak.Array(theta_mc[simID][booll]),
+           ak.Array(phi_mc[simID][booll]),
+           -np.log(np.tan((ak.Array(theta_mc[simID][booll]))/2))]
+RC_list = [ak.Array(momentum_rc[recID][booll]),
+           ak.Array(theta_rc[recID][booll]),
+           ak.Array(phi_rc[recID][booll]),
+           -np.log(np.tan((ak.Array(theta_rc[recID][booll]))/2))]
+title_list = ['Momentum','Theta','Phi','Eta']
+for i in range(len(MC_list)):
+    X = MC_list[i]
+    Y = RC_list[i]
+
+    X_len = ak.count(X,axis=None)
+    Y_len = ak.count(Y,axis=None)
+
+    if X_len > Y_len: 
+        Y_boolean = np.ones_like(Y) == 1
+        X_F = X[Y_boolean] 
+        Y_F = Y[Y_boolean] 
+    else: 
+        X_boolean = np.ones_like(X) == 1
+        X_F = X[X_boolean] 
+        Y_F = Y[X_boolean] 
+
+    X_s = np.array(ak.flatten(X_F)) 
+    Y_s = np.array(ak.flatten(Y_F)) 
+    
+    if i == 3:
+        h, xedges, yedges, image = plt.hist2d(x=X_s,y= Y_s,  bins = 11,range = [[-5,5],[-5,5]]) 
+        plt.close()
+    else:
+        h, xedges, yedges, image = plt.hist2d(x=X_s,y= Y_s, bins = 10) 
+        plt.close()
+    
+    col_sum = ak.sum(h,axis=-1) #number of events in each (verticle) column 
+    norm_h = [] #norm_h is the normalized matrix
+    norm_h_text = [] #display labels matrix
+    for j in range(len(col_sum)):
+        if col_sum[j] != 0:
+            norm_c = h[j]/col_sum[j] #normalized column = column values divide by sum of the column
+        else:
+            norm_c = h[j]
+        norm_h.append(norm_c)
+        norm_c_text = [ '%.3f' % elem for elem in norm_c ] #display value to 3 dp
+        norm_h_text.append(norm_c_text)
+    
+    fig, axs = plt.subplots(1, 2, figsize=(20, 10))
+    if i == 3:
+        axs[0].hist2d(x=X_s,y=Y_s, bins = 11,range = [[-5,5],[-5,5]]) 
+    else:
+        axs[0].hist2d(x=X_s,y=Y_s, bins = 10) 
+    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])
+    fig.suptitle('(%g x %g)GeV  %gGeV minQ2  %s events'%(k,p,minq2,Nevents))
+    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[1].set_title('%s Correlation'%(title_list[i]))
+    plt.show()
+
+
+# # Momentum
+
+# In[30]:
+
+
+booll = (PDG_mc[simID])==(PDG_rc[recID])
+X1 = ak.Array(momentum_mc[simID][booll])
+Y1 = ak.Array(momentum_rc[recID][booll])
+
+boolean_pion = np.logical_or(ak.Array(PDG_mc[simID][booll])==-211, ak.Array(PDG_mc[simID][booll])==+211)
+boolean_proton = np.logical_or(ak.Array(PDG_mc[simID][booll])==-2212, ak.Array(PDG_mc[simID][booll])==+2212)
+boolean_electron = ak.Array(PDG_mc[simID][booll])==11
+boolean_neutron = ak.Array(PDG_mc[simID][booll])==2112
+boolean_photon = ak.Array(PDG_mc[simID][booll])==22
+
+X_list = [ak.Array(X1),
+          ak.Array(X1[boolean_pion]),
+          ak.Array(X1[boolean_proton]),
+          ak.Array(X1[boolean_electron]),
+          ak.Array(X1[boolean_neutron]),
+          ak.Array(X1[boolean_photon])]
+Y_list = [ak.Array(Y1),
+          ak.Array(Y1[boolean_pion]),
+          ak.Array(Y1[boolean_proton]),
+          ak.Array(Y1[boolean_electron]),
+          ak.Array(Y1[boolean_neutron]),
+          ak.Array(Y1[boolean_photon])]
+
+X_plot = [0,0,0,0,0,0]
+Y_plot = [0,0,0,0,0,0]
+
+for i in range(len(X_list)):
+    X = X_list[i]
+    Y = Y_list[i]
+    X_len = ak.count(X,axis=None)
+    Y_len = ak.count(Y,axis=None)
+
+    if X_len > Y_len: 
+        Y_boolean = np.ones_like(Y) == 1
+        X_F = X[Y_boolean] 
+        Y_F = Y[Y_boolean] 
+    else: 
+        X_boolean = np.ones_like(X) == 1
+        X_F = X[X_boolean] 
+        Y_F = Y[X_boolean] 
+
+    X_s = np.array(ak.flatten(X_F)) 
+    Y_s = np.array(ak.flatten(Y_F))
+
+    ratio = np.array((ak.Array(Y_s)/ak.Array(X_s)))
+    X_plot[i] = X_s
+    Y_plot[i] = 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('')
+ax1.scatter(-X_plot[0], Y_plot[0], s = 0.01)
+ax2.scatter(-X_plot[1], Y_plot[1], s = 0.01)
+ax3.scatter(-X_plot[2], Y_plot[2], s = 0.01)
+ax4.scatter(-X_plot[3], Y_plot[3], s = 0.01)
+ax5.scatter(-X_plot[4], Y_plot[4], s = 0.01)
+ax6.scatter(-X_plot[5], Y_plot[5], s = 0.01)
+ax1.set_yscale('log')
+ax1.set_xscale('log')
+ax2.set_yscale('log')
+ax2.set_xscale('log')
+ax3.set_yscale('log')
+ax3.set_xscale('log')
+ax4.set_yscale('log')
+ax4.set_xscale('log')
+ax5.set_yscale('log')
+ax5.set_xscale('log')
+ax6.set_yscale('log')
+ax6.set_xscale('log')
+ax1.set_title('Momentum Ratio   %s - %s   %s events'%(minE,maxE,tevents))
+ax2.set_title('Pions')
+ax3.set_title('Protons')
+ax4.set_title('Electrons')
+ax5.set_title('Neutrons')
+ax6.set_title('Photons')
+ax1.set_ylabel('rc/mc')
+ax3.set_ylabel('rc/mc')
+ax5.set_ylabel('rc/mc')
+ax5.set_xlabel('Momentum')
+ax6.set_xlabel('Momentum')
+fig.set_figwidth(16)
+fig.set_figheight(10)
+plt.show()
+
+
+# # Theta
+
+# In[21]:
+
+
+n = 'Theta'
+n_s = 'theta'
+booll = (PDG_mc[simID])==(PDG_rc[recID])
+X1 = ak.Array(theta_mc[simID][booll])
+Y1 = ak.Array(theta_rc[recID][booll])
+
+boolean_pion = np.logical_or(ak.Array(PDG_mc[simID][booll])==-211, ak.Array(PDG_mc[simID][booll])==+211)
+boolean_proton = np.logical_or(ak.Array(PDG_mc[simID][booll])==-2212, ak.Array(PDG_mc[simID][booll])==+2212)
+boolean_electron = ak.Array(PDG_mc[simID][booll])==11
+boolean_neutron = ak.Array(PDG_mc[simID][booll])==2112
+boolean_photon = ak.Array(PDG_mc[simID][booll])==22
+
+X_list = [ak.Array(X1),
+          ak.Array(X1[boolean_pion]),
+          ak.Array(X1[boolean_proton]),
+          ak.Array(X1[boolean_electron]),
+          ak.Array(X1[boolean_neutron]),
+          ak.Array(X1[boolean_photon])]
+Y_list = [ak.Array(Y1),
+          ak.Array(Y1[boolean_pion]),
+          ak.Array(Y1[boolean_proton]),
+          ak.Array(Y1[boolean_electron]),
+          ak.Array(Y1[boolean_neutron]),
+          ak.Array(Y1[boolean_photon])]
+
+X_plot = [0,0,0,0,0,0]
+Y_plot = [0,0,0,0,0,0]
+
+for i in range(len(X_list)):
+    X = X_list[i]
+    Y = Y_list[i]
+    X_len = ak.count(X,axis=None)
+    Y_len = ak.count(Y,axis=None)
+
+    if X_len > Y_len: 
+        Y_boolean = np.ones_like(Y) == 1
+        X_F = X[Y_boolean] 
+        Y_F = Y[Y_boolean] 
+    else: 
+        X_boolean = np.ones_like(X) == 1
+        X_F = X[X_boolean] 
+        Y_F = Y[X_boolean] 
+
+    X_s = np.array(ak.flatten(X_F)) 
+    Y_s = np.array(ak.flatten(Y_F))
+
+    ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
+    X_plot[i] = X_s
+    Y_plot[i] = 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('')
+ax1.scatter(-X_plot[0], Y_plot[0], s = 0.01)
+ax2.scatter(-X_plot[1], Y_plot[1], s = 0.01)
+ax3.scatter(-X_plot[2], Y_plot[2], s = 0.01)
+ax4.scatter(-X_plot[3], Y_plot[3], s = 0.01)
+ax5.scatter(-X_plot[4], Y_plot[4], s = 0.01)
+ax6.scatter(-X_plot[5], Y_plot[5], s = 0.01)
+ax1.set_title('%s Difference   %s - %s   %s events'%(n,minE,maxE,tevents))
+ax2.set_title('Pions')
+ax3.set_title('Protons')
+ax4.set_title('Electrons')
+ax5.set_title('Neutrons')
+ax6.set_title('Photons')
+ax1.set_ylabel('rc-mc')
+ax3.set_ylabel('rc-mc')
+ax5.set_ylabel('rc-mc')
+ax5.set_xlabel('%s'%n)
+ax6.set_xlabel('%s'%n)
+fig.set_figwidth(16)
+fig.set_figheight(10)
+plt.show()
+
+
+# In[22]:
+
+
+booll = (PDG_mc[simID])==(PDG_rc[recID])
+X1 = ak.Array(theta_mc[simID][booll])
+Y1 = ak.Array(theta_rc[recID][booll])
+# X1 = np.log(np.tan(X1/2))
+# Y1 = np.log(np.tan(Y1/2))
+
+X2 = ak.Array(phi_mc[simID][booll])
+Y2 = ak.Array(phi_rc[recID][booll])
+
+boolean_photon = ak.Array(PDG_mc[simID][booll])==22
+
+X1 = ak.Array(X1[boolean_photon])
+Y1 = ak.Array(Y1[boolean_photon])
+X2 = ak.Array(X2[boolean_photon])
+Y2 = ak.Array(Y2[boolean_photon])
+X_len = ak.count(X1,axis=None)
+Y_len = ak.count(Y1,axis=None)
+
+if X_len > Y_len: 
+    Y_boolean = np.ones_like(Y1) == 1
+    X_F = X1[Y_boolean] 
+    Y_F = Y1[Y_boolean]
+    x_f = X2[Y_boolean]
+    y_f = Y2[Y_boolean]
+else: 
+    X_boolean = np.ones_like(X1) == 1
+    X_F = X1[X_boolean] 
+    Y_F = Y1[X_boolean] 
+    x_f = X2[X_boolean]
+    y_f = Y2[X_boolean]
+
+X_s = np.array(ak.flatten(X_F)) 
+Y_s = np.array(ak.flatten(Y_F))
+x_s = np.array(ak.flatten(x_f))
+y_s = np.array(ak.flatten(y_f))
+
+ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
+boolean_range = np.absolute(ratio) > 0
+X_plot = X_s[boolean_range]
+Y_plot = ratio[boolean_range]
+x_plot = x_s[boolean_range]
+y_plot = y_s[boolean_range]
+
+fig = plt.figure()
+plt.scatter(-X_plot, Y_plot, s = 0.01)
+plt.xlabel('theta mc')
+plt.ylabel('rc-mc')
+plt.title('%s - %s   %s events\nPhotons Theta Difference vs Theta MC'%(minE,maxE,tevents))
+fig.set_figwidth(16)
+fig.set_figheight(6)
+plt.show()
+
+X_plot = Y_s[boolean_range]
+fig = plt.figure()
+plt.scatter(-X_plot, Y_plot, s = 0.01)
+plt.xlabel('theta rc')
+plt.ylabel('rc-mc')
+plt.title('%s - %s   %s events\nPhotons Theta Difference vs Theta RC'%(minE,maxE,tevents))
+fig.set_figwidth(16)
+fig.set_figheight(6)
+plt.show()
+
+X_plot = X_s[boolean_range]
+fig = plt.figure()
+plt.scatter(-X_plot, x_plot, s = 0.01)
+plt.xlabel('theta mc')
+plt.ylabel('phi mc')
+plt.title('')
+fig.set_figwidth(16)
+fig.set_figheight(6)
+plt.show()
+
+X_plot = Y_s[boolean_range]
+fig = plt.figure()
+plt.scatter(-X_plot, y_plot, s = 0.01)
+plt.xlabel('theta rc')
+plt.ylabel('phi rc')
+plt.title('%s - %s   %s events\nPhotons Phi RC vs Theta RC'%(minE,maxE,tevents))
+# plt.title('')
+fig.set_figwidth(16)
+fig.set_figheight(6)
+plt.show()
+
+
+X_plot = Y_s[boolean_range]
+fig = plt.figure()
+plt.scatter(y_plot,-X_plot, s = 0.01)
+            # y_plot, s = 0.01)
+plt.xlabel('phi rc')
+plt.ylabel('theta rc')
+plt.title('%s - %s   %s events\nPhotons Theta RC vs Phi RC'%(minE,maxE,tevents))
+# plt.title('')
+fig.set_figwidth(16)
+fig.set_figheight(6)
+plt.show()
+
+
+# In[23]:
+
+
+M_list = [ak.Array(momentum_mc[simID][booll]),
+          ak.Array(momentum_mc[simID][booll][boolean_pion]),
+          ak.Array(momentum_mc[simID][booll][boolean_proton]),
+          ak.Array(momentum_mc[simID][booll][boolean_electron]),
+          ak.Array(momentum_mc[simID][booll][boolean_neutron]),
+          ak.Array(momentum_mc[simID][booll][boolean_photon])]
+          
+X_plot = [0,0,0,0,0,0]
+Y_plot = [0,0,0,0,0,0]
+
+for i in range(len(X_list)):
+    X = X_list[i]
+    Y = Y_list[i]
+    M_mc = M_list[i]
+    boolean_M = np.ones_like(M_mc) == 1
+    X_F = X[boolean_M] 
+    Y_F = Y[boolean_M] 
+    X_s = np.array(ak.flatten(X_F)) 
+    Y_s = np.array(ak.flatten(Y_F)) 
+    M_s = np.array(ak.flatten(M_mc))
+    ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
+    X_plot[i] = M_s
+    Y_plot[i] = 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('')
+ax1.scatter(X_plot[0], Y_plot[0], s = 0.01)
+ax2.scatter(X_plot[1], Y_plot[1], s = 0.01)
+ax3.scatter(X_plot[2], Y_plot[2], s = 0.01)
+ax4.scatter(X_plot[3], Y_plot[3], s = 0.01)
+ax5.scatter(X_plot[4], Y_plot[4], s = 0.01)
+ax6.scatter(X_plot[5], Y_plot[5], s = 0.01)
+ax1.set_title('%s Difference Vs Momentum   %s - %s   %s events'%(n,minE,maxE,tevents))
+ax2.set_title('Pions')
+ax3.set_title('Protons')
+ax4.set_title('Electrons')
+ax5.set_title('Neutrons')
+ax6.set_title('Photons')
+ax1.set_ylabel('rc-mc')
+ax3.set_ylabel('rc-mc')
+ax5.set_ylabel('rc-mc')
+ax5.set_xlabel('Momentum')
+ax6.set_xlabel('Momentum')
+fig.set_figwidth(16)
+fig.set_figheight(10)
+ax1.set_xscale('log')
+ax2.set_xscale('log')
+ax3.set_xscale('log')
+ax4.set_xscale('log')
+ax5.set_xscale('log')
+ax6.set_xscale('log')
+plt.show()
+
+
+# # Phi
+
+# In[32]:
+
+
+n = 'Phi'
+n_s = 'phi'
+booll = (PDG_mc[simID])==(PDG_rc[recID])
+X1 = ak.Array(phi_mc[simID][booll])
+Y1 = ak.Array(phi_rc[recID][booll])
+
+boolean_pion = np.logical_or(ak.Array(PDG_mc[simID][booll])==-211, ak.Array(PDG_mc[simID][booll])==+211)
+boolean_proton = np.logical_or(ak.Array(PDG_mc[simID][booll])==-2212, ak.Array(PDG_mc[simID][booll])==+2212)
+boolean_electron = ak.Array(PDG_mc[simID][booll])==11
+boolean_neutron = ak.Array(PDG_mc[simID][booll])==2112
+boolean_photon = ak.Array(PDG_mc[simID][booll])==22
+
+X_list = [ak.Array(X1),
+          ak.Array(X1[boolean_pion]),
+          ak.Array(X1[boolean_proton]),
+          ak.Array(X1[boolean_electron]),
+          ak.Array(X1[boolean_neutron]),
+          ak.Array(X1[boolean_photon])]
+Y_list = [ak.Array(Y1),
+          ak.Array(Y1[boolean_pion]),
+          ak.Array(Y1[boolean_proton]),
+          ak.Array(Y1[boolean_electron]),
+          ak.Array(Y1[boolean_neutron]),
+          ak.Array(Y1[boolean_photon])]
+X_plot = [0,0,0,0,0,0]
+Y_plot = [0,0,0,0,0,0]
+
+for i in range(len(X_list)):
+    X = X_list[i]
+    Y = Y_list[i]
+    X_len = ak.count(X,axis=None)
+    Y_len = ak.count(Y,axis=None)
+
+    if X_len > Y_len: 
+        Y_boolean = np.ones_like(Y) == 1
+        X_F = X[Y_boolean] 
+        Y_F = Y[Y_boolean] 
+    else: 
+        X_boolean = np.ones_like(X) == 1
+        X_F = X[X_boolean] 
+        Y_F = Y[X_boolean] 
+
+    X_s = np.array(ak.flatten(X_F)) 
+    Y_s = np.array(ak.flatten(Y_F))
+
+    ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
+    X_plot[i] = X_s
+    Y_plot[i] = 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('')
+ax1.scatter(X_plot[0], Y_plot[0], s = 0.01)
+ax2.scatter(X_plot[1], Y_plot[1], s = 0.01)
+ax3.scatter(X_plot[2], Y_plot[2], s = 0.01)
+ax4.scatter(X_plot[3], Y_plot[3], s = 0.01)
+ax5.scatter(X_plot[4], Y_plot[4], s = 0.01)
+ax6.scatter(X_plot[5], Y_plot[5], s = 0.01)
+ax1.set_title('%s Difference   %s - %s   %s events'%(n,minE,maxE,tevents))
+ax2.set_title('Pions')
+ax3.set_title('Protons')
+ax4.set_title('Electrons')
+ax5.set_title('Neutrons')
+ax6.set_title('Photons')
+ax1.set_ylabel('rc-mc')
+ax3.set_ylabel('rc-mc')
+ax5.set_ylabel('rc-mc')
+ax5.set_xlabel('%s'%n)
+ax6.set_xlabel('%s'%n)
+fig.set_figwidth(16)
+fig.set_figheight(10)
+plt.show()
+
+
+# In[26]:
+
+
+M_list = [ak.Array(momentum_mc[simID][booll]),
+          ak.Array(momentum_mc[simID][booll][boolean_pion]),
+          ak.Array(momentum_mc[simID][booll][boolean_proton]),
+          ak.Array(momentum_mc[simID][booll][boolean_electron]),
+          ak.Array(momentum_mc[simID][booll][boolean_neutron]),
+          ak.Array(momentum_mc[simID][booll][boolean_photon])]
+          
+X_plot = [0,0,0,0,0,0]
+Y_plot = [0,0,0,0,0,0]
+
+for i in range(len(X_list)):
+    X = X_list[i]
+    Y = Y_list[i]
+    M_mc = M_list[i]
+    boolean_M = np.ones_like(M_mc) == 1
+    X_F = X[boolean_M] 
+    Y_F = Y[boolean_M] 
+    X_s = np.array(ak.flatten(X_F)) 
+    Y_s = np.array(ak.flatten(Y_F)) 
+    M_s = np.array(ak.flatten(M_mc))
+    ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
+    X_plot[i] = M_s
+    Y_plot[i] = 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('')
+ax1.scatter(X_plot[0], Y_plot[0], s = 0.01)
+ax2.scatter(X_plot[1], Y_plot[1], s = 0.01)
+ax3.scatter(X_plot[2], Y_plot[2], s = 0.01)
+ax4.scatter(X_plot[3], Y_plot[3], s = 0.01)
+ax5.scatter(X_plot[4], Y_plot[4], s = 0.01)
+ax6.scatter(X_plot[5], Y_plot[5], s = 0.01)
+ax1.set_title('%s Difference Vs Momentum   %s - %s   %s events'%(n,minE,maxE,tevents))
+ax2.set_title('Pions')
+ax3.set_title('Protons')
+ax4.set_title('Electrons')
+ax5.set_title('Neutrons')
+ax6.set_title('Photons')
+ax1.set_ylabel('rc-mc')
+ax3.set_ylabel('rc-mc')
+ax5.set_ylabel('rc-mc')
+ax5.set_xlabel('Momentum')
+ax6.set_xlabel('Momentum')
+fig.set_figwidth(16)
+fig.set_figheight(10)
+ax1.set_xscale('log')
+ax2.set_xscale('log')
+ax3.set_xscale('log')
+ax4.set_xscale('log')
+ax5.set_xscale('log')
+ax6.set_xscale('log')
+plt.show()
+
+
+# # Eta
+
+# In[28]:
+
+
+n = 'Eta'
+n_s = 'eta'
+booll = (PDG_mc[simID])==(PDG_rc[recID])
+X1 = ak.Array(theta_mc[simID][booll])
+Y1 = ak.Array(theta_rc[recID][booll])
+X1 = -np.log(np.tan(X1/2))
+Y1 = -np.log(np.tan(Y1/2))
+
+boolean_pion = np.logical_or(ak.Array(PDG_mc[simID][booll])==-211, ak.Array(PDG_mc[simID][booll])==+211)
+boolean_proton = np.logical_or(ak.Array(PDG_mc[simID][booll])==-2212, ak.Array(PDG_mc[simID][booll])==+2212)
+boolean_electron = ak.Array(PDG_mc[simID][booll])==11
+boolean_neutron = ak.Array(PDG_mc[simID][booll])==2112
+boolean_photon = ak.Array(PDG_mc[simID][booll])==22
+
+X_list = [ak.Array(X1),
+          ak.Array(X1[boolean_pion]),
+          ak.Array(X1[boolean_proton]),
+          ak.Array(X1[boolean_electron]),
+          ak.Array(X1[boolean_neutron]),
+          ak.Array(X1[boolean_photon])]
+Y_list = [ak.Array(Y1),
+          ak.Array(Y1[boolean_pion]),
+          ak.Array(Y1[boolean_proton]),
+          ak.Array(Y1[boolean_electron]),
+          ak.Array(Y1[boolean_neutron]),
+          ak.Array(Y1[boolean_photon])]
+X_plot = [0,0,0,0,0,0]
+Y_plot = [0,0,0,0,0,0]
+
+for i in range(len(X_list)):
+    X = X_list[i]
+    Y = Y_list[i]
+    X_len = ak.count(X,axis=None)
+    Y_len = ak.count(Y,axis=None)
+
+    if X_len > Y_len: 
+        Y_boolean = np.ones_like(Y) == 1
+        X_F = X[Y_boolean] 
+        Y_F = Y[Y_boolean] 
+    else: 
+        X_boolean = np.ones_like(X) == 1
+        X_F = X[X_boolean] 
+        Y_F = Y[X_boolean] 
+
+    X_s = np.array(ak.flatten(X_F)) 
+    Y_s = np.array(ak.flatten(Y_F))
+
+    ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
+    X_plot[i] = X_s
+    Y_plot[i] = 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('')
+ax1.scatter(X_plot[0], Y_plot[0], s = 0.01)
+ax2.scatter(X_plot[1], Y_plot[1], s = 0.01)
+ax3.scatter(X_plot[2], Y_plot[2], s = 0.01)
+ax4.scatter(X_plot[3], Y_plot[3], s = 0.01)
+ax5.scatter(X_plot[4], Y_plot[4], s = 0.01)
+ax6.scatter(X_plot[5], Y_plot[5], s = 0.01)
+ax1.set_title('%s Difference %s - %s   %s events'%(n,minE,maxE,tevents))
+ax2.set_title('Pions')
+ax3.set_title('Protons')
+ax4.set_title('Electrons')
+ax5.set_title('Neutrons')
+ax6.set_title('Photons')
+ax1.set_ylabel('rc-mc')
+ax3.set_ylabel('rc-mc')
+ax5.set_ylabel('rc-mc')
+ax5.set_xlabel('%s'%n)
+ax6.set_xlabel('%s'%n)
+fig.set_figwidth(16)
+fig.set_figheight(10)
+plt.show()
+
+
+# In[29]:
+
+
+M_list = [ak.Array(momentum_mc[simID][booll]),
+          ak.Array(momentum_mc[simID][booll][boolean_pion]),
+          ak.Array(momentum_mc[simID][booll][boolean_proton]),
+          ak.Array(momentum_mc[simID][booll][boolean_electron]),
+          ak.Array(momentum_mc[simID][booll][boolean_neutron]),
+          ak.Array(momentum_mc[simID][booll][boolean_photon])]
+          
+X_plot = [0,0,0,0,0,0]
+Y_plot = [0,0,0,0,0,0]
+
+for i in range(len(X_list)):
+    X = X_list[i]
+    Y = Y_list[i]
+    M_mc = M_list[i]
+    boolean_M = np.ones_like(M_mc) == 1
+    X_F = X[boolean_M] 
+    Y_F = Y[boolean_M] 
+    X_s = np.array(ak.flatten(X_F)) 
+    Y_s = np.array(ak.flatten(Y_F)) 
+    M_s = np.array(ak.flatten(M_mc))
+    ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
+    X_plot[i] = M_s
+    Y_plot[i] = 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('')
+ax1.scatter(X_plot[0], Y_plot[0], s = 0.01)
+ax2.scatter(X_plot[1], Y_plot[1], s = 0.01)
+ax3.scatter(X_plot[2], Y_plot[2], s = 0.01)
+ax4.scatter(X_plot[3], Y_plot[3], s = 0.01)
+ax5.scatter(X_plot[4], Y_plot[4], s = 0.01)
+ax6.scatter(X_plot[5], Y_plot[5], s = 0.01)
+ax1.set_title('%s Difference Vs Momentum %s - %s   %s events'%(n,minE,maxE,tevents))
+ax2.set_title('Pions')
+ax3.set_title('Protons')
+ax4.set_title('Electrons')
+ax5.set_title('Neutrons')
+ax6.set_title('Photons')
+ax1.set_ylabel('rc-mc')
+ax3.set_ylabel('rc-mc')
+ax5.set_ylabel('rc-mc')
+ax5.set_xlabel('Momentum')
+ax6.set_xlabel('Momentum')
+fig.set_figwidth(16)
+fig.set_figheight(10)
+ax1.set_xscale('log')
+ax2.set_xscale('log')
+ax3.set_xscale('log')
+ax4.set_xscale('log')
+ax5.set_xscale('log')
+ax6.set_xscale('log')
+plt.show()
+