Page MenuHomeHEPForge
Files F19244571

jets.cc
No OneTemporary
Actions

Size
29 KB
Referenced Files
None
Subscribers
None

jets.cc
View Options

/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include <algorithm>
#include <cmath>
#include "HEJ/Constants.hh"
namespace {
// short hand for math functions
using std::abs;
using std::conj;
using std::pow;
using std::sqrt;
double metric(const size_t mu, const size_t nu) {
if(mu != nu) return 0.;
return (mu == 0)?1.:-1.;
}
} // Anonymous Namespace
namespace HEJ {
namespace currents {
// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
// @TODO: this is not a current and should be moved somewhere else
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(
HLV const & pout,
HLV const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
CCurrent CCurrent::operator+(const CCurrent& other) const
{
COM result_c0=c0 + other.c0;
COM result_c1=c1 + other.c1;
COM result_c2=c2 + other.c2;
COM result_c3=c3 + other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator-(const CCurrent& other) const
{
COM result_c0=c0 - other.c0;
COM result_c1=c1 - other.c1;
COM result_c2=c2 - other.c2;
COM result_c3=c3 - other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const double x) const
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const double x) const
{
COM result_c0=CCurrent::c0/x;
COM result_c1=CCurrent::c1/x;
COM result_c2=CCurrent::c2/x;
COM result_c3=CCurrent::c3/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const COM x) const
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const COM x) const
{
COM result_c0=(CCurrent::c0)/x;
COM result_c1=(CCurrent::c1)/x;
COM result_c2=(CCurrent::c2)/x;
COM result_c3=(CCurrent::c3)/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
{
os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
return os;
}
CCurrent operator * ( double x, CCurrent const & m)
{
return m*x;
}
CCurrent operator * ( COM const & x, CCurrent const & m)
{
return m*x;
}
CCurrent operator / ( double x, CCurrent const & m)
{
return m/x;
}
CCurrent operator / ( COM const & x, CCurrent const & m)
{
return m/x;
}
COM CCurrent::dot(HLV const & p1) const
{
// Current goes (E,px,py,pz)
// Vector goes (px,py,pz,E)
return p1[3]*c0-p1[0]*c1-p1[1]*c2-p1[2]*c3;
}
COM CCurrent::dot(CCurrent const & p1) const
{
return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3;
}
//Current Functions
void joi(HLV const & pout, bool helout,
HLV const & pin, bool helin, current &cur
){
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
const double sqpop = sqrt(abs(pout.plus()));
const double sqpom = sqrt(abs(pout.minus()));
// Allow for "jii" format
const COM poperp = (pout.x()==0 && pout.y() ==0) ? -1 : pout.x()+COM(0,1)*pout.y();
if (helout != helin) {
throw std::invalid_argument{"Non-matching helicities"};
} else if (helout == false) { // negative helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(abs(pin.plus()));
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * poperp / abs(poperp);
cur[2] = -COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(abs(pin.minus()));
cur[0] = -sqpom * sqpim * poperp / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = COM(0,1) * cur[1];
cur[3] = -cur[0];
}
} else { // positive helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(abs(pin.plus()));
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
cur[2] = COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
double sqpim = sqrt(abs(pin.minus()));
cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = -COM(0,1) * cur[1];
cur[3] = -cur[0];
}
}
}
CCurrent joi(HLV const & pout, bool helout, HLV const & pin, bool helin)
{
current cur;
joi(pout, helout, pin, helin, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void jio(HLV const & pin, bool helin,
HLV const & pout, bool helout, current &cur
) {
joi(pout, !helout, pin, !helin, cur);
}
CCurrent jio(HLV const & pin, bool helin, HLV const & pout, bool helout){
current cur;
jio(pin, helin, pout, helout, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {
// Zero our current
cur[0] = 0.0;
cur[1] = 0.0;
cur[2] = 0.0;
cur[3] = 0.0;
if (heli!=helj) {
throw std::invalid_argument{"Non-matching helicities"};
} else if ( heli == true ) { // If positive helicity swap momenta
std::swap(pi,pj);
}
const double sqpjp = sqrt(abs(pj.plus() ));
const double sqpjm = sqrt(abs(pj.minus()));
const double sqpip = sqrt(abs(pi.plus() ));
const double sqpim = sqrt(abs(pi.minus()));
// Allow for "jii" format
const COM piperp = (pi.x()==0 && pi.y() ==0) ? -1 : pi.x()+COM(0,1)*pi.y();
const COM pjperp = (pj.x()==0 && pj.y() ==0) ? -1 : pj.x()+COM(0,1)*pj.y();
const COM phasei = piperp / abs(piperp);
const COM phasej = pjperp / abs(pjperp);
cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
}
CCurrent joo(HLV const & pi, bool heli, HLV const & pj, bool helj)
{
current cur;
joo(pi, heli, pj, helj, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
namespace {
//@{
/**
* @brief Pure Jet FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1.
* @param p1in Incoming Particle 1.
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
*
* Calculates j_\mu j^\mu.
* Handles all possible incoming states. Helicity doesn't matter since we sum
* over all of them.
*/
double j_j(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
HLV const q1=p1in-p1out;
HLV const q2=-(p2in-p2out);
current mj1m,mj1p,mj2m,mj2p;
// Note need to flip helicities in anti-quark case.
joi(p1out, false, p1in, false, mj1p);
joi(p1out, true, p1in, true, mj1m);
joi(p2out, false, p2in, false, mj2p);
joi(p2out, true, p2in, true, mj2m);
COM const Mmp=cdot(mj1m,mj2p);
COM const Mmm=cdot(mj1m,mj2m);
COM const Mpp=cdot(mj1p,mj2p);
COM const Mpm=cdot(mj1p,mj2m);
double const sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);
// Multiply by Cf^2
return C_F*C_F*(sst)/(q1.m2()*q2.m2());
}
} // Anonymous Namespace
double ME_qQ(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qbarQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/C_F;
}
double ME_qbarg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/(C_F);
}
double ME_gg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p1out, p1in)*K_g(p2out, p2in)/(C_F*C_F);
}
//@}
namespace {
double juno_j(HLV const & pg, HLV const & p1out,
HLV const & p1in, HLV const & p2out, HLV const & p2in
){
// This construction is taking rapidity order: pg > p1out >> p2out
HLV q1=p1in-p1out; // Top End
HLV q2=-(p2in-p2out); // Bottom End
HLV qg=p1in-p1out-pg; // Extra bit post-gluon
// Note <p1|eps|pa> current split into two by gauge choice.
// See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
CCurrent mj1p=joi(p1out, false, p1in, false);
CCurrent mj1m=joi(p1out, true, p1in, true);
CCurrent jgap=joi(pg, false, p1in, false);
CCurrent jgam=joi(pg, true, p1in, true);
// Note for function joo(): <p1+|pg+> = <pg-|p1->.
CCurrent j2gp=joo(p1out, false, pg, false);
CCurrent j2gm=joo(p1out, true, pg, true);
CCurrent mj2p=joi(p2out, false, p2in, false);
CCurrent mj2m=joi(p2out, true, p2in, true);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p);
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
CCurrent p1o(p1out),p2o(p2out),p2i(p2in),qsum(q1+qg),p1i(p1in);
CCurrent Lmm=(qsum*(Mmm)+(-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
CCurrent Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
CCurrent Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
CCurrent Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
CCurrent U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
CCurrent U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
CCurrent U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
CCurrent U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
CCurrent U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
CCurrent U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
CCurrent U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
CCurrent U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
constexpr double cf=C_F;
double amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
double amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
//Divide by t-channels
ampsq/=q2.m2()*qg.m2();
ampsq/=16.;
// Factor of (Cf/Ca) for each quark to match j_j.
ampsq*=(C_F*C_F)/(C_A*C_A);
return ampsq;
}
} // Anonymous Namespace
//Unordered bits for pure jet
double ME_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qg( HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F;
}
double ME_unob_qbarg(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F;
}
namespace {
void eps(HLV const & refmom, HLV const & kb, bool hel, current &ep){
//Positive helicity eps has negative helicity choices for spinors and vice versa
joi(refmom,hel,kb,hel,ep);
double norm;
if(kb.z()<0.) norm = sqrt(2.*refmom.plus()*kb.minus());
else norm = sqrt(2.*refmom.minus()*kb.plus());
// Normalise
std::for_each(begin(ep), end(ep), [&,norm](auto & val){val/=norm;});
}
COM qggm1(HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2,
HLV const & p3, bool helchain, bool heltop, bool helb, HLV const & refmom
){
// Since everything is defined with currents, need to use compeleness relation
// to expand p slash. i.e. pslash = |p><p|. Only one helicity 'survives' as
// defined by the helicities of the spinors at the end of the chain.
current cur33, cur23, curb3, cur2b, cur1a, ep;
joo(p3, helchain, p3, helchain,cur33);
joo(p2,helchain,p3,helchain,cur23);
jio(pb,helchain,p3,helchain,curb3);
joi(p2,helchain,pb,helchain,cur2b);
joi(p1, heltop, pa, heltop,cur1a);
const double t2 = (p3-pb)*(p3-pb);
//Calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
COM v1[4][4];
for(int u=0; u<4;++u){
for(int v=0; v<4;++v){
v1[u][v]=(cur23[u]*cur33[v]-cur2b[u]*curb3[v])/t2*(-1.);
}
}
//Dot in current and eps
//eps
eps(refmom,pb,helb, ep);
COM M1=0.;
// Perform Contraction: g^{ik} j_{1a}_k * v1_i^j eps^l g_lj
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M1+= metric(i,i) *cur1a[i]*(v1[i][j])*ep[j]*metric(j,j);
}
}
return M1;
}
COM qggm2(HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2,
HLV const & p3, bool helchain, bool heltop, bool helb, HLV const & refmom
){
// Since everything is defined with currents, need to use completeness relation
// to expand p slash. i.e. pslash = |p><p|. Only one helicity 'survives' as
// defined by the helicities of the spinors at the end of the chain.
current cur22, cur23, curb3, cur2b, cur1a, ep;
joo(p2, helchain, p2, helchain, cur22);
joo(p2, helchain, p3, helchain, cur23);
jio(pb, helchain, p3, helchain, curb3);
joi(p2, helchain, pb, helchain, cur2b);
joi(p1, heltop, pa, heltop, cur1a);
const double t2t = (p2-pb)*(p2-pb);
//Calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
COM v2[4][4]={};
for(int u=0; u<4;++u){
for(int v=0; v<4; ++v){
v2[u][v]=(cur22[v]*cur23[u]-cur2b[v]*curb3[u])/t2t;
}
}
//Dot in current and eps
//eps
eps(refmom,pb,helb, ep);
COM M2=0.;
// Perform Contraction: g^{ik} j_{1a}_k * v2_i^j eps^l g_lj
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M2+= metric(i,i)*cur1a[i]*(v2[i][j])*ep[j]*metric(j,j);
}
}
return M2;
}
COM qggm3(HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2,
HLV const & p3, bool helchain, bool heltop, bool helb, HLV const & refmom
){
current qqcur,ep,cur1a;
const double s23 = (p2+p3)*(p2+p3);
joo(p2,helchain,p3,helchain,qqcur);
joi(p1, heltop, pa, heltop,cur1a);
//Redefine relevant momenta as currents - for ease of calling correct part of vector
const current kb{pb.e(), pb.x(), pb.y(), pb.z()};
const current k2{p2.e(), p2.x(), p2.y(), p2.z()};
const current k3{p3.e(), p3.x(), p3.y(), p3.z()};
//Calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
COM V3g[4][4]={};
const COM kbqq=kb[0]*qqcur[0] -kb[1]*qqcur[1] -kb[2]*qqcur[2] -kb[3]*qqcur[3];
for(int u=0;u<4;++u){
for(int v=0;v<4;++v){
V3g[u][v] += 2.*COM(0.,1.)*(((k2[v]+k3[v])*qqcur[u] - (kb[u])*qqcur[v])+
kbqq*metric(u,v))/s23;
}
}
eps(refmom,pb,helb, ep);
COM M3=0.;
// Perform Contraction: g^{ik} j_{1a}_k * (v2_i^j) eps^l g_lj
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M3+= metric(i,i)*cur1a[i]*(V3g[i][j])*ep[j]*metric(j,j);
}
}
return M3;
}
//Now the function to give helicity/colour sum/average
double MqgtqQQ(HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2, HLV const & p3
){
// 4 indepedent helicity choices (complex conjugation related).
//Need to evalute each independent hel configuration and store that result somewhere
const COM Mmmm1 = qggm1(pa,pb,p1,p2,p3,false,false,false, pa);
const COM Mmmm2 = qggm2(pa,pb,p1,p2,p3,false,false,false, pa);
const COM Mmmm3 = qggm3(pa,pb,p1,p2,p3,false,false,false, pa);
const COM Mmmp1 = qggm1(pa,pb,p1,p2,p3,false,true, false, pa);
const COM Mmmp2 = qggm2(pa,pb,p1,p2,p3,false,true, false, pa);
const COM Mmmp3 = qggm3(pa,pb,p1,p2,p3,false,true, false, pa);
const COM Mpmm1 = qggm1(pa,pb,p1,p2,p3,false,false,true, pa);
const COM Mpmm2 = qggm2(pa,pb,p1,p2,p3,false,false,true, pa);
const COM Mpmm3 = qggm3(pa,pb,p1,p2,p3,false,false,true, pa);
const COM Mpmp1 = qggm1(pa,pb,p1,p2,p3,false,true, true, pa);
const COM Mpmp2 = qggm2(pa,pb,p1,p2,p3,false,true, true, pa);
const COM Mpmp3 = qggm3(pa,pb,p1,p2,p3,false,true, true, pa);
//Colour factors:
const COM cm1m1 = 8./3.;
const COM cm2m2 = 8./3.;
const COM cm3m3 = 6.;
const COM cm1m2 = -1./3.;
const COM cm1m3 = -3.*COM(0.,1.);
const COM cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
const double Mmmm = real(cm1m1*pow(abs(Mmmm1),2)+cm2m2*pow(abs(Mmmm2),2)+
cm3m3*pow(abs(Mmmm3),2)+2.*real(cm1m2*Mmmm1*conj(Mmmm2))+
2.*real(cm1m3*Mmmm1*conj(Mmmm3))+2.*real(cm2m3*Mmmm2*conj(Mmmm3)));
const double Mmmp = real(cm1m1*pow(abs(Mmmp1),2)+cm2m2*pow(abs(Mmmp2),2)+
cm3m3*pow(abs(Mmmp3),2)+2.*real(cm1m2*Mmmp1*conj(Mmmp2))+
2.*real(cm1m3*Mmmp1*conj(Mmmp3))+2.*real(cm2m3*Mmmp2*conj(Mmmp3)));
const double Mpmm = real(cm1m1*pow(abs(Mpmm1),2)+cm2m2*pow(abs(Mpmm2),2)+
cm3m3*pow(abs(Mpmm3),2)+2.*real(cm1m2*Mpmm1*conj(Mpmm2))+
2.*real(cm1m3*Mpmm1*conj(Mpmm3))+2.*real(cm2m3*Mpmm2*conj(Mpmm3)));
const double Mpmp = real(cm1m1*pow(abs(Mpmp1),2)+cm2m2*pow(abs(Mpmp2),2)+
cm3m3*pow(abs(Mpmp3),2)+2.*real(cm1m2*Mpmp1*conj(Mpmp2))+
2.*real(cm1m3*Mpmp1*conj(Mpmp3))+2.*real(cm2m3*Mpmp2*conj(Mpmp3)));
// Factor of 2 for the helicity for conjugate configurations
return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/3.)/(pa-p1).m2()/(p2+p3-pb).m2();
}
} // Anonymous Namespace
// Extremal qqx
double ME_Exqqx_qbarqQ(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout);
}
double ME_Exqqx_qqbarQ(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout);
}
double ME_Exqqx_qbarqg(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout)*K_g(p2out,p2in)/C_F;
}
double ME_Exqqx_qqbarg(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout)*K_g(p2out,p2in)/C_F;
}
namespace {
void CurrentMatrix(current const & j1, current const & j2,
COM array[4][4]
){
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
array[i][j]=j1[i]*j2[j];
}
}
}
//qqbar produced in the middle
COM m1(current const & jtop, current const & jbot, bool hel2,
HLV const & ka, HLV const & kb, HLV const & k2, HLV const & k3,
std::vector<HLV> const & partons, unsigned int nabove
){
const double s23 = 2.*(k2*k3);
const double sa2 = 2.*(ka*k2);
const double sa3 = 2.*(ka*k3);
const double s12 = 2.*(partons.front()*k2);
const double s13 = 2.*(partons.front()*k3);
const double sb2 = 2.*(kb*k2);
const double sb3 = 2.*(kb*k3);
const double s42 = 2.*(partons.back()*k2);
const double s43 = 2.*(partons.back()*k3);
HLV q1=ka-partons.front();
for(unsigned int i=1;i<nabove+1;++i)
q1-=partons.at(i);
const HLV q2=q1-partons.at(nabove+2)-partons.at(nabove+1);
const double t1 = q1.m2();
const double t3 = q2.m2();
current cur23;
joo(k2,hel2,k3,hel2,cur23);
const current curka{ka.e(),ka.px(),ka.py(),ka.pz()};
const current curkb{kb.e(),kb.px(),kb.py(),kb.pz()};
const current curk1{partons.front().e(),partons.front().px(),
partons.front().py(),partons.front().pz()};
const current curk2{k2.e(),k2.px(),k2.py(),k2.pz()};
const current curk3{k3.e(),k3.px(),k3.py(),k3.pz()};
const current curk4{partons.back().e(),partons.back().px(),
partons.back().py(),partons.back().pz()};
const current qc1{q1.e(),q1.px(),q1.py(),q1.pz()};
const current qc2{q2.e(),q2.px(),q2.py(),q2.pz()};
//Create the two bits of this vertex
COM veik[4][4],v3g[4][4];
for(int i=0;i<4;++i) {
for(int j=0;j<4;++j){
veik[i][j] = (cdot(cur23,curka)*(t1/(sa2+sa3))+cdot(cur23,curk1)*
(t1/(s12+s13))-cdot(cur23,curkb)*(t3/(sb2+sb3))-
cdot(cur23,curk4)*(t3/(s42+s43)))*metric(i,j);
}
}
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
v3g[i][j] = qc1[j]*cur23[i]+curk2[j]*cur23[i]+curk3[j]*cur23[i]+
qc2[i]*cur23[j]-curk2[i]*cur23[j]-curk3[i]*cur23[j]-
(cdot(qc1,cur23)+cdot(qc2,cur23))*metric(i,j);
}
}
//Now dot in the currents - potential problem here with Lorentz
//indicies, so check this
COM M1=0;
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M1+= metric(i,i)*jtop[i]*(veik[i][j]+v3g[i][j])*jbot[j]*metric(j,j);
}
}
M1/=s23;
return M1;
}
COM m2(current const & jtop, current const & jbot, bool hel2,
HLV const & ka, HLV k2, HLV const & k3,
std::vector<HLV> const & partons, unsigned int nabove
){
//In order to get correct momentum dependence in the vertex, forst
//have to work with CCurrent objects and then convert to 'current'
current cur22,cur23,cur2q,curq3;
COM qarray[4][4]={};
COM temp[4][4]={};
joo(k2,hel2,k2,hel2,cur22);
joo(k2,hel2,k3,hel2,cur23);
joi(k2,hel2,ka,hel2,cur2q);
jio(ka,hel2,k3,hel2,curq3);
CurrentMatrix(cur2q, curq3, qarray);
for(unsigned int i =0; i<nabove+1; ++i){
joo(k2,hel2,partons.at(i),hel2,cur2q);
joo(partons.at(i),hel2,k3,hel2,curq3);
CurrentMatrix(cur2q, curq3, temp);
for(int ii=0;ii<4;++ii){
for(int jj=0;jj<4;++jj){
qarray[ii][jj]=qarray[ii][jj]-temp[ii][jj];
}
}
}
HLV qt=ka-k2;
for(unsigned int i=0; i<nabove+1;++i){
qt-=partons.at(i);
}
const double t2=qt*qt;
COM tempv[4][4];
for(int i=0; i<4;++i){
for(int j=0;j<4;++j){
tempv[i][j] = COM(0.,1.)*(qarray[i][j]-cur22[i]*cur23[j]);
}
}
COM M2=0.;
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M2+= metric(i,i)*jtop[i]*(tempv[i][j])*jbot[j]*metric(j,j);
}
}
M2/=t2;
return M2;
}
COM m3(current const & jtop, current const & jbot, bool hel2,
HLV const & ka, HLV const & k2, HLV const & k3,
std::vector<HLV> const & partons, unsigned int nabove
){
COM M3=0.;
current cur23,cur33,cur2q,curq3;
COM qarray[4][4]={};
COM temp[4][4]={};
joo(k3,hel2,k3,hel2,cur33);
joo(k2,hel2,k3,hel2,cur23);
joi(k2,hel2,ka,hel2,cur2q);
jio(ka,hel2,k3,hel2,curq3);
CurrentMatrix(cur2q, curq3, qarray);
for(unsigned int i =0; i<nabove+1; ++i){
joo(k2,hel2,partons.at(i),hel2,cur2q);
joo(partons.at(i),hel2,k3,hel2,curq3);
CurrentMatrix(cur2q, curq3, temp);
for(int ii=0;ii<4;++ii){
for(int jj=0;jj<4;++jj){
qarray[ii][jj]=qarray[ii][jj]-temp[ii][jj];
}
}
}
HLV qt=ka-k3;
for(unsigned int i=0; i<nabove+1;++i){
qt-=partons.at(i);
}
const double t2t=qt*qt;
COM tempv[4][4];
for(int i=0; i<4;++i){
for(int j=0;j<4;++j){
tempv[i][j] = COM(0.,-1.)*(qarray[j][i]-cur23[j]*cur33[i]);
}
}
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M3+= metric(i,i)*jtop[i]*(tempv[i][j])*jbot[j]*metric(j,j);
}
}
M3/= t2t;
return M3;
}
} // Anonymous Namespace
double ME_Cenqqx_qq(HLV const & ka, HLV const & kb,
std::vector<HLV> const & partons, bool aqlinepa,
bool aqlinepb, bool qqxmarker, int nabove
){
//Get all the possible outer currents
current j1p,j1m,j4p,j4m;
if(!(aqlinepa)){
joi(partons.front(),true,ka,true,j1p);
joi(partons.front(),false,ka,false,j1m);
}
if(aqlinepa){
jio(ka,true,partons.front(),true,j1p);
jio(ka,false,partons.front(),false,j1m);
}
if(!(aqlinepb)){
joi(partons.back(),true,kb,true,j4p);
joi(partons.back(),false,kb,false,j4m);
}
if(aqlinepb){
jio(kb,true,partons.back(),true,j4p);
jio(kb,false,partons.back(),false,j4m);
}
HLV k2,k3;
if(!(qqxmarker)){
k2=partons.at(nabove+1);
k3=partons.at(nabove+2);
}
else{
k2=partons.at(nabove+2);
k3=partons.at(nabove+1);
}
//8 helicity choices we can make, but only 4 indepedent ones
//(complex conjugation related).
const COM Mmmm1 = m1(j1m,j4m,false,ka,kb,k2,k3,partons,nabove);
const COM Mmmm2 = m2(j1m,j4m,false,ka, k2,k3,partons,nabove);
const COM Mmmm3 = m3(j1m,j4m,false,ka, k2,k3,partons,nabove);
const COM Mmmp1 = m1(j1m,j4m,true, ka,kb,k2,k3,partons,nabove);
const COM Mmmp2 = m2(j1m,j4m,true, ka, k2,k3,partons,nabove);
const COM Mmmp3 = m3(j1m,j4m,true, ka, k2,k3,partons,nabove);
const COM Mpmm1 = m1(j1p,j4m,false,ka,kb,k2,k3,partons,nabove);
const COM Mpmm2 = m2(j1p,j4m,false,ka, k2,k3,partons,nabove);
const COM Mpmm3 = m3(j1p,j4m,false,ka, k2,k3,partons,nabove);
const COM Mpmp1 = m1(j1p,j4m,true, ka,kb,k2,k3,partons,nabove);
const COM Mpmp2 = m2(j1p,j4m,true, ka, k2,k3,partons,nabove);
const COM Mpmp3 = m3(j1p,j4m,true, ka, k2,k3,partons,nabove);
//Colour factors:
const COM cm1m1=3.;
const COM cm2m2=4./3.;
const COM cm3m3=4./3.;
const COM cm1m2 =3./2.*COM(0.,1.);
const COM cm1m3 = -3./2.*COM(0.,1.);
const COM cm2m3 = -1./6.;
//Square and sum for each helicity config:
const double Mmmm = real(cm1m1*pow(abs(Mmmm1),2)+cm2m2*pow(abs(Mmmm2),2)+
cm3m3*pow(abs(Mmmm3),2)+2.*real(cm1m2*Mmmm1*conj(Mmmm2))+
2.*real(cm1m3*Mmmm1*conj(Mmmm3))+2.*real(cm2m3*Mmmm2*conj(Mmmm3)));
const double Mmmp = real(cm1m1*pow(abs(Mmmp1),2)+cm2m2*pow(abs(Mmmp2),2)+
cm3m3*pow(abs(Mmmp3),2)+2.*real(cm1m2*Mmmp1*conj(Mmmp2))+
2.*real(cm1m3*Mmmp1*conj(Mmmp3))+2.*real(cm2m3*Mmmp2*conj(Mmmp3)));
const double Mpmm = real(cm1m1*pow(abs(Mpmm1),2)+cm2m2*pow(abs(Mpmm2),2)+
cm3m3*pow(abs(Mpmm3),2)+2.*real(cm1m2*Mpmm1*conj(Mpmm2))+
2.*real(cm1m3*Mpmm1*conj(Mpmm3))+2.*real(cm2m3*Mpmm2*conj(Mpmm3)));
const double Mpmp = real(cm1m1*pow(abs(Mpmp1),2)+cm2m2*pow(abs(Mpmp2),2)+
cm3m3*pow(abs(Mpmp3),2)+2.*real(cm1m2*Mpmp1*conj(Mpmp2))+
2.*real(cm1m3*Mpmp1*conj(Mpmp3))+2.*real(cm2m3*Mpmp2*conj(Mpmp3)));
//Result (averaged, without coupling or t-channel props). Factor of
//2 for the 4 helicity configurations I didn't work out explicitly
HLV prop1 = ka;
for(int i=0; i<=nabove; ++i){
prop1 -= partons[i];
}
const HLV prop2 = prop1 - k2 - k3;
return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/9./4.) /
((ka-partons.front()).m2()*(kb-partons.back()).m2()*prop1.m2()*prop2.m2());
}
} // namespace currents
} // namespace HEJ

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 4:43 AM (1 d, 6 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6480479
Default Alt Text
jets.cc (29 KB)

Event Timeline

HEPForge.org · .