No OneTemporary
Actions

Size

56 KB

Subscribers

None

View Options

	diff --git a/Changes.md b/Changes.md
	index ed7de34..758d3e1 100644
	--- a/Changes.md
	+++ b/Changes.md
	@@ -1,60 +1,61 @@
	# Changelog

	This is the log for changes to the HEJ program. Further changes to the HEJ API
	are documented in `Changes-API.md`. If you are using HEJ as a library, please
	also read the changes there.

	## Version 2.X

	### 2.X.0

	* Resummation for W bosons with jets
	- New subleading processes `extremal qqx` & `central qqx` for a quark and
	anti-quark in the final state, e.g. `g g => u d_bar Wm g` (the other
	subleading processes also work with W's)
	- `HEJFOG` can generate mutliple jets together with a (off-shell) W bosons
	decaying into lepton & neutrino
	* Resummation can now be performed on `unordered` subleading processes
	in pure jets.
	* Allow multiplication and division of multiple scale functions e.g.
	`H_T/2*m_j1j2`
	* Print cross sections at end of run
	* Follow HepMC convention for particle Status codes: incoming = 11,
	decaying = 2, outgoing = 1 (unchanged)
	* Partons now have a Colour charge
	- Colours are read from and written to LHE files
	- For reweighted events the colours are created according to leading colour in
	the FKL limit
	+* Grouped `event treatment` for subleading channels together in runcard
	+ - Rename `non-HEJ` processes to `non-resummable`
	* Allow changing the regulator lambda in input (`regulator parameter`, only for
	advanced users)
	* Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`)
	* Added support to read `hdf5` event files suggested in
	[arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs
	[HighFive](https://github.com/BlueBrain/HighFive))
	* Support input with avarage weight equal to the cross section (`IDWTUP=1 or 4`)
	-* Rename `non-HEJ` Processes to `non-resummable`
	* Dropped support for HepMC 3.0.0, either HepMC version 2 or >3.1 is required
	- It is now possible to write out both HepMC 2 and HepMC 3 events at the same
	time

	## 2.0.5

	* Fixed event classification for input not ordered in rapidity

	### 2.0.4

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

	### 2.0.3

	* Fixed parsing of (numerical factor) * (base scale) in configuration
	* Don't change scale names, but sanitise Rivet output file names instead

	### 2.0.2

	* Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
	`"/"` and `"*"` before)

	### 2.0.1

	* Fixed name of fixed-order generator in error message.
	diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc
	index 48aeaff..f0d6197 100644
	--- a/FixedOrderGen/src/config.cc
	+++ b/FixedOrderGen/src/config.cc
	@@ -1,409 +1,408 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "config.hh"

	#include <cctype>

	#include "Subleading.hh"

	#include "HEJ/config.hh"
	#include "HEJ/YAMLreader.hh"

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

	namespace{
	//! 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",
	"event output", "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", "Wp", "W+", "Wm", "W-", "Z"}){
	for(auto && particle_opt: {"mass", "width"}){
	supported["particle properties"][particle_type][particle_opt] = "";
	}
	supported["particle properties"][particle_type]["decays"]["into"] = "";
	supported["particle properties"][particle_type]["decays"]["branching ratio"] = "";
	}
	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, end = 0; end != str.npos;){
	begin = str.find_first_not_of(delims, end);
	if(begin == str.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];
	if(pidsum == +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;
	}
	if(pidsum == -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;
	}
	throw HEJ::not_implemented{
	"final state with leptons "+HEJ::name(ids[0])+" and "+HEJ::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::channels;
	set_from_yaml(name, yaml);
	if(name == "none")
	return none;
	if(name == "all")
	return all;
	if(name == "unordered" \|\| name == "uno")
	return uno;
	if(name == "qqx")
	return qqx;
	throw HEJ::unknown_option("Unknown subleading channel '"+name+"'");
	-
	}

	unsigned int get_subleading_channels(YAML::Node const & node){
	using YAML::NodeType;
	using namespace HEJFOG::channels;
	// all channels allowed by default
	if(!node) return all;
	switch(node.Type()){
	case NodeType::Undefined:
	return all;
	case NodeType::Null:
	return none;
	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:
	unsigned int channels = HEJFOG::Subleading::none;
	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){
	Decay decay;
	set_from_yaml(decay.products, node, "into");
	decay.branching_ratio=1;
	set_from_yaml_if_defined(decay.branching_ratio, node, "branching ratio");
	return decay;
	}

	std::vector<Decay> get_decays(YAML::Node const & node){
	using YAML::NodeType;
	if(!node) return {};
	switch(node.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)};
	case NodeType::Sequence:
	std::vector<Decay> result;
	for(auto && decay_str: node){
	result.emplace_back(get_decay(decay_str));
	}
	return result;
	}
	throw std::logic_error{"unreachable"};
	}

	ParticleProperties get_particle_properties(
	YAML::Node const & node, std::string const & entry
	){
	ParticleProperties result;
	set_from_yaml(result.mass, node, entry, "mass");
	set_from_yaml(result.width, node, entry, "width");
	try{
	result.decays = get_decays(node[entry]["decays"]);
	}
	catch(HEJ::missing_option const & ex){
	throw HEJ::missing_option{entry + ": decays: " + ex.what()};
	}
	catch(HEJ::invalid_type const & ex){
	throw HEJ::invalid_type{entry + ": decays: " + ex.what()};
	}
	return result;
	}

	ParticlesPropMap get_all_particles_properties(YAML::Node const & node){
	ParticlesPropMap result;
	for(auto const & entry: node) {
	const auto name = entry.first.as<std::string>();
	const auto id = HEJ::to_ParticleID(name);
	result.emplace(id, get_particle_properties(node,name));
	}
	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){
	const auto & 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_fraction > 1){
	throw std::invalid_argument{
	"subleading fraction has to be between 0 and 1"
	};
	}
	if(config.subleading_fraction == 0)
	config.subleading_channels = Subleading::none;
	else
	config.subleading_channels = get_subleading_channels(yaml["subleading channels"]);

	if(yaml["particle properties"]){
	config.particles_properties = get_all_particles_properties(
	yaml["particle properties"]);
	}
	if(config.process.boson
	&& config.particles_properties.find(*(config.process.boson))
	== config.particles_properties.end())
	throw HEJ::missing_option("Process wants to generate boson "
	+std::to_string(*(config.process.boson))+", but particle properties are missing");
	set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
	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"]) config.unweight = get_unweight(yaml, "unweight");
	return config;
	}

	} // namespace anonymous

	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;
	}
	}
	}
	diff --git a/config.yml b/config.yml
	index b4fa674..a0fabd4 100644
	--- a/config.yml
	+++ b/config.yml
	@@ -1,99 +1,100 @@
	# number of attempted resummation phase space points for each input event
	trials: 10

	min extparton pt: 30 # minimum transverse momentum of extremal partons

	# maximum soft transverse momentum fraction in extremal jets
	#
	# max ext soft pt fraction: 0.1

	resummation jets: # resummation jet properties
	min pt: 35 # minimum jet transverse momentum
	algorithm: antikt # jet clustering algorithm
	R: 0.4 # jet R parameter

	fixed order jets: # properties of input jets
	min pt: 30
	# by default, algorithm and R are like for resummation jets

	# treatment of he various event classes
	# the supported settings are: reweight, keep, discard
	# non-resummable events cannot be reweighted
	-FKL: reweight
	-unordered: keep
	-extremal qqx: keep
	-central qqx: keep
	-non-resummable: keep
	+event treatment:
	+ FKL: reweight
	+ unordered: keep
	+ extremal qqx: keep
	+ central qqx: keep
	+ non-resummable: keep

	# central scale choice or choices
	#
	# scales: [125, max jet pperp, H_T/2, 2*jet invariant mass, m_j1j2]
	scales: 91.188

	# factors by which the central scales should be multiplied
	# renormalisation and factorisation scales are varied independently
	#
	# scale factors: [0.5, 0.7071, 1, 1.41421, 2]

	# maximum ratio between renormalisation and factorisation scale
	#
	# max scale ratio: 2.0001

	# import scale setting functions
	#
	# import scales:
	# lib_my_scales.so: [scale0,scale1]

	log correction: false # whether or not to include higher order logs

	# event output files
	#
	# the supported formats are
	# - Les Houches (suffix .lhe)
	# - HepMC2 (suffix .hepmc)
	# - HepMC3 (suffix .hepmc3 or .hepmc)
	# TODO: - ROOT ntuples (suffix .root)
	#
	# An output file's format is deduced either automatically from the suffix
	# or from an explicit specification, e.g.
	# - Les Houches: outfile

	event output:
	- HEJ.lhe
	# - HEJ_events.hepmc

	# to use a rivet analysis
	#
	# analysis:
	# rivet: MC_XS # rivet analysis name
	# output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added
	#
	# to use a custom analysis
	#
	# analysis:
	# plugin: /path/to/libmyanalysis.so
	# my analysis parameter: some value

	# selection of random number generator and seed
	# the choices are
	# - mixmax (seed is an integer)
	# - ranlux64 (seed is a filename containing parameters)
	random generator:
	name: mixmax
	# seed: 1

	# parameters for Higgs-gluon couplings
	# this requires compilation with qcdloop
	#
	# Higgs coupling:
	# use impact factors: false
	# mt: 174
	# include bottom: true
	# mb: 4.7

	## ---------------------------------------------------------------------- ##
	## The following settings are only intended for advances users. ##
	## Please DO NOT SET them unless you know exactly what you are doing! ##
	## ---------------------------------------------------------------------- ##
	#
	# regulator parameter: 0.2 # The regulator lambda for the subtraction terms
	diff --git a/doc/sphinx/HEJ.rst b/doc/sphinx/HEJ.rst
	index de02a72..cc6171d 100644
	--- a/doc/sphinx/HEJ.rst
	+++ b/doc/sphinx/HEJ.rst
	@@ -1,325 +1,326 @@
	.. _`Running HEJ 2`:

	Running HEJ 2
	=============

	Quick start
	-----------

	In order to run HEJ 2, you need a configuration file and a file
	containing fixed-order events. A sample configuration is given by the
	:file:`config.yml` file distributed together with HEJ 2. Events in the
	Les Houches Event File format can be generated with standard Monte Carlo
	generators like `MadGraph5_aMC@NLO <https://launchpad.net/mg5amcnlo>`_
	or `Sherpa <https://sherpa.hepforge.org/trac/wiki>`_. If HEJ 2 was
	compiled with `HDF5 <https://www.hdfgroup.org/>`_ support, it can also
	read event files in the format suggested in
	`arXiv:1905.05120 <https://arxiv.org/abs/1905.05120>`_.
	HEJ 2 assumes that the cross section is given by the sum of the event
	weights. Depending on the fixed-order generator it may be necessary to
	adjust the weights in the Les Houches Event File accordingly.

	The processes supported by HEJ 2 are

	- Pure multijet production
	- Production of a Higgs boson with jets
	- Production of a W boson with jets

	..
	- TODO Production of a Z boson or photon with jets

	where at least two jets are required in each case. For the time being,
	only leading-order events are supported.

	After generating an event file :file:`events.lhe` adjust the parameters
	under the `fixed order jets`_ setting in :file:`config.yml` to the
	settings in the fixed-order generation. Resummation can then be added by
	running::

	HEJ config.yml events.lhe

	Using the default settings, this will produce an output event file
	:file:`HEJ.lhe` with events including high-energy resummation.

	When using the `Docker image <https://hub.docker.com/r/hejdock/hej>`_,
	HEJ can be run with

	.. code-block:: bash

	docker run -v $PWD:$PWD -w $PWD hejdock/hej HEJ config.yml events.lhe

	.. _`HEJ 2 settings`:

	Settings
	--------

	HEJ 2 configuration files follow the `YAML <http://yaml.org/>`_
	format. The following configuration parameters are supported:

	.. _`trials`:

	trials
	High-energy resummation is performed by generating a number of
	resummation phase space configurations corresponding to an input
	fixed-order event. This parameter specifies how many such
	configurations HEJ 2 should try to generate for each input
	event. Typical values vary between 10 and 100.

	.. _`min extparton pt`:

	min extparton pt
	Specifies the minimum transverse momentum in GeV of the most forward
	and the most backward parton. This setting is needed to regulate an
	otherwise uncancelled divergence. Its value should be slightly below
	the minimum transverse momentum of jets specified by `resummation
	jets: min pt`_. See also the `max ext soft pt fraction`_ setting.

	.. _`max ext soft pt fraction`:

	max ext soft pt fraction
	Specifies the maximum fraction that soft radiation can contribute to
	the transverse momentum of each the most forward and the most backward
	jet. Values between around 0.05 and 0.1 are recommended. See also the
	`min extparton pt`_ setting.

	.. _`fixed order jets`:

	fixed order jets
	This tag collects a number of settings specifying the jet definition
	in the event input. The settings should correspond to the ones used in
	the fixed-order Monte Carlo that generated the input events.

	.. _`fixed order jets: min pt`:

	min pt
	Minimum transverse momentum in GeV of fixed-order jets.

	.. _`fixed order jets: algorithm`:

	algorithm
	The algorithm used to define jets. Allowed settings are
	:code:`kt`, :code:`cambridge`, :code:`antikt`,
	:code:`cambridge for passive`. See the `FastJet
	<http://fastjet.fr/>`_ documentation for a description of these
	algorithms.

	.. _`fixed order jets: R`:

	R
	The R parameter used in the jet algorithm, roughly corresponding
	to the jet radius in the plane spanned by the rapidity and the
	azimuthal angle.

	.. _`resummation jets`:

	resummation jets
	This tag collects a number of settings specifying the jet definition
	in the observed, i.e. resummed events. These settings are optional, by
	default the same values as for the `fixed order jets`_ are assumed.

	.. _`resummation jets: min pt`:

	min pt
	Minimum transverse momentum in GeV of resummation jets. This
	should be between 25% and 50% larger than the minimum transverse
	momentum of fixed order jets set by `fixed order jets: min pt`_.

	.. _`resummation jets: algorithm`:

	algorithm
	The algorithm used to define jets. The HEJ 2 approach to
	resummation relies on properties of :code:`antikt` jets, so this
	value is strongly recommended. For a list of possible other
	values, see the `fixed order jets: algorithm`_ setting.

	.. _`resummation jets: R`:

	R
	The R parameter used in the jet algorithm.

	-.. _`FKL`:
	-
	-FKL
	- Specifies how to treat events respecting FKL rapidity ordering. These
	- configurations are dominant in the high-energy limit. The possible
	- values are :code:`reweight` to enable resummation, :code:`keep` to
	- keep the events as they are up to a possible change of
	- renormalisation and factorisation scale, and :code:`discard` to
	- discard these events.
	-
	-.. _`unordered`:
	-
	-unordered
	-
	- Specifies how to treat events with one emission that does not respect FKL
	- ordering, e.g. :code:`u d => g u d`. In the high-energy limit, such
	- configurations are logarithmically suppressed compared to FKL configurations.
	- The possible values are the same as for the `FKL`_ setting.
	-
	-.. _`extremal qqx`:
	-
	-extremal qqx
	- Specifies how to treat events with a quark-antiquark pair as extremal partons
	- in rapidity, e.g. :code:`g d => u u_bar d`. In the high-energy limit, such
	- configurations are logarithmically suppressed compared to FKL configurations.
	- The possible values are the same as for the `FKL`_ setting, but
	- :code:`reweight` is currently only supported for W boson plus jets
	- production.
	-
	-.. _`central qqx`:
	-
	-central qqx
	- Specifies how to treat events with a quark-antiquark pair central in
	- rapidity, e.g. :code:`g g => g u u_bar g`. In the high-energy limit, such
	- configurations are logarithmically suppressed compared to FKL configurations.
	- The possible values are the same as for the `FKL`_ setting, but
	- :code:`reweight` is currently only supported for W boson plus jets
	- production.
	-
	-.. _`non-resummable`:
	-
	-non-resummable
	- Specifies how to treat events where no resummation is possible. The
	- allowed values are :code:`keep` to keep the events as they are up to
	- a possible change of renormalisation and factorisation scale and
	- :code:`discard` to discard these events.
	-
	-.. _`scales`:
	+.. _`event treatment`:
	+
	+event treatment
	+ Specify how to treat different event types. The different event types
	+ contribute to different orders in the high-energy limit. The possible values
	+ are :code:`reweight` to enable resummation, :code:`keep` to keep the events as
	+ they are up to a possible change of renormalisation and factorisation scale,
	+ and :code:`discard` to discard these events.
	+
	+ .. _`FKL`:
	+
	+ FKL
	+ Specifies how to treat events respecting FKL rapidity ordering. These
	+ configurations are dominant in the high-energy limit.
	+
	+ .. _`unordered`:
	+
	+ unordered
	+ Specifies how to treat events with one emission that does not respect FKL
	+ ordering, e.g. :code:`u d => g u d`. In the high-energy limit, such
	+ configurations are logarithmically suppressed compared to FKL
	+ configurations.
	+
	+ .. _`extremal qqx`:
	+
	+ extremal qqx
	+ Specifies how to treat events with a quark-antiquark pair as extremal
	+ partons in rapidity, e.g. :code:`g d => u u_bar d`. In the high-energy
	+ limit, such configurations are logarithmically suppressed compared to FKL
	+ configurations. :code:`reweight` is currently only supported for W boson
	+ plus jets production.
	+
	+ .. _`central qqx`:
	+
	+ central qqx
	+ Specifies how to treat events with a quark-antiquark pair central in
	+ rapidity, e.g. :code:`g g => g u u_bar g`. In the high-energy limit, such
	+ configurations are logarithmically suppressed compared to FKL
	+ configurations. :code:`reweight` is currently only supported for W boson
	+ plus jets production.
	+
	+ .. _`non-resummable`:
	+
	+ non-resummable
	+ Specifies how to treat events where no resummation is possible. Only
	+ :code:`keep` or :code:`discard` are valid options, not :code:`reweight`
	+ for obvious reasons.
	+
	+ .. _`scales`:

	scales
	Specifies the renormalisation and factorisation scales for the output
	events. This can either be a single entry or a list :code:`[scale1,
	scale2, ...]`. For the case of a list the first entry defines the
	central scale. Possible values are fixed numbers to set the scale in
	GeV or the following:

	- :code:`H_T`: The sum of the scalar transverse momenta of all
	final-state particles
	- :code:`max jet pperp`: The maximum transverse momentum of all jets
	- :code:`jet invariant mass`: Sum of the invariant masses of all jets
	- :code:`m_j1j2`: Invariant mass between the two hardest jets.

	Scales can be multiplied or divided by overall factors, e.g. :code:`H_T/2`.

	It is also possible to import scales from an external library, see
	:ref:`Custom scales`

	.. _`scale factors`:

	scale factors
	A list of numeric factors by which each of the `scales`_ should be
	multiplied. Renormalisation and factorisation scales are varied
	independently. For example, a list with entries :code:`[0.5, 2]`
	would give the four scale choices (0.5μ\ :sub:`r`, 0.5μ\ :sub:`f`);
	(0.5μ\ :sub:`r`, 2μ\ :sub:`f`); (2μ\ :sub:`r`, 0.5μ\ :sub:`f`); (2μ\
	:sub:`r`, 2μ\ :sub:`f`) in this order. The ordering corresponds to
	the order of the final event weights.

	.. _`max scale ratio`:

	max scale ratio
	Specifies the maximum factor by which renormalisation and
	factorisation scales may difer. For a value of :code:`2` and the
	example given for the `scale factors`_ the scale choices
	(0.5μ\ :sub:`r`, 2μ\ :sub:`f`) and (2μ\ :sub:`r`, 0.5μ\ :sub:`f`)
	will be discarded.

	.. _`log correction`:

	log correction
	Whether to include corrections due to the evolution of the strong
	coupling constant in the virtual corrections. Allowed values are
	:code:`true` and :code:`false`.

	.. _`event output`:

	event output
	Specifies the name of a single event output file or a list of such
	files. The file format is either specified explicitly or derived from
	the suffix. For example, :code:`events.lhe` or, equivalently
	:code:`Les Houches: events.lhe` generates an output event file
	:code:`events.lhe` in the Les Houches format. The supported formats
	are

	- :code:`file.lhe` or :code:`Les Houches: file`: The Les Houches
	event file format.
	- :code:`file.hepmc2` or :code:`HepMC2: file`: HepMC format version 2.
	- :code:`file.hepmc3` or :code:`HepMC3: file`: HepMC format version 3.
	- :code:`file.hepmc` or :code:`HepMC: file`: The latest supported
	version of the HepMC format, currently version 3.

	.. _`random generator`:

	random generator
	Sets parameters for random number generation.

	.. _`random generator: name`:

	name
	Which random number generator to use. Currently, :code:`mixmax`
	and :code:`ranlux64` are supported. Mixmax is recommended. See
	the `CLHEP documentation
	<http://proj-clhep.web.cern.ch/proj-clhep/index.html#docu>`_ for
	details on the generators.

	.. _`random generator: seed`:

	seed
	The seed for random generation. This should be a single number for
	:code:`mixmax` and the name of a state file for :code:`ranlux64`.

	.. _`analysis`:

	analysis
	Name and Setting for the event analyses; either a custom
	analysis plugin or Rivet. For the first the :code:`plugin` sub-entry
	should be set to the analysis file path. All further entries are passed on
	to the analysis. To use Rivet a list of Rivet analyses have to be
	given in :code:`rivet` and prefix for the yoda file has to be set
	through :code:`output`. See :ref:`Writing custom analyses` for details.

	.. _`Higgs coupling`:

	Higgs coupling
	This collects a number of settings concerning the effective coupling
	of the Higgs boson to gluons. This is only relevant for the
	production process of a Higgs boson with jets and only supported if
	HEJ 2 was compiled with `QCDLoop
	<https://github.com/scarrazza/qcdloop>`_ support.

	.. _`Higgs coupling: use impact factors`:

	use impact factors
	Whether to use impact factors for the coupling to the most forward
	and most backward partons. Impact factors imply the infinite
	top-quark mass limit.

	.. _`Higgs coupling: mt`:

	mt
	The value of the top-quark mass in GeV. If this is not specified,
	the limit of an infinite mass is taken.

	.. _`Higgs coupling: include bottom`:

	include bottom
	Whether to include the Higgs coupling to bottom quarks.

	.. _`Higgs coupling: mb`:

	mb
	The value of the bottom-quark mass in GeV. Only used for the Higgs
	coupling, external bottom-quarks are always assumed to be massless.

	Advanced Settings
	~~~~~~~~~~~~~~~~~

	All of the following settings are optional. Please do not set any of the
	following options, unless you know exactly what you are doing. The default
	behaviour gives the most reliable results for a wide range of observables.

	.. _`regulator parameter`:

	regulator parameter
	Slicing parameter to regularise the subtraction term, called :math:`\lambda`
	in `arxiv:1706.01002 <https://arxiv.org/abs/1706.01002>`_. Default is 0.2
	diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
	index 4298e54..96caa71 100644
	--- a/src/YAMLreader.cc
	+++ b/src/YAMLreader.cc
	@@ -1,478 +1,480 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/YAMLreader.hh"

	#include <algorithm>
	#include <iostream>
	#include <limits>
	#include <map>
	#include <string>
	#include <unordered_map>
	#include <vector>

	#include <dlfcn.h>

	#include "HEJ/ScaleFunction.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/output_formats.hh"
	#include "HEJ/Constants.hh"

	namespace HEJ{
	class Event;

	namespace{
	//! 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 = {
	"trials", "min extparton pt", "max ext soft pt fraction",
	- "FKL", "unordered", "extremal qqx", "central qqx", "non-resummable",
	"scales", "scale factors", "max scale ratio", "import scales",
	"log correction", "event output", "analysis", "regulator parameter"
	};
	// add subnodes to "supported" - the assigned value is irrelevant
	for(auto && opt: opts) supported[opt] = "";
	for(auto && jet_opt: {"min pt", "algorithm", "R"}){
	supported["resummation jets"][jet_opt] = "";
	supported["fixed order jets"][jet_opt] = "";
	}
	for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
	supported["Higgs coupling"][opt] = "";
	}
	for(auto && opt: {"name", "seed"}){
	supported["random generator"][opt] = "";
	}
	+ for(auto && opt: {"FKL", "unordered", "extremal qqx", "central qqx", "non-resummable"}){
	+ supported["event treatment"][opt] = "";
	+ }
	return supported;
	}();
	return supported;
	}

	fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){
	using namespace fastjet;
	static const std::map<std::string, fastjet::JetAlgorithm> known = {
	{"kt", kt_algorithm},
	{"cambridge", cambridge_algorithm},
	{"antikt", antikt_algorithm},
	{"cambridge for passive", cambridge_for_passive_algorithm},
	{"plugin", plugin_algorithm}
	};
	const auto res = known.find(algo);
	if(res == known.end()){
	- throw std::invalid_argument("Unknown jet algorithm " + algo);
	+ throw std::invalid_argument("Unknown jet algorithm \"" + algo + "\"");
	}
	return res->second;
	}

	EventTreatment to_EventTreatment(std::string const & name){
	static const std::map<std::string, EventTreatment> known = {
	{"reweight", EventTreatment::reweight},
	{"keep", EventTreatment::keep},
	{"discard", EventTreatment::discard}
	};
	const auto res = known.find(name);
	if(res == known.end()){
	- throw std::invalid_argument("Unknown event treatment " + name);
	+ throw std::invalid_argument("Unknown event treatment \"" + name + "\"");
	}
	return res->second;
	}

	} // namespace anonymous

	namespace detail{
	void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){
	setting = to_JetAlgorithm(yaml.as<std::string>());
	}

	void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){
	setting = to_EventTreatment(yaml.as<std::string>());
	}

	void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){
	setting = to_ParticleID(yaml.as<std::string>());
	}
	} // namespace detail

	JetParameters get_jet_parameters(
	YAML::Node const & node,
	std::string const & entry
	){
	assert(node);
	JetParameters result;
	fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm;
	double R;
	set_from_yaml_if_defined(jet_algo, node, entry, "algorithm");
	set_from_yaml(R, node, entry, "R");
	result.def = fastjet::JetDefinition{jet_algo, R};
	set_from_yaml(result.min_pt, node, entry, "min pt");
	return result;
	}

	RNGConfig to_RNGConfig(
	YAML::Node const & node,
	std::string const & entry
	){
	assert(node);
	RNGConfig result;
	set_from_yaml(result.name, node, entry, "name");
	set_from_yaml_if_defined(result.seed, node, entry, "seed");
	return result;
	}

	HiggsCouplingSettings get_Higgs_coupling(
	YAML::Node const & node,
	std::string const & entry
	){
	assert(node);
	static constexpr double mt_max = 2e4;
	#ifndef HEJ_BUILD_WITH_QCDLOOP
	if(node[entry]){
	throw std::invalid_argument{
	"Higgs coupling settings require building HEJ 2 "
	"with QCDloop support"
	};
	}
	#endif
	HiggsCouplingSettings settings;
	set_from_yaml_if_defined(settings.mt, node, entry, "mt");
	set_from_yaml_if_defined(settings.mb, node, entry, "mb");
	set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom");
	set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors");
	if(settings.use_impact_factors){
	if(settings.mt != std::numeric_limits<double>::infinity()){
	throw std::invalid_argument{
	"Conflicting settings: "
	"impact factors may only be used in the infinite top mass limit"
	};
	}
	}
	else{
	// huge values of the top mass are numerically unstable
	settings.mt = std::min(settings.mt, mt_max);
	}
	return settings;
	}

	FileFormat to_FileFormat(std::string const & name){
	static const std::map<std::string, FileFormat> known = {
	{"Les Houches", FileFormat::Les_Houches},
	{"HepMC", FileFormat::HepMC},
	{"HepMC2", FileFormat::HepMC2},
	{"HepMC3", FileFormat::HepMC3}
	};
	const auto res = known.find(name);
	if(res == known.end()){
	- throw std::invalid_argument("Unknown file format " + name);
	+ throw std::invalid_argument("Unknown file format \"" + name + "\"");
	}
	return res->second;
	}

	std::string extract_suffix(std::string const & filename){
	size_t separator = filename.rfind('.');
	if(separator == filename.npos) return {};
	return filename.substr(separator + 1);
	}

	FileFormat format_from_suffix(std::string const & filename){
	const std::string suffix = extract_suffix(filename);
	if(suffix == "lhe") return FileFormat::Les_Houches;
	if(suffix == "hepmc") return FileFormat::HepMC;
	if(suffix == "hepmc3") return FileFormat::HepMC3;
	if(suffix == "hepmc2") return FileFormat::HepMC2;
	throw std::invalid_argument{
	- "Can't determine format for output file " + filename
	+ "Can't determine format for output file \"" + filename + "\""
	};
	}

	void assert_all_options_known(
	YAML::Node const & conf, YAML::Node const & supported
	){
	if(!conf.IsMap()) return;
	if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"};
	for(auto const & entry: conf){
	const auto name = entry.first.as<std::string>();
	if(! supported[name]) throw unknown_option{name};
	/* check sub-options, e.g. 'resummation jets: min pt'
	* we don't check analysis sub-options
	* those depend on the analysis being used and should be checked there
	* similar for "import scales"
	*/
	if(name != "analysis" && name != "import scales"){
	try{
	assert_all_options_known(conf[name], supported[name]);
	}
	catch(unknown_option const & ex){
	throw unknown_option{name + ": " + ex.what()};
	}
	catch(invalid_type const & ex){
	throw invalid_type{name + ": " + ex.what()};
	}
	}
	}
	}

	} // namespace HEJ

	namespace YAML {

	Node convert<HEJ::OutputFile>::encode(HEJ::OutputFile const & outfile) {
	Node node;
	node[to_string(outfile.format)] = outfile.name;
	return node;
	};

	bool convert<HEJ::OutputFile>::decode(Node const & node, HEJ::OutputFile & out) {
	switch(node.Type()){
	case NodeType::Map: {
	YAML::const_iterator it = node.begin();
	out.format = HEJ::to_FileFormat(it->first.as<std::string>());
	out.name = it->second.as<std::string>();
	return true;
	}
	case NodeType::Scalar:
	out.name = node.as<std::string>();
	out.format = HEJ::format_from_suffix(out.name);
	return true;
	default:
	return false;
	}
	}
	} // namespace YAML

	namespace HEJ{

	namespace detail{
	void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){
	setting = yaml.as<OutputFile>();
	}
	}

	namespace{
	void update_fixed_order_jet_parameters(
	JetParameters & fixed_order_jets, YAML::Node const & yaml
	){
	if(!yaml["fixed order jets"]) return;
	set_from_yaml_if_defined(
	fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt"
	);
	fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm();
	set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm");
	double R = fixed_order_jets.def.R();
	set_from_yaml_if_defined(R, yaml, "fixed order jets", "R");
	fixed_order_jets.def = fastjet::JetDefinition{algo, R};
	}

	// like std::stod, but throw if not the whole string can be converted
	double to_double(std::string const & str){
	std::size_t pos;
	const double result = std::stod(str, &pos);
	if(pos < str.size()){
	throw std::invalid_argument(str + " is not a valid double value");
	}
	return result;
	}

	using EventScale = double (*)(Event const &);

	void import_scale_functions(
	std::string const & file,
	std::vector<std::string> const & scale_names,
	std::unordered_map<std::string, EventScale> & known
	) {
	auto handle = dlopen(file.c_str(), RTLD_NOW);
	char * error = dlerror();
	if(error != nullptr) throw std::runtime_error{error};

	for(auto const & scale: scale_names) {
	void * sym = dlsym(handle, scale.c_str());
	error = dlerror();
	if(error != nullptr) throw std::runtime_error{error};
	known.emplace(scale, reinterpret_cast<EventScale>(sym));
	}
	}

	auto get_scale_map(
	YAML::Node const & yaml
	) {
	std::unordered_map<std::string, EventScale> scale_map;
	scale_map.emplace("H_T", H_T);
	scale_map.emplace("max jet pperp", max_jet_pt);
	scale_map.emplace("jet invariant mass", jet_invariant_mass);
	scale_map.emplace("m_j1j2", m_j1j2);
	if(yaml["import scales"]) {
	if(! yaml["import scales"].IsMap()) {
	throw invalid_type{"Entry 'import scales' is not a map"};
	}
	for(auto const & import: yaml["import scales"]) {
	const auto file = import.first.as<std::string>();
	const auto scale_names =
	import.second.IsSequence()
	?import.second.as<std::vector<std::string>>()
	:std::vector<std::string>{import.second.as<std::string>()};
	import_scale_functions(file, scale_names, scale_map);
	}
	}
	return scale_map;
	}

	// simple (as in non-composite) scale functions
	/**
	* An example for a simple scale function would be H_T,
	* H_T/2 is then composite (take H_T and then divide by 2)
	*/
	ScaleFunction parse_simple_ScaleFunction(
	std::string const & scale_fun,
	std::unordered_map<std::string, EventScale> const & known
	) {
	assert(
	scale_fun.empty() \|\|
	(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
	);
	const auto it = known.find(scale_fun);
	if(it != end(known)) return {it->first, it->second};
	try{
	const double scale = to_double(scale_fun);
	return {scale_fun, FixedScale{scale}};
	} catch(std::invalid_argument const &){}
	- throw std::invalid_argument{"Unknown scale choice: " + scale_fun};
	+ throw std::invalid_argument{"Unknown scale choice: \"" + scale_fun + "\""};
	}

	std::string trim_front(std::string const & str){
	const auto new_begin = std::find_if(
	begin(str), end(str), [](char c){ return ! std::isspace(c); }
	);
	return std::string(new_begin, end(str));
	}

	std::string trim_back(std::string str){
	size_t pos = str.size() - 1;
	// use guaranteed wrap-around behaviour to check whether we have
	// traversed the whole string
	for(; pos < str.size() && std::isspace(str[pos]); --pos) {}
	str.resize(pos + 1); // note that pos + 1 can be 0
	return str;
	}

	ScaleFunction parse_ScaleFunction(
	std::string const & scale_fun,
	std::unordered_map<std::string, EventScale> const & known
	){
	assert(
	scale_fun.empty() \|\|
	(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
	);
	// parse from right to left => a/b/c gives (a/b)/c
	const size_t delim = scale_fun.find_last_of("*/");
	if(delim == scale_fun.npos){
	return parse_simple_ScaleFunction(scale_fun, known);
	}
	const std::string first = trim_back(std::string{scale_fun, 0, delim});
	const std::string second = trim_front(std::string{scale_fun, delim+1});
	if(scale_fun[delim] == '/'){
	return parse_ScaleFunction(first, known)
	/ parse_ScaleFunction(second, known);
	}
	else{
	assert(scale_fun[delim] == '*');
	return parse_ScaleFunction(first, known)
	* parse_ScaleFunction(second, known);
	}
	}

	EventTreatMap get_event_treatment(
	- YAML::Node const & yaml
	+ YAML::Node const & node, std::string const & entry
	){
	using namespace event_type;
	EventTreatMap treat {
	- {no_2_jets, EventTreatment::discard},
	- {bad_final_state, EventTreatment::discard},
	- {FKL, EventTreatment::reweight},
	- {unob, EventTreatment::keep},
	- {unof, EventTreatment::keep},
	- {qqxexb, EventTreatment::keep},
	- {qqxexf, EventTreatment::keep},
	- {qqxmid, EventTreatment::keep},
	- {non_resummable, EventTreatment::keep}
	- };
	- set_from_yaml(treat.at(FKL), yaml, "FKL");
	- set_from_yaml(treat.at(unob), yaml, "unordered");
	+ {no_2_jets, EventTreatment::discard},
	+ {bad_final_state, EventTreatment::discard},
	+ {FKL, EventTreatment::discard},
	+ {unob, EventTreatment::discard},
	+ {unof, EventTreatment::discard},
	+ {qqxexb, EventTreatment::discard},
	+ {qqxexf, EventTreatment::discard},
	+ {qqxmid, EventTreatment::discard},
	+ {non_resummable, EventTreatment::discard}
	+ };
	+ set_from_yaml(treat.at(FKL), node, entry, "FKL");
	+ set_from_yaml(treat.at(unob), node, entry, "unordered");
	treat.at(unof) = treat.at(unob);
	- set_from_yaml(treat.at(qqxexb), yaml, "extremal qqx");
	- set_from_yaml(treat.at(qqxexf), yaml, "extremal qqx");
	- set_from_yaml(treat.at(qqxmid), yaml, "central qqx");
	- set_from_yaml(treat.at(non_resummable), yaml, "non-resummable");
	+ set_from_yaml(treat.at(qqxexb), node, entry, "extremal qqx");
	+ treat.at(qqxexf) = treat.at(qqxexb);
	+ set_from_yaml(treat.at(qqxmid), node, entry, "central qqx");
	+ set_from_yaml(treat.at(non_resummable), node, entry, "non-resummable");
	if(treat[non_resummable] == EventTreatment::reweight){
	throw std::invalid_argument{"Cannot reweight Fixed Order events"};
	}
	return treat;
	}

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

	Config config;
	config.resummation_jets = get_jet_parameters(yaml, "resummation jets");
	config.fixed_order_jets = config.resummation_jets;
	update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
	set_from_yaml(config.min_extparton_pt, yaml, "min extparton pt");
	// Sets the standard value, then changes this if defined
	config.regulator_lambda=CLAMBDA;
	set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter");
	config.max_ext_soft_pt_fraction = 0.1;
	set_from_yaml_if_defined(
	config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
	);
	set_from_yaml(config.trials, yaml, "trials");
	set_from_yaml(config.log_correction, yaml, "log correction");
	- config.treat = get_event_treatment(yaml);
	+ config.treat = get_event_treatment(yaml, "event treatment");
	set_from_yaml_if_defined(config.output, yaml, "event output");
	config.rng = to_RNGConfig(yaml, "random generator");
	set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
	config.scales = to_ScaleConfig(yaml);
	config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
	return config;
	}

	} // namespace anonymous

	ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
	ScaleConfig config;
	auto scale_funs = get_scale_map(yaml);
	std::vector<std::string> scales;
	set_from_yaml(scales, yaml, "scales");
	config.base.reserve(scales.size());
	std::transform(
	begin(scales), end(scales), std::back_inserter(config.base),
	[scale_funs](auto const & entry){
	return parse_ScaleFunction(entry, scale_funs);
	}
	);
	set_from_yaml_if_defined(config.factors, yaml, "scale factors");
	config.max_ratio = std::numeric_limits<double>::infinity();
	set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio");
	return config;
	}

	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 HEJ
	diff --git a/t/ME_data/config_mt.yml b/t/ME_data/config_mt.yml
	index 3e81c97..fd871f3 100644
	--- a/t/ME_data/config_mt.yml
	+++ b/t/ME_data/config_mt.yml
	@@ -1,30 +1,31 @@
	trials: 1

	min extparton pt: 30

	resummation jets:
	min pt: 30
	algorithm: antikt
	R: 0.4

	fixed order jets:
	min pt: 30

	-FKL: reweight
	-unordered: reweight
	-extremal qqx: discard
	-central qqx: discard
	-non-resummable: discard
	+event treatment:
	+ FKL: reweight
	+ unordered: reweight
	+ extremal qqx: discard
	+ central qqx: discard
	+ non-resummable: discard

	scales: 125

	log correction: false

	random generator:
	name: mixmax
	seed: 1

	Higgs coupling:
	use impact factors: false
	mt: 174
	include bottom: false
	diff --git a/t/ME_data/config_mtinf.yml b/t/ME_data/config_mtinf.yml
	index e99541a..0400893 100644
	--- a/t/ME_data/config_mtinf.yml
	+++ b/t/ME_data/config_mtinf.yml
	@@ -1,25 +1,26 @@
	trials: 1

	min extparton pt: 30

	resummation jets:
	min pt: 30
	algorithm: antikt
	R: 0.4

	fixed order jets:
	min pt: 30

	-FKL: reweight
	-unordered: reweight
	-extremal qqx: discard
	-central qqx: discard
	-non-resummable: discard
	+event treatment:
	+ FKL: reweight
	+ unordered: reweight
	+ extremal qqx: discard
	+ central qqx: discard
	+ non-resummable: discard

	scales: 125

	log correction: false

	random generator:
	name: mixmax
	seed: 1
	diff --git a/t/ME_data/config_mtmb.yml b/t/ME_data/config_mtmb.yml
	index 26a0760..3c7ce17 100644
	--- a/t/ME_data/config_mtmb.yml
	+++ b/t/ME_data/config_mtmb.yml
	@@ -1,31 +1,32 @@
	trials: 1

	min extparton pt: 30

	resummation jets:
	min pt: 30
	algorithm: antikt
	R: 0.4

	fixed order jets:
	min pt: 30

	-FKL: reweight
	-unordered: reweight
	-extremal qqx: discard
	-central qqx: discard
	-non-resummable: discard
	+event treatment:
	+ FKL: reweight
	+ unordered: reweight
	+ extremal qqx: discard
	+ central qqx: discard
	+ non-resummable: discard

	scales: 125

	log correction: false

	random generator:
	name: mixmax
	seed: 1

	Higgs coupling:
	use impact factors: false
	mt: 174
	include bottom: true
	mb: 4.7
	diff --git a/t/ME_data/config_w_ME.yml b/t/ME_data/config_w_ME.yml
	index b605ec1..c389201 100644
	--- a/t/ME_data/config_w_ME.yml
	+++ b/t/ME_data/config_w_ME.yml
	@@ -1,25 +1,26 @@
	trials: 1

	min extparton pt: 30

	resummation jets:
	min pt: 30
	algorithm: antikt
	R: 0.4

	fixed order jets:
	min pt: 30

	-FKL: reweight
	-unordered: reweight
	-extremal qqx: reweight
	-central qqx: reweight
	-non-resummable: discard
	+event treatment:
	+ FKL: reweight
	+ unordered: reweight
	+ extremal qqx: reweight
	+ central qqx: reweight
	+ non-resummable: discard

	scales: 125

	log correction: false

	random generator:
	name: mixmax
	seed: 1
	diff --git a/t/jet_config.yml b/t/jet_config.yml
	index 0a021f9..9309e04 100644
	--- a/t/jet_config.yml
	+++ b/t/jet_config.yml
	@@ -1,27 +1,28 @@
	trials: 10

	min extparton pt: 30

	resummation jets:
	min pt: 35
	algorithm: antikt
	R: 0.4

	fixed order jets:
	min pt: 30

	-FKL: reweight
	-unordered: keep
	-extremal qqx: discard
	-central qqx: keep
	-non-resummable: keep
	+event treatment:
	+ FKL: reweight
	+ unordered: keep
	+ extremal qqx: discard
	+ central qqx: keep
	+ non-resummable: keep

	log correction: false

	scales: 91.188

	random generator:
	name: ranlux64

	event output:
	- tst.lhe
	diff --git a/t/jet_config_with_import.yml b/t/jet_config_with_import.yml
	index 1d3e04f..924cad3 100644
	--- a/t/jet_config_with_import.yml
	+++ b/t/jet_config_with_import.yml
	@@ -1,30 +1,31 @@
	trials: 10

	min extparton pt: 30

	resummation jets:
	min pt: 35
	algorithm: antikt
	R: 0.4

	fixed order jets:
	min pt: 30

	-FKL: reweight
	-unordered: discard
	-extremal qqx: discard
	-central qqx: discard
	-non-resummable: discard
	+event treatment:
	+ FKL: reweight
	+ unordered: discard
	+ extremal qqx: discard
	+ central qqx: discard
	+ non-resummable: discard

	log correction: false

	scales: softest_jet_pt

	event output:
	- tst.lhe

	random generator:
	name: ranlux64

	import scales:
	./libscales.so: softest_jet_pt

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:27 PM (16 h, 5 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486594
Default Alt Text: (56 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions