No OneTemporary
Actions

Size

38 KB

Subscribers

None

View Options

	diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
	index 3be02f3..461202a 100644
	--- a/FixedOrderGen/src/main.cc
	+++ b/FixedOrderGen/src/main.cc
	@@ -1,233 +1,233 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#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);
	}
	}

	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()];
	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()];
	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));
	if(! unweighted) continue;
	ev = std::move(unweighted);
	}
	events.emplace_back(std::move(*ev));
	++progress;
	}
	}
	std::cout << std::endl;

	HEJ::CrossSectionAccumulator xs;
	for(auto & ev: events){
	ev.parameters() *= 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";
	}
	-
	+ return EXIT_SUCCESS;
	}
	diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
	index 4733eb4..eb089aa 100644
	--- a/src/bin/HEJ.cc
	+++ b/src/bin/HEJ.cc
	@@ -1,239 +1,240 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <chrono>
	#include <iostream>
	#include <limits>
	#include <memory>
	#include <numeric>

	#include "yaml-cpp/yaml.h"

	#include "fastjet/ClusterSequence.hh"

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

	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;
	}
	fastjet::ClusterSequence::print_banner();

	// read configuration
	const HEJ::Config config = load_config(argv[1]);
	auto reader = HEJ::make_reader(argv[2]);
	assert(reader);

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

	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();

	analysis->initialise(heprup);
	HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};

	double global_reweight = 1.;
	int max_events = std::numeric_limits<int>::max();
	// if we need the event number:
	if(std::abs(heprup.IDWTUP) == 4 \|\| std::abs(heprup.IDWTUP) == 1 \|\| argn > 3){
	// try to read from LHE head
	auto input_events{reader->number_events()};
	if(!input_events) {
	// else count manually
	auto t_reader = HEJ::make_reader(argv[2]);
	input_events = 0;
	while(t_reader->read_event()) ++(*input_events);
	}
	if(std::abs(heprup.IDWTUP) == 4 \|\| std::abs(heprup.IDWTUP) == 1){
	// IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs
	std::cout << "Found IDWTUP " << heprup.IDWTUP << ": "
	<< "assuming \"cross section = average weight\".\n"
	<< "converting to \"cross section = sum of weights\" ";
	global_reweight /= *input_events;
	}
	if(argn > 3){
	// maximal number of events given
	max_events = std::stoi(argv[3]);
	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
	};

	// status infos & eye candy
	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;
	std::map<HEJ::StatusCode, int> status_counter;
	size_t total_trials = 0;

	// Loop over the events in the input file
	while(reader->read_event() && nevent < max_events){
	++nevent;

	// reweight events so that the total cross section is conserved
	auto hepeup = reader->hepeup();
	hepeup.setWeight(0, global_reweight * hepeup.weight());

	HEJ::Event::EventData event_data{hepeup};
	event_data.reconstruct_intermediate();

	// calculate HEJ weight
	HEJ::Event FO_event{
	std::move(event_data).cluster(
	config.fixed_order_jets.def, config.fixed_order_jets.min_pt
	)
	};
	if(FO_event.central().weight == 0) {
	static const bool warned_once = [argv,nevent](){
	std::cerr
	<< "WARNING: event number " << nevent
	<< " in " << argv[2] << " has zero weight. "
	"Ignoring this and all further events with vanishing weight.\n";
	return true;
	}();
	(void) warned_once; // shut up compiler warnings
	continue;
	}

	auto resummed_events = hej.reweight(FO_event, config.trials);
	for(auto const & s: hej.status())
	++status_counter[s];
	total_trials+=hej.status().size();
	++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';
	std::cout << '\t' << name(EventType::first_type) << ": "
	<< nevent_type[EventType::first_type]
	<< ", failed to reconstruct " << nfailed_type[EventType::first_type]
	<< '\n';
	for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){
	std::cout << '\t' << name(static_cast<EventType>(i)) << ": "
	<< nevent_type[i]
	<< ", failed to reconstruct " << nfailed_type[i]
	<< '\n';
	}

	std::cout << '\n' << xs << '\n';

	std::cout << "Generation statistic: "
	<< status_counter[HEJ::StatusCode::good] << "/" << total_trials
	<< " trials successful.\n";
	for(auto && entry: status_counter){
	const double fraction = static_cast<double>(entry.second)/total_trials;
	const int percent = std::round(100*fraction);
	std::cout << std::left << std::setw(17) << (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 << "] " <<std::setw(2)<<std::right<< percent << "%\n";
	}

	std::chrono::duration<double> run_time = (clock::now() - start_time);
	std::cout << "\nFinished " << 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";

	+ return EXIT_SUCCESS;
	}
	diff --git a/t/check_hepmc.cc b/t/check_hepmc.cc
	index 0f824ac..e26be89 100644
	--- a/t/check_hepmc.cc
	+++ b/t/check_hepmc.cc
	@@ -1,26 +1,27 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <iostream>
	#include <stdexcept>

	#include "HepMC3/ReaderAscii.h"

	static constexpr double ep = 1e-3;

	int main(int argn, char** argv) {
	if(argn != 2){
	std::cerr << "Usage: check_hepmc hepmc_file\n";
	return EXIT_FAILURE;
	}

	HepMC3::ReaderAscii input{argv[1]};
	if(input.failed()) throw std::runtime_error{"failed to open HepMC file"};
	while(true){
	HepMC3::GenEvent ev{};
	if ( !input.read_event(ev) \|\| ev.event_number() == 0 ) break;
	if(input.failed()) throw std::runtime_error{"failed to read HepMC event"};
	}
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/check_lhe.cc b/t/check_lhe.cc
	index aafc0be..b770227 100644
	--- a/t/check_lhe.cc
	+++ b/t/check_lhe.cc
	@@ -1,40 +1,41 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <iostream>
	#include <unordered_map>

	#include "HEJ/stream.hh"
	#include "LHEF/LHEF.h"

	static constexpr double ep = 1e-3;

	int main(int argn, char** argv) {
	if(argn != 2){
	std::cerr << "Usage: check_lhe lhe_file\n";
	return EXIT_FAILURE;
	}

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

	std::unordered_map<int, double> xsec_ref;
	for(int i=0; i < reader.heprup.NPRUP; ++i)
	xsec_ref[reader.heprup.LPRUP[i]] = 0.;

	while(reader.readEvent()){
	xsec_ref[reader.hepeup.IDPRUP] += reader.hepeup.weight();
	}

	for(size_t i = 0; i < xsec_ref.size(); ++i){
	double const ref = xsec_ref[reader.heprup.LPRUP[i]];
	double const calc = reader.heprup.XSECUP[i];
	std::cout << ref << '\t' << calc << '\n';
	if(std::abs(calc-ref) > ep*calc){
	std::cerr << "Cross sections deviate substantially";
	return EXIT_FAILURE;
	}
	}
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/check_res.cc b/t/check_res.cc
	index 0d04ea9..274427b 100644
	--- a/t/check_res.cc
	+++ b/t/check_res.cc
	@@ -1,165 +1,166 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <iostream>
	#include <math.h>

	#include "LHEF/LHEF.h"

	#include "HEJ/Event.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/Mixmax.hh"
	#include "HEJ/stream.hh"

	#define ASSERT(x) if(!(x)) { \
	std::cerr << "Assertion '" #x "' failed.\n"; \
	return EXIT_FAILURE; \
	}

	namespace{
	const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
	const fastjet::JetDefinition Born_jet_def{jet_def};
	constexpr double Born_jetptmin = 30;
	constexpr double extpartonptmin = 30;
	constexpr double max_ext_soft_pt_fraction = 0.1;
	constexpr double jetptmin = 35;
	constexpr bool log_corr = false;
	using EventTreatment = HEJ::EventTreatment;
	using namespace HEJ::event_type;
	HEJ::EventTreatMap treat{
	{no_2_jets, EventTreatment::discard},
	{bad_final_state, EventTreatment::discard},
	{non_resummable, EventTreatment::discard},
	{unof, EventTreatment::discard},
	{unob, EventTreatment::discard},
	{qqxexb, EventTreatment::discard},
	{qqxexf, EventTreatment::discard},
	{qqxmid, EventTreatment::discard},
	{FKL, EventTreatment::reweight}
	};

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

	int main(int argn, char** argv) {
	if(argn == 5 && std::string(argv[4]) == "unof"){
	--argn;
	treat[unof] = EventTreatment::reweight;
	treat[unob] = EventTreatment::discard;
	treat[FKL] = EventTreatment::discard;
	}
	if(argn == 5 && std::string(argv[4]) == "unob"){
	--argn;
	treat[unof] = EventTreatment::discard;
	treat[unob] = EventTreatment::reweight;
	treat[FKL] = EventTreatment::discard;
	}
	else if(argn == 5 && std::string(argv[4]) == "splitf"){
	--argn;
	treat[qqxexb] = EventTreatment::discard;
	treat[qqxexf] = EventTreatment::reweight;
	treat[FKL] = EventTreatment::discard;
	}
	else if(argn == 5 && std::string(argv[4]) == "splitb"){
	--argn;
	treat[qqxexb] = EventTreatment::reweight;
	treat[qqxexf] = EventTreatment::discard;
	treat[FKL] = EventTreatment::discard;
	}
	else if(argn == 5 && std::string(argv[4]) == "qqxmid"){
	--argn;
	treat[qqxmid] = EventTreatment::reweight;
	treat[FKL] = EventTreatment::discard;
	}
	if(argn != 4){
	std::cerr << "Usage: check_res eventfile xsection tolerance [uno]";
	return EXIT_FAILURE;
	}

	const double xsec_ref = std::stod(argv[2]);
	const double tolerance = std::stod(argv[3]);

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

	HEJ::PhaseSpacePointConfig psp_conf;
	psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin};
	psp_conf.min_extparton_pt = extpartonptmin;
	psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
	HEJ::MatrixElementConfig ME_conf;
	ME_conf.log_correction = log_corr;
	ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{};
	HEJ::EventReweighterConfig conf;
	conf.psp_config = std::move(psp_conf);
	conf.ME_config = std::move(ME_conf);
	conf.jet_param = psp_conf.jet_param;
	conf.treat = treat;

	reader.readEvent();
	const bool has_Higgs = std::find(
	begin(reader.hepeup.IDUP),
	end(reader.hepeup.IDUP),
	25
	) != end(reader.hepeup.IDUP);
	const double mu = has_Higgs?125.:91.188;
	HEJ::ScaleGenerator scale_gen{
	{{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1.
	};
	HEJ::Mixmax ran{};
	HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran};

	double xsec = 0.;
	double xsec_err = 0.;
	do{
	auto ev_data = HEJ::Event::EventData{reader.hepeup};
	ev_data.reconstruct_intermediate();
	HEJ::Event ev{
	ev_data.cluster(
	Born_jet_def, Born_jetptmin
	)
	};
	auto resummed_events = hej.reweight(ev, 100);
	for(auto const & ev: resummed_events) {
	ASSERT(correct_colour(ev));
	ASSERT(isfinite(ev.central().weight));
	xsec += ev.central().weight;
	xsec_err += ev.central().weight*ev.central().weight;
	}
	} while(reader.readEvent());
	xsec_err = std::sqrt(xsec_err);
	const double significance =
	std::abs(xsec - xsec_ref) / std::sqrt( xsec_errxsec_err + tolerancetolerance );
	std::cout << xsec_ref << " +/- " << tolerance << " ~ "
	<< xsec << " +- " << xsec_err << " => " << significance << " sigma\n";

	if(significance > 3.){
	std::cerr << "Cross section is off by over 3 sigma!\n";
	return EXIT_FAILURE;
	}
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/cmp_events.cc b/t/cmp_events.cc
	index 93f7e9f..0a3c833 100644
	--- a/t/cmp_events.cc
	+++ b/t/cmp_events.cc
	@@ -1,53 +1,53 @@
	#include "HEJ/stream.hh"
	#include "HEJ/Event.hh"

	#include "LHEF/LHEF.h"

	int main(int argn, char** argv) {
	if(argn != 3){
	std::cerr << "Usage: cmp_events eventfile eventfile_no_boson ";
	return EXIT_FAILURE;
	}
	HEJ::istream in{argv[1]};
	LHEF::Reader reader{in};
	HEJ::istream in_no_boson{argv[2]};
	LHEF::Reader reader_no_boson{in_no_boson};

	while(reader.readEvent()) {
	if(! reader_no_boson.readEvent()) {
	std::cerr << "wrong number of events in " << argv[2] << '\n';
	return EXIT_FAILURE;
	}
	const auto is_AWZH = [](HEJ::Particle const & p) {
	return is_AWZH_boson(p);
	};
	const HEJ::Event::EventData event_data{reader.hepeup};
	const auto boson = std::find_if(
	begin(event_data.outgoing), end(event_data.outgoing), is_AWZH
	);
	if(boson == end(event_data.outgoing)) {
	std::cerr << "no boson in event\n";
	return EXIT_FAILURE;
	}
	HEJ::Event::EventData data_no_boson{reader_no_boson.hepeup};
	data_no_boson.reconstruct_intermediate();
	const auto new_boson = std::find_if(
	begin(event_data.outgoing), end(event_data.outgoing), is_AWZH
	);
	if(new_boson == end(data_no_boson.outgoing)) {
	std::cerr << "failed to find reconstructed boson\n";
	return EXIT_FAILURE;
	}

	if(new_boson->type != boson->type) {
	std::cerr << "reconstructed wrong boson type\n";
	return EXIT_FAILURE;
	}

	if(!HEJ::nearby(new_boson->p, boson->p)) {
	std::cerr << "reconstructed boson momentum is off\n";
	return EXIT_FAILURE;
	}
	-
	}
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/test_classify_ref.cc b/t/test_classify_ref.cc
	index a53a0a3..7a8edb4 100644
	--- a/t/test_classify_ref.cc
	+++ b/t/test_classify_ref.cc
	@@ -1,95 +1,95 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <fstream>
	#include <random>
	#include <algorithm>

	#include "LHEF/LHEF.h"

	#include "HEJ/YAMLreader.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/EventReader.hh"

	namespace{
	// this is deliberately chosen bigger than in the generation,
	// to cluster multiple partons in one jet
	constexpr double min_jet_pt = 40.;
	const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};

	void shuffle_particles(HEJ::Event::EventData & ev) {
	static std::mt19937_64 ran{0};
	// incoming
	std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
	// outgoing (through index)
	auto old_outgoing = std::move(ev.outgoing);
	std::vector<size_t> idx(old_outgoing.size());
	std::iota(idx.begin(), idx.end(), 0);
	std::shuffle(begin(idx), end(idx), ran);
	ev.outgoing.clear();
	ev.outgoing.reserve(old_outgoing.size());
	for(size_t i: idx) {
	ev.outgoing.emplace_back(std::move(old_outgoing[i]));
	}
	// find decays again
	if(!ev.decays.empty()){
	auto old_decays = std::move(ev.decays);
	ev.decays.clear();
	for(size_t i=0; i<idx.size(); ++i) {
	auto decay = old_decays.find(idx[i]);
	if(decay != old_decays.end())
	ev.decays.emplace(i, std::move(decay->second));
	}
	}
	}
	}

	int main(int argn, char** argv) {
	if(argn != 3 && argn != 4){
	std::cerr << "Usage: " << argv[0]
	<< " reference_classification input_file.lhe\n";
	return EXIT_FAILURE;
	}
	bool OUTPUT_MODE = false;
	if(argn == 4 && std::string("OUTPUT")==std::string(argv[3]))
	OUTPUT_MODE = true;

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

	auto reader{ HEJ::make_reader(argv[2]) };

	while(reader->read_event()){
	HEJ::Event::EventData data{ reader->hepeup() };
	shuffle_particles(data);
	const HEJ::Event ev{
	data.cluster(
	jet_def, min_jet_pt
	)
	};
	if ( OUTPUT_MODE ) {
	ref_file << ev.type() << std::endl;
	} else {
	std::string line;
	if(!std::getline(ref_file,line)) break;
	const auto expected{static_cast<HEJ::event_type::EventType>(std::stoi(line))};

	if(ev.type() != expected){
	using HEJ::event_type::name;
	std::cerr << "wrong classification of event:\n" << ev
	<< "classified as " << name(ev.type())
	<< ", is " << name(expected) << '\n';
	return EXIT_FAILURE;
	}
	}
	}
	-
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/test_descriptions.cc b/t/test_descriptions.cc
	index 2d3c9d2..4ac462a 100644
	--- a/t/test_descriptions.cc
	+++ b/t/test_descriptions.cc
	@@ -1,67 +1,68 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <iostream>
	#include <cstddef>

	#include "HEJ/Event.hh"
	#include "HEJ/EventReweighter.hh"
	#include "HEJ/ScaleFunction.hh"

	#define ASSERT(x) if(!(x)) { \
	std::cerr << "Assertion '" #x "' failed.\n"; \
	return EXIT_FAILURE; \
	}

	int main() {
	constexpr double mu = 125.;
	HEJ::ScaleFunction fun{"125", HEJ::FixedScale{mu}};
	ASSERT(fun.name() == "125");

	HEJ::ScaleGenerator scale_gen{
	{std::move(fun)}, {0.5, 1, 2.}, 2.1
	};
	HEJ::Event::EventData tmp;
	tmp.outgoing.push_back(
	{HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.), {}}
	);
	tmp.outgoing.push_back(
	{HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.), {}}
	);
	HEJ::Event ev{
	tmp.cluster(
	fastjet::JetDefinition{fastjet::kt_algorithm, 0.4}, 20.
	)
	};

	auto rescaled = scale_gen(std::move(ev));
	ASSERT(rescaled.central().description->scale_name == "125");
	for(auto const & var: rescaled.variations()) {
	ASSERT(var.description->scale_name == "125");
	}
	ASSERT(rescaled.central().description->mur_factor == 1.);
	ASSERT(rescaled.central().description->muf_factor == 1.);

	ASSERT(rescaled.variations(0).description->mur_factor == 1.);
	ASSERT(rescaled.variations(0).description->muf_factor == 1.);

	ASSERT(rescaled.variations(1).description->mur_factor == 0.5);
	ASSERT(rescaled.variations(1).description->muf_factor == 0.5);

	ASSERT(rescaled.variations(2).description->mur_factor == 0.5);
	ASSERT(rescaled.variations(2).description->muf_factor == 1.);

	ASSERT(rescaled.variations(3).description->mur_factor == 1.);
	ASSERT(rescaled.variations(3).description->muf_factor == 0.5);

	ASSERT(rescaled.variations(4).description->mur_factor == 1.);
	ASSERT(rescaled.variations(4).description->muf_factor == 2.);

	ASSERT(rescaled.variations(5).description->mur_factor == 2.);
	ASSERT(rescaled.variations(5).description->muf_factor == 1.);

	ASSERT(rescaled.variations(6).description->mur_factor == 2.);
	ASSERT(rescaled.variations(6).description->muf_factor == 2.);
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/test_hdf5.cc b/t/test_hdf5.cc
	index 08f7c10..3f370e3 100644
	--- a/t/test_hdf5.cc
	+++ b/t/test_hdf5.cc
	@@ -1,45 +1,46 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <iostream>

	#include "HEJ/EventReader.hh"

	int main(int argc, char** argv) {
	if(argc != 2) {
	std::cerr << "Usage: " << argv[0] << " file.hdf5\n";
	return EXIT_FAILURE;
	}
	auto reader = HEJ::make_reader(argv[1]);
	if(
	reader->heprup().EBMUP != std::make_pair(7000., 7000.)
	\|\| reader->heprup().PDFSUP != std::make_pair(13000, 13000)
	) {
	std::cerr << "Read incorrect init parameters\n";
	return EXIT_FAILURE;
	}
	int nevent = 0;
	while(reader->read_event()) {
	++nevent;
	if(reader->hepeup().NUP != 13) {
	std::cerr << "Read wrong number of particles: "
	<< reader->hepeup().NUP << " != 13 in event " << nevent;
	return EXIT_FAILURE;
	}
	for(size_t i = 0; i < 2; ++i) {
	for(size_t j = 0; j < 2; ++j) {
	if(reader->hepeup().PUP[i][j] != 0) {
	std::cerr << "Non-vanishing transverse momentum in incoming particle"
	" in event " << nevent;
	return EXIT_FAILURE;
	}
	}
	}
	}
	if(nevent != 51200) {
	std::cerr << "Wrong number of events " << nevent << " != 51200\n";
	return EXIT_FAILURE;
	}
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/test_psp.cc b/t/test_psp.cc
	index c232813..92b6a2b 100644
	--- a/t/test_psp.cc
	+++ b/t/test_psp.cc
	@@ -1,68 +1,68 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "LHEF/LHEF.h"
	#include "HEJ/stream.hh"
	#include "HEJ/config.hh"
	#include "HEJ/event_types.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/PhaseSpacePoint.hh"
	#include "HEJ/Ranlux64.hh"

	namespace{
	constexpr int max_trials = 100;
	constexpr double extpartonptmin = 45.;
	constexpr double max_ext_soft_pt_fraction =
	std::numeric_limits<double>::infinity();
	const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
	constexpr double min_jet_pt = 50;
	};

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

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

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

	HEJ::Ranlux64 ran{};

	while(reader.readEvent()){
	const HEJ::Event ev{
	HEJ::Event::EventData{reader.hepeup}( jet_def, min_jet_pt )
	};
	for(int trial = 0; trial < max_trials; ++trial){
	HEJ::PhaseSpacePoint psp{ev, conf, ran};
	if(psp.weight() != 0){
	HEJ::Event::EventData tmp_ev;
	tmp_ev.incoming = psp.incoming();
	tmp_ev.outgoing = psp.outgoing();
	tmp_ev.parameters.central = {0,0,0};
	HEJ::Event out_ev{ tmp_ev(jet_def, min_jet_pt) };
	if(out_ev.type() != ev.type()){
	using HEJ::event_type::name;
	std::cerr << "Wrong class of phase space point:\n"
	"original event of class " << name(ev.type()) << ":\n";
	writer.hepeup = reader.hepeup;
	writer.writeEvent();
	std::cerr << "Phase space point of class " << name(out_ev.type()) << ":\n";
	writer.hepeup = to_HEPEUP(out_ev, &writer.heprup);
	writer.writeEvent();
	return EXIT_FAILURE;
	}
	}
	}
	}
	-
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/test_scale_arithmetics.cc b/t/test_scale_arithmetics.cc
	index 3a66bab..933971e 100644
	--- a/t/test_scale_arithmetics.cc
	+++ b/t/test_scale_arithmetics.cc
	@@ -1,120 +1,120 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <fstream>
	#include <random>
	#include <algorithm>

	#include "LHEF/LHEF.h"

	#include "HEJ/EventReweighter.hh"
	#include "HEJ/make_RNG.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/YAMLreader.hh"
	#include "HEJ/stream.hh"

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

	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;
	}
	}

	void shuffle_particles(HEJ::Event::EventData & ev) {
	static std::mt19937_64 ran{0};
	// incoming
	std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
	// outgoing (through index)
	auto old_outgoing = std::move(ev.outgoing);
	std::vector<size_t> idx(old_outgoing.size());
	std::iota(idx.begin(), idx.end(), 0);
	std::shuffle(begin(idx), end(idx), ran);
	ev.outgoing.clear();
	ev.outgoing.reserve(old_outgoing.size());
	for(size_t i: idx) {
	ev.outgoing.emplace_back(std::move(old_outgoing[i]));
	}
	// find decays again
	if(!ev.decays.empty()){
	auto old_decays = std::move(ev.decays);
	ev.decays.clear();
	for(size_t i=0; i<idx.size(); ++i) {
	auto decay = old_decays.find(idx[i]);
	if(decay != old_decays.end())
	ev.decays.emplace(i, std::move(decay->second));
	}
	}
	}

	int main(int argn, char** argv){
	if(argn != 3){
	std::cerr << "\n# Usage:\n."<< argv[0] <<" config.yml input_file.lhe\n\n";
	return EXIT_FAILURE;
	}
	HEJ::Config config = HEJ::load_config(argv[1]);
	config.scales = HEJ::to_ScaleConfig(
	YAML::Load("scales: [H_T, 1 * H_T, 2/2 * H_T, 2H_T/2, H_T/22, H_T/2/24, H_TH_T/H_T]")
	);

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

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

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

	HEJ::EventReweighter resum{
	reader.heprup,
	std::move(scale_gen),
	to_EventReweighterConfig(config),
	*ran
	};

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

	HEJ::Event::EventData data{reader.hepeup};
	shuffle_particles(data);

	HEJ::Event event{
	data.cluster(
	config.resummation_jets.def,
	config.resummation_jets.min_pt
	)
	};

	auto resummed = resum.reweight(event, config.trials);
	for(auto && ev: resummed) {
	for(auto &&var: ev.variations()) {
	if(std::abs(var.muf - ev.central().muf) > ep) {
	std::cerr
	<< std::setprecision(15)
	<< "unequal scales: " << var.muf
	<< " != " << ev.central().muf << '\n'
	<< "in resummed event:\n";
	dump(ev);
	std::cerr << "\noriginal event:\n";
	dump(event);
	return EXIT_FAILURE;
	}
	}
	}
	}
	-
	+ return EXIT_SUCCESS;
	}
	diff --git a/t/test_scale_import.cc b/t/test_scale_import.cc
	index 36b72e9..2fc7a28 100644
	--- a/t/test_scale_import.cc
	+++ b/t/test_scale_import.cc
	@@ -1,34 +1,35 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include <stdexcept>
	#include <iostream>

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

	int main(int argc, char** argv) {
	constexpr double ep = 1e-7;
	if (argc != 2) {
	throw std::logic_error{"wrong number of args"};
	}
	const HEJ::Config config = HEJ::load_config(argv[1]);

	HEJ::Event::EventData tmp;
	tmp.outgoing.push_back(
	{HEJ::ParticleID::gluon, fastjet::PtYPhiM(50., -1., 0.3, 0.), {}}
	);
	tmp.outgoing.push_back(
	{HEJ::ParticleID::gluon, fastjet::PtYPhiM(30., 1., -0.3, 0.), {}}
	);
	HEJ::Event ev{
	tmp(fastjet::JetDefinition{fastjet::kt_algorithm, 0.4}, 20.)
	};

	const double softest_pt = config.scales.base[0](ev);
	if(std::abs(softest_pt-30.) > ep){
	throw std::logic_error{"wrong softest pt"};
	}
	+ return EXIT_SUCCESS;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:27 PM (11 h, 41 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486589
Default Alt Text: (38 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions