diff --git a/src/Hjets.cc b/src/Hjets.cc index 328b840..dff0705 100644 --- a/src/Hjets.cc +++ b/src/Hjets.cc @@ -1,430 +1,387 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \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; } - /** - * @brief Higgs+Jets FKL Contributions, function to handle all incoming types. - * @param p1out Outgoing Particle 1 - * @param p1in Incoming Particle 1 - * @param p2out Outgoing Particle 2 - * @param p2in Incoming Particle 2 - * @param qH1 t-channel momenta into higgs vertex - * @param qH2 t-channel momenta out of higgs vertex - * @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 j^\mu H j_\mu. FKL with higgs vertex somewhere in the FKL chain. - * Handles all possible incoming states. - */ - - // } - } // 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]); return 2.*C_F*std::real((l-u1)*std::conj(l+u2)) + 2.*C_F*C_F/3.*std::norm(u1+u2) ; } +} // namespace -//@{ -/** - * @brief Higgs+Jets Unordered Contributions, function to handle all incoming types. - * @param pg Unordered Gluon momenta - * @param p1out Outgoing Particle 1 - * @param p1in Incoming Particle 1 - * @param p2out Outgoing Particle 2 - * @param p2in Incoming Particle 2 - * @param qH1 t-channel momenta into higgs vertex - * @param qH2 t-channel momenta out of higgs vertex - * @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) - * - * Calculates j_{uno}^\mu H j_\mu. Unordered with higgs vertex - * somewhere in the FKL chain. Handles all possible incoming states. - */ - double juno_h_j( + 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 incBot, const double mb, const double vev + 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, incBot, mb, vev); + 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); + // 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 / 16.; } -} // 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, double mt, bool include_bottom, double mb, - double vev -){ - return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); -} //@} } // namespace currents } // namespace HEJ