Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-07-15 08:15:00 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!