diff --git a/FixedOrderGen/t/CMakeLists.txt b/FixedOrderGen/t/CMakeLists.txt index 98b7c44..5de2656 100644 --- a/FixedOrderGen/t/CMakeLists.txt +++ b/FixedOrderGen/t/CMakeLists.txt @@ -1,166 +1,182 @@ set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}") set(runcard_dir ${tst_dir}/runcards) 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_test( NAME xs_h5j # calculated with HEJ revision dbde2ffbb3b383ae6709b2424d8f0f9d5658270b COMMAND xs_gen ${runcard_dir}/h5j.yml 0.0112504 0.000199633 ) 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_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_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 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_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_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_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_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 +) diff --git a/FixedOrderGen/t/PSP_subl.cc b/FixedOrderGen/t/PSP_subl.cc new file mode 100644 index 0000000..56b9cab --- /dev/null +++ b/FixedOrderGen/t/PSP_subl.cc @@ -0,0 +1,214 @@ +/** + * 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 +#include +#include +#include + +#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 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 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){ + std::cerr << "Usage: " << argv[0] << " config.yaml\n"; + return EXIT_FAILURE; + } + std::cout <(); + + // 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; icentral().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; icentral().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 << " " <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/5j.yml b/FixedOrderGen/t/runcards/5j.yml new file mode 100644 index 0000000..666455e --- /dev/null +++ b/FixedOrderGen/t/runcards/5j.yml @@ -0,0 +1,39 @@ +events: 500000 + +jets: + min pt: 30 + R: 0.4 + algorithm: antikt + max rapidity: 6 + +beam: + energy: 7000 + particles: [p, p] + +pdf: 11000 + +process: p p => 5j + +subleading fraction: 0.1 + +scales: 125 + +random generator: + name: mixmax + +particle properties: + Higgs: + mass: 125 + width: 0 + W: + mass: 80.385 + width: 2.085 + Z: + mass: 91.187 + width: 2.495 + +vev: 246.2196508 + +unweight: + sample size: 100 + max deviation: -10 diff --git a/FixedOrderGen/t/runcards/Wm5j.yml b/FixedOrderGen/t/runcards/Wm5j.yml new file mode 100644 index 0000000..7bd9e17 --- /dev/null +++ b/FixedOrderGen/t/runcards/Wm5j.yml @@ -0,0 +1,38 @@ +events: 1000000 + +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/include/HEJ/CrossSectionAccumulator.hh b/include/HEJ/CrossSectionAccumulator.hh index a8589ea..42bdd63 100644 --- a/include/HEJ/CrossSectionAccumulator.hh +++ b/include/HEJ/CrossSectionAccumulator.hh @@ -1,99 +1,107 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include "HEJ/event_types.hh" namespace HEJ { class Event; //! Collection of Cross Section with its uncertainty template struct XSWithError { T value = T{}; //!< Cross Section - T error = T{}; //!< Error + T error = T{}; //!< Squared error (Variance) }; /** * @brief Sum of Cross Section for different subproccess */ class CrossSectionAccumulator { public: //! Fill with single event //! @note for multiple resummation events use fill_correlated() instead void fill(Event const & ev); //! Fill by weight and type //! @note for multiple resummation events use fill_correlated() instead void fill(double wt, event_type::EventType type); //! Fill by weight, error and type //! @note The error will be _added_ to the current error void fill(double wt, double err, event_type::EventType type); /** * @brief Fill with multiple correlated weights * @details This should be used to fill multiple reweighted events, * coming from the same fixed order point. * Total error for \f$N\f$ fixed order points each giving \f$M_i\f$ * resummation events is: * \f[ * \delta^2=\sum_i \left(\sum_j w_{i,j}\right)^2 * +\sum_{i,j} \left(w_{i,j}\right)^2, * \f] * @note This is equivalent to fill() for only one reweighted event * coming from each fixed order point (\f$M_i=1\f$) */ void fill_correlated(std::vector const & evts); //! iterator implementation of fill_correlated() template void fill_correlated(ConstIt begin, ConstIt end); //! explicit version of fill_correlated() by giving sum(wt) and sum(wt^2) void fill_correlated(double sum, double sum2, event_type::EventType type); //! begin of Cross Section and error for subprocesses auto begin() const { return std::begin(xs_); } //! end of Cross Section and error for subprocesses auto end() const { return std::end(xs_); } //! total Cross Section and error XSWithError const & total() const { return total_; } + //! Cross Section and error of specific type + XSWithError const & at(event_type::EventType type) const { + return xs_.at(type); + } + //! Cross Section and error of specific type + XSWithError const & operator[](event_type::EventType type) const { + return at(type); + } private: std::map> xs_; XSWithError total_; }; //! Print CrossSectionAccumulator to stream std::ostream& operator<<(std::ostream& os, const CrossSectionAccumulator& xs); // ------------ Implementation ------------ template void CrossSectionAccumulator::fill_correlated(ConstIt begin, ConstIt end){ if(std::distance(begin, end) < 2){ // only one event fill(*begin); return; } double sum = 0.; double sum2 = 0.; const auto type = begin->type(); for(; begin != end; ++begin){ double const wt = begin->central().weight; sum += wt; sum2 += wt*wt; } fill_correlated(sum, sum2, type); } } // namespace HEJ