diff --git a/include/HEJ/Hjets.hh b/include/HEJ/Hjets.hh index 1486a70..f36e1d0 100644 --- a/include/HEJ/Hjets.hh +++ b/include/HEJ/Hjets.hh @@ -1,114 +1,111 @@ /** * \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 H+Jets. * * This file contains all the H+Jet specific components to compute - * the current contractions for valid HEJ processes, to form a full - * H+Jets ME, currently one would have to use functions from the - * jets.hh header also. We have FKL and also unordered components for - * H+Jets. + * the current contractions for valid HEJ processes. */ #pragma once #include "CLHEP/Vector/LorentzVector.h" namespace HEJ { namespace currents { using HLV = CLHEP::HepLorentzVector; /** * @brief Square of gq->Hq current contraction * @param ph Outgoing Higgs boson momentum * @param pg Incoming gluon momentum * @param pn Momentum of outgoing particle in FKL current * @param pb Momentum of incoming particle in FKL current * @param mt top mass (inf or value) * @param inc_bottom whether to include bottom mass effects (true) or not (false) * @param mb bottom mass (value) * @param vev Higgs vacuum expectation value * * Calculates helicity-averaged || \epsilon_\mu V_H^{\mu\nu} j_\nu ||^2. * See eq:S_gf_Hf in developer manual */ double ME_H_gq( HLV const & ph, HLV const & pa, HLV const & pn, HLV const & pb, double mt, bool inc_bottom, double mb, double vev ); /** * @brief Square of qg -> gqH current contraction * @param pg Outgoing (unordered) gluon momentum * @param ph Outgoing quark momentum * @param pa Incoming quark momentum * @param ph Outgoing Higgs boson momentum * @param pb Incoming gluon momentum * @param mt top mass (inf or value) * @param inc_bottom whether to include bottom mass effects (true) or not (false) * @param mb bottom mass (value) * @param vev Higgs vacuum expectation value * * Calculates helicity-averaged || j_{uno \mu} V_H^{\mu\nu} \epsilon_\nu ||^2. * See eq:S_gf_Hf in developer manual */ double ME_juno_jgH( HLV const & pg, HLV const & p1, HLV const & pa, HLV const & ph, HLV const & pb, double mt, bool inc_bottom, double mb, double vev ); //! Square of qQ->qHQ Higgs+Jets Scattering Current /** * @param p1out Momentum of final state quark * @param p1in Momentum of initial state quark * @param p2out Momentum of final state quark * @param p2in Momentum of intial state quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for qQ->qQ * * q~p1 Q~p2 (i.e. ALWAYS p1 for quark, p2 for quark) * should be called with qH1 meant to be contracted with p2 in first part of * vertex (i.e. if Q is backward, qH1 is forward) */ double ME_H_qQ(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qQ->qQg Higgs+Jets Unordered b Scattering Current /** * @param pg Momentum of unordered b gluon * @param p1out Momentum of final state quark * @param p1in Momentum of initial state quark * @param p2out Momentum of final state quark * @param p2in Momentum of intial state quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for qQ->qQg * * This construction is taking rapidity order: p1out >> p2out > pg */ double ME_H_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! @} } // namespace currents } // namespace HEJ diff --git a/include/HEJ/Wjets.hh b/include/HEJ/Wjets.hh index 95333c7..41e7218 100644 --- a/include/HEJ/Wjets.hh +++ b/include/HEJ/Wjets.hh @@ -1,185 +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, to form a full - * W+Jets ME, currently one would have to use functions from the - * jets.hh header also. We can calculate all subleading processes for - * W+Jets. + * 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 /** * @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 * * This returns the square of the current contractions in qQg->qQg 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); //! 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 * * This returns the square of the current contractions in gqQ->gqQ 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); //! 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/include/HEJ/jets.hh b/include/HEJ/jets.hh index bf03ea7..7cceff3 100644 --- a/include/HEJ/jets.hh +++ b/include/HEJ/jets.hh @@ -1,97 +1,93 @@ /** * \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 pure jets. * * This file contains all the necessary functions to compute the * current contractions for all valid pure jet HEJ processes, which * so far is FKL and unordered processes. It will also contain some * pure jet ME components used in other process ME calculations * */ #pragma once #include #include #include #include #include #include "CLHEP/Vector/LorentzVector.h" namespace HEJ { namespace currents { using HLV = CLHEP::HepLorentzVector; //! Helicity-summed current contraction /** * @param p1out Final-state momentum of first current * @param p1in Initial-state momentum of first current * @param p2out Final-state momentum of second current * @param p2in Initial-state momentum of second current * @returns Current contractions * * \internal * See eq:S in the developer manual */ double ME_qq( HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ); //! Helicity-summed current contraction with unordered gluon emission /** * @param pg Momentum of gluon emitted off first current * @param p1out Final-state momentum of first current * @param p1in Initial-state momentum of first current * @param p2out Final-state momentum of second current * @param p2in Initial-state momentum of second current * @returns Current contraction * * \internal * See eq:S_uno in the developer manual */ double ME_unob_qq( HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ); //! Helicity-summed current contraction with extremal quark/anti-quark emission /** * @param pa Initial-state momentum of first current * @param pb Initial-state momentum of gluon splitting into qqbar * @param p1 Final-state momentum of first current * @param p2 Less extremal (anti-)quark from splitting momentum * @param p3 More extremal (anti-)quark from splitting momentum * @returns Current contraction */ double ME_qqbar_qg( HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2, HLV const & p3 ); //! Square of qq->qQQbarq Pure Jets Central qqbar Scattering Current /** * @param pa Momentum of incoming leg a * @param pb Momentum of incoming leg b * @param partons vector of outgoing partons * @param qqbarmarker Is anti-quark further back in rapidity than quark * (qqbar pair) * @param nabove Number of gluons emitted above qqbar pair (back in rap) * @returns Square of the current contractions for qq->q+QQbar+q */ double ME_Cenqqbar_qq( HLV const & pa, HLV const & pb, std::vector const & partons, bool qqbarmarker, std::size_t nabove ); - - //! @TODO These are not currents and should be moved elsewhere. - double K_g(double p1minus, double paminus); - double K_g(HLV const & pout, HLV const & pin); } // namespace currents } // namespace HEJ diff --git a/src/Hjets.cc b/src/Hjets.cc index d70b079..db2bbd0 100644 --- a/src/Hjets.cc +++ b/src/Hjets.cc @@ -1,389 +1,388 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2022 * \copyright GPLv2 or later */ -#include "HEJ/jets.hh" #include "HEJ/Hjets.hh" #include #include #include #include #include #include "HEJ/ConfigFlags.hh" #include "HEJ/Constants.hh" #include "HEJ/LorentzVector.hh" #include "HEJ/utility.hh" // generated headers #include "HEJ/currents/j_h_j.hh" #include "HEJ/currents/jgh_j.hh" #include "HEJ/currents/juno_h_j.hh" #include "HEJ/currents/juno_jgh.hh" #ifdef HEJ_BUILD_WITH_QCDLOOP #include "qcdloop/qcdloop.h" #else #include "HEJ/exceptions.hh" #endif namespace HEJ { namespace currents { namespace { using COM = std::complex; constexpr double infinity = std::numeric_limits::infinity(); // NOLINT // Loop integrals #ifdef HEJ_BUILD_WITH_QCDLOOP const COM LOOPRWFACTOR = (COM(0.,1.)*M_PI*M_PI)/std::pow((2.*M_PI),4); COM B0DD(HLV const & q, double mq) { static std::vector> result(3); static auto ql_B0 = [](){ ql::Bubble,double,double> ql_B0; ql_B0.setCacheSize(100); return ql_B0; }(); static std::vector masses(2); static std::vector momenta(1); for(auto & m: masses) m = mq*mq; momenta.front() = q.m2(); ql_B0.integral(result, 1, masses, momenta); return result[0]; } COM C0DD(HLV const & q1, HLV const & q2, double mq) { static std::vector> result(3); static auto ql_C0 = [](){ ql::Triangle,double,double> ql_C0; ql_C0.setCacheSize(100); return ql_C0; }(); static std::vector masses(3); static std::vector momenta(3); for(auto & m: masses) m = mq*mq; momenta[0] = q1.m2(); momenta[1] = q2.m2(); momenta[2] = (q1+q2).m2(); ql_C0.integral(result, 1, masses, momenta); return result[0]; } // Kallen lambda functions, see eq:lambda in developer manual double lambda(const double s1, const double s2, const double s3) { return s1*s1 + s2*s2 + s3*s3 - 2*s1*s2 - 2*s1*s3 - 2*s2*s3; } // eq: T_1 in developer manual COM T1(HLV const & q1, HLV const & q2, const double m) { const double q12 = q1.m2(); const double q22 = q2.m2(); const HLV ph = q1 - q2; const double ph2 = ph.m2(); const double lam = lambda(q12, q22, ph2); assert(m > 0.); const double m2 = m*m; return - C0DD(q1, -q2, m)*(2.*m2 + 1./2.*(q12 + q22 - ph2) + 2.*q12*q22*ph2/lam) - (B0DD(q2, m) - B0DD(ph, m))*(q22 - q12 - ph2)*q22/lam - (B0DD(q1, m) - B0DD(ph, m))*(q12 - q22 - ph2)*q12/lam - 1.; } // eq: T_2 in developer manual COM T2(HLV const & q1, HLV const & q2, const double m) { const double q12 = q1.m2(); const double q22 = q2.m2(); const HLV ph = q1 - q2; const double ph2 = ph.m2(); const double lam = lambda(q12, q22, ph2); assert(m > 0.); const double m2 = m*m; return C0DD(q1, -q2, m)*( 4.*m2/lam*(ph2 - q12 - q22) - 1. - 4.*q12*q22/lam*( 1 + 3.*ph2*(q12 + q22 - ph2)/lam ) ) - (B0DD(q2, m) - B0DD(ph, m))*(1. + 6.*q12/lam*(q22 - q12 + ph2))*2.*q22/lam - (B0DD(q1, m) - B0DD(ph, m))*(1. + 6.*q22/lam*(q12 - q22 + ph2))*2.*q12/lam - 2.*(q12 + q22 - ph2)/lam; } #else // no QCDloop COM T1(HLV const & /*q1*/, HLV const & /*q2*/, double /*mt*/){ throw std::logic_error{"T1 called without QCDloop support"}; } COM T2(HLV const & /*q1*/, HLV const & /*q2*/, double /*mt*/){ throw std::logic_error{"T2 called without QCDloop support"}; } #endif // prefactors of g^{\mu \nu} and q_2^\mu q_1^\nu in Higgs boson emission vertex // see eq:VH in developer manual, but *without* global factor \alpha_s std::array TT( HLV const & qH1, HLV const & qH2, const double mt, const bool inc_bottom, const double mb, const double vev ) { if(mt == infinity) { std::array res = {qH1.dot(qH2), 1.}; for(auto & tt: res) tt /= (3.*M_PI*vev); return res; } std::array res = {T1(qH1, qH2, mt), T2(qH1, qH2, mt)}; for(auto & tt: res) tt *= mt*mt; if(inc_bottom) { res[0] += mb*mb*T1(qH1, qH2, mb); res[1] += mb*mb*T2(qH1, qH2, mb); } for(auto & tt: res) tt /= M_PI*vev; return res; } } // namespace double ME_H_qQ( HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, const double mt, const bool inc_bottom, const double mb, const double vev ){ using helicity::plus; using helicity::minus; const auto qqH1 = split_into_lightlike(qH1); const HLV qH11 = qqH1.first; const HLV qH12 = -qqH1.second; const auto qqH2 = split_into_lightlike(qH2); const HLV qH21 = qqH2.first; const HLV qH22 = -qqH2.second; // since qH1.m2(), qH2.m2() < 0 the following assertions are always true assert(qH11.e() > 0); assert(qH12.e() > 0); assert(qH21.e() > 0); assert(qH22.e() > 0); const auto T_pref = TT(qH1, qH2, mt, inc_bottom, mb, vev); const COM amp_mm = HEJ::j_h_j( p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1] ); const COM amp_mp = HEJ::j_h_j( p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1] ); const COM amp_pm = HEJ::j_h_j( p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1] ); const COM amp_pp = HEJ::j_h_j( p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1] ); // square of amplitudes, averaged over helicities return (norm(amp_mm) + norm(amp_mp) + norm(amp_pm) + norm(amp_pp)); } //@} double ME_H_gq( HLV const & ph, HLV const & pa, HLV const & pn, HLV const & pb, const double mt, const bool inc_bottom, const double mb, const double vev ){ using helicity::plus; using helicity::minus; const auto pH = split_into_lightlike(ph); const HLV ph1 = pH.first; const HLV ph2 = pH.second; // since pH.m2() > 0 the following assertions are always true assert(ph1.e() > 0); assert(ph2.e() > 0); const auto T_pref = TT(pa, pa-ph, mt, inc_bottom, mb, vev); // only distinguish between same and different helicities, // the other two combinations just add a factor of 2 const COM amp_mm = HEJ::jgh_j( pa, pb, pn, ph1, ph2, T_pref[0], T_pref[1] ); const COM amp_mp = HEJ::jgh_j( pa, pb, pn, ph1, ph2, T_pref[0], T_pref[1] ); constexpr double hel_factor = 2.; // sum over squares of helicity amplitudes return hel_factor*(norm(amp_mm) + norm(amp_mp)); } namespace { template double amp_juno_jgh( HLV const & pg, HLV const & p1, HLV const & pa, HLV const & ph1, HLV const & ph2, HLV const & pb, std::array const & T_pref ) { // TODO: code duplication with Wjets and pure jets const COM u1 = U1_jgh(pg, p1, pa, ph1, ph2, pb, T_pref[0], T_pref[1]); const COM u2 = U2_jgh(pg, p1, pa, ph1, ph2, pb, T_pref[0], T_pref[1]); const COM l = L_jgh(pg, p1, pa, ph1, ph2, pb, T_pref[0], T_pref[1]); return C_F*std::norm(u1+u2) - C_A*std::real((u1-l)*std::conj(u2+l)); } } // namespace double ME_juno_jgH( HLV const & pg, HLV const & p1, HLV const & pa, HLV const & ph, HLV const & pb, const double mt, const bool inc_bottom, const double mb, const double vev ) { using Helicity::plus; using Helicity::minus; const auto T_pref = TT(pb, pb-ph, mt, inc_bottom, mb, vev); const auto pH = split_into_lightlike(ph); const HLV ph1 = pH.first; const HLV ph2 = pH.second; // since pH.m2() > 0 the following assertions are always true assert(ph1.e() > 0); assert(ph2.e() > 0); // only 4 out of the 8 helicity amplitudes are independent // we still compute all of them for better numerical stability (mirror check) MultiArray amp; // NOLINTNEXTLINE #define ASSIGN_HEL(RES, J, H1, H2, HG) \ RES[H1][H2][HG] = J( \ pg, p1, pa, ph1, ph2, pb, T_pref \ ) ASSIGN_HEL(amp, amp_juno_jgh, minus, minus, minus); ASSIGN_HEL(amp, amp_juno_jgh, minus, minus, plus); ASSIGN_HEL(amp, amp_juno_jgh, minus, plus, minus); ASSIGN_HEL(amp, amp_juno_jgh, minus, plus, plus); ASSIGN_HEL(amp, amp_juno_jgh, plus, minus, minus); ASSIGN_HEL(amp, amp_juno_jgh, plus, minus, plus); ASSIGN_HEL(amp, amp_juno_jgh, plus, plus, minus); ASSIGN_HEL(amp, amp_juno_jgh, plus, plus, plus); #undef ASSIGN_HEL double ampsq = 0.; for(Helicity h1: {minus, plus}) { for(Helicity h2: {minus, plus}) { for(Helicity hg: {minus, plus}) { ampsq += amp[h1][h2][hg]; } } } return ampsq; } namespace { template double amp_juno_h_j( HLV const & pa, HLV const & pb, HLV const & pg, HLV const & p1, HLV const & p2, HLV const & qH11, HLV const & qH12, HLV const & qH21, HLV const & qH22, std::array const & T_pref ) { // TODO: code duplication with Wjets and pure jets const COM u1 = U1_h_j(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]); const COM u2 = U2_h_j(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]); const COM l = L_h_j(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]); const COM x = u1 - l; const COM y = u2 + l; return C_F*norm(x + y) - C_A*(x*conj(y)).real(); } } // namespace double ME_H_unob_qQ( HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, const double mt, const bool include_bottom, const double mb, const double vev ){ using helicity::plus; using helicity::minus; const auto qqH1 = split_into_lightlike(qH1); const HLV qH11 = qqH1.first; const HLV qH12 = -qqH1.second; const auto qqH2 = split_into_lightlike(qH2); const HLV qH21 = qqH2.first; const HLV qH22 = -qqH2.second; // since qH1.m2(), qH2.m2() < 0 the following assertions are always true assert(qH11.e() > 0); assert(qH12.e() > 0); assert(qH21.e() > 0); assert(qH22.e() > 0); const auto T_pref = TT(qH1, qH2, mt, include_bottom, mb, vev); // only 4 out of the 8 helicity amplitudes are independent // we still compute all of them for better numerical stability (mirror check) MultiArray amp{}; // NOLINTNEXTLINE #define ASSIGN_HEL(RES, J, H1, H2, HG) \ RES[H1][H2][HG] = J( \ p1in, p2in, pg, p1out, p2out, qH11, qH12, qH21, qH22, T_pref \ ) ASSIGN_HEL(amp, amp_juno_h_j, minus, minus, minus); ASSIGN_HEL(amp, amp_juno_h_j, minus, minus, plus); ASSIGN_HEL(amp, amp_juno_h_j, minus, plus, minus); ASSIGN_HEL(amp, amp_juno_h_j, minus, plus, plus); ASSIGN_HEL(amp, amp_juno_h_j, plus, minus, minus); ASSIGN_HEL(amp, amp_juno_h_j, plus, minus, plus); ASSIGN_HEL(amp, amp_juno_h_j, plus, plus, minus); ASSIGN_HEL(amp, amp_juno_h_j, plus, plus, plus); #undef ASSIGN_HEL double ampsq = 0.; for(Helicity h1: {minus, plus}) { for(Helicity h2: {minus, plus}) { for(Helicity hg: {minus, plus}) { ampsq += amp[h1][h2][hg]; } } } return ampsq; } //@} } // namespace currents } // namespace HEJ diff --git a/src/WWjets.cc b/src/WWjets.cc index d0928ef..cb50cc6 100644 --- a/src/WWjets.cc +++ b/src/WWjets.cc @@ -1,135 +1,134 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2020-2022 * \copyright GPLv2 or later */ #include "HEJ/WWjets.hh" #include #include #include #include #include #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/LorentzVector.hh" #include "HEJ/exceptions.hh" -#include "HEJ/jets.hh" // generated headers #include "HEJ/currents/jV_jV.hh" namespace HEJ { namespace currents { namespace { using std::conj; using COM = std::complex; // W Propagator 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; } //! WW+Jets FKL Contributions /** * @brief WW+Jets Currents * @param p1out Outgoing Particle 1 (W1 emission) * @param pl1bar Outgoing anti-lepton 1 momenta * @param pl1 Outgoing lepton 1 momenta * @param p1in Incoming Particle 1 (W1 emission) * @param p2out Outgoing Particle 2 (W2 emission) * @param pl2bar Outgoing anti-lepton 2 momenta * @param pl2 Outgoing lepton 2 momenta * @param p2in Incoming Particle 2 (W2 emission) * @param aqlineb Bool. Is Backwards quark line an anti-quark line? * @param aqlinef Bool. Is Forwards quark line an anti-quark line? * * Calculates jW^\mu jW_\mu */ template std::vector jW_jW( const HLV & p1out, const HLV & pl1bar, const HLV & pl1, const HLV & p1in, const HLV & p2out, const HLV & pl2bar, const HLV & pl2, const HLV & p2in, ParticleProperties const & wprop ){ using helicity::minus; const COM amp_top = jV_jV(p1in,p1out,p2in,p2out,pl1,pl1bar,pl2,pl2bar); const COM amp_bot = jV_jV(p1in,p1out,p2in,p2out,pl2,pl2bar,pl1,pl1bar); const double res_top = norm(amp_top); const double res_bot = norm(amp_bot); const double res_mix = 2.*real(amp_top*conj(amp_bot)); const double pref = WProp(pl1bar, pl1, wprop) * WProp(pl2bar, pl2, wprop); return {res_top*pref, res_bot*pref, res_mix*pref}; } } // namespace std::vector ME_WW_qQ( const HLV & p1out, const HLV& pl1bar, const HLV& pl1, const HLV & p1in, const HLV & p2out, const HLV& pl2bar, const HLV& pl2, const HLV & p2in, ParticleProperties const & wprop ){ using helicity::minus; return jW_jW( p1out, pl1bar, pl1, p1in, p2out, pl2bar, pl2, p2in, wprop ); } std::vector ME_WW_qbarQ( const HLV & p1out, const HLV& pl1bar, const HLV& pl1, const HLV & p1in, const HLV & p2out, const HLV& pl2bar, const HLV& pl2, const HLV & p2in, ParticleProperties const & wprop ){ using helicity::minus; using helicity::plus; return jW_jW( p1out, pl1bar, pl1, p1in, p2out, pl2bar, pl2, p2in, wprop ); } std::vector ME_WW_qQbar( const HLV & p1out, const HLV& pl1bar, const HLV& pl1, const HLV & p1in, const HLV & p2out, const HLV& pl2bar, const HLV& pl2, const HLV & p2in, ParticleProperties const & wprop ){ using helicity::minus; using helicity::plus; return jW_jW( p1out, pl1bar, pl1, p1in, p2out, pl2bar, pl2, p2in, wprop ); } std::vector ME_WW_qbarQbar( const HLV & p1out, const HLV& pl1bar, const HLV& pl1, const HLV & p1in, const HLV & p2out, const HLV& pl2bar, const HLV& pl2, const HLV & p2in, ParticleProperties const & wprop ){ using helicity::plus; return jW_jW( p1out, pl1bar, pl1, p1in, p2out, pl2bar, pl2, p2in, wprop ); } } // namespace currents } // namespace HEJ diff --git a/src/Wjets.cc b/src/Wjets.cc index cd73bfb..dadae22 100644 --- a/src/Wjets.cc +++ b/src/Wjets.cc @@ -1,521 +1,520 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2020 + * \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" -#include "HEJ/jets.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; } /** * @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 ){ 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; } } // 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 { // 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 { // 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 double ME2h1pp = jM2Wuno( pg, p1, plbar, pl, pa, p2, pb, wprop ); const double ME2h1pm = jM2Wuno( pg, p1, plbar, pl, pa, p2, pb, wprop ); 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 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; 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)) ; } // 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 diff --git a/src/Zjets.cc b/src/Zjets.cc index 51317d2..7ee4278 100644 --- a/src/Zjets.cc +++ b/src/Zjets.cc @@ -1,456 +1,456 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) - * \date 2020 + * \date 2020-2022 * \copyright GPLv2 or later */ #include #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/PDG_codes.hh" -#include "HEJ/jets.hh" #include "HEJ/utility.hh" +#include "HEJ/Zjets.hh" // generated headers #include "HEJ/currents/jV_j.hh" #include "HEJ/currents/jV_juno.hh" #include "HEJ/currents/jVuno_j.hh" namespace HEJ { namespace currents { namespace { using COM = std::complex; // Z propagator COM ZProp(const double q, ParticleProperties const & zprop){ return 1. / (q - zprop.mass*zprop.mass + COM(0.,1.)*zprop.width*zprop.mass); } // Photon propagator COM GProp(const double q) { return 1. / q; } // Weak charge template double Zq(ParticleID PID, double stw2, double ctw); // Weak charge - Positive Spin template<> double Zq( const ParticleID PID, const double stw2, const double ctw ) { using namespace pid; // quarks if (PID == d || PID == s || PID == b) return (+ 1.0 * stw2 / 3.0) / ctw; if (PID == u || PID == c) return (- 2.0 * stw2 / 3.0) / ctw; // antiquarks if (PID == d_bar || PID == s_bar || PID == b_bar) return (+ 0.5 - 1.0 * stw2 / 3.0) / ctw; if (PID == u_bar || PID == c_bar) return (- 0.5 + 2.0 * stw2 / 3.0) / ctw; // electron if (PID == electron) return stw2 / ctw; throw std::logic_error("ERROR! No weak charge found"); } // Weak charge - Negative Spin template<> double Zq( const ParticleID PID, const double stw2, const double ctw ) { using namespace pid; // quarks if (PID == d || PID == s || PID == b) return (- 0.5 + 1.0 * stw2 / 3.0) / ctw; if (PID == u || PID == c) return (+ 0.5 - 2.0 * stw2 / 3.0) / ctw; // antiquarks if (PID == d_bar || PID == s_bar || PID == b_bar) return (- 1.0 * stw2 / 3.0) / ctw; if (PID == u_bar || PID == c_bar) return (+ 2.0 * stw2 / 3.0) / ctw; // electron if (PID == electron) return (-1.0 / 2.0 + stw2) / ctw; throw std::logic_error("ERROR! No weak charge found"); } // Electric charge double Gq (const ParticleID PID) { using namespace pid; if (PID == d || PID == s || PID == b) return -1./3.; if (PID == u || PID == c) return +2./3.; if (PID == d_bar || PID == s_bar || PID == b_bar) return +1./3.; if (PID == u_bar || PID == c_bar) return -2./3.; throw std::logic_error("ERROR! No electric charge found"); } //! Prefactor for Z+Jets Contributions /** * @brief Z+Jets Contributions Prefactor * @param aptype Incoming Particle 1 type (Z emission) * @param propZ Z Propagator * @param propG Photon Propagator * @param stw2 Value of sin(theta_w)^2 * @param ctw Value of cos(theta_w) * @returns Prefactors for Z+Jets for all helicity combinations * (includes couplings and propagators) */ MultiArray Z_amp_pref( const ParticleID aptype, COM const & propZ, COM const & propG, const double stw2, const double ctw ){ using helicity::plus; using helicity::minus; const double zq_a_p = Zq(aptype, stw2, ctw); const double zq_a_m = Zq(aptype, stw2, ctw); const double ze_p = Zq(pid::electron, stw2, ctw); const double ze_m = Zq(pid::electron, stw2, ctw); const double gq_a = Gq(aptype); MultiArray res; res[ plus][ plus] = -2.*(zq_a_p * ze_p * propZ - gq_a * propG * stw2); res[ plus][minus] = -2.*(zq_a_p * ze_m * propZ - gq_a * propG * stw2); res[minus][minus] = -2.*(zq_a_m * ze_m * propZ - gq_a * propG * stw2); res[minus][ plus] = -2.*(zq_a_m * ze_p * propZ - gq_a * propG * stw2); return res; } //! Z+Jets FKL Contribution /** * @brief Z+Jets FKL Contribution * @param pa Incoming Particle 1 (Z emission) * @param pb Incoming Particle 2 * @param p1 Outgoing Particle 1 (Z emission) * @param p2 Outgoing Particle 2 * @param pep Outgoing positron * @param pem Outgoing electron * @returns j_Z^\mu j_\mu for all helicities h1, hl, h2 */ MultiArray jZ_j( const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2, const HLV & pep, const HLV & pem ){ using helicity::plus; using helicity::minus; MultiArray res; // NOLINTNEXTLINE #define ASSIGN_HEL(RES, J, H1, HL, H2) \ RES[H1][HL][H2] = J(pa, p1, pb, p2, pem, pep) ASSIGN_HEL(res, jV_j, plus, minus, minus); ASSIGN_HEL(res, jV_j, plus, minus, plus); ASSIGN_HEL(res, jV_j, plus, plus, minus); ASSIGN_HEL(res, jV_j, plus, plus, plus); #undef ASSIGN_HEL for(auto hl: {minus, plus}) { for(auto h2: {minus, plus}) { res[minus][hl][h2] = std::conj(res[plus][flip(hl)][flip(h2)]); } } return res; } // X and Y as used in contractions with unordered currents struct XY { COM X; COM Y; }; /** * @brief Z+Jets Unordered Contribution, unordered on Z side * @tparam h1 Helicity of line 1 (Z emission line) * @tparam hl Lepton Helicity * @tparam h2 Helicity of line 2 * @tparam hg Helicity of unordered gluon * @param pa Incoming Particle 1 (Z and Uno emission) * @param pb Incoming Particle 2 * @param pg Unordered Gluon * @param p1 Outgoing Particle 1 (Z and Uno emission) * @param p2 Outgoing Particle 2 * @param pep Outgoing positron * @param pem Outgoing electron * @returns X: (U1-L), Y: (U2+l) * * Calculates j_Z_{uno}^\mu j_\mu. Ie, unordered with Z emission same side. */ template XY amp_jZuno_j( const HLV & pa, const HLV & pb, const HLV & pg, const HLV & p1, const HLV & p2, const HLV & pep, const HLV & pem ){ const COM u1 = U1(p1, p2, pa, pb, pg, pem, pep); const COM u2 = U2(p1, p2, pa, pb, pg, pem, pep); const COM l = L (p1, p2, pa, pb, pg, pem, pep); return {u1 - l, u2 + l}; } MultiArray jZuno_j( const HLV & pa, const HLV & pb, const HLV & pg, const HLV & p1, const HLV & p2, const HLV & pep, const HLV & pem ){ using helicity::plus; using helicity::minus; MultiArray xy; // NOLINTNEXTLINE #define ASSIGN_HEL(XY, J, H1, H2, H3, H4) \ XY[H1][H2][H3][H4] = J(pa, pb, pg, p1, p2, pep, pem) // NOLINT ASSIGN_HEL(xy, amp_jZuno_j, minus, minus, minus, minus); ASSIGN_HEL(xy, amp_jZuno_j, minus, minus, minus, plus); ASSIGN_HEL(xy, amp_jZuno_j, minus, minus, plus, minus); ASSIGN_HEL(xy, amp_jZuno_j, minus, minus, plus, plus); ASSIGN_HEL(xy, amp_jZuno_j, minus, plus, minus, minus); ASSIGN_HEL(xy, amp_jZuno_j, minus, plus, minus, plus); ASSIGN_HEL(xy, amp_jZuno_j, minus, plus, plus, minus); ASSIGN_HEL(xy, amp_jZuno_j, minus, plus, plus, plus); ASSIGN_HEL(xy, amp_jZuno_j, plus, minus, minus, minus); ASSIGN_HEL(xy, amp_jZuno_j, plus, minus, minus, plus); ASSIGN_HEL(xy, amp_jZuno_j, plus, minus, plus, minus); ASSIGN_HEL(xy, amp_jZuno_j, plus, minus, plus, plus); ASSIGN_HEL(xy, amp_jZuno_j, plus, plus, minus, minus); ASSIGN_HEL(xy, amp_jZuno_j, plus, plus, minus, plus); ASSIGN_HEL(xy, amp_jZuno_j, plus, plus, plus, minus); ASSIGN_HEL(xy, amp_jZuno_j, plus, plus, plus, plus); #undef ASSIGN_HEL return xy; } /** * @brief Z+Jets Unordered Contribution, unordered opposite to Z side * @tparam h1 Helicity of line 1 (Z emission) * @tparam hl Lepton Helicity * @tparam h2 Helicity of line 2 (unordered emission) * @tparam hg Helicity of unordered gluon * @param pa Incoming Particle 1 (Z emission) * @param pb Incoming Particle 2 (unordered emission) * @param p1 Outgoing Particle 1 (Z emission) * @param p2 Outgoing Particle 2 (unordered emission) * @param pg Unordered Gluon * @param pep Outgoing positron * @param pem Outgoing electron * @returns X: (U1-L), Y: (U2+l) * * Calculates j_Z^\mu j_{uno}_\mu. Ie, unordered with Z emission opposite side. */ template XY amp_jZ_juno( const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2, const HLV & pg, const HLV & pep, const HLV & pem ){ const COM u1 = U1_jV(pa, p1, pb, p2, pg, pem, pep); const COM u2 = U2_jV(pa, p1, pb, p2, pg, pem, pep); const COM l = L_jV (pa, p1, pb, p2, pg, pem, pep); return {u1 - l, u2 + l}; } MultiArray jZ_juno( const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2, const HLV & pg, const HLV & pep, const HLV & pem ){ using helicity::plus; using helicity::minus; MultiArray xy; // NOLINTNEXTLINE #define ASSIGN_HEL(XY, J, H1, H2, H3, H4) \ XY[H1][H2][H3][H4] = J(pa, pb, p1, p2, pg, pep, pem) ASSIGN_HEL(xy, amp_jZ_juno, minus, minus, minus, minus); ASSIGN_HEL(xy, amp_jZ_juno, minus, minus, minus, plus); ASSIGN_HEL(xy, amp_jZ_juno, minus, minus, plus, minus); ASSIGN_HEL(xy, amp_jZ_juno, minus, minus, plus, plus); ASSIGN_HEL(xy, amp_jZ_juno, minus, plus, minus, minus); ASSIGN_HEL(xy, amp_jZ_juno, minus, plus, minus, plus); ASSIGN_HEL(xy, amp_jZ_juno, minus, plus, plus, minus); ASSIGN_HEL(xy, amp_jZ_juno, minus, plus, plus, plus); ASSIGN_HEL(xy, amp_jZ_juno, plus, minus, minus, minus); ASSIGN_HEL(xy, amp_jZ_juno, plus, minus, minus, plus); ASSIGN_HEL(xy, amp_jZ_juno, plus, minus, plus, minus); ASSIGN_HEL(xy, amp_jZ_juno, plus, minus, plus, plus); ASSIGN_HEL(xy, amp_jZ_juno, plus, plus, minus, minus); ASSIGN_HEL(xy, amp_jZ_juno, plus, plus, minus, plus); ASSIGN_HEL(xy, amp_jZ_juno, plus, plus, plus, minus); ASSIGN_HEL(xy, amp_jZ_juno, plus, plus, plus, plus); #undef ASSIGN_HEL return xy; } } // Anonymous Namespace std::vector ME_Z_qQ(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2, const HLV & pep, const HLV & pem, const ParticleID aptype, const ParticleID bptype, ParticleProperties const & zprop, const double stw2, const double ctw ){ using helicity::minus; using helicity::plus; const HLV pZ = pep + pem; const COM propZ = ZProp(pZ.m2(), zprop); const COM propG = GProp(pZ.m2()); MultiArray pref_top = Z_amp_pref(aptype, propZ, propG, stw2, ctw); MultiArray pref_bot = Z_amp_pref(bptype, propZ, propG, stw2, ctw); MultiArray coeff_top = jZ_j(pa, pb, p1, p2, pep, pem); MultiArray coeff_bot = jZ_j(pb, pa, p2, p1, pep, pem); double sum_top=0.; double sum_bot=0.; double sum_mix=0.; for(auto h1: {minus, plus}){ for(auto hl: {minus, plus}){ for(auto h2: {minus, plus}){ const COM res_top = pref_top[h1][hl] * coeff_top[h1][hl][h2]; const COM res_bot = pref_bot[h2][hl] * coeff_bot[h2][hl][h1]; sum_top += norm(res_top); sum_bot += norm(res_bot); sum_mix += 2.0 * real(res_top * conj(res_bot)); } } } return {sum_top, sum_bot, sum_mix}; } double ME_Z_qg(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2, const HLV & pep, const HLV & pem, const ParticleID aptype, const ParticleID /*bptype*/, ParticleProperties const & zprop, const double stw2, const double ctw ){ using helicity::minus; using helicity::plus; const HLV pZ = pep + pem; const COM propZ = ZProp(pZ.m2(), zprop); const COM propG = GProp(pZ.m2()); MultiArray pref = Z_amp_pref(aptype, propZ, propG, stw2, ctw); MultiArray coeff = jZ_j(pa, pb, p1, p2, pep, pem); double sum = 0.; for(auto h1: {minus, plus}){ for(auto hl: {minus, plus}){ for(auto h2: {minus, plus}){ sum += norm(pref[h1][hl] * coeff[h1][hl][h2]); } } } return sum; } std::vector ME_Zuno_qQ(const HLV & pa, const HLV & pb, const HLV & pg, const HLV & p1, const HLV & p2, const HLV & pep, const HLV & pem, const ParticleID aptype, const ParticleID bptype, ParticleProperties const & zprop, const double stw2, const double ctw ){ using helicity::minus; using helicity::plus; const HLV pZ = pep + pem; const COM propZ = ZProp(pZ.m2(), zprop); const COM propG = GProp(pZ.m2()); MultiArray prefact_top = Z_amp_pref(aptype, propZ, propG, stw2, ctw); MultiArray prefact_bot = Z_amp_pref(bptype, propZ, propG, stw2, ctw); const MultiArray coeff_top = jZuno_j(pa, pb, pg, p1, p2, pep, pem); const MultiArray coeff_bot = jZ_juno(pb, pa, p2, p1, pg, pep, pem); double sum_top=0.; double sum_bot=0.; double sum_mix=0.; for(auto h1: {minus, plus}){ for(auto hl: {minus, plus}){ for(auto h2: {minus, plus}){ for(auto hg: {minus, plus}){ const COM pref_top = prefact_top[h1][hl]; const COM x_top = coeff_top[h1][hl][h2][hg].X; const COM y_top = coeff_top[h1][hl][h2][hg].Y; const COM pref_bot = prefact_bot[h2][hl]; const COM x_bot = coeff_bot[h2][hl][h1][hg].X; const COM y_bot = coeff_bot[h2][hl][h1][hg].Y; sum_top += norm(pref_top) * (C_A*C_F*C_F/2.*(norm(x_top)+norm(y_top)) - C_F/2.*(x_top*conj(y_top)).real()); sum_bot += norm(pref_bot) * (C_A*C_F*C_F/2.*(norm(x_bot)+norm(y_bot)) - C_F/2.*(x_bot*conj(y_bot)).real()); const COM xx = C_A*C_F*C_F/2. * pref_top * x_top * conj(pref_bot * x_bot); const COM yy = C_A*C_F*C_F/2. * pref_top * y_top * conj(pref_bot * y_bot); const COM xy = -C_F/2. * (pref_top * x_top * conj(pref_bot * y_bot) + pref_top * y_top * conj(pref_bot * x_bot)); sum_mix += 2.0 * real(xx + yy + xy); } } } } //Helicity sum and average over initial states const double pref = 1./(4.*C_A*C_A); return {sum_top*pref, sum_bot*pref, sum_mix*pref}; } double ME_Zuno_qg(const HLV & pa, const HLV & pb, const HLV & pg, const HLV & p1, const HLV & p2, const HLV & pep, const HLV & pem, const ParticleID aptype, const ParticleID /*bptype*/, ParticleProperties const & zprop, const double stw2, const double ctw ){ using helicity::minus; using helicity::plus; const HLV pZ = pep + pem; const COM propZ = ZProp(pZ.m2(), zprop); const COM propG = GProp(pZ.m2()); MultiArray pref = Z_amp_pref(aptype, propZ, propG, stw2, ctw); const auto coeff = jZuno_j(pa, pb, pg, p1, p2, pep, pem); double sum = 0.; for(auto h1: {minus, plus}){ for(auto hl: {minus, plus}){ for(auto h2: {minus, plus}){ for(auto hg: {minus, plus}){ const COM X = coeff[h1][hl][h2][hg].X; const COM Y = coeff[h1][hl][h2][hg].Y; sum += norm(pref[h1][hl]) * (C_A*C_F*C_F/2.*(norm(X)+norm(Y)) - C_F/2.*(X*conj(Y)).real()); } } } } //Helicity sum and average over initial states return sum / (4.*C_A*C_A); } } // namespace currents } // namespace HEJ diff --git a/src/jets.cc b/src/jets.cc index fa136eb..b006d8b 100644 --- a/src/jets.cc +++ b/src/jets.cc @@ -1,244 +1,230 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2022 * \copyright GPLv2 or later */ #include "HEJ/jets.hh" #include #include #include #include "HEJ/Constants.hh" // generated headers #include "HEJ/currents/j_j.hh" #include "HEJ/currents/j_jqqbar.hh" #include "HEJ/currents/j_qqbar_j.hh" #include "HEJ/currents/juno_j.hh" namespace { using COM = std::complex; } // Anonymous Namespace namespace HEJ { namespace currents { - - // Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113 - // @TODO: code multiplication with MatrixElement.cc - double K_g(double p1minus, double paminus) { - return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A; - } - double K_g( - HLV const & pout, - HLV const & pin - ) { - if(pin.z() > 0) return K_g(pout.plus(), pin.plus()); - return K_g(pout.minus(), pin.minus()); - } - double ME_qq( HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ using helicity::plus; using helicity::minus; // we only have to distinguish between the two possibilities of // contracting same-helicity or opposite-helicity currents. COM const Mp = HEJ::j_j(p1in, p1out, p2in, p2out); COM const Mm = HEJ::j_j(p1in, p1out, p2in, p2out); double const sst = std::norm(Mm) + std::norm(Mp); // Multiply by factor 2 (helicities) return 2.*sst; } namespace { template double amp_juno_j( HLV const & pa, HLV const & pb, HLV const & pg, HLV const & p1, HLV const & p2 ) { // TODO: code duplication with Wjets const COM u1 = U1_j(pa,p1,pb,p2,pg); const COM u2 = U2_j(pa,p1,pb,p2,pg); const COM l = L_j(pa,p1,pb,p2,pg); 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_unob_qq( HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in ){ using helicity::minus; using helicity::plus; // only calculate half of the helicity amplitudes, // the remaining ones follow from CP symmetry const double amm = amp_juno_j(p1in, p2in, pg, p1out, p2out); const double amp = amp_juno_j(p1in, p2in, pg, p1out, p2out); const double apm = amp_juno_j< plus, minus>(p1in, p2in, pg, p1out, p2out); const double app = amp_juno_j< plus, plus>(p1in, p2in, pg, p1out, p2out); constexpr double hel_fac = 2.; return hel_fac * (amm + amp + apm + app); } namespace { // helicity amplitudes for j jqqbar contraction template double amp_j_jqqbar( HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2, HLV const & p3 ) { // TODO: code duplication with Wjets.cc // TODO: symbolic colour factors 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 = j_qggm1(pb,p2,p3,pa,p1); const COM m2 = j_qggm2(pb,p2,p3,pa,p1); const COM m3 = j_qggm3(pb,p2,p3,pa,p1); 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 double ME_qqbar_qg( HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2, HLV const & p3 ){ using helicity::minus; using helicity::plus; const double Mmmm = amp_j_jqqbar(pa, pb, p1, p2, p3); const double Mmmp = amp_j_jqqbar(pa, pb, p1, p2, p3); const double Mpmm = amp_j_jqqbar< plus, minus>(pa, pb, p1, p2, p3); const double Mpmp = amp_j_jqqbar< plus, plus>(pa, pb, p1, p2, p3); // Factor of 2 for the helicity for conjugate configurations return 2.*(Mmmm+Mmmp+Mpmm+Mpmp)/3./C_F; } namespace { // helicity amplitudes for matrix element with central qqbar // TODO: update to arxiv:2012.10310, around eq. (3.33) template double amp_Cenqqbar_qq( HLV const & pa, HLV const & p1, HLV const & pb, HLV const & p4, HLV const & pq, HLV const & pqbar, HLV const & q11, HLV const & q12 ){ using std::norm; const COM sym = M_sym( pa, p1, pb, p4, pq, pqbar, q11, q12 ); const COM cross = M_cross( pa, p1, pb, p4, pq, pqbar, q11, q12 ); const COM uncross = M_qbar( pa, p1, pb, p4, pq, pqbar, 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 double ME_Cenqqbar_qq( HLV const & pa, HLV const & pb, std::vector const & partons, const bool qqbarmarker, const std::size_t nabove ){ using helicity::plus; using helicity::minus; CLHEP::HepLorentzVector q1 = pa; for(std::size_t 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 = partons[nabove+1]; HLV pqbar = partons[nabove+2]; if(qqbarmarker) std::swap(pq, pqbar); const HLV p1 = partons.front(); const HLV pn = partons.back(); // 8 helicity choices, but only 4 independent ones //(complex conjugation related). // In principle, the proper helicity labeling depends on // whether we have antiquark lines (aqlinea and aqlineb). // However, this only corresponds to a relabeling, // which doesn't change the sum over all helicities below. const double amp_mm = amp_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second ); const double amp_mp = amp_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second ); const double amp_pm = amp_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second ); const double amp_pp = amp_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, qq.first, -qq.second ); //Result (averaged, without coupling or t-channel props). Factor of //2 for the 4 helicity configurations I didn't work out explicitly const HLV q3 = q1 - pq - pqbar; constexpr double hel_fac = 2.; // TODO: explain the 1/64 return hel_fac/64.*(amp_mm+amp_mp+amp_pm+amp_pp) / (q1.m2()*q3.m2()); } } // namespace currents } // namespace HEJ