diff --git a/current_generator/include/currents.frm b/current_generator/include/currents.frm
index afe5196..8ece75d 100644
--- a/current_generator/include/currents.frm
+++ b/current_generator/include/currents.frm
@@ -1,83 +1,83 @@
*/**
* \brief Predefined current functions and polarisation vectors
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020-2025
* \copyright GPLv2 or later
*/
#ifndef `$HEJCurrentsIncluded'
#$HEJCurrentsIncluded = 1;
* FKL currents for leptonically decaying vector boson and stable photon
cfunction JV,Jgamma;
* Polarisation vectors
cfunction Eps;
* Polarisation vectors contracted with Higgs emission vertex
cfunction EPSH;
* Higgs emission vertex
cfunction VH;
cfunction m2,plus,minus,perphat,sqrt,conj,dot;
s T1,T2;
* internal symbols are all uppercase
symbols H, H1, HL, HPHOTON, M;
indices MU1,...,MU50;
vectors PA,P1,PH,PLBAR,PL,PG,PPHOTON,PR,Q1,Q2;
#procedure InsSpecialCurrents
#call InsJV
#call InsJgamma
#call InsEps
#call InsVH
#endprocedure
* Replace leptonically decaying vector boson FKL current
#procedure InsJV()
* we use a loop here to introduce a new dummy index for each occurence
repeat;
once JV(H1?, HL?, MU1?, PA?, P1?, PLBAR?, PL?) = Current(HL, PL, MU2, PLBAR)*(
+ Current(H1, P1, MU2, P1+PLBAR+PL, MU1, PA)/m2(P1+PLBAR+PL)
+ Current(H1, P1, MU1, PA-PLBAR-PL, MU2, PA)/m2(PA-PLBAR-PL)
);
sum MU2;
endrepeat;
#endprocedure
* Replace photon FKL current
* see eq:J_photon in developer manual
#procedure InsJgamma()
* we use a loop here to introduce a new dummy index for each occurence
repeat;
once Jgamma(H1?, HPHOTON?, MU1?, PA?, P1?, PPHOTON?, PR?) = (
- + Current(H1, P1, MU2, P1+PPHOTON, MU1, PA)
- + Current(H1, P1, MU1, PA-PPHOTON, MU2, PA)
+ + Current(H1, P1, MU2, P1+PPHOTON, MU1, PA)/m2(P1+PPHOTON)
+ + Current(H1, P1, MU1, PA-PPHOTON, MU2, PA)/m2(PA-PPHOTON)
)*Eps(HPHOTON, MU2, PPHOTON, PR);
sum MU2;
endrepeat;
#endprocedure
* Replace polarisation vectors by currents
#procedure InsEps()
id Conjugate(Eps(H1?, MU1?, PG?, PR?)) = - Eps(-H1, MU1, PG, PR);
id Eps(-1, MU1?, PG?, PR?) = sqrt(2)/2*SpinorChain(PR, MU1, PG)/AngleChain(PG,PR);
id Eps(+1, MU1?, PG?, PR?) = sqrt(2)/2*SpinorChain(PG, MU1, PR)/SquareChain(PG,PR);
#endprocedure
* Higgs emission vertex
* see eq:VH in developer manual, but without overall prefactor
#procedure InsVH()
id VH(MU1?, MU2?, Q1?, Q2?) = d_(MU1, MU2)*T1 - Q2(MU1)*Q1(MU2)*T2;
#endprocedure
#endif
diff --git a/doc/developer_manual/currents.tex b/doc/developer_manual/currents.tex
index 6f84656..4ea4e73 100644
--- a/doc/developer_manual/currents.tex
+++ b/doc/developer_manual/currents.tex
@@ -1,721 +1,721 @@
\section{Currents}
\label{sec:currents_impl}
The following section contains a list of all the currents implemented in \HEJ.
Clean up of the code structure is ongoing. Each implemented current has its own
separate source file (e.g. \texttt{src/<PROCESS>.cc}), and associated header file
(e.g. \texttt{include/HEJ/<PROCESS>.hh}). The processes (and their filename) that are implemented are:
Pure jets (\texttt{jets}),
$W$+jets (\texttt{Wjets}),
$Z/\gamma$+jets (\texttt{Zjets}),
$h$+jets (\texttt{Hjets}),
$W^\pm W^\pm $+jets (\texttt{WWjets}).
The naming convention for the current contraction $\left\|S_{f_1 f_2\to f_1
f_2}\right\|^2$ is \lstinline!ME_[Boson]_[subleading-type]_[incoming]!. For
example \lstinline!ME_W_unob_qq! corresponds to the contraction $j_W^\mu
j_{\text{uno}, \mu}$ ($qQ\to \bar{q}WQg$). For bosons on the same side as the
subleading we drop the connecting underscore, e.g. \lstinline!ME_Wuno_qq!
gives $j_{W,\text{uno}}^\mu j_\mu$ ($qQ\to g\bar{q}WQ$).
\subsection{Pure Jets}
\subsubsection{Quarks and Gluons}
\label{sec:current_quark}
The basic quark current is
\begin{align}
\label{eq:j}
j^\pm_\mu(p_i,p_j)=\bar{u}^\pm(p_i)\gamma_\mu u^\pm(p_j),
\end{align}
see also equation~(\ref{eq:j_gen}). The contraction of two currents is
\begin{equation}
\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\,.
\end{equation}
Gluons use the same currents as quarks. The differences in the
tree-level matrix elements are absorbed into the colour acceleration
multipliers defined in equations (\ref{eq:K_g}) and (\ref{eq:K_q}).
\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}
using the short-hand notation $j_{ij} = j^{\lambda_1}(p_i, p_j)$.
$j^{{\rm uno}\; \mu}$ is currently not calculated as a separate
current, but always as needed for the ME (i.e. in
\lstinline!ME_unob_XX!). The contraction with a basic quark/gluon current is
\begin{equation}
\label{eq:S_uno}
\left\|S^{\text{uno}}_{q f_2\to g q f_2}\right\|^2 = 2\sum_{\substack{\lambda_a =
+,-\\\lambda_g = +,-}} \left[C_F |X + Y|^2 - C_A \Re(XY^*)\right]
\end{equation}
with
\begin{align}
\label{eq:X_munu}
X^{\mu\nu} ={}& U_1^{\mu\nu} - L^{\mu\nu}, \\
\label{eq:Y_munu}
Y^{\mu\nu} ={}& U_2^{\mu\nu} + L^{\mu\nu}, \\
\label{eq:X}
X ={}& j_\mu^-(p_b, p_2) \epsilon_\nu^{\lambda_g}(p_g, p_r) X^{\mu\nu},\\
\label{eq:Y}
Y ={}& j_\mu^-(p_b, p_2) \epsilon_\nu^{\lambda_g}(p_g, p_r) Y^{\mu\nu}.
\end{align}
In the current implementation we choose $p_r = p_1$.
\subsubsection{Extremal \texorpdfstring{$q\bar{q}$}{qqbar}}
In Pure jets we also include the subleading process which arises when
an incoming gluon splits into a $q\bar{q}$ pair. This splitting impact
factor is related to the unordered current by simple means of a
crossing symmetry. Namely, starting from the unordered current
$q(p_a) \to g(p_g) q(p_1) g^*(\tilde{q})$, we make the following replacements:
\begin{equation}
\label{eq:crossing}
p_a \to - p_{\bar{q}}\,, \qquad p_g \to - p_a\,, \qquad p_1 \to p_q\,.
\end{equation}
After crossing we end up with negative-energy spinors such as $| - p_a \rangle$.
Using analytic continuation these can be rewritten as positive-energy ones by pulling
out a phase factor $i$, e.g. $| - p_a \rangle \to i \, | p_a \rangle$~\cite{Gunion:1985vca}.
\subsubsection{Central \texorpdfstring{$q\bar{q}$}{qqbar}}
\label{sec:central_qqbar}
The final subleading process type in the Pure Jets case is Central
$q\bar{q}$. In this process type, we have two currents scattering off
of each other, but this time, via an effective vertex, which connects
together two FKL chains. Each FKL chain t-channel gluon splits into a
$q\bar{q}$ and this results in a quark and anti-quark in between the
most forward and backward jets. One can see an example of such a
process in Figure \ref{fig:qqbarcen_example}.
\begin{figure}[ht]
\centering
\includegraphics[scale=0.8]{Cenqqbar_jx}
\caption{Momentum labeling for a central $q\bar{q}$ process.}
\label{fig:qqbarcen_example}
\end{figure}
As the new central $q\bar{q}$ piece contains the quark propagator, we will treat
this as part of the skeleton process. This means that we do not impose strong ordering
between the $q\bar{q}$-pair taking
\begin{align}
\label{eq:cenqqbarraporder}
y_1 \ll y_q,y_{\bar{q}} \ll y_n.
\end{align}
The skeleton function is obtained by summing all the leading-order diagrams
contributing to this process, shown in Figure~\ref{fig:qq_qQQq_graphs}, and
dropping terms kinematically suppressed in the QMRK limit to arrive at a
factorised form. It reads~\cite{Andersen:2020yax}
\footnote{In~\cite{Andersen:2020yax} the two colour factors are rewritten as $C_A C_F^2/4=C_F$ and $C_F/4=1/C_A$ using $N_c=3$. Note the plus sign for the second term (typo in eq.~(3.39) of arXiv:2012.10310v2).\label{footnote_qqbar}}
\todo{Update footnote when corrected version is on arxiv}
\begin{equation}
\label{eq:Scen}
\left\|S_{f_a f_b \to f_a q\bar{q} f_b}^{\rm cen}\right\|^2 =
\frac{1}{2} \frac{1}{q_1^2q_3^2}\sum_{\{h_i\}}\ \left[\frac{C_A C_F^2}{4}\left(|V|^2 + |W|^2\right) +
\frac{C_F}{4} {\rm Re}\left(V W^*\right) \right],
\end{equation}
where
\begin{align}
\label{eq:Vcen}
V &= j_{\mu}(p_a,p_1) \left( X_s^{\mu\nu} + X^{\mu \nu}_{6} \right) j_{\nu}(p_b,p_n)\\
\label{eq:Wcen}
W &= j_{\mu}(p_a,p_1) \left( X_s^{\mu\nu} + X^{\mu \nu}_{7} \right) j_{\nu}(p_b,p_n)
\end{align}
with
\begin{align}
\label{eq:Xs}
X_s^{\mu\nu} &= \frac{1}{s_{q\bar{q}}}\left(\eta^{\mu \nu} X_{sym}^\sigma + X^{\mu \nu \sigma}_{5}
\right) \left[ \bar{u}_{q} \gamma_{\sigma} v_{\bar{q}}\right], \\
\label{eq:Xsym}
X_{sym}^\sigma &= q_1^2 \left(\frac{p_a^\sigma}{s_{aq}+s_{a\bar{q}}} + \frac{p_1^\sigma}{s_{1q}+s_{1\bar{q}}} \right)
- q_3^2 \left(\frac{p_b^\sigma}{s_{bq}+s_{b\bar{q}}} + \frac{p_n^\sigma}{s_{nq}+s_{n\bar{q}}} \right), \\
\label{eq:X5}
X_{5}^{\mu \nu \sigma} &= (q_1 + p_q + p_{\bar{q}})^\nu \eta^{\mu \sigma} + (q_3 -
p_q - p_{\bar{q}})^\mu \eta^{\nu \sigma} - (q_1 + q_3)^\sigma \eta^{\mu \nu},
\end{align}
\begin{align}
\label{eq:X6}
X_{6\;\mu \nu} &= \frac{\bar{u}_q \gamma_\mu (\slashed{q}_1 -
\slashed{p}_q)\gamma_\nu v_{\bar{q}}}{(q_1-p_q)^2}, \\
\label{eq:X7}
X_{7\;\mu\nu } &= \frac{\bar{u}_q \gamma_\nu (\slashed{q}_1 - \slashed{p}_{\bar{q}})\gamma_\mu v_{\bar{q}}}{(q_1-p_{\bar{q}})^2} \,,
\end{align}
and $q_3 = q_1 - p_q - p_{\bar{q}}$. In these expressions $X_{sym}$ corresponds
to the sum of the first 4 diagrams in Figure~\ref{fig:qq_qQQq_graphs} after
taking the QMRK limit, while $X_5$, $X_6$ and $X_7$ correspond to the 5th, 6th
and 7th diagrams respectively (see~\cite{Andersen:2020yax} for more details).
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen1}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen2}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen3}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen4}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen5}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen6}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{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}
\subsection{Higgs}
% They individual pieces are listed in section~\ref{sec:currents_impl} and the contraction is done in the current generator, section~\ref{sec:cur_gen}.
\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{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$.
This vertex is used in the current contractions for central and extremal Higgs boson emission:
\begin{align}
\label{eq:S_ff_fHf}
\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\,,\\
\label{eq:S_gf_Hf}
\left\|S_{g f_2 \to H f_2}\right\|^2 ={}& \sum_{
\substack{
\lambda_{g} = +,-\\
\lambda_{b} = +,-
}}
\left|\epsilon_\mu^{\lambda_g}(p_g, p_r)\ V_H^{\mu\nu}(p_g, p_g-p_H)\ j_\nu^{\lambda_b}(p_n, p_b)\right|^2\,,
\end{align}
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/Hjets.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}
\subsection{Vector Boson + Jets}
\label{sec:currents_W}
\subsubsection{Quark+ Vector Boson}
\begin{figure}
\centering
\begin{minipage}[b]{0.2\textwidth}
\includegraphics[width=\textwidth]{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]{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]{Wbits3.pdf}
\end{minipage}
\caption{The $j_V$ current is constructed from the two diagrams which
contribute to the emission of a vector boson from a given quark line.}
\label{fig:jV}
\end{figure}
For a $W, Z$, or photon emission we require a fermion. The current is actually a sum of
two separate contributions, see figure~\ref{fig:jV}, one with a vector boson emission
from the initial state, and one with the vector boson emission from the final state.
This can be seen as the following two
terms, given for the example of a $W$ emission~\cite{Andersen:2012gk,Andersen:2020yax}
\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.
The coupling and propagator factor in eq.~(\ref{eq:Weffcur1}) have to
be adapted depending on the emitted boson. The remaining product of
currents
\begin{equation}
\label{eq:J_V}
J_{\text{V}}^\mu(p_2,p_l,p_{\bar{l}},p_3)=\left( \frac{ \bar{u}_2 \gamma^\nu (\slashed{p}_2 +
\slashed{p}_l +
\slashed{p}_{\bar{l}}) \gamma^\mu u_3}{s_{2l\bar{l}}} - \frac{\bar u_2
\gamma^\mu(\slashed{p}_3 + \slashed{p}_l + \slashed{p}_{\bar{l}}) \gamma^\nu
u_3}{s_{3l\bar{l}}} \right) [\bar{u}_l \gamma_\nu u_{\bar{l}}]
\end{equation}
with $s_{il\bar{l}} = (p_i + p_l +p_{\bar{l}})^2$ is universal. The
implementation is in \texttt{include/currents.frm} inside the
\texttt{current\_generator} (see section~\ref{sec:cur_gen}). To use it
inside \FORM use the place-holder
\lstinline!JV(h1, hl, mu, pa, p1, plbar, pl)!, where \lstinline!h1! is
the helicity of the quark line and \lstinline!hl! the helicity of the
lepton line.
-For a final-state photon the propagator denominators and the leptonic
-current are replaced by the photon polarisation vector, viz.
+For a final-state photon the leptonic current is replaced by the
+photon polarisation vector, viz.
\begin{equation}
\label{eq:J_photon}
- j_\gamma^\mu(p_a,p_\gamma,p_1) = \bar{u}(p_1)\left[ \slashed{\varepsilon}(p_\gamma) (\slashed{p}_1 +
- p_\gamma) \gamma^\mu +
- \gamma^\mu(\slashed{p}_a - p_\gamma) \slashed{\varepsilon}(p_\gamma)
- \right]u(p_a)
+ j_\gamma^\mu(p_a,p_\gamma,p_1) = \bar{u}(p_1)\left( \frac{\slashed{\varepsilon}(p_\gamma) (\slashed{p}_1 +
+ p_\gamma) \gamma^\mu}{(p_1 + p_\gamma)^2} +
+ \frac{\gamma^\mu(\slashed{p}_a - p_\gamma) \slashed{\varepsilon}(p_\gamma)}{(p_a - p_\gamma)^2}
+ \right)u(p_a).
\end{equation}
\subsubsection{Vector boson with unordered emission}
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{vuno1}
\caption{}
\label{fig:U1diags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{vuno2}
\caption{}
\label{fig:U2diags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{vuno3}
\caption{}
\label{fig:Cdiags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{vuno4}
\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 vector
boson are emitted from the same quark line and the gluon comes after the
$t$-channel propagator. In (b), the gluon and vector 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:Vunodiags}
\end{figure}
It is necessary to include subleading processes in vector boson + jets
also. Similarly to the pure jet case, the unordered currents are not
calculated separately, and only in the ME functions when required in
the \texttt{src/Wjets.cc} or \texttt{src/Zjets.cc} file. The following shows the derivation of the calculation of
this ME within HEJ. We start with a contraction of two currents:
\begin{equation}
\label{eq:SabsVuno}
S_{qQ\to Vgq^\prime Q} =
j_{V{\rm uno}\,\mu}^d(p_a,p_1,p_2,p_\ell,p_{\bar\ell})\ g^{\mu
\nu}\ T^d_{3b}\ j^{h_b,h_3}_{\nu}(p_b,p_{3}),
\end{equation}
where $j_{V,{\rm uno}}$ is our new unordered current which is 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:Vunodiags}). Using
$T^m_{ij}$ represent fundamental colour matrices between quark state
$i$ and $j$ with adjoint index $m$ we find
\begin{align}\label{eq:vunocurrent}
\begin{split}
j^{d\,\mu}_{\rm V,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}_V)\rho |a\rangle }{s_{12}t_{aV}} + \frac{\langle 2|\nu (\slashed{p}_2+ \slashed{p}_1)\rho (\slashed{p}_2+\slashed{p}_1 + \slashed{p}_V)\mu |a\rangle }{s_{12}s_{12V}} \\
&+ \frac{\langle 2|\rho (\slashed{p}_2+ \slashed{p}_V) \nu
(\slashed{p}_1 + \slashed{p}_2+\slashed{p}_V)\mu |a\rangle}{s_{2V}s_{12V}}\,,
\end{split}\\
\label{eq:U2tensor}
\begin{split}
\tilde U_2^{\nu\mu\rho} ={}&\frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_V-\slashed{p}_1)\nu (\slashed{p}_a - \slashed{p}_V)\rho |a\rangle }{t_{aV1}t_{aV}} + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_V- \slashed{p}_1)\rho (\slashed{p}_a-\slashed{p}_1) \nu |a\rangle }{t_{a1V}t_{a1}} \\
&+ \frac{\langle 2|\rho (\slashed{p}_2+ \slashed{p}_V) \mu
(\slashed{p}_a-\slashed{p}_1)\nu |a\rangle}{s_{2V}t_{a1}}\,,
\end{split}\\
\label{eq:Ltensor}
\begin{split}
\tilde L^{\nu\mu\rho} ={}& \frac{1}{t_{aV2}}\left[
\frac{\langle 2 | \sigma (\slashed{p}_a-\slashed{p}_V)\rho|a\rangle}{t_{aV}}
+\frac{\langle 2 | \rho (\slashed{p}_2+\slashed{p}_V)\sigma|a\rangle}{s_{2V}}
\right]\\
&\times \left\{\left(\frac{p_b^\nu}{s_{1b}} + \frac{p_3^\nu}{s_{13}}\right)(q_1-p_1)^2g^{\mu\sigma}+(2q_1-p_1)^\nu g^{\mu\sigma} - 2p_1^\mu g^{\nu\sigma} + (2p_1-q_1)^\sigma g^{\mu\nu} \right\}\,,
\end{split}
\end{align}
where $s_{ij\dots} = (p_i + p_j + \dots)^2, t_{ij\dots} = (p_i - p_j -
\dots)^2$ and $q_1 = p_a-p_2-p_V$. This is actually calculated and
used in the code in a much cleaner way as follows:
\begin{align}\label{eq:spinorstringVuno}
S_{qQ\to Vgq^\prime Q} = i\varepsilon_\nu (p_g) \bar{u}(p_\ell) & \gamma_\rho v(p_{\bar \ell})\times T^d_{3b} \bar{u}^{h_3}(p_3)\gamma_\mu u^{h_b}(p_b) \times \nonumber \\
&\left( T^1_{2i}T^d_{ia} \left( \tilde{U}_1^{\nu\mu\rho}-\tilde{L}^{\nu\mu\rho}\right)+T^d_{2i}T^1_{ia}\left(\tilde{U}_2^{\nu\mu\rho}+\tilde{L}^{\nu\mu\rho}\right) \right)
\end{align}
If we define the objects:
\begin{align}\label{eq:VunoX}
X &= \varepsilon_\nu(p_g) \left[ \bar{u}(p_\ell) \gamma_\rho v(p_{\bar \ell}) \right] \left[ \bar{u}^{h_3}(p_3)\gamma_\mu u^{h_b}(p_b)\right] \left( \tilde{U}_1^{\nu\mu\rho}-\tilde{L}^{\nu\mu\rho}\right)\\
Y &= \varepsilon_\nu(p_g) \left[ \bar{u}(p_\ell) \gamma_\rho v(p_{\bar \ell}) \right] \left[ \bar{u}^{h_3}(p_3)\gamma_\mu u^{h_b}(p_b)\right] \left( \tilde{U}_2^{\nu\mu\rho}+\tilde{L}^{\nu\mu\rho}\right)
\label{eq:VunoY}
\end{align}
then we can rewrite Equation \eqref{eq:spinorstringVuno} in the much simpler form:
\begin{equation}
S_{qQ\to Vgq^\prime Q} = iT^d_{3b} \left( T^{1}_{2i}T^d_{ia} X + T^d_{2i}T^1_{ia} Y \right)
\end{equation}
then, by using:
\begin{align}
\sum_{\text{all indices}}& T^d_{3b}T^e_{b3}T^1_{2i}T^d_{ia}T^e_{ai}T^1_{i2} = \frac{1}{2}C_F^2C_A \\
\sum_{\text{all indices}}& T^d_{3b}T^e_{b3}T^1_{2i}T^d_{ia}T^1_{ai}T^e_{i2} = \frac{1}{2}C_F^2C_A - \frac{1}{4}C_A^2C_F = -\frac{1}{4}C_F
\end{align}
giving then, the spin summed and helicity averaged spinor string as:
\begin{equation}\label{eq:VunoSumAveS}
||\;\bar{S}_{qQ\to Vgq^\prime Q}\;|| = \frac{1}{4N_C^2} \left( \frac{1}{2}C_F^2C_A\left(|X|^2+|Y|^2\right)-\frac{1}{4}C_F\times2\mathrm{Re}\left(XY^*\right)\right)
\end{equation}
In the case of $Z$ boson emission and if the two incoming particles are (anti-)quarks one has to take into account the possibility that the boson can be emitted from either the backward or the forward leg. The current contractions needed in the second, third and fourth lines in eq.~(\ref{eq:ME_Z}) are:
\begin{align}
||\;\bar{S}_{\rm top}\;|| &= \frac{1}{4N_C^2} \left( \frac{1}{2}C_F^2C_A\left(|X_{\rm top}|^2+|Y_{\rm top}|^2\right)-\frac{1}{4}C_F\times2\mathrm{Re}\left(X_{\rm top} Y_{\rm top}^*\right)\right) \label{eq:Z_uno_top} \\
||\;\bar{S}_{\rm bot}\;|| &= \frac{1}{4N_C^2} \left( \frac{1}{2}C_F^2C_A\left(|X_{\rm bot}|^2+|Y_{\rm bot}|^2\right)-\frac{1}{4}C_F\times2\mathrm{Re}\left(X_{\rm bot} Y_{\rm bot}^*\right)\right) \label{eq:Z_uno_bot} \\
||\;\bar{S}_{\rm int.}\;|| &= \frac{1}{4N_C^2} \left( \frac{1}{2}C_F^2C_A\left(X_{\rm top} X_{\rm bot}^* + Y_{\rm top} Y_{\rm bot}^*\right)-\frac{1}{4}C_F\left(X_{\rm top} Y_{\rm bot}^* + Y_{\rm top} X_{\rm bot}^*\right)\right) \label{eq:Z_uno_int}
\end{align}
where $X_{\rm top}, Y_{\rm top}$ correspond to the case where the boson is emitted from the unordered current~(\ref{eq:VunoX}, \ref{eq:VunoY}) while $X_{\rm bot}, Y_{\rm bot}$ correspond to boson emission from the quark leg opposite to the unordered current and are obtained by contracting an unordered current~(\ref{eq:juno}) with a boson emission current (see eq.(\ref{eq:Weffcur1}) with $W \to Z$).
\subsubsection{\texorpdfstring{$V$}{V}+Extremal \texorpdfstring{$\mathbf{q\bar{q}}$}{qqbar}}
The $V$+Jet sub-leading processes containing an extremal $q\bar{q}$ are
related by crossing symmetry to the $V$+Jet unordered processes. This
means that one can simply perform a crossing symmetry argument on
eq.~(\ref{eq:vunocurrent}) 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 $V$-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 V q \bar{q}$ by evaluating
the current for $V$ 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 V,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
\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}_V)\rho |1\rangle }{t_{a2}s_{1V}} + \frac{\langle 2|\nu (\slashed{p}_a- \slashed{p}_2)\rho (\slashed{p}_a-\slashed{p}_2 - \slashed{p}_V)\mu |1\rangle }{t_{a2}t_{a2V}} \\
&- \frac{\langle 2|\rho (\slashed{p}_2+ \slashed{p}_V) \nu
(\slashed{p}_a - \slashed{p}_2-\slashed{p}_V)\mu
|1\rangle}{s_{2V}t_{a2V}}\,,
\end{split}\\
\label{eq:U2tensorX}
\begin{split}
\tilde U_{2,X}^{\nu\mu\rho} =&-\frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_V-\slashed{p}_1)\nu (\slashed{p}_1 + \slashed{p}_V)\rho |1\rangle }{t_{aV1}s_{1V}} + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_V- \slashed{p}_1)\rho (\slashed{p}_a-\slashed{p}_1) \nu |1\rangle }{t_{a1V}t_{a1}} \\
&+ \frac{\langle 2|\rho (\slashed{p}_2+ \slashed{p}_V) \mu
(\slashed{p}_a-\slashed{p}_1)\nu |1\rangle}{s_{2V}t_{a1}}\,,
\end{split}\\
\label{eq:LtensorX}
\begin{split}
\tilde L^{\nu\mu\rho}_X &= \frac{1}{s_{V12}}\left[-\frac{\langle 2 |\sigma (\slashed{p}_1 + \slashed{p}_V) \rho | 1\rangle}{s_{1V}} + \frac{\langle 2 |\rho (\slashed{p}_2 + \slashed{p}_V) \sigma | 1\rangle }{s_{2V}} \right] \\
&\vphantom{+\frac{1}{t_{aV2}}}\quad\cdot \left( -\left(
\frac{p_b^\nu}{s_{ab}} + \frac{p_n^\nu}{s_{an}} \right) (q_1+p_a)^2 g^{\sigma\mu}+ g^{\sigma \mu} (2q_1 +p_a)^\nu - g^{\mu \nu}(2p_a+q_1)^\sigma+ 2g^{\nu \sigma}p_a^\mu \right)\,,
\end{split}
\end{align}
where $q_1=-(p_1+p_2+p_V)$. Notice in particular the similarity to the $V$+uno scenario (from which
this has been derived).
\subsubsection{Central \texorpdfstring{$\mathbf{q\bar{q}}$}{qqbar} Vertex}
\label{sec:central_qqbar_V}
The final subleading process in the $V$+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 \FORM inside
\texttt{current\_generator/j\_Vqqbar\_j.frm}, which is detailed in
section~\ref{sec:contr_calc}.
If the boson is a $W$, it can be emitted from the central effective vertex (scenario
dealt with by the function \lstinline!ME_WCenqqbar_qq! in the file
\texttt{src/Wjets.cc}); or from either of the external quark legs
(scenario dealt with by \lstinline!ME_W_Cenqqbar_qq! in same file).
In the latter case we need to contract the pure jets central $q\bar{q}$
vertex with a quark current on one side and a boson emission current
on the other side. In practice this means that we can use the expressions
in section~\ref{sec:central_qqbar}, replacing $j_{\nu}(p_b,p_n)$ in
eqs.~(\ref{eq:Vcen}), (\ref{eq:Wcen}) with a boson current~(\ref{eq:Weffcur1}).
This contraction is implemented in \texttt{current\_generator/j\_qqbar\_jV.frm}.
To deal with the case of $V$ emission from the central $q\bar{q}$ pair we need
to compute the corresponding effective vertex. Since the boson emission doesn't
affect colour factors, the calculation can be organised in a similar way as in
the pure jets case, see section~\ref{sec:central_qqbar}. Each pure QCD diagram
in Figure~\ref{fig:qq_qQQq_graphs} will give rise to a subset of new diagrams
corresponding to all the possible ways to emit the boson. One can write the
skeleton function for this process as~\cite{Andersen:2020yax}
\footnote{See footnote~\ref{footnote_qqbar} for more details on the colour factors.}
\begin{equation}
\label{eq:Scen_V}
\left\|S_{f_a f_b\to f_a V q\bar{q}' f_b}^{\rm cen}\right\|^2 =
\frac{1}{2} \frac{1}{q_1^2q_3^2}\sum_{\{h_i\}}\ \left[\frac{C_A C_F^2}{4}\left(|\widetilde{V}|^2 + |\widetilde{W}|^2\right) +
\frac{C_F}{4} {\rm Re}\left(\widetilde{V} \widetilde{W}^*\right) \right],
\end{equation}
where
\begin{align}
\label{eq:Vcen_V}
\widetilde{V} &= j_{\mu}(p_a,p_1) \left( \tilde{X}_s^{\mu\nu} + \tilde{X}^{\mu \nu}_{6} \right)
j_{\nu}(p_b,p_n) \\
\label{eq:Wcen_V}
\widetilde{W} &= j_{\mu}(p_a,p_1) \left( \tilde{X}_s^{\mu\nu} + \tilde{X}^{\mu
\nu}_{7} \right) j_{\nu}(p_b,p_n)
\end{align}
with
\begin{equation}
\label{eq:Xs_V}
\tilde X_s^{\mu\nu} = \frac{1}{s_{q\bar{q}l\bar{l}}}\left(\tilde{X}_{\rm
sym}^{\sigma} \eta^{\mu\nu} + \tilde X_5^{\sigma\mu\nu}\right)J_{\text{V}
\sigma}(p_q,p_l,p_{\bar{l}},p_{\bar{q}}),
\end{equation}
\begin{align}
\label{eq:Xsym_V}
\begin{split}
\tilde{X}_{\rm sym}^{\sigma} =& q_1^2 \left(
\frac{p^\sigma_a}{s_{aq}+s_{a\bar{q}}+s_{a\ell}+s_{a\bar\ell}} +
\frac{p^\sigma_1}{s_{1q}+s_{1\bar{q}}+s_{1\ell}+s_{1\bar\ell}} \right) \\
& - q_3^2 \left(
\frac{p^\sigma_b}{s_{bq}+s_{b\bar{q}}+s_{b\ell}+s_{b\bar\ell}} +
\frac{p^\sigma_n}{s_{nq}+s_{n\bar{q}}+s_{n\ell}+s_{n\bar\ell}}\right),
\end{split}
\end{align}
\begin{equation}
\label{eq:X5_V}
\tilde X_5^{\sigma\mu\nu} = \left(q_1+p_q+p_{\bar{q}}+p_l+p_{\bar{l}}\right)^\nu\eta^{\mu\sigma}
+\left(q_3-p_q-p_{\bar{q}}-p_l-p_{\bar{l}}\right)^\mu\eta^{\sigma\nu}
-\left(q_1+q_3\right)^\sigma\eta^{\mu\nu},
\end{equation}
\begin{equation}
\label{eq:JVcen}
J_V^\sigma(p_q,p_\ell,p_{\bar{\ell}},p_{\bar{q}})=K_V \
[\bar{u}_\ell \gamma_\rho u_{\bar{\ell}}] \left( \frac{ [\bar{u}_q \gamma^\rho (\slashed{p}_q +
\slashed{p}_l +
\slashed{p}_{\bar{l}}) \gamma^\sigma u_{\bar{q}}]}{s_{q\ell\bar{\ell}}} -
\frac{[\bar u_q
\gamma^\sigma(\slashed{p}_{\bar{q}} + \slashed{p}_l + \slashed{p}_{\bar{l}}) \gamma^\rho
u_{\bar{q}}]}{s_{\bar{q}\ell\bar{\ell}}} \right),
\end{equation}
\begin{align}
\label{eq:X6_V}
\tilde X^{\mu\nu}_{6}=-K_V[\bar{u}_\ell \gamma_\sigma v_{\bar{\ell}}] \times \bar{u}_q \Biggl[-\frac{
\gamma^\sigma(\slashed{p}_q+\slashed{p}_l+\slashed{p}_{\bar{l}})\gamma^\mu
(\slashed{q}_3+ \slashed{p}_{\bar{q}})\gamma^\nu}{(p_q + p_l +
p_{\bar{l}})^2(q_3 +p_{\bar{q}})^2} & +\frac{\gamma^\mu(\slashed{q}_1-\slashed{p}_q)\gamma^\sigma(\slashed{q}_3+\slashed{p}_{\bar{q}})\gamma^\nu}{(q_1 - p_q)^2(q_3+p_{\bar{q}})^2} \nonumber \\
& +\frac{\gamma^\mu(\slashed{q}_1-\slashed{p}_q)\gamma^\nu(\slashed{p}_{\bar{q}}+\slashed{p}_l+\slashed{p}_{\bar{l}})\gamma^\sigma
}{(q_1 - p_q)^2(p_{\bar{q}} + p_l + p_{\bar{l}})^2}\Biggr]v_{\bar{q}},
\end{align}
and
\begin{align}
\label{eq:X7_V}
\tilde X^{\mu\nu}_{7}=K_V[\bar{u}_\ell
\gamma_\sigma v_{\bar{\ell}}] \times \bar{u}_q \Biggl[-\frac{
\gamma^\nu(\slashed{q}_3+\slashed{p}_q)\gamma^\mu
(\slashed{p}_{\bar{q}}+\slashed{p}_l+\slashed{p}_{\bar{l}})\gamma^\sigma}{(q_3
+p_q)^2(p_{\bar{q}}+p_l+p_{\bar{l}})^2} & +\frac{\gamma^\nu (\slashed{q}_3+\slashed{p}_q) \gamma^\sigma
(\slashed{q}_1-\slashed{p}_{\bar{q}}) \gamma^\mu}{(q_3 + p_q)^2(q_1 - p_{\bar{q}})^2} \nonumber \\
& + \frac{\gamma^\sigma(\slashed{p}_q+\slashed{p}_l+\slashed{p}_{\bar{l}})\gamma^\nu(\slashed{q}_1-\slashed{p}_{\bar{q}})\gamma^\mu
}{(p_q+p_l+p_{\bar{l}})^2(q_1 - p_{\bar{q}})^2}\Biggr]v_{\bar{q}}.
\end{align}
If for example the emitted boson is a $W$ we have $K_V=K_W$ with
\begin{equation}
\label{eq:KW}
K_W= \left(\frac{ig_W}{\sqrt{2}}\right)^2\frac{1}{\left(p_W^2 - M_W^2 + i \Gamma_W M_W\right)}.
\end{equation}
The numbering of the terms follows the same logic as in section~\ref{sec:central_qqbar},
e.g. $\tilde X_5$ corresponds to the subset of diagrams obtained by adding
a $V$ emission where allowed in diagram 5 of Figure~\ref{fig:qq_qQQq_graphs}.
If the emitted boson is a $Z$ the number of contributions at the amplitude level
depends on the incoming state. In the case of two incoming gluons the boson can
only be emitted from the $q\bar{q}$ pair and thus the skeleton function has the
same form as~(\ref{eq:Scen_V}) (taking into account the fact that the $Z/\gamma^*$
couples to all helicities). This is dealt with by the function
\lstinline!ME_ZCenqqbar_gg! in \texttt{src/Zjets.cc}. On the other hand, if the
two incoming particles are quarks, we have to take into account emission from
the upper quark leg, the $q\bar{q}$ pair and the lower quark leg. The 6 current
contractions entering the generalized version of eq.~(\ref{eq:ME_Z}) to the case
$n=3$ are:
\begin{align}
\left\|S_{\rm top}\right\|^2 &=
\frac{1}{2} \frac{1}{q_{1, {\rm top}}^2q_{3, {\rm top}}^2}\sum_{\{h_i\}}\ \left[\frac{C_A C_F^2}{4}\left(|\widetilde{V}_{\rm top}|^2 + |\widetilde{W}_{\rm top}|^2\right) +
\frac{C_F}{4} {\rm Re}\left(\widetilde{V}_{\rm top} \widetilde{W}_{\rm top}^*\right) \right] \label{eq:Z_cen_top} \\
\left\|S_{\rm mid}\right\|^2 &=
\frac{1}{2} \frac{1}{q_{1, {\rm mid}}^2q_{3, {\rm mid}}^2}\sum_{\{h_i\}}\ \left[\frac{C_A C_F^2}{4}\left(|\widetilde{V}_{\rm mid}|^2 + |\widetilde{W}_{\rm mid}|^2\right) +
\frac{C_F}{4} {\rm Re}\left(\widetilde{V}_{\rm mid} \widetilde{W}_{\rm mid}^*\right) \right] \label{eq:Z_cen_mid} \\
\left\|S_{\rm bot}\right\|^2 &=
\frac{1}{2} \frac{1}{q_{1, {\rm bot}}^2q_{3, {\rm bot}}^2}\sum_{\{h_i\}}\ \left[\frac{C_A C_F^2}{4}\left(|\widetilde{V}_{\rm bot}|^2 + |\widetilde{W}_{\rm bot}|^2\right) +
\frac{C_F}{4} {\rm Re}\left(\widetilde{V}_{\rm bot} \widetilde{W}_{\rm bot}^*\right) \right] \label{eq:Z_cen_bot} \\
\left\|S_{\rm top, mid}\right\|^2 &=
\frac{1}{2} \frac{1}{|q_{1, {\rm top}}||q_{3, {\rm top}}||q_{1, {\rm mid}}||q_{3, {\rm mid}}|}\sum_{\{h_i\}}\ \left[\frac{C_A C_F^2}{4}\left(\widetilde{V}_{\rm top}\widetilde{V}_{\rm mid}^* + \widetilde{W}_{\rm top}\widetilde{W}_{\rm mid}^*\right) \right. \nonumber \\
& \hspace{6cm} \left. + \frac{C_F}{8} \left(\widetilde{V}_{\rm top} \widetilde{W}_{\rm mid}^* + \widetilde{W}_{\rm top} \widetilde{V}_{\rm mid}^*\right) \right] \label{eq:Z_cen_top_mid} \\
\left\|S_{\rm top, bot}\right\|^2 &=
\frac{1}{2} \frac{1}{|q_{1, {\rm top}}||q_{3, {\rm top}}||q_{1, {\rm bot}}||q_{3, {\rm bot}}|}\sum_{\{h_i\}}\ \left[\frac{C_A C_F^2}{4}\left(\widetilde{V}_{\rm top}\widetilde{V}_{\rm bot}^* + \widetilde{W}_{\rm top}\widetilde{W}_{\rm bot}^*\right) \right. \nonumber \\
& \hspace{6cm} \left. + \frac{C_F}{8} \left(\widetilde{V}_{\rm top} \widetilde{W}_{\rm bot}^* + \widetilde{W}_{\rm top} \widetilde{V}_{\rm bot}^*\right) \right] \label{eq:Z_cen_top_bot} \\
\left\|S_{\rm mid, bot}\right\|^2 &=
\frac{1}{2} \frac{1}{|q_{1, {\rm mid}}||q_{3, {\rm mid}}||q_{1, {\rm bot}}||q_{3, {\rm bot}}|}\sum_{\{h_i\}}\ \left[\frac{C_A C_F^2}{4}\left(\widetilde{V}_{\rm mid}\widetilde{V}_{\rm bot}^* + \widetilde{W}_{\rm mid}\widetilde{W}_{\rm bot}^*\right) \right. \nonumber \\
& \hspace{6cm} \left. + \frac{C_F}{8} \left(\widetilde{V}_{\rm mid} \widetilde{W}_{\rm bot}^* + \widetilde{W}_{\rm mid} \widetilde{V}_{\rm bot}^*\right) \right] \label{eq:Z_cen_mid_bot}
\end{align}
In these expressions the subscripts ``top'', ``mid'' and ``bot'' are used to
indicate that the boson is emitted from the upper quark leg, the central $q\bar{q}$
pair or the lower quark leg respectively. One can note that
$q_{1, {\rm bot}}=q_{1, {\rm mid}}$ and $q_{3, {\rm mid}}=q_{3, {\rm top}}$.
The function \lstinline!ME_ZCenqqbar_qq! in \texttt{src/Zjets.cc} returns the
6 values obtained by evaluating eqs.~(\ref{eq:Z_cen_top})-(\ref{eq:Z_cen_mid_bot}),
arranged in the same order. If one of the incoming partons is a gluon the boson
cannot be emitted from this leg so we can drop all the terms involving the
corresponding subscript in the previous expressions. We are thus left with two
squares and one interference term, as with other processes having two contributions
at the amplitude level. This is dealt with by \lstinline!ME_ZCenqqbar_gq! in
\texttt{src/Zjets.cc}.
\subsubsection{$W^+ W^+$ + jets}
The production of same-sign $WW$ + jets is implemented through the contraction
of two vector boson currents; see eqn~\ref{eq:J_V} and more generally
subsection~\ref{sec:currents_W}. This contraction is available in FORM under the
header \texttt{jV\_jV.hh}, while the matrix elements can be found in
\texttt{src/WWjets.cc}.
There are two distinct contributions to this process which correspond to the
possible pairings of the $W$-bosons and incoming legs. Labelling the bosons as
\texttt{W1} and \texttt{W2}, the two configurations correspond to
\texttt{W1}-forward with \texttt{W2}-backward, and \texttt{W2}-forward with
\texttt{W1}-backward. The interference between these contributions is included.
For same-flavour decays, \texttt{reconstruct\_intermediate()} will compute the
difference between reconstruced mass and reference mass in possible pairings and
select the pairing which minimises this quantity.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "developer_manual"
%%% End: