diff --git a/.clang-tidy b/.clang-tidy index 0dc056c..f96b147 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -1,38 +1,39 @@ Checks: ' *, -cert-dcl21-cpp, -clang-analyzer-optin.cplusplus.VirtualCall, -cppcoreguidelines-avoid-magic-numbers, -cppcoreguidelines-pro-bounds-array-to-pointer-decay, -cppcoreguidelines-pro-bounds-constant-array-index, -fuchsia*, -google-build-using-namespace, -google-explicit-constructor, -google-readability*, -google-runtime-int, -google-runtime-references, -hicpp*, -llvm-header-guard, -modernize-use-trailing-return-type, -readability-braces-around-statements, -readability-uppercase-literal-suffix, + -bugprone-macro-parentheses, ' CheckOptions: - key: cppcoreguidelines-special-member-functions.AllowSoleDefaultDtor value: '1' - key: llvm-namespace-comment.ShortNamespaceLines value: '10' - key: misc-non-private-member-variables-in-classes.IgnoreClassesWithAllMemberVariablesBeingPublic value: '1' - key: performance-move-const-arg.CheckTriviallyCopyableMove value: '0' - key: readability-identifier-naming.GlobalConstantCase value: UPPER_CASE - key: readability-identifier-naming.GlobalVariableCase value: UPPER_CASE - key: readability-identifier-naming.PrivateMemberSuffix value: '_' - key: readability-implicit-bool-conversion.AllowIntegerConditions value: '1' - key: readability-magic-numbers.IgnoreAllFloatingPointValues value: '1' diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc index 2e3fea7..3cdacfe 100644 --- a/FixedOrderGen/src/config.cc +++ b/FixedOrderGen/src/config.cc @@ -1,441 +1,419 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "config.hh" #include #include #include #include #include #include #include #include #include "HEJ/PDG_codes.hh" #include "HEJ/YAMLreader.hh" #include "HEJ/exceptions.hh" namespace HEJFOG { using HEJ::set_from_yaml; using HEJ::set_from_yaml_if_defined; using HEJ::pid::ParticleID; namespace { //! Get YAML tree of supported options /** * The configuration file is checked against this tree of options * in assert_all_options_known. */ YAML::Node const & get_supported_options(){ const static YAML::Node supported = [](){ YAML::Node supported; static const auto opts = { "process", "events", "subleading fraction","subleading channels", "scales", "scale factors", "max scale ratio", "pdf", "vev", "event output", "analyses", "analysis", "import scales" }; // add subnodes to "supported" - the assigned value is irrelevant for(auto && opt: opts) supported[opt] = ""; for(auto && jet_opt: {"min pt", "peak pt", "algorithm", "R", "max rapidity"}){ supported["jets"][jet_opt] = ""; } for(auto && particle_type: {"Higgs", "W", "Z"}){ for(auto && particle_opt: {"mass", "width"}){ supported["particle properties"][particle_type][particle_opt] = ""; } } for(auto && particle_type: {"Higgs", "Wp", "W+", "Wm", "W-", "Z"}){ for(auto && particle_opt: {"into", "branching ratio"}){ supported["decays"][particle_type][particle_opt] = ""; } } for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){ supported["Higgs coupling"][opt] = ""; } for(auto && beam_opt: {"energy", "particles"}){ supported["beam"][beam_opt] = ""; } for(auto && unweight_opt: {"sample size", "max deviation"}){ supported["unweight"][unweight_opt] = ""; } for(auto && opt: {"name", "seed"}){ supported["random generator"][opt] = ""; } return supported; }(); return supported; } JetParameters get_jet_parameters( YAML::Node const & node, std::string const & entry ){ const auto p = HEJ::get_jet_parameters(node, entry); JetParameters result; result.def = p.def; result.min_pt = p.min_pt; set_from_yaml(result.max_y, node, entry, "max rapidity"); set_from_yaml_if_defined(result.peak_pt, node, entry, "peak pt"); if(result.peak_pt && *result.peak_pt <= result.min_pt) throw std::invalid_argument{ "Value of option 'peak pt' has to be larger than 'min pt'." }; return result; } Beam get_Beam( YAML::Node const & node, std::string const & entry ){ Beam beam; std::vector particles; set_from_yaml(beam.energy, node, entry, "energy"); set_from_yaml_if_defined(particles, node, entry, "particles"); if(! particles.empty()){ for(HEJ::ParticleID particle: particles){ if(particle != HEJ::pid::p && particle != HEJ::pid::p_bar){ throw std::invalid_argument{ "Unsupported value in option " + entry + ": particles:" " only proton ('p') and antiproton ('p_bar') beams are supported" }; } } if(particles.size() != 2){ throw std::invalid_argument{"Not exactly two beam particles"}; } beam.particles.front() = particles.front(); beam.particles.back() = particles.back(); } return beam; } std::vector split( std::string const & str, std::string const & delims ){ std::vector result; for(size_t begin = 0, end = 0; end != std::string::npos;){ begin = str.find_first_not_of(delims, end); if(begin == std::string::npos) break; end = str.find_first_of(delims, begin + 1); result.emplace_back(str.substr(begin, end - begin)); } return result; } std::invalid_argument invalid_incoming(std::string const & what){ return std::invalid_argument{ "Incoming particle type " + what + " not supported," " incoming particles have to be 'p', 'p_bar' or partons" }; } std::invalid_argument invalid_outgoing(std::string const & what){ return std::invalid_argument{ "Outgoing particle type " + what + " not supported," " outgoing particles have to be 'j', 'photon', 'H', 'Wm', 'Wp', 'e-', 'e+', 'nu_e', 'nu_e_bar'" }; } HEJ::ParticleID reconstruct_boson_id( std::vector const & ids ){ assert(ids.size()==2); const int pidsum = ids[0] + ids[1]; if(pidsum == +1) { assert(HEJ::is_antilepton(ids[0])); if(HEJ::is_antineutrino(ids[0])) { throw HEJ::not_implemented{"lepton-flavour violating final state"}; } assert(HEJ::is_neutrino(ids[1])); // charged antilepton + neutrino means we had a W+ return HEJ::pid::Wp; } if(pidsum == -1) { assert(HEJ::is_antilepton(ids[0])); if(HEJ::is_neutrino(ids[1])) { throw HEJ::not_implemented{"lepton-flavour violating final state"}; } assert(HEJ::is_antineutrino(ids[0])); // charged lepton + antineutrino means we had a W- return HEJ::pid::Wm; } throw HEJ::not_implemented{ "final state with leptons "+name(ids[0])+" and "+name(ids[1]) +" not supported" }; } Process get_process( YAML::Node const & node, std::string const & entry ){ Process result; std::string process_string; set_from_yaml(process_string, node, entry); assert(! process_string.empty()); const auto particles = split(process_string, " \n\t\v=>"); if(particles.size() < 3){ throw std::invalid_argument{ "Bad format in option process: '" + process_string + "', expected format is 'in1 in2 => out1 ...'" }; } result.incoming.front() = HEJ::to_ParticleID(particles[0]); result.incoming.back() = HEJ::to_ParticleID(particles[1]); for(size_t i = 0; i < result.incoming.size(); ++i){ const HEJ::ParticleID in = result.incoming[i]; if( in != HEJ::pid::proton && in != HEJ::pid::p_bar && !HEJ::is_parton(in) ){ throw invalid_incoming(particles[i]); } } result.njets = 0; for(size_t i = result.incoming.size(); i < particles.size(); ++i){ assert(! particles[i].empty()); if(particles[i] == "j") ++result.njets; else if(std::isdigit(particles[i].front()) && particles[i].back() == 'j') result.njets += std::stoi(particles[i]); else{ const auto pid = HEJ::to_ParticleID(particles[i]); if(pid==HEJ::pid::Higgs || pid==HEJ::pid::Wp || pid==HEJ::pid::Wm){ if(result.boson) throw std::invalid_argument{ "More than one outgoing boson is not supported" }; if(!result.boson_decay.empty()) throw std::invalid_argument{ "Production of a boson together with a lepton is not supported" }; result.boson = pid; } else if (HEJ::is_anylepton(pid)){ // Do not accept more leptons, if two leptons are already mentioned if( result.boson_decay.size()>=2 ) throw std::invalid_argument{"Too many leptons provided"}; if(result.boson) throw std::invalid_argument{ "Production of a lepton together with a boson is not supported" }; result.boson_decay.emplace_back(pid); } else { throw invalid_outgoing(particles[i]); } } } if(result.njets < 2){ throw std::invalid_argument{ "Process has to include at least two jets ('j')" }; } if(!result.boson_decay.empty()){ std::sort(begin(result.boson_decay),end(result.boson_decay)); assert(!result.boson); result.boson = reconstruct_boson_id(result.boson_decay); } return result; } HEJFOG::Subleading to_subleading_channel(YAML::Node const & yaml){ std::string name; using namespace HEJFOG::subleading; set_from_yaml(name, yaml); if(name == "none") return NONE; if(name == "all") return ALL; Subleading channel = NONE; if(name == "unordered" || name == "uno") return channel.set(uno); if(name == "central qqx" || name == "cqqx") return channel.set(cqqx); if(name == "extremal qqx" || name == "eqqx") return channel.set(eqqx); throw HEJ::unknown_option("Unknown subleading channel '"+name+"'"); } Subleading get_subleading_channels(YAML::Node const & node){ using YAML::NodeType; using namespace HEJFOG::subleading; // all channels allowed by default if(!node) return ALL; switch(node.Type()){ case NodeType::Undefined: return ALL; case NodeType::Null: return NONE; case NodeType::Scalar: return to_subleading_channel(node); case NodeType::Map: throw HEJ::invalid_type{"map is not a valid option for subleading channels"}; case NodeType::Sequence: Subleading channels; for(auto && channel_node: node){ channels |= get_subleading_channels(channel_node); } return channels; } throw std::logic_error{"unreachable"}; } Decay get_decay(YAML::Node const & node, std::string const & entry, std::string const & boson ){ Decay decay; set_from_yaml(decay.products, node, entry, boson, "into"); decay.branching_ratio=1; set_from_yaml_if_defined(decay.branching_ratio, node, entry, boson, "branching ratio"); return decay; } std::vector get_decays(YAML::Node const & node, std::string const & entry, std::string const & boson ){ using YAML::NodeType; if(!node[entry][boson]) return {}; switch(node[entry][boson].Type()){ case NodeType::Null: case NodeType::Undefined: return {}; case NodeType::Scalar: throw HEJ::invalid_type{"value is not a list of decays"}; case NodeType::Map: return {get_decay(node, entry, boson)}; case NodeType::Sequence: std::vector result; for(auto && decay_str: node[entry][boson]){ result.emplace_back(get_decay(decay_str, entry, boson)); } return result; } throw std::logic_error{"unreachable"}; } ParticlesDecayMap get_all_decays(YAML::Node const & node, std::string const & entry ){ if(!node[entry]) return {}; if(!node[entry].IsMap()) throw HEJ::invalid_type{entry + " have to be a map"}; ParticlesDecayMap result; for(auto const & sub_node: node[entry]) { const auto boson = sub_node.first.as(); const auto id = HEJ::to_ParticleID(boson); result.emplace(id, get_decays(node, entry, boson)); } return result; } - HEJ::ParticleProperties get_particle_properties( - YAML::Node const & node, std::string const & entry, - std::string const & boson - ){ - HEJ::ParticleProperties result{}; - set_from_yaml(result.mass, node, entry, boson, "mass"); - set_from_yaml(result.width, node, entry, boson, "width"); - return result; - } - - HEJ::EWConstants get_ew_parameters(YAML::Node const & node){ - HEJ::EWConstants result; - double vev = NAN; - set_from_yaml(vev, node, "vev"); - result.set_vevWZH(vev, - get_particle_properties(node, "particle properties", "W"), - get_particle_properties(node, "particle properties", "Z"), - get_particle_properties(node, "particle properties", "Higgs") - ); - return result; - } - UnweightSettings get_unweight( YAML::Node const & node, std::string const & entry ){ UnweightSettings result{}; set_from_yaml(result.sample_size, node, entry, "sample size"); if(result.sample_size <= 0){ throw std::invalid_argument{ "negative sample size " + std::to_string(result.sample_size) }; } set_from_yaml(result.max_dev, node, entry, "max deviation"); return result; } Config to_Config(YAML::Node const & yaml){ try{ HEJ::assert_all_options_known(yaml, get_supported_options()); } catch(HEJ::unknown_option const & ex){ throw HEJ::unknown_option{std::string{"Unknown option '"} + ex.what() + "'"}; } Config config; config.process = get_process(yaml, "process"); set_from_yaml(config.events, yaml, "events"); config.jets = get_jet_parameters(yaml, "jets"); config.beam = get_Beam(yaml, "beam"); for(size_t i = 0; i < config.process.incoming.size(); ++i){ auto const & in = config.process.incoming[i]; using namespace HEJ::pid; if( (in == p || in == p_bar) && in != config.beam.particles[i]){ throw std::invalid_argument{ "Particle type of beam " + std::to_string(i+1) + " incompatible" + " with type of incoming particle " + std::to_string(i+1) }; } } set_from_yaml(config.pdf_id, yaml, "pdf"); set_from_yaml(config.subleading_fraction, yaml, "subleading fraction"); if(config.subleading_fraction == 0) config.subleading_channels.reset(); else config.subleading_channels = get_subleading_channels(yaml["subleading channels"]); - config.ew_parameters = get_ew_parameters(yaml); + config.ew_parameters = HEJ::get_ew_parameters(yaml); config.particle_decays = get_all_decays(yaml, "decays"); if(config.process.boson // check that Ws always decay && std::abs(*config.process.boson) == HEJ::ParticleID::Wp && config.process.boson_decay.empty() ){ auto const & decay = config.particle_decays.find(*config.process.boson); if(decay == config.particle_decays.end() || decay->second.empty()) throw std::invalid_argument{ "Decay for "+name(*config.process.boson)+" is required"}; } set_from_yaml_if_defined(config.analyses_parameters, yaml, "analyses"); if(yaml["analysis"].IsDefined()){ std::cerr << "WARNING: Configuration entry 'analysis' is deprecated. " " Use 'analyses' instead.\n"; YAML::Node ana; set_from_yaml(ana, yaml, "analysis"); if(!ana.IsNull()){ config.analyses_parameters.push_back(ana); } } config.scales = HEJ::to_ScaleConfig(yaml); set_from_yaml_if_defined(config.output, yaml, "event output"); config.rng = HEJ::to_RNGConfig(yaml, "random generator"); config.Higgs_coupling = HEJ::get_Higgs_coupling(yaml, "Higgs coupling"); if(yaml["unweight"].IsDefined()) config.unweight = get_unweight(yaml, "unweight"); return config; } } // namespace Config load_config(std::string const & config_file){ try{ return to_Config(YAML::LoadFile(config_file)); } catch(...){ std::cerr << "Error reading " << config_file << ":\n "; throw; } } } // namespace HEJFOG 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/include/HEJ/PDF.hh b/include/HEJ/PDF.hh index 44a09b2..a1c2e39 100644 --- a/include/HEJ/PDF.hh +++ b/include/HEJ/PDF.hh @@ -1,81 +1,81 @@ /** \file * * \brief Contains all the necessary classes and functions for interaction with PDFs. * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "HEJ/PDG_codes.hh" namespace LHAPDF { class PDF; } namespace HEJ { //! Class for interaction with a PDF set class PDF { public: /** * \brief PDF Constructor * @param id Particle ID according to PDG * @param beam1 Particle ID of particle in beam 1 * @param beam2 Particle ID of particle in beam 2 */ PDF(int id, ParticleID beam1, ParticleID beam2); /** * \brief Calculate the pdf value x*f(x, q) * @param beam_idx Beam number (0 or 1) * @param x Momentum fraction * @param q Energy scale * @param id PDG particle id * @returns x*f(x, q) * * Returns 0 if x or q are outside the range covered by the PDF set */ double pdfpt(std::size_t beam_idx, double x, double q, ParticleID id) const; /** * \brief Value of the strong coupling \f$\alpha_s(q)\f$ at the given scale * @param q Renormalisation scale * @returns Value of the strong coupling constant */ double Halphas(double q) const; //! Check if the energy scale is within the range covered by the PDF set /** * @param q Energy Scale * @returns true if q is within the covered range, false otherwise */ bool inRangeQ(double q) const; //! Check if the momentum fraction is within the range covered by the PDF set /** * @param x Momentum Fraction * @returns true if x is within the covered range, false otherwise */ bool inRangeX(double x) const; //! PDF id of the current set int id() const; PDF(PDF const & other) = delete; //!< Copy not possible PDF & operator=(PDF const & other) = delete; //!< Copy not possible - PDF(PDF && other); //!< Moving allowed - PDF & operator=(PDF && other); //!< Moving allowed + PDF(PDF && other) noexcept; //!< Moving allowed + PDF & operator=(PDF && other) noexcept; //!< Moving allowed ~PDF(); private: //! @internal unique_ptr does not allow copy std::unique_ptr pdf_; - std::array beamtype_; + std::array beamtype_{}; }; } // namespace HEJ diff --git a/include/HEJ/Particle.hh b/include/HEJ/Particle.hh index d8f0aa2..21cf292 100644 --- a/include/HEJ/Particle.hh +++ b/include/HEJ/Particle.hh @@ -1,225 +1,227 @@ /** * \file Particle.hh * \brief Contains the particle struct * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include "fastjet/PseudoJet.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/optional.hh" namespace HEJ { using Colour = std::pair; //! Class representing a particle struct Particle { //! particle type ParticleID type = pid::unspecified; //! particle momentum fastjet::PseudoJet p; //! (optional) colour & anti-colour optional colour; //! get rapidity double rapidity() const{ return p.rapidity(); } //! get transverse momentum double perp() const{ return p.perp(); } //! get transverse momentum double pt() const{ return perp(); } //! get momentum in x direction double px() const{ return p.px(); } //! get momentum in y direction double py() const{ return p.py(); } //! get momentum in z direction double pz() const{ return p.pz(); } //! get energy double E() const{ return p.E(); } //! get mass double m() const{ return p.m(); } }; //! Functor to compare rapidities /** * This can be used whenever a rapidity comparison function is needed, * for example in many standard library functions. * * @see pz_less */ struct rapidity_less{ template bool operator()(FourVector const & p1, FourVector const & p2){ return p1.rapidity() < p2.rapidity(); } }; //! Functor to compare momenta in z direction /** * This can be used whenever a pz comparison function is needed, * for example in many standard library functions. * * @see rapidity_less */ struct pz_less{ template bool operator()(FourVector const & p1, FourVector const & p2){ return p1.pz() < p2.pz(); } }; //! Convert a vector of Particles to a vector of particle momenta inline std::vector to_PseudoJet( std::vector const & v ){ std::vector result; result.reserve(v.size()); for(auto && sp: v) result.emplace_back(sp.p); return result; } //! Check if a particle is a parton, i.e. quark, antiquark, or gluon inline - bool is_parton(Particle const & p){ + constexpr bool is_parton(Particle const & p){ return is_parton(p.type); } //! Check if a particle is a quark inline - bool is_quark(Particle const & p){ + constexpr bool is_quark(Particle const & p){ return is_quark(p.type); } //! Check if a particle is an anti-quark inline - bool is_antiquark(Particle const & p){ + constexpr bool is_antiquark(Particle const & p){ return is_antiquark(p.type); } //! Check if a particle is a quark or anit-quark inline - bool is_anyquark(Particle const & p){ + constexpr bool is_anyquark(Particle const & p){ return is_anyquark(p.type); } /** * \brief Function to determine if particle is a lepton * @param p the particle * @returns true if the particle is a lepton, false otherwise */ inline constexpr bool is_lepton(Particle const & p){ return is_lepton(p.type); } /** * \brief Function to determine if particle is an antilepton * @param p the particle * @returns true if the particle is an antilepton, false otherwise */ inline constexpr bool is_antilepton(Particle const & p){ return is_antilepton(p.type); } /** * \brief Function to determine if particle is an (anti-)lepton * @param p the particle * @returns true if the particle is a lepton or antilepton, false otherwise */ inline constexpr bool is_anylepton(Particle const & p){ return is_anylepton(p.type); } /** * \brief Function to determine if particle is a neutrino * @param p the particle * @returns true if the particle is a neutrino, false otherwise */ inline constexpr bool is_neutrino(Particle const & p){ return is_neutrino(p.type); } /** * \brief Function to determine if particle is an antineutrino * @param p the particle * @returns true if the particle is an antineutrino, false otherwise */ inline constexpr bool is_antineutrino(Particle const & p){ return is_antineutrino(p.type); } /** * \brief Function to determine if particle is an (anti-)neutrino * @param p the particle * @returns true if the particle is a neutrino or antineutrino, false otherwise */ inline constexpr bool is_anyneutrino(Particle const & p){ return is_anyneutrino(p.type); } //! Check if a particle is a photon, W or Z boson - inline bool is_AWZ_boson(Particle const & particle){ + inline + constexpr bool is_AWZ_boson(Particle const & particle){ return is_AWZ_boson(particle.type); } //! Check if a particle is a photon, W, Z, or Higgs boson - inline bool is_AWZH_boson(Particle const & particle){ + inline + constexpr bool is_AWZH_boson(Particle const & particle){ return is_AWZH_boson(particle.type); } //! Extract all partons from a vector of particles inline std::vector filter_partons( std::vector const & v ){ std::vector result; result.reserve(v.size()); std::copy_if( begin(v), end(v), std::back_inserter(result), [](Particle const & p){ return is_parton(p); } ); return result; } //! Extract all AWZH bosons from a vector of particles inline std::vector filter_AWZH_bosons( std::vector const & v ){ std::vector result; std::copy_if( begin(v), end(v), std::back_inserter(result), [](Particle const & p){ return is_AWZH_boson(p); } ); return result; } } // namespace HEJ diff --git a/include/HEJ/event_types.hh b/include/HEJ/event_types.hh index 9723543..b53ba36 100644 --- a/include/HEJ/event_types.hh +++ b/include/HEJ/event_types.hh @@ -1,108 +1,109 @@ /** \file * \brief Define different types of events. * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include "HEJ/exceptions.hh" namespace HEJ { //! Namespace for event types namespace event_type { //! Possible event types enum EventType: std::size_t { non_resummable = 0, //!< event configuration not covered by All Order resummation bad_final_state = 1, //!< event with an unsupported final state no_2_jets = 2, //!< event with less than two jets FKL = 4, //!< FKL-type event unordered_backward = 8, //!< event with unordered backward emission unordered_forward = 16, //!< event with unordered forward emission extremal_qqxb = 32, //!< event with a backward extremal qqbar extremal_qqxf = 64, //!< event with a forward extremal qqbar central_qqx = 128, //!< event with a central qqbar unob = unordered_backward, //!< alias for unordered_backward unof = unordered_forward, //!< alias for unordered_forward qqxexb = extremal_qqxb, //!< alias for extremal_qqxb qqxexf = extremal_qqxf, //!< alias for extremal_qqxf qqxmid = central_qqx, //!< alias for central_qqx first_type = non_resummable, //!< alias for non_resummable last_type = central_qqx //!< alias for central_qqx }; //! Event type names /** * For example, name(FKL) is the string "FKL" */ - inline std::string name(EventType type) { + inline + std::string name(EventType type) { switch(type) { case FKL: return "FKL"; case unordered_backward: return "unordered backward"; case unordered_forward: return "unordered forward"; case extremal_qqxb: return "extremal qqbar backward"; case extremal_qqxf: return "extremal qqbar forward"; case central_qqx: return "central qqbar"; case non_resummable: return "non-resummable"; case no_2_jets: return "no 2 jets"; case bad_final_state: return "bad final state"; default: throw std::logic_error{"Unreachable"}; } } //! Returns True for a HEJ \ref event_type::EventType "EventType" inline - bool is_resummable(EventType type) { + constexpr bool is_resummable(EventType type) { switch(type) { case FKL: case unordered_backward: case unordered_forward: case extremal_qqxb: case extremal_qqxf: case central_qqx: return true; default: return false; } } //! Returns True for an unordered \ref event_type::EventType "EventType" inline - bool is_uno(EventType type) { + constexpr bool is_uno(EventType type) { return type == unordered_backward || type == unordered_forward; } //! Returns True for an extremal_qqx \ref event_type::EventType "EventType" inline - bool is_ex_qqx(EventType type) { + constexpr bool is_ex_qqx(EventType type) { return type == extremal_qqxb || type == extremal_qqxf; } //! Returns True for an central_qqx \ref event_type::EventType "EventType" inline - bool is_mid_qqx(EventType type) { + constexpr bool is_mid_qqx(EventType type) { return type == central_qqx; } //! Returns True for any qqx event \ref event_type::EventType "EventType" inline - bool is_qqx(EventType type) { + constexpr bool is_qqx(EventType type) { return is_ex_qqx(type) || is_mid_qqx(type); } } // namespace event_type } // namespace HEJ diff --git a/include/HEJ/utility.hh b/include/HEJ/utility.hh index 6f94205..997dce8 100644 --- a/include/HEJ/utility.hh +++ b/include/HEJ/utility.hh @@ -1,103 +1,104 @@ /** * \file * \brief Contains various utilities * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "boost/core/demangle.hpp" #include "fastjet/PseudoJet.hh" namespace HEJ { inline std::string join( std::string const & /* delim */ ){ return ""; } inline std::string join( std::string const & /* delim */, std::string const & str ){ return str; } //! Join strings with a delimiter /** * @param delim Delimiter to be put between consecutive strings * @param first First string * @param second Second string * @param rest Remaining strings */ template std::string join( std::string const & delim, std::string const & first, std::string const & second, Strings&&... rest ){ return join(delim, first + delim + second, std::forward(rest)...); } //! Return the name of the argument's type template std::string type_string(T&& /*unused*/){ return boost::core::demangle(typeid(T).name()); } //! Eliminate compiler warnings for unused variables template constexpr void ignore(T&&... /*unused*/) {} //! Check whether two doubles are closer than ep > 0 to each other inline - bool nearby_ep(double a, double b, double ep){ + constexpr bool nearby_ep(double a, double b, double ep){ assert(ep > 0); return std::abs(a-b) < ep; } //! Check whether all components of two PseudoJets are closer than ep to each other inline bool nearby_ep( fastjet::PseudoJet const & pa, fastjet::PseudoJet const & pb, double ep ){ assert(ep > 0); for(size_t i = 0; i < 4; ++i){ if(!nearby_ep(pa[i], pb[i], ep)) return false; } return true; } inline bool nearby( - fastjet::PseudoJet const & pa, fastjet::PseudoJet const & pb, double const norm = 1. + fastjet::PseudoJet const & pa, fastjet::PseudoJet const & pb, + double const norm = 1. ){ return nearby_ep(pa, pb, 1e-7*norm); } namespace detail { template struct ArrayTag{ using type = typename ArrayTag, Ns...>::type; }; template struct ArrayTag { using type = std::array; }; } // helper for multidimensional std::array, for example // MultiArray = std::array, N2> template using MultiArray = typename detail::ArrayTag::type; } // namespace HEJ diff --git a/src/PDF.cc b/src/PDF.cc index 4cdd837..34f2e10 100644 --- a/src/PDF.cc +++ b/src/PDF.cc @@ -1,76 +1,76 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/PDF.hh" #include #include #include #include #include #include "LHAPDF/LHAPDF.h" namespace HEJ { namespace { int to_beam(ParticleID id){ if(std::abs(id) == pid::proton){ return (id > 0)?1:-1; } throw std::invalid_argument( "unknown beam type: " + std::to_string(id) ); } } // namespace PDF::PDF(int id, ParticleID beam1, ParticleID beam2): pdf_{LHAPDF::mkPDF(id)}, beamtype_{{to_beam(beam1), to_beam(beam2)}} {} double PDF::pdfpt(size_t beam_idx, double x, double q, ParticleID id) const{ if(!(inRangeQ(q) && inRangeX(x))) return 0.; if(id == pid::gluon){ return pdf_->xfxQ(pid::gluon,x,q); } if(std::abs(id) <= pid::top){ return pdf_->xfxQ(id*beamtype_[beam_idx],x,q); } std::cerr << "particle type unknown: "<< id << std::endl; return 0.0; } double PDF::Halphas(double q) const{ double as = pdf_->alphasQ(q); if (std::isnan(as) || as > 0.5) { as = 0.5; } return as; } int PDF::id() const{ return pdf_->lhapdfID(); } bool PDF::inRangeQ(double q) const{ return pdf_->inRangeQ(q); } bool PDF::inRangeX(double x) const{ return pdf_->inRangeX(x); } // can't be in header since we forward declare LHAPDF::PDF PDF::~PDF() = default; - PDF::PDF(PDF && other) = default; - PDF & PDF::operator=(PDF && other) = default; + PDF::PDF(PDF && other) noexcept = default; + PDF & PDF::operator=(PDF && other) noexcept = default; } // namespace HEJ diff --git a/src/jets.cc b/src/jets.cc index 265ea9b..88efa02 100644 --- a/src/jets.cc +++ b/src/jets.cc @@ -1,357 +1,355 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/jets.hh" #include #include #include #include "HEJ/Constants.hh" // generated headers #include "HEJ/currents/j_j.hh" #include "HEJ/currents/j_jqqbar.hh" #include "HEJ/currents/j_qqbar_j.hh" #include "HEJ/currents/juno_j.hh" namespace { // short hand for math functions - using std::abs; - using std::conj; - using std::sqrt; -} // Anonymous Namespace + using std::conj; + } // Anonymous Namespace namespace HEJ { namespace currents { // Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113 // @TODO: this is not a current and should be moved somewhere else double K_g(double p1minus, double paminus) { return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A; } double K_g( HLV const & pout, HLV const & pin ) { if(pin.z() > 0) return K_g(pout.plus(), pin.plus()); return K_g(pout.minus(), pin.minus()); } namespace { //@{ /** * @brief Pure Jet FKL Contributions, function to handle all incoming types. * @param p1out Outgoing Particle 1. * @param p1in Incoming Particle 1. * @param p2out Outgoing Particle 2 * @param p2in Incoming Particle 2 * * Calculates j_\mu j^\mu. * Handles all possible incoming states. Helicity doesn't matter since we sum * over all of them. In addition, we only have to distinguish between the two * possibilities of contracting same-helicity or opposite-helicity currents. */ double j_j(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ using helicity::plus; using helicity::minus; COM const Mp = HEJ::j_j(p1in, p1out, p2in, p2out); COM const Mm = HEJ::j_j(p1in, p1out, p2in, p2out); double const sst = std::norm(Mm) + std::norm(Mp); HLV const q1=p1in-p1out; HLV const q2=-(p2in-p2out); // Multiply by Cf^2 (colour) * 2 (helicities) return 2.*C_F*C_F*(sst)/(q1.m2()*q2.m2()); } } // Anonymous Namespace double ME_qQ(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return j_j(p1out, p1in, p2out, p2in); } double ME_qQbar(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return j_j(p1out, p1in, p2out, p2in); } double ME_qbarQbar(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return j_j(p1out, p1in, p2out, p2in); } double ME_qg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/C_F; } double ME_qbarg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/(C_F); } double ME_gg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return j_j(p1out, p1in, p2out, p2in)*K_g(p1out, p1in)*K_g(p2out, p2in)/(C_F*C_F); } //@} namespace { template double amp_juno_j( HLV const & pa, HLV const & pb, HLV const & pg, HLV const & p1, HLV const & p2 ) { // TODO: code duplication with Wjets const COM u1 = U1_j(pa,p1,pb,p2,pg); const COM u2 = U2_j(pa,p1,pb,p2,pg); const COM l = L_j(pa,p1,pb,p2,pg); const COM x = u1 - l; const COM y = u2 + l; return C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real(); } double juno_j(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ using helicity::minus; using helicity::plus; const HLV q2 = p2out - p2in; // Bottom End const HLV qg = p1in - p1out - pg; // Extra bit post-gluon // only calculate half of the helicity amplitudes, // the remaining ones follow from CP symmetry const double amm = amp_juno_j(p1in, p2in, pg, p1out, p2out); const double amp = amp_juno_j(p1in, p2in, pg, p1out, p2out); const double apm = amp_juno_j< plus, minus>(p1in, p2in, pg, p1out, p2out); const double app = amp_juno_j< plus, plus>(p1in, p2in, pg, p1out, p2out); constexpr double hel_fac = 2.; return hel_fac/(4.*C_A*C_A)*(amm + amp + apm + app)/(q2.m2()*qg.m2()); } } // Anonymous Namespace //Unordered bits for pure jet double ME_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return juno_j(pg, p1out, p1in, p2out, p2in); } double ME_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return juno_j(pg, p1out, p1in, p2out, p2in); } double ME_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return juno_j(pg, p1out, p1in, p2out, p2in); } double ME_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return juno_j(pg, p1out, p1in, p2out, p2in); } double ME_unob_qg( HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F; } double ME_unob_qbarg(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F; } namespace { // helicity amplitudes for j jqqbar contraction template double amp_j_jqqbar( HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2, HLV const & p3 ) { // TODO: code duplication with Wjets.cc using std::norm; static constexpr double cm1m1 = 8./3.; static constexpr double cm2m2 = 8./3.; static constexpr double cm3m3 = 6.; static constexpr double cm1m2 =-1./3.; static constexpr COM cm1m3 = COM{0.,-3.}; static constexpr COM cm2m3 = COM{0.,3.}; const COM m1 = j_qggm1(pb,p2,p3,pa,p1); const COM m2 = j_qggm2(pb,p2,p3,pa,p1); const COM m3 = j_qggm3(pb,p2,p3,pa,p1); return + cm1m1*norm(m1) + cm2m2*norm(m2) + cm3m3*norm(m3) + 2.*real(cm1m2*m1*conj(m2)) + 2.*real(cm1m3*m1*conj(m3)) + 2.*real(cm2m3*m2*conj(m3)) ; } //Now the function to give helicity/colour sum/average double MqgtqQQ(HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2, HLV const & p3 ){ using helicity::minus; using helicity::plus; const double Mmmm = amp_j_jqqbar(pa, pb, p1, p2, p3); const double Mmmp = amp_j_jqqbar(pa, pb, p1, p2, p3); const double Mpmm = amp_j_jqqbar< plus, minus>(pa, pb, p1, p2, p3); const double Mpmp = amp_j_jqqbar< plus, plus>(pa, pb, p1, p2, p3); // Factor of 2 for the helicity for conjugate configurations return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/3.)/(pa-p1).m2()/(p2+p3-pb).m2(); } } // Anonymous Namespace // Extremal qqx double ME_Exqqx_qbarqQ(HLV const & pgin, HLV const & pqout, HLV const & pqbarout, HLV const & p2out, HLV const & p2in ){ return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout); } double ME_Exqqx_qqbarQ(HLV const & pgin, HLV const & pqout, HLV const & pqbarout, HLV const & p2out, HLV const & p2in ){ return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout); } double ME_Exqqx_qbarqg(HLV const & pgin, HLV const & pqout, HLV const & pqbarout, HLV const & p2out, HLV const & p2in ){ return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout)*K_g(p2out,p2in)/C_F; } double ME_Exqqx_qqbarg(HLV const & pgin, HLV const & pqout, HLV const & pqbarout, HLV const & p2out, HLV const & p2in ){ return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout)*K_g(p2out,p2in)/C_F; } namespace { // helicity amplitudes for matrix element with central qqbar template double amp_Cenqqx_qq( HLV const & pa, HLV const & p1, HLV const & pb, HLV const & p4, HLV const & pq, HLV const & pqbar, HLV const & q11, HLV const & q12 ){ using std::norm; const COM sym = M_sym( pa, p1, pb, p4, pq, pqbar, q11, q12 ); const COM cross = M_cross( pa, p1, pb, p4, pq, pqbar, q11, q12 ); const COM uncross = M_qbar( pa, p1, pb, p4, pq, pqbar, q11, q12 ); // Colour factors static constexpr double cmsms = 3.; static constexpr double cmumu = 4./3.; static constexpr double cmcmc = 4./3.; static constexpr COM cmsmu = COM{0., 3./2.}; static constexpr COM cmsmc = COM{0.,-3./2.}; static constexpr double cmumc = -1./6.; return +cmsms*norm(sym) +cmumu*norm(uncross) +cmcmc*norm(cross) +2.*real(cmsmu*sym*conj(uncross)) +2.*real(cmsmc*sym*conj(cross)) +2.*real(cmumc*uncross*conj(cross)) ; } } // Anonymous Namespace double ME_Cenqqx_qq( HLV const & pa, HLV const & pb, std::vector const & partons, bool /* aqlinepa */, bool /* aqlinepb */, const bool qqxmarker, const std::size_t nabove ){ using helicity::plus; using helicity::minus; CLHEP::HepLorentzVector q1 = pa; for(std::size_t i = 0; i <= nabove; ++i){ q1 -= partons[i]; } const auto qq = split_into_lightlike(q1); // since q1.m2() < 0 the following assertion is always true // see documentation for split_into_lightlike assert(qq.second.e() < 0); HLV pq = partons[nabove+1]; HLV pqbar = partons[nabove+2]; if(qqxmarker) std::swap(pq, pqbar); const HLV p1 = partons.front(); const HLV pn = partons.back(); // 8 helicity choices, but only 4 independent ones //(complex conjugation related). // In principle, the proper helicity labeling depends on // whether we have antiquark lines (aqlinea and aqlineb). // However, this only corresponds to a relabeling, // which doesn't change the sum over all helicities below. const double amp_mm = amp_Cenqqx_qq( pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second ); const double amp_mp = amp_Cenqqx_qq( pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second ); const double amp_pm = amp_Cenqqx_qq( pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second ); const double amp_pp = amp_Cenqqx_qq( pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second ); //Result (averaged, without coupling or t-channel props). Factor of //2 for the 4 helicity configurations I didn't work out explicitly const HLV q3 = q1 - pq - pqbar; return (2.*(amp_mm+amp_mp+amp_pm+amp_pp)/9./4.) / ((pa-p1).m2()*(pb-pn).m2()*q1.m2()*q3.m2()); } } // namespace currents } // 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; +}