Page MenuHomeHEPForge
Files F11222055

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/Wjets.hh b/include/HEJ/Wjets.hh
index 41e7218..6c1df54 100644
--- a/include/HEJ/Wjets.hh
+++ b/include/HEJ/Wjets.hh
@@ -1,182 +1,182 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2022
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in W+Jets.
*
* This file contains all the W+Jet specific components to compute
* the current contractions for valid HEJ processes.
*/
#pragma once
#include <vector>
#include "CLHEP/Vector/LorentzVector.h"
namespace HEJ {
struct ParticleProperties;
namespace currents {
using HLV = CLHEP::HepLorentzVector;
//! Square of qQ->qenuQ W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qQ->qenuQ
*
* This returns the square of the current contractions in qQ->qenuQ scattering
* with an emission of a W Boson off the quark q.
*/
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);
- //! qQg Wjets Unordered backwards opposite leg to W
+ //! W+uno same leg
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
- * @returns Square of the current contractions for qQ->qQg Scattering
+ * @returns Square of the current contractions for qQ->qQg
*
- * This returns the square of the current contractions in qQg->qQg scattering
+ * This returns the square of the current contractions in gqQ->gqQ scattering
* with an emission of a W Boson.
*/
- 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);
+ 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);
- //! W+uno same leg
+ //! qQg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param wprop Mass and width of the W boson
- * @returns Square of the current contractions for qQ->qQg
+ * @returns Square of the current contractions for qQ->qQg Scattering
*
- * This returns the square of the current contractions in gqQ->gqQ scattering
+ * This returns the square of the current contractions in qQg->qQg scattering
* with an emission of a W Boson.
*/
- 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);
+ 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);
//! W+Extremal qqbar. qqbar+Q
/**
* @param pgin Momentum of initial state gluon
* @param pqbarout Momentum of final state anti-quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqout Momentum of final state quark a
* @param p2out Momentum of initial state anti-quark b
* @param p2in Momentum of final state gluon b
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for ->qqbarQ
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqbar_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);
//! W+Extremal qqbar. gg->qqbarg. qqbar on forwards leg, W emission backwards leg.
/**
* @param pa Momentum of initial state (anti-)quark
* @param pb Momentum of initial state gluon
* @param p1 Momentum of final state (anti-)quark (after W emission)
* @param p2 Momentum of final state anti-quark
* @param p3 Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param aqlinepa Is opposite extremal leg to qqbar a quark or antiquark line
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for gq->qqbarqW
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated via current contraction of existing
* currents. Assumes qqbar split from forwards leg, W emission from backwards
* leg. Switch input (pa<->pb, p1<->pn) if calculating forwards qqbar.
*/
double ME_W_Exqqbar_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);
//! W+Jets qqbarCentral. qqbar W emission.
/**
* @param pa Momentum of initial state particle a
* @param pb Momentum of initial state particle b
* @param pl Momentum of final state lepton
* @param plbar Momentum of final state anti-lepton
* @param partons Vector of outgoing parton momenta
* @param aqlinepa True= pa is anti-quark
* @param aqlinepb True= pb is anti-quark
* @param qqbar_marker Ordering of the qqbar pair produced (q-qbar vs qbar-q)
* @param nabove Number of lipatov vertices "above" qqbar pair
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qq>qQQbarWq
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WCenqqbar_qq(HLV const & pa, HLV const & pb,
HLV const & pl, HLV const & plbar,
std::vector<HLV> const & partons,
bool aqlinepa, bool aqlinepb, bool qqbar_marker,
int nabove,
ParticleProperties const & wprop);
//! W+Jets qqbarCentral. W emission from backwards leg.
/**
* @param pa Momentum of initial state particle a
* @param pb Momentum of initial state particle b
* @param pl Momentum of final state lepton
* @param plbar Momentum of final state anti-lepton
* @param partons outgoing parton momenta
* @param aqlinepa True= pa is anti-quark
* @param aqlinepb True= pb is anti-quark
* @param qqbar_marker Ordering of the qqbar pair produced (q-qbar vs qbar-q)
* @param nabove Number of lipatov vertices "above" qqbar pair
* @param nbelow Number of lipatov vertices "below" qqbar pair
* @param forwards Swap to emission off front leg
* TODO: remove so args can be const
* @param wprop Mass and width of the W boson
* @returns Square of the current contractions for qq>qQQbarWq
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_W_Cenqqbar_qq(HLV pa, HLV pb,
HLV pl, HLV plbar,
std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqbar_marker,
int nabove, int nbelow, bool forwards,
ParticleProperties const & wprop);
} // namespace currents
} // namespace HEJ
diff --git a/src/Wjets.cc b/src/Wjets.cc
index dadae22..cb6a9b5 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,520 +1,522 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020-2022
* \copyright GPLv2 or later
*/
#include "HEJ/Wjets.hh"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <iostream>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/exceptions.hh"
// generated headers
#include "HEJ/currents/jV_j.hh"
#include "HEJ/currents/jV_juno.hh"
#include "HEJ/currents/jVuno_j.hh"
#include "HEJ/currents/jW_jqqbar.hh"
#include "HEJ/currents/jW_qqbar_j.hh"
#include "HEJ/currents/jWqqbar_j.hh"
#include "HEJ/currents/j_Wqqbar_j.hh"
namespace HEJ {
namespace currents {
namespace {
using COM = std::complex<double>;
// --- 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;
}
+ } // Anonymous Namespace
+
+ //! W+Jets FKL Contributions
+ 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
+ ){
+ using helicity::minus;
+ using helicity::plus;
+ 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);
+
+ return WProp(plbar, pl, wprop) * Msqr;
+ }
+
+ namespace {
/**
* @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;
// eq:VunoSumAveS in developer manual
// TODO: use same form as for other uno currents,
// extracting at least a factor of K_q = C_F
const double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
return WProp(plbar, pl, wprop) * 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
+ // 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 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 double ME2h1pp = jM2Wuno<h1, plus, plus>(
+ pg, p1, plbar, pl, pa, p2, pb, wprop
+ );
- const COM x = u1Xcontr - lXcontr;
- const COM y = u2Xcontr + lXcontr;
+ const double ME2h1pm = jM2Wuno<h1, plus, minus>(
+ pg, p1, plbar, pl, pa, p2, pb, wprop
+ );
- const double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
- return WProp(plbar, pl, wprop) * amp;
+ 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;
}
} // Anonymous Namespace
- //! W+Jets FKL Contributions
- 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
+ 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
){
- using helicity::minus;
- using helicity::plus;
-
- 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);
-
- return WProp(plbar, pl, wprop) * Msqr;
+ return jWuno_j_helsum<helicity::minus>(
+ pg, p1out, plbar, pl, p1in, p2out, p2in, wprop
+ )/(4.*C_A*C_A);
}
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);
const COM x = u1 - l;
const COM y = u2 + l;
return C_F*norm(x + y) - C_A*(x*conj(y)).real();
}
} // 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
){
using helicity::minus;
using helicity::plus;
// helicity assignments assume quarks
// in the end, this is irrelevant since we sum over all helicities
const double ampsq =
+ ampsq_jW_juno<minus, minus>(p2in,p2out,p1in,p1out,pg,pl,plbar)
+ ampsq_jW_juno<minus, plus> (p2in,p2out,p1in,p1out,pg,pl,plbar)
+ ampsq_jW_juno<plus, minus> (p2in,p2out,p1in,p1out,pg,pl,plbar)
+ ampsq_jW_juno<plus, plus> (p2in,p2out,p1in,p1out,pg,pl,plbar)
;
return WProp(plbar, pl, wprop)*ampsq;
}
namespace {
+ /**
+ * @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
+ ){
- // helicity sum helper for jWuno_j(...)
+ 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;
+
+ const double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real();
+ return WProp(plbar, pl, wprop) * amp;
+ }
+
+ // helicity sum helper for jWqqbar_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
+ double jWqqbar_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 = jM2Wuno<h1, plus, plus>(
- pg, p1, plbar, pl, pa, p2, pb, wprop
+ const double amp_h1pp = jM2_W_extremal_qqbar<h1, plus, plus>(
+ pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
-
- const double ME2h1pm = jM2Wuno<h1, plus, minus>(
- pg, p1, plbar, pl, pa, p2, pb, wprop
+ const double amp_h1pm = jM2_W_extremal_qqbar<h1, plus, minus>(
+ pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
-
- const double ME2h1mp = jM2Wuno<h1, minus, plus>(
- pg, p1, plbar, pl, pa, p2, pb, wprop
+ const double amp_h1mp = jM2_W_extremal_qqbar<h1, minus, plus>(
+ pg, pq, plbar, pl, pqbar, p3, pb, wprop
);
-
- const double ME2h1mm = jM2Wuno<h1, minus, minus>(
- pg, p1, plbar, pl, pa, p2, pb, wprop
+ 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;
}
-
} // 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_helsum<helicity::minus>(
- pg, p1out, plbar, pl, p1in, p2out, p2in, wprop
- )/(4.*C_A*C_A);
- }
-
- // helicity sum helper for jWqqbar_j(...)
- template<Helicity h1>
- double jWqqbar_j_helsum(
- HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl,
- HLV const & pqbar, HLV const & p3, HLV const & pb,
- ParticleProperties const & wprop
- ){
- using helicity::minus;
- using helicity::plus;
-
- const double amp_h1pp = jM2_W_extremal_qqbar<h1, plus, plus>(
- pg, pq, plbar, pl, pqbar, p3, pb, wprop
- );
- const double amp_h1pm = jM2_W_extremal_qqbar<h1, plus, minus>(
- pg, pq, plbar, pl, pqbar, p3, pb, wprop
- );
- const double amp_h1mp = jM2_W_extremal_qqbar<h1, minus, plus>(
- pg, pq, plbar, pl, pqbar, p3, pb, wprop
- );
- const double amp_h1mm = jM2_W_extremal_qqbar<h1, minus, minus>(
- pg, pq, plbar, pl, pqbar, p3, pb, wprop
- );
-
- return amp_h1pp + amp_h1pm + amp_h1mp + amp_h1mm;
- }
-
double ME_WExqqbar_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
){
//Helicity sum and average over initial states.
double ME2 = jWqqbar_j_helsum<helicity::plus>(
pgin, pqbarout, plbar, pl, pqout, p2out, p2in, wprop
)/(4.*C_A*C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
namespace {
template<Helicity h1, Helicity hg>
double jW_jqqbar(
HLV const & pb,
HLV const & p2,
HLV const & p3,
HLV const & pa,
HLV const & p1,
HLV const & pl,
HLV const & plbar
){
using std::norm;
static constexpr double cm1m1 = 8./3.;
static constexpr double cm2m2 = 8./3.;
static constexpr double cm3m3 = 6.;
static constexpr double cm1m2 =-1./3.;
static constexpr COM cm1m3 = COM{0.,-3.};
static constexpr COM cm2m3 = COM{0.,3.};
const COM m1 = jW_qggm1<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
const COM m2 = jW_qggm2<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
const COM m3 = jW_qggm3<h1,hg>(pb,p2,p3,pa,p1,pl,plbar);
return
+ cm1m1*norm(m1)
+ cm2m2*norm(m2)
+ cm3m3*norm(m3)
+ 2.*real(cm1m2*m1*conj(m2))
+ 2.*real(cm1m3*m1*conj(m3))
+ 2.*real(cm2m3*m2*conj(m3)) ;
}
} // Anonymous Namespace
// contraction of W current with extremal qqbar on the other side
double ME_W_Exqqbar_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.;
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_WCenqqbar_qq(
HLV const & pa, HLV const & p1,
HLV const & pb, HLV const & p4,
HLV const & pq, HLV const & pqbar,
HLV const & pl, HLV const & plbar,
HLV const & q11, HLV const & q12
){
using std::norm;
const COM sym = M_sym_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM cross = M_cross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
const COM uncross = M_uncross_W<h1a, h4b>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12
);
// Colour factors
static constexpr double cmsms = 3.;
static constexpr double cmumu = 4./3.;
static constexpr double cmcmc = 4./3.;
static constexpr COM cmsmu = COM{0., 3./2.};
static constexpr COM cmsmc = COM{0.,-3./2.};
static constexpr double cmumc = -1./6.;
return
+cmsms*norm(sym)
+cmumu*norm(uncross)
+cmcmc*norm(cross)
+2.*real(cmsmu*sym*conj(uncross))
+2.*real(cmsmc*sym*conj(cross))
+2.*real(cmumc*uncross*conj(cross))
;
}
} // Anonymous Namespace
// matrix element for W + Jets with W emission off central qqbar
double ME_WCenqqbar_qq(
HLV const & pa, HLV const & pb, HLV const & pl, HLV const & plbar,
std::vector<HLV> const & partons, bool /* aqlinepa */, bool /* aqlinepb */,
bool qqbar_marker, int nabove, ParticleProperties const & wprop
){
using helicity::plus;
using helicity::minus;
CLHEP::HepLorentzVector q1 = pa;
for(int i = 0; i <= nabove; ++i){
q1 -= partons[i];
}
const auto qq = split_into_lightlike(q1);
// since q1.m2() < 0 the following assertion is always true
// see documentation for split_into_lightlike
assert(qq.second.e() < 0);
HLV pq;
HLV pqbar;
if (qqbar_marker){
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_WCenqqbar_qq<minus, minus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_mp = amp_WCenqqbar_qq<minus, plus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pm = amp_WCenqqbar_qq<plus, minus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pp = amp_WCenqqbar_qq<plus, plus>(
pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second
);
// Divide by t channels, extremal + adjacent central vertex
const double ta = (pa-p1).m2();
const double t1 = q1.m2();
const double t3 = (q1-pl-plbar-pq-pqbar).m2();
const double tb = (p4-pb).m2();
const double amp = WProp(plbar, pl, wprop)*(
amp_mm+amp_mp+amp_pm+amp_pp
)/(9.*4.*ta*t1*t3*tb);
return amp;
}
- // helper function for matrix element for W + Jets with central qqbar
- // W emitted off extremal parton
- // @TODO: code duplication with amp_WCenqqbar_qq
- template<Helicity h1a, Helicity hqqbar>
- double amp_W_Cenqqbar_qq(
- HLV const & pa, HLV const & p1,
- HLV const & pb, HLV const & pn,
- HLV const & pq, HLV const & pqbar,
- HLV const & pl, HLV const & plbar,
- HLV const & q11, HLV const & q12
- ){
- using std::norm;
+ namespace {
+ // helper function for matrix element for W + Jets with central qqbar
+ // W emitted off extremal parton
+ // @TODO: code duplication with amp_WCenqqbar_qq
+ template<Helicity h1a, Helicity hqqbar>
+ double amp_W_Cenqqbar_qq(
+ HLV const & pa, HLV const & p1,
+ HLV const & pb, HLV const & pn,
+ HLV const & pq, HLV const & pqbar,
+ HLV const & pl, HLV const & plbar,
+ HLV const & q11, HLV const & q12
+ ){
+ using std::norm;
- const COM crossed = M_cross<h1a, hqqbar>(
- pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
- );
- const COM uncrossed = M_qbar<h1a, hqqbar>(
- pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
- );
- const COM sym = M_sym<h1a, hqqbar>(
- pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12
- );
+ 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 double cmsms = 3.;
- static constexpr double cmumu = 4./3.;
- static constexpr double cmcmc = 4./3.;
- static constexpr COM cmsmu = COM{0.,3./2.};
- static constexpr COM cmsmc = COM{0.,-3./2.};
- static constexpr double cmumc = -1./6.;
-
- return
- +cmsms*norm(sym)
- +cmumu*norm(uncrossed)
- +cmcmc*norm(crossed)
- +2.*real(cmsmu*sym*conj(uncrossed))
- +2.*real(cmsmc*sym*conj(crossed))
- +2.*real(cmumc*uncrossed*conj(crossed))
- ;
- }
+ //Colour factors:
+ static constexpr double cmsms = 3.;
+ static constexpr double cmumu = 4./3.;
+ static constexpr double cmcmc = 4./3.;
+ static constexpr COM cmsmu = COM{0.,3./2.};
+ static constexpr COM cmsmc = COM{0.,-3./2.};
+ static constexpr double cmumc = -1./6.;
+
+ return
+ +cmsms*norm(sym)
+ +cmumu*norm(uncrossed)
+ +cmcmc*norm(crossed)
+ +2.*real(cmsmu*sym*conj(uncrossed))
+ +2.*real(cmsmc*sym*conj(crossed))
+ +2.*real(cmumc*uncrossed*conj(crossed))
+ ;
+ }
+ } // Anonymous Namespace
// matrix element for W + Jets with W emission *not* off central qqbar
double ME_W_Cenqqbar_qq(
HLV pa, HLV pb, HLV pl,HLV plbar,
std::vector<HLV> partons, bool aqlinepa, bool aqlinepb,
bool qqbar_marker, 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);
qqbar_marker = !qqbar_marker;
std::swap(nabove,nbelow);
}
HLV q1=pa;
for(int i=0;i<nabove+1;++i){
q1-=partons.at(i);
}
const auto qq = split_into_lightlike(q1);
HLV pq;
HLV pqbar;
if (qqbar_marker){
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_Cenqqbar_qq<minus, minus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_mp = amp_W_Cenqqbar_qq<minus, plus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pm = amp_W_Cenqqbar_qq<plus, minus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
const double amp_pp = amp_W_Cenqqbar_qq<plus, plus>(
pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second
);
// Divide by t channels, extremal + adjacent central vertex
const double ta = (pa-p1).m2();
const double t1 = q1.m2();
const double t3 = (q1 - pq - pqbar).m2();
const double tb = (pn+pl+plbar-pb).m2();
const double amp= WProp(plbar, pl, wprop)*(
amp_mm+amp_mp+amp_pm+amp_pp
)/(9.*4.*ta*t1*t3*tb);
return amp;
}
} // namespace currents
} // namespace HEJ

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 11:20 AM (17 h, 14 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5047585
Default Alt Text
(31 KB)

Event Timeline