diff --git a/current_generator/jW_j.frm b/current_generator/jW_j.frm
new file mode 100644
index 0000000..60d3930
--- /dev/null
+++ b/current_generator/jW_j.frm
@@ -0,0 +1,35 @@
+*/**
+* \brief Contraction of W current with FKL current
+*
+* \authors The HEJ collaboration (see AUTHORS for details)
+* \date 2020
+* \copyright GPLv2 or later
+*/
+#include- include/helspin.frm
+#include- include/write.frm
+
+v pb,p2,pa,p1,pg,pl,plbar,pW;
+i mu,nu,sigma,alpha;
+cf jW,epsg,m2,sqrt;
+s h;
+
+#do h2={+,-}
+ l [jW_j `h2'] = jW(-1, pa, p1, plbar, pl)*Current(`h2'1, p2, mu, pb);
+#enddo
+
+* eq:Weffcur1 in developer manual
+id jW(h?, pa?, p1?, plbar?, pl?) = Current(-1, pl, alpha, plbar)*(
+ + Current(h,p1,alpha,p1+pW,mu,pa)/m2(p1+pW)
+ + Current(h,p1,mu,pa-pW,alpha,pa)/m2(pa-pW)
+);
+
+multiply replace_(pW,pl+plbar);
+.sort
+#call ContractCurrents
+.sort
+format O4;
+format c;
+#call WriteHeader(`OUTPUT')
+#call WriteOptimised(`OUTPUT',jW_j,1,pa,p1,pb,p2,pl,plbar)
+#call WriteFooter(`OUTPUT')
+.end
diff --git a/src/Wjets.cc b/src/Wjets.cc
index 06f0c9f..5ecc097 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,782 +1,736 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <array>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/jets.hh"
#include "HEJ/Tensor.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/exceptions.hh"
// generated headers
+#include "HEJ/currents/jW_j.hh"
#include "HEJ/currents/jWuno_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/jW_juno.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 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 HEJ::C_A;
using HEJ::C_F;
const COM U1 = HEJ::U1<h1, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM U2 = HEJ::U2<h1, h2, pol>(p1, p2, pa, pb, pg, pl, plbar);
const COM L = HEJ::L<h1, 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)) - 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;
}
/**
* @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
){
using HEJ::C_A;
using HEJ::C_F;
const COM U1Xcontr = HEJ::U1X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM U2Xcontr = HEJ::U2X<h1, h2, hg>(pg, pq, plbar, pl, pqbar, p3, pb);
const COM LXcontr = HEJ::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)) - HEJ::C_F/2.*(X*conj(Y)).real();
amp *= WProp(plbar, pl, wprop);
// what do I have to put here?
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 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));
- }
+ // 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 = HEJ::jW_j<minus>(p1in,p1out,p2in,p2out,pl,plbar);
+ const COM ampp = HEJ::jW_j<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 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{
// 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
){
const COM U1 = HEJ::U1_jW<h2,hg>(pa,p1,pb,p2,pg,pl,plbar);
const COM U2 = HEJ::U2_jW<h2,hg>(pa,p1,pb,p2,pg,pl,plbar);
const COM L = HEJ::L_jW<h2,hg>(pa,p1,pb,p2,pg,pl,plbar);
return 2.*HEJ::C_F*std::real((L-U1)*std::conj(L+U2))
+ 2.*HEJ::C_F*HEJ::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.
*/
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;
// 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()));
}
}
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{
// 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 pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
HLV p2out, HLV 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.*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;
}
// 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>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double ME2h1pm = 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>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
const double ME2h1mm = jM2_W_extremal_qqbar<h1, minus, minus>(
pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm;
}
/**
* @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
){
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.*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 {
template<Helicity h1, Helicity hg>
double jW_jqqbar(
HLV const & pb,
HLV const & p2,
HLV const & p3,
HLV const & pa,
HLV const & p1,
HLV const & pl,
HLV const & plbar
) {
using helicity::minus;
using helicity::plus;
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 = HEJ::jW_qggm1<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
const COM m2 = HEJ::jW_qggm2<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
const COM m3 = HEJ::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)) ;
}
}
// contraction of W current with extremal qqbar on the other side
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;
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<HEJ::Helicity h1a, HEJ::Helicity h4b>
double amp_WCenqqx_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p4,
HLV const & pq, HLV const & pqbar,
HLV const & pl, HLV const & plbar,
HLV const & q11, HLV const & q12
) {
const COM sym = HEJ::M_sym_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM cross = HEJ::M_cross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM uncross = HEJ::M_uncross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
// Colour factors
static constexpr COM cmsms = 3.;
static constexpr COM cmumu = 4./3.;
static constexpr COM cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0., 3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
static constexpr COM cmumc = -1./6.;
return real(
cmsms*pow(abs(sym),2)
+cmumu*pow(abs(uncross),2)
+cmcmc*pow(abs(cross),2)
)
+2.*real(cmsmu*sym*conj(uncross))
+2.*real(cmsmc*sym*conj(cross))
+2.*real(cmumc*uncross*conj(cross))
;
}
}
// matrix element for W + Jets with W emission off central qqbar
double ME_WCenqqx_qq(
HLV pa, HLV pb, HLV pl, HLV plbar, std::vector<HLV> partons,
bool /* aqlinepa */, bool /* aqlinepb */, bool qqxmarker, int nabove,
ParticleProperties const & wprop
) {
using helicity::plus;
using helicity::minus;
CLHEP::HepLorentzVector q1 = pa;
for(int i = 0; i <= nabove; ++i){
q1 -= partons[i];
}
const auto qq = HEJ::split_into_lightlike(q1);
// since q1.m2() < 0 the following assertion is always true
// see documentation for split_into_lightlike
assert(qq.second.e() < 0);
HLV pq, pqbar;
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<HEJ::Helicity h1a, HEJ::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 = HEJ::M_cross<h1a, hqqbar>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
);
const COM uncrossed = HEJ::M_qbar<h1a, hqqbar>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
);
const COM sym = HEJ::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 = HEJ::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;
}