Page MenuHomeHEPForge
Files F7878576

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/jets.hh b/include/HEJ/jets.hh
index 178a631..93a5381 100644
--- a/include/HEJ/jets.hh
+++ b/include/HEJ/jets.hh
@@ -1,339 +1,341 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in pure jets.
*
* This file contains all the necessary functions to compute the
* current contractions for all valid pure jet HEJ processes, which
* so far is FKL and unordered processes. It will also contain some
* pure jet ME components used in other process ME calculations
*
* @TODO add a namespace
*/
#pragma once
#include <complex>
#include <ostream>
#include <CLHEP/Vector/LorentzVector.h>
typedef std::complex<double> COM;
typedef COM current[4];
typedef CLHEP::HepLorentzVector HLV;
//! Square of qQ->qQ Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qQ Scattering
*/
double ME_qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qQbar Scattering
*
* @note this can be used for qbarQ->qbarQ Scattering by inputting arguments
* appropriately.
*/
double ME_qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*/
double ME_qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qg->qg Scattering
*
* @note this can be used for gq->gq Scattering by inputting arguments
* appropriately.
*/
double ME_qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* @note this can be used for gqbar->gqbar Scattering by inputting arguments
* appropriately.
*/
double ME_qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of gg->gg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for gg->gg Scattering
*/
double ME_gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//Unordered Backwards contributions:
//! Square of qQ->qQ Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qQ Scattering
*/
double ME_unob_qQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQ->qbarQ Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qbarQ->qbarQ Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qbarQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qQbar Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qQbar Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarQbar Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qbarQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qg Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qg->qg Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarg Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qbarg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
/** \class CCurrent jets.hh "include/HEJ/jets.hh"
* \brief This is the a new class structure for currents.
*/
class CCurrent
{
public:
CCurrent(COM sc0, COM sc1, COM sc2, COM sc3)
:c0(sc0),c1(sc1),c2(sc2),c3(sc3)
{};
CCurrent(const HLV p)
{
c0=p.e();
c1=p.px();
c2=p.py();
c3=p.pz();
};
CCurrent()
{};
CCurrent operator+(const CCurrent& other);
CCurrent operator-(const CCurrent& other);
CCurrent operator*(const double x);
CCurrent operator*(const COM x);
CCurrent operator/(const double x);
CCurrent operator/(const COM x);
friend std::ostream& operator<<(std::ostream& os, const CCurrent& cur);
COM dot(HLV p1);
COM dot(CCurrent p1);
COM c0,c1,c2,c3;
};
/* std::ostream& operator <<(std::ostream& os, const CCurrent& cur); */
CCurrent operator * ( double x, CCurrent& m);
CCurrent operator * ( COM x, CCurrent& m);
CCurrent operator / ( double x, CCurrent& m);
CCurrent operator / ( COM x, CCurrent& m);
//! Current <incoming state | mu | outgoing state>
/**
* This is a wrapper function around \see joi() note helicity flip to
* give same answer.
*/
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur);
//! Current <outgoing state | mu | outgoing state>
/**
* @param pi bra state momentum
* @param heli helicity of pi
* @param pj ket state momentum
* @param helj helicity of pj. (must be same as heli)
* @param cur reference to current which is saved.
*
* This function is for building <i (out)| mu |j (out)> currents. It
* must be called with pi as the bra, and pj as the ket.
*
* @TODO Remove heli/helj and just have helicity of current as argument.
*/
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* @param pout bra state momentum
* @param helout helicity of pout
* @param pin ket state momentum
* @param helin helicity of pin. (must be same as helout)
* @param cur reference to current which is saved.
*
* This function is for building <out| mu |in> currents. It must be
* called with pout as the bra, and pin as the ket. jio calls this
* with flipped helicity
*
* @TODO Remove helout/helin and just have helicity of current as argument.
*/
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* This is a wrapper function around the void function of the same name. \see joi
+ *
+ * @TODO This is never used
*/
CCurrent joi (HLV pout, bool helout, HLV pin, bool helin);
//! Current <incoming state | mu | outgoing state>
/**
* This is a wrapper function around the void function of the same name. \see jio
*/
CCurrent jio (HLV pout, bool helout, HLV pin, bool helin);
//! Current <outgoing state | mu | outgoing state>
/**
* This is a wrapper function around the void function of the same name. \see joo
*/
CCurrent joo (HLV pout, bool helout, HLV pin, bool helin);
inline COM cdot(const current & j1, const current & j2)
{
return j1[0]*j2[0]-j1[1]*j2[1]-j1[2]*j2[2]-j1[3]*j2[3];
}
inline COM cdot(const HLV & p, const current & j1) {
return j1[0]*p.e()-j1[1]*p.x()-j1[2]*p.y()-j1[3]*p.z();
}
inline void cmult(const COM & factor, const current & j1, current &cur)
{
cur[0]=factor*j1[0];
cur[1]=factor*j1[1];
cur[2]=factor*j1[2];
cur[3]=factor*j1[3];
}
// WHY!?!
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, const current & j5, current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0]+j4[0]+j5[0];
sum[1]=j1[1]+j2[1]+j3[1]+j4[1]+j5[1];
sum[2]=j1[2]+j2[2]+j3[2]+j4[2]+j5[2];
sum[3]=j1[3]+j2[3]+j3[3]+j4[3]+j5[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, current &sum) {
sum[0] = j1[0] + j2[0] + j3[0] + j4[0];
sum[1] = j1[1] + j2[1] + j3[1] + j4[1];
sum[2] = j1[2] + j2[2] + j3[2] + j4[2];
sum[3] = j1[3] + j2[3] + j3[3] + j4[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0];
sum[1]=j1[1]+j2[1]+j3[1];
sum[2]=j1[2]+j2[2]+j3[2];
sum[3]=j1[3]+j2[3]+j3[3];
}
inline void cadd(const current & j1, const current & j2, current &sum)
{
sum[0]=j1[0]+j2[0];
sum[1]=j1[1]+j2[1];
sum[2]=j1[2]+j2[2];
sum[3]=j1[3]+j2[3];
}
inline double abs2(const COM & a)
{
return (a*conj(a)).real();
}
inline double vabs2(const CCurrent & cur)
{
return abs2(cur.c0)-abs2(cur.c1)-abs2(cur.c2)-abs2(cur.c3);
}
inline double vre(const CCurrent & a, const CCurrent & b)
{
return real(a.c0*conj(b.c0)-a.c1*conj(b.c1)-a.c2*conj(b.c2)-a.c3*conj(b.c3));
}
//! @TODO These are not currents and should be moved elsewhere.
double K_g(double p1minus, double paminus);
double K_g(HLV const & pout, HLV const & pin);
diff --git a/src/Hjets.cc b/src/Hjets.cc
index 0f8912b..e2fa700 100644
--- a/src/Hjets.cc
+++ b/src/Hjets.cc
@@ -1,1085 +1,1084 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Hjets.hh"
#include <assert.h>
#include <limits>
#include "HEJ/Constants.hh"
#ifdef HEJ_BUILD_WITH_QCDLOOP
#include "qcdloop/qcdloop.h"
#endif
const COM looprwfactor = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4);
constexpr double infinity = std::numeric_limits<double>::infinity();
namespace {
// Loop integrals
#ifdef HEJ_BUILD_WITH_QCDLOOP
COM B0DD(HLV q, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_B0 = [](){
ql::Bubble<std::complex<double>,double,double> ql_B0;
ql_B0.setCacheSize(100);
return ql_B0;
}();
static std::vector<double> masses(2);
static std::vector<double> momenta(1);
for(auto & m: masses) m = mq*mq;
momenta.front() = q.m2();
ql_B0.integral(result, 1, masses, momenta);
return result[0];
}
COM C0DD(HLV q1, HLV q2, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_C0 = [](){
ql::Triangle<std::complex<double>,double,double> ql_C0;
ql_C0.setCacheSize(100);
return ql_C0;
}();
static std::vector<double> masses(3);
static std::vector<double> momenta(3);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = (q1+q2).m2();
ql_C0.integral(result, 1, masses, momenta);
return result[0];
}
COM D0DD(HLV q1,HLV q2, HLV q3, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_D0 = [](){
ql::Box<std::complex<double>,double,double> ql_D0;
ql_D0.setCacheSize(100);
return ql_D0;
}();
static std::vector<double> masses(4);
static std::vector<double> momenta(6);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = q3.m2();
momenta[3] = (q1+q2+q3).m2();
momenta[4] = (q1+q2).m2();
momenta[5] = (q2+q3).m2();
ql_D0.integral(result, 1, masses, momenta);
return result[0];
}
COM A1(HLV q1, HLV q2, double mt)
// As given in Eq. (B.2) of VDD
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
assert(mt > 0.);
mt2=mt*mt;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 4.*mt2/Delta3*(Q2-q12-q22)
-1.-4.*q12*q22/Delta3-12.*q12*q22*Q2/Delta3/Delta3*(q12+q22-Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q2,mt)-B0DD(Q,mt) )
* ( 2.*q22/Delta3+12.*q12*q22/Delta3/Delta3*(q22-q12+Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q1,mt)-B0DD(Q,mt) )
* ( 2.*q12/Delta3+12.*q12*q22/Delta3/Delta3*(q12-q22+Q2) )
- 2./Delta3/16/M_PI/M_PI*(q12+q22-Q2);
return ans;
}
COM A2(HLV q1, HLV q2, double mt)
// As given in Eq. (B.2) of VDD, but with high energy limit
// of invariants taken.
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
assert(mt > 0.);
mt2=mt*mt;
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 2.*mt2+1./2.*(q12+q22-Q2)
+2.*q12*q22*Q2/Delta3 )
+looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))
*q22*(q22-q12-Q2)/Delta3
+looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))
*q12*(q12-q22-Q2)/Delta3+1./16/M_PI/M_PI;
return ans;
}
#else // no QCDloop
COM A1(HLV, HLV, double) {
throw std::logic_error{"A1 called without QCDloop support"};
}
COM A2(HLV, HLV, double) {
throw std::logic_error{"A2 called without QCDloop support"};
}
#endif
void to_current(const HLV & q, current & ret){
ret[0]=q.e();
ret[1]=q.x();
ret[2]=q.y();
ret[3]=q.z();
}
/**
* @brief Higgs vertex contracted with current @param C1 and @param C2
*/
COM cHdot(const current & C1, const current & C2, const current & q1,
const current & q2, double mt, bool incBot, double mb, double vev)
{
if (mt == infinity) {
return (cdot(C1,C2)*cdot(q1,q2)-cdot(C1,q2)*cdot(C2,q1))/(3*M_PI*vev);
}
else {
HLV vq1,vq2;
vq1.set(q1[1].real(),q1[2].real(),q1[3].real(),q1[0].real());
vq2.set(q2[1].real(),q2[2].real(),q2[3].real(),q2[0].real());
// first minus sign obtained because of q1-difference to VDD
// Factor is because 4 mt^2 g^2/vev A1 -> 16 pi mt^2/vev alphas,
if(!(incBot))
return 16.*M_PI*mt*mt/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)
-cdot(C1,C2)*A2(-vq1,vq2,mt));
else
return 16.*M_PI*mt*mt/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)
-cdot(C1,C2)*A2(-vq1,vq2,mt))
+ 16.*M_PI*mb*mb/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)
-cdot(C1,C2)*A2(-vq1,vq2,mb));
}
}
//@{
/**
* @brief Higgs+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param q1 t-channel momenta into higgs vertex
* @param q2 t-channel momenta out of higgs vertex
* @param mt top mass (inf or value)
* @param incBot Bool, to include bottom mass (true) or not (false)?
* @param mb bottom mass (value)
* @param pg Unordered Gluon momenta
- * @param aqlineb Bool. Is Backwards quark line an anti-quark line?
- * @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j^\mu H j_\mu. FKL with higgs vertex somewhere in the FKL chain.
* Handles all possible incoming states.
*/
- double j_h_j(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV q1, HLV q2, double mt, bool incBot,
- double mb, double vev, bool aqlineb, bool aqlinef){
+ double j_h_j(HLV const & p1out, HLV const & p1in, HLV const & p2out,
+ HLV const & p2in, HLV const & q1, HLV const & q2,
+ double mt, bool incBot, double mb, double vev
+ ){
current j1p,j1m,j2p,j2m, q1v, q2v;
// Note need to flip helicities in anti-quark case.
- joi(p1out,!aqlineb, p1in,!aqlineb, j1p);
- joi(p1out, aqlineb, p1in, aqlineb, j1m);
- joi(p2out,!aqlinef, p2in,!aqlinef, j2p);
- joi(p2out, aqlinef, p2in, aqlinef, j2m);
+ joi(p1out, false, p1in, false, j1p);
+ joi(p1out, true, p1in, true, j1m);
+ joi(p2out, false, p2in, false, j2p);
+ joi(p2out, true, p2in, true, j2m);
to_current(q1, q1v);
to_current(q2, q2v);
COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb, vev);
COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb, vev);
COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb, vev);
COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb, vev);
// average over helicities
const double sst=(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm))/4.;
return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
}
} // namespace anonymous
double ME_H_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev, false, false);
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev, false, true);
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev, true, false);
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev, true, true);
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev, false, false)
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/HEJ::C_A;
}
double ME_H_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev, true, false)
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/HEJ::C_A;
}
double ME_H_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb, double vev){
- return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev, false, false)
+ return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/HEJ::C_A * K_g(p1out,p1in)/HEJ::C_A;
}
//@}
namespace {
//@{
/// @brief Higgs vertex contracted with one current
- CCurrent jH(HLV pout, bool helout, HLV pin,
- bool helin, HLV q1, HLV q2,
+ CCurrent jH(HLV const & pout, bool helout, HLV const & pin,
+ bool helin, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev)
{
CCurrent j2 = joi(pout,helout,pin,helin);
CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());
if(mt == infinity)
return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*vev);
else
{
if(incBot)
return (-16.*M_PI*mb*mb/vev*j2.dot(q1)*jq2*A1(-q1,q2,mb)
-16.*M_PI*mb*mb/vev*j2*A2(-q1,q2,mb))
+ (-16.*M_PI*mt*mt/vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)
-16.*M_PI*mt*mt/vev*j2*A2(-q1,q2,mt));
else
return (-16.*M_PI*mt*mt/vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)
-16.*M_PI*mt*mt/vev*j2*A2(-q1,q2,mt));
}
}
//@}
//@{
/**
* @brief Higgs+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (W emission)
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param q1 t-channel momenta into higgs vertex
* @param q2 t-channel momenta out of higgs vertex
* @param mt top mass (inf or value)
* @param incBot Bool, to include bottom mass (true) or not (false)?
* @param mb bottom mass (value)
- * @param aqlineb Bool. Is Backwards quark line an anti-quark line?
- * @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_{uno}^\mu H j_\mu. Unordered with higgs vertex
* somewhere in the FKL chain. Handles all possible incoming states.
*/
- double juno_h_j(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1, HLV qH2,
- double mt, bool incBot, double mb, double vev, bool aqlineb
- ){
+ double juno_h_j(HLV const & pg, HLV const & p1out, HLV const & p1in,
+ HLV const & p2out, HLV const & p2in,
+ HLV const & qH1, HLV const & qH2,
+ double mt, bool incBot, double mb, double vev
+ ){
// 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,!aqlineb, p1in,!aqlineb);
- CCurrent mj1m=joi(p1out, aqlineb, p1in, aqlineb);
- CCurrent jgap=joi(pg, !aqlineb, p1in,!aqlineb);
- CCurrent jgam=joi(pg, aqlineb, p1in, aqlineb);
+ 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, !aqlineb, pg, !aqlineb);
- CCurrent j2gm=joo(p1out, aqlineb, pg, aqlineb);
+ CCurrent j2gp=joo(p1out, false, pg, false);
+ CCurrent j2gm=joo(p1out, true, pg, true);
CCurrent mjH2p=jH(p2out, true,p2in, true,qH1,qH2,mt,incBot,mb, vev);
CCurrent mjH2m=jH(p2out,false,p2in,false,qH1,qH2,mt,incBot,mb, vev);
// Dot products of these which occur again and again
COM MHmp=mj1m.dot(mjH2p);
COM MHmm=mj1m.dot(mjH2m);
COM MHpp=mj1p.dot(mjH2p);
COM MHpm=mj1p.dot(mjH2m);
CCurrent p2o(p2out), p2i(p2in), p1o(p1out), p1i(p1in), qsum(q1+qg);
CCurrent Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmm/2.) )/q1.m2();
CCurrent Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmp/2.) )/q1.m2();
CCurrent Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpm/2.) )/q1.m2();
CCurrent Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpp/2.) )/q1.m2();
CCurrent U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2();
CCurrent U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2();
CCurrent U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2();
CCurrent U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2();
CCurrent U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2();
CCurrent U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
CCurrent U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
CCurrent U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(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)/(q2.m2()*qH2.m2());
// Now add the t-channels for the Higgs
ampsq/=qH1.m2()*qg.m2();
ampsq/=16.;
// Factor of (Cf/Ca) for each quark to match ME_H_qQ.
ampsq*=HEJ::C_F*HEJ::C_F/HEJ::C_A/HEJ::C_A;
return ampsq;
}
} // namespace anonymous
double ME_H_unob_qQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1,
HLV qH2, double mt, bool incBot, double mb, double vev){
- return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev, false);
+ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qbarQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
- return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev, true);
+ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
- return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev, false);
+ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qbarQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
- return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev, true);
+ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_gQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
- return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev, false)
+ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev)
*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_H_unob_gQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb, double vev){
- return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev, true)
+ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev)
*K_g(p2out,p2in)/HEJ::C_F;
}
//@}
// Begin finite mass stuff
#ifdef HEJ_BUILD_WITH_QCDLOOP
namespace {
// All the stuff needed for the box functions in qg->qgH now...
COM E1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2=-(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) +
S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) +
2.*(s34 + S1)*(s34 + S1)/Delta +
S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
2.*(s34 + S1)/Delta +
S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) -
C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) -
C0DD(k1 + k2, q2, mq)*2.*s34/
S1 - (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 +
S1)/(S1*Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(2.*s34*(s34 +
S1)*(S1 - S2)/(Delta*Sigma) +
2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S1)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM F1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-S2*D0DD(k1, k2, q2,
mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) -
s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
S2*pow((s34 + S2),2)/Delta/Delta)
- (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12)
- C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2,
mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 +
S2)*(S2 - S1)/(Delta*Sigma) + (B0DD(
k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S2)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM G1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor*(S2*D0DD(k1, q2, k2,
mq)*(Delta/s12/s12 - 4.*mq*mq/s12) -
S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) -
C0DD(k1, q2, mq))*4.*mq*mq/
s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./
s12 + (B0DD(k1 + q2, mq) -
B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1)));
}
COM E4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k2, k1, q2,
mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) -
s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S1),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta +
2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) +
s12*(S1 - S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))
*((2.*s12*(2.*s12*s34 - S1*(S1 + S2) + s34*(S2 - S1)))/(Delta*Sigma)));
}
COM F4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k1, k2, q2,
mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) +
s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S2),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 +
S2)/Delta + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 -
S1*(S1 + S2) +
s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))
*(2.*s12*(2.*s12*s34 - S2*(S1 + S2) + s34*(S1 - S2))/(Delta*Sigma)));
}
COM G4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2,
mq)*(Delta/s12 + (s12 + S1)/2. -
4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD(
k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD(
k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*2./(s12 + S2));
}
COM E10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta +
12.*mq*mq*S1*(s34 + S1)/Delta/Delta -
4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) +
C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) -
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) +
8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s34*(4.*s12 + 3.*S1 +
S2)/(Delta*Sigma) +
8.*s12*s34*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2,
mq))*(12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*S1*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM F10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (s12*D0DD(k1, k2, q2,
mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta +
12.*mq*mq*s34*(s12 + S1)/Delta/Delta -
4.*s12*pow((s34 + S2),2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*(s34 + S2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) -
C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) +
4.*s12*s34*(S2 - S1)/(Delta*Sigma) +
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) - (B0DD(
k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) +
8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(-12*s34*(2*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s12*s34*s34/(S2*Delta*Delta) +
4.*s34*S1/(Delta*Sigma) -
4.*s34*(s12*s34*(2.*s12 + S2) -
S1*S1*(2.*s12 +
S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*(2.*s34 + S1)/(Delta*Sigma) -
8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) +
s12*(S2 -
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM G10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. +
4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2,
k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*4.*(s34 +
S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) -
B0DD(k2 + q2, mq))*4.*s34/(Delta*S2));
}
COM H1(HLV k1, HLV k2, HLV kh, double mq){
return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq);
}
COM H4(HLV k1, HLV k2, HLV kh, double mq){
return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq);
}
COM H10(HLV k1, HLV k2, HLV kh, double mq){
return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq);
}
COM H2(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H1(k2,k1,kh,mq);
}
COM H5(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H4(k2,k1,kh,mq);
}
COM H12(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H10(k2,k1,kh,mq);
}
// FL and FT functions
COM FL(HLV q1, HLV q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*((2.-
3.*q1.m2()*q2.dot(Q)/detQ2)*(B0DD(q1, mq) -
B0DD(Q, mq)) + (2. -
3.*q2.m2()*q1.dot(Q)/detQ2)*(B0DD(q2, mq) -
B0DD(Q, mq)) - (4.*mq*mq + q1.m2() + q2.m2() +
Q.m2() - 3.*q1.m2()*q2.m2()*Q.m2()/detQ2)*C0DD(
q1, q2, mq) - 2.);
}
COM FT(HLV q1, HLV q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*(Q.m2()*(B0DD(q1, mq) + B0DD(q2, mq) - 2.*B0DD(Q, mq) -
2.*q1.dot(q2)*C0DD(q1, q2, mq)) + (q1.m2() -
q2.m2()) *(B0DD(q1, mq) - B0DD(q2, mq))) -
q1.dot(q2)*FL(q1, q2, mq);
}
HLV ParityFlip(HLV p){
HLV flippedVector;
flippedVector.setE(p.e());
flippedVector.setX(-p.x());
flippedVector.setY(-p.y());
flippedVector.setZ(-p.z());
return flippedVector;
}
/// @brief HC amp for qg->qgH with finite top (i.e. j^{++}_H)
void g_gH_HC(HLV pa, HLV p1,
HLV pH, double mq, double vev, current &retAns)
{
current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1,
epsHapart2,conjepsH1part1,conjepsH1part2;
COM ang1a,sqa1;
const double F = 4.*mq*mq/vev;
// Easier to have the whole thing as current object so I can use cdot functionality.
// Means I need to write pa,p1 as current objects
to_current(pa, pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
//HEJ gauge
jio(pa,false,p1,false,cura1);
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
} else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(-1./sqrt(2)/ang1a,cura1,conjeps1);
cmult(1./sqrt(2)/sqa1,cura1,epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM h4 = H4(p1,-pa,pH,mq);
const COM h5 = H5(p1,-pa,pH,mq);
const COM h10 = H10(p1,-pa,pH,mq);
const COM h12 = H12(p1,-pa,pH,mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
const COM aH1 = cdot(pHcur, cura1);
current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/ang1a, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/ang1a, epsHa, T2);
}
cmult(sqrt(2.)/ang1a*aH1, epsHa, T3);
cmult(sqrt(2.)/sqa1*aH1, conjepsH1, T4);
cmult(-sqrt(2.)*Fta*pa.dot(p1)*aH1/sqa1, conjeps1, T5);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/ang1a, epsa, T6);
cmult(-aH1/sqrt(2.)/sqa1*h4*8.*COM(0.,1.)*M_PI*M_PI, conjeps1, T7);
cmult(aH1/sqrt(2.)/ang1a*h5*8.*COM(0.,1.)*M_PI*M_PI, epsa, T8);
cmult(aH1*aH1/2./ang1a/sqa1*h10*8.*COM(0.,1.)*M_PI*M_PI, pacur, T9);
cmult(-aH1*aH1/2./ang1a/sqa1*h12*8.*COM(0.,1.)*M_PI*M_PI, p1cur, T10);
current ans;
for(int i=0;i<4;i++)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5[i]+T6[i]+T7[i]+T8[i]+T9[i]+T10[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
/// @brief HNC amp for qg->qgH with finite top (i.e. j^{+-}_H)
void g_gH_HNC(HLV pa, HLV p1, HLV pH, double mq, double vev, current &retAns)
{
const double F = 4.*mq*mq/vev;
COM ang1a,sqa1;
current conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur,
p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1,
conjepsH1part2;
// Find here if pa, meaning the gluon, is forward or backward
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
jio(pa,true,p1,true,cura1);
joi(p1,true,pa,true,cur1a);
to_current(pa,pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
to_current(pa+p1,paplusp1cur);
to_current(p1-pa,p1minuspacur);
const COM aH1 = cdot(pHcur,cura1);
const COM oneHa = std::conj(aH1); // = cdot(pHcur,cur1a)
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
}
else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(1./sqrt(2)/sqa1, cur1a, epsa);
cmult(-1./sqrt(2)/sqa1, cura1, conjeps1);
const COM phase = cdot(conjeps1, epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM Falpha = FT(p1-pa,pa-p1-pH,mq);
const COM Fbeta = FL(p1-pa,pa-p1-pH,mq);
const COM h1 = H1(p1,-pa, pH, mq);
const COM h2 = H2(p1,-pa, pH, mq);
const COM h4 = H4(p1,-pa, pH, mq);
const COM h5 = H5(p1,-pa, pH, mq);
const COM h10 = H10(p1,-pa, pH, mq);
const COM h12 = H12(p1,-pa, pH, mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
current T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a,
T11b,T12a,T12b,T13;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/sqa1, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())*((p1.x()+COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, epsHa, T2);
}
const COM boxdiagFact = 8.*COM(0.,1.)*M_PI*M_PI;
cmult(aH1*sqrt(2.)/sqa1, epsHa, T3);
cmult(oneHa*sqrt(2.)/sqa1, conjepsH1, T4);
cmult(-2.*phase*Fta*pa.dot(pH), p1cur, T5a);
cmult(2.*phase*Ft1*p1.dot(pH), pacur, T5b);
cmult(-sqrt(2.)*Fta*p1.dot(pa)*oneHa/sqa1, conjeps1, T6);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/sqa1, epsa, T7);
cmult(-boxdiagFact*phase*h2, pacur, T8a);
cmult(boxdiagFact*phase*h1, p1cur, T8b);
cmult(boxdiagFact*aH1/sqrt(2.)/sqa1*h5, epsa, T9);
cmult(-boxdiagFact*oneHa/sqrt(2.)/sqa1*h4, conjeps1, T10);
cmult(boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h10, pacur, T11a);
cmult(-boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h12, p1cur, T11b);
cmult(-phase/(pa-p1).m2()*Falpha*(p1-pa).dot(pa-p1-pH), paplusp1cur, T12a);
cmult(phase/(pa-p1).m2()*Falpha*(pa+p1).dot(pa-p1-pH), p1minuspacur, T12b);
cmult(-phase*Fbeta*(pa-p1-pH).m2(), paplusp1cur, T13);
current ans;
for(int i=0;i<4;i++)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]
+T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
} // namespace anonymous
// JDC - new amplitude with Higgs emitted close to gluon with full mt effects.
// Keep usual HEJ-style function call
double ME_Houtside_gq(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pH,
double mq, bool includeBottom, double mq2, double vev
){
current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip;
current retAns,retAnsb;
joi(p2out,true,p2in,true,cur2bplus);
joi(p2out,false,p2in,false,cur2bminus);
joi(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip);
joi(ParityFlip(p2out),false,ParityFlip(p2in),false,cur2bminusFlip);
COM app1,app2,apm1,apm2;
COM app3, app4, apm3, apm4;
if(!includeBottom)
{
g_gH_HC(p1in,p1out,pH,mq,vev,retAns);
app1=cdot(retAns,cur2bplus);
app2=cdot(retAns,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
app3=cdot(retAns,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,vev,retAns);
apm1=cdot(retAns,cur2bplus);
apm2=cdot(retAns,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
apm3=cdot(retAns,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip);
} else {
g_gH_HC(p1in,p1out,pH,mq,vev,retAns);
g_gH_HC(p1in,p1out,pH,mq2,vev,retAnsb);
app1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
app2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,vev,retAnsb);
app3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,vev,retAns);
g_gH_HNC(p1in,p1out,pH,mq2,vev,retAnsb);
apm1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
apm2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,vev,retAnsb);
apm3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
}
return abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1)
+ abs2(apm2) + abs2(apm3) + abs2(apm4);
}
#endif // HEJ_BUILD_WITH_QCDLOOP
double C2gHgm(HLV p2, HLV p1, HLV pH, double vev)
{
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
double s12,p1p,p2p;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p1p=p1.plus();
p2p=p2.plus();
} else { // opposite case
p1p=p1.minus();
p2p=p2.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=COM(0,1)*A/(2.*s12)*(p2p/p1p*conj(p1perp)*p3perp+p1p/p2p*p1perp*conj(p3perp));
temp=temp*conj(temp);
return temp.real();
}
double C2gHgp(HLV p2, HLV p1, HLV pH, double vev)
{
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.23) in hep-ph/0301013
double s12,php,p1p,phm;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
php=pH.plus();
phm=pH.minus();
p1p=p1.plus();
} else { // opposite case
php=pH.minus();
phm=pH.plus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=-COM(0,1)*A/(2.*s12)*( conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p)
+s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm)
-pow(conj(p3perp)
+(1.+php/p1p)*conj(p1perp),2)/((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH)))) );
temp=temp*conj(temp);
return temp.real();
}
double C2qHqm(HLV p2, HLV p1, HLV pH, double vev)
{
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.22) in hep-ph/0301013
double s12,p2p,p1p;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p2p=p2.plus();
p1p=p1.plus();
} else { // opposite case
p2p=p2.minus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=A/(2.*s12)*( sqrt(p2p/p1p)*p3perp*conj(p1perp)
+sqrt(p1p/p2p)*p1perp*conj(p3perp) );
temp=temp*conj(temp);
return temp.real();
}
diff --git a/src/jets.cc b/src/jets.cc
index c97ca19..4a07c63 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,369 +1,372 @@
/**
* \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
- * @param aqlineb Bool. Is Backwards quark line an anti-quark line?
- * @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_\mu j^\mu.
- * Handles all possible incoming states.
+ * Handles all possible incoming states. Helicity doesn't matter since we sum
+ * over all of them.
*/
- double j_j(HLV p1out, HLV p1in, HLV p2out, HLV p2in, bool aqlineb, bool aqlinef){
- HLV q1=p1in-p1out;
- HLV q2=-(p2in-p2out);
+ 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,!aqlineb, p1in,!aqlineb, mj1p);
- joi(p1out, aqlineb, p1in, aqlineb, mj1m);
- joi(p2out,!aqlinef, p2in,!aqlinef, mj2p);
- joi(p2out, aqlinef, p2in, aqlinef, mj2m);
+ 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 Mmp=cdot(mj1m,mj2p);
- COM Mmm=cdot(mj1m,mj2m);
- COM Mpp=cdot(mj1p,mj2p);
- COM Mpm=cdot(mj1p,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 sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);
+ 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, false, false);
+ 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, false, true);
+ 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, true, true);
+ 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, false, false)*K_g(p2out, p2in)/HEJ::C_F;
+ 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, true, false)*K_g(p2out, p2in)/(HEJ::C_F);
+ 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, false, false)*K_g(p1out, p1in)*K_g(p2out, p2in)/(HEJ::C_F*HEJ::C_F);
+ 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 pg, HLV p1out, HLV p1in, HLV p2out,
- HLV p2in, bool aqlinea, bool aqlineb){
+ 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,!aqlinea,p1in,!aqlinea);
- CCurrent mj1m=joi(p1out,aqlinea,p1in,aqlinea);
- CCurrent jgap=joi(pg,!aqlinea,p1in,!aqlinea);
- CCurrent jgam=joi(pg,aqlinea,p1in,aqlinea);
+ 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,!aqlinea,pg,!aqlinea);
- CCurrent j2gm=joo(p1out,aqlinea,pg,aqlinea);
+ CCurrent j2gp=joo(p1out, false, pg, false);
+ CCurrent j2gm=joo(p1out, true, pg, true);
- CCurrent mj2p=joi(p2out,!aqlineb,p2in,!aqlineb);
- CCurrent mj2m=joi(p2out,aqlineb,p2in,aqlineb);
+ 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, false, false);
+ 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, true, false);
+ 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, false, true);
+ 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, true, true);
+ 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, false, false)*K_g(p2out,p2in)/HEJ::C_F;
+ 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, true, false)*K_g(p2out,p2in)/HEJ::C_F;
+ return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/HEJ::C_F;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:22 PM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3802744
Default Alt Text
(68 KB)

Event Timeline