Merge pull request #239 from mirguest/master

WIP: Fix the inconsistent between G4 and EDM4hep due to some particles's generator status is not 1

Analysis/CMakeLists.txt

 
 add_subdirectory(TotalInvMass)
 add_subdirectory(TrackInspect)
+add_subdirectory(DumpEvent)

Analysis/DumpEvent/CMakeLists.txt

 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

Analysis/DumpEvent/src/DumpSimHitAlg.cpp

+/*
+ * 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;
+}

Examples/options/dump_simhit.py

+#!/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
+)

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;
 }
 

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());

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,

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

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 << ")"

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

Simulation/DetSimInterface/include/DetSimInterface/CommonUserEventInfo.hh

+#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

Simulation/DetSimInterface/src/CommonUserEventInfo.cc

+#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;
+    }
+}