diff --git a/src/Zjets.cc b/src/Zjets.cc
index bf369bf..e58a124 100644
--- a/src/Zjets.cc
+++ b/src/Zjets.cc
@@ -1,169 +1,170 @@
/**
* \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/jets.hh"
// generated headers
#include "HEJ/currents/jZ_j.hh"
using HEJ::Helicity;
using HEJ::ParticleProperties;
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");
}
void Z_amp(HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, int aptype, COM PZ, COM PG,
double stw2, double ctw, COM (&res)[2][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);
const COM pref_pp = -2.*(Zq_a_p * Ze_p * PZ - Gq_a * PG * stw2);
const COM pref_pm = -2.*(Zq_a_p * Ze_m * PZ - Gq_a * PG * stw2);
const COM pref_mm = -2.*(Zq_a_m * Ze_m * PZ - Gq_a * PG * stw2);
const COM pref_mp = -2.*(Zq_a_m * Ze_p * PZ - Gq_a * PG * stw2);
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] = (pref_pp * jZ_j_ppp);
res[1][1][0] = (pref_pp * jZ_j_ppm);
res[1][0][1] = (pref_pm * jZ_j_pmp);
res[1][0][0] = (pref_pm * jZ_j_pmm);
res[0][0][0] = (pref_mm * conj(jZ_j_ppp));
res[0][0][1] = (pref_mm * conj(jZ_j_ppm));
res[0][1][0] = (pref_mp * conj(jZ_j_pmp));
res[0][1][1] = (pref_mp * conj(jZ_j_pmm));
}
} // Anonymous Namespace
std::vector <double> ME_Z_qQ(HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem,
int aptype, int 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 Coeff_top[2][2][2], Coeff_bot[2][2][2];
Z_amp(pa, pb, p1, p2, pep, pem, aptype, PZ, PG, stw2, ctw, Coeff_top);
Z_amp(pb, pa, p2, p1, pep, pem, bptype, PZ, PG, stw2, ctw, 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++){
sum_top += norm(Coeff_top[h1][hl][h2]);
sum_bot += norm(Coeff_bot[h1][hl][h2]);
sum_mix += 2.0*real(Coeff_top[h1][hl][h2]*conj(Coeff_bot[h2][hl][h1]));
}
}
}
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);
- // Division by colour and Helicity average (Nc2-1)(4)
- // Multiply by Cf^2
- const double pref = HEJ::C_F*HEJ::C_F / ((HEJ::N_C*HEJ::N_C - 1)*4);
+ // 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,
int aptype, int /*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 Coeff[2][2][2];
Z_amp(pa, pb, p1, p2, pep, pem, aptype, PZ, PG, stw2, ctw, 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(Coeff[h1][hl][h2]);
}
}
}
sum /= (pa-p1-pZ).m2()*(pb-p2).m2();
- // Division by colour and Helicity average (Nc2-1)(4)
- // Multiply by Cf
- sum *= HEJ::C_F / ((HEJ::N_C*HEJ::N_C - 1)*4);
+ // 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);
}