No OneTemporary
Actions

Size

45 KB

Subscribers

None

View Options

	diff --git a/Changes.md b/Changes.md
	index aea4b72..b30bb04 100644
	--- a/Changes.md
	+++ b/Changes.md
	@@ -1,33 +1,34 @@
	# 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
	+* Removed default constructor for `Event`

	## 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/EventGenerator.hh b/FixedOrderGen/include/EventGenerator.hh
	index a039144..e2064c7 100644
	--- a/FixedOrderGen/include/EventGenerator.hh
	+++ b/FixedOrderGen/include/EventGenerator.hh
	@@ -1,56 +1,57 @@
	#pragma once

	-#include "HEJ/PDF.hh"
	#include "HEJ/MatrixElement.hh"
	+#include "HEJ/optional.hh"
	+#include "HEJ/PDF.hh"
	#include "HEJ/RNG.hh"

	+#include "Beam.hh"
	#include "JetParameters.hh"
	+#include "ParticleProperties.hh"
	#include "Process.hh"
	-#include "Beam.hh"
	#include "Status.hh"
	-#include "ParticleProperties.hh"

	namespace HEJ{
	class Event;
	class HiggsCouplingSettings;
	class ScaleGenerator;
	}

	//! Namespace for HEJ Fixed Order Generator
	namespace HEJFOG{
	class EventGenerator{
	public:
	EventGenerator(
	Process process,
	Beam beam,
	HEJ::ScaleGenerator scale_gen,
	JetParameters jets,
	int pdf_id,
	double subl_change,
	unsigned int subl_channels,
	ParticlesPropMap particles_properties,
	HEJ::HiggsCouplingSettings Higgs_coupling,
	HEJ::RNG & ran
	);

	- HEJ::Event gen_event();
	+ HEJ::optional<HEJ::Event> gen_event();

	Status status() const {
	return status_;
	}

	private:
	HEJ::PDF pdf_;
	HEJ::MatrixElement ME_;
	HEJ::ScaleGenerator scale_gen_;
	Process process_;
	JetParameters jets_;
	Beam beam_;
	Status status_;
	double subl_change_;
	unsigned int subl_channels_;
	ParticlesPropMap particles_properties_;
	std::reference_wrapper<HEJ::RNG> ran_;
	};

	}
	diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
	index 3a27b74..be6b61e 100644
	--- a/FixedOrderGen/src/EventGenerator.cc
	+++ b/FixedOrderGen/src/EventGenerator.cc
	@@ -1,79 +1,79 @@
	#include "EventGenerator.hh"

	#include "Process.hh"
	#include "Beam.hh"
	#include "JetParameters.hh"
	#include "PhaseSpacePoint.hh"

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

	namespace HEJFOG{
	EventGenerator::EventGenerator(
	Process process,
	Beam beam,
	HEJ::ScaleGenerator scale_gen,
	JetParameters jets,
	int pdf_id,
	double subl_change,
	unsigned int subl_channels,
	ParticlesPropMap particles_properties,
	HEJ::HiggsCouplingSettings Higgs_coupling,
	HEJ::RNG & ran
	):
	pdf_{pdf_id, beam.particles[0], beam.particles[1]},
	ME_{
	[this](double mu){ return pdf_.Halphas(mu); },
	HEJ::MatrixElementConfig{
	false,
	std::move(Higgs_coupling)
	}
	},
	scale_gen_{std::move(scale_gen)},
	process_{std::move(process)},
	jets_{std::move(jets)},
	beam_{std::move(beam)},
	subl_change_{subl_change},
	subl_channels_{subl_channels},
	particles_properties_{std::move(particles_properties)},
	ran_{ran}
	{
	}

	- HEJ::Event EventGenerator::gen_event(){
	+ HEJ::optional<HEJ::Event> EventGenerator::gen_event(){
	HEJFOG::PhaseSpacePoint psp{
	process_,
	jets_,
	pdf_, beam_.energy,
	subl_change_, subl_channels_,
	particles_properties_,
	ran_
	};
	status_ = psp.status();
	if(status_ != good) return {};

	HEJ::Event ev = scale_gen_(
	HEJ::Event{
	to_EventData( std::move(psp) ).cluster( jets_.def, jets_.min_pt)
	}
	);

	const double shat = HEJ::shat(ev);
	const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy);
	const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy);

	// evaluate matrix element on this point
	const auto ME_weights = ME_.tree(ev);
	ev.central().weight = ME_weights.central/(shatshat);
	ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
	ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
	for(size_t i = 0; i < ev.variations().size(); ++i){
	auto & var = ev.variations(i);
	var.weight = ME_weights.variations[i]/(shatshat);
	var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type);
	var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type);
	}
	return ev;
	}

	}
	diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
	index 6d55edb..7f4b583 100644
	--- a/FixedOrderGen/src/main.cc
	+++ b/FixedOrderGen/src/main.cc
	@@ -1,233 +1,235 @@
	/**
	* 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 "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));
	+ assert( (generator.status() == HEJFOG::good) == bool(ev) );
	+ 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)) {
	+ assert( (generator.status() == HEJFOG::good) == bool(ev) );
	+ if(generator.status() == HEJFOG::good && analysis->pass_cuts(ev, ev)) {
	if(unweighter) {
	- auto unweighted = unweighter->unweight(std::move(ev));
	+ auto unweighted = unweighter->unweight(std::move(*ev));
	if(! unweighted) continue;
	- ev = std::move(*unweighted);
	+ ev = std::move(unweighted);
	}
	- events.emplace_back(std::move(ev));
	+ events.emplace_back(std::move(*ev));
	++progress;
	}
	}
	std::cout << std::endl;

	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 << xs << '\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/FixedOrderGen/t/2j.cc b/FixedOrderGen/t/2j.cc
	index 80bdff3..ea68d04 100644
	--- a/FixedOrderGen/t/2j.cc
	+++ b/FixedOrderGen/t/2j.cc
	@@ -1,59 +1,60 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif

	#include <algorithm>
	#include <cmath>
	#include <cassert>
	#include <iostream>

	#include "config.hh"
	#include "EventGenerator.hh"
	#include "HEJ/Mixmax.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/MatrixElement.hh"

	using namespace HEJFOG;

	int main(){
	constexpr double invGeV2_to_pb = 389379292.;
	constexpr double xs_ref = 86.42031848*1e6; //calculated with "combined" HEJ svn r3480

	auto config = load_config("config_2j.yml");

	HEJ::Mixmax ran{};
	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	HEJ::ScaleGenerator{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	},
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	ran
	};

	double xs = 0., xs_err = 0.;
	for (int trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	if(generator.status() != good) continue;
	- ev.central().weight *= invGeV2_to_pb;
	- ev.central().weight /= config.events;
	+ assert(ev);
	+ ev->central().weight *= invGeV2_to_pb;
	+ ev->central().weight /= config.events;

	- xs += ev.central().weight;
	- xs_err += ev.central().weight*ev.central().weight;
	+ xs += ev->central().weight;
	+ xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';

	assert(std::abs(xs - xs_ref) < 3*xs_err);
	assert(xs_err < 0.01*xs);
	}
	diff --git a/FixedOrderGen/t/4j.cc b/FixedOrderGen/t/4j.cc
	index 5b4e1b4..a597242 100644
	--- a/FixedOrderGen/t/4j.cc
	+++ b/FixedOrderGen/t/4j.cc
	@@ -1,60 +1,61 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif

	#include <algorithm>
	#include <cmath>
	#include <cassert>
	#include <iostream>

	#include "config.hh"
	#include "EventGenerator.hh"
	#include "HEJ/Mixmax.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/MatrixElement.hh"

	using namespace HEJFOG;

	int main(){
	constexpr double invGeV2_to_pb = 389379292.;
	constexpr double xs_ref = 0.81063619*1e6; //calculated with "combined" HEJ svn r3480

	auto config = load_config("config_2j.yml");
	config.process.njets = 4;

	HEJ::Mixmax ran{};
	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	HEJ::ScaleGenerator{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	},
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	ran
	};

	double xs = 0., xs_err = 0.;
	for (int trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	if(generator.status() != good) continue;
	- ev.central().weight *= invGeV2_to_pb;
	- ev.central().weight /= config.events;
	+ assert(ev);
	+ ev->central().weight *= invGeV2_to_pb;
	+ ev->central().weight /= config.events;

	- xs += ev.central().weight;
	- xs_err += ev.central().weight*ev.central().weight;
	+ xs += ev->central().weight;
	+ xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';

	assert(std::abs(xs - xs_ref) < 3*xs_err);
	assert(xs_err < 0.03*xs);
	}
	diff --git a/FixedOrderGen/t/h_2j.cc b/FixedOrderGen/t/h_2j.cc
	index 3e3135e..4ad315d 100644
	--- a/FixedOrderGen/t/h_2j.cc
	+++ b/FixedOrderGen/t/h_2j.cc
	@@ -1,67 +1,68 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif

	#include <algorithm>
	#include <cmath>
	#include <cassert>
	#include <iostream>

	#include "config.hh"
	#include "EventGenerator.hh"
	#include "HEJ/Ranlux64.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/MatrixElement.hh"

	using namespace HEJFOG;

	int main(){
	constexpr double invGeV2_to_pb = 389379292.;
	constexpr double xs_ref = 2.04928; // +- 0.00377252
	//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20

	auto config = load_config("config_h_2j.yml");

	HEJ::Ranlux64 ran{};
	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	HEJ::ScaleGenerator{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	},
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	ran
	};

	double xs = 0., xs_err = 0.;
	for (int trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	if(generator.status() != good) continue;
	- ev.central().weight *= invGeV2_to_pb;
	- ev.central().weight /= config.events;
	+ assert(ev);
	+ ev->central().weight *= invGeV2_to_pb;
	+ ev->central().weight /= config.events;

	const auto the_Higgs = std::find_if(
	- begin(ev.outgoing()), end(ev.outgoing()),
	+ begin(ev->outgoing()), end(ev->outgoing()),
	[](HEJ::Particle const & p){ return p.type == HEJ::ParticleID::h; }
	);
	- assert(the_Higgs != end(ev.outgoing()));
	+ assert(the_Higgs != end(ev->outgoing()));
	if(std::abs(the_Higgs->rapidity()) > 5.) continue;

	- xs += ev.central().weight;
	- xs_err += ev.central().weight*ev.central().weight;
	+ xs += ev->central().weight;
	+ xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl;

	assert(std::abs(xs - xs_ref) < 3*xs_err);
	assert(xs_err < 0.01*xs);
	}
	diff --git a/FixedOrderGen/t/h_2j_decay.cc b/FixedOrderGen/t/h_2j_decay.cc
	index 8454b7b..e404707 100644
	--- a/FixedOrderGen/t/h_2j_decay.cc
	+++ b/FixedOrderGen/t/h_2j_decay.cc
	@@ -1,86 +1,87 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif

	#include <algorithm>
	#include <cmath>
	#include <cassert>
	#include <iostream>

	#include "config.hh"
	#include "EventGenerator.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/MatrixElement.hh"
	#include "HEJ/Particle.hh"
	#include "HEJ/PDF.hh"
	#include "HEJ/Ranlux64.hh"
	#include "HEJ/utility.hh"

	using namespace HEJFOG;

	bool pass_dR_cut(
	std::vector<fastjet::PseudoJet> const & jets,
	std::vector<HEJ::Particle> const & photons
	){
	constexpr double delta_R_min = 0.7;
	for(auto const & jet: jets){
	for(auto const & photon: photons){
	if(jet.delta_R(photon.p) < delta_R_min) return false;
	}
	}
	return true;
	}

	int main(){
	constexpr double invGeV2_to_pb = 389379292.;
	constexpr double xs_ref = 0.00429198; // +- 1.0488e-05
	//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20

	auto config = load_config("config_h_2j_decay.yml");

	HEJ::Ranlux64 ran{};
	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	HEJ::ScaleGenerator{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	},
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	ran
	};

	double xs = 0., xs_err = 0.;
	for (int trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	if(generator.status() != good) continue;
	- assert(ev.decays().size() == 1);
	- const auto decay = begin(ev.decays());
	- assert(ev.outgoing().size() > decay->first);
	- const auto & the_Higgs = ev.outgoing()[decay->first];
	+ assert(ev);
	+ assert(ev->decays().size() == 1);
	+ const auto decay = begin(ev->decays());
	+ assert(ev->outgoing().size() > decay->first);
	+ const auto & the_Higgs = ev->outgoing()[decay->first];
	assert(the_Higgs.type == HEJ::pid::Higgs);
	assert(decay->second.size() == 2);
	auto const & gamma = decay->second;
	assert(gamma[0].type == HEJ::pid::photon);
	assert(gamma[1].type == HEJ::pid::photon);
	assert(HEJ::nearby_ep(gamma[0].p + gamma[1].p, the_Higgs.p, 1e-6));
	- if(!pass_dR_cut(ev.jets(), gamma)) continue;
	- ev.central().weight *= invGeV2_to_pb;
	- ev.central().weight /= config.events;
	+ if(!pass_dR_cut(ev->jets(), gamma)) continue;
	+ ev->central().weight *= invGeV2_to_pb;
	+ ev->central().weight /= config.events;

	- xs += ev.central().weight;
	- xs_err += ev.central().weight*ev.central().weight;
	+ xs += ev->central().weight;
	+ xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl;

	assert(std::abs(xs - xs_ref) < 3*xs_err);
	assert(xs_err < 0.012*xs);
	}
	diff --git a/FixedOrderGen/t/h_3j.cc b/FixedOrderGen/t/h_3j.cc
	index cad089a..bb2baff 100644
	--- a/FixedOrderGen/t/h_3j.cc
	+++ b/FixedOrderGen/t/h_3j.cc
	@@ -1,68 +1,69 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif

	#include <algorithm>
	#include <cmath>
	#include <cassert>
	#include <iostream>

	#include "config.hh"
	#include "EventGenerator.hh"
	#include "HEJ/Ranlux64.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/MatrixElement.hh"
	#include "HEJ/PDF.hh"

	using namespace HEJFOG;

	int main(){
	constexpr double invGeV2_to_pb = 389379292.;
	constexpr double xs_ref = 1.07807; // +- 0.0071
	//calculated with HEJ revision 93efdc851b02a907a6fcc63956387f9f4c1111c2 +1

	auto config = load_config("config_h_2j.yml");
	config.process.njets = 3;

	HEJ::Ranlux64 ran{};
	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	HEJ::ScaleGenerator{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	},
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	ran
	};

	double xs = 0., xs_err = 0.;
	for (int trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	if(generator.status() != good) continue;
	- ev.central().weight *= invGeV2_to_pb;
	- ev.central().weight /= config.events;
	+ assert(ev);
	+ ev->central().weight *= invGeV2_to_pb;
	+ ev->central().weight /= config.events;

	const auto the_Higgs = std::find_if(
	- begin(ev.outgoing()), end(ev.outgoing()),
	+ begin(ev->outgoing()), end(ev->outgoing()),
	[](HEJ::Particle const & p){ return p.type == HEJ::ParticleID::h; }
	);
	- assert(the_Higgs != end(ev.outgoing()));
	+ assert(the_Higgs != end(ev->outgoing()));
	if(std::abs(the_Higgs->rapidity()) > 5.) continue;

	- xs += ev.central().weight;
	- xs_err += ev.central().weight*ev.central().weight;
	+ xs += ev->central().weight;
	+ xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl;

	assert(std::abs(xs - xs_ref) < 3*xs_err);
	assert(xs_err < 0.02*xs);
	}
	diff --git a/FixedOrderGen/t/h_3j_uno1.cc b/FixedOrderGen/t/h_3j_uno1.cc
	index 4c0f6f3..c100dcb 100644
	--- a/FixedOrderGen/t/h_3j_uno1.cc
	+++ b/FixedOrderGen/t/h_3j_uno1.cc
	@@ -1,72 +1,73 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif

	// check that adding uno emissions doesn't change the FKL cross section

	#include <algorithm>
	#include <cassert>
	#include <cmath>
	#include <iostream>

	#include "config.hh"
	#include "EventGenerator.hh"
	#include "HEJ/Ranlux64.hh"
	#include "Subleading.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/MatrixElement.hh"
	#include "HEJ/PDF.hh"

	using namespace HEJFOG;

	int main(){
	constexpr double invGeV2_to_pb = 389379292.;
	constexpr double xs_ref = 0.0243548; // +- 0.000119862
	//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20

	auto config = load_config("config_h_2j.yml");
	config.process.njets = 3;
	config.process.incoming = {HEJ::pid::u, HEJ::pid::u};
	config.subleading_channels = HEJFOG::Subleading::uno;

	HEJ::Ranlux64 ran{};
	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	HEJ::ScaleGenerator{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	},
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	ran
	};

	double xs = 0., xs_err = 0.;
	int uno_found = 0;
	for (int trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	if(generator.status() != good) continue;
	- if(ev.type() != HEJ::event_type::FKL){
	+ assert(ev);
	+ if(ev->type() != HEJ::event_type::FKL){
	++uno_found;
	continue;
	}
	- ev.central().weight *= invGeV2_to_pb;
	- ev.central().weight /= config.events;
	+ ev->central().weight *= invGeV2_to_pb;
	+ ev->central().weight /= config.events;

	- xs += ev.central().weight;
	- xs_err += ev.central().weight*ev.central().weight;
	+ xs += ev->central().weight;
	+ xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
	std::cout << uno_found << " events with unordered emission" << std::endl;
	assert(uno_found > 0);
	assert(std::abs(xs - xs_ref) < 3*xs_err);
	assert(xs_err < 0.05*xs);
	}
	diff --git a/FixedOrderGen/t/h_3j_uno2.cc b/FixedOrderGen/t/h_3j_uno2.cc
	index 5c6c409..4dc4770 100644
	--- a/FixedOrderGen/t/h_3j_uno2.cc
	+++ b/FixedOrderGen/t/h_3j_uno2.cc
	@@ -1,67 +1,68 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif

	// check uno cross section

	#include <algorithm>
	#include <cassert>
	#include <cmath>
	#include <iostream>

	#include "config.hh"
	#include "EventGenerator.hh"
	#include "HEJ/Ranlux64.hh"
	#include "Subleading.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/MatrixElement.hh"
	#include "HEJ/PDF.hh"

	using namespace HEJFOG;

	int main(){
	constexpr double invGeV2_to_pb = 389379292.;
	constexpr double xs_ref = 0.00347538; // +- 3.85875e-05
	//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20

	auto config = load_config("config_h_2j.yml");
	config.process.njets = 3;
	config.process.incoming = {HEJ::pid::u, HEJ::pid::u};
	config.subleading_fraction = 1.;
	config.subleading_channels = HEJFOG::Subleading::uno;

	HEJ::Ranlux64 ran{};
	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	HEJ::ScaleGenerator{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	},
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	ran
	};

	double xs = 0., xs_err = 0.;
	for (int trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	- if(ev.type() == HEJ::event_type::FKL) continue;
	if(generator.status() != good) continue;
	- ev.central().weight *= invGeV2_to_pb;
	- ev.central().weight /= config.events;
	+ assert(ev);
	+ if(ev->type() == HEJ::event_type::FKL) continue;
	+ ev->central().weight *= invGeV2_to_pb;
	+ ev->central().weight /= config.events;

	- xs += ev.central().weight;
	- xs_err += ev.central().weight*ev.central().weight;
	+ xs += ev->central().weight;
	+ xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl;
	assert(std::abs(xs - xs_ref) < 3*xs_err);
	assert(xs_err < 0.05*xs);
	}
	diff --git a/FixedOrderGen/t/h_5j.cc b/FixedOrderGen/t/h_5j.cc
	index 6aa3d40..a1b60e8 100644
	--- a/FixedOrderGen/t/h_5j.cc
	+++ b/FixedOrderGen/t/h_5j.cc
	@@ -1,63 +1,64 @@
	#ifdef NDEBUG
	#undef NDEBUG
	#endif
	// This is a regression test
	// the reference cross section has not been checked against any other program

	#include <algorithm>
	#include <cassert>
	#include <cmath>
	#include <iostream>

	#include "config.hh"
	#include "EventGenerator.hh"
	#include "HEJ/Ranlux64.hh"

	#include "HEJ/Event.hh"
	#include "HEJ/MatrixElement.hh"
	#include "HEJ/PDF.hh"

	using namespace HEJFOG;

	int main(){
	constexpr double invGeV2_to_pb = 389379292.;
	constexpr double xs_ref = 0.252273; // +- 0.00657742
	//calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20

	auto config = load_config("config_h_2j.yml");
	config.process.njets = 5;

	HEJ::Ranlux64 ran{};
	HEJFOG::EventGenerator generator{
	config.process,
	config.beam,
	HEJ::ScaleGenerator{
	config.scales.base,
	config.scales.factors,
	config.scales.max_ratio
	},
	config.jets,
	config.pdf_id,
	config.subleading_fraction,
	config.subleading_channels,
	config.particles_properties,
	config.Higgs_coupling,
	ran
	};

	double xs = 0., xs_err = 0.;
	for (int trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	if(generator.status() != good) continue;
	- ev.central().weight *= invGeV2_to_pb;
	- ev.central().weight /= config.events;
	+ assert(ev);
	+ ev->central().weight *= invGeV2_to_pb;
	+ ev->central().weight /= config.events;

	- xs += ev.central().weight;
	- xs_err += ev.central().weight*ev.central().weight;
	+ xs += ev->central().weight;
	+ xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << std::endl;

	assert(std::abs(xs - xs_ref) < 3*xs_err);
	assert(xs_err < 0.06*xs);
	}
	diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
	index 7ba72fa..2d11440 100644
	--- a/include/HEJ/Event.hh
	+++ b/include/HEJ/Event.hh
	@@ -1,342 +1,342 @@
	/** \file
	* \brief Declares the Event class and helpers
	*
	* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
	* \date 2019
	* \copyright GPLv2 or later
	*/

	#pragma once

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

	#include "HEJ/event_types.hh"
	#include "HEJ/Particle.hh"

	#include "fastjet/ClusterSequence.hh"

	namespace LHEF{
	class HEPEUP;
	class HEPRUP;
	}

	namespace fastjet{
	class JetDefinition;
	}

	namespace HEJ{


	struct ParameterDescription;

	//! Event parameters
	struct EventParameters{
	double mur; /*< Value of the Renormalisation Scale /
	double muf; /*< Value of the Factorisation Scale /
	double weight; /*< Event Weight /
	//! Optional description
	std::shared_ptr<ParameterDescription> description = nullptr;
	};

	//! Description of event parameters
	struct ParameterDescription {
	//! Name of central scale choice (e.g. "H_T/2")
	std::string scale_name;
	//! Actual renormalisation scale divided by central scale
	double mur_factor;
	//! Actual factorisation scale divided by central scale
	double muf_factor;

	ParameterDescription() = default;
	ParameterDescription(
	std::string scale_name, double mur_factor, double muf_factor
	):
	scale_name{scale_name}, mur_factor{mur_factor}, muf_factor{muf_factor}
	{};
	};

	struct UnclusteredEvent;

	/** An event with clustered jets
	*
	* This is the main HEJ 2 event class.
	* It contains kinematic information including jet clustering,
	* parameter (e.g. scale) settings and the event weight.
	*/
	class Event{
	public:
	class EventData;
	//! Default Event Constructor
	- Event() = default;
	+ Event() = delete;
	//! Event Constructor adding jet clustering to an unclustered event
	//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.3.0
	[[deprecated("UnclusteredEvent will be replaced by EventData")]]
	Event(
	UnclusteredEvent const & ev,
	fastjet::JetDefinition const & jet_def, double min_jet_pt
	);


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

	//! The corresponding event before jet clustering
	// [[deprecated]]
	// UnclusteredEvent unclustered() const {
	// return UnclusteredEvent(ev_);
	// TODO what to do with this?
	// }

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

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

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

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

	//! Indices of the jets the outgoing partons belong to
	/**
	* @param jets Jets to be tested
	* @returns A vector containing, for each outgoing parton,
	* the index in the vector of jets the considered parton
	* belongs to. If the parton is not inside any of the
	* passed jets, the corresponding index is set to -1.
	*/
	std::vector<int> particle_jet_indices(
	std::vector<fastjet::PseudoJet> const & jets
	) const{
	return cs_.particle_jet_indices(jets);
	}

	//! Jet definition used for clustering
	fastjet::JetDefinition const & jet_def() const{
	return cs_.jet_def();
	}

	//! Minimum jet transverse momentum
	double min_jet_pt() const{
	return min_jet_pt_;
	}

	//! Event type
	event_type::EventType type() const{
	return type_;
	}

	private:
	//! \internal Event Constructor adding jet clustering to an bare, unclustered event
	Event(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	std::unordered_map<size_t, std::vector<Particle>> const & decays,
	EventParameters const & central,
	std::vector<EventParameters> const & variations,
	fastjet::JetDefinition const & jet_def,
	double const min_jet_pt
	): incoming_{incoming},
	outgoing_{outgoing},
	decays_{decays},
	central_{central},
	variations_{variations},
	cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
	min_jet_pt_{min_jet_pt}
	{};

	//! \internal sort particles
	void sort();

	std::array<Particle, 2> incoming_;
	std::vector<Particle> outgoing_;
	std::unordered_map<size_t, std::vector<Particle>> decays_;
	//! @TODO replace this by "ParameterVariations"
	EventParameters central_;
	//! @TODO replace this by "ParameterVariations"
	std::vector<EventParameters> variations_;
	fastjet::ClusterSequence cs_;
	double min_jet_pt_;
	event_type::EventType type_;
	}; // end class Event

	//! Class to store general Event setup, i.e. Phase space and weights
	class Event::EventData{
	public:
	//! Default Constructor
	EventData() = default;
	//! Constructor from LesHouches event information
	EventData(LHEF::HEPEUP const & hepeup);
	//! Constructor with all values given
	EventData(
	std::array<Particle, 2> const & incoming,
	std::vector<Particle> const & outgoing,
	std::unordered_map<size_t, std::vector<Particle>> const & decays,
	EventParameters const & central,
	std::vector<EventParameters> const & variations
	):
	incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
	decays_(std::move(decays)),
	central_(std::move(central)), variations_(std::move(variations))
	{};
	//! Move Constructor with all values given
	EventData(
	std::array<Particle, 2> && incoming,
	std::vector<Particle> && outgoing,
	std::unordered_map<size_t, std::vector<Particle>> && decays,
	EventParameters && central,
	std::vector<EventParameters> && variations
	):
	incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
	decays_(std::move(decays)),
	central_(std::move(central)), variations_(std::move(variations))
	{};

	//! Generate an Event from the stored EventData.
	/**
	* @details Do jet clustering and classification.
	* Use this to generate an Event.
	*
	* @param jet_def Jet definition
	* @param min_jet_pt minimal \f$p_T\f$ for each jet
	*
	* @returns Full clustered and classified event.
	*/
	Event cluster(
	fastjet::JetDefinition const & jet_def, double const min_jet_pt) const;

	//! Alias for cluster()
	Event operator()(
	fastjet::JetDefinition const & jet_def, double const min_jet_pt) const{
	return cluster(jet_def, min_jet_pt);
	};

	//! Get Incoming Particles
	std::array<Particle, 2> const & get_incoming() const{
	return incoming_;
	}
	//! Set Incoming Particles
	void set_incoming(std::array<Particle, 2> const & in){
	incoming_ = in;
	}

	//! Get Outgoing Particles
	std::vector<Particle> const & get_outgoing() const{
	return outgoing_;
	}
	//! Set Outgoing Particles
	void set_outgoing(std::vector<Particle> const & out){
	outgoing_ = out;
	}

	//! Get Particle decays in the format {outgoing index, decay products}
	std::unordered_map<size_t, std::vector<Particle>> const & get_decays() const{
	return decays_;
	}
	//! Set Particle decays in the format {outgoing index, decay products}
	void set_decays(
	std::unordered_map<size_t, std::vector<Particle>> const & decays
	){
	decays_ = decays;
	}

	//! Get Central parameter (e.g. scale) choice
	EventParameters const & get_central() const{
	return central_;
	}
	//! Set Central parameter (e.g. scale) choice
	void set_central(EventParameters const & central){
	central_ = central;
	}

	//! Get parameter variation
	std::vector<EventParameters> const & get_variations() const{
	return variations_;
	}
	//! Set parameter variation
	void set_variations(std::vector<EventParameters> const & variations){
	variations_ = variations;
	}

	private:
	std::array<Particle, 2> incoming_;
	std::vector<Particle> outgoing_;
	std::unordered_map<size_t, std::vector<Particle>> decays_;
	//! @TODO replace this by "ParameterVariations"
	EventParameters central_;
	//! @TODO replace this by "ParameterVariations"
	std::vector<EventParameters> variations_;
	}; // end class EventData

	//! Square of the partonic centre-of-mass energy \f$\hat{s}\f$
	double shat(Event const & ev);

	//! Convert an event to a LHEF::HEPEUP
	LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *);

	// put deprecated warning at the end, so don't get the warning inside Event.hh,
	// additionally doxygen can not identify [[deprecated]] correctly
	struct [[deprecated("UnclusteredEvent will be replaced by EventData")]]
	UnclusteredEvent;
	//! An event before jet clustering
	//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.3.0
	struct UnclusteredEvent{
	//! Default Constructor
	UnclusteredEvent() = default;
	//! Constructor from LesHouches event information
	UnclusteredEvent(LHEF::HEPEUP const & hepeup);

	std::array<Particle, 2> incoming; /*< Incoming Particles /
	std::vector<Particle> outgoing; /*< Outgoing Particles /
	//! Particle decays in the format {outgoing index, decay products}
	std::unordered_map<size_t, std::vector<Particle>> decays;
	//! Central parameter (e.g. scale) choice
	EventParameters central;
	std::vector<EventParameters> variations; /*< For parameter variation /
	};

	}
	diff --git a/t/test_ME_generic.cc b/t/test_ME_generic.cc
	index 48c072f..86e9d12 100644
	--- a/t/test_ME_generic.cc
	+++ b/t/test_ME_generic.cc
	@@ -1,105 +1,108 @@

	// Generic tester for the ME for a given set of PSP
	// reference weights and PSP (as LHE file) have to be given as _individual_ files

	#include <fstream>

	#include "LHEF/LHEF.h"

	#include "HEJ/MatrixElement.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/YAMLreader.hh"
	#include "HEJ/stream.hh"

	constexpr double alpha_s = 0.118;
	constexpr double ep = 1e-6;

	void dump(HEJ::Event const & ev){
	{
	LHEF::Writer writer{std::cout};
	std::cout << std::setprecision(6);
	writer.hepeup = to_HEPEUP(std::move(ev), nullptr);
	writer.writeEvent();
	}
	std::cout << "Rapidity ordering:\n";
	for(const auto & part: ev.outgoing()){
	std::cout << std::setw(2) << part.type << ": "<< std::setw(7) << part.rapidity() << std::endl;
	}
	}

	int main(int argn, char** argv){
	if(argn != 4 && argn != 5){
	std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml ME_weights input_file.lhe\n\n";
	return EXIT_FAILURE;
	}
	bool OUTPUT_MODE = false;
	if(argn == 5 && std::string("OUTPUT")==std::string(argv[4]))
	OUTPUT_MODE = true;
	const HEJ::Config config = HEJ::load_config(argv[1]);

	std::fstream wgt_file;
	if ( OUTPUT_MODE ) {
	std::cout << "_______________________USING OUTPUT MODE!_______________________" << std::endl;
	wgt_file.open(argv[2], std::fstream::out);
	wgt_file.precision(10);
	} else {
	wgt_file.open(argv[2], std::fstream::in);
	}

	HEJ::istream in{argv[3]};
	LHEF::Reader reader{in};

	HEJ::MatrixElement ME{
	[](double){ return alpha_s; },
	HEJ::to_MatrixElementConfig(config)
	};
	double max_ratio = 0.;
	size_t idx_max_ratio = 0;
	- HEJ::Event ev_max_ratio;
	+ HEJ::Event ev_max_ratio(HEJ::UnclusteredEvent(),
	+ config.resummation_jets.def,
	+ config.resummation_jets.min_pt
	+ );
	double av_ratio = 0;

	size_t i = 0;
	while(reader.readEvent()){
	++i;

	HEJ::Event event{
	HEJ::Event::EventData{reader.hepeup}(
	config.resummation_jets.def,
	config.resummation_jets.min_pt
	)
	};
	const double our_ME = ME.tree(event).central;

	if ( OUTPUT_MODE ) {
	wgt_file << our_ME << std::endl;
	} else {
	std::string line;
	if(!std::getline(wgt_file,line)) break;
	const double ref_ME = std::stod(line);
	const double diff = std::abs(our_ME/ref_ME-1.);
	av_ratio+=diff;
	if( diff > max_ratio ) {
	max_ratio = diff;
	idx_max_ratio = i;
	ev_max_ratio = event;
	}
	if( diff > ep ){
	size_t precision(std::cout.precision());
	std::cout.precision(16);
	std::cout<< "Large difference in PSP " << i << "\nis: "<<our_ME << " should: " << ref_ME << " => difference: " << diff << std::endl;
	std::cout.precision(precision);
	dump(event);
	return EXIT_FAILURE;
	}
	}
	}
	wgt_file.close();
	if ( !OUTPUT_MODE ) {
	size_t precision(std::cout.precision());
	std::cout.precision(16);
	std::cout << "Avg ratio after " << i << " PSP: " << av_ratio/i << std::endl;
	std::cout << "maximal ratio at " << idx_max_ratio << ": " << max_ratio << std::endl;
	std::cout.precision(precision);
	}
	return EXIT_SUCCESS;
	}

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions