Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7877773
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
20 KB
Subscribers
None
View Options
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
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 4:23 PM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3798070
Default Alt Text
(20 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment