diff --git a/src/Wjets.cc b/src/Wjets.cc
index cacc376..63bc609 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,799 +1,800 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/jets.hh"
// generated headers
#include "HEJ/currents/jV_j.hh"
#include "HEJ/currents/jV_juno.hh"
#include "HEJ/currents/jVuno_j.hh"
#include "HEJ/currents/jW_jqqbar.hh"
#include "HEJ/currents/jW_qqbar_j.hh"
#include "HEJ/currents/jWqqbar_j.hh"
#include "HEJ/currents/j_Wqqbar_j.hh"
namespace HEJ {
namespace currents {
namespace {
// short hands
- using std::abs;
using std::conj;
- using std::pow;
// --- Helper Functions ---
// FKL W Helper Functions
double WProp(const HLV & plbar, const HLV & pl,
ParticleProperties const & wprop
){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
+ COM(0.,1.)*wprop.mass*wprop.width);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
/**
* @brief Contraction of W + unordered jet current with FKL current
*
* See eq:wunocurrent in the developer manual for the definition
* of the W + unordered jet current
*/
template<Helicity h1, Helicity h2, Helicity pol>
double jM2Wuno(
HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
HLV const & pa, HLV const & p2, HLV const & pb,
ParticleProperties const & wprop
){
using helicity::minus;
const COM u1 = U1<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM u2 = U2<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM l = L<h1, minus, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM x = u1 - l;
const COM y = u2 + l;
double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
amp *= WProp(plbar, pl, wprop);
const HLV q1g = pa-pl-plbar-p1-pg;
const HLV q2 = p2-pb;
const double t1 = q1g.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
/**
* @brief Contraction of W + extremal qqbar current with FKL current
*
* See eq:crossed in the developer manual for the definition
* of the W + extremal qqbar current.
*
*/
template<Helicity h1, Helicity h2, Helicity hg>
double jM2_W_extremal_qqbar(
HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
HLV const & pqbar, HLV const & p3, HLV const & pb,
ParticleProperties const & wprop
){
const COM u1Xcontr = U1X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM u2Xcontr = U2X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM lXcontr = LX<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM x = u1Xcontr - lXcontr;
const COM y = u2Xcontr + lXcontr;
double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
amp *= WProp(plbar, pl, wprop);
const HLV q1 = pg-pl-plbar-pq-pqbar;
const HLV q2 = p3-pb;
const double t1 = q1.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
//! 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 const & p1out, HLV plbar, HLV pl,
HLV const & p1in, HLV const & p2out, HLV const & p2in,
bool aqlineb, bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out);
const double wPropfact = WProp(plbar, pl, wprop);
// we assume that the W is emitted off a quark line
// if this is not the case, we have to apply CP conjugation,
// which is equivalent to swapping lepton and antilepton momenta
if(aqlineb) std::swap(pl, plbar);
const COM ampm = jV_j<minus, minus, minus>(p1in,p1out,p2in,p2out,pl,plbar);
const COM ampp = jV_j<minus, minus, plus>(p1in,p1out,p2in,p2out,pl,plbar);
const double Msqr = std::norm(ampm) + std::norm(ampp);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return C_F*C_F*wPropfact*Msqr/(q1.m2()*q2.m2()*(N_C*N_C - 1)*4);
}
} // Anonymous Namespace
double ME_W_qQ(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop);
}
double ME_W_qQbar(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true, wprop);
}
double ME_W_qbarQ(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop);
}
double ME_W_qbarQbar(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true, wprop);
}
double ME_W_qg(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop)
*K_g(p2out, p2in)/C_F;
}
double ME_W_qbarg(
HLV const & p1out,
HLV const & plbar, HLV const & pl,
HLV const & p1in,
HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop)
*K_g(p2out, p2in)/C_F;
}
namespace {
// helicity amplitude squares for contraction of W current with unordered
// current
template<Helicity h2, Helicity hg>
double ampsq_jW_juno(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p2,
HLV const & pg,
HLV const & pl, HLV const & plbar
){
using helicity::minus;
const COM u1 = U1_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
const COM u2 = U2_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
const COM l = L_jV<minus, minus, h2, hg>(pa,p1,pb,p2,pg,pl,plbar);
return 2.*C_F*std::real((l-u1)*std::conj(l+u2))
+ 2.*C_F*C_F/3.*std::norm(u1+u2)
;
}
/**
* @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.
*
* @TODO use const & for pl plbar
*/
double jW_juno(
HLV const & p1out, HLV plbar, HLV pl, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pg, bool aqlineb,
bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
// we assume that the W is emitted off a quark line
// if this is not the case, we have to apply CP conjugation,
// which is equivalent to swapping lepton and antilepton momenta
if(aqlineb) std::swap(pl, plbar);
// helicity assignments assume aqlinef == false
// in the end, this is irrelevant since we sum over all helicities
const double ampsq =
+ ampsq_jW_juno<minus, minus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
+ ampsq_jW_juno<minus, plus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
+ ampsq_jW_juno<plus, minus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
+ ampsq_jW_juno<plus, plus>(p1in,p1out,p2in,p2out,pg,pl,plbar)
;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out-pg);
//! @TODO explain magic 16
return WProp(plbar, pl, wprop)*ampsq/(16.*q2.m2()*q1.m2());
}
} // Anonymous Namespace
double ME_W_unob_qQ(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl, ParticleProperties
const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
}
double ME_W_unob_qQbar(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true, wprop);
}
double ME_W_unob_qbarQ(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false, wprop);
}
double ME_W_unob_qbarQbar(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true, wprop);
}
namespace {
// helicity sum helper for jWuno_j(...)
template<Helicity h1>
double jWuno_j_helsum(
HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
HLV const & pa, HLV const & p2, HLV const & pb,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const double ME2h1pp = jM2Wuno<h1, plus, plus>(
pg, p1, plbar, pl, pa, p2, pb, wprop
);
const double ME2h1pm = jM2Wuno<h1, plus, minus>(
pg, p1, plbar, pl, pa, p2, pb, wprop
);
const double ME2h1mp = jM2Wuno<h1, minus, plus>(
pg, p1, plbar, pl, pa, p2, pb, wprop
);
const double ME2h1mm = jM2Wuno<h1, minus, minus>(
pg, p1, plbar, pl, pa, p2, pb, wprop
);
return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm;
}
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming
* @types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same
* side. Handles all possible incoming states. Note this handles both forward
* and back- -ward Wuno emission. For forward, ensure p1out is the uno and W
* emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean
up * ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(
HLV const & pg, HLV const & p1out, HLV const & plbar, HLV const & pl,
HLV const & p1in, HLV const & p2out, HLV const & p2in, bool aqlineb,
ParticleProperties const & wprop
){
const auto helsum = aqlineb?
jWuno_j_helsum<helicity::plus>:
jWuno_j_helsum<helicity::minus>;
//Helicity sum and average over initial states
return helsum(pg, p1out,plbar,pl,p1in,p2out,p2in, wprop)/(4.*C_A*C_A);
}
} // Anonymous Namespace
double ME_Wuno_qQ(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qQbar(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qbarQ(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qbarQbar(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qg(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop)
*K_g(p2out, p2in)/C_F;
}
double ME_Wuno_qbarg(
HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in,
HLV const & pg, HLV const & plbar, HLV const & pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop)
*K_g(p2out, p2in)/C_F;
}
// helicity sum helper for jWqqx_j(...)
template<Helicity h1>
double jWqqx_j_helsum(
HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
HLV const & pqbar, HLV const & p3, HLV const & pb,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const double amp_h1pp = jM2_W_extremal_qqbar<h1, plus, plus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double amp_h1pm = jM2_W_extremal_qqbar<h1, plus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double amp_h1mp = jM2_W_extremal_qqbar<h1, minus, plus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double amp_h1mm = jM2_W_extremal_qqbar<h1, minus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
return amp_h1pp + amp_h1pm + amp_h1mp + amp_h1mm;
}
/**
* @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 const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in, bool aqlinef,
ParticleProperties const & wprop
){
const auto helsum = aqlinef?
jWqqx_j_helsum<helicity::plus>:
jWqqx_j_helsum<helicity::minus>;
//Helicity sum and average over initial states.
double ME2 = helsum(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, wprop)/
(4.*C_A*C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(
HLV const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
}
double ME_WExqqx_qqbarQ(
HLV const & pgin, HLV const & pqbarout, HLV const & plbar, HLV const & pl,
HLV const & pqout, HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
}
double ME_WExqqx_qbarqg(
HLV const & pgin, HLV const & pqout, HLV const & plbar, HLV const & pl,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
*K_g(p2out,p2in)/C_F;
}
double ME_WExqqx_qqbarg(
HLV const & pgin, HLV const & pqbarout, HLV const & plbar, HLV const & pl,
HLV const & pqout, HLV const & p2out, HLV const & p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
*K_g(p2out,p2in)/C_F;
}
namespace {
template<Helicity h1, Helicity hg>
double jW_jqqbar(
HLV const & pb,
HLV const & p2,
HLV const & p3,
HLV const & pa,
HLV const & p1,
HLV const & pl,
HLV const & plbar
){
+ using std::norm;
- static constexpr COM cm1m1 = 8./3.;
- static constexpr COM cm2m2 = 8./3.;
- static constexpr COM cm3m3 = 6.;
- static constexpr COM cm1m2 =-1./3.;
+ static constexpr double cm1m1 = 8./3.;
+ static constexpr double cm2m2 = 8./3.;
+ static constexpr double cm3m3 = 6.;
+ static constexpr double cm1m2 =-1./3.;
static constexpr COM cm1m3 = COM{0.,-3.};
static constexpr COM cm2m3 = COM{0.,3.};
const COM m1 = jW_qggm1<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
const COM m2 = jW_qggm2<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
const COM m3 = jW_qggm3<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
- return real(
- cm1m1*pow(abs(m1),2) + cm2m2*pow(abs(m2),2)
- + cm3m3*pow(abs(m3),2) + 2.*real(cm1m2*m1*conj(m2))
- )
+ return
+ + cm1m1*norm(m1)
+ + cm2m2*norm(m2)
+ + cm3m3*norm(m3)
+ + 2.*real(cm1m2*m1*conj(m2))
+ 2.*real(cm1m3*m1*conj(m3))
+ 2.*real(cm2m3*m2*conj(m3)) ;
}
} // Anonymous Namespace
// contraction of W current with extremal qqbar on the other side
double ME_W_Exqqx_QQq(
HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2,
HLV const & p3, HLV const & plbar, HLV const & pl, bool aqlinepa,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const double wPropfact = WProp(plbar, pl, wprop);
const double prefactor = 2.*wPropfact/24./4./(
(pa-p1-pl-plbar).m2()*(p2+p3-pb).m2()
);
if(aqlinepa) {
return prefactor*(
jW_jqqbar<plus, minus>(pb,p2,p3,pa,p1,pl,plbar)
+ jW_jqqbar<plus, plus>(pb,p2,p3,pa,p1,pl,plbar)
);
}
return prefactor*(
jW_jqqbar<minus, minus>(pb,p2,p3,pa,p1,pl,plbar)
+ jW_jqqbar<minus, plus>(pb,p2,p3,pa,p1,pl,plbar)
);
}
namespace {
// helper function for matrix element for W + Jets with central qqbar
template<Helicity h1a, Helicity h4b>
double amp_WCenqqx_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p4,
HLV const & pq, HLV const & pqbar,
HLV const & pl, HLV const & plbar,
HLV const & q11, HLV const & q12
){
+ using std::norm;
+
const COM sym = M_sym_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM cross = M_cross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM uncross = M_uncross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
// Colour factors
- static constexpr COM cmsms = 3.;
- static constexpr COM cmumu = 4./3.;
- static constexpr COM cmcmc = 4./3.;
+ static constexpr double cmsms = 3.;
+ static constexpr double cmumu = 4./3.;
+ static constexpr double cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0., 3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
- static constexpr COM cmumc = -1./6.;
+ static constexpr double cmumc = -1./6.;
- return real(
- cmsms*pow(abs(sym),2)
- +cmumu*pow(abs(uncross),2)
- +cmcmc*pow(abs(cross),2)
- )
+ return
+ +cmsms*norm(sym)
+ +cmumu*norm(uncross)
+ +cmcmc*norm(cross)
+2.*real(cmsmu*sym*conj(uncross))
+2.*real(cmsmc*sym*conj(cross))
+2.*real(cmumc*uncross*conj(cross))
;
}
} // Anonymous Namespace
// matrix element for W + Jets with W emission off central qqbar
double ME_WCenqqx_qq(
HLV const & pa, HLV const & pb, HLV const & pl, HLV const & plbar,
std::vector<HLV> const & partons, bool /* aqlinepa */, bool /* aqlinepb */, bool
qqxmarker, int nabove, ParticleProperties const & wprop
){
using helicity::plus;
using helicity::minus;
CLHEP::HepLorentzVector q1 = pa;
for(int i = 0; i <= nabove; ++i){
q1 -= partons[i];
}
const auto qq = split_into_lightlike(q1);
// since q1.m2() < 0 the following assertion is always true
// see documentation for split_into_lightlike
assert(qq.second.e() < 0);
HLV pq;
HLV pqbar;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];
}
const HLV p1 = partons.front();
const HLV p4 = partons.back();
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
// the first helicity label is for aqlinepa == true,
// the second one for aqlinepb == true
// In principle the corresponding helicity should be flipped
// if either of them is false, but we sum up all helicities,
// so this has no effect in the end.
const double amp_mm = amp_WCenqqx_qq<minus, minus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_mp = amp_WCenqqx_qq<minus, plus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pm = amp_WCenqqx_qq<plus, minus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pp = amp_WCenqqx_qq<plus, plus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double t1 = q1.m2();
const double t3 = (q1-pl-plbar-pq-pqbar).m2();
const double amp = WProp(plbar, pl, wprop)*(
amp_mm+amp_mp+amp_pm+amp_pp
)/(9.*4.*t1*t1*t3*t3);
return amp;
}
// helper function for matrix element for W + Jets with central qqbar
// W emitted off extremal parton
// @TODO: code duplication with amp_WCenqqx_qq
template<Helicity h1a, Helicity hqqbar>
double amp_W_Cenqqx_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & pn,
HLV const & pq, HLV const & pqbar,
HLV const & pl, HLV const & plbar,
HLV const & q11, HLV const & q12
){
+ using std::norm;
const COM crossed = M_cross<h1a, hqqbar>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
);
const COM uncrossed = M_qbar<h1a, hqqbar>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
);
const COM sym = M_sym<h1a, hqqbar>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
);
//Colour factors:
- static constexpr COM cmsms = 3.;
- static constexpr COM cmumu = 4./3.;
- static constexpr COM cmcmc = 4./3.;
+ static constexpr double cmsms = 3.;
+ static constexpr double cmumu = 4./3.;
+ static constexpr double cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0.,3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
- static constexpr COM cmumc = -1./6.;
+ static constexpr double cmumc = -1./6.;
- return real(
- cmsms*pow(abs(sym),2)
- +cmumu*pow(abs(uncrossed),2)
- +cmcmc*pow(abs(crossed),2)
- )
+ return
+ +cmsms*norm(sym)
+ +cmumu*norm(uncrossed)
+ +cmcmc*norm(crossed)
+2.*real(cmsmu*sym*conj(uncrossed))
+2.*real(cmsmc*sym*conj(crossed))
+2.*real(cmumc*uncrossed*conj(crossed))
;
}
// matrix element for W + Jets with W emission *not* off central qqbar
double ME_W_Cenqqx_qq(
HLV pa, HLV pb, HLV pl,HLV plbar,
std::vector<HLV> partons, bool aqlinepa, bool aqlinepb,
bool qqxmarker, int nabove, int nbelow, bool forwards,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
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 q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
const auto qq = split_into_lightlike(q1);
HLV pq;
HLV pqbar;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
// we assume that the W is emitted off a quark line
// if this is not the case, we have to apply CP conjugation,
// which is equivalent to swapping lepton and antilepton momenta
if(aqlinepb) std::swap(pl, plbar);
const HLV p1 = partons.front();
const HLV pn = partons.back();
// helicity labels are for aqlinepa == aqlineb == false,
// In principle the helicities should be flipped
// if either of them is true, but we sum up all helicities,
// so this has no effect in the end.
const double amp_mm = amp_W_Cenqqx_qq<minus, minus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_mp = amp_W_Cenqqx_qq<minus, plus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pm = amp_W_Cenqqx_qq<plus, minus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pp = amp_W_Cenqqx_qq<plus, plus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double t1 = q1.m2();
const double t3 = (q1 - pq - pqbar).m2();
const double amp= WProp(plbar, pl, wprop)*(
amp_mm+amp_mp+amp_pm+amp_pp
)/(9.*4.*t1*t1*t3*t3);
return amp;
}
} // namespace currents
} // namespace HEJ
diff --git a/src/jets.cc b/src/jets.cc
index ca7c962..b61a436 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,592 +1,592 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include "HEJ/Constants.hh"
// generated headers
#include "HEJ/currents/j_j.hh"
#include "HEJ/currents/j_qqbar_j.hh"
#include "HEJ/currents/j_jqqbar.hh"
using HEJ::Helicity;
namespace helicity = HEJ::helicity;
namespace {
// short hand for math functions
using std::abs;
using std::conj;
- using std::pow;
using std::sqrt;
} // Anonymous Namespace
namespace HEJ {
namespace currents {
// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
// @TODO: this is not a current and should be moved somewhere else
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
HLV const & pout,
HLV const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
CCurrent CCurrent::operator+(const CCurrent& other) const
{
COM result_c0=c0 + other.c0;
COM result_c1=c1 + other.c1;
COM result_c2=c2 + other.c2;
COM result_c3=c3 + other.c3;
return {result_c0,result_c1,result_c2,result_c3};
}
CCurrent CCurrent::operator-(const CCurrent& other) const
{
COM result_c0=c0 - other.c0;
COM result_c1=c1 - other.c1;
COM result_c2=c2 - other.c2;
COM result_c3=c3 - other.c3;
return {result_c0,result_c1,result_c2,result_c3};
}
CCurrent CCurrent::operator*(const double x) const
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return {result_c0,result_c1,result_c2,result_c3};
}
CCurrent CCurrent::operator/(const double x) const
{
COM result_c0=CCurrent::c0/x;
COM result_c1=CCurrent::c1/x;
COM result_c2=CCurrent::c2/x;
COM result_c3=CCurrent::c3/x;
return {result_c0,result_c1,result_c2,result_c3};
}
CCurrent CCurrent::operator*(const COM x) const
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return {result_c0,result_c1,result_c2,result_c3};
}
CCurrent CCurrent::operator/(const COM x) const
{
COM result_c0=(CCurrent::c0)/x;
COM result_c1=(CCurrent::c1)/x;
COM result_c2=(CCurrent::c2)/x;
COM result_c3=(CCurrent::c3)/x;
return {result_c0,result_c1,result_c2,result_c3};
}
std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
{
os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
return os;
}
CCurrent operator * ( double x, CCurrent const & m)
{
return m*x;
}
CCurrent operator * ( COM const & x, CCurrent const & m)
{
return m*x;
}
CCurrent operator / ( double x, CCurrent const & m)
{
return m/x;
}
CCurrent operator / ( COM const & x, CCurrent const & m)
{
return m/x;
}
COM CCurrent::dot(HLV const & p1) const
{
// Current goes (E,px,py,pz)
// Vector goes (px,py,pz,E)
return p1[3]*c0-p1[0]*c1-p1[1]*c2-p1[2]*c3;
}
COM CCurrent::dot(CCurrent const & p1) const
{
return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3;
}
//Current Functions
void joi(HLV const & pout, bool helout,
HLV const & pin, bool helin, current &cur
){
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
const double sqpop = sqrt(abs(pout.plus()));
const double sqpom = sqrt(abs(pout.minus()));
// Allow for "jii" format
const COM poperp = (pout.x()==0 && pout.y() ==0) ? -1 : pout.x()+COM(0,1)*pout.y();
if (helout != helin) {
throw std::invalid_argument{"Non-matching helicities"};
}
if (!helout) { // negative helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(abs(pin.plus()));
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * poperp / abs(poperp);
cur[2] = -COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(abs(pin.minus()));
cur[0] = -sqpom * sqpim * poperp / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = COM(0,1) * cur[1];
cur[3] = -cur[0];
}
} else { // positive helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(abs(pin.plus()));
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
cur[2] = COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
double sqpim = sqrt(abs(pin.minus()));
cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = -COM(0,1) * cur[1];
cur[3] = -cur[0];
}
}
}
CCurrent joi(HLV const & pout, bool helout, HLV const & pin, bool helin)
{
current cur;
joi(pout, helout, pin, helin, cur);
return {cur[0],cur[1],cur[2],cur[3]};
}
void jio(HLV const & pin, bool helin,
HLV const & pout, bool helout, current &cur
) {
joi(pout, !helout, pin, !helin, cur);
}
CCurrent jio(HLV const & pin, bool helin, HLV const & pout, bool helout){
current cur;
jio(pin, helin, pout, helout, cur);
return {cur[0],cur[1],cur[2],cur[3]};
}
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {
// Zero our current
cur[0] = 0.0;
cur[1] = 0.0;
cur[2] = 0.0;
cur[3] = 0.0;
if (heli!=helj) {
throw std::invalid_argument{"Non-matching helicities"};
}
if ( heli ) { // If positive helicity swap momenta
std::swap(pi,pj);
}
const double sqpjp = sqrt(abs(pj.plus() ));
const double sqpjm = sqrt(abs(pj.minus()));
const double sqpip = sqrt(abs(pi.plus() ));
const double sqpim = sqrt(abs(pi.minus()));
// Allow for "jii" format
const COM piperp = (pi.x()==0 && pi.y() ==0) ? -1 : pi.x()+COM(0,1)*pi.y();
const COM pjperp = (pj.x()==0 && pj.y() ==0) ? -1 : pj.x()+COM(0,1)*pj.y();
const COM phasei = piperp / abs(piperp);
const COM phasej = pjperp / abs(pjperp);
cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
}
CCurrent joo(HLV const & pi, bool heli, HLV const & pj, bool helj)
{
current cur;
joo(pi, heli, pj, helj, cur);
return {cur[0],cur[1],cur[2],cur[3]};
}
namespace {
//@{
/**
* @brief Pure Jet FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1.
* @param p1in Incoming Particle 1.
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
*
* Calculates j_\mu j^\mu.
* Handles all possible incoming states. Helicity doesn't matter since we sum
* over all of them. In addition, we only have to distinguish between the two
* possibilities of contracting same-helicity or opposite-helicity currents.
*/
double j_j(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
COM const Mp = HEJ::j_j<helicity::plus>(p1in, p1out, p2in, p2out);
COM const Mm = HEJ::j_j<helicity::minus>(p1in, p1out, p2in, p2out);
double const sst=abs2(Mm)+abs2(Mp);
HLV const q1=p1in-p1out;
HLV const q2=-(p2in-p2out);
// Multiply by Cf^2 (colour) * 2 (helicities)
return 2.*C_F*C_F*(sst)/(q1.m2()*q2.m2());
}
} // Anonymous Namespace
double ME_qQ(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qbarQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/C_F;
}
double ME_qbarg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/(C_F);
}
double ME_gg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p1out, p1in)*K_g(p2out, p2in)/(C_F*C_F);
}
//@}
namespace {
double juno_j(HLV const & pg, HLV const & p1out,
HLV const & p1in, HLV const & p2out, HLV const & p2in
){
// 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
// Note <p1|eps|pa> current split into two by gauge choice.
// See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
CCurrent mj1p=joi(p1out, false, p1in, false);
CCurrent mj1m=joi(p1out, true, p1in, true);
CCurrent jgap=joi(pg, false, p1in, false);
CCurrent jgam=joi(pg, true, p1in, true);
// Note for function joo(): <p1+|pg+> = <pg-|p1->.
CCurrent j2gp=joo(p1out, false, pg, false);
CCurrent j2gm=joo(p1out, true, pg, true);
CCurrent mj2p=joi(p2out, false, p2in, false);
CCurrent mj2m=joi(p2out, true, p2in, true);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p);
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
CCurrent p1o(p1out);
CCurrent p2o(p2out);
CCurrent p2i(p2in);
CCurrent qsum(q1+qg);
CCurrent p1i(p1in);
CCurrent Lmm=(qsum*(Mmm)+(-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
CCurrent Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
CCurrent Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
CCurrent Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
CCurrent U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
CCurrent U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
CCurrent U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
CCurrent U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
CCurrent U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
CCurrent U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
CCurrent U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
CCurrent U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
constexpr double cf=C_F;
double amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
double amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
//Divide by t-channels
ampsq/=q2.m2()*qg.m2();
ampsq/=16.;
// Factor of (Cf/Ca) for each quark to match j_j.
ampsq*=(C_F*C_F)/(C_A*C_A);
return ampsq;
}
} // Anonymous Namespace
//Unordered bits for pure jet
double ME_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qg( HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F;
}
double ME_unob_qbarg(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F;
}
namespace {
// helicity amplitudes for j jqqbar contraction
template<Helicity h1, Helicity hg>
double amp_j_jqqbar(
HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2, HLV const & p3
) {
// TODO: code duplication with Wjets.cc
using std::norm;
static constexpr double cm1m1 = 8./3.;
static constexpr double cm2m2 = 8./3.;
static constexpr double cm3m3 = 6.;
static constexpr double cm1m2 =-1./3.;
static constexpr COM cm1m3 = COM{0.,-3.};
static constexpr COM cm2m3 = COM{0.,3.};
const COM m1 = j_qggm1<h1,hg>(pb,p2,p3,pa,p1);
const COM m2 = j_qggm2<h1,hg>(pb,p2,p3,pa,p1);
const COM m3 = j_qggm3<h1,hg>(pb,p2,p3,pa,p1);
return
+ cm1m1*norm(m1)
+ cm2m2*norm(m2)
+ cm3m3*norm(m3)
+ 2.*real(cm1m2*m1*conj(m2))
+ 2.*real(cm1m3*m1*conj(m3))
+ 2.*real(cm2m3*m2*conj(m3))
;
}
//Now the function to give helicity/colour sum/average
double MqgtqQQ(HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2, HLV const & p3
){
using helicity::minus;
using helicity::plus;
const double Mmmm = amp_j_jqqbar<minus, minus>(pa, pb, p1, p2, p3);
const double Mmmp = amp_j_jqqbar<minus, plus>(pa, pb, p1, p2, p3);
const double Mpmm = amp_j_jqqbar< plus, minus>(pa, pb, p1, p2, p3);
const double Mpmp = amp_j_jqqbar< plus, plus>(pa, pb, p1, p2, p3);
// Factor of 2 for the helicity for conjugate configurations
return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/3.)/(pa-p1).m2()/(p2+p3-pb).m2();
}
} // Anonymous Namespace
// Extremal qqx
double ME_Exqqx_qbarqQ(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout);
}
double ME_Exqqx_qqbarQ(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout);
}
double ME_Exqqx_qbarqg(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout)*K_g(p2out,p2in)/C_F;
}
double ME_Exqqx_qqbarg(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout)*K_g(p2out,p2in)/C_F;
}
namespace {
// helicity amplitudes for matrix element with central qqbar
template<Helicity h1a, Helicity h4b>
double amp_Cenqqx_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p4,
HLV const & pq, HLV const & pqbar,
HLV const & q11, HLV const & q12
){
+ using std::norm;
+
const COM sym = M_sym<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, q11, q12
);
const COM cross = M_cross<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, q11, q12
);
const COM uncross = M_qbar<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, q11, q12
);
// Colour factors
- static constexpr COM cmsms = 3.;
- static constexpr COM cmumu = 4./3.;
- static constexpr COM cmcmc = 4./3.;
+ static constexpr double cmsms = 3.;
+ static constexpr double cmumu = 4./3.;
+ static constexpr double cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0., 3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
- static constexpr COM cmumc = -1./6.;
+ static constexpr double cmumc = -1./6.;
- return real(
- cmsms*pow(abs(sym),2)
- +cmumu*pow(abs(uncross),2)
- +cmcmc*pow(abs(cross),2)
- )
+ return
+ +cmsms*norm(sym)
+ +cmumu*norm(uncross)
+ +cmcmc*norm(cross)
+2.*real(cmsmu*sym*conj(uncross))
+2.*real(cmsmc*sym*conj(cross))
+2.*real(cmumc*uncross*conj(cross))
;
}
} // Anonymous Namespace
double ME_Cenqqx_qq(
HLV const & pa, HLV const & pb,
std::vector<HLV> const & partons,
bool /* aqlinepa */, bool /* aqlinepb */,
const bool qqxmarker, const std::size_t nabove
){
using helicity::plus;
using helicity::minus;
CLHEP::HepLorentzVector q1 = pa;
for(std::size_t i = 0; i <= nabove; ++i){
q1 -= partons[i];
}
const auto qq = split_into_lightlike(q1);
// since q1.m2() < 0 the following assertion is always true
// see documentation for split_into_lightlike
assert(qq.second.e() < 0);
HLV pq = partons[nabove+1];
HLV pqbar = partons[nabove+2];
if(qqxmarker) std::swap(pq, pqbar);
const HLV p1 = partons.front();
const HLV pn = partons.back();
// 8 helicity choices, but only 4 independent ones
//(complex conjugation related).
// In principle, the proper helicity labeling depends on
// whether we have antiquark lines (aqlinea and aqlineb).
// However, this only corresponds to a relabeling,
// which doesn't change the sum over all helicities below.
const double amp_mm = amp_Cenqqx_qq<minus, minus>(
pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
);
const double amp_mp = amp_Cenqqx_qq<minus, plus>(
pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
);
const double amp_pm = amp_Cenqqx_qq<plus, minus>(
pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
);
const double amp_pp = amp_Cenqqx_qq<plus, plus>(
pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second
);
//Result (averaged, without coupling or t-channel props). Factor of
//2 for the 4 helicity configurations I didn't work out explicitly
const HLV q3 = q1 - pq - pqbar;
return (2.*(amp_mm+amp_mp+amp_pm+amp_pp)/9./4.) /
((pa-p1).m2()*(pb-pn).m2()*q1.m2()*q3.m2());
}
} // namespace currents
} // namespace HEJ