No OneTemporary
Actions

Size

58 KB

Subscribers

None

View Options

	diff --git a/FixedOrderGen/.clang-tidy b/FixedOrderGen/.clang-tidy
	index c9f60aa..d2b65fd 100644
	--- a/FixedOrderGen/.clang-tidy
	+++ b/FixedOrderGen/.clang-tidy
	@@ -1,46 +1,44 @@
	Checks: '
	*,
	- -clang-analyzer-*,
	- -cppcoreguidelines-*,
	-google-*,
	-llvm-*,
	-misc-*,
	-modernize-*,
	-performance-*,
	-readability-*,
	-cert-dcl21-cpp,
	-clang-analyzer-optin.cplusplus.VirtualCall,
	-cppcoreguidelines-avoid-magic-numbers,
	-cppcoreguidelines-pro-bounds-array-to-pointer-decay,
	-cppcoreguidelines-pro-bounds-constant-array-index,
	-fuchsia*,
	-google-build-using-namespace,
	-google-explicit-constructor,
	-google-readability*,
	-google-runtime-int,
	-google-runtime-references,
	-hicpp*,
	-llvm-header-guard,
	-modernize-use-trailing-return-type,
	-readability-braces-around-statements,
	-readability-uppercase-literal-suffix,
	'
	CheckOptions:
	- key: cppcoreguidelines-special-member-functions.AllowSoleDefaultDtor
	value: '1'
	- key: llvm-namespace-comment.ShortNamespaceLines
	value: '10'
	- key: misc-non-private-member-variables-in-classes.IgnoreClassesWithAllMemberVariablesBeingPublic
	value: '1'
	- key: performance-move-const-arg.CheckTriviallyCopyableMove
	value: '0'
	- key: readability-identifier-naming.GlobalConstantCase
	value: UPPER_CASE
	- key: readability-identifier-naming.GlobalVariableCase
	value: UPPER_CASE
	- key: readability-identifier-naming.PrivateMemberSuffix
	value: '_'
	- key: readability-implicit-bool-conversion.AllowIntegerConditions
	value: '1'
	- key: readability-magic-numbers.IgnoreAllFloatingPointValues
	value: '1'
	diff --git a/FixedOrderGen/include/Beam.hh b/FixedOrderGen/include/Beam.hh
	index 4054083..57b55fe 100644
	--- a/FixedOrderGen/include/Beam.hh
	+++ b/FixedOrderGen/include/Beam.hh
	@@ -1,19 +1,19 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <array>

	#include "HEJ/PDG_codes.hh"

	namespace HEJFOG {
	struct Beam{
	- double energy;
	+ double energy{};
	std::array<HEJ::ParticleID, 2> particles{{
	HEJ::pid::proton, HEJ::pid::proton
	}};
	};
	}
	diff --git a/FixedOrderGen/include/Decay.hh b/FixedOrderGen/include/Decay.hh
	index 2c92262..aee1c34 100644
	--- a/FixedOrderGen/include/Decay.hh
	+++ b/FixedOrderGen/include/Decay.hh
	@@ -1,27 +1,27 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <unordered_map>
	#include <vector>

	#include "HEJ/PDG_codes.hh"

	namespace HEJFOG {
	struct Decay{
	std::vector<HEJ::ParticleID> products;
	- double branching_ratio;
	+ double branching_ratio{};
	};
	#if !defined(__clang__) && defined(__GNUC__) && (__GNUC__ < 6)
	// gcc version < 6 explicitly needs hash function for enum
	// see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=60970
	using ParticlesDecayMap
	= std::unordered_map<HEJ::ParticleID, std::vector<Decay>, std::hash<int>>;
	#else
	using ParticlesDecayMap
	= std::unordered_map<HEJ::ParticleID, std::vector<Decay>>;
	#endif
	}
	diff --git a/FixedOrderGen/include/JetParameters.hh b/FixedOrderGen/include/JetParameters.hh
	index 6ea2711..f648d2c 100644
	--- a/FixedOrderGen/include/JetParameters.hh
	+++ b/FixedOrderGen/include/JetParameters.hh
	@@ -1,19 +1,19 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include "fastjet/JetDefinition.hh"

	#include "HEJ/optional.hh"

	namespace HEJFOG {
	struct JetParameters{
	fastjet::JetDefinition def;
	- double min_pt;
	- double max_y;
	+ double min_pt{};
	+ double max_y{};
	HEJ::optional<double> peak_pt;
	};
	}
	diff --git a/FixedOrderGen/include/Process.hh b/FixedOrderGen/include/Process.hh
	index 5ae39e9..e6c950c 100644
	--- a/FixedOrderGen/include/Process.hh
	+++ b/FixedOrderGen/include/Process.hh
	@@ -1,23 +1,23 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <array>
	#include <vector>

	#include "HEJ/PDG_codes.hh"
	#include "HEJ/optional.hh"

	namespace HEJFOG {
	struct Process{
	//! @internal pair to match Event::incoming
	- std::array<HEJ::ParticleID, 2> incoming;
	- std::size_t njets;
	+ std::array<HEJ::ParticleID, 2> incoming{};
	+ std::size_t njets{};
	HEJ::optional<HEJ::ParticleID> boson;
	std::vector<HEJ::ParticleID> boson_decay;
	};

	}
	diff --git a/FixedOrderGen/include/config.hh b/FixedOrderGen/include/config.hh
	index 69c62e4..36f1259 100644
	--- a/FixedOrderGen/include/config.hh
	+++ b/FixedOrderGen/include/config.hh
	@@ -1,50 +1,50 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <cstddef>
	#include <string>
	#include <vector>

	#include "HEJ/Config.hh"
	#include "HEJ/EWConstants.hh"
	#include "HEJ/Fraction.hh"
	#include "HEJ/HiggsCouplingSettings.hh"
	#include "HEJ/optional.hh"
	#include "HEJ/output_formats.hh"

	#include "yaml-cpp/yaml.h"

	#include "Beam.hh"
	#include "Decay.hh"
	#include "JetParameters.hh"
	#include "Process.hh"
	#include "Subleading.hh"
	#include "UnweightSettings.hh"

	namespace HEJFOG {

	struct Config{
	Process process;
	- std::size_t events;
	+ std::size_t events{};
	JetParameters jets;
	Beam beam;
	- int pdf_id;
	- HEJ::Fraction<double> subleading_fraction;
	+ int pdf_id{};
	+ HEJ::Fraction<double> subleading_fraction{};
	Subleading subleading_channels; //! < see HEJFOG::Subleading
	ParticlesDecayMap particle_decays;
	std::vector<YAML::Node> analyses_parameters;
	HEJ::ScaleConfig scales;
	std::vector<HEJ::OutputFile> output;
	HEJ::RNGConfig rng;
	HEJ::HiggsCouplingSettings Higgs_coupling;
	HEJ::EWConstants ew_parameters;
	HEJ::optional<UnweightSettings> unweight;
	};

	Config load_config(std::string const & config_file);

	}
	diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
	index 935b1eb..13654c0 100644
	--- a/FixedOrderGen/src/EventGenerator.cc
	+++ b/FixedOrderGen/src/EventGenerator.cc
	@@ -1,100 +1,101 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "EventGenerator.hh"

	#include <cstddef>
	#include <utility>

	#include "HEJ/Config.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/EWConstants.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/HiggsCouplingSettings.hh"

	#include "Process.hh"
	#include "Beam.hh"
	#include "JetParameters.hh"
	#include "PhaseSpacePoint.hh"

	namespace HEJFOG {
	EventGenerator::EventGenerator(
	Process process,
	Beam beam,
	HEJ::ScaleGenerator scale_gen,
	JetParameters jets,
	int pdf_id,
	double subl_chance,
	Subleading subl_channels,
	ParticlesDecayMap particle_decays,
	HEJ::HiggsCouplingSettings Higgs_coupling,
	HEJ::EWConstants ew_parameters,
	std::shared_ptr<HEJ::RNG> ran
	):
	pdf_{pdf_id, beam.particles[0], beam.particles[1]},
	ME_{
	[this](double mu){ return pdf_.Halphas(mu); },
	HEJ::MatrixElementConfig{
	false,
	std::move(Higgs_coupling),
	ew_parameters
	}
	},
	scale_gen_{std::move(scale_gen)},
	process_{std::move(process)},
	jets_{std::move(jets)},
	beam_{std::move(beam)},
	+ status_{Status::good},
	subl_chance_{subl_chance},
	subl_channels_{subl_channels},
	particle_decays_{std::move(particle_decays)},
	ew_parameters_{ew_parameters},
	ran_{std::move(ran)}
	{
	}

	HEJ::optional<HEJ::Event> EventGenerator::gen_event(){
	HEJFOG::PhaseSpacePoint psp{
	process_,
	jets_,
	pdf_, beam_.energy,
	subl_chance_, subl_channels_,
	particle_decays_,
	ew_parameters_,
	*ran_
	};
	status_ = psp.status();
	if(status_ != Status::good) return {};

	HEJ::Event ev = scale_gen_(
	HEJ::Event{
	to_EventData( std::move(psp) ).cluster( jets_.def, jets_.min_pt)
	}
	);
	if(!is_resummable(ev.type())){
	throw HEJ::not_implemented("Tried to generate a event type, "
	"which is not yet implemented in HEJ.");
	}

	ev.generate_colours(*ran_);

	const double shat = HEJ::shat(ev);
	const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy);
	const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy);

	// evaluate matrix element
	ev.parameters() = ME_.tree(ev)/(shatshat);
	// and PDFs
	ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
	ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
	for(std::size_t i = 0; i < ev.variations().size(); ++i){
	auto & var = ev.variations(i);
	var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type);
	var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type);
	}
	return ev;
	}

	}
	diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
	index 0d4d772..dee306e 100644
	--- a/FixedOrderGen/src/PhaseSpacePoint.cc
	+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
	@@ -1,974 +1,974 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "PhaseSpacePoint.hh"

	#include <algorithm>
	#include <cassert>
	#include <cmath>
	#include <cstdlib>
	#include <iostream>
	#include <iterator>
	#include <limits>
	#include <tuple>
	#include <type_traits>
	#include <utility>

	#include "fastjet/ClusterSequence.hh"

	#include "HEJ/Constants.hh"
	#include "HEJ/EWConstants.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/kinematics.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/RNG.hh"
	#include "HEJ/utility.hh"

	#include "JetParameters.hh"
	#include "Process.hh"

	namespace {
	static_assert(
	std::numeric_limits<double>::has_quiet_NaN,
	"no quiet NaN for double"
	);
	constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
	} // namespace anonymous

	namespace HEJFOG {
	HEJ::Event::EventData to_EventData(PhaseSpacePoint psp){
	//! @TODO Same function already in HEJ
	HEJ::Event::EventData result;
	result.incoming = std::move(psp).incoming_; // NOLINT(bugprone-use-after-move)
	result.outgoing = std::move(psp).outgoing_; // NOLINT(bugprone-use-after-move)
	// technically Event::EventData doesn't have to be sorted,
	// but PhaseSpacePoint should be anyway
	assert(
	std::is_sorted(
	begin(result.outgoing), end(result.outgoing),
	HEJ::rapidity_less{}
	)
	);
	assert(result.outgoing.size() >= 2);
	result.decays = std::move(psp).decays_; // NOLINT(bugprone-use-after-move)
	result.parameters.central = {NaN, NaN, psp.weight()}; // NOLINT(bugprone-use-after-move)
	return result;
	}

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

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

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

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

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

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

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

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

	namespace {
	bool can_swap_to_uno(
	HEJ::Particle const & p1, HEJ::Particle const & p2
	){
	assert(is_parton(p1) && is_parton(p2));
	return p1.type != HEJ::pid::gluon
	&& p2.type == HEJ::pid::gluon;
	}

	size_t count_gluons(PhaseSpacePoint::ConstPartonIterator first,
	PhaseSpacePoint::ConstPartonIterator last
	){
	return std::count_if(first, last, [](HEJ::Particle const & p)
	{return p.type == HEJ::pid::gluon;});
	}

	/** assumes FKL configurations between first and last,
	* else there can be a quark in a non-extreme position
	* e.g. uno configuration gqg would pass
	*/
	Subleading possible_qqx(
	PhaseSpacePoint::ConstPartonIterator first,
	PhaseSpacePoint::ConstReversePartonIterator last
	){
	using namespace subleading;
	assert( std::distance( first,last.base() )>2 );
	Subleading channels = ALL;
	channels.reset(eqqx);
	channels.reset(cqqx);
	auto const ngluon = count_gluons(first,last.base());
	if(ngluon < 2) return channels;
	if(first->type==HEJ::pid::gluon \|\| last->type==HEJ::pid::gluon){
	channels.set(eqqx);
	}
	if(std::distance(first,last.base())>=4){
	channels.set(cqqx);
	}
	return channels;
	}

	template<class PartonIt, class OutIt>
	bool uno_possible(
	PartonIt first_parton, OutIt first_out
	){
	using namespace HEJ;
	// Special case: Higgs can not be outside of uno
	if(first_out->type == pid::Higgs
	\|\| std::next(first_out)->type==pid::Higgs){
	return false;
	}
	// decide what kind of subleading process is allowed
	return can_swap_to_uno( first_parton, std::next(first_parton) );
	}

	bool is_AWZ_proccess(Process const & proc){
	return proc.boson && HEJ::is_AWZ_boson(*proc.boson);
	}

	bool is_up_type(HEJ::Particle const & part){
	return is_anyquark(part) && !(std::abs(part.type)%2);
	}
	bool is_down_type(HEJ::Particle const & part){
	return is_anyquark(part) && std::abs(part.type)%2;
	}
	bool can_couple_to_W(
	HEJ::Particle const & part, HEJ::pid::ParticleID const W_id
	){
	const int W_charge = W_id>0?1:-1;
	return std::abs(part.type)<5
	&& ( (W_charge*part.type > 0 && is_up_type(part))
	\|\| (W_charge*part.type < 0 && is_down_type(part)) );
	}

	Subleading ensure_AWZ(
	double & subl_chance, bool & allow_strange,
	HEJ::ParticleID const boson,
	PhaseSpacePoint::ConstPartonIterator first_parton,
	PhaseSpacePoint::ConstPartonIterator last_parton
	){
	auto channels = subleading::ALL;
	if(std::none_of(first_parton, last_parton,
	[&boson](HEJ::Particle const & p){
	return can_couple_to_W(p, boson);})) {
	// enforce qqx if A/W/Z can't couple somewhere else
	// this is ensured to work through filter_partons in reconstruct_incoming
	channels.reset(subleading::uno);
	assert(channels.any());
	subl_chance = 1.;
	// strange not allowed for W
	if(std::abs(boson)== HEJ::pid::Wp) allow_strange = false;
	}
	return channels;
	}
	}

	void PhaseSpacePoint::turn_to_subl(
	Subleading const channels,
	bool const can_be_uno_backward, bool const can_be_uno_forward,
	bool const allow_strange,
	HEJ::RNG & ran
	){
	double const nchannels = channels.count();
	double const step = 1./nchannels;
	weight_*=nchannels;
	unsigned selected = subleading::first;
	double rnd = nchannels>1?ran.flat():0.;
	// @TODO optimise this sampling
	for(; selected<=subleading::last; ++selected){
	assert(rnd>=0);
	if(channels[selected]){
	if(rnd<step) break;
	rnd-=step;
	}
	}
	switch(selected){
	case subleading::uno:
	return turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
	case subleading::cqqx:
	return turn_to_cqqx(allow_strange, ran);
	case subleading::eqqx:
	return turn_to_eqqx(allow_strange, ran);
	default:
	throw std::logic_error{"unreachable"};
	}
	}

	void PhaseSpacePoint::maybe_turn_to_subl(
	double chance, Subleading channels, Process const & proc, HEJ::RNG & ran
	){
	using namespace HEJ;
	if(proc.njets <= 2) return;
	assert(outgoing_.size() >= 2);

	// decide what kind of subleading process is allowed
	bool const can_be_uno_backward = uno_possible(cbegin_partons(),
	outgoing_.cbegin());
	bool const can_be_uno_forward = uno_possible(crbegin_partons(),
	outgoing_.crbegin());
	if(channels[subleading::uno]){
	channels.set(subleading::uno, can_be_uno_backward \|\| can_be_uno_forward);
	}
	channels &= possible_qqx(cbegin_partons(), crbegin_partons());

	bool allow_strange = true;
	if(is_AWZ_proccess(proc)) {
	channels &= ensure_AWZ(chance, allow_strange, *proc.boson,
	cbegin_partons(), cend_partons());
	}

	std::size_t const nchannels = channels.count();
	// no subleading
	if(nchannels==0) return;
	if(ran.flat() >= chance){
	weight_ /= 1 - chance;
	return;
	}
	weight_ /= chance;

	// select channel
	return turn_to_subl( channels, can_be_uno_backward, can_be_uno_forward,
	allow_strange, ran);
	}

	void PhaseSpacePoint::turn_to_uno(
	const bool can_be_uno_backward, const bool can_be_uno_forward,
	HEJ::RNG & ran
	){
	if(!can_be_uno_backward && !can_be_uno_forward) return;
	if(can_be_uno_backward && can_be_uno_forward){
	weight_ *= 2.;
	if(ran.flat() < 0.5){
	return std::swap(begin_partons()->type, std::next(begin_partons())->type);
	}
	return std::swap(rbegin_partons()->type, std::next(rbegin_partons())->type);
	}
	if(can_be_uno_backward){
	return std::swap(begin_partons()->type, std::next(begin_partons())->type);
	}
	assert(can_be_uno_forward);
	std::swap(rbegin_partons()->type, std::next(rbegin_partons())->type);
	}

	//! select flavour of quark
	HEJ::ParticleID PhaseSpacePoint::select_qqx_flavour(
	const bool allow_strange, HEJ::RNG & ran
	){
	const double r1 = 2.*ran.flat()-1.;
	const double max_flavour = allow_strange?HEJ::N_F:HEJ::N_F-1;
	weight_ = max_flavour2;
	double const flavour = HEJ::pid::down + std::floor(std::abs(r1)*max_flavour);
	return static_cast<HEJ::ParticleID>(flavour*(r1<0.?-1:1));
	}

	void PhaseSpacePoint::turn_to_cqqx(const bool allow_strange, HEJ::RNG & ran){
	// we assume all FKL partons to be gluons
	auto first = std::next(begin_partons());
	auto last = std::next(rbegin_partons());
	auto const ng = std::distance(first, last.base());
	if(ng < 2)
	throw std::logic_error("not enough gluons to create qqx");
	auto flavour = select_qqx_flavour(allow_strange, ran);
	// select gluon for switch
	if(ng!=2){
	const double steps = 1./(ng-1.);
	weight_ /= steps;
	for(auto rnd = ran.flat(); rnd>steps; ++first){
	rnd-=steps;
	}
	}
	first->type = flavour;
	std::next(first)->type = anti(flavour);
	}

	void PhaseSpacePoint::turn_to_eqqx(const bool allow_strange, HEJ::RNG & ran){
	/// find first and last gluon in FKL chain
	auto first = begin_partons();
	const bool can_forward = !is_anyquark(*first);
	auto last = rbegin_partons();
	const bool can_backward = !is_anyquark(*last);
	if(std::distance(first, last.base()) < 2)
	throw std::logic_error("not enough gluons to create qqx");
	auto flavour = select_qqx_flavour(allow_strange, ran);
	// select gluon for switch
	if(can_forward && !can_backward){
	first->type = flavour;
	std::next(first)->type = anti(flavour);
	return;
	}
	if(!can_forward && can_backward){
	last->type = flavour;
	std::next(last)->type = anti(flavour);
	return;
	}
	assert(can_forward && can_backward);
	weight_*=2.;
	if(ran.flat()>0.5){
	first->type = flavour;
	std::next(first)->type = anti(flavour);
	return;
	}
	last->type = flavour;
	std::next(last)->type = anti(flavour);
	}

	template<class ParticleMomenta>
	fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum(
	ParticleMomenta const & other_momenta,
	const double mass_square, const double y
	) const {
	std::array<double,2> pt{0.,0.};
	for (auto const & p: other_momenta) {
	pt[0]-= p.px();
	pt[1]-= p.py();
	}

	const double mperp = std::sqrt(pt[0]pt[0]+pt[1]pt[1]+mass_square);
	const double pz=mperp*std::sinh(y);
	const double E=mperp*std::cosh(y);

	return {pt[0], pt[1], pz, E};
	}

	namespace {
	//! adds a particle to target (in correct rapidity ordering)
	//! @returns positon of insertion
	auto insert_particle(std::vector<HEJ::Particle> & target,
	HEJ::Particle && particle
	){
	const auto pos = std::upper_bound(
	begin(target),end(target),particle,HEJ::rapidity_less{}
	);
	target.insert(pos, std::move(particle));
	return pos;
	}
	}

	PhaseSpacePoint::PhaseSpacePoint(
	Process const & proc,
	JetParameters const & jet_param,
	HEJ::PDF & pdf, double E_beam,
	double const subl_chance,
	Subleading subl_channels,
	ParticlesDecayMap const & particle_decays,
	HEJ::EWConstants const & ew_parameters,
	HEJ::RNG & ran
	){
	assert(proc.njets >= 2);
	status_ = Status::good;
	weight_ = 1;

	// ensure that all setting are consistent
	if(subl_chance == 0.)
	subl_channels.reset();

	const std::size_t nout = proc.njets + (proc.boson?1:0)
	+ proc.boson_decay.size();
	outgoing_.reserve(nout);

	// generate parton momenta
	const bool is_pure_jets = (nout == proc.njets);
	auto partons = gen_LO_partons(
	proc.njets, is_pure_jets, jet_param, E_beam, ran
	);
	// pre fill flavour with gluons
	for(auto&& p_out: partons) {
	outgoing_.emplace_back(HEJ::Particle{HEJ::pid::gluon, std::move(p_out), {}});
	}
	if(status_ != Status::good) return;

	if(proc.boson){ // decay boson
	auto const & boson_prop = ew_parameters.prop(*proc.boson) ;
	auto boson{ gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran) };
	const auto pos{insert_particle(outgoing_, std::move(boson))};
	const size_t boson_idx = std::distance(begin(outgoing_), pos);
	auto const & boson_decay = particle_decays.find(*proc.boson);
	if( !proc.boson_decay.empty() ){ // decay given in proc
	decays_.emplace(
	boson_idx,
	decay_boson(outgoing_[boson_idx], proc.boson_decay, ran)
	);
	} else if( boson_decay != particle_decays.end()
	&& !boson_decay->second.empty() ){ // decay given explicitly
	decays_.emplace(
	boson_idx,
	decay_boson(outgoing_[boson_idx], boson_decay->second, ran)
	);
	}
	}
	// normalisation of momentum-conserving delta function
	weight_ = std::pow(2M_PI, 4);

	/** @TODO
	* uf (jet_param.min_pt) doesn't correspond to our final scale choice.
	* The HEJ scale generators currently expect a full event as input,
	* so fixing this is not completely trivial
	*/
	reconstruct_incoming(proc, subl_chance, subl_channels, pdf, E_beam, jet_param.min_pt, ran);
	if(status_ != Status::good) return;
	// set outgoing states
	begin_partons()->type = incoming_[0].type;
	rbegin_partons()->type = incoming_[1].type;

	maybe_turn_to_subl(subl_chance, subl_channels, proc, ran);

	if(proc.boson) couple_boson(*proc.boson, ran);
	}

	// pt generation, see eq:pt_sampling in developer manual
	double PhaseSpacePoint::gen_hard_pt(
	const int np , const double ptmin, const double ptmax, const double /* y */,
	HEJ::RNG & ran
	){
	// heuristic parameter for pt sampling, see eq:pt_par in developer manual
	const double ptpar = ptmin + np/5.;

	const double arctan = std::atan((ptmax - ptmin)/ptpar);
	const double xpt = ran.flat();
	const double pt = ptmin + ptparstd::tan(xptarctan);
	const double cosine = std::cos(xpt*arctan);
	weight_ = ptptpararctan/(cosinecosine);
	return pt;
	}

	double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, HEJ::RNG & ran) {
	constexpr double ptpar = 4.;

	const double r = ran.flat();
	const double pt = max_pt + ptpar/np*std::log(r);
	weight_ = ptptpar/(np*r);
	return pt;
	}

	double PhaseSpacePoint::gen_parton_pt(
	int count, JetParameters const & jet_param, double max_pt, double y,
	HEJ::RNG & ran
	) {
	constexpr double p_small_pt = 0.02;

	if(! jet_param.peak_pt) {
	return gen_hard_pt(count, jet_param.min_pt, max_pt, y, ran);
	}
	const double r = ran.flat();
	if(r > p_small_pt) {
	weight_ /= 1. - p_small_pt;
	return gen_hard_pt(count, *jet_param.peak_pt, max_pt, y, ran);
	}
	weight_ /= p_small_pt;
	const double pt = gen_soft_pt(count, *jet_param.peak_pt, ran);
	if(pt < jet_param.min_pt) {
	weight_=0.0;
	status_ = Status::not_enough_jets;
	return jet_param.min_pt;
	}
	return pt;
	}

	std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_LO_partons(
	int np, bool is_pure_jets,
	JetParameters const & jet_param,
	double max_pt,
	HEJ::RNG & ran
	){
	if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};

	weight_ /= std::pow(16.*std::pow(M_PI,3),np);
	weight_ /= std::tgamma(np+1); //remove rapidity ordering
	std::vector<fastjet::PseudoJet> partons;
	partons.reserve(np);
	for(int i = 0; i < np; ++i){
	const double y = -jet_param.max_y + 2jet_param.max_yran.flat();
	weight_ = 2jet_param.max_y;

	const bool is_last_parton = i+1 == np;
	if(is_pure_jets && is_last_parton) {
	constexpr double parton_mass_sq = 0.;
	partons.emplace_back(gen_last_momentum(partons, parton_mass_sq, y));
	break;
	}

	const double phi = 2M_PIran.flat();
	weight_ = 2.0M_PI;

	const double pt = gen_parton_pt(np, jet_param, max_pt, y, ran);
	if(weight_ == 0.0) return {};

	partons.emplace_back(fastjet::PtYPhiM(pt, y, phi));
	assert(jet_param.min_pt <= partons[i].pt());
	assert(partons[i].pt() <= max_pt+1e-5);
	}
	// Need to check that at LO, the number of jets = number of partons;
	fastjet::ClusterSequence cs(partons, jet_param.def);
	auto cluster_jets=cs.inclusive_jets(jet_param.min_pt);
	if (cluster_jets.size()!=unsigned(np)){
	weight_=0.0;
	status_ = Status::not_enough_jets;
	return {};
	}

	std::sort(begin(partons), end(partons), HEJ::rapidity_less{});

	return partons;
	}

	HEJ::Particle PhaseSpacePoint::gen_boson(
	HEJ::ParticleID bosonid, double mass, double width,
	HEJ::RNG & ran
	){
	// Usual phase space measure
	weight_ /= 16.*std::pow(M_PI, 3);

	// Generate a y Gaussian distributed around 0
	/// @TODO check magic numbers for different boson Higgs
	/// @TODO better sampling for W
	const double stddev_y = 1.6;
	const double y = random_normal(stddev_y, ran);
	const double r1 = ran.flat();
	const double s_boson = mass*(
	mass + widthstd::tan(M_PI/2.r1 + (r1-1.)*std::atan(mass/width))
	);
	// off-shell s_boson sampling, compensates for Breit-Wigner
	/// @TODO use a flag instead
	if(std::abs(bosonid) == HEJ::pid::Wp){
	weight_/=M_PIM_PI16.; //Corrects B-W factors, see git issue 132
	weight_= masswidth( M_PI+2.std::atan(mass/width) )
	/ ( 1. + std::cos( M_PIr1 + 2.(r1-1.)*std::atan(mass/width) ) );
	}

	auto p = gen_last_momentum(outgoing_, s_boson, y);

	return HEJ::Particle{bosonid, std::move(p), {}};
	}

	namespace {
	/// partons are ordered: even = anti, 0 = gluon
	HEJ::ParticleID index_to_pid(size_t i){
	if(!i) return HEJ::pid::gluon;
	return static_cast<HEJ::ParticleID>( i%2 ? (i+1)/2 : -i/2 );
	}

	/// partons are ordered: even = anti, 0 = gluon
	size_t pid_to_index(HEJ::ParticleID id){
	if(id==HEJ::pid::gluon) return 0;
	return id>0 ? id2-1 : std::abs(id)2;
	}

	using part_mask = std::bitset<11>; //!< Selection mask for partons

	part_mask init_allowed(HEJ::ParticleID const id){
	if(std::abs(id) == HEJ::pid::proton)
	return ~0;
	part_mask out{0};
	if(HEJ::is_parton(id))
	out[pid_to_index(id)] = 1;
	return out;
	}

	/// decides which "index" (see index_to_pid) are allowed for process
	part_mask allowed_quarks(HEJ::ParticleID const boson){
	if(std::abs(boson) != HEJ::pid::Wp){
	return ~1; // not a gluon
	}
	// special case W:
	// Wp: anti-down or up-type quark, no b/t
	// Wm: down or anti-up-type quark, no b/t
	return boson>0?0b00011001100
	:0b00100110010;
	}
	}

	std::array<part_mask,2> PhaseSpacePoint::incoming_AWZ(
	Process const & proc, Subleading const subl_channels,
	std::array<part_mask,2> allowed_partons,
	HEJ::RNG & ran
	){
	assert(proc.boson);
	auto couple_parton = allowed_quarks(*proc.boson);
	if( // coupling possible through input
	allowed_partons[0] == (couple_parton&allowed_partons[0])
	\|\| allowed_partons[1] == (couple_parton&allowed_partons[1])
	// cqqx possible
	\|\| (proc.njets >= 4 && subl_channels[subleading::cqqx])
	){
	return allowed_partons;
	}
	// eqqx only possible if one incoming is a gluon
	if(proc.njets >= 3 && subl_channels[subleading::eqqx]){
	couple_parton.set(pid_to_index(HEJ::ParticleID::gluon));
	}
	// only first can couple
	if( (allowed_partons[0]&couple_parton).any()
	&&(allowed_partons[1]&couple_parton).none()
	){
	return {allowed_partons[0]&couple_parton, allowed_partons[1]};
	}
	// only second can couple
	if( (allowed_partons[0]&couple_parton).none()
	&& (allowed_partons[1]&couple_parton).any()
	){
	return {allowed_partons[0], allowed_partons[1]&couple_parton};
	}
	// both can couple
	if( (allowed_partons[0]&couple_parton).any()
	&& (allowed_partons[1]&couple_parton).any()
	){
	double rnd = ran.flat();
	weight_*=3.;
	if(rnd<1./3.){
	return {
	allowed_partons[0] & couple_parton,
	allowed_partons[1] & ~couple_parton
	};
	}
	if(rnd<2./3.){
	return {
	allowed_partons[0] & ~couple_parton,
	allowed_partons[1] & couple_parton
	};
	}
	return {
	allowed_partons[0] & couple_parton,
	allowed_partons[1] & couple_parton
	};
	}
	throw std::invalid_argument{"Incoming state not allowed."};
	}

	std::array<part_mask,2> PhaseSpacePoint::incoming_eqqx(
	std::array<part_mask,2> allowed_partons, HEJ::RNG & ran
	){
	auto const gluon_idx = pid_to_index(HEJ::pid::gluon);
	auto & first_beam = allowed_partons[0];
	auto & second_beam = allowed_partons[1];
	if(first_beam[gluon_idx] && !second_beam[gluon_idx]){
	first_beam.reset();
	first_beam.set(gluon_idx);
	return allowed_partons;
	}
	if(!first_beam[gluon_idx] && second_beam[gluon_idx]) {
	second_beam.reset();
	second_beam.set(gluon_idx);
	return allowed_partons;
	}
	if(first_beam[gluon_idx] && second_beam[gluon_idx]) {
	// both beams can be gluons
	// if one can't be a quark everything is good
	auto first_quarks = first_beam;
	first_quarks.reset(gluon_idx);
	auto second_quarks = second_beam;
	second_quarks.reset(gluon_idx);
	if(first_quarks.none() \|\| second_quarks.none()){
	return allowed_partons;
	}
	// else choose one to be a gluon
	double rnd = ran.flat();
	weight_*=3.;
	if(rnd<1./3.){
	allowed_partons[0].reset();
	allowed_partons[0].set(gluon_idx);
	allowed_partons[1].reset(gluon_idx);
	} else if(rnd<2./3.){
	allowed_partons[1].reset();
	allowed_partons[1].set(gluon_idx);
	allowed_partons[0].reset(gluon_idx);
	} else {
	allowed_partons[0].reset();
	allowed_partons[0].set(gluon_idx);
	allowed_partons[1].reset();
	allowed_partons[1].set(gluon_idx);
	}
	return allowed_partons;
	}
	throw std::invalid_argument{
	"Incoming state not allowed for pure extremal qqx."};
	}

	std::array<part_mask,2> PhaseSpacePoint::incoming_uno(
	std::array<part_mask,2> allowed_partons, HEJ::RNG & ran
	){
	auto const gluon_idx = pid_to_index(HEJ::pid::gluon);
	auto & first_beam = allowed_partons[0];
	auto first_quarks = first_beam;
	first_quarks.reset(gluon_idx);
	auto & second_beam = allowed_partons[1];
	auto second_quarks = second_beam;
	second_quarks.reset(gluon_idx);
	if(first_quarks.any() && second_quarks.none()){
	first_beam.reset(gluon_idx);
	return allowed_partons;
	}
	if(first_quarks.none() && second_quarks.any()) {
	second_beam.reset(gluon_idx);
	return allowed_partons;
	}
	if(first_quarks.any() && second_quarks.any()) {
	// both beams can be quarks
	// if one can't be gluon everything is good
	if(!first_beam[gluon_idx] \|\| !second_beam[gluon_idx]){
	return allowed_partons;
	}
	// else choose one to be a quark
	double rnd = ran.flat();
	weight_*=3.;
	if(rnd<1./3.){
	allowed_partons[0].reset(gluon_idx);
	allowed_partons[1].reset();
	allowed_partons[1].set(gluon_idx);
	} else if(rnd<2./3.){
	allowed_partons[1].reset(gluon_idx);
	allowed_partons[0].reset();
	allowed_partons[0].set(gluon_idx);
	} else {
	allowed_partons[0].reset(gluon_idx);
	allowed_partons[1].reset(gluon_idx);
	}
	return allowed_partons;
	}
	throw std::invalid_argument{
	"Incoming state not allowed for pure unordered."};
	}

	/**
	* @brief Returns list of all allowed initial states partons
	*
	* checks which partons are allowed as initial state:
	* 1. only allow what is given in the Runcard (p -> all)
	* 2. A/W/Z require something to couple to
	* a) no qqx => no incoming gluon
	* b) 2j => no incoming gluon
	* c) >3j => can couple OR is gluon => 2 gluons become qqx later
	* 3. pure eqqx requires at least one gluon
	* 4. pure uno requires at least one quark
	*/
	std::array<part_mask,2> PhaseSpacePoint::allowed_incoming(
	Process const & proc,
	double const subl_chance, Subleading const subl_channels,
	HEJ::RNG & ran
	){
	// all possible incoming states
	std::array<part_mask,2> allowed_partons{
	init_allowed(proc.incoming[0]),
	init_allowed(proc.incoming[1])
	};
	// special case A/W/Z
	if(proc.boson && is_AWZ_proccess(proc)){
	allowed_partons = incoming_AWZ(proc, subl_channels, allowed_partons, ran);
	}
	// special case: pure subleading
	if(subl_chance!=1.){
	return allowed_partons;
	}
	auto other_channels = subl_channels;
	// pure eqqx
	other_channels.reset(subleading::eqqx);
	if(other_channels.none()){
	return incoming_eqqx(allowed_partons, ran);
	}
	other_channels = subl_channels;
	// pure uno
	other_channels.reset(subleading::uno);
	if(other_channels.none()){
	return incoming_uno(allowed_partons, ran);
	}
	return allowed_partons;
	}


	void PhaseSpacePoint::reconstruct_incoming(
	Process const & proc,
	double const subl_chance, Subleading const subl_channels,
	HEJ::PDF & pdf, double E_beam,
	double uf,
	HEJ::RNG & ran
	){
	std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
	// calculate xa, xb
	const double sqrts=2*E_beam;
	const double xa=(incoming_[0].E()-incoming_[0].pz())/sqrts;
	const double xb=(incoming_[1].E()+incoming_[1].pz())/sqrts;
	// abort if phase space point is outside of collider energy reach
	if (xa>1. \|\| xb>1.){
	weight_=0;
	status_ = Status::too_much_energy;
	return;
	}
	auto const & ids = proc.incoming;
	std::array<part_mask,2> allowed_partons
	= allowed_incoming(proc, subl_chance, subl_channels, ran);
	for(size_t i = 0; i < 2; ++i){
	if(ids[i] == HEJ::pid::proton \|\| ids[i] == HEJ::pid::p_bar){
	// pick ids according to pdfs
	incoming_[i].type =
	generate_incoming_id(i, i?xb:xa, uf, pdf, allowed_partons[i], ran);
	} else {
	assert(allowed_partons[i][pid_to_index(ids[i])]);
	incoming_[i].type = ids[i];
	}
	}
	assert(momentum_conserved(1e-7));
	}

	HEJ::ParticleID PhaseSpacePoint::generate_incoming_id(
	size_t const beam_idx, double const x, double const uf,
	HEJ::PDF & pdf, part_mask allowed_partons, HEJ::RNG & ran
	){
	- std::array<double,11> pdf_wt;
	+ std::array<double,11> pdf_wt{};
	pdf_wt[0] = allowed_partons[0]?
	std::fabs(pdf.pdfpt(beam_idx,x,uf,HEJ::pid::gluon)):0.;
	double pdftot = pdf_wt[0];
	for(size_t i = 1; i < pdf_wt.size(); ++i){
	pdf_wt[i] = allowed_partons[i]?
	4./9.*std::fabs(pdf.pdfpt(beam_idx,x,uf,index_to_pid(i))):0;
	pdftot += pdf_wt[i];
	}
	const double r1 = pdftot * ran.flat();
	double sum = 0;
	for(size_t i=0; i < pdf_wt.size(); ++i){
	if (r1 < (sum+=pdf_wt[i])){
	weight_*= pdftot/pdf_wt[i];
	return index_to_pid(i);
	}
	}
	std::cerr << "Error in choosing incoming parton: "<<x<<" "<<uf<<" "
	<<sum<<" "<<pdftot<<" "<<r1<<std::endl;
	throw std::logic_error{"Failed to choose parton flavour"};
	}

	void PhaseSpacePoint::couple_boson(
	HEJ::ParticleID const boson, HEJ::RNG & ran
	){
	if(std::abs(boson) != HEJ::pid::Wp) return; // only matters for W
	/// @TODO this could be use to sanity check gamma and Z

	// find all possible quarks
	std::vector<PartonIterator> allowed_parts;
	for(auto part_it=begin_partons(); part_it!=end_partons(); ++part_it){
	// Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom
	if ( can_couple_to_W(*part_it, boson) )
	allowed_parts.push_back(part_it);
	}
	if(allowed_parts.size() == 0){
	throw std::logic_error{"Found no parton for coupling with boson"};
	}

	// select one and flip it
	size_t idx = 0;
	if(allowed_parts.size() > 1){
	/// @TODO more efficient sampling
	/// old code: probability[i] = exp(parton[i].y - W.y)
	idx = std::floor(ran.flat()*allowed_parts.size());
	weight_ *= allowed_parts.size();
	}
	const int W_charge = boson>0?1:-1;
	allowed_parts[idx]->type =
	static_cast<HEJ::ParticleID>( allowed_parts[idx]->type - W_charge );
	}

	double PhaseSpacePoint::random_normal( double stddev, HEJ::RNG & ran ){
	const double r1 = ran.flat();
	const double r2 = ran.flat();
	const double lninvr1 = -std::log(r1);
	const double result = stddevstd::sqrt(2.lninvr1)std::cos(2.M_PI*r2);
	weight_ = exp(resultresult/(2stddevstddev))std::sqrt(2.M_PI)*stddev;
	return result;
	}

	bool PhaseSpacePoint::momentum_conserved(double ep) const{
	fastjet::PseudoJet diff;
	for(auto const & in: incoming()) diff += in.p;
	for(auto const & out: outgoing()) diff -= out.p;
	return HEJ::nearby_ep(diff, fastjet::PseudoJet{}, ep);
	}

	Decay PhaseSpacePoint::select_decay_channel(
	std::vector<Decay> const & decays,
	HEJ::RNG & ran
	){
	double br_total = 0.;
	for(auto const & decay: decays) br_total += decay.branching_ratio;
	// adjust weight
	// this is given by (channel branching ratio)/(chance to pick channel)
	// where (chance to pick channel) =
	// (channel branching ratio)/(total branching ratio)
	weight_ *= br_total;
	if(decays.size()==1) return decays.front();
	const double r1 = br_total*ran.flat();
	double br_sum = 0.;
	for(auto const & decay: decays){
	br_sum += decay.branching_ratio;
	if(r1 < br_sum) return decay;
	}
	throw std::logic_error{"unreachable"};
	}

	std::vector<HEJ::Particle> PhaseSpacePoint::decay_boson(
	HEJ::Particle const & parent,
	std::vector<Decay> const & decays,
	HEJ::RNG & ran
	){
	const auto channel = select_decay_channel(decays, ran);
	return decay_boson(parent, channel.products, ran);
	}

	std::vector<HEJ::Particle> PhaseSpacePoint::decay_boson(
	HEJ::Particle const & parent,
	std::vector<HEJ::ParticleID> const & decays,
	HEJ::RNG & ran
	){
	if(decays.size() != 2){
	throw HEJ::not_implemented{
	"only decays into two particles are implemented"
	};
	}
	std::vector<HEJ::Particle> decay_products(decays.size());
	for(size_t i = 0; i < decays.size(); ++i){
	decay_products[i].type = decays[i];
	}
	// choose polar and azimuth angle in parent rest frame
	const double E = parent.m()/2;
	const double theta = 2.M_PIran.flat();
	const double cos_phi = 2.*ran.flat()-1.; // Jacobian Factors for W in line 418
	const double sin_phi = std::sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
	const double px = Estd::cos(theta)sin_phi;
	const double py = Estd::sin(theta)sin_phi;
	const double pz = E*cos_phi;
	decay_products[0].p.reset(px, py, pz, E);
	decay_products[1].p.reset(-px, -py, -pz, E);
	for(auto & particle: decay_products) particle.p.boost(parent.p);
	return decay_products;
	}
	}
	diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc
	index 91f75cb..76be58f 100644
	--- a/FixedOrderGen/src/config.cc
	+++ b/FixedOrderGen/src/config.cc
	@@ -1,440 +1,442 @@
	/**
	* \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/exceptions.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/YAMLreader.hh"
	+#include "math.h"

	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", "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, end = 0; end != str.npos;){
	+ for(size_t begin = 0, 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 "+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 qqx" \|\| name == "cqqx")
	return channel.set(cqqx);
	if(name == "extremal qqx" \|\| name == "eqqx")
	return channel.set(eqqx);
	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;
	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:
	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;
	}

	HEJ::ParticleProperties get_particle_properties(
	YAML::Node const & node, std::string const & entry,
	std::string const & boson
	){
	- HEJ::ParticleProperties result;
	+ HEJ::ParticleProperties result{};
	set_from_yaml(result.mass, node, entry, boson, "mass");
	set_from_yaml(result.width, node, entry, boson, "width");
	return result;
	}

	HEJ::EWConstants get_ew_parameters(YAML::Node const & node){
	HEJ::EWConstants result;
	- double vev;
	+ double vev = NAN;
	set_from_yaml(vev, node, "vev");
	result.set_vevWZH(vev,
	get_particle_properties(node, "particle properties", "W"),
	get_particle_properties(node, "particle properties", "Z"),
	get_particle_properties(node, "particle properties", "Higgs")
	);
	return result;
	}

	UnweightSettings get_unweight(
	YAML::Node const & node, std::string const & entry
	){
	- UnweightSettings result;
	+ 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"]);

	config.ew_parameters = 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"]){
	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"]) 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;
	}
	}
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:21 PM (1 d, 16 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804937
Default Alt Text: (58 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions