diff --git a/include/HEJ/LesHouchesReader.hh b/include/HEJ/LesHouchesReader.hh index d610484..fffe553 100644 --- a/include/HEJ/LesHouchesReader.hh +++ b/include/HEJ/LesHouchesReader.hh @@ -1,123 +1,93 @@ /** \file * \brief Header file for reading events in the Les Houches Event File format. * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2022 * \copyright GPLv2 or later */ #pragma once #include #include #include "LHEF/LHEF.h" #include "HEJ/EventReader.hh" #include "HEJ/stream.hh" namespace HEJ { //! Class for reading events from a file in the Les Houches Event File format class LesHouchesReader : public EventReader{ public: //! Contruct object reading from the given file - explicit LesHouchesReader(std::string const & filename): - stream_{filename}, - reader_{stream_} - { - // always use the newest LHE version - reader_.heprup.version = LHEF::HEPRUP().version; - } + explicit LesHouchesReader(std::string const & filename); //! Read an event - bool read_event() override { - return reader_.readEvent(); - } + bool read_event() override; //! Access header text std::string const & header() const override { return reader_.headerBlock; } //! Access run information LHEF::HEPRUP const & heprup() const override { return reader_.heprup; } //! Access last read event LHEF::HEPEUP const & hepeup() const override { return reader_.hepeup; } - private: - istream stream_; - LHEF::Reader reader_; - - protected: - //! Underlying reader - LHEF::Reader & reader(){ - return reader_; - } - }; - - /** - * @brief Les Houches Event file reader for LHE files created by Sherpa - * @details In Sherpa the cross section is given by - * sum(weights)/(number of trials). This EventReader converts the - * weights such that cross section=sum(weights) - * @note Reading from a pipe is not possible! - */ - class HEJLHEReader : public LesHouchesReader{ - public: - //! Inialise Reader for a Sherpa LHE file - explicit HEJLHEReader(std::string const & filename); - - bool read_event() override; - std::optional number_events() const override { if (generator_==Generator::Sherpa) { return num_events_; } else { std::size_t start = header().rfind("Number of Events"); start = header().find_first_of("123456789", start); if(start == std::string::npos) { return {}; } const std::size_t end = header().find_first_not_of("0123456789", start); return std::stoi(header().substr(start, end - start)); } } double scalefactor() override { if (generator_==Generator::Sherpa) { - int num_trials=std::accumulate(num_trials_.begin(),num_trials_.end(),0); - return 1./double(num_trials); + const std::size_t num_trials = std::accumulate( + num_trials_.begin(), num_trials_.end(), 0 + ); + return 1./static_cast(num_trials); } else - return 1; + return 1.; } private: enum class Generator{ HEJ, HEJFOG, Sherpa, MG, unknown }; bool calculate_XSECUP() { return calculate_XSECUP_; } Generator get_generator( LHEF::HEPRUP const & heprup, std::string const & header ); + istream stream_; + LHEF::Reader reader_; std::vector num_trials_; size_t num_events_; Generator generator_; bool calculate_XSECUP_; - - }; + }; } // namespace HEJ diff --git a/src/EventReader.cc b/src/EventReader.cc index 5015b74..6dd5189 100644 --- a/src/EventReader.cc +++ b/src/EventReader.cc @@ -1,34 +1,32 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2022 * \copyright GPLv2 or later */ #include "HEJ/EventReader.hh" #include #include #include "LHEF/LHEF.h" #include "HEJ/ConfigFlags.hh" #include "HEJ/HDF5Reader.hh" #include "HEJ/LesHouchesReader.hh" #include "HEJ/Version.hh" namespace HEJ { std::unique_ptr make_reader(std::string const & filename) { try { - auto reader{ std::make_unique(filename) }; + auto reader{ std::make_unique(filename) }; return reader; } catch(std::runtime_error&) { #ifdef HEJ_BUILD_WITH_HDF5 return std::make_unique(filename); #else throw; #endif } } - - - } // namespace HEJ +} // namespace HEJ diff --git a/src/LesHouchesReader.cc b/src/LesHouchesReader.cc index 856a7cc..54c46d7 100644 --- a/src/LesHouchesReader.cc +++ b/src/LesHouchesReader.cc @@ -1,82 +1,85 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2022 * \copyright GPLv2 or later */ #include "HEJ/LesHouchesReader.hh" #include "HEJ/Version.hh" #include namespace HEJ { - HEJLHEReader::Generator HEJLHEReader::get_generator( + LesHouchesReader::Generator LesHouchesReader::get_generator( LHEF::HEPRUP const & heprup, std::string const & header ){ // try getting generator name from specific tag if(!heprup.generators.empty()){ std::string const & gen_name = heprup.generators.back().name; if(gen_name == "HEJ" || gen_name == HEJ::Version::String()) return Generator::HEJ; if(gen_name == "HEJ Fixed Order Generation") return Generator::HEJFOG; if(gen_name == "MadGraph5_aMC@NLO") return Generator::MG; std::cerr << "Unknown LHE Generator " << gen_name << " using default LHE interface.\n"; return Generator::unknown; } // The generator tag is not always used -> check by hand if(header.find("generated with HEJ")!=std::string::npos) return Generator::HEJ; if(header.find("# created by SHERPA")!=std::string::npos) return Generator::Sherpa; if(header.find("")!=std::string::npos) return Generator::MG; std::cerr<<"Could not determine LHE Generator using default LHE interface.\n"; return Generator::unknown; } - HEJLHEReader::HEJLHEReader(std::string const & filename): - LesHouchesReader{filename}, - num_trials_(10,0), + LesHouchesReader::LesHouchesReader( + std::string const & filename + ): + stream_{filename}, + reader_{stream_}, + num_trials_(reader_.heprup.NPRUP, 0), num_events_(0), - calculate_XSECUP_{true} + generator_{get_generator(reader_.heprup, reader_.headerBlock)}, + // calculate cross section if it was not set to > 0. + calculate_XSECUP_{reader_.heprup.XSECUP.front() <= 0.} { - generator_=get_generator(reader().heprup, reader().headerBlock); - int nprup=reader().heprup.NPRUP; - num_trials_.resize(nprup); - reader().heprup.XSECUP.resize(nprup); - // check if XSECUP is already set - if (reader().heprup.XSECUP.front()>0) {// if it is set - calculate_XSECUP_=false; - } + // always use the newest LHE version + reader_.heprup.version = LHEF::HEPRUP().version; + + reader_.heprup.XSECUP.resize(reader_.heprup.NPRUP); // For IDWTUP == 1 or 4 we assume avg(weight)=xs // With the modified weights we have in Sherpa sum(weight)=xs // -> overwrite IDWTUP to "something neutral" - reader().heprup.IDWTUP = reader().heprup.IDWTUP>0?3:-3; + if(generator_ == Generator::Sherpa) { + reader_.heprup.IDWTUP = reader_.heprup.IDWTUP>0?3:-3; + } } - bool HEJLHEReader::read_event() { - if(!LesHouchesReader::read_event()) return false; + bool LesHouchesReader::read_event() { + if(!reader_.readEvent()) return false; - int nprup=reader().heprup.NPRUP; + int nprup=reader_.heprup.NPRUP; int position=std::max(1,std::min(hepeup().IDPRUP,nprup))-1; // 0..nprup-1 if (generator_ == Generator::Sherpa) { ++num_events_; - reader().heprup.XSECUP.at(position) *= num_trials_.at(position); + reader_.heprup.XSECUP.at(position) *= num_trials_.at(position); num_trials_.at(position) += std::stod(hepeup().attributes.at("trials")); - reader().heprup.XSECUP.at(position) += hepeup().XWGTUP; - reader().heprup.XSECUP.at(position) /= + reader_.heprup.XSECUP.at(position) += hepeup().XWGTUP; + reader_.heprup.XSECUP.at(position) /= num_trials_.at(position); } else { ++num_events_; - if (HEJLHEReader::calculate_XSECUP()) { - reader().heprup.XSECUP.at(position) += hepeup().XWGTUP; + if (LesHouchesReader::calculate_XSECUP()) { + reader_.heprup.XSECUP.at(position) += hepeup().XWGTUP; } } return true; } } // namespace HEJ