No OneTemporary
Actions

Size

64 KB

Subscribers

None

View Options

	diff --git a/Changes-API.md b/Changes-API.md
	index e4bbd44..7afe670 100644
	--- a/Changes-API.md
	+++ b/Changes-API.md
	@@ -1,87 +1,90 @@
	# Changelog for HEJ API

	This log lists only changes on the HEJ API. These are primarily code changes
	relevant for calling HEJ as an API. This file should only be read as an addition
	to [`Changes.md`](Changes.md), where the main features are documented.

	## Version 2.1

	### 2.1.0

	#### Changes to Event class
	* Restructured `Event` class
	- `Event` can now only be build from a (new) `Event::EventData` class
	- Removed default constructor for `Event`
	- `Event::EventData` replaces the old `UnclusteredEvent` struct.
	- `UnclusteredEvent` is now deprecated, and will be removed in HEJ Version
	2.3.0
	- Removed `Event.unclustered()` function
	- Added new member function `Events.parameters()`, to directly access
	(underlying) `Parameters<EventParameters>`
	- New member functions `begin_partons`, `end_partons` (`rbegin_partons`,
	`rend_partons`) with aliases `cbegin_partons`, `cend_partons`
	(`crbegin_partons`, `crend_partons`) for constant (reversed) iterators over
	outgoing partons.
	* New function `Event::EventData.reconstruct_intermediate()` to reconstruct
	bosons from decays, e.g. `positron + nu_e => Wp`
	* Added optional Colour charges to particles (`Particle.colour`)
	- Colour connection in the HEJ limit can be generated via
	`Event.generate_colours` (automatically done in the resummation)
	- Colour configuration of input events can be checked with
	`Event.is_leading_colour`
	* Added function `Event.valid_hej_state` to check if event _could have_ been
	produced by `HEJ` according to the `max ext soft pt fraction` cut on the jets

	#### Updated functions
	* Renamed `EventType::nonHEJ` to `EventType::non_resummable` and `is_HEJ()`
	to `is_resummable()` such that run card is consistent with internal workings
	* Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public
	* New `EventReweighter` member function `treatment` to query the
	treatment with respect to resummation for the various event types.
	* Added auxiliary functions to obtain a `std::string` from and
	`EventDescription` (`to_string` for human readable, and `to_simple_string`
	for easy parsable string)

	#### New classes
	* New class `Unweighter` to do unweighting
	* New class `CrossSectionAccumulator` to keep track of Cross Section of the
	different subprocess
	* New template struct `Parameters` similar to old `Weights`
	- `Weights` are now an alias for `Parameters<double>`. Calling `Weights` did
	not change
	- `Weights.hh` was replaced by `Parameters.hh`. The old `Weights.hh` header
	will be removed in HEJ Version 2.2.0
	* Function to multiplication and division of `EventParameters.weight` by double
	- This can be combined with `Parameters`, e.g.
	`Parameters<EventParameters>*Weights`, see also `Events.parameters()`
	- Moved `EventParameters` to `Parameters.hh` header
	* new class `EWConstants` replaces previously hard coded `vev`
	- - `EWConstants` have to be set in the general `config` and the
	+ - `EWConstants` have to be set in the general `Config` and the
	`MatrixElementConfig`

	#### Input/Output
	* New abstract `EventReader` class, as base for reading events from files
	- Moved LHE file reader to `HEJ::LesHouchesReader`
	* Support reading (`HDF5Reader`) and writing (`HDF5Writer`) `hdf5` files
	* New `BufferedEventReader` class that allows to put events back into
	the reader.
	* `get_analysis` now requires `YAML::Node` and `LHEF::HEPRUP` as arguments
	* Replaced `HepMCInterface` and `HepMCWriter` by `HepMCInterfaceX` and
	`HepMCWriterX` respectively, with `X` being the major version of HepMC (2 or
	3)
	- Renamed `HepMCInterfaceX::init_kinematics` to `HepMCInterfaceX::init_event`
	and protected it, use `HepMCInterfaceX::operator()` instead
	- Removed redundant `HepMCInterfaceX::add_variation` function

	#### Linking
	* Export cmake target `HEJ::HEJ` to link directly against `libHEJ`

	+#### Miscellaneous
	+* Capitalisation of `Config.hh` now matches class `Config` (was `config.hh`)
	+
	## Version 2.0

	### 2.0.4

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

	### 2.0.0

	* First release
	diff --git a/FixedOrderGen/include/config.hh b/FixedOrderGen/include/config.hh
	index 45e4c16..6dfff7b 100644
	--- a/FixedOrderGen/include/config.hh
	+++ b/FixedOrderGen/include/config.hh
	@@ -1,45 +1,45 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include "yaml-cpp/yaml.h"

	#include "HEJ/HiggsCouplingSettings.hh"
	#include "HEJ/optional.hh"
	-#include "HEJ/config.hh"
	+#include "HEJ/Config.hh"
	#include "HEJ/output_formats.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/EWConstants.hh"

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

	namespace HEJFOG{

	struct Config{
	Process process;
	size_t events;
	JetParameters jets;
	Beam beam;
	int pdf_id;
	double subleading_fraction;
	unsigned int subleading_channels; //! < see HEJFOG::Subleading
	ParticlesDecayMap particle_decays;
	YAML::Node analysis_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 dd53fd4..c92661e 100644
	--- a/FixedOrderGen/src/EventGenerator.cc
	+++ b/FixedOrderGen/src/EventGenerator.cc
	@@ -1,93 +1,93 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "EventGenerator.hh"

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

	-#include "HEJ/config.hh"
	+#include "HEJ/Config.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/EWConstants.hh"

	namespace HEJFOG{
	EventGenerator::EventGenerator(
	Process process,
	Beam beam,
	HEJ::ScaleGenerator scale_gen,
	JetParameters jets,
	int pdf_id,
	double subl_change,
	unsigned int subl_channels,
	ParticlesDecayMap particle_decays,
	HEJ::HiggsCouplingSettings Higgs_coupling,
	HEJ::EWConstants ew_parameters,
	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)},
	subl_change_{subl_change},
	subl_channels_{subl_channels},
	particle_decays_{std::move(particle_decays)},
	ew_parameters_{ew_parameters},
	ran_{ran}
	{
	}

	HEJ::optional<HEJ::Event> EventGenerator::gen_event(){
	HEJFOG::PhaseSpacePoint psp{
	process_,
	jets_,
	pdf_, beam_.energy,
	subl_change_, subl_channels_,
	particle_decays_,
	ew_parameters_,
	ran_
	};
	status_ = psp.status();
	if(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(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/config.cc b/FixedOrderGen/src/config.cc
	index f068ff9..e241857 100644
	--- a/FixedOrderGen/src/config.cc
	+++ b/FixedOrderGen/src/config.cc
	@@ -1,427 +1,427 @@
	/**
	* \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/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", "vev",
	"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", "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;){
	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,
	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;
	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;
	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;
	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"]);

	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 "+HEJ::name(*config.process.boson)+" is required"};
	}
	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/include/HEJ/config.hh b/include/HEJ/Config.hh
	similarity index 100%
	rename from include/HEJ/config.hh
	rename to include/HEJ/Config.hh
	diff --git a/include/HEJ/EventReweighter.hh b/include/HEJ/EventReweighter.hh
	index 0e1af41..ed5a953 100644
	--- a/include/HEJ/EventReweighter.hh
	+++ b/include/HEJ/EventReweighter.hh
	@@ -1,201 +1,201 @@
	/** \file
	* \brief Declares the EventReweighter class
	*
	* EventReweighter is the main class used within HEJ 2. It reweights the
	* resummation events.
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <array>
	#include <functional>
	#include <utility>
	#include <vector>

	-#include "HEJ/config.hh"
	+#include "HEJ/Config.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/MatrixElement.hh"
	#include "HEJ/Parameters.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/RNG.hh"
	#include "HEJ/ScaleFunction.hh"
	#include "HEJ/StatusCode.hh"

	namespace LHEF {
	class HEPRUP;
	}

	namespace HEJ {
	class Event;

	//! Beam parameters
	/**
	* Currently, only symmetric beams are supported,
	* so there is a single beam energy.
	*/
	struct Beam{
	double E; /*< Beam energy /
	std::array<ParticleID, 2> type; /*< Beam particles /
	};

	//! Main class for reweighting events in HEJ.
	class EventReweighter{
	using EventType = event_type::EventType;

	public:
	EventReweighter(
	Beam beam, /*< Beam Energy /
	int pdf_id, /*< PDF ID /
	ScaleGenerator scale_gen, /*< Scale settings /
	EventReweighterConfig conf, /*< Configuration parameters /
	HEJ::RNG & ran /*< Random number generator /
	);

	EventReweighter(
	LHEF::HEPRUP const & heprup, /*< LHEF event header /
	ScaleGenerator scale_gen, /*< Scale settings /
	EventReweighterConfig conf, /*< Configuration parameters /
	HEJ::RNG & ran /*< Random number generator /
	);

	//! Get the used pdf
	PDF const & pdf() const;

	//! Get event treatment
	EventTreatment treatment(EventType type) const;

	//! Generate resummation events for a given fixed-order event
	/**
	* @param ev Fixed-order event corresponding
	* to the resummation events
	* @param num_events Number of trial resummation configurations.
	* @returns A vector of resummation events.
	*
	* The result vector depends on the type of the input event and the
	* treatment of different types as specified in the constructor:
	*
	* \ref reweight The result vector contains between
	* 0 and num_events resummation events.
	*
	* \ref keep If the input event passes the resummation jet cuts
	* the result vector contains one event. Otherwise it is empty.
	*
	* \ref discard The result vector is empty
	*/
	std::vector<Event> reweight(
	Event const & ev,
	int num_events
	);

	//! Gives all StatusCodes of the last reweight()
	/**
	* Each StatusCode corresponds to one tried generation. Only good
	* StatusCodes generated an event.
	*/
	std::vector<StatusCode> const & status() const {
	return status_;
	}

	private:
	template<typename... T>
	PDF const & pdf(T&& ...);

	/** \internal
	* \brief main generation/reweighting function:
	* generate phase space points and divide out Born factors
	*/
	std::vector<Event> gen_res_events(
	Event const & ev, size_t num_events
	);
	std::vector<Event> rescale(
	Event const & Born_ev, std::vector<Event> events
	) const;

	/** \internal
	* \brief Do the Jets pass the resummation Cuts?
	*
	* @param ev Event in Question
	* @returns 0 or 1 depending on if ev passes Jet Cuts
	*/
	bool jets_pass_resummation_cuts(Event const & ev) const;

	/** \internal
	* \brief pdf_factors Function
	*
	* @param ev Event in Question
	* @returns EventFactor due to PDFs
	*
	* Calculates the Central value and the variation due
	* to the PDF choice made.
	*/
	Weights pdf_factors(Event const & ev) const;

	/** \internal
	* \brief matrix_elements Function
	*
	* @param ev Event in question
	* @returns EventFactor due to MatrixElements
	*
	* Calculates the Central value and the variation due
	* to the Matrix Element.
	*/
	Weights matrix_elements(Event const & ev) const;

	/** \internal
	* \brief Scale-dependent part of fixed-order matrix element
	*
	* @param ev Event in question
	* @returns EventFactor scale variation due to FO-ME.
	*
	* This is only called to compute the scale variation for events where
	* we don't do resummation (e.g. non-FKL).
	* Since at tree level the scale dependence is just due to alpha_s,
	* it is enough to return the alpha_s(mur) factors in the matrix element.
	* The rest drops out in the ratio of (output event ME)/(input event ME),
	* so we never have to compute it.
	*/
	Weights fixed_order_scale_ME(Event const & ev) const;

	/** \internal
	* \brief Computes the tree level matrix element
	*
	* @param ev Event in Question
	* @returns HEJ approximation to Tree level Matrix Element
	*
	* This computes the HEJ approximation to the tree level FO
	* Matrix element which is used within the LO weighting process.
	*/
	double tree_matrix_element(Event const & ev) const;

	//! \internal General parameters
	EventReweighterConfig param_;

	//! \internal Beam energy
	double E_beam_;

	//! \internal PDF
	PDF pdf_;

	//! \internal Object to calculate the square of the matrix element
	MatrixElement MEt2_;
	//! \internal Object to calculate event renormalisation and factorisation scales
	ScaleGenerator scale_gen_;
	/** \internal random number generator
	*
	* \note We use a reference_wrapper so that EventReweighter objects can
	* still be copied (which would be impossible with a reference).
	*/
	std::reference_wrapper<HEJ::RNG> ran_;
	std::vector<StatusCode> status_;
	};

	template<typename... T>
	PDF const & EventReweighter::pdf(T&&... t){
	return pdf_ = PDF{std::forward<T>(t)...};
	}

	} // namespace HEJ
	diff --git a/include/HEJ/MatrixElement.hh b/include/HEJ/MatrixElement.hh
	index 0f25944..8991863 100644
	--- a/include/HEJ/MatrixElement.hh
	+++ b/include/HEJ/MatrixElement.hh
	@@ -1,186 +1,186 @@
	/** \file
	* \brief Contains the MatrixElement Class
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <functional>
	#include <vector>

	#include "fastjet/PseudoJet.hh"

	-#include "HEJ/config.hh"
	+#include "HEJ/Config.hh"
	#include "HEJ/Parameters.hh"
	#include "HEJ/PDG_codes.hh"

	namespace CLHEP {
	class HepLorentzVector;
	}

	namespace HEJ{
	class Event;
	class Particle;

	//! Class to calculate the squares of matrix elements
	class MatrixElement{
	public:
	/** \brief MatrixElement Constructor
	* @param alpha_s Function taking the renormalisation scale
	* and returning the strong coupling constant
	* @param conf General matrix element settings
	*/
	MatrixElement(
	std::function<double (double)> alpha_s,
	MatrixElementConfig conf
	);

	/**
	* \brief squares of regulated HEJ matrix elements
	* @param event The event for which to calculate matrix elements
	* @returns The squares of HEJ matrix elements including virtual corrections
	*
	* This function returns one value for the central parameter choice
	* and one additional value for each entry in \ref Event.variations().
	* See eq. (22) in \cite Andersen:2011hs for the definition of the squared
	* matrix element.
	*
	* \internal Relation to standard HEJ Met2: MatrixElement = Met2shat^2/(pdftapdftb)
	*/
	Weights operator()(Event const & event) const;

	//! Squares of HEJ tree-level matrix elements
	/**
	* @param event The event for which to calculate matrix elements
	* @returns The squares of HEJ matrix elements without virtual corrections
	*
	* cf. eq. (22) in \cite Andersen:2011hs
	*/
	Weights tree(Event const & event) const;

	/**
	* \brief Virtual corrections to matrix element squares
	* @param event The event for which to calculate matrix elements
	* @returns The virtual corrections to the squares of the matrix elements
	*
	* The all order virtual corrections to LL in the MRK limit is
	* given by replacing 1/t in the scattering amplitude according to the
	* lipatov ansatz.
	*
	* cf. second-to-last line of eq. (22) in \cite Andersen:2011hs
	* note that indices are off by one, i.e. out[0].p corresponds to p_1
	*/
	Weights virtual_corrections(Event const & event) const;

	/**
	* \brief Scale-dependent part of tree-level matrix element squares
	* @param event The event for which to calculate matrix elements
	* @returns The scale-dependent part of the squares of the
	* tree-level matrix elements
	*
	* The tree-level matrix elements factorises into a renormalisation-scale
	* dependent part, given by the strong coupling to some power, and a
	* scale-independent remainder. This function only returns the former parts
	* for the central scale choice and all \ref Event.variations().
	*
	* @see tree, tree_kin
	*/
	Weights tree_param(Event const & event) const;

	/**
	* \brief Kinematic part of tree-level matrix element squares
	* @param event The event for which to calculate matrix elements
	* @returns The kinematic part of the squares of the
	* tree-level matrix elements
	*
	* The tree-level matrix elements factorises into a renormalisation-scale
	* dependent part, given by the strong coupling to some power, and a
	* scale-independent remainder. This function only returns the latter part.
	* Since it does not depend on the parameter variations, only a single value
	* is returned.
	*
	* @see tree, tree_param
	*/
	double tree_kin(Event const & event) const;

	private:
	double tree_param(
	Event const & event,
	double mur
	) const;

	double virtual_corrections_W(
	Event const & event,
	double mur,
	Particle const & WBoson
	) const;
	double virtual_corrections(
	Event const & event,
	double mur
	) const;

	//! \internal cf. last line of eq. (22) in \cite Andersen:2011hs
	double omega0(
	double alpha_s, double mur,
	fastjet::PseudoJet const & q_j
	) const;

	double tree_kin_jets(
	Event const & ev
	) const;
	double tree_kin_W(
	Event const & ev
	) const;
	double tree_kin_Higgs(
	Event const & ev
	) const;
	double tree_kin_Higgs_first(
	Event const & ev
	) const;
	double tree_kin_Higgs_last(
	Event const & ev
	) const;

	/**
	* \internal
	* \brief Higgs inbetween extremal partons.
	*
	* Note that in the case of unordered emission, the Higgs is always
	* treated as if in between the extremal (FKL) partons, even if its
	* rapidity is outside the extremal parton rapidities
	*/
	double tree_kin_Higgs_between(
	Event const & ev
	) const;

	double tree_param_partons(
	double alpha_s, double mur,
	std::vector<Particle> const & partons
	) const;

	std::vector<int> in_extremal_jet_indices(
	std::vector<fastjet::PseudoJet> const & partons
	) const;

	std::vector<Particle> tag_extremal_jet_partons(
	Event const & ev
	) const;

	double MH2_forwardH(
	CLHEP::HepLorentzVector const & p1out,
	CLHEP::HepLorentzVector const & p1in,
	pid::ParticleID type2,
	CLHEP::HepLorentzVector const & p2out,
	CLHEP::HepLorentzVector const & p2in,
	CLHEP::HepLorentzVector const & pH,
	double t1, double t2
	) const;

	std::function<double (double)> alpha_s_;

	MatrixElementConfig param_;
	};

	} // namespace HEJ
	diff --git a/include/HEJ/PhaseSpacePoint.hh b/include/HEJ/PhaseSpacePoint.hh
	index f0cc3e1..19ad82e 100644
	--- a/include/HEJ/PhaseSpacePoint.hh
	+++ b/include/HEJ/PhaseSpacePoint.hh
	@@ -1,197 +1,197 @@
	/** \file
	* \brief Contains the PhaseSpacePoint Class
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <array>
	#include <functional>
	#include <unordered_map>
	#include <vector>

	-#include "HEJ/config.hh"
	+#include "HEJ/Config.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/RNG.hh"
	#include "HEJ/StatusCode.hh"

	namespace HEJ{
	class Event;

	//! A point in resummation phase space
	class PhaseSpacePoint{
	public:
	//! No default PhaseSpacePoint Constructor
	PhaseSpacePoint() = delete;

	//! PhaseSpacePoint Constructor
	/**
	* @param ev Clustered Jet Event
	* @param conf Configuration parameters
	* @param ran Random number generator
	*/
	PhaseSpacePoint(
	Event const & ev,
	PhaseSpacePointConfig conf,
	RNG & ran
	);

	//! Get phase space point weight
	double weight() const{
	return weight_;
	}

	//! Access incoming particles
	std::array<Particle, 2> const & incoming() const{
	return incoming_;
	}

	//! Access outgoing particles
	std::vector<Particle> const & outgoing() const{
	return outgoing_;
	}

	//! Particle decays
	/**
	* The key in the returned map corresponds to the index in the
	* vector returned by outgoing()
	*/
	std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
	return decays_;
	}

	//! Status code of generation
	StatusCode status() const{
	return status_;
	}

	static constexpr int ng_max = 1000; //!< maximum number of extra gluons

	private:
	//! /internal returns the clustered jets sorted in rapidity
	std::vector<fastjet::PseudoJet> cluster_jets(
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	bool pass_resummation_cuts(
	std::vector<fastjet::PseudoJet> const & jets
	) const;
	bool pass_extremal_cuts(
	fastjet::PseudoJet const & ext_parton,
	fastjet::PseudoJet const & jet
	) const;
	int sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets);
	int sample_ng_jets(int ng, std::vector<fastjet::PseudoJet> const & Born_jets);
	double probability_in_jet(
	std::vector<fastjet::PseudoJet> const & Born_jets
	) const;
	std::vector<fastjet::PseudoJet> gen_non_jet(
	int ng_non_jet,
	double ptmin, double ptmax
	);
	void rescale_qqx_rapidities(
	std::vector<fastjet::PseudoJet> & out_partons,
	std::vector<fastjet::PseudoJet> const & jets,
	const double ymin1, const double ymax2,
	const int qqxbackjet
	);
	void rescale_rapidities(
	std::vector<fastjet::PseudoJet> & partons,
	double ymin, double ymax
	);
	std::vector<fastjet::PseudoJet> reshuffle(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	fastjet::PseudoJet const & q
	);
	/** \interal final jet test
	* - number of jets must match Born kinematics
	* - no partons designated as nonjet may end up inside jets
	* - all other outgoing partons must end up inside jets
	* - the extremal (in rapidity) partons must be inside the extremal jets
	* - rapidities must be the same (by construction)
	*/
	bool jets_ok(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	void reconstruct_incoming(std::array<Particle, 2> const & Born_incoming);
	double phase_space_normalisation(
	int num_Born_jets,
	int num_res_partons
	) const;
	/** \interal Distribute gluon in jet
	* @param jets jets to distribute gluon in
	* @param ng_jets number of gluons
	* @param qqxbackjet position of first (backwards) qqx jet
	*
	* relies on JetSplitter
	*/
	std::vector<fastjet::PseudoJet> split(
	std::vector<fastjet::PseudoJet> const & jets,
	int ng_jets, size_t qqxbackjet
	);
	std::vector<int> distribute_jet_partons(
	int ng_jets, std::vector<fastjet::PseudoJet> const & jets
	);
	std::vector<fastjet::PseudoJet> split(
	std::vector<fastjet::PseudoJet> const & jets,
	std::vector<int> const & np_in_jet,
	size_t qqxbackjet
	);
	bool split_preserved_jets(
	std::vector<fastjet::PseudoJet> const & jets,
	std::vector<fastjet::PseudoJet> const & jet_partons
	) const;
	template<class Particle>
	Particle const & most_backward_FKL(
	std::vector<Particle> const & partons
	) const;
	template<class Particle>
	Particle const & most_forward_FKL(
	std::vector<Particle> const & partons
	) const;
	template<class Particle>
	Particle & most_backward_FKL(std::vector<Particle> & partons) const;
	template<class Particle>
	Particle & most_forward_FKL(std::vector<Particle> & partons) const;
	bool extremal_ok(
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	/** \internal
	* Assigns PDG IDs to outgoing partons, i.e. labels them as quarks
	*
	* \note This function assumes outgoing_ to be pure partonic when called,
	* i.e. A/W/Z/h bosons should _not be set_ at this stage
	*/
	void label_quarks(Event const & event);
	/** \internal
	* This function will label the qqx pair in a qqx event back to their
	* original types from the input event.
	*/
	void label_qqx(Event const & event);
	void copy_AWZH_boson_from(Event const & event);

	bool momentum_conserved() const;

	bool contains_idx(
	fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
	) const;

	bool unob_, unof_, qqxb_, qqxf_, qqxmid_;

	double weight_;

	PhaseSpacePointConfig param_;

	std::array<Particle, 2> incoming_;
	std::vector<Particle> outgoing_;
	//! \internal Particle decays in the format {outgoing index, decay products}
	std::unordered_map<size_t, std::vector<Particle>> decays_;

	std::reference_wrapper<HEJ::RNG> ran_;

	StatusCode status_;
	};

	} // namespace HEJ
	diff --git a/include/HEJ/YAMLreader.hh b/include/HEJ/YAMLreader.hh
	index 26927e9..993fbc7 100644
	--- a/include/HEJ/YAMLreader.hh
	+++ b/include/HEJ/YAMLreader.hh
	@@ -1,253 +1,253 @@
	/** \file
	* \brief The file which handles the configuration file parameters
	*
	* The configuration files parameters are read and then stored
	* within this objects.
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <string>
	#include <utility>
	#include <vector>

	#include "yaml-cpp/yaml.h"

	#include "fastjet/JetDefinition.hh"

	-#include "HEJ/config.hh"
	+#include "HEJ/Config.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/optional.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/utility.hh"

	namespace HEJ{
	class OutputFile;
	//! Load configuration from file
	/**
	* @param config_file Name of the YAML configuration file
	* @returns The HEJ 2 configuration
	*/
	Config load_config(std::string const & config_file);

	//! Set option using the corresponding YAML entry
	/**
	* @param setting Option variable to be set
	* @param yaml Root of the YAML configuration
	* @param names Name of the entry
	*
	* If the entry does not exist or has the wrong type or format
	* an exception is thrown.
	*
	* For example
	* @code
	* set_from_yaml(foobar, yaml, "foo", "bar")
	* @endcode
	* is equivalent to
	* @code
	* foobar = yaml["foo"]["bar"].as<decltype(foobar)>()
	* @endcode
	* with improved diagnostics on errors.
	*
	* @see set_from_yaml_if_defined
	*/
	template<typename T, typename... YamlNames>
	void set_from_yaml(
	T & setting,
	YAML::Node const & yaml, YamlNames const & ... names
	);

	//! Set option using the corresponding YAML entry, if present
	/**
	* @param setting Option variable to be set
	* @param yaml Root of the YAML configuration
	* @param names Name of the entry
	*
	* This function works similar to set_from_yaml, but does not
	* throw any exception if the requested YAML entry does not exist.
	*
	* @see set_from_yaml
	*/
	template<typename T, typename... YamlNames>
	void set_from_yaml_if_defined(
	T & setting,
	YAML::Node const & yaml, YamlNames const & ... names
	);

	//! Extract jet parameters from YAML configuration
	JetParameters get_jet_parameters(
	YAML::Node const & node, std::string const & entry
	);

	//! Extract Higgs coupling settings from YAML configuration
	HiggsCouplingSettings get_Higgs_coupling(
	YAML::Node const & node, std::string const & entry
	);

	//! Extract scale setting parameters from YAML configuration
	ScaleConfig to_ScaleConfig(YAML::Node const & yaml);

	//! Extract random number generator settings from YAML configuration
	RNGConfig to_RNGConfig(YAML::Node const & node, std::string const & entry);

	//! Check whether all options in configuration are supported
	/**
	* @param conf Configuration to be checked
	* @param supported Tree of supported options
	*
	* If conf contains an entry that does not appear in supported
	* an unknown_option exception is thrown. Sub-entries of "analysis"
	* (if present) are not checked.
	*
	* @see unknown_option
	*/
	void assert_all_options_known(
	YAML::Node const & conf, YAML::Node const & supported
	);

	namespace detail{
	void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml);
	void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml);
	void set_from_yaml(ParticleID & setting, YAML::Node const & yaml);
	void set_from_yaml(OutputFile & setting, YAML::Node const & yaml);
	void set_from_yaml(WeightType & setting, YAML::Node const & yaml);

	inline
	void set_from_yaml(YAML::Node & setting, YAML::Node const & yaml){
	setting = yaml;
	}

	template<typename Scalar>
	void set_from_yaml(Scalar & setting, YAML::Node const & yaml){
	assert(yaml);
	if(!yaml.IsScalar()){
	throw invalid_type{"value is not a scalar"};
	}
	try{
	setting = yaml.as<Scalar>();
	}
	catch(...){
	throw invalid_type{
	"value " + yaml.as<std::string>()
	+ " cannot be converted to a " + type_string(setting)
	};
	}
	}

	template<typename T>
	void set_from_yaml(optional<T> & setting, YAML::Node const & yaml){
	T tmp;
	set_from_yaml(tmp, yaml);
	setting = tmp;
	}

	template<typename T>
	void set_from_yaml(std::vector<T> & setting, YAML::Node const & yaml){
	assert(yaml);
	// special case: treat a single value like a vector with one element
	if(yaml.IsScalar()){
	setting.resize(1);
	return set_from_yaml(setting.front(), yaml);
	}
	if(yaml.IsSequence()){
	setting.resize(yaml.size());
	for(size_t i = 0; i < setting.size(); ++i){
	set_from_yaml(setting[i], yaml[i]);
	}
	return;
	}
	throw invalid_type{""};
	}

	template<typename T, typename FirstName, typename... YamlNames>
	void set_from_yaml(
	T & setting,
	YAML::Node const & yaml, FirstName const & name,
	YamlNames && ... names
	){
	if(!yaml[name]) throw missing_option{""};
	set_from_yaml(
	setting,
	yaml[name], std::forward<YamlNames>(names)...
	);
	}

	template<typename T>
	void set_from_yaml_if_defined(T & setting, YAML::Node const & yaml){
	return set_from_yaml(setting, yaml);
	}

	template<typename T, typename FirstName, typename... YamlNames>
	void set_from_yaml_if_defined(
	T & setting,
	YAML::Node const & yaml, FirstName const & name,
	YamlNames && ... names
	){
	if(!yaml[name]) return;
	set_from_yaml_if_defined(
	setting,
	yaml[name], std::forward<YamlNames>(names)...
	);
	}
	} // namespace detail

	template<typename T, typename... YamlNames>
	void set_from_yaml(
	T & setting,
	YAML::Node const & yaml, YamlNames const & ... names
	){
	try{
	detail::set_from_yaml(setting, yaml, names...);
	}
	catch(invalid_type const & ex){
	throw invalid_type{
	"In option " + join(": ", names...) + ": " + ex.what()
	};
	}
	catch(missing_option const &){
	throw missing_option{
	"No entry for mandatory option " + join(": ", names...)
	};
	}
	catch(std::invalid_argument const & ex){
	throw missing_option{
	"In option " + join(": ", names...) + ":"
	" invalid value " + ex.what()
	};
	}
	}

	template<typename T, typename... YamlNames>
	void set_from_yaml_if_defined(
	T & setting,
	YAML::Node const & yaml, YamlNames const & ... names
	){
	try{
	detail::set_from_yaml_if_defined(setting, yaml, names...);
	}
	catch(invalid_type const & ex){
	throw invalid_type{
	"In option " + join(": ", names...) + ": " + ex.what()
	};
	}
	catch(std::invalid_argument const & ex){
	throw missing_option{
	"In option " + join(": ", names...) + ":"
	" invalid value " + ex.what()
	};
	}
	}

	} // namespace HEJ

	namespace YAML {

	template<>
	struct convert<HEJ::OutputFile> {
	static Node encode(HEJ::OutputFile const & outfile);
	static bool decode(Node const & node, HEJ::OutputFile & out);
	};
	}
	diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
	index fe28c5d..d2b8bb6 100644
	--- a/src/bin/HEJ.cc
	+++ b/src/bin/HEJ.cc
	@@ -1,369 +1,369 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <chrono>
	#include <iostream>
	#include <limits>
	#include <memory>
	#include <numeric>

	#include "yaml-cpp/yaml.h"

	#include "fastjet/ClusterSequence.hh"

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

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

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

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

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

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

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

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

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

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


	void train(
	HEJ::Unweighter & unweighter,
	HEJ::BufferedEventReader & reader,
	HEJ::EventReweighter & reweighter,
	const int total_trials,
	const double max_dev,
	double reweight_factor,
	HEJ::JetParameters const & fixed_order_jets
	) {
	std::cout << "Reading up to " << total_trials << " training events...\n";
	auto FO_events = peek_events(reader, total_trials);
	const int trials = total_trials/FO_events.size();
	// adjust reweight factor so that the overall normalisation
	// is the same as in the full run
	reweight_factor *= trials;
	for(auto & hepeup: FO_events) {
	hepeup.setWeight(0, reweight_factor * hepeup.weight());
	}
	std::cout << "Training unweighter with "
	<< trials << '*' << FO_events.size() << " events\n";

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	// analysis
	for(auto & ev: resummed_events){
	//TODO: move pass_cuts to after phase space point generation
	if(analysis->pass_cuts(ev, FO_event)){
	analysis->fill(ev, FO_event);
	writer.write(ev);
	} else {
	ev.parameters()*=0; // do not use discarded events afterwards
	}
	}
	xs.fill_correlated(resummed_events);
	total_resum += resummed_events.size();
	if(progress) progress->increment(FO_event.central().weight);
	} // main event loop
	std::cout << '\n';
	analysis->finalise();

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

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

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

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

	return EXIT_SUCCESS;
	}
	diff --git a/t/test_psp.cc b/t/test_psp.cc
	index 1c502aa..c62e0c4 100644
	--- a/t/test_psp.cc
	+++ b/t/test_psp.cc
	@@ -1,68 +1,68 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "LHEF/LHEF.h"
	#include "HEJ/stream.hh"
	-#include "HEJ/config.hh"
	+#include "HEJ/Config.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/PhaseSpacePoint.hh"
	#include "HEJ/Ranlux64.hh"

	#include "hej_test.hh"

	namespace{
	constexpr int max_trials = 100;
	constexpr double max_ext_soft_pt_fraction = 0.1;
	const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
	constexpr double min_jet_pt = 50;
	}

	int main(int argn, char** argv) {
	if(argn != 2){
	std::cerr << "Usage: " << argv[0] << " eventfile";
	return EXIT_FAILURE;
	}

	HEJ::istream in{argv[1]};
	LHEF::Reader reader{in};
	LHEF::Writer writer{std::cerr};
	writer.heprup = reader.heprup;

	HEJ::PhaseSpacePointConfig conf;
	conf.jet_param = HEJ::JetParameters{jet_def, min_jet_pt};
	conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;

	HEJ::Ranlux64 ran{};

	while(reader.readEvent()){
	HEJ::Event::EventData ev_data{ reader.hepeup };
	shuffle_particles(ev_data);
	const HEJ::Event ev{ ev_data( jet_def, min_jet_pt ) };
	for(int trial = 0; trial < max_trials; ++trial){
	HEJ::PhaseSpacePoint psp{ev, conf, ran};
	if(psp.weight() != 0){
	HEJ::Event::EventData tmp_ev;
	tmp_ev.incoming = psp.incoming();
	tmp_ev.outgoing = psp.outgoing();
	tmp_ev.parameters.central = {0,0,0};
	shuffle_particles(tmp_ev);
	HEJ::Event out_ev{ tmp_ev(jet_def, min_jet_pt) };
	if(out_ev.type() != ev.type()){
	using HEJ::event_type::name;
	std::cerr << "Wrong class of phase space point:\n"
	"original event of class " << name(ev.type()) << ":\n";
	writer.hepeup = reader.hepeup;
	writer.writeEvent();
	std::cerr << "Phase space point of class " << name(out_ev.type()) << ":\n";
	writer.hepeup = to_HEPEUP(out_ev, &writer.heprup);
	writer.writeEvent();
	return EXIT_FAILURE;
	}
	}
	}
	}
	return EXIT_SUCCESS;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 3:50 PM (1 d, 20 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023347
Default Alt Text: (64 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions