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 #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 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 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 #include #include #include #include #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; // --- 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(p1in,p1out,p2in,p2out,pl,plbar); + const COM ampp = jV_j(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 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(p1, p2, pa, pb, pg, pl, plbar); const COM u2 = U2(p1, p2, pa, pb, pg, pl, plbar); const COM l = L(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 - 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 + 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(pg, pq, plbar, pl, pqbar, p3, pb); - const COM u2Xcontr = U2X(pg, pq, plbar, pl, pqbar, p3, pb); - const COM lXcontr = LX(pg, pq, plbar, pl, pqbar, p3, pb); + const double ME2h1pp = jM2Wuno( + pg, p1, plbar, pl, pa, p2, pb, wprop + ); - const COM x = u1Xcontr - lXcontr; - const COM y = u2Xcontr + lXcontr; + const double ME2h1pm = jM2Wuno( + 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( + pg, p1, plbar, pl, pa, p2, pb, wprop + ); + + const double ME2h1mm = jM2Wuno( + 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(p1in,p1out,p2in,p2out,pl,plbar); - const COM ampp = jV_j(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( + 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 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(pa,p1,pb,p2,pg,pl,plbar); const COM u2 = U2_jV(pa,p1,pb,p2,pg,pl,plbar); const COM l = L_jV(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(p2in,p2out,p1in,p1out,pg,pl,plbar) + ampsq_jW_juno (p2in,p2out,p1in,p1out,pg,pl,plbar) + ampsq_jW_juno (p2in,p2out,p1in,p1out,pg,pl,plbar) + ampsq_jW_juno (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 + 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(pg, pq, plbar, pl, pqbar, p3, pb); + const COM u2Xcontr = U2X(pg, pq, plbar, pl, pqbar, p3, pb); + const COM lXcontr = LX(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 - 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( - pg, p1, plbar, pl, pa, p2, pb, wprop + const double amp_h1pp = jM2_W_extremal_qqbar( + pg, pq, plbar, pl, pqbar, p3, pb, wprop ); - - const double ME2h1pm = jM2Wuno( - pg, p1, plbar, pl, pa, p2, pb, wprop + const double amp_h1pm = jM2_W_extremal_qqbar( + pg, pq, plbar, pl, pqbar, p3, pb, wprop ); - - const double ME2h1mp = jM2Wuno( - pg, p1, plbar, pl, pa, p2, pb, wprop + const double amp_h1mp = jM2_W_extremal_qqbar( + pg, pq, plbar, pl, pqbar, p3, pb, wprop ); - - const double ME2h1mm = jM2Wuno( - pg, p1, plbar, pl, pa, p2, pb, wprop + const double amp_h1mm = jM2_W_extremal_qqbar( + 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( - pg, p1out, plbar, pl, p1in, p2out, p2in, wprop - )/(4.*C_A*C_A); - } - - // helicity sum helper for jWqqbar_j(...) - template - 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( - pg, pq, plbar, pl, pqbar, p3, pb, wprop - ); - const double amp_h1pm = jM2_W_extremal_qqbar( - pg, pq, plbar, pl, pqbar, p3, pb, wprop - ); - const double amp_h1mp = jM2_W_extremal_qqbar( - pg, pq, plbar, pl, pqbar, p3, pb, wprop - ); - const double amp_h1mm = jM2_W_extremal_qqbar( - 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( 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 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(pb,p2,p3,pa,p1,pl,plbar); const COM m2 = jW_qggm2(pb,p2,p3,pa,p1,pl,plbar); const COM m3 = jW_qggm3(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(pb,p2,p3,pa,p1,pl,plbar) + jW_jqqbar(pb,p2,p3,pa,p1,pl,plbar) ); } return prefactor*( jW_jqqbar(pb,p2,p3,pa,p1,pl,plbar) + jW_jqqbar(pb,p2,p3,pa,p1,pl,plbar) ); } namespace { // helper function for matrix element for W + Jets with central qqbar template 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( pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12 ); const COM cross = M_cross_W( pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12 ); const COM uncross = M_uncross_W( 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 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( pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_mp = amp_WCenqqbar_qq( pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_pm = amp_WCenqqbar_qq( pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_pp = amp_WCenqqbar_qq( 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 - 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 + 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( - pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 - ); - const COM uncrossed = M_qbar( - pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 - ); - const COM sym = M_sym( - pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 - ); + const COM crossed = M_cross( + pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 + ); + const COM uncrossed = M_qbar( + pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 + ); + const COM sym = M_sym( + 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 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( pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_mp = amp_W_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_pm = amp_W_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_pp = amp_W_Cenqqbar_qq( 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