No OneTemporary
Actions

Size

47 KB

Subscribers

None

View Options

	diff --git a/CMakeLists.txt b/CMakeLists.txt
	index 063b160..d572ae2 100644
	--- a/CMakeLists.txt
	+++ b/CMakeLists.txt
	@@ -1,396 +1,405 @@
	cmake_minimum_required(VERSION 3.1 FATAL_ERROR)

	set(CMAKE_LEGACY_CYGWIN_WIN32 0)
	set(CMAKE_EXPORT_COMPILE_COMMANDS ON)

	project("HEJ" VERSION 2.0.5 LANGUAGES C CXX)

	# Set a default build type if none was specified
	set(default_build_type "RelWithDebInfo")

	if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
	message(STATUS "Setting build type to '${default_build_type}' as none was specified.")
	set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
	STRING "Choose the type of build." FORCE)
	# Set the possible values of build type for cmake-gui
	set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
	"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
	endif()

	## Flags for the compiler. No warning allowed.
	if (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra")
	elseif (MSVC)
	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4 /WX /EHsc")
	endif()
	set(CMAKE_CXX_STANDARD_REQUIRED ON)
	set(CMAKE_CXX_STANDARD 14)

	## Create Version
	set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")

	# Get the latest abbreviated commit hash of the working branch
	execute_process(
	COMMAND git rev-parse HEAD
	WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
	OUTPUT_VARIABLE PROJECT_GIT_REVISION
	OUTPUT_STRIP_TRAILING_WHITESPACE
	)
	# Get the current working branch
	execute_process(
	COMMAND git rev-parse --abbrev-ref HEAD
	WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
	OUTPUT_VARIABLE PROJECT_GIT_BRANCH
	OUTPUT_STRIP_TRAILING_WHITESPACE
	)

	## target directories for install
	set(INSTALL_INCLUDE_DIR_BASE include)
	set(INSTALL_INCLUDE_DIR ${INSTALL_INCLUDE_DIR_BASE}/HEJ)
	set(INSTALL_BIN_DIR bin)
	set(INSTALL_LIB_DIR lib)
	set(INSTALL_CONFIG_DIR lib/cmake/HEJ)

	## Template files
	configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in
	${PROJECT_BINARY_DIR}/include/HEJ/Version.hh @ONLY )

	configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/HEJ-config.cc.in
	${PROJECT_BINARY_DIR}/src/bin/HEJ-config.cc @ONLY )

	# Generate CMake config file
	include(CMakePackageConfigHelpers)
	configure_package_config_file(
	cmake/Templates/hej-config.cmake.in
	${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config.cmake
	INSTALL_DESTINATION ${INSTALL_CONFIG_DIR}
	PATH_VARS INSTALL_INCLUDE_DIR_BASE INSTALL_LIB_DIR
	)
	write_basic_package_version_file(
	${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config-version.cmake
	COMPATIBILITY SameMajorVersion
	)
	install(
	FILES ${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config.cmake
	${CMAKE_CURRENT_BINARY_DIR}/${INSTALL_CONFIG_DIR}/hej-config-version.cmake
	DESTINATION ${INSTALL_CONFIG_DIR})

	## Add directories and find dependences
	include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include ${PROJECT_BINARY_DIR}/include)

	find_package(fastjet REQUIRED)
	include_directories(${fastjet_INCLUDE_DIRS})
	find_package(CLHEP 2.3 REQUIRED)
	include_directories(${CLHEP_INCLUDE_DIRS})
	find_package(LHAPDF REQUIRED)
	include_directories(${LHAPDF_INCLUDE_DIRS})
	## Amend unintuitive behaviour of FindBoost.cmake
	## Priority of BOOST_ROOT over BOOSTROOT matches FindBoost.cmake
	## at least for cmake 3.12
	if(DEFINED BOOST_ROOT)
	message("BOOST_ROOT set - only looking for Boost in ${BOOST_ROOT}")
	set(Boost_NO_BOOST_CMAKE ON)
	elseif(DEFINED BOOSTROOT)
	message("BOOSTROOT set - only looking for Boost in ${BOOSTROOT}")
	set(Boost_NO_BOOST_CMAKE ON)
	endif()
	find_package(Boost REQUIRED COMPONENTS iostreams)
	include_directories(${Boost_INCLUDE_DIRS})
	find_package(yaml-cpp) # requiring yaml does not work with fedora
	include_directories(${YAML_CPP_INCLUDE_DIR})
	if(${EXCLUDE_HepMC})
	message(STATUS "Skipping HepMC")
	# avoid "unused variable" warning if EXCLUDE_rivet is set by user
	set(EXCLUDE_rivet TRUE)
	else()
	find_package(HepMC 2)
	endif()
	if(${HepMC_FOUND})
	set(
	CMAKE_CXX_FLAGS
	"${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_HepMC_VERSION=${HepMC_VERSION_MAJOR}"
	)
	include_directories(${HepMC_INCLUDE_DIRS})
	if(${EXCLUDE_rivet})
	message(STATUS "Skipping rivet")
	else()
	find_package(rivet)
	endif()
	if(${rivet_FOUND})
	include_directories(${rivet_INCLUDE_DIRS})
	set(
	CMAKE_CXX_FLAGS
	"${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_RIVET"
	)
	endif()
	endif()
	if(${EXCLUDE_QCDloop})
	message(STATUS "Skipping QCDloop")
	else()
	find_package(QCDloop 2)
	endif()
	if(${QCDloop_FOUND})
	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_QCDLOOP")
	include_directories(SYSTEM ${QCDloop_INCLUDE_DIRS})
	endif()
	if(${EXCLUDE_HighFive})
	message(STATUS "Skipping HighFive")
	else()
	find_package(HighFive QUIET)
	find_package_handle_standard_args( HighFive CONFIG_MODE )
	endif()
	if(${HighFive_FOUND})
	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHEJ_BUILD_WITH_HDF5")
	include_directories($<TARGET_PROPERTY:HighFive,INTERFACE_INCLUDE_DIRECTORIES>)
	endif()

	add_subdirectory(src)

	## define executable
	add_executable(HEJ src/bin/HEJ.cc)

	## link libraries
	target_link_libraries(HEJ hejlib)

	add_executable(HEJ-config src/bin/HEJ-config.cc)

	file(GLOB hej_headers
	${CMAKE_CURRENT_SOURCE_DIR}/include/HEJ/*.hh
	${PROJECT_BINARY_DIR}/include/HEJ/*.hh
	)
	file(GLOB hej_detail_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/HEJ/detail/*.hh)
	file(GLOB lhef_headers ${CMAKE_CURRENT_SOURCE_DIR}/include/LHEF/*.h)
	install(FILES ${hej_headers} DESTINATION ${INSTALL_INCLUDE_DIR})
	install(FILES ${hej_detail_headers} DESTINATION ${INSTALL_INCLUDE_DIR}/detail/)
	install(FILES ${lhef_headers} DESTINATION ${INSTALL_INCLUDE_DIR_BASE}/LHEF/)
	install(TARGETS HEJ HEJ-config DESTINATION ${INSTALL_BIN_DIR})

	## tests
	enable_testing()
	set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}/t")
	set(tst_ME_data_dir "${tst_dir}/ME_data")

	# test event classification
	add_executable(test_classify ${tst_dir}/test_classify.cc)
	target_link_libraries(test_classify hejlib)
	add_test(
	NAME t_classify
	COMMAND test_classify
	)
	add_executable(test_classify_ref ${tst_dir}/test_classify_ref.cc)
	target_link_libraries(test_classify_ref hejlib)
	add_test(
	NAME t_classify_ref
	COMMAND test_classify_ref ${tst_dir}/classify_ref ${tst_dir}/classify.lhe.gz
	)
	add_test(
	NAME t_classify_ref_4j
	COMMAND test_classify_ref ${tst_dir}/classify_ref_4j ${tst_dir}/4j.lhe.gz
	)

	# test phase space point
	add_executable(test_psp ${tst_dir}/test_psp.cc)
	target_link_libraries(test_psp hejlib)
	add_test(
	NAME t_psp
	COMMAND test_psp ${tst_dir}/psp_gen.lhe.gz
	)

	# test importing scales
	add_library(scales SHARED ${tst_dir}/scales.cc)
	add_executable(test_scale_import ${tst_dir}/test_scale_import)
	target_link_libraries(test_scale_import hejlib)
	add_test(
	NAME t_scale_import
	COMMAND test_scale_import ${tst_dir}/jet_config_with_import.yml
	)

	# test scale arithmetic (e.g. 2*H_T/4)
	add_executable(test_scale_arithmetics ${tst_dir}/test_scale_arithmetics)
	target_link_libraries(test_scale_arithmetics hejlib)
	add_test(
	NAME t_scale_arithmetics
	COMMAND test_scale_arithmetics ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz
	)

	# test "ParameterDescription"
	add_executable(test_descriptions ${tst_dir}/test_descriptions)
	target_link_libraries(test_descriptions hejlib)
	add_test(
	NAME t_descriptions
	COMMAND test_descriptions
	)

	# test "EventParameters*Weight"
	add_executable(test_parameters ${tst_dir}/test_parameters)
	target_link_libraries(test_parameters hejlib)
	add_test(
	NAME test_parameters
	COMMAND test_parameters
	)

	# test colour generation
	add_executable(test_colours ${tst_dir}/test_colours)
	target_link_libraries(test_colours hejlib)
	add_test(
	NAME t_colour_flow
	COMMAND test_colours
	)

	# test matrix elements
	add_executable(test_ME_generic ${tst_dir}/test_ME_generic.cc)
	target_link_libraries(test_ME_generic hejlib)
	add_test(
	NAME t_ME_j
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz
	)
	add_test(
	NAME t_ME_h
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
	)
	if(${QCDloop_FOUND})
	add_test(
	NAME t_ME_h_mt
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_h_mt_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
	)
	add_test(
	NAME t_ME_h_mtmb
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_h_mtmb_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz
	)
	endif()
	add_test(
	NAME t_ME_w_FKL
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_tree.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz
	)
	add_test(
	NAME t_ME_w_FKL_virt
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_virt.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz
	)
	add_test(
	NAME t_ME_Wp
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wp.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz
	)
	add_test(
	NAME t_ME_Wm
	COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wm.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz
	)

	# test main executable
	if(${HepMC_FOUND})
	file(READ "${tst_dir}/jet_config.yml" config)
	file(WRITE "${tst_dir}/jet_config_withHepMC.yml" "${config} - tst.hepmc")
	if(${rivet_FOUND})
	file(READ "${tst_dir}/jet_config_withHepMC.yml" config)
	file(WRITE "${tst_dir}/jet_config_withRivet.yml" "${config}\n\nanalysis:\n rivet: MC_XS\n output: tst")
	add_test(
	NAME t_main
	COMMAND HEJ ${tst_dir}/jet_config_withRivet.yml ${tst_dir}/2j.lhe.gz
	)
	else() # no rivet
	add_test(
	NAME t_main
	COMMAND HEJ ${tst_dir}/jet_config_withHepMC.yml ${tst_dir}/2j.lhe.gz
	)
	endif()
	if(${HepMC_VERSION_MAJOR} GREATER 2)
	add_executable(check_hepmc ${tst_dir}/check_hepmc.cc)
	target_link_libraries(check_hepmc hejlib)
	# check that HepMC output is correct
	add_test(
	NAME t_hepmc
	COMMAND check_hepmc tst.hepmc
	)
	endif()
	else() # no HepMC
	add_test(
	NAME t_main
	COMMAND HEJ ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz
	)
	endif()
	# check that LHEF output is correct
	add_executable(check_lhe ${tst_dir}/check_lhe.cc)
	target_link_libraries(check_lhe hejlib)
	add_test(
	NAME t_lhe
	COMMAND check_lhe tst.lhe
	)
	# check HDF5 reader
	if(${HighFive_FOUND})
	add_executable(test_hdf5 ${tst_dir}/test_hdf5.cc)
	target_link_libraries(test_hdf5 hejlib)
	add_test(
	NAME t_hdf5
	COMMAND test_hdf5 ${tst_dir}/Wm9-g4-repack.hdf5
	)
	endif()
	+# check rivet interface
	+if(${RIVET_FOUND})
	+ add_executable(check_rivet ${tst_dir}/check_rivet.cc)
	+ target_link_libraries(check_rivet hejlib)
	+ add_test(
	+ NAME t_rivet
	+ COMMAND check_rivet
	+ )
	+endif()

	# test boson reconstruction
	add_executable(cmp_events ${tst_dir}/cmp_events.cc)
	target_link_libraries(cmp_events hejlib)
	add_test(
	NAME t_epnu_2j_noW
	COMMAND cmp_events ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz ${tst_dir}/epnu2jLOFKL_unweight_noW.lhe.gz
	)

	# test resummed result
	add_executable(check_res ${tst_dir}/check_res.cc)
	target_link_libraries(check_res hejlib)

	if(${TEST_ALL}) # deactivate long tests by default
	add_test(
	NAME t_2j
	COMMAND check_res ${tst_dir}/2j.lhe.gz 3.49391e+07 419684
	)
	add_test(
	NAME t_3j
	COMMAND check_res ${tst_dir}/3j.lhe.gz 2.37902e+06 25746.6
	)
	add_test(
	NAME t_4j
	COMMAND check_res ${tst_dir}/4j.lhe.gz 603713 72822.6
	)
	add_test(
	NAME t_h_3j
	COMMAND check_res ${tst_dir}/h_3j.lhe.gz 0.821622 0.0220334
	)
	add_test(
	NAME t_h_3j_unof
	COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0127362 0.000271555 unof
	)
	add_test(
	NAME t_h_3j_unob
	COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0130615 0.000224793 unob
	)
	add_test(
	NAME t_epnu_2j
	COMMAND check_res ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz 262.7 3
	)
	add_test(
	NAME t_MGepnu_3j
	COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.tar.gz 38.9512 1
	)
	add_test(
	NAME t_MGemnubar_3j
	COMMAND check_res ${tst_dir}/MGemnubar3j_unweighted.tar.gz 24.1575 1
	)
	add_test(
	NAME t_MGepnu_3j_splitf
	COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.tar.gz 2.91995 0.0463182 splitf
	)
	add_test(
	NAME t_MGepnu_3j_splitb
	COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.tar.gz 3.40708 0.0550975 splitb
	)
	add_test(
	NAME t_MGepnu_4j
	COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.tar.gz 10.2542 0.135106
	)
	add_test(
	NAME t_MGemnubar_4j
	COMMAND check_res ${tst_dir}/MGemnubar4j_unweighted.tar.gz 5.57909 0.0300496
	)
	add_test(
	NAME t_MGepnu_4j_qqxmid
	COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.tar.gz 0.732084 0.005 qqxmid
	)
	endif()
	diff --git a/Changes-API.md b/Changes-API.md
	index c3c285c..de1da09 100644
	--- a/Changes-API.md
	+++ b/Changes-API.md
	@@ -1,64 +1,66 @@
	# 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`, where the main features are documented.

	## Version 2.X

	### 2.X.0

	* Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public
	* New class `CrossSectionAccumulator` to keep track of Cross Section of the
	different subproccess
	* 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 HEJ Version 2.3.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
	* 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 HEJ Version
	2.3.0
	- Removed `Event.unclustered()` function
	- Added new member function `Events.parameters()`, to directly access
	(underlying) `Parameters<EventParameters>`
	- New member functions `begin_partons`, `end_partons` with aliases
	`cbegin_partons`, `cend_partons` for constant 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)
	* New abstact `EventReader` class, as base for reading events from files
	- Moved LHE file reader to `HEJ::LesHouchesReader`
	- New `HEJ::HDF5Reader` to read `hdf5` files
	+* New function `Analysis.initialise(LHEF::HEPRUP const &)` to pass `HEPRUP` to
	+ the analysis.

	## Version 2.0

	### 2.0.5

	* no further changes to API

	### 2.0.4

	* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`

	### 2.0.3

	* no further changes to API

	### 2.0.2

	* no further changes to API

	### 2.0.1

	* no further changes to API
	diff --git a/include/HEJ/Analysis.hh b/include/HEJ/Analysis.hh
	index f00d7fc..999702a 100644
	--- a/include/HEJ/Analysis.hh
	+++ b/include/HEJ/Analysis.hh
	@@ -1,47 +1,57 @@
	/** \file
	* \brief Header file for the Analysis interface
	*
	* This header contains declarations that faciliate creating custom analyses
	* to be used with HEJ 2.
	* \todo link to user documentation
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	+namespace LHEF {
	+ class HEPRUP;
	+}
	+
	//! Main HEJ 2 Namespace
	namespace HEJ{
	class Event;

	//! Analysis base class
	/**
	* This is the interface that all analyses should implement,
	* i.e. all custom analyses have to be derived from this struct.
	*/
	struct Analysis{
	+ //! Optional function to read general setup from LHEF::HEPRUP
	+ /**
	+ * This function should be called at the start, directly after the
	+ * constructor.
	+ */
	+ virtual void initialise(LHEF::HEPRUP const &){};
	//! Fill event into analysis (e.g. to histograms)
	/**
	* @param res_event The event in resummation phase space
	* @param FO_event The original fixed-order event
	*/
	virtual void fill(Event const & res_event, Event const & FO_event) = 0;
	//! Decide whether an event passes the cuts
	/**
	* @param res_event The event in resummation phase space
	* @param FO_event The original fixed-order event
	* @returns Whether the event passes all cuts
	*/
	virtual bool pass_cuts(Event const & res_event, Event const & FO_event) = 0;
	//! Finalise analysis
	/**
	* This function is called after all events have been processed and
	* can be used for example to print out or save the results.
	*/
	virtual void finalise() = 0;

	virtual ~Analysis() = default;
	};

	}
	diff --git a/include/HEJ/HepMCInterface.hh b/include/HEJ/HepMCInterface.hh
	index 94c788e..d32bebd 100644
	--- a/include/HEJ/HepMCInterface.hh
	+++ b/include/HEJ/HepMCInterface.hh
	@@ -1,72 +1,81 @@
	/** \file
	* \brief Header file for the HepMCInterface
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	+#include <array>
	#include <sys/types.h>
	#include <vector>

	+#include "HEJ/PDG_codes.hh"
	+
	namespace HepMC{
	class GenCrossSection;
	class GenEvent;
	}

	+namespace LHEF{
	+ class HEPRUP;
	+}
	+
	namespace HEJ{
	class Event;
	class EventParameters;
	//! 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 HepMCInterface{
	public:
	- HepMCInterface();
	+ HepMCInterface(LHEF::HEPRUP const &);
	/**
	* \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, ssize_t weight_index = -1);
	/**
	* \brief initialise the event kinematics (everything but the weights)
	*
	* \param event Event to convert
	* \param weight_index optional selection of specific weight
	* (negative value gives central weight)
	*/
	HepMC::GenEvent init_kinematics(Event const & event);
	/**
	* \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")
	*/
	void set_central(HepMC::GenEvent & out_ev, Event const & event,
	ssize_t weight_index = -1);
	/**
	* \brief Add the event \p variations to \p out_ev
	*/
	void add_variation(HepMC::GenEvent & out_ev,
	std::vector<EventParameters> const & variations);
	private:
	+ std::array<ParticleID,2> beam_particle_;
	+ std::array<double,2> beam_energy_;
	size_t event_count_;
	double tot_weight_;
	double tot_weight2_;
	HepMC::GenCrossSection cross_section() const;

	};
	}
	diff --git a/include/HEJ/RivetAnalysis.hh b/include/HEJ/RivetAnalysis.hh
	index cac5738..07e22e5 100644
	--- a/include/HEJ/RivetAnalysis.hh
	+++ b/include/HEJ/RivetAnalysis.hh
	@@ -1,68 +1,74 @@
	/** \file
	* \brief HEJ 2 interface to rivet analyses
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <memory>
	#include <string>
	#include <vector>

	+#include "LHEF/LHEF.h"
	+
	#include "HEJ/Analysis.hh"
	#include "HEJ/HepMCInterface.hh"

	namespace Rivet {
	class AnalysisHandler;
	}
	namespace YAML {
	class Node;
	}

	namespace HEJ {
	/**
	* @brief Class representing a Rivet analysis
	*
	* This class inherits from Analysis and can therefore be used
	* like any other HEJ 2 analysis.
	*/
	class RivetAnalysis: public HEJ::Analysis {
	public:
	static std::unique_ptr<Analysis> create(YAML::Node const & config);

	//! Constructor
	/**
	* @param config Configuration parameters
	*
	- * config["analysis"] should be the name of a single Rivet analysis or
	+ * config["rivet"] should be the name of a single Rivet analysis or
	* a list of Rivet analyses. config["output"] is the prefix for
	* the .yoda output files.
	*/
	RivetAnalysis(YAML::Node const & config);
	+ //! Initialise common informations, given in LHEF::HEPRUP
	+ void initialise(LHEF::HEPRUP const &) override;
	//! Pass an event to the underlying Rivet analysis
	void fill(HEJ::Event const & event, HEJ::Event const &) override;
	bool pass_cuts(HEJ::Event const &, HEJ::Event const &) override
	{return true;} //< no additional cuts are applied
	void finalise() override;
	private:
	std::vector<std::string> analyses_names_;
	const std::string output_name_;
	+ LHEF::HEPRUP heprup_;

	/// struct to organise the infos per rivet run/scale setting
	struct rivet_info {
	std::unique_ptr<Rivet::AnalysisHandler> handler;
	std::string name;
	HEJ::HepMCInterface hepmc;
	};
	std::vector<rivet_info> rivet_runs_;

	/**
	* \internal
	- * @brief Calculates the scale variation from the first event for the output file
	+ * @brief Calculates the scale variation from the first event for the output
	+ * file
	*/
	void init(HEJ::Event const & event);
	bool first_event_;
	};
	}
	diff --git a/src/HepMCInterface.cc b/src/HepMCInterface.cc
	index 7e5b7fe..01cb081 100644
	--- a/src/HepMCInterface.cc
	+++ b/src/HepMCInterface.cc
	@@ -1,176 +1,208 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/HepMCInterface.hh"

	#include "HEJ/exceptions.hh"

	#ifdef HEJ_BUILD_WITH_HepMC_VERSION

	#include <math.h>
	#include <utility>

	#include "HEJ/Event.hh"
	#include "HEJ/Particle.hh"

	+#include "LHEF/LHEF.h"
	+
	#include "HepMC/GenCrossSection.h"
	#include "HepMC/GenEvent.h"
	#include "HepMC/GenParticle.h"
	#include "HepMC/GenVertex.h"

	namespace HEJ{

	namespace {
	HepMC::FourVector to_FourVector(Particle const & sp){
	return {sp.px(), sp.py(), sp.pz(), sp.E()};
	}
	+ constexpr int status_beam = 4;
	constexpr int status_in = 11;
	constexpr int status_decayed = 2;
	constexpr int status_out = 1;

	template<class HepMCClass, typename... Args>
	auto make_ptr(Args&&... args){
	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	return HepMC::make_shared<HepMCClass>(std::forward<Args>(args)...);
	#else
	return new HepMCClass(std::forward<Args>(args)...);
	#endif
	}
	} // namespace anonymous

	- HepMCInterface::HepMCInterface():
	+ HepMCInterface::HepMCInterface(LHEF::HEPRUP const & heprup):
	event_count_(0.), tot_weight_(0.), tot_weight2_(0.)
	- {}
	+ {
	+ beam_particle_[0] = static_cast<ParticleID>(heprup.IDBMUP.first);
	+ beam_particle_[1] = static_cast<ParticleID>(heprup.IDBMUP.second);
	+ beam_energy_[0] = heprup.EBMUP.first;
	+ beam_energy_[1] = heprup.EBMUP.second;
	+ }

	HepMC::GenCrossSection HepMCInterface::cross_section() const {
	HepMC::GenCrossSection xs;
	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	xs.set_cross_section(tot_weight_, sqrt(tot_weight2_), event_count_);
	/// @TODO add number of attempted events
	#else // HepMC 2
	xs.set_cross_section(tot_weight_, sqrt(tot_weight2_));
	#endif
	return xs;
	}

	HepMC::GenEvent HepMCInterface::init_kinematics(Event const & event) {
	+ #if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	+ using GenParticlePtr = HepMC::GenParticlePtr;
	+ #else
	+ using GenParticlePtr = HepMC::GenParticle*;
	+ #endif
	HepMC::GenEvent out_ev{HepMC::Units::GEV, HepMC::Units::MM};
	+
	+ //! @TODO don't hardcode beam energy
	+ #if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	+ static
	+ #endif
	+ const std::array<GenParticlePtr,2> beam {
	+ make_ptr<HepMC::GenParticle>(
	+ HepMC::FourVector(0,0,-beam_energy_[0],beam_energy_[0]),
	+ beam_particle_[0],status_beam ),
	+ make_ptr<HepMC::GenParticle>(
	+ HepMC::FourVector(0,0, beam_energy_[1],beam_energy_[1]),
	+ beam_particle_[1],status_beam )
	+ };
	+ out_ev.set_beam_particles(beam[0],beam[1]);
	+
	auto vx = make_ptr<HepMC::GenVertex>();
	- for(auto const & in: event.incoming()){
	- vx->add_particle_in(
	- make_ptr<HepMC::GenParticle>(
	- to_FourVector(in), static_cast<int>(in.type), status_in
	- )
	- );
	+ for(size_t i=0; i<event.incoming().size(); ++i){
	+ auto particle{ make_ptr<HepMC::GenParticle>(
	+ to_FourVector(event.incoming()[i]),
	+ static_cast<int>(event.incoming()[i].type), status_in
	+ ) };
	+ auto vx_beam{ make_ptr<HepMC::GenVertex>() };
	+ vx_beam->add_particle_in(beam[i]);
	+ vx_beam->add_particle_out(particle);
	+ out_ev.add_vertex(vx_beam);
	+ vx->add_particle_in(particle);
	}
	for(size_t i=0; i < event.outgoing().size(); ++i){
	auto const & out = event.outgoing()[i];
	auto particle = make_ptr<HepMC::GenParticle>(
	to_FourVector(out), static_cast<int>(out.type), status_out
	);
	const int status = event.decays().count(i)?status_decayed:status_out;
	particle->set_status(status);
	if( status == status_decayed){
	auto vx_decay = make_ptr<HepMC::GenVertex>();
	vx_decay->add_particle_in(particle);
	for( auto const & out: event.decays().at(i)){
	vx_decay->add_particle_out(
	make_ptr<HepMC::GenParticle>(
	to_FourVector(out), static_cast<int>(out.type), status_out
	)
	);
	}
	out_ev.add_vertex(vx_decay);
	}
	vx->add_particle_out(particle);
	}
	out_ev.add_vertex(vx);

	return out_ev;
	}

	void HepMCInterface::set_central(HepMC::GenEvent & out_ev, Event const & event,
	ssize_t const weight_index
	) {
	EventParameters event_param;
	if(weight_index < 0)
	event_param = event.central();
	else if ( (size_t) weight_index < event.variations().size())
	event_param = event.variations(weight_index);
	else
	throw std::invalid_argument{
	"HepMCInterface tried to access a weight outside of the variation range."
	};
	const double wt = event_param.weight;

	tot_weight_ += wt;
	tot_weight2_ += wt * wt;

	if(out_ev.weights().size() == 0){
	out_ev.weights().push_back(wt);
	} else { // central always on first
	out_ev.weights()[0] = wt;
	}

	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	out_ev.set_cross_section(
	HepMC::make_shared<HepMC::GenCrossSection>(cross_section()) );
	#else // HepMC 2
	out_ev.set_cross_section( cross_section() );
	out_ev.set_signal_process_id(event.type()+1); // "+1": conistent with lhe
	out_ev.set_event_scale(event_param.mur);
	#endif

	++event_count_;
	out_ev.set_event_number(event_count_);

	/// @TODO add alphaQCD (need function) and alphaQED
	/// @TODO output pdf (currently not avaiable from event alone)

	}

	void HepMCInterface::add_variation(HepMC::GenEvent & out_ev,
	std::vector<EventParameters> const & varis
	) {
	for(auto const & var: varis){
	out_ev.weights().push_back(var.weight);
	}
	/// @TODO add name list for weights
	}

	HepMC::GenEvent HepMCInterface::operator()(Event const & event,
	ssize_t const weight_index
	) {
	HepMC::GenEvent out_ev(init_kinematics(event));
	set_central(out_ev, event, weight_index);
	add_variation(out_ev, event.variations());
	return out_ev;
	}

	}
	#else // no HepMC => empty class
	namespace HepMC {
	class GenEvent {};
	class GenCrossSection {};
	}
	namespace HEJ{
	- HepMCInterface::HepMCInterface(){
	+ HepMCInterface::HepMCInterface(LHEF::HEPRUP const &){
	throw std::invalid_argument(
	"Failed to create HepMCInterface: "
	"HEJ 2 was built without HepMC support"
	);
	}

	HepMC::GenEvent HepMCInterface::operator()(Event const &, ssize_t)
	{return HepMC::GenEvent();}
	HepMC::GenEvent HepMCInterface::init_kinematics(Event const &)
	{return HepMC::GenEvent();}
	void HepMCInterface::add_variation(HepMC::GenEvent &,
	std::vector<EventParameters> const &){}
	void HepMCInterface::set_central(HepMC::GenEvent &, Event const &, ssize_t) {}
	HepMC::GenCrossSection HepMCInterface::cross_section() const
	{return HepMC::GenCrossSection();}
	}
	#endif
	diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc
	index c23014b..938a558 100644
	--- a/src/HepMCWriter.cc
	+++ b/src/HepMCWriter.cc
	@@ -1,148 +1,148 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/HepMCWriter.hh"

	#include <cassert>

	#include "LHEF/LHEF.h"

	#ifdef HEJ_BUILD_WITH_HepMC_VERSION

	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	#include "HepMC/LHEFAttributes.h"
	#include "HepMC/WriterAscii.h"

	#include "HEJ/Version.hh"

	#else
	#include "HepMC/IO_GenEvent.h"
	#endif

	#include <utility>

	#include "HepMC/GenParticle.h"
	#include "HepMC/GenVertex.h"

	#include "HEJ/Event.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/HepMCInterface.hh"

	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	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.};
	}
	HepMC::shared_ptr<HepMC::GenRunInfo> init_runinfo(LHEF::HEPRUP && heprup){
	reset_weight_info(heprup);
	HepMC::GenRunInfo runinfo;

	auto hepr = HepMC::make_shared<HepMC::HEPRUPAttribute>();
	hepr->heprup = heprup;
	runinfo.add_attribute(std::string("HEPRUP"), hepr);

	for (int i = 0, N = hepr->heprup.generators.size(); i < N; ++i ){
	HepMC::GenRunInfo::ToolInfo tool;
	tool.name = hepr->heprup.generators[i].name;
	tool.version = hepr->heprup.generators[i].version;
	tool.description = hepr->heprup.generators[i].contents;
	runinfo.tools().push_back(tool);
	}
	return HepMC::make_shared<HepMC::GenRunInfo>(runinfo);
	}

	} // namespace anonymous
	#endif // HepMC 3

	namespace HEJ{

	struct HepMCWriter::HepMCWriterImpl{
	HepMCInterface hepmc_;

	HepMCWriterImpl & operator=(HepMCWriterImpl const & other) = delete;
	HepMCWriterImpl(HepMCWriterImpl const & other) = delete;
	HepMCWriterImpl & operator=(HepMCWriterImpl && other) = delete;
	HepMCWriterImpl(HepMCWriterImpl && other) = delete;

	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3

	HepMC::WriterAscii writer_;

	HepMCWriterImpl(
	std::string const & file, LHEF::HEPRUP && heprup
	):
	- hepmc_(),
	+ hepmc_(heprup),
	writer_{file, init_runinfo(std::move(heprup))}
	{}

	~HepMCWriterImpl(){
	writer_.close();
	}

	#else // HepMC 2

	HepMC::IO_GenEvent writer_;

	HepMCWriterImpl(
	- std::string const & file, LHEF::HEPRUP &&
	+ std::string const & file, LHEF::HEPRUP && heprup
	):
	- hepmc_(),
	+ hepmc_(heprup),
	writer_{file}
	{}

	#endif

	void write(Event const & ev){
	auto out_ev = hepmc_(ev);
	#if HEJ_BUILD_WITH_HepMC_VERSION >= 3
	writer_.write_event(out_ev);
	#else // HepMC 2
	writer_.write_event(&out_ev);
	#endif
	}
	};

	void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
	delete p;
	}

	HepMCWriter::HepMCWriter(std::string const & file, LHEF::HEPRUP heprup):
	impl_{std::unique_ptr<HepMCWriterImpl, HepMCWriterImplDeleter>{
	new HepMCWriterImpl(file, std::move(heprup))
	}}
	{}

	void HepMCWriter::write(Event const & ev){
	impl_->write(ev);
	}
	} // namespace HEJ

	#else // no HepMC

	namespace HEJ{

	class HepMCWriter::HepMCWriterImpl{};

	HepMCWriter::HepMCWriter(std::string const &, LHEF::HEPRUP){
	throw std::invalid_argument(
	"Failed to create HepMC writer: "
	"HEJ 2 was built without HepMC support"
	);
	}

	void HepMCWriter::write(Event const &){
	assert(false);
	}

	void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
	delete p;
	}
	}
	#endif
	diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc
	index cb07fd9..440d47f 100644
	--- a/src/RivetAnalysis.cc
	+++ b/src/RivetAnalysis.cc
	@@ -1,143 +1,150 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/RivetAnalysis.hh"

	#ifdef HEJ_BUILD_WITH_RIVET

	#include <ostream>
	#include <stddef.h>

	#include "yaml-cpp/yaml.h"

	#include "Rivet/AnalysisHandler.hh"

	#include "HepMC/GenEvent.h"

	#include "HEJ/Event.hh"
	#include "HEJ/exceptions.hh"

	#endif

	namespace HEJ{
	std::unique_ptr<Analysis> RivetAnalysis::create(YAML::Node const & config){
	return std::unique_ptr<Analysis>{new RivetAnalysis{config}};
	}
	}

	#ifdef HEJ_BUILD_WITH_RIVET

	namespace HEJ {

	RivetAnalysis::RivetAnalysis(YAML::Node const & config):
	output_name_{config["output"].as<std::string>()},
	first_event_(true)
	{
	// read in analyses
	const auto & name_node = config["rivet"];
	switch(name_node.Type()){
	case YAML::NodeType::Scalar:
	analyses_names_.push_back(name_node.as<std::string>());
	break;
	case YAML::NodeType::Sequence:
	for(YAML::const_iterator it = name_node.begin(); it != name_node.end(); ++it){
	analyses_names_.push_back(it->as<std::string>());
	}
	break;
	default:
	throw std::invalid_argument{
	"No Analysis was provided to rivet. "
	"Either give an analysis or deactivate rivet."
	};
	}
	}

	namespace {
	void replace(
	std::string & str,
	char to_replace, std::string const & replacement
	) {
	for(
	auto pos = str.find(to_replace);
	pos != str.npos;
	pos = str.find(to_replace, pos)
	) {
	str.replace(pos, 1, replacement);
	pos += replacement.size();
	}
	}

	// remove "special" characters from scale name
	// so that we can more easily use it as part of a file name
	std::string sanitise_scalename(std::string scalename) {
	replace(scalename, '/', "_over_");
	replace(scalename, '*', "_times_");
	return scalename;
	}
	}

	+ void RivetAnalysis::initialise(LHEF::HEPRUP const & heprup){
	+ heprup_ = heprup;
	+ }
	+
	+ // it is a bit ugly that we can't do this directly in `initialise`
	void RivetAnalysis::init(Event const & event){
	+ rivet_runs_.reserve(event.variations().size()+1);
	rivet_runs_.push_back(
	- {std::make_unique<Rivet::AnalysisHandler>(), "", HepMCInterface()}
	+ {std::make_unique<Rivet::AnalysisHandler>(), "", HepMCInterface(heprup_)}
	);
	rivet_runs_.back().handler->addAnalyses(analyses_names_);
	if( !event.variations().empty() ){
	- rivet_runs_.reserve(event.variations().size()+1);
	for(auto const & vari : event.variations()){
	std::ostringstream name;
	name << ".Scale" << sanitise_scalename(vari.description->scale_name)
	<< "_MuR" << vari.description->mur_factor
	<< "_MuF" << vari.description->muf_factor;
	rivet_runs_.push_back(
	- {std::make_unique<Rivet::AnalysisHandler>(), name.str(), HepMCInterface()}
	+ {std::make_unique<Rivet::AnalysisHandler>(), name.str(),
	+ HepMCInterface(heprup_)}
	);
	rivet_runs_.back().handler->addAnalyses(analyses_names_);
	}
	}
	}

	void RivetAnalysis::fill(Event const & event, Event const &){
	if(first_event_){
	first_event_=false;
	init(event);
	}
	HepMC::GenEvent hepmc_kin(rivet_runs_[0].hepmc.init_kinematics(event));
	for(size_t i = 0; i < rivet_runs_.size(); ++i){
	auto & run = rivet_runs_[i];
	run.hepmc.set_central(hepmc_kin, event, i-1); // -1: first = central
	run.handler->analyze(hepmc_kin);
	}
	}

	void RivetAnalysis::finalise(){
	for(auto const & run: rivet_runs_){
	run.handler->finalize();
	run.handler->writeData(output_name_+run.name+std::string(".yoda"));
	}
	}
	} // namespace HEJ

	#else // no rivet => create empty analysis

	namespace Rivet {
	class AnalysisHandler {};
	}

	namespace HEJ {

	RivetAnalysis::RivetAnalysis(YAML::Node const &)
	{
	throw std::invalid_argument(
	"Failed to create RivetAnalysis: "
	"HEJ 2 was built without rivet support"
	);
	}

	+ void RivetAnalysis::initialise(LHEF::HEPRUP const &){}
	void RivetAnalysis::init(Event const &){}
	void RivetAnalysis::fill(Event const &, Event const &){}
	void RivetAnalysis::finalise(){}
	} // namespace HEJ

	#endif
	diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
	index 6f0a1a4..b4f1a12 100644
	--- a/src/bin/HEJ.cc
	+++ b/src/bin/HEJ.cc
	@@ -1,228 +1,230 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <chrono>
	#include <iostream>
	#include <limits>
	#include <memory>
	#include <numeric>

	#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/EventReweighter.hh"
	#include "HEJ/get_analysis.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/ProgressBar.hh"
	#include "HEJ/stream.hh"
	#include "HEJ/Version.hh"
	#include "HEJ/YAMLreader.hh"

	int event_number(std::string const & record){
	size_t start = record.rfind("Number of Events");
	start = record.find_first_of("123456789", start);
	if(start == std::string::npos) {
	throw std::invalid_argument("no event number record found");
	}
	const size_t end = record.find_first_not_of("0123456789", start);
	return std::stoi(record.substr(start, end - start));
	}

	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::unique_ptr<HEJ::Analysis> get_analysis(
	YAML::Node const & parameters
	){
	try{
	return HEJ::get_analysis(parameters);
	}
	catch(std::exception const & exc){
	std::cerr << "Failed to load analysis: " << exc.what() << '\n';
	std::exit(EXIT_FAILURE);
	}
	}

	// unique_ptr is a workaround:
	// HEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0
	std::unique_ptr<HEJ::ProgressBar<double>> make_progress_bar(
	std::vector<double> const & xs
	) {
	if(xs.empty()) return {};
	const double Born_xs = std::accumulate(begin(xs), end(xs), 0.);
	return std::make_unique<HEJ::ProgressBar<double>>(std::cout, Born_xs);
	}

	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;
	}

	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);

	std::unique_ptr<HEJ::Analysis> analysis = get_analysis(
	config.analysis_parameters
	);
	assert(analysis != nullptr);

	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();
	+
	+ analysis->initialise(heprup);
	HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};

	double global_reweight = 1.;
	int max_events = std::numeric_limits<int>::max();
	if(argn > 3){
	max_events = std::stoi(argv[3]);
	const int input_events = event_number(reader->header());
	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
	};
	auto 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
	};

	// status infos & eye candy
	int nevent = 0;
	std::array<int, HEJ::event_type::last_type + 1>
	nevent_type{0}, nfailed_type{0};
	auto progress = make_progress_bar(reader->heprup().XSECUP);
	HEJ::CrossSectionAccumulator xs;
	std::map<HEJ::StatusCode, int> status_counter;
	size_t total_trials = 0;

	// Loop over the events in the input file
	while(reader->read_event()){
	// reweight events so that the total cross section is conserved
	auto hepeup = reader->hepeup();
	hepeup.setWeight(0, global_reweight * hepeup.weight());

	if(nevent == max_events) break;
	++nevent;

	HEJ::Event::EventData event_data{hepeup};
	event_data.reconstruct_intermediate();

	// calculate HEJ weight
	HEJ::Event FO_event{
	std::move(event_data).cluster(
	config.fixed_order_jets.def, config.fixed_order_jets.min_pt
	)
	};
	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);
	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()];

	for(auto const & ev: resummed_events){
	//TODO: move pass_cuts to after phase space point generation
	if(analysis->pass_cuts(ev, FO_event)){
	analysis->fill(ev, FO_event);
	writer.write(ev);
	xs.fill(ev);
	}
	}
	if(progress) progress->increment(FO_event.central().weight);
	} // main event loop
	std::cout << '\n';
	analysis->finalise();

	using namespace HEJ::event_type;
	std::cout<< "Events processed: " << nevent << '\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";

	}
	diff --git a/t/check_rivet.cc b/t/check_rivet.cc
	new file mode 100644
	index 0000000..4d54eb4
	--- /dev/null
	+++ b/t/check_rivet.cc
	@@ -0,0 +1,73 @@
	+/**
	+ * \authors The HEJ collaboration (see AUTHORS for details)
	+ * \date 2019
	+ * \copyright GPLv2 or later
	+ */
	+/**
	+ * small test to see if we can rivet interface, and if rivet accepts our beams
	+ */
	+#include <array>
	+#include <iostream>
	+#include <string>
	+#include <vector>
	+
	+#include "HEJ/Event.hh"
	+#include "HEJ/exceptions.hh"
	+#include "HEJ/PDG_codes.hh"
	+#include "HEJ/RivetAnalysis.hh"
	+
	+#include "LHEF/LHEF.h"
	+
	+#include "Rivet/AnalysisHandler.hh"
	+
	+#include "yaml-cpp/yaml.h"
	+
	+namespace {
	+ const fastjet::JetDefinition jet_def{fastjet::JetAlgorithm::antikt_algorithm, 0.4};
	+ const double min_jet_pt{30.};
	+
	+ void test_analysis( std::vector<std::string> analysis,
	+ std::array<double,2> energy, std::array<HEJ::ParticleID,2> beam
	+ ){
	+ HEJ::Event::EventData ev;
	+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -78, -54, 95}, {}});
	+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -64, -44, 81}, {}});
	+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, -24, 113}, {}});
	+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, -8, 93}, {}});
	+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, 44, 109}, {}});
	+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -52, 36, 65}, {}});
	+ ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -84, 97, 137}, {}});
	+ ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -323, 323}, {}};
	+ ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 370, 370}, {}};
	+ auto const event{ ev.cluster(jet_def, min_jet_pt) };
	+
	+ std::string yaml_setup{ "{output: rivet_test, rivet: [" };
	+ for(auto const & ana: analysis)
	+ yaml_setup+=ana+",";
	+ yaml_setup.pop_back();
	+ yaml_setup+="]}";
	+ auto ana_conf{ YAML::Load(yaml_setup) };
	+
	+ LHEF::HEPRUP heprup;
	+ heprup.IDBMUP.first = beam[0];
	+ heprup.IDBMUP.second = beam[1];
	+ heprup.EBMUP.first = energy[0];
	+ heprup.EBMUP.second = energy[1];
	+
	+ HEJ::RivetAnalysis rivet{ ana_conf };
	+ rivet.initialise(heprup);
	+ if(!rivet.pass_cuts(event, event))
	+ throw std::logic_error{"Rivet should not induce additional cuts!"};
	+ rivet.fill(event, event);
	+ rivet.finalise();
	+ }
	+
	+}
	+
	+int main(){
	+ using namespace HEJ;
	+ test_analysis({"MC_XS"},{6500,6500},{ParticleID::e_bar,ParticleID::e}); // MC_XS accepts everything
	+ test_analysis({"MC_XS","ATLAS_2016_I1419652"},{6500,6500},{ParticleID::p,ParticleID::p});
	+ test_analysis({"ATLAS_2014_I1319490"},{3500,3500},{ParticleID::p,ParticleID::p});
	+ return EXIT_SUCCESS;
	+}

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 11:56 AM (1 h, 36 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5111564
Default Alt Text: (47 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions