diff --git a/Changes-API.md b/Changes-API.md index 9f1ad56..921cc3c 100644 --- a/Changes-API.md +++ b/Changes-API.md @@ -1,133 +1,134 @@ # Changelog for HEJ API This log lists only changes on the HEJ API. These are primarily code changes relevant for calling HEJ as an API. This file should only be read as an addition to [`Changes.md`](Changes.md), where the main features are documented. ## Version 2.2 ### 2.2.0 #### New and updated functions * New functions to decide whether particles are massive, charged, charged leptons, antiparticles. * Added `to_string` function for `Event`. * Updated function `Event::EventData.reconstruct_intermediate()` - Now requires an `EWConstants` argument - Added support for WpWp and WmWm events. - In the case of WW same-flavour the reconstruction will minimise the total off-shell momentum. * Replaced `no_2_jets` and `bad_final_state` enumerators in `EventType` with `invalid` and `unknown`. * Added helper functions `is_backward_g_to_h`, `is_forward_g_to_h` to detect incoming gluon to outgoing Higgs transitions. * Updated function `implemented_types()` to now be in HEJ namespace. - +* `EventWriter` has a new `set_xs_scale` member function to set the + ratio between the sum of event weights and the cross section. ## Version 2.1 ### 2.1.0 #### Changes to Event class * Restructured `Event` class - `Event` can now only be build from a (new) `Event::EventData` class - Removed default constructor for `Event` - `Event::EventData` replaces the old `UnclusteredEvent` struct. - `UnclusteredEvent` is now **deprecated**, and will be removed in version 2.2.0 - Removed `Event.unclustered()` function - Added new member function `Events.parameters()`, to directly access (underlying) `Parameters` - New member functions `begin_partons`, `end_partons` (`rbegin_partons`, `rend_partons`) with aliases `cbegin_partons`, `cend_partons` (`crbegin_partons`, `crend_partons`) for constant (reversed) iterators over outgoing partons. * New function `Event::EventData.reconstruct_intermediate()` to reconstruct bosons from decays, e.g. `positron + nu_e => Wp` * Added optional Colour charges to particles (`Particle.colour`) - Colour connection in the HEJ limit can be generated via `Event.generate_colours` (automatically done in the resummation) - Colour configuration of input events can be checked with `Event.is_leading_colour` * Added function `Event.valid_hej_state` to check if event _could have_ been produced by `HEJ` according to the `soft pt regulator` cut on the jets #### New and updated functions * Renamed `EventType::nonHEJ` to `EventType::non_resummable` and `is_HEJ()` to `is_resummable()` such that run card is consistent with internal workings * Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public * New `EventReweighter` member function `treatment` to query the treatment with respect to resummation for the various event types. * Added auxiliary functions to obtain a `std::string` from `EventDescription` (`to_string` for human readable, and `to_simple_string` for easy parsable string) * New `get_analyses` function to read in multiple `HEJ::Analysis` at once, similar to `get_analysis` * New `get_ew_parameters` function to extract electroweak parameters from YAML configuration. #### New classes * New class `Unweighter` to do unweighting * New class `CrossSectionAccumulator` to keep track of Cross Section of the different subprocess * New template struct `Parameters` similar to old `Weights` - `Weights` are now an alias for `Parameters`. Calling `Weights` did not change - `Weights.hh` was replaced by `Parameters.hh`. The old `Weights.hh` header will be removed in version 2.2.0 * Function to multiplication and division of `EventParameters.weight` by double - This can be combined with `Parameters`, e.g. `Parameters*Weights`, see also `Events.parameters()` - Moved `EventParameters` to `Parameters.hh` header * new class `EWConstants` replaces previously hard coded `vev` - `EWConstants` have to be set in the general `Config` and the `MatrixElementConfig` #### Input/Output * New abstract `EventReader` class, as base for reading events from files - Moved LHE file reader to `HEJ::LesHouchesReader` * New (optional) function `finish()` in abstract class `EventWriter`. `finish()` is called _after_ all events are written. * Support reading (`HDF5Reader`) and writing (`HDF5Writer`) `hdf5` files * New `BufferedEventReader` class that allows to put events back into the reader. * New `SherpaLHEReader` to read Sherpa LHE files with correct weights * `get_analysis` now requires `YAML::Node` and `LHEF::HEPRUP` as arguments * Replaced `HepMCInterface` and `HepMCWriter` by `HepMCInterfaceX` and `HepMCWriterX` respectively, with `X` being the major version of HepMC (2 or 3) - Renamed `HepMCInterfaceX::init_kinematics` to `HepMCInterfaceX::init_event` and protected it, use `HepMCInterfaceX::operator()` instead - Removed redundant `HepMCInterfaceX::add_variation` function #### Linking * Export cmake target `HEJ::HEJ` to link directly against `libHEJ` * Preprocessor flags (`HEJ_BUILD_WITH_XYZ`) for enabled/disabled dependencies are now written to `ConfigFlags.hh` * Provide links to version specific object files, e.g. `libHEJ.so -> libHEJ.so.2.1 (soname) -> libHEJ.so.2.1.0` * Removed `LHAPDF` from public interface #### Miscellaneous * Capitalisation of `Config.hh` now matches class `Config` (was `config.hh`) * Renamed `Config::max_ext_soft_pt_fraction` to `Config::soft_pt_regulator`. The new `Config::soft_pt_regulator` parameter is optional and the old `Config::max_ext_soft_pt_fraction` is **deprecated**. * Replaced redundant member `EventReweighterConfig::jet_param` with getter function `EventReweighter.jet_param()` (`JetParameters` are already in `PhaseSpacePointConfig`) * Avoid storing reference to the Random Number Generator inside classes - Constructors of `EventReweighter` now expect `std::shared_ptr` (was reference) - Moved reference to `HEJ::RNG` from constructor of `JetSplitter` to `JetSplitter.split` ## Version 2.0 ### 2.0.4 * Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake` ### 2.0.0 * First release diff --git a/include/HEJ/CombinedEventWriter.hh b/include/HEJ/CombinedEventWriter.hh index 66ff8e6..3933b65 100644 --- a/include/HEJ/CombinedEventWriter.hh +++ b/include/HEJ/CombinedEventWriter.hh @@ -1,50 +1,53 @@ /** \file * \brief Declares the CombinedEventWriter class * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2023 * \copyright GPLv2 or later */ #pragma once #include #include #include "HEJ/EventWriter.hh" namespace LHEF { class HEPRUP; } namespace HEJ { class Event; struct OutputFile; //! Write event output to zero or more output files. class CombinedEventWriter: public EventWriter { public: //! Constructor /** * @param outfiles Specifies files output should be written to. * Each entry in the vector contains a file name * and output format. * @param heprup General process information */ CombinedEventWriter( std::vector const & outfiles, LHEF::HEPRUP const & heprup ); ~CombinedEventWriter() override; //! Write one event to all output files void write(Event const & /*ev*/) override; + //! Set the ratio (cross section) / (sum of event weights) + void set_xs_scale(double scale) override; + //! Finish writing, e.g. write out general informations at end of run //! @warning using the writer after finishing is not guaranteed to work void finish() override; private: std::vector> writers_; }; } // namespace HEJ diff --git a/include/HEJ/EventWriter.hh b/include/HEJ/EventWriter.hh index c4b437e..9e636d1 100644 --- a/include/HEJ/EventWriter.hh +++ b/include/HEJ/EventWriter.hh @@ -1,41 +1,45 @@ /** \file * \brief Header file for the EventWriter interface. * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2023 * \copyright GPLv2 or later */ #pragma once #include namespace HEJ { class Event; //! Pure abstract base class for event writers class EventWriter { public: //! Write an event virtual void write(Event const &) = 0; + + //! Set the ratio (cross section) / (sum of event weights) + virtual void set_xs_scale(double scale) = 0; + //! Finish writing, e.g. write out general informations at end of run //! @warning using the writer after finishing is not guaranteed to work virtual void finish(){finished_=true;}; virtual ~EventWriter() = default; protected: //! @brief If writer is not finished run finish() and abort on error //! //! @param writer writer that is expected to finish, i.e. `this` //! @param name name for error message //! //! Used in the destructor of inherited EventWriter to force finish() //! @warning When finish() failed this will abort the whole program! void finish_or_abort( EventWriter* writer, std::string const & name ) const noexcept; bool finished() const {return finished_;}; private: bool finished_ = false; }; } // namespace HEJ diff --git a/include/HEJ/HDF5Writer.hh b/include/HEJ/HDF5Writer.hh index 2843eca..afddd7e 100644 --- a/include/HEJ/HDF5Writer.hh +++ b/include/HEJ/HDF5Writer.hh @@ -1,57 +1,60 @@ /** \file * \brief Contains the EventWriter for HDF5 Output. * * The output format is specified in arXiv:1905.05120. * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2023 * \copyright GPLv2 or later */ #pragma once #include #include #include "HEJ/EventWriter.hh" namespace LHEF { class HEPRUP; } namespace HEJ { class Event; //! This is an event writer specifically for HDF5 output. /** * \internal Implementation note: This uses the pimpl ("pointer to * implementation") idiom. HDF5 support is optional. Without pimpl, * we would have to specify whether HDF5 is available via the * preprocessor whenever this header is included. We don't want to * burden users of the HEJ library (for example the HEJ fixed-order * generator) with those details */ class HDF5Writer: public EventWriter{ public: //! Constructor /** * @param file name of the output file * @param heprup general process information */ HDF5Writer(std::string const & file, LHEF::HEPRUP heprup); HDF5Writer() = delete; //! Write an event to the output file void write(Event const & ev) override; + //! Set the ratio (cross section) / (sum of event weights) + void set_xs_scale(double scale) override; + //! Finish writing void finish() override; ~HDF5Writer() override; private: struct HDF5WriterImpl; std::unique_ptr impl_; }; } // namespace HEJ diff --git a/include/HEJ/HepMC2Interface.hh b/include/HEJ/HepMC2Interface.hh index 40c918f..f896e58 100644 --- a/include/HEJ/HepMC2Interface.hh +++ b/include/HEJ/HepMC2Interface.hh @@ -1,77 +1,81 @@ /** \file * \brief Header file for the HepMC2Interface * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2023 * \copyright GPLv2 or later */ #pragma once #include #include #include "HEJ/PDG_codes.hh" namespace HepMC { class GenCrossSection; class GenEvent; } namespace LHEF { class HEPRUP; } namespace HEJ { class Event; //! This class converts the Events into HepMC::GenEvents /** * \details The output is depended on the HepMC version HEJ is compiled with, * both HepMC 2 and HepMC 3 are supported. If HEJ 2 is compiled * without HepMC calling this interface will throw an error. * * This interface will also keep track of the cross section of all the events that * being fed into it. */ class HepMC2Interface{ public: HepMC2Interface(LHEF::HEPRUP const & /*heprup*/); /** * \brief main function to convert an event into HepMC::GenEvent * * \param event Event to convert * \param weight_index optional selection of specific weight * (negative value gives central weight) */ HepMC::GenEvent operator()(Event const & event, int weight_index = -1); /** * \brief Sets the central value from \p event to \p out_ev * * \param out_ev HepMC::GenEvent to write to * \param event Event to convert * \param weight_index optional selection of specific weight * (negative value gives "central") * * This overwrites the central value of \p out_ev. */ void set_central(HepMC::GenEvent & out_ev, Event const & event, int weight_index = -1); + //! Set the ratio (cross section) / (sum of event weights) + void set_xs_scale(double scale); + protected: /** * \brief initialise generic event infomrations (not central weights) * * \param event Event to convert */ HepMC::GenEvent init_event(Event const & event) const; private: std::array beam_particle_{}; std::array beam_energy_{}; std::size_t event_count_; double tot_weight_; double tot_weight2_; + double xs_scale_ = 1.; HepMC::GenCrossSection cross_section() const; }; } // namespace HEJ diff --git a/include/HEJ/HepMC2Writer.hh b/include/HEJ/HepMC2Writer.hh index 004c08a..ba4b354 100644 --- a/include/HEJ/HepMC2Writer.hh +++ b/include/HEJ/HepMC2Writer.hh @@ -1,53 +1,56 @@ /** \file * \brief Contains the EventWriter for HepMC Output. * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2023 * \copyright GPLv2 or later */ #pragma once #include #include #include "HEJ/EventWriter.hh" namespace LHEF { class HEPRUP; } namespace HEJ { class Event; //! This is an event writer specifically for HepMC output. /** * \internal Implementation note: * This uses the pimpl ("pointer to implementation") idiom. * HepMC support is optional and the implementation depends on the * HepMC version. Without pimpl, we would have to specify the HepMC version * via the preprocessor whenever this header is included. We don't want to * burden users of the HEJ library (for example the HEJ fixed-order generator) * with those details */ class HepMC2Writer: public EventWriter{ public: //! Constructor /** * @param file name of the output file * @param heprup general process information */ HepMC2Writer(std::string const & file, LHEF::HEPRUP heprup); HepMC2Writer() = delete; //! Write an event to the output file void write(Event const & ev) override; + //! Set the ratio (cross section) / (sum of event weights) + void set_xs_scale(double scale) override; + ~HepMC2Writer() override; private: struct HepMC2WriterImpl; std::unique_ptr impl_; }; } // namespace HEJ diff --git a/include/HEJ/HepMC3Interface.hh b/include/HEJ/HepMC3Interface.hh index 94f821f..83e9d6e 100644 --- a/include/HEJ/HepMC3Interface.hh +++ b/include/HEJ/HepMC3Interface.hh @@ -1,89 +1,93 @@ /** \file * \brief Header file for the HepMC3Interface * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2023 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "HEJ/PDG_codes.hh" namespace HepMC3 { class GenCrossSection; class GenEvent; class GenRunInfo; } namespace LHEF { class HEPRUP; } namespace HEJ { class Event; //! This class converts the Events into HepMC3::GenEvents /** * \details The output is depended on the HepMC3 version HEJ is compiled with, * both HepMC3 2 and HepMC3 3 are supported. If HEJ 2 is compiled * without HepMC3 calling this interface will throw an error. * * This interface will also keep track of the cross section of all the events that * being fed into it. */ class HepMC3Interface{ public: HepMC3Interface(LHEF::HEPRUP heprup); ~HepMC3Interface(); HepMC3Interface & operator=(HepMC3Interface const & other) = default; HepMC3Interface(HepMC3Interface const & other) = default; HepMC3Interface & operator=(HepMC3Interface && other) = default; HepMC3Interface(HepMC3Interface && other) = default; /** * \brief main function to convert an event into HepMC3::GenEvent * * \param event Event to convert * \param weight_index optional selection of specific weight * (negative value gives central weight) */ HepMC3::GenEvent operator()(Event const & event, int weight_index = -1); /** * \brief Sets the central value from \p event to \p out_ev * * \param out_ev HepMC3::GenEvent to write to * \param event Event to convert * \param weight_index optional selection of specific weight * (negative value gives "central") * * This overwrites the central value of \p out_ev. */ void set_central(HepMC3::GenEvent & out_ev, Event const & event, int weight_index = -1); //! Pointer to generic run informations std::shared_ptr run_info(){ return run_info_; } + //! Set the ratio (cross section) / (sum of event weights) + void set_xs_scale(double scale); + protected: /** * \brief initialise generic event infomrations (not central weights) * * \param event Event to convert */ HepMC3::GenEvent init_event(Event const & event) const; private: std::array beam_particle_{}; std::array beam_energy_{}; std::shared_ptr run_info_; std::size_t event_count_; double tot_weight_; double tot_weight2_; + double xs_scale_ = 1.; std::shared_ptr xs_; }; } // namespace HEJ diff --git a/include/HEJ/HepMC3Writer.hh b/include/HEJ/HepMC3Writer.hh index 05a6b88..b628520 100644 --- a/include/HEJ/HepMC3Writer.hh +++ b/include/HEJ/HepMC3Writer.hh @@ -1,55 +1,58 @@ /** \file * \brief Contains the EventWriter for HepMC3 Output. * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2023 * \copyright GPLv2 or later */ #pragma once #include #include #include "HEJ/EventWriter.hh" namespace LHEF { class HEPRUP; } namespace HEJ { class Event; //! This is an event writer specifically for HepMC3 output. /** * \internal Implementation note: * This uses the pimpl ("pointer to implementation") idiom. * HepMC3 support is optional and the implementation depends on the * HepMC3 version. Without pimpl, we would have to specify the HepMC3 version * via the preprocessor whenever this header is included. We don't want to * burden users of the HEJ library (for example the HEJ fixed-order generator) * with those details */ class HepMC3Writer: public EventWriter{ public: //! Constructor /** * @param file name of the output file * @param heprup general process information */ HepMC3Writer(std::string const & file, LHEF::HEPRUP heprup); HepMC3Writer() = delete; //! Write an event to the output file void write(Event const & ev) override; + //! Set the ratio (cross section) / (sum of event weights) + void set_xs_scale(double scale) override; + //! Finish writing void finish() override; ~HepMC3Writer() override; private: struct HepMC3WriterImpl; std::unique_ptr impl_; }; } // namespace HEJ diff --git a/include/HEJ/LesHouchesWriter.hh b/include/HEJ/LesHouchesWriter.hh index 077ec71..60c7652 100644 --- a/include/HEJ/LesHouchesWriter.hh +++ b/include/HEJ/LesHouchesWriter.hh @@ -1,61 +1,65 @@ /** \file * \brief Contains the writer for LesHouches output * * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2019-2020 + * \date 2019-2023 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "LHEF/LHEF.h" #include "HEJ/EventWriter.hh" namespace HEJ { class Event; //! Class for writing events to a file in the Les Houches Event File format class LesHouchesWriter : public EventWriter{ public: //! Constructor /** * @param file Name of output file * @param heprup General process information */ LesHouchesWriter(std::string const & file, LHEF::HEPRUP heprup); LesHouchesWriter(LesHouchesWriter const & other) = delete; LesHouchesWriter & operator=(LesHouchesWriter const & other) = delete; /** @TODO in principle, this class should be movable * but that somehow(?) breaks the write member function */ LesHouchesWriter(LesHouchesWriter && other) = delete; LesHouchesWriter & operator=(LesHouchesWriter && other) = delete; ~LesHouchesWriter() override; //! Write an event to the file specified in the constructor void write(Event const & ev) override; + //! Set the ratio (cross section) / (sum of event weights) + void set_xs_scale(double scale) override; + //! Dump general information & finish writing void finish() override; private: void write_init(){ writer_->init(); } LHEF::HEPRUP & heprup(){ return writer_->heprup; } LHEF::HEPEUP & hepeup(){ return writer_->hepeup; } std::fstream out_; std::unique_ptr writer_; + double xs_scale_ = 1.; }; } // namespace HEJ diff --git a/src/CombinedEventWriter.cc b/src/CombinedEventWriter.cc index 3139d8b..4951b45 100644 --- a/src/CombinedEventWriter.cc +++ b/src/CombinedEventWriter.cc @@ -1,41 +1,47 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/CombinedEventWriter.hh" #include #include "HEJ/exceptions.hh" #include "HEJ/make_writer.hh" #include "HEJ/output_formats.hh" namespace HEJ { CombinedEventWriter::CombinedEventWriter( std::vector const & outfiles, LHEF::HEPRUP const & heprup ){ writers_.reserve(outfiles.size()); for(OutputFile const & outfile: outfiles){ writers_.emplace_back( make_format_writer(outfile.format, outfile.name, heprup) ); } } void CombinedEventWriter::write(Event const & ev){ for(auto & writer: writers_) writer->write(ev); } + void CombinedEventWriter::set_xs_scale(const double scale) { + for(auto & writer: writers_) { + writer->set_xs_scale(scale); + } + } + void CombinedEventWriter::finish(){ EventWriter::finish(); for(auto & writer: writers_) writer->finish(); } CombinedEventWriter::~CombinedEventWriter(){ finish_or_abort(this, "CombinedEventWriter"); } } // namespace HEJ diff --git a/src/HDF5Writer.cc b/src/HDF5Writer.cc index ea90b58..06508c1 100644 --- a/src/HDF5Writer.cc +++ b/src/HDF5Writer.cc @@ -1,435 +1,451 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/HDF5Writer.hh" #include #include "LHEF/LHEF.h" #include "HEJ/ConfigFlags.hh" #ifdef HEJ_BUILD_WITH_HDF5 #include #include #include #include #include #include #include "highfive/H5File.hpp" #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #else #include "HEJ/exceptions.hh" #endif #ifdef HEJ_BUILD_WITH_HDF5 namespace HEJ { using HighFive::File; using HighFive::DataSpace; namespace { using std::size_t; constexpr size_t CHUNK_SIZE = 1000; constexpr unsigned COMPRESSION_LEVEL = 3; size_t to_index(event_type::EventType const type){ return type==0?0:std::floor(std::log2(static_cast(type)))+1; } template void write_dataset(HighFive::Group & group, std::string const & name, T val) { using data_t = std::decay_t; group.createDataSet(name, DataSpace::From(val)).write(val); } template void write_dataset( HighFive::Group & group, std::string const & name, std::vector const & val ) { using data_t = std::decay_t; group.createDataSet(name, DataSpace::From(val)).write(val); } struct Cache { explicit Cache(size_t capacity): capacity{capacity} { nparticles.reserve(capacity); start.reserve(capacity); pid.reserve(capacity); weight.reserve(capacity); scale.reserve(capacity); fscale.reserve(capacity); rscale.reserve(capacity); aqed.reserve(capacity); aqcd.reserve(capacity); trials.reserve(capacity); npLO.reserve(capacity); npNLO.reserve(capacity); } void fill(Event const & ev) { const auto hepeup = to_HEPEUP(ev, nullptr); // HEJ event to get nice wrapper const auto num_partons = std::distance(ev.cbegin_partons(), ev.cend_partons()); assert(num_partons>0); // Number of partons for LO matching, HEJ requires at least 2 partons npLO.emplace_back(num_partons>1?num_partons-2:num_partons); // Number of real emissions in NLO, HEJ is LO -> -1 npNLO.emplace_back(-1); fill_event_params(hepeup); fill_event_particles(hepeup); } void fill_event_params(LHEF::HEPEUP const & ev) { nparticles.emplace_back(ev.NUP); start.emplace_back(particle_pos); pid.emplace_back(ev.IDPRUP); weight.emplace_back(ev.XWGTUP); scale.emplace_back(ev.SCALUP); fscale.emplace_back(ev.scales.muf); rscale.emplace_back(ev.scales.mur); aqed.emplace_back(ev.AQEDUP); aqcd.emplace_back(ev.AQCDUP); // set first trial=1 for first event // -> sum(trials) = 1 -> xs=sum(weights)/sum(trials) as in Sherpa if(particle_pos == 0){ trials.emplace_back(1.); } else { trials.emplace_back(0.); } particle_pos += ev.NUP; } void fill_event_particles(LHEF::HEPEUP const & ev) { id.insert(end(id), begin(ev.IDUP), end(ev.IDUP)); status.insert(end(status), begin(ev.ISTUP), end(ev.ISTUP)); lifetime.insert(end(lifetime), begin(ev.VTIMUP), end(ev.VTIMUP)); spin.insert(end(spin), begin(ev.SPINUP), end(ev.SPINUP)); for(int i = 0; i < ev.NUP; ++i) { mother1.emplace_back(ev.MOTHUP[i].first); mother2.emplace_back(ev.MOTHUP[i].second); color1.emplace_back(ev.ICOLUP[i].first); color2.emplace_back(ev.ICOLUP[i].second); px.emplace_back(ev.PUP[i][0]); py.emplace_back(ev.PUP[i][1]); pz.emplace_back(ev.PUP[i][2]); e.emplace_back(ev.PUP[i][3]); m.emplace_back(ev.PUP[i][4]); } } bool is_full() const { return nparticles.size() >= capacity; } void clear() { nparticles.clear(); start.clear(); pid.clear(); id.clear(); status.clear(); mother1.clear(); mother2.clear(); color1.clear(); color2.clear(); weight.clear(); scale.clear(); fscale.clear(); rscale.clear(); aqed.clear(); aqcd.clear(); trials.clear(); npLO.clear(); npNLO.clear(); px.clear(); py.clear(); pz.clear(); e.clear(); m.clear(); lifetime.clear(); spin.clear(); } size_t capacity; std::vector nparticles, start, pid, id, status, mother1, mother2, color1, color2, npLO, npNLO; std::vector weight, scale, fscale, rscale, aqed, aqcd, trials, px, py, pz, e, m, lifetime, spin; size_t particle_pos = 0; }; } // namespace struct HDF5Writer::HDF5WriterImpl: EventWriter{ File file; LHEF::HEPRUP heprup; Cache cache{CHUNK_SIZE}; size_t event_idx = 0; size_t particle_idx = 0; + double xs_scale = 1.; HDF5WriterImpl(std::string const & filename, LHEF::HEPRUP && hepr): file{filename, File::ReadWrite | File::Create | File::Truncate}, heprup{std::forward(hepr)} { // TODO: code duplication with Les Houches Writer const size_t max_number_types = to_index(event_type::last_type)+1; heprup.NPRUP = max_number_types; // ids of event types heprup.LPRUP.clear(); heprup.LPRUP.reserve(max_number_types); heprup.LPRUP.emplace_back(0); for(size_t i=event_type::first_type+1; i<=event_type::last_type; i*=2) { heprup.LPRUP.emplace_back(i); } heprup.XSECUP = std::vector(max_number_types, 0.); heprup.XERRUP = std::vector(max_number_types, 0.); heprup.XMAXUP = std::vector(max_number_types, 0.); write_init(); create_event_group(); create_particle_group(); } void write_init() { auto init = file.createGroup("init"); write_dataset(init, "PDFgroupA" , heprup.PDFGUP.first); write_dataset(init, "PDFgroupB" , heprup.PDFGUP.second); write_dataset(init, "PDFsetA" , heprup.PDFSUP.first); write_dataset(init, "PDFsetB" , heprup.PDFSUP.second); write_dataset(init, "beamA" , heprup.IDBMUP.first); write_dataset(init, "beamB" , heprup.IDBMUP.second); write_dataset(init, "energyA" , heprup.EBMUP.first); write_dataset(init, "energyB" , heprup.EBMUP.second); write_dataset(init, "numProcesses" , heprup.NPRUP); write_dataset(init, "weightingStrategy", heprup.IDWTUP); auto proc_info = file.createGroup("procInfo"); write_dataset(proc_info, "procId", heprup.LPRUP); } static HighFive::DataSetCreateProps const & hdf5_chunk() { static const auto props = [](){ HighFive::DataSetCreateProps props; props.add(HighFive::Chunking({CHUNK_SIZE})); props.add(HighFive::Deflate(COMPRESSION_LEVEL)); return props; }(); return props; } void create_event_group() { static const auto dim = DataSpace({0}, {DataSpace::UNLIMITED}); auto events = file.createGroup("event"); events.createDataSet("nparticles", dim, hdf5_chunk()); events.createDataSet("start", dim, hdf5_chunk()); events.createDataSet("pid", dim, hdf5_chunk()); events.createDataSet("weight", dim, hdf5_chunk()); events.createDataSet("scale", dim, hdf5_chunk()); events.createDataSet("fscale", dim, hdf5_chunk()); events.createDataSet("rscale", dim, hdf5_chunk()); events.createDataSet("aqed", dim, hdf5_chunk()); events.createDataSet("aqcd", dim, hdf5_chunk()); events.createDataSet("trials", dim, hdf5_chunk()); events.createDataSet("npLO", dim, hdf5_chunk()); events.createDataSet("npNLO", dim, hdf5_chunk()); } void resize_event_group(size_t new_size) { auto events = file.getGroup("event"); events.getDataSet("nparticles").resize({new_size}); events.getDataSet("start").resize({new_size}); events.getDataSet("pid").resize({new_size}); events.getDataSet("weight").resize({new_size}); events.getDataSet("scale").resize({new_size}); events.getDataSet("fscale").resize({new_size}); events.getDataSet("rscale").resize({new_size}); events.getDataSet("aqed").resize({new_size}); events.getDataSet("aqcd").resize({new_size}); events.getDataSet("trials").resize({new_size}); events.getDataSet("npLO").resize({new_size}); events.getDataSet("npNLO").resize({new_size}); } void create_particle_group() { static const auto dim = DataSpace({0}, {DataSpace::UNLIMITED}); auto particles = file.createGroup("particle"); particles.createDataSet("id", dim, hdf5_chunk()); particles.createDataSet("status", dim, hdf5_chunk()); particles.createDataSet("mother1", dim, hdf5_chunk()); particles.createDataSet("mother2", dim, hdf5_chunk()); particles.createDataSet("color1", dim, hdf5_chunk()); particles.createDataSet("color2", dim, hdf5_chunk()); particles.createDataSet("px", dim, hdf5_chunk()); particles.createDataSet("py", dim, hdf5_chunk()); particles.createDataSet("pz", dim, hdf5_chunk()); particles.createDataSet("e", dim, hdf5_chunk()); particles.createDataSet("m", dim, hdf5_chunk()); particles.createDataSet("lifetime", dim, hdf5_chunk()); particles.createDataSet("spin", dim, hdf5_chunk()); } void resize_particle_group(size_t new_size) { auto particles = file.getGroup("particle"); particles.getDataSet("id").resize({new_size}); particles.getDataSet("status").resize({new_size}); particles.getDataSet("mother1").resize({new_size}); particles.getDataSet("mother2").resize({new_size}); particles.getDataSet("color1").resize({new_size}); particles.getDataSet("color2").resize({new_size}); particles.getDataSet("px").resize({new_size}); particles.getDataSet("py").resize({new_size}); particles.getDataSet("pz").resize({new_size}); particles.getDataSet("e").resize({new_size}); particles.getDataSet("m").resize({new_size}); particles.getDataSet("lifetime").resize({new_size}); particles.getDataSet("spin").resize({new_size}); } void write(Event const & ev) override { cache.fill(ev); if(cache.is_full()) { dump_cache(); } const double wt = ev.central().weight; const size_t idx = to_index(ev.type()); heprup.XSECUP[idx] += wt; heprup.XERRUP[idx] += wt*wt; if(wt > heprup.XMAXUP[idx]){ heprup.XMAXUP[idx] = wt; } } + void set_xs_scale(const double scale) override { + xs_scale = scale; + } + void dump_cache() { write_event_params(); write_event_particles(); cache.clear(); } void write_event_params() { auto events = file.getGroup("event"); // choose arbitrary dataset to find size const auto dataset = events.getDataSet("nparticles"); const size_t size = dataset.getSpace().getDimensions().front(); resize_event_group(size + cache.nparticles.size()); // NOLINTNEXTLINE #define WRITE_FROM_CACHE(GROUP, PROPERTY) \ GROUP.getDataSet(#PROPERTY).select({size}, {cache.PROPERTY.size()}).write(cache.PROPERTY) WRITE_FROM_CACHE(events, nparticles); WRITE_FROM_CACHE(events, start); WRITE_FROM_CACHE(events, pid); WRITE_FROM_CACHE(events, weight); WRITE_FROM_CACHE(events, scale); WRITE_FROM_CACHE(events, fscale); WRITE_FROM_CACHE(events, rscale); WRITE_FROM_CACHE(events, aqed); WRITE_FROM_CACHE(events, aqcd); WRITE_FROM_CACHE(events, trials); WRITE_FROM_CACHE(events, npLO); WRITE_FROM_CACHE(events, npNLO); } void write_event_particles() { auto particles = file.getGroup("particle"); // choose arbitrary dataset to find size const auto dataset = particles.getDataSet("id"); const size_t size = dataset.getSpace().getDimensions().front(); resize_particle_group(size + cache.id.size()); WRITE_FROM_CACHE(particles, id); WRITE_FROM_CACHE(particles, status); WRITE_FROM_CACHE(particles, lifetime); WRITE_FROM_CACHE(particles, spin); WRITE_FROM_CACHE(particles, mother1); WRITE_FROM_CACHE(particles, mother2); WRITE_FROM_CACHE(particles, color1); WRITE_FROM_CACHE(particles, color2); WRITE_FROM_CACHE(particles, px); WRITE_FROM_CACHE(particles, py); WRITE_FROM_CACHE(particles, pz); WRITE_FROM_CACHE(particles, e); WRITE_FROM_CACHE(particles, m); } #undef WRITE_FROM_CACHE void finish() override { if(finished()) throw std::ios_base::failure("HDF5Writer writer already finished."); EventWriter::finish(); dump_cache(); auto proc_info = file.getGroup("procInfo"); + for(auto & xs: heprup.XSECUP) { + xs *= xs_scale; + } write_dataset(proc_info, "xSection", heprup.XSECUP); for(auto & xs_err: heprup.XERRUP) { - xs_err = std::sqrt(xs_err); + xs_err = xs_scale*std::sqrt(xs_err); } write_dataset(proc_info, "error", heprup.XERRUP); write_dataset(proc_info, "unitWeight", heprup.XMAXUP); file.flush(); } ~HDF5WriterImpl() override { finish_or_abort(this, "HDF5Writer"); } }; HDF5Writer::HDF5Writer(std::string const & filename, LHEF::HEPRUP heprup): impl_{std::make_unique( filename, std::move(heprup))} {} void HDF5Writer::write(Event const & ev){ impl_->write(ev); } + void HDF5Writer::set_xs_scale(const double scale) { + impl_->set_xs_scale(scale); + } + void HDF5Writer::finish(){ impl_->finish(); } } // namespace HEJ #else // no HDF5 support namespace HEJ { struct HDF5Writer::HDF5WriterImpl{}; HDF5Writer::HDF5Writer(std::string const & /*file*/, LHEF::HEPRUP /*heprup*/){ throw std::invalid_argument{ "Failed to create HDF5 writer: " "HEJ 2 was built without HDF5 support" }; } void HDF5Writer::write(Event const & /*ev*/){ assert(false); } + void HDF5Writer::set_xs_scale(double /* scale */) { + assert(false); + } + void HDF5Writer::finish(){ assert(false); } } #endif namespace HEJ { HDF5Writer::~HDF5Writer() = default; } diff --git a/src/HepMC2Interface.cc b/src/HepMC2Interface.cc index edba099..38307ec 100644 --- a/src/HepMC2Interface.cc +++ b/src/HepMC2Interface.cc @@ -1,165 +1,172 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/HepMC2Interface.hh" #include "HEJ/ConfigFlags.hh" #include "HEJ/exceptions.hh" #ifdef HEJ_BUILD_WITH_HepMC2 #include #include #include #include "LHEF/LHEF.h" #include "HepMC/GenCrossSection.h" #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/GenVertex.h" #include "HepMC/SimpleVector.h" #include "HepMC/Units.h" #include "HEJ/Event.hh" #include "HEJ/Particle.hh" #include "HEJ/detail/HepMCInterface_common.hh" #else // no HepMC2 #include "HEJ/utility.hh" #endif #ifdef HEJ_BUILD_WITH_HepMC2 namespace HEJ { namespace detail_HepMC { template<> struct HepMCVersion<2> { using GenEvent = HepMC::GenEvent; using Beam = std::array; }; template<> auto make_particle_ptr<2> ( Particle const & sp, int status ) { return new HepMC::GenParticle( to_FourVector(sp), static_cast (sp.type), status ); } template<> auto make_vx_ptr<2>() { return new HepMC::GenVertex(); } } // namespace detail_HepMC HepMC2Interface::HepMC2Interface(LHEF::HEPRUP const & heprup): beam_particle_{static_cast(heprup.IDBMUP.first), static_cast(heprup.IDBMUP.second)}, beam_energy_{heprup.EBMUP.first, heprup.EBMUP.second}, event_count_(0.), tot_weight_(0.), tot_weight2_(0.) {} HepMC::GenCrossSection HepMC2Interface::cross_section() const { HepMC::GenCrossSection xs; - xs.set_cross_section(tot_weight_, std::sqrt(tot_weight2_)); + xs.set_cross_section( + xs_scale_ * tot_weight_, + xs_scale_ * std::sqrt(tot_weight2_) + ); return xs; } + void HepMC2Interface::set_xs_scale(const double scale) { + xs_scale_ = scale; + } + HepMC::GenEvent HepMC2Interface::init_event(Event const & event) const { const std::array beam { new HepMC::GenParticle( HepMC::FourVector(0,0,-beam_energy_[0],beam_energy_[0]), beam_particle_[0], detail_HepMC::Status::beam ), new HepMC::GenParticle( HepMC::FourVector(0,0, beam_energy_[1],beam_energy_[1]), beam_particle_[1], detail_HepMC::Status::beam ) }; auto hepmc_ev{ detail_HepMC::HepMC_init_kinematics<2>( event, beam, HepMC::GenEvent{ HepMC::Units::GEV, HepMC::Units::MM } ) }; hepmc_ev.weights().push_back( event.central().weight ); for(auto const & var: event.variations()){ hepmc_ev.weights().push_back( var.weight ); // no weight name for HepMC2 since rivet3 seem to mix them up // (could be added via hepmc_ev.weights()[name]=weight) } return hepmc_ev; } void HepMC2Interface::set_central( HepMC::GenEvent & out_ev, Event const & event, int const weight_index ){ EventParameters event_param; if(weight_index < 0) event_param = event.central(); else if ( static_cast(weight_index) < event.variations().size()) event_param = event.variations(weight_index); else throw std::invalid_argument{ "HepMC2Interface tried to access a weight outside of the variation range." }; const double wt = event_param.weight; tot_weight_ += wt; tot_weight2_ += wt * wt; ++event_count_; // central always on first assert(out_ev.weights().size() == event.variations().size()+1); out_ev.weights()[0] = wt; out_ev.set_cross_section( cross_section() ); out_ev.set_signal_process_id(event.type()); out_ev.set_event_scale(event_param.mur); out_ev.set_event_number(event_count_); /// @TODO add alphaQCD (need function) and alphaQED /// @TODO output pdf (currently not avaiable from event alone) } HepMC::GenEvent HepMC2Interface::operator()(Event const & event, int const weight_index ){ HepMC::GenEvent out_ev(init_event(event)); set_central(out_ev, event, weight_index); return out_ev; } } // namespace HEJ #else // no HepMC2 => empty class namespace HepMC { class GenEvent {}; class GenCrossSection {}; } namespace HEJ { HepMC2Interface::HepMC2Interface(LHEF::HEPRUP const & /*heprup*/){ ignore(beam_particle_,beam_energy_,event_count_,tot_weight_,tot_weight2_); throw std::invalid_argument( "Failed to create HepMC2Interface: " "HEJ 2 was built without HepMC2 support" ); } HepMC::GenEvent HepMC2Interface::operator()( Event const & /*event*/, int /*weight_index*/ ){return HepMC::GenEvent();} HepMC::GenEvent HepMC2Interface::init_event(Event const & /*event*/) const {return HepMC::GenEvent();} void HepMC2Interface::set_central( HepMC::GenEvent & /*out_ev*/, Event const & /*event*/, int /*weight_index*/ ){} HepMC::GenCrossSection HepMC2Interface::cross_section() const {return HepMC::GenCrossSection();} } #endif diff --git a/src/HepMC2Writer.cc b/src/HepMC2Writer.cc index 04f77e0..d001305 100644 --- a/src/HepMC2Writer.cc +++ b/src/HepMC2Writer.cc @@ -1,88 +1,101 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/HepMC2Writer.hh" #include "LHEF/LHEF.h" #include "HEJ/ConfigFlags.hh" #ifdef HEJ_BUILD_WITH_HepMC2 #include #include "HepMC/IO_GenEvent.h" #include "HEJ/HepMC2Interface.hh" #else #include #include "HEJ/exceptions.hh" #endif #ifdef HEJ_BUILD_WITH_HepMC2 namespace HEJ { struct HepMC2Writer::HepMC2WriterImpl{ HepMC2Interface hepmc_; HepMC2WriterImpl & operator=(HepMC2WriterImpl const & other) = delete; HepMC2WriterImpl(HepMC2WriterImpl const & other) = delete; HepMC2WriterImpl & operator=(HepMC2WriterImpl && other) = delete; HepMC2WriterImpl(HepMC2WriterImpl && other) = delete; ~HepMC2WriterImpl() = default; HepMC::IO_GenEvent writer_; HepMC2WriterImpl( std::string const & file, LHEF::HEPRUP && heprup ): hepmc_(heprup), writer_{file} {} void write(Event const & ev){ auto out_ev = hepmc_(ev); writer_.write_event(&out_ev); } + + void set_xs_scale(const double scale) { + hepmc_.set_xs_scale(scale); + } + }; HepMC2Writer::HepMC2Writer(std::string const & file, LHEF::HEPRUP heprup): impl_{std::make_unique(file, std::move(heprup))} {} + void HepMC2Writer::set_xs_scale(const double scale) { + impl_->set_xs_scale(scale); + } + void HepMC2Writer::write(Event const & ev){ impl_->write(ev); } } // namespace HEJ #else // no HepMC2 namespace HEJ { struct HepMC2Writer::HepMC2WriterImpl{}; HepMC2Writer::HepMC2Writer( std::string const & /*file*/, LHEF::HEPRUP /*heprup*/ ){ throw std::invalid_argument( "Failed to create HepMC writer: " "HEJ 2 was built without HepMC2 support" ); } void HepMC2Writer::write(Event const & /*ev*/){ assert(false); } + void HepMC2Writer::set_xs_scale(double /* scale */) { + assert(false); + } + } #endif namespace HEJ { HepMC2Writer::~HepMC2Writer() = default; } diff --git a/src/HepMC3Interface.cc b/src/HepMC3Interface.cc index d6c2714..48a7dd2 100644 --- a/src/HepMC3Interface.cc +++ b/src/HepMC3Interface.cc @@ -1,226 +1,233 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/HepMC3Interface.hh" #include "HEJ/ConfigFlags.hh" #include "HEJ/exceptions.hh" // include before HepMC3 to override include of "HepMC3/LHEF.h" // (theoretically both files should be the same) #include "LHEF/LHEF.h" #ifdef HEJ_BUILD_WITH_HepMC3 #include #include #include #include #include #include "HepMC3/Attribute.h" #include "HepMC3/FourVector.h" #include "HepMC3/GenCrossSection.h" #include "HepMC3/GenCrossSection_fwd.h" #include "HepMC3/GenEvent.h" #include "HepMC3/GenParticle.h" #include "HepMC3/GenParticle_fwd.h" #include "HepMC3/GenRunInfo.h" #include "HepMC3/GenVertex.h" #include "HepMC3/LHEFAttributes.h" #include "HepMC3/Units.h" #include "HEJ/Event.hh" #include "HEJ/Parameters.hh" #include "HEJ/Particle.hh" #include "HEJ/detail/HepMCInterface_common.hh" #else // no HepMC3 #include "HEJ/utility.hh" #endif #ifdef HEJ_BUILD_WITH_HepMC3 namespace HEJ { namespace detail_HepMC { template<> struct HepMCVersion<3> { using GenEvent = HepMC3::GenEvent; using Beam = std::array; }; template<> auto make_particle_ptr<3> ( Particle const & sp, int status ) { return HepMC3::make_shared( to_FourVector(sp), static_cast (sp.type), status ); } template<> auto make_vx_ptr<3>() { return HepMC3::make_shared(); } } // namespace detail_HepMC namespace { void reset_weight_info(LHEF::HEPRUP & heprup){ heprup.IDWTUP = 2; // use placeholders for unknown init block values // we can overwrite them after processing all events heprup.XSECUP = {0.}; heprup.XERRUP = {0.}; heprup.XMAXUP = {0.}; } HepMC3::shared_ptr init_runinfo(LHEF::HEPRUP heprup){ reset_weight_info(heprup); auto runinfo{ HepMC3::make_shared() }; auto hepr{ HepMC3::make_shared() }; hepr->heprup = std::move(heprup); runinfo->add_attribute(std::string("HEPRUP"), hepr); for(auto const & gen: hepr->heprup.generators){ runinfo->tools().emplace_back( HepMC3::GenRunInfo::ToolInfo{gen.name, gen.version, gen.contents} ); } return runinfo; } std::vector get_weight_names(Event const & ev){ std::vector names; names.reserve(ev.variations().size()+1); // +1 from central names.emplace_back(""); // rivet assumes central band to have no name for(auto const & var : ev.variations()){ if(var.description){ names.emplace_back( to_simple_string(*var.description) ); } else { names.emplace_back( "" ); } } assert(names.size() == ev.variations().size()+1); return names; } } // namespace HepMC3Interface::HepMC3Interface(LHEF::HEPRUP heprup): beam_particle_{static_cast(heprup.IDBMUP.first), static_cast(heprup.IDBMUP.second)}, beam_energy_{heprup.EBMUP.first, heprup.EBMUP.second}, run_info_{ init_runinfo(std::move(heprup)) }, event_count_(0.), tot_weight_(0.), tot_weight2_(0.), xs_{std::make_shared()} { } HepMC3::GenEvent HepMC3Interface::init_event(Event const & event) const { const std::array beam { HepMC3::make_shared( HepMC3::FourVector(0,0,-beam_energy_[0],beam_energy_[0]), beam_particle_[0], detail_HepMC::Status::beam ), HepMC3::make_shared( HepMC3::FourVector(0,0, beam_energy_[1],beam_energy_[1]), beam_particle_[1], detail_HepMC::Status::beam ) }; auto hepmc_ev{ detail_HepMC::HepMC_init_kinematics<3>( event, beam, HepMC3::GenEvent{ HepMC3::Units::GEV, HepMC3::Units::MM } ) }; // set up run specific informations if( run_info_->weight_names().size() != event.variations().size()+1 ){ run_info_->set_weight_names( get_weight_names(event) ); } // order matters: weights in hepmc_ev initialised when registering run_info hepmc_ev.set_run_info(run_info_); assert(hepmc_ev.weights().size() == event.variations().size()+1); for(size_t i=0; i(weight_index) < event.variations().size()) event_param = event.variations(weight_index); else throw std::invalid_argument{ "HepMC3Interface tried to access a weight outside of the variation range." }; const double wt = event_param.weight; tot_weight_ += wt; tot_weight2_ += wt * wt; ++event_count_; // central always on first assert(out_ev.weights().size() == event.variations().size()+1); out_ev.weights()[0] = wt; // out_ev can be setup with a different central scale -> save xs manually out_ev.set_cross_section(xs_); assert(out_ev.cross_section() && out_ev.cross_section() == xs_); // overwrites all previously set xs ... - xs_->set_cross_section(tot_weight_,std::sqrt(tot_weight2_)); + xs_->set_cross_section( + xs_scale_ * tot_weight_, + xs_scale_ * std::sqrt(tot_weight2_) + ); out_ev.set_event_number(event_count_); /// @TODO add number of attempted events xs_->set_accepted_events(event_count_); /// @TODO add alphaQCD (need function) and alphaQED /// @TODO output pdf (currently not avaiable from event alone) } HepMC3::GenEvent HepMC3Interface::operator()( Event const & event, int const weight_index ){ HepMC3::GenEvent out_ev(init_event(event)); set_central(out_ev, event, weight_index); return out_ev; } + void HepMC3Interface::set_xs_scale(const double scale) { + xs_scale_ = scale; + } + } // namespace HEJ #else // no HepMC3 => empty class namespace HepMC3 { class GenEvent {}; class GenCrossSection {}; class GenRunInfo {}; } namespace HEJ { HepMC3Interface::HepMC3Interface(LHEF::HEPRUP /*heprup*/){ ignore(beam_particle_,beam_energy_,event_count_,tot_weight_,tot_weight2_); throw std::invalid_argument( "Failed to create HepMC3Interface: " "HEJ 2 was built without HepMC3 support" ); } HepMC3::GenEvent HepMC3Interface::operator()( Event const & /*event*/, int /*weight_index*/ ){return HepMC3::GenEvent();} HepMC3::GenEvent HepMC3Interface::init_event(Event const & /*event*/) const {return HepMC3::GenEvent();} void HepMC3Interface::set_central( HepMC3::GenEvent & /*out_ev*/, Event const & /*event*/, int /*weight_index*/ ){} } #endif namespace HEJ { HepMC3Interface::~HepMC3Interface() = default; } diff --git a/src/HepMC3Writer.cc b/src/HepMC3Writer.cc index 5059cbf..18b7f9c 100644 --- a/src/HepMC3Writer.cc +++ b/src/HepMC3Writer.cc @@ -1,111 +1,124 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/HepMC3Writer.hh" #include "LHEF/LHEF.h" #include "HEJ/ConfigFlags.hh" #include "HEJ/exceptions.hh" #ifdef HEJ_BUILD_WITH_HepMC3 #include #include "HepMC3/Attribute.h" #include "HepMC3/WriterAscii.h" #include "HEJ/HepMC3Interface.hh" #else #include #endif #ifdef HEJ_BUILD_WITH_HepMC3 namespace HEJ { struct HepMC3Writer::HepMC3WriterImpl: EventWriter{ HepMC3Interface HepMC3_; std::unique_ptr writer_; std::string const file_; HepMC3WriterImpl( std::string file, LHEF::HEPRUP && heprup ): HepMC3_{std::forward(heprup)}, file_{std::move(file)} {} void init_writer(){ writer_ = std::make_unique(file_, HepMC3_.run_info()); } void finish() override { if(finished()) throw std::ios_base::failure("HepMC3 writer already finished."); EventWriter::finish(); if(!writer_) // make sure that we always write something init_writer(); writer_->close(); } void write(Event const & ev) override { auto out_ev = HepMC3_(ev); //! weight names are only available after first event if(!writer_) init_writer(); writer_->write_event(out_ev); } + void set_xs_scale(const double scale) override { + HepMC3_.set_xs_scale(scale); + } + ~HepMC3WriterImpl() override { finish_or_abort(this, "HepMC3Writer"); } }; HepMC3Writer::HepMC3Writer(std::string const & file, LHEF::HEPRUP heprup): impl_{ std::make_unique(file, std::move(heprup)) } {} void HepMC3Writer::write(Event const & ev){ impl_->write(ev); } + void HepMC3Writer::set_xs_scale(const double scale) { + impl_->set_xs_scale(scale); + } + void HepMC3Writer::finish(){ impl_->finish(); } } // namespace HEJ #else // no HepMC3 namespace HEJ { struct HepMC3Writer::HepMC3WriterImpl{}; HepMC3Writer::HepMC3Writer( std::string const & /*file*/, LHEF::HEPRUP /*heprup*/ ){ throw std::invalid_argument( "Failed to create HepMC3 writer: " "HEJ 2 was built without HepMC3 support" ); } void HepMC3Writer::write(Event const & /*ev*/){ assert(false); } + + void HepMC3Writer::set_xs_scale(double /* scale */) { + assert(false); + } + void HepMC3Writer::finish(){ assert(false); } } #endif namespace HEJ { HepMC3Writer::~HepMC3Writer() = default; } diff --git a/src/LesHouchesWriter.cc b/src/LesHouchesWriter.cc index 8046ede..a9b08f7 100644 --- a/src/LesHouchesWriter.cc +++ b/src/LesHouchesWriter.cc @@ -1,130 +1,134 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/LesHouchesWriter.hh" #include #include #include #include #include #include "HEJ/Event.hh" #include "HEJ/Parameters.hh" #include "HEJ/event_types.hh" #include "HEJ/exceptions.hh" #include "HEJ/utility.hh" namespace HEJ { namespace { std::size_t to_index(event_type::EventType const type){ return type==0?0:1+std::floor(std::log2(static_cast(type))); } } LesHouchesWriter::LesHouchesWriter( std::string const & file, LHEF::HEPRUP heprup ): out_{file, std::fstream::in | std::fstream::out | std::fstream::trunc}, writer_{std::make_unique(out_)} { if(! out_.is_open()){ throw std::ios_base::failure("Failed to open " + file); }; // scientific style is needed to allow rewriting the init block out_ << std::scientific; writer_->heprup = std::move(heprup); // lhe Standard: IDWTUP (negative => weights = +/-) // IDWTUP: HEJ -> SHG/Pythia/next program // 1: weighted->unweighted, xs = mean(weight), XMAXUP given // 2: weighted->unweighted, xs = XSECUP, XMAXUP given // 3: unweighted (weight=+1)->unweighted, no additional information // 4: weighted->weighted, xs = mean(weight) // // None of these codes actually match what we want: // 1 and 4 require xs = mean(weight), which is impossible until after generation // 2 tells the SHG to unweight our events, which is wasteful // 3 claims we produce unweighted events, which is both wasteful _and_ // impossible until after generation (we don't know the maximum weight before) // // For the time being, we choose -3. If the consumer (like Pythia) assumes // weight=+1, the final weights have to be corrected by multiplying with // the original weight we provided. We are also often use NLO-PDFs which can // give negative weights, hence the native IDWTUP. // writer_->heprup.IDWTUP = -3; // always use the newest LHE version // Pythia only saves weights (hepeup.XWGTUP) for version >=2 writer_->heprup.version = LHEF::HEPRUP().version; const std::size_t max_number_types = to_index(event_type::last_type)+1; writer_->heprup.NPRUP = max_number_types; // ids of event types writer_->heprup.LPRUP.clear(); writer_->heprup.LPRUP.reserve(max_number_types); writer_->heprup.LPRUP.emplace_back(0); for(auto i=event_type::first_type+1; i<=event_type::last_type; i*=2) writer_->heprup.LPRUP.emplace_back(i); // use placeholders for unknown init block values // we can overwrite them after processing all events writer_->heprup.XSECUP = std::vector(max_number_types, 0.); writer_->heprup.XERRUP = std::vector(max_number_types, 0.); writer_->heprup.XMAXUP = std::vector(max_number_types, 0.); write_init(); } void LesHouchesWriter::write(Event const & ev){ if(finished()) { throw std::logic_error{ "Tried to write a Les Houches event after calling `finish`" }; } assert(writer_ && out_.is_open()); const double wt = ev.central().weight; writer_->hepeup = to_HEPEUP(ev, &heprup()); writer_->writeEvent(); assert(heprup().XSECUP.size() > to_index(ev.type())); heprup().XSECUP[to_index(ev.type())] += wt; heprup().XERRUP[to_index(ev.type())] += wt*wt; if(wt > heprup().XMAXUP[to_index(ev.type())]){ heprup().XMAXUP[to_index(ev.type())] = wt; } } + void LesHouchesWriter::set_xs_scale(const double scale) { + xs_scale_ = scale; + } + // this function is called after overwritting the Les Houches init block // assert that we have overwritten *exactly* the init block, // i.e. we are at the end of the file or an intact event block is next void assert_next_event_intact(std::istream & out){ ignore(out); // suppress compiler warnings if not in debug mode #ifndef NDEBUG std::string line; getline(out, line); assert(out.eof() || line.rfind("