HepMC3FileReader.cpp

Go to the documentation of this file.

//==========================================================================
//  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 "HepMC3EventReader.h"
 
#include <DDG4/EventParameters.h>
#include <DDG4/RunParameters.h>
 
#include <HepMC3/ReaderFactory.h>
#include <HepMC3/Version.h>
 
namespace dd4hep  {
 
  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;
          }
          if(auto int_attr = std::dynamic_pointer_cast<HepMC3::IntAttribute>(inAttr.second)) {
            m_intValues[strstr.str()] = {int_attr->value()};
          } else if(auto flt_attr = std::dynamic_pointer_cast<HepMC3::FloatAttribute>(inAttr.second)) {
            m_fltValues[strstr.str()] = {flt_attr->value()};
          } else if(auto dbl_attr = std::dynamic_pointer_cast<HepMC3::DoubleAttribute>(inAttr.second)) {
            m_dblValues[strstr.str()] = {dbl_attr->value()};
          } else { // anything else
            m_strValues[strstr.str()] = {inAttr.second->unparsed_string()};
          }
        }
      }
 
      if (const auto& weights = genEvent.weights(); !weights.empty()) {
        m_dblValues["EventWeights"] = weights;
      }
      // Not using the GenEven::weight_names here because that checks for
      // emptyness and throws in that case. We don't care if we get an empty
      // vector at this point
      if (genEvent.run_info() && !genEvent.run_info()->weight_names().empty()) {
        m_strValues["EventWeightNames"] = genEvent.weight_names();
      }
    }
 
    template <class T=HepMC3::GenRunInfo> void RunParameters::ingestParameters(T const& runInfo) {
      // This attributes is not the same return type as for GenEvent!
      for(auto const& attr: runInfo.attributes()){
        if(auto int_attr = std::dynamic_pointer_cast<HepMC3::IntAttribute>(attr.second)) {
          m_intValues[attr.first] = {int_attr->value()};
        } else if(auto flt_attr = std::dynamic_pointer_cast<HepMC3::FloatAttribute>(attr.second)) {
          m_fltValues[attr.first] = {flt_attr->value()};
        } else if(auto dbl_attr = std::dynamic_pointer_cast<HepMC3::DoubleAttribute>(attr.second)) {
          m_dblValues[attr.first] = {dbl_attr->value()};
        } else { // anything else
          m_strValues[attr.first] = {attr.second->unparsed_string()};
        }
      }
    }
 
 
    class HEPMC3FileReader : public HEPMC3EventReader  {
    protected:
      std::shared_ptr<HepMC3::Reader> m_reader;
    public:
      HEPMC3FileReader(const std::string& nam);
      virtual ~HEPMC3FileReader() = default;
 
      virtual EventReaderStatus readGenEvent(int event_number, HepMC3::GenEvent& genEvent);
      virtual EventReaderStatus moveToEvent(int event_number);
      //virtual EventReaderStatus skipEvent() { return EVENT_READER_OK; }
      virtual EventReaderStatus setParameters(std::map< std::string, std::string >& parameters);
      virtual void registerRunParameters();
 
    };
  }
}
 
#include <DD4hep/Printout.h>
#include <DDG4/Factories.h>
 
using dd4hep::sim::HEPMC3FileReader;
using dd4hep::sim::Geant4EventReader;
 
// Factory entry
DECLARE_GEANT4_EVENT_READER_NS(dd4hep::sim,HEPMC3FileReader)
 
HEPMC3FileReader::HEPMC3FileReader(const std::string& nam)
: HEPMC3EventReader(nam)
{
  printout(INFO,"HEPMC3FileReader","Created file reader. Try to open input %s", nam.c_str());
  m_reader = HepMC3::deduce_reader(nam);
#if HEPMC3_VERSION_CODE >= 3002006
  // to get the runInfo in the Ascii reader we have to force HepMC to read the first event
  HepMC3::GenEvent dummy;
  m_reader->read_event(dummy);
  // then we get the run info (shared pointer)
  auto runInfo = m_reader->run_info();
  // and deallocate the reader
  m_reader.reset();
  // so we can open the file again from the start
  m_reader = HepMC3::deduce_reader(nam);
  // and set the run info object now
  m_reader->set_run_info(std::move(runInfo));
#endif
  m_directAccess = false;
}
 
void HEPMC3FileReader::registerRunParameters() {
  try {
    auto *parameters = new RunParameters();
    parameters->ingestParameters(*(m_reader->run_info()));
    context()->run().addExtension<RunParameters>(parameters);
  } catch(std::exception &e) {
    printout(ERROR,"HEPMC3FileReader::registerRunParameters","Failed to register run parameters: %s", e.what());
  }
}
 
Geant4EventReader::EventReaderStatus
HEPMC3FileReader::moveToEvent(int event_number) {
  printout(INFO,"HEPMC3FileReader::moveToEvent","Skipping the first %d events ", event_number);
  while( m_currEvent != event_number) {
    printout(INFO,"HEPMC3FileReader::moveToEvent","Event number before skipping: %d", m_currEvent );
    HepMC3::GenEvent genEvent;
    auto status_OK = m_reader->skip(event_number);
    if(not status_OK) {
      return EVENT_READER_IO_ERROR;
    }
    m_currEvent = event_number;
    printout(INFO,"HEPMC3FileReader::moveToEvent","Event number after skipping: %d", m_currEvent );
  }
  return EVENT_READER_OK;
}
 
Geant4EventReader::EventReaderStatus
HEPMC3FileReader::readGenEvent(int /*event_number*/, HepMC3::GenEvent& genEvent)  {
  auto status_OK = m_reader->read_event(genEvent);
    if(not status_OK) {
      return EVENT_READER_IO_ERROR;
    }
  ++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;
}
 
Geant4EventReader::EventReaderStatus
HEPMC3FileReader::setParameters( std::map< std::string, std::string > & parameters ) {
  _getParameterValue(parameters, "Flow1", m_flow1, std::string("flow1"));
  _getParameterValue(parameters, "Flow2", m_flow2, std::string("flow2"));
  return EVENT_READER_OK;
}