Page MenuHomeHEPForge
Files F8309399

No OneTemporary
Actions

View Options

diff --git a/doc/sphinx/HEJFOG.rst b/doc/sphinx/HEJFOG.rst
index 1a72676..70046d7 100644
--- a/doc/sphinx/HEJFOG.rst
+++ b/doc/sphinx/HEJFOG.rst
@@ -1,309 +1,311 @@
The HEJ Fixed Order Generator
=============================
For high jet multiplicities event generation with standard fixed-order
generators becomes increasingly cumbersome. For example, the
leading-order production of a Higgs Boson with five or more jets is
computationally prohibitively expensive.
To this end, HEJ 2 provides the ``HEJFOG`` fixed-order generator
that allows to generate events with high jet multiplicities. To
facilitate the computation the limit of Multi-Regge Kinematics with
large invariant masses between all outgoing particles is assumed in the
matrix elements. The typical use of the ``HEJFOG`` is to supplement
low-multiplicity events from standard generators with high-multiplicity
events before using the HEJ 2 program to add high-energy
resummation.
Installation
------------
The ``HEJFOG`` comes bundled together with HEJ 2 and the
installation is very similar. After downloading HEJ 2 and
installing the prerequisites as described in :ref:`Installation` the
``HEJFOG`` can be installed with::
cmake /path/to/FixedOrderGen -DCMAKE_INSTALL_PREFIX=target/directory
make install
where :file:`/path/to/FixedOrderGen` refers to the :file:`FixedOrderGen`
subdirectory in the HEJ 2 directory. If HEJ 2 was
installed to a non-standard location, it may be necessary to specify the
directory containing :file:`HEJ-config.cmake`. If the base installation
directory is :file:`/path/to/HEJ`, :file:`HEJ-config.cmake` should be
found in :file:`/path/to/HEJ/lib/cmake/HEJ` and the commands for
installing the ``HEJFOG`` would read::
cmake /path/to/FixedOrderGen -DHEJ_DIR=/path/to/HEJ/lib/cmake/HEJ -DCMAKE_INSTALL_PREFIX=target/directory
make install
The installation can be tested with::
make test
provided that the CT10nlo PDF set is installed.
Running the fixed-order generator
---------------------------------
After installing the ``HEJFOG`` you can modify the provided
configuration file :file:`configFO.yml` and run the generator with::
HEJFOG configFO.yml
The resulting event file, by default named :file:`HEJFO.lhe`, can then be
fed into HEJ 2 like any event file generated from a standard
fixed-order generator, see :ref:`Running HEJ 2`.
Settings
--------
Similar to HEJ 2, the ``HEJFOG`` uses a `YAML
<http://yaml.org/>`_ configuration file. The settings are
.. _`process`:
**process**
The scattering process for which events are being generated. The
format is
:code:`in1 in2 => out1 out2 ...`
The incoming particles, :code:`in1`, :code:`in2` can be
- quarks: :code:`u`, :code:`d`, :code:`u_bar`, and so on
- gluons: :code:`g`
- protons :code:`p` or antiprotons :code:`p_bar`
The outgoing particles, can be
- jets: :code:`j`, multiple jets can be grouped together e.g. :code:`2j`
- bosons: Any gauge boson, they require `particle properties`_
- - leptons: if they can be reconstructed to a :code:`W+` or :code:`W-`, e.g. :code:`e+ nu_e`
+ - leptons: if they can be reconstructed to a :code:`W+` or :code:`W-`,
+ e.g. :code:`e+ nu_e`
At most one of the outgoing particles can be a boson, the rest has to be
partonic. At the moment only Higgs :code:`h`, :code:`W+` and :code:`W-`
bosons are supported. All other outgoing particles are jets. There have to be
at least two jets. Instead of the leptons decays of the bosons can be added
through the :ref:`particle properties<particle properties: particle:
decays>`. The latter is required to decay a Higgs boson. So :code:`p p => Wp
j j` is the same as :code:`p p => e+ nu_e 2j`, if the decay :code:`Wp => e+
nu_e` is specified in `particle properties`_.
.. _`events`:
**events**
Specifies the number of events to generate.
.. _`jets`:
**jets**
Defines the properties of the generated jets.
.. _`jets: min pt`:
**min pt**
Minimum jet transverse momentum in GeV.
.. _`jets: peak pt`:
**peak pt**
Optional setting to specify the dominant jet transverse momentum
in GeV. If the generated events are used as input for HEJ
resummation, this should be set to the minimum transverse momentum
of the resummation jets. In all cases it has to be larger than
:code:`min pt`. The effect is that only a small fraction of jets
will be generated with a transverse momentum below the value of
this setting.
.. _`jets: algorithm`:
**algorithm**
The algorithm used to define jets. Allowed settings are
:code:`kt`, :code:`cambridge`, :code:`antikt`,
:code:`cambridge for passive`. See the `FastJet
<http://fastjet.fr/>`_ documentation for a description of these
algorithms.
.. _`jets: R`:
**R**
The R parameter used in the jet algorithm.
.. _`jets: max rapidity`:
**max rapidity**
Maximum absolute value of the jet rapidity.
.. _`beam`:
**beam**
Defines various properties of the collider beam.
.. _`beam: energy`:
**energy**
The beam energy in GeV. For example, the 13
TeV LHC corresponds to a value of 6500.
.. _`beam: particles`:
**particles**
A list :code:`[p1, p2]` of two beam particles. The only supported
entries are protons :code:`p` and antiprotons :code:`p_bar`.
.. _`pdf`:
**pdf**
The `LHAPDF number <https://lhapdf.hepforge.org/pdfsets>`_ of the PDF set.
For example, 230000 corresponds to an NNPDF 2.3 NLO PDF set.
.. _`subleading fraction`:
**subleading fraction**
This setting is related to the fraction of events that are not a FKL
configuration. All possible subleading process are listed in
:ref:`subleading channels<subleading channels>`. This value must be
positive and not larger than 1. It should typically be chosen between
0.01 and 0.1. Note that while this parameter influences the chance of
generating subleading configurations, it generally does not
correspond to the actual fraction of subleading events.
.. _`subleading channels`:
**subleading channels**
Optional parameters to select a specific subleading process, multiple
channels can be selected at once. If multiple subleading configurations
are possible one will be selected at random for each event, thus each
final state will include at most one subleading process, e.g. either
:code:`unordered` or :code:`qqx`. Only has an effect if
:code:`subleading fraction` is non-zero. The following values are allowed:
- :code:`all`: All channels allowed. This is the default.
- :code:`none`: No subleading contribution, only the leading (FKL)
configurations are allowed. This is equivalent to
:code:`subleading fraction: 0`.
- :code:`unordered`: Unordered emission allowed.
Unordered emission are any rapidity ordering where exactly one gluon is
emitted outside the rapidity ordering required in FKL events. More
precisely, if at least one of the incoming particles is a quark or
antiquark and there are more than two jets in the final state,
:code:`subleading fraction` states the probability that the flavours
of the outgoing particles are assigned in such a way that an unordered
configuration arises. Unordered emissions are currently not implemented
for pure jet production.
- :code:`qqx`: Production of a quark-antiquark pair. In the leading
(FKL) configurations all partons except for the most backward and
forward ones are gluons. If the :code:`qqx` channel is allowed,
:code:`subleading fraction` gives the probability of emitting a
quark-antiquark pair in place of two gluons that are adjacent in
rapidity. Quark-antiquark pairs are currently only implemented for
W boson production.
.. _`unweight`:
**unweight**
This setting defines the parameters for the partial unweighting of
events. You can disable unweighting by removing this entry from the
configuration file.
.. _`unweight: sample size`:
**sample size**
The number of weighted events used to calibrate the unweighting.
A good default is to set this to the number of target
`events`_. If the number of `events`_ is large this can
lead to significant memory consumption and a lower value should be
chosen. Contrarily, for large multiplicities the unweighting
efficiency becomes worse and the sample size should be increased.
.. _`unweight: max deviation`:
**max deviation**
Controls the range of events to which unweighting is applied. A
larger value means that a larger fraction of events are unweighted.
Typical values are between -1 and 1.
.. _`particle properties`:
**particle properties**
Specifies various properties of the different particles (Higgs, W or Z).
This is only relevant if the chosen `process`_ is the production of the
corresponding particles with jets. E.g. for the `process`_
:code:`p p => h 2j` the :code:`mass`, :code:`width` and (optionally)
:code:`decays` of the :code:`Higgs` boson are required, while all other
particle properties will be ignored.
.. _`particle properties: particle`:
**Higgs, W+, W- or Z**
The particle (Higgs, |W+|, |W-|, Z) for which the following
properties are defined.
.. |W+| replace:: W\ :sup:`+`
.. |W-| replace:: W\ :sup:`-`
.. _`particle properties: particle: mass`:
**mass**
The mass of the particle in GeV.
.. _`particle properties: particle: width`:
**width**
The total decay width of the particle in GeV.
.. _`particle properties: particle: decays`:
**decays**
Optional setting specifying the decays of the particle. Only the decay
into two particles is implemented. Each decay has the form
:code:`{into: [p1,p2], branching ratio: r}`
where :code:`p1` and :code:`p2` are the particle names of the
decay product (e.g. :code:`photon`) and :code:`r` is the branching
ratio. The branching ratio is optional (:code:`1` by default).
Decays of a Higgs boson are treated as the production and subsequent
decay of an on-shell Higgs boson, so decays into e.g. Z bosons are not
- supported. In contrast W bosons are decayed off-shell, thus the branching ratio should not be changed or set to :code:`1`.
+ supported. In contrast W bosons are decayed off-shell, thus the
+ branching ratio should not be changed or set to :code:`1`.
.. _`scales`:
**scales**
Specifies the renormalisation and factorisation scales for the output
events. For details, see the corresponding entry in the HEJ 2
:ref:`HEJ 2 settings`. Note that this should usually be a
single value, as the weights resulting from additional scale choices
will simply be ignored when adding high-energy resummation with HEJ 2.
.. _`event output`:
**event output**
Specifies the name of a single event output file or a list of such
files. See the corresponding entry in the HEJ 2
:ref:`HEJ 2 settings` for details.
.. _`RanLux init`:
.. _`random generator`:
**random generator**
Sets parameters for random number generation. See the HEJ 2
:ref:`HEJ 2 settings` for details.
.. _`analysis`:
**analysis**
Specifies the name and settings for a custom analysis library. This
can be useful to specify cuts at the fixed-order level. See the
corresponding entry in the HEJ 2 :ref:`HEJ 2 settings`
for details.
.. _`Higgs coupling`:
**Higgs coupling**
This collects a number of settings concerning the effective coupling
of the Higgs boson to gluons. See the corresponding entry in the
HEJ 2 :ref:`HEJ 2 settings` for details.
diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst
index 7e4cc5b..24c5f27 100644
--- a/doc/sphinx/installation.rst
+++ b/doc/sphinx/installation.rst
@@ -1,94 +1,97 @@
.. _Installation:
Installation
============
Download
--------
A tar archive of the HEJ 2 source code can be downloaded and
decompressed with the command::
curl https://hej.web.cern.ch/HEJ/downloads/HEJ_2.0.tar.gz | tar -xz
To obtain the latest stable HEJ version, `HEJ_2.0.tar.gz` should be
replaced by `HEJ.tar.gz`.
Alternatively, the HEJ source code can be obtained by installing the
`git version control system <https://git-scm.com/>`_. and running::
git clone https://phab.hepforge.org/source/hej.git
We also provide a `Docker image <https://hub.docker.com/r/hejdock/hej>`_
containing a HEJ 2 installation on. This image can be pulled with::
docker pull hejdock/hej
When using the Docker image the remaining installation steps can be
skipped.
Prerequisites
-------------
Before installing HEJ 2, you need the following programs and
libraries:
- `CMake <https://cmake.org/>`_ version 3.1
-- A compiler supporting the C++14 standard, for example `gcc <https://gcc.gnu.org/>`_ 5 or later
+- A compiler supporting the C++14 standard, for example
+ `gcc <https://gcc.gnu.org/>`_ 5 or later
- `FastJet <http://fastjet.fr/>`_
- `CLHEP <https://gitlab.cern.ch/CLHEP/CLHEP>`_
- `LHAPDF <https://lhapdf.hepforge.org/>`_
- The `IOStreams` and `uBLAS` `boost <https://www.boost.org>`_ libraries
- `yaml-cpp <https://github.com/jbeder/yaml-cpp>`_
In addition, some optional features have additional dependencies:
- `Version 2 of QCDLoop <https://github.com/scarrazza/qcdloop>`_ is
required to include finite top mass corrections in Higgs boson + jets
production.
- `HepMC versions 2 and 3 <https://hepmc.web.cern.ch/hepmc/>`_ enable
event output in the respective format.
- `Rivet <https://rivet.hepforge.org/>`_ together with HepMC 2 allow
using Rivet analyses.
- `HighFive <https://github.com/BlueBrain/HighFive>`_ has to be
installed in order to read event files in the
`HDF5 <https://www.hdfgroup.org/>`_-based format suggested in
`arXiv:1905.05120 <https://arxiv.org/abs/1905.05120>`_.
Compilation
-----------
To compile and install HEJ 2 run::
cmake source/directory -DCMAKE_INSTALL_PREFIX=target/directory
make install
:file:`source/directory` is the directory containing the file
:file:`CMakeLists.txt`. If you omit
:code:`-DCMAKE_INSTALL_PREFIX=target/directory` HEJ 2 will be
installed to some default location.
In case some of the aforementioned prerequisites are not found by
:code:`cmake` you can give a hint by adding an additional argument
:code:`-Dlibname_ROOT_DIR=/directory/with/library`, where
:code:`libname` should be replaced by the name of the library in
question.
To not include specific packages one can add
:code:`-DEXCLUDE_packagename=TRUE` to :code:`cmake`, e.g. by
setting :code:`-DEXCLUDE_rivet=TRUE` HEJ 2 will not be interfaced
to `rivet` even if rivet is available on the system.
Testing
-------
To test your installation, download the NNPDF 2.3 PDF set with::
lhapdf install NNPDF23_nlo_as_0119
and run::
make test
The test data of HEJ are sorted in a
-`Git Large File Storage <https://git-lfs.github.com/>`_ format. :code:`git clone` therefore requires :code:`git-lfs` to download the data correctly.
+`Git Large File Storage <https://git-lfs.github.com/>`_ format.
+:code:`git clone` therefore requires :code:`git-lfs` to download the data
+correctly.
diff --git a/include/HEJ/Constants.hh b/include/HEJ/Constants.hh
index 30856b5..19be946 100644
--- a/include/HEJ/Constants.hh
+++ b/include/HEJ/Constants.hh
@@ -1,39 +1,41 @@
/** \file
* \brief Header file defining all global constants used for HEJ
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
namespace HEJ{
/// @name QCD parameters
//@{
constexpr double N_C = 3.; //!< number of Colours
constexpr double C_A = N_C; //!< \f$C_A\f$
constexpr double C_F = (N_C*N_C - 1.)/(2.*N_C); //!< \f$C_F\f$
constexpr double t_f = 0.5; //!< \f$t_f\f$
constexpr double n_f = 5.; //!< number light flavours
constexpr double beta0 = 11./3.*C_A - 4./3.*t_f*n_f; //!< \f$\beta_0\f$
//@}
/// @name QFT parameters
//@{
constexpr double vev = 246.2196508; //!< Higgs vacuum expectation value in GeV
constexpr double gw = 0.653233; //!< elector-weak coupling
constexpr double MW = 80.419; //!< The W mass in GeV/c^2
constexpr double GammaW = 2.0476; //!< the W width in GeV/c^2
//@}
/// @name Generation Parameters
//@{
- //! Default scale for virtual correction, \f$\lambda\f$ cf. eq. (20) in \cite Andersen:2011hs
+ //! @brief Default scale for virtual correction
+ //! \f$\lambda\f$ cf. eq. (20) in \cite Andersen:2011hs
constexpr double CLAMBDA = 0.2;
constexpr double CMINPT = 0.2; //!< minimal \f$p_t\f$ of all partons
//@}
/// @name Conventional Parameters
//@{
- //! Value of first colour for colour dressing, according to LHE convention \cite Boos:2001cv
+ //! Value of first colour for colour dressing, according to LHE convention
+ //! \cite Boos:2001cv
constexpr int COLOUR_OFFSET = 501;
//@}
}
diff --git a/include/HEJ/CrossSectionAccumulator.hh b/include/HEJ/CrossSectionAccumulator.hh
index f66bb7d..a839d0c 100644
--- a/include/HEJ/CrossSectionAccumulator.hh
+++ b/include/HEJ/CrossSectionAccumulator.hh
@@ -1,74 +1,75 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <map>
#include <ostream>
#include <iomanip>
#include <string>
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
namespace HEJ {
template<typename T>
struct XSWithError {
T value = T{};
T error = T{};
};
/**
* @brief Sum of Cross Section for different subproccess
*/
class CrossSectionAccumulator {
public:
void fill(HEJ::Event const & ev) {
const double wt = ev.central().weight;
auto & entry = xs_[ev.type()];
entry.value += wt;
entry.error += wt*wt;
total_.value += wt;
total_.error += wt*wt;
}
auto begin() const {
return std::begin(xs_);
}
auto end() const {
return std::end(xs_);
}
//! total Cross Section and error
XSWithError<double> total() const {
return total_;
}
private:
std::map<HEJ::event_type::EventType, XSWithError<double>> xs_;
XSWithError<double> total_;
};
std::ostream& operator<<(std::ostream& os, const CrossSectionAccumulator& xs){
const std::streamsize orig_prec = os.precision();
os << std::scientific << std::setprecision(3)
<< " " << std::left << std::setw(25)
<< "Cross section: " << xs.total().value
<< " +- " << std::sqrt(xs.total().error) << " (pb)\n";
for(auto const & xs_type: xs) {
os << " " << std::left << std::scientific <<std::setw(25)
<< (event_type::name(xs_type.first) + std::string(": "));
os << xs_type.second.value << " +- "
<< std::sqrt(xs_type.second.error) << " (pb) "
<< std::fixed << std::setprecision(3)
- << "[" <<std::setw(6) <<std::right<< ((xs_type.second.value)/xs.total().value)*100 << "%]"
+ << "[" <<std::setw(6) <<std::right
+ << ((xs_type.second.value)/xs.total().value)*100 << "%]"
<< std::endl;
}
os << std::defaultfloat;
os.precision(orig_prec);
return os;
}
}
diff --git a/include/HEJ/event_types.hh b/include/HEJ/event_types.hh
index e5be748..3ea75da 100644
--- a/include/HEJ/event_types.hh
+++ b/include/HEJ/event_types.hh
@@ -1,106 +1,106 @@
/** \file
* \brief Define different types of events.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include "HEJ/utility.hh"
namespace HEJ{
//! Namespace for event types
namespace event_type{
//! Possible event types
enum EventType: size_t{
- non_resummable=0, /**< event configuration not covered by All Order resummation */
- bad_final_state=1, /**< event with an unsupported final state */
- no_2_jets=2, /**< event with less than two jets */
- FKL=4, /**< FKL-type event */
- unordered_backward=8, /**< event with unordered backward emission */
- unordered_forward=16, /**< event with unordered forward emission */
- extremal_qqxb=32, /**< event with a backward extremal qqbar */
- extremal_qqxf=64, /**< event with a forward extremal qqbar */
- central_qqx=128, /**< event with a central qqbar */
- unob = unordered_backward,
- unof = unordered_forward,
- qqxexb = extremal_qqxb,
- qqxexf = extremal_qqxf,
- qqxmid = central_qqx,
- first_type = non_resummable,
- last_type = central_qqx
+ non_resummable=0, //!< event configuration not covered by All Order resummation
+ bad_final_state=1, //!< event with an unsupported final state
+ no_2_jets=2, //!< event with less than two jets
+ FKL=4, //!< FKL-type event
+ unordered_backward=8, //!< event with unordered backward emission
+ unordered_forward=16, //!< event with unordered forward emission
+ extremal_qqxb=32, //!< event with a backward extremal qqbar
+ extremal_qqxf=64, //!< event with a forward extremal qqbar
+ central_qqx=128, //!< event with a central qqbar
+ unob = unordered_backward, //!< alias for unordered_backward
+ unof = unordered_forward, //!< alias for unordered_forward
+ qqxexb = extremal_qqxb, //!< alias for extremal_qqxb
+ qqxexf = extremal_qqxf, //!< alias for extremal_qqxf
+ qqxmid = central_qqx, //!< alias for central_qqx
+ first_type = non_resummable, //!< alias for non_resummable
+ last_type = central_qqx //!< alias for central_qqx
};
//! Event type names
/**
* For example, name(FKL) is the string "FKL"
*/
inline std::string name(EventType type) {
switch(type) {
case FKL:
return "FKL";
case unordered_backward:
return "unordered backward";
case unordered_forward:
return "unordered forward";
case extremal_qqxb:
return "extremal qqbar backward";
case extremal_qqxf:
return "extremal qqbar forward";
case central_qqx:
return "central qqbar";
case non_resummable:
return "non-resummable";
case no_2_jets:
return "no 2 jets";
case bad_final_state:
return "bad final state";
default:
throw std::logic_error{"Unreachable"};
}
};
//! Returns True for a HEJ \ref event_type::EventType "EventType"
inline
bool is_resummable(EventType type) {
switch(type) {
case FKL:
case unordered_backward:
case unordered_forward:
case extremal_qqxb:
case extremal_qqxf:
case central_qqx:
return true;
default:
return false;
}
}
//! Returns True for an unordered \ref event_type::EventType "EventType"
inline
bool is_uno(EventType type) {
return type == unordered_backward || type == unordered_forward;
}
//! Returns True for an extremal_qqx \ref event_type::EventType "EventType"
inline
bool is_ex_qqx(EventType type) {
return type == extremal_qqxb || type == extremal_qqxf;
}
//! Returns True for an central_qqx \ref event_type::EventType "EventType"
inline
bool is_mid_qqx(EventType type) {
return type == central_qqx;
}
//! Returns True for any qqx event \ref event_type::EventType "EventType"
inline
bool is_qqx(EventType type) {
return is_ex_qqx(type) || is_mid_qqx(type);
}
} // namespace event_type
} // namespace HEJ
diff --git a/src/Hjets.cc b/src/Hjets.cc
index 82fa2e2..d00baab 100644
--- a/src/Hjets.cc
+++ b/src/Hjets.cc
@@ -1,1104 +1,1117 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Hjets.hh"
#include <assert.h>
#include <limits>
#include "HEJ/Constants.hh"
#ifdef HEJ_BUILD_WITH_QCDLOOP
#include "qcdloop/qcdloop.h"
#endif
const COM looprwfactor = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4);
constexpr double infinity = std::numeric_limits<double>::infinity();
namespace {
// Loop integrals
#ifdef HEJ_BUILD_WITH_QCDLOOP
COM B0DD(HLV q, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_B0 = [](){
ql::Bubble<std::complex<double>,double,double> ql_B0;
ql_B0.setCacheSize(100);
return ql_B0;
}();
static std::vector<double> masses(2);
static std::vector<double> momenta(1);
for(auto & m: masses) m = mq*mq;
momenta.front() = q.m2();
ql_B0.integral(result, 1, masses, momenta);
return result[0];
}
COM C0DD(HLV q1, HLV q2, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_C0 = [](){
ql::Triangle<std::complex<double>,double,double> ql_C0;
ql_C0.setCacheSize(100);
return ql_C0;
}();
static std::vector<double> masses(3);
static std::vector<double> momenta(3);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = (q1+q2).m2();
ql_C0.integral(result, 1, masses, momenta);
return result[0];
}
COM D0DD(HLV q1,HLV q2, HLV q3, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_D0 = [](){
ql::Box<std::complex<double>,double,double> ql_D0;
ql_D0.setCacheSize(100);
return ql_D0;
}();
static std::vector<double> masses(4);
static std::vector<double> momenta(6);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = q3.m2();
momenta[3] = (q1+q2+q3).m2();
momenta[4] = (q1+q2).m2();
momenta[5] = (q2+q3).m2();
ql_D0.integral(result, 1, masses, momenta);
return result[0];
}
COM A1(HLV q1, HLV q2, double mt)
// As given in Eq. (B.2) of VDD
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
assert(mt > 0.);
mt2=mt*mt;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 4.*mt2/Delta3*(Q2-q12-q22)
-1.-4.*q12*q22/Delta3-12.*q12*q22*Q2/Delta3/Delta3*(q12+q22-Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q2,mt)-B0DD(Q,mt) )
* ( 2.*q22/Delta3+12.*q12*q22/Delta3/Delta3*(q22-q12+Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q1,mt)-B0DD(Q,mt) )
* ( 2.*q12/Delta3+12.*q12*q22/Delta3/Delta3*(q12-q22+Q2) )
- 2./Delta3/16/M_PI/M_PI*(q12+q22-Q2);
return ans;
}
COM A2(HLV q1, HLV q2, double mt)
// As given in Eq. (B.2) of VDD, but with high energy limit
// of invariants taken.
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
assert(mt > 0.);
mt2=mt*mt;
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 2.*mt2+1./2.*(q12+q22-Q2)
+2.*q12*q22*Q2/Delta3 )
+looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))
*q22*(q22-q12-Q2)/Delta3
+looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))
*q12*(q12-q22-Q2)/Delta3+1./16/M_PI/M_PI;
return ans;
}
#else // no QCDloop
COM A1(HLV, HLV, double) {
throw std::logic_error{"A1 called without QCDloop support"};
}
COM A2(HLV, HLV, double) {
throw std::logic_error{"A2 called without QCDloop support"};
}
#endif
void to_current(const HLV & q, current & ret){
ret[0]=q.e();
ret[1]=q.x();
ret[2]=q.y();
ret[3]=q.z();
}
/**
* @brief Higgs vertex contracted with current @param C1 and @param C2
*/
COM cHdot(const current & C1, const current & C2, const current & q1,
const current & q2, double mt, bool incBot, double mb)
{
if (mt == infinity) {
return (cdot(C1,C2)*cdot(q1,q2)-cdot(C1,q2)*cdot(C2,q1))/(3*M_PI*HEJ::vev);
}
else {
HLV vq1,vq2;
vq1.set(q1[1].real(),q1[2].real(),q1[3].real(),q1[0].real());
vq2.set(q2[1].real(),q2[2].real(),q2[3].real(),q2[0].real());
// first minus sign obtained because of q1-difference to VDD
// Factor is because 4 mt^2 g^2/HEJ::vev A1 -> 16 pi mt^2/HEJ::vev alphas,
if(!(incBot))
- return 16.*M_PI*mt*mt/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt));
+ return 16.*M_PI*mt*mt/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)
+ -cdot(C1,C2)*A2(-vq1,vq2,mt));
else
- return 16.*M_PI*mt*mt/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt))
- + 16.*M_PI*mb*mb/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)-cdot(C1,C2)*A2(-vq1,vq2,mb));
+ return 16.*M_PI*mt*mt/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)
+ -cdot(C1,C2)*A2(-vq1,vq2,mt))
+ + 16.*M_PI*mb*mb/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)
+ -cdot(C1,C2)*A2(-vq1,vq2,mb));
}
}
//@{
/**
* @brief Higgs+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @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 q1 t-channel momenta into higgs vertex
* @param q2 t-channel momenta out of higgs vertex
* @param mt top mass (inf or value)
* @param incBot Bool, to include bottom mass (true) or not (false)?
* @param mb bottom mass (value)
* @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^\mu H j_\mu. FKL with higgs vertex somewhere in the FKL chain.
* Handles all possible incoming states.
*/
double j_h_j(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV q1, HLV q2, double mt, bool incBot,
double mb, bool aqlineb, bool aqlinef){
current j1p,j1m,j2p,j2m, q1v, q2v;
// Note need to flip helicities in anti-quark case.
joi(p1out,!aqlineb, p1in,!aqlineb, j1p);
joi(p1out, aqlineb, p1in, aqlineb, j1m);
joi(p2out,!aqlinef, p2in,!aqlinef, j2p);
joi(p2out, aqlinef, p2in, aqlinef, j2m);
to_current(q1, q1v);
to_current(q2, q2v);
COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);
// average over helicities
const double sst=(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm))/4.;
return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
}
} // namespace anonymous
double ME_H_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, false);
}
double ME_H_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, true);
}
double ME_H_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, true, false);
}
double ME_H_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, true, true);
}
double ME_H_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, false)*K_g(p2out,p2in)/HEJ::C_A;
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, false)
+ * K_g(p2out,p2in)/HEJ::C_A;
}
double ME_H_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, true, false)*K_g(p2out,p2in)/HEJ::C_A;
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, true, false)
+ * K_g(p2out,p2in)/HEJ::C_A;
}
double ME_H_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, false)*(K_g(p2out,p2in)/HEJ::C_A)*(K_g(p1out,p1in)/HEJ::C_A);
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, false)
+ * K_g(p2out,p2in)/HEJ::C_A * K_g(p1out,p1in)/HEJ::C_A;
}
//@}
namespace {
//@{
/// @brief Higgs vertex contracted with one current
CCurrent jH(HLV pout, bool helout, HLV pin,
bool helin, HLV q1, HLV q2,
double mt, bool incBot, double mb)
{
CCurrent j2 = joi(pout,helout,pin,helin);
CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());
if(mt == infinity)
return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*HEJ::vev);
else
{
if(incBot)
- return (-16.*M_PI*mb*mb/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mb)-16.*M_PI*mb*mb/HEJ::vev*j2*A2(-q1,q2,mb))
- + (-16.*M_PI*mt*mt/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/HEJ::vev*j2*A2(-q1,q2,mt));
+ return (-16.*M_PI*mb*mb/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mb)
+ -16.*M_PI*mb*mb/HEJ::vev*j2*A2(-q1,q2,mb))
+ + (-16.*M_PI*mt*mt/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)
+ -16.*M_PI*mt*mt/HEJ::vev*j2*A2(-q1,q2,mt));
else
- return (-16.*M_PI*mt*mt/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/HEJ::vev*j2*A2(-q1,q2,mt));
+ return (-16.*M_PI*mt*mt/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)
+ -16.*M_PI*mt*mt/HEJ::vev*j2*A2(-q1,q2,mt));
}
}
//@}
//@{
/**
* @brief Higgs+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (W emission)
* @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 q1 t-channel momenta into higgs vertex
* @param q2 t-channel momenta out of higgs vertex
* @param mt top mass (inf or value)
* @param incBot Bool, to include bottom mass (true) or not (false)?
* @param mb bottom mass (value)
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_{uno}^\mu H j_\mu. Unordered with higgs vertex
* somewhere in the FKL chain. Handles all possible incoming states.
*/
double j_h_juno(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1, HLV qH2,
double mt, bool incBot, double mb, bool aqlineb
){
// This construction is taking rapidity order: pg > p1out >> p2out
HLV q1=p1in-p1out; // Top End
HLV q2=-(p2in-p2out); // Bottom End
HLV qg=p1in-p1out-pg; // Extra bit post-gluon
CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p,jgam,jgap,j2gm,j2gp;
CCurrent test1, test2, test3, test4;
// Note <p1|eps|pa> current split into two by gauge choice.
// See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
mj1p=joi(p1out,!aqlineb, p1in,!aqlineb);
mj1m=joi(p1out, aqlineb, p1in, aqlineb);
jgap=joi(pg, !aqlineb, p1in,!aqlineb);
jgam=joi(pg, aqlineb, p1in, aqlineb);
// Note for function joo(): <p1+|pg+> = <pg-|p1->.
j2gp=joo(p1out, !aqlineb, pg, !aqlineb);
j2gm=joo(p1out, aqlineb, pg, aqlineb);
mjH2p=jH(p2out, true,p2in, true,qH1,qH2,mt,incBot,mb);
mjH2m=jH(p2out,false,p2in,false,qH1,qH2,mt,incBot,mb);
// Dot products of these which occur again and again
COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
COM MHmm=mj1m.dot(mjH2m);
COM MHpp=mj1p.dot(mjH2p);
COM MHpm=mj1p.dot(mjH2m);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm;
CCurrent p2o(p2out), p2i(p2in), p1o(p1out), p1i(p1in), qsum(q1+qg);
Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmm/2.) )/q1.m2();
Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmp/2.) )/q1.m2();
Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpm/2.) )/q1.m2();
Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpp/2.) )/q1.m2();
U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2();
U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2();
U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2();
U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2();
U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2();
double amm,amp,apm,app;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
app=HEJ::C_F*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());
// Now add the t-channels for the Higgs
double th=qH1.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
// Factor of (Cf/Ca) for each quark to match ME_H_qQ.
ampsq*=HEJ::C_F*HEJ::C_F/HEJ::C_A/HEJ::C_A;
return ampsq;
}
} // namespace anonymous
double ME_H_unof_qQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1,
HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, false);
}
double ME_H_unof_qbarQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, true);
}
double ME_H_unof_qQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, false);
}
double ME_H_unof_qbarQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, true);
}
double ME_H_unof_qg(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, false)*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_H_unof_qbarg(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, true)*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_H_unob_qQ(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, false);
}
double ME_H_unob_qbarQ(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, false);
}
double ME_H_unob_qQbar(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, true);
}
double ME_H_unob_qbarQbar(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, true);
}
double ME_H_unob_gQ(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, false)*K_g(p1out,p1in)/HEJ::C_F;
}
double ME_H_unob_gQbar(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, true)*K_g(p1out,p1in)/HEJ::C_F;
}
//@}
// Begin finite mass stuff
#ifdef HEJ_BUILD_WITH_QCDLOOP
namespace {
// All the stuff needed for the box functions in qg->qgH now...
COM E1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2=-(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) +
S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) +
2.*(s34 + S1)*(s34 + S1)/Delta +
S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
2.*(s34 + S1)/Delta +
S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) -
C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) -
C0DD(k1 + k2, q2, mq)*2.*s34/
S1 - (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 +
S1)/(S1*Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(2.*s34*(s34 +
S1)*(S1 - S2)/(Delta*Sigma) +
2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S1)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM F1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-S2*D0DD(k1, k2, q2,
mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) -
s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
- S2*pow((s34 + S2),2)/Delta/Delta) - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12) -
- C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) -
+ S2*pow((s34 + S2),2)/Delta/Delta)
+ - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12)
+ - C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2,
mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 +
S2)*(S2 - S1)/(Delta*Sigma) + (B0DD(
k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S2)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM G1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor*(S2*D0DD(k1, q2, k2,
mq)*(Delta/s12/s12 - 4.*mq*mq/s12) -
S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) -
C0DD(k1, q2, mq))*4.*mq*mq/
s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./
s12 + (B0DD(k1 + q2, mq) -
B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1)));
}
COM E4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k2, k1, q2,
mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) -
s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S1),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta +
2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) +
- s12*(S1 -
- S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) -
- B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*((2.*s12*(2.*s12*s34 - S1*(S1 + S2) +
- s34*(S2 - S1)))/(Delta*Sigma)));
+ s12*(S1 - S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) -
+ B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))
+ *((2.*s12*(2.*s12*s34 - S1*(S1 + S2) + s34*(S2 - S1)))/(Delta*Sigma)));
}
COM F4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k1, k2, q2,
mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) +
s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S2),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 +
S2)/Delta + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 -
S1*(S1 + S2) +
s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) -
- B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(2.*s12*(2.*s12*s34 - S2*(S1 + S2) +
- s34*(S1 - S2))/(Delta*Sigma)));
+ B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))
+ *(2.*s12*(2.*s12*s34 - S2*(S1 + S2) + s34*(S1 - S2))/(Delta*Sigma)));
}
COM G4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2,
mq)*(Delta/s12 + (s12 + S1)/2. -
4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD(
k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD(
k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*2./(s12 + S2));
}
COM E10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta +
12.*mq*mq*S1*(s34 + S1)/Delta/Delta -
4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) +
C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) -
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) +
8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s34*(4.*s12 + 3.*S1 +
S2)/(Delta*Sigma) +
8.*s12*s34*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2,
mq))*(12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*S1*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM F10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (s12*D0DD(k1, k2, q2,
mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta +
12.*mq*mq*s34*(s12 + S1)/Delta/Delta -
4.*s12*pow((s34 + S2),2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*(s34 + S2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) -
C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) +
4.*s12*s34*(S2 - S1)/(Delta*Sigma) +
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) - (B0DD(
k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) +
8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(-12*s34*(2*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s12*s34*s34/(S2*Delta*Delta) +
4.*s34*S1/(Delta*Sigma) -
4.*s34*(s12*s34*(2.*s12 + S2) -
S1*S1*(2.*s12 +
S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*(2.*s34 + S1)/(Delta*Sigma) -
8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) +
s12*(S2 -
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM G10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. +
4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2,
k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*4.*(s34 +
S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) -
B0DD(k2 + q2, mq))*4.*s34/(Delta*S2));
}
COM H1(HLV k1, HLV k2, HLV kh, double mq){
return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq);
}
COM H4(HLV k1, HLV k2, HLV kh, double mq){
return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq);
}
COM H10(HLV k1, HLV k2, HLV kh, double mq){
return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq);
}
COM H2(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H1(k2,k1,kh,mq);
}
COM H5(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H4(k2,k1,kh,mq);
}
COM H12(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H10(k2,k1,kh,mq);
}
// FL and FT functions
COM FL(HLV q1, HLV q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*((2.-
3.*q1.m2()*q2.dot(Q)/detQ2)*(B0DD(q1, mq) -
B0DD(Q, mq)) + (2. -
3.*q2.m2()*q1.dot(Q)/detQ2)*(B0DD(q2, mq) -
B0DD(Q, mq)) - (4.*mq*mq + q1.m2() + q2.m2() +
Q.m2() - 3.*q1.m2()*q2.m2()*Q.m2()/detQ2)*C0DD(
q1, q2, mq) - 2.);
}
COM FT(HLV q1, HLV q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*(Q.m2()*(B0DD(q1, mq) + B0DD(q2, mq) - 2.*B0DD(Q, mq) -
2.*q1.dot(q2)*C0DD(q1, q2, mq)) + (q1.m2() -
q2.m2()) *(B0DD(q1, mq) - B0DD(q2, mq))) -
q1.dot(q2)*FL(q1, q2, mq);
}
HLV ParityFlip(HLV p){
HLV flippedVector;
flippedVector.setE(p.e());
flippedVector.setX(-p.x());
flippedVector.setY(-p.y());
flippedVector.setZ(-p.z());
return flippedVector;
}
/// @brief HC amp for qg->qgH with finite top (i.e. j^{++}_H)
void g_gH_HC(HLV pa, HLV p1,
HLV pH, double mq, current &retAns)
{
current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1,
epsHapart2,conjepsH1part1,conjepsH1part2;
COM ang1a,sqa1;
const double F = 4.*mq*mq/HEJ::vev;
// Easier to have the whole thing as current object so I can use cdot functionality.
// Means I need to write pa,p1 as current objects
to_current(pa, pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
//HEJ gauge
jio(pa,false,p1,false,cura1);
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
} else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(-1./sqrt(2)/ang1a,cura1,conjeps1);
cmult(1./sqrt(2)/sqa1,cura1,epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM h4 = H4(p1,-pa,pH,mq);
const COM h5 = H5(p1,-pa,pH,mq);
const COM h10 = H10(p1,-pa,pH,mq);
const COM h12 = H12(p1,-pa,pH,mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
const COM aH1 = cdot(pHcur, cura1);
current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/ang1a, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/ang1a, epsHa, T2);
}
cmult(sqrt(2.)/ang1a*aH1, epsHa, T3);
cmult(sqrt(2.)/sqa1*aH1, conjepsH1, T4);
cmult(-sqrt(2.)*Fta*pa.dot(p1)*aH1/sqa1, conjeps1, T5);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/ang1a, epsa, T6);
cmult(-aH1/sqrt(2.)/sqa1*h4*8.*COM(0.,1.)*M_PI*M_PI, conjeps1, T7);
cmult(aH1/sqrt(2.)/ang1a*h5*8.*COM(0.,1.)*M_PI*M_PI, epsa, T8);
cmult(aH1*aH1/2./ang1a/sqa1*h10*8.*COM(0.,1.)*M_PI*M_PI, pacur, T9);
cmult(-aH1*aH1/2./ang1a/sqa1*h12*8.*COM(0.,1.)*M_PI*M_PI, p1cur, T10);
current ans;
for(int i=0;i<4;i++)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5[i]+T6[i]+T7[i]+T8[i]+T9[i]+T10[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
/// @brief HNC amp for qg->qgH with finite top (i.e. j^{+-}_H)
void g_gH_HNC(HLV pa, HLV p1, HLV pH, double mq, current &retAns)
{
const double F = 4.*mq*mq/HEJ::vev;
COM ang1a,sqa1;
current conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur,
p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1,
conjepsH1part2;
// Find here if pa, meaning the gluon, is forward or backward
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
jio(pa,true,p1,true,cura1);
joi(p1,true,pa,true,cur1a);
to_current(pa,pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
to_current(pa+p1,paplusp1cur);
to_current(p1-pa,p1minuspacur);
const COM aH1 = cdot(pHcur,cura1);
const COM oneHa = std::conj(aH1); // = cdot(pHcur,cur1a)
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
}
else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(1./sqrt(2)/sqa1, cur1a, epsa);
cmult(-1./sqrt(2)/sqa1, cura1, conjeps1);
const COM phase = cdot(conjeps1, epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM Falpha = FT(p1-pa,pa-p1-pH,mq);
const COM Fbeta = FL(p1-pa,pa-p1-pH,mq);
const COM h1 = H1(p1,-pa, pH, mq);
const COM h2 = H2(p1,-pa, pH, mq);
const COM h4 = H4(p1,-pa, pH, mq);
const COM h5 = H5(p1,-pa, pH, mq);
const COM h10 = H10(p1,-pa, pH, mq);
const COM h12 = H12(p1,-pa, pH, mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
current T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a,
T11b,T12a,T12b,T13;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/sqa1, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())*((p1.x()+COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, epsHa, T2);
}
const COM boxdiagFact = 8.*COM(0.,1.)*M_PI*M_PI;
cmult(aH1*sqrt(2.)/sqa1, epsHa, T3);
cmult(oneHa*sqrt(2.)/sqa1, conjepsH1, T4);
cmult(-2.*phase*Fta*pa.dot(pH), p1cur, T5a);
cmult(2.*phase*Ft1*p1.dot(pH), pacur, T5b);
cmult(-sqrt(2.)*Fta*p1.dot(pa)*oneHa/sqa1, conjeps1, T6);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/sqa1, epsa, T7);
cmult(-boxdiagFact*phase*h2, pacur, T8a);
cmult(boxdiagFact*phase*h1, p1cur, T8b);
cmult(boxdiagFact*aH1/sqrt(2.)/sqa1*h5, epsa, T9);
cmult(-boxdiagFact*oneHa/sqrt(2.)/sqa1*h4, conjeps1, T10);
cmult(boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h10, pacur, T11a);
cmult(-boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h12, p1cur, T11b);
cmult(-phase/(pa-p1).m2()*Falpha*(p1-pa).dot(pa-p1-pH), paplusp1cur, T12a);
cmult(phase/(pa-p1).m2()*Falpha*(pa+p1).dot(pa-p1-pH), p1minuspacur, T12b);
cmult(-phase*Fbeta*(pa-p1-pH).m2(), paplusp1cur, T13);
current ans;
for(int i=0;i<4;i++)
{
- ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]+T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i];
+ ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]
+ +T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
} // namespace anonymous
-// JDC - new amplitude with Higgs emitted close to gluon with full mt effects. Keep usual HEJ-style function call
-double ME_Houtside_gq(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pH, double mq, bool includeBottom, double mq2)
-{
+
+// JDC - new amplitude with Higgs emitted close to gluon with full mt effects.
+// Keep usual HEJ-style function call
+double ME_Houtside_gq(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pH,
+ double mq, bool includeBottom, double mq2
+){
current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip;
current retAns,retAnsb;
joi(p2out,true,p2in,true,cur2bplus);
joi(p2out,false,p2in,false,cur2bminus);
joi(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip);
joi(ParityFlip(p2out),false,ParityFlip(p2in),false,cur2bminusFlip);
COM app1,app2,apm1,apm2;
COM app3, app4, apm3, apm4;
if(!includeBottom)
{
g_gH_HC(p1in,p1out,pH,mq,retAns);
app1=cdot(retAns,cur2bplus);
app2=cdot(retAns,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
app3=cdot(retAns,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,retAns);
apm1=cdot(retAns,cur2bplus);
apm2=cdot(retAns,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
apm3=cdot(retAns,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip);
} else {
g_gH_HC(p1in,p1out,pH,mq,retAns);
g_gH_HC(p1in,p1out,pH,mq2,retAnsb);
app1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
app2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb);
app3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,retAns);
g_gH_HNC(p1in,p1out,pH,mq2,retAnsb);
apm1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
apm2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb);
apm3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
}
return abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1)
+ abs2(apm2) + abs2(apm3) + abs2(apm4);
}
#endif // HEJ_BUILD_WITH_QCDLOOP
double C2gHgm(HLV p2, HLV p1, HLV pH)
{
static double A=1./(3.*M_PI*HEJ::vev);
// Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
double s12,p1p,p2p;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p1p=p1.plus();
p2p=p2.plus();
} else { // opposite case
p1p=p1.minus();
p2p=p2.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=COM(0,1)*A/(2.*s12)*(p2p/p1p*conj(p1perp)*p3perp+p1p/p2p*p1perp*conj(p3perp));
temp=temp*conj(temp);
return temp.real();
}
double C2gHgp(HLV p2, HLV p1, HLV pH)
{
static double A=1./(3.*M_PI*HEJ::vev);
// Implements Eq. (4.23) in hep-ph/0301013
double s12,php,p1p,phm;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
php=pH.plus();
phm=pH.minus();
p1p=p1.plus();
} else { // opposite case
php=pH.minus();
phm=pH.plus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=-COM(0,1)*A/(2.*s12)*( conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p)
+s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm)
-pow(conj(p3perp)
+(1.+php/p1p)*conj(p1perp),2)/((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH)))) );
temp=temp*conj(temp);
return temp.real();
}
double C2qHqm(HLV p2, HLV p1, HLV pH)
{
static double A=1./(3.*M_PI*HEJ::vev);
// Implements Eq. (4.22) in hep-ph/0301013
double s12,p2p,p1p;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p2p=p2.plus();
p1p=p1.plus();
} else { // opposite case
p2p=p2.minus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=A/(2.*s12)*( sqrt(p2p/p1p)*p3perp*conj(p1perp)
+sqrt(p1p/p2p)*p1perp*conj(p3perp) );
temp=temp*conj(temp);
return temp.real();
}
diff --git a/src/Wjets.cc b/src/Wjets.cc
index 8226de6..cde969d 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1213 +1,1217 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Tensor.hh"
#include "HEJ/Constants.hh"
#include <array>
#include <iostream>
using HEJ::Tensor;
using HEJ::init_sigma_index;
using HEJ::metric;
using HEJ::rank3_current;
using HEJ::rank5_current;
using HEJ::eps;
using HEJ::Construct1Tensor;
using HEJ::Helicity;
namespace helicity = HEJ::helicity;
namespace { // Helper Functions
// FKL W Helper Functions
double WProp (const HLV & plbar, const HLV & pl){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - HEJ::MW*HEJ::MW + COM(0.,1.)*HEJ::MW*HEJ::GammaW);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
CCurrent jW (
HLV pout, Helicity helout,
HLV plbar, Helicity hellbar,
HLV pl, Helicity hell,
HLV pin, Helicity helin
){
COM cur[4];
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
CCurrent sum(0.,0.,0.,0.);
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is
// anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hellbar==hell) {
HLV qa=pout+plbar+pl;
HLV qb=pin-plbar-pl;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pl,hell,plbar,hellbar);
CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
COM brac615=t65.dot(vout);
COM brac645=t65.dot(vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
temp2 = joo(pout,helout,pl,helout);
COM prod1665=temp2.dot(t65);
temp3 = joi(plbar,helin,pin,helin);
COM prod5465=temp3.dot(t65);
temp2=joo(pout,helout,plbar,helout);
temp3=joi(pl,hell,pin,helin);
temp5=joi(pout,helout,pin,helin);
CCurrent term1,term2,term3;
term1=(2.*brac615/ta+2.*brac645/tb)*temp5;
term2=(prod1665/ta)*temp3;
term3=(-prod5465/tb)*temp2;
sum=term1+term2+term3;
}
return sum;
}
CCurrent jWbar (
HLV pout, Helicity helout,
HLV plbar, Helicity hellbar,
HLV pl, Helicity hell,
HLV pin, Helicity helin
){
COM cur[4];
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
CCurrent sum(0.,0.,0.,0.);
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is
// anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hellbar==hell) {
HLV qa=pout+plbar+pl;
HLV qb=pin-plbar-pl;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pl,hell,plbar,hellbar);
CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
COM brac615=t65.dot(vout);
COM brac645=t65.dot(vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
temp2 = joo(plbar,helout,pout,helout); // temp2 is <5|alpha|1>
COM prod5165=temp2.dot(t65);
temp3 = jio(pin,helin,pl,helin); // temp3 is <4|alpha|6>
COM prod4665=temp3.dot(t65);
temp2=joo(pl,helout,pout,helout); // temp2 is now <6|mu|1>
temp3=jio(pin,helin,plbar,helin); // temp3 is now <4|mu|5>
temp5=jio(pin,helin,pout,helout); // temp5 is <4|mu|1>
CCurrent term1,term2,term3;
term1 =(-2.*brac615/ta-2.*brac645/tb)*temp5;
term2 =(-prod5165/ta)*temp3;
term3 =(prod4665/tb)*temp2;
sum = term1 + term2 + term3;
}
return sum;
}
// Extremal quark current with W emission.
// Using Tensor class rather than CCurrent
Tensor <1> jW(HLV pin, HLV pout, HLV plbar, HLV pl, bool aqline){
// Build the external quark line W Emmision
Tensor<1> ABCurr = HEJ::current(pl, plbar, helicity::minus);
Tensor<3> J4bBlank;
if (aqline){
J4bBlank = rank3_current(pin,pout,helicity::minus);
}
else{
J4bBlank = rank3_current(pout,pin,helicity::minus);
}
double t4AB = (pout+pl+plbar).m2();
double tbAB = (pin-pl-plbar).m2();
Tensor<2> J4b1 = (J4bBlank.contract(pout+pl+plbar,2))/t4AB;
Tensor<2> J4b2 = (J4bBlank.contract(pin-pl-plbar,2))/tbAB;
Tensor<2> T4bmMom(0.);
if (aqline){
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
T4bmMom(mu, nu) = (J4b1(nu,mu) + J4b2(mu,nu))*COM(0,-1);
}
}
}
else{
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
T4bmMom(nu,mu) = (J4b1(nu,mu) + J4b2(mu,nu))*COM(0,1);
}
}
}
Tensor<1> T4bm = T4bmMom.contract(ABCurr,1);
return T4bm;
}
/**
* @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
){
//@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 tb2 = (pb-p2).m2();
const double tb2g = (pb-p2-pg).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+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 = Construct1Tensor((q1+q2)/taW1 + (pb/pb.dot(pg)
+p2/p2.dot(pg)) * tb2/(2*tb2g));
Tensor<3> J31a = rank3_current(p1, pa, h1);
Tensor<2> J2_qaW =J31a.contract((pa-pW)/taW, 2);
Tensor<2> J2_p1W =J31a.contract((p1+pW)/s1W, 2);
Tensor<3> L1a = outer(Tq1q2, J2_qaW);
Tensor<3> L1b = outer(Tq1q2, J2_p1W);
Tensor<3> L2a = outer(-pg-q1,J2_qaW)/taW1;
Tensor<3> L2b = outer(-pg-q1, J2_p1W)/taW1;
Tensor<3> L3 = outer(metric(), J2_qaW.contract(pg-q2,1)+J2_p1W.contract(pg-q2,2))/taW1;
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();
double WPropfact = WProp(plbar, pl);
//Divide by WProp
amp*=WPropfact;
//Divide by t-channels
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
){
CCurrent mj1m,mj2p,mj2m;
HLV q1=p1in-p1out-plbar-pl;
HLV q2=-(p2in-p2out);
- if(aqlineb) mj1m=jWbar(p1out,helicity::minus,plbar,helicity::minus,pl,helicity::minus,p1in,helicity::minus);
- else mj1m=jW(p1out,helicity::minus,plbar,helicity::minus,pl,helicity::minus,p1in,helicity::minus);
+ if(aqlineb) mj1m=jWbar(p1out,helicity::minus, plbar,helicity::minus,
+ pl,helicity::minus, p1in,helicity::minus);
+ else mj1m=jW(p1out,helicity::minus, plbar,helicity::minus,
+ pl,helicity::minus, p1in,helicity::minus);
mj2p=joi(p2out,!aqlinef, p2in,!aqlinef);
mj2m=joi(p2out, aqlinef, p2in, aqlinef);
COM Mmp=mj1m.dot(mj2p);
COM Mmm=mj1m.dot(mj2m);
// sum of spinor strings ||^2
double a2Mmp=abs2(Mmp);
double a2Mmm=abs2(Mmm);
double WPropfact = WProp(plbar, pl);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(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){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false);
}
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true);
}
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false);
}
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true);
}
double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false)*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false)*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){
CCurrent mj1m,mj2p,mj2m, jgbm,jgbp,j2gm,j2gp;
HLV q1=p1in-p1out-plbar-pl;
HLV q2=-(p2in-p2out-pg);
HLV q3=-(p2in-p2out);
- if(aqlineb) mj1m=jWbar(p1out,helicity::minus,plbar,helicity::minus,pl,helicity::minus,p1in,helicity::minus);
- else mj1m=jW(p1out,helicity::minus,plbar,helicity::minus,pl,helicity::minus,p1in,helicity::minus);
+ if(aqlineb) mj1m=jWbar(p1out,helicity::minus, plbar,helicity::minus,
+ pl,helicity::minus, p1in,helicity::minus);
+ else mj1m=jW(p1out,helicity::minus, plbar,helicity::minus,
+ pl,helicity::minus, p1in,helicity::minus);
// Note <p1|eps|pa> current split into two by gauge choice.
// See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
mj2p=joi(p2out,!aqlinef,p2in,!aqlinef);
mj2m=joi(p2out,aqlinef,p2in,aqlinef);
jgbp=joi(pg,!aqlinef,p2in,!aqlinef);
jgbm=joi(pg,aqlinef,p2in,aqlinef);
// Note for function joo(): <p1+|pg+> = <pg-|p1->.
j2gp=joo(p2out,!aqlinef,pg,!aqlinef);
j2gm=joo(p2out,aqlinef,pg,aqlinef);
// Dot products of these which occur again and again
COM MWmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent qsum(q2+q3);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
CCurrent p2o(p2out);
CCurrent p2i(p2in);
Lmm=( (-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m + 2.*mj2m.dot(pg)*mj1m
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmm/2. ) )/q3.m2();
Lmp=( (-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p + 2.*mj2p.dot(pg)*mj1m
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmp/2. ) )/q3.m2();
U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
U2mp=((-1.)*j2gp.dot(mj1m)*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.*vabs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
ampsq*=WProp(plbar, pl);
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
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true);
}
double ME_W_unof_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, false, false);
}
double ME_W_unof_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, false, true);
}
double ME_W_unof_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, true, false);
}
double ME_W_unof_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, true, true);
}
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
){
//Calculate different Helicity choices
const Helicity h = aqlineb?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::plus,helicity::plus);
double ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::plus,helicity::minus);
double ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::minus,helicity::plus);
double ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::minus,helicity::minus);
//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
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false)*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
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true)*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){
//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);
double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::plus,helicity::minus);
double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::minus,helicity::plus);
double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::minus,helicity::minus);
//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){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false);
}
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true);
}
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false)*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){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true)*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
){
init_sigma_index();
// 2 independent helicity choices (complex conjugation related).
Tensor<1> TMmmm1 = qggm1(pb,p2,p3,helicity::minus,helicity::minus, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,helicity::minus,helicity::minus, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,helicity::minus,helicity::minus, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,helicity::minus,helicity::plus, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,helicity::minus,helicity::plus, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,helicity::minus,helicity::plus, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa);
//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
double WPropfact = WProp(plbar, pl);
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
){
init_sigma_index();
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, 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);
// 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
double WPropfact = WProp(plbar, pl);
amp*=WPropfact;
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
){
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, 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);}
Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, helicity::minus, nabove);
Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, helicity::minus, nabove);
Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, helicity::minus, nabove);
Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, helicity::plus, nabove);
Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, helicity::plus, nabove);
Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, helicity::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
double WPropfact = WProp(plbar, pl);
amp*=WPropfact;
return amp;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 3:14 PM (1 d, 6 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023250
Default Alt Text
(110 KB)

Event Timeline