Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F9501436
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
40 KB
Subscribers
None
View Options
diff --git a/src/Hjets.cc b/src/Hjets.cc
index 4ae5158..0f8912b 100644
--- a/src/Hjets.cc
+++ b/src/Hjets.cc
@@ -1,1083 +1,1085 @@
/**
* \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){
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);
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);
}
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);
}
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);
}
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);
}
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)
* 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)
* 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)
* 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,
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
){
// 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);
// Note for function joo(): <p1+|pg+> = <pg-|p1->.
CCurrent j2gp=joo(p1out, !aqlineb, pg, !aqlineb);
CCurrent j2gm=joo(p1out, aqlineb, pg, aqlineb);
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);
}
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);
}
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);
}
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);
}
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)*K_g(p2out,p2in)/HEJ::C_F;
+ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev, false)
+ *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)*K_g(p2out,p2in)/HEJ::C_F;
+ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev, true)
+ *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();
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:27 PM (8 h, 44 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4484519
Default Alt Text
(40 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment