diff --git a/include/HEJ/EventReader.hh b/include/HEJ/EventReader.hh index 8f44b0f..f92f26a 100644 --- a/include/HEJ/EventReader.hh +++ b/include/HEJ/EventReader.hh @@ -1,73 +1,63 @@ /** \file * \brief Header file for event reader interface * * This header defines an abstract base class for reading events from files. * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2022 * \copyright GPLv2 or later */ #pragma once #include #include #include #include namespace LHEF { class HEPEUP; class HEPRUP; } namespace HEJ { - enum class generator{ - HEJ, - HEJFOG, - Sherpa, - MG, - unknown - }; - //! Abstract base class for reading events from files struct EventReader { - generator Generator = generator::unknown; - //! Read an event virtual bool read_event() = 0; //! Access header text virtual std::string const & header() const = 0; //! Access run information virtual LHEF::HEPRUP const & heprup() const = 0; //! Access last read event virtual LHEF::HEPEUP const & hepeup() const = 0; //! Guess number of events from header virtual std::optional number_events() const { 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)); } virtual double scalefactor() { return 1.; }; virtual ~EventReader() = default; }; //! Factory function for event readers /** * @param filename The name of the input file * @returns A pointer to an instance of an EventReader * for the input file */ std::unique_ptr make_reader(std::string const & filename); } // namespace HEJ diff --git a/include/HEJ/LesHouchesReader.hh b/include/HEJ/LesHouchesReader.hh index 45fc8f0..5b4197c 100644 --- a/include/HEJ/LesHouchesReader.hh +++ b/include/HEJ/LesHouchesReader.hh @@ -1,112 +1,123 @@ /** \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 * \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; } //! Read an event bool read_event() override { return reader_.readEvent(); } //! 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==HEJ::generator::Sherpa) { + if (Generator==generator::Sherpa) { int num_trials=std::accumulate(num_trials_.begin(),num_trials_.end(),0); return 1./double(num_trials); } else 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 + ); std::vector num_trials_; size_t num_events_; generator Generator; bool calculate_XSECUP_; }; } // namespace HEJ diff --git a/src/LesHouchesReader.cc b/src/LesHouchesReader.cc index 87abbc1..1c56832 100644 --- a/src/LesHouchesReader.cc +++ b/src/LesHouchesReader.cc @@ -1,82 +1,82 @@ /** * \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 { - HEJ::generator get_generator( + HEJLHEReader::generator HEJLHEReader::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 HEJ::generator::HEJ; + return generator::HEJ; if(gen_name == "HEJ Fixed Order Generation") - return HEJ::generator::HEJFOG; - if(gen_name == "MadGraph5_aMC@NLO") return HEJ::generator::MG; + 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 HEJ::generator::unknown; + return generator::unknown; } // The generator tag is not always used -> check by hand if(header.find("generated with HEJ")!=std::string::npos) - return HEJ::generator::HEJ; + return generator::HEJ; if(header.find("# created by SHERPA")!=std::string::npos) - return HEJ::generator::Sherpa; + return generator::Sherpa; if(header.find("")!=std::string::npos) - return HEJ::generator::MG; + return generator::MG; std::cerr<<"Could not determine LHE Generator using default LHE interface.\n"; - return HEJ::generator::unknown; + return generator::unknown; } HEJLHEReader::HEJLHEReader(std::string const & filename): LesHouchesReader{filename}, num_trials_(10,0), num_events_(0), calculate_XSECUP_{true} { 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; } // 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; } bool HEJLHEReader::read_event() { if(!LesHouchesReader::read_event()) return false; int nprup=reader().heprup.NPRUP; int position=std::max(1,std::min(hepeup().IDPRUP,nprup))-1; // 0..nprup-1 - if (Generator == HEJ::generator::Sherpa) { + if (Generator == generator::Sherpa) { ++num_events_; 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) /= num_trials_.at(position); } else { ++num_events_; if (HEJLHEReader::calculate_XSECUP()) { reader().heprup.XSECUP.at(position) += hepeup().XWGTUP; } } return true; } } // namespace HEJ