Warning, /eic-spack/packages/dd4hep/revert-Geant4Output2EDM4hep-dd4hep-1-25-1-to-1-23.patch is written in an unsupported language. File is not indexed.
0001 diff --git a/DDG4/edm4hep/Geant4Output2EDM4hep.cpp b/DDG4/edm4hep/Geant4Output2EDM4hep.cpp
0002 --- a/DDG4/edm4hep/Geant4Output2EDM4hep.cpp
0003 +++ b/DDG4/edm4hep/Geant4Output2EDM4hep.cpp
0004 @@ -10,19 +10,40 @@
0005 // Author : F.Gaede, DESY
0006 //
0007 //==========================================================================
0008 +
0009 #ifndef DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
0010 #define DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
0011
0012 -/// Framework include files
0013 -#include <DD4hep/Detector.h>
0014 -#include <DDG4/EventParameters.h>
0015 -#include <DDG4/Geant4OutputAction.h>
0016 -
0017 -/// edm4hep include files
0018 -#include <edm4hep/MCParticleCollection.h>
0019 -/// podio include files
0020 -#include <podio/Frame.h>
0021 -#include <podio/ROOTFrameWriter.h>
0022 +// Framework include files
0023 +#include "DD4hep/Detector.h"
0024 +#include "DD4hep/VolumeManager.h"
0025 +#include "DDG4/Geant4HitCollection.h"
0026 +#include "DDG4/Geant4OutputAction.h"
0027 +#include "DDG4/Geant4SensDetAction.h"
0028 +#include "DDG4/Geant4DataConversion.h"
0029 +#include "DDG4/EventParameters.h"
0030 +
0031 +// Geant4 headers
0032 +#include "G4Threading.hh"
0033 +#include "G4AutoLock.hh"
0034 +#include "G4Version.hh"
0035 +
0036 +// use the Geant4 units in namespace CLHEP
0037 +#include "CLHEP/Units/SystemOfUnits.h"
0038 +
0039 +// edm4hep include files
0040 +#include "edm4hep/EventHeaderCollection.h"
0041 +#include "edm4hep/MCParticleCollection.h"
0042 +#include "edm4hep/SimTrackerHitCollection.h"
0043 +#include "edm4hep/CaloHitContributionCollection.h"
0044 +#include "edm4hep/SimCalorimeterHitCollection.h"
0045 +#include "podio/EventStore.h"
0046 +#include "podio/ROOTWriter.h"
0047 +
0048 +#include <typeinfo>
0049 +#include <iostream>
0050 +#include <ctime>
0051 +#include <unordered_map>
0052 #include <podio/podioVersion.h>
0053
0054 /// Namespace for the AIDA detector description toolkit
0055 @@ -33,6 +54,20 @@ namespace dd4hep {
0056 /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0057 namespace sim {
0058
0059 + template <class T=podio::EventStore> void EventParameters::extractParameters(T& event){
0060 + auto& lcparameters = event.getEventMetaData();
0061 +
0062 + for(auto const& ival: this->intParameters()) {
0063 + lcparameters.setValue(ival.first, ival.second);
0064 + }
0065 + for(auto const& ival: this->fltParameters()) {
0066 + lcparameters.setValue(ival.first, ival.second);
0067 + }
0068 + for(auto const& ival: this->strParameters()) {
0069 + lcparameters.setValue(ival.first, ival.second);
0070 + }
0071 + }
0072 +
0073 class Geant4ParticleMap;
0074
0075 /// Base class to output Geant4 event data to EDM4hep
0076 @@ -43,22 +78,19 @@ namespace dd4hep {
0077 */
0078 class Geant4Output2EDM4hep : public Geant4OutputAction {
0079 protected:
0080 - using writer_t = podio::ROOTFrameWriter;
0081 - using stringmap_t = std::map< std::string, std::string >;
0082 - std::unique_ptr<writer_t> m_file { };
0083 - podio::Frame m_frame { };
0084 - edm4hep::MCParticleCollection m_particles { };
0085 - stringmap_t m_runHeader;
0086 - stringmap_t m_eventParametersInt;
0087 - stringmap_t m_eventParametersFloat;
0088 - stringmap_t m_eventParametersString;
0089 - std::string m_section_name { "events" };
0090 - int m_runNo { 0 };
0091 - int m_runNumberOffset { 0 };
0092 - int m_eventNo { 0 };
0093 - int m_eventNumberOffset { 0 };
0094 - bool m_filesByRun { false };
0095 -
0096 + podio::EventStore* m_store;
0097 + podio::ROOTWriter* m_file;
0098 + int m_runNo;
0099 + int m_runNumberOffset;
0100 + int m_eventNumberOffset;
0101 + std::map< std::string, std::string > m_runHeader;
0102 + std::map< std::string, std::string > m_eventParametersInt;
0103 + std::map< std::string, std::string > m_eventParametersFloat;
0104 + std::map< std::string, std::string > m_eventParametersString;
0105 + bool m_FirstEvent = true ;
0106 +
0107 + std::unordered_map<std::string, podio::CollectionBase*> m_collections;
0108 +
0109 /// create the podio collections for the particles and hits
0110 void createCollections(OutputContext<G4Event>& ctxt) ;
0111 /// Data conversion interface for MC particles to EDM4hep format
0112 @@ -87,35 +119,33 @@ namespace dd4hep {
0113 protected:
0114 /// Fill event parameters in EDM4hep event
0115 template <typename T>
0116 - void saveEventParameters(const std::map<std::string, std::string >& parameters) {
0117 - for(const auto& p : parameters) {
0118 - info("Saving event parameter: %-32s = %s", p.first.c_str(), p.second.c_str());
0119 - m_frame.putParameter(p.first, p.second);
0120 - }
0121 - }
0122 + void saveEventParameters(const std::map<std::string, std::string >& parameters);
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 + /// Fill event parameters in EDM4hep event
0138 + template <typename T>
0139 + inline void Geant4Output2EDM4hep::saveEventParameters(const std::map<std::string, std::string >& parameters) {
0140 + for(std::map<std::string, std::string >::const_iterator iter = parameters.begin(), endIter = parameters.end() ; iter != endIter ; ++iter) {
0141 + T parameter;
0142 + std::istringstream iss(iter->second);
0143 + if ( (iss >> parameter).fail() ) {
0144 + printout(FATAL,"saveEventParameters","+++ Event parameter %s: FAILED to convert to type :%s",iter->first.c_str(),typeid(T).name());
0145 + continue;
0146 + }
0147 + auto& evtMD = m_store->getEventMetaData();
0148 + evtMD.setValue(iter->first,parameter);
0149 }
0150 -#if podio_VERSION_MAJOR > 0 || podio_VERSION_MINOR > 16 || podio_VERSION_PATCH > 2
0151 - // This functionality is only present in podio > 0.16.2
0152 - for (auto const& p: this->dblParameters()) {
0153 - printout(DEBUG, "Geant4OutputEDM4hep", "Saving event parameter: %s", p.first.c_str());
0154 - frame.putParameter(p.first, p.second);
0155 + }
0156 +
0157 + /// Fill event parameters in EDM4hep event - std::string specialization
0158 + template <>
0159 + inline void Geant4Output2EDM4hep::saveEventParameters<std::string>(const std::map<std::string, std::string >& parameters) {
0160 + for(std::map<std::string, std::string >::const_iterator iter = parameters.begin(), endIter = parameters.end() ; iter != endIter ; ++iter) {
0161 + auto& evtMD = m_store->getEventMetaData();
0162 + evtMD.setValue(iter->first,iter->second);
0163 }
0164 -#endif
0165 }
0166 +
0167 } // End namespace sim
0168 } // End namespace dd4hep
0169 #endif // DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
0170 @@ -133,38 +163,26 @@ namespace dd4hep {
0171 //
0172 //==========================================================================
0173
0174 -/// Framework include files
0175 -#include <DD4hep/InstanceCount.h>
0176 -#include <DD4hep/VolumeManager.h>
0177 -
0178 -#include <DDG4/Geant4HitCollection.h>
0179 -#include <DDG4/Geant4DataConversion.h>
0180 -#include <DDG4/Geant4SensDetAction.h>
0181 -#include <DDG4/Geant4Context.h>
0182 -#include <DDG4/Geant4Particle.h>
0183 -#include <DDG4/Geant4Data.h>
0184 -
0185 -///#include <DDG4/Geant4Output2EDM4hep.h>
0186 -/// Geant4 headers
0187 -#include <G4Threading.hh>
0188 -#include <G4AutoLock.hh>
0189 -#include <G4Version.hh>
0190 -#include <G4ParticleDefinition.hh>
0191 -#include <G4VProcess.hh>
0192 -#include <G4Event.hh>
0193 -#include <G4Run.hh>
0194 -/// use the Geant4 units in namespace CLHEP
0195 -#include <CLHEP/Units/SystemOfUnits.h>
0196 -
0197 -/// edm4hep include files
0198 -#include <edm4hep/EventHeaderCollection.h>
0199 -#include <edm4hep/SimTrackerHitCollection.h>
0200 -#include <edm4hep/CaloHitContributionCollection.h>
0201 -#include <edm4hep/SimCalorimeterHitCollection.h>
0202 +// Framework include files
0203 +#include "DD4hep/InstanceCount.h"
0204 +#include "DD4hep/Detector.h"
0205 +#include "DDG4/Geant4HitCollection.h"
0206 +#include "DDG4/Geant4DataConversion.h"
0207 +#include "DDG4/Geant4Context.h"
0208 +#include "DDG4/Geant4Particle.h"
0209 +#include "DDG4/Geant4Data.h"
0210 +#include "DDG4/Geant4Action.h"
0211 +
0212 +//#include "DDG4/Geant4Output2EDM4hep.h"
0213 +#include "G4ParticleDefinition.hh"
0214 +#include "G4VProcess.hh"
0215 +#include "G4Event.hh"
0216 +#include "G4Run.hh"
0217 +
0218
0219 using namespace dd4hep::sim;
0220 using namespace dd4hep;
0221 -
0222 +using namespace std;
0223 namespace {
0224 G4Mutex action_mutex=G4MUTEX_INITIALIZER;
0225 }
0226 @@ -173,65 +191,55 @@ namespace {
0227 DECLARE_GEANT4ACTION(Geant4Output2EDM4hep)
0228
0229 /// Standard constructor
0230 -Geant4Output2EDM4hep::Geant4Output2EDM4hep(Geant4Context* ctxt, const std::string& nam)
0231 -: Geant4OutputAction(ctxt,nam), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
0232 +Geant4Output2EDM4hep::Geant4Output2EDM4hep(Geant4Context* ctxt, const string& nam)
0233 +: Geant4OutputAction(ctxt,nam), m_store(0), m_file(0), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
0234 {
0235 declareProperty("RunHeader", m_runHeader);
0236 declareProperty("EventParametersInt", m_eventParametersInt);
0237 declareProperty("EventParametersFloat", m_eventParametersFloat);
0238 declareProperty("EventParametersString", m_eventParametersString);
0239 - declareProperty("RunNumberOffset", m_runNumberOffset);
0240 - declareProperty("EventNumberOffset", m_eventNumberOffset);
0241 - declareProperty("SectionName", m_section_name);
0242 - declareProperty("FilesByRun", m_filesByRun);
0243 - info("Writer is now instantiated ..." );
0244 + declareProperty("RunNumberOffset", m_runNumberOffset);
0245 + declareProperty("EventNumberOffset", m_eventNumberOffset);
0246 + printout( INFO, "Geant4Output2EDM4hep" ," instantiated ..." ) ;
0247 InstanceCount::increment(this);
0248 }
0249
0250 /// Default destructor
0251 Geant4Output2EDM4hep::~Geant4Output2EDM4hep() {
0252 G4AutoLock protection_lock(&action_mutex);
0253 - m_file.reset();
0254 + if ( m_file ) {
0255 + m_file->finish();
0256 + detail::deletePtr(m_file);
0257 + }
0258 + if (nullptr != m_store) {
0259 + delete m_store;
0260 + }
0261 InstanceCount::decrement(this);
0262 }
0263
0264 // Callback to store the Geant4 run information
0265 void Geant4Output2EDM4hep::beginRun(const G4Run* run) {
0266 G4AutoLock protection_lock(&action_mutex);
0267 - std::string fname = m_output;
0268 - m_runNo = run->GetRunID();
0269 - if ( m_filesByRun ) {
0270 - std::size_t idx = m_output.rfind(".");
0271 - if ( idx != std::string::npos ) {
0272 - fname = m_output.substr(0, idx) + _toString(m_runNo, ".run%08d") + m_output.substr(idx);
0273 - }
0274 - }
0275 - if ( !m_file && !fname.empty() ) {
0276 - m_file = std::make_unique<podio::ROOTFrameWriter>(fname);
0277 - if ( !m_file ) {
0278 - fatal("+++ Failed to open output file: %s", fname.c_str());
0279 - }
0280 - printout( INFO, "Geant4Output2EDM4hep" ,"Opened %s for output", fname.c_str() ) ;
0281 + if ( 0 == m_file && !m_output.empty() ) {
0282 + m_store = new podio::EventStore ;
0283 + m_file = new podio::ROOTWriter(m_output, m_store);
0284 +
0285 + printout( INFO, "Geant4Output2EDM4hep" ," opened %s for output", m_output.c_str() ) ;
0286 }
0287 +
0288 + saveRun(run);
0289 }
0290
0291 /// Callback to store the Geant4 run information
0292 -void Geant4Output2EDM4hep::endRun(const G4Run* run) {
0293 - saveRun(run);
0294 - if ( m_file ) {
0295 - m_file->finish();
0296 - m_file.reset();
0297 - }
0298 +void Geant4Output2EDM4hep::endRun(const G4Run* /*run*/) {
0299 }
0300
0301 /// Commit data at end of filling procedure
0302 void Geant4Output2EDM4hep::commit( OutputContext<G4Event>& /* ctxt */) {
0303 if ( m_file ) {
0304 G4AutoLock protection_lock(&action_mutex);
0305 - m_frame.put( std::move(m_particles), "MCParticles");
0306 - m_file->writeFrame(m_frame, m_section_name);
0307 - m_particles.clear();
0308 - m_frame = {};
0309 + m_file->writeEvent();
0310 + m_store->clearCollections();
0311 return;
0312 }
0313 except("+++ Failed to write output file. [Stream is not open]");
0314 @@ -239,31 +247,12 @@ void Geant4Output2EDM4hep::commit( OutputContext<G4Event>& /* ctxt */) {
0315
0316 /// Callback to store the Geant4 run information
0317 void Geant4Output2EDM4hep::saveRun(const G4Run* run) {
0318 - G4AutoLock protection_lock(&action_mutex);
0319 - // --- write an edm4hep::RunHeader ---------
0320 - // Runs are just Frames with different contents in EDM4hep / podio. We simply
0321 - // store everything as parameters for now
0322 - auto runHeader = podio::Frame();
0323 -
0324 - for (const auto& [key, value] : m_runHeader) {
0325 - runHeader.putParameter(key, value);
0326 - }
0327 + //G4AutoLock protection_lock(&action_mutex);
0328
0329 - m_runNo = m_runNumberOffset > 0 ? m_runNumberOffset + run->GetRunID() : run->GetRunID();
0330 - runHeader.putParameter("runNumber", m_runNo);
0331 - runHeader.putParameter("GEANT4Version", G4Version);
0332 - runHeader.putParameter("DD4hepVersion", versionString());
0333 - runHeader.putParameter("detectorName", context()->detectorDescription().header().name());
0334 -
0335 - m_file->writeFrame(runHeader, "runs");
0336 + warning("saveRun: RunHeader not implemented in EDM4hep, nothing written for run: %d", run->GetRunID());
0337 }
0338
0339 -void Geant4Output2EDM4hep::begin(const G4Event* event) {
0340 - /// Create event frame object
0341 - m_eventNo = event->GetEventID();
0342 - m_particles.clear();
0343 - m_particles = {};
0344 - m_frame = {};
0345 +void Geant4Output2EDM4hep::begin(const G4Event* /* event */) {
0346 }
0347
0348 /// Data conversion interface for MC particles to EDM4hep format
0349 @@ -271,13 +260,13 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
0350 typedef detail::ReferenceBitMask<const int> PropertyMask;
0351 typedef Geant4ParticleMap::ParticleMap ParticleMap;
0352 const ParticleMap& pm = particles->particleMap;
0353 + auto mcpc = static_cast<edm4hep::MCParticleCollection*>(m_collections["MCParticles"]);
0354
0355 - m_particles.clear();
0356 if ( pm.size() > 0 ) {
0357 size_t cnt = 0;
0358 // Mapping of ids in the ParticleMap to indices in the MCParticle collection
0359 - std::map<int,int> p_ids;
0360 - std::vector<const Geant4Particle*> p_part;
0361 + map<int,int> p_ids;
0362 + vector<const Geant4Particle*> p_part;
0363 p_part.reserve(pm.size());
0364 // First create the particles
0365 for (const auto& iParticle : pm) {
0366 @@ -286,7 +275,7 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
0367 PropertyMask mask(p->status);
0368 // std::cout << " ********** mcp status : 0x" << std::hex << p->status << ", mask.isSet(G4PARTICLE_GEN_STABLE) x" << std::dec << mask.isSet(G4PARTICLE_GEN_STABLE) <<std::endl ;
0369 const G4ParticleDefinition* def = p.definition();
0370 - auto mcp = m_particles.create();
0371 + auto mcp = mcpc->create();
0372 mcp.setPDG(p->pdgID);
0373
0374 float ps_fa[3] = {float(p->psx/CLHEP::GeV),float(p->psy/CLHEP::GeV),float(p->psz/CLHEP::GeV)};
0375 @@ -308,13 +297,13 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
0376 // Set generator status
0377 mcp.setGeneratorStatus(0);
0378 if( p->genStatus ) {
0379 - mcp.setGeneratorStatus( p->genStatus ) ;
0380 + mcp.setGeneratorStatus( p->genStatus ) ;
0381 } else {
0382 - if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) mcp.setGeneratorStatus(1);
0383 - else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) ) mcp.setGeneratorStatus(2);
0384 - else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) mcp.setGeneratorStatus(3);
0385 - else if ( mask.isSet(G4PARTICLE_GEN_BEAM) ) mcp.setGeneratorStatus(4);
0386 - else if ( mask.isSet(G4PARTICLE_GEN_OTHER) ) mcp.setGeneratorStatus(9);
0387 + if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) mcp.setGeneratorStatus(1);
0388 + else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) ) mcp.setGeneratorStatus(2);
0389 + else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) mcp.setGeneratorStatus(3);
0390 + else if ( mask.isSet(G4PARTICLE_GEN_BEAM) ) mcp.setGeneratorStatus(4);
0391 + else if ( mask.isSet(G4PARTICLE_GEN_OTHER) ) mcp.setGeneratorStatus(9);
0392 }
0393
0394 // Set simulation status
0395 @@ -329,7 +318,7 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
0396
0397 //fg: if simstatus !=0 we have to set the generator status to 0:
0398 if( mcp.isCreatedInSimulation() )
0399 - mcp.setGeneratorStatus( 0 ) ;
0400 + mcp.setGeneratorStatus( 0 ) ;
0401
0402 mcp.setSpin(p->spin);
0403 mcp.setColorFlow(p->colorFlow);
0404 @@ -341,30 +330,30 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
0405 // Now establish parent-daughter relationships
0406 for(size_t i=0; i < p_ids.size(); ++i) {
0407 const Geant4Particle* p = p_part[i];
0408 - auto q = m_particles[i];
0409 + auto q = (*mcpc)[i];
0410
0411 for (const auto& idau : p->daughters) {
0412 - const auto k = p_ids.find(idau);
0413 - if (k == p_ids.end()) {
0414 - fatal("+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
0415 - continue;
0416 - }
0417 - int iqdau = (*k).second;
0418 - auto qdau = m_particles[iqdau];
0419 - q.addToDaughters(qdau);
0420 + const auto k = p_ids.find(idau);
0421 + if (k == p_ids.end()) {
0422 + printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
0423 + continue;
0424 + }
0425 + int iqdau = (*k).second;
0426 + auto qdau = (*mcpc)[iqdau];
0427 + q.addToDaughters(qdau);
0428 }
0429
0430 for (const auto& ipar : p->parents) {
0431 - if (ipar >= 0) { // A parent ID of -1 means NO parent, because a base of 0 is perfectly legal
0432 - const auto k = p_ids.find(ipar);
0433 - if (k == p_ids.end()) {
0434 - fatal("+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
0435 - continue;
0436 - }
0437 - int iqpar = (*k).second;
0438 - auto qpar = m_particles[iqpar];
0439 - q.addToParents(qpar);
0440 - }
0441 + if (ipar >= 0) { // A parent ID of -1 means NO parent, because a base of 0 is perfectly legal
0442 + const auto k = p_ids.find(ipar);
0443 + if (k == p_ids.end()) {
0444 + printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
0445 + continue;
0446 + }
0447 + int iqpar = (*k).second;
0448 + auto qpar = (*mcpc)[iqpar];
0449 + q.addToParents(qpar);
0450 + }
0451 }
0452 }
0453 }
0454 @@ -372,6 +361,12 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
0455
0456 /// Callback to store the Geant4 event
0457 void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt) {
0458 +
0459 + if( m_FirstEvent ){
0460 + createCollections( ctxt ) ;
0461 + m_FirstEvent = false ;
0462 + }
0463 +
0464 EventParameters* parameters = context()->event().extension<EventParameters>(false);
0465 int runNumber(0), eventNumber(0);
0466 const int eventNumberOffset(m_eventNumberOffset > 0 ? m_eventNumberOffset : 0);
0467 @@ -381,11 +376,7 @@ void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt) {
0468 if ( parameters ) {
0469 runNumber = parameters->runNumber() + runNumberOffset;
0470 eventNumber = parameters->eventNumber() + eventNumberOffset;
0471 - parameters->extractParameters(m_frame);
0472 -#if podio_VERSION_MAJOR > 0 || podio_VERSION_MINOR > 16 || podio_VERSION_PATCH > 2
0473 - // This functionality is only present in podio > 0.16.2
0474 - eventWeight = m_frame.getParameter<double>("EventWeights");
0475 -#endif
0476 + parameters->extractParameters(*m_store);
0477 } else { // ... or from DD4hep framework
0478 runNumber = m_runNo + runNumberOffset;
0479 eventNumber = ctxt.context->GetEventID() + eventNumberOffset;
0480 @@ -393,14 +384,15 @@ void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt) {
0481 printout(INFO,"Geant4Output2EDM4hep","+++ Saving EDM4hep event %d run %d.", eventNumber, runNumber);
0482
0483 // this does not compile as create() is we only get a const ref - need to review PODIO EventStore API
0484 - edm4hep::EventHeaderCollection header_collection;
0485 - auto header = header_collection.create();
0486 - header.setRunNumber(runNumber);
0487 - header.setEventNumber(eventNumber);
0488 - header.setWeight(eventWeight);
0489 - //not implemented in EDM4hep ? header.setDetectorName(context()->detectorDescription().header().name());
0490 - header.setTimeStamp( std::time(nullptr) ) ;
0491 - m_frame.put( std::move(header_collection), "EventHeader");
0492 + // auto& evtHCol = m_store->get<edm4hep::EventHeaderCollection>("EventHeader") ;
0493 +// auto evtHdr = evtHCol.create() ;
0494 + auto* evtHCol = static_cast<edm4hep::EventHeaderCollection*>(m_collections["EventHeader"]);
0495 + auto evtHdr = evtHCol->create() ;
0496 +
0497 + evtHdr.setRunNumber(runNumber);
0498 + evtHdr.setEventNumber(eventNumber);
0499 +//not implemented in EDM4hep ? evtHdr.setDetectorName(context()->detectorDescription().header().name());
0500 + evtHdr.setTimeStamp( std::time(nullptr) ) ;
0501
0502 saveEventParameters<int>(m_eventParametersInt);
0503 saveEventParameters<float>(m_eventParametersFloat);
0504 @@ -417,34 +409,49 @@ void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt) {
0505
0506 /// Callback to store each Geant4 hit collection
0507 void Geant4Output2EDM4hep::saveCollection(OutputContext<G4Event>& /*ctxt*/, G4VHitsCollection* collection) {
0508 - Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(collection);
0509 +
0510 + size_t nhits = collection->GetSize();
0511 std::string colName = collection->GetName();
0512 +
0513 + printout(DEBUG,"Geant4Output2EDM4hep","+++ Saving EDM4hep collection %s with %d entries.",
0514 + colName.c_str(),int(nhits));
0515 +
0516 + Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(collection);
0517 if( coll == nullptr ){
0518 - error(" no Geant4HitCollection: %s ", colName.c_str());
0519 + printout(ERROR, "Geant4Output2EDM4hep" , " no Geant4HitCollection: %s ", colName.c_str() );
0520 return ;
0521 }
0522 - size_t nhits = collection->GetSize();
0523 +
0524 Geant4ParticleMap* pm = context()->event().extension<Geant4ParticleMap>(false);
0525 - debug("+++ Saving EDM4hep collection %s with %d entries.", colName.c_str(), int(nhits));
0526 +
0527 + auto* mcpc = static_cast<edm4hep::MCParticleCollection*>(m_collections["MCParticles"]);
0528 +
0529 //-------------------------------------------------------------------
0530 if( typeid( Geant4Tracker::Hit ) == coll->type().type() ){
0531 - edm4hep::SimTrackerHitCollection sthc;
0532 +
0533 + auto* sthc = static_cast<edm4hep::SimTrackerHitCollection*>(m_collections[colName]);
0534 +
0535 for(unsigned i=0 ; i < nhits ; ++i){
0536 - auto sth = sthc->create();
0537 + auto sth = sthc->create() ;
0538 +
0539 const Geant4Tracker::Hit* hit = coll->hit(i);
0540 const Geant4Tracker::Hit::Contribution& t = hit->truth;
0541 - int trackID = pm->particleID(t.trackID);
0542 - auto mcp = m_particles.at(trackID);
0543 - const auto& mom = hit->momentum;
0544 - const auto& pos = hit->position;
0545 - edm4hep::Vector3f();
0546 + int trackID = pm->particleID(t.trackID);
0547 +
0548 + auto mcp = mcpc->at(trackID);
0549 +
0550 sth.setCellID( hit->cellID ) ;
0551 sth.setEDep(hit->energyDeposit/CLHEP::GeV);
0552 sth.setPathLength(hit->length/CLHEP::mm);
0553 sth.setTime(hit->truth.time/CLHEP::ns);
0554 sth.setMCParticle(mcp);
0555 - sth.setPosition( {pos.x()/CLHEP::mm, pos.y()/CLHEP::mm, pos.z()/CLHEP::mm} );
0556 - sth.setMomentum( {float(mom.x()/CLHEP::GeV),float(mom.y()/CLHEP::GeV),float(mom.z()/CLHEP::GeV)} );
0557 + sth.setPosition({hit->position.x()/CLHEP::mm,
0558 + hit->position.y()/CLHEP::mm,
0559 + hit->position.z()/CLHEP::mm});
0560 + sth.setMomentum(edm4hep::Vector3f(hit->momentum.x()/CLHEP::GeV,
0561 + hit->momentum.y()/CLHEP::GeV,
0562 + hit->momentum.z()/CLHEP::GeV ));
0563 +
0564 auto particleIt = pm->particles().find(trackID);
0565 if( ( particleIt != pm->particles().end()) ){
0566 // if the original track ID of the particle is not the same as the
0567 @@ -453,49 +460,139 @@ void Geant4Output2EDM4hep::saveCollection(OutputContext<G4Event>& /*ctxt*/, G4VH
0568 sth.setProducedBySecondary( (particleIt->second->originalG4ID != t.trackID) );
0569 }
0570 }
0571 - m_frame.put(std::move(sthc), colName);
0572 - //-------------------------------------------------------------------
0573 + //-------------------------------------------------------------------
0574 }
0575 else if( typeid( Geant4Calorimeter::Hit ) == coll->type().type() ){
0576 - Geant4Sensitive* sd = coll->sensitive();
0577 +
0578 + Geant4Sensitive* sd = coll->sensitive();
0579 int hit_creation_mode = sd->hitCreationMode();
0580
0581 - edm4hep::SimCalorimeterHitCollection sCaloHitColl;
0582 - edm4hep::CaloHitContributionCollection sCaloHitContColl;
0583 + auto* sCaloHitColl = static_cast<edm4hep::SimCalorimeterHitCollection*>(m_collections[colName]);
0584 +
0585 + colName += "Contributions" ;
0586 +
0587 + auto* sCaloHitContColl = static_cast<edm4hep::CaloHitContributionCollection*>(m_collections[colName]);
0588 +
0589 +
0590 for(unsigned i=0 ; i < nhits ; ++i){
0591 auto sch = sCaloHitColl->create() ;
0592 +
0593 const Geant4Calorimeter::Hit* hit = coll->hit(i);
0594 - const auto& pos = hit->position;
0595 +
0596 sch.setCellID( hit->cellID );
0597 - sch.setPosition({float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm)});
0598 + sch.setPosition({
0599 + float(hit->position.x()/CLHEP::mm),
0600 + float(hit->position.y()/CLHEP::mm),
0601 + float(hit->position.z()/CLHEP::mm)});
0602 sch.setEnergy( hit->energyDeposit/CLHEP::GeV );
0603
0604 -
0605 // now add the individual step contributions
0606 - for(auto ci=hit->truth.begin(); ci != hit->truth.end(); ++ci){
0607 + for(Geant4HitData::Contributions::const_iterator ci=hit->truth.begin();
0608 + ci!=hit->truth.end(); ++ci){
0609
0610 auto sCaloHitCont = sCaloHitContColl->create();
0611 sch.addToContributions( sCaloHitCont );
0612
0613 const Geant4HitData::Contribution& c = *ci;
0614 int trackID = pm->particleID(c.trackID);
0615 - auto mcp = m_particles.at(trackID);
0616 + auto mcp = mcpc->at(trackID);
0617 +
0618 sCaloHitCont.setEnergy( c.deposit/CLHEP::GeV );
0619 sCaloHitCont.setTime( c.time/CLHEP::ns );
0620 sCaloHitCont.setParticle( mcp );
0621
0622 if ( hit_creation_mode == Geant4Sensitive::DETAILED_MODE ) {
0623 - edm4hep::Vector3f p(c.x/CLHEP::mm, c.y/CLHEP::mm, c.z/CLHEP::mm);
0624 sCaloHitCont.setPDG( c.pdgID );
0625 - sCaloHitCont.setStepPosition( p );
0626 + sCaloHitCont.setStepPosition( edm4hep::Vector3f(
0627 + c.x/CLHEP::mm,
0628 + c.y/CLHEP::mm,
0629 + c.z/CLHEP::mm) );
0630 }
0631 }
0632 }
0633 - m_frame.put(std::move(sCaloHitColl), colName);
0634 - m_frame.put(std::move(sCaloHitContColl), colName + "Contributions");
0635 - //-------------------------------------------------------------------
0636 + //-------------------------------------------------------------------
0637 } else {
0638 - error("+++ unknown type in Geant4HitCollection %s ", coll->type().type().name());
0639 +
0640 + printout(ERROR, "Geant4Output2EDM4hep" , " unknown type in Geant4HitCollection %s ",
0641 + coll->type().type().name() );
0642 }
0643 }
0644
0645 +void Geant4Output2EDM4hep::createCollections(OutputContext<G4Event>& ctxt){
0646 +
0647 + auto* mcColl = new edm4hep::MCParticleCollection();
0648 + m_collections.emplace("MCParticles", mcColl);
0649 + m_store->registerCollection("MCParticles", mcColl);
0650 + m_file->registerForWrite("MCParticles");
0651 + printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection MCParticles" );
0652 +
0653 + auto* evtHeader = new edm4hep::EventHeaderCollection();
0654 + m_collections.emplace("EventHeader", evtHeader);
0655 + m_store->registerCollection("EventHeader", evtHeader);
0656 + m_file->registerForWrite("EventHeader");
0657 + printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection EventHeader" );
0658 +
0659 +
0660 + const G4Event* evt = ctxt.context ;
0661 + G4HCofThisEvent* hce = evt->GetHCofThisEvent();
0662 + int nCol = hce->GetNumberOfCollections();
0663 +
0664 + for (int i = 0; i < nCol; ++i) {
0665 + G4VHitsCollection* hc = hce->GetHC(i);
0666 + std::string colName = hc->GetName() ;
0667 + Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(hc);
0668 + if( coll == nullptr ){
0669 + printout(WARNING, "Geant4Output2EDM4hep" , " no Geant4HitCollection: %s ", colName.c_str() );
0670 + continue ;
0671 + }
0672 +
0673 + Geant4Sensitive* sd = coll->sensitive();
0674 + string sd_enc = dd4hep::sim::Geant4ConversionHelper::encoding(sd->sensitiveDetector());
0675 +
0676 + if( typeid( Geant4Tracker::Hit ) == coll->type().type() ){
0677 +
0678 + auto* sthc = new edm4hep::SimTrackerHitCollection();
0679 + auto pairColCheck = m_collections.emplace(colName, sthc);
0680 + if ( pairColCheck.second ) {
0681 + m_store->registerCollection(colName, sthc);
0682 + m_file->registerForWrite(colName);
0683 + auto& sthc_md = m_store->getCollectionMetaData( sthc->getID() );
0684 + sthc_md.setValue("CellIDEncodingString", sd_enc);
0685 + printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );
0686 + } else {
0687 + delete sthc;
0688 + }
0689 + printout(DEBUG,"Geant4Output2EDM4hep","+++ re-using collection %s",colName.c_str() );
0690 +
0691 + }
0692 + else if( typeid( Geant4Calorimeter::Hit ) == coll->type().type() ){
0693 +
0694 + auto* schc = new edm4hep::SimCalorimeterHitCollection();
0695 + auto pairColCheck = m_collections.emplace(colName, schc);
0696 + if ( pairColCheck.second ) {
0697 + m_store->registerCollection(colName, schc);
0698 + m_file->registerForWrite(colName);
0699 + auto& schc_md = m_store->getCollectionMetaData( schc->getID() );
0700 + schc_md.setValue("CellIDEncodingString", sd_enc);
0701 + printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );
0702 +
0703 + colName += "Contributions" ;
0704 + auto* chContribColl = new edm4hep::CaloHitContributionCollection();
0705 + m_collections.emplace(colName, chContribColl);
0706 + m_store->registerCollection(colName, chContribColl);
0707 + m_file->registerForWrite(colName);
0708 + printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );
0709 + } else {
0710 + delete schc;
0711 + printout(DEBUG,"Geant4Output2EDM4hep","+++ re-using collection %s",colName.c_str() );
0712 + }
0713 +
0714 + } else {
0715 +
0716 + printout(WARNING, "Geant4Output2EDM4hep" ,
0717 + " unknown type in Geant4HitCollection %s ", coll->type().type().name() );
0718 + }
0719 + }
0720 +
0721 +
0722 +}
