From 188c27298ec730e5b866b22cbf1b43c86b8cd274 Mon Sep 17 00:00:00 2001 From: Andre Sailer Date: Tue, 17 Dec 2024 16:49:36 +0100 Subject: [PATCH] ReadEDM4hep: protect against missing g4context, minor cleanups --- DDG4/edm4hep/EDM4hepFileReader.cpp | 355 ++++++++++++++--------------- 1 file changed, 177 insertions(+), 178 deletions(-) diff --git a/DDG4/edm4hep/EDM4hepFileReader.cpp b/DDG4/edm4hep/EDM4hepFileReader.cpp index 2e7da7b72..9e79eb739 100644 --- a/DDG4/edm4hep/EDM4hepFileReader.cpp +++ b/DDG4/edm4hep/EDM4hepFileReader.cpp @@ -43,7 +43,8 @@ namespace dd4hep::sim { // we will not support MCParticles from different collections using MCPARTICLE_MAP=std::map; - /// get the parameters from the GenericParameters of the input EDM4hep frame and store them in the EventParameters extension + /// get the parameters from the GenericParameters of the input EDM4hep frame and store them in the EventParameters + /// extension template void EventParameters::ingestParameters(T const& source) { auto const& intKeys = source.template getKeys(); for(auto const& key: intKeys) { @@ -64,6 +65,8 @@ namespace dd4hep::sim { } } + /// get the parameters from the GenericParameters of the input EDM4hep run frame and store them in the RunParameters + /// extension template void RunParameters::ingestParameters(T const& source) { auto const& intKeys = source.template getKeys(); for(auto const& key: intKeys) { @@ -84,8 +87,7 @@ namespace dd4hep::sim { } } - - /// Base class to read EDM4hep files + /// Class to read EDM4hep files /** * \version 1.0 * \ingroup DD4HEP_SIMULATION @@ -94,213 +96,210 @@ namespace dd4hep::sim { protected: /// Reference to reader object podio::ROOTReader m_reader {}; + /// Name of the MCParticle collection to read std::string m_collectionName; + /// Name of the EventHeader collection to read std::string m_eventHeaderCollectionName; public: /// Initializing constructor EDM4hepFileReader(const std::string& nam); - /// Read an event and fill a vector of MCParticles. + /// Read parameters set for this action virtual EventReaderStatus setParameters(std::map< std::string, std::string >& parameters); /// Read an event and fill a vector of MCParticles. virtual EventReaderStatus readParticles(int event_number, Vertices& vertices, std::vector& particles); - /// Read an event and return a LCCollectionVec of MCParticles. - + /// Call to register the run parameters void registerRunParameters(); }; - - - -/// Initializing constructor -dd4hep::sim::EDM4hepFileReader::EDM4hepFileReader(const std::string& nam) -: Geant4EventReader(nam) -, m_collectionName("MCParticles") -, m_eventHeaderCollectionName("EventHeader") -{ - printout(INFO,"EDM4hepFileReader","Created file reader. Try to open input %s",nam.c_str()); - m_reader.openFile(nam); - m_directAccess = true; -} - - -void EDM4hepFileReader::registerRunParameters() { - try { - auto *parameters = new RunParameters(); - podio::Frame runFrame = m_reader.readEntry("runs", 0); - parameters->ingestParameters(runFrame.getParameters()); - podio::Frame metaFrame = m_reader.readEntry("metadata", 0); - parameters->ingestParameters(metaFrame.getParameters()); - context()->run().addExtension(parameters); - - } catch(std::exception &e) { - printout(ERROR,"EDM4hepFileReader::registerRunParameters","Failed to register run parameters: %s", e.what()); + /// Initializing constructor + dd4hep::sim::EDM4hepFileReader::EDM4hepFileReader(const std::string& nam) + : Geant4EventReader(nam) + , m_collectionName("MCParticles") + , m_eventHeaderCollectionName("EventHeader") + { + printout(INFO,"EDM4hepFileReader","Created file reader. Try to open input %s",nam.c_str()); + m_reader.openFile(nam); + m_directAccess = true; } -} + void EDM4hepFileReader::registerRunParameters() { + try { + auto *parameters = new RunParameters(); + podio::Frame runFrame = m_reader.readEntry("runs", 0); + parameters->ingestParameters(runFrame.getParameters()); + podio::Frame metaFrame = m_reader.readEntry("metadata", 0); + parameters->ingestParameters(metaFrame.getParameters()); + context()->run().addExtension(parameters); + } catch(std::exception &e) { + printout(ERROR,"EDM4hepFileReader::registerRunParameters","Failed to register run parameters: %s", e.what()); + } + } -namespace { - /// Helper function to look up MCParticles from mapping - inline int GET_ENTRY(MCPARTICLE_MAP const& mcparts, int partID) { - MCPARTICLE_MAP::const_iterator ip = mcparts.find(partID); - if ( ip == mcparts.end() ) { - throw std::runtime_error("Unknown particle identifier look-up!"); + namespace { + /// Helper function to look up MCParticles from mapping + inline int GET_ENTRY(MCPARTICLE_MAP const& mcparts, int partID) { + MCPARTICLE_MAP::const_iterator ip = mcparts.find(partID); + if ( ip == mcparts.end() ) { + throw std::runtime_error("Unknown particle identifier look-up!"); + } + return (*ip).second; } - return (*ip).second; } -} - - -/// Read an event and fill a vector of MCParticles. -EDM4hepFileReader::EventReaderStatus -EDM4hepFileReader::readParticles(int event_number, Vertices& vertices, std::vector& particles) { - - podio::Frame frame = m_reader.readEntry("events", event_number); - const auto& primaries = frame.get(m_collectionName); - if (primaries.isValid()) { - //Read the event header collection and add it to the context as an extension - const auto& eventHeaderCollection = frame.get(m_eventHeaderCollectionName); - if(eventHeaderCollection.isValid() && eventHeaderCollection.size() == 1){ - const auto& eh = eventHeaderCollection.at(0); - auto eventNumber = eh.getEventNumber(); - auto runNumber = eh.getRunNumber(); - Geant4Context* ctx = context(); - auto clonedEh = new edm4hep::MutableEventHeader(eh.clone()); - ctx->event().addExtension(clonedEh); + /// Read an event and fill a vector of MCParticles + EDM4hepFileReader::EventReaderStatus + EDM4hepFileReader::readParticles(int event_number, Vertices& vertices, std::vector& particles) { + m_currEvent = event_number; + podio::Frame frame = m_reader.readEntry("events", event_number); + const auto& primaries = frame.get(m_collectionName); + if (primaries.isValid()) { printout(INFO,"EDM4hepFileReader","read collection %s from event %d in run %d ", m_collectionName.c_str(), eventNumber, runNumber); - // Create input event parameters context - try { - EventParameters *parameters = new EventParameters(); - parameters->setRunNumber(runNumber); - parameters->setEventNumber(eventNumber); - parameters->ingestParameters(frame.getParameters()); - ctx->event().addExtension(parameters); - } catch(std::exception &) {} + //Read the event header collection and add it to the context as an extension + const auto& eventHeaderCollection = frame.get(m_eventHeaderCollectionName); + if(eventHeaderCollection.isValid() && eventHeaderCollection.size() == 1){ + const auto& eh = eventHeaderCollection.at(0); + auto eventNumber = eh.getEventNumber(); + auto runNumber = eh.getRunNumber(); + try { + Geant4Context* ctx = context(); + ctx->event().addExtension(new edm4hep::MutableEventHeader(eh.clone())); + } catch(std::exception& ) {} + // Create input event parameters context + try { + Geant4Context* ctx = context(); + EventParameters *parameters = new EventParameters(); + parameters->setRunNumber(runNumber); + parameters->setEventNumber(eventNumber); + parameters->ingestParameters(frame.getParameters()); + ctx->event().addExtension(parameters); + } catch(std::exception &) {} + } else { + printout(WARNING,"EDM4hepFileReader","No EventHeader collection found or too many event headers!"); + try { + Geant4Context* ctx = context(); + EventParameters *parameters = new EventParameters(); + parameters->setRunNumber(0); + parameters->setEventNumber(event_number); + parameters->ingestParameters(frame.getParameters()); + ctx->event().addExtension(parameters); + } catch(std::exception &) {} + } } else { - printout(WARNING,"EDM4hepFileReader","No EventHeader collection found"); - try { - Geant4Context* ctx = context(); - EventParameters *parameters = new EventParameters(); - parameters->setRunNumber(0); - parameters->setEventNumber(event_number); - parameters->ingestParameters(frame.getParameters()); - ctx->event().addExtension(parameters); - } catch(std::exception &) {} + return EVENT_READER_EOF; } - } else { - return EVENT_READER_EOF; - } - - printout(INFO,"EDM4hepFileReader", "We read the particle collection"); - int NHEP = primaries.size(); - // check if there is at least one particle - if ( NHEP == 0 ) { - printout(WARNING,"EDM4hepFileReader", "We have no particles"); - return EVENT_READER_NO_PRIMARIES; - } - - MCPARTICLE_MAP mcparts; - std::vector mcpcoll; - mcpcoll.resize(NHEP); - for(int i=0; ipdgID = pdg; - p->charge = int(mcp.getCharge()*3.0); - p->psx = mom[0]*CLHEP::GeV; - p->psy = mom[1]*CLHEP::GeV; - p->psz = mom[2]*CLHEP::GeV; - p->time = mcp.getTime()*CLHEP::ns; - p->properTime = mcp.getTime()*CLHEP::ns; - p->vsx = vsx[0]*CLHEP::mm; - p->vsy = vsx[1]*CLHEP::mm; - p->vsz = vsx[2]*CLHEP::mm; - p->vex = vex[0]*CLHEP::mm; - p->vey = vex[1]*CLHEP::mm; - p->vez = vex[2]*CLHEP::mm; - p->process = 0; - p->spin[0] = spin[0]; - p->spin[1] = spin[1]; - p->spin[2] = spin[2]; - p->colorFlow[0] = color[0]; - p->colorFlow[1] = color[1]; - p->mass = mcp.getMass()*CLHEP::GeV; - const auto par = mcp.getParents(), &dau=mcp.getDaughters(); - for(int num=dau.size(),k=0; kdaughters.insert(GET_ENTRY(mcparts,dau_k.getObjectID().index)); - } - for(int num=par.size(),k=0; kparents.insert(GET_ENTRY(mcparts, par_k.getObjectID().index)); + printout(INFO,"EDM4hepFileReader", "We read the particle collection"); + int NHEP = primaries.size(); + // check if there is at least one particle + if ( NHEP == 0 ) { + printout(WARNING,"EDM4hepFileReader", "We have no particles"); + return EVENT_READER_NO_PRIMARIES; } - PropertyMask status(p->status); - int genStatus = mcp.getGeneratorStatus(); - // Copy raw generator status - p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK; - m_inputAction->setGeneratorStatus(genStatus, status); - - //fg: we simply add all particles without parents as with their own vertex. - // This might include the incoming beam particles, e.g. in - // the case of lcio files written with Whizard2, which is slightly odd, - // however should be treated correctly in Geant4InputHandling.cpp. - // We no longer make an attempt to identify the incoming particles - // based on the generator status, as this varies widely with different - // generators. - - if ( p->parents.size() == 0 ) { - - Geant4Vertex* vtx = new Geant4Vertex ; - vertices.emplace_back( vtx ); - vtx->x = p->vsx; - vtx->y = p->vsy; - vtx->z = p->vsz; - vtx->time = p->time; - - vtx->out.insert(p->id) ; + MCPARTICLE_MAP mcparts; + std::vector mcpcoll; + mcpcoll.resize(NHEP); + for(int i=0; ipdgID = pdg; + p->charge = int(mcp.getCharge()*3.0); + p->psx = mom[0]*CLHEP::GeV; + p->psy = mom[1]*CLHEP::GeV; + p->psz = mom[2]*CLHEP::GeV; + p->time = mcp.getTime()*CLHEP::ns; + p->properTime = mcp.getTime()*CLHEP::ns; + p->vsx = vsx[0]*CLHEP::mm; + p->vsy = vsx[1]*CLHEP::mm; + p->vsz = vsx[2]*CLHEP::mm; + p->vex = vex[0]*CLHEP::mm; + p->vey = vex[1]*CLHEP::mm; + p->vez = vex[2]*CLHEP::mm; + p->process = 0; + p->spin[0] = spin[0]; + p->spin[1] = spin[1]; + p->spin[2] = spin[2]; + p->colorFlow[0] = color[0]; + p->colorFlow[1] = color[1]; + p->mass = mcp.getMass()*CLHEP::GeV; + const auto par = mcp.getParents(), &dau=mcp.getDaughters(); + for(int num=dau.size(),k=0; kdaughters.insert(GET_ENTRY(mcparts,dau_k.getObjectID().index)); + } + for(int num=par.size(),k=0; kparents.insert(GET_ENTRY(mcparts, par_k.getObjectID().index)); + } + + PropertyMask status(p->status); + int genStatus = mcp.getGeneratorStatus(); + // Copy raw generator status + p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK; + if(m_inputAction) { + // in some tests we do not set up the inputAction + m_inputAction->setGeneratorStatus(genStatus, status); + } + + //fg: we simply add all particles without parents as with their own vertex. + // This might include the incoming beam particles, e.g. in + // the case of lcio files written with Whizard2, which is slightly odd, + // however should be treated correctly in Geant4InputHandling.cpp. + // We no longer make an attempt to identify the incoming particles + // based on the generator status, as this varies widely with different + // generators. + + if ( p->parents.size() == 0 ) { + + Geant4Vertex* vtx = new Geant4Vertex ; + vertices.emplace_back( vtx ); + vtx->x = p->vsx; + vtx->y = p->vsy; + vtx->z = p->vsz; + vtx->time = p->time; + + vtx->out.insert(p->id) ; + } + + if ( mcp.isCreatedInSimulation() ) status.set(G4PARTICLE_SIM_CREATED); + if ( mcp.isBackscatter() ) status.set(G4PARTICLE_SIM_BACKSCATTER); + if ( mcp.vertexIsNotEndpointOfParent() ) status.set(G4PARTICLE_SIM_PARENT_RADIATED); + if ( mcp.isDecayedInTracker() ) status.set(G4PARTICLE_SIM_DECAY_TRACKER); + if ( mcp.isDecayedInCalorimeter() ) status.set(G4PARTICLE_SIM_DECAY_CALO); + if ( mcp.hasLeftDetector() ) status.set(G4PARTICLE_SIM_LEFT_DETECTOR); + if ( mcp.isStopped() ) status.set(G4PARTICLE_SIM_STOPPED); + if ( mcp.isOverlay() ) status.set(G4PARTICLE_SIM_OVERLAY); + particles.emplace_back(p); + } + return EVENT_READER_OK; } - return EVENT_READER_OK; -} - -/// Set the parameters for the class, in particular the name of the MCParticle -/// list -Geant4EventReader::EventReaderStatus -dd4hep::sim::EDM4hepFileReader::setParameters( std::map< std::string, std::string > & parameters ) { - _getParameterValue( parameters, "MCParticleCollectionName", m_collectionName, m_collectionName); - return EVENT_READER_OK; -} - + /// Set the parameters for the class, in particular the name of the MCParticle + /// list + Geant4EventReader::EventReaderStatus + dd4hep::sim::EDM4hepFileReader::setParameters( std::map< std::string, std::string > & parameters ) { + _getParameterValue(parameters, "MCParticleCollectionName", m_collectionName, m_collectionName); + _getParameterValue(parameters, "EventHeaderCollectionName", m_eventHeaderCollectionName, m_eventHeaderCollectionName); + return EVENT_READER_OK; + } } //end dd4hep::sim