diff --git a/FixedOrderGen/t/CMakeLists.txt b/FixedOrderGen/t/CMakeLists.txt
index cd3bb57..c627280 100644
--- a/FixedOrderGen/t/CMakeLists.txt
+++ b/FixedOrderGen/t/CMakeLists.txt
@@ -1,209 +1,213 @@
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_example
COMMAND HEJFOG ${PROJECT_SOURCE_DIR}/configFO.yml
)
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
)
# special case where coupling not trivial
add_test(
NAME main_W4j_uno+eqqx
COMMAND HEJFOG ${runcard_dir}/Wp4j_uno+eqqx.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(
+add_long_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 667eb768fbefa99148bf6fe9ffb6fcd16c0f976e
COMMAND xs_gen ${runcard_dir}/Wp3j_qqx.yml 2.267785e+01 3.707774e-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 667eb768fbefa99148bf6fe9ffb6fcd16c0f976e
COMMAND xs_gen ${runcard_dir}/Wp4j_qqx.yml 1.159881e-01 3.633528e-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(
+add_long_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(
+add_long_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(
+add_long_test(
NAME channel_W2j
COMMAND PSP_channel ${runcard_dir}/Wm2j.yml
)
-add_test(
+add_long_test(
NAME channel_W3j_uno
COMMAND PSP_channel ${runcard_dir}/Wp3j_uno.yml
)
-add_test(
+add_long_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)
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_W3j
+ COMMAND PSP_subl ${runcard_dir}/Wp3j_qqx.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 dfde2b0..29ae1ca 100644
--- a/FixedOrderGen/t/PSP_channel.cc
+++ b/FixedOrderGen/t/PSP_channel.cc
@@ -1,222 +1,226 @@
/**
* 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 && 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;
+ config.events/=2.;
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]){
+ && config.subleading_fraction>0.
+ && config.subleading_channels[HEJFOG::subleading::cqqx]){
// this will force central qqx
+ } else if(config.process.njets>2
+ && config.subleading_fraction>0.
+ && config.subleading_channels[HEJFOG::subleading::eqqx]){
+ // this will force extremal 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_channels*err_channels+tot_err*tot_err));
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/PSP_subl.cc b/FixedOrderGen/t/PSP_subl.cc
index e6de902..d41bf08 100644
--- a/FixedOrderGen/t/PSP_subl.cc
+++ b/FixedOrderGen/t/PSP_subl.cc
@@ -1,218 +1,220 @@
/**
* 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 && 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
){
+ if(config.process.njets < 4 && channel == HEJFOG::subleading::cqqx)
+ continue;
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 56f13e7..77b3cd5 100644
--- a/FixedOrderGen/t/runcards/Wm5j.yml
+++ b/FixedOrderGen/t/runcards/Wm5j.yml
@@ -1,38 +1,38 @@
-events: 2000000
+events: 4000000
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