diff --git a/DDG4/hepmc/HepMC3EventReader.cpp b/DDG4/hepmc/HepMC3EventReader.cpp index 4f993de5880c7347538e3b5deadabe3bf058c948..790c78cc79a251a2fdc055819e6eb03fe483a2d4 100644 --- a/DDG4/hepmc/HepMC3EventReader.cpp +++ b/DDG4/hepmc/HepMC3EventReader.cpp @@ -27,6 +27,7 @@ #include "G4Event.hh" + using dd4hep::sim::HEPMC3EventReader; using SGenParticle = std::shared_ptr<HepMC3::GenParticle>; using PropertyMask = dd4hep::detail::ReferenceBitMask<int>; @@ -65,6 +66,19 @@ HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles mcpcoll[i] = p; } + //treat event attributes, flow[12] + std::vector<std::map<int, int>> colorFlow(2, std::map<int, int>()); + std::map<std::string, std::string> eventAttributes {}; + for(auto const& attr: genEvent.attributes()){ + for(auto const& inAttr: attr.second){ + if(attr.first == m_flow1){ + colorFlow[0][inAttr.first] = std::atoi(inAttr.second->unparsed_string().c_str()); + } else if(attr.first == m_flow1){ + colorFlow[1][inAttr.first] = std::atoi(inAttr.second->unparsed_string().c_str()); + } + } + } + double mom_unit = (genEvent.momentum_unit() == HepMC3::Units::MomentumUnit::GEV) ? CLHEP::GeV : CLHEP::MeV; double len_unit = (genEvent.length_unit() == HepMC3::Units::LengthUnit::MM) ? CLHEP::mm : CLHEP::cm; // build collection of MCParticles @@ -72,11 +86,10 @@ HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles auto const& mcp = mcpcoll[i]; Geant4ParticleHandle p(new Particle(i)); auto const& mom = mcp->momentum(); - auto vsx = mcp->production_vertex()->position(); - auto vex = (mcp->end_vertex()) ? mcp->end_vertex()->position() : HepMC3::FourVector(); - + auto const& vsx = mcp->production_vertex()->position(); + auto const& vex = (mcp->end_vertex()) ? mcp->end_vertex()->position() : HepMC3::FourVector(); const float spin[] = {0.0, 0.0, 0.0}; //mcp->getSpin(); // FIXME - const int color[] = {0, 0}; //mcp->getColorFlow(); // FIXME + const int color[] = {colorFlow[0][mcp->id()], colorFlow[1][mcp->id()]}; const int pdg = mcp->pid(); PropertyMask status(p->status); p->pdgID = pdg; @@ -84,8 +97,8 @@ HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles p->psx = mom.get_component(0) * mom_unit; p->psy = mom.get_component(1) * mom_unit; p->psz = mom.get_component(2) * mom_unit; - p->time = vsx.get_component(3) * len_unit / CLHEP::c_light; // FIXME - p->properTime = vsx.get_component(3) * len_unit / CLHEP::c_light; // FIXME + p->time = vsx.get_component(3) * len_unit / CLHEP::c_light; + p->properTime = vsx.get_component(3) * len_unit / CLHEP::c_light; p->vsx = vsx.get_component(0) * len_unit; p->vsy = vsx.get_component(1) * len_unit; p->vsz = vsx.get_component(2) * len_unit; @@ -97,8 +110,8 @@ HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles p->spin[1] = spin[1]; p->spin[2] = spin[2]; p->colorFlow[0] = color[0]; - p->colorFlow[0] = color[1]; - p->mass = mcp->generated_mass()*mom_unit; // FIXME + p->colorFlow[1] = color[1]; + p->mass = mcp->generated_mass() * mom_unit; auto const &par = mcp->parents(), &dau=mcp->children(); for(int num=dau.size(), k=0; k<num; ++k) p->daughters.insert(GET_ENTRY(mcparts,dau[k])); @@ -109,7 +122,6 @@ HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY); else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE); else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED); - else if ( genStatus == 11 ) status.set(G4PARTICLE_GEN_DECAYED); // FIXME else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION); else if ( genStatus == 4 ) status.set(G4PARTICLE_GEN_BEAM); else diff --git a/DDG4/hepmc/HepMC3EventReader.h b/DDG4/hepmc/HepMC3EventReader.h index cb1f04b03ccc1f39ef8c3f3aed17d159f5ae979d..da6f3a514a325773f573e839cd4af02eebb7025b 100644 --- a/DDG4/hepmc/HepMC3EventReader.h +++ b/DDG4/hepmc/HepMC3EventReader.h @@ -15,7 +15,7 @@ // Framework include files #include "DDG4/Geant4InputAction.h" -namespace HepMC3{ struct GenEvent; } +namespace HepMC3{ class GenEvent; } /// Namespace for the AIDA detector description toolkit namespace dd4hep { @@ -27,20 +27,27 @@ namespace dd4hep { /** * \version 1.0 * \ingroup DD4HEP_SIMULATION + * @property: Parameters.Flow1 + * @property: Parameters.Flow2 */ class HEPMC3EventReader : public Geant4EventReader { public: /// Initializing constructor - HEPMC3EventReader(const std::string& nam); + explicit HEPMC3EventReader(const std::string& fileName); /// Default destructor - virtual ~HEPMC3EventReader(); + virtual ~HEPMC3EventReader() = default; - /// Read an event and fill a vector of MCParticles. - virtual EventReaderStatus readParticles(int event_number, - Vertices& vertices, - std::vector<Particle*>& particles); - /// Read an event and return a MCParticles. + /// Read an event and fill a vector of Particles and vertices + virtual EventReaderStatus readParticles(int event_number, Vertices& vertices, Particles& particles); + /// Read an event virtual EventReaderStatus readGenEvent(int event_number, HepMC3::GenEvent& genEvent) = 0; + + protected: + /// name of the GenEvent Attribute storing the color flow1 + std::string m_flow1 = "flow1"; + /// name of the GenEvent Attribute storing the color flow2 + std::string m_flow2 = "flow2"; + }; } /* End namespace sim */ diff --git a/DDG4/hepmc/HepMC3FileReader.cpp b/DDG4/hepmc/HepMC3FileReader.cpp index 4c65435941012cc682fc04905079cfbdfe7d0b1c..67fa3fd0fe4d0e19037dce59d3fbf4b58d7b5acd 100644 --- a/DDG4/hepmc/HepMC3FileReader.cpp +++ b/DDG4/hepmc/HepMC3FileReader.cpp @@ -14,7 +14,7 @@ * @{ \package HepMC3FileReader - * \brief Plugin to read HepMC3 Ascii files + * \brief Plugin to read HepMC3 files * And here we can put a longer description, parameters, examples... * @@ -23,6 +23,8 @@ #include "HepMC3EventReader.h" +#include "DDG4/EventParameters.h" + #include <HepMC3/ReaderFactory.h> /// Namespace for the AIDA detector description toolkit @@ -31,6 +33,29 @@ namespace dd4hep { /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit namespace sim { + template <class T=HepMC3::GenEvent> void EventParameters::ingestParameters(T const& genEvent) { + for(auto const& attr: genEvent.attributes()){ + for(auto const& inAttr: attr.second){ + std::stringstream strstr; + strstr << attr.first; + // if more than one entry for given key, or if counter not 0 append "_counterValue" + if(attr.second.size() > 1 or inAttr.first != 0){ + strstr << "_" << inAttr.first; + } + auto attr_as_int = std::dynamic_pointer_cast<HepMC3::IntAttribute>(inAttr.second); + auto attr_as_flt = std::dynamic_pointer_cast<HepMC3::FloatAttribute>(inAttr.second); + if(attr_as_int) { + m_intValues[strstr.str()] = {attr_as_int->value()}; + } else if(attr_as_flt) { + m_fltValues[strstr.str()] = {attr_as_flt->value()}; + } else { // anything else + m_strValues[strstr.str()] = {inAttr.second->unparsed_string()}; + } + } + } + } + + /// Base class to read hepmc3 event files /** * \version 1.0 @@ -95,6 +120,17 @@ HEPMC3FileReader::readGenEvent(int /*event_number*/, HepMC3::GenEvent& genEvent) ++m_currEvent; if (genEvent.particles().size()) { printout(INFO,"HEPMC3FileReader","Read event from file"); + // Create input event parameters context + try { + Geant4Context* ctx = context(); + EventParameters *parameters = new EventParameters(); + parameters->setEventNumber(genEvent.event_number()); + parameters->ingestParameters(genEvent); + ctx->event().addExtension<EventParameters>(parameters); + } + catch(std::exception &) + { + } return EVENT_READER_OK; } return EVENT_READER_EOF; diff --git a/DDG4/include/DDG4/EventParameters.h b/DDG4/include/DDG4/EventParameters.h new file mode 100644 index 0000000000000000000000000000000000000000..8a05e9e2ddce33110550d52f93a5f4d90521577d --- /dev/null +++ b/DDG4/include/DDG4/EventParameters.h @@ -0,0 +1,69 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// +//========================================================================== +#ifndef DD4HEP_DDG4_EVENTPARAMETERS_H +#define DD4HEP_DDG4_EVENTPARAMETERS_H + +#include <map> +#include <string> +#include <vector> + + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim { + + /// Event extension to pass input event data to output event + /** + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class EventParameters { + protected: + int m_runNumber = -1; + int m_eventNumber = -1; + std::map<std::string, std::vector<int>> m_intValues {}; + std::map<std::string, std::vector<float>> m_fltValues {}; + std::map<std::string, std::vector<std::string>> m_strValues {}; + + public: + /// Initializing constructor + EventParameters() = default; + /// Default destructor + ~EventParameters() = default; + + /// Set the event parameters + void setRunNumber(int runNumber); + void setEventNumber(int eventNumber); + /// Get the run number + int runNumber() const { return m_runNumber; } + /// Get the event number + int eventNumber() const { return m_eventNumber; } + + /// Copy the parameters from source + template <class T> void ingestParameters(T const& source); + /// Put parameters into destination + template <class T> void extractParameters(T& destination); + + /// Get the int event parameters + auto const& intParameters() const { return m_intValues; } + /// Get the float event parameters + auto const& fltParameters() const { return m_fltValues; } + /// Get the string event parameters + auto const& strParameters() const { return m_strValues; } + + }; + + } /* End namespace sim */ +} /* End namespace dd4hep */ +#endif /* DD4HEP_DDG4_EVENTPARAMETERS_H */ diff --git a/DDG4/lcio/Geant4Output2LCIO.cpp b/DDG4/lcio/Geant4Output2LCIO.cpp index d0727c29582a5d0a2032ab835bb932725ab1442f..02119ca7482b31adadcd4696144251a388d6a553 100644 --- a/DDG4/lcio/Geant4Output2LCIO.cpp +++ b/DDG4/lcio/Geant4Output2LCIO.cpp @@ -17,7 +17,8 @@ // Framework include files #include "DD4hep/VolumeManager.h" #include "DDG4/Geant4OutputAction.h" -#include "LCIOEventParameters.h" + +#include "DDG4/EventParameters.h" // Geant4 headers #include "G4Threading.hh" #include "G4AutoLock.hh" @@ -41,6 +42,21 @@ namespace dd4hep { /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit namespace sim { + template <class T=lcio::LCEventImpl> void EventParameters::extractParameters(T& event){ + auto& lcparameters = event.parameters(); + event.setRunNumber(this->runNumber()); + event.setEventNumber(this->eventNumber()); + 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; @@ -351,18 +367,16 @@ lcio::LCCollectionVec* Geant4Output2LCIO::saveParticles(Geant4ParticleMap* parti /// Callback to store the Geant4 event void Geant4Output2LCIO::saveEvent(OutputContext<G4Event>& ctxt) { lcio::LCEventImpl* e = context()->event().extension<lcio::LCEventImpl>(); - LCIOEventParameters* parameters = context()->event().extension<LCIOEventParameters>(false); + EventParameters* parameters = context()->event().extension<EventParameters>(false); int runNumber(0), eventNumber(0); const int eventNumberOffset(m_eventNumberOffset > 0 ? m_eventNumberOffset : 0); const int runNumberOffset(m_runNumberOffset > 0 ? m_runNumberOffset : 0); // Get event number, run number and parameters from extension ... - if ( parameters ) { + if (parameters) { runNumber = parameters->runNumber() + runNumberOffset; eventNumber = parameters->eventNumber() + eventNumberOffset; - LCIOEventParameters::copyLCParameters(parameters->eventParameters(),e->parameters()); - } - // ... or from DD4hep framework - else { + parameters->extractParameters(*e); + } else { // ... or from DD4hep framework runNumber = m_runNo + runNumberOffset; eventNumber = ctxt.context->GetEventID() + eventNumberOffset; } diff --git a/DDG4/src/EventParameters.cpp b/DDG4/src/EventParameters.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0596e5a56fb992ec047d3c760f10fe384ba069ba --- /dev/null +++ b/DDG4/src/EventParameters.cpp @@ -0,0 +1,23 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// +//==================================================================== + +#include "DDG4/EventParameters.h" + +#include <map> +#include <string> + +void dd4hep::sim::EventParameters::setRunNumber(int runNumber) { + m_runNumber = runNumber; +} +void dd4hep::sim::EventParameters::setEventNumber(int eventNumber) { + m_eventNumber = eventNumber; +}