Page MenuHomeHEPForge
Files F9501184

No OneTemporary
Actions

View Options

diff --git a/current_generator/W_central_qqbar.frm b/current_generator/W_central_qqbar.frm
index 19fe09b..18b588b 100644
--- a/current_generator/W_central_qqbar.frm
+++ b/current_generator/W_central_qqbar.frm
@@ -1,116 +1,116 @@
*/**
-* \brief Calculation of W + jets with central qqbar emission, W emitted from extremal leg
+* \brief Calculation of W + jets with central qqbar emission, W emitted from central leg
*
* TODO: unify conventions with developer manual
* the current dictionary is as follows:
*
* code | manual
* pq | p_2
* pqbar | p_3
* pl | p_A
* plbar | p_B
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include- include/helspin.frm
#include- include/write.frm
i mu1,...,mu20;
i mu,nu,rho,sigma;
cf T1a, T4b, Xunc, Xcro, X1a, X4b, V3g, JV;
v p1,pa,p4,pb,pq,pqbar,pl,plbar,q1,q3,q11,q12;
s h;
s tunc1,tunc2,tcro1,tcro2,s2AB,s3AB,s23AB;
s t1,t3,s12,s13,s1A,s1B,sa2,sa3,saA,saB,s42,s43,s4A,s4B,sb2,sb3,sbA,sbB;
cf m2;
#do h1={+,-}
#do h2={+,-}
l [M_uncross_W `h1'`h2'] = T1a(`h1'1, mu)*T4b(`h2'1, nu)*Xunc(mu, nu);
l [M_cross_W `h1'`h2'] = T1a(`h1'1, mu)*T4b(`h2'1, nu)*Xcro(mu, nu);
l [M_sym_W `h1'`h2'] = T1a(`h1'1, mu)*T4b(`h2'1, nu)*(
X1a(mu, nu, sigma) - X4b(mu, nu, sigma) + V3g(mu, nu, sigma)
)/s23AB * JV(sigma);
#enddo
#enddo
id T1a(h?, mu?) = Current(h, p1, mu, pa);
id T4b(h?, mu?) = Current(h, p4, mu, pb);
* see eq:X_Unc in developer manual, multiplied by -t1*t3
id Xunc(mu, nu) = -Current(-1, pl, sigma, plbar)*(
- Current(-1, pq, sigma, pq + pl + plbar, mu, q3 + pqbar, nu, pqbar)/s2AB/tunc2
+ Current(-1, pq, mu, q1 - pq, sigma, q3 + pqbar, nu, pqbar)/tunc1/tunc2
+ Current(-1, pq, mu, q1 - pq, nu, pqbar + pl + plbar, sigma, pqbar)/tunc1/s3AB
);
* see eq:X_Cro in developer manual, multiplied by -t1*t3
id Xcro(mu, nu) = -Current(-1, pl, sigma, plbar)*(
- Current(-1, pq, nu, q3 + pq, mu, pqbar + pl + plbar, sigma, pqbar)/tcro1/s3AB
+ Current(-1, pq, nu, q3 + pq, sigma, q1 - pqbar, mu, pqbar)/tcro1/tcro2
+ Current(-1, pq, sigma, pq + pl + plbar, nu, q1 - pqbar, mu, pqbar)/s2AB/tcro2
);
* see eq:X_1a in developer manual, multiplied by -t1*t3
id X1a(mu, nu, sigma) = t1*d_(mu, nu)*(
+ p1(sigma)/(s12 + s13 + s1A + s1B)
+ pa(sigma)/(sa2 + sa3 + saA + saB)
);
* see eq:X_4b in developer manual, multiplied by t1*t3
id X4b(mu, nu, sigma) = t3*d_(mu, nu)*(
+ p4(sigma)/(s42 + s43 + s4A + s4B)
+ pb(sigma)/(sb2 + sb3 + sbA + sbB)
);
* see eq:3GluonWEmit in developer manual, multiplied by t1*s23AB*t3
id V3g(mu, nu, sigma) = (
+ (q1(nu) + pq(nu) + pqbar(nu) + pl(nu) + plbar(nu))*d_(mu, sigma)
+ (q3(mu) - pq(mu) - pqbar(mu) - pl(mu) - plbar(mu))*d_(nu, sigma)
- (q1(sigma) + q3(sigma))*d_(mu, nu)
);
* see eq:J_V in developer manual
id JV(sigma?) = -i_*(
+ Current(-1, pq, rho, pq + pl + plbar, sigma, pqbar)/s2AB
- Current(-1, pq, sigma, pqbar + pl + plbar, rho, pqbar)/s3AB
)*Current(-1, pl, rho, plbar);
multiply replace_(
s2AB, m2(pl+plbar+pq),
s3AB, m2(pl+plbar+pqbar),
s23AB, m2(pl+plbar+pq+pqbar),
tunc1, m2(q1-pq),
tunc2, m2(q3+pqbar),
tcro1, m2(q3+pq),
tcro2, m2(q1-pqbar),
t1, m2(q1),
t3, m2(q3),
s12, m2(p1+pq),
s13, m2(p1+pqbar),
s1A, m2(p1+pl),
s1B, m2(p1+plbar),
sa2, m2(pa+pq),
sa3, m2(pa+pqbar),
saA, m2(pa+pl),
saB, m2(pa+plbar),
s42, m2(p4+pq),
s43, m2(p4+pqbar),
s4A, m2(p4+pl),
s4B, m2(p4+plbar),
sb2, m2(pb+pq),
sb3, m2(pb+pqbar),
sbA, m2(pb+pl),
sbB, m2(pb+plbar)
);
multiply replace_(q3, q1-pl-plbar-pq-pqbar);
* replace q1 by sum of lightlike momenta
* @TODO: choose q11 as collinear to some other vector
multiply replace_(q1, q11-q12);
.sort
#call ContractCurrents()
.sort
format O4;
format c;
#call WriteHeader(`OUTPUT')
#call WriteOptimised(`OUTPUT',M_uncross_W,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
#call WriteOptimised(`OUTPUT',M_cross_W,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
#call WriteOptimised(`OUTPUT',M_sym_W,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
#call WriteFooter(`OUTPUT')
.end
diff --git a/current_generator/W_extremal_qqbar.frm b/current_generator/W_extremal_qqbar.frm
index 856fb9f..7305403 100644
--- a/current_generator/W_extremal_qqbar.frm
+++ b/current_generator/W_extremal_qqbar.frm
@@ -1,102 +1,105 @@
*/**
* \brief Calculation of W + extremal qqbar current contraction
*
* TODO: unify conventions with developer manual
* the current dictionary is as follows:
*
* code | manual
* pq | p_2
* pqbar | p_1
* pg | p_a
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include- include/helspin.frm
#include- include/write.frm
s h,sgb,sgn,sqqbarW;
v pq,pqbar,pg,pl,plbar,pn,pb,pW,pr,q1,q2;
i mu,nu,rho,sigma;
i mu1,...,mu20;
cf epsg,m2,m2inv,sqrt;
#do h1={+,-}
* eq:U1tensorX in developer manual, up to factors 1/sij, 1/tij
l [U1X `h1'] = (
+ Current(`h1'1, pq, nu, pg-pq, mu, pqbar+pW, rho, pqbar)
+ Current(`h1'1, pq, nu, pg-pq, rho, pg-pq-pW, mu, pqbar)
- Current(`h1'1, pq, rho, pq+pW, nu, pg-pq-pW, mu, pqbar)
);
* eq:U2tensorX in developer manual, up to factors 1/sij, 1/tij
l [U2X `h1'] = (
- Current(`h1'1, pq, mu, pg-pqbar-pW, nu, pqbar+pW, rho, pqbar)
+ Current(`h1'1, pq, mu, pg-pqbar-pW, rho, pg-pqbar, nu, pqbar)
+ Current(`h1'1, pq, rho, pq+pW, mu, pg-pqbar, nu, pqbar)
);
* eq:LtensorX in developer manual, up to factors 1/sij, 1/tij
l [LX `h1'] = (
- Current(`h1'1, pq, sigma, pqbar+pW, rho, pqbar) +
+ Current(`h1'1, pq, rho, pq+pW, sigma, pqbar)
)*(
+ (-(pb(nu)/sgb + pn(nu)/sgn)*m2(q1 + pg) + 2*q1(nu) - pg(nu))*d_(mu, sigma)
- (2*pg(sigma) + q1(sigma))*d_(mu, nu)
+ 2*pg(mu)*d_(nu, sigma)
)/sqqbarW;
#enddo
.sort
* restore kinematic factors
id Current(h?, pq, mu?, q1?, nu?, q2?, rho?, pqbar) = (
Current(h, pq, mu, q1, nu, q2, rho, pqbar)*m2inv(q1)*m2inv(q2)
);
id Current(h?, pq, mu?, q1?, nu?, pqbar) = (
Current(h, pq, mu, q1, nu, pqbar)*m2inv(q1)
);
.sort
drop;
* multiply with polarisation vector and other currents
#do h1={+,-}
#do h2={+,-}
#do hg={+,-}
#do TENSOR={U1X,U2X,LX}
l [`TENSOR' `h1'`h2'`hg'] = (
[`TENSOR' `h1']
*epsg(`hg'1, nu)
*Current(`h2'1, pn, mu, pb)
*Current(-1, pl, rho, plbar)
);
#enddo
#enddo
#enddo
#enddo
-* TODO: choice of best reference vector (pn or pb)
-id epsg(h?, mu?) = epsg(h, mu, pn);
+* choice of reference vector (pn or pb)
+* this is not exactly the optimal choice for all helicities,
+* but it's pretty close and the logic is simple
+id epsg(h?, nu?)*Current(h?, pn, mu?, pb) = epsg(h, nu, pn)*Current(h, pn, mu, pb);
+also epsg(h?, nu?) = epsg(h, nu, pb);
id epsg(-1, mu?, pr?) = sqrt(2)/2*SpinorChain(pr, mu, pg)/AngleChain(pg,pr);
id epsg(+1, mu?, pr?) = sqrt(2)/2*SpinorChain(pg, mu, pr)/SquareChain(pg,pr);
multiply replace_(q1,-(pq+pqbar+pW));
multiply replace_(pW,pl+plbar);
.sort
#call ContractCurrents
multiply replace_(
sgn,m2(pg+pn),
sgb,m2(pg+pb),
sqqbarW,m2(pq+pqbar+pW)
);
id m2inv(q1?) = 1/m2(q1);
multiply replace_(pW,pl+plbar);
.sort
format O4;
format c;
#call WriteHeader(`OUTPUT')
#call WriteOptimised(`OUTPUT',U1X,3,pg,pq,plbar,pl,pqbar,pn,pb)
#call WriteOptimised(`OUTPUT',U2X,3,pg,pq,plbar,pl,pqbar,pn,pb)
#call WriteOptimised(`OUTPUT',LX,3,pg,pq,plbar,pl,pqbar,pn,pb)
#call WriteFooter(`OUTPUT')
.end
diff --git a/current_generator/extremal_W_central_qqbar.frm b/current_generator/extremal_W_central_qqbar.frm
new file mode 100644
index 0000000..a3898d9
--- /dev/null
+++ b/current_generator/extremal_W_central_qqbar.frm
@@ -0,0 +1,82 @@
+*/**
+* \brief Calculation of W + jets with central qqbar emission, W emitted from extremal leg
+*
+* TODO: document in developer manual!
+*
+* \authors The HEJ collaboration (see AUTHORS for details)
+* \date 2020
+* \copyright GPLv2 or later
+*/
+#include- include/helspin.frm
+#include- include/write.frm
+
+cf Xsym,V3g,jW;
+i alpha,mu,nu,sigma;
+s h,sqqbar,saq,saqbar,s1q,s1qbar,sbq,sbqbar,s4q,s4qbar;
+cf m2;
+v p,pa,p1,pb,p4,pq,pqbar,pl,plbar,q1,q3,q11,q12,pW;
+
+* cf. eq:Xcen in developer manual
+#do h1={+,-}
+ #do h2={+,-}
+* eq:Xcrossed in developer manual, normalisation added later on
+ l [M_cross `h1'`h2'] = -Current(`h1'1, p1, mu, pa)*Current(
+ `h2'1, pq, nu, q1 - pqbar, mu, pqbar
+ );
+* eq:Xqprop in developer manual, normalisation added later on
+ l [M_qbar `h1'`h2'] = Current(`h1'1, p1, mu, pa)*Current(
+ `h2'1, pq, mu, q1 - pq, nu, pqbar
+ );
+ l [M_sym `h1'`h2'] = -i_*Current(`h1'1, p1, mu, pa)*(
+ d_(mu,nu)*Xsym(sigma) + V3g(mu, nu, sigma)
+ )/sqqbar*Current(`h2'1, pq, sigma, pqbar);
+ #enddo
+#enddo
+multiply jW(-1, pb, p4, plbar, pl);
+* eq:Weffcur1 in developer manual
+id jW(-1, pa?, p1?, plbar?, pl?) = Current(-1, pl, alpha, plbar)*(
+ + Current(-1,p1,alpha,p1+pW,nu,pa)
+ + Current(-1,p1,nu,pa-pW,alpha,pa)
+);
+
+id Current(h?, pq?, mu?, p?, nu?, pqbar?) = Current(h, pq, mu, p, nu, pqbar)/m2(p);
+
+* eq:Xsym in developer manual
+id Xsym(sigma?) = m2(q1)*(
+ pa(sigma)/(saq + saqbar) + p1(sigma)/(s1q + s1qbar)
+) - m2(q3)*(
+ pb(sigma)/(sbq + sbqbar) + p4(sigma)/(s4q + s4qbar)
+);
+* eq:V3g in developer manual
+id V3g(mu?, nu?, sigma?) = (
+ (q1(nu) + pq(nu) + pqbar(nu))*d_(mu, sigma)
+ + (q3(mu) - pq(mu) - pqbar(mu))*d_(nu, sigma)
+ - (q1(sigma) + q3(sigma))*d_(mu, nu)
+);
+
+multiply replace_(pW, pl+plbar);
+multiply replace_(q3, q1-pq-pqbar);
+* replace q1 by sum of lightlike momenta
+multiply replace_(q1, q11-q12);
+.sort
+#call ContractCurrents()
+multiply replace_(
+ sqqbar, m2(pq+pqbar),
+ saq, m2(pa+pq),
+ saqbar, m2(pa+pqbar),
+ s1q, m2(p1+pq),
+ s1qbar, m2(p1+pqbar),
+ sbq, m2(pb+pq),
+ sbqbar, m2(pb+pqbar),
+ s4q, m2(p4+pq),
+ s4qbar, m2(p4+pqbar)
+);
+.sort
+format O4;
+format c;
+#call WriteHeader(`OUTPUT')
+#call WriteOptimised(`OUTPUT',M_qbar,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
+#call WriteOptimised(`OUTPUT',M_cross,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
+#call WriteOptimised(`OUTPUT',M_sym,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
+#call WriteFooter(`OUTPUT')
+.end
diff --git a/doc/developer_manual/currents.tex b/doc/developer_manual/currents.tex
index d6725ae..a5aa485 100644
--- a/doc/developer_manual/currents.tex
+++ b/doc/developer_manual/currents.tex
@@ -1,672 +1,693 @@
\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. All $W$+Jet currents
are located in \texttt{src/Wjets.cc}, all Higgs+Jets currents are
defined in \texttt{src/Hjets.cc}, and pure jet currents are defined in
in \texttt{src/jets.cc}. All of these have their own separate header
files: \texttt{include/HEJ/Wjets.hh}, \texttt{include/HEJ/Hjets.hh} and
\texttt{include/HEJ/jets.hh} respectively.
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{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!ME_unob_XX!).
\subsubsection{Extremal $q\bar{q}$}
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.
\subsubsection{Central $q\bar{q}$}
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}.
\\
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_2,y_3 \ll y_4.
\end{align}
The HEJ Matrix element for this process can be calculated as:
\begin{align}
\label{eq:Mcentral}
i\mathcal{M} &= g_s^4 T^d_{1a} T^e_{4b}\ \frac{j_{\mu}(p_a,p_1)\ X^{ab\, \mu
\nu}_{{\rm cen}}(p_2,p_3,q_1,q_3)\
j_{\nu}(p_b,p_4)}{t_{a1}t_{b4}}.
\end{align}
where $X^{\mu \nu}_{\rm cen}$ is given by:
-\begin{align}
- X^{\mu \nu}_{\rm cen} &=\frac{f^{ced}T^c_{23}}{s_{23}}
+\begin{equation}
+ \label{eq:Xcen}
+ \begin{split}
+ X^{\mu \nu}_{\rm cen} ={}&\frac{f^{ced}T^c_{23}}{s_{23}}
\left(\eta^{\mu \nu} X_{sym}^\sigma + V^{\mu \nu \sigma}_{3g} \right)
\bar{u}(p_2) \gamma^\sigma u(p_3) \\ & \qquad + \frac{i
- T^d_{2q}T^e_{q3}}{(q_1-p_2)^2} X^{\mu\nu}_{qprop} - \frac{i
- T^e_{2q}T^d_{q3}}{(q_1-p_3)^2} X^{\mu\nu}_{crossed},
+ T^d_{2q}T^e_{q3}}{(q_1-p_2)^2} X^{\mu\nu}_{\text{qprop}} - \frac{i
+ T^e_{2q}T^d_{q3}}{(q_1-p_3)^2} X^{\mu\nu}_{\text{crossed}}\,,
+ \end{split}
+\end{equation}
+with
+\begin{align}
+ \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_4^\sigma}{s_{4q} + s_{4\bar{q}}}
+ \right)\,,\\
+ \label{eq:V3g}
+ V_{3g}^{\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}\,,\\
+ \label{eq:Xqprop}
+ X^{\mu\nu}_{\text{qprop}} ={}& \frac{\langle p_q | \mu (q_1-p_q) \nu | p_{\bar{q}}\rangle}{(q_1-p_q)^2}\,,\\
+ \label{eq:Xcrossed}
+ X^{\mu\nu}_{\text{crossed}} ={}& \frac{\langle p_q | \nu (q_1-p_{\bar{q}}) \mu | p_{\bar{q}}\rangle}{(q_1-p_{\bar{q}})^2}\,,
\end{align}
+and $q_3 = q_1 - p_q - p_{\bar{q}}$.
\begin{figure}[t]
\centering
\includegraphics[]{figures/Cenqqbar_jx}
\caption{An example of a Central $q\bar{q}$ process}
\label{fig:qqbarcen_example}
\end{figure}
\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{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/Hjets.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/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}
\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/Hjets.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/Hjets.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/Hjets.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/Hjets.cc}, each function
calculates the square of eq. (4.23) and (4.22) from~\cite{DelDuca:2003ba} respectively.
\subsection{$W$+Jets}
\label{sec:currents_W}
\subsubsection{Quark+$W$}
\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_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~\cite{Andersen:2012gk}\todo{cite W subleading paper}:
\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 FKL $W$ current is
\todo{Text. The following is what's currently implemented in \lstinline!jW!}
\begin{align}
\label{eq:jW-}
j^-_{W}(p_a, p_1, p_{\bar{l}}, p_{l}) ={}&
- \frac{2 \spa 1.l}{(p_1+p_l+p_{\bar{l}})^2}\*\Big(
- \spa 1.{\bar{l}} j^-(p_1, p_a) + \spb l.{\bar{l}}
+ \frac{2 \spa 1.l}{(p_1+p_W)^2}\*\Big(
+ \spb 1.{\bar{l}} j^-(p_1, p_a) + \spb l.{\bar{l}}
j^-(p_{l}, p_a)\Big)\notag\\
& + \frac{2\spb a.{\bar{l}}}{(p_a-p_W)^2}\Big(
- \spa 1.l j^-(p_l, p_a) + \spa a.l
+ \spa l.{\bar{l}} j^-(p_1, p_{\bar{l}}) + \spa a.l
j^-(p_1, p_a)
\Big)\,,\\
\label{eq:jW+}
j^+_{W}(p_a, p_1, p_{\bar{l}}, p_{l}) =& \big[j^-_{W}(p_a, p_1, p_l, p_{\bar{l}})\big]^*\,,
\end{align}
where the negative-helicity current is used for emission off a quark
line and the positive-helicity current for emissions off antiquark.
\subsubsection{$W$+uno}
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{wuno1}
\caption{}
\label{fig:U1diags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{wuno2}
\caption{}
\label{fig:U2diags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{wuno3}
\caption{}
\label{fig:Cdiags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{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}\\
\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}\\
\label{eq:Ltensor}
\begin{split}
\tilde L^{\nu\mu\rho} ={}& \frac{1}{t_{aW2}}\left[
\frac{\langle 2 | \sigma (\slashed{p}_a-\slashed{p}_W)\rho|a\rangle}{t_{aW}}
+\frac{\langle 2 | \rho (\slashed{p}_2+\slashed{p}_W)\sigma|a\rangle}{s_{2W}}
\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_W$.
\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{1}{s_{W12}}\left[-\frac{\langle 2 |\sigma (\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) \sigma | 1\rangle }{s_{2W}} \right] \\
&\vphantom{+\frac{1}{t_{aW2}}}\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_W)$. 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 \lstinline!ME_WCenqqx_qq! in the file
\texttt{src/Wjets.cc}); or from either of the external quark legs
(scenario dealt with by \lstinline!ME_W_Cenqqx_qq! 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{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}
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 and
\todo{Different convention for $J_V$}
\begin{equation}
\label{eq:J_V}
J_V^\mu(p_2,p_l,p_{\bar{l}},p_3)=-i\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}
The following tensor structures correspond to groupings of diagrams in
figure~\ref{fig:qq_qQQq_graphs}.
\begin{eqnarray}
\label{eq:X_1a}
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:X_4b}
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}
V^{\mu\nu\sigma}_{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 \lstinline!M_sym_W! that is part of the header
\texttt{W\_central\_qqbar.hh} generated by the current generator
(section~\ref{sec:cur_gen}) returns a contraction of the above tensor
containing the information from these 5 groupings of contributions (30
diagrams in total).
\begin{align}
\label{eq:X_Unc}
\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 current
contraction can be obtained in the code with the function
\lstinline!M_uncross_W! in the autogenerated header \texttt{W\_central\_qqbar.hh}.
\begin{align}
\begin{split}
\label{eq:X_Cro}
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 current
contraction can be obtained in the code with the function
\lstinline!M_cross_W! in the generated header
\texttt{W\_central\_qqbar.hh}.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "developer_manual"
%%% End:
diff --git a/src/LorentzVector.cc b/src/LorentzVector.cc
index 745bf6b..926381c 100644
--- a/src/LorentzVector.cc
+++ b/src/LorentzVector.cc
@@ -1,33 +1,33 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include "HEJ/LorentzVector.hh"
#include <cmath>
#include <utility>
namespace HEJ {
// see eq:P_massive, eq:P_massive_p, eq:P_massive_q, , eq:P_massive_plus,
// eq:P_massive_minus in the developer manual
std::pair<CLHEP::HepLorentzVector, CLHEP::HepLorentzVector>
split_into_lightlike(CLHEP::HepLorentzVector const & P) {
const double pt = P.perp();
if(P.plus() > 0){
const double y = std::log(P.plus()/pt);
const double E = P.m2()/(2.*P.plus());
return std::make_pair(
CLHEP::HepLorentzVector{P.x(), P.y(), pt*std::sinh(y), pt*std::cosh(y)},
- CLHEP::HepLorentzVector{0., 0., E, -E}
+ CLHEP::HepLorentzVector{0., 0., -E, E}
);
} else {
const double y = std::log(pt/P.minus());
const double E = P.m2()/(2.*P.minus());
return std::make_pair(
CLHEP::HepLorentzVector{P.x(), P.y(), pt*std::sinh(y), pt*std::cosh(y)},
- CLHEP::HepLorentzVector{0., 0., E, +E}
+ CLHEP::HepLorentzVector{0., 0., +E, E}
);
}
}
}
diff --git a/src/Wjets.cc b/src/Wjets.cc
index 7bd2f54..34861c5 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1100 +1,830 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <array>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/jets.hh"
#include "HEJ/Tensor.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/exceptions.hh"
// generated headers
#include "HEJ/currents/jW_uno.hh"
#include "HEJ/currents/W_central_qqbar.hh"
#include "HEJ/currents/W_extremal_qqbar.hh"
+#include "HEJ/currents/extremal_W_central_qqbar.hh"
using HEJ::Tensor;
using HEJ::init_sigma_index;
using HEJ::metric;
using HEJ::rank3_current;
using HEJ::rank5_current;
using HEJ::eps;
using HEJ::to_tensor;
using HEJ::Helicity;
using HEJ::angle;
using HEJ::square;
using HEJ::flip;
using HEJ::ParticleProperties;
namespace helicity = HEJ::helicity;
namespace { // Helper Functions
// FKL W Helper Functions
double WProp (const HLV & plbar, const HLV & pl, ParticleProperties const & wprop){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
+ COM(0.,1.)*wprop.mass*wprop.width);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
namespace {
// FKL current including W emission off negative helicities
// See eq. (87) {eq:jW-} in developer manual
// Note that the terms are rearranged
Tensor<1> jW_minus(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl
){
using HEJ::helicity::minus;
const double tWin = (pa-pl-plbar).m2();
const double tWout = (p1+pl+plbar).m2();
// C++ arithmetic operators are evaluated left-to-right,
// so the following first computes complex scalar coefficients,
// which then multiply a current, reducing the number
// of multiplications
return 2.*(
+ angle(p1, pl)*square(p1, plbar)/tWout
+ square(pa, plbar)*angle(pa, pl)/tWin
)*HEJ::current(p1, pa, helicity::minus)
+ 2.*angle(p1, pl)*square(pl, plbar)/tWout
*HEJ::current(pl, pa, helicity::minus)
+ 2.*square(pa, plbar)*angle(pl, plbar)/tWin
*HEJ::current(p1, plbar, helicity::minus);
}
}
// FKL current including W emission
// see eqs. (87), (88) {eq:jW-}, {eq:jW+} in developer manual
Tensor<1> jW(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl,
Helicity h
){
if(h == helicity::minus) {
return jW_minus(pa, p1, plbar, pl);
}
return jW_minus(pa, p1, pl, plbar).complex_conj();
}
/**
- * @brief Contraction of the \tilde{U}_1 tensor
- *
- * This is the contraction of the tensor defined in eq:U1tensor
- * in the developer manual with the uno gluon polarisation vector,
- * the leptonic and the partonic current (see eq:wunocurrent) in the
- * developer manual)
- */
- COM U1contr(
- HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
- HLV const & pa, Helicity h1,
- HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
- ) {
- using helicity::plus;
- using helicity::minus;
-
- if(h1 == plus) {
- if(h2 == plus) {
- if(hg == plus) {
- return HEJ::U1<plus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::U1<plus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- else {
- if(hg == helicity::plus) {
- return HEJ::U1<plus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::U1<plus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- }
- else {
- if(h2 == helicity::plus) {
- if(hg == helicity::plus) {
- return HEJ::U1<minus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::U1<minus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- else {
- if(hg == helicity::plus) {
- return HEJ::U1<minus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::U1<minus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- }
- }
-
- /**
- * @brief Contraction of the \tilde{U}_1 tensor
- *
- * This is the contraction of the tensor defined in eq:U2tensor
- * in the developer manual with the uno gluon polarisation vector,
- * the leptonic and the partonic current (see eq:wunocurrent) in the
- * developer manual)
- */
- COM U2contr(
- HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
- HLV const & pa, Helicity h1,
- HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
- ) {
- using helicity::plus;
- using helicity::minus;
-
- if(h1 == plus) {
- if(h2 == plus) {
- if(hg == plus) {
- return HEJ::U2<plus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::U2<plus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- else {
- if(hg == helicity::plus) {
- return HEJ::U2<plus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::U2<plus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- }
- else {
- if(h2 == helicity::plus) {
- if(hg == helicity::plus) {
- return HEJ::U2<minus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::U2<minus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- else {
- if(hg == helicity::plus) {
- return HEJ::U2<minus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::U2<minus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- }
- }
-
- /**
- * @brief Contraction of the \tilde{U}_1 tensor
- *
- * This is the contraction of the tensor defined in eq:Ltensor
- * in the developer manual with the uno gluon polarisation vector,
- * the leptonic and the partonic current (see eq:wunocurrent) in the
- * developer manual)
- */
- COM Lcontr(
- HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
- HLV const & pa, Helicity h1,
- HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
- ) {
- using helicity::plus;
- using helicity::minus;
-
- if(h1 == plus) {
- if(h2 == plus) {
- if(hg == plus) {
- return HEJ::L<plus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::L<plus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- else {
- if(hg == helicity::plus) {
- return HEJ::L<plus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::L<plus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- }
- else {
- if(h2 == helicity::plus) {
- if(hg == helicity::plus) {
- return HEJ::L<minus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::L<minus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- else {
- if(hg == helicity::plus) {
- return HEJ::L<minus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- else {
- return HEJ::L<minus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
- }
- }
- }
- }
-
- /**
* @brief New implementation of unordered W+Jets current
*
* See eq:wunocurrent in the developer manual for the definition
* of this current
*/
+ template<Helicity h1, Helicity h2, Helicity pol>
double jM2Wuno(
HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
- HLV const & pa, Helicity const h1,
- HLV const & p2, HLV const & pb, Helicity const h2, Helicity const pol,
+ HLV const & pa, HLV const & p2, HLV const & pb,
ParticleProperties const & wprop
){
using HEJ::C_A;
using HEJ::C_F;
- const COM U1 = U1contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
- const COM U2 = U2contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
- const COM L = Lcontr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
+ const COM U1 = HEJ::U1<h1, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
+ const COM U2 = HEJ::U2<h1, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
+ const COM L = HEJ::L<h1, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM X = U1 - L;
const COM Y = U2 + L;
double amp = C_A*C_F*C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
amp *= WProp(plbar, pl, wprop);
const HLV q1g = pa-pl-plbar-p1-pg;
const HLV q2 = p2-pb;
const double t1 = q1g.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
/**
* @brief Contraction of W + extremal qqbar current
*
* See eq:crossed in the developer manual for the definition
* of the W + extremal qqbar current.
*
*/
template<Helicity h1, Helicity h2, Helicity hg>
double jM2_W_extremal_qqbar(
HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
HLV const & pqbar, HLV const & p3, HLV const & pb,
ParticleProperties const & wprop
){
using HEJ::C_A;
using HEJ::C_F;
const COM U1Xcontr = HEJ::U1X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM U2Xcontr = HEJ::U2X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM LXcontr = HEJ::LX<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM X = U1Xcontr - LXcontr;
const COM Y = U2Xcontr + LXcontr;
double amp = C_A*C_F*C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
amp *= WProp(plbar, pl, wprop);
// what do I have to put here?
const HLV q1 = pg-pl-plbar-pq-pqbar;
const HLV q2 = p3-pb;
const double t1 = q1.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
- Tensor <2> MCross(HLV pa, HLV pq, HLV pqbar, std::vector<HLV> partons,
- Helicity hq, int nabove
- ){
- //@TODO Simplify the calculation below Maybe combine with MCrossW?
- HLV q1;
- q1=pa;
- for(int i=0;i<nabove+1;++i){
- q1-=partons.at(i);
- }
-
- double t2=(q1-pqbar).m2();
-
- //Blank 3 gamma Current
- Tensor<3> J323 = rank3_current(pq,pqbar,hq);
-
- // 2 gamma current (with 1 contraction already).
- Tensor<2> XCroCont = J323.contract((q1-pqbar),2)/(t2);
-
- //Initialise the Crossed Vertex
- Tensor<2> Xcro(0.);
-
- for(int mu=0; mu<4;++mu){
- for(int nu=0;nu<4;++nu){
- Xcro(mu,nu) = XCroCont(nu,mu);
- }
- }
-
- return Xcro;
- }
-
- // Helper Functions Calculate the Uncrossed Contribution
- Tensor <2> MUncross(HLV pa, HLV pq,HLV pqbar, std::vector<HLV> partons,
- Helicity hq, int nabove
- ){
- //@TODO Simplify the calculation below Maybe combine with MUncrossW?
- HLV q1;
- q1=pa;
- for(int i=0;i<nabove+1;++i){
- q1-=partons.at(i);
- }
- double t2 = (q1-pq).m2();
-
- //Blank 3 gamma Current
- Tensor<3> J323 = rank3_current(pq,pqbar,hq);
-
- // 2 gamma currents (with 1 contraction already).
- Tensor<2> XUncCont = J323.contract((q1-pq),2)/t2;
-
- //Initialise the Uncrossed Vertex
- Tensor<2> Xunc(0.);
-
- for(int mu=0; mu<4;++mu){
- for(int nu=0;nu<4;++nu){
- Xunc(mu,nu) = -XUncCont(mu,nu);
- }
- }
-
- return Xunc;
- }
-
- // Helper Functions Calculate the Eikonal Contributions
- Tensor <2> MSym(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
- std::vector<HLV> partons, Helicity hq, int nabove
- ){
- //@TODO Simplify the calculation below Maybe combine with MsymW?
- HLV q1, q3;
- q1=pa;
- for(int i=0;i<nabove+1;++i){
- q1-=partons.at(i);
- }
- q3 = q1-pq-pqbar;
- double t1 = (q1).m2();
- double t3 = (q3).m2();
-
- double s23 = (pq+pqbar).m2();
- double sa2 = (pa+pq).m2();
- double sa3 = (pa+pqbar).m2();
- double s12 = (p1+pq).m2();
- double s13 = (p1+pqbar).m2();
- double sb2 = (pb+pq).m2();
- double sb3 = (pb+pqbar).m2();
- double s42 = (p4+pq).m2();
- double s43 = (p4+pqbar).m2();
-
- Tensor<1> qqxCur = HEJ::current(pq, pqbar, hq);
-
- // // 1a gluon emisson Contribution
- Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13))+ pa*(t1/(sa2+sa3)));
- Tensor<2> X1aCont = X1a.contract(qqxCur,3);
-
- // //4b gluon emission Contribution
- Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43)) + pb*(t3/(sb2+sb3)));
- Tensor<2> X4bCont = X4b.contract(qqxCur,3);
-
- // New Formulation Corresponding to New Analytics
- Tensor<3> X3g1 = outer(q1+pq+pqbar, metric());
- Tensor<3> X3g2 = outer(q3-pq-pqbar, metric());
- Tensor<3> X3g3 = outer(q1+q3, metric());
-
- // Note the contraction of indices changes term by term
- Tensor<2> X3g1Cont = X3g1.contract(qqxCur,3);
- Tensor<2> X3g2Cont = X3g2.contract(qqxCur,2);
- Tensor<2> X3g3Cont = X3g3.contract(qqxCur,1);
-
- Tensor<2>Xsym(0.);
-
- for(int mu=0; mu<4;++mu){
- for(int nu=0;nu<4;++nu){
- Xsym(mu, nu) = COM(0,1) * ( (X3g1Cont(nu,mu) + X3g2Cont(mu,nu)
- - X3g3Cont(nu,mu)) + (X1aCont(mu,nu) - X4bCont(mu,nu)) );
- }
- }
- return Xsym/s23;
- }
-
//! W+Jets FKL Contributions
/**
* @brief W+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_\mu.
* Handles all possible incoming states.
*/
double jW_j( HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
bool aqlineb, bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out);
const double WPropfact = WProp(plbar, pl, wprop);
const auto j_W = COM{0,-1}*jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
double Msqr = 0.;
for(const auto h: {plus, minus}) {
const auto j = HEJ::current(p2out, p2in, h);
Msqr += abs2(j_W.contract(j, 1));
}
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*Msqr/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
} // Anonymous Namespace
double ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop);
}
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true, wprop);
}
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop);
}
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true, wprop);
}
double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param pg Unordered Gluon momenta
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_{uno}_\mu. Ie, unordered with W emission opposite side.
* Handles all possible incoming states.
*/
double jW_juno(HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out,
HLV p2in, HLV pg, bool aqlineb, bool aqlinef,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out-pg);
const HLV q3=-(p2in-p2out);
const Helicity fhel = aqlinef?plus:minus;
const auto j_W = jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
const auto mj2p = HEJ::current(p2out, p2in, flip(fhel));
const auto mj2m = HEJ::current(p2out, p2in, fhel);
const auto jgbp = HEJ::current(pg, p2in, flip(fhel));
const auto jgbm = HEJ::current(pg, p2in, fhel);
const auto j2gp = HEJ::current(p2out, pg, flip(fhel));
const auto j2gm = HEJ::current(p2out, pg, fhel);
// Dot products of these which occur again and again
COM MWmp=j_W.dot(mj2p); // And now for the Higgs ones
COM MWmm=j_W.dot(mj2m);
const auto qsum = to_tensor(q2+q3);
const auto p1o = to_tensor(p1out);
const auto p1i = to_tensor(p1in);
const auto p2o = to_tensor(p2out);
const auto p2i = to_tensor(p2in);
const auto Lmm=( (-1.)*qsum*(MWmm) + (-2.*COM{j_W.dot(pg)})*mj2m + 2.*COM{mj2m.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmm/2. ) )/q3.m2();
const auto Lmp=( (-1.)*qsum*(MWmp) + (-2.*COM{j_W.dot(pg)})*mj2p + 2.*COM{mj2p.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmp/2. ) )/q3.m2();
const auto U1mm=(COM{jgbm.dot(j_W)}*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
const auto U1mp=(COM{jgbp.dot(j_W)}*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
const auto U2mm=((-1.)*COM{j2gm.dot(j_W)}*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
const auto U2mp=((-1.)*COM{j2gp.dot(j_W)}*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
double amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
ampsq*=WProp(plbar, pl, wprop);
return ampsq/((16)*(q2.m2()*q1.m2()));
}
}
double ME_W_unob_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true, wprop);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false, wprop);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true, wprop);
}
namespace{
+
+ // helicity sum helper for jWuno_j(...)
+ template<Helicity h1>
+ double jWuno_j_helsum(
+ HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
+ HLV const & pa, HLV const & p2, HLV const & pb,
+ ParticleProperties const & wprop
+ ){
+ using helicity::minus;
+ using helicity::plus;
+
+ const double ME2h1pp = jM2Wuno<h1, plus, plus>(
+ pg, p1, plbar, pl, pa, p2, pb, wprop
+ );
+
+ const double ME2h1pm = jM2Wuno<h1, plus, minus>(
+ pg, p1, plbar, pl, pa, p2, pb, wprop
+ );
+
+ const double ME2h1mp = jM2Wuno<h1, minus, plus>(
+ pg, p1, plbar, pl, pa, p2, pb, wprop
+ );
+
+ const double ME2h1mm = jM2Wuno<h1, minus, minus>(
+ pg, p1, plbar, pl, pa, p2, pb, wprop
+ );
+
+ return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm;
+ }
+
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same side.
* Handles all possible incoming states. Note this handles both forward and back-
* -ward Wuno emission. For forward, ensure p1out is the uno and W emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean up
* ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(HLV pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
HLV p2out, HLV p2in, bool aqlineb,
ParticleProperties const & wprop
){
- //Calculate different Helicity choices
- const Helicity h = aqlineb?helicity::plus:helicity::minus;
- double ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
- helicity::plus,helicity::plus, wprop);
- double ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
- helicity::plus,helicity::minus, wprop);
- double ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
- helicity::minus,helicity::plus, wprop);
- double ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
- helicity::minus,helicity::minus, wprop);
+ const auto helsum = aqlineb?
+ jWuno_j_helsum<helicity::plus>:
+ jWuno_j_helsum<helicity::minus>;
//Helicity sum and average over initial states
- return (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
+ return helsum(pg, p1out,plbar,pl,p1in,p2out,p2in, wprop)/(4.*HEJ::C_A*HEJ::C_A);
}
}
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
// helicity sum helper for jWqqx_j(...)
template<Helicity h1>
double jWqqx_j_helsum(
HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
HLV const & pqbar, HLV const & p3, HLV const & pb,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const double ME2h1pp = jM2_W_extremal_qqbar<h1, plus, plus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double ME2h1pm = jM2_W_extremal_qqbar<h1, plus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double ME2h1mp = jM2_W_extremal_qqbar<h1, minus, plus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double ME2h1mm = jM2_W_extremal_qqbar<h1, minus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm;
}
/**
* @brief W+Jets Extremal qqx Contributions, function to handle all incoming types.
* @param pgin Incoming gluon which will split into qqx.
* @param pqout Quark of extremal qqx outgoing (W-Emission).
* @param plbar Outgoing anti-lepton momenta
* @param pl Outgoing lepton momenta
* @param pqbarout Anti-quark of extremal qqx pair. (W-Emission)
* @param pout Outgoing Particle 2 (end of FKL chain)
* @param p2in Incoming Particle 2
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W_{qqx} ^\mu j_\mu. Ie, Ex-QQX with W emission same side.
* Handles all possible incoming states. Calculated via crossing symmetry from jWuno_j.
*/
double jWqqx_j(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef,
ParticleProperties const & wprop
){
const auto helsum = aqlinef?
jWqqx_j_helsum<helicity::plus>:
jWqqx_j_helsum<helicity::minus>;
//Helicity sum and average over initial states.
double ME2 = helsum(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, wprop)/
(4.*HEJ::C_A*HEJ::C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
}
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
}
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm1(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below. (Less Tensor class use?)
double t1 = (p3-pb)*(p3-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p3-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,2)/t1;
return gqqCur*(-1);
}
//Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm2(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double t1 = (p2-pb)*(p2-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb,refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p2-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm3(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s23 = (p2+p3)*(p2+p3);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank1 = outer(p2+p3, metric())/s23;
Tensor<3> qqCurBlank2 = outer(pb, metric())/s23;
Tensor<1> Cur23 = HEJ::current(p2, p3,hel2);
Tensor<2> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2,
HLV p3,HLV plbar, HLV pl, bool aqlinepa,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
// 2 independent helicity choices (complex conjugation related).
Tensor<1> TMmmm1 = qggm1(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,minus,minus, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,minus,plus, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa?plus:minus);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1);
COM Mmmm2 = TMmmm2.contract(cur1a,1);
COM Mmmm3 = TMmmm3.contract(cur1a,1);
COM Mpmm1 = TMpmm1.contract(cur1a,1);
COM Mpmm2 = TMpmm2.contract(cur1a,1);
COM Mpmm3 = TMpmm3.contract(cur1a,1);
//Colour factors:
COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
cm1m1=8./3.;
cm2m2=8./3.;
cm3m3=6.;
cm1m2 =-1./3.;
cm1m3 = -3.*COM(0.,1.);
cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
double Mmmm = real( cm1m1*pow(abs(Mmmm1),2) + cm2m2*pow(abs(Mmmm2),2)
+ cm3m3*pow(abs(Mmmm3),2) + 2.*real(cm1m2*Mmmm1*conj(Mmmm2))
+ 2.*real(cm1m3*Mmmm1*conj(Mmmm3))
+ 2.*real(cm2m3*Mmmm2*conj(Mmmm3)) );
double Mpmm = real( cm1m1*pow(abs(Mpmm1),2) + cm2m2*pow(abs(Mpmm2),2)
+ cm3m3*pow(abs(Mpmm3),2) + 2.*real(cm1m2*Mpmm1*conj(Mpmm2))
+ 2.*real(cm1m3*Mpmm1*conj(Mpmm3))
+ 2.*real(cm2m3*Mpmm2*conj(Mpmm3)) );
// Divide by WProp
const double WPropfact = WProp(plbar, pl, wprop);
return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
}
namespace {
// helper function for matrix element for W + Jets with central qqbar
template<HEJ::Helicity h1a, HEJ::Helicity h4b>
double amp_WCenqqx_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p4,
HLV const & pq, HLV const & pqbar,
HLV const & pl, HLV const & plbar,
HLV const & q11, HLV const & q12
) {
const COM sym = HEJ::M_sym_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM cross = HEJ::M_cross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM uncross = HEJ::M_uncross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
// Colour factors
static constexpr COM cmsms = 3.;
static constexpr COM cmumu = 4./3.;
static constexpr COM cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0., 3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
static constexpr COM cmumc = -1./6.;
return real(
cmsms*pow(abs(sym),2)
+cmumu*pow(abs(uncross),2)
+cmcmc*pow(abs(cross),2)
)
+2.*real(cmsmu*sym*conj(uncross))
+2.*real(cmsmc*sym*conj(cross))
+2.*real(cmumc*uncross*conj(cross))
;
}
}
-// matrix element for W + Jets with central qqbar
+// matrix element for W + Jets with W emission off central qqbar
double ME_WCenqqx_qq(
HLV pa, HLV pb, HLV pl, HLV plbar, std::vector<HLV> partons,
bool /* aqlinepa */, bool /* aqlinepb */, bool qqxmarker, int nabove,
ParticleProperties const & wprop
) {
using helicity::plus;
using helicity::minus;
CLHEP::HepLorentzVector q1 = pa;
for(int i = 0; i <= nabove; ++i){
q1 -= partons[i];
}
const auto qq = HEJ::split_into_lightlike(q1);
// since q1.m2() < 0 the following assertion is always true
// see documentation for split_into_lightlike
assert(qq.second.e() < 0);
HLV pq, pqbar;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];
}
const HLV p1 = partons.front();
const HLV p4 = partons.back();
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
// the first helicity label is for aqlinepa == true,
// the second one for aqlinepb == true
// In principle the corresponding helicity should be flipped
// if either of them is false, but we sum up all helicities,
// so this has no effect in the end.
const double amp_mm = amp_WCenqqx_qq<minus, minus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_mp = amp_WCenqqx_qq<minus, plus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pm = amp_WCenqqx_qq<plus, minus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pp = amp_WCenqqx_qq<plus, plus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double t1 = q1.m2();
const double t3 = (q1-pl-plbar-pq-pqbar).m2();
const double amp = WProp(plbar, pl, wprop)*(
amp_mm+amp_mp+amp_pm+amp_pp
)/(9.*4.*t1*t1*t3*t3);
return amp;
}
-// no wqq emission
-double ME_W_Cenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
- bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
- int nbelow, bool forwards, ParticleProperties const & wprop
+// helper function for matrix element for W + Jets with central qqbar
+// W emitted off extremal parton
+// @TODO: code duplication with amp_WCenqqx_qq
+template<HEJ::Helicity h1a, HEJ::Helicity hqqbar>
+double amp_W_Cenqqx_qq(
+ HLV const & pa, HLV const & p1,
+ HLV const & pb, HLV const & p4,
+ HLV const & pq, HLV const & pqbar,
+ HLV const & pl, HLV const & plbar,
+ HLV const & q11, HLV const & q12
+ ) {
+
+ const COM crossed = HEJ::M_cross<h1a, hqqbar>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
+ );
+ const COM uncrossed = HEJ::M_qbar<h1a, hqqbar>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
+ );
+ const COM sym = HEJ::M_sym<h1a, hqqbar>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
+ );
+
+ //Colour factors:
+ static constexpr COM cmsms = 3.;
+ static constexpr COM cmumu = 4./3.;
+ static constexpr COM cmcmc = 4./3.;
+ static constexpr COM cmsmu = COM{0.,3./2.};
+ static constexpr COM cmsmc = COM{0.,-3./2.};
+ static constexpr COM cmumc = -1./6.;
+
+ return real(
+ cmsms*pow(abs(sym),2)
+ +cmumu*pow(abs(uncrossed),2)
+ +cmcmc*pow(abs(crossed),2)
+ )
+ +2.*real(cmsmu*sym*conj(uncrossed))
+ +2.*real(cmsmc*sym*conj(crossed))
+ +2.*real(cmumc*uncrossed*conj(crossed))
+ ;
+}
+
+// matrix element for W + Jets with W emission *not* off central qqbar
+double ME_W_Cenqqx_qq(
+ HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
+ bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
+ int nbelow, bool forwards, ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
- init_sigma_index();
-
if (!forwards){ //If Emission from Leg a instead, flip process.
std::swap(pa, pb);
std::reverse(partons.begin(),partons.end());
std::swap(aqlinepa, aqlinepb);
qqxmarker = !qqxmarker;
std::swap(nabove,nbelow);
}
- HLV pq, pqbar, p1,p4;
+ HLV q1=pa;
+ for(int i=0;i<nabove+1;++i){
+ q1-=partons.at(i);
+ }
+ const auto qq = HEJ::split_into_lightlike(q1);
+
+ HLV pq, pqbar;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
- p1 = partons.front();
- p4 = partons.back();
-
- Tensor<1> T1am(0.), T1ap(0.);
- if(!(aqlinepa)){
- T1ap = HEJ::current(p1, pa, plus);
- T1am = HEJ::current(p1, pa, minus);}
- else if(aqlinepa){
- T1ap = HEJ::current(pa, p1, plus);
- T1am = HEJ::current(pa, p1, minus);}
-
- Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb?plus:minus);
-
- // Calculate the 3 separate contributions to the effective vertex
- Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, minus, nabove);
- Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, minus, nabove);
- Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, minus, nabove);
-
- Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, plus, nabove);
- Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, plus, nabove);
- Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, plus, nabove);
-
- // (- - hel choice)
- COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1));
- COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1));
- COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1));
- // (- + hel choice)
- COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1));
- COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1));
- COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1));
- // (+ - hel choice)
- COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1));
- COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1));
- COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1));
- // (+ + hel choice)
- COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1));
- COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1));
- COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1));
+ // we assume that the W is emitted off a quark line
+ // if this is not the case, we have to apply CP conjugation,
+ // which is equivalent to swapping lepton and antilepton momenta
+ if(aqlinepb) std::swap(pl, plbar);
- //Colour factors:
- COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
- cmsms=3.;
- cmumu=4./3.;
- cmcmc=4./3.;
- cmsmu =3./2.*COM(0.,1.);
- cmsmc = -3./2.*COM(0.,1.);
- cmumc = -1./6.;
-
- // Work Out Interference in each case of helicity:
- double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
- +cmumu*pow(abs(M_mmUnc),2)
- +cmcmc*pow(abs(M_mmCro),2)
- +2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
- +2.*real(cmsmc*M_mmSym*conj(M_mmCro))
- +2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
-
- double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
- +cmumu*pow(abs(M_mpUnc),2)
- +cmcmc*pow(abs(M_mpCro),2)
- +2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
- +2.*real(cmsmc*M_mpSym*conj(M_mpCro))
- +2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
-
- double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
- +cmumu*pow(abs(M_pmUnc),2)
- +cmcmc*pow(abs(M_pmCro),2)
- +2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
- +2.*real(cmsmc*M_pmSym*conj(M_pmCro))
- +2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
-
- double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
- +cmumu*pow(abs(M_ppUnc),2)
- +cmcmc*pow(abs(M_ppCro),2)
- +2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
- +2.*real(cmsmc*M_ppSym*conj(M_ppCro))
- +2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
-
- double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
-
- HLV q1,q3;
- q1=pa;
- for(int i=0;i<nabove+1;++i){
- q1-=partons.at(i);
- }
- q3 = q1 - pq - pqbar;
+ const HLV p1 = partons.front();
+ const HLV p4 = partons.back();
- double t1 = (q1).m2();
- double t3 = (q3).m2();
+ // helicity labels are for aqlinepa == aqlineb == false,
+ // In principle the helicities should be flipped
+ // if either of them is true, but we sum up all helicities,
+ // so this has no effect in the end.
+ const double amp_mm = amp_W_Cenqqx_qq<minus, minus>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
+ );
+ const double amp_mp = amp_W_Cenqqx_qq<minus, plus>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
+ );
+ const double amp_pm = amp_W_Cenqqx_qq<plus, minus>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
+ );
+ const double amp_pp = amp_W_Cenqqx_qq<plus, plus>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
+ );
- //Divide by t-channels
- amp/=(t1*t1*t3*t3);
+ const double t1 = q1.m2();
+ const double t3 = (q1 - pq - pqbar).m2();
- //Divide by WProp
- amp*=WProp(plbar, pl, wprop);
+ const double amp= WProp(plbar, pl, wprop)*(
+ amp_mm+amp_mp+amp_pm+amp_pp
+ )/(9.*4.*t1*t1*t3*t3);
return amp;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:01 PM (1 h, 48 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486434
Default Alt Text
(85 KB)

Event Timeline