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.
Currently, most 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
+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>.frm}
inside \texttt{current\_generator} it generates a corresponding header
\texttt{include/HEJ/currents/<j>.hh} inside the build directory. The
header can be included with
\begin{lstlisting}[language=C++,caption={}]
#include "HEJ/currents/<j>.hh"
\end{lstlisting}
\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}
\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
+$\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