No OneTemporary
Actions

Size

24 KB

Referenced Files

None

Subscribers

None

View Options

	diff --git a/FixedOrderGen/t/CMakeLists.txt b/FixedOrderGen/t/CMakeLists.txt
	index abb527f9..0cc08e5f 100644
	--- a/FixedOrderGen/t/CMakeLists.txt
	+++ b/FixedOrderGen/t/CMakeLists.txt
	@@ -1,185 +1,200 @@
	set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}")
	set(runcard_dir ${tst_dir}/runcards)
	+
	+function(add_long_test)
	+ set(options NONE)
	+ set(oneValueArgs NAME)
	+ set(multiValueArgs COMMAND)
	+ cmake_parse_arguments(TEST "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN})
	+ add_test(
	+ NAME ${TEST_NAME}_short
	+ COMMAND ${TEST_COMMAND} short
	+ )
	+ if(TEST_ALL)
	+ add_test(
	+ NAME ${TEST_NAME}
	+ COMMAND ${TEST_COMMAND}
	+ )
	+ endif()
	+endfunction()
	+
	foreach(tst W_reconstruct_enu W_2j_classify W_nj_classify)
	add_executable(test_${tst} ${tst_dir}/${tst}.cc)
	target_link_libraries(test_${tst} hejfog_lib)
	add_test(NAME ${tst}
	COMMAND test_${tst}
	WORKING_DIRECTORY ${runcard_dir}
	)
	endforeach()
	# this only tests if the runcard actually works, not if the result is correct
	add_test(
	NAME main_2j
	COMMAND HEJFOG ${runcard_dir}/2j.yml
	)
	add_test(
	NAME main_h2j
	COMMAND HEJFOG ${runcard_dir}/h2j.yml
	)
	add_test(
	NAME main_h2j_decay
	COMMAND HEJFOG ${runcard_dir}/h2j_decay.yml
	)
	add_test(
	NAME peakpt
	COMMAND HEJFOG ${runcard_dir}/2j_peak.yml
	)

	# check that uno emission doesn't change FKL xs
	add_executable(FKL_uno FKL_uno.cc)
	target_link_libraries(FKL_uno hejfog_lib)
	add_test(
	NAME FKL_uno
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND FKL_uno ${runcard_dir}/h3j_uno.yml 0.0243548 0.000119862
	)
	# xs tests
	add_executable(xs_gen xs_gen.cc)
	target_link_libraries(xs_gen hejfog_lib)
	## Higgs
	add_test(
	NAME xs_hqQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/hqQ.yml 1.612e-02 1.26303e-04
	)
	add_test(
	NAME xs_h2j
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND xs_gen ${runcard_dir}/h2j.yml 2.04928 0.00956022
	)
	add_test(
	NAME xs_h3j
	# calculated with HEJ revision bd4388fe55cbc3f5a7b6139096456c551294aa31
	COMMAND xs_gen ${runcard_dir}/h3j.yml 1.07807 0.0132409
	)
	+add_long_test(
	+ NAME xs_h3j_uno
	+ # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	+ COMMAND xs_gen ${runcard_dir}/h3j_uno.yml 0.00347538 3.85875e-05
	+)
	add_test(
	NAME xs_h5j
	# calculated with HEJ revision dbde2ffbb3b383ae6709b2424d8f0f9d5658270b
	COMMAND xs_gen ${runcard_dir}/h5j.yml 0.0112504 0.000199633
	)
	add_test(
	NAME xs_h2j_decay
	# calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	COMMAND xs_gen ${runcard_dir}/h2j_decay.yml 0.00466994 2.20995e-05
	)
	## pure jets
	add_test(
	NAME xs_qQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/qQ.yml 7.354e+05 1.905e+03
	)
	add_test(
	NAME xs_2j
	# calculated with "combined" HEJ svn r3480
	COMMAND xs_gen ${runcard_dir}/2j.yml 86.42031848e06 590570
	)
	add_test(
	NAME xs_3j_uno
	# calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1
	COMMAND xs_gen ${runcard_dir}/3j_uno.yml 900 14.3131
	)
	add_test(
	NAME xs_3j_qqx
	# calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1
	COMMAND xs_gen ${runcard_dir}/3j_qqx.yml 62040 1005
	)
	add_test(
	NAME xs_4j_qqx
	# calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1
	COMMAND xs_gen ${runcard_dir}/4j_qqx.yml 28936 550
	)
	add_test(
	NAME xs_4j
	# calculated with HEJ revision 13207b5f67a5f40a2141aa7ee515b022bd4efb65
	COMMAND xs_gen ${runcard_dir}/4j.yml 915072 15402.4
	)
	## W
	add_test(
	NAME xs_WqQ
	# calculated with Sherpa see #132
	COMMAND xs_gen ${runcard_dir}/WpqQ.yml 3.086e+00 4.511e-02
	)
	add_test(
	NAME xs_W2j
	# calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	COMMAND xs_gen ${runcard_dir}/Wm2j.yml 4.177443e+02 7.446928e+00
	)
	+add_long_test(
	+ NAME xs_W3j_eqqx
	+ # calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	+ COMMAND xs_gen ${runcard_dir}/Wp3j_qqx.yml 8.710751e+00 1.245725e-01
	+)
	+add_long_test(
	+ NAME xs_W3j_uno
	+ # calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	+ COMMAND xs_gen ${runcard_dir}/Wp3j_uno.yml 3.000955e-01 5.831799e-03
	+)
	+add_long_test(
	+ NAME xs_W4j_qqx
	+ # calculated with HEJ revision be065dc9a21e5965ce57583f6c0a3d953664b82b
	+ COMMAND xs_gen ${runcard_dir}/Wp4j_qqx.yml 9.274718e-02 4.875742e-03
	+)

	# Test that sum of partons == proton
	add_executable(PSP_channel PSP_channel.cc)
	target_link_libraries(PSP_channel hejfog_lib)
	# pure jets
	add_test(
	NAME channel_2j
	COMMAND PSP_channel ${runcard_dir}/2j.yml
	)
	add_test(
	NAME channel_3j_qqx
	COMMAND PSP_channel ${runcard_dir}/3j_qqx.yml
	)
	add_test(
	NAME channel_3j_uno
	COMMAND PSP_channel ${runcard_dir}/3j_uno.yml
	)
	+add_long_test(
	+ NAME channel_4j_qqx
	+ COMMAND PSP_channel ${runcard_dir}/4j_qqx.yml
	+)
	# W+jets
	# pure jets
	add_test(
	NAME channel_W2j
	COMMAND PSP_channel ${runcard_dir}/Wm2j.yml
	)
	add_test(
	NAME channel_W3j_uno
	COMMAND PSP_channel ${runcard_dir}/Wp3j_uno.yml
	)
	add_test(
	NAME channel_W3j_eqqx
	COMMAND PSP_channel ${runcard_dir}/Wp3j_qqx.yml
	)
	+add_long_test(
	+ NAME channel_W4j_qqx
	+ COMMAND PSP_channel ${runcard_dir}/Wp4j_qqx.yml
	+)

	# Test that each subleading channel is generated consistently
	add_executable(PSP_subl PSP_subl.cc)
	target_link_libraries(PSP_subl hejfog_lib)
	-
	-if(TEST_ALL) # deactivate long tests by default
	- add_test(
	- NAME xs_W3j_eqqx
	- # calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	- COMMAND xs_gen ${runcard_dir}/Wp3j_qqx.yml 8.710751e+00 1.245725e-01
	- )
	- add_test(
	- NAME xs_h3j_uno
	- # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20
	- COMMAND xs_gen ${runcard_dir}/h3j_uno.yml 0.00347538 3.85875e-05
	- )
	- add_test(
	- NAME xs_W4j_qqx
	- # calculated with HEJ revision be065dc9a21e5965ce57583f6c0a3d953664b82b
	- COMMAND xs_gen ${runcard_dir}/Wp4j_qqx.yml 9.274718e-02 4.875742e-03
	- )
	- add_test(
	- NAME xs_W3j_uno
	- # calculated with HEJ revision 449f2f6b597020e9c9e35699568edc05c827fc11+1
	- COMMAND xs_gen ${runcard_dir}/Wp3j_uno.yml 3.000955e-01 5.831799e-03
	- )
	- add_test(
	- NAME channel_4j_qqx
	- COMMAND PSP_channel ${runcard_dir}/4j_qqx.yml
	- )
	- add_test(
	- NAME channel_W4j_qqx
	- COMMAND PSP_channel ${runcard_dir}/Wp4j_qqx.yml
	- )
	- add_test(
	- NAME subl_5j
	- COMMAND PSP_subl ${runcard_dir}/5j.yml
	- )
	- add_test(
	- NAME subl_h5j
	- COMMAND PSP_subl ${runcard_dir}/h5j.yml
	- )
	- add_test(
	- NAME subl_W5j
	- COMMAND PSP_subl ${runcard_dir}/Wm5j.yml
	- )
	-endif()
	+add_long_test(
	+ NAME subl_5j
	+ COMMAND PSP_subl ${runcard_dir}/5j.yml
	+)
	+add_long_test(
	+ NAME subl_h5j
	+ COMMAND PSP_subl ${runcard_dir}/h5j.yml
	+)
	+add_long_test(
	+ NAME subl_W5j
	+ COMMAND PSP_subl ${runcard_dir}/Wm5j.yml
	+)
	diff --git a/FixedOrderGen/t/PSP_channel.cc b/FixedOrderGen/t/PSP_channel.cc
	index 3d9539a8..dfde2b09 100644
	--- a/FixedOrderGen/t/PSP_channel.cc
	+++ b/FixedOrderGen/t/PSP_channel.cc
	@@ -1,218 +1,222 @@
	/**
	* check that the sum of all possible quarks is the same as a proton
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2020
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <cmath>
	#include <iostream>
	#include <memory>

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

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

	//! throw error if condition not fulfilled
	#define ASSERT(x) if(!(x)) { \
	throw std::logic_error("Assertion '" #x "' failed."); \
	}

	//! throw error if condition not fulfilled
	#define ASSERT_THROW(x, exception) try { \
	x; \
	std::cerr << "'" #x "' did not throw an exception.\n"; \
	throw; \
	} catch(exception const & ex){ \
	std::cout << "throw triggered: " << ex.what() << std::endl; \
	} \
	catch (...) { \
	std::cerr << "Unexpected exception thrown for '" #x "'.\n"; \
	throw; \
	}

	namespace {

	HEJFOG::EventGenerator make_generator(
	HEJFOG::Config const & config,
	std::shared_ptr<HEJ::RNG> ran
	){
	return HEJFOG::EventGenerator{
	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.particle_decays,
	config.Higgs_coupling,
	config.ew_parameters,
	ran
	};
	}

	const static std::array<HEJ::ParticleID, 11> PARTONS{
	HEJ::ParticleID::g,
	HEJ::ParticleID::d, HEJ::ParticleID::u,
	HEJ::ParticleID::s, HEJ::ParticleID::c,
	HEJ::ParticleID::b,
	HEJ::ParticleID::d_bar, HEJ::ParticleID::u_bar,
	HEJ::ParticleID::s_bar, HEJ::ParticleID::c_bar,
	HEJ::ParticleID::b_bar
	};

	bool only_channel(
	HEJFOG::Subleading channels,
	HEJFOG::subleading::Channels selected
	){
	return channels[selected] && channels.reset(selected).none();
	}

	void generate_till_throw(HEJFOG::EventGenerator & gen ){
	for(std::size_t i=0; i<100; ++i){
	auto ev = gen.gen_event();
	if(gen.status()==HEJFOG::Status::good){
	std::cerr << "Wrongfully generated valid event!" << *ev <<"\n";
	throw;
	}
	}
	std::cerr << "Unable to generate a proper throw!\n";
	throw;
	}

	bool can_couple_to_W(HEJ::ParticleID const id, HEJ::ParticleID const boson){
	using namespace HEJ;
	if(!is_anyquark(id)){
	return false;
	}
	if(std::abs(id)==HEJ::ParticleID::b \|\| std::abs(id)==HEJ::ParticleID::top)
	return false;
	// Wp:
	if(boson>0){
	// anti-down
	if(id%2==0){
	return is_quark(id);
	}
	// or up-type quark
	return is_antiquark(id);
	}
	// Wm:
	// down
	if(id%2==0){
	return is_antiquark(id);
	}
	// or anti-up-type quark
	return is_quark(id);
	}
	}

	int main(int argc, char const *argv[])
	{
	- if(argc != 2){
	+ if(argc != 2 && argc != 3){
	std::cerr << "Usage: " << argv[0] << " config.yaml\n";
	return EXIT_FAILURE;
	}
	+ const bool short_only = argc==3;
	std::cout <<std::scientific;

	auto config = HEJFOG::load_config(argv[1]);
	config.events = 100000;
	auto ran = std::make_shared<HEJ::Mixmax>();

	+ if(short_only)
	+ config.events/=10.;
	+
	double tot_weight = 0.;
	double tot_err = 0.;
	config.process.incoming[0] = config.process.incoming[1] = HEJ::ParticleID::proton;
	{
	HEJFOG::EventGenerator gen = make_generator(config, ran);
	for(std::size_t i=0; i<config.events; ++i){
	auto const ev = gen.gen_event();
	if(ev){
	double const wgt = ev->central().weight/config.events;
	tot_weight += wgt;
	tot_err += wgt*wgt;
	}
	}
	ASSERT(tot_weight!=0.);
	}
	tot_err = std::sqrt(tot_err);


	config.events /= PARTONS.size();
	double tot_channels = 0.;
	double err_channels = 0.;
	for(auto b1: PARTONS){
	for(auto b2: PARTONS){
	double wgt_channel = 0.;
	double err_channel = 0.;
	config.process.incoming[0] = b1;
	config.process.incoming[1] = b2;
	HEJFOG::EventGenerator gen = make_generator(config, ran);
	std::cout << name(b1) << "+" << name(b2) << " " <<std::flush;
	// for pure subleading configurations some setups states should throw
	if(config.subleading_fraction == 1.){
	if( only_channel(config.subleading_channels, HEJFOG::subleading::uno) ){
	if(HEJ::is_gluon(b1) && HEJ::is_gluon(b2)){
	std::cout << "uno" << std::endl;
	ASSERT_THROW(generate_till_throw(gen), std::invalid_argument);
	continue;
	}
	}
	if( only_channel(config.subleading_channels, HEJFOG::subleading::eqqx) ){
	if(HEJ::is_anyquark(b1) && HEJ::is_anyquark(b2)){
	std::cout << "eqqx" << std::endl;
	ASSERT_THROW(generate_till_throw(gen), std::invalid_argument);
	continue;
	}
	}
	}
	// for W+jets we need the correct quarks
	if(config.process.boson
	&& std::abs(*config.process.boson) == HEJ::ParticleID::Wp
	){
	if(config.process.njets>3
	&& config.subleading_fraction>0.
	&& config.subleading_channels[HEJFOG::subleading::cqqx]){
	// this will force central qqx
	} else {
	auto const boson = *config.process.boson;
	if(!can_couple_to_W(b1, boson)
	&& !can_couple_to_W(b2, boson)
	){
	std::cout << "bad " << name(boson) << std::endl;
	ASSERT_THROW(generate_till_throw(gen), std::invalid_argument);
	continue;
	}
	}
	}
	for(std::size_t i=0; i<config.events; ++i){
	// everything else should work
	auto const ev = gen.gen_event();
	if(ev){
	double const wgt = ev->central().weight/config.events;
	wgt_channel += wgt;
	err_channel += wgt*wgt;
	}
	}
	ASSERT(wgt_channel!=0.);
	tot_channels += wgt_channel;
	err_channels += err_channel;
	err_channel = std::sqrt(err_channel);
	std::cout << "=> " << wgt_channel << " +/- " << err_channel << std::endl;
	}
	}
	err_channels = std::sqrt(err_channels);
	std::cout << tot_weight << " +/- " << tot_err
	<< " vs. " << tot_channels << " +/- " << err_channels << std::endl;

	ASSERT(std::abs(tot_weight - tot_channels)
	< 2.sqrt(err_channelserr_channels+tot_err*tot_err));

	return EXIT_SUCCESS;
	}
	diff --git a/FixedOrderGen/t/PSP_subl.cc b/FixedOrderGen/t/PSP_subl.cc
	index 56b9caba..e6de9025 100644
	--- a/FixedOrderGen/t/PSP_subl.cc
	+++ b/FixedOrderGen/t/PSP_subl.cc
	@@ -1,214 +1,218 @@
	/**
	* check that each subprocess xs is correctly generated when only selecting one
	*
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2020
	* \copyright GPLv2 or later
	*/
	#include <array>
	#include <cmath>
	#include <iostream>
	#include <memory>

	#include "HEJ/Mixmax.hh"
	#include "HEJ/Event.hh"
	#include "HEJ/CrossSectionAccumulator.hh"

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

	//! throw error if condition not fulfilled
	#define ASSERT(x) if(!(x)) { \
	throw std::logic_error("Assertion '" #x "' failed."); \
	}

	//! throw error if condition not fulfilled
	#define ASSERT_THROW(x, exception) try { \
	x; \
	std::cerr << "'" #x "' did not throw an exception.\n"; \
	throw; \
	} catch(exception const & ex){ \
	std::cout << "throw triggered: " << ex.what() << std::endl; \
	} \
	catch (...) { \
	std::cerr << "Unexpected exception thrown for '" #x "'.\n"; \
	throw; \
	}

	namespace {
	const static std::array<HEJ::event_type::EventType, 6> ALL_TYPES{
	HEJ::event_type::FKL,
	HEJ::event_type::unordered_backward,
	HEJ::event_type::unordered_forward,
	HEJ::event_type::extremal_qqxb,
	HEJ::event_type::extremal_qqxf,
	HEJ::event_type::central_qqx
	};

	HEJFOG::EventGenerator make_generator(
	HEJFOG::Config const & config,
	std::shared_ptr<HEJ::RNG> ran
	){
	return HEJFOG::EventGenerator{
	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.particle_decays,
	config.Higgs_coupling,
	config.ew_parameters,
	ran
	};
	}

	bool valid_type(
	HEJFOG::Subleading const channels, HEJ::event_type::EventType type
	){
	using namespace HEJ::event_type;
	using c = HEJFOG::subleading::Channels;
	if(channels.none())
	return type == FKL;
	if(channels.count()!=1)
	return false;
	if(channels[c::uno])
	return type == unob \|\| type == unof;
	if(channels[c::eqqx])
	return type == qqxexb \|\| type == qqxexf;
	if(channels[c::cqqx])
	return type == qqxmid;
	throw HEJ::unknown_option{"wrong channel"};
	}

	void generate_till_throw(HEJFOG::EventGenerator & gen ){
	for(std::size_t i=0; i<100; ++i){
	auto ev = gen.gen_event();
	if(gen.status()==HEJFOG::Status::good){
	std::cerr << "Wrongfully generated valid event!" << *ev <<"\n";
	throw;
	}
	}
	std::cerr << "Unable to generate a proper throw!\n";
	throw;
	}
	}

	int main(int argc, char const *argv[])
	{
	- if(argc != 2){
	+ if(argc != 2 && argc != 3){
	std::cerr << "Usage: " << argv[0] << " config.yaml\n";
	return EXIT_FAILURE;
	}
	+ const bool short_only = argc==3;
	std::cout <<std::scientific;

	auto config = HEJFOG::load_config(argv[1]);
	auto ran = std::make_shared<HEJ::Mixmax>();

	+ if(short_only)
	+ config.events/=10.;
	+
	// mixes sample
	HEJ::CrossSectionAccumulator xs_tot;
	config.process.incoming[0] = config.process.incoming[1] = HEJ::ParticleID::proton;
	config.subleading_fraction=0.8;
	if(config.process.boson && *config.process.boson == HEJ::ParticleID::Higgs){
	config.subleading_channels.reset();
	config.subleading_channels.set(HEJFOG::subleading::uno);
	} else {
	config.subleading_channels.set();
	}
	{
	HEJFOG::EventGenerator gen = make_generator(config, ran);
	for(std::size_t i=0; i<config.events; ++i){
	auto ev = gen.gen_event();
	if(ev){
	ev->central().weight /= config.events;
	xs_tot.fill(*ev);
	}
	}
	ASSERT(xs_tot.total().value!=0.);
	}

	// config.events /= HEJFOG::subleading::last+1;

	// pure FKL
	HEJ::CrossSectionAccumulator xs_subl;
	{
	config.subleading_fraction = 0.;
	HEJFOG::EventGenerator gen = make_generator(config, ran);
	for(std::size_t i=0; i<config.events; ++i){
	auto ev = gen.gen_event();
	if(ev){
	ev->central().weight /= config.events;
	xs_subl.fill(*ev);
	ASSERT(valid_type(0,ev->type()));
	}
	}
	ASSERT(xs_subl.total().value!=0.);
	std::cout << "=>\n" << xs_subl << std::endl;
	}

	// pure subleading
	config.subleading_fraction = 1.;
	for(unsigned channel = HEJFOG::subleading::first;
	channel <= HEJFOG::subleading::last; ++channel
	){
	config.subleading_channels.reset();
	config.subleading_channels.set(channel);
	HEJFOG::EventGenerator gen = make_generator(config, ran);
	HEJ::CrossSectionAccumulator xs_channel;
	std::cout << config.subleading_channels << " " <<std::flush;
	// Higgs+jets can only generate uno events
	if(config.process.boson
	&& *config.process.boson == HEJ::ParticleID::Higgs
	&& channel != HEJFOG::subleading::uno
	){
	ASSERT_THROW(generate_till_throw(gen), HEJ::not_implemented);
	continue;
	}
	for(std::size_t i=0; i<config.events; ++i){
	// everything else should work
	auto ev = gen.gen_event();
	if(ev){
	//! @FIXME for Higgs+jets we do actually generate FKL events ...
	if(config.process.boson
	&& *config.process.boson == HEJ::ParticleID::Higgs
	&& ev->type()==HEJ::event_type::FKL)
	continue;
	ev->central().weight /= config.events;
	xs_subl.fill(*ev);
	xs_channel.fill(*ev);
	ASSERT(valid_type(config.subleading_channels,ev->type()));
	}
	}
	ASSERT(xs_subl.total().value!=0.);
	std::cout << "=>\n" << xs_channel << std::endl;
	}
	std::cout << "Total:\n" << xs_tot << " vs.\n" << xs_subl << std::endl;

	for(auto type: ALL_TYPES){
	double diff = 0.;
	double err = 1.;
	try {
	auto const & tot = xs_tot[type];
	auto const & subl = xs_subl[type];
	diff = tot.value - subl.value;
	err = sqrt(tot.error+subl.error);
	} catch (std::out_of_range const &){
	std::cout << name(type) << " not set" << std::endl;
	ASSERT_THROW(xs_tot[type], std::out_of_range);
	ASSERT_THROW(xs_subl[type], std::out_of_range);
	continue;
	}
	if(std::abs(diff) > 2.*err){
	std::cerr << "Large difference in " << name(type)
	<< " (" << (diff/err) << " sigma)\n";
	return EXIT_FAILURE;
	}
	}

	return EXIT_SUCCESS;
	}
	diff --git a/FixedOrderGen/t/runcards/Wm5j.yml b/FixedOrderGen/t/runcards/Wm5j.yml
	index 7bd9e177..56f13e7c 100644
	--- a/FixedOrderGen/t/runcards/Wm5j.yml
	+++ b/FixedOrderGen/t/runcards/Wm5j.yml
	@@ -1,38 +1,38 @@
	-events: 1000000
	+events: 2000000

	jets:
	min pt: 50
	R: 0.4
	algorithm: antikt
	max rapidity: 5

	beam:
	energy: 15000
	particles: [p, p]

	pdf: 11000

	process: p p => W- 5j

	subleading fraction: 1.

	scales: 91

	particle properties:
	Higgs:
	mass: 125
	width: 0
	W:
	mass: 80.385
	width: 2.085
	Z:
	mass: 91.187
	width: 2.495

	decays:
	Wm: {into: [e-, electron_antineutrino], branching ratio: 1}

	random generator:
	name: mixmax

	vev: 246.2196508
	diff --git a/FixedOrderGen/t/xs_gen.cc b/FixedOrderGen/t/xs_gen.cc
	index 7566d25e..1b02581a 100644
	--- a/FixedOrderGen/t/xs_gen.cc
	+++ b/FixedOrderGen/t/xs_gen.cc
	@@ -1,108 +1,114 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2020
	* \copyright GPLv2 or later
	*/
	#include <algorithm>
	#include <cmath>
	#include <cstdlib>
	#include <iostream>
	#include <memory>

	#include "HEJ/Event.hh"
	#include "HEJ/Mixmax.hh"
	#include "HEJ/PDG_codes.hh"
	#include "HEJ/ScaleFunction.hh"
	#include "HEJ/utility.hh"

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

	//! throw error if condition not fulfilled
	#define ASSERT(x) if(!(x)) { \
	throw std::logic_error("Assertion '" #x "' failed."); \
	}

	namespace {
	constexpr double invGeV2_to_pb = 389379292.;
	}

	int main(int argn, char const *argv[]){
	using namespace HEJFOG;
	- if(argn != 4){
	+ if(argn != 4 && argn != 5){
	std::cerr << "Usage: " << argv[0] << " config.yml xs xs_err";
	return EXIT_FAILURE;
	}

	+ const bool short_only = argn==5;
	+
	const double xs_ref = std::stod(argv[2]);
	const double err_ref = std::stod(argv[3]);

	- const auto config = load_config(argv[1]);
	+ auto config = load_config(argv[1]);

	std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Mixmax>()};
	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.particle_decays,
	config.Higgs_coupling,
	config.ew_parameters,
	ran
	};
	+ if(short_only)
	+ config.events/=10.;

	double xs = 0., xs_err = 0.;
	for (std::size_t trials = 0; trials < config.events; ++trials){
	auto ev = generator.gen_event();
	if(generator.status() != Status::good) continue;
	ASSERT(ev);

	if(config.process.boson){
	const auto boson = std::find_if(
	begin(ev->outgoing()), end(ev->outgoing()),
	[&](HEJ::Particle const & p){ return p.type == *config.process.boson; }
	);
	ASSERT(boson != end(ev->outgoing()));

	if(!config.process.boson_decay.empty()){
	ASSERT(ev->decays().size() == 1);
	const auto decay = ev->decays().begin();
	ASSERT(ev->outgoing().size() > decay->first);
	ASSERT(decay->first == static_cast<std::size_t>(
	std::distance(ev->outgoing().begin(), boson)));
	ASSERT(decay->second.size() == 2);
	auto const & decay_part = decay->second;
	ASSERT(decay_part[0].type == config.process.boson_decay[0]);
	ASSERT(decay_part[1].type == config.process.boson_decay[1]);
	ASSERT(HEJ::nearby_ep(decay_part[0].p + decay_part[1].p, boson->p, 1e-6));
	}
	}

	xs += ev->central().weight;
	xs_err += ev->central().weight*ev->central().weight;
	}
	xs_err = std::sqrt(xs_err);
	xs *= invGeV2_to_pb/config.events;
	xs_err *= invGeV2_to_pb/config.events;
	std::cout <<std::scientific
	<< xs_ref << " +- " << err_ref
	<<std::fixed<< " (" << err_ref/xs_ref*100. << "%)"
	<<std::scientific<< " ~ " << xs << " +- " << xs_err
	<<std::fixed<< " (" << xs_err/xs*100. << "%)" << std::endl;
	std::cout << "=> "
	<< std::abs(xs - xs_ref)/sqrt(xs_errxs_err+err_referr_ref)
	<< " sigma" << std::endl;

	ASSERT(std::abs(xs - xs_ref) < 2.sqrt(xs_errxs_err+err_ref*err_ref));
	- ASSERT(std::abs(err_ref - xs_err) < 0.1*xs_err);
	+ if(!short_only){
	+ ASSERT(std::abs(err_ref - xs_err) < 0.1*xs_err);
	+ }

	return EXIT_SUCCESS;
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Sep 30, 6:06 AM (1 d, 22 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6566423
Default Alt Text: (24 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions