Page MenuHomeHEPForge
Files F19251870

jets.cc
No OneTemporary
Actions

Size
11 KB
Referenced Files
None
Subscribers
None

jets.cc
View Options

/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/jets.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(pout.plus());
const double sqpom = sqrt(pout.minus());
const COM poperp = 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(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(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(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
const double sqpim = sqrt(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 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(pj.plus());
const double sqpjm = sqrt(pj.minus());
const double sqpip = sqrt(pi.plus());
const double sqpim = sqrt(pi.minus());
const COM piperp = pi.x() + COM(0,1) * pi.y();
const COM pjperp = 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 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;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:12 AM (20 h, 56 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6490754
Default Alt Text
jets.cc (11 KB)

Event Timeline

HEPForge.org · .