diff --git a/FixedOrderGen/include/Subleading.hh b/FixedOrderGen/include/Subleading.hh index 9f785e1..69093c4 100644 --- a/FixedOrderGen/include/Subleading.hh +++ b/FixedOrderGen/include/Subleading.hh @@ -1,34 +1,36 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019 * \copyright GPLv2 or later */ #pragma once #include +#include namespace HEJFOG { namespace subleading { //!< @TODO confusing name with capital Subleading /** * Bit position of different subleading channels * e.g. (unsigned int) 1 => only unordered */ enum Channels: unsigned { uno, unordered = uno, cqqbar, central_qqbar = cqqbar, eqqbar, extremal_qqbar = eqqbar, first = uno, last = eqqbar, }; } // namespace subleading + std::string_view name(subleading::Channels c); using Subleading = std::bitset; namespace subleading { static constexpr Subleading ALL{~0u}; static constexpr Subleading NONE{0u}; } } // namespace HEJFOG diff --git a/FixedOrderGen/include/Utility.hh b/FixedOrderGen/include/Utility.hh new file mode 100644 index 0000000..9a4439e --- /dev/null +++ b/FixedOrderGen/include/Utility.hh @@ -0,0 +1,29 @@ +/** + * \authors The HEJ collaboration (see AUTHORS for details) + * \date 2022 + * \copyright GPLv2 or later + */ +#pragma once + +#include + +namespace HEJFOG { + class Process; + + bool vector_boson_can_couple_to( + HEJ::pid::ParticleID vector_boson, + HEJ::pid::ParticleID particle + ); + + bool parton_can_couple_to_W( + HEJ::Particle const & part, + HEJ::pid::ParticleID W_id + ); + + bool parton_can_couple_to_W( + HEJ::pid::ParticleID part, + HEJ::pid::ParticleID W_id + ); + + bool is_AWZ_process(Process const & proc); +} diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc index 796c75f..3b9bde9 100644 --- a/FixedOrderGen/src/EventGenerator.cc +++ b/FixedOrderGen/src/EventGenerator.cc @@ -1,127 +1,246 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "EventGenerator.hh" +#include #include +#include +#include #include #include "HEJ/Config.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Event.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/event_types.hh" #include "HEJ/exceptions.hh" #include "HEJ/make_RNG.hh" +#include "HEJ/PDG_codes.hh" #include "Beam.hh" #include "JetParameters.hh" #include "PhaseSpacePoint.hh" #include "Process.hh" +#include "Subleading.hh" +#include "Utility.hh" #include "config.hh" namespace HEJFOG { + namespace { + + bool allows_gluon(const HEJ::ParticleID particle_id) { + return + std::abs(particle_id) == HEJ::pid::proton + || particle_id == HEJ::pid::gluon; + } + + bool allows_anyquark(const HEJ::ParticleID particle_id) { + return + std::abs(particle_id) == HEJ::pid::proton + || HEJ::is_anyquark(particle_id); + } + + Subleading possible_subl_channels(Process const &process) { + std::size_t nhiggs_bosons = + (process.boson && *process.boson == HEJ::pid::Higgs) ? 1 : 0; + if(process.njets + nhiggs_bosons <= 2){ + return subleading::NONE; + } + auto allowed = subleading::ALL; + if(process.njets + nhiggs_bosons < 4){ + allowed.reset(subleading::central_qqbar); + } + if( + !allows_gluon(process.incoming[0]) + && !allows_gluon(process.incoming[1]) + ){ + allowed.reset(subleading::extremal_qqbar); + } + + if( + !allows_anyquark(process.incoming[0]) + && !allows_anyquark(process.incoming[1]) + ){ + allowed.reset(subleading::unordered); + } + + if( + process.boson && std::abs(*process.boson) == HEJ::pid::Wp + && !vector_boson_can_couple_to(process.incoming[0], *process.boson) + && !vector_boson_can_couple_to(process.incoming[1], *process.boson) + ){ + allowed.reset(subleading::unordered); + } + return allowed; + } + + bool can_generate_FKL(Process const &process) { + std::size_t nhiggs_bosons = + (process.boson && *process.boson == HEJ::pid::Higgs) ? 1 : 0; + if(process.njets + nhiggs_bosons < 2){ + return false; + } + if( + is_AWZ_process(process) + && !vector_boson_can_couple_to(process.incoming[0], *process.boson) + && !vector_boson_can_couple_to(process.incoming[1], *process.boson) + ) return false; + return true; + } + + // TODO: maybe we don't even want to accept inconsistent settings here + void ensure_valid_subl( + Process const &process, + double &subl_chance, + Subleading &subl_channels + ){ + if(subl_chance == 0. && subl_channels.any()) { + std::cerr << "WARNING: Subleading fraction is set to 0," + " disabling all subleading channels\n"; + subl_channels = subleading::NONE; + } + const Subleading possible_channels = possible_subl_channels(process); + for ( + unsigned channel = subleading::Channels::first; + channel <= subleading::Channels::last; + ++channel + ){ + if(subl_channels.test(channel) && !possible_channels.test(channel)){ + std::cerr + << "WARNING: selected subleading channel " + << name(static_cast(channel)) + << " cannot be generated for the chosen process.\n"; + subl_channels.reset(channel); + } + } + if(subl_channels == subleading::NONE && subl_chance > 0.) { + if(subl_chance == 1.) { + throw std::invalid_argument{ + "Cannot generate subleading event configurations" + }; + } + std::cerr << + "WARNING: Positive fraction of subleading configurations requested," + " but cannot generate any.\n"; + subl_chance = 0.; + } + if(!can_generate_FKL(process) && subl_chance < 1.) { + if(subl_chance > 0.) { + std::cerr << + "WARNING: Can only generate subleading configurations.\n"; + subl_chance = 1.; + } else { + throw std::invalid_argument{ + "Cannot generate events for given process" + }; + } + } + } + } + EventGenerator::EventGenerator(Config const & config): EventGenerator{ config.process, config.beam, HEJ::ScaleGenerator{ config.scales.base, config.scales.factors, config.scales.max_ratio }, config.jets, config.pdf_id, config.subleading_fraction, config.subleading_channels, config.particle_decays, config.Higgs_coupling, config.ew_parameters, config.nlo, HEJ::make_RNG(config.rng.name, config.rng.seed) } {} EventGenerator::EventGenerator( Process process, Beam beam, HEJ::ScaleGenerator scale_gen, JetParameters jets, int pdf_id, double subl_chance, Subleading subl_channels, ParticlesDecayMap particle_decays, HEJ::HiggsCouplingSettings Higgs_coupling, HEJ::EWConstants ew_parameters, HEJ::NLOConfig nlo, std::shared_ptr ran ): pdf_{pdf_id, beam.particles[0], beam.particles[1]}, ME_{ [this](double mu){ return pdf_.Halphas(mu); }, HEJ::MatrixElementConfig{ false, std::move(Higgs_coupling), ew_parameters, nlo } }, scale_gen_{std::move(scale_gen)}, process_{std::move(process)}, jets_{std::move(jets)}, beam_{std::move(beam)}, status_{Status::good}, subl_chance_{subl_chance}, subl_channels_{subl_channels}, particle_decays_{std::move(particle_decays)}, ew_parameters_{ew_parameters}, ran_{std::move(ran)} { + ensure_valid_subl(process_, subl_chance_, subl_channels_); assert(ran_ != nullptr); } std::optional EventGenerator::gen_event(){ HEJFOG::PhaseSpacePoint psp{ process_, jets_, pdf_, beam_.energy, subl_chance_, subl_channels_, particle_decays_, ew_parameters_, *ran_ }; status_ = psp.status(); if(status_ != Status::good) return {}; HEJ::Event ev = scale_gen_( HEJ::Event{ to_EventData( std::move(psp) ).cluster( jets_.def, jets_.min_pt) } ); if(!is_resummable(ev.type())){ throw HEJ::not_implemented("Tried to generate a event type, " "which is not yet implemented in HEJ."); } ev.generate_colours(*ran_); const double shat = HEJ::shat(ev); const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy); const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy); // evaluate matrix element ev.parameters() *= ME_.tree(ev)/(shat*shat); // and PDFs ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type); ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type); for(std::size_t i = 0; i < ev.variations().size(); ++i){ auto & var = ev.variations(i); var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type); var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type); } return ev; } } // namespace HEJFOG diff --git a/FixedOrderGen/src/Subleading.cc b/FixedOrderGen/src/Subleading.cc new file mode 100644 index 0000000..053bfed --- /dev/null +++ b/FixedOrderGen/src/Subleading.cc @@ -0,0 +1,19 @@ +#include "Subleading.hh" +#include + +namespace HEJFOG { + using namespace subleading; + using namespace std::string_view_literals; + + std::string_view name(Channels const c){ + switch(c) { + case unordered: + return "unordered"sv; + case central_qqbar: + return "central qqbar"sv; + case extremal_qqbar: + return "extremal qqbar"sv; + } + throw std::logic_error{"Channel not included"}; + } +} diff --git a/FixedOrderGen/src/Utility.cc b/FixedOrderGen/src/Utility.cc new file mode 100644 index 0000000..09a8975 --- /dev/null +++ b/FixedOrderGen/src/Utility.cc @@ -0,0 +1,71 @@ +#include "Utility.hh" + +#include +#include +#include +#include + +#include "Process.hh" + +namespace HEJFOG { + using HEJ::is_anyquark; + + namespace { + bool is_up_type(HEJ::pid::ParticleID const part) { + return HEJ::is_anyquark(part) && ((std::abs(part) % 2) == 0); + } + bool is_down_type(HEJ::pid::ParticleID const part){ + return HEJ::is_anyquark(part) && ((std::abs(part)%2) != 0); + } + } + + bool parton_can_couple_to_W( + HEJ::pid::ParticleID const part, HEJ::pid::ParticleID const W_id + ){ + const int W_charge = W_id>0?1:-1; + return std::abs(part) 0 && is_up_type(part)) + || (W_charge*part < 0 && is_down_type(part)) ); + } + + bool parton_can_couple_to_W( + HEJ::Particle const & part, HEJ::pid::ParticleID const W_id + ){ + return parton_can_couple_to_W(part.type, W_id); + } + + bool is_AWZ_process(Process const & proc){ + return proc.boson && HEJ::is_AWZ_boson(*proc.boson); + } + + bool vector_boson_can_couple_to( + HEJ::pid::ParticleID boson, + HEJ::pid::ParticleID particle + ) { + if(!HEJ::is_AWZ_boson(boson)) { + if(HEJ::is_AWZ_boson(particle)) { + std::swap(boson, particle); + } else { + throw std::invalid_argument{"One argument has to be a vector boson"}; + } + } + switch(particle){ + case HEJ::pid::proton: + case HEJ::pid::antiproton: + return true; + case HEJ::pid::gluon: + return false; + default: + if (!is_anyquark(particle)) { + throw HEJ::not_implemented{ + "Test if boson couples to " + name(particle) + }; + } + if(std::abs(boson) == HEJ::pid::ParticleID::Wp) { + return parton_can_couple_to_W(particle, boson); + } + assert(boson == HEJ::pid::Z_photon_mix); + return true; + } + } +} // namespace HEJFOG diff --git a/FixedOrderGen/t/PSP_channel.cc b/FixedOrderGen/t/PSP_channel.cc index 2e9f658..fdec315 100644 --- a/FixedOrderGen/t/PSP_channel.cc +++ b/FixedOrderGen/t/PSP_channel.cc @@ -1,250 +1,238 @@ /** * check that the sum of all possible quarks is the same as a proton * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2020 * \copyright GPLv2 or later */ #include #include #include #include #include "HEJ/Mixmax.hh" #include "HEJ/Event.hh" #include "config.hh" #include "EventGenerator.hh" namespace { constexpr double max_sigma = 3.6; } //! throw error if condition not fulfilled #define ASSERT(x) if(!(x)) { \ throw std::logic_error("Assertion '" #x "' failed."); \ } //! throw error if condition not fulfilled #define ASSERT_THROW(x, exception) try { \ x; \ std::cerr << "'" #x "' did not throw an exception.\n"; \ throw; \ } catch(exception const & ex){ \ std::cout << "throw triggered: " << ex.what() << std::endl; \ } \ catch (...) { \ std::cerr << "Unexpected exception thrown for '" #x "'.\n"; \ throw; \ } namespace { const static std::array PARTONS{ HEJ::ParticleID::g, HEJ::ParticleID::d, HEJ::ParticleID::u, HEJ::ParticleID::s, HEJ::ParticleID::c, HEJ::ParticleID::b, HEJ::ParticleID::d_bar, HEJ::ParticleID::u_bar, HEJ::ParticleID::s_bar, HEJ::ParticleID::c_bar, HEJ::ParticleID::b_bar }; bool only_channel( HEJFOG::Subleading channels, HEJFOG::subleading::Channels selected ){ return channels[selected] && channels.reset(selected).none(); } - void generate_till_throw(HEJFOG::EventGenerator & gen ){ - for(std::size_t i=0; i<100; ++i){ - auto ev = gen.gen_event(); - if(gen.status()==HEJFOG::Status::good){ - std::cerr << "Wrongfully generated valid event!" << *ev <<"\n"; - throw; - } - } - std::cerr << "Unable to generate a proper throw!\n"; - throw; - } - bool can_couple_to_W(HEJ::ParticleID const id, HEJ::ParticleID const boson){ using namespace HEJ; if(!is_anyquark(id)){ return false; } if(std::abs(id)==HEJ::ParticleID::b || std::abs(id)==HEJ::ParticleID::top) return false; // Wp: if(boson>0){ // anti-down if(id%2==0){ return is_quark(id); } // or up-type quark return is_antiquark(id); } // Wm: // down if(id%2==0){ return is_antiquark(id); } // or anti-up-type quark return is_quark(id); } bool can_couple_to_Z(HEJ::ParticleID const id){ using namespace HEJ; return is_anyquark(id); } } int main(int argc, char const *argv[]) { if(argc != 2 && argc != 3){ std::cerr << "Usage: " << argv[0] << " config.yaml\n"; return EXIT_FAILURE; } const bool short_only = argc==3; std::cout <outgoing()), end(ev->outgoing()), [&](HEJ::Particle const & p){ return p.type == *config.process.boson; } ); if(boson->m() < 81. || boson->m() > 101.){ continue; } } double const wgt = ev->central().weight/config.events; tot_weight += wgt; tot_err += wgt*wgt; } } ASSERT(tot_weight!=0.); } tot_err = std::sqrt(tot_err); config.events /= PARTONS.size(); double tot_channels = 0.; double err_channels = 0.; for(auto b1: PARTONS){ for(auto b2: PARTONS){ double wgt_channel = 0.; double err_channel = 0.; config.process.incoming[0] = b1; config.process.incoming[1] = b2; - HEJFOG::EventGenerator gen{config}; std::cout << name(b1) << "+" << name(b2) << " " <3 && config.subleading_fraction>0. && config.subleading_channels[HEJFOG::subleading::cqqbar]){ // this will force central qqbar } else if(config.process.njets>2 && config.subleading_fraction>0. && config.subleading_channels[HEJFOG::subleading::eqqbar]){ // this will force extremal qqbar } else { auto const boson = *config.process.boson; if(!can_couple_to_W(b1, boson) && !can_couple_to_W(b2, boson) ){ std::cout << "bad " << name(boson) << std::endl; - ASSERT_THROW(generate_till_throw(gen), std::invalid_argument); + ASSERT_THROW(HEJFOG::EventGenerator{config}, std::invalid_argument); continue; } } } // for Z+jets we need quarks if(config.process.boson && *config.process.boson == HEJ::ParticleID::Z_photon_mix){ if(config.process.njets>3 && config.subleading_fraction>0. && config.subleading_channels[HEJFOG::subleading::cqqbar]){ // this will force central qqbar } else if(config.process.njets>2 && config.subleading_fraction>0. && config.subleading_channels[HEJFOG::subleading::eqqbar]){ // this will force extremal qqbar } else { auto const boson = *config.process.boson; if(!can_couple_to_Z(b1) && !can_couple_to_Z(b2) ){ std::cout << "bad " << name(boson) << std::endl; - ASSERT_THROW(generate_till_throw(gen), std::invalid_argument); + ASSERT_THROW(HEJFOG::EventGenerator{config}, std::invalid_argument); continue; } } } + // everything else should work + HEJFOG::EventGenerator gen{config}; for(std::size_t i=0; ioutgoing()), end(ev->outgoing()), [&](HEJ::Particle const & p){ return p.type == *config.process.boson; } ); if(boson->m() < 81. || boson->m() > 101.){ continue; } } double const wgt = ev->central().weight/config.events; wgt_channel += wgt; err_channel += wgt*wgt; } } ASSERT(wgt_channel!=0.); tot_channels += wgt_channel; err_channels += err_channel; err_channel = std::sqrt(err_channel); std::cout << "=> " << wgt_channel << " +/- " << err_channel << std::endl; } } err_channels = std::sqrt(err_channels); std::cout << tot_weight << " +/- " << tot_err << " vs. " << tot_channels << " +/- " << err_channels << std::endl; const double deviation = std::abs(tot_weight - tot_channels) / sqrt(err_channels*err_channels+tot_err*tot_err); std::cout << "(" << deviation << " sigma)" << std::endl; ASSERT(deviation < max_sigma); return EXIT_SUCCESS; }