DDG4/edm4hep/Geant4Output2EDM4hep.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : F.Gaede, DESY
0011 //
0012 //==========================================================================
0013 #ifndef DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
0014 #define DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
0015
0016 ///  Framework include files
0017 #include <DD4hep/Detector.h>
0018 #include <DDG4/EventParameters.h>
0019 #include <DDG4/FileParameters.h>
0020 #include <DDG4/Geant4OutputAction.h>
0021 #include <DDG4/RunParameters.h>
0022
0023 /// edm4hep include files
0024 #include <edm4hep/MCParticleCollection.h>
0025 #include <edm4hep/SimTrackerHitCollection.h>
0026 #include <edm4hep/CaloHitContributionCollection.h>
0027 #include <edm4hep/SimCalorimeterHitCollection.h>
0028 #include <edm4hep/EDM4hepVersion.h>
0029 #include <edm4hep/Constants.h>
0030 #if EDM4HEP_BUILD_VERSION < EDM4HEP_VERSION(0, 99, 0)
0031   using edm4hep::CellIDEncoding;
0032 #else
0033   using edm4hep::labels::CellIDEncoding;
0034 #endif
0035 /// podio include files
0036 #include <podio/CollectionBase.h>
0037 #include <podio/podioVersion.h>
0038 #include <podio/Frame.h>
0039 #include <podio/FrameCategories.h>
0040 #if PODIO_BUILD_VERSION >= PODIO_VERSION(0, 99, 0)
0041 #include <podio/ROOTWriter.h>
0042 #else
0043 #include <podio/ROOTFrameWriter.h>
0044 namespace podio {
0045   using ROOTWriter = podio::ROOTFrameWriter;
0046 }
0047 #endif
0048
0049 /// Namespace for the AIDA detector description toolkit
0050 namespace dd4hep {
0051
0052   class ComponentCast;
0053
0054   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0055   namespace sim {
0056
0057     class  Geant4ParticleMap;
0058
0059     /// Base class to output Geant4 event data to EDM4hep
0060     /**
0061      *  \author  F.Gaede
0062      *  \version 1.0
0063      *  \ingroup DD4HEP_SIMULATION
0064      */
0065     class Geant4Output2EDM4hep : public Geant4OutputAction  {
0066     protected:
0067       using writer_t = podio::ROOTWriter;
0068       using stringmap_t = std::map< std::string, std::string >;
0069       using trackermap_t = std::map< std::string, edm4hep::SimTrackerHitCollection >;
0070       using calorimeterpair_t = std::pair< edm4hep::SimCalorimeterHitCollection, edm4hep::CaloHitContributionCollection >;
0071       using calorimetermap_t = std::map< std::string, calorimeterpair_t >;
0072       std::unique_ptr<writer_t>     m_file  { };
0073       podio::Frame                  m_frame { };
0074       edm4hep::MCParticleCollection m_particles { };
0075       trackermap_t                  m_trackerHits;
0076       calorimetermap_t              m_calorimeterHits;
0077       stringmap_t                   m_runHeader;
0078       stringmap_t                   m_eventParametersInt;
0079       stringmap_t                   m_eventParametersFloat;
0080       stringmap_t                   m_eventParametersString;
0081       stringmap_t                   m_cellIDEncodingStrings{};
0082       std::string                   m_section_name      { "events" };
0083       int                           m_runNo             { 0 };
0084       int                           m_runNumberOffset   { 0 };
0085       int                           m_eventNo           { 0 };
0086       int                           m_eventNumberOffset { 0 };
0087       bool                          m_filesByRun        { false };
0088
0089       /// Data conversion interface for MC particles to EDM4hep format
0090       void saveParticles(Geant4ParticleMap* particles);
0091       /// Store the metadata frame with e.g. the cellID encoding strings
0092       void saveFileMetaData();
0093     public:
0094       /// Standard constructor
0095       Geant4Output2EDM4hep(Geant4Context* ctxt, const std::string& nam);
0096       /// Default destructor
0097       virtual ~Geant4Output2EDM4hep();
0098       /// Callback to store the Geant4 run information
0099       virtual void beginRun(const G4Run* run);
0100       /// Callback to store the Geant4 run information
0101       virtual void endRun(const G4Run* run);
0102
0103       /// Callback to store the Geant4 run information
0104       virtual void saveRun(const G4Run* run);
0105       /// Callback to store the Geant4 event
0106       virtual void saveEvent( OutputContext<G4Event>& ctxt);
0107       /// Callback to store each Geant4 hit collection
0108       virtual void saveCollection( OutputContext<G4Event>& ctxt, G4VHitsCollection* collection);
0109       /// Commit data at end of filling procedure
0110       virtual void commit( OutputContext<G4Event>& ctxt);
0111
0112       /// begin-of-event callback - creates EDM4hep event and adds it to the event context
0113       virtual void begin(const G4Event* event);
0114     protected:
0115       /// Fill event parameters in EDM4hep event
0116       template <typename T>
0117       void saveEventParameters(const std::map<std::string, std::string >& parameters)   {
0118         for(const auto& p : parameters)   {
0119           info("Saving event parameter: %-32s = %s", p.first.c_str(), p.second.c_str());
0120           m_frame.putParameter(p.first, p.second);
0121         }
0122       }
0123     };
0124
0125     template <> void EventParameters::extractParameters(podio::Frame& frame)   {
0126       for(auto const& p: this->intParameters()) {
0127         printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str());
0128         frame.putParameter(p.first, p.second);
0129       }
0130       for(auto const& p: this->fltParameters()) {
0131         printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str());
0132         frame.putParameter(p.first, p.second);
0133       }
0134       for(auto const& p: this->strParameters()) {
0135         printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str());
0136         frame.putParameter(p.first, p.second);
0137       }
0138       // This functionality is only present in podio > 0.16.2
0139       for (auto const& p: this->dblParameters()) {
0140         printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str());
0141         frame.putParameter(p.first, p.second);
0142       }
0143     }
0144
0145     template <> void RunParameters::extractParameters(podio::Frame& frame)   {
0146       for(auto const& p: this->intParameters()) {
0147         printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str());
0148         frame.putParameter(p.first, p.second);
0149       }
0150       for(auto const& p: this->fltParameters()) {
0151         printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str());
0152         frame.putParameter(p.first, p.second);
0153       }
0154       for(auto const& p: this->strParameters()) {
0155         printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str());
0156         frame.putParameter(p.first, p.second);
0157       }
0158       // This functionality is only present in podio > 0.16.2
0159       for (auto const& p: this->dblParameters()) {
0160         printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str());
0161         frame.putParameter(p.first, p.second);
0162       }
0163     }
0164     template <> void FileParameters::extractParameters(podio::Frame& frame)   {
0165       for(auto const& p: this->intParameters()) {
0166         printout(DEBUG, "Geant4OutputEDM4hep", "Saving meta parameter: %s", p.first.c_str());
0167         frame.putParameter(p.first, p.second);
0168       }
0169       for(auto const& p: this->fltParameters()) {
0170         printout(DEBUG, "Geant4OutputEDM4hep", "Saving meta parameter: %s", p.first.c_str());
0171         frame.putParameter(p.first, p.second);
0172       }
0173       for(auto const& p: this->strParameters()) {
0174         printout(DEBUG, "Geant4OutputEDM4hep", "Saving meta parameter: %s", p.first.c_str());
0175         frame.putParameter(p.first, p.second);
0176       }
0177       // This functionality is only present in podio > 0.16.2
0178       for (auto const& p: this->dblParameters()) {
0179         printout(DEBUG, "Geant4OutputEDM4hep", "Saving meta parameter: %s", p.first.c_str());
0180         frame.putParameter(p.first, p.second);
0181       }
0182     }
0183
0184   }    // End namespace sim
0185 }      // End namespace dd4hep
0186 #endif // DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
0187
0188 //==========================================================================
0189 //  AIDA Detector description implementation
0190 //--------------------------------------------------------------------------
0191 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0192 // All rights reserved.
0193 //
0194 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0195 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0196 //
0197 // Author     : F.Gaede, DESY
0198 //
0199 //==========================================================================
0200
0201 /// Framework include files
0202 #include <DD4hep/InstanceCount.h>
0203 #include <DD4hep/VolumeManager.h>
0204
0205 #include <DDG4/Geant4HitCollection.h>
0206 #include <DDG4/Geant4DataConversion.h>
0207 #include <DDG4/Geant4SensDetAction.h>
0208 #include <DDG4/Geant4Context.h>
0209 #include <DDG4/Geant4Particle.h>
0210 #include <DDG4/Geant4Data.h>
0211
0212 ///#include <DDG4/Geant4Output2EDM4hep.h>
0213 /// Geant4 headers
0214 #include <G4Threading.hh>
0215 #include <G4AutoLock.hh>
0216 #include <G4Version.hh>
0217 #include <G4ParticleDefinition.hh>
0218 #include <G4VProcess.hh>
0219 #include <G4Event.hh>
0220 #include <G4Run.hh>
0221 /// use the Geant4 units in namespace CLHEP
0222 #include <CLHEP/Units/SystemOfUnits.h>
0223
0224 /// edm4hep include files
0225 #include <edm4hep/EventHeaderCollection.h>
0226
0227 using namespace dd4hep::sim;
0228 using namespace dd4hep;
0229
0230 namespace {
0231   G4Mutex action_mutex = G4MUTEX_INITIALIZER;
0232 }
0233
0234 #include <DDG4/Factories.h>
0235 DECLARE_GEANT4ACTION(Geant4Output2EDM4hep)
0236
0237 /// Standard constructor
0238 Geant4Output2EDM4hep::Geant4Output2EDM4hep(Geant4Context* ctxt, const std::string& nam)
0239 : Geant4OutputAction(ctxt,nam), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
0240 {
0241   declareProperty("RunHeader",             m_runHeader);
0242   declareProperty("EventParametersInt",    m_eventParametersInt);
0243   declareProperty("EventParametersFloat",  m_eventParametersFloat);
0244   declareProperty("EventParametersString", m_eventParametersString);
0245   declareProperty("RunNumberOffset",       m_runNumberOffset);
0246   declareProperty("EventNumberOffset",     m_eventNumberOffset);
0247   declareProperty("SectionName",           m_section_name);
0248   declareProperty("FilesByRun",            m_filesByRun);
0249   info("Writer is now instantiated ..." );
0250   InstanceCount::increment(this);
0251 }
0252
0253 /// Default destructor
0254 Geant4Output2EDM4hep::~Geant4Output2EDM4hep()  {
0255   G4AutoLock protection_lock(&action_mutex);
0256   m_file.reset();
0257   InstanceCount::decrement(this);
0258 }
0259
0260 // Callback to store the Geant4 run information
0261 void Geant4Output2EDM4hep::beginRun(const G4Run* run)  {
0262   G4AutoLock protection_lock(&action_mutex);
0263   std::string fname = m_output;
0264   m_runNo = run->GetRunID();
0265   if ( m_filesByRun )    {
0266     std::size_t idx = m_output.rfind(".");
0267     if ( idx != std::string::npos )   {
0268       fname = m_output.substr(0, idx) + _toString(m_runNo, ".run%08d") + m_output.substr(idx);
0269     }
0270   }
0271   if ( !fname.empty() )   {
0272     m_file = std::make_unique<podio::ROOTWriter>(fname);
0273     if ( !m_file )   {
0274       fatal("+++ Failed to open output file: %s", fname.c_str());
0275     }
0276     printout( INFO, "Geant4Output2EDM4hep" ,"Opened %s for output", fname.c_str() ) ;
0277   }
0278 }
0279
0280 /// Callback to store the Geant4 run information
0281 void Geant4Output2EDM4hep::endRun(const G4Run* run)  {
0282   saveRun(run);
0283   saveFileMetaData();
0284   if ( m_file )   {
0285     m_file->finish();
0286     m_file.reset();
0287   }
0288 }
0289
0290 void Geant4Output2EDM4hep::saveFileMetaData() {
0291   podio::Frame metaFrame{};
0292   for (const auto& [name, encodingStr] : m_cellIDEncodingStrings) {
0293     metaFrame.putParameter(podio::collMetadataParamName(name, CellIDEncoding), encodingStr);
0294   }
0295
0296   m_file->writeFrame(metaFrame, "metadata");
0297 }
0298
0299 /// Commit data at end of filling procedure
0300 void Geant4Output2EDM4hep::commit( OutputContext<G4Event>& /* ctxt */)   {
0301   if ( m_file )   {
0302     G4AutoLock protection_lock(&action_mutex);
0303     m_frame.put( std::move(m_particles), "MCParticles");
0304     for (auto it = m_trackerHits.begin(); it != m_trackerHits.end(); ++it)   {
0305       m_frame.put( std::move(it->second), it->first);
0306     }
0307     for (auto& [colName, calorimeterHits] : m_calorimeterHits) {
0308       m_frame.put( std::move(calorimeterHits.first), colName);
0309       m_frame.put( std::move(calorimeterHits.second), colName + "Contributions");
0310     }
0311     m_file->writeFrame(m_frame, m_section_name);
0312     m_particles.clear();
0313     m_trackerHits.clear();
0314     m_calorimeterHits.clear();
0315     m_frame = {};
0316     return;
0317   }
0318   except("+++ Failed to write output file. [Stream is not open]");
0319 }
0320
0321 /// Callback to store the Geant4 run information
0322 void Geant4Output2EDM4hep::saveRun(const G4Run* run)   {
0323   G4AutoLock protection_lock(&action_mutex);
0324   // --- write an edm4hep::RunHeader ---------
0325   // Runs are just Frames with different contents in EDM4hep / podio. We simply
0326   // store everything as parameters for now
0327   podio::Frame runHeader {};
0328   for (const auto& [key, value] : m_runHeader)
0329     runHeader.putParameter(key, value);
0330
0331   m_runNo = m_runNumberOffset > 0 ? m_runNumberOffset + run->GetRunID() : run->GetRunID();
0332   runHeader.putParameter("runNumber", m_runNo);
0333   runHeader.putParameter("GEANT4Version", G4Version);
0334   runHeader.putParameter("DD4hepVersion", versionString());
0335   runHeader.putParameter("detectorName", context()->detectorDescription().header().name());
0336   {
0337     RunParameters* parameters = context()->run().extension<RunParameters>(false);
0338     if ( parameters ) {
0339       parameters->extractParameters(runHeader);
0340     }
0341     m_file->writeFrame(runHeader, "runs");
0342   }
0343   {
0344     podio::Frame metaFrame {};
0345     FileParameters* parameters = context()->run().extension<FileParameters>(false);
0346     if ( parameters ) {
0347       parameters->extractParameters(metaFrame);
0348     }
0349     m_file->writeFrame(metaFrame, "meta");
0350   }
0351 }
0352
0353 void Geant4Output2EDM4hep::begin(const G4Event* event)  {
0354   /// Create event frame object
0355   m_eventNo = event->GetEventID();
0356   m_frame = {};
0357   m_particles = {};
0358   m_trackerHits.clear();
0359   m_calorimeterHits.clear();
0360 }
0361
0362 /// Data conversion interface for MC particles to EDM4hep format
0363 void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles)    {
0364   typedef detail::ReferenceBitMask<const int> PropertyMask;
0365   typedef Geant4ParticleMap::ParticleMap ParticleMap;
0366   const ParticleMap& pm = particles->particleMap;
0367
0368   m_particles.clear();
0369   if ( pm.size() > 0 )  {
0370     size_t cnt = 0;
0371     // Mapping of ids in the ParticleMap to indices in the MCParticle collection
0372     std::map<int,int> p_ids;
0373     std::vector<const Geant4Particle*> p_part;
0374     p_part.reserve(pm.size());
0375     // First create the particles
0376     for (const auto& iParticle : pm) {
0377       int id = iParticle.first;
0378       const Geant4ParticleHandle p = iParticle.second;
0379       PropertyMask mask(p->status);
0380       //      std::cout << " ********** mcp status : 0x" << std::hex << p->status << ", mask.isSet(G4PARTICLE_GEN_STABLE) x" << std::dec << mask.isSet(G4PARTICLE_GEN_STABLE)  <<std::endl ;
0381       const G4ParticleDefinition* def = p.definition();
0382       auto mcp = m_particles.create();
0383       mcp.setPDG(p->pdgID);
0384       // Because EDM4hep is switching between vector3f[loat] and vector3d[ouble]
0385       using MT = decltype(std::declval<edm4hep::MCParticle>().getMomentum().x);
0386       mcp.setMomentum( {MT(p->psx/CLHEP::GeV),MT(p->psy/CLHEP::GeV),MT(p->psz/CLHEP::GeV)} );
0387       mcp.setMomentumAtEndpoint( {MT(p->pex/CLHEP::GeV),MT(p->pey/CLHEP::GeV),MT(p->pez/CLHEP::GeV)} );
0388
0389       double vs_fa[3] = { p->vsx/CLHEP::mm, p->vsy/CLHEP::mm, p->vsz/CLHEP::mm } ;
0390       mcp.setVertex( vs_fa );
0391
0392       double ve_fa[3] = { p->vex/CLHEP::mm, p->vey/CLHEP::mm, p->vez/CLHEP::mm } ;
0393       mcp.setEndpoint( ve_fa );
0394
0395       mcp.setTime(p->time/CLHEP::ns);
0396       mcp.setMass(p->mass/CLHEP::GeV);
0397       mcp.setCharge(def ? def->GetPDGCharge() : 0); // Charge(e+) = 1 !
0398
0399       // Set generator status
0400       mcp.setGeneratorStatus(0);
0401       if( p->genStatus ) {
0402         mcp.setGeneratorStatus( p->genStatus ) ;
0403       } else {
0404         if ( mask.isSet(G4PARTICLE_GEN_STABLE) )             mcp.setGeneratorStatus(1);
0405         else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) )       mcp.setGeneratorStatus(2);
0406         else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) mcp.setGeneratorStatus(3);
0407         else if ( mask.isSet(G4PARTICLE_GEN_BEAM) )          mcp.setGeneratorStatus(4);
0408         else if ( mask.isSet(G4PARTICLE_GEN_OTHER) )         mcp.setGeneratorStatus(9);
0409       }
0410
0411       // Set simulation status
0412       mcp.setCreatedInSimulation(         mask.isSet(G4PARTICLE_SIM_CREATED) );
0413       mcp.setBackscatter(                 mask.isSet(G4PARTICLE_SIM_BACKSCATTER) );
0414       mcp.setVertexIsNotEndpointOfParent( mask.isSet(G4PARTICLE_SIM_PARENT_RADIATED) );
0415       mcp.setDecayedInTracker(            mask.isSet(G4PARTICLE_SIM_DECAY_TRACKER) );
0416       mcp.setDecayedInCalorimeter(        mask.isSet(G4PARTICLE_SIM_DECAY_CALO) );
0417       mcp.setHasLeftDetector(             mask.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) );
0418       mcp.setStopped(                     mask.isSet(G4PARTICLE_SIM_STOPPED) );
0419       mcp.setOverlay(                     false );
0420
0421       //fg: if simstatus !=0 we have to set the generator status to 0:
0422       if( mcp.isCreatedInSimulation() )
0423         mcp.setGeneratorStatus( 0 )  ;
0424
0425       mcp.setSpin(p->spin);
0426
0427       p_ids[id] = cnt++;
0428       p_part.push_back(p);
0429     }
0430
0431     // Now establish parent-daughter relationships
0432     for(size_t i=0; i < p_ids.size(); ++i)   {
0433       const Geant4Particle* p = p_part[i];
0434       auto q = m_particles[i];
0435
0436       for (const auto& idau : p->daughters) {
0437         const auto k = p_ids.find(idau);
0438         if (k == p_ids.end()) {
0439           fatal("+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
0440           continue;
0441         }
0442         int iqdau = (*k).second;
0443         auto qdau = m_particles[iqdau];
0444         q.addToDaughters(qdau);
0445       }
0446
0447       for (const auto& ipar : p->parents) {
0448         if (ipar >= 0) { // A parent ID of -1 means NO parent, because a base of 0 is perfectly legal
0449           const auto k = p_ids.find(ipar);
0450           if (k == p_ids.end()) {
0451             fatal("+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
0452             continue;
0453           }
0454           int iqpar = (*k).second;
0455           auto qpar = m_particles[iqpar];
0456           q.addToParents(qpar);
0457         }
0458       }
0459     }
0460   }
0461 }
0462
0463 /// Callback to store the Geant4 event
0464 void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt)  {
0465   EventParameters* parameters = context()->event().extension<EventParameters>(false);
0466   int runNumber(0), eventNumber(0);
0467   const int eventNumberOffset(m_eventNumberOffset > 0 ? m_eventNumberOffset : 0);
0468   const int runNumberOffset(m_runNumberOffset > 0 ? m_runNumberOffset : 0);
0469   double eventWeight{0};
0470   // Get event number, run number and parameters from extension ...
0471   if ( parameters ) {
0472     runNumber = parameters->runNumber() + runNumberOffset;
0473     eventNumber = parameters->eventNumber() + eventNumberOffset;
0474     parameters->extractParameters(m_frame);
0475 #if PODIO_BUILD_VERSION > PODIO_VERSION(0, 99, 0)
0476     eventWeight = m_frame.getParameter<double>("EventWeights").value_or(0.0);
0477 #else
0478     eventWeight = m_frame.getParameter<double>("EventWeights");
0479 #endif
0480   } else { // ... or from DD4hep framework
0481     runNumber = m_runNo + runNumberOffset;
0482     eventNumber = ctxt.context->GetEventID() + eventNumberOffset;
0483   }
0484   printout(INFO,"Geant4Output2EDM4hep","+++ Saving EDM4hep event %d run %d.", eventNumber, runNumber);
0485
0486   // this does not compile as create() is we only get a const ref - need to review PODIO EventStore API
0487   edm4hep::EventHeaderCollection header_collection;
0488
0489   auto header = header_collection.create();
0490   header.setRunNumber(runNumber);
0491   header.setEventNumber(eventNumber);
0492   header.setWeight(eventWeight);
0493   //not implemented in EDM4hep ?  header.setDetectorName(context()->detectorDescription().header().name());
0494   header.setTimeStamp(std::time(nullptr));
0495
0496   // extract event header, in case we come from edm4hep input
0497   auto* meh = context()->event().extension<edm4hep::MutableEventHeader>(false);
0498   if(meh) {
0499     header.setTimeStamp(meh->getTimeStamp());
0500 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
0501     for (auto const& weight: meh->getWeights()) {
0502       header.addToWeights(weight);
0503     }
0504 #endif
0505   }
0506
0507   m_frame.put(std::move(header_collection), "EventHeader");
0508   saveEventParameters<int>(m_eventParametersInt);
0509   saveEventParameters<float>(m_eventParametersFloat);
0510   saveEventParameters<std::string>(m_eventParametersString);
0511
0512   Geant4ParticleMap* part_map = context()->event().extension<Geant4ParticleMap>(false);
0513   if ( part_map )   {
0514     print("+++ Saving %d EDM4hep particles....",int(part_map->particleMap.size()));
0515     if ( part_map->particleMap.size() > 0 )  {
0516       saveParticles(part_map);
0517     }
0518   }
0519 }
0520
0521 /**
0522  * Helper struct that can be used together with map::try_emplace to construct
0523  * the encoding only once per collection (name).
0524  */
0525 struct LazyEncodingExtraction {
0526   /// Constructor that does effectively nothing. This will be called in every
0527   /// try_emplace call
0528   LazyEncodingExtraction(Geant4HitCollection* coll) : m_coll(coll) {}
0529   /// Defer the real work to the implicit conversion to std::string that will
0530   /// only be called if the value is actually emplaced into the map
0531   operator std::string() const {
0532     const auto* sd = m_coll->sensitive();
0533     return dd4hep::sim::Geant4ConversionHelper::encoding(sd->sensitiveDetector());
0534   }
0535 private:
0536   Geant4HitCollection* m_coll{nullptr};
0537 };
0538
0539
0540 /// Callback to store each Geant4 hit collection
0541 void Geant4Output2EDM4hep::saveCollection(OutputContext<G4Event>& /*ctxt*/, G4VHitsCollection* collection)  {
0542   Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(collection);
0543   std::string colName = collection->GetName();
0544   if( coll == nullptr ){
0545     error(" no Geant4HitCollection:  %s ", colName.c_str());
0546     return ;
0547   }
0548   size_t nhits = collection->GetSize();
0549   Geant4ParticleMap* pm = context()->event().extension<Geant4ParticleMap>(false);
0550   debug("+++ Saving EDM4hep collection %s with %d entries.", colName.c_str(), int(nhits));
0551
0552   // Using try_emplace here to only fill this the first time we come across
0553   m_cellIDEncodingStrings.try_emplace(colName, LazyEncodingExtraction{coll});
0554
0555   //-------------------------------------------------------------------
0556   if( typeid( Geant4Tracker::Hit ) == coll->type().type()  ){
0557     // Create the hit container even if there are no entries!
0558     auto& hits = m_trackerHits[colName];
0559     for(unsigned i=0 ; i < nhits ; ++i){
0560       auto sth = hits->create();
0561       const Geant4Tracker::Hit* hit = coll->hit(i);
0562       const Geant4Tracker::Hit::Contribution& t = hit->truth;
0563       int   trackID   = pm->particleID(t.trackID);
0564       auto  mcp       = m_particles.at(trackID);
0565       const auto& mom = hit->momentum;
0566       const auto& pos = hit->position;
0567       edm4hep::Vector3f();
0568       sth.setCellID( hit->cellID ) ;
0569       sth.setEDep(hit->energyDeposit/CLHEP::GeV);
0570       sth.setPathLength(hit->length/CLHEP::mm);
0571       sth.setTime(hit->truth.time/CLHEP::ns);
0572 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 10, 99)
0573       sth.setParticle(mcp);
0574 #else
0575       sth.setMCParticle(mcp);
0576 #endif
0577       sth.setPosition( {pos.x()/CLHEP::mm, pos.y()/CLHEP::mm, pos.z()/CLHEP::mm} );
0578       sth.setMomentum( {float(mom.x()/CLHEP::GeV),float(mom.y()/CLHEP::GeV),float(mom.z()/CLHEP::GeV)} );
0579       auto particleIt = pm->particles().find(trackID);
0580       if( ( particleIt != pm->particles().end()) ){
0581         // if the original track ID of the particle is not the same as the
0582         // original track ID of the hit it was produced by an MCParticle that
0583         // is no longer stored
0584         sth.setProducedBySecondary( (particleIt->second->originalG4ID != t.trackID) );
0585       }
0586     }
0587     //-------------------------------------------------------------------
0588   }
0589   else if( typeid( Geant4Calorimeter::Hit ) == coll->type().type() ){
0590     Geant4Sensitive* sd = coll->sensitive();
0591     int hit_creation_mode = sd->hitCreationMode();
0592     // Create the hit container even if there are no entries!
0593     auto& hits = m_calorimeterHits[colName];
0594     for(unsigned i=0 ; i < nhits ; ++i){
0595       auto sch = hits.first->create();
0596       const Geant4Calorimeter::Hit* hit = coll->hit(i);
0597       const auto& pos = hit->position;
0598       sch.setCellID( hit->cellID );
0599       sch.setPosition({float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm)});
0600       sch.setEnergy( hit->energyDeposit/CLHEP::GeV );
0601
0602
0603       // now add the individual step contributions
0604       for(auto ci=hit->truth.begin(); ci != hit->truth.end(); ++ci){
0605
0606         auto sCaloHitCont = hits.second->create();
0607         sch.addToContributions( sCaloHitCont );
0608
0609         const Geant4HitData::Contribution& c = *ci;
0610         int trackID = pm->particleID(c.trackID);
0611         auto mcp = m_particles.at(trackID);
0612         sCaloHitCont.setEnergy( c.deposit/CLHEP::GeV );
0613         sCaloHitCont.setTime( c.time/CLHEP::ns );
0614         sCaloHitCont.setParticle( mcp );
0615
0616         if ( hit_creation_mode == Geant4Sensitive::DETAILED_MODE )     {
0617           edm4hep::Vector3f p(c.x/CLHEP::mm, c.y/CLHEP::mm, c.z/CLHEP::mm);
0618           sCaloHitCont.setPDG( c.pdgID );
0619           sCaloHitCont.setStepPosition( p );
0620         }
0621       }
0622     }
0623     //-------------------------------------------------------------------
0624   } else {
0625     error("+++ unknown type in Geant4HitCollection %s ", coll->type().type().name());
0626   }
0627 }