currents.cc
No OneTemporary
Actions

Size

83 KB

Referenced Files

None

Subscribers

None

currents.cc
View Options

	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/currents.hh"

	#include <iostream>
	#include <limits>
	#include <utility>
	#include <vector>

	#ifdef HEJ_BUILD_WITH_QCDLOOP
	#include "qcdloop/qcdloop.h"
	#endif

	#include "HEJ/Constants.hh"
	#include "HEJ/exceptions.hh"
	#include "HEJ/PDG_codes.hh"

	const COM looprwfactor = (COM(0.,1.)M_PIM_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=q12q12+q22q22+Q2Q2-2q12q22-2q12Q2-2q22*Q2;
	if (mt < 0.)
	std::cerr<<"Problem in A1! mt = "<<mt<<std::endl;
	mt2=mt*mt;

	ans=looprwfactorCOM(0,-1)C0DD(q1,q2,mt)( 4.mt2/Delta3*(Q2-q12-q22)
	-1.-4.q12q22/Delta3-12.q12q22Q2/Delta3/Delta3(q12+q22-Q2) )
	- looprwfactorCOM(0,-1)( B0DD(q2,mt)-B0DD(Q,mt) )
	* ( 2.q22/Delta3+12.q12q22/Delta3/Delta3(q22-q12+Q2) )
	- looprwfactorCOM(0,-1)( B0DD(q1,mt)-B0DD(Q,mt) )
	* ( 2.q12/Delta3+12.q12q22/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.));

	if (mt < 0.)
	std::cerr<<"Problem in A2! mt = "<<mt<<std::endl;
	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=q12q12+q22q22+Q2Q2-2q12q22-2q12Q2-2q22*Q2;
	ans=looprwfactorCOM(0,-1)C0DD(q1,q2,mt)( 2.mt2+1./2.*(q12+q22-Q2)
	+2.q12q22*Q2/Delta3 )
	+looprwfactorCOM(0,-1)(B0DD(q2,mt)-B0DD(Q,mt))
	q22(q22-q12-Q2)/Delta3
	+looprwfactorCOM(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();
	}
	} // namespace anonymous

	// 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.c0c0-p1.c1c1-p1.c2c2-p1.c3c3;
	}

	//Current Functions

	// Current for <outgoing state \| mu \| incoming state>
	/// @TODO always use this instead of "j"
	/// @TODO isn't this jio with flipt helicities?
	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]);
	}

	// Current for <incoming state \| mu \| outgoing state>
	void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) {

	cur[0] = 0.0;
	cur[1] = 0.0;
	cur[2] = 0.0;
	cur[3] = 0.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 * 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];
	}
	} else { // positive 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];
	}
	}
	}

	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]);
	}

	// Current for <outgoing state \| mu \| outgoing state>
	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]);
	}

	double jM2qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in)
	{

	HLV q1=p1in-p1out;
	HLV q2=-(p2in-p2out);
	current mj1m,mj1p,mj2m,mj2p;
	joi(p1out,true,p1in,true,mj1p);
	joi(p1out,false,p1in,false,mj1m);
	joi(p2out,true,p2in,true,mj2p);
	joi(p2out,false,p2in,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	double sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);

	// Multiply by Cf^2
	return HEJ::C_FHEJ::C_F(sst)/(q1.m2()*q2.m2());
	}

	double jM2qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in)
	{
	HLV q1=p1in-p1out;
	HLV q2=-(p2in-p2out);
	current mj1m,mj1p,mj2m,mj2p;
	joi(p1out,true,p1in,true,mj1p);
	joi(p1out,false,p1in,false,mj1m);
	jio(p2in,true,p2out,true,mj2p);
	jio(p2in,false,p2out,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	double sumsq=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);

	// Multiply by Cf^2
	return HEJ::C_FHEJ::C_F(sumsq)/(q1.m2()*q2.m2());
	}

	double jM2qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in)
	{
	HLV q1=p1in-p1out;
	HLV q2=-(p2in-p2out);
	current mj1m,mj1p,mj2m,mj2p;
	jio(p1in,true,p1out,true,mj1p);
	jio(p1in,false,p1out,false,mj1m);
	jio(p2in,true,p2out,true,mj2p);
	jio(p2in,false,p2out,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	double sumsq=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);

	// Multiply by Cf^2
	return HEJ::C_FHEJ::C_F(sumsq)/(q1.m2()*q2.m2());
	}

	double jM2qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in)
	// Calculates the square of the current contractions for qg scattering
	// p1: quark
	// p2: gluon
	{
	HLV q1=p1in-p1out;
	HLV q2=-(p2in-p2out);

	current mj1m,mj1p,mj2m,mj2p;
	joi(p1out,true,p1in,true,mj1p);
	joi(p1out,false,p1in,false,mj1m);
	joi(p2out,true,p2in,true,mj2p);
	joi(p2out,false,p2in,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	const double K = K_g(p2out, p2in);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);
	double a2Mpp=abs2(Mpp);
	double a2Mpm=abs2(Mpm);
	double sst = K/HEJ::C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);

	return HEJ::C_FHEJ::C_Asst/(q1.m2()*q2.m2());

	}

	double jM2qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in)
	// Calculates the square of the current contractions for qg scattering
	// p1: quark
	// p2: gluon
	{
	HLV q1=p1in-p1out;
	HLV q2=-(p2in-p2out);

	current mj1m,mj1p,mj2m,mj2p;
	jio(p1in,true,p1out,true,mj1p);
	jio(p1in,false,p1out,false,mj1m);
	joi(p2out,true,p2in,true,mj2p);
	joi(p2out,false,p2in,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	const double K = K_g(p2out, p2in);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);
	double a2Mpp=abs2(Mpp);
	double a2Mpm=abs2(Mpm);
	double sst = K/HEJ::C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);

	// Cf*Ca=4
	return HEJ::C_FHEJ::C_Asst/(q1.m2()*q2.m2());

	}

	double jM2gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in)
	// Calculates the square of the current contractions for gg scattering
	// p1: gluon
	// p2: gluon
	{
	HLV q1=p1in-p1out;
	HLV q2=-(p2in-p2out);

	current mj1m,mj1p,mj2m,mj2p;
	joi(p1out,true,p1in,true,mj1p);
	joi(p1out,false,p1in,false,mj1m);
	joi(p2out,true,p2in,true,mj2p);
	joi(p2out,false,p2in,false,mj2m);

	COM Mmp=cdot(mj1m,mj2p);
	COM Mmm=cdot(mj1m,mj2m);
	COM Mpp=cdot(mj1p,mj2p);
	COM Mpm=cdot(mj1p,mj2m);

	const double K_g1 = K_g(p1out, p1in);
	const double K_g2 = K_g(p2out, p2in);

	// sum of spinor strings \|\|^2
	double a2Mmp=abs2(Mmp);
	double a2Mmm=abs2(Mmm);
	double a2Mpp=abs2(Mpp);
	double a2Mpm=abs2(Mpm);
	double sst = K_g1/HEJ::C_AK_g2/HEJ::C_A(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
	// Ca*Ca=9
	return HEJ::C_AHEJ::C_Asst/(q1.m2()*q2.m2());
	}

	namespace {
	/**
	* @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)
	{
	if (mt == infinity) {
	return (cdot(C1,C2)cdot(q1,q2)-cdot(C1,q2)cdot(C2,q1))/(6M_PIHEJ::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/HEJ::vev A1 -> 16 pi mt^2/HEJ::vev alphas,
	// and we divide by a factor 4 at the amp sqaured level later
	// which I absorb here (i.e. I divide by 2)
	/// @TODO move factor 1/2 from S to \|ME\|^2 => consistent with general notation
	if(!(incBot))
	return 8.M_PImtmt/HEJ::vev(-cdot(C1,q2)cdot(C2,q1)A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt));
	else
	return 8.M_PImtmt/HEJ::vev(-cdot(C1,q2)cdot(C2,q1)A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt))
	+ 8.M_PImbmb/HEJ::vev(-cdot(C1,q2)cdot(C2,q1)A1(-vq1,vq2,mb)-cdot(C1,C2)*A2(-vq1,vq2,mb));
	}
	}
	} // namespace anonymous

	double MH2qQ (HLV p1out, HLV p1in,
	HLV p2out, HLV p2in,
	HLV q1, HLV q2,
	double mt, bool incBot, double mb)
	{
	current j1p,j1m,j2p,j2m, q1v, q2v;

	joi (p1out,true,p1in,true,j1p);
	joi (p1out,false,p1in,false,j1m);

	joi (p2out,true,p2in,true,j2p);
	joi (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);

	double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
	double mt, bool incBot, double mb
	){
	current j1p,j1m,j2p,j2m,q1v,q2v;

	joi (p1out,true,p1in,true,j1p);
	joi (p1out,false,p1in,false,j1m);
	jio (p2in,true,p2out,true,j2p);

	jio (p2in,false,p2out,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);

	double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qbarQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
	double mt, bool incBot, double mb
	){
	current j1p,j1m,j2p,j2m,q1v,q2v;

	jio (p1in,true,p1out,true,j1p);
	jio (p1in,false,p1out,false,j1m);

	joi (p2out,true,p2in,true,j2p);
	joi (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);

	double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
	double mt, bool incBot, double mb
	){
	current j1p,j1m,j2p,j2m,q1v,q2v;

	jio (p1in,true,p1out,true,j1p);
	jio (p1in,false,p1out,false,j1m);

	jio (p2in,true,p2out,true,j2p);
	jio (p2in,false,p2out,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);

	double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
	double mt, bool incBot, double mb
	){
	// q~p1 g~p2 (i.e. ALWAYS p1 for quark, p2 for gluon)
	// should be called with q1 meant to be contracted with p2 in first part of vertex
	// (i.e. if g is backward, q1 is forward)
	current j1p,j1m,j2p,j2m,q1v,q2v;

	joi (p1out,true,p1in,true,j1p);
	joi (p1out,false,p1in,false,j1m);

	joi (p2out,true,p2in,true,j2p);
	joi (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	// First, calculate the non-flipping amplitudes:

	COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);
	COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);

	const double K = K_g(p2out, p2in);

	double sst=K/HEJ::C_A*(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm));

	return sst/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
	double mt, bool incBot, double mb
	){
	// qbar~p1 g~p2 (i.e. ALWAYS p1 for anti-quark, p2 for gluon)
	// should be called with q1 meant to be contracted with p2 in first part of vertex
	// (i.e. if g is backward, q1 is forward)
	current j1p,j1m,j2p,j2m,q1v,q2v;

	jio (p1in,true,p1out,true,j1p);
	jio (p1in,false,p1out,false,j1m);

	joi (p2out,true,p2in,true,j2p);
	joi (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	// First, calculate the non-flipping amplitudes:

	COM amp,amm,apm,app;
	app=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	apm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);
	amp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	amm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);

	double MH2sum = abs2(app)+abs2(amm)+abs2(apm)+abs2(amp);

	const double K = K_g(p2out, p2in);
	MH2sum*=K/HEJ::C_A;

	return MH2sum/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	double MH2gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
	double mt, bool incBot, double mb
	){
	// g~p1 g~p2
	// should be called with q1 meant to be contracted with p2 in first part of vertex
	// (i.e. if g is backward, q1 is forward)
	current j1p,j1m,j2p,j2m,q1v,q2v;

	joi (p1out,true,p1in,true,j1p);
	joi (p1out,false,p1in,false,j1m);

	joi (p2out,true,p2in,true,j2p);
	joi (p2out,false,p2in,false,j2m);

	to_current(q1, q1v);
	to_current(q2, q2v);

	// First, calculate the non-flipping amplitudes:

	COM amp,amm,apm,app;
	app=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
	apm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);
	amp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
	amm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);

	double MH2sum = abs2(app)+abs2(amm)+abs2(apm)+abs2(amp);

	const double K_g1 = K_g(p1out, p1in);
	const double K_g2 = K_g(p2out, p2in);

	MH2sum=K_g1/HEJ::C_AK_g2/HEJ::C_A;

	return MH2sum/((p1in-p1out).m2()(p2in-p2out).m2()q1.m2()*q2.m2());
	}

	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)
	{

	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)/(3M_PIHEJ::vev);
	else
	{
	if(incBot)
	return (-16.M_PImbmb/HEJ::vevj2.dot(q1)jq2A1(-q1,q2,mb)-16.M_PImbmb/HEJ::vevj2*A2(-q1,q2,mb))
	+ (-16.M_PImtmt/HEJ::vevj2.dot(q1)jq2A1(-q1,q2,mt)-16.M_PImtmt/HEJ::vevj2*A2(-q1,q2,mt));
	else
	return (-16.M_PImtmt/HEJ::vevj2.dot(q1)jq2A1(-q1,q2,mt)-16.M_PImtmt/HEJ::vevj2*A2(-q1,q2,mt));
	}
	}

	CCurrent jioH (HLV pin, bool helin, HLV pout,
	bool helout, HLV q1, HLV q2,
	double mt, bool incBot, double mb)
	{

	CCurrent j2 = jio(pin,helin,pout,helout);
	CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());

	if(mt == infinity)
	return ((q1.dot(q2))j2 - j2.dot(q1)jq2)/(3M_PIHEJ::vev);
	else
	{
	if(incBot)
	return (-16.M_PImbmb/HEJ::vevj2.dot(q1)jq2A1(-q1,q2,mb)-16.M_PImbmb/HEJ::vevj2*A2(-q1,q2,mb))
	+ (-16.M_PImtmt/HEJ::vevj2.dot(q1)jq2A1(-q1,q2,mt)-16.M_PImtmt/HEJ::vevj2*A2(-q1,q2,mt));
	else
	return (-16.M_PImtmt/HEJ::vevj2.dot(q1)jq2A1(-q1,q2,mt)-16.M_PImtmt/HEJ::vevj2*A2(-q1,q2,mt));
	}
	}

	CCurrent jHtop (HLV pout, bool helout, HLV pin,
	bool helin, HLV q1, HLV q2,
	double mt, bool incBot, double mb)
	{

	CCurrent j1 = joi(pout,helout,pin,helin);
	CCurrent jq1(q1.e(),q1.px(),q1.py(),q1.pz());

	if(mt == infinity)
	return ((q1.dot(q2))j1 - j1.dot(q2)jq1)/(3M_PIHEJ::vev);
	else
	{
	if(incBot)
	return (-16.M_PImbmb/HEJ::vevj1.dot(q2)jq1A1(-q1,q2,mb)-16.M_PImbmb/HEJ::vevj1*A2(-q1,q2,mb))
	+ (-16.M_PImtmt/HEJ::vevj1.dot(q2)jq1A1(-q1,q2,mt)-16.M_PImtmt/HEJ::vevj1*A2(-q1,q2,mt));
	else
	return (-16.M_PImtmt/HEJ::vevj1.dot(q2)jq1A1(-q1,q2,mt)-16.M_PImtmt/HEJ::vevj1*A2(-q1,q2,mt));
	}
	}

	CCurrent jioHtop (HLV pin, bool helin, HLV pout,
	bool helout, HLV q1, HLV q2,
	double mt, bool incBot, double mb)
	{

	CCurrent j1 = jio(pin,helin,pout,helout);
	CCurrent jq1(q1.e(),q1.px(),q1.py(),q1.pz());

	if(mt == infinity)
	return ((q1.dot(q2))j1 - j1.dot(q2)jq1)/(3M_PIHEJ::vev);
	else
	{
	if(incBot)
	return (-16.M_PImbmb/HEJ::vevj1.dot(q2)jq1A1(-q1,q2,mb)-16.M_PImbmb/HEJ::vevj1*A2(-q1,q2,mb))
	+ (-16.M_PImtmt/HEJ::vevj1.dot(q2)jq1A1(-q1,q2,mt)-16.M_PImtmt/HEJ::vevj1*A2(-q1,q2,mt));
	else
	return (-16.M_PImtmt/HEJ::vevj1.dot(q2)jq1A1(-q1,q2,mt)-16.M_PImtmt/HEJ::vevj1*A2(-q1,q2,mt));
	}
	}
	//@}
	} // namespace anonymous

	double jM2unogqHQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1,
	HLV qH2, double mt, bool incBot, double mb
	){
	// 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

	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=joi(p1out,true,p1in,true);
	mj1m=joi(p1out,false,p1in,false);
	mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(p1out,true,pg,true);
	j2gm=joo(p1out,false,pg,false);
	jgap=joi(pg,true,p1in,true);
	jgam=joi(pg,false,p1in,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	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();
	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();
	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();
	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();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	// Factor of (Cf/Ca) for each quark to match MH2qQ.
	ampsq=HEJ::C_FHEJ::C_F/HEJ::C_A/HEJ::C_A;

	return ampsq;
	}

	double jM2unogqbarHQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){
	// 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

	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=jio(p1in,true,p1out,true);
	mj1m=jio(p1in,false,p1out,false);
	mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(pg,true,p1out,true);
	j2gm=joo(pg,false,p1out,false);
	jgap=jio(p1in,true,pg,true);
	jgam=jio(p1in,false,pg,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	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();
	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();
	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();
	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();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}

	double jM2unogqHQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){
	// 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

	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=joi(p1out,true,p1in,true);
	mj1m=joi(p1out,false,p1in,false);
	mjH2p=jioH(p2in,true,p2out,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jioH(p2in,false,p2out,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(p1out,true,pg,true);
	j2gm=joo(p1out,false,pg,false);
	jgap=joi(pg,true,p1in,true);
	jgam=joi(pg,false,p1in,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	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();
	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();
	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();
	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();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}

	double jM2unogqbarHQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){
	// 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

	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=jio(p1in,true,p1out,true);
	mj1m=jio(p1in,false,p1out,false);
	mjH2p=jioH(p2in,true,p2out,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jioH(p2in,false,p2out,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(pg,true,p1out,true);
	j2gm=joo(pg,false,p1out,false);
	jgap=jio(p1in,true,pg,true);
	jgam=jio(p1in,false,pg,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	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();
	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();
	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();
	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();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}

	double jM2unogqHg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){
	// 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

	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=joi(p1out,true,p1in,true);
	mj1m=joi(p1out,false,p1in,false);
	mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(p1out,true,pg,true);
	j2gm=joo(p1out,false,pg,false);
	jgap=joi(pg,true,p1in,true);
	jgam=joi(pg,false,p1in,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	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();
	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();
	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();
	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();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4./9.4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.
	// here we need 2 to match with the normalization
	// gq is 9./4. times the qQ

	const double K = K_g(p2out, p2in);

	return ampsqK/HEJ::C_A9./4.; //ca/cf = 9/4
	}

	double jM2unogqbarHg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){
	// 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

	CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p;
	mj1p=jio(p1in,true,p1out,true);
	mj1m=jio(p1in,false,p1out,false);
	mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
	mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);

	// Dot products of these which occur again and again
	COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
	COM MHmm=mj1m.dot(mjH2m);
	COM MHpp=mj1p.dot(mjH2p);
	COM MHpm=mj1p.dot(mjH2m);

	// Currents with pg
	CCurrent jgam,jgap,j2gm,j2gp;
	j2gp=joo(pg,true,p1out,true);
	j2gm=joo(pg,false,p1out,false);
	jgap=jio(p1in,true,pg,true);
	jgam=jio(p1in,false,pg,false);

	CCurrent qsum(q1+qg);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
	CCurrent p1o(p1out);
	CCurrent p1i(p1in);

	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();
	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();
	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();
	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();

	U1mm=(jgam.dot(mjH2m)j2gm+2.p1o*MHmm)/(p1out+pg).m2();
	U1mp=(jgam.dot(mjH2p)j2gm+2.p1o*MHmp)/(p1out+pg).m2();
	U1pm=(jgap.dot(mjH2m)j2gp+2.p1o*MHpm)/(p1out+pg).m2();
	U1pp=(jgap.dot(mjH2p)j2gp+2.p1o*MHpp)/(p1out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH2m)jgam+2.p1iMHmm)/(p1in-pg).m2();
	U2mp=((-1.)j2gm.dot(mjH2p)jgam+2.p1iMHmp)/(p1in-pg).m2();
	U2pm=((-1.)j2gp.dot(mjH2m)jgap+2.p1iMHpm)/(p1in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH2p)jgap+2.p1iMHpp)/(p1in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());

	// Now add the t-channels for the Higgs
	double th=qH1.m2()*qg.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4./9.4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.
	// here we need 2 to match with the normalization
	// gq is 9./4. times the qQ

	const double K = K_g(p2out, p2in);

	return ampsq*K/HEJ::C_F;
	}

	double jM2unobqHQg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){

	HLV q1=p1in-p1out; // Top End
	HLV q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	HLV q3=-(p2in-p2out); // Bottom End

	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2, mt, incBot, mb);
	mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2, mt, incBot, mb);
	mj2p=joi(p2out,true,p2in,true);
	mj2m=joi(p2out,false,p2in,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);

	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(p2out,true,pg,true);
	j2gm=joo(p2out,false,pg,false);
	jgbp=joi(pg,true,p2in,true);
	jgbm=joi(pg,false,p2in,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	Lmm=((-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m+2.mj2m.dot(pg)mjH1m
	+ (p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmm/2.))/q3.m2();
	Lmp=((-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p+2.mj2p.dot(pg)mjH1m
	+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHmp/2.))/q3.m2();
	Lpm=((-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m+2.mj2m.dot(pg)mjH1p
	+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpm/2.))/q3.m2();
	Lpp=((-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p+2.mj2p.dot(pg)mjH1p
	+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))(q2.m2()MHpp/2.))/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	// 1/3. = 1/HEJ::C_A ?
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// Now add the t-channels for the Higgs
	const double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=HEJ::C_FHEJ::C_F/(HEJ::C_A*HEJ::C_A); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}

	double jM2unobqbarHQg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){

	HLV q1=p1in-p1out; // Top End
	HLV q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	HLV q3=-(p2in-p2out); // Bottom End

	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jioHtop(p1in,true,p1out,true,qH1,qH2, mt, incBot, mb);
	mjH1m=jioHtop(p1in,false,p1out,false,qH1,qH2, mt, incBot, mb);
	mj2p=joi(p2out,true,p2in,true);
	mj2m=joi(p2out,false,p2in,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);

	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(p2out,true,pg,true);
	j2gm=joo(p2out,false,pg,false);
	jgbp=joi(pg,true,p2in,true);
	jgbm=joi(pg,false,p2in,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);

	Lmm=( (-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmm/2. ) )/q3.m2();
	Lmp=( (-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmp/2. ) )/q3.m2();
	Lpm=( (-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpm/2. ) )/q3.m2();
	Lpp=( (-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpp/2. ) )/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}

	double jM2unobqHQbarg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){

	HLV q1=p1in-p1out; // Top End
	HLV q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	HLV q3=-(p2in-p2out); // Bottom End

	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb);
	mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb);
	mj2p=jio(p2in,true,p2out,true);
	mj2m=jio(p2in,false,p2out,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);

	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(pg,true,p2out,true);
	j2gm=joo(pg,false,p2out,false);
	jgbp=jio(p2in,true,pg,true);
	jgbm=jio(p2in,false,pg,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);

	Lmm=( (-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmm/2. ) )/q3.m2();
	Lmp=( (-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmp/2. ) )/q3.m2();
	Lpm=( (-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpm/2. ) )/q3.m2();
	Lpp=( (-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpp/2. ) )/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}

	double jM2unobqbarHQbarg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){

	HLV q1=p1in-p1out; // Top End
	HLV q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	HLV q3=-(p2in-p2out); // Bottom End

	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jioHtop(p1in,true,p1out,true,qH1,qH2,mt, incBot, mb);
	mjH1m=jioHtop(p1in,false,p1out,false,qH1,qH2,mt, incBot, mb);
	mj2p=jio(p2in,true,p2out,true);
	mj2m=jio(p2in,false,p2out,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);

	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(pg,true,p2out,true);
	j2gm=joo(pg,false,p2out,false);
	jgbp=jio(p2in,true,pg,true);
	jgbm=jio(p2in,false,pg,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);

	Lmm=( (-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmm/2. ) )/q3.m2();
	Lmp=( (-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmp/2. ) )/q3.m2();
	Lpm=( (-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpm/2. ) )/q3.m2();
	Lpp=( (-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpp/2. ) )/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4.4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.

	return ampsq;
	}

	double jM2unobgHQg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){

	HLV q1=p1in-p1out; // Top End
	HLV q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	HLV q3=-(p2in-p2out); // Bottom End

	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb);
	mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb);
	mj2p=joi(p2out,true,p2in,true);
	mj2m=joi(p2out,false,p2in,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);

	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(p2out,true,pg,true);
	j2gm=joo(p2out,false,pg,false);
	jgbp=joi(pg,true,p2in,true);
	jgbm=joi(pg,false,p2in,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	// COM test=pminus.dot(p1in);

	Lmm=( (-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmm/2. ) )/q3.m2();
	Lmp=( (-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmp/2. ) )/q3.m2();
	Lpm=( (-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpm/2. ) )/q3.m2();
	Lpp=( (-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpp/2. ) )/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4./9.4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.
	// need twice to match the normalization

	const double K = K_g(p1out, p1in);

	return ampsq*K/HEJ::C_F;
	}

	double jM2unobgHQbarg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
	HLV qH1, HLV qH2, double mt, bool incBot, double mb
	){

	HLV q1=p1in-p1out; // Top End
	HLV q2=-(p2in-p2out-pg); // Extra bit pre-gluon
	HLV q3=-(p2in-p2out); // Bottom End

	CCurrent mjH1m,mjH1p,mj2m,mj2p;
	mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb);
	mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb);
	mj2p=jio(p2in,true,p2out,true);
	mj2m=jio(p2in,false,p2out,false);

	// Dot products of these which occur again and again
	COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones
	COM MHmm=mjH1m.dot(mj2m);
	COM MHpp=mjH1p.dot(mj2p);
	COM MHpm=mjH1p.dot(mj2m);

	// Currents with pg
	CCurrent jgbm,jgbp,j2gm,j2gp;
	j2gp=joo(pg,true,p2out,true);
	j2gm=joo(pg,false,p2out,false);
	jgbp=jio(p2in,true,pg,true);
	jgbm=jio(p2in,false,pg,false);

	CCurrent qsum(q2+q3);

	CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
	CCurrent p2o(p2out);
	CCurrent p2i(p2in);

	CCurrent pplus((p1in+p1out)/2.);
	CCurrent pminus((p2in+p2out)/2.);

	Lmm=( (-1.)qsum(MHmm) + (-2.mjH1m.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmm/2. ) )/q3.m2();
	Lmp=( (-1.)qsum(MHmp) + (-2.mjH1m.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1m
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHmp/2. ) )/q3.m2();
	Lpm=( (-1.)qsum(MHpm) + (-2.mjH1p.dot(pg))mj2m + 2.mj2m.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpm/2. ) )/q3.m2();
	Lpp=( (-1.)qsum(MHpp) + (-2.mjH1p.dot(pg))mj2p + 2.mj2p.dot(pg)mjH1p
	+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )( q2.m2()MHpp/2. ) )/q3.m2();
	U1mm=(jgbm.dot(mjH1m)j2gm+2.p2o*MHmm)/(p2out+pg).m2();
	U1mp=(jgbp.dot(mjH1m)j2gp+2.p2o*MHmp)/(p2out+pg).m2();
	U1pm=(jgbm.dot(mjH1p)j2gm+2.p2o*MHpm)/(p2out+pg).m2();
	U1pp=(jgbp.dot(mjH1p)j2gp+2.p2o*MHpp)/(p2out+pg).m2();
	U2mm=((-1.)j2gm.dot(mjH1m)jgbm+2.p2iMHmm)/(p2in-pg).m2();
	U2mp=((-1.)j2gp.dot(mjH1m)jgbp+2.p2iMHmp)/(p2in-pg).m2();
	U2pm=((-1.)j2gm.dot(mjH1p)jgbm+2.p2iMHpm)/(p2in-pg).m2();
	U2pp=((-1.)j2gp.dot(mjH1p)jgbp+2.p2iMHpp)/(p2in-pg).m2();

	const double cf=HEJ::C_F;
	double amm,amp,apm,app;

	amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2());

	// Now add the t-channels for the Higgs
	double th=qH2.m2()*q2.m2();
	ampsq/=th;
	ampsq/=16.;
	ampsq=4./9.4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.

	const double K = K_g(p1out, p1in);

	return ampsq*K/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 = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor(-s12D0DD(k2, k1, q2, mq)(1 - 8.mqmq/s12 + S2/(2.s12) +
	S2(s12 - 8.mqmq)(s34 + S1)/(2.s12Delta) +
	2.(s34 + S1)(s34 + S1)/Delta +
	S2pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)C0DD(k2,
	k1 + q2, mq) -
	s12C0DD(k1, k2, mq) + (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S1*C0DD(k1, q2,
	mq))(S2(s12 - 4.mqmq)/(2.s12Delta) +
	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)/(DeltaSigma) +
	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))/(DeltaSigma));
	}

	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 = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor(-S2D0DD(k1, k2, q2,
	mq)(0.5 - (s12 - 8.mqmq)(s34 + S2)/(2.*Delta) -
	s12pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)C0DD(k1,
	k2 + q2, mq) -
	s12C0DD(k1, k2, mq) - (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S2*C0DD(k2, q2,
	mq))(S2(s12 - 4.mqmq)/(2.s12Delta) +
	S2pow((s34 + S2),2)/Delta/Delta) - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))(1. - 4.mqmq/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) +
	s12C0DD(k1 + k2, q2, mq))2.s34(s34 +
	S2)(S2 - S1)/(DeltaSigma) + (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))/(DeltaSigma));
	}

	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 = s12s34 - S1S2;

	return looprwfactor(S2D0DD(k1, q2, k2,
	mq)(Delta/s12/s12 - 4.mq*mq/s12) -
	S2((s12 + S1)C0DD(k1, k2 + q2, mq) -
	S1C0DD(k1, q2, mq))(1./
	s12/s12 - (s12 - 4.mqmq)/(2.s12Delta)) -
	S2((s12 + S2)C0DD(k1 + q2, k2, mq) -
	S2C0DD(k2, q2, mq))(1./
	s12/s12 + (s12 - 4.mqmq)/(2.s12Delta)) -
	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 = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor* (-s12*D0DD(k2, k1, q2,
	mq)(0.5 - (S1 - 8.mqmq)(s34 + S1)/(2.*Delta) -
	s12pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)C0DD(k2,
	k1 + q2, mq) -
	s12C0DD(k1, k2, mq) + (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S1C0DD(k1, q2, mq))((S1 - 4.mqmq)/(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.s12s34 - 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.s12s34 - S1*(S1 + S2) +
	s34(S2 - S1)))/(DeltaSigma)));
	}

	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 = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor* (-s12*D0DD(k1, k2, q2,
	mq)(0.5 + (S1 - 8.mqmq)(s34 + S2)/(2.*Delta) +
	s12pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)C0DD(k1,
	k2 + q2, mq) -
	s12C0DD(k1, k2, mq) - (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S2C0DD(k2, q2, mq))((S1 - 4.mqmq)/(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) +
	s12C0DD(k1 + k2, q2, mq))2.s34(2.s12s34 -
	S1*(S1 + S2) +
	s12(S2 - S1))/(DeltaSigma) - (B0DD(k1 + k2, mq) -
	B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)C0DD(k1 + k2, q2, mq))(2.s12(2.s12s34 - S2*(S1 + S2) +
	s34(S1 - S2))/(DeltaSigma)));
	}

	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 = s12s34 - S1S2;

	return looprwfactor* (-D0DD(k1, q2, k2,
	mq)*(Delta/s12 + (s12 + S1)/2. -
	4.mqmq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
	S1C0DD(k1, q2, mq))(1./
	s12 - (S1 - 4.mqmq)/(2.Delta)) + ((s12 + S2)C0DD(
	k1 + q2, k2, mq) -
	S2C0DD(k2, q2, mq))(1./
	s12 + (S1 - 4.mqmq)/(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 = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor(-s12D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta +
	12.mqmqS1(s34 + S1)/Delta/Delta -
	4.s12S1pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)C0DD(k2, k1 + q2, mq) -
	s12C0DD(k1, k2, mq) + (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S1C0DD(k1, q2, mq))(1./Delta +
	4.mqmq*S1/Delta/Delta -
	4.s12S1*pow((s34 + S1),2)/Delta/Delta/Delta) +
	C0DD(k1 + k2, q2, mq)(4.s12s34(S1 - S2)/(Delta*Sigma) -
	4.*(s12 -
	2.mqmq)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma)) + (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) +
	s12C0DD(k1 + k2, q2, mq))(12.s34(2.*s12 + S1 +
	S2)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma*Sigma) -
	4.s34(4.s12 + 3.S1 +
	S2)/(Delta*Sigma) +
	8.s12s34(s34(s12 + S2) -
	S1*(s34 +
	S1))/(DeltaDeltaSigma)) + (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))/(DeltaSigma*Sigma) +
	8.s12S1(s34(s12 + S2) -
	S1*(s34 +
	S1))/(DeltaDeltaSigma))) + (COM(0.,1.)/(4.M_PIM_PI))((2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma));
	}

	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 = s12s34 - S1S2;
	Sigma = 4.s12s34 - pow(S1+S2,2);

	return looprwfactor* (s12*D0DD(k1, k2, q2,
	mq)((s34 + S2)/Delta - 4.mq*mq/Delta +
	12.mqmqs34(s12 + S1)/Delta/Delta -
	4.s12pow((s34 + S2),2)/Delta/Delta -
	4.s12S1pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)C0DD(k1, k2 + q2, mq) -
	s12C0DD(k1, k2, mq) - (S1 - S2)C0DD(k1 + k2, q2, mq) -
	S2C0DD(k2, q2, mq))(1./Delta +
	4.mqmq*S1/Delta/Delta -
	4.s12(s34 + S2)/Delta/Delta -
	4.s12S1*pow((s34 + S2),2)/Delta/Delta/Delta) -
	C0DD(k1 + k2, q2, mq)(4.s12s34/(S2Delta) +
	4.s12s34(S2 - S1)/(DeltaSigma) +
	4.*(s12 -
	2.mqmq)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma)) - (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))(-12s34(2s12 + S1 +
	S2)(2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma*Sigma) -
	4.s12s34s34/(S2Delta*Delta) +
	4.s34S1/(Delta*Sigma) -
	4.s34(s12s34(2.*s12 + S2) -
	S1S1(2.*s12 +
	S1))/(DeltaDeltaSigma)) - (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))/(DeltaSigma*Sigma) +
	8.s12(2.s34 + S1)/(DeltaSigma) -
	8.s12s34(2.s12s34 - S1(S1 + S2) +
	s12*(S2 -
	S1))/(DeltaDeltaSigma))) + (COM(0.,1.)/(4.M_PIM_PI))((2.s12*s34 -
	S1(S1 + S2))/(DeltaSigma));

	}

	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 = s12s34 - S1S2;

	return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. +
	4.S1mqmq/Delta) + ((s12 + S1)C0DD(k1,
	k2 + q2, mq) -
	S1C0DD(k1, q2, mq))(1./Delta +
	4.S1mqmq/Delta/Delta) - ((s12 + S2)C0DD(k1 + q2,
	k2, mq) - S2C0DD(k2, q2, mq))(1./Delta +
	4.S1mq*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.mqmq + 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, current &retAns)
	{
	current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1,
	epsHapart2,conjepsH1part1,conjepsH1part2;
	COM ang1a,sqa1;

	const double F = 4.mqmq/HEJ::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.Ftacdot(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.)Ftapa.dot(p1)*aH1/sqa1, conjeps1, T5);
	cmult(-sqrt(2.)Ft1pa.dot(p1)*aH1/ang1a, epsa, T6);

	cmult(-aH1/sqrt(2.)/sqa1h48.COM(0.,1.)M_PI*M_PI, conjeps1, T7);
	cmult(aH1/sqrt(2.)/ang1ah58.COM(0.,1.)M_PI*M_PI, epsa, T8);
	cmult(aH1aH1/2./ang1a/sqa1h108.COM(0.,1.)M_PIM_PI, pacur, T9);
	cmult(-aH1aH1/2./ang1a/sqa1h128.COM(0.,1.)M_PIM_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, current &retAns)
	{
	const double F = 4.mqmq/HEJ::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.Ftacdot(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.phaseFta*pa.dot(pH), p1cur, T5a);
	cmult(2.phaseFt1*p1.dot(pH), pacur, T5b);
	cmult(-sqrt(2.)Ftap1.dot(pa)*oneHa/sqa1, conjeps1, T6);
	cmult(-sqrt(2.)Ft1pa.dot(p1)*aH1/sqa1, epsa, T7);

	cmult(-boxdiagFactphaseh2, pacur, T8a);
	cmult(boxdiagFactphaseh1, p1cur, T8b);
	cmult(boxdiagFactaH1/sqrt(2.)/sqa1h5, epsa, T9);
	cmult(-boxdiagFactoneHa/sqrt(2.)/sqa1h4, conjeps1, T10);
	cmult(boxdiagFactaH1oneHa/2./sqa1/sqa1*h10, pacur, T11a);
	cmult(-boxdiagFactaH1oneHa/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(-phaseFbeta(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 MH2gq_outsideH(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pH, double mq, bool includeBottom, double mq2)
	{

	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,retAns);
	app1=cdot(retAns,cur2bplus);
	app2=cdot(retAns,cur2bminus);

	g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
	app3=cdot(retAns,cur2bplusFlip);
	app4=cdot(retAns,cur2bminusFlip);

	// And non-conserving bits
	g_gH_HNC(p1in,p1out,pH,mq,retAns);
	apm1=cdot(retAns,cur2bplus);
	apm2=cdot(retAns,cur2bminus);

	g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
	apm3=cdot(retAns,cur2bplusFlip);
	apm4=cdot(retAns,cur2bminusFlip);
	} else {
	g_gH_HC(p1in,p1out,pH,mq,retAns);
	g_gH_HC(p1in,p1out,pH,mq2,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,retAns);
	g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,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,retAns);
	g_gH_HNC(p1in,p1out,pH,mq2,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,retAns);
	g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,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)
	{
	static double A=1./(3.M_PIHEJ::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 og 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/p1pconj(p1perp)p3perp+p1p/p2pp1perp*conj(p3perp));
	temp=temp*conj(temp);
	return temp.real();
	}

	double C2gHgp(HLV p2, HLV p1, HLV pH)
	{
	static double A=1./(3.M_PIHEJ::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(p1perpp3perp)*pow(php/p1p,2)/(1.+php/p1p)
	+s12(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1pphm)
	-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)
	{
	static double A=1./(3.M_PIHEJ::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)p3perpconj(p1perp)
	+sqrt(p1p/p2p)p1perpconj(p3perp) );
	temp=temp*conj(temp);
	return temp.real();
	}

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 6:09 AM (1 d, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6542481
Default Alt Text: currents.cc (83 KB)

currents.ccNo OneTemporaryActions

currents.ccView Options

File Metadata

Event Timeline

currents.cc
No OneTemporary
Actions

currents.cc
View Options