diff --git a/FixedOrderGen/t/CMakeLists.txt b/FixedOrderGen/t/CMakeLists.txt index 9503a66..aba4113 100644 --- a/FixedOrderGen/t/CMakeLists.txt +++ b/FixedOrderGen/t/CMakeLists.txt @@ -1,128 +1,147 @@ set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}") 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 ${tst_dir}) endforeach() # this only tests if the runcard actually works, not if the result is correct add_test( NAME main_2j COMMAND HEJFOG ${tst_dir}/config_2j.yml ) add_test( NAME main_h2j COMMAND HEJFOG ${tst_dir}/config_h2j.yml ) add_test( NAME main_h2j_decay COMMAND HEJFOG ${tst_dir}/config_h2j_decay.yml ) add_test( NAME peakpt COMMAND HEJFOG ${tst_dir}/config_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 ${tst_dir}/config_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_h2j # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 COMMAND xs_gen ${tst_dir}/config_h2j.yml 2.04928 0.00956022 ) add_test( NAME xs_h3j # calculated with HEJ revision bd4388fe55cbc3f5a7b6139096456c551294aa31 COMMAND xs_gen ${tst_dir}/config_h3j.yml 1.07807 0.0132409 ) add_test( NAME xs_h5j # calculated with HEJ revision dbde2ffbb3b383ae6709b2424d8f0f9d5658270b COMMAND xs_gen ${tst_dir}/config_h5j.yml 0.0112504 0.000199633 ) add_test( NAME xs_h3j_uno # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 COMMAND xs_gen ${tst_dir}/config_h3j_uno.yml 0.00347538 3.85875e-05 ) add_test( NAME xs_h2j_decay # calculated with HEJ revision 9570e3809613272ac4b8bf3236279ba23cf64d20 COMMAND xs_gen ${tst_dir}/config_h2j_decay.yml 0.00466994 2.20995e-05 ) ## pure jets add_test( NAME xs_2j # calculated with "combined" HEJ svn r3480 COMMAND xs_gen ${tst_dir}/config_2j.yml 86.42031848e06 590570 ) add_test( NAME xs_3j_uno # calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1 COMMAND xs_gen ${tst_dir}/config_3j_uno.yml 900 14.3131 ) add_test( NAME xs_3j_qqx # calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1 COMMAND xs_gen ${tst_dir}/config_3j_qqx.yml 62040 1005 ) add_test( NAME xs_4j_qqx # calculated with HEJ revision 9401196fba75b5d16bc33f2a309175fecaca00f1 COMMAND xs_gen ${tst_dir}/config_4j_qqx.yml 28936 550 ) add_test( NAME xs_4j # calculated with HEJ revision 13207b5f67a5f40a2141aa7ee515b022bd4efb65 COMMAND xs_gen ${tst_dir}/config_4j.yml 915072 15402.4 ) ## W add_test( NAME xs_W2j # calculated with HEJ revision 987bb30a7985a24a7961f98cfbbc74d3a7992970 COMMAND xs_gen ${tst_dir}/config_Wm2j.yml 231.404 3.43798 ) add_test( NAME xs_W3j_uno # calculated with HEJ revision 996f728a56ed67b980fccbe79948fe56914aaccd+1 COMMAND xs_gen ${tst_dir}/config_Wp3j_uno.yml 0.295275 0.0056667 ) add_test( NAME xs_W3j_eqqx # calculated with HEJ revision 987bb30a7985a24a7961f98cfbbc74d3a7992970 COMMAND xs_gen ${tst_dir}/config_Wp3j_qqx.yml 8.45323 0.103449 ) add_test( NAME xs_W4j_qqx # calculated with HEJ revision 987bb30a7985a24a7961f98cfbbc74d3a7992970 COMMAND xs_gen ${tst_dir}/config_Wp4j_qqx.yml 0.0851908 0.00300447 ) -# PSP tests +# 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 ${tst_dir}/config_2j.yml ) add_test( NAME channel_3j_qqx COMMAND PSP_channel ${tst_dir}/config_3j_qqx.yml ) add_test( NAME channel_3j_uno COMMAND PSP_channel ${tst_dir}/config_3j_uno.yml ) add_test( NAME channel_4j_qqx COMMAND PSP_channel ${tst_dir}/config_4j_qqx.yml ) +# W+jets +# pure jets +add_test( + NAME channel_W2j + COMMAND PSP_channel ${tst_dir}/config_Wm2j.yml +) +add_test( + NAME channel_W3j_uno + COMMAND PSP_channel ${tst_dir}/config_Wp3j_uno.yml +) +add_test( + NAME channel_W3j_eqqx + COMMAND PSP_channel ${tst_dir}/config_Wp3j_qqx.yml +) +add_test( + NAME channel_W4j_qqx + COMMAND PSP_channel ${tst_dir}/config_Wp4j_qqx.yml +) diff --git a/FixedOrderGen/t/PSP_channel.cc b/FixedOrderGen/t/PSP_channel.cc index 2c2e494..3d9539a 100644 --- a/FixedOrderGen/t/PSP_channel.cc +++ b/FixedOrderGen/t/PSP_channel.cc @@ -1,172 +1,218 @@ /** * 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 +#include #include #include #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 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 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){ std::cerr << "Usage: " << argv[0] << " config.yaml\n"; return EXIT_FAILURE; } std::cout <(); 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; icentral().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) << " " <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; icentral().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; }