diff --git a/benchmarks/dis/analysis/truth_reconstruction.py b/benchmarks/dis/analysis/truth_reconstruction.py
index 3089488c5a10fa24d385f8663488497eb7b83c5c..131185617f02303b5f2e61a39701b04b8a30741d 100644
--- a/benchmarks/dis/analysis/truth_reconstruction.py
+++ b/benchmarks/dis/analysis/truth_reconstruction.py
@@ -16,75 +16,86 @@ args = parser.parse_args()
kwargs = vars(args)
rec_file = args.rec_file
-config = args.config
Nevents = int(args.nevents)
-r_path = args.results_path + '/truth_reconstruction/'
+r_path = args.results_path + '/truth_reconstruction/' #Path for output figures and file.
+Dconfig = 'epic' + args.config.split('_epic')[1].strip() #Detector config
+config = args.config.split('_epic')[0].strip()
for array in ur.iterate(rec_file + ':events',['MCParticles/MCParticles.generatorStatus',
'MCParticles/MCParticles.PDG',
'MCParticles/MCParticles.momentum.x',
'MCParticles/MCParticles.momentum.y',
'MCParticles/MCParticles.momentum.z',
- 'ReconstructedParticles/ReconstructedParticles.PDG',
- 'ReconstructedParticles/ReconstructedParticles.momentum.x',
- 'ReconstructedParticles/ReconstructedParticles.momentum.y',
- 'ReconstructedParticles/ReconstructedParticles.momentum.z',
- 'ReconstructedParticlesAssoc/ReconstructedParticlesAssoc.simID',
- 'ReconstructedParticlesAssoc/ReconstructedParticlesAssoc.recID',],step_size=Nevents):
- PDG_mc = array['MCParticles/MCParticles.PDG']
+ 'ReconstructedChargedParticles/ReconstructedChargedParticles.PDG',
+ 'ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.x',
+ 'ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.y',
+ 'ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.z',
+ 'ReconstructedChargedParticlesAssoc/ReconstructedChargedParticlesAssoc.simID',
+ 'ReconstructedChargedParticlesAssoc/ReconstructedChargedParticlesAssoc.recID'],step_size=Nevents):
+ PDG_mc = array['MCParticles/MCParticles.PDG'] #Monte Carlo (MC) particle numbering scheme.
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']
-
+ PDG_rc = array['ReconstructedChargedParticles/ReconstructedChargedParticles.PDG']
+ px_rc = array['ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.x']
+ py_rc = array['ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.y']
+ pz_rc = array['ReconstructedChargedParticles/ReconstructedChargedParticles.momentum.z']
+ simID = array['ReconstructedChargedParticlesAssoc/ReconstructedChargedParticlesAssoc.simID']
+ recID = array['ReconstructedChargedParticlesAssoc/ReconstructedChargedParticlesAssoc.recID']
+ #SimID and recID contain the indices of the MCParticles and ReconstructedParticles entry for that event.
+
+### MCParticles Variables
momentum_mc = np.sqrt(((px_mc**2)+(py_mc**2)+(pz_mc**2)))
theta_mc = np.arctan2(np.sqrt(px_mc**2+py_mc**2), pz_mc)
phi_mc = np.arctan2(py_mc, px_mc)
-
+### ReconstructedParticles Variables
momentum_rc = np.sqrt(((px_rc**2)+(py_rc**2)+(pz_rc**2)))
theta_rc = np.arctan2(np.sqrt(px_rc**2+py_rc**2), pz_rc)
phi_rc = np.arctan2(py_rc, px_rc)
-booll = (PDG_mc[simID])==(PDG_rc[recID])
-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
-
-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))]
+booll = (PDG_mc[simID])==(PDG_rc[recID]) #boolean that allows events where the same particle is reconstructed
+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
+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
+boolean_electron = ak.Array(PDG_mc[simID][booll])==11 #boolean that allows events involving electrons
+boolean_neutron = ak.Array(PDG_mc[simID][booll])==2112 #boolean that allows events involving neutrons
+boolean_photon = ak.Array(PDG_mc[simID][booll])==22 #boolean that allows events involving photons
+
+### MCParticles variables list
+MC_list = [ak.Array(momentum_mc[simID][booll]), #Momentum
+ ak.Array(theta_mc[simID][booll]), #Theta
+ ak.Array(phi_mc[simID][booll]), #Phi
+ -np.log(np.tan((ak.Array(theta_mc[simID][booll]))/2))] #Eta
+### ReconstructedParticles variables list
+RC_list = [ak.Array(momentum_rc[recID][booll]), #Momentum
+ ak.Array(theta_rc[recID][booll]), #Theta
+ ak.Array(phi_rc[recID][booll]), #Phi
+ -np.log(np.tan((ak.Array(theta_rc[recID][booll]))/2))] #Eta
+title_list = ['Momentum','Theta','Phi','Eta']
+### MC Momentum for different particles list
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])]
-title_list = ['Momentum','Theta','Phi','Eta']
+#Marker Size in plots
if Nevents == 100:
ssize = 1
else:
ssize = 0.01
+#Particle type for Single events
particle = config.split('-')[0].strip()
particle_dict = {'e':[boolean_electron,'Electrons'],'pi':[boolean_pion,'Pions']}
-###################
-for i in range(len(MC_list)):
- X1 = MC_list[i]
- Y1 = RC_list[i]
+####################################################################################################
+ #Ratio
+####################################################################################################
+
+for i in range(len(MC_list)): #Repeat the following steps for each variable (momentum,theta,phi,eta)
+ X1 = MC_list[i] #MCParticles events to be plotted on x-axis
+ Y1 = RC_list[i] #ReconstructedParticles events
X_list = [ak.Array(X1),
ak.Array(X1[boolean_pion]),
ak.Array(X1[boolean_proton]),
@@ -97,11 +108,10 @@ for i in range(len(MC_list)):
ak.Array(Y1[boolean_electron]),
ak.Array(Y1[boolean_neutron]),
ak.Array(Y1[boolean_photon])]
-
X_plot = list(np.zeros(len(X_list)))
Y_plot = list(np.zeros(len(X_list)))
- for j in range(len(X_list)):
+ for j in range(len(X_list)): #Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
X = X_list[j]
Y = Y_list[j]
X_len = ak.count(X,axis=None)
@@ -110,11 +120,11 @@ for i in range(len(MC_list)):
F_boolean = np.ones_like(Y) == 1
else:
F_boolean = np.ones_like(X) == 1
- X_s = np.array(ak.flatten(X[F_boolean]))
+ X_s = np.array(ak.flatten(X[F_boolean])) #Filtered lists
Y_s = np.array(ak.flatten(Y[F_boolean]))
if i == 0: #Momentum
ratio = np.array((ak.Array(Y_s)/ak.Array(X_s)))
- else:
+ else: #Angle difference
ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
X_plot[j] = X_s
Y_plot[j] = ratio
@@ -123,7 +133,7 @@ for i in range(len(MC_list)):
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: #theta
+ 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)
@@ -132,17 +142,15 @@ for i in range(len(MC_list)):
ax5.scatter(X_plot[4], Y_plot[4], s = ssize)
ax6.scatter(X_plot[5], Y_plot[5], s = ssize)
- ax_list = [ax1,ax2,ax3,ax4,ax5]
- if i == 0:
- ax1.set_ylabel('rc/mc')
+ if i == 0: # for momentum
+ ax1.set_ylabel('rc/mc') #ratio
ax3.set_ylabel('rc/mc')
ax5.set_ylabel('rc/mc')
title ='ratio'
- for ax in ax_list:
- ax.set_yscale('log')
- ax.set_xscale('log')
- else:
- ax1.set_ylabel('rc-mc')
+ 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'
@@ -165,12 +173,57 @@ for i in range(len(MC_list)):
x_range = list(ax1.get_xlim())
fig.set_figwidth(20)
fig.set_figheight(10)
- ax1.set_title('%s %s %s %s events'%(title_list[i],title,config,Nevents))
+ 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[i],title,config)))
plt.close()
- ############
-
- #Correlation
+
+###################################################################################################
+ #Ratio vs momentum
+###################################################################################################
+
+ if i > 0: #for each variable theta, phi, and eta
+ for j in range(len(M_list)): #Repeat the following steps for each particle (pions,protons,electrons,neutrons,photons)
+ X = X_list[j]
+ Y = Y_list[j]
+ M_mc = M_list[j]
+ boolean_M = np.ones_like(M_mc) == 1
+ X_s = np.array(ak.flatten(X[boolean_M]))
+ Y_s = np.array(ak.flatten(Y[boolean_M]))
+ M_s = np.array(ak.flatten(M_mc))
+ ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
+ 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)
+ 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[i],config)))
+
+###################################################################################################
+ #Correlation
+###################################################################################################
+
+ #Repeat the following steps for each variable (momentum,theta,phi,eta)
X_len = ak.count(X1,axis=None)
Y_len = ak.count(Y1,axis=None)
if X_len > Y_len:
@@ -179,7 +232,8 @@ for i in range(len(MC_list)):
F_boolean = np.ones_like(X1) == 1
X_s = np.array(ak.flatten(X1[F_boolean]))
Y_s = np.array(ak.flatten(Y1[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)
else:
@@ -188,7 +242,7 @@ for i in range(len(MC_list)):
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
+ norm_h_text = []
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
@@ -211,53 +265,15 @@ for i in range(len(MC_list)):
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]))
- fig.suptitle('%s %s events'%(config,Nevents))
+ fig.suptitle('%s %s events\n DETECTOR_CONFIG: %s'%(config,Nevents,Dconfig))
plt.savefig(os.path.join(r_path, '%s_correlation_%s.png' % (title_list[i],config)))
-###############
-
- if i > 0:
- for j in range(len(X_list)):
- X = X_list[j]
- Y = Y_list[j]
- M_mc = M_list[j]
- boolean_M = np.ones_like(M_mc) == 1
- X_s = np.array(ak.flatten(X[boolean_M]))
- Y_s = np.array(ak.flatten(Y[boolean_M]))
- M_s = np.array(ak.flatten(M_mc))
- ratio = np.array((ak.Array(Y_s)-(ak.Array(X_s))))
- 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)
- # fig.suptitle('')
- 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)
- ax1.set_title('%s Difference Vs Momentum %s %s events'%(title_list[i],config,Nevents))
- plt.savefig(os.path.join(r_path, '%s_difference_vs_momentum_%s.png' % (title_list[i],config)))
-
-################
+###################################################################################################
+ #Phi vs Theta plots
+###################################################################################################
def particle_plots(boolean_particle):
+ #filtered lists w.r.t the particle
theta_mc_fil = ak.Array(theta_mc[simID][booll])[boolean_particle]
theta_rc_fil = ak.Array(theta_rc[recID][booll])[boolean_particle]
phi_mc_fil = ak.Array(phi_mc[simID][booll])[boolean_particle]
@@ -269,6 +285,7 @@ def particle_plots(boolean_particle):
F_boolean = np.ones_like(theta_rc_fil) == 1
else:
F_boolean = np.ones_like(theta_mc_fil) == 1
+ #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]))
@@ -298,7 +315,7 @@ if particle in particle_dict.keys():
particle_name = particle_dict[particle][1]
particle_plots(boolean_particle)
- plt.suptitle('%s in %s %s events'%(particle_name,config,Nevents))
+ 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,config)))
else:
for i in [[boolean_photon,'Photons'],[boolean_electron,'Electrons'],[boolean_pion,'Pions']]:
@@ -306,7 +323,7 @@ else:
particle_name = i[1]
particle_plots(boolean_particle)
- plt.suptitle('%s in %s %s events'%(particle_name,config,Nevents))
+ 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,config)))
diff --git a/benchmarks/dis/dis.sh b/benchmarks/dis/dis.sh
index 2c25719a64a05c61a3c3e3b7cead3597dfc56d10..7b5d952bbcaa697e80fb2551e81e57a257627882 100755
--- a/benchmarks/dis/dis.sh
+++ b/benchmarks/dis/dis.sh
@@ -128,7 +128,7 @@ if [[ "$?" -ne "0" ]] ; then
exit 1
fi
-python benchmarks/dis/analysis/truth_reconstruction.py --rec_file ${REC_FILE} --config ${PLOT_TAG} --results_path ${RESULTS_PATH} --nevents ${JUGGLER_N_EVENTS}
+python benchmarks/dis/analysis/truth_reconstruction.py --rec_file ${REC_FILE} --config ${PLOT_TAG}_${DETECTOR_CONFIG} --results_path ${RESULTS_PATH} --nevents ${JUGGLER_N_EVENTS}
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running truth_reconstruction script"
exit 1
diff --git a/benchmarks/single/analysis/analyze.cxx b/benchmarks/single/analysis/analyze.cxx
index 4070868e5e475f4b33a38c8b163b2f201ca96d3b..8151c94a5a0437ef0890d9d415981bcad49b82c4 100644
--- a/benchmarks/single/analysis/analyze.cxx
+++ b/benchmarks/single/analysis/analyze.cxx
@@ -8,7 +8,7 @@
int analyze(std::string file)
{
// open dataframe
- ROOT::RDataFrame df("events", file, {"GeneratedParticles", "ReconstructedParticles"});
+ ROOT::RDataFrame df("events", file, {"GeneratedParticles", "ReconstructedChargedParticles"});
// count total events
auto count = df.Count();
@@ -21,7 +21,7 @@ int analyze(std::string file)
auto d = df
.Define("n_tracks_gen", n_tracks, {"GeneratedParticles"})
- .Define("n_tracks_rec", n_tracks, {"ReconstructedParticles"})
+ .Define("n_tracks_rec", n_tracks, {"ReconstructedChargedParticles"})
;
auto stats_n_tracks_gen = d.Stats("n_tracks_gen");
diff --git a/benchmarks/single/analyze.sh b/benchmarks/single/analyze.sh
index 6271888dfc00173af9d9bdf661041970d6415e7e..681e32ccc91705d5b46a3fc6483b04b101d442ac 100644
--- a/benchmarks/single/analyze.sh
+++ b/benchmarks/single/analyze.sh
@@ -10,7 +10,7 @@ if [[ "$?" -ne "0" ]] ; then
exit 1
fi
-python benchmarks/dis/analysis/truth_reconstruction.py --rec_file ${JUGGLER_REC_FILE} --config ${JUGGLER_FILE_NAME_TAG} --results_path ${RESULTS_PATH} --nevents ${JUGGLER_N_EVENTS}
+python benchmarks/dis/analysis/truth_reconstruction.py --rec_file ${JUGGLER_REC_FILE} --config ${JUGGLER_FILE_NAME_TAG}_${DETECTOR_CONFIG} --results_path ${RESULTS_PATH} --nevents ${JUGGLER_N_EVENTS}
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running truth_reconstruction script"
exit 1
diff --git a/benchmarks/single/pi-_1GeV_3to45deg.steer b/benchmarks/single/pi-_1GeV_3to45deg.steer
index f547e2e116d6b7a7bb1896dae75ee09090195fb5..40d767e5e0037f91e293e53dd3fc7ad5c468a2f6 100644
--- a/benchmarks/single/pi-_1GeV_3to45deg.steer
+++ b/benchmarks/single/pi-_1GeV_3to45deg.steer
@@ -3,7 +3,7 @@ from g4units import mm, GeV, MeV, degree
SIM = DD4hepSimulation()
SIM.gun.energy = 1*GeV
-SIM.gun.particle = "e-"
+SIM.gun.particle = "pi-"
SIM.gun.position = (0.0, 0.0, 0.0)
SIM.gun.direction = (0.0, 0.0, 1.0)
SIM.gun.distribution = "cos(theta)"