diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp index 9de89e9447aec654c51a1567f9dd482cfcce00fe..ca3bf7b63705b2b7ecf841a30d9e2b089c31cfff 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.cpp @@ -30,6 +30,13 @@ void Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) { msg() << "Event " << anEvent->GetEventID() << endmsg; + // event info + m_userinfo = nullptr; + if (anEvent->GetUserInformation()) { + m_userinfo = dynamic_cast<CommonUserEventInfo*>(anEvent->GetUserInformation()); + m_userinfo->dumpIdxG4Track2Edm4hep(); + } + auto mcGenCol = m_mcParGenCol.get(); mcCol = m_mcParCol.createAndPut(); @@ -50,7 +57,7 @@ Edm4hepWriterAnaElemTool::BeginOfEventAction(const G4Event* anEvent) { newparticle.setColorFlow (mcGenParticle.getColorFlow()); } - msg() << "mcCol size: " << mcCol->size() << endmsg; + msg() << "mcCol size (original) : " << mcCol->size() << endmsg; // reset m_track2primary.clear(); @@ -257,7 +264,12 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { pritrkid = 1; } - edm_trk_hit.setMCParticle(mcCol->at(pritrkid-1)); + if (m_userinfo) { + auto idxedm4hep = m_userinfo->idxG4Track2Edm4hep(pritrkid); + if (idxedm4hep != -1) { + edm_trk_hit.setMCParticle(mcCol->at(idxedm4hep)); + } + } if (pritrkid != trackID) { // If the track is a secondary, then the primary track id and current track id is different @@ -297,7 +309,12 @@ Edm4hepWriterAnaElemTool::EndOfEventAction(const G4Event* anEvent) { pritrkid = 1; } - edm_calo_contrib.setParticle(mcCol->at(pritrkid-1)); // todo + if (m_userinfo) { + auto idxedm4hep = m_userinfo->idxG4Track2Edm4hep(pritrkid); + if (idxedm4hep != -1) { + edm_calo_contrib.setParticle(mcCol->at(idxedm4hep)); // todo + } + } edm_calo_hit.addToContributions(edm_calo_contrib); } } @@ -345,17 +362,27 @@ void Edm4hepWriterAnaElemTool::PostUserTrackingAction(const G4Track* track) { int curtrkid = track->GetTrackID(); // starts from 1 int curparid = track->GetParentID(); + int idxedm4hep = -1; - if (curparid == 0) { + if (m_userinfo) { + idxedm4hep = m_userinfo->idxG4Track2Edm4hep(curtrkid); + } + + if (curparid == 0) { // Primary Track // select the primary tracks (parentID == 0) // auto mcCol = m_mcParCol.get(); - if (curtrkid-1>=mcCol->size()) { + if (idxedm4hep == -1) { + error () << "Failed to get idxedm4hep according to the g4track id " << curtrkid << endmsg; + return; + } + + if (idxedm4hep>=mcCol->size()) { error() << "out of range: curtrkid is " << curtrkid << " while the MCParticle size is " << mcCol->size() << endmsg; return; } - auto primary_particle = mcCol->at(curtrkid-1); + auto primary_particle = mcCol->at(idxedm4hep); const G4ThreeVector& stop_pos = track->GetPosition(); edm4hep::Vector3d endpoint(stop_pos.x()/CLHEP::mm, diff --git a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h index 937869b4a7a1f19d0c26c7828976934f12e8e8c2..b290c40bd8f10545ce16d1c8588eb12f5aa54998 100644 --- a/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h +++ b/Simulation/DetSimAna/src/Edm4hepWriterAnaElemTool.h @@ -6,6 +6,7 @@ #include "GaudiKernel/AlgTool.h" #include "k4FWCore/DataHandle.h" #include "DetSimInterface/IAnaElemTool.h" +#include "DetSimInterface/CommonUserEventInfo.hh" #include "edm4hep/MCParticleCollection.h" #include "edm4hep/SimTrackerHitCollection.h" @@ -140,6 +141,7 @@ private: std::map<int, int> m_track2primary; + CommonUserEventInfo* m_userinfo = nullptr; }; #endif diff --git a/Simulation/DetSimCore/src/G4PrimaryCnvTool.cpp b/Simulation/DetSimCore/src/G4PrimaryCnvTool.cpp index b03a403e6d2eb013c0b0bb8d88dfebe7623eb3f8..35711ba17a37a365cf22c1232790c552f6433456 100644 --- a/Simulation/DetSimCore/src/G4PrimaryCnvTool.cpp +++ b/Simulation/DetSimCore/src/G4PrimaryCnvTool.cpp @@ -5,10 +5,19 @@ #include "G4IonTable.hh" #include "G4ParticleDefinition.hh" +#include "DetSimInterface/CommonUserEventInfo.hh" + DECLARE_COMPONENT(G4PrimaryCnvTool) bool G4PrimaryCnvTool::mutate(G4Event* anEvent) { + // create a event info + auto eventinfo = new CommonUserEventInfo(); + anEvent->SetUserInformation(eventinfo); + + int idxG4 = 0; // valid: [1, N+1) + int idxEdm4hep = -1; // valid: [0, N) + auto mcCol = m_mcParCol.get(); info() << "Start a new event: " << endmsg; for ( auto p : *mcCol ) { @@ -20,11 +29,19 @@ bool G4PrimaryCnvTool::mutate(G4Event* anEvent) { } info() << " ]); " << endmsg; + // idx in mc particle collection + ++idxEdm4hep; + // only the GeneratorStatus == 1 is used. if (p.getGeneratorStatus() != 1) { continue; } + // idx in g4 collection + ++idxG4; + + eventinfo->setIdxG4Track2Edm4hep(idxG4, idxEdm4hep); + // vertex const edm4hep::Vector3d& vertex = p.getVertex(); double t = p.getTime()*CLHEP::ns; diff --git a/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh b/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh index 733c66669b02fcb9c39ad4b0ae38342e8821583c..3319598a270d0d4cb38607d87b6c4fd3b2a85ad4 100644 --- a/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh +++ b/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh @@ -14,6 +14,7 @@ */ #include "G4VUserEventInformation.hh" +#include <map> class CommonUserEventInfo: public G4VUserEventInformation { public: @@ -24,6 +25,15 @@ public: public: virtual void Print() const; + // set the relationship between idx in geant4 and idx in mc particle collection. + // idxG4: G4 track ID (starts from 1) + // idxEdm4hep: index in MC Particle collection (starts from 0) + bool setIdxG4Track2Edm4hep(int idxG4, int idxEdm4hep); + int idxG4Track2Edm4hep(int idxG4) const; + void dumpIdxG4Track2Edm4hep() const; + +private: + std::map<int, int> m_g4track_to_edm4hep; }; #endif diff --git a/Simulation/DetSimInterface/src/CommonUserEventInfo.cc b/Simulation/DetSimInterface/src/CommonUserEventInfo.cc index 3f7d2af760f7fa5fd2071e4a38378f35a0187afb..7992c69c897967e728040b3602a472663f604c6f 100644 --- a/Simulation/DetSimInterface/src/CommonUserEventInfo.cc +++ b/Simulation/DetSimInterface/src/CommonUserEventInfo.cc @@ -1,4 +1,5 @@ #include "DetSimInterface/CommonUserEventInfo.hh" +#include <iostream> CommonUserEventInfo::CommonUserEventInfo() { @@ -12,3 +13,41 @@ void CommonUserEventInfo::Print() const { } + +bool +CommonUserEventInfo::setIdxG4Track2Edm4hep(int idxG4, int idxEdm4hep) { + auto it = m_g4track_to_edm4hep.find(idxG4); + + // if already exists, return false + if (it != m_g4track_to_edm4hep.end()) { + return false; + } + + m_g4track_to_edm4hep[idxG4] = idxEdm4hep; + + return true; +} + +int +CommonUserEventInfo::idxG4Track2Edm4hep(int idxG4) const { + int ret = -1; + + auto it = m_g4track_to_edm4hep.find(idxG4); + + // if found + if (it != m_g4track_to_edm4hep.end()) { + ret = it->second; + } + + return ret; +} + +void +CommonUserEventInfo::dumpIdxG4Track2Edm4hep() const { + std::cout << "---- Dumping IdxG4Track2Edm4hep: " + << " total size: " << m_g4track_to_edm4hep.size() + << std::endl; + for (auto [idxG4, idxE4]: m_g4track_to_edm4hep) { + std::cout << " - " << idxG4 << " -> " << idxE4 << std::endl; + } +}