Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7877650
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
115 KB
Subscribers
None
View Options
diff --git a/include/HEJ/Zjets.hh b/include/HEJ/Zjets.hh
index 5c15606..f3e6cf9 100644
--- a/include/HEJ/Zjets.hh
+++ b/include/HEJ/Zjets.hh
@@ -1,115 +1,123 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in Z+Jets.
*/
#pragma once
#include <vector>
#include "CLHEP/Vector/LorentzVector.h"
#include "HEJ/PDG_codes.hh"
typedef CLHEP::HepLorentzVector HLV;
namespace HEJ {
class ParticleProperties;
}
//! Square of qQ->qQe+e- Z+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 pep Momentum of final state positron
* @param pem Momentum of final state electron
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @param zprop Mass and width of the Z boson
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns Square of the current contractions for qQ->qQe+e- Scattering
*
* This returns the square of the current contractions in qQ->qQ scattering
* with an emission of a Z Boson.
*/
-std::vector <double> ME_Z_qQ(HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem,
- HEJ::ParticleID aptype, HEJ::ParticleID bptype,
+std::vector <double> ME_Z_qQ(const HLV & pa, const HLV & pb,
+ const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem,
+ const HEJ::ParticleID aptype, const HEJ::ParticleID bptype,
HEJ::ParticleProperties const & zprop,
- double stw2, double ctw);
+ const double stw2, const double ctw);
//! Square of qg->qge+e- Z+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 pep Momentum of final state positron
* @param pem Momentum of final state electron
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @param zprop Mass and width of the Z boson
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns Square of the current contractions for qg->qge+e- Scattering
*
* This returns the square of the current contractions in qg->qg scattering
* with an emission of a Z Boson.
*/
-double ME_Z_qg(HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem,
- HEJ::ParticleID aptype, HEJ::ParticleID bptype,
+double ME_Z_qg(const HLV & pa, const HLV & pb,
+ const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem,
+ const HEJ::ParticleID aptype, const HEJ::ParticleID bptype,
HEJ::ParticleProperties const & zprop,
- double stw2, double ctw);
+ const double stw2, const double ctw);
//! Square of qQ->gqQe+e- Z+Jets Unordered Current
/**
* @param pa Momentum of initial state quark a
* @param pb Momentum of initial state quark b
* @param pg Momentum of final state unordered gluon
* @param p1 Momentum of final state quark a
* @param p2 Momentum of final state quark b
* @param pep Momentum of final state positron
* @param pem Momentum of final state electron
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @param zprop Mass and width of the Z boson
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns Square of the current contractions for qQ->gqQe+e- Scattering
*
* This returns the square of the current contractions in qQ->gqQ scattering
* with an emission of a Z Boson.
*/
-std::vector <double> ME_Zuno_qQ(HLV pa, HLV pb, HLV pg, HLV p1, HLV p2, HLV pep, HLV pem,
- HEJ::ParticleID aptype, HEJ::ParticleID bptype,
+std::vector <double> ME_Zuno_qQ(const HLV & pa, const HLV & pb,
+ const HLV & pg, const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem,
+ const HEJ::ParticleID aptype, const HEJ::ParticleID bptype,
HEJ::ParticleProperties const & zprop,
- double stw2, double ctw);
+ const double stw2, const double ctw);
//! Square of qg->gqge+e- Z+Jets Unordered Current
/**
* @param pa Momentum of initial state quark
* @param pb Momentum of initial state gluon
* @param pg Momentum of final state unordered gluon
* @param p1 Momentum of final state quark
* @param p2 Momentum of final state gluon
* @param pep Momentum of final state positron
* @param pem Momentum of final state electron
* @param aptype Initial particle 1 type
* @param bptype Initial particle 2 type
* @param zprop Mass and width of the Z boson
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns Square of the current contractions for qg->gqge+e- Scattering
*
* This returns the square of the current contractions in qg->gqg scattering
* with an emission of a Z Boson.
*/
-double ME_Zuno_qg(HLV pa, HLV pb, HLV pg, HLV p1, HLV p2, HLV pep, HLV pem,
- HEJ::ParticleID aptype, HEJ::ParticleID bptype,
+double ME_Zuno_qg(const HLV & pa, const HLV & pb,
+ const HLV & pg, const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem,
+ const HEJ::ParticleID aptype, const HEJ::ParticleID bptype,
HEJ::ParticleProperties const & zprop,
- double stw2, double ctw);
\ No newline at end of file
+ const double stw2, const double ctw);
\ No newline at end of file
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index ed6acf7..e4be323 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,2203 +1,2204 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/jets.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/utility.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Zjets.hh"
namespace HEJ{
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*log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
return (
1. + alpha_s/(4.*M_PI)*beta0*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(size_t i=0; i<tree_kin_part.size(); ++i) {
sum += tree_kin_part.at(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(! is_resummable(event.type())) {
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;
}
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;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
#ifndef NDEBUG
bool wc = true;
#endif
bool wqq = false;
// With extremal qqx 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::qqxexb) {
q -= partons[1].p;
++first_idx;
if (abs(partons[0].type) != abs(partons[1].type)){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else {
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
if (in[0].type != partons[0].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(abs(backquark->type) != abs((backquark+1)->type)) {
wqq=true;
#ifndef NDEBUG
wc=false;
#endif
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqx) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
#ifndef NDEBUG
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
#endif
return exp(exponent);
}
std::vector <double> MatrixElement::virtual_corrections_Z_qq(
Event const & event,
- double mur,
+ const double mur,
Particle const & ZBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p;
fastjet::PseudoJet q_b = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q_t -= partons[1].p;
q_b -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum_top=0., sum_bot=0., sum_mix=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
- double dy = partons[j+1].rapidity() - partons[j].rapidity();
- double tmp_top = omega0(alpha_s, mur, q_t)*dy;
- double tmp_bot = omega0(alpha_s, mur, q_b)*dy;
+ const double dy = partons[j+1].rapidity() - partons[j].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 -= partons[j+1].p;
q_b -= partons[j+1].p;
}
return {exp(sum_top), exp(sum_bot), exp(sum_mix)};
}
double MatrixElement::virtual_corrections_Z_qg(
Event const & event,
- double mur,
+ const double mur,
Particle const & ZBoson,
- bool is_gq_event
+ const bool is_gq_event
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
// If this is a gq event, don't subtract the Z momentum from first q
fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p);
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity()
- partons[j].rapidity());
q -= partons[j+1].p;
}
return exp(sum);
}
std::vector<double> MatrixElement::virtual_corrections(
Event const & event,
const double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
const auto AWZH_boson = std::find_if(
begin(out), end(out),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson != end(out) && abs(AWZH_boson->type) == pid::Wp){
return {virtual_corrections_W(event, mur, *AWZH_boson)};
}
if(AWZH_boson != end(out) && AWZH_boson->type == pid::Z_photon_mix){
if(is_gluon(in.back().type)){
// This is a qg event
return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, false)};
}else if(is_gluon(in.front().type)){
// This is a gq event
return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, true)};
}else{
// This is a qq event
return virtual_corrections_Z_qq(event, mur, *AWZH_boson);
}
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
size_t first_idx = 0;
size_t last_idx = out.size() - 1;
// if there is a Higgs boson, extremal qqx 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)
|| event.type() == event_type::unob
|| event.type() == event_type::qqxexb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return {0.};
last_idx = std::distance(begin(out), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqx) q -= out[last_idx+1].p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return {exp(exponent)};
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
const CLHEP::HepLorentzVector temptrans=-(qav+qbv);
const CLHEP::HepLorentzVector p5=qav-qbv;
const CLHEP::HepLorentzVector CL=temptrans
+ 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.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);
}
//! 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 temptrans=-(qav+qbv);
const CLHEP::HepLorentzVector p5=qav-qbv;
const CLHEP::HepLorentzVector CL=temptrans
+ 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.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_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
- HLV temptrans_b = -(qav_b + qbv_b);
- HLV temptrans_t = -(qav_t + qbv_t);
+ const HLV temptrans_b = -(qav_b + qbv_b);
+ const HLV temptrans_t = -(qav_t + qbv_t);
- HLV p5_b = qav_b - qbv_b;
- HLV p5_t = qav_t - qbv_t;
+ const HLV p5_b = qav_b - qbv_b;
+ const HLV p5_t = qav_t - qbv_t;
- HLV CL_b = -(qav_b + qbv_b)
- + (qav_b.m2() / p5_b.dot(p1) + 2.0 * p5_b.dot(p2) / p1.dot(p2)) * p1
- - p2 * (qbv_b.m2() / p5_b.dot(p2) + 2.0 * p5_b.dot(p1) / p1.dot(p2));
- HLV CL_t = -(qav_t + qbv_t)
- + (qav_t.m2() / p5_t.dot(p1) + 2.0 * p5_t.dot(p2) / p1.dot(p2)) * p1
- - p2 * (qbv_t.m2() / p5_t.dot(p2) + 2.0 * p5_t.dot(p1) / p1.dot(p2));
+ const HLV CL_b = -(qav_b + qbv_b)
+ + (qav_b.m2() / p5_b.dot(p1) + 2.0 * p5_b.dot(p2) / p1.dot(p2)) * p1
+ - p2 * (qbv_b.m2() / p5_b.dot(p2) + 2.0 * p5_b.dot(p1) / p1.dot(p2));
+ const HLV CL_t = -(qav_t + qbv_t)
+ + (qav_t.m2() / p5_t.dot(p1) + 2.0 * p5_t.dot(p2) / p1.dot(p2)) * p1
+ - p2 * (qbv_t.m2() / p5_t.dot(p2) + 2.0 * p5_t.dot(p1) / p1.dot(p2));
return -CL_b.dot(CL_t);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
- double lambda
+ const double lambda
) {
- double kperp = (qav_b - qbv_b).perp();
+ const double kperp = (qav_b - qbv_b).perp();
if (kperp > lambda) {
return C2Lipatov_Mix(qav_b, qbv_b, qav_t, qbv_t, p1, p2)
/ sqrt(qav_b.m2() * qbv_b.m2() * qav_t.m2() * qbv_t.m2());
} else {
- double Cls = C2Lipatov_Mix(qav_b, qbv_b, qav_t, qbv_t, p1, p2)
- / sqrt(qav_b.m2() * qbv_b.m2() * qav_t.m2() * qbv_t.m2());
+ const double Cls = C2Lipatov_Mix(qav_b, qbv_b, qav_t, qbv_t, p1, p2)
+ / sqrt(qav_b.m2() * qbv_b.m2() * qav_t.m2() * qbv_t.m2());
return Cls - 4.0 / (kperp * kperp);
}
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
) {
- HLV temptrans_b = -(qav_b + qbv_b);
- HLV temptrans_t = -(qav_t + qbv_t);
-
- HLV p5_b = qav_b - qbv_b;
- HLV p5_t = qav_t - qbv_t;
-
- HLV CL_b = temptrans_b
- + qav_b.m2()*(1./p5_b.dot(pip)*pip+1./p5_b.dot(pop)*pop)/2.
- - qbv_b.m2()*(1./p5_b.dot(pim)*pim+1./p5_b.dot(pom)*pom)/2.
- + (pip*(p5_b.dot(pim)/pip.dot(pim)+p5_b.dot(pom)/pip.dot(pom))
- +pop*(p5_b.dot(pim)/pop.dot(pim)+p5_b.dot(pom)/pop.dot(pom))
- -pim*(p5_b.dot(pip)/pip.dot(pim) + p5_b.dot(pop)/pop.dot(pim))
- -pom*(p5_b.dot(pip)/pip.dot(pom) + p5_b.dot(pop)/pop.dot(pom)))/2.;
-
- HLV CL_t = temptrans_t
- + qav_t.m2()*(1./p5_t.dot(pip)*pip+1./p5_t.dot(pop)*pop)/2.
- - qbv_t.m2()*(1./p5_t.dot(pim)*pim+1./p5_t.dot(pom)*pom)/2.
- + (pip*(p5_t.dot(pim)/pip.dot(pim)+p5_t.dot(pom)/pip.dot(pom))
- +pop*(p5_t.dot(pim)/pop.dot(pim)+p5_t.dot(pom)/pop.dot(pom))
- -pim*(p5_t.dot(pip)/pip.dot(pim) + p5_t.dot(pop)/pop.dot(pim))
- -pom*(p5_t.dot(pip)/pip.dot(pom) + p5_t.dot(pop)/pop.dot(pom)))/2.;
+ const HLV temptrans_b = -(qav_b + qbv_b);
+ const HLV temptrans_t = -(qav_t + qbv_t);
+
+ const HLV p5_b = qav_b - qbv_b;
+ const HLV p5_t = qav_t - qbv_t;
+
+ const HLV CL_b = temptrans_b
+ + qav_b.m2()*(1./p5_b.dot(pip)*pip+1./p5_b.dot(pop)*pop)/2.
+ - qbv_b.m2()*(1./p5_b.dot(pim)*pim+1./p5_b.dot(pom)*pom)/2.
+ + (pip*(p5_b.dot(pim)/pip.dot(pim)+p5_b.dot(pom)/pip.dot(pom))
+ +pop*(p5_b.dot(pim)/pop.dot(pim)+p5_b.dot(pom)/pop.dot(pom))
+ -pim*(p5_b.dot(pip)/pip.dot(pim) + p5_b.dot(pop)/pop.dot(pim))
+ -pom*(p5_b.dot(pip)/pip.dot(pom) + p5_b.dot(pop)/pop.dot(pom)))/2.;
+
+ const HLV CL_t = temptrans_t
+ + qav_t.m2()*(1./p5_t.dot(pip)*pip+1./p5_t.dot(pop)*pop)/2.
+ - qbv_t.m2()*(1./p5_t.dot(pim)*pim+1./p5_t.dot(pom)*pom)/2.
+ + (pip*(p5_t.dot(pim)/pip.dot(pim)+p5_t.dot(pom)/pip.dot(pom))
+ +pop*(p5_t.dot(pim)/pop.dot(pim)+p5_t.dot(pom)/pop.dot(pom))
+ -pim*(p5_t.dot(pip)/pip.dot(pim) + p5_t.dot(pop)/pop.dot(pim))
+ -pom*(p5_t.dot(pip)/pip.dot(pom) + p5_t.dot(pop)/pop.dot(pom)))/2.;
return -CL_b.dot(CL_t);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
- double lambda
+ const double lambda
) {
- double kperp = (qav_b - qbv_b).perp();
+ const double kperp = (qav_b - qbv_b).perp();
if (kperp > lambda) {
return C2Lipatov_Mix(qav_b, qbv_b, qav_t, qbv_t, pa, pb, p1, p2)
/ sqrt(qav_b.m2() * qbv_b.m2() * qav_t.m2() * qbv_t.m2());
} else {
- double Cls = C2Lipatov_Mix(qav_b, qbv_b, qav_t, qbv_t, pa, pb, p1, p2)
- / sqrt(qav_b.m2() * qbv_b.m2() * qav_t.m2() * qbv_t.m2());
- double temp = Cls - 4.0 / (kperp * kperp);
+ const double Cls = C2Lipatov_Mix(qav_b, qbv_b, qav_t, qbv_t, pa, pb, p1, p2)
+ / sqrt(qav_b.m2() * qbv_b.m2() * qav_t.m2() * qbv_t.m2());
+ const double temp = Cls - 4.0 / (kperp * kperp);
return temp;
}
}
/** 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(
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
if (bptype==pid::gluon) {
if (is_quark(aptype))
return ME_unob_qg(pg,p1,pa,pn,pb);
else
return ME_unob_qbarg(pg,p1,pa,pn,pb);
}
else if (is_antiquark(bptype)) {
if (is_quark(aptype))
return ME_unob_qQbar(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQbar(pg,p1,pa,pn,pb);
}
else { //bptype == quark
if (is_quark(aptype))
return ME_unob_qQ(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQ(pg,p1,pa,pn,pb);
}
throw std::logic_error("unreachable");
}
/** Matrix element squared for tree-level current-current scattering
* @param bptype Particle b PDG ID
* @param pgin Incoming gluon momentum
* @param pq Quark from splitting Momentum
* @param pqbar Anti-quark from splitting Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param swap_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal.
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The qqxf contribution can be calculated by reversing the argument ordering.
*/
double ME_qqx_current(
ParticleID bptype,
CLHEP::HepLorentzVector const & pgin,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
bool const swap_q_qx
){
if (bptype==pid::gluon) {
if (swap_q_qx) // pq extremal
return ME_Exqqx_qqbarg(pgin,pq,pqbar,pn,pb);
else // pqbar extremal
return ME_Exqqx_qbarqg(pgin,pq,pqbar,pn,pb);
}
else { // b leg quark line
if (swap_q_qx) //extremal pq
return ME_Exqqx_qqbarQ(pgin,pq,pqbar,pn,pb);
else
return ME_Exqqx_qbarqQ(pgin,pq,pqbar,pn,pb);
}
throw std::logic_error("unreachable");
}
/* \brief Matrix element squared for central qqx 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 qqxpair
* @param nbelow Number of gluons emitted after central qqxpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_qqxmid_current(
ParticleID aptype, ParticleID bptype, int nabove,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<HLV> const & partons){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/HEJ::C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/HEJ::C_F;
return wt*ME_Cenqqx_qq(pa, pb, partons,is_antiquark(bptype),is_antiquark(aptype), swap_q_qx, nabove);
}
/** 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
){
if (aptype==pid::gluon && bptype==pid::gluon) {
return ME_gg(pn,pb,p1,pa);
} else if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_qg(pn,pb,p1,pa);
else
return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_qg(p1,pa,pn,pb);
else
return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_qQ(pn,pb,p1,pa);
else
return ME_qQbar(pn,pb,p1,pa);
}
else {
if (is_quark(aptype))
return ME_qQbar(p1,pa,pn,pb);
else
return ME_qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unreachable");
}
/** 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
){
// We know it cannot be gg incoming.
assert(!(aptype==pid::gluon && bptype==pid::gluon));
if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
}
else if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
else {
if (is_quark(aptype))
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
}
else{ // emission off a, (first argument paout)
if (is_quark(aptype)) {
if (is_quark(bptype))
return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
else { // a is anti-quark
if (is_quark(bptype))
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
}
}
throw std::logic_error("unreachable");
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_W_uno_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc, ParticleProperties const & Wprop
){
// we know they are not both gluons
assert(bptype != pid::gluon || aptype != pid::gluon);
if (bptype == pid::gluon && aptype != pid::gluon) {
// b gluon => W emission off a
if (is_quark(aptype))
return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else { // they are both quark
if (wc) {// emission off b, i.e. b is first current
if (is_quark(bptype)){
if (is_quark(aptype))
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else{
if (is_quark(aptype))
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (is_quark(aptype)) {
if (is_quark(bptype)) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else //qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else { // a is anti-quark
if (is_quark(bptype)) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else //qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
}
}
throw std::logic_error("unreachable");
}
/** \brief Matrix element squared for backward qqx tree-level current-current
* scattering With W+Jets
*
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param swap_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal.
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*
* @note calculate forwards qqx contribution by reversing argument ordering.
*/
double ME_W_qqx_current(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const swap_q_qx, bool const wc,
ParticleProperties const & Wprop
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
const double CFbackward = K_g( (swap_q_qx)?pq:pqbar ,pa)/HEJ::C_F;
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==pid::gluon && bptype==pid::gluon) {
//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission
// Site.
if (swap_q_qx)
return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)*CFbackward;
else
return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)*CFbackward;
}
else {
assert(aptype==pid::gluon && bptype!=pid::gluon );
//a gluon => W emission off b leg or qqx
if (!wc){ // W Emitted from backwards qqx
if (swap_q_qx)
return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)*CFbackward;
else
return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)*CFbackward;
}
else { // W Must be emitted from forwards leg.
if (swap_q_qx)
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, is_antiquark(bptype), Wprop)*CFbackward;
else
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop)*CFbackward;
}
}
throw std::logic_error("unreachable");
}
/* \brief Matrix element squared for central qqx 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 qqxpair
* @param nbelow Number of gluons emitted after central qqxpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqx
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_W_qqxmid_current(
ParticleID aptype, ParticleID bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<HLV> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/HEJ::C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/HEJ::C_F;
if(wqq)
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons,(is_antiquark(bptype)),(is_antiquark(aptype)),
swap_q_qx, nabove, Wprop);
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, (is_antiquark(bptype)), (is_antiquark(aptype)),
swap_q_qx, 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(
- ParticleID aptype, ParticleID bptype,
+ const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
- double stw2, double ctw
+ const double stw2, const double ctw
){
// 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
- double current_factor=ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
+ const double current_factor=ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
return { current_factor };
}else if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
- double current_factor=ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw);
+ const double current_factor=ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw);
return { current_factor };
}else if(is_anyquark(aptype) && is_anyquark(bptype)){
// This is a qq event
return ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}else{
// It cannot be gg incoming
throw std::logic_error("Bad particle types");
}
}
/** 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(
- ParticleID aptype, ParticleID bptype,
+ const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
ParticleProperties const & Zprop,
- double stw2, double ctw
+ const double stw2, const double ctw
){
// 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
- double current_factor=ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
+ const double current_factor=ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
return { current_factor };
} else if (is_gluon(aptype) && is_anyquark(bptype)) {
// This is a gq event
- double current_factor=ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw);
+ const double current_factor=ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw);
return { current_factor };
} else if (is_anyquark(aptype) && is_anyquark(bptype)) {
// This is a qq event
return ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}
throw std::logic_error("unreachable");
}
/** \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
){
if (aptype==pid::gluon && bptype==pid::gluon)
// gg initial state
return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev);
else if (aptype==pid::gluon&&bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
else
return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
}
else if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
else
return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
}
else { // they are both quark
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
else {
if (is_quark(aptype))
return ME_H_qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
}
throw std::logic_error("unreachable");
}
/** \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
){
if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
else { // they are both quark
if (is_quark(aptype)) {
if (is_quark(bptype))
return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
else {
if (is_quark(bptype))
return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
}
throw std::logic_error("unreachable");
}
CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(HEJ::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 anonymous
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(! is_resummable(ev.type())) return {0.};
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return {tree_kin_jets(ev)};
switch(AWZH_boson->type){
case pid::Higgs:
return {tree_kin_Higgs(ev)};
case pid::Wp:
case pid::Wm:
return {tree_kin_W(ev)};
case pid::Z_photon_mix:
return tree_kin_Z(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
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,
- double lambda
+ 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};
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
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;
}
const auto & jets = ev.jets();
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission or qqx
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
assert(jets.size() >= 3);
++most_backward;
}
else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = ev.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(HEJ::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;
}
namespace {
double tree_kin_jets_qqxmid(
ParticleID aptype, ParticleID bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
double lambda
){
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
const int nabove = std::distance(begin_ladder, backmidquark);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_qqxmid_current(
aptype, bptype, nabove, pa, pb,
pq, pqbar, partonsHLV
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqxt, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_qqx(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const bool swap_q_qx = is_quark(*BeginPart);
const auto pgin = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon.
const double current_factor = ME_qqx_current(
(EndIn-1)->type, pgin, pq, pqbar, pn, pb, swap_q_qx
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pgin-pq-pqbar, pgin, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const double current_factor = ME_uno_current(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pa-p1-pg, pa, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
}
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()==HEJ::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
);
}
else if (ev.type()==HEJ::event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::extremal_qqxb){
return tree_kin_jets_qqx(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::extremal_qqxf){
return tree_kin_jets_qqx(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::central_qqx){
return tree_kin_jets_qqxmid(incoming[0].type, incoming[1].type,
to_HepLorentzVector(incoming[0]),
to_HepLorentzVector(incoming[1]),
partons, param_.regulator_lambda);
}
else {
throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets");
}
}
namespace{
double tree_kin_W_FKL(
ParticleID aptype, ParticleID bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda, ParticleProperties const & Wprop
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_uno(InIter BeginIn, partIter BeginPart,
partIter EndPart, const HLV & plbar, const HLV & pl,
double lambda, ParticleProperties const & Wprop){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w
auto q0 = pa - p1 - pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_uno_current(
(BeginIn)->type, (BeginIn+1)->type, pn, pb,
p1, pa, pg, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_qqx(InIter BeginIn, partIter BeginPart,
partIter EndPart, const HLV & plbar, const HLV & pl,
double lambda, ParticleProperties const & Wprop){
const bool swap_q_qx=is_quark(*BeginPart);
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqx_current(
(BeginIn)->type, (BeginIn+1)->type, pa, pb,
pq, pqbar, pn, plbar, pl, swap_q_qx, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxmid(
ParticleID aptype, ParticleID bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda, ParticleProperties const & Wprop
){
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons.front());
auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqx emit W
bool wc = !wqq & (aptype==partons.front().type); //leg b emits w
assert(!wqq || (wqq && !wc));
if(wqq){ // emission from qqx
qqxt -= pl + plbar;
} else if(!wc) { // emission from leg a
q0 -= pl + plbar;
qqxt -= pl + plbar;
}
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
const int nabove = std::distance(begin_ladder, backmidquark);
const int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqxt, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
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
HLV plbar, 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_qqxb){
return tree_kin_W_qqx(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqx(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
assert(ev.type() == central_qqx);
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Wprop());
}
namespace{
std::vector <double> tree_kin_Z_FKL(
- ParticleID aptype, ParticleID bptype, HLV pa, HLV pb,
+ const ParticleID aptype, const ParticleID bptype,
+ HLV const & pa, HLV const & pb,
std::vector<Particle> const & partons,
- HLV plbar, HLV pl,
- double lambda, ParticleProperties const & Zprop,
- double stw2, double ctw
+ HLV const & plbar, HLV const & pl,
+ const double lambda, ParticleProperties const & Zprop,
+ const double stw2, const double ctw
){
- auto p1 = to_HepLorentzVector(partons[0]);
- auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
+ 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;
- std::vector <double> current_factor = ME_Z_current(
+ const std::vector <double> current_factor = ME_Z_current(
aptype, bptype, pn, pb, p1, pa,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
- auto q0 = pa-p1-plbar-pl;
+ 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
- auto q0 = pa-p1;
+ 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
- auto q0 = pa-p1-plbar-pl;
- auto q1 = pa-p1;
+ const auto q0 = pa-p1-plbar-pl;
+ const auto q1 = pa-p1;
ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); i++){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
template<class InIter, class partIter>
std::vector <double> tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart,
const HLV & plbar, const HLV & pl,
- double lambda, ParticleProperties const & Zprop,
- double stw2, double ctw){
+ const double lambda, ParticleProperties const & Zprop,
+ const double stw2, const double ctw){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const ParticleID aptype = (BeginIn)->type;
const ParticleID bptype = (BeginIn+1)->type;
- std::vector <double> current_factor = ME_Z_uno_current(
+ const std::vector <double> current_factor = ME_Z_uno_current(
aptype, bptype, pn, pb, p1, pa, pg,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
- auto q0 = pa-pg-p1-plbar-pl;
+ const auto q0 = pa-pg-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else if(is_gluon(aptype)){
// This is a gq event
- auto q0 = pa-pg-p1;
+ const auto q0 = pa-pg-p1;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
- auto q0 = pa-pg-p1-plbar-pl;
- auto q1 = pa-pg-p1;
+ const auto q0 = pa-pg-p1-plbar-pl;
+ const auto q1 = pa-pg-p1;
ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); i++){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
} // namespace anonymous
std::vector<double> MatrixElement::tree_kin_Z(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
// find decay products of Z
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0))
&& decay.at(0).type==-decay.at(1).type);
// get leptons
HLV plbar, 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);
- double stw2 = param_.ew_parameters.sin2_tw();
- double ctw = param_.ew_parameters.cos_tw();
+ const double stw2 = param_.ew_parameters.sin2_tw();
+ const double ctw = param_.ew_parameters.cos_tw();
if(ev.type() == FKL){
return tree_kin_Z_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_backward){
return tree_kin_Z_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_forward){
return tree_kin_Z_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets");
}
double MatrixElement::tree_kin_Higgs(Event const & ev) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev);
}
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);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
#ifdef HEJ_BUILD_WITH_QCDLOOP
// TODO: code duplication with jets.cc
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
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());
}
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == pid::gluon) return K_g(pout, pin);
return C_F;
}
#endif
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == pid::gluon)?C_A:C_F;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector const & p1out,
CLHEP::HepLorentzVector const & p1in,
ParticleID type2,
CLHEP::HepLorentzVector const & p2out,
CLHEP::HepLorentzVector const & p2in,
CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
const double vev = param_.ew_parameters.vev();
// gluon case
#ifdef HEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*ME_Houtside_gq(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH,vev) + C2gHgm(p1in,p1out,pH,vev)
)/(t1*t2);
}
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(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
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());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
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(outgoing[outgoing.size()-2].type != pid::gluon) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.back());
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.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
return MH2_forwardH(
pn, pb, incoming[0].type, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, const HLV & qH, const HLV & qHp1,
double mt, bool inc_bot, double mb, double vev){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb, vev
);
}
}
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;
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 = HEJ::C_A*HEJ::C_A/2*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 = HEJ::C_A*HEJ::C_A/2*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*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) {
const auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](auto const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())) return 1.;
switch(AWZH_boson->type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return alpha_w*alpha_w;
case pid::Z_photon_mix:
return alpha_w*alpha_w;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
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*log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ
diff --git a/src/Zjets.cc b/src/Zjets.cc
index 42d5243..76dcc40 100644
--- a/src/Zjets.cc
+++ b/src/Zjets.cc
@@ -1,509 +1,517 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include <vector>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/jets.hh"
// generated headers
#include "HEJ/currents/jZ_j.hh"
#include "HEJ/currents/jZuno_j.hh"
#include "HEJ/currents/jZ_juno.hh"
using HEJ::Helicity;
using HEJ::ParticleProperties;
using HEJ::ParticleID;
namespace helicity = HEJ::helicity;
namespace {
// Z propagator
- COM ZProp(double q, ParticleProperties const & zprop){
+ COM ZProp(const double q, ParticleProperties const & zprop){
return 1. / (q - zprop.mass*zprop.mass + COM(0.,1.)*zprop.width*zprop.mass);
}
// Photon propagator
- COM GProp(double q) {
+ COM GProp(const double q) {
return 1. / q;
}
// Weak charge
- double Zq (ParticleID PID, Helicity h, double stw2, double ctw) {
+ double Zq (const ParticleID PID, const Helicity h, const double stw2, const double ctw) {
using namespace HEJ::pid;
// Positive Spin
if (h == helicity::plus) {
// quarks
if (PID == d || PID == s || PID == b) return (+ 1.0 * stw2 / 3.0) / ctw;
if (PID == u || PID == c) return (- 2.0 * stw2 / 3.0) / ctw;
// antiquarks
if (PID == d_bar || PID == s_bar || PID == b_bar) return (+ 0.5 - 1.0 * stw2 / 3.0) / ctw;
if (PID == u_bar || PID == c_bar) return (- 0.5 + 2.0 * stw2 / 3.0) / ctw;
// electron
if (PID == electron) return stw2 / ctw;
}
// Negative Spin
else {
// quarks
if (PID == d || PID == s || PID == b) return (- 0.5 + 1.0 * stw2 / 3.0) / ctw;
if (PID == u || PID == c) return (+ 0.5 - 2.0 * stw2 / 3.0) / ctw;
// antiquarks
if (PID == d_bar || PID == s_bar || PID == b_bar) return (- 1.0 * stw2 / 3.0) / ctw;
if (PID == u_bar || PID == c_bar) return (+ 2.0 * stw2 / 3.0) / ctw;
// electron
if (PID == electron) return (-1.0 / 2.0 + stw2) / ctw;
}
throw std::logic_error("ERROR! No weak charge found");
}
// Electric charge
- double Gq (ParticleID PID) {
+ double Gq (const ParticleID PID) {
using namespace HEJ::pid;
if (PID == d || PID == s || PID == b) return -1./3.;
if (PID == u || PID == c) return +2./3.;
if (PID == d_bar || PID == s_bar || PID == b_bar) return +1./3.;
if (PID == u_bar || PID == c_bar) return -2./3.;
throw std::logic_error("ERROR! No electric charge found");
}
//! Prefactor for Z+Jets Contributions
/**
* @brief Z+Jets Contributions Prefactor
* @param aptype Incoming Particle 1 type (Z emission)
* @param PZ Z Propagator
* @param PG Photon Propagator
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @param res Ouput: pref[h1][hl]
*
* Calculates prefactor for Z+Jets (includes couplings and propagators)
*/
- void Z_amp_pref(ParticleID aptype, COM PZ, COM PG, double stw2, double ctw, COM (&res)[2][2]
+ void Z_amp_pref(const ParticleID aptype, const COM PZ, const COM PG,
+ const double stw2, const double ctw, COM (&res)[2][2]
){
using helicity::plus;
using helicity::minus;
const double Zq_a_p = Zq(aptype, plus, stw2, ctw);
const double Zq_a_m = Zq(aptype, minus, stw2, ctw);
const double Ze_p = Zq(HEJ::pid::electron, plus, stw2, ctw);
const double Ze_m = Zq(HEJ::pid::electron, minus, stw2, ctw);
const double Gq_a = Gq(aptype);
res[1][1] = -2.*(Zq_a_p * Ze_p * PZ - Gq_a * PG * stw2);
res[1][0] = -2.*(Zq_a_p * Ze_m * PZ - Gq_a * PG * stw2);
res[0][0] = -2.*(Zq_a_m * Ze_m * PZ - Gq_a * PG * stw2);
res[0][1] = -2.*(Zq_a_m * Ze_p * PZ - Gq_a * PG * stw2);
}
//! Z+Jets FKL Contribution
/**
* @brief Z+Jets FKL Contribution
* @param pa Incoming Particle 1 (Z emission)
* @param pb Incoming Particle 2
* @param p1 Outgoing Particle 1 (Z emission)
* @param p2 Outgoing Particle 2
* @param pep Outgoing positron
* @param pem Outgoing electron
* @param res Ouput: jZ_j[h1][hl][h2]
*
* Calculates j_Z^\mu j_\mu.
*/
- void jZ_j(HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, COM (&res)[2][2][2]
+ void jZ_j(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem, COM (&res)[2][2][2]
){
using helicity::plus;
using helicity::minus;
res[1][1][1] = HEJ::jZ_j<plus, plus, plus> (pa, p1, pb, p2, pem, pep);
res[1][1][0] = HEJ::jZ_j<plus, plus, minus>(pa, p1, pb, p2, pem, pep);
res[1][0][1] = HEJ::jZ_j<plus, minus, plus> (pa, p1, pb, p2, pem, pep);
res[1][0][0] = HEJ::jZ_j<plus, minus, minus>(pa, p1, pb, p2, pem, pep);
res[0][0][0] = conj(res[1][1][1]);
res[0][0][1] = conj(res[1][1][0]);
res[0][1][0] = conj(res[1][0][1]);
res[0][1][1] = conj(res[1][0][0]);
}
/**
* @brief Z+Jets Unordered Contribution, unordered on Z side
* @param pa Incoming Particle 1 (Z and Uno emission)
* @param pb Incoming Particle 2
* @param pg Unordered Gluon
* @param p1 Outgoing Particle 1 (Z and Uno emission)
* @param p2 Outgoing Particle 2
* @param pep Outgoing positron
* @param pem Outgoing electron
* @param X Ouput: (U1-L)[h1][hl][h2][hg]
* @param Y Ouput: (U2+L)[h1][hl][h2][hg]
*
* Calculates j_Z_{uno}^\mu j_\mu. Ie, unordered with Z emission same side.
*/
- void jZuno_j(HLV pa, HLV pb, HLV pg, HLV p1, HLV p2, HLV pep, HLV pem,
- COM (&X)[2][2][2][2], COM (&Y)[2][2][2][2]
+ void jZuno_j(const HLV & pa, const HLV & pb, const HLV & pg, const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem, COM (&X)[2][2][2][2], COM (&Y)[2][2][2][2]
){
using helicity::plus;
using helicity::minus;
COM U1[2][2][2][2];
U1[1][1][1][1] = HEJ::U1<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
U1[1][1][1][0] = HEJ::U1<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
U1[1][1][0][1] = HEJ::U1<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
U1[1][1][0][0] = HEJ::U1<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
U1[1][0][1][1] = HEJ::U1<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
U1[1][0][1][0] = HEJ::U1<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
U1[1][0][0][1] = HEJ::U1<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
U1[1][0][0][0] = HEJ::U1<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
U1[0][1][1][1] = HEJ::U1<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
U1[0][1][1][0] = HEJ::U1<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
U1[0][1][0][1] = HEJ::U1<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
U1[0][1][0][0] = HEJ::U1<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
U1[0][0][1][1] = HEJ::U1<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
U1[0][0][1][0] = HEJ::U1<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
U1[0][0][0][1] = HEJ::U1<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
U1[0][0][0][0] = HEJ::U1<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
COM U2[2][2][2][2];
U2[1][1][1][1] = HEJ::U2<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
U2[1][1][1][0] = HEJ::U2<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
U2[1][1][0][1] = HEJ::U2<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
U2[1][1][0][0] = HEJ::U2<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
U2[1][0][1][1] = HEJ::U2<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
U2[1][0][1][0] = HEJ::U2<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
U2[1][0][0][1] = HEJ::U2<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
U2[1][0][0][0] = HEJ::U2<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
U2[0][1][1][1] = HEJ::U2<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
U2[0][1][1][0] = HEJ::U2<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
U2[0][1][0][1] = HEJ::U2<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
U2[0][1][0][0] = HEJ::U2<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
U2[0][0][1][1] = HEJ::U2<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
U2[0][0][1][0] = HEJ::U2<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
U2[0][0][0][1] = HEJ::U2<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
U2[0][0][0][0] = HEJ::U2<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
COM L[2][2][2][2];
L[1][1][1][1] = HEJ::L<plus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
L[1][1][1][0] = HEJ::L<plus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
L[1][1][0][1] = HEJ::L<plus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
L[1][1][0][0] = HEJ::L<plus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
L[1][0][1][1] = HEJ::L<plus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
L[1][0][1][0] = HEJ::L<plus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
L[1][0][0][1] = HEJ::L<plus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
L[1][0][0][0] = HEJ::L<plus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
L[0][1][1][1] = HEJ::L<minus, plus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
L[0][1][1][0] = HEJ::L<minus, plus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
L[0][1][0][1] = HEJ::L<minus, plus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
L[0][1][0][0] = HEJ::L<minus, plus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
L[0][0][1][1] = HEJ::L<minus, minus, plus, plus> (p1, p2, pa, pb, pg, pem, pep);
L[0][0][1][0] = HEJ::L<minus, minus, plus, minus>(p1, p2, pa, pb, pg, pem, pep);
L[0][0][0][1] = HEJ::L<minus, minus, minus, plus> (p1, p2, pa, pb, pg, pem, pep);
L[0][0][0][0] = HEJ::L<minus, minus, minus, minus>(p1, p2, pa, pb, pg, pem, pep);
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
for(int hg=0; hg<2; hg++){
X[h1][hl][h2][hg] = U1[h1][hl][h2][hg] - L[h1][hl][h2][hg];
Y[h1][hl][h2][hg] = U2[h1][hl][h2][hg] + L[h1][hl][h2][hg];
}
}
}
}
}
/**
* @brief Z+Jets Unordered Contribution, unordered opposite to Z side
* @param pa Incoming Particle 1 (Z emission)
* @param pb Incoming Particle 2 (unordered emission)
* @param p1 Outgoing Particle 1 (Z emission)
* @param p2 Outgoing Particle 2 (unordered emission)
* @param pg Unordered Gluon
* @param pep Outgoing positron
* @param pem Outgoing electron
* @param X Ouput: (U1-L)[h1][hl][h2][hg]
* @param Y Ouput: (U2+L)[h1][hl][h2][hg]
*
* Calculates j_Z^\mu j_{uno}_\mu. Ie, unordered with Z emission opposite side.
*/
- void jZ_juno(HLV pa, HLV pb, HLV p1, HLV p2, HLV pg, HLV pep, HLV pem,
- COM (&X)[2][2][2][2], COM (&Y)[2][2][2][2]
+ void jZ_juno(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2, const HLV & pg,
+ const HLV & pep, const HLV & pem, COM (&X)[2][2][2][2], COM (&Y)[2][2][2][2]
){
using helicity::plus;
using helicity::minus;
COM U1[2][2][2][2];
U1[1][1][1][1] = HEJ::U1_jZ<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
U1[1][1][1][0] = HEJ::U1_jZ<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
U1[1][1][0][1] = HEJ::U1_jZ<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
U1[1][1][0][0] = HEJ::U1_jZ<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
U1[1][0][1][1] = HEJ::U1_jZ<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
U1[1][0][1][0] = HEJ::U1_jZ<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
U1[1][0][0][1] = HEJ::U1_jZ<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
U1[1][0][0][0] = HEJ::U1_jZ<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
U1[0][1][1][1] = HEJ::U1_jZ<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
U1[0][1][1][0] = HEJ::U1_jZ<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
U1[0][1][0][1] = HEJ::U1_jZ<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
U1[0][1][0][0] = HEJ::U1_jZ<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
U1[0][0][1][1] = HEJ::U1_jZ<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
U1[0][0][1][0] = HEJ::U1_jZ<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
U1[0][0][0][1] = HEJ::U1_jZ<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
U1[0][0][0][0] = HEJ::U1_jZ<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
COM U2[2][2][2][2];
U2[1][1][1][1] = HEJ::U2_jZ<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
U2[1][1][1][0] = HEJ::U2_jZ<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
U2[1][1][0][1] = HEJ::U2_jZ<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
U2[1][1][0][0] = HEJ::U2_jZ<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
U2[1][0][1][1] = HEJ::U2_jZ<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
U2[1][0][1][0] = HEJ::U2_jZ<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
U2[1][0][0][1] = HEJ::U2_jZ<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
U2[1][0][0][0] = HEJ::U2_jZ<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
U2[0][1][1][1] = HEJ::U2_jZ<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
U2[0][1][1][0] = HEJ::U2_jZ<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
U2[0][1][0][1] = HEJ::U2_jZ<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
U2[0][1][0][0] = HEJ::U2_jZ<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
U2[0][0][1][1] = HEJ::U2_jZ<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
U2[0][0][1][0] = HEJ::U2_jZ<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
U2[0][0][0][1] = HEJ::U2_jZ<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
U2[0][0][0][0] = HEJ::U2_jZ<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
COM L[2][2][2][2];
L[1][1][1][1] = HEJ::L_jZ<plus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
L[1][1][1][0] = HEJ::L_jZ<plus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
L[1][1][0][1] = HEJ::L_jZ<plus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
L[1][1][0][0] = HEJ::L_jZ<plus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
L[1][0][1][1] = HEJ::L_jZ<plus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
L[1][0][1][0] = HEJ::L_jZ<plus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
L[1][0][0][1] = HEJ::L_jZ<plus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
L[1][0][0][0] = HEJ::L_jZ<plus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
L[0][1][1][1] = HEJ::L_jZ<minus, plus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
L[0][1][1][0] = HEJ::L_jZ<minus, plus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
L[0][1][0][1] = HEJ::L_jZ<minus, plus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
L[0][1][0][0] = HEJ::L_jZ<minus, plus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
L[0][0][1][1] = HEJ::L_jZ<minus, minus, plus, plus> (pa, p1, pb, p2, pg, pem, pep);
L[0][0][1][0] = HEJ::L_jZ<minus, minus, plus, minus>(pa, p1, pb, p2, pg, pem, pep);
L[0][0][0][1] = HEJ::L_jZ<minus, minus, minus, plus> (pa, p1, pb, p2, pg, pem, pep);
L[0][0][0][0] = HEJ::L_jZ<minus, minus, minus, minus>(pa, p1, pb, p2, pg, pem, pep);
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
for(int hg=0; hg<2; hg++){
X[h1][hl][h2][hg] = U1[h1][hl][h2][hg] - L[h1][hl][h2][hg];
Y[h1][hl][h2][hg] = U2[h1][hl][h2][hg] + L[h1][hl][h2][hg];
}
}
}
}
}
} // Anonymous Namespace
-std::vector <double> ME_Z_qQ(HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem,
- ParticleID aptype, ParticleID bptype,
+std::vector <double> ME_Z_qQ(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem,
+ const ParticleID aptype, const ParticleID bptype,
ParticleProperties const & zprop,
- double stw2, double ctw
+ const double stw2, const double ctw
){
const HLV pZ = pep + pem;
const COM PZ = ZProp(pZ.m2(), zprop);
const COM PG = GProp(pZ.m2());
COM pref_top[2][2], pref_bot[2][2];
Z_amp_pref(aptype, PZ, PG, stw2, ctw, pref_top);
Z_amp_pref(bptype, PZ, PG, stw2, ctw, pref_bot);
COM Coeff_top[2][2][2], Coeff_bot[2][2][2];
jZ_j(pa, pb, p1, p2, pep, pem, Coeff_top);
jZ_j(pb, pa, p2, p1, pep, pem, Coeff_bot);
double sum_top=0., sum_bot=0., sum_mix=0.;
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
- COM res_top = pref_top[h1][hl] * Coeff_top[h1][hl][h2];
- COM res_bot = pref_bot[h2][hl] * Coeff_bot[h2][hl][h1];
+ const COM res_top = pref_top[h1][hl] * Coeff_top[h1][hl][h2];
+ const COM res_bot = pref_bot[h2][hl] * Coeff_bot[h2][hl][h1];
sum_top += norm(res_top);
sum_bot += norm(res_bot);
sum_mix += 2.0 * real(res_top * conj(res_bot));
}
}
}
const double t1_top = (pa-p1-pZ).m2();
const double t2_top = (pb-p2 ).m2();
const double t1_bot = (pa-p1 ).m2();
const double t2_bot = (pb-p2-pZ).m2();
sum_top /= t1_top * t2_top;
sum_bot /= t1_bot * t2_bot;
sum_mix /= sqrt(t1_top * t2_top * t1_bot * t2_bot);
// Colour factor: (CF*CA)/2
// Colour and helicity average: 1/(4*Nc^2)
const double pref = (HEJ::C_F*HEJ::C_A) / (8*HEJ::N_C*HEJ::N_C);
return {sum_top*pref, sum_bot*pref, sum_mix*pref};
}
-double ME_Z_qg(HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem,
- ParticleID aptype, ParticleID /*bptype*/,
+double ME_Z_qg(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem,
+ const ParticleID aptype, const ParticleID /*bptype*/,
ParticleProperties const & zprop,
- double stw2, double ctw
+ const double stw2, const double ctw
){
const HLV pZ = pep + pem;
const COM PZ = ZProp(pZ.m2(), zprop);
const COM PG = GProp(pZ.m2());
COM pref[2][2], Coeff[2][2][2];
Z_amp_pref(aptype, PZ, PG, stw2, ctw, pref);
jZ_j(pa, pb, p1, p2, pep, pem, Coeff);
double sum = 0.;
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
sum += norm(pref[h1][hl] * Coeff[h1][hl][h2]);
}
}
}
sum /= (pa-p1-pZ).m2()*(pb-p2).m2();
// Colour factor: (CF*CA)/2
// Colour and helicity average: 1/(4*Nc^2)
// Divide by CF because of gluon (K_g -> CA)
sum *= HEJ::C_A / (8*HEJ::N_C*HEJ::N_C);
// Multiply by CAM
return sum * K_g(p2, pb);
}
-std::vector <double> ME_Zuno_qQ(HLV pa, HLV pb, HLV pg, HLV p1, HLV p2, HLV pep, HLV pem,
- ParticleID aptype, ParticleID bptype,
+std::vector <double> ME_Zuno_qQ(const HLV & pa, const HLV & pb,
+ const HLV & pg, const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem,
+ const ParticleID aptype, const ParticleID bptype,
ParticleProperties const & zprop,
- double stw2, double ctw
+ const double stw2, const double ctw
){
using HEJ::C_A;
using HEJ::C_F;
const HLV pZ = pep + pem;
const COM PZ = ZProp(pZ.m2(), zprop);
const COM PG = GProp(pZ.m2());
COM prefact_top[2][2], prefact_bot[2][2];
Z_amp_pref(aptype, PZ, PG, stw2, ctw, prefact_top);
Z_amp_pref(bptype, PZ, PG, stw2, ctw, prefact_bot);
COM CoeffX_top[2][2][2][2], CoeffY_top[2][2][2][2];
jZuno_j(pa, pb, pg, p1, p2, pep, pem, CoeffX_top, CoeffY_top);
COM CoeffX_bot[2][2][2][2], CoeffY_bot[2][2][2][2];
jZ_juno(pb, pa, p2, p1, pg, pep, pem, CoeffX_bot, CoeffY_bot);
double sum_top=0., sum_bot=0., sum_mix=0.;
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
for(int hg=0; hg<2; hg++){
- COM pref_top = prefact_top[h1][hl];
- COM X_top = CoeffX_top[h1][hl][h2][hg];
- COM Y_top = CoeffY_top[h1][hl][h2][hg];
+ const COM pref_top = prefact_top[h1][hl];
+ const COM X_top = CoeffX_top[h1][hl][h2][hg];
+ const COM Y_top = CoeffY_top[h1][hl][h2][hg];
- COM pref_bot = prefact_bot[h2][hl];
- COM X_bot = CoeffX_bot[h2][hl][h1][hg];
- COM Y_bot = CoeffY_bot[h2][hl][h1][hg];
+ const COM pref_bot = prefact_bot[h2][hl];
+ const COM X_bot = CoeffX_bot[h2][hl][h1][hg];
+ const COM Y_bot = CoeffY_bot[h2][hl][h1][hg];
sum_top += norm(pref_top) * (C_A*C_F*C_F/2.*(norm(X_top)+norm(Y_top))
- C_F/2.*(X_top*conj(Y_top)).real());
sum_bot += norm(pref_bot) * (C_A*C_F*C_F/2.*(norm(X_bot)+norm(Y_bot))
- C_F/2.*(X_bot*conj(Y_bot)).real());
- COM XX = C_A*C_F*C_F/2. * pref_top * X_top * conj(pref_bot * X_bot);
- COM YY = C_A*C_F*C_F/2. * pref_top * Y_top * conj(pref_bot * Y_bot);
- COM XY = -C_F/2. * (pref_top * X_top * conj(pref_bot * Y_bot)
- + pref_top * Y_top * conj(pref_bot * X_bot));
+ const COM XX = C_A*C_F*C_F/2. * pref_top * X_top * conj(pref_bot * X_bot);
+ const COM YY = C_A*C_F*C_F/2. * pref_top * Y_top * conj(pref_bot * Y_bot);
+ const COM XY = -C_F/2. * (pref_top * X_top * conj(pref_bot * Y_bot)
+ + pref_top * Y_top * conj(pref_bot * X_bot));
sum_mix += 2.0 * real(XX + YY + XY);
}
}
}
}
const double t1_top = (pa-pg-p1-pZ).m2();
const double t2_top = (pb-p2 ).m2();
const double t1_bot = (pa-pg-p1).m2();
const double t2_bot = (pb-p2-pZ).m2();
sum_top /= t1_top * t2_top;
sum_bot /= t1_bot * t2_bot;
sum_mix /= sqrt(t1_top * t2_top * t1_bot * t2_bot);
//Helicity sum and average over initial states
const double pref = 1./(4.*C_A*C_A);
return {sum_top*pref, sum_bot*pref, sum_mix*pref};
}
-double ME_Zuno_qg(HLV pa, HLV pb, HLV pg, HLV p1, HLV p2, HLV pep, HLV pem,
- ParticleID aptype, ParticleID /*bptype*/,
+double ME_Zuno_qg(const HLV & pa, const HLV & pb,
+ const HLV & pg, const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem,
+ const ParticleID aptype, const ParticleID /*bptype*/,
ParticleProperties const & zprop,
- double stw2, double ctw
+ const double stw2, const double ctw
){
using HEJ::C_A;
using HEJ::C_F;
const HLV pZ = pep + pem;
const COM PZ = ZProp(pZ.m2(), zprop);
const COM PG = GProp(pZ.m2());
COM pref[2][2], CoeffX[2][2][2][2], CoeffY[2][2][2][2];
Z_amp_pref(aptype, PZ, PG, stw2, ctw, pref);
jZuno_j(pa, pb, pg, p1, p2, pep, pem, CoeffX, CoeffY);
double sum = 0.;
for(int h1=0; h1<2; h1++){
for(int hl=0; hl<2; hl++){
for(int h2=0; h2<2; h2++){
for(int hg=0; hg<2; hg++){
- COM X = CoeffX[h1][hl][h2][hg];
- COM Y = CoeffY[h1][hl][h2][hg];
+ const COM X = CoeffX[h1][hl][h2][hg];
+ const COM Y = CoeffY[h1][hl][h2][hg];
sum += norm(pref[h1][hl]) * (C_A*C_F*C_F/2.*(norm(X)+norm(Y))
- C_F/2.*(X*conj(Y)).real());
}
}
}
}
sum /= (pa-pg-p1-pZ).m2()*(p2-pb).m2();
//Helicity sum and average over initial states
sum /= (4.*C_A*C_A);
// Multiply by CAM
return sum * (K_g(p2, pb) / C_F);
}
\ No newline at end of file
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 4:01 PM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805029
Default Alt Text
(115 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment