diff --git a/src/Zjets.cc b/src/Zjets.cc
index 86c591d..61148d1 100644
--- a/src/Zjets.cc
+++ b/src/Zjets.cc
@@ -1,509 +1,509 @@
/**
* \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){
return 1. / (q - zprop.mass*zprop.mass + COM(0.,1.)*zprop.width*zprop.mass);
}
// Photon propagator
COM GProp(double q) {
return 1. / q;
}
// Weak charge
double Zq (int PID, Helicity h, double stw2, double ctw) {
// Positive Spin
if (h == helicity::plus) {
// quarks
if (PID == 1 || PID == 3 || PID == 5) return (+ 1.0 * stw2 / 3.0) / ctw;
if (PID == 2 || PID == 4) return (- 2.0 * stw2 / 3.0) / ctw;
if (PID == -1 || PID == -3 || PID == -5) return (- 1.0 * stw2 / 3.0) / ctw;
if (PID == -2 || PID == -4) return (+ 2.0 * stw2 / 3.0) / ctw;
// electron
if (PID == 11) return stw2 / ctw;
}
// Negative Spin
else {
// quarks
if (PID == 1 || PID == 3 || PID == 5) return (-0.5 + 1.0 * stw2 / 3.0) / ctw;
if (PID == 2 || PID == 4) return ( 0.5 - 2.0 * stw2 / 3.0) / ctw;
if (PID == -1 || PID == -3 || PID == -5) return ( 0.5 - 1.0 * stw2 / 3.0) / ctw;
if (PID == -2 || PID == -4) return (-0.5 + 2.0 * stw2 / 3.0) / ctw;
// electron
if (PID == 11) return (-1.0 / 2.0 + stw2) / ctw;
}
throw std::logic_error("ERROR! No weak charge found");
}
// Electric charge
double Gq (int PID) {
if (PID == 1 || PID == 3 || PID == 5) return -1.0 / 3.0;
if (PID == 2 || PID == 4) return 2.0 / 3.0;
if (PID == -1 || PID == -3 || PID == -5) return 1.0 / 3.0;
if (PID == -2 || PID == -4) return -2.0 / 3.0;
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]
){
using helicity::plus;
using helicity::minus;
const double Zq_a_p = Zq(aptype, (aptype>0 ? plus : minus), stw2, ctw);
const double Zq_a_m = Zq(aptype, (aptype>0 ? minus : plus), 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]
){
using helicity::plus;
using helicity::minus;
const COM jZ_j_ppp = HEJ::jZ_j<plus, plus, plus> (pa, p1, pb, p2, pem, pep);
const COM jZ_j_ppm = HEJ::jZ_j<plus, plus, minus>(pa, p1, pb, p2, pem, pep);
const COM jZ_j_pmp = HEJ::jZ_j<plus, minus, plus> (pa, p1, pb, p2, pem, pep);
const COM jZ_j_pmm = HEJ::jZ_j<plus, minus, minus>(pa, p1, pb, p2, pem, pep);
res[1][1][1] = jZ_j_ppp;
res[1][1][0] = jZ_j_ppm;
res[1][0][1] = jZ_j_pmp;
res[1][0][0] = jZ_j_pmm;
res[0][0][0] = conj(jZ_j_ppp);
res[0][0][1] = conj(jZ_j_ppm);
res[0][1][0] = conj(jZ_j_pmp);
res[0][1][1] = conj(jZ_j_pmm);
}
/**
* @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]
){
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]
){
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,
ParticleProperties const & zprop,
double stw2, 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];
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*/,
ParticleProperties const & zprop,
double stw2, 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,
ParticleProperties const & zprop,
double stw2, 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];
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];
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));
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,
+ ParticleID aptype, ParticleID /*bptype*/,
ParticleProperties const & zprop,
double stw2, 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];
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