From d7240a6fe40aec70e310861396f7e8b3b3f7b457 Mon Sep 17 00:00:00 2001
From: Markus Frank <Markus.Frank@cern.ch>
Date: Mon, 16 Jan 2023 14:53:26 +0100
Subject: [PATCH] Move DDG4 edm4hep writer to use ROOTFrameWriters
---
DDG4/edm4hep/Geant4Output2EDM4hep.cpp | 507 +++++++-----------
examples/ClientTests/scripts/DDG4TestSetup.py | 9 +-
.../ClientTests/scripts/MiniTelGenerate.py | 6 +-
examples/ClientTests/scripts/MiniTelSetup.py | 4 +
4 files changed, 207 insertions(+), 319 deletions(-)
diff --git a/DDG4/edm4hep/Geant4Output2EDM4hep.cpp b/DDG4/edm4hep/Geant4Output2EDM4hep.cpp
index 2bafbdbd1..950f3e209 100644
--- a/DDG4/edm4hep/Geant4Output2EDM4hep.cpp
+++ b/DDG4/edm4hep/Geant4Output2EDM4hep.cpp
@@ -10,40 +10,19 @@
// Author : F.Gaede, DESY
//
//==========================================================================
-
#ifndef DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
#define DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
-// Framework include files
-#include "DD4hep/Detector.h"
-#include "DD4hep/VolumeManager.h"
-#include "DDG4/Geant4HitCollection.h"
-#include "DDG4/Geant4OutputAction.h"
-#include "DDG4/Geant4SensDetAction.h"
-#include "DDG4/Geant4DataConversion.h"
-#include "DDG4/EventParameters.h"
-
-// Geant4 headers
-#include "G4Threading.hh"
-#include "G4AutoLock.hh"
-#include "G4Version.hh"
-
-// use the Geant4 units in namespace CLHEP
-#include "CLHEP/Units/SystemOfUnits.h"
-
-// edm4hep include files
-#include "edm4hep/EventHeaderCollection.h"
-#include "edm4hep/MCParticleCollection.h"
-#include "edm4hep/SimTrackerHitCollection.h"
-#include "edm4hep/CaloHitContributionCollection.h"
-#include "edm4hep/SimCalorimeterHitCollection.h"
-#include "podio/EventStore.h"
-#include "podio/ROOTWriter.h"
-
-#include <typeinfo>
-#include <iostream>
-#include <ctime>
-#include <unordered_map>
+/// Framework include files
+#include <DD4hep/Detector.h>
+#include <DDG4/EventParameters.h>
+#include <DDG4/Geant4OutputAction.h>
+
+/// edm4hep include files
+#include <edm4hep/MCParticleCollection.h>
+/// podio include files
+#include <podio/Frame.h>
+#include <podio/ROOTFrameWriter.h>
/// Namespace for the AIDA detector description toolkit
namespace dd4hep {
@@ -53,20 +32,6 @@ namespace dd4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace sim {
- template <class T=podio::EventStore> void EventParameters::extractParameters(T& event){
- auto& lcparameters = event.getEventMetaData();
-
- for(auto const& ival: this->intParameters()) {
- lcparameters.setValues(ival.first, ival.second);
- }
- for(auto const& ival: this->fltParameters()) {
- lcparameters.setValues(ival.first, ival.second);
- }
- for(auto const& ival: this->strParameters()) {
- lcparameters.setValues(ival.first, ival.second);
- }
- }
-
class Geant4ParticleMap;
/// Base class to output Geant4 event data to EDM4hep
@@ -77,19 +42,22 @@ namespace dd4hep {
*/
class Geant4Output2EDM4hep : public Geant4OutputAction {
protected:
- podio::EventStore* m_store;
- podio::ROOTWriter* m_file;
- int m_runNo;
- int m_runNumberOffset;
- int m_eventNumberOffset;
- std::map< std::string, std::string > m_runHeader;
- std::map< std::string, std::string > m_eventParametersInt;
- std::map< std::string, std::string > m_eventParametersFloat;
- std::map< std::string, std::string > m_eventParametersString;
- bool m_FirstEvent = true ;
-
- std::unordered_map<std::string, podio::CollectionBase*> m_collections;
-
+ using writer_t = podio::ROOTFrameWriter;
+ using stringmap_t = std::map< std::string, std::string >;
+ std::unique_ptr<writer_t> m_file { };
+ podio::Frame m_frame { };
+ edm4hep::MCParticleCollection m_particles { };
+ stringmap_t m_runHeader;
+ stringmap_t m_eventParametersInt;
+ stringmap_t m_eventParametersFloat;
+ stringmap_t m_eventParametersString;
+ std::string m_section_name { "events" };
+ int m_runNo { 0 };
+ int m_runNumberOffset { 0 };
+ int m_eventNo { 0 };
+ int m_eventNumberOffset { 0 };
+ bool m_filesByRun { false };
+
/// create the podio collections for the particles and hits
void createCollections(OutputContext<G4Event>& ctxt) ;
/// Data conversion interface for MC particles to EDM4hep format
@@ -118,33 +86,28 @@ namespace dd4hep {
protected:
/// Fill event parameters in EDM4hep event
template <typename T>
- void saveEventParameters(const std::map<std::string, std::string >& parameters);
+ void saveEventParameters(const std::map<std::string, std::string >& parameters) {
+ for(const auto& p : parameters) {
+ info("Saving event parameter: %-32s = %s", p.first.c_str(), p.second.c_str());
+ m_frame.putParameter(p.first, p.second);
+ }
+ }
};
- /// Fill event parameters in EDM4hep event
- template <typename T>
- inline void Geant4Output2EDM4hep::saveEventParameters(const std::map<std::string, std::string >& parameters) {
- for(std::map<std::string, std::string >::const_iterator iter = parameters.begin(), endIter = parameters.end() ; iter != endIter ; ++iter) {
- T parameter;
- std::istringstream iss(iter->second);
- if ( (iss >> parameter).fail() ) {
- printout(FATAL,"saveEventParameters","+++ Event parameter %s: FAILED to convert to type :%s",iter->first.c_str(),typeid(T).name());
- continue;
- }
- auto& evtMD = m_store->getEventMetaData();
- evtMD.setValue(iter->first,parameter);
+ template <> void EventParameters::extractParameters(podio::Frame& frame) {
+ for(auto const& p: this->intParameters()) {
+ std::cout << "Saving event parameter: " << p.first << std::endl;
+ frame.putParameter(p.first, p.second);
}
- }
-
- /// Fill event parameters in EDM4hep event - std::string specialization
- template <>
- inline void Geant4Output2EDM4hep::saveEventParameters<std::string>(const std::map<std::string, std::string >& parameters) {
- for(std::map<std::string, std::string >::const_iterator iter = parameters.begin(), endIter = parameters.end() ; iter != endIter ; ++iter) {
- auto& evtMD = m_store->getEventMetaData();
- evtMD.setValue(iter->first,iter->second);
+ for(auto const& p: this->fltParameters()) {
+ std::cout << "Saving event parameter: " << p.first << std::endl;
+ frame.putParameter(p.first, p.second);
+ }
+ for(auto const& p: this->strParameters()) {
+ std::cout << "Saving event parameter: " << p.first << std::endl;
+ frame.putParameter(p.first, p.second);
}
}
-
} // End namespace sim
} // End namespace dd4hep
#endif // DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
@@ -162,111 +125,136 @@ namespace dd4hep {
//
//==========================================================================
-// Framework include files
-#include "DD4hep/InstanceCount.h"
-#include "DD4hep/Detector.h"
-#include "DDG4/Geant4HitCollection.h"
-#include "DDG4/Geant4DataConversion.h"
-#include "DDG4/Geant4Context.h"
-#include "DDG4/Geant4Particle.h"
-#include "DDG4/Geant4Data.h"
-#include "DDG4/Geant4Action.h"
-
-//#include "DDG4/Geant4Output2EDM4hep.h"
-#include "G4ParticleDefinition.hh"
-#include "G4VProcess.hh"
-#include "G4Event.hh"
-#include "G4Run.hh"
-
+/// Framework include files
+#include <DD4hep/InstanceCount.h>
+#include <DD4hep/VolumeManager.h>
+
+#include <DDG4/Geant4HitCollection.h>
+#include <DDG4/Geant4DataConversion.h>
+#include <DDG4/Geant4SensDetAction.h>
+#include <DDG4/Geant4Context.h>
+#include <DDG4/Geant4Particle.h>
+#include <DDG4/Geant4Data.h>
+
+///#include <DDG4/Geant4Output2EDM4hep.h>
+/// Geant4 headers
+#include <G4Threading.hh>
+#include <G4AutoLock.hh>
+#include <G4Version.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4VProcess.hh>
+#include <G4Event.hh>
+#include <G4Run.hh>
+/// use the Geant4 units in namespace CLHEP
+#include <CLHEP/Units/SystemOfUnits.h>
+
+/// edm4hep include files
+#include <edm4hep/EventHeaderCollection.h>
+#include <edm4hep/SimTrackerHitCollection.h>
+#include <edm4hep/CaloHitContributionCollection.h>
+#include <edm4hep/SimCalorimeterHitCollection.h>
using namespace dd4hep::sim;
using namespace dd4hep;
-using namespace std;
+
namespace {
G4Mutex action_mutex=G4MUTEX_INITIALIZER;
}
-#include "DDG4/Factories.h"
+#include <DDG4/Factories.h>
DECLARE_GEANT4ACTION(Geant4Output2EDM4hep)
/// Standard constructor
-Geant4Output2EDM4hep::Geant4Output2EDM4hep(Geant4Context* ctxt, const string& nam)
-: Geant4OutputAction(ctxt,nam), m_store(0), m_file(0), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
+Geant4Output2EDM4hep::Geant4Output2EDM4hep(Geant4Context* ctxt, const std::string& nam)
+: Geant4OutputAction(ctxt,nam), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
{
- declareProperty("RunHeader", m_runHeader);
+ declareProperty("RunHeader", m_runHeader);
declareProperty("EventParametersInt", m_eventParametersInt);
declareProperty("EventParametersFloat", m_eventParametersFloat);
declareProperty("EventParametersString", m_eventParametersString);
- declareProperty("RunNumberOffset", m_runNumberOffset);
- declareProperty("EventNumberOffset", m_eventNumberOffset);
- printout( INFO, "Geant4Output2EDM4hep" ," instantiated ..." ) ;
+ declareProperty("RunNumberOffset", m_runNumberOffset);
+ declareProperty("EventNumberOffset", m_eventNumberOffset);
+ declareProperty("SectionName", m_section_name);
+ declareProperty("FilesByRun", m_filesByRun);
+ info("Writer is now instantiated ..." );
InstanceCount::increment(this);
}
/// Default destructor
Geant4Output2EDM4hep::~Geant4Output2EDM4hep() {
G4AutoLock protection_lock(&action_mutex);
- if ( m_file ) {
- m_file->finish();
- detail::deletePtr(m_file);
- }
- if (nullptr != m_store) {
- delete m_store;
- }
+ m_file.reset();
InstanceCount::decrement(this);
}
// Callback to store the Geant4 run information
void Geant4Output2EDM4hep::beginRun(const G4Run* run) {
G4AutoLock protection_lock(&action_mutex);
- if ( 0 == m_file && !m_output.empty() ) {
- m_store = new podio::EventStore ;
- m_file = new podio::ROOTWriter(m_output, m_store);
-
- printout( INFO, "Geant4Output2EDM4hep" ," opened %s for output", m_output.c_str() ) ;
+ std::string fname = m_output;
+ m_runNo = run->GetRunID();
+ if ( m_filesByRun ) {
+ std::size_t idx = m_output.rfind(".");
+ if ( idx != std::string::npos ) {
+ fname = m_output.substr(0, idx) + _toString(m_runNo, ".run%08d") + m_output.substr(idx);
+ }
+ }
+ if ( !m_file && !fname.empty() ) {
+ m_file = std::make_unique<podio::ROOTFrameWriter>(fname);
+ if ( !m_file ) {
+ fatal("+++ Failed to open output file: %s", fname.c_str());
+ }
+ printout( INFO, "Geant4Output2EDM4hep" ,"Opened %s for output", fname.c_str() ) ;
}
-
- saveRun(run);
}
/// Callback to store the Geant4 run information
-void Geant4Output2EDM4hep::endRun(const G4Run* /*run*/) {
- // saveRun(run);
+void Geant4Output2EDM4hep::endRun(const G4Run* run) {
+ saveRun(run);
+ if ( m_file ) {
+ m_file->finish();
+ m_file.reset();
+ }
}
/// Commit data at end of filling procedure
void Geant4Output2EDM4hep::commit( OutputContext<G4Event>& /* ctxt */) {
if ( m_file ) {
G4AutoLock protection_lock(&action_mutex);
- m_file->writeEvent();
- m_store->clearCollections();
+ m_frame.put( std::move(m_particles), "MCParticles");
+ m_file->writeFrame(m_frame, m_section_name);
+ m_particles.clear();
+ m_frame = {};
return;
}
except("+++ Failed to write output file. [Stream is not open]");
}
/// Callback to store the Geant4 run information
-void Geant4Output2EDM4hep::saveRun(const G4Run* /*run*/) {
+void Geant4Output2EDM4hep::saveRun(const G4Run* run) {
//G4AutoLock protection_lock(&action_mutex);
- printout( WARNING, "Geant4Output2EDM4hep" ,"saveRun(): RunHeader not implemented in EDM4hep, nothing written ..." ) ;
+ warning("saveRun: RunHeader not implemented in EDM4hep, nothing written for run: %d", run->GetRunID());
// --- write an edm4hep::RunHeader ---------
//FIXME: need a suitable Runheader object in EDM4hep
-
-// edm4hep::LCRunHeaderImpl* rh = new edm4hep::LCRunHeaderImpl;
-// for (std::map< std::string, std::string >::iterator it = m_runHeader.begin(); it != m_runHeader.end(); ++it) {
-// rh->parameters().setValue( it->first, it->second );
-// }
-// m_runNo = m_runNumberOffset > 0 ? m_runNumberOffset + run->GetRunID() : run->GetRunID();
-// rh->parameters().setValue("GEANT4Version", G4Version);
-// rh->parameters().setValue("DD4HEPVersion", versionString());
-// rh->setRunNumber(m_runNo);
-// rh->setDetectorName(context()->detectorDescription().header().name());
-// m_file->writeRunHeader(rh);
+ // edm4hep::LCRunHeaderImpl* rh = new edm4hep::LCRunHeaderImpl;
+ // for (stringmap_t::iterator it = m_runHeader.begin(); it != m_runHeader.end(); ++it) {
+ // rh->parameters().setValue( it->first, it->second );
+ // }
+ // m_runNo = m_runNumberOffset > 0 ? m_runNumberOffset + run->GetRunID() : run->GetRunID();
+ // rh->parameters().setValue("GEANT4Version", G4Version);
+ // rh->parameters().setValue("DD4HEPVersion", versionString());
+ // rh->setRunNumber(m_runNo);
+ // rh->setDetectorName(context()->detectorDescription().header().name());
+ // m_file->writeRunHeader(rh);
}
-void Geant4Output2EDM4hep::begin(const G4Event* /* event */) {
+void Geant4Output2EDM4hep::begin(const G4Event* event) {
+ /// Create event frame object
+ m_eventNo = event->GetEventID();
+ m_particles.clear();
+ m_particles = {};
+ m_frame = {};
}
/// Data conversion interface for MC particles to EDM4hep format
@@ -274,13 +262,13 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
typedef detail::ReferenceBitMask<const int> PropertyMask;
typedef Geant4ParticleMap::ParticleMap ParticleMap;
const ParticleMap& pm = particles->particleMap;
- auto mcpc = static_cast<edm4hep::MCParticleCollection*>(m_collections["MCParticles"]);
+ m_particles.clear();
if ( pm.size() > 0 ) {
size_t cnt = 0;
// Mapping of ids in the ParticleMap to indices in the MCParticle collection
- map<int,int> p_ids;
- vector<const Geant4Particle*> p_part;
+ std::map<int,int> p_ids;
+ std::vector<const Geant4Particle*> p_part;
p_part.reserve(pm.size());
// First create the particles
for (const auto& iParticle : pm) {
@@ -289,7 +277,7 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
PropertyMask mask(p->status);
// std::cout << " ********** mcp status : 0x" << std::hex << p->status << ", mask.isSet(G4PARTICLE_GEN_STABLE) x" << std::dec << mask.isSet(G4PARTICLE_GEN_STABLE) <<std::endl ;
const G4ParticleDefinition* def = p.definition();
- auto mcp = mcpc->create();
+ auto mcp = m_particles.create();
mcp.setPDG(p->pdgID);
float ps_fa[3] = {float(p->psx/CLHEP::GeV),float(p->psy/CLHEP::GeV),float(p->psz/CLHEP::GeV)};
@@ -311,13 +299,13 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
// Set generator status
mcp.setGeneratorStatus(0);
if( p->genStatus ) {
- mcp.setGeneratorStatus( p->genStatus ) ;
+ mcp.setGeneratorStatus( p->genStatus ) ;
} else {
- if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) mcp.setGeneratorStatus(1);
- else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) ) mcp.setGeneratorStatus(2);
- else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) mcp.setGeneratorStatus(3);
- else if ( mask.isSet(G4PARTICLE_GEN_BEAM) ) mcp.setGeneratorStatus(4);
- else if ( mask.isSet(G4PARTICLE_GEN_OTHER) ) mcp.setGeneratorStatus(9);
+ if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) mcp.setGeneratorStatus(1);
+ else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) ) mcp.setGeneratorStatus(2);
+ else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) mcp.setGeneratorStatus(3);
+ else if ( mask.isSet(G4PARTICLE_GEN_BEAM) ) mcp.setGeneratorStatus(4);
+ else if ( mask.isSet(G4PARTICLE_GEN_OTHER) ) mcp.setGeneratorStatus(9);
}
// Set simulation status
@@ -332,7 +320,7 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
//fg: if simstatus !=0 we have to set the generator status to 0:
if( mcp.isCreatedInSimulation() )
- mcp.setGeneratorStatus( 0 ) ;
+ mcp.setGeneratorStatus( 0 ) ;
mcp.setSpin(p->spin);
mcp.setColorFlow(p->colorFlow);
@@ -344,30 +332,30 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
// Now establish parent-daughter relationships
for(size_t i=0; i < p_ids.size(); ++i) {
const Geant4Particle* p = p_part[i];
- auto q = (*mcpc)[i];
+ auto q = m_particles[i];
for (const auto& idau : p->daughters) {
- const auto k = p_ids.find(idau);
- if (k == p_ids.end()) {
- printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
- continue;
- }
- int iqdau = (*k).second;
- auto qdau = (*mcpc)[iqdau];
- q.addToDaughters(qdau);
+ const auto k = p_ids.find(idau);
+ if (k == p_ids.end()) {
+ fatal("+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
+ continue;
+ }
+ int iqdau = (*k).second;
+ auto qdau = m_particles[iqdau];
+ q.addToDaughters(qdau);
}
for (const auto& ipar : p->parents) {
- if (ipar >= 0) { // A parent ID of -1 means NO parent, because a base of 0 is perfectly legal
- const auto k = p_ids.find(ipar);
- if (k == p_ids.end()) {
- printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
- continue;
- }
- int iqpar = (*k).second;
- auto qpar = (*mcpc)[iqpar];
- q.addToParents(qpar);
- }
+ if (ipar >= 0) { // A parent ID of -1 means NO parent, because a base of 0 is perfectly legal
+ const auto k = p_ids.find(ipar);
+ if (k == p_ids.end()) {
+ fatal("+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
+ continue;
+ }
+ int iqpar = (*k).second;
+ auto qpar = m_particles[iqpar];
+ q.addToParents(qpar);
+ }
}
}
}
@@ -375,12 +363,6 @@ void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles) {
/// Callback to store the Geant4 event
void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt) {
-
- if( m_FirstEvent ){
- createCollections( ctxt ) ;
- m_FirstEvent = false ;
- }
-
EventParameters* parameters = context()->event().extension<EventParameters>(false);
int runNumber(0), eventNumber(0);
const int eventNumberOffset(m_eventNumberOffset > 0 ? m_eventNumberOffset : 0);
@@ -389,7 +371,7 @@ void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt) {
if ( parameters ) {
runNumber = parameters->runNumber() + runNumberOffset;
eventNumber = parameters->eventNumber() + eventNumberOffset;
- parameters->extractParameters(*m_store);
+ parameters->extractParameters(m_frame);
} else { // ... or from DD4hep framework
runNumber = m_runNo + runNumberOffset;
eventNumber = ctxt.context->GetEventID() + eventNumberOffset;
@@ -397,15 +379,13 @@ void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt) {
printout(INFO,"Geant4Output2EDM4hep","+++ Saving EDM4hep event %d run %d.", eventNumber, runNumber);
// this does not compile as create() is we only get a const ref - need to review PODIO EventStore API
- // auto& evtHCol = m_store->get<edm4hep::EventHeaderCollection>("EventHeader") ;
-// auto evtHdr = evtHCol.create() ;
- auto* evtHCol = static_cast<edm4hep::EventHeaderCollection*>(m_collections["EventHeader"]);
- auto evtHdr = evtHCol->create() ;
-
- evtHdr.setRunNumber(runNumber);
- evtHdr.setEventNumber(eventNumber);
-//not implemented in EDM4hep ? evtHdr.setDetectorName(context()->detectorDescription().header().name());
- evtHdr.setTimeStamp( std::time(nullptr) ) ;
+ edm4hep::EventHeaderCollection header_collection;
+ auto header = header_collection.create();
+ header.setRunNumber(runNumber);
+ header.setEventNumber(eventNumber);
+ //not implemented in EDM4hep ? header.setDetectorName(context()->detectorDescription().header().name());
+ header.setTimeStamp( std::time(nullptr) ) ;
+ m_frame.put( std::move(header_collection), "EventHeader");
saveEventParameters<int>(m_eventParametersInt);
saveEventParameters<float>(m_eventParametersFloat);
@@ -422,49 +402,34 @@ void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt) {
/// Callback to store each Geant4 hit collection
void Geant4Output2EDM4hep::saveCollection(OutputContext<G4Event>& /*ctxt*/, G4VHitsCollection* collection) {
-
- size_t nhits = collection->GetSize();
- std::string colName = collection->GetName();
-
- printout(DEBUG,"Geant4Output2EDM4hep","+++ Saving EDM4hep collection %s with %d entries.",
- colName.c_str(),int(nhits));
-
Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(collection);
+ std::string colName = collection->GetName();
if( coll == nullptr ){
- printout(ERROR, "Geant4Output2EDM4hep" , " no Geant4HitCollection: %s ", colName.c_str() );
+ error(" no Geant4HitCollection: %s ", colName.c_str());
return ;
}
-
+ size_t nhits = collection->GetSize();
Geant4ParticleMap* pm = context()->event().extension<Geant4ParticleMap>(false);
-
- auto* mcpc = static_cast<edm4hep::MCParticleCollection*>(m_collections["MCParticles"]);
-
+ debug("+++ Saving EDM4hep collection %s with %d entries.", colName.c_str(), int(nhits));
//-------------------------------------------------------------------
if( typeid( Geant4Tracker::Hit ) == coll->type().type() ){
-
- auto* sthc = static_cast<edm4hep::SimTrackerHitCollection*>(m_collections[colName]);
-
+ edm4hep::SimTrackerHitCollection sthc;
for(unsigned i=0 ; i < nhits ; ++i){
- auto sth = sthc->create() ;
-
+ auto sth = sthc->create();
const Geant4Tracker::Hit* hit = coll->hit(i);
const Geant4Tracker::Hit::Contribution& t = hit->truth;
- int trackID = pm->particleID(t.trackID);
-
- auto mcp = mcpc->at(trackID);
-
+ int trackID = pm->particleID(t.trackID);
+ auto mcp = m_particles.at(trackID);
+ const auto& mom = hit->momentum;
+ const auto& pos = hit->position;
+ edm4hep::Vector3f();
sth.setCellID( hit->cellID ) ;
sth.setEDep(hit->energyDeposit/CLHEP::GeV);
sth.setPathLength(hit->length/CLHEP::mm);
sth.setTime(hit->truth.time/CLHEP::ns);
sth.setMCParticle(mcp);
- sth.setPosition({hit->position.x()/CLHEP::mm,
- hit->position.y()/CLHEP::mm,
- hit->position.z()/CLHEP::mm});
- sth.setMomentum(edm4hep::Vector3f(hit->momentum.x()/CLHEP::GeV,
- hit->momentum.y()/CLHEP::GeV,
- hit->momentum.z()/CLHEP::GeV ));
-
+ sth.setPosition( {pos.x()/CLHEP::mm, pos.y()/CLHEP::mm, pos.z()/CLHEP::mm} );
+ sth.setMomentum( {float(mom.x()/CLHEP::GeV),float(mom.y()/CLHEP::GeV),float(mom.z()/CLHEP::GeV)} );
auto particleIt = pm->particles().find(trackID);
if( ( particleIt != pm->particles().end()) ){
// if the original track ID of the particle is not the same as the
@@ -473,139 +438,49 @@ void Geant4Output2EDM4hep::saveCollection(OutputContext<G4Event>& /*ctxt*/, G4VH
sth.setProducedBySecondary( (particleIt->second->originalG4ID != t.trackID) );
}
}
- //-------------------------------------------------------------------
+ m_frame.put(std::move(sthc), colName);
+ //-------------------------------------------------------------------
}
else if( typeid( Geant4Calorimeter::Hit ) == coll->type().type() ){
-
- Geant4Sensitive* sd = coll->sensitive();
+ Geant4Sensitive* sd = coll->sensitive();
int hit_creation_mode = sd->hitCreationMode();
- auto* sCaloHitColl = static_cast<edm4hep::SimCalorimeterHitCollection*>(m_collections[colName]);
-
- colName += "Contributions" ;
-
- auto* sCaloHitContColl = static_cast<edm4hep::CaloHitContributionCollection*>(m_collections[colName]);
-
-
+ edm4hep::SimCalorimeterHitCollection sCaloHitColl;
+ edm4hep::CaloHitContributionCollection sCaloHitContColl;
for(unsigned i=0 ; i < nhits ; ++i){
auto sch = sCaloHitColl->create() ;
-
const Geant4Calorimeter::Hit* hit = coll->hit(i);
-
+ const auto& pos = hit->position;
sch.setCellID( hit->cellID );
- sch.setPosition({
- float(hit->position.x()/CLHEP::mm),
- float(hit->position.y()/CLHEP::mm),
- float(hit->position.z()/CLHEP::mm)});
+ sch.setPosition({float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm)});
sch.setEnergy( hit->energyDeposit/CLHEP::GeV );
+
// now add the individual step contributions
- for(Geant4HitData::Contributions::const_iterator ci=hit->truth.begin();
- ci!=hit->truth.end(); ++ci){
+ for(auto ci=hit->truth.begin(); ci != hit->truth.end(); ++ci){
auto sCaloHitCont = sCaloHitContColl->create();
sch.addToContributions( sCaloHitCont );
const Geant4HitData::Contribution& c = *ci;
int trackID = pm->particleID(c.trackID);
- auto mcp = mcpc->at(trackID);
-
+ auto mcp = m_particles.at(trackID);
sCaloHitCont.setEnergy( c.deposit/CLHEP::GeV );
sCaloHitCont.setTime( c.time/CLHEP::ns );
sCaloHitCont.setParticle( mcp );
if ( hit_creation_mode == Geant4Sensitive::DETAILED_MODE ) {
+ edm4hep::Vector3f p(c.x/CLHEP::mm, c.y/CLHEP::mm, c.z/CLHEP::mm);
sCaloHitCont.setPDG( c.pdgID );
- sCaloHitCont.setStepPosition( edm4hep::Vector3f(
- c.x/CLHEP::mm,
- c.y/CLHEP::mm,
- c.z/CLHEP::mm) );
+ sCaloHitCont.setStepPosition( p );
}
}
}
- //-------------------------------------------------------------------
+ m_frame.put(std::move(sCaloHitColl), colName);
+ m_frame.put(std::move(sCaloHitContColl), colName + "Contributions");
+ //-------------------------------------------------------------------
} else {
-
- printout(ERROR, "Geant4Output2EDM4hep" , " unknown type in Geant4HitCollection %s ",
- coll->type().type().name() );
+ error("+++ unknown type in Geant4HitCollection %s ", coll->type().type().name());
}
}
-void Geant4Output2EDM4hep::createCollections(OutputContext<G4Event>& ctxt){
-
- auto* mcColl = new edm4hep::MCParticleCollection();
- m_collections.emplace("MCParticles", mcColl);
- m_store->registerCollection("MCParticles", mcColl);
- m_file->registerForWrite("MCParticles");
- printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection MCParticles" );
-
- auto* evtHeader = new edm4hep::EventHeaderCollection();
- m_collections.emplace("EventHeader", evtHeader);
- m_store->registerCollection("EventHeader", evtHeader);
- m_file->registerForWrite("EventHeader");
- printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection EventHeader" );
-
-
- const G4Event* evt = ctxt.context ;
- G4HCofThisEvent* hce = evt->GetHCofThisEvent();
- int nCol = hce->GetNumberOfCollections();
-
- for (int i = 0; i < nCol; ++i) {
- G4VHitsCollection* hc = hce->GetHC(i);
- std::string colName = hc->GetName() ;
- Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(hc);
- if( coll == nullptr ){
- printout(WARNING, "Geant4Output2EDM4hep" , " no Geant4HitCollection: %s ", colName.c_str() );
- continue ;
- }
-
- Geant4Sensitive* sd = coll->sensitive();
- string sd_enc = dd4hep::sim::Geant4ConversionHelper::encoding(sd->sensitiveDetector());
-
- if( typeid( Geant4Tracker::Hit ) == coll->type().type() ){
-
- auto* sthc = new edm4hep::SimTrackerHitCollection();
- auto pairColCheck = m_collections.emplace(colName, sthc);
- if ( pairColCheck.second ) {
- m_store->registerCollection(colName, sthc);
- m_file->registerForWrite(colName);
- auto& sthc_md = m_store->getCollectionMetaData( sthc->getID() );
- sthc_md.setValue("CellIDEncodingString", sd_enc);
- printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );
- } else {
- delete sthc;
- }
- printout(DEBUG,"Geant4Output2EDM4hep","+++ re-using collection %s",colName.c_str() );
-
- }
- else if( typeid( Geant4Calorimeter::Hit ) == coll->type().type() ){
-
- auto* schc = new edm4hep::SimCalorimeterHitCollection();
- auto pairColCheck = m_collections.emplace(colName, schc);
- if ( pairColCheck.second ) {
- m_store->registerCollection(colName, schc);
- m_file->registerForWrite(colName);
- auto& schc_md = m_store->getCollectionMetaData( schc->getID() );
- schc_md.setValue("CellIDEncodingString", sd_enc);
- printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );
-
- colName += "Contributions" ;
- auto* chContribColl = new edm4hep::CaloHitContributionCollection();
- m_collections.emplace(colName, chContribColl);
- m_store->registerCollection(colName, chContribColl);
- m_file->registerForWrite(colName);
- printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );
- } else {
- delete schc;
- printout(DEBUG,"Geant4Output2EDM4hep","+++ re-using collection %s",colName.c_str() );
- }
-
- } else {
-
- printout(WARNING, "Geant4Output2EDM4hep" ,
- " unknown type in Geant4HitCollection %s ", coll->type().type().name() );
- }
- }
-
-
-}
diff --git a/examples/ClientTests/scripts/DDG4TestSetup.py b/examples/ClientTests/scripts/DDG4TestSetup.py
index e25039056..4f1664777 100644
--- a/examples/ClientTests/scripts/DDG4TestSetup.py
+++ b/examples/ClientTests/scripts/DDG4TestSetup.py
@@ -46,8 +46,13 @@ class Setup:
def defineOutput(self, output):
# Configure I/O
- evt_root = self.geant4.setupROOTOutput('RootOutput', output, mc_truth=True)
- return evt_root
+ evt_write = self.geant4.setupROOTOutput('RootOutput', output, mc_truth=True)
+ return evt_write
+
+ def defineEdm4hepOutput(self, output):
+ # Configure I/O
+ evt_write = self.geant4.setupEDM4hepOutput('Edm4hepOutput', output)
+ return evt_write
def setupGun(self, name="Gun", particle='pi-', energy=100 * GeV, multiplicity=1):
# Setup particle gun
diff --git a/examples/ClientTests/scripts/MiniTelGenerate.py b/examples/ClientTests/scripts/MiniTelGenerate.py
index 4c7242a61..c177f8ba7 100644
--- a/examples/ClientTests/scripts/MiniTelGenerate.py
+++ b/examples/ClientTests/scripts/MiniTelGenerate.py
@@ -38,8 +38,12 @@ def run():
cmds.append('/ddg4/UI/terminate')
m.ui.Commands = cmds
m.configure()
- wr = m.defineOutput(output='MiniTel')
+ if args.output and str(args.output).lower().find('edm4hep') >= 0:
+ wr = m.defineEdm4hepOutput(args.output)
+ else:
+ wr = m.defineOutput(output='MiniTel')
wr.FilesByRun = True
+
gen = DDG4.GeneratorAction(kernel, "Geant4GeneratorActionInit/GenerationInit")
kernel.generatorAction().adopt(gen)
diff --git a/examples/ClientTests/scripts/MiniTelSetup.py b/examples/ClientTests/scripts/MiniTelSetup.py
index df88967d3..f1c5a9b16 100644
--- a/examples/ClientTests/scripts/MiniTelSetup.py
+++ b/examples/ClientTests/scripts/MiniTelSetup.py
@@ -41,3 +41,7 @@ class Setup(DDG4TestSetup.Setup):
def defineOutput(self, output='MiniTel_' + time.strftime('%Y-%m-%d_%H-%M')):
return DDG4TestSetup.Setup.defineOutput(self, output)
+
+ def defineEdm4hepOutput(self, output='MiniTel_' + time.strftime('%Y-%m-%d_%H-%M')):
+ return DDG4TestSetup.Setup.defineEdm4hepOutput(self, output)
+
--
GitLab