No OneTemporary
Actions

Size

76 KB

Subscribers

None

View Options

	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<EventParameters>`
	- 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<double>`. 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<EventParameters>*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<HEJ::RNG>`
	(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 <memory>
	#include <vector>

	#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<OutputFile> 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<std::unique_ptr<EventWriter>> 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 <string>

	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 <memory>
	#include <string>

	#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<HDF5WriterImpl> 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 <array>
	#include <cstddef>

	#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<ParticleID,2> beam_particle_{};
	std::array<double,2> 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 <memory>
	#include <string>

	#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<HepMC2WriterImpl> 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 <array>
	#include <cstddef>
	#include <memory>

	#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<HepMC3::GenRunInfo> 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<ParticleID,2> beam_particle_{};
	std::array<double,2> beam_energy_{};
	std::shared_ptr<HepMC3::GenRunInfo> run_info_;
	std::size_t event_count_;
	double tot_weight_;
	double tot_weight2_;
	+ double xs_scale_ = 1.;
	std::shared_ptr<HepMC3::GenCrossSection> 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 <memory>
	#include <string>

	#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<HepMC3WriterImpl> impl_;
	};

	} // namespace HEJ
	diff --git a/include/HEJ/LesHouchesWriter.hh b/include/HEJ/LesHouchesWriter.hh
	index 077ec71..70c4d82 100644
	--- a/include/HEJ/LesHouchesWriter.hh
	+++ b/include/HEJ/LesHouchesWriter.hh
	@@ -1,61 +1,67 @@
	/** \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 <fstream>
	#include <memory>
	#include <string>

	#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();
	}

	+ void rewrite_init_and_close();
	+
	LHEF::HEPRUP & heprup(){
	return writer_->heprup;
	}

	LHEF::HEPEUP & hepeup(){
	return writer_->hepeup;
	}

	std::fstream out_;
	std::unique_ptr<LHEF::Writer> 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 <iostream>

	#include "HEJ/exceptions.hh"
	#include "HEJ/make_writer.hh"
	#include "HEJ/output_formats.hh"

	namespace HEJ {

	CombinedEventWriter::CombinedEventWriter(
	std::vector<OutputFile> 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 0e0b33b..06508c1 100644
	--- a/src/HDF5Writer.cc
	+++ b/src/HDF5Writer.cc
	@@ -1,432 +1,451 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "HEJ/HDF5Writer.hh"

	#include <cassert>

	#include "LHEF/LHEF.h"

	#include "HEJ/ConfigFlags.hh"

	#ifdef HEJ_BUILD_WITH_HDF5

	#include <cmath>
	#include <cstddef>
	#include <iterator>
	#include <type_traits>
	#include <utility>
	#include <vector>

	#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<size_t>(type)))+1;
	}

	template<typename T>
	void write_dataset(HighFive::Group & group, std::string const & name, T val) {
	using data_t = std::decay_t<T>;
	group.createDataSet<data_t>(name, DataSpace::From(val)).write(val);
	}

	template<typename T>
	void write_dataset(
	HighFive::Group & group, std::string const & name,
	std::vector<T> const & val
	) {
	using data_t = std::decay_t<T>;
	group.createDataSet<data_t>(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<int> nparticles, start, pid, id, status,
	mother1, mother2, color1, color2, npLO, npNLO;
	std::vector<double> 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<LHEF::HEPRUP>(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<double>(max_number_types, 0.);
	heprup.XERRUP = std::vector<double>(max_number_types, 0.);
	heprup.XMAXUP = std::vector<double>(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<int>("nparticles", dim, hdf5_chunk());
	events.createDataSet<int>("start", dim, hdf5_chunk());
	events.createDataSet<int>("pid", dim, hdf5_chunk());
	events.createDataSet<double>("weight", dim, hdf5_chunk());
	events.createDataSet<double>("scale", dim, hdf5_chunk());
	events.createDataSet<double>("fscale", dim, hdf5_chunk());
	events.createDataSet<double>("rscale", dim, hdf5_chunk());
	events.createDataSet<double>("aqed", dim, hdf5_chunk());
	events.createDataSet<double>("aqcd", dim, hdf5_chunk());
	events.createDataSet<double>("trials", dim, hdf5_chunk());
	events.createDataSet<int>("npLO", dim, hdf5_chunk());
	events.createDataSet<int>("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<int>("id", dim, hdf5_chunk());
	particles.createDataSet<int>("status", dim, hdf5_chunk());
	particles.createDataSet<int>("mother1", dim, hdf5_chunk());
	particles.createDataSet<int>("mother2", dim, hdf5_chunk());
	particles.createDataSet<int>("color1", dim, hdf5_chunk());
	particles.createDataSet<int>("color2", dim, hdf5_chunk());
	particles.createDataSet<double>("px", dim, hdf5_chunk());
	particles.createDataSet<double>("py", dim, hdf5_chunk());
	particles.createDataSet<double>("pz", dim, hdf5_chunk());
	particles.createDataSet<double>("e", dim, hdf5_chunk());
	particles.createDataSet<double>("m", dim, hdf5_chunk());
	particles.createDataSet<double>("lifetime", dim, hdf5_chunk());
	particles.createDataSet<double>("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 = 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<HDF5Writer::HDF5WriterImpl>(
	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 <cassert>
	#include <cmath>
	#include <utility>

	#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<HepMC::GenParticle*,2>;
	};

	template<>
	auto make_particle_ptr<2> (
	Particle const & sp, int status
	) {
	return new HepMC::GenParticle(
	to_FourVector<HepMC::FourVector>(sp),
	static_cast<int> (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<ParticleID>(heprup.IDBMUP.first),
	static_cast<ParticleID>(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<HepMC::GenParticle*,2> 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<std::size_t>(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 <utility>

	#include "HepMC/IO_GenEvent.h"

	#include "HEJ/HepMC2Interface.hh"

	#else

	#include <cassert>

	#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<HepMC2WriterImpl>(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 <cassert>
	#include <cmath>
	#include <string>
	#include <utility>
	#include <vector>

	#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<HepMC3::GenParticlePtr,2>;
	};

	template<>
	auto make_particle_ptr<3> (
	Particle const & sp, int status
	) {
	return HepMC3::make_shared<HepMC3::GenParticle>(
	to_FourVector<HepMC3::FourVector>(sp),
	static_cast<int> (sp.type),
	status
	);
	}

	template<>
	auto make_vx_ptr<3>() {
	return HepMC3::make_shared<HepMC3::GenVertex>();
	}
	} // 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<HepMC3::GenRunInfo> init_runinfo(LHEF::HEPRUP heprup){
	reset_weight_info(heprup);
	auto runinfo{ HepMC3::make_shared<HepMC3::GenRunInfo>() };

	auto hepr{ HepMC3::make_shared<HepMC3::HEPRUPAttribute>() };
	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<std::string> get_weight_names(Event const & ev){
	std::vector<std::string> 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<ParticleID>(heprup.IDBMUP.first),
	static_cast<ParticleID>(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::GenCrossSection>()}
	{
	}

	HepMC3::GenEvent HepMC3Interface::init_event(Event const & event) const {

	const std::array<HepMC3::GenParticlePtr,2> beam {
	HepMC3::make_shared<HepMC3::GenParticle>(
	HepMC3::FourVector(0,0,-beam_energy_[0],beam_energy_[0]),
	beam_particle_[0], detail_HepMC::Status::beam ),
	HepMC3::make_shared<HepMC3::GenParticle>(
	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<event.variations().size(); ++i){
	hepmc_ev.weights()[i+1] = event.variations()[i].weight;
	//! @TODO set variation specific cross section
	//! the problem is that set_cross_section overwrites everything
	}
	return hepmc_ev;
	}

	void HepMC3Interface::set_central(
	HepMC3::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<size_t>(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 <utility>

	#include "HepMC3/Attribute.h"
	#include "HepMC3/WriterAscii.h"

	#include "HEJ/HepMC3Interface.hh"

	#else

	#include <cassert>

	#endif

	#ifdef HEJ_BUILD_WITH_HepMC3

	namespace HEJ {

	struct HepMC3Writer::HepMC3WriterImpl: EventWriter{
	HepMC3Interface HepMC3_;

	std::unique_ptr<HepMC3::WriterAscii> writer_;
	std::string const file_;

	HepMC3WriterImpl(
	std::string file, LHEF::HEPRUP && heprup
	):
	HepMC3_{std::forward<LHEF::HEPRUP>(heprup)},
	file_{std::move(file)}
	{}

	void init_writer(){
	writer_ = std::make_unique<HepMC3::WriterAscii>(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<HepMC3WriterImpl>(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..5631e41 100644
	--- a/src/LesHouchesWriter.cc
	+++ b/src/LesHouchesWriter.cc
	@@ -1,130 +1,154 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "HEJ/LesHouchesWriter.hh"

	#include <cassert>
	#include <cmath>
	#include <cstddef>
	#include <utility>
	#include <vector>

	#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<std::size_t>(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<LHEF::Writer>(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<double>(max_number_types, 0.);
	writer_->heprup.XERRUP = std::vector<double>(max_number_types, 0.);
	writer_->heprup.XMAXUP = std::vector<double>(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("<event", 0) == 0);
	#endif
	}

	+ void LesHouchesWriter::rewrite_init_and_close(){
	+ assert(writer_ && out_.is_open());
	+
	+ for(auto & xs: heprup().XSECUP) {
	+ xs *= xs_scale_;
	+ }
	+ for(auto & xs_err: heprup().XERRUP) {
	+ xs_err = xs_scale_ * std::sqrt(xs_err);
	+ }
	+
	+ // replace placeholder entries
	+ out_.seekp(0);
	+ if(! out_.fail()) {
	+ write_init();
	+ assert_next_event_intact(out_);
	+ } else {
	+ std::cerr << "Failed to rewrite init block in Les Houches output\n";
	+ }
	+
	+ // close the stream before the writer destructor is called
	+ // so that it can't write a superfluous end tag
	+ out_.close();
	+ writer_.reset();
	+ }
	+
	void LesHouchesWriter::finish(){
	if(finished())
	throw std::ios_base::failure("Les Houches Output already finished");
	EventWriter::finish();

	assert(writer_);
	- for(auto & xs_err: heprup().XERRUP)
	- {
	- xs_err = std::sqrt(xs_err);
	- }
	- writer_.reset();
	-
	- out_.close();
	+ out_ << "</LesHouchesEvents>" << std::endl;
	+ rewrite_init_and_close();
	}

	LesHouchesWriter::~LesHouchesWriter(){
	finish_or_abort(this, "LesHouchesWriter");
	}

	} // namespace HEJ
	diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
	index d902f75..383320f 100644
	--- a/src/bin/HEJ.cc
	+++ b/src/bin/HEJ.cc
	@@ -1,447 +1,448 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2022
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <chrono>
	#include <iostream>
	#include <limits>
	#include <memory>
	#include <numeric>
	#include <optional>


	#include "yaml-cpp/yaml.h"

	#include "fastjet/ClusterSequence.hh"

	#include "HEJ/CombinedEventWriter.hh"
	#include "HEJ/Config.hh"
	#include "HEJ/CrossSectionAccumulator.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/EventReader.hh"
	#include "HEJ/BufferedEventReader.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/filesystem.hh"
	#include "HEJ/get_analysis.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/ProgressBar.hh"
	#include "HEJ/stream.hh"
	#include "HEJ/Unweighter.hh"
	#include "HEJ/Version.hh"
	#include "HEJ/YAMLreader.hh"

	HEJ::Config load_config(char const * filename){
	try{
	return HEJ::load_config(filename);
	}
	catch(std::exception const & exc){
	std::cerr << "Error: " << exc.what() << '\n';
	std::exit(EXIT_FAILURE);
	}
	}

	std::vector<std::unique_ptr<HEJ::Analysis>> get_analyses(
	std::vector<YAML::Node> const & parameters, LHEF::HEPRUP const & heprup
	){
	try{
	return HEJ::get_analyses(parameters, heprup);
	}
	catch(std::exception const & exc){
	std::cerr << "Failed to load analysis: " << exc.what() << '\n';
	std::exit(EXIT_FAILURE);
	}
	}

	std::string time_to_string(const time_t time){
	char s[30];
	struct tm * p = localtime(&time);
	strftime(s, 30, "%a %b %d %Y %H:%M:%S", p);
	return s;
	}

	void check_one_parton_per_jet(HEJ::Event const & ev) {
	const std::size_t npartons = std::count_if(
	ev.outgoing().begin(), ev.outgoing().end(),
	[](HEJ::Particle const & p) { return is_parton(p); }
	);
	if(ev.jets().size() < npartons) {
	throw std::invalid_argument{
	"Number of partons (" + std::to_string(npartons) + ") in input event"
	" does not match number of jets (" + std::to_string(ev.jets().size())+ ")"
	};
	}
	}

	HEJ::Event to_event(
	LHEF::HEPEUP const & hepeup,
	HEJ::JetParameters const & fixed_order_jets,
	HEJ::EWConstants const & ew_parameters,
	const double off_shell_tolerance
	) {
	HEJ::Event::EventData event_data{hepeup};
	event_data.reconstruct_intermediate(ew_parameters);
	event_data.repair_momenta(off_shell_tolerance);

	HEJ::Event ev{
	std::move(event_data).cluster(
	fixed_order_jets.def, fixed_order_jets.min_pt
	)
	};
	check_one_parton_per_jet(ev);
	return ev;
	}

	void unweight(
	HEJ::Unweighter & unweighter,
	HEJ::WeightType weight_type,
	std::vector<HEJ::Event> & events,
	HEJ::RNG & ran
	) {
	if(weight_type == HEJ::WeightType::unweighted_resum){
	unweighter.set_cut_to_maxwt(events);
	}
	events.erase(
	unweighter.unweight(begin(events), end(events), ran),
	end(events)
	);
	}

	// peek up to nevents events from reader
	std::vector<LHEF::HEPEUP> peek_events(
	HEJ::BufferedEventReader & reader,
	const int nevents
	) {
	std::vector<LHEF::HEPEUP> events;
	while(
	static_cast<int>(events.size()) < nevents
	&& reader.read_event()
	) {
	events.emplace_back(reader.hepeup());
	}
	// put everything back into the reader
	for(auto it = rbegin(events); it != rend(events); ++it) {
	reader.emplace(*it);
	}
	return events;
	}

	void append_resummed_events(
	std::vector<HEJ::Event> & resummation_events,
	HEJ::EventReweighter & reweighter,
	LHEF::HEPEUP const & hepeup,
	const size_t trials,
	HEJ::JetParameters const & fixed_order_jets,
	HEJ::EWConstants const & ew_parameters,
	const double off_shell_tolerance
	) {
	const HEJ::Event FO_event = to_event(
	hepeup,
	fixed_order_jets,
	ew_parameters,
	off_shell_tolerance
	);
	if(reweighter.treatment(FO_event.type()) != HEJ::EventTreatment::reweight) {
	return;
	}
	const auto resummed = reweighter.reweight(FO_event, trials);
	resummation_events.insert(
	end(resummation_events),
	begin(resummed), end(resummed)
	);
	}

	std::optional<HEJ::ProgressBar<std::size_t>> make_progress_bar(
	std::optional<std::size_t> n_input_events,
	char const * event_file
	) {
	if(n_input_events) {
	return HEJ::ProgressBar{std::cout, *n_input_events};
	}
	if(HEJ::is_regular(event_file)) {
	auto t_reader = HEJ::make_reader(event_file);
	std::size_t n_input_events = 0;
	while(t_reader->read_event()) ++n_input_events;
	return HEJ::ProgressBar{std::cout, n_input_events};
	}
	return {};
	}

	void train(
	HEJ::Unweighter & unweighter,
	HEJ::BufferedEventReader & reader,
	HEJ::EventReweighter & reweighter,
	const size_t total_trials,
	const double max_dev,
	double reweight_factor,
	HEJ::JetParameters const & fixed_order_jets,
	HEJ::EWConstants const & ew_parameters,
	const double off_shell_tolerance
	) {
	std::cout << "Reading up to " << total_trials << " training events...\n";
	auto FO_events = peek_events(reader, total_trials);
	if(FO_events.empty()) {
	throw std::runtime_error{
	"No events generated to calibrate the unweighting weight!"
	"Please increase the number \"trials\" or deactivate the unweighting."
	};
	}
	const size_t trials = total_trials/FO_events.size();
	// adjust reweight factor so that the overall normalisation
	// is the same as in the full run
	reweight_factor *= trials;
	for(auto & hepeup: FO_events) {
	hepeup.XWGTUP *= reweight_factor;
	}
	std::cout << "Training unweighter with "
	<< trials << '*' << FO_events.size() << " events\n";

	auto progress = HEJ::ProgressBar<int>{
	std::cout, static_cast<int>(FO_events.size())
	};
	std::vector<HEJ::Event> resummation_events;
	for(auto const & hepeup: FO_events) {
	append_resummed_events(
	resummation_events,
	reweighter, hepeup, trials, fixed_order_jets,
	ew_parameters,
	off_shell_tolerance
	);
	++progress;
	}

	unweighter.set_cut_to_peakwt(resummation_events, max_dev);
	std::cout << "\nUnweighting events with weight up to "
	<< unweighter.get_cut() << '\n';
	}

	int main(int argn, char** argv) {
	using clock = std::chrono::system_clock;

	if (argn != 3) {
	std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n";
	return EXIT_FAILURE;
	}

	const auto start_time = clock::now();
	{
	std::cout << "Starting " << HEJ::Version::package_name_full()
	<< ", revision " << HEJ::Version::revision() << " ("
	<< time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
	}
	fastjet::ClusterSequence::print_banner();

	// read configuration
	const HEJ::Config config = load_config(argv[1]);
	auto reader = HEJ::make_reader(argv[2]);
	assert(reader);

	auto heprup{ reader->heprup() };
	heprup.generators.emplace_back(LHEF::XMLTag{});
	heprup.generators.back().name = HEJ::Version::package_name();
	heprup.generators.back().version = HEJ::Version::String();

	auto analyses = get_analyses( config.analyses_parameters, heprup );
	assert(analyses.empty() \|\| analyses.front() != nullptr);

	HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};

	double global_reweight = 1.;
	auto const & max_events = config.max_events;

	// if we need the event number:
	if(std::abs(heprup.IDWTUP) == 4 \|\| std::abs(heprup.IDWTUP) == 1 \|\| max_events){
	// try to read from LHE head
	auto input_events{reader->number_events()};
	if(!input_events) {
	// else count manually
	auto t_reader = HEJ::make_reader(argv[2]);
	input_events = 0;
	while(t_reader->read_event()) ++(*input_events);
	}
	if(std::abs(heprup.IDWTUP) == 4 \|\| std::abs(heprup.IDWTUP) == 1){
	// IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs
	std::cout << "Found IDWTUP " << heprup.IDWTUP << ": "
	<< "assuming \"cross section = average weight\".\n"
	<< "converting to \"cross section = sum of weights\" ";
	global_reweight /= *input_events;
	}
	if(max_events && (input_events > max_events)){
	// maximal number of events given
	global_reweight = input_events/static_cast<double>(*max_events);
	std::cout << "Processing " << *max_events
	<< " out of " << *input_events << " events\n";
	}
	}

	HEJ::ScaleGenerator scale_gen{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	};
	std::shared_ptr<HEJ::RNG> ran{
	HEJ::make_RNG(config.rng.name, config.rng.seed)};
	assert(ran != nullptr);
	HEJ::EventReweighter hej{
	reader->heprup(),
	std::move(scale_gen),
	to_EventReweighterConfig(config),
	ran
	};

	std::optional<HEJ::Unweighter> unweighter{};
	if(config.weight_type != HEJ::WeightType::weighted) {
	unweighter = HEJ::Unweighter{};
	}
	if(config.weight_type == HEJ::WeightType::partially_unweighted) {
	HEJ::BufferedEventReader buffered_reader{std::move(reader)};
	assert(config.unweight_config);
	train(
	*unweighter,
	buffered_reader,
	hej,
	config.unweight_config->trials,
	config.unweight_config->max_dev,
	global_reweight/config.trials,
	config.fixed_order_jets,
	config.ew_parameters,
	config.off_shell_tolerance
	);
	reader = std::make_unique<HEJ::BufferedEventReader>(
	std::move(buffered_reader)
	);
	}

	// status infos & eye candy
	if (config.nlo.enabled){
	std::cout << "HEJ@NLO Truncation Enabled for NLO order: "
	<< config.nlo.nj << std::endl;
	}
	size_t nevent = 0;
	std::array<int, HEJ::event_type::last_type + 1>
	nevent_type{0}, nfailed_type{0};

	auto progress = make_progress_bar(reader->number_events(), argv[2]);
	if(!progress) {
	std::cout
	<< "Cannot determine number of input events. "
	"Disabling progress bar"
	<< std::endl;
	}

	HEJ::CrossSectionAccumulator xs;
	std::map<HEJ::StatusCode, int> status_counter;
	size_t total_trials = 0;
	size_t total_resum = 0;

	// Loop over the events in the input file
	while(reader->read_event() && (!max_events \|\| nevent < *max_events) ){
	++nevent;

	// reweight events so that the total cross section is conserved
	auto hepeup = reader->hepeup();
	hepeup.XWGTUP *= global_reweight;

	const auto FO_event = to_event(
	hepeup,
	config.fixed_order_jets,
	config.ew_parameters,
	config.off_shell_tolerance
	);
	if(FO_event.central().weight == 0) {
	static const bool warned_once = [argv,nevent](){
	std::cerr
	<< "WARNING: event number " << nevent
	<< " in " << argv[2] << " has zero weight. "
	"Ignoring this and all further events with vanishing weight.\n";
	return true;
	}();
	(void) warned_once; // shut up compiler warnings
	continue;
	}

	auto resummed_events{ hej.reweight(FO_event, config.trials) };

	// some bookkeeping
	for(auto const & s: hej.status())
	++status_counter[s];
	total_trials+=hej.status().size();
	++nevent_type[FO_event.type()];

	if(resummed_events.empty()) ++nfailed_type[FO_event.type()];

	if(unweighter) {
	unweight(unweighter, config.weight_type, resummed_events, ran);
	}

	// analysis
	for(auto & ev: resummed_events){
	//TODO: move pass_cuts to after phase space point generation
	bool passed = analyses.empty();
	for(auto const & analysis: analyses){
	if(analysis->pass_cuts(ev, FO_event)){
	passed = true;
	analysis->fill(ev, FO_event);
	};
	}
	if(passed){
	+ writer.set_xs_scale(reader->scalefactor());
	writer.write(ev);
	} else {
	ev.parameters()*=0; // do not use discarded events afterwards
	}
	}
	xs.fill_correlated(resummed_events);
	total_resum += resummed_events.size();
	if(progress) ++*progress;
	} // main event loop
	std::cout << '\n';
	// scale to account for num_trials of sherpa
	if(const double scalefactor = reader->scalefactor(); scalefactor != 1.) {
	for(auto const & analysis: analyses){
	analysis->scale(scalefactor);
	}
	xs.scale(scalefactor);
	}

	for(auto const & analysis: analyses){
	analysis->finalise();
	}
	writer.finish();

	using namespace HEJ::event_type;
	std::cout<< "Events processed: " << nevent << " (" << total_resum << " resummed)"<< '\n';
	std::cout << '\t' << name(EventType::first_type) << ": "
	<< nevent_type[EventType::first_type]
	<< ", failed to reconstruct " << nfailed_type[EventType::first_type]
	<< '\n';
	for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){
	std::cout << '\t' << name(static_cast<EventType>(i)) << ": "
	<< nevent_type[i]
	<< ", failed to reconstruct " << nfailed_type[i]
	<< '\n';
	}

	std::cout << '\n' << xs << '\n';

	std::cout << "Generation statistic: "
	<< status_counter[HEJ::StatusCode::good] << "/" << total_trials
	<< " trials successful.\n";
	for(auto && entry: status_counter){
	const double fraction = static_cast<double>(entry.second)/total_trials;
	const int percent = std::round(100*fraction);
	std::cout << std::left << std::setw(17) << (to_string(entry.first) + ":")
	<< " [";
	for(int i = 0; i < percent/2; ++i) std::cout << '#';
	for(int i = percent/2; i < 50; ++i) std::cout << ' ';
	std::cout << "] " <<std::setw(2)<<std::right<< percent << "%\n";
	}

	std::chrono::duration<double> run_time = (clock::now() - start_time);
	std::cout << "\nFinished " << HEJ::Version::package_name() << " at "
	<< time_to_string(clock::to_time_t(clock::now()))
	<< "\n=> Runtime: " << run_time.count() << " sec ("
	<< nevent/run_time.count() << " Events/sec).\n";

	return EXIT_SUCCESS;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 6:38 PM (5 h, 25 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4006916
Default Alt Text: (76 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions