diff --git a/Analysis/CMakeLists.txt b/Analysis/CMakeLists.txt
index 7e666e5bb77e848d9218a4ee66b5018fcd871a94..58590ac4f319bc8f4bc40f2fd0f093d97927a4a4 100644
--- a/Analysis/CMakeLists.txt
+++ b/Analysis/CMakeLists.txt
@@ -1,3 +1,4 @@
add_subdirectory(TotalInvMass)
add_subdirectory(TrackInspect)
+add_subdirectory(DumpEvent)
diff --git a/Analysis/DumpEvent/CMakeLists.txt b/Analysis/DumpEvent/CMakeLists.txt
index 3c5f697d2e4355b72a6ee88c16e271edba076353..d432b646a874c1a640a2459ccbbec64d92f53527 100644
--- a/Analysis/DumpEvent/CMakeLists.txt
+++ b/Analysis/DumpEvent/CMakeLists.txt
@@ -1,5 +1,6 @@
gaudi_add_module(DumpEvent
SOURCES src/DumpMCParticleAlg.cpp
+ src/DumpSimHitAlg.cpp
#src/DumpHitAlg.cpp
src/DumpTrackAlg.cpp
#src/DumpCalorimeterAlg.cpp
@@ -10,6 +11,7 @@ gaudi_add_module(DumpEvent
${CLHEP_LIBRARIES}
${DD4hep_COMPONENT_LIBRARIES}
DetInterface
+ k4FWCore::k4FWCore
)
install(TARGETS DumpEvent
diff --git a/Analysis/DumpEvent/src/DumpSimHitAlg.cpp b/Analysis/DumpEvent/src/DumpSimHitAlg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5c9a3d1f22f99b75f7d4e987e4befa9e5066e30e
--- /dev/null
+++ b/Analysis/DumpEvent/src/DumpSimHitAlg.cpp
@@ -0,0 +1,76 @@
+/*
+ * Description:
+ * Dump the simulated information.
+ *
+ * Author:
+ * Tao Lin <lintao AT ihep.ac.cn>
+ */
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/Algorithm.h"
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/CaloHitContributionCollection.h"
+
+#include "GaudiKernel/NTuple.h"
+
+class DumpSimHitAlg: public Algorithm {
+public:
+
+ DumpSimHitAlg(const std::string& name, ISvcLocator* pSvcLocator);
+
+ // Three mandatory member functions of any algorithm
+ StatusCode initialize() override;
+ StatusCode execute() override;
+ StatusCode finalize() override;
+
+private:
+ // - collection MCParticleG4: the simulated particles in Geant4
+ DataHandle<edm4hep::MCParticleCollection> m_mcParCol{"MCParticle",
+ Gaudi::DataHandle::Reader, this};
+
+ // Dedicated collections for CEPC
+ DataHandle<edm4hep::SimTrackerHitCollection> m_VXDCol{"VXDCollection",
+ Gaudi::DataHandle::Reader, this};
+
+};
+
+DECLARE_COMPONENT( DumpSimHitAlg )
+
+DumpSimHitAlg::DumpSimHitAlg(const std::string& name, ISvcLocator* pSvcLocator)
+: Algorithm(name, pSvcLocator) {
+
+}
+
+StatusCode DumpSimHitAlg::initialize() {
+ return StatusCode::SUCCESS;
+}
+
+StatusCode DumpSimHitAlg::execute() {
+ auto mcCol = m_mcParCol.get();
+
+ auto vxdCol = m_VXDCol.get();
+
+ for (auto hit: *vxdCol) {
+ auto mcparticle = hit.getMCParticle();
+
+ if (mcparticle.getGeneratorStatus() != 1) {
+ error() << "Found generator status is not 1 for hit. " << endmsg;
+ }
+
+ info() << "hit -> "
+ << " mcparticle ("
+ << " ID: " << mcparticle.getObjectID().index << "; "
+ << " generator status: " << mcparticle.getGeneratorStatus() << "; "
+ << " simulator status: " << mcparticle.getSimulatorStatus() << ") "
+ << endmsg;
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode DumpSimHitAlg::finalize() {
+ return StatusCode::SUCCESS;
+}
diff --git a/Examples/options/dump_simhit.py b/Examples/options/dump_simhit.py
new file mode 100644
index 0000000000000000000000000000000000000000..1c699b7a4f98ca5860e6f7a77763f2c91e16eb07
--- /dev/null
+++ b/Examples/options/dump_simhit.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+
+from Gaudi.Configuration import *
+
+##############################################################################
+# Geometry Svc
+##############################################################################
+
+# geometry_option = "CepC_v4-onlyTracker.xml"
+# geometry_option = "CepC_v4-onlyVXD.xml"
+geometry_option = "CepC_v4.xml"
+
+if not os.getenv("DETCEPCV4ROOT"):
+ print("Can't find the geometry. Please setup envvar DETCEPCV4ROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCEPCV4ROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+##############################################################################
+# Event Data Svc
+##############################################################################
+
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc", input="test-detsim10.root")
+
+##############################################################################
+# NTuple Svc
+##############################################################################
+
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = ["MyTuples DATAFILE='result.root' OPT='NEW' TYP='ROOT'"]
+
+##############################################################################
+# DumpAlg
+##############################################################################
+
+from Configurables import DumpSimHitAlg
+alg = DumpSimHitAlg("DumpSimHitAlg")
+
+from Configurables import PodioInput
+podioinput = PodioInput("PodioReader", collections=[
+ "MCParticle",
+ "VXDCollection"
+ ])
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+ApplicationMgr( TopAlg = [podioinput, alg],
+ EvtSel = 'NONE',
+ EvtMax = -1,
+ ExtSvc = [dsvc, ntsvc],
+ HistogramPersistency = "ROOT",
+ # OutputLevel=DEBUG
+)
diff --git a/Generator/src/GenPrinter.cpp b/Generator/src/GenPrinter.cpp
index 48c6d499451ec24b1eb9b0874b971cdc0678f3b2..4840bee179b367531ab4c2f32cf14ec6d3c7eebe 100644
--- a/Generator/src/GenPrinter.cpp
+++ b/Generator/src/GenPrinter.cpp
@@ -6,26 +6,39 @@ DECLARE_COMPONENT(GenPrinter)
bool GenPrinter::mutate(MyHepMC::GenEvent& event){
std::cout << "print mc info for event "<< event.getID() << ", mc size ="<< event.m_mc_vec.size() << std::endl;
for ( int i =0; i < event.m_mc_vec.size(); i++ ) {
- auto p = event.m_mc_vec.at(i);
- std::cout<< "PDG :"<< p.getPDG ()<<std::endl
- << "id :"<< p.id ()<<std::endl
- << "ID :"<< p.getObjectID().index <<std::endl
- << "GeneratorStatus :"<< p.getGeneratorStatus ()<<std::endl
- << "SimulatorStatus :"<< p.getSimulatorStatus ()<<std::endl
- << "Charge :"<< p.getCharge ()<<std::endl
- << "Time :"<< p.getTime ()<<std::endl
- << "Mass :"<< p.getMass ()<<std::endl
- << "Vertex :"<< p.getVertex ()<<std::endl
- << "Endpoint :"<< p.getEndpoint ()<<std::endl
- << "Momentum :"<< p.getMomentum ()<<std::endl
- << "MomentumAtEndpoint:"<< p.getMomentumAtEndpoint()<<std::endl
- << "Spin :"<< p.getSpin ()<<std::endl
- << "ColorFlow :"<< p.getColorFlow ()<<std::endl
- << "Parent size :"<< p.parents_size ()<<std::endl
- << "Daughter size :"<< p.daughters_size ()<<std::endl;
- //for(unsigned int j=0; j<p.parents_size(); j++) std::cout << " for parent: "<< j << ",PDG="<< p.getParents(j).getPDG() << ",id=:"<< p.getParents(j).id()<<std::endl;
- for (auto it = p.parents_begin(), end = p.parents_end(); it != end ; ++it ) std::cout << " for parent ,PDG="<< it->getPDG() << ",id=:"<< it->getObjectID().index<<std::endl;
+ auto p = event.m_mc_vec.at(i);
+ std::cout<< "PDG :"<< p.getPDG ()<<std::endl
+ << "id :"<< p.id ()<<std::endl
+ << "ID :"<< p.getObjectID().index <<std::endl
+ << "GeneratorStatus :"<< p.getGeneratorStatus ()<<std::endl
+ << "SimulatorStatus :"<< p.getSimulatorStatus ()<<std::endl
+ << "Charge :"<< p.getCharge ()<<std::endl
+ << "Time :"<< p.getTime ()<<std::endl
+ << "Mass :"<< p.getMass ()<<std::endl
+ << "Vertex :"<< p.getVertex ()<<std::endl
+ << "Endpoint :"<< p.getEndpoint ()<<std::endl
+ << "Momentum :"<< p.getMomentum ()<<std::endl
+ << "MomentumAtEndpoint:"<< p.getMomentumAtEndpoint()<<std::endl
+ << "Spin :"<< p.getSpin ()<<std::endl
+ << "ColorFlow :"<< p.getColorFlow ()<<std::endl
+ << "Parent size :"<< p.parents_size ()<<std::endl
+ << "Daughter size :"<< p.daughters_size ()<<std::endl;
+ //for(unsigned int j=0; j<p.parents_size(); j++) std::cout << " for parent: "<< j << ",PDG="<< p.getParents(j).getPDG() << ",id=:"<< p.getParents(j).id()<<std::endl;
+ for (auto it = p.parents_begin(), end = p.parents_end(); it != end ; ++it ) {
+ std::cout << " - parent, PDG=" << it->getPDG()
+ << ", id=" << it->id()
+ << ", ID=" << it->getObjectID().index
+ << std::endl;
+ }
+
+ for (auto it = p.daughters_begin(), end = p.daughters_end(); it != end ; ++it ) {
+ std::cout << " - daughter, PDG=" << it->getPDG()
+ << ", id=" << it->id()
+ << ", ID=" << it->getObjectID().index
+ << std::endl;
+ }
}
+
return true;
}
diff --git a/Generator/src/StdHepRdr.cpp b/Generator/src/StdHepRdr.cpp
index b56342b0bfb25de4eb4ceee6b6d6362004bcb1bd..4255e4bd314f3549f0bea3f56be3d4e470221c5a 100644
--- a/Generator/src/StdHepRdr.cpp
+++ b/Generator/src/StdHepRdr.cpp
@@ -36,13 +36,13 @@ bool StdHepRdr::mutate(MyHepMC::GenEvent& event){
m_processed_event ++;
int n_mc = mc_vec->getNumberOfElements();
//std::cout<<"Debug: Read event :"<< m_processed_event <<", mc size :"<< n_mc <<std::endl;
- std::map<int, int> pmcid_lmcid;
+ std::map<int, int> pmcid_lmcid; // mapping between the obj idx in edm4hep and the idx in stdhep
for (int i=0; i < n_mc; i++){
MCParticleImpl* mc = (MCParticleImpl*) mc_vec->getElementAt(i);
//std::cout<<"At mc :"<< i <<std::endl;
auto mcp = event.m_mc_vec.create();
pmcid_lmcid.insert(std::pair<int, int>(mc->id(),i));
- //std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl;
+ // std::cout<<"map<id,i>:"<<mc->id()<<","<< i <<std::endl;
mcp.setPDG (mc->getPDG());
mcp.setGeneratorStatus (mc->getGeneratorStatus());
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 c756412f62dccf71d5628191b10f4a4e6f555e26..35711ba17a37a365cf22c1232790c552f6433456 100644
--- a/Simulation/DetSimCore/src/G4PrimaryCnvTool.cpp
+++ b/Simulation/DetSimCore/src/G4PrimaryCnvTool.cpp
@@ -5,24 +5,43 @@
#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 ) {
- info() << p.getObjectID().index << " : [";
+ info() << " gen track: " << p.getObjectID().index
+ << " : (status: " << p.getGeneratorStatus() << ")"
+ << " : (daughters: [";
for ( auto it = p.daughters_begin(), end = p.daughters_end(); it != end; ++it ) {
info() << " " << it->getObjectID().index;
}
- info() << " ]; " << endmsg;
+ 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;
@@ -31,7 +50,7 @@ bool G4PrimaryCnvTool::mutate(G4Event* anEvent) {
vertex.z*CLHEP::mm,
t);
- info() << "Geant4 Primary Vertex: ("
+ info() << "--> Creating Geant4 Primary Vertex: ("
<< vertex.x*CLHEP::mm << ","
<< vertex.y*CLHEP::mm << ","
<< vertex.z*CLHEP::mm << ")"
diff --git a/Simulation/DetSimInterface/CMakeLists.txt b/Simulation/DetSimInterface/CMakeLists.txt
index ed99eead6c48943bbe639731e2705a8e70423d31..b5fd14616a4bf34b98f3d4545d6c398954423bda 100644
--- a/Simulation/DetSimInterface/CMakeLists.txt
+++ b/Simulation/DetSimInterface/CMakeLists.txt
@@ -4,7 +4,9 @@
gaudi_add_library(DetSimInterface
SOURCES src/IDetSimSvc.cpp
+ src/CommonUserEventInfo.cc
LINK Gaudi::GaudiKernel
+ ${Geant4_LIBRARIES}
)
install(TARGETS DetSimInterface
diff --git a/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh b/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3319598a270d0d4cb38607d87b6c4fd3b2a85ad4
--- /dev/null
+++ b/Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh
@@ -0,0 +1,39 @@
+#ifndef CommonUserEventInfo_hh
+#define CommonUserEventInfo_hh
+
+/*
+ * Description:
+ * This class is a part of simulation framework to allow users to extend the G4Event.
+ *
+ * For example, when G4 converts the EDM4hep/G4 primary vertex/particle to G4 track,
+ * the relationship between the EDM4hep track and G4 track is missing.
+ * So a map is used as a bookkeeping.
+ *
+ * Author:
+ * Tao Lin <lintao AT ihep.ac.cn>
+ */
+
+#include "G4VUserEventInformation.hh"
+#include <map>
+
+class CommonUserEventInfo: public G4VUserEventInformation {
+public:
+
+ CommonUserEventInfo();
+ virtual ~CommonUserEventInfo();
+
+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
new file mode 100644
index 0000000000000000000000000000000000000000..7992c69c897967e728040b3602a472663f604c6f
--- /dev/null
+++ b/Simulation/DetSimInterface/src/CommonUserEventInfo.cc
@@ -0,0 +1,53 @@
+#include "DetSimInterface/CommonUserEventInfo.hh"
+#include <iostream>
+
+CommonUserEventInfo::CommonUserEventInfo() {
+
+}
+
+CommonUserEventInfo::~CommonUserEventInfo() {
+
+}
+
+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;
+ }
+}