+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{Processes with interference}
+\label{sec:interference}
+
+For some processes several contributions can produce the same final state and it is necessary to take into account the interference between them. As an example let us consider the FKL squared matrix element for $Z/\gamma^*$ emission with two incoming quarks. It reads~\cite{Andersen:2016vkp}
+\todo{This equation is copied from the paper}
+\begin{align}
+\label{eq:ME_Z}
+\begin{split}
+|\mathcal{M}^{HEJ-{\rm reg}}_{qQ\to Z q(n-2)gQ}|^2 &=\ g_s^2 \frac{C_F}{8N_c}\ ( g_s^2
+The first two terms correspond to a boson emission from either the top or bottom quark line in both the amplitude and the complex conjugate, while the last term is the interference between the top and bottom emissions. To deal with this structure the intermediate objects are kept as vectors. For example, \lstinline!MatrixElement::FKL_ladder_weight_mix! returns the ladder factors for top, bottom and mixed emission, the latter one making use of \lstinline!C2Lipatovots_Mix! which implements the contraction $V^\mu(q_{ai},q_{a(i+1)})V_\mu(q_{bi},q_{b(i+1)})$. \lstinline!MatrixElement::operator()! and \lstinline!MatrixElement::tree! then use these vectors to perform the sum in eq.~(\ref{eq:ME_Z}).
+
\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}_\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)}
Currently, some of the expressions for the current contractions are hard-coded. Work is ongoing to use the current generator described in section~\ref{sec:cur_gen} instead.
\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. Naive unweighting can be
performed by using the \lstinline!set_cut_to_maxwt! member function of the
\lstinline!Unweighter! on the events and then call the
\lstinline!unweight! member function. It can be enabled for the
resummation events as explained in the user documentation.
To ameliorate the problem of naive unweighting, we also implement
partial unweighting. That is, we perform unweighting only for events
with sufficiently small weights. When using the \lstinline!Unweighter!
member function \lstinline!set_cut_to_peakwt! 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; events with a larger weight are not
touched. The actual unweighting is again done in the
\lstinline!Unweighter::unweight! function.
To estimate the peak weight we employ the following heuristic
algorithm. For a calibration sample of $n$ events, create a histogram
with $b=\sqrt{n}$ equal-sized bins. The histogram ranges from $
\log(\min |w_i|)$ to $\log(|\max w_i|)$, where $w_i$ are the event weights. For
each event, add $|w_i|$ to the corresponding bin. We then prune the
histogram by setting all bins containing less than $c=n/b$
events to zero. This effectively removes statistical outliers. The
logarithm of the peak weight is then the centre of the highest bin in
the histogram. In principle, the number of bins $b$ and the pruning
parameter $c$ could be tuned further.
To illustrate the principle, here is a weight-weight histogram
filled with a sample of 100000 event weights before the pruning:
\begin{center}
\includegraphics[width=0.7\linewidth]{wtwt}
\end{center}
The peaks to the right are clearly outliers caused by single
events. After pruning we get the following histogram:
\begin{center}
\includegraphics[width=0.7\linewidth]{wtwt_cut}
\end{center}
The actual peak weight probably lies above the cut, and the algorithm
can certainly be improved. Still, the estimate we get from the pruned
histogram is already good enough to eliminate about $99\%$ of the
low-weight events.
\section{The current generator}
\label{sec:cur_gen}
The current generator in the \texttt{current\_generator} directory is
automatically invoked when building \HEJ 2. Its task is to compute and
generate C++ code for the current contractions listed in
section~\ref{sec:currents}. For each source file \texttt{<j\_j>.frm}
inside \texttt{current\_generator} it generates a corresponding header
\texttt{include/HEJ/currents/<j\_j>.hh} inside the build directory. The
header can be included with
\begin{lstlisting}[language=C++,caption={}]
#include "HEJ/currents/<j>.hh"
\end{lstlisting}
The naming scheme is
\begin{itemize}
\item \texttt{j1\_j2.frm} for the contraction of current \texttt{j1} with current \texttt{j2}.
\item \texttt{j1\_vx\_j2.frm} for the contraction of current \texttt{j1}
with the vertex \texttt{vx} and the current \texttt{j2}.
\end{itemize}
For instance, \texttt{juno\_qqbarW\_j.frm} would indicate the
contraction of an unordered current with a (central) quark-antiquark-W
emission vertex and a standard FKL current.
\subsection{Implementing new current contractions}
\label{sec:cur_gen_new}
To implement a new current contraction \lstinline!jcontr! create a new
file \texttt{jcontr.frm} inside the \texttt{current\_generator}
directory. \texttt{jcontr.frm} should contain \FORM code, see
\url{https://www.nikhef.nl/~form/} for more information on \FORM. Here
is a small example:
\begin{lstlisting}[caption={}]
* FORM comments are started by an asterisk * at the beginning of a line
* First include the relevant headers
#include- include/helspin.frm
#include- include/write.frm
* Define the symbols that appear.
* UPPERCASE symbols are reserved for internal use
vectors p1,...,p10;
indices mu1,...,mu10;
* Define local expressions of the form [NAME HELICITIES]
* for the current contractions for all relevant helicity configurations
#do HELICITY1={+,-}
#do HELICITY2={+,-}
* We use the Current function
* Current(h, p1, mu1, ..., muX, p2) =
* u(p1) \gamma_{mu1} ... \gamma_{muX} u(p2)
* where h=+1 or h=-1 is the spinor helicity.
* All momenta appearing as arguments have to be *lightlike*
local [jcontr `HELICITY1'`HELICITY2'] =
Current(`HELICITY1'1, p1, mu1, p2, mu2, p3)
*Current(`HELICITY2'1, p4, mu2, p2, mu1, p1);
#enddo
#enddo
.sort
* Main procedure that calculates the contraction
#call ContractCurrents
.sort
* Optimise expression
format O4;
* Format in a (mostly) c compatible way
format c;
* Write start of C++ header file
#call WriteHeader(`OUTPUT')
* Write a template function jcontr
* taking as template arguments two helicities
* and as arguments the momenta p1,...,p4
* returning the contractions [jcontr HELICITIES] defined above
to split up currents that are contracted with momenta.
\item Apply the Fierz transformation~\eqref{eq:Fierz} to eliminate
contractions between vector currents.
\item Write the arguments of the antisymmetric angle and scalar products in canonical order, see equations~\eqref{eq:angle_product} ,\eqref{eq:square_product}.
\end{enumerate}
The corresponding \lstinline!ContractCurrents! procedure is implemented in
\texttt{include/helspin.fm}.
\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