Page MenuHomeHEPForge
Files F8723909

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/CMakeLists.txt b/FixedOrderGen/CMakeLists.txt
index 48d00ec..d824bb7 100644
--- a/FixedOrderGen/CMakeLists.txt
+++ b/FixedOrderGen/CMakeLists.txt
@@ -1,112 +1,113 @@
cmake_minimum_required(VERSION 3.1 FATAL_ERROR)
set(CMAKE_LEGACY_CYGWIN_WIN32 0)
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
project("HEJ Fixed Order Generation" VERSION 2.0.5 LANGUAGES C CXX)
# Set a default build type if none was specified
set(default_build_type "RelWithDebInfo")
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "Setting build type to '${default_build_type}' as none was specified.")
set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
STRING "Choose the type of build." FORCE)
# Set the possible values of build type for cmake-gui
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
## Flags for the compiler. No warning allowed.
if (CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra")
elseif (MSVC)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /W4 /WX /EHsc")
endif()
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_STANDARD 14)
## Create Version
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")
# Get the latest abbreviated commit hash of the working branch
execute_process(
COMMAND git rev-parse HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
OUTPUT_VARIABLE PROJECT_GIT_REVISION
OUTPUT_STRIP_TRAILING_WHITESPACE
)
# Get the current working branch
execute_process(
COMMAND git rev-parse --abbrev-ref HEAD
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
OUTPUT_VARIABLE PROJECT_GIT_BRANCH
OUTPUT_STRIP_TRAILING_WHITESPACE
)
CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/Version.hh.in
${PROJECT_BINARY_DIR}/include/Version.hh @ONLY )
## Add directories and find dependences
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include ${PROJECT_BINARY_DIR}/include)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/../cmake/Modules/")
find_package(HEJ 2 REQUIRED)
include_directories(${HEJ_INCLUDE_DIR})
find_package(fastjet REQUIRED)
include_directories(${fastjet_INCLUDE_DIRS})
find_package(CLHEP 2.3 REQUIRED)
include_directories(${CLHEP_INCLUDE_DIRS})
find_package(LHAPDF REQUIRED)
include_directories(${LHAPDF_INCLUDE_DIRS})
## Amend unintuitive behaviour of FindBoost.cmake
## Priority of BOOST_ROOT over BOOSTROOT matches FindBoost.cmake
## at least for cmake 3.12
if(DEFINED BOOST_ROOT)
message("BOOST_ROOT set - only looking for Boost in ${BOOST_ROOT}")
set(Boost_NO_BOOST_CMAKE ON)
elseif(DEFINED BOOSTROOT)
message("BOOSTROOT set - only looking for Boost in ${BOOSTROOT}")
set(Boost_NO_BOOST_CMAKE ON)
endif()
find_package(Boost REQUIRED COMPONENTS iostreams)
include_directories(${Boost_INCLUDE_DIRS})
find_package(yaml-cpp) # requiring yaml does not work with fedora
include_directories(${YAML_CPP_INCLUDE_DIR})
## define executable
file(GLOB HEJFOG_source ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cc)
list(REMOVE_ITEM HEJFOG_source ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cc)
add_library(hejfog STATIC ${HEJFOG_source})
add_executable(HEJFOG ${CMAKE_CURRENT_SOURCE_DIR}/src/main.cc)
## link libraries
set(libraries ${CMAKE_DL_LIBS} ${LHAPDF_LIBRARIES} ${CLHEP_LIBRARIES}
${fastjet_LIBRARIES} ${Boost_LIBRARIES} ${YAML_CPP_LIBRARIES} ${HEJ_LIBRARIES})
target_link_libraries(hejfog ${libraries})
target_link_libraries(HEJFOG hejfog)
install(TARGETS HEJFOG DESTINATION bin)
## tests
enable_testing()
set(tst_dir "${CMAKE_CURRENT_SOURCE_DIR}/t")
-foreach(tst h_2j h_3j h_5j h_3j_uno1 h_3j_uno2 h_2j_decay 2j 4j)# W_2j_classify W_nj_classify)
+foreach(tst h_2j h_3j h_5j h_3j_uno1 h_3j_uno2 h_2j_decay 2j 4j W_reconstruct_enu
+ W_2j_classify W_nj_classify)
add_executable(test_${tst} ${tst_dir}/${tst}.cc)
target_link_libraries(test_${tst} hejfog)
add_test(NAME t_${tst} COMMAND test_${tst} WORKING_DIRECTORY ${tst_dir})
endforeach()
add_test(
NAME t_main_2j
COMMAND HEJFOG ${tst_dir}/config_2j.yml
)
add_test(
NAME t_main_h2j
COMMAND HEJFOG ${tst_dir}/config_h_2j.yml
)
add_test(
NAME t_main_h2j_decay
COMMAND HEJFOG ${tst_dir}/config_h_2j_decay.yml
)
diff --git a/FixedOrderGen/include/Decay.hh b/FixedOrderGen/include/Decay.hh
index 64cae8d..05d9c75 100644
--- a/FixedOrderGen/include/Decay.hh
+++ b/FixedOrderGen/include/Decay.hh
@@ -1,16 +1,16 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include "HEJ/PDG_codes.hh"
#include <vector>
namespace HEJFOG{
struct Decay{
- std::vector<HEJ::pid::ParticleID> products;
+ std::vector<HEJ::ParticleID> products;
double branching_ratio;
};
}
diff --git a/FixedOrderGen/include/PhaseSpacePoint.hh b/FixedOrderGen/include/PhaseSpacePoint.hh
index 0227f7c..3fc11b9 100644
--- a/FixedOrderGen/include/PhaseSpacePoint.hh
+++ b/FixedOrderGen/include/PhaseSpacePoint.hh
@@ -1,221 +1,225 @@
/** \file PhaseSpacePoint.hh
* \brief Contains the PhaseSpacePoint Class
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <bitset>
#include <vector>
#include "HEJ/Event.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDF.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/RNG.hh"
#include "JetParameters.hh"
#include "ParticleProperties.hh"
#include "Status.hh"
namespace HEJFOG{
class Process;
using HEJ::Particle;
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! Default PhaseSpacePoint Constructor
PhaseSpacePoint() = default;
//! PhaseSpacePoint Constructor
/**
* @param proc The process to generate
* @param jet_properties Jet defintion & cuts
* @param pdf The pdf set (used for sampling)
* @param E_beam Energie of the beam
* @param subl_chance Chance to turn a potentially unordered
* emission into an actual one
* @param subl_channels Possible subleading channels.
* see HEJFOG::Subleading
* @param particle_properties Properties of producted boson
*
* Initially, only FKL phase space points are generated. subl_chance gives
* the change of turning one emissions into a subleading configuration,
* i.e. either unordered or central quark/anti-quark pair. Unordered
* emissions require that the most extremal emission in any direction is
* a quark or anti-quark and the next emission is a gluon. Quark/anti-quark
* pairs are only generated for W processes. At most one subleading
* emission will be generated in this way.
*/
PhaseSpacePoint(
Process const & proc,
JetParameters const & jet_properties,
HEJ::PDF & pdf, double E_beam,
double subl_chance,
unsigned int subl_channels,
ParticlesPropMap const & particles_properties,
HEJ::RNG & ran
);
//! Get Weight Function
/**
* @returns Weight of Event
*/
double weight() const{
return weight_;
}
Status status() const{
return status_;
}
//! Get Incoming Function
/**
* @returns Incoming Particles
*/
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Get Outgoing Function
/**
* @returns Outgoing Particles
*/
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
return decays_;
}
private:
/**
* @internal
* @brief Generate LO parton momentum
*
* @param count Number of partons to generate
* @param is_pure_jets If true ensures momentum conservation in x and y
* @param jet_param Jet properties to fulfil
* @param max_pt max allowed pt for a parton (typically E_CMS)
* @param ran Random Number Generator
*
* @returns Momentum of partons
*
* Ensures that each parton is in its own jet.
* Generation is independent of parton flavour. Output is sorted in rapidity.
*/
std::vector<fastjet::PseudoJet> gen_LO_partons(
int count, bool is_pure_jets,
JetParameters const & jet_param,
double max_pt,
HEJ::RNG & ran
);
- std::vector<Particle> gen_enu(
- std::vector<HEJ::pid::ParticleID> const & pair, HEJ::RNG & ran
- );
Particle gen_boson(
HEJ::ParticleID bosonid, double mass, double width,
HEJ::RNG & ran
);
template<class ParticleMomenta>
fastjet::PseudoJet gen_last_momentum(
ParticleMomenta const & other_momenta,
double mass_square, double y
) const;
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
/**
* @internal
* @brief Generate incoming partons according to the PDF
*
* @param uf Scale used in the PDF
*/
void reconstruct_incoming(
Process const & proc, unsigned int subl_channels,
HEJ::PDF & pdf, double E_beam,
double uf,
HEJ::RNG & ran
);
/**
* @internal
* @brief Returns list of all allowed initial states partons
*/
std::array<std::bitset<11>,2> filter_partons(
Process const & proc, unsigned int const subl_channels,
HEJ::RNG & ran
);
HEJ::ParticleID generate_incoming_id(
size_t beam_idx, double x, double uf, HEJ::PDF & pdf,
std::bitset<11> allowed_partons, HEJ::RNG & ran
);
bool momentum_conserved(double ep) const;
HEJ::Particle const & most_backward_FKL(
std::vector<HEJ::Particle> const & partons
) const;
HEJ::Particle const & most_forward_FKL(
std::vector<HEJ::Particle> const & partons
) const;
HEJ::Particle & most_backward_FKL(std::vector<HEJ::Particle> & partons) const;
HEJ::Particle & most_forward_FKL(std::vector<HEJ::Particle> & partons) const;
bool extremal_FKL_ok(
std::vector<fastjet::PseudoJet> const & partons
) const;
double random_normal(double stddev, HEJ::RNG & ran);
/**
* @internal
* @brief Turns a FKL configuration into a subleading one
*
* @param chance Change to switch to subleading configuration
* @param channels Allowed channels for subleading process
* @param proc Process to decide which subleading
* configurations are allowed
*
* With a chance of "chance" the FKL configuration is either turned into
* a unordered configuration or, for A/W/Z bosons, a configuration with
* a central quark/anti-quark pair.
*/
void maybe_turn_to_subl(double chance, unsigned int channels,
Process const & proc, HEJ::RNG & ran);
void turn_to_uno(bool can_be_uno_backward, bool can_be_uno_forward, HEJ::RNG & ran);
void turn_to_qqx(bool allow_strange, HEJ::RNG & ran);
+ //! decay where we select the decay channel
std::vector<Particle> decay_boson(
HEJ::Particle const & parent,
std::vector<Decay> const & decays,
HEJ::RNG & ran
);
+ //! generate decay products of a boson
+ std::vector<Particle> decay_boson(
+ HEJ::Particle const & parent,
+ std::vector<HEJ::ParticleID> const & decays,
+ HEJ::RNG & ran
+ );
/// @brief setup outgoing partons to ensure correct coupling to boson
void couple_boson(HEJ::ParticleID boson, HEJ::RNG & ran);
Decay select_decay_channel(
std::vector<Decay> const & decays,
HEJ::RNG & ran
);
double gen_hard_pt(
int np, double ptmin, double ptmax, double y,
HEJ::RNG & ran
);
double gen_soft_pt(int np, double ptmax, HEJ::RNG & ran);
double gen_parton_pt(
int count, JetParameters const & jet_param, double ptmax, double y,
HEJ::RNG & ran
);
double weight_;
Status status_;
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<size_t, std::vector<Particle>> decays_;
};
HEJ::Event::EventData to_EventData(PhaseSpacePoint const & psp);
}
diff --git a/FixedOrderGen/include/Process.hh b/FixedOrderGen/include/Process.hh
index d2b0beb..aef9426 100644
--- a/FixedOrderGen/include/Process.hh
+++ b/FixedOrderGen/include/Process.hh
@@ -1,22 +1,22 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <vector>
#include "HEJ/PDG_codes.hh"
#include "HEJ/optional.hh"
namespace HEJFOG{
struct Process{
std::array<HEJ::ParticleID, 2> incoming;
int njets;
HEJ::optional<HEJ::ParticleID> boson;
- std::vector <HEJ::ParticleID> leptons;
+ std::vector<HEJ::ParticleID> boson_decay;
};
}
diff --git a/FixedOrderGen/src/PhaseSpacePoint.cc b/FixedOrderGen/src/PhaseSpacePoint.cc
index d5d80d5..02880b6 100644
--- a/FixedOrderGen/src/PhaseSpacePoint.cc
+++ b/FixedOrderGen/src/PhaseSpacePoint.cc
@@ -1,742 +1,695 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "PhaseSpacePoint.hh"
#include <algorithm>
#include "CLHEP/Vector/LorentzVector.h"
#include "HEJ/Constants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/kinematics.hh"
#include "HEJ/Particle.hh"
#include "HEJ/utility.hh"
#include "Process.hh"
#include "Subleading.hh"
using namespace HEJ;
namespace HEJFOG{
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
HEJ::Event::EventData to_EventData(PhaseSpacePoint const & psp){
HEJ::Event::EventData result;
result.incoming = psp.incoming();
assert(result.incoming.size() == 2);
result.outgoing=psp.outgoing();
// technically Event::EventData doesn't have to be sorted,
// but PhaseSpacePoint should be anyway
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
HEJ::rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.decays=psp.decays();
result.parameters.central= {NaN, NaN, psp.weight() };
return result;
}
namespace{
bool can_swap_to_uno(
HEJ::Particle const & p1, HEJ::Particle const & p2
){
return is_parton(p1)
&& p1.type != pid::gluon
&& p2.type == pid::gluon;
}
size_t count_gluons(std::vector<Particle>::const_iterator first,
std::vector<Particle>::const_iterator last){
return std::count_if(first, last, [](Particle const & p)
{return p.type == pid::gluon;});
}
/** assumes FKL configurations between first and last,
* else there can be a quark in a non-extreme position
* e.g. uno configuration gqg would pass
*/
bool can_change_to_qqx(
std::vector<Particle>::const_iterator first,
std::vector<Particle>::const_iterator last){
return 1 < count_gluons(first,last);
}
bool is_AWZ_proccess(Process const & proc){
return proc.boson && is_AWZ_boson(*proc.boson);
}
bool is_up_type(Particle const & part){
return HEJ::is_anyquark(part) && !(abs(part.type)%2);
}
bool is_down_type(Particle const & part){
return HEJ::is_anyquark(part) && abs(part.type)%2;
}
bool can_couple_to_W(Particle const & part, pid::ParticleID const W_id){
const int W_charge = W_id>0?1:-1;
return abs(part.type)<5
&& ( (W_charge*part.type > 0 && is_up_type(part))
|| (W_charge*part.type < 0 && is_down_type(part)) );
}
}
void PhaseSpacePoint::maybe_turn_to_subl(
double chance,
unsigned int const channels,
Process const & proc,
HEJ::RNG & ran
){
if(proc.njets <= 2) return;
assert(outgoing_.size() >= 2);
// decide what kind of subleading process is allowed
bool allow_uno = false;
bool allow_strange = true;
const size_t nout = outgoing_.size();
const bool can_be_uno_backward = (channels&Subleading::uno)
&& can_swap_to_uno(outgoing_[0], outgoing_[1]);
const bool can_be_uno_forward = (channels&Subleading::uno)
&& can_swap_to_uno(outgoing_[nout-1], outgoing_[nout-2]);
allow_uno = can_be_uno_backward || can_be_uno_forward;
bool allow_qqx = false;
if(is_AWZ_proccess(proc)) {
allow_qqx = (channels&Subleading::qqx)
&& can_change_to_qqx(outgoing_.cbegin(), outgoing_.cend());
if(std::none_of(outgoing_.cbegin(), outgoing_.cend(),
[&proc](Particle const & p){ return can_couple_to_W(p, *proc.boson);})) {
// enforce qqx if A/W/Z can't couple somewhere else
assert(allow_qqx);
allow_uno = false;
chance = 1.;
// strange not allowed for W
if(abs(*proc.boson)== pid::Wp) allow_strange = false;
}
}
if(!allow_uno && !allow_qqx) return;
if(ran.flat() < chance){
weight_ /= chance;
if(allow_uno && !allow_qqx){
turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
} else if (!allow_uno && allow_qqx) {
turn_to_qqx(allow_strange, ran);
} else {
assert( allow_uno && allow_qqx);
if(ran.flat() < 0.5) turn_to_uno(can_be_uno_backward, can_be_uno_forward, ran);
else turn_to_qqx(allow_strange, ran);
weight_ *= 2.;
}
} else weight_ /= 1 - chance;
}
void PhaseSpacePoint::turn_to_uno(
const bool can_be_uno_backward, const bool can_be_uno_forward,
HEJ::RNG & ran
){
if(!can_be_uno_backward && !can_be_uno_forward) return;
const size_t nout = outgoing_.size();
if(can_be_uno_backward && can_be_uno_forward){
if(ran.flat() < 0.5){
std::swap(outgoing_[0].type, outgoing_[1].type);
} else {
std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
}
weight_ *= 2.;
} else if(can_be_uno_backward){
std::swap(outgoing_[0].type, outgoing_[1].type);
} else {
assert(can_be_uno_forward);
std::swap(outgoing_[nout-1].type, outgoing_[nout-2].type);
}
}
void PhaseSpacePoint::turn_to_qqx(const bool allow_strange, HEJ::RNG & ran){
/// find first and last gluon in FKL chain
auto first = std::find_if(outgoing_.begin(), outgoing_.end(),
[](Particle const & p){return p.type == pid::gluon;});
std::vector<Particle*> FKL_gluons;
for(auto p = first; p!=outgoing_.end(); ++p){
if(p->type == pid::gluon) FKL_gluons.push_back(&*p);
else if(is_anyquark(*p)) break;
}
const size_t ng = FKL_gluons.size();
if(ng < 2)
throw std::logic_error("not enough gluons to create qqx");
// select flavour of quark
const double r1 = 2.*ran.flat()-1.;
const double max_flavour = allow_strange?n_f:n_f-1;
weight_ *= max_flavour*2;
int flavour = pid::down + std::floor(std::abs(r1)*max_flavour);
flavour*=r1<0.?-1:1;
// select gluon for switch
const size_t idx = floor((ng-1) * ran.flat());
weight_ *= (ng-1);
FKL_gluons[idx]->type = ParticleID(flavour);
FKL_gluons[idx+1]->type = ParticleID(-flavour);
}
template<class ParticleMomenta>
fastjet::PseudoJet PhaseSpacePoint::gen_last_momentum(
ParticleMomenta const & other_momenta,
const double mass_square, const double y
) const {
std::array<double,2> pt{0.,0.};
for (auto const & p: other_momenta) {
pt[0]-= p.px();
pt[1]-= p.py();
}
const double mperp = sqrt(pt[0]*pt[0]+pt[1]*pt[1]+mass_square);
const double pz=mperp*sinh(y);
const double E=mperp*cosh(y);
return {pt[0], pt[1], pz, E};
}
namespace {
//! adds a particle to target (in correct rapidity ordering)
//! @returns positon of insertion
auto insert_particle(std::vector<HEJ::Particle> & target,
HEJ::Particle && particle
){
const auto pos = std::upper_bound(
begin(target),end(target),particle,rapidity_less{}
);
target.insert(pos, std::move(particle));
return pos;
}
}
PhaseSpacePoint::PhaseSpacePoint(
- Process const & proc,
+ Process const & proc_in,
JetParameters const & jet_param,
HEJ::PDF & pdf, double E_beam,
double const subl_chance,
unsigned int const subl_channels,
- ParticlesPropMap const & particles_properties,
+ ParticlesPropMap const & particles_properties_in,
HEJ::RNG & ran
)
{
+ auto proc{ proc_in };
+ auto particles_properties{ particles_properties_in };
assert(proc.njets >= 2);
if(proc.boson
&& particles_properties.find(*(proc.boson))
== particles_properties.end())
throw HEJ::missing_option("Boson "
+std::to_string(*(proc.boson))+" can't be generated: missing properties");
status_ = good;
weight_ = 1;
- const int nout = proc.njets + (proc.boson?1:0) + proc.leptons.size();
+ const int nout = proc.njets + (proc.boson?1:0) + proc.boson_decay.size();
outgoing_.reserve(nout);
// generate parton momenta
const bool is_pure_jets = (nout == proc.njets);
auto partons = gen_LO_partons(
proc.njets, is_pure_jets, jet_param, E_beam, ran
);
// pre fill flavour with gluons
for(auto&& p_out: partons) {
outgoing_.emplace_back(Particle{pid::gluon, std::move(p_out), {}});
}
if(status_ != good) return;
- // create phase space point for (lepton,neutrino)-pair of W+/W-
- // @TODO this have massive overlap with the "create boson" section
- if(proc.leptons.size()==2){
- auto pair(gen_enu(proc.leptons,ran));
- for(auto&& p: pair)
- insert_particle(outgoing_, std::move(p));
- }
-
- // create boson
- if(proc.boson){
+ if(proc.boson){ // decay boson
const auto & boson_prop = particles_properties.at(*proc.boson);
- auto boson(gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran));
+ auto boson{ gen_boson(*proc.boson, boson_prop.mass, boson_prop.width, ran) };
const auto pos{insert_particle(outgoing_, std::move(boson))};
- if(! boson_prop.decays.empty()){
- const size_t boson_idx = std::distance(begin(outgoing_), pos);
+ const size_t boson_idx = std::distance(begin(outgoing_), pos);
+ if( !proc.boson_decay.empty() ){
+ decays_.emplace(
+ boson_idx,
+ decay_boson(outgoing_[boson_idx], proc.boson_decay, ran)
+ );
+ } else if( !boson_prop.decays.empty() ){
decays_.emplace(
boson_idx,
decay_boson(outgoing_[boson_idx], boson_prop.decays, ran)
);
}
}
// normalisation of momentum-conserving delta function
weight_ *= pow(2*M_PI, 4);
/** @TODO
* uf (jet_param.min_pt) doesn't correspond to our final scale choice.
* The HEJ scale generators currently expect a full event as input,
* so fixing this is not completely trivial
*/
reconstruct_incoming(proc, subl_channels, pdf, E_beam, jet_param.min_pt, ran);
if(status_ != good) return;
// set outgoing states
most_backward_FKL(outgoing_).type = incoming_[0].type;
most_forward_FKL(outgoing_).type = incoming_[1].type;
maybe_turn_to_subl(subl_chance, subl_channels, proc, ran);
if(proc.boson) couple_boson(*proc.boson, ran);
}
double PhaseSpacePoint::gen_hard_pt(
int np , double ptmin, double ptmax, double y,
HEJ::RNG & ran
) {
// heuristic parameters for pt sampling
const double ptpar = ptmin + np/5.;
const double arg_small_y = atan((ptmax - ptmin)/ptpar);
const double y_cut = 3.;
const double r1 = ran.flat();
if(y < y_cut){
const double pt = ptmin + ptpar*tan(r1*arg_small_y);
const double temp = cos(r1*arg_small_y);
weight_ *= pt*ptpar*arg_small_y/(temp*temp);
return pt;
}
const double ptpar2 = ptpar/(1 + 5*(y-y_cut));
const double temp = 1. - std::exp((ptmin-ptmax)/ptpar2);
const double pt = ptmin - ptpar2*std::log(1-r1*temp);
weight_ *= pt*ptpar2*temp/(1-r1*temp);
return pt;
}
double PhaseSpacePoint::gen_soft_pt(int np, double max_pt, HEJ::RNG & ran) {
constexpr double ptpar = 4.;
const double r = ran.flat();
const double pt = max_pt + ptpar/np*std::log(r);
weight_ *= pt*ptpar/(np*r);
return pt;
}
double PhaseSpacePoint::gen_parton_pt(
int count, JetParameters const & jet_param, double max_pt, double y,
HEJ::RNG & ran
) {
constexpr double p_small_pt = 0.02;
if(! jet_param.peak_pt) {
return gen_hard_pt(count, jet_param.min_pt, max_pt, y, ran);
}
const double r = ran.flat();
if(r > p_small_pt) {
weight_ /= 1. - p_small_pt;
return gen_hard_pt(count, *jet_param.peak_pt, max_pt, y, ran);
}
weight_ /= p_small_pt;
const double pt = gen_soft_pt(count, *jet_param.peak_pt, ran);
if(pt < jet_param.min_pt) {
weight_=0.0;
status_ = not_enough_jets;
return jet_param.min_pt;
}
return pt;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_LO_partons(
int np, bool is_pure_jets,
JetParameters const & jet_param,
double max_pt,
HEJ::RNG & ran
){
if (np<2) throw std::invalid_argument{"Not enough partons in gen_LO_partons"};
weight_ /= pow(16.*pow(M_PI,3),np);
weight_ /= std::tgamma(np+1); //remove rapidity ordering
std::vector<fastjet::PseudoJet> partons;
partons.reserve(np);
for(int i = 0; i < np; ++i){
const double y = -jet_param.max_y + 2*jet_param.max_y*ran.flat();
weight_ *= 2*jet_param.max_y;
const bool is_last_parton = i+1 == np;
if(is_pure_jets && is_last_parton) {
constexpr double parton_mass_sq = 0.;
partons.emplace_back(gen_last_momentum(partons, parton_mass_sq, y));
break;
}
const double phi = 2*M_PI*ran.flat();
weight_ *= 2.0*M_PI;
const double pt = gen_parton_pt(np, jet_param, max_pt, y, ran);
if(weight_ == 0.0) return {};
partons.emplace_back(fastjet::PtYPhiM(pt, y, phi));
assert(jet_param.min_pt <= partons[i].pt());
assert(partons[i].pt() <= max_pt+1e-5);
}
// Need to check that at LO, the number of jets = number of partons;
fastjet::ClusterSequence cs(partons, jet_param.def);
auto cluster_jets=cs.inclusive_jets(jet_param.min_pt);
if (cluster_jets.size()!=unsigned(np)){
weight_=0.0;
status_ = not_enough_jets;
return {};
}
std::sort(begin(partons), end(partons), rapidity_less{});
return partons;
}
- std::vector<Particle> PhaseSpacePoint::gen_enu(
- std::vector<HEJ::pid::ParticleID> const & pair, HEJ::RNG & ran
- ){
- assert(pair.size() == 2);
- // Now we know the transverse momentum of the W, given by the array kt
- // double pW[4];
- // Choose its mass according to Breit-Wigner
- double r1=ran.flat();
- double sW=HEJ::MW*(HEJ::MW + HEJ::GammaW*tan((M_PI*r1)/2. + (-1. + r1)*atan(HEJ::MW/HEJ::GammaW)));
- // Multiply by derivate (d sap)/(d r)
- {
- static double temp=atan(HEJ::MW/HEJ::GammaW);
- weight_*=(HEJ::GammaW*HEJ::MW*(M_PI+2.*temp))/(1.+cos(M_PI*r1+2.*(-1.+r1)*temp));
- }
-
- // Generate a yW Gaussian distributed around 0
- double yW;
- {
- double lninvr1,r1,r2,temp,a;
- r1=ran.flat();
- r2=ran.flat();
- lninvr1=-log(r1);
- a=0.7; // tuned number
- temp=a*sqrt(2.*lninvr1)*cos(2.*M_PI*r2);
- yW=temp;
- weight_=weight_*(exp(temp*temp/2./a/a))*sqrt(2.*M_PI)*a;
- }
- auto p = gen_last_momentum(outgoing_, sW, yW);
- CLHEP::HepLorentzVector pWv(p.px(),p.py(),p.pz(),p.e());
-
- double ppW[4],ppWs;
- CLHEP::HepLorentzVector ppWv,appWv,pd1,pd2;
- ppWs=sqrt(sW)/2.;
- // Choose theta and phi
- double pptheta=2.*M_PI*ran.flat();
- double cosppphi=2.*ran.flat()-1.;
- weight_*=2.*M_PI*2.;
- double sinphi=sqrt(1.-cosppphi*cosppphi); // Know 0 < phi < pi
- weight_*=1./(pow(2.*M_PI,3)*2.)/(4.); // Divide by 8, *2, see WPhaseSpace.tex @TODO where?!?
- // (phi -> cosphi now so no sinphi in J)
- // construct 4-vector in W rest frame
- ppW[0]=ppWs;
- ppW[1]=ppWs*cos(pptheta)*sinphi;
- ppW[2]=ppWs*sin(pptheta)*sinphi;
- ppW[3]=ppWs*cosppphi;
- ppWv.set(ppW[1],ppW[2],ppW[3],ppW[0]);
- appWv.set(-ppW[1],-ppW[2],-ppW[3],ppW[0]);
-
- // translate to lab frame
- pd1=ppWv.boost(pWv.boostVector()); //particle
- pd2=appWv.boost(pWv.boostVector()); //anti-particle
-
- // assign particle and anti-particle momentum
- // find the particle, assuming the pair is ordered (anti-lepton,lepton)
- HEJ::Particle alepton{pair[0],
- fastjet::PseudoJet{pd2.px(),pd2.py(),pd2.pz(),pd2.e()}, {}};
- HEJ::Particle lepton{pair[1],
- fastjet::PseudoJet{pd1.px(),pd1.py(),pd1.pz(),pd1.e()}, {}};
-
- return std::vector<HEJ::Particle>{alepton,lepton};
- }
-
Particle PhaseSpacePoint::gen_boson(
HEJ::ParticleID bosonid, double mass, double width,
HEJ::RNG & ran
){
// Usual phase space measure
weight_ /= 16.*pow(M_PI, 3);
// Generate a y Gaussian distributed around 0
- /// @TODO: magic number only for Higgs
+ /// @TODO check magic numbers for different boson Higgs
/// @TODO better sampling for W
- const double y = random_normal(1.6, ran);
+ const double stddev_y = 1.6;
+ const double y = random_normal(stddev_y, ran);
const double r1 = ran.flat();
- const double sH = mass*(
+ const double s_boson = mass*(
mass + width*tan(M_PI/2.*r1 + (r1-1.)*atan(mass/width))
);
+ // off-shell s_boson sampling, compensates for Breit-Wigner
+ /// @TODO use a flag instead
+ if(abs(bosonid) == pid::Wp){
+ weight_/=M_PI*M_PI*8.;
+ weight_*= mass*width*( M_PI+2.*atan(mass/width) )
+ / ( 1. + cos( M_PI*r1 + 2.*(r1-1.)*atan(mass/width) ) );
+ }
- auto p = gen_last_momentum(outgoing_, sH, y);
+ auto p = gen_last_momentum(outgoing_, s_boson, y);
return Particle{bosonid, std::move(p), {}};
}
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
if(!HEJ::is_parton(partons[0])) return partons[1];
return partons[0];
}
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t last_idx = partons.size() - 1;
if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
if(!HEJ::is_parton(partons[0])) return partons[1];
return partons[0];
}
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t last_idx = partons.size() - 1;
if(!HEJ::is_parton(partons[last_idx])) return partons[last_idx-1];
return partons[last_idx];
}
namespace {
/// partons are ordered: even = anti, 0 = gluon
ParticleID index_to_pid(size_t i){
if(!i) return pid::gluon;
return static_cast<ParticleID>(i%2?(i+1)/2:-i/2);
}
/// partons are ordered: even = anti, 0 = gluon
size_t pid_to_index(ParticleID id){
if(id==pid::gluon) return 0;
return id>0?id*2-1:abs(id)*2;
}
std::bitset<11> init_allowed(ParticleID const id){
if(abs(id) == pid::proton)
return ~0;
std::bitset<11> out = 0;
if(is_parton(id))
out[pid_to_index(id)] = 1;
return out;
}
/// decides which "index" (see index_to_pid) are allowed for process
std::bitset<11> allowed_quarks(ParticleID const boson){
std::bitset<11> allowed = ~0;
if(abs(boson) == pid::Wp){
// special case W:
// Wp: anti-down or up-type quark, no b/t -> 0001100110(1) = 205
// Wm: down or anti-up-type quark, no b/t -> 0010011001(1) = 307
allowed = boson>0?205:307;
}
return allowed;
}
}
/**
* checks which partons are allowed as initial state:
* 1. only allow what is given in the Runcard (p -> all)
* 2. A/W/Z require something to couple to
* a) no qqx => no incoming gluon
* b) 2j => no incoming gluon
* c) 3j => can couple OR is gluon => 2 gluons become qqx later
*/
std::array<std::bitset<11>,2> PhaseSpacePoint::filter_partons(
Process const & proc, unsigned int const subl_channels, HEJ::RNG & ran
){
std::array<std::bitset<11>,2> allowed_partons{
init_allowed(proc.incoming[0]),
init_allowed(proc.incoming[1])
};
bool const allow_qqx = subl_channels&Subleading::qqx;
// special case A/W/Z
if(is_AWZ_proccess(proc) && ((proc.njets < 4) || !allow_qqx)){
// all possible incoming states
auto allowed(allowed_quarks(*proc.boson));
if(proc.njets == 2 || !allow_qqx) allowed[0]=0;
// possible states per leg
std::array<std::bitset<11>,2> const maybe_partons{
allowed_partons[0]&allowed, allowed_partons[1]&allowed};
if(maybe_partons[0].any() && maybe_partons[1].any()){
// two options to get allowed initial state => choose one at random
const size_t idx = ran.flat() < 0.5;
allowed_partons[idx] = maybe_partons[idx];
// else choose the possible
} else if(maybe_partons[0].any()) {
allowed_partons[0] = maybe_partons[0];
} else if(maybe_partons[1].any()) {
allowed_partons[1] = maybe_partons[1];
} else{
throw std::invalid_argument{"Incoming state not allowed."};
}
}
return allowed_partons;
}
void PhaseSpacePoint::reconstruct_incoming(
Process const & proc, unsigned int const subl_channels,
HEJ::PDF & pdf, double E_beam,
double uf,
HEJ::RNG & ran
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
// calculate xa, xb
const double sqrts=2*E_beam;
const double xa=(incoming_[0].p.e()-incoming_[0].p.pz())/sqrts;
const double xb=(incoming_[1].p.e()+incoming_[1].p.pz())/sqrts;
// abort if phase space point is outside of collider energy reach
if (xa>1. || xb>1.){
weight_=0;
status_ = too_much_energy;
return;
}
auto const & ids = proc.incoming;
std::array<std::bitset<11>,2> allowed_partons(
filter_partons(proc, subl_channels, ran));
for(size_t i = 0; i < 2; ++i){
if(ids[i] == pid::proton || ids[i] == pid::p_bar){
// pick ids according to pdfs
incoming_[i].type =
generate_incoming_id(i, i?xb:xa, uf, pdf, allowed_partons[i], ran);
} else {
assert(allowed_partons[i][pid_to_index(ids[i])]);
incoming_[i].type = ids[i];
}
}
assert(momentum_conserved(1e-7));
}
HEJ::ParticleID PhaseSpacePoint::generate_incoming_id(
size_t const beam_idx, double const x, double const uf,
HEJ::PDF & pdf, std::bitset<11> allowed_partons, HEJ::RNG & ran
){
std::array<double,11> pdf_wt;
pdf_wt[0] = allowed_partons[0]?fabs(pdf.pdfpt(beam_idx,x,uf,pid::gluon)):0.;
double pdftot = pdf_wt[0];
for(size_t i = 1; i < pdf_wt.size(); ++i){
pdf_wt[i] = allowed_partons[i]?4./9.*fabs(pdf.pdfpt(beam_idx,x,uf,index_to_pid(i))):0;
pdftot += pdf_wt[i];
}
const double r1 = pdftot * ran.flat();
double sum = 0;
for(size_t i=0; i < pdf_wt.size(); ++i){
if (r1 < (sum+=pdf_wt[i])){
weight_*= pdftot/pdf_wt[i];
return index_to_pid(i);
}
}
std::cerr << "Error in choosing incoming parton: "<<x<<" "<<uf<<" "
<<sum<<" "<<pdftot<<" "<<r1<<std::endl;
throw std::logic_error{"Failed to choose parton flavour"};
}
void PhaseSpacePoint::couple_boson(
HEJ::ParticleID const boson, HEJ::RNG & ran
){
if(abs(boson) != pid::Wp) return; // only matters for W
/// @TODO this could be use to sanity check gamma and Z
// find all possible quarks
std::vector<Particle*> allowed_parts;
for(auto & part: outgoing_){
// Wp -> up OR anti-down, Wm -> anti-up OR down, no bottom
if ( can_couple_to_W(part, boson) )
allowed_parts.push_back(&part);
}
if(allowed_parts.size() == 0){
throw std::logic_error{"Found no parton for coupling with boson"};
}
// select one and flip it
size_t idx = 0;
if(allowed_parts.size() > 1){
/// @TODO more efficient sampling
/// old code: probability[i] = exp(parton[i].y - W.y)
idx = floor(ran.flat()*allowed_parts.size());
weight_ *= allowed_parts.size();
}
const int W_charge = boson>0?1:-1;
allowed_parts[idx]->type =
static_cast<ParticleID>( allowed_parts[idx]->type - W_charge );
}
double PhaseSpacePoint::random_normal(
double stddev,
HEJ::RNG & ran
){
const double r1 = ran.flat();
const double r2 = ran.flat();
const double lninvr1 = -log(r1);
const double result = stddev*sqrt(2.*lninvr1)*cos(2.*M_PI*r2);
weight_ *= exp(result*result/(2*stddev*stddev))*sqrt(2.*M_PI)*stddev;
return result;
}
bool PhaseSpacePoint::momentum_conserved(double ep) const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
for(auto const & out: outgoing()) diff -= out.p;
return nearby_ep(diff, fastjet::PseudoJet{}, ep);
}
Decay PhaseSpacePoint::select_decay_channel(
std::vector<Decay> const & decays,
HEJ::RNG & ran
){
double br_total = 0.;
for(auto const & decay: decays) br_total += decay.branching_ratio;
// adjust weight
// this is given by (channel branching ratio)/(chance to pick channel)
// where (chance to pick channel) =
// (channel branching ratio)/(total branching ratio)
weight_ *= br_total;
+ if(decays.size()==1) return decays.front();
const double r1 = br_total*ran.flat();
double br_sum = 0.;
for(auto const & decay: decays){
br_sum += decay.branching_ratio;
if(r1 < br_sum) return decay;
}
throw std::logic_error{"unreachable"};
}
std::vector<Particle> PhaseSpacePoint::decay_boson(
HEJ::Particle const & parent,
std::vector<Decay> const & decays,
HEJ::RNG & ran
){
const auto channel = select_decay_channel(decays, ran);
- if(channel.products.size() != 2){
+ return decay_boson(parent, channel.products, ran);
+ }
+
+ std::vector<Particle> PhaseSpacePoint::decay_boson(
+ HEJ::Particle const & parent,
+ std::vector<HEJ::ParticleID> const & decays,
+ HEJ::RNG & ran
+ ){
+ if(decays.size() != 2){
throw HEJ::not_implemented{
"only decays into two particles are implemented"
};
}
- std::vector<Particle> decay_products(channel.products.size());
- for(size_t i = 0; i < channel.products.size(); ++i){
- decay_products[i].type = channel.products[i];
+ std::vector<Particle> decay_products(decays.size());
+ for(size_t i = 0; i < decays.size(); ++i){
+ decay_products[i].type = decays[i];
}
// choose polar and azimuth angle in parent rest frame
const double E = parent.m()/2;
const double theta = 2.*M_PI*ran.flat();
const double cos_phi = 2.*ran.flat()-1.;
const double sin_phi = sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi
const double px = E*cos(theta)*sin_phi;
const double py = E*sin(theta)*sin_phi;
const double pz = E*cos_phi;
decay_products[0].p.reset(px, py, pz, E);
decay_products[1].p.reset(-px, -py, -pz, E);
for(auto & particle: decay_products) particle.p.boost(parent.p);
return decay_products;
}
}
diff --git a/FixedOrderGen/src/config.cc b/FixedOrderGen/src/config.cc
index cf43551..2a87eb6 100644
--- a/FixedOrderGen/src/config.cc
+++ b/FixedOrderGen/src/config.cc
@@ -1,393 +1,411 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "config.hh"
#include <cctype>
#include "Subleading.hh"
#include "HEJ/config.hh"
#include "HEJ/YAMLreader.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",
"event output", "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", "Wp", "W+", "Wm", "W-", "Z"}){
for(auto && particle_opt: {"mass", "width"}){
supported["particle properties"][particle_type][particle_opt] = "";
}
supported["particle properties"][particle_type]["decays"]["into"] = "";
supported["particle properties"][particle_type]["decays"]["branching ratio"] = "";
}
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<HEJ::ParticleID> 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<std::string> split(
std::string const & str, std::string const & delims
){
std::vector<std::string> result;
for(size_t begin, end = 0; end != str.npos;){
begin = str.find_first_not_of(delims, end);
if(begin == str.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', 'e-', 'e+', 'nu_e', 'nu_e_bar'"
+ " outgoing particles have to be 'j', 'photon', 'H', 'Wm', 'Wp', 'e-', 'e+', 'nu_e', 'nu_e_bar'"
};
}
- bool leptons_flavours_OK(std::vector<ParticleID> const & pid) {
- assert(pid.size()==2);
- assert(is_sorted(begin(pid),end(pid)));
- // Check it is a W+/W-
- const int pidsum = pid[0]+pid[1];
- if (HEJ::is_antineutrino(pid[0])){
- if (pidsum!=-1) return false;
- } else {
- if (pidsum!=1) return false;
+ HEJ::ParticleID reconstruct_boson_id(
+ std::vector<HEJ::ParticleID> 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;
}
- return true;
+ 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 "+HEJ::name(ids[0])+" and "+HEJ::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){
- if(result.boson){
+ 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.leptons.size()==2||result.boson){
- throw std::invalid_argument {
- "Too many leptons required "
- };
- }
- result.leptons.emplace_back(pid);
+ 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]);
}
}
}
- std::sort(std::begin(result.leptons),std::end(result.leptons));
if(result.njets < 2){
throw std::invalid_argument{
"Process has to include at least two jets ('j')"
};
}
- if((result.leptons.size()>0) && (!leptons_flavours_OK(result.leptons))){
- throw std::invalid_argument{
- "Requested process has unsupported combinations of leptons"
- };
+ 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 HEJFOG::Subleading;
set_from_yaml(name, yaml);
if(name == "none")
return none;
if(name == "all")
return all;
if(name == "unordered" || name == "uno")
return uno;
if(name == "qqx")
return qqx;
throw HEJ::unknown_option("Unknown subleading channel '"+name+"'");
}
unsigned int get_subleading_channels(YAML::Node const & node){
using YAML::NodeType;
using 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:
unsigned int channels = HEJFOG::Subleading::none;
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){
Decay decay;
set_from_yaml(decay.products, node, "into");
- set_from_yaml(decay.branching_ratio, node, "branching ratio");
+ decay.branching_ratio=1;
+ set_from_yaml_if_defined(decay.branching_ratio, node, "branching ratio");
return decay;
}
std::vector<Decay> get_decays(YAML::Node const & node){
using YAML::NodeType;
if(!node) return {};
switch(node.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)};
case NodeType::Sequence:
std::vector<Decay> result;
for(auto && decay_str: node){
- result.emplace_back();
- set_from_yaml(result.back().products, decay_str, "into");
- set_from_yaml(result.back().branching_ratio, decay_str, "branching ratio");
+ result.emplace_back(get_decay(decay_str));
}
return result;
}
throw std::logic_error{"unreachable"};
}
ParticleProperties get_particle_properties(
YAML::Node const & node, std::string const & entry
){
ParticleProperties result;
set_from_yaml(result.mass, node, entry, "mass");
set_from_yaml(result.width, node, entry, "width");
try{
result.decays = get_decays(node[entry]["decays"]);
}
catch(HEJ::missing_option const & ex){
throw HEJ::missing_option{entry + ": decays: " + ex.what()};
}
catch(HEJ::invalid_type const & ex){
throw HEJ::invalid_type{entry + ": decays: " + ex.what()};
}
return result;
}
ParticlesPropMap get_all_particles_properties(YAML::Node const & node){
ParticlesPropMap result;
for(auto const & entry: node) {
const auto name = entry.first.as<std::string>();
const auto id = HEJ::to_ParticleID(name);
result.emplace(id, get_particle_properties(node,name));
}
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){
const auto & 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_fraction > 1){
throw std::invalid_argument{
"subleading fraction has to be between 0 and 1"
};
}
if(config.subleading_fraction == 0)
config.subleading_channels = Subleading::none;
else
config.subleading_channels = get_subleading_channels(yaml["subleading channels"]);
if(!config.process.boson && config.subleading_channels != Subleading::none)
throw HEJ::not_implemented("Subleading processes for pure Jet production not implemented yet");
if(yaml["particle properties"]){
config.particles_properties = get_all_particles_properties(
yaml["particle properties"]);
}
if(config.process.boson
&& config.particles_properties.find(*(config.process.boson))
== config.particles_properties.end())
throw HEJ::missing_option("Process wants to generate boson "
+std::to_string(*(config.process.boson))+", but particle properties are missing");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
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"]) config.unweight = get_unweight(yaml, "unweight");
return config;
}
} // namespace anonymous
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;
}
}
}
diff --git a/FixedOrderGen/t/W_2j_classify.cc b/FixedOrderGen/t/W_2j_classify.cc
index 5dfc5c8..bf1d22b 100644
--- a/FixedOrderGen/t/W_2j_classify.cc
+++ b/FixedOrderGen/t/W_2j_classify.cc
@@ -1,152 +1,152 @@
/**
* \brief check that the PSP generates only "valid" W + 2 jets events
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include "JetParameters.hh"
#include "ParticleProperties.hh"
#include "PhaseSpacePoint.hh"
#include "Process.hh"
#include "Subleading.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/PDF.hh"
#include "HEJ/utility.hh"
using namespace HEJFOG;
using namespace HEJ;
namespace {
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) && !(abs(part.type)%2);
}
bool is_down_type(Particle const & part){
return HEJ::is_anyquark(part) && abs(part.type)%2;
}
bool check_W2j(PhaseSpacePoint const & psp, ParticleID const W_type){
bool found_quark = false;
bool found_anti = false;
std::vector<Particle> out_partons;
std::vector<Particle> 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(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(abs(in.type)>4 || 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 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_change = 0.5;
constexpr auto subl_channels = Subleading::all;
const ParticlesPropMap boson_prop{
{pid::Wp, {91.1876, 2.085, {Decay{ {pid::e_bar, pid::nu_e}, 1.}} }},
{pid::Wm, {91.1876, 2.085, {Decay{ {pid::e, pid::nu_e_bar}, 1.}} }}
};
HEJ::Mixmax ran{};
// Wp2j
- Process proc {{pid::proton,pid::proton}, 2, pid::Wp};
+ Process proc {{pid::proton,pid::proton}, 2, pid::Wp, {}};
size_t n_psp = n_psp_base;
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
check_W2j(psp, *proc.boson);
} else { // bad process -> 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};
+ proc = Process{{pid::proton,pid::proton}, 2, pid::Wm, {}};
n_psp = n_psp_base;
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
check_W2j(psp, *proc.boson);
} else { // bad process -> 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 9e832ab..bb95029 100644
--- a/FixedOrderGen/t/W_nj_classify.cc
+++ b/FixedOrderGen/t/W_nj_classify.cc
@@ -1,184 +1,184 @@
/**
* \brief check that the PSP generates the all W+jet subleading processes
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include "JetParameters.hh"
#include "ParticleProperties.hh"
#include "PhaseSpacePoint.hh"
#include "Process.hh"
#include "Subleading.hh"
#include "HEJ/Event.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/PDF.hh"
#include "HEJ/utility.hh"
using namespace HEJFOG;
using namespace HEJ;
namespace {
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};
}
}
int main(){
constexpr 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_change = 0.8;
const ParticlesPropMap boson_prop{
{pid::Wp, {91.1876, 2.085, {Decay{ {pid::e_bar, pid::nu_e}, 1.}} }},
{pid::Wm, {91.1876, 2.085, {Decay{ {pid::e, pid::nu_e_bar}, 1.}} }}
};
HEJ::Mixmax ran{};
auto subl_channels = Subleading::all;
std::vector<event_type::EventType> 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};
+ Process proc {{pid::proton,pid::proton}, 3, pid::Wp, {}};
size_t n_psp = n_psp_base;
std::unordered_map<event_type::EventType, size_t> type_counter;
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt) };
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
bail_out(psp, "Found not allowed event of type "
- +std::string(event_type::names[ev.type()]));
+ +std::string(event_type::name(ev.type())));
}
} else { // bad process -> 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<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
- << (event_type::names[entry.first] + std::string(":"))
+ << (event_type::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 " << event_type::names[t] << std::endl;
+ std::cerr << "Less than 5% of the events are of type " << event_type::name(t) << std::endl;
return EXIT_FAILURE;
}
}
// Wm3j - only uno
- proc = Process{{pid::proton,pid::proton}, 3, pid::Wm};
+ proc = Process{{pid::proton,pid::proton}, 3, pid::Wm, {}};
n_psp = n_psp_base;
subl_channels = Subleading::uno;
allowed_types = {event_type::FKL, event_type::unob, event_type::unof};
type_counter.clear();
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt) };
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
bail_out(psp, "Found not allowed event of type "
- +std::string(event_type::names[ev.type()]));
+ +std::string(event_type::name(ev.type())));
}
} else { // bad process -> 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<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
- << (event_type::names[entry.first] + std::string(":"))
+ << (event_type::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 " << event_type::names[t] << std::endl;
+ std::cerr << "Less than 5% of the events are of type " << event_type::name(t) << std::endl;
return EXIT_FAILURE;
}
}
// Wm4j
- proc = Process{{pid::proton,pid::proton}, 4, pid::Wm};
+ proc = Process{{pid::proton,pid::proton}, 4, pid::Wm, {}};
n_psp = n_psp_base;
subl_channels = Subleading::all;
allowed_types = {event_type::FKL,
event_type::unob, event_type::unof, event_type::qqxexb, event_type::qqxexf,
event_type::qqxmid};
type_counter.clear();
for( size_t i = 0; i<n_psp; ++i){
const PhaseSpacePoint psp{proc,jet_para,pdf,E_cms, subl_change,subl_channels,
boson_prop, ran};
if(psp.status()==good){
const Event ev{ to_EventData(psp).cluster(jet_para.def, jet_para.min_pt)};
++type_counter[ev.type()];
if( std::find(allowed_types.cbegin(), allowed_types.cend(), ev.type())
== allowed_types.cend()) {
bail_out(psp, "Found not allowed event of type "
- +std::string(event_type::names[ev.type()]));
+ +std::string(event_type::name(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<double>(entry.second)/n_psp_base;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(25)
- << (event_type::names[entry.first] + std::string(":"))
+ << (event_type::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 " << event_type::names[t] << std::endl;
+ std::cerr << "Less than 3% of the events are of type " << event_type::name(t) << std::endl;
return EXIT_FAILURE;
}
}
std::cout << "All processes passed." << std::endl;
return EXIT_SUCCESS;
}
diff --git a/FixedOrderGen/t/W_reconstruct_enu.cc b/FixedOrderGen/t/W_reconstruct_enu.cc
new file mode 100644
index 0000000..1db956a
--- /dev/null
+++ b/FixedOrderGen/t/W_reconstruct_enu.cc
@@ -0,0 +1,71 @@
+/**
+ * \brief that the reconstruction of the W works
+ *
+ * \authors The HEJ collaboration (see AUTHORS for details)
+ * \date 2019
+ * \copyright GPLv2 or later
+ */
+#ifdef NDEBUG
+#undef NDEBUG
+#endif
+
+#include <algorithm>
+
+#include "config.hh"
+#include "EventGenerator.hh"
+
+#include "HEJ/Event.hh"
+#include "HEJ/Mixmax.hh"
+
+using namespace HEJFOG;
+using namespace HEJ;
+namespace {
+ constexpr size_t num_events = 1000;
+ constexpr double invGeV2_to_pb = 389379292.;
+}
+
+double get_xs(std::string config_name){
+ auto config { load_config(config_name) };
+ config.events = num_events;
+ HEJ::Mixmax ran{1};
+ HEJFOG::EventGenerator generator{
+ 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.particles_properties,
+ config.Higgs_coupling,
+ ran
+ };
+
+ double xs = 0.;
+ for (int trials = 0; trials < config.events; ++trials){
+ auto ev = generator.gen_event();
+ if(generator.status() != good) continue;
+ assert(ev);
+ ev->central().weight *= invGeV2_to_pb;
+ ev->central().weight /= config.events;
+ xs += ev->central().weight;
+ }
+ return xs;
+}
+
+int main(){
+
+ double xs_W{ get_xs("config_Wp_2j.yml")};
+ double xs_enu{ get_xs("config_Wp_2j_decay.yml")};
+ if(std::abs(xs_W/xs_enu-1.)>1e-6){
+ std::cerr << "Reconstructing the W in the runcard gave a different results ("
+ << xs_W << " vs. "<< xs_enu << " -> " << std::abs(xs_W/xs_enu-1.)*100 << "%)\n";
+ return EXIT_FAILURE;
+ }
+
+ return EXIT_SUCCESS;
+}
diff --git a/FixedOrderGen/t/config_Wp_2j.yml b/FixedOrderGen/t/config_Wp_2j.yml
new file mode 100644
index 0000000..1f3874d
--- /dev/null
+++ b/FixedOrderGen/t/config_Wp_2j.yml
@@ -0,0 +1,27 @@
+events: 200000
+
+jets:
+ min pt: 30
+ R: 0.4
+ algorithm: antikt
+ max rapidity: 5
+
+beam:
+ energy: 6500
+ particles: [p, p]
+
+pdf: 11000
+
+process: p p => e+ nu_e 2j
+
+subleading fraction: 0.37
+
+scales: 125
+
+particle properties:
+ Wp:
+ mass: 80.385
+ width: 2.085
+
+random generator:
+ name: mixmax
diff --git a/FixedOrderGen/t/config_Wp_2j_decay.yml b/FixedOrderGen/t/config_Wp_2j_decay.yml
new file mode 100644
index 0000000..7979293
--- /dev/null
+++ b/FixedOrderGen/t/config_Wp_2j_decay.yml
@@ -0,0 +1,28 @@
+events: 200000
+
+jets:
+ min pt: 30
+ R: 0.4
+ algorithm: antikt
+ max rapidity: 5
+
+beam:
+ energy: 6500
+ particles: [p, p]
+
+pdf: 11000
+
+process: p p => Wp 2j
+
+subleading fraction: 0.37
+
+scales: 125
+
+particle properties:
+ Wp:
+ mass: 80.385
+ width: 2.085
+ decays: {into: [e+, nu_e], branching ratio: 1}
+
+random generator:
+ name: mixmax
diff --git a/doc/sphinx/HEJFOG.rst b/doc/sphinx/HEJFOG.rst
index bd0eddc..1a72676 100644
--- a/doc/sphinx/HEJFOG.rst
+++ b/doc/sphinx/HEJFOG.rst
@@ -1,301 +1,309 @@
The HEJ Fixed Order Generator
=============================
For high jet multiplicities event generation with standard fixed-order
generators becomes increasingly cumbersome. For example, the
leading-order production of a Higgs Boson with five or more jets is
computationally prohibitively expensive.
To this end, HEJ 2 provides the ``HEJFOG`` fixed-order generator
that allows to generate events with high jet multiplicities. To
facilitate the computation the limit of Multi-Regge Kinematics with
large invariant masses between all outgoing particles is assumed in the
matrix elements. The typical use of the ``HEJFOG`` is to supplement
low-multiplicity events from standard generators with high-multiplicity
events before using the HEJ 2 program to add high-energy
resummation.
Installation
------------
The ``HEJFOG`` comes bundled together with HEJ 2 and the
installation is very similar. After downloading HEJ 2 and
installing the prerequisites as described in :ref:`Installation` the
``HEJFOG`` can be installed with::
cmake /path/to/FixedOrderGen -DCMAKE_INSTALL_PREFIX=target/directory
make install
where :file:`/path/to/FixedOrderGen` refers to the :file:`FixedOrderGen`
subdirectory in the HEJ 2 directory. If HEJ 2 was
installed to a non-standard location, it may be necessary to specify the
directory containing :file:`HEJ-config.cmake`. If the base installation
directory is :file:`/path/to/HEJ`, :file:`HEJ-config.cmake` should be
found in :file:`/path/to/HEJ/lib/cmake/HEJ` and the commands for
installing the ``HEJFOG`` would read::
cmake /path/to/FixedOrderGen -DHEJ_DIR=/path/to/HEJ/lib/cmake/HEJ -DCMAKE_INSTALL_PREFIX=target/directory
make install
The installation can be tested with::
make test
provided that the CT10nlo PDF set is installed.
Running the fixed-order generator
---------------------------------
After installing the ``HEJFOG`` you can modify the provided
configuration file :file:`configFO.yml` and run the generator with::
HEJFOG configFO.yml
The resulting event file, by default named :file:`HEJFO.lhe`, can then be
fed into HEJ 2 like any event file generated from a standard
fixed-order generator, see :ref:`Running HEJ 2`.
Settings
--------
Similar to HEJ 2, the ``HEJFOG`` uses a `YAML
<http://yaml.org/>`_ configuration file. The settings are
.. _`process`:
**process**
The scattering process for which events are being generated. The
format is
:code:`in1 in2 => out1 out2 ...`
The incoming particles, :code:`in1`, :code:`in2` can be
- quarks: :code:`u`, :code:`d`, :code:`u_bar`, and so on
- gluons: :code:`g`
- protons :code:`p` or antiprotons :code:`p_bar`
+ The outgoing particles, can be
+
+ - jets: :code:`j`, multiple jets can be grouped together e.g. :code:`2j`
+ - bosons: Any gauge boson, they require `particle properties`_
+ - leptons: if they can be reconstructed to a :code:`W+` or :code:`W-`, e.g. :code:`e+ nu_e`
+
At most one of the outgoing particles can be a boson, the rest has to be
- partonic. At the moment only the Higgs boson :code:`h` is supported. All
- other outgoing particles are jets. Multiple jets can be grouped together, so
- :code:`p p => h j j` is the same as :code:`p p => h 2j`. There have to be at
- least two jets. Further decays of the boson can be added through the
- :ref:`particle properties<particle properties: particle: decays>`.
+ partonic. At the moment only Higgs :code:`h`, :code:`W+` and :code:`W-`
+ bosons are supported. All other outgoing particles are jets. There have to be
+ at least two jets. Instead of the leptons decays of the bosons can be added
+ through the :ref:`particle properties<particle properties: particle:
+ decays>`. The latter is required to decay a Higgs boson. So :code:`p p => Wp
+ j j` is the same as :code:`p p => e+ nu_e 2j`, if the decay :code:`Wp => e+
+ nu_e` is specified in `particle properties`_.
.. _`events`:
**events**
Specifies the number of events to generate.
.. _`jets`:
**jets**
Defines the properties of the generated jets.
.. _`jets: min pt`:
**min pt**
Minimum jet transverse momentum in GeV.
.. _`jets: peak pt`:
**peak pt**
Optional setting to specify the dominant jet transverse momentum
in GeV. If the generated events are used as input for HEJ
resummation, this should be set to the minimum transverse momentum
of the resummation jets. In all cases it has to be larger than
:code:`min pt`. The effect is that only a small fraction of jets
will be generated with a transverse momentum below the value of
this setting.
.. _`jets: algorithm`:
**algorithm**
The algorithm used to define jets. Allowed settings are
:code:`kt`, :code:`cambridge`, :code:`antikt`,
:code:`cambridge for passive`. See the `FastJet
<http://fastjet.fr/>`_ documentation for a description of these
algorithms.
.. _`jets: R`:
**R**
The R parameter used in the jet algorithm.
.. _`jets: max rapidity`:
**max rapidity**
Maximum absolute value of the jet rapidity.
.. _`beam`:
**beam**
Defines various properties of the collider beam.
.. _`beam: energy`:
**energy**
The beam energy in GeV. For example, the 13
TeV LHC corresponds to a value of 6500.
.. _`beam: particles`:
**particles**
A list :code:`[p1, p2]` of two beam particles. The only supported
entries are protons :code:`p` and antiprotons :code:`p_bar`.
.. _`pdf`:
**pdf**
The `LHAPDF number <https://lhapdf.hepforge.org/pdfsets>`_ of the PDF set.
For example, 230000 corresponds to an NNPDF 2.3 NLO PDF set.
.. _`subleading fraction`:
**subleading fraction**
This setting is related to the fraction of events that are not a FKL
configuration. All possible subleading process are listed in
:ref:`subleading channels<subleading channels>`. This value must be
positive and not larger than 1. It should typically be chosen between
0.01 and 0.1. Note that while this parameter influences the chance of
generating subleading configurations, it generally does not
correspond to the actual fraction of subleading events.
.. _`subleading channels`:
**subleading channels**
Optional parameters to select a specific subleading process, multiple
channels can be selected at once. If multiple subleading configurations
are possible one will be selected at random for each event, thus each
final state will include at most one subleading process, e.g. either
:code:`unordered` or :code:`qqx`. Only has an effect if
:code:`subleading fraction` is non-zero. The following values are allowed:
- :code:`all`: All channels allowed. This is the default.
- :code:`none`: No subleading contribution, only the leading (FKL)
configurations are allowed. This is equivalent to
:code:`subleading fraction: 0`.
- :code:`unordered`: Unordered emission allowed.
Unordered emission are any rapidity ordering where exactly one gluon is
emitted outside the rapidity ordering required in FKL events. More
precisely, if at least one of the incoming particles is a quark or
antiquark and there are more than two jets in the final state,
:code:`subleading fraction` states the probability that the flavours
of the outgoing particles are assigned in such a way that an unordered
configuration arises. Unordered emissions are currently not implemented
for pure jet production.
- :code:`qqx`: Production of a quark-antiquark pair. In the leading
(FKL) configurations all partons except for the most backward and
forward ones are gluons. If the :code:`qqx` channel is allowed,
:code:`subleading fraction` gives the probability of emitting a
quark-antiquark pair in place of two gluons that are adjacent in
rapidity. Quark-antiquark pairs are currently only implemented for
W boson production.
.. _`unweight`:
**unweight**
This setting defines the parameters for the partial unweighting of
events. You can disable unweighting by removing this entry from the
configuration file.
.. _`unweight: sample size`:
**sample size**
The number of weighted events used to calibrate the unweighting.
A good default is to set this to the number of target
`events`_. If the number of `events`_ is large this can
lead to significant memory consumption and a lower value should be
chosen. Contrarily, for large multiplicities the unweighting
efficiency becomes worse and the sample size should be increased.
.. _`unweight: max deviation`:
**max deviation**
Controls the range of events to which unweighting is applied. A
larger value means that a larger fraction of events are unweighted.
Typical values are between -1 and 1.
.. _`particle properties`:
**particle properties**
Specifies various properties of the different particles (Higgs, W or Z).
This is only relevant if the chosen `process`_ is the production of the
corresponding particles with jets. E.g. for the `process`_
:code:`p p => h 2j` the :code:`mass`, :code:`width` and (optionally)
:code:`decays` of the :code:`Higgs` boson are required, while all other
particle properties will be ignored.
.. _`particle properties: particle`:
**Higgs, W+, W- or Z**
The particle (Higgs, |W+|, |W-|, Z) for which the following
properties are defined.
.. |W+| replace:: W\ :sup:`+`
.. |W-| replace:: W\ :sup:`-`
.. _`particle properties: particle: mass`:
**mass**
The mass of the particle in GeV.
.. _`particle properties: particle: width`:
**width**
The total decay width of the particle in GeV.
.. _`particle properties: particle: decays`:
**decays**
Optional setting specifying the decays of the particle. Only the decay
into two particles is implemented. Each decay has the form
:code:`{into: [p1,p2], branching ratio: r}`
where :code:`p1` and :code:`p2` are the particle names of the
decay product (e.g. :code:`photon`) and :code:`r` is the branching
- ratio.
+ ratio. The branching ratio is optional (:code:`1` by default).
Decays of a Higgs boson are treated as the production and subsequent
decay of an on-shell Higgs boson, so decays into e.g. Z bosons are not
- supported.
+ supported. In contrast W bosons are decayed off-shell, thus the branching ratio should not be changed or set to :code:`1`.
.. _`scales`:
**scales**
Specifies the renormalisation and factorisation scales for the output
events. For details, see the corresponding entry in the HEJ 2
:ref:`HEJ 2 settings`. Note that this should usually be a
single value, as the weights resulting from additional scale choices
will simply be ignored when adding high-energy resummation with HEJ 2.
.. _`event output`:
**event output**
Specifies the name of a single event output file or a list of such
files. See the corresponding entry in the HEJ 2
:ref:`HEJ 2 settings` for details.
.. _`RanLux init`:
.. _`random generator`:
**random generator**
Sets parameters for random number generation. See the HEJ 2
:ref:`HEJ 2 settings` for details.
.. _`analysis`:
**analysis**
Specifies the name and settings for a custom analysis library. This
can be useful to specify cuts at the fixed-order level. See the
corresponding entry in the HEJ 2 :ref:`HEJ 2 settings`
for details.
.. _`Higgs coupling`:
**Higgs coupling**
This collects a number of settings concerning the effective coupling
of the Higgs boson to gluons. See the corresponding entry in the
HEJ 2 :ref:`HEJ 2 settings` for details.
diff --git a/include/HEJ/Constants.hh b/include/HEJ/Constants.hh
index 7b43915..30856b5 100644
--- a/include/HEJ/Constants.hh
+++ b/include/HEJ/Constants.hh
@@ -1,39 +1,39 @@
/** \file
* \brief Header file defining all global constants used for HEJ
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
namespace HEJ{
/// @name QCD parameters
//@{
constexpr double N_C = 3.; //!< number of Colours
constexpr double C_A = N_C; //!< \f$C_A\f$
constexpr double C_F = (N_C*N_C - 1.)/(2.*N_C); //!< \f$C_F\f$
constexpr double t_f = 0.5; //!< \f$t_f\f$
constexpr double n_f = 5.; //!< number light flavours
constexpr double beta0 = 11./3.*C_A - 4./3.*t_f*n_f; //!< \f$\beta_0\f$
//@}
/// @name QFT parameters
//@{
constexpr double vev = 246.2196508; //!< Higgs vacuum expectation value in GeV
- constexpr double gw = 0.653233;
- constexpr double MW = 80.419; // The W mass in GeV/c^2
- constexpr double GammaW = 2.0476; // the W width in GeV/c^2
+ constexpr double gw = 0.653233; //!< elector-weak coupling
+ constexpr double MW = 80.419; //!< The W mass in GeV/c^2
+ constexpr double GammaW = 2.0476; //!< the W width in GeV/c^2
//@}
/// @name Generation Parameters
//@{
//! Default scale for virtual correction, \f$\lambda\f$ cf. eq. (20) in \cite Andersen:2011hs
constexpr double CLAMBDA = 0.2;
constexpr double CMINPT = 0.2; //!< minimal \f$p_t\f$ of all partons
//@}
/// @name Conventional Parameters
//@{
//! Value of first colour for colour dressing, according to LHE convention \cite Boos:2001cv
constexpr int COLOUR_OFFSET = 501;
//@}
}
diff --git a/src/Event.cc b/src/Event.cc
index 63ea509..7166357 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,771 +1,785 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Event.hh"
#include <algorithm>
#include <assert.h>
#include <numeric>
#include <utility>
#include "LHEF/LHEF.h"
#include "fastjet/JetDefinition.hh"
#include "HEJ/Constants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/PDG_codes.hh"
namespace HEJ{
namespace {
constexpr int status_in = -1;
constexpr int status_decayed = 2;
constexpr int status_out = 1;
/// @name helper functions to determine event type
//@{
/**
* \brief check if final state valid for HEJ
*
* check if there is at most one photon, W, H, Z in the final state
* and all the rest are quarks or gluons
*/
bool final_state_ok(std::vector<Particle> const & outgoing){
bool has_AWZH_boson = false;
for(auto const & out: outgoing){
if(is_AWZH_boson(out.type)){
if(has_AWZH_boson) return false;
has_AWZH_boson = true;
}
else if(! is_parton(out.type)) return false;
}
return true;
}
/**
* \brief function which determines if type change is consistent with Wp emission.
* @param in incoming Particle id
* @param out outgoing Particle id
* @param qqx Current both incoming/both outgoing?
*
* \see is_Wm_Change
*/
bool is_Wp_Change(ParticleID in, ParticleID out, bool qqx){
if(!qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out== (in-1);
if( qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out==-(in+1);
return false;
}
/**
* \brief function which determines if type change is consistent with Wm emission.
* @param in incoming Particle id
* @param out outgoing Particle id
* @param qqx Current both incoming/both outgoing?
*
* Ensures that change type of quark line is possible by a flavour changing
* Wm emission. Allows checking of qqx currents also.
*/
bool is_Wm_Change(ParticleID in, ParticleID out, bool qqx){
if(!qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out== (in+1);
if( qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out==-(in-1);
return false;
}
/**
* \brief checks if particle type remains same from incoming to outgoing
* @param in incoming Particle
* @param out outgoing Particle
* @param qqx Current both incoming/outgoing?
*/
bool no_flavour_change(ParticleID in, ParticleID out, bool qqx){
const int qqxCurrent = qqx?-1:1;
if(abs(in)<=6 || in==pid::gluon) return (in==out*qqxCurrent);
else return false;
}
bool has_2_jets(Event const & event){
return event.jets().size() >= 2;
}
/**
* \brief check if we have a valid Impact factor
* @param in incoming Particle
* @param out outgoing Particle
* @param qqx Current both incoming/outgoing?
* @param qqx returns +1 if Wp, -1 if Wm, else 0
*/
bool is_valid_impact_factor(
ParticleID in, ParticleID out, bool qqx, int & W_change
){
if( no_flavour_change(in, out, qqx) ){
return true;
}
if( is_Wp_Change(in, out, qqx) ) {
W_change+=1;
return true;
}
if( is_Wm_Change(in, out, qqx) ) {
W_change-=1;
return true;
}
return false;
}
//! Returns all possible classifications from the impact factors
// the beginning points are changed s.t. after the the classification they
// point to the beginning of the (potential) FKL chain
// sets W_change: + if Wp change
// 0 if no change
// - if Wm change
// This function can be used with forward & backwards iterators
template<class OutIterator, class IndexIterator>
size_t possible_impact_factors(
ParticleID incoming_id, // incoming
OutIterator & begin_out, OutIterator const & end_out, // outgoing
IndexIterator & begin_idx, // jet indices
int & W_change, std::vector<Particle> const & boson,
bool const backward // backward?
){
using event_type::EventType;
assert(boson.size() < 2);
// keep track of all states that we don't test
size_t not_tested = EventType::qqxmid;
if(backward)
not_tested |= EventType::unof | EventType::qqxexf;
else
not_tested |= EventType::unob | EventType::qqxexb;
// Is this LL current?
if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){
++begin_out;
++begin_idx;
return not_tested | EventType::FKL;
}
// or NLL current?
// -> needs two partons in two different jets
if( std::distance(begin_out, end_out)>=2
&& *begin_idx>=0 && *(begin_idx+1)>=0 && *begin_idx!=*(begin_idx+1)
){
// Is this unordered emisson?
if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){
if( is_valid_impact_factor(
incoming_id, (begin_out+1)->type, false, W_change )
){
// veto Higgs inside uno
assert((begin_out+1)<end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < (begin_out+1)->rapidity())
||(!backward && boson.front().rapidity() > (begin_out+1)->rapidity()))
return EventType::FixedOrder;
}
begin_out+=2;
begin_idx+=2;
return not_tested | (backward?EventType::unob:EventType::unof);
}
}
// Is this QQbar?
else if( incoming_id==pid::gluon ){
if( is_valid_impact_factor(
begin_out->type, (begin_out+1)->type, true, W_change )
){
// veto Higgs inside qqx
assert((begin_out+1)<end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < (begin_out+1)->rapidity())
||(!backward && boson.front().rapidity() > (begin_out+1)->rapidity()))
return EventType::FixedOrder;
}
begin_out+=2;
begin_idx+=2;
return not_tested | (backward?EventType::qqxexb:EventType::qqxexf);
}
}
}
return EventType::FixedOrder;
}
//! Returns all possible classifications from central emissions
// the beginning points are changed s.t. after the the classification they
// point to the end of the emission chain
// sets W_change: + if Wp change
// 0 if no change
// - if Wm change
template<class OutIterator, class IndexIterator>
size_t possible_central(
OutIterator & begin_out, OutIterator const & end_out,
IndexIterator & begin_idx,
int & W_change, std::vector<Particle> const & boson
){
using event_type::EventType;
assert(boson.size() < 2);
// if we already passed the central chain,
// then it is not a valid all-order state
if(std::distance(begin_out, end_out) < 0) return EventType::FixedOrder;
// keep track of all states that we don't test
size_t possible = EventType::unob | EventType::unof
| EventType::qqxexb | EventType::qqxexf;
// Find the first non-gluon/non-FKL
while( (begin_out->type==pid::gluon) && (begin_out<end_out) ){
++begin_out;
++begin_idx;
}
// end of chain -> FKL
if( begin_out==end_out ){
return possible | EventType::FKL;
}
// is this a qqbar-pair?
// needs two partons in two separate jets
if( is_valid_impact_factor(
begin_out->type, (begin_out+1)->type, true, W_change )
&& *begin_idx>=0 && *(begin_idx+1)>=0 && *begin_idx!=*(begin_idx+1)
){
// veto Higgs inside qqx
if( !boson.empty() && boson.front().type == ParticleID::h
&& boson.front().rapidity() > begin_out->rapidity()
&& boson.front().rapidity() < (begin_out+1)->rapidity()
){
return EventType::FixedOrder;
}
begin_out+=2;
begin_idx+=2;
// remaining chain should be pure gluon/FKL
for(; begin_out<end_out; ++begin_out){
if(begin_out->type != pid::gluon) return EventType::FixedOrder;
++begin_idx;
}
return possible | EventType::qqxmid;
}
return EventType::FixedOrder;
}
/**
* \brief Checks for all event types
* @param ev Event
* @returns Event Type
*
*/
event_type::EventType classify(Event const & ev){
using event_type::EventType;
if(! has_2_jets(ev))
return EventType::no_2_jets;
// currently we can't handle multiple boson states in the ME. So they are
// considered "bad_final_state" even though the "classify" could work with
// them.
if(! final_state_ok(ev.outgoing()))
return EventType::bad_final_state;
// initialise variables
auto const & in = ev.incoming();
auto const & out = filter_partons(ev.outgoing());
auto indices{ev.particle_jet_indices({ev.jets()})};
assert(std::distance(begin(in), end(in)) == 2);
assert(out.size() >= 2);
assert(std::distance(begin(out), end(out)) >= 2);
assert(std::is_sorted(begin(out), end(out), rapidity_less{}));
auto const boson{ filter_AWZH_bosons(ev.outgoing()) };
// we only allow one boson through final_state_ok
assert(boson.size()<=1);
// keep track of potential W couplings, at the end the sum should be 0
int remaining_Wp = 0;
int remaining_Wm = 0;
if(!boson.empty() && abs(boson.front().type) == ParticleID::Wp ){
if(boson.front().type>0) ++remaining_Wp;
else ++remaining_Wm;
}
int W_change = 0;
// range for current checks
auto begin_out{out.cbegin()};
auto end_out{out.crbegin()};
auto begin_idx{indices.cbegin()};
auto end_idx{indices.crbegin()};
size_t final_type = ~(EventType::no_2_jets | EventType::bad_final_state);
// check forward impact factor
final_type &= possible_impact_factors(
in.front().type,
begin_out, end_out.base(), begin_idx,
W_change, boson, true );
assert(std::distance(begin_out, end_out.base())
== std::distance(begin_idx, end_idx.base()));
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
W_change = 0;
// check backward impact factor
final_type &= possible_impact_factors(
in.back().type,
end_out, std::make_reverse_iterator(begin_out), end_idx,
W_change, boson, false );
assert(std::distance(begin_out, end_out.base())
== std::distance(begin_idx, end_idx.base()));
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
W_change = 0;
// check central emissions
final_type &= possible_central(
begin_out, end_out.base(), begin_idx, W_change, boson );
assert(std::distance(begin_out, end_out.base())
== std::distance(begin_idx, end_idx.base()));
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
// Check whether the right number of Ws are present
if( remaining_Wp != 0 || remaining_Wm != 0 ) return EventType::FixedOrder;
// result has to be unique
if( (final_type & (final_type-1)) != 0) return EventType::FixedOrder;
return static_cast<EventType>(final_type);
}
//@}
Particle extract_particle(LHEF::HEPEUP const & hepeup, int i){
const ParticleID id = static_cast<ParticleID>(hepeup.IDUP[i]);
const fastjet::PseudoJet momentum{
hepeup.PUP[i][0], hepeup.PUP[i][1],
hepeup.PUP[i][2], hepeup.PUP[i][3]
};
if(is_parton(id))
return Particle{ id, std::move(momentum), hepeup.ICOLUP[i] };
return Particle{ id, std::move(momentum), {} };
}
bool is_decay_product(std::pair<int, int> const & mothers){
if(mothers.first == 0) return false;
return mothers.second == 0 || mothers.first == mothers.second;
}
} // namespace anonymous
Event::EventData::EventData(LHEF::HEPEUP const & hepeup){
parameters.central = EventParameters{
hepeup.scales.mur, hepeup.scales.muf, hepeup.weight()
};
size_t in_idx = 0;
for (int i = 0; i < hepeup.NUP; ++i) {
// skip decay products
// we will add them later on, but we have to ensure that
// the decayed particle is added before
if(is_decay_product(hepeup.MOTHUP[i])) continue;
auto particle = extract_particle(hepeup, i);
// needed to identify mother particles for decay products
particle.p.set_user_index(i+1);
if(hepeup.ISTUP[i] == status_in){
if(in_idx > incoming.size()) {
throw std::invalid_argument{
"Event has too many incoming particles"
};
}
incoming[in_idx++] = std::move(particle);
}
else outgoing.emplace_back(std::move(particle));
}
// add decay products
for (int i = 0; i < hepeup.NUP; ++i) {
if(!is_decay_product(hepeup.MOTHUP[i])) continue;
const int mother_id = hepeup.MOTHUP[i].first;
const auto mother = std::find_if(
begin(outgoing), end(outgoing),
[mother_id](Particle const & particle){
return particle.p.user_index() == mother_id;
}
);
if(mother == end(outgoing)){
throw std::invalid_argument{"invalid decay product parent"};
}
const int mother_idx = std::distance(begin(outgoing), mother);
assert(mother_idx >= 0);
decays[mother_idx].emplace_back(extract_particle(hepeup, i));
}
}
Event::Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double const min_jet_pt
):
Event( Event::EventData{
ev.incoming, ev.outgoing, ev.decays,
Parameters<EventParameters>{ev.central, ev.variations}
}.cluster(jet_def, min_jet_pt) )
{}
//! @TODO remove in HEJ 2.2.0
UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){
Event::EventData const evData{hepeup};
incoming = evData.incoming;
outgoing = evData.outgoing;
decays = evData.decays;
central = evData.parameters.central;
variations = evData.parameters.variations;
}
void Event::EventData::sort(){
// sort particles
std::sort(
begin(incoming), end(incoming),
[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
);
auto old_outgoing = std::move(outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){
return old_outgoing[i].rapidity() < old_outgoing[j].rapidity();
});
outgoing.clear();
outgoing.reserve(old_outgoing.size());
for(size_t i: idx) {
outgoing.emplace_back(std::move(old_outgoing[i]));
}
// find decays again
if(!decays.empty()){
auto old_decays = std::move(decays);
decays.clear();
for(size_t i=0; i<idx.size(); ++i) {
auto decay = old_decays.find(idx[i]);
if(decay != old_decays.end())
decays.emplace(i, std::move(decay->second));
}
assert(old_decays.size() == decays.size());
}
}
namespace {
Particle reconstruct_boson(std::vector<Particle> const & leptons) {
Particle decayed_boson;
decayed_boson.p = leptons[0].p + leptons[1].p;
const int pidsum = leptons[0].type + leptons[1].type;
if(pidsum == +1) {
assert(is_antilepton(leptons[0]));
if(is_antineutrino(leptons[0])) {
throw not_implemented{"lepton-flavour violating final state"};
}
assert(is_neutrino(leptons[1]));
// charged antilepton + neutrino means we had a W+
decayed_boson.type = pid::Wp;
}
else if(pidsum == -1) {
assert(is_antilepton(leptons[0]));
if(is_neutrino(leptons[1])) {
throw not_implemented{"lepton-flavour violating final state"};
}
assert(is_antineutrino(leptons[0]));
// charged lepton + antineutrino means we had a W-
decayed_boson.type = pid::Wm;
}
else {
throw not_implemented{
"final state with leptons "
+ name(leptons[0].type)
+ " and "
+ name(leptons[1].type)
};
}
return decayed_boson;
}
}
void Event::EventData::reconstruct_intermediate() {
const auto begin_leptons = std::partition(
begin(outgoing), end(outgoing),
[](Particle const & p) {return !is_anylepton(p);}
);
if(begin_leptons == end(outgoing)) return;
assert(is_anylepton(*begin_leptons));
std::vector<Particle> leptons(begin_leptons, end(outgoing));
outgoing.erase(begin_leptons, end(outgoing));
if(leptons.size() != 2) {
throw not_implemented{"Final states with one or more than two leptons"};
}
std::sort(
begin(leptons), end(leptons),
[](Particle const & p0, Particle const & p1) {
return p0.type < p1.type;
}
);
outgoing.emplace_back(reconstruct_boson(leptons));
decays.emplace(outgoing.size()-1, std::move(leptons));
}
Event Event::EventData::cluster(
fastjet::JetDefinition const & jet_def, double const min_jet_pt
){
sort();
Event ev{ std::move(incoming), std::move(outgoing), std::move(decays),
std::move(parameters),
jet_def, min_jet_pt
};
assert(std::is_sorted(begin(ev.outgoing_), end(ev.outgoing_),
rapidity_less{}));
ev.type_ = classify(ev);
return ev;
}
Event::Event(
std::array<Particle, 2> && incoming,
std::vector<Particle> && outgoing,
std::unordered_map<size_t, std::vector<Particle>> && decays,
Parameters<EventParameters> && parameters,
fastjet::JetDefinition const & jet_def,
double const min_jet_pt
): incoming_{std::move(incoming)},
outgoing_{std::move(outgoing)},
decays_{std::move(decays)},
parameters_{std::move(parameters)},
cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
min_jet_pt_{min_jet_pt}
{
jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_));
}
namespace {
void connect_incoming(Particle & in, int & colour, int & anti_colour){
in.colour = std::make_pair(anti_colour, colour);
// gluon
if(in.type == pid::gluon)
return;
if(in.type > 0){
// quark
assert(is_quark(in));
in.colour->second = 0;
colour*=-1;
return;
}
// anti-quark
assert(is_antiquark(in));
in.colour->first = 0;
anti_colour*=-1;
return;
}
}
bool Event::generate_colours(RNG & ran){
// generate only for HEJ events
if(!event_type::is_HEJ(type()))
return false;
assert(std::is_sorted(
begin(outgoing()), end(outgoing()), rapidity_less{}));
assert(incoming()[0].pz() < incoming()[1].pz());
// positive (anti-)colour -> can connect
// negative (anti-)colour -> not available/used up by (anti-)quark
int colour = COLOUR_OFFSET;
int anti_colour = colour+1;
// initialise first
connect_incoming(incoming_[0], colour, anti_colour);
for(auto & part: outgoing_){
assert(colour>0 || anti_colour>0);
if(part.type == ParticleID::gluon){
// gluon
if(colour>0 && anti_colour>0){
// on g line => connect to colour OR anti-colour (random)
if(ran.flat() < 0.5){
part.colour = std::make_pair(colour+2,colour);
colour+=2;
} else {
part.colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(colour > 0){
// on q line => connect to available colour
part.colour = std::make_pair(colour+2, colour);
colour+=2;
} else {
assert(colour<0 && anti_colour>0);
// on qx line => connect to available anti-colour
part.colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(is_quark(part)) {
// quark
assert(anti_colour>0);
if(colour>0){
// on g line => connect and remove anti-colour
part.colour = std::make_pair(anti_colour, 0);
anti_colour+=2;
anti_colour*=-1;
} else {
// on qx line => new colour
colour*=-1;
part.colour = std::make_pair(colour, 0);
}
} else if(is_antiquark(part)) {
// anti-quark
assert(colour>0);
if(anti_colour>0){
// on g line => connect and remove colour
part.colour = std::make_pair(0, colour);
colour+=2;
colour*=-1;
} else {
// on q line => new anti-colour
anti_colour*=-1;
part.colour = std::make_pair(0, anti_colour);
}
}
// else not a parton
}
// Connect last
connect_incoming(incoming_[1], anti_colour, colour);
return true;
} // generate_colours
Event::ConstPartonIterator Event::begin_partons() const {
return cbegin_partons();
};
Event::ConstPartonIterator Event::cbegin_partons() const {
return boost::make_filter_iterator(
static_cast<bool (*)(Particle const &)>(is_parton),
cbegin(outgoing()),
cend(outgoing())
);
};
Event::ConstPartonIterator Event::end_partons() const {
return cend_partons();
};
Event::ConstPartonIterator Event::cend_partons() const {
return boost::make_filter_iterator(
static_cast<bool (*)(Particle const &)>(is_parton),
cend(outgoing()),
cend(outgoing())
);
};
namespace {
void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){
const std::streamsize orig_prec = os.precision();
os <<std::scientific<<std::setprecision(6) << "["
<<std::setw(13)<<std::right<< part.px() << ", "
<<std::setw(13)<<std::right<< part.py() << ", "
<<std::setw(13)<<std::right<< part.pz() << ", "
<<std::setw(13)<<std::right<< part.E() << "]"<< std::fixed;
os.precision(orig_prec);
}
}
std::ostream& operator<<(std::ostream & os, Event const & ev){
const std::streamsize orig_prec = os.precision();
os <<std::setprecision(4)<<std::fixed;
std::cout << "########## " << event_type::name(ev.type()) << " ##########" << std::endl;
std::cout << "Incoming particles:\n";
for(auto const & in: ev.incoming()){
std::cout <<std::setw(3)<< in.type << ": ";
print_momentum(os, in.p);
std::cout << std::endl;
}
std::cout << "\nOutgoing particles: " << ev.outgoing().size() << "\n";
for(auto const & out: ev.outgoing()){
std::cout <<std::setw(3)<< out.type << ": ";
print_momentum(os, out.p);
std::cout << " => rapidity="
<<std::setw(7)<<std::right<< out.rapidity() << std::endl;
}
std::cout << "\nForming Jets: " << ev.jets().size() << "\n";
for(auto const & jet: ev.jets()){
print_momentum(os, jet);
std::cout << " => rapidity="
<<std::setw(7)<<std::right<< jet.rapidity() << std::endl;
}
+ if(ev.decays().size() > 0 ){
+ std::cout << "\nDecays: " << ev.decays().size() << "\n";
+ for(auto const & decay: ev.decays()){
+ std::cout <<std::setw(3)<< ev.outgoing()[decay.first].type
+ << " (" << decay.first << ") to:\n";
+ for(auto const & out: decay.second){
+ std::cout <<" "<<std::setw(3)<< out.type << ": ";
+ print_momentum(os, out.p);
+ std::cout << " => rapidity="
+ <<std::setw(7)<<std::right<< out.rapidity() << std::endl;
+ }
+ }
+
+ }
os << std::defaultfloat;
os.precision(orig_prec);
return os;
}
double shat(Event const & ev){
return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
}
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
LHEF::HEPEUP result;
result.heprup = heprup;
result.weights = {{event.central().weight, nullptr}};
for(auto const & var: event.variations()){
result.weights.emplace_back(var.weight, nullptr);
}
size_t num_particles = event.incoming().size() + event.outgoing().size();
for(auto const & decay: event.decays()) num_particles += decay.second.size();
result.NUP = num_particles;
// the following entries are pretty much meaningless
result.IDPRUP = event.type(); // event type
result.AQEDUP = 1./128.; // alpha_EW
//result.AQCDUP = 0.118 // alpha_QCD
// end meaningless part
result.XWGTUP = event.central().weight;
result.SCALUP = event.central().muf;
result.scales.muf = event.central().muf;
result.scales.mur = event.central().mur;
result.scales.SCALUP = event.central().muf;
result.pdfinfo.p1 = event.incoming().front().type;
result.pdfinfo.p2 = event.incoming().back().type;
result.pdfinfo.scale = event.central().muf;
result.IDUP.reserve(num_particles); // PID
result.ISTUP.reserve(num_particles); // status (in, out, decay)
result.PUP.reserve(num_particles); // momentum
result.MOTHUP.reserve(num_particles); // index mother particle
result.ICOLUP.reserve(num_particles); // colour
// incoming
for(Particle const & in: event.incoming()){
result.IDUP.emplace_back(in.type);
result.ISTUP.emplace_back(status_in);
result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
result.MOTHUP.emplace_back(0, 0);
assert(in.colour);
result.ICOLUP.emplace_back(*in.colour);
}
// outgoing
for(size_t i = 0; i < event.outgoing().size(); ++i){
Particle const & out = event.outgoing()[i];
result.IDUP.emplace_back(out.type);
const int status = event.decays().count(i)?status_decayed:status_out;
result.ISTUP.emplace_back(status);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
result.MOTHUP.emplace_back(1, 2);
if(out.colour)
result.ICOLUP.emplace_back(*out.colour);
else{
assert(is_AWZH_boson(out));
result.ICOLUP.emplace_back(std::make_pair(0,0));
}
}
// decays
for(auto const & decay: event.decays()){
for(auto const out: decay.second){
result.IDUP.emplace_back(out.type);
result.ISTUP.emplace_back(status_out);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
const size_t mother_idx = 1 + event.incoming().size() + decay.first;
result.MOTHUP.emplace_back(mother_idx, mother_idx);
result.ICOLUP.emplace_back(0,0);
}
}
assert(result.ICOLUP.size() == num_particles);
static constexpr double unknown_spin = 9.; //per Les Houches accord
result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
result.SPINUP = result.VTIMUP;
return result;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:51 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242554
Default Alt Text
(113 KB)

Event Timeline