Page MenuHomeHEPForge
Files F9514797

No OneTemporary
Actions

View Options

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 <HEJ/Particle.hh>
#include <cstddef>
+#include <iostream>
+#include <stdexcept>
#include <utility>
#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<subleading::Channels>(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<HEJ::RNG> 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<HEJ::Event> 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/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 <array>
#include <cmath>
#include <iostream>
#include <memory>
#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<HEJ::ParticleID, 11> 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 <<std::scientific;
auto config = HEJFOG::load_config(argv[1]);
config.events/=2.;
if(short_only)
config.events/=20.;
double tot_weight = 0.;
double tot_err = 0.;
config.process.incoming[0] = config.process.incoming[1] = HEJ::ParticleID::proton;
{
HEJFOG::EventGenerator gen{config};
for(std::size_t i=0; i<config.events; ++i){
auto const ev = gen.gen_event();
if(ev){
if(config.process.boson && *config.process.boson == HEJ::ParticleID::Z_photon_mix){
const auto boson = std::find_if(
begin(ev->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) << " " <<std::flush;
// for pure subleading configurations some setups states should throw
if(config.subleading_fraction == 1.){
if( only_channel(config.subleading_channels, HEJFOG::subleading::uno) ){
if(HEJ::is_gluon(b1) && HEJ::is_gluon(b2)){
std::cout << "uno" << std::endl;
- ASSERT_THROW(generate_till_throw(gen), std::invalid_argument);
+ ASSERT_THROW(HEJFOG::EventGenerator{config}, std::invalid_argument);
continue;
}
}
if( only_channel(config.subleading_channels, HEJFOG::subleading::eqqbar) ){
if(HEJ::is_anyquark(b1) && HEJ::is_anyquark(b2)){
std::cout << "eqqbar" << std::endl;
- ASSERT_THROW(generate_till_throw(gen), std::invalid_argument);
+ ASSERT_THROW(HEJFOG::EventGenerator{config}, std::invalid_argument);
continue;
}
}
}
// for W+jets we need the correct quarks
if(config.process.boson
&& std::abs(*config.process.boson) == HEJ::ParticleID::Wp
){
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_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; i<config.events; ++i){
- // everything else should work
auto const ev = gen.gen_event();
if(ev){
if(config.process.boson && *config.process.boson == HEJ::ParticleID::Z_photon_mix){
const auto boson = std::find_if(
begin(ev->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;
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;
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Feb 24, 6:45 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4493371
Default Alt Text
(16 KB)

Event Timeline