No OneTemporary
Actions

Size

49 KB

Subscribers

None

View Options

	diff --git a/src/Event.cc b/src/Event.cc
	index f690529..bdce1d9 100644
	--- a/src/Event.cc
	+++ b/src/Event.cc
	@@ -1,1242 +1,1242 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \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/Constants.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;
	}

	/**
	* \brief function which determines if type change is consistent with W emission.
	* @param in incoming Particle
	* @param out outgoing Particle
	*
	* Ensures that change type of quark line is possible by a flavour changing
	* W emission.
	*/
	bool is_W_Current(ParticleID in, ParticleID out){
	if((in==1 && out==2)\|\|(in==2 && out==1)){
	return true;
	}
	else if((in==-1 && out==-2)\|\|(in==-2 && out==-1)){
	return true;
	}
	else if((in==3 && out==4)\|\|(in==4 && out==3)){
	return true;
	}
	else if((in==-3 && out==-4)\|\|(in==-4 && out==-3)){
	return true;
	}
	else{
	return false;
	}
	}

	/**
	* \brief checks if particle type remains same from incoming to outgoing
	* @param in incoming Particle
	* @param out outgoing Particle
	*/
	bool is_Pure_Current(ParticleID in, ParticleID out){
	if(abs(in)<=6 \|\| in==21) return (in==out);
	else return false;
	}


	// @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);

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

	template<class ConstIterator, class Iterator>
	bool is_W_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);

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

	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)));

	const auto WEmit = std::find_if(
	begin(outgoing), end(outgoing),
	[](Particle const & s){ return abs(s.type) == pid::Wp; }
	);

	if (WEmit != end(outgoing) && abs(WEmit->type) == pid::Wp){
	return is_W_FKL(
	begin(incoming), end(incoming),
	begin(outgoing), end(outgoing)
	);
	}
	else{
	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 = 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 = ev.jets();
	const auto indices = ev.particle_jet_indices({jets.back()});
	return indices.back() >= 0 && indices[indices.size()-2] == -1;
	}


	/**
	* \brief Checks for a forward extremal qqx
	* @param ev Event
	* @returns Is the event a forward extremal qqx event
	*
	* Checks there is 3 or more than 3 constituents in the final state
	* Checks there is 3 or more than 3 jets
	* Checks most forwards incoming is gluon
	* Checks most extremal particle is not a Higgs (either direction)
	* Checks the second most forwards particle is not Higgs boson
	* Checks remaining chain is pure gluon
	* Checks the most forwards parton is a either quark or anti-quark.
	* Checks the second most forwards parton is anti-quark or quark.
	* Checks the qqbar pair form 2 separate jets.
	*/
	bool is_Ex_qqxf(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)));

	int fkl_start=is_AWZ_boson(out.front());
	int fkl_end=2;

	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::Higgs \|\| out.front().type == pid::Higgs
	\|\| out.rbegin()[1].type == pid::Higgs) return false;

	// if extremal AWZ
	if(is_AWZ_boson(out.back())){ // if extremal AWZ
	++fkl_end;
	if (is_quark(out.rbegin()[1])){ //if second quark
	if (!(is_antiquark(out.rbegin()[2]))) return false;// third must be anti-quark
	}
	else if (is_antiquark(out.rbegin()[1])){ //if second anti-quark
	if (!(is_quark(out.rbegin()[2]))) return false;// third must be quark
	}
	else return false;
	}
	else if (is_quark(out.rbegin()[0])){ //if extremal quark
	if(is_AWZ_boson(out.rbegin()[1])){ // if second AWZ
	++fkl_end;
	if (!(is_antiquark(out.rbegin()[2]))) return false;// third must be anti-quark
	}
	else if (!(is_antiquark(out.rbegin()[1]))) return false;// second must be anti-quark
	}
	else if (is_antiquark(out.rbegin()[0])){ //if extremal anti-quark
	if(is_AWZ_boson(out.rbegin()[1])){ // if second AWZ
	++fkl_end;
	if (!(is_quark(out.rbegin()[2]))) return false;// third must be quark
	}
	else if (!(is_quark(out.rbegin()[1]))) return false;// second must be quark
	}
	else return false;

	// When skipping the qqbar
	// New last outgoing particle must not be a Higgs
	if (out.rbegin()[fkl_end].type == pid::Higgs) return false;

	// no other quark emission (first doesn't matter)
	if(std::any_of(
	out.begin()+1+fkl_start, out.end()-fkl_end,
	[](Particle const & p){ return is_anyquark(p); })
	) return false;

	const auto & jets = ev.jets();
	const auto indices = ev.particle_jet_indices({jets});

	// Ensure qqbar pair are in separate jets
	if(indices[indices.size()-2] != indices[indices.size()-1]-1) return false;

	// Opposite current should be logical to process
	if (is_AWZ_boson(out.front().type)){
	return (is_Pure_Current(in.front().type, out[1].type)
	\|\| is_W_Current(in.front().type,out[1].type));
	}
	return (is_Pure_Current(in.front().type, out[0].type)
	\|\| is_W_Current(in.front().type,out[0].type));
	}

	/**
	* \brief Checks for a backward extremal qqx
	* @param ev Event
	* @returns Is the event a backward extremal qqx event
	*
	* Checks there is 3 or more than 3 constituents in the final state
	* Checks there is 3 or more than 3 jets
	* Checks most backwards incoming is gluon
	* Checks most extremal particle is not a Higgs (either direction) y
	* Checks the second most backwards particle is not Higgs boson y
	* Checks remaining chain is pure gluon
	* Checks the most backwards parton is a either quark or anti-quark. y
	* Checks the second most backwards parton is anti-quark or quark. y
	* Checks the qqbar pair form 2 separate jets.
	*/
	bool is_Ex_qqxb(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)));

	int fkl_start=2;
	int fkl_end=is_AWZ_boson(out.back());

	if(out.size() < 3) return false;
	if(ev.jets().size() < 3) return false;
	if(in.front().type != pid::gluon) return false;
	if(out.back().type == pid::Higgs \|\| out.front().type == pid::Higgs
	\|\| out[1].type == pid::Higgs) return false;

	// @TODO why don't we test to the Higgs as well? i.e. is_AWZH_boson
	if(is_AWZ_boson(out.front())){ // if extremal AWZ
	++fkl_start;
	if (is_quark(out[1])){ //if second quark
	if (!(is_antiquark(out[2]))) return false;// third must be anti-quark
	}
	else if (is_antiquark(out[1])){ //if second anti-quark
	if (!(is_quark(out[2]))) return false;// third must be quark
	}
	else return false;
	}
	else if (is_quark(out[0])){ // if extremal quark
	if(is_AWZ_boson(out[1])){ // if second AWZ
	++fkl_start;
	if (!(is_antiquark(out[2]))) return false;// third must be anti-quark
	}
	else if (!(is_antiquark(out[1]))) return false;// second must be anti-quark
	}
	else if (is_antiquark(out[0])){ //if extremal anti-quark
	if(is_AWZ_boson(out[1])){ // if second AWZ
	++fkl_start;
	if (!(is_quark(out[2]))) return false;// third must be quark
	}
	else if (!(is_quark(out[1]))) return false;// second must be quark
	}
	else return false;

	// When skipping the qqbar
	// New last outgoing particle must not be a Higgs.
	if (out[fkl_start].type == pid::Higgs) return false;

	// no other quark emission (last doesn't matter)
	if(std::any_of(
	out.cbegin()+fkl_start, out.cend()-1-fkl_end,
	[](Particle const & p){ return is_anyquark(p); })
	) return false;

	const auto & jets = ev.jets();
	const auto indices = ev.particle_jet_indices({jets});

	// Ensure qqbar pair form separate jets.
	if(indices[0] != indices[1]-1) return false;

	// Other current should be logical to process
	if (is_AWZ_boson(out.back())){
	return (is_Pure_Current(in.back().type, out.rbegin()[1].type)
	\|\| is_W_Current(in.back().type,out.rbegin()[1].type));
	}
	else
	return (is_Pure_Current(in.back().type, out.rbegin()[0].type)
	\|\| is_W_Current(in.back().type, out.rbegin()[0].type));
	}


	/**
	* \brief Checks for a central qqx
	* @param ev Event
	* @returns Is the event a central extremal qqx event
	*
	* Checks there is 4 or more than 4 constuents in the final state
	* Checks there is 4 or more than 4 jets
	* Checks most extremal particle is not a Higgs (either direction) y
	* Checks for a central quark in the outgoing states
	* Checks for adjacent anti-quark parton. (allowing for AWZ boson emission between)
	* Check that other partons are only gluons in chain
	* Checks external currents are logically sound.
	*/
	bool is_Mid_qqx(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() < 4) return false;
	if(ev.jets().size() < 4) return false;

	if(out.back().type == pid::Higgs \|\| out.front().type == pid::Higgs)
	return false;

	size_t start_FKL=0;
	size_t end_FKL=0;

	if (is_AWZ_boson(out.back())){
	++end_FKL;
	}
	if (is_AWZ_boson(out.front())){
	++start_FKL;
	}

	if ((is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type)
	&& is_Pure_Current(in.front().type,out[start_FKL].type))){
	//nothing to do
	}
	else if (is_W_Current(in.back().type,out.rbegin()[end_FKL].type)
	&& is_Pure_Current(in.front().type,out[start_FKL].type)){
	//nothing to do
	}
	else if (!(is_Pure_Current(in.back().type,out.rbegin()[end_FKL].type)
	&& is_W_Current(in.front().type,out[start_FKL].type))){
	return false;
	}

	const auto & jets = ev.jets();
	const auto indices = ev.particle_jet_indices({jets});
	auto const out_partons = filter_partons(out);

	// search for qqx pair
	for (size_t i = 1; i<out_partons.size()-2; ++i){
	if ((is_quark(out_partons[i]) && (is_antiquark(out_partons[i+1])))
	\|\| (is_antiquark(out_partons[i]) && (is_quark(out_partons[i+1])))
	){
	// no quarks before qqx (beside first)
	if(std::any_of(
	out_partons.cbegin()+1, out_partons.cbegin()+i,
	[](Particle const & p){ return is_anyquark(p); })
	) return false;
	// no quarks after qqx (beside last)
	if(std::any_of(
	out_partons.cbegin()+i+2, out_partons.cend()-1,
	[](Particle const & p){ return is_anyquark(p); })
	) return false;
	// should be in seperate jets
	return (indices[i+1] == indices[i]+1 && indices[i] != -1);
	}
	}
	return false;
	}

	using event_type::EventType;

	/**
	* \brief Checks for all event types
	* @param ev Event
	* @returns Event Type
	*
	*/
	EventType classify_new(Event const & ev){
	bool IFfuno(false),IFfgqq(false), IFfLL(false),
	IFbuno(false), IFbgqq(false), IFbLL(false);
	bool IFfunoWp(false),IFfgqqWp(false), IFfLLWp(false),
	IFbunoWp(false), IFbgqqWp(false), IFbLLWp(false);
	bool IFfunoWm(false),IFfgqqWm(false), IFfLLWm(false),
	IFbunoWm(false), IFbgqqWm(false), IFbLLWm(false);
	bool Midqqbar(false), MidqqbarWp(false), MidqqbarWm(false);

	#ifndef NDEBUG
	std::cerr <<"got started\n"<<std::endl;
	#endif
	if(! final_state_ok(ev.outgoing()))
	return EventType::bad_final_state;
	if(! has_2_jets(ev))
	return EventType::no_2_jets;
	// now check event parton by parton
	// check first the impact factors at either end, and loop over the middle partons
	auto const & in = ev.incoming();
	auto const & out = filter_partons(ev.outgoing());
	auto outN=out.size();
	assert(std::distance(begin(in), end(in)) == 2);
	assert(outN >= 2);
	assert(std::distance(begin(out), end(out)) >= 2);
	assert(std::is_sorted(begin(out), end(out), rapidity_less{}));

	long unsigned int thispartonN=0;
	long unsigned int nextpartonN=1;
	if (in.front().type == out.at(thispartonN).type) { // standard jet LL vertex
	IFbLL=true;
	++thispartonN;
	++nextpartonN;
	} else if ((in.front().type==pid::u&&out.at(thispartonN).type==pid::d)
	\|\|(in.front().type==pid::c&&out.at(thispartonN).type==pid::s)
	\|\|(in.front().type==pid::d_bar&&out.at(thispartonN).type==pid::u_bar)
	\|\|(in.front().type==pid::s_bar&&out.at(thispartonN).type==pid::c_bar)){
	// if in SU(2)-up and out SU(2)-down then LLWp-type
	IFbLLWp=true;
	++thispartonN;
	++nextpartonN;
	} else if ((in.front().type==pid::u_bar&&out.at(thispartonN).type==pid::d_bar)
	\|\|(in.front().type==pid::c_bar&&out.at(thispartonN).type==pid::s_bar)
	\|\|(in.front().type==pid::d&&out.at(thispartonN).type==pid::u)
	\|\|(in.front().type==pid::s&&out.at(thispartonN).type==pid::c)){
	// if in SU(2)-down and out SU(2)-up then LLWm-type
	IFbLLWm=true;
	std::cerr <<"here\n";
	++thispartonN;
	++nextpartonN;
	} else if (outN>=3) {
	// Possibility of unordered gluon emission or gluon splitting
	if ((in.front().type!=pid::gluon)&&out.at(thispartonN).type==pid::gluon) {
	// Possibility of unordered gluon emission
	if (in.front().type == out.at(nextpartonN).type) { // unordered jet vertex
	IFbuno=true;
	thispartonN+=2;
	nextpartonN+=2;
	}
	else if ((in.front().type==pid::u&&out.at(nextpartonN).type==pid::d)
	\|\|(in.front().type==pid::c&&out.at(nextpartonN).type==pid::s)
	\|\|(in.front().type==pid::d_bar&&out.at(nextpartonN).type==pid::u_bar)
	\|\|(in.front().type==pid::s_bar&&out.at(nextpartonN).type==pid::c_bar)){
	// if in SU(2)-up and out SU(2)-down then unoWp-type
	IFbunoWp=true;
	thispartonN+=2;
	nextpartonN+=2;
	}
	else if ((in.front().type==pid::u_bar&&out.at(nextpartonN).type==pid::d_bar)
	\|\|(in.front().type==pid::c_bar&&out.at(nextpartonN).type==pid::s_bar)
	\|\|(in.front().type==pid::d&&out.at(nextpartonN).type==pid::u)
	\|\|(in.front().type==pid::s&&out.at(nextpartonN).type==pid::c)){
	// if in SU(2)-down and out SU(2)-up then unoWm-type
	IFbunoWm=true;
	thispartonN+=2;
	nextpartonN+=2;
	}
	} else if (in.front().type==pid::gluon) {
	// possibility of backward gluon splitting
	if (out.at(thispartonN).type==-out.at(nextpartonN).type) {
	IFbgqq=true;
	thispartonN+=2;
	nextpartonN+=2;
	} else if ((out.at(thispartonN).type==pid::u_bar&&out.at(nextpartonN).type==pid::d)
	\|\|(out.at(thispartonN).type==pid::c_bar&&out.at(nextpartonN).type==pid::s)
	\|\|(out.at(thispartonN).type==pid::d&&out.at(nextpartonN).type==pid::u_bar)
	\|\|(out.at(thispartonN).type==pid::s&&out.at(nextpartonN).type==pid::c_bar)){
	// if in SU(2)-up and out SU(2)-down then gqqWp-type
	IFbgqqWp=true;
	thispartonN+=2;
	nextpartonN+=2;
	} else if ((out.at(thispartonN).type==pid::u&&out.at(nextpartonN).type==pid::d_bar)
	\|\|(out.at(thispartonN).type==pid::c&&out.at(nextpartonN).type==pid::s_bar)
	\|\|(out.at(thispartonN).type==pid::d_bar&&out.at(nextpartonN).type==pid::u)
	\|\|(out.at(thispartonN).type==pid::s_bar&&out.at(nextpartonN).type==pid::c)){
	// if in SU(2)-down and out SU(2)-up then gqqWm-type
	IFbgqqWm=true;
	thispartonN+=2;
	nextpartonN+=2;
	}
	}
	}
	// if we reached this far and thisparton==begin(out)
	// then it is not a valid all-order state
	if (thispartonN==0) return EventType::FixedOrder;
	#ifndef NDEBUG
	std::cerr <<"got this far\n"<<std::endl;
	std::cerr << "IFbLLWp : "<<IFbLLWp<<std::endl;
	std::cerr << "IFbunoWp : "<<IFbunoWp<<std::endl;
	std::cerr << "IFbLL : "<<IFbLL<<std::endl;
	#endif
	// else check forward particles.
	// Need to ensure the same particles don't count for both sub-leading
	// impact factors
	long unsigned int thispartonfN=outN-1;
	long unsigned int prevpartonfN=outN-2;

	if (in.back().type == out.at(thispartonfN).type) {
	// standard forward jet LL vertex
	IFfLL=true;
	--thispartonfN;
	--prevpartonfN;
	} else if ((in.back().type==pid::u&&out.at(thispartonfN).type==pid::d)
	\|\|(in.back().type==pid::c&&out.at(thispartonfN).type==pid::s)
	\|\|(in.back().type==pid::d_bar&&out.at(thispartonfN).type==pid::u_bar)
	\|\|(in.back().type==pid::s_bar&&out.at(thispartonfN).type==pid::c_bar)){
	// if in SU(2)-up and out SU(2)-down then LLWp-type
	IFfLLWp=true;
	--thispartonfN;
	--prevpartonfN;
	} else if ((in.back().type==pid::u_bar&&out.at(thispartonfN).type==pid::d_bar)
	\|\|(in.back().type==pid::c_bar&&out.at(thispartonfN).type==pid::s_bar)
	\|\|(in.back().type==pid::d&&out.at(thispartonfN).type==pid::u)
	\|\|(in.back().type==pid::s&&out.at(thispartonfN).type==pid::c)){
	// if in SU(2)-down and out SU(2)-up then LLWm-type
	IFfLLWm=true;
	--thispartonfN;
	--prevpartonfN;
	} else if (outN>=3&&(prevpartonfN!=(thispartonN-1))) { // if the prevparton is not
	// the last parton included in backward IF
	// Possibility of unordered gluon emission or gluon splitting
	// unless we already used up all the partons
	if ((in.back().type!=pid::gluon)&&out.at(thispartonfN).type==pid::gluon) {
	#ifndef NDEBUG
	std::cerr <<"got this far2\n"<<std::endl;
	std::cerr << "thispartonfN, prevpartonfN : "<<thispartonfN <<", "<<prevpartonfN<<std::endl;
	std::cerr <<"out.at(prevpartonfN).type : "<<out.at(prevpartonfN).type<<std::endl;
	std::cerr <<"out.at(thispartonfN).type : "<<out.at(thispartonfN).type<<std::endl;
	#endif
	// possibility of unordered emissions
	if (in.back().type == out.at(prevpartonfN).type) { // unordered jet vertex
	IFfuno=true;
	thispartonfN-=2;
	prevpartonfN-=2;
	}
	else if ((in.back().type==pid::u&&out.at(prevpartonfN).type==pid::d)
	\|\|(in.back().type==pid::c&&out.at(prevpartonfN).type==pid::s)
	\|\|(in.back().type==pid::d_bar&&out.at(prevpartonfN).type==pid::u_bar)
	\|\|(in.back().type==pid::s_bar&&out.at(prevpartonfN).type==pid::c_bar)){
	// if in SU(2)-up and out SU(2)-down then unoWp-type
	IFfunoWp=true;
	thispartonfN-=2;
	prevpartonfN-=2;
	}
	else if ((in.back().type==pid::u_bar&&out.at(prevpartonfN).type==pid::d_bar)
	\|\|(in.back().type==pid::c_bar&&out.at(prevpartonfN).type==pid::s_bar)
	\|\|(in.back().type==pid::d&&out.at(prevpartonfN).type==pid::u)
	\|\|(in.back().type==pid::s&&out.at(prevpartonfN).type==pid::c)){
	// if in SU(2)-down and out SU(2)-up then unoWm-type
	IFfunoWm=true;
	thispartonfN-=2;
	prevpartonfN-=2;
	}
	} else if (in.back().type==pid::gluon) {
	// possibility of forward gluon splitting
	if (out.at(thispartonfN).type==-out.at(prevpartonfN).type) {
	IFfgqq=true;
	thispartonfN-=2;
	prevpartonfN-=2;
	} else if ((out.at(thispartonfN).type==pid::u_bar&&out.at(prevpartonfN).type==pid::d)
	\|\|(out.at(thispartonfN).type==pid::c_bar&&out.at(prevpartonfN).type==pid::s)
	\|\|(out.at(thispartonfN).type==pid::d&&out.at(prevpartonfN).type==pid::u_bar)
	\|\|(out.at(thispartonfN).type==pid::s&&out.at(prevpartonfN).type==pid::c_bar)){
	// if in SU(2)-up and out SU(2)-down then gqqWp-type
	IFfgqqWp=true;
	thispartonfN-=2;
	prevpartonfN-=2;
	} else if ((out.at(thispartonfN).type==pid::u&&out.at(prevpartonfN).type==pid::d_bar)
	\|\|(out.at(thispartonfN).type==pid::c&&out.at(prevpartonfN).type==pid::s_bar)
	\|\|(out.at(thispartonfN).type==pid::d_bar&&out.at(prevpartonfN).type==pid::u)
	\|\|(out.at(thispartonfN).type==pid::s_bar&&out.at(prevpartonfN).type==pid::c)){
	// if in SU(2)-down and out SU(2)-up then gqqWm-type
	IFfgqqWm=true;
	thispartonfN-=2;
	prevpartonfN-=2;
	}
	}
	} //impact factors done
	// if we reached this far and thispartonf==end(out),
	// then it is not a valid all-order state
	#ifndef NDEBUG
	std::cerr << "IFfLL : "<<IFfLL<<std::endl;
	std::cerr << "IFfuno : "<<IFfuno<<std::endl;
	#endif
	if (thispartonfN==outN-1) return EventType::FixedOrder;
	// Now check the remaining middle partons
	// the first (in rapidity) parton in the forward impact factor
	long unsigned int firstforwardIFpartonN = thispartonfN+1;
	// allow only one qqbar(+Wp/Wm) insertion
	while ((out.at(thispartonN).type==pid::gluon)&&(thispartonN!=firstforwardIFpartonN)) {
	++thispartonN;
	++nextpartonN;
	#ifndef NDEBUG
	std::cerr<<"one gluon found\n";
	#endif
	}
	#ifndef NDEBUG
	std::cerr<<"thisparton : "<<out.at(thispartonN).type<<std::endl<<"firstforwardIFpartonN : "<<out.at(firstforwardIFpartonN).type<<std::endl;
	std::cerr<< "(thispartonN!=firstforwardIFpartonN) : "<<(thispartonN!=firstforwardIFpartonN) <<std::endl;
	// std::cerr<< "(thisparton==firstforwardIFparton) : "<<((thisparton)==(firstforwardIFparton)) <<std::endl;
	#endif
	if (thispartonN!=firstforwardIFpartonN) {
	// is this a qqbar-pair?
	if (out.at(thispartonN).type==-out.at(nextpartonN).type) {
	Midqqbar=true;
	thispartonN+=2;
	nextpartonN+=2;
	} else if ((out.at(thispartonN).type==pid::u_bar&&out.at(nextpartonN).type==pid::d)
	\|\|(out.at(thispartonN).type==pid::c_bar&&out.at(nextpartonN).type==pid::s)
	\|\|(out.at(thispartonN).type==pid::d&&out.at(nextpartonN).type==pid::u_bar)
	\|\|(out.at(thispartonN).type==pid::s&&out.at(nextpartonN).type==pid::c_bar)) {
	// if in SU(2)-up and out SU(2)-down then gqqWp-type
	MidqqbarWp=true;
	thispartonN+=2;
	nextpartonN+=2;
	} else if ((out.at(thispartonN).type==pid::u&&out.at(nextpartonN).type==pid::d_bar)
	\|\|(out.at(thispartonN).type==pid::c&&out.at(nextpartonN).type==pid::s_bar)
	\|\|(out.at(thispartonN).type==pid::d_bar&&out.at(nextpartonN).type==pid::u)
	\|\|(out.at(thispartonN).type==pid::s_bar&&out.at(nextpartonN).type==pid::c)) {
	// if in SU(2)-down and out SU(2)-up then gqqWm-type
	MidqqbarWm=true;
	thispartonN+=2;
	nextpartonN+=2;
	}

	//continue checking
	while ((out.at(thispartonN).type==pid::gluon)&&(thispartonN!=firstforwardIFpartonN)) {
	++thispartonN;
	++nextpartonN;
	}
	}
	// are we at the end of the partons? if not, this is not resummable
	if (thispartonN!=firstforwardIFpartonN) return EventType::FixedOrder;

	#ifndef NDEBUG
	std::cerr<<!(Midqqbar\|\|MidqqbarWm\|\|MidqqbarWp)<<std::endl;
	#endif
	// Checks using the non-partonic momenta
	// auto const boson { std::find_if(out.cbegin(), out.cend(), [](Particle const & p){return is_AZHW(p);}) }
	int WpN(0),WmN(0),HN(0),ZN(0);
	// if(boson!= out.cend(){ // if some bosons were found
	// // count number of W+, W-, H, Z/gamma
	// for (auto bitr=boson.begin(); bitr!=boson.end(); bitr++) {
	// if (bitr->type==24) ++WpN;
	// else if (bitr->type==-24) ++WmN;
	// else if (bitr->type==25) ++HN;
	// else if (bitr->type==23) ++ZN;
	// }
	// }

	// veto Higgs+uno with Higgs mixing in uno particles
	// veto Higgs between central qqbar pair


	// Check whether the right number of Ws are present
	// count the numbers of Wp Wm required by the identification of partons in the events
	- int NWpR=(IFbLLWp?true:1)+(IFfLLWp?true:1)+(MidqqbarWp?true:1)+(IFbunoWp?true:1)+(IFfunoWp?true:1);
	- int NWmR=(IFbLLWm?true:1)+(IFfLLWm?true:1)+(MidqqbarWm?true:1)+(IFbunoWm?true:1)+(IFfunoWm?true:1);
	+ int NWpR=(IFbLLWp?1)+(IFfLLWp?1)+(MidqqbarWp?1)+(IFbunoWp?1)+(IFfunoWp?1);
	+ int NWmR=(IFbLLWm?1)+(IFfLLWm?1)+(MidqqbarWm?1)+(IFbunoWm?1)+(IFfunoWm?1);

	if (NWpR!=WpN) return EventType::FixedOrder;
	if (NWmR!=WmN) return EventType::FixedOrder;
	if (WmN>1\|\|WpN>1\|\|HN>1\|\|ZN>1) return EventType::FixedOrder;

	// return the event type
	if (!(Midqqbar\|\|MidqqbarWm\|\|MidqqbarWp)) { // if there are no mid qqbar
	if ((IFbLL\|\|IFbLLWm\|\|IFbLLWp)&&IFfLL) return EventType::FKL;
	if ((IFfLLWm\|\|IFfLLWp)&&IFbLL) return EventType::FKL;
	if ((IFbuno\|\|IFbunoWm\|\|IFbunoWp)&&IFfLL) return EventType::unordered_backward;
	if ((IFbuno)&&(IFfLLWp\|\|IFfLLWm)) return EventType::unordered_backward;
	if ((IFfuno\|\|IFfunoWm\|\|IFfunoWp)&&IFbLL) return EventType::unordered_forward;
	if ((IFfuno)&&(IFbLLWp\|\|IFbLLWm)) return EventType::unordered_forward;
	if ((IFbgqq\|\|IFbgqqWp\|\|IFbgqqWm)&&IFfLL) return EventType::extremal_qqxb;
	if ((IFbgqq)&&(IFfLLWp\|\|IFfLLWm)) return EventType::extremal_qqxb;
	if ((IFfgqq\|\|IFfgqqWp\|\|IFfgqqWm)&&IFbLL) return EventType::extremal_qqxf;
	if ((IFfgqq)&&(IFbLLWp\|\|IFbLLWm)) return EventType::extremal_qqxf;
	} else { // there is a middle qqbar pair
	if (IFbLL)
	if (IFfLL\|\|IFfLLWp\|\|IFfLLWm) return EventType::central_qqx;
	if (IFfLL)
	if (IFbLL\|\|IFbLLWp\|\|IFbLLWm) return EventType::central_qqx;
	}

	return EventType::FixedOrder;
	}

	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;
	if(is_Ex_qqxb(ev))
	return EventType::extremal_qqxb;
	if(is_Ex_qqxf(ev))
	return EventType::extremal_qqxf;
	if(is_Mid_qqx(ev))
	return EventType::central_qqx;
	return EventType::FixedOrder;
	}
	//@}

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

	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

	Event::EventData::EventData(LHEF::HEPEUP const & hepeup){
	parameters.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));
	}
	}

	Event::Event(
	UnclusteredEvent const & ev,
	fastjet::JetDefinition const & jet_def, double const min_jet_pt
	):
	Event( Event::EventData{
	ev.incoming, ev.outgoing, ev.decays,
	Parameters<EventParameters>{ev.central, ev.variations}
	}.cluster(jet_def, min_jet_pt) )
	{}

	//! @TODO remove in HEJ 2.2.0
	UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){
	Event::EventData const evData{hepeup};
	incoming = evData.incoming;
	outgoing = evData.outgoing;
	decays = evData.decays;
	central = evData.parameters.central;
	variations = evData.parameters.variations;
	}

	void Event::EventData::sort(){
	// sort particles
	std::sort(
	begin(incoming), end(incoming),
	[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
	);

	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();
	});
	outgoing.clear();
	outgoing.reserve(old_outgoing.size());
	for(size_t i: idx) {
	outgoing.emplace_back(std::move(old_outgoing[i]));
	}

	// find decays again
	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())
	decays.emplace(i, std::move(decay->second));
	}
	assert(old_decays.size() == decays.size());
	}
	}

	namespace {
	Particle reconstruct_boson(std::vector<Particle> const & leptons) {
	HEJ::Particle decayed_boson;
	decayed_boson.p = leptons[0].p + leptons[1].p;
	const int pidsum = leptons[0].type + leptons[1].type;
	if(pidsum == +1) {
	assert(is_antilepton(leptons[0]));
	if(is_antineutrino(leptons[0])) {
	throw HEJ::not_implemented{"lepton-flavour violating final state"};
	}
	assert(is_neutrino(leptons[1]));
	// charged antilepton + neutrino means we had a W+
	decayed_boson.type = HEJ::pid::Wp;
	}
	else if(pidsum == -1) {
	assert(is_antilepton(leptons[0]));
	if(is_neutrino(leptons[1])) {
	throw HEJ::not_implemented{"lepton-flavour violating final state"};
	}
	assert(is_antineutrino(leptons[0]));
	// charged lepton + antineutrino means we had a W-
	decayed_boson.type = HEJ::pid::Wm;
	}
	else {
	throw HEJ::not_implemented{
	"final state with leptons "
	+ HEJ::name(leptons[0].type)
	+ " and "
	+ HEJ::name(leptons[1].type)
	};
	}
	return decayed_boson;
	}
	}

	void HEJ::Event::EventData::reconstruct_intermediate() {
	const auto begin_leptons = std::partition(
	begin(outgoing), end(outgoing),
	[](HEJ::Particle const & p) {return !HEJ::is_anylepton(p);}
	);
	if(begin_leptons == end(outgoing)) return;
	assert(is_anylepton(*begin_leptons));
	std::vector<HEJ::Particle> leptons(begin_leptons, end(outgoing));
	outgoing.erase(begin_leptons, end(outgoing));
	if(leptons.size() != 2) {
	throw HEJ::not_implemented{"Final states with one or more than two leptons"};
	}
	std::sort(
	begin(leptons), end(leptons),
	[](HEJ::Particle const & p0, HEJ::Particle const & p1) {
	return p0.type < p1.type;
	}
	);
	outgoing.emplace_back(reconstruct_boson(leptons));
	decays.emplace(outgoing.size()-1, std::move(leptons));
	}

	Event Event::EventData::cluster(
	fastjet::JetDefinition const & jet_def, double const min_jet_pt
	){
	sort();
	Event ev{ std::move(incoming), std::move(outgoing), std::move(decays),
	std::move(parameters),
	jet_def, min_jet_pt
	};
	assert(std::is_sorted(begin(ev.outgoing_), end(ev.outgoing_),
	rapidity_less{}));
	ev.type_ = classify(ev);
	if (classify_new(ev)!=ev.type_) {
	std::cout<< classify_new(ev) << " "<<ev.type_<<std::endl;
	std::cout<<ev;
	}
	assert (classify_new(ev)==ev.type_);
	return ev;
	}

	Event::Event(
	std::array<Particle, 2> && incoming,
	std::vector<Particle> && outgoing,
	std::unordered_map<size_t, std::vector<Particle>> && decays,
	Parameters<EventParameters> && parameters,
	fastjet::JetDefinition const & jet_def,
	double const min_jet_pt
	): incoming_{std::move(incoming)},
	outgoing_{std::move(outgoing)},
	decays_{std::move(decays)},
	parameters_{std::move(parameters)},
	cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
	min_jet_pt_{min_jet_pt}
	{
	jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_));
	}

	namespace {
	void connect_incoming(Particle & in, int & colour, int & anti_colour){
	in.colour = std::make_pair(anti_colour, colour);
	// gluon
	if(in.type == pid::gluon)
	return;
	if(in.type > 0){
	// quark
	assert(is_quark(in));
	in.colour->second = 0;
	colour*=-1;
	return;
	}
	// anti-quark
	assert(is_antiquark(in));
	in.colour->first = 0;
	anti_colour*=-1;
	return;
	}
	}

	bool Event::generate_colours(RNG & ran){
	// generate only for HEJ events
	if(!event_type::is_HEJ(type()))
	return false;
	assert(std::is_sorted(
	begin(outgoing()), end(outgoing()), rapidity_less{}));
	assert(incoming()[0].pz() < incoming()[1].pz());

	// positive (anti-)colour -> can connect
	// negative (anti-)colour -> not available/used up by (anti-)quark
	int colour = COLOUR_OFFSET;
	int anti_colour = colour+1;
	// initialise first
	connect_incoming(incoming_[0], colour, anti_colour);

	for(auto & part: outgoing_){
	assert(colour>0 \|\| anti_colour>0);
	if(part.type == ParticleID::gluon){
	// gluon
	if(colour>0 && anti_colour>0){
	// on g line => connect to colour OR anti-colour (random)
	if(ran.flat() < 0.5){
	part.colour = std::make_pair(colour+2,colour);
	colour+=2;
	} else {
	part.colour = std::make_pair(anti_colour, anti_colour+2);
	anti_colour+=2;
	}
	} else if(colour > 0){
	// on q line => connect to available colour
	part.colour = std::make_pair(colour+2, colour);
	colour+=2;
	} else {
	assert(colour<0 && anti_colour>0);
	// on qx line => connect to available anti-colour
	part.colour = std::make_pair(anti_colour, anti_colour+2);
	anti_colour+=2;
	}
	} else if(is_quark(part)) {
	// quark
	assert(anti_colour>0);
	if(colour>0){
	// on g line => connect and remove anti-colour
	part.colour = std::make_pair(anti_colour, 0);
	anti_colour+=2;
	anti_colour*=-1;
	} else {
	// on qx line => new colour
	colour*=-1;
	part.colour = std::make_pair(colour, 0);
	}
	} else if(is_antiquark(part)) {
	// anti-quark
	assert(colour>0);
	if(anti_colour>0){
	// on g line => connect and remove colour
	part.colour = std::make_pair(0, colour);
	colour+=2;
	colour*=-1;
	} else {
	// on q line => new anti-colour
	anti_colour*=-1;
	part.colour = std::make_pair(0, anti_colour);
	}
	}
	// else not a parton
	}
	// Connect last
	connect_incoming(incoming_[1], anti_colour, colour);
	return true;
	} // generate_colours

	namespace {
	void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){
	const std::streamsize orig_prec = os.precision();
	os <<std::scientific<<std::setprecision(6) << "["
	<<std::setw(13)<<std::right<< part.px() << ", "
	<<std::setw(13)<<std::right<< part.py() << ", "
	<<std::setw(13)<<std::right<< part.pz() << ", "
	<<std::setw(13)<<std::right<< part.E() << "]"<< std::fixed;
	os.precision(orig_prec);
	}
	}

	std::ostream& operator<<(std::ostream & os, Event const & ev){
	const std::streamsize orig_prec = os.precision();
	os <<std::setprecision(4)<<std::fixed;
	std::cout << "########## " << event_type::name(ev.type()) << " ##########" << std::endl;
	std::cout << "Incoming particles:\n";
	for(auto const & in: ev.incoming()){
	std::cout <<std::setw(3)<< in.type << ": ";
	print_momentum(os, in.p);
	std::cout << std::endl;
	}
	std::cout << "\nOutgoing particles: " << ev.outgoing().size() << "\n";
	for(auto const & out: ev.outgoing()){
	std::cout <<std::setw(3)<< out.type << ": ";
	print_momentum(os, out.p);
	std::cout << " => rapidity="
	<<std::setw(7)<<std::right<< out.rapidity() << std::endl;
	}
	std::cout << "\nForming Jets: " << ev.jets().size() << "\n";
	for(auto const & jet: ev.jets()){
	print_momentum(os, jet);
	std::cout << " => rapidity="
	<<std::setw(7)<<std::right<< jet.rapidity() << std::endl;
	}
	os << std::defaultfloat;
	os.precision(orig_prec);
	return os;
	}

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

	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;
	result.IDPRUP = event.type(); // event type
	// the following entries are pretty much meaningless
	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;

	result.IDUP.reserve(num_particles); // PID
	result.ISTUP.reserve(num_particles); // status (in, out, decay)
	result.PUP.reserve(num_particles); // momentum
	result.MOTHUP.reserve(num_particles); // index mother particle
	result.ICOLUP.reserve(num_particles); // colour
	// incoming
	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);
	assert(in.colour);
	result.ICOLUP.emplace_back(*in.colour);
	}
	// outgoing
	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);
	if(out.colour)
	result.ICOLUP.emplace_back(*out.colour);
	else{
	assert(is_AWZH_boson(out));
	result.ICOLUP.emplace_back(std::make_pair(0,0));
	}
	}
	// decays
	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;
	}

	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:17 PM (1 d, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3798535
Default Alt Text: (49 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions