diff --git a/FixedOrderGen/t/W_2j_classify.cc b/FixedOrderGen/t/W_2j_classify.cc index 414cb31..7b7ee17 100644 --- a/FixedOrderGen/t/W_2j_classify.cc +++ b/FixedOrderGen/t/W_2j_classify.cc @@ -1,165 +1,163 @@ /** * \brief check that the PSP generates only "valid" W + 2 jets events * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include #include #include #include #include "fastjet/JetDefinition.hh" #include "HEJ/EWConstants.hh" #include "HEJ/exceptions.hh" #include "HEJ/Mixmax.hh" #include "HEJ/Particle.hh" #include "HEJ/PDF.hh" #include "HEJ/PDG_codes.hh" #include "Decay.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" #include "Process.hh" #include "Status.hh" #include "Subleading.hh" namespace { using namespace HEJFOG; using namespace HEJ; void print_psp(PhaseSpacePoint const & psp){ std::cerr << "Process:\n" << psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> "; for(auto const & out: psp.outgoing()){ std::cerr << out.type << " "; } std::cerr << "\n"; } void bail_out(PhaseSpacePoint const & psp, std::string msg){ print_psp(psp); throw std::logic_error{msg}; } bool is_up_type(Particle const & part){ return HEJ::is_anyquark(part) && !(std::abs(part.type)%2); } bool is_down_type(Particle const & part){ return HEJ::is_anyquark(part) && std::abs(part.type)%2; } bool check_W2j(PhaseSpacePoint const & psp, ParticleID const W_type){ bool found_quark = false; bool found_anti = false; std::vector out_partons; std::vector Wp; for(auto const & p: psp.outgoing()){ if(p.type == W_type) Wp.push_back(p); else if(is_parton(p)) out_partons.push_back(p); else bail_out(psp, "Found particle with is not " +std::to_string(int(W_type))+" or parton"); } if(Wp.size() != 1 || out_partons.size() != 2){ bail_out(psp, "Found wrong number of outgoing partons"); } for(std::size_t j=0; j<2; ++j){ auto const & in = psp.incoming()[j]; auto const & out = out_partons[j]; if(is_quark(in) || is_antiquark(in)) { found_quark = true; if(in.type != out.type) { // switch in quark type -> Wp couples to it if(found_anti){ // already found qq for coupling to W bail_out(psp, "Found second up/down pair"); } else if(std::abs(in.type)>4 || std::abs(out.type)>4){ bail_out(psp, "Found bottom/top pair"); } found_anti = true; if( is_up_type(in)) { // "up" in if(W_type > 0){ // -> only allowed u -> Wp + d if(in.type < 0 || is_up_type(out) || out.type < 0) bail_out(psp, "u -/> Wp + d"); } else { // -> only allowed ux -> Wm + dx if(in.type > 0 || is_up_type(out) || out.type > 0) bail_out(psp, "ux -/> Wm + dx"); } } else { // "down" in if(W_type > 0){ // -> only allowed dx -> Wp + ux if(in.type > 0 || is_down_type(out) || out.type > 0) bail_out(psp, "dx -/> Wp + ux"); } else { // -> only allowed d -> Wm + u if(in.type < 0 || is_down_type(out) || out.type < 0) bail_out(psp, "d -/> Wm + u"); } } } } } if(!found_quark) { bail_out(psp, "Found no initial quarks"); } else if(!found_anti){ bail_out(psp, "Found no up/down pair"); } return true; } } int main(){ constexpr std::size_t n_psp_base = 1337; const JetParameters jet_para{ fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30, 5, 30}; PDF pdf(11000, pid::proton, pid::proton); constexpr double E_cms = 13000.; constexpr double subl_chance = 0.5; const Subleading subl_channels = subleading::ALL; const ParticlesDecayMap boson_decays{ {pid::Wp, {Decay{ {pid::e_bar, pid::nu_e}, 1.} }}, {pid::Wm, {Decay{ {pid::e, pid::nu_e_bar}, 1.} }} }; - const EWConstants ew_constants{246.2196508, - ParticleProperties{80.385, 2.085}, - ParticleProperties{91.187, 2.495}, - ParticleProperties{125, 0.004165} - }; + constexpr EWConstants ew_constants{246.2196508, {80.385, 2.085}, + {91.187, 2.495}, + {125, 0.004165} }; HEJ::Mixmax ran{}; // Wp2j Process proc {{pid::proton,pid::proton}, 2, pid::Wp, {}}; std::size_t n_psp = n_psp_base; for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wp+2j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; // Wm2j proc = Process{{pid::proton,pid::proton}, 2, pid::Wm, {}}; n_psp = n_psp_base; for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wm+2j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "All processes passed." << std::endl; return EXIT_SUCCESS; } diff --git a/FixedOrderGen/t/W_nj_classify.cc b/FixedOrderGen/t/W_nj_classify.cc index 2cabb8e..185bc91 100644 --- a/FixedOrderGen/t/W_nj_classify.cc +++ b/FixedOrderGen/t/W_nj_classify.cc @@ -1,227 +1,226 @@ /** * \brief check that the PSP generates the all W+jet subleading processes * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include #include #include #include #include #include #include #include #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/EWConstants.hh" #include "HEJ/exceptions.hh" #include "HEJ/Mixmax.hh" #include "HEJ/PDF.hh" #include "HEJ/PDG_codes.hh" #include "fastjet/JetDefinition.hh" #include "Decay.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" #include "Process.hh" #include "Status.hh" #include "Subleading.hh" namespace { using namespace HEJFOG; using namespace HEJ; void print_psp(PhaseSpacePoint const & psp){ std::cerr << "Process:\n" << psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> "; for(auto const & out: psp.outgoing()){ std::cerr << out.type << " "; } std::cerr << "\n"; } void bail_out(PhaseSpacePoint const & psp, std::string msg){ print_psp(psp); throw std::logic_error{msg}; } #if !defined(__clang__) && defined(__GNUC__) && (__GNUC__ < 6) // gcc version < 6 explicitly needs hash function for enum // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=60970 using type_map = std::unordered_map>; #else using type_map = std::unordered_map; #endif } int main(){ constexpr std::size_t n_psp_base = 10375; const JetParameters jet_para{ fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30, 5, 30}; PDF pdf(11000, pid::proton, pid::proton); constexpr double E_cms = 13000.; constexpr double subl_chance = 0.8; const ParticlesDecayMap boson_decays{ {pid::Wp, {Decay{ {pid::e_bar, pid::nu_e}, 1.} }}, {pid::Wm, {Decay{ {pid::e, pid::nu_e_bar}, 1.} }} }; - const EWConstants ew_constants{246.2196508, - ParticleProperties{80.385, 2.085}, - ParticleProperties{91.187, 2.495}, - ParticleProperties{125, 0.004165} + constexpr EWConstants ew_constants{246.2196508, {80.385, 2.085}, + {91.187, 2.495}, + {125, 0.004165} }; HEJ::Mixmax ran{}; Subleading subl_channels = subleading::ALL; std::vector allowed_types{event_type::FKL, event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf}; std::cout << "Wp3j" << std::endl; // Wp3j Process proc {{pid::proton,pid::proton}, 3, pid::Wp, {}}; std::size_t n_psp = n_psp_base; type_map type_counter; for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wp+3j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.05 * n_psp_base){ std::cerr << "Less than 5% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } // Wm3j - only uno proc = Process{{pid::proton,pid::proton}, 3, pid::Wm, {}}; n_psp = n_psp_base; subl_channels.reset(); subl_channels.set(subleading::uno); allowed_types = {event_type::FKL, event_type::unob, event_type::unof}; type_counter.clear(); for( std::size_t i = 0; i try again ++n_psp; } } std::cout << "Wm+3j (only uno): Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.05 * n_psp_base){ std::cerr << "Less than 5% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } // Wm4j proc = Process{{pid::proton,pid::proton}, 4, pid::Wm, {}}; n_psp = n_psp_base; subl_channels.set(); allowed_types = {event_type::FKL, event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf, event_type::qqxmid}; type_counter.clear(); std::array wpos_counter; // position of the W boson (back, central, forward) for( std::size_t i = 0; itype==pid::Wm){ ++wpos_counter[0][ev.type()]; } else if(ev.outgoing().rbegin()->type==pid::Wm){ ++wpos_counter[2][ev.type()]; } else { ++wpos_counter[1][ev.type()]; } } else { // bad process -> try again ++n_psp; } } std::cout << "Wm+4j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "States by classification:\n"; for(auto const & entry: type_counter){ const double fraction = static_cast(entry.second)/n_psp_base; const int percent = std::round(100*fraction); std::cout << std::left << std::setw(25) << (name(entry.first) + std::string(":")) << entry.second << " (" << percent << "%)\n"; } for(auto const & t: allowed_types){ if(type_counter[t] < 0.03 * n_psp_base){ std::cerr << "Less than 3% of the events are of type " << name(t) << std::endl; return EXIT_FAILURE; } } std::cout << "Stats by Wpos:\n"; for(std::size_t i=0; i +#include #include "HEJ/PDG_codes.hh" #include "HEJ/exceptions.hh" namespace HEJ { //! collection of basic particle properties struct ParticleProperties { double mass; //!< Mass double width; //!< Decay width }; //! Collection of electro-weak constants class EWConstants { public: + //! Default constructor, values have to be initialised separately EWConstants() = default; //! initialise by Vacuum expectation value & boson properties - EWConstants( + constexpr EWConstants( double vev, //!< vacuum expectation value - ParticleProperties const & Wprop, //!< W boson mass & width - ParticleProperties const & Zprop, //!< Z boson mass & width - ParticleProperties const & Hprop //!< Higgs boson mass & width - ): set_{true}, vev_{vev}, Wprop_{Wprop}, Zprop_{Zprop}, Hprop_{Hprop} {} + ParticleProperties Wprop, //!< W boson mass & width + ParticleProperties Zprop, //!< Z boson mass & width + ParticleProperties Hprop //!< Higgs boson mass & width + ): set_{true}, vev_{vev}, Wprop_{std::move(Wprop)}, + Zprop_{std::move(Zprop)}, Hprop_{std::move(Hprop)} + {} //! set constants by Vacuum expectation value & boson properties void set_vevWZH( double vev, //!< vacuum expectation value - ParticleProperties const & Wprop, //!< W boson mass & width - ParticleProperties const & Zprop, //!< Z boson mass & width - ParticleProperties const & Hprop //!< Higgs boson mass & width + ParticleProperties Wprop, //!< W boson mass & width + ParticleProperties Zprop, //!< Z boson mass & width + ParticleProperties Hprop //!< Higgs boson mass & width ){ - set_ = true; vev_= vev; Wprop_= Wprop; Zprop_= Zprop; Hprop_= Hprop; + set_ = true; vev_= vev; Wprop_= std::move(Wprop); + Zprop_= std::move(Zprop); Hprop_= std::move(Hprop); } //! vacuum expectation value double vev() const {check_set(); return vev_;} //! Properties of the W boson ParticleProperties const & Wprop() const {check_set(); return Wprop_;} //! Properties of the Z boson ParticleProperties const & Zprop() const {check_set(); return Zprop_;} //! Properties of the Higgs boson ParticleProperties const & Hprop() const {check_set(); return Hprop_;} //! access Properties by boson id ParticleProperties const & prop(ParticleID const id) const { using namespace pid; switch(id){ case Wp: case Wm: return Wprop(); case Z: return Zprop(); case h: return Hprop(); default: throw std::invalid_argument("No properties available for particle "+name(id)); } } //! cosine of Weinberg angle double cos_tw() const {return Wprop().mass/Zprop().mass;} //! cosine square of Weinberg angle double cos2_tw() const {return cos_tw()*cos_tw();} //! sinus Weinberg angle double sin_tw() const {return std::sqrt(sin2_tw());} //! sinus square of Weinberg angle double sin2_tw() const {return 1. - cos2_tw();} //! elector magnetic coupling double alpha_em() const {return e2()/4./M_PI;} //! weak coupling double alpha_w() const {return gw2()/2.;} private: double gw2() const {return 4*Wprop().mass/vev()*Wprop().mass/vev();} double e2() const {return gw2()*sin2_tw();} void check_set() const { if(!set_) throw std::invalid_argument("EW constants not specified"); } bool set_{false}; double vev_{}; ParticleProperties Wprop_{}; ParticleProperties Zprop_{}; ParticleProperties Hprop_{}; }; } // namespace HEJ diff --git a/t/CMakeLists.txt b/t/CMakeLists.txt index f0375ef..60c8e55 100644 --- a/t/CMakeLists.txt +++ b/t/CMakeLists.txt @@ -1,497 +1,505 @@ set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}") set(tst_ME_data_dir "${tst_dir}/ME_data") # small library for common test functions add_library(hej_test SHARED hej_test.cc) target_include_directories(hej_test PUBLIC ${tst_dir}) target_link_libraries(hej_test HEJ) # test event classification # test explicit configurations add_executable(test_classify ${tst_dir}/test_classify.cc) target_compile_options(test_classify PRIVATE "-O0") # avoid compiler optimisation target_link_libraries(test_classify HEJ hej_test) add_test( NAME classify COMMAND test_classify ) # test against reference data add_executable(test_classify_ref ${tst_dir}/test_classify_ref.cc) target_link_libraries(test_classify_ref HEJ hej_test) add_test( NAME classify_ref COMMAND test_classify_ref ${tst_dir}/classify_ref ${tst_dir}/classify.lhe.gz ) add_test( NAME classify_ref_4j COMMAND test_classify_ref ${tst_dir}/classify_ref_4j ${tst_dir}/4j.lhe.gz ) add_test( NAME classify_ref_W4j COMMAND test_classify_ref ${tst_dir}/classify_ref_W4j ${tst_dir}/MGemnubar4j_unweighted.lhe.gz ) add_test( NAME t_classify_ref_Z4j COMMAND test_classify_ref ${tst_dir}/classify_ref_Z4j ${tst_dir}/Z_4j.lhe.gz ) # test for valid W decays add_executable(test_decay ${tst_dir}/test_decay.cc) target_link_libraries(test_decay HEJ hej_test) add_test( NAME valid_decay COMMAND test_decay ) # test valid jet cuts on tagging jets add_executable(test_jet_cuts ${tst_dir}/test_jet_cuts.cc) target_link_libraries(test_jet_cuts HEJ hej_test) add_test( NAME jet_cuts COMMAND test_jet_cuts ) # test phase space point add_executable(test_psp ${tst_dir}/test_psp.cc) target_link_libraries(test_psp HEJ hej_test) add_test( NAME PhaseSpace COMMAND test_psp ${tst_dir}/psp_gen.lhe.gz ) # test importing analyses file(COPY "${tst_dir}/analysis_config.yml" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}") set(test_config "${CMAKE_CURRENT_BINARY_DIR}/analysis_config.yml") get_target_property(ANALYSIS_PATH AnalysisTemplate_lib BINARY_DIR) get_target_property(ANALYSIS_LIB AnalysisTemplate_lib OUTPUT_NAME) set(ANALYSIS_PARAMETERS "") file(READ ${test_config} config) file(WRITE ${test_config} "${config} - plugin: ${ANALYSIS_PATH}/lib${ANALYSIS_LIB}.so\n ${ANALYSIS_PARAMETERS}\n") configure_file( ${tst_dir}/analysis_config_single.yml.in analysis_config_simple.yml @ONLY ) add_test( NAME analysis_simple COMMAND $ analysis_config_simple.yml ${tst_dir}/2j.lhe.gz ) get_target_property(ANALYSIS_PATH AnalysisPrint_lib BINARY_DIR) get_target_property(ANALYSIS_LIB AnalysisPrint_lib OUTPUT_NAME) set(ANALYSIS_PARAMETERS " output: ana_output") file(READ ${test_config} config) file(WRITE ${test_config} "${config} - plugin: ${ANALYSIS_PATH}/lib${ANALYSIS_LIB}.so\n ${ANALYSIS_PARAMETERS}\n") configure_file( ${tst_dir}/analysis_config_single.yml.in analysis_config_print.yml @ONLY ) add_test( NAME analysis_print COMMAND $ analysis_config_print.yml ${tst_dir}/2j.lhe.gz ) if(ROOT_FOUND) get_target_property(ANALYSIS_PATH AnalysisROOT_lib BINARY_DIR) get_target_property(ANALYSIS_LIB AnalysisROOT_lib OUTPUT_NAME) set(ANALYSIS_PARAMETERS "") file(READ ${test_config} config) file(WRITE ${test_config} "${config} - plugin: ${ANALYSIS_PATH}/lib${ANALYSIS_LIB}.so\n ${ANALYSIS_PARAMETERS}\n") configure_file( ${tst_dir}/analysis_config_single.yml.in analysis_config_root.yml @ONLY ) add_test( NAME analysis_root COMMAND $ analysis_config_root.yml ${tst_dir}/2j.lhe.gz ) endif() if(RIVET_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config} - rivet: MC_XS\n output: ana_rivet\n") endif() add_test( NAME analysis_all COMMAND $ ${test_config} ${tst_dir}/2j.lhe.gz ) # test importing scales (from examples/softestptScale) add_executable(test_scale_import ${tst_dir}/test_scale_import) target_link_libraries(test_scale_import HEJ) get_target_property(SCALE_PATH softestptScale_lib BINARY_DIR) get_target_property(SCALE_LIB softestptScale_lib OUTPUT_NAME) set(SCALE_NAME "softest_jet_pt") configure_file( ${tst_dir}/jet_config_with_import.yml.in jet_config_with_import.yml @ONLY ) add_test( NAME scale_import COMMAND test_scale_import jet_config_with_import.yml ) # test scale arithmetic (e.g. 2*H_T/4) add_executable(test_scale_arithmetics ${tst_dir}/test_scale_arithmetics) target_link_libraries(test_scale_arithmetics HEJ hej_test) add_test( NAME scale_arithmetics COMMAND test_scale_arithmetics ${tst_dir}/jet_config.yml ${tst_dir}/2j.lhe.gz ) # test "ParameterDescription" add_executable(test_descriptions ${tst_dir}/test_descriptions) target_link_libraries(test_descriptions HEJ hej_test) add_test( NAME descriptions COMMAND test_descriptions ) # test "EventParameters*Weight" add_executable(test_parameters ${tst_dir}/test_parameters) target_link_libraries(test_parameters HEJ hej_test) add_test( NAME parameters COMMAND test_parameters ) # test Particle IDs add_executable(test_pid ${tst_dir}/test_pid) target_link_libraries(test_pid HEJ hej_test) add_test( NAME pid COMMAND test_pid ) +# test EW constants +add_executable(test_EWConstants ${tst_dir}/test_EWConstants) +target_link_libraries(test_EWConstants HEJ hej_test) +add_test( + NAME EWConstants + COMMAND test_EWConstants + ) + # test unweighting add_executable(test_unweighter ${tst_dir}/test_unweighter) target_link_libraries(test_unweighter HEJ hej_test) add_test( NAME unweighter COMMAND test_unweighter ${tst_dir}/4j.lhe.gz ) # test colour generation add_executable(test_colours ${tst_dir}/test_colours) target_link_libraries(test_colours HEJ hej_test) add_test( NAME colour_flow COMMAND test_colours ) add_executable(test_colours2 ${tst_dir}/test_colours2) target_link_libraries(test_colours2 HEJ hej_test) add_test( NAME colour_flow2 COMMAND test_colours2 ) # test matrix elements add_executable(test_ME_generic ${tst_dir}/test_ME_generic.cc) target_link_libraries(test_ME_generic HEJ hej_test) add_test( NAME ME_j COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets_tree.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz ) add_test( NAME ME_j_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_jets_tree_virt.dat ${tst_ME_data_dir}/PSP_jets.lhe.gz ) add_test( NAME ME_h COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) add_test( NAME ME_h_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtinf.yml ${tst_ME_data_dir}/ME_h_mtinf_virt.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) if(QCDloop_FOUND) add_test( NAME ME_h_mt COMMAND test_ME_generic ${tst_ME_data_dir}/config_mt.yml ${tst_ME_data_dir}/ME_h_mt_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) add_test( NAME ME_h_mtmb COMMAND test_ME_generic ${tst_ME_data_dir}/config_mtmb.yml ${tst_ME_data_dir}/ME_h_mtmb_tree.dat ${tst_ME_data_dir}/PSP_h.lhe.gz ) endif() add_test( NAME ME_j_subl COMMAND test_ME_generic ${tst_ME_data_dir}/config_pure.yml ${tst_ME_data_dir}/ME_jets_subl_tree.dat ${tst_ME_data_dir}/PSP_jets_subl.lhe.gz ) add_test( NAME ME_j_subl_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_pure.yml ${tst_ME_data_dir}/ME_jets_subl_tree_virt.dat ${tst_ME_data_dir}/PSP_jets_subl.lhe.gz ) add_test( NAME ME_j_subl_4j COMMAND test_ME_generic ${tst_ME_data_dir}/config_pure.yml ${tst_ME_data_dir}/ME_jets_subl_tree_new.dat ${tst_dir}/4j.lhe.gz ) add_test( NAME ME_j_subl_4j_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_pure.yml ${tst_ME_data_dir}/ME_jets_subl_tree_new_virt.dat ${tst_dir}/4j.lhe.gz ) add_test( NAME ME_w_FKL COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_tree.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz ) add_test( NAME ME_w_FKL_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_w_FKL_virt.dat ${tst_ME_data_dir}/PSP_w_FKL.lhe.gz ) add_test( NAME ME_Wp COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wp.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz ) add_test( NAME ME_Wp_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wp_virt.dat ${tst_ME_data_dir}/PSP_Wp.lhe.gz ) add_test( NAME ME_Wm COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wm.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz ) add_test( NAME ME_Wm_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_w_ME.yml ${tst_ME_data_dir}/ME_Wm_virt.dat ${tst_ME_data_dir}/PSP_Wm.lhe.gz ) add_test( NAME t_ME_Z_FKL COMMAND test_ME_generic ${tst_ME_data_dir}/config_Z_ME.yml ${tst_ME_data_dir}/ME_Z_FKL_tree.dat ${tst_ME_data_dir}/PSP_Z_FKL.lhe.gz ) add_test( NAME t_ME_Z_FKL_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_Z_ME.yml ${tst_ME_data_dir}/ME_Z_FKL_virt.dat ${tst_ME_data_dir}/PSP_Z_FKL.lhe.gz ) add_test( NAME t_ME_Z_uno COMMAND test_ME_generic ${tst_ME_data_dir}/config_Z_ME.yml ${tst_ME_data_dir}/ME_Z_uno_tree.dat ${tst_ME_data_dir}/PSP_Z_uno.lhe.gz ) add_test( NAME t_ME_Z_uno_virt COMMAND test_ME_generic ${tst_ME_data_dir}/config_Z_ME.yml ${tst_ME_data_dir}/ME_Z_uno_virt.dat ${tst_ME_data_dir}/PSP_Z_uno.lhe.gz ) # test main executable file(COPY "${tst_dir}/jet_config.yml" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}") set(test_config "${CMAKE_CURRENT_BINARY_DIR}/jet_config.yml") if(HighFive_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config} - tst.hdf5\n") endif() if(HepMC3_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config} - tst.hepmc\n") endif() if(HepMC_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config} - tst.hepmc2\n") endif() if(rivet_FOUND) file(READ ${test_config} config) file(WRITE ${test_config} "${config}\nanalysis:\n rivet: MC_XS\n output: tst\n") endif() set(test_cmd_main "$\\\;${test_config}\\\;${tst_dir}/2j.lhe.gz") # check that HepMC3 output is correct if(HepMC3_FOUND) add_executable(check_hepmc ${tst_dir}/check_hepmc.cc) target_link_libraries(check_hepmc HEJ ${HEPMC3_LIBRARIES}) target_include_directories(check_hepmc PRIVATE ${HEPMC3_INCLUDE_DIR}) set(test_cmd_hepmc "$\\\;tst.hepmc") else() set(test_cmd_hepmc "") endif() # check that LHEF output is correct add_executable(check_lhe ${tst_dir}/check_lhe.cc) target_link_libraries(check_lhe HEJ hej_test) set(test_cmd_lhe "$\\\;tst.lhe") # check that rivet interface is consistent with naive rivet if(rivet_FOUND) # this assumes "rivet" and "yodadiff" are found in PATH if(rivet_USE_HEPMC3) set(hepmc_file "tst.hepmc") else() set(hepmc_file "tst.hepmc2") endif() if(rivet_USE_HEPMC3 OR (rivet_VERSION VERSION_LESS 3)) set(histo_exclude "") else() # rivet 3 with HepMC 2 is inconsistent in order of weights # -> interface != direct call (by permutation) # REQUIRES Yoda 1.7.5 set(histo_exclude "-M\\\;\\\\d") endif() set(test_cmd_rivet "rivet\\\;-a\\\;MC_XS\\\;${hepmc_file}\\\;-o\\\;tst_direct.yoda\ \;yodadiff\\\;${histo_exclude}\\\;tst.yoda\\\;tst_direct.yoda") else() set(test_cmd_rivet "") endif() # Run dependent tests in one command to ensure correct execution order # Note: The commands are concatenated with "\;" to escape CMake lists. # Thus arguments have to be escaped twice "\\\;". # e.g. "cmd1\;cmd2\\\;arg1\\\;arg2" is executed like "cmd1 && cmd2 arg1 arg2" add_test( NAME main COMMAND ${CMAKE_COMMAND} -DCMDS=${test_cmd_main}\;${test_cmd_hepmc}\;${test_cmd_lhe}\;${test_cmd_rivet} -P ${PROJECT_SOURCE_DIR}/cmake/run_multiple_tests.cmake ) add_test( NAME main_example COMMAND $ ${PROJECT_SOURCE_DIR}/config.yml ${tst_dir}/h_3j.lhe.gz ) # check that Sherpas LHE input can be read add_executable(check_lhe_sherpa ${tst_dir}/check_lhe_sherpa.cc) target_link_libraries(check_lhe_sherpa HEJ hej_test) add_test( NAME sherpa_reader COMMAND check_lhe_sherpa ${tst_dir}/SherpaLHE.lhe 1.62624e+08 ) # check HDF5 reader & writer if(HighFive_FOUND) add_executable(test_hdf5 ${tst_dir}/test_hdf5.cc) target_link_libraries(test_hdf5 HEJ) add_test( NAME hdf5_read COMMAND test_hdf5 ${tst_dir}/Wm9-g4-repack.hdf5 ) add_executable(test_hdf5_write ${tst_dir}/test_hdf5_write.cc) target_link_libraries(test_hdf5_write HEJ hej_test) add_test( NAME hdf5_write COMMAND test_hdf5_write ${tst_dir}/Wm9-g4-repack.hdf5 output.hdf5 ) endif() # check rivet interface if(RIVET_FOUND) add_executable(check_rivet ${tst_dir}/check_rivet.cc) target_link_libraries(check_rivet HEJ rivet::rivet) add_test( NAME rivet COMMAND check_rivet ) endif() # test boson reconstruction add_executable(cmp_events ${tst_dir}/cmp_events.cc) target_link_libraries(cmp_events HEJ) add_test( NAME reconstruct_W COMMAND cmp_events ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz ${tst_dir}/epnu2jLOFKL_unweight_noW.lhe.gz ) # test resummed result add_executable(check_res ${tst_dir}/check_res.cc) target_link_libraries(check_res HEJ hej_test) if(TEST_ALL) # deactivate long tests by default add_test( NAME xs_2j COMMAND check_res ${tst_dir}/2j.lhe.gz 3.49391e+07 419684 ) add_test( NAME xs_3j COMMAND check_res ${tst_dir}/3j.lhe.gz 2.37902e+06 25746.6 ) add_test( NAME xs_3j_unof COMMAND check_res ${tst_dir}/3j.lhe.gz 133399 4688.83 unof ) add_test( NAME xs_3j_unob COMMAND check_res ${tst_dir}/3j.lhe.gz 105247 3449.45 unob ) add_test( NAME xs_3j_splitf COMMAND check_res ${tst_dir}/3j.lhe.gz 97659.9 2748.65 splitf ) add_test( NAME xs_3j_splitb COMMAND check_res ${tst_dir}/3j.lhe.gz 107150 2799.8 splitb ) add_test( NAME xs_4j COMMAND check_res ${tst_dir}/4j.lhe.gz 603713 72822.6 ) add_test( NAME xs_4j_qqxmid COMMAND check_res ${tst_dir}/4j.lhe.gz 21866.7 1716.96 qqxmid ) add_test( NAME xs_h_3j COMMAND check_res ${tst_dir}/h_3j.lhe.gz 0.821622 0.0220334 ) add_test( NAME xs_h_3j_unof COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0127362 0.000271555 unof ) add_test( NAME xs_h_3j_unob COMMAND check_res ${tst_dir}/h_3j_uno.lhe.gz 0.0130615 0.000224793 unob ) add_test( NAME xs_epnu_2j COMMAND check_res ${tst_dir}/epnu2jLOFKL_unweight.lhe.tar.gz 262.7 3 ) add_test( NAME xs_epnu_3j COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 38.9512 1 ) add_test( NAME xs_emnubar_3j COMMAND check_res ${tst_dir}/MGemnubar3j_unweighted.lhe.gz 24.1575 1 ) add_test( NAME xs_epnu_3j_unof COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 9.63702 0.128355 unof ) add_test( NAME xs_epnu_3j_unob COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 9.70119 0.108436 unob ) add_test( NAME xs_epnu_3j_splitf COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 2.91995 0.0463182 splitf ) add_test( NAME xs_epnu_3j_splitb COMMAND check_res ${tst_dir}/MGepnu3j_unweighted.lhe.gz 3.40708 0.0550975 splitb ) add_test( NAME xs_epnu_4j COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.lhe.gz 10.2542 0.135106 ) add_test( NAME xs_emnubar_4j COMMAND check_res ${tst_dir}/MGemnubar4j_unweighted.lhe.gz 5.73282 0.0583738 ) add_test( NAME xs_epnu_4j_qqxmid COMMAND check_res ${tst_dir}/MGepnu4j_unweighted.lhe.gz 0.802603 0.0126908 qqxmid ) endif() diff --git a/t/check_res.cc b/t/check_res.cc index 11837c6..101a4fa 100644 --- a/t/check_res.cc +++ b/t/check_res.cc @@ -1,166 +1,166 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ - #include "hej_test.hh" +#include "hej_test.hh" #include #include #include #include #include #include #include #include #include "HEJ/Config.hh" #include "HEJ/CrossSectionAccumulator.hh" #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/EventReweighter.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Fraction.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/Mixmax.hh" #include "HEJ/Parameters.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/stream.hh" #include "fastjet/JetDefinition.hh" #include "LHEF/LHEF.h" namespace HEJ { struct RNG; } namespace { const fastjet::JetDefinition JET_DEF{fastjet::kt_algorithm, 0.4}; const fastjet::JetDefinition BORN_JET_DEF{JET_DEF}; constexpr double BORN_JETPTMIN = 30; constexpr double MAX_EXT_SOFT_PT_FRACTION = 0.1; constexpr double JETPTMIN = 35; constexpr bool LOG_CORR = false; constexpr std::size_t NUM_TRIES = 100; - const HEJ::ParticleProperties WPROP{80.385, 2.085}; - const HEJ::ParticleProperties ZPROP{91.187, 2.495}; - const HEJ::ParticleProperties HPROP{125, 0.004165}; + constexpr HEJ::ParticleProperties WPROP{80.385, 2.085}; + constexpr HEJ::ParticleProperties ZPROP{91.187, 2.495}; + constexpr HEJ::ParticleProperties HPROP{125, 0.004165}; constexpr double VEV = 246.2196508; using EventTreatment = HEJ::EventTreatment; using namespace HEJ::event_type; HEJ::EventTreatMap TREAT{ {no_2_jets, EventTreatment::discard}, {bad_final_state, EventTreatment::discard}, {non_resummable, EventTreatment::discard}, {unof, EventTreatment::discard}, {unob, EventTreatment::discard}, {qqxexb, EventTreatment::discard}, {qqxexf, EventTreatment::discard}, {qqxmid, EventTreatment::discard}, {FKL, EventTreatment::reweight} }; bool correct_colour(HEJ::Event const & ev){ if(!HEJ::event_type::is_resummable(ev.type())) return true; return ev.is_leading_colour(); } } // namespace int main(int argn, char** argv) { if(argn == 5 && std::string(argv[4]) == "unof"){ --argn; TREAT[unof] = EventTreatment::reweight; TREAT[unob] = EventTreatment::discard; TREAT[FKL] = EventTreatment::discard; } if(argn == 5 && std::string(argv[4]) == "unob"){ --argn; TREAT[unof] = EventTreatment::discard; TREAT[unob] = EventTreatment::reweight; TREAT[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitf"){ --argn; TREAT[qqxexb] = EventTreatment::discard; TREAT[qqxexf] = EventTreatment::reweight; TREAT[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "splitb"){ --argn; TREAT[qqxexb] = EventTreatment::reweight; TREAT[qqxexf] = EventTreatment::discard; TREAT[FKL] = EventTreatment::discard; } else if(argn == 5 && std::string(argv[4]) == "qqxmid"){ --argn; TREAT[qqxmid] = EventTreatment::reweight; TREAT[FKL] = EventTreatment::discard; } if(argn != 4){ std::cerr << "Usage: check_res eventfile xsection tolerance [uno]"; return EXIT_FAILURE; } const double xsec_ref = std::stod(argv[2]); const double tolerance = std::stod(argv[3]); HEJ::istream in{argv[1]}; LHEF::Reader reader{in}; HEJ::PhaseSpacePointConfig psp_conf; psp_conf.jet_param = HEJ::JetParameters{JET_DEF, JETPTMIN}; psp_conf.max_ext_soft_pt_fraction = MAX_EXT_SOFT_PT_FRACTION; HEJ::MatrixElementConfig ME_conf; ME_conf.log_correction = LOG_CORR; ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{}; ME_conf.ew_parameters.set_vevWZH(VEV, WPROP, ZPROP, HPROP); HEJ::EventReweighterConfig conf; conf.psp_config = std::move(psp_conf); conf.ME_config = std::move(ME_conf); conf.treat = TREAT; reader.readEvent(); const bool has_Higgs = std::find( begin(reader.hepeup.IDUP), end(reader.hepeup.IDUP), 25 ) != end(reader.hepeup.IDUP); const double mu = has_Higgs?125.:91.188; HEJ::ScaleGenerator scale_gen{ {{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1. }; std::shared_ptr ran{std::make_shared()}; HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran}; HEJ::CrossSectionAccumulator xs; do{ auto ev_data = HEJ::Event::EventData{reader.hepeup}; shuffle_particles(ev_data); ev_data.reconstruct_intermediate(); HEJ::Event ev{ ev_data.cluster( BORN_JET_DEF, BORN_JETPTMIN ) }; auto resummed_events = hej.reweight(ev, NUM_TRIES); for(auto const & res_ev: resummed_events) { ASSERT(correct_colour(res_ev)); ASSERT(std::isfinite(res_ev.central().weight)); // we fill the xs uncorrelated since we only want to test the uncertainty // of the resummation xs.fill(res_ev); } } while(reader.readEvent()); const double xsec = xs.total().value; const double xsec_err = std::sqrt(xs.total().error); const double significance = std::abs(xsec - xsec_ref) / std::sqrt( xsec_err*xsec_err + tolerance*tolerance ); std::cout << xsec_ref << " +/- " << tolerance << " ~ " << xsec << " +- " << xsec_err << " => " << significance << " sigma\n"; if(significance > 3.){ std::cerr << "Cross section is off by over 3 sigma!\n"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/t/test_EWConstants.cc b/t/test_EWConstants.cc new file mode 100644 index 0000000..4b4b723 --- /dev/null +++ b/t/test_EWConstants.cc @@ -0,0 +1,51 @@ +/** + * \authors The HEJ collaboration (see AUTHORS for details) + * \date 2020 + * \copyright GPLv2 or later + */ +#include "hej_test.hh" + +#include "HEJ/EWConstants.hh" + +namespace { + constexpr double VEV = 38.; + constexpr HEJ::ParticleProperties WPROP{VEV, 11.}; + constexpr HEJ::ParticleProperties ZPROP{VEV*2., 13.}; + constexpr HEJ::ParticleProperties HPROP{19., 17.}; + constexpr double AWEAK = 2.*WPROP.mass/VEV*WPROP.mass/VEV; + constexpr double COS = WPROP.mass/ZPROP.mass; + inline bool operator==( + HEJ::ParticleProperties const & p1, HEJ::ParticleProperties const & p2 + ){ + return p1.mass==p2.mass && p1.width == p2.width; + } + inline bool operator==( + HEJ::EWConstants const & p1, HEJ::EWConstants const & p2 + ){ + return p1.vev()==p2.vev() + && p1.Wprop() == p2.Wprop() + && p1.Zprop() == p2.Zprop() + && p1.Hprop() == p2.Hprop(); + } +} + +int main() { + constexpr HEJ::EWConstants EWCONST{VEV, WPROP, ZPROP, HPROP}; + // check that everything is copied correctly + ASSERT(EWCONST.vev() == VEV); + ASSERT(EWCONST.prop(HEJ::ParticleID::Wp)==EWCONST.Wprop()); + ASSERT(EWCONST.prop(HEJ::ParticleID::Wm)==EWCONST.Wprop()); + ASSERT(EWCONST.prop(HEJ::ParticleID::Z) ==EWCONST.Zprop()); + ASSERT(EWCONST.prop(HEJ::ParticleID::h) ==EWCONST.Hprop()); + ASSERT(EWCONST.alpha_w() == AWEAK); + ASSERT(EWCONST.cos_tw() == COS); + ASSERT(EWCONST.cos2_tw()+EWCONST.sin2_tw() == 1.); + + // check that default is not set + HEJ::EWConstants ew2; + ASSERT_THROW(ew2.vev(), std::invalid_argument); + ew2.set_vevWZH(EWCONST.vev(), EWCONST.Wprop(), EWCONST.Zprop(), EWCONST.Hprop()); + ASSERT(ew2 == EWCONST); + + return EXIT_SUCCESS; +}