No OneTemporary
Actions

Size

147 KB

Subscribers

None

View Options

	diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
	index 60cce8e..e7e4b10 100644
	--- a/FixedOrderGen/src/EventGenerator.cc
	+++ b/FixedOrderGen/src/EventGenerator.cc
	@@ -1,79 +1,81 @@
	#include "EventGenerator.hh"

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

	#include "RHEJ/Event.hh"
	#include "RHEJ/config.hh"

	namespace HEJFOG{
	EventGenerator::EventGenerator(
	Process process,
	Beam beam,
	RHEJ::ScaleGenerator & scale_gen,
	JetParameters jets,
	int pdf_id,
	double uno_chance,
	HiggsProperties higgs_properties,
	RHEJ::HiggsCouplingSettings Higgs_coupling,
	CLHEP::Ranlux64Engine & ran
	):
	pdf_{pdf_id, beam.particles[0], beam.particles[1]},
	ME_{
	[this](double mu){ return pdf_.Halphas(mu); },
	- jets.def, jets.min_pt,
	- false,
	- std::move(Higgs_coupling)
	+ RHEJ::MatrixElementConfig{
	+ RHEJ::JetParameters{jets.def, jets.min_pt},
	+ false,
	+ std::move(Higgs_coupling)
	+ }
	},
	scale_gen_{scale_gen},
	process_{std::move(process)},
	jets_{std::move(jets)},
	beam_{std::move(beam)},
	uno_chance_{uno_chance},
	higgs_properties_{std::move(higgs_properties)},
	ran_{ran}
	{
	}

	RHEJ::Event EventGenerator::gen_event(){
	HEJFOG::PhaseSpacePoint psp{
	process_,
	jets_.def, jets_.min_pt, jets_.max_y,
	pdf_, beam_.energy,
	uno_chance_,
	higgs_properties_,
	ran_
	};
	status_ = psp.status();
	if(status_ != good) return {};

	RHEJ::Event ev = scale_gen_(
	RHEJ::Event{
	to_UnclusteredEvent(std::move(psp)),
	jets_.def, jets_.min_pt
	}
	);

	const double shat = RHEJ::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 on this point
	ev.central().weight *= ME_.tree(
	ev.central().mur, ev.incoming(), ev.outgoing(), false
	)/(shat*shat);
	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(auto & var: ev.variations()){
	var.weight *= ME_.tree(
	var.mur, ev.incoming(), ev.outgoing(), false
	)/(shat*shat);
	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 cb0dd26..22a6d7b 100644
	--- a/FixedOrderGen/src/config.cc
	+++ b/FixedOrderGen/src/config.cc
	@@ -1,282 +1,283 @@
	#include "config.hh"

	#include <cctype>

	#include "RHEJ/config.hh"
	+#include "RHEJ/YAMLreader.hh"

	namespace HEJFOG{
	using RHEJ::set_from_yaml;
	using RHEJ::set_from_yaml_if_defined;

	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", "trials", "unordered fraction", "scales", "scale factors",
	"max scale ratio", "pdf", "event output", "analysis", "RanLux init",
	};
	// add subnodes to "supported" - the assigned value is irrelevant
	for(auto && opt: opts) supported[opt] = "";
	for(auto && jet_opt: {"min pt", "algorithm", "R", "max rapidity"}){
	supported["jets"][jet_opt] = "";
	}
	for(auto && Higgs_opt: {"mass", "width"}){
	supported["Higgs properties"][Higgs_opt] = "";
	}
	supported["Higgs properties"]["decays"]["into"] = "";
	supported["Higgs properties"]["decays"]["branching ratio"] = "";
	for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
	supported["Higgs coupling"][opt] = "";
	}
	for(auto && beam_opt: {"energy", "particles"}){
	supported["beam"][beam_opt] = "";
	}
	for(auto && unweight_opt: {"sample size", "max deviation"}){
	supported["unweight"][unweight_opt] = "";
	}
	return supported;
	}();
	return supported;
	}

	JetParameters get_jet_parameters(
	YAML::Node const & node, std::string const & entry
	){
	const auto p = RHEJ::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");
	return result;
	}

	Beam get_Beam(
	YAML::Node const & node, std::string const & entry
	){
	Beam beam;
	std::vector<RHEJ::ParticleID> particles;
	set_from_yaml(beam.energy, node, entry, "energy");
	set_from_yaml_if_defined(particles, node, entry, "particles");
	if(! particles.empty()){
	for(RHEJ::ParticleID particle: particles){
	if(particle != RHEJ::pid::p && particle != RHEJ::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', 'W+', 'W-', 'Z', 'H'"
	};
	}

	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() = RHEJ::to_ParticleID(particles[0]);
	result.incoming.back() = RHEJ::to_ParticleID(particles[1]);

	for(size_t i = 0; i < result.incoming.size(); ++i){
	const RHEJ::ParticleID in = result.incoming[i];
	if(
	in != RHEJ::pid::proton && in != RHEJ::pid::p_bar
	&& !RHEJ::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())){
	if(particles[i].back() != 'j'){
	throw invalid_outgoing(particles[i]);
	}
	result.njets += std::stoi(particles[i]);
	}
	else{
	const auto pid = RHEJ::to_ParticleID(particles[i]);
	if(!RHEJ::is_AWZH_boson(pid)){
	throw invalid_outgoing(particles[i]);
	}
	if(result.boson){
	throw std::invalid_argument{
	"More than one outgoing boson is not supported"
	};
	}
	result.boson = pid;
	}
	}
	if(result.njets < 2){
	throw std::invalid_argument{
	"Process has to include at least two jets ('j')"
	};
	}
	return result;
	}

	Decay get_decay(YAML::Node const & node){
	Decay decay;
	set_from_yaml(decay.products, node, "into");
	set_from_yaml(decay.branching_ratio, node, "branching ratio");
	return decay;
	}

	std::vector<Decay> get_decays(YAML::Node const & node){
	using YAML::NodeType;
	if(!node) return {};
	switch(node.Type()){
	case NodeType::Null:
	case NodeType::Undefined:
	return {};
	case NodeType::Scalar:
	throw RHEJ::invalid_type{"value is not a list of decays"};
	case NodeType::Map:
	return {get_decay(node)};
	case NodeType::Sequence:
	std::vector<Decay> result;
	for(auto && decay_str: node){
	result.emplace_back();
	set_from_yaml(result.back().products, decay_str, "into");
	set_from_yaml(result.back().branching_ratio, decay_str, "branching ratio");
	}
	return result;
	}
	throw std::logic_error{"unreachable"};
	}

	HiggsProperties get_higgs_properties(
	YAML::Node const & node, std::string const & entry
	){
	HiggsProperties result;
	set_from_yaml(result.mass, node, entry, "mass");
	set_from_yaml(result.width, node, entry, "width");
	try{
	result.decays = get_decays(node[entry]["decays"]);
	}
	catch(RHEJ::missing_option const & ex){
	throw RHEJ::missing_option{entry + ": decays: " + ex.what()};
	}
	catch(RHEJ::invalid_type const & ex){
	throw RHEJ::invalid_type{entry + ": decays: " + ex.what()};
	}
	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{
	RHEJ::assert_all_options_known(yaml, get_supported_options());
	}
	catch(RHEJ::unknown_option const & ex){
	throw RHEJ::unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
	}

	Config config;
	config.process = get_process(yaml, "process");
	set_from_yaml(config.trials, yaml, "trials");
	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 RHEJ::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.unordered_fraction, yaml, "unordered fraction");
	if(config.unordered_fraction < 0 \|\| config.unordered_fraction > 1){
	throw std::invalid_argument{
	"unordered fraction has to be between 0 and 1"
	};
	}
	config.Higgs_properties = get_higgs_properties(yaml, "Higgs properties");
	set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
	config.scale_gen = RHEJ::ScaleGenerator{RHEJ::to_ScaleConfig(yaml)};
	set_from_yaml_if_defined(config.output, yaml, "event output");
	set_from_yaml_if_defined(config.RanLux_init, yaml, "RanLux init");
	config.Higgs_coupling = RHEJ::get_Higgs_coupling(yaml, "Higgs coupling");
	if(yaml["unweight"]) config.unweight = get_unweight(yaml, "unweight");
	return config;
	}

	}

	Config load_config(std::string const & config_file){
	try{
	return to_Config(YAML::LoadFile(config_file));
	}
	catch(...){
	std::cerr << "Error reading " << config_file << ":\n ";
	throw;
	}
	}
	}
	diff --git a/include/RHEJ/CombinedEventWriter.hh b/include/RHEJ/CombinedEventWriter.hh
	index 44308a0..a1f1646 100644
	--- a/include/RHEJ/CombinedEventWriter.hh
	+++ b/include/RHEJ/CombinedEventWriter.hh
	@@ -1,53 +1,52 @@
	/** \file CombinedEventWriter.hh
	* \brief Details the CombinedEventWriter class
	*
	* CombinedEventWriter Class which handles all of the EventWriters being used by a single use
	* of RHEJ. It calls all other EventWriters in its list to use their write function and output.
	*/

	#pragma once

	#include <memory>
	#include <vector>

	#include "RHEJ/EventWriter.hh"
	#include "RHEJ/make_writer.hh"

	namespace LHEF{
	/** \struct HEPRUP CombinedEventWriter.hh "include/RHEJ/CombinedEventWriter.hh"
	* \brief HEPRUP struct which contains input event information.
	*
	* This struct contains important information from the LHEF header such as the beam energy.
	*/
	struct HEPRUP;
	}

	namespace RHEJ{
	- struct config;

	/** \class CombinedEventWriter CombinedEventWriter.hh "include/RHEJ/CombinedEventWriter.hh"
	* \brief Intended as an EventWriter to Multiple Output types.
	*
	* Inherits from EventWriter, and then uses its write function which calls all of the event
	* writers write function to output to their specific output files in their specific format.
	* This handles this complexity of having multiple EventWriters.
	*/
	class CombinedEventWriter: public EventWriter{

	public:
	//! CombinedEventWriter Constructor
	CombinedEventWriter(
	- /** Contains all configurations for running rHEJ */
	- Config const & conf,
	+ /** Files to which event output should be written */
	+ std::vector<OutputFile> const & outfiles,
	/** Contains important information from the Input */
	LHEF::HEPRUP const & heprup
	);

	void write(Event const &) override;

	private:
	//! A vector of the EventWriters
	std::vector<std::unique_ptr<EventWriter>> writers_;
	};

	}
	diff --git a/include/RHEJ/EventReweighter.hh b/include/RHEJ/EventReweighter.hh
	index 22229f8..f94e745 100644
	--- a/include/RHEJ/EventReweighter.hh
	+++ b/include/RHEJ/EventReweighter.hh
	@@ -1,175 +1,175 @@
	/** \file EventReweighter.hh
	* \brief Main class for reweighting events according to the reversed HEJ method
	*
	* EventReweighter Class is the main class used within RHEJ. It reweights the
	* resummation events.
	*/
	#pragma once

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

	#include "LHEF/LHEF.h"

	#include "RHEJ/config.hh"
	#include "RHEJ/PDF.hh"
	#include "RHEJ/utility.hh"
	#include "RHEJ/Event.hh"
	#include "RHEJ/event_types.hh"
	#include "RHEJ/PhaseSpacePoint.hh"
	#include "RHEJ/MatrixElement.hh"

	namespace RHEJ{

	/** \struct Beam EventReweighter.hh "include/RHEJ/EventReweighter.hh"
	* \brief Beam Struct. Contains beam energy and incoming particle content.
	*
	* Contains the value of the energy of the beam and the identity of the two incoming Particles
	*/
	struct Beam{
	double E;
	std::array<ParticleID, 2> type;
	};

	/** \class EventReweighter EventReweighter.hh "include/RHEJ/EventReweighter.hh"
	* \brief EventReweighter Main class for reweighting events in RHEJ.
	*
	* This is the class which reweights the resummation phase space points once they have been
	* generated by PhaseSpacePoint. This class also handles how many phase space points are generated
	* for every fixed order event which is processed.
	*/
	class EventReweighter{
	using EventType = event_type::EventType;
	public:

	EventReweighter(
	Beam beam, /*< Beam Energy /
	int pdf_id, /*< PDF ID /
	- Config const & conf /*< Configuration parameters /
	+ EventReweighterConfig conf /*< Configuration parameters /
	);

	EventReweighter(
	LHEF::HEPRUP const & heprup, /*< LHEF event header /
	- Config const & conf /*< Configuration parameters /
	+ EventReweighterConfig conf /*< Configuration parameters /
	);

	PDF const & pdf() const; /*< PDF Used /


	std::vector<Event> reweight(
	Event const & ev, /*< Event For Reweighing /
	int num_events /*< Number of Events /
	);


	std::vector<Event> reweight(
	Event const & ev, /*< Event For Reweighing /
	int num_events, /*< Number of Events /
	ScaleGenerator const & gen_scales /*< Generated Scales see config.hh /
	);

	private:

	/** \struct EventFactors EventReweighter.hh "include/RHEJ/EventReweighter.hh"
	* \brief EventFactors Struct containing central value and variations
	*
	* Contains the Central Value and Variation around
	*/
	struct EventFactors{
	double central; /*< Central Value for the Factor /
	std::vector<double> variations; /*< Vector of Variations from central value /
	};

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

	std::vector<Event> gen_res_events(
	Event const & ev, int num_events
	);
	std::vector<Event> rescale(
	Event const & Born_ev, std::vector<Event> events
	) const;

	/**
	* \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;

	/**
	* \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.
	*/
	EventFactors pdf_factors(Event const & ev) const;

	/**
	* \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.
	*/
	EventFactors matrix_elements(Event const & ev) const;

	/**
	* \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.
	*/
	EventFactors fixed_order_scale_ME(Event const & ev) const;

	/**
	* \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;

	/**
	* \brief Generation of Resummation Phase Space Points
	*
	* @param ev Event In Question
	* @returns Resummation PhaseSpacePoint
	*
	* Where is this even called? I can't find it with a grep -r in top directory???
	*/
	PhaseSpacePoint gen_PhaseSpacePoint(Event const & ev) const;


	- EventReweighterParameters param_;
	+ EventReweighterConfig param_;

	double E_beam_; /*< Energy of Beam /

	PDF pdf_; /*< Relevant PDF /

	MatrixElement MEt2_; /*< Matrix Element Object /

	};

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

	}
	diff --git a/include/RHEJ/HepMCWriter.hh b/include/RHEJ/HepMCWriter.hh
	index c6d2f0d..ce61c00 100644
	--- a/include/RHEJ/HepMCWriter.hh
	+++ b/include/RHEJ/HepMCWriter.hh
	@@ -1,47 +1,46 @@
	/** \file HepMCWriter.hh
	* \brief Contains the EventWriter which is necessary for HepMC Output.
	*/

	#pragma once

	#include <string>

	#include "RHEJ/EventWriter.hh"

	#include "LHEF/LHEF.h"

	namespace RHEJ{
	class Event;
	- struct Config;

	/** \class HepMCWriter HepMCWriter.hh "include/RHEJ/HepMCWriter.hh"
	* \brief This is an event writer specifically for HepMC output.
	* as such it inherits everything from the EventWriter class
	* and also includes HepMCWriter constructors
	*
	* Implementation note:
	* This uses the pimpl ("pointer to implementation") idiom.
	* HepMC support is optional and the implementation depends on the
	* HepMC version. Without pimpl, we would have to specify the HepMC version
	* via the preprocessor whenever this header is included. We don't want to
	* burden users of the rHEJ library (for example the HEJ fixed-order generator)
	* with those details
	*
	*/
	class HepMCWriter: public EventWriter{
	public:
	- HepMCWriter(std::string const & file, Config const & conf, LHEF::HEPRUP heprup);
	+ HepMCWriter(std::string const & file, LHEF::HEPRUP heprup);
	~HepMCWriter() override = default;

	void write(Event const & ev) override;

	private:
	struct HepMCWriterImpl;
	struct HepMCWriterImplDeleter {
	void operator()(HepMCWriterImpl* p);
	};

	std::unique_ptr<HepMCWriterImpl, HepMCWriterImplDeleter> impl_;
	};

	}
	diff --git a/include/RHEJ/LesHouchesWriter.hh b/include/RHEJ/LesHouchesWriter.hh
	index f8666f7..a312644 100644
	--- a/include/RHEJ/LesHouchesWriter.hh
	+++ b/include/RHEJ/LesHouchesWriter.hh
	@@ -1,57 +1,55 @@
	/** \file LesHouchesWriter.hh
	* \brief Contains the writer for LesHouches Output
	*/

	#pragma once

	#include <fstream>

	#include "RHEJ/EventWriter.hh"

	#include "LHEF/LHEF.h"

	namespace RHEJ{
	class Event;
	- struct Config;

	/** \class LesHouchesWriter LesHouchesWriter.hh "include/RHEJ/LesHouchesWriter.hh"
	* \brief EventWriter Class specifically for LesHouches Output
	*
	- * This is an event writer specifically for Les Houches Output. This
	+ * This is an event writer specifically for Les Houches Output. This
	* inherits from the EventWriter Class and contains its own
	- * contructors. TODO: in principle, this class should be movable but that
	+ * contructors. TODO: in principle, this class should be movable but that
	* somehow(?) breaks the write member function.
	**/
	class LesHouchesWriter : public EventWriter{
	public:
	- LesHouchesWriter(std::string const & file, Config const & conf, LHEF::HEPRUP heprup);
	+ LesHouchesWriter(std::string const & file, LHEF::HEPRUP heprup);
	LesHouchesWriter(LesHouchesWriter const & other) = delete;
	LesHouchesWriter & operator=(LesHouchesWriter const & other) = delete;
	// TODO: in principle, this class should be movable
	// but that somehow(?) breaks the write member function
	LesHouchesWriter(LesHouchesWriter && other) = delete;
	LesHouchesWriter & operator=(LesHouchesWriter && other) = delete;
	~LesHouchesWriter() override;

	void write(Event const & ev) override;

	private:
	void write_init(){
	writer_->init();
	}

	void rewrite_init();

	LHEF::HEPRUP & heprup(){
	return writer_->heprup;
	}

	LHEF::HEPEUP & hepeup(){
	return writer_->hepeup;
	}

	- Config const & config_;
	std::fstream out_;
	std::unique_ptr<LHEF::Writer> writer_;
	};
	}
	diff --git a/include/RHEJ/MatrixElement.hh b/include/RHEJ/MatrixElement.hh
	index ced2cb4..20302f4 100644
	--- a/include/RHEJ/MatrixElement.hh
	+++ b/include/RHEJ/MatrixElement.hh
	@@ -1,216 +1,216 @@
	/** \file MatrixElement.hh
	* \brief The header file which contains the MatrixElement Class
	*
	* This contains the MatrixElement Class which contains many functions
	* used to calculate many different MatrixElements and their components.
	*/
	#pragma once

	#include <functional>

	#include "RHEJ/config.hh"
	#include "RHEJ/utility.hh"
	#include "RHEJ/HiggsCouplingSettings.hh"

	#include "CLHEP/Vector/LorentzVector.h"


	namespace RHEJ{

	/** \class MatrixElement MatrixElement.hh "include/RHEJ/MatrixElement.hh
	* \brief MatrixElement class. Functions for obtaining various ME and components.
	*/
	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,
	- Config const & conf
	+ MatrixElementConfig conf
	);

	/**
	* \brief regulated HEJ matrix element
	* @param mur Value of the renormalisation scale
	* @param incoming Incoming particles
	* @param partons Outgoing particles
	* @param check_momenta Special treatment for partons inside extremal jets
	* @returns The HEJ matrix element including virtual corrections
	*
	* cf. eq. (22) in \ref Andersen:2011hs
	* Incoming particles should be ordered by ascending z momentum.
	* Outgoing particles should be ordered by ascending rapidity.
	*
	* \internal Relation to standard HEJ Met2: MatrixElement = Met2shat^2/(pdftapdftb)
	*/
	double operator()(
	double mur,
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const;

	//! HEJ tree-level matrix element
	/**
	* @param mur Value of the renormalisation scale
	* @param incoming Incoming particles
	* @param outgoing Outgoing particles
	* @param check_momenta Special treatment for partons inside extremal jets
	* @returns The HEJ matrix element without virtual corrections
	*
	* cf. eq. (22) in \ref Andersen:2011hs
	* Incoming particles should be ordered by ascending z momentum.
	* Outgoing particles should be ordered by ascending rapidity.
	*/
	double tree(
	double mur,
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const;

	//! HEJ tree-level matrix element - parametric part
	/**
	* @param mur Value of the renormalisation scale
	* @param outgoing Outgoing particles
	* @returns The parametric part of the tree matrix element
	*
	* cf. eq. (22) in \ref Andersen:2011hs
	*
	* The tree level matrix element factorises into a parametric part
	* which depends on the theory parameters (alpha_s and scale)
	* and a kinematic part comprising the dependence on the particle momenta
	* and colour factors. This function returns the former.
	*/
	double tree_param(
	double mur,
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing
	) const;

	//! HEJ tree-level matrix element - kinematic part
	/**
	* @param incoming Incoming particles
	* @param partons Outgoing particles
	* @param check_momenta Special treatment for partons inside extremal jets
	* @returns The kinematic part of the tree matrix element
	*
	* cf. eq. (22) in \ref Andersen:2011hs
	* Incoming particles should be ordered by ascending z momentum.
	* Outgoing particles should be ordered by ascending rapidity.
	*
	* The tree level matrix element factorises into a parametric part
	* which depends on the theory parameters (alpha_s and scale)
	* and a kinematic part comprising the dependence on the particle momenta
	* and colour factors. This function returns the latter.
	*/
	double tree_kin(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const;

	/**
	* \brief Calculates the Virtual Corrections
	* @param mur Value of the renormalisation scale
	* @param in Incoming particles
	* @param out Outgoing particles
	* @returns The Virtual Corrections of the Matrix Element
	*
	* Incoming particles should be ordered by ascending z momentum.
	* Outgoing particles should be ordered by ascending rapidity.
	*
	* 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 \ref Andersen:2011hs
	* note that indices are off by one, i.e. out[0].p corresponds to p_1
	*/
	double virtual_corrections(
	double mur,
	std::array<Sparticle, 2> const & in,
	std::vector<Sparticle> const & out
	) const;

	private:
	/**
	* \brief cf. last line of eq. (22) in \ref Andersen:2011hs
	* @param mur Value of Renormalisation Scale
	* @param q_j ???
	* @param lambda ???
	*/
	double omega0(
	double alpha_s, double mur,
	fastjet::PseudoJet const & q_j, double lambda
	) const;

	double tree_kin_jets(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> partons,
	bool check_momenta
	) const;
	double tree_kin_Higgs(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const;
	double tree_kin_Higgs_first(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const;
	double tree_kin_Higgs_last(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const;

	/**
	* \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(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const;


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


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


	std::vector<Sparticle> tag_extremal_jet_partons(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> out_partons, bool check_momenta
	) const;

	double MH2_forwardH(
	RHEJ::ParticleID id,
	CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector pH,
	double t1, double t2
	) const;

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

	- MatrixElementParameters param_;
	+ MatrixElementConfig param_;
	};


	}
	diff --git a/include/RHEJ/PhaseSpacePoint.hh b/include/RHEJ/PhaseSpacePoint.hh
	index 59364c7..b050bdb 100644
	--- a/include/RHEJ/PhaseSpacePoint.hh
	+++ b/include/RHEJ/PhaseSpacePoint.hh
	@@ -1,157 +1,157 @@
	/** \file PhaseSpacePoint.hh
	* \brief Contains the PhaseSpacePoint Class
	*/

	#pragma once

	#include <vector>

	#include <CLHEP/Random/Randomize.h>
	#include <CLHEP/Random/RanluxEngine.h>

	#include "RHEJ/utility.hh"
	#include "RHEJ/config.hh"
	#include "RHEJ/Event.hh"
	#include "RHEJ/Splitter.hh"

	namespace RHEJ{

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

	//! PhaseSpacePoint Constructor
	/**
	* @param ev Clustered Jet Event
	* @param conf Configuration parameters
	*/
	PhaseSpacePoint(
	Event const & ev,
	- Config const & conf
	+ PhaseSpacePointConfig conf
	);

	//! Get Weight Function
	/**
	* @returns Weight of Event
	*/
	double weight() const{
	return weight_;
	}


	//! Get Incoming Function
	/**
	* @returns Incoming Particles
	*/
	std::array<Sparticle, 2> const & incoming() const{
	return incoming_;
	}

	//! Get Outgoing Function
	/**
	* @returns Outgoing Particles
	*/
	std::vector<Sparticle> const & outgoing() const{
	return outgoing_;
	}

	std::unordered_map<int, std::vector<Sparticle>> const & decays() const{
	return decays_;
	}

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


	static void reset_ranlux(std::string const & init_file);
	static void reset_ranlux(char const * init_file);

	private:

	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_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
	);
	bool jets_ok(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	void reconstruct_incoming(std::array<Sparticle, 2> const & Born_incoming);
	double phase_space_normalisation(
	int num_Born_jets,
	int num_res_partons
	) const;
	std::vector<fastjet::PseudoJet> split(
	std::vector<fastjet::PseudoJet> const & jets, int ng_jets
	);
	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
	);
	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;
	void copy_AWZH_boson_from(Event const & event);

	bool momentum_conserved() const;

	bool unob_, unof_;

	double weight_;

	- PhaseSpacePointParameters param_;
	+ PhaseSpacePointConfig param_;


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

	static CLHEP::Ranlux64Engine ran_;
	HejSplit splitter_;
	};

	}
	diff --git a/include/RHEJ/config.hh b/include/RHEJ/config.hh
	index 26d031c..dfbbd69 100644
	--- a/include/RHEJ/config.hh
	+++ b/include/RHEJ/config.hh
	@@ -1,194 +1,185 @@
	/** \file config.hh
	* \brief The file which handles the configuration file parameters
	*
	* Contains the JetParameters Struct and ScaleConfig Struct. Also
	* contains the ScaleGenerator Class and the EventTreatment class.
	* Config struct is also defined here. The configuration files
	* parameters are read and then stored within this objects.
	*/

	#pragma once

	#include <string>
	#include <memory>

	#include "fastjet/JetDefinition.hh"
	#include "yaml-cpp/yaml.h"

	#include "RHEJ/event_types.hh"
	#include "RHEJ/Event.hh"
	#include "RHEJ/HiggsCouplingSettings.hh"
	#include "RHEJ/optional.hh"
	#include "RHEJ/output_formats.hh"
	#include "RHEJ/ScaleFunction.hh"
	#include "RHEJ/utility.hh"


	namespace RHEJ{

	/**! \struct JetParameters config.hh "include/RHEJ/config.hh"
	* \brief Contains Jet Definition and Minimum Pt to be a Jet.
	*
	* Contains the Fastjet definition of a Jet and a threshold (minimum) value for pt
	*/
	struct JetParameters{
	fastjet::JetDefinition def; /*< Jet Definition /
	double min_pt; /*< Minimum Jet Pt /
	};

	/**! \struct ScaleConfig config.hh "include/RHEJ/config.hh"
	* \brief Sets up the possible choices for scale variation?
	*
	* There is a vector which contains the possible scale choices and
	* also vector of doubles with the scale factors along with a max
	* scale ratio.
	*/
	struct ScaleConfig{
	std::vector<std::shared_ptr<ScaleFunction>> scales; /*< Vector of possible Scale choices /
	std::vector<double> scale_factors; /*< Vector of possible Scale Factors /
	double max_scale_ratio; /*< Maximum ratio for the scales /
	};

	/**! \class ScaleGenerator config.hh "include/RHEJ/config.hh"
	* \brief A Class used to generate scales
	*
	* This class has a clustered event class and scale config struct
	* within it.
	*/
	class ScaleGenerator{
	public:
	ScaleGenerator() = default;
	explicit ScaleGenerator(ScaleConfig && settings):
	settings_{std::move(settings)}
	{}

	Event operator()(Event ev) const;
	ScaleConfig const & settings() const;

	private:
	ScaleConfig settings_;
	};


	/**! \enum EventTreatment config.hh "include/RHEJ/config.hh"
	* \brief Enumeration which deals with how to treat events.
	*
	* The program will decide on how to treat an event based on
	* the value of this enumeration.
	*/
	enum class EventTreatment{
	reweight, /*< Reweigh the event /
	keep, /*< Keep the event /
	discard, /*< Discard the event /
	};

	/**
	* \brief An ordered Map with EventType keys and mapped values EventTreatment
	*
	* If called upon with EventTreatMap.at(EventType) it will return the treatment which has
	* been specified for that EventType.
	*/
	using EventTreatMap = std::map<event_type::EventType, EventTreatment>;

	/**! \struct Config config.hh "include/RHEJ/config.hh"
	* \brief Config Struct for user input parameters.
	*
	* The struct which handles the input card given by the user.
	*
	* \internal To add a new option:
	* 1. Add a member to the Config struct.
	* 2. Inside "src/YAMLreader.cc":
	* - Add the option name to the "supported" Node in
	* get_supported_options.
	* - Initialise the new Config member in to_Config.
	* The functions set_from_yaml (for mandatory options) and
	* set_from_yaml_if_defined (non-mandatory) may be helpful.
	* 3. Add a new entry (with short description) to config.yaml
	*/
	struct Config {
	ScaleGenerator scale_gen; /*< Scale /
	JetParameters resummation_jets; /*< Resummation Jets /
	JetParameters fixed_order_jets; /*< Fixed-Order Jets /
	double min_extparton_pt; /*< Min External Parton Pt/
	double max_ext_soft_pt_fraction; /*< Min External Parton Pt Fraction /
	int trials; /*< Number of resummation points to generate per FO /
	//! Log Correct from running of \f$\alpha_s$\f On or Off
	bool log_correction;
	bool unweight; /*< Unweight On or Off /
	std::vector<OutputFile> output; /*< Output Files Type /
	optional<std::string> RanLux_init; /*< Ranlux File Choice/
	//! Map to decide what to do for differnt EventTypes
	EventTreatMap treat;
	YAML::Node analysis_parameters; /*< Analysis Parameters /
	//! Settings for Effective Higgs-Gluon Coupling
	HiggsCouplingSettings Higgs_coupling;
	};

	- class Parameters{
	- protected:
	- const Config & config_;
	- public:
	- explicit Parameters(const Config & conf): config_(conf){}
	- const Config & get_Config() const {return config_;}
	+ //! Configuration options for the PhaseSpacePoint class
	+ struct PhaseSpacePointConfig {
	+ JetParameters jet_param;
	+ double min_extparton_pt;
	+ double max_ext_soft_pt_fraction;
	};

	- class PhaseSpacePointParameters: public Parameters {
	- public:
	- explicit PhaseSpacePointParameters(const Config & conf): Parameters(conf){}
	- ///@{
	- /// @name getter functions
	- const JetParameters & jet_parameters() const
	- {return config_.resummation_jets;}
	- double jet_ptmin() const
	- {return jet_parameters().min_pt;}
	- const fastjet::JetDefinition & jet_def() const
	- {return jet_parameters().def;}
	- double jet_R() const
	- {return jet_def().R();}
	-
	- double min_ext_pt() const
	- {return config_.min_extparton_pt;}
	- double max_ext_soft_pt_fraction() const
	- {return config_.max_ext_soft_pt_fraction;}
	- ///@}
	+ //! Configuration options for the MatrixElement class
	+ struct MatrixElementConfig {
	+ JetParameters jet_param;
	+ bool log_correction;
	+ HiggsCouplingSettings Higgs_coupling;
	};

	- class MatrixElementParameters: public Parameters {
	- public:
	- explicit MatrixElementParameters(const Config & conf):
	- Parameters(conf){}
	- ///@{
	- /// @name getter functions
	- const JetParameters & jet_parameters() const
	- {return config_.resummation_jets;}
	- double jet_ptmin() const
	- {return jet_parameters().min_pt;}
	- const fastjet::JetDefinition & jet_def() const
	- {return jet_parameters().def;}
	-
	- bool log_corr() const
	- {return config_.log_correction;}
	-
	- const HiggsCouplingSettings & Hgg_setting() const
	- {return config_.Higgs_coupling;}
	- ///@}
	+ //! Configuration options for the EventReweighter class
	+ struct EventReweighterConfig {
	+ PhaseSpacePointConfig psp_config;
	+ MatrixElementConfig ME_config;
	+ JetParameters jet_param;
	+ EventTreatMap treat;
	};

	- class EventReweighterParameters: public Parameters {
	- public:
	- explicit EventReweighterParameters(const Config & conf):
	- Parameters(conf){}
	- ///@{
	- /// @name getter functions
	- const JetParameters & jet_parameters() const
	- {return config_.resummation_jets;}
	- double jet_ptmin() const
	- {return jet_parameters().min_pt;}
	- const fastjet::JetDefinition & jet_def() const
	- {return jet_parameters().def;}
	-
	- const EventTreatMap & treat() const
	- {return config_.treat;}
	- ///@}
	- };
	+ /**! Extract PhaseSpacePointConfig from Config
	+ *
	+ * We do not provide a PhaseSpacePointConfig constructor from Config
	+ * so that PhaseSpacePointConfig remains an aggregate.
	+ * This faciliates writing client code (e.g. the HEJ fixed-order generator)
	+ * that creates a PhaseSpacePointConfig without a Config object.
	+ *
	+ * @see to_MatrixElementConfig, to_EventReweighterConfig
	+ */
	+ inline
	+ PhaseSpacePointConfig to_PhaseSpacePointConfig(Config const & conf) {
	+ return {conf.resummation_jets, conf.min_extparton_pt, conf.max_ext_soft_pt_fraction};
	+ }
	+
	+ /**! Extract MatrixElementConfig from Config
	+ *
	+ * @see to_PhaseSpacePointConfig, to_EventReweighterConfig
	+ */
	+ inline
	+ MatrixElementConfig to_MatrixElementConfig(Config const & conf) {
	+ return {conf.resummation_jets, conf.log_correction, conf.Higgs_coupling};
	+ }
	+
	+ /**! Extract EventReweighterConfig from Config
	+ *
	+ * @see to_PhaseSpacePointConfig, to_MatrixElementConfig
	+ */
	+ inline
	+ EventReweighterConfig to_EventReweighterConfig(Config const & conf) {
	+ return {
	+ to_PhaseSpacePointConfig(conf),
	+ to_MatrixElementConfig(conf),
	+ conf.resummation_jets, conf.treat
	+ };
	+ }
	+
	} // namespace RHEJ
	diff --git a/include/RHEJ/make_writer.hh b/include/RHEJ/make_writer.hh
	index 8fb1492..c27f2a7 100644
	--- a/include/RHEJ/make_writer.hh
	+++ b/include/RHEJ/make_writer.hh
	@@ -1,32 +1,30 @@
	/** \file make_writer.hh
	* \brief Contains make_format_writer which returns the EventWriter which is asked for.
	*/

	#pragma once

	#include <string>
	#include <memory>

	#include "RHEJ/output_formats.hh"
	#include "RHEJ/EventWriter.hh"

	namespace LHEF{
	struct HEPRUP;
	}

	namespace RHEJ{
	- struct Config;

	/**
	* \brief returns the EventWriter which is asked for.
	*
	* Cases which are possible: Les_Houches and HepMC
	*/
	std::unique_ptr<EventWriter> make_format_writer(
	FileFormat format, /*< Output File Type /
	std::string const & outfile, /*< Output File Name /
	- Config const & conf, /*< Config for running HEJ /
	LHEF::HEPRUP const & heprup /*< Contains important information from the input LHEF /
	);

	}
	diff --git a/src/CombinedEventWriter.cc b/src/CombinedEventWriter.cc
	index 883958d..5b1cf38 100644
	--- a/src/CombinedEventWriter.cc
	+++ b/src/CombinedEventWriter.cc
	@@ -1,24 +1,21 @@
	#include "RHEJ/CombinedEventWriter.hh"

	-#include "RHEJ/config.hh"
	-
	namespace RHEJ{

	CombinedEventWriter::CombinedEventWriter(
	- Config const & conf,
	+ std::vector<OutputFile> const & outfiles,
	LHEF::HEPRUP const & heprup
	){
	- std::vector<OutputFile> const & outfiles = conf.output;
	writers_.reserve(outfiles.size());
	for(OutputFile const & outfile: outfiles){
	writers_.emplace_back(
	- make_format_writer(outfile.format, outfile.name, conf, heprup)
	+ make_format_writer(outfile.format, outfile.name, heprup)
	);
	}
	}

	void CombinedEventWriter::write(Event const & ev){
	for(auto & writer: writers_) writer->write(ev);
	}

	}
	diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
	index 3f9672d..1d7a8d1 100644
	--- a/src/EventReweighter.cc
	+++ b/src/EventReweighter.cc
	@@ -1,302 +1,303 @@
	#include "RHEJ/EventReweighter.hh"

	#include <string>
	#include <unordered_map>

	#include "RHEJ/PhaseSpacePoint.hh"
	#include "RHEJ/PDG_codes.hh"
	#include "RHEJ/debug.hh"

	namespace RHEJ{
	using EventType = event_type::EventType;

	namespace {

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

	UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
	UnclusteredEvent result;
	result.incoming = psp.incoming();
	std::sort(
	begin(result.incoming), end(result.incoming),
	[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
	);
	assert(result.incoming.size() == 2);
	result.outgoing = psp.outgoing();
	assert(
	std::is_sorted(
	begin(result.outgoing), end(result.outgoing),
	rapidity_less{}
	)
	);
	assert(result.outgoing.size() >= 2);
	result.decays = psp.decays();
	result.central.mur = NaN;
	result.central.muf = NaN;
	result.central.weight = psp.weight();
	return result;
	}

	} // namespace anonymous

	EventReweighter::EventReweighter(
	LHEF::HEPRUP const & heprup,
	- Config const & conf
	+ EventReweighterConfig conf
	):
	EventReweighter{
	RHEJ::Beam{
	heprup.EBMUP.first,
	{{
	static_cast<RHEJ::ParticleID>(heprup.IDBMUP.first),
	static_cast<RHEJ::ParticleID>(heprup.IDBMUP.second)
	}}
	},
	heprup.PDFSUP.first,
	- conf
	+ std::move(conf)
	}
	{
	if(heprup.EBMUP.second != E_beam_){
	throw std::invalid_argument(
	"asymmetric beam: " + std::to_string(E_beam_)
	+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
	);
	};
	if(heprup.PDFSUP.second != pdf_.id()){
	throw std::invalid_argument(
	"conflicting PDF ids: " + std::to_string(pdf_.id())
	+ " vs. " + std::to_string(heprup.PDFSUP.second)
	);
	}
	}

	EventReweighter::EventReweighter(
	Beam beam,
	int pdf_id,
	- Config const & conf
	+ EventReweighterConfig conf
	):
	- param_(conf),
	+ param_{std::move(conf)},
	E_beam_{beam.E},
	pdf_{pdf_id, beam.type.front(), beam.type.back()},
	MEt2_{
	[this](double mu){ return pdf_.Halphas(mu); },
	- param_.get_Config()
	+ param_.ME_config
	}
	{}

	PDF const & EventReweighter::pdf() const{
	return pdf_;
	}

	namespace{
	static_assert(
	std::numeric_limits<double>::has_infinity, "infinity not supported"
	);

	// there is no simple way to create a vector of move-only objects
	// this is a clumsy work-around
	ScaleGenerator make_identity(){
	std::vector<std::shared_ptr<ScaleFunction>> scale_fun;
	scale_fun.emplace_back(new InputScales);
	return ScaleGenerator{
	ScaleConfig{
	std::move(scale_fun), {}, std::numeric_limits<double>::infinity()
	}
	};
	}

	const ScaleGenerator identity{make_identity()};
	}

	std::vector<Event> EventReweighter::reweight(
	Event const & input_ev, int num_events
	){
	return reweight(input_ev, num_events, identity);
	}

	std::vector<Event> EventReweighter::reweight(
	Event const & input_ev, int num_events,
	ScaleGenerator const & gen_scales
	){
	auto res_events = gen_res_events(input_ev, num_events);
	if(res_events.empty()) return {};
	for(auto & event: res_events) event = gen_scales(event);
	return rescale(input_ev, std::move(res_events));
	}

	/**
	* \brief main generation/reweighting function:
	* generate phase space points and divide out Born factors
	*/
	std::vector<Event> EventReweighter::gen_res_events(
	Event const & ev,
	int phase_space_points
	){
	assert(ev.variations().empty());

	- switch(param_.treat().at(ev.type())){
	+ switch(param_.treat.at(ev.type())){
	case EventTreatment::discard: return {};
	case EventTreatment::keep:
	if(! jets_pass_resummation_cuts(ev)) return {};
	else return {ev};
	default:;
	}
	const double Born_shat = shat(ev);

	std::vector<Event> resummation_events;
	for(int psp_number = 0; psp_number < phase_space_points; ++psp_number){
	PhaseSpacePoint psp{
	ev,
	- param_.get_Config()
	+ param_.psp_config
	};
	if(psp.weight() == 0.) continue;
	if(psp.incoming()[0].E() > E_beam_ \|\| psp.incoming()[1].E() > E_beam_) continue;

	resummation_events.emplace_back(
	- to_UnclusteredEvent(std::move(psp)), param_.jet_def(), param_.jet_ptmin()
	+ to_UnclusteredEvent(std::move(psp)),
	+ param_.jet_param.def, param_.jet_param.min_pt
	);
	auto & new_event = resummation_events.back();
	assert(new_event.variations().empty());
	new_event.central().mur = ev.central().mur;
	new_event.central().muf = ev.central().muf;
	const double resum_shat = shat(new_event);
	new_event.central().weight = ev.central().weightBorn_shat*Born_shat/
	(phase_space_pointsresum_shatresum_shat);
	}
	return resummation_events;
	}

	std::vector<Event> EventReweighter::rescale(
	Event const & Born_ev,
	std::vector<Event> events
	) const{
	const double Born_pdf = pdf_factors(Born_ev).central;
	const double Born_ME = tree_matrix_element(Born_ev);

	for(auto & cur_event: events){
	const auto pdf = pdf_factors(cur_event);
	assert(pdf.variations.size() == cur_event.variations().size());
	const auto ME = matrix_elements(cur_event);
	assert(ME.variations.size() == cur_event.variations().size());
	cur_event.central().weight = pdf.centralME.central/(Born_pdf*Born_ME);
	for(size_t i = 0; i < cur_event.variations().size(); ++i){
	cur_event.variations(i).weight *=
	pdf.variations[i]ME.variations[i]/(Born_pdfBorn_ME);
	}
	}
	return events;
	};

	bool EventReweighter::jets_pass_resummation_cuts(
	Event const & ev
	) const{
	const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
	- fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_def()};
	- return cs.inclusive_jets(param_.jet_ptmin()).size() == ev.jets().size();
	+ fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_param.def};
	+ return cs.inclusive_jets(param_.jet_param.min_pt).size() == ev.jets().size();
	}



	EventReweighter::EventFactors
	EventReweighter::pdf_factors(Event const & ev) const{
	auto const & a = ev.incoming().front();
	auto const & b = ev.incoming().back();
	const double xa = a.p.e()/E_beam_;
	const double xb = b.p.e()/E_beam_;

	EventFactors result;
	std::unordered_map<double, double> known_pdf;
	result.central =
	pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
	pdf_.pdfpt(1,xb,ev.central().muf,b.type);
	known_pdf.emplace(ev.central().muf, result.central);

	result.variations.reserve(ev.variations().size());
	for(auto const & ev_param: ev.variations()){
	const double muf = ev_param.muf;
	auto cur_pdf = known_pdf.find(muf);
	if(cur_pdf == known_pdf.end()){
	cur_pdf = known_pdf.emplace(
	muf,
	pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
	).first;
	}
	result.variations.emplace_back(cur_pdf->second);
	}
	assert(result.variations.size() == ev.variations().size());
	return result;
	}




	EventReweighter::EventFactors
	EventReweighter::matrix_elements(Event const & ev) const{
	- assert(param_.treat().count(ev.type()) > 0);
	- if(param_.treat().find(ev.type())->second == EventTreatment::keep){
	+ assert(param_.treat.count(ev.type()) > 0);
	+ if(param_.treat.find(ev.type())->second == EventTreatment::keep){
	return fixed_order_scale_ME(ev);
	}

	// precompute overall kinematic factor
	const double ME_kin = MEt2_.tree_kin(ev.incoming(), ev.outgoing(), true);

	EventFactors result;
	std::unordered_map<double, double> known_ME;
	result.central = MEt2_(
	ev.central().mur,
	ev.incoming(), ev.outgoing(),
	true
	);
	known_ME.emplace(ev.central().mur, result.central);

	result.variations.reserve(ev.variations().size());
	- for(auto const & param_: ev.variations()){
	- const double mur = param_.mur;
	+ for(auto const & param: ev.variations()){
	+ const double mur = param.mur;
	auto cur_ME = known_ME.find(mur);
	if(cur_ME == known_ME.end()){
	const double ME = MEt2_.tree_param(
	mur, ev.incoming(), ev.outgoing()
	)ME_kinMEt2_.virtual_corrections(
	mur, ev.incoming(), ev.outgoing()
	);
	cur_ME = known_ME.emplace(mur, ME).first;
	}
	result.variations.emplace_back(cur_ME->second);
	}
	assert(result.variations.size() == ev.variations().size());
	return result;
	}

	double EventReweighter::tree_matrix_element(Event const & ev) const{
	assert(ev.variations().empty());
	- assert(param_.treat().count(ev.type()) > 0);
	- if(param_.treat().find(ev.type())->second == EventTreatment::keep){
	+ assert(param_.treat.count(ev.type()) > 0);
	+ if(param_.treat.find(ev.type())->second == EventTreatment::keep){
	return fixed_order_scale_ME(ev).central;
	}
	return MEt2_.tree(
	ev.central().mur,
	ev.incoming(), ev.outgoing(),
	false
	);
	}

	EventReweighter::EventFactors
	EventReweighter::fixed_order_scale_ME(Event const & ev) const{
	const int alpha_s_power = std::count_if(
	begin(ev.outgoing()), end(ev.outgoing()),
	[](Sparticle const & p){ return is_parton(p); }
	);
	EventFactors result;
	result.central = pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
	for(auto const & var: ev.variations()){
	result.variations.emplace_back(
	pow(pdf_.Halphas(var.mur), alpha_s_power)
	);
	}
	return result;
	}

	}
	diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc
	index 4e9353b..163f298 100644
	--- a/src/HepMCWriter.cc
	+++ b/src/HepMCWriter.cc
	@@ -1,256 +1,251 @@
	#include "RHEJ/HepMCWriter.hh"

	#include <cassert>

	#ifdef RHEJ_BUILD_WITH_HepMC_VERSION

	#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
	#include "HepMC/WriterAscii.h"
	#include "HepMC/LHEFAttributes.h"
	#else
	#include "HepMC/IO_GenEvent.h"
	#endif

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

	-#include "RHEJ/config.hh"
	#include "RHEJ/utility.hh"
	#include "RHEJ/Event.hh"
	#include "RHEJ/exceptions.hh"
	#include "RHEJ/Version.hh"

	namespace RHEJ{

	#if RHEJ_BUILD_WITH_HepMC_VERSION >= 3
	namespace {
	void add_generator_tag(LHEF::HEPRUP & heprup){
	heprup.generators.emplace_back(LHEF::XMLTag{});
	heprup.generators.back().name = Version::package_name();
	heprup.generators.back().version = Version::String();
	}
	void reset_weight_info(LHEF::HEPRUP & heprup){
	heprup.IDWTUP = -4;
	// use placeholders for unknown init block values
	// we can overwrite them after processing all events
	heprup.XSECUP = {0.};
	heprup.XERRUP = {0.};
	heprup.XMAXUP = {0.};
	}
	HepMC::FourVector to_FourVector(Sparticle const & sp){
	return {sp.px(), sp.py(), sp.pz(), sp.E()};
	}

	HepMC::shared_ptr<HepMC::GenRunInfo> & add_HEPRUP(
	HepMC::shared_ptr<HepMC::GenRunInfo> & runinfo,
	LHEF::HEPRUP heprup
	){
	add_generator_tag(heprup);
	reset_weight_info(heprup);

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

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

	- struct HepMCWriter::HepMCWriterImpl: public EventWriter{
	+ struct HepMCWriter::HepMCWriterImpl{

	- Config config_;
	HepMC::shared_ptr<HepMC::GenRunInfo> runinfo_;
	HepMC::WriterAscii writer_;

	HepMCWriterImpl(
	- std::string const & file, Config const & conf, LHEF::HEPRUP && heprup
	+ std::string const & file, LHEF::HEPRUP && heprup
	):
	- config_(conf),
	runinfo_{HepMC::make_shared<HepMC::GenRunInfo>()},
	writer_{file, add_HEPRUP(runinfo_, std::move(heprup))}
	{}

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

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

	void write(Event const & ev){
	if(!ev.decays().empty()){
	throw not_implemented{"HepMC output for events with decays"};
	}
	auto vx = HepMC::make_shared<HepMC::GenVertex>();
	for(auto const & in: ev.incoming()){
	vx->add_particle_in(
	HepMC::make_shared<HepMC::GenParticle>(
	to_FourVector(in), static_cast<int>(in.type), -1
	)
	);
	}
	for(auto const & out: ev.outgoing()){
	vx->add_particle_out(
	HepMC::make_shared<HepMC::GenParticle>(
	to_FourVector(out), static_cast<int>(out.type), +1
	)
	);
	}
	HepMC::GenEvent out_ev{runinfo_, HepMC::Units::GEV, HepMC::Units::MM};
	out_ev.add_vertex(vx);
	out_ev.weights().reserve(ev.variations().size() + 1u);
	out_ev.weights().emplace_back(ev.central().weight);
	for(auto const & var: ev.variations()){
	out_ev.weights().emplace_back(var.weight);
	}
	writer_.write_event(out_ev);
	}
	};

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

	#else // HepMC 2

	namespace{
	HepMC::FourVector to_FourVector(Sparticle const & sp){
	return {sp.px(), sp.py(), sp.pz(), sp.E()};
	}
	}

	- struct HepMCWriter::HepMCWriterImpl: public EventWriter{
	+ struct HepMCWriter::HepMCWriterImpl{

	struct HepMCInfo{
	size_t event_count = 0.;
	double tot_weight = 0.;
	double tot_weight2 = 0.;

	void add_event(Event const & ev){
	double wt = ev.central().weight;
	tot_weight += wt;
	tot_weight2 += wt * wt;
	++event_count;
	}
	};

	- Config config_;
	HepMC::IO_GenEvent writer_;
	HepMCInfo genInfo_;

	HepMCWriterImpl(
	- std::string const & file, Config const & conf, LHEF::HEPRUP &&
	+ std::string const & file, LHEF::HEPRUP &&
	):
	- config_(conf),
	writer_{file}
	{}

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

	void write(Event const & ev){
	auto vx = new HepMC::GenVertex();
	for(auto const & in: ev.incoming()){
	vx->add_particle_in(
	new HepMC::GenParticle{
	to_FourVector(in), static_cast<int>(in.type), -1
	}
	);
	}
	for(int i=0; i<ev.outgoing().size(); ++i){
	auto const & out = ev.outgoing()[i];
	bool decays = false;
	for(auto const & decay: ev.decays()){
	if( i == decay.first ){ // if particle decays replace it with decay products
	for( auto const & out: decay.second){
	vx->add_particle_out(
	new HepMC::GenParticle{
	to_FourVector(out), static_cast<int>(out.type), +1
	}
	);
	}
	decays = true;
	break;
	}
	}
	if(!decays)
	vx->add_particle_out(
	new HepMC::GenParticle{
	to_FourVector(out), static_cast<int>(out.type), +1
	}
	);
	}
	HepMC::GenEvent out_ev{HepMC::Units::GEV, HepMC::Units::MM};
	out_ev.add_vertex(vx);
	/// weights
	out_ev.weights().push_back(ev.central().weight);
	for(auto const & var: ev.variations()){
	out_ev.weights().push_back(var.weight);
	}
	/// @TODO add name list
	/// general informations
	genInfo_.add_event(ev);
	out_ev.set_signal_process_id(ev.type()+1); /// "+1": conistent with lhe
	out_ev.set_event_scale(ev.central().mur); /// event scale
	out_ev.set_event_number(genInfo_.event_count); /// event number
	/// cross section
	HepMC::GenCrossSection t_xs;
	t_xs.set_cross_section( genInfo_.tot_weight , sqrt(genInfo_.tot_weight2) );
	out_ev.set_cross_section( t_xs );

	/// @TODO add alphaQCD (need function) and alphaQED
	/// @TODO output pdf (currently not avaiable from event alone)
	writer_.write_event(&out_ev);
	}
	};

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

	#endif // HepMC 2

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

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

	#else // no HepMC

	namespace RHEJ{

	class HepMCWriter::HepMCWriterImpl{};

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

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

	void HepMCWriter::HepMCWriterImplDeleter::operator()(HepMCWriterImpl* p) {
	delete p;
	}
	}
	#endif
	diff --git a/src/LesHouchesWriter.cc b/src/LesHouchesWriter.cc
	index 7561dad..5f0d6c6 100644
	--- a/src/LesHouchesWriter.cc
	+++ b/src/LesHouchesWriter.cc
	@@ -1,92 +1,89 @@
	#include <stdexcept>
	#include <memory>
	#include <cassert>

	#include "RHEJ/LesHouchesWriter.hh"
	#include "RHEJ/event_types.hh"

	-#include "RHEJ/config.hh"
	#include "RHEJ/Event.hh"
	-#include "RHEJ/Version.hh"

	namespace RHEJ{
	namespace{
	template<class T, class... Args>
	std::unique_ptr<T> make_unique(Args&&... a){
	return std::unique_ptr<T>{new T{std::forward<Args>(a)...}};
	}
	}

	LesHouchesWriter::LesHouchesWriter(
	- std::string const & file, Config const & conf, LHEF::HEPRUP heprup
	+ std::string const & file, LHEF::HEPRUP heprup
	):
	- config_(conf),
	out_{file, std::fstream::in \| std::fstream::out \| std::fstream::trunc},
	writer_{RHEJ::make_unique<LHEF::Writer>(out_)}
	{
	if(! out_.is_open()){
	throw std::ios_base::failure("Failed to open " + file);
	};
	writer_->heprup = std::move(heprup);
	// lhe Stardard: IDWTUP (negative => weights = +/-)
	// 3: weight=+/-, xs given in head (same as default MG)
	// 4: weight=+/-, xs = avg(weights)
	writer_->heprup.IDWTUP = -3;
	writer_->heprup.generators.emplace_back(LHEF::XMLTag{});
	- writer_->heprup.generators.back().name = Version::package_name();
	- writer_->heprup.generators.back().version = Version::String();
	+ writer_->heprup.generators.back().name = "HEJ";
	+ writer_->heprup.generators.back().version = "0.0.1";
	// use placeholders for unknown init block values
	// we can overwrite them after processing all events
	writer_->heprup.XSECUP = std::vector<double>(event_type::last_type+1, 0.);
	writer_->heprup.XERRUP = std::vector<double>(event_type::last_type+1, 0.);
	writer_->heprup.XMAXUP = std::vector<double>(event_type::last_type+1, 0.);
	write_init();
	}

	void LesHouchesWriter::write(Event const & ev){
	assert(writer_ && out_.is_open());

	const double wt = ev.central().weight;
	writer_->hepeup = RHEJ::to_HEPEUP(std::move(ev), &heprup());
	writer_->writeEvent();
	heprup().XSECUP[ev.type()] += wt;
	heprup().XERRUP[ev.type()] += wt*wt;
	if(wt > heprup().XMAXUP[ev.type()]){
	heprup().XMAXUP[ev.type()] = wt;
	}
	}

	// this function is called after overwritting the Les Houches init block
	// assert that we have overwritten exactly the init block,
	// i.e. an intact event block is next
	void assert_next_event_intact(std::istream & out){
	(void) out; // suppress compiler warnings if not in debug mode
	#ifndef NDEBUG
	std::string line;
	getline(out, line);
	assert(line == "<event>");
	#endif
	}

	void LesHouchesWriter::rewrite_init(){
	assert(writer_ && out_.is_open());

	// replace placeholder entries
	const auto pos = out_.tellp();
	out_.seekp(0);
	writer_->init();
	assert_next_event_intact(out_);
	out_.seekp(pos);
	}

	LesHouchesWriter::~LesHouchesWriter(){
	assert(writer_ && out_.is_open());

	for(auto & xs_err: heprup().XERRUP)
	{
	xs_err = sqrt(xs_err);
	}
	rewrite_init();
	}

	}
	diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
	index ed5cc21..84815eb 100644
	--- a/src/MatrixElement.cc
	+++ b/src/MatrixElement.cc
	@@ -1,749 +1,752 @@
	#include "RHEJ/MatrixElement.hh"

	#include <CLHEP/Random/Randomize.h>
	#include <CLHEP/Random/RanluxEngine.h>

	#include "RHEJ/Constants.hh"
	#include "RHEJ/currents.hh"
	#include "RHEJ/PDG_codes.hh"
	#include "RHEJ/uno.hh"
	#include "RHEJ/debug.hh"

	namespace RHEJ{
	//cf. last line of eq. (22) in \ref Andersen:2011hs
	double MatrixElement::omega0(
	double alpha_s, double mur,
	fastjet::PseudoJet const & q_j, double lambda
	) const {
	const double result = - alpha_sN_C/M_PIlog(q_j.perp2()/(lambda*lambda));
	- if(! param_.log_corr()) return result;
	+ if(! param_.log_correction) return result;
	// use alpha_s(sqrt(q_j*lambda)), evolved to mur
	return (
	1. + alpha_s/(4.M_PI)beta0log(murmur/(q_j.perp()*lambda))
	)*result;
	}

	double MatrixElement::virtual_corrections(
	double mur,
	std::array<Sparticle, 2> const & in,
	std::vector<Sparticle> const & out
	) const{
	fastjet::PseudoJet const & pa = in.front().p;
	#ifndef NDEBUG
	fastjet::PseudoJet const & pb = in.back().p;
	double const norm = (in.front().p + in.back().p).E();
	#endif

	assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
	assert(out.size() >= 2);
	assert(pa.pz() < pb.pz());

	fastjet::PseudoJet q = pa - out[0].p;
	size_t first_idx = 0;
	size_t last_idx = out.size() - 1;
	// if there is a Higgs or unordered gluon outside the extremal partons
	// then it is not part of the FKL ladder and does not contribute
	// to the virtual corrections
	if(out.front().type == pid::Higgs \|\| has_unob_gluon(in, out)){
	q -= out[1].p;
	++first_idx;
	}
	if(out.back().type == pid::Higgs \|\| has_unof_gluon(in, out)){
	--last_idx;
	}

	double exponent = 0;
	const double alpha_s = alpha_s_(mur);
	for(size_t j = first_idx; j < last_idx; ++j){
	exponent += omega0(alpha_s, mur, q, CLAMBDA)*(
	out[j+1].rapidity() - out[j].rapidity()
	);
	q -= out[j+1].p;
	}
	assert(
	nearby(q, -1*pb, norm)
	\|\| out.back().type == pid::Higgs
	\|\| has_unof_gluon(in, out)
	);
	return exp(exponent);
	}
	} // namespace RHEJ

	namespace {
	//! Lipatov vertex for partons emitted into extremal jets
	double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
	{
	CLHEP::HepLorentzVector temptrans=-(qav+qbv);
	CLHEP::HepLorentzVector p5=qav-qbv;
	CLHEP::HepLorentzVector CL=temptrans+(qav.m2()/p5.dot(p1)+2.p5.dot(p2)/p1.dot(p2))p1-p2(qbv.m2()/p5.dot(p2)+2.p5.dot(p1)/p1.dot(p2));

	#if printoutput
	cout << "#Fadin qa : "<<qav<<endl;
	cout << "#Fadin qb : "<<qbv<<endl;
	cout << "#Fadin p1 : "<<p1<<endl;
	cout << "#Fadin p2 : "<<p2<<endl;
	cout << "#Fadin p5 : "<<p5<<endl;
	cout << "#Fadin Gauge Check : "<< CL.dot(p5)<<endl;
	cout << "#Fadin C2L : "<< -CL.dot(CL)<<" "<<-CL.dot(CL)/(qav.m2()*qbv.m2())/(4./p5.perp2())<<endl;
	#endif
	#if 0
	if (-CL.dot(CL)<0.)
	// if (fabs(CL.dot(p5))>fabs(CL.dot(CL))) // not sufficient!
	return 0.;
	else
	#endif
	return -CL.dot(CL);
	}

	//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
	double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
	{
	double kperp=(qav-qbv).perp();
	if (kperp>RHEJ::CLAMBDA)
	return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
	else {
	double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
	double temp=Cls-4./(kperp*kperp);

	// std::cout <<kperp <<" "<<temp<<" "<<4./(kperpkperp)<<" "<<(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()qbv.m2()))<<std::endl;
	return temp;
	}
	}

	//! Lipatov vertex
	double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
	{
	CLHEP::HepLorentzVector temptrans=-(qav+qbv);
	CLHEP::HepLorentzVector p5=qav-qbv;
	CLHEP::HepLorentzVector CL=temptrans+qav.m2()(1./p5.dot(pip)pip+1./p5.dot(pop)pop)/2.-qbv.m2()(1./p5.dot(pim)pim+1./p5.dot(pom)pom)/2.+(pip(p5.dot(pim)/pip.dot(pim)+p5.dot(pom)/pip.dot(pom)) + pop(p5.dot(pim)/pop.dot(pim)+p5.dot(pom)/pop.dot(pom)) - pim(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim)) - pom(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;

	return -CL.dot(CL);
	}

	//! Lipatov vertex with soft subtraction
	double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
	{
	double kperp=(qav-qbv).perp();
	if (kperp>RHEJ::CLAMBDA)
	return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
	else {
	double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
	double temp=Cls-4./(kperp*kperp);
	return temp;
	}
	}

	/** Matrix element squared for tree-level current-current scattering
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @returns ME Squared for Tree-Level Current-Current Scattering
	*/
	double ME_current(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa
	){
	if (aptype==21&&bptype==21) {
	return jM2gg(pn,pb,p1,pa);
	} else if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return jM2qg(pn,pb,p1,pa);
	else
	return jM2qbarg(pn,pb,p1,pa);
	}
	else if (bptype==21&&aptype!=21) { // ----- \|\| -----
	if (aptype > 0)
	return jM2qg(p1,pa,pn,pb);
	else
	return jM2qbarg(p1,pa,pn,pb);
	}
	else { // they are both quark
	if (bptype>0) {
	if (aptype>0)
	return jM2qQ(pn,pb,p1,pa);
	else
	return jM2qQbar(pn,pb,p1,pa);
	}
	else {
	if (aptype>0)
	return jM2qQbar(p1,pa,pn,pb);
	else
	return jM2qbarQbar(pn,pb,p1,pa);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** \brief Matrix element squared for tree-level current-current scattering with Higgs
	* @param aptype Particle a PDG ID
	* @param bptype Particle b PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param qH t-channel momentum before Higgs
	* @param qHp1 t-channel momentum after Higgs
	* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
	*/
	double ME_Higgs_current(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
	CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
	double mt, bool include_bottom, double mb
	){
	if (aptype==21&&bptype==21) // gg initial state
	return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
	else
	return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
	}
	else if (bptype==21&&aptype!=21) {
	if (aptype > 0)
	return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
	else
	return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
	}
	else { // they are both quark
	if (bptype>0) {
	if (aptype>0)
	return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)4.4./(9.*9.);
	else
	return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)4.4./(9.*9.);
	}
	else {
	if (aptype>0)
	return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)4.4./(9.*9.);
	else
	return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)4.4./(9.*9.);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** \brief Current matrix element squared with Higgs and unordered forward emission
	* @param aptype Particle A PDG ID
	* @param bptype Particle B PDG ID
	* @param punof Unordered Particle Momentum
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param qH t-channel momentum before Higgs
	* @param qHp1 t-channel momentum after Higgs
	* @returns ME Squared with Higgs and unordered forward emission
	*/
	double ME_Higgs_current_unof(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & punof,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
	CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
	double mt, bool include_bottom, double mb
	){
	if (aptype==21&&bptype!=21) {
	if (bptype > 0)
	return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	else { // they are both quark
	if (bptype>0) {
	if (aptype>0)
	return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	else {
	if (aptype>0)
	return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	/** \brief Current matrix element squared with Higgs and unordered backward emission
	* @param aptype Particle A PDG ID
	* @param bptype Particle B PDG ID
	* @param pn Particle n Momentum
	* @param pb Particle b Momentum
	* @param punob Unordered back Particle Momentum
	* @param p1 Particle 1 Momentum
	* @param pa Particle a Momentum
	* @param qH t-channel momentum before Higgs
	* @param qHp1 t-channel momentum after Higgs
	* @returns ME Squared with Higgs and unordered backward emission
	*/
	double ME_Higgs_current_unob(
	int aptype, int bptype,
	CLHEP::HepLorentzVector const & pn,
	CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & punob,
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & pa,
	CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
	CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
	double mt, bool include_bottom, double mb
	){
	if (bptype==21&&aptype!=21) {
	if (aptype > 0)
	return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	else { // they are both quark
	if (aptype>0) {
	if (bptype>0)
	return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	else {
	if (bptype>0)
	return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	else
	return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
	}
	}
	throw std::logic_error("unknown particle types");
	}

	CLHEP::HepLorentzVector to_HepLorentzVector(RHEJ::Sparticle const & particle){
	return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
	}
	} // namespace anonymous

	namespace RHEJ{
	MatrixElement::MatrixElement(
	std::function<double (double)> alpha_s,
	- Config const & conf
	+ MatrixElementConfig conf
	):
	alpha_s_{std::move(alpha_s)},
	- param_(conf)
	+ param_{std::move(conf)}
	{}

	double MatrixElement::operator()(
	double mur,
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const {
	return tree(
	mur,
	incoming, outgoing,
	check_momenta
	)*virtual_corrections(
	mur,
	incoming, outgoing
	);
	}

	double MatrixElement::tree_kin(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const {
	assert(
	std::is_sorted(
	incoming.begin(), incoming.end(),
	[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
	)
	);
	assert(std::is_sorted(outgoing.begin(), outgoing.end(), rapidity_less{}));

	auto AWZH_boson = std::find_if(
	begin(outgoing), end(outgoing),
	[](Sparticle const & p){return is_AWZH_boson(p);}
	);

	if(AWZH_boson == end(outgoing)){
	return tree_kin_jets(incoming, outgoing, check_momenta);
	}

	switch(AWZH_boson->type){
	case pid::Higgs: {
	static constexpr double mH = 125.;
	const double alpha_s_mH = alpha_s_(mH);
	return alpha_s_mHalpha_s_mH/(256.pow(M_PI, 5))*tree_kin_Higgs(
	incoming, outgoing, check_momenta
	);
	}
	// TODO
	case pid::Wp:
	case pid::Wm:
	case pid::photon:
	case pid::Z:
	default:
	throw std::logic_error("Emission of boson of unsupported type");
	}
	}

	namespace{
	constexpr int extremal_jet_idx = 1;
	constexpr int no_extremal_jet_idx = 0;

	bool treat_as_extremal(Sparticle const & parton){
	return parton.p.user_index() == extremal_jet_idx;
	}

	template<class InputIterator>
	double FKL_ladder_weight(
	InputIterator begin_gluon, InputIterator end_gluon,
	CLHEP::HepLorentzVector const & q0,
	CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
	CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn
	){
	double wt = 1;
	auto qi = q0;
	for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
	assert(gluon_it->type == pid::gluon);
	const auto g = to_HepLorentzVector(*gluon_it);
	const auto qip1 = qi - g;

	if(treat_as_extremal(*gluon_it)){
	wt = C2Lipatovots(qip1, qi, pa, pb)C_A;
	} else{
	wt = C2Lipatovots(qip1, qi, pa, pb, p1, pn)C_A;
	}

	qi = qip1;
	}
	return wt;
	}

	} // namespace anonymous

	std::vector<Sparticle> MatrixElement::tag_extremal_jet_partons(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> out_partons, bool check_momenta
	) const{
	if(!check_momenta){
	for(auto & parton: out_partons){
	parton.p.set_user_index(no_extremal_jet_idx);
	}
	return out_partons;
	}
	- fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param_.jet_def());
	- const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_ptmin()));
	+ fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param_.jet_param.def);
	+ const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
	assert(jets.size() >= 2);
	auto most_backward = begin(jets);
	auto most_forward = end(jets) - 1;
	// skip jets caused by unordered emission
	if(has_unob_gluon(incoming, out_partons)){
	assert(jets.size() >= 3);
	++most_backward;
	}
	else if(has_unof_gluon(incoming, out_partons)){
	assert(jets.size() >= 3);
	--most_forward;
	}
	const auto extremal_jet_indices = cs.particle_jet_indices(
	{most_backward, most_forward}
	);
	assert(extremal_jet_indices.size() == out_partons.size());
	for(size_t i = 0; i < out_partons.size(); ++i){
	assert(RHEJ::is_parton(out_partons[i]));
	const int idx = (extremal_jet_indices[i]>=0)?
	extremal_jet_idx:
	no_extremal_jet_idx;
	out_partons[i].p.set_user_index(idx);
	}
	return out_partons;
	}

	double MatrixElement::tree_kin_jets(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> partons,
	bool check_momenta
	) const {
	partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
	if(has_unob_gluon(incoming, partons) \|\| has_unof_gluon(incoming, partons)){
	throw std::logic_error("unordered emission not implemented for pure jets");
	}

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	const auto p1 = to_HepLorentzVector(partons.front());
	const auto pn = to_HepLorentzVector(partons.back());

	return ME_current(
	incoming[0].type, incoming[1].type,
	pn, pb, p1, pa
	)/(4(N_CN_C - 1))*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	pa - p1, pa, pb, p1, pn
	);
	}

	double MatrixElement::tree_kin_Higgs(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const {
	if(has_uno_gluon(incoming, outgoing)){
	return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
	}
	if(outgoing.front().type == pid::Higgs){
	return tree_kin_Higgs_first(incoming, outgoing, check_momenta);
	}
	if(outgoing.back().type == pid::Higgs){
	return tree_kin_Higgs_last(incoming, outgoing, check_momenta);
	}
	return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
	}

	double MatrixElement::MH2_forwardH(
	RHEJ::ParticleID id,
	CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
	CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
	CLHEP::HepLorentzVector pH,
	double t1, double t2
	) const{
	ignore(p2out, p2in);
	const double shat = p1in.invariantMass2(p2in);
	assert(RHEJ::is_parton(id));
	if(id != RHEJ::pid::gluon){
	return 9./2.shatshatC2qHqm(p1in,p1out,pH)/(t1t2);
	}
	// gluon case
	#ifdef RHEJ_BUILD_WITH_QCDLOOP
	- if(!param_.Hgg_setting().use_impact_factors){
	+ if(!param_.Higgs_coupling.use_impact_factors){
	return C_A/C_F1./(16M_PIM_PI)t1/t2*MH2gq_outsideH(
	p1out, p1in, p2out, p2in, pH,
	- param_.Hgg_setting().mt, param_.Hgg_setting().include_bottom,
	- param_.Hgg_setting().mb
	+ param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
	+ param_.Higgs_coupling.mb
	);
	}
	#endif
	return 9./2.shatshat*(
	C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
	)/(t1*t2);
	}

	double MatrixElement::tree_kin_Higgs_first(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const {
	assert(outgoing.front().type == pid::Higgs);
	const auto pH = to_HepLorentzVector(outgoing.front());
	const auto partons = tag_extremal_jet_partons(
	incoming,
	std::vector<Sparticle>(begin(outgoing) + 1, end(outgoing)),
	check_momenta
	);

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	const auto p1 = to_HepLorentzVector(partons.front());
	const auto pn = to_HepLorentzVector(partons.back());

	const auto q0 = pa - p1 - pH;

	const double t1 = q0.m2();
	const double t2 = (pn - pb).m2();

	double wt = MH2_forwardH(
	incoming[0].type, p1, pa, pn, pb, pH,
	t1, t2
	)*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	q0, pa, pb, p1, pn
	);

	for(auto const & inc: incoming){
	if(inc.type != pid::gluon) wt *= C_F/C_A;
	}
	return wt;
	}

	double MatrixElement::tree_kin_Higgs_last(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const {
	assert(outgoing.back().type == pid::Higgs);
	const auto pH = to_HepLorentzVector(outgoing.back());
	const auto partons = tag_extremal_jet_partons(
	incoming,
	std::vector<Sparticle>(begin(outgoing), end(outgoing) - 1),
	check_momenta
	);

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	auto p1 = to_HepLorentzVector(partons.front());
	const auto pn = to_HepLorentzVector(partons.back());

	auto q0 = pa - p1;

	const double t1 = q0.m2();
	const double t2 = (pn + pH - pb).m2();

	double wt = MH2_forwardH(
	incoming[1].type, pn, pb, p1, pa, pH,
	t2, t1
	)*FKL_ladder_weight(
	begin(partons) + 1, end(partons) - 1,
	q0, pa, pb, p1, pn
	);

	for(auto const & inc: incoming){
	if(inc.type != pid::gluon) wt *= C_F/C_A;
	}
	return wt;
	}


	double MatrixElement::tree_kin_Higgs_between(
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const {
	const auto the_Higgs = std::find_if(
	begin(outgoing), end(outgoing),
	[](Sparticle const & s){ return s.type == pid::Higgs; }
	);
	assert(the_Higgs != end(outgoing));
	const auto pH = to_HepLorentzVector(*the_Higgs);
	std::vector<Sparticle> partons(begin(outgoing), the_Higgs);
	partons.insert(end(partons), the_Higgs + 1, end(outgoing));
	partons = tag_extremal_jet_partons(incoming, partons, check_momenta);

	const auto pa = to_HepLorentzVector(incoming[0]);
	const auto pb = to_HepLorentzVector(incoming[1]);

	auto p1 = to_HepLorentzVector(
	partons[has_unob_gluon(incoming, outgoing)?1:0]
	);
	auto pn = to_HepLorentzVector(
	partons[partons.size() - (has_unof_gluon(incoming, outgoing)?2:1)]
	);

	auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
	assert(
	(first_after_Higgs == end(partons) && has_unob_gluon(incoming, outgoing))
	\|\| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
	);
	assert(
	(first_after_Higgs == begin(partons) && has_unof_gluon(incoming, outgoing))
	\|\| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
	);
	// always treat the Higgs as if it were in between the extremal FKL partons
	if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
	else if(first_after_Higgs == end(partons)) --first_after_Higgs;

	// t-channel momentum before Higgs
	auto qH = pa;
	for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
	qH -= to_HepLorentzVector(*parton_it);
	}

	auto q0 = pa - p1;
	auto begin_ladder = begin(partons) + 1;
	auto end_ladder = end(partons) - 1;

	double current_factor;
	if(has_unob_gluon(incoming, outgoing)){
	current_factor = 9./2.*ME_Higgs_current_unob(
	incoming[0].type, incoming[1].type,
	pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH,
	- param_.Hgg_setting().mt, param_.Hgg_setting().include_bottom, param_.Hgg_setting().mb
	+ param_.Higgs_coupling.mt,
	+ param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
	);
	const auto p_unob = to_HepLorentzVector(partons.front());
	q0 -= p_unob;
	p1 += p_unob;
	++begin_ladder;
	}
	else if(has_unof_gluon(incoming, outgoing)){
	current_factor = 9./2.*ME_Higgs_current_unof(
	incoming[0].type, incoming[1].type,
	to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH,
	- param_.Hgg_setting().mt, param_.Hgg_setting().include_bottom, param_.Hgg_setting().mb
	+ param_.Higgs_coupling.mt,
	+ param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
	);
	pn += to_HepLorentzVector(partons.back());
	--end_ladder;
	}
	else{
	current_factor = ME_Higgs_current(
	incoming[0].type, incoming[1].type,
	pn, pb, p1, pa, qH, qH - pH,
	- param_.Hgg_setting().mt, param_.Hgg_setting().include_bottom, param_.Hgg_setting().mb
	+ param_.Higgs_coupling.mt,
	+ param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
	);
	}

	const double ladder_factor = FKL_ladder_weight(
	begin_ladder, first_after_Higgs,
	q0, pa, pb, p1, pn
	)*FKL_ladder_weight(
	first_after_Higgs, end_ladder,
	qH - pH, pa, pb, p1, pn
	);
	return current_factor9./8.ladder_factor;
	}

	double MatrixElement::tree_param_partons(
	double alpha_s, double mur,
	std::vector<Sparticle> const & partons
	) const{
	const double gs2 = 4.M_PIalpha_s;
	double wt = std::pow(gs2, partons.size());
	- if(param_.log_corr()){
	+ if(param_.log_correction){
	// use alpha_s(q_perp), evolved to mur
	assert(partons.size() >= 2);
	for(size_t i = 1; i < partons.size()-1; ++i){
	wt = 1 + alpha_s/(2M_PI)beta0log(mur/partons[i].p.perp());
	}
	}
	return wt;
	}

	double MatrixElement::tree_param(
	double mur,
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing
	) const{
	const double alpha_s = alpha_s_(mur);
	if(has_unob_gluon(incoming, outgoing)){
	return 4M_PIalpha_s*tree_param_partons(
	alpha_s, mur, filter_partons({begin(outgoing) + 1, end(outgoing)})
	);
	}
	if(has_unof_gluon(incoming, outgoing)){
	return 4M_PIalpha_s*tree_param_partons(
	alpha_s, mur, filter_partons({begin(outgoing), end(outgoing) - 1})
	);
	}
	return tree_param_partons(alpha_s, mur, filter_partons(outgoing));
	}

	double MatrixElement::tree(
	double mur,
	std::array<Sparticle, 2> const & incoming,
	std::vector<Sparticle> const & outgoing,
	bool check_momenta
	) const {
	return tree_param(mur, incoming, outgoing)*tree_kin(
	incoming, outgoing, check_momenta
	);
	}
	} // namespace RHEJ
	diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
	index 160ef6f..46bd4bc 100644
	--- a/src/PhaseSpacePoint.cc
	+++ b/src/PhaseSpacePoint.cc
	@@ -1,579 +1,579 @@
	#include "RHEJ/PhaseSpacePoint.hh"

	#include <random>

	#include "RHEJ/Constants.hh"
	#include "RHEJ/resummation_jet_momenta.hh"
	#include "RHEJ/Jacobian.hh"
	#include "RHEJ/RNGWrapper.hh"
	#include "RHEJ/uno.hh"
	#include "RHEJ/debug.hh"
	#include "RHEJ/kinematics.hh"

	namespace RHEJ{

	namespace {
	//generate Ranlux64Engine with fixed, predefined state
	/*
	* some (all?) of the Ranlux64Engine constructors leave fields
	* uninitialised, invoking undefined behaviour. This can be
	* circumvented by restoring the state from a file
	*/
	CLHEP::Ranlux64Engine gen_Ranlux64Engine(){
	static const std::string state =
	"9876\n"
	"0.91280703978419097666\n"
	"0.41606065829518357191\n"
	"0.99156342622341142601\n"
	"0.030922955274050423213\n"
	"0.16206278421638486975\n"
	"0.76151768001958330956\n"
	"0.43765760066092695979\n"
	"0.42904698253748563275\n"
	"0.11476317525663759511\n"
	"0.026620053590963976831\n"
	"0.65953715764414511114\n"
	"0.30136722624439826745\n"
	"3.5527136788005009294e-15 4\n"
	"1 202\n";
	const std::string file = std::tmpnam(nullptr);
	{
	std::ofstream out{file};
	out << state;
	}
	CLHEP::Ranlux64Engine result;
	result.restoreStatus(file.c_str());
	return result;
	}
	}

	CLHEP::Ranlux64Engine PhaseSpacePoint::ran_{gen_Ranlux64Engine()};


	void PhaseSpacePoint::reset_ranlux(std::string const & init_file){
	reset_ranlux(init_file.c_str());
	}

	void PhaseSpacePoint::reset_ranlux(char const * init_file){
	ran_.restoreStatus(init_file);
	}

	namespace {
	constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;

	bool is_nonjet_parton(fastjet::PseudoJet const & parton){
	assert(parton.user_index() != -1);
	return parton.user_index() > max_jet_user_idx;
	}

	bool is_jet_parton(fastjet::PseudoJet const & parton){
	assert(parton.user_index() != -1);
	return parton.user_index() <= max_jet_user_idx;
	}

	// user indices for partons with extremal rapidity
	constexpr int unob_idx = -5;
	constexpr int unof_idx = -4;
	constexpr int backward_FKL_idx = -3;
	constexpr int forward_FKL_idx = -2;
	}

	namespace {
	double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
	const double delta_y =
	Born_jets.back().rapidity() - Born_jets.front().rapidity();
	assert(delta_y > 0);
	// Formula derived from fit in reversed HEJ intro paper
	return 0.975052*delta_y;
	}
	}

	std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
	std::vector<fastjet::PseudoJet> const & partons
	) const{
	- fastjet::ClusterSequence cs(partons, param_.jet_def());
	- return cs.inclusive_jets(param_.jet_ptmin());
	+ fastjet::ClusterSequence cs(partons, param_.jet_param.def);
	+ return cs.inclusive_jets(param_.jet_param.min_pt);
	}

	bool PhaseSpacePoint::pass_resummation_cuts(
	std::vector<fastjet::PseudoJet> const & jets
	) const{
	return cluster_jets(jets).size() == jets.size();
	}

	int PhaseSpacePoint::sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets){
	const double ng_mean = estimate_ng_mean(Born_jets);
	std::poisson_distribution<int> dist(ng_mean);
	RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
	const int ng = dist(rng);
	assert(ng >= 0);
	assert(ng < ng_max);
	weight_ = std::tgamma(ng + 1)std::exp(ng_mean)*std::pow(ng_mean, -ng);
	return ng;
	}

	void PhaseSpacePoint::copy_AWZH_boson_from(Event const & event){
	auto const & from = event.outgoing();

	const auto AWZH_boson = std::find_if(
	begin(from), end(from),
	[](Sparticle const & p){ return is_AWZH_boson(p); }
	);
	if(AWZH_boson == end(from)) return;
	auto insertion_point = std::lower_bound(
	begin(outgoing_), end(outgoing_), *AWZH_boson, rapidity_less{}
	);
	outgoing_.insert(insertion_point, *AWZH_boson);

	// copy decay products
	const int idx = std::distance(begin(from), AWZH_boson);
	const auto decay_it = event.decays().find(idx);
	if(decay_it != end(event.decays())){
	const int new_idx = std::distance(begin(outgoing_), insertion_point);
	assert(outgoing_[new_idx].type == AWZH_boson->type);
	decays_.emplace(new_idx, decay_it->second);
	}

	assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
	}

	- PhaseSpacePoint::PhaseSpacePoint(Event const & ev, Config const & conf):
	+ PhaseSpacePoint::PhaseSpacePoint(Event const & ev, PhaseSpacePointConfig conf):
	unob_{has_unob_gluon(ev.incoming(), ev.outgoing())},
	unof_{!unob_ && has_unof_gluon(ev.incoming(), ev.outgoing())},
	- param_(conf),
	- splitter_{param_.jet_R(), param_.jet_def(), param_.jet_ptmin(), ran_}
	+ param_{std::move(conf)},
	+ splitter_{param_.jet_param.def.R(), param_.jet_param.def, param_.jet_param.min_pt, ran_}
	{
	weight_ = 1;
	const auto Born_jets = sorted_by_rapidity(ev.jets());
	const int ng = sample_ng(Born_jets);
	weight_ /= std::tgamma(ng + 1);
	const int ng_jets = sample_ng_jets(ng, Born_jets);
	std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
	- ng - ng_jets, CMINPT, param_.jet_ptmin()
	+ ng - ng_jets, CMINPT, param_.jet_param.min_pt
	);
	{
	const auto qperp = std::accumulate(
	begin(out_partons), end(out_partons),
	fastjet::PseudoJet{}
	);
	const auto jets = reshuffle(Born_jets, qperp);
	if(weight_ == 0.) return;
	if(! pass_resummation_cuts(jets)){
	weight_ = 0.;
	return;
	}
	std::vector<fastjet::PseudoJet> jet_partons = split(jets, ng_jets);
	if(weight_ == 0.) return;
	rescale_rapidities(
	out_partons,
	most_backward_FKL(jet_partons).rapidity(),
	most_forward_FKL(jet_partons).rapidity()
	);
	if(! cluster_jets(out_partons).empty()){
	weight_ = 0.;
	return;
	}
	std::sort(begin(out_partons), end(out_partons), rapidity_less{});
	assert(
	std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{})
	);
	const auto first_jet_parton = out_partons.insert(
	end(out_partons), begin(jet_partons), end(jet_partons)
	);
	std::inplace_merge(
	begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{}
	);
	}
	if(! jets_ok(Born_jets, out_partons)){
	weight_ = 0.;
	return;
	}
	weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size());

	outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H
	for(auto & p: out_partons){
	outgoing_.emplace_back(Sparticle{pid::gluon, std::move(p)});
	}
	most_backward_FKL(outgoing_).type = ev.incoming().front().type;
	most_forward_FKL(outgoing_).type = ev.incoming().back().type;
	copy_AWZH_boson_from(ev);

	assert(!outgoing_.empty());
	reconstruct_incoming(ev.incoming());
	}

	std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
	int count, double ptmin, double ptmax
	){
	// heuristic parameters for pt sampling
	const double ptpar = 1.3 + count/5.;
	const double temp1 = atan((ptmax - ptmin)/ptpar);

	std::vector<fastjet::PseudoJet> partons(count);
	for(size_t i = 0; i < (size_t) count; ++i){
	const double r1 = ran_.flat();
	const double pt = ptmin + ptpartan(r1temp1);
	const double temp2 = cos(r1*temp1);
	const double phi = 2M_PIran_.flat();
	weight_ = 2.0M_PIptptpartemp1/(temp2temp2);
	// we don't know the allowed rapidity span yet,
	// set a random value to be rescaled later on
	const double y = ran_.flat();
	partons[i].reset_PtYPhiM(pt, y, phi);
	// Set user index higher than any jet-parton index
	// in order to assert that these are not inside jets
	partons[i].set_user_index(i + 1 + RHEJ::PhaseSpacePoint::ng_max);

	assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5);
	}
	assert(std::all_of(partons.cbegin(), partons.cend(), is_nonjet_parton));
	return partons;
	}

	void PhaseSpacePoint::rescale_rapidities(
	std::vector<fastjet::PseudoJet> & partons,
	double ymin, double ymax
	){
	constexpr double ep = 1e-7;
	for(auto & parton: partons){
	assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
	const double dy = ymax - ymin - 2*ep;
	const double y = ymin + ep + dy*parton.rapidity();
	parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
	weight_ *= dy;
	assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
	}
	}

	namespace {
	template<typename T, typename... Rest>
	auto min(T const & a, T const & b, Rest&&... r) {
	using std::min;
	return min(a, min(b, std::forward<Rest>(r)...));
	}
	}

	double PhaseSpacePoint::probability_in_jet(
	std::vector<fastjet::PseudoJet> const & Born_jets
	) const{
	assert(std::is_sorted(begin(Born_jets), end(Born_jets), rapidity_less{}));
	assert(Born_jets.size() >= 2);

	const double dy =
	Born_jets.back().rapidity() - Born_jets.front().rapidity();
	- const double R = param_.jet_R();
	+ const double R = param_.jet_param.def.R();
	const int njets = Born_jets.size();
	const double p_J_y_large = (njets-1)RR/(2.*dy);
	const double p_J_y0 = njets*R/M_PI;
	return min(p_J_y_large, p_J_y0, 1.);
	}

	int PhaseSpacePoint::sample_ng_jets(
	int ng, std::vector<fastjet::PseudoJet> const & Born_jets
	){
	RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
	const double p_J = probability_in_jet(Born_jets);
	std::binomial_distribution<> bin_dist(ng, p_J);
	const int ng_J = bin_dist(rng);
	weight_ = std::pow(p_J, -ng_J)std::pow(1 - p_J, ng_J - ng);
	return ng_J;
	}

	std::vector<fastjet::PseudoJet>
	PhaseSpacePoint::reshuffle(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	fastjet::PseudoJet const & q
	){
	if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
	std::vector<fastjet::PseudoJet> jets = resummation_jet_momenta(Born_jets, q);
	if(jets.empty()){
	weight_ = 0;
	return {};
	}
	// transform delta functions to integration over resummation momenta
	weight_ /= Jacobian(jets, q);
	return jets;
	}

	std::vector<int> PhaseSpacePoint::distribute_jet_partons(
	int ng_jets, std::vector<fastjet::PseudoJet> const & jets
	){
	size_t first_valid_jet = 0;
	size_t num_valid_jets = jets.size();
	- const double R_eff = 5./3.*param_.jet_R();
	+ const double R_eff = 5./3.*param_.jet_param.def.R();
	// if there is an unordered jet too far away from the FKL jets
	// then extra gluon constituents of the unordered jet would
	// violate the FKL rapidity ordering
	if(unob_ && jets[0].delta_R(jets[1]) > R_eff){
	++first_valid_jet;
	--num_valid_jets;
	}
	else if(unof_ && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
	--num_valid_jets;
	}
	std::vector<int> np(jets.size(), 1);
	for(int i = 0; i < ng_jets; ++i){
	++np[first_valid_jet + ran_.flat() * num_valid_jets];
	}
	weight_ *= std::pow(num_valid_jets, ng_jets);
	return np;
	}

	#ifndef NDEBUG
	namespace{
	bool tagged_FKL_backward(
	std::vector<fastjet::PseudoJet> const & jet_partons
	){
	return std::find_if(
	begin(jet_partons), end(jet_partons),
	[](fastjet::PseudoJet const & p){
	return p.user_index() == backward_FKL_idx;
	}
	) != end(jet_partons);
	}

	bool tagged_FKL_forward(
	std::vector<fastjet::PseudoJet> const & jet_partons
	){
	// the most forward FKL parton is most likely near the end of jet_partons;
	// start search from there
	return std::find_if(
	jet_partons.rbegin(), jet_partons.rend(),
	[](fastjet::PseudoJet const & p){
	return p.user_index() == forward_FKL_idx;
	}
	) != jet_partons.rend();
	}

	bool tagged_FKL_extremal(
	std::vector<fastjet::PseudoJet> const & jet_partons
	){
	return tagged_FKL_backward(jet_partons) && tagged_FKL_forward(jet_partons);
	}

	} // namespace anonymous
	#endif

	std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
	std::vector<fastjet::PseudoJet> const & jets,
	int ng_jets
	){
	return split(jets, distribute_jet_partons(ng_jets, jets));
	}

	bool PhaseSpacePoint::pass_extremal_cuts(
	fastjet::PseudoJet const & ext_parton,
	fastjet::PseudoJet const & jet
	) const{
	- if(ext_parton.pt() < param_.min_ext_pt()) return false;
	- return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction();
	+ if(ext_parton.pt() < param_.min_extparton_pt) return false;
	+ return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction;
	}

	std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
	std::vector<fastjet::PseudoJet> const & jets,
	std::vector<int> const & np
	){
	assert(! jets.empty());
	assert(jets.size() == np.size());
	assert(pass_resummation_cuts(jets));

	const size_t most_backward_FKL_idx = 0 + unob_;
	const size_t most_forward_FKL_idx = jets.size() - 1 - unof_;

	std::vector<fastjet::PseudoJet> jet_partons;
	// randomly distribute jet gluons among jets
	for(size_t i = 0; i < jets.size(); ++i){
	weight_ *= splitter_.Split(jets[i], np[i]);
	if(weight_ == 0) return {};
	assert(
	std::all_of(
	begin(splitter_.get_jcons()), end(splitter_.get_jcons()),
	is_jet_parton
	)
	);
	const auto first_new_parton = jet_partons.insert(
	end(jet_partons),
	begin(splitter_.get_jcons()), end(splitter_.get_jcons())
	);
	// mark uno and extremal FKL emissions here so we can check
	// their position once all emissions are generated
	auto extremal = end(jet_partons);
	if((unob_ && i == 0) \|\| i == most_backward_FKL_idx){
	// unordered or FKL backward emission
	extremal = std::min_element(
	first_new_parton, end(jet_partons), rapidity_less{}
	);
	extremal->set_user_index(
	(i == most_backward_FKL_idx)?backward_FKL_idx:unob_idx
	);
	}
	else if((unof_ && i == jets.size() - 1) \|\| i == most_forward_FKL_idx){
	// unordered or FKL forward emission
	extremal = std::max_element(
	first_new_parton, end(jet_partons), rapidity_less{}
	);
	extremal->set_user_index(
	(i == most_forward_FKL_idx)?forward_FKL_idx:unof_idx
	);
	}
	if(
	extremal != end(jet_partons)
	&& !pass_extremal_cuts(*extremal, jets[i])
	){
	weight_ = 0;
	return {};
	}
	}
	assert(tagged_FKL_extremal(jet_partons));
	std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
	if(
	!extremal_ok(jet_partons)
	\|\| !split_preserved_jets(jets, jet_partons)
	){
	weight_ = 0.;
	return {};
	}
	return jet_partons;
	}

	bool PhaseSpacePoint::extremal_ok(
	std::vector<fastjet::PseudoJet> const & partons
	) const{
	assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
	if(unob_ && partons.front().user_index() != unob_idx) return false;
	if(unof_ && partons.back().user_index() != unof_idx) return false;
	return
	most_backward_FKL(partons).user_index() == backward_FKL_idx
	&& most_forward_FKL(partons).user_index() == forward_FKL_idx;
	}

	bool PhaseSpacePoint::split_preserved_jets(
	std::vector<fastjet::PseudoJet> const & jets,
	std::vector<fastjet::PseudoJet> const & jet_partons
	) const{
	assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));

	const auto split_jets = sorted_by_rapidity(cluster_jets(jet_partons));
	// this can happen if two overlapping jets
	// are both split into more than one parton
	if(split_jets.size() != jets.size()) return false;
	for(size_t i = 0; i < split_jets.size(); ++i){
	// this can happen if there are two overlapping jets
	// and a parton is assigned to the "wrong" jet
	if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
	return false;
	}
	}
	return true;
	}

	template<class Particle>
	Particle const & PhaseSpacePoint::most_backward_FKL(
	std::vector<Particle> const & partons
	) const{
	return partons[0 + unob_];
	}

	template<class Particle>
	Particle const & PhaseSpacePoint::most_forward_FKL(
	std::vector<Particle> const & partons
	) const{
	const size_t idx = partons.size() - 1 - unof_;
	assert(idx < partons.size());
	return partons[idx];
	}

	template<class Particle>
	Particle & PhaseSpacePoint::most_backward_FKL(
	std::vector<Particle> & partons
	) const{
	return partons[0 + unob_];
	}

	template<class Particle>
	Particle & PhaseSpacePoint::most_forward_FKL(
	std::vector<Particle> & partons
	) const{
	const size_t idx = partons.size() - 1 - unof_;
	assert(idx < partons.size());
	return partons[idx];
	}

	namespace {
	bool contains_idx(
	fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
	){
	auto const & constituents = jet.constituents();
	const int idx = parton.user_index();
	return std::find_if(
	begin(constituents), end(constituents),
	[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
	) != end(constituents);
	}
	}

	/**
	* 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 PhaseSpacePoint::jets_ok(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	std::vector<fastjet::PseudoJet> const & partons
	) const{
	- fastjet::ClusterSequence cs(partons, param_.jet_def());
	- const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_ptmin()));
	+ fastjet::ClusterSequence cs(partons, param_.jet_param.def);
	+ const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
	if(jets.size() != Born_jets.size()) return false;
	int in_jet = 0;
	for(size_t i = 0; i < jets.size(); ++i){
	assert(jets[i].has_constituents());
	for(auto && parton: jets[i].constituents()){
	if(is_nonjet_parton(parton)) return false;
	}
	in_jet += jets[i].constituents().size();
	}
	const int expect_in_jet = std::count_if(
	partons.cbegin(), partons.cend(), is_jet_parton
	);
	if(in_jet != expect_in_jet) return false;
	// note that PseudoJet::contains does not work here
	if(! (
	contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
	&& contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
	)) return false;
	if(unob_ && !contains_idx(jets.front(), partons.front())) return false;
	if(unof_ && !contains_idx(jets.back(), partons.back())) return false;
	for(size_t i = 0; i < jets.size(); ++i){
	assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
	}
	return true;
	}

	void PhaseSpacePoint::reconstruct_incoming(
	std::array<Sparticle, 2> const & Born_incoming
	){
	std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
	for(size_t i = 0; i < incoming_.size(); ++i){
	incoming_[i].type = Born_incoming[i].type;
	}
	assert(momentum_conserved());
	}

	double PhaseSpacePoint::phase_space_normalisation(
	int num_Born_jets, int num_out_partons
	) const{
	return pow(16*pow(M_PI,3), num_Born_jets - num_out_partons);
	}

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

	} //namespace RHEJ
	diff --git a/src/main.cc b/src/main.cc
	index 4f034da..e2dfafd 100644
	--- a/src/main.cc
	+++ b/src/main.cc
	@@ -1,160 +1,160 @@
	/**
	* Name: main.cc
	* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
	*
	*/
	#include <fstream>
	#include <algorithm>
	#include <memory>
	#include <chrono>
	#include <iostream>

	#include "yaml-cpp/yaml.h"

	#include "RHEJ/CombinedEventWriter.hh"
	#include "RHEJ/config.hh"
	#include "RHEJ/EventReweighter.hh"
	#include "RHEJ/get_analysis.hh"
	#include "LHEF/LHEF.h"
	#include "RHEJ/utility.hh"
	#include "RHEJ/Version.hh"
	#include "RHEJ/stream.hh"
	#include "RHEJ/YAMLreader.hh"

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

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

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

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

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

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

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

	// read configuration
	const RHEJ::Config config = load_config(argv[1]);

	RHEJ::istream in{argv[2]};

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

	LHEF::Reader reader{in};
	- RHEJ::CombinedEventWriter writer{config, reader.heprup};
	+ RHEJ::CombinedEventWriter writer{config.output, reader.heprup};

	double global_reweight = 1.;
	int max_events = std::numeric_limits<int>::max();
	if(argn > 3){
	max_events = std::stoi(argv[3]);
	const int input_events = event_number(reader.headerBlock);
	global_reweight = input_events/static_cast<double>(max_events);
	std::cout << "Processing " << max_events
	<< " out of " << input_events << " events\n";
	}
	RHEJ::EventReweighter rhej{
	reader.heprup,
	- config
	+ to_EventReweighterConfig(config)
	};
	if(config.RanLux_init){
	RHEJ::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
	}

	int nevent = 0;
	std::array<int, RHEJ::event_type::last_type + 1>
	nevent_type{0}, nfailed_type{0};

	// Loop over the events in the inputfile
	while(reader.readEvent()){
	// reweight events so that the total cross section is conserved
	reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());

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

	if (nevent % 10000 == 0){
	std::cout << "Passed " << nevent << " events ("
	<< time_to_string(clock::to_time_t(clock::now())) << ")"<< std::endl;
	}

	// calculate rHEJ weight
	RHEJ::Event FO_event{
	RHEJ::UnclusteredEvent{reader.hepeup},
	config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
	};
	auto resummed_events = rhej.reweight(
	FO_event,
	config.trials, config.scale_gen
	);
	++nevent_type[FO_event.type()];

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

	for(auto const & ev: resummed_events){
	//TODO: move pass_cuts to after phase space point generation
	if(analysis->pass_cuts(ev, FO_event)){
	analysis->fill(ev, FO_event);
	writer.write(ev);
	}
	}
	} // main event loop

	using namespace RHEJ::event_type;
	std::cout<< "Events processed: " << nevent << '\n';
	for(size_t ev_type = first_type; ev_type <= last_type; ++ev_type){
	std::cout << '\t' << names[ev_type] << ": " << nevent_type[ev_type]
	<< ", failed to reconstruct " << nfailed_type[ev_type]
	<< '\n';
	}
	std::chrono::duration<double> run_time = (clock::now() - start_time);
	std::cout << "Finished " << RHEJ::Version::package_name() << " at "
	<< time_to_string(clock::to_time_t(clock::now()))
	<< "\n=> Runtime: " << run_time.count() << " sec ("
	<< nevent/run_time.count() << " Events/sec).\n";

	}
	diff --git a/src/make_writer.cc b/src/make_writer.cc
	index 3235dad..430863f 100644
	--- a/src/make_writer.cc
	+++ b/src/make_writer.cc
	@@ -1,26 +1,25 @@
	#include "RHEJ/make_writer.hh"

	#include "RHEJ/LesHouchesWriter.hh"
	#include "RHEJ/HepMCWriter.hh"

	namespace RHEJ{
	std::unique_ptr<EventWriter> make_format_writer(
	FileFormat format, std::string const & outfile,
	- Config const & conf,
	LHEF::HEPRUP const & heprup
	){
	switch(format){
	case Les_Houches:
	return std::unique_ptr<EventWriter>{
	- new LesHouchesWriter{outfile, conf, heprup}
	+ new LesHouchesWriter{outfile, heprup}
	};
	case HepMC:
	return std::unique_ptr<EventWriter>{
	- new HepMCWriter{outfile, conf, heprup}
	+ new HepMCWriter{outfile, heprup}
	};
	default:
	throw std::logic_error("unhandled file format");
	}
	}

	}
	diff --git a/t/check_res.cc b/t/check_res.cc
	index c8f6c8a..9bfab88 100644
	--- a/t/check_res.cc
	+++ b/t/check_res.cc
	@@ -1,72 +1,78 @@
	#include <iostream>

	#include "LHEF/LHEF.h"
	#include "RHEJ/stream.hh"
	#include "RHEJ/EventReweighter.hh"

	namespace{
	const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
	const fastjet::JetDefinition Born_jet_def{jet_def};
	constexpr double Born_jetptmin = 30;
	constexpr double extpartonptmin = 30;
	constexpr double max_ext_soft_pt_fraction =
	std::numeric_limits<double>::infinity();
	constexpr double jetptmin = 35;
	constexpr bool log_corr = false;
	using EventTreatment = RHEJ::EventTreatment;
	using namespace RHEJ::event_type;
	RHEJ::EventTreatMap treat{
	{no_2_jets, EventTreatment::discard},
	{bad_final_state, EventTreatment::discard},
	{nonFKL, EventTreatment::discard},
	{unof, EventTreatment::discard},
	{unob, EventTreatment::discard},
	{FKL, EventTreatment::reweight}
	};
	};

	int main(int argn, char** argv) {
	if(argn == 5 && std::string(argv[4]) == "uno"){
	--argn;
	treat[unof] = EventTreatment::reweight;
	treat[unob] = EventTreatment::reweight;
	treat[FKL] = EventTreatment::discard;
	}
	if(argn != 4){
	std::cerr << "Usage: check_res eventfile xsection tolerance [uno]";
	return EXIT_FAILURE;
	}

	const double xsec_ref = std::stod(argv[2]);
	const double tolerance = std::stod(argv[3]);

	RHEJ::istream in{argv[1]};
	LHEF::Reader reader{in};

	- RHEJ::Config conf;
	- conf.resummation_jets = RHEJ::JetParameters{jet_def, jetptmin};
	+ RHEJ::PhaseSpacePointConfig psp_conf;
	+ psp_conf.jet_param = RHEJ::JetParameters{jet_def, jetptmin};
	+ psp_conf.min_extparton_pt = extpartonptmin;
	+ psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
	+ RHEJ::MatrixElementConfig ME_conf;
	+ ME_conf.jet_param = psp_conf.jet_param;
	+ ME_conf.log_correction = log_corr;
	+ ME_conf.Higgs_coupling = RHEJ::HiggsCouplingSettings{};
	+ RHEJ::EventReweighterConfig conf;
	+ conf.psp_config = std::move(psp_conf);
	+ conf.ME_config = std::move(ME_conf);
	+ conf.jet_param = psp_conf.jet_param;
	conf.treat = treat;
	- conf.min_extparton_pt = extpartonptmin;
	- conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
	- conf.log_correction = log_corr;
	- conf.Higgs_coupling = RHEJ::HiggsCouplingSettings{};

	RHEJ::EventReweighter rhej(reader.heprup, conf);

	double xsec = 0.;
	while(reader.readEvent()){
	RHEJ::Event ev{
	RHEJ::UnclusteredEvent{reader.hepeup},
	Born_jet_def, Born_jetptmin
	};
	auto resummed_events = rhej.reweight(ev, 10);
	for(auto const & ev: resummed_events) xsec += ev.central().weight;
	}

	if(std::abs(xsec - xsec_ref) > tolerance){
	std::cerr << "cross section is off: "
	<< xsec << " != " << xsec_ref
	<< " +- " << tolerance << '\n';
	return EXIT_FAILURE;
	}
	}
	diff --git a/t/test_ME_h_3j.cc b/t/test_ME_h_3j.cc
	index c8fad15..ab3a28d 100644
	--- a/t/test_ME_h_3j.cc
	+++ b/t/test_ME_h_3j.cc
	@@ -1,99 +1,99 @@
	#include "LHEF/LHEF.h"
	#include "RHEJ/config.hh"
	#include "RHEJ/MatrixElement.hh"
	#include "RHEJ/Event.hh"

	struct Event{
	std::array<RHEJ::Sparticle, 2> incoming;
	std::vector<RHEJ::Sparticle> outgoing;

	Event(
	RHEJ::Sparticle in1, RHEJ::Sparticle in2,
	std::initializer_list<RHEJ::Sparticle> out
	):
	incoming{std::move(in1), std::move(in2)},
	outgoing{out}
	{}
	};

	constexpr
	RHEJ::pid::ParticleID PDG(int id){
	return static_cast<RHEJ::pid::ParticleID>(id);
	}

	constexpr double alpha_s = 0.113559;
	constexpr double mu = 125.;
	constexpr double R = 0.4;
	constexpr double min_jet_pt = 50.;
	constexpr auto jet_def = fastjet::antikt_algorithm;
	constexpr double ep = 1e-4;

	void dump(Event const & ev){
	LHEF::Writer writer{std::cerr};
	std::cerr << std::setprecision(15);
	RHEJ::UnclusteredEvent tmp;
	tmp.incoming = ev.incoming;
	tmp.outgoing = ev.outgoing;
	tmp.central = {mu, mu, 0.};
	RHEJ::Event out_ev{std::move(tmp), {jet_def, R}, min_jet_pt};
	writer.hepeup = to_HEPEUP(std::move(out_ev), nullptr);
	writer.writeEvent();
	}

	int main(){
	std::vector<Event> events = {
	/*
	* reference matrix elements obtained from traditional HEJ (svn r3355)
	* weights correspond to shat2wt/(pdfPtr->pdftapdfPtr->pdftb)
	* at the end of MEt2 in HJets.cxx
	*/
	#include "ME_test_events.dat"
	};
	const std::vector<double> ref_weights{
	1.20840366448943e-08,3.41482815561959e-18,4.09113550770465e-12,4.62704699562129e-15,4.48439126356004e-15,1.87972263060254e-12,1.44306589954128e-14,9.10660553817556e-13,6.78185181417658e-13,2.2313219279205e-14,5.04225136522914e-10,1.15445795400542e-13,4.25305992011228e-15,6.34464495221583e-11,8.1504882635896e-09,8.13050604628519e-10,1.16718768515638e-12,2.10641257506331e-08,6.62448031524987e-23,4.3319311109536e-07,1.20944494921349e-15,5.92539274973405e-05,4.4796169513938e-17,1.49273221119393e-15,2.43311569643869e-18,2.402863950606e-12,1.0620402696538e-11,5.35907540952987e-14,1.23031328706173e-12,2.79984692016006e-15,9.15681008462036e-19,6.82192590709917e-19,7.89237067022333e-13,1.86125787460036e-14,6.46979689455398e-07,3.04911830916725e-15,8.13163686285418e-15,1.53766124466997e-14,8.32276518653951e-14,8.88400099785712e-11,8.91102728936822e-06,1.80107644544759e-15,2.10582866365414e-09,3.10947796493188e-17,3.91243813845236e-15,3.26654480787734e-17,9.5447871679842e-14,2.59793669360488e-15,4.96134344054012e-13,5.81162371703576e-14,1.111167877034e-09,2.79109572058797e-16,4.46160513513727e-16,6.75218397576417e-15,2.68260561114305e-11,4.28989454505788e-16,5.8329868340234e-16,2.2024897389957e-09,2.57397955688743e-15,2.44654345522128e-14,2.30866752912176e-14,1.41827747056589e-07,1.08893147410797e-14,4.03382910318587e-11,3.1389869623674e-10,3.90505557430021e-13,1.08295237971538e-12,7.454591488958e-15,7.54483306433453e-14,1.02908259302562e-15,5.63531574394432e-14,6.30249429292698e-11,4.53762691808614e-13,1.89024041595359e-16,6.8158270448436e-15,1.14399844221906e-12,8.90256332512462e-11,6.06548092223796e-17,9.98362270779565e-12,4.59859716748231e-15,2.15099414771066e-10,2.11591109618363e-13,4.9742178019335e-09,1.33845681137089e-10,2.71186491627684e-16,8.86739836718398e-13,1.8249669781697e-14,1.92215499935624e-13,5.63538105112205e-14,4.14171568963499e-12,2.2717979957039e-15,5.55611156772101e-15,1.28364738716082e-16,5.3851134765245e-12,1.24211103260246e-08,4.83232530709888e-17,0.000521569751168655,1.6806219324588e-15,6.59977348355555e-11,9.1627970683281e-14,2.76820118228537e-12,4.0198982918336e-11,1.48426883756103e-18,1.40114232815435e-18,1.7965991343524e-11,7.6977628175923e-10,4.62196685302415e-14,2.69246069544232e-08,1.51516950568948e-13,6.99496684443045e-13,5.18573546306991e-16,0.0434630460883225,1.9837330600509e-23,2.67017601462925e-12,9.8878856049861e-16,8.72580364190897e-11,9.78094461016429e-06,8.71781003862465e-17,1.1559862791241e-08,4.84083438123941e-11,0.0139742563571754,1.99847601444777e-13,3.04915596353023e-17,5.93151845056274e-16,2.17435618150856e-14,5.33225221812805e-08,3.87773807827851e-14,4.50810202343348e-08,1.21607462589437e-08,3.1594598104258e-13,3.28656273761685e-16,2.50210924939855e-21,2.38126296013289e-13,1.25433376653286e-11,2.14322725557091e-14,1.49146107213853e-09,1.242585653879e-15,3.55954759815053e-17,3.58323190858796e-10,4.51037144634267e-14,7.55931651143187e-09,5.37698824199559e-14,7.38868781004529e-14,1.69357650230295e-11,5.80639629823124e-10,7.4837350747944e-17,6.96956083393151e-16,1.81834207579059e-10,2.82540028744915e-12,4.0395726561383e-13,1.79694539029015e-12,1.31041890677814e-16,1.30438189219159e-07,1.43221371664219e-10,9.114047729275e-16,7.34577355123687e-09,3.22593068239571e-15,4.15681922353836e-15,3.36408379142654e-16,1.44480854726148e-11,1.01363611458694e-13,1.00219358430676e-12,0.00338137229726986,4.29548017717404e-15,2.60395354258443e-14,4.1476296627427e-14,1.06855911096978e-13,8.22439039784712e-11,2.73629652512136e-10,5.15659349602634e-11,1.01584718933351e-13,6.94380398977803e-16,2.18931000347143e-18,3.07933496617785e-09,1.66390942837878e-16,3.52525907142114e-13,3.78080647432571e-14,1.17556440172277e-12,6.59751630435329e-17,1.21755902521401e-08,2.55450721249666e-11,1.53256550428156e-14,1.08380893991568e-12,4.16098986326189e-13,2.90205083341855e-14,5.511460594056e-18,1.57699398911427e-16,7.00104401760928e-13,3.58046798793545e-12,2.01229361268661e-12,3.73614518690553e-07,2.14710426845508e-14,3.70105600162546e-13,7.98071394035332e-15,7.37834963791502e-16,1.60529982053429e-09,1.08256973923745e-12,1.0334353499219e-13,9.49543099778791e-11,1.9115591364306e-12,4.0448657562838e-16,8.71422438179998e-15,1.25631778664749e-16,3.14479096095934e-19,2.64428535684163e-16,1.12107615349188e-11,1.23625298812219e-13,4.01428912903254e-13,2.70797308155832e-11,6.56279497820654e-09,7.13734059806893e-16,0.000782868291989126,3.76614026148548e-16,1.26257000484765e-13,9.58631950869197e-13,8.31708875806124e-11,2.34133301418612e-13,9.56863619261736e-14,4.64712032868663e-12,2.86185199150703e-12,1.33576465609689e-11,1.10079981208014e-15,1.0127254462106e-12,1.1462054990501e-12,6.6945137657857e-11,5.58956247731427e-17,9.55921087729154e-12,2.34217482968982e-13,9.06915368555449e-19,5.71005344841003e-15,7.30755001164554e-09,6.21879071656432e-07,3.18226338142145e-17,7.85187718829052e-15,1.12784019414768e-13,1.26620458266236e-12,3.69177699852588e-14,5.15075992359208e-08,4.92072565553936e-15,1.85781482577606e-18,4.45444190278706e-16,1.26266755594178e-13,2.4288203970545e-16,3.05176309462854e-13,1.10490541122231e-16,8.47749453725764e-16,5.15877353528789e-11,5.24254149207253e-17,1.91355494316452e-10,5.00556344079427e-17,2.17154626918891e-15,8.84358943960976e-15,3.61414423657153e-11,4.82735747825927e-09,1.29064865560832e-13,8.42127739585314e-10,3.92906216459766e-16,6.04446459041361e-17,5.89524177334148e-12,5.46194519536606e-19,1.34088171433487e-11,5.52860041827642e-15,1.75477788501097e-15,1.73667838135065e-16,1.65397026290053e-13,1.12185521855334e-09,6.08386662697057e-17,4.59141650497559e-07,1.5209325219688e-11,2.23037744763897e-13,2.01952877105987e-13,6.8597135743949e-17,1.21998302502671e-11,3.21171910404898e-17,5.26021980482437e-15,7.19705632175135e-19,1.20963419308456e-09,6.88528237531901e-13,1.98848279147693e-11,4.01262000581642e-12,5.20247739003617e-13,1.30428205622828e-14,4.11298291694491e-13,5.52043003142038e-10,1.23093415613032e-16,7.36616975809533e-17,3.45364902177034e-16,4.34985181424009e-08,5.38804276057675e-14,9.48300760869968e-13,1.22337104642822e-09,8.35787040119371e-14,1.55423397182344e-15,1.42203125672685e-13,3.18903277676011e-13,3.88652741986258e-15,7.04896809011702e-18,1.32646740227354e-21,7.99231856556843e-13,5.95802234870006e-09,1.65149306543534e-12,1.33064116256099e-11,1.56832702886549e-12,1.10299934631458e-08,3.67898989067657e-16,7.84682358054423e-17,2.57971628718788e-09,8.33182695406546e-17,1.18568992167241e-10,1.84355449893326e-12,2.55066898328171e-13,2.83673294248314e-14,1.78606190612042e-12,3.58348766614931e-12,4.59943073492537e-16,9.51522624162865e-13,8.0286195911246e-14,1.36695977212813e-16,2.66996900743536e-13,1.18874419459462e-15,4.20697505402443e-07,1.44157427471053e-12,3.88352207261609e-12,1.10352682706254e-15,9.16382540998255e-16,1.6918487768822e-14,6.27001020277801e-13,1.19982339734628e-15,2.08524585760556e-15,3.98635216491517e-11,2.29164410091261e-16,6.59993812570474e-14,9.09720299660866e-19,1.00791143935241e-13,3.48282488442434e-15,3.82462701684145e-11,1.82912920125976e-14,2.35934504263123e-12,3.07876467756239e-11,3.01958860729077e-05,1.34442938440306e-12,1.98172170171169e-11,2.54803366934403e-14,4.63525528427102e-13,1.98623090357032e-11,1.17532063246858e-12,1.62645374532443e-16,2.23093836031326e-11,5.23718977936426e-13,6.30254161980384e-15,5.64457070654369e-16,5.8598582297477e-15,1.54121478101166e-11,2.54757492800786e-11,7.24784320481052e-13,7.10590291356503e-11,1.90012029508776e-12,5.38421849395334e-05,9.2027530054952e-13,0.00802523961864568,1.69971085315301e-12,4.17128853468654e-15,2.58531014487319e-15,6.72871249369096e-12,3.68392430578061e-16,1.56192266889718e-14,8.09758960836169e-16,1.33670176619002e-07,3.40402526854793e-10,7.14029639915786e-15,1.15907463202726e-14
	};
	const std::vector<double> ref_virt{
	0.000407437434887122,0.0616303680115081,0.0619338272462471,0.00113600815259533,0.0611805205486735,0.000565698082480416,0.516976611025519,0.0361140019339901,0.00145895360905874,0.134071980986426,0.000248975936203567,0.199294625870235,0.24358305091454,4.74361172738652e-05,5.89666196640687e-05,0.00531027816302915,0.000200243383490462,0.0805261711344379,0.00120027603203917,0.00201705178179895,0.339181495419129,0.000229863490455131,0.0523929260778993,0.000263057002724373,0.00244241108460007,0.0261467924657004,5.9657372131903e-05,0.200512468910312,0.00810099141722299,0.013973033333956,0.0108154079325531,0.158574821336736,0.00070998197888614,0.000498608284914822,3.73054997210496e-05,0.0589374030166269,0.00230622721548646,0.0339477696166811,0.000912463632412847,0.000736710687806361,0.00409015383244036,0.245951477936874,0.0837392661445892,0.222078993106222,0.00532590926455614,0.625630979237125,0.0496003388518398,0.10116468200925,0.0661786669164834,0.0210544045965245,0.0244625178289293,0.0146510511459362,0.101257018174142,0.00316763508016474,0.00143243489499069,0.333592621341247,0.00370907687545401,0.00213872056006003,0.01900030889465,0.167144065343861,2.6100961817368e-05,4.91550277446228e-05,0.0166499896061438,0.00011398434602661,0.000359768751250008,0.00225147558163185,0.0317147645964709,0.0842595514695342,0.35407612826878,0.00010526694382132,0.0430608610020923,3.5715013356257e-05,0.000838215519554185,0.00192018265054966,0.0713712284182796,0.000485232329561933,0.00152471141515176,0.000311862603200811,0.166619420201183,0.0199021050517686,1.41779057397457e-05,0.0017318628596011,0.00053760812759845,0.0686125677803084,0.0532620439395034,0.00166666380089786,0.0113533421860016,0.000109349820296255,0.0115667016173654,0.00173679164546385,0.0256607191277514,0.00120155627964714,0.0111815333349668,0.00349052851306921,0.00531735228893806,0.00152975285273546,4.66037248418355e-05,0.0425850611539657,0.00117697776465851,1.23387751381228e-05,0.0038237470220905,0.0905644410585002,0.02820127090728,0.00024367136851271,0.147636957219027,0.000593136877076053,0.0618471045067722,9.07687985864349e-06,0.000158960279301719,0.0887358461047651,0.0332353914977947,7.4341716367417e-06,0.000892008177580303,0.000741022113249675,0.00237741099974918,0.00354813689304541,0.00074173961384706,0.00467141577324516,0.00842296482521057,0.00821661392007147,4.83792120557372e-05,0.00425252516494632,0.000432204419963046,0.00586736518573752,0.000956530707046858,0.000236851602228008,0.0892427215701044,1.38830962250345e-05,0.0080267826500828,0.284615916615109,0.000228095910528196,0.000135342214864163,7.55341699407512e-06,0.000571884401074523,0.0101768016413147,0.0171951107578717,0.00189669876455143,0.00394930354845142,0.00328294862282636,0.00123193340562429,0.00331783483858544,0.00319051646342221,0.00155129659005506,8.87358953340508e-05,0.27623170363077,0.000223827470760779,0.0568015505265356,0.000286151801691539,0.0193226057191209,0.00123755083758374,0.000465248541258071,0.0760142560624272,0.00112490607994365,0.000144842352216183,0.197405681090406,0.00129883979140284,0.0930965339402953,0.0010308437485118,0.670013359364299,0.00432667098624779,0.0454915871138651,0.00251703872973342,0.00423441033466031,0.000345157767746845,0.00161415356368555,0.00416069958683314,0.00164939718639514,0.0111406419016027,0.00655421252542946,0.00851532238518754,0.00348709627275328,0.496059905302338,0.0371917367183897,0.000211531938236732,0.000258942974591912,0.17759337052515,0.00418622722335427,0.166120697713642,0.00360502919950296,0.000407414278464457,0.0360738279541873,0.0288893174873495,0.00047259015965458,1.56547248986137e-05,0.834956154090444,0.00732703650642937,0.00474657800944094,0.00273096554923287,0.257207154892208,0.00469808491585188,0.000270576144489275,0.000136075839471612,0.097205919234904,0.000607126031374257,0.000136619984636834,0.000160087925098193,0.0073858227644151,0.538982397271589,0.00336284303209108,0.00152656630103797,0.367113302445716,7.48038626150403e-05,2.60037367946894e-05,4.68016275421979e-05,0.000158257654442151,0.000139074963089943,0.000665698879408135,0.000343142682314476,0.0686327476602143,0.000242720484154181,0.0343708294206813,0.00114889968298977,0.0586454267607763,0.0681515891243611,0.00825429790468652,0.023992046337532,0.0176721744671367,0.000118812692493978,0.152442537631382,0.017153140472422,0.0064172004743378,3.15311961394686e-05,0.0176719977165195,0.00523267370956132,0.00696013339412051,0.0283129688131134,0.000387153677372906,0.00691199037138631,0.000146275451696295,0.701230156661611,0.000192846626056869,0.000164329434653947,0.00161262788918379,0.138649850457724,0.0298984433917698,0.0318122273046375,3.04827880239274e-05,0.000447390422746664,0.00113026494701542,0.00113285595566244,2.89098339981399e-05,0.00113856357833898,0.0539144571449436,0.0708152617112735,0.000146823107858066,0.0226524231026193,0.00337461891944858,0.000951319767163979,0.00173528759879105,0.027553837556441,0.000644790987241654,0.00202675597720372,0.0904812026317863,0.000121859990107361,0.0778263959001868,0.0078676495423855,0.00106041104696191,0.0460181415526085,0.111418229601818,0.00485578450515022,0.000301615731939099,0.0991074788208953,0.000545043001929007,0.00699096247420379,0.203553031240112,0.00129391506623192,0.0163036481873179,9.84784360850182e-06,0.0011930992876943,0.000233138513444116,0.00327134762224621,0.00934510369593952,0.0257690022271133,0.0586620958023085,0.0044269686967688,0.00458204111385776,0.000575173870997244,0.0196091953126155,0.0172446254918132,0.0023874923763825,0.00203444177442346,0.00590189106361814,0.160798676736245,0.000218465854877618,0.0449130791723244,0.0221907455403962,0.00577930006864582,4.68129848124974e-05,4.39088532565896e-05,0.0957291068496414,0.000222196152405895,0.00244492299475677,0.0649283175892222,0.00246575641284538,0.00395895850035338,0.452609135498808,2.18247701431579e-05,0.000103691495491759,0.0300254483476907,0.000997825085990438,0.00120168803041095,0.000488310402876861,5.02375648302248e-05,0.0889218155230037,0.886196353445691,0.537472442381394,0.000989158503397272,0.00585452117968649,0.000166646041974894,0.0559913557856292,0.271624614149917,1.9604193016588e-05,0.00400455993655262,0.000719953658022318,0.00769221521719457,0.0122222590022294,0.00869539856191144,0.00665334452394116,0.00269456106141176,0.0203144580365219,1.06269256412071e-05,0.00266305904353431,0.000866750257077464,0.0118678506617259,0.00774140571086026,0.472668086097713,0.00449314655439557,0.000363054262255041,0.000123958362353204,0.0593128481039698,6.49758497533422e-05,0.0147733538718899,0.000181422709422017,0.124358363237549,0.0149418567595413,0.000496790567694211,0.00319609174497007,0.171371757877329,0.000132314177013318,0.000118421517677613,0.0384411149022143,0.0560023977555578,0.000545951922843868,0.00923849055867683,7.42750022847054e-05,0.0034657141222543,0.0814951536826408,0.000247436489650107,0.00279135741337402,0.128217103654628,0.000239656955359697,0.0286081920788208,0.0264861644616376,0.000700223039958121,0.000224261722502375,0.206999869456071,0.0245477336638465,0.00020080171595025,0.000534370708560753,0.00114853261836608,0.0140588914482069,0.0491695416950957,0.0219392429988939,0.000187573794241001,2.76073824918907e-05,0.0177905351398601,0.0227102533101537,0.0156338659659345,0.0045130866637909,0.0118789100672146,0.0684121072806861
	};
	- RHEJ::Config config;
	- config.resummation_jets.def = fastjet::JetDefinition(jet_def, R);
	- config.resummation_jets.min_pt = min_jet_pt;
	+ RHEJ::MatrixElementConfig config;
	+ config.jet_param.def = fastjet::JetDefinition(jet_def, R);
	+ config.jet_param.min_pt = min_jet_pt;
	config.log_correction = false;
	config.Higgs_coupling = RHEJ::HiggsCouplingSettings{};

	RHEJ::MatrixElement ME{
	[](double){ return alpha_s; },
	config
	};
	assert(events.size() == ref_weights.size());
	for(size_t i = 0; i < events.size(); ++i){
	std::sort(
	begin(events[i].outgoing), end(events[i].outgoing),
	RHEJ::rapidity_less{}
	);
	const double our_ME = ME.tree(
	mu, events[i].incoming, events[i].outgoing, true
	);
	const double deviation = our_ME/ref_weights[i] - 1.;
	if(std::abs(deviation) > ep){
	std::cerr
	<< "Matrix element deviates by factor of " << 1+deviation << ":\n"
	"Is " << our_ME << ", should be " << ref_weights[i] << "\n"
	"Event:\n";
	dump(events[i]);
	return EXIT_FAILURE;
	}
	const double our_virt = ME.virtual_corrections(
	mu, events[i].incoming, events[i].outgoing
	);
	const double virt_deviation = our_virt/ref_virt[i] - 1.;
	if(std::abs(virt_deviation) > ep){
	std::cerr
	<< "Virtual corrections deviate by factor of " << 1+virt_deviation << ":\n"
	"Is " << our_virt << ", should be " << ref_virt[i] << "\n"
	"Event:\n";
	dump(events[i]);
	return EXIT_FAILURE;
	}
	}
	}
	diff --git a/t/test_ME_hjets_mt174.cc b/t/test_ME_hjets_mt174.cc
	index 4c54a60..45591c1 100644
	--- a/t/test_ME_hjets_mt174.cc
	+++ b/t/test_ME_hjets_mt174.cc
	@@ -1,90 +1,90 @@
	#include "LHEF/LHEF.h"
	#include "RHEJ/MatrixElement.hh"
	#include "RHEJ/Event.hh"

	struct Event{
	double weight;
	std::array<RHEJ::Sparticle, 2> incoming;
	std::vector<RHEJ::Sparticle> outgoing;

	Event(
	double wt,
	RHEJ::Sparticle in1, RHEJ::Sparticle in2,
	std::initializer_list<RHEJ::Sparticle> out
	):
	weight{wt},
	incoming{std::move(in1), std::move(in2)},
	outgoing{out}
	{}
	};

	constexpr
	RHEJ::pid::ParticleID PDG(int id){
	return static_cast<RHEJ::pid::ParticleID>(id);
	}

	constexpr double alpha_s = 0.112654;
	constexpr double mu = 125.;
	constexpr double R = 0.4;
	constexpr double min_jet_pt = 35.;
	constexpr auto jet_def = fastjet::antikt_algorithm;
	constexpr double ep = 1e-4;
	constexpr double alpha_s_mH_RHEJ = 0.113559;

	void dump(Event const & ev){
	LHEF::Writer writer{std::cerr};
	std::cerr << std::setprecision(15);
	RHEJ::UnclusteredEvent tmp;
	tmp.incoming = ev.incoming;
	tmp.outgoing = ev.outgoing;
	tmp.central = {mu, mu, 0.};
	RHEJ::Event out_ev{std::move(tmp), {jet_def, R}, min_jet_pt};
	writer.hepeup = to_HEPEUP(std::move(out_ev), nullptr);
	writer.writeEvent();
	}

	int main(){
	std::vector<Event> events = {
	/*
	* reference matrix elements obtained from traditional HEJ (svn r3355)
	* weights correspond to shat2wt/(pdfPtr->pdftapdfPtr->pdftb)
	* at the end of MEt2 in HJets.cxx
	*/
	#include "ME_hjets_mt174.dat"
	};

	- RHEJ::Config config;
	- config.resummation_jets.def = fastjet::JetDefinition(jet_def, R);
	- config.resummation_jets.min_pt = min_jet_pt;
	+ RHEJ::MatrixElementConfig config;
	+ config.jet_param.def = fastjet::JetDefinition(jet_def, R);
	+ config.jet_param.min_pt = min_jet_pt;
	config.log_correction = false;

	RHEJ::HiggsCouplingSettings settings;
	settings.mt = 174.;
	settings.use_impact_factors = false;
	settings.mb = 5.;
	settings.include_bottom = true;
	config.Higgs_coupling = settings;

	RHEJ::MatrixElement ME{
	[](double){ return alpha_s; },
	config
	};
	for(size_t i = 0; i < events.size(); ++i){
	std::sort(
	begin(events[i].outgoing), end(events[i].outgoing),
	RHEJ::rapidity_less{}
	);
	const double our_ME = ME.tree(
	mu, events[i].incoming, events[i].outgoing, true
	);
	const double deviation = our_ME/events[i].weight - 1.;
	if(std::abs(deviation) > ep){
	std::cerr
	<< "Matrix element deviates by factor of " << 1+deviation << ":\n"
	"Is " << our_ME << ", should be " << events[i].weight << "\n"
	"Event:\n";
	dump(events[i]);
	return EXIT_FAILURE;
	}
	}
	}
	diff --git a/t/test_psp.cc b/t/test_psp.cc
	index 434d33f..fe3660e 100644
	--- a/t/test_psp.cc
	+++ b/t/test_psp.cc
	@@ -1,66 +1,66 @@
	#include "LHEF/LHEF.h"
	#include "RHEJ/stream.hh"
	#include "RHEJ/config.hh"
	#include "RHEJ/event_types.hh"
	#include "RHEJ/Event.hh"
	#include "RHEJ/PhaseSpacePoint.hh"

	namespace{
	constexpr int max_trials = 100;
	constexpr double extpartonptmin = 45.;
	constexpr double max_ext_soft_pt_fraction =
	std::numeric_limits<double>::infinity();
	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: test_psp eventfile";
	return EXIT_FAILURE;
	}

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


	- RHEJ::Config conf;
	- conf.resummation_jets = RHEJ::JetParameters{jet_def, min_jet_pt};
	+ RHEJ::PhaseSpacePointConfig conf;
	+ conf.jet_param = RHEJ::JetParameters{jet_def, min_jet_pt};
	conf.min_extparton_pt = extpartonptmin;
	conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;

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

	}

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 9:10 PM (23 h, 20 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4242384
Default Alt Text: (147 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions