Page MenuHomeHEPForge
Files F9501395

No OneTemporary
Actions

View Options

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/Changes-API.md b/Changes-API.md
index de1da09..ad032e0 100644
--- a/Changes-API.md
+++ b/Changes-API.md
@@ -1,66 +1,68 @@
# Changelog for HEJ API
This log lists only changes on the HEJ API. These are primarily code changes
relevant for calling HEJ as an API. This file should only be read as an addition
to `Changes.md`, where the main features are documented.
## Version 2.X
### 2.X.0
* Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public
* New class `CrossSectionAccumulator` to keep track of Cross Section of the
different subproccess
* New template struct `Parameters` similar to old `Weights`
- `Weights` are now an alias for `Parameters<double>`. Calling `Weights` did
not change
- `Weights.hh` was replaced by `Parameters.hh`. The old `Weights.hh` header
- will be removed in HEJ Version 2.3.0
+ will be removed in HEJ Version 2.2.0
* Function to multiplication and division of `EventParameters.weight` by double
- This can be combined with `Parameters`, e.g.
`Parameters<EventParameters>*Weights`, see also `Events.parameters()`
- Moved `EventParameters` to `Parameters.hh` header
* Restructured `Event` class
- `Event` can now only be build from a (new) `Event::EventData` class
- Removed default constructor for `Event`
- `Event::EventData` replaces the old `UnclusteredEvent` struct.
- `UnclusteredEvent` is now deprecated, and will be removed in HEJ Version
2.3.0
- Removed `Event.unclustered()` function
- Added new member function `Events.parameters()`, to directly access
(underlying) `Parameters<EventParameters>`
- New member functions `begin_partons`, `end_partons` with aliases
`cbegin_partons`, `cend_partons` for constant iterators over
outgoing partons.
* New function `Event::EventData.reconstruct_intermediate()` to reconstruct
bosons from decays, e.g. `positron + nu_e => Wp`
* Added optional Colour charges to particles (`Particle.colour`)
- Colour connection in the HEJ limit can be generated via
`Event.generate_colours` (automatically done in the resummation)
* New abstact `EventReader` class, as base for reading events from files
- Moved LHE file reader to `HEJ::LesHouchesReader`
- New `HEJ::HDF5Reader` to read `hdf5` files
* New function `Analysis.initialise(LHEF::HEPRUP const &)` to pass `HEPRUP` to
the analysis.
+* Renamed `EventType::nonHEJ` to `EventType::non_resummable` and `is_HEJ()`
+ to `is_resummable()` such that Run card is consistent with internal workings
## Version 2.0
### 2.0.5
* no further changes to API
### 2.0.4
* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
### 2.0.3
* no further changes to API
### 2.0.2
* no further changes to API
### 2.0.1
* no further changes to API
diff --git a/Changes.md b/Changes.md
index 176f662..0aecee0 100644
--- a/Changes.md
+++ b/Changes.md
@@ -1,54 +1,55 @@
# Changelog
This is the log for changes to the HEJ program. Further changes to the HEJ API
are documented in `Changes-API.md`. If you are using HEJ as a library, please
also read the changes there.
## Version 2.X
### 2.X.0
* Resummation for W bosons with jets
- New subleading processes `extremal qqx` & `central qqx` for a quark and
anti-quark in the final state, e.g. `g g => u d_bar Wm g` (the other
subleading processes also work with W's)
- `HEJFOG` can generate mutliple jets together with a (off-shell) W bosons
decaying into lepton & neutrino
* Allow multiplication and division of multiple scale functions e.g.
`H_T/2*m_j1j2`
* Print cross sections at end of run
* Follow HepMC convention for particle Status codes: incoming = 11,
decaying = 2, outgoing = 1 (unchanged)
* Partons now have a Colour charge
- Colours are read from and written to LHE files
- For reweighted events the colours are created according to leading colour in
the FKL limit
* Allow changing the regulator lambda in input (`regulator parameter`, only for
advanced users)
* Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`)
* Added support to read `hdf5` event files suggested in
[arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs
[HighFive](https://github.com/BlueBrain/HighFive))
* Support input with avarage weight equal to the cross section (`IDWTUP=1 or 4`)
+* Rename `non-HEJ` Processes to `non-resummable`
## 2.0.5
* Fixed event classification for input not ordered in rapidity
### 2.0.4
* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
### 2.0.3
* Fixed parsing of (numerical factor) * (base scale) in configuration
* Don't change scale names, but sanitise Rivet output file names instead
### 2.0.2
* Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
`"/"` and `"*"` before)
### 2.0.1
* Fixed name of fixed-order generator in error message.
diff --git a/config.yml b/config.yml
index f5db005..6885bca 100644
--- a/config.yml
+++ b/config.yml
@@ -1,98 +1,98 @@
# number of attempted resummation phase space points for each input event
trials: 10
min extparton pt: 30 # minimum transverse momentum of extremal partons
# maximum soft transverse momentum fraction in extremal jets
#
# max ext soft pt fraction: 0.1
resummation jets: # resummation jet properties
min pt: 35 # minimum jet transverse momentum
algorithm: antikt # jet clustering algorithm
R: 0.4 # jet R parameter
fixed order jets: # properties of input jets
min pt: 30
# by default, algorithm and R are like for resummation jets
# treatment of he various event classes
# the supported settings are: reweight, keep, discard
-# non-HEJ events cannot be reweighted
+# non-resummable events cannot be reweighted
FKL: reweight
unordered: keep
extremal qqx: keep
central qqx: keep
-non-HEJ: keep
+non-resummable: keep
# central scale choice or choices
#
# scales: [125, max jet pperp, H_T/2, 2*jet invariant mass, m_j1j2]
scales: 91.188
# factors by which the central scales should be multiplied
# renormalisation and factorisation scales are varied independently
#
# scale factors: [0.5, 0.7071, 1, 1.41421, 2]
# maximum ratio between renormalisation and factorisation scale
#
# max scale ratio: 2.0001
# import scale setting functions
#
# import scales:
# lib_my_scales.so: [scale0,scale1]
log correction: false # whether or not to include higher order logs
# event output files
#
# the supported formats are
# - Les Houches (suffix .lhe)
# - HepMC (suffix .hepmc3)
# TODO: - ROOT ntuples (suffix .root)
#
# An output file's format is deduced either automatically from the suffix
# or from an explicit specification, e.g.
# - Les Houches: outfile
event output:
- HEJ.lhe
# - HEJ_events.hepmc
# to use a rivet analysis
#
# analysis:
# rivet: MC_XS # rivet analysis name
# output: HEJ # name of the yoda files, ".yoda" and scale suffix will be added
#
# to use a custom analysis
#
# analysis:
# plugin: /path/to/libmyanalysis.so
# my analysis parameter: some value
# selection of random number generator and seed
# the choices are
# - mixmax (seed is an integer)
# - ranlux64 (seed is a filename containing parameters)
random generator:
name: mixmax
# seed: 1
# parameters for Higgs-gluon couplings
# this requires compilation with qcdloop
#
# Higgs coupling:
# use impact factors: false
# mt: 174
# include bottom: true
# mb: 4.7
## ---------------------------------------------------------------------- ##
## The following settings are only intended for advances users. ##
## Please DO NOT SET them unless you know exactly what you are doing! ##
## ---------------------------------------------------------------------- ##
#
# regulator parameter: 0.2 # The regulator lambda for the subtraction terms
diff --git a/doc/developer_manual/developer_manual.tex b/doc/developer_manual/developer_manual.tex
index 96b33fc..3c5f98c 100644
--- a/doc/developer_manual/developer_manual.tex
+++ b/doc/developer_manual/developer_manual.tex
@@ -1,1579 +1,1579 @@
\documentclass[a4paper,11pt]{article}
\usepackage{fourier}
\usepackage[T1]{fontenc}
\usepackage{microtype}
\usepackage{geometry}
\usepackage{enumitem}
\setlist[description]{leftmargin=\parindent,labelindent=\parindent}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage[utf8x]{inputenc}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{todonotes}
\usepackage{listings}
\usepackage{xspace}
\usepackage{tikz}
\usepackage{slashed}
\usepackage{subcaption}
\usetikzlibrary{arrows.meta}
\usetikzlibrary{shapes}
\usetikzlibrary{calc}
\usepackage[colorlinks,linkcolor={blue!50!black}]{hyperref}
\graphicspath{{build/figures/}{figures/}}
\emergencystretch \hsize
\newcommand{\HEJ}{{\tt HEJ}\xspace}
\newcommand{\HIGHEJ}{\emph{High Energy Jets}\xspace}
\newcommand{\cmake}{\href{https://cmake.org/}{cmake}\xspace}
\newcommand{\html}{\href{https://www.w3.org/html/}{html}\xspace}
\newcommand{\YAML}{\href{http://yaml.org/}{YAML}\xspace}
\newcommand{\QCDloop}{\href{https://github.com/scarrazza/qcdloop}{QCDloop}\xspace}
\newcommand\matel[4][]{\mathinner{\langle#2\vert#3\vert#4\rangle}_{#1}}
\newcommand{\as}{\alpha_s}
\DeclareRobustCommand{\mathgraphics}[1]{\vcenter{\hbox{\includegraphics{#1}}}}
\def\spa#1.#2{\left\langle#1\,#2\right\rangle}
\def\spb#1.#2{\left[#1\,#2\right]} \def\spaa#1.#2.#3{\langle\mskip-1mu{#1} |
#2 | {#3}\mskip-1mu\rangle} \def\spbb#1.#2.#3{[\mskip-1mu{#1} | #2 |
{#3}\mskip-1mu]} \def\spab#1.#2.#3{\langle\mskip-1mu{#1} | #2 |
{#3}\mskip-1mu\rangle} \def\spba#1.#2.#3{\langle\mskip-1mu{#1}^+ | #2 |
{#3}^+\mskip-1mu\rangle} \def\spav#1.#2.#3{\|\mskip-1mu{#1} | #2 |
{#3}\mskip-1mu\|^2} \def\jc#1.#2.#3{j^{#1}_{#2#3}}
\definecolor{darkgreen}{rgb}{0,0.4,0}
\lstset{ %
backgroundcolor=\color{lightgray}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}
basicstyle=\footnotesize\usefont{T1}{DejaVuSansMono-TLF}{m}{n}, % the size of the fonts that are used for the code
breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace
breaklines=false, % sets automatic line breaking
captionpos=t, % sets the caption-position to bottom
commentstyle=\color{red}, % comment style
deletekeywords={...}, % if you want to delete keywords from the given language
escapeinside={\%*}{*)}, % if you want to add LaTeX within your code
extendedchars=true, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8
frame=false, % adds a frame around the code
keepspaces=true, % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible)
keywordstyle=\color{blue}, % keyword style
otherkeywords={}, % if you want to add more keywords to the set
numbers=none, % where to put the line-numbers; possible values are (none, left, right)
numbersep=5pt, % how far the line-numbers are from the code
rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here))
showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces'
showstringspaces=false, % underline spaces within strings only
showtabs=false, % show tabs within strings adding particular underscores
stepnumber=2, % the step between two line-numbers. If it's 1, each line will be numbered
stringstyle=\color{gray}, % string literal style
tabsize=2, % sets default tabsize to 2 spaces
title=\lstname,
emph={},
emphstyle=\color{darkgreen}
}
\begin{document}
\tikzstyle{mynode}=[rectangle split,rectangle split parts=2, draw,rectangle split part fill={lightgray, none}]
\title{HEJ 2 developer manual}
\author{}
\maketitle
\tableofcontents
\newpage
\section{Overview}
\label{sec:overview}
HEJ 2 is a C++ program and library implementing an algorithm to
apply \HIGHEJ resummation~\cite{Andersen:2008ue,Andersen:2008gc} to
pre-generated fixed-order events. This document is intended to give an
overview over the concepts and structure of this implementation.
\subsection{Project structure}
\label{sec:project}
HEJ 2 is developed under the \href{https://git-scm.com/}{git}
version control system. The main repository is on the IPPP
\href{https://gitlab.com/}{gitlab} server under
\url{https://gitlab.dur.scotgrid.ac.uk/hej/hej}. To get a local
copy, get an account on the gitlab server and use
\begin{lstlisting}[language=sh,caption={}]
git clone git@gitlab.dur.scotgrid.ac.uk:hej/hej.git
\end{lstlisting}
This should create a directory \texttt{hej} with the following
contents:
\begin{description}
\item[doc:] Contains additional documentation, see section~\ref{sec:doc}.
\item[include:] Contains the C++ header files.
\item[src:] Contains the C++ source files.
\item[t:] Contains the source code for the automated tests.
\item[CMakeLists.txt:] Configuration file for the \cmake build
system. See section~\ref{sec:cmake}.
\item[cmake:] Auxiliary files for \cmake. This includes modules for
finding installed software in \texttt{cmake/Modules} and templates for
code generation during the build process in \texttt{cmake/Templates}.
\item[config.yml:] Sample configuration file for running HEJ 2.
\item[FixedOrderGen:] Contains the code for the fixed-order generator,
see section~\ref{sec:HEJFOG}.
\end{description}
In the following all paths are given relative to the
\texttt{hej} directory.
\subsection{Documentation}
\label{sec:doc}
The \texttt{doc} directory contains user documentation in
\texttt{doc/sphinx} and the configuration to generate source code
documentation in \texttt{doc/doxygen}.
The user documentation explains how to install and run HEJ 2. The
format is
\href{http://docutils.sourceforge.net/rst.html}{reStructuredText}, which
is mostly human-readable. Other formats, like \html, can be generated with the
\href{http://www.sphinx-doc.org/en/master/}{sphinx} generator with
\begin{lstlisting}[language=sh,caption={}]
make html
\end{lstlisting}
To document the source code we use
\href{https://www.stack.nl/~dimitri/doxygen/}{doxygen}. To generate
\html documentation, use the command
\begin{lstlisting}[language=sh,caption={}]
doxygen Doxyfile
\end{lstlisting}
in the \texttt{doc/doxygen} directory.
\subsection{Build system}
\label{sec:cmake}
For the most part, HEJ 2 is a library providing classes and
functions that can be used to add resummation to fixed-order events. In
addition, there is a relatively small executable program leveraging this
library to read in events from an input file and produce resummation
events. Both the library and the program are built and installed with
the help of \cmake.
Debug information can be turned on by using
\begin{lstlisting}[language=sh,caption={}]
cmake base/directory -DCMAKE_BUILD_TYPE=Debug
make install
\end{lstlisting}
This facilitates the use of debuggers like \href{https://www.gnu.org/software/gdb/}{gdb}.
The main \cmake configuration file is \texttt{CMakeLists.txt}. It defines the
compiler flags, software prerequisites, header and source files used to
build HEJ 2, and the automated tests.
\texttt{cmake/Modules} contains module files that help with the
detection of the software prerequisites and \texttt{cmake/Templates}
template files for the automatic generation of header and
source files. For example, this allows to only keep the version
information in one central location (\texttt{CMakeLists.txt}) and
automatically generate a header file from the template \texttt{Version.hh.in} to propagate this to the C++ code.
\subsection{General coding guidelines}
\label{sec:notes}
The goal is to make the HEJ 2 code well-structured and
readable. Here are a number of guidelines to this end.
\begin{description}
\item[Observe the boy scout rule.] Always leave the code cleaner
than how you found it. Ugly hacks can be useful for testing, but
shouldn't make their way into the main branch.
\item[Ask if something is unclear.] Often there is a good reason why
code is written the way it is. Sometimes that reason is only obvious to
the original author (use \lstinline!git blame! to find them), in which
case they should be poked to add a comment. Sometimes there is no good
reason, but nobody has had the time to come up with something better,
yet. In some places the code might just be bad.
\item[Don't break tests.] There are a number of tests in the \texttt{t}
directory, which can be run with \lstinline!make test!. Ideally, all
tests should run successfully in each git revision. If your latest
commit broke a test and you haven't pushed to the central repository
yet, you can fix it with \lstinline!git commit --amend!. If an earlier
local commit broke a test, you can use \lstinline!git rebase -i! if
you feel confident. Additionally each \lstinline!git push! is also
automatically tested via the GitLab CI (see appendix~\ref{sec:CI}).
\item[Test your new code.] When you add some new functionality, also add an
automated test. This can be useful even if you don't know the
``correct'' result because it prevents the code from changing its behaviour
silently in the future. \href{http://www.valgrind.org/}{valgrind} is a
very useful tool to detect potential memory leaks.
\item[Stick to the coding style.] It is somewhat easier to read code
that has a uniform coding and indentation style. We don't have a
strict style, but it helps if your code looks similar to what is
already there.
\end{description}
\section{Program flow}
\label{sec:flow}
A run of the HEJ 2 program has three stages: initialisation,
event processing, and cleanup. The following sections outline these
stages and their relations to the various classes and functions in the
code. Unless denoted otherwise, all classes and functions are part of
the \lstinline!HEJ! namespace. The code for the HEJ 2 program is
in \texttt{src/bin/HEJ.cc}, all other code comprises the HEJ 2
library. Classes and free functions are usually implemented in header
and source files with a corresponding name, i.e. the code for
\lstinline!MyClass! can usually be found in
\texttt{include/HEJ/MyClass.hh} and \texttt{src/MyClass.cc}.
\subsection{Initialisation}
\label{sec:init}
The first step is to load and parse the \YAML configuration file. The
entry point for this is the \lstinline!load_config! function and the
related code can be found in \texttt{include/HEJ/YAMLreader.hh},
\texttt{include/HEJ/config.hh} and the corresponding \texttt{.cc} files
in the \texttt{src} directory. The implementation is based on the
\href{https://github.com/jbeder/yaml-cpp}{yaml-cpp} library.
The \lstinline!load_config! function returns a \lstinline!Config! object
containing all settings. To detect potential mistakes as early as
possible, we throw an exception whenever one of the following errors
occurs:
\begin{itemize}
\item There is an unknown option in the \YAML file.
\item A setting is invalid, for example a string is given where a number
would be expected.
\item An option value is not set.
\end{itemize}
The third rule is sometimes relaxed for ``advanced'' settings with an
obvious default, like for importing custom scales or analyses.
The information stored in the \lstinline!Config! object is then used to
initialise various objects required for the event processing stage
described in section~\ref{sec:processing}. First, the
\lstinline!get_analysis! function creates an object that inherits from
the \lstinline!Analysis! interface.\footnote{In the context of C++ the
proper technical expression is ``pure abstract class''.} Using an
interface allows us to decide the concrete type of the analysis at run
time instead of having to make a compile-time decision. Depending on the
settings, \lstinline!get_analysis! creates either a user-defined
analysis loaded from an external library (see the user documentation
\url{https://hej.web.cern.ch/HEJ/doc/current/user/}) or the default \lstinline!EmptyAnalysis!, which does
nothing.
Together with a number of further objects, whose roles are described in
section~\ref{sec:processing}, we also initialise the global random
number generator. We again use an interface to defer deciding the
concrete type until the program is actually run. Currently, we support the
\href{https://mixmax.hepforge.org/}{MIXMAX}
(\texttt{include/HEJ/Mixmax.hh}) and Ranlux64
(\texttt{include/HEJ/Ranlux64.hh}) random number generators, both are provided
by \href{http://proj-clhep.web.cern.ch/}{CLHEP}.
We also set up a \lstinline!HEJ::EventReader! object for reading events
either in the the Les Houches event file format~\cite{Alwall:2006yp} or
an \href{https://www.hdfgroup.org/}{HDF5}-based
format~\cite{Hoeche:2019rti}. To allow making the decision at run time,
\lstinline!HEJ::EventReader! is an abstract base class defined in
\texttt{include/HEJ/EventReader.hh} and the implementations of the
derived classes are in \texttt{include/HEJ/LesHouchesReader.hh},
\texttt{include/HEJ/HDF5Reader.hh} and the corresponding \texttt{.cc}
source files in the \texttt{src} directory. The
\lstinline!LesHouchesReader! leverages
\href{http://home.thep.lu.se/~leif/LHEF/}{\texttt{include/LHEF/LHEF.h}}. A
small wrapper around the
\href{https://www.boost.org/doc/libs/1_67_0/libs/iostreams/doc/index.html}{boost
iostreams} library allows us to also read event files compressed with
\href{https://www.gnu.org/software/gzip/}{gzip}. The wrapper code is in
\texttt{include/HEJ/stream.hh} and the \texttt{src/stream.cc}.
\subsection{Event processing}
\label{sec:processing}
In the second stage events are continously read from the event
file. After jet clustering, a number of corresponding resummation events
are generated for each input event and fed into the analysis and a
number of output files. The roles of various classes and functions are
illustrated in the following flow chart:
\begin{center}
\begin{tikzpicture}[node distance=2cm and 5mm]
\node (reader) [mynode]
{\lstinline!EventReader::read_event!\nodepart{second}{read event}};
\node
(data) [mynode,below=of reader]
{\lstinline!Event::EventData! constructor\nodepart{second}{convert to \HEJ object}};
\node
(cluster) [mynode,below=of data]
{\lstinline!Event::EventData::cluster!\nodepart{second}{cluster jets \&
classify \lstinline!EventType!}};
\node
(resum) [mynode,below=of cluster]
{\lstinline!EventReweighter::reweight!\nodepart{second}{perform resummation}};
\node
(cut) [mynode,below=of resum]
{\lstinline!Analysis::pass_cuts!\nodepart{second}{apply cuts}};
\node
(fill) [mynode,below left=of cut]
{\lstinline!Analysis::fill!\nodepart{second}{analyse event}};
\node
(write) [mynode,below right=of cut]
{\lstinline!CombinedEventWriter::write!\nodepart{second}{write out event}};
\node
(control) [below=of cut] {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(reader.south) -- node[left] {\lstinline!LHEF::HEPEUP!} (data.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(data.south) -- node[left] {\lstinline!Event::EventData!} (cluster.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(cluster.south) -- node[left] {\lstinline!Event!} (resum.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(resum.south) -- (cut.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(resum.south)+(7mm, 0cm)$) -- ($(cut.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(resum.south)-(7mm, 0cm)$) -- node[left] {\lstinline!Event!} ($(cut.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)-(3mm,0mm)$) .. controls ($(control)-(3mm,0mm)$) ..node[left] {\lstinline!Event!} (fill.east);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)-(3mm,0mm)$) .. controls ($(control)-(3mm,0mm)$) .. (write.west);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)+(3mm,0mm)$) .. controls ($(control)+(3mm,0mm)$) .. (fill.east);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)+(3mm,0mm)$) .. controls ($(control)+(3mm,0mm)$) ..node[right] {\lstinline!Event!} (write.west);
\end{tikzpicture}
\end{center}
\lstinline!EventData! is an intermediate container, its members are completely
accessible. In contrast after jet clustering and classification the phase space
inside \lstinline!Event! can not be changed any more
(\href{https://wikipedia.org/wiki/Builder_pattern}{Builder design pattern}). The
resummation is performed by the \lstinline!EventReweighter! class, which is
described in more detail in section~\ref{sec:resum}. The
\lstinline!CombinedEventWriter! writes events to zero or more output files. To
this end, it contains a number of objects implementing the
\lstinline!EventWriter! interface. These event writers typically write the
events to a file in a given format. We currently have the
\lstinline!LesHouchesWriter! for event files in the Les Houches Event File
format and the \lstinline!HepMCWriter! for the
\href{https://hepmc.web.cern.ch/hepmc/}{HepMC} format (Version 2 and 3).
\subsection{Resummation}
\label{sec:resum}
In the \lstinline!EventReweighter::reweight! member function, we first
-classify the input fixed-order event (FKL, unordered, non-HEJ, \dots)
+classify the input fixed-order event (FKL, unordered, non-resummable, \dots)
and decide according to the user settings whether to discard, keep, or
resum the event. If we perform resummation for the given event, we
generate a number of trial \lstinline!PhaseSpacePoint! objects. Phase
space generation is discussed in more detail in
section~\ref{sec:pspgen}. We then perform jet clustering according to
the settings for the resummation jets on each
\lstinline!PhaseSpacePoint!, update the factorisation and
renormalisation scale in the resulting \lstinline!Event! and reweight it
according to the ratio of pdf factors and \HEJ matrix elements between
resummation and original fixed-order event:
\begin{center}
\begin{tikzpicture}[node distance=1.5cm and 5mm]
\node (in) {};
\node (treat) [diamond,draw,below=of in,minimum size=3.5cm,
label={[anchor=west, inner sep=8pt]west:discard},
label={[anchor=east, inner sep=14pt]east:keep},
label={[anchor=south, inner sep=20pt]south:reweight}
] {};
\draw (treat.north west) -- (treat.south east);
\draw (treat.north east) -- (treat.south west);
\node
(psp) [mynode,below=of treat]
{\lstinline!PhaseSpacePoint! constructor};
\node
(cluster) [mynode,below=of psp]
{\lstinline!Event::EventData::cluster!\nodepart{second}{cluster jets}};
\node
(colour) [mynode,below=of cluster]
{\lstinline!Event::generate_colours()!\nodepart{second}{generate particle colour}};
\node
(gen_scales) [mynode,below=of colour]
{\lstinline!ScaleGenerator::operator()!\nodepart{second}{update scales}};
\node
(rescale) [mynode,below=of gen_scales]
{\lstinline!PDF::pdfpt!,
\lstinline!MatrixElement!\nodepart{second}{reweight}};
\node (out) [below of=rescale] {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(in.south) -- node[left] {\lstinline!Event!} (treat.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(treat.south) -- node[left] {\lstinline!Event!} (psp.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(psp.south) -- (cluster.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(psp.south)+(7mm, 0cm)$) -- ($(cluster.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(psp.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!PhaseSpacePoint!} ($(cluster.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(cluster.south) -- (colour.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cluster.south)+(7mm, 0cm)$) -- ($(colour.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cluster.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(colour.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(colour.south) -- (gen_scales.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(colour.south)+(7mm, 0cm)$) -- ($(gen_scales.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(colour.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(gen_scales.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(gen_scales.south) -- (rescale.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(gen_scales.south)+(7mm, 0cm)$) -- ($(rescale.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(gen_scales.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(rescale.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(rescale.south) -- (out.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(rescale.south)+(7mm, 0cm)$) -- ($(out.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(rescale.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(out.north)-(7mm, 0cm)$);
\node (helper) at ($(treat.east) + (15mm,0cm)$) {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(treat.east) -- ($(treat.east) + (15mm,0cm)$)
-- node[left] {\lstinline!Event!} (helper |- gen_scales.east) -- (gen_scales.east)
;
\end{tikzpicture}
\end{center}
\subsection{Phase space point generation}
\label{sec:pspgen}
The resummed and matched \HEJ cross section for pure jet production of
FKL configurations is given by (cf. eq. (3) of~\cite{Andersen:2018tnm})
\begin{align}
\label{eq:resumdijetFKLmatched2}
% \begin{split}
\sigma&_{2j}^\mathrm{resum, match}=\sum_{f_1, f_2}\ \sum_m
\prod_{j=1}^m\left(
\int_{p_{j\perp}^B=0}^{p_{j\perp}^B=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{j\perp}^B}{(2\pi)^3}\ \int
\frac{\mathrm{d} y_j^B}{2} \right) \
(2\pi)^4\ \delta^{(2)}\!\!\left(\sum_{k=1}^{m}
\mathbf{p}_{k\perp}^B\right)\nonumber\\
&\times\ x_a^B\ f_{a, f_1}(x_a^B, Q_a^B)\ x_b^B\ f_{b, f_2}(x_b^B, Q_b^B)\
\frac{\overline{\left|\mathcal{M}_\text{LO}^{f_1f_2\to f_1g\cdots
gf_2}\big(\big\{p^B_j\big\}\big)\right|}^2}{(\hat {s}^B)^2}\nonumber\\
& \times (2\pi)^{-4+3m}\ 2^m \nonumber\\
&\times\ \sum_{n=2}^\infty\
\int_{p_{1\perp}=p_{\perp,\mathrm{min}} }^{p_{1\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{1\perp}}{(2\pi)^3}\
\int_{p_{n\perp}=p_{\perp,\mathrm{min}}}^{p_{n\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{n\perp}}{(2\pi)^3}\
\prod_{i=2}^{n-1}\int_{p_{i\perp}=\lambda}^{p_{i\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{i\perp}}{(2\pi)^3}\ (2\pi)^4\ \delta^{(2)}\!\!\left(\sum_{k=1}^n
\mathbf{p}_{k\perp}\right )\\
&\times \ \mathbf{T}_y \prod_{i=1}^n
\left(\int \frac{\mathrm{d} y_i}{2}\right)\
\mathcal{O}_{mj}^e\
\left(\prod_{l=1}^{m-1}\delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l}\perp}^B -
\mathbf{j}_{l\perp})\right)\
\left(\prod_{l=1}^m\delta(y^B_{\mathcal{J}_l}-y_{\mathcal{J}_l})\right)
\ \mathcal{O}_{2j}(\{p_i\})\nonumber\\
&\times \frac{(\hat{s}^B)^2}{\hat{s}^2}\ \frac{x_a f_{a,f_1}(x_a, Q_a)\ x_b f_{b,f_2}(x_b, Q_b)}{x_a^B\ f_{a,f_1}(x_a^B, Q_a^B)\ x_b^B\ f_{b,f_2}(x_b^B, Q_b^B)}\ \frac{\overline{\left|\mathcal{M}_{\mathrm{HEJ}}^{f_1 f_2\to f_1 g\cdots
gf_2}(\{ p_i\})\right|}^2}{\overline{\left|\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to f_1g\cdots
gf_2}\big(\big\{p^B_j\big\}\big)\right|}^{2}} \,.\nonumber
% \end{split}
\end{align}
The first two lines correspond to the generation of the fixed-order
input events with incoming partons $f_1, f_2$ and outgoing momenta
$p_j^B$, where $\mathbf{p}_{j\perp}^B$ and $y_j^B$ denote the respective
transverse momentum and rapidity. Note that, at leading order, these
coincide with the fixed-order jet momenta $p_{\mathcal{J}_j}^B$.
$f_{a,f_1}(x_a, Q_a),f_{b,f_2}(x_b, Q_b)$ are the pdf factors for the incoming partons with
momentum fractions $x_a$ and $x_b$. The square of the partonic
centre-of-mass energy is denoted by $\hat{s}^B$ and
$\mathcal{M}_\text{LO}^{f_1f_2\to f_1g\cdots gf_2}$ is the
leading-order matrix element.
The third line is a factor accounting for the different multiplicities
between fixed-order and resummation events. Lines four and five are
the integration over the resummation phase space described in this
section. $p_i$ are the momenta of the outgoing partons in resummation
phase space. $\mathbf{T}_y$ denotes rapidity
ordering and $\mathcal{O}_{mj}^e$ projects out the exclusive $m$-jet
component. The relation between resummation and fixed-order momenta is
fixed by the $\delta$ functions. The first sets each transverse fixed-order jet
momentum to some function $\mathbf{j_{l\perp}}$ of the resummation
momenta. The exact form is described in section~\ref{sec:ptj_res}. The second
$\delta$ forces the rapidities of resummation and fixed-order jets to be
the same. Finally, the last line is the reweighting of pdf and matrix
element factors already shown in section~\ref{sec:resum}.
There are two kinds of cut-off in the integration over the resummation
partons. $\lambda$ is a technical cut-off connected to the cancellation
of infrared divergencies between real and virtual corrections. Its
numerical value is set in
\texttt{include/HEJ/Constants.h}. $p_{\perp,\mathrm{min}}$ regulates
and \emph{uncancelled} divergence in the extremal parton momenta. Its
size is set by the user configuration \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJ.html#settings}.
It is straightforward to generalise eq.~(\ref{eq:resumdijetFKLmatched2})
to unordered configurations and processes with additional colourless
emissions, for example a Higgs or electroweak boson. In the latter case only
the fixed-order integration and the matrix elements change.
\subsubsection{Gluon Multiplicity}
\label{sec:psp_ng}
The first step in evaluating the resummation phase space in
eq.~(\ref{eq:resumdijetFKLmatched2}) is to randomly pick terms in the
sum over the number of emissions. This sampling of the gluon
multiplicity is done in the \lstinline!PhaseSpacePoint::sample_ng!
function in \texttt{src/PhaseSpacePoint.cc}.
The typical number of extra emissions depends strongly on the rapidity
span of the underlying fixed-order event. Let us, for example, consider
a fixed-order FKL-type multi-jet configuration with rapidities
$y_{j_f},\,y_{j_b}$ of the most forward and backward jets,
respectively. By eq.~(\ref{eq:resumdijetFKLmatched2}), the jet
multiplicity and the rapidity of each jet are conserved when adding
resummation. This implies that additional hard radiation is restricted
to rapidities $y$ within a region $y_{j_b} \lesssim y \lesssim
y_{j_f}$. Within \HEJ, we require the most forward and most backward
emissions to be hard \todo{specify how hard} in order to avoid divergences, so this constraint
in fact applies to \emph{all} additional radiation.
To simplify the remaining discussion, let us remove the FKL rapidity
ordering
\begin{equation}
\label{eq:remove_y_order}
\mathbf{T}_y \prod_{i=1}^n\int \frac{\mathrm{d}y_i}{2} =
\frac{1}{n!}\prod_{i=1}^n\int
\frac{\mathrm{d}y_i}{2}\,,
\end{equation}
where all rapidity integrals now cover a region which is approximately
bounded by $y_{j_b}$ and $y_{j_f}$. Each of the $m$ jets has to contain at least
one parton; selecting random emissions we can rewrite the phase space
integrals as
\begin{equation}
\label{eq:select_jets}
\frac{1}{n!}\prod_{i=1}^n\int [\mathrm{d}p_i] =
\left(\prod_{i=1}^{m}\int [\mathrm{d}p_i]\ {\cal J}_i(p_i)\right)
\frac{1}{n_g!}\prod_{i=m+1}^{m+n_g}\int [\mathrm{d}p_i]
\end{equation}
with jet selection functions
\begin{equation}
\label{eq:def_jet_selection}
{\cal J}_i(p) =
\begin{cases}
1 &p\text{ clustered into jet }i\\
0 & \text{otherwise}
\end{cases}
\end{equation}
and $n_g \equiv n - m$. Here and in the following we use the short-hand
notation $[\mathrm{d}p_i]$ to denote the phase-space measure for parton
$i$. As is evident from eq.~\eqref{eq:select_jets}, adding an extra emission
$n_g+1$ introduces a suppression factor $\tfrac{1}{n_g+1}$. However, the
additional phase space integral also results in an enhancement proportional
to $\Delta y_{j_f j_b} = y_{j_f} - y_{j_b}$. This is a result of the
rapidity-independence of the MRK limit of the integrand, consisting of the
matrix elements divided by the flux factor. Indeed, we observe that the
typical number of gluon emissions is to a good approximation proportional to
the rapidity separation and the phase space integral is dominated by events
with $n_g \approx \Delta y_{j_f j_b}$.
For the actual phase space sampling, we assume a Poisson distribution
and extract the mean number of gluon emissions in different rapidity
bins and fit the results to a linear function in $\Delta y_{j_f j_b}$,
finding a coefficient of $0.975$ for the inclusive production of a Higgs
boson with two jets. Here are the observed and fitted average gluon
multiplicities as a function of $\Delta y_{j_f j_b}$:
\begin{center}
\includegraphics[width=.75\textwidth]{ng_mean}
\end{center}
As shown for two rapidity slices the assumption of a Poisson
distribution is also a good approximation:
\begin{center}
\includegraphics[width=.49\textwidth]{{ng_1.5}.pdf}\hfill
\includegraphics[width=.49\textwidth]{{ng_5.5}.pdf}
\end{center}
\subsubsection{Number of Gluons inside Jets}
\label{sec:psp_ng_jet}
For each of the $n_g$ gluon emissions we can split the phase-space
integral into a (disconnected) region inside the jets and a remainder:
\begin{equation}
\label{eq:psp_split}
\int [\mathrm{d}p_i] = \int [\mathrm{d}p_i]\,
\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg) + \int [\mathrm{d}p_i]\,
\bigg[1-\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg)\bigg]\,.
\end{equation}
The next step is to decide how many of the gluons will form part of a
jet. This is done in the \lstinline!PhaseSpacePoint::sample_ng_jets!
function.
We choose an importance sampling which is flat in the plane
spanned by the azimuthal angle $\phi$ and the rapidity $y$. This is
observed in BFKL and valid in the limit of Multi-Regge-Kinematics
(MRK). Furthermore, we assume anti-$k_t$ jets, which cover an area of
$\pi R^2$.
In principle, the total accessible area in the $y$-$\phi$ plane is given
by $2\pi \Delta y_{fb}$, where $\Delta y_{fb}\geq \Delta y_{j_f j_b}$ is
the a priori unknown rapidity separation between the most forward and
backward partons. In most cases the extremal jets consist of single
partons, so that $\Delta y_{fb} = \Delta y_{j_f j_b}$. For the less common
case of two partons forming a jet we observe a maximum distance of $R$
between the constituents and the jet centre. In rare cases jets have
more than two constituents. Empirically, they are always within a
distance of $\tfrac{5}{3}R$ to the centre of the jet, so
$\Delta y_{fb} \leq \Delta y_{j_f j_b} + \tfrac{10}{3} R$. In practice, the
extremal partons are required to carry a large fraction of the jet
transverse momentum and will therefore be much closer to the jet axis.
In summary, for sufficiently large rapidity separations we can use the
approximation $\Delta y_{fb} \approx \Delta y_{j_f j_b}$. This scenario
is depicted here:
\begin{center}
\includegraphics[width=0.5\linewidth]{ps_large_y}
\end{center}
If there is no overlap between jets, the probability $p_{\cal J, >}$ for
an extra gluon to end up inside a jet is then given by
\begin{equation}
\label{eq:p_J_large}
p_{\cal J, >} = \frac{(m - 1)\*R^2}{2\Delta y_{j_f j_b}}\,.
\end{equation}
For a very small rapidity separation, eq.~\eqref{eq:p_J_large}
obviously overestimates the true probability. The maximum phase space
covered by jets in the limit of a vanishing rapidity distance between
all partons is $2mR \Delta y_{fb}$:
\begin{center}
\includegraphics[width=0.5\linewidth]{ps_small_y}
\end{center}
We therefore estimate the probability for a parton to end up inside a jet as
\begin{equation}
\label{eq:p_J}
p_{\cal J} = \min\bigg(\frac{(m - 1)\*R^2}{2\Delta y_{j_f j_b}}, \frac{mR}{\pi}\bigg)\,.
\end{equation}
Here we compare this estimate with the actually observed
fraction of additional emissions into jets as a function of the rapidity
separation:
\begin{center}
\includegraphics[width=0.75\linewidth]{pJ}
\end{center}
\subsubsection{Gluons outside Jets}
\label{sec:gluons_nonjet}
Using our estimate for the probability of a gluon to be a jet
constituent, we choose a number $n_{g,{\cal J}}$ of gluons inside
jets, which also fixes the number $n_g - n_{g,{\cal J}}$ of gluons
outside jets. As explained later on, we need to generate the momenta of
the gluons outside jets first. This is done in
\lstinline!PhaseSpacePoint::gen_non_jet!.
The azimuthal angle $\phi$ is generated flat within $0\leq \phi \leq 2
\pi$. The allowed rapidity interval is set by the most forward and
backward partons, which are necessarily inside jets. Since these parton
rapidities are not known at this point, we also have to postpone the
rapidity generation for the gluons outside jets. For the scalar
transverse momentum $p_\perp = |\mathbf{p}_\perp|$ of a gluon outside
jets we use the parametrisation
\begin{equation}
\label{eq:p_nonjet}
p_\perp = \lambda + \tilde{p}_\perp\*\tan(\tau\*r)\,, \qquad
\tau = \arctan\bigg(\frac{p_{\perp{\cal J}_\text{min}} - \lambda}{\tilde{p}_\perp}\bigg)\,.
\end{equation}
For $r \in [0,1)$, $p_\perp$ is always less than the minimum momentum
$p_{\perp{\cal J}_\text{min}}$ required for a jet. $\tilde{p}_\perp$ is
a free parameter, a good empirical value is $\tilde{p}_\perp = [1.3 +
0.2\*(n_g - n_{g,\cal J})]\,$GeV
\subsubsection{Resummation jet momenta}
\label{sec:ptj_res}
On the one hand, each jet momentum is given by the sum of its
constituent momenta. On the other hand, the resummation jet momenta are
fixed by the constraints in line five of the master
equation~\eqref{eq:resumdijetFKLmatched2}. We therefore have to
calculate the resummation jet momenta from these constraints before
generating the momenta of the gluons inside jets. This is done in
\lstinline!PhaseSpacePoint::reshuffle! and in the free
\lstinline!resummation_jet_momenta! function (declared in \texttt{resummation\_jet.hh}).
The resummation jet momenta are determined by the $\delta$ functions in
line five of eq.~(\ref{eq:resumdijetFKLmatched2}). The rapidities are
fixed to the rapidities of the jets in the input fixed-order events, so
that the FKL ordering is guaranteed to be preserved.
In traditional \HEJ reshuffling the transverse momentum are given through
\begin{equation}
\label{eq:ptreassign_old}
\mathbf{p}^B_{\mathcal{J}_{l\perp}} = \mathbf{j}_{l\perp} \equiv \mathbf{p}_{\mathcal{J}_{l}\perp}
+ \mathbf{q}_\perp \,\frac{|\mathbf{p}_{\mathcal{J}_{l}\perp}|}{P_\perp},
\end{equation}
where $\mathbf{q}_\perp = \sum_{j=1}^n \mathbf{p}_{i\perp}
\bigg[1-\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg)\bigg] $ is the
total transverse momentum of all partons \emph{outside} jets and
$P_\perp = \sum_{j=1}^m |\mathbf{p}_{\mathcal{J}_{j}\perp}|$. Since the
total transverse momentum of an event vanishes, we can also use
$\mathbf{q}_\perp = - \sum_{j=1}^m
\mathbf{p}_{\mathcal{J}_{j}\perp}$. Eq.~(\ref{eq:ptreassign}) is a
non-linear system of equations in the resummation jet momenta
$\mathbf{p}_{\mathcal{J}_{l}\perp}$. Hence we would have to solve
\begin{equation}
\label{eq:ptreassign_eq}
\mathbf{p}_{\mathcal{J}_{l}\perp}=\mathbf{j}^B_{l\perp} \equiv\mathbf{j}_{l\perp}^{-1}
\left(\mathbf{p}^B_{\mathcal{J}_{l\perp}}\right)
\end{equation}
numerically.
Since solving such a system is computationally expensive, we instead
change the reshuffling around to be linear in the resummation jet
momenta. Hence~\eqref{eq:ptreassign_eq} gets replaces by
\begin{equation}
\label{eq:ptreassign}
\mathbf{p}_{\mathcal{J}_{l\perp}} = \mathbf{j}^B_{l\perp} \equiv \mathbf{p}^B_{\mathcal{J}_{l}\perp}
- \mathbf{q}_\perp \,\frac{|\mathbf{p}^B_{\mathcal{J}_{l}\perp}|}{P^B_\perp},
\end{equation}
which is linear in the resummation momentum. Consequently the equivalent
of~\eqref{eq:ptreassign_old} is non-linear in the Born momentum. However
the exact form of~\eqref{eq:ptreassign_old} is not relevant for the resummation.
Both methods have been tested for two and three jets with the \textsc{rivet}
standard analysis \texttt{MC\_JETS}. They didn't show any differences even
after $10^9$ events.
The reshuffling relation~\eqref{eq:ptreassign} allows the transverse
momenta $p^B_{\mathcal{J}_{l\perp}}$ of the fixed-order jets to be
somewhat below the minimum transverse momentum of resummation jets. It
is crucial that this difference does not become too large, as the
fixed-order cross section diverges for vanishing transverse momenta. In
the production of a Higgs boson with resummation jets above $30\,$GeV we observe
that the contribution from fixed-order events with jets softer than
about $20\,$GeV can be safely neglected. This is shown in the following
plot of the differential cross section over the transverse momentum of
the softest fixed-order jet:
\begin{center}
\includegraphics[width=.75\textwidth]{ptBMin}
\end{center}
Finally, we have to account for the fact that the reshuffling
relation~\eqref{eq:ptreassign} is non-linear in the Born momenta. To
arrive at the master formula~\eqref{eq:resumdijetFKLmatched2} for the
cross section, we have introduced unity in the form of an integral over
the Born momenta with $\delta$ functions in the integrand, that is
\begin{equation}
\label{eq:delta_intro}
1 = \int_{p_{j\perp}^B=0}^{p_{j\perp}^B=\infty}
\mathrm{d}^2\mathbf{p}_{j\perp}^B\delta^{(2)}(\mathbf{p}_{\mathcal{J}_{j\perp}}^B -
\mathbf{j}_{j\perp})\,.
\end{equation}
If the arguments of the $\delta$ functions are not linear in the Born
momenta, we have to compensate with additional Jacobians as
factors. Explicitly, for the reshuffling relation~\eqref{eq:ptreassign}
we have
\begin{equation}
\label{eq:delta_rewrite}
\prod_{l=1}^m \delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l\perp}}^B -
\mathbf{j}_{l\perp}) = \Delta \prod_{l=1}^m \delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l\perp}} -
\mathbf{j}_{l\perp}^B)\,,
\end{equation}
where $\mathbf{j}_{l\perp}^B$ is given by~\eqref{eq:ptreassign_eq} and only
depends on the Born momenta. We have extended the product to run to $m$
instead of $m-1$ by eliminating the last $\delta$ function
$\delta^{(2)}\!\!\left(\sum_{k=1}^n \mathbf{p}_{k\perp}\right )$.
The Jacobian $\Delta$ is the determinant of a $2m \times 2m$ matrix with $l, l' = 1,\dots,m$
and $X, X' = x,y$.
\begin{equation}
\label{eq:jacobian}
\Delta = \left|\frac{\partial\,\mathbf{j}^B_{l'\perp}}{\partial\, \mathbf{p}^B_{{\cal J}_l \perp}} \right|
= \left| \delta_{l l'} \delta_{X X'} - \frac{q_X\, p^B_{{\cal
J}_{l'}X'}}{\left|\mathbf{p}^B_{{\cal J}_{l'} \perp}\right| P^B_\perp}\left(\delta_{l l'}
- \frac{\left|\mathbf{p}^B_{{\cal J}_l \perp}\right|}{P^B_\perp}\right)\right|\,.
\end{equation}
The determinant is calculated in \lstinline!resummation_jet_weight!,
again coming from the \texttt{resummation\_jet.hh} header.
Having to introduce this Jacobian is not a disadvantage specific to the new
reshuffling. If we instead use the old reshuffling
relation~\eqref{eq:ptreassign_old} we \emph{also} have to introduce a
similar Jacobian since we actually want to integrate over the
resummation phase space and need to transform the argument of the
$\delta$ function to be linear in the resummation momenta for this.
\subsubsection{Gluons inside Jets}
\label{sec:gluons_jet}
After the steps outlined in section~\ref{sec:psp_ng_jet}, we have a
total number of $m + n_{g,{\cal J}}$ constituents. In
\lstinline!PhaseSpacePoint::distribute_jet_partons! we distribute them
randomly among the jets such that each jet has at least one
constituent. We then generate their momenta in
\lstinline!PhaseSpacePoint::split! using the \lstinline!Splitter! class.
The phase space integral for a jet ${\cal J}$ is given by
\begin{equation}
\label{eq:ps_jetparton} \prod_{i\text{ in }{\cal J}} \bigg(\int
\mathrm{d}\mathbf{p}_{i\perp}\ \int \mathrm{d} y_i
\bigg)\delta^{(2)}\Big(\sum_{i\text{ in }{\cal J}} \mathbf{p}_{i\perp} -
\mathbf{j}_{\perp}^B\Big)\delta(y_{\mathcal{J}}-y^B_{\mathcal{J}})\,.
\end{equation}
For jets with a single constituent, the parton momentum is obiously equal to the
jet momentum. In the case of two constituents, we observe that the
partons are always inside the jet cone with radius $R$ and often very
close to the jet centre. The following plots show the typical relative
distance $\Delta R/R$ for this scenario:
\begin{center}
\includegraphics[width=0.45\linewidth]{dR_2}
\includegraphics[width=0.45\linewidth]{dR_2_small}
\end{center}
According to this preference for small values of $\Delta R$, we
parametrise the $\Delta R$ integrals as
\begin{equation}
\label{eq:dR_sampling}
\frac{\Delta R}{R} =
\begin{cases}
0.25\,x_R & x_R < 0.4 \\
1.5\,x_R - 0.5 & x_R \geq 0.4
\end{cases}\,.
\end{equation}
Next, we generate $\Theta_1 \equiv \Theta$ and use the constraint $\Theta_2 = \Theta
\pm \pi$. The transverse momentum of the first parton is then given by
\begin{equation}
\label{eq:delta_constraints}
p_{1\perp} =
\frac{p_{\mathcal{J} y} - \tan(\phi_2) p_{\mathcal{J} x}}{\sin(\phi_1)
- \tan(\phi_2)\cos(\phi_1)}\,.
\end{equation}
We get $p_{2\perp}$ by exchanging $1 \leftrightarrow 2$ in the
indices. To obtain the Jacobian of the transformation, we start from the
single jet phase space eq.~(\ref{eq:ps_jetparton}) with the rapidity
delta function already rewritten to be linear in the rapidity of the
last parton, i.e.
\begin{equation}
\label{eq:jet_2p}
\prod_{i=1,2} \bigg(\int
\mathrm{d}\mathbf{p}_{i\perp}\ \int \mathrm{d} y_i
\bigg)\delta^{(2)}\Big(\mathbf{p}_{1\perp} + \mathbf{p}_{2\perp} -
\mathbf{j}_{\perp}^B\Big)\delta(y_2- \dots)\,.
\end{equation}
The integral over the second parton momentum is now trivial; we can just replace
the integral over $y_2$ with the equivalent constraint
\begin{equation}
\label{eq:R2}
\int \mathrm{d}R_2 \ \delta\bigg(R_2 - \bigg[\phi_{\cal J} - \arctan
\bigg(\frac{p_{{\cal J}y} - p_{1y}}{p_{{\cal J}x} -
p_{1x}}\bigg)\bigg]/\cos \Theta\bigg) \,.
\end{equation}
In order to fix the integral over $p_{1\perp}$ instead, we rewrite this
$\delta$ function. This introduces the Jacobian
\begin{equation}
\label{eq:jac_pt1}
\bigg|\frac{\partial p_{1\perp}}{\partial R_2} \bigg| =
\frac{\cos(\Theta)\mathbf{p}_{2\perp}^2}{p_{{\cal J}\perp}\sin(\phi_{\cal J}-\phi_1)}\,.
\end{equation}
The final form of the integral over the two parton momenta is then
\begin{equation}
\label{eq:ps_jet_2p}
\int \mathrm{d}R_1\ R_1 \int \mathrm{d}R_2 \int \mathrm{d}x_\Theta\ 2\pi \int
\mathrm{d}p_{1\perp}\ p_{1\perp} \int \mathrm{d}p_{2\perp}
\ \bigg|\frac{\partial p_{1\perp}}{\partial R_2} \bigg|\delta(p_{1\perp}
-\dots) \delta(p_{2\perp} - \dots)\,.
\end{equation}
As is evident from section~\ref{sec:psp_ng_jet}, jets with three or more
constituents are rare and an efficient phase-space sampling is less
important. For such jets, we exploit the observation that partons with a
distance larger than $R_{\text{max}} = \tfrac{5}{3} R$ to
the jet centre are never clustered into the jet. Assuming $N$
constituents, we generate all components
for the first $N-1$ partons and fix the remaining parton with the
$\delta$-functional. In order to end up inside the jet, we use the
parametrisation
\begin{align}
\label{eq:ps_jet_param}
\phi_i ={}& \phi_{\cal J} + \Delta \phi_i\,, & \Delta \phi_i ={}& \Delta
R_i
\cos(\Theta_i)\,, \\
y_i ={}& y_{\cal J} + \Delta y_i\,, & \Delta y_i ={}& \Delta
R_i
\sin(\Theta_i)\,,
\end{align}
and generate $\Theta_i$ and $\Delta R_i$ randomly with $\Delta R_i \leq
R_{\text{max}}$ and the empiric value $R_{\text{max}} = 5\*R/3$. We can
then write the phase space integral for a single parton as $(p_\perp = |\mathbf{p}_\perp|)$
\begin{equation}
\label{eq:ps_jetparton_x}
\int \mathrm{d}\mathbf{p}_{\perp}\ \int
\mathrm{d} y \approx \int_{\Box} \mathrm{d}x_{\perp}
\mathrm{d}x_{ R}
\mathrm{d}x_{\theta}\
2\*\pi\,\*R_{\text{max}}^2\,\*x_{R}\,\*p_{\perp}\,\*(p_{\perp,\text{max}}
- p_{\perp,\text{min}})
\end{equation}
with
\begin{align}
\label{eq:ps_jetparton_parameters}
\Delta \phi ={}& R_{\text{max}}\*x_{R}\*\cos(2\*\pi\*x_\theta)\,,&
\Delta y ={}& R_{\text{max}}\*x_{R}\*\sin(2\*\pi\*x_\theta)\,, \\
p_{\perp} ={}& (p_{\perp,\text{max}} - p_{\perp,\text{min}})\*x_\perp +
p_{\perp,\text{min}}\,.
\end{align}
$p_{\perp,\text{max}}$ is determined from the requirement that the total
contribution from the first $n-1$ partons --- i.e. the projection onto the
jet $p_{\perp}$ axis --- must never exceed the jet $p_\perp$. This gives
\todo{This bound is too high}
\begin{equation}
\label{eq:pt_max}
p_{i\perp,\text{max}} = \frac{p_{{\cal J}\perp} - \sum_{j<i} p_{j\perp}
\cos \Delta
\phi_j}{\cos \Delta
\phi_i}\,.
\end{equation}
The $x$ and $y$ components of the last parton follow immediately from
the first $\delta$ function. The last rapidity is fixed by the condition that
the jet rapidity is kept fixed by the reshuffling, i.e.
\begin{equation}
\label{eq:yJ_delta}
y^B_{\cal J} = y_{\cal J} = \frac 1 2 \ln \frac{\sum_{i=1}^n E_i+ p_{iz}}{\sum_{i=1}^n E_i - p_{iz}}\,.
\end{equation}
With $E_n \pm p_{nz} = p_{n\perp}\exp(\pm y_n)$ this can be rewritten to
\begin{equation}
\label{eq:yn_quad_eq}
\exp(2y_{\cal J}) = \frac{\sum_{i=1}^{n-1} E_i+ p_{iz}+p_{n\perp} \exp(y_n)}{\sum_{i=1}^{n-1} E_i - p_{iz}+p_{n\perp} \exp(-y_n)}\,,
\end{equation}
which is a quadratic equation in $\exp(y_n)$. The physical solution is
\begin{align}
\label{eq:yn}
y_n ={}& \log\Big(-b + \sqrt{b^2 + \exp(2y_{\cal J})}\,\Big)\,,\\
b ={}& \bigg(\sum_{i=1}^{n-1} E_i + p_{iz} - \exp(2y_{\cal J})
\sum_{i=1}^{n-1} E_i - p_{iz}\bigg)/(2 p_{n\perp})\,.
\end{align}
\todo{what's wrong with the following?} To eliminate the remaining rapidity
integral, we transform the $\delta$ function to be linear in the
rapidity $y$ of the last parton. The corresponding Jacobian is
\begin{equation}
\label{eq:jacobian_y}
\bigg|\frac{\partial y_{\cal J}}{\partial y_n}\bigg|^{-1} = 2 \bigg( \frac{E_n +
p_{nz}}{E_{\cal J} + p_{{\cal J}z}} + \frac{E_n - p_{nz}}{E_{\cal J} -
p_{{\cal J}z}}\bigg)^{-1}\,.
\end{equation}
Finally, we check that all designated constituents are actually
clustered into the considered jet.
\subsubsection{Final steps}
\label{sec:final}
Knowing the rapidity span covered by the extremal partons, we can now
generate the rapdities for the partons outside jets. We perform jet
clustering on all partons and check in
\lstinline!PhaseSpacePoint::jets_ok! that all the following criteria are
fulfilled:
\begin{itemize}
\item The number of resummation jets must match the number of
fixed-order jets.
\item No partons designated to be outside jets may end up inside jets.
\item All other outgoing partons \emph{must} end up inside jets.
\item The extremal (in rapidity) partons must be inside the extremal
jets. If there is, for example, an unordered forward emission, the
most forward parton must end up inside the most forward jet and the
next parton must end up inside second jet.
\item The rapidities of fixed-order and resummation jets must match.
\end{itemize}
After this, we adjust the phase-space normalisation according to the
third line of eq.~(\ref{eq:resumdijetFKLmatched2}), determine the
flavours of the outgoing partons, and adopt any additional colourless
bosons from the fixed-order input event. Finally, we use momentum
conservation to reconstruct the momenta of the incoming partons.
\subsection{Colour connection}
\label{sec:Colour}
\begin{figure}
\input{src/ColourConnect.tex}
\caption{Left: Non-crossing colour flow dominating in the MRK limit. The
crossing of the colour line connecting to particle 2 can be resolved by
writing particle 2 on the left. Right: A colour flow with a (manifest)
colour-crossing. The crossing can only be resolved if one breaks the
rapidities order, e.g. switching particles 2 and 3. From~\cite{Andersen:2017sht}.}
\label{fig:Colour_crossing}
\end{figure}
After the phase space for the resummation event is generated, we can construct
the colour for each particle. To generate the colour flow one has to call
\lstinline!Event::generate_colours! on any \HEJ configuration. For non-\HEJ
event we do not change the colour, and assume it is provided by the user (e.g.
through the LHE file input). The colour connection is done in the large $N_c$
(infinite number of colour) limit with leading colour in
MRK~\cite{Andersen:2008ue, Andersen:2017sht}. The idea is to allow only
$t$-channel colour exchange, without any crossing colour lines. For example the
colour crossing in the colour connection on the left of
figure~\ref{fig:Colour_crossing} can be resolved by switching \textit{particle
2} to the left.
We can write down the colour connections by following the colour flow from
\textit{gluon a} to \textit{gluon b} and back to \textit{gluon a}, e.g.
figure~\ref{fig:Colour_gleft} corresponds to $a123ba$. In such an expression any
valid, non-crossing colour flow will connect all external legs while respecting
the rapidity ordering. Thus configurations like the left of
figure~\ref{fig:Colour_crossing} are allowed ($a134b2a$), but the right of the
same figures breaks the rapidity ordering between 2 and 3 ($a1324ba$). Note that
connections between $b$ and $a$ are in inverse order, e.g. $ab321a$ corresponds to~\ref{fig:Colour_gright} ($a123ba$) just with colour and anti-colour swapped.
\begin{figure}
\centering
\subcaptionbox{$a123ba$\label{fig:Colour_gright}}{
\includegraphics[height=0.25\textwidth]{figures/colour_gright.jpg}}
\subcaptionbox{$a13b2a$\label{fig:Colour_gleft}}{
\includegraphics[height=0.25\textwidth]{figures/colour_gleft.jpg}}
\subcaptionbox{$a\_123ba$\label{fig:Colour_qx}}{
\includegraphics[height=0.25\textwidth]{figures/colour_qx.jpg}}
\subcaptionbox{$a\_23b1a$\label{fig:Colour_uno}}{
\includegraphics[height=0.25\textwidth]{figures/colour_uno.jpg}}
\subcaptionbox{$a14b3\_2a$\label{fig:Colour_qqx}}{
\includegraphics[height=0.25\textwidth]{figures/colour_centralqqx.jpg}}
\caption{Different colour non-crossing colour connections. Both incoming
particles are drawn at the top or bottom and the outgoing left or right.
The Feynman diagram is shown in black and the colour flow in blue.}
%TODO Maybe make these plots nicer (in Latex/asy)
\end{figure}
If we replace two gluons with a quark, (anti-)quark pair we break one of the
colour connections. Still the basic concept from before holds, we just have to
treat the connection between two (anti-)quarks like an unmovable (anti-)colour.
We denote such a connection by a underscore (e.g. $1\_a$). For example the
equivalent of~\ref{fig:Colour_gright} ($a123ba$) with an incoming antiquark
is~\ref{fig:Colour_qx} ($a\_123ba$). As said this also holds for other
subleading configurations like unordered emission~\ref{fig:Colour_uno} or
central quark-antiquark pairs~\ref{fig:Colour_qqx} \footnote{Obviously this can
not be guaranteed for non-\HEJ configurations, e.g. $qQ\to Qq$ requires a
$u$-channel exchange.}.
Some rapidity ordering can have multiple possible colour connections,
e.g.~\ref{fig:Colour_gright} and~\ref{fig:Colour_gleft}. This is always the case
if a gluon radiates off a gluon line. In that case we randomly connect the gluon
to either the colour or anti-colour. Thus in the generation we keep track
whether we are on a quark or gluon line, and act accordingly.
\subsection{The matrix element }
\label{sec:ME}
The derivation of the \HEJ matrix element is explained in some detail
in~\cite{Andersen:2017kfc}, where also results for leading and
subleading matrix elements for pure multijet production and production
of a Higgs boson with at least two associated jets are listed. Matrix
elements for $Z/\gamma^*$ production together with jets are
given in~\cite{Andersen:2016vkp}, but not yet included.
A full list of all implemented currents is given in
section~\ref{sec:currents_impl}.
The matrix elements are implemented in the \lstinline!MatrixElement!
class. To discuss the structure, let us consider the squared matrix
element for FKL multijet production with $n$ final-state partons:
\begin{align}
\label{eq:ME}
\begin{split}
\overline{\left|\mathcal{M}_\text{HEJ}^{f_1 f_2 \to f_1
g\cdots g f_2}\right|}^2 = \ &\frac {(4\pi\alpha_s)^n} {4\ (N_c^2-1)}
\cdot\ \textcolor{blue}{\frac {K_{f_1}(p_1^-, p_a^-)} {t_1}\ \cdot\ \frac{K_{f_2}(p_n^+, p_b^+)}{t_{n-1}}\ \cdot\ \left\|S_{f_1 f_2\to f_1 f_2}\right\|^2}\\
& \cdot \prod_{i=1}^{n-2} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
& \cdot \prod_{j=1}^{n-1} \textcolor{red}{\exp\left[\omega^0(q_{j\perp})(y_{j+1}-y_j)\right]}.
\end{split}
\end{align}
The structure and momentum assignment of the unsquared matrix element is
as illustrated here:
\begin{center}
\includegraphics{HEJ_amplitude}
\end{center}
The square
of the complete matrix element as given in eq.~\eqref{eq:ME} is
calculated by \lstinline!MatrixElement::operator()!. The \textcolor{red}{last line} of
eq.~\eqref{eq:ME} constitutes the all-order virtual correction,
implemented in
\lstinline!MatrixElement::virtual_corrections!.
$\omega^0$ is the
\textit{regularised Regge trajectory}
\begin{equation}
\label{eq:omega_0}
\omega^0(q_\perp) = - C_A \frac{\alpha_s}{\pi} \log \left(\frac{q_\perp^2}{\lambda^2}\right)\,,
\end{equation}
where $\lambda$ is the slicing parameter limiting the softness of real
gluon emissions, cf. eq.~\eqref{eq:resumdijetFKLmatched2}. $\lambda$ can be
changed at runtime by setting \lstinline!regulator parameter! in
\lstinline!conifg.yml!.
The remaining parts, which correspond to the square of the leading-order
HEJ matrix element $\overline{\left|\mathcal{M}_\text{LO,
HEJ}^{f_1f_2\to f_1g\cdots
gf_2}\big(\big\{p^B_j\big\}\big)\right|}^{2}$, are computed in
\lstinline!MatrixElement::tree!. We can further factor off the
scale-dependent ``parametric'' part
\lstinline!MatrixElement::tree_param! containing all factors of the
strong coupling $4\pi\alpha_s$. Using this function saves some CPU time
when adjusting the renormalisation scale, see
section~\ref{sec:resum}. The remaining ``kinematic'' factors are
calculated in \lstinline!MatrixElement::kin!.
\subsubsection{Matrix elements for Higgs plus dijet}
\label{sec:ME_h_jets}
In the production of a Higgs boson together with jets the parametric
parts and the virtual corrections only require minor changes in the
respective functions. However, in the ``kinematic'' parts we have to
distinguish between several cases, which is done in
\lstinline!MatrixElement::tree_kin_Higgs!. The Higgs boson can be
\emph{central}, i.e. inside the rapidity range spanned by the extremal
partons (\lstinline!MatrixElement::tree_kin_Higgs_central!) or
\emph{peripheral} and outside this range
(\lstinline!MatrixElement::tree_kin_Higgs_first! or
\lstinline!MatrixElement::tree_kin_Higgs_last!). Currently the current for an
unordered emission with an Higgs on the same side it not implemented
\footnote{In principle emitting a Higgs boson \textit{on the other
side} of the unordered gluon is possible by contracting an unordered and
external Higgs current. Obviously this would not cover all possible
configurations, e.g. $qQ\to HgqQ$ requires contraction of the standard $Q\to Q$
current with an (unknown) $q\to Hgq$ one.}.
If a Higgs boson with momentum $p_H$ is emitted centrally, after parton
$j$ in rapidity, the matrix element reads
\begin{equation}
\label{eq:ME_h_jets_central}
\begin{split}
\overline{\left|\mathcal{M}_\text{HEJ}^{f_1 f_2 \to f_1 g\cdot H
\cdot g f_2}\right|}^2 = \ &\frac {\alpha_s^2 (4\pi\alpha_s)^n} {4\ (N_c^2-1)}
\cdot\ \textcolor{blue}{\frac {K_{f_1}(p_1^-, p_a^-)} {t_1}\
\cdot\ \frac{1}{t_j t_{j+1}} \cdot\ \frac{K_{f_2}(p_n^+, p_b^+)}{t_{n}}\ \cdot\ \left\|S_{f_1
f_2\to f_1 H f_2}\right\|^2}\\
& \cdot \prod_{\substack{i=1\\i \neq j}}^{n-1} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
& \cdot \textcolor{red}{\prod_{i=1}^{n-1}
\exp\left[\omega^0(q_{i\perp})\Delta y_i\right]}
\end{split}
\end{equation}
with the momentum definitions
\begin{center}
\includegraphics{HEJ_central_Higgs_amplitude}
\end{center}
$q_i$ is the $i$th $t$-channel momentum and $\Delta y_i$ the rapidity
gap between outgoing \emph{particles} (not partons) $i$ and $i+1$ in
rapidity ordering.
For \emph{peripheral} emission in the backward direction
(\lstinline!MatrixElement::tree_kin_Higgs_first!) we first check whether
the most backward parton is a gluon or an (anti-)quark. In the latter
case the leading contribution to the matrix element arises through
emission off the $t$-channel gluons and we can use the same formula
eq.~(\ref{eq:ME_h_jets_central}) as for central emission. If the most
backward parton is a gluon, the square of the matrix element can be
written as
\begin{equation}
\label{eq:ME_h_jets_peripheral}
\begin{split}
\overline{\left|\mathcal{M}_\text{HEJ}^{g f_2 \to H g\cdot g f_2}\right|}^2 = \ &\frac {\alpha_s^2 (4\pi\alpha_s)^n} {\textcolor{blue}{4\ (N_c^2-1)}}
\textcolor{blue}{\cdot\ K_{H}\
\frac{K_{f_2}(p_n^+, p_b^+)}{t_{n-1}}\ \cdot\ \left\|S_{g
f_2\to H g f_2}\right\|^2}\\
& \cdot \prod_{\substack{i=1}}^{n-2} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
& \cdot \textcolor{red}{\prod_{i=1}^{n-1}
\exp\left[\omega^0(q_{i\perp}) (y_{i+1} - y_i)\right]}
\end{split}
\end{equation}
with the momenta as follows:
\begin{center}
\includegraphics{HEJ_peripheral_Higgs_amplitude}
\end{center}
The \textcolor{blue}{blue part} is implemented in
\lstinline!MatrixElement::MH2_forwardH!. All other building blocks are
already available.\todo{Impact factors} The actual current contraction
is calculated in \lstinline!MH2gq_outsideH! inside
\lstinline!currents.cc!, which corresponds to $\tfrac{16 \pi^2}{t_1} \left\|S_{g
f_2\to H g f_2}\right\|^2$.\todo{Fix this insane normalisation}
The forward emission of a Higgs boson is completely analogous. We can
use the same function \lstinline!MatrixElement::MH2_forwardH!, swapping
$p_1 \leftrightarrow p_n,\,p_a \leftrightarrow p_b$.
\subsubsection{FKL ladder and Lipatov vertices}
\label{sec:FKL_ladder}
The ``FKL ladder'' is the product
\begin{equation}
\label{eq:FKL_ladder}
\prod_{i=1}^{n-2} \left( \frac{-C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)
\end{equation}
appearing in the square of the matrix element for $n$ parton production,
cf. eq.~(\ref{eq:ME}), and implemented in
\lstinline!MatrixElement::FKL_ladder_weight!. The Lipatov vertex contraction
$V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1})$ is implemented \lstinline!C2Lipatovots!.
It is given by \todo{equation} \todo{mention difference between the two versions
of \lstinline!C2Lipatovots!, maybe even get rid of one}.
\subsubsection{Currents}
\label{sec:currents}
The current factors $\frac{K_{f_1}K_{f_2}}{t_1 t_{n-1}}\left\|S_{f_1
f_2\to f_1 f_2}\right\|^2$ and their extensions for unordered and Higgs
boson emissions are implemented in the \lstinline!jM2!$\dots$ functions
of \texttt{src/currents.cc}. \todo{Only $\|S\|^2$ should be in currents}
\footnote{The current implementation for
Higgs production in \texttt{src/currents.cc} includes the $1/4$ factor
inside $S$, opposing to~\eqref{eq:ME}. Thus the overall normalisation is
unaffected.} The ``colour acceleration multiplier'' (CAM) $K_{f}$
for a parton $f\in\{g,q,\bar{q}\}$ is defined as
\begin{align}
\label{eq:K_g}
K_g(p_1^-, p_a^-) ={}& \frac{1}{2}\left(\frac{p_1^-}{p_a^-} + \frac{p_a^-}{p_1^-}\right)\left(C_A -
\frac{1}{C_A}\right)+\frac{1}{C_A}\\
\label{eq:K_q}
K_q(p_1^-, p_a^-) ={}&K_{\bar{q}}(p_1^-, p_a^-) = C_F\,.
\end{align}
The Higgs current CAM used in eq.~(\ref{eq:ME_h_jets_peripheral}) is
\begin{equation}
\label{eq:K_H}
K_H = C_A\,.
\end{equation}
The current contractions are given by\todo{check all this
carefully!}
\begin{align}
\label{eq:S}
\left\|S_{f_1 f_2\to f_1 f_2}\right\|^2 ={}& \sum_{\substack{\lambda_a =
+,-\\\lambda_b = +,-}} \left|j^{\lambda_a}_\mu(p_1, p_a)\
j^{\lambda_b\,\mu}(p_n, p_b)\right|^2 = 2\sum_{\lambda =
+,-} \left|j^{-}_\mu(p_1, p_a)\ j^{\lambda\,\mu}(p_n, p_b)\right|^2\,,\\
\left\|S_{f_1 f_2\to f_1 H f_2}\right\|^2 ={}& \sum_{\substack{\lambda_a =
+,-\\\lambda_b = +,-}} \left|j^{\lambda_a}_\mu(p_1, p_a)V_H^{\mu\nu}(q_j, q_{j+1})\
j^{\lambda_b}_\nu(p_n, p_b)\right|^2\,,\\
\left\|S_{g f_2 \to H g f_2}\right\|^2 ={}& \sum_{
\substack{
\lambda_{a} = +,-\\
\lambda_{1} =+,-\\
\lambda_{b} = +,-
}}
\left|j^{\lambda_a\lambda_1}_{H\,\mu}(p_1, p_a, p_H)\ j^{\lambda_b\,\mu}(p_n, p_b)\right|^2\,.
\end{align}
The ``basic'' currents $j$ are independent of the parton flavour and read
\begin{equation}
\label{eq:j}
j^\pm_\mu(p, q) = u^{\pm,\dagger}(p)\ \sigma^\pm_\mu\ u^{\pm}(q)\,,
\end{equation}
where $\sigma_\mu^\pm = (1, \pm \sigma_i)$ and $\sigma_i$ are the Pauli
matrices
\begin{equation}
\label{eq:Pauli_matrices}
\sigma_1 =
\begin{pmatrix}
0 & 1\\ 1 & 0
\end{pmatrix}
\,,
\qquad \sigma_2 =
\begin{pmatrix}
0 & -i\\ i & 0
\end{pmatrix}
\,,
\qquad \sigma_3 =
\begin{pmatrix}
1 & 0\\ 0 & -1
\end{pmatrix}
\,.
\end{equation}
The two-component chiral spinors are given by
\begin{align}
\label{eq:u_plus}
u^+(p)={}& \left(\sqrt{p^+}, \sqrt{p^-} \hat{p}_\perp \right) \,,\\
\label{eq:u_minus}
u^-(p)={}& \left(\sqrt{p^-} \hat{p}^*_\perp, -\sqrt{p^+}\right)\,,
\end{align}
with $p^\pm = E\pm p_z,\, \hat{p}_\perp = \tfrac{p_\perp}{|p_\perp|},\,
p_\perp = p_x + i p_y$. The spinors for vanishing transverse momentum
are obtained by replacing $\hat{p_\perp} \to -1$.
Explicitly, the currents read
\begin{align}
\label{eq:j-_explicit}
j^-_\mu(p, q) ={}&
\begin{pmatrix}
\sqrt{p^+\,q^+} + \sqrt{p^-\,q^-} \hat{p}_{\perp} \hat{q}_{\perp}^*\\
\sqrt{p^-\,q^+}\, \hat{p}_{\perp} + \sqrt{p^+\,q^-}\,\hat{q}_{\perp}^*\\
-i \sqrt{p^-\,q^+}\, \hat{p}_{\perp} + i \sqrt{p^+\,q^-}\, \hat{q}_{\perp}^*\\
\sqrt{p^+\,q^+} - \sqrt{p^-\,q^-}\, \hat{p}_{\perp}\, \hat{q}_{\perp}^*
\end{pmatrix}\\
j^+_\mu(p, q) ={}&\big(j^-_\mu(p, q)\big)^*
\end{align}
If $q= p_{\text{in}}$ is the momentum of an incoming parton, we have
$\hat{p}_{\text{in} \perp} = -1$ and either $p_{\text{in}}^+ = 0$ or
$p_{\text{in}}^- = 0$. The current simplifies further:\todo{Helicities flipped w.r.t code}
\begin{align}
\label{eq:j_explicit}
j^-_\mu(p_{\text{out}}, p_{\text{in}}) ={}&
\begin{pmatrix}
\sqrt{p_{\text{in}}^+\,p_{\text{out}}^+}\\
\sqrt{p_{\text{in}}^+\,p_{\text{out}}^-} \ \hat{p}_{\text{out}\,\perp}\\
-i\,j^-_1\\
j^-_0
\end{pmatrix}
& p_{\text{in}\,z} > 0\,,\\
j^-_\mu(p_{\text{out}}, p_{\text{in}}) ={}&
\begin{pmatrix}
-\sqrt{p_{\text{in}}^-\,p_{\text{out}}^{-\phantom{+}}} \ \hat{p}_{\text{out}\,\perp}\\
- \sqrt{p_{\text{in}}^-\,p_{\text{out}}^+}\\
i\,j^-_1\\
-j^-_0
\end{pmatrix} & p_{\text{in}\,z} < 0\,.
\end{align}
\section{The fixed-order generator}
\label{sec:HEJFOG}
Even at leading order, standard fixed-order generators can only generate
events with a limited number of final-state particles within reasonable
CPU time. The purpose of the fixed-order generator is to supplement this
with high-multiplicity input events according to the first two lines of
eq.~\eqref{eq:resumdijetFKLmatched2} with the \HEJ approximation
$\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to f_1g\cdots gf_2}$ instead of the
full fixed-order matrix element $\mathcal{M}_\text{LO}^{f_1f_2\to
f_1g\cdots gf_2}$. Its usage is described in the user
documentation \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJFOG.html}.
\subsection{File structure}
\label{sec:HEJFOG_structure}
The code for the fixed-order generator is in the \texttt{FixedOrderGen}
directory, which contains the following:
\begin{description}
\item[include:] Contains the C++ header files.
\item[src:] Contains the C++ source files.
\item[t:] Contains the source code for the automated tests.
\item[CMakeLists.txt:] Configuration file for the \cmake build system.
\item[configFO.yml:] Sample configuration file for the fixed-order generator.
\end{description}
The code is generally in the \lstinline!HEJFOG! namespace. Functions and
classes \lstinline!MyClass! are usually declared in
\texttt{include/MyClass.hh} and implemented in \texttt{src/MyClass.cc}.
\subsection{Program flow}
\label{sec:prog_flow}
A single run of the fixed-order generator consists of three or four
stages.
First, we perform initialisation similar to HEJ 2, see
section~\ref{sec:init}. Since there is a lot of overlap we frequently
reuse classes and functions from HEJ 2, i.e. from the
\lstinline!HEJ! namespace. The code for parsing the configuration file
is in \texttt{include/config.hh} and implemented in
\texttt{src/config.cc}.
If partial unweighting is requested in the user settings \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJFOG.html#settings},
the initialisation is followed by a calibration phase. We use a
\lstinline!EventGenerator! to produce a number of trial
events. We use these to calibrate the \lstinline!Unweighter! in
its constructor and produce a first batch of partially unweighted
events. This also allows us to estimate our unweighting efficiency.
In the next step, we continue to generate events and potentially
unweight them. Once the user-defined target number of events is reached,
we adjust their weights according to the number of required trials. As
in HEJ 2 (see section~\ref{sec:processing}), we pass the final
events to a \lstinline!HEJ::Analysis! and a
\lstinline!HEJ::CombinedEventWriter!.
\subsection{Event generation}
\label{sec:evgen}
Event generation is performed by the
\lstinline!EventGenerator::gen_event! member function. We begin by generating a
\lstinline!PhaseSpacePoint!. This is not to be confused with
the resummation phase space points represented by
\lstinline!HEJ::PhaseSpacePoint!! After jet clustering, we compute the
leading-order matrix element (see section~\ref{sec:ME}) and pdf factors.
The phase space point generation is performed in the
\lstinline!PhaseSpacePoint! constructor. We first construct the
user-defined number of $n_p$ partons (by default gluons) in
\lstinline!PhaseSpacePoint::gen_LO_partons!. We use flat sampling in
rapidity and azimuthal angle. For the scalar transverse momenta, we
distinguish between two cases. By default, they are generated based on a
random variable $x_{p_\perp}$ according to
\begin{equation}
\label{eq:pt_sampling}
p_\perp = p_{\perp,\text{min}} +
\begin{cases}
p_{\perp,\text{par}}
\tan\left(
x_{p_\perp}
\arctan\left(
\frac{p_{\perp,\text{max}} - p_{\perp,\text{min}}}{p_{\perp,\text{par}}}
\right)
\right)
& y < y_\text{cut}
\\
- \tilde{p}_{\perp,\text{par}}\log\left(1 - x_{p_\perp}\left[1 -
\exp\left(\frac{p_{\perp,\text{min}} -
p_{\perp,\text{max}}}{\tilde{p}_{\perp,\text{par}}}\right)\right]\right)
& y \geq y_\text{cut}
\end{cases}\,,
\end{equation}
where $p_{\perp,\text{min}}$ is the minimum jet transverse momentum,
$p_{\perp,\text{max}}$ is the maximum transverse parton momentum,
tentatively set to the beam energy, and $y_\text{cut}$, $p_{\perp,\text{par}}$
and $\tilde{p}_{\perp,\text{par}}$ are generation parameters set to
heuristically determined values of
\begin{align}
y_\text{cut}&=3,\\
p_{\perp,\text{par}}&=p_{\perp,\min}+\frac{n_p}{5}, \\
\tilde{p}_{\perp,\text{par}}&=\frac{p_{\perp,\text{par}}}{1 +
5(y-y_\text{cut})}.
\end{align}
The problem with this generation is that the transverse momenta peak at
the minimum transverse momentum required for fixed-order jets. However,
if we use the generated events as input for \HEJ resummation, events
with such soft transverse momenta hardly contribute, see
section~\ref{sec:ptj_res}. To generate efficient input for resummation,
there is the user option \texttt{peak pt}, which specifies the
dominant transverse momentum for resummation jets. If this option is
set, most jets will be generated as above, but with
$p_{\perp,\text{min}}$ set to the peak transverse momentum $p_{\perp,
\text{peak}}$. In addition, there is a small chance of around $2\%$ to
generate softer jets. The heuristic ansatz for the transverse momentum
distribution in the ``soft'' region is
\begin{equation}
\label{FO_pt_soft}
\frac{\partial \sigma}{\partial p_\perp} \propto e^{n_p\frac{p_\perp- p_{\perp,
\text{peak}}}{\bar{p}_\perp}}\,,
\end{equation}
where $n_p$ is the number of partons and $\bar{p}_\perp \approx
4\,$GeV. To achieve this distribution, we use
\begin{equation}
\label{eq:FO_pt_soft_sampling}
p_\perp = p_{\perp, \text{peak}} + \bar{p}_\perp \frac{\log x_{p_\perp}}{n_p}
\end{equation}
and discard the phase space point if the parton is too soft, i.e. below the threshold for
fixed-order jets.
After ensuring that all partons form separate jets, we generate any
potential colourless emissions. We then determine the incoming momenta
and flavours in \lstinline!PhaseSpacePoint::reconstruct_incoming! and
adjust the outgoing flavours to ensure an FKL configuration. Finally, we
may reassign outgoing flavours to generate suppressed (for example
unordered) configurations.
\subsection{Unweighting}
\label{sec:unweight}
Straightforward event generation tends to produce many events with small
weights. Those events have a negligible contribution to the final
observables, but can take up considerable storage space and CPU time in
later processing stages. This problem can be addressed by unweighting.
For naive unweighting, one would determine the maximum weight
$w_\text{max}$ of all events, discard each event with weight $w$ with a
probability $p=w/w_\text{max}$, and set the weights of all remaining
events to $w_\text{max}$. The downside to this procedure is that it also
eliminates a sizeable fraction of events with moderate weight, so that
the statistical convergence deteriorates.
To ameliorate this problem, we perform unweighting only for events with
sufficiently small weights. This is done by the
\lstinline!Unweighter! class. In the constructor we estimate the
mean and width of the weight-weight distribution from a sample of
events. We use these estimates to determine the maximum weight below
which unweighting is performed. The actual unweighting is the done in
the \lstinline!Unweighter::unweight! function.
\input{currents}
\input{tensor}
\appendix
\section{Continuous Integration}
\label{sec:CI}
If you are implementing something new or fixed a bug please also add a test to
the \texttt{t/} folder and add it to the main \lstinline!CMakeLists.txt! via
\lstinline!add_test!. These test can be triggered by running
\lstinline!make test! or \lstinline!ctest! after compiling. A typical test
should be at most a few seconds, so it can be potentially run on each commit
change by each developer. If you require a longer, more careful test, preferably
on top of a small one, surround it with
\begin{lstlisting}[caption={}]
if(${TEST_ALL})
add_test(
NAME t_feature
COMMAND really_long_test
)
endif()
\end{lstlisting}
Afterwards you can execute the longer tests with\footnote{No recompiling is
needed, as long as only the \lstinline!add_test! command is guarded, not the
compiling commands itself.}
\begin{lstlisting}[language=sh,caption={}]
cmake base/directory -DTEST_ALL=TRUE
make test
\end{lstlisting}
On top of that you should add
\href{https://en.cppreference.com/w/cpp/error/assert}{\lstinline!assert!s} in
the code itself. They are only executed when compiled with
\lstinline!CMAKE_BUILD_TYPE=Debug!, without slowing down release code. So you
can use them everywhere to test \textit{expected} or \textit{assumed} behaviour,
e.g. requiring a Higgs boson or relying on rapidity ordering.
GitLab provides ways to directly test code via \textit{Continuous integrations}.
The CI is controlled by \texttt{.gitlab-ci.yml}. For all options for the YAML
file see \href{https://docs.gitlab.com/ee/ci/yaml/}{docs.gitlab.com/ee/ci/yaml/}.
GitLab also provides a small tool to check that YAML syntax is correct under
\lstinline!CI/CD > Pipelines > CI Lint! or
\href{https://gitlab.dur.scotgrid.ac.uk/hej/HEJ/-/ci/lint}{gitlab.dur.scotgrid.ac.uk/hej/HEJ/-/ci/lint}.
Currently the CI is configured to trigger a \textit{Pipeline} on each
\lstinline!git push!. The corresponding \textit{GitLab runners} are configured
under \lstinline!CI/CD Settings>Runners! in the GitLab UI. All runners use a
\href{https://www.docker.com/}{docker} image as virtual environments\footnote{To
use only Docker runners set the \lstinline!docker! tag in
\texttt{.gitlab-ci.yml}.}. The specific docker images maintained separately. If
you add a new dependences, please also provide a docker image for the CI. The
goal to be able to test \HEJ with all possible configurations.
Each pipeline contains multiple stages (see \lstinline!stages! in
\texttt{.gitlab-ci.yml}) which are executed in order from top to bottom.
Additionally each stage contains multiple jobs. For example the stage
\textit{build} contains the jobs \lstinline!build:basic!,
\lstinline!build:qcdloop!, \lstinline!build:rivet!, etc., which compile \HEJ for
different environments and dependences, by using different in the Docker images.
Jobs staring with an dot are ignored by the Runner, e.g. \lstinline!.HEJ_build!
is only used as a template but never executed directly. Only after all jobs of
the previous stage was executed without any error the next stage will start.
To pass information between one stage and the next we use \lstinline!artifacts!.
The runner will automatically load all artifacts form all
\lstinline!dependencies! for each job\footnote{If no dependencies are defined
\textit{all} artifacts from all previous jobs are downloaded. Thus please
specify an empty dependence if you do not want to load any artifacts.}. For
example the compiled \HEJ code from \lstinline!build:basic! gets loaded in
\lstinline!test:basic! and \lstinline!FOG:build:basic!, without recompiling \HEJ
again. Additionally artifacts can be downloaded from the GitLab web page, which
could be handy for debugging.
We also trigger some jobs \lstinline!only! on specific events. For example we
only push the code to
\href{https://phab.hepforge.org/source/hej/repository/v2.0/}{HepForge} on
release branches (e.g. v2.0). Also we only execute the \textit{long} tests for
merge requests, on a release or the \lstinline!master! branch, or when triggered
manually from the GitLab web page.
The actual commands are given in the \lstinline!before_script!,
\lstinline!script! and \lstinline!after_script!
\footnote{\lstinline!after_script! is always executed} sections, and are
standard Linux shell commands (dependent on the docker image). Any failed
command, i.e. returning not zero, stops the job and making the pipeline fail
entirely. Most tests are just running \lstinline!make test! or are based on it.
Thus, to emphasise it again, write tests for your code in \lstinline!cmake!. The
CI is only intended to make automated testing in different environments easier.
\bibliographystyle{JHEP}
\bibliography{biblio}
\end{document}
diff --git a/doc/sphinx/HEJ.rst b/doc/sphinx/HEJ.rst
index bf8e1fe..c70cfe7 100644
--- a/doc/sphinx/HEJ.rst
+++ b/doc/sphinx/HEJ.rst
@@ -1,324 +1,324 @@
.. _`Running HEJ 2`:
Running HEJ 2
=============
Quick start
-----------
In order to run HEJ 2, you need a configuration file and a file
containing fixed-order events. A sample configuration is given by the
:file:`config.yml` file distributed together with HEJ 2. Events in the
Les Houches Event File format can be generated with standard Monte Carlo
generators like `MadGraph5_aMC@NLO <https://launchpad.net/mg5amcnlo>`_
or `Sherpa <https://sherpa.hepforge.org/trac/wiki>`_. If HEJ 2 was
compiled with `HDF5 <https://www.hdfgroup.org/>`_ support, it can also
read event files in the format suggested in
`arXiv:1905.05120 <https://arxiv.org/abs/1905.05120>`_.
HEJ 2 assumes that the cross section is given by the sum of the event
weights. Depending on the fixed-order generator it may be necessary to
adjust the weights in the Les Houches Event File accordingly.
The processes supported by HEJ 2 are
- Pure multijet production
- Production of a Higgs boson with jets
- Production of a W boson with jets
..
- *TODO* Production of a Z boson or photon with jets
where at least two jets are required in each case. For the time being,
only leading-order events are supported.
After generating an event file :file:`events.lhe` adjust the parameters
under the `fixed order jets`_ setting in :file:`config.yml` to the
settings in the fixed-order generation. Resummation can then be added by
running::
HEJ config.yml events.lhe
Using the default settings, this will produce an output event file
:file:`HEJ.lhe` with events including high-energy resummation.
When using the `Docker image <https://hub.docker.com/r/hejdock/hej>`_,
HEJ can be run with
.. code-block:: bash
docker run -v $PWD:$PWD -w $PWD hejdock/hej HEJ config.yml events.lhe
.. _`HEJ 2 settings`:
Settings
--------
HEJ 2 configuration files follow the `YAML <http://yaml.org/>`_
format. The following configuration parameters are supported:
.. _`trials`:
**trials**
High-energy resummation is performed by generating a number of
resummation phase space configurations corresponding to an input
fixed-order event. This parameter specifies how many such
configurations HEJ 2 should try to generate for each input
event. Typical values vary between 10 and 100.
.. _`min extparton pt`:
**min extparton pt**
Specifies the minimum transverse momentum in GeV of the most forward
and the most backward parton. This setting is needed to regulate an
otherwise uncancelled divergence. Its value should be slightly below
the minimum transverse momentum of jets specified by `resummation
jets: min pt`_. See also the `max ext soft pt fraction`_ setting.
.. _`max ext soft pt fraction`:
**max ext soft pt fraction**
Specifies the maximum fraction that soft radiation can contribute to
the transverse momentum of each the most forward and the most backward
jet. Values between around 0.05 and 0.1 are recommended. See also the
`min extparton pt`_ setting.
.. _`fixed order jets`:
**fixed order jets**
This tag collects a number of settings specifying the jet definition
in the event input. The settings should correspond to the ones used in
the fixed-order Monte Carlo that generated the input events.
.. _`fixed order jets: min pt`:
**min pt**
Minimum transverse momentum in GeV of fixed-order jets.
.. _`fixed order 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.
.. _`fixed order jets: R`:
**R**
The R parameter used in the jet algorithm, roughly corresponding
to the jet radius in the plane spanned by the rapidity and the
azimuthal angle.
.. _`resummation jets`:
**resummation jets**
This tag collects a number of settings specifying the jet definition
in the observed, i.e. resummed events. These settings are optional, by
default the same values as for the `fixed order jets`_ are assumed.
.. _`resummation jets: min pt`:
**min pt**
Minimum transverse momentum in GeV of resummation jets. This
should be between 25% and 50% larger than the minimum transverse
momentum of fixed order jets set by `fixed order jets: min pt`_.
.. _`resummation jets: algorithm`:
**algorithm**
The algorithm used to define jets. The HEJ 2 approach to
resummation relies on properties of :code:`antikt` jets, so this
value is strongly recommended. For a list of possible other
values, see the `fixed order jets: algorithm`_ setting.
.. _`resummation jets: R`:
**R**
The R parameter used in the jet algorithm.
.. _`FKL`:
**FKL**
Specifies how to treat events respecting FKL rapidity ordering. These
configurations are dominant in the high-energy limit. The possible
values are :code:`reweight` to enable resummation, :code:`keep` to
keep the events as they are up to a possible change of
renormalisation and factorisation scale, and :code:`discard` to
discard these events.
.. _`unordered`:
**unordered**
Specifies how to treat events with one emission that does not respect FKL
ordering, e.g. :code:`u d => g u d`. In the high-energy limit, such
configurations are logarithmically suppressed compared to FKL configurations.
The possible values are the same as for the `FKL`_ setting, but
:code:`reweight` is currently only supported for Higgs or W bosons plus jets
production.
-.. _`ex_qqx`:
+.. _`extremal qqx`:
**extremal qqx**
Specifies how to treat events with a quark-antiquark pair as extremal partons
in rapidity, e.g. :code:`g d => u u_bar d`. In the high-energy limit, such
configurations are logarithmically suppressed compared to FKL configurations.
The possible values are the same as for the `FKL`_ setting, but
:code:`reweight` is currently only supported for W boson plus jets
production.
-.. _`min_qqx`:
+.. _`central qqx`:
**central qqx**
Specifies how to treat events with a quark-antiquark pair central in
rapidity, e.g. :code:`g g => g u u_bar g`. In the high-energy limit, such
configurations are logarithmically suppressed compared to FKL configurations.
The possible values are the same as for the `FKL`_ setting, but
:code:`reweight` is currently only supported for W boson plus jets
production.
-.. _`non-HEJ`:
+.. _`non-resummable`:
-**non-HEJ**
+**non-resummable**
Specifies how to treat events where no resummation is possible. The
allowed values are :code:`keep` to keep the events as they are up to
a possible change of renormalisation and factorisation scale and
:code:`discard` to discard these events.
.. _`scales`:
**scales**
Specifies the renormalisation and factorisation scales for the output
events. This can either be a single entry or a list :code:`[scale1,
scale2, ...]`. For the case of a list the first entry defines the
central scale. Possible values are fixed numbers to set the scale in
GeV or the following:
- :code:`H_T`: The sum of the scalar transverse momenta of all
final-state particles
- :code:`max jet pperp`: The maximum transverse momentum of all jets
- :code:`jet invariant mass`: Sum of the invariant masses of all jets
- :code:`m_j1j2`: Invariant mass between the two hardest jets.
Scales can be multiplied or divided by overall factors, e.g. :code:`H_T/2`.
It is also possible to import scales from an external library, see
:ref:`Custom scales`
.. _`scale factors`:
**scale factors**
A list of numeric factors by which each of the `scales`_ should be
multiplied. Renormalisation and factorisation scales are varied
independently. For example, a list with entries :code:`[0.5, 2]`
would give the four scale choices (0.5μ\ :sub:`r`, 0.5μ\ :sub:`f`);
(0.5μ\ :sub:`r`, 2μ\ :sub:`f`); (2μ\ :sub:`r`, 0.5μ\ :sub:`f`); (2μ\
:sub:`r`, 2μ\ :sub:`f`) in this order. The ordering corresponds to
the order of the final event weights.
.. _`max scale ratio`:
**max scale ratio**
Specifies the maximum factor by which renormalisation and
factorisation scales may difer. For a value of :code:`2` and the
example given for the `scale factors`_ the scale choices
(0.5μ\ :sub:`r`, 2μ\ :sub:`f`) and (2μ\ :sub:`r`, 0.5μ\ :sub:`f`)
will be discarded.
.. _`log correction`:
**log correction**
Whether to include corrections due to the evolution of the strong
coupling constant in the virtual corrections. Allowed values are
:code:`true` and :code:`false`.
.. _`event output`:
**event output**
Specifies the name of a single event output file or a list of such
files. The file format is either specified explicitly or derived from
the suffix. For example, :code:`events.lhe` or, equivalently
:code:`Les Houches: events.lhe` generates an output event file
:code:`events.lhe` in the Les Houches format. The supported formats
are
- :code:`file.lhe` or :code:`Les Houches: file`: The Les Houches
event file format.
- :code:`file.hepmc` or :code:`HepMC: file`: The HepMC format.
.. _`random generator`:
**random generator**
Sets parameters for random number generation.
.. _`random generator: name`:
**name**
Which random number generator to use. Currently, :code:`mixmax`
and :code:`ranlux64` are supported. Mixmax is recommended. See
the `CLHEP documentation
<http://proj-clhep.web.cern.ch/proj-clhep/index.html#docu>`_ for
details on the generators.
.. _`random generator: seed`:
**seed**
The seed for random generation. This should be a single number for
:code:`mixmax` and the name of a state file for :code:`ranlux64`.
.. _`analysis`:
**analysis**
Name and Setting for the event analyses; either a custom
analysis plugin or Rivet. For the first the :code:`plugin` sub-entry
should be set to the analysis file path. All further entries are passed on
to the analysis. To use Rivet a list of Rivet analyses have to be
given in :code:`rivet` and prefix for the yoda file has to be set
through :code:`output`. See :ref:`Writing custom analyses` for details.
.. _`Higgs coupling`:
**Higgs coupling**
This collects a number of settings concerning the effective coupling
of the Higgs boson to gluons. This is only relevant for the
production process of a Higgs boson with jets and only supported if
HEJ 2 was compiled with `QCDLoop
<https://github.com/scarrazza/qcdloop>`_ support.
.. _`Higgs coupling: use impact factors`:
**use impact factors**
Whether to use impact factors for the coupling to the most forward
and most backward partons. Impact factors imply the infinite
top-quark mass limit.
.. _`Higgs coupling: mt`:
**mt**
The value of the top-quark mass in GeV. If this is not specified,
the limit of an infinite mass is taken.
.. _`Higgs coupling: include bottom`:
**include bottom**
Whether to include the Higgs coupling to bottom quarks.
.. _`Higgs coupling: mb`:
**mb**
The value of the bottom-quark mass in GeV. Only used for the Higgs
coupling, external bottom-quarks are always assumed to be massless.
Advanced Settings
~~~~~~~~~~~~~~~~~
All of the following settings are optional. Please **do not set** any of the
following options, unless you know exactly what you are doing. The default
behaviour gives the most reliable results for a wide range of observables.
.. _`regulator parameter`:
**regulator parameter**
Slicing parameter to regularise the subtraction term, called :math:`\lambda`
in `arxiv:1706.01002 <https://arxiv.org/abs/1706.01002>`_. Default is 0.2
diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index c0a93a3..4bc711a 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,298 +1,298 @@
/** \file
* \brief Declares the Event class and helpers
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <memory>
#include <string>
#include <unordered_map>
#include <vector>
#include <boost/iterator/filter_iterator.hpp>
#include "HEJ/event_types.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/Particle.hh"
#include "HEJ/RNG.hh"
#include "fastjet/ClusterSequence.hh"
namespace LHEF{
class HEPEUP;
class HEPRUP;
}
namespace fastjet{
class JetDefinition;
}
namespace HEJ{
struct UnclusteredEvent;
/** @brief An event with clustered jets
*
* This is the main HEJ 2 event class.
* It contains kinematic information including jet clustering,
* parameter (e.g. scale) settings and the event weight.
*/
class Event{
public:
class EventData;
using ConstPartonIterator = boost::filter_iterator<
bool (*)(Particle const &),
std::vector<Particle>::const_iterator
>;
//! No default Constructor
Event() = delete;
//! Event Constructor adding jet clustering to an unclustered event
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
[[deprecated("UnclusteredEvent will be replaced by EventData")]]
Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
);
//! Incoming particles
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Outgoing particles
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
//! Iterator to the first outgoing parton
ConstPartonIterator begin_partons() const;
//! Iterator to the first outgoing parton
ConstPartonIterator cbegin_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator end_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator cend_partons() const;
//! Particle decays
/**
* The key in the returned map corresponds to the index in the
* vector returned by outgoing()
*/
std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
return decays_;
}
//! The jets formed by the outgoing partons, sorted in rapidity
std::vector<fastjet::PseudoJet> const & jets() const{
return jets_;
}
//! All chosen parameter, i.e. scale choices (const version)
Parameters<EventParameters> const & parameters() const{
return parameters_;
}
//! All chosen parameter, i.e. scale choices
Parameters<EventParameters> & parameters(){
return parameters_;
}
//! Central parameter choice (const version)
EventParameters const & central() const{
return parameters_.central;
}
//! Central parameter choice
EventParameters & central(){
return parameters_.central;
}
//! Parameter (scale) variations (const version)
std::vector<EventParameters> const & variations() const{
return parameters_.variations;
}
//! Parameter (scale) variations
std::vector<EventParameters> & variations(){
return parameters_.variations;
}
//! Parameter (scale) variation (const version)
/**
* @param i Index of the requested variation
*/
EventParameters const & variations(size_t i) const{
return parameters_.variations[i];
}
//! Parameter (scale) variation
/**
* @param i Index of the requested variation
*/
EventParameters & variations(size_t i){
return parameters_.variations[i];
}
//! Indices of the jets the outgoing partons belong to
/**
* @param jets Jets to be tested
* @returns A vector containing, for each outgoing parton,
* the index in the vector of jets the considered parton
* belongs to. If the parton is not inside any of the
* passed jets, the corresponding index is set to -1.
*/
std::vector<int> particle_jet_indices(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cs_.particle_jet_indices(jets);
}
//! Jet definition used for clustering
fastjet::JetDefinition const & jet_def() const{
return cs_.jet_def();
}
//! Minimum jet transverse momentum
double min_jet_pt() const{
return min_jet_pt_;
}
//! Event type
event_type::EventType type() const{
return type_;
}
//! Give colours to each particle
/**
- * @returns true if new colours are generated, i.e. same as is_HEJ()
+ * @returns true if new colours are generated, i.e. same as is_resummable()
* @details Colour ordering is done according to leading colour in the MRK
* limit, see \cite Andersen:2011zd. This only affects \ref
- * is_HEJ() "HEJ" configurations, all other \ref event_type
+ * is_resummable() "HEJ" configurations, all other \ref event_type
* "EventTypes" will be ignored.
* @note This overwrites all previously set colours.
*/
bool generate_colours(HEJ::RNG &);
private:
//! \internal
//! @brief Construct Event explicitly from input.
/** This is only intended to be called from EventData.
*
* \warning The input is taken _as is_, sorting and classification has to be
* done externally, i.e. by EventData
*/
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
);
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
std::unordered_map<size_t, std::vector<Particle>> decays_;
std::vector<fastjet::PseudoJet> jets_;
Parameters<EventParameters> parameters_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
event_type::EventType type_;
}; // end class Event
//! Class to store general Event setup, i.e. Phase space and weights
class Event::EventData{
public:
//! Default Constructor
EventData() = default;
//! Constructor from LesHouches event information
EventData(LHEF::HEPEUP const & hepeup);
//! Constructor with all values given
EventData(
std::array<Particle, 2> const & incoming_,
std::vector<Particle> const & outgoing_,
std::unordered_map<size_t, std::vector<Particle>> const & decays_,
Parameters<EventParameters> const & parameters_
):
incoming(incoming_), outgoing(outgoing_),
decays(decays_), parameters(parameters_)
{};
//! Move Constructor with all values given
EventData(
std::array<Particle, 2> && incoming_,
std::vector<Particle> && outgoing_,
std::unordered_map<size_t, std::vector<Particle>> && decays_,
Parameters<EventParameters> && parameters_
):
incoming(std::move(incoming_)), outgoing(std::move(outgoing_)),
decays(std::move(decays_)), parameters(std::move(parameters_))
{};
//! Generate an Event from the stored EventData.
/**
* @details Do jet clustering and classification.
* Use this to generate an Event.
*
* @note Calling this function destroys EventData
*
* @param jet_def Jet definition
* @param min_jet_pt minimal \f$p_T\f$ for each jet
*
* @returns Full clustered and classified event.
*/
Event cluster(
fastjet::JetDefinition const & jet_def, double const min_jet_pt);
//! Alias for cluster()
Event operator()(
fastjet::JetDefinition const & jet_def, double const min_jet_pt){
return cluster(jet_def, min_jet_pt);
};
//! Sort particles in rapidity
void sort();
//! Reconstruct intermediate particles from final-state leptons
/**
* Final-state leptons are created from virtual photons, W, or Z bosons.
* This function tries to reconstruct such intermediate bosons if they
* are not part of the event record.
*/
void reconstruct_intermediate();
std::array<Particle, 2> incoming;
std::vector<Particle> outgoing;
std::unordered_map<size_t, std::vector<Particle>> decays;
Parameters<EventParameters> parameters;
}; // end class EventData
//! Print Event
std::ostream& operator<<(std::ostream & os, Event const & ev);
//! Square of the partonic centre-of-mass energy \f$\hat{s}\f$
double shat(Event const & ev);
//! Convert an event to a LHEF::HEPEUP
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *);
// put deprecated warning at the end, so don't get the warning inside Event.hh,
// additionally doxygen can not identify [[deprecated]] correctly
struct [[deprecated("UnclusteredEvent will be replaced by EventData")]]
UnclusteredEvent;
//! An event before jet clustering
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
struct UnclusteredEvent{
//! Default Constructor
UnclusteredEvent() = default;
//! Constructor from LesHouches event information
UnclusteredEvent(LHEF::HEPEUP const & hepeup);
std::array<Particle, 2> incoming; /**< Incoming Particles */
std::vector<Particle> outgoing; /**< Outgoing Particles */
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<size_t, std::vector<Particle>> decays;
//! Central parameter (e.g. scale) choice
EventParameters central;
std::vector<EventParameters> variations; /**< For parameter variation */
};
}
diff --git a/include/HEJ/event_types.hh b/include/HEJ/event_types.hh
index 849b29a..e5be748 100644
--- a/include/HEJ/event_types.hh
+++ b/include/HEJ/event_types.hh
@@ -1,106 +1,106 @@
/** \file
* \brief Define different types of events.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include "HEJ/utility.hh"
namespace HEJ{
//! Namespace for event types
namespace event_type{
//! Possible event types
enum EventType: size_t{
- FixedOrder=0, /**< event configuration not covered by All Order resummation */
+ non_resummable=0, /**< event configuration not covered by All Order resummation */
bad_final_state=1, /**< event with an unsupported final state */
no_2_jets=2, /**< event with less than two jets */
FKL=4, /**< FKL-type event */
unordered_backward=8, /**< event with unordered backward emission */
unordered_forward=16, /**< event with unordered forward emission */
extremal_qqxb=32, /**< event with a backward extremal qqbar */
extremal_qqxf=64, /**< event with a forward extremal qqbar */
central_qqx=128, /**< event with a central qqbar */
unob = unordered_backward,
unof = unordered_forward,
qqxexb = extremal_qqxb,
qqxexf = extremal_qqxf,
qqxmid = central_qqx,
- first_type = FixedOrder,
+ first_type = non_resummable,
last_type = central_qqx
};
//! Event type names
/**
* For example, name(FKL) is the string "FKL"
*/
inline std::string name(EventType type) {
switch(type) {
case FKL:
return "FKL";
case unordered_backward:
return "unordered backward";
case unordered_forward:
return "unordered forward";
case extremal_qqxb:
return "extremal qqbar backward";
case extremal_qqxf:
return "extremal qqbar forward";
case central_qqx:
return "central qqbar";
- case FixedOrder:
- return "FixedOrder";
+ case non_resummable:
+ return "non-resummable";
case no_2_jets:
return "no 2 jets";
case bad_final_state:
return "bad final state";
default:
throw std::logic_error{"Unreachable"};
}
};
//! Returns True for a HEJ \ref event_type::EventType "EventType"
inline
- bool is_HEJ(EventType type) {
+ bool is_resummable(EventType type) {
switch(type) {
case FKL:
case unordered_backward:
case unordered_forward:
case extremal_qqxb:
case extremal_qqxf:
case central_qqx:
return true;
default:
return false;
}
}
//! Returns True for an unordered \ref event_type::EventType "EventType"
inline
bool is_uno(EventType type) {
return type == unordered_backward || type == unordered_forward;
}
//! Returns True for an extremal_qqx \ref event_type::EventType "EventType"
inline
bool is_ex_qqx(EventType type) {
return type == extremal_qqxb || type == extremal_qqxf;
}
//! Returns True for an central_qqx \ref event_type::EventType "EventType"
inline
bool is_mid_qqx(EventType type) {
return type == central_qqx;
}
//! Returns True for any qqx event \ref event_type::EventType "EventType"
inline
bool is_qqx(EventType type) {
return is_ex_qqx(type) || is_mid_qqx(type);
}
} // namespace event_type
} // namespace HEJ
diff --git a/src/Event.cc b/src/Event.cc
index 7166357..ed1ad8d 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,785 +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;
+ return EventType::non_resummable;
}
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;
+ return EventType::non_resummable;
}
begin_out+=2;
begin_idx+=2;
return not_tested | (backward?EventType::qqxexb:EventType::qqxexf);
}
}
}
- return EventType::FixedOrder;
+ return EventType::non_resummable;
}
//! 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;
+ if(std::distance(begin_out, end_out) < 0) return EventType::non_resummable;
// 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;
+ return EventType::non_resummable;
}
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;
+ if(begin_out->type != pid::gluon) return EventType::non_resummable;
++begin_idx;
}
return possible | EventType::qqxmid;
}
- return EventType::FixedOrder;
+ return EventType::non_resummable;
}
/**
* \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;
+ if( remaining_Wp != 0 || remaining_Wm != 0 ) return EventType::non_resummable;
// result has to be unique
- if( (final_type & (final_type-1)) != 0) return EventType::FixedOrder;
+ if( (final_type & (final_type-1)) != 0) return EventType::non_resummable;
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()))
+ if(!event_type::is_resummable(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;
}
}
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 5d35ee0..1396d29 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1755 +1,1755 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/currents.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/utility.hh"
namespace HEJ{
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
// use alpha_s(sqrt(q_j*lambda)), evolved to mur
return (
1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(
Event const & event
) const {
return tree(event)*virtual_corrections(event);
}
Weights MatrixElement::tree(
Event const & event
) const {
return tree_param(event)*tree_kin(event);
}
Weights MatrixElement::tree_param(
Event const & event
) const {
- if(! is_HEJ(event.type())) {
+ if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
Weights MatrixElement::virtual_corrections(
Event const & event
) const {
- if(! is_HEJ(event.type())) {
+ if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = virtual_corrections(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = virtual_corrections(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
bool wc = true;
bool wqq = false;
// With extremal qqx or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::FKL) {
if (in[0].type != partons[0].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::qqxexb) {
q -= partons[1].p;
++first_idx;
if (abs(partons[0].type) != abs(partons[1].type)){
q -= WBoson.p;
wc = false;
}
}
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(abs(backquark->type) != abs((backquark+1)->type)) {
wqq=true;
wc=false;
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqx) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
double MatrixElement::virtual_corrections(
Event const & event,
double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
const auto AWZH_boson = std::find_if(
begin(out), end(out),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson != end(out) && abs(AWZH_boson->type) == pid::Wp){
return virtual_corrections_W(event, mur, *AWZH_boson);
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
size_t first_idx = 0;
size_t last_idx = out.size() - 1;
// if there is a Higgs boson, extremal qqx or unordered gluon
// outside the extremal partons then it is not part of the FKL
// ladder and does not contribute to the virtual corrections
if((out.front().type == pid::Higgs)
|| event.type() == event_type::unob
|| event.type() == event_type::qqxexb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return 0;
last_idx = std::distance(begin(out), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqx) q -= out[last_idx+1].p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
} // namespace HEJ
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector qav,
CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1,
CLHEP::HepLorentzVector p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
}
//! Lipatov vertex
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,
CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector qav,
CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pa,
CLHEP::HepLorentzVector pb,
CLHEP::HepLorentzVector p1,
CLHEP::HepLorentzVector p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
return temp;
}
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
if (aptype==21&&bptype==21) {
return ME_gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_qg(pn,pb,p1,pa);
else
return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_qg(p1,pa,pn,pb);
else
return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_qQ(pn,pb,p1,pa);
else
return ME_qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return ME_qQbar(p1,pa,pn,pb);
else
return ME_qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// We know it cannot be gg incoming.
assert(!(aptype==21 && bptype==21));
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_W_qg(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_W_qg(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qbarg(p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (bptype>0) {
if (aptype>0)
return ME_W_qQ(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa);
}
else {
if (aptype>0)
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return ME_W_qQ(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0)
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb);
}
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*/
double ME_W_unob_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (bptype == 21 && aptype != 21) { // b gluon => W emission off a
if (aptype > 0)
return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else{
if (aptype>0)
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl);
else //qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else //qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
}
}
/** Matrix element squared for uno forward tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unof Tree-Level Current-Current Scattering
*/
double ME_W_unof_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (aptype==21 && bptype!=21) {//a gluon => W emission off b
if (bptype > 0)
return ME_Wuno_qg(pn, pb, p1, pa, pg, plbar, pl);
else
return ME_Wuno_qbarg(pn, pb, p1, pa, pg, plbar, pl);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return ME_Wuno_qQ(pn,pb,p1,pa,pg,plbar,pl);
else
return ME_Wuno_qQbar(pn,pb,p1,pa,pg,plbar,pl);
}
else{
if (aptype>0)
return ME_Wuno_qbarQ(pn,pb,p1,pa,pg,plbar,pl);
else
return ME_Wuno_qbarQbar(pn,pb,p1,pa,pg,plbar,pl);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return ME_W_unof_qQ(p1,pa,pn,pb,pg,plbar,pl);
// return ME_W_unof_qQ(pn,pb,p1,pa,pg,plbar,pl);
else //qqbar
return ME_W_unof_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return ME_W_unof_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else //qbarqbar
return ME_W_unof_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
}
}
/** \brief Matrix element squared for backward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*/
double ME_W_qqxb_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFbackward;
if (pqbar.rapidity() > pq.rapidity()){
swapQuarkAntiquark=true;
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pq.minus())+(pq.minus())/pa.minus())+1./3.)*3./4.;
}
else{
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pqbar.minus())+(pqbar.minus())/pa.minus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark){
return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn,pb)*CFbackward;}
else {
return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn,pb)*CFbackward;}
}
else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
if (wc!=1){ // W Emitted from backwards qqx
if (swapQuarkAntiquark){
return ME_WExqqx_qqbarQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
else{
return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
}
else { // W Must be emitted from forwards leg.
if(bptype > 0){
if (swapQuarkAntiquark){
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;}
else{
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;}
} else {
if (swapQuarkAntiquark){
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;}
else{
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, true)*CFbackward;}
}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxb");
}
}
/* \brief Matrix element squared for forward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param p1 Final state 1 Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxf Tree-Level Current-Current Scattering
*/
double ME_W_qqxf_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFforward;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pq.plus())+(pq.plus())/pb.plus())+1./3.)*3./4.;
}
else{
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pqbar.plus())+(pqbar.plus())/pb.plus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark){
return ME_WExqqx_qqbarg(pb, pqbar, plbar, pl, pq, p1,pa)*CFforward;}
else {
return ME_WExqqx_qbarqg(pb, pq, plbar, pl, pqbar, p1,pa)*CFforward;}
}
else if (bptype==21&&aptype!=21) {// b gluon => W emission off a or qqx
if (wc==1){ // W Emitted from forwards qqx
if (swapQuarkAntiquark){
return ME_WExqqx_qbarqQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
else {
return ME_WExqqx_qqbarQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
}
// W Must be emitted from backwards leg.
if (aptype > 0){
if (swapQuarkAntiquark){
return ME_W_Exqqx_QQq(pa,pb, p1, pqbar, pq, plbar, pl, false)*CFforward;}
else{
return ME_W_Exqqx_QQq(pa,pb, p1, pq, pqbar, plbar, pl, false)*CFforward;}
} else
{
if (swapQuarkAntiquark){
return ME_W_Exqqx_QQq(pa,pb, p1, pqbar, pq, plbar, pl, true)*CFforward;}
else{
return ME_W_Exqqx_QQq(pa,pb, p1, pq, pqbar, plbar, pl, true)*CFforward;}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxf");
}
}
/* \brief Matrix element squared for central qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqxpair
* @param nbelow Number of gluons emitted after central qqxpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqx
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_W_qqxmid_current(
int aptype, int bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<HLV> partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
bool swapQuarkAntiquark=false;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
}
double wt=1.;
if (aptype==21) wt*=K_g(partons.front(),pa)/HEJ::C_F;
if (bptype==21) wt*=K_g(partons.back(),pb)/HEJ::C_F;
if (aptype <=0 && bptype <=0){ // Both External AntiQuark
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl,plbar, partons,true,true,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true,true,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true,true,
swapQuarkAntiquark, nabove, nbelow, false);
}
} // end both antiquark
else if (aptype<=0){ // a is antiquark
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl,plbar, partons, false, true,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons,false,true,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, false, true,
swapQuarkAntiquark, nabove, nbelow, false);
}
} // end a is antiquark
else if (bptype<=0){ // b is antiquark
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl,plbar, partons, true, false,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true, false,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true, false,
swapQuarkAntiquark, nabove, nbelow, false);
}
} //end b is antiquark
else{ //Both Quark or gluon
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons, false, false,
swapQuarkAntiquark, nabove);}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, false, false,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, false, false,
swapQuarkAntiquark, nabove, nbelow, false);
}
}
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype==21) // gg initial state
return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return ME_H_qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered forward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param punof Unordered Particle Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered forward emission
*/
double ME_Higgs_current_unof(
int aptype, int bptype,
CLHEP::HepLorentzVector const & punof,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_H_unof_qg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unof_qbarg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_H_unof_qQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unof_qQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (aptype>0)
return ME_H_unof_qbarQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unof_qbarQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param punob Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*/
double ME_Higgs_current_unob(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & punob,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_unob_gQ(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unob_gQbar(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return ME_H_unob_qQ(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unob_qbarQ(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (bptype>0)
return ME_H_unob_qQbar(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return ME_H_unob_qbarQbar(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(HEJ::MatrixElementConfig const & config) {
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace anonymous
namespace HEJ{
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
double MatrixElement::tree_kin(
Event const & ev
) const {
- if(! is_HEJ(ev.type())) return 0.;
+ if(! is_resummable(ev.type())) return 0.;
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return tree_kin_jets(ev);
switch(AWZH_boson->type){
case pid::Higgs:
return tree_kin_Higgs(ev);
case pid::Wp:
case pid::Wm:
return tree_kin_W(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
// TODO: avoid reclustering
fastjet::ClusterSequence cs(to_PseudoJet(out_partons), ev.jet_def());
const auto jets = sorted_by_rapidity(cs.inclusive_jets(ev.min_jet_pt()));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission or qqx
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
assert(jets.size() >= 3);
++most_backward;
}
else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = cs.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(HEJ::is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
extremal_jet_idx:
no_extremal_jet_idx;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
double MatrixElement::tree_kin_jets(
Event const & ev
) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if(is_uno(ev.type())){
throw not_implemented("unordered emission not implemented for pure jets");
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn,
param_.regulator_lambda
);
}
namespace{
double tree_kin_W_FKL(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
double tree_kin_W_unob(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto pg = to_HepLorentzVector(partons[0]);
auto p1 = to_HepLorentzVector(partons[1]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - p1 -pg;
if (aptype!=partons[1].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_unob_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_unof(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 2]);
auto pg = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_unof_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxb(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
HLV pq,pqbar;
if(is_quark(partons[0])){
pq = to_HepLorentzVector(partons[0]);
pqbar = to_HepLorentzVector(partons[1]);
}
else{
pq = to_HepLorentzVector(partons[1]);
pqbar = to_HepLorentzVector(partons[0]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - pq - pqbar;
if (partons[1].type!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_qqxb_current(
aptype, bptype, pa, pb,
pq, pqbar, pn, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxf(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
HLV pq,pqbar;
if(is_quark(partons[partons.size() - 1])){
pq = to_HepLorentzVector(partons[partons.size() - 1]);
pqbar = to_HepLorentzVector(partons[partons.size() - 2]);
}
else{
pq = to_HepLorentzVector(partons[partons.size() - 2]);
pqbar = to_HepLorentzVector(partons[partons.size() - 1]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_qqxf_current(
aptype, bptype, pa, pb,
pq, pqbar, p1, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxmid(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wc, wqq;
if (backmidquark->type == -(backmidquark+1)->type){ // Central qqx does not emit
wqq=false;
if (aptype==partons[0].type) {
wc = true;
}
else{
wc = false;
q0-=pl+plbar;
}
}
else{
wqq = true;
wc = false;
qqxt-=pl+plbar;
}
auto begin_ladder = begin(partons) + 1;
auto end_ladder_1 = (backmidquark);
auto begin_ladder_2 = (backmidquark+2);
auto end_ladder = end(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark);
int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqxt, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const & decays(ev.decays());
HLV plbar, pl;
for (auto& x: decays) {
if (x.second.at(0).type < 0){
plbar = to_HepLorentzVector(x.second.at(0));
pl = to_HepLorentzVector(x.second.at(1));
}
else{
pl = to_HepLorentzVector(x.second.at(0));
plbar = to_HepLorentzVector(x.second.at(1));
}
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == unordered_backward){
return tree_kin_W_unob(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == unordered_forward){
return tree_kin_W_unof(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqxb(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqxf(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == central_qqx){
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
double MatrixElement::tree_kin_Higgs(
Event const & ev
) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev);
}
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
#ifdef HEJ_BUILD_WITH_QCDLOOP
// TODO: code duplication with currents.cc
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == ParticleID::gluon) return K_g(pout, pin);
return C_F;
}
#endif
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == ParticleID::gluon)?C_A:C_F;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
ParticleID type2,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
// gluon case
#ifdef HEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*ME_Houtside_gq(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(
Event const & ev
) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_last(
Event const & ev
) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(outgoing[outgoing.size()-2].type != pid::gluon) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
return MH2_forwardH(
pn, pb, incoming[0].type, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_between(
Event const & ev
) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor;
if(ev.type() == unob){
current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s) {
const auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](auto const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())) return 1.;
switch(AWZH_boson->type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return gw*gw*gw*gw/4.;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
double MatrixElement::tree_param(
Event const & ev,
double mur
) const{
- assert(is_HEJ(ev.type()));
+ assert(is_resummable(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*beta0*log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s)*res;
}
} // namespace HEJ
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index 282c5ef..59b5332 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,474 +1,474 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/YAMLreader.hh"
#include <algorithm>
#include <iostream>
#include <limits>
#include <map>
#include <string>
#include <unordered_map>
#include <vector>
#include <dlfcn.h>
#include "HEJ/ScaleFunction.hh"
#include "HEJ/event_types.hh"
#include "HEJ/output_formats.hh"
#include "HEJ/Constants.hh"
namespace HEJ{
class Event;
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 = {
"trials", "min extparton pt", "max ext soft pt fraction",
- "FKL", "unordered", "extremal qqx", "central qqx", "non-HEJ",
+ "FKL", "unordered", "extremal qqx", "central qqx", "non-resummable",
"scales", "scale factors", "max scale ratio", "import scales",
"log correction", "event output", "analysis", "regulator parameter"
};
// add subnodes to "supported" - the assigned value is irrelevant
for(auto && opt: opts) supported[opt] = "";
for(auto && jet_opt: {"min pt", "algorithm", "R"}){
supported["resummation jets"][jet_opt] = "";
supported["fixed order jets"][jet_opt] = "";
}
for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
supported["Higgs coupling"][opt] = "";
}
for(auto && opt: {"name", "seed"}){
supported["random generator"][opt] = "";
}
return supported;
}();
return supported;
}
fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){
using namespace fastjet;
static const std::map<std::string, fastjet::JetAlgorithm> known = {
{"kt", kt_algorithm},
{"cambridge", cambridge_algorithm},
{"antikt", antikt_algorithm},
{"cambridge for passive", cambridge_for_passive_algorithm},
{"plugin", plugin_algorithm}
};
const auto res = known.find(algo);
if(res == known.end()){
throw std::invalid_argument("Unknown jet algorithm " + algo);
}
return res->second;
}
EventTreatment to_EventTreatment(std::string const & name){
static const std::map<std::string, EventTreatment> known = {
{"reweight", EventTreatment::reweight},
{"keep", EventTreatment::keep},
{"discard", EventTreatment::discard}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown event treatment " + name);
}
return res->second;
}
} // namespace anonymous
namespace detail{
void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){
setting = to_JetAlgorithm(yaml.as<std::string>());
}
void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){
setting = to_EventTreatment(yaml.as<std::string>());
}
void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){
setting = to_ParticleID(yaml.as<std::string>());
}
} // namespace detail
JetParameters get_jet_parameters(
YAML::Node const & node,
std::string const & entry
){
assert(node);
JetParameters result;
fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm;
double R;
set_from_yaml_if_defined(jet_algo, node, entry, "algorithm");
set_from_yaml(R, node, entry, "R");
result.def = fastjet::JetDefinition{jet_algo, R};
set_from_yaml(result.min_pt, node, entry, "min pt");
return result;
}
RNGConfig to_RNGConfig(
YAML::Node const & node,
std::string const & entry
){
assert(node);
RNGConfig result;
set_from_yaml(result.name, node, entry, "name");
set_from_yaml_if_defined(result.seed, node, entry, "seed");
return result;
}
HiggsCouplingSettings get_Higgs_coupling(
YAML::Node const & node,
std::string const & entry
){
assert(node);
static constexpr double mt_max = 2e4;
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(node[entry]){
throw std::invalid_argument{
"Higgs coupling settings require building HEJ 2 "
"with QCDloop support"
};
}
#endif
HiggsCouplingSettings settings;
set_from_yaml_if_defined(settings.mt, node, entry, "mt");
set_from_yaml_if_defined(settings.mb, node, entry, "mb");
set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom");
set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors");
if(settings.use_impact_factors){
if(settings.mt != std::numeric_limits<double>::infinity()){
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
else{
// huge values of the top mass are numerically unstable
settings.mt = std::min(settings.mt, mt_max);
}
return settings;
}
FileFormat to_FileFormat(std::string const & name){
static const std::map<std::string, FileFormat> known = {
{"Les Houches", FileFormat::Les_Houches},
{"HepMC", FileFormat::HepMC}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown file format " + name);
}
return res->second;
}
std::string extract_suffix(std::string const & filename){
size_t separator = filename.rfind('.');
if(separator == filename.npos) return {};
return filename.substr(separator + 1);
}
FileFormat format_from_suffix(std::string const & filename){
const std::string suffix = extract_suffix(filename);
if(suffix == "lhe") return FileFormat::Les_Houches;
if(suffix == "hepmc") return FileFormat::HepMC;
throw std::invalid_argument{
"Can't determine format for output file " + filename
};
}
void assert_all_options_known(
YAML::Node const & conf, YAML::Node const & supported
){
if(!conf.IsMap()) return;
if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"};
for(auto const & entry: conf){
const auto name = entry.first.as<std::string>();
if(! supported[name]) throw unknown_option{name};
/* check sub-options, e.g. 'resummation jets: min pt'
* we don't check analysis sub-options
* those depend on the analysis being used and should be checked there
* similar for "import scales"
*/
if(name != "analysis" && name != "import scales"){
try{
assert_all_options_known(conf[name], supported[name]);
}
catch(unknown_option const & ex){
throw unknown_option{name + ": " + ex.what()};
}
catch(invalid_type const & ex){
throw invalid_type{name + ": " + ex.what()};
}
}
}
}
} // namespace HEJ
namespace YAML {
Node convert<HEJ::OutputFile>::encode(HEJ::OutputFile const & outfile) {
Node node;
node[to_string(outfile.format)] = outfile.name;
return node;
};
bool convert<HEJ::OutputFile>::decode(Node const & node, HEJ::OutputFile & out) {
switch(node.Type()){
case NodeType::Map: {
YAML::const_iterator it = node.begin();
out.format = HEJ::to_FileFormat(it->first.as<std::string>());
out.name = it->second.as<std::string>();
return true;
}
case NodeType::Scalar:
out.name = node.as<std::string>();
out.format = HEJ::format_from_suffix(out.name);
return true;
default:
return false;
}
}
} // namespace YAML
namespace HEJ{
namespace detail{
void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){
setting = yaml.as<OutputFile>();
}
}
namespace{
void update_fixed_order_jet_parameters(
JetParameters & fixed_order_jets, YAML::Node const & yaml
){
if(!yaml["fixed order jets"]) return;
set_from_yaml_if_defined(
fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt"
);
fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm();
set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm");
double R = fixed_order_jets.def.R();
set_from_yaml_if_defined(R, yaml, "fixed order jets", "R");
fixed_order_jets.def = fastjet::JetDefinition{algo, R};
}
// like std::stod, but throw if not the whole string can be converted
double to_double(std::string const & str){
std::size_t pos;
const double result = std::stod(str, &pos);
if(pos < str.size()){
throw std::invalid_argument(str + " is not a valid double value");
}
return result;
}
using EventScale = double (*)(Event const &);
void import_scale_functions(
std::string const & file,
std::vector<std::string> const & scale_names,
std::unordered_map<std::string, EventScale> & known
) {
auto handle = dlopen(file.c_str(), RTLD_NOW);
char * error = dlerror();
if(error != nullptr) throw std::runtime_error{error};
for(auto const & scale: scale_names) {
void * sym = dlsym(handle, scale.c_str());
error = dlerror();
if(error != nullptr) throw std::runtime_error{error};
known.emplace(scale, reinterpret_cast<EventScale>(sym));
}
}
auto get_scale_map(
YAML::Node const & yaml
) {
std::unordered_map<std::string, EventScale> scale_map;
scale_map.emplace("H_T", H_T);
scale_map.emplace("max jet pperp", max_jet_pt);
scale_map.emplace("jet invariant mass", jet_invariant_mass);
scale_map.emplace("m_j1j2", m_j1j2);
if(yaml["import scales"]) {
if(! yaml["import scales"].IsMap()) {
throw invalid_type{"Entry 'import scales' is not a map"};
}
for(auto const & import: yaml["import scales"]) {
const auto file = import.first.as<std::string>();
const auto scale_names =
import.second.IsSequence()
?import.second.as<std::vector<std::string>>()
:std::vector<std::string>{import.second.as<std::string>()};
import_scale_functions(file, scale_names, scale_map);
}
}
return scale_map;
}
// simple (as in non-composite) scale functions
/**
* An example for a simple scale function would be H_T,
* H_T/2 is then composite (take H_T and then divide by 2)
*/
ScaleFunction parse_simple_ScaleFunction(
std::string const & scale_fun,
std::unordered_map<std::string, EventScale> const & known
) {
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
const auto it = known.find(scale_fun);
if(it != end(known)) return {it->first, it->second};
try{
const double scale = to_double(scale_fun);
return {scale_fun, FixedScale{scale}};
} catch(std::invalid_argument const &){}
throw std::invalid_argument{"Unknown scale choice: " + scale_fun};
}
std::string trim_front(std::string const & str){
const auto new_begin = std::find_if(
begin(str), end(str), [](char c){ return ! std::isspace(c); }
);
return std::string(new_begin, end(str));
}
std::string trim_back(std::string str){
size_t pos = str.size() - 1;
// use guaranteed wrap-around behaviour to check whether we have
// traversed the whole string
for(; pos < str.size() && std::isspace(str[pos]); --pos) {}
str.resize(pos + 1); // note that pos + 1 can be 0
return str;
}
ScaleFunction parse_ScaleFunction(
std::string const & scale_fun,
std::unordered_map<std::string, EventScale> const & known
){
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
// parse from right to left => a/b/c gives (a/b)/c
const size_t delim = scale_fun.find_last_of("*/");
if(delim == scale_fun.npos){
return parse_simple_ScaleFunction(scale_fun, known);
}
const std::string first = trim_back(std::string{scale_fun, 0, delim});
const std::string second = trim_front(std::string{scale_fun, delim+1});
if(scale_fun[delim] == '/'){
return parse_ScaleFunction(first, known)
/ parse_ScaleFunction(second, known);
}
else{
assert(scale_fun[delim] == '*');
return parse_ScaleFunction(first, known)
* parse_ScaleFunction(second, known);
}
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
using namespace event_type;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{FKL, EventTreatment::reweight},
{unob, EventTreatment::keep},
{unof, EventTreatment::keep},
{qqxexb, EventTreatment::keep},
{qqxexf, EventTreatment::keep},
{qqxmid, EventTreatment::keep},
- {FixedOrder, EventTreatment::keep}
+ {non_resummable, EventTreatment::keep}
};
set_from_yaml(treat.at(FKL), yaml, "FKL");
set_from_yaml(treat.at(unob), yaml, "unordered");
treat.at(unof) = treat.at(unob);
set_from_yaml(treat.at(qqxexb), yaml, "extremal qqx");
set_from_yaml(treat.at(qqxexf), yaml, "extremal qqx");
set_from_yaml(treat.at(qqxmid), yaml, "central qqx");
- set_from_yaml(treat.at(FixedOrder), yaml, "non-HEJ");
- if(treat[FixedOrder] == EventTreatment::reweight){
+ set_from_yaml(treat.at(non_resummable), yaml, "non-resummable");
+ if(treat[non_resummable] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight Fixed Order events"};
}
return treat;
}
Config to_Config(YAML::Node const & yaml){
try{
assert_all_options_known(yaml, get_supported_options());
}
catch(unknown_option const & ex){
throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
}
Config config;
config.resummation_jets = get_jet_parameters(yaml, "resummation jets");
config.fixed_order_jets = config.resummation_jets;
update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
set_from_yaml(config.min_extparton_pt, yaml, "min extparton pt");
// Sets the standard value, then changes this if defined
config.regulator_lambda=CLAMBDA;
set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter");
config.max_ext_soft_pt_fraction = 0.1;
set_from_yaml_if_defined(
config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
);
set_from_yaml(config.trials, yaml, "trials");
set_from_yaml(config.log_correction, yaml, "log correction");
config.treat = get_event_treatment(yaml);
set_from_yaml_if_defined(config.output, yaml, "event output");
config.rng = to_RNGConfig(yaml, "random generator");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
config.scales = to_ScaleConfig(yaml);
config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
} // namespace anonymous
ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
ScaleConfig config;
auto scale_funs = get_scale_map(yaml);
std::vector<std::string> scales;
set_from_yaml(scales, yaml, "scales");
config.base.reserve(scales.size());
std::transform(
begin(scales), end(scales), std::back_inserter(config.base),
[scale_funs](auto const & entry){
return parse_ScaleFunction(entry, scale_funs);
}
);
set_from_yaml_if_defined(config.factors, yaml, "scale factors");
config.max_ratio = std::numeric_limits<double>::infinity();
set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio");
return config;
}
Config load_config(std::string const & config_file){
try{
return to_Config(YAML::LoadFile(config_file));
}
catch(...){
std::cerr << "Error reading " << config_file << ":\n ";
throw;
}
}
} // namespace HEJ
diff --git a/t/ME_data/config_mt.yml b/t/ME_data/config_mt.yml
index a88b482..3e81c97 100644
--- a/t/ME_data/config_mt.yml
+++ b/t/ME_data/config_mt.yml
@@ -1,30 +1,30 @@
trials: 1
min extparton pt: 30
resummation jets:
min pt: 30
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
FKL: reweight
unordered: reweight
extremal qqx: discard
central qqx: discard
-non-HEJ: discard
+non-resummable: discard
scales: 125
log correction: false
random generator:
name: mixmax
seed: 1
Higgs coupling:
use impact factors: false
mt: 174
include bottom: false
diff --git a/t/ME_data/config_mtinf.yml b/t/ME_data/config_mtinf.yml
index e74df97..e99541a 100644
--- a/t/ME_data/config_mtinf.yml
+++ b/t/ME_data/config_mtinf.yml
@@ -1,25 +1,25 @@
trials: 1
min extparton pt: 30
resummation jets:
min pt: 30
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
FKL: reweight
unordered: reweight
extremal qqx: discard
central qqx: discard
-non-HEJ: discard
+non-resummable: discard
scales: 125
log correction: false
random generator:
name: mixmax
seed: 1
diff --git a/t/ME_data/config_mtmb.yml b/t/ME_data/config_mtmb.yml
index 8cd06b8..26a0760 100644
--- a/t/ME_data/config_mtmb.yml
+++ b/t/ME_data/config_mtmb.yml
@@ -1,31 +1,31 @@
trials: 1
min extparton pt: 30
resummation jets:
min pt: 30
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
FKL: reweight
unordered: reweight
extremal qqx: discard
central qqx: discard
-non-HEJ: discard
+non-resummable: discard
scales: 125
log correction: false
random generator:
name: mixmax
seed: 1
Higgs coupling:
use impact factors: false
mt: 174
include bottom: true
mb: 4.7
diff --git a/t/ME_data/config_w_ME.yml b/t/ME_data/config_w_ME.yml
index 3a35832..b605ec1 100644
--- a/t/ME_data/config_w_ME.yml
+++ b/t/ME_data/config_w_ME.yml
@@ -1,25 +1,25 @@
trials: 1
min extparton pt: 30
resummation jets:
min pt: 30
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
FKL: reweight
unordered: reweight
extremal qqx: reweight
central qqx: reweight
-non-HEJ: discard
+non-resummable: discard
scales: 125
log correction: false
random generator:
name: mixmax
seed: 1
diff --git a/t/check_res.cc b/t/check_res.cc
index a8a2546..0d04ea9 100644
--- a/t/check_res.cc
+++ b/t/check_res.cc
@@ -1,165 +1,165 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <iostream>
#include <math.h>
#include "LHEF/LHEF.h"
#include "HEJ/Event.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/stream.hh"
#define ASSERT(x) if(!(x)) { \
std::cerr << "Assertion '" #x "' failed.\n"; \
return EXIT_FAILURE; \
}
namespace{
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
const fastjet::JetDefinition Born_jet_def{jet_def};
constexpr double Born_jetptmin = 30;
constexpr double extpartonptmin = 30;
constexpr double max_ext_soft_pt_fraction = 0.1;
constexpr double jetptmin = 35;
constexpr bool log_corr = false;
using EventTreatment = HEJ::EventTreatment;
using namespace HEJ::event_type;
HEJ::EventTreatMap treat{
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
- {FixedOrder, EventTreatment::discard},
+ {non_resummable, EventTreatment::discard},
{unof, EventTreatment::discard},
{unob, EventTreatment::discard},
{qqxexb, EventTreatment::discard},
{qqxexf, EventTreatment::discard},
{qqxmid, EventTreatment::discard},
{FKL, EventTreatment::reweight}
};
/// true if colour is allowed for particle
bool correct_colour(HEJ::Particle const & part){
if(HEJ::is_AWZH_boson(part) && !part.colour) return true;
if(!part.colour) return false;
int const colour = part.colour->first;
int const anti_colour = part.colour->second;
if(part.type == HEJ::ParticleID::gluon)
return colour != anti_colour && colour > 0 && anti_colour > 0;
if(HEJ::is_quark(part))
return anti_colour == 0 && colour > 0;
return colour == 0 && anti_colour > 0;
}
bool correct_colour(HEJ::Event const & ev){
- if(!HEJ::event_type::is_HEJ(ev.type()))
+ if(!HEJ::event_type::is_resummable(ev.type()))
return true;
for(auto const & part: ev.incoming()){
if(!correct_colour(part))
return false;
}
for(auto const & part: ev.outgoing()){
if(!correct_colour(part))
return false;
}
return true;
}
};
int main(int argn, char** argv) {
if(argn == 5 && std::string(argv[4]) == "unof"){
--argn;
treat[unof] = EventTreatment::reweight;
treat[unob] = EventTreatment::discard;
treat[FKL] = EventTreatment::discard;
}
if(argn == 5 && std::string(argv[4]) == "unob"){
--argn;
treat[unof] = EventTreatment::discard;
treat[unob] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "splitf"){
--argn;
treat[qqxexb] = EventTreatment::discard;
treat[qqxexf] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "splitb"){
--argn;
treat[qqxexb] = EventTreatment::reweight;
treat[qqxexf] = EventTreatment::discard;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "qqxmid"){
--argn;
treat[qqxmid] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
if(argn != 4){
std::cerr << "Usage: check_res eventfile xsection tolerance [uno]";
return EXIT_FAILURE;
}
const double xsec_ref = std::stod(argv[2]);
const double tolerance = std::stod(argv[3]);
HEJ::istream in{argv[1]};
LHEF::Reader reader{in};
HEJ::PhaseSpacePointConfig psp_conf;
psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin};
psp_conf.min_extparton_pt = extpartonptmin;
psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
HEJ::MatrixElementConfig ME_conf;
ME_conf.log_correction = log_corr;
ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{};
HEJ::EventReweighterConfig conf;
conf.psp_config = std::move(psp_conf);
conf.ME_config = std::move(ME_conf);
conf.jet_param = psp_conf.jet_param;
conf.treat = treat;
reader.readEvent();
const bool has_Higgs = std::find(
begin(reader.hepeup.IDUP),
end(reader.hepeup.IDUP),
25
) != end(reader.hepeup.IDUP);
const double mu = has_Higgs?125.:91.188;
HEJ::ScaleGenerator scale_gen{
{{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1.
};
HEJ::Mixmax ran{};
HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran};
double xsec = 0.;
double xsec_err = 0.;
do{
auto ev_data = HEJ::Event::EventData{reader.hepeup};
ev_data.reconstruct_intermediate();
HEJ::Event ev{
ev_data.cluster(
Born_jet_def, Born_jetptmin
)
};
auto resummed_events = hej.reweight(ev, 100);
for(auto const & ev: resummed_events) {
ASSERT(correct_colour(ev));
ASSERT(isfinite(ev.central().weight));
xsec += ev.central().weight;
xsec_err += ev.central().weight*ev.central().weight;
}
} while(reader.readEvent());
xsec_err = std::sqrt(xsec_err);
const double significance =
std::abs(xsec - xsec_ref) / std::sqrt( xsec_err*xsec_err + tolerance*tolerance );
std::cout << xsec_ref << " +/- " << tolerance << " ~ "
<< xsec << " +- " << xsec_err << " => " << significance << " sigma\n";
if(significance > 3.){
std::cerr << "Cross section is off by over 3 sigma!\n";
return EXIT_FAILURE;
}
}
diff --git a/t/jet_config.yml b/t/jet_config.yml
index 480bcfd..0a021f9 100644
--- a/t/jet_config.yml
+++ b/t/jet_config.yml
@@ -1,27 +1,27 @@
trials: 10
min extparton pt: 30
resummation jets:
min pt: 35
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
FKL: reweight
unordered: keep
extremal qqx: discard
central qqx: keep
-non-HEJ: keep
+non-resummable: keep
log correction: false
scales: 91.188
random generator:
name: ranlux64
event output:
- tst.lhe
diff --git a/t/jet_config_with_import.yml b/t/jet_config_with_import.yml
index b3d41c6..1d3e04f 100644
--- a/t/jet_config_with_import.yml
+++ b/t/jet_config_with_import.yml
@@ -1,30 +1,30 @@
trials: 10
min extparton pt: 30
resummation jets:
min pt: 35
algorithm: antikt
R: 0.4
fixed order jets:
min pt: 30
FKL: reweight
unordered: discard
extremal qqx: discard
central qqx: discard
-non-HEJ: discard
+non-resummable: discard
log correction: false
scales: softest_jet_pt
event output:
- tst.lhe
random generator:
name: ranlux64
import scales:
./libscales.so: softest_jet_pt
diff --git a/t/test_classify.cc b/t/test_classify.cc
index f1be154..2351e3c 100644
--- a/t/test_classify.cc
+++ b/t/test_classify.cc
@@ -1,835 +1,835 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <iostream>
#include <random>
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#define ASSERT(x) if(!(x)) { \
throw std::logic_error("Assertion '" #x "' failed."); \
}
namespace {
const fastjet::JetDefinition jet_def{fastjet::JetAlgorithm::antikt_algorithm, 0.4};
const double min_jet_pt{30.};
const std::array<std::string, 6> all_quarks{"-4","-1","1","2","3","4"};
const std::array<std::string, 7> all_partons{"g","-2","-1","1","2","3","4"};
const std::array<std::string, 3> all_bosons{"h", "Wp", "Wm"};
static std::mt19937_64 ran{0};
void shuffle_particles(HEJ::Event::EventData & ev) {
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
std::shuffle(begin(ev.outgoing), end(ev.outgoing), ran);
}
// if pos_boson = -1 (or not implemented) -> no boson
// njet==7 is special: has less jets, i.e. multiple parton in one jet,
// pos_boson < 0 to select process (see list for details)
HEJ::Event::EventData get_process(int const njet, int const pos_boson){
HEJ::Event::EventData ev;
if(njet == 2){
switch(pos_boson){
case 0:
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 198, 33, -170, 291}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-154, 68, 44, 174}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -44, -101, 88, 141}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -322, 322}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 284, 284}, {}};
return ev;
case 1:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -6, 82, -159, 179}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 195, -106, 74, 265}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-189, 24, 108, 219}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -320, 320}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 343, 343}, {}};
return ev;
case 2:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -80, -80, -140, 180}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -60, -32, 0, 68}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 140, 112, 177, 281}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -246, 246}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 283, 283}, {}};
return ev;
default:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -72, 24, 18, 78}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 72, -24, 74, 106}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -46, 46}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 138, 138}, {}};
return ev;
}
}
if(njet == 3){
switch(pos_boson){
case 0:
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 152, -117, -88, 245}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-146, 62, -11, 159}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 126, -72, 96, 174}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-132, 127, 144, 233}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -335, 335}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 476, 476}, {}};
return ev;
case 1:
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-191, 188, -128, 297}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 184, -172, -8, 252}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-192, -88, 54, 218}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -591, 591}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 433, 433}, {}};
return ev;
case 2:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -42, 18, -49, 67}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, -54, -28, 62}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 99, 32, -16, 163}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -45, 4, 72, 85}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -199, 199}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 178, 178}, {}};
return ev;
case 3:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -65, -32, -76, 105}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -22, 31, -34, 51}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, -67, -36, 77}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 99, 68, -4, 173}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -278, 278}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 128, 128}, {}};
return ev;
default:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -90, -135, 30, 165}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-108, 198, 76, 238}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, -63, 126, 243}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -207, 207}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 439, 439}, {}};
return ev;
}
}
if(njet == 4){
switch(pos_boson){
case 0:
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-200, -155, -64, 261}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, 194, 57, 283}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 1, 32, 8, 33}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-198, -143, 186, 307}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -515, 515}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 626, 626}, {}};
return ev;
case 1:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, 61, -162, 263}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 72, -76, 257}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-200, 135, 144, 281}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-198, -186, 171, 321}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 1, -82, 122, 147}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -535, 535}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 734, 734}, {}};
return ev;
case 2:
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-180, -27, -164, 245}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-108, 78, -36, 138}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 196, -189, 68, 307}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-107, 136, 76, 189}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 199, 2, 178, 267}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -512, 512}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 634, 634}, {}};
return ev;
case 3:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, -30, -84, 90}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -72, 22, -96, 122}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 0, -51, 85}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 64, 72, -81, 177}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -64, 84, 116}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -409, 409}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 181, 181}, {}};
return ev;
case 4:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -72, -49, -72, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, 0, -36, 60}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 54, -36, 66}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, -77, -56, 117}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 64, 72, -81, 177}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -407, 407}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 126, 126}, {}};
return ev;
default:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 248, -56, -122, 282}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 249, 30, -10, 251}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-249, -18, 26, 251}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-248, 44, 199, 321}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -506, 506}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 599, 599}, {}};
return ev;
}
}
if(njet == 6){
switch(pos_boson){
case 0:
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 349, 330, -94, 505}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-315, -300, 0, 435}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 347, 306, 18, 463}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-249, -342, 162, 453}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 345, 312, 284, 545}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-324, -126, 292, 454}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-153, -180, 304, 385}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -1137, 1137}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 2103, 2103}, {}};
return ev;
case 1:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 242, 241, -182, 387}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 243, 238, -190, 409}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-218, -215, -74, 315}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-224, -224, 112, 336}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 241, 182, 154, 339}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -53, -234, 126, 271}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-231, 12, 156, 279}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -1117, 1117}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 1219, 1219}, {}};
return ev;
case 2:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 151, 102, -42, 187}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -86, -46, -17, 99}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 152, 153, 0, 249}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -60, -135, 64, 161}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 150, 123, 110, 223}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-154, -49, 98, 189}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-153, -148, 144, 257}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -504, 504}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 861, 861}, {}};
return ev;
case 3:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, 197, -66, 287}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-198, -189, -54, 279}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-200, -64, 2, 210}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 158, 6, 283}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-199, -184, 172, 321}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 196, 168, 177, 313}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 4, -86, 92, 126}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -745, 745}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 1074, 1074}, {}};
return ev;
case 4:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 151, 102, -42, 187}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -86, -133, -14, 159}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-154, -104, -8, 186}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -60, 11, 0, 61}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 152, 153, 0, 249}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 150, 125, 90, 215}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-153, -154, 126, 251}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -578, 578}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 730, 730}, {}};
return ev;
case 5:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -90, -94, 131}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 82, -74, 111}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -80, -64, 105}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -25, -36, 65}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 68, 92, -18, 170}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -78, 54, 95}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -513, 513}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 265, 265}, {}};
return ev;
case 6:
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 198, 197, -66, 287}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 4, -84, -18, 86}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-198, -60, -36, 210}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 196, -78, -36, 214}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-200, 45, 0, 205}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-199, -178, 2, 267}, {}});
ev.outgoing.push_back({HEJ::ParticleID::higgs, { 199, 158, 6, 283}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -850, 850}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 702, 702}, {}};
return ev;
default:
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-350, -112, -280, 462}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 347, 266, -322, 543}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-349, -314, -38, 471}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 349, 348, 12, 493}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, {-342, -54, 23, 347}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 345, -134, 138, 395}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -1589, 1589}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 1122, 1122}, {}};
return ev;
}
}
if(njet == 7){
switch(pos_boson){
case -1: // jet idx: -1 0 1 2 3 4 5
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, -18, -42, 47}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, 26, -18, 35}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, -24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -54, -6, 59}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -44, 8, 45}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -96, 44, 116}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, 88, 133}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -249, 249}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 299, 299}, {}};
return ev;
case -2: // jet idx: 0 1 2 3 4 -1 -1
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -80, -64, 105}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -84, -12, 85}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, 88, 133}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -24, 62, 82}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, -18, 42, 47}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, 20, 60, 65}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -215, 215}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 415, 415}, {}};
return ev;
case -3: // jet idx: 0 0 1 2 3 4 5
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -86, -70, 111}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -52, -36, 65}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -44, 109}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -60, 5, 77}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, 8, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -80, 64, 105}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -361, 361}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 312, 312}, {}};
return ev;
case -4: // jet idx: 0 1 2 3 4 5 5
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -40, -56, 69}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -88, 133}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -84, -72, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, 8, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -58, 30, 67}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -96, 96, 144}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -395, 395}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 337, 337}, {}};
return ev;
case -5: // jet idx: 0 1 -1 -1 2 3 4
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -64, -80, 105}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, 20, 0, 25}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, -10, 2, 15}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -72, 54, 102}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -60, 48, 77}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -253, 253}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 285, 285}, {}};
return ev;
case -6: // jet idx: 0 0 0 1 2 2 3
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -60, -48, 77}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -52, -36, 65}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -96, -58, 122}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, -24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, 8, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -70, 14, 75}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -403, 403}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 243, 243}, {}};
return ev;
case -7: // jet idx: 0 1 1 2 2 3 4
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -46, -46, 69}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -90, -70, 115}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -78, -54, 95}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 88, -28, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, -16, 101}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -60, 5, 77}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -424, 424}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 239, 239}, {}};
return ev;
case -8: // jet idx: 0 1 2 2 2 3 4
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 23, -84, -84, 121}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -11, 92, -8, 93}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -15, -36, 0, 39}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -5, -62, 10, 63}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -48, -96, 19, 109}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { 68, 87, 24, 113}, {}});
ev.outgoing.push_back({HEJ::ParticleID::gluon, { -12, 99, 88, 133}, {}});
ev.incoming[0] = {HEJ::ParticleID::gluon, { 0, 0, -311, 311}, {}};
ev.incoming[1] = {HEJ::ParticleID::gluon, { 0, 0, 360, 360}, {}};
return ev;
}
}
throw HEJ::unknown_option{"unkown process"};
}
bool couple_quark(std::string const & boson, std::string & quark){
if(abs(HEJ::to_ParticleID(boson)) == HEJ::ParticleID::Wp){
auto qflav{ HEJ::to_ParticleID(quark) };
if(!HEJ::is_anyquark(qflav)) return false;
const int W_charge = HEJ::to_ParticleID(boson)>0?1:-1;
if(W_charge*qflav < 0 && !(abs(qflav)%2)) return false; // not anti-down
if(W_charge*qflav > 0 && (abs(qflav)%2)) return false; // not up
quark=std::to_string(qflav-W_charge);
}
return true;
}
// create event corresponding from given Configuration
// overwrite_boson to force a specific boson position, indepentent from input
// (useful for njet == 7)
HEJ::Event parse_configuration(
std::array<std::string,2> const & in, std::vector<std::string> const & out,
int const overwrite_boson = 0
){
auto boson = std::find_if(out.cbegin(), out.cend(),
[](std::string id){ return !HEJ::is_parton(HEJ::to_ParticleID(id)); });
int const pos_boson = (overwrite_boson!=0)?overwrite_boson:
((boson == out.cend())?-1:std::distance(out.cbegin(), boson));
size_t njets = out.size();
if(boson != out.cend()) --njets;
HEJ::Event::EventData ev{get_process(njets, pos_boson)};
ASSERT((pos_boson<0) || (ev.outgoing[pos_boson].type == HEJ::ParticleID::higgs));
for(size_t i=0; i<out.size(); ++i){
ev.outgoing[i].type = HEJ::to_ParticleID(out[i]);
}
for(size_t i=0; i<in.size(); ++i){
ev.incoming[i].type = HEJ::to_ParticleID(in[i]);
}
shuffle_particles(ev);
return ev.cluster(jet_def, min_jet_pt);
}
bool match_expectation( HEJ::event_type::EventType expected,
std::array<std::string,2> const & in, std::vector<std::string> const & out,
int const overwrite_boson = 0
){
HEJ::Event ev{parse_configuration(in,out,overwrite_boson)};
if(ev.type() != expected){
std::cerr << "Expected type " << HEJ::event_type::name(expected)
<< " but found " << HEJ::event_type::name(ev.type()) << "\n" << ev;
auto jet_idx{ ev.particle_jet_indices(ev.jets()) };
std::cout << "Particle Jet indices: ";
for(int const i: jet_idx)
std::cout << i << " ";
std::cout << std::endl;
return false;
}
return true;
}
bool check_fkl(){
using namespace HEJ;
for(std::string const & first: all_partons) // all quark flavours
for(std::string const & last: all_partons){
for(int njet=2; njet<=6; ++njet){ // all multiplicities
if(njet==5) continue;
std::array<std::string,2> base_in{first,last};
std::vector<std::string> base_out(njet, "g");
base_out.front() = first;
base_out.back() = last;
if(!match_expectation(event_type::FKL, base_in, base_out))
return false;
for(auto const & boson: all_bosons) // any boson
for(int pos=0; pos<=njet; ++pos){ // at any position
auto in{base_in};
auto out{base_out};
// change quark flavours for W
int couple_idx{ std::uniform_int_distribution<int>{0,1}(ran) };
if(!couple_quark(boson, couple_idx?out.back():out.front()))
continue;
out.insert(out.begin()+pos, boson);
if(!match_expectation(event_type::FKL, in, out))
return false;
}
}
for(int i=-1; i>-9;--i){ // jet setup doesn't matter for FKL
std::array<std::string,2> base_in{first,last};
std::vector<std::string> base_out(7, "g");
base_out.front() = first;
base_out.back() = last;
if(!match_expectation(event_type::FKL, base_in, base_out, i))
return false;
}
}
return true;
}
bool check_uno(){
using namespace HEJ;
auto const b{ event_type::unob };
auto const f{ event_type::unof };
for(std::string const & uno: all_quarks) // all quark flavours
for(std::string const & fkl: all_partons){
for(int njet=3; njet<=6; ++njet){ // all multiplicities >2
if(njet==5) continue;
for(int i=0; i<2; ++i){ // forward & backwards
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
const int uno_pos = i?1:(njet-2);
const int fkl_pos = i?(njet-1):0;
base_in[i?0:1] = uno;
base_in[i?1:0] = fkl;
base_out[uno_pos] = uno;
base_out[fkl_pos] = fkl;
auto expectation{ i?b:f };
if( !match_expectation(expectation, base_in, base_out) )
return false;
for(auto const & boson: all_bosons){ // any boson
// at any position (higgs only inside FKL chain)
int start = 0;
int end = njet;
if(to_ParticleID(boson) == pid::higgs){
start = i?(uno_pos+1):fkl_pos;
end = i?(fkl_pos+1):uno_pos;
}
for(int pos=start; pos<=end; ++pos){
auto in{base_in};
auto out{base_out};
// change quark flavours for W
bool couple_idx{ std::uniform_int_distribution<int>{0,1}(ran) };
if(!couple_quark(boson, couple_idx?out[fkl_pos]:out[uno_pos]))
continue;
out.insert(out.begin()+pos, boson);
if(!match_expectation(expectation, in, out))
return false;
}
}
}
}
// test allowed jet configurations
if(!(match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -1)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -2)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -3)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -4)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -5)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -5)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -6)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -7)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -7)
&& match_expectation(f,{fkl,uno},{fkl,"g","g","g","g",uno,"g"}, -8)
&& match_expectation(b,{uno,fkl},{"g",uno,"g","g","g","g",fkl}, -8)
))
return false;
}
return true;
}
bool check_extremal_qqx(){
using namespace HEJ;
auto const b{ event_type::qqxexb };
auto const f{ event_type::qqxexf };
for(std::string const & qqx: all_quarks) // all quark flavours
for(std::string const & fkl: all_partons){
std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) };
for(int njet=3; njet<=6; ++njet){ // all multiplicities >2
if(njet==5) continue;
for(int i=0; i<2; ++i){ // forward & backwards
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
const int qqx_pos = i?0:(njet-2);
const int fkl_pos = i?(njet-1):0;
base_in[i?0:1] = "g";
base_in[i?1:0] = fkl;
base_out[fkl_pos] = fkl;
base_out[qqx_pos] = qqx;
base_out[qqx_pos+1] = qqx2;
auto expectation{ i?b:f };
if( !match_expectation(expectation, base_in, base_out) )
return false;
for(auto const & boson: all_bosons){ // any boson
// at any position (higgs only inside FKL chain)
int start = 0;
int end = njet;
if(to_ParticleID(boson) == pid::higgs){
start = i?(qqx_pos+2):fkl_pos;
end = i?(fkl_pos+1):qqx_pos;
}
for(int pos=start; pos<=end; ++pos){
auto in{base_in};
auto out{base_out};
// change quark flavours for W
bool couple_idx{ std::uniform_int_distribution<int>{0,1}(ran) };
if(couple_idx || !couple_quark(boson, out[fkl_pos]) ){
// (randomly) try couple to FKL, else fall-back to qqx
if(!couple_quark(boson, out[qqx_pos]))
couple_quark(boson, out[qqx_pos+1]);
}
out.insert(out.begin()+pos, boson);
if(!match_expectation(expectation, in, out))
return false;
}
}
}
}
// test allowed jet configurations
if(!(match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -1)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -2)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -3)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -4)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -5)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -5)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -6)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -7)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -7)
&& match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -8)
&& match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -8)
))
return false;
}
return true;
}
bool check_central_qqx(){
using namespace HEJ;
auto const t{ event_type::qqxmid };
for(std::string const & qqx: all_quarks) // all quark flavours
for(std::string const & fkl1: all_partons)
for(std::string const & fkl2: all_partons){
std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) };
for(int njet=4; njet<=6; ++njet){ // all multiplicities >3
if(njet==5) continue;
for(int qqx_pos=1; qqx_pos<njet-2; ++qqx_pos){ // any qqx position
std::array<std::string,2> base_in;
std::vector<std::string> base_out(njet, "g");
base_in[0] = fkl1;
base_in[1] = fkl2;
base_out.front() = fkl1;
base_out.back() = fkl2;
base_out[qqx_pos] = qqx;
base_out[qqx_pos+1] = qqx2;
if( !match_expectation(t, base_in, base_out) )
return false;
for(auto const & boson: all_bosons) // any boson
for(int pos=0; pos<=njet; ++pos){ // at any position
if( to_ParticleID(boson) == pid::higgs
&& (pos==qqx_pos || pos==qqx_pos+1) )
continue;
auto in{base_in};
auto out{base_out};
// change quark flavours for W
int couple_idx{ std::uniform_int_distribution<int>{0,2}(ran) };
// (randomly) try couple to FKL, else fall-back to qqx
if( couple_idx == 0 && couple_quark(boson, out.front()) ){}
else if( couple_idx == 1 && couple_quark(boson, out.back()) ){}
else {
if(!couple_quark(boson, out[qqx_pos]))
couple_quark(boson, out[qqx_pos+1]);
}
out.insert(out.begin()+pos, boson);
if(!match_expectation(t, in, out))
return false;
}
}
}
// test allowed jet configurations (non exhaustive collection)
if(!(match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -1)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -2)
&& match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -3)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -3)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -4)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -4)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -5)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -6)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g",qqx,qqx2,"g",fkl2}, -6)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g",qqx,qqx2,"g","g",fkl2}, -7)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -7)
&& match_expectation(t,{fkl1,fkl2},{fkl1,qqx,qqx2,"g","g","g",fkl2}, -8)
&& match_expectation(t,{fkl1,fkl2},{fkl1,"g","g","g",qqx,qqx2,fkl2}, -8)
))
return false;
}
return true;
}
- // this checks a (non excessive) list of Non-HEJ states
- bool check_non_HEJ(){
- auto type{ HEJ::event_type::FixedOrder};
+ // this checks a (non excessive) list of non-resummable states
+ bool check_non_resummable(){
+ auto type{ HEJ::event_type::non_resummable};
return
// 2j - crossing lines
match_expectation(type, {"g","2"}, {"2","g"})
&& match_expectation(type, {"-1","g"}, {"g","-1"})
&& match_expectation(type, {"1","-1"}, {"-1","1"})
&& match_expectation(type, {"g","2"}, {"2","g","h"})
&& match_expectation(type, {"1","2"}, {"2","h","1"})
&& match_expectation(type, {"1","-1"}, {"h","-1","1"})
&& match_expectation(type, {"g","2"}, {"Wp","1","g"})
&& match_expectation(type, {"1","-1"}, {"-2","Wp","1"})
&& match_expectation(type, {"4","g"}, {"g","3","Wp"})
&& match_expectation(type, {"1","-2"}, {"-1","Wm","1"})
&& match_expectation(type, {"g","3"}, {"4","g","Wm"})
&& match_expectation(type, {"1","3"}, {"Wm","4","1"})
// 2j - qqx
&& match_expectation(type, {"g","g"}, {"1","-1"})
&& match_expectation(type, {"g","g"}, {"-2","2","h"})
&& match_expectation(type, {"g","g"}, {"-4","Wp","3"})
&& match_expectation(type, {"g","g"}, {"Wm","-1","2"})
// 3j - crossing lines
&& match_expectation(type, {"g","4"}, {"4","g","g"})
&& match_expectation(type, {"-1","g"}, {"g","g","-1"})
&& match_expectation(type, {"1","3"}, {"3","g","1"})
&& match_expectation(type, {"-2","2"}, {"2","g","-2","h"})
&& match_expectation(type, {"-3","g"}, {"g","g","Wp","-4"})
&& match_expectation(type, {"1","-2"}, {"Wm","-1","g","1"})
&& match_expectation(type, {"-1","g"}, {"1","-1","-1"})
// higgs inside uno
&& match_expectation(type, {"-1","g"}, {"g","h","-1","g"})
&& match_expectation(type, {"-1","1"}, {"g","h","-1","1"})
&& match_expectation(type, {"g","2"}, {"g","2","h","g"})
&& match_expectation(type, {"-1","1"}, {"-1","1","h","g"})
// higgs outside uno
&& match_expectation(type, {"-1","g"}, {"h","g","-1","g"})
&& match_expectation(type, {"-1","1"}, {"-1","1","g","h"})
// higgs inside qqx
&& match_expectation(type, {"g","g"}, {"-1","h","1","g","g"})
&& match_expectation(type, {"g","g"}, {"g","-1","h","1","g"})
&& match_expectation(type, {"g","g"}, {"g","g","2","h","-2"})
// higgs outside qqx
&& match_expectation(type, {"g","g"}, {"h","-1","1","g","g"})
&& match_expectation(type, {"g","g"}, {"g","g","2","-2","h"})
// 4j - two uno
&& match_expectation(type, {"-2","2"}, {"g","-2","2","g"})
&& match_expectation(type, {"1","3"}, {"g","1","h","3","g"})
&& match_expectation(type, {"1","2"}, {"g","1","3","Wp","g"})
&& match_expectation(type, {"1","-2"}, {"g","Wm","1","-1","g"})
// 4j - two gluon outside
&& match_expectation(type, {"g","4"}, {"g","4","g","g"})
&& match_expectation(type, {"1","3"}, {"1","3","h","g","g"})
&& match_expectation(type, {"1","2"}, {"1","3","g","Wp","g"})
&& match_expectation(type, {"1","-2"}, {"1","Wm","-1","g","g"})
&& match_expectation(type, {"-1","g"}, {"g","g","-1","g"})
&& match_expectation(type, {"1","3"}, {"g","g","1","3","h"})
&& match_expectation(type, {"1","2"}, {"g","g","1","Wp","3"})
&& match_expectation(type, {"1","-2"}, {"Wm","g","g","1","-1"})
// 4j - ggx+uno
&& match_expectation(type, {"g","4"}, {"1","-1","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-3","3"})
&& match_expectation(type, {"g","4"}, {"1","-1","h","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-3","3","h"})
&& match_expectation(type, {"g","4"}, {"Wp","1","-1","3","g"})
&& match_expectation(type, {"2","g"}, {"g","2","-4","Wp","3"})
&& match_expectation(type, {"g","4"}, {"2","Wm","-1","4","g"})
&& match_expectation(type, {"2","g"}, {"g","2","Wp","-3","4"})
// 3j - crossing+uno
&& match_expectation(type, {"1","4"}, {"g","4","1"})
&& match_expectation(type, {"1","4"}, {"4","1","g"})
&& match_expectation(type, {"1","4"}, {"g","h","4","1"})
&& match_expectation(type, {"-1","-3"},{"Wm","g","-4","-1"})
&& match_expectation(type, {"1","4"}, {"3","1","Wp","g"})
&& match_expectation(type, {"1","4"}, {"3","1","g","h"})
// 3j - crossing+qqx
&& match_expectation(type, {"1","g"}, {"-1","1","g","1"})
&& match_expectation(type, {"1","g"}, {"-1","1","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","1","-1"})
&& match_expectation(type, {"g","1"}, {"g","1","1","-1"})
&& match_expectation(type, {"1","g"}, {"2","-2","g","1"})
&& match_expectation(type, {"1","g"}, {"2","-2","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","-2","2"})
&& match_expectation(type, {"g","1"}, {"g","1","-2","2"})
&& match_expectation(type, {"1","g"}, {"-1","1","h","g","1"})
&& match_expectation(type, {"1","g"}, {"-1","h","1","1","g"})
&& match_expectation(type, {"g","1"}, {"1","g","1","h","-1"})
&& match_expectation(type, {"g","1"}, {"h","g","1","1","-1"})
&& match_expectation(type, {"1","g"}, {"2","-2","1","g","h"})
&& match_expectation(type, {"g","1"}, {"g","h","1","-2","2"})
&& match_expectation(type, {"1","g"}, {"Wp","3","-4","g","1"})
&& match_expectation(type, {"3","g"}, {"-2","Wm","1","3","g"})
&& match_expectation(type, {"g","1"}, {"1","g","Wm","-3","4"})
&& match_expectation(type, {"g","-3"}, {"g","-3","-1","Wp","2"})
// 4j- gluon in qqx
&& match_expectation(type, {"g","1"}, {"1","g","-1","1"})
&& match_expectation(type, {"1","g"}, {"1","-1","g","1"})
&& match_expectation(type, {"g","1"}, {"1","g","Wm","-2","1"})
&& match_expectation(type, {"2","g"}, {"2","-2","g","Wp","1"})
&& match_expectation(type, {"g","g"}, {"Wp","3","g","-4","g"})
&& match_expectation(type, {"1","g"}, {"1","h","-1","g","1"})
// 6j - two qqx
&& match_expectation(type, {"g","g"}, {"1","-1","g","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"1","-1","g","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","1","-1","-1","g"})
&& match_expectation(type, {"g","g"}, {"h","1","-1","g","g","1","-1"})
&& match_expectation(type, {"g","g"}, {"1","Wp","-2","g","1","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","Wp","-1","g","1","-2"})
&& match_expectation(type, {"g","g"}, {"g","1","-1","Wm","2","-1","g"})
&& match_expectation(type, {"g","g"}, {"g","1","2","-1","Wm","-1","g"})
// random stuff (can be non-physical)
&& match_expectation(type, {"g","g"}, {"1","-2","2","-1"}) // != 2 qqx
&& match_expectation(type, {"g","g"}, {"1","-2","2","g"}) // could be qqx
&& match_expectation(type, {"e+","e-"},{"1","-1"}) // bad initial state
&& match_expectation(type, {"1","e-"}, {"g","1","Wm"}) // bad initial state
&& match_expectation(type, {"h","g"}, {"g","g"}) // bad initial state
&& match_expectation(type, {"-1","g"}, {"-1","1","1"}) // bad qqx
&& match_expectation(type, {"-1","g"}, {"1","1","-1"}) // crossing in bad qqx
&& match_expectation(type, {"-1","g"}, {"-2","1","1","Wp"}) // bad qqx
&& match_expectation(type, {"1","2"}, {"1","-1","g","g","g","2"}) // bad qqx
&& match_expectation(type, {"1","2"}, {"1","-1","-2","g","g","2"}) // gluon in bad qqx
&& match_expectation(type, {"g","g"}, {"-1","2","g","g"}) // wrong back qqx
&& match_expectation(type, {"g","g"}, {"g","g","2","1"}) // wrong forward qqx
&& match_expectation(type, {"g","g"}, {"g","-2","1","g"}) // wrong central qqx
;
}
// Events failing the jet requirements, e.g. qqx inside one jet
//! @FIXME they are currently allowed
bool check_illegal_jets(){
- auto type{ HEJ::event_type::FixedOrder};
+ auto type{ HEJ::event_type::non_resummable};
return true
// uno backward not in jet
&& match_expectation(type, {"1","g"}, {"g","1","g","g","g","g","g"}, -1)
// & also legal uno on other side
&& match_expectation(type, {"1","1"}, {"g","1","g","g","g","1","g"}, -1)
// qqx backward not in jet
&& match_expectation(type, {"g","2"}, {"-1","1","g","g","g","g","2"}, -1)
// uno forward not in jet
&& match_expectation(type, {"3","3"}, {"3","g","g","g","g","3","g"}, -2)
// qqx backward not in jet
&& match_expectation(type, {"g","g"}, {"g","g","g","g","g","-2","2"}, -2)
// uno backward in same jet
&& match_expectation(type, {"1","g"}, {"g","1","g","g","g","g","g"}, -3)
// & also legal uno on other side
&& match_expectation(type, {"1","1"}, {"g","1","g","g","g","1","g"}, -3)
// qqx backward in same jet
&& match_expectation(type, {"g","2"}, {"-4","4","g","g","g","g","2"}, -3)
// uno forward in same jet
&& match_expectation(type, {"3","2"}, {"3","g","g","g","g","2","g"}, -4)
// qqx backward in same jet
&& match_expectation(type, {"g","g"}, {"g","g","g","g","g","-2","2"}, -4)
// central qqx not in jet
&& match_expectation(type, {"1","2"}, {"1","g","-1","1","g","g","2"}, -5)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","-1","1","g","g","g","2"}, -6)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","g","g","g","2","-2","2"}, -6)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","-1","1","g","g","g","2"}, -7)
// central qqx in same jet
&& match_expectation(type, {"1","2"}, {"1","g","g","3","-3","g","2"}, -7)
// central qqx in same jet
&& match_expectation(type, {"g","3"}, {"g","g","-2","2","g","g","3"}, -8)
// central qqx in same jet
&& match_expectation(type, {"g","-2"}, {"g","g","g","2","-2","g","-2"}, -8)
;
}
// Two boson states, that are currently not implemented
bool check_bad_FS(){
auto type{ HEJ::event_type::bad_final_state};
return
match_expectation(type, {"g","g"}, {"g","h","h","g"})
&& match_expectation(type, {"g","g"}, {"h","g","h","g"})
&& match_expectation(type, {"g","-1"}, {"g","h","Wp","-2"})
&& match_expectation(type, {"-3","-1"},{"-4","g","Wp","Wp","-2"})
&& match_expectation(type, {"-4","-1"},{"-3","Wp","g","Wm","-2"})
&& match_expectation(type, {"-4","-1"},{"g","-3","Wp","Wm","-2"})
;
}
}
int main() {
// tests for "no false negatives"
// i.e. all HEJ-configurations get classified correctly
if(!check_fkl()) return EXIT_FAILURE;
if(!check_uno()) return EXIT_FAILURE;
if(!check_extremal_qqx()) return EXIT_FAILURE;
if(!check_central_qqx()) return EXIT_FAILURE;
// test for "no false positive"
- // i.e. non-HEJ gives non-HEJ
- if(!check_non_HEJ()) return EXIT_FAILURE;
+ // i.e. non-resummable gives non-resummable
+ if(!check_non_resummable()) return EXIT_FAILURE;
if(!check_illegal_jets()) return EXIT_FAILURE;
if(!check_bad_FS()) return EXIT_FAILURE;
//! @TODO missing:
//! checks for partons sharing a jet
return EXIT_SUCCESS;
}
diff --git a/t/test_colours.cc b/t/test_colours.cc
index 60dc566..d47b225 100644
--- a/t/test_colours.cc
+++ b/t/test_colours.cc
@@ -1,281 +1,281 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <random>
#include <stdexcept>
#include <utility>
#include "HEJ/Event.hh"
#include "HEJ/RNG.hh"
#define ASSERT(x) if(!(x)) { \
throw std::logic_error("Assertion '" #x "' failed."); \
}
/// biased RNG to connect always to colour
class dum_rnd: public HEJ::DefaultRNG {
public:
dum_rnd() = default;
double flat() override {
return 0.;
};
};
void shuffle_particles(HEJ::Event::EventData & ev) {
static std::mt19937_64 ran{0};
std::shuffle(begin(ev.incoming), end(ev.incoming), ran);
std::shuffle(begin(ev.outgoing), end(ev.outgoing), ran);
}
void dump_event(HEJ::Event const & ev){
for(auto const & in: ev.incoming()){
std::cerr << "in type=" << in.type
<< ", colour={" << (*in.colour).first
<< ", " << (*in.colour).second << "}\n";
}
for(auto const & out: ev.outgoing()){
std::cerr << "out type=" << out.type << ", colour={";
if(out.colour)
std::cerr << (*out.colour).first << ", " << (*out.colour).second;
else
std::cerr << "non, non";
std::cerr << "}\n";
}
}
/// true if colour is allowed for particle
bool correct_colour(HEJ::Particle const & part){
if(HEJ::is_AWZH_boson(part) && !part.colour) return true;
if(!part.colour) return false;
int const colour = part.colour->first;
int const anti_colour = part.colour->second;
if(part.type == HEJ::ParticleID::gluon)
return colour != anti_colour && colour > 0 && anti_colour > 0;
if(HEJ::is_quark(part))
return anti_colour == 0 && colour > 0;
return colour == 0 && anti_colour > 0;
}
bool correct_colour(HEJ::Event const & ev){
for(auto const & part: ev.incoming()){
if(!correct_colour(part))
return false;
}
for(auto const & part: ev.outgoing()){
if(!correct_colour(part))
return false;
}
return true;
}
bool match_expected(
HEJ::Event const & ev,
std::vector<HEJ::Colour> const & expected
){
ASSERT(ev.outgoing().size()+2==expected.size());
for(size_t i=0; i<ev.incoming().size(); ++i){
ASSERT(ev.incoming()[i].colour);
if( *ev.incoming()[i].colour != expected[i])
return false;
}
for(size_t i=2; i<ev.outgoing().size()+2; ++i){
if( ev.outgoing()[i-2].colour ){
if( *ev.outgoing()[i-2].colour != expected[i] )
return false;
} else if( expected[i].first != 0 || expected[i].second != 0)
return false;
}
return true;
}
void check_event(
HEJ::Event::EventData unc_ev, std::vector<HEJ::Colour> const & expected_colours
){
shuffle_particles(unc_ev); // make sure incoming order doesn't matter
HEJ::Event ev{unc_ev.cluster(
fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30.)
};
- ASSERT(HEJ::event_type::is_HEJ(ev.type()));
+ ASSERT(HEJ::event_type::is_resummable(ev.type()));
dum_rnd rng;
ASSERT(ev.generate_colours(rng));
if(!correct_colour(ev)){
std::cerr << "Found illegal colours for event\n";
dump_event(ev);
throw std::invalid_argument("Illegal colour set");
}
if(!match_expected(ev, expected_colours)){
std::cerr << "Colours didn't match expectation. Found\n";
dump_event(ev);
std::cerr << "but expected\n";
for(auto const & col: expected_colours){
std::cerr << "colour={" << col.first << ", " << col.second << "}\n";
}
throw std::logic_error("Colours did not match expectation");
}
}
int main() {
HEJ::Event::EventData ev;
std::vector<HEJ::Colour> expected_colours(7);
/// pure gluon
ev.incoming[0] = { HEJ::ParticleID::gluon, { 0, 0,-427, 427}, {}};
ev.incoming[1] = { HEJ::ParticleID::gluon, { 0, 0, 851, 851}, {}};
ev.outgoing.push_back({ HEJ::ParticleID::gluon, { 196, 124, -82, 246}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-167,-184, 16, 249}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::higgs, { 197, 180, 168, 339}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-190, -57, 126, 235}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, { -36, -63, 196, 209}, {}});
expected_colours[0] = {502, 501};
expected_colours[1] = {509, 502};
expected_colours[2] = {503, 501};
expected_colours[3] = {505, 503};
expected_colours[4] = { 0, 0};
expected_colours[5] = {507, 505};
expected_colours[6] = {509, 507};
check_event(ev, expected_colours);
/// last g to Qx (=> gQx -> g ... Qx)
ev.incoming[1].type = HEJ::ParticleID::d_bar;
ev.outgoing[4].type = HEJ::ParticleID::d_bar;
// => only end changes
expected_colours[1].first = 0;
expected_colours[6].first = 0;
check_event(ev, expected_colours);
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno forward (=> gQx -> g ... Qx g)
std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
// => uno quarks eats colour and gluon connects to anti-colour
new_expected[5] = {0, expected_colours[3].first};
new_expected[6] = {expected_colours[0].first, expected_colours[0].first+2};
new_expected[1].second += 2; // one more anti-colour in line
check_event(new_ev, new_expected);
}
/// swap Qx <-> Q (=> gQ -> g ... Q)
ev.incoming[1].type = HEJ::ParticleID::d;
ev.outgoing[4].type = HEJ::ParticleID::d;
// => swap: colour<->anti && inital<->final
std::swap(expected_colours[1], expected_colours[6]);
std::swap(expected_colours[1].first, expected_colours[1].second);
std::swap(expected_colours[6].first, expected_colours[6].second);
check_event(ev, expected_colours);
/// first g to qx (=> qxQ -> qx ... Q)
ev.incoming[0].type = HEJ::ParticleID::u_bar;
ev.outgoing[0].type = HEJ::ParticleID::u_bar;
expected_colours[0] = { 0, 501};
// => shift anti-colour index one up
expected_colours[1].first -= 2;
expected_colours[5] = expected_colours[3];
expected_colours[3] = expected_colours[2];
expected_colours[2] = { 0, 502};
check_event(ev, expected_colours);
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno backward (=> qxQ -> g qx ... Q)
std::swap(new_ev.outgoing[0].type, new_ev.outgoing[1].type);
// => uno gluon connects to quark colour
new_expected[3] = expected_colours[2];
new_expected[2] = {expected_colours[0].second+2, expected_colours[0].second};
check_event(new_ev, new_expected);
/// swap qx <-> q (=> qQ -> g q ... Q)
new_ev.incoming[0].type = HEJ::ParticleID::u;
new_ev.outgoing[1].type = HEJ::ParticleID::u;
// => swap: colour<->anti && inital<->final
std::swap(new_expected[0], new_expected[3]);
std::swap(new_expected[0].first, new_expected[0].second);
std::swap(new_expected[3].first, new_expected[3].second);
// => & connect first gluon with remaining anti-colour
new_expected[2] = {new_expected[0].first, new_expected[0].first+2};
// shift colour line one down
new_expected[1].first-=2;
new_expected[5].first-=2;
new_expected[5].second-=2;
// shift anti-colour line one up
new_expected[6].first+=2;
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno forward (=> qxQ -> qx ... Q g)
std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
// => uno gluon connects to remaining colour
new_expected[5] = expected_colours[6];
new_expected[6] = {expected_colours[3].first+2, expected_colours[3].first};
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// qqx backward (=> gQ -> qx q ... Q)
// => swap: incoming q <-> outgoing gluon
std::swap(new_ev.incoming[0].type, new_ev.outgoing[1].type);
new_ev.outgoing[1].type=static_cast<HEJ::ParticleID>(
-1*new_ev.outgoing[1].type);
// incoming q -> outgoing q (colour<->anti)
std::swap(new_expected[0], new_expected[3]);
std::swap(new_expected[3].first, new_expected[3].second);
new_expected[3].first+=2;
new_expected[0].first-=1; // skip one index
// couple first in to first out
new_expected[2].second=new_expected[0].second;
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// qqx forward (=> qx g -> qx ... Qx Q)
// => swap: incoming Q <-> outgoing gluon
std::swap(new_ev.incoming[1].type, new_ev.outgoing[3].type);
new_ev.outgoing[3].type=static_cast<HEJ::ParticleID>(
-1*new_ev.outgoing[3].type);
// incoming q -> outgoing q (colour<->anti)
std::swap(new_expected[1], new_expected[5]);
std::swap(new_expected[5].first, new_expected[5].second);
new_expected[5].second-=2;
new_expected[1].second-=1; // skip one index
// couple last in to last out
new_expected[6].first=new_expected[1].first;
check_event(new_ev, new_expected);
// move Higgs to position 1 (=> qx g -> qx h g Qx Q)
std::swap(new_ev.outgoing[1].type, new_ev.outgoing[2].type);
std::swap(new_expected[3], new_expected[4]); // trivial
check_event(new_ev, new_expected);
// central qqx (=> qx g -> qx h Q Qx g)
// => swap: Q <-> g
std::swap(new_ev.outgoing[2].type, new_ev.outgoing[4].type);
std::swap(new_expected[4], new_expected[6]);
// gluon was connected on left side, i.e. doesn't matter for QQx
// => couple Q to out qx
new_expected[4].first = new_expected[2].second;
// Qx next in line
new_expected[5].second = new_expected[4].first+2;
// incoming g shifted by one position in line
new_expected[1].first-=2;
new_expected[1].second+=2;
check_event(new_ev, new_expected);
}
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:25 PM (15 h, 58 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486564
Default Alt Text
(281 KB)

Event Timeline