add MCParticle collection to EDM4hep output

include/SimpleCaloEDM4hepEventAction.hh

@@ -42,6 +42,7 @@ namespace podio{
 
 namespace edm4hep{
   class SimCalorimeterHitCollection ;
+  class MCParticleCollection ;
 }
 
 
@@ -54,13 +55,14 @@ public:
   virtual void   EndOfEventAction(const G4Event*);
   
 private:
-  G4int    fNevent;
-  G4double fDtime;
-  G4int    fCalorimeterCollectionId;
-  G4Timer  fTimerIntern;
+  G4int    fNevent=0;
+  G4double fDtime=0.;
+  G4int    fCalorimeterCollectionId=-1;
+  G4Timer  fTimerIntern={};
   std::string fFileName = {"SimpleCaloEDM4hep.root"} ;
   podio::EventStore* fStore={} ;
   podio::ROOTWriter* fWriter={} ;
+  edm4hep::MCParticleCollection* fMCP={} ;
   edm4hep::SimCalorimeterHitCollection* fSHC={} ;
 
 };

src/SimpleCaloEDM4hepEventAction.cc

@@ -48,14 +48,14 @@ using namespace std;
 #include "podio/ROOTWriter.h"
 
 // edm4hep 
+#include "edm4hep/MCParticleCollection.h"
 #include "edm4hep/SimCalorimeterHitCollection.h"
 
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 SimpleCaloEDM4hepEventAction::SimpleCaloEDM4hepEventAction()
- : G4UserEventAction(),
-   fNevent(0),fDtime(0.0),fCalorimeterCollectionId(-1)
+  : G4UserEventAction()
 {
 
   fStore =  new podio::EventStore  ;
@@ -63,8 +63,8 @@ SimpleCaloEDM4hepEventAction::SimpleCaloEDM4hepEventAction()
 
   // create the collections to be written
 
-  // auto& mcps  = store.create<edm4hep::MCParticleCollection>("MCParticles");
-  // writer.registerForWrite("MCParticles");
+  fMCP = & fStore->create<edm4hep::MCParticleCollection>("MCParticles");
+  fWriter->registerForWrite("MCParticles");
 
   fSHC = & fStore->create<edm4hep::SimCalorimeterHitCollection>("SimCalorimeterHits");
   fWriter->registerForWrite("SimCalorimeterHits");
@@ -126,17 +126,24 @@ void SimpleCaloEDM4hepEventAction::EndOfEventAction(const G4Event *evt)
     /// Hits in sensitive Detector
     int n_hit = THC->entries();
     
-    // FIXME: create MCParticle here ...
-    // G4cout<<"  " << n_hit<< " hits are stored in SimpleCaloHitsCollection "<<G4endl;
-    // G4PrimaryVertex* pvertex=evt->GetPrimaryVertex();   
-    // ///Computing (x,y,z) of vertex of initial particles  
-    // G4ThreeVector vtx=pvertex->GetPosition();
-    // G4PrimaryParticle* pparticle=pvertex->GetPrimary();
-    // // direction of the Shower
-    // G4ThreeVector mom=pparticle->GetMomentum()/pparticle->GetMomentum().mag();
-    
+    // create MCParticle
+
+    G4PrimaryVertex* pvertex=evt->GetPrimaryVertex();
+    G4ThreeVector vtx=pvertex->GetPosition();
+    G4PrimaryParticle* pparticle=pvertex->GetPrimary();
+    G4ThreeVector mom=pparticle->GetMomentum() ;
+
+    auto mcp = fMCP->create() ;
+
+    mcp.setMass(  pparticle->GetMass() / GeV ) ;
+    mcp.setPDG( pparticle->GetPDGcode() ) ;
+    mcp.setMomentum( { float( mom[0] / GeV ), float( mom[1] / GeV ), float( mom[2] / GeV ) } ) ;
+    mcp.setVertex( { float( vtx[0] / GeV ), float( vtx[1] / GeV ), float( vtx[2] / GeV ) } ) ;
+    mcp.setGeneratorStatus( 1 ) ;
+    mcp.setCharge( pparticle->GetCharge() ) ;
+
     
-    /// For now simple convert all Hits in sensitive detector
+    /// For now simply convert all Hits in sensitive detector
     for (int i=0;i<n_hit;i++)
     {
       G4double estep = (*THC)[i]->GetEdep()/GeV;