DDG4/hepmc/HepMC3FileReader.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 //
0011 //==========================================================================
0012
0013 /** \addtogroup Geant4EventReader
0014  *
0015  @{
0016   \package HepMC3FileReader
0017  * \brief Plugin to read HepMC3 files
0018  *
0019  * The HepMC3 File reader uses the HepMC3 Reader Factory to open any HepMC3 compatible file.
0020  *
0021  * However, it was only tested on HepMC3 ascii files produced by Whizard 2.8.2, so please use it carefully, and validate
0022  * its conversion carefully before extensive usage.
0023  *
0024 @}
0025  */
0026
0027 #include "HepMC3EventReader.h"
0028
0029 #include <DDG4/EventParameters.h>
0030 #include <DDG4/RunParameters.h>
0031
0032 #include <HepMC3/ReaderFactory.h>
0033 #include <HepMC3/Version.h>
0034
0035 /// Namespace for the AIDA detector description toolkit
0036 namespace dd4hep  {
0037
0038   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0039   namespace sim  {
0040
0041     template <class T=HepMC3::GenEvent> void EventParameters::ingestParameters(T const& genEvent) {
0042       for(auto const& attr: genEvent.attributes()){
0043         for(auto const& inAttr: attr.second){
0044           std::stringstream strstr;
0045           strstr << attr.first;
0046           // if more than one entry for given key, or if counter not 0 append "_counterValue"
0047           if(attr.second.size() > 1 or inAttr.first != 0){
0048             strstr << "_" << inAttr.first;
0049           }
0050           if(auto int_attr = std::dynamic_pointer_cast<HepMC3::IntAttribute>(inAttr.second)) {
0051             m_intValues[strstr.str()] = {int_attr->value()};
0052           } else if(auto flt_attr = std::dynamic_pointer_cast<HepMC3::FloatAttribute>(inAttr.second)) {
0053             m_fltValues[strstr.str()] = {flt_attr->value()};
0054           } else if(auto dbl_attr = std::dynamic_pointer_cast<HepMC3::DoubleAttribute>(inAttr.second)) {
0055             m_dblValues[strstr.str()] = {dbl_attr->value()};
0056           } else { // anything else
0057             m_strValues[strstr.str()] = {inAttr.second->unparsed_string()};
0058           }
0059         }
0060       }
0061
0062       if (const auto& weights = genEvent.weights(); !weights.empty()) {
0063         m_dblValues["EventWeights"] = weights;
0064       }
0065       // Not using the GenEven::weight_names here because that checks for
0066       // emptyness and throws in that case. We don't care if we get an empty
0067       // vector at this point
0068       if (genEvent.run_info() && !genEvent.run_info()->weight_names().empty()) {
0069         m_strValues["EventWeightNames"] = genEvent.weight_names();
0070       }
0071     }
0072
0073     template <class T=HepMC3::GenRunInfo> void RunParameters::ingestParameters(T const& runInfo) {
0074       // This attributes is not the same return type as for GenEvent!
0075       for(auto const& attr: runInfo.attributes()){
0076         if(auto int_attr = std::dynamic_pointer_cast<HepMC3::IntAttribute>(attr.second)) {
0077           m_intValues[attr.first] = {int_attr->value()};
0078         } else if(auto flt_attr = std::dynamic_pointer_cast<HepMC3::FloatAttribute>(attr.second)) {
0079           m_fltValues[attr.first] = {flt_attr->value()};
0080         } else if(auto dbl_attr = std::dynamic_pointer_cast<HepMC3::DoubleAttribute>(attr.second)) {
0081           m_dblValues[attr.first] = {dbl_attr->value()};
0082         } else { // anything else
0083           m_strValues[attr.first] = {attr.second->unparsed_string()};
0084         }
0085       }
0086     }
0087
0088     /// Base class to read hepmc3 event files
0089     /**
0090      *  \version 1.0
0091      *  \ingroup DD4HEP_SIMULATION
0092      */
0093     class HEPMC3FileReader : public HEPMC3EventReader  {
0094     protected:
0095       /// Reference to reader object
0096       std::shared_ptr<HepMC3::Reader> m_reader;
0097     public:
0098       /// Initializing constructor
0099       HEPMC3FileReader(const std::string& nam);
0100       /// Default destructor
0101       virtual ~HEPMC3FileReader() = default;
0102
0103       /// Read an event and fill a vector of MCParticles.
0104       virtual EventReaderStatus readGenEvent(int event_number, HepMC3::GenEvent& genEvent);
0105       /// skip to event with event_number
0106       virtual EventReaderStatus moveToEvent(int event_number);
0107       //virtual EventReaderStatus skipEvent() { return EVENT_READER_OK; }
0108       virtual EventReaderStatus setParameters(std::map< std::string, std::string >& parameters);
0109       /// register the run parameters into an extension for the run context
0110       virtual void registerRunParameters();
0111
0112     };
0113   }
0114 }
0115
0116 #include <DD4hep/Printout.h>
0117 #include <DDG4/Factories.h>
0118
0119 using dd4hep::sim::HEPMC3FileReader;
0120 using dd4hep::sim::Geant4EventReader;
0121
0122 // Factory entry
0123 DECLARE_GEANT4_EVENT_READER_NS(dd4hep::sim,HEPMC3FileReader)
0124
0125 /// Initializing constructor
0126 HEPMC3FileReader::HEPMC3FileReader(const std::string& nam)
0127 : HEPMC3EventReader(nam)
0128 {
0129   printout(INFO,"HEPMC3FileReader","Created file reader. Try to open input %s", nam.c_str());
0130   m_reader = HepMC3::deduce_reader(nam);
0131 #if HEPMC3_VERSION_CODE >= 3002006
0132   // to get the runInfo in the Ascii reader we have to force HepMC to read the first event
0133   HepMC3::GenEvent dummy;
0134   m_reader->read_event(dummy);
0135   // then we get the run info (shared pointer)
0136   auto runInfo = m_reader->run_info();
0137   // and deallocate the reader
0138   m_reader.reset();
0139   // so we can open the file again from the start
0140   m_reader = HepMC3::deduce_reader(nam);
0141   // and set the run info object now
0142   m_reader->set_run_info(std::move(runInfo));
0143 #endif
0144   m_directAccess = false;
0145 }
0146
0147 void HEPMC3FileReader::registerRunParameters() {
0148   try {
0149     // get RunParameters or create new if not existent yet
0150     auto *parameters = context()->run().extension<RunParameters>(false);
0151     if (!parameters) {
0152       parameters = new RunParameters();
0153       context()->run().addExtension<RunParameters>(parameters);
0154     }
0155     parameters->ingestParameters(*(m_reader->run_info()));
0156   } catch(std::exception &e) {
0157     printout(ERROR,"HEPMC3FileReader::registerRunParameters","Failed to register run parameters: %s", e.what());
0158   }
0159 }
0160
0161 /// moveToSpecifiedEvent, a.k.a. skipNEvents
0162 Geant4EventReader::EventReaderStatus
0163 HEPMC3FileReader::moveToEvent(int event_number) {
0164   printout(INFO,"HEPMC3FileReader::moveToEvent","Skipping the first %d events ", event_number);
0165   while( m_currEvent != event_number) {
0166     printout(INFO,"HEPMC3FileReader::moveToEvent","Event number before skipping: %d", m_currEvent );
0167     HepMC3::GenEvent genEvent;
0168     auto status_OK = m_reader->skip(event_number);
0169     if(not status_OK) {
0170       return EVENT_READER_IO_ERROR;
0171     }
0172     m_currEvent = event_number;
0173     printout(INFO,"HEPMC3FileReader::moveToEvent","Event number after skipping: %d", m_currEvent );
0174   }
0175   return EVENT_READER_OK;
0176 }
0177
0178 /// Read an event and fill a vector of MCParticles.
0179 Geant4EventReader::EventReaderStatus
0180 HEPMC3FileReader::readGenEvent(int /*event_number*/, HepMC3::GenEvent& genEvent)  {
0181   auto status_OK = m_reader->read_event(genEvent);
0182     if(not status_OK) {
0183       return EVENT_READER_IO_ERROR;
0184     }
0185   ++m_currEvent;
0186   if (genEvent.particles().size()) {
0187     printout(INFO,"HEPMC3FileReader","Read event from file");
0188     // Create input event parameters context
0189     try {
0190       // get EventParameters or create new if not existent yet
0191       Geant4Context* ctx = context();
0192       auto* parameters = ctx->event().extension<EventParameters>(false);
0193       if (!parameters) {
0194         parameters = new EventParameters();
0195         ctx->event().addExtension<EventParameters>(parameters);
0196       }
0197       parameters->setEventNumber(genEvent.event_number());
0198       parameters->ingestParameters(genEvent);
0199     }
0200     catch(std::exception &)
0201     {
0202     }
0203     return EVENT_READER_OK;
0204   }
0205   return EVENT_READER_EOF;
0206 }
0207
0208 /// Set the parameters for the class
0209 Geant4EventReader::EventReaderStatus
0210 HEPMC3FileReader::setParameters( std::map< std::string, std::string > & parameters ) {
0211   _getParameterValue(parameters, "Flow1", m_flow1, std::string("flow1"));
0212   _getParameterValue(parameters, "Flow2", m_flow2, std::string("flow2"));
0213   return EVENT_READER_OK;
0214 }