diff --git a/doc/developer_manual/.gitattributes b/doc/developer_manual/.gitattributes
new file mode 100644
index 0000000..591ce85
--- /dev/null
+++ b/doc/developer_manual/.gitattributes
@@ -0,0 +1 @@
+figures/* filter=lfs diff=lfs merge=lfs -text
diff --git a/doc/developer_manual/currents.tex b/doc/developer_manual/currents.tex
index 5756461..6a3c7a5 100644
--- a/doc/developer_manual/currents.tex
+++ b/doc/developer_manual/currents.tex
@@ -1,226 +1,580 @@
 \section{Currents}
 \label{sec:currents}
-In the following is a list of all the currents currently implemented in HEJ. All of them are defined inside \texttt{include/currents.hh}.
-
+The following section contains a list of all the currents implemented
+in \HEJ. Clean up of the code structure is ongoing. All $W$+Jet currents
+are located in \texttt{src/Wjets.cc}. All other currents are defined
+in \texttt{src/currents.cc}. There is a common header between these
+files located at \texttt{include/HEJ/currents.hh}.
 \subsection{Pure Jets}
 
 \subsubsection{Quark}
 \label{sec:current_quark}
 \begin{align}
 j_\mu(p_i,p_j)=\bar{u}(p_i)\gamma_\mu u(p_j)
 \end{align}
 The exact for depends on the helicity and direction (forward/backwards) for the quarks. Currently all different contractions of incoming and outgoing states are defined in \lstinline!joi!, \lstinline!jio! and \lstinline!joo!.
 
 \subsubsection{Gluon}
 In \HEJ the currents for gluons and quarks are the same, up to a colour factor $K_g/C_F$, where
 \begin{align}
  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}.
 \end{align}
 Thus we can just reuse the results from sec.~\ref{sec:current_quark}.
 
 \subsubsection{Single unordered gluon}
 Configuration $q(p_a) \to g(p_1) q(p_2) g^*(\tilde{q}_2)$~\cite{Andersen:2017kfc}
 \begin{align}
 \label{eq:juno}
 \begin{split}
 &j^{{\rm uno}\; \mu\ cd}(p_2,p_1,p_a) = i \varepsilon_{1\nu} \left(  T_{2i}^{c}T_{ia}^d\
 \left(U_1^{\mu\nu}-L^{\mu\nu} \right) + T_{2i}^{d}T_{ia}^c\ \left(U_2^{\mu\nu} +
 L^{\mu\nu} \right) \right). \\
 U_1^{\mu\nu} &= \frac1{s_{21}} \left( j_{21}^\nu j_{1a}^\mu + 2 p_2^\nu
 j_{2a}^\mu \right) \qquad  \qquad U_2^{\mu\nu} = \frac1{t_{a1}} \left( 2
 j_{2a}^\mu p_a^\nu - j_{21}^\mu  j_{1a}^\nu \right) \\
 L^{\mu\nu} &= \frac1{t_{a2}} \left(-2p_1^\mu j_{2a}^\nu+2p_1.j_{2a}
 g^{\mu\nu} + (\tilde{q}_1+\tilde{q}_2)^\nu j_{2a}^\mu + \frac{t_{b2}}{2} j_{2a}^\mu \left(
 \frac{p_2^\nu}{p_1.p_2} + \frac{p_b^\nu}{p_1.p_b} \right) \right) ,
 \end{split}
 \end{align}
 $j^{{\rm uno}\; \mu}$ is currently not calculated as a separate current, but always as needed for the ME (i.e. in \lstinline!jM2unoXXX!).
 
 \subsection{Higgs}
 Different rapidity orderings \todo{give name of functions}
 \begin{enumerate}
 \item $qQ\to HqQ/qHQ/qQH$ (any rapidity order, full LO ME) $\Rightarrow$ see~\ref{sec:V_H}
 \item $qg\to Hqg$ (Higgs outside quark) $\Rightarrow$ see~\ref{sec:V_H}
 \item $qg\to qHg$ (central Higgs) $\Rightarrow$ see~\ref{sec:V_H}
 \item $qg\to qgH$ (Higgs outside gluon) $\Rightarrow$ see~\ref{sec:jH_mt}
 \item $gg\to gHg$ (central Higgs) $\Rightarrow$ see~\ref{sec:V_H}
 \item $gg\to ggH$ (Higgs outside gluon) $\Rightarrow$ see~\ref{sec:jH_mt}
 \end{enumerate}
 
 \subsubsection{Higgs gluon vertex}
 \label{sec:V_H}
 The coupling of the Higgs boson to gluons via a virtual quark loop can be written as
 \begin{align}
 \label{eq:VH}
 V^{\mu\nu}_H(q_1, q_2) = \mathgraphics{build/figures/V_H.pdf} &=
 \frac{\alpha_s m^2}{\pi v}\big[
 g^{\mu\nu} T_1(q_1, q_2) - q_2^\mu q_1^\nu T_2(q_1, q_2)
 \big]\, \\
 &\xrightarrow{m \to \infty} \frac{\alpha_s}{3\pi
     v} \left(g^{\mu\nu} q_1\cdot q_2 - q_2^\mu q_1^\nu\right).
 \end{align}
 The outgoing momentum of the Higgs boson is $p_H = q_1 - q_2$.
 As a contraction with two currents this by implemented in \lstinline!cHdot! inside \texttt{src/currents.cc}.
 The form factors $T_1$ and $T_2$ are then given by~\cite{DelDuca:2003ba}
 \begin{align}
 \label{eq:T_1}
 T_1(q_1, q_2) ={}& -C_0(q_1, q_2)\*\left[2\*m^2+\frac{1}{2}\*\left(q_1^2+q_2^2-p_H^2\right)+\frac{2\*q_1^2\*q_2^2\*p_H^2}{\lambda}\right]\notag\\
 &-\left[B_0(q_2)-B_0(p_H)\right]\*\frac{q_2^2}{\lambda}\*\left(q_2^2-q_1^2-p_H^2\right)\notag\\
 &-\left[B_0(q_1)-B_0(p_H)\right]\*\frac{q_1^2}{\lambda}\*\left(q_1^2-q_2^2-p_H^2\right)-1\,,\displaybreak[0]\\
 \label{eq:T_2}
 T_2(q_1, q_2) ={}& C_0(q_1,
 q_2)\*\left[\frac{4\*m^2}{\lambda}\*\left(p_H^2-q_1^2-q_2^2\right)-1-\frac{4\*q_1^2\*q_2^2}{\lambda}
 -
 \frac{12\*q_1^2\*q_2^2\*p_H^2}{\lambda^2}\*\left(q_1^2+q_2^2-p_H^2\right)\right]\notag\\
 &-\left[B_0(q_2)-B_0(p_H)\right]\*\left[\frac{2\*q_2^2}{\lambda}+\frac{12\*q_1^2\*q_2^2}{\lambda^2}\*\left(q_2^2-q_1^2+p_H^2\right)\right]\notag\\
 &-\left[B_0(q_1)-B_0(p_H)\right]\*\left[\frac{2\*q_1^2}{\lambda}+\frac{12\*q_1^2\*q_2^2}{\lambda^2}\*\left(q_1^2-q_2^2+p_H^2\right)\right]\notag\\
 &-\frac{2}{\lambda}\*\left(q_1^2+q_2^2-p_H^2\right)\,,
 \end{align}
 where we have used the scalar bubble and triangle integrals
 \begin{align}
 \label{eq:B0}
 B_0\left(p\right) ={}& \int \frac{d^dl}{i\pi^{\frac{d}{2}}}
 \frac{1}{\left(l^2-m^2\right)\left((l+p)^2-m^2\right)}\,,\\
 \label{eq:C0}
 C_0\left(p,q\right) ={}& \int \frac{d^dl}{i\pi^{\frac{d}{2}}} \frac{1}{\left(l^2-m^2\right)\left((l+p)^2-m^2\right)\left((l+p-q)^2-m^2\right)}\,,
 \end{align}
 and the K\"{a}ll\'{e}n function
 \begin{equation}
 \label{eq:lambda}
 \lambda = q_1^4 + q_2^4 + p_H^4 - 2\*q_1^2\*q_2^2 - 2\*q_1^2\*p_H^2- 2\*q_2^2\*p_H^2\,.
 \end{equation}
 The Integrals as such are provided by \QCDloop{} (see wrapper functions \lstinline!B0DD! and \lstinline!C0DD! in \texttt{src/currents.cc}).
 In the code we are sticking to the convention of~\cite{DelDuca:2003ba}, thus instead of the $T_{1/2}$ we implement (in the functions \lstinline!A1! and \lstinline!A2!)
 \begin{align}
 \label{eq:A_1}
 A_1(q_1, q_2) ={}& \frac{i}{16\pi^2}\*T_2(-q_1, q_2)\,,\\
 \label{eq:A_2}
 A_2(q_1, q_2) ={}& -\frac{i}{16\pi^2}\*T_1(-q_1, q_2)\,.
 \end{align}
 
 \subsubsection{Peripheral Higgs emission - Finite quark mass}
 \label{sec:jH_mt}
 
 We describe the emission of a peripheral Higgs boson close to a
 scattering gluon with an effective current. In the following we consider
 a lightcone decomposition of the gluon momenta, i.e. $p^\pm = E \pm p_z$
 and $p_\perp = p_x + i p_y$. The incoming gluon momentum $p_a$ defines
 the $-$ direction, so that $p_a^+ = p_{a\perp} = 0$. The outgoing
 momenta are $p_1$ for the gluon and $p_H$ for the Higgs boson. We choose
 the following polarisation vectors:
 \begin{equation}
   \label{eq:pol_vectors}
   \epsilon_\mu^\pm(p_a) = \frac{j_\mu^\pm(p_1, p_a)}{\sqrt{2}
     \bar{u}^\pm(p_a)u^\mp(p_1)}\,, \quad \epsilon_\mu^{\pm,*}(p_1) = -\frac{j_\mu^\pm(p_1, p_a)}{\sqrt{2}
     \bar{u}^\mp(p_1)u^\pm(p_a)}\,.
 \end{equation}
 Following~\cite{DelDuca:2001fn}, we introduce effective polarisation
 vectors to describe the contraction with the Higgs-boson production
 vertex eq.~\eqref{eq:VH}:
 \begin{align}
   \label{eq:eps_H}
   \epsilon_{H,\mu}(p_a) = \frac{T_2(p_a, p_a-p_H)}{(p_a-p_H)^2}\big[p_a\cdot
   p_H\epsilon_\mu(p_a) - p_H\cdot\epsilon(p_a) p_{a,\mu}\big]\,,\\
   \epsilon_{H,\mu}^*(p_1) = -\frac{T_2(p_1+p_H, p_1)}{(p_1+p_H)^2}\big[p_1\cdot
   p_H\epsilon_\mu^*(p_1) - p_H\cdot\epsilon^*(p_1) p_{1,\mu}\big]\,,
 \end{align}
 We also employ the usual short-hand notation
 \begin{equation}
   \label{eq:spinor_helicity}
   \spa i.j = \bar{u}^-(p_i)u^+(p_j)\,,\qquad \spb i.j =
   \bar{u}^+(p_i)u^-(p_j)\,, \qquad[ i | H | j\rangle = j_\mu^+(p_i, p_j)p_H^\mu\,.
 \end{equation}
 Without loss of generality, we consider only the case where the incoming
 gluon has positive helicity. The remaining helicity configurations can
 be obtained through parity transformation.
 
 Labelling the effective current by the helicities of the gluons we obtain
 for the same-helicity case
 \begin{equation}
   \label{eq:jH_same_helicity}
   \begin{split}
     j_{H,\mu}^{++}{}&(p_1,p_a,p_H) =
     \frac{m^2}{\pi v}\bigg[\\
       &-\sqrt{\frac{2p_1^-}{p_a^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spb a.1}\epsilon^{+,*}_{H,\mu}(p_1)
       +\sqrt{\frac{2p_a^-}{p_1^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spa 1.a}\epsilon^{+}_{H,\mu}(p_a)\\
       &+ [1|H|a\rangle \bigg(
         \frac{\sqrt{2}}{\spa 1.a}\epsilon^{+}_{H,\mu}(p_a)
         +         \frac{\sqrt{2}}{\spb a.1}\epsilon^{+,*}_{H,\mu}(p_1)-\frac{\spa 1.a T_2(p_a, p_a-p_H)}{\sqrt{2}(p_a-p_H)^2}\epsilon^{+,*}_{\mu}(p_1)\\
         &
         \qquad
         -\frac{\spb a.1 T_2(p_1+p_H,
           p_1)}{\sqrt{2}(p_1+p_H)^2}\epsilon^{+}_{\mu}(p_a)-\frac{RH_4}{\sqrt{2}\spb a.1}\epsilon^{+,*}_{\mu}(p_1)+\frac{RH_5}{\sqrt{2}\spa 1.a}\epsilon^{+}_{\mu}(p_a)
         \bigg)\\
         &
         - \frac{[1|H|a\rangle^2}{2 t_1}(p_{a,\mu} RH_{10} - p_{1,\mu} RH_{12})\bigg]
   \end{split}
 \end{equation}
 with $t_1 = (p_a-p_1)^2$, $t_2 = (p_a-p_1-p_H)^2$ and $R = 8 \pi^2$. Eq.~\eqref{eq:jH_same_helicity}
 is implemented by \lstinline!g_gH_HC! in \texttt{src/currents.cc}
 \footnote{\lstinline!g_gH_HC! and \lstinline!g_gH_HNC! includes an additional
 $1/t_2$ factor, which should be in the Matrix element instead.}.
 
 The currents with a helicity flip is given through
 \begin{equation}
   \label{eq:jH_helicity_flip}
   \begin{split}
     j_{H,\mu}^{+-}{}&(p_1,p_a,p_H) =
     \frac{m^2}{\pi v}\bigg[\\
     &-\sqrt{\frac{2p_1^-}{p_a^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spb
       a.1}\epsilon^{-,*}_{H,\mu}(p_1)
     +\sqrt{\frac{2p_a^-}{p_1^-}}\frac{p_{1\perp}}{|p_{1\perp}|}\frac{t_2}{\spb a.1}\epsilon^{+}_{H,\mu}(p_a)\\
     &+ [1|H|a\rangle \left(
     \frac{\sqrt{2}}{\spb a.1} \epsilon^{-,*}_{H,\mu}(p_1)
     -\frac{\spa 1.a T_2(p_a, p_a-p_H)}{\sqrt{2}(p_a-p_H)^2}\epsilon^{-,*}_{\mu}(p_1)
     - \frac{RH_4}{\sqrt{2}\spb a.1}\epsilon^{-,*}_{\mu}(p_1)\right)
     \\
     &+ [a|H|1\rangle \left(
       \frac{\sqrt{2}}{\spb a.1}\epsilon^{+}_{H,\mu}(p_a)
       -\frac{\spa 1.a
         T_2(p_1+p_H,p_1)}{\sqrt{2}(p_1+p_H)^2}\epsilon^{+}_{\mu}(p_a)
       +\frac{RH_5}{\sqrt{2}\spb a.1}\epsilon^{+}_{\mu}(p_a)
     \right)\\
     & - \frac{[1|H|a\rangle [a|H|1\rangle}{2 \spb a.1 ^2}(p_{a,\mu} RH_{10} - p_{1,\mu}
     RH_{12})\\
     &+ \frac{\spa 1.a}{\spb a.1}\bigg(RH_1p_{1,\mu}-RH_2p_{a,\mu}+2
     p_1\cdot p_H \frac{T_2(p_1+p_H, p_1)}{(p_1+p_H)^2} p_{a,\mu}
     \\
     &
     \qquad- 2p_a \cdot p_H \frac{T_2(p_a, p_a-p_H)}{(p_a-p_H)^2} p_{1,\mu}+ T_1(p_a-p_1, p_a-p_1-p_H)\frac{(p_1 + p_a)_\mu}{t_1}\\
     &\qquad-\frac{(p_1+p_a)\cdot p_H}{t_1} T_2(p_a-p_1, p_a-p_1-p_H)(p_1 - p_a)_\mu
     \bigg)
     \bigg]\,,
   \end{split}
 \end{equation}
 and implemented by \lstinline!g_gH_HNC! again in \texttt{src/currents.cc}.
 If we instead choose the gluon momentum in the $+$ direction, so that
 $p_a^- = p_{a\perp} = 0$, the corresponding currents are obtained by
 replacing $p_1^- \to p_1^+, p_a^- \to p_a^+,
 \frac{p_{1\perp}}{|p_{1\perp}|} \to -1$ in the second line of
 eq.~\eqref{eq:jH_same_helicity} and eq.~\eqref{eq:jH_helicity_flip} (see variables \lstinline!ang1a! and \lstinline!sqa1! in the implementation).
 
  The form factors $H_1,H_2,H_4,H_5, H_{10}, H_{12}$ are given in~\cite{DelDuca:2003ba}, and are implemented under their name in \texttt{src/currents.cc}. They reduce down to fundamental QCD integrals, which are again provided by \QCDloop.
 
 \subsubsection{Peripheral Higgs emission - Infinite top mass}
 \label{sec:jH_eff}
 
 To get the result with infinite top mass we could either take the limit $m_t\to \infty$ in~\eqref{eq:jH_helicity_flip} and~\eqref{eq:jH_same_helicity}, or use the \textit{impact factors} as given in~\cite{DelDuca:2003ba}. Both methods are equivalent, and lead to the same result. For the first one would find
 \begin{align}
 \lim_{m_t\to\infty} m_t^2 H_1    &= i \frac{1}{24 \pi^2}\\
 \lim_{m_t\to\infty} m_t^2 H_2    &=-i \frac{1}{24 \pi^2}\\
 \lim_{m_t\to\infty} m_t^2 H_4    &= i \frac{1}{24 \pi^2}\\
 \lim_{m_t\to\infty} m_t^2 H_5    &=-i \frac{1}{24 \pi^2}\\
 \lim_{m_t\to\infty} m_t^2 H_{10} &= 0 \\
 \lim_{m_t\to\infty} m_t^2 H_{12} &= 0.
 \end{align}
 \todo{double check this, see James thesis eq. 4.33}
 However only the second method is implemented in the code through \lstinline!C2gHgp!
 and \lstinline!C2gHgm! inside \texttt{src/currents.cc}, each function
 calculates the square of eq. (4.23) and (4.22) from~\cite{DelDuca:2003ba} respectively.
 
+\subsection{$W$+Jets}
+\subsubsection{Quark+$W$}
+
+\begin{figure}
+  \centering
+  \begin{minipage}[b]{0.2\textwidth}
+    \includegraphics[width=\textwidth]{figures/Wbits.pdf}
+  \end{minipage}
+  \begin{minipage}[b]{0.1\textwidth}
+    \centering{=}
+    \vspace{0.7cm}
+  \end{minipage}
+  \begin{minipage}[b]{0.2\textwidth}
+    \includegraphics[width=\textwidth]{figures/Wbits2.pdf}
+  \end{minipage}
+  \begin{minipage}[b]{0.1\textwidth}
+    \centering{+}
+    \vspace{0.7cm}
+  \end{minipage}
+  \begin{minipage}[b]{0.2\textwidth}
+    \includegraphics[width=\textwidth]{figures/Wbits3.pdf}
+  \end{minipage}
+  \caption{The $j_W$ current is constructed from the two diagrams which
+    contribute to the emission of a $W$-boson from a given quark line.}
+  \label{fig:jW}
+\end{figure}
+
+For a $W$ emission we require a fermion. The $j_W$ current is actually a sum of
+two separate contributions, see figure~\ref{fig:jW}, one with a $W$-emission
+from the initial state, and one with the $W$-emission from the final state.
+Mathematically this can be seen as the following two terms:
+
+\begin{align}
+  \label{eq:Weffcur1}
+  j_W^\mu(p_a,p_{\ell},p_{\bar{\ell}}, p_1) =&\ \frac{g_W^2}{2}\
+     \frac1{p_W^2-M_W^2+i\ \Gamma_W M_W}\ \bar{u}^-(p_\ell) \gamma_\alpha
+                                               v^-(p_{\bar\ell})\nonumber \\
+& \cdot \left( \frac{ \bar{u}^-(p_1) \gamma^\alpha (\slashed{p}_W +
+  \slashed{p}_1)\gamma^\mu u^-(p_a)}{(p_W+p_1)^2} +
+\frac{ \bar{u}^-(p_1)\gamma^\mu (\slashed{p}_a + \slashed{p}_W)\gamma^\alpha u^-(p_a)}{(p_a-p_W)^2} \right).
+\end{align}
+
+There are a couple of subtleties here. There is a minus sign
+distinction between the quark-anti-quark cases due to the fermion flow
+of the propagator in the current. Note that the type of $W$ emission
+(+ or -) will depend on the quark flavour, and that the handedness of
+the quark-line is given by whether its a quark or anti-quark.
+
+\subsubsection{$W$+uno}
+\begin{figure}
+  \centering
+  \begin{subfigure}{0.45\textwidth}
+  \centering
+  \includegraphics{figures/wuno1}
+  \caption{}
+  \label{fig:U1diags}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+  \centering
+\includegraphics{figures/wuno2}
+\caption{}
+\label{fig:U2diags}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+  \centering
+  \includegraphics{figures/wuno3}
+\caption{}
+\label{fig:Cdiags}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{figures/wuno4}
+\caption{}
+ \label{fig:Ddiags}
+ \end{subfigure}
+  \vspace{0.4cm}
+  \caption{Examples of each of the four categories of Feynman diagram which
+    contribute to at leading-order; there are twelve in total.  (a) is an example where the gluon and $W$
+    boson are emitted from the same quark line and the gluon comes after the
+    $t$-channel propagator.  In (b), the gluon and $W$ boson are emitted from
+    the same quark line and the gluon comes before the $t$-channel proagator.
+    In (c) the gluon is emitted from the $t$-channel gluon and in (d) the gluon
+    is emitted from the $b$--$3$ quark line.}
+  \label{fig:Wunodiags}
+\end{figure}
+It is necessary to include subleading processes in $W$+Jets also. All of
+these currents have been built for the \lstinline!Tensor! Class detailed in
+section~\ref{sec:tensor}. Similarly to the pure jet case, the uno currents are
+not calculated separately, and only in the ME functions when required
+in the \texttt{src/Wjets.cc} file. For unordered emissions a new
+current is required, $j_{W,{\rm uno}}$, it is only non-zero for
+$h_a=h_1=-$ and hence we have suppressed its helicity indices.  It is
+derived from the 12 leading-order Feynman diagrams in the QMRK
+limit (see figure~\ref{fig:Wunodiags}). Using $T^m_{ij}$ represent fundamental
+colour matrices between quark state $i$ and $j$ with adjoint index $m$ we find
+\begin{align}\label{eq:wunocurrent}
+\begin{split}
+  j^{d\,\mu}_{\rm W,uno}(p_a,p_1,p_2,p_\ell,p_{\bar{\ell}}) =& \ i \varepsilon_{\nu}(p_1)\
+  \bar{u}^-(p_\ell) \gamma_\rho v^-(p_{\bar \ell}) \\ & \quad \times\
+  \left(T^1_{2i} T^d_{ia} (\tilde U_1^{\nu\mu\rho}-\tilde L^{\nu\mu\rho}) +
+    T^d_{2i} T^1_{ia} (\tilde U_2^{\nu\mu\rho}+\tilde L^{\nu\mu\rho}) \right),
+  \end{split}
+\end{align}
+where expressions for $\tilde U_{1,2}^{\nu\mu\rho}$ and $\tilde L^{\nu\mu\rho}$
+are given as:
+\begin{align}
+  \label{eq:U1tensor}
+  \begin{split}
+    \tilde U_1^{\nu\mu\rho} =&\frac{\langle 2|\nu (\slashed{p}_2+ \slashed{p}_1)\mu (\slashed{p}_a - \slashed{p}_W)\rho P_L |a\rangle }{s_{12}t_{aW}}  + \frac{\langle 2|\nu (\slashed{p}_2+ \slashed{p}_1)\rho P_L (\slashed{p}_2+\slashed{p}_1 + \slashed{p}_W)\mu |a\rangle }{s_{12}s_{12W}} \\
+    &+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \nu
+      (\slashed{p}_1 + \slashed{p}_2+\slashed{p}_W)\mu |a\rangle}{s_{2W}s_{12W}}.
+  \end{split}
+\end{align}
+\begin{align}
+  \label{eq:U2tensor}
+  \begin{split}
+    \tilde U_2^{\nu\mu\rho} =&\frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W-\slashed{p}_1)\nu (\slashed{p}_a - \slashed{p}_W)\rho P_L |a\rangle }{t_{aW1}t_{aW}}  + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W- \slashed{p}_1)\rho P_L (\slashed{p}_a-\slashed{p}_1) \nu |a\rangle }{t_{a1W}t_{a1}}  \\
+    &+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \mu
+      (\slashed{p}_a-\slashed{p}_1)\nu |a\rangle}{s_{2W}t_{a1}}.
+  \end{split}
+\end{align}
+\begin{align}
+\label{eq:Ltensor}
+ \tilde L^{\nu\mu\rho} &= \frac{q_2^2}{2t_{aW2}}  \left[\frac{\langle 2 |\mu (\slashed{p}_a - \slashed{p}_W) \rho P_L | a\rangle}{t_{aW}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \mu | a\rangle }{s_{2W}} \right]
+ \cdot \left( \frac{p_b^{\nu}}{p_b\cdot p_1} + \frac{p_3^{\nu}}{p_3\cdot p_1} \right) \nonumber \\
+ &\quad+\frac{1}{t_{aW2}}\left[\frac{\langle 2 |\sigma (\slashed{p}_a - \slashed{p}_W) \rho P_L | a\rangle}{t_{aW}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \sigma | a\rangle }{s_{2W}} \right] \nonumber \\
+ &\vphantom{+\frac{1}{t_{aW2}}}\quad\cdot \left( g^{\sigma \mu} (q_1 +q_2)^\nu + g^{\mu \nu}(-q_2 +p_1)^\sigma+ g^{\nu \sigma}(-p_1 -q_1)^\mu \right).
+\end{align}
+
+\subsubsection{$W$+Extremal $\mathbf{q\bar{q}}$}
+The $W$+Jet sub-leading processes containing an extremal $q\bar{q}$ are
+related by crossing symmetry to the $W$+Jet unordered processes. This
+means that one can simply perform a crossing symmetry argument on
+eq.~\ref{eq:wunocurrent} to arrive at the extremal $q\bar{q}$ current
+required.We show the basic structure of the extremal $q\bar{q}$
+current in figure~\ref{fig:qgimp}, neglecting the $W$-emission for
+simplicity.
+
+\begin{figure}
+\centering
+\includegraphics[width=0.3\textwidth]{{qqbarex_schem}}
+\caption{Schematic structure of the $gq \to \bar{Q}Qq$ amplitude in the limit
+  $y_1 \sim y_2 \ll y_3$}
+\label{fig:qgimp}
+\end{figure}
+
+\begin{figure}
+\centering
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{qqbarex1}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{qqbarex2}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{qqbarex4}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{qqbarex5}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{qqbarex3}
+\end{subfigure}
+\caption{The five tree-level graphs which contribute to the process $gq \to \bar{Q}Qq$.}
+\label{fig:qg_qQQ_graphs}
+\end{figure}
+We can obtain the current for $g\rightarrow W q \bar{q}$ by evaluating
+the current for $W$ plus unordered emissions with the normal arguments
+$p_a \leftrightarrow -p_1 $ interchanged.  This is a non-trivial
+statement: due to the minimality of the approximations made, the
+crossing symmetry normally present in the full amplitude may be
+extended to the factorised current.
+
+We must again note that swapping $p_a \leftrightarrow -p_1$ will lead
+to $u$-spinors with momenta with negative energy.  These are identical
+to $v$-spinors with momenta with positive energy, up to an overall
+phase which is common to all terms, and can therefore be neglected.
+
+Mathematically, this is given by:
+
+\begin{align}\label{eq:crossedJ}
+  j^\mu_{\rm W,g\to Q\bar{Q}}(p_a,p_1,p_2,p_\ell,p_{\bar{\ell}}) =i \varepsilon_{g\nu}
+  \langle \ell | \rho | \bar \ell \rangle_L
+\left(T^1_{2i} T^d_{ia} (\tilde U_{1,X}^{\nu\mu\rho}-\tilde L^{\nu\mu\rho}_X) + T^d_{2i} T^1_{ia} (\tilde U_{2,X}^{\nu\mu\rho}+\tilde L_X^{\nu\mu\rho}) \right),
+\end{align}
+where the components are now given by
+\begin{align}
+  \label{eq:U1tensorX}
+  \begin{split}
+    \tilde U_{1,X}^{\nu\mu\rho} =&\frac{\langle 2|\nu (\slashed{p}_a- \slashed{p}_2)\mu (\slashed{p}_1 + \slashed{p}_W)\rho P_L |1\rangle }{t_{a2}s_{1W}}  + \frac{\langle 2|\nu (\slashed{p}_a- \slashed{p}_2)\rho P_L (\slashed{p}_a-\slashed{p}_2 - \slashed{p}_W)\mu |1\rangle }{t_{a2}t_{a2W}} \\
+    &- \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \nu
+      (\slashed{p}_a - \slashed{p}_2-\slashed{p}_W)\mu
+      |1\rangle}{s_{2W}t_{a2W}}, \\
+    \tilde U_{2,X}^{\nu\mu\rho} =&-\frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W-\slashed{p}_1)\nu (\slashed{p}_1 + \slashed{p}_W)\rho P_L |1\rangle }{t_{aW1}s_{1W}}  + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W- \slashed{p}_1)\rho P_L (\slashed{p}_a-\slashed{p}_1) \nu |1\rangle }{t_{a1W}t_{a1}}  \\
+    &+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \mu
+      (\slashed{p}_a-\slashed{p}_1)\nu |1\rangle}{s_{2W}t_{a1}}, \\
+    \tilde L^{\nu\mu\rho}_X &= \frac{q_2^2}{2s_{1W2}}  \left[\frac{\langle 2 |\mu (\slashed{p}_1 + \slashed{p}_W) \rho P_L | 1\rangle}{s_{1W}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \mu | 1\rangle }{s_{2W}} \right]
+    \cdot \left( \frac{p_b^{\nu}}{p_a\cdot p_b} + \frac{p_3^{\nu}}{p_a\cdot p_3} \right) \\
+    &\quad+\frac{1}{s_{W12}}\left[-\frac{\langle 2 |\sigma (\slashed{p}_1 + \slashed{p}_W) \rho P_L | a\rangle}{t_{aW}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \sigma | 1\rangle }{s_{2W}} \right] \\
+    &\vphantom{+\frac{1}{t_{aW2}}}\quad\cdot \left( g^{\sigma \mu} (q_1 +q_2)^\nu + g^{\mu \nu}(-q_2 +p_1)^\sigma+ g^{\nu \sigma}(-p_1 -q_1)^\mu \right).
+  \end{split}
+\end{align}
+Notice in particular the similarity to the $W$+uno scenario (from which
+this has been derived).
+\subsubsection{Central $\mathbf{q\bar{q}}$ Vertex}
+The final subleading process in the $W$+Jet case is the Central
+$q\bar{q}$ vertex. This subleading process does not require an altered
+current, but an effective vertex which is contracted with two regular
+\HEJ currents. This complexity is dealt with nicely by the \lstinline!Tensor!
+class, which is detailed in section~\ref{sec:tensor}.
+
+The $W$-emission can be from the central effective vertex (scenario
+dealt with by the function \texttt{jM2WqqtoqQQq()} in the file
+\texttt{src/Wjets.cc}); or from either of the external quark legs
+(scenario dealt with by \texttt{jM2WqqtoqQQqW()} in same file).  In
+the pure jets case, there are 7 separate diagrams which contribute to
+this, which can be seen in figure~\ref{fig:qq_qQQq_graphs}. In the $W$+Jets
+case, there are then 45 separate contributions.
+
+\begin{figure}
+\centering
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{figures/qqbarcen1}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{figures/qqbarcen2}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{figures/qqbarcen3}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{figures/qqbarcen4}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{figures/qqbarcen5}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{figures/qqbarcen6}
+\end{subfigure}
+\begin{subfigure}{0.45\textwidth}
+\centering
+\includegraphics{figures/qqbarcen7}
+\end{subfigure}
+\caption{All Feynman diagrams which contribute to $qq' \to qQ\bar{Q}q'$ at
+  leading order.}
+\label{fig:qq_qQQq_graphs}
+\end{figure}
+
+The end result is of the effective vertex, after derivation, is:
+
+\begin{align}
+  \label{eq:EffectiveVertexqqbar}
+  \begin{split}
+  V^{\mu\nu}_{\text{Eff}}=&
+  \frac{C_1}{s_{23AB}}\left(X^{\mu\nu\sigma}_{1a}\hat{t_1} + X^{\mu\nu\sigma}_{4b}\hat{t_3} +V^{\mu\nu\sigma}_{3g}\right)J_{\text{V} \sigma}(p_2,p_A,p_B,p_3)
+  \\
+  &\quad +iC_2X^{\mu\nu}_{Unc}+iC_3X^{\mu\nu}_{Cro},
+  \end{split}
+\end{align}
+
+where:
+
+\begin{eqnarray}
+  \begin{split}
+  C_1=&T^e_{1q}T^g_{qa}T^e_{23}T^g_{4b} -
+  T^g_{1q}T^e_{qa}T^e_{23}T^g_{4b} = f^{egc}T^c_{1a}T^e_{23}T^g_{4b},
+  \\
+  C_2=&T^g_{1a}T^g_{2q}T^{g'}_{q3}T^{g'}_{4b}
+  \\
+  C_3=&T^g_{1a}T^{g'}_{2q}T^g_{q3}T^{g'}_{4b}
+  \end{split}
+\end{eqnarray}
+are the colour factors of different contributions. The following
+tensor structures correspond to groupings of diagrams in
+figure~\ref{fig:qq_qQQq_graphs}.
+\begin{eqnarray}
+  \label{eq:1aFixed}
+  X_{1a}^{\mu\nu\sigma} &=
+\frac{-g^{\mu\nu}}{s_{23AB}\hat{t_3}}\left(\frac{p^\sigma_a}{(s_{a23AB})} +
+                    \frac{p^\sigma_1}{(s_{123AB})}\right)
+  \\
+  \label{eq:4bFixed}
+  X_{4b}^{\mu\nu\sigma}
+&=\frac{g^{\mu\nu}}{s_{23AB}\hat{t_1}}\left(\frac{p^\sigma_b}{(s_{23bAB})}+\frac{p^\sigma_4}{(s_{234AB}}\right)
+\end{eqnarray}
+correspond to the first and second row of diagrams in figure~\ref{fig:qq_qQQq_graphs}.
+\begin{align}
+  \label{eq:3GluonWEmit}
+  \begin{split}
+  X^{\mu\nu}_{3g}=\frac{1}{
+    \hat{t}_1s_{23AB}\,\hat{t}_3}
+  \bigg[&\left(q_1+p_2+p_3+p_A+p_B\right)^\nu
+  g^{\mu\sigma}+
+  \\
+  &\quad\left(q_3-p_2-p_3-p_A-p_B\right)^\mu g^{\sigma\nu}-
+  \\
+  & \qquad\left(q_1+q_3\right)^\sigma g^{\mu\nu}\bigg]J_{\text{V} \sigma}(p_2,p_A,p_B,p_3)
+  \end{split}
+\end{align}
+corresponds to the left diagram on the third row in
+figure~\ref{fig:qq_qQQq_graphs}. One notes that all of these contributions have
+the same colour factor, and as such we can group them together nicely
+before summing over helicities etc. As such, the function
+\texttt{MSymW()} in \texttt{src/Wjets.cc} returns a \lstinline!Tensor!
+containing the information from these 5 groupings of contributions (30 diagrams
+in total).
+\begin{align}
+  \begin{split}
+      X^{\mu\nu}_{Unc}=\frac{\langle A|\sigma P_L|B\rangle}{\hat{t_1}\hat{t_3}} \bar{u}_2&\left[-\frac{
+  \gamma^\sigma P_L(\slashed{p}_2+\slashed{p}_A+\slashed{p}_B)\gamma^\mu
+  (\slashed{q}_3+ \slashed{p}_3)\gamma^\nu}{(s_{2AB})(t_{unc_{2}})}\right.+
+  \\
+    &\qquad\left. \frac{\gamma^\mu(\slashed{q}_1-\slashed{p}_2)\gamma^\sigma
+        P_L(\slashed{q}_3+\slashed{p}_3)\gamma^\nu}{(t_{unc_{1}})(t_{unc_{2}})}\right. +
+  \\
+    &\qquad\qquad\left. \frac{\gamma^\mu(\slashed{q}_1-\slashed{p}_2)\gamma^\nu(\slashed{p}_3+\slashed{p}_A+\slashed{p}_B)\gamma^\sigma P_L
+  }{(t_{unc_1})(s_{3AB})}\right]v_3
+  \end{split}
+\end{align}
+corresponds to the diagram on the right of row three in
+figure~\ref{fig:qq_qQQq_graphs}. This contribution to the effective vertex can
+be obtained in the code with the function \texttt{MUncW()} in file
+\texttt{src/Wjets.cc}.
+\begin{align}
+  \begin{split}
+    X^{\mu\nu}_{Cro}=\frac{\langle
+A|\sigma P_L|B\rangle}{\hat{t_1}\hat{t_3}} \bar{u}_2&\left[-\frac{
+\gamma^\nu(\slashed{q}_3+\slashed{p}_2)\gamma^\mu
+(\slashed{p}_3+\slashed{p}_A+\slashed{p}_B)\gamma^\sigma P_L}{(t_{cro_1})(s_{3AB})}\right.+
+\\
+&\qquad\left. \frac{\gamma^\nu(\slashed{q}_3+\slashed{p}_2)\gamma^\sigma
+    P_L(\slashed{q}_1-\slashed{p}_3)\gamma^\mu}{(t_{cro_{1}})(t_{cro_{2}})}\right.+
+\\ &\qquad\qquad\left
+. \frac{\gamma^\sigma P_L(\slashed{p}_2+\slashed{p}_A+\slashed{p}_B)\gamma^\nu(\slashed{q}_1-\slashed{p}_3)\gamma^\mu
+}{(s_{2AB})(t_{cro_2})}\right]v_3
+  \end{split}
+\end{align}
+corresponds to the last diagram in figure~\ref{fig:qq_qQQq_graphs}. This
+contribution to the effective vertex can be obtained in the code with
+the function \texttt{MCroW()} in file \texttt{src/Wjets.cc}.
 
 %%% Local Variables:
 %%% mode: latex
 %%% TeX-master: "developer_manual"
 %%% End:
diff --git a/doc/developer_manual/developer_manual.tex b/doc/developer_manual/developer_manual.tex
index 2729cce..387e65a 100644
--- a/doc/developer_manual/developer_manual.tex
+++ b/doc/developer_manual/developer_manual.tex
@@ -1,1582 +1,1584 @@
 \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)
 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 $W$ and $Z/\gamma^*$ production together with jets are
 given in~\cite{Andersen:2012gk,Andersen:2016vkp}, but not yet included.
 
 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!). In case there is also
 an unordered gluon emission the matrix element is already suppressed by one
 logarithm $\log(s/t)$. Therefore at NLL the Higgs boson can only be emitted
 centrally\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/developer_manual/figures/ColourConnects_FKL.pdf b/doc/developer_manual/figures/ColourConnects_FKL.pdf
index f011d62..0db27cf 100755
Binary files a/doc/developer_manual/figures/ColourConnects_FKL.pdf and b/doc/developer_manual/figures/ColourConnects_FKL.pdf differ
diff --git a/doc/developer_manual/figures/ColourConnects_crossed.pdf b/doc/developer_manual/figures/ColourConnects_crossed.pdf
index e161928..d360dce 100755
Binary files a/doc/developer_manual/figures/ColourConnects_crossed.pdf and b/doc/developer_manual/figures/ColourConnects_crossed.pdf differ
diff --git a/doc/developer_manual/figures/Wbits.pdf b/doc/developer_manual/figures/Wbits.pdf
new file mode 100644
index 0000000..fb61e4d
--- /dev/null
+++ b/doc/developer_manual/figures/Wbits.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4b2935ead78945daf81d185eaec4ff04279f3e949c43353f31fcd4ce736d721e
+size 3173
diff --git a/doc/developer_manual/figures/Wbits2.pdf b/doc/developer_manual/figures/Wbits2.pdf
new file mode 100644
index 0000000..9929ee1
--- /dev/null
+++ b/doc/developer_manual/figures/Wbits2.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9db642b30f29c8598f456c74c8a16af3f149e24e0a98c09832105c6f877598be
+size 3161
diff --git a/doc/developer_manual/figures/Wbits3.pdf b/doc/developer_manual/figures/Wbits3.pdf
new file mode 100644
index 0000000..79a20fc
--- /dev/null
+++ b/doc/developer_manual/figures/Wbits3.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c97e6c7e2565194c45ad38395a5e098fa62298bfc1c8445be5df42477fbce2cd
+size 3482
diff --git a/doc/developer_manual/figures/colour_centralqqx.jpg b/doc/developer_manual/figures/colour_centralqqx.jpg
index 84ea16a..194a571 100644
Binary files a/doc/developer_manual/figures/colour_centralqqx.jpg and b/doc/developer_manual/figures/colour_centralqqx.jpg differ
diff --git a/doc/developer_manual/figures/colour_gleft.jpg b/doc/developer_manual/figures/colour_gleft.jpg
index 9a24558..7b38379 100644
Binary files a/doc/developer_manual/figures/colour_gleft.jpg and b/doc/developer_manual/figures/colour_gleft.jpg differ
diff --git a/doc/developer_manual/figures/colour_gright.jpg b/doc/developer_manual/figures/colour_gright.jpg
index 32e5179..7c65f02 100644
Binary files a/doc/developer_manual/figures/colour_gright.jpg and b/doc/developer_manual/figures/colour_gright.jpg differ
diff --git a/doc/developer_manual/figures/colour_qx.jpg b/doc/developer_manual/figures/colour_qx.jpg
index 344bbc0..8372340 100644
Binary files a/doc/developer_manual/figures/colour_qx.jpg and b/doc/developer_manual/figures/colour_qx.jpg differ
diff --git a/doc/developer_manual/figures/colour_uno.jpg b/doc/developer_manual/figures/colour_uno.jpg
index 362686c..2e50c63 100644
Binary files a/doc/developer_manual/figures/colour_uno.jpg and b/doc/developer_manual/figures/colour_uno.jpg differ
diff --git a/doc/developer_manual/figures/dR_2.pdf b/doc/developer_manual/figures/dR_2.pdf
index 5a8745c..0eb9f05 100644
Binary files a/doc/developer_manual/figures/dR_2.pdf and b/doc/developer_manual/figures/dR_2.pdf differ
diff --git a/doc/developer_manual/figures/dR_2_small.pdf b/doc/developer_manual/figures/dR_2_small.pdf
index 4880eea..00bcb44 100644
Binary files a/doc/developer_manual/figures/dR_2_small.pdf and b/doc/developer_manual/figures/dR_2_small.pdf differ
diff --git a/doc/developer_manual/figures/ng_1.5.pdf b/doc/developer_manual/figures/ng_1.5.pdf
index d405ea8..2ce30de 100644
Binary files a/doc/developer_manual/figures/ng_1.5.pdf and b/doc/developer_manual/figures/ng_1.5.pdf differ
diff --git a/doc/developer_manual/figures/ng_5.5.pdf b/doc/developer_manual/figures/ng_5.5.pdf
index 0435138..dbf676e 100644
Binary files a/doc/developer_manual/figures/ng_5.5.pdf and b/doc/developer_manual/figures/ng_5.5.pdf differ
diff --git a/doc/developer_manual/figures/ng_mean.pdf b/doc/developer_manual/figures/ng_mean.pdf
index b4c4097..bc8eeed 100644
Binary files a/doc/developer_manual/figures/ng_mean.pdf and b/doc/developer_manual/figures/ng_mean.pdf differ
diff --git a/doc/developer_manual/figures/pJ.pdf b/doc/developer_manual/figures/pJ.pdf
index 4db6d25..0089a57 100644
Binary files a/doc/developer_manual/figures/pJ.pdf and b/doc/developer_manual/figures/pJ.pdf differ
diff --git a/doc/developer_manual/figures/ptBMin.pdf b/doc/developer_manual/figures/ptBMin.pdf
index 607cc20..6397d77 100644
Binary files a/doc/developer_manual/figures/ptBMin.pdf and b/doc/developer_manual/figures/ptBMin.pdf differ
diff --git a/doc/developer_manual/figures/qqbarcen1.pdf b/doc/developer_manual/figures/qqbarcen1.pdf
new file mode 100644
index 0000000..b180327
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarcen1.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:782c11e4a4eb6b491751af1f6897d7a055f3d4a1df8e4421de0716b81c55a780
+size 5169
diff --git a/doc/developer_manual/figures/qqbarcen2.pdf b/doc/developer_manual/figures/qqbarcen2.pdf
new file mode 100644
index 0000000..3573495
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarcen2.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fb394c8ddf110f0353d2f8c165a413208e1b0755da5d4147ca1437d188cca83c
+size 5322
diff --git a/doc/developer_manual/figures/qqbarcen3.pdf b/doc/developer_manual/figures/qqbarcen3.pdf
new file mode 100644
index 0000000..1bd8c12
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarcen3.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7dca47fbd26056dbd09eb186f3912c0c7a39399c4a624fb6fac7ccc728cc36f4
+size 5162
diff --git a/doc/developer_manual/figures/qqbarcen4.pdf b/doc/developer_manual/figures/qqbarcen4.pdf
new file mode 100644
index 0000000..32d020a
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarcen4.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6315caf833c4e19d041ec29a2f0e38d43ad367a56130777132a408ff6079bae8
+size 5329
diff --git a/doc/developer_manual/figures/qqbarcen5.pdf b/doc/developer_manual/figures/qqbarcen5.pdf
new file mode 100644
index 0000000..93d176f
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarcen5.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:20be4cdf7b79395e811cee155d933cd8b5b18e232dbb32cda71984d06595a40c
+size 4976
diff --git a/doc/developer_manual/figures/qqbarcen6.pdf b/doc/developer_manual/figures/qqbarcen6.pdf
new file mode 100644
index 0000000..0d24cc5
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarcen6.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3a1cb40f4cc2d875f4198db6ceea33f770a0ea6b2204ef53d80c6191276c56be
+size 4834
diff --git a/doc/developer_manual/figures/qqbarcen7.pdf b/doc/developer_manual/figures/qqbarcen7.pdf
new file mode 100644
index 0000000..7586171
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarcen7.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:bc813739ba3e07832f474dfd07c0cc28f924e1b707434d18e06ebffd4d9bafbf
+size 4872
diff --git a/doc/developer_manual/figures/qqbarex1.pdf b/doc/developer_manual/figures/qqbarex1.pdf
new file mode 100644
index 0000000..41a778a
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarex1.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2230151e2eeb50f773c19218fc24a83d9d85d0e59e4fcc9ac26202abb478fda8
+size 5101
diff --git a/doc/developer_manual/figures/qqbarex2.pdf b/doc/developer_manual/figures/qqbarex2.pdf
new file mode 100644
index 0000000..fdcf011
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarex2.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ffd952ae26d5fff4f4888b001fbe198dcdef2a543a97c3ccc4171dd9d27569c0
+size 5285
diff --git a/doc/developer_manual/figures/qqbarex3.pdf b/doc/developer_manual/figures/qqbarex3.pdf
new file mode 100644
index 0000000..ee0526c
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarex3.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ad7a129e3820d2428e0bd4bdbad26222ce77e8f1a44885c8126fc0025c8b86b4
+size 5014
diff --git a/doc/developer_manual/figures/qqbarex4.pdf b/doc/developer_manual/figures/qqbarex4.pdf
new file mode 100644
index 0000000..adb1011
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarex4.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a7c1e11dff2270760423a87e149babdffc708bebb0793a55a6e529e3f2966371
+size 4840
diff --git a/doc/developer_manual/figures/qqbarex5.pdf b/doc/developer_manual/figures/qqbarex5.pdf
new file mode 100644
index 0000000..845a46a
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarex5.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6bd05346c376d77e44b5ac4022424632e2dfa365e20e475031a08a5a4819040d
+size 4876
diff --git a/doc/developer_manual/figures/qqbarex_schem.pdf b/doc/developer_manual/figures/qqbarex_schem.pdf
new file mode 100644
index 0000000..0dce5f6
--- /dev/null
+++ b/doc/developer_manual/figures/qqbarex_schem.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4a275254f30d614606c1ae85dd886c8fd87819b518fa2c61b657d9f699cd8444
+size 6814
diff --git a/doc/developer_manual/figures/wuno1.pdf b/doc/developer_manual/figures/wuno1.pdf
new file mode 100644
index 0000000..1110954
--- /dev/null
+++ b/doc/developer_manual/figures/wuno1.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2a13a926a7c6c3468fb09d6a968a4e0cf321cf26db28f276c1e8e2531bc4c0f6
+size 5247
diff --git a/doc/developer_manual/figures/wuno2.pdf b/doc/developer_manual/figures/wuno2.pdf
new file mode 100644
index 0000000..5b1d742
--- /dev/null
+++ b/doc/developer_manual/figures/wuno2.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:38fefdfb069a43916677f110a35aa6e65e49820d1ef664067a6499a8ae33c68b
+size 5489
diff --git a/doc/developer_manual/figures/wuno3.pdf b/doc/developer_manual/figures/wuno3.pdf
new file mode 100644
index 0000000..4aad276
--- /dev/null
+++ b/doc/developer_manual/figures/wuno3.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:76b070e9ae86f2b71f3f8eb90c16814123d5656fbfa8caad5e6496163cccafc5
+size 5602
diff --git a/doc/developer_manual/figures/wuno4.pdf b/doc/developer_manual/figures/wuno4.pdf
new file mode 100644
index 0000000..7fb6e3f
--- /dev/null
+++ b/doc/developer_manual/figures/wuno4.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5c7984ffecd1a52ad84bd4a1b89661e724d82996d49945ad3b6de4a70ebae136
+size 5594