Page MenuHomeHEPForge
Files F11221591

No OneTemporary
Actions

View Options

diff --git a/current_generator/j_jqqbar.frm b/current_generator/j_jqqbar.frm
deleted file mode 100644
index dde640a..0000000
--- a/current_generator/j_jqqbar.frm
+++ /dev/null
@@ -1,48 +0,0 @@
-*/**
-* \brief Contraction of FKL current with extremal qqbar emission current
-*
-* TODO: transcribe formulas to developer manual
-*
-* \authors The HEJ collaboration (see AUTHORS for details)
-* \date 2020
-* \copyright GPLv2 or later
-*/
-#include- include/helspin.frm
-#include- include/write.frm
-
-v pb,p2,p3,pa,p1,p;
-i mu,nu,sigma;
-s h;
-
-* Equation 3.23 in James Cockburn's Thesis.
-l qggm1 = Current(-1, p2, mu, p3-pb, nu, p3);
-l qggm2 = -Current(-1, p2, nu, p2-pb, mu, p3);
-l qggm3 = -2*i_/m2(p2+p3)*(
- (p2(nu) + p3(nu))*d_(mu, sigma)
- - pb(mu)*d_(nu, sigma)
- + pb(sigma)*d_(mu, nu)
-)*Current(-1, p2, sigma, p3);
-.sort
-drop;
-
-#do i=1,3
- #do h1={+,-}
- #do hg={+,-}
- l [j_qggm`i' `h1'`hg'] = qggm`i'*Eps(`hg'1, nu, pa, pb)*Current(
- `h1'1, pa, mu, p1
- );
- #enddo
- #enddo
-#enddo
-id Current(h?, p1?, mu?, p?, nu?, p2?) = Current(h, p1, mu, p, nu, p2)/m2(p);
-
-#call ContractCurrents
-.sort
-format O4;
-format c;
-#call WriteHeader(`OUTPUT')
-#call WriteOptimised(`OUTPUT',j_qggm1,2,pb,p2,p3,pa,p1)
-#call WriteOptimised(`OUTPUT',j_qggm2,2,pb,p2,p3,pa,p1)
-#call WriteOptimised(`OUTPUT',j_qggm3,2,pb,p2,p3,pa,p1)
-#call WriteFooter(`OUTPUT')
-.end
diff --git a/current_generator/jqqbar_j.frm b/current_generator/jqqbar_j.frm
new file mode 100644
index 0000000..626885a
--- /dev/null
+++ b/current_generator/jqqbar_j.frm
@@ -0,0 +1,66 @@
+*/**
+* \brief Contraction of FKL current with extremal qqbar current
+*
+* \authors The HEJ collaboration (see AUTHORS for details)
+* \date 2020
+* \copyright GPLv2 or later
+*/
+#include- include/helspin.frm
+#include- include/write.frm
+
+* Obtained from juno_j.frm using crossing symmetry (see eq:crossing in developer manual)
+
+v pb,p2,pqbar,pq,pa,q2,q3,p,pr;
+i mu,nu,sigma;
+s h;
+
+#do h1={+,-}
+ l [U1X_j `h1'] = (
+ + i_*Current(`h1'1, pq, nu, pa)*i_*i_*Current(`h1'1, pa, mu, pqbar)
+ + 2*pq(nu)*i_*Current(`h1'1, pq, mu, pqbar)
+ )/m2(pq - pa);
+ l [U2X_j `h1'] = (
+ - 2*i_*Current(`h1'1, pq, mu, pqbar)*pqbar(nu)
+ - i_*Current(`h1'1, pq, mu, pa)*i_*i_*Current(`h1'1, pa, nu, pqbar)
+ )/m2(-pqbar + pa);
+ l [LX_j `h1'] = (
+ - (q2(nu) + q3(nu))*i_*Current(`h1'1, pq, mu, pqbar)
+ + 2*pa(mu)*i_*Current(`h1'1, pq, nu, pqbar)
+ - 2*i_*Current(`h1'1, pq, pa, pqbar)*d_(mu, nu)
+ + m2(q2)*i_*Current(`h1'1, pq, mu, pqbar)*(p2(nu)/m2(p2-pa) + pb(nu)/m2(pb-pa))
+ )/m2(q3);
+#enddo
+.sort
+drop;
+
+* multiply with polarisation vector and other currents
+#do h1={+,-}
+ #do hg={+,-}
+ #do TENSOR={U1X_j,U2X_j,LX_j}
+ l [`TENSOR' `h1'`hg'] = (
+ [`TENSOR' `h1']
+ *Current(-1, pb, mu, p2)
+ *Eps(`hg'1, nu, pr)
+ );
+ #enddo
+ #enddo
+#enddo
+
+* choice of auxiliary vector
+id Eps(h?, mu?, pr) = Eps(h, mu, pa, pq);
+
+multiply replace_(
+ q2,pq-pa+pqbar,
+ q3,pq+pqbar
+);
+.sort
+#call ContractCurrents
+.sort
+format O4;
+format c;
+#call WriteHeader(`OUTPUT')
+#call WriteOptimised(`OUTPUT',U1X_j,2,pa,pb,pq,pqbar,p2)
+#call WriteOptimised(`OUTPUT',U2X_j,2,pa,pb,pq,pqbar,p2)
+#call WriteOptimised(`OUTPUT',LX_j,2,pa,pb,pq,pqbar,p2)
+#call WriteFooter(`OUTPUT')
+.end
diff --git a/doc/developer_manual/biblio.bib b/doc/developer_manual/biblio.bib
index 71aae00..1e9f72b 100644
--- a/doc/developer_manual/biblio.bib
+++ b/doc/developer_manual/biblio.bib
@@ -1,194 +1,204 @@
@Article{Andersen:2008gc,
author = "Andersen, Jeppe R. and Del Duca, Vittorio and White, Chris
D.",
title = "{Higgs Boson Production in Association with Multiple Hard
Jets}",
journal = "JHEP",
volume = "02",
year = "2009",
pages = "015",
eprint = "0808.3696",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
doi = "10.1088/1126-6708/2009/02/015",
SLACcitation = "%%CITATION = 0808.3696;%%"
}
@Article{Andersen:2008ue,
author = "Andersen, Jeppe R. and White, Chris D.",
title = "{A New Framework for Multijet Predictions and its
application to Higgs Boson production at the LHC}",
journal = "Phys. Rev.",
volume = "D78",
year = "2008",
pages = "051501",
eprint = "0802.2858",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
doi = "10.1103/PhysRevD.78.051501",
SLACcitation = "%%CITATION = 0802.2858;%%"
}
@article{Alwall:2006yp,
author = "Alwall, Johan and others",
title = "{A Standard format for Les Houches event files}",
booktitle = "{Monte Carlos for the LHC: A Workshop on the Tools for
LHC Event Simulation (MC4LHC) Geneva, Switzerland, July
17-16, 2006}",
journal = "Comput. Phys. Commun.",
volume = "176",
year = "2007",
pages = "300-304",
doi = "10.1016/j.cpc.2006.11.010",
eprint = "hep-ph/0609017",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "FERMILAB-PUB-06-337-T, CERN-LCGAPP-2006-03",
SLACcitation = "%%CITATION = HEP-PH/0609017;%%"
}
@article{Andersen:2017kfc,
author = "Andersen, Jeppe R. and Hapola, Tuomas and Maier, Andreas
and Smillie, Jennifer M.",
title = "{Higgs Boson Plus Dijets: Higher Order Corrections}",
journal = "JHEP",
volume = "09",
year = "2017",
pages = "065",
doi = "10.1007/JHEP09(2017)065",
eprint = "1706.01002",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "EDINBURGH-2017-11, IPPP-17-33, DCPT-17-66, MCNET-17-9",
SLACcitation = "%%CITATION = ARXIV:1706.01002;%%"
}
@article{Andersen:2018tnm,
author = "Andersen, Jeppe R. and Hapola, Tuomas and Heil, Marian
and Maier, Andreas and Smillie, Jennifer M.",
title = "{Higgs-boson plus Dijets: Higher-Order Matching for
High-Energy Predictions}",
year = "2018",
eprint = "1805.04446",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "DCPT/18/66, IPPP/18/33, MCnet-18-10, DCPT-18-66,
IPPP-18-33, MCNET-18-10",
SLACcitation = "%%CITATION = ARXIV:1805.04446;%%"
}
@article{Andersen:2012gk,
author = "Andersen, Jeppe R. and Hapola, Tuomas and Smillie,
Jennifer M.",
title = "{W Plus Multiple Jets at the LHC with High Energy Jets}",
journal = "JHEP",
volume = "09",
year = "2012",
pages = "047",
doi = "10.1007/JHEP09(2012)047",
eprint = "1206.6763",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "EDINBURGH-2012-13, IPPP-12-45, DCPT-12-90,
CP3-ORIGINS-2012-017, DIAS-2012-18, --DIAS-2012-18",
SLACcitation = "%%CITATION = ARXIV:1206.6763;%%"
}
@article{Andersen:2016vkp,
author = "Andersen, Jeppe R. and Medley, Jack J. and Smillie,
Jennifer M.",
title = "{$Z/\gamma^{∗}$ plus multiple hard jets in high energy
collisions}",
journal = "JHEP",
volume = "05",
year = "2016",
pages = "136",
doi = "10.1007/JHEP05(2016)136",
eprint = "1603.05460",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "EDINBURGH-2016-03, IPPP-16-19, DCPT-16-38, MCNET-16-08",
SLACcitation = "%%CITATION = ARXIV:1603.05460;%%"
}
@article{DelDuca:2003ba,
author = "Del Duca, V. and Kilgore, W. and Oleari, C. and Schmidt,
C.R. and Zeppenfeld, D.",
title = "{Kinematical limits on Higgs boson production via gluon
fusion in association with jets}",
journal = "Phys.Rev.",
volume = "D67",
pages = "073003",
doi = "10.1103/PhysRevD.67.073003",
year = "2003",
eprint = "hep-ph/0301013",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "DCPT-02-148, DFTT-20-2002, IPPP-02-74, MADPH-02-1276,
MSUHEP-20620",
SLACcitation = "%%CITATION = HEP-PH/0301013;%%",
}
@article{DelDuca:2001fn,
author = "Del Duca, V. and Kilgore, W. and Oleari, C. and Schmidt,
C. and Zeppenfeld, D.",
title = "{Gluon fusion contributions to H + 2 jet production}",
journal = "Nucl.Phys.",
volume = "B616",
pages = "367-399",
doi = "10.1016/S0550-3213(01)00446-1",
year = "2001",
eprint = "hep-ph/0108030",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "MADPH-01-1235, BNL-HET-01-28, MSUHEP-10709, DFTT-19-2001",
SLACcitation = "%%CITATION = HEP-PH/0108030;%%",
}
@article{Andersen:2011zd,
author = "Andersen, Jeppe R. and Lonnblad, Leif and Smillie,
Jennifer M.",
title = "{A Parton Shower for High Energy Jets}",
journal = "JHEP",
volume = "07",
year = "2011",
pages = "110",
doi = "10.1007/JHEP07(2011)110",
eprint = "1104.1316",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "CERN-PH-TH-2011-072, CP3-ORIGINS-2011-14,
EDINBURGH-2011-16, LU-TP-11-15, MCNET-11-12, LU-TP
--11-15",
SLACcitation = "%%CITATION = ARXIV:1104.1316;%%"
}
@article{Andersen:2017sht,
author = "Andersen, Jeppe R. and Brooks, Helen M. and Lönnblad,
Leif",
title = "{Merging High Energy with Soft and Collinear Logarithms
using HEJ and PYTHIA}",
journal = "JHEP",
volume = "09",
year = "2018",
pages = "074",
doi = "10.1007/JHEP09(2018)074",
eprint = "1712.00178",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "CPT/17/184, IPPP/17/92, LU-TP 17-38, MCnet-17-22,
CoEPP-MN-17-22, CPT-17-184, IPPP-17-92, LU-TP-17-38,
MCNET-17-22, COEPP-MN-17-22",
SLACcitation = "%%CITATION = ARXIV:1712.00178;%%"
}
@article{Hoeche:2019rti,
author = "Höche, Stefan and Prestel, Stefan and Schulz, Holger",
title = "{Simulation of vector boson plus many jet final states at
the high luminosity LHC}",
year = "2019",
eprint = "1905.05120",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "FERMILAB-PUB-19-192-T, LU-TP 19-14, MCNET-19-09",
SLACcitation = "%%CITATION = ARXIV:1905.05120;%%"
}
@article{Andersen:2020yax,
author = "Andersen, Jeppe R. and Black, James A. and Brooks, Helen M. and Byrne, Emmet P. and Maier, Andreas and Smillie, Jennifer M.",
title = "{Combined subleading high-energy logarithms and NLO accuracy for W production in association with multiple jets}",
eprint = "2012.10310",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "DESY-20-233, DCPT/20/134, IPPP/20/67, MCnet-20-27, SAGEX-20-28",
month = "12",
year = "2020"
}
+@article{Gunion:1985vca,
+ author = "Gunion, J. F. and Kunszt, Z.",
+ title = "{Improved Analytic Techniques for Tree Graph Calculations and the G g q anti-q Lepton anti-Lepton Subprocess}",
+ reportNumber = "OITS-296, UCD-85/3",
+ doi = "10.1016/0370-2693(85)90774-9",
+ journal = "Phys. Lett. B",
+ volume = "161",
+ pages = "333",
+ year = "1985"
+}
diff --git a/doc/developer_manual/currents.tex b/doc/developer_manual/currents.tex
index 9305029..f79529f 100644
--- a/doc/developer_manual/currents.tex
+++ b/doc/developer_manual/currents.tex
@@ -1,653 +1,660 @@
\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.\todo[inline]{Current implementation still uses
-different conventions (James Cockburn's thesis)}
+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}}
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[]{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 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_{nb}\ \frac{j_{\mu}(p_a,p_1)\ X^{ab\, \mu
\nu}_{{\rm cen}}(p_q,p_{\bar{q}},q_1,q_3)\
j_{\nu}(p_b,p_n)}{t_{a1}t_{bn}}.
\end{align}
where $X^{\mu \nu}_{\rm cen}$ is given by:
\begin{equation}
\label{eq:Xcen}
\begin{split}
X^{\mu \nu}_{\rm cen} ={}&\frac{f^{ced}T^c_{q\bar{q}}}{s_{q\bar{q}}}
\left(\eta^{\mu \nu} X_{sym}^\sigma + V^{\mu \nu \sigma}_{\bar{q}g} \right)
\bar{u}(p_q) \gamma^\sigma u(p_{\bar{q}}) \\ & \qquad + \frac{i
T^d_{qj}T^e_{j\bar{q}}}{(q_1-p_q)^2} X^{\mu\nu}_{\text{qprop}} - \frac{i
T^e_{qj}T^d_{j\bar{q}}}{(q_1-p_{\bar{q}})^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_n^\sigma}{s_{nq} + s_{n\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}}$.
\todo[inline]{Document formula for current contraction, like equation~(\ref{eq:S_uno})}
\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.
\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:wunocurrent}
\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 P_L |a\rangle }{s_{12}t_{aV}} + \frac{\langle 2|\nu (\slashed{p}_2+ \slashed{p}_1)\rho P_L (\slashed{p}_2+\slashed{p}_1 + \slashed{p}_V)\mu |a\rangle }{s_{12}s_{12V}} \\
&+ \frac{\langle 2|\rho P_L (\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 P_L |a\rangle }{t_{aV1}t_{aV}} + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_V- \slashed{p}_1)\rho P_L (\slashed{p}_a-\slashed{p}_1) \nu |a\rangle }{t_{a1V}t_{a1}} \\
&+ \frac{\langle 2|\rho P_L (\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_2)&\gamma_\rho\nu(p_{\bar{q}})\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_2)\gamma_\rho\nu(p_{\bar{q}})\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_2)\gamma_\rho\nu(p_{\bar{q}})\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:WunoY}
\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}
\subsubsection{\texorpdfstring{$W$}{W}+Extremal \texorpdfstring{$\mathbf{q\bar{q}}$}{qqbar}}
\todo{Update when included in $Z$ + jets}
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}}\,,
\end{split}\\
\label{eq:U2tensorX}
\begin{split}
\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}}\,,
\end{split}\\
\label{eq:LtensorX}
\begin{split}
\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 \texorpdfstring{$\mathbf{q\bar{q}}$}{qqbar} 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 \FORM inside the
\texttt{current\_generator/j\_Wqqbar\_j.frm}, which is detailed in
section~\ref{sec:contr_calc}.
The $W$-emission can be 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 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}}=&
- i \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{align}
\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{align}
are the colour factors of different contributions and $J_\text{V}$ is
given in equation~(\ref{eq:J_V}).
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_{a2} + s_{a3} + s_{aA} + s_{aB}} +
\frac{p^\sigma_1}{s_{12} + s_{13} + s_{1A} + s_{1B}}\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_{b2} + s_{b3} + s_{bA} + s_{bB}}+
\frac{p^\sigma_4}{s_{42} + s_{43} + s_{4A} + s_{4B}}\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! returns a contraction of the above tensor containing the
information from these 5 groupings of contributions (30 diagrams in total). It
is available through the generated header \texttt{j\_Wqqbar\_j.hh} (see
section~\ref{sec:cur_gen}).
\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!.
\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!.
\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:
diff --git a/include/HEJ/jets.hh b/include/HEJ/jets.hh
index 7cceff3..c7f9ee1 100644
--- a/include/HEJ/jets.hh
+++ b/include/HEJ/jets.hh
@@ -1,93 +1,93 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2022
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in pure jets.
*
* This file contains all the necessary functions to compute the
* current contractions for all valid pure jet HEJ processes, which
* so far is FKL and unordered processes. It will also contain some
* pure jet ME components used in other process ME calculations
*
*/
#pragma once
#include <array>
#include <cmath>
#include <complex>
#include <ostream>
#include <vector>
#include "CLHEP/Vector/LorentzVector.h"
namespace HEJ {
namespace currents {
using HLV = CLHEP::HepLorentzVector;
//! Helicity-summed current contraction
/**
* @param p1out Final-state momentum of first current
* @param p1in Initial-state momentum of first current
* @param p2out Final-state momentum of second current
* @param p2in Initial-state momentum of second current
* @returns Current contractions
*
* \internal
* See eq:S in the developer manual
*/
double ME_qq(
HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
);
//! Helicity-summed current contraction with unordered gluon emission
/**
* @param pg Momentum of gluon emitted off first current
* @param p1out Final-state momentum of first current
* @param p1in Initial-state momentum of first current
* @param p2out Final-state momentum of second current
* @param p2in Initial-state momentum of second current
* @returns Current contraction
*
* \internal
* See eq:S_uno in the developer manual
*/
double ME_unob_qq(
HLV const & pg,
HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
);
//! Helicity-summed current contraction with extremal quark/anti-quark emission
/**
- * @param pa Initial-state momentum of first current
- * @param pb Initial-state momentum of gluon splitting into qqbar
- * @param p1 Final-state momentum of first current
+ * @param pa Initial-state momentum of gluon splitting into qqbar
+ * @param pb Initial-state momentum of second current
+ * @param p1 More extremal (anti-)quark from splitting momentum
* @param p2 Less extremal (anti-)quark from splitting momentum
- * @param p3 More extremal (anti-)quark from splitting momentum
+ * @param p3 Final-state momentum of second current
* @returns Current contraction
*/
double ME_qqbar_qg(
HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2, HLV const & p3
);
//! Square of qq->qQQbarq Pure Jets Central qqbar Scattering Current
/**
* @param pa Momentum of incoming leg a
* @param pb Momentum of incoming leg b
* @param partons vector of outgoing partons
* @param qqbarmarker Is anti-quark further back in rapidity than quark
* (qqbar pair)
* @param nabove Number of gluons emitted above qqbar pair (back in rap)
* @returns Square of the current contractions for qq->q+QQbar+q
*/
double ME_Cenqqbar_qq(
HLV const & pa, HLV const & pb,
std::vector<HLV> const & partons,
bool qqbarmarker, std::size_t nabove
);
} // namespace currents
} // namespace HEJ
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 6d7fcdb..5bc18b8 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,2475 +1,2475 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2022
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <iterator>
#include <limits>
#include <unordered_map>
#include <utility>
#include "fastjet/PseudoJet.hh"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/Event.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "HEJ/WWjets.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Zjets.hh"
#include "HEJ/event_types.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/jets.hh"
#include "HEJ/utility.hh"
namespace HEJ {
namespace {
// Colour acceleration multiplier for gluons
// see eq:K_g in developer manual
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
// Colour acceleration multiplier for quarks
// see eq:K_q in developer manual
constexpr double K_q = C_F;
// Colour acceleration multipliers
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
){
if(type == pid::gluon) return K_g(pout, pin);
return K_q;
}
} // namespace
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*std::log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
return (
1. + alpha_s/(4.*M_PI)*BETA0*std::log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
std::vector <Weights> virtual_part=virtual_corrections(event);
if(tree_kin_part.size() != virtual_part.size()) {
throw std::logic_error("tree and virtuals have different sizes");
}
Weights sum = Weights{0., std::vector<double>(event.variations().size(), 0.)};
for(size_t i=0; i<tree_kin_part.size(); ++i) {
sum += tree_kin_part.at(i)*virtual_part.at(i);
}
return tree_param(event)*sum;
}
Weights MatrixElement::tree(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
double sum = 0.;
for(double i : tree_kin_part) {
sum += i;
}
return tree_param(event)*sum;
}
Weights MatrixElement::tree_param(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
std::vector<Weights> MatrixElement::virtual_corrections(Event const & event) const {
if(!event.valid_hej_state(param_.soft_pt_regulator)) {
return {Weights{0., std::vector<double>(event.variations().size(), 0.)}};
}
// only compute once for each renormalisation scale
std::unordered_map<double, std::vector<double> > known_vec;
std::vector<double> central_vec=virtual_corrections(event, event.central().mur);
known_vec.emplace(event.central().mur, central_vec);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
if(ME_it == end(known_vec)) {
known_vec.emplace(var.mur, virtual_corrections(event, var.mur));
}
}
// At this stage known_vec contains one vector of virtual corrections for each mur value
// Now put this into a vector of Weights
std::vector<Weights> result_vec;
for(size_t i=0; i<central_vec.size(); ++i) {
Weights result;
result.central = central_vec.at(i);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
result.variations.emplace_back(ME_it->second.at(i));
}
result_vec.emplace_back(result);
}
return result_vec;
}
template<class InputIterator>
std::vector <double> MatrixElement::virtual_corrections_interference(
InputIterator begin_parton, InputIterator end_parton,
fastjet::PseudoJet const & q0_t,
fastjet::PseudoJet const & q0_b,
const double mur
) const{
const double alpha_s = alpha_s_(mur);
auto qi_t = q0_t;
auto qi_b = q0_b;
double sum_top = 0.;
double sum_bot = 0.;
double sum_mix = 0.;
for(auto parton_it = begin_parton; parton_it != end_parton; ++parton_it){
Particle parton = *parton_it;
Particle parton_next = *(parton_it+1);
const double dy = parton_next.rapidity() - parton.rapidity();
const double tmp_top = omega0(alpha_s, mur, qi_t)*dy;
const double tmp_bot = omega0(alpha_s, mur, qi_b)*dy;
sum_top += tmp_top;
sum_bot += tmp_bot;
sum_mix += (tmp_top + tmp_bot) / 2.;
qi_t -= parton_next.p;
qi_b -= parton_next.p;
}
if (param_.nlo.enabled){
return {(sum_top), (sum_bot), (sum_mix)};
}
return {exp(sum_top), exp(sum_bot), exp(sum_mix)};
}
double MatrixElement::virtual_corrections_W(
Event const & event,
const double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
std::size_t first_idx = 0;
std::size_t last_idx = partons.size() - 1;
#ifndef NDEBUG
bool wc = true;
#endif
bool wqq = false;
// With extremal qqbar or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else if (event.type() == event_type::qqbar_exb) {
q -= partons[1].p;
++first_idx;
if (std::abs(partons[0].type) != std::abs(partons[1].type)){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else {
if(event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf){
--last_idx;
}
if (in[0].type != partons[0].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
std::size_t first_idx_qqbar = last_idx;
std::size_t last_idx_qqbar = last_idx;
//if qqbarMid event, virtual correction do not occur between qqbar pair.
if(event.type() == event_type::qqbar_mid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(std::abs(backquark->type) != std::abs((backquark+1)->type)) {
wqq=true;
#ifndef NDEBUG
wc=false;
#endif
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqbar = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(std::size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqbar) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
#ifndef NDEBUG
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf
);
#endif
if (param_.nlo.enabled){
double nlo_virtual =1.;
//Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) nlo_virtual+=exponent;
return nlo_virtual;
}
return std::exp(exponent);
}
std::vector <double> MatrixElement::virtual_corrections_WW(
Event const & event,
const double mur
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
assert(event.decays().size() == 2);
std::vector<fastjet::PseudoJet> plbar;
std::vector<fastjet::PseudoJet> pl;
for(auto const & decay_pair : event.decays()) {
auto const decay = decay_pair.second;
if(decay.at(0).type < 0) {
plbar.emplace_back(decay.at(0).p);
pl .emplace_back(decay.at(1).p);
}
else {
pl .emplace_back(decay.at(0).p);
plbar.emplace_back(decay.at(1).p);
}
}
fastjet::PseudoJet q_t = pa - partons[0].p - pl[0] - plbar[0];
fastjet::PseudoJet q_b = pa - partons[0].p - pl[1] - plbar[1];
auto const begin_parton = cbegin(partons);
auto const end_parton = cend(partons) - 1;
if (param_.nlo.enabled){
std::vector<double> virt_corrections_nlo {1.0,1.0,1.0}; //set virtual corrections to 1.
// Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) {
std::vector <double> virt_corrections_nlo_interference =
virtual_corrections_interference(
begin_parton, end_parton, q_t, q_b, mur
);
assert(
virt_corrections_nlo_interference.size()
== virt_corrections_nlo.size()
);
for(std::size_t i = 0; i < virt_corrections_nlo.size(); ++i){
virt_corrections_nlo[i] += virt_corrections_nlo_interference[i];
}
}
return virt_corrections_nlo;
}
return virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur);
}
std::vector <double> MatrixElement::virtual_corrections_Z_qq(
Event const & event,
const double mur,
Particle const & ZBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p;
fastjet::PseudoJet q_b = pa - partons[0].p;
auto begin_parton = cbegin(partons);
auto end_parton = cend(partons) - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q_t -= partons[1].p;
q_b -= partons[1].p;
++begin_parton;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--end_parton;
}
if (param_.nlo.enabled){
//set virtual corrections to 1.
std::vector<double> virt_corrections_nlo {1.0,1.0,1.0};
// Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) {
std::vector <double> virt_corrections_nlo_interference =
virtual_corrections_interference(
begin_parton, end_parton, q_t, q_b, mur
);
assert(
virt_corrections_nlo.size()
== virt_corrections_nlo_interference.size()
);
for(std::size_t i = 0; i < virt_corrections_nlo.size(); ++i){
virt_corrections_nlo[i] += virt_corrections_nlo_interference[i];
}
}
return virt_corrections_nlo;
}
return virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur);
}
double MatrixElement::virtual_corrections_Z_qg(
Event const & event,
const double mur,
Particle const & ZBoson,
const bool is_gq_event
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
// If this is a gq event, don't subtract the Z momentum from first q
fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p);
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity()
- partons[j].rapidity());
q -= partons[j+1].p;
}
if (param_.nlo.enabled){
double nlo_virtual =1.;
//Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) nlo_virtual+=sum;
return nlo_virtual;
}
return exp(sum);
}
std::vector<double> MatrixElement::virtual_corrections(
Event const & event,
const double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
std::vector<Particle> bosons = filter_AWZH_bosons(out);
if(event.jets().size() != param_.nlo.nj && param_.nlo.enabled) {
throw std::logic_error{
"Input event has number of jets different to stated NLO "
"input in config file: " + std::to_string(event.jets().size())
+ " vs " +std::to_string(param_.nlo.nj) + "\n"
};
}
if(bosons.size() > 2) {
throw not_implemented("Emission of >2 bosons is unsupported");
}
if(bosons.size() == 2) {
if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp) {
return virtual_corrections_WW(event, mur);
}
else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm) {
return virtual_corrections_WW(event, mur);
}
throw not_implemented("Emission of bosons of unsupported type");
}
if(bosons.size() == 1) {
const auto AWZH_boson = bosons[0];
if(std::abs(AWZH_boson.type) == pid::Wp){
return {virtual_corrections_W(event, mur, AWZH_boson)};
}
if(AWZH_boson.type == pid::Z_photon_mix){
if(is_gluon(in.back().type)){
// This is a qg event
return {virtual_corrections_Z_qg(event, mur, AWZH_boson, false)};
}
if(is_gluon(in.front().type)){
// This is a gq event
return {virtual_corrections_Z_qg(event, mur, AWZH_boson, true)};
}
// This is a qq event
return virtual_corrections_Z_qq(event, mur, AWZH_boson);
}
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
std::size_t first_idx = 0;
std::size_t last_idx = out.size() - 1;
// if there is a Higgs boson _not_ emitted off an incoming gluon,
// extremal qqbar or unordered gluon outside the extremal partons
// then it is not part of the FKL ladder
// and does not contribute to the virtual corrections
if((out.front().type == pid::Higgs && in.front().type != pid::gluon)
|| event.type() == event_type::unob
|| event.type() == event_type::qqbar_exb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs && in.back().type != pid::gluon)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf){
--last_idx;
}
std::size_t first_idx_qqbar = last_idx;
std::size_t last_idx_qqbar = last_idx;
//if central qqbar event, virtual correction do not occur between q-qbar.
if(event.type() == event_type::qqbar_mid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return {0.};
last_idx = std::distance(begin(out), backquark);
first_idx_qqbar = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(std::size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqbar) q -= out[last_idx+1].p;
for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| (out.back().type == pid::Higgs && in.back().type != pid::gluon)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf
);
if (param_.nlo.enabled){
const auto partons = filter_partons(event.outgoing());
double nlo_virtual =1.;
//Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) nlo_virtual+=exponent;
return {nlo_virtual};
}
return {std::exp(exponent)};
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return CL;
}
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, p1, p2);
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, p1, p2);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, p1, p2);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda){
return Cls;
}
return Cls - 4.0 / (kperp * kperp);
}
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
){
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return CL;
}
//! Lipatov vertex
double C2Lipatov( // B
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim,
CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom,
CLHEP::HepLorentzVector const & pop
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, pim, pip, pom, pop);
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, pim, pip, pom, pop);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, pim, pip, pom, pop);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, pa, pb, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda) {
return Cls;
}
return Cls - 4.0 / (kperp * kperp);
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pg Unordered gluon momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_uno_current(
[[maybe_unused]] ParticleID aptype,
ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
assert(aptype!=pid::gluon); // aptype cannot be gluon
const double t1 = (pa - p1 - pg).m2();
const double t2 = (pb - pn).m2();
return K_q * K(bptype, pn, pb)*currents::ME_unob_qq(pg, p1, pa, pn, pb) / (t1 * t2);
}
/** Matrix element squared for tree-level current-current scattering
* @param bptype Particle b PDG ID
* @param pgin Incoming gluon momentum
* @param p1 More backward (anti-)quark from splitting Momentum
* @param p2 Less backward (anti-)quark from splitting Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The forward qqbar contribution can be calculated by reversing the
* argument ordering.
*/
double ME_qqbar_current(
ParticleID bptype,
CLHEP::HepLorentzVector const & pgin,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb
){
const double t1 = (pgin - p1 - p2).m2();
const double t2 = (pn - pb).m2();
- return K(bptype, pn, pb)*currents::ME_qqbar_qg(pb, pgin, pn, p2, p1) / (t1 * t2);
+ return K_q * K(bptype, pn, pb)*currents::ME_qqbar_qg(pgin, pb, p1, p2, pn) / (t1 * t2);
}
/* \brief Matrix element squared for central qqbar tree-level current-current
* scattering
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqbarpair
* @param nbelow Number of gluons emitted after central qqbarpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering
*/
double ME_qqbar_mid_current(
ParticleID aptype, ParticleID bptype, int nabove,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons
){
using namespace currents;
// CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards)
const bool swap_qqbar=pqbar.rapidity() < pq.rapidity();
CLHEP::HepLorentzVector const & p1 = partons.front();
CLHEP::HepLorentzVector const & pn = partons.back();
const double t1 = (p1 - pa).m2();
const double t2 = (pb - pn).m2();
return K(aptype, p1, pa)
*K(bptype, pn, pb)
*ME_Cenqqbar_qq(
pa, pb, partons,
swap_qqbar, nabove
) / (t1 * t2);
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
const double t1 = (p1 - pa).m2();
const double t2 = (pb - pn).m2();
return K(aptype, p1, pa)
* K(bptype, pn, pb)
* currents::ME_qq(p1, pa, pn, pb)/(t1 * t2);
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc, ParticleProperties const & Wprop
){
using namespace currents;
assert(!(aptype==pid::gluon && bptype==pid::gluon));
if(aptype == pid::gluon || wc) {
// swap currents to ensure that the W is emitted off the first one
return ME_W_current(bptype, aptype, p1, pa, pn, pb, plbar, pl, false, Wprop);
}
// 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
const double current_contr = is_quark(aptype)?
ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop):
ME_W_qQ(p1,pl,plbar,pa,pn,pb,Wprop);
const double t1 = (pa - p1 - pl - plbar).m2();
const double tn = (pn - pb).m2();
return K(aptype, p1, pa)
* K(bptype, pn, pb)
* current_contr/(4.*(N_C*N_C - 1) * t1 * tn);
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_W_uno_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector plbar,
CLHEP::HepLorentzVector pl,
bool const wc, ParticleProperties const & Wprop
){
using namespace currents;
assert(bptype != pid::gluon || aptype != pid::gluon);
if(aptype == pid::gluon || wc) {
// emission off pb -- pn line
// 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(is_antiquark(bptype)) std::swap(plbar, pl);
const double t1 = (pa - p1 - pg).m2();
const double tn = (pb - pn - plbar - pl).m2();
return K_q*K_q*ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop)/(4.*(N_C*N_C - 1) * t1 * tn);
}
// emission off pa -- p1 line
if(is_antiquark(aptype)) std::swap(plbar, pl);
const double t1 = (pa - p1 - pg - plbar - pl).m2();
const double tn = (pb - pn).m2();
return K(bptype, pn, pb)/C_F*ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop)/(t1 * tn);
}
/** \brief Matrix element squared for backward qqbar tree-level current-current
* scattering With W+Jets
*
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqbarb Tree-Level Current-Current Scattering
*
* @note calculate forwards qqbar contribution by reversing argument ordering.
*/
double ME_W_qqbar_current(
ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector pq,
CLHEP::HepLorentzVector pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
if(is_anyquark(bptype) && wc) {
// W Must be emitted from forwards leg.
const double t1 = (pa - pq - pqbar).m2();
const double tn = (pb - pn - pl - plbar).m2();
return ME_W_Exqqbar_QQq(
pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop
)/(t1 * tn);
}
const double t1 = (pa - pl - plbar - pq - pqbar).m2();
const double tn = (pb - pn).m2();
return K(bptype, pn, pb)/C_F
* ME_WExqqbar_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop) / (t1 * tn);
}
/* \brief Matrix element squared for central qqbar tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqbarpair
* @param nbelow Number of gluons emitted after central qqbarpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqbar
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering
*/
double ME_W_qqbar_mid_current(
ParticleID aptype, ParticleID bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
// CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards)
const bool swap_qqbar=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;
if(wqq)
return wt*ME_WCenqqbar_qq(pa, pb, pl, plbar, partons,
is_antiquark(aptype),is_antiquark(bptype),
swap_qqbar, nabove, Wprop);
return wt*ME_W_Cenqqbar_qq(pa, pb, pl, plbar, partons,
is_antiquark(aptype), is_antiquark(bptype),
swap_qqbar, nabove, nbelow, wc, Wprop);
}
/** Matrix element squared for tree-level current-current scattering With Z+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
std::vector<double> ME_Z_current(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
const auto pZ = pl + plbar;
const double pref = K(aptype, p1, pa) * K(bptype, pn, pb)/(4.*(N_C*N_C-1));
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if(is_anyquark(aptype) && is_gluon(bptype)){
// This is a qg event
const double t1 = (pa-p1-pZ).m2();
const double tn = (pb-pn ).m2();
return { pref*ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw)/(t1 * tn) };
}
if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
const double t1 = (pa-p1 ).m2();
const double tn = (pb-pn-pZ).m2();
return { pref*ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw)/(t1 * tn) };
}
// This is a qq event
assert(is_anyquark(aptype) && is_anyquark(bptype));
const double t1_top = (pa-p1-pZ).m2();
const double t2_top = (pb-pn ).m2();
const double t1_bot = (pa-p1 ).m2();
const double t2_bot = (pb-pn-pZ).m2();
std::vector<double> res = ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
assert(res.size() == 3);
res[0] *= pref/(t1_top * t2_top);
res[1] *= pref/(t1_bot * t2_bot);
res[2] *= pref/sqrt(t1_top * t2_top * t1_bot * t2_bot);
return res;
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With Z+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
std::vector<double> ME_Z_uno_current(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
const auto pZ = pl + plbar;
const double pref = K(aptype, p1, pa)/C_F * K(bptype, pn, pb)/C_F;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if (is_anyquark(aptype) && is_gluon(bptype)) {
// This is a qg event
const double t1 = (pa-p1-pg-pZ).m2();
const double tn = (pb-pn ).m2();
return { pref*ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw)/(t1 * tn) };
}
if (is_gluon(aptype) && is_anyquark(bptype)) {
// This is a gq event
const double t1 = (pa-p1 ).m2();
const double tn = (pb-pn-pg-pZ).m2();
return { pref*ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw)/(t1 * tn) };
}
// This is a qq event
assert(is_anyquark(aptype) && is_anyquark(bptype));
const double t1_top = (pa-pg-p1-pZ).m2();
const double t2_top = (pb-pn ).m2();
const double t1_bot = (pa-pg-p1).m2();
const double t2_bot = (pb-pn-pZ).m2();
std::vector<double> res = ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
assert(res.size() == 3);
res[0] *= pref/(t1_top * t2_top);
res[1] *= pref/(t1_bot * t2_bot);
res[2] *= pref/sqrt(t1_top * t2_top * t1_bot * t2_bot);
return res;
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb, double vev
){
const double t1 = (pa - p1).m2();
const double tn = (pb - pn).m2();
return
K(aptype, p1, pa)
*K(bptype, pn, pb)
*currents::ME_H_qQ(
pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev
) / (t1 * tn * qH.m2() * qHp1.m2());
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param pg Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*
* @note This function assumes unordered gluon backwards from pa-p1 current.
* For unof, reverse call order
*/
double ME_Higgs_current_uno(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb, double vev
){
const double t1 = (pa - p1 - pg).m2();
const double tn = (pn - pb).m2();
return
K(aptype, p1, pa)
*K(bptype, pn, pb)
*currents::ME_H_unob_qQ(
pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev
) / (t1 * qH.m2() * qHp1.m2() * tn);
}
/** Matrix element squared for tree-level scattering with WW+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pl1bar Particle l1bar Momentum
* @param pl1 Particle l1 Momentum
* @param pl2bar Particle l2bar Momentum
* @param pl2 Particle l2 Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
std::vector <double> ME_WW_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pl1bar,
CLHEP::HepLorentzVector const & pl1,
CLHEP::HepLorentzVector const & pl2bar,
CLHEP::HepLorentzVector const & pl2,
ParticleProperties const & Wprop
){
using namespace currents;
if (aptype > 0 && bptype > 0)
return ME_WW_qQ(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
if (aptype < 0 && bptype > 0)
return ME_WW_qbarQ(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
if (aptype > 0 && bptype < 0)
return ME_WW_qQbar(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
if (aptype < 0 && bptype < 0)
return ME_WW_qbarQbar(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
throw std::logic_error("unreachable");
}
void validate(MatrixElementConfig const & config) {
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
std::vector<double> MatrixElement::tree_kin(
Event const & ev
) const {
if(!ev.valid_hej_state(param_.soft_pt_regulator)) return {0.};
if(!ev.valid_incoming()){
throw std::invalid_argument{
"Invalid momentum for one or more incoming particles. "
"Incoming momenta must have vanishing mass and transverse momentum."
};
}
std::vector<Particle> bosons = filter_AWZH_bosons(ev.outgoing());
if(bosons.empty()) {
return {tree_kin_jets(ev)};
}
if(bosons.size() == 1) {
switch(bosons[0].type){
case pid::Higgs:
return {tree_kin_Higgs(ev)};
case pid::Wp:
case pid::Wm:
return {tree_kin_W(ev)};
case pid::Z_photon_mix:
return tree_kin_Z(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
if(bosons.size() == 2) {
if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp){
return tree_kin_WW(ev);
}
else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm){
return tree_kin_WW(ev);
}
throw not_implemented("Emission of bosons of unsupported type");
}
throw not_implemented("Emission of >2 bosons is unsupported");
}
namespace {
constexpr int EXTREMAL_JET_IDX = 1;
constexpr int NO_EXTREMAL_JET_IDX = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == EXTREMAL_JET_IDX;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
template<class InputIterator>
std::vector <double> FKL_ladder_weight_mix(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0_t, CLHEP::HepLorentzVector const & q0_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
const double lambda
){
double wt_top = 1;
double wt_bot = 1;
double wt_mix = 1;
auto qi_t = q0_t;
auto qi_b = q0_b;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1_t = qi_t - g;
const auto qip1_b = qi_b - g;
if(treat_as_extremal(*gluon_it)){
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, lambda)*C_A;
} else{
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, p1, pn, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
}
qi_t = qip1_t;
qi_b = qip1_b;
}
return {wt_top, wt_bot, wt_mix};
}
std::vector<Particle> tag_extremal_jet_partons( Event const & ev ){
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(NO_EXTREMAL_JET_IDX);
}
return out_partons;
}
auto const & jets = ev.jets();
std::vector<fastjet::PseudoJet> extremal_jets;
if(! is_backward_g_to_h(ev)) {
auto most_backward = begin(jets);
// skip jets caused by unordered emission or qqbar
if(ev.type() == event_type::unob || ev.type() == event_type::qqbar_exb){
assert(jets.size() >= 2);
++most_backward;
}
extremal_jets.emplace_back(*most_backward);
}
if(! is_forward_g_to_h(ev)) {
auto most_forward = end(jets) - 1;
if(ev.type() == event_type::unof || ev.type() == event_type::qqbar_exf){
assert(jets.size() >= 2);
--most_forward;
}
extremal_jets.emplace_back(*most_forward);
}
const auto extremal_jet_indices = ev.particle_jet_indices(
extremal_jets
);
assert(extremal_jet_indices.size() == out_partons.size());
for(std::size_t i = 0; i < out_partons.size(); ++i){
assert(is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
EXTREMAL_JET_IDX:
NO_EXTREMAL_JET_IDX;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
double tree_kin_jets_qqbar_mid(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
double lambda
){
CLHEP::HepLorentzVector pq;
CLHEP::HepLorentzVector pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqbar
auto qqbart = q0;
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqbart -= to_HepLorentzVector(*parton_it);
}
const int nabove = std::distance(begin_ladder, backmidquark);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_qqbar_mid_current(
aptype, bptype, nabove, pa, pb,
pq, pqbar, partonsHLV
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqbart, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_qqbar(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const auto pgin = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto p1 = to_HepLorentzVector(*BeginPart);
const auto p2 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon.
const double current_factor = ME_qqbar_current(
(EndIn-1)->type, pgin, p1, p2, pn, pb
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pgin-p1-p2, pgin, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const double current_factor = ME_uno_current(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pa-p1-pg, pa, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
} // namespace
double MatrixElement::tree_kin_jets(Event const & ev) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if (ev.type()==event_type::FKL){
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn,
param_.regulator_lambda
);
}
if (ev.type()==event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
if (ev.type()==event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
if (ev.type()==event_type::extremal_qqbar_backward){
return tree_kin_jets_qqbar(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
if (ev.type()==event_type::extremal_qqbar_forward){
return tree_kin_jets_qqbar(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
if (ev.type()==event_type::central_qqbar){
return tree_kin_jets_qqbar_mid(incoming[0].type, incoming[1].type,
to_HepLorentzVector(incoming[0]),
to_HepLorentzVector(incoming[1]),
partons, param_.regulator_lambda);
}
throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets");
}
namespace {
double tree_kin_W_FKL(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
double lambda, ParticleProperties const & Wprop
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_uno(InIter BeginIn, partIter BeginPart,
partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w
auto q0 = pa - p1 - pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_uno_current(
(BeginIn)->type, (BeginIn+1)->type, pn, pb,
p1, pa, pg, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_qqbar(InIter BeginIn, partIter BeginPart,
partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const bool swap_qqbar=is_quark(*BeginPart);
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_qqbar?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_qqbar?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqbar_current(
(BeginIn+1)->type, pa, pb,
pq, pqbar, pn, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
double tree_kin_W_qqbar_mid(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
double lambda, ParticleProperties const & Wprop
){
CLHEP::HepLorentzVector pq;
CLHEP::HepLorentzVector pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons.front());
auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
// t-channel momentum after qqbar
auto qqbart = q0;
bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqbar emit W
bool wc = !wqq && (aptype==partons.front().type); //leg b emits w
assert(!wqq || !wc);
if(wqq){ // emission from qqbar
qqbart -= pl + plbar;
} else if(!wc) { // emission from leg a
q0 -= pl + plbar;
qqbart -= pl + plbar;
}
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqbart -= to_HepLorentzVector(*parton_it);
}
const int nabove = std::distance(begin_ladder, backmidquark);
const int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqbar_mid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqbart, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
} // namespace
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
#ifndef NDEBUG
// assert that there is exactly one decay corresponding to the W
assert(ev.decays().size() == 1);
auto const & w_boson{
std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(),
[] (Particle const & p) -> bool {
return std::abs(p.type) == ParticleID::Wp;
}) };
assert(w_boson != ev.outgoing().cend());
assert( static_cast<long int>(ev.decays().cbegin()->first)
== std::distance(ev.outgoing().cbegin(), w_boson) );
#endif
// find decay products of W
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) )
|| ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) );
// get lepton & neutrino
CLHEP::HepLorentzVector plbar;
CLHEP::HepLorentzVector pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.at(1));
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == FKL){
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_backward){
return tree_kin_W_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_forward){
return tree_kin_W_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqbar_backward){
return tree_kin_W_qqbar(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqbar_forward){
return tree_kin_W_qqbar(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
assert(ev.type() == central_qqbar);
return tree_kin_W_qqbar_mid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
namespace /* WW */ {
std::vector <double> tree_kin_WW_FKL(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & pl1bar, CLHEP::HepLorentzVector const & pl1,
CLHEP::HepLorentzVector const & pl2bar, CLHEP::HepLorentzVector const & pl2,
double lambda, ParticleProperties const & Wprop
){
assert(is_anyquark(aptype));
assert(is_anyquark(bptype));
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const std::vector <double> current_factor = ME_WW_current(
aptype, bptype, pn, pb, p1, pa,
pl1bar, pl1, pl2bar, pl2,
Wprop
);
auto const begin_ladder = cbegin(partons) + 1;
auto const end_ladder = cend(partons) - 1;
// pa -> W1 p1, pb -> W2 + pn
const auto q0 = pa - p1 - (pl1 + pl1bar);
// pa -> W2 p1, pb -> W1 + pn
const auto q1 = pa - p1 - (pl2 + pl2bar);
const std::vector <double> ladder_factor = FKL_ladder_weight_mix(
begin_ladder, end_ladder,
q0, q1, pa, pb, p1, pn,
lambda
);
assert(current_factor.size() == 3);
assert(current_factor.size() == ladder_factor.size());
std::vector<double> result(current_factor.size());
for(size_t i=0; i<current_factor.size(); ++i){
result[i] = K_q*K_q/(4.*(N_C*N_C - 1.))
*current_factor.at(i)
*ladder_factor.at(i);
}
const double t1_top = q0.m2();
const double t2_top = (pb-pn-pl2bar-pl2).m2();
const double t1_bot = q1.m2();
const double t2_bot = (pb-pn-pl1bar-pl1).m2();
result[0] /= t1_top * t2_top;
result[1] /= t1_bot * t2_bot;
result[2] /= sqrt(t1_top * t2_top * t1_bot * t2_bot);
return result;
}
} // namespace
std::vector <double> MatrixElement::tree_kin_WW(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const pa = to_HepLorentzVector(incoming[0]);
auto const pb = to_HepLorentzVector(incoming[1]);
auto const partons = tag_extremal_jet_partons(ev);
// W1 & W2
assert(ev.decays().size() == 2);
std::vector<CLHEP::HepLorentzVector> plbar;
std::vector<CLHEP::HepLorentzVector> pl;
for(auto const & decay_pair : ev.decays()) {
auto const decay = decay_pair.second;
// TODO: how to label W1, W2
if(decay.at(0).type < 0) {
plbar.emplace_back(to_HepLorentzVector(decay.at(0)));
pl .emplace_back(to_HepLorentzVector(decay.at(1)));
}
else {
pl .emplace_back(to_HepLorentzVector(decay.at(0)));
plbar.emplace_back(to_HepLorentzVector(decay.at(1)));
}
}
if(ev.type() == FKL) {
return tree_kin_WW_FKL(
incoming[0].type, incoming[1].type,
pa, pb, partons,
plbar[0], pl[0], plbar[1], pl[1],
param_.regulator_lambda,
param_.ew_parameters.Wprop()
);
}
throw std::logic_error("Can only reweight FKL events in WW");
}
namespace{
std::vector <double> tree_kin_Z_FKL(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw
){
const auto p1 = to_HepLorentzVector(partons[0]);
const auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
const std::vector <double> current_factor = ME_Z_current(
aptype, bptype, pn, pb, p1, pa,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
} else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-p1;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
} else {
// This is a qq event
const auto q0 = pa-p1-plbar-pl;
const auto q1 = pa-p1;
ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
template<class InIter, class partIter>
std::vector <double> tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const ParticleID aptype = (BeginIn)->type;
const ParticleID bptype = (BeginIn+1)->type;
const std::vector <double> current_factor = ME_Z_uno_current(
aptype, bptype, pn, pb, p1, pa, pg,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-pg-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-pg-p1;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q0 = pa-pg-p1-plbar-pl;
const auto q1 = pa-pg-p1;
ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
} // namespace
std::vector<double> MatrixElement::tree_kin_Z(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
// find decay products of Z
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0))
&& decay.at(0).type==-decay.at(1).type);
// get leptons
CLHEP::HepLorentzVector plbar;
CLHEP::HepLorentzVector pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.at(1));
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
const double stw2 = param_.ew_parameters.sin2_tw();
const double ctw = param_.ew_parameters.cos_tw();
if(ev.type() == FKL){
return tree_kin_Z_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_backward){
return tree_kin_Z_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_forward){
return tree_kin_Z_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets");
}
double MatrixElement::tree_kin_Higgs(Event const & ev) const {
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
// kinetic matrix element square for backward Higgs emission
// cf. eq:ME_h_jets_peripheral in developer manual,
// but without factors \alpha_s and the FKL ladder
double MatrixElement::MH2_backwardH(
const ParticleID type_forward,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pH,
CLHEP::HepLorentzVector const & pn
) const {
using namespace currents;
const double vev = param_.ew_parameters.vev();
return K(type_forward, pn, pb)*ME_H_gq(
pH, pa, pn, pb,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(4.*(N_C*N_C - 1)*(pa-pH).m2()*(pb-pn).m2());
}
// kinetic matrix element square for backward unordered emission
// and forward g -> Higgs
double MatrixElement::MH2_unob_forwardH(
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pH
) const {
using namespace currents;
const double vev = param_.ew_parameters.vev();
constexpr double K_f1 = K_q;
constexpr double nhel = 4.;
return K_f1*ME_juno_jgH(
pg, p1, pa, pH, pb,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(nhel*(N_C*N_C - 1)*(pa - p1 - pg).m2()*(pb - pH).m2());
}
double MatrixElement::tree_kin_Higgs_first(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(is_anyquark(incoming.front())) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming.front());
const auto pb = to_HepLorentzVector(incoming.back());
const auto pH = to_HepLorentzVector(outgoing.front());
const auto end_ladder = end(partons) - ((ev.type() == event_type::unof)?2:1);
const auto pn = to_HepLorentzVector(*end_ladder);
const double ladder = FKL_ladder_weight(
begin(partons), end_ladder,
pa - pH, pa, pb, pa, pn,
param_.regulator_lambda
);
if(ev.type() == event_type::unof) {
const auto pg = to_HepLorentzVector(outgoing.back());
return MH2_unob_forwardH(
pb, pa, pg, pn, pH
)*ladder;
}
return MH2_backwardH(
incoming.back().type,
pa, pb, pH, pn
)*ladder;
}
double MatrixElement::tree_kin_Higgs_last(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(is_anyquark(incoming.back())) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming.front());
const auto pb = to_HepLorentzVector(incoming.back());
const auto pH = to_HepLorentzVector(outgoing.back());
auto begin_ladder = begin(partons) + 1;
auto q0 = pa - to_HepLorentzVector(partons.front());
if(ev.type() == event_type::unob) {
q0 -= to_HepLorentzVector(*begin_ladder);
++begin_ladder;
}
const auto p1 = to_HepLorentzVector(*(begin_ladder - 1));
const double ladder = FKL_ladder_weight(
begin_ladder, end(partons),
q0, pa, pb, p1, pb,
param_.regulator_lambda
);
if(ev.type() == event_type::unob) {
const auto pg = to_HepLorentzVector(outgoing.front());
return MH2_unob_forwardH(
pa, pb, pg, p1, pH
)*ladder;
}
return MH2_backwardH(
incoming.front().type,
pb, pa, pH, p1
)*ladder;
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart,
CLHEP::HepLorentzVector const & qH,
CLHEP::HepLorentzVector const & qHp1,
double mt, bool inc_bot, double mb, double vev
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb, vev
);
}
} // namespace
double MatrixElement::tree_kin_Higgs_between(Event const & ev) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor = NAN;
if(ev.type() == FKL){
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
}
else if(ev.type() == unob){
current_factor = tree_kin_Higgs_uno(
begin(incoming), end(incoming), begin(partons),
end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = tree_kin_Higgs_uno(
rbegin(incoming), rend(incoming), rbegin(partons),
rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
throw std::logic_error("Can only reweight FKL or uno processes in H+Jets");
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor/(4.*(N_C*N_C-1.))*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) {
std::vector<Particle> bosons = filter_AWZH_bosons(ev.outgoing());
if(bosons.empty()) {
return 1.;
}
if(bosons.size() == 1) {
switch(bosons[0].type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return alpha_w*alpha_w;
case pid::Z_photon_mix:
return alpha_w*alpha_w;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
if(bosons.size() == 2) {
if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp) {
return alpha_w*alpha_w*alpha_w*alpha_w;
}
else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm) {
return alpha_w*alpha_w*alpha_w*alpha_w;
}
throw not_implemented("Emission of bosons of unsupported type");
}
throw not_implemented("Emission of >2 bosons is unsupported");
}
} // namespace
double MatrixElement::tree_param(Event const & ev, double mur) const {
assert(is_resummable(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*BETA0*std::log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ
diff --git a/src/jets.cc b/src/jets.cc
index b006d8b..5ba1f8b 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,230 +1,222 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2022
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include "HEJ/Constants.hh"
// generated headers
#include "HEJ/currents/j_j.hh"
-#include "HEJ/currents/j_jqqbar.hh"
+#include "HEJ/currents/jqqbar_j.hh"
#include "HEJ/currents/j_qqbar_j.hh"
#include "HEJ/currents/juno_j.hh"
namespace {
using COM = std::complex<double>;
} // Anonymous Namespace
namespace HEJ {
namespace currents {
double ME_qq(
HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
using helicity::plus;
using helicity::minus;
// we only have to distinguish between the two possibilities of
// contracting same-helicity or opposite-helicity currents.
COM const Mp = HEJ::j_j<minus, plus>(p1in, p1out, p2in, p2out);
COM const Mm = HEJ::j_j<minus, minus>(p1in, p1out, p2in, p2out);
double const sst = std::norm(Mm) + std::norm(Mp);
// Multiply by factor 2 (helicities)
return 2.*sst;
}
namespace {
template<Helicity h1, Helicity hg>
double amp_juno_j(
HLV const & pa, HLV const & pb,
HLV const & pg, HLV const & p1, HLV const & p2
) {
// TODO: code duplication with Wjets
const COM u1 = U1_j<h1, hg>(pa,p1,pb,p2,pg);
const COM u2 = U2_j<h1, hg>(pa,p1,pb,p2,pg);
const COM l = L_j<h1, hg>(pa,p1,pb,p2,pg);
const COM x = u1 - l;
const COM y = u2 + l;
return C_F*norm(x + y) - C_A*(x*conj(y)).real();
}
} // Anonymous Namespace
double ME_unob_qq(
HLV const & pg,
HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
using helicity::minus;
using helicity::plus;
// only calculate half of the helicity amplitudes,
// the remaining ones follow from CP symmetry
const double amm = amp_juno_j<minus, minus>(p1in, p2in, pg, p1out, p2out);
const double amp = amp_juno_j<minus, plus>(p1in, p2in, pg, p1out, p2out);
const double apm = amp_juno_j< plus, minus>(p1in, p2in, pg, p1out, p2out);
const double app = amp_juno_j< plus, plus>(p1in, p2in, pg, p1out, p2out);
constexpr double hel_fac = 2.;
return hel_fac * (amm + amp + apm + app);
}
namespace {
- // helicity amplitudes for j jqqbar contraction
+ // helicity amplitudes for jqqbar j contraction
template<Helicity h1, Helicity hg>
- double amp_j_jqqbar(
+ double amp_jqqbar_j(
HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2, HLV const & p3
) {
- // TODO: code duplication with Wjets.cc
- // TODO: symbolic colour factors
+ const COM u1 = U1X_j<h1, hg>(pa,pb,p1,p2,p3);
+ const COM u2 = U2X_j<h1, hg>(pa,pb,p1,p2,p3);
+ const COM l = LX_j <h1, hg>(pa,pb,p1,p2,p3);
- using std::norm;
-
- static constexpr double cm1m1 = 8./3.;
- static constexpr double cm2m2 = 8./3.;
- static constexpr double cm3m3 = 6.;
- static constexpr double cm1m2 =-1./3.;
- static constexpr COM cm1m3 = COM{0.,-3.};
- static constexpr COM cm2m3 = COM{0.,3.};
-
- const COM m1 = j_qggm1<h1,hg>(pb,p2,p3,pa,p1);
- const COM m2 = j_qggm2<h1,hg>(pb,p2,p3,pa,p1);
- const COM m3 = j_qggm3<h1,hg>(pb,p2,p3,pa,p1);
+ const COM x = u1 - l;
+ const COM y = u2 + l;
- return
- + cm1m1*norm(m1)
- + cm2m2*norm(m2)
- + cm3m3*norm(m3)
- + 2.*real(cm1m2*m1*conj(m2))
- + 2.*real(cm1m3*m1*conj(m3))
- + 2.*real(cm2m3*m2*conj(m3))
- ;
+ return C_F*norm(x + y) - C_A*(x*conj(y)).real();
}
} // Anonymous Namespace
double ME_qqbar_qg(
HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2, HLV const & p3
){
using helicity::minus;
using helicity::plus;
- const double Mmmm = amp_j_jqqbar<minus, minus>(pa, pb, p1, p2, p3);
- const double Mmmp = amp_j_jqqbar<minus, plus>(pa, pb, p1, p2, p3);
- const double Mpmm = amp_j_jqqbar< plus, minus>(pa, pb, p1, p2, p3);
- const double Mpmp = amp_j_jqqbar< plus, plus>(pa, pb, p1, p2, p3);
+ // only calculate half of the helicity amplitudes,
+ // the remaining ones follow from CP symmetry
+ const double ampsq =
+ + amp_jqqbar_j<minus, minus>(pa, pb, p1, p2, p3)
+ + amp_jqqbar_j<minus, plus> (pa, pb, p1, p2, p3)
+ + amp_jqqbar_j<plus, minus>(pa, pb, p1, p2, p3)
+ + amp_jqqbar_j<plus, plus> (pa, pb, p1, p2, p3)
+ ;
+
+ constexpr double hel_fac = 2.;
+
+ // correct colour averaging after crossing:
+ constexpr double avg_fac = N_C / (N_C*N_C - 1.);
- // Factor of 2 for the helicity for conjugate configurations
- return 2.*(Mmmm+Mmmp+Mpmm+Mpmp)/3./C_F;
+ return avg_fac * hel_fac * ampsq;
}
namespace {
// helicity amplitudes for matrix element with central qqbar
// TODO: update to arxiv:2012.10310, around eq. (3.33)
template<Helicity h1a, Helicity h4b>
double amp_Cenqqbar_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p4,
HLV const & pq, HLV const & pqbar,
HLV const & q11, HLV const & q12
){
using std::norm;
const COM sym = M_sym<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, q11, q12
);
const COM cross = M_cross<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, q11, q12
);
const COM uncross = M_qbar<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, q11, q12
);
// Colour factors
static constexpr double cmsms = 3.;
static constexpr double cmumu = 4./3.;
static constexpr double cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0., 3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
static constexpr double cmumc = -1./6.;
return
+cmsms*norm(sym)
+cmumu*norm(uncross)
+cmcmc*norm(cross)
+2.*real(cmsmu*sym*conj(uncross))
+2.*real(cmsmc*sym*conj(cross))
+2.*real(cmumc*uncross*conj(cross))
;
}
} // Anonymous Namespace
double ME_Cenqqbar_qq(
HLV const & pa, HLV const & pb,
std::vector<HLV> const & partons,
const bool qqbarmarker, const std::size_t nabove
){
using helicity::plus;
using helicity::minus;
CLHEP::HepLorentzVector q1 = pa;
for(std::size_t i = 0; i <= nabove; ++i){
q1 -= partons[i];
}
const auto qq = 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 = partons[nabove+1];
HLV pqbar = partons[nabove+2];
if(qqbarmarker) std::swap(pq, pqbar);
const HLV p1 = partons.front();
const HLV pn = partons.back();
// 8 helicity choices, but only 4 independent ones
//(complex conjugation related).
// In principle, the proper helicity labeling depends on
// whether we have antiquark lines (aqlinea and aqlineb).
// However, this only corresponds to a relabeling,
// which doesn't change the sum over all helicities below.
const double amp_mm = amp_Cenqqbar_qq<minus, minus>(
pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
);
const double amp_mp = amp_Cenqqbar_qq<minus, plus>(
pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
);
const double amp_pm = amp_Cenqqbar_qq<plus, minus>(
pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
);
const double amp_pp = amp_Cenqqbar_qq<plus, plus>(
pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
);
//Result (averaged, without coupling or t-channel props). Factor of
//2 for the 4 helicity configurations I didn't work out explicitly
const HLV q3 = q1 - pq - pqbar;
constexpr double hel_fac = 2.;
// TODO: explain the 1/64
return hel_fac/64.*(amp_mm+amp_mp+amp_pm+amp_pp) / (q1.m2()*q3.m2());
}
} // namespace currents
} // namespace HEJ
diff --git a/t/ME_data/ME_jets_subl_tree.dat b/t/ME_data/ME_jets_subl_tree.dat
index 9661e27..a672850 100644
--- a/t/ME_data/ME_jets_subl_tree.dat
+++ b/t/ME_data/ME_jets_subl_tree.dat
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:5dd7a8c9890012143e74f5da58cb697d0d0eee571d0fff6421cd37b251734a8f
+oid sha256:ae98eb4606999f889461cca3807664a6f8afe85c2c1ae4a99efc92759c101abd
size 126752
diff --git a/t/ME_data/ME_jets_subl_tree_new.dat b/t/ME_data/ME_jets_subl_tree_new.dat
index 81332ba..50b7aaa 100644
--- a/t/ME_data/ME_jets_subl_tree_new.dat
+++ b/t/ME_data/ME_jets_subl_tree_new.dat
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:a0459af6f2937220b6872e827fa69dfadaaef5551f19346948eee9e8927f52fb
+oid sha256:08550c187d24baf3e2dc0a50f9c0cccc0449abe34d21aa642726373f238b952e
size 52080

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 10:35 AM (1 d, 11 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5093171
Default Alt Text
(143 KB)

Event Timeline