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