No OneTemporary
Actions

Size

97 KB

Subscribers

None

View Options

	diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
	index 2f3d62c..87b7185 100644
	--- a/.gitlab-ci.yml
	+++ b/.gitlab-ci.yml
	@@ -1,156 +1,296 @@
	-.tags_template:
	- tags: &tags_def
	- - docker
	+# ---------------------------------
	+# - General Setup -
	+# ---------------------------------

	stages:
	- - minimal_compile
	- - different_compile
	+ - build
	+ - test
	+ - FOG:build
	+ - FOG:test
	- clean_code
	- publish

	variables:
	- HEJ_Install_DIR: tmp_HEJ/HEJ_installed
	- HEJ_build_DIR: tmp_HEJ/build_HEJ
	- FOG_build_DIR: tmp_HEJ/build_FOG
	+ # build directories
	+ HEJ_BUILD_DIR: tmp_HEJ/HEJ_build
	+ HEJ_INSTALL_DIR: tmp_HEJ/HEJ_installed
	+ FOG_BUILD_DIR: tmp_HEJ/FOG_build
	+ FOG_INSTALL_DIR: ${HEJ_INSTALL_DIR}
	+ # docker images
	+ DOCKER_BASIC: hejdock/hepenv
	+ DOCKER_HEPMC3: hejdock/hepmc3env
	+ DOCKER_QCDLOOP: hejdock/qcdloopenv
	+ DOCKER_RIVET: hejdock/rivetenv
	+
	+# ----------- Macros -----------

	after_script:
	- date

	-.HEJ_compile:
	+.tags_template:
	+ tags: &tags_def
	+ - docker
	+
	+# save complete history of failed tests
	+.save_failure:
	+ artifacts: &artifacts_failed
	+ when: on_failure
	+ untracked: true
	+ expire_in: 3d
	+
	+# ---------------------------------
	+# - Script Templates -
	+# ---------------------------------
	+
	+# ----------- Build -----------
	+
	+.HEJ_build:
	tags: *tags_def
	+ stage: build
	before_script:
	- date
	- source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh \|\| exit 1
	- - export t_HEJ_Install_DIR=${PWD}/$HEJ_Install_DIR
	- - export t_HEJ_build_DIR=${PWD}/$HEJ_build_DIR
	- - export t_HEJ_DIR=${PWD}
	- - mkdir -p ${t_HEJ_build_DIR}
	- - cd ${t_HEJ_build_DIR}
	- - cmake ${t_HEJ_DIR} -DCMAKE_BUILD_TYPE=DEBUG -DCMAKE_INSTALL_PREFIX=${t_HEJ_Install_DIR}
	- - make -j 8
	-
	-Minimal_setup:
	- stage: minimal_compile
	- extends: .HEJ_compile
	- image: hejdock/hepenv
	+ # prepare build
	+ - t_HEJ_DIR=${PWD}
	+ - t_HEJ_INSTALL_DIR=${PWD}/${HEJ_INSTALL_DIR}
	+ - t_HEJ_BUILD_DIR=${PWD}/${HEJ_BUILD_DIR}
	+ - mkdir -p ${t_HEJ_BUILD_DIR}
	+ - cd ${t_HEJ_BUILD_DIR}
	+ - cmake ${t_HEJ_DIR} -DCMAKE_BUILD_TYPE=DEBUG
	+ -DCMAKE_INSTALL_PREFIX=${t_HEJ_INSTALL_DIR}
	script:
	- - make test
	+ - make -j $(nproc --ignore=1)
	+ - make install
	+ artifacts:
	+ # save build and installed folder
	+ name: build
	+ expire_in: 1d
	+ paths:
	+ - ${HEJ_INSTALL_DIR}
	+ - ${HEJ_BUILD_DIR}

	-FOG:
	- stage: minimal_compile
	- extends: .HEJ_compile
	- image: hejdock/hepenv
	+# ----------- Test -----------
	+
	+.HEJ_test:
	+ tags: *tags_def
	+ stage: test
	+ before_script:
	+ - date
	+ - source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh \|\| exit 1
	+ # load HEJ
	+ - t_HEJ_DIR=${PWD}
	+ - t_HEJ_INSTALL_DIR=${PWD}/${HEJ_INSTALL_DIR}
	+ - export LD_LIBRARY_PATH=${t_HEJ_INSTALL_DIR}/lib:${LD_LIBRARY_PATH}
	+ - export PATH=${t_HEJ_INSTALL_DIR}/bin:${PATH}
	+ - cd ${HEJ_BUILD_DIR}
	+ - cmake ${t_HEJ_DIR} # rerun cmake to create all test configure files
	script:
	+ - ctest --output-on-failure
	+ artifacts: *artifacts_failed
	+
	+## ----------- FOG build -----------
	+
	+.FOG_build:
	+ tags: *tags_def
	+ stage: FOG:build
	+ before_script:
	- date
	- # install HEJ
	- - make install
	- # setup env
	- - export LD_LIBRARY_PATH=${t_HEJ_Install_DIR}/lib:${LD_LIBRARY_PATH}
	- - export PATH=${t_HEJ_Install_DIR}/bin:${PATH}
	- - export t_FOG_build_DIR=${PWD}/$FOG_build_DIR
	- - export t_FOG_DIR=${t_HEJ_DIR}/FixedOrderGen
	- # compile
	- - mkdir -p ${t_FOG_build_DIR}
	- - cd ${t_FOG_build_DIR}
	+ - source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh \|\| exit 1
	+ # load HEJ
	+ - t_HEJ_INSTALL_DIR=${PWD}/${HEJ_INSTALL_DIR}
	+ - export LD_LIBRARY_PATH=${t_HEJ_INSTALL_DIR}/lib:${LD_LIBRARY_PATH}
	+ - export PATH=${t_HEJ_INSTALL_DIR}/bin:${PATH}
	+ # prepare build
	+ - t_FOG_DIR=${PWD}/FixedOrderGen
	+ - t_FOG_INSTALL_DIR=${PWD}/${FOG_INSTALL_DIR}
	+ - t_FOG_BUILD_DIR=${PWD}/${FOG_BUILD_DIR}
	+ - mkdir -p ${t_FOG_BUILD_DIR}
	+ - cd ${t_FOG_BUILD_DIR}
	- cmake ${t_FOG_DIR} -DCMAKE_BUILD_TYPE=DEBUG
	- - make -j 8
	- - make test
	+ -DCMAKE_INSTALL_PREFIX=${t_FOG_INSTALL_DIR}
	+ script:
	+ - make -j $(nproc --ignore=1)
	+ - make install
	+ artifacts:
	+ # save build and installed folder
	+ name: FOG_build
	+ expire_in: 1d
	+ paths:
	+ - ${HEJ_INSTALL_DIR}
	+ - ${FOG_INSTALL_DIR}
	+ - ${FOG_BUILD_DIR}

	+## ----------- FOG test -----------

	-rivet:
	- stage: different_compile
	- extends: .HEJ_compile
	- image: hejdock/rivetenv
	- variables:
	- HEJ_build_DIR: tmp_HEJ/build_HEJ_rivet
	+.FOG_test:
	+ tags: *tags_def
	+ stage: FOG:test
	+ before_script:
	+ - date
	+ - source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh \|\| exit 1
	+ # load HEJ
	+ - t_FOG_DIR=${PWD}/FixedOrderGen
	+ - t_HEJ_INSTALL_DIR=${PWD}/${HEJ_INSTALL_DIR}
	+ - t_FOG_INSTALL_DIR=${PWD}/${FOG_INSTALL_DIR}
	+ - export LD_LIBRARY_PATH=${t_HEJ_INSTALL_DIR}/lib:${LD_LIBRARY_PATH}
	+ - export PATH=${t_HEJ_INSTALL_DIR}/bin:${t_FOG_INSTALL_DIR}/bin:${PATH}
	+ - t_FOG_BUILD_DIR=${PWD}/${FOG_BUILD_DIR}
	+ - cd ${t_FOG_BUILD_DIR}
	+ - cmake ${t_FOG_DIR} # rerun cmake to create all test configure files
	script:
	- make test
	+ artifacts: *artifacts_failed
	+
	+# ---------------------------------
	+# - Build & Test -
	+# ---------------------------------
	+
	+# ----------- basic -----------
	+
	+build:basic:
	+ image: ${DOCKER_BASIC}
	+ extends: .HEJ_build
	+
	+test:basic:
	+ image: ${DOCKER_BASIC}
	+ extends: .HEJ_test
	+ dependencies:
	+ - build:basic
	+
	+FOG:build:basic:
	+ image: ${DOCKER_BASIC}
	+ extends: .FOG_build
	+ dependencies:
	+ - build:basic
	+
	+FOG:test:basic:
	+ image: ${DOCKER_BASIC}
	+ extends: .FOG_test
	+ dependencies:
	+ - FOG:build:basic
	+
	+# ----------- HepMC 3 -----------
	+
	+build:hepmc3:
	+ image: ${DOCKER_HEPMC3}
	+ extends: .HEJ_build
	+
	+test:hepmc3:
	+ image: ${DOCKER_HEPMC3}
	+ extends: .HEJ_test
	+ dependencies:
	+ - build:hepmc3
	+
	+
	+# ----------- QCDloop -----------
	+
	+build:qcdloop:
	+ image: ${DOCKER_QCDLOOP}
	+ extends: .HEJ_build
	+
	+test:qcdloop:
	+ image: ${DOCKER_QCDLOOP}
	+ extends: .HEJ_test
	+ dependencies:
	+ - build:qcdloop
	+
	+# ----------- rivet -----------
	+
	+build:rivet:
	+ image: ${DOCKER_RIVET}
	+ extends: .HEJ_build
	+
	+test:rivet:
	+ image: ${DOCKER_RIVET}
	+ extends: .HEJ_test
	+ dependencies:
	+ - build:rivet
	+ script:
	+ - ctest --output-on-failure
	- bash -c '[ -f tst.yoda ]' && echo "found rivet output"
	- rivet-cmphistos tst.yoda
	- bash -c '[ -f MC_XS_XS.dat ]' && echo "yoda not empty"

	-QCDloop:
	- stage: different_compile
	- extends: .HEJ_compile
	- image: hejdock/qcdloopenv
	- variables:
	- HEJ_build_DIR: tmp_HEJ/build_HEJ_qcd
	- script:
	- - make test
	-
	-HepMC3:
	- stage: different_compile
	- extends: .HEJ_compile
	- image: hejdock/hepmc3env
	- variables:
	- HEJ_build_DIR: tmp_HEJ/build_HEJ_HepMC3
	- script:
	- - make test
	+# ---------------------------------
	+# - Clean Code -
	+# ---------------------------------

	No_tabs:
	stage: clean_code
	tags: *tags_def
	image: hejdock/git
	+ dependencies: []
	script:
	- date
	- check_tabs

	-No_warnings:
	+# ----------- No gcc warnings -----------
	+
	+.Warning_build:
	+ extends: .HEJ_build
	stage: clean_code
	- tags: *tags_def
	- image: hejdock/hepenv
	- before_script:
	- - date
	- - source /cvmfs/pheno.egi.eu/HEJ/HEJ_env.sh \|\| exit 1
	- - export t_HEJ_Install_DIR=${PWD}/$HEJ_Install_DIR
	- - export t_HEJ_build_DIR=${PWD}/$HEJ_build_DIR
	- - export t_HEJ_DIR=${PWD}
	- - export LD_LIBRARY_PATH=${t_HEJ_Install_DIR}/lib:${LD_LIBRARY_PATH}
	- - export PATH=${t_HEJ_Install_DIR}/bin:${PATH}
	- - export t_FOG_build_DIR=${PWD}/$FOG_build_DIR
	- - export t_FOG_DIR=${t_HEJ_DIR}/FixedOrderGen
	- - mkdir -p ${t_HEJ_build_DIR}
	- - mkdir -p ${t_FOG_build_DIR}
	+ dependencies: []
	script:
	+ - cd ${t_HEJ_DIR}
	+ # suppress warnings from side packages
	- sed -i 's/include_directories(${LHAPDF/include_directories(SYSTEM ${LHAPDF/g' CMakeLists.txt
	- sed -i 's/include_directories(${fastjet/include_directories(SYSTEM ${fastjet/g' CMakeLists.txt
	- sed -i 's/include_directories(${Boost/include_directories(SYSTEM ${Boost/g' CMakeLists.txt
	- - cd ${t_HEJ_build_DIR}
	- - cmake ${t_HEJ_DIR} -DCMAKE_CXX_FLAGS="-Werror" -DCMAKE_INSTALL_PREFIX=${t_HEJ_Install_DIR}
	- - make -j 8
	- - make install
	+ - cd ${t_HEJ_BUILD_DIR}
	+ - cmake ${t_HEJ_DIR} -DCMAKE_CXX_FLAGS="-Werror"
	+ - make -j $(nproc --ignore=1)
	+ artifacts: # don't save anything
	+
	+.Warning_FOG:
	+ extends: .FOG_build
	+ stage: clean_code
	+ script:
	- cd ${t_FOG_DIR}
	+ # suppress warnings from side packages
	- sed -i 's/include_directories(${LHAPDF/include_directories(SYSTEM ${LHAPDF/g' CMakeLists.txt
	- sed -i 's/include_directories(${fastjet/include_directories(SYSTEM ${fastjet/g' CMakeLists.txt
	- sed -i 's/include_directories(${Boost/include_directories(SYSTEM ${Boost/g' CMakeLists.txt
	- - cd ${t_FOG_build_DIR}
	+ - cd ${t_FOG_BUILD_DIR}
	- cmake ${t_FOG_DIR} -DCMAKE_CXX_FLAGS="-Werror"
	- - make -j 8
	+ - make -j $(nproc --ignore=1)
	+ artifacts: # don't save anything
	+
	+No_Warning:basic:
	+ image: ${DOCKER_BASIC}
	+ extends: .Warning_build
	+
	+No_Warning:basic:FOG:
	+ image: ${DOCKER_BASIC}
	+ extends: .Warning_FOG
	+ dependencies:
	+ - build:basic

	-.Publish:
	+# ---------------------------------
	+# - Publish -
	+# ---------------------------------
	+
	+Publish_version:
	stage: publish
	tags: *tags_def
	image: hejdock/git
	+ dependencies: []
	+ only:
	+ - /^v\d+\.\d+$/
	+ except:
	+ - tags
	+ when: on_success
	before_script:
	- mkdir -p .ssh/ && echo "${SSH_KEY}" > .ssh/id_rsa && chmod 0600 .ssh/id_rsa
	- rm -rf ~/.ssh/id_rsa; mkdir -p ~/.ssh/
	- ln -s $PWD/.ssh/id_rsa ~/.ssh/id_rsa && chmod 0600 ~/.ssh/id_rsa
	- ssh -T ${PUBLIC_GIT_PREFIX} -o "StrictHostKeyChecking no" \|\| echo "added ssh"
	script:
	- git remote add public ${PUBLIC_GIT_PREFIX}${PUBLIC_GIT_POSTFIX}
	- git checkout $CI_COMMIT_REF_NAME
	- git branch
	- git pull
	- git push public
	- git push public --tags
	after_script:
	- rm -f /root/.ssh/id_rsa && rm -fr .ssh
	- git remote rm public
	-
	-Publish_version:
	- extends: .Publish
	- only:
	- - /^v\d+\.\d+$/
	- except:
	- - tags
	- when: on_success
	diff --git a/Changes.md b/Changes.md
	index d1f98d3..aea4b72 100644
	--- a/Changes.md
	+++ b/Changes.md
	@@ -1,31 +1,33 @@
	# Version 2.0

	## 2.X.0

	* Made `MatrixElement::tree_kin` and `MatrixElement::tree_param` public
	* Allow multiplication and division of multiple scale functions e.g.
	`H_T/2*m_j1j2`
	* Follow HepMC convention for particle Status codes: incoming = 11,
	decaying = 2, outgoing = 1 (unchanged)
	+* New class `CrossSectionAccumulator` to keep track of Cross Section of the
	+ different subproccess

	## 2.0.4

	* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
	* Correctly include rivet headers
	* New `EXCLUDE_package` variable in `cmake` to not interface to specific
	packages
	* Consistent search and include for side packages in `cmake`

	## 2.0.3

	* Fixed parsing of (numerical factor) * (base scale) in configuration
	* Don't change scale names, but sanitise Rivet output file names instead

	## 2.0.2

	* Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
	`"/"` and `"*"` before)

	## 2.0.1

	* Fixed name of fixed-order generator in error message.
	diff --git a/FixedOrderGen/include/CrossSectionAccumulator.hh b/FixedOrderGen/include/CrossSectionAccumulator.hh
	deleted file mode 100644
	index 00f20c6..0000000
	--- a/FixedOrderGen/include/CrossSectionAccumulator.hh
	+++ /dev/null
	@@ -1,43 +0,0 @@
	-#pragma once
	-
	-#include <map>
	-
	-#include "HEJ/Event.hh"
	-#include "HEJ/event_types.hh"
	-
	-namespace HEJFOG {
	- template<typename T>
	- struct WithError {
	- T value = T{};
	- T error = T{};
	- };
	-
	- class CrossSectionAccumulator {
	- public:
	- void fill(HEJ::Event const & ev) {
	- const double wt = ev.central().weight;
	- auto & entry = xs_[ev.type()];
	- entry.value += wt;
	- entry.error += wt*wt;
	- total_.value += wt;
	- total_.error += wt*wt;
	- }
	-
	- auto begin() const {
	- return std::begin(xs_);
	- }
	-
	- auto end() const {
	- return std::end(xs_);
	- }
	-
	- WithError<double> total() const {
	- return total_;
	- }
	-
	- private:
	- std::map<HEJ::event_type::EventType, WithError<double>> xs_;
	- WithError<double> total_;
	- };
	-
	-}
	diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
	index 7e0cdf3..6d55edb 100644
	--- a/FixedOrderGen/src/main.cc
	+++ b/FixedOrderGen/src/main.cc
	@@ -1,246 +1,233 @@
	/**
	* Name: main.cc
	* Authors: Jeppe R. Andersen
	*/

	#include <algorithm>
	#include <chrono>
	#include <fstream>
	#include <iostream>
	#include <map>
	#include <memory>

	#include "yaml-cpp/yaml.h"

	#include "LHEF/LHEF.h"

	#include "HEJ/CombinedEventWriter.hh"
	+#include "HEJ/CrossSectionAccumulator.hh"
	#include "HEJ/get_analysis.hh"
	#include "HEJ/LesHouchesWriter.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/ProgressBar.hh"
	#include "HEJ/stream.hh"

	#include "config.hh"
	-#include "CrossSectionAccumulator.hh"
	#include "EventGenerator.hh"
	#include "PhaseSpacePoint.hh"
	#include "Unweighter.hh"
	#include "Version.hh"

	namespace{
	constexpr auto banner =
	" __ ___ __ ______ __ __ \n"
	" / / / (_)___ _/ /_ / ____/___ ___ _________ ___ __ / /__ / /______ \n"
	" / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _ \\/ __/ ___/ \n"
	" / __ / / /_/ / / / / / /___/ / / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) \n"
	" /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\____/\\___/\\__/____/ \n"
	" ____///__/ __ ____ ///__//____/ ______ __ \n"
	" / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / ____/__ ____ ___ _________ _/ /_____ _____\n"
	" / /_ / / \|/_/ _ \\/ __ / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ __/ __ \\/ ___/\n"
	" / __/ / /> </ __/ /_/ / / /_/ / / / /_/ / __/ / / /_/ / __/ / / / __/ / / /_/ / /_/ /_/ / / \n"
	" /_/ /_/_/\|_\|\\___/\\__,_/ \\____/_/ \\__,_/\\___/_/ \\____/\\___/_/ /_/\\___/_/ \\__,_/\\__/\\____/_/ \n"
	;

	constexpr double invGeV2_to_pb = 389379292.;
	constexpr long long max_warmup_events = 10000;
	}

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

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

	void rescale_weights(HEJ::Event & ev, double factor) {
	ev.central().weight *= factor;
	for(auto & var: ev.variations()){
	var.weight *= factor;
	}
	}

	int main(int argn, char** argv) {
	using namespace std::string_literals;
	if (argn < 2) {
	std::cerr << "\n# Usage:\n." << argv[0] << " config_file\n";
	return EXIT_FAILURE;
	}
	std::cout << banner;
	std::cout << "Version " << HEJFOG::Version::String()
	<< ", revision " << HEJFOG::Version::revision() << std::endl;
	fastjet::ClusterSequence::print_banner();
	using clock = std::chrono::system_clock;

	const auto start_time = clock::now();

	// read configuration
	auto config = load_config(argv[1]);

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

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

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

	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	std::move(scale_gen),
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	*ran
	};

	LHEF::HEPRUP heprup;
	heprup.IDBMUP=std::pair<long,long>(config.beam.particles[0], config.beam.particles[1]);
	heprup.EBMUP=std::make_pair(config.beam.energy, config.beam.energy);
	heprup.PDFGUP=std::make_pair(0,0);
	heprup.PDFSUP=std::make_pair(config.pdf_id,config.pdf_id);
	heprup.NPRUP=1;
	heprup.XSECUP=std::vector<double>(1.);
	heprup.XERRUP=std::vector<double>(1.);
	heprup.LPRUP=std::vector<int>{1};
	heprup.generators.emplace_back(LHEF::XMLTag{});
	heprup.generators.back().name = HEJFOG::Version::package_name();
	heprup.generators.back().version = HEJFOG::Version::String();

	HEJ::CombinedEventWriter writer{config.output, heprup};

	HEJ::optional<HEJFOG::Unweighter> unweighter{};
	std::map<HEJFOG::Status, int> status_counter;

	std::vector<HEJ::Event> events;
	int trials = 0;
	// warm-up phase to train unweighter
	if(config.unweight) {
	std::cout << "Calibrating unweighting ...\n";
	const auto warmup_start = clock::now();
	const size_t warmup_events = config.unweight->sample_size;
	HEJ::ProgressBar<size_t> warmup_progress{std::cout, warmup_events};
	for(; events.size() < warmup_events; ++trials){
	auto ev = generator.gen_event();
	++status_counter[generator.status()];
	if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
	events.emplace_back(std::move(ev));
	++warmup_progress;
	}
	}
	std::cout << std::endl;
	unweighter = HEJFOG::Unweighter{
	begin(events), end(events), config.unweight->max_dev, *ran,
	config.jets.peak_pt?(*config.jets.peak_pt):0.
	};
	std::vector<HEJ::Event> unweighted_events;
	for(auto && ev: events) {
	auto unweighted = unweighter->unweight(std::move(ev));
	if(unweighted) {
	unweighted_events.emplace_back(std::move(*unweighted));
	}
	}
	events = std::move(unweighted_events);
	if(events.empty()) {
	std::cerr <<
	"Failed to generate events. Please increase \"unweight: sample size\""
	" or reduce \"unweight: max deviation\"\n";
	return EXIT_FAILURE;
	}
	const auto warmup_end = clock::now();
	const double completion = static_cast<double>(events.size())/config.events;
	const std::chrono::duration<double> remaining_time =
	(warmup_end- warmup_start)*(1./completion - 1);
	const auto finish = clock::to_time_t(
	std::chrono::time_point_cast<std::chrono::seconds>(warmup_end + remaining_time)
	);
	std::cout
	<< "Generated " << events.size() << "/" << config.events << " events ("
	<< static_cast<int>(std::round(100*completion)) << "%)\n"
	<< "Estimated remaining generation time: "
	<< remaining_time.count() << " seconds ("
	<< std::put_time(std::localtime(&finish), "%c") << ")\n\n";
	}

	HEJ::ProgressBar<long long> progress{std::cout, config.events};
	progress.increment(events.size());
	events.reserve(config.events);
	for(; events.size() < static_cast<size_t>(config.events); ++trials){
	auto ev = generator.gen_event();
	++status_counter[generator.status()];
	if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
	if(unweighter) {
	auto unweighted = unweighter->unweight(std::move(ev));
	if(! unweighted) continue;
	ev = std::move(*unweighted);
	}
	events.emplace_back(std::move(ev));
	++progress;
	}
	}
	std::cout << std::endl;

	- HEJFOG::CrossSectionAccumulator xs;
	+ HEJ::CrossSectionAccumulator xs;
	for(auto & ev: events){
	rescale_weights(ev, invGeV2_to_pb/trials);
	analysis->fill(ev, ev);
	writer.write(ev);
	xs.fill(ev);
	}
	analysis->finalise();

	const std::chrono::duration<double> run_time = (clock::now() - start_time);
	std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n\n";

	- std::cout.precision(10);
	- std::cout.setf(std::ios::fixed);
	- std::cout
	- << " " << std::left << std::setw(20)
	- << "Cross section: " << xs.total().value
	- << " +- " << std::sqrt(xs.total().error) << " (pb)\n";
	- for(auto const & xs_type: xs) {
	- std::cout
	- << " " << std::left << std::setw(20)
	- << (HEJ::event_type::names[xs_type.first] + ": "s)
	- << xs_type.second.value << " +- "
	- << std::sqrt(xs_type.second.error) << " (pb)\n";
	- }
	+ std::cout << xs << '\n';

	- std::cout << '\n';
	for(auto && entry: status_counter){
	const double fraction = static_cast<double>(entry.second)/trials;
	const int percent = std::round(100*fraction);
	std::cout << "status "
	<< std::left << std::setw(16) << (to_string(entry.first) + ":")
	<< " [";
	for(int i = 0; i < percent/2; ++i) std::cout << '#';
	for(int i = percent/2; i < 50; ++i) std::cout << ' ';
	std::cout << "] " << percent << "%\n";
	}

	}
	diff --git a/doc/developer_manual/developer_manual.tex b/doc/developer_manual/developer_manual.tex
	index aa0435f..de04414 100644
	--- a/doc/developer_manual/developer_manual.tex
	+++ b/doc/developer_manual/developer_manual.tex
	@@ -1,1385 +1,1436 @@
	\documentclass[a4paper,11pt]{article}

	\usepackage{fourier}
	\usepackage[T1]{fontenc}
	\usepackage{microtype}
	\usepackage{geometry}
	\usepackage{enumitem}
	\setlist[description]{leftmargin=\parindent,labelindent=\parindent}
	\usepackage{amsmath}
	\usepackage{amssymb}
	\usepackage[utf8x]{inputenc}
	\usepackage{graphicx}
	\usepackage{xcolor}
	\usepackage{todonotes}
	\usepackage{listings}
	\usepackage{xspace}
	\usepackage{tikz}
	\usetikzlibrary{arrows.meta}
	\usetikzlibrary{shapes}
	\usetikzlibrary{calc}
	\usepackage[colorlinks,linkcolor={blue!50!black}]{hyperref}
	\graphicspath{{build/figures/}{figures/}}
	\emergencystretch \hsize

	\newcommand{\HEJ}{{\tt HEJ}\xspace}
	\newcommand{\HIGHEJ}{\emph{High Energy Jets}\xspace}
	\newcommand{\cmake}{\href{https://cmake.org/}{cmake}\xspace}
	\newcommand{\html}{\href{https://www.w3.org/html/}{html}\xspace}
	\newcommand{\YAML}{\href{http://yaml.org/}{YAML}\xspace}
	\newcommand{\QCDloop}{\href{https://github.com/scarrazza/qcdloop}{QCDloop}\xspace}

	\newcommand{\as}{\alpha_s}
	\DeclareRobustCommand{\mathgraphics}[1]{\vcenter{\hbox{\includegraphics{#1}}}}
	\def\spa#1.#2{\left\langle#1\,#2\right\rangle}
	\def\spb#1.#2{\left[#1\,#2\right]} \def\spaa#1.#2.#3{\langle\mskip-1mu{#1} \|
	#2 \| {#3}\mskip-1mu\rangle} \def\spbb#1.#2.#3{[\mskip-1mu{#1} \| #2 \|
	{#3}\mskip-1mu]} \def\spab#1.#2.#3{\langle\mskip-1mu{#1} \| #2 \|
	{#3}\mskip-1mu\rangle} \def\spba#1.#2.#3{\langle\mskip-1mu{#1}^+ \| #2 \|
	{#3}^+\mskip-1mu\rangle} \def\spav#1.#2.#3{\\|\mskip-1mu{#1} \| #2 \|
	{#3}\mskip-1mu\\|^2} \def\jc#1.#2.#3{j^{#1}_{#2#3}}

	\definecolor{darkgreen}{rgb}{0,0.4,0}
	\lstset{ %
	backgroundcolor=\color{lightgray}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}
	basicstyle=\footnotesize\usefont{T1}{DejaVuSansMono-TLF}{m}{n}, % the size of the fonts that are used for the code
	breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace
	breaklines=false, % sets automatic line breaking
	captionpos=t, % sets the caption-position to bottom
	commentstyle=\color{red}, % comment style
	deletekeywords={...}, % if you want to delete keywords from the given language
	escapeinside={\%}{)}, % if you want to add LaTeX within your code
	extendedchars=true, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8
	frame=false, % adds a frame around the code
	keepspaces=true, % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible)
	keywordstyle=\color{blue}, % keyword style
	otherkeywords={}, % if you want to add more keywords to the set
	numbers=none, % where to put the line-numbers; possible values are (none, left, right)
	numbersep=5pt, % how far the line-numbers are from the code
	rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here))
	showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces'
	showstringspaces=false, % underline spaces within strings only
	showtabs=false, % show tabs within strings adding particular underscores
	stepnumber=2, % the step between two line-numbers. If it's 1, each line will be numbered
	stringstyle=\color{gray}, % string literal style
	tabsize=2, % sets default tabsize to 2 spaces
	title=\lstname,
	emph={},
	emphstyle=\color{darkgreen}
	}

	\begin{document}

	\tikzstyle{mynode}=[rectangle split,rectangle split parts=2, draw,rectangle split part fill={lightgray, none}]


	\title{HEJ 2 developer manual}
	\author{}
	\maketitle
	\tableofcontents

	\newpage

	\section{Overview}
	\label{sec:overview}

	HEJ 2 is a C++ program and library implementing an algorithm to
	apply \HIGHEJ resummation~\cite{Andersen:2008ue,Andersen:2008gc} to
	pre-generated fixed-order events. This document is intended to give an
	overview over the concepts and structure of this implementation.

	\subsection{Project structure}
	\label{sec:project}

	HEJ 2 is developed under the \href{https://git-scm.com/}{git}
	version control system. The main repository is on the IPPP
	\href{https://gitlab.com/}{gitlab} server under
	\url{https://gitlab.dur.scotgrid.ac.uk/hej/hej}. To get a local
	copy, get an account on the gitlab server and use
	\begin{lstlisting}[language=sh,caption={}]
	git clone git@gitlab.dur.scotgrid.ac.uk:hej/hej.git
	\end{lstlisting}
	This should create a directory \texttt{hej} with the following
	contents:
	\begin{description}
	\item[doc:] Contains additional documentation, see section~\ref{sec:doc}.
	\item[include:] Contains the C++ header files.
	\item[src:] Contains the C++ source files.
	\item[t:] Contains the source code for the automated tests.
	\item[CMakeLists.txt:] Configuration file for the \cmake build
	system. See section~\ref{sec:cmake}.
	\item[cmake:] Auxiliary files for \cmake. This includes modules for
	finding installed software in \texttt{cmake/Modules} and templates for
	code generation during the build process in \texttt{cmake/Templates}.
	\item[config.yml:] Sample configuration file for running HEJ 2.
	\item[FixedOrderGen:] Contains the code for the fixed-order generator,
	see section~\ref{sec:HEJFOG}.
	\end{description}
	In the following all paths are given relative to the
	\texttt{hej} directory.

	\subsection{Documentation}
	\label{sec:doc}

	The \texttt{doc} directory contains user documentation in
	\texttt{doc/sphinx} and the configuration to generate source code
	documentation in \texttt{doc/doxygen}.

	The user documentation explains how to install and run HEJ 2. The
	format is
	\href{http://docutils.sourceforge.net/rst.html}{reStructuredText}, which
	is mostly human-readable. Other formats, like \html, can be generated with the
	\href{http://www.sphinx-doc.org/en/master/}{sphinx} generator with
	\begin{lstlisting}[language=sh,caption={}]
	make html
	\end{lstlisting}


	To document the source code we use
	\href{https://www.stack.nl/~dimitri/doxygen/}{doxygen}. To generate
	\html documentation, use the command
	\begin{lstlisting}[language=sh,caption={}]
	doxygen Doxyfile
	\end{lstlisting}
	in the \texttt{doc/doxygen} directory.

	\subsection{Build system}
	\label{sec:cmake}

	For the most part, HEJ 2 is a library providing classes and
	functions that can be used to add resummation to fixed-order events. In
	addition, there is a relatively small executable program leveraging this
	library to read in events from an input file and produce resummation
	events. Both the library and the program are built and installed with
	the help of \cmake.
	Debug information can be turned on by using
	\begin{lstlisting}[language=sh,caption={}]
	cmake base/directory -DCMAKE_BUILD_TYPE=Debug
	make install
	\end{lstlisting}
	This facilitates the use of debuggers like \href{https://www.gnu.org/software/gdb/}{gdb}.

	The main \cmake configuration file is \texttt{CMakeLists.txt}. It defines the
	compiler flags, software prerequisites, header and source files used to
	build HEJ 2, and the automated tests.
	\texttt{cmake/Modules} contains module files that help with the
	detection of the software prerequisites and \texttt{cmake/Templates}
	template files for the automatic generation of header and
	source files. For example, this allows to only keep the version
	information in one central location (\texttt{CMakeLists.txt}) and
	automatically generate a header file from the template \texttt{Version.hh.in} to propagate this to the C++ code.

	\subsection{General coding guidelines}
	\label{sec:notes}

	The goal is to make the HEJ 2 code well-structured and
	readable. Here are a number of guidelines to this end.

	\begin{description}
	\item[Observe the boy scout rule.] Always leave the code cleaner
	than how you found it. Ugly hacks can be useful for testing, but
	shouldn't make their way into the main branch.
	\item[Ask if something is unclear.] Often there is a good reason why
	code is written the way it is. Sometimes that reason is only obvious to
	the original author (use \lstinline!git blame! to find them), in which
	case they should be poked to add a comment. Sometimes there is no good
	reason, but nobody has had the time to come up with something better,
	yet. In some places the code might just be bad.
	\item[Don't break tests.] There are a number of tests in the \texttt{t}
	directory, which can be run with \lstinline!make test!. Ideally, all
	tests should run successfully in each git revision. If your latest
	commit broke a test and you haven't pushed to the central repository
	yet, you can fix it with \lstinline!git commit --amend!. If an earlier
	local commit broke a test, you can use \lstinline!git rebase -i! if
	- you feel confident.
	+ you feel confident. Additionally each \lstinline!git push! is also
	+ automatically tested via the GitLab CI (see appendix~\ref{sec:gitlabCI}).
	\item[Test your new code.] When you add some new functionality, also add an
	automated test. This can be useful even if you don't know the
	``correct'' result because it prevents the code from changing its behaviour
	silently in the future. \href{http://www.valgrind.org/}{valgrind} is a
	very useful tool to detect potential memory leaks.
	\item[Stick to the coding style.] It is somewhat easier to read code
	that has a uniform coding and indentation style. We don't have a
	strict style, but it helps if your code looks similar to what is
	already there.
	\end{description}

	\section{Program flow}
	\label{sec:flow}

	A run of the HEJ 2 program has three stages: initialisation,
	event processing, and cleanup. The following sections outline these
	stages and their relations to the various classes and functions in the
	code. Unless denoted otherwise, all classes and functions are part of
	the \lstinline!HEJ! namespace. The code for the HEJ 2 program is
	in \texttt{src/bin/HEJ.cc}, all other code comprises the HEJ 2
	library. Classes and free functions are usually implemented in header
	and source files with a corresponding name, i.e. the code for
	\lstinline!MyClass! can usually be found in
	\texttt{include/HEJ/MyClass.hh} and \texttt{src/MyClass.cc}.

	\subsection{Initialisation}
	\label{sec:init}

	The first step is to load and parse the \YAML configuration file. The
	entry point for this is the \lstinline!load_config! function and the
	related code can be found in \texttt{include/HEJ/YAMLreader.hh},
	\texttt{include/HEJ/config.hh} and the corresponding \texttt{.cc} files
	in the \texttt{src} directory. The implementation is based on the
	\href{https://github.com/jbeder/yaml-cpp}{yaml-cpp} library.
	The \lstinline!load_config! function returns a \lstinline!Config! object
	containing all settings. To detect potential mistakes as early as
	possible, we throw an exception whenever one of the following errors
	occurs:
	\begin{itemize}
	\item There is an unknown option in the \YAML file.
	\item A setting is invalid, for example a string is given where a number
	would be expected.
	\item An option value is not set.
	\end{itemize}
	The third rule is sometimes relaxed for ``advanced'' settings with an
	obvious default, like for importing custom scales or analyses.

	The information stored in the \lstinline!Config! object is then used to
	initialise various objects required for the event processing stage
	described in section~\ref{sec:processing}. First, the
	\lstinline!get_analysis! function creates an object that inherits from
	the \lstinline!Analysis! interface.\footnote{In the context of C++ the
	proper technical expression is ``pure abstract class''.} Using an
	interface allows us to decide the concrete type of the analysis at run
	time instead of having to make a compile-time decision. Depending on the
	settings, \lstinline!get_analysis! creates either a user-defined
	analysis loaded from an external library (see the user documentation
	\url{https://hej.web.cern.ch/HEJ/doc/current/user/}) or the default \lstinline!EmptyAnalysis!, which does
	nothing.

	Together with a number of further objects, whose roles are described in
	section~\ref{sec:processing}, we also initialise the global random
	number generator. We again use an interface to defer deciding the
	concrete type until the program is actually run. Currently, we support the
	\href{https://mixmax.hepforge.org/}{MIXMAX}
	(\texttt{include/HEJ/Mixmax.hh}) and Ranlux64
	(\texttt{include/HEJ/Ranlux64.hh}) random number generators, both are provided
	by \href{http://proj-clhep.web.cern.ch/}{CLHEP}.

	We also set up a \lstinline!LHEF::Reader! object (see
	\href{http://home.thep.lu.se/~leif/LHEF/}{\texttt{include/LHEF/LHEF.h}}) for
	reading events from a file in the Les
	Houches event file format~\cite{Alwall:2006yp}. A small wrapper around
	the
	\href{https://www.boost.org/doc/libs/1_67_0/libs/iostreams/doc/index.html}{boost
	iostreams} library allows us to also read event files compressed with
	\href{https://www.gnu.org/software/gzip/}{gzip}. The wrapper code is in
	\texttt{include/HEJ/stream.hh} and the \texttt{src/stream.cc}.

	\subsection{Event processing}
	\label{sec:processing}

	In the second stage events are continously read from the event
	file. After jet clustering, a number of corresponding resummation events
	are generated for each input event and fed into the analysis and a
	number of output files. The roles of various classes and functions are
	illustrated in the following flow chart:
	\begin{center}
	\begin{tikzpicture}[node distance=2cm and 5mm]
	\node (reader) [mynode]
	{\lstinline!LHEF::Reader::readEvent!\nodepart{second}{read event}};
	\node
	(cluster) [mynode,below=of reader]
	{\lstinline!Event! constructor\nodepart{second}{cluster jets}};
	\node
	(resum) [mynode,below=of cluster]
	{\lstinline!EventReweighter::reweight!\nodepart{second}{perform resummation}};
	\node
	(cut) [mynode,below=of resum]
	{\lstinline!Analysis::pass_cuts!\nodepart{second}{apply cuts}};
	\node
	(fill) [mynode,below left=of cut]
	{\lstinline!Analysis::fill!\nodepart{second}{analyse event}};
	\node
	(write) [mynode,below right=of cut]
	{\lstinline!CombinedEventWriter::write!\nodepart{second}{write out event}};
	\node
	(control) [below=of cut] {};
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(reader.south) -- node[left] {\lstinline!LHEF::HEPEUP!} (cluster.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(cluster.south) -- node[left] {\lstinline!Event!} (resum.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(resum.south) -- (cut.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(resum.south)+(7mm, 0cm)$) -- ($(cut.north)+(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(resum.south)-(7mm, 0cm)$) -- node[left] {\lstinline!Event!} ($(cut.north)-(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(cut.south)-(3mm,0mm)$) .. controls ($(control)-(3mm,0mm)$) ..node[left] {\lstinline!Event!} (fill.east);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(cut.south)-(3mm,0mm)$) .. controls ($(control)-(3mm,0mm)$) .. (write.west);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(cut.south)+(3mm,0mm)$) .. controls ($(control)+(3mm,0mm)$) .. (fill.east);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(cut.south)+(3mm,0mm)$) .. controls ($(control)+(3mm,0mm)$) ..node[right] {\lstinline!Event!} (write.west);
	\end{tikzpicture}
	\end{center}
	The resummation is performed by the \lstinline!EventReweighter! class,
	which is described in more detail in section~\ref{sec:resum}. The
	\lstinline!CombinedEventWriter! writes events to zero or more output
	files. To this end, it contains a number of objects implementing the
	\lstinline!EventWriter! interface. These event writers typically write
	the events to a file in a given format. We currently have the
	\lstinline!LesHouchesWriter! for event files in the Les Houches Event
	File format and the \lstinline!HepMCWriter! for the
	\href{https://hepmc.web.cern.ch/hepmc/}{HepMC} format (Version 2 and 3).

	\subsection{Resummation}
	\label{sec:resum}

	In the \lstinline!EventReweighter::reweight! member function, we first
	classify the input fixed-order event (FKL, unordered, non-HEJ, \dots)
	and decide according to the user settings whether to discard, keep, or
	resum the event. If we perform resummation for the given event, we
	generate a number of trial \lstinline!PhaseSpacePoint! objects. Phase
	space generation is discussed in more detail in
	section~\ref{sec:pspgen}. We then perform jet clustering according to
	the settings for the resummation jets on each
	\lstinline!PhaseSpacePoint!, update the factorisation and
	renormalisation scale in the resulting \lstinline!Event! and reweight it
	according to the ratio of pdf factors and \HEJ matrix elements between
	resummation and original fixed-order event:
	\begin{center}
	\begin{tikzpicture}[node distance=2cm and 5mm]
	\node (in) {};
	\node (treat) [diamond,draw,below=of in,minimum size=3.5cm,
	label={[anchor=west, inner sep=8pt]west:discard},
	label={[anchor=east, inner sep=14pt]east:keep},
	label={[anchor=south, inner sep=20pt]south:reweight}
	] {};
	\draw (treat.north west) -- (treat.south east);
	\draw (treat.north east) -- (treat.south west);
	\node
	(psp) [mynode,below=of treat]
	{\lstinline!PhaseSpacePoint! constructor};
	\node
	(cluster) [mynode,below=of psp]
	{\lstinline!Event! constructor\nodepart{second}{cluster jets}};
	\node
	(gen_scales) [mynode,below=of cluster]
	{\lstinline!ScaleGenerator::operator()!\nodepart{second}{update scales}};
	\node
	(rescale) [mynode,below=of gen_scales]
	{\lstinline!PDF::pdfpt!,
	\lstinline!MatrixElement!\nodepart{second}{reweight}};
	\node (out) [below of=rescale] {};
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(in.south) -- node[left] {\lstinline!Event!} (treat.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(treat.south) -- node[left] {\lstinline!Event!} (psp.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(psp.south) -- (cluster.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(psp.south)+(7mm, 0cm)$) -- ($(cluster.north)+(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(psp.south)-(7mm, 0cm)$) -- node[left]
	{\lstinline!PhaseSpacePoint!} ($(cluster.north)-(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(cluster.south) -- (gen_scales.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(cluster.south)+(7mm, 0cm)$) -- ($(gen_scales.north)+(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(cluster.south)-(7mm, 0cm)$) -- node[left]
	{\lstinline!Event!} ($(gen_scales.north)-(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(gen_scales.south) -- (rescale.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(gen_scales.south)+(7mm, 0cm)$) -- ($(rescale.north)+(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(gen_scales.south)-(7mm, 0cm)$) -- node[left]
	{\lstinline!Event!} ($(rescale.north)-(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(rescale.south) -- (out.north);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(rescale.south)+(7mm, 0cm)$) -- ($(out.north)+(7mm, 0cm)$);
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	($(rescale.south)-(7mm, 0cm)$) -- node[left]
	{\lstinline!Event!} ($(out.north)-(7mm, 0cm)$);
	\node (helper) at ($(treat.east) + (15mm,0cm)$) {};
	\draw[-{Latex[length=3mm, width=1.5mm]}]
	(treat.east) -- ($(treat.east) + (15mm,0cm)$)
	-- node[left] {\lstinline!Event!} (helper \|- gen_scales.east) -- (gen_scales.east)
	;
	\end{tikzpicture}
	\end{center}


	\subsection{Phase space point generation}
	\label{sec:pspgen}

	The resummed and matched \HEJ cross section for pure jet production of
	FKL configurations is given by (cf. eq. (3) of~\cite{Andersen:2018tnm})
	\begin{align}
	\label{eq:resumdijetFKLmatched2}
	% \begin{split}
	\sigma&_{2j}^\mathrm{resum, match}=\sum_{f_1, f_2}\ \sum_m
	\prod_{j=1}^m\left(
	\int_{p_{j\perp}^B=0}^{p_{j\perp}^B=\infty}
	\frac{\mathrm{d}^2\mathbf{p}_{j\perp}^B}{(2\pi)^3}\ \int
	\frac{\mathrm{d} y_j^B}{2} \right) \
	(2\pi)^4\ \delta^{(2)}\!\!\left(\sum_{k=1}^{m}
	\mathbf{p}_{k\perp}^B\right)\nonumber\\
	&\times\ x_a^B\ f_{a, f_1}(x_a^B, Q_a^B)\ x_b^B\ f_{b, f_2}(x_b^B, Q_b^B)\
	\frac{\overline{\left\|\mathcal{M}_\text{LO}^{f_1f_2\to f_1g\cdots
	gf_2}\big(\big\{p^B_j\big\}\big)\right\|}^2}{(\hat {s}^B)^2}\nonumber\\
	& \times (2\pi)^{-4+3m}\ 2^m \nonumber\\
	&\times\ \sum_{n=2}^\infty\
	\int_{p_{1\perp}=p_{\perp,\mathrm{min}} }^{p_{1\perp}=\infty}
	\frac{\mathrm{d}^2\mathbf{p}_{1\perp}}{(2\pi)^3}\
	\int_{p_{n\perp}=p_{\perp,\mathrm{min}}}^{p_{n\perp}=\infty}
	\frac{\mathrm{d}^2\mathbf{p}_{n\perp}}{(2\pi)^3}\
	\prod_{i=2}^{n-1}\int_{p_{i\perp}=\lambda}^{p_{i\perp}=\infty}
	\frac{\mathrm{d}^2\mathbf{p}_{i\perp}}{(2\pi)^3}\ (2\pi)^4\ \delta^{(2)}\!\!\left(\sum_{k=1}^n
	\mathbf{p}_{k\perp}\right )\\
	&\times \ \mathbf{T}_y \prod_{i=1}^n
	\left(\int \frac{\mathrm{d} y_i}{2}\right)\
	\mathcal{O}_{mj}^e\
	\left(\prod_{l=1}^{m-1}\delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l}\perp}^B -
	\mathbf{j}_{l\perp})\right)\
	\left(\prod_{l=1}^m\delta(y^B_{\mathcal{J}_l}-y_{\mathcal{J}_l})\right)
	\ \mathcal{O}_{2j}(\{p_i\})\nonumber\\
	&\times \frac{(\hat{s}^B)^2}{\hat{s}^2}\ \frac{x_a f_{a,f_1}(x_a, Q_a)\ x_b f_{b,f_2}(x_b, Q_b)}{x_a^B\ f_{a,f_1}(x_a^B, Q_a^B)\ x_b^B\ f_{b,f_2}(x_b^B, Q_b^B)}\ \frac{\overline{\left\|\mathcal{M}_{\mathrm{HEJ}}^{f_1 f_2\to f_1 g\cdots
	gf_2}(\{ p_i\})\right\|}^2}{\overline{\left\|\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to f_1g\cdots
	gf_2}\big(\big\{p^B_j\big\}\big)\right\|}^{2}} \,.\nonumber
	% \end{split}
	\end{align}
	The first two lines correspond to the generation of the fixed-order
	input events with incoming partons $f_1, f_2$ and outgoing momenta
	$p_j^B$, where $\mathbf{p}_{j\perp}^B$ and $y_j^B$ denote the respective
	transverse momentum and rapidity. Note that, at leading order, these
	coincide with the fixed-order jet momenta $p_{\mathcal{J}_j}^B$.
	$f_{a,f_1}(x_a, Q_a),f_{b,f_2}(x_b, Q_b)$ are the pdf factors for the incoming partons with
	momentum fractions $x_a$ and $x_b$. The square of the partonic
	centre-of-mass energy is denoted by $\hat{s}^B$ and
	$\mathcal{M}_\text{LO}^{f_1f_2\to f_1g\cdots gf_2}$ is the
	leading-order matrix element.

	The third line is a factor accounting for the different multiplicities
	between fixed-order and resummation events. Lines four and five are
	the integration over the resummation phase space described in this
	section. $p_i$ are the momenta of the outgoing partons in resummation
	phase space. $\mathbf{T}_y$ denotes rapidity
	ordering and $\mathcal{O}_{mj}^e$ projects out the exclusive $m$-jet
	component. The relation between resummation and fixed-order momenta is
	fixed by the $\delta$ functions. The first sets each transverse fixed-order jet
	momentum to some function $\mathbf{j_{l\perp}}$ of the resummation
	momenta. The exact form is described in section~\ref{sec:ptj_res}. The second
	$\delta$ forces the rapidities of resummation and fixed-order jets to be
	the same. Finally, the last line is the reweighting of pdf and matrix
	element factors already shown in section~\ref{sec:resum}.

	There are two kinds of cut-off in the integration over the resummation
	partons. $\lambda$ is a technical cut-off connected to the cancellation
	of infrared divergencies between real and virtual corrections. Its
	numerical value is set in
	\texttt{include/HEJ/Constants.h}. $p_{\perp,\mathrm{min}}$ regulates
	and \emph{uncancelled} divergence in the extremal parton momenta. Its
	size is set by the user configuration \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJ.html#settings}.

	It is straightforward to generalise eq.~(\ref{eq:resumdijetFKLmatched2})
	to unordered configurations and processes with additional colourless
	emissions, for example a Higgs or electroweak boson. In the latter case only
	the fixed-order integration and the matrix elements change.

	\subsubsection{Gluon Multiplicity}
	\label{sec:psp_ng}

	The first step in evaluating the resummation phase space in
	eq.~(\ref{eq:resumdijetFKLmatched2}) is to randomly pick terms in the
	sum over the number of emissions. This sampling of the gluon
	multiplicity is done in the \lstinline!PhaseSpacePoint::sample_ng!
	function in \texttt{src/PhaseSpacePoint.cc}.

	The typical number of extra emissions depends strongly on the rapidity
	span of the underlying fixed-order event. Let us, for example, consider
	a fixed-order FKL-type multi-jet configuration with rapidities
	$y_{j_f},\,y_{j_b}$ of the most forward and backward jets,
	respectively. By eq.~(\ref{eq:resumdijetFKLmatched2}), the jet
	multiplicity and the rapidity of each jet are conserved when adding
	resummation. This implies that additional hard radiation is restricted
	to rapidities $y$ within a region $y_{j_b} \lesssim y \lesssim
	y_{j_f}$. Within \HEJ, we require the most forward and most backward
	emissions to be hard \todo{specify how hard} in order to avoid divergences, so this constraint
	in fact applies to \emph{all} additional radiation.

	To simplify the remaining discussion, let us remove the FKL rapidity
	ordering
	\begin{equation}
	\label{eq:remove_y_order}
	\mathbf{T}_y \prod_{i=1}^n\int \frac{\mathrm{d}y_i}{2} =
	\frac{1}{n!}\prod_{i=1}^n\int
	\frac{\mathrm{d}y_i}{2}\,,
	\end{equation}
	where all rapidity integrals now cover a region which is approximately
	bounded by $y_{j_b}$ and $y_{j_f}$. Each of the $m$ jets has to contain at least
	one parton; selecting random emissions we can rewrite the phase space
	integrals as
	\begin{equation}
	\label{eq:select_jets}
	\frac{1}{n!}\prod_{i=1}^n\int [\mathrm{d}p_i] =
	\left(\prod_{i=1}^{m}\int [\mathrm{d}p_i]\ {\cal J}_i(p_i)\right)
	\frac{1}{n_g!}\prod_{i=m+1}^{m+n_g}\int [\mathrm{d}p_i]
	\end{equation}
	with jet selection functions
	\begin{equation}
	\label{eq:def_jet_selection}
	{\cal J}_i(p) =
	\begin{cases}
	1 &p\text{ clustered into jet }i\\
	0 & \text{otherwise}
	\end{cases}
	\end{equation}
	and $n_g \equiv n - m$. Here and in the following we use the short-hand
	notation $[\mathrm{d}p_i]$ to denote the phase-space measure for parton
	$i$. As is evident from eq.~\eqref{eq:select_jets}, adding an extra emission
	$n_g+1$ introduces a suppression factor $\tfrac{1}{n_g+1}$. However, the
	additional phase space integral also results in an enhancement proportional
	to $\Delta y_{j_f j_b} = y_{j_f} - y_{j_b}$. This is a result of the
	rapidity-independence of the MRK limit of the integrand, consisting of the
	matrix elements divided by the flux factor. Indeed, we observe that the
	typical number of gluon emissions is to a good approximation proportional to
	the rapidity separation and the phase space integral is dominated by events
	with $n_g \approx \Delta y_{j_f j_b}$.

	For the actual phase space sampling, we assume a Poisson distribution
	and extract the mean number of gluon emissions in different rapidity
	bins and fit the results to a linear function in $\Delta y_{j_f j_b}$,
	finding a coefficient of $0.975$ for the inclusive production of a Higgs
	boson with two jets. Here are the observed and fitted average gluon
	multiplicities as a function of $\Delta y_{j_f j_b}$:
	\begin{center}
	\includegraphics[width=.75\textwidth]{ng_mean}
	\end{center}
	As shown for two rapidity slices the assumption of a Poisson
	distribution is also a good approximation:
	\begin{center}
	\includegraphics[width=.49\textwidth]{{ng_1.5}.pdf}\hfill
	\includegraphics[width=.49\textwidth]{{ng_5.5}.pdf}
	\end{center}

	\subsubsection{Number of Gluons inside Jets}
	\label{sec:psp_ng_jet}

	For each of the $n_g$ gluon emissions we can split the phase-space
	integral into a (disconnected) region inside the jets and a remainder:
	\begin{equation}
	\label{eq:psp_split}
	\int [\mathrm{d}p_i] = \int [\mathrm{d}p_i]\,
	\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg) + \int [\mathrm{d}p_i]\,
	\bigg[1-\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg)\bigg]\,.
	\end{equation}
	The next step is to decide how many of the gluons will form part of a
	jet. This is done in the \lstinline!PhaseSpacePoint::sample_ng_jets!
	function.

	We choose an importance sampling which is flat in the plane
	spanned by the azimuthal angle $\phi$ and the rapidity $y$. This is
	observed in BFKL and valid in the limit of Multi-Regge-Kinematics
	(MRK). Furthermore, we assume anti-$k_t$ jets, which cover an area of
	$\pi R^2$.

	In principle, the total accessible area in the $y$-$\phi$ plane is given
	by $2\pi \Delta y_{fb}$, where $\Delta y_{fb}\geq \Delta y_{j_f j_b}$ is
	the a priori unknown rapidity separation between the most forward and
	backward partons. In most cases the extremal jets consist of single
	partons, so that $\Delta y_{fb} = \Delta y_{j_f j_b}$. For the less common
	case of two partons forming a jet we observe a maximum distance of $R$
	between the constituents and the jet centre. In rare cases jets have
	more than two constituents. Empirically, they are always within a
	distance of $\tfrac{5}{3}R$ to the centre of the jet, so
	$\Delta y_{fb} \leq \Delta y_{j_f j_b} + \tfrac{10}{3} R$. In practice, the
	extremal partons are required to carry a large fraction of the jet
	transverse momentum and will therefore be much closer to the jet axis.

	In summary, for sufficiently large rapidity separations we can use the
	approximation $\Delta y_{fb} \approx \Delta y_{j_f j_b}$. This scenario
	is depicted here:
	\begin{center}
	\includegraphics[width=0.5\linewidth]{ps_large_y}
	\end{center}
	If there is no overlap between jets, the probability $p_{\cal J, >}$ for
	an extra gluon to end up inside a jet is then given by
	\begin{equation}
	\label{eq:p_J_large}
	p_{\cal J, >} = \frac{(m - 1)\*R^2}{2\Delta y_{j_f j_b}}\,.
	\end{equation}
	For a very small rapidity separation, eq.~\eqref{eq:p_J_large}
	obviously overestimates the true probability. The maximum phase space
	covered by jets in the limit of a vanishing rapidity distance between
	all partons is $2mR \Delta y_{fb}$:
	\begin{center}
	\includegraphics[width=0.5\linewidth]{ps_small_y}
	\end{center}
	We therefore estimate the probability for a parton to end up inside a jet as
	\begin{equation}
	\label{eq:p_J}
	p_{\cal J} = \min\bigg(\frac{(m - 1)\*R^2}{2\Delta y_{j_f j_b}}, \frac{mR}{\pi}\bigg)\,.
	\end{equation}
	Here we compare this estimate with the actually observed
	fraction of additional emissions into jets as a function of the rapidity
	separation:
	\begin{center}
	\includegraphics[width=0.75\linewidth]{pJ}
	\end{center}

	\subsubsection{Gluons outside Jets}
	\label{sec:gluons_nonjet}

	Using our estimate for the probability of a gluon to be a jet
	constituent, we choose a number $n_{g,{\cal J}}$ of gluons inside
	jets, which also fixes the number $n_g - n_{g,{\cal J}}$ of gluons
	outside jets. As explained later on, we need to generate the momenta of
	the gluons outside jets first. This is done in
	\lstinline!PhaseSpacePoint::gen_non_jet!.

	The azimuthal angle $\phi$ is generated flat within $0\leq \phi \leq 2
	\pi$. The allowed rapidity interval is set by the most forward and
	backward partons, which are necessarily inside jets. Since these parton
	rapidities are not known at this point, we also have to postpone the
	rapidity generation for the gluons outside jets. For the scalar
	transverse momentum $p_\perp = \|\mathbf{p}_\perp\|$ of a gluon outside
	jets we use the parametrisation
	\begin{equation}
	\label{eq:p_nonjet}
	p_\perp = \lambda + \tilde{p}_\perp\\tan(\tau\r)\,, \qquad
	\tau = \arctan\bigg(\frac{p_{\perp{\cal J}_\text{min}} - \lambda}{\tilde{p}_\perp}\bigg)\,.
	\end{equation}
	For $r \in [0,1)$, $p_\perp$ is always less than the minimum momentum
	$p_{\perp{\cal J}_\text{min}}$ required for a jet. $\tilde{p}_\perp$ is
	a free parameter, a good empirical value is $\tilde{p}_\perp = [1.3 +
	0.2\*(n_g - n_{g,\cal J})]\,$GeV

	\subsubsection{Resummation jet momenta}
	\label{sec:ptj_res}

	On the one hand, each jet momentum is given by the sum of its
	constituent momenta. On the other hand, the resummation jet momenta are
	fixed by the constraints in line five of the master
	equation~\eqref{eq:resumdijetFKLmatched2}. We therefore have to
	calculate the resummation jet momenta from these constraints before
	generating the momenta of the gluons inside jets. This is done in
	\lstinline!PhaseSpacePoint::reshuffle! and in the free
	\lstinline!resummation_jet_momenta! function (declared in \texttt{resummation\_jet.hh}).

	The resummation jet momenta are determined by the $\delta$ functions in
	line five of eq.~(\ref{eq:resumdijetFKLmatched2}). The rapidities are
	fixed to the rapidities of the jets in the input fixed-order events, so
	that the FKL ordering is guaranteed to be preserved.

	In traditional \HEJ reshuffling the transverse momentum are given through
	\begin{equation}
	\label{eq:ptreassign_old}
	\mathbf{p}^B_{\mathcal{J}_{l\perp}} = \mathbf{j}_{l\perp} \equiv \mathbf{p}_{\mathcal{J}_{l}\perp}
	+ \mathbf{q}_\perp \,\frac{\|\mathbf{p}_{\mathcal{J}_{l}\perp}\|}{P_\perp},
	\end{equation}
	where $\mathbf{q}_\perp = \sum_{j=1}^n \mathbf{p}_{i\perp}
	\bigg[1-\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg)\bigg] $ is the
	total transverse momentum of all partons \emph{outside} jets and
	$P_\perp = \sum_{j=1}^m \|\mathbf{p}_{\mathcal{J}_{j}\perp}\|$. Since the
	total transverse momentum of an event vanishes, we can also use
	$\mathbf{q}_\perp = - \sum_{j=1}^m
	\mathbf{p}_{\mathcal{J}_{j}\perp}$. Eq.~(\ref{eq:ptreassign}) is a
	non-linear system of equations in the resummation jet momenta
	$\mathbf{p}_{\mathcal{J}_{l}\perp}$. Hence we would have to solve
	\begin{equation}
	\label{eq:ptreassign_eq}
	\mathbf{p}_{\mathcal{J}_{l}\perp}=\mathbf{j}^B_{l\perp} \equiv\mathbf{j}_{l\perp}^{-1}
	\left(\mathbf{p}^B_{\mathcal{J}_{l\perp}}\right)
	\end{equation}
	numerically.

	Since solving such a system is computationally expensive, we instead
	change the reshuffling around to be linear in the resummation jet
	momenta. Hence~\eqref{eq:ptreassign_eq} gets replaces by
	\begin{equation}
	\label{eq:ptreassign}
	\mathbf{p}_{\mathcal{J}_{l\perp}} = \mathbf{j}^B_{l\perp} \equiv \mathbf{p}^B_{\mathcal{J}_{l}\perp}
	- \mathbf{q}_\perp \,\frac{\|\mathbf{p}^B_{\mathcal{J}_{l}\perp}\|}{P^B_\perp},
	\end{equation}
	which is linear in the resummation momentum. Consequently the equivalent
	of~\eqref{eq:ptreassign_old} is non-linear in the Born momentum. However
	the exact form of~\eqref{eq:ptreassign_old} is not relevant for the resummation.
	Both methods have been tested for two and three jets with the \textsc{rivet}
	standard analysis \texttt{MC\_JETS}. They didn't show any differences even
	after $10^9$ events.

	The reshuffling relation~\eqref{eq:ptreassign} allows the transverse
	momenta $p^B_{\mathcal{J}_{l\perp}}$ of the fixed-order jets to be
	somewhat below the minimum transverse momentum of resummation jets. It
	is crucial that this difference does not become too large, as the
	fixed-order cross section diverges for vanishing transverse momenta. In
	the production of a Higgs boson with resummation jets above $30\,$GeV we observe
	that the contribution from fixed-order events with jets softer than
	about $20\,$GeV can be safely neglected. This is shown in the following
	plot of the differential cross section over the transverse momentum of
	the softest fixed-order jet:
	\begin{center}
	\includegraphics[width=.75\textwidth]{ptBMin}
	\end{center}

	Finally, we have to account for the fact that the reshuffling
	relation~\eqref{eq:ptreassign} is non-linear in the Born momenta. To
	arrive at the master formula~\eqref{eq:resumdijetFKLmatched2} for the
	cross section, we have introduced unity in the form of an integral over
	the Born momenta with $\delta$ functions in the integrand, that is
	\begin{equation}
	\label{eq:delta_intro}
	1 = \int_{p_{j\perp}^B=0}^{p_{j\perp}^B=\infty}
	\mathrm{d}^2\mathbf{p}_{j\perp}^B\delta^{(2)}(\mathbf{p}_{\mathcal{J}_{j\perp}}^B -
	\mathbf{j}_{j\perp})\,.
	\end{equation}
	If the arguments of the $\delta$ functions are not linear in the Born
	momenta, we have to compensate with additional Jacobians as
	factors. Explicitly, for the reshuffling relation~\eqref{eq:ptreassign}
	we have
	\begin{equation}
	\label{eq:delta_rewrite}
	\prod_{l=1}^m \delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l\perp}}^B -
	\mathbf{j}_{l\perp}) = \Delta \prod_{l=1}^m \delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l\perp}} -
	\mathbf{j}_{l\perp}^B)\,,
	\end{equation}
	where $\mathbf{j}_{l\perp}^B$ is given by~\eqref{eq:ptreassign_eq} and only
	depends on the Born momenta. We have extended the product to run to $m$
	instead of $m-1$ by eliminating the last $\delta$ function
	$\delta^{(2)}\!\!\left(\sum_{k=1}^n \mathbf{p}_{k\perp}\right )$.
	The Jacobian $\Delta$ is the determinant of a $2m \times 2m$ matrix with $l, l' = 1,\dots,m$
	and $X, X' = x,y$.
	\begin{equation}
	\label{eq:jacobian}
	\Delta = \left\|\frac{\partial\,\mathbf{j}^B_{l'\perp}}{\partial\, \mathbf{p}^B_{{\cal J}_l \perp}} \right\|
	= \left\| \delta_{l l'} \delta_{X X'} - \frac{q_X\, p^B_{{\cal
	J}_{l'}X'}}{\left\|\mathbf{p}^B_{{\cal J}_{l'} \perp}\right\| P^B_\perp}\left(\delta_{l l'}
	- \frac{\left\|\mathbf{p}^B_{{\cal J}_l \perp}\right\|}{P^B_\perp}\right)\right\|\,.
	\end{equation}
	The determinant is calculated in \lstinline!resummation_jet_weight!,
	again coming from the \texttt{resummation\_jet.hh} header.
	Having to introduce this Jacobian is not a disadvantage specific to the new
	reshuffling. If we instead use the old reshuffling
	relation~\eqref{eq:ptreassign_old} we \emph{also} have to introduce a
	similar Jacobian since we actually want to integrate over the
	resummation phase space and need to transform the argument of the
	$\delta$ function to be linear in the resummation momenta for this.

	\subsubsection{Gluons inside Jets}
	\label{sec:gluons_jet}

	After the steps outlined in section~\ref{sec:psp_ng_jet}, we have a
	total number of $m + n_{g,{\cal J}}$ constituents. In
	\lstinline!PhaseSpacePoint::distribute_jet_partons! we distribute them
	randomly among the jets such that each jet has at least one
	constituent. We then generate their momenta in
	\lstinline!PhaseSpacePoint::split! using the \lstinline!Splitter! class.

	The phase space integral for a jet ${\cal J}$ is given by
	\begin{equation}
	\label{eq:ps_jetparton} \prod_{i\text{ in }{\cal J}} \bigg(\int
	\mathrm{d}\mathbf{p}_{i\perp}\ \int \mathrm{d} y_i
	\bigg)\delta^{(2)}\Big(\sum_{i\text{ in }{\cal J}} \mathbf{p}_{i\perp} -
	\mathbf{j}_{\perp}^B\Big)\delta(y_{\mathcal{J}}-y^B_{\mathcal{J}})\,.
	\end{equation}
	For jets with a single constituent, the parton momentum is obiously equal to the
	jet momentum. In the case of two constituents, we observe that the
	partons are always inside the jet cone with radius $R$ and often very
	close to the jet centre. The following plots show the typical relative
	distance $\Delta R/R$ for this scenario:
	\begin{center}
	\includegraphics[width=0.45\linewidth]{dR_2}
	\includegraphics[width=0.45\linewidth]{dR_2_small}
	\end{center}
	According to this preference for small values of $\Delta R$, we
	parametrise the $\Delta R$ integrals as
	\begin{equation}
	\label{eq:dR_sampling}
	\frac{\Delta R}{R} =
	\begin{cases}
	0.25\,x_R & x_R < 0.4 \\
	1.5\,x_R - 0.5 & x_R \geq 0.4
	\end{cases}\,.
	\end{equation}
	Next, we generate $\Theta_1 \equiv \Theta$ and use the constraint $\Theta_2 = \Theta
	\pm \pi$. The transverse momentum of the first parton is then given by
	\begin{equation}
	\label{eq:delta_constraints}
	p_{1\perp} =
	\frac{p_{\mathcal{J} y} - \tan(\phi_2) p_{\mathcal{J} x}}{\sin(\phi_1)
	- \tan(\phi_2)\cos(\phi_1)}\,.
	\end{equation}
	We get $p_{2\perp}$ by exchanging $1 \leftrightarrow 2$ in the
	indices. To obtain the Jacobian of the transformation, we start from the
	single jet phase space eq.~(\ref{eq:ps_jetparton}) with the rapidity
	delta function already rewritten to be linear in the rapidity of the
	last parton, i.e.
	\begin{equation}
	\label{eq:jet_2p}
	\prod_{i=1,2} \bigg(\int
	\mathrm{d}\mathbf{p}_{i\perp}\ \int \mathrm{d} y_i
	\bigg)\delta^{(2)}\Big(\mathbf{p}_{1\perp} + \mathbf{p}_{2\perp} -
	\mathbf{j}_{\perp}^B\Big)\delta(y_2- \dots)\,.
	\end{equation}
	The integral over the second parton momentum is now trivial; we can just replace
	the integral over $y_2$ with the equivalent constraint
	\begin{equation}
	\label{eq:R2}
	\int \mathrm{d}R_2 \ \delta\bigg(R_2 - \bigg[\phi_{\cal J} - \arctan
	\bigg(\frac{p_{{\cal J}y} - p_{1y}}{p_{{\cal J}x} -
	p_{1x}}\bigg)\bigg]/\cos \Theta\bigg) \,.
	\end{equation}
	In order to fix the integral over $p_{1\perp}$ instead, we rewrite this
	$\delta$ function. This introduces the Jacobian
	\begin{equation}
	\label{eq:jac_pt1}
	\bigg\|\frac{\partial p_{1\perp}}{\partial R_2} \bigg\| =
	\frac{\cos(\Theta)\mathbf{p}_{2\perp}^2}{p_{{\cal J}\perp}\sin(\phi_{\cal J}-\phi_1)}\,.
	\end{equation}
	The final form of the integral over the two parton momenta is then
	\begin{equation}
	\label{eq:ps_jet_2p}
	\int \mathrm{d}R_1\ R_1 \int \mathrm{d}R_2 \int \mathrm{d}x_\Theta\ 2\pi \int
	\mathrm{d}p_{1\perp}\ p_{1\perp} \int \mathrm{d}p_{2\perp}
	\ \bigg\|\frac{\partial p_{1\perp}}{\partial R_2} \bigg\|\delta(p_{1\perp}
	-\dots) \delta(p_{2\perp} - \dots)\,.
	\end{equation}

	As is evident from section~\ref{sec:psp_ng_jet}, jets with three or more
	constituents are rare and an efficient phase-space sampling is less
	important. For such jets, we exploit the observation that partons with a
	distance larger than $R_{\text{max}} = \tfrac{5}{3} R$ to
	the jet centre are never clustered into the jet. Assuming $N$
	constituents, we generate all components
	for the first $N-1$ partons and fix the remaining parton with the
	$\delta$-functional. In order to end up inside the jet, we use the
	parametrisation
	\begin{align}
	\label{eq:ps_jet_param}
	\phi_i ={}& \phi_{\cal J} + \Delta \phi_i\,, & \Delta \phi_i ={}& \Delta
	R_i
	\cos(\Theta_i)\,, \\
	y_i ={}& y_{\cal J} + \Delta y_i\,, & \Delta y_i ={}& \Delta
	R_i
	\sin(\Theta_i)\,,
	\end{align}
	and generate $\Theta_i$ and $\Delta R_i$ randomly with $\Delta R_i \leq
	R_{\text{max}}$ and the empiric value $R_{\text{max}} = 5\*R/3$. We can
	then write the phase space integral for a single parton as $(p_\perp = \|\mathbf{p}_\perp\|)$
	\begin{equation}
	\label{eq:ps_jetparton_x}
	\int \mathrm{d}\mathbf{p}_{\perp}\ \int
	\mathrm{d} y \approx \int_{\Box} \mathrm{d}x_{\perp}
	\mathrm{d}x_{ R}
	\mathrm{d}x_{\theta}\
	2\\pi\,\R_{\text{max}}^2\,\x_{R}\,\p_{\perp}\,\*(p_{\perp,\text{max}}
	- p_{\perp,\text{min}})
	\end{equation}
	with
	\begin{align}
	\label{eq:ps_jetparton_parameters}
	\Delta \phi ={}& R_{\text{max}}\x_{R}\\cos(2\\pi\x_\theta)\,,&
	\Delta y ={}& R_{\text{max}}\x_{R}\\sin(2\\pi\x_\theta)\,, \\
	p_{\perp} ={}& (p_{\perp,\text{max}} - p_{\perp,\text{min}})\*x_\perp +
	p_{\perp,\text{min}}\,.
	\end{align}
	$p_{\perp,\text{max}}$ is determined from the requirement that the total
	contribution from the first $n-1$ partons --- i.e. the projection onto the
	jet $p_{\perp}$ axis --- must never exceed the jet $p_\perp$. This gives
	\todo{This bound is too high}
	\begin{equation}
	\label{eq:pt_max}
	p_{i\perp,\text{max}} = \frac{p_{{\cal J}\perp} - \sum_{j<i} p_{j\perp}
	\cos \Delta
	\phi_j}{\cos \Delta
	\phi_i}\,.
	\end{equation}
	The $x$ and $y$ components of the last parton follow immediately from
	the first $\delta$ function. The last rapidity is fixed by the condition that
	the jet rapidity is kept fixed by the reshuffling, i.e.
	\begin{equation}
	\label{eq:yJ_delta}
	y^B_{\cal J} = y_{\cal J} = \frac 1 2 \ln \frac{\sum_{i=1}^n E_i+ p_{iz}}{\sum_{i=1}^n E_i - p_{iz}}\,.
	\end{equation}
	With $E_n \pm p_{nz} = p_{n\perp}\exp(\pm y_n)$ this can be rewritten to
	\begin{equation}
	\label{eq:yn_quad_eq}
	\exp(2y_{\cal J}) = \frac{\sum_{i=1}^{n-1} E_i+ p_{iz}+p_{n\perp} \exp(y_n)}{\sum_{i=1}^{n-1} E_i - p_{iz}+p_{n\perp} \exp(-y_n)}\,,
	\end{equation}
	which is a quadratic equation in $\exp(y_n)$. The physical solution is
	\begin{align}
	\label{eq:yn}
	y_n ={}& \log\Big(-b + \sqrt{b^2 + \exp(2y_{\cal J})}\,\Big)\,,\\
	b ={}& \bigg(\sum_{i=1}^{n-1} E_i + p_{iz} - \exp(2y_{\cal J})
	\sum_{i=1}^{n-1} E_i - p_{iz}\bigg)/(2 p_{n\perp})\,.
	\end{align}
	\todo{what's wrong with the following?} To eliminate the remaining rapidity
	integral, we transform the $\delta$ function to be linear in the
	rapidity $y$ of the last parton. The corresponding Jacobian is
	\begin{equation}
	\label{eq:jacobian_y}
	\bigg\|\frac{\partial y_{\cal J}}{\partial y_n}\bigg\|^{-1} = 2 \bigg( \frac{E_n +
	p_{nz}}{E_{\cal J} + p_{{\cal J}z}} + \frac{E_n - p_{nz}}{E_{\cal J} -
	p_{{\cal J}z}}\bigg)^{-1}\,.
	\end{equation}
	Finally, we check that all designated constituents are actually
	clustered into the considered jet.

	\subsubsection{Final steps}
	\label{sec:final}

	Knowing the rapidity span covered by the extremal partons, we can now
	generate the rapdities for the partons outside jets. We perform jet
	clustering on all partons and check in
	\lstinline!PhaseSpacePoint::jets_ok! that all the following criteria are
	fulfilled:
	\begin{itemize}
	\item The number of resummation jets must match the number of
	fixed-order jets.
	\item No partons designated to be outside jets may end up inside jets.
	\item All other outgoing partons \emph{must} end up inside jets.
	\item The extremal (in rapidity) partons must be inside the extremal
	jets. If there is, for example, an unordered forward emission, the
	most forward parton must end up inside the most forward jet and the
	next parton must end up inside second jet.
	\item The rapidities of fixed-order and resummation jets must match.
	\end{itemize}
	After this, we adjust the phase-space normalisation according to the
	third line of eq.~(\ref{eq:resumdijetFKLmatched2}), determine the
	flavours of the outgoing partons, and adopt any additional colourless
	bosons from the fixed-order input event. Finally, we use momentum
	conservation to reconstruct the momenta of the incoming partons.

	\subsection{The matrix element }
	\label{sec:ME}

	The derivation of the \HEJ matrix element is explained in some detail
	in~\cite{Andersen:2017kfc}, where also results for leading and
	subleading matrix elements for pure multijet production and production
	of a Higgs boson with at least two associated jets are listed. Matrix
	elements for $W$ and $Z/\gamma^*$ production together with jets are
	given in~\cite{Andersen:2012gk,Andersen:2016vkp}, but not yet included.

	The matrix elements are implemented in the \lstinline!MatrixElement!
	class. To discuss the structure, let us consider the squared matrix
	element for FKL multijet production with $n$ final-state partons:
	\begin{align}
	\label{eq:ME}
	\begin{split}
	\overline{\left\|\mathcal{M}_\text{HEJ}^{f_1 f_2 \to f_1
	g\cdots g f_2}\right\|}^2 = \ &\frac {(4\pi\alpha_s)^n} {4\ (N_c^2-1)}
	\cdot\ \textcolor{blue}{\frac {K_{f_1}(p_1^-, p_a^-)} {t_1}\ \cdot\ \frac{K_{f_2}(p_n^+, p_b^+)}{t_{n-1}}\ \cdot\ \left\\|S_{f_1 f_2\to f_1 f_2}\right\\|^2}\\
	& \cdot \prod_{i=1}^{n-2} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
	V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
	& \cdot \prod_{j=1}^{n-1} \textcolor{red}{\exp\left[\omega^0(q_{j\perp})(y_{j+1}-y_j)\right]}.
	\end{split}
	\end{align}
	The structure and momentum assignment of the unsquared matrix element is
	as illustrated here:
	\begin{center}
	\includegraphics{HEJ_amplitude}
	\end{center}
	The square
	of the complete matrix element as given in eq.~(\ref{eq:ME}) is
	calculated by \lstinline!MatrixElement::operator()!. The \textcolor{red}{last line} of
	eq.~\eqref{eq:ME} constitutes the all-order virtual correction,
	implemented in
	\lstinline!MatrixElement::virtual_corrections!.
	$\omega^0$ is the
	\textit{regularised Regge trajectory}
	\begin{equation}
	\label{eq:omega_0}
	\omega^0(q_\perp) = - C_A \frac{\alpha_s}{\pi} \log \left(\frac{q_\perp^2}{\lambda^2}\right)\,,
	\end{equation}
	where $\lambda$ is the slicing parameter limiting the softness of real
	gluon emissions, cf. eq.~(\ref{eq:resumdijetFKLmatched2}).

	The remaining parts, which correspond to the square of the leading-order
	HEJ matrix element $\overline{\left\|\mathcal{M}_\text{LO,
	HEJ}^{f_1f_2\to f_1g\cdots
	gf_2}\big(\big\{p^B_j\big\}\big)\right\|}^{2}$, are computed in
	\lstinline!MatrixElement::tree!. We can further factor off the
	scale-dependent ``parametric'' part
	\lstinline!MatrixElement::tree_param! containing all factors of the
	strong coupling $4\pi\alpha_s$. Using this function saves some CPU time
	when adjusting the renormalisation scale, see
	section~\ref{sec:resum}. The remaining ``kinematic'' factors are
	calculated in \lstinline!MatrixElement::kin!.

	\subsubsection{Matrix elements for Higgs plus dijet}
	\label{sec:ME_h_jets}

	In the production of a Higgs boson together with jets the parametric
	parts and the virtual corrections only require minor changes in the
	respective functions. However, in the ``kinematic'' parts we have to
	distinguish between several cases, which is done in
	\lstinline!MatrixElement::tree_kin_Higgs!. The Higgs boson can be
	\emph{central}, i.e. inside the rapidity range spanned by the extremal
	partons (\lstinline!MatrixElement::tree_kin_Higgs_central!) or
	\emph{peripheral} and outside this range
	(\lstinline!MatrixElement::tree_kin_Higgs_first! or
	\lstinline!MatrixElement::tree_kin_Higgs_last!). In case there is also
	an unordered gluon emission the matrix element is already suppressed by one
	logarithm $\log(s/t)$. Therefore at NLL the Higgs boson can only be emitted
	centrally\footnote{In principle emitting a Higgs boson \textit{on the other
	side} of the unordered gluon is possible by contracting an unordered and
	external Higgs current. Obviously this would not cover all possible
	configurations, e.g. $qQ\to HgqQ$ requires contraction of the standard $Q\to Q$
	current with an (unknown) $q\to Hgq$ one.}.

	If a Higgs boson with momentum $p_H$ is emitted centrally, after parton
	$j$ in rapidity, the matrix element reads
	\begin{equation}
	\label{eq:ME_h_jets_central}
	\begin{split}
	\overline{\left\|\mathcal{M}_\text{HEJ}^{f_1 f_2 \to f_1 g\cdot H
	\cdot g f_2}\right\|}^2 = \ &\frac {\alpha_s^2 (4\pi\alpha_s)^n} {4\ (N_c^2-1)}
	\cdot\ \textcolor{blue}{\frac {K_{f_1}(p_1^-, p_a^-)} {t_1}\
	\cdot\ \frac{1}{t_j t_{j+1}} \cdot\ \frac{K_{f_2}(p_n^+, p_b^+)}{t_{n}}\ \cdot\ \left\\|S_{f_1
	f_2\to f_1 H f_2}\right\\|^2}\\
	& \cdot \prod_{\substack{i=1\\i \neq j}}^{n-1} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
	V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
	& \cdot \textcolor{red}{\prod_{i=1}^{n-1}
	\exp\left[\omega^0(q_{i\perp})\Delta y_i\right]}
	\end{split}
	\end{equation}
	with the momentum definitions
	\begin{center}
	\includegraphics{HEJ_central_Higgs_amplitude}
	\end{center}
	$q_i$ is the $i$th $t$-channel momentum and $\Delta y_i$ the rapidity
	gap between outgoing \emph{particles} (not partons) $i$ and $i+1$ in
	rapidity ordering.

	For \emph{peripheral} emission in the backward direction
	(\lstinline!MatrixElement::tree_kin_Higgs_first!) we first check whether
	the most backward parton is a gluon or an (anti-)quark. In the latter
	case the leading contribution to the matrix element arises through
	emission off the $t$-channel gluons and we can use the same formula
	eq.~(\ref{eq:ME_h_jets_central}) as for central emission. If the most
	backward parton is a gluon, the square of the matrix element can be
	written as
	\begin{equation}
	\label{eq:ME_h_jets_peripheral}
	\begin{split}
	\overline{\left\|\mathcal{M}_\text{HEJ}^{g f_2 \to H g\cdot g f_2}\right\|}^2 = \ &\frac {\alpha_s^2 (4\pi\alpha_s)^n} {\textcolor{blue}{4\ (N_c^2-1)}}
	\textcolor{blue}{\cdot\ K_{H}\
	\frac{K_{f_2}(p_n^+, p_b^+)}{t_{n-1}}\ \cdot\ \left\\|S_{g
	f_2\to H g f_2}\right\\|^2}\\
	& \cdot \prod_{\substack{i=1}}^{n-2} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
	V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
	& \cdot \textcolor{red}{\prod_{i=1}^{n-1}
	\exp\left[\omega^0(q_{i\perp}) (y_{i+1} - y_i)\right]}
	\end{split}
	\end{equation}
	with the momenta as follows:
	\begin{center}
	\includegraphics{HEJ_peripheral_Higgs_amplitude}
	\end{center}
	The \textcolor{blue}{blue part} is implemented in
	\lstinline!MatrixElement::MH2_forwardH!. All other building blocks are
	already available.\todo{Impact factors} The actual current contraction
	is calculated in \lstinline!MH2gq_outsideH! inside
	\lstinline!currents.cc!, which corresponds to $\tfrac{16 \pi^2}{t_1} \left\\|S_{g
	f_2\to H g f_2}\right\\|^2$.\todo{Fix this insane normalisation}

	The forward emission of a Higgs boson is completely analogous. We can
	use the same function \lstinline!MatrixElement::MH2_forwardH!, swapping
	$p_1 \leftrightarrow p_n,\,p_a \leftrightarrow p_b$.

	\subsubsection{FKL ladder and Lipatov vertices}
	\label{sec:FKL_ladder}

	The ``FKL ladder'' is the product
	\begin{equation}
	\label{eq:FKL_ladder}
	\prod_{i=1}^{n-2} \left( \frac{-C_A}{t_it_{i+1}}\
	V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)
	\end{equation}
	appearing in the square of the matrix element for $n$ parton production,
	cf. eq.~(\ref{eq:ME}), and implemented in
	\lstinline!MatrixElement::FKL_ladder_weight!. The Lipatov vertex contraction
	$V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1})$ is implemented \lstinline!C2Lipatovots!.
	It is given by \todo{equation} \todo{mention difference between the two versions
	of \lstinline!C2Lipatovots!, maybe even get rid of one}.

	\subsubsection{Currents}
	\label{sec:currents}

	The current factors $\frac{K_{f_1}K_{f_2}}{t_1 t_{n-1}}\left\\|S_{f_1
	f_2\to f_1 f_2}\right\\|^2$ and their extensions for unordered and Higgs
	boson emissions are implemented in the \lstinline!jM2!$\dots$ functions
	of \texttt{src/currents.cc}. \todo{Only $\\|S\\|^2$ should be in currents}
	\footnote{The current implementation for
	Higgs production in \texttt{src/currents.cc} includes the $1/4$ factor
	inside $S$, opposing to~\eqref{eq:ME}. Thus the overall normalisation is
	unaffected.} The ``colour acceleration multiplier'' (CAM) $K_{f}$
	for a parton $f\in\{g,q,\bar{q}\}$ is defined as
	\begin{align}
	\label{eq:K_g}
	K_g(p_1^-, p_a^-) ={}& \frac{1}{2}\left(\frac{p_1^-}{p_a^-} + \frac{p_a^-}{p_1^-}\right)\left(C_A -
	\frac{1}{C_A}\right)+\frac{1}{C_A}\\
	\label{eq:K_q}
	K_q(p_1^-, p_a^-) ={}&K_{\bar{q}}(p_1^-, p_a^-) = C_F\,.
	\end{align}
	The Higgs current CAM used in eq.~(\ref{eq:ME_h_jets_peripheral}) is
	\begin{equation}
	\label{eq:K_H}
	K_H = C_A\,.
	\end{equation}
	The current contractions are given by\todo{check all this
	carefully!}
	\begin{align}
	\label{eq:S}
	\left\\|S_{f_1 f_2\to f_1 f_2}\right\\|^2 ={}& \sum_{\substack{\lambda_a =
	+,-\\\lambda_b = +,-}} \left\|j^{\lambda_a}_\mu(p_1, p_a)\
	j^{\lambda_b\,\mu}(p_n, p_b)\right\|^2 = 2\sum_{\lambda =
	+,-} \left\|j^{-}_\mu(p_1, p_a)\ j^{\lambda\,\mu}(p_n, p_b)\right\|^2\,,\\
	\left\\|S_{f_1 f_2\to f_1 H f_2}\right\\|^2 ={}& \sum_{\substack{\lambda_a =
	+,-\\\lambda_b = +,-}} \left\|j^{\lambda_a}_\mu(p_1, p_a)V_H^{\mu\nu}(q_j, q_{j+1})\
	j^{\lambda_b}_\nu(p_n, p_b)\right\|^2\,,\\
	\left\\|S_{g f_2 \to H g f_2}\right\\|^2 ={}& \sum_{
	\substack{
	\lambda_{a} = +,-\\
	\lambda_{1} =+,-\\
	\lambda_{b} = +,-
	}}
	\left\|j^{\lambda_a\lambda_1}_{H\,\mu}(p_1, p_a, p_H)\ j^{\lambda_b\,\mu}(p_n, p_b)\right\|^2\,.
	\end{align}
	The ``basic'' currents $j$ are independent of the parton flavour and read
	\begin{equation}
	\label{eq:j}
	j^\pm_\mu(p, q) = u^{\pm,\dagger}(p)\ \sigma^\pm_\mu\ u^{\pm}(q)\,,
	\end{equation}
	where $\sigma_\mu^\pm = (1, \pm \sigma_i)$ and $\sigma_i$ are the Pauli
	matrices
	\begin{equation}
	\label{eq:Pauli_matrices}
	\sigma_1 =
	\begin{pmatrix}
	0 & 1\\ 1 & 0
	\end{pmatrix}
	\,,
	\qquad \sigma_2 =
	\begin{pmatrix}
	0 & -i\\ i & 0
	\end{pmatrix}
	\,,
	\qquad \sigma_3 =
	\begin{pmatrix}
	1 & 0\\ 0 & -1
	\end{pmatrix}
	\,.
	\end{equation}
	The two-component chiral spinors are given by
	\begin{align}
	\label{eq:u_plus}
	u^+(p)={}& \left(\sqrt{p^+}, \sqrt{p^-} \hat{p}_\perp \right) \,,\\
	\label{eq:u_minus}
	u^-(p)={}& \left(\sqrt{p^-} \hat{p}^*_\perp, -\sqrt{p^+}\right)\,,
	\end{align}
	with $p^\pm = E\pm p_z,\, \hat{p}_\perp = \tfrac{p_\perp}{\|p_\perp\|},\,
	p_\perp = p_x + i p_y$. The spinors for vanishing transverse momentum
	are obtained by replacing $\hat{p_\perp} \to -1$.

	Explicitly, the currents read
	\begin{align}
	\label{eq:j-_explicit}
	j^-_\mu(p, q) ={}&
	\begin{pmatrix}
	\sqrt{p^+\,q^+} + \sqrt{p^-\,q^-} \hat{p}_{\perp} \hat{q}_{\perp}^*\\
	\sqrt{p^-\,q^+}\, \hat{p}_{\perp} + \sqrt{p^+\,q^-}\,\hat{q}_{\perp}^*\\
	-i \sqrt{p^-\,q^+}\, \hat{p}_{\perp} + i \sqrt{p^+\,q^-}\, \hat{q}_{\perp}^*\\
	\sqrt{p^+\,q^+} - \sqrt{p^-\,q^-}\, \hat{p}_{\perp}\, \hat{q}_{\perp}^*
	\end{pmatrix}\\
	j^+_\mu(p, q) ={}&\big(j^-_\mu(p, q)\big)^*
	\end{align}
	If $q= p_{\text{in}}$ is the momentum of an incoming parton, we have
	$\hat{p}_{\text{in} \perp} = -1$ and either $p_{\text{in}}^+ = 0$ or
	$p_{\text{in}}^- = 0$. The current simplifies further:\todo{Helicities flipped w.r.t code}
	\begin{align}
	\label{eq:j_explicit}
	j^-_\mu(p_{\text{out}}, p_{\text{in}}) ={}&
	\begin{pmatrix}
	\sqrt{p_{\text{in}}^+\,p_{\text{out}}^+}\\
	\sqrt{p_{\text{in}}^+\,p_{\text{out}}^-} \ \hat{p}_{\text{out}\,\perp}\\
	-i\,j^-_1\\
	j^-_0
	\end{pmatrix}
	& p_{\text{in}\,z} > 0\,,\\
	j^-_\mu(p_{\text{out}}, p_{\text{in}}) ={}&
	\begin{pmatrix}
	-\sqrt{p_{\text{in}}^-\,p_{\text{out}}^{-\phantom{+}}} \ \hat{p}_{\text{out}\,\perp}\\
	- \sqrt{p_{\text{in}}^-\,p_{\text{out}}^+}\\
	i\,j^-_1\\
	-j^-_0
	\end{pmatrix} & p_{\text{in}\,z} < 0\,.
	\end{align}

	\section{The fixed-order generator}
	\label{sec:HEJFOG}

	Even at leading order, standard fixed-order generators can only generate
	events with a limited number of final-state particles within reasonable
	CPU time. The purpose of the fixed-order generator is to supplement this
	with high-multiplicity input events according to the first two lines of
	eq.~\eqref{eq:resumdijetFKLmatched2} with the \HEJ approximation
	$\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to f_1g\cdots gf_2}$ instead of the
	full fixed-order matrix element $\mathcal{M}_\text{LO}^{f_1f_2\to
	f_1g\cdots gf_2}$. Its usage is described in the user
	documentation \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJFOG.html}.

	\subsection{File structure}
	\label{sec:HEJFOG_structure}

	The code for the fixed-order generator is in the \texttt{FixedOrderGen}
	directory, which contains the following:
	\begin{description}
	\item[include:] Contains the C++ header files.
	\item[src:] Contains the C++ source files.
	\item[t:] Contains the source code for the automated tests.
	\item[CMakeLists.txt:] Configuration file for the \cmake build system.
	\item[configFO.yml:] Sample configuration file for the fixed-order generator.
	\end{description}
	The code is generally in the \lstinline!HEJFOG! namespace. Functions and
	classes \lstinline!MyClass! are usually declared in
	\texttt{include/MyClass.hh} and implemented in \texttt{src/MyClass.cc}.

	\subsection{Program flow}
	\label{sec:prog_flow}

	A single run of the fixed-order generator consists of three or four
	stages.

	First, we perform initialisation similar to HEJ 2, see
	section~\ref{sec:init}. Since there is a lot of overlap we frequently
	reuse classes and functions from HEJ 2, i.e. from the
	\lstinline!HEJ! namespace. The code for parsing the configuration file
	is in \texttt{include/config.hh} and implemented in
	\texttt{src/config.cc}.

	If partial unweighting is requested in the user settings \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJFOG.html#settings},
	the initialisation is followed by a calibration phase. We use a
	\lstinline!EventGenerator! to produce a number of trial
	events. We use these to calibrate the \lstinline!Unweighter! in
	its constructor and produce a first batch of partially unweighted
	events. This also allows us to estimate our unweighting efficiency.

	In the next step, we continue to generate events and potentially
	unweight them. Once the user-defined target number of events is reached,
	we adjust their weights according to the number of required trials. As
	in HEJ 2 (see section~\ref{sec:processing}), we pass the final
	events to a \lstinline!HEJ::Analysis! and a
	\lstinline!HEJ::CombinedEventWriter!.

	\subsection{Event generation}
	\label{sec:evgen}

	Event generation is performed by the
	\lstinline!EventGenerator::gen_event! member function. We begin by generating a
	\lstinline!PhaseSpacePoint!. This is not to be confused with
	the resummation phase space points represented by
	\lstinline!HEJ::PhaseSpacePoint!! After jet clustering, we compute the
	leading-order matrix element (see section~\ref{sec:ME}) and pdf factors.

	The phase space point generation is performed in the
	\lstinline!PhaseSpacePoint! constructor. We first construct the
	user-defined number of $n_p$ partons (by default gluons) in
	\lstinline!PhaseSpacePoint::gen_LO_partons!. We use flat sampling in
	rapidity and azimuthal angle. For the scalar transverse momenta, we
	distinguish between two cases. By default, they are generated based on a
	random variable $x_{p_\perp}$ according to
	\begin{equation}
	\label{eq:pt_sampling}
	p_\perp = p_{\perp,\text{min}} +
	\begin{cases}
	p_{\perp,\text{par}}
	\tan\left(
	x_{p_\perp}
	\arctan\left(
	\frac{p_{\perp,\text{max}} - p_{\perp,\text{min}}}{p_{\perp,\text{par}}}
	\right)
	\right)
	& y < y_\text{cut}
	\\
	- \tilde{p}_{\perp,\text{par}}\log\left(1 - x_{p_\perp}\left[1 -
	\exp\left(\frac{p_{\perp,\text{min}} -
	p_{\perp,\text{max}}}{\tilde{p}_{\perp,\text{par}}}\right)\right]\right)
	& y \geq y_\text{cut}
	\end{cases}\,,
	\end{equation}
	where $p_{\perp,\text{min}}$ is the minimum jet transverse momentum,
	$p_{\perp,\text{max}}$ is the maximum transverse parton momentum,
	tentatively set to the beam energy, and $y_\text{cut}$, $p_{\perp,\text{par}}$
	and $\tilde{p}_{\perp,\text{par}}$ are generation parameters set to
	heuristically determined values of
	\begin{align}
	y_\text{cut}&=3,\\
	p_{\perp,\text{par}}&=p_{\perp,\min}+\frac{n_p}{5}, \\
	\tilde{p}_{\perp,\text{par}}&=\frac{p_{\perp,\text{par}}}{1 +
	5(y-y_\text{cut})}.
	\end{align}
	The problem with this generation is that the transverse momenta peak at
	the minimum transverse momentum required for fixed-order jets. However,
	if we use the generated events as input for \HEJ resummation, events
	with such soft transverse momenta hardly contribute, see
	section~\ref{sec:ptj_res}. To generate efficient input for resummation,
	there is the user option \texttt{peak pt}, which specifies the
	dominant transverse momentum for resummation jets. If this option is
	set, most jets will be generated as above, but with
	$p_{\perp,\text{min}}$ set to the peak transverse momentum $p_{\perp,
	\text{peak}}$. In addition, there is a small chance of around $2\%$ to
	generate softer jets. The heuristic ansatz for the transverse momentum
	distribution in the ``soft'' region is
	\begin{equation}
	\label{FO_pt_soft}
	\frac{\partial \sigma}{\partial p_\perp} \propto e^{n_p\frac{p_\perp- p_{\perp,
	\text{peak}}}{\bar{p}_\perp}}\,,
	\end{equation}
	where $n_p$ is the number of partons and $\bar{p}_\perp \approx
	4\,$GeV. To achieve this distribution, we use
	\begin{equation}
	\label{eq:FO_pt_soft_sampling}
	p_\perp = p_{\perp, \text{peak}} + \bar{p}_\perp \frac{\log x_{p_\perp}}{n_p}
	\end{equation}
	and discard the phase space point if the parton is too soft, i.e. below the threshold for
	fixed-order jets.

	After ensuring that all partons form separate jets, we generate any
	potential colourless emissions. We then determine the incoming momenta
	and flavours in \lstinline!PhaseSpacePoint::reconstruct_incoming! and
	adjust the outgoing flavours to ensure an FKL configuration. Finally, we
	may reassign outgoing flavours to generate suppressed (for example
	unordered) configurations.

	\subsection{Unweighting}
	\label{sec:unweight}

	Straightforward event generation tends to produce many events with small
	weights. Those events have a negligible contribution to the final
	observables, but can take up considerable storage space and CPU time in
	later processing stages. This problem can be addressed by unweighting.

	For naive unweighting, one would determine the maximum weight
	$w_\text{max}$ of all events, discard each event with weight $w$ with a
	probability $p=w/w_\text{max}$, and set the weights of all remaining
	events to $w_\text{max}$. The downside to this procedure is that it also
	eliminates a sizeable fraction of events with moderate weight, so that
	the statistical convergence deteriorates.

	To ameliorate this problem, we perform unweighting only for events with
	sufficiently small weights. This is done by the
	\lstinline!Unweighter! class. In the constructor we estimate the
	mean and width of the weight-weight distribution from a sample of
	events. We use these estimates to determine the maximum weight below
	which unweighting is performed. The actual unweighting is the done in
	the \lstinline!Unweighter::unweight! function.

	\input{currents}

	+\appendix
	+\section{Continuous Integration}
	+\label{sec:gitlabCI}
	+
	+GitLab provides ways to directly test code via \textit{Continuous integrations}.
	+The CI is controlled by \texttt{.gitlab-ci.yml}. For all options for the YAML
	+file see \href{https://docs.gitlab.com/ee/ci/yaml/}{docs.gitlab.com/ee/ci/yaml/}.
	+GitLab also provides a small tool to check that YAML syntax is correct under
	+\lstinline!CI/CD > Pipelines > CI Lint! or
	+\href{https://gitlab.dur.scotgrid.ac.uk/hej/HEJ/-/ci/lint}{gitlab.dur.scotgrid.ac.uk/hej/HEJ/-/ci/lint}.
	+
	+Currently the CI is configured to trigger a \textit{Pipeline} on each
	+\lstinline!git push!. The corresponding \textit{GitLab runners} are configured
	+under \lstinline!CI/CD Settings>Runners! in the GitLab UI. All runners use a
	+\href{https://www.docker.com/}{docker} image as virtual environments\footnote{To
	+use only Docker runners set the \lstinline!docker! tag in
	+\texttt{.gitlab-ci.yml}.}. The specific docker images maintained separately. If
	+you add a new dependences, please also provide a docker image for the CI. The
	+goal to be able to test \HEJ with all possible configurations.
	+
	+Each pipeline contains multiple stages (see \lstinline!stages! in
	+\texttt{.gitlab-ci.yml}) which are executed in order from top to bottom.
	+Additionally each stage contains multiple jobs. For example the stage
	+\textit{build} contains the jobs \lstinline!build:basic!,
	+\lstinline!build:qcdloop!, \lstinline!build:rivet!, etc., which compile \HEJ for
	+different environments and dependences, by using different in the Docker images.
	+Jobs staring with an dot are ignored by the Runner, e.g. \lstinline!.HEJ_build!
	+is only used as a template but never executed directly. Only after all jobs of
	+the previous stage was executed without any error the next stage will start.
	+
	+To pass information between one stage and the next we use \lstinline!artifacts!.
	+The runner will automatically load all artifacts form all
	+\lstinline!dependencies! for each job\footnote{If no dependencies are defined
	+\textit{all} artifacts from all previous jobs are downloaded. Thus please
	+specify an empty dependence if you do not want to load any artifacts.}. For
	+example the compiled \HEJ code from \lstinline!build:basic! gets loaded in
	+\lstinline!test:basic! and \lstinline!FOG:build:basic!, without recompiling \HEJ
	+again. Additionally artifacts can be downloaded from the GitLab web page, which
	+could be handy for debugging.
	+
	+The actual commands are given in the \lstinline!before_script!,
	+\lstinline!script! and \lstinline!after_script!
	+\footnote{\lstinline!after_script! is always executed} sections, and are
	+standard Linux shell commands (dependent on the docker image). Any failed
	+command, i.e. returning not zero, stops the job and making the pipeline fail
	+entirely. Most tests are just running \lstinline!make test! or are based on it.
	+Thus, to emphasise it again, write tests for your code in \lstinline!cmake!. The
	+CI is only intended to make automated testing in different environments easier.
	+
	+
	\bibliographystyle{JHEP}
	\bibliography{biblio}


	\end{document}
	diff --git a/include/HEJ/CrossSectionAccumulator.hh b/include/HEJ/CrossSectionAccumulator.hh
	new file mode 100644
	index 0000000..9be4d9f
	--- /dev/null
	+++ b/include/HEJ/CrossSectionAccumulator.hh
	@@ -0,0 +1,69 @@
	+#pragma once
	+
	+#include <map>
	+#include <ostream>
	+#include <iomanip>
	+#include <string>
	+
	+#include "HEJ/Event.hh"
	+#include "HEJ/event_types.hh"
	+
	+namespace HEJ {
	+ template<typename T>
	+ struct XSWithError {
	+ T value = T{};
	+ T error = T{};
	+ };
	+
	+ /**
	+ * @brief Sum of Cross Section for different subproccess
	+ */
	+ class CrossSectionAccumulator {
	+ public:
	+ void fill(HEJ::Event const & ev) {
	+ const double wt = ev.central().weight;
	+ auto & entry = xs_[ev.type()];
	+ entry.value += wt;
	+ entry.error += wt*wt;
	+ total_.value += wt;
	+ total_.error += wt*wt;
	+ }
	+
	+ auto begin() const {
	+ return std::begin(xs_);
	+ }
	+
	+ auto end() const {
	+ return std::end(xs_);
	+ }
	+
	+ //! total Cross Section and error
	+ XSWithError<double> total() const {
	+ return total_;
	+ }
	+
	+
	+ private:
	+ std::map<HEJ::event_type::EventType, XSWithError<double>> xs_;
	+ XSWithError<double> total_;
	+ };
	+
	+ std::ostream& operator<<(std::ostream& os, const CrossSectionAccumulator& xs){
	+ const std::streamsize orig_prec = os.precision();
	+ os << std::scientific << std::setprecision(10)
	+ << " " << std::left << std::setw(20)
	+ << "Cross section: " << xs.total().value
	+ << " +- " << std::sqrt(xs.total().error) << " (pb)\n";
	+ for(auto const & xs_type: xs) {
	+ os
	+ << " " << std::left << std::setw(20)
	+ << (HEJ::event_type::names[xs_type.first] + std::string(": "))
	+ << xs_type.second.value << " +- "
	+ << std::sqrt(xs_type.second.error) << " (pb)\n";
	+ }
	+ os << std::defaultfloat;
	+ os.precision(orig_prec);
	+ return os;
	+ }
	+
	+}
	diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
	index ca81ec7..f0c672a 100644
	--- a/src/bin/HEJ.cc
	+++ b/src/bin/HEJ.cc
	@@ -1,185 +1,191 @@
	/**
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <chrono>
	#include <iostream>
	#include <limits>
	#include <memory>
	#include <numeric>

	#include "yaml-cpp/yaml.h"

	#include "LHEF/LHEF.h"

	#include "fastjet/ClusterSequence.hh"

	#include "HEJ/CombinedEventWriter.hh"
	#include "HEJ/config.hh"
	+#include "HEJ/CrossSectionAccumulator.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/get_analysis.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/ProgressBar.hh"
	#include "HEJ/stream.hh"
	#include "HEJ/Version.hh"
	#include "HEJ/YAMLreader.hh"

	int event_number(std::string const & record){
	size_t start = record.rfind("Number of Events");
	start = record.find_first_of("123456789", start);
	if(start == std::string::npos) {
	throw std::invalid_argument("no event number record found");
	}
	const size_t end = record.find_first_not_of("0123456789", start);
	return std::stoi(record.substr(start, end - start));
	}

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

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

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

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

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

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


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

	fastjet::ClusterSequence::print_banner();

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

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

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

	LHEF::Reader reader{in};
	auto heprup = reader.heprup;
	heprup.generators.emplace_back(LHEF::XMLTag{});
	heprup.generators.back().name = HEJ::Version::package_name();
	heprup.generators.back().version = HEJ::Version::String();
	HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};

	double global_reweight = 1.;
	int max_events = std::numeric_limits<int>::max();
	if(argn > 3){
	max_events = std::stoi(argv[3]);
	const int input_events = event_number(reader.headerBlock);
	global_reweight = input_events/static_cast<double>(max_events);
	std::cout << "Processing " << max_events
	<< " out of " << input_events << " events\n";
	}
	HEJ::ScaleGenerator scale_gen{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	};
	auto ran = HEJ::make_RNG(config.rng.name, config.rng.seed);
	assert(ran != nullptr);
	HEJ::EventReweighter hej{
	reader.heprup,
	std::move(scale_gen),
	to_EventReweighterConfig(config),
	*ran
	};

	int nevent = 0;
	std::array<int, HEJ::event_type::last_type + 1>
	nevent_type{0}, nfailed_type{0};
	auto progress = make_progress_bar(reader.heprup.XSECUP);
	+ HEJ::CrossSectionAccumulator xs;
	// Loop over the events in the inputfile
	while(reader.readEvent()){
	// reweight events so that the total cross section is conserved
	reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());

	if(nevent == max_events) break;
	++nevent;

	// calculate HEJ weight
	HEJ::Event FO_event{
	HEJ::UnclusteredEvent{reader.hepeup},
	config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
	};
	auto resummed_events = hej.reweight(FO_event, config.trials);
	++nevent_type[FO_event.type()];

	if(resummed_events.empty()) ++nfailed_type[FO_event.type()];

	for(auto const & ev: resummed_events){
	//TODO: move pass_cuts to after phase space point generation
	if(analysis->pass_cuts(ev, FO_event)){
	analysis->fill(ev, FO_event);
	writer.write(ev);
	+ xs.fill(ev);
	}
	}
	if(progress) progress->increment(FO_event.central().weight);
	} // main event loop
	std::cout << '\n';
	analysis->finalise();

	using namespace HEJ::event_type;
	std::cout<< "Events processed: " << nevent << '\n';
	for(size_t ev_type = first_type; ev_type <= last_type; ++ev_type){
	std::cout << '\t' << names[ev_type] << ": " << nevent_type[ev_type]
	<< ", failed to reconstruct " << nfailed_type[ev_type]
	<< '\n';
	}
	+
	+ std::cout << '\n' << xs << '\n';
	+
	std::chrono::duration<double> run_time = (clock::now() - start_time);
	std::cout << "Finished " << HEJ::Version::package_name() << " at "
	<< time_to_string(clock::to_time_t(clock::now()))
	<< "\n=> Runtime: " << run_time.count() << " sec ("
	<< nevent/run_time.count() << " Events/sec).\n";

	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 6:54 PM (1 d, 14 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805685
Default Alt Text: (97 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions