Page MenuHomeHEPForge
Files F7877316

No OneTemporary
Actions

View Options

diff --git a/src/Wjets.cc b/src/Wjets.cc
index 6abe7dd..5458ed6 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,796 +1,796 @@
/**
* \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/jets.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/exceptions.hh"
// generated headers
#include "HEJ/currents/jV_j.hh"
#include "HEJ/currents/jVuno_j.hh"
#include "HEJ/currents/jWqqbar_j.hh"
#include "HEJ/currents/jW_qqbar_j.hh"
#include "HEJ/currents/j_Wqqbar_j.hh"
#include "HEJ/currents/jW_jqqbar.hh"
#include "HEJ/currents/jV_juno.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);
+ 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);
+ 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);
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 ME2h1pp = jM2_W_extremal_qqbar<h1, plus, plus>(
+ const double amp_h1pp = jM2_W_extremal_qqbar<h1, plus, plus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
- const double ME2h1pm = jM2_W_extremal_qqbar<h1, plus, minus>(
+ const double amp_h1pm = jM2_W_extremal_qqbar<h1, plus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
- const double ME2h1mp = jM2_W_extremal_qqbar<h1, minus, plus>(
+ const double amp_h1mp = jM2_W_extremal_qqbar<h1, minus, plus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
- const double ME2h1mm = jM2_W_extremal_qqbar<h1, minus, minus>(
+ const double amp_h1mm = jM2_W_extremal_qqbar<h1, minus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
- return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm;
+ 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
){
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 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))
)
+ 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./(
+ 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
){
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 COM cmsmu = COM{0., 3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
static constexpr COM cmumc = -1./6.;
return real(
cmsms*pow(abs(sym),2)
+cmumu*pow(abs(uncross),2)
+cmcmc*pow(abs(cross),2)
)
+2.*real(cmsmu*sym*conj(uncross))
+2.*real(cmsmc*sym*conj(cross))
+2.*real(cmumc*uncross*conj(cross))
;
}
} // 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, 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
){
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 COM cmsmu = COM{0.,3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
static constexpr COM cmumc = -1./6.;
return real(
cmsms*pow(abs(sym),2)
+cmumu*pow(abs(uncrossed),2)
+cmcmc*pow(abs(crossed),2)
)
+2.*real(cmsmu*sym*conj(uncrossed))
+2.*real(cmsmc*sym*conj(crossed))
+2.*real(cmumc*uncrossed*conj(crossed))
;
}
// matrix element for W + Jets with W emission *not* off central qqbar
double ME_W_Cenqqx_qq(
HLV pa, HLV pb, HLV pl,HLV plbar,
std::vector<HLV> partons, bool aqlinepa, bool aqlinepb,
bool qqxmarker, int nabove, int nbelow, bool forwards,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
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, 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/Zjets.cc b/src/Zjets.cc
index 3b32cbd..77be9e8 100644
--- a/src/Zjets.cc
+++ b/src/Zjets.cc
@@ -1,510 +1,511 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include <vector>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/jets.hh"
// generated headers
#include "HEJ/currents/jV_j.hh"
#include "HEJ/currents/jVuno_j.hh"
#include "HEJ/currents/jV_juno.hh"
namespace HEJ {
namespace currents {
namespace {
// Z propagator
COM ZProp(const double q, ParticleProperties const & zprop){
return 1. / (q - zprop.mass*zprop.mass + COM(0.,1.)*zprop.width*zprop.mass);
}
// Photon propagator
COM GProp(const double q) {
return 1. / q;
}
// Weak charge
double Zq (const ParticleID PID, const Helicity h, const double stw2, const double ctw) {
using namespace pid;
// Positive Spin
if (h == helicity::plus) {
// quarks
if (PID == d || PID == s || PID == b) return (+ 1.0 * stw2 / 3.0) / ctw;
if (PID == u || PID == c) return (- 2.0 * stw2 / 3.0) / ctw;
// antiquarks
if (PID == d_bar || PID == s_bar || PID == b_bar) return (+ 0.5 - 1.0 * stw2 / 3.0) / ctw;
if (PID == u_bar || PID == c_bar) return (- 0.5 + 2.0 * stw2 / 3.0) / ctw;
// electron
if (PID == electron) return stw2 / ctw;
}
// Negative Spin
else {
// quarks
if (PID == d || PID == s || PID == b) return (- 0.5 + 1.0 * stw2 / 3.0) / ctw;
if (PID == u || PID == c) return (+ 0.5 - 2.0 * stw2 / 3.0) / ctw;
// antiquarks
if (PID == d_bar || PID == s_bar || PID == b_bar) return (- 1.0 * stw2 / 3.0) / ctw;
if (PID == u_bar || PID == c_bar) return (+ 2.0 * stw2 / 3.0) / ctw;
// electron
if (PID == electron) return (-1.0 / 2.0 + stw2) / ctw;
}
throw std::logic_error("ERROR! No weak charge found");
}
// Electric charge
double Gq (const ParticleID PID) {
using namespace pid;
if (PID == d || PID == s || PID == b) return -1./3.;
if (PID == u || PID == c) return +2./3.;
if (PID == d_bar || PID == s_bar || PID == b_bar) return +1./3.;
if (PID == u_bar || PID == c_bar) return -2./3.;
throw std::logic_error("ERROR! No electric charge found");
}
//! Prefactor for Z+Jets Contributions
/**
* @brief Z+Jets Contributions Prefactor
* @param aptype Incoming Particle 1 type (Z emission)
- * @param PZ Z Propagator
- * @param PG Photon Propagator
+ * @param propZ Z Propagator
+ * @param propG Photon Propagator
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @param res Ouput: pref[h1][hl]
*
* Calculates prefactor for Z+Jets (includes couplings and propagators)
*/
- void Z_amp_pref(const ParticleID aptype, const COM PZ, const COM PG,
+ void Z_amp_pref(const ParticleID aptype,
+ COM const & propZ, COM const & propG,
const double stw2, const double ctw, COM (&res)[2][2]
){
using helicity::plus;
using helicity::minus;
- const double Zq_a_p = Zq(aptype, plus, stw2, ctw);
- const double Zq_a_m = Zq(aptype, minus, stw2, ctw);
- const double Ze_p = Zq(pid::electron, plus, stw2, ctw);
- const double Ze_m = Zq(pid::electron, minus, stw2, ctw);
- const double Gq_a = Gq(aptype);
+ const double zq_a_p = Zq(aptype, plus, stw2, ctw);
+ const double zq_a_m = Zq(aptype, minus, stw2, ctw);
+ const double ze_p = Zq(pid::electron, plus, stw2, ctw);
+ const double ze_m = Zq(pid::electron, minus, stw2, ctw);
+ const double gq_a = Gq(aptype);
- res[1][1] = -2.*(Zq_a_p * Ze_p * PZ - Gq_a * PG * stw2);
- res[1][0] = -2.*(Zq_a_p * Ze_m * PZ - Gq_a * PG * stw2);
- res[0][0] = -2.*(Zq_a_m * Ze_m * PZ - Gq_a * PG * stw2);
- res[0][1] = -2.*(Zq_a_m * Ze_p * PZ - Gq_a * PG * stw2);
+ res[1][1] = -2.*(zq_a_p * ze_p * propZ - gq_a * propG * stw2);
+ res[1][0] = -2.*(zq_a_p * ze_m * propZ - gq_a * propG * stw2);
+ res[0][0] = -2.*(zq_a_m * ze_m * propZ - gq_a * propG * stw2);
+ res[0][1] = -2.*(zq_a_m * ze_p * propZ - gq_a * propG * stw2);
}
//! Z+Jets FKL Contribution
/**
* @brief Z+Jets FKL Contribution
* @param pa Incoming Particle 1 (Z emission)
* @param pb Incoming Particle 2
* @param p1 Outgoing Particle 1 (Z emission)
* @param p2 Outgoing Particle 2
* @param pep Outgoing positron
* @param pem Outgoing electron
* @param res Ouput: jZ_j[h1][hl][h2]
*
* Calculates j_Z^\mu j_\mu.
*/
void jZ_j(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem, COM (&res)[2][2][2]
){
using helicity::plus;
using helicity::minus;
res[1][1][1] = jV_j<plus, plus, plus> (pa, p1, pb, p2, pem, pep);
res[1][1][0] = jV_j<plus, plus, minus>(pa, p1, pb, p2, pem, pep);
res[1][0][1] = jV_j<plus, minus, plus> (pa, p1, pb, p2, pem, pep);
res[1][0][0] = jV_j<plus, minus, minus>(pa, p1, pb, p2, pem, pep);
res[0][0][0] = conj(res[1][1][1]);
res[0][0][1] = conj(res[1][1][0]);
res[0][1][0] = conj(res[1][0][1]);
res[0][1][1] = conj(res[1][0][0]);
}
/**
* @brief Z+Jets Unordered Contribution, unordered on Z side
* @param pa Incoming Particle 1 (Z and Uno emission)
* @param pb Incoming Particle 2
* @param pg Unordered Gluon
* @param p1 Outgoing Particle 1 (Z and Uno emission)
* @param p2 Outgoing Particle 2
* @param pep Outgoing positron
* @param pem Outgoing electron
- * @param X Ouput: (U1-L)[h1][hl][h2][hg]
- * @param Y Ouput: (U2+L)[h1][hl][h2][hg]
+ * @param x Ouput: (U1-L)[h1][hl][h2][hg]
+ * @param y Ouput: (U2+L)[h1][hl][h2][hg]
*
* Calculates j_Z_{uno}^\mu j_\mu. Ie, unordered with Z emission same side.
*/
void jZuno_j(const HLV & pa, const HLV & pb, const HLV & pg, const HLV & p1, const HLV & p2,
- const HLV & pep, const HLV & pem, COM (&X)[2][2][2][2], COM (&Y)[2][2][2][2]
+ const HLV & pep, const HLV & pem, COM (&x)[2][2][2][2], COM (&y)[2][2][2][2]
){
using helicity::plus;
using helicity::minus;
- COM U1[2][2][2][2];
- U1[1][1][1][1] = HEJ::U1<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U1[1][1][1][0] = HEJ::U1<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- U1[1][1][0][1] = HEJ::U1<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U1[1][1][0][0] = HEJ::U1<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- U1[1][0][1][1] = HEJ::U1<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U1[1][0][1][0] = HEJ::U1<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- U1[1][0][0][1] = HEJ::U1<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U1[1][0][0][0] = HEJ::U1<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- U1[0][1][1][1] = HEJ::U1<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U1[0][1][1][0] = HEJ::U1<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- U1[0][1][0][1] = HEJ::U1<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U1[0][1][0][0] = HEJ::U1<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- U1[0][0][1][1] = HEJ::U1<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U1[0][0][1][0] = HEJ::U1<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- U1[0][0][0][1] = HEJ::U1<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U1[0][0][0][0] = HEJ::U1<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- COM U2[2][2][2][2];
- U2[1][1][1][1] = HEJ::U2<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U2[1][1][1][0] = HEJ::U2<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- U2[1][1][0][1] = HEJ::U2<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U2[1][1][0][0] = HEJ::U2<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- U2[1][0][1][1] = HEJ::U2<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U2[1][0][1][0] = HEJ::U2<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- U2[1][0][0][1] = HEJ::U2<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U2[1][0][0][0] = HEJ::U2<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- U2[0][1][1][1] = HEJ::U2<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U2[0][1][1][0] = HEJ::U2<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- U2[0][1][0][1] = HEJ::U2<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U2[0][1][0][0] = HEJ::U2<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- U2[0][0][1][1] = HEJ::U2<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U2[0][0][1][0] = HEJ::U2<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- U2[0][0][0][1] = HEJ::U2<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- U2[0][0][0][0] = HEJ::U2<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- COM L[2][2][2][2];
- L[1][1][1][1] = HEJ::L<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- L[1][1][1][0] = HEJ::L<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- L[1][1][0][1] = HEJ::L<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- L[1][1][0][0] = HEJ::L<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- L[1][0][1][1] = HEJ::L<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- L[1][0][1][0] = HEJ::L<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- L[1][0][0][1] = HEJ::L<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- L[1][0][0][0] = HEJ::L<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- L[0][1][1][1] = HEJ::L<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- L[0][1][1][0] = HEJ::L<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- L[0][1][0][1] = HEJ::L<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- L[0][1][0][0] = HEJ::L<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
-
- L[0][0][1][1] = HEJ::L<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
- L[0][0][1][0] = HEJ::L<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
- L[0][0][0][1] = HEJ::L<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
- L[0][0][0][0] = HEJ::L<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ COM u1[2][2][2][2];
+ u1[1][1][1][1] = U1<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u1[1][1][1][0] = U1<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ u1[1][1][0][1] = U1<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u1[1][1][0][0] = U1<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ u1[1][0][1][1] = U1<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u1[1][0][1][0] = U1<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ u1[1][0][0][1] = U1<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u1[1][0][0][0] = U1<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ u1[0][1][1][1] = U1<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u1[0][1][1][0] = U1<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ u1[0][1][0][1] = U1<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u1[0][1][0][0] = U1<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ u1[0][0][1][1] = U1<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u1[0][0][1][0] = U1<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ u1[0][0][0][1] = U1<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u1[0][0][0][0] = U1<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ COM u2[2][2][2][2];
+ u2[1][1][1][1] = U2<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u2[1][1][1][0] = U2<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ u2[1][1][0][1] = U2<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u2[1][1][0][0] = U2<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ u2[1][0][1][1] = U2<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u2[1][0][1][0] = U2<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ u2[1][0][0][1] = U2<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u2[1][0][0][0] = U2<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ u2[0][1][1][1] = U2<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u2[0][1][1][0] = U2<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ u2[0][1][0][1] = U2<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u2[0][1][0][0] = U2<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ u2[0][0][1][1] = U2<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u2[0][0][1][0] = U2<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ u2[0][0][0][1] = U2<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ u2[0][0][0][0] = U2<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ COM l[2][2][2][2];
+ l[1][1][1][1] = L<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ l[1][1][1][0] = L<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ l[1][1][0][1] = L<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ l[1][1][0][0] = L<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ l[1][0][1][1] = L<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ l[1][0][1][0] = L<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ l[1][0][0][1] = L<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ l[1][0][0][0] = L<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ l[0][1][1][1] = L<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ l[0][1][1][0] = L<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ l[0][1][0][1] = L<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ l[0][1][0][0] = L<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
+
+ l[0][0][1][1] = L<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ l[0][0][1][0] = L<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
+ l[0][0][0][1] = L<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
+ l[0][0][0][0] = L<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
for(int hg=0; hg<2; hg++){
- X[h1][hl][h2][hg] = U1[h1][hl][h2][hg] - L[h1][hl][h2][hg];
- Y[h1][hl][h2][hg] = U2[h1][hl][h2][hg] + L[h1][hl][h2][hg];
+ x[h1][hl][h2][hg] = u1[h1][hl][h2][hg] - l[h1][hl][h2][hg];
+ y[h1][hl][h2][hg] = u2[h1][hl][h2][hg] + l[h1][hl][h2][hg];
}
}
}
}
}
/**
* @brief Z+Jets Unordered Contribution, unordered opposite to Z side
* @param pa Incoming Particle 1 (Z emission)
* @param pb Incoming Particle 2 (unordered emission)
* @param p1 Outgoing Particle 1 (Z emission)
* @param p2 Outgoing Particle 2 (unordered emission)
* @param pg Unordered Gluon
* @param pep Outgoing positron
* @param pem Outgoing electron
- * @param X Ouput: (U1-L)[h1][hl][h2][hg]
- * @param Y Ouput: (U2+L)[h1][hl][h2][hg]
+ * @param x Ouput: (U1-L)[h1][hl][h2][hg]
+ * @param y Ouput: (U2+L)[h1][hl][h2][hg]
*
* Calculates j_Z^\mu j_{uno}_\mu. Ie, unordered with Z emission opposite side.
*/
void jZ_juno(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2, const HLV & pg,
- const HLV & pep, const HLV & pem, COM (&X)[2][2][2][2], COM (&Y)[2][2][2][2]
+ const HLV & pep, const HLV & pem, COM (&x)[2][2][2][2], COM (&y)[2][2][2][2]
){
using helicity::plus;
using helicity::minus;
- COM U1[2][2][2][2];
- U1[1][1][1][1] = HEJ::U1_jV<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U1[1][1][1][0] = HEJ::U1_jV<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- U1[1][1][0][1] = HEJ::U1_jV<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U1[1][1][0][0] = HEJ::U1_jV<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- U1[1][0][1][1] = HEJ::U1_jV<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U1[1][0][1][0] = HEJ::U1_jV<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- U1[1][0][0][1] = HEJ::U1_jV<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U1[1][0][0][0] = HEJ::U1_jV<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- U1[0][1][1][1] = HEJ::U1_jV<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U1[0][1][1][0] = HEJ::U1_jV<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- U1[0][1][0][1] = HEJ::U1_jV<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U1[0][1][0][0] = HEJ::U1_jV<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- U1[0][0][1][1] = HEJ::U1_jV<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U1[0][0][1][0] = HEJ::U1_jV<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- U1[0][0][0][1] = HEJ::U1_jV<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U1[0][0][0][0] = HEJ::U1_jV<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- COM U2[2][2][2][2];
- U2[1][1][1][1] = HEJ::U2_jV<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U2[1][1][1][0] = HEJ::U2_jV<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- U2[1][1][0][1] = HEJ::U2_jV<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U2[1][1][0][0] = HEJ::U2_jV<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- U2[1][0][1][1] = HEJ::U2_jV<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U2[1][0][1][0] = HEJ::U2_jV<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- U2[1][0][0][1] = HEJ::U2_jV<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U2[1][0][0][0] = HEJ::U2_jV<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- U2[0][1][1][1] = HEJ::U2_jV<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U2[0][1][1][0] = HEJ::U2_jV<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- U2[0][1][0][1] = HEJ::U2_jV<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U2[0][1][0][0] = HEJ::U2_jV<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- U2[0][0][1][1] = HEJ::U2_jV<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U2[0][0][1][0] = HEJ::U2_jV<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- U2[0][0][0][1] = HEJ::U2_jV<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- U2[0][0][0][0] = HEJ::U2_jV<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- COM L[2][2][2][2];
- L[1][1][1][1] = HEJ::L_jV<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- L[1][1][1][0] = HEJ::L_jV<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- L[1][1][0][1] = HEJ::L_jV<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- L[1][1][0][0] = HEJ::L_jV<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- L[1][0][1][1] = HEJ::L_jV<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- L[1][0][1][0] = HEJ::L_jV<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- L[1][0][0][1] = HEJ::L_jV<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- L[1][0][0][0] = HEJ::L_jV<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- L[0][1][1][1] = HEJ::L_jV<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- L[0][1][1][0] = HEJ::L_jV<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- L[0][1][0][1] = HEJ::L_jV<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- L[0][1][0][0] = HEJ::L_jV<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
-
- L[0][0][1][1] = HEJ::L_jV<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
- L[0][0][1][0] = HEJ::L_jV<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
- L[0][0][0][1] = HEJ::L_jV<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
- L[0][0][0][0] = HEJ::L_jV<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ COM u1[2][2][2][2];
+ u1[1][1][1][1] = U1_jV<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u1[1][1][1][0] = U1_jV<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ u1[1][1][0][1] = U1_jV<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u1[1][1][0][0] = U1_jV<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ u1[1][0][1][1] = U1_jV<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u1[1][0][1][0] = U1_jV<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ u1[1][0][0][1] = U1_jV<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u1[1][0][0][0] = U1_jV<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ u1[0][1][1][1] = U1_jV<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u1[0][1][1][0] = U1_jV<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ u1[0][1][0][1] = U1_jV<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u1[0][1][0][0] = U1_jV<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ u1[0][0][1][1] = U1_jV<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u1[0][0][1][0] = U1_jV<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ u1[0][0][0][1] = U1_jV<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u1[0][0][0][0] = U1_jV<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ COM u2[2][2][2][2];
+ u2[1][1][1][1] = U2_jV<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u2[1][1][1][0] = U2_jV<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ u2[1][1][0][1] = U2_jV<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u2[1][1][0][0] = U2_jV<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ u2[1][0][1][1] = U2_jV<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u2[1][0][1][0] = U2_jV<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ u2[1][0][0][1] = U2_jV<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u2[1][0][0][0] = U2_jV<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ u2[0][1][1][1] = U2_jV<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u2[0][1][1][0] = U2_jV<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ u2[0][1][0][1] = U2_jV<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u2[0][1][0][0] = U2_jV<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ u2[0][0][1][1] = U2_jV<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u2[0][0][1][0] = U2_jV<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ u2[0][0][0][1] = U2_jV<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ u2[0][0][0][0] = U2_jV<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ COM l[2][2][2][2];
+ l[1][1][1][1] = L_jV<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ l[1][1][1][0] = L_jV<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ l[1][1][0][1] = L_jV<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ l[1][1][0][0] = L_jV<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ l[1][0][1][1] = L_jV<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ l[1][0][1][0] = L_jV<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ l[1][0][0][1] = L_jV<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ l[1][0][0][0] = L_jV<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ l[0][1][1][1] = L_jV<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ l[0][1][1][0] = L_jV<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ l[0][1][0][1] = L_jV<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ l[0][1][0][0] = L_jV<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
+
+ l[0][0][1][1] = L_jV<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ l[0][0][1][0] = L_jV<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
+ l[0][0][0][1] = L_jV<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
+ l[0][0][0][0] = L_jV<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
for(int hg=0; hg<2; hg++){
- X[h1][hl][h2][hg] = U1[h1][hl][h2][hg] - L[h1][hl][h2][hg];
- Y[h1][hl][h2][hg] = U2[h1][hl][h2][hg] + L[h1][hl][h2][hg];
+ x[h1][hl][h2][hg] = u1[h1][hl][h2][hg] - l[h1][hl][h2][hg];
+ y[h1][hl][h2][hg] = u2[h1][hl][h2][hg] + l[h1][hl][h2][hg];
}
}
}
}
}
} // Anonymous Namespace
std::vector <double> ME_Z_qQ(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem,
const ParticleID aptype, const ParticleID bptype,
ParticleProperties const & zprop,
const double stw2, const double ctw
){
const HLV pZ = pep + pem;
- const COM PZ = ZProp(pZ.m2(), zprop);
- const COM PG = GProp(pZ.m2());
+ const COM propZ = ZProp(pZ.m2(), zprop);
+ const COM propG = GProp(pZ.m2());
COM pref_top[2][2], pref_bot[2][2];
- Z_amp_pref(aptype, PZ, PG, stw2, ctw, pref_top);
- Z_amp_pref(bptype, PZ, PG, stw2, ctw, pref_bot);
+ Z_amp_pref(aptype, propZ, propG, stw2, ctw, pref_top);
+ Z_amp_pref(bptype, propZ, propG, stw2, ctw, pref_bot);
- COM Coeff_top[2][2][2], Coeff_bot[2][2][2];
- jZ_j(pa, pb, p1, p2, pep, pem, Coeff_top);
- jZ_j(pb, pa, p2, p1, pep, pem, Coeff_bot);
+ COM coeff_top[2][2][2], coeff_bot[2][2][2];
+ jZ_j(pa, pb, p1, p2, pep, pem, coeff_top);
+ jZ_j(pb, pa, p2, p1, pep, pem, coeff_bot);
double sum_top=0., sum_bot=0., sum_mix=0.;
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
- const COM res_top = pref_top[h1][hl] * Coeff_top[h1][hl][h2];
- const COM res_bot = pref_bot[h2][hl] * Coeff_bot[h2][hl][h1];
+ const COM res_top = pref_top[h1][hl] * coeff_top[h1][hl][h2];
+ const COM res_bot = pref_bot[h2][hl] * coeff_bot[h2][hl][h1];
sum_top += norm(res_top);
sum_bot += norm(res_bot);
sum_mix += 2.0 * real(res_top * conj(res_bot));
}
}
}
const double t1_top = (pa-p1-pZ).m2();
const double t2_top = (pb-p2 ).m2();
const double t1_bot = (pa-p1 ).m2();
const double t2_bot = (pb-p2-pZ).m2();
sum_top /= t1_top * t2_top;
sum_bot /= t1_bot * t2_bot;
sum_mix /= sqrt(t1_top * t2_top * t1_bot * t2_bot);
// Colour factor: (CF*CA)/2
// Colour and helicity average: 1/(4*Nc^2)
const double pref = (C_F*C_A) / (8*N_C*N_C);
return {sum_top*pref, sum_bot*pref, sum_mix*pref};
}
double ME_Z_qg(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem,
const ParticleID aptype, const ParticleID /*bptype*/,
ParticleProperties const & zprop,
const double stw2, const double ctw
){
const HLV pZ = pep + pem;
- const COM PZ = ZProp(pZ.m2(), zprop);
- const COM PG = GProp(pZ.m2());
+ const COM propZ = ZProp(pZ.m2(), zprop);
+ const COM propG = GProp(pZ.m2());
- COM pref[2][2], Coeff[2][2][2];
- Z_amp_pref(aptype, PZ, PG, stw2, ctw, pref);
- jZ_j(pa, pb, p1, p2, pep, pem, Coeff);
+ COM pref[2][2], coeff[2][2][2];
+ Z_amp_pref(aptype, propZ, propG, stw2, ctw, pref);
+ jZ_j(pa, pb, p1, p2, pep, pem, coeff);
double sum = 0.;
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
- sum += norm(pref[h1][hl] * Coeff[h1][hl][h2]);
+ sum += norm(pref[h1][hl] * coeff[h1][hl][h2]);
}
}
}
sum /= (pa-p1-pZ).m2()*(pb-p2).m2();
// Colour factor: (CF*CA)/2
// Colour and helicity average: 1/(4*Nc^2)
// Divide by CF because of gluon (K_g -> CA)
sum *= C_A / (8*N_C*N_C);
// Multiply by CAM
return sum * K_g(p2, pb);
}
std::vector <double> ME_Zuno_qQ(const HLV & pa, const HLV & pb,
const HLV & pg, const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem,
const ParticleID aptype, const ParticleID bptype,
ParticleProperties const & zprop,
const double stw2, const double ctw
){
const HLV pZ = pep + pem;
- const COM PZ = ZProp(pZ.m2(), zprop);
- const COM PG = GProp(pZ.m2());
+ const COM propZ = ZProp(pZ.m2(), zprop);
+ const COM propG = GProp(pZ.m2());
COM prefact_top[2][2], prefact_bot[2][2];
- Z_amp_pref(aptype, PZ, PG, stw2, ctw, prefact_top);
- Z_amp_pref(bptype, PZ, PG, stw2, ctw, prefact_bot);
+ Z_amp_pref(aptype, propZ, propG, stw2, ctw, prefact_top);
+ Z_amp_pref(bptype, propZ, propG, stw2, ctw, prefact_bot);
- COM CoeffX_top[2][2][2][2], CoeffY_top[2][2][2][2];
- jZuno_j(pa, pb, pg, p1, p2, pep, pem, CoeffX_top, CoeffY_top);
+ COM coeffX_top[2][2][2][2], coeffY_top[2][2][2][2];
+ jZuno_j(pa, pb, pg, p1, p2, pep, pem, coeffX_top, coeffY_top);
- COM CoeffX_bot[2][2][2][2], CoeffY_bot[2][2][2][2];
- jZ_juno(pb, pa, p2, p1, pg, pep, pem, CoeffX_bot, CoeffY_bot);
+ COM coeffX_bot[2][2][2][2], coeffY_bot[2][2][2][2];
+ jZ_juno(pb, pa, p2, p1, pg, pep, pem, coeffX_bot, coeffY_bot);
double sum_top=0., sum_bot=0., sum_mix=0.;
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
for(int hg=0; hg<2; hg++){
const COM pref_top = prefact_top[h1][hl];
- const COM X_top = CoeffX_top[h1][hl][h2][hg];
- const COM Y_top = CoeffY_top[h1][hl][h2][hg];
+ const COM x_top = coeffX_top[h1][hl][h2][hg];
+ const COM y_top = coeffY_top[h1][hl][h2][hg];
const COM pref_bot = prefact_bot[h2][hl];
- const COM X_bot = CoeffX_bot[h2][hl][h1][hg];
- const COM Y_bot = CoeffY_bot[h2][hl][h1][hg];
-
- sum_top += norm(pref_top) * (C_A*C_F*C_F/2.*(norm(X_top)+norm(Y_top))
- - C_F/2.*(X_top*conj(Y_top)).real());
- sum_bot += norm(pref_bot) * (C_A*C_F*C_F/2.*(norm(X_bot)+norm(Y_bot))
- - C_F/2.*(X_bot*conj(Y_bot)).real());
-
- const COM XX = C_A*C_F*C_F/2. * pref_top * X_top * conj(pref_bot * X_bot);
- const COM YY = C_A*C_F*C_F/2. * pref_top * Y_top * conj(pref_bot * Y_bot);
- const COM XY = -C_F/2. * (pref_top * X_top * conj(pref_bot * Y_bot)
- + pref_top * Y_top * conj(pref_bot * X_bot));
- sum_mix += 2.0 * real(XX + YY + XY);
+ const COM x_bot = coeffX_bot[h2][hl][h1][hg];
+ const COM y_bot = coeffY_bot[h2][hl][h1][hg];
+
+ sum_top += norm(pref_top) * (C_A*C_F*C_F/2.*(norm(x_top)+norm(y_top))
+ - C_F/2.*(x_top*conj(y_top)).real());
+ sum_bot += norm(pref_bot) * (C_A*C_F*C_F/2.*(norm(x_bot)+norm(y_bot))
+ - C_F/2.*(x_bot*conj(y_bot)).real());
+
+ const COM xx = C_A*C_F*C_F/2. * pref_top * x_top * conj(pref_bot * x_bot);
+ const COM yy = C_A*C_F*C_F/2. * pref_top * y_top * conj(pref_bot * y_bot);
+ const COM xy = -C_F/2. * (pref_top * x_top * conj(pref_bot * y_bot)
+ + pref_top * y_top * conj(pref_bot * x_bot));
+ sum_mix += 2.0 * real(xx + yy + xy);
}
}
}
}
const double t1_top = (pa-pg-p1-pZ).m2();
const double t2_top = (pb-p2 ).m2();
const double t1_bot = (pa-pg-p1).m2();
const double t2_bot = (pb-p2-pZ).m2();
sum_top /= t1_top * t2_top;
sum_bot /= t1_bot * t2_bot;
sum_mix /= sqrt(t1_top * t2_top * t1_bot * t2_bot);
//Helicity sum and average over initial states
const double pref = 1./(4.*C_A*C_A);
return {sum_top*pref, sum_bot*pref, sum_mix*pref};
}
double ME_Zuno_qg(const HLV & pa, const HLV & pb,
const HLV & pg, const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem,
const ParticleID aptype, const ParticleID /*bptype*/,
ParticleProperties const & zprop,
const double stw2, const double ctw
){
const HLV pZ = pep + pem;
- const COM PZ = ZProp(pZ.m2(), zprop);
- const COM PG = GProp(pZ.m2());
+ const COM propZ = ZProp(pZ.m2(), zprop);
+ const COM propG = GProp(pZ.m2());
- COM pref[2][2], CoeffX[2][2][2][2], CoeffY[2][2][2][2];
- Z_amp_pref(aptype, PZ, PG, stw2, ctw, pref);
- jZuno_j(pa, pb, pg, p1, p2, pep, pem, CoeffX, CoeffY);
+ COM pref[2][2], coeffX[2][2][2][2], coeffY[2][2][2][2];
+ Z_amp_pref(aptype, propZ, propG, stw2, ctw, pref);
+ jZuno_j(pa, pb, pg, p1, p2, pep, pem, coeffX, coeffY);
double sum = 0.;
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
for(int hg=0; hg<2; hg++){
- const COM X = CoeffX[h1][hl][h2][hg];
- const COM Y = CoeffY[h1][hl][h2][hg];
- sum += norm(pref[h1][hl]) * (C_A*C_F*C_F/2.*(norm(X)+norm(Y))
- - C_F/2.*(X*conj(Y)).real());
+ const COM x = coeffX[h1][hl][h2][hg];
+ const COM y = coeffY[h1][hl][h2][hg];
+ sum += norm(pref[h1][hl]) * (C_A*C_F*C_F/2.*(norm(x)+norm(y))
+ - C_F/2.*(x*conj(y)).real());
}
}
}
}
sum /= (pa-pg-p1-pZ).m2()*(p2-pb).m2();
//Helicity sum and average over initial states
sum /= (4.*C_A*C_A);
// Multiply by CAM
return sum * (K_g(p2, pb) / C_F);
}
} // namespace currents
} // namespace HEJ

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:09 PM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804916
Default Alt Text
(63 KB)

Event Timeline