diff --git a/include/HEJ/jets.hh b/include/HEJ/jets.hh
index 4ffac16..0013971 100644
--- a/include/HEJ/jets.hh
+++ b/include/HEJ/jets.hh
@@ -1,395 +1,339 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in pure jets.
*
* This file contains all the necessary functions to compute the
* current contractions for all valid pure jet HEJ processes. It will
* also contain some pure jet ME components used in other process ME
* calculations
*
* @TODO add a namespace
*/
#pragma once
#include <complex>
#include <ostream>
#include <CLHEP/Vector/LorentzVector.h>
typedef std::complex<double> COM;
typedef COM current[4];
typedef CLHEP::HepLorentzVector HLV;
//! Square of qQ->qQ Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qQ Scattering
*/
double ME_qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qQbar Scattering
*
* @note this can be used for qbarQ->qbarQ Scattering by inputting arguments
* appropriately.
*/
double ME_qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*/
double ME_qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qg->qg Scattering
*
* @note this can be used for gq->gq Scattering by inputting arguments
* appropriately.
*/
double ME_qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* @note this can be used for gqbar->gqbar Scattering by inputting arguments
* appropriately.
*/
double ME_qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of gg->gg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for gg->gg Scattering
*/
double ME_gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//Unordered Backwards contributions:
//! Square of qQ->qQ Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qQ Scattering
*/
double ME_unob_qQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQ->qbarQ Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qbarQ->qbarQ Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
*/
double ME_unob_qbarQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qQbar Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qQbar Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
*/
double ME_unob_qQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarQbar Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
*/
double ME_unob_qbarQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qg Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qg->qg Scattering
+ *
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
*/
double ME_unob_qg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarg Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarg->qbarg Scattering
- */
-double ME_unob_qbarg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
-
-//Unordered Forwards contributions:
-//! Square of qQ->qQ Pure Jets Scattering Current
-/**
- * @param p1out Momentum of final state quark
- * @param p1in Momentum of initial state quark
- * @param pg Momentum of unordered gluon
- * @param p2out Momentum of final state quark
- * @param p2in Momentum of intial state quark
- * @returns Square of the current contractions for qQ->qQ Scattering
- */
-double ME_unof_qQ (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in);
-
-//! Square of qbarQ->qbarQ Pure Jets Unordered forwards Scattering Current
-/**
- * @param p1out Momentum of final state anti-quark
- * @param p1in Momentum of initial state anti-quark
- * @param pg Momentum of unordered gluon
- * @param p2out Momentum of final state quark
- * @param p2in Momentum of intial state quark
- * @returns Square of the current contractions for qbarQ->qbarQ Scattering
- */
-double ME_unof_qbarQ (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in);
-
-//! Square of qQbar->qQbar Pure Jets Unordered forwards Scattering Current
-/**
- * @param p1out Momentum of final state quark
- * @param p1in Momentum of initial state quark
- * @param pg Momentum of unordered gluon
- * @param p2out Momentum of final state anti-quark
- * @param p2in Momentum of intial state anti-quark
- * @returns Square of the current contractions for qQbar->qQbar Scattering
*
- * @note this can be used for qbarQ->qbarQ Scattering by inputting arguments
- * appropriately.
- */
-double ME_unof_qQbar (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in);
-
-//! Square of qbarQbar->qbarQbar Pure Jets Unordered forwards Scattering Current
-/**
- * @param p1out Momentum of final state anti-quark
- * @param p1in Momentum of initial state anti-quark
- * @param pg Momentum of unordered gluon
- * @param p2out Momentum of final state anti-quark
- * @param p2in Momentum of intial state anti-quark
- * @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
+ * @note this can be used for unof contributions by inputting
+ * arguments appropriately.
*/
-double ME_unof_qbarQbar (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in);
-
-//! Square of gQ->gQ Pure Jets Unordered forwards Scattering Current
-/**
- * @param p1out Momentum of final state gluon
- * @param p1in Momentum of initial state gluon
- * @param pg Momentum of unordered gluon
- * @param p2out Momentum of final state quark
- * @param p2in Momentum of intial state quark
- * @returns Square of the current contractions for gQ->gQ Scattering
- */
-double ME_unof_gQ (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in);
-
-//! Square of gQbar->gQbar Pure Jets Unordered forwards Scattering Current
-/**
- * @param p1out Momentum of final state gluon
- * @param p1in Momentum of initial state gluon
- * @param pg Momentum of unordered gluon
- * @param p2out Momentum of final state anti-quark
- * @param p2in Momentum of intial state anti-quark
- * @returns Square of the current contractions for gQbar->gQbar Scattering
- */
-double ME_unof_gQbar(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in);
-
+double ME_unob_qbarg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
/** \class CCurrent jets.hh "include/HEJ/jets.hh"
* \brief This is the a new class structure for currents.
*/
class CCurrent
{
public:
CCurrent(COM sc0, COM sc1, COM sc2, COM sc3)
:c0(sc0),c1(sc1),c2(sc2),c3(sc3)
{};
CCurrent(const HLV p)
{
c0=p.e();
c1=p.px();
c2=p.py();
c3=p.pz();
};
CCurrent()
{};
CCurrent operator+(const CCurrent& other);
CCurrent operator-(const CCurrent& other);
CCurrent operator*(const double x);
CCurrent operator*(const COM x);
CCurrent operator/(const double x);
CCurrent operator/(const COM x);
friend std::ostream& operator<<(std::ostream& os, const CCurrent& cur);
COM dot(HLV p1);
COM dot(CCurrent p1);
COM c0,c1,c2,c3;
};
/* std::ostream& operator <<(std::ostream& os, const CCurrent& cur); */
CCurrent operator * ( double x, CCurrent& m);
CCurrent operator * ( COM x, CCurrent& m);
CCurrent operator / ( double x, CCurrent& m);
CCurrent operator / ( COM x, CCurrent& m);
//! Current <incoming state | mu | outgoing state>
/**
* This is a wrapper function around \see joi() note helicity flip to
* give same answer.
*/
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur);
//! Current <outgoing state | mu | outgoing state>
/**
* @param pi bra state momentum
* @param heli helicity of pi
* @param pj ket state momentum
* @param helj helicity of pj. (must be same as heli)
* @param cur reference to current which is saved.
*
* This function is for building <i (out)| mu |j (out)> currents. It
* must be called with pi as the bra, and pj as the ket.
*
* @TODO Remove heli/helj and just have helicity of current as argument.
*/
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* @param pout bra state momentum
* @param helout helicity of pout
* @param pin ket state momentum
* @param helin helicity of pin. (must be same as helout)
* @param cur reference to current which is saved.
*
* This function is for building <out| mu |in> currents. It must be
* called with pout as the bra, and pin as the ket. jio calls this
* with flipped helicity
*
* @TODO Remove helout/helin and just have helicity of current as argument.
*/
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* This is a wrapper function around the void function of the same name. \see joi
*/
CCurrent joi (HLV pout, bool helout, HLV pin, bool helin);
//! Current <incoming state | mu | outgoing state>
/**
* This is a wrapper function around the void function of the same name. \see jio
*/
CCurrent jio (HLV pout, bool helout, HLV pin, bool helin);
//! Current <outgoing state | mu | outgoing state>
/**
* This is a wrapper function around the void function of the same name. \see joo
*/
CCurrent joo (HLV pout, bool helout, HLV pin, bool helin);
inline COM cdot(const current & j1, const current & j2)
{
return j1[0]*j2[0]-j1[1]*j2[1]-j1[2]*j2[2]-j1[3]*j2[3];
}
inline COM cdot(const HLV & p, const current & j1) {
return j1[0]*p.e()-j1[1]*p.x()-j1[2]*p.y()-j1[3]*p.z();
}
inline void cmult(const COM & factor, const current & j1, current &cur)
{
cur[0]=factor*j1[0];
cur[1]=factor*j1[1];
cur[2]=factor*j1[2];
cur[3]=factor*j1[3];
}
// WHY!?!
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, const current & j5, current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0]+j4[0]+j5[0];
sum[1]=j1[1]+j2[1]+j3[1]+j4[1]+j5[1];
sum[2]=j1[2]+j2[2]+j3[2]+j4[2]+j5[2];
sum[3]=j1[3]+j2[3]+j3[3]+j4[3]+j5[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, current &sum) {
sum[0] = j1[0] + j2[0] + j3[0] + j4[0];
sum[1] = j1[1] + j2[1] + j3[1] + j4[1];
sum[2] = j1[2] + j2[2] + j3[2] + j4[2];
sum[3] = j1[3] + j2[3] + j3[3] + j4[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0];
sum[1]=j1[1]+j2[1]+j3[1];
sum[2]=j1[2]+j2[2]+j3[2];
sum[3]=j1[3]+j2[3]+j3[3];
}
inline void cadd(const current & j1, const current & j2, current &sum)
{
sum[0]=j1[0]+j2[0];
sum[1]=j1[1]+j2[1];
sum[2]=j1[2]+j2[2];
sum[3]=j1[3]+j2[3];
}
inline double abs2(const COM & a)
{
return (a*conj(a)).real();
}
inline double vabs2(const CCurrent & cur)
{
return abs2(cur.c0)-abs2(cur.c1)-abs2(cur.c2)-abs2(cur.c3);
}
inline double vre(const CCurrent & a, const CCurrent & b)
{
return real(a.c0*conj(b.c0)-a.c1*conj(b.c1)-a.c2*conj(b.c2)-a.c3*conj(b.c3));
}
//! @TODO These are not currents and should be moved elsewhere.
double K_g(double p1minus, double paminus);
double K_g(HLV const & pout, HLV const & pin);
diff --git a/src/jets.cc b/src/jets.cc
index 028aa7d..c9e3d95 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,1245 +1,737 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Constants.hh"
// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
// @TODO: this is not a current and should be moved somewhere else
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(HEJ::C_A - 1./HEJ::C_A) + 1./HEJ::C_A;
}
double K_g(
HLV const & pout,
HLV const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
CCurrent CCurrent::operator+(const CCurrent& other)
{
COM result_c0=c0 + other.c0;
COM result_c1=c1 + other.c1;
COM result_c2=c2 + other.c2;
COM result_c3=c3 + other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator-(const CCurrent& other)
{
COM result_c0=c0 - other.c0;
COM result_c1=c1 - other.c1;
COM result_c2=c2 - other.c2;
COM result_c3=c3 - other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const double x)
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const double x)
{
COM result_c0=CCurrent::c0/x;
COM result_c1=CCurrent::c1/x;
COM result_c2=CCurrent::c2/x;
COM result_c3=CCurrent::c3/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const COM x)
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const COM x)
{
COM result_c0=(CCurrent::c0)/x;
COM result_c1=(CCurrent::c1)/x;
COM result_c2=(CCurrent::c2)/x;
COM result_c3=(CCurrent::c3)/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
{
os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
return os;
}
CCurrent operator * ( double x, CCurrent& m)
{
return m*x;
}
CCurrent operator * ( COM x, CCurrent& m)
{
return m*x;
}
CCurrent operator / ( double x, CCurrent& m)
{
return m/x;
}
CCurrent operator / ( COM x, CCurrent& m)
{
return m/x;
}
COM CCurrent::dot(HLV p1)
{
// Current goes (E,px,py,pz)
// Vector goes (px,py,pz,E)
return p1[3]*c0-p1[0]*c1-p1[1]*c2-p1[2]*c3;
}
COM CCurrent::dot(CCurrent p1)
{
return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3;
}
//Current Functions
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
const double sqpop = sqrt(pout.plus());
const double sqpom = sqrt(pout.minus());
const COM poperp = pout.x() + COM(0, 1) * pout.y();
if (helout != helin) {
throw std::invalid_argument{"Non-matching helicities"};
} else if (helout == false) { // negative helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(pin.plus());
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * poperp / abs(poperp);
cur[2] = -COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(pin.minus());
cur[0] = -sqpom * sqpim * poperp / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = COM(0,1) * cur[1];
cur[3] = -cur[0];
}
} else { // positive helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(pin.plus());
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
cur[2] = COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(pin.minus());
cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = -COM(0,1) * cur[1];
cur[3] = -cur[0];
}
}
}
CCurrent joi (HLV pout, bool helout, HLV pin, bool helin)
{
current cur;
joi(pout, helout, pin, helin, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) {
joi(pout, !helout, pin, !helin, cur);
}
CCurrent jio (HLV pin, bool helin, HLV pout, bool helout)
{
current cur;
jio(pin, helin, pout, helout, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {
// Zero our current
cur[0] = 0.0;
cur[1] = 0.0;
cur[2] = 0.0;
cur[3] = 0.0;
if (heli!=helj) {
throw std::invalid_argument{"Non-matching helicities"};
} else if ( heli == true ) { // If positive helicity swap momenta
std::swap(pi,pj);
}
const double sqpjp = sqrt(pj.plus());
const double sqpjm = sqrt(pj.minus());
const double sqpip = sqrt(pi.plus());
const double sqpim = sqrt(pi.minus());
const COM piperp = pi.x() + COM(0,1) * pi.y();
const COM pjperp = pj.x() + COM(0,1) * pj.y();
const COM phasei = piperp / abs(piperp);
const COM phasej = pjperp / abs(pjperp);
cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
}
CCurrent joo (HLV pi, bool heli, HLV pj, bool helj)
{
current cur;
joo(pi, heli, pj, helj, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
namespace{
//@{
/**
* @brief Pure Jet FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1.
* @param p1in Incoming Particle 1.
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_\mu j^\mu.
* Handles all possible incoming states.
*/
double j_j(HLV p1out, HLV p1in, HLV p2out, HLV p2in, bool aqlineb, bool aqlinef){
HLV q1=p1in-p1out;
HLV q2=-(p2in-p2out);
current mj1m,mj1p,mj2m,mj2p;
// Note need to flip helicities in anti-quark case.
joi(p1out,!aqlineb, p1in,!aqlineb, mj1p);
joi(p1out, aqlineb, p1in, aqlineb, mj1m);
joi(p2out,!aqlinef, p2in,!aqlinef, mj2p);
joi(p2out, aqlinef, p2in, aqlinef, mj2m);
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_F*HEJ::C_F*(sst)/(q1.m2()*q2.m2());
}
} //anonymous namespace
double ME_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false);
}
double ME_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, true);
}
double ME_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, true, true);
}
double ME_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false)*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, true, false)*K_g(p2out, p2in)/(HEJ::C_F);
}
double ME_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false)*K_g(p1out, p1in)*K_g(p2out, p2in)/(HEJ::C_F*HEJ::C_F);
}
//@}
//Unordered bits for pure jet
double ME_unob_qQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
// This construction is taking rapidity order: pg > p1out >> p2out
// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;
CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon
// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
CCurrent mj1m,mj1p,mj2m,mj2p;
mj1p=joi(p1out,true,p1in,true);
mj1m=joi(p1out,false,p1in,false);
mj2p=joi(p2out,true,p2in,true);
mj2m=joi(p2out,false,p2in,false);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
// std::cout<< p1out.rapidity() << " " << p2out.rapidity()<< " " << qH1 << " " << qH2 << "\n" <<MHmm << " " << MHmp << " " << MHpm << " " << MHpp << std::endl;
// 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*(Mmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
double cf=4./3.;
double amm,amp,apm,app;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// Now add the t-channels for the Higgs
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
//Higgs coupling is included in Hjets.C
return ampsq;
}
double ME_unob_qbarQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
// This construction is taking rapidity order: pg > p1out >> p2out
// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;
CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon
// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
CCurrent mj1m,mj1p,mj2m,mj2p;
mj1p=jio(p1in,true,p1out,true);
mj1m=jio(p1in,false,p1out,false);
mj2p=joi(p2out,true,p2in,true);
mj2m=joi(p2out,false,p2in,false);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
// 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*(Mmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
double cf=4./3.;
double amm,amp,apm,app;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// Now add the t-channels for the Higgs
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
//Higgs coupling is included in Hjets.C
return ampsq;
}
double ME_unob_qQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
// This construction is taking rapidity order: pg > p1out >> p2out
// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;
CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon
// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
CCurrent mj1m,mj1p,mj2m,mj2p;
mj1p=joi(p1out,true,p1in,true);
mj1m=joi(p1out,false,p1in,false);
mj2p=jio(p2in,true,p2out,true);
mj2m=jio(p2in,false,p2out,false);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
// 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*(Mmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
double cf=4./3.;
double amm,amp,apm,app;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// Now add the t-channels for the Higgs
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
//Higgs coupling is included in Hjets.C
return ampsq;
}
double ME_unob_qbarQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
// This construction is taking rapidity order: pg > p1out >> p2out
// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;
CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon
// std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
CCurrent mj1m,mj1p,mj2m,mj2p;
mj1p=jio(p1in,true,p1out,true);
mj1m=jio(p1in,false,p1out,false);
mj2p=jio(p2in,true,p2out,true);
mj2m=jio(p2in,false,p2out,false);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
// 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*(Mmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
double cf=4./3.;
double amm,amp,apm,app;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// Now add the t-channels for the Higgs
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
//Higgs coupling is included in Hjets.C
ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
return ampsq;
}
double ME_unob_qg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
// This construction is taking rapidity order: pg > p1out >> p2out
// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;
CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon
CCurrent mj1m,mj1p,mj2m,mj2p;
mj1p=joi(p1out,true,p1in,true);
mj1m=joi(p1out,false,p1in,false);
mj2p=joi(p2out,true,p2in,true);
mj2m=joi(p2out,false,p2in,false);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
// 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*(Mmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
double cf=4./3.;
double amm,amp,apm,app;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// Now add the t-channels for the Higgs
double th=q2.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
//Higgs coupling is included in Hjets.C
double ratio; // p2-/pb- in the notes
// if (p2in.plus()>0) // if the gluon is the positive
if (p2in.pz()>0) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.;
return ampsq*nonflipcolourmult*9./4.; //ca/cf = 9/4
}
double ME_unob_qbarg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
// This construction is taking rapidity order: pg > p1out >> p2out
// std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl;
CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End
CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon
CCurrent mj1m,mj1p,mj2m,mj2p;
mj1p=jio(p1in,true,p1out,true);
mj1m=jio(p1in,false,p1out,false);
mj2p=joi(p2out,true,p2in,true);
mj2m=joi(p2out,false,p2in,false);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p); // And now for the iggs ones
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
// 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*(Mmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
double cf=4./3.;
double amm,amp,apm,app;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
// std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
// Now add the t-channels for the Higgs
double th=q2.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
//Higgs coupling is included in Hjets.C
double ratio; // p2-/pb- in the notes
// if (p2in.plus()>0) // if the gluon is the positive
if (p2in.pz()>0) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.;
return ampsq*nonflipcolourmult*9./4.; //ca/cf = 9/4
}
-
-double ME_unof_qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
-{
- // std::cout << "####################\n";
- // std::cout << "# p1in : "<<p1in<< " "<<p1in.plus()<<" "<<p1in.minus()<<std::endl;
- // std::cout << "# p2in : "<<p2in<< " "<<p2in.plus()<<" "<<p2in.minus()<<std::endl;
- // std::cout << "# p1out : "<<p1out<< " "<<p1out.rapidity()<<std::endl;
- // std::cout << "# (qH1-qH2) : "<<(qH1-qH2)<< " "<<(qH1-qH2).rapidity()<<std::endl;
- // std::cout << "# pg : "<<pg<< " "<<pg.rapidity()<<std::endl;
- // std::cout << "# p2out : "<<p2out<< " "<<p2out.rapidity()<<std::endl;
- // std::cout << "####################\n";
-
- CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
- CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
- CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End
-
- // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
- CCurrent mj1m,mj1p,mj2m,mj2p;
- mj1p=joi(p1out,true,p1in,true);
- mj1m=joi(p1out,false,p1in,false);
- mj2p=joi(p2out,true,p2in,true);
- mj2m=joi(p2out,false,p2in,false);
-
- // Dot products of these which occur again and again
- COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
- COM Mmm=mj1m.dot(mj2m);
- COM Mpp=mj1p.dot(mj2p);
- COM Mpm=mj1p.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*(Mmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmm/2.))/q3.m2();
- Lmp=((-1.)*qsum*(Mmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmp/2.))/q3.m2();
- Lpm=((-1.)*qsum*(Mpm) + (-2.*mj1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpm/2.))/q3.m2();
- Lpp=((-1.)*qsum*(Mpp) + (-2.*mj1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpp/2.))/q3.m2();
- U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*Mmm)/(p2out+pg).m2();
- U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*Mmp)/(p2out+pg).m2();
- U1pm=(jgbm.dot(mj1p)*j2gm+2.*p2o*Mpm)/(p2out+pg).m2();
- U1pp=(jgbp.dot(mj1p)*j2gp+2.*p2o*Mpp)/(p2out+pg).m2();
- U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*Mmm)/(p2in-pg).m2();
- U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*Mmp)/(p2in-pg).m2();
- U2pm=((-1.)*j2gm.dot(mj1p)*jgbm+2.*p2i*Mpm)/(p2in-pg).m2();
- U2pp=((-1.)*j2gp.dot(mj1p)*jgbp+2.*p2i*Mpp)/(p2in-pg).m2();
-
- double cf=4./3.;
- double amm,amp,apm,app;
-
- amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
- amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
- apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
- app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
- double ampsq=-(amm+amp+apm+app);
-
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
-
- // Now add the t-channels for the Higgs
- double th=q1.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 ME_unof_qbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
-{
-
- CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
- CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
- CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End
-
- // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
- CCurrent mj1m,mj1p,mj2m,mj2p;
- mj1p=jio(p1in,true,p1out,true);
- mj1m=jio(p1in,false,p1out,false);
- mj2p=joi(p2out,true,p2in,true);
- mj2m=joi(p2out,false,p2in,false);
-
- // Dot products of these which occur again and again
- COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
- COM Mmm=mj1m.dot(mj2m);
- COM Mpp=mj1p.dot(mj2p);
- COM Mpm=mj1p.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*(Mmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmm/2.))/q3.m2();
- Lmp=((-1.)*qsum*(Mmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmp/2.))/q3.m2();
- Lpm=((-1.)*qsum*(Mpm) + (-2.*mj1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpm/2.))/q3.m2();
- Lpp=((-1.)*qsum*(Mpp) + (-2.*mj1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpp/2.))/q3.m2();
- U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*Mmm)/(p2out+pg).m2();
- U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*Mmp)/(p2out+pg).m2();
- U1pm=(jgbm.dot(mj1p)*j2gm+2.*p2o*Mpm)/(p2out+pg).m2();
- U1pp=(jgbp.dot(mj1p)*j2gp+2.*p2o*Mpp)/(p2out+pg).m2();
- U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*Mmm)/(p2in-pg).m2();
- U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*Mmp)/(p2in-pg).m2();
- U2pm=((-1.)*j2gm.dot(mj1p)*jgbm+2.*p2i*Mpm)/(p2in-pg).m2();
- U2pp=((-1.)*j2gp.dot(mj1p)*jgbp+2.*p2i*Mpp)/(p2in-pg).m2();
-
-
- double cf=4./3.;
- double amm,amp,apm,app;
-
- amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
- amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
- apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
- app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
- double ampsq=-(amm+amp+apm+app);
-
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
-
- // Now add the t-channels for the Higgs
- double th=q1.m2()*q2.m2();
- ampsq/=th;
- ampsq/=16.;
- ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
- //Higgs coupling is included in Hjets.C
-
- return ampsq;
-}
-
-
-
-double ME_unof_qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
-{
-
- CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
- CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
- CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End
-
- // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
- CCurrent mj1m,mj1p,mj2m,mj2p;
- mj1p=joi(p1out,true,p1in,true);
- mj1m=joi(p1out,false,p1in,false);
- mj2p=jio(p2in,true,p2out,true);
- mj2m=jio(p2in,false,p2out,false);
-
- // Dot products of these which occur again and again
- COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
- COM Mmm=mj1m.dot(mj2m);
- COM Mpp=mj1p.dot(mj2p);
- COM Mpm=mj1p.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*(Mmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmm/2.))/q3.m2();
- Lmp=((-1.)*qsum*(Mmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmp/2.))/q3.m2();
- Lpm=((-1.)*qsum*(Mpm) + (-2.*mj1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpm/2.))/q3.m2();
- Lpp=((-1.)*qsum*(Mpp) + (-2.*mj1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpp/2.))/q3.m2();
- U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*Mmm)/(p2out+pg).m2();
- U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*Mmp)/(p2out+pg).m2();
- U1pm=(jgbm.dot(mj1p)*j2gm+2.*p2o*Mpm)/(p2out+pg).m2();
- U1pp=(jgbp.dot(mj1p)*j2gp+2.*p2o*Mpp)/(p2out+pg).m2();
- U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*Mmm)/(p2in-pg).m2();
- U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*Mmp)/(p2in-pg).m2();
- U2pm=((-1.)*j2gm.dot(mj1p)*jgbm+2.*p2i*Mpm)/(p2in-pg).m2();
- U2pp=((-1.)*j2gp.dot(mj1p)*jgbp+2.*p2i*Mpp)/(p2in-pg).m2();
-
-
-
- double cf=4./3.;
- double amm,amp,apm,app;
-
- amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
- amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
- apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
- app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
- double ampsq=-(amm+amp+apm+app);
-
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
-
- // Now add the t-channels for the Higgs
- double th=q1.m2()*q2.m2();
- ampsq/=th;
- ampsq/=16.;
- ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
- //Higgs coupling is included in Hjets.C
-
- return ampsq;
-}
-
-
-
-double ME_unof_qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
-{
-
- CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
- CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
- CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End
-
- // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
- CCurrent mj1m,mj1p,mj2m,mj2p;
- mj1p=jio(p1in,true,p1out,true);
- mj1m=jio(p1in,false,p1out,false);
- mj2p=jio(p2in,true,p2out,true);
- mj2m=jio(p2in,false,p2out,false);
-
- // Dot products of these which occur again and again
- COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
- COM Mmm=mj1m.dot(mj2m);
- COM Mpp=mj1p.dot(mj2p);
- COM Mpm=mj1p.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*(Mmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmm/2.))/q3.m2();
- Lmp=((-1.)*qsum*(Mmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmp/2.))/q3.m2();
- Lpm=((-1.)*qsum*(Mpm) + (-2.*mj1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpm/2.))/q3.m2();
- Lpp=((-1.)*qsum*(Mpp) + (-2.*mj1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpp/2.))/q3.m2();
- U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*Mmm)/(p2out+pg).m2();
- U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*Mmp)/(p2out+pg).m2();
- U1pm=(jgbm.dot(mj1p)*j2gm+2.*p2o*Mpm)/(p2out+pg).m2();
- U1pp=(jgbp.dot(mj1p)*j2gp+2.*p2o*Mpp)/(p2out+pg).m2();
- U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*Mmm)/(p2in-pg).m2();
- U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*Mmp)/(p2in-pg).m2();
- U2pm=((-1.)*j2gm.dot(mj1p)*jgbm+2.*p2i*Mpm)/(p2in-pg).m2();
- U2pp=((-1.)*j2gp.dot(mj1p)*jgbp+2.*p2i*Mpp)/(p2in-pg).m2();
-
-
-
- double cf=4./3.;
- double amm,amp,apm,app;
-
- amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
- amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
- apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
- app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
- double ampsq=-(amm+amp+apm+app);
-
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
-
- // Now add the t-channels for the Higgs
- double th=q1.m2()*q2.m2();
- ampsq/=th;
- ampsq/=16.;
- ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ.
- //Higgs coupling is included in Hjets.C
-
-
-
- return ampsq;
-}
-
-double ME_unof_gQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
-{
- // std::cout << "####################\n";
- // std::cout << "# p1in : "<<p1in<< " "<<p1in.plus()<<" "<<p1in.minus()<<std::endl;
- // std::cout << "# p2in : "<<p2in<< " "<<p2in.plus()<<" "<<p2in.minus()<<std::endl;
- // std::cout << "# p1out : "<<p1out<< " "<<p1out.rapidity()<<std::endl;
- // std::cout << "# (qH1-qH2) : "<<(qH1-qH2)<< " "<<(qH1-qH2).rapidity()<<std::endl;
- // std::cout << "# pg : "<<pg<< " "<<pg.rapidity()<<std::endl;
- // std::cout << "# p2out : "<<p2out<< " "<<p2out.rapidity()<<std::endl;
- // std::cout << "####################\n";
-
- CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
- CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
- CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End
-
- // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
- CCurrent mj1m,mj1p,mj2m,mj2p;
- mj1p=joi(p1out,true,p1in,true);
- mj1m=joi(p1out,false,p1in,false);
- mj2p=joi(p2out,true,p2in,true);
- mj2m=joi(p2out,false,p2in,false);
-
- // Dot products of these which occur again and again
- COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
- COM Mmm=mj1m.dot(mj2m);
- COM Mpp=mj1p.dot(mj2p);
- COM Mpm=mj1p.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*(Mmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmm/2.))/q3.m2();
- Lmp=((-1.)*qsum*(Mmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmp/2.))/q3.m2();
- Lpm=((-1.)*qsum*(Mpm) + (-2.*mj1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpm/2.))/q3.m2();
- Lpp=((-1.)*qsum*(Mpp) + (-2.*mj1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpp/2.))/q3.m2();
- U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*Mmm)/(p2out+pg).m2();
- U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*Mmp)/(p2out+pg).m2();
- U1pm=(jgbm.dot(mj1p)*j2gm+2.*p2o*Mpm)/(p2out+pg).m2();
- U1pp=(jgbp.dot(mj1p)*j2gp+2.*p2o*Mpp)/(p2out+pg).m2();
- U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*Mmm)/(p2in-pg).m2();
- U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*Mmp)/(p2in-pg).m2();
- U2pm=((-1.)*j2gm.dot(mj1p)*jgbm+2.*p2i*Mpm)/(p2in-pg).m2();
- U2pp=((-1.)*j2gp.dot(mj1p)*jgbp+2.*p2i*Mpp)/(p2in-pg).m2();
-
- double cf=4./3.;
- double amm,amp,apm,app;
-
- amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
- amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
- apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
- app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
- double ampsq=-(amm+amp+apm+app);
-
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
-
- // Now add the t-channels for the Higgs
- double th=q1.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
- //Higgs coupling is included in Hjets.C
-
- double ratio; // p2-/pb- in the notes
- // if (p2in.plus()>0) // if the gluon is the positive
- if (p1in.pz()>0) // if the gluon is the positive
- ratio=p1out.plus()/p1in.plus();
- else // the gluon is the negative
- ratio=p1out.minus()/p1in.minus();
-
- double nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.;
-
- return ampsq*nonflipcolourmult*9./4.; //ca/cf = 9/4
-}
-
-double ME_unof_gQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
-{
-
- CLHEP::HepLorentzVector q1=p1in-p1out; // Top End
- CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon
- CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End
-
- // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl;
- CCurrent mj1m,mj1p,mj2m,mj2p;
- mj1p=joi(p1out,true,p1in,true);
- mj1m=joi(p1out,false,p1in,false);
- mj2p=jio(p2in,true,p2out,true);
- mj2m=jio(p2in,false,p2out,false);
-
- // Dot products of these which occur again and again
- COM Mmp=mj1m.dot(mj2p); // And now for the Higgs ones
- COM Mmm=mj1m.dot(mj2m);
- COM Mpp=mj1p.dot(mj2p);
- COM Mpm=mj1p.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*(Mmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmm/2.))/q3.m2();
- Lmp=((-1.)*qsum*(Mmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mmp/2.))/q3.m2();
- Lpm=((-1.)*qsum*(Mpm) + (-2.*mj1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpm/2.))/q3.m2();
- Lpp=((-1.)*qsum*(Mpp) + (-2.*mj1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*Mpp/2.))/q3.m2();
- U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*Mmm)/(p2out+pg).m2();
- U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*Mmp)/(p2out+pg).m2();
- U1pm=(jgbm.dot(mj1p)*j2gm+2.*p2o*Mpm)/(p2out+pg).m2();
- U1pp=(jgbp.dot(mj1p)*j2gp+2.*p2o*Mpp)/(p2out+pg).m2();
- U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*Mmm)/(p2in-pg).m2();
- U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*Mmp)/(p2in-pg).m2();
- U2pm=((-1.)*j2gm.dot(mj1p)*jgbm+2.*p2i*Mpm)/(p2in-pg).m2();
- U2pp=((-1.)*j2gp.dot(mj1p)*jgbp+2.*p2i*Mpp)/(p2in-pg).m2();
-
-
-
- double cf=4./3.;
- double amm,amp,apm,app;
-
- amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
- amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
- apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
- app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
- double ampsq=-(amm+amp+apm+app);
-
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
- // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999)
- // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl;
-
- // Now add the t-channels for the Higgs
- double th=q1.m2()*q2.m2();
- ampsq/=th;
- ampsq/=16.;
- ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ.
- //Higgs coupling is included in Hjets.C
-
- double ratio; // p2-/pb- in the notes
- // if (p2in.plus()>0) // if the gluon is the positive
- if (p1in.pz()>0) // if the gluon is the positive
- ratio=p1out.plus()/p1in.plus();
- else // the gluon is the negative
- ratio=p1out.minus()/p1in.minus();
-
- double nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.;
-
- return ampsq*nonflipcolourmult*9./4.; //ca/cf = 9/4
-
-}