No OneTemporary
Actions

Size

33 KB

Subscribers

None

View Options

	diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
	index a110609..a037a82 100644
	--- a/include/HEJ/Event.hh
	+++ b/include/HEJ/Event.hh
	@@ -1,203 +1,210 @@
	/** \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 "HEJ/RNG.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}
	{};
	};

	//! An event before jet clustering & classification
	struct UnclusteredEvent{
	//! Default Constructor
	UnclusteredEvent() = default;
	//! Constructor from LesHouches event information
	UnclusteredEvent(LHEF::HEPEUP const & hepeup);

	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 /
	- //! @brief Generate the particle colour leading in the MRK limit, see \cite Andersen:2011zd
	- //! @note This will overwrite all existing colours
	- void generate_colour_flow(HEJ::RNG &);
	- void sort(); /< Sort Particles in rapidity /
	};

	/** 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:
	//! Default Event Constructor
	Event() = default;
	//! Event Constructor adding jet clustering to an unclustered event
	Event(
	UnclusteredEvent ev,
	fastjet::JetDefinition const & jet_def, double min_jet_pt
	);

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

	//! The corresponding event before jet clustering
	UnclusteredEvent const & unclustered() const {
	return ev_;
	}

	//! Central parameter choice
	EventParameters const & central() const{
	return ev_.central;
	}

	//! Central parameter choice
	EventParameters & central(){
	return ev_.central;
	}

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

	//! Outgoing particles
	std::vector<Particle> const & outgoing() const{
	return ev_.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;
	}

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

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

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

	//! Parameter (scale) variation
	/**
	* @param i Index of the requested variation
	*/
	EventParameters & variations(size_t i){
	return ev_.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_;
	}

	+ //! Assign colours to each particle
	+ /**
	+ * @details Colour ordering is done according to leading colour in the MRK
	+ * limit, see \cite Andersen:2011zd. This only affects \ref
	+ * is_HEJ() "HEJ" configurations, all other \ref event_type
	+ * "EventTypes" will be ignored.
	+ * @note This overwrites all previous set colours.
	+ */
	+ void generate_colour_flow(HEJ::RNG &);
	+
	private:
	+ void sort(); /< Sort Particles in rapidity /
	UnclusteredEvent ev_;
	fastjet::ClusterSequence cs_;
	double min_jet_pt_;
	event_type::EventType type_;
	};

	//! 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 *);

	}
	diff --git a/include/HEJ/event_types.hh b/include/HEJ/event_types.hh
	index 1b145a8..df60728 100644
	--- a/include/HEJ/event_types.hh
	+++ b/include/HEJ/event_types.hh
	@@ -1,63 +1,65 @@
	/** \file
	* \brief Define different types of events.
	*
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/

	#pragma once

	#include "HEJ/utility.hh"

	namespace HEJ{

	//! Namespace for event types
	namespace event_type{
	//! Possible event types
	enum EventType: size_t{
	FKL, /*< FKL-type event /
	unordered_backward, /*< event with unordered backward emission /
	unordered_forward, /*< event with unordered forward emission /
	nonHEJ, /*< event configuration not covered by HEJ /
	no_2_jets, /*< event with less than two jets /
	bad_final_state, /*< event with an unsupported final state /
	unob = unordered_backward,
	unof = unordered_forward,
	first_type = FKL,
	last_type = bad_final_state
	};

	//! Event type names
	/**
	* For example, names[FKL] is the string "FKL"
	*/
	static constexpr auto names = make_array(
	"FKL",
	"unordered backward",
	"unordered forward",
	"nonHEJ",
	"no 2 jets",
	"bad final state"
	);

	+ //! Returns True for a HEJ \ref event_type::EventType "EventType"
	inline
	bool is_HEJ(EventType type) {
	switch(type) {
	case FKL:
	case unordered_backward:
	case unordered_forward:
	return true;
	default:
	return false;
	}
	}

	+ //! Returns True for an unordered \ref event_type::EventType "EventType"
	inline
	bool is_uno(EventType type) {
	return type == unordered_backward \|\| type == unordered_forward;
	}

	}

	}
	diff --git a/src/Event.cc b/src/Event.cc
	index 7df221a..31245b4 100644
	--- a/src/Event.cc
	+++ b/src/Event.cc
	@@ -1,486 +1,490 @@
	/**
	* \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/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;
	}

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

	UnclusteredEvent::UnclusteredEvent(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));
	}
	}

	namespace {
	int connect_incoming(Particle & part, int & colour, int & anti_colour){
	- // TODO something is wrong here
	part.colour = std::make_pair(anti_colour, colour);
	// gluon
	if(part.type == pid::gluon) {
	return 0;
	}
	// q
	if(part.type > 0){
	part.colour->second = 0;
	colour*=-1;
	return -1;
	}
	// qx
	part.colour->first = 0;
	anti_colour*=-1;
	return +1;
	}
	}

	- void UnclusteredEvent::generate_colour_flow(RNG & ran){
	- sort();
	+ void Event::generate_colour_flow(RNG & ran){
	+ // only sort HEJ events
	+ if(!event_type::is_HEJ(type()))
	+ return;
	+
	+ assert(std::is_sorted(
	+ begin(outgoing()), end(outgoing()), rapidity_less{}));
	+ assert(incoming()[0].pz() < incoming()[1].pz());

	// initialise first
	int colour = COLOUR_OFFSET;
	int anti_colour = colour+1;
	// -1 anti, 0 no, +1 colour
	- int is_quark_line = connect_incoming(incoming[0], colour, anti_colour);
	+ int is_quark_line = connect_incoming(ev_.incoming[0], colour, anti_colour);

	- for(auto & part: outgoing){
	+ for(auto & part: ev_.outgoing){
	if(part.type == ParticleID::gluon){
	// gluon
	if(is_quark_line == 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(is_quark_line > 0){
	// on q line => connect to available colour
	part.colour = std::make_pair(colour+2, colour);
	colour+=2;
	} else {
	// 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
	if(is_quark_line == 0){
	// on g line => connect and remove colour
	is_quark_line = 1;
	part.colour = std::make_pair(anti_colour, 0);
	anti_colour+=2;
	anti_colour*=-1;
	} else if(is_quark_line < 0){
	// on qx line => new colour
	is_quark_line = 0;
	colour*=-1;
	part.colour = std::make_pair(colour,0);

	} else {
	throw std::logic_error("Colour connection impossible: Expected anti-quark but found quark.");
	}
	} else if(is_antiquark(part)) {
	// anti-quark
	if(is_quark_line == 0){
	// on g line => connect and remove anti-colour
	is_quark_line = -1;
	part.colour = std::make_pair(0, colour);
	colour+=2;
	colour*=-1;
	} else if(is_quark_line > 0){
	// on q line => new anti-colour
	is_quark_line = 0;
	anti_colour*=-1;
	part.colour = std::make_pair(0,anti_colour);
	} else {
	throw std::logic_error("Colour connection impossible: Expected quark but found anti-quark.");
	}
	}
	}
	// Connect last
	- connect_incoming(incoming[1], anti_colour, colour);
	+ connect_incoming(ev_.incoming[1], anti_colour, colour);

	}

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

	//sort outgoing
	- if(std::is_sorted(begin(outgoing), end(outgoing), rapidity_less{}))
	+ if(std::is_sorted(begin(ev_.outgoing), end(ev_.outgoing), rapidity_less{}))
	return;
	- auto old_outgoing = std::move(outgoing);
	+ auto old_outgoing = std::move(ev_.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());
	+ ev_.outgoing.clear();
	for(size_t i: idx) {
	- outgoing.emplace_back(std::move(old_outgoing[i]));
	+ ev_.outgoing.emplace_back(std::move(old_outgoing[i]));
	}

	// find decays again
	- if(!decays.empty()){
	- auto old_decays = std::move(decays);
	- decays.clear();
	+ if(!ev_.decays.empty()){
	+ auto old_decays = std::move(ev_.decays);
	+ ev_.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));
	+ ev_.decays.emplace(i, std::move(decay->second));
	}
	- assert(old_decays.size() == decays.size());
	+ assert(old_decays.size() == ev_.decays.size());
	}

	}

	Event::Event(
	UnclusteredEvent ev,
	fastjet::JetDefinition const & jet_def, double min_jet_pt
	):
	ev_{std::move(ev)},
	cs_{to_PseudoJet(filter_partons(ev_.outgoing)), jet_def},
	min_jet_pt_{min_jet_pt}
	{
	// sort particles
	- ev_.sort();
	+ sort();

	// 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/t/test_colour_flow.cc b/t/test_colour_flow.cc
	index e76cd1f..44156c6 100644
	--- a/t/test_colour_flow.cc
	+++ b/t/test_colour_flow.cc
	@@ -1,205 +1,209 @@
	#include <stdexcept>
	#include <utility>

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

	/// biased RNG to connect always to colour
	class dum_rnd: public HEJ::DefaultRNG {
	public:
	dum_rnd() = default;
	double flat() override {
	return 0.;
	};
	};

	-void dump_event(HEJ::UnclusteredEvent const & ev){
	- for(auto const & in: ev.incoming){
	+void dump_event(HEJ::Event const & ev){
	+ for(auto const & in: ev.incoming()){
	std::cerr << "in type=" << in.type
	<< ", colour={" << (*in.colour).first
	<< ", " << (*in.colour).second << "}\n";
	}
	- for(auto const & out: ev.outgoing){
	+ for(auto const & out: ev.outgoing()){
	std::cerr << "out type=" << out.type << ", colour={";
	if(out.colour)
	std::cerr << (out.colour).first << ", " << (out.colour).second;
	else
	std::cerr << "non, non";
	std::cerr << "}\n";
	}
	}

	/// true if colour is allowed for particle
	bool correct_colour(HEJ::Particle const & part){
	if(HEJ::is_AWZH_boson(part) && !part.colour) return true;
	if(!part.colour) return false;
	int const colour = part.colour->first;
	int const anti_colour = part.colour->second;
	if(part.type == HEJ::ParticleID::gluon)
	return colour != anti_colour && colour > 0 && anti_colour > 0;
	if(HEJ::is_quark(part))
	return anti_colour == 0 && colour > 0;
	return colour == 0 && anti_colour > 0;
	}

	-bool correct_colour(HEJ::UnclusteredEvent const & ev){
	- for(auto const & part: ev.incoming){
	+bool correct_colour(HEJ::Event const & ev){
	+ for(auto const & part: ev.incoming()){
	if(!correct_colour(part))
	return false;
	}
	- for(auto const & part: ev.outgoing){
	+ for(auto const & part: ev.outgoing()){
	if(!correct_colour(part))
	return false;
	}
	return true;
	}

	bool match_expected(
	- HEJ::UnclusteredEvent const & ev,
	+ HEJ::Event const & ev,
	std::vector<HEJ::Colour> const & expected
	){
	- assert(ev.outgoing.size()+2==expected.size());
	- for(size_t i=0; i<ev.incoming.size(); ++i){
	- assert(ev.incoming[i].colour);
	- if( *ev.incoming[i].colour != expected[i])
	+ assert(ev.outgoing().size()+2==expected.size());
	+ for(size_t i=0; i<ev.incoming().size(); ++i){
	+ assert(ev.incoming()[i].colour);
	+ if( *ev.incoming()[i].colour != expected[i])
	return false;
	}
	- for(size_t i=2; i<ev.outgoing.size()+2; ++i){
	- if( ev.outgoing[i-2].colour ){
	- if( *ev.outgoing[i-2].colour != expected[i] )
	+ for(size_t i=2; i<ev.outgoing().size()+2; ++i){
	+ if( ev.outgoing()[i-2].colour ){
	+ if( *ev.outgoing()[i-2].colour != expected[i] )
	return false;
	} else if( expected[i].first != 0 \|\| expected[i].second != 0)
	return false;
	}
	return true;
	}

	void check_event(
	- HEJ::UnclusteredEvent ev, std::vector<HEJ::Colour> const & expected_colours
	+ HEJ::UnclusteredEvent unc_ev, std::vector<HEJ::Colour> const & expected_colours
	){
	+ HEJ::Event ev(
	+ unc_ev,
	+ fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30.);
	+ assert(HEJ::event_type::is_HEJ(ev.type()));
	dum_rnd rng;
	ev.generate_colour_flow(rng);
	if(!correct_colour(ev)){
	std::cerr << "Found illegal colours for event\n";
	dump_event(ev);
	throw std::invalid_argument("Illegal colour set");
	}
	if(!match_expected(ev, expected_colours)){
	std::cerr << "Colours didn't match expectation. Found\n";
	dump_event(ev);
	std::cerr << "but expected\n";
	for(auto const & col: expected_colours){
	std::cerr << "colour={" << col.first << ", " << col.second << "}\n";
	}
	throw std::logic_error("Colours did not match expectation");
	}
	}

	int main() {
	HEJ::UnclusteredEvent ev;
	std::vector<HEJ::Colour> expected_colours(7);

	/// pure gluon
	ev.incoming[0] = { HEJ::ParticleID::gluon, { 0, 0,-427, 427}, {}};
	ev.incoming[1] = { HEJ::ParticleID::gluon, { 0, 0, 851, 851}, {}};
	ev.outgoing.push_back({ HEJ::ParticleID::gluon, { 196, 124, -82, 246}, {}});
	ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-167,-184, 16, 249}, {}});
	ev.outgoing.push_back({ HEJ::ParticleID::higgs, { 197, 180, 168, 339}, {}});
	ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-190, -57, 126, 235}, {}});
	ev.outgoing.push_back({ HEJ::ParticleID::gluon, { -36, -63, 196, 209}, {}});

	expected_colours[0] = {502, 501};
	expected_colours[1] = {509, 502};
	expected_colours[2] = {503, 501};
	expected_colours[3] = {505, 503};
	expected_colours[4] = { 0, 0};
	expected_colours[5] = {507, 505};
	expected_colours[6] = {509, 507};
	check_event(ev, expected_colours);

	/// last g to Qx (=> gQx -> g ... Qx)
	ev.incoming[1].type = HEJ::ParticleID::d_bar;
	ev.outgoing[4].type = HEJ::ParticleID::d_bar;
	// => only end changes
	expected_colours[1].first = 0;
	expected_colours[6].first = 0;
	check_event(ev, expected_colours);

	{
	// don't overwrite
	auto new_expected = expected_colours;
	auto new_ev = ev;
	/// uno forward (=> gQx -> g ... Qx g)
	std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
	// => uno quarks eats colour and gluon connects to anti-colour
	new_expected[5] = {0, expected_colours[3].first};
	new_expected[6] = {expected_colours[0].first, expected_colours[0].first+2};
	new_expected[1].second += 2; // one more anti-colour in line
	check_event(new_ev, new_expected);
	}

	/// swap Qx <-> Q (=> gQx -> g ... Qx)
	ev.incoming[1].type = HEJ::ParticleID::d;
	ev.outgoing[4].type = HEJ::ParticleID::d;
	// => swap: colour<->anti && inital<->final
	std::swap(expected_colours[1], expected_colours[6]);
	std::swap(expected_colours[1].first, expected_colours[1].second);
	std::swap(expected_colours[6].first, expected_colours[6].second);
	check_event(ev, expected_colours);

	/// first g to q (=> qxQ -> qx ... Q)
	ev.incoming[0].type = HEJ::ParticleID::u_bar;
	ev.outgoing[0].type = HEJ::ParticleID::u_bar;
	expected_colours[0] = { 0, 501};
	// => shift anti-colour index one up
	expected_colours[1].first -= 2;
	expected_colours[5] = expected_colours[3];
	expected_colours[3] = expected_colours[2];
	expected_colours[2] = { 0, 502};
	check_event(ev, expected_colours);

	{
	// don't overwrite
	auto new_expected = expected_colours;
	auto new_ev = ev;
	/// uno backward (=> qxQ -> g qx ... Q)
	std::swap(new_ev.outgoing[0].type, new_ev.outgoing[1].type);
	// => uno gluon connects to quark colour
	new_expected[3] = expected_colours[2];
	new_expected[2] = {expected_colours[0].second+2, expected_colours[0].second};
	check_event(new_ev, new_expected);

	/// swap qx <-> q (=> qQ -> g q ... Q)
	new_ev.incoming[0].type = HEJ::ParticleID::u;
	new_ev.outgoing[1].type = HEJ::ParticleID::u;
	// => swap: colour<->anti && inital<->final
	std::swap(new_expected[0], new_expected[3]);
	std::swap(new_expected[0].first, new_expected[0].second);
	std::swap(new_expected[3].first, new_expected[3].second);
	// => & connect first gluon with remaining anti-colour
	new_expected[2] = {new_expected[0].first, new_expected[0].first+2};
	// shift colour line one down
	new_expected[1].first-=2;
	new_expected[5].first-=2;
	new_expected[5].second-=2;
	// shift anti-colour line one up
	new_expected[6].first+=2;
	check_event(new_ev, new_expected);
	}

	{
	// don't overwrite
	auto new_expected = expected_colours;
	auto new_ev = ev;
	/// uno forward (=> qxQ -> qx ... Q g)
	std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
	// => uno gluon connects to remaining colour
	new_expected[5] = expected_colours[6];
	new_expected[6] = {expected_colours[3].first+2, expected_colours[3].first};
	check_event(new_ev, new_expected);
	}

	/// @TODO add qqx test when implemented (it should work)

	return EXIT_SUCCESS;
	}

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions