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