No OneTemporary
Actions

Size

104 KB

Subscribers

None

View Options

	diff --git a/Changes-API.md b/Changes-API.md
	index 4019fed..86be069 100644
	--- a/Changes-API.md
	+++ b/Changes-API.md
	@@ -1,129 +1,130 @@
	# 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
	* 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.
	* Renamed `no_2_jets` error flag to `not_enough_jets`.
	* 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.


	## 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/FixedOrderGen/configFO.yml b/FixedOrderGen/configFO.yml
	index d2472dc..2f7eca5 100644
	--- a/FixedOrderGen/configFO.yml
	+++ b/FixedOrderGen/configFO.yml
	@@ -1,95 +1,93 @@
	## Number of generated events
	events: 200

	jets:
	min pt: 20 # minimal jet pt, should be slightly below analysis cut
	peak pt: 30 # peak pt of jets, should be at analysis cut
	algorithm: antikt # jet clustering algorithm
	R: 0.4 # jet R parameter
	max rapidity: 5 # maximum jet rapidity

	## Particle beam
	beam:
	energy: 6500
	particles: [p, p]

	## PDF ID according to LHAPDF
	pdf: 230000

	## Scattering process
	process: p p => Wp 4j

	## Particle decays (multiple decays are allowed)
	decays:
	Higgs: {into: [photon, photon], branching ratio: 0.0023568762400521404}
	Wp: {into: [e+, nu_e], branching ratio: 1}
	Wm: {into: [e-, nu_e_bar]} # equivalent to "branching ratio: 1"

	## Fraction of events with two extremal emissions in one direction
	## that contain an subleading emission e.g. unordered emission
	subleading fraction: 0.05

	## Allow different subleading configurations.
	-## By default all processes are allowed.
	-## This does not check if the processes are implemented in HEJ!
	-#
	-# subleading channels:
	+## By default all subleading channels are disabled.
	+subleading channels: none
	# - unordered
	# - central qqbar
	# - extremal qqbar

	## Unweighting parameters
	## remove to obtain weighted events
	# unweight:
	# sample size: 200 # should be similar to "events:", but not more than ~10000
	# max deviation: 0

	## --- The following settings are equivalent to HEJ, --- ##
	## --- see HEJ for reference. --- ##

	## Central scale choice
	scales: max jet pperp

	## Event output files
	event output:
	- HEJFO.lhe
	# - HEJFO_events.hepmc

	## Analyses
	#
	# analyses:
	## Rivet analysis
	# - rivet: MC_XS # rivet analysis name
	# output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added
	## Custom analysis
	# - plugin: /path/to/libmyanalysis.so
	# my analysis parameter: some value

	## Selection of random number generator and seed
	random generator:
	name: mixmax
	# seed: 1

	## Vacuum expectation value
	vev: 246.2196508

	## Properties of the weak gauge bosons
	particle properties:
	Higgs:
	mass: 125
	width: 0.004165
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	## Parameters for Higgs-gluon couplings
	## this requires compilation with QCDloop
	#
	# Higgs coupling:
	# use impact factors: false
	# mt: 174
	# include bottom: true
	# mb: 4.7
	diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc
	index 104dab1..3cdff3b 100644
	--- a/FixedOrderGen/src/config.cc
	+++ b/FixedOrderGen/src/config.cc
	@@ -1,427 +1,488 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "config.hh"

	#include <algorithm>
	#include <cassert>
	#include <cctype>
	#include <cmath>
	#include <cstdlib>
	#include <iostream>
	#include <iterator>
	#include <stdexcept>

	#include "HEJ/PDG_codes.hh"
	#include "HEJ/YAMLreader.hh"
	#include "HEJ/exceptions.hh"
	+#include "HEJ/Event.hh"
	+

	namespace HEJFOG {
	using HEJ::set_from_yaml;
	using HEJ::set_from_yaml_if_defined;
	using HEJ::pid::ParticleID;

	namespace {
	+
	+ //! Convert name of channel in config to name of subleading process in HEJ
	+ size_t channel_to_subleading(
	+ unsigned const channel
	+ ){
	+ using namespace HEJ::event_type;
	+ using namespace HEJFOG::subleading;
	+ switch (channel){
	+ case uno:
	+ return unob;
	+ case cqqbar:
	+ return qqbar_mid;
	+ case eqqbar:
	+ return qqbar_exb;
	+ default:
	+ throw std::logic_error{"Invalid subleading channel"};
	+ }
	+ }
	+
	+ //! Check subleading processes implemented in HEJ2
	+ bool is_valid_subleading(
	+ std::optional<HEJ::ParticleID> const boson,
	+ Subleading const subleading_channels
	+ ){
	+
	+ std::vector<HEJ::Particle> bosons;
	+ if(boson){
	+ HEJ::Particle indicator_boson;
	+ indicator_boson.type = *boson;
	+ bosons.push_back(indicator_boson);
	+ }
	+ size_t const valid_types = HEJ::implemented_types(bosons);
	+
	+ for (
	+ unsigned channel = subleading::first;
	+ channel<=subleading::last;
	+ ++channel
	+ ){
	+ if (subleading_channels[channel]){
	+ if ( !(channel_to_subleading(channel) & valid_types) ){
	+ return false;
	+ }
	+ }
	+ }
	+ return true;
	+ }
	+
	+
	+
	//! Get YAML tree of supported options
	/**
	* The configuration file is checked against this tree of options
	* in assert_all_options_known.
	*/
	YAML::Node const & get_supported_options(){
	const static YAML::Node supported = [](){
	YAML::Node supported;
	static const auto opts = {
	"process", "events", "subleading fraction","subleading channels",
	"scales", "scale factors", "max scale ratio", "pdf", "vev",
	"event output", "analyses", "analysis", "import scales"
	};
	// add subnodes to "supported" - the assigned value is irrelevant
	for(auto && opt: opts) supported[opt] = "";
	for(auto && jet_opt: {"min pt", "peak pt", "algorithm", "R", "max rapidity"}){
	supported["jets"][jet_opt] = "";
	}
	for(auto && particle_type: {"Higgs", "W", "Z"}){
	for(auto && particle_opt: {"mass", "width"}){
	supported["particle properties"][particle_type][particle_opt] = "";
	}
	}
	for(auto && particle_type: {"Higgs", "Wp", "W+", "Wm", "W-", "Z"}){
	for(auto && particle_opt: {"into", "branching ratio"}){
	supported["decays"][particle_type][particle_opt] = "";
	}
	}
	for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
	supported["Higgs coupling"][opt] = "";
	}
	for(auto && beam_opt: {"energy", "particles"}){
	supported["beam"][beam_opt] = "";
	}
	for(auto && unweight_opt: {"sample size", "max deviation"}){
	supported["unweight"][unweight_opt] = "";
	}
	for(auto && opt: {"name", "seed"}){
	supported["random generator"][opt] = "";
	}
	return supported;
	}();
	return supported;
	}

	JetParameters get_jet_parameters(
	YAML::Node const & node, std::string const & entry
	){
	const auto p = HEJ::get_jet_parameters(node, entry);
	JetParameters result;
	result.def = p.def;
	result.min_pt = p.min_pt;
	set_from_yaml(result.max_y, node, entry, "max rapidity");
	set_from_yaml_if_defined(result.peak_pt, node, entry, "peak pt");
	if(result.peak_pt && *result.peak_pt <= result.min_pt)
	throw std::invalid_argument{
	"Value of option 'peak pt' has to be larger than 'min pt'."
	};
	return result;
	}

	Beam get_Beam(
	YAML::Node const & node, std::string const & entry
	){
	Beam beam;
	std::vector<HEJ::ParticleID> particles;
	set_from_yaml(beam.energy, node, entry, "energy");
	set_from_yaml_if_defined(particles, node, entry, "particles");
	if(! particles.empty()){
	for(HEJ::ParticleID particle: particles){
	if(particle != HEJ::pid::p && particle != HEJ::pid::p_bar){
	throw std::invalid_argument{
	"Unsupported value in option " + entry + ": particles:"
	" only proton ('p') and antiproton ('p_bar') beams are supported"
	};
	}
	}
	if(particles.size() != 2){
	throw std::invalid_argument{"Not exactly two beam particles"};
	}
	beam.particles.front() = particles.front();
	beam.particles.back() = particles.back();
	}
	return beam;
	}

	std::vector<std::string> split(
	std::string const & str, std::string const & delims
	){
	std::vector<std::string> result;
	for(size_t begin = 0, end = 0; end != std::string::npos;){
	begin = str.find_first_not_of(delims, end);
	if(begin == std::string::npos) break;
	end = str.find_first_of(delims, begin + 1);
	result.emplace_back(str.substr(begin, end - begin));
	}
	return result;
	}

	std::invalid_argument invalid_incoming(std::string const & what){
	return std::invalid_argument{
	"Incoming particle type " + what + " not supported,"
	" incoming particles have to be 'p', 'p_bar' or partons"
	};
	}

	std::invalid_argument invalid_outgoing(std::string const & what){
	return std::invalid_argument{
	"Outgoing particle type " + what + " not supported,"
	" outgoing particles have to be 'j', 'photon', 'H', 'Wm', 'Wp', 'e-', 'e+', 'nu_e', 'nu_e_bar'"
	};
	}

	HEJ::ParticleID reconstruct_boson_id(
	std::vector<HEJ::ParticleID> const & ids
	){
	assert(ids.size()==2);
	const int pidsum = ids[0] + ids[1];
	switch(pidsum){
	case +1:
	assert(HEJ::is_antilepton(ids[0]));
	if(HEJ::is_antineutrino(ids[0])) {
	throw HEJ::not_implemented{"lepton-flavour violating final state"};
	}
	assert(HEJ::is_neutrino(ids[1]));
	// charged antilepton + neutrino means we had a W+
	return HEJ::pid::Wp;
	case -1:
	assert(HEJ::is_antilepton(ids[0]));
	if(HEJ::is_neutrino(ids[1])) {
	throw HEJ::not_implemented{"lepton-flavour violating final state"};
	}
	assert(HEJ::is_antineutrino(ids[0]));
	// charged lepton + antineutrino means we had a W-
	return HEJ::pid::Wm;
	case 0:
	assert(HEJ::is_anylepton(ids[0]));
	if(HEJ::is_anyneutrino(ids[0])) {
	throw HEJ::not_implemented{"final state with two neutrinos"};
	}
	// charged lepton-antilepton pair means we had a Z/photon
	return HEJ::pid::Z_photon_mix;
	default:
	throw HEJ::not_implemented{
	"final state with leptons "+name(ids[0])+" and "+name(ids[1])
	+" not supported"
	};
	}
	}

	Process get_process(
	YAML::Node const & node, std::string const & entry
	){
	Process result;

	std::string process_string;
	set_from_yaml(process_string, node, entry);
	assert(! process_string.empty());
	const auto particles = split(process_string, " \n\t\v=>");
	if(particles.size() < 3){
	throw std::invalid_argument{
	"Bad format in option process: '" + process_string
	+ "', expected format is 'in1 in2 => out1 ...'"
	};
	}
	result.incoming.front() = HEJ::to_ParticleID(particles[0]);
	result.incoming.back() = HEJ::to_ParticleID(particles[1]);

	for(size_t i = 0; i < result.incoming.size(); ++i){
	const HEJ::ParticleID in = result.incoming[i];
	if(
	in != HEJ::pid::proton && in != HEJ::pid::p_bar
	&& !HEJ::is_parton(in)
	){
	throw invalid_incoming(particles[i]);
	}
	}
	result.njets = 0;
	for(size_t i = result.incoming.size(); i < particles.size(); ++i){
	assert(! particles[i].empty());
	if(particles[i] == "j") ++result.njets;
	else if(std::isdigit(particles[i].front())
	&& particles[i].back() == 'j')
	result.njets += std::stoi(particles[i]);
	else{
	const auto pid = HEJ::to_ParticleID(particles[i]);
	if(pid==HEJ::pid::Higgs \|\| pid==HEJ::pid::Wp \|\| pid==HEJ::pid::Wm){
	if(result.boson)
	throw std::invalid_argument{
	"More than one outgoing boson is not supported"
	};
	if(!result.boson_decay.empty())
	throw std::invalid_argument{
	"Production of a boson together with a lepton is not supported"
	};
	result.boson = pid;
	} else if (HEJ::is_anylepton(pid)){
	// Do not accept more leptons, if two leptons are already mentioned
	if( result.boson_decay.size()>=2 )
	throw std::invalid_argument{"Too many leptons provided"};
	if(result.boson)
	throw std::invalid_argument{
	"Production of a lepton together with a boson is not supported"
	};
	result.boson_decay.emplace_back(pid);
	} else {
	throw invalid_outgoing(particles[i]);
	}
	}
	}
	if(result.njets < 2){
	throw std::invalid_argument{
	"Process has to include at least two jets ('j')"
	};
	}
	if(!result.boson_decay.empty()){
	std::sort(begin(result.boson_decay),end(result.boson_decay));
	assert(!result.boson);
	result.boson = reconstruct_boson_id(result.boson_decay);
	}
	return result;
	}

	HEJFOG::Subleading to_subleading_channel(YAML::Node const & yaml){
	std::string name;
	using namespace HEJFOG::subleading;
	set_from_yaml(name, yaml);
	if(name == "none")
	return NONE;
	if(name == "all")
	return ALL;
	Subleading channel = NONE;
	if(name == "unordered" \|\| name == "uno")
	return channel.set(uno);
	if(name == "central qqbar" \|\| name == "cqqbar")
	return channel.set(cqqbar);
	if(name == "extremal qqbar" \|\| name == "eqqbar")
	return channel.set(eqqbar);
	throw HEJ::unknown_option("Unknown subleading channel '"+name+"'");
	}

	Subleading get_subleading_channels(YAML::Node const & node){
	using YAML::NodeType;
	using namespace HEJFOG::subleading;
	- // all channels allowed by default
	- if(!node) return ALL;
	+ if(!node)
	+ throw std::invalid_argument{"Subleading channels must be specified"};
	switch(node.Type()){
	case NodeType::Undefined:
	- return ALL;
	case NodeType::Null:
	- return NONE;
	+ throw std::invalid_argument("Subleading channels must be specified");
	case NodeType::Scalar:
	return to_subleading_channel(node);
	case NodeType::Map:
	throw HEJ::invalid_type{"map is not a valid option for subleading channels"};
	case NodeType::Sequence:
	Subleading channels;
	for(auto && channel_node: node){
	channels \|= get_subleading_channels(channel_node);
	}
	return channels;
	}
	throw std::logic_error{"unreachable"};
	}

	Decay get_decay(YAML::Node const & node,
	std::string const & entry, std::string const & boson
	){
	Decay decay;
	set_from_yaml(decay.products, node, entry, boson, "into");
	decay.branching_ratio=1;
	set_from_yaml_if_defined(decay.branching_ratio, node, entry, boson,
	"branching ratio");
	return decay;
	}

	std::vector<Decay> get_decays(YAML::Node const & node,
	std::string const & entry, std::string const & boson
	){
	using YAML::NodeType;
	if(!node[entry][boson]) return {};
	switch(node[entry][boson].Type()){
	case NodeType::Null:
	case NodeType::Undefined:
	return {};
	case NodeType::Scalar:
	throw HEJ::invalid_type{"value is not a list of decays"};
	case NodeType::Map:
	return {get_decay(node, entry, boson)};
	case NodeType::Sequence:
	std::vector<Decay> result;
	for(auto && decay_str: node[entry][boson]){
	result.emplace_back(get_decay(decay_str, entry, boson));
	}
	return result;
	}
	throw std::logic_error{"unreachable"};
	}

	ParticlesDecayMap get_all_decays(YAML::Node const & node,
	std::string const & entry
	){
	if(!node[entry]) return {};
	if(!node[entry].IsMap())
	throw HEJ::invalid_type{entry + " have to be a map"};
	ParticlesDecayMap result;
	for(auto const & sub_node: node[entry]) {
	const auto boson = sub_node.first.as<std::string>();
	const auto id = HEJ::to_ParticleID(boson);
	result.emplace(id, get_decays(node, entry, boson));
	}
	return result;
	}

	UnweightSettings get_unweight(
	YAML::Node const & node, std::string const & entry
	){
	UnweightSettings result{};
	set_from_yaml(result.sample_size, node, entry, "sample size");
	if(result.sample_size <= 0){
	throw std::invalid_argument{
	"negative sample size " + std::to_string(result.sample_size)
	};
	}
	set_from_yaml(result.max_dev, node, entry, "max deviation");
	return result;
	}

	Config to_Config(YAML::Node const & yaml){
	try{
	HEJ::assert_all_options_known(yaml, get_supported_options());
	}
	catch(HEJ::unknown_option const & ex){
	throw HEJ::unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
	}

	Config config;
	config.process = get_process(yaml, "process");
	set_from_yaml(config.events, yaml, "events");
	config.jets = get_jet_parameters(yaml, "jets");
	config.beam = get_Beam(yaml, "beam");
	for(size_t i = 0; i < config.process.incoming.size(); ++i){
	auto const & in = config.process.incoming[i];
	using namespace HEJ::pid;
	if( (in == p \|\| in == p_bar) && in != config.beam.particles[i]){
	throw std::invalid_argument{
	"Particle type of beam " + std::to_string(i+1) + " incompatible"
	+ " with type of incoming particle " + std::to_string(i+1)
	};
	}
	}
	set_from_yaml(config.pdf_id, yaml, "pdf");

	set_from_yaml(config.subleading_fraction, yaml, "subleading fraction");
	if(config.subleading_fraction == 0)
	config.subleading_channels.reset();
	else
	config.subleading_channels = get_subleading_channels(yaml["subleading channels"]);

	+ //check to see if subleading channels allowed for given process
	+ if (!is_valid_subleading(
	+ config.process.boson,
	+ config.subleading_channels
	+ )
	+ ){
	+ throw HEJ::not_implemented{
	+ "Invalid subleading channels for given process"
	+ };
	+ }
	+
	config.ew_parameters = HEJ::get_ew_parameters(yaml);
	config.particle_decays = get_all_decays(yaml, "decays");
	if(config.process.boson // check that Ws always decay
	&& std::abs(*config.process.boson) == HEJ::ParticleID::Wp
	&& config.process.boson_decay.empty()
	){
	auto const & decay = config.particle_decays.find(*config.process.boson);
	if(decay == config.particle_decays.end() \|\| decay->second.empty())
	throw std::invalid_argument{
	"Decay for "+name(*config.process.boson)+" is required"};
	}
	set_from_yaml_if_defined(config.analyses_parameters, yaml, "analyses");
	if(yaml["analysis"].IsDefined()){
	std::cerr <<
	"WARNING: Configuration entry 'analysis' is deprecated. "
	" Use 'analyses' instead.\n";
	YAML::Node ana;
	set_from_yaml(ana, yaml, "analysis");
	if(!ana.IsNull()){
	config.analyses_parameters.push_back(ana);
	}
	}
	config.scales = HEJ::to_ScaleConfig(yaml);
	set_from_yaml_if_defined(config.output, yaml, "event output");
	config.rng = HEJ::to_RNGConfig(yaml, "random generator");
	config.Higgs_coupling = HEJ::get_Higgs_coupling(yaml, "Higgs coupling");
	if(yaml["unweight"].IsDefined()) config.unweight = get_unweight(yaml, "unweight");
	return config;
	}

	} // namespace

	Config load_config(std::string const & config_file){
	try{
	return to_Config(YAML::LoadFile(config_file));
	}
	catch(...){
	std::cerr << "Error reading " << config_file << ":\n ";
	throw;
	}
	}
	} // namespace HEJFOG
	diff --git a/FixedOrderGen/t/CMakeLists.txt b/FixedOrderGen/t/CMakeLists.txt
	index bd76158..1ba17e9 100644
	--- a/FixedOrderGen/t/CMakeLists.txt
	+++ b/FixedOrderGen/t/CMakeLists.txt
	@@ -1,241 +1,283 @@
	set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}")
	set(runcard_dir ${tst_dir}/runcards)

	function(add_long_test)
	set(options NONE)
	set(oneValueArgs NAME)
	set(multiValueArgs COMMAND)
	cmake_parse_arguments(TEST "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN})
	add_test(
	NAME ${TEST_NAME}_short
	COMMAND ${TEST_COMMAND} short
	)
	if(TEST_ALL)
	add_test(
	NAME ${TEST_NAME}
	COMMAND ${TEST_COMMAND}
	)
	endif()
	endfunction()

	foreach(tst W_reconstruct_enu W_2j_classify W_nj_classify)
	add_executable(test_${tst} ${tst_dir}/${tst}.cc)
	target_link_libraries(test_${tst} hejfog_lib)
	add_test(NAME ${tst}
	COMMAND test_${tst}
	WORKING_DIRECTORY ${runcard_dir}
	)
	endforeach()
	foreach(tst Z_2j_classify Z_nj_classify)
	add_executable(test_${tst} ${tst_dir}/${tst}.cc)
	target_link_libraries(test_${tst} hejfog_lib)
	add_test(NAME ${tst}
	COMMAND test_${tst}
	WORKING_DIRECTORY ${runcard_dir}
	)
	endforeach()
	# this only tests if the runcard actually works, not if the result is correct
	add_test(
	NAME main_example
	COMMAND HEJFOG ${PROJECT_SOURCE_DIR}/configFO.yml
	)
	add_test(
	NAME main_2j
	COMMAND HEJFOG ${runcard_dir}/2j.yml
	)
	add_test(
	NAME main_h2j
	COMMAND HEJFOG ${runcard_dir}/h2j.yml
	)
	add_test(
	NAME main_h2j_decay
	COMMAND HEJFOG ${runcard_dir}/h2j_decay.yml
	)
	add_test(
	NAME peakpt
	COMMAND HEJFOG ${runcard_dir}/2j_peak.yml
	)
	# special case where coupling not trivial
	add_test(
	NAME main_W4j_uno+eqqbar
	COMMAND HEJFOG ${runcard_dir}/Wp4j_uno+eqqbar.yml
	)

	# check that uno emission doesn't change FKL xs
	add_executable(FKL_uno FKL_uno.cc)
	target_link_libraries(FKL_uno hejfog_lib)
	add_test(
	NAME FKL_uno
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND FKL_uno ${runcard_dir}/h3j_uno.yml 0.0243548 0.000119862
	)
	# xs tests
	add_executable(xs_gen xs_gen.cc)
	target_link_libraries(xs_gen hejfog_lib)
	## Higgs
	add_test(
	NAME xs_hqQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/hqQ.yml 1.612e-02 1.26303e-04
	)
	add_test(
	NAME xs_h2j
	COMMAND xs_gen ${runcard_dir}/h2j.yml 2.718938e 0.01428279
	)
	add_test(
	NAME xs_h3j
	COMMAND xs_gen ${runcard_dir}/h3j.yml 1.249224e 0.01264546
	)
	add_long_test(
	NAME xs_h3j_uno
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND xs_gen ${runcard_dir}/h3j_uno.yml 0.00347538 3.85875e-05
	)
	add_test(
	NAME xs_h5j
	COMMAND xs_gen ${runcard_dir}/h5j.yml 0.01227874 0.0002203041
	)
	add_test(
	NAME xs_h2j_decay
	COMMAND xs_gen ${runcard_dir}/h2j_decay.yml 6.215591e-03 3.285611e-05
	)
	## pure jets
	add_test(
	NAME xs_qQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/qQ.yml 7.354e+05 1.905e+03
	)
	add_test(
	NAME xs_2j
	# calculated with "combined" HEJ svn r3480
	COMMAND xs_gen ${runcard_dir}/2j.yml 86.42031848e06 590570
	)
	add_test(
	NAME xs_3j_uno
	# calculated with HEJ revision 59b8452accfe1462796406fa8e008c8b16b9657b
	COMMAND xs_gen ${runcard_dir}/3j_uno.yml 1.449280e+05 2.318716e+03
	)
	add_test(
	NAME xs_3j_qqbar
	# calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1
	COMMAND xs_gen ${runcard_dir}/3j_qqbar.yml 62040 1005
	)
	add_long_test(
	NAME xs_4j_qqbar
	# calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1
	COMMAND xs_gen ${runcard_dir}/4j_qqbar.yml 28936 550
	)
	add_test(
	NAME xs_4j
	# calculated with HEJ revision 59b8452accfe1462796406fa8e008c8b16b9657b
	COMMAND xs_gen ${runcard_dir}/4j.yml 1.042808e+06 2.137257e+04
	)
	## W
	add_test(
	NAME xs_WqQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/WpqQ.yml 3.086e+00 4.511e-02
	)
	add_test(
	NAME xs_W2j
	# calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	COMMAND xs_gen ${runcard_dir}/Wm2j.yml 4.177443e+02 7.446928e+00
	)
	add_long_test(
	NAME xs_W3j_eqqbar
	# calculated with HEJ revision 667eb768fbefa99148bf6fe9ffb6fcd16c0f976e
	COMMAND xs_gen ${runcard_dir}/Wp3j_qqbar.yml 2.267785e+01 3.707774e-01
	)
	add_long_test(
	NAME xs_W3j_uno
	# calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	COMMAND xs_gen ${runcard_dir}/Wp3j_uno.yml 3.000955e-01 5.831799e-03
	)
	add_long_test(
	NAME xs_W4j_qqbar
	# calculated with HEJ revision e8929bea8dd53ea52f0e8efd9e235e0b10142aee
	COMMAND xs_gen ${runcard_dir}/Wp4j_qqbar.yml 1.083712e-01 3.297341e-03
	)
	## Z
	add_long_test(
	NAME xs_Z2j
	# calculated with HEJ revision c6cf722b9ed20531422f0723f311088e0d19c396
	COMMAND xs_gen ${runcard_dir}/Z2j.yml 1.018429e+02 6.063405e-01
	)
	add_long_test(
	NAME xs_Z3j_uno
	# calculated with HEJ revision c6cf722b9ed20531422f0723f311088e0d19c396
	COMMAND xs_gen ${runcard_dir}/Z3j_uno.yml 2.743604e-02 5.593320e-04
	)

	# Test that sum of partons == proton
	add_executable(PSP_channel PSP_channel.cc)
	target_link_libraries(PSP_channel hejfog_lib)
	# pure jets
	add_long_test(
	NAME channel_2j
	COMMAND PSP_channel ${runcard_dir}/2j.yml
	)
	add_test(
	NAME channel_3j_qqbar
	COMMAND PSP_channel ${runcard_dir}/3j_qqbar.yml
	)
	add_long_test(
	NAME channel_3j_uno
	COMMAND PSP_channel ${runcard_dir}/3j_uno.yml
	)
	add_long_test(
	NAME channel_4j_qqbar
	COMMAND PSP_channel ${runcard_dir}/4j_qqbar.yml
	)
	# W+jets
	# pure jets
	add_long_test(
	NAME channel_W2j
	COMMAND PSP_channel ${runcard_dir}/Wm2j.yml
	)
	add_long_test(
	NAME channel_W3j_uno
	COMMAND PSP_channel ${runcard_dir}/Wp3j_uno.yml
	)
	add_long_test(
	NAME channel_W3j_eqqbar
	COMMAND PSP_channel ${runcard_dir}/Wp3j_qqbar.yml
	)
	add_long_test(
	NAME channel_W4j_qqbar
	COMMAND PSP_channel ${runcard_dir}/Wp4j_qqbar.yml
	)
	# Z+jets
	add_long_test(
	NAME channel_Z2j
	COMMAND PSP_channel ${runcard_dir}/Z2j.yml
	)
	add_long_test(
	NAME channel_Z3j_uno
	COMMAND PSP_channel ${runcard_dir}/Z3j_uno.yml
	)

	# Test that each subleading channel is generated consistently
	add_executable(PSP_subl PSP_subl.cc)
	target_link_libraries(PSP_subl hejfog_lib)
	add_long_test(
	NAME subl_5j
	COMMAND PSP_subl ${runcard_dir}/5j.yml
	)
	add_long_test(
	NAME subl_h5j
	COMMAND PSP_subl ${runcard_dir}/h5j.yml
	)
	add_long_test(
	NAME subl_W3j
	COMMAND PSP_subl ${runcard_dir}/Wp3j_qqbar.yml
	)
	add_long_test(
	NAME subl_W5j
	COMMAND PSP_subl ${runcard_dir}/Wm5j.yml
	)
	add_long_test(
	NAME subl_Z3j
	COMMAND PSP_subl ${runcard_dir}/Z3j_uno.yml
	)
	+
	+# Test that invalid choices for subleading processes cause the correct errors
	+set(BOSONTYPE h)
	+set(SUBLEADINGCHANNEL cqqbar)
	+configure_file( ${runcard_dir}/unsupported.yml.in
	+ check_unsupported_h_cqqbar.yml @ONLY )
	+set(SUBLEADINGCHANNEL eqqbar)
	+configure_file( ${runcard_dir}/unsupported.yml.in
	+ check_unsupported_h_eqqbar.yml @ONLY )
	+set(BOSONTYPE "e+ e-")
	+set(SUBLEADINGCHANNEL cqqbar)
	+configure_file( ${runcard_dir}/unsupported.yml.in
	+ check_unsupported_z_cqqbar.yml @ONLY )
	+set(SUBLEADINGCHANNEL eqqbar)
	+configure_file( ${runcard_dir}/unsupported.yml.in
	+ check_unsupported_z_eqqbar.yml @ONLY )
	+
	+add_test(
	+ NAME check_unsupported_h_cqqbar
	+ COMMAND HEJFOG check_unsupported_h_cqqbar.yml
	+)
	+add_test(
	+ NAME check_unsupported_h_eqqbar
	+ COMMAND HEJFOG check_unsupported_h_eqqbar.yml
	+)
	+add_test(
	+ NAME check_unsupported_z_cqqbar
	+ COMMAND HEJFOG check_unsupported_z_cqqbar.yml
	+)
	+add_test(
	+ NAME check_unsupported_z_eqqbar
	+ COMMAND HEJFOG check_unsupported_z_eqqbar.yml
	+)
	+set_tests_properties(
	+ check_unsupported_h_cqqbar
	+ check_unsupported_h_eqqbar
	+ check_unsupported_z_cqqbar
	+ check_unsupported_z_eqqbar
	+ PROPERTIES
	+ PASS_REGULAR_EXPRESSION ".*Invalid subleading channels for given process"
	+)
	+
	diff --git a/FixedOrderGen/t/runcards/5j.yml b/FixedOrderGen/t/runcards/5j.yml
	index a1af6a9..69d325c 100644
	--- a/FixedOrderGen/t/runcards/5j.yml
	+++ b/FixedOrderGen/t/runcards/5j.yml
	@@ -1,40 +1,44 @@
	events: 500000

	jets:
	min pt: 30
	R: 0.4
	algorithm: antikt
	max rapidity: 6

	beam:
	energy: 7000
	particles: [p, p]

	pdf: 11000

	process: p p => 5j

	subleading fraction: 0.1
	+subleading channels:
	+ - cqqbar
	+ - eqqbar
	+ - uno

	scales: 125

	random generator:
	name: mixmax
	seed: 1

	particle properties:
	Higgs:
	mass: 125
	width: 0
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	vev: 246.2196508

	unweight:
	sample size: 100
	max deviation: -10
	diff --git a/FixedOrderGen/t/runcards/Wm5j.yml b/FixedOrderGen/t/runcards/Wm5j.yml
	index 4975ca6..c22edf3 100644
	--- a/FixedOrderGen/t/runcards/Wm5j.yml
	+++ b/FixedOrderGen/t/runcards/Wm5j.yml
	@@ -1,39 +1,43 @@
	events: 4000000

	jets:
	min pt: 50
	R: 0.4
	algorithm: antikt
	max rapidity: 5

	beam:
	energy: 15000
	particles: [p, p]

	pdf: 11000

	process: p p => W- 5j

	subleading fraction: 1.
	+subleading channels:
	+ - cqqbar
	+ - eqqbar
	+ - uno

	scales: 91

	particle properties:
	Higgs:
	mass: 125
	width: 0
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	decays:
	Wm: {into: [e-, electron_antineutrino], branching ratio: 1}

	random generator:
	name: mixmax
	seed: 1

	vev: 246.2196508
	diff --git a/FixedOrderGen/t/runcards/Wp_2j.yml b/FixedOrderGen/t/runcards/Wp_2j.yml
	index 0a9dd09..1ecf4ca 100644
	--- a/FixedOrderGen/t/runcards/Wp_2j.yml
	+++ b/FixedOrderGen/t/runcards/Wp_2j.yml
	@@ -1,36 +1,37 @@
	events: 200000

	jets:
	min pt: 30
	R: 0.4
	algorithm: antikt
	max rapidity: 5

	beam:
	energy: 6500
	particles: [p, p]

	pdf: 11000

	process: p p => e+ nu_e 2j

	subleading fraction: 0.37
	+subleading channels: none

	scales: 125

	particle properties:
	Higgs:
	mass: 125
	width: 0.004165
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	random generator:
	name: mixmax
	seed: 1

	vev: 246.2196508
	diff --git a/FixedOrderGen/t/runcards/Wp_2j_decay.yml b/FixedOrderGen/t/runcards/Wp_2j_decay.yml
	index cdef493..4f577eb 100644
	--- a/FixedOrderGen/t/runcards/Wp_2j_decay.yml
	+++ b/FixedOrderGen/t/runcards/Wp_2j_decay.yml
	@@ -1,39 +1,40 @@
	events: 200000

	jets:
	min pt: 30
	R: 0.4
	algorithm: antikt
	max rapidity: 5

	beam:
	energy: 6500
	particles: [p, p]

	pdf: 11000

	process: p p => Wp 2j

	subleading fraction: 0.37
	+subleading channels: none

	scales: 125

	particle properties:
	Higgs:
	mass: 125
	width: 0.004165
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	decays:
	Wp: {into: [e+, nu_e], branching ratio: 1}

	random generator:
	name: mixmax
	seed: 1

	vev: 246.2196508
	diff --git a/FixedOrderGen/t/runcards/WpqQ.yml b/FixedOrderGen/t/runcards/WpqQ.yml
	index 6903eca..58bdcfd 100644
	--- a/FixedOrderGen/t/runcards/WpqQ.yml
	+++ b/FixedOrderGen/t/runcards/WpqQ.yml
	@@ -1,36 +1,37 @@
	events: 400000

	jets:
	min pt: 20
	algorithm: antikt
	R: 0.4
	max rapidity: 6

	beam:
	energy: 6500
	particles: [p, p]

	pdf: 11000

	process: u s => e+ nu_e 2j

	subleading fraction: .4
	+subleading channels: none

	scales: 125.

	random generator:
	name: mixmax
	seed: 1

	vev: 246.2196508

	particle properties:
	Higgs:
	mass: 125
	width: 0.004165
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.1876
	width: 2.495
	diff --git a/FixedOrderGen/t/runcards/h3j.yml b/FixedOrderGen/t/runcards/h3j.yml
	index a462045..593bfc3 100644
	--- a/FixedOrderGen/t/runcards/h3j.yml
	+++ b/FixedOrderGen/t/runcards/h3j.yml
	@@ -1,37 +1,37 @@
	events: 200000

	jets:
	min pt: 30
	R: 0.4
	algorithm: antikt
	max rapidity: 5

	beam:
	energy: 6500
	particles: [p, p]

	pdf: 11000

	process: p p => h 3j

	subleading fraction: 0.37
	-subleading channels:
	+subleading channels: none

	scales: 125

	particle properties:
	Higgs:
	mass: 125
	width: 0
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	random generator:
	name: mixmax
	seed: 1

	vev: 246.2196508
	diff --git a/FixedOrderGen/t/runcards/h5j.yml b/FixedOrderGen/t/runcards/h5j.yml
	index be69386..4db5faf 100644
	--- a/FixedOrderGen/t/runcards/h5j.yml
	+++ b/FixedOrderGen/t/runcards/h5j.yml
	@@ -1,37 +1,38 @@
	events: 1000000

	jets:
	min pt: 100
	R: 0.4
	algorithm: antikt
	max rapidity: 5

	beam:
	energy: 10000
	particles: [p, p]

	pdf: 11000

	process: p p => h 5j

	subleading fraction: 0.37
	-subleading channels: none
	+subleading channels:
	+ - unordered

	scales: 125

	particle properties:
	Higgs:
	mass: 125
	width: 0
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	random generator:
	name: mixmax
	seed: 1

	vev: 246.2196508
	diff --git a/FixedOrderGen/t/runcards/hqQ.yml b/FixedOrderGen/t/runcards/hqQ.yml
	index e1b7fb9..502b4c0 100644
	--- a/FixedOrderGen/t/runcards/hqQ.yml
	+++ b/FixedOrderGen/t/runcards/hqQ.yml
	@@ -1,41 +1,42 @@
	events: 400000

	jets:
	min pt: 20
	algorithm: antikt
	R: 0.4
	max rapidity: 6

	beam:
	energy: 6500
	particles: [p, p]

	pdf: 11000

	process: u d => h 2j

	subleading fraction: .4
	+subleading channels: none

	scales: 125.

	random generator:
	name: mixmax
	seed: 1

	decays:
	Higgs:
	branching ratio: 1
	into: [photon, photon]

	vev: 246.2196508

	particle properties:
	Higgs:
	mass: 125
	width: 0.004165
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.1876
	width: 2.495
	diff --git a/FixedOrderGen/t/runcards/qQ.yml b/FixedOrderGen/t/runcards/qQ.yml
	index 27b4b3a..7050b69 100644
	--- a/FixedOrderGen/t/runcards/qQ.yml
	+++ b/FixedOrderGen/t/runcards/qQ.yml
	@@ -1,36 +1,37 @@
	events: 400000

	jets:
	min pt: 20
	algorithm: antikt
	R: 0.4
	max rapidity: 6

	beam:
	energy: 7000
	particles: [p, p]

	pdf: 11000

	process: u d => 2j

	subleading fraction: .4
	+subleading channels: none

	scales: 125.

	random generator:
	name: mixmax
	seed: 1

	vev: 246.2196508

	particle properties:
	Higgs:
	mass: 125
	width: 0.004165
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.1876
	width: 2.495
	diff --git a/FixedOrderGen/t/runcards/h3j.yml b/FixedOrderGen/t/runcards/unsupported.yml.in
	similarity index 85%
	copy from FixedOrderGen/t/runcards/h3j.yml
	copy to FixedOrderGen/t/runcards/unsupported.yml.in
	index a462045..ac6c062 100644
	--- a/FixedOrderGen/t/runcards/h3j.yml
	+++ b/FixedOrderGen/t/runcards/unsupported.yml.in
	@@ -1,37 +1,38 @@
	-events: 200000
	+events: 100

	jets:
	min pt: 30
	R: 0.4
	algorithm: antikt
	max rapidity: 5

	beam:
	energy: 6500
	particles: [p, p]

	pdf: 11000

	-process: p p => h 3j
	+process: p p => @BOSONTYPE@ 3j

	subleading fraction: 0.37
	subleading channels:
	+ - @SUBLEADINGCHANNEL@

	scales: 125

	particle properties:
	Higgs:
	mass: 125
	width: 0
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	random generator:
	name: mixmax
	seed: 1

	vev: 246.2196508
	diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
	index dae110a..e4e6186 100644
	--- a/include/HEJ/Event.hh
	+++ b/include/HEJ/Event.hh
	@@ -1,413 +1,416 @@
	/** \file
	* \brief Declares the Event class and helpers
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <array>
	#include <cstddef>
	#include <iosfwd>
	#include <iterator>
	#include <unordered_map>
	#include <utility>
	#include <vector>

	#include "boost/iterator/filter_iterator.hpp"

	#include "fastjet/ClusterSequence.hh"
	#include "fastjet/PseudoJet.hh"

	#include "HEJ/Constants.hh"
	#include "HEJ/Parameters.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/event_types.hh"

	namespace LHEF {
	class HEPEUP;
	class HEPRUP;
	}

	namespace fastjet {
	class JetDefinition;
	}

	namespace HEJ {
	class EWConstants;
	struct RNG;
	struct UnclusteredEvent;

	+ //! Method for accessing implemented types
	+ size_t implemented_types(std::vector<Particle> const & bosons);
	+
	/** @brief An event with clustered jets
	*
	* This is the main HEJ 2 event class.
	* It contains kinematic information including jet clustering,
	* parameter (e.g. scale) settings and the event weight.
	*/
	class Event {
	public:
	class EventData;

	//! Iterator over partons
	using ConstPartonIterator = boost::filter_iterator<
	bool (*)(Particle const &),
	std::vector<Particle>::const_iterator
	>;
	//! Reverse Iterator over partons
	using ConstReversePartonIterator = std::reverse_iterator<
	ConstPartonIterator>;
	//! No default Constructor
	Event() = delete;
	//! Event Constructor adding jet clustering to an unclustered event
	//! @deprecated UnclusteredEvent got superseded by EventData
	//! and will be removed in HEJ 2.2.0
	[[deprecated("UnclusteredEvent got superseded by EventData")]]
	Event(
	UnclusteredEvent const & ev,
	fastjet::JetDefinition const & jet_def, double min_jet_pt
	);

	//! @name Particle Access
	//! @{

	//! Incoming particles
	std::array<Particle, 2> const & incoming() const{
	return incoming_;
	}
	//! Outgoing particles
	std::vector<Particle> const & outgoing() const{
	return outgoing_;
	}
	//! Iterator to the first outgoing parton
	ConstPartonIterator begin_partons() const;
	//! Iterator to the first outgoing parton
	ConstPartonIterator cbegin_partons() const;

	//! Iterator to the end of the outgoing partons
	ConstPartonIterator end_partons() const;
	//! Iterator to the end of the outgoing partons
	ConstPartonIterator cend_partons() const;

	//! Reverse Iterator to the first outgoing parton
	ConstReversePartonIterator rbegin_partons() const;
	//! Reverse Iterator to the first outgoing parton
	ConstReversePartonIterator crbegin_partons() const;
	//! Reverse Iterator to the first outgoing parton
	ConstReversePartonIterator rend_partons() const;
	//! Reverse Iterator to the first outgoing parton
	ConstReversePartonIterator crend_partons() const;

	//! Particle decays
	/**
	* The key in the returned map corresponds to the index in the
	* vector returned by outgoing()
	*/
	std::unordered_map<std::size_t, std::vector<Particle>> const & decays()
	const {
	return decays_;
	}
	//! The jets formed by the outgoing partons, sorted in rapidity
	std::vector<fastjet::PseudoJet> const & jets() const{
	return jets_;
	}
	//! @}

	//! @name Weight variations
	//! @{

	//! All chosen parameter, i.e. scale choices (const version)
	Parameters<EventParameters> const & parameters() const{
	return parameters_;
	}
	//! All chosen parameter, i.e. scale choices
	Parameters<EventParameters> & parameters(){
	return parameters_;
	}

	//! Central parameter choice (const version)
	EventParameters const & central() const{
	return parameters_.central;
	}
	//! Central parameter choice
	EventParameters & central(){
	return parameters_.central;
	}

	//! Parameter (scale) variations (const version)
	std::vector<EventParameters> const & variations() const{
	return parameters_.variations;
	}
	//! Parameter (scale) variations
	std::vector<EventParameters> & variations(){
	return parameters_.variations;
	}

	//! Parameter (scale) variation (const version)
	/**
	* @param i Index of the requested variation
	*/
	EventParameters const & variations(std::size_t i) const{
	return parameters_.variations.at(i);
	}
	//! Parameter (scale) variation
	/**
	* @param i Index of the requested variation
	*/
	EventParameters & variations(std::size_t i){
	return parameters_.variations.at(i);
	}
	//! @}

	//! Indices of the jets the outgoing partons belong to
	/**
	* @param jets Jets to be tested
	* @returns A vector containing, for each outgoing parton,
	* the index in the vector of jets the considered parton
	* belongs to. If the parton is not inside any of the
	* passed jets, the corresponding index is set to -1.
	*/
	std::vector<int> particle_jet_indices(
	std::vector<fastjet::PseudoJet> const & jets
	) const {
	return cs_.particle_jet_indices(jets);
	}
	//! particle_jet_indices() of the Event jets()
	std::vector<int> particle_jet_indices() const {
	return particle_jet_indices(jets());
	}

	//! Jet definition used for clustering
	fastjet::JetDefinition const & jet_def() const{
	return cs_.jet_def();
	}

	//! Minimum jet transverse momentum
	double min_jet_pt() const{
	return min_jet_pt_;
	}

	//! Event type
	event_type::EventType type() const{
	return type_;
	}

	//! Give colours to each particle
	/**
	* @returns true if new colours are generated, i.e. same as is_resummable()
	* @details Colour ordering is done according to leading colour in the MRK
	* limit, see \cite Andersen:2011zd. This only affects \ref
	* is_resummable() "HEJ" configurations, all other \ref event_type
	* "EventTypes" will be ignored.
	* @note This overwrites all previously set colours.
	*/
	bool generate_colours(RNG & /ran/);

	//! Check that current colours are leading in the high energy limit
	/**
	* @details Checks that the colour configuration can be split up in
	* multiple, rapidity ordered, non-overlapping ladders. Such
	* configurations are leading in the MRK limit, see
	* \cite Andersen:2011zd
	*
	* @note This is _not_ to be confused with \ref is_resummable(), however
	* for all resummable states it is possible to create a leading colour
	* configuration, see generate_colours()
	*/
	bool is_leading_colour() const;

	/**
	* @brief Check if given event could have been produced by HEJ
	* @details A HEJ state has to fulfil:
	* 1. type() has to be \ref is_resummable() "resummable"
	* 2. Soft radiation in the tagging jets contributes at most to
	* `soft_pt_regulator` of the total jet \f$ p_\perp \f$
	*
	* @note This is true for any resummed stated produced by the
	* EventReweighter or any \ref is_resummable() "resummable" Leading
	* Order state.
	*
	* @param soft_pt_regulator Maximum transverse momentum fraction from soft
	* radiation in tagging jets
	* @param min_pt Absolute minimal \f$ p_\perp \f$,
	* \b deprecated use soft_pt_regulator instead
	* @return True if this state could have been produced by HEJ
	*/
	bool valid_hej_state(
	double soft_pt_regulator = DEFAULT_SOFT_PT_REGULATOR,
	double min_pt = 0.
	) const;

	//! Check that the incoming momenta are valid
	/**
	* @details Checks that the incoming parton momenta are on-shell and have
	* vanishing transverse components.
	*
	*/
	bool valid_incoming() const;

	private:
	//! \internal
	//! @brief Construct Event explicitly from input.
	/** This is only intended to be called from EventData.
	*
	* \warning The input is taken _as is_, sorting and classification has to be
	* done externally, i.e. by EventData
	*/
	Event(
	std::array<Particle, 2> && incoming,
	std::vector<Particle> && outgoing,
	std::unordered_map<std::size_t, std::vector<Particle>> && decays,
	Parameters<EventParameters> && parameters,
	fastjet::JetDefinition const & jet_def,
	double min_jet_pt
	);

	//! Iterator over partons (non-const)
	using PartonIterator = boost::filter_iterator<
	bool (*)(Particle const &),
	std::vector<Particle>::iterator
	>;
	//! Reverse Iterator over partons (non-const)
	using ReversePartonIterator = std::reverse_iterator<PartonIterator>;

	//! Iterator to the first outgoing parton (non-const)
	PartonIterator begin_partons();
	//! Iterator to the end of the outgoing partons (non-const)
	PartonIterator end_partons();

	//! Reverse Iterator to the first outgoing parton (non-const)
	ReversePartonIterator rbegin_partons();
	//! Reverse Iterator to the first outgoing parton (non-const)
	ReversePartonIterator rend_partons();

	std::array<Particle, 2> incoming_;
	std::vector<Particle> outgoing_;
	std::unordered_map<std::size_t, std::vector<Particle>> decays_;
	std::vector<fastjet::PseudoJet> jets_;
	Parameters<EventParameters> parameters_;
	fastjet::ClusterSequence cs_;
	double min_jet_pt_;
	event_type::EventType type_;
	}; // end class Event


	//! Detect if a backward incoming gluon turns into a backward outgoing Higgs boson
	inline
	bool is_backward_g_to_h(Event const & ev) {
	return ev.outgoing().front().type == pid::higgs
	&& ev.incoming().front().type == pid::gluon;
	}

	//! Detect if a forward incoming gluon turns into a forward outgoing Higgs boson
	inline
	bool is_forward_g_to_h(Event const & ev) {
	return ev.outgoing().back().type == pid::higgs
	&& ev.incoming().back().type == pid::gluon;
	}

	//! Class to store general Event setup, i.e. Phase space and weights
	class Event::EventData {
	public:
	//! Default Constructor
	EventData() = default;
	//! Constructor from LesHouches event information
	EventData(LHEF::HEPEUP const & hepeup);
	//! Constructor with all values given
	EventData(
	std::array<Particle, 2> incoming,
	std::vector<Particle> outgoing,
	std::unordered_map<std::size_t, std::vector<Particle>> decays,
	Parameters<EventParameters> parameters
	):
	incoming(std::move(incoming)), outgoing(std::move(outgoing)),
	decays(std::move(decays)), parameters(std::move(parameters))
	{}

	//! Generate an Event from the stored EventData.
	/**
	* @details Do jet clustering and classification.
	* Use this to generate an Event.
	*
	* @note Calling this function destroys EventData
	*
	* @param jet_def Jet definition
	* @param min_jet_pt minimal \f$p_T\f$ for each jet
	*
	* @returns Full clustered and classified event.
	*/
	Event cluster(
	fastjet::JetDefinition const & jet_def, double min_jet_pt);

	//! Alias for cluster()
	Event operator()(
	fastjet::JetDefinition const & jet_def, double const min_jet_pt){
	return cluster(jet_def, min_jet_pt);
	}

	//! Sort particles in rapidity
	void sort();

	//! Reconstruct intermediate particles from final-state leptons
	/**
	* Final-state leptons are created from virtual photons, W, or Z bosons.
	* This function tries to reconstruct such intermediate bosons if they
	* are not part of the event record.
	*/
	void reconstruct_intermediate(EWConstants const & /ew_parameters/);

	//! Repair momenta of massless particles
	/**
	* This function changes the momenta of massless particles as follows:
	* - Close-to-zero incoming transverse momenta are set to zero.
	* - Nearly on-shell momenta are made lightlike by rescaling energy
	* and spatial components. This rescaling is applied to both incoming
	* and outgoing particles, including decay products.
	*/
	void repair_momenta(double tolerance);

	//! Incoming particles
	std::array<Particle, 2> incoming;
	//! Outcoing particles
	std::vector<Particle> outgoing;
	//! Particle decays in the format {outgoing index, decay products}
	std::unordered_map<std::size_t, std::vector<Particle>> decays;
	//! Parameters, e.g. scale or inital weight
	Parameters<EventParameters> parameters;
	}; // end class EventData

	//! Print Event
	std::ostream& operator<<(std::ostream & os, Event const & ev);

	//! Square of the partonic centre-of-mass energy \f$\hat{s}\f$
	double shat(Event const & ev);

	//! Tolerance parameter for validity check on incoming momenta
	static constexpr double TOL = 1e-6;

	//! Convert an event to a LHEF::HEPEUP
	LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * /heprup/);

	// put deprecated warning at the end, so don't get the warning inside Event.hh,
	// additionally doxygen can not identify [[deprecated]] correctly
	struct [[deprecated("UnclusteredEvent will be replaced by EventData")]]
	UnclusteredEvent;
	//! An event before jet clustering
	//! @deprecated UnclusteredEvent got superseded by EventData
	//! and will be removed in HEJ 2.2.0
	struct UnclusteredEvent{
	//! Default Constructor
	UnclusteredEvent() = default;
	//! Constructor from LesHouches event information
	UnclusteredEvent(LHEF::HEPEUP const & hepeup);

	std::array<Particle, 2> incoming; /*< Incoming Particles /
	std::vector<Particle> outgoing; /*< Outgoing Particles /
	//! Particle decays in the format {outgoing index, decay products}
	std::unordered_map<std::size_t, std::vector<Particle>> decays;
	//! Central parameter (e.g. scale) choice
	EventParameters central;
	std::vector<EventParameters> variations; /*< For parameter variation /
	};

	} // namespace HEJ
	diff --git a/src/Event.cc b/src/Event.cc
	index 7e935ac..8319cf4 100644
	--- a/src/Event.cc
	+++ b/src/Event.cc
	@@ -1,1394 +1,1396 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "HEJ/Event.hh"

	#include <algorithm>
	#include <cassert>
	#include <cstdlib>
	#include <iomanip>
	#include <iterator>
	#include <memory>
	#include <numeric>
	#include <optional>
	#include <ostream>
	#include <string>
	#include <utility>

	#include "fastjet/ClusterSequence.hh"
	#include "fastjet/JetDefinition.hh"
	#include "fastjet/PseudoJet.hh"

	#include "LHEF/LHEF.h"

	#include "HEJ/Constants.hh"
	#include "HEJ/EWConstants.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/RNG.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/LorentzVector.hh"
	#include "HEJ/utility.hh"

	namespace HEJ {

	+ /**
	+ * returns all EventTypes implemented in HEJ
	+ */
	+ size_t implemented_types(std::vector<Particle> const & bosons){
	+ using namespace event_type;
	+ // no bosons
	+ if(bosons.empty()) return FKL \| unob \| unof \| qqbar_exb \| qqbar_exf \| qqbar_mid;
	+ // 1 boson
	+ if(bosons.size()== 1) {
	+ switch (bosons[0].type) {
	+ case ParticleID::Wp:
	+ case ParticleID::Wm:
	+ return FKL \| unob \| unof \| qqbar_exb \| qqbar_exf \| qqbar_mid;
	+ case ParticleID::Z_photon_mix:
	+ return FKL \| unob \| unof;
	+ case ParticleID::h:
	+ return FKL \| unob \| unof;
	+ default:
	+ return non_resummable;
	+ }
	+ }
	+ // 2 bosons
	+ if(bosons.size() == 2) {
	+ // Resum only samesign W events
	+ if(bosons[0].type == ParticleID::Wp && bosons[1].type == ParticleID::Wp) {
	+ return FKL;
	+ }
	+ else if(bosons[0].type == ParticleID::Wm && bosons[1].type == ParticleID::Wm) {
	+ return FKL;
	+ }
	+ }
	+
	+ return non_resummable;
	+ }
	+
	namespace {
	using std::size_t;

	//! LHE status codes
	namespace lhe_status {
	enum Status: int {
	in = -1,
	decay = 2,
	out = 1,
	};
	}
	using LHE_Status = lhe_status::Status;

	//! true if leptonic W decay
	bool valid_W_decay( int const w_type, // sign of W
	std::vector<Particle> const & decays
	){
	if(decays.size() != 2) // no 1->2 decay
	return false;
	const int pidsum = decays[0].type + decays[1].type;
	if( std::abs(pidsum) != 1 \|\| pidsum != w_type ) // correct charge
	return false;
	// leptonic decay (only check first, second follows from pidsum)
	if( w_type == 1 ) // W+
	return is_antilepton(decays[0]) \|\| is_neutrino(decays[0]);
	// W-
	return is_lepton(decays[0]) \|\| is_antineutrino(decays[0]);
	}

	//! true for Z decay to charged leptons
	bool valid_Z_decay(std::vector<Particle> const & decays){
	if(decays.size() != 2) // no 1->2 decay
	return false;
	const int pidsum = decays[0].type + decays[1].type;
	if( std::abs(pidsum) != 0 ) // correct charge
	return false;
	// leptonic decay (only check first, second follows from pidsum)
	return is_anylepton(decays[0]) && !is_anyneutrino(decays[0]);
	}

	//! true if supported decay
	bool valid_decay(std::vector<Particle> const & decays){
	return valid_W_decay(+1, decays) \|\| // Wp
	valid_W_decay(-1, decays) \|\| // Wm
	valid_Z_decay( decays) // Z/gamma
	;
	}

	/// @name helper functions to determine event type
	//@{

	/**
	* \brief check if final state valid for HEJ
	*
	* check final state has the expected number of valid decays for bosons
	* and all the rest are quarks or gluons
	*/
	bool final_state_ok(Event const & ev){
	size_t invalid_decays = ev.decays().size();

	std::vector<Particle> const & outgoing = ev.outgoing();
	for( size_t i=0; i<outgoing.size(); ++i ){
	auto const & out{ outgoing[i] };
	if(is_AWZH_boson(out.type)){
	auto const decay = ev.decays().find(i);

	// Momentum Conservating Decays
	if(decay != ev.decays().cend()) {
	auto const progeny = decay -> second;

	fastjet::PseudoJet res = std::accumulate(
	progeny.cbegin(), progeny.cend(), fastjet::PseudoJet(),
	[](fastjet::PseudoJet & sum, Particle const & p) -> fastjet::PseudoJet {
	return std::move(sum) + p.p;
	}
	);

	if(!nearby(out.p, res, out.E())){ return false; }
	}

	// W decays (required)
	if(std::abs(out.type) == ParticleID::Wp){
	if( decay != ev.decays().cend() &&
	valid_W_decay(out.type>0?+1:-1, decay->second)
	){
	--invalid_decays;
	}
	else return false;
	}

	// Higgs decays (optional)
	else if(out.type == ParticleID::h){
	if(decay != ev.decays().cend()) --invalid_decays;
	}

	// Z decays (required)
	else if(out.type == ParticleID::Z_photon_mix){
	if( decay != ev.decays().cend() &&
	valid_Z_decay(decay->second)
	){
	--invalid_decays;
	}
	else return false;
	}
	}
	else if(! is_parton(out.type)) return false;
	}
	// any invalid decays?
	return invalid_decays == 0;
	}

	- /**
	- * returns all EventTypes implemented in HEJ
	- */
	- size_t implemented_types(std::vector<Particle> const & bosons){
	- using namespace event_type;
	- // no bosons
	- if(bosons.empty()) return FKL \| unob \| unof \| qqbar_exb \| qqbar_exf \| qqbar_mid;
	- // 1 boson
	- if(bosons.size()== 1) {
	- switch (bosons[0].type) {
	- case ParticleID::Wp:
	- case ParticleID::Wm:
	- return FKL \| unob \| unof \| qqbar_exb \| qqbar_exf \| qqbar_mid;
	- case ParticleID::Z_photon_mix:
	- return FKL \| unob \| unof;
	- case ParticleID::h:
	- return FKL \| unob \| unof;
	- default:
	- return non_resummable;
	- }
	- }
	- // 2 bosons
	- if(bosons.size() == 2) {
	- // Resum only samesign W events
	- if(bosons[0].type == ParticleID::Wp && bosons[1].type == ParticleID::Wp) {
	- return FKL;
	- }
	- else if(bosons[0].type == ParticleID::Wm && bosons[1].type == ParticleID::Wm) {
	- return FKL;
	- }
	- }

	- return non_resummable;
	- }

	/**
	* \brief function which determines if type change is consistent with Wp emission.
	* @param in incoming Particle id
	* @param out outgoing Particle id
	* @param is_qqbar Current both incoming/both outgoing?
	*
	* \see is_Wm_Change
	*/
	bool is_Wp_Change(ParticleID in, ParticleID out, bool is_qqbar){
	using namespace pid;
	if(!is_qqbar && (in==d_bar \|\| in==u \|\| in==s_bar \|\| in==c))
	return out == (in-1);
	if( is_qqbar && (in==d \|\| in==u_bar \|\| in==s \|\| in==c_bar))
	return out == -(in+1);
	return false;
	}

	/**
	* \brief function which determines if type change is consistent with Wm emission.
	* @param in incoming Particle id
	* @param out outgoing Particle id
	* @param is_qqbar Current both incoming/both outgoing?
	*
	* Ensures that change type of quark line is possible by a flavour changing
	* Wm emission. Allows checking of is_qqbar currents also.
	*/
	bool is_Wm_Change(ParticleID in, ParticleID out, bool is_qqbar){
	using namespace pid;
	if(!is_qqbar && (in==d \|\| in==u_bar \|\| in==s \|\| in==c_bar))
	return out == (in+1);
	if( is_qqbar && (in==d_bar \|\| in==u \|\| in==s_bar \|\| in==c))
	return out == -(in-1);
	return false;
	}

	/**
	* \brief checks if particle type remains same from incoming to outgoing
	* @param in incoming Particle
	* @param out outgoing Particle
	* @param is_qqbar Current both incoming/outgoing?
	*/
	bool no_flavour_change(ParticleID in, ParticleID out, bool is_qqbar){
	const int qqbarCurrent = is_qqbar?-1:1;
	if(std::abs(in)<=pid::top \|\| in==pid::gluon)
	return (in==out*qqbarCurrent);
	return false;
	}

	bool has_enough_jets(Event const & event){
	if(event.jets().size() >= 2) return true;
	if(event.jets().empty()) return false;
	// check for h+jet
	const auto the_higgs = std::find_if(
	begin(event.outgoing()), end(event.outgoing()),
	[](const auto & particle) { return particle.type == pid::higgs; }
	);
	return the_higgs != end(event.outgoing());
	}

	bool is_gluon_to_Higgs(const ParticleID in, const ParticleID out) {
	return in == pid::gluon && out == pid::Higgs;
	}

	/**
	* \brief check if we have a valid Impact factor
	* @param in incoming Particle
	* @param out outgoing Particle
	* @param is_qqbar Current both incoming/outgoing?
	* @param W_change returns +1 if Wp, -1 if Wm, else 0
	*/
	bool is_valid_impact_factor(
	ParticleID in, ParticleID out, bool is_qqbar, int & W_change
	){
	if( no_flavour_change(in, out, is_qqbar) \|\| is_gluon_to_Higgs(in, out)) {
	return true;
	}
	if( is_Wp_Change(in, out, is_qqbar) ) {
	W_change+=1;
	return true;
	}
	if( is_Wm_Change(in, out, is_qqbar) ) {
	W_change-=1;
	return true;
	}
	return false;
	}

	bool is_extremal_higgs_off_quark(
	const ParticleID in,
	const ParticleID extremal_out,
	const ParticleID out
	) {
	return in == out && extremal_out == pid::higgs && is_anyquark(in);
	}

	//! Returns all possible classifications from the impact factors
	// the beginning points are changed s.t. after the the classification they
	// point to the beginning of the (potential) FKL chain
	// sets W_change: + if Wp change
	// 0 if no change
	// - if Wm change
	// This function can be used with forward & backwards iterators
	template<class OutIterator>
	size_t possible_impact_factors(
	ParticleID incoming_id, // incoming
	OutIterator & begin_out, OutIterator const & end_out, // outgoing
	int & W_change, std::vector<Particle> const & boson,
	bool const backward // backward?
	){
	using namespace event_type;


	if(begin_out == end_out) return non_resummable;
	// keep track of all states that we don't test
	size_t not_tested = qqbar_mid;
	if(backward)
	not_tested \|= unof \| qqbar_exf;
	else
	not_tested \|= unob \| qqbar_exb;

	// Is this LL current?
	if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){
	++begin_out;
	return not_tested \| FKL;
	}
	// q -> H q and qbar -> H qbar are technically not LL,
	// but we treat them as such anyway
	const auto next = std::next(begin_out);
	if(
	// first ensure that the next particle is not part of the other impact factor
	next != end_out
	&& is_extremal_higgs_off_quark(incoming_id, begin_out->type, next->type)
	) {
	std::advance(begin_out, 2);
	return not_tested \| FKL;
	}

	// or NLL current?
	// -> needs two partons in two different jets
	if( std::distance(begin_out, end_out)>=2
	){
	auto next = std::next(begin_out);
	// Is this unordered emisson?
	if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){
	if( is_valid_impact_factor(
	incoming_id, next->type, false, W_change )
	){
	// veto Higgs inside uno
	assert(next!=end_out);
	if( !boson.empty() && boson.front().type == ParticleID::h
	){
	if( (backward && boson.front().rapidity() < next->rapidity())
	\|\|(!backward && boson.front().rapidity() > next->rapidity()))
	return non_resummable;
	}
	begin_out = std::next(next);
	return not_tested \| (backward?unob:unof);
	}
	}
	// Is this QQbar?
	else if( incoming_id==pid::gluon ){
	if( is_valid_impact_factor(
	begin_out->type, next->type, true, W_change )
	){
	// veto Higgs inside qqbar
	assert(next!=end_out);
	if( !boson.empty() && boson.front().type == ParticleID::h
	){
	if( (backward && boson.front().rapidity() < next->rapidity())
	\|\|(!backward && boson.front().rapidity() > next->rapidity()))
	return non_resummable;
	}
	begin_out = std::next(next);
	return not_tested \| (backward?qqbar_exb:qqbar_exf);
	}
	}
	}
	return non_resummable;
	}

	//! Returns all possible classifications from central emissions
	// the beginning points are changed s.t. after the the classification they
	// point to the end of the emission chain
	// sets W_change: + if Wp change
	// 0 if no change
	// - if Wm change
	template<class OutIterator>
	size_t possible_central(
	OutIterator & begin_out, OutIterator const & end_out,
	int & W_change, std::vector<Particle> const & boson
	){
	using namespace event_type;
	// if we already passed the central chain,
	// then it is not a valid all-order state
	if(std::distance(begin_out, end_out) < 0) return non_resummable;
	// keep track of all states that we don't test
	size_t possible = unob \| unof
	\| qqbar_exb \| qqbar_exf;

	// Find the first quark or antiquark emission
	begin_out = std::find_if(
	begin_out, end_out,
	[](Particle const & p) { return is_anyquark(p); }
	);
	// end of chain -> FKL
	if( begin_out==end_out ){
	return possible \| FKL;
	}

	// is this a qqbar-pair?
	// needs two partons in two separate jets
	auto next = std::next(begin_out);
	if(
	next != end_out
	&& is_valid_impact_factor(begin_out->type, next->type, true, W_change)
	){
	// veto Higgs inside qqbar
	if( !boson.empty() && boson.front().type == ParticleID::h
	&& boson.front().rapidity() > begin_out->rapidity()
	&& boson.front().rapidity() < next->rapidity()
	){
	return non_resummable;
	}
	begin_out = std::next(next);
	// remaining chain should be pure FKL (gluon or higgs)
	if(std::any_of(
	begin_out, end_out,
	[](Particle const & p) { return is_anyquark(p); }
	)) {
	return non_resummable;
	}
	return possible \| qqbar_mid;
	}
	return non_resummable;
	}

	namespace {
	bool is_parton_or_higgs(Particle const & p) {
	return is_parton(p) \|\| p.type == pid::higgs;
	}
	}

	/**
	* \brief Checks for all event types
	* @param ev Event
	* @returns Event Type
	*
	*/
	event_type::EventType classify(Event const & ev){
	using namespace event_type;
	if(! has_enough_jets(ev))
	return not_enough_jets;
	// currently we can't handle multiple boson states in the ME. So they are
	// considered "bad_final_state" even though the "classify" could work with
	// them.
	if(! final_state_ok(ev))
	return bad_final_state;

	// initialise variables
	auto const & in = ev.incoming();

	// range for current checks
	auto begin_out = boost::make_filter_iterator(
	is_parton_or_higgs, cbegin(ev.outgoing()), cend(ev.outgoing())
	);
	auto rbegin_out = std::make_reverse_iterator(
	boost::make_filter_iterator(
	is_parton_or_higgs, cend(ev.outgoing()), cend(ev.outgoing())
	)
	);

	assert(std::distance(begin(in), end(in)) == 2);
	assert(std::distance(begin_out, rbegin_out.base()) >= 2);
	assert(std::is_sorted(begin_out, rbegin_out.base(), rapidity_less{}));

	auto const bosons{ filter_AWZH_bosons(ev.outgoing()) };

	// keep track of potential W couplings, at the end the sum should be 0
	int remaining_Wp = 0;
	int remaining_Wm = 0;
	for(auto const & boson : bosons){
	if(boson.type == ParticleID::Wp) ++remaining_Wp;
	else if(boson.type == ParticleID::Wm) ++remaining_Wm;
	}

	size_t final_type = ~(not_enough_jets \| bad_final_state);

	// check forward impact factor
	int W_change = 0;
	final_type &= possible_impact_factors(
	in.front().type,
	begin_out, rbegin_out.base(),
	W_change, bosons, true );
	if( final_type == non_resummable )
	return non_resummable;
	if(W_change>0) remaining_Wp-=W_change;
	else if(W_change<0) remaining_Wm+=W_change;

	// check backward impact factor
	W_change = 0;
	final_type &= possible_impact_factors(
	in.back().type,
	rbegin_out, std::make_reverse_iterator(begin_out),
	W_change, bosons, false );
	if( final_type == non_resummable )
	return non_resummable;
	if(W_change>0) remaining_Wp-=W_change;
	else if(W_change<0) remaining_Wm+=W_change;

	// check central emissions
	W_change = 0;
	final_type &= possible_central(
	begin_out, rbegin_out.base(), W_change, bosons );
	if( final_type == non_resummable )
	return non_resummable;
	if(W_change>0) remaining_Wp-=W_change;
	else if(W_change<0) remaining_Wm+=W_change;

	// Check whether the right number of Ws are present
	if( remaining_Wp != 0 \|\| remaining_Wm != 0 ) return non_resummable;

	// result has to be unique
	if( (final_type & (final_type-1)) != 0) return non_resummable;

	// check that each sub processes is implemented
	// (has to be done at the end)
	if( (final_type & ~implemented_types(bosons)) != 0 )
	return non_resummable;

	return static_cast<EventType>(final_type);
	}
	//@}

	Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){
	auto id = static_cast<ParticleID>(hepeup.IDUP[i]);
	auto colour = is_parton(id)?hepeup.ICOLUP[i]:std::optional<Colour>();
	return { id,
	{ hepeup.PUP[i][0], hepeup.PUP[i][1],
	hepeup.PUP[i][2], hepeup.PUP[i][3] },
	colour
	};
	}

	bool is_decay_product(std::pair<int, int> const & mothers){
	if(mothers.first == 0) return false;
	return mothers.second == 0 \|\| mothers.first == mothers.second;
	}

	} // namespace

	Event::EventData::EventData(LHEF::HEPEUP const & hepeup){
	parameters.central = EventParameters{
	hepeup.scales.mur, hepeup.scales.muf, hepeup.XWGTUP
	};
	size_t in_idx = 0;
	for (int i = 0; i < hepeup.NUP; ++i) {
	// skip decay products
	// we will add them later on, but we have to ensure that
	// the decayed particle is added before
	if(is_decay_product(hepeup.MOTHUP[i])) continue;

	auto particle = extract_particle(hepeup, i);
	// needed to identify mother particles for decay products
	particle.p.set_user_index(i+1);

	if(hepeup.ISTUP[i] == LHE_Status::in){
	if(in_idx > incoming.size()) {
	throw std::invalid_argument{
	"Event has too many incoming particles"
	};
	}
	incoming[in_idx++] = std::move(particle);
	}
	else outgoing.emplace_back(std::move(particle));
	}

	// add decay products
	for (int i = 0; i < hepeup.NUP; ++i) {
	if(!is_decay_product(hepeup.MOTHUP[i])) continue;
	const int mother_id = hepeup.MOTHUP[i].first;
	const auto mother = std::find_if(
	begin(outgoing), end(outgoing),
	[mother_id](Particle const & particle){
	return particle.p.user_index() == mother_id;
	}
	);
	if(mother == end(outgoing)){
	throw std::invalid_argument{"invalid decay product parent"};
	}
	const int mother_idx = std::distance(begin(outgoing), mother);
	assert(mother_idx >= 0);
	decays[mother_idx].emplace_back(extract_particle(hepeup, i));
	}
	}

	Event::Event(
	UnclusteredEvent const & ev,
	fastjet::JetDefinition const & jet_def, double const min_jet_pt
	):
	Event( Event::EventData{
	ev.incoming, ev.outgoing, ev.decays,
	Parameters<EventParameters>{ev.central, ev.variations}
	}.cluster(jet_def, min_jet_pt) )
	{}

	//! @TODO remove in HEJ 2.2.0
	UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){
	Event::EventData const evData{hepeup};
	incoming = evData.incoming;
	outgoing = evData.outgoing;
	decays = evData.decays;
	central = evData.parameters.central;
	variations = evData.parameters.variations;
	}

	void Event::EventData::sort(){
	// sort particles
	std::sort(
	begin(incoming), end(incoming),
	[](Particle const & o1, Particle const & o2){return o1.p.pz()<o2.p.pz();}
	);

	auto old_outgoing = std::move(outgoing);
	std::vector<size_t> idx(old_outgoing.size());
	std::iota(idx.begin(), idx.end(), 0);
	std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){
	return old_outgoing[i].rapidity() < old_outgoing[j].rapidity();
	});
	outgoing.clear();
	outgoing.reserve(old_outgoing.size());
	for(size_t i: idx) {
	outgoing.emplace_back(std::move(old_outgoing[i]));
	}

	// find decays again
	if(!decays.empty()){
	auto old_decays = std::move(decays);
	decays.clear();
	for(size_t i=0; i<idx.size(); ++i) {
	auto decay = old_decays.find(idx[i]);
	if(decay != old_decays.end())
	decays.emplace(i, std::move(decay->second));
	}
	assert(old_decays.size() == decays.size());
	}
	}

	namespace {
	// use valid_X_decay to determine boson type
	ParticleID reconstruct_type(std::vector<Particle> const & progeny) {
	if(valid_W_decay(+1, progeny)) { return ParticleID::Wp; }
	if(valid_W_decay(-1, progeny)) { return ParticleID::Wm; }
	if(valid_Z_decay(progeny)) { return ParticleID::Z_photon_mix; }

	throw not_implemented{
	"final state with decay X -> "
	+ name(progeny[0].type)
	+ " + "
	+ name(progeny[1].type)
	};
	}

	// reconstruct particle with explicit ParticleID
	Particle reconstruct_boson(
	std::vector<Particle> const & progeny,
	ParticleID const & type
	) {
	Particle progenitor;
	progenitor.p = progeny[0].p + progeny[1].p;
	progenitor.type = type;
	return progenitor;
	}

	// reconstruct via call to reconstruct_type
	Particle reconstruct_boson(std::vector<Particle> const & progeny) {
	Particle progenitor {reconstruct_boson(progeny, reconstruct_type(progeny))};
	assert(is_AWZH_boson(progenitor));
	return progenitor;
	}

	using GroupedParticles = std::vector<std::vector<Particle> >;
	using Decay = std::pair<Particle, std::vector<Particle> >;
	using Decays = std::vector<Decay>;

	// return groups of reconstructable progeny
	std::vector<GroupedParticles> group_progeny(std::vector<Particle> & leptons) {
	/**
	Warning: The partition in to charged/neutral leptons is valid ONLY for WW.
	**/
	assert(leptons.size() == 4);
	auto const begin_neutrino = std::partition(
	begin(leptons), end(leptons),
	[](Particle const & p) {return !is_anyneutrino(p);}
	);

	std::vector<Particle> neutrinos (begin_neutrino, end(leptons));
	leptons.erase(begin_neutrino, end(leptons));

	if(leptons.size() != 2 \|\| neutrinos.size() != 2) { return {}; }
	assert(leptons.size() == 2 && neutrinos.size() == 2);

	std::vector<GroupedParticles> valid_groupings;

	GroupedParticles candidate_grouping{{leptons[0], neutrinos[0]}, {leptons[1], neutrinos[1]}};
	if(valid_decay(candidate_grouping.front()) && valid_decay(candidate_grouping.back())) {
	valid_groupings.emplace_back(std::move(candidate_grouping));
	}

	candidate_grouping = {{leptons[1], neutrinos[0]}, {leptons[0], neutrinos[1]}};
	if(valid_decay(candidate_grouping.front()) && valid_decay(candidate_grouping.back())) {
	valid_groupings.emplace_back(std::move(candidate_grouping));
	}

	return valid_groupings;
	}

	// 'best' decay ordering measure
	double decay_measure(const Particle& reconstructed, EWConstants const & params) {
	ParticleProperties ref = params.prop(reconstructed.type);
	return std::abs(reconstructed.p.m() - ref.mass);
	}

	// decay_measure accumulated over decays
	double decay_measure(Decays const & decays, EWConstants const & params) {
	return
	std::accumulate(
	cbegin(decays), cend(decays), 0.,
	[&params] (double dm, Decay const & decay) -> double {
	return dm + decay_measure(decay.first, params);
	}
	);
	}

	// select best combination of decays for the event
	Decays select_decays
	(
	std::vector<Particle> & leptons,
	EWConstants const & ew_parameters
	) {
	std::vector<GroupedParticles> groupings = group_progeny(leptons);

	std::vector<Decays> valid_decays;
	valid_decays.reserve(groupings.size());

	// Reconstruct all groupings
	for(GroupedParticles const & group : groupings) {
	Decays decays;
	for(auto const & progeny : group) {
	decays.emplace_back(make_pair(reconstruct_boson(progeny), progeny));
	}
	valid_decays.emplace_back(decays);
	}
	if (valid_decays.empty()) {
	throw not_implemented{"No supported intermediate reconstruction available"};
	}
	if (valid_decays.size() == 1) {
	return valid_decays[0];
	}
	// select decay with smallest decay_measure
	auto selected = std::min_element(cbegin(valid_decays), cend(valid_decays),
	[&ew_parameters] (auto const & d1, auto const & d2) -> bool {
	return decay_measure(d1, ew_parameters) < decay_measure(d2, ew_parameters);
	}
	);
	return *selected;
	}
	} // namespace

	void Event::EventData::reconstruct_intermediate(EWConstants const & ew_parameters) {
	auto const begin_leptons = std::partition(
	begin(outgoing), end(outgoing),
	[](Particle const & p) {return !is_anylepton(p);}
	);
	std::vector<Particle> leptons(begin_leptons, end(outgoing));
	outgoing.erase(begin_leptons, end(outgoing));
	if(leptons.empty()) { return; } // nothing to do
	if(leptons.size() == 2) {
	outgoing.emplace_back(reconstruct_boson(leptons));
	std::sort(begin(leptons), end(leptons), type_less{});
	decays.emplace(outgoing.size()-1, std::move(leptons));
	}
	else if(leptons.size() == 4) {
	Decays valid_decays = select_decays(leptons, ew_parameters);
	for(auto &decay : valid_decays) {
	outgoing.emplace_back(decay.first);
	std::sort(begin(decay.second), end(decay.second), type_less{});
	decays.emplace(outgoing.size()-1, std::move(decay.second));
	}
	}
	else {
	throw not_implemented {
	std::to_string(leptons.size())
	+ " leptons in the final state"
	};
	}
	}

	namespace {
	void repair_momentum(fastjet::PseudoJet & p, const double tolerance) {
	if(p.e() > 0. && p.m2() != 0. && (p.m2() < tolerance * tolerance)) {
	const double rescale = std::sqrt(p.modp() / p.e());
	const double e = p.e() * rescale;
	const double px = p.px() / rescale;
	const double py = p.py() / rescale;
	const double pz = p.pz() / rescale;
	p.reset(px, py, pz, e);
	}
	}
	}

	void Event::EventData::repair_momenta(const double tolerance) {
	for(auto & in: incoming) {
	if(is_massless(in)) {
	const double px = (std::abs(in.px()) < tolerance)?0.:in.px();
	const double py = (std::abs(in.py()) < tolerance)?0.:in.py();
	in.p.reset(px, py, in.p.pz(), in.p.e());
	repair_momentum(in.p, tolerance);
	}
	}
	for(auto & out: outgoing) {
	if(is_massless(out)) repair_momentum(out.p, tolerance);
	}
	for(auto & decay: decays) {
	for(auto & out: decay.second) {
	if(is_massless(out)) repair_momentum(out.p, tolerance);
	}
	}
	}

	Event Event::EventData::cluster(
	fastjet::JetDefinition const & jet_def, double const min_jet_pt
	){
	sort();
	return Event{ std::move(incoming), std::move(outgoing), std::move(decays),
	std::move(parameters),
	jet_def, min_jet_pt
	};
	}

	Event::Event(
	std::array<Particle, 2> && incoming,
	std::vector<Particle> && outgoing,
	std::unordered_map<size_t, std::vector<Particle>> && decays,
	Parameters<EventParameters> && parameters,
	fastjet::JetDefinition const & jet_def,
	double const min_jet_pt
	): incoming_{std::move(incoming)},
	outgoing_{std::move(outgoing)},
	decays_{std::move(decays)},
	parameters_{std::move(parameters)},
	cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
	min_jet_pt_{min_jet_pt}
	{
	jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_));
	assert(std::is_sorted(begin(outgoing_), end(outgoing_),
	rapidity_less{}));
	type_ = classify(*this);
	}

	namespace {
	//! check that Particles have a reasonable colour
	bool correct_colour(Particle const & part){
	ParticleID id{ part.type };
	if(!is_parton(id))
	return !part.colour;

	if(!part.colour)
	return false;

	Colour const & col{ *part.colour };
	if(is_quark(id))
	return col.first != 0 && col.second == 0;
	if(is_antiquark(id))
	return col.first == 0 && col.second != 0;
	assert(id==ParticleID::gluon);
	return col.first != 0 && col.second != 0 && col.first != col.second;
	}

	//! Connect parton to t-channel colour line & update the line
	//! returns false if connection not possible
	template<class OutIterator>
	bool try_connect_t(OutIterator const & it_part, Colour & line_colour){
	if( line_colour.first == it_part->colour->second ){
	line_colour.first = it_part->colour->first;
	return true;
	}
	if( line_colour.second == it_part->colour->first ){
	line_colour.second = it_part->colour->second;
	return true;
	}
	return false;
	}

	//! Connect parton to u-channel colour line & update the line
	//! returns false if connection not possible
	template<class OutIterator>
	bool try_connect_u(OutIterator & it_part, Colour & line_colour){
	auto it_next = std::next(it_part);
	if( try_connect_t(it_next, line_colour)
	&& try_connect_t(it_part, line_colour)
	){
	it_part=it_next;
	return true;
	}
	return false;
	}
	} // namespace

	bool Event::is_leading_colour() const {
	if( !correct_colour(incoming()[0]) \|\| !correct_colour(incoming()[1]) )
	return false;

	Colour line_colour = *incoming()[0].colour;
	std::swap(line_colour.first, line_colour.second);

	// reasonable colour
	if(!std::all_of(outgoing().cbegin(), outgoing().cend(), correct_colour))
	return false;

	for(auto it_part = cbegin_partons(); it_part!=cend_partons(); ++it_part){

	switch (type()) {
	case event_type::FKL:
	if( !try_connect_t(it_part, line_colour) )
	return false;
	break;
	case event_type::unob:
	case event_type::qqbar_exb: {
	if( !try_connect_t(it_part, line_colour)
	// u-channel only allowed at impact factor
	&& (std::distance(cbegin_partons(), it_part)!=0
	\|\| !try_connect_u(it_part, line_colour)))
	return false;
	break;
	}
	case event_type::unof:
	case event_type::qqbar_exf: {
	if( !try_connect_t(it_part, line_colour)
	// u-channel only allowed at impact factor
	&& (std::distance(it_part, cend_partons())!=2
	\|\| !try_connect_u(it_part, line_colour)))
	return false;
	break;
	}
	case event_type::qqbar_mid:{
	auto it_next = std::next(it_part);
	if( !try_connect_t(it_part, line_colour)
	// u-channel only allowed at q-qbar/qbar-q pair
	&& ( ( !(is_quark(it_part) && is_antiquark(it_next))
	&& !(is_antiquark(it_part) && is_quark(it_next)))
	\|\| !try_connect_u(it_part, line_colour))
	)
	return false;
	break;
	}
	default:
	throw std::logic_error{"unreachable"};
	}

	// no colour singlet exchange/disconnected diagram
	if(line_colour.first == line_colour.second)
	return false;
	}

	return (incoming()[1].colour->first == line_colour.first)
	&& (incoming()[1].colour->second == line_colour.second);
	}

	namespace {
	//! connect incoming Particle to colour flow
	void connect_incoming(Particle & in, int & colour, int & anti_colour){
	in.colour = std::make_pair(anti_colour, colour);
	// gluon
	if(in.type == pid::gluon)
	return;
	if(in.type > 0){
	// quark
	assert(is_quark(in));
	in.colour->second = 0;
	colour*=-1;
	return;
	}
	// anti-quark
	assert(is_antiquark(in));
	in.colour->first = 0;
	anti_colour*=-1;
	}

	//! connect outgoing Particle to t-channel colour flow
	template<class OutIterator>
	void connect_tchannel(
	OutIterator & it_part, int & colour, int & anti_colour, RNG & ran
	){
	assert(colour>0 \|\| anti_colour>0);
	if(it_part->type == ParticleID::gluon){
	// gluon
	if(colour>0 && anti_colour>0){
	// on g line => connect to colour OR anti-colour (random)
	if(ran.flat() < 0.5){
	it_part->colour = std::make_pair(colour+2,colour);
	colour+=2;
	} else {
	it_part->colour = std::make_pair(anti_colour, anti_colour+2);
	anti_colour+=2;
	}
	} else if(colour > 0){
	// on q line => connect to available colour
	it_part->colour = std::make_pair(colour+2, colour);
	colour+=2;
	} else {
	assert(colour<0 && anti_colour>0);
	// on qbar line => connect to available anti-colour
	it_part->colour = std::make_pair(anti_colour, anti_colour+2);
	anti_colour+=2;
	}
	} else if(is_quark(*it_part)) {
	// quark
	assert(anti_colour>0);
	if(colour>0){
	// on g line => connect and remove anti-colour
	it_part->colour = std::make_pair(anti_colour, 0);
	anti_colour+=2;
	anti_colour*=-1;
	} else {
	// on qbar line => new colour
	colour*=-1;
	it_part->colour = std::make_pair(colour, 0);
	}
	} else if(is_antiquark(*it_part)) {
	// anti-quark
	assert(colour>0);
	if(anti_colour>0){
	// on g line => connect and remove colour
	it_part->colour = std::make_pair(0, colour);
	colour+=2;
	colour*=-1;
	} else {
	// on q line => new anti-colour
	anti_colour*=-1;
	it_part->colour = std::make_pair(0, anti_colour);
	}
	} else { // not a parton
	assert(!is_parton(*it_part));
	it_part->colour = {};
	}
	}

	//! connect to t- or u-channel colour flow
	template<class OutIterator>
	void connect_utchannel(
	OutIterator & it_part, int & colour, int & anti_colour, RNG & ran
	){
	OutIterator it_first = it_part++;
	if(ran.flat()<.5) {// t-channel
	connect_tchannel(it_first, colour, anti_colour, ran);
	connect_tchannel(it_part, colour, anti_colour, ran);
	}
	else { // u-channel
	connect_tchannel(it_part, colour, anti_colour, ran);
	connect_tchannel(it_first, colour, anti_colour, ran);
	}
	}
	} // namespace

	bool Event::generate_colours(RNG & ran){
	// generate only for HEJ events
	if(!event_type::is_resummable(type()))
	return false;
	assert(std::is_sorted(
	begin(outgoing()), end(outgoing()), rapidity_less{}));
	assert(incoming()[0].pz() < incoming()[1].pz());

	// positive (anti-)colour -> can connect
	// negative (anti-)colour -> not available/used up by (anti-)quark
	int colour = COLOUR_OFFSET;
	int anti_colour = colour+1;
	// initialise first
	connect_incoming(incoming_[0], colour, anti_colour);

	// reset outgoing colours
	std::for_each(outgoing_.begin(), outgoing_.end(),
	[](Particle & part){ part.colour = {};});

	for(auto it_part = begin_partons(); it_part!=end_partons(); ++it_part){
	switch (type()) {
	// subleading can connect to t- or u-channel
	case event_type::unob:
	case event_type::qqbar_exb: {
	if( std::distance(begin_partons(), it_part)==0)
	connect_utchannel(it_part, colour, anti_colour, ran);
	else
	connect_tchannel(it_part, colour, anti_colour, ran);
	break;
	}
	case event_type::unof:
	case event_type::qqbar_exf: {
	if( std::distance(it_part, end_partons())==2)
	connect_utchannel(it_part, colour, anti_colour, ran);
	else
	connect_tchannel(it_part, colour, anti_colour, ran);
	break;
	}
	case event_type::qqbar_mid:{
	auto it_next = std::next(it_part);
	if( std::distance(begin_partons(), it_part)>0
	&& std::distance(it_part, end_partons())>2
	&& ( (is_quark(it_part) && is_antiquark(it_next))
	\|\| (is_antiquark(it_part) && is_quark(it_next)) )
	)
	connect_utchannel(it_part, colour, anti_colour, ran);
	else
	connect_tchannel(it_part, colour, anti_colour, ran);
	break;
	}
	default: // rest has to be t-channel
	connect_tchannel(it_part, colour, anti_colour, ran);
	}
	}
	// Connect last
	connect_incoming(incoming_[1], anti_colour, colour);
	assert(is_leading_colour());
	return true;
	} // generate_colours

	namespace {
	bool valid_parton(
	std::vector<fastjet::PseudoJet> const & jets,
	Particle const & parton, int const idx,
	double const soft_pt_regulator, double const min_extparton_pt
	){
	// TODO code overlap with PhaseSpacePoint::pass_extremal_cuts
	if(min_extparton_pt > parton.pt()) return false;
	if(idx<0) return false;
	assert(static_cast<int>(jets.size())>=idx);
	auto const & jet{ jets[idx] };
	return (parton.p - jet).pt()/jet.pt() <= soft_pt_regulator;
	}

	} // namespace

	// this should work with multiple types
	bool Event::valid_hej_state(double const soft_pt_regulator,
	double const min_pt
	) const {
	using namespace event_type;
	if(!is_resummable(type()))
	return false;

	auto const & jet_idx{ particle_jet_indices() };
	auto idx_begin{ jet_idx.cbegin() };
	auto idx_end{ jet_idx.crbegin() };

	auto part_begin{ cbegin_partons() };
	auto part_end{ crbegin_partons() };

	// always seperate extremal jets
	if(!is_backward_g_to_h(*this)) {
	if(! valid_parton(jets(), part_begin, idx_begin, soft_pt_regulator, min_pt)) {
	return false;
	}
	++part_begin;
	++idx_begin;
	// unob -> second parton in own jet
	if( type() & (unob \| qqbar_exb) ){
	if( !valid_parton(jets(), part_begin, idx_begin,
	soft_pt_regulator, min_pt) )
	return false;
	++part_begin;
	++idx_begin;
	}
	}
	if(!is_forward_g_to_h(*this)) {
	if(!valid_parton(jets(), part_end, idx_end, soft_pt_regulator, min_pt)) {
	return false;
	}
	++part_end;
	++idx_end;
	if( type() & (unof \| qqbar_exf) ){
	if( !valid_parton(jets(), part_end, idx_end,
	soft_pt_regulator, min_pt) )
	return false;
	++part_end;
	// ++idx_end; // last check, we don't need idx_end afterwards
	}
	}

	if( type() & qqbar_mid ){
	// find qqbar pair
	auto begin_qqbar{ std::find_if( part_begin, part_end.base(),
	[](Particle const & part) -> bool {
	return part.type != ParticleID::gluon;
	}
	)};
	assert(begin_qqbar != part_end.base());
	long int qqbar_pos{ std::distance(part_begin, begin_qqbar) };
	assert(qqbar_pos >= 0);
	idx_begin+=qqbar_pos;
	if( !( valid_parton(jets(), begin_qqbar, idx_begin,
	soft_pt_regulator, min_pt)
	&& valid_parton(jets(), std::next(begin_qqbar), std::next(idx_begin),
	soft_pt_regulator, min_pt)
	))
	return false;
	}
	return true;
	}

	bool Event::valid_incoming() const{
	for(std::size_t i=0; i < incoming_.size(); ++i){
	if(!(HEJ::nearby_ep(std::abs(incoming_[i].pz()), incoming_[i].E(), TOL*incoming_[i].E())
	&& (incoming_[i].pt()==0.)))
	return false;
	}
	return true;
	}

	Event::ConstPartonIterator Event::begin_partons() const {
	return cbegin_partons();
	}
	Event::ConstPartonIterator Event::cbegin_partons() const {
	return {HEJ::is_parton, cbegin(outgoing()), cend(outgoing())};
	}

	Event::ConstPartonIterator Event::end_partons() const {
	return cend_partons();
	}
	Event::ConstPartonIterator Event::cend_partons() const {
	return {HEJ::is_parton, cend(outgoing()), cend(outgoing())};
	}

	Event::ConstReversePartonIterator Event::rbegin_partons() const {
	return crbegin_partons();
	}
	Event::ConstReversePartonIterator Event::crbegin_partons() const {
	return std::reverse_iterator<ConstPartonIterator>( cend_partons() );
	}

	Event::ConstReversePartonIterator Event::rend_partons() const {
	return crend_partons();
	}
	Event::ConstReversePartonIterator Event::crend_partons() const {
	return std::reverse_iterator<ConstPartonIterator>( cbegin_partons() );
	}

	Event::PartonIterator Event::begin_partons() {
	return {HEJ::is_parton, begin(outgoing_), end(outgoing_)};
	}

	Event::PartonIterator Event::end_partons() {
	return {HEJ::is_parton, end(outgoing_), end(outgoing_)};
	}

	Event::ReversePartonIterator Event::rbegin_partons() {
	return std::reverse_iterator<PartonIterator>( end_partons() );
	}

	Event::ReversePartonIterator Event::rend_partons() {
	return std::reverse_iterator<PartonIterator>( begin_partons() );
	}

	namespace {
	void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){
	constexpr int prec = 6;
	const std::streamsize orig_prec = os.precision();
	os <<std::scientific<<std::setprecision(prec) << "["
	<<std::setw(2*prec+1)<<std::right<< part.px() << ", "
	<<std::setw(2*prec+1)<<std::right<< part.py() << ", "
	<<std::setw(2*prec+1)<<std::right<< part.pz() << ", "
	<<std::setw(2*prec+1)<<std::right<< part.E() << "]"<< std::fixed;
	os.precision(orig_prec);
	}

	void print_colour(std::ostream & os, std::optional<Colour> const & col){
	constexpr int width = 3;
	if(!col)
	os << "(no color)"; // American spelling for better alignment
	else
	os << "(" <<std::setw(width)<<std::right<< col->first
	<< ", " <<std::setw(width)<<std::right<< col->second << ")";
	}
	} // namespace

	std::ostream& operator<<(std::ostream & os, Event const & ev){
	constexpr int prec = 4;
	constexpr int wtype = 3; // width for types
	const std::streamsize orig_prec = os.precision();
	os <<std::setprecision(prec)<<std::fixed;
	os << "########## " << name(ev.type()) << " ##########" << std::endl;
	os << "Incoming particles:\n";
	for(auto const & in: ev.incoming()){
	os <<std::setw(wtype)<< in.type << ": ";
	print_colour(os, in.colour);
	os << " ";
	print_momentum(os, in.p);
	os << std::endl;
	}
	os << "\nOutgoing particles: " << ev.outgoing().size() << "\n";
	for(auto const & out: ev.outgoing()){
	os <<std::setw(wtype)<< out.type << ": ";
	print_colour(os, out.colour);
	os << " ";
	print_momentum(os, out.p);
	os << " => rapidity="
	<<std::setw(2*prec-1)<<std::right<< out.rapidity() << std::endl;
	}
	os << "\nForming Jets: " << ev.jets().size() << "\n";
	for(auto const & jet: ev.jets()){
	print_momentum(os, jet);
	os << " => rapidity="
	<<std::setw(2*prec-1)<<std::right<< jet.rapidity() << std::endl;
	}
	if(!ev.decays().empty() ){
	os << "\nDecays: " << ev.decays().size() << "\n";
	for(auto const & decay: ev.decays()){
	os <<std::setw(wtype)<< ev.outgoing()[decay.first].type
	<< " (" << decay.first << ") to:\n";
	for(auto const & out: decay.second){
	os <<" "<<std::setw(wtype)<< out.type << ": ";
	print_momentum(os, out.p);
	os << " => rapidity="
	<<std::setw(2*prec-1)<<std::right<< out.rapidity() << std::endl;
	}
	}

	}
	os << std::defaultfloat;
	os.precision(orig_prec);
	return os;
	}

	double shat(Event const & ev){
	return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
	}

	LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
	LHEF::HEPEUP result;
	result.heprup = heprup;
	result.weights = {{event.central().weight, nullptr}};
	for(auto const & var: event.variations()){
	result.weights.emplace_back(var.weight, nullptr);
	}
	size_t num_particles = event.incoming().size() + event.outgoing().size();
	for(auto const & decay: event.decays()) num_particles += decay.second.size();
	result.NUP = num_particles;
	// the following entries are pretty much meaningless
	result.IDPRUP = event.type(); // event type
	result.AQEDUP = 1./128.; // alpha_EW
	//result.AQCDUP = 0.118 // alpha_QCD
	// end meaningless part
	result.XWGTUP = event.central().weight;
	result.SCALUP = event.central().muf;
	result.scales.muf = event.central().muf;
	result.scales.mur = event.central().mur;
	result.scales.SCALUP = event.central().muf;
	result.pdfinfo.p1 = event.incoming().front().type;
	result.pdfinfo.p2 = event.incoming().back().type;
	result.pdfinfo.scale = event.central().muf;

	result.IDUP.reserve(num_particles); // PID
	result.ISTUP.reserve(num_particles); // status (in, out, decay)
	result.PUP.reserve(num_particles); // momentum
	result.MOTHUP.reserve(num_particles); // index mother particle
	result.ICOLUP.reserve(num_particles); // colour
	// incoming
	std::array<Particle, 2> incoming{ event.incoming() };
	// First incoming should be positive pz according to LHE standard
	// (or at least most (everyone?) do it this way, and Pythia assumes it)
	if(incoming[0].pz() < incoming[1].pz())
	std::swap(incoming[0], incoming[1]);
	for(Particle const & in: incoming){
	result.IDUP.emplace_back(in.type);
	result.ISTUP.emplace_back(LHE_Status::in);
	result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
	result.MOTHUP.emplace_back(0, 0);
	assert(in.colour);
	result.ICOLUP.emplace_back(*in.colour);
	}
	// outgoing
	for(size_t i = 0; i < event.outgoing().size(); ++i){
	Particle const & out = event.outgoing()[i];
	result.IDUP.emplace_back(out.type);
	const int status = event.decays().count(i) != 0u
	?LHE_Status::decay
	:LHE_Status::out;
	result.ISTUP.emplace_back(status);
	result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
	result.MOTHUP.emplace_back(1, 2);
	if(out.colour)
	result.ICOLUP.emplace_back(*out.colour);
	else{
	result.ICOLUP.emplace_back(std::make_pair(0,0));
	}
	}
	// decays
	for(auto const & decay: event.decays()){
	for(auto const & out: decay.second){
	result.IDUP.emplace_back(out.type);
	result.ISTUP.emplace_back(LHE_Status::out);
	result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
	const size_t mother_idx = 1 + event.incoming().size() + decay.first;
	result.MOTHUP.emplace_back(mother_idx, mother_idx);
	result.ICOLUP.emplace_back(0,0);
	}
	}

	assert(result.ICOLUP.size() == num_particles);
	static constexpr double unknown_spin = 9.; //per Les Houches accord
	result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
	result.SPINUP = result.VTIMUP;
	return result;
	}

	} // namespace HEJ

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:42 PM (1 d, 7 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3798176
Default Alt Text: (104 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions