Page MenuHomeHEPForge
Files F19244387

Wjets.cc
No OneTemporary
Actions

Size
65 KB
Referenced Files
None
Subscribers
None

Wjets.cc
View Options

/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <array>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/jets.hh"
#include "HEJ/Tensor.hh"
using HEJ::Tensor;
using HEJ::init_sigma_index;
using HEJ::metric;
using HEJ::rank3_current;
using HEJ::rank5_current;
using HEJ::eps;
using HEJ::to_tensor;
using HEJ::Helicity;
using HEJ::angle;
using HEJ::square;
using HEJ::flip;
using HEJ::ParticleProperties;
namespace helicity = HEJ::helicity;
namespace { // Helper Functions
// FKL W Helper Functions
double WProp (const HLV & plbar, const HLV & pl, ParticleProperties const & wprop){
COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
+ COM(0.,1.)*wprop.mass*wprop.width);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
namespace {
// FKL current including W emission off negative helicities
// See eq. (87) {eq:jW-} in developer manual
// Note that the terms are rearranged
Tensor<1> jW_minus(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl
){
using HEJ::helicity::minus;
const double tWin = (pa-pl-plbar).m2();
const double tWout = (p1+pl+plbar).m2();
// C++ arithmetic operators are evaluated left-to-right,
// so the following first computes complex scalar coefficients,
// which then multiply a current, reducing the number
// of multiplications
return 2.*(
+ angle(p1, pl)*square(p1, plbar)/tWout
+ square(pa, plbar)*angle(pa, pl)/tWin
)*HEJ::current(p1, pa, helicity::minus)
+ 2.*angle(p1, pl)*square(pl, plbar)/tWout
*HEJ::current(pl, pa, helicity::minus)
+ 2.*square(pa, plbar)*angle(pl, plbar)/tWin
*HEJ::current(p1, plbar, helicity::minus);
}
}
// FKL current including W emission
// see eqs. (87), (88) {eq:jW-}, {eq:jW+} in developer manual
Tensor<1> jW(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl,
Helicity h
){
if(h == helicity::minus) {
return jW_minus(pa, p1, plbar, pl);
}
return jW_minus(pa, p1, pl, plbar).complex_conj();
}
/**
* @brief Contraction of the \tilde{U}_1 tensor
*
* This is the contraction of the tensor defined in eq:U1tensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM U1contr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
const HLV pW = pl+plbar;
const HLV q1g = pa-pW-p1-pg;
const HLV q1 = pa-p1-pW;
const HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
if(h1 == helicity::plus) {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity +++
return (4.*sqrt(2.)*(-1.*taW*angle(pa,pb)*square(p1,plbar)*
(s1W*angle(p1,pg)*angle(pg,pl)*square(p1,p2)*square(p2,pg) -
1.*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
((s1g + s1W)*angle(p1,pg)*square(p1,p2) +
s1g*(angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar)))) +
s1gW*s1W*angle(p1,pg)*angle(pa,pl)*pow(square(p1,p2),2.)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*square(p2,pg));
}
else {
// helicity ++-
return (sqrt(2.)*(taW*angle(pa,pb)*(4.*s1g*square(p1,plbar)*
(angle(p1,pl)*square(p1,pg)*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar)) +
angle(pl,plbar)*(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg) +
angle(pb,pl)*square(p2,pl) + angle(pb,plbar)*square(p2,plbar))*square(pg,plbar)) +
square(p1,pg)*(4.*angle(p1,pb)*square(p1,plbar)*
((s1g + s1W)*angle(p1,pl)*square(p1,p2) - 1.*s1W*angle(pg,pl)*square(p2,pg) +
s1W*angle(pl,plbar)*square(p2,plbar)) -
4.*s1W*angle(pb,pg)*(angle(p1,pl)*square(p1,p2) - 1.*angle(pg,pl)*square(p2,pg) +
angle(pl,plbar)*square(p2,plbar))*square(pg,plbar))) +
4.*s1gW*s1W*angle(pa,pl)*square(p1,pg)*
(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg))*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity +-+
return (sqrt(2.)*(4.*taW*angle(p2,pa)*square(p1,plbar)*
(s1W*angle(p1,pg)*angle(pg,pl)*square(p1,pb)*square(pb,pg) -
1.*(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar))*
((s1g + s1W)*angle(p1,pg)*square(p1,pb) +
s1g*(angle(pg,pl)*square(pb,pl) + angle(pg,plbar)*square(pb,plbar)))) -
4.*s1gW*s1W*angle(p1,pg)*angle(pa,pl)*pow(square(p1,pb),2.)*
(angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*square(pb,pg));
}
else {
// helicity +--
return (-1.*sqrt(2.)*(taW*angle(p2,pa)*(4.*s1g*square(p1,plbar)*
(angle(p1,pl)*square(p1,pg)*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar)) +
angle(pl,plbar)*(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg) +
angle(p2,pl)*square(pb,pl) + angle(p2,plbar)*square(pb,plbar))*square(pg,plbar)) +
square(p1,pg)*(4.*angle(p1,p2)*square(p1,plbar)*
((s1g + s1W)*angle(p1,pl)*square(p1,pb) - 1.*s1W*angle(pg,pl)*square(pb,pg) +
s1W*angle(pl,plbar)*square(pb,plbar)) -
4.*s1W*angle(p2,pg)*(angle(p1,pl)*square(p1,pb) - 1.*angle(pg,pl)*square(pb,pg) +
angle(pl,plbar)*square(pb,plbar))*square(pg,plbar))) +
4.*s1gW*s1W*angle(pa,pl)*square(p1,pg)*
(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg))*
(angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*angle(p2,pg));
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity -++
return (
sqrt(2.)*(-4.*s1gW*s1W*angle(p1,pg)*(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg))*
(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*square(pa,plbar) +
taW*angle(pg,pl)*square(p2,pa)*(4.*s1g*angle(p1,pl)*
(angle(p1,pb)*square(p1,p2) + angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar))*square(pl,plbar) +
4.*s1W*angle(p1,pg)*square(p2,pg)*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pg)*square(pg,plbar) -
1.*angle(pb,pl)*square(pl,plbar))) -
4.*taW*angle(p1,pg)*angle(p1,pl)*square(p2,pa)*
(square(p1,plbar)*((s1g + s1W)*angle(p1,pb)*square(p1,p2) +
s1g*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar))) -
1.*s1W*square(p1,p2)*(angle(pb,pg)*square(pg,plbar) + angle(pb,pl)*square(pl,plbar))))
)/(s1g*s1gW*s1W*taW*square(p2,pg));
}
else {
// helicity -+-
return (-1.*sqrt(2.)*(4.*s1W*angle(p1,pb)*square(p1,pg)*
(s1gW*angle(p1,pb)*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
square(pa,plbar) - 1.*taW*angle(p1,pl)*angle(pb,pg)*square(p2,pa)*square(pg,plbar)) +
taW*angle(p1,pl)*square(p2,pa)*((s1g + s1W)*angle(p1,pb)*square(p1,pg) +
s1g*(angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar)))*
(4.*angle(p1,pb)*square(p1,plbar) - 4.*angle(pb,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity --+
return (sqrt(2.)*(4.*s1gW*s1W*angle(p1,pg)*square(pa,plbar)*
(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg))*
(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) +
taW*square(pa,pb)*(-1.*angle(pg,pl)*
(4.*s1g*angle(p1,pl)*(angle(p1,p2)*square(p1,pb) + angle(p2,pg)*square(pb,pg) +
angle(p2,pl)*square(pb,pl) + angle(p2,plbar)*square(pb,plbar))*square(pl,plbar) +
4.*s1W*angle(p1,pg)*square(pb,pg)*
(angle(p1,p2)*square(p1,plbar) - 1.*angle(p2,pg)*square(pg,plbar) -
1.*angle(p2,pl)*square(pl,plbar))) +
4.*angle(p1,pg)*angle(p1,pl)*(square(p1,plbar)*
((s1g + s1W)*angle(p1,p2)*square(p1,pb) +
s1g*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar))) -
1.*s1W*square(p1,pb)*(angle(p2,pg)*square(pg,plbar) + angle(p2,pl)*square(pl,plbar)))))
)/(s1g*s1gW*s1W*taW*square(pb,pg));
}
else {
// helicity ---
return (sqrt(2.)*(s1W*angle(p1,p2)*square(p1,pg)*
(4.*s1gW*angle(p1,p2)*square(pa,plbar)*
(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) -
4.*taW*angle(p1,pl)*angle(p2,pg)*square(pa,pb)*square(pg,plbar)) +
taW*angle(p1,pl)*square(pa,pb)*((s1g + s1W)*angle(p1,p2)*square(p1,pg) +
s1g*(angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))*
(4.*angle(p1,p2)*square(p1,plbar) - 4.*angle(p2,pl)*square(pl,plbar))))/
(s1g*s1gW*s1W*taW*angle(p2,pg));
}
}
}
}
/**
* @brief Contraction of the \tilde{U}_2 tensor
*
* This is the contraction of the tensor defined in eq:U2tensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM U2contr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
const HLV pW = pl+plbar;
const HLV q1g=pa-pW-p1-pg;
const HLV q1 = pa-p1-pW;
const HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double s1W = (p1+pW).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
if(h1 == helicity::plus) {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity +++
return (-4.*sqrt(2.)*(taW*angle(pa,pg)*(taWg*square(p1,plbar)*
(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg))*
(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) -
1.*s1W*angle(pg,pl)*square(p1,p2)*square(p2,pg)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pg)*square(pg,plbar) +
angle(pb,pl)*square(pl,plbar))) +
s1W*angle(pa,pl)*square(p1,p2)*(tag*
(angle(pa,pg)*(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar))*square(pa,plbar) +
angle(pg,pl)*(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg) +
angle(pb,pl)*square(p2,pl) + angle(pb,plbar)*square(p2,plbar))*square(pl,plbar)) +
angle(pa,pg)*square(p2,pa)*((tag + taW)*angle(pa,pb)*square(pa,plbar) +
taW*(angle(pb,pg)*square(pg,plbar) + angle(pb,pl)*square(pl,plbar))))))/
(s1W*tag*taW*taWg*square(p2,pg));
}
else {
// helicity ++-
return (-4.*sqrt(2.)*(s1W*tag*angle(pa,pl)*square(p1,p2)*
(angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar))*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) -
1.*angle(pa,pb)*square(pa,pg)*(taW*taWg*angle(pa,pb)*square(p1,plbar)*
(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
s1W*angle(pa,pl)*square(p1,p2)*
(taW*angle(pb,pg)*square(pg,plbar) +
(tag + taW)*(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))))))/
(s1W*tag*taW*taWg*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity +-+
return (-4.*sqrt(2.)*(taW*angle(pa,pg)*(taWg*square(p1,plbar)*
(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg))*
(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar)) -
1.*s1W*square(p1,pb)*(angle(pa,pl)*square(pa,pb) + angle(pg,pl)*square(pb,pg))*
(angle(p2,pg)*square(pg,plbar) + angle(p2,pl)*square(pl,plbar))) +
s1W*square(p1,pb)*(tag*angle(pa,pl)*
(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar))*
(angle(pa,pg)*square(pa,plbar) + angle(pg,pl)*square(pl,plbar)) +
angle(p2,pa)*(angle(pa,pg)*square(pa,plbar)*
((tag + taW)*angle(pa,pl)*square(pa,pb) + taW*angle(pg,pl)*square(pb,pg)) +
tag*angle(pa,pl)*angle(pg,pl)*square(pa,pb)*square(pl,plbar)))))/
(s1W*tag*taW*taWg*square(pb,pg));
}
else {
// helicity +--
return (sqrt(2.)*(4.*taW*angle(p2,pa)*square(pa,pg)*
(taWg*angle(p2,pa)*square(p1,plbar)*
(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar)) -
1.*s1W*angle(p2,pg)*angle(pa,pl)*square(p1,pb)*square(pg,plbar)) +
s1W*angle(pa,pl)*square(p1,pb)*((tag + taW)*angle(p2,pa)*square(pa,pg) +
tag*(angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))*
(4.*angle(p2,pa)*square(pa,plbar) - 4.*angle(p2,pl)*square(pl,plbar))))/
(s1W*tag*taW*taWg*angle(p2,pg));
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity -++
return (sqrt(2.)*(-4.*s1W*angle(p1,pb)*(taW*angle(pa,pg)*angle(pg,pl)*square(p2,pa)*square(p2,pg) -
1.*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
((tag + taW)*angle(pa,pg)*square(p2,pa) +
tag*(angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))))*square(pa,plbar)
+ 4.*taW*taWg*angle(p1,pl)*angle(pa,pg)*pow(square(p2,pa),2.)*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar))))/
(s1W*tag*taW*taWg*square(p2,pg));
}
else {
// helicity -+-
return (sqrt(2.)*(4.*s1W*angle(p1,pb)*(tag*angle(pl,plbar)*
(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar))*square(pa,plbar)*square(pg,plbar) +
taW*angle(pg,pl)*square(p2,pg)*square(pa,pg)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pg)*square(pg,plbar)) -
1.*square(pa,pg)*((tag*angle(pa,pl)*
(angle(pb,pg)*square(p2,pg) + angle(pb,pl)*square(p2,pl) +
angle(pb,plbar)*square(p2,plbar)) +
angle(pa,pb)*((tag + taW)*angle(pa,pl)*square(p2,pa) +
taW*angle(pl,plbar)*square(p2,plbar)))*square(pa,plbar) +
taW*angle(pb,pg)*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
square(pg,plbar))) - 4.*taW*taWg*angle(p1,pl)*
(angle(pa,pb)*square(p2,pa) + angle(pb,pg)*square(p2,pg))*square(pa,pg)*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar))))/
(s1W*tag*taW*taWg*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity --+
return (4.*sqrt(2.)*(angle(p1,p2)*(taW*taWg*angle(p1,pl)*angle(pa,pg)*pow(square(pa,pb),2.)*
square(p1,plbar) + s1W*square(pa,plbar)*
(tag*(angle(pg,pl)*square(pb,pl) + angle(pg,plbar)*square(pb,plbar))*
(-1.*angle(pa,pl)*square(pa,pb) + angle(pl,plbar)*square(pb,plbar)) +
angle(pa,pg)*square(pa,pb)*((tag + taW)*angle(pa,pl)*square(pa,pb) +
taW*angle(pg,pl)*square(pb,pg) - 1.*(tag + taW)*angle(pl,plbar)*square(pb,plbar))))
- 1.*taW*taWg*angle(p1,pl)*angle(p2,pl)*angle(pa,pg)*pow(square(pa,pb),2.)*square(pl,plbar))
)/(s1W*tag*taW*taWg*square(pb,pg));
}
else {
// helicity ---
return (sqrt(2.)*(s1W*angle(p1,p2)*(-4.*square(pa,pg)*square(pa,plbar)*
(tag*angle(pa,pl)*(angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar)) +
angle(p2,pa)*((tag + taW)*angle(pa,pl)*square(pa,pb) +
taW*angle(pg,pl)*square(pb,pg) - 1.*taW*angle(pl,plbar)*square(pb,plbar))) +
4.*tag*angle(pl,plbar)*square(pa,plbar)*
(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg) + angle(p2,pl)*square(pb,pl) +
angle(p2,plbar)*square(pb,plbar))*square(pg,plbar) +
4.*taW*angle(p2,pg)*square(pa,pg)*
(angle(pa,pl)*square(pa,pb) + angle(pg,pl)*square(pb,pg) -
1.*angle(pl,plbar)*square(pb,plbar))*square(pg,plbar)) -
4.*taW*taWg*angle(p1,pl)*square(pa,pg)*
(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg))*
(angle(p1,p2)*square(p1,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
(s1W*tag*taW*taWg*angle(p2,pg));
}
}
}
}
/**
* @brief Contraction of the \tilde{L} tensor
*
* This is the contraction of the tensor defined in eq:Ltensor
* in the developer manual with the uno gluon polarisation vector,
* the leptonic and the partonic current (see eq:wunocurrent) in the
* developer manual)
*/
COM Lcontr(
HLV const & pg, HLV const & p1,HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity h1,
HLV const & p2, HLV const & pb, Helicity h2, Helicity hg
) {
const HLV pW = pl+plbar;
const HLV q1g=pa-pW-p1-pg;
const HLV q1 = pa-p1-pW;
const HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double s1W = (p1+pW).m2();
const double s2g = (p2+pg).m2();
const double sbg = (pb+pg).m2();
if(h1 == helicity::plus) {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity +++
return (-2.*sqrt(2.)*(angle(pb,pg)*q1g.m2()*square(p2,pb)*
(taW*angle(pa,pb)*square(p1,plbar)*
(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
s1W*angle(pa,pl)*square(p1,p2)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar))) +
2.*sbg*(taW*square(p1,plbar)*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
(angle(pa,pg)*angle(pb,pg)*square(p2,pg) +
angle(pa,pb)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))) +
s1W*angle(pa,pl)*square(p1,p2)*
((angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) + angle(pg,pl)*square(p2,pl) +
angle(pg,plbar)*square(p2,plbar))*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) +
angle(pb,pg)*square(p2,pg)*(angle(pa,pg)*square(pa,plbar) +
angle(pg,pl)*square(pl,plbar))))))/(s1W*sbg*taW*taW1*square(p2,pg));
}
else {
// helicity ++-
return (-1.*sqrt(2.)*(s2g*taW*angle(pa,pb)*square(p1,plbar)*
(4.*(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar))*
(angle(p1,pb)*square(p1,pg) - 1.*angle(pa,pb)*square(pa,pg) +
angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar)) +
4.*angle(pb,pg)*square(p2,pg)*(angle(p1,pl)*square(p1,pg) +
angle(pl,plbar)*square(pg,plbar))) +
s1W*s2g*angle(pa,pl)*(4.*angle(pb,pg)*square(p1,pg)*square(p2,pg) -
4.*angle(pa,pb)*square(p1,p2)*square(pa,pg) +
4.*square(p1,p2)*(angle(p1,pb)*square(p1,pg) + angle(pb,pl)*square(pg,pl) +
angle(pb,plbar)*square(pg,plbar)))*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)) -
2.*angle(p2,pb)*q1g.m2()*square(p2,pg)*
(taW*angle(pa,pb)*square(p1,plbar)*
(angle(p1,pl)*square(p1,p2) + angle(pl,plbar)*square(p2,plbar)) +
s1W*angle(pa,pl)*square(p1,p2)*
(angle(pa,pb)*square(pa,plbar) + angle(pb,pl)*square(pl,plbar)))))/
(s1W*s2g*taW*taW1*angle(pb,pg));
}
}
else {
if(hg == helicity::plus) {
// helicity +-+
return (-1.*sqrt(2.)*(2.*taW*square(p1,plbar)*(angle(p1,pl)*square(p1,pb) +
angle(pl,plbar)*square(pb,plbar))*
(2.*s2g*angle(pa,pg)*(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg)) +
angle(p2,pa)*(angle(p2,pg)*q1g.m2()*square(p2,pb) -
2.*s2g*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar)))) +
s1W*angle(pa,pl)*square(p1,pb)*(4.*s2g*square(pa,plbar)*
(angle(pa,pg)*(angle(p2,pa)*square(pa,pb) + angle(p2,pg)*square(pb,pg)) -
1.*angle(p2,pa)*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar))) +
s2g*(-4.*angle(p2,pl)*angle(pa,pg)*square(pa,pb) +
4.*angle(p2,pg)*angle(pg,pl)*square(pb,pg) +
4.*angle(p2,pl)*(angle(p1,pg)*square(p1,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar)))*square(pl,plbar) +
2.*angle(p2,pg)*q1g.m2()*square(p2,pb)*
(angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar)))))/
(s1W*s2g*taW*taW1*square(pb,pg));
}
else {
// helicity +--
return (sqrt(2.)*(2.*taW*angle(p2,pa)*square(p1,plbar)*
(2.*sbg*angle(p2,pg)*square(pb,pg)*
(angle(p1,pl)*square(p1,pg) + angle(pl,plbar)*square(pg,plbar)) +
(angle(p1,pl)*square(p1,pb) + angle(pl,plbar)*square(pb,plbar))*
(angle(p2,pb)*q1g.m2()*square(pb,pg) +
2.*sbg*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))) +
2.*s1W*angle(pa,pl)*(2.*sbg*angle(p1,p2)*square(p1,pb)*square(p1,pg) +
2.*sbg*angle(p2,pa)*square(p1,pb)*square(pa,pg) +
angle(p2,pb)*q1g.m2()*square(p1,pb)*square(pb,pg) +
2.*sbg*angle(p2,pg)*square(p1,pg)*square(pb,pg) +
2.*sbg*angle(p2,pl)*square(p1,pb)*square(pg,pl) +
2.*sbg*angle(p2,plbar)*square(p1,pb)*square(pg,plbar))*
(angle(p2,pa)*square(pa,plbar) - 1.*angle(p2,pl)*square(pl,plbar))))/
(s1W*sbg*taW*taW1*angle(p2,pg));
}
}
}
else {
if(h2 == helicity::plus) {
if(hg == helicity::plus) {
// helicity -++
return (sqrt(2.)*(2.*s1W*(angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
(2.*sbg*angle(p1,pg)*angle(pb,pg)*square(p2,pg) +
angle(p1,pb)*(angle(pb,pg)*q1g.m2()*square(p2,pb) +
2.*sbg*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))))*square(pa,plbar) +
taW*angle(p1,pl)*square(p2,pa)*(4.*sbg*square(p1,plbar)*
(angle(p1,pg)*angle(pb,pg)*square(p2,pg) +
angle(p1,pb)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar))) -
4.*sbg*(angle(pb,pg)*angle(pg,pl)*square(p2,pg) +
angle(pb,pl)*(angle(p1,pg)*square(p1,p2) + angle(pa,pg)*square(p2,pa) +
angle(pg,pl)*square(p2,pl) + angle(pg,plbar)*square(p2,plbar)))*square(pl,plbar) +
2.*angle(pb,pg)*q1g.m2()*square(p2,pb)*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar)))))/
(s1W*sbg*taW*taW1*square(p2,pg));
}
else {
// helicity -+-
return (sqrt(2.)*((angle(p1,pb)*square(pa,plbar)*
(((angle(pa,pl)*square(p2,pa) + angle(pl,plbar)*square(p2,plbar))*
(4.*angle(p1,pb)*square(p1,pg) - (2.*angle(p2,pb)*q1g.m2()*square(p2,pg))/s2g -
4.*angle(pa,pb)*square(pa,pg) + 4.*angle(pb,pl)*square(pg,pl) +
4.*angle(pb,plbar)*square(pg,plbar)))/angle(pb,pg) +
square(p2,pg)*(-4.*angle(pa,pl)*square(pa,pg) + 4.*angle(pl,plbar)*square(pg,plbar))))/
taW + (2.*angle(p1,pl)*(2.*s2g*angle(p1,pb)*square(p1,pg)*square(p2,pa) -
1.*angle(p2,pb)*q1g.m2()*square(p2,pa)*square(p2,pg) +
2.*s2g*(-1.*angle(pa,pb)*square(p2,pa)*square(pa,pg) -
1.*angle(pb,pg)*square(p2,pg)*square(pa,pg) +
square(p2,pa)*(angle(pb,pl)*square(pg,pl) + angle(pb,plbar)*square(pg,plbar))))*
(angle(p1,pb)*square(p1,plbar) - 1.*angle(pb,pl)*square(pl,plbar)))/(s1W*s2g*angle(pb,pg))
))/taW1;
}
}
else {
if(hg == helicity::plus) {
// helicity --+
return (2.*sqrt(2.)*(-1.*angle(p1,p2)*(2.*s2g*angle(p1,pg)*square(p1,pb) -
1.*angle(p2,pg)*q1g.m2()*square(p2,pb) +
2.*s2g*(-1.*angle(pa,pg)*square(pa,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar)))*
(taW*angle(p1,pl)*square(p1,plbar)*square(pa,pb) +
s1W*square(pa,plbar)*(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)))
+ taW*angle(p1,pl)*square(pa,pb)*
(angle(p2,pg)*(-1.*angle(p2,pl)*q1g.m2()*square(p2,pb) +
2.*s2g*angle(pg,pl)*square(pb,pg)) +
2.*s2g*angle(p2,pl)*(-1.*angle(pa,pg)*square(pa,pb) + angle(pg,pl)*square(pb,pl) +
angle(pg,plbar)*square(pb,plbar)))*square(pl,plbar) -
2.*s2g*angle(p1,pg)*(s1W*angle(p2,pg)*square(pa,plbar)*square(pb,pg)*
(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar)) +
taW*angle(p1,pl)*square(pa,pb)*
(angle(p2,pg)*square(p1,plbar)*square(pb,pg) -
1.*angle(p2,pl)*square(p1,pb)*square(pl,plbar)))))/(s1W*s2g*taW*taW1*square(pb,pg));
}
else {
// helicity ---
return (2.*sqrt(2.)*(-1.*angle(p1,p2)*(angle(p2,pb)*q1g.m2()*square(pb,pg)*
(taW*angle(p1,pl)*square(p1,plbar)*square(pa,pb) +
s1W*square(pa,plbar)*(angle(pa,pl)*square(pa,pb) - 1.*angle(pl,plbar)*square(pb,plbar))
) + 2.*sbg*(taW*angle(p1,pl)*square(p1,plbar) + s1W*angle(pa,pl)*square(pa,plbar))*
(angle(p2,pg)*square(pa,pg)*square(pb,pg) +
square(pa,pb)*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar))) -
2.*s1W*sbg*angle(pl,plbar)*square(pa,plbar)*
(angle(p2,pg)*square(pb,pg)*square(pg,plbar) +
square(pb,plbar)*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar)))) +
taW*angle(p1,pl)*angle(p2,pl)*(2.*sbg*angle(p2,pg)*square(pa,pg)*square(pb,pg) +
square(pa,pb)*(angle(p2,pb)*q1g.m2()*square(pb,pg) +
2.*sbg*(angle(p1,p2)*square(p1,pg) + angle(p2,pa)*square(pa,pg) +
angle(p2,pl)*square(pg,pl) + angle(p2,plbar)*square(pg,plbar))))*square(pl,plbar)))/
(s1W*sbg*taW*taW1*angle(p2,pg));
}
}
}
}
/**
* @brief W+Jets Unordered Contribution Helper Functions
* @returns result of equation (4.1.28) in Helen's Thesis (p.100)
*/
double jM2Wuno(HLV pg, HLV p1,HLV plbar, HLV pl, HLV pa, Helicity h1,
HLV p2, HLV pb, Helicity h2, Helicity pol,
ParticleProperties const & wprop
){
//@TODO Simplify the below (less Tensor class?)
init_sigma_index();
HLV pW = pl+plbar;
HLV q1g=pa-pW-p1-pg;
HLV q1 = pa-p1-pW;
HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
const double s2g = (p2+pg).m2();
const double sbg = (pb+pg).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
//use p1 as ref vec in pol tensor
Tensor<1> epsg = eps(pg,p2,pol);
Tensor<1> epsW = HEJ::current(pl,plbar,helicity::minus);
Tensor<1> j2b = HEJ::current(p2,pb,h2);
Tensor<1> Tq1q2 = to_tensor(
(pb/sbg + p2/s2g)*(q1 - pg).m2() + 2*q1 - pg
);
Tensor<3> J31a = rank3_current(p1, pa, h1);
Tensor<2> J2_qaW =J31a.contract((pa-pW)/taW/taW1, 2);
Tensor<2> J2_p1W =J31a.contract((p1+pW)/s1W/taW1, 2);
Tensor<3> L1a = outer(Tq1q2, J2_qaW);
Tensor<3> L1b = outer(Tq1q2, J2_p1W);
Tensor<3> L2a = outer(-2*pg, J2_qaW);
Tensor<3> L2b = outer(-2*pg, J2_p1W);
Tensor<3> L3 = outer(metric(), J2_qaW.contract(2*pg-q1,1)+J2_p1W.contract(2*pg-q1,2));
Tensor<3> L(0.);
Tensor<5> J51a = rank5_current(p1, pa, h1);
Tensor<4> J_qaW = J51a.contract((pa-pW),4);
Tensor<4> J_qag = J51a.contract(pa-pg,4);
Tensor<4> J_p1gW = J51a.contract(p1+pg+pW,4);
Tensor<3> U1a = J_qaW.contract(p1+pg,2);
Tensor<3> U1b = J_p1gW.contract(p1+pg,2);
Tensor<3> U1c = J_p1gW.contract(p1+pW,2);
Tensor<3> U1(0.);
Tensor<3> U2a = J_qaW.contract(pa-pg-pW,2);
Tensor<3> U2b = J_qag.contract(pa-pg-pW,2);
Tensor<3> U2c = J_qag.contract(p1+pW,2);
Tensor<3> U2(0.);
for(int nu=0; nu<4;++nu){
for(int mu=0;mu<4;++mu){
for(int rho=0;rho<4;++rho){
L(nu, mu, rho) = L1a(nu,mu,rho) + L1b(nu,rho,mu)
+ L2a(mu,nu,rho) + L2b(mu,rho,nu) + L3(mu,nu,rho);
U1(nu, mu, rho) = U1a(nu, mu, rho) / (s1g*taW)
+ U1b(nu,rho,mu) / (s1g*s1gW) + U1c(rho,nu,mu) / (s1W*s1gW);
U2(nu,mu,rho) = U2a(mu,nu,rho) / (taWg*taW)
+ U2b(mu,rho,nu) / (taWg*tag) + U2c(rho,mu,nu) / (s1W*tag);
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
double amp = HEJ::C_A*HEJ::C_F*HEJ::C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
double t1 = q1g.m2();
double t2 = q2.m2();
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
//Divide by t-channels
amp/=(t1*t2);
return amp;
}
/**
* @brief New implementation of unordered W+Jets current
*
* See eq:wunocurrent in the developer manual for the definition
* of this current
*
* The aim is to replace jM2Wuno and eventually all uses of the Tensor class.
* Explicit tensor contractions are rather slow and the initialisation code
* in Tensor.cc is not very transparent.
*
* At the moment, this function _only_ works for positive-energy spinors,
* so jM2Wuno has to be kept for qqbar emission.
*/
double jM2Wuno_pos_energy(
HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl,
HLV const & pa, Helicity const h1,
HLV const & p2, HLV const & pb, Helicity const h2, Helicity const pol,
ParticleProperties const & wprop
){
using HEJ::C_A;
using HEJ::C_F;
const COM U1 = U1contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM U2 = U2contr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM L = Lcontr(pg, p1, plbar, pl, pa, h1, p2, pb, h2, pol);
const COM X = U1 - L;
const COM Y = U2 + L;
double amp = C_A*C_F*C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
amp *= WProp(plbar, pl, wprop);
const HLV q1g = pa-pl-plbar-p1-pg;
const HLV q2 = p2-pb;
const double t1 = q1g.m2();
const double t2 = q2.m2();
amp /= (t1*t2);
return amp;
}
// Relevant Wqqx Helper Functions.
//g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
Tensor <1> gtqqxW(HLV pq,HLV pqbar,HLV pl,HLV plbar){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//blank 3 Gamma Current
Tensor<3> JV23 = rank3_current(pq,pqbar,helicity::minus);
// Components of g->qqW before W Contraction
Tensor<2> JV1 = JV23.contract((pq + pl + plbar),2)/(s2AB);
Tensor<2> JV2 = JV23.contract((pqbar + pl + plbar),2)/(s3AB);
// g->qqW Current. Note Minus between terms due to momentum flow.
// Also note: (-I)^2 from W vert. (I) from Quark prop.
Tensor<1> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
return JVCur;
}
// Helper Functions Calculate the Crossed Contribution
Tensor <2> MCrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar, HLV pl,
HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCross?
// Useful propagator factors
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
for(int i=0; i<nabove+1;++i){
q1=q1-partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double tcro1=(q3+pq).m2();
double tcro2=(q1-pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar,helicity::minus);
//Blank 5 gamma Current
Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_q3q = J523.contract((q3 + pq),2);
Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
// Components of Crossed Vertex Contribution
Tensor<3> Xcro1 = J_q3q.contract((pqbar + pl + plbar),3);
Tensor<3> Xcro2 = J_q3q.contract((q1 - pqbar),3);
Tensor<3> Xcro3 = J_2AB.contract((q1 - pqbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
Tensor<2> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
Tensor<2> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
//Initialise the Crossed Vertex Object
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xcro(mu,nu) = -(-Xcro1Cont(nu,mu)+Xcro2Cont(nu,mu)+Xcro3Cont(nu,mu));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar,
HLV pl, HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncross?
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
for(int i=0; i<nabove+1;++i){
q1=q1-partons.at(i);
}
q3 = q1 - pl - plbar - pq - pqbar;
double tunc1 = (q1-pq).m2();
double tunc2 = (q3+pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//Blank 5 gamma Current
Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
Tensor<4> J_q1q = J523.contract((q1 - pq),2);
// 2 Contractions taken care of.
Tensor<3> Xunc1 = J_2AB.contract((q3 + pqbar),3);
Tensor<3> Xunc2 = J_q1q.contract((q3 + pqbar),3);
Tensor<3> Xunc3 = J_q1q.contract((pqbar + pl + plbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
Tensor<2> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
Tensor<2> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
//Initialise the Uncrossed Vertex Object
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xunc(mu,nu) = -(- Xunc1Cont(mu,nu)+Xunc2Cont(mu,nu) +Xunc3Cont(mu,nu));
}
}
return Xunc;
}
// Helper Functions Calculate the g->qqxW (Eikonal) Contributions
Tensor <2> MSymW(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
HLV pl,HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MSym?
double sa2=(pa+pq).m2();
double s12=(p1+pq).m2();
double sa3=(pa+pqbar).m2();
double s13=(p1+pqbar).m2();
double saA=(pa+pl).m2();
double s1A=(p1+pl).m2();
double saB=(pa+plbar).m2();
double s1B=(p1+plbar).m2();
double sb2=(pb+pq).m2();
double s42=(p4+pq).m2();
double sb3=(pb+pqbar).m2();
double s43=(p4+pqbar).m2();
double sbA=(pb+pl).m2();
double s4A=(p4+pl).m2();
double sbB=(pb+plbar).m2();
double s4B=(p4+plbar).m2();
double s23AB=(pl+plbar+pq+pqbar).m2();
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3=q1-pq-pqbar-plbar-pl;
double t1 = (q1).m2();
double t3 = (q3).m2();
// g->qqW Current (Factors of sqrt2 dealt with in this function.)
Tensor<1> JV = gtqqxW(pq,pqbar,pl,plbar);
// 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13+s1A+s1B))
+ pa*(t1/(sa2+sa3+saA+saB)) );
Tensor<2> X1aCont = X1a.contract(JV,3);
//4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43+s4A+s4B))
+ pb*(t3/(sb2+sb3+sbA+sbB)) );
Tensor<2> X4bCont = X4b.contract(JV,3);
//Set up each term of 3G diagram.
Tensor<3> X3g1 = outer(q1+pq+pqbar+pl+plbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar-pl-plbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(JV,3);
Tensor<2> X3g2Cont = X3g2.contract(JV,2);
Tensor<2> X3g3Cont = X3g3.contract(JV,1);
// XSym is an amalgamation of x1a, X4b and X3g.
// Makes sense from a colour factor point of view.
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xsym(mu,nu) = (X3g1Cont(nu,mu) + X3g2Cont(mu,nu) - X3g3Cont(nu,mu))
+ (X1aCont(mu,nu) - X4bCont(mu,nu));
}
}
return Xsym/s23AB;
}
Tensor <2> MCross(HLV pa, HLV pq, HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2> XCroCont = J323.contract((q1-pqbar),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xcro(mu,nu) = XCroCont(nu,mu);
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncross(HLV pa, HLV pq,HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2> XUncCont = J323.contract((q1-pq),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xunc(mu,nu) = -XUncCont(mu,nu);
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2> MSym(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
std::vector<HLV> partons, Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MsymW?
HLV q1, q3;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3 = q1-pq-pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
double s23 = (pq+pqbar).m2();
double sa2 = (pa+pq).m2();
double sa3 = (pa+pqbar).m2();
double s12 = (p1+pq).m2();
double s13 = (p1+pqbar).m2();
double sb2 = (pb+pq).m2();
double sb3 = (pb+pqbar).m2();
double s42 = (p4+pq).m2();
double s43 = (p4+pqbar).m2();
Tensor<1> qqxCur = HEJ::current(pq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13))+ pa*(t1/(sa2+sa3)));
Tensor<2> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43)) + pb*(t3/(sb2+sb3)));
Tensor<2> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3> X3g1 = outer(q1+pq+pqbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;++mu){
for(int nu=0;nu<4;++nu){
Xsym(mu, nu) = COM(0,1) * ( (X3g1Cont(nu,mu) + X3g2Cont(mu,nu)
- X3g3Cont(nu,mu)) + (X1aCont(mu,nu) - X4bCont(mu,nu)) );
}
}
return Xsym/s23;
}
//! W+Jets FKL Contributions
/**
* @brief W+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_\mu.
* Handles all possible incoming states.
*/
double jW_j( HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
bool aqlineb, bool /* aqlinef */,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out);
const double WPropfact = WProp(plbar, pl, wprop);
const auto j_W = COM{0,-1}*jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
double Msqr = 0.;
for(const auto h: {plus, minus}) {
const auto j = HEJ::current(p2out, p2in, h);
Msqr += abs2(j_W.contract(j, 1));
}
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*Msqr/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
} // Anonymous Namespace
double ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop);
}
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true, wprop);
}
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop);
}
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true, wprop);
}
double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param pg Unordered Gluon momenta
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_{uno}_\mu. Ie, unordered with W emission opposite side.
* Handles all possible incoming states.
*/
double jW_juno(HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out,
HLV p2in, HLV pg, bool aqlineb, bool aqlinef,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out-pg);
const HLV q3=-(p2in-p2out);
const Helicity fhel = aqlinef?plus:minus;
const auto j_W = jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
const auto mj2p = HEJ::current(p2out, p2in, flip(fhel));
const auto mj2m = HEJ::current(p2out, p2in, fhel);
const auto jgbp = HEJ::current(pg, p2in, flip(fhel));
const auto jgbm = HEJ::current(pg, p2in, fhel);
const auto j2gp = HEJ::current(p2out, pg, flip(fhel));
const auto j2gm = HEJ::current(p2out, pg, fhel);
// Dot products of these which occur again and again
COM MWmp=j_W.dot(mj2p); // And now for the Higgs ones
COM MWmm=j_W.dot(mj2m);
const auto qsum = to_tensor(q2+q3);
const auto p1o = to_tensor(p1out);
const auto p1i = to_tensor(p1in);
const auto p2o = to_tensor(p2out);
const auto p2i = to_tensor(p2in);
const auto Lmm=( (-1.)*qsum*(MWmm) + (-2.*COM{j_W.dot(pg)})*mj2m + 2.*COM{mj2m.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmm/2. ) )/q3.m2();
const auto Lmp=( (-1.)*qsum*(MWmp) + (-2.*COM{j_W.dot(pg)})*mj2p + 2.*COM{mj2p.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmp/2. ) )/q3.m2();
const auto U1mm=(COM{jgbm.dot(j_W)}*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
const auto U1mp=(COM{jgbp.dot(j_W)}*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
const auto U2mm=((-1.)*COM{j2gm.dot(j_W)}*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
const auto U2mp=((-1.)*COM{j2gp.dot(j_W)}*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
double amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
ampsq*=WProp(plbar, pl, wprop);
return ampsq/((16)*(q2.m2()*q1.m2()));
}
}
double ME_W_unob_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true, wprop);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false, wprop);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true, wprop);
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same side.
* Handles all possible incoming states. Note this handles both forward and back-
* -ward Wuno emission. For forward, ensure p1out is the uno and W emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean up
* ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(HLV pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
HLV p2out, HLV p2in, bool aqlineb,
ParticleProperties const & wprop
){
//Calculate different Helicity choices
const Helicity h = aqlineb?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::plus,helicity::plus, wprop);
double ME2mpm = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::plus,helicity::minus, wprop);
double ME2mmp = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::minus,helicity::plus, wprop);
double ME2mmm = jM2Wuno_pos_energy(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states
return (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
}
}
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl,
ParticleProperties const & wprop
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop)
*K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Extremal qqx Contributions, function to handle all incoming types.
* @param pgin Incoming gluon which will split into qqx.
* @param pqout Quark of extremal qqx outgoing (W-Emission).
* @param plbar Outgoing anti-lepton momenta
* @param pl Outgoing lepton momenta
* @param pqbarout Anti-quark of extremal qqx pair. (W-Emission)
* @param pout Outgoing Particle 2 (end of FKL chain)
* @param p2in Incoming Particle 2
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W_{qqx} ^\mu j_\mu. Ie, Ex-QQX with W emission same side.
* Handles all possible incoming states. Calculated via crossing symmetry from jWuno_j.
*/
double jWqqx_j(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef,
ParticleProperties const & wprop
){
//Calculate Different Helicity Configurations.
const Helicity h = aqlinef?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::plus,helicity::plus, wprop);
double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::plus,helicity::minus, wprop);
double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::minus,helicity::plus, wprop);
double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states.
double ME2 = (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
}
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
}
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout, HLV plbar, HLV pl,
HLV pqbarout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in,
ParticleProperties const & wprop
){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm1(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below. (Less Tensor class use?)
double t1 = (p3-pb)*(p3-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p3-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,2)/t1;
return gqqCur*(-1);
}
//Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm2(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double t1 = (p2-pb)*(p2-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb,refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p2-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm3(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s23 = (p2+p3)*(p2+p3);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank1 = outer(p2+p3, metric())/s23;
Tensor<3> qqCurBlank2 = outer(pb, metric())/s23;
Tensor<1> Cur23 = HEJ::current(p2, p3,hel2);
Tensor<2> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2,
HLV p3,HLV plbar, HLV pl, bool aqlinepa,
ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
// 2 independent helicity choices (complex conjugation related).
Tensor<1> TMmmm1 = qggm1(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,minus,minus, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,minus,plus, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa?plus:minus);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1);
COM Mmmm2 = TMmmm2.contract(cur1a,1);
COM Mmmm3 = TMmmm3.contract(cur1a,1);
COM Mpmm1 = TMpmm1.contract(cur1a,1);
COM Mpmm2 = TMpmm2.contract(cur1a,1);
COM Mpmm3 = TMpmm3.contract(cur1a,1);
//Colour factors:
COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
cm1m1=8./3.;
cm2m2=8./3.;
cm3m3=6.;
cm1m2 =-1./3.;
cm1m3 = -3.*COM(0.,1.);
cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
double Mmmm = real( cm1m1*pow(abs(Mmmm1),2) + cm2m2*pow(abs(Mmmm2),2)
+ cm3m3*pow(abs(Mmmm3),2) + 2.*real(cm1m2*Mmmm1*conj(Mmmm2))
+ 2.*real(cm1m3*Mmmm1*conj(Mmmm3))
+ 2.*real(cm2m3*Mmmm2*conj(Mmmm3)) );
double Mpmm = real( cm1m1*pow(abs(Mpmm1),2) + cm2m2*pow(abs(Mpmm2),2)
+ cm3m3*pow(abs(Mpmm3),2) + 2.*real(cm1m2*Mpmm1*conj(Mpmm2))
+ 2.*real(cm1m3*Mpmm1*conj(Mpmm3))
+ 2.*real(cm2m3*Mpmm2*conj(Mpmm3)) );
// Divide by WProp
const double WPropfact = WProp(plbar, pl, wprop);
return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
}
// W+Jets qqxCentral
double ME_WCenqqx_qq(HLV pa, HLV pb,HLV pl, HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
ParticleProperties const & wprop
){
init_sigma_index();
HLV pq, pqbar, p1, p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am, T4bm, T1ap, T4bp;
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, helicity::plus);
T1am = HEJ::current(p1, pa, helicity::minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, helicity::plus);
T1am = HEJ::current(pa, p1, helicity::minus);}
if(!(aqlinepb)){
T4bp = HEJ::current(p4, pb, helicity::plus);
T4bm = HEJ::current(p4, pb, helicity::minus);}
else if(aqlinepb){
T4bp = HEJ::current(pb, p4, helicity::plus);
T4bm = HEJ::current(pb, p4, helicity::minus);}
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xsym = MSymW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
// (- - hel choice)
COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1));
COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1));
COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1));
COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1));
COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
return amp;
}
// no wqq emission
double ME_W_Cenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
int nbelow, bool forwards, ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
if (!forwards){ //If Emission from Leg a instead, flip process.
std::swap(pa, pb);
std::reverse(partons.begin(),partons.end());
std::swap(aqlinepa, aqlinepb);
qqxmarker = !qqxmarker;
std::swap(nabove,nbelow);
}
HLV pq, pqbar, p1,p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, plus);
T1am = HEJ::current(p1, pa, minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, plus);
T1am = HEJ::current(pa, p1, minus);}
Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb?plus:minus);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, minus, nabove);
Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, plus, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
amp*=WProp(plbar, pl, wprop);
return amp;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 4:41 AM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6486586
Default Alt Text
Wjets.cc (65 KB)

Event Timeline

HEPForge.org · .