diff --git a/current_generator/jphoton_juno.frm b/current_generator/jphoton_juno.frm
index 8fe7fb2..98339ba 100644
--- a/current_generator/jphoton_juno.frm
+++ b/current_generator/jphoton_juno.frm
@@ -1,74 +1,76 @@
*/**
* \brief Contraction of vector boson current (photon) with unordered current
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2025
* \copyright GPLv2 or later
*
* Current obtained from jphoton_j.frm and jV_juno.frm
*
*/
#include- include/helspin.frm
#include- include/write.frm
v pb,p2,pa,p1,pphoton,pr,pg,pr1,q2,q3;
i mu,nu;
s h;
#do h2={+,-}
* eq:juno in developer manual with pa -> pb, p1 -> pg
l [U1_jphoton `h2'] = (
+ Current(`h2'1, p2, nu, pg)*Current(`h2'1, pg, mu, pb)
+ 2*p2(nu)*Current(`h2'1, p2, mu, pb)
)/m2(p2 + pg);
l [U2_jphoton `h2'] = (
+ 2*Current(`h2'1, p2, mu, pb)*pb(nu)
- Current(`h2'1, p2, mu, pg)*Current(`h2'1, pg, nu, pb)
)/m2(pb - pg);
l [L_jphoton `h2'] = (
- (q2(nu) + q3(nu))*Current(`h2'1, p2, mu, pb)
- 2*pg(mu)*Current(`h2'1, p2, nu, pb)
+ 2*Current(`h2'1, p2, pg, pb)*d_(mu, nu)
+ m2(q2)*Current(`h2'1, p2, mu, pb)*(p1(nu)/m2(p1+pg) + pa(nu)/m2(pa+pg))
)/m2(q3);
#enddo
.sort
drop;
* multiply with polarisation vector and other currents
#do h1={+,-}
#do hphoton={+,-}
#do h2={+,-}
#do hg={+,-}
#do TENSOR={U1_jphoton,U2_jphoton,L_jphoton}
l [`TENSOR' `h1'`hphoton'`h2'`hg'] = (
[`TENSOR' `h2']
*Jgamma(`h1'1, `hphoton'1, mu, pa, p1, pphoton, pr)
*Eps(`hg'1, nu, pr1)
);
#enddo
#enddo
#enddo
#enddo
#enddo
*choice of auxiliary vector
id Eps(h?, mu?, pr1)*Current(h?, ?a) = Eps(h, mu, pg, p2)*Current(h, ?a);
also Eps(h?, mu?, pr1) = Eps(h, mu, pg, pb);
+id pa?.pb? = dot(pa, pb);
+
multiply replace_(
pr, p1,
q2, p2+pg-pb,
q3, p2-pb
);
.sort
#call ContractCurrents
.sort
format O4;
format c;
#call WriteHeader(`OUTPUT')
#call WriteOptimised(`OUTPUT',U1_jphoton,4,pa,pb,pphoton,p1,p2,pg)
#call WriteOptimised(`OUTPUT',U2_jphoton,4,pa,pb,pphoton,p1,p2,pg)
#call WriteOptimised(`OUTPUT',L_jphoton,4,pa,pb,pphoton,p1,p2,pg)
#call WriteFooter(`OUTPUT')
.end
diff --git a/current_generator/jphotonuno_j.frm b/current_generator/jphotonuno_j.frm
index e56ed5f..306a9f5 100644
--- a/current_generator/jphotonuno_j.frm
+++ b/current_generator/jphotonuno_j.frm
@@ -1,110 +1,110 @@
*/**
* \brief Contraction of vector boson (photon) unordered current with FKL current
*
* TODO: unify conventions with developer manual
* the current dictionary is as follows:
*
* code | manual
* pg | p_1
* p1 | p_2
* pa | p_a
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*
* Current contraction derived from jVuno_j.frm
*
*/
#include- include/helspin.frm
#include- include/write.frm
s h,s2g,sbg,taV1;
v p,p1,p2,pa,pb,pg,pV,pr,q1,q1g,q2,pr1,pr2;
i mu,nu,rho,sigma;
cf m2inv;
#do h1={+,-}
* eq:U1tensor in developer manual, up to factors 1/sij, 1/tij
l [U1_jphotonuno `h1'] = (
+ Current(`h1'1, p1, nu, p1+pg, mu, pa-pV, rho, pa)
+ Current(`h1'1, p1, nu, p1+pg, rho, p1+pg+pV, mu, pa)
+ Current(`h1'1, p1, rho, p1+pV, nu, p1+pg+pV, mu, pa)
);
* eq:U2tensor in developer manual, up to factors 1/sij, 1/tij
l [U2_jphotonuno `h1'] = (
+ Current(`h1'1, p1, mu, pa - pV - pg, nu, pa - pV, rho, pa)
+ Current(`h1'1, p1, mu, pa - pV - pg, rho, pa - pg, nu, pa)
+ Current(`h1'1, p1, rho, p1 + pV, mu, pa - pg, nu, pa)
);
* eq:Ltensor in developer manual, up to factors 1/sij, 1/tij
l [L_jphotonuno `h1'] = (
Current(`h1'1, p1, sigma, pa - pV, rho, pa) +
Current(`h1'1, p1, rho, p1 + pV, sigma, pa)
)*(
((pb(nu)/sbg + p2(nu)/s2g)*m2(q1g) + 2*q1(nu) - pg(nu))*d_(mu, sigma)
- 2*pg(mu)*d_(nu, sigma)
+ (2*pg(sigma) - q1(sigma))*d_(mu, nu)
)/taV1;
#enddo
.sort
* restore kinematic factors
id Current(h?, p1, mu?, q1?, nu?, q2?, rho?, pa) = (
Current(h, p1, mu, q1, nu, q2, rho, pa)*m2inv(q1)*m2inv(q2)
);
id Current(h?, p1, mu?, q1?, nu?, pa) = (
Current(h, p1, mu, q1, nu, pa)*m2inv(q1)
);
.sort
drop;
* multiply with polarisation vector and other currents
#do h1={+,-}
#do hp={+,-}
#do h2={+,-}
#do hg={+,-}
#do TENSOR={U1_jphotonuno,U2_jphotonuno,L_jphotonuno}
l [`TENSOR' `h1'`hp'`h2'`hg'] = (
[`TENSOR' `h1']
*Eps(`hg'1, nu, pr1)
*Current(`h2'1, p2, mu, pb)
*Eps(`hp'1, rho, pr2)
);
#enddo
#enddo
#enddo
#enddo
#enddo
* choice of best reference vector (p2 or pb)
id Eps(h?, nu?, pr1)*Current(h?, p2, mu?, pb) = Eps(h, nu, pg, p2)*Current(h, p2, mu, pb);
also Eps(h?, nu?, pr1) = Eps(h, nu, pg, pb);
* choice of best reference vector (p2 or pb)
-id Eps(h?, nu?, pr2)*Current(h?, p2, mu?, pb) = Eps(h, nu, pg, p2)*Current(h, p2, mu, pb);
-also Eps(h?, nu?, pr2) = Eps(h, nu, pg, pb);
+id Eps(h?, nu?, pr2)*Current(h?, p2, mu?, pb) = Eps(h, nu, pV, p1)*Current(h, p2, mu, pb);
+also Eps(h?, nu?, pr2) = Eps(h, nu, pV, pa);
multiply replace_(q1g,q1-pg);
multiply replace_(q1,pa-p1-pV);
.sort
#call ContractCurrents
multiply replace_(
s2g,m2(p2+pg),
sbg,m2(pb+pg),
taV1,m2(pa-pV-p1)
);
id m2inv(q1?) = 1/m2(q1);
.sort
format O4;
format c;
#call WriteHeader(`OUTPUT')
#call WriteOptimised(`OUTPUT',U1_jphotonuno,4,p1,p2,pa,pb,pg,pV)
#call WriteOptimised(`OUTPUT',U2_jphotonuno,4,p1,p2,pa,pb,pg,pV)
#call WriteOptimised(`OUTPUT',L_jphotonuno,4,p1,p2,pa,pb,pg,pV)
#call WriteFooter(`OUTPUT')
.end
\ No newline at end of file
diff --git a/include/HEJ/Photonjets.hh b/include/HEJ/Photonjets.hh
index 12dfa78..ff82555 100644
--- a/include/HEJ/Photonjets.hh
+++ b/include/HEJ/Photonjets.hh
@@ -1,85 +1,108 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020-2023
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in Photon+Jets.
*/
#pragma once
#include <vector>
#include "CLHEP/Vector/LorentzVector.h"
#include "HEJ/PDG_codes.hh"
namespace HEJ::currents {
using HLV = CLHEP::HepLorentzVector;
//! Square of qQ->qQphoton Photon+Jets Scattering Current
/**
* @param pa Momentum of initial state quark
* @param pb Momentum of initial state quark
* @param p1 Momentum of final state quark
* @param p2 Momentum of final state quark
* @param pphoton Momentum of final state photon
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @returns Square of the current contractions for qQ->qQphoton Scattering
*
* This returns the square of the current contractions in qQ->qQ scattering
* with an emission of a photon.
*/
std::vector <double> ME_photon_qQ(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pphoton,
ParticleID aptype, ParticleID bptype
);
//! Square of qg->qgphoton Photon+Jets Scattering Current
/**
* @param pa Momentum of initial state quark
* @param pb Momentum of initial state gluon
* @param p1 Momentum of final state quark
* @param p2 Momentum of final state gluon
* @param pphoton Momentum of final state photon
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @returns Square of the current contractions for qg->qgphoton Scattering
*
* This returns the square of the current contractions in qg->qg scattering
* with an emission of a photon.
*/
double ME_photon_qg(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pp,
ParticleID aptype, ParticleID bptype
);
//! Photon+Jets UNO Contributions
/**
* @brief Photon+Jets Unordered Contribution, unordered opposite to Photon side
* @param pa Incoming Particle 1 (photon emission)
* @param pb Incoming Particle 2 (unorderd emission)
* @param p1 Outgoing Particle 1 (photon emission)
* @param p2 Outgoing Particle 2 (unordered emission)
* @param pg Unordered Gluon
* @param pp Outgoing Photon
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
*
* @returns Square of the current contractions for for qQ -> photon qQ g
*
* This returns the square of the current contractions in qQ->qQg scattering
* with an emission of a photon.
*/
std::vector <double> ME_photon_uno_qQ(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pp, const HLV & pg,
ParticleID aptype, ParticleID bptype
);
+
+ /**
+ * @brief Photon+Jets Unordered Contribution, unordered same to Photon side
+ * @param pa Incoming Particle 1 (photon and unordered emission)
+ * @param pb Incoming Particle 2
+ * @param p1 Outgoing Particle 1 (photon emission and unordered)
+ * @param p2 Outgoing Particle 2
+ * @param pg Unordered Gluon
+ * @param pp Outgoing Photon
+ * @param aptype Initial particle 1 type
+ * @param bptype Initial particle 2 type
+ *
+ * @returns Square of the current contractions for for qg -> photon g qg
+ *
+ * This returns the square of the current contractions in qQ->qQg scattering
+ * with an emission of a photon.
+ */
+ double ME_photon_uno_qg(
+ const HLV & pa, const HLV & pb, const HLV & pg,
+ const HLV & p1, const HLV & p2, const HLV & pp,
+ const ParticleID aptype, const ParticleID /*bptype*/
+ );
+
} // namespace HEJ::currents
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 93c2d81..14fdce9 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,3252 +1,3261 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2025
* \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/Photonjets.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>
double MatrixElement::virtual_corrections_no_interference(
InputIterator begin_parton, InputIterator end_parton,
fastjet::PseudoJet & q,
const double mur
) const{
const double alpha_s = alpha_s_(mur);
double exponent = 0.;
for(auto parton_it = begin_parton; parton_it != end_parton; ++parton_it){
exponent += omega0(alpha_s, mur, q)*(
(parton_it+1)->rapidity() - parton_it->rapidity()
);
q -= (parton_it+1)->p;
}
return exponent;
}
template<class InputIterator>
std::vector <double> MatrixElement::virtual_corrections_interference(
InputIterator begin_parton, InputIterator end_parton,
fastjet::PseudoJet & q_t,
fastjet::PseudoJet & q_b,
const double mur
) const{
const double alpha_s = alpha_s_(mur);
double sum_top = 0.;
double sum_bot = 0.;
double sum_mix = 0.;
for(auto parton_it = begin_parton; parton_it != end_parton; ++parton_it){
const double dy = (parton_it+1)->rapidity() - parton_it->rapidity();
const double tmp_top = omega0(alpha_s, mur, q_t) * dy;
const double tmp_bot = omega0(alpha_s, mur, q_b) * dy;
sum_top += tmp_top;
sum_bot += tmp_bot;
sum_mix += (tmp_top + tmp_bot) / 2.;
q_t -= (parton_it+1)->p;
q_b -= (parton_it+1)->p;
}
return {sum_top, sum_bot, 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;
auto first_idx = cbegin(partons);
auto last_idx = cend(partons) - 1;
#ifndef NDEBUG
bool wc = true;
#endif
// 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
}
}
auto first_idx_qqbar = last_idx;
auto last_idx_qqbar = last_idx;
bool wqq = false;
//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); }
);
assert(backquark!=end(partons)-1 && (backquark+1)->type!=pid::gluon);
if(std::abs(backquark->type) != std::abs((backquark+1)->type)) {
wqq=true;
#ifndef NDEBUG
wc=false;
#endif
}
last_idx = backquark;
first_idx_qqbar = last_idx+1;
}
double exponent = virtual_corrections_no_interference(first_idx, last_idx, q, mur);
if(last_idx != first_idx_qqbar) {
q -= (last_idx+1)->p;
if(wqq) q -= WBoson.p;
exponent += virtual_corrections_no_interference(first_idx_qqbar, last_idx_qqbar, q, mur);
}
#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;
auto res = virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur);
if(param_.nlo.enabled) {
if(partons.size() > param_.nlo.nj) return std::vector(res.size(), 1.);
assert(partons.size() == param_.nlo.nj);
for(double & virt: res) virt += 1.;
} else {
for(double & virt: res) virt = exp(virt);
}
return res;
}
// Virtual corrections for Z or photon emission processes without
// interference (only one contribution)
double MatrixElement::virtual_corrections_Zphoton_single(
Event const & event,
const double mur,
Particle const & boson,
const bool emit_fwd
) 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 the boson is emitted from forward leg or central qqbar
// don't subtract the boson momentum from first q
fastjet::PseudoJet q;
if(emit_fwd || event.type()==event_type::central_qqbar) {
q = pa - partons[0].p;
} else {
q = pa - partons[0].p - boson.p;
}
auto first_idx = cbegin(partons);
auto last_idx = cend(partons) - 1;
// for uno/exqqbar the two extremal partons do not contribute to the virtual corrections
if (event.type() == event_type::unob
|| event.type() == event_type::qqbar_exb) {
// unordered gluon or first quark is partons[0] and is already subtracted
// partons[1] is the backward quark (uno case) or the second quark (exqqbar case)
q -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf) {
// end sum at next to last parton
--last_idx;
}
// if central qqbar event, virtual correction do not occur between qqbar pair
auto first_idx_qqbar = last_idx;
auto last_idx_qqbar = last_idx;
if(event.type() == event_type::central_qqbar){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
assert(backquark!=end(partons)-1 && (backquark+1)->type!=pid::gluon);
last_idx = backquark;
first_idx_qqbar = last_idx+1;
}
double exponent = virtual_corrections_no_interference(first_idx, last_idx, q, mur);
if(last_idx != first_idx_qqbar) {
q -= (last_idx+1)->p;
q -= boson.p;
exponent += virtual_corrections_no_interference(first_idx_qqbar, last_idx_qqbar, q, mur);
}
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 exp(exponent);
}
// Virtual corrections for Z or photon processes with 2 interfering
// contributions
// Returns 3 values (2 squares + 1 interference term)
std::vector<double> MatrixElement::virtual_corrections_Zphoton_double(
Event const & event,
const double mur,
Particle const & boson
) 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;
if(event.type() == event_type::central_qqbar && is_gluon(in.front().type)) {
// gq initiated central qqbar event -> "top" emission from qqbar
q_t = pa - partons[0].p;
} else {
// top emission from top quark leg
q_t = pa - partons[0].p - boson.p;
}
fastjet::PseudoJet q_b = pa - partons[0].p;
auto first_idx = cbegin(partons);
auto last_idx = cend(partons) - 1;
// for uno/exqqbar the two extremal partons do not contribute to the virtual corrections
if (event.type() == event_type::unob
|| event.type() == event_type::qqbar_exb) {
// unordered gluon or first quark is partons[0] and is already subtracted
// partons[1] is the backward quark (uno case) or the second quark (exqqbar case)
q_t -= partons[1].p;
q_b -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof
|| event.type() == event_type::qqbar_exf) {
// end sum at next to last parton
--last_idx;
}
// if central qqbar event, virtual correction do not occur between qqbar pair
auto first_idx_qqbar = last_idx;
auto last_idx_qqbar = last_idx;
if(event.type() == event_type::central_qqbar){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
assert(backquark!=end(partons)-1 && (backquark+1)->type!=pid::gluon);
last_idx = backquark;
first_idx_qqbar = last_idx+1;
}
auto res = virtual_corrections_interference(first_idx, last_idx, q_t, q_b, mur);
if(last_idx != first_idx_qqbar) {
q_t -= (last_idx+1)->p;
q_b -= (last_idx+1)->p;
if(is_gluon(in.front().type)) {
// "top" emission from qqbar, "bot" emission from bottom quark leg
q_t -= boson.p;
} else {
// "top" emission from top quark leg, "bot" emission from qqbar
q_b -= boson.p;
}
auto res_2 = virtual_corrections_interference(first_idx_qqbar, last_idx_qqbar, q_t, q_b, mur);
for(std::size_t i=0; i<res.size(); ++i){
res[i] += res_2[i];
}
}
if(param_.nlo.enabled) {
if(partons.size() > param_.nlo.nj) return std::vector(res.size(), 1.);
assert(partons.size() == param_.nlo.nj);
for(double & virt: res) virt += 1.;
} else {
for(double & virt: res) virt = exp(virt);
}
return res;
}
// Virtual corrections for Z or photon processes with 3 interfering
// contributions
// Returns 6 values (3 squares + 3 interference terms)
std::vector <double> MatrixElement::virtual_corrections_Zphoton_triple(
Event const & event,
const double mur,
Particle const & boson
) const{
assert(event.type() == event_type::central_qqbar);
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());
enum emission_site {top, mid, bot};
fastjet::PseudoJet q_t = pa - partons[0].p - boson.p;
fastjet::PseudoJet q_b = pa - partons[0].p;
auto first_idx = cbegin(partons);
auto last_idx = cend(partons) - 1;
auto first_idx_qqbar = last_idx;
auto last_idx_qqbar = last_idx;
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
assert(backquark!=end(partons)-1 && (backquark+1)->type!=pid::gluon);
last_idx = backquark;
first_idx_qqbar = last_idx+1;
// here q[bot] = q[mid], thus we only need to evaluate virtual_corrections_interference[top][mid]
const auto res_1_top_mid = virtual_corrections_interference(first_idx, last_idx, q_t, q_b, mur);
double res_tmp[3][3];
res_tmp[top][top] = res_1_top_mid.at(0);
res_tmp[mid][mid] = res_1_top_mid.at(1);
res_tmp[bot][bot] = res_tmp[mid][mid];
res_tmp[top][mid] = res_1_top_mid.at(2);
res_tmp[top][bot] = res_tmp[top][mid];
res_tmp[mid][bot] = res_tmp[mid][mid];
if(last_idx != first_idx_qqbar) {
q_t -= (last_idx+1)->p;
q_b -= (last_idx+1)->p;
// here q[mid] = q[top], thus we only need to evaluate virtual_corrections_interference[top][bot]
const auto res_2_top_bot = virtual_corrections_interference(first_idx_qqbar, last_idx_qqbar,
q_t, q_b, mur);
res_tmp[top][top] += res_2_top_bot.at(0);
res_tmp[mid][mid] += res_2_top_bot.at(0); // same as [top][top]
res_tmp[bot][bot] += res_2_top_bot.at(1);
res_tmp[top][mid] += res_2_top_bot.at(0); // same as [top][top]
res_tmp[top][bot] += res_2_top_bot.at(2);
res_tmp[mid][bot] += res_2_top_bot.at(2); // same as [top][bot]
}
std::vector<double> res = {res_tmp[top][top], res_tmp[mid][mid], res_tmp[bot][bot],
res_tmp[top][mid], res_tmp[top][bot], res_tmp[mid][bot]};
if(param_.nlo.enabled) {
if(partons.size() > param_.nlo.nj) return std::vector(res.size(), 1.);
assert(partons.size() == param_.nlo.nj);
for(double & virt: res) virt += 1.;
} else {
for(double & virt: res) virt = exp(virt);
}
return res;
}
std::vector<double> MatrixElement::virtual_corrections_Zphoton(
Event const & event,
const double mur,
Particle const & boson
) const{
using namespace event_type;
auto const & in = event.incoming();
if(event.type() == central_qqbar) {
if(is_gluon(in.front().type)) {
if(is_gluon(in.back().type)) {
// gg event
return {virtual_corrections_Zphoton_single(event, mur, boson, false)};
} else {
// gq event
return virtual_corrections_Zphoton_double(event, mur, boson);
}
} else {
if(is_gluon(in.back().type)) {
// qg event
return virtual_corrections_Zphoton_double(event, mur, boson);
} else {
// qq event
return virtual_corrections_Zphoton_triple(event, mur, boson);
}
}
}
if(event.type() == qqbar_exb && is_gluon(in.back().type)) {
// gg event, emission from backward leg
return {virtual_corrections_Zphoton_single(event, mur, boson, false)};
}
if(event.type() == qqbar_exf && is_gluon(in.front().type)) {
// gg event, emission from forward leg
return {virtual_corrections_Zphoton_single(event, mur, boson, true)};
}
if(event.type() == FKL || event.type() == unob || event.type() == unof) {
if(is_gluon(in.back().type)) {
// qg event, emission from backward leg
return {virtual_corrections_Zphoton_single(event, mur, boson, false)};
}
if(is_gluon(in.front().type)) {
// gq event, emission from forward leg
return {virtual_corrections_Zphoton_single(event, mur, boson, true)};
}
}
// qq initiated FKL/uno or qg/gq initiated exqqbar
return virtual_corrections_Zphoton_double(event, mur, boson);
}
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
|| AWZH_boson.type == pid::Z
|| AWZH_boson.type == pid::photon
){
return virtual_corrections_Zphoton(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;
auto first_idx = cbegin(out);
auto last_idx = cend(out) - 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;
}
auto first_idx_qqbar = last_idx;
auto 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)); }
);
assert(backquark!=end(out)-1 && (backquark+1)->type!=pid::gluon);
last_idx = backquark;
first_idx_qqbar = last_idx+1;
}
double exponent = virtual_corrections_no_interference(first_idx, last_idx, q, mur);
if(last_idx != first_idx_qqbar) {
q -= (last_idx+1)->p;
exponent += virtual_corrections_no_interference(first_idx_qqbar, last_idx_qqbar, q, mur);
}
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_q * K(bptype, pn, pb)*currents::ME_qqbar_qg(pgin, pb, p1, p2, pn) / (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 partons Vector of all outgoing partons
* @param qbar_first Ordering of the qqbar pair (true: qbar-q, false: q-qbar)
* @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,
std::vector<CLHEP::HepLorentzVector> const & partons,
bool const qbar_first
){
using namespace currents;
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,
qbar_first, nabove
) / (t1 * t2 * (4.*(N_C*N_C - 1)));
}
/** 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;
assert(!(aptype==pid::gluon && bptype==pid::gluon));
if(aptype == pid::gluon || wc) {
// swap currents to ensure that the W is emitted off the first one
return ME_W_current(bptype, aptype, p1, pa, pn, pb, plbar, pl, false, Wprop);
}
// we assume that the W is emitted off a quark line
// if this is not the case, we have to apply CP conjugation,
// which is equivalent to swapping lepton and antilepton momenta
const double current_contr = is_quark(aptype)?
ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop):
ME_W_qQ(p1,pl,plbar,pa,pn,pb,Wprop);
const double t1 = (pa - p1 - pl - plbar).m2();
const double tn = (pn - pb).m2();
return K(aptype, p1, pa)
* K(bptype, pn, pb)
* current_contr/(4.*(N_C*N_C - 1) * t1 * tn);
}
/** 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 plbar,
CLHEP::HepLorentzVector pl,
bool const wc, ParticleProperties const & Wprop
){
using namespace currents;
assert(bptype != pid::gluon || aptype != pid::gluon);
if(aptype == pid::gluon || wc) {
// emission off pb -- pn line
// we assume that the W is emitted off a quark line
// if this is not the case, we have to apply CP conjugation,
// which is equivalent to swapping lepton and antilepton momenta
if(is_antiquark(bptype)) std::swap(plbar, pl);
const double t1 = (pa - p1 - pg).m2();
const double tn = (pb - pn - plbar - pl).m2();
return K_q*K_q*ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop)/(4.*(N_C*N_C - 1) * t1 * tn);
}
// emission off pa -- p1 line
if(is_antiquark(aptype)) std::swap(plbar, pl);
const double t1 = (pa - p1 - pg - plbar - pl).m2();
const double tn = (pb - pn).m2();
return K(bptype, pn, pb)/C_F*ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop)/(t1 * tn);
}
/** \brief Matrix element squared for backward qqbar tree-level current-current
* scattering With W+Jets
*
* @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 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 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 wc,
ParticleProperties const & Wprop
){
using namespace currents;
if(is_anyquark(bptype) && wc) {
// W Must be emitted from forwards leg.
const double t1 = (pa - pq - pqbar).m2();
const double tn = (pb - pn - pl - plbar).m2();
return K_q * K_q * ME_W_Exqqbar_QQq(
pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop
) / (4.*(N_C*N_C - 1) * t1 * tn);
}
const double t1 = (pa - pl - plbar - pq - pqbar).m2();
const double tn = (pb - pn).m2();
return K(bptype, pn, pb)/C_F
* ME_WExqqbar_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop) / (t1 * tn);
}
/* \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 partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param qbar_first Ordering of the qqbar pair (true: qbar-q, false: q-qbar)
* @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,
std::vector<CLHEP::HepLorentzVector> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const qbar_first, bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
const double wt = K(aptype, partons.front(), pa) * K(bptype, partons.back(), pb) / (4.*(N_C*N_C - 1));
if(wqq)
return wt*ME_WCenqqbar_qq(pa, pb, pl, plbar, partons,
is_antiquark(aptype),is_antiquark(bptype),
qbar_first, nabove, Wprop);
return wt*ME_W_Cenqqbar_qq(pa, pb, pl, plbar, partons,
is_antiquark(aptype), is_antiquark(bptype),
qbar_first, 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,
const ParticleID ltype,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
const auto pZ = pl + plbar;
const double pref = K(aptype, p1, pa) * K(bptype, pn, pb)/(4.*(N_C*N_C-1));
// 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
const double t1 = (pa-p1-pZ).m2();
const double tn = (pb-pn ).m2();
return { pref*ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,ltype,Zprop,stw2,ctw)/(t1 * tn) };
}
if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
const double t1 = (pa-p1 ).m2();
const double tn = (pb-pn-pZ).m2();
return { pref*ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,ltype,Zprop,stw2,ctw)/(t1 * tn) };
}
// This is a qq event
assert(is_anyquark(aptype) && is_anyquark(bptype));
const double t1_top = (pa-p1-pZ).m2();
const double t2_top = (pb-pn ).m2();
const double t1_bot = (pa-p1 ).m2();
const double t2_bot = (pb-pn-pZ).m2();
std::vector<double> res = ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,ltype,Zprop,stw2,ctw);
assert(res.size() == 3);
res[0] *= pref/(t1_top * t2_top);
res[1] *= pref/(t1_bot * t2_bot);
res[2] *= pref/sqrt(t1_top * t2_top * t1_bot * t2_bot);
return res;
}
/** 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,
const ParticleID ltype,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
const auto pZ = pl + plbar;
const double pref = K(aptype, p1, pa)/C_F * K(bptype, pn, pb)/C_F;
// we only evaluate unordered backward -> first incoming particle must be a quark
assert(is_anyquark(aptype));
const double t1_top = (pa-pg-p1-pZ).m2();
const double t2_top = (pb-pn ).m2();
if (is_gluon(bptype)) {
// This is a qg event -> Z emission from top leg
return {
pref*ME_Zuno_qg(
pa,pb,pg,p1,pn,plbar,pl,
aptype,bptype,ltype,
Zprop,stw2,ctw
)/(t1_top * t2_top)
};
}
// This is a qq event
assert(is_anyquark(bptype));
const double t1_bot = (pa-pg-p1).m2();
const double t2_bot = (pb-pn-pZ).m2();
std::vector<double> res = ME_Zuno_qQ(
pa,pb,pg,p1,pn,plbar,pl,
aptype,bptype,ltype,
Zprop,stw2,ctw
);
assert(res.size() == 3);
res[0] *= pref/(t1_top * t2_top);
res[1] *= pref/(t1_bot * t2_bot);
res[2] *= pref/sqrt(t1_top * t2_top * t1_bot * t2_bot);
return res;
}
/** Matrix element squared for backwards extremal qqbar tree-level current-current
* scattering With Z+Jets
*
* @param qptype PDG ID of final state quark in qqbar pair
* @param bptype Incoming particle b PDG ID
* @param pa Incoming particle a momentum
* @param pb Incoming particle b momentum
* @param pq Final state quark momentum
* @param pqbar Final state anti-quark momentum
* @param pn Final state particle n 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 qqbar_exb Tree-Level Current-Current Scattering
*
* @note The qqbar_exf contribution can be calculated by reversing the argument ordering.
*/
std::vector<double> ME_Z_exqqbar_current(
const ParticleID qptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
const ParticleID ltype,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
const auto pZ = pl + plbar;
const double t1_top = (pa-pq-pqbar-pZ).m2();
const double t2_top = (pb-pn).m2();
if (is_gluon(bptype)) {
// This is a gg event -> Z emission from top leg
return { K(bptype, pn, pb)/C_F
* ME_ZExqqbar_gg(pa,pb,pq,pqbar,pn,plbar,pl,qptype,ltype,Zprop,stw2,ctw)
/ (t1_top * t2_top) };
}
// This is a gq event
assert(is_anyquark(bptype));
const double t1_bot = (pa-pq-pqbar).m2();
const double t2_bot = (pb-pn-pZ).m2();
std::vector<double> res = ME_ZExqqbar_gq(
pa,pb,pq,pqbar,pn,plbar,pl,
qptype,bptype,ltype,
Zprop,stw2,ctw
);
assert(res.size() == 3);
res[0] /= (t1_top * t2_top);
res[1] /= (t1_bot * t2_bot);
res[2] /= sqrt(t1_top * t2_top * t1_bot * t2_bot);
return res;
}
/** Matrix element squared for central qqbar tree-level current-current
* scattering With Z+Jets
*
* @param aptype Incoming particle a PDG ID
* @param bptype Incoming particle b PDG ID
* @param qptype PDG ID of final state quark in qqbar pair
* @param qbar_first Ordering of the qqbar pair (true: qbar-q, false: q-qbar)
* @param nabove Number of gluons emitted before central qqbarpair
* @param pa Incoming particle a momentum
* @param pb Incoming particle b momentum
* @param partons Vector of all outgoing partons
* @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 central qqbar Tree-Level Current-Current Scattering
*/
std::vector<double> ME_Z_qqbar_mid_current(
const ParticleID aptype, const ParticleID bptype,
const ParticleID qptype,
const bool qbar_first, const int nabove,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
std::vector<CLHEP::HepLorentzVector> const & partons,
const ParticleID ltype,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
std::vector<double> res;
if(is_gluon(aptype)) {
if(is_gluon(bptype)) {
// gg event
res = {
ME_ZCenqqbar_gg(
pa,pb,plbar,pl,
partons,qbar_first,nabove,
qptype,ltype,
Zprop,stw2,ctw
)
};
} else {
// gq event
res = ME_ZCenqqbar_gq(
pa,pb,plbar,pl,
partons,qbar_first,nabove,
bptype,qptype,ltype,
Zprop,stw2,ctw
);
}
} else {
assert(is_anyquark(bptype));
// qq event
res = ME_ZCenqqbar_qq(
pa,pb,plbar,pl,
partons,qbar_first,nabove,
aptype,bptype,qptype,ltype,
Zprop,stw2,ctw
);
}
const double wt = K(aptype, partons.front(), pa) * K(bptype, partons.back(), pb) / (4.*(N_C*N_C - 1));
for(double & me: res) me *= wt;
return res;
}
/** Matrix element squared for tree-level current-current scattering With photon+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 p_photon Final state photon momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
std::vector<double> ME_photon_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 & p_photon
){
using namespace currents;
const double pref = K(aptype, p1, pa) * K(bptype, pn, pb)/(4.*(N_C*N_C-1));
// 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
const double t1 = (pa-p1-p_photon).m2();
const double tn = (pb-pn ).m2();
return { pref*ME_photon_qg(pa,pb,p1,pn,p_photon,aptype,bptype)/(t1 * tn) };
}
if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
const double t1 = (pa-p1 ).m2();
const double tn = (pb-pn-p_photon).m2();
return { pref*ME_photon_qg(pb,pa,pn,p1,p_photon,bptype,aptype)/(t1 * tn) };
}
// This is a qq event
assert(is_anyquark(aptype) && is_anyquark(bptype));
const double t1_top = (pa-p1-p_photon).m2();
const double t2_top = (pb-pn ).m2();
const double t1_bot = (pa-p1 ).m2();
const double t2_bot = (pb-pn-p_photon).m2();
std::vector<double> res = ME_photon_qQ(pa,pb,p1,pn,p_photon,aptype,bptype);
assert(res.size() == 3);
res[0] *= pref/(t1_top * t2_top);
res[1] *= pref/(t1_bot * t2_bot);
res[2] *= pref/sqrt(t1_top * t2_top * t1_bot * t2_bot);
return res;
}
/** Matrix element squared for tree-level current-current scattering With photon+Jets (UNO)
* @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 p_photon Final state photon momentum
*
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
std::vector<double> ME_photon_uno_current(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & p_photon,
CLHEP::HepLorentzVector const & pg
){
using namespace currents;
const double pref = K(aptype, p1, pa)/C_F * K(bptype, pn, pb)/C_F;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
// we only evaluate unordered backward -> first incoming particle must be a quark
assert(is_anyquark(aptype));
if(is_anyquark(aptype) && is_gluon(bptype)){
// This is a qg event
- const double t1 = (pa-p1-p_photon).m2();
+ const double t1 = (pa-p1-p_photon-pg).m2();
const double tn = (pb-pn ).m2();
- return { pref/**ME_photon_qg(pa,pb,p1,pn,p_photon,aptype,bptype)/(t1 * tn)*/ };
+ return { pref*ME_photon_uno_qg(pa,pb,pg,p1,pn,p_photon,aptype,bptype)/(t1 * tn) };
}
if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
const double t1 = (pa-p1 ).m2();
- const double tn = (pb-pn-p_photon).m2();
- return { pref/**ME_photon_qg(pb,pa,pn,p1,p_photon,bptype,aptype)/(t1 * tn)*/ };
+ const double tn = (pb-pn-p_photon-pg).m2();
+ return { pref*ME_photon_uno_qg(pb,pa,pg,pn,p1,p_photon,bptype,aptype)/(t1 * tn)};
}
// This is a qq event
assert(is_anyquark(aptype) && is_anyquark(bptype));
const double t1_top = (pa-p1-pg-p_photon).m2();
const double t2_top = (pb-pn).m2();
const double t1_bot = (pa-p1-pg).m2();
const double t2_bot = (pb-pn-p_photon).m2();
std::vector<double> res = ME_photon_uno_qQ(pa,pb,p1,pn,p_photon,pg,aptype,bptype);
assert(res.size() == 3);
res[0] *= pref/(t1_top * t2_top);
res[1] *= pref/(t1_bot * t2_bot);
res[2] *= pref/sqrt(t1_top * t2_top * t1_bot * t2_bot);
return res;
}
/** \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:
case pid::Z_photon_mix:
return tree_kin_Z(ev);
case pid::photon:
return tree_kin_photon(ev);
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
){
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);
const bool qbar_first = is_antiquark(backmidquark->type);
const auto pq = to_HepLorentzVector(*(backmidquark+(qbar_first?1:0)));
const auto pqbar = to_HepLorentzVector(*(backmidquark+(qbar_first?0:1)));
const auto p1 = to_HepLorentzVector(partons[0]);
const auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder_1 = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder_2 = cend(partons) - 1;
const int nabove = std::distance(begin_ladder_1, end_ladder_1);
const auto q0 = pa - p1;
// t-channel momentum after qqbar
auto q3 = q0;
for(auto parton_it = begin_ladder_1; parton_it != begin_ladder_2; ++parton_it){
q3 -= to_HepLorentzVector(*parton_it);
}
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for(Particle const & parton: partons) {
partonsHLV.push_back(to_HepLorentzVector(parton));
}
const double current_factor = ME_qqbar_mid_current(
aptype, bptype, nabove, pa, pb,
partonsHLV, qbar_first
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder_1, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder_2,
q3, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_qqbar(InIter begin_in, InIter end_in, partIter begin_part,
partIter end_part, double lambda){
const auto pgin = to_HepLorentzVector(*begin_in);
const auto pb = to_HepLorentzVector(*(end_in-1));
const auto p1 = to_HepLorentzVector(*begin_part);
const auto p2 = to_HepLorentzVector(*(begin_part+1));
const auto pn = to_HepLorentzVector(*(end_part-1));
assert((begin_in)->type==pid::gluon); // Incoming a must be gluon.
const double current_factor = ME_qqbar_current(
(end_in-1)->type, pgin, p1, p2, pn, pb
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(begin_part+2), (end_part-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 begin_in, InIter end_in, partIter begin_part,
partIter end_part, double lambda
){
const auto pa = to_HepLorentzVector(*begin_in);
const auto pb = to_HepLorentzVector(*(end_in-1));
const auto pg = to_HepLorentzVector(*begin_part);
const auto p1 = to_HepLorentzVector(*(begin_part+1));
const auto pn = to_HepLorentzVector(*(end_part-1));
const double current_factor = ME_uno_current(
(begin_in)->type, (end_in-1)->type, pg, pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(begin_part+2), (end_part-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 begin_in, partIter begin_part,
partIter end_part,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const auto pa = to_HepLorentzVector(*begin_in);
const auto pb = to_HepLorentzVector(*(begin_in+1));
const auto pg = to_HepLorentzVector(*begin_part);
const auto p1 = to_HepLorentzVector(*(begin_part+1));
const auto pn = to_HepLorentzVector(*(end_part-1));
bool wc = (begin_in)->type==(begin_part+1)->type; //leg b emits w
auto q0 = pa - p1 - pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_uno_current(
(begin_in)->type, (begin_in+1)->type, pn, pb,
p1, pa, pg, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_part+2, end_part-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_qqbar(InIter begin_in, partIter begin_part,
partIter end_part,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
double lambda, ParticleProperties const & Wprop
){
const bool qbar_first=is_quark(*begin_part);
const auto pa = to_HepLorentzVector(*begin_in);
const auto pb = to_HepLorentzVector(*(begin_in+1));
const auto pq = to_HepLorentzVector(*(begin_part+(qbar_first?0:1)));
const auto pqbar = to_HepLorentzVector(*(begin_part+(qbar_first?1:0)));
const auto p1 = to_HepLorentzVector(*(begin_part));
const auto pn = to_HepLorentzVector(*(end_part-1));
const bool wc = (begin_in+1)->type!=(end_part-1)->type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqbar_current(
(begin_in+1)->type, pa, pb,
pq, pqbar, pn, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_part+2, end_part-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
){
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);
const bool qbar_first = is_antiquark(backmidquark->type);
const auto pq = to_HepLorentzVector(*(backmidquark+(qbar_first?1:0)));
const auto pqbar = to_HepLorentzVector(*(backmidquark+(qbar_first?0:1)));
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto begin_ladder_1 = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder_2 = cend(partons) - 1;
const int nabove = std::distance(begin_ladder_1, end_ladder_1);
const int nbelow = std::distance(begin_ladder_2, end_ladder_2);
const bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqbar emit W
const bool wc = !wqq && (aptype==partons.front().type); //leg b emits w
assert(!wqq || !wc);
auto q0 = pa - p1;
// t-channel momentum after qqbar
auto q3 = q0;
if(wqq){ // emission from qqbar
q3 -= pl + plbar;
} else if(!wc) { // emission from leg a
q0 -= pl + plbar;
q3 -= pl + plbar;
}
for(auto parton_it = begin_ladder_1; parton_it != begin_ladder_2; ++parton_it){
q3 -= to_HepLorentzVector(*parton_it);
}
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for(Particle const & parton: partons) {
partonsHLV.push_back(to_HepLorentzVector(parton));
}
const double current_factor = ME_W_qqbar_mid_current(
aptype, bptype, nabove, nbelow, pa, pb, partonsHLV,
plbar, pl, qbar_first, wqq, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder_1, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder_2,
q3, 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
){
assert(is_anyquark(aptype));
assert(is_anyquark(bptype));
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
);
assert(current_factor.size() == 3);
assert(current_factor.size() == ladder_factor.size());
std::vector<double> result(current_factor.size());
for(size_t i=0; i<current_factor.size(); ++i){
result[i] = K_q*K_q/(4.*(N_C*N_C - 1.))
*current_factor.at(i)
*ladder_factor.at(i);
}
const double t1_top = q0.m2();
const double t2_top = (pb-pn-pl2bar-pl2).m2();
const double t1_bot = q1.m2();
const double t2_bot = (pb-pn-pl1bar-pl1).m2();
result[0] /= t1_top * t2_top;
result[1] /= t1_bot * t2_bot;
result[2] /= sqrt(t1_top * t2_top * t1_bot * t2_bot);
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,
const ParticleID ltype,
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,
ltype, 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 begin_in, partIter begin_part, partIter end_part,
const ParticleID ltype,
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(*begin_in);
const auto pb = to_HepLorentzVector(*(begin_in+1));
const auto pg = to_HepLorentzVector(*begin_part);
const auto p1 = to_HepLorentzVector(*(begin_part+1));
const auto pn = to_HepLorentzVector(*(end_part-1));
const ParticleID aptype = (begin_in)->type;
const ParticleID bptype = (begin_in+1)->type;
// we only evaluate unordered backward -> first incoming particle must be a quark
assert(is_anyquark(aptype));
const std::vector <double> current_factor = ME_Z_uno_current(
aptype, bptype, pn, pb, p1, pa, pg,
ltype, plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
const auto q0 = pa-pg-p1-plbar-pl;
if(is_gluon(bptype)){
// This is a qg event
ladder_factor.push_back(FKL_ladder_weight(begin_part+2, end_part-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q1 = pa-pg-p1;
ladder_factor=FKL_ladder_weight_mix(begin_part+2, end_part-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;
}
template<class InIter, class partIter>
std::vector <double> tree_kin_Z_qqbar(
InIter begin_in, partIter begin_part, partIter end_part,
const ParticleID ltype,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw
){
const bool qbar_first = is_antiquark(*begin_part);
const auto pa = to_HepLorentzVector(*begin_in);
const auto pb = to_HepLorentzVector(*(begin_in+1));
const auto pq = to_HepLorentzVector(*(begin_part+(qbar_first?1:0)));
const auto pqbar = to_HepLorentzVector(*(begin_part+(qbar_first?0:1)));
const auto p1 = to_HepLorentzVector(*(begin_part));
const auto pn = to_HepLorentzVector(*(end_part-1));
const ParticleID qptype = (begin_part+(qbar_first?1:0))->type;
const ParticleID bptype = (begin_in+1)->type;
const std::vector <double> current_factor = ME_Z_exqqbar_current(
qptype, bptype, pa, pb, pq, pqbar, pn,
ltype, plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
const auto q0 = pa-pq-pqbar-plbar-pl;
if(is_gluon(bptype)){
// This is a gg event
ladder_factor.push_back(FKL_ladder_weight(begin_part+2, end_part-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a gq event
const auto q1 = pa-pq-pqbar;
ladder_factor=FKL_ladder_weight_mix(begin_part+2, end_part-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;
}
template<class InIter, class partIter>
std::vector<double> tree_kin_Z_qqbar_mid(
InIter begin_in, partIter begin_part, partIter end_part,
const ParticleID ltype,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw
){
const ParticleID aptype = (begin_in)->type;
const ParticleID bptype = (begin_in+1)->type;
const auto pa = to_HepLorentzVector(*begin_in);
const auto pb = to_HepLorentzVector(*(begin_in+1));
const auto backmidquark = std::find_if(
begin_part+1, end_part-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end_part-1);
const bool qbar_first = is_antiquark(backmidquark->type);
const auto pq = to_HepLorentzVector(*(backmidquark+(qbar_first?1:0)));
const auto pqbar = to_HepLorentzVector(*(backmidquark+(qbar_first?0:1)));
const ParticleID qptype = (backmidquark+(qbar_first?1:0))->type;
const auto p1 = to_HepLorentzVector(*(begin_part));
const auto pn = to_HepLorentzVector(*(end_part-1));
const auto begin_ladder_1 = begin_part + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder_2 = end_part - 1;
const int nabove = std::distance(begin_ladder_1, end_ladder_1);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(std::distance(begin_part, end_part));
for(auto parton_it = begin_part; parton_it != end_part; ++parton_it){
partonsHLV.push_back(to_HepLorentzVector(*parton_it));
}
const std::vector<double> current_factor = ME_Z_qqbar_mid_current(
aptype, bptype, qptype, qbar_first, nabove, pa, pb,
partonsHLV, ltype, plbar, pl, Zprop, stw2, ctw
);
std::vector<double> ladder_factor;
if(is_gluon(aptype)) {
if(is_gluon(bptype)) {
// gg event -> Z emitted from central qqbar
// first t-channel momentum
const auto q0 = pa - p1;
// t-channel momentum after qqbar
auto q3 = q0 - pl - plbar;
for(auto parton_it = begin_ladder_1; parton_it != begin_ladder_2; ++parton_it){
q3 -= to_HepLorentzVector(*parton_it);
}
ladder_factor.push_back(
FKL_ladder_weight(begin_ladder_1, end_ladder_1, q0, pa, pb, p1, pn, lambda)
* FKL_ladder_weight(begin_ladder_2, end_ladder_2, q3, pa, pb, p1, pn, lambda)
);
} else {
// gq event -> Z emitted from central qqbar or bottom leg
auto q_bot = pa - p1;
// q0_mid = q0_bot -> no need to evaluate FKL_ladder_weight_mix for the first ladder
const double ladder_1 = FKL_ladder_weight(begin_ladder_1, end_ladder_1,
q_bot, pa, pb, p1, pn, lambda);
for(auto parton_it = begin_ladder_1; parton_it != begin_ladder_2; ++parton_it) {
q_bot -= to_HepLorentzVector(*parton_it);
}
auto q_mid = q_bot - pl - plbar;
const auto ladder_2 = FKL_ladder_weight_mix(begin_ladder_2, end_ladder_2,
q_mid, q_bot, pa, pb, p1, pn, lambda);
for(size_t i=0; i<ladder_2.size(); ++i){
ladder_factor.push_back(ladder_1 * ladder_2.at(i));
}
}
} else {
assert(is_anyquark(bptype));
// qq event -> 3 contributions
auto q_t = pa - p1 - pl - plbar;
auto q_b = pa - p1;
// q0_bot = q0_mid -> only need to evaluate FKL_ladder_weight_mix[top][mid]
const auto ladder_1_top_mid = FKL_ladder_weight_mix(begin_ladder_1, end_ladder_1,
q_t, q_b, pa, pb, p1, pn, lambda);
enum emission_site {top, mid, bot};
double ladder[3][3];
ladder[top][top] = ladder_1_top_mid.at(0);
ladder[mid][mid] = ladder_1_top_mid.at(1);
ladder[bot][bot] = ladder[mid][mid];
ladder[top][mid] = ladder_1_top_mid.at(2);
ladder[top][bot] = ladder[top][mid];
ladder[mid][bot] = ladder[mid][mid];
for(auto parton_it = begin_ladder_1; parton_it != begin_ladder_2; ++parton_it) {
q_b -= to_HepLorentzVector(*parton_it);
}
q_t = q_b - pl - plbar;
// q3_mid = q3_top -> only need to evaluate FKL_ladder_weight_mix[top][bot]
const auto ladder_2_top_bot = FKL_ladder_weight_mix(begin_ladder_2, end_ladder_2,
q_t, q_b, pa, pb, p1, pn, lambda);
ladder[top][top] *= ladder_2_top_bot.at(0);
ladder[mid][mid] *= ladder_2_top_bot.at(0); // same as [top][top]
ladder[bot][bot] *= ladder_2_top_bot.at(1);
ladder[top][mid] *= ladder_2_top_bot.at(0); // same as [top][top]
ladder[top][bot] *= ladder_2_top_bot.at(2);
ladder[mid][bot] *= ladder_2_top_bot.at(2); // same as [top][bot]
ladder_factor = {ladder[top][top], ladder[mid][mid], ladder[bot][bot],
ladder[top][mid], ladder[top][bot], ladder[mid][bot]};
}
assert(current_factor.size() == ladder_factor.size());
std::vector<double> res;
for(size_t i=0; i<current_factor.size(); ++i) {
res.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return res;
}
} // 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)) && decay.at(0).type==-decay.at(1).type);
// get leptons
CLHEP::HepLorentzVector plbar;
CLHEP::HepLorentzVector pl;
ParticleID ltype;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
ltype = decay.at(1).type;
}
else{
pl = to_HepLorentzVector(decay.at(0));
ltype = decay.at(0).type;
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, ltype, 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), ltype, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_forward){
auto kin_rev = tree_kin_Z_uno(crbegin(incoming), crbegin(partons),
crend(partons), ltype, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
if(!is_gluon(incoming[0].type)){
// qq unordered forward: reorder contributions such that first/second
// value corresponds to Z emission from top/bottom leg respectively
std::swap(kin_rev[0], kin_rev[1]);
}
return kin_rev;
}
if(ev.type() == extremal_qqbar_backward){
return tree_kin_Z_qqbar(cbegin(incoming), cbegin(partons),
cend(partons), ltype, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == extremal_qqbar_forward){
auto kin_rev = tree_kin_Z_qqbar(crbegin(incoming), crbegin(partons),
crend(partons), ltype, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
if(!is_gluon(incoming[0].type)){
// qqbar forward with incoming qg: reorder contributions such that first/second
// value corresponds to Z emission from top/bottom leg respectively
std::swap(kin_rev[0], kin_rev[1]);
}
return kin_rev;
}
assert(ev.type() == central_qqbar);
if(!is_gluon(incoming[0].type) && is_gluon(incoming[1].type)){
// qg initiated central qqbar: reverse event to evaluate as gq
auto kin_rev = tree_kin_Z_qqbar_mid(crbegin(incoming), crbegin(partons),
crend(partons), ltype, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
std::swap(kin_rev[0], kin_rev[1]);
return kin_rev;
}
return tree_kin_Z_qqbar_mid(cbegin(incoming), cbegin(partons),
cend(partons), ltype, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
namespace {
std::vector <double> tree_kin_photon_FKL(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & p_photon,
const double lambda
){
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_photon_current(
aptype, bptype, pn, pb, p1, pa, p_photon
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-p1-p_photon;
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-p_photon;
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_photon_uno(
InIter begin_in, partIter begin_part, partIter end_part,
const CLHEP::HepLorentzVector & pp,
const double lambda
){
const auto pa = to_HepLorentzVector(*begin_in);
const auto pb = to_HepLorentzVector(*(begin_in+1));
const auto pg = to_HepLorentzVector(*begin_part);
const auto p1 = to_HepLorentzVector(*(begin_part+1));
const auto pn = to_HepLorentzVector(*(end_part-1));
const ParticleID aptype = (begin_in)->type;
const ParticleID bptype = (begin_in+1)->type;
// we only evaluate unordered backward -> second incoming particle must be a quark
assert(is_anyquark(aptype));
const std::vector <double> current_factor = ME_photon_uno_current(
aptype, bptype, pa, pb, p1, pn, pp, pg
);
std::vector <double> ladder_factor;
const auto q0 = pa-pg-p1-pp;
if(is_gluon(bptype)){
// This is a qg event
ladder_factor.push_back(FKL_ladder_weight(begin_part+2, end_part-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q1 = pa-pg-p1;
ladder_factor=FKL_ladder_weight_mix(begin_part+2, end_part-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;
}
}
std::vector<double> MatrixElement::tree_kin_photon(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const & outgoing(ev.outgoing());
assert(ev.decays().empty());
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
const auto the_photon = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::photon; }
);
assert(the_photon != end(outgoing));
const auto p_photon = to_HepLorentzVector(*the_photon);
switch(ev.type()) {
- case FKL:
+ case FKL:{
return tree_kin_photon_FKL(
incoming[0].type, incoming[1].type,
pa, pb, partons, p_photon,
param_.regulator_lambda
);
- case unob:
- return tree_kin_photon_uno(
+ }
+ case unob:{
+ return tree_kin_photon_uno(
cbegin(incoming), cbegin(partons), cend(partons),
p_photon, param_.regulator_lambda
);
- case unof:
- return tree_kin_photon_uno(
+ }
+ case unof:{
+ std::vector <double> kin_rev = tree_kin_photon_uno(
crbegin(incoming), crbegin(partons), crend(partons),
p_photon, param_.regulator_lambda
);
+ if(!is_gluon(incoming[0].type)){
+ // qq unordered forward: reorder contributions such that first/second
+ // value corresponds to Z emission from top/bottom leg respectively
+ std::swap(kin_rev[0], kin_rev[1]);
+ }
+ return kin_rev;
+ }
default:
// TODO: rest
throw HEJ::not_implemented{"photon + jets matrix elements beyond LL"};
}
}
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 begin_in, InIter end_in, partIter begin_part,
partIter end_part,
CLHEP::HepLorentzVector const & qH,
CLHEP::HepLorentzVector const & qHp1,
double mt, bool inc_bot, double mb, double vev
){
const auto pa = to_HepLorentzVector(*begin_in);
const auto pb = to_HepLorentzVector(*(end_in-1));
const auto pg = to_HepLorentzVector(*begin_part);
const auto p1 = to_HepLorentzVector(*(begin_part+1));
const auto pn = to_HepLorentzVector(*(end_part-1));
return ME_Higgs_current_uno(
(begin_in)->type, (end_in-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 = 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 = 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,
EWConstants const & ew
) {
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:
case pid::Z:
case pid::Z_photon_mix: {
const double alpha_w = ew.alpha_w();
return alpha_w*alpha_w;
}
case pid::photon:
// only one power of \alpha
// in contrast to W/Z we don't have any coupling to decay products
return 4.*M_PI*ew.alpha_em();
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
if(bosons.size() == 2) {
const double alpha_w = ew.alpha_w();
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)*res;
}
} // namespace HEJ
diff --git a/src/Photonjets.cc b/src/Photonjets.cc
index 0dabbb9..15dfd10 100644
--- a/src/Photonjets.cc
+++ b/src/Photonjets.cc
@@ -1,308 +1,339 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2024
* \copyright GPLv2 or later
*/
#include <vector>
#include "HEJ/Constants.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/Photonjets.hh"
// generated headers
#include "HEJ/currents/jphoton_j.hh"
#include "HEJ/currents/jphoton_juno.hh"
#include "HEJ/currents/jphotonuno_j.hh"
namespace HEJ::currents {
namespace {
using COM = std::complex<double>;
// X and Y as used in contractions with unordered currents
struct XY {
COM X;
COM Y;
};
//! Photon+Jets FKL Contribution
/**
* @brief Z+Jets FKL Contribution
* @param pa Incoming Particle 1 (photon emission)
* @param pb Incoming Particle 2
* @param p1 Outgoing Particle 1 (photon emission)
* @param p2 Outgoing Particle 2
* @param pphoton Outgoing photon
* @returns j_\gamma^\mu j_\mu for all helicities h1, h\gamma, h2
*/
MultiArray<COM, 2, 2, 2> jgamma_j(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pphoton
){
using helicity::plus;
using helicity::minus;
MultiArray<COM, 2, 2, 2> res;
// NOLINTNEXTLINE
#define ASSIGN_HEL(RES, J, H1, HL, H2) \
RES[H1][HL][H2] = J<H1, HL, H2>(pa, p1, pb, p2, pphoton)
ASSIGN_HEL(res, jphoton_j, plus, minus, minus);
ASSIGN_HEL(res, jphoton_j, plus, minus, plus);
ASSIGN_HEL(res, jphoton_j, plus, plus, minus);
ASSIGN_HEL(res, jphoton_j, plus, plus, plus);
#undef ASSIGN_HEL
for(auto hphoton: {minus, plus}) {
for(auto h2: {minus, plus}) {
res[minus][hphoton][h2] = std::conj(res[plus][flip(hphoton)][flip(h2)]);
}
}
return res;
}
} // Anonymous Namespace
std::vector <double> ME_photon_qQ(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pphoton,
const ParticleID aptype, const ParticleID bptype
) {
using helicity::minus;
using helicity::plus;
MultiArray<COM, 2, 2, 2> hel_amp_top = jgamma_j(pa, pb, p1, p2, pphoton);
MultiArray<COM, 2, 2, 2> hel_amp_bot = jgamma_j(pb, pa, p2, p1, pphoton);
const auto charge_top = boost::rational_cast<double>(charge(aptype));
const auto charge_bot = boost::rational_cast<double>(charge(bptype));
double sum_top=0.;
double sum_bot=0.;
double sum_mix=0.;
for(auto h1: {minus, plus}){
for(auto hl: {minus, plus}){
for(auto h2: {minus, plus}){
const COM res_top = charge_top * hel_amp_top[h1][hl][h2];
const COM res_bot = charge_bot * hel_amp_bot[h2][hl][h1];
sum_top += norm(res_top);
sum_bot += norm(res_bot);
sum_mix += 2.0 * real(res_top * conj(res_bot));
}
}
}
return {sum_top, sum_bot, sum_mix};
}
double ME_photon_qg(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pphoton,
const ParticleID aptype, const ParticleID /*bptype*/
){
using helicity::minus;
using helicity::plus;
MultiArray<COM, 2, 2, 2> hel_amp = jgamma_j(pa, pb, p1, p2, pphoton);
const auto q_charge = boost::rational_cast<double>(charge(aptype));
double sum = 0.;
for(auto h1: {minus, plus}){
for(auto hphoton: {minus, plus}){
for(auto h2: {minus, plus}){
sum += norm(q_charge * hel_amp[h1][hphoton][h2]);
}
}
}
return sum;
}
namespace {
//! Photon+Jets UNO Contributions
/**
* @brief Photon+Jets Unordered Contribution, unordered on Photon side
* @tparam h1 Helicity of line 1 (Photon emission line)
* @tparam hp Photon Helicity
* @tparam h2 Helicity of line 2
* @tparam hg Helicity of unordered gluon
* @param pa Incoming Particle 1 (Photon and Uno emission)
* @param pb Incoming Particle 2
* @param pg Unordered Gluon
* @param p1 Outgoing Particle 1 (Photon and Uno emission)
* @param p2 Outgoing Particle 2
* @param pp Outgoing Photon
* @returns X: (U1-L), Y: (U2+l)
*
* Calculates j_Gamma_{uno}^\mu j_\mu. Ie, unordered with Photon emission same side.
*/
template<Helicity h1, Helicity hp, Helicity h2, Helicity hg>
XY amp_jgammauno_j(
const HLV & pa, const HLV & pb, const HLV & pg,
const HLV & p1, const HLV & p2, const HLV & pp
){
const COM u1 = U1_jphotonuno<h1, hp, h2, hg>(p1, p2, pa, pb, pg, pp);
const COM u2 = U2_jphotonuno<h1, hp, h2, hg>(p1, p2, pa, pb, pg, pp);
const COM l = L_jphotonuno <h1, hp, h2, hg>(p1, p2, pa, pb, pg, pp);
return {u1 - l, u2 + l};
}
MultiArray<XY, 2, 2, 2, 2> jgammauno_j(
const HLV & pa, const HLV & pb, const HLV & pg,
const HLV & p1, const HLV & p2, const HLV & pp
){
using helicity::plus;
using helicity::minus;
MultiArray<XY, 2, 2, 2, 2> xy;
// NOLINTNEXTLINE
#define ASSIGN_HEL(XY, J, H1, H2, H3, H4) \
XY[H1][H2][H3][H4] = J<H1, H2, H3, H4>(pa, pb, pg, p1, p2, pp) // NOLINT
ASSIGN_HEL(xy, amp_jgammauno_j, minus, minus, minus, minus);
ASSIGN_HEL(xy, amp_jgammauno_j, minus, minus, minus, plus);
ASSIGN_HEL(xy, amp_jgammauno_j, minus, minus, plus, minus);
ASSIGN_HEL(xy, amp_jgammauno_j, minus, minus, plus, plus);
ASSIGN_HEL(xy, amp_jgammauno_j, minus, plus, minus, minus);
ASSIGN_HEL(xy, amp_jgammauno_j, minus, plus, minus, plus);
ASSIGN_HEL(xy, amp_jgammauno_j, minus, plus, plus, minus);
ASSIGN_HEL(xy, amp_jgammauno_j, minus, plus, plus, plus);
ASSIGN_HEL(xy, amp_jgammauno_j, plus, minus, minus, minus);
ASSIGN_HEL(xy, amp_jgammauno_j, plus, minus, minus, plus);
ASSIGN_HEL(xy, amp_jgammauno_j, plus, minus, plus, minus);
ASSIGN_HEL(xy, amp_jgammauno_j, plus, minus, plus, plus);
ASSIGN_HEL(xy, amp_jgammauno_j, plus, plus, minus, minus);
ASSIGN_HEL(xy, amp_jgammauno_j, plus, plus, minus, plus);
ASSIGN_HEL(xy, amp_jgammauno_j, plus, plus, plus, minus);
ASSIGN_HEL(xy, amp_jgammauno_j, plus, plus, plus, plus);
#undef ASSIGN_HEL
return xy;
}
/**
* @brief Photon+Jets Unordered Contribution, unordered opposite to Photon side
* @tparam h1 Helicity of line 1 (photon emission)
* @tparam hp Helicity of the outgoing photon
* @tparam h2 Helicity of line 2 (unordered emission)
* @tparam hg Helicity of unordered gluon
* @param pa Incoming Particle 1 (photon emission)
* @param pb Incoming Particle 2 (unorderd emission)
* @param p1 Outgoing Particle 1 (photon emission)
* @param p2 Outgoing Particle 2 (unordered emission)
* @param pg Unordered Gluon
* @param pp Outgoing Photon
*
* @returns X: (U1-L), Y: (U2+l)
*
* Calculates j_Gamma^\mu j_{uno}_\mu. Ie, unordered with Photon emission opposite side.
*/
template<Helicity h1, Helicity hp, Helicity h2, Helicity hg>
XY amp_jgamma_juno(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pp, const HLV & pg
){
const COM u1 = U1_jphoton<h1, hp, h2, hg>(pa, pb, pp, p1, p2, pg);
const COM u2 = U2_jphoton<h1, hp, h2, hg>(pa, pb, pp, p1, p2, pg);
const COM l = L_jphoton <h1, hp, h2, hg>(pa, pb, pp, p1, p2, pg);
return {u1 - l, u2 + l};
}
MultiArray<XY, 2, 2, 2, 2> jgamma_juno(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pp, const HLV & pg
){
using helicity::plus;
using helicity::minus;
MultiArray<XY, 2, 2, 2, 2> xy;
// NOLINTNEXTLINE
#define ASSIGN_HEL(XY, J, H1, H2, H3, H4) \
XY[H1][H2][H3][H4] = J<H1, H2, H3, H4>(pa, pb, p1, p2, pp, pg)
ASSIGN_HEL(xy, amp_jgamma_juno, minus, minus, minus, minus);
ASSIGN_HEL(xy, amp_jgamma_juno, minus, minus, minus, plus);
ASSIGN_HEL(xy, amp_jgamma_juno, minus, minus, plus, minus);
ASSIGN_HEL(xy, amp_jgamma_juno, minus, minus, plus, plus);
ASSIGN_HEL(xy, amp_jgamma_juno, minus, plus, minus, minus);
ASSIGN_HEL(xy, amp_jgamma_juno, minus, plus, minus, plus);
ASSIGN_HEL(xy, amp_jgamma_juno, minus, plus, plus, minus);
ASSIGN_HEL(xy, amp_jgamma_juno, minus, plus, plus, plus);
ASSIGN_HEL(xy, amp_jgamma_juno, plus, minus, minus, minus);
ASSIGN_HEL(xy, amp_jgamma_juno, plus, minus, minus, plus);
ASSIGN_HEL(xy, amp_jgamma_juno, plus, minus, plus, minus);
ASSIGN_HEL(xy, amp_jgamma_juno, plus, minus, plus, plus);
ASSIGN_HEL(xy, amp_jgamma_juno, plus, plus, minus, minus);
ASSIGN_HEL(xy, amp_jgamma_juno, plus, plus, minus, plus);
ASSIGN_HEL(xy, amp_jgamma_juno, plus, plus, plus, minus);
ASSIGN_HEL(xy, amp_jgamma_juno, plus, plus, plus, plus);
#undef ASSIGN_HEL
return xy;
}
} // Anonymous Namespace
std::vector <double> ME_photon_uno_qQ(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2,
const HLV & pp, const HLV & pg,
const ParticleID aptype, const ParticleID bptype
){
using namespace std;
using helicity::minus;
using helicity::plus;
MultiArray<XY, 2, 2, 2, 2> hel_amp_top = jgammauno_j(pa, pb, pg, p1, p2, pp);
MultiArray<XY, 2, 2, 2, 2> hel_amp_bot = jgamma_juno(pb, pa, p2, p1, pp, pg);
const auto charge_top = boost::rational_cast<double>(charge(aptype));
const auto charge_bot = boost::rational_cast<double>(charge(bptype));
double sum_top = 0.;
double sum_bot = 0.;
double sum_mix = 0.;
for(auto h1: {minus, plus}){
for(auto hp: {minus, plus}){
for(auto h2: {minus, plus}){
for(auto hg: {minus, plus}){
const COM X_top = hel_amp_top[h1][hp][h2][hg].X;
const COM Y_top = hel_amp_top[h1][hp][h2][hg].Y;
const COM X_bot = hel_amp_bot[h2][hp][h1][hg].X;
const COM Y_bot = hel_amp_bot[h2][hp][h1][hg].Y;
// see eq:Z_uno_top, eq:Z_uno_bot & eq:Z_uno_int in developer manual
sum_top += norm(charge_top) * (C_A*C_F*C_F/2.*(norm(X_top)+norm(Y_top))
- C_F/2.*(X_top*conj(Y_top)).real());
sum_bot += norm(charge_bot) *(C_A*C_F*C_F/2.*(norm(X_bot)+norm(Y_bot))
- C_F/2.*(X_bot*conj(Y_bot)).real());
const COM xx = C_A*C_F*C_F/2. * charge_top * X_top * conj(charge_bot * X_bot);
const COM yy = C_A*C_F*C_F/2. * charge_top * Y_top * conj(charge_bot * Y_bot);
const COM xy = -C_F/4. * (charge_top * X_top * conj(charge_bot * Y_bot)
+ charge_top * Y_top * conj(charge_bot * X_bot));
sum_mix += 2.0 * real(xx + yy + xy);
}
}
}
}
//Helicity sum and average over initial states
const double pref = 1./(4.*C_A*C_A);
return {sum_top*pref, sum_bot*pref, sum_mix*pref};
}
+
+ double ME_photon_uno_qg(
+ const HLV & pa, const HLV & pb, const HLV & pg,
+ const HLV & p1, const HLV & p2, const HLV & pp,
+ const ParticleID aptype, const ParticleID /*bptype*/
+ ){
+ using namespace std;
+ using helicity::minus;
+ using helicity::plus;
+
+ MultiArray<XY, 2, 2, 2, 2> hel_amp = jgammauno_j(pa, pb, pg, p1, p2, pp);
+ const auto q_charge = boost::rational_cast<double>(charge(aptype));
+
+ double sum = 0.;
+ for(auto h1: {minus, plus}){
+ for(auto hp: {minus, plus}){
+ for(auto h2: {minus, plus}){
+ for(auto hg: {minus, plus}){
+ const COM X = hel_amp[h1][hp][h2][hg].X;
+ const COM Y = hel_amp[h1][hp][h2][hg].Y;
+ sum += norm(q_charge) * (C_A*C_F*C_F/2.*(norm(X)+norm(Y))
+ - C_F/2.*(X*conj(Y)).real());
+ }
+ }
+ }
+ }
+
+ //Helicity sum and average over initial states
+ return sum / (4.*C_A*C_A);
+ }
+
} // namespace HEJ::currents