Page MenuHomeHEPForge
Files F8309879

No OneTemporary
Actions

View Options

diff --git a/src/Zjets.cc b/src/Zjets.cc
index 0100208..2e22a47 100644
--- a/src/Zjets.cc
+++ b/src/Zjets.cc
@@ -1,480 +1,488 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include <array>
#include <vector>
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/jets.hh"
// generated headers
#include "HEJ/currents/jV_j.hh"
#include "HEJ/currents/jVuno_j.hh"
#include "HEJ/currents/jV_juno.hh"
namespace HEJ {
namespace currents {
// helper for multidimensional std::array, for example
// Array<T, N1, N2> = std::array<std::array<T, N1>, N2>
template<typename T, std::size_t N, std::size_t... Ns>
struct ArrayTag{
using type = typename ArrayTag<std::array<T, N>, Ns...>::type;
};
template<typename T, std::size_t N>
struct ArrayTag<T, N> {
using type = std::array<T, N>;
};
template<typename T, std::size_t N, std::size_t... Ns>
using Array = typename ArrayTag<T, N, Ns...>::type;
namespace {
// Z propagator
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(const double q) {
return 1. / q;
}
// Weak charge
template<Helicity h>
double Zq(ParticleID PID, double stw2, double ctw);
// Positive Spin
template<>
double Zq<helicity::plus>(
const ParticleID PID,
const double stw2,
const double ctw
) {
using namespace pid;
// 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;
throw std::logic_error("ERROR! No weak charge found");
}
// Negative Spin
template<>
double Zq<helicity::minus>(
const ParticleID PID,
const double stw2,
const double ctw
) {
using namespace pid;
// 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 (const ParticleID PID) {
using namespace 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 propZ Z Propagator
* @param propG 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)
*/
Array<COM, 2, 2> Z_amp_pref(
const ParticleID aptype,
COM const & propZ, COM const & propG,
const double stw2, const double ctw
){
using helicity::plus;
using helicity::minus;
const double zq_a_p = Zq<plus>(aptype, stw2, ctw);
const double zq_a_m = Zq<minus>(aptype, stw2, ctw);
const double ze_p = Zq<plus>(pid::electron, stw2, ctw);
const double ze_m = Zq<minus>(pid::electron, stw2, ctw);
const double gq_a = Gq(aptype);
Array<COM, 2, 2> res;
res[ plus][ plus] = -2.*(zq_a_p * ze_p * propZ - gq_a * propG * stw2);
res[ plus][minus] = -2.*(zq_a_p * ze_m * propZ - gq_a * propG * stw2);
res[minus][minus] = -2.*(zq_a_m * ze_m * propZ - gq_a * propG * stw2);
res[minus][ plus] = -2.*(zq_a_m * ze_p * propZ - gq_a * propG * stw2);
return res;
}
//! 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(const HLV & pa, const HLV & pb, const HLV & p1, const HLV & p2,
- const HLV & pep, const HLV & pem, COM (&res)[2][2][2]
+ Array<COM, 2, 2, 2> jZ_j(
+ const HLV & pa, const HLV & pb,
+ const HLV & p1, const HLV & p2,
+ const HLV & pep, const HLV & pem
){
using helicity::plus;
using helicity::minus;
- res[1][1][1] = jV_j<plus, plus, plus> (pa, p1, pb, p2, pem, pep);
- res[1][1][0] = jV_j<plus, plus, minus>(pa, p1, pb, p2, pem, pep);
- res[1][0][1] = jV_j<plus, minus, plus> (pa, p1, pb, p2, pem, pep);
- res[1][0][0] = jV_j<plus, minus, minus>(pa, p1, pb, p2, pem, pep);
+ Array<COM, 2, 2, 2> res;
+#define ASSIGN_HEL(RES, J, H1, HL, H2) \
+ RES[H1][HL][H2] = J<H1, HL, H2>(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]);
+ ASSIGN_HEL(res, jV_j, plus, minus, minus);
+ ASSIGN_HEL(res, jV_j, plus, minus, plus);
+ ASSIGN_HEL(res, jV_j, plus, plus, minus);
+ ASSIGN_HEL(res, jV_j, plus, plus, plus);
+
+#undef ASSIGN_HEL
+
+ for(auto hl: {minus, plus}) {
+ for(auto h2: {minus, plus}) {
+ res[minus][hl][h2] = std::conj(res[plus][flip(hl)][flip(h2)]);
+ }
+ }
+ return res;
}
// X and Y as used in contractions with unordered currents
struct XY {
COM X;
COM Y;
};
/**
* @brief Z+Jets Unordered Contribution, unordered on Z side
* @tparam h1 Helicity of line 1 (Z emission line)
* @tparam hl Lepton Helicity
* @tparam h2 Helicity of line 2
* @tparam hg Helicity of unordered gluon
* @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
* @returns X: (U1-L), Y: (U2+l)
*
* Calculates j_Z_{uno}^\mu j_\mu. Ie, unordered with Z emission same side.
*/
template<Helicity h1, Helicity hl, Helicity h2, Helicity hg>
XY amp_jZuno_j(
const HLV & pa, const HLV & pb, const HLV & pg,
const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem
){
const COM u1 = U1<h1, hl, h2, hg>(p1, p2, pa, pb, pg, pem, pep);
const COM u2 = U2<h1, hl, h2, hg>(p1, p2, pa, pb, pg, pem, pep);
const COM l = L <h1, hl, h2, hg>(p1, p2, pa, pb, pg, pem, pep);
return {u1 - l, u2 + l};
}
Array<XY, 2, 2, 2, 2> jZuno_j(
const HLV & pa, const HLV & pb, const HLV & pg,
const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem
){
using helicity::plus;
using helicity::minus;
Array<XY, 2, 2, 2, 2> xy;
#define ASSIGN_HEL(XY, J, H1, H2, H3, H4) \
XY[H1][H2][H3][H4] = J<H1, H2, H3, H4>(pa, pb, pg, p1, p2, pep, pem)
ASSIGN_HEL(xy, amp_jZuno_j, minus, minus, minus, minus);
ASSIGN_HEL(xy, amp_jZuno_j, minus, minus, minus, plus);
ASSIGN_HEL(xy, amp_jZuno_j, minus, minus, plus, minus);
ASSIGN_HEL(xy, amp_jZuno_j, minus, minus, plus, plus);
ASSIGN_HEL(xy, amp_jZuno_j, minus, plus, minus, minus);
ASSIGN_HEL(xy, amp_jZuno_j, minus, plus, minus, plus);
ASSIGN_HEL(xy, amp_jZuno_j, minus, plus, plus, minus);
ASSIGN_HEL(xy, amp_jZuno_j, minus, plus, plus, plus);
ASSIGN_HEL(xy, amp_jZuno_j, plus, minus, minus, minus);
ASSIGN_HEL(xy, amp_jZuno_j, plus, minus, minus, plus);
ASSIGN_HEL(xy, amp_jZuno_j, plus, minus, plus, minus);
ASSIGN_HEL(xy, amp_jZuno_j, plus, minus, plus, plus);
ASSIGN_HEL(xy, amp_jZuno_j, plus, plus, minus, minus);
ASSIGN_HEL(xy, amp_jZuno_j, plus, plus, minus, plus);
ASSIGN_HEL(xy, amp_jZuno_j, plus, plus, plus, minus);
ASSIGN_HEL(xy, amp_jZuno_j, plus, plus, plus, plus);
#undef ASSIGN_HEL
return xy;
}
/**
* @brief Z+Jets Unordered Contribution, unordered opposite to Z side
* @tparam h1 Helicity of line 1 (Z emission)
* @tparam hl Lepton Helicity
* @tparam h2 Helicity of line 2 (unordered emission)
* @tparam hg Helicity of unordered gluon
* @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
* @returns X: (U1-L), Y: (U2+l)
*
* Calculates j_Z^\mu j_{uno}_\mu. Ie, unordered with Z emission opposite side.
*/
template<Helicity h1, Helicity hl, Helicity h2, Helicity hg>
XY amp_jZ_juno(
const HLV & pa, const HLV & pb,
const HLV & p1, const HLV & p2, const HLV & pg,
const HLV & pep, const HLV & pem
){
const COM u1 = U1_jV<h1, hl, h2, hg>(pa, p1, pb, p2, pg, pem, pep);
const COM u2 = U2_jV<h1, hl, h2, hg>(pa, p1, pb, p2, pg, pem, pep);
const COM l = L_jV <h1, hl, h2, hg>(pa, p1, pb, p2, pg, pem, pep);
return {u1 - l, u2 + l};
}
Array<XY, 2, 2, 2, 2> jZ_juno(
const HLV & pa, const HLV & pb, const HLV & pg,
const HLV & p1, const HLV & p2,
const HLV & pep, const HLV & pem
){
using helicity::plus;
using helicity::minus;
Array<XY, 2, 2, 2, 2> xy;
#define ASSIGN_HEL(XY, J, H1, H2, H3, H4) \
XY[H1][H2][H3][H4] = J<H1, H2, H3, H4>(pa, pb, pg, p1, p2, pep, pem)
ASSIGN_HEL(xy, amp_jZ_juno, minus, minus, minus, minus);
ASSIGN_HEL(xy, amp_jZ_juno, minus, minus, minus, plus);
ASSIGN_HEL(xy, amp_jZ_juno, minus, minus, plus, minus);
ASSIGN_HEL(xy, amp_jZ_juno, minus, minus, plus, plus);
ASSIGN_HEL(xy, amp_jZ_juno, minus, plus, minus, minus);
ASSIGN_HEL(xy, amp_jZ_juno, minus, plus, minus, plus);
ASSIGN_HEL(xy, amp_jZ_juno, minus, plus, plus, minus);
ASSIGN_HEL(xy, amp_jZ_juno, minus, plus, plus, plus);
ASSIGN_HEL(xy, amp_jZ_juno, plus, minus, minus, minus);
ASSIGN_HEL(xy, amp_jZ_juno, plus, minus, minus, plus);
ASSIGN_HEL(xy, amp_jZ_juno, plus, minus, plus, minus);
ASSIGN_HEL(xy, amp_jZ_juno, plus, minus, plus, plus);
ASSIGN_HEL(xy, amp_jZ_juno, plus, plus, minus, minus);
ASSIGN_HEL(xy, amp_jZ_juno, plus, plus, minus, plus);
ASSIGN_HEL(xy, amp_jZ_juno, plus, plus, plus, minus);
ASSIGN_HEL(xy, amp_jZ_juno, plus, plus, plus, plus);
#undef ASSIGN_HEL
return xy;
}
} // Anonymous Namespace
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,
const double stw2, const double ctw
){
const HLV pZ = pep + pem;
const COM propZ = ZProp(pZ.m2(), zprop);
const COM propG = GProp(pZ.m2());
Array<COM, 2, 2> pref_top = Z_amp_pref(aptype, propZ, propG, stw2, ctw);
Array<COM, 2, 2> pref_bot = Z_amp_pref(bptype, propZ, propG, stw2, ctw);
- 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);
+ Array<COM, 2, 2, 2> coeff_top = jZ_j(pa, pb, p1, p2, pep, pem);
+ Array<COM, 2, 2, 2> coeff_bot = jZ_j(pb, pa, p2, p1, pep, pem);
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){
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 = (C_F*C_A) / (8*N_C*N_C);
return {sum_top*pref, sum_bot*pref, sum_mix*pref};
}
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,
const double stw2, const double ctw
){
const HLV pZ = pep + pem;
const COM propZ = ZProp(pZ.m2(), zprop);
const COM propG = GProp(pZ.m2());
- COM coeff[2][2][2];
Array<COM, 2, 2> pref = Z_amp_pref(aptype, propZ, propG, stw2, ctw);
- jZ_j(pa, pb, p1, p2, pep, pem, coeff);
+ Array<COM, 2, 2, 2> coeff = jZ_j(pa, pb, p1, p2, pep, pem);
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 *= C_A / (8*N_C*N_C);
// Multiply by CAM
return sum * K_g(p2, pb);
}
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,
const double stw2, const double ctw
){
const HLV pZ = pep + pem;
const COM propZ = ZProp(pZ.m2(), zprop);
const COM propG = GProp(pZ.m2());
Array<COM, 2, 2> prefact_top = Z_amp_pref(aptype, propZ, propG, stw2, ctw);
Array<COM, 2, 2> prefact_bot = Z_amp_pref(bptype, propZ, propG, stw2, ctw);
const Array<XY, 2, 2, 2, 2> coeff_top = jZuno_j(pa, pb, pg, p1, p2, pep, pem);
const Array<XY, 2, 2, 2, 2> coeff_bot = jZ_juno(pb, pa, p2, p1, pg, pep, pem);
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){
const COM pref_top = prefact_top[h1][hl];
const COM x_top = coeff_top[h1][hl][h2][hg].X;
const COM y_top = coeff_top[h1][hl][h2][hg].Y;
const COM pref_bot = prefact_bot[h2][hl];
const COM x_bot = coeff_bot[h2][hl][h1][hg].X;
const COM y_bot = coeff_bot[h2][hl][h1][hg].Y;
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());
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(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,
const double stw2, const double ctw
){
const HLV pZ = pep + pem;
const COM propZ = ZProp(pZ.m2(), zprop);
const COM propG = GProp(pZ.m2());
Array<COM, 2, 2> pref = Z_amp_pref(aptype, propZ, propG, stw2, ctw);
const auto coeff = jZuno_j(pa, pb, pg, p1, p2, pep, pem);
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++){
const COM X = coeff[h1][hl][h2][hg].X;
const COM Y = coeff[h1][hl][h2][hg].Y;
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);
}
} // namespace currents
} // namespace HEJ

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 4:47 PM (21 h, 37 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023512
Default Alt Text
(19 KB)

Event Timeline