No OneTemporary
Actions

Size

81 KB

Subscribers

None

View Options

	diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh
	index 2635bf8..6ffdfec 100644
	--- a/FixedOrderGen/include/PhaseSpacePoint.hh
	+++ b/FixedOrderGen/include/PhaseSpacePoint.hh
	@@ -1,187 +1,187 @@
	/** \file PhaseSpacePoint.hh
	* \brief Contains the PhaseSpacePoint Class
	*/

	#pragma once

	#include <vector>

	#include "HEJ/Event.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/RNG.hh"

	#include "JetParameters.hh"
	#include "ParticleProperties.hh"
	#include "Status.hh"

	namespace HEJFOG{
	class Process;

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

	//! PhaseSpacePoint Constructor
	/**
	* @param proc The process to generate
	* @param jet_properties Jet defintion & cuts
	* @param pdf The pdf set (used for sampling)
	* @param E_beam Energie of the beam
	* @param subl_chance Chance to turn a potentially unordered
	* emission into an actual one
	* @param subl_channels Possible subleading channels.
	* see HEJFOG::Subleading
	* @param particle_properties Properties of producted boson
	*
	* Initially, only FKL phase space points are generated. If the most
	* extremal emission in any direction is a quark or anti-quark and the
	* next emission is a gluon, subl_chance is the chance to turn this into
	* an unordered emission event by swapping the two flavours. At most one
	* unordered emission will be generated in this way.
	*/
	PhaseSpacePoint(
	Process const & proc,
	JetParameters const & jet_properties,
	HEJ::PDF & pdf, double E_beam,
	double subl_chance,
	unsigned int subl_channels,
	ParticlesPropMap const & particles_properties,
	HEJ::RNG & ran
	);

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

	Status status() const{
	return status_;
	}

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

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

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

	private:
	/**
	* \internal
	* \brief Generate LO parton momentum
	*
	* @param count Number of partons to generate
	* @param is_pure_jets If true ensures momentum conservation in x and y
	* @param jet_param Jet properties to fulfil
	* @param max_pt max allowed pt for a parton (typically E_CMS)
	* @param ran Random Number Generator
	*
	* @returns Momentum of partons
	*
	* Ensures that each parton is in its own jet.
	* Generation is independent of parton flavour. Output is sorted in rapidity.
	*/
	std::vector<fastjet::PseudoJet> gen_LO_partons(
	int count, bool is_pure_jets,
	JetParameters const & jet_param,
	double max_pt,
	HEJ::RNG & ran
	);
	Particle gen_boson(
	HEJ::ParticleID bosonid, double mass, double width,
	HEJ::RNG & ran
	);
	template<class ParticleMomenta>
	fastjet::PseudoJet gen_last_momentum(
	ParticleMomenta const & other_momenta,
	double mass_square, double y
	) const;

	bool jets_ok(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	/**
	* \internal
	* \brief Generate incoming partons according to the PDF
	*
	* @param uf Scale used in the PDF
	*/
	void reconstruct_incoming(
	std::array<HEJ::ParticleID, 2> const & ids,
	HEJ::PDF & pdf, double E_beam,
	double uf,
	HEJ::RNG & ran
	);
	HEJ::ParticleID generate_incoming_id(
	size_t beam_idx, double x, double uf, HEJ::PDF & pdf,
	HEJ::RNG & ran
	);

	bool momentum_conserved(double ep) const;

	HEJ::Particle const & most_backward_FKL(
	std::vector<HEJ::Particle> const & partons
	) const;
	HEJ::Particle const & most_forward_FKL(
	std::vector<HEJ::Particle> const & partons
	) const;
	HEJ::Particle & most_backward_FKL(std::vector<HEJ::Particle> & partons) const;
	HEJ::Particle & most_forward_FKL(std::vector<HEJ::Particle> & partons) const;
	bool extremal_FKL_ok(
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	double random_normal(double stddev, HEJ::RNG & ran);
	void maybe_turn_to_uno(double chance, HEJ::RNG & ran);
	std::vector<Particle> decay_boson(
	HEJ::Particle const & parent,
	std::vector<Decay> const & decays,
	HEJ::RNG & ran
	);
	Decay select_decay_channel(
	std::vector<Decay> const & decays,
	HEJ::RNG & ran
	);
	double gen_hard_pt(
	int np, double ptmin, double ptmax, double y,
	HEJ::RNG & ran
	);
	double gen_soft_pt(int np, double ptmax, HEJ::RNG & ran);
	double gen_parton_pt(
	int count, JetParameters const & jet_param, double ptmax, double y,
	HEJ::RNG & ran
	);


	double weight_;

	Status status_;

	std::array<Particle, 2> incoming_;
	std::vector<Particle> outgoing_;
	//! Particle decays in the format {outgoing index, decay products}
	std::unordered_map<size_t, std::vector<Particle>> decays_;
	};
	- HEJ::EventData to_EventData(PhaseSpacePoint const & psp);
	+ HEJ::Event::EventData to_EventData(PhaseSpacePoint const & psp);
	}
	diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
	index 961f6aa..fd9e165 100644
	--- a/FixedOrderGen/src/PhaseSpacePoint.cc
	+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
	@@ -1,461 +1,457 @@
	#include "PhaseSpacePoint.hh"

	#include <random>

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

	#include "HEJ/exceptions.hh"
	#include "HEJ/kinematics.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/utility.hh"

	#include "Process.hh"
	#include "Subleading.hh"


	using namespace HEJ;

	namespace HEJFOG{

	static_assert(
	std::numeric_limits<double>::has_quiet_NaN,
	"no quiet NaN for double"
	);
	constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
	- HEJ::EventData to_EventData(PhaseSpacePoint const & psp){
	- HEJ::EventData result;
	- result.incoming() = psp.incoming();
	- std::sort(
	- begin(result.incoming()), end(result.incoming()),
	- [](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
	- );
	- assert(result.incoming().size() == 2);
	- result.outgoing() = psp.outgoing();
	+ HEJ::Event::EventData to_EventData(PhaseSpacePoint const & psp){
	+ HEJ::Event::EventData result;
	+ result.set_incoming(psp.incoming());
	+ assert(result.get_incoming().size() == 2);
	+ result.set_outgoing(psp.outgoing());
	+ // technically Event::EventData doesn't have to be sorted,
	+ // but PhaseSpacePoint should be anyway
	assert(
	std::is_sorted(
	- begin(result.outgoing()), end(result.outgoing()),
	+ begin(result.get_outgoing()), end(result.get_outgoing()),
	HEJ::rapidity_less{}
	)
	);
	- assert(result.outgoing().size() >= 2);
	- result.decays() = psp.decays();
	- result.central().mur = NaN;
	- result.central().muf = NaN;
	- result.central().weight = psp.weight();
	+ assert(result.get_outgoing().size() >= 2);
	+ result.set_decays(psp.decays());
	+ result.set_central( {NaN, NaN, psp.weight() });
	return result;
	}

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

	void PhaseSpacePoint::maybe_turn_to_uno(
	double chance,
	HEJ::RNG & ran
	){
	assert(outgoing_.size() >= 2);
	const size_t nout = outgoing_.size();
	const bool can_be_uno_backward = can_swap_to_uno(
	outgoing_[0], outgoing_[1]
	);
	const bool can_be_uno_forward = can_swap_to_uno(
	outgoing_[nout-1], outgoing_[nout-2]
	);
	if(!can_be_uno_backward && !can_be_uno_forward) return;
	if(ran.flat() < chance){
	weight_ /= chance;
	if(can_be_uno_backward && can_be_uno_forward){
	if(ran.flat() < 0.5){
	std::swap(outgoing_[0].type, outgoing_[1].type);
	}
	else{
	std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
	}
	weight_ *= 2.;
	}
	else if(can_be_uno_backward){
	std::swap(outgoing_[0].type, outgoing_[1].type);
	}
	else{
	assert(can_be_uno_forward);
	std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
	}
	}
	else weight_ /= 1 - chance;
	}

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

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

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

	PhaseSpacePoint::PhaseSpacePoint(
	Process const & proc,
	JetParameters const & jet_param,
	HEJ::PDF & pdf, double E_beam,
	double subl_change,
	unsigned int subl_channels,
	ParticlesPropMap const & particles_properties,
	HEJ::RNG & ran
	)
	{
	assert(proc.njets >= 2);
	if(proc.boson
	&& particles_properties.find(*(proc.boson))
	== particles_properties.end())
	throw HEJ::missing_option("Boson "
	+std::to_string(*(proc.boson))+" can't be generated: missing properties");
	status_ = good;
	weight_ = 1;

	const int nout = proc.njets + (proc.boson?1:0);
	outgoing_.reserve(nout);

	const bool is_pure_jets = !proc.boson;
	auto partons = gen_LO_partons(
	proc.njets, is_pure_jets, jet_param, E_beam, ran
	);
	for(auto&& p_out: partons) {
	outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out)});
	}
	if(status_ != good) return;

	// create boson
	if(proc.boson){
	const auto & boson_prop = particles_properties.at(*proc.boson);
	auto boson(gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran));
	const auto pos = std::upper_bound(
	begin(outgoing_),end(outgoing_),boson,rapidity_less{}
	);
	outgoing_.insert(pos, std::move(boson));
	if(! boson_prop.decays.empty()){
	const size_t boson_idx = std::distance(begin(outgoing_), pos);
	decays_.emplace(
	boson_idx,
	decay_boson(outgoing_[boson_idx], boson_prop.decays, ran)
	);
	}
	}
	// normalisation of momentum-conserving delta function
	weight_ = pow(2M_PI, 4);

	reconstruct_incoming(proc.incoming, pdf, E_beam, jet_param.min_pt, ran);
	if(status_ != good) return;
	// set outgoing states
	most_backward_FKL(outgoing_).type = incoming_[0].type;
	most_forward_FKL(outgoing_).type = incoming_[1].type;

	if(proc.njets > 2 && subl_channels & Subleading::uno) maybe_turn_to_uno(subl_change, ran);
	}

	double PhaseSpacePoint::gen_hard_pt(
	int np , double ptmin, double ptmax, double y,
	HEJ::RNG & ran
	) {
	// heuristic parameters for pt sampling
	const double ptpar = ptmin + np/5.;
	const double arg_small_y = atan((ptmax - ptmin)/ptpar);
	const double y_cut = 3.;

	const double r1 = ran.flat();
	if(y < y_cut){
	const double pt = ptmin + ptpartan(r1arg_small_y);
	const double temp = cos(r1*arg_small_y);
	weight_ = ptptpararg_small_y/(temptemp);
	return pt;
	}

	const double ptpar2 = ptpar/(1 + 5*(y-y_cut));
	const double temp = 1. - std::exp((ptmin-ptmax)/ptpar2);
	const double pt = ptmin - ptpar2std::log(1-r1temp);
	weight_ = ptptpar2temp/(1-r1temp);
	return pt;
	}

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

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

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

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

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

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

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

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

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

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

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

	return partons;
	}

	Particle PhaseSpacePoint::gen_boson(
	HEJ::ParticleID bosonid, double mass, double width,
	HEJ::RNG & ran
	){
	if(bosonid != HEJ::pid::Higgs)
	throw HEJ::not_implemented("Production of boson "+std::to_string(bosonid)
	+" not implemented yet.");
	// Usual phase space measure
	weight_ /= 16.*pow(M_PI, 3);

	// Generate a y Gaussian distributed around 0
	// TODO: magic number only for Higgs
	const double y = random_normal(1.6, ran);
	const double r1 = ran.flat();
	const double sH = mass*(
	mass + widthtan(M_PI/2.r1 + (r1-1.)*atan(mass/width))
	);

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

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

	Particle const & PhaseSpacePoint::most_backward_FKL(
	std::vector<Particle> const & partons
	) const{
	if(!HEJ::is_parton(partons[0])) return partons[1];
	return partons[0];
	}

	Particle const & PhaseSpacePoint::most_forward_FKL(
	std::vector<Particle> const & partons
	) const{
	const size_t last_idx = partons.size() - 1;
	if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
	return partons[last_idx];
	}

	Particle & PhaseSpacePoint::most_backward_FKL(
	std::vector<Particle> & partons
	) const{
	if(!HEJ::is_parton(partons[0])) return partons[1];
	return partons[0];
	}

	Particle & PhaseSpacePoint::most_forward_FKL(
	std::vector<Particle> & partons
	) const{
	const size_t last_idx = partons.size() - 1;
	if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
	return partons[last_idx];
	}

	void PhaseSpacePoint::reconstruct_incoming(
	std::array<HEJ::ParticleID, 2> const & ids,
	HEJ::PDF & pdf, double E_beam,
	double uf,
	HEJ::RNG & ran
	){
	std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
	// calculate xa, xb
	const double sqrts=2*E_beam;
	const double xa=(incoming_[0].p.e()-incoming_[0].p.pz())/sqrts;
	const double xb=(incoming_[1].p.e()+incoming_[1].p.pz())/sqrts;
	// abort if phase space point is outside of collider energy reach
	if (xa>1. \|\| xb>1.){
	weight_=0;
	status_ = too_much_energy;
	return;
	}
	// pick pdfs
	/** @TODO:
	* ufa, ufb don't correspond to our final scale choice.
	* The HEJ scale generators currently expect a full event as input,
	* so fixing this is not completely trivial
	*/
	for(size_t i = 0; i < 2; ++i){
	if(ids[i] == pid::proton \|\| ids[i] == pid::p_bar){
	incoming_[i].type = generate_incoming_id(i, i?xb:xa, uf, pdf, ran);
	}
	else {
	incoming_[i].type = ids[i];
	}
	}
	assert(momentum_conserved(1e-7));
	}

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

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

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

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

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

	std::vector<Particle> PhaseSpacePoint::decay_boson(
	HEJ::Particle const & parent,
	std::vector<Decay> const & decays,
	HEJ::RNG & ran
	){
	const auto channel = select_decay_channel(decays, ran);
	if(channel.products.size() != 2){
	throw HEJ::not_implemented{
	"only decays into two particles are implemented"
	};
	}
	std::vector<Particle> decay_products(channel.products.size());
	for(size_t i = 0; i < channel.products.size(); ++i){
	decay_products[i].type = channel.products[i];
	}
	// choose polar and azimuth angle in parent rest frame
	const double E = parent.m()/2;
	const double theta = 2.M_PIran.flat();
	const double cos_phi = 2.*ran.flat()-1.;
	const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
	const double px = Ecos(theta)sin_phi;
	const double py = Esin(theta)sin_phi;
	const double pz = E*cos_phi;
	decay_products[0].p.reset(px, py, pz, E);
	decay_products[1].p.reset(-px, -py, -pz, E);
	for(auto & particle: decay_products) particle.p.boost(parent.p);
	return decay_products;
	}
	}
	diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
	index 28b3e7a..d5e222c 100644
	--- a/include/HEJ/Event.hh
	+++ b/include/HEJ/Event.hh
	@@ -1,307 +1,310 @@
	/** \file
	* \brief Declares the Event class and helpers
	*
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/

	#pragma once

	#include <array>
	#include <memory>
	#include <string>
	#include <unordered_map>
	#include <vector>

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

	#include "fastjet/ClusterSequence.hh"

	namespace LHEF{
	class HEPEUP;
	class HEPRUP;
	}

	namespace fastjet{
	class JetDefinition;
	}

	namespace HEJ{


	struct ParameterDescription;

	//! Event parameters
	struct EventParameters{
	double mur; /*< Value of the Renormalisation Scale /
	double muf; /*< Value of the Factorisation Scale /
	double weight; /*< Event Weight /
	//! Optional description
	std::shared_ptr<ParameterDescription> description = nullptr;
	};

	//! Description of event parameters
	struct ParameterDescription {
	//! Name of central scale choice (e.g. "H_T/2")
	std::string scale_name;
	//! Actual renormalisation scale divided by central scale
	double mur_factor;
	//! Actual factorisation scale divided by central scale
	double muf_factor;

	ParameterDescription() = default;
	ParameterDescription(
	std::string scale_name, double mur_factor, double muf_factor
	):
	scale_name{scale_name}, mur_factor{mur_factor}, muf_factor{muf_factor}
	{};
	};

	- //! Class to store general Event setup, i.e. Phase space and weights
	- class EventData{
	- public:
	- //! Default Constructor
	- EventData() = default;
	- //! Constructor from LesHouches event information
	- EventData(LHEF::HEPEUP const & hepeup);
	- //! Constructor with all values given
	- EventData(
	- std::array<Particle, 2> const & incoming,
	- std::vector<Particle> const & outgoing,
	- std::unordered_map<size_t, std::vector<Particle>> const & decays,
	- EventParameters const & central,
	- std::vector<EventParameters> const & variations
	- ):
	- incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
	- decays_(std::move(decays)),
	- central_(std::move(central)), variations_(std::move(variations))
	- {};
	- //! Move Constructor with all values given
	- EventData(
	- std::array<Particle, 2> && incoming,
	- std::vector<Particle> && outgoing,
	- std::unordered_map<size_t, std::vector<Particle>> && decays,
	- EventParameters && central,
	- std::vector<EventParameters> && variations
	- ):
	- incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
	- decays_(std::move(decays)),
	- central_(std::move(central)), variations_(std::move(variations))
	- {};
	-
	- //! Incoming Particles
	- std::array<Particle, 2> const & incoming() const{
	- return incoming_;
	- }
	- //! Incoming Particles
	- std::array<Particle, 2> & incoming(){
	- return incoming_;
	- }
	-
	- //! Outgoing Particles
	- std::vector<Particle> const & outgoing() const{
	- return outgoing_;
	- }
	- //! Outgoing Particles
	- std::vector<Particle> & outgoing(){
	- return outgoing_;
	- }
	-
	- //! Particle decays in the format {outgoing index, decay products}
	- std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
	- return decays_;
	- }
	- //! Particle decays in the format {outgoing index, decay products}
	- std::unordered_map<size_t, std::vector<Particle>> & decays(){
	- return decays_;
	- }
	-
	- //! Central parameter (e.g. scale) choice
	- EventParameters const & central() const{
	- return central_;
	- }
	- //! Central parameter (e.g. scale) choice
	- EventParameters & central(){
	- return central_;
	- }
	-
	- //! For parameter variation
	- std::vector<EventParameters> const & variations() const{
	- return variations_;
	- }
	- //! For parameter variation
	- std::vector<EventParameters> & variations(){
	- return variations_;
	- }
	-
	- private:
	- std::array<Particle, 2> incoming_;
	- std::vector<Particle> outgoing_;
	- std::unordered_map<size_t, std::vector<Particle>> decays_;
	- //! @TODO replace this by "ParameterVariations"
	- EventParameters central_;
	- //! @TODO replace this by "ParameterVariations"
	- std::vector<EventParameters> variations_;
	- };
	-
	- //! An event before jet clustering
	- //! @TODO remove in HEJ 2.3.0
	- struct [[deprecated("Use EventData instead")]] UnclusteredEvent{
	- //! Default Constructor
	- UnclusteredEvent() = default;
	- //! Constructor from LesHouches event information
	- UnclusteredEvent(LHEF::HEPEUP const & hepeup):
	- UnclusteredEvent(EventData(hepeup)){}
	- //! Constructor from EventData
	- UnclusteredEvent(EventData const & evData);
	-
	- std::array<Particle, 2> incoming; /*< Incoming Particles /
	- std::vector<Particle> outgoing; /*< Outgoing Particles /
	- //! Particle decays in the format {outgoing index, decay products}
	- std::unordered_map<size_t, std::vector<Particle>> decays;
	- //! Central parameter (e.g. scale) choice
	- EventParameters central;
	- std::vector<EventParameters> variations; /*< For parameter variation /
	- };
	+ struct UnclusteredEvent;

	/** An event with clustered jets
	*
	* This is the main HEJ 2 event class.
	* It contains kinematic information including jet clustering,
	* parameter (e.g. scale) settings and the event weight.
	*/
	class Event{
	public:
	+ class EventData;
	//! Default Event Constructor
	Event() = default;
	//! Event Constructor adding jet clustering to an bare, unclustered event
	Event(
	- EventData ev,
	+ Event::EventData const & ev,
	fastjet::JetDefinition const & jet_def, double min_jet_pt
	);
	//! Event Constructor adding jet clustering to an unclustered event
	//! @TODO remove in HEJ 2.3.0
	[[deprecated("Use constructor with EventData instead")]]
	Event(
	- UnclusteredEvent ev,
	+ UnclusteredEvent const & ev,
	fastjet::JetDefinition const & jet_def, double min_jet_pt
	- ):Event(
	- EventData{ev.incoming,ev.outgoing,ev.decays,ev.central,ev.variations},
	- jet_def, min_jet_pt
	- ){};
	+ );
	+

	//! The jets formed by the outgoing partons
	std::vector<fastjet::PseudoJet> jets() const;

	//! The corresponding event before jet clustering
	//! @TODO remove in HEJ 2.3.0
	- [[deprecated("Use data() instead")]]
	- UnclusteredEvent unclustered() const {
	- return UnclusteredEvent(ev_);
	- }
	-
	- //! The corresponding bare event before jet clustering
	- EventData const & data() const {
	- return ev_;
	- }
	-
	- //! Central parameter choice
	- EventParameters const & central() const{
	- return ev_.central();
	- }
	-
	- //! Central parameter choice
	- EventParameters & central(){
	- return ev_.central();
	- }
	+ // [[deprecated]]
	+ // UnclusteredEvent unclustered() const {
	+ // return UnclusteredEvent(ev_);
	+ // TODO what to do with this?
	+ // }

	//! Incoming particles
	std::array<Particle, 2> const & incoming() const{
	- return ev_.incoming();
	+ return incoming_;
	}
	-
	//! Outgoing particles
	std::vector<Particle> const & outgoing() const{
	- return ev_.outgoing();
	+ return outgoing_;
	}
	-
	//! Particle decays
	/**
	* The key in the returned map corresponds to the index in the
	* vector returned by outgoing()
	*/
	std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
	- return ev_.decays();
	+ return decays_;
	+ }
	+
	+ //! Central parameter choice
	+ EventParameters const & central() const{
	+ return central_;
	+ }
	+ //! Central parameter choice
	+ EventParameters & central(){
	+ return central_;
	}

	//! Parameter (scale) variations
	std::vector<EventParameters> const & variations() const{
	- return ev_.variations();
	+ return variations_;
	}
	-
	//! Parameter (scale) variations
	std::vector<EventParameters> & variations(){
	- return ev_.variations();
	+ return variations_;
	}

	//! Parameter (scale) variation
	/**
	* @param i Index of the requested variation
	*/
	EventParameters const & variations(size_t i) const{
	- return ev_.variations()[i];
	+ return variations_[i];
	}
	-
	//! Parameter (scale) variation
	/**
	* @param i Index of the requested variation
	*/
	EventParameters & variations(size_t i){
	- return ev_.variations()[i];
	+ return variations_[i];
	}

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

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

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

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

	private:
	- EventData ev_;
	+ std::array<Particle, 2> incoming_;
	+ std::vector<Particle> outgoing_;
	+ std::unordered_map<size_t, std::vector<Particle>> decays_;
	+ //! @TODO replace this by "ParameterVariations"
	+ EventParameters central_;
	+ //! @TODO replace this by "ParameterVariations"
	+ std::vector<EventParameters> variations_;
	fastjet::ClusterSequence cs_;
	double min_jet_pt_;
	event_type::EventType type_;
	- };
	+ }; // end class Event
	+
	+ //! Class to store general Event setup, i.e. Phase space and weights
	+ class Event::EventData{
	+ public:
	+ //! Default Constructor
	+ EventData() = default;
	+ //! Constructor from LesHouches event information
	+ EventData(LHEF::HEPEUP const & hepeup);
	+ //! Constructor with all values given
	+ EventData(
	+ std::array<Particle, 2> const & incoming,
	+ std::vector<Particle> const & outgoing,
	+ std::unordered_map<size_t, std::vector<Particle>> const & decays,
	+ EventParameters const & central,
	+ std::vector<EventParameters> const & variations
	+ ):
	+ incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
	+ decays_(std::move(decays)),
	+ central_(std::move(central)), variations_(std::move(variations))
	+ {};
	+ //! Move Constructor with all values given
	+ EventData(
	+ std::array<Particle, 2> && incoming,
	+ std::vector<Particle> && outgoing,
	+ std::unordered_map<size_t, std::vector<Particle>> && decays,
	+ EventParameters && central,
	+ std::vector<EventParameters> && variations
	+ ):
	+ incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
	+ decays_(std::move(decays)),
	+ central_(std::move(central)), variations_(std::move(variations))
	+ {};
	+
	+ //! Get Incoming Particles
	+ std::array<Particle, 2> const & get_incoming() const{
	+ return incoming_;
	+ }
	+ //! Set Incoming Particles
	+ void set_incoming(std::array<Particle, 2> const & in){
	+ incoming_ = in;
	+ }
	+
	+ //! Get Outgoing Particles
	+ std::vector<Particle> const & get_outgoing() const{
	+ return outgoing_;
	+ }
	+ //! Set Outgoing Particles
	+ void set_outgoing(std::vector<Particle> const & out){
	+ outgoing_ = out;
	+ }
	+
	+ //! Get Particle decays in the format {outgoing index, decay products}
	+ std::unordered_map<size_t, std::vector<Particle>> const & get_decays() const{
	+ return decays_;
	+ }
	+ //! Set Particle decays in the format {outgoing index, decay products}
	+ void set_decays(
	+ std::unordered_map<size_t, std::vector<Particle>> const & decays
	+ ){
	+ decays_ = decays;
	+ }
	+
	+ //! Get Central parameter (e.g. scale) choice
	+ EventParameters const & get_central() const{
	+ return central_;
	+ }
	+ //! Set Central parameter (e.g. scale) choice
	+ void set_central(EventParameters const & central){
	+ central_ = central;
	+ }
	+
	+ //! Get parameter variation
	+ std::vector<EventParameters> const & get_variations() const{
	+ return variations_;
	+ }
	+ //! Set parameter variation
	+ void set_variations(std::vector<EventParameters> const & variations){
	+ variations_ = variations;
	+ }
	+
	+ private:
	+ std::array<Particle, 2> incoming_;
	+ std::vector<Particle> outgoing_;
	+ std::unordered_map<size_t, std::vector<Particle>> decays_;
	+ //! @TODO replace this by "ParameterVariations"
	+ EventParameters central_;
	+ //! @TODO replace this by "ParameterVariations"
	+ std::vector<EventParameters> variations_;
	+ }; // end class EventData

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

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

	+ //! An event before jet clustering
	+ //! @TODO remove in HEJ 2.3.0
	+ struct [[deprecated("Use EventData instead")]] UnclusteredEvent{
	+ //! Default Constructor
	+ UnclusteredEvent() = default;
	+ //! Constructor from LesHouches event information
	+ UnclusteredEvent(LHEF::HEPEUP const & hepeup):
	+ UnclusteredEvent(Event::EventData(hepeup)){};
	+ //! Constructor from EventData
	+ UnclusteredEvent(Event::EventData const & evData):
	+ incoming{evData.get_incoming()}, outgoing{evData.get_outgoing()},
	+ decays{evData.get_decays()},
	+ central{evData.get_central()}, variations{evData.get_variations()} {};
	+
	+ std::array<Particle, 2> incoming; /*< Incoming Particles /
	+ std::vector<Particle> outgoing; /*< Outgoing Particles /
	+ //! Particle decays in the format {outgoing index, decay products}
	+ std::unordered_map<size_t, std::vector<Particle>> decays;
	+ //! Central parameter (e.g. scale) choice
	+ EventParameters central;
	+ std::vector<EventParameters> variations; /*< For parameter variation /
	+ };
	+
	}
	diff --git a/src/Event.cc b/src/Event.cc
	index b2cb0c5..cca568f 100644
	--- a/src/Event.cc
	+++ b/src/Event.cc
	@@ -1,392 +1,397 @@
	/**
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/Event.hh"

	#include <algorithm>
	#include <assert.h>
	#include <numeric>
	#include <utility>

	#include "LHEF/LHEF.h"

	#include "fastjet/JetDefinition.hh"

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

	namespace HEJ{

	namespace{
	constexpr int status_in = -1;
	constexpr int status_decayed = 2;
	constexpr int status_out = 1;

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

	/**
	* \brief check if final state valid for HEJ
	*
	* check if there is at most one photon, W, H, Z in the final state
	* and all the rest are quarks or gluons
	*/
	bool final_state_ok(std::vector<Particle> const & outgoing){
	bool has_AWZH_boson = false;
	for(auto const & out: outgoing){
	if(is_AWZH_boson(out.type)){
	if(has_AWZH_boson) return false;
	has_AWZH_boson = true;
	}
	else if(! is_parton(out.type)) return false;
	}
	return true;
	}

	template<class Iterator>
	Iterator remove_AWZH(Iterator begin, Iterator end){
	return std::remove_if(
	begin, end, [](Particle const & p){return is_AWZH_boson(p);}
	);
	}

	template<class Iterator>
	bool valid_outgoing(Iterator begin, Iterator end){
	return std::distance(begin, end) >= 2
	&& std::is_sorted(begin, end, rapidity_less{})
	&& std::count_if(
	begin, end, [](Particle const & s){return is_AWZH_boson(s);}
	) < 2;
	}

	/// @note that this changes the outgoing range!
	template<class ConstIterator, class Iterator>
	bool is_FKL(
	ConstIterator begin_incoming, ConstIterator end_incoming,
	Iterator begin_outgoing, Iterator end_outgoing
	){
	assert(std::distance(begin_incoming, end_incoming) == 2);
	assert(std::distance(begin_outgoing, end_outgoing) >= 2);

	// One photon, W, H, Z in the final state is allowed.
	// Remove it for remaining tests,
	end_outgoing = remove_AWZH(begin_outgoing, end_outgoing);

	// Test if this is a standard FKL configuration.
	return
	(begin_incoming->type == begin_outgoing->type)
	&& ((end_incoming-1)->type == (end_outgoing-1)->type)
	&& std::all_of(
	begin_outgoing + 1, end_outgoing - 1,
	[](Particle const & p){ return p.type == pid::gluon; }
	);
	}

	bool is_FKL(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> outgoing
	){
	assert(std::is_sorted(begin(incoming), end(incoming), pz_less{}));
	assert(valid_outgoing(begin(outgoing), end(outgoing)));

	return is_FKL(
	begin(incoming), end(incoming),
	begin(outgoing), end(outgoing)
	);
	}

	bool has_2_jets(Event const & event){
	return event.jets().size() >= 2;
	}

	/**
	* \brief Checks whether event is unordered backwards
	* @param ev Event
	* @returns Is Event Unordered Backwards
	*
	* - Checks there is more than 3 constuents in the final state
	* - Checks there is more than 3 jets
	* - Checks the most backwards parton is a gluon
	* - Checks the most forwards jet is not a gluon
	* - Checks the rest of the event is FKL
	* - Checks the second most backwards is not a different boson
	* - Checks the unordered gluon actually forms a jet
	*/
	bool is_unordered_backward(Event const & ev){
	auto const & in = ev.incoming();
	auto const & out = ev.outgoing();
	assert(std::is_sorted(begin(in), end(in), pz_less{}));
	assert(valid_outgoing(begin(out), end(out)));

	if(out.size() < 3) return false;
	if(ev.jets().size() < 3) return false;
	if(in.front().type == pid::gluon) return false;
	if(out.front().type != pid::gluon) return false;
	// When skipping the unordered emission
	// the remainder should be a regular FKL event,
	// except that the (new) first outgoing particle must not be a A,W,Z,H.
	const auto FKL_begin = next(begin(out));
	if(is_AWZH_boson(*FKL_begin)) return false;
	if(!is_FKL(in, {FKL_begin, end(out)})) return false;
	// check that the unordered gluon forms an extra jet
	const auto jets = sorted_by_rapidity(ev.jets());
	const auto indices = ev.particle_jet_indices({jets.front()});
	return indices[0] >= 0 && indices[1] == -1;
	}

	/**
	* \brief Checks for a forward unordered gluon emission
	* @param ev Event
	* @returns Is the event a forward unordered emission
	*
	* \see is_unordered_backward
	*/
	bool is_unordered_forward(Event const & ev){
	auto const & in = ev.incoming();
	auto const & out = ev.outgoing();
	assert(std::is_sorted(begin(in), end(in), pz_less{}));
	assert(valid_outgoing(begin(out), end(out)));

	if(out.size() < 3) return false;
	if(ev.jets().size() < 3) return false;
	if(in.back().type == pid::gluon) return false;
	if(out.back().type != pid::gluon) return false;
	// When skipping the unordered emission
	// the remainder should be a regular FKL event,
	// except that the (new) last outgoing particle must not be a A,W,Z,H.
	const auto FKL_end = prev(end(out));
	if(is_AWZH_boson(*prev(FKL_end))) return false;
	if(!is_FKL(in, {begin(out), FKL_end})) return false;
	// check that the unordered gluon forms an extra jet
	const auto jets = sorted_by_rapidity(ev.jets());
	const auto indices = ev.particle_jet_indices({jets.back()});
	return indices.back() >= 0 && indices[indices.size()-2] == -1;
	}

	using event_type::EventType;

	EventType classify(Event const & ev){
	if(! final_state_ok(ev.outgoing())) return EventType::bad_final_state;
	if(! has_2_jets(ev)) return EventType::no_2_jets;
	if(is_FKL(ev.incoming(), ev.outgoing())) return EventType::FKL;
	if(is_unordered_backward(ev)){
	return EventType::unordered_backward;
	}
	if(is_unordered_forward(ev)){
	return EventType::unordered_forward;
	}
	return EventType::nonHEJ;
	}
	//@}

	Particle extract_particle(LHEF::HEPEUP const & hepeup, int i){
	return Particle{
	static_cast<ParticleID>(hepeup.IDUP[i]),
	fastjet::PseudoJet{
	hepeup.PUP[i][0], hepeup.PUP[i][1],
	hepeup.PUP[i][2], hepeup.PUP[i][3]
	}
	};
	}

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

	} // namespace anonymous

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

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

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

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

	+ //! @TODO remove in HEJ 2.3.0
	+ Event::Event(
	+ UnclusteredEvent const & ev,
	+ fastjet::JetDefinition const & jet_def, double min_jet_pt
	+ ):Event(
	+ Event::EventData{
	+ ev.incoming, ev.outgoing, ev.decays, ev.central, ev.variations},
	+ jet_def, min_jet_pt
	+ ){}
	+
	Event::Event(
	- EventData ev,
	- fastjet::JetDefinition const & jet_def, double min_jet_pt
	- ):
	- ev_{std::move(ev)},
	- cs_{to_PseudoJet(filter_partons(ev_.outgoing())), jet_def},
	+ Event::EventData const & ev,
	+ fastjet::JetDefinition const & jet_def, double const min_jet_pt
	+ ): incoming_{ev.get_incoming()},
	+ outgoing_{ev.get_outgoing()},
	+ decays_{ev.get_decays()},
	+ central_{ev.get_central()},
	+ variations_{ev.get_variations()},
	+ cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
	min_jet_pt_{min_jet_pt}
	{
	// sort particles
	std::sort(
	- begin(ev_.incoming()), end(ev_.incoming()),
	+ begin(incoming_), end(incoming_),
	[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
	);

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

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

	// classify event
	type_ = classify(*this);

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

	std::vector<fastjet::PseudoJet> Event::jets() const{
	return cs_.inclusive_jets(min_jet_pt_);
	}

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

	namespace{
	// colour flow according to Les Houches standard
	// TODO: stub
	std::vector<std::pair<int, int>> colour_flow(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing
	){
	std::vector<std::pair<int, int>> result(
	incoming.size() + outgoing.size()
	);
	for(auto & col: result){
	col = std::make_pair(-1, -1);
	}
	return result;
	}
	}

	LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
	LHEF::HEPEUP result;
	result.heprup = heprup;
	result.weights = {{event.central().weight, nullptr}};
	for(auto const & var: event.variations()){
	result.weights.emplace_back(var.weight, nullptr);
	}
	size_t num_particles = event.incoming().size() + event.outgoing().size();
	for(auto const & decay: event.decays()) num_particles += decay.second.size();
	result.NUP = num_particles;
	// the following entries are pretty much meaningless
	result.IDPRUP = event.type()+1; // event ID
	result.AQEDUP = 1./128.; // alpha_EW
	//result.AQCDUP = 0.118 // alpha_QCD
	// end meaningless part
	result.XWGTUP = event.central().weight;
	result.SCALUP = event.central().muf;
	result.scales.muf = event.central().muf;
	result.scales.mur = event.central().mur;
	result.scales.SCALUP = event.central().muf;
	result.pdfinfo.p1 = event.incoming().front().type;
	result.pdfinfo.p2 = event.incoming().back().type;
	result.pdfinfo.scale = event.central().muf;
	for(Particle const & in: event.incoming()){
	result.IDUP.emplace_back(in.type);
	result.ISTUP.emplace_back(status_in);
	result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
	result.MOTHUP.emplace_back(0, 0);
	}
	for(size_t i = 0; i < event.outgoing().size(); ++i){
	Particle const & out = event.outgoing()[i];
	result.IDUP.emplace_back(out.type);
	const int status = event.decays().count(i)?status_decayed:status_out;
	result.ISTUP.emplace_back(status);
	result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
	result.MOTHUP.emplace_back(1, 2);
	}
	result.ICOLUP = colour_flow(
	event.incoming(), filter_partons(event.outgoing())
	);
	if(result.ICOLUP.size() < num_particles){
	const size_t AWZH_boson_idx = std::find_if(
	begin(event.outgoing()), end(event.outgoing()),
	[](Particle const & s){ return is_AWZH_boson(s); }
	) - begin(event.outgoing()) + event.incoming().size();
	assert(AWZH_boson_idx <= result.ICOLUP.size());
	result.ICOLUP.insert(
	begin(result.ICOLUP) + AWZH_boson_idx,
	std::make_pair(0,0)
	);
	}
	for(auto const & decay: event.decays()){
	for(auto const out: decay.second){
	result.IDUP.emplace_back(out.type);
	result.ISTUP.emplace_back(status_out);
	result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
	const size_t mother_idx = 1 + event.incoming().size() + decay.first;
	result.MOTHUP.emplace_back(mother_idx, mother_idx);
	result.ICOLUP.emplace_back(0,0);
	}
	}
	assert(result.ICOLUP.size() == num_particles);
	static constexpr double unknown_spin = 9.; //per Les Houches accord
	result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
	result.SPINUP = result.VTIMUP;
	return result;
	}

	}
	diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
	index a23923e..7c79c7e 100644
	--- a/src/EventReweighter.cc
	+++ b/src/EventReweighter.cc
	@@ -1,274 +1,270 @@
	/**
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/EventReweighter.hh"

	#include <algorithm>
	#include <assert.h>
	#include <limits>
	#include <math.h>
	#include <stddef.h>
	#include <string>
	#include <unordered_map>

	#include "fastjet/ClusterSequence.hh"

	#include "LHEF/LHEF.h"

	#include "HEJ/Event.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/PhaseSpacePoint.hh"

	namespace HEJ{
	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();

	- EventData to_EventData(PhaseSpacePoint const & psp){
	- EventData result;
	- result.incoming() = psp.incoming();
	- std::sort(
	- begin(result.incoming()), end(result.incoming()),
	- [](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
	- );
	- assert(result.incoming().size() == 2);
	- result.outgoing() = psp.outgoing();
	+ Event::EventData to_EventData(PhaseSpacePoint const & psp){
	+ Event::EventData result;
	+ result.set_incoming(psp.incoming());
	+ assert(result.get_incoming().size() == 2);
	+ result.set_outgoing(psp.outgoing());
	+ // technically Event::EventData doesn't have to be sorted,
	+ // but PhaseSpacePoint should be anyway
	assert(
	std::is_sorted(
	- begin(result.outgoing()), end(result.outgoing()),
	+ begin(result.get_outgoing()), end(result.get_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();
	+ assert(result.get_outgoing().size() >= 2);
	+ result.set_decays( psp.decays() );
	+ result.set_central({NaN, NaN, psp.weight()});
	return result;
	}

	} // namespace anonymous

	EventReweighter::EventReweighter(
	LHEF::HEPRUP const & heprup,
	ScaleGenerator scale_gen,
	EventReweighterConfig conf,
	HEJ::RNG & ran
	):
	EventReweighter{
	HEJ::Beam{
	heprup.EBMUP.first,
	{{
	static_cast<HEJ::ParticleID>(heprup.IDBMUP.first),
	static_cast<HEJ::ParticleID>(heprup.IDBMUP.second)
	}}
	},
	heprup.PDFSUP.first,
	std::move(scale_gen),
	std::move(conf),
	ran
	}
	{
	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,
	ScaleGenerator scale_gen,
	EventReweighterConfig conf,
	HEJ::RNG & ran
	):
	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_.ME_config
	},
	scale_gen_(std::move(scale_gen)),
	ran_{ran}
	{}

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

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

	std::vector<Event> EventReweighter::gen_res_events(
	Event const & ev,
	int phase_space_points
	){
	assert(ev.variations().empty());

	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_.psp_config, ran_};
	if(psp.weight() == 0.) continue;
	if(psp.incoming()[0].E() > E_beam_ \|\| psp.incoming()[1].E() > E_beam_) continue;

	resummation_events.emplace_back(
	to_EventData(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_param.def};
	return cs.inclusive_jets(param_.jet_param.min_pt).size() == ev.jets().size();
	}

	Weights
	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_;

	Weights 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;
	}

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

	return MEt2_(ev);
	}

	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){
	return fixed_order_scale_ME(ev).central;
	}
	return MEt2_.tree(ev).central;
	}

	Weights
	EventReweighter::fixed_order_scale_ME(Event const & ev) const{
	int alpha_s_power = 0;
	for(auto const & part: ev.outgoing()){
	if(is_parton(part))
	++alpha_s_power;
	else {
	switch(part.type){
	case pid::Higgs: {
	alpha_s_power += 2;
	break;
	}
	// TODO
	case pid::Wp:
	case pid::Wm:
	case pid::photon:
	case pid::Z:
	default:
	throw not_implemented("Emission of boson of unsupported type");
	}
	}
	}

	Weights 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;
	}

	} // namespace HEJ
	diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
	index dc60d30..3c91df0 100644
	--- a/src/bin/HEJ.cc
	+++ b/src/bin/HEJ.cc
	@@ -1,188 +1,188 @@
	/**
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <chrono>
	#include <iostream>
	#include <limits>
	#include <memory>
	#include <numeric>

	#include "yaml-cpp/yaml.h"

	#include "LHEF/LHEF.h"

	#include "fastjet/ClusterSequence.hh"

	#include "HEJ/CombinedEventWriter.hh"
	#include "HEJ/config.hh"
	#include "HEJ/CrossSectionAccumulator.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/get_analysis.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/ProgressBar.hh"
	#include "HEJ/stream.hh"
	#include "HEJ/Version.hh"
	#include "HEJ/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));
	}

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

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

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

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

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

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

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

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

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

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

	LHEF::Reader reader{in};
	auto heprup = reader.heprup;
	heprup.generators.emplace_back(LHEF::XMLTag{});
	heprup.generators.back().name = HEJ::Version::package_name();
	heprup.generators.back().version = HEJ::Version::String();
	HEJ::CombinedEventWriter writer{config.output, std::move(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";
	}
	HEJ::ScaleGenerator scale_gen{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	};
	auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
	assert(ran != nullptr);
	HEJ::EventReweighter hej{
	reader.heprup,
	std::move(scale_gen),
	to_EventReweighterConfig(config),
	*ran
	};

	int nevent = 0;
	std::array<int, HEJ::event_type::last_type + 1>
	nevent_type{0}, nfailed_type{0};
	auto progress = make_progress_bar(reader.heprup.XSECUP);
	HEJ::CrossSectionAccumulator xs;
	// 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;

	// calculate HEJ weight
	HEJ::Event FO_event{
	- HEJ::EventData{reader.hepeup},
	+ HEJ::Event::EventData{reader.hepeup},
	config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
	};
	auto resummed_events = hej.reweight(FO_event, config.trials);
	++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);
	xs.fill(ev);
	}
	}
	if(progress) progress->increment(FO_event.central().weight);
	} // main event loop
	std::cout << '\n';
	analysis->finalise();

	using namespace HEJ::event_type;
	std::cout<< "Events processed: " << nevent << '\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::cout << '\n' << xs << '\n';

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

	}
	diff --git a/t/check_res.cc b/t/check_res.cc
	index 108b3ee..bd40e1f 100644
	--- a/t/check_res.cc
	+++ b/t/check_res.cc
	@@ -1,98 +1,98 @@
	#include <iostream>

	#include "LHEF/LHEF.h"

	#include "HEJ/Event.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/Mixmax.hh"
	#include "HEJ/stream.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 = HEJ::EventTreatment;
	using namespace HEJ::event_type;
	HEJ::EventTreatMap treat{
	{no_2_jets, EventTreatment::discard},
	{bad_final_state, EventTreatment::discard},
	{nonHEJ, 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]);

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

	HEJ::PhaseSpacePointConfig psp_conf;
	psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin};
	psp_conf.min_extparton_pt = extpartonptmin;
	psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
	HEJ::MatrixElementConfig ME_conf;
	ME_conf.log_correction = log_corr;
	ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{};
	HEJ::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;

	reader.readEvent();
	const bool has_Higgs = std::find(
	begin(reader.hepeup.IDUP),
	end(reader.hepeup.IDUP),
	25
	) != end(reader.hepeup.IDUP);
	const double mu = has_Higgs?125.:91.188;
	HEJ::ScaleGenerator scale_gen{
	{{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1.
	};
	HEJ::Mixmax ran{};
	HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran};

	double xsec = 0.;
	double xsec_err = 0.;
	do{
	HEJ::Event ev{
	- HEJ::EventData{reader.hepeup},
	+ HEJ::Event::EventData{reader.hepeup},
	Born_jet_def, Born_jetptmin
	};
	auto resummed_events = hej.reweight(ev, 20);
	for(auto const & ev: resummed_events) {
	xsec += ev.central().weight;
	xsec_err += ev.central().weight*ev.central().weight;
	}
	} while(reader.readEvent());
	xsec_err = std::sqrt(xsec_err);
	const double significance =
	std::abs(xsec - xsec_ref) / std::sqrt( xsec_errxsec_err + tolerancetolerance );
	std::cout << xsec_ref << " +/- " << tolerance << " ~ "
	<< xsec << " +- " << xsec_err << " => " << significance << " sigma\n";

	if(significance > 3.){
	std::cerr << "Cross section is off by over 3 sigma!\n";
	return EXIT_FAILURE;
	}
	}
	diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc
	index be54a75..e41a330 100644
	--- a/t/test_ME_generic.cc
	+++ b/t/test_ME_generic.cc
	@@ -1,104 +1,104 @@

	// Generic tester for the ME for a given set of PSP
	// reference weights and PSP (as LHE file) have to be given as _individual_ files

	#include <fstream>

	#include "LHEF/LHEF.h"

	#include "HEJ/MatrixElement.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/YAMLreader.hh"
	#include "HEJ/stream.hh"

	constexpr double alpha_s = 0.118;
	constexpr double ep = 1e-6;

	void dump(HEJ::Event const & ev){
	{
	LHEF::Writer writer{std::cout};
	std::cout << std::setprecision(6);
	writer.hepeup = to_HEPEUP(std::move(ev), nullptr);
	writer.writeEvent();
	}
	std::cout << "Rapidity ordering:\n";
	for(const auto & part: ev.outgoing()){
	std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl;
	}
	}

	int main(int argn, char** argv){
	if(argn != 4 && argn != 5){
	std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml ME_weights input_file.lhe\n\n";
	return EXIT_FAILURE;
	}
	bool OUTPUT_MODE = false;
	if(argn == 5 && std::string("OUTPUT")==std::string(argv[4]))
	OUTPUT_MODE = true;
	const HEJ::Config config = HEJ::load_config(argv[1]);

	std::fstream wgt_file;
	if ( OUTPUT_MODE ) {
	std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl;
	wgt_file.open(argv[2], std::fstream::out);
	wgt_file.precision(10);
	} else {
	wgt_file.open(argv[2], std::fstream::in);
	}

	HEJ::istream in{argv[3]};
	LHEF::Reader reader{in};

	HEJ::MatrixElement ME{
	[](double){ return alpha_s; },
	HEJ::to_MatrixElementConfig(config)
	};
	double max_ratio = 0.;
	size_t idx_max_ratio = 0;
	HEJ::Event ev_max_ratio;
	double av_ratio = 0;

	size_t i = 0;
	while(reader.readEvent()){
	++i;

	HEJ::Event event{
	- HEJ::EventData{reader.hepeup},
	+ HEJ::Event::EventData{reader.hepeup},
	config.resummation_jets.def,
	config.resummation_jets.min_pt
	};
	const double our_ME = ME.tree(event).central;

	if ( OUTPUT_MODE ) {
	wgt_file << our_ME << std::endl;
	} else {
	std::string line;
	if(!std::getline(wgt_file,line)) break;
	const double ref_ME = std::stod(line);
	const double diff = std::abs(our_ME/ref_ME-1.);
	av_ratio+=diff;
	if( diff > max_ratio ) {
	max_ratio = diff;
	idx_max_ratio = i;
	ev_max_ratio = event;
	}
	if( diff > ep ){
	size_t precision(std::cout.precision());
	std::cout.precision(16);
	std::cout<< "Large difference in PSP " << i << "\nis: "<<our_ME << " should: " << ref_ME << " => difference: " << diff << std::endl;
	std::cout.precision(precision);
	dump(event);
	return EXIT_FAILURE;
	}
	}
	}
	wgt_file.close();
	if ( !OUTPUT_MODE ) {
	size_t precision(std::cout.precision());
	std::cout.precision(16);
	std::cout << "Avg ratio after " << i << " PSP: " << av_ratio/i << std::endl;
	std::cout << "maximal ratio at " << idx_max_ratio << ": " << max_ratio << std::endl;
	std::cout.precision(precision);
	}
	return EXIT_SUCCESS;
	}
	diff --git a/t/test_classify.cc b/t/test_classify.cc
	index e0b58b4..e10aa42 100644
	--- a/t/test_classify.cc
	+++ b/t/test_classify.cc
	@@ -1,50 +1,50 @@
	#include "LHEF/LHEF.h"
	#include "HEJ/stream.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/Event.hh"

	namespace{
	constexpr double min_jet_pt = 30.;
	const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
	using namespace HEJ::event_type;
	static const std::vector<EventType> results{
	unob,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,unob,FKL,FKL,FKL,unof,FKL,unob,FKL,
	FKL,unob,unob,FKL,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unof,
	FKL,FKL,unof,FKL,FKL,FKL,FKL,FKL,unof,FKL,FKL,FKL,unof,FKL,FKL,unob,unof,
	FKL,unof,FKL,unob,FKL,FKL,unob,FKL,unob,unof,unob,unof,FKL,FKL,FKL,FKL,FKL,
	FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,
	FKL,FKL,FKL,unof,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,
	FKL,FKL,FKL,FKL,unob,FKL,unob,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unof,unob,FKL
	};
	}

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

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

	for(auto const & expected: results){
	reader.readEvent();
	const HEJ::Event ev{
	- HEJ::EventData{reader.hepeup},
	+ HEJ::Event::EventData{reader.hepeup},
	jet_def, min_jet_pt
	};

	if(ev.type() != expected){
	using HEJ::event_type::names;
	writer.hepeup = reader.hepeup;
	std::cerr << "wrong classification of event:\n";
	writer.writeEvent();
	std::cerr << "classified as " << names[ev.type()]
	<< ", is " << names[expected] << '\n';
	return EXIT_FAILURE;
	}
	}

	}
	diff --git a/t/test_descriptions.cc b/t/test_descriptions.cc
	index f796b02..ed61be2 100644
	--- a/t/test_descriptions.cc
	+++ b/t/test_descriptions.cc
	@@ -1,62 +1,61 @@
	#include <iostream>
	#include <cstddef>

	#include "HEJ/Event.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/ScaleFunction.hh"

	#define ASSERT(x) if(!(x)) { \
	std::cerr << "Assertion '" #x "' failed.\n"; \
	return EXIT_FAILURE; \
	}

	int main() {
	constexpr double mu = 125.;
	HEJ::ScaleFunction fun{"125", HEJ::FixedScale{mu}};
	ASSERT(fun.name() == "125");

	HEJ::ScaleGenerator scale_gen{
	{std::move(fun)}, {0.5, 1, 2.}, 2.1
	};
	- HEJ::EventData tmp;
	- tmp.outgoing().push_back(
	- {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)}
	- );
	- tmp.outgoing().push_back(
	+ HEJ::Event::EventData tmp;
	+ tmp.set_outgoing({
	+ {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)},
	{HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.)}
	+ }
	);
	HEJ::Event ev{
	std::move(tmp),
	fastjet::JetDefinition{fastjet::kt_algorithm, 0.4},
	20.
	};

	auto rescaled = scale_gen(std::move(ev));
	ASSERT(rescaled.central().description->scale_name == "125");
	for(auto const & var: rescaled.variations()) {
	ASSERT(var.description->scale_name == "125");
	}
	ASSERT(rescaled.central().description->mur_factor == 1.);
	ASSERT(rescaled.central().description->muf_factor == 1.);

	ASSERT(rescaled.variations(0).description->mur_factor == 1.);
	ASSERT(rescaled.variations(0).description->muf_factor == 1.);

	ASSERT(rescaled.variations(1).description->mur_factor == 0.5);
	ASSERT(rescaled.variations(1).description->muf_factor == 0.5);

	ASSERT(rescaled.variations(2).description->mur_factor == 0.5);
	ASSERT(rescaled.variations(2).description->muf_factor == 1.);

	ASSERT(rescaled.variations(3).description->mur_factor == 1.);
	ASSERT(rescaled.variations(3).description->muf_factor == 0.5);

	ASSERT(rescaled.variations(4).description->mur_factor == 1.);
	ASSERT(rescaled.variations(4).description->muf_factor == 2.);

	ASSERT(rescaled.variations(5).description->mur_factor == 2.);
	ASSERT(rescaled.variations(5).description->muf_factor == 1.);

	ASSERT(rescaled.variations(6).description->mur_factor == 2.);
	ASSERT(rescaled.variations(6).description->muf_factor == 2.);
	}
	diff --git a/t/test_psp.cc b/t/test_psp.cc
	index 35ce3e7..507319b 100644
	--- a/t/test_psp.cc
	+++ b/t/test_psp.cc
	@@ -1,69 +1,69 @@
	#include "LHEF/LHEF.h"
	#include "HEJ/stream.hh"
	#include "HEJ/config.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/PhaseSpacePoint.hh"
	#include "HEJ/Ranlux64.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";
	+ std::cerr << "Usage: " << argv[0] << " eventfile";
	return EXIT_FAILURE;
	}

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


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

	HEJ::Ranlux64 ran{};

	while(reader.readEvent()){
	const HEJ::Event ev{
	- HEJ::EventData{reader.hepeup},
	+ HEJ::Event::EventData{reader.hepeup},
	jet_def, min_jet_pt
	};
	for(int trial = 0; trial < max_trials; ++trial){
	HEJ::PhaseSpacePoint psp{ev, conf, ran};
	if(psp.weight() != 0){
	- HEJ::EventData tmp_ev;
	- tmp_ev.incoming() = psp.incoming();
	- tmp_ev.outgoing() = psp.outgoing();
	- tmp_ev.central() = {0,0,0};
	+ HEJ::Event::EventData tmp_ev;
	+ tmp_ev.set_incoming( psp.incoming() );
	+ tmp_ev.set_outgoing( psp.outgoing() );
	+ tmp_ev.set_central( {0,0,0} );
	HEJ::Event out_ev{
	std::move(tmp_ev),
	jet_def, min_jet_pt
	};
	if(out_ev.type() != ev.type()){
	using HEJ::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;
	}
	}
	}
	}

	}
	diff --git a/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc
	index 484357a..079a398 100644
	--- a/t/test_scale_arithmetics.cc
	+++ b/t/test_scale_arithmetics.cc
	@@ -1,87 +1,87 @@

	// Generic tester for the ME for a given set of PSP
	// reference weights and PSP (as LHE file) have to be given as _individual_ files

	#include <fstream>

	#include "LHEF/LHEF.h"

	#include "HEJ/EventReweighter.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/YAMLreader.hh"
	#include "HEJ/stream.hh"

	constexpr double alpha_s = 0.118;
	constexpr double ep = 1e-13;

	void dump(HEJ::Event const & ev){
	{
	LHEF::Writer writer{std::cout};
	std::cout << std::setprecision(6);
	writer.hepeup = to_HEPEUP(std::move(ev), nullptr);
	writer.writeEvent();
	}
	std::cout << "Rapidity ordering:\n";
	for(const auto & part: ev.outgoing()){
	std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl;
	}
	}

	int main(int argn, char** argv){
	if(argn != 3){
	std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n";
	return EXIT_FAILURE;
	}
	HEJ::Config config = HEJ::load_config(argv[1]);
	config.scales = HEJ::to_ScaleConfig(
	YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2H_T/2, H_T/22, H_T/2/24, H_TH_T/H_T]")
	);

	HEJ::istream in{argv[2]};
	LHEF::Reader reader{in};

	auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);

	HEJ::ScaleGenerator scale_gen{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	};

	HEJ::EventReweighter resum{
	reader.heprup,
	std::move(scale_gen),
	to_EventReweighterConfig(config),
	*ran
	};

	size_t i = 0;
	while(reader.readEvent()){
	++i;

	HEJ::Event event{
	- HEJ::EventData{reader.hepeup},
	+ HEJ::Event::EventData{reader.hepeup},
	config.resummation_jets.def,
	config.resummation_jets.min_pt
	};

	auto resummed = resum.reweight(event, config.trials);
	for(auto && ev: resummed) {
	for(auto &&var: ev.variations()) {
	if(std::abs(var.muf - ev.central().muf) > ep) {
	std::cerr
	<< std::setprecision(15)
	<< "unequal scales: " << var.muf
	<< " != " << ev.central().muf << '\n'
	<< "in resummed event:\n";
	dump(ev);
	std::cerr << "\noriginal event:\n";
	dump(event);
	return EXIT_FAILURE;
	}
	}
	}
	}

	}
	diff --git a/t/test_scale_import.cc b/t/test_scale_import.cc
	index 24ac7f3..e26ac1f 100644
	--- a/t/test_scale_import.cc
	+++ b/t/test_scale_import.cc
	@@ -1,31 +1,31 @@
	#include <stdexcept>
	#include <iostream>

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

	int main(int argc, char** argv) {
	constexpr double ep = 1e-7;
	if (argc != 2) {
	throw std::logic_error{"wrong number of args"};
	}
	const HEJ::Config config = HEJ::load_config(argv[1]);

	- HEJ::EventData tmp;
	- tmp.outgoing().push_back(
	- {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)}
	- );
	- tmp.outgoing().push_back(
	+ HEJ::Event::EventData tmp;
	+ tmp.set_outgoing(
	+ {
	+ {HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.)},
	{HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.)}
	+ }
	);
	HEJ::Event ev{
	std::move(tmp),
	fastjet::JetDefinition{fastjet::kt_algorithm, 0.4},
	20.
	};

	const double softest_pt = config.scales.base[0](ev);
	if(std::abs(softest_pt-30.) > ep){
	throw std::logic_error{"wrong softest pt"};
	}
	}

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 11:33 PM (1 d, 6 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4229940
Default Alt Text: (81 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions