From 5ed37b550d17d185a0c416c63df9538ca99bd461 Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Tue, 17 Aug 2021 14:19:44 +0000
Subject: [PATCH] Resolve "Add true decaying particle info onto event display"
---
.../imaging_ecal/options/imaging_2dcluster.py | 9 +++++-
benchmarks/imaging_ecal/run_emcal_barrel.sh | 2 +-
benchmarks/imaging_ecal/run_imcal_pion0.sh | 2 +-
.../imaging_ecal/scripts/draw_cluster.py | 26 ++++++++++++++++-
benchmarks/imaging_ecal/scripts/utils.py | 29 +++++++++++++++++++
5 files changed, 64 insertions(+), 4 deletions(-)
diff --git a/benchmarks/imaging_ecal/options/imaging_2dcluster.py b/benchmarks/imaging_ecal/options/imaging_2dcluster.py
index 6802e07d..d2a038ca 100644
--- a/benchmarks/imaging_ecal/options/imaging_2dcluster.py
+++ b/benchmarks/imaging_ecal/options/imaging_2dcluster.py
@@ -6,6 +6,7 @@ from GaudiKernel.SystemOfUnits import mm, MeV, rad, ns
from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
from Configurables import Jug__Digi__CalorimeterHitDigi as CalHitDigi
from Configurables import Jug__Reco__CalorimeterHitReco as CalHitReco
from Configurables import Jug__Reco__CalorimeterHitsEtaPhiProjector as CalHitsProj
@@ -37,6 +38,12 @@ podioevent = EICDataSvc('EventDataSvc', inputs=kwargs['input'].split(','), Outpu
podin = PodioInput('PodioReader', collections=['mcparticles', 'EcalBarrelHits'], OutputLevel=DEBUG)
podout = PodioOutput('out', filename=kwargs['output'])
+
+mccopier = MCCopier('MCCopier',
+ # OutputLevel=DEBUG,
+ inputCollection='mcparticles',
+ outputCollection='mcparticles2')
+
# use the same daq_setting for digi/reco pair
imcal_barrel_daq = dict(
dynamicRangeADC=3.*MeV,
@@ -86,7 +93,7 @@ imcal_barrel_clreco = RecoCoG('imcal_barrel_clreco',
podout.outputCommands = ['keep *']
ApplicationMgr(
- TopAlg=[podin,
+ TopAlg=[podin, mccopier,
imcal_barrel_digi, imcal_barrel_reco,
imcal_barrel_merger, imcal_barrel_cl, imcal_barrel_clreco,
podout],
diff --git a/benchmarks/imaging_ecal/run_emcal_barrel.sh b/benchmarks/imaging_ecal/run_emcal_barrel.sh
index 73bae02c..0dca31e5 100644
--- a/benchmarks/imaging_ecal/run_emcal_barrel.sh
+++ b/benchmarks/imaging_ecal/run_emcal_barrel.sh
@@ -87,7 +87,7 @@ ls -lh ${CB_EMCAL_GEN_FILE}
# Run geant4 simulations
npsim --runType batch \
-v WARNING \
- --part.minimalKineticEnergy "1*TeV" \
+ --part.minimalKineticEnergy "0.5*MeV" \
--numberOfEvents ${CB_EMCAL_NUMEV} \
--compactFile ${CB_EMCAL_COMPACT_PATH} \
--inputFiles ${CB_EMCAL_GEN_FILE} \
diff --git a/benchmarks/imaging_ecal/run_imcal_pion0.sh b/benchmarks/imaging_ecal/run_imcal_pion0.sh
index 760de223..9506c117 100644
--- a/benchmarks/imaging_ecal/run_imcal_pion0.sh
+++ b/benchmarks/imaging_ecal/run_imcal_pion0.sh
@@ -83,7 +83,7 @@ ls -lh ${CB_EMCAL_GEN_FILE}
# Run geant4 simulations
npsim --runType batch \
-v WARNING \
- --part.minimalKineticEnergy "1*TeV" \
+ --part.minimalKineticEnergy "0.5*MeV" \
--numberOfEvents ${CB_EMCAL_NUMEV} \
--compactFile ${CB_EMCAL_COMPACT_PATH} \
--inputFiles ${CB_EMCAL_GEN_FILE} \
diff --git a/benchmarks/imaging_ecal/scripts/draw_cluster.py b/benchmarks/imaging_ecal/scripts/draw_cluster.py
index 3c896f2c..f810b6a9 100644
--- a/benchmarks/imaging_ecal/scripts/draw_cluster.py
+++ b/benchmarks/imaging_ecal/scripts/draw_cluster.py
@@ -5,6 +5,10 @@
Author: Chao Peng (ANL)
Date: 04/30/2021
+
+ Added true decaying particles on eta-phi plane projection plot
+ Author: Jihee Kim (ANL)
+ Data: 08/06/2021
'''
import os
@@ -143,9 +147,26 @@ if __name__ == '__main__':
df['theta'] = np.arccos(df['z'].values/df['r'].values)*1000.
df['eta'] = -np.log(np.tan(df['theta'].values/1000./2.))
+ # Read all mc particles
+ dfallmcp = get_all_mcp(args.file, args.iev, 'mcparticles2')
+ pdgbase = ROOT.TDatabasePDG()
+ # Select decaying particles
+ dftemp = dfallmcp[dfallmcp['g4Parent'] == 1.0]
+ if len(dftemp) > 0:
+ dfdecaymcp = dftemp.copy()
+ for iptl in [0, len(dfdecaymcp) - 1]:
+ infoptl = pdgbase.GetParticle(int(dfdecaymcp['pid'].iloc[iptl]))
+ print("{} Decaying particle = {}, pdgcode = {}, charge = {}, mass = {}"\
+ .format(iptl, infoptl.GetName(), infoptl.PdgCode(), infoptl.Charge(), infoptl.Mass()))
+ # Calculate geometric variables of decaying particles
+ dfdecaymcp['r'] = np.sqrt(dfdecaymcp['vex'].values**2 + dfdecaymcp['vey'].values**2 + dfdecaymcp['vez'].values**2)
+ dfdecaymcp['phi'] = np.arctan2(dfdecaymcp['vey'].values, dfdecaymcp['vex'].values)*1000.
+ dfdecaymcp['theta'] = np.arccos(dfdecaymcp['vez'].values/dfdecaymcp['r'].values)*1000.
+ dfdecaymcp['eta'] = -np.log(np.tan(dfdecaymcp['theta'].values/1000./2.))
+
# truth
dfmcp = get_mcp_simple(args.file, args.iev, 'mcparticles2').iloc[0]
- pdgbase = ROOT.TDatabasePDG()
+ #pdgbase = ROOT.TDatabasePDG()
inpart = pdgbase.GetParticle(int(dfmcp['pid']))
print("Incoming particle = {}, pdgcode = {}, charge = {}, mass = {}"\
.format(inpart.GetName(), inpart.PdgCode(), inpart.Charge(), inpart.Mass()))
@@ -216,6 +237,9 @@ if __name__ == '__main__':
ax, sm = draw_heatmap(axs[0], df['eta'].values, df['phi'].values, weights=df['energy'].values,
bins=(np.arange(*eta_rg, step=args.topo_size/1000.), np.arange(*phi_rg, step=args.topo_size)),
cmap=cmap, cmin=0., pc_kw=dict(alpha=0.8, edgecolor='k'))
+ # draw true decaying particle position
+ if len(dftemp) > 0:
+ ax.scatter(dfdecaymcp['eta'].values, dfdecaymcp['phi'].values, marker='x', color='red', s=22**2, linewidth=5.0)
ax.set_ylabel(r'$\phi$ (mrad)', fontsize=32)
ax.set_xlabel(r'$\eta$', fontsize=32)
diff --git a/benchmarks/imaging_ecal/scripts/utils.py b/benchmarks/imaging_ecal/scripts/utils.py
index ca58634a..19f3722d 100644
--- a/benchmarks/imaging_ecal/scripts/utils.py
+++ b/benchmarks/imaging_ecal/scripts/utils.py
@@ -3,6 +3,10 @@
Author: Chao Peng (ANL)
Date: 04/30/2021
+
+ Added all mc particle info to obtain decaying particles
+ Author: Jihee Kim (ANL)
+ Date: 08/06/2021
'''
import os
@@ -79,6 +83,31 @@ def get_mcp_simple(path, evnums=None, branch='mcparticles2'):
idb += 1
return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status'])
+#######################################
+# read all mc particles from root file
+#######################################
+def get_all_mcp(path, evnums=None, branch='mcparticles2'):
+ f = ROOT.TFile(path)
+ events = f.events
+ if evnums is None:
+ evnums = np.arange(events.GetEntries())
+ elif isinstance(evnums, int):
+ evnums = [evnums]
+
+ dbuf = np.zeros(shape=(2000*len(evnums), 10))
+ idb = 0
+ for iev in evnums:
+ if iev >= events.GetEntries():
+ print('Error: event {:d} is out of range (0 - {:d})'.format(iev, events.GetEntries() - 1))
+ continue
+
+ events.GetEntry(iev)
+ # extract mc particle data
+ for ptl in getattr(events, branch):
+ dbuf[idb] = (iev, ptl.ps.x, ptl.ps.y, ptl.ps.z, ptl.pdgID, ptl.status, ptl.g4Parent, ptl.ve.x, ptl.ve.y, ptl.ve.z)
+ idb += 1
+
+ return pd.DataFrame(data=dbuf[:idb], columns=['event', 'px', 'py', 'pz', 'pid', 'status', 'g4Parent', 'vex', 'vey', 'vez'])
# read hits data from root file
def get_hits_data(path, evnums=None, branch='RecoEcalBarrelHits'):
--
GitLab