Page MenuHomeHEPForge

No OneTemporary

diff --git a/Changes.md b/Changes.md
index 6b1c3a3..b417c6b 100644
--- a/Changes.md
+++ b/Changes.md
@@ -1,126 +1,129 @@
# Changelog
This is the log for changes to the HEJ program. Further changes to the HEJ API
are documented in [`Changes-API.md`](Changes-API.md). If you are using HEJ as a
library, please also read the changes there.
## Version 2.1
Implemented W boson with jets for complete set of first subleading processes
(unordered gluon, central and extremal quark-antiquark pair), see
[arXiv:TODO](https://arxiv.org/abs/TODO). Ported unordered gluon emissions for
pure jets from HEJ 1. This release include many changes to the code, which
affect users of HEJ as a library (see [`Changes-API.md`](Changes-API.md)).
### 2.1.0
#### New Processes
* Resummation for W bosons with jets
- New subleading processes `extremal qqx` & `central qqx` for a quark and
anti-quark in the final state, e.g. `g g => u d_bar Wm g` (the other
subleading processes also work with W's)
- `HEJFOG` can generate multiple jets together with a (off-shell) W bosons
decaying into lepton & neutrino
* Resummation can now be performed on all subleading processes within pure
jets also. This includes `unordered`, `extremal qqbar` and `central
qqbar` processes.
#### More Physics implementation
* Partons now have a Colour charge
- Colours are read from and written to LHE files
- For reweighted events the colours are created according to leading colour in
the FKL limit
* Use relative fraction for soft transverse momentum in extremal jets (`max ext
soft pt fraction`). This supersedes `min extparton pt`, which is now
deprecated and will be removed in future versions.
#### Updates to Runcard
* Allow multiplication and division of multiple scale functions e.g.
`H_T/2*m_j1j2`
* Grouped `event treatment` for subleading channels together in runcard
- Rename `non-HEJ` processes to `non-resummable`
* Read electro-weak constants from input
- new mandatory setting `vev` to change vacuum expectation value
- new mandatory settings `particle properties` to specify mass & width of
bosons
- FOG: decays are now specified in `decays` setting (previously under
`particle properties`)
* Allow loading multiple analyses with `analyses`. The old `analysis` (with "i")
is marked deprecated.
* Optional setting to specify maximal number of Fixed Order events (`max
events`, default is all)
* Allow changing the regulator lambda in input (`regulator parameter`, only for
advanced users)
#### Changes to Input/Output
* Added support to read & write `hdf5` event files suggested in
[arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs
[HighFive](https://github.com/BlueBrain/HighFive))
* Support input with average weight equal to the cross section (`IDWTUP=1 or 4`)
* Support unmodified Les Houches Event Files written by Sherpa with
`cross section = sum(weights)/sum(trials)`
* Analyses now get general run information (`LHEF::HEPRUP`) in the
constructor. **This might break previously written, external analyses!**
- external analyses should now be created with
`make_analysis(YAML::Node const & config, LHEF::HEPRUP const & heprup)`
* Support `rivet` version 3 with both `HepMC` version 2 and 3
- Multiple weights with `rivet 3` will only create one `.yoda` file (instead
of one per weight/scale)
* Added option to unweight only resummation events
(`unweight: {type: resummation}`)
* Added option for partially unweighting resummation events, similar to
the fixed-order generator.
* Improved unweighting algorithm.
* Follow HepMC convention for particle Status codes: incoming = 11,
decaying = 2, outgoing = 1 (unchanged)
#### Miscellaneous
+* Currents are now generated with [`FORM`](https://github.com/vermaseren/form)
+ - `FORM` is included as a `git submodule`, use `git submodule update --init`
+ to download `FORM`
* Print cross sections at end of run
* Added example analysis & scale to `examples/`. Everything in `examples/` will
be build when the flag `-DBUILD_EXAMPLES=TRUE` is set in `cmake`.
* Use `git-lfs` for raw data in test (`make test` now requires `git-lfs`)
* Dropped support for HepMC 3.0.0, either HepMC version 2 or >3.1 is required
- It is now possible to write out both HepMC 2 and HepMC 3 events at the same
time
* Create [Sphinx](http://sphinx-doc.org/) and [Doxygen](http://doxygen.org/)
documentation by `make sphinx` or `make doxygen` in your `build/` folder
## Version 2.0
First release of HEJ 2. Complete code rewrite compared to HEJ 1. Improved
matching to Fixed Order ([arXiv:1805.04446](https://arxiv.org/abs/1805.04446)).
Implemented processes: Higgs boson with jets (FKL and unordered gluon emission,
with finite quark mass loop,
[arXiv:1812.08072](https://arxiv.org/abs/1812.08072)), and pure jets (only FKL).
See [arXiv:1902.08430](https://arxiv.org/abs/1902.08430)
### 2.0.6
* Fixed compiling rivet when YODA headers are _outside_ of rivet directory
### 2.0.5
* Fixed event classification for input not ordered in rapidity
### 2.0.4
* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
### 2.0.3
* Fixed parsing of (numerical factor) * (base scale) in configuration
* Don't change scale names, but sanitise Rivet output file names instead
### 2.0.2
* Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
`"/"` and `"*"` before)
### 2.0.1
* Fixed name of fixed-order generator in error message.
### 2.0.0
* First release
diff --git a/current_generator/W_central_qqbar.frm b/current_generator/W_central_qqbar.frm
new file mode 100644
index 0000000..19fe09b
--- /dev/null
+++ b/current_generator/W_central_qqbar.frm
@@ -0,0 +1,116 @@
+*/**
+* \brief Calculation of W + jets with central qqbar emission, W emitted from extremal leg
+*
+* TODO: unify conventions with developer manual
+* the current dictionary is as follows:
+*
+* code | manual
+* pq | p_2
+* pqbar | p_3
+* pl | p_A
+* plbar | p_B
+*
+* \authors The HEJ collaboration (see AUTHORS for details)
+* \date 2020
+* \copyright GPLv2 or later
+*/
+#include- include/helspin.frm
+#include- include/write.frm
+
+i mu1,...,mu20;
+i mu,nu,rho,sigma;
+cf T1a, T4b, Xunc, Xcro, X1a, X4b, V3g, JV;
+v p1,pa,p4,pb,pq,pqbar,pl,plbar,q1,q3,q11,q12;
+s h;
+s tunc1,tunc2,tcro1,tcro2,s2AB,s3AB,s23AB;
+s t1,t3,s12,s13,s1A,s1B,sa2,sa3,saA,saB,s42,s43,s4A,s4B,sb2,sb3,sbA,sbB;
+cf m2;
+
+#do h1={+,-}
+ #do h2={+,-}
+ l [M_uncross_W `h1'`h2'] = T1a(`h1'1, mu)*T4b(`h2'1, nu)*Xunc(mu, nu);
+ l [M_cross_W `h1'`h2'] = T1a(`h1'1, mu)*T4b(`h2'1, nu)*Xcro(mu, nu);
+ l [M_sym_W `h1'`h2'] = T1a(`h1'1, mu)*T4b(`h2'1, nu)*(
+ X1a(mu, nu, sigma) - X4b(mu, nu, sigma) + V3g(mu, nu, sigma)
+ )/s23AB * JV(sigma);
+ #enddo
+#enddo
+
+id T1a(h?, mu?) = Current(h, p1, mu, pa);
+id T4b(h?, mu?) = Current(h, p4, mu, pb);
+* see eq:X_Unc in developer manual, multiplied by -t1*t3
+id Xunc(mu, nu) = -Current(-1, pl, sigma, plbar)*(
+ - Current(-1, pq, sigma, pq + pl + plbar, mu, q3 + pqbar, nu, pqbar)/s2AB/tunc2
+ + Current(-1, pq, mu, q1 - pq, sigma, q3 + pqbar, nu, pqbar)/tunc1/tunc2
+ + Current(-1, pq, mu, q1 - pq, nu, pqbar + pl + plbar, sigma, pqbar)/tunc1/s3AB
+);
+* see eq:X_Cro in developer manual, multiplied by -t1*t3
+id Xcro(mu, nu) = -Current(-1, pl, sigma, plbar)*(
+ - Current(-1, pq, nu, q3 + pq, mu, pqbar + pl + plbar, sigma, pqbar)/tcro1/s3AB
+ + Current(-1, pq, nu, q3 + pq, sigma, q1 - pqbar, mu, pqbar)/tcro1/tcro2
+ + Current(-1, pq, sigma, pq + pl + plbar, nu, q1 - pqbar, mu, pqbar)/s2AB/tcro2
+);
+* see eq:X_1a in developer manual, multiplied by -t1*t3
+id X1a(mu, nu, sigma) = t1*d_(mu, nu)*(
+ + p1(sigma)/(s12 + s13 + s1A + s1B)
+ + pa(sigma)/(sa2 + sa3 + saA + saB)
+);
+* see eq:X_4b in developer manual, multiplied by t1*t3
+id X4b(mu, nu, sigma) = t3*d_(mu, nu)*(
+ + p4(sigma)/(s42 + s43 + s4A + s4B)
+ + pb(sigma)/(sb2 + sb3 + sbA + sbB)
+);
+* see eq:3GluonWEmit in developer manual, multiplied by t1*s23AB*t3
+id V3g(mu, nu, sigma) = (
+ + (q1(nu) + pq(nu) + pqbar(nu) + pl(nu) + plbar(nu))*d_(mu, sigma)
+ + (q3(mu) - pq(mu) - pqbar(mu) - pl(mu) - plbar(mu))*d_(nu, sigma)
+ - (q1(sigma) + q3(sigma))*d_(mu, nu)
+);
+* see eq:J_V in developer manual
+id JV(sigma?) = -i_*(
+ + Current(-1, pq, rho, pq + pl + plbar, sigma, pqbar)/s2AB
+ - Current(-1, pq, sigma, pqbar + pl + plbar, rho, pqbar)/s3AB
+)*Current(-1, pl, rho, plbar);
+
+multiply replace_(
+ s2AB, m2(pl+plbar+pq),
+ s3AB, m2(pl+plbar+pqbar),
+ s23AB, m2(pl+plbar+pq+pqbar),
+ tunc1, m2(q1-pq),
+ tunc2, m2(q3+pqbar),
+ tcro1, m2(q3+pq),
+ tcro2, m2(q1-pqbar),
+ t1, m2(q1),
+ t3, m2(q3),
+ s12, m2(p1+pq),
+ s13, m2(p1+pqbar),
+ s1A, m2(p1+pl),
+ s1B, m2(p1+plbar),
+ sa2, m2(pa+pq),
+ sa3, m2(pa+pqbar),
+ saA, m2(pa+pl),
+ saB, m2(pa+plbar),
+ s42, m2(p4+pq),
+ s43, m2(p4+pqbar),
+ s4A, m2(p4+pl),
+ s4B, m2(p4+plbar),
+ sb2, m2(pb+pq),
+ sb3, m2(pb+pqbar),
+ sbA, m2(pb+pl),
+ sbB, m2(pb+plbar)
+);
+multiply replace_(q3, q1-pl-plbar-pq-pqbar);
+* replace q1 by sum of lightlike momenta
+* @TODO: choose q11 as collinear to some other vector
+multiply replace_(q1, q11-q12);
+.sort
+#call ContractCurrents()
+.sort
+format O4;
+format c;
+#call WriteHeader(`OUTPUT')
+#call WriteOptimised(`OUTPUT',M_uncross_W,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
+#call WriteOptimised(`OUTPUT',M_cross_W,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
+#call WriteOptimised(`OUTPUT',M_sym_W,2,pa,p1,pb,p4,pq,pqbar,pl,plbar,q11,q12)
+#call WriteFooter(`OUTPUT')
+.end
diff --git a/current_generator/include/write.frm b/current_generator/include/write.frm
index 5741f5f..a9023d3 100644
--- a/current_generator/include/write.frm
+++ b/current_generator/include/write.frm
@@ -1,129 +1,133 @@
* Write start of C++ header to `OUTPUT'
#procedure WriteHeader(OUTPUT)
#write <`OUTPUT'> "#pragma once"
#write <`OUTPUT'> ""
#write <`OUTPUT'> "#include <complex>"
#write <`OUTPUT'> "#include \"HEJ/helicity.hh\""
#write <`OUTPUT'> "#include \"HEJ/LorentzVector.hh\""
+ #write <`OUTPUT'> "#include \"HEJ/Tensor.hh\""
#write <`OUTPUT'> ""
#write <`OUTPUT'> "namespace HEJ {"
#endprocedure
* Write optimised expression to C++ header `OUTPUT'
#procedure WriteOptimised(OUTPUT,EXPRESSION,NUMHELICITIES,?MOMENTA)
#write <`OUTPUT'> " template<%"
#define FIRST "1"
#do i=1,`NUMHELICITIES'
#ifdef `FIRST'
#undefine FIRST
#write <`OUTPUT'> "Helicity%"
#else
#write <`OUTPUT'> ", Helicity%"
#endif
#enddo
#write <`OUTPUT'> ">"
#write <`OUTPUT'> " std::complex<double> `EXPRESSION'("
#call WriteMomenta(`?MOMENTA')
#write <`OUTPUT'> "\n );\n"
#call WriteOptimisedHelper(`OUTPUT',`EXPRESSION',`NUMHELICITIES',`?MOMENTA')
#endprocedure
*INTERNAL PROCEDURE
#procedure WriteOptimisedHelper(OUTPUT,EXPRESSION,NUMHELICITIES,?REST)
#if `NUMHELICITIES' > 0
#do H={+,-}
#call WriteOptimisedHelper(`OUTPUT',`EXPRESSION',{`NUMHELICITIES'-1},`?REST',`H')
#enddo
#else
#define HELSTRING ""
#define TEMPLATEARGS ""
#define MOMENTA ""
#do ARG={`?REST'}
#if ("`ARG'" == "+") || ("`ARG'" == "-")
* arguments that define helicities
#redefine HELSTRING "`HELSTRING'`ARG'"
#if "`TEMPLATEARGS'" != ""
#redefine TEMPLATEARGS "`TEMPLATEARGS',"
#endif
#if "`ARG'" == "+"
#redefine TEMPLATEARGS "`TEMPLATEARGS'Helicity::plus"
#else
#redefine TEMPLATEARGS "`TEMPLATEARGS'Helicity::minus"
#endif
#else
* arguments that define momenta
#if "`MOMENTA'" != ""
#redefine MOMENTA "`MOMENTA',"
#endif
#redefine MOMENTA "`MOMENTA'`ARG'"
#endif
#enddo
#optimize [`EXPRESSION'`HELSTRING']
+ #write "operations in [`EXPRESSION'`HELSTRING']: `optimvalue_'"
#write <`OUTPUT'> " template<>"
#write <`OUTPUT'> " inline std::complex<double> `EXPRESSION'<%"
* we use a loop here because otherwise FORM will insert line breaks
* if the string is too large
#define FIRST "1"
#do TEMPLATEARG={`TEMPLATEARGS'}
#ifdef `FIRST'
#undefine FIRST
#else
#write <`OUTPUT'> ", %"
#endif
#write <`OUTPUT'> "`TEMPLATEARG'%"
#enddo
#write <`OUTPUT'> ">("
#call WriteMomenta(`MOMENTA')
#write <`OUTPUT'> "\n ) {"
+ #write <`OUTPUT'> " static constexpr std::complex<double> i_{0., 1.};"
+ #write <`OUTPUT'> " (void) i_; //potentially unused"
#if `optimmaxvar_' > 0
#write <`OUTPUT'> " std::complex<double> %"
#define FIRST "1"
#do i=1,`optimmaxvar_'
#ifdef `FIRST'
#undefine FIRST
#write <`OUTPUT'> "Z`i'_%"
#else
#write <`OUTPUT'> ", Z`i'_%"
#endif
#enddo
#write <`OUTPUT'> ";"
#endif
#write <`OUTPUT'> " %O"
#write <`OUTPUT'> " return %E;" [`EXPRESSION'`HELSTRING']
#write <`OUTPUT'> " }\n" [`EXPRESSION'`HELSTRING']
#clearoptimize
#endif
#endprocedure
*INTERNAL PROCEDURE
* Write momenta as C++ function arguments to `OUTPUT'
#procedure WriteMomenta(?MOMENTA)
#define FIRST "1"
#do P={`?MOMENTA'}
#ifdef `FIRST'
#undefine FIRST
#write <`OUTPUT'> " CLHEP::HepLorentzVector const & `P'%"
#else
#write <`OUTPUT'> ",\n CLHEP::HepLorentzVector const & `P'%"
#endif
#enddo
#endprocedure
* Write end of C++ header to `OUTPUT'
#procedure WriteFooter(OUTPUT)
#write <`OUTPUT'> "}"
#endprocedure
diff --git a/doc/developer_manual/currents.tex b/doc/developer_manual/currents.tex
index ed1fc98..d6725ae 100644
--- a/doc/developer_manual/currents.tex
+++ b/doc/developer_manual/currents.tex
@@ -1,657 +1,672 @@
\section{Currents}
\label{sec:currents_impl}
The following section contains a list of all the currents implemented
in \HEJ. Clean up of the code structure is ongoing. All $W$+Jet currents
are located in \texttt{src/Wjets.cc}, all Higgs+Jets currents are
defined in \texttt{src/Hjets.cc}, and pure jet currents are defined in
in \texttt{src/jets.cc}. All of these have their own separate header
files: \texttt{include/HEJ/Wjets.hh}, \texttt{include/HEJ/Hjets.hh} and
\texttt{include/HEJ/jets.hh} respectively.
The naming convention for the current contraction $\left\|S_{f_1 f_2\to f_1
f_2}\right\|^2$ is \lstinline!ME_[Boson]_[subleading-type]_[incoming]!. For
example \lstinline!ME_W_unob_qq! corresponds to the contraction $j_W^\mu
j_{\text{uno}, \mu}$ ($qQ\to \bar{q}WQg$). For bosons on the same side as the
subleading we drop the connecting underscore, e.g. \lstinline!ME_Wuno_qq!
gives $j_{W,\text{uno}}^\mu j_\mu$ ($qQ\to g\bar{q}WQ$).
\subsection{Pure Jets}
\subsubsection{Quark}
\label{sec:current_quark}
\begin{align}
j_\mu(p_i,p_j)=\bar{u}(p_i)\gamma_\mu u(p_j)
\end{align}
The exact for depends on the helicity and direction (forward/backwards) for the quarks. Currently all different contractions of incoming and outgoing states are defined in \lstinline!joi!, \lstinline!jio! and \lstinline!joo!.
\subsubsection{Gluon}
In \HEJ the currents for gluons and quarks are the same, up to a colour factor $K_g/C_F$, where
\begin{align}
K_g(p_1^-, p_a^-) = \frac{1}{2}\left(\frac{p_1^-}{p_a^-} + \frac{p_a^-}{p_1^-}\right)\left(C_A -
\frac{1}{C_A}\right)+\frac{1}{C_A}.
\end{align}
Thus we can just reuse the results from sec.~\ref{sec:current_quark}.
\subsubsection{Single unordered gluon}
Configuration $q(p_a) \to g(p_1) q(p_2) g^*(\tilde{q}_2)$~\cite{Andersen:2017kfc}
\begin{align}
\label{eq:juno}
\begin{split}
&j^{{\rm uno}\; \mu\ cd}(p_2,p_1,p_a) = i \varepsilon_{1\nu} \left( T_{2i}^{c}T_{ia}^d\
\left(U_1^{\mu\nu}-L^{\mu\nu} \right) + T_{2i}^{d}T_{ia}^c\ \left(U_2^{\mu\nu} +
L^{\mu\nu} \right) \right). \\
U_1^{\mu\nu} &= \frac1{s_{21}} \left( j_{21}^\nu j_{1a}^\mu + 2 p_2^\nu
j_{2a}^\mu \right) \qquad \qquad U_2^{\mu\nu} = \frac1{t_{a1}} \left( 2
j_{2a}^\mu p_a^\nu - j_{21}^\mu j_{1a}^\nu \right) \\
L^{\mu\nu} &= \frac1{t_{a2}} \left(-2p_1^\mu j_{2a}^\nu+2p_1.j_{2a}
g^{\mu\nu} + (\tilde{q}_1+\tilde{q}_2)^\nu j_{2a}^\mu + \frac{t_{b2}}{2} j_{2a}^\mu \left(
\frac{p_2^\nu}{p_1.p_2} + \frac{p_b^\nu}{p_1.p_b} \right) \right) ,
\end{split}
\end{align}
$j^{{\rm uno}\; \mu}$ is currently not calculated as a separate current, but always as needed for the ME (i.e. in \lstinline!ME_unob_XX!).
\subsubsection{Extremal $q\bar{q}$}
In Pure jets we also include the subleading process which arises when
an incoming gluon splits into a $q\bar{q}$ pair. This splitting impact
factor is related to the unordered current by simple means of a
crossing symmetry.
\subsubsection{Central $q\bar{q}$}
The final subleading process type in the Pure Jets case is Central
$q\bar{q}$. In this process type, we have two currents scattering off
of each other, but this time, via an effective vertex, which connects
together two FKL chains. Each FKL chain t-channel gluon splits into a
$q\bar{q}$ and this results in a quark and anti-quark in between the
most forward and backward jets. One can see an example of such a
process in Figure \ref{fig:qqbarcen_example}.
\\
As the new central $q\bar{q}$ piece contains the quark propagator, we will treat
this as part of the skeleton process. This means that we do not impose strong ordering
between the $q\bar{q}$-pair taking
\begin{align}
\label{eq:cenqqbarraporder}
y_1 \ll y_2,y_3 \ll y_4.
\end{align}
The HEJ Matrix element for this process can be calculated as:
\begin{align}
\label{eq:Mcentral}
i\mathcal{M} &= g_s^4 T^d_{1a} T^e_{4b}\ \frac{j_{\mu}(p_a,p_1)\ X^{ab\, \mu
\nu}_{{\rm cen}}(p_2,p_3,q_1,q_3)\
j_{\nu}(p_b,p_4)}{t_{a1}t_{b4}}.
\end{align}
where $X^{\mu \nu}_{\rm cen}$ is given by:
\begin{align}
X^{\mu \nu}_{\rm cen} &=\frac{f^{ced}T^c_{23}}{s_{23}}
\left(\eta^{\mu \nu} X_{sym}^\sigma + V^{\mu \nu \sigma}_{3g} \right)
\bar{u}(p_2) \gamma^\sigma u(p_3) \\ & \qquad + \frac{i
T^d_{2q}T^e_{q3}}{(q_1-p_2)^2} X^{\mu\nu}_{qprop} - \frac{i
T^e_{2q}T^d_{q3}}{(q_1-p_3)^2} X^{\mu\nu}_{crossed},
\end{align}
\begin{figure}[t]
\centering
\includegraphics[]{figures/Cenqqbar_jx}
\caption{An example of a Central $q\bar{q}$ process}
\label{fig:qqbarcen_example}
\end{figure}
\subsection{Higgs}
Different rapidity orderings \todo{give name of functions}
\begin{enumerate}
\item $qQ\to HqQ/qHQ/qQH$ (any rapidity order, full LO ME) $\Rightarrow$ see~\ref{sec:V_H}
\item $qg\to Hqg$ (Higgs outside quark) $\Rightarrow$ see~\ref{sec:V_H}
\item $qg\to qHg$ (central Higgs) $\Rightarrow$ see~\ref{sec:V_H}
\item $qg\to qgH$ (Higgs outside gluon) $\Rightarrow$ see~\ref{sec:jH_mt}
\item $gg\to gHg$ (central Higgs) $\Rightarrow$ see~\ref{sec:V_H}
\item $gg\to ggH$ (Higgs outside gluon) $\Rightarrow$ see~\ref{sec:jH_mt}
\end{enumerate}
\subsubsection{Higgs gluon vertex}
\label{sec:V_H}
The coupling of the Higgs boson to gluons via a virtual quark loop can be written as
\begin{align}
\label{eq:VH}
V^{\mu\nu}_H(q_1, q_2) = \mathgraphics{V_H.pdf} &=
\frac{\alpha_s m^2}{\pi v}\big[
g^{\mu\nu} T_1(q_1, q_2) - q_2^\mu q_1^\nu T_2(q_1, q_2)
\big]\, \\
&\xrightarrow{m \to \infty} \frac{\alpha_s}{3\pi
v} \left(g^{\mu\nu} q_1\cdot q_2 - q_2^\mu q_1^\nu\right).
\end{align}
The outgoing momentum of the Higgs boson is $p_H = q_1 - q_2$.
As a contraction with two currents this by implemented in \lstinline!cHdot! inside \texttt{src/Hjets.cc}.
The form factors $T_1$ and $T_2$ are then given by~\cite{DelDuca:2003ba}
\begin{align}
\label{eq:T_1}
T_1(q_1, q_2) ={}& -C_0(q_1, q_2)\*\left[2\*m^2+\frac{1}{2}\*\left(q_1^2+q_2^2-p_H^2\right)+\frac{2\*q_1^2\*q_2^2\*p_H^2}{\lambda}\right]\notag\\
&-\left[B_0(q_2)-B_0(p_H)\right]\*\frac{q_2^2}{\lambda}\*\left(q_2^2-q_1^2-p_H^2\right)\notag\\
&-\left[B_0(q_1)-B_0(p_H)\right]\*\frac{q_1^2}{\lambda}\*\left(q_1^2-q_2^2-p_H^2\right)-1\,,\displaybreak[0]\\
\label{eq:T_2}
T_2(q_1, q_2) ={}& C_0(q_1,
q_2)\*\left[\frac{4\*m^2}{\lambda}\*\left(p_H^2-q_1^2-q_2^2\right)-1-\frac{4\*q_1^2\*q_2^2}{\lambda}
-
\frac{12\*q_1^2\*q_2^2\*p_H^2}{\lambda^2}\*\left(q_1^2+q_2^2-p_H^2\right)\right]\notag\\
&-\left[B_0(q_2)-B_0(p_H)\right]\*\left[\frac{2\*q_2^2}{\lambda}+\frac{12\*q_1^2\*q_2^2}{\lambda^2}\*\left(q_2^2-q_1^2+p_H^2\right)\right]\notag\\
&-\left[B_0(q_1)-B_0(p_H)\right]\*\left[\frac{2\*q_1^2}{\lambda}+\frac{12\*q_1^2\*q_2^2}{\lambda^2}\*\left(q_1^2-q_2^2+p_H^2\right)\right]\notag\\
&-\frac{2}{\lambda}\*\left(q_1^2+q_2^2-p_H^2\right)\,,
\end{align}
where we have used the scalar bubble and triangle integrals
\begin{align}
\label{eq:B0}
B_0\left(p\right) ={}& \int \frac{d^dl}{i\pi^{\frac{d}{2}}}
\frac{1}{\left(l^2-m^2\right)\left((l+p)^2-m^2\right)}\,,\\
\label{eq:C0}
C_0\left(p,q\right) ={}& \int \frac{d^dl}{i\pi^{\frac{d}{2}}} \frac{1}{\left(l^2-m^2\right)\left((l+p)^2-m^2\right)\left((l+p-q)^2-m^2\right)}\,,
\end{align}
and the K\"{a}ll\'{e}n function
\begin{equation}
\label{eq:lambda}
\lambda = q_1^4 + q_2^4 + p_H^4 - 2\*q_1^2\*q_2^2 - 2\*q_1^2\*p_H^2- 2\*q_2^2\*p_H^2\,.
\end{equation}
The Integrals as such are provided by \QCDloop{} (see wrapper functions \lstinline!B0DD! and \lstinline!C0DD! in \texttt{src/Hjets.cc}).
In the code we are sticking to the convention of~\cite{DelDuca:2003ba}, thus instead of the $T_{1/2}$ we implement (in the functions \lstinline!A1! and \lstinline!A2!)
\begin{align}
\label{eq:A_1}
A_1(q_1, q_2) ={}& \frac{i}{16\pi^2}\*T_2(-q_1, q_2)\,,\\
\label{eq:A_2}
A_2(q_1, q_2) ={}& -\frac{i}{16\pi^2}\*T_1(-q_1, q_2)\,.
\end{align}
\subsubsection{Peripheral Higgs emission - Finite quark mass}
\label{sec:jH_mt}
We describe the emission of a peripheral Higgs boson close to a
scattering gluon with an effective current. In the following we consider
a lightcone decomposition of the gluon momenta, i.e. $p^\pm = E \pm p_z$
and $p_\perp = p_x + i p_y$. The incoming gluon momentum $p_a$ defines
the $-$ direction, so that $p_a^+ = p_{a\perp} = 0$. The outgoing
momenta are $p_1$ for the gluon and $p_H$ for the Higgs boson. We choose
the following polarisation vectors:
\begin{equation}
\label{eq:pol_vectors}
\epsilon_\mu^\pm(p_a) = \frac{j_\mu^\pm(p_1, p_a)}{\sqrt{2}
\bar{u}^\pm(p_a)u^\mp(p_1)}\,, \quad \epsilon_\mu^{\pm,*}(p_1) = -\frac{j_\mu^\pm(p_1, p_a)}{\sqrt{2}
\bar{u}^\mp(p_1)u^\pm(p_a)}\,.
\end{equation}
Following~\cite{DelDuca:2001fn}, we introduce effective polarisation
vectors to describe the contraction with the Higgs-boson production
vertex eq.~\eqref{eq:VH}:
\begin{align}
\label{eq:eps_H}
\epsilon_{H,\mu}(p_a) = \frac{T_2(p_a, p_a-p_H)}{(p_a-p_H)^2}\big[p_a\cdot
p_H\epsilon_\mu(p_a) - p_H\cdot\epsilon(p_a) p_{a,\mu}\big]\,,\\
\epsilon_{H,\mu}^*(p_1) = -\frac{T_2(p_1+p_H, p_1)}{(p_1+p_H)^2}\big[p_1\cdot
p_H\epsilon_\mu^*(p_1) - p_H\cdot\epsilon^*(p_1) p_{1,\mu}\big]\,,
\end{align}
We also employ the usual short-hand notation
\begin{equation}
\label{eq:spinor_helicity}
\spa i.j = \bar{u}^-(p_i)u^+(p_j)\,,\qquad \spb i.j =
\bar{u}^+(p_i)u^-(p_j)\,, \qquad[ i | H | j\rangle = j_\mu^+(p_i, p_j)p_H^\mu\,.
\end{equation}
Without loss of generality, we consider only the case where the incoming
gluon has positive helicity. The remaining helicity configurations can
be obtained through parity transformation.
Labelling the effective current by the helicities of the gluons we obtain
for the same-helicity case
\begin{equation}
\label{eq:jH_same_helicity}
\begin{split}
j_{H,\mu}^{++}{}&(p_1,p_a,p_H) =
\frac{m^2}{\pi v}\bigg[\\
&-\sqrt{\frac{2p_1^-}{p_a^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spb a.1}\epsilon^{+,*}_{H,\mu}(p_1)
+\sqrt{\frac{2p_a^-}{p_1^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spa 1.a}\epsilon^{+}_{H,\mu}(p_a)\\
&+ [1|H|a\rangle \bigg(
\frac{\sqrt{2}}{\spa 1.a}\epsilon^{+}_{H,\mu}(p_a)
+ \frac{\sqrt{2}}{\spb a.1}\epsilon^{+,*}_{H,\mu}(p_1)-\frac{\spa 1.a T_2(p_a, p_a-p_H)}{\sqrt{2}(p_a-p_H)^2}\epsilon^{+,*}_{\mu}(p_1)\\
&
\qquad
-\frac{\spb a.1 T_2(p_1+p_H,
p_1)}{\sqrt{2}(p_1+p_H)^2}\epsilon^{+}_{\mu}(p_a)-\frac{RH_4}{\sqrt{2}\spb a.1}\epsilon^{+,*}_{\mu}(p_1)+\frac{RH_5}{\sqrt{2}\spa 1.a}\epsilon^{+}_{\mu}(p_a)
\bigg)\\
&
- \frac{[1|H|a\rangle^2}{2 t_1}(p_{a,\mu} RH_{10} - p_{1,\mu} RH_{12})\bigg]
\end{split}
\end{equation}
with $t_1 = (p_a-p_1)^2$, $t_2 = (p_a-p_1-p_H)^2$ and $R = 8 \pi^2$. Eq.~\eqref{eq:jH_same_helicity}
is implemented by \lstinline!g_gH_HC! in \texttt{src/Hjets.cc}
\footnote{\lstinline!g_gH_HC! and \lstinline!g_gH_HNC! includes an additional
$1/t_2$ factor, which should be in the Matrix element instead.}.
The currents with a helicity flip is given through
\begin{equation}
\label{eq:jH_helicity_flip}
\begin{split}
j_{H,\mu}^{+-}{}&(p_1,p_a,p_H) =
\frac{m^2}{\pi v}\bigg[\\
&-\sqrt{\frac{2p_1^-}{p_a^-}}\frac{p_{1\perp}^*}{|p_{1\perp}|}\frac{t_2}{\spb
a.1}\epsilon^{-,*}_{H,\mu}(p_1)
+\sqrt{\frac{2p_a^-}{p_1^-}}\frac{p_{1\perp}}{|p_{1\perp}|}\frac{t_2}{\spb a.1}\epsilon^{+}_{H,\mu}(p_a)\\
&+ [1|H|a\rangle \left(
\frac{\sqrt{2}}{\spb a.1} \epsilon^{-,*}_{H,\mu}(p_1)
-\frac{\spa 1.a T_2(p_a, p_a-p_H)}{\sqrt{2}(p_a-p_H)^2}\epsilon^{-,*}_{\mu}(p_1)
- \frac{RH_4}{\sqrt{2}\spb a.1}\epsilon^{-,*}_{\mu}(p_1)\right)
\\
&+ [a|H|1\rangle \left(
\frac{\sqrt{2}}{\spb a.1}\epsilon^{+}_{H,\mu}(p_a)
-\frac{\spa 1.a
T_2(p_1+p_H,p_1)}{\sqrt{2}(p_1+p_H)^2}\epsilon^{+}_{\mu}(p_a)
+\frac{RH_5}{\sqrt{2}\spb a.1}\epsilon^{+}_{\mu}(p_a)
\right)\\
& - \frac{[1|H|a\rangle [a|H|1\rangle}{2 \spb a.1 ^2}(p_{a,\mu} RH_{10} - p_{1,\mu}
RH_{12})\\
&+ \frac{\spa 1.a}{\spb a.1}\bigg(RH_1p_{1,\mu}-RH_2p_{a,\mu}+2
p_1\cdot p_H \frac{T_2(p_1+p_H, p_1)}{(p_1+p_H)^2} p_{a,\mu}
\\
&
\qquad- 2p_a \cdot p_H \frac{T_2(p_a, p_a-p_H)}{(p_a-p_H)^2} p_{1,\mu}+ T_1(p_a-p_1, p_a-p_1-p_H)\frac{(p_1 + p_a)_\mu}{t_1}\\
&\qquad-\frac{(p_1+p_a)\cdot p_H}{t_1} T_2(p_a-p_1, p_a-p_1-p_H)(p_1 - p_a)_\mu
\bigg)
\bigg]\,,
\end{split}
\end{equation}
and implemented by \lstinline!g_gH_HNC! again in \texttt{src/Hjets.cc}.
If we instead choose the gluon momentum in the $+$ direction, so that
$p_a^- = p_{a\perp} = 0$, the corresponding currents are obtained by
replacing $p_1^- \to p_1^+, p_a^- \to p_a^+,
\frac{p_{1\perp}}{|p_{1\perp}|} \to -1$ in the second line of
eq.~\eqref{eq:jH_same_helicity} and eq.~\eqref{eq:jH_helicity_flip} (see variables \lstinline!ang1a! and \lstinline!sqa1! in the implementation).
The form factors $H_1,H_2,H_4,H_5, H_{10}, H_{12}$ are given in~\cite{DelDuca:2003ba}, and are implemented under their name in \texttt{src/Hjets.cc}. They reduce down to fundamental QCD integrals, which are again provided by \QCDloop.
\subsubsection{Peripheral Higgs emission - Infinite top mass}
\label{sec:jH_eff}
To get the result with infinite top mass we could either take the limit $m_t\to \infty$ in~\eqref{eq:jH_helicity_flip} and~\eqref{eq:jH_same_helicity}, or use the \textit{impact factors} as given in~\cite{DelDuca:2003ba}. Both methods are equivalent, and lead to the same result. For the first one would find
\begin{align}
\lim_{m_t\to\infty} m_t^2 H_1 &= i \frac{1}{24 \pi^2}\\
\lim_{m_t\to\infty} m_t^2 H_2 &=-i \frac{1}{24 \pi^2}\\
\lim_{m_t\to\infty} m_t^2 H_4 &= i \frac{1}{24 \pi^2}\\
\lim_{m_t\to\infty} m_t^2 H_5 &=-i \frac{1}{24 \pi^2}\\
\lim_{m_t\to\infty} m_t^2 H_{10} &= 0 \\
\lim_{m_t\to\infty} m_t^2 H_{12} &= 0.
\end{align}
\todo{double check this, see James thesis eq. 4.33}
However only the second method is implemented in the code through \lstinline!C2gHgp!
and \lstinline!C2gHgm! inside \texttt{src/Hjets.cc}, each function
calculates the square of eq. (4.23) and (4.22) from~\cite{DelDuca:2003ba} respectively.
\subsection{$W$+Jets}
\label{sec:currents_W}
\subsubsection{Quark+$W$}
\begin{figure}
\centering
\begin{minipage}[b]{0.2\textwidth}
\includegraphics[width=\textwidth]{Wbits.pdf}
\end{minipage}
\begin{minipage}[b]{0.1\textwidth}
\centering{=}
\vspace{0.7cm}
\end{minipage}
\begin{minipage}[b]{0.2\textwidth}
\includegraphics[width=\textwidth]{Wbits2.pdf}
\end{minipage}
\begin{minipage}[b]{0.1\textwidth}
\centering{+}
\vspace{0.7cm}
\end{minipage}
\begin{minipage}[b]{0.2\textwidth}
\includegraphics[width=\textwidth]{Wbits3.pdf}
\end{minipage}
\caption{The $j_W$ current is constructed from the two diagrams which
contribute to the emission of a $W$-boson from a given quark line.}
\label{fig:jW}
\end{figure}
For a $W$ emission we require a fermion. The $j_W$ current is actually a sum of
two separate contributions, see figure~\ref{fig:jW}, one with a $W$-emission
from the initial state, and one with the $W$-emission from the final state.
Mathematically this can be seen as the following two
terms~\cite{Andersen:2012gk}\todo{cite W subleading paper}:
\begin{align}
\label{eq:Weffcur1}
j_W^\mu(p_a,p_{\ell},p_{\bar{\ell}}, p_1) =&\ \frac{g_W^2}{2}\
\frac1{p_W^2-M_W^2+i\ \Gamma_W M_W}\ \bar{u}^-(p_\ell) \gamma_\alpha
v^-(p_{\bar\ell})\nonumber \\
& \cdot \left( \frac{ \bar{u}^-(p_1) \gamma^\alpha (\slashed{p}_W +
\slashed{p}_1)\gamma^\mu u^-(p_a)}{(p_W+p_1)^2} +
\frac{ \bar{u}^-(p_1)\gamma^\mu (\slashed{p}_a + \slashed{p}_W)\gamma^\alpha u^-(p_a)}{(p_a-p_W)^2} \right).
\end{align}
There are a couple of subtleties here. There is a minus sign
distinction between the quark-anti-quark cases due to the fermion flow
of the propagator in the current. Note that the type of $W$ emission
(+ or -) will depend on the quark flavour, and that the handedness of
the quark-line is given by whether its a quark or anti-quark.
The FKL $W$ current is
\todo{Text. The following is what's currently implemented in \lstinline!jW!}
\begin{align}
\label{eq:jW-}
j^-_{W}(p_a, p_1, p_{\bar{l}}, p_{l}) ={}&
\frac{2 \spa 1.l}{(p_1+p_l+p_{\bar{l}})^2}\*\Big(
\spa 1.{\bar{l}} j^-(p_1, p_a) + \spb l.{\bar{l}}
j^-(p_{l}, p_a)\Big)\notag\\
& + \frac{2\spb a.{\bar{l}}}{(p_a-p_W)^2}\Big(
\spa 1.l j^-(p_l, p_a) + \spa a.l
j^-(p_1, p_a)
\Big)\,,\\
\label{eq:jW+}
j^+_{W}(p_a, p_1, p_{\bar{l}}, p_{l}) =& \big[j^-_{W}(p_a, p_1, p_l, p_{\bar{l}})\big]^*\,,
\end{align}
where the negative-helicity current is used for emission off a quark
line and the positive-helicity current for emissions off antiquark.
\subsubsection{$W$+uno}
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{wuno1}
\caption{}
\label{fig:U1diags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{wuno2}
\caption{}
\label{fig:U2diags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{wuno3}
\caption{}
\label{fig:Cdiags}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{wuno4}
\caption{}
\label{fig:Ddiags}
\end{subfigure}
\vspace{0.4cm}
\caption{Examples of each of the four categories of Feynman diagram which
contribute to at leading-order; there are twelve in total. (a) is an example where the gluon and $W$
boson are emitted from the same quark line and the gluon comes after the
$t$-channel propagator. In (b), the gluon and $W$ boson are emitted from
the same quark line and the gluon comes before the $t$-channel proagator.
In (c) the gluon is emitted from the $t$-channel gluon and in (d) the gluon
is emitted from the $b$--$3$ quark line.}
\label{fig:Wunodiags}
\end{figure}
It is necessary to include subleading processes in $W$+Jets also. All of
these currents have been built for the \lstinline!Tensor! Class detailed in
section~\ref{sec:tensor}. Similarly to the pure jet case, the uno currents are
not calculated separately, and only in the ME functions when required
in the \texttt{src/Wjets.cc} file. For unordered emissions a new
current is required, $j_{W,{\rm uno}}$, it is only non-zero for
$h_a=h_1=-$ and hence we have suppressed its helicity indices. It is
derived from the 12 leading-order Feynman diagrams in the QMRK
limit (see figure~\ref{fig:Wunodiags}). Using $T^m_{ij}$ represent fundamental
colour matrices between quark state $i$ and $j$ with adjoint index $m$ we find
\begin{align}\label{eq:wunocurrent}
\begin{split}
j^{d\,\mu}_{\rm W,uno}(p_a,p_1,p_2,p_\ell,p_{\bar{\ell}}) =& \ i \varepsilon_{\nu}(p_1)\
\bar{u}^-(p_\ell) \gamma_\rho v^-(p_{\bar \ell}) \\ & \quad \times\
\left(T^1_{2i} T^d_{ia} (\tilde U_1^{\nu\mu\rho}-\tilde L^{\nu\mu\rho}) +
T^d_{2i} T^1_{ia} (\tilde U_2^{\nu\mu\rho}+\tilde L^{\nu\mu\rho}) \right),
\end{split}
\end{align}
where expressions for $\tilde U_{1,2}^{\nu\mu\rho}$ and $\tilde L^{\nu\mu\rho}$
are given as:
\begin{align}
\label{eq:U1tensor}
\begin{split}
\tilde U_1^{\nu\mu\rho} ={}&\frac{\langle 2|\nu (\slashed{p}_2+ \slashed{p}_1)\mu (\slashed{p}_a - \slashed{p}_W)\rho P_L |a\rangle }{s_{12}t_{aW}} + \frac{\langle 2|\nu (\slashed{p}_2+ \slashed{p}_1)\rho P_L (\slashed{p}_2+\slashed{p}_1 + \slashed{p}_W)\mu |a\rangle }{s_{12}s_{12W}} \\
&+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \nu
(\slashed{p}_1 + \slashed{p}_2+\slashed{p}_W)\mu |a\rangle}{s_{2W}s_{12W}}\,,
\end{split}\\
\label{eq:U2tensor}
\begin{split}
\tilde U_2^{\nu\mu\rho} ={}&\frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W-\slashed{p}_1)\nu (\slashed{p}_a - \slashed{p}_W)\rho P_L |a\rangle }{t_{aW1}t_{aW}} + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W- \slashed{p}_1)\rho P_L (\slashed{p}_a-\slashed{p}_1) \nu |a\rangle }{t_{a1W}t_{a1}} \\
&+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \mu
(\slashed{p}_a-\slashed{p}_1)\nu |a\rangle}{s_{2W}t_{a1}}\,,
\end{split}\\
\label{eq:Ltensor}
\begin{split}
\tilde L^{\nu\mu\rho} ={}& \frac{1}{t_{aW2}}\left[
\frac{\langle 2 | \sigma (\slashed{p}_a-\slashed{p}_W)\rho|a\rangle}{t_{aW}}
+\frac{\langle 2 | \rho (\slashed{p}_2+\slashed{p}_W)\sigma|a\rangle}{s_{2W}}
\right]\\
&\times \left\{\left(\frac{p_b^\nu}{s_{1b}} + \frac{p_3^\nu}{s_{13}}\right)(q_1-p_1)^2g^{\mu\sigma}+(2q_1-p_1)^\nu g^{\mu\sigma} - 2p_1^\mu g^{\nu\sigma} + (2p_1-q_1)^\sigma g^{\mu\nu} \right\}\,,
\end{split}
\end{align}
where $s_{ij\dots} = (p_i + p_j + \dots)^2, t_{ij\dots} = (p_i - p_j - \dots)^2$ and $q_1 = p_a-p_2-p_W$.
\subsubsection{$W$+Extremal $\mathbf{q\bar{q}}$}
The $W$+Jet sub-leading processes containing an extremal $q\bar{q}$ are
related by crossing symmetry to the $W$+Jet unordered processes. This
means that one can simply perform a crossing symmetry argument on
eq.~\ref{eq:wunocurrent} to arrive at the extremal $q\bar{q}$ current
required.We show the basic structure of the extremal $q\bar{q}$
current in figure~\ref{fig:qgimp}, neglecting the $W$-emission for
simplicity.
\begin{figure}
\centering
\includegraphics[width=0.3\textwidth]{{qqbarex_schem}}
\caption{Schematic structure of the $gq \to \bar{Q}Qq$ amplitude in the limit
$y_1 \sim y_2 \ll y_3$}
\label{fig:qgimp}
\end{figure}
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex1}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex2}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex4}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex5}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarex3}
\end{subfigure}
\caption{The five tree-level graphs which contribute to the process $gq \to \bar{Q}Qq$.}
\label{fig:qg_qQQ_graphs}
\end{figure}
We can obtain the current for $g\rightarrow W q \bar{q}$ by evaluating
the current for $W$ plus unordered emissions with the normal arguments
$p_a \leftrightarrow -p_1 $ interchanged. This is a non-trivial
statement: due to the minimality of the approximations made, the
crossing symmetry normally present in the full amplitude may be
extended to the factorised current.
We must again note that swapping $p_a \leftrightarrow -p_1$ will lead
to $u$-spinors with momenta with negative energy. These are identical
to $v$-spinors with momenta with positive energy, up to an overall
phase which is common to all terms, and can therefore be neglected.
Mathematically, this is given by:
\begin{align}\label{eq:crossedJ}
j^\mu_{\rm W,g\to Q\bar{Q}}(p_a,p_1,p_2,p_\ell,p_{\bar{\ell}}) =i \varepsilon_{g\nu}
\langle \ell | \rho | \bar \ell \rangle_L
\left(T^1_{2i} T^d_{ia} (\tilde U_{1,X}^{\nu\mu\rho}-\tilde L^{\nu\mu\rho}_X) + T^d_{2i} T^1_{ia} (\tilde U_{2,X}^{\nu\mu\rho}+\tilde L_X^{\nu\mu\rho}) \right),
\end{align}
where the components are now given by
\begin{align}
\label{eq:U1tensorX}
\begin{split}
\tilde U_{1,X}^{\nu\mu\rho} =&\frac{\langle 2|\nu (\slashed{p}_a- \slashed{p}_2)\mu (\slashed{p}_1 + \slashed{p}_W)\rho P_L |1\rangle }{t_{a2}s_{1W}} + \frac{\langle 2|\nu (\slashed{p}_a- \slashed{p}_2)\rho P_L (\slashed{p}_a-\slashed{p}_2 - \slashed{p}_W)\mu |1\rangle }{t_{a2}t_{a2W}} \\
&- \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \nu
(\slashed{p}_a - \slashed{p}_2-\slashed{p}_W)\mu
|1\rangle}{s_{2W}t_{a2W}}, \\
\tilde U_{2,X}^{\nu\mu\rho} =&-\frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W-\slashed{p}_1)\nu (\slashed{p}_1 + \slashed{p}_W)\rho P_L |1\rangle }{t_{aW1}s_{1W}} + \frac{\langle 2|\mu (\slashed{p}_a-\slashed{p}_W- \slashed{p}_1)\rho P_L (\slashed{p}_a-\slashed{p}_1) \nu |1\rangle }{t_{a1W}t_{a1}} \\
&+ \frac{\langle 2|\rho P_L (\slashed{p}_2+ \slashed{p}_W) \mu
(\slashed{p}_a-\slashed{p}_1)\nu |1\rangle}{s_{2W}t_{a1}}, \\
\tilde L^{\nu\mu\rho}_X &= \frac{1}{s_{W12}}\left[-\frac{\langle 2 |\sigma (\slashed{p}_1 + \slashed{p}_W) \rho P_L | 1\rangle}{s_{1W}} + \frac{\langle 2 |\rho P_L (\slashed{p}_2 + \slashed{p}_W) \sigma | 1\rangle }{s_{2W}} \right] \\
&\vphantom{+\frac{1}{t_{aW2}}}\quad\cdot \left( -\left(
\frac{p_b^\nu}{s_{ab}} + \frac{p_n^\nu}{s_{an}} \right) (q_1+p_a)^2 g^{\sigma\mu}+ g^{\sigma \mu} (2q_1 +p_a)^\nu - g^{\mu \nu}(2p_a+q_1)^\sigma+ 2g^{\nu \sigma}p_a^\mu \right),
\end{split}
\end{align}
where $q_1=-(p_1+p_2+p_W)$. Notice in particular the similarity to the $W$+uno scenario (from which
this has been derived).
\subsubsection{Central $\mathbf{q\bar{q}}$ Vertex}
The final subleading process in the $W$+Jet case is the Central
$q\bar{q}$ vertex. This subleading process does not require an altered
current, but an effective vertex which is contracted with two regular
\HEJ currents. This complexity is dealt with nicely by the \lstinline!Tensor!
class, which is detailed in section~\ref{sec:tensor}.
The $W$-emission can be from the central effective vertex (scenario
-dealt with by the function \texttt{jM2WqqtoqQQq()} in the file
+dealt with by the function \lstinline!ME_WCenqqx_qq! in the file
\texttt{src/Wjets.cc}); or from either of the external quark legs
-(scenario dealt with by \texttt{jM2WqqtoqQQqW()} in same file). In
+(scenario dealt with by \lstinline!ME_W_Cenqqx_qq! in same file). In
the pure jets case, there are 7 separate diagrams which contribute to
this, which can be seen in figure~\ref{fig:qq_qQQq_graphs}. In the $W$+Jets
case, there are then 45 separate contributions.
\begin{figure}
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen1}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen2}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen3}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen4}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen5}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen6}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics{qqbarcen7}
\end{subfigure}
\caption{All Feynman diagrams which contribute to $qq' \to qQ\bar{Q}q'$ at
leading order.}
\label{fig:qq_qQQq_graphs}
\end{figure}
The end result is of the effective vertex, after derivation, is:
\begin{align}
\label{eq:EffectiveVertexqqbar}
\begin{split}
V^{\mu\nu}_{\text{Eff}}=&
\frac{C_1}{s_{23AB}}\left(X^{\mu\nu\sigma}_{1a}\hat{t_1} + X^{\mu\nu\sigma}_{4b}\hat{t_3} +V^{\mu\nu\sigma}_{3g}\right)J_{\text{V} \sigma}(p_2,p_A,p_B,p_3)
\\
&\quad +iC_2X^{\mu\nu}_{Unc}+iC_3X^{\mu\nu}_{Cro},
\end{split}
\end{align}
where:
\begin{eqnarray}
\begin{split}
C_1=&T^e_{1q}T^g_{qa}T^e_{23}T^g_{4b} -
T^g_{1q}T^e_{qa}T^e_{23}T^g_{4b} = f^{egc}T^c_{1a}T^e_{23}T^g_{4b},
\\
C_2=&T^g_{1a}T^g_{2q}T^{g'}_{q3}T^{g'}_{4b}
\\
C_3=&T^g_{1a}T^{g'}_{2q}T^g_{q3}T^{g'}_{4b}
\end{split}
\end{eqnarray}
-are the colour factors of different contributions. The following
-tensor structures correspond to groupings of diagrams in
+are the colour factors of different contributions and
+\todo{Different convention for $J_V$}
+\begin{equation}
+ \label{eq:J_V}
+ J_V^\mu(p_2,p_l,p_{\bar{l}},p_3)=-i\left( \frac{ \bar{u}_2 \gamma^\nu (\slashed{p}_2 +
+ \slashed{p}_l +
+ \slashed{p}_{\bar{l}}) \gamma^\mu u_3}{s_{2l\bar{l}}} - \frac{\bar u_2
+ \gamma^\mu(\slashed{p}_3 + \slashed{p}_l + \slashed{p}_{\bar{l}}) \gamma^\nu
+ u_3}{s_{3l\bar{l}}} \right) [\bar{u}_l \gamma_\nu u_{\bar{l}}].
+\end{equation}
+The following tensor structures correspond to groupings of diagrams in
figure~\ref{fig:qq_qQQq_graphs}.
\begin{eqnarray}
- \label{eq:1aFixed}
+ \label{eq:X_1a}
X_{1a}^{\mu\nu\sigma} &=
\frac{-g^{\mu\nu}}{s_{23AB}\hat{t_3}}\left(\frac{p^\sigma_a}{(s_{a23AB})} +
\frac{p^\sigma_1}{(s_{123AB})}\right)
\\
- \label{eq:4bFixed}
+ \label{eq:X_4b}
X_{4b}^{\mu\nu\sigma}
&=\frac{g^{\mu\nu}}{s_{23AB}\hat{t_1}}\left(\frac{p^\sigma_b}{(s_{23bAB})}+\frac{p^\sigma_4}{(s_{234AB}}\right)
\end{eqnarray}
correspond to the first and second row of diagrams in figure~\ref{fig:qq_qQQq_graphs}.
\begin{align}
\label{eq:3GluonWEmit}
\begin{split}
- X^{\mu\nu}_{3g}=\frac{1}{
+ 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
-\texttt{MSymW()} in \texttt{src/Wjets.cc} returns a \lstinline!Tensor!
-containing the information from these 5 groupings of contributions (30 diagrams
-in total).
+figure~\ref{fig:qq_qQQq_graphs}. One notes that all of these
+contributions have the same colour factor, and as such we can group
+them together nicely before summing over helicities etc. As such, the
+function \lstinline!M_sym_W! that is part of the header
+\texttt{W\_central\_qqbar.hh} generated by the current generator
+(section~\ref{sec:cur_gen}) returns a contraction of the above tensor
+containing the information from these 5 groupings of contributions (30
+diagrams in total).
\begin{align}
+ \label{eq:X_Unc}
\begin{split}
X^{\mu\nu}_{Unc}=\frac{\langle A|\sigma P_L|B\rangle}{\hat{t_1}\hat{t_3}} \bar{u}_2&\left[-\frac{
\gamma^\sigma P_L(\slashed{p}_2+\slashed{p}_A+\slashed{p}_B)\gamma^\mu
(\slashed{q}_3+ \slashed{p}_3)\gamma^\nu}{(s_{2AB})(t_{unc_{2}})}\right.+
\\
&\qquad\left. \frac{\gamma^\mu(\slashed{q}_1-\slashed{p}_2)\gamma^\sigma
P_L(\slashed{q}_3+\slashed{p}_3)\gamma^\nu}{(t_{unc_{1}})(t_{unc_{2}})}\right. +
\\
&\qquad\qquad\left. \frac{\gamma^\mu(\slashed{q}_1-\slashed{p}_2)\gamma^\nu(\slashed{p}_3+\slashed{p}_A+\slashed{p}_B)\gamma^\sigma P_L
}{(t_{unc_1})(s_{3AB})}\right]v_3
\end{split}
\end{align}
corresponds to the diagram on the right of row three in
-figure~\ref{fig:qq_qQQq_graphs}. This contribution to the effective vertex can
-be obtained in the code with the function \texttt{MUncW()} in file
-\texttt{src/Wjets.cc}.
+figure~\ref{fig:qq_qQQq_graphs}. This contribution to the current
+contraction can be obtained in the code with the function
+\lstinline!M_uncross_W! in the autogenerated header \texttt{W\_central\_qqbar.hh}.
\begin{align}
\begin{split}
+ \label{eq:X_Cro}
X^{\mu\nu}_{Cro}=\frac{\langle
A|\sigma P_L|B\rangle}{\hat{t_1}\hat{t_3}} \bar{u}_2&\left[-\frac{
\gamma^\nu(\slashed{q}_3+\slashed{p}_2)\gamma^\mu
(\slashed{p}_3+\slashed{p}_A+\slashed{p}_B)\gamma^\sigma P_L}{(t_{cro_1})(s_{3AB})}\right.+
\\
&\qquad\left. \frac{\gamma^\nu(\slashed{q}_3+\slashed{p}_2)\gamma^\sigma
P_L(\slashed{q}_1-\slashed{p}_3)\gamma^\mu}{(t_{cro_{1}})(t_{cro_{2}})}\right.+
\\ &\qquad\qquad\left
. \frac{\gamma^\sigma P_L(\slashed{p}_2+\slashed{p}_A+\slashed{p}_B)\gamma^\nu(\slashed{q}_1-\slashed{p}_3)\gamma^\mu
}{(s_{2AB})(t_{cro_2})}\right]v_3
\end{split}
\end{align}
-corresponds to the last diagram in figure~\ref{fig:qq_qQQq_graphs}. This
-contribution to the effective vertex can be obtained in the code with
-the function \texttt{MCroW()} in file \texttt{src/Wjets.cc}.
+corresponds to the last diagram in
+figure~\ref{fig:qq_qQQq_graphs}. This contribution to the current
+contraction can be obtained in the code with the function
+\lstinline!M_cross_W! in the generated header
+\texttt{W\_central\_qqbar.hh}.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "developer_manual"
%%% End:
diff --git a/doc/developer_manual/developer_manual.tex b/doc/developer_manual/developer_manual.tex
index 03deb07..80293c9 100644
--- a/doc/developer_manual/developer_manual.tex
+++ b/doc/developer_manual/developer_manual.tex
@@ -1,1932 +1,1965 @@
\documentclass[a4paper,11pt]{article}
\usepackage{fourier}
\usepackage[T1]{fontenc}
\usepackage{microtype}
\usepackage{geometry}
\usepackage{enumitem}
\setlist[description]{leftmargin=\parindent,labelindent=\parindent}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage[utf8x]{inputenc}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{todonotes}
\usepackage{listings}
\usepackage{xspace}
\usepackage{tikz}
\usepackage{slashed}
\usepackage{subcaption}
\usetikzlibrary{arrows.meta}
\usetikzlibrary{shapes}
\usetikzlibrary{calc}
\usepackage[colorlinks,linkcolor={blue!50!black}]{hyperref}
\graphicspath{{build/figures/}{figures/}}
\usepackage[left]{showlabels}
\renewcommand{\showlabelfont}{\footnotesize\color{darkgreen}}
\emergencystretch \hsize
\newcommand{\HEJ}{{\tt HEJ}\xspace}
\newcommand{\HIGHEJ}{\emph{High Energy Jets}\xspace}
\newcommand{\cmake}{\href{https://cmake.org/}{cmake}\xspace}
\newcommand{\html}{\href{https://www.w3.org/html/}{html}\xspace}
\newcommand{\YAML}{\href{http://yaml.org/}{YAML}\xspace}
\newcommand{\QCDloop}{\href{https://github.com/scarrazza/qcdloop}{QCDloop}\xspace}
\newcommand\matel[4][]{\mathinner{\langle#2\vert#3\vert#4\rangle}_{#1}}
\newcommand{\as}{\alpha_s}
\DeclareRobustCommand{\mathgraphics}[1]{\vcenter{\hbox{\includegraphics{#1}}}}
\def\spa#1.#2{\left\langle#1\,#2\right\rangle}
\def\spb#1.#2{\left[#1\,#2\right]} \def\spaa#1.#2.#3{\langle\mskip-1mu{#1} |
#2 | {#3}\mskip-1mu\rangle} \def\spbb#1.#2.#3{[\mskip-1mu{#1} | #2 |
{#3}\mskip-1mu]} \def\spab#1.#2.#3{\langle\mskip-1mu{#1} | #2 |
{#3}\mskip-1mu\rangle} \def\spba#1.#2.#3{\langle\mskip-1mu{#1}^+ | #2 |
{#3}^+\mskip-1mu\rangle} \def\spav#1.#2.#3{\|\mskip-1mu{#1} | #2 |
{#3}\mskip-1mu\|^2} \def\jc#1.#2.#3{j^{#1}_{#2#3}}
% expectation value
\newcommand{\ev}[1]{\text{E}\left[#1\right]}
\definecolor{darkgreen}{rgb}{0,0.4,0}
\lstset{ %
backgroundcolor=\color{lightgray}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}
basicstyle=\footnotesize\usefont{T1}{DejaVuSansMono-TLF}{m}{n}, % the size of the fonts that are used for the code
breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace
breaklines=false, % sets automatic line breaking
captionpos=t, % sets the caption-position to bottom
commentstyle=\color{red}, % comment style
deletekeywords={...}, % if you want to delete keywords from the given language
escapeinside={\%*}{*)}, % if you want to add LaTeX within your code
extendedchars=true, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8
frame=false, % adds a frame around the code
keepspaces=true, % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible)
keywordstyle=\color{blue}, % keyword style
otherkeywords={}, % if you want to add more keywords to the set
numbers=none, % where to put the line-numbers; possible values are (none, left, right)
numbersep=5pt, % how far the line-numbers are from the code
rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here))
showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces'
showstringspaces=false, % underline spaces within strings only
showtabs=false, % show tabs within strings adding particular underscores
stepnumber=2, % the step between two line-numbers. If it's 1, each line will be numbered
stringstyle=\color{gray}, % string literal style
tabsize=2, % sets default tabsize to 2 spaces
title=\lstname,
emph={},
emphstyle=\color{darkgreen}
}
\begin{document}
\tikzstyle{mynode}=[rectangle split,rectangle split parts=2, draw,rectangle split part fill={lightgray, none}]
\title{\HEJ 2 developer manual}
\author{}
\maketitle
\tableofcontents
\newpage
\section{Overview}
\label{sec:overview}
\HEJ 2 is a C++ program and library implementing an algorithm to
apply \HIGHEJ resummation~\cite{Andersen:2008ue,Andersen:2008gc} to
pre-generated fixed-order events. This document is intended to give an
overview over the concepts and structure of this implementation.
\subsection{Project structure}
\label{sec:project}
\HEJ 2 is developed under the \href{https://git-scm.com/}{git}
version control system. The main repository is on the IPPP
\href{https://gitlab.com/}{gitlab} server under
\url{https://gitlab.dur.scotgrid.ac.uk/hej/hej}. To get a local
copy, get an account on the gitlab server and use
\begin{lstlisting}[language=sh,caption={}]
git clone git@gitlab.dur.scotgrid.ac.uk:hej/hej.git
\end{lstlisting}
This should create a directory \texttt{hej} with the following
contents:
\begin{description}
\item[README.md:] Basic information concerning \HEJ.
\item[doc:] Contains additional documentation, see section~\ref{sec:doc}.
\item[include:] Contains the C++ header files.
\item[src:] Contains the C++ source files.
\item[examples:] Example analyses and scale setting code.
\item[t:] Contains the source code for the automated tests.
\item[CMakeLists.txt:] Configuration file for the \cmake build
system. See section~\ref{sec:cmake}.
\item[cmake:] Auxiliary files for \cmake. This includes modules for
finding installed software in \texttt{cmake/Modules} and templates for
code generation during the build process in \texttt{cmake/Templates}.
\item[config.yml:] Sample configuration file for running \HEJ 2.
\item[current\_generator:] Contains the code for the current generator,
see section~\ref{sec:cur_gen}.
\item[FixedOrderGen:] Contains the code for the fixed-order generator,
see section~\ref{sec:HEJFOG}.
\item[COPYING:] License file.
\item[AUTHORS:] List of \HEJ authors.
\item[Changes-API.md:] Change history for the \HEJ API (application programming interface).
\item[Changes.md:] History of changes relevant for \HEJ users.
\end{description}
In the following all paths are given relative to the
\texttt{hej} directory.
\subsection{Documentation}
\label{sec:doc}
The \texttt{doc} directory contains user documentation in
\texttt{doc/sphinx} and the configuration to generate source code
documentation in \texttt{doc/doxygen}.
The user documentation explains how to install and run \HEJ 2. The
format is
\href{http://docutils.sourceforge.net/rst.html}{reStructuredText}, which
is mostly human-readable. Other formats, like \html, can be generated with the
\href{http://www.sphinx-doc.org/en/master/}{sphinx} generator with
\begin{lstlisting}[language=sh,caption={}]
make sphinx
\end{lstlisting}
To document the source code we use
\href{https://www.stack.nl/~dimitri/doxygen/}{doxygen}. To generate
\html documentation, use the command
\begin{lstlisting}[language=sh,caption={}]
make doxygen
\end{lstlisting}
in your \texttt{build/} directory.
\subsection{Build system}
\label{sec:cmake}
For the most part, \HEJ 2 is a library providing classes and
functions that can be used to add resummation to fixed-order events. In
addition, there is a relatively small executable program leveraging this
library to read in events from an input file and produce resummation
events. Both the library and the program are built and installed with
the help of \cmake.
Debug information can be turned on by using
\begin{lstlisting}[language=sh,caption={}]
cmake base/directory -DCMAKE_BUILD_TYPE=Debug
make install
\end{lstlisting}
This facilitates the use of debuggers like \href{https://www.gnu.org/software/gdb/}{gdb}.
The main \cmake configuration file is \texttt{CMakeLists.txt}. It defines the
compiler flags, software prerequisites, header and source files used to
build \HEJ 2, and the automated tests.
\texttt{cmake/Modules} contains module files that help with the
detection of the software prerequisites and \texttt{cmake/Templates}
template files for the automatic generation of header and
source files. For example, this allows to only keep the version
information in one central location (\texttt{CMakeLists.txt}) and
automatically generate a header file from the template \texttt{Version.hh.in} to propagate this to the C++ code.
\subsection{General coding guidelines}
\label{sec:notes}
The goal is to make the \HEJ 2 code well-structured and
readable. Here are a number of guidelines to this end.
\begin{description}
\item[Observe the boy scout rule.] Always leave the code cleaner
than how you found it. Ugly hacks can be useful for testing, but
shouldn't make their way into the main branch.
\item[Ask if something is unclear.] Often there is a good reason why
code is written the way it is. Sometimes that reason is only obvious to
the original author (use \lstinline!git blame! to find them), in which
case they should be poked to add a comment. Sometimes there is no good
reason, but nobody has had the time to come up with something better,
yet. In some places the code might just be bad.
\item[Don't break tests.] There are a number of tests in the \texttt{t}
directory, which can be run with \lstinline!make test!. Ideally, all
tests should run successfully in each git revision. If your latest
commit broke a test and you haven't pushed to the central repository
yet, you can fix it with \lstinline!git commit --amend!. If an earlier
local commit broke a test, you can use \lstinline!git rebase -i! if
you feel confident. Additionally each \lstinline!git push! is also
automatically tested via the GitLab CI (see appendix~\ref{sec:CI}).
\item[Test your new code.] When you add some new functionality, also add an
automated test. This can be useful even if you don't know the
``correct'' result because it prevents the code from changing its behaviour
silently in the future. \href{http://www.valgrind.org/}{valgrind} is a
very useful tool to detect potential memory leaks. The code coverage of all
tests can be generated with \href{https://gcovr.com/en/stable/}{gcovr}.
Therefore add the flag \lstinline!-DTEST_COVERAGE=True! to cmake and run
\lstinline!make ctest_coverage!.
\item[Stick to the coding style.] It is somewhat easier to read code
that has a uniform coding and indentation style. We don't have a
strict style, but it helps if your code looks similar to what is
already there.
\end{description}
\section{Program flow}
\label{sec:flow}
A run of the \HEJ 2 program has three stages: initialisation,
event processing, and cleanup. The following sections outline these
stages and their relations to the various classes and functions in the
code. Unless denoted otherwise, all classes and functions are part of
the \lstinline!HEJ! namespace. The code for the \HEJ 2 program is
in \texttt{src/bin/HEJ.cc}, all other code comprises the \HEJ 2
library. Classes and free functions are usually implemented in header
and source files with a corresponding name, i.e. the code for
\lstinline!MyClass! can usually be found in
\texttt{include/HEJ/MyClass.hh} and \texttt{src/MyClass.cc}.
\subsection{Initialisation}
\label{sec:init}
The first step is to load and parse the \YAML configuration file. The
entry point for this is the \lstinline!load_config! function and the
related code can be found in \texttt{include/HEJ/YAMLreader.hh},
\texttt{include/HEJ/config.hh} and the corresponding \texttt{.cc} files
in the \texttt{src} directory. The implementation is based on the
\href{https://github.com/jbeder/yaml-cpp}{yaml-cpp} library.
The \lstinline!load_config! function returns a \lstinline!Config! object
containing all settings. To detect potential mistakes as early as
possible, we throw an exception whenever one of the following errors
occurs:
\begin{itemize}
\item There is an unknown option in the \YAML file.
\item A setting is invalid, for example a string is given where a number
would be expected.
\item An option value is not set.
\end{itemize}
The third rule is sometimes relaxed for ``advanced'' settings with an
obvious default, like for importing custom scales or analyses.
The information stored in the \lstinline!Config! object is then used to
initialise various objects required for the event processing stage described in
section~\ref{sec:processing}. First, the \lstinline!get_analyses! function
creates a number objects that inherits from the \lstinline!Analysis!
interface.\footnote{In the context of C++ the proper technical expression is
``pure abstract class''.} Using an interface allows us to decide the concrete
type of the analysis at run time instead of having to make a compile-time
decision. Depending on the settings, \lstinline!get_analyses! creates a number
of user-defined analyses loaded from one or more external libraries and
\textsc{rivet} analyses (see the user documentation
\url{https://hej.web.cern.ch/HEJ/doc/current/user/})
Together with a number of further objects, whose roles are described in
section~\ref{sec:processing}, we also initialise the global random
number generator. We again use an interface to defer deciding the
concrete type until the program is actually run. Currently, we support the
\href{https://mixmax.hepforge.org/}{MIXMAX}
(\texttt{include/HEJ/Mixmax.hh}) and Ranlux64
(\texttt{include/HEJ/Ranlux64.hh}) random number generators, both are provided
by \href{http://proj-clhep.web.cern.ch/}{CLHEP}.
We also set up a \lstinline!HEJ::EventReader! object for reading events
either in the the Les Houches event file format~\cite{Alwall:2006yp} or
an \href{https://www.hdfgroup.org/}{HDF5}-based
format~\cite{Hoeche:2019rti}. To allow making the decision at run time,
\lstinline!HEJ::EventReader! is an abstract base class defined in
\texttt{include/HEJ/EventReader.hh} and the implementations of the
derived classes are in \texttt{include/HEJ/LesHouchesReader.hh},
\texttt{include/HEJ/HDF5Reader.hh} and the corresponding \texttt{.cc}
source files in the \texttt{src} directory. The
\lstinline!LesHouchesReader! leverages
\href{http://home.thep.lu.se/~leif/LHEF/}{\texttt{include/LHEF/LHEF.h}}. A
small wrapper around the
\href{https://www.boost.org/doc/libs/1_67_0/libs/iostreams/doc/index.html}{boost
iostreams} library allows us to also read event files compressed with
\href{https://www.gnu.org/software/gzip/}{gzip}. The wrapper code is in
\texttt{include/HEJ/stream.hh} and the \texttt{src/stream.cc}.
If unweighting is enabled, we also initialise an unweighter as defined
in \texttt{include/HEJ/Unweighter.hh}. The unweighting strategies are
explained in section~\ref{sec:unweight}.
\subsection{Event processing}
\label{sec:processing}
In the second stage events are continuously read from the event
file. After jet clustering, a number of corresponding resummation events
are generated for each input event and fed into the analyses and a
number of output files. The roles of various classes and functions are
illustrated in the following flow chart:
\begin{center}
\begin{tikzpicture}[node distance=2cm and 5mm]
\node (reader) [mynode]
{\lstinline!EventReader::read_event!\nodepart{second}{read event}};
\node
(data) [mynode,below=of reader]
{\lstinline!Event::EventData! constructor\nodepart{second}{convert to \HEJ object}};
\node
(cluster) [mynode,below=of data]
{\lstinline!Event::EventData::cluster!\nodepart{second}{cluster jets \&
classify \lstinline!EventType!}};
\node
(resum) [mynode,below=of cluster]
{\lstinline!EventReweighter::reweight!\nodepart{second}{perform resummation}};
\node
(cut) [mynode,below=of resum]
{\lstinline!Analysis::pass_cuts!\nodepart{second}{apply cuts}};
\node
(cut) [mynode,below=of resum]
{\lstinline!Analysis::pass_cuts!\nodepart{second}{apply cuts}};
\node
(unweight) [mynode,below=of cut]
{\lstinline!Unweighter::unweight!\nodepart{second}{unweight (optional)}};
\node
(fill) [mynode,below left=of unweight]
{\lstinline!Analysis::fill!\nodepart{second}{analyse event}};
\node
(write) [mynode,below right=of unweight]
{\lstinline!CombinedEventWriter::write!\nodepart{second}{write out event}};
\node
(control) [below=of unweight] {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(reader.south) -- node[left] {\lstinline!LHEF::HEPEUP!} (data.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(data.south) -- node[left] {\lstinline!Event::EventData!} (cluster.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(cluster.south) -- node[left] {\lstinline!Event!} (resum.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(resum.south) -- (cut.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(resum.south)+(10mm, 0cm)$) -- ($(cut.north)+(10mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(resum.south)+(5mm, 0cm)$) -- ($(cut.north)+(5mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(resum.south)-(5mm, 0cm)$) -- ($(cut.north)-(5mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(resum.south)-(10mm, 0cm)$) -- node[left] {\lstinline!Event!} ($(cut.north)-(10mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(cut.south) -- (unweight.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)+(7mm, 0cm)$) -- ($(unweight.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cut.south)-(7mm, 0cm)$) -- node[left] {\lstinline!Event!} ($(unweight.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(unweight.south)-(3mm,0mm)$) .. controls ($(control)-(3mm,0mm)$) ..node[left] {\lstinline!Event!} (fill.east);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(unweight.south)-(3mm,0mm)$) .. controls ($(control)-(3mm,0mm)$) .. (write.west);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(unweight.south)+(3mm,0mm)$) .. controls ($(control)+(3mm,0mm)$) .. (fill.east);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(unweight.south)+(3mm,0mm)$) .. controls ($(control)+(3mm,0mm)$) ..node[right] {\lstinline!Event!} (write.west);
\end{tikzpicture}
\end{center}
\lstinline!EventData! is an intermediate container, its members are completely
accessible. In contrast after jet clustering and classification the phase space
inside \lstinline!Event! can not be changed any more
(\href{https://wikipedia.org/wiki/Builder_pattern}{Builder design pattern}). The
resummation is performed by the \lstinline!EventReweighter! class, which is
described in more detail in section~\ref{sec:resum}. The
\lstinline!CombinedEventWriter! writes events to zero or more output files. To
this end, it contains a number of objects implementing the
\lstinline!EventWriter! interface. These event writers typically write the
events to a file in a given format. We currently have the
\lstinline!LesHouchesWriter! for event files in the Les Houches Event File
format, the \lstinline!HDF5Writer! for
\href{https://www.hdfgroup.org/}{HDF5}~\cite{Hoeche:2019rti} and the
\lstinline!HepMC2Writer! or \lstinline!HepMC3Writer! for the
\href{https://hepmc.web.cern.ch/hepmc/}{HepMC} format (Version 2 and
3).
\subsection{Resummation}
\label{sec:resum}
In the \lstinline!EventReweighter::reweight! member function, we first
classify the input fixed-order event (FKL, unordered, non-resummable, \dots)
and decide according to the user settings whether to discard, keep, or
resum the event. If we perform resummation for the given event, we
generate a number of trial \lstinline!PhaseSpacePoint! objects. Phase
space generation is discussed in more detail in
section~\ref{sec:pspgen}. We then perform jet clustering according to
the settings for the resummation jets on each
\lstinline!PhaseSpacePoint!, update the factorisation and
renormalisation scale in the resulting \lstinline!Event! and reweight it
according to the ratio of pdf factors and \HEJ matrix elements between
resummation and original fixed-order event:
\begin{center}
\begin{tikzpicture}[node distance=1.5cm and 5mm]
\node (in) {};
\node (treat) [diamond,draw,below=of in,minimum size=3.5cm,
label={[anchor=west, inner sep=8pt]west:discard},
label={[anchor=east, inner sep=14pt]east:keep},
label={[anchor=south, inner sep=20pt]south:reweight}
] {};
\draw (treat.north west) -- (treat.south east);
\draw (treat.north east) -- (treat.south west);
\node
(psp) [mynode,below=of treat]
{\lstinline!PhaseSpacePoint! constructor};
\node
(cluster) [mynode,below=of psp]
{\lstinline!Event::EventData::cluster!\nodepart{second}{cluster jets}};
\node
(colour) [mynode,below=of cluster]
{\lstinline!Event::generate_colours()!\nodepart{second}{generate particle colour}};
\node
(gen_scales) [mynode,below=of colour]
{\lstinline!ScaleGenerator::operator()!\nodepart{second}{update scales}};
\node
(rescale) [mynode,below=of gen_scales]
{\lstinline!PDF::pdfpt!,
\lstinline!MatrixElement!\nodepart{second}{reweight}};
\node (out) [below of=rescale] {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(in.south) -- node[left] {\lstinline!Event!} (treat.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(treat.south) -- node[left] {\lstinline!Event!} (psp.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(psp.south) -- (cluster.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(psp.south)+(7mm, 0cm)$) -- ($(cluster.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(psp.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!PhaseSpacePoint!} ($(cluster.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(cluster.south) -- (colour.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cluster.south)+(7mm, 0cm)$) -- ($(colour.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(cluster.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(colour.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(colour.south) -- (gen_scales.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(colour.south)+(7mm, 0cm)$) -- ($(gen_scales.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(colour.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(gen_scales.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(gen_scales.south) -- (rescale.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(gen_scales.south)+(7mm, 0cm)$) -- ($(rescale.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(gen_scales.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(rescale.north)-(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
(rescale.south) -- (out.north);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(rescale.south)+(7mm, 0cm)$) -- ($(out.north)+(7mm, 0cm)$);
\draw[-{Latex[length=3mm, width=1.5mm]}]
($(rescale.south)-(7mm, 0cm)$) -- node[left]
{\lstinline!Event!} ($(out.north)-(7mm, 0cm)$);
\node (helper) at ($(treat.east) + (15mm,0cm)$) {};
\draw[-{Latex[length=3mm, width=1.5mm]}]
(treat.east) -- ($(treat.east) + (15mm,0cm)$)
-- node[left] {\lstinline!Event!} (helper |- gen_scales.east) -- (gen_scales.east)
;
\end{tikzpicture}
\end{center}
\subsection{Phase space point generation}
\label{sec:pspgen}
The resummed and matched \HEJ cross section for pure jet production of
FKL configurations is given by (cf. eq. (3) of~\cite{Andersen:2018tnm})
\begin{align}
\label{eq:resumdijetFKLmatched2}
% \begin{split}
\sigma&_{2j}^\mathrm{resum, match}=\sum_{f_1, f_2}\ \sum_m
\prod_{j=1}^m\left(
\int_{p_{j\perp}^B=0}^{p_{j\perp}^B=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{j\perp}^B}{(2\pi)^3}\ \int
\frac{\mathrm{d} y_j^B}{2} \right) \
(2\pi)^4\ \delta^{(2)}\!\!\left(\sum_{k=1}^{m}
\mathbf{p}_{k\perp}^B\right)\nonumber\\
&\times\ x_a^B\ f_{a, f_1}(x_a^B, Q_a^B)\ x_b^B\ f_{b, f_2}(x_b^B, Q_b^B)\
\frac{\overline{\left|\mathcal{M}_\text{LO}^{f_1f_2\to f_1g\cdots
gf_2}\big(\big\{p^B_j\big\}\big)\right|}^2}{(\hat {s}^B)^2}\nonumber\\
& \times (2\pi)^{-4+3m}\ 2^m \nonumber\\
&\times\ \sum_{n=2}^\infty\
\int_{p_{1\perp}=p_{\perp,\mathrm{min}} }^{p_{1\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{1\perp}}{(2\pi)^3}\
\int_{p_{n\perp}=p_{\perp,\mathrm{min}}}^{p_{n\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{n\perp}}{(2\pi)^3}\
\prod_{i=2}^{n-1}\int_{p_{i\perp}=\lambda}^{p_{i\perp}=\infty}
\frac{\mathrm{d}^2\mathbf{p}_{i\perp}}{(2\pi)^3}\ (2\pi)^4\ \delta^{(2)}\!\!\left(\sum_{k=1}^n
\mathbf{p}_{k\perp}\right )\\
&\times \ \mathbf{T}_y \prod_{i=1}^n
\left(\int \frac{\mathrm{d} y_i}{2}\right)\
\mathcal{O}_{mj}^e\
\left(\prod_{l=1}^{m-1}\delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l}\perp}^B -
\mathbf{j}_{l\perp})\right)\
\left(\prod_{l=1}^m\delta(y^B_{\mathcal{J}_l}-y_{\mathcal{J}_l})\right)
\ \mathcal{O}_{2j}(\{p_i\})\nonumber\\
&\times \frac{(\hat{s}^B)^2}{\hat{s}^2}\ \frac{x_a f_{a,f_1}(x_a, Q_a)\ x_b f_{b,f_2}(x_b, Q_b)}{x_a^B\ f_{a,f_1}(x_a^B, Q_a^B)\ x_b^B\ f_{b,f_2}(x_b^B, Q_b^B)}\ \frac{\overline{\left|\mathcal{M}_{\mathrm{HEJ}}^{f_1 f_2\to f_1 g\cdots
gf_2}(\{ p_i\})\right|}^2}{\overline{\left|\mathcal{M}_\text{LO, HEJ}^{f_1f_2\to f_1g\cdots
gf_2}\big(\big\{p^B_j\big\}\big)\right|}^{2}} \,.\nonumber
% \end{split}
\end{align}
The first two lines correspond to the generation of the fixed-order
input events with incoming partons $f_1, f_2$ and outgoing momenta
$p_j^B$, where $\mathbf{p}_{j\perp}^B$ and $y_j^B$ denote the respective
transverse momentum and rapidity. Note that, at leading order, these
coincide with the fixed-order jet momenta $p_{\mathcal{J}_j}^B$.
$f_{a,f_1}(x_a, Q_a),f_{b,f_2}(x_b, Q_b)$ are the pdf factors for the incoming partons with
momentum fractions $x_a$ and $x_b$. The square of the partonic
centre-of-mass energy is denoted by $\hat{s}^B$ and
$\mathcal{M}_\text{LO}^{f_1f_2\to f_1g\cdots gf_2}$ is the
leading-order matrix element.
The third line is a factor accounting for the different multiplicities
between fixed-order and resummation events. Lines four and five are
the integration over the resummation phase space described in this
section. $p_i$ are the momenta of the outgoing partons in resummation
phase space. $\mathbf{T}_y$ denotes rapidity
ordering and $\mathcal{O}_{mj}^e$ projects out the exclusive $m$-jet
component. The relation between resummation and fixed-order momenta is
fixed by the $\delta$ functions. The first sets each transverse fixed-order jet
momentum to some function $\mathbf{j_{l\perp}}$ of the resummation
momenta. The exact form is described in section~\ref{sec:ptj_res}. The second
$\delta$ forces the rapidities of resummation and fixed-order jets to be
the same. Finally, the last line is the reweighting of pdf and matrix
element factors already shown in section~\ref{sec:resum}.
There are two kinds of cut-off in the integration over the resummation
partons. $\lambda$ is a technical cut-off connected to the cancellation
of infrared divergencies between real and virtual corrections. Its
numerical value is set in
\texttt{include/HEJ/Constants.h}. $p_{\perp,\mathrm{min}}$ regulates
and \emph{uncancelled} divergence in the extremal parton momenta. Its
size is set by the user configuration \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJ.html#settings}.
It is straightforward to generalise eq.~(\ref{eq:resumdijetFKLmatched2})
to unordered configurations and processes with additional colourless
emissions, for example a Higgs or electroweak boson. In the latter case only
the fixed-order integration and the matrix elements change.
\subsubsection{Gluon Multiplicity}
\label{sec:psp_ng}
The first step in evaluating the resummation phase space in
eq.~(\ref{eq:resumdijetFKLmatched2}) is to randomly pick terms in the
sum over the number of emissions. This sampling of the gluon
multiplicity is done in the \lstinline!PhaseSpacePoint::sample_ng!
function in \texttt{src/PhaseSpacePoint.cc}.
The typical number of extra emissions depends strongly on the rapidity
span of the underlying fixed-order event. Let us, for example, consider
a fixed-order FKL-type multi-jet configuration with rapidities
$y_{j_f},\,y_{j_b}$ of the most forward and backward jets,
respectively. By eq.~(\ref{eq:resumdijetFKLmatched2}), the jet
multiplicity and the rapidity of each jet are conserved when adding
resummation. This implies that additional hard radiation is restricted
to rapidities $y$ within a region $y_{j_b} \lesssim y \lesssim
y_{j_f}$. Within \HEJ, we require the most forward and most backward
emissions to be hard \todo{specify how hard} in order to avoid divergences, so this constraint
in fact applies to \emph{all} additional radiation.
To simplify the remaining discussion, let us remove the FKL rapidity
ordering
\begin{equation}
\label{eq:remove_y_order}
\mathbf{T}_y \prod_{i=1}^n\int \frac{\mathrm{d}y_i}{2} =
\frac{1}{n!}\prod_{i=1}^n\int
\frac{\mathrm{d}y_i}{2}\,,
\end{equation}
where all rapidity integrals now cover a region which is approximately
bounded by $y_{j_b}$ and $y_{j_f}$. Each of the $m$ jets has to contain at least
one parton; selecting random emissions we can rewrite the phase space
integrals as
\begin{equation}
\label{eq:select_jets}
\frac{1}{n!}\prod_{i=1}^n\int [\mathrm{d}p_i] =
\left(\prod_{i=1}^{m}\int [\mathrm{d}p_i]\ {\cal J}_i(p_i)\right)
\frac{1}{n_g!}\prod_{i=m+1}^{m+n_g}\int [\mathrm{d}p_i]
\end{equation}
with jet selection functions
\begin{equation}
\label{eq:def_jet_selection}
{\cal J}_i(p) =
\begin{cases}
1 &p\text{ clustered into jet }i\\
0 & \text{otherwise}
\end{cases}
\end{equation}
and $n_g \equiv n - m$. Here and in the following we use the short-hand
notation $[\mathrm{d}p_i]$ to denote the phase-space measure for parton
$i$. As is evident from eq.~\eqref{eq:select_jets}, adding an extra emission
$n_g+1$ introduces a suppression factor $\tfrac{1}{n_g+1}$. However, the
additional phase space integral also results in an enhancement proportional
to $\Delta y_{j_f j_b} = y_{j_f} - y_{j_b}$. This is a result of the
rapidity-independence of the MRK limit of the integrand, consisting of the
matrix elements divided by the flux factor. Indeed, we observe that the
typical number of gluon emissions is to a good approximation proportional to
the rapidity separation and the phase space integral is dominated by events
with $n_g \approx \Delta y_{j_f j_b}$.
For the actual phase space sampling, we assume a Poisson distribution
and extract the mean number of gluon emissions in different rapidity
bins and fit the results to a linear function in $\Delta y_{j_f j_b}$,
finding a coefficient of $0.975$ for the inclusive production of a Higgs
boson with two jets. Here are the observed and fitted average gluon
multiplicities as a function of $\Delta y_{j_f j_b}$:
\begin{center}
\includegraphics[width=.75\textwidth]{ng_mean}
\end{center}
As shown for two rapidity slices the assumption of a Poisson
distribution is also a good approximation:
\begin{center}
\includegraphics[width=.49\textwidth]{{ng_1.5}.pdf}\hfill
\includegraphics[width=.49\textwidth]{{ng_5.5}.pdf}
\end{center}
\subsubsection{Number of Gluons inside Jets}
\label{sec:psp_ng_jet}
For each of the $n_g$ gluon emissions we can split the phase-space
integral into a (disconnected) region inside the jets and a remainder:
\begin{equation}
\label{eq:psp_split}
\int [\mathrm{d}p_i] = \int [\mathrm{d}p_i]\,
\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg) + \int [\mathrm{d}p_i]\,
\bigg[1-\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg)\bigg]\,.
\end{equation}
The next step is to decide how many of the gluons will form part of a
jet. This is done in the \lstinline!PhaseSpacePoint::sample_ng_jets!
function.
We choose an importance sampling which is flat in the plane
spanned by the azimuthal angle $\phi$ and the rapidity $y$. This is
observed in BFKL and valid in the limit of Multi-Regge-Kinematics
(MRK). Furthermore, we assume anti-$k_t$ jets, which cover an area of
$\pi R^2$.
In principle, the total accessible area in the $y$-$\phi$ plane is given
by $2\pi \Delta y_{fb}$, where $\Delta y_{fb}\geq \Delta y_{j_f j_b}$ is
the a priori unknown rapidity separation between the most forward and
backward partons. In most cases the extremal jets consist of single
partons, so that $\Delta y_{fb} = \Delta y_{j_f j_b}$. For the less common
case of two partons forming a jet we observe a maximum distance of $R$
between the constituents and the jet centre. In rare cases jets have
more than two constituents. Empirically, they are always within a
distance of $\tfrac{5}{3}R$ to the centre of the jet, so
$\Delta y_{fb} \leq \Delta y_{j_f j_b} + \tfrac{10}{3} R$. In practice, the
extremal partons are required to carry a large fraction of the jet
transverse momentum and will therefore be much closer to the jet axis.
In summary, for sufficiently large rapidity separations we can use the
approximation $\Delta y_{fb} \approx \Delta y_{j_f j_b}$. This scenario
is depicted here:
\begin{center}
\includegraphics[width=0.5\linewidth]{ps_large_y}
\end{center}
If there is no overlap between jets, the probability $p_{\cal J, >}$ for
an extra gluon to end up inside a jet is then given by
\begin{equation}
\label{eq:p_J_large}
p_{\cal J, >} = \frac{(m - 1)\*R^2}{2\Delta y_{j_f j_b}}\,.
\end{equation}
For a very small rapidity separation, eq.~\eqref{eq:p_J_large}
obviously overestimates the true probability. The maximum phase space
covered by jets in the limit of a vanishing rapidity distance between
all partons is $2mR \Delta y_{fb}$:
\begin{center}
\includegraphics[width=0.5\linewidth]{ps_small_y}
\end{center}
We therefore estimate the probability for a parton to end up inside a jet as
\begin{equation}
\label{eq:p_J}
p_{\cal J} = \min\bigg(\frac{(m - 1)\*R^2}{2\Delta y_{j_f j_b}}, \frac{mR}{\pi}\bigg)\,.
\end{equation}
Here we compare this estimate with the actually observed
fraction of additional emissions into jets as a function of the rapidity
separation:
\begin{center}
\includegraphics[width=0.75\linewidth]{pJ}
\end{center}
\subsubsection{Gluons outside Jets}
\label{sec:gluons_nonjet}
Using our estimate for the probability of a gluon to be a jet
constituent, we choose a number $n_{g,{\cal J}}$ of gluons inside
jets, which also fixes the number $n_g - n_{g,{\cal J}}$ of gluons
outside jets. As explained later on, we need to generate the momenta of
the gluons outside jets first. This is done in
\lstinline!PhaseSpacePoint::gen_non_jet!.
The azimuthal angle $\phi$ is generated flat within $0\leq \phi \leq 2
\pi$. The allowed rapidity interval is set by the most forward and
backward partons, which are necessarily inside jets. Since these parton
rapidities are not known at this point, we also have to postpone the
rapidity generation for the gluons outside jets. For the scalar
transverse momentum $p_\perp = |\mathbf{p}_\perp|$ of a gluon outside
jets we use the parametrisation
\begin{equation}
\label{eq:p_nonjet}
p_\perp = \lambda + \tilde{p}_\perp\*\tan(\tau\*r)\,, \qquad
\tau = \arctan\bigg(\frac{p_{\perp{\cal J}_\text{min}} - \lambda}{\tilde{p}_\perp}\bigg)\,.
\end{equation}
For $r \in [0,1)$, $p_\perp$ is always less than the minimum momentum
$p_{\perp{\cal J}_\text{min}}$ required for a jet. $\tilde{p}_\perp$ is
a free parameter, a good empirical value is $\tilde{p}_\perp = [1.3 +
0.2\*(n_g - n_{g,\cal J})]\,$GeV
\subsubsection{Resummation jet momenta}
\label{sec:ptj_res}
On the one hand, each jet momentum is given by the sum of its
constituent momenta. On the other hand, the resummation jet momenta are
fixed by the constraints in line five of the master
equation~\eqref{eq:resumdijetFKLmatched2}. We therefore have to
calculate the resummation jet momenta from these constraints before
generating the momenta of the gluons inside jets. This is done in
\lstinline!PhaseSpacePoint::reshuffle! and in the free
\lstinline!resummation_jet_momenta! function (declared in \texttt{resummation\_jet.hh}).
The resummation jet momenta are determined by the $\delta$ functions in
line five of eq.~(\ref{eq:resumdijetFKLmatched2}). The rapidities are
fixed to the rapidities of the jets in the input fixed-order events, so
that the FKL ordering is guaranteed to be preserved.
In traditional \HEJ reshuffling the transverse momentum are given through
\begin{equation}
\label{eq:ptreassign_old}
\mathbf{p}^B_{\mathcal{J}_{l\perp}} = \mathbf{j}_{l\perp} \equiv \mathbf{p}_{\mathcal{J}_{l}\perp}
+ \mathbf{q}_\perp \,\frac{|\mathbf{p}_{\mathcal{J}_{l}\perp}|}{P_\perp},
\end{equation}
where $\mathbf{q}_\perp = \sum_{j=1}^n \mathbf{p}_{i\perp}
\bigg[1-\theta\bigg(\sum_{j=1}^{m}{\cal J}_j(p_i)\bigg)\bigg] $ is the
total transverse momentum of all partons \emph{outside} jets and
$P_\perp = \sum_{j=1}^m |\mathbf{p}_{\mathcal{J}_{j}\perp}|$. Since the
total transverse momentum of an event vanishes, we can also use
$\mathbf{q}_\perp = - \sum_{j=1}^m
\mathbf{p}_{\mathcal{J}_{j}\perp}$. Eq.~(\ref{eq:ptreassign}) is a
non-linear system of equations in the resummation jet momenta
$\mathbf{p}_{\mathcal{J}_{l}\perp}$. Hence we would have to solve
\begin{equation}
\label{eq:ptreassign_eq}
\mathbf{p}_{\mathcal{J}_{l}\perp}=\mathbf{j}^B_{l\perp} \equiv\mathbf{j}_{l\perp}^{-1}
\left(\mathbf{p}^B_{\mathcal{J}_{l\perp}}\right)
\end{equation}
numerically.
Since solving such a system is computationally expensive, we instead
change the reshuffling around to be linear in the resummation jet
momenta. Hence~\eqref{eq:ptreassign_eq} gets replaces by
\begin{equation}
\label{eq:ptreassign}
\mathbf{p}_{\mathcal{J}_{l\perp}} = \mathbf{j}^B_{l\perp} \equiv \mathbf{p}^B_{\mathcal{J}_{l}\perp}
- \mathbf{q}_\perp \,\frac{|\mathbf{p}^B_{\mathcal{J}_{l}\perp}|}{P^B_\perp},
\end{equation}
which is linear in the resummation momentum. Consequently the equivalent
of~\eqref{eq:ptreassign_old} is non-linear in the Born momentum. However
the exact form of~\eqref{eq:ptreassign_old} is not relevant for the resummation.
Both methods have been tested for two and three jets with the \textsc{rivet}
standard analysis \texttt{MC\_JETS}. They didn't show any differences even
after $10^9$ events.
The reshuffling relation~\eqref{eq:ptreassign} allows the transverse
momenta $p^B_{\mathcal{J}_{l\perp}}$ of the fixed-order jets to be
somewhat below the minimum transverse momentum of resummation jets. It
is crucial that this difference does not become too large, as the
fixed-order cross section diverges for vanishing transverse momenta. In
the production of a Higgs boson with resummation jets above $30\,$GeV we observe
that the contribution from fixed-order events with jets softer than
about $20\,$GeV can be safely neglected. This is shown in the following
plot of the differential cross section over the transverse momentum of
the softest fixed-order jet:
\begin{center}
\includegraphics[width=.75\textwidth]{ptBMin}
\end{center}
Finally, we have to account for the fact that the reshuffling
relation~\eqref{eq:ptreassign} is non-linear in the Born momenta. To
arrive at the master formula~\eqref{eq:resumdijetFKLmatched2} for the
cross section, we have introduced unity in the form of an integral over
the Born momenta with $\delta$ functions in the integrand, that is
\begin{equation}
\label{eq:delta_intro}
1 = \int_{p_{j\perp}^B=0}^{p_{j\perp}^B=\infty}
\mathrm{d}^2\mathbf{p}_{j\perp}^B\delta^{(2)}(\mathbf{p}_{\mathcal{J}_{j\perp}}^B -
\mathbf{j}_{j\perp})\,.
\end{equation}
If the arguments of the $\delta$ functions are not linear in the Born
momenta, we have to compensate with additional Jacobians as
factors. Explicitly, for the reshuffling relation~\eqref{eq:ptreassign}
we have
\begin{equation}
\label{eq:delta_rewrite}
\prod_{l=1}^m \delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l\perp}}^B -
\mathbf{j}_{l\perp}) = \Delta \prod_{l=1}^m \delta^{(2)}(\mathbf{p}_{\mathcal{J}_{l\perp}} -
\mathbf{j}_{l\perp}^B)\,,
\end{equation}
where $\mathbf{j}_{l\perp}^B$ is given by~\eqref{eq:ptreassign_eq} and only
depends on the Born momenta. We have extended the product to run to $m$
instead of $m-1$ by eliminating the last $\delta$ function
$\delta^{(2)}\!\!\left(\sum_{k=1}^n \mathbf{p}_{k\perp}\right )$.
The Jacobian $\Delta$ is the determinant of a $2m \times 2m$ matrix with $l, l' = 1,\dots,m$
and $X, X' = x,y$.
\begin{equation}
\label{eq:jacobian}
\Delta = \left|\frac{\partial\,\mathbf{j}^B_{l'\perp}}{\partial\, \mathbf{p}^B_{{\cal J}_l \perp}} \right|
= \left| \delta_{l l'} \delta_{X X'} - \frac{q_X\, p^B_{{\cal
J}_{l'}X'}}{\left|\mathbf{p}^B_{{\cal J}_{l'} \perp}\right| P^B_\perp}\left(\delta_{l l'}
- \frac{\left|\mathbf{p}^B_{{\cal J}_l \perp}\right|}{P^B_\perp}\right)\right|\,.
\end{equation}
The determinant is calculated in \lstinline!resummation_jet_weight!,
again coming from the \texttt{resummation\_jet.hh} header.
Having to introduce this Jacobian is not a disadvantage specific to the new
reshuffling. If we instead use the old reshuffling
relation~\eqref{eq:ptreassign_old} we \emph{also} have to introduce a
similar Jacobian since we actually want to integrate over the
resummation phase space and need to transform the argument of the
$\delta$ function to be linear in the resummation momenta for this.
\subsubsection{Gluons inside Jets}
\label{sec:gluons_jet}
After the steps outlined in section~\ref{sec:psp_ng_jet}, we have a
total number of $m + n_{g,{\cal J}}$ constituents. In
\lstinline!PhaseSpacePoint::distribute_jet_partons! we distribute them
randomly among the jets such that each jet has at least one
constituent. We then generate their momenta in
\lstinline!PhaseSpacePoint::split! using the \lstinline!Splitter! class.
The phase space integral for a jet ${\cal J}$ is given by
\begin{equation}
\label{eq:ps_jetparton} \prod_{i\text{ in }{\cal J}} \bigg(\int
\mathrm{d}\mathbf{p}_{i\perp}\ \int \mathrm{d} y_i
\bigg)\delta^{(2)}\Big(\sum_{i\text{ in }{\cal J}} \mathbf{p}_{i\perp} -
\mathbf{j}_{\perp}^B\Big)\delta(y_{\mathcal{J}}-y^B_{\mathcal{J}})\,.
\end{equation}
For jets with a single constituent, the parton momentum is obiously equal to the
jet momentum. In the case of two constituents, we observe that the
partons are always inside the jet cone with radius $R$ and often very
close to the jet centre. The following plots show the typical relative
distance $\Delta R/R$ for this scenario:
\begin{center}
\includegraphics[width=0.45\linewidth]{dR_2}
\includegraphics[width=0.45\linewidth]{dR_2_small}
\end{center}
According to this preference for small values of $\Delta R$, we
parametrise the $\Delta R$ integrals as
\begin{equation}
\label{eq:dR_sampling}
\frac{\Delta R}{R} =
\begin{cases}
0.25\,x_R & x_R < 0.4 \\
1.5\,x_R - 0.5 & x_R \geq 0.4
\end{cases}\,.
\end{equation}
Next, we generate $\Theta_1 \equiv \Theta$ and use the constraint $\Theta_2 = \Theta
\pm \pi$. The transverse momentum of the first parton is then given by
\begin{equation}
\label{eq:delta_constraints}
p_{1\perp} =
\frac{p_{\mathcal{J} y} - \tan(\phi_2) p_{\mathcal{J} x}}{\sin(\phi_1)
- \tan(\phi_2)\cos(\phi_1)}\,.
\end{equation}
We get $p_{2\perp}$ by exchanging $1 \leftrightarrow 2$ in the
indices. To obtain the Jacobian of the transformation, we start from the
single jet phase space eq.~(\ref{eq:ps_jetparton}) with the rapidity
delta function already rewritten to be linear in the rapidity of the
last parton, i.e.
\begin{equation}
\label{eq:jet_2p}
\prod_{i=1,2} \bigg(\int
\mathrm{d}\mathbf{p}_{i\perp}\ \int \mathrm{d} y_i
\bigg)\delta^{(2)}\Big(\mathbf{p}_{1\perp} + \mathbf{p}_{2\perp} -
\mathbf{j}_{\perp}^B\Big)\delta(y_2- \dots)\,.
\end{equation}
The integral over the second parton momentum is now trivial; we can just replace
the integral over $y_2$ with the equivalent constraint
\begin{equation}
\label{eq:R2}
\int \mathrm{d}R_2 \ \delta\bigg(R_2 - \bigg[\phi_{\cal J} - \arctan
\bigg(\frac{p_{{\cal J}y} - p_{1y}}{p_{{\cal J}x} -
p_{1x}}\bigg)\bigg]/\cos \Theta\bigg) \,.
\end{equation}
In order to fix the integral over $p_{1\perp}$ instead, we rewrite this
$\delta$ function. This introduces the Jacobian
\begin{equation}
\label{eq:jac_pt1}
\bigg|\frac{\partial p_{1\perp}}{\partial R_2} \bigg| =
\frac{\cos(\Theta)\mathbf{p}_{2\perp}^2}{p_{{\cal J}\perp}\sin(\phi_{\cal J}-\phi_1)}\,.
\end{equation}
The final form of the integral over the two parton momenta is then
\begin{equation}
\label{eq:ps_jet_2p}
\int \mathrm{d}R_1\ R_1 \int \mathrm{d}R_2 \int \mathrm{d}x_\Theta\ 2\pi \int
\mathrm{d}p_{1\perp}\ p_{1\perp} \int \mathrm{d}p_{2\perp}
\ \bigg|\frac{\partial p_{1\perp}}{\partial R_2} \bigg|\delta(p_{1\perp}
-\dots) \delta(p_{2\perp} - \dots)\,.
\end{equation}
As is evident from section~\ref{sec:psp_ng_jet}, jets with three or more
constituents are rare and an efficient phase-space sampling is less
important. For such jets, we exploit the observation that partons with a
distance larger than $R_{\text{max}} = \tfrac{5}{3} R$ to
the jet centre are never clustered into the jet. Assuming $N$
constituents, we generate all components
for the first $N-1$ partons and fix the remaining parton with the
$\delta$-functional. In order to end up inside the jet, we use the
parametrisation
\begin{align}
\label{eq:ps_jet_param}
\phi_i ={}& \phi_{\cal J} + \Delta \phi_i\,, & \Delta \phi_i ={}& \Delta
R_i
\cos(\Theta_i)\,, \\
y_i ={}& y_{\cal J} + \Delta y_i\,, & \Delta y_i ={}& \Delta
R_i
\sin(\Theta_i)\,,
\end{align}
and generate $\Theta_i$ and $\Delta R_i$ randomly with $\Delta R_i \leq
R_{\text{max}}$ and the empiric value $R_{\text{max}} = 5\*R/3$. We can
then write the phase space integral for a single parton as $(p_\perp = |\mathbf{p}_\perp|)$
\begin{equation}
\label{eq:ps_jetparton_x}
\int \mathrm{d}\mathbf{p}_{\perp}\ \int
\mathrm{d} y \approx \int_{\Box} \mathrm{d}x_{\perp}
\mathrm{d}x_{ R}
\mathrm{d}x_{\theta}\
2\*\pi\,\*R_{\text{max}}^2\,\*x_{R}\,\*p_{\perp}\,\*(p_{\perp,\text{max}}
- p_{\perp,\text{min}})
\end{equation}
with
\begin{align}
\label{eq:ps_jetparton_parameters}
\Delta \phi ={}& R_{\text{max}}\*x_{R}\*\cos(2\*\pi\*x_\theta)\,,&
\Delta y ={}& R_{\text{max}}\*x_{R}\*\sin(2\*\pi\*x_\theta)\,, \\
p_{\perp} ={}& (p_{\perp,\text{max}} - p_{\perp,\text{min}})\*x_\perp +
p_{\perp,\text{min}}\,.
\end{align}
$p_{\perp,\text{max}}$ is determined from the requirement that the total
contribution from the first $n-1$ partons --- i.e. the projection onto the
jet $p_{\perp}$ axis --- must never exceed the jet $p_\perp$. This gives
\todo{This bound is too high}
\begin{equation}
\label{eq:pt_max}
p_{i\perp,\text{max}} = \frac{p_{{\cal J}\perp} - \sum_{j<i} p_{j\perp}
\cos \Delta
\phi_j}{\cos \Delta
\phi_i}\,.
\end{equation}
The $x$ and $y$ components of the last parton follow immediately from
the first $\delta$ function. The last rapidity is fixed by the condition that
the jet rapidity is kept fixed by the reshuffling, i.e.
\begin{equation}
\label{eq:yJ_delta}
y^B_{\cal J} = y_{\cal J} = \frac 1 2 \ln \frac{\sum_{i=1}^n E_i+ p_{iz}}{\sum_{i=1}^n E_i - p_{iz}}\,.
\end{equation}
With $E_n \pm p_{nz} = p_{n\perp}\exp(\pm y_n)$ this can be rewritten to
\begin{equation}
\label{eq:yn_quad_eq}
\exp(2y_{\cal J}) = \frac{\sum_{i=1}^{n-1} E_i+ p_{iz}+p_{n\perp} \exp(y_n)}{\sum_{i=1}^{n-1} E_i - p_{iz}+p_{n\perp} \exp(-y_n)}\,,
\end{equation}
which is a quadratic equation in $\exp(y_n)$. The physical solution is
\begin{align}
\label{eq:yn}
y_n ={}& \log\Big(-b + \sqrt{b^2 + \exp(2y_{\cal J})}\,\Big)\,,\\
b ={}& \bigg(\sum_{i=1}^{n-1} E_i + p_{iz} - \exp(2y_{\cal J})
\sum_{i=1}^{n-1} E_i - p_{iz}\bigg)/(2 p_{n\perp})\,.
\end{align}
\todo{what's wrong with the following?} To eliminate the remaining rapidity
integral, we transform the $\delta$ function to be linear in the
rapidity $y$ of the last parton. The corresponding Jacobian is
\begin{equation}
\label{eq:jacobian_y}
\bigg|\frac{\partial y_{\cal J}}{\partial y_n}\bigg|^{-1} = 2 \bigg( \frac{E_n +
p_{nz}}{E_{\cal J} + p_{{\cal J}z}} + \frac{E_n - p_{nz}}{E_{\cal J} -
p_{{\cal J}z}}\bigg)^{-1}\,.
\end{equation}
Finally, we check that all designated constituents are actually
clustered into the considered jet.
\subsubsection{Final steps}
\label{sec:final}
Knowing the rapidity span covered by the extremal partons, we can now
generate the rapdities for the partons outside jets. We perform jet
clustering on all partons and check in
\lstinline!PhaseSpacePoint::jets_ok! that all the following criteria are
fulfilled:
\begin{itemize}
\item The number of resummation jets must match the number of
fixed-order jets.
\item No partons designated to be outside jets may end up inside jets.
\item All other outgoing partons \emph{must} end up inside jets.
\item The extremal (in rapidity) partons must be inside the extremal
jets. If there is, for example, an unordered forward emission, the
most forward parton must end up inside the most forward jet and the
next parton must end up inside second jet.
\item The rapidities of fixed-order and resummation jets must match.
\end{itemize}
After this, we adjust the phase-space normalisation according to the
third line of eq.~(\ref{eq:resumdijetFKLmatched2}), determine the
flavours of the outgoing partons, and adopt any additional colourless
bosons from the fixed-order input event. Finally, we use momentum
conservation to reconstruct the momenta of the incoming partons.
\subsection{Colour connection}
\label{sec:Colour}
\begin{figure}
\input{src/ColourConnect.tex}
\caption{Left: Non-crossing colour flow dominating in the MRK limit. The
crossing of the colour line connecting to particle 2 can be resolved by
writing particle 2 on the left. Right: A colour flow with a (manifest)
colour-crossing. The crossing can only be resolved if one breaks the
rapidities order, e.g. switching particles 2 and 3. From~\cite{Andersen:2017sht}.}
\label{fig:Colour_crossing}
\end{figure}
After the phase space for the resummation event is generated, we can construct
the colour for each particle. To generate the colour flow one has to call
\lstinline!Event::generate_colours! on any \HEJ configuration. For non-\HEJ
event we do not change the colour, and assume it is provided by the user (e.g.
through the LHE file input). The colour connection is done in the large $N_c$
(infinite number of colour) limit with leading colour in
MRK~\cite{Andersen:2008ue, Andersen:2017sht}. The idea is to allow only
$t$-channel colour exchange, without any crossing colour lines. For example the
colour crossing in the colour connection on the left of
figure~\ref{fig:Colour_crossing} can be resolved by switching \textit{particle
2} to the left.
We can write down the colour connections by following the colour flow from
\textit{gluon a} to \textit{gluon b} and back to \textit{gluon a}, e.g.
figure~\ref{fig:Colour_gleft} corresponds to $a123ba$. In such an expression any
valid, non-crossing colour flow will connect all external legs while respecting
the rapidity ordering. Thus configurations like the left of
figure~\ref{fig:Colour_crossing} are allowed ($a134b2a$), but the right of the
same figures breaks the rapidity ordering between 2 and 3 ($a1324ba$). Note that
connections between $b$ and $a$ are in inverse order, e.g. $ab321a$ corresponds to~\ref{fig:Colour_gright} ($a123ba$) just with colour and anti-colour swapped.
\begin{figure}
\centering
\subcaptionbox{$a123ba$\label{fig:Colour_gright}}{
\includegraphics[height=0.25\textwidth]{colour_gright.jpg}}
\subcaptionbox{$a13b2a$\label{fig:Colour_gleft}}{
\includegraphics[height=0.25\textwidth]{colour_gleft.jpg}}
\subcaptionbox{$a\_123ba$\label{fig:Colour_qx}}{
\includegraphics[height=0.25\textwidth]{colour_qx.jpg}}
\subcaptionbox{$a\_23b1a$\label{fig:Colour_uno}}{
\includegraphics[height=0.25\textwidth]{colour_uno.jpg}}
\subcaptionbox{$a14b3\_2a$\label{fig:Colour_qqx}}{
\includegraphics[height=0.25\textwidth]{colour_centralqqx.jpg}}
\caption{Different colour non-crossing colour connections. Both incoming
particles are drawn at the top or bottom and the outgoing left or right.
The Feynman diagram is shown in black and the colour flow in blue.}
%TODO Maybe make these plots nicer (in Latex/asy)
\end{figure}
If we replace two gluons with a quark, (anti-)quark pair we break one of the
colour connections. Still the basic concept from before holds, we just have to
treat the connection between two (anti-)quarks like an unmovable (anti-)colour.
We denote such a connection by a underscore (e.g. $1\_a$). For example the
equivalent of~\ref{fig:Colour_gright} ($a123ba$) with an incoming antiquark
is~\ref{fig:Colour_qx} ($a\_123ba$). As said this also holds for other
subleading configurations like unordered emission~\ref{fig:Colour_uno} or
central quark-antiquark pairs~\ref{fig:Colour_qqx} \footnote{Obviously this can
not be guaranteed for non-\HEJ configurations, e.g. $qQ\to Qq$ requires a
$u$-channel exchange.}.
Some rapidity ordering can have multiple possible colour connections,
e.g.~\ref{fig:Colour_gright} and~\ref{fig:Colour_gleft}. This is always the case
if a gluon radiates off a gluon line. In that case we randomly connect the gluon
to either the colour or anti-colour. Thus in the generation we keep track
whether we are on a quark or gluon line, and act accordingly.
\subsection{The matrix element }
\label{sec:ME}
The derivation of the \HEJ matrix element is explained in some detail
in~\cite{Andersen:2017kfc}, where also results for leading and
subleading matrix elements for pure multijet production and production
of a Higgs boson with at least two associated jets are listed. Matrix
elements for $Z/\gamma^*$ production together with jets are
given in~\cite{Andersen:2016vkp}, but not yet included.
A full list of all implemented currents is given in
section~\ref{sec:currents_impl}.
The matrix elements are implemented in the \lstinline!MatrixElement!
class. To discuss the structure, let us consider the squared matrix
element for FKL multijet production with $n$ final-state partons:
\begin{align}
\label{eq:ME}
\begin{split}
\overline{\left|\mathcal{M}_\HEJ^{f_1 f_2 \to f_1
g\cdots g f_2}\right|}^2 = \ &\frac {(4\pi\alpha_s)^n} {4\ (N_c^2-1)}
\cdot\ \textcolor{blue}{\frac {K_{f_1}(p_1^-, p_a^-)} {t_1}\ \cdot\ \frac{K_{f_2}(p_n^+, p_b^+)}{t_{n-1}}\ \cdot\ \left\|S_{f_1 f_2\to f_1 f_2}\right\|^2}\\
& \cdot \prod_{i=1}^{n-2} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
& \cdot \prod_{j=1}^{n-1} \textcolor{red}{\exp\left[\omega^0(q_{j\perp})(y_{j+1}-y_j)\right]}.
\end{split}
\end{align}
The structure and momentum assignment of the unsquared matrix element is
as illustrated here:
\begin{center}
\includegraphics{HEJ_amplitude}
\end{center}
The square
of the complete matrix element as given in eq.~\eqref{eq:ME} is
calculated by \lstinline!MatrixElement::operator()!. The \textcolor{red}{last line} of
eq.~\eqref{eq:ME} constitutes the all-order virtual correction,
implemented in
\lstinline!MatrixElement::virtual_corrections!.
$\omega^0$ is the
\textit{regularised Regge trajectory}
\begin{equation}
\label{eq:omega_0}
\omega^0(q_\perp) = - C_A \frac{\alpha_s}{\pi} \log \left(\frac{q_\perp^2}{\lambda^2}\right)\,,
\end{equation}
where $\lambda$ is the slicing parameter limiting the softness of real
gluon emissions, cf. eq.~\eqref{eq:resumdijetFKLmatched2}. $\lambda$ can be
changed at runtime by setting \lstinline!regulator parameter! in
\lstinline!conifg.yml!.
The remaining parts, which correspond to the square of the leading-order
\HEJ matrix element $\overline{\left|\mathcal{M}_\text{LO,
\HEJ}^{f_1f_2\to f_1g\cdots
gf_2}\big(\big\{p^B_j\big\}\big)\right|}^{2}$, are computed in
\lstinline!MatrixElement::tree!. We can further factor off the
scale-dependent ``parametric'' part
\lstinline!MatrixElement::tree_param! containing all factors of the
strong coupling $4\pi\alpha_s$. Using this function saves some CPU time
when adjusting the renormalisation scale, see
section~\ref{sec:resum}. The remaining ``kinematic'' factors are
calculated in \lstinline!MatrixElement::kin!.
\subsubsection{Matrix elements for Higgs plus dijet}
\label{sec:ME_h_jets}
In the production of a Higgs boson together with jets the parametric
parts and the virtual corrections only require minor changes in the
respective functions. However, in the ``kinematic'' parts we have to
distinguish between several cases, which is done in
\lstinline!MatrixElement::tree_kin_Higgs!. The Higgs boson can be
\emph{central}, i.e. inside the rapidity range spanned by the extremal
partons (\lstinline!MatrixElement::tree_kin_Higgs_central!) or
\emph{peripheral} and outside this range
(\lstinline!MatrixElement::tree_kin_Higgs_first! or
\lstinline!MatrixElement::tree_kin_Higgs_last!). Currently the current for an
unordered emission with an Higgs on the same side it not implemented
\footnote{In principle emitting a Higgs boson \textit{on the other
side} of the unordered gluon is possible by contracting an unordered and
external Higgs current. Obviously this would not cover all possible
configurations, e.g. $qQ\to HgqQ$ requires contraction of the standard $Q\to Q$
current with an (unknown) $q\to Hgq$ one.}.
If a Higgs boson with momentum $p_H$ is emitted centrally, after parton
$j$ in rapidity, the matrix element reads
\begin{equation}
\label{eq:ME_h_jets_central}
\begin{split}
\overline{\left|\mathcal{M}_\HEJ^{f_1 f_2 \to f_1 g\cdot H
\cdot g f_2}\right|}^2 = \ &\frac {\alpha_s^2 (4\pi\alpha_s)^n} {4\ (N_c^2-1)}
\cdot\ \textcolor{blue}{\frac {K_{f_1}(p_1^-, p_a^-)} {t_1}\
\cdot\ \frac{1}{t_j t_{j+1}} \cdot\ \frac{K_{f_2}(p_n^+, p_b^+)}{t_{n}}\ \cdot\ \left\|S_{f_1
f_2\to f_1 H f_2}\right\|^2}\\
& \cdot \prod_{\substack{i=1\\i \neq j}}^{n-1} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
& \cdot \textcolor{red}{\prod_{i=1}^{n-1}
\exp\left[\omega^0(q_{i\perp})\Delta y_i\right]}
\end{split}
\end{equation}
with the momentum definitions
\begin{center}
\includegraphics{HEJ_central_Higgs_amplitude}
\end{center}
$q_i$ is the $i$th $t$-channel momentum and $\Delta y_i$ the rapidity
gap between outgoing \emph{particles} (not partons) $i$ and $i+1$ in
rapidity ordering.
For \emph{peripheral} emission in the backward direction
(\lstinline!MatrixElement::tree_kin_Higgs_first!) we first check whether
the most backward parton is a gluon or an (anti-)quark. In the latter
case the leading contribution to the matrix element arises through
emission off the $t$-channel gluons and we can use the same formula
eq.~(\ref{eq:ME_h_jets_central}) as for central emission. If the most
backward parton is a gluon, the square of the matrix element can be
written as
\begin{equation}
\label{eq:ME_h_jets_peripheral}
\begin{split}
\overline{\left|\mathcal{M}_\HEJ^{g f_2 \to H g\cdot g f_2}\right|}^2 = \ &\frac {\alpha_s^2 (4\pi\alpha_s)^n} {\textcolor{blue}{4\ (N_c^2-1)}}
\textcolor{blue}{\cdot\ K_{H}\
\frac{K_{f_2}(p_n^+, p_b^+)}{t_{n-1}}\ \cdot\ \left\|S_{g
f_2\to H g f_2}\right\|^2}\\
& \cdot \prod_{\substack{i=1}}^{n-2} \textcolor{gray}{\left( \frac{-C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)}\\
& \cdot \textcolor{red}{\prod_{i=1}^{n-1}
\exp\left[\omega^0(q_{i\perp}) (y_{i+1} - y_i)\right]}
\end{split}
\end{equation}
with the momenta as follows:
\begin{center}
\includegraphics{HEJ_peripheral_Higgs_amplitude}
\end{center}
The \textcolor{blue}{blue part} is implemented in
\lstinline!MatrixElement::MH2_forwardH!. All other building blocks are
already available.\todo{Impact factors} The actual current contraction
is calculated in \lstinline!MH2gq_outsideH! inside
\lstinline!src/Hjets.cc!, which corresponds to $\tfrac{16 \pi^2}{t_1} \left\|S_{g
f_2\to H g f_2}\right\|^2$.\todo{Fix this insane normalisation}
The forward emission of a Higgs boson is completely analogous. We can
use the same function \lstinline!MatrixElement::MH2_forwardH!, swapping
$p_1 \leftrightarrow p_n,\,p_a \leftrightarrow p_b$.
\subsubsection{FKL ladder and Lipatov vertices}
\label{sec:FKL_ladder}
The ``FKL ladder'' is the product
\begin{equation}
\label{eq:FKL_ladder}
\prod_{i=1}^{n-2} \left( \frac{-C_A}{t_it_{i+1}}\
V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1}) \right)
\end{equation}
appearing in the square of the matrix element for $n$ parton production,
cf. eq.~(\ref{eq:ME}), and implemented in
\lstinline!MatrixElement::FKL_ladder_weight!. The Lipatov vertex contraction
$V^\mu(q_i,q_{i+1})V_\mu(q_i,q_{i+1})$ is implemented \lstinline!C2Lipatovots!.
It is given by \todo{equation} \todo{mention difference between the two versions
of \lstinline!C2Lipatovots!, maybe even get rid of one}.
\subsubsection{Currents}
\label{sec:currents}
The current factors $\frac{K_{f_1}K_{f_2}}{t_1 t_{n-1}}\left\|S_{f_1
f_2\to f_1 f_2}\right\|^2$ and their extensions for unordered and Higgs
boson emissions are implemented in the \lstinline!jM2!$\dots$ functions
of \texttt{src/Hjets.cc}. \todo{Only $\|S\|^2$ should be in currents}
\footnote{The current implementation for
Higgs production in \texttt{src/Hjets.cc} includes the $1/4$ factor
inside $S$, opposing to~\eqref{eq:ME}. Thus the overall normalisation is
unaffected.} The ``colour acceleration multiplier'' (CAM) $K_{f}$
for a parton $f\in\{g,q,\bar{q}\}$ is defined as
\begin{align}
\label{eq:K_g}
K_g(p_1^-, p_a^-) ={}& \frac{1}{2}\left(\frac{p_1^-}{p_a^-} + \frac{p_a^-}{p_1^-}\right)\left(C_A -
\frac{1}{C_A}\right)+\frac{1}{C_A}\\
\label{eq:K_q}
K_q(p_1^-, p_a^-) ={}&K_{\bar{q}}(p_1^-, p_a^-) = C_F\,.
\end{align}
The Higgs current CAM used in eq.~(\ref{eq:ME_h_jets_peripheral}) is
\begin{equation}
\label{eq:K_H}
K_H = C_A\,.
\end{equation}
The current contractions are given by\todo{check all this
carefully!}
\begin{align}
\label{eq:S}
\left\|S_{f_1 f_2\to f_1 f_2}\right\|^2 ={}& \sum_{\substack{\lambda_a =
+,-\\\lambda_b = +,-}} \left|j^{\lambda_a}_\mu(p_1, p_a)\
j^{\lambda_b\,\mu}(p_n, p_b)\right|^2 = 2\sum_{\lambda =
+,-} \left|j^{-}_\mu(p_1, p_a)\ j^{\lambda\,\mu}(p_n, p_b)\right|^2\,,\\
\left\|S_{f_1 f_2\to f_1 H f_2}\right\|^2 ={}& \sum_{\substack{\lambda_a =
+,-\\\lambda_b = +,-}} \left|j^{\lambda_a}_\mu(p_1, p_a)V_H^{\mu\nu}(q_j, q_{j+1})\
j^{\lambda_b}_\nu(p_n, p_b)\right|^2\,,\\
\left\|S_{g f_2 \to H g f_2}\right\|^2 ={}& \sum_{
\substack{
\lambda_{a} = +,-\\
\lambda_{1} =+,-\\
\lambda_{b} = +,-
}}
\left|j^{\lambda_a\lambda_1}_{H\,\mu}(p_1, p_a, p_H)\ j^{\lambda_b\,\mu}(p_n, p_b)\right|^2\,.
\end{align}
Currently, most of the expressions for the current contractions are hard-coded. Work is ongoing to use the current generator described in section~\ref{sec:cur_gen} instead.
\subsection{Unweighting}
\label{sec:unweight}
Straightforward event generation tends to produce many events with small
weights. Those events have a negligible contribution to the final
observables, but can take up considerable storage space and CPU time in
later processing stages. This problem can be addressed by unweighting.
For naive unweighting, one would determine the maximum weight
$w_\text{max}$ of all events, discard each event with weight $w$ with a
probability $p=w/w_\text{max}$, and set the weights of all remaining
events to $w_\text{max}$. The downside to this procedure is that it also
eliminates a sizeable fraction of events with moderate weight, so that
the statistical convergence deteriorates. Naive unweighting can be
performed by using the \lstinline!set_cut_to_maxwt! member function of the
\lstinline!Unweighter! on the events and then call the
\lstinline!unweight! member function. It can be enabled for the
resummation events as explained in the user documentation.
To ameliorate the problem of naive unweighting, we also implement
partial unweighting. That is, we perform unweighting only for events
with sufficiently small weights. When using the \lstinline!Unweighter!
member function \lstinline!set_cut_to_peakwt! we estimate the mean and
width of the weight-weight distribution from a sample of events. We
use these estimates to determine the maximum weight below which
unweighting is performed; events with a larger weight are not
touched. The actual unweighting is again done in the
\lstinline!Unweighter::unweight! function.
To estimate the peak weight we employ the following heuristic
algorithm. For a calibration sample of $n$ events, create a histogram
with $b=\sqrt{n}$ equal-sized bins. The histogram ranges from $
\log(\min |w_i|)$ to $\log(|\max w_i|)$, where $w_i$ are the event weights. For
each event, add $|w_i|$ to the corresponding bin. We then prune the
histogram by setting all bins containing less than $c=n/b$
events to zero. This effectively removes statistical outliers. The
logarithm of the peak weight is then the centre of the highest bin in
the histogram. In principle, the number of bins $b$ and the pruning
parameter $c$ could be tuned further.
To illustrate the principle, here is a weight-weight histogram
filled with a sample of 100000 event weights before the pruning:
\begin{center}
\includegraphics[width=0.7\linewidth]{wtwt}
\end{center}
The peaks to the right are clearly outliers caused by single
events. After pruning we get the following histogram:
\begin{center}
\includegraphics[width=0.7\linewidth]{wtwt_cut}
\end{center}
The actual peak weight probably lies above the cut, and the algorithm
can certainly be improved. Still, the estimate we get from the pruned
histogram is already good enough to eliminate about $99\%$ of the
low-weight events.
\section{The current generator}
\label{sec:cur_gen}
The current generator in the \texttt{current\_generator} directory is
automatically invoked when building \HEJ 2. Its task is to compute and
generate C++ code for the current contractions listed in
section~\ref{sec:currents}. For each source file \texttt{<j>.frm}
inside \texttt{current\_generator} it generates a corresponding header
\texttt{include/HEJ/currents/<j>.hh} inside the build directory. The
header can be included with
\begin{lstlisting}[language=C++,caption={}]
#include "HEJ/currents/<j>.hh"
\end{lstlisting}
\subsection{Implementing new current contractions}
\label{sec:cur_gen_new}
To implement a new current contraction \lstinline!jcontr! create a new
file \texttt{jcontr.frm} inside the \texttt{current\_generator}
directory. \texttt{jcontr.frm} should contain FORM code, see
\url{https://www.nikhef.nl/~form/} for more information on FORM. Here
is a small example:
\begin{lstlisting}[caption={}]
* FORM comments are started by an asterisk * at the beginning of a line
* First include the relevant headers
#include- include/helspin.frm
#include- include/write.frm
* Define the symbols that appear.
* UPPERCASE symbols are reserved for internal use
vectors p1,...,p10;
indices mu1,...,mu10;
* Define local expressions of the form [NAME HELICITIES]
* for the current contractions for all relevant helicity configurations
#do HELICITY1={+,-}
#do HELICITY2={+,-}
* We use the Current function
* Current(h, p1, mu1, ..., muX, p2) =
* u(p1) \gamma_{mu1} ... \gamma_{muX} u(p2)
* where h=+1 or h=-1 is the spinor helicity.
* All momenta appearing as arguments have to be *lightlike*
local [jcontr `HELICITY1'`HELICITY2'] =
Current(`HELICITY1'1, p1, mu1, p2, mu2, p3)
*Current(`HELICITY2'1, p4, mu2, p2, mu1, p1);
#enddo
#enddo
.sort
* Main procedure that calculates the contraction
#call ContractCurrents
.sort
* Optimise expression
format O4;
* Format in a (mostly) c compatible way
format c;
* Write start of C++ header file
#call WriteHeader(`OUTPUT')
* Write a template function jcontr
* taking as template arguments two helicities
* and as arguments the momenta p1,...,p4
* returning the contractions [jcontr HELICITIES] defined above
#call WriteOptimised(`OUTPUT',jcontr,2,p1,p2,p3,p4)
* Wrap up
#call WriteFooter(`OUTPUT')
.end
\end{lstlisting}
\subsection{Calculation of contractions}
\label{sec:contr_calc}
In order to describe the algorithm for the calculation of current
contractions we first have to define the currents and establish some
useful relations.
+\subsubsection{Massive momenta}
+\label{sec:p_massive}
+
+We want to use relations for lightlike momenta. Momenta $P$ that are
+\emph{not} lightlike can be written as the sum of two lightlike
+momenta:
+\begin{equation}
+ \label{eq:P_massive}
+ P^\mu = p^\mu + q^\mu\,, \qquad p^2 = q^2 = 0 \,.
+\end{equation}
+This decomposition is not unique. If we impose the arbitrary
+constraint $q_\perp = 0$ and require real-valued momentum
+components we can use the ansatz
+\begin{align}
+ \label{eq:P_massive_p}
+ p_\mu ={}& P_\perp\*(\cosh y, \cos \phi, \sin \phi, \sinh y)\,,\\
+ \label{eq:P_massive_q}
+ q_\mu ={}& (E, 0, 0, s\,E)\,,
+\end{align}
+where $P_\perp$ is the transverse momentum of $P$ and $\phi$ the corresponding azimuthal angle. For the remaining parameters we obtain
+\begin{align}
+ \label{eq:P_massive_plus}
+P^+ > 0:& & y ={}& \log \frac{P^+}{P_\perp}\,,\qquad E = \frac{P^2}{2P^+}\,,\qquad s = -1\,,\\
+ \label{eq:P_massive_minus}
+P^- > 0:& & y ={}& \log \frac{P_\perp}{P^-}\,,\qquad E = \frac{P^2}{2P^-}\,,\qquad s = +1\,.
+\end{align}
+
+\subsubsection{Currents and current relations}
+\label{sec:current_relations}
+
Our starting point are generalised currents
\begin{equation}
\label{eq:j_gen}
j^{\pm}(p, \mu_1,\dots,\mu_{2N-1},q) = \bar{u}^{\pm}(p)\gamma_{\mu_1} \dots \gamma_{\mu_{2N-1}} u^\pm(q)\,.
\end{equation}
Since there are no masses, we can consider two-component chiral spinors
\begin{align}
\label{eq:u_plus}
u^+(p)={}& \left(\sqrt{p^+}, \sqrt{p^-} \hat{p}_\perp \right) \,,\\
\label{eq:u_minus}
u^-(p)={}& \left(\sqrt{p^-} \hat{p}^*_\perp, -\sqrt{p^+}\right)\,,
\end{align}
with $p^\pm = E\pm p_z,\, \hat{p}_\perp = \tfrac{p_\perp}{|p_\perp|},\,
p_\perp = p_x + i p_y$. The spinors for vanishing transverse momentum
are obtained by replacing $\hat{p}_\perp \to -1$.
This gives
\begin{equation}
\label{eq:jminus_gen}
j^-(p,\mu_1,\dots,\mu_{2N-1},q) = u^{-,\dagger}(p)\ \sigma^-_{\mu_1}\ \sigma^+_{\mu_2}\dots\sigma^-_{\mu_{2N-1}}\ u^{-}(q)\,.
\end{equation}
where $\sigma_\mu^\pm = (1, \pm \sigma_i)$ and $\sigma_i$ are the Pauli
matrices
\begin{equation}
\label{eq:Pauli_matrices}
\sigma_1 =
\begin{pmatrix}
0 & 1\\ 1 & 0
\end{pmatrix}
\,,
\qquad \sigma_2 =
\begin{pmatrix}
0 & -i\\ i & 0
\end{pmatrix}
\,,
\qquad \sigma_3 =
\begin{pmatrix}
1 & 0\\ 0 & -1
\end{pmatrix}
\,.
\end{equation}
For positive-helicity currents we can either flip all helicities in
eq.~(\ref{eq:jminus_gen}) or reverse the order of the arguments, i.e.
\begin{equation}
\label{eq:jplus_gen}
j^+(p,\mu_1,\dots,\mu_{2N-1},q) = \big(j^-(p,\mu_1,\dots,\mu_{2N-1},q)\big)^* = j^-(q,\mu_{2N-1},\dots,\mu_1,p) \,.
\end{equation}
Using the standard spinor-helicity notation we have
\begin{gather}
\label{eq:spinors_spinhel}
u^+(p) = | p \rangle\,, \qquad u^-(p) = | p ]\,, \qquad u^{+,\dagger}(p) = [ p |\,, \qquad u^{-,\dagger}(p) = \langle p |\,,\\
\label{eq:current_spinhel}
j^-(p,\mu_1,\dots,\mu_{2N-1},q) = \langle p |\ \mu_1\ \dots\ \mu_{2N-1}\ | q ] \,.\\
\label{eq:contraction_spinhel}
P_{\mu_i} j^-(p,\mu_1,\dots,\mu_{2N-1},q) = \langle p |\ \mu_1\ \dots\ \mu_{i-1}\ P\ \mu_{i+1}\ \dots\ \mu_{2N-1}\ | q ] \,.
\end{gather}
Lightlike momenta $p$ can be decomposed into spinor products:
\begin{equation}
\label{eq:p_decomp}
\slashed{p} = |p\rangle [p| + |p] \langle p |\,.
\end{equation}
Taking into account helicity conservation this gives the following relations:
\begingroup
\addtolength{\jot}{1em}
\begin{align}
\label{eq:p_in_current}
\langle p |\ \mu_1\ \dots\ \mu_i\ P\ \mu_{i+1}\ \dots\ \mu_{2N-1}\ | q ] ={}&
\begin{cases}
\langle p |\ \mu_1\ \dots\ \mu_i\ |P]\ \langle P|\ \mu_{i+1}\ \dots\ \mu_{2N-1}\ | q ]& i \text{ even}\\
\langle p |\ \mu_1\ \dots\ \mu_i\ |P\rangle\ [ P|\ \mu_{i+1}\ \dots\ \mu_{2N-1}\ | q ]& i \text{ odd}
\end{cases}\,,\\
\label{eq:p_in_angle}
\langle p |\ \mu_1\ \dots\ \mu_i\ P\ \mu_{i+1}\ \dots\ \mu_{2N}\ | q \rangle ={}&
\begin{cases}
\langle p |\ \mu_1\ \dots\ \mu_i\ |P]\ \langle P| \mu_{i+1}\ \dots\ \mu_{2N}\ | q \rangle & i \text{ even}\\
\langle p |\ \mu_1\ \dots\ \mu_i\ |P\rangle\ \big(\langle P| \mu_{i+1}\ \dots\ \mu_{2N}\ | q ]\big)^* & i \text{ odd}
\end{cases}\,,\\
\label{eq:p_in_square}
[ p |\ \mu_1\ \dots\ \mu_i\ P\ \mu_{i+1}\ \dots\ \mu_{2N}\ | q ] ={}&
\begin{cases}
\big(\langle p |\ \mu_1\ \dots\ \mu_i\ |P]\big)^* \ [ P| \mu_{i+1}\ \dots\ \mu_{2N}\ | q ] & i \text{ even}\\
[ p |\ \mu_1\ \dots\ \mu_i\ |P]\ \langle P| \mu_{i+1}\ \dots\ \mu_{2N}\ | q ] & i \text{ odd}
\end{cases}\,.
\end{align}
\endgroup
For contractions of vector currents we can use the Fierz identity
\begin{equation}
\label{eq:Fierz}
\langle p|\ \mu\ |q]\ \langle k|\ \mu\ |l] = 2 \spa p.k \spb l.q\,.
\end{equation}
The scalar angle and square products are given by
\begin{align}
\label{eq:angle_product}
\spa p.q ={}& {\big(u^-(p)\big)}^\dagger u^+(q) =
\sqrt{p^-q^+}\hat{p}_{i,\perp} - \sqrt{p^+q^-}\hat{q}_{j,\perp} = - \spa q.p\,,\\
\label{eq:square_product}
\spb p.q ={}& {\big(u^+(p)\big)}^\dagger u^-(q) = -\spa p.q ^* = - \spb q.p\,.
\end{align}
+\subsubsection{Contraction algorithm}
+\label{sec:contr_calc_algo}
+
The contractions are now calculated as follows:
\begin{enumerate}
\item Use equations \eqref{eq:jplus_gen}, \eqref{eq:current_spinhel} to write all currents in a canonical form.
\item Assume that all momenta are lightlike and use the relations
\eqref{eq:p_in_current}, \eqref{eq:p_in_angle}, \eqref{eq:p_in_square}
to split up currents that are contracted with momenta.
\item Apply the Fierz transformation~\eqref{eq:Fierz} to eliminate
contractions between vector currents.
\item Write the arguments of the antisymmetric angle and scalar products in canonical order, see equations~\eqref{eq:angle_product} ,\eqref{eq:square_product}.
\end{enumerate}
\section{The fixed-order generator}
\label{sec:HEJFOG}
Even at leading order, standard fixed-order generators can only generate
events with a limited number of final-state particles within reasonable
CPU time. The purpose of the fixed-order generator is to supplement this
with high-multiplicity input events according to the first two lines of
eq.~\eqref{eq:resumdijetFKLmatched2} with the \HEJ approximation
$\mathcal{M}_\text{LO, \HEJ}^{f_1f_2\to f_1g\cdots gf_2}$ instead of the
full fixed-order matrix element $\mathcal{M}_\text{LO}^{f_1f_2\to
f_1g\cdots gf_2}$. Its usage is described in the user
documentation \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJFOG.html}.
\subsection{File structure}
\label{sec:HEJFOG_structure}
The code for the fixed-order generator is in the \texttt{FixedOrderGen}
directory, which contains the following:
\begin{description}
\item[include:] Contains the C++ header files.
\item[src:] Contains the C++ source files.
\item[t:] Contains the source code for the automated tests.
\item[CMakeLists.txt:] Configuration file for the \cmake build system.
\item[configFO.yml:] Sample configuration file for the fixed-order generator.
\end{description}
The code is generally in the \lstinline!HEJFOG! namespace. Functions and
classes \lstinline!MyClass! are usually declared in
\texttt{include/MyClass.hh} and implemented in \texttt{src/MyClass.cc}.
\subsection{Program flow}
\label{sec:prog_flow}
A single run of the fixed-order generator consists of three or four
stages.
First, we perform initialisation similar to \HEJ 2, see
section~\ref{sec:init}. Since there is a lot of overlap we frequently
reuse classes and functions from \HEJ 2, i.e. from the
\lstinline!HEJ! namespace. The code for parsing the configuration file
is in \texttt{include/config.hh} and implemented in
\texttt{src/config.cc}.
If partial unweighting is requested in the user settings \url{https://hej.web.cern.ch/HEJ/doc/current/user/HEJFOG.html#settings},
the initialisation is followed by a calibration phase. We use a
\lstinline!EventGenerator! to produce a number of trial
events. We use these to calibrate the \lstinline!Unweighter! in
its constructor and produce a first batch of partially unweighted
events. This also allows us to estimate our unweighting efficiency.
In the next step, we continue to generate events and potentially
unweight them. Once the user-defined target number of events is reached,
we adjust their weights according to the number of required trials. As
in \HEJ 2 (see section~\ref{sec:processing}), we pass the final
events to a number of \lstinline!HEJ::Analysis! objects and a
\lstinline!HEJ::CombinedEventWriter!.
\subsection{Event generation}
\label{sec:evgen}
Event generation is performed by the
\lstinline!EventGenerator::gen_event! member function. We begin by generating a
\lstinline!PhaseSpacePoint!. This is not to be confused with
the resummation phase space points represented by
\lstinline!HEJ::PhaseSpacePoint!! After jet clustering, we compute the
leading-order matrix element (see section~\ref{sec:ME}) and pdf factors.
The phase space point generation is performed in the
\lstinline!PhaseSpacePoint! constructor. We first construct the
user-defined number of $n_p$ partons (by default gluons) in
\lstinline!PhaseSpacePoint::gen_LO_partons!. We use flat sampling in
rapidity and azimuthal angle. For the scalar transverse momenta, we
distinguish between two cases. By default, they are generated based on a
random variable $x_{p_\perp}$ according to
\begin{equation}
\label{eq:pt_sampling}
p_\perp = p_{\perp,\text{min}} +
\begin{cases}
p_{\perp,\text{par}}
\tan\left(
x_{p_\perp}
\arctan\left(
\frac{p_{\perp,\text{max}} - p_{\perp,\text{min}}}{p_{\perp,\text{par}}}
\right)
\right)
& y < y_\text{cut}
\\
- \tilde{p}_{\perp,\text{par}}\log\left(1 - x_{p_\perp}\left[1 -
\exp\left(\frac{p_{\perp,\text{min}} -
p_{\perp,\text{max}}}{\tilde{p}_{\perp,\text{par}}}\right)\right]\right)
& y \geq y_\text{cut}
\end{cases}\,,
\end{equation}
where $p_{\perp,\text{min}}$ is the minimum jet transverse momentum,
$p_{\perp,\text{max}}$ is the maximum transverse parton momentum,
tentatively set to the beam energy, and $y_\text{cut}$, $p_{\perp,\text{par}}$
and $\tilde{p}_{\perp,\text{par}}$ are generation parameters set to
heuristically determined values of
\begin{align}
y_\text{cut}&=3,\\
p_{\perp,\text{par}}&=p_{\perp,\min}+\frac{n_p}{5}, \\
\tilde{p}_{\perp,\text{par}}&=\frac{p_{\perp,\text{par}}}{1 +
5(y-y_\text{cut})}.
\end{align}
The problem with this generation is that the transverse momenta peak at
the minimum transverse momentum required for fixed-order jets. However,
if we use the generated events as input for \HEJ resummation, events
with such soft transverse momenta hardly contribute, see
section~\ref{sec:ptj_res}. To generate efficient input for resummation,
there is the user option \texttt{peak pt}, which specifies the
dominant transverse momentum for resummation jets. If this option is
set, most jets will be generated as above, but with
$p_{\perp,\text{min}}$ set to the peak transverse momentum $p_{\perp,
\text{peak}}$. In addition, there is a small chance of around $2\%$ to
generate softer jets. The heuristic ansatz for the transverse momentum
distribution in the ``soft'' region is
\begin{equation}
\label{FO_pt_soft}
\frac{\partial \sigma}{\partial p_\perp} \propto e^{n_p\frac{p_\perp- p_{\perp,
\text{peak}}}{\bar{p}_\perp}}\,,
\end{equation}
where $n_p$ is the number of partons and $\bar{p}_\perp \approx
4\,$GeV. To achieve this distribution, we use
\begin{equation}
\label{eq:FO_pt_soft_sampling}
p_\perp = p_{\perp, \text{peak}} + \bar{p}_\perp \frac{\log x_{p_\perp}}{n_p}
\end{equation}
and discard the phase space point if the parton is too soft, i.e. below the threshold for
fixed-order jets.
After ensuring that all partons form separate jets, we generate any
potential colourless emissions. We then determine the incoming momenta
and flavours in \lstinline!PhaseSpacePoint::reconstruct_incoming! and
adjust the outgoing flavours to ensure an FKL configuration. Finally, we
may reassign outgoing flavours to generate suppressed (for example
unordered) configurations.
\input{currents}
\input{tensor}
\appendix
\section{Continuous Integration}
\label{sec:CI}
Whenever you are implementing something new or fixed a bug, please also add a
test for the new behaviour to \texttt{t/CMakeLists.txt} via
\lstinline!add_test!. These test can be triggered by running
\lstinline!make test! or \lstinline!ctest! after compiling. A typical test
should be at most a few seconds, so it can be potentially run on each commit
change by each developer. If you require a longer, more careful test, preferably
on top of a small one, surround it with
\begin{lstlisting}[caption={}]
if(${TEST_ALL})
add_test(
NAME t_feature
COMMAND really_long_test
)
endif()
\end{lstlisting}
Afterwards you can execute the longer tests with\footnote{No recompiling is
needed, as long as only the \lstinline!add_test! command is guarded, not the
compiling commands itself.}
\begin{lstlisting}[language=sh,caption={}]
cmake base/directory -DTEST_ALL=TRUE
make test
\end{lstlisting}
On top of that you should add
\href{https://en.cppreference.com/w/cpp/error/assert}{\lstinline!assert!s} in
the code itself. They are only executed when compiled with
\lstinline!CMAKE_BUILD_TYPE=Debug!, without slowing down release code. So you
can use them everywhere to test \textit{expected} or \textit{assumed} behaviour,
e.g. requiring a Higgs boson or relying on rapidity ordering.
GitLab provides ways to directly test code via \textit{Continuous integrations}.
The CI is controlled by \texttt{.gitlab-ci.yml}. For all options for the YAML
file see \href{https://docs.gitlab.com/ee/ci/yaml/}{docs.gitlab.com/ee/ci/yaml/}.https://gitlab.dur.scotgrid.ac.uk/hej/docold/tree/master/Theses
GitLab also provides a small tool to check that YAML syntax is correct under
\lstinline!CI/CD > Pipelines > CI Lint! or
\href{https://gitlab.dur.scotgrid.ac.uk/hej/HEJ/-/ci/lint}{gitlab.dur.scotgrid.ac.uk/hej/HEJ/-/ci/lint}.
Currently the CI is configured to trigger a \textit{Pipeline} on each
\lstinline!git push!. The corresponding \textit{GitLab runners} are configured
under \lstinline!CI/CD Settings>Runners! in the GitLab UI. All runners use a
\href{https://www.docker.com/}{docker} image as virtual environments\footnote{To
use only Docker runners set the \lstinline!docker! tag in
\texttt{.gitlab-ci.yml}.}. The specific docker images maintained separately. If
you add a new dependencies, please also provide a docker image for the CI. The
goal to be able to test \HEJ with all possible configurations.
Each pipeline contains multiple stages (see \lstinline!stages! in
\texttt{.gitlab-ci.yml}) which are executed in order from top to bottom.
Additionally each stage contains multiple jobs. For example the stage
\lstinline!build! contains the jobs \lstinline!build:basic!,
\lstinline!build:qcdloop!, \lstinline!build:rivet!, etc., which compile \HEJ for
different environments and dependencies, by using different Docker images. Jobs
starting with an dot are ignored by the Runner, e.g. \lstinline!.HEJ_build! is
only used as a template, but never executed directly. Only after all jobs of the
previous stage was executed without any error the next stage will start.
To pass information between multiple stages we use \lstinline!artifacts!. The
runner will automatically load all artifacts form all \lstinline!dependencies!
for each job\footnote{If no dependencies are defined \textit{all} artifacts from
all previous jobs are downloaded. Thus please specify an empty dependence if you
do not want to load any artifacts.}. For example the compiled \HEJ code from
\lstinline!build:basic! gets loaded in \lstinline!test:basic! and
\lstinline!FOG:build:basic!, without recompiling \HEJ again. Additionally
artifacts can be downloaded from the GitLab web page, which could be handy for
debugging.
We also trigger some jobs \lstinline!only! on specific events. For example we
only push the code to
\href{https://phab.hepforge.org/source/hej/repository/v2.0/}{HepForge} on
release branches (e.g. v2.0). Also we only execute the \textit{long} tests for
merge requests, on pushes for any release or the \lstinline!master! branch, or
when triggered manually from the GitLab web page.
The actual commands are given in the \lstinline!before_script!,
\lstinline!script! and \lstinline!after_script!
\footnote{\lstinline!after_script! is always executed} sections, and are
standard Linux shell commands (dependent on the docker image). Any failed
command, i.e. returning not zero, stops the job and making the pipeline fail
entirely. Most tests are just running \lstinline!make test! or are based on it.
Thus, to emphasise it again, write tests for your code in \lstinline!cmake!. The
CI is only intended to make automated testing in different environments easier.
\section{Monte Carlo uncertainty}
\label{sec:MC_err}
Since \HEJ is reweighting each Fixed Order point with multiple resummation
events, the Monte Carlo uncertainty of \HEJ is a little bit more complicated
then usual. We start by defining the \HEJ cross section after $N$ FO points
\begin{align}
\sigma_N:=\sum_{i}^N x_i \sum_{j}^{M_i} y_{i,j}=:\sum_i^N\sum_j^{M_i} w_{i,j},
\end{align}
where $x_i$ are the FO weights\footnote{In this definition $x_i$ can be zero,
see the discussion in the next section.}, $y_{i,j}$ are the reweighting weights
, and $M_i$ the number of resummation points. We can set $M=M_i \forall i$ by
potentially adding some points with $y_{i,j}=0$, i.e. $M$ correspond to the
\lstinline!trials! in \lstinline!EventReweighter!. $w_{i,j}$ are the weights as
written out by \HEJ. The expectation value of $\sigma$ is then
\begin{align}
\ev{\sigma_N}= \sum_i \ev{x_i}\sum_j\ev{y_{i,j}}=M \mu_x\sum_i\mu_{y_i},\label{eq:true_sigma}
\end{align}
with $\mu_{x/y}$ being the (true) mean value of $x$ or $y$, i.e.
\begin{align}
\mu_{x}:=\ev{\bar{x}}=\ev{\frac{\sum_i x_i}{N}}=\ev{x}.
\end{align}
The true underlying standard derivation on $\sigma_N$, assuming $\delta_{x}$
and $\delta_{y_i}$ are the standard derivations of $x$ and $y_i$ is
\begin{align}
\delta_{\sigma_N}^2&=M^2 \delta_{x}^2 \sum_i \mu_{y_i}^2
+M \mu_x^2 \sum_i \delta_{y_i}^2. \label{eq:true_err}
\end{align}
Notice that each point $i$ can have an different expectation for $y_i$.
Since we do not know the true distribution of $x$ and $y$ we need to estimate
it. We use the standard derivation
\begin{align}
\tilde{\delta}_{x_i}^2&:=\left(x_i-\bar x\right)^2
=\left(\frac{N-1}{N} x_i - \frac{\sum_{j\neq i} x_j}{N}\right)^2
\label{eq:err_x}\\
\tilde{\delta}_{y_{i,j}}^2&:=\left(y_{i,j}-\bar y_i\right)^2 \label{eq:err_y},
\end{align}
and the mean values $\bar x$ and $\bar y$, to get an estimator for
$\delta_{\sigma_N}$
\begin{align}
\tilde\delta_{\sigma_N}^2&=M^2 \sum_i \tilde\delta_{x_i}^2 \bar{y_i}^2
+\sum_{i,j} x_i^2\tilde\delta_{y_{i,j}}^2. \label{eq:esti_err}
\end{align}
Trough error propagation we can connect the estimated uncertainties back to the
fundamental ones
\begin{align}
\delta_{\tilde{\delta}_{x_i}}^2=\frac{N-1}{N} \delta_x^2.
\end{align}
Together with $\delta_x^2=\ev{x^2}-\ev{x}^2$ and $\ev{\tilde\delta}=0$ this
leads to
\begin{align}
\ev{\tilde{\delta}_{x_i}^2 \bar y_i^2}&=\ev{\tilde{\delta}_{x_i} \bar y_i}^2
+\delta_{\tilde{\delta}_{x_i}}^2 \mu_{y_i}^2
+\delta_{y_i}^2 \mu_{\tilde\delta}^2 \\
&=\frac{N-1}{N} \delta_x^2\mu_{y_i}^2,
\end{align}
and a similar results for $y$. Therefore
\begin{align}
\ev{\delta_{\sigma_N}}=\frac{N-1}{N} M^2 \delta_{x}^2 \sum_i \mu_{y_i}^2
+\frac{M-1}{M} M \mu_x^2 \sum_i \delta_{y_i}^2,
\end{align}
where we can compensate for the additional factors compared to~\eqref{eq:true_err}, by replacing
\begin{align}
\tilde\delta_x&\to\frac{N}{N-1}\tilde\delta_x \label{eq:xcom_bias}\\
\tilde\delta_{y_i}&\to\frac{M}{M-1}\tilde\delta_{y_i}. \label{eq:ycom_bias}
\end{align}
Thus~\eqref{eq:esti_err} is an unbiased estimator of $\delta_{\sigma_N}$.
\subsection{Number of events vs. number of trials}
Even though the above calculation is completely valid, it is unpractical. Both
$x_i$ and $y_{ij}$ could be zero, but zero weight events are typically not
written out. In that sense $N$ and $M$ are the \textit{number of trials} it took
to generate $N'$ and $M'$ (non-zero) events. We can not naively replace all $N$
and $M$ with $N'$ and $M'$ in the above equations, since this would also change
the definition of the average $\bar x$ and $\bar y$.
For illustration let us consider unweighted events, with all weights equal to
$x'$, without changing the cross section $\sum_i^N x_i=\sum_i^{N'} x'_i=N' x'$.
Then the average trial weight is unequal to the average event weight
\begin{align}
\bar x = \frac{\sum_i^{N} x_i}{N} = \frac{\sum_i^{N'} x'}{N}=x'\frac{N'}{N}
\neq x'=\frac{\sum_i^{N'} x'}{N'}.
\end{align}
$N=N'$ would correspond to an $100\%$ efficient unweighting, i.e. a perfect
sampling, where we know the analytical results. In particular using $N'$ instead
of $N$ in the standard derivation gives
\begin{align}
\sum_i \left(x_i-\frac{\sum_i^{N} x_i}{N'}\right)^2=\sum_i \left(x'-x' \frac{\sum_i^{N'}}{N'}\right)^2=0,
\end{align}
which is obviously not true in general for $\tilde\delta^2_x$.
Hence we would have to use the number of trials $N$ everywhere. This would
require an additional parameter to be passed with each events, which is not
always available in practice\footnote{ \texttt{Sherpa} gives the number of
trials, as an \lstinline!attribute::trials! of \lstinline!HEPEUP! in the
\texttt{LHE} file, or similarly as a data member in the HDF5 format
\cite{Hoeche:2019rti}. The \texttt{LHE} standard itself provides the
variable \lstinline!ntries! per event (see
\href{https://phystev.cnrs.fr/wiki/2017:groups:tools:lhe}{this proposal}),
though I have not seen this used anywhere.}. Instead we use
\begin{align}
\tilde\delta_{x}'^2:=\sum_i^{N} x_i^2\geq\tilde\delta_x^2, \label{eq:err_prac}
\end{align}
where the bias of $\delta_x'^2$ vanishes for large $N$. Thus we can use the sum
of weight squares~\eqref{eq:err_prac} instead of~\eqref{eq:err_x}
and~\eqref{eq:err_y}, without worrying about the difference between trials and
generated events. The total error~\eqref{eq:esti_err} becomes
\begin{align}
\tilde\delta_{\sigma_N}^2=\sum_i \left(\sum_j w_{i,j}\right)^2+\sum_{i,j} \left(w_{i,j}\right)^2,
\end{align}
which (conveniently) only dependent on the \HEJ weights $w_{i,j}$.
\section{Explicit formulas for vector currents}
\label{sec:j_vec}
Using eqs.~\eqref{eq:u_plus}\eqref{eq:u_minus}\eqref{eq:jminus_gen}\eqref{eq:Pauli_matrices}\eqref{eq:jplus_gen}, the vector currents read
\begin{align}
\label{eq:j-_explicit}
j^-_\mu(p, q) ={}&
\begin{pmatrix}
\sqrt{p^+\,q^+} + \sqrt{p^-\,q^-} \hat{p}_{\perp} \hat{q}_{\perp}^*\\
\sqrt{p^-\,q^+}\, \hat{p}_{\perp} + \sqrt{p^+\,q^-}\,\hat{q}_{\perp}^*\\
-i \sqrt{p^-\,q^+}\, \hat{p}_{\perp} + i \sqrt{p^+\,q^-}\, \hat{q}_{\perp}^*\\
\sqrt{p^+\,q^+} - \sqrt{p^-\,q^-}\, \hat{p}_{\perp}\, \hat{q}_{\perp}^*
\end{pmatrix}\,,\\
j^+_\mu(p, q) ={}&\big(j^-_\mu(p, q)\big)^*\,,\\
j^\pm_\mu(q, p) ={}&\big(j^\pm_\mu(p, q)\big)^*\,.
\end{align}
If $q= p_{\text{in}}$ is the momentum of an incoming parton, we have
$\hat{p}_{\text{in} \perp} = -1$ and either $p_{\text{in}}^+ = 0$ or
$p_{\text{in}}^- = 0$. The current simplifies further:\todo{Helicities flipped w.r.t code}
\begin{align}
\label{eq:j_explicit}
j^-_\mu(p_{\text{out}}, p_{\text{in}}) ={}&
\begin{pmatrix}
\sqrt{p_{\text{in}}^+\,p_{\text{out}}^+}\\
\sqrt{p_{\text{in}}^+\,p_{\text{out}}^-} \ \hat{p}_{\text{out}\,\perp}\\
-i\,j^-_1\\
j^-_0
\end{pmatrix}
& p_{\text{in}\,z} > 0\,,\\
j^-_\mu(p_{\text{out}}, p_{\text{in}}) ={}&
\begin{pmatrix}
-\sqrt{p_{\text{in}}^-\,p_{\text{out}}^{-\phantom{+}}} \ \hat{p}_{\text{out}\,\perp}\\
- \sqrt{p_{\text{in}}^-\,p_{\text{out}}^+}\\
i\,j^-_1\\
-j^-_0
\end{pmatrix} & p_{\text{in}\,z} < 0\,.
\end{align}
We also employ the usual short-hand notation
For the gluon polarisation vectors with gluon momentum $p_g$ and auxiliary
reference vector $p_r$ we use
\begin{equation}
\label{eq:pol_vector}
\epsilon_\mu^+(p_g, p_r) = \frac{j_\mu^+(p_r, p_g)}{\sqrt{2}\spb g.r}\,,\qquad\epsilon_\mu^-(p_g, p_r) = \frac{j_\mu^-(p_r, p_g)}{\sqrt{2}\spa g.r}\,.
\end{equation}
\bibliographystyle{JHEP}
\bibliography{biblio}
\end{document}
diff --git a/include/HEJ/LorentzVector.hh b/include/HEJ/LorentzVector.hh
index 0b5ed8e..df42d12 100644
--- a/include/HEJ/LorentzVector.hh
+++ b/include/HEJ/LorentzVector.hh
@@ -1,16 +1,39 @@
+/** \file
+ * \brief Auxiliary functions for Lorentz vectors
+ *
+ * \authors The HEJ collaboration (see AUTHORS for details)
+ * \date 2019
+ * \copyright GPLv2 or later
+ */
+#pragma once
+
#include "CLHEP/Vector/LorentzVector.h"
namespace HEJ {
inline
auto dot(
CLHEP::HepLorentzVector const & h1,
CLHEP::HepLorentzVector const & h2
) {
return h1.dot(h2);
}
inline
auto m2(CLHEP::HepLorentzVector const & h1) {
return h1.m2();
}
+
+ //! Split a single Lorentz vector into two lightlike Lorentz vectors
+ /**
+ * @param P Lorentz vector to be split
+ * @returns A pair (p, q) of Lorentz vectors with P = p + q and p^2 = q^2 = 0
+ *
+ * P.perp() has to be positive.
+ *
+ * p.e() is guaranteed to be positive.
+ * In addition, if either of P.plus() or P.minus() is positive,
+ * q.e() has the same sign as P.m2()
+ */
+ std::pair<CLHEP::HepLorentzVector, CLHEP::HepLorentzVector>
+ split_into_lightlike(CLHEP::HepLorentzVector const & P);
}
diff --git a/src/LorentzVector.cc b/src/LorentzVector.cc
new file mode 100644
index 0000000..745bf6b
--- /dev/null
+++ b/src/LorentzVector.cc
@@ -0,0 +1,33 @@
+/**
+ * \authors The HEJ collaboration (see AUTHORS for details)
+ * \date 2020
+ * \copyright GPLv2 or later
+ */
+#include "HEJ/LorentzVector.hh"
+
+#include <cmath>
+#include <utility>
+
+namespace HEJ {
+ // see eq:P_massive, eq:P_massive_p, eq:P_massive_q, , eq:P_massive_plus,
+ // eq:P_massive_minus in the developer manual
+ std::pair<CLHEP::HepLorentzVector, CLHEP::HepLorentzVector>
+ split_into_lightlike(CLHEP::HepLorentzVector const & P) {
+ const double pt = P.perp();
+ if(P.plus() > 0){
+ const double y = std::log(P.plus()/pt);
+ const double E = P.m2()/(2.*P.plus());
+ return std::make_pair(
+ CLHEP::HepLorentzVector{P.x(), P.y(), pt*std::sinh(y), pt*std::cosh(y)},
+ CLHEP::HepLorentzVector{0., 0., E, -E}
+ );
+ } else {
+ const double y = std::log(pt/P.minus());
+ const double E = P.m2()/(2.*P.minus());
+ return std::make_pair(
+ CLHEP::HepLorentzVector{P.x(), P.y(), pt*std::sinh(y), pt*std::cosh(y)},
+ CLHEP::HepLorentzVector{0., 0., E, +E}
+ );
+ }
+ }
+}
diff --git a/src/Wjets.cc b/src/Wjets.cc
index 55cd707..2d36888 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1345 +1,1140 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <array>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/jets.hh"
#include "HEJ/Tensor.hh"
+#include "HEJ/LorentzVector.hh"
+#include "HEJ/exceptions.hh"
// generated headers
#include "HEJ/currents/jW_uno.hh"
+#include "HEJ/currents/W_central_qqbar.hh"
using HEJ::Tensor;
using HEJ::init_sigma_index;
using HEJ::metric;
using HEJ::rank3_current;
using HEJ::rank5_current;
using HEJ::eps;
using HEJ::to_tensor;
using HEJ::Helicity;
using HEJ::angle;
using HEJ::square;
using HEJ::flip;
using HEJ::ParticleProperties;
namespace helicity = HEJ::helicity;
namespace { // Helper Functions
// FKL W Helper Functions
double WProp (const HLV & plbar, const HLV & pl, ParticleProperties const & wprop){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
+ COM(0.,1.)*wprop.mass*wprop.width);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
namespace {
// FKL current including W emission off negative helicities
// See eq. (87) {eq:jW-} in developer manual
// Note that the terms are rearranged
Tensor<1> jW_minus(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl
){
using HEJ::helicity::minus;
const double tWin = (pa-pl-plbar).m2();
const double tWout = (p1+pl+plbar).m2();
// C++ arithmetic operators are evaluated left-to-right,
// so the following first computes complex scalar coefficients,
// which then multiply a current, reducing the number
// of multiplications
return 2.*(
+ angle(p1, pl)*square(p1, plbar)/tWout
+ square(pa, plbar)*angle(pa, pl)/tWin
)*HEJ::current(p1, pa, helicity::minus)
+ 2.*angle(p1, pl)*square(pl, plbar)/tWout
*HEJ::current(pl, pa, helicity::minus)
+ 2.*square(pa, plbar)*angle(pl, plbar)/tWin
*HEJ::current(p1, plbar, helicity::minus);
}
}
// FKL current including W emission
// see eqs. (87), (88) {eq:jW-}, {eq:jW+} in developer manual
Tensor<1> jW(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl,
Helicity h
){
if(h == helicity::minus) {
return jW_minus(pa, p1, plbar, pl);
}
return jW_minus(pa, p1, pl, plbar).complex_conj();
}
/**
* @brief Contraction of the \tilde{U}_1 tensor
*
* This is the contraction of the tensor defined in eq:U1tensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM U1contr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
using helicity::plus;
using helicity::minus;
if(h1 == plus) {
if(h2 == plus) {
if(hg == plus) {
return HEJ::U1<plus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::U1<plus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
else {
if(hg == helicity::plus) {
return HEJ::U1<plus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::U1<plus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
return HEJ::U1<minus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::U1<minus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
else {
if(hg == helicity::plus) {
return HEJ::U1<minus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::U1<minus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
}
}
/**
* @brief Contraction of the \tilde{U}_1 tensor
*
* This is the contraction of the tensor defined in eq:U2tensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM U2contr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
using helicity::plus;
using helicity::minus;
if(h1 == plus) {
if(h2 == plus) {
if(hg == plus) {
return HEJ::U2<plus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::U2<plus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
else {
if(hg == helicity::plus) {
return HEJ::U2<plus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::U2<plus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
return HEJ::U2<minus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::U2<minus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
else {
if(hg == helicity::plus) {
return HEJ::U2<minus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::U2<minus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
}
}
/**
* @brief Contraction of the \tilde{U}_1 tensor
*
* This is the contraction of the tensor defined in eq:Ltensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM Lcontr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
using helicity::plus;
using helicity::minus;
if(h1 == plus) {
if(h2 == plus) {
if(hg == plus) {
return HEJ::L<plus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::L<plus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
else {
if(hg == helicity::plus) {
return HEJ::L<plus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::L<plus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
return HEJ::L<minus,plus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::L<minus,plus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
else {
if(hg == helicity::plus) {
return HEJ::L<minus,minus,plus>(p1,p2,pa,pb,pg,pl,plbar);
}
else {
return HEJ::L<minus,minus,minus>(p1,p2,pa,pb,pg,pl,plbar);
}
}
}
}
/**
* @brief W+Jets Unordered Contribution Helper Functions
* @returns result of equation (4.1.28) in Helen's Thesis (p.100)
*/
double jM2Wuno(HLV pg, HLV p1,HLV plbar, HLV pl, HLV pa, Helicity h1,
HLV p2, HLV pb, Helicity h2, Helicity pol,
ParticleProperties const & wprop
){
//@TODO Simplify the below (less Tensor class?)
init_sigma_index();
HLV pW = pl+plbar;
HLV q1g=pa-pW-p1-pg;
HLV q1 = pa-p1-pW;
HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
const double s2g = (p2+pg).m2();
const double sbg = (pb+pg).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
//use p1 as ref vec in pol tensor
Tensor<1> epsg = eps(pg,p2,pol);
Tensor<1> epsW = HEJ::current(pl,plbar,helicity::minus);
Tensor<1> j2b = HEJ::current(p2,pb,h2);
Tensor<1> Tq1q2 = to_tensor(
(pb/sbg + p2/s2g)*(q1 - pg).m2() + 2*q1 - pg
);
Tensor<3> J31a = rank3_current(p1, pa, h1);
Tensor<2> J2_qaW =J31a.contract((pa-pW)/taW/taW1, 2);
Tensor<2> J2_p1W =J31a.contract((p1+pW)/s1W/taW1, 2);
Tensor<3> L1a = outer(Tq1q2, J2_qaW);
Tensor<3> L1b = outer(Tq1q2, J2_p1W);
Tensor<3> L2a = outer(-2*pg, J2_qaW);
Tensor<3> L2b = outer(-2*pg, J2_p1W);
Tensor<3> L3 = outer(metric(), J2_qaW.contract(2*pg-q1,1)+J2_p1W.contract(2*pg-q1,2));
Tensor<3> L(0.);
Tensor<5> J51a = rank5_current(p1, pa, h1);
Tensor<4> J_qaW = J51a.contract((pa-pW),4);
Tensor<4> J_qag = J51a.contract(pa-pg,4);
Tensor<4> J_p1gW = J51a.contract(p1+pg+pW,4);
Tensor<3> U1a = J_qaW.contract(p1+pg,2);
Tensor<3> U1b = J_p1gW.contract(p1+pg,2);
Tensor<3> U1c = J_p1gW.contract(p1+pW,2);
Tensor<3> U1(0.);
Tensor<3> U2a = J_qaW.contract(pa-pg-pW,2);
Tensor<3> U2b = J_qag.contract(pa-pg-pW,2);
Tensor<3> U2c = J_qag.contract(p1+pW,2);
Tensor<3> U2(0.);
for(int nu=0; nu<4;++nu){
for(int mu=0;mu<4;++mu){
for(int rho=0;rho<4;++rho){
L(nu, mu, rho) = L1a(nu,mu,rho) + L1b(nu,rho,mu)
+ L2a(mu,nu,rho) + L2b(mu,rho,nu) + L3(mu,nu,rho);
U1(nu, mu, rho) = U1a(nu, mu, rho) / (s1g*taW)
+ U1b(nu,rho,mu) / (s1g*s1gW) + U1c(rho,nu,mu) / (s1W*s1gW);
U2(nu,mu,rho) = U2a(mu,nu,rho) / (taWg*taW)
+ U2b(mu,rho,nu) / (taWg*tag) + U2c(rho,mu,nu) / (s1W*tag);
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
double amp = HEJ::C_A*HEJ::C_F*HEJ::C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
double t1 = q1g.m2();
double t2 = q2.m2();
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
//Divide by t-channels
amp/=(t1*t2);
return amp;
}
/**
* @brief New implementation of unordered W+Jets current
*
* See eq:wunocurrent in the developer manual for the definition
* of this current
*
* The aim is to replace jM2Wuno and eventually all uses of the Tensor class.
* Explicit tensor contractions are rather slow and the initialisation code
* in Tensor.cc is not very transparent.
*
* At the moment, this function _only_ works for positive-energy spinors,
* so jM2Wuno has to be kept for qqbar emission.
*/
double jM2Wuno_pos_energy(
HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity const h1,
HLV const & p2, HLV const & pb, Helicity const h2, Helicity const pol,
ParticleProperties const & wprop
){
using HEJ::C_A;
using HEJ::C_F;
const COM U1 = U1contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM U2 = U2contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM L = Lcontr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM X = U1 - L;
const COM Y = U2 + L;
double amp = C_A*C_F*C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
amp *= WProp(plbar, pl, wprop);
const HLV q1g = pa-pl-plbar-p1-pg;
const HLV q2 = p2-pb;
const double t1 = q1g.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
- // Relevant Wqqx Helper Functions.
- //g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
- Tensor <1> gtqqxW(HLV pq,HLV pqbar,HLV pl,HLV plbar){
- //@TODO Simplify the calculation below (Less Tensor class use?)
- double s2AB=(pl+plbar+pq).m2();
- double s3AB=(pl+plbar+pqbar).m2();
-
- // Define llx current.
- Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
-
- //blank 3 Gamma Current
- Tensor<3> JV23 = rank3_current(pq,pqbar,helicity::minus);
-
- // Components of g->qqW before W Contraction
- Tensor<2> JV1 = JV23.contract((pq + pl + plbar),2)/(s2AB);
- Tensor<2> JV2 = JV23.contract((pqbar + pl + plbar),2)/(s3AB);
-
- // g->qqW Current. Note Minus between terms due to momentum flow.
- // Also note: (-I)^2 from W vert. (I) from Quark prop.
- Tensor<1> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
-
- return JVCur;
- }
-
- // Helper Functions Calculate the Crossed Contribution
- Tensor <2> MCrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar, HLV pl,
- HLV plbar, std::vector<HLV> partons, int nabove
- ){
- //@TODO Simplify the calculation below Maybe combine with MCross?
- // Useful propagator factors
- double s2AB=(pl+plbar+pq).m2();
- double s3AB=(pl+plbar+pqbar).m2();
-
- HLV q1, q3;
- q1=pa;
- for(int i=0; i<nabove+1;++i){
- q1=q1-partons.at(i);
- }
- q3 = q1 - pq - pqbar - pl - plbar;
-
- double tcro1=(q3+pq).m2();
- double tcro2=(q1-pqbar).m2();
-
- // Define llx current.
- Tensor<1> ABCur = HEJ::current(pl, plbar,helicity::minus);
-
- //Blank 5 gamma Current
- Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
-
- // 4 gamma currents (with 1 contraction already).
- Tensor<4> J_q3q = J523.contract((q3 + pq),2);
- Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
-
- // Components of Crossed Vertex Contribution
- Tensor<3> Xcro1 = J_q3q.contract((pqbar + pl + plbar),3);
- Tensor<3> Xcro2 = J_q3q.contract((q1 - pqbar),3);
- Tensor<3> Xcro3 = J_2AB.contract((q1 - pqbar),3);
-
- // Term Denominators Taken Care of at this stage
- Tensor<2> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
- Tensor<2> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
- Tensor<2> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
-
- //Initialise the Crossed Vertex Object
- Tensor<2> Xcro(0.);
-
- for(int mu=0; mu<4;++mu){
- for(int nu=0;nu<4;++nu){
- Xcro(mu,nu) = -(-Xcro1Cont(nu,mu)+Xcro2Cont(nu,mu)+Xcro3Cont(nu,mu));
- }
- }
-
- return Xcro;
- }
-
- // Helper Functions Calculate the Uncrossed Contribution
- Tensor <2> MUncrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar,
- HLV pl, HLV plbar, std::vector<HLV> partons, int nabove
- ){
- //@TODO Simplify the calculation below Maybe combine with MUncross?
- double s2AB=(pl+plbar+pq).m2();
- double s3AB=(pl+plbar+pqbar).m2();
-
- HLV q1, q3;
- q1=pa;
- for(int i=0; i<nabove+1;++i){
- q1=q1-partons.at(i);
- }
- q3 = q1 - pl - plbar - pq - pqbar;
- double tunc1 = (q1-pq).m2();
- double tunc2 = (q3+pqbar).m2();
-
- // Define llx current.
- Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
-
- //Blank 5 gamma Current
- Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
-
- // 4 gamma currents (with 1 contraction already).
- Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
- Tensor<4> J_q1q = J523.contract((q1 - pq),2);
-
- // 2 Contractions taken care of.
- Tensor<3> Xunc1 = J_2AB.contract((q3 + pqbar),3);
- Tensor<3> Xunc2 = J_q1q.contract((q3 + pqbar),3);
- Tensor<3> Xunc3 = J_q1q.contract((pqbar + pl + plbar),3);
-
- // Term Denominators Taken Care of at this stage
- Tensor<2> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
- Tensor<2> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
- Tensor<2> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
-
- //Initialise the Uncrossed Vertex Object
- Tensor<2> Xunc(0.);
-
- for(int mu=0; mu<4;++mu){
- for(int nu=0;nu<4;++nu){
- Xunc(mu,nu) = -(- Xunc1Cont(mu,nu)+Xunc2Cont(mu,nu) +Xunc3Cont(mu,nu));
- }
- }
-
- return Xunc;
- }
-
- // Helper Functions Calculate the g->qqxW (Eikonal) Contributions
- Tensor <2> MSymW(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
- HLV pl,HLV plbar, std::vector<HLV> partons, int nabove
- ){
- //@TODO Simplify the calculation below Maybe combine with MSym?
- double sa2=(pa+pq).m2();
- double s12=(p1+pq).m2();
- double sa3=(pa+pqbar).m2();
- double s13=(p1+pqbar).m2();
- double saA=(pa+pl).m2();
- double s1A=(p1+pl).m2();
- double saB=(pa+plbar).m2();
- double s1B=(p1+plbar).m2();
- double sb2=(pb+pq).m2();
- double s42=(p4+pq).m2();
- double sb3=(pb+pqbar).m2();
- double s43=(p4+pqbar).m2();
- double sbA=(pb+pl).m2();
- double s4A=(p4+pl).m2();
- double sbB=(pb+plbar).m2();
- double s4B=(p4+plbar).m2();
- double s23AB=(pl+plbar+pq+pqbar).m2();
-
- HLV q1,q3;
- q1=pa;
- for(int i=0;i<nabove+1;++i){
- q1-=partons.at(i);
- }
- q3=q1-pq-pqbar-plbar-pl;
- double t1 = (q1).m2();
- double t3 = (q3).m2();
-
- // g->qqW Current (Factors of sqrt2 dealt with in this function.)
- Tensor<1> JV = gtqqxW(pq,pqbar,pl,plbar);
-
- // 1a gluon emisson Contribution
- Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13+s1A+s1B))
- + pa*(t1/(sa2+sa3+saA+saB)) );
- Tensor<2> X1aCont = X1a.contract(JV,3);
-
- //4b gluon emission Contribution
- Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43+s4A+s4B))
- + pb*(t3/(sb2+sb3+sbA+sbB)) );
- Tensor<2> X4bCont = X4b.contract(JV,3);
-
- //Set up each term of 3G diagram.
- Tensor<3> X3g1 = outer(q1+pq+pqbar+pl+plbar, metric());
- Tensor<3> X3g2 = outer(q3-pq-pqbar-pl-plbar, metric());
- Tensor<3> X3g3 = outer(q1+q3, metric());
-
- // Note the contraction of indices changes term by term
- Tensor<2> X3g1Cont = X3g1.contract(JV,3);
- Tensor<2> X3g2Cont = X3g2.contract(JV,2);
- Tensor<2> X3g3Cont = X3g3.contract(JV,1);
-
- // XSym is an amalgamation of x1a, X4b and X3g.
- // Makes sense from a colour factor point of view.
- Tensor<2>Xsym(0.);
-
- for(int mu=0; mu<4;++mu){
- for(int nu=0;nu<4;++nu){
- Xsym(mu,nu) = (X3g1Cont(nu,mu) + X3g2Cont(mu,nu) - X3g3Cont(nu,mu))
- + (X1aCont(mu,nu) - X4bCont(mu,nu));
- }
- }
- return Xsym/s23AB;
- }
-
Tensor <2> MCross(HLV pa, HLV pq, HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2> XCroCont = J323.contract((q1-pqbar),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xcro(mu,nu) = XCroCont(nu,mu);
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncross(HLV pa, HLV pq,HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2> XUncCont = J323.contract((q1-pq),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xunc(mu,nu) = -XUncCont(mu,nu);
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2> MSym(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
std::vector<HLV> partons, Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MsymW?
HLV q1, q3;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3 = q1-pq-pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
double s23 = (pq+pqbar).m2();
double sa2 = (pa+pq).m2();
double sa3 = (pa+pqbar).m2();
double s12 = (p1+pq).m2();
double s13 = (p1+pqbar).m2();
double sb2 = (pb+pq).m2();
double sb3 = (pb+pqbar).m2();
double s42 = (p4+pq).m2();
double s43 = (p4+pqbar).m2();
Tensor<1> qqxCur = HEJ::current(pq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13))+ pa*(t1/(sa2+sa3)));
Tensor<2> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43)) + pb*(t3/(sb2+sb3)));
Tensor<2> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3> X3g1 = outer(q1+pq+pqbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xsym(mu, nu) = COM(0,1) * ( (X3g1Cont(nu,mu) + X3g2Cont(mu,nu)
- X3g3Cont(nu,mu)) + (X1aCont(mu,nu) - X4bCont(mu,nu)) );
}
}
return Xsym/s23;
}
//! W+Jets FKL Contributions
/**
* @brief W+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_\mu.
* Handles all possible incoming states.
*/
double jW_j( HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
bool aqlineb, bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out);
const double WPropfact = WProp(plbar, pl, wprop);
const auto j_W = COM{0,-1}*jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
double Msqr = 0.;
for(const auto h: {plus, minus}) {
const auto j = HEJ::current(p2out, p2in, h);
Msqr += abs2(j_W.contract(j, 1));
}
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*Msqr/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
} // Anonymous Namespace
double ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop);
}
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true, wprop);
}
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop);
}
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true, wprop);
}
double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param pg Unordered Gluon momenta
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_{uno}_\mu. Ie, unordered with W emission opposite side.
* Handles all possible incoming states.
*/
double jW_juno(HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out,
HLV p2in, HLV pg, bool aqlineb, bool aqlinef,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out-pg);
const HLV q3=-(p2in-p2out);
const Helicity fhel = aqlinef?plus:minus;
const auto j_W = jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
const auto mj2p = HEJ::current(p2out, p2in, flip(fhel));
const auto mj2m = HEJ::current(p2out, p2in, fhel);
const auto jgbp = HEJ::current(pg, p2in, flip(fhel));
const auto jgbm = HEJ::current(pg, p2in, fhel);
const auto j2gp = HEJ::current(p2out, pg, flip(fhel));
const auto j2gm = HEJ::current(p2out, pg, fhel);
// Dot products of these which occur again and again
COM MWmp=j_W.dot(mj2p); // And now for the Higgs ones
COM MWmm=j_W.dot(mj2m);
const auto qsum = to_tensor(q2+q3);
const auto p1o = to_tensor(p1out);
const auto p1i = to_tensor(p1in);
const auto p2o = to_tensor(p2out);
const auto p2i = to_tensor(p2in);
const auto Lmm=( (-1.)*qsum*(MWmm) + (-2.*COM{j_W.dot(pg)})*mj2m + 2.*COM{mj2m.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmm/2. ) )/q3.m2();
const auto Lmp=( (-1.)*qsum*(MWmp) + (-2.*COM{j_W.dot(pg)})*mj2p + 2.*COM{mj2p.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmp/2. ) )/q3.m2();
const auto U1mm=(COM{jgbm.dot(j_W)}*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
const auto U1mp=(COM{jgbp.dot(j_W)}*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
const auto U2mm=((-1.)*COM{j2gm.dot(j_W)}*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
const auto U2mp=((-1.)*COM{j2gp.dot(j_W)}*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
double amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
ampsq*=WProp(plbar, pl, wprop);
return ampsq/((16)*(q2.m2()*q1.m2()));
}
}
double ME_W_unob_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true, wprop);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false, wprop);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true, wprop);
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same side.
* Handles all possible incoming states. Note this handles both forward and back-
* -ward Wuno emission. For forward, ensure p1out is the uno and W emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean up
* ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(HLV pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
HLV p2out, HLV p2in, bool aqlineb,
ParticleProperties const & wprop
){
//Calculate different Helicity choices
const Helicity h = aqlineb?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::plus,helicity::plus, wprop);
double ME2mpm = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::plus,helicity::minus, wprop);
double ME2mmp = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::minus,helicity::plus, wprop);
double ME2mmm = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states
return (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
}
}
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Extremal qqx Contributions, function to handle all incoming types.
* @param pgin Incoming gluon which will split into qqx.
* @param pqout Quark of extremal qqx outgoing (W-Emission).
* @param plbar Outgoing anti-lepton momenta
* @param pl Outgoing lepton momenta
* @param pqbarout Anti-quark of extremal qqx pair. (W-Emission)
* @param pout Outgoing Particle 2 (end of FKL chain)
* @param p2in Incoming Particle 2
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W_{qqx} ^\mu j_\mu. Ie, Ex-QQX with W emission same side.
* Handles all possible incoming states. Calculated via crossing symmetry from jWuno_j.
*/
double jWqqx_j(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef,
ParticleProperties const & wprop
){
//Calculate Different Helicity Configurations.
const Helicity h = aqlinef?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::plus,helicity::plus, wprop);
double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::plus,helicity::minus, wprop);
double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::minus,helicity::plus, wprop);
double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states.
double ME2 = (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
}
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
}
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm1(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below. (Less Tensor class use?)
double t1 = (p3-pb)*(p3-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p3-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,2)/t1;
return gqqCur*(-1);
}
//Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm2(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double t1 = (p2-pb)*(p2-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb,refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p2-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm3(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s23 = (p2+p3)*(p2+p3);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank1 = outer(p2+p3, metric())/s23;
Tensor<3> qqCurBlank2 = outer(pb, metric())/s23;
Tensor<1> Cur23 = HEJ::current(p2, p3,hel2);
Tensor<2> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2,
HLV p3,HLV plbar, HLV pl, bool aqlinepa,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
// 2 independent helicity choices (complex conjugation related).
Tensor<1> TMmmm1 = qggm1(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,minus,minus, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,minus,plus, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa?plus:minus);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1);
COM Mmmm2 = TMmmm2.contract(cur1a,1);
COM Mmmm3 = TMmmm3.contract(cur1a,1);
COM Mpmm1 = TMpmm1.contract(cur1a,1);
COM Mpmm2 = TMpmm2.contract(cur1a,1);
COM Mpmm3 = TMpmm3.contract(cur1a,1);
//Colour factors:
COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
cm1m1=8./3.;
cm2m2=8./3.;
cm3m3=6.;
cm1m2 =-1./3.;
cm1m3 = -3.*COM(0.,1.);
cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
double Mmmm = real( cm1m1*pow(abs(Mmmm1),2) + cm2m2*pow(abs(Mmmm2),2)
+ cm3m3*pow(abs(Mmmm3),2) + 2.*real(cm1m2*Mmmm1*conj(Mmmm2))
+ 2.*real(cm1m3*Mmmm1*conj(Mmmm3))
+ 2.*real(cm2m3*Mmmm2*conj(Mmmm3)) );
double Mpmm = real( cm1m1*pow(abs(Mpmm1),2) + cm2m2*pow(abs(Mpmm2),2)
+ cm3m3*pow(abs(Mpmm3),2) + 2.*real(cm1m2*Mpmm1*conj(Mpmm2))
+ 2.*real(cm1m3*Mpmm1*conj(Mpmm3))
+ 2.*real(cm2m3*Mpmm2*conj(Mpmm3)) );
// Divide by WProp
const double WPropfact = WProp(plbar, pl, wprop);
return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
}
-// W+Jets qqxCentral
-double ME_WCenqqx_qq(HLV pa, HLV pb,HLV pl, HLV plbar, std::vector<HLV> partons,
- bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
- ParticleProperties const & wprop
-){
- init_sigma_index();
+namespace {
+ // helper function for matrix element for W + Jets with central qqbar
+ template<HEJ::Helicity h1a, HEJ::Helicity h4b>
+ double amp_WCenqqx_qq(
+ HLV const & pa, HLV const & p1,
+ HLV const & pb, HLV const & p4,
+ HLV const & pq, HLV const & pqbar,
+ HLV const & pl, HLV const & plbar,
+ HLV const & q11, HLV const & q12
+ ) {
+ const COM sym = HEJ::M_sym_W<h1a, h4b>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
+ );
+ const COM cross = HEJ::M_cross_W<h1a, h4b>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
+ );
+ const COM uncross = HEJ::M_uncross_W<h1a, h4b>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
+ );
+
+ // Colour factors
+ static constexpr COM cmsms = 3.;
+ static constexpr COM cmumu = 4./3.;
+ static constexpr COM cmcmc = 4./3.;
+ static constexpr COM cmsmu = COM{0., 3./2.};
+ static constexpr COM cmsmc = COM{0.,-3./2.};
+ static constexpr COM cmumc = -1./6.;
+
+ return real(
+ cmsms*pow(abs(sym),2)
+ +cmumu*pow(abs(uncross),2)
+ +cmcmc*pow(abs(cross),2)
+ )
+ +2.*real(cmsmu*sym*conj(uncross))
+ +2.*real(cmsmc*sym*conj(cross))
+ +2.*real(cmumc*uncross*conj(cross))
+ ;
+ }
+}
+
+// matrix element for W + Jets with central qqbar
+double ME_WCenqqx_qq(
+ HLV pa, HLV pb, HLV pl, HLV plbar, std::vector<HLV> partons,
+ bool /* aqlinepa */, bool /* aqlinepb */, bool qqxmarker, int nabove,
+ ParticleProperties const & wprop
+) {
+ using helicity::plus;
+ using helicity::minus;
+
+ CLHEP::HepLorentzVector q1 = pa;
+ for(int i = 0; i <= nabove; ++i){
+ q1 -= partons[i];
+ }
+ const auto qq = HEJ::split_into_lightlike(q1);
+ // since q1.m2() < 0 the following assertion is always true
+ // see documentation for split_into_lightlike
+ assert(qq.second.e() < 0);
- HLV pq, pqbar, p1, p4;
+ HLV pq, pqbar;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
- pqbar = partons[nabove+2];}
-
- p1 = partons.front();
- p4 = partons.back();
-
- Tensor<1> T1am, T4bm, T1ap, T4bp;
- if(!(aqlinepa)){
- T1ap = HEJ::current(p1, pa, helicity::plus);
- T1am = HEJ::current(p1, pa, helicity::minus);}
- else if(aqlinepa){
- T1ap = HEJ::current(pa, p1, helicity::plus);
- T1am = HEJ::current(pa, p1, helicity::minus);}
- if(!(aqlinepb)){
- T4bp = HEJ::current(p4, pb, helicity::plus);
- T4bm = HEJ::current(p4, pb, helicity::minus);}
- else if(aqlinepb){
- T4bp = HEJ::current(pb, p4, helicity::plus);
- T4bm = HEJ::current(pb, p4, helicity::minus);}
-
- // Calculate the 3 separate contributions to the effective vertex
- Tensor<2> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
- Tensor<2> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
- Tensor<2> Xsym = MSymW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
+ pqbar = partons[nabove+2];
+ }
+ const HLV p1 = partons.front();
+ const HLV p4 = partons.back();
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
- // (- - hel choice)
- COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1));
- COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1));
- COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1));
- // (- + hel choice)
- COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1));
- COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1));
- COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1));
- // (+ - hel choice)
- COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1));
- COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1));
- COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1));
- // (+ + hel choice)
- COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1));
- COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1));
- COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1));
-
- //Colour factors:
- COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
- cmsms=3.;
- cmumu=4./3.;
- cmcmc=4./3.;
- cmsmu =3./2.*COM(0.,1.);
- cmsmc = -3./2.*COM(0.,1.);
- cmumc = -1./6.;
-
- // Work Out Interference in each case of helicity:
- double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
- +cmumu*pow(abs(M_mmUnc),2)
- +cmcmc*pow(abs(M_mmCro),2)
- +2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
- +2.*real(cmsmc*M_mmSym*conj(M_mmCro))
- +2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
-
- double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
- +cmumu*pow(abs(M_mpUnc),2)
- +cmcmc*pow(abs(M_mpCro),2)
- +2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
- +2.*real(cmsmc*M_mpSym*conj(M_mpCro))
- +2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
-
- double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
- +cmumu*pow(abs(M_pmUnc),2)
- +cmcmc*pow(abs(M_pmCro),2)
- +2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
- +2.*real(cmsmc*M_pmSym*conj(M_pmCro))
- +2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
-
- double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
- +cmumu*pow(abs(M_ppUnc),2)
- +cmcmc*pow(abs(M_ppCro),2)
- +2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
- +2.*real(cmsmc*M_ppSym*conj(M_ppCro))
- +2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
-
- double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
-
- HLV q1,q3;
- q1=pa;
- for(int i=0;i<nabove+1;++i){
- q1-=partons.at(i);
- }
- q3 = q1 - pq - pqbar - pl - plbar;
-
- double t1 = (q1).m2();
- double t3 = (q3).m2();
-
- //Divide by t-channels
- amp/=(t1*t1*t3*t3);
-
- //Divide by WProp
- amp*=WProp(plbar, pl, wprop);
+ // the first helicity label is for aqlinepa == true,
+ // the second one for aqlinepb == true
+ // In principle the corresponding helicity should be flipped
+ // if either of them is false, but we sum up all helicities,
+ // so this has no effect in the end.
+ const double amp_mm = amp_WCenqqx_qq<minus, minus>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
+ );
+ const double amp_mp = amp_WCenqqx_qq<minus, plus>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
+ );
+ const double amp_pm = amp_WCenqqx_qq<plus, minus>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
+ );
+ const double amp_pp = amp_WCenqqx_qq<plus, plus>(
+ pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
+ );
+
+ const double t1 = q1.m2();
+ const double t3 = (q1-pl-plbar-pq-pqbar).m2();
+
+ const double amp = WProp(plbar, pl, wprop)*(
+ amp_mm+amp_mp+amp_pm+amp_pp
+ )/(9.*4.*t1*t1*t3*t3);
return amp;
}
// no wqq emission
double ME_W_Cenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
int nbelow, bool forwards, ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
if (!forwards){ //If Emission from Leg a instead, flip process.
std::swap(pa, pb);
std::reverse(partons.begin(),partons.end());
std::swap(aqlinepa, aqlinepb);
qqxmarker = !qqxmarker;
std::swap(nabove,nbelow);
}
HLV pq, pqbar, p1,p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, plus);
T1am = HEJ::current(p1, pa, minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, plus);
T1am = HEJ::current(pa, p1, minus);}
Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb?plus:minus);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, minus, nabove);
Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, plus, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
return amp;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 3:01 PM (5 h, 33 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486801
Default Alt Text
(191 KB)

Event Timeline