Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F10664332
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
28 KB
Subscribers
None
View Options
diff --git a/src/jets.cc b/src/jets.cc
index 4d41515..2ec1465 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,842 +1,830 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Tensor.hh"
#include "HEJ/Constants.hh"
// 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)*(HEJ::C_A - 1./HEJ::C_A) + 1./HEJ::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)
{
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)
{
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)
{
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)
{
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)
{
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)
{
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& m)
{
return m*x;
}
CCurrent operator * ( COM x, CCurrent& m)
{
return m*x;
}
CCurrent operator / ( double x, CCurrent& m)
{
return m/x;
}
CCurrent operator / ( COM x, CCurrent& m)
{
return m/x;
}
COM CCurrent::dot(HLV p1)
{
// 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 p1)
{
return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3;
}
//Current Functions
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
const double sqpop = sqrt(std::abs(pout.plus()));
const double sqpom = sqrt(std::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(std::abs(pin.plus()));
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * poperp / std::abs(poperp);
cur[2] = -COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(std::abs(pin.minus()));
cur[0] = -sqpom * sqpim * poperp / std::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(std::abs(pin.plus()));
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * conj(poperp) / std::abs(poperp);
cur[2] = COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
double sqpim = sqrt(std::abs(pin.minus()));
cur[0] = -sqpom * sqpim * conj(poperp) / std::abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = -COM(0,1) * cur[1];
cur[3] = -cur[0];
}
}
}
CCurrent joi (HLV pout, bool helout, HLV pin, bool helin)
{
current cur;
joi(pout, helout, pin, helin, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) {
joi(pout, !helout, pin, !helin, cur);
}
CCurrent jio (HLV pin, bool helin, HLV 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(std::abs(pj.plus() ));
const double sqpjm = sqrt(std::abs(pj.minus()));
const double sqpip = sqrt(std::abs(pi.plus() ));
const double sqpim = sqrt(std::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 / std::abs(piperp);
const COM phasej = pjperp / std::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 pi, bool heli, HLV 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 HEJ::C_F*HEJ::C_F*(sst)/(q1.m2()*q2.m2());
}
} //anonymous namespace
double ME_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/(HEJ::C_F);
}
double ME_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in)*K_g(p1out, p1in)*K_g(p2out, p2in)/(HEJ::C_F*HEJ::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=HEJ::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*=(HEJ::C_F*HEJ::C_F)/(HEJ::C_A*HEJ::C_A);
return ampsq;
}
}
//Unordered bits for pure jet
double ME_unob_qQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qbarQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qbarQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_unob_qbarg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
void eps(HLV refmom, HLV 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 pa, HLV pb, HLV p1, HLV p2, HLV p3, bool helchain,
bool heltop, bool helb,HLV 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
//Metric tensor
auto eta = HEJ::metric();
//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+= eta(i,i) *cur1a[i]*(v1[i][j])*ep[j]*eta(j,j);
}
}
return M1;
}
COM qggm2(HLV pa, HLV pb, HLV p1, HLV p2, HLV p3, bool helchain, bool heltop,
bool helb,HLV 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
auto eta=HEJ::metric();
//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+= eta(i,i)*cur1a[i]*(v2[i][j])*ep[j]*eta(j,j);
}
}
return M2;
}
COM qggm3(HLV pa, HLV pb, HLV p1, HLV p2, HLV p3, bool helchain, bool heltop,
bool helb,HLV 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()};
auto eta=HEJ::metric();
//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*eta(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+= eta(i,i)*cur1a[i]*(V3g[i][j])*ep[j]*eta(j,j);
}
}
return M3;
}
//Now the function to give helicity/colour sum/average
double MqgtqQQ(HLV pa, HLV pb, HLV p1, HLV p2, HLV 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();
}
}
// Extremal qqx
double ME_Exqqx_qbarqQ(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout);
}
double ME_Exqqx_qqbarQ(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout);
}
double ME_Exqqx_qbarqg(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout)*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_Exqqx_qqbarg(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout)*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
void CurrentMatrix(current j1, current 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 jtop, current jbot, bool hel2, HLV ka, HLV kb, HLV k2, HLV k3,
std::vector<HLV> partons, unsigned int nabove){
//Define a load of invaraints I need
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();
//Easier to have everything in terms of 'currents'
//(E,px,py,pz). To make the distinction between actual currents of
//the form ubar gamma u and 4-vectors being placed under the
//'current' class (to make dot products work out), all actual
//currents will have either a j or 'cur' in their variable name.
current cur23,curka,curkb,curk1,curk2,curk3,curk4,qc1,qc2;
//From what I gather, joo is what I use for two outgoing momenta,
//jio with one outgoing and one incoming (in that order) and j for
//lines when pa,pb are on the right of the spinor product
joo(k2,hel2,k3,hel2,cur23);
curka[0]=ka.e();
curka[1]=ka.px();
curka[2]=ka.py();
curka[3]=ka.pz();
curkb[0]=kb.e();
curkb[1]=kb.px();
curkb[2]=kb.py();
curkb[3]=kb.pz();
curk1[0]=partons.front().e();
curk1[1]=partons.front().px();
curk1[2]=partons.front().py();
curk1[3]=partons.front().pz();
curk2[0]=k2.e();
curk2[1]=k2.px();
curk2[2]=k2.py();
curk2[3]=k2.pz();
curk3[0]=k3.e();
curk3[1]=k3.px();
curk3[2]=k3.py();
curk3[3]=k3.pz();
curk4[0]=partons.back().e();
curk4[1]=partons.back().px();
curk4[2]=partons.back().py();
curk4[3]=partons.back().pz();
qc1[0]=q1.e();
qc1[1]=q1.px();
qc1[2]=q1.py();
qc1[3]=q1.pz();
qc2[0]=q2.e();
qc2[1]=q2.px();
qc2[2]=q2.py();
qc2[3]=q2.pz();
//Metric tensor
auto eta=HEJ::metric();
//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)))*eta(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))*eta(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++){
- for(int k=0; k<4; k++){
- for(int l=0; l<4;l++){
- M1+= eta(i,k)*jtop[k]*(veik[i][j]+v3g[i][j])*jbot[l]*eta(l,j);
- }
- }
+ M1+= eta(i,i)*jtop[i]*(veik[i][j]+v3g[i][j])*jbot[j]*eta(j,j);
}
}
M1/=s23;
return M1;
}
COM m2 (current jtop, current jbot, bool hel2, HLV ka, HLV k2,
HLV k3, std::vector<HLV> 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;
//Metric tensor
auto eta=HEJ::metric();
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++){
- for(int k=0; k<4; k++){
- for(int l=0; l<4;l++){
- M2+= eta(i,k)*jtop[k]*(tempv[i][j])*jbot[l]*eta(l,j);
- }
- }
+ M2+= eta(i,i)*jtop[i]*(tempv[i][j])*jbot[j]*eta(j,j);
}
}
M2/=t2;
return M2;
}
COM m3 (current jtop, current jbot, bool hel2, HLV ka, HLV k2,
HLV k3, std::vector<HLV> 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;
//Metric tensor
auto eta=HEJ::metric();
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++){
- for(int k=0; k<4; k++){
- for(int l=0; l<4;l++){
- M3+= eta(i,k)*jtop[k]*(tempv[i][j])*jbot[l]*eta(l,j);
- }
- }
+ M3+= eta(i,i)*jtop[i]*(tempv[i][j])*jbot[j]*eta(j,j);
}
}
M3/= t2t;
return M3;
}
}
double ME_Cenqqx_qq(HLV ka, HLV kb, std::vector<HLV> partons, bool aqlinepa,
bool aqlinepb, bool qqxmarker, int nabove){
//Partons in the 'wrong' ordering, so reverse it. Just put ka =
//forward and kb = backward into the function call
std::reverse(partons.begin(),partons.end());
//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());
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:36 AM (1 d, 15 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4887774
Default Alt Text
(28 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment