Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7879845
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
100 KB
Subscribers
None
View Options
diff --git a/src/Hjets.cc b/src/Hjets.cc
index c0d5ed6..d70b079 100644
--- a/src/Hjets.cc
+++ b/src/Hjets.cc
@@ -1,387 +1,389 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2022
* \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;
}
} // 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)
- ;
+
+ 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<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);
#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.;
+ return ampsq;
}
//@}
} // namespace currents
} // namespace HEJ
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 2d8b8d4..ab2af6e 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,2439 +1,2439 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2022
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <iterator>
#include <limits>
#include <unordered_map>
#include <utility>
#include "fastjet/PseudoJet.hh"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/Event.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "HEJ/WWjets.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Zjets.hh"
#include "HEJ/event_types.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/jets.hh"
#include "HEJ/utility.hh"
namespace HEJ {
namespace {
// Colour acceleration multiplier for gluons
// see eq:K_g in developer manual
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(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
// Colour acceleration multiplier for quarks
// see eq:K_q in developer manual
constexpr double K_q = C_F;
// Colour acceleration multipliers
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
){
if(type == pid::gluon) return K_g(pout, pin);
return K_q;
}
} // namespace
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*std::log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
return (
1. + alpha_s/(4.*M_PI)*BETA0*std::log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
std::vector <Weights> virtual_part=virtual_corrections(event);
if(tree_kin_part.size() != virtual_part.size()) {
throw std::logic_error("tree and virtuals have different sizes");
}
Weights sum = Weights{0., std::vector<double>(event.variations().size(), 0.)};
for(size_t i=0; i<tree_kin_part.size(); ++i) {
sum += tree_kin_part.at(i)*virtual_part.at(i);
}
return tree_param(event)*sum;
}
Weights MatrixElement::tree(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
double sum = 0.;
for(double i : tree_kin_part) {
sum += i;
}
return tree_param(event)*sum;
}
Weights MatrixElement::tree_param(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
std::vector<Weights> MatrixElement::virtual_corrections(Event const & event) const {
if(!event.valid_hej_state(param_.soft_pt_regulator)) {
return {Weights{0., std::vector<double>(event.variations().size(), 0.)}};
}
// only compute once for each renormalisation scale
std::unordered_map<double, std::vector<double> > known_vec;
std::vector<double> central_vec=virtual_corrections(event, event.central().mur);
known_vec.emplace(event.central().mur, central_vec);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
if(ME_it == end(known_vec)) {
known_vec.emplace(var.mur, virtual_corrections(event, var.mur));
}
}
// At this stage known_vec contains one vector of virtual corrections for each mur value
// Now put this into a vector of Weights
std::vector<Weights> result_vec;
for(size_t i=0; i<central_vec.size(); ++i) {
Weights result;
result.central = central_vec.at(i);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
result.variations.emplace_back(ME_it->second.at(i));
}
result_vec.emplace_back(result);
}
return result_vec;
}
template<class InputIterator>
std::vector <double> MatrixElement::virtual_corrections_interference(
InputIterator begin_parton, InputIterator end_parton,
fastjet::PseudoJet const & q0_t,
fastjet::PseudoJet const & q0_b,
const double mur
) const{
const double alpha_s = alpha_s_(mur);
auto qi_t = q0_t;
auto qi_b = q0_b;
double sum_top = 0.;
double sum_bot = 0.;
double sum_mix = 0.;
for(auto parton_it = begin_parton; parton_it != end_parton; ++parton_it){
Particle parton = *parton_it;
Particle parton_next = *(parton_it+1);
const double dy = parton_next.rapidity() - parton.rapidity();
const double tmp_top = omega0(alpha_s, mur, qi_t)*dy;
const double tmp_bot = omega0(alpha_s, mur, qi_b)*dy;
sum_top += tmp_top;
sum_bot += tmp_bot;
sum_mix += (tmp_top + tmp_bot) / 2.;
qi_t -= parton_next.p;
qi_b -= parton_next.p;
}
if (param_.nlo.enabled){
return {(sum_top), (sum_bot), (sum_mix)};
}
return {exp(sum_top), exp(sum_bot), exp(sum_mix)};
}
double MatrixElement::virtual_corrections_W(
Event const & event,
const double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
std::size_t first_idx = 0;
std::size_t last_idx = partons.size() - 1;
#ifndef NDEBUG
bool wc = true;
#endif
bool wqq = false;
// With extremal qqbar or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else if (event.type() == event_type::qqbar_exb) {
q -= partons[1].p;
++first_idx;
if (std::abs(partons[0].type) != std::abs(partons[1].type)){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else {
if(event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf){
--last_idx;
}
if (in[0].type != partons[0].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
std::size_t first_idx_qqbar = last_idx;
std::size_t last_idx_qqbar = last_idx;
//if qqbarMid event, virtual correction do not occur between qqbar pair.
if(event.type() == event_type::qqbar_mid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(std::abs(backquark->type) != std::abs((backquark+1)->type)) {
wqq=true;
#ifndef NDEBUG
wc=false;
#endif
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqbar = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(std::size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqbar) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
#ifndef NDEBUG
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf
);
#endif
if (param_.nlo.enabled){
double nlo_virtual =1.;
//Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) nlo_virtual+=exponent;
return nlo_virtual;
}
return std::exp(exponent);
}
std::vector <double> MatrixElement::virtual_corrections_WW(
Event const & event,
const double mur
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
assert(event.decays().size() == 2);
std::vector<fastjet::PseudoJet> plbar;
std::vector<fastjet::PseudoJet> pl;
for(auto const & decay_pair : event.decays()) {
auto const decay = decay_pair.second;
if(decay.at(0).type < 0) {
plbar.emplace_back(decay.at(0).p);
pl .emplace_back(decay.at(1).p);
}
else {
pl .emplace_back(decay.at(0).p);
plbar.emplace_back(decay.at(1).p);
}
}
fastjet::PseudoJet q_t = pa - partons[0].p - pl[0] - plbar[0];
fastjet::PseudoJet q_b = pa - partons[0].p - pl[1] - plbar[1];
auto const begin_parton = cbegin(partons);
auto const end_parton = cend(partons) - 1;
if (param_.nlo.enabled){
std::vector<double> virt_corrections_nlo {1.0,1.0,1.0}; //set virtual corrections to 1.
// Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) {
std::vector <double> virt_corrections_nlo_interference =
virtual_corrections_interference(
begin_parton, end_parton, q_t, q_b, mur
);
assert(
virt_corrections_nlo_interference.size()
== virt_corrections_nlo.size()
);
for(std::size_t i = 0; i < virt_corrections_nlo.size(); ++i){
virt_corrections_nlo[i] += virt_corrections_nlo_interference[i];
}
}
return virt_corrections_nlo;
}
return virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur);
}
std::vector <double> MatrixElement::virtual_corrections_Z_qq(
Event const & event,
const double mur,
Particle const & ZBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p;
fastjet::PseudoJet q_b = pa - partons[0].p;
auto begin_parton = cbegin(partons);
auto end_parton = cend(partons) - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q_t -= partons[1].p;
q_b -= partons[1].p;
++begin_parton;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--end_parton;
}
if (param_.nlo.enabled){
//set virtual corrections to 1.
std::vector<double> virt_corrections_nlo {1.0,1.0,1.0};
// Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) {
std::vector <double> virt_corrections_nlo_interference =
virtual_corrections_interference(
begin_parton, end_parton, q_t, q_b, mur
);
assert(
virt_corrections_nlo.size()
== virt_corrections_nlo_interference.size()
);
for(std::size_t i = 0; i < virt_corrections_nlo.size(); ++i){
virt_corrections_nlo[i] += virt_corrections_nlo_interference[i];
}
}
return virt_corrections_nlo;
}
return virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur);
}
double MatrixElement::virtual_corrections_Z_qg(
Event const & event,
const double mur,
Particle const & ZBoson,
const bool is_gq_event
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
// If this is a gq event, don't subtract the Z momentum from first q
fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p);
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity()
- partons[j].rapidity());
q -= partons[j+1].p;
}
if (param_.nlo.enabled){
double nlo_virtual =1.;
//Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) nlo_virtual+=sum;
return nlo_virtual;
}
return exp(sum);
}
std::vector<double> MatrixElement::virtual_corrections(
Event const & event,
const double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
std::vector<Particle> bosons = filter_AWZH_bosons(out);
if(event.jets().size() != param_.nlo.nj && param_.nlo.enabled) {
throw std::logic_error{
"Input event has number of jets different to stated NLO "
"input in config file: " + std::to_string(event.jets().size())
+ " vs " +std::to_string(param_.nlo.nj) + "\n"
};
}
if(bosons.size() > 2) {
throw not_implemented("Emission of >2 bosons is unsupported");
}
if(bosons.size() == 2) {
if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp) {
return virtual_corrections_WW(event, mur);
}
else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm) {
return virtual_corrections_WW(event, mur);
}
throw not_implemented("Emission of bosons of unsupported type");
}
if(bosons.size() == 1) {
const auto AWZH_boson = bosons[0];
if(std::abs(AWZH_boson.type) == pid::Wp){
return {virtual_corrections_W(event, mur, AWZH_boson)};
}
if(AWZH_boson.type == pid::Z_photon_mix){
if(is_gluon(in.back().type)){
// This is a qg event
return {virtual_corrections_Z_qg(event, mur, AWZH_boson, false)};
}
if(is_gluon(in.front().type)){
// This is a gq event
return {virtual_corrections_Z_qg(event, mur, AWZH_boson, true)};
}
// This is a qq event
return virtual_corrections_Z_qq(event, mur, AWZH_boson);
}
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
std::size_t first_idx = 0;
std::size_t last_idx = out.size() - 1;
// if there is a Higgs boson _not_ emitted off an incoming gluon,
// extremal qqbar or unordered gluon outside the extremal partons
// then it is not part of the FKL ladder
// and does not contribute to the virtual corrections
if((out.front().type == pid::Higgs && in.front().type != pid::gluon)
|| event.type() == event_type::unob
|| event.type() == event_type::qqbar_exb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs && in.back().type != pid::gluon)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf){
--last_idx;
}
std::size_t first_idx_qqbar = last_idx;
std::size_t last_idx_qqbar = last_idx;
//if central qqbar event, virtual correction do not occur between q-qbar.
if(event.type() == event_type::qqbar_mid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return {0.};
last_idx = std::distance(begin(out), backquark);
first_idx_qqbar = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(std::size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqbar) q -= out[last_idx+1].p;
for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| (out.back().type == pid::Higgs && in.back().type != pid::gluon)
|| event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf
);
if (param_.nlo.enabled){
const auto partons = filter_partons(event.outgoing());
double nlo_virtual =1.;
//Only apply virtual corrections to a nlo order event.
if (partons.size() == param_.nlo.nj) nlo_virtual+=exponent;
return {nlo_virtual};
}
return {std::exp(exponent)};
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return CL;
}
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, p1, p2);
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, p1, p2);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, p1, p2);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda){
return Cls;
}
return Cls - 4.0 / (kperp * kperp);
}
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
){
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return CL;
}
//! Lipatov vertex
double C2Lipatov( // B
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim,
CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom,
CLHEP::HepLorentzVector const & pop
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, pim, pip, pom, pop);
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, pim, pip, pom, pop);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, pim, pip, pom, pop);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, pa, pb, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda) {
return Cls;
}
return Cls - 4.0 / (kperp * kperp);
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pg Unordered gluon momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_uno_current(
[[maybe_unused]] ParticleID aptype,
ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
assert(aptype!=pid::gluon); // aptype cannot be gluon
const double t1 = (pa - p1 - pg).m2();
const double t2 = (pb - pn).m2();
return K_q * K(bptype, pn, pb)*currents::ME_unob_qq(pg, p1, pa, pn, pb) / (t1 * t2);
}
/** Matrix element squared for tree-level current-current scattering
* @param bptype Particle b PDG ID
* @param pgin Incoming gluon momentum
* @param p1 More backward (anti-)quark from splitting Momentum
* @param p2 Less backward (anti-)quark from splitting Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The forward qqbar contribution can be calculated by reversing the
* argument ordering.
*/
double ME_qqbar_current(
ParticleID bptype,
CLHEP::HepLorentzVector const & pgin,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb
){
const double t1 = (pgin - p1 - p2).m2();
const double t2 = (pn - pb).m2();
return K(bptype, pn, pb)*currents::ME_qqbar_qg(pb, pgin, pn, p2, p1) / (t1 * t2);
}
/* \brief Matrix element squared for central qqbar tree-level current-current
* scattering
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqbarpair
* @param nbelow Number of gluons emitted after central qqbarpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering
*/
double ME_qqbar_mid_current(
ParticleID aptype, ParticleID bptype, int nabove,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons
){
using namespace currents;
// CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards)
const bool swap_qqbar=pqbar.rapidity() < pq.rapidity();
CLHEP::HepLorentzVector const & p1 = partons.front();
CLHEP::HepLorentzVector const & pn = partons.back();
const double t1 = (p1 - pa).m2();
const double t2 = (pb - pn).m2();
return K(aptype, p1, pa)
*K(bptype, pn, pb)
*ME_Cenqqbar_qq(
pa, pb, partons,
swap_qqbar, nabove
) / (t1 * t2);
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
const double t1 = (p1 - pa).m2();
const double t2 = (pb - pn).m2();
return K(aptype, p1, pa)
* K(bptype, pn, pb)
* currents::ME_qq(p1, pa, pn, pb)/(t1 * t2);
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc, ParticleProperties const & Wprop
){
using namespace currents;
// We know it cannot be gg incoming.
assert(!(aptype==pid::gluon && bptype==pid::gluon));
if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
}
if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop);
}
// they are both quark
if (wc){ // emission off b, (first argument pbout)
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
if (is_quark(aptype))
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
// emission off a, (first argument paout)
if (is_quark(aptype)) {
if (is_quark(bptype))
return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
// a is anti-quark
if (is_quark(bptype))
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_W_uno_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc, ParticleProperties const & Wprop
){
using namespace currents;
// we know they are not both gluons
assert(bptype != pid::gluon || aptype != pid::gluon);
if (bptype == pid::gluon && aptype != pid::gluon) {
// b gluon => W emission off a
if (is_quark(aptype))
return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
// they are both quark
if (wc) {// emission off b, i.e. b is first current
if (is_quark(bptype)){
if (is_quark(aptype))
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
if (is_quark(aptype))
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
// wc == false, emission off a, i.e. a is first current
if (is_quark(aptype)) {
if (is_quark(bptype)) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
//qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
// a is anti-quark
if (is_quark(bptype)) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
//qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
/** \brief Matrix element squared for backward qqbar tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param swap_qqbar Ordering of qqbar pair. False: pqbar extremal.
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqbarb Tree-Level Current-Current Scattering
*
* @note calculate forwards qqbar contribution by reversing argument ordering.
*/
double ME_W_qqbar_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const swap_qqbar, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==pid::gluon && bptype==pid::gluon) {
//a gluon, b gluon gg->qqbarWg
// This will be a wqqbar emission as there is no other possible W Emission
// Site.
if (swap_qqbar)
return ME_WExqqbar_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop);
return ME_WExqqbar_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop);
}
assert(aptype==pid::gluon && bptype!=pid::gluon );
//a gluon => W emission off b leg or qqbar
if (!wc){ // W Emitted from backwards qqbar
if (swap_qqbar)
return ME_WExqqbar_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop);
return ME_WExqqbar_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop);
}
// W Must be emitted from forwards leg.
return ME_W_Exqqbar_QQq(pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop);
}
/* \brief Matrix element squared for central qqbar tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqbarpair
* @param nbelow Number of gluons emitted after central qqbarpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqbar
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering
*/
double ME_W_qqbar_mid_current(
ParticleID aptype, ParticleID bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
// CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards)
const bool swap_qqbar=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;
if(wqq)
return wt*ME_WCenqqbar_qq(pa, pb, pl, plbar, partons,
is_antiquark(aptype),is_antiquark(bptype),
swap_qqbar, nabove, Wprop);
return wt*ME_W_Cenqqbar_qq(pa, pb, pl, plbar, partons,
is_antiquark(aptype), is_antiquark(bptype),
swap_qqbar, nabove, nbelow, wc, Wprop);
}
/** Matrix element squared for tree-level current-current scattering With Z+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
std::vector<double> ME_Z_current(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if(is_anyquark(aptype) && is_gluon(bptype)){
// This is a qg event
return { ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw) };
}
if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
return { ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw) };
}
assert(is_anyquark(aptype) && is_anyquark(bptype));
// This is a qq event
return ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With Z+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
std::vector<double> ME_Z_uno_current(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if (is_anyquark(aptype) && is_gluon(bptype)) {
// This is a qg event
return { ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw) };
}
if (is_gluon(aptype) && is_anyquark(bptype)) {
// This is a gq event
return { ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw) };
}
assert(is_anyquark(aptype) && is_anyquark(bptype));
// This is a qq event
return ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb, double vev
){
const double t1 = (pa - p1).m2();
const double tn = (pb - pn).m2();
return
K(aptype, p1, pa)
*K(bptype, pn, pb)
*currents::ME_H_qQ(
pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev
) / (t1 * tn * qH.m2() * qHp1.m2());
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param pg Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*
* @note This function assumes unordered gluon backwards from pa-p1 current.
* For unof, reverse call order
*/
double ME_Higgs_current_uno(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb, double vev
){
const double t1 = (pa - p1 - pg).m2();
const double tn = (pn - pb).m2();
return
K(aptype, p1, pa)
*K(bptype, pn, pb)
*currents::ME_H_unob_qQ(
pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev
) / (t1 * qH.m2() * qHp1.m2() * tn);
}
/** Matrix element squared for tree-level scattering with WW+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pl1bar Particle l1bar Momentum
* @param pl1 Particle l1 Momentum
* @param pl2bar Particle l2bar Momentum
* @param pl2 Particle l2 Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
std::vector <double> ME_WW_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pl1bar,
CLHEP::HepLorentzVector const & pl1,
CLHEP::HepLorentzVector const & pl2bar,
CLHEP::HepLorentzVector const & pl2,
ParticleProperties const & Wprop
){
using namespace currents;
if (aptype > 0 && bptype > 0)
return ME_WW_qQ(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
if (aptype < 0 && bptype > 0)
return ME_WW_qbarQ(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
if (aptype > 0 && bptype < 0)
return ME_WW_qQbar(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
if (aptype < 0 && bptype < 0)
return ME_WW_qbarQbar(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop);
throw std::logic_error("unreachable");
}
void validate(MatrixElementConfig const & config) {
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
std::vector<double> MatrixElement::tree_kin(
Event const & ev
) const {
if(!ev.valid_hej_state(param_.soft_pt_regulator)) return {0.};
if(!ev.valid_incoming()){
throw std::invalid_argument{
"Invalid momentum for one or more incoming particles. "
"Incoming momenta must have vanishing mass and transverse momentum."
};
}
std::vector<Particle> bosons = filter_AWZH_bosons(ev.outgoing());
if(bosons.empty()) {
return {tree_kin_jets(ev)};
}
if(bosons.size() == 1) {
switch(bosons[0].type){
case pid::Higgs:
return {tree_kin_Higgs(ev)};
case pid::Wp:
case pid::Wm:
return {tree_kin_W(ev)};
case pid::Z_photon_mix:
return tree_kin_Z(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
if(bosons.size() == 2) {
if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp){
return tree_kin_WW(ev);
}
else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm){
return tree_kin_WW(ev);
}
throw not_implemented("Emission of bosons of unsupported type");
}
throw not_implemented("Emission of >2 bosons is unsupported");
}
namespace {
constexpr int EXTREMAL_JET_IDX = 1;
constexpr int NO_EXTREMAL_JET_IDX = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == EXTREMAL_JET_IDX;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
template<class InputIterator>
std::vector <double> FKL_ladder_weight_mix(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0_t, CLHEP::HepLorentzVector const & q0_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
const double lambda
){
double wt_top = 1;
double wt_bot = 1;
double wt_mix = 1;
auto qi_t = q0_t;
auto qi_b = q0_b;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1_t = qi_t - g;
const auto qip1_b = qi_b - g;
if(treat_as_extremal(*gluon_it)){
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, lambda)*C_A;
} else{
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, p1, pn, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
}
qi_t = qip1_t;
qi_b = qip1_b;
}
return {wt_top, wt_bot, wt_mix};
}
std::vector<Particle> tag_extremal_jet_partons( Event const & ev ){
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(NO_EXTREMAL_JET_IDX);
}
return out_partons;
}
auto const & jets = ev.jets();
std::vector<fastjet::PseudoJet> extremal_jets;
if(! is_backward_g_to_h(ev)) {
auto most_backward = begin(jets);
// skip jets caused by unordered emission or qqbar
if(ev.type() == event_type::unob || ev.type() == event_type::qqbar_exb){
assert(jets.size() >= 2);
++most_backward;
}
extremal_jets.emplace_back(*most_backward);
}
if(! is_forward_g_to_h(ev)) {
auto most_forward = end(jets) - 1;
if(ev.type() == event_type::unof || ev.type() == event_type::qqbar_exf){
assert(jets.size() >= 2);
--most_forward;
}
extremal_jets.emplace_back(*most_forward);
}
const auto extremal_jet_indices = ev.particle_jet_indices(
extremal_jets
);
assert(extremal_jet_indices.size() == out_partons.size());
for(std::size_t i = 0; i < out_partons.size(); ++i){
assert(is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
EXTREMAL_JET_IDX:
NO_EXTREMAL_JET_IDX;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
double tree_kin_jets_qqbar_mid(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
double lambda
){
CLHEP::HepLorentzVector pq;
CLHEP::HepLorentzVector pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqbar
auto qqbart = q0;
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqbart -= to_HepLorentzVector(*parton_it);
}
const int nabove = std::distance(begin_ladder, backmidquark);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_qqbar_mid_current(
aptype, bptype, nabove, pa, pb,
pq, pqbar, partonsHLV
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqbart, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_qqbar(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const auto pgin = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto p1 = to_HepLorentzVector(*BeginPart);
const auto p2 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon.
const double current_factor = ME_qqbar_current(
(EndIn-1)->type, pgin, p1, p2, pn, pb
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pgin-p1-p2, pgin, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const double current_factor = ME_uno_current(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pa-p1-pg, pa, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
} // namespace
double MatrixElement::tree_kin_jets(Event const & ev) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if (ev.type()==event_type::FKL){
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn,
param_.regulator_lambda
);
}
if (ev.type()==event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
if (ev.type()==event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
if (ev.type()==event_type::extremal_qqbar_backward){
return tree_kin_jets_qqbar(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
if (ev.type()==event_type::extremal_qqbar_forward){
return tree_kin_jets_qqbar(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
if (ev.type()==event_type::central_qqbar){
return tree_kin_jets_qqbar_mid(incoming[0].type, incoming[1].type,
to_HepLorentzVector(incoming[0]),
to_HepLorentzVector(incoming[1]),
partons, param_.regulator_lambda);
}
throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets");
}
namespace {
double tree_kin_W_FKL(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
double lambda, ParticleProperties const & Wprop
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_uno(InIter BeginIn, partIter BeginPart,
partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w
auto q0 = pa - p1 - pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_uno_current(
(BeginIn)->type, (BeginIn+1)->type, pn, pb,
p1, pa, pg, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_qqbar(InIter BeginIn, partIter BeginPart,
partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const bool swap_qqbar=is_quark(*BeginPart);
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_qqbar?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_qqbar?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqbar_current(
(BeginIn)->type, (BeginIn+1)->type, pa, pb,
pq, pqbar, pn, plbar, pl, swap_qqbar, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
double tree_kin_W_qqbar_mid(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
double lambda, ParticleProperties const & Wprop
){
CLHEP::HepLorentzVector pq;
CLHEP::HepLorentzVector pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons.front());
auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
// t-channel momentum after qqbar
auto qqbart = q0;
bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqbar emit W
bool wc = !wqq && (aptype==partons.front().type); //leg b emits w
assert(!wqq || !wc);
if(wqq){ // emission from qqbar
qqbart -= pl + plbar;
} else if(!wc) { // emission from leg a
q0 -= pl + plbar;
qqbart -= pl + plbar;
}
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqbart -= to_HepLorentzVector(*parton_it);
}
const int nabove = std::distance(begin_ladder, backmidquark);
const int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqbar_mid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqbart, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
} // namespace
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
#ifndef NDEBUG
// assert that there is exactly one decay corresponding to the W
assert(ev.decays().size() == 1);
auto const & w_boson{
std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(),
[] (Particle const & p) -> bool {
return std::abs(p.type) == ParticleID::Wp;
}) };
assert(w_boson != ev.outgoing().cend());
assert( static_cast<long int>(ev.decays().cbegin()->first)
== std::distance(ev.outgoing().cbegin(), w_boson) );
#endif
// find decay products of W
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) )
|| ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) );
// get lepton & neutrino
CLHEP::HepLorentzVector plbar;
CLHEP::HepLorentzVector pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.at(1));
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == FKL){
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_backward){
return tree_kin_W_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_forward){
return tree_kin_W_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqbar_backward){
return tree_kin_W_qqbar(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqbar_forward){
return tree_kin_W_qqbar(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
assert(ev.type() == central_qqbar);
return tree_kin_W_qqbar_mid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
namespace /* WW */ {
std::vector <double> tree_kin_WW_FKL(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & pl1bar, CLHEP::HepLorentzVector const & pl1,
CLHEP::HepLorentzVector const & pl2bar, CLHEP::HepLorentzVector const & pl2,
double lambda, ParticleProperties const & Wprop
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const std::vector <double> current_factor = ME_WW_current(
aptype, bptype, pn, pb, p1, pa,
pl1bar, pl1, pl2bar, pl2,
Wprop
);
auto const begin_ladder = cbegin(partons) + 1;
auto const end_ladder = cend(partons) - 1;
// pa -> W1 p1, pb -> W2 + pn
const auto q0 = pa - p1 - (pl1 + pl1bar);
// pa -> W2 p1, pb -> W1 + pn
const auto q1 = pa - p1 - (pl2 + pl2bar);
const std::vector <double> ladder_factor = FKL_ladder_weight_mix(
begin_ladder, end_ladder,
q0, q1, pa, pb, p1, pn,
lambda
);
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
} // namespace
std::vector <double> MatrixElement::tree_kin_WW(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const pa = to_HepLorentzVector(incoming[0]);
auto const pb = to_HepLorentzVector(incoming[1]);
auto const partons = tag_extremal_jet_partons(ev);
// W1 & W2
assert(ev.decays().size() == 2);
std::vector<CLHEP::HepLorentzVector> plbar;
std::vector<CLHEP::HepLorentzVector> pl;
for(auto const & decay_pair : ev.decays()) {
auto const decay = decay_pair.second;
// TODO: how to label W1, W2
if(decay.at(0).type < 0) {
plbar.emplace_back(to_HepLorentzVector(decay.at(0)));
pl .emplace_back(to_HepLorentzVector(decay.at(1)));
}
else {
pl .emplace_back(to_HepLorentzVector(decay.at(0)));
plbar.emplace_back(to_HepLorentzVector(decay.at(1)));
}
}
if(ev.type() == FKL) {
return tree_kin_WW_FKL(
incoming[0].type, incoming[1].type,
pa, pb, partons,
plbar[0], pl[0], plbar[1], pl[1],
param_.regulator_lambda,
param_.ew_parameters.Wprop()
);
}
throw std::logic_error("Can only reweight FKL events in WW");
}
namespace{
std::vector <double> tree_kin_Z_FKL(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw
){
const auto p1 = to_HepLorentzVector(partons[0]);
const auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
const std::vector <double> current_factor = ME_Z_current(
aptype, bptype, pn, pb, p1, pa,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
} else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-p1;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
} else {
// This is a qq event
const auto q0 = pa-p1-plbar-pl;
const auto q1 = pa-p1;
ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
template<class InIter, class partIter>
std::vector <double> tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const ParticleID aptype = (BeginIn)->type;
const ParticleID bptype = (BeginIn+1)->type;
const std::vector <double> current_factor = ME_Z_uno_current(
aptype, bptype, pn, pb, p1, pa, pg,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-pg-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-pg-p1;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q0 = pa-pg-p1-plbar-pl;
const auto q1 = pa-pg-p1;
ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
} // namespace
std::vector<double> MatrixElement::tree_kin_Z(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
// find decay products of Z
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0))
&& decay.at(0).type==-decay.at(1).type);
// get leptons
CLHEP::HepLorentzVector plbar;
CLHEP::HepLorentzVector pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.at(1));
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
const double stw2 = param_.ew_parameters.sin2_tw();
const double ctw = param_.ew_parameters.cos_tw();
if(ev.type() == FKL){
return tree_kin_Z_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_backward){
return tree_kin_Z_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_forward){
return tree_kin_Z_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets");
}
double MatrixElement::tree_kin_Higgs(Event const & ev) const {
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
// kinetic matrix element square for backward Higgs emission
// cf. eq:ME_h_jets_peripheral in developer manual,
// but without factors \alpha_s and the FKL ladder
double MatrixElement::MH2_backwardH(
const ParticleID type_forward,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pH,
CLHEP::HepLorentzVector const & pn
) const {
using namespace currents;
const double vev = param_.ew_parameters.vev();
return K(type_forward, pn, pb)*ME_H_gq(
pH, pa, pn, pb,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(4.*(N_C*N_C - 1)*(pa-pH).m2()*(pb-pn).m2());
}
// kinetic matrix element square for backward unordered emission
// and forward g -> Higgs
double MatrixElement::MH2_unob_forwardH(
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pH
) const {
using namespace currents;
const double vev = param_.ew_parameters.vev();
constexpr double K_f1 = K_q;
constexpr double nhel = 4.;
return K_f1*ME_juno_jgH(
pg, p1, pa, pH, pb,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(nhel*(N_C*N_C - 1)*(pa - p1 - pg).m2()*(pb - pH).m2());
}
double MatrixElement::tree_kin_Higgs_first(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(is_anyquark(incoming.front())) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming.front());
const auto pb = to_HepLorentzVector(incoming.back());
const auto pH = to_HepLorentzVector(outgoing.front());
const auto end_ladder = end(partons) - ((ev.type() == event_type::unof)?2:1);
const auto pn = to_HepLorentzVector(*end_ladder);
const double ladder = FKL_ladder_weight(
begin(partons), end_ladder,
pa - pH, pa, pb, pa, pn,
param_.regulator_lambda
);
if(ev.type() == event_type::unof) {
const auto pg = to_HepLorentzVector(outgoing.back());
return MH2_unob_forwardH(
pb, pa, pg, pn, pH
)*ladder;
}
return MH2_backwardH(
incoming.back().type,
pa, pb, pH, pn
)*ladder;
}
double MatrixElement::tree_kin_Higgs_last(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(is_anyquark(incoming.back())) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming.front());
const auto pb = to_HepLorentzVector(incoming.back());
const auto pH = to_HepLorentzVector(outgoing.back());
auto begin_ladder = begin(partons) + 1;
auto q0 = pa - to_HepLorentzVector(partons.front());
if(ev.type() == event_type::unob) {
q0 -= to_HepLorentzVector(*begin_ladder);
++begin_ladder;
}
const auto p1 = to_HepLorentzVector(*(begin_ladder - 1));
const double ladder = FKL_ladder_weight(
begin_ladder, end(partons),
q0, pa, pb, p1, pb,
param_.regulator_lambda
);
if(ev.type() == event_type::unob) {
const auto pg = to_HepLorentzVector(outgoing.front());
return MH2_unob_forwardH(
pa, pb, pg, p1, pH
)*ladder;
}
return MH2_backwardH(
incoming.front().type,
pb, pa, pH, p1
)*ladder;
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart,
CLHEP::HepLorentzVector const & qH,
CLHEP::HepLorentzVector const & qHp1,
double mt, bool inc_bot, double mb, double vev
){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb, vev
);
}
} // namespace
double MatrixElement::tree_kin_Higgs_between(Event const & ev) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor = NAN;
if(ev.type() == FKL){
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
}
else if(ev.type() == unob){
- current_factor = 4*C_A*C_A/2*tree_kin_Higgs_uno(
+ current_factor = tree_kin_Higgs_uno(
begin(incoming), end(incoming), begin(partons),
end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
- current_factor = 4*C_A*C_A/2*tree_kin_Higgs_uno(
+ current_factor = tree_kin_Higgs_uno(
rbegin(incoming), rend(incoming), rbegin(partons),
rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
throw std::logic_error("Can only reweight FKL or uno processes in H+Jets");
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor/(4.*(N_C*N_C-1.))*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) {
std::vector<Particle> bosons = filter_AWZH_bosons(ev.outgoing());
if(bosons.empty()) {
return 1.;
}
if(bosons.size() == 1) {
switch(bosons[0].type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return alpha_w*alpha_w;
case pid::Z_photon_mix:
return alpha_w*alpha_w;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
if(bosons.size() == 2) {
if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp) {
return alpha_w*alpha_w*alpha_w*alpha_w;
}
else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm) {
return alpha_w*alpha_w*alpha_w*alpha_w;
}
throw not_implemented("Emission of bosons of unsupported type");
}
throw not_implemented("Emission of >2 bosons is unsupported");
}
} // namespace
double MatrixElement::tree_param(Event const & ev, double mur) const {
assert(is_resummable(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*BETA0*std::log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 9:16 PM (1 d, 30 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806227
Default Alt Text
(100 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment