13 #ifndef DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
14 #define DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
24 #include <edm4hep/MCParticleCollection.h>
25 #include <edm4hep/SimTrackerHitCollection.h>
26 #include <edm4hep/CaloHitContributionCollection.h>
27 #include <edm4hep/SimCalorimeterHitCollection.h>
28 #include <edm4hep/EDM4hepVersion.h>
29 #include <edm4hep/Constants.h>
30 #if EDM4HEP_BUILD_VERSION < EDM4HEP_VERSION(0, 99, 0)
31 using edm4hep::CellIDEncoding;
33 using edm4hep::labels::CellIDEncoding;
35 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 3)
36 #include <edm4hep/GeneratorEventParametersCollection.h>
40 #include <podio/CollectionBase.h>
41 #include <podio/podioVersion.h>
42 #include <podio/Frame.h>
43 #include <podio/FrameCategories.h>
44 #if PODIO_BUILD_VERSION >= PODIO_VERSION(1, 0, 0)
45 #include <podio/Writer.h>
47 #if PODIO_BUILD_VERSION >= PODIO_VERSION(0, 99, 0)
48 #include <podio/ROOTWriter.h>
50 #include <podio/ROOTFrameWriter.h>
52 using ROOTWriter = podio::ROOTFrameWriter;
77 #if PODIO_BUILD_VERSION >= PODIO_VERSION(1, 0, 0)
85 using trackermap_t = std::map< std::string, edm4hep::SimTrackerHitCollection >;
86 using calorimeterpair_t = std::pair< edm4hep::SimCalorimeterHitCollection, edm4hep::CaloHitContributionCollection >;
88 std::unique_ptr<writer_t>
m_file { };
120 virtual void beginRun(
const G4Run* run);
122 virtual void endRun(
const G4Run* run);
125 virtual void saveRun(
const G4Run* run);
127 virtual void saveEvent( OutputContext<G4Event>& ctxt);
131 virtual void commit( OutputContext<G4Event>& ctxt);
134 virtual void begin(
const G4Event* event);
137 template <
typename T>
139 for(
const auto& p : parameters) {
140 std::stringstream output;
141 output <<
"Saving event parameter: "
142 << std::setw(32) << p.first
143 << std::setw(20) << p.second;
144 info(output.str().c_str());
145 m_frame.putParameter(p.first, p.second);
152 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving event parameter: %s", p.first.c_str());
153 frame.putParameter(p.first, p.second);
156 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving event parameter: %s", p.first.c_str());
157 frame.putParameter(p.first, p.second);
160 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving event parameter: %s", p.first.c_str());
161 frame.putParameter(p.first, p.second);
165 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving event parameter: %s", p.first.c_str());
166 frame.putParameter(p.first, p.second);
172 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving run parameter: %s", p.first.c_str());
173 frame.putParameter(p.first, p.second);
176 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving run parameter: %s", p.first.c_str());
177 frame.putParameter(p.first, p.second);
180 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving run parameter: %s", p.first.c_str());
181 frame.putParameter(p.first, p.second);
185 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving run parameter: %s", p.first.c_str());
186 frame.putParameter(p.first, p.second);
191 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving meta parameter: %s", p.first.c_str());
192 frame.putParameter(p.first, p.second);
195 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving meta parameter: %s", p.first.c_str());
196 frame.putParameter(p.first, p.second);
199 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving meta parameter: %s", p.first.c_str());
200 frame.putParameter(p.first, p.second);
204 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving meta parameter: %s", p.first.c_str());
205 frame.putParameter(p.first, p.second);
211 #endif // DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
239 #include <G4Threading.hh>
240 #include <G4AutoLock.hh>
241 #include <G4Version.hh>
242 #include <G4ParticleDefinition.hh>
243 #include <G4VProcess.hh>
244 #include <G4Event.hh>
247 #include <CLHEP/Units/SystemOfUnits.h>
250 #include <edm4hep/EventHeaderCollection.h>
256 G4Mutex action_mutex = G4MUTEX_INITIALIZER;
264 :
Geant4OutputAction(ctxt,nam), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
266 declareProperty(
"RunHeader", m_runHeader);
267 declareProperty(
"EventParametersInt", m_eventParametersInt);
268 declareProperty(
"EventParametersFloat", m_eventParametersFloat);
269 declareProperty(
"EventParametersString", m_eventParametersString);
270 declareProperty(
"RunParametersInt", m_runParametersInt);
271 declareProperty(
"RunParametersFloat", m_runParametersFloat);
272 declareProperty(
"RunParametersString", m_runParametersString);
273 declareProperty(
"RunNumberOffset", m_runNumberOffset);
274 declareProperty(
"EventNumberOffset", m_eventNumberOffset);
275 declareProperty(
"SectionName", m_section_name);
276 declareProperty(
"FilesByRun", m_filesByRun);
277 declareProperty(
"RNTuple", m_rntuple);
279 info(
"Writer is now instantiated ..." );
285 G4AutoLock protection_lock(&action_mutex);
291 G4AutoLock protection_lock(&action_mutex);
295 std::size_t idx =
m_output.rfind(
".");
296 if ( idx != std::string::npos ) {
301 if ( !fname.empty() && !
m_file ) {
302 #if PODIO_BUILD_VERSION >= PODIO_VERSION(1, 0, 0)
303 m_file = std::make_unique<podio::Writer>(podio::makeWriter(fname,
m_rntuple ?
"rntuple" :
"default"));
305 m_file = std::make_unique<podio::ROOTWriter>(fname);
308 fatal(
"+++ Failed to open output file: %s", fname.c_str());
310 printout( INFO,
"Geant4Output2EDM4hep" ,
"Opened %s for output", fname.c_str() ) ;
323 G4AutoLock protection_lock(&action_mutex);
332 podio::Frame metaFrame{};
334 metaFrame.putParameter(podio::collMetadataParamName(
name, CellIDEncoding), encodingStr);
336 G4AutoLock protection_lock(&action_mutex);
337 m_file->writeFrame(metaFrame,
"metadata");
343 G4AutoLock protection_lock(&action_mutex);
346 m_frame.put( std::move(it->second), it->first);
349 m_frame.put( std::move(calorimeterHits.first), colName);
350 m_frame.put( std::move(calorimeterHits.second), colName +
"Contributions");
359 except(
"+++ Failed to write output file. [Stream is not open]");
364 G4AutoLock protection_lock(&action_mutex);
368 podio::Frame runHeader {};
370 runHeader.putParameter(
key, value);
373 runHeader.putParameter(
key, value);
376 runHeader.putParameter(
key, value);
379 runHeader.putParameter(
key, value);
382 runHeader.putParameter(
"runNumber",
m_runNo);
383 runHeader.putParameter(
"GEANT4Version", G4Version);
385 runHeader.putParameter(
"detectorName",
context()->detectorDescription().header().
name());
388 if (
context()->runPtr() !=
nullptr) {
393 m_file->writeFrame(runHeader,
"runs");
399 podio::Frame metaFrame {};
404 m_file->writeFrame(metaFrame,
"meta");
420 typedef detail::ReferenceBitMask<const int>
PropertyMask;
425 if ( pm.size() > 0 ) {
428 std::map<int,int> p_ids;
429 std::vector<const Geant4Particle*> p_part;
430 p_part.reserve(pm.size());
432 for (
const auto& iParticle : pm) {
433 int id = iParticle.first;
437 const G4ParticleDefinition* def = p.
definition();
439 mcp.setPDG(p->
pdgID);
441 using MT = decltype(std::declval<edm4hep::MCParticle>().getMomentum().x);
442 mcp.setMomentum( {MT(p->
psx/CLHEP::GeV),MT(p->
psy/CLHEP::GeV),MT(p->
psz/CLHEP::GeV)} );
443 mcp.setMomentumAtEndpoint( {MT(p->
pex/CLHEP::GeV),MT(p->
pey/CLHEP::GeV),MT(p->
pez/CLHEP::GeV)} );
445 double vs_fa[3] = { p->
vsx/CLHEP::mm, p->
vsy/CLHEP::mm, p->
vsz/CLHEP::mm } ;
446 mcp.setVertex( vs_fa );
448 double ve_fa[3] = { p->
vex/CLHEP::mm, p->
vey/CLHEP::mm, p->
vez/CLHEP::mm } ;
449 mcp.setEndpoint( ve_fa );
451 mcp.setTime(p->
time/CLHEP::ns);
452 mcp.setMass(p->
mass/CLHEP::GeV);
453 mcp.setCharge(def ? def->GetPDGCharge() : 0);
456 mcp.setGeneratorStatus(0);
475 mcp.setOverlay(
false );
478 if( mcp.isCreatedInSimulation() )
479 mcp.setGeneratorStatus( 0 ) ;
481 #if EDM4HEP_MCPARTICLE_HAS_HELICITY
482 mcp.setHelicity(p->
spin[2]);
484 mcp.setSpin(p->
spin);
492 for(
size_t i=0; i < p_ids.size(); ++i) {
497 const auto k = p_ids.find(idau);
498 if (k == p_ids.end()) {
499 fatal(
"+++ Particle %d: FAILED to find daughter with ID:%d",p->
id,idau);
502 int iqdau = (*k).second;
504 q.addToDaughters(qdau);
507 for (
const auto& ipar : p->
parents) {
509 const auto k = p_ids.find(ipar);
510 if (k == p_ids.end()) {
511 fatal(
"+++ Particle %d: FAILED to find parent with ID:%d",p->
id,ipar);
514 int iqpar = (*k).second;
516 q.addToParents(qpar);
526 int runNumber(0), eventNumber(0);
529 double eventWeight{0};
532 runNumber = parameters->
runNumber() + runNumberOffset;
533 eventNumber = parameters->
eventNumber() + eventNumberOffset;
535 #if PODIO_BUILD_VERSION > PODIO_VERSION(0, 99, 0)
536 eventWeight =
m_frame.getParameter<
double>(
"EventWeights").value_or(0.0);
538 eventWeight =
m_frame.getParameter<
double>(
"EventWeights");
541 runNumber =
m_runNo + runNumberOffset;
542 eventNumber = ctxt.
context->GetEventID() + eventNumberOffset;
544 printout(INFO,
"Geant4Output2EDM4hep",
"+++ Saving EDM4hep event %d run %d.", eventNumber, runNumber);
547 edm4hep::EventHeaderCollection header_collection;
549 auto header = header_collection.create();
550 header.setRunNumber(runNumber);
551 header.setEventNumber(eventNumber);
552 header.setWeight(eventWeight);
554 header.setTimeStamp(std::time(
nullptr));
559 header.setTimeStamp(meh->getTimeStamp());
560 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
561 for (
auto const& weight: meh->getWeights()) {
562 header.addToWeights(weight);
567 m_frame.put(std::move(header_collection),
"EventHeader");
569 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 3)
573 edm4hep::GeneratorEventParametersCollection genEvtParamsColl{};
574 genEvtParamsColl.push_back(*genEvtParams);
575 m_frame.put(std::move(genEvtParamsColl), edm4hep::labels::GeneratorEventParameters);
585 print(
"+++ Saving %d EDM4hep particles....",
int(part_map->
particleMap.size()));
602 operator std::string()
const {
603 const auto* sd = m_coll->sensitive();
614 std::string colName = collection->GetName();
615 if( coll ==
nullptr ){
616 error(
" no Geant4HitCollection: %s ", colName.c_str());
619 size_t nhits = collection->GetSize();
621 debug(
"+++ Saving EDM4hep collection %s with %d entries.", colName.c_str(),
int(nhits));
630 for(
unsigned i=0 ; i < nhits ; ++i){
631 auto sth = hits.create();
639 sth.setCellID( hit->
cellID ) ;
641 sth.setPathLength(hit->
length/CLHEP::mm);
643 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 10, 99)
644 sth.setParticle(mcp);
646 sth.setMCParticle(mcp);
648 sth.setPosition( {pos.x()/CLHEP::mm, pos.y()/CLHEP::mm, pos.z()/CLHEP::mm} );
649 sth.setMomentum( {float(mom.x()/CLHEP::GeV),float(mom.y()/CLHEP::GeV),float(mom.z()/CLHEP::GeV)} );
650 auto particleIt = pm->
particles().find(trackID);
651 if( ( particleIt != pm->
particles().end()) ){
655 sth.setProducedBySecondary( (particleIt->second->originalG4ID != t.
trackID) );
665 for(
unsigned i=0 ; i < nhits ; ++i){
666 auto sch = hits.first.create();
669 sch.setCellID( hit->
cellID );
670 sch.setPosition({float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm)});
675 for(
auto ci=hit->
truth.begin(); ci != hit->
truth.end(); ++ci){
677 auto sCaloHitCont = hits.second.create();
678 sch.addToContributions( sCaloHitCont );
683 sCaloHitCont.setEnergy( c.
deposit/CLHEP::GeV );
684 sCaloHitCont.setTime( c.
time/CLHEP::ns );
685 sCaloHitCont.setParticle( mcp );
688 edm4hep::Vector3f p(c.
x/CLHEP::mm, c.
y/CLHEP::mm, c.
z/CLHEP::mm);
689 sCaloHitCont.setPDG( c.
pdgID );
690 sCaloHitCont.setStepPosition( p );
691 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 3)
692 sCaloHitCont.setStepLength(c.
length / CLHEP::mm);
699 error(
"+++ unknown type in Geant4HitCollection %s ", coll->
type().type().name());