Page MenuHomeHEPForge
Files F10882039

No OneTemporary
Actions

View Options

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 <array>
#include <cassert>
#include <cmath>
#include <complex>
#include <limits>
#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<double>;
constexpr double infinity = std::numeric_limits<double>::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<std::complex<double>> result(3);
static auto ql_B0 = [](){
ql::Bubble<std::complex<double>,double,double> ql_B0;
ql_B0.setCacheSize(100);
return ql_B0;
}();
static std::vector<double> masses(2);
static std::vector<double> 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<std::complex<double>> result(3);
static auto ql_C0 = [](){
ql::Triangle<std::complex<double>,double,double> ql_C0;
ql_C0.setCacheSize(100);
return ql_C0;
}();
static std::vector<double> masses(3);
static std::vector<double> 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<COM, 2> 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<COM, 2> res = {qH1.dot(qH2), 1.};
for(auto & tt: res) tt /= (3.*M_PI*vev);
return res;
}
std::array<COM, 2> 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<minus, minus>(
p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
);
const COM amp_mp = HEJ::j_h_j<minus, plus>(
p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
);
const COM amp_pm = HEJ::j_h_j<plus, minus>(
p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
);
const COM amp_pp = HEJ::j_h_j<plus, plus>(
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<minus, minus>(
pa, pb, pn, ph1, ph2, T_pref[0], T_pref[1]
);
const COM amp_mp = HEJ::jgh_j<minus, plus>(
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<Helicity h1, Helicity h2, Helicity hg>
double amp_juno_jgh(
HLV const & pg, HLV const & p1, HLV const & pa,
HLV const & ph1, HLV const & ph2, HLV const & pb,
std::array<COM, 2> const & T_pref
) {
// TODO: code duplication with Wjets and pure jets
const COM u1 = U1_jgh<h1, h2, hg>(pg, p1, pa, ph1, ph2, pb, T_pref[0], T_pref[1]);
const COM u2 = U2_jgh<h1, h2, hg>(pg, p1, pa, ph1, ph2, pb, T_pref[0], T_pref[1]);
const COM l = L_jgh<h1, h2, hg>(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<double, 2, 2, 2> amp;
// NOLINTNEXTLINE
#define ASSIGN_HEL(RES, J, H1, H2, HG) \
RES[H1][H2][HG] = J<H1, H2, HG>( \
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<Helicity h1, Helicity h2, Helicity hg>
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<COM, 2> const & T_pref
) {
// TODO: code duplication with Wjets and pure jets
const COM u1 = U1_h_j<h1, h2, hg>(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]);
const COM u2 = U2_h_j<h1, h2, hg>(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]);
const COM l = L_h_j<h1, h2, hg>(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<double, 2, 2, 2> amp{};
-// NOLINTNEXTLINE
-#define ASSIGN_HEL(RES, J, H1, H2, HG) \
- RES[H1][H2][HG] = J<H1, H2, HG>( \
- 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<H1, H2, HG>( \
+ 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

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 7:01 AM (6 h, 35 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983255
Default Alt Text
(15 KB)

Event Timeline