No OneTemporary
Actions

Size

61 KB

Subscribers

None

View Options

	diff --git a/Changes-API.md b/Changes-API.md
	index 02cbab5..60baa9e 100644
	--- a/Changes-API.md
	+++ b/Changes-API.md
	@@ -1,98 +1,99 @@
	# Changelog for HEJ API

	This log lists only changes on the HEJ API. These are primarily code changes
	relevant for calling HEJ as an API. This file should only be read as an addition
	to [`Changes.md`](Changes.md), where the main features are documented.

	## Version 2.1

	### 2.1.0

	#### Changes to Event class
	* Restructured `Event` class
	- `Event` can now only be build from a (new) `Event::EventData` class
	- Removed default constructor for `Event`
	- `Event::EventData` replaces the old `UnclusteredEvent` struct.
	- `UnclusteredEvent` is now deprecated, and will be removed in HEJ Version
	2.3.0
	- Removed `Event.unclustered()` function
	- Added new member function `Events.parameters()`, to directly access
	(underlying) `Parameters<EventParameters>`
	- New member functions `begin_partons`, `end_partons` (`rbegin_partons`,
	`rend_partons`) with aliases `cbegin_partons`, `cend_partons`
	(`crbegin_partons`, `crend_partons`) for constant (reversed) iterators over
	outgoing partons.
	* New function `Event::EventData.reconstruct_intermediate()` to reconstruct
	bosons from decays, e.g. `positron + nu_e => Wp`
	* Added optional Colour charges to particles (`Particle.colour`)
	- Colour connection in the HEJ limit can be generated via
	`Event.generate_colours` (automatically done in the resummation)
	- Colour configuration of input events can be checked with
	`Event.is_leading_colour`
	* Added function `Event.valid_hej_state` to check if event _could have_ been
	produced by `HEJ` according to the `max ext soft pt fraction` cut on the jets

	#### Updated functions
	* Renamed `EventType::nonHEJ` to `EventType::non_resummable` and `is_HEJ()`
	to `is_resummable()` such that run card is consistent with internal workings
	* Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public
	* New `EventReweighter` member function `treatment` to query the
	treatment with respect to resummation for the various event types.
	* Added auxiliary functions to obtain a `std::string` from and
	`EventDescription` (`to_string` for human readable, and `to_simple_string`
	for easy parsable string)

	#### New classes
	* New class `Unweighter` to do unweighting
	* New class `CrossSectionAccumulator` to keep track of Cross Section of the
	different subprocess
	* New template struct `Parameters` similar to old `Weights`
	- `Weights` are now an alias for `Parameters<double>`. Calling `Weights` did
	not change
	- `Weights.hh` was replaced by `Parameters.hh`. The old `Weights.hh` header
	will be removed in HEJ Version 2.2.0
	* Function to multiplication and division of `EventParameters.weight` by double
	- This can be combined with `Parameters`, e.g.
	`Parameters<EventParameters>*Weights`, see also `Events.parameters()`
	- Moved `EventParameters` to `Parameters.hh` header
	* new class `EWConstants` replaces previously hard coded `vev`
	- `EWConstants` have to be set in the general `Config` and the
	`MatrixElementConfig`

	#### Input/Output
	* New abstract `EventReader` class, as base for reading events from files
	- Moved LHE file reader to `HEJ::LesHouchesReader`
	* Support reading (`HDF5Reader`) and writing (`HDF5Writer`) `hdf5` files
	* New `BufferedEventReader` class that allows to put events back into
	the reader.
	* New `SherpaLHEReader` to read Sherpa LHE files with correct weights
	* `get_analysis` now requires `YAML::Node` and `LHEF::HEPRUP` as arguments
	* Replaced `HepMCInterface` and `HepMCWriter` by `HepMCInterfaceX` and
	`HepMCWriterX` respectively, with `X` being the major version of HepMC (2 or
	3)
	- Renamed `HepMCInterfaceX::init_kinematics` to `HepMCInterfaceX::init_event`
	and protected it, use `HepMCInterfaceX::operator()` instead
	- Removed redundant `HepMCInterfaceX::add_variation` function

	#### Linking
	* Export cmake target `HEJ::HEJ` to link directly against `libHEJ`

	#### Miscellaneous
	* Capitalisation of `Config.hh` now matches class `Config` (was `config.hh`)
	* Replaced redundant member `EventReweighterConfig::jet_param` getter function
	`EventReweighter.jet_param()` (`JetParameters` are already in
	`PhaseSpacePointConfig`)
	-* The Random Number Generator is now stored as `std::shared_ptr<HEJ::RNG>`
	- inside classes (was `std::reference_wrapper`)
	- - Constructors of `EventReweighter`, `JetSplitter` and `PhaseSpacePoint` now
	- expect pointers (was reference)
	+* Avoid storing reference to the Random Number Generator inside classes
	+ - Constructors of `EventReweighter` now expect `std::shared_ptr<HEJ::RNG>`
	+ (was reference)
	+ - Moved reference to `HEJ::RNG` from constructor of `JetSplitter` to
	+ `JetSplitter.split`

	## Version 2.0

	### 2.0.4

	* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`

	### 2.0.0

	* First release
	diff --git a/include/HEJ/JetSplitter.hh b/include/HEJ/JetSplitter.hh
	index 38dbafc..bf0fabd 100644
	--- a/include/HEJ/JetSplitter.hh
	+++ b/include/HEJ/JetSplitter.hh
	@@ -1,72 +1,70 @@
	/**
	* \file
	* \brief Declaration of the JetSplitter class
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

	#include <memory>
	#include <vector>

	#include "fastjet/JetDefinition.hh"

	namespace fastjet {
	class PseudoJet;
	}

	namespace HEJ {
	class RNG;

	//! Class to split jets into their constituents
	class JetSplitter {
	public:
	struct SplitResult {
	std::vector<fastjet::PseudoJet> constituents;
	double weight;
	};

	//! Constructor
	/**
	- * @param jet_def Jet definition
	- * @param min_pt Minimum jet transverse momentum
	- * @param ran Random number generator
	+ * @param jet_def Jet definition
	+ * @param min_pt Minimum jet transverse momentum
	*/
	- JetSplitter(
	- fastjet::JetDefinition jet_def, double min_pt,
	- std::shared_ptr<RNG> ran
	- );
	+ JetSplitter(fastjet::JetDefinition jet_def, double min_pt);

	//! Split a get into constituents
	/**
	- * @param j2split Jet to be split
	- * @param ncons Number of constituents
	- * @returns The constituent momenta
	- * together with the associated weight
	+ * @param j2split Jet to be split
	+ * @param ncons Number of constituents
	+ * @param ran Random number generator
	+ * @returns The constituent momenta together with the associated
	+ * weight
	*/
	- SplitResult split(fastjet::PseudoJet const & j2split, int ncons) const;
	+ SplitResult split(
	+ fastjet::PseudoJet const & j2split, int ncons, RNG & ran
	+ ) const;

	//! Maximum distance of constituents to jet axis
	static constexpr double R_factor = 5./3.;

	private:
	//! \internal split jet into two partons
	- SplitResult Split2(fastjet::PseudoJet const & j2split) const;
	+ SplitResult Split2(fastjet::PseudoJet const & j2split, RNG & ran) const;

	/** \internal
	* @brief sample y-phi distance to jet pt axis for a jet splitting into two
	* partons
	*
	* @param wt Multiplied by the weight of the sampling point
	+ * @param ran Random number generator
	* @returns The distance in units of the jet radius
	*/
	- double sample_distance_2p(double & wt) const;
	+ double sample_distance_2p(double & wt, RNG & ran) const;

	- double R_;
	double min_jet_pt_;
	fastjet::JetDefinition jet_def_;
	- std::shared_ptr<RNG> ran_;
	};
	} // namespace HEJ
	diff --git a/include/HEJ/PhaseSpacePoint.hh b/include/HEJ/PhaseSpacePoint.hh
	index 9c36cc3..389f803 100644
	--- a/include/HEJ/PhaseSpacePoint.hh
	+++ b/include/HEJ/PhaseSpacePoint.hh
	@@ -1,197 +1,198 @@
	/** \file
	* \brief Contains the PhaseSpacePoint Class
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#pragma once

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

	#include "HEJ/Config.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/StatusCode.hh"

	namespace HEJ{
	class Event;
	class RNG;

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

	//! PhaseSpacePoint Constructor
	/**
	* @param ev Clustered Jet Event
	* @param conf Configuration parameters
	* @param ran Random number generator
	*/
	PhaseSpacePoint(
	Event const & ev,
	PhaseSpacePointConfig conf,
	- std::shared_ptr<RNG> ran
	+ RNG & ran
	);

	//! Get phase space point weight
	double weight() const{
	return weight_;
	}

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

	//! Access outgoing particles
	std::vector<Particle> const & outgoing() const{
	return outgoing_;
	}

	//! Particle decays
	/**
	* The key in the returned map corresponds to the index in the
	* vector returned by outgoing()
	*/
	std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
	return decays_;
	}

	//! Status code of generation
	StatusCode status() const{
	return status_;
	}

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

	private:
	//! /internal returns the clustered jets sorted in rapidity
	std::vector<fastjet::PseudoJet> cluster_jets(
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	bool pass_resummation_cuts(
	std::vector<fastjet::PseudoJet> const & jets
	) const;
	bool pass_extremal_cuts(
	fastjet::PseudoJet const & ext_parton,
	fastjet::PseudoJet const & jet
	) const;
	- int sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets);
	- int sample_ng_jets(int ng, std::vector<fastjet::PseudoJet> const & Born_jets);
	+ int sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran);
	+ int sample_ng_jets(
	+ int ng, std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran
	+ );
	double probability_in_jet(
	std::vector<fastjet::PseudoJet> const & Born_jets
	) const;
	std::vector<fastjet::PseudoJet> gen_non_jet(
	int ng_non_jet,
	- double ptmin, double ptmax
	+ double ptmin, double ptmax,
	+ RNG & ran
	);
	void rescale_qqx_rapidities(
	std::vector<fastjet::PseudoJet> & out_partons,
	std::vector<fastjet::PseudoJet> const & jets,
	const double ymin1, const double ymax2,
	const int qqxbackjet
	);
	void rescale_rapidities(
	std::vector<fastjet::PseudoJet> & partons,
	double ymin, double ymax
	);
	std::vector<fastjet::PseudoJet> reshuffle(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	fastjet::PseudoJet const & q
	);
	/** \interal final jet test
	* - number of jets must match Born kinematics
	* - no partons designated as nonjet may end up inside jets
	* - all other outgoing partons must end up inside jets
	* - the extremal (in rapidity) partons must be inside the extremal jets
	* - rapidities must be the same (by construction)
	*/
	bool jets_ok(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	void reconstruct_incoming(std::array<Particle, 2> const & Born_incoming);
	double phase_space_normalisation(
	int num_Born_jets,
	int num_res_partons
	) const;
	/** \interal Distribute gluon in jet
	* @param jets jets to distribute gluon in
	* @param ng_jets number of gluons
	* @param qqxbackjet position of first (backwards) qqx jet
	*
	* relies on JetSplitter
	*/
	std::vector<fastjet::PseudoJet> split(
	std::vector<fastjet::PseudoJet> const & jets,
	- int ng_jets, size_t qqxbackjet
	+ int ng_jets, size_t qqxbackjet, RNG & ran
	);
	std::vector<int> distribute_jet_partons(
	- int ng_jets, std::vector<fastjet::PseudoJet> const & jets
	+ int ng_jets, std::vector<fastjet::PseudoJet> const & jets, RNG & ran
	);
	std::vector<fastjet::PseudoJet> split(
	std::vector<fastjet::PseudoJet> const & jets,
	std::vector<int> const & np_in_jet,
	- size_t qqxbackjet
	+ size_t qqxbackjet,
	+ RNG & ran
	);
	bool split_preserved_jets(
	std::vector<fastjet::PseudoJet> const & jets,
	std::vector<fastjet::PseudoJet> const & jet_partons
	) const;
	template<class Particle>
	Particle const & most_backward_FKL(
	std::vector<Particle> const & partons
	) const;
	template<class Particle>
	Particle const & most_forward_FKL(
	std::vector<Particle> const & partons
	) const;
	template<class Particle>
	Particle & most_backward_FKL(std::vector<Particle> & partons) const;
	template<class Particle>
	Particle & most_forward_FKL(std::vector<Particle> & partons) const;
	bool extremal_ok(
	std::vector<fastjet::PseudoJet> const & partons
	) const;
	/** \internal
	* Assigns PDG IDs to outgoing partons, i.e. labels them as quarks
	*
	* \note This function assumes outgoing_ to be pure partonic when called,
	* i.e. A/W/Z/h bosons should _not be set_ at this stage
	*/
	void label_quarks(Event const & event);
	/** \internal
	* This function will label the qqx pair in a qqx event back to their
	* original types from the input event.
	*/
	void label_qqx(Event const & event);
	void copy_AWZH_boson_from(Event const & event);

	bool momentum_conserved() const;

	bool contains_idx(
	fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
	) const;

	bool unob_, unof_, qqxb_, qqxf_, qqxmid_;

	double weight_;

	PhaseSpacePointConfig param_;

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

	- std::shared_ptr<RNG> ran_;
	-
	StatusCode status_;
	};

	} // namespace HEJ
	diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
	index c27d033..bf28520 100644
	--- a/src/EventReweighter.cc
	+++ b/src/EventReweighter.cc
	@@ -1,284 +1,284 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/EventReweighter.hh"

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

	#include "fastjet/ClusterSequence.hh"

	#include "LHEF/LHEF.h"

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

	namespace HEJ{
	namespace {

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

	Event::EventData to_EventData(PhaseSpacePoint const & psp){
	Event::EventData result;
	result.incoming=psp.incoming();
	assert(result.incoming.size() == 2);
	result.outgoing=psp.outgoing();
	// technically Event::EventData doesn't have to be sorted,
	// but PhaseSpacePoint should be anyway
	assert(
	std::is_sorted(
	begin(result.outgoing), end(result.outgoing),
	rapidity_less{}
	)
	);
	assert(result.outgoing.size() >= 2);
	result.decays = psp.decays();
	result.parameters.central = {NaN, NaN, psp.weight()};
	return result;
	}

	} // namespace anonymous

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

	EventReweighter::EventReweighter(
	Beam beam,
	int pdf_id,
	ScaleGenerator scale_gen,
	EventReweighterConfig conf,
	std::shared_ptr<RNG> ran
	):
	param_{std::move(conf)},
	E_beam_{beam.E},
	pdf_{pdf_id, beam.type.front(), beam.type.back()},
	MEt2_{
	[this](double mu){ return pdf_.Halphas(mu); },
	param_.ME_config
	},
	scale_gen_{std::move(scale_gen)},
	ran_{std::move(ran)}
	{
	assert(ran_);
	}

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

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

	EventTreatment EventReweighter::treatment(EventType type) const {
	return param_.treat.at(type);
	}

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

	switch(treatment(ev.type())){
	case EventTreatment::discard: {
	status_.emplace_back(StatusCode::discard);
	return {};
	}
	case EventTreatment::keep:
	if(! jets_pass_resummation_cuts(ev)) {
	status_.emplace_back(StatusCode::failed_resummation_cuts);
	return {};
	}
	else {
	status_.emplace_back(StatusCode::good);
	return {ev};
	}
	default:;
	}
	const double Born_shat = shat(ev);

	std::vector<Event> resummation_events;
	status_.reserve(phase_space_points);
	for(size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){
	- PhaseSpacePoint psp{ev, param_.psp_config, ran_};
	+ PhaseSpacePoint psp{ev, param_.psp_config, *ran_};
	status_.emplace_back(psp.status());
	assert(psp.status() != StatusCode::unspecified);
	if(psp.status() != StatusCode::good) continue;
	assert(psp.weight() != 0.);
	if(psp.incoming()[0].E() > E_beam_ \|\| psp.incoming()[1].E() > E_beam_) {
	status_.back() = StatusCode::too_much_energy;
	continue;
	}

	resummation_events.emplace_back(
	to_EventData( std::move(psp) ).cluster(
	param_.jet_param().def, param_.jet_param().min_pt
	)
	);
	auto & new_event = resummation_events.back();
	assert( new_event.valid_hej_state(
	param_.psp_config.max_ext_soft_pt_fraction,
	param_.psp_config.min_extparton_pt ) );
	if( new_event.type() != ev.type() )
	throw std::logic_error{"Resummation Event does not match Born event"};
	new_event.generate_colours(*ran_);
	assert(new_event.variations().empty());
	new_event.central().mur = ev.central().mur;
	new_event.central().muf = ev.central().muf;
	const double resum_shat = shat(new_event);
	new_event.central().weight = ev.central().weightBorn_shat*Born_shat/
	(phase_space_pointsresum_shatresum_shat);
	}
	return resummation_events;
	}

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

	for(auto & cur_event: events){
	const auto pdf = pdf_factors(cur_event);
	assert(pdf.variations.size() == cur_event.variations().size());
	const auto ME = matrix_elements(cur_event);
	assert(ME.variations.size() == cur_event.variations().size());
	cur_event.parameters() = pdfME/(Born_pdf*Born_ME);
	}
	return events;
	}

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

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

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

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

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

	return MEt2_(ev);
	}

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

	Weights
	EventReweighter::fixed_order_scale_ME(Event const & ev) const{
	int alpha_s_power = 0;
	for(auto const & part: ev.outgoing()){
	if(is_parton(part))
	++alpha_s_power;
	else if(part.type == pid::Higgs) {
	alpha_s_power += 2;
	}
	// nothing to do for other uncoloured particles
	}

	Weights result;
	result.central = pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
	for(auto const & var: ev.variations()){
	result.variations.emplace_back(
	pow(pdf_.Halphas(var.mur), alpha_s_power)
	);
	}
	return result;
	}

	} // namespace HEJ
	diff --git a/src/JetSplitter.cc b/src/JetSplitter.cc
	index bf85e2f..96a9162 100644
	--- a/src/JetSplitter.cc
	+++ b/src/JetSplitter.cc
	@@ -1,191 +1,188 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/JetSplitter.hh"

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

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

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

	namespace HEJ {
	namespace{
	constexpr double ccut=HEJ::CMINPT; // min parton pt

	template<class Iterator>
	bool same_pt_and_rapidity(
	Iterator begin, Iterator end,
	fastjet::PseudoJet const & jet
	){
	constexpr double ep = 1e-2;
	const fastjet::PseudoJet reconstructed_jet = std::accumulate(
	begin, end, fastjet::PseudoJet{}
	);
	return
	(std::abs(reconstructed_jet.pt() - jet.pt()) < ep)
	&& (std::abs(reconstructed_jet.rapidity() - jet.rapidity()) < ep)
	;
	}

	bool all_in_one_jet(
	std::vector<fastjet::PseudoJet> const & partons,
	fastjet::JetDefinition jet_def, double min_jet_pt
	){
	fastjet::ClusterSequence ev(partons, jet_def);
	const std::vector<fastjet::PseudoJet> testjet = ev.inclusive_jets(min_jet_pt);
	return testjet.size() == 1u
	&& testjet[0].constituents().size() == partons.size();
	}
	}
	JetSplitter::JetSplitter(
	- fastjet::JetDefinition jet_def, double min_pt,
	- std::shared_ptr<RNG> ran
	+ fastjet::JetDefinition jet_def, double min_pt
	):
	- R_{jet_def.R()},
	min_jet_pt_{min_pt},
	- jet_def_{jet_def},
	- ran_{std::move(ran)}
	- {
	- assert(ran_);
	- }
	+ jet_def_{std::move(jet_def)}
	+ {}

	using SplitResult = JetSplitter::SplitResult;

	SplitResult JetSplitter::split(
	- fastjet::PseudoJet const & j2split, int ncons
	+ fastjet::PseudoJet const & j2split, int ncons, RNG & ran
	) const{
	if(ncons <= 0) {
	throw std::invalid_argument{
	"number of requested jet constituents less than 1"
	};
	}
	double swt = 1.;

	std::vector<fastjet::PseudoJet> jcons;
	if(ncons == 1){
	jcons.emplace_back(j2split);
	jcons.back().set_user_index(0);
	return {jcons, swt};
	}
	if(ncons == 2){
	- return Split2(j2split);
	+ return Split2(j2split, ran);
	}
	- const double R_max = R_factor*R_;
	+ const double R_max = R_factor*jet_def_.R();
	assert(R_max < M_PI);

	double pt_remaining = j2split.pt();
	const double phi_jet = j2split.phi();
	const double y_jet = j2split.rapidity();
	for(int i = 0; i < ncons - 1; ++i){
	/**
	* Generate rapidity and azimuthal angle with a distance
	* R = sqrt(delta_y^2 + delta_phi^2) < R_max
	* from the jet centre
	*/
	- const double R = R_max*ran_->flat();
	- const double theta = 2M_PIran_->flat();
	+ const double R = R_max*ran.flat();
	+ const double theta = 2M_PIran.flat();
	const double delta_phi = R*cos(theta);
	const double delta_y = R*sin(theta);

	/**
	* Generate pt such that the total contribution of all partons
	* along the jet pt axis does not exceed the jet pt
	*/
	const double pt_max = pt_remaining/cos(delta_phi);
	assert(pt_max > 0);
	if(pt_max < ccut) return {}; // no pt remaining for this parton
	- const double pt = (pt_max - ccut)*ran_->flat() + ccut;
	+ const double pt = (pt_max - ccut)*ran.flat() + ccut;
	pt_remaining -= pt*cos(delta_phi);

	jcons.emplace_back(
	ptcos(phi_jet + delta_phi), ptsin(phi_jet + delta_phi),
	ptsinh(y_jet + delta_y), ptcosh(y_jet + delta_y)
	);
	jcons.back().set_user_index(i);
	swt = 2M_PIRR_maxpt(pt_max - ccut);
	}

	const fastjet::PseudoJet p_total = std::accumulate(
	jcons.begin(), jcons.end(), fastjet::PseudoJet{}
	);

	// Calculate the pt of the last parton
	const double last_px = j2split.px() - p_total.px();
	const double last_py = j2split.py() - p_total.py();
	const double last_pt = sqrt(last_pxlast_px + last_pylast_py);
	if(last_pt < ccut) return {};

	// Calculate the rapidity of the last parton using the requirement that the
	// new jet must have the same rapidity as the LO jet.
	const double exp_2y_jet = (j2split.e() + j2split.pz())/(j2split.e() - j2split.pz());
	const double bb = (p_total.e()+p_total.pz()) - exp_2y_jet*(p_total.e()-p_total.pz());
	const double lasty = log((-bb+sqrt(bbbb+4.exp_2y_jetlast_ptlast_pt))/(2.*last_pt));

	jcons.emplace_back(
	last_px, last_py, last_ptsinh(lasty), last_ptcosh(lasty)
	);
	jcons.back().set_user_index(ncons-1);
	assert(same_pt_and_rapidity(begin(jcons), end(jcons), j2split));

	// Test that the last parton is not too far away from the jet centre.
	if (jcons.back().delta_R(j2split) > R_max) return {};

	if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};

	return {jcons, swt};
	}

	- double JetSplitter::sample_distance_2p(double & wt) const{
	+ double JetSplitter::sample_distance_2p(double & wt, RNG & ran) const{
	static constexpr double x_small = 0.1;
	static constexpr double p_small = 0.4;

	- const double pR = ran_->flat();
	+ const double pR = ran.flat();
	if(pR < p_small){
	wt *= x_small/p_small;
	return x_small/p_small*pR;
	}
	wt *= (1-x_small)/(1-p_small);
	return (1-x_small)/(1-p_small)*(pR-p_small) + x_small;
	}

	- SplitResult JetSplitter::Split2(fastjet::PseudoJet const & j2split) const{
	+ SplitResult JetSplitter::Split2(
	+ fastjet::PseudoJet const & j2split, RNG & ran
	+ ) const{
	static constexpr size_t ncons = 2;
	std::vector<fastjet::PseudoJet> jcons(ncons);
	std::array<double, ncons> R, phi, y, pt;
	double wt = 1;

	- const double theta = 2M_PIran_->flat(); // angle in y-phi plane
	+ const double theta = 2M_PIran.flat(); // angle in y-phi plane
	// empiric observation: we are always within the jet radius
	- R[0] = sample_distance_2p(wt)*R_;
	- R[1] = -sample_distance_2p(wt)*R_;
	+ R[0] = sample_distance_2p(wt, ran)*jet_def_.R();
	+ R[1] = -sample_distance_2p(wt, ran)*jet_def_.R();
	for(size_t i = 0; i <= 1; ++i){
	phi[i] = j2split.phi() + R[i]*cos(theta);
	y[i] = j2split.rapidity() + R[i]*sin(theta);
	}
	for(size_t i = 0; i <= 1; ++i){
	pt[i] = (j2split.py() - tan(phi[1-i])*j2split.px())/
	(sin(phi[i]) - tan(phi[1-i])*cos(phi[i]));
	if(pt[i] < ccut) return {};
	jcons[i].reset_PtYPhiM(pt[i], y[i], phi[i]);
	jcons[i].set_user_index(i);
	}
	assert(same_pt_and_rapidity(begin(jcons), end(jcons), j2split));

	if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
	- wt = 2M_PIpt[0]R[0]R_R_;
	+ wt = 2M_PIpt[0]R[0]jet_def_.R()jet_def_.R();
	// from transformation of delta(R[1] - ...) to delta(pt[0] - ...)
	const double dphi0 = phi[0] - j2split.phi();
	const double ptJ = j2split.pt();
	const double jacobian = cos(theta)pt[1]pt[1]/(ptJ*sin(dphi0));
	return {jcons, jacobian*wt};
	}
	}
	diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
	index aef2b4e..9c5b669 100644
	--- a/src/PhaseSpacePoint.cc
	+++ b/src/PhaseSpacePoint.cc
	@@ -1,811 +1,814 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/PhaseSpacePoint.hh"

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

	#include "fastjet/ClusterSequence.hh"

	#include "HEJ/Constants.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/JetSplitter.hh"
	#include "HEJ/kinematics.hh"
	#include "HEJ/resummation_jet.hh"
	#include "HEJ/utility.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/event_types.hh"

	namespace HEJ{
	namespace {
	constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;

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

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

	// user indices for partons with extremal rapidity
	constexpr int qqxmid1_idx = -9;
	constexpr int qqxmid2_idx = -8;
	constexpr int qqxb_idx = -7;
	constexpr int qqxf_idx = -6;
	constexpr int unob_idx = -5;
	constexpr int unof_idx = -4;
	constexpr int backward_FKL_idx = -3;
	constexpr int forward_FKL_idx = -2;
	}

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

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

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

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

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

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

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

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

	namespace {

	auto get_first_anyquark_emission(Event const & ev) {
	// find born quarks (ignore extremal partons)
	auto const firstquark = std::find_if(
	std::next(ev.begin_partons()), std::prev(ev.end_partons(), 2),
	[](Particle const & s){ return (is_anyquark(s)); }
	);
	// assert that it is a q-q_bar pair.
	assert(std::distance(firstquark, ev.end_partons()) != 2);
	assert(
	( is_quark(firstquark) && is_antiquark(std::next(firstquark)) )
	\|\| ( is_antiquark(firstquark) && is_quark(std::next(firstquark)) )
	);
	return firstquark;
	}

	//! returns index of most backward q-qbar jet
	template<class Iterator>
	int get_back_quark_jet(Event const & ev, Iterator firstquark){
	// find jets at FO corresponding to the quarks
	// technically this isn't necessary for LO
	std::vector<int> const born_indices{ ev.particle_jet_indices() };
	const auto firstquark_idx = std::distance(ev.begin_partons(), firstquark);
	int const firstjet_idx = born_indices[firstquark_idx];
	assert(firstjet_idx>0);
	assert( born_indices[firstquark_idx+1] == firstjet_idx+1 );

	return firstjet_idx;
	}

	//! returns index of most backward q-qbar jet
	int getBackQuarkJet(Event const & ev){
	const auto firstquark = get_first_anyquark_emission(ev);
	return get_back_quark_jet(ev, firstquark);
	}

	template<class ConstIterator, class Iterator>
	void label_extremal_qqx(
	ConstIterator born_begin, ConstIterator born_end,
	Iterator first_out
	){
	// find born quarks
	const auto firstquark = std::find_if(
	born_begin, born_end-1,
	[](Particle const & s){ return (is_anyquark(s)); }
	);
	assert(firstquark != born_end-1);
	const auto secondquark = std::find_if(
	firstquark+1, born_end,
	[](Particle const & s){ return (is_anyquark(s)); }
	);
	assert(secondquark != born_end);
	assert( ( is_quark(firstquark) && is_antiquark(secondquark) )
	\|\| ( is_antiquark(firstquark) && is_quark(secondquark) ));
	assert(first_out->type == ParticleID::gluon);
	assert((first_out+1)->type == ParticleID::gluon);

	// copy type from born
	first_out->type = firstquark->type;
	(first_out+1)->type = secondquark->type;
	}
	}

	void PhaseSpacePoint::label_qqx(Event const & event){
	assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
	assert(filter_partons(outgoing_).size() == outgoing_.size());
	if(qqxb_){
	label_extremal_qqx( event.outgoing().cbegin(), event.outgoing().cend(),
	outgoing_.begin()
	);
	return;
	}
	if(qqxf_){ // same as qqxb with reversed order
	label_extremal_qqx( event.outgoing().crbegin(), event.outgoing().crend(),
	outgoing_.rbegin()
	);
	return;
	}
	// central qqx
	const auto firstquark = get_first_anyquark_emission(event);

	// find jets at FO corresponding to the quarks
	// technically this isn't necessary for LO
	const auto firstjet_idx = get_back_quark_jet(event, firstquark);

	// find corresponding jets after resummation
	fastjet::ClusterSequence cs{to_PseudoJet(outgoing_), param_.jet_param.def};
	auto const jets = fastjet::sorted_by_rapidity(
	cs.inclusive_jets( param_.jet_param.min_pt ));
	std::vector<int> const resum_indices{ cs.particle_jet_indices({jets}) };

	// assert that jets didn't move
	assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx )->rapidity(),
	jets[ firstjet_idx ].rapidity(), 1e-2) );
	assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx+1 )->rapidity(),
	jets[ firstjet_idx+1 ].rapidity(), 1e-2) );

	// find last partons in first (central) jet
	size_t idx_out = 0;
	for(size_t i=resum_indices.size()-2; i>0; --i)
	if(resum_indices[i] == firstjet_idx){
	idx_out = i;
	break;
	}
	assert(idx_out != 0);

	// check that there is sufficient pt in jets from the quarks
	const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction;
	if (outgoing_[idx_out].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx )->pt()){
	weight_=0.;
	status_ = StatusCode::wrong_jets;
	return;
	}

	if (outgoing_[idx_out+1].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx+1 )->pt()){
	weight_=0.;
	status_ = StatusCode::wrong_jets;
	return;
	}
	// check that no additional emission between jets
	// such configurations are possible if we have an gluon gets generated
	// inside the rapidities of the qqx chain, but clusted to a
	// differnet/outside jet. Changing this is non trivial
	if(resum_indices[idx_out+1] != resum_indices[idx_out]+1){
	weight_=0.;
	status_ = StatusCode::gluon_in_qqx;
	return;
	}
	outgoing_[idx_out].type = firstquark->type;
	outgoing_[idx_out+1].type = std::next(firstquark)->type;
	}

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

	if (WEmit != end(ev.outgoing())){
	if(!qqxb_) {
	const size_t backward_FKL_idx = unob_?1:0;
	const auto backward_FKL = std::next(ev.begin_partons(), backward_FKL_idx);
	outgoing_[backward_FKL_idx].type = backward_FKL->type;
	}
	if(!qqxf_) {
	const size_t forward_FKL_idx = unof_?1:0;
	const auto forward_FKL = std::prev(ev.end_partons(), 1+forward_FKL_idx);
	outgoing_.rbegin()[unof_].type = forward_FKL->type;
	}
	} else {
	most_backward_FKL(outgoing_).type = ev.incoming().front().type;
	most_forward_FKL(outgoing_).type = ev.incoming().back().type;
	}

	if(qqxmid_\|\|qqxb_\|\|qqxf_){
	label_qqx(ev);
	}
	}

	PhaseSpacePoint::PhaseSpacePoint(
	- Event const & ev, PhaseSpacePointConfig conf, std::shared_ptr<RNG> ran
	+ Event const & ev, PhaseSpacePointConfig conf, RNG & ran
	):
	unob_{ev.type() == event_type::unob},
	unof_{ev.type() == event_type::unof},
	qqxb_{ev.type() == event_type::qqxexb},
	qqxf_{ev.type() == event_type::qqxexf},
	qqxmid_{ev.type() == event_type::qqxmid},
	param_{std::move(conf)},
	- ran_{std::move(ran)},
	status_{unspecified}
	{
	- assert(ran_);
	weight_ = 1;
	const auto & Born_jets = ev.jets();
	- const int ng = sample_ng(Born_jets);
	+ const int ng = sample_ng(Born_jets, ran);
	weight_ /= std::tgamma(ng + 1);
	- const int ng_jets = sample_ng_jets(ng, Born_jets);
	+ const int ng_jets = sample_ng_jets(ng, Born_jets, ran);
	std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
	- ng - ng_jets, CMINPT, param_.jet_param.min_pt
	+ ng - ng_jets, CMINPT, param_.jet_param.min_pt, ran
	);

	int qqxbackjet(-1);
	if(qqxmid_){
	qqxbackjet = getBackQuarkJet(ev);
	}

	const auto qperp = std::accumulate(
	begin(out_partons), end(out_partons),
	fastjet::PseudoJet{}
	);
	const auto jets = reshuffle(Born_jets, qperp);
	if(weight_ == 0.) {
	status_ = failed_reshuffle;
	return;
	}
	if(! pass_resummation_cuts(jets)){
	status_ = failed_resummation_cuts;
	weight_ = 0.;
	return;
	}
	- std::vector<fastjet::PseudoJet> jet_partons = split(jets, ng_jets, qqxbackjet);
	+ std::vector<fastjet::PseudoJet> jet_partons = split(
	+ jets, ng_jets, qqxbackjet, ran
	+ );
	if(weight_ == 0.) {
	status_ = StatusCode::failed_split;
	return;
	}
	if(qqxmid_){
	rescale_qqx_rapidities(
	out_partons, jets,
	most_backward_FKL(jet_partons).rapidity(),
	most_forward_FKL(jet_partons).rapidity(),
	qqxbackjet
	);
	}
	else{
	rescale_rapidities(
	out_partons,
	most_backward_FKL(jet_partons).rapidity(),
	most_forward_FKL(jet_partons).rapidity()
	- );
	+ );
	}
	if(! cluster_jets(out_partons).empty()){
	weight_ = 0.;
	status_ = StatusCode::empty_jets;
	return;
	}
	std::sort(begin(out_partons), end(out_partons), rapidity_less{});
	assert(
	std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{})
	);
	const auto first_jet_parton = out_partons.insert(
	end(out_partons), begin(jet_partons), end(jet_partons)
	);
	std::inplace_merge(
	begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{}
	);

	if(! jets_ok(Born_jets, out_partons)){
	weight_ = 0.;
	status_ = StatusCode::wrong_jets;
	return;
	}

	weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size());

	outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H
	for(auto & p: out_partons){
	outgoing_.emplace_back(Particle{pid::gluon, std::move(p), {}});
	}
	assert(!outgoing_.empty());

	label_quarks(ev);
	if(weight_ == 0.) {
	//! @TODO optimise s.t. this is not possible
	// status is handled internally
	return;
	}

	copy_AWZH_boson_from(ev);

	reconstruct_incoming(ev.incoming());
	status_ = StatusCode::good;
	}

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

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

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

	void PhaseSpacePoint::rescale_qqx_rapidities(
	std::vector<fastjet::PseudoJet> & out_partons,
	std::vector<fastjet::PseudoJet> const & jets,
	const double ymin1, const double ymax2,
	const int qqxbackjet
	){
	const double ymax1 = jets[qqxbackjet].rapidity();
	const double ymin2 = jets[qqxbackjet+1].rapidity();
	constexpr double ep = 1e-7;
	const double tot_y = ymax1 - ymin1 + ymax2 - ymin2;

	std::vector<std::reference_wrapper<fastjet::PseudoJet>> refpart(
	out_partons.begin(), out_partons.end());

	double ratio = (ymax1 - ymin1)/tot_y;

	const auto gap{ std::find_if(refpart.begin(), refpart.end(),
	[ratio](fastjet::PseudoJet p){
	return (p.rapidity()>=ratio);} ) };

	double ymin = ymin1;
	double ymax = ymax1;
	double dy = ymax - ymin - 2*ep;
	double offset = 0.;
	for(auto it_part=refpart.begin(); it_part<refpart.end(); ++it_part){
	if(it_part == gap){
	ymin = ymin2;
	ymax = ymax2;
	dy = ymax - ymin - 2*ep;
	offset = ratio;
	ratio = 1-ratio;
	}
	fastjet::PseudoJet & part = *it_part;
	assert(offset <= part.rapidity() && part.rapidity() < ratio+offset);
	const double y = ymin + ep + dy*((part.rapidity()-offset)/ratio);
	part.reset_momentum_PtYPhiM(part.pt(), y, part.phi());
	weight_ = tot_y-4.ep;
	assert(ymin <= part.rapidity() && part.rapidity() <= ymax);
	}
	assert(is_sorted(begin(out_partons), end(out_partons), rapidity_less{}));
	}

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

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

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

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

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

	std::vector<fastjet::PseudoJet> PhaseSpacePoint::reshuffle(
	std::vector<fastjet::PseudoJet> const & Born_jets,
	fastjet::PseudoJet const & q
	){
	if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
	const auto jets = resummation_jet_momenta(Born_jets, q);
	if(jets.empty()){
	weight_ = 0;
	return {};
	}

	// additional Jacobian to ensure Born integration over delta gives 1
	weight_ *= resummation_jet_weight(Born_jets, q);
	return jets;
	}

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

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

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

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

	} // namespace anonymous
	#endif

	std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
	std::vector<fastjet::PseudoJet> const & jets,
	- int ng_jets,
	- size_t qqxbackjet
	+ int ng_jets, size_t qqxbackjet, RNG & ran
	){
	- return split(jets, distribute_jet_partons(ng_jets, jets), qqxbackjet);
	+ return split(
	+ jets, distribute_jet_partons(ng_jets, jets, ran), qqxbackjet, ran);
	}

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

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

	const size_t most_backward_FKL_idx = 0 + unob_ + qqxb_;
	const size_t most_forward_FKL_idx = jets.size() - 1 - unof_ - qqxf_;
	const auto & jet = param_.jet_param;
	- const JetSplitter jet_splitter{jet.def, jet.min_pt, ran_};
	+ const JetSplitter jet_splitter{jet.def, jet.min_pt};

	std::vector<fastjet::PseudoJet> jet_partons;
	// randomly distribute jet gluons among jets
	for(size_t i = 0; i < jets.size(); ++i){
	- auto split_res = jet_splitter.split(jets[i], np[i]);
	+ auto split_res = jet_splitter.split(jets[i], np[i], ran);
	weight_ *= split_res.weight;
	if(weight_ == 0) return {};
	assert(
	std::all_of(
	begin(split_res.constituents), end(split_res.constituents),
	is_jet_parton
	)
	);
	const auto first_new_parton = jet_partons.insert(
	end(jet_partons),
	begin(split_res.constituents), end(split_res.constituents)
	);
	// mark uno and extremal FKL emissions here so we can check
	// their position once all emissions are generated
	// also mark qqxmid partons, and apply appropriate pt cut.
	auto extremal = end(jet_partons);
	if (i == most_backward_FKL_idx){ //FKL backward emission
	extremal = std::min_element(
	first_new_parton, end(jet_partons), rapidity_less{}
	);
	extremal->set_user_index(backward_FKL_idx);
	}
	else if(((unob_ \|\| qqxb_) && i == 0)){
	// unordered/qqxb
	extremal = std::min_element(
	first_new_parton, end(jet_partons), rapidity_less{}
	);
	extremal->set_user_index((unob_)?unob_idx:qqxb_idx);
	}

	else if (i == most_forward_FKL_idx){
	extremal = std::max_element(
	first_new_parton, end(jet_partons), rapidity_less{}
	);
	extremal->set_user_index(forward_FKL_idx);
	}
	else if(((unof_ \|\| qqxf_) && i == jets.size() - 1)){
	// unordered/qqxf
	extremal = std::max_element(
	first_new_parton, end(jet_partons), rapidity_less{}
	);
	extremal->set_user_index((unof_)?unof_idx:qqxf_idx);
	}
	else if((qqxmid_ && i == qqxbackjet)){
	extremal = std::max_element(
	first_new_parton, end(jet_partons), rapidity_less{}
	);
	extremal->set_user_index(qqxmid1_idx);
	}
	else if((qqxmid_ && i == qqxbackjet+1)){
	extremal = std::min_element(
	first_new_parton, end(jet_partons), rapidity_less{}
	);
	extremal->set_user_index(qqxmid2_idx);
	}
	if(
	extremal != end(jet_partons)
	&& !pass_extremal_cuts(*extremal, jets[i])
	){
	weight_ = 0;
	return {};
	}
	}
	assert(tagged_FKL_extremal(jet_partons));
	std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
	if(
	!extremal_ok(jet_partons)
	\|\| !split_preserved_jets(jets, jet_partons)
	){
	weight_ = 0.;
	return {};
	}
	return jet_partons;
	}

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

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

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

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

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

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

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

	bool PhaseSpacePoint::contains_idx(
	fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
	) const {
	auto const & constituents = jet.constituents();
	const int idx = parton.user_index();
	const bool injet = std::find_if(
	begin(constituents), end(constituents),
	[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
	) != end(constituents);
	const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction;
	return ((parton.pt()>minpartonjetpt*jet.pt())&&injet);
	}

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

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

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

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

	} //namespace HEJ
	diff --git a/t/test_psp.cc b/t/test_psp.cc
	index b3283b2..c62e0c4 100644
	--- a/t/test_psp.cc
	+++ b/t/test_psp.cc
	@@ -1,68 +1,68 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "LHEF/LHEF.h"
	#include "HEJ/stream.hh"
	#include "HEJ/Config.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/PhaseSpacePoint.hh"
	#include "HEJ/Ranlux64.hh"

	#include "hej_test.hh"

	namespace{
	constexpr int max_trials = 100;
	constexpr double max_ext_soft_pt_fraction = 0.1;
	const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
	constexpr double min_jet_pt = 50;
	}

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

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

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

	- std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Ranlux64>()};
	+ HEJ::Ranlux64 ran{};

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

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:38 PM (1 d, 2 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806095
Default Alt Text: (61 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions