From db51dba6904b8ee4850cebdf23b1e56a6fb70dd1 Mon Sep 17 00:00:00 2001
From: "lintao@ihep.ac.cn" <lintao@ihep.ac.cn>
Date: Wed, 26 Mar 2025 09:52:50 +0000
Subject: [PATCH] WIP: improve the MCTruth handling in detector simulation
---
Detector/DetCRD/scripts/TDR_o1_v01/sim.py | 4 +++
.../src/Edm4hepWriterAnaElemTool.cpp | 35 ++++++++++++++++---
.../DetSimAna/src/Edm4hepWriterAnaElemTool.h | 7 ++--
3 files changed, 39 insertions(+), 7 deletions(-)
diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/sim.py b/Detector/DetCRD/scripts/TDR_o1_v01/sim.py
index 8bcb3eeb..89159f9d 100644
--- a/Detector/DetCRD/scripts/TDR_o1_v01/sim.py
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/sim.py
@@ -93,6 +93,10 @@ detsimalg.AnaElems = [
]
detsimalg.RootDetElem = "WorldDetElemTool"
+from Configurables import Edm4hepWriterAnaElemTool
+detsim_anatool = Edm4hepWriterAnaElemTool("Edm4hepWriterAnaElemTool")
+detsim_anatool.IsTrk2Primary = False # True: primary; False: ancestor
+
from Configurables import TimeProjectionChamberSensDetTool
tpc_sensdettool = TimeProjectionChamberSensDetTool("TimeProjectionChamberSensDetTool")
tpc_sensdettool.TypeOption = 1
diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp
index 5f20cc77..3139e2ef 100644
--- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp
+++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp
@@ -282,7 +282,7 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) {
int pritrkid = m_track2primary[trackID];
if (pritrkid <= 0) {
- error() << "Failed to find the primary track for trackID #" << trackID << endmsg;
+ error() << "Failed to find the primary (or ancestor) track for trackID #" << trackID << endmsg;
pritrkid = 1;
}
@@ -363,6 +363,16 @@ Edm4hepWriterAnaElemTool::PreUserTrackingAction(const G4Track* track) {
int curparid = track->GetParentID();
int pritrkid = curparid;
+ // if it is ancestor mode, then we need to check whether current track is associated
+ // with MCParticle or not.
+ // If it is associated, then we will record this info.
+ if (not m_istrk2primary.value()) {
+ auto trackinfo = dynamic_cast<CommonUserTrackInfo*>(track->GetUserInformation());
+ if (trackinfo) {
+ curparid = 0; // force to zero, so that it will be set as primary
+ }
+ }
+
// try to find the primary track id from the parent track id.
if (curparid) {
auto it = m_track2primary.find(curparid);
@@ -429,6 +439,11 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
auto trackinfo = dynamic_cast<CommonUserTrackInfo*>(track->GetUserInformation());
if (trackinfo) {
idxedm4hep = trackinfo->idxEdm4hep();
+ // if the current track is a secondary track and it is associated
+ // with a MCParticle, we need to update the event level user info
+ if (idxedm4hep) {
+ m_userinfo->setIdxG4Track2Edm4hep(curtrkid, idxedm4hep);
+ }
}
}
@@ -475,6 +490,7 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
// flags
bool is_decay = false;
+ G4ThreeVector decay_pos; // only valid when is decay
size_t nSeco = secondaries->size();
for (size_t i = 0; i < nSeco; ++i) {
@@ -506,7 +522,9 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
<< " momentum: " << sectrk->GetMomentum() //
<< endmsg;
- if (creatorProcess==proc_decay) is_decay = true;
+ if (creatorProcess==proc_decay) {
+ is_decay = true;
+ }
// create secondaries in MC particles
// todo: convert the const collection to non-const
@@ -524,6 +542,11 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
double y=sec_init_pos.y()/CLHEP::mm;
double z=sec_init_pos.z()/CLHEP::mm;
+ // if the particle decay, record the position
+ if (is_decay) {
+ decay_pos = sec_init_pos;
+ }
+
const G4ThreeVector& sec_init_mom = sectrk->GetMomentum();
double px=sec_init_mom.x()/CLHEP::GeV;
double py=sec_init_mom.y()/CLHEP::GeV;
@@ -548,8 +571,12 @@ Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) {
// now update
if (is_decay) {
- primary_particle.setDecayedInTracker(is_decay);
- primary_particle.setDecayedInCalorimeter(is_decay);
+ bool is_in_tracker = decay_pos.perp() < R && std::fabs(decay_pos.z()) < Z;
+ if (is_in_tracker) {
+ primary_particle.setDecayedInTracker(is_decay);
+ } else {
+ primary_particle.setDecayedInCalorimeter(is_decay);
+ }
}
}
diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
index 7b9b6e71..b386b122 100644
--- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
+++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h
@@ -77,11 +77,12 @@ private:
Gaudi::Property<double> m_sectrk_Ek{this, "SecTrackEk", 100., "Ek (MeV) threshold to record a secondary track"};
Gaudi::Property<double> m_sectrk_rho{this, "SecTrackRho", 1830., "rho (mm) threshold to record a secondary track"};
Gaudi::Property<double> m_sectrk_z{this, "SecTrackZ", 2900., "+/- z (mm) threshold to record a secondary track"};
-
+
+ Gaudi::Property<bool> m_istrk2primary{this, "IsTrk2Primary", true, "For m_track2primary, the value is primary or ancestor"};
private:
- // in order to associate the hit contribution with the primary track,
+ // in order to associate the hit contribution with the primary track or ancestor track,
// we have a bookkeeping of every track.
- // The primary track will assign the same key/value.
+ // The primary or ancestor track will assign the same key/value.
// Following is an example:
// 1 -> 1,
--
GitLab