Event.cc
No OneTemporary
Actions

Size

39 KB

Referenced Files

None

Subscribers

None

Event.cc
View Options

	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019-2020
	* \copyright GPLv2 or later
	*/
	#include "HEJ/Event.hh"

	#include <algorithm>
	#include <cassert>
	#include <cstdlib>
	#include <iomanip>
	#include <iterator>
	#include <memory>
	#include <numeric>
	#include <ostream>
	#include <string>
	#include <utility>

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

	#include "LHEF/LHEF.h"

	#include "HEJ/Constants.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/RNG.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/optional.hh"

	namespace HEJ {

	namespace {
	using std::size_t;

	//! LHE status codes
	namespace lhe_status {
	enum Status: int {
	in = -1,
	decay = 2,
	out = 1,
	};
	}
	using LHE_Status = lhe_status::Status;

	//! true if leptonic W decay
	bool valid_W_decay( int const w_type, // sign of W
	std::vector<Particle> const & decays
	){
	if(decays.size() != 2) // no 1->2 decay
	return false;
	const int pidsum = decays[0].type + decays[1].type;
	if( std::abs(pidsum) != 1 \|\| pidsum != w_type ) // correct charge
	return false;
	// leptonic decay (only check first, second follows from pidsum)
	if( w_type == 1 ) // W+
	return is_antilepton(decays[0]) \|\| is_neutrino(decays[0]);
	// W-
	return is_lepton(decays[0]) \|\| is_antineutrino(decays[0]);
	}

	//! true for Z decay to charged leptons
	bool valid_Z_decay(std::vector<Particle> const & decays){
	if(decays.size() != 2) // no 1->2 decay
	return false;
	const int pidsum = decays[0].type + decays[1].type;
	if( std::abs(pidsum) != 0 ) // correct charge
	return false;
	// leptonic decay (only check first, second follows from pidsum)
	return is_anylepton(decays[0]) && !is_anyneutrino(decays[0]);
	}

	/// @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(Event const & ev){
	std::vector<Particle> const & outgoing = ev.outgoing();
	if(ev.decays().size() > 1) // at most one decay
	return false;
	bool has_AWZH_boson = false;
	for( size_t i=0; i<outgoing.size(); ++i ){
	auto const & out{ outgoing[i] };
	if(is_AWZH_boson(out.type)){
	// at most one boson
	if(has_AWZH_boson) return false;
	has_AWZH_boson = true;

	// valid decay for W
	if(std::abs(out.type) == ParticleID::Wp){
	// exactly 1 decay of W
	if( ev.decays().size() != 1 \|\| ev.decays().cbegin()->first != i )
	return false;
	if( !valid_W_decay(out.type>0?+1:-1, ev.decays().cbegin()->second) )
	return false;
	}

	// valid decay for Z
	if(out.type == ParticleID::Z_photon_mix){
	// exactly 1 decay
	if( ev.decays().size() != 1 \|\| ev.decays().cbegin()->first != i )
	return false;
	if( !valid_Z_decay(ev.decays().cbegin()->second) )
	return false;
	}
	}
	else if(! is_parton(out.type)) return false;
	}
	return true;
	}

	/**
	* returns all EventTypes implemented in HEJ
	*/
	size_t implemented_types(std::vector<Particle> const & bosons){
	using namespace event_type;
	if(bosons.empty()) return FKL \| unob \| unof \| qqxexb \| qqxexf \| qqxmid;
	if(bosons.size()>1) return non_resummable; // multi boson
	switch (bosons[0].type) {
	case ParticleID::Wp:
	case ParticleID::Wm:
	return FKL \| unob \| unof \| qqxexb \| qqxexf \| qqxmid;
	case ParticleID::Z_photon_mix:
	return FKL \| unob \| unof;
	case ParticleID::h:
	return FKL \| unob \| unof;
	default:
	return non_resummable;
	}
	}

	/**
	* \brief function which determines if type change is consistent with Wp emission.
	* @param in incoming Particle id
	* @param out outgoing Particle id
	* @param qqx Current both incoming/both outgoing?
	*
	* \see is_Wm_Change
	*/
	bool is_Wp_Change(ParticleID in, ParticleID out, bool qqx){
	using namespace pid;
	if(!qqx && (in==d_bar \|\| in==u \|\| in==s_bar \|\| in==c))
	return out == (in-1);
	if( qqx && (in==d \|\| in==u_bar \|\| in==s \|\| in==c_bar))
	return out == -(in+1);
	return false;
	}

	/**
	* \brief function which determines if type change is consistent with Wm emission.
	* @param in incoming Particle id
	* @param out outgoing Particle id
	* @param qqx Current both incoming/both outgoing?
	*
	* Ensures that change type of quark line is possible by a flavour changing
	* Wm emission. Allows checking of qqx currents also.
	*/
	bool is_Wm_Change(ParticleID in, ParticleID out, bool qqx){
	using namespace pid;
	if(!qqx && (in==d \|\| in==u_bar \|\| in==s \|\| in==c_bar))
	return out == (in+1);
	if( qqx && (in==d_bar \|\| in==u \|\| in==s_bar \|\| in==c))
	return out == -(in-1);
	return false;
	}

	/**
	* \brief checks if particle type remains same from incoming to outgoing
	* @param in incoming Particle
	* @param out outgoing Particle
	* @param qqx Current both incoming/outgoing?
	*/
	bool no_flavour_change(ParticleID in, ParticleID out, bool qqx){
	const int qqxCurrent = qqx?-1:1;
	if(std::abs(in)<=pid::top \|\| in==pid::gluon)
	return (in==out*qqxCurrent);
	return false;
	}

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

	/**
	* \brief check if we have a valid Impact factor
	* @param in incoming Particle
	* @param out outgoing Particle
	* @param qqx Current both incoming/outgoing?
	* @param W_change returns +1 if Wp, -1 if Wm, else 0
	*/
	bool is_valid_impact_factor(
	ParticleID in, ParticleID out, bool qqx, int & W_change
	){
	if( no_flavour_change(in, out, qqx) ){
	return true;
	}
	if( is_Wp_Change(in, out, qqx) ) {
	W_change+=1;
	return true;
	}
	if( is_Wm_Change(in, out, qqx) ) {
	W_change-=1;
	return true;
	}
	return false;
	}

	//! Returns all possible classifications from the impact factors
	// the beginning points are changed s.t. after the the classification they
	// point to the beginning of the (potential) FKL chain
	// sets W_change: + if Wp change
	// 0 if no change
	// - if Wm change
	// This function can be used with forward & backwards iterators
	template<class OutIterator>
	size_t possible_impact_factors(
	ParticleID incoming_id, // incoming
	OutIterator & begin_out, OutIterator const & end_out, // outgoing
	int & W_change, std::vector<Particle> const & boson,
	bool const backward // backward?
	){
	using namespace event_type;
	assert(boson.size() < 2);
	// keep track of all states that we don't test
	size_t not_tested = qqxmid;
	if(backward)
	not_tested \|= unof \| qqxexf;
	else
	not_tested \|= unob \| qqxexb;

	// Is this LL current?
	if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){
	++begin_out;
	return not_tested \| FKL;
	}

	// or NLL current?
	// -> needs two partons in two different jets
	if( std::distance(begin_out, end_out)>=2
	){
	auto next = std::next(begin_out);
	// Is this unordered emisson?
	if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){
	if( is_valid_impact_factor(
	incoming_id, next->type, false, W_change )
	){
	// veto Higgs inside uno
	assert(next!=end_out);
	if( !boson.empty() && boson.front().type == ParticleID::h
	){
	if( (backward && boson.front().rapidity() < next->rapidity())
	\|\|(!backward && boson.front().rapidity() > next->rapidity()))
	return non_resummable;
	}
	begin_out = std::next(next);
	return not_tested \| (backward?unob:unof);
	}
	}
	// Is this QQbar?
	else if( incoming_id==pid::gluon ){
	if( is_valid_impact_factor(
	begin_out->type, next->type, true, W_change )
	){
	// veto Higgs inside qqx
	assert(next!=end_out);
	if( !boson.empty() && boson.front().type == ParticleID::h
	){
	if( (backward && boson.front().rapidity() < next->rapidity())
	\|\|(!backward && boson.front().rapidity() > next->rapidity()))
	return non_resummable;
	}
	begin_out = std::next(next);
	return not_tested \| (backward?qqxexb:qqxexf);
	}
	}
	}
	return non_resummable;
	}

	//! Returns all possible classifications from central emissions
	// the beginning points are changed s.t. after the the classification they
	// point to the end of the emission chain
	// sets W_change: + if Wp change
	// 0 if no change
	// - if Wm change
	template<class OutIterator>
	size_t possible_central(
	OutIterator & begin_out, OutIterator const & end_out,
	int & W_change, std::vector<Particle> const & boson
	){
	using namespace event_type;
	assert(boson.size() < 2);
	// if we already passed the central chain,
	// then it is not a valid all-order state
	if(std::distance(begin_out, end_out) < 0) return non_resummable;
	// keep track of all states that we don't test
	size_t possible = unob \| unof
	\| qqxexb \| qqxexf;

	// Find the first non-gluon/non-FKL
	while( (begin_out->type==pid::gluon) && (begin_out!=end_out) ){
	++begin_out;
	}
	// end of chain -> FKL
	if( begin_out==end_out ){
	return possible \| FKL;
	}

	// is this a qqbar-pair?
	// needs two partons in two separate jets
	auto next = std::next(begin_out);
	if( is_valid_impact_factor(
	begin_out->type, next->type, true, W_change )
	){
	// veto Higgs inside qqx
	if( !boson.empty() && boson.front().type == ParticleID::h
	&& boson.front().rapidity() > begin_out->rapidity()
	&& boson.front().rapidity() < next->rapidity()
	){
	return non_resummable;
	}
	begin_out = std::next(next);
	// remaining chain should be pure gluon/FKL
	for(; begin_out!=end_out; ++begin_out){
	if(begin_out->type != pid::gluon) return non_resummable;
	}
	return possible \| qqxmid;
	}
	return non_resummable;
	}

	/**
	* \brief Checks for all event types
	* @param ev Event
	* @returns Event Type
	*
	*/
	event_type::EventType classify(Event const & ev){
	using namespace event_type;
	if(! has_2_jets(ev))
	return no_2_jets;
	// currently we can't handle multiple boson states in the ME. So they are
	// considered "bad_final_state" even though the "classify" could work with
	// them.
	if(! final_state_ok(ev))
	return bad_final_state;

	// initialise variables
	auto const & in = ev.incoming();

	// range for current checks
	auto begin_out{ev.cbegin_partons()};
	auto end_out{ev.crbegin_partons()};

	assert(std::distance(begin(in), end(in)) == 2);
	assert(std::distance(begin_out, end_out.base()) >= 2);
	assert(std::is_sorted(begin_out, end_out.base(), rapidity_less{}));

	auto const boson{ filter_AWZH_bosons(ev.outgoing()) };
	// we only allow one boson through final_state_ok
	assert(boson.size()<=1);

	// keep track of potential W couplings, at the end the sum should be 0
	int remaining_Wp = 0;
	int remaining_Wm = 0;
	if(!boson.empty() && std::abs(boson.front().type) == ParticleID::Wp ){
	if(boson.front().type>0) ++remaining_Wp;
	else ++remaining_Wm;
	}
	int W_change = 0;

	size_t final_type = ~(no_2_jets \| bad_final_state);

	// check forward impact factor
	final_type &= possible_impact_factors(
	in.front().type,
	begin_out, end_out.base(),
	W_change, boson, true );
	if( final_type == non_resummable )
	return non_resummable;
	if(W_change>0) remaining_Wp-=W_change;
	else if(W_change<0) remaining_Wm+=W_change;
	W_change = 0;

	// check backward impact factor
	final_type &= possible_impact_factors(
	in.back().type,
	end_out, std::make_reverse_iterator(begin_out),
	W_change, boson, false );
	if( final_type == non_resummable )
	return non_resummable;
	if(W_change>0) remaining_Wp-=W_change;
	else if(W_change<0) remaining_Wm+=W_change;
	W_change = 0;

	// check central emissions
	final_type &= possible_central(
	begin_out, end_out.base(), W_change, boson );
	if( final_type == non_resummable )
	return non_resummable;
	if(W_change>0) remaining_Wp-=W_change;
	else if(W_change<0) remaining_Wm+=W_change;

	// Check whether the right number of Ws are present
	if( remaining_Wp != 0 \|\| remaining_Wm != 0 ) return non_resummable;

	// result has to be unique
	if( (final_type & (final_type-1)) != 0) return non_resummable;

	// check that each sub processes is implemented
	// (has to be done at the end)
	if( (final_type & ~implemented_types(boson)) != 0 )
	return non_resummable;

	return static_cast<EventType>(final_type);
	}
	//@}

	Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){
	auto id = static_cast<ParticleID>(hepeup.IDUP[i]);
	auto colour = is_parton(id)?hepeup.ICOLUP[i]:optional<Colour>();
	return { id,
	{ hepeup.PUP[i][0], hepeup.PUP[i][1],
	hepeup.PUP[i][2], hepeup.PUP[i][3] },
	colour
	};
	}

	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

	Event::EventData::EventData(LHEF::HEPEUP const & hepeup){
	parameters.central = EventParameters{
	hepeup.scales.mur, hepeup.scales.muf, hepeup.XWGTUP
	};
	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] == LHE_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 const & o1, Particle const & 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) {
	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 not_implemented{"lepton-flavour violating final state"};
	}
	assert(is_neutrino(leptons[1]));
	// charged antilepton + neutrino means we had a W+
	decayed_boson.type = pid::Wp;
	}
	else if(pidsum == -1) {
	assert(is_antilepton(leptons[0]));
	if(is_neutrino(leptons[1])) {
	throw not_implemented{"lepton-flavour violating final state"};
	}
	assert(is_antineutrino(leptons[0]));
	// charged lepton + antineutrino means we had a W-
	decayed_boson.type = pid::Wm;
	}
	else if(pidsum == 0) {
	assert(is_anylepton(leptons[0]));
	if(is_anyneutrino(leptons[0])) {
	throw not_implemented{"final state with two neutrinos"};
	}
	// charged lepton-antilepton pair means we had a Z/photon
	decayed_boson.type = pid::Z_photon_mix;
	}
	else {
	throw not_implemented{
	"final state with leptons "
	+ name(leptons[0].type)
	+ " and "
	+ name(leptons[1].type)
	};
	}
	return decayed_boson;
	}
	} // namespace

	void Event::EventData::reconstruct_intermediate() {
	const auto begin_leptons = std::partition(
	begin(outgoing), end(outgoing),
	[](Particle const & p) {return !is_anylepton(p);}
	);
	// We can only reconstruct FS with 2 leptons
	if(std::distance(begin_leptons, end(outgoing)) != 2) return;
	std::vector<Particle> leptons(begin_leptons, end(outgoing));
	std::sort(
	begin(leptons), end(leptons),
	[](Particle const & p0, Particle const & p1) {
	assert(is_anylepton(p0) && is_anylepton(p1));
	return p0.type < p1.type;
	}
	);
	// `reconstruct_boson` can throw, it should therefore be called before
	// changing `outgoing` to allow the user to recover the original EventData
	auto boson = reconstruct_boson(leptons);
	outgoing.erase(begin_leptons, end(outgoing));
	outgoing.emplace_back(std::move(boson));
	decays.emplace(outgoing.size()-1, std::move(leptons));
	}

	Event Event::EventData::cluster(
	fastjet::JetDefinition const & jet_def, double const min_jet_pt
	){
	sort();
	return Event{ std::move(incoming), std::move(outgoing), std::move(decays),
	std::move(parameters),
	jet_def, min_jet_pt
	};
	}

	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_));
	assert(std::is_sorted(begin(outgoing_), end(outgoing_),
	rapidity_less{}));
	type_ = classify(*this);
	}

	namespace {
	//! check that Particles have a reasonable colour
	bool correct_colour(Particle const & part){
	ParticleID id{ part.type };
	if(!is_parton(id))
	return !part.colour;

	if(!part.colour)
	return false;

	Colour const & col{ *part.colour };
	if(is_quark(id))
	return col.first != 0 && col.second == 0;
	if(is_antiquark(id))
	return col.first == 0 && col.second != 0;
	assert(id==ParticleID::gluon);
	return col.first != 0 && col.second != 0 && col.first != col.second;
	}

	//! Connect parton to t-channel colour line & update the line
	//! returns false if connection not possible
	template<class OutIterator>
	bool try_connect_t(OutIterator const & it_part, Colour & line_colour){
	if( line_colour.first == it_part->colour->second ){
	line_colour.first = it_part->colour->first;
	return true;
	}
	if( line_colour.second == it_part->colour->first ){
	line_colour.second = it_part->colour->second;
	return true;
	}
	return false;
	}

	//! Connect parton to u-channel colour line & update the line
	//! returns false if connection not possible
	template<class OutIterator>
	bool try_connect_u(OutIterator & it_part, Colour & line_colour){
	auto it_next = std::next(it_part);
	if( try_connect_t(it_next, line_colour)
	&& try_connect_t(it_part, line_colour)
	){
	it_part=it_next;
	return true;
	}
	return false;
	}
	} // namespace

	bool Event::is_leading_colour() const {
	if( !correct_colour(incoming()[0]) \|\| !correct_colour(incoming()[1]) )
	return false;

	Colour line_colour = *incoming()[0].colour;
	std::swap(line_colour.first, line_colour.second);

	// reasonable colour
	if(!std::all_of(outgoing().cbegin(), outgoing().cend(), correct_colour))
	return false;

	for(auto it_part = cbegin_partons(); it_part!=cend_partons(); ++it_part){

	switch (type()) {
	case event_type::FKL:
	if( !try_connect_t(it_part, line_colour) )
	return false;
	break;
	case event_type::unob:
	case event_type::qqxexb: {
	if( !try_connect_t(it_part, line_colour)
	// u-channel only allowed at impact factor
	&& (std::distance(cbegin_partons(), it_part)!=0
	\|\| !try_connect_u(it_part, line_colour)))
	return false;
	break;
	}
	case event_type::unof:
	case event_type::qqxexf: {
	if( !try_connect_t(it_part, line_colour)
	// u-channel only allowed at impact factor
	&& (std::distance(it_part, cend_partons())!=2
	\|\| !try_connect_u(it_part, line_colour)))
	return false;
	break;
	}
	case event_type::qqxmid:{
	auto it_next = std::next(it_part);
	if( !try_connect_t(it_part, line_colour)
	// u-channel only allowed at qqx/qxq pair
	&& ( ( !(is_quark(it_part) && is_antiquark(it_next))
	&& !(is_antiquark(it_part) && is_quark(it_next)))
	\|\| !try_connect_u(it_part, line_colour))
	)
	return false;
	break;
	}
	default:
	throw std::logic_error{"unreachable"};
	}

	// no colour singlet exchange/disconnected diagram
	if(line_colour.first == line_colour.second)
	return false;
	}

	return (incoming()[1].colour->first == line_colour.first)
	&& (incoming()[1].colour->second == line_colour.second);
	}

	namespace {
	//! connect incoming Particle to colour flow
	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;
	}

	//! connect outgoing Particle to t-channel colour flow
	template<class OutIterator>
	void connect_tchannel(
	OutIterator & it_part, int & colour, int & anti_colour, RNG & ran
	){
	assert(colour>0 \|\| anti_colour>0);
	if(it_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){
	it_part->colour = std::make_pair(colour+2,colour);
	colour+=2;
	} else {
	it_part->colour = std::make_pair(anti_colour, anti_colour+2);
	anti_colour+=2;
	}
	} else if(colour > 0){
	// on q line => connect to available colour
	it_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
	it_part->colour = std::make_pair(anti_colour, anti_colour+2);
	anti_colour+=2;
	}
	} else if(is_quark(*it_part)) {
	// quark
	assert(anti_colour>0);
	if(colour>0){
	// on g line => connect and remove anti-colour
	it_part->colour = std::make_pair(anti_colour, 0);
	anti_colour+=2;
	anti_colour*=-1;
	} else {
	// on qx line => new colour
	colour*=-1;
	it_part->colour = std::make_pair(colour, 0);
	}
	} else if(is_antiquark(*it_part)) {
	// anti-quark
	assert(colour>0);
	if(anti_colour>0){
	// on g line => connect and remove colour
	it_part->colour = std::make_pair(0, colour);
	colour+=2;
	colour*=-1;
	} else {
	// on q line => new anti-colour
	anti_colour*=-1;
	it_part->colour = std::make_pair(0, anti_colour);
	}
	} else { // not a parton
	assert(!is_parton(*it_part));
	it_part->colour = {};
	}
	}

	//! connect to t- or u-channel colour flow
	template<class OutIterator>
	void connect_utchannel(
	OutIterator & it_part, int & colour, int & anti_colour, RNG & ran
	){
	OutIterator it_first = it_part++;
	if(ran.flat()<.5) {// t-channel
	connect_tchannel(it_first, colour, anti_colour, ran);
	connect_tchannel(it_part, colour, anti_colour, ran);
	}
	else { // u-channel
	connect_tchannel(it_part, colour, anti_colour, ran);
	connect_tchannel(it_first, colour, anti_colour, ran);
	}
	}
	} // namespace

	bool Event::generate_colours(RNG & ran){
	// generate only for HEJ events
	if(!event_type::is_resummable(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);

	// reset outgoing colours
	std::for_each(outgoing_.begin(), outgoing_.end(),
	[](Particle & part){ part.colour = {};});

	for(auto it_part = begin_partons(); it_part!=end_partons(); ++it_part){
	switch (type()) {
	// subleading can connect to t- or u-channel
	case event_type::unob:
	case event_type::qqxexb: {
	if( std::distance(begin_partons(), it_part)==0)
	connect_utchannel(it_part, colour, anti_colour, ran);
	else
	connect_tchannel(it_part, colour, anti_colour, ran);
	break;
	}
	case event_type::unof:
	case event_type::qqxexf: {
	if( std::distance(it_part, end_partons())==2)
	connect_utchannel(it_part, colour, anti_colour, ran);
	else
	connect_tchannel(it_part, colour, anti_colour, ran);
	break;
	}
	case event_type::qqxmid:{
	auto it_next = std::next(it_part);
	if( std::distance(begin_partons(), it_part)>0
	&& std::distance(it_part, end_partons())>2
	&& ( (is_quark(it_part) && is_antiquark(it_next))
	\|\| (is_antiquark(it_part) && is_quark(it_next)) )
	)
	connect_utchannel(it_part, colour, anti_colour, ran);
	else
	connect_tchannel(it_part, colour, anti_colour, ran);
	break;
	}
	default: // rest has to be t-channel
	connect_tchannel(it_part, colour, anti_colour, ran);
	}
	}
	// Connect last
	connect_incoming(incoming_[1], anti_colour, colour);
	assert(is_leading_colour());
	return true;
	} // generate_colours

	namespace {
	bool valid_parton(
	std::vector<fastjet::PseudoJet> const & jets,
	Particle const & parton, int const idx,
	double const soft_pt_regulator, double const min_extparton_pt
	){
	// TODO code overlap with PhaseSpacePoint::pass_extremal_cuts
	if(min_extparton_pt > parton.pt()) return false;
	if(idx<0) return false;
	assert(static_cast<int>(jets.size())>=idx);
	auto const & jet{ jets[idx] };
	return (parton.p - jet).pt()/jet.pt() <= soft_pt_regulator;
	}
	} // namespace

	// this should work with multiple types
	bool Event::valid_hej_state(double const soft_pt_regulator,
	double const min_pt
	) const {
	using namespace event_type;
	if(!is_resummable(type()))
	return false;

	auto const & jet_idx{ particle_jet_indices() };
	auto idx_begin{ jet_idx.cbegin() };
	auto idx_end{ jet_idx.crbegin() };

	auto part_begin{ cbegin_partons() };
	auto part_end{ crbegin_partons() };

	// always seperate extremal jets
	if( !valid_parton(jets(), part_begin, idx_begin,
	soft_pt_regulator, min_pt) )
	return false;
	++part_begin;
	++idx_begin;
	if( !valid_parton(jets(), part_end, idx_end,
	soft_pt_regulator, min_pt) )
	return false;
	++part_end;
	++idx_end;

	// unob -> second parton in own jet
	if( type() & (unob \| qqxexb) ){
	if( !valid_parton(jets(), part_begin, idx_begin,
	soft_pt_regulator, min_pt) )
	return false;
	++part_begin;
	++idx_begin;
	}

	if( type() & (unof \| qqxexf) ){
	if( !valid_parton(jets(), part_end, idx_end,
	soft_pt_regulator, min_pt) )
	return false;
	++part_end;
	// ++idx_end; // last check, we don't need idx_end afterwards
	}

	if( type() & qqxmid ){
	// find qqx pair
	auto begin_qqx{ std::find_if( part_begin, part_end.base(),
	[](Particle const & part) -> bool {
	return part.type != ParticleID::gluon;
	}
	)};
	assert(begin_qqx != part_end.base());
	long int qqx_pos{ std::distance(part_begin, begin_qqx) };
	assert(qqx_pos >= 0);
	idx_begin+=qqx_pos;
	if( !( valid_parton(jets(), begin_qqx, idx_begin,
	soft_pt_regulator, min_pt)
	&& valid_parton(jets(), std::next(begin_qqx), std::next(idx_begin),
	soft_pt_regulator, min_pt)
	))
	return false;
	}
	return true;
	}

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

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

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

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

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

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

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

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

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

	void print_colour(std::ostream & os, optional<Colour> const & col){
	constexpr int width = 3;
	if(!col)
	os << "(no color)"; // American spelling for better alignment
	else
	os << "(" <<std::setw(width)<<std::right<< col->first
	<< ", " <<std::setw(width)<<std::right<< col->second << ")";
	}
	} // namespace

	std::ostream& operator<<(std::ostream & os, Event const & ev){
	constexpr int prec = 4;
	constexpr int wtype = 3; // width for types
	const std::streamsize orig_prec = os.precision();
	os <<std::setprecision(prec)<<std::fixed;
	os << "########## " << name(ev.type()) << " ##########" << std::endl;
	os << "Incoming particles:\n";
	for(auto const & in: ev.incoming()){
	os <<std::setw(wtype)<< in.type << ": ";
	print_colour(os, in.colour);
	os << " ";
	print_momentum(os, in.p);
	os << std::endl;
	}
	os << "\nOutgoing particles: " << ev.outgoing().size() << "\n";
	for(auto const & out: ev.outgoing()){
	os <<std::setw(wtype)<< out.type << ": ";
	print_colour(os, out.colour);
	os << " ";
	print_momentum(os, out.p);
	os << " => rapidity="
	<<std::setw(2*prec-1)<<std::right<< out.rapidity() << std::endl;
	}
	os << "\nForming Jets: " << ev.jets().size() << "\n";
	for(auto const & jet: ev.jets()){
	print_momentum(os, jet);
	os << " => rapidity="
	<<std::setw(2*prec-1)<<std::right<< jet.rapidity() << std::endl;
	}
	if(!ev.decays().empty() ){
	os << "\nDecays: " << ev.decays().size() << "\n";
	for(auto const & decay: ev.decays()){
	os <<std::setw(wtype)<< ev.outgoing()[decay.first].type
	<< " (" << decay.first << ") to:\n";
	for(auto const & out: decay.second){
	os <<" "<<std::setw(wtype)<< out.type << ": ";
	print_momentum(os, out.p);
	os << " => rapidity="
	<<std::setw(2*prec-1)<<std::right<< out.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;
	// the following entries are pretty much meaningless
	result.IDPRUP = event.type(); // event type
	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
	std::array<Particle, 2> incoming{ event.incoming() };
	// First incoming should be positive pz according to LHE standard
	// (or at least most (everyone?) do it this way, and Pythia assumes it)
	if(incoming[0].pz() < incoming[1].pz())
	std::swap(incoming[0], incoming[1]);
	for(Particle const & in: incoming){
	result.IDUP.emplace_back(in.type);
	result.ISTUP.emplace_back(LHE_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) != 0u
	?LHE_Status::decay
	:LHE_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{
	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(LHE_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;
	}

	} // namespace HEJ

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 6:11 AM (1 d, 1 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6543391
Default Alt Text: Event.cc (39 KB)

Event.ccNo OneTemporaryActions

Event.ccView Options

File Metadata

Event Timeline

Event.cc
No OneTemporary
Actions

Event.cc
View Options