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 #include <podio/CollectionBase.h>
37 #include <podio/podioVersion.h>
38 #include <podio/Frame.h>
39 #include <podio/FrameCategories.h>
40 #if PODIO_BUILD_VERSION >= PODIO_VERSION(0, 99, 0)
41 #include <podio/ROOTWriter.h>
43 #include <podio/ROOTFrameWriter.h>
45 using ROOTWriter = podio::ROOTFrameWriter;
73 using trackermap_t = std::map< std::string, edm4hep::SimTrackerHitCollection >;
74 using calorimeterpair_t = std::pair< edm4hep::SimCalorimeterHitCollection, edm4hep::CaloHitContributionCollection >;
76 std::unique_ptr<writer_t>
m_file { };
107 virtual void beginRun(
const G4Run* run);
109 virtual void endRun(
const G4Run* run);
112 virtual void saveRun(
const G4Run* run);
114 virtual void saveEvent( OutputContext<G4Event>& ctxt);
118 virtual void commit( OutputContext<G4Event>& ctxt);
121 virtual void begin(
const G4Event* event);
124 template <
typename T>
126 for(
const auto& p : parameters) {
127 std::stringstream output;
128 output <<
"Saving event parameter: "
129 << std::setw(32) << p.first
130 << std::setw(20) << p.second;
131 info(output.str().c_str());
132 m_frame.putParameter(p.first, p.second);
139 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving event parameter: %s", p.first.c_str());
140 frame.putParameter(p.first, p.second);
143 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving event parameter: %s", p.first.c_str());
144 frame.putParameter(p.first, p.second);
147 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving event parameter: %s", p.first.c_str());
148 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);
159 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving run parameter: %s", p.first.c_str());
160 frame.putParameter(p.first, p.second);
163 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving run parameter: %s", p.first.c_str());
164 frame.putParameter(p.first, p.second);
167 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving run parameter: %s", p.first.c_str());
168 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);
178 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving meta parameter: %s", p.first.c_str());
179 frame.putParameter(p.first, p.second);
182 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving meta parameter: %s", p.first.c_str());
183 frame.putParameter(p.first, p.second);
186 printout(DEBUG,
"Geant4OutputEDM4hep",
"Saving meta parameter: %s", p.first.c_str());
187 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);
198 #endif // DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
226 #include <G4Threading.hh>
227 #include <G4AutoLock.hh>
228 #include <G4Version.hh>
229 #include <G4ParticleDefinition.hh>
230 #include <G4VProcess.hh>
231 #include <G4Event.hh>
234 #include <CLHEP/Units/SystemOfUnits.h>
237 #include <edm4hep/EventHeaderCollection.h>
243 G4Mutex action_mutex = G4MUTEX_INITIALIZER;
251 :
Geant4OutputAction(ctxt,nam), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
253 declareProperty(
"RunHeader", m_runHeader);
254 declareProperty(
"EventParametersInt", m_eventParametersInt);
255 declareProperty(
"EventParametersFloat", m_eventParametersFloat);
256 declareProperty(
"EventParametersString", m_eventParametersString);
257 declareProperty(
"RunParametersInt", m_runParametersInt);
258 declareProperty(
"RunParametersFloat", m_runParametersFloat);
259 declareProperty(
"RunParametersString", m_runParametersString);
260 declareProperty(
"RunNumberOffset", m_runNumberOffset);
261 declareProperty(
"EventNumberOffset", m_eventNumberOffset);
262 declareProperty(
"SectionName", m_section_name);
263 declareProperty(
"FilesByRun", m_filesByRun);
264 info(
"Writer is now instantiated ..." );
270 G4AutoLock protection_lock(&action_mutex);
276 G4AutoLock protection_lock(&action_mutex);
280 std::size_t idx =
m_output.rfind(
".");
281 if ( idx != std::string::npos ) {
286 if ( !fname.empty() && !
m_file ) {
287 m_file = std::make_unique<podio::ROOTWriter>(fname);
289 fatal(
"+++ Failed to open output file: %s", fname.c_str());
291 printout( INFO,
"Geant4Output2EDM4hep" ,
"Opened %s for output", fname.c_str() ) ;
304 G4AutoLock protection_lock(&action_mutex);
313 podio::Frame metaFrame{};
315 metaFrame.putParameter(podio::collMetadataParamName(
name, CellIDEncoding), encodingStr);
317 G4AutoLock protection_lock(&action_mutex);
318 m_file->writeFrame(metaFrame,
"metadata");
324 G4AutoLock protection_lock(&action_mutex);
327 m_frame.put( std::move(it->second), it->first);
330 m_frame.put( std::move(calorimeterHits.first), colName);
331 m_frame.put( std::move(calorimeterHits.second), colName +
"Contributions");
340 except(
"+++ Failed to write output file. [Stream is not open]");
345 G4AutoLock protection_lock(&action_mutex);
349 podio::Frame runHeader {};
351 runHeader.putParameter(
key, value);
354 runHeader.putParameter(
key, value);
357 runHeader.putParameter(
key, value);
360 runHeader.putParameter(
key, value);
363 runHeader.putParameter(
"runNumber",
m_runNo);
364 runHeader.putParameter(
"GEANT4Version", G4Version);
366 runHeader.putParameter(
"detectorName",
context()->detectorDescription().header().
name());
369 if (
context()->runPtr() !=
nullptr) {
374 m_file->writeFrame(runHeader,
"runs");
380 podio::Frame metaFrame {};
385 m_file->writeFrame(metaFrame,
"meta");
401 typedef detail::ReferenceBitMask<const int>
PropertyMask;
406 if ( pm.size() > 0 ) {
409 std::map<int,int> p_ids;
410 std::vector<const Geant4Particle*> p_part;
411 p_part.reserve(pm.size());
413 for (
const auto& iParticle : pm) {
414 int id = iParticle.first;
418 const G4ParticleDefinition* def = p.
definition();
420 mcp.setPDG(p->
pdgID);
422 using MT = decltype(std::declval<edm4hep::MCParticle>().getMomentum().x);
423 mcp.setMomentum( {MT(p->
psx/CLHEP::GeV),MT(p->
psy/CLHEP::GeV),MT(p->
psz/CLHEP::GeV)} );
424 mcp.setMomentumAtEndpoint( {MT(p->
pex/CLHEP::GeV),MT(p->
pey/CLHEP::GeV),MT(p->
pez/CLHEP::GeV)} );
426 double vs_fa[3] = { p->
vsx/CLHEP::mm, p->
vsy/CLHEP::mm, p->
vsz/CLHEP::mm } ;
427 mcp.setVertex( vs_fa );
429 double ve_fa[3] = { p->
vex/CLHEP::mm, p->
vey/CLHEP::mm, p->
vez/CLHEP::mm } ;
430 mcp.setEndpoint( ve_fa );
432 mcp.setTime(p->
time/CLHEP::ns);
433 mcp.setMass(p->
mass/CLHEP::GeV);
434 mcp.setCharge(def ? def->GetPDGCharge() : 0);
437 mcp.setGeneratorStatus(0);
456 mcp.setOverlay(
false );
459 if( mcp.isCreatedInSimulation() )
460 mcp.setGeneratorStatus( 0 ) ;
462 #if EDM4HEP_MCPARTICLE_HAS_HELICITY
463 mcp.setHelicity(p->
spin[2]);
465 mcp.setSpin(p->
spin);
473 for(
size_t i=0; i < p_ids.size(); ++i) {
478 const auto k = p_ids.find(idau);
479 if (k == p_ids.end()) {
480 fatal(
"+++ Particle %d: FAILED to find daughter with ID:%d",p->
id,idau);
483 int iqdau = (*k).second;
485 q.addToDaughters(qdau);
488 for (
const auto& ipar : p->
parents) {
490 const auto k = p_ids.find(ipar);
491 if (k == p_ids.end()) {
492 fatal(
"+++ Particle %d: FAILED to find parent with ID:%d",p->
id,ipar);
495 int iqpar = (*k).second;
497 q.addToParents(qpar);
507 int runNumber(0), eventNumber(0);
510 double eventWeight{0};
513 runNumber = parameters->
runNumber() + runNumberOffset;
514 eventNumber = parameters->
eventNumber() + eventNumberOffset;
516 #if PODIO_BUILD_VERSION > PODIO_VERSION(0, 99, 0)
517 eventWeight =
m_frame.getParameter<
double>(
"EventWeights").value_or(0.0);
519 eventWeight =
m_frame.getParameter<
double>(
"EventWeights");
522 runNumber =
m_runNo + runNumberOffset;
523 eventNumber = ctxt.
context->GetEventID() + eventNumberOffset;
525 printout(INFO,
"Geant4Output2EDM4hep",
"+++ Saving EDM4hep event %d run %d.", eventNumber, runNumber);
528 edm4hep::EventHeaderCollection header_collection;
530 auto header = header_collection.create();
531 header.setRunNumber(runNumber);
532 header.setEventNumber(eventNumber);
533 header.setWeight(eventWeight);
535 header.setTimeStamp(std::time(
nullptr));
540 header.setTimeStamp(meh->getTimeStamp());
541 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
542 for (
auto const& weight: meh->getWeights()) {
543 header.addToWeights(weight);
548 m_frame.put(std::move(header_collection),
"EventHeader");
555 print(
"+++ Saving %d EDM4hep particles....",
int(part_map->
particleMap.size()));
572 operator std::string()
const {
573 const auto* sd = m_coll->sensitive();
584 std::string colName = collection->GetName();
585 if( coll ==
nullptr ){
586 error(
" no Geant4HitCollection: %s ", colName.c_str());
589 size_t nhits = collection->GetSize();
591 debug(
"+++ Saving EDM4hep collection %s with %d entries.", colName.c_str(),
int(nhits));
600 for(
unsigned i=0 ; i < nhits ; ++i){
601 auto sth = hits.create();
609 sth.setCellID( hit->
cellID ) ;
611 sth.setPathLength(hit->
length/CLHEP::mm);
613 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 10, 99)
614 sth.setParticle(mcp);
616 sth.setMCParticle(mcp);
618 sth.setPosition( {pos.x()/CLHEP::mm, pos.y()/CLHEP::mm, pos.z()/CLHEP::mm} );
619 sth.setMomentum( {float(mom.x()/CLHEP::GeV),float(mom.y()/CLHEP::GeV),float(mom.z()/CLHEP::GeV)} );
620 auto particleIt = pm->
particles().find(trackID);
621 if( ( particleIt != pm->
particles().end()) ){
625 sth.setProducedBySecondary( (particleIt->second->originalG4ID != t.
trackID) );
635 for(
unsigned i=0 ; i < nhits ; ++i){
636 auto sch = hits.first.create();
639 sch.setCellID( hit->
cellID );
640 sch.setPosition({float(pos.x()/CLHEP::mm), float(pos.y()/CLHEP::mm), float(pos.z()/CLHEP::mm)});
645 for(
auto ci=hit->
truth.begin(); ci != hit->
truth.end(); ++ci){
647 auto sCaloHitCont = hits.second.create();
648 sch.addToContributions( sCaloHitCont );
653 sCaloHitCont.setEnergy( c.
deposit/CLHEP::GeV );
654 sCaloHitCont.setTime( c.
time/CLHEP::ns );
655 sCaloHitCont.setParticle( mcp );
658 edm4hep::Vector3f p(c.
x/CLHEP::mm, c.
y/CLHEP::mm, c.
z/CLHEP::mm);
659 sCaloHitCont.setPDG( c.
pdgID );
660 sCaloHitCont.setStepPosition( p );
661 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 3)
662 sCaloHitCont.setStepLength(c.
length / CLHEP::mm);
669 error(
"+++ unknown type in Geant4HitCollection %s ", coll->
type().type().name());