diff --git a/src/currents.cc b/src/currents.cc
index 779e1ca..7e917bf 100644
--- a/src/currents.cc
+++ b/src/currents.cc
@@ -1,3825 +1,3746 @@
 //////////////////////////////////////////////////
 //////////////////////////////////////////////////
 // This source code is Copyright (2012) of      //
 //  Jeppe R. Andersen and Jennifer M. Smillie   //
 // and is distributed under the                 //
 // Gnu Public License version 2                 //
 // http://www.gnu.org/licenses/gpl-2.0.html     //
 // You are allowed to distribute and alter the  //
 // source under the conditions of the GPLv2     //
 // as long as this copyright notice             //
 // is unaltered and distributed with the source //
 // Any use should comply with the               //
 //             MCNET GUIDELINES                 //
 //    for Event Generator Authors and Users     //
 // as distributed with this source code         //
 //////////////////////////////////////////////////
 //////////////////////////////////////////////////
 #include "RHEJ/currents.hh"
 //#include "ZJets/Flags.h"
 #include "RHEJ/utility.hh"
+#include "RHEJ/PDG_codes.hh"
 
 const COM looprwfactor = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4);
 //const double HVE = 246.21845810181637;
 
 #ifdef RHEJ_BUILD_WITH_QCDLOOP
 #include "qcdloop/qcdloop.h"
 #endif
 
 #include <iostream>
 
 namespace {
   // Loop integrals
   #ifdef RHEJ_BUILD_WITH_QCDLOOP
 
   COM B0DD(CLHEP::HepLorentzVector 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(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector 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(CLHEP::HepLorentzVector q1,CLHEP::HepLorentzVector q2, CLHEP::HepLorentzVector 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(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt)
   // As given in Eq. (B.2) of VDD
   {
     double q12,q22,Q2;
     CLHEP::HepLorentzVector 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();
     // std::cout<<"Higgs mass? : "<<sqrt(Q2)<<std::endl;
 
     Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
     if (mt < 0.)
       std::cerr<<"Problem in A1! mt = "<<mt<<std::endl;
     mt2=mt*mt;
 
     ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*(4.*mt2/Delta3*(Q2-q12-q22)-1.-4.*q12*q22/Delta3-12.*q12*q22*Q2/Delta3/Delta3*(q12+q22-Q2));
     ans=ans-looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))*(2.*q22/Delta3+12.*q12*q22/Delta3/Delta3*(q22-q12+Q2));
     ans=ans-looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))*(2.*q12/Delta3+12.*q12*q22/Delta3/Delta3*(q12-q22+Q2));
     ans=ans-2./Delta3/16/M_PI/M_PI*(q12+q22-Q2);
 
     //cout << "q12, q22= "<<q12<<" "<<q22<<" "<<endl;
 
     return ans;
 
   }
 
 
   COM A2(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt)
   // As given in Eq. (B.2) of VDD, but with high energy limit
   // of invariants taken.
   {
     double q12,q22,Q2;
     CLHEP::HepLorentzVector 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();
     // std::cout<<"Higgs mass Square? : "<<Q2<<std::endl;
 
     Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
     ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*(2.*mt2+1./2.*(q12+q22-Q2)+2.*q12*q22*Q2/Delta3);
     ans=ans+looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))*q22*(q22-q12-Q2)/Delta3;
     ans=ans+looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))*q12*(q12-q22-Q2)/Delta3+1./16/M_PI/M_PI;
 
     return ans;
   }
 
 #else // no QCDloop
 
   COM A1(CLHEP::HepLorentzVector, CLHEP::HepLorentzVector, double) {
     throw std::logic_error{"A1 called without QCDloop support"};
   }
 
   COM A2(CLHEP::HepLorentzVector, CLHEP::HepLorentzVector, double) {
     throw std::logic_error{"A2 called without QCDloop support"};
   }
 
 #endif
 
   void to_current(const CLHEP::HepLorentzVector & q, current & ret){
     ret[0]=q.e();
     ret[1]=q.x();
     ret[2]=q.y();
     ret[3]=q.z();
   }
+
+  constexpr double C_A = 3.;
+  constexpr double C_F = 4./3.;
+  using ParticleID = RHEJ::pid::ParticleID;
+
+  // 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)*(C_A - 1/C_A) + 1/C_A;
+  }
+  double K_g(
+      CLHEP::HepLorentzVector const & pout,
+      CLHEP::HepLorentzVector const & pin
+  ) {
+    if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
+    return K_g(pout.minus(), pin.minus());
+  }
+
 } // namespace anonymous
 
 /// @TODO move all of these functions to header? clean up in general
 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];
 }
 
 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();
 }
 
 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!?!
 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];
 }
 
 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];
 }
 
 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];
 }
 
 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];
 }
 
 double abs2(const COM & a)
 {
     return (a*conj(a)).real();
 }
 
 double vabs2(const CCurrent & cur)
 {
     return abs2(cur.c0)-abs2(cur.c1)-abs2(cur.c2)-abs2(cur.c3);
 }
 
 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));
 }
 
 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(CLHEP::HepLorentzVector p1)
 {
     //  Current goes (E,px,py,pz)
     //  std::cout<<"current = ("<<c0<<","<<c1<<","<<c2<<","<<c3<<")\n";
     //  Vector goes (px,py,pz,E)
     //  std::cout<<"vector = ("<<p1[0]<<","<<p1[1]<<","<<p1[2]<<","<<p1[3]<<")\n";
     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;
 }
 
 
 
 void j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin,current &cur) {
 
   cur[0]=0.;
   cur[1]=0.;
   cur[2]=0.;
   cur[3]=0.;
 
   double sqpop=sqrt(pout.plus());
   double sqpom=sqrt(pout.minus());
   COM poperp=pout.x()+COM(0,1)*pout.y();
 
   if (helout!=helin) {
     std::cerr<< "void j : Non-matching helicities at line " << __LINE__ << std::endl;
   } else if (helout==false) { // negative helicity
     if (pin.plus()>pin.minus()) { // if forward
       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
       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
       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
       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 j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin)
 {
     COM cur[4];
 
     cur[0]=0.;
     cur[1]=0.;
     cur[2]=0.;
     cur[3]=0.;
 
     double sqpop=sqrt(pout.plus());
     double sqpom=sqrt(pout.minus());
     COM poperp=pout.x()+COM(0,1)*pout.y();
 
     if (helout!=helin) {
         std::cerr<< "void j : Non-matching helicities\n";
     } else if (helout==false) { // negative helicity
         if (pin.plus()>pin.minus()) { // if forward
             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
             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
             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
             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 temp(cur[0],cur[1],cur[2],cur[3]);
     return temp;
 }
 
 CCurrent jio (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout)
 {
     COM cur[4];
 
     cur[0]=0.;
     cur[1]=0.;
     cur[2]=0.;
     cur[3]=0.;
 
     double sqpop=sqrt(pout.plus());
     double sqpom=sqrt(pout.minus());
     COM poperp=pout.x()+COM(0,1)*pout.y();
 
     if (helout!=helin) {
         std::cerr<< "void j : Non-matching helicities\n";
     } else if (helout==false) { // negative helicity
         if (pin.plus()>pin.minus()) { // if forward
             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
             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
             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
             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 temp(cur[0],cur[1],cur[2],cur[3]);
     return temp;
 }
 
 
 // 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;
   if(helin!=helout){
     std::cout<<__LINE__<<" "<<__FILE__<<std::endl;
   }
   double sqpop  = sqrt(pout.plus());
   double sqpom  = sqrt(pout.minus());
   COM    poperp = pout.x() + COM(0, 1) * pout.y();
 
   if (helout == false) {
 
     if (pin.plus() > pin.minus()) { // if forward
       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 {
       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 {
     if (pin.plus() > pin.minus()) { // if forward
       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 {
       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];
     }
   }
 }
 
 
 // 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(helj){
     std::cout<<__LINE__<<" "<<__FILE__<<std::endl;
   }
   // If positive helicity swap momenta
   if (heli == true) {
     HLV dummy;
     dummy = pi;
     pi = pj;
     pj = dummy;
   }
 
   double sqpjp = sqrt(pj.plus());
   double sqpjm = sqrt(pj.minus());
   double sqpip = sqrt(pi.plus());
   double sqpim = sqrt(pi.minus());
 
   COM piperp = pi.x() + COM(0,1) * pi.y();
   COM pjperp = pj.x() + COM(0,1) * pj.y();
   COM phasei = piperp / abs(piperp);
   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 (CLHEP::HepLorentzVector pi, bool heli, CLHEP::HepLorentzVector pj, bool helj)
 {
     COM cur[4];
 
     if (heli!=helj) {
         std::cerr<< "void j : Non-matching helicities\n";
     } else if (heli==true) { // negative helicity
         CLHEP::HepLorentzVector dummy;
         dummy=pi;
         pi=pj;
         pj=dummy;
     }
     double sqpjp=sqrt(pj.plus());
     double sqpjm=sqrt(pj.minus());
     double sqpip=sqrt(pi.plus());
     double sqpim=sqrt(pi.minus());
     COM piperp=pi.x()+COM(0,1)*pi.y();
     COM pjperp=pj.x()+COM(0,1)*pj.y();
     COM phasei=piperp/abs(piperp);
     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 temp(cur[0],cur[1],cur[2],cur[3]);
     return temp;
 }
 
 
 // Current Functions
 // Current for <outgoing state | mu | incoming state>
 void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
 
   cur[0] = 0.0;
   cur[1] = 0.0;
   cur[2] = 0.0;
   cur[3] = 0.0;
   if(helin){
     std::cout<<__LINE__<<" "<<__FILE__<<std::endl;
   }
   double sqpop  = sqrt(pout.plus());
   double sqpom  = sqrt(pout.minus());
   COM    poperp = pout.x() + COM(0, 1) * pout.y();
 
   if (helout == false) {
 
     if (pin.plus() > pin.minus()) { // if forward
       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 {
       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 {
     if (pin.plus() > pin.minus()) { // if forward
       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 {
       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];
     }
   }
 }
 
 namespace {
   /// @TODO unused function
   // double jM2 (CLHEP::HepLorentzVector p1out, bool hel1out, CLHEP::HepLorentzVector p1in, bool hel1in, CLHEP::HepLorentzVector p2out, bool hel2out, CLHEP::HepLorentzVector p2in, bool hel2in)
   // {
   //   CLHEP::HepLorentzVector q1=p1in-p1out;
   //   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   //   current C1,C2;
   //   j (p1out,hel1out,p1in,hel1in, C1);
   //   j (p2out,hel2out,p2in,hel2in, C2);
 
   //   std::cout << "# From Currents, C1 : ("<<C1[0]<<","<<C1[1]<<","<<C1[2]<<","<<C1[3]<<"\n";
   //   std::cout << "# From Currents, C2 : ("<<C2[0]<<","<<C2[1]<<","<<C2[2]<<","<<C2[3]<<"\n";
 
   //   COM M=cdot(C1,C2);
 
   //   return (M*conj(M)).real()/(q1.m2()*q2.m2());
   // }
 
   void jW (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin, current cur)
   {
     // NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
     // Need to swap e and nu for events with W- --> e- nubar!
     if (helin==helout && hele==helnu) {
       CLHEP::HepLorentzVector qa=pout+pe+pnu;
       CLHEP::HepLorentzVector qb=pin-pe-pnu;
       double ta(qa.m2()),tb(qb.m2());
 
       current t65,vout,vin,temp2,temp3,temp5;
       joo(pnu,helnu,pe,hele,t65);
       vout[0]=pout.e();
       vout[1]=pout.x();
       vout[2]=pout.y();
       vout[3]=pout.z();
       vin[0]=pin.e();
       vin[1]=pin.x();
       vin[2]=pin.y();
       vin[3]=pin.z();
 
       COM brac615=cdot(t65,vout);
       COM brac645=cdot(t65,vin);
 
       // prod1565 and prod6465 are zero for Ws (not Zs)!!
       // noalias(temp)=prod(trans(CurrentOutOut(pout,helout,pnu,helout)),metric);
       joo(pout,helout,pnu,helout,temp2);
       // noalias(temp2)=prod(temp,ctemp);
       COM prod1665=cdot(temp2,t65);
       // noalias(temp)=prod(trans(Current(pe,helin,pin,helin)),metric);
       // noalias(temp2)=prod(temp,ctemp);
       j(pe,helin,pin,helin,temp3);
       COM prod5465=cdot(temp3,t65);
       // noalias(temp)=prod(trans(Current(pnu,helin,pin,helin)),metric);
       // noalias(temp2)=prod(temp,ctemp);
 
       joo(pout,helout,pe,helout,temp2);
       j(pnu,helnu,pin,helin,temp3);
       j(pout,helout,pin,helin,temp5);
 
       current term1,term2,term3,sum;
       cmult(2.*brac615/ta+2.*brac645/tb,temp5,term1);
       cmult(prod1665/ta,temp3,term2);
       cmult(-prod5465/tb,temp2,term3);
 
      //    cur=((2.*brac615*Current(pout,helout,pin,helin)+prod1565*Current(pe,helin,pin,helin)+prod1665*Current(pnu,helin,pin,helin))/ta + (2.*brac645*Current(pout,helout,pin,helin)-prod5465*CurrentOutOut(pout,helout,pe,helout)-prod6465*CurrentOutOut(pout,helout,pnu,helout))/tb);
       //    cur=((2.*brac615*temp5+prod1565*temp3+prod1665*temp4)/ta + (2.*brac645*temp5-prod5465*temp1-prod6465*temp2)/tb);
       cadd(term1,term2,term3,sum);
       //    std::cout<<"sum: ("<<sum[0]<<","<<sum[1]<<","<<sum[2]<<","<<sum[3]<<")\n";
       cur[0]=sum[0];
       cur[1]=sum[1];
       cur[2]=sum[2];
       cur[3]=sum[3];
     }
   }
 
   void jWbar (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin, current cur)
   {
     // NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
     // Need to swap e and nu for events with W- --> e- nubar!
     if (helin==helout && hele==helnu) {
       CLHEP::HepLorentzVector qa=pout+pe+pnu;
       CLHEP::HepLorentzVector qb=pin-pe-pnu;
       double ta(qa.m2()),tb(qb.m2());
 
       current t65,vout,vin,temp2,temp3,temp5;
       joo(pnu,helnu,pe,hele,t65);
       vout[0]=pout.e();
       vout[1]=pout.x();
       vout[2]=pout.y();
       vout[3]=pout.z();
       vin[0]=pin.e();
       vin[1]=pin.x();
       vin[2]=pin.y();
       vin[3]=pin.z();
 
       COM brac615=cdot(t65,vout);
       COM brac645=cdot(t65,vin);
 
       // prod1565 and prod6465 are zero for Ws (not Zs)!!
       joo(pe,helout,pout,helout,temp2);  //  temp2 is <5|alpha|1>
       COM prod5165=cdot(temp2,t65);
       jio(pin,helin,pnu,helin,temp3);      // temp3 is <4|alpha|6>
       COM prod4665=cdot(temp3,t65);
 
       joo(pnu,helout,pout,helout,temp2);  // temp2 is now <6|mu|1>
       jio(pin,helin,pe,helin,temp3);        // temp3 is now <4|mu|5>
       jio(pin,helin,pout,helout,temp5);     //  temp5 is <4|mu|1>
 
       current term1,term2,term3,sum;
       cmult(-2.*brac615/ta-2.*brac645/tb,temp5,term1);
       cmult(-prod5165/ta,temp3,term2);
       cmult(prod4665/tb,temp2,term3);
 
      //    cur=((2.*brac615*Current(pout,helout,pin,helin)+prod1565*Current(pe,helin,pin,helin)+prod1665*Current(pnu,helin,pin,helin))/ta + (2.*brac645*Current(pout,helout,pin,helin)-prod5465*CurrentOutOut(pout,helout,pe,helout)-prod6465*CurrentOutOut(pout,helout,pnu,helout))/tb);
       //    cur=((2.*brac615*temp5+prod1565*temp3+prod1665*temp4)/ta + (2.*brac645*temp5-prod5465*temp1-prod6465*temp2)/tb);
       cadd(term1,term2,term3,sum);
       //    std::cout<<"term1: ("<<temp5[0]<<"  "<<temp5[1]<<"  "<<temp5[2]<<"  "<<temp5[3]<<")"<<std::endl;
       //    std::cout<<"sum: ("<<sum[0]<<","<<sum[1]<<","<<sum[2]<<","<<sum[3]<<")\n";
       cur[0]=sum[0];
       cur[1]=sum[1];
       cur[2]=sum[2];
       cur[3]=sum[3];
     }
   }
 } // namespace anonymous
 
 double jMWqQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   current mj1m,mj2p,mj2m;
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
 
   jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
   j(p2out,true,p2in,true,mj2p);
   j(p2out,false,p2in,false,mj2m);
 
   // std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n";
   // std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n";
   // std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n";
 
   // mj1m.mj2p
 
   COM Mmp=cdot(mj1m,mj2p);
 
   // mj1m.mj2m
   COM Mmm=cdot(mj1m,mj2m);
 
   // sum of spinor strings ||^2
   double a2Mmp=abs2(Mmp);
   double a2Mmm=abs2(Mmm);
 
 //   // Leave division by colour and Helicity avg until Tree files
   // Leave multi. of couplings to later
   // Multiply by Cf^2
   return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());
 
 }
 
 double jMWqQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   current mj1m,mj2p,mj2m;
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
 
   jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
   jio(p2in,true,p2out,true,mj2p);
   jio(p2in,false,p2out,false,mj2m);
 
   // std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n";
   // std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n";
   // std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n";
 
   // mj1m.mj2p
 
   COM Mmp=cdot(mj1m,mj2p);
 
   // mj1m.mj2m
   COM Mmm=cdot(mj1m,mj2m);
 
   // sum of spinor strings ||^2
   double a2Mmp=abs2(Mmp);
   double a2Mmm=abs2(Mmm);
 
 //   // Leave division by colour and Helicity avg until Tree files
   // Leave multi. of couplings to later
   // Multiply by Cf^2
   return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());
 
 }
 
 double jMWqbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   current mj1m,mj2p,mj2m;
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
 
   jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
   j(p2out,true,p2in,true,mj2p);
   j(p2out,false,p2in,false,mj2m);
 
   // std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n";
   // std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n";
   // std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n";
 
   // mj1m.mj2p
 
   COM Mmp=cdot(mj1m,mj2p);
 
   // mj1m.mj2m
   COM Mmm=cdot(mj1m,mj2m);
 
   // sum of spinor strings ||^2
   double a2Mmp=abs2(Mmp);
   double a2Mmm=abs2(Mmm);
 
 //   // Leave division by colour and Helicity avg until Tree files
   // Leave multi. of couplings to later
   // Multiply by Cf^2
   return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());
 
 }
 
 double jMWqbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   current mj1m,mj2p,mj2m;
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
 
   jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
   jio(p2in,true,p2out,true,mj2p);
   jio(p2in,false,p2out,false,mj2m);
 
   // std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n";
   // std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n";
   // std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n";
 
   // mj1m.mj2p
 
   COM Mmp=cdot(mj1m,mj2p);
 
   // mj1m.mj2m
   COM Mmm=cdot(mj1m,mj2m);
 
   // sum of spinor strings ||^2
   double a2Mmp=abs2(Mmp);
   double a2Mmm=abs2(Mmm);
 
 //   // Leave division by colour and Helicity avg until Tree files
   // Leave multi. of couplings to later
   // Multiply by Cf^2
   return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());
 
 }
 
 double jMWqg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qg->qenug scattering
 // p1: quark
 // p2: gluon
 {
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   current mj1m,mj2p,mj2m;
 
   jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
 
   j(p2out,true,p2in,true,mj2p);
   j(p2out,false,p2in,false,mj2m);
 
   // mj1m.mj2p
 
   COM Mmp=cdot(mj1m,mj2p);
 
   // mj1m.mj2m
   COM Mmm=cdot(mj1m,mj2m);
 
-  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.;
+  const double K = K_g(p2out, p2in);
 
   // sum of spinor strings ||^2
   double a2Mmp=abs2(Mmp);
   double a2Mmm=abs2(Mmm);
-  double sst = nonflipcolourmult*(a2Mmp+a2Mmm);
+  double sst = K/C_A*(a2Mmp+a2Mmm);
   // double sstsave=sst;
 
 //   // Leave division by colour and Helicity avg until Tree files
   // Leave multi. of couplings to later
   // Multiply by Cf*Ca=4
   return 4.*sst/(q1.m2()*q2.m2());
 
 }
 
 double jMWqbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qg->qenug scattering
 // p1: quark
 // p2: gluon
 {
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   current mj1m,mj2p,mj2m;
 
   jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m);
 
   j(p2out,true,p2in,true,mj2p);
   j(p2out,false,p2in,false,mj2m);
 
   // mj1m.mj2p
 
   COM Mmp=cdot(mj1m,mj2p);
 
   // mj1m.mj2m
   COM Mmm=cdot(mj1m,mj2m);
 
-  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.;
+  const double K = K_g(p2out, p2in);
 
   // sum of spinor strings ||^2
   double a2Mmp=abs2(Mmp);
   double a2Mmm=abs2(Mmm);
-  double sst = nonflipcolourmult*(a2Mmp+a2Mmm);
+  double sst = K/C_A*(a2Mmp+a2Mmm);
   // double sstsave=sst;
 
 //   // Leave division by colour and Helicity avg until Tree files
   // Leave multi. of couplings to later
   // Multiply by Cf*Ca=4
   return 4.*sst/(q1.m2()*q2.m2());
 
 }
 
 double jM2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //  std::cerr<<"Current: "<<p1out<<"  "<<p1in<<"  "<<p2out<<"  "<<p2in<<std::endl;
 
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   //  std::cerr<<"Current:  "<<q1.m2()<<"  "<<q2.m2()<<std::endl;
   current mj1m,mj1p,mj2m,mj2p;
   j(p1out,true,p1in,true,mj1p);
   j(p1out,false,p1in,false,mj1m);
   j(p2out,true,p2in,true,mj2p);
   j(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 (4./3.)*(4./3.)*(sst)/(q1.m2()*q2.m2());
 }
 
 double jM2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   current mj1m,mj1p,mj2m,mj2p;
   j(p1out,true,p1in,true,mj1p);
   j(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 (4./3.)*(4./3.)*(sumsq)/(q1.m2()*q2.m2());
 }
 
 double jM2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector 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 (4./3.)*(4./3.)*(sumsq)/(q1.m2()*q2.m2());
 }
 
 double jM2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qg scattering
 // p1: quark
 // p2: gluon
 {
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
 
   current mj1m,mj1p,mj2m,mj2p;
   j(p1out,true,p1in,true,mj1p);
   j(p1out,false,p1in,false,mj1m);
   j(p2out,true,p2in,true,mj2p);
   j(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 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.;
+  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 = nonflipcolourmult*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
+  double sst = K/C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
   // double sstsave=sst;
 
    // std::cout <<"ratio: "<<sst/sstsave<<std::endl;
   // Cf*Ca=4
   return 4.*sst/(q1.m2()*q2.m2());
 
 }
 
 double jM2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qg scattering
 // p1: quark
 // p2: gluon
 {
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
 
   current mj1m,mj1p,mj2m,mj2p;
   jio(p1in,true,p1out,true,mj1p);
   jio(p1in,false,p1out,false,mj1m);
   j(p2out,true,p2in,true,mj2p);
   j(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 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.;
+  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 = nonflipcolourmult*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
+  double sst = K/C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
   // double sstsave=sst;
 
    // std::cout <<"ratio: "<<sst/sstsave<<std::endl;
   // Cf*Ca=4
   return 4.*sst/(q1.m2()*q2.m2());
 
 }
 
 double jM2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for gg scattering
 // p1: gluon
 // p2: gluon
 {
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector q2=-(p2in-p2out);
 
   current mj1m,mj1p,mj2m,mj2p;
   j(p1out,true,p1in,true,mj1p);
   j(p1out,false,p1in,false,mj1m);
   j(p2out,true,p2in,true,mj2p);
   j(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 ratio1; // p1-/pa- in the notes
-  // if (p1in.plus()>1.) // if the gluon is the positive. Should have been a
-  //          // test against 0, but 1. is better
-  if (p1in.pz()>0.) // a much better test
-    ratio1=p1out.plus()/p1in.plus();
-  else // the gluon is the negative
-    ratio1=p1out.minus()/p1in.minus();
-  double nonflipcolourmult1=(1.-1./9.)/2.*(ratio1+1./ratio1)+1./9.;
-
-  double ratio2; // p2-/pb- in the notes
-  if (p2in.pz()>0.) // a much better test
-    ratio2=p2out.plus()/p2in.plus();
-  else // the gluon is the negative
-    ratio2=p2out.minus()/p2in.minus();
-  double nonflipcolourmult2=(1.-1./9.)/2.*(ratio2+1./ratio2)+1./9.;
+  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 = nonflipcolourmult1*nonflipcolourmult2*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
+  double sst = K_g1/C_A*K_g2/C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm);
   // double sstsave=sst;
    // std::cout <<"ratio: "<<sst/sstsave<<std::endl;
   // Ca*Ca=9
   return 9.*sst/(q1.m2()*q2.m2());
 
 }
 
 namespace {
   /// @TODO what was this intended to do?
   // double MH2helper(current C1, current C2, current q1, current q2)
   // {
   //   COM M;
   //   COM temp1,temp2;
   //   // First the C1.q2 * C2.q1 - part
   //   temp1=cdot(C1,q2);
   //   temp2=cdot(C2,q1);
   //   M=temp1*temp2;
 
   //   // Then the C1.C2 * q1.q2
   //   temp1=cdot(C1,C2);
   //   temp2=cdot(q1,q2);
   //   M-=temp1*temp2;
 
   //   return (M*conj(M)).real();
   // }
 
   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))/(6*M_PI*v);
     }
     else {
       CLHEP::HepLorentzVector 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
       // std::cout<<"A1 : " << A1(-vq1,vq2)<<std::endl;
       // std::cout<<"A2 : " << A2(-vq1,vq2)<<std::endl;
       if(!(incBot)) // Factor is because 4 mt^2 g^2/v A1 -> 16 pi mt^2/v alphas, and we divide by a factor 4 at the amp sqaured level later which I absorb here (i.e. I divide by 2)
         return 8.*M_PI*mt*mt/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt));
       else
         return 8.*M_PI*mt*mt/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt))
              + 8.*M_PI*mb*mb/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)-cdot(C1,C2)*A2(-vq1,vq2,mb));
     }
   }
 } // namespace anonymous
 
 double MH2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
               CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
               CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2,
               double mt, bool incBot, double mb)
 {
 //   CLHEP::HepLorentzVector q1=p1in-p1out;
 //   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   current j1p,j1m,j2p,j2m, q1v, q2v;
 
   j (p1out,true,p1in,true,j1p);
   j (p1out,false,p1in,false,j1m);
 
   j (p2out,true,p2in,true,j2p);
   j (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 (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
   return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
 }
 
 double MH2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
 {
 //   CLHEP::HepLorentzVector q1=p1in-p1out;
 //   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   current j1p,j1m,j2p,j2m,q1v,q2v;
 
   j (p1out,true,p1in,true,j1p);
   j (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 (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
   return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
 }
 
 double MH2qbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
 {
 //   CLHEP::HepLorentzVector q1=p1in-p1out;
 //   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   current j1p,j1m,j2p,j2m,q1v,q2v;
 
   jio (p1in,true,p1out,true,j1p);
   jio (p1in,false,p1out,false,j1m);
 
   j (p2out,true,p2in,true,j2p);
   j (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 (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
   return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
 }
 
 double MH2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
 {
 //   CLHEP::HepLorentzVector q1=p1in-p1out;
 //   CLHEP::HepLorentzVector q2=-(p2in-p2out);
   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 (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
   return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
 }
 
 double MH2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector 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;
 
   j (p1out,true,p1in,true,j1p);
   j (p1out,false,p1in,false,j1m);
 
   j (p2out,true,p2in,true,j2p);
   j (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);
 
   //cout << "Bits in MH2qg: " << Mpp << " " << Mpm << " " << Mmp << " " << Mmm << endl;
 
-  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();
+  const double K = K_g(p2out, p2in);
 
-  double nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.;
-  double sst=nonflipcolourmult*(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm));
+  double sst=K/C_A*(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm));
 
   // Cf*Ca=4
   // return 4.*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
   return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
 }
 
 double MH2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector 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);
 
   j (p2out,true,p2in,true,j2p);
   j (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);
 
-  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.;
-  MH2sum*=nonflipcolourmult;
+  const double K = K_g(p2out, p2in);
+  MH2sum*=K/C_A;
 
   // Cf*Ca=4
   // return 4.*MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
   return MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
 }
 
 double MH2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector 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;
 
   j (p1out,true,p1in,true,j1p);
   j (p1out,false,p1in,false,j1m);
 
   j (p2out,true,p2in,true,j2p);
   j (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);
 
-  double ratio1; // p1-/pa- in the notes
-  // if (p1in.plus()>0) // if the gluon is the positive
-  if (p1in.pz()>0.) // if the gluon is the positive
-    ratio1=p1out.plus()/p1in.plus();
-  else // the gluon is the negative
-    ratio1=p1out.minus()/p1in.minus();
-
-  double nonflipcolourmult1=(1.-1./9.)/2.*(ratio1+1./ratio1)+1./9.;
-
-  double ratio2; // 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
-    ratio2=p2out.plus()/p2in.plus();
-  else // the gluon is the negative
-    ratio2=p2out.minus()/p2in.minus();
+  const double K_g1 = K_g(p1out, p1in);
+  const double K_g2 = K_g(p2out, p2in);
 
-  double nonflipcolourmult2=(1.-1./9.)/2.*(ratio2+1./ratio2)+1./9.;
-  MH2sum*=nonflipcolourmult1*nonflipcolourmult2;
+  MH2sum*=K_g1/C_A*K_g2/C_A;
 
   // Ca*Ca=9
   // return 9.*MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
   return MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
 }
 
 // // Z's stuff
 // void jZ(HLV pin, HLV pout, HLV pem, HLV pep, bool HelPartons, bool HelLeptons, current cur) {
 
 //  // Init current to zero
 //  cur[0] = 0.0;
 //  cur[1] = 0.0;
 //  cur[2] = 0.0;
 //  cur[3] = 0.0;
 
 //  // Temporary variables
 //  COM temp;
 //  current Term_1, Term_2, Term_3, Term_4, J_temp, TempCur1, TempCur2;
 
 //  // Momentum of virtual gluons aroun weak boson emission site
 //  HLV qa = pout + pep + pem;
 //  HLV qb = pin  - pep - pem;
 
 //  double ta = qa.m2();
 //  double tb = qb.m2();
 
 //  // Out-Out currents:
 //  current Em_Ep, Out_Em, Out_Ep;
 
 //  // Other currents:
 //  current Out_In, Em_In, Ep_In;
 
 //  joi(pout, HelPartons, pin, HelPartons, Out_In);
 //  joi(pem,  HelLeptons, pin, HelPartons, Em_In);
 //  joi(pep,  HelLeptons, pin, HelPartons, Ep_In);
 
 //  joo(pem,  HelLeptons, pep, HelLeptons, Em_Ep);
 //  joo(pout, HelPartons, pem, HelLeptons, Out_Em);
 //  joo(pout, HelPartons, pep, HelLeptons, Out_Ep);
 
 //  if (HelLeptons == HelPartons) {
 
 //    temp = 2.0 * cdot(pout, Em_Ep);
 //    cmult(temp / ta, Out_In, Term_1);
 
 //    temp = cdot(Out_Em, Em_Ep);
 //    cmult(temp / ta , Em_In, Term_2);
 
 //    temp = 2.0 * cdot(pin, Em_Ep);
 //    cmult(temp / tb, Out_In, Term_3);
 
 //    temp = -cdot(Ep_In, Em_Ep);
 //    cmult(temp / tb, Out_Ep, Term_4);
 
 //    cadd(Term_1, Term_2, Term_3, Term_4, J_temp);
 
 //    cur[0] = J_temp[0];
 //    cur[1] = J_temp[1];
 //    cur[2] = J_temp[2];
 //    cur[3] = J_temp[3];
 //  }
 
 //  else {
 //    if (HelPartons == true) {
 //      temp = 2.0 * cdot(pout, Em_Ep);
 //      cmult(temp / ta, Out_In, Term_1);
 
 //      joo(pout, true, pep, true, TempCur1);
 //      joi(pep,  true, pin, true, TempCur2);
 
 //      temp = cdot(TempCur1, Em_Ep);
 //      cmult(temp / ta , TempCur2, Term_2);
 
 //      temp = 2.0 * cdot(pin, Em_Ep);
 //      cmult(temp / tb, Out_In, Term_3);
 
 //      joo(pout, true, pem, true, TempCur1);
 //      joi(pem,  true, pin, true, TempCur2);
 
 //      temp = -cdot(TempCur2, Em_Ep);
 //      cmult(temp / tb, TempCur1, Term_4);
 
 //      cadd(Term_1, Term_2, Term_3, Term_4, J_temp);
 
 //      cur[0] = J_temp[0];
 //      cur[1] = J_temp[1];
 //      cur[2] = J_temp[2];
 //      cur[3] = J_temp[3];
 //    }
 
 //    else {
 //      temp = 2.0 * cdot(pout, Em_Ep);
 //      cmult(temp / ta, Out_In, Term_1);
 
 //      joo(pout, false, pep, false, TempCur1);
 //      joi(pep,  false, pin, false, TempCur2);
 
 //      temp = cdot(TempCur1, Em_Ep);
 //      cmult(temp / ta, TempCur2, Term_2);
 
 //      temp = 2.0 * cdot(pin, Em_Ep);
 //      cmult(temp / tb, Out_In, Term_3);
 
 //      joo(pout, false, pem, false, TempCur1);
 //      joi(pem,  false, pin, false, TempCur2);
 
 //      temp = -cdot(TempCur2, Em_Ep);
 //      cmult(temp / tb, TempCur1, Term_4);
 
 //      cadd(Term_1, Term_2, Term_3, Term_4, J_temp);
 
 //      cur[0] = J_temp[0];
 //      cur[1] = J_temp[1];
 //      cur[2] = J_temp[2];
 //      cur[3] = J_temp[3];
 //    }
 
 //  }
 // }
 
 // void jZbar(HLV pin, HLV pout, HLV pem, HLV pep, bool HelPartons, bool HelLeptons, current cur) {
 
 //  // Init current to zero
 //  cur[0] = 0.0;
 //  cur[1] = 0.0;
 //  cur[2] = 0.0;
 //  cur[3] = 0.0;
 
 //  // Temporary variables
 //  COM temp;
 //  current Term_1, Term_2, Term_3, Term_4, J_temp, TempCur1, TempCur2;
 
 //  // Transfered 4-momenta
 //  HLV qa = pout + pep + pem;
 //  HLV qb = pin  - pep - pem;
 
 //  // The square of the transfered 4-momenta
 //  double ta = qa.m2();
 //  double tb = qb.m2();
 
 //  // Out-Out currents:
 //  current Em_Ep, Em_Out, Ep_Out;
 
 //  // In-Out currents:
 //  current In_Out, In_Em, In_Ep;
 
 //  // Safe to use the currents since helicity structure is ok
 //  if (HelPartons == HelLeptons) {
 //    jio(pin, HelPartons, pout, HelPartons, In_Out);
 //    joo(pem, HelLeptons, pep,  HelLeptons, Em_Ep);
 //    jio(pin, HelPartons, pem,  HelLeptons, In_Em);
 //    jio(pin, HelPartons, pep,  HelLeptons, In_Ep);
 //    joo(pem, HelLeptons, pout, HelPartons, Em_Out);
 //    joo(pep, HelLeptons, pout, HelPartons, Ep_Out);
 //  }
 
 //  else {
 //    jio(pin, HelPartons, pout, HelPartons, In_Out);
 //    joo(pem, HelLeptons, pep,  HelLeptons, Em_Ep);
 
 //    In_Em[0] = 0.0;
 //    In_Em[1] = 0.0;
 //    In_Em[2] = 0.0;
 //    In_Em[3] = 0.0;
 
 //    In_Ep[0] = 0.0;
 //    In_Ep[1] = 0.0;
 //    In_Ep[2] = 0.0;
 //    In_Ep[3] = 0.0;
 
 //    Em_Out[0] = 0.0;
 //    Em_Out[1] = 0.0;
 //    Em_Out[2] = 0.0;
 //    Em_Out[3] = 0.0;
 
 //    Ep_Out[0] = 0.0;
 //    Ep_Out[1] = 0.0;
 //    Ep_Out[2] = 0.0;
 //    Ep_Out[3] = 0.0;
 //  }
 
 //  if (HelLeptons == HelPartons) {
 
 //    temp = 2.0 * cdot(pout, Em_Ep);
 //    cmult(temp / ta, In_Out, Term_1);
 
 //    temp = cdot(Ep_Out, Em_Ep);
 //    cmult(temp / ta, In_Ep, Term_2);
 
 //    temp = 2.0 * cdot(pin, Em_Ep);
 //    cmult(temp / tb, In_Out, Term_3);
 
 //    temp = - cdot(In_Em, Em_Ep);
 //    cmult(temp / tb, Em_Out, Term_4);
 
 //    cadd(Term_1, Term_2, Term_3, Term_4, J_temp);
 
 //    cur[0] = J_temp[0];
 //    cur[1] = J_temp[1];
 //    cur[2] = J_temp[2];
 //    cur[3] = J_temp[3];
 //  }
 
 //  else {
 //    if (HelPartons == true) {
 
 //      temp = 2.0 * cdot(pout, Em_Ep);
 //      cmult(temp / ta, In_Out, Term_1);
 
 //      joo(pem, true, pout, true, TempCur1);
 //      jio(pin, true, pem,  true, TempCur2);
 
 //      temp = cdot(TempCur1, Em_Ep);
 //      cmult(temp / ta , TempCur2, Term_2);
 
 //      temp = 2.0 * cdot(pin, Em_Ep);
 //      cmult(temp / tb, In_Out, Term_3);
 
 //      joo(pep, true, pout, true, TempCur1);
 //      jio(pin, true, pep,  true, TempCur2);
 
 //      temp = - cdot(TempCur2, Em_Ep);
 //      cmult(temp / tb, TempCur1, Term_4);
 
 //      cadd(Term_1, Term_2, Term_3, Term_4, J_temp);
 
 //      cur[0] = J_temp[0];
 //      cur[1] = J_temp[1];
 //      cur[2] = J_temp[2];
 //      cur[3] = J_temp[3];
 //    }
 
 //    else {
 
 //      temp = 2.0 * cdot(pout, Em_Ep);
 //      cmult(temp / ta, In_Out, Term_1);
 
 //      joo(pem, false, pout, false, TempCur1);
 //      jio(pin, false, pem,  false, TempCur2);
 
 //      temp = cdot(TempCur1, Em_Ep);
 //      cmult(temp / ta , TempCur2, Term_2);
 
 //      temp = 2.0 * cdot(pin, Em_Ep);
 //      cmult(temp / tb, In_Out, Term_3);
 
 //      joo(pep, false, pout, false, TempCur1);
 //      jio(pin, false, pep,  false, TempCur2);
 
 //      temp = - cdot(TempCur2, Em_Ep);
 //      cmult(temp / tb, TempCur1, Term_4);
 
 //      cadd(Term_1, Term_2, Term_3, Term_4, J_temp);
 
 //      cur[0] = J_temp[0];
 //      cur[1] = J_temp[1];
 //      cur[2] = J_temp[2];
 //      cur[3] = J_temp[3];
 //    }
 //  }
 // }
 
 // // Progagators
 // COM PZ(double s) {
 
 //  double MZ, GammaZ;
 
 //  MZ       = 9.118800e+01;  // Mass of the mediating gauge boson
 //  GammaZ   = 2.441404e+00;  // Z peak width
 
 //  // Return Z Prop value
 //  return 1.0 / (s - MZ * MZ + COM(0.0, 1.0) * GammaZ * MZ);
 // }
 
 // COM PG(double s) {
 //  return 1.0 / s;
 // }
 
 // // Non-gluonic with pa emitting
 // std::vector <double> jMZqQ (HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, std::vector <double> VProducts, std::vector < std::vector <double> > Virtuals, int aptype, int bptype, bool UseVirtuals, bool BottomLineEmit) {
 
 //  std::vector <double> ScaledWeights;
 
 //  double Sum;
 
 //  // Propagator factors
 //  COM PZs = PZ((pep + pem).m2());
 //  COM PGs = PG((pep + pem).m2());
 
 
 //  // Emitting current initialisation
 //  current j1pptop, j1pmtop; // Emission from top line
 //  current j1ppbot, j1pmbot; // Emission from bottom line
 
 //  // Non-emitting current initialisation
 //  current j2ptop, j2mtop;   // Emission from top line
 //  current j2pbot, j2mbot;   // Emission from bottom line
 
 //  // Currents for top emission
 //  // Upper current calculations
 //  // if a is a quark
 //  if (aptype > 0) {
 //    jZ(pa, p1, pem, pep, true,  true,  j1pptop);
 //    jZ(pa, p1, pem, pep, true,  false, j1pmtop);
 //  }
 //  // if a is an antiquark
 //  else {
 //    jZbar(pa, p1, pem, pep, true,  true,  j1pptop);
 //    jZbar(pa, p1, pem, pep, true,  false, j1pmtop);
 //  }
 
 //  // Lower current calculations
 //  // if b is a quark
 //  if (bptype > 0) {
 //    joi(p2, true,  pb, true,  j2ptop);
 //    joi(p2, false, pb, false, j2mtop);
 //  }
 //  // if b is an antiquark
 //  else {
 //    jio(pb, true,  p2, true,  j2ptop);
 //    jio(pb, false, p2, false, j2mtop);
 //  }
 
 //  // Currents for bottom emission
 //  // Lower current calculations
 //  if (bptype > 0) {
 //    jZ(pb, p2, pem, pep, true,  true,  j1ppbot);
 //    jZ(pb, p2, pem, pep, true,  false, j1pmbot);
 //  }
 //  else {
 //    jZbar(pb, p2, pem, pep, true,  true,  j1ppbot);
 //    jZbar(pb, p2, pem, pep, true,  false, j1pmbot);
 //  }
 
 //  // Upper current calculations
 //  if (aptype > 0) {
 //    joi(p1, true,  pa, true,  j2pbot);
 //    joi(p1, false, pa, false, j2mbot);
 //  }
 //  else {
 //    jio(pa, true,  p1, true,  j2pbot);
 //    jio(pa, false, p1, false, j2mbot);
 //  }
 
 //  COM Coeff[2][8];
 
 //  if (!Interference) {
 
 //    double ZCharge_a_P = Zq(aptype, true);
 //    double ZCharge_a_M = Zq(aptype, false);
 //    double ZCharge_b_P = Zq(bptype, true);
 //    double ZCharge_b_M = Zq(bptype, false);
 
 //    if (BottomLineEmit) {
 //      // Emission from top-line quark (pa/p1 line)
 //      Coeff[0][0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pptop, j2ptop);
 //      Coeff[0][1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pptop, j2mtop);
 //      Coeff[0][2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pmtop, j2ptop);
 //      Coeff[0][3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pmtop, j2mtop);
 //      Coeff[0][4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop));
 //      Coeff[0][5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop));
 //      Coeff[0][6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop));
 //      Coeff[0][7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop));
 //    }
 
 //    else {
 //      // Emission from bottom-line quark (pb/p2 line)
 //      Coeff[1][0] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) *      cdot(j1ppbot, j2pbot);
 //      Coeff[1][7] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) *      cdot(j1ppbot, j2mbot);
 //      Coeff[1][2] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) *      cdot(j1pmbot, j2pbot);
 //      Coeff[1][5] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) *      cdot(j1pmbot, j2mbot);
 //      Coeff[1][4] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2pbot));
 //      Coeff[1][3] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2mbot));
 //      Coeff[1][6] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2pbot));
 //      Coeff[1][1] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2mbot));
 //    }
 //  }
 
 //  // Else calculate all the possiblities
 //  else {
 
 //    double ZCharge_a_P = Zq(aptype, true);
 //    double ZCharge_a_M = Zq(aptype, false);
 //    double ZCharge_b_P = Zq(bptype, true);
 //    double ZCharge_b_M = Zq(bptype, false);
 
 //    // Emission from top-line quark (pa/p1 line)
 //    Coeff[0][0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pptop, j2ptop);
 //    Coeff[0][1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pptop, j2mtop);
 //    Coeff[0][2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pmtop, j2ptop);
 //    Coeff[0][3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pmtop, j2mtop);
 //    Coeff[0][4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop));
 //    Coeff[0][5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop));
 //    Coeff[0][6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop));
 //    Coeff[0][7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop));
 
 //    // Emission from bottom-line quark (pb/p2 line)
 //    Coeff[1][0] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) *      cdot(j1ppbot, j2pbot);
 //    Coeff[1][7] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) *      cdot(j1ppbot, j2mbot);
 //    Coeff[1][2] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) *      cdot(j1pmbot, j2pbot);
 //    Coeff[1][5] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) *      cdot(j1pmbot, j2mbot);
 //    Coeff[1][4] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2pbot));
 //    Coeff[1][3] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2mbot));
 //    Coeff[1][6] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2pbot));
 //    Coeff[1][1] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2mbot));
 //  }
 
 //  // Find the numbers of scales
 //  int ScaleCount;
 //  #if calcscaleunc
 //    ScaleCount = 20;
 //  #else
 //    ScaleCount = 1;
 //  #endif
 
 //  // For each scale...
 //  for (int j = 0; j < ScaleCount; j++) {
 
 //    Sum = 0.0;
 
 //    // If we want to compare back to the W's code only emit from one quark and only couple to left handed particles
 //    // virtuals arent here since they are calculated and included in weight() call.
 //    if (!Interference) {
 
 //      if (BottomLineEmit) for (int i = 0; i < 8; i++) Sum += abs2(Coeff[1][i]) * VProducts.at(1);
 //      else                for (int i = 0; i < 8; i++) Sum += abs2(Coeff[0][i]) * VProducts.at(0);
 //    }
 
 //    // Else work out the full interference
 //    else {
 
 //      // For the full calculation...
 //      if (UseVirtuals) {
 //        for (int i = 0; i < 8; i++) {
 //          Sum += abs2(Coeff[0][i])  * VProducts.at(0) * Virtuals.at(j).at(0)
 //               + abs2(Coeff[1][i])  * VProducts.at(1) * Virtuals.at(j).at(1)
 //                     + 2.0 * real(Coeff[0][i] * conj(Coeff[1][i])) * VProducts.at(2) * Virtuals.at(j).at(2);
 //        }
 //      }
 
 //      // For the tree level calculation...
 //      else {
 //        for (int i = 0; i < 8; i++) {
 //          Sum += abs2(Coeff[0][i])  * VProducts.at(0)
 //               + abs2(Coeff[1][i])  * VProducts.at(1)
 //        + 2.0 * real(Coeff[0][i] * conj(Coeff[1][i])) * VProducts.at(2);
 //        }
 //      }
 //    }
 
 //    // Add this to the vector to be returned with the other factors of C_A and the helicity sum/average factors.
 //    ScaledWeights.push_back(Sum / 18.0);
 //  }
 
 //  // Return all the scale values
 //  return ScaledWeights;
 // }
 
 // // Semi-gluonic with pa emitting
 // std::vector <double> jMZqg (HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, std::vector <double> VProducts, std::vector < std::vector <double> > Virtuals, int aptype, int bptype, bool UseVirtuals, bool BottomLineEmit) {
 
 //  COM Coeff[8];
 
 //  double Sum;
 
 //  std::vector <double> ScaledWeights;
 
 //  COM PZs = PZ((pep + pem).m2());
 //  COM PGs = PG((pep + pem).m2());
 
 //  // Emitting current initialisation - Emission from top line
 //  current j1pptop, j1pmtop;
 
 //  // Non-emitting current initialisation - Emission from top line
 //  current j2ptop, j2mtop;
 
 //  // Currents for top emission
 //  // Upper current calculations
 //  if (aptype > 0) {
 //    jZ   (pa, p1, pem, pep, true,  true,  j1pptop);
 //    jZ   (pa, p1, pem, pep, true,  false, j1pmtop);
 //  }
 //  else {
 //    jZbar(pa, p1, pem, pep, true,  true,  j1pptop);
 //    jZbar(pa, p1, pem, pep, true,  false, j1pmtop);
 //  }
 
 //  // Lower current calculations
 //  joi(p2, true,  pb, true,  j2ptop);
 //  joi(p2, false, pb, false, j2mtop);
 
 //  // Calculate all the possiblities
 //  double ZCharge_a_P = Zq(aptype, true);
 //  double ZCharge_a_M = Zq(aptype, false);
 
 //  // Emission from top-line quark (pa/p1 line)
 //  Coeff[0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pptop, j2ptop);
 //  Coeff[1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pptop, j2mtop);
 //  Coeff[2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pmtop, j2ptop);
 //  Coeff[3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) *      cdot(j1pmtop, j2mtop);
 //  Coeff[4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop));
 //  Coeff[5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop));
 //  Coeff[6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop));
 //  Coeff[7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop));
 
 //  // Calculate gluon colour accelerated factor
 //  double CAMFactor, z;
 
 //  // If b is a forward moving gluon define z (C.F. multiple jets papers)
 //  if (pb.pz() > 0) z = p2.plus()  / pb.plus();
 //  else             z = p2.minus() / pb.minus();
 
 //  CAMFactor = (1.0 - 1.0 / 9.0) / 2.0 * (z + 1.0 / z) + 1.0 / 9.0;
 
 //  // Find the numbers of scales
 //  int ScaleCount;
 //  #if calcscaleunc
 //    ScaleCount = 20;
 //  #else
 //    ScaleCount = 1;
 //  #endif
 
 //  // For each scale...
 //  for (int j = 0; j < ScaleCount; j++) {
 
 //    Sum = 0.0;
 
 //    // If we dont want the interference
 //    if (!Interference) for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0);
 
 //    // Else work out the full interference
 //    else {
 //      if (UseVirtuals) {
 //        for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0) * Virtuals.at(j).at(0);
 //      }
 //      else {
 //        for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0);
 //      }
 //    }
 
 //    // Add this to the vector to be returned with the other factors of C_A, the colour accelerated factor and the helicity sum/average factors.: (4/3)*3/32
 //    ScaledWeights.push_back(CAMFactor * Sum / 8.0);
 //  }
 
 //  return ScaledWeights;
 // }
 
 // // Electroweak Charge Functions
 // double Zq (int PID, bool Helcitiy) {
 
 //  double temp;
 
 //  // Positive Spin
 //  if (Helcitiy == true) {
 //    if (PID ==   1 || PID ==  3 || PID ==  5) temp = (+ 1.0 * stw2 / 3.0) / ctw;
 //    if (PID ==   2 || PID ==  4)              temp = (- 2.0 * stw2 / 3.0) / ctw;
 //    if (PID ==  -1 || PID == -3 || PID == -5) temp = (- 1.0 * stw2 / 3.0) / ctw;
 //    if (PID ==  -2 || PID == -4)              temp = (+ 2.0 * stw2 / 3.0) / ctw;
 
 //    // If electron or positron
 //    if (PID ==  7 || PID == -7) temp = Zep;
 //  }
 
 //  // Negative Spin
 //  else {
 //    if (PID ==  1 || PID ==  3 || PID ==  5) temp = (-0.5 + 1.0 * stw2 / 3.0) / ctw;
 //    if (PID ==  2 || PID ==  4)              temp = ( 0.5 - 2.0 * stw2 / 3.0) / ctw;
 //    if (PID == -1 || PID == -3 || PID == -5) temp = ( 0.5 - 1.0 * stw2 / 3.0) / ctw;
 //    if (PID == -2 || PID == -4)              temp = (-0.5 + 2.0 * stw2 / 3.0) / ctw;
 
 //    // If electron or positron
 //    if (PID ==  7 || PID == -7) temp = Zem;
 
 //  }
 
 //  return temp;
 // }
 
 // double Gq (int PID) {
 
 //  if (!VirtualPhoton) return 0.0;
 
 //  if (PID == -1) return  1.0 * ee / 3.0;
 //  if (PID == -2) return -2.0 * ee / 3.0;
 //  if (PID == -3) return  1.0 * ee / 3.0;
 //  if (PID == -4) return -2.0 * ee / 3.0;
 //  if (PID == -5) return  1.0 * ee / 3.0;
 //  if (PID ==  1) return -1.0 * ee / 3.0;
 //  if (PID ==  2) return  2.0 * ee / 3.0;
 //  if (PID ==  3) return -1.0 * ee / 3.0;
 //  if (PID ==  4) return  2.0 * ee / 3.0;
 //  if (PID ==  5) return -1.0 * ee / 3.0;
 
 //  std::cout << "ERROR! No Electroweak Charge Found at line " << __LINE__ << "..." << std::endl;
 //  return 0.0;
 // }
 
 
 namespace {
   CCurrent jH (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
   {
 
       CCurrent j2 = j(pout,helout,pin,helin);
       CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());
 
       if(mt == infinity)
         return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*v);
       else
       {
         if(incBot)
           return (-16.*M_PI*mb*mb/v*j2.dot(q1)*jq2*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j2*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt));
         else
           return (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt));
       }
   }
 
   CCurrent jioH (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector 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)/(3*M_PI*v);
       else
       {
         if(incBot)
           return (-16.*M_PI*mb*mb/v*j2.dot(q1)*jq2*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j2*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt));
         else
           return (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt));
       }
   }
 
   CCurrent jHtop (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb)
   {
 
       CCurrent j1 = j(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)/(3*M_PI*v);
       else
       {
         if(incBot)
           return (-16.*M_PI*mb*mb/v*j1.dot(q2)*jq1*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j1*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt));
         else
           return (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt));
       }
   }
 
 
   CCurrent jioHtop (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector 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)/(3*M_PI*v);
       else
       {
         if(incBot)
           return (-16.*M_PI*mb*mb/v*j1.dot(q2)*jq1*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j1*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt));
         else
           return (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt));
       }
   }
 } // namespace anonymous
 
 double jM2unogqHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
     //  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,mjH2m,mjH2p;
     mj1p=j(p1out,true,p1in,true);
     mj1m=j(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);
 
     //  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=j(pg,true,p1in,true);
     jgam=j(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.*p1i*MHmm)/(p1in-pg).m2();
     U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
     U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(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)/(q2.m2()*qH2.m2());
 
   // 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=qH1.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 jM2unogqbarHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
     //  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,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.*p1i*MHmm)/(p1in-pg).m2();
     U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
     U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(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)/(q2.m2()*qH2.m2());
 
   // 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=qH1.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 jM2unogqHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
     //  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,mjH2m,mjH2p;
     mj1p=j(p1out,true,p1in,true);
     mj1m=j(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=j(pg,true,p1in,true);
     jgam=j(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.*p1i*MHmm)/(p1in-pg).m2();
     U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
     U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(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)/(q2.m2()*qH2.m2());
 
   // 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=qH1.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 jM2unogqbarHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
     //  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,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.*p1i*MHmm)/(p1in-pg).m2();
     U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
     U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(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)/(q2.m2()*qH2.m2());
 
   // 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=qH1.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 jM2unogqHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
     //  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,mjH2m,mjH2p;
     mj1p=j(p1out,true,p1in,true);
     mj1m=j(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=j(pg,true,p1in,true);
     jgam=j(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.*p1i*MHmm)/(p1in-pg).m2();
     U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
     U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(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)/(q2.m2()*qH2.m2());
 
   // 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=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
     //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();
+    const double K = K_g(p2out, p2in);
 
-    double nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.;
-
-    return ampsq*nonflipcolourmult*9./4.; //ca/cf = 9/4
+    return ampsq*K/C_A*9./4.; //ca/cf = 9/4
 }
 
 double jM2unogqbarHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
     //  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,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.*p1i*MHmm)/(p1in-pg).m2();
     U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
     U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(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)/(q2.m2()*qH2.m2());
 
   // 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=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
     //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.;
+    const double K = K_g(p2out, p2in);
 
-    return ampsq*nonflipcolourmult*9./4.; //ca/cf = 9/4
+    return ampsq*K/C_F;
 }
 
 double jM2unobqHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
     // 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 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=j(p2out,true,p2in,true);
     mj2m=j(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=j(pg,true,p2in,true);
     jgbm=j(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.*p2i*MHmm)/(p2in-pg).m2();
     U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2();
     U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(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)/(q1.m2()*qH1.m2());
 
   // 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=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 jM2unobqbarHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
 
     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 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=j(p2out,true,p2in,true);
     mj2m=j(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=j(pg,true,p2in,true);
     jgbm=j(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.*p2i*MHmm)/(p2in-pg).m2();
     U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2();
     U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(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)/(q1.m2()*qH1.m2());
 
   // 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=qH2.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 jM2unobqHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
 
     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 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.*p2i*MHmm)/(p2in-pg).m2();
     U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2();
     U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(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)/(q1.m2()*qH1.m2());
 
   // 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=qH2.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 jM2unobqbarHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
 
     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 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.*p2i*MHmm)/(p2in-pg).m2();
     U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2();
     U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(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)/(q1.m2()*qH1.m2());
 
   // 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=qH2.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 jM2unobgHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
     // 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 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=j(p2out,true,p2in,true);
     mj2m=j(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=j(pg,true,p2in,true);
     jgbm=j(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.*p2i*MHmm)/(p2in-pg).m2();
     U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2();
     U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(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)/(q1.m2()*qH1.m2());
 
   // 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=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
     //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();
+    const double K = K_g(p1out, p1in);
 
-    double nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.;
-
-    return ampsq*nonflipcolourmult*9./4.; //ca/cf = 9/4
+    return ampsq*K/C_F;
 }
 
 
 
 double jM2unobgHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb)
 {
 
     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 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.*p2i*MHmm)/(p2in-pg).m2();
     U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2();
     U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2();
     U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(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)/(q1.m2()*qH1.m2());
 
   // 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=qH2.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.;
+    const double K = K_g(p1out, p1in);
 
-    return ampsq*nonflipcolourmult*9./4.; //ca/cf = 9/4
+    return ampsq*K/C_F; //ca/cf = 9/4
 
 }
 
 // Begin finite mass stuff
 #ifdef RHEJ_BUILD_WITH_QCDLOOP
 namespace {
 
 
   // All the stuff needed for the box functions in qg->qgH now...
 
   //COM E1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM E1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
     Sigma = 4.*s12*s34 - pow(S1+S2,2);
 
     return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) +
           S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) +
           2.*(s34 + S1)*(s34 + S1)/Delta +
           S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2,
             k1 + q2, mq) -
           s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
           S1*C0DD(k1, q2,
             mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
           2.*(s34 + S1)/Delta +
           S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) -
           C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) -
        C0DD(k1 + k2, q2, mq)*2.*s34/
          S1 - (B0DD(k1 + q2, mq) -
           B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 +
            S1)/(S1*Delta) + (B0DD(q2, mq) -
           B0DD(k1 + k2 + q2, mq) +
           s12*C0DD(k1 + k2, q2,
             mq))*(2.*s34*(s34 +
              S1)*(S1 - S2)/(Delta*Sigma) +
           2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) -
           B0DD(k1 + k2 + q2,
            mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
           S1)*(2.*s12*s34 -
            S2*(S1 + S2))/(Delta*Sigma));
   }
   //COM F1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM F1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
     Sigma = 4.*s12*s34 - pow(S1+S2,2);
 
     return looprwfactor*(-S2*D0DD(k1, k2, q2,
          mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) -
           s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1,
             k2 + q2, mq) -
           s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
           S2*C0DD(k2, q2,
             mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
           S2*pow((s34 + S2),2)/Delta/Delta) - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12) -
        C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) -
           B0DD(k1 + k2 + q2,
            mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD(
            q2, mq) - B0DD(k1 + k2 + q2, mq) +
           s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 +
           S2)*(S2 - S1)/(Delta*Sigma) + (B0DD(
            k1 + k2, mq) -
           B0DD(k1 + k2 + q2,
            mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
           S2)*(2.*s12*s34 -
            S2*(S1 + S2))/(Delta*Sigma));
   }
   //COM G1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM G1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
 
     return looprwfactor*(S2*D0DD(k1, q2, k2,
          mq)*(Delta/s12/s12 - 4.*mq*mq/s12) -
        S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) -
           S1*C0DD(k1, q2, mq))*(1./
            s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
        S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) -
           S2*C0DD(k2, q2, mq))*(1./
            s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
        C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) -
           C0DD(k1, q2, mq))*4.*mq*mq/
          s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./
          s12 + (B0DD(k1 + q2, mq) -
           B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) -
           B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1)));
   }
   //COM E4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM E4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
     Sigma = 4.*s12*s34 - pow(S1+S2,2);
 
     return looprwfactor* (-s12*D0DD(k2, k1, q2,
          mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) -
           s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2,
             k1 + q2, mq) -
           s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
           S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
            s12*pow((s34 + S1),2)/Delta/Delta) -
        C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) -
           B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta +
           2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD(
            q2, mq) - B0DD(k1 + k2 + q2, mq) +
           s12*C0DD(k1 + k2, q2,
             mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) +
               s12*(S1 -
                  S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) -
           B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*((2.*s12*(2.*s12*s34 - S1*(S1 + S2) +
               s34*(S2 - S1)))/(Delta*Sigma)));
   }
   //COM F4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM F4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
     Sigma = 4.*s12*s34 - pow(S1+S2,2);
 
     return looprwfactor* (-s12*D0DD(k1, k2, q2,
          mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) +
           s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1,
             k2 + q2, mq) -
           s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
           S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
            s12*pow((s34 + S2),2)/Delta/Delta) -
        C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) -
           B0DD(k1 + k2 + q2, mq))*2.*(s34 +
            S2)/Delta + (B0DD(q2, mq) -
           B0DD(k1 + k2 + q2, mq) +
           s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 -
            S1*(S1 + S2) +
            s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) -
           B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(2.*s12*(2.*s12*s34 - S2*(S1 + S2) +
              s34*(S1 - S2))/(Delta*Sigma)));
   }
   //COM G4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM G4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
 
     return looprwfactor* (-D0DD(k1, q2, k2,
           mq)*(Delta/s12 + (s12 + S1)/2. -
           4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
           S1*C0DD(k1, q2, mq))*(1./
            s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD(
             k1 + q2, k2, mq) -
           S2*C0DD(k2, q2, mq))*(1./
            s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD(
            k1 + k2 + q2, mq) -
           B0DD(k1 + q2, mq))*2./(s12 + S2));
   }
   //COM E10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM E10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
     Sigma = 4.*s12*s34 - pow(S1+S2,2);
 
     return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta +
            12.*mq*mq*S1*(s34 + S1)/Delta/Delta -
            4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) -
            s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
             S1*C0DD(k1, q2, mq))*(1./Delta +
            4.*mq*mq*S1/Delta/Delta -
            4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) +
         C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) -
            4.*(s12 -
               2.*mq*mq)*(2.*s12*s34 -
                S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) -
            B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) +
            8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) -
            B0DD(k1 + k2 + q2, mq) +
            s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 +
               S2)*(2.*s12*s34 -
                S1*(S1 + S2))/(Delta*Sigma*Sigma) -
            4.*s34*(4.*s12 + 3.*S1 +
                S2)/(Delta*Sigma) +
            8.*s12*s34*(s34*(s12 + S2) -
                S1*(s34 +
                   S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) -
            B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2,
              mq))*(12.*s12*(2.*s34 + S1 +
               S2)*(2.*s12*s34 -
                S1*(S1 + S2))/(Delta*Sigma*Sigma) +
            8.*s12*S1*(s34*(s12 + S2) -
                S1*(s34 +
                   S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
           S1*(S1 + S2))/(Delta*Sigma));
   }
 
   //COM F10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM F10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
     Sigma = 4.*s12*s34 - pow(S1+S2,2);
 
     return looprwfactor* (s12*D0DD(k1, k2, q2,
           mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta +
            12.*mq*mq*s34*(s12 + S1)/Delta/Delta -
            4.*s12*pow((s34 + S2),2)/Delta/Delta -
            4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
            s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
             S2*C0DD(k2, q2, mq))*(1./Delta +
            4.*mq*mq*S1/Delta/Delta -
            4.*s12*(s34 + S2)/Delta/Delta -
            4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) -
         C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) +
            4.*s12*s34*(S2 - S1)/(Delta*Sigma) +
            4.*(s12 -
               2.*mq*mq)*(2.*s12*s34 -
                S1*(S1 + S2))/(Delta*Sigma)) - (B0DD(
             k2 + q2, mq) -
            B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) +
            8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) -
            B0DD(k1 + k2 + q2, mq) +
            s12*C0DD(k1 + k2, q2,
              mq))*(-12*s34*(2*s12 + S1 +
               S2)*(2.*s12*s34 -
                S1*(S1 + S2))/(Delta*Sigma*Sigma) -
            4.*s12*s34*s34/(S2*Delta*Delta) +
            4.*s34*S1/(Delta*Sigma) -
            4.*s34*(s12*s34*(2.*s12 + S2) -
                S1*S1*(2.*s12 +
                   S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) -
            B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 +
               S2)*(2.*s12*s34 -
                S1*(S1 + S2))/(Delta*Sigma*Sigma) +
            8.*s12*(2.*s34 + S1)/(Delta*Sigma) -
            8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) +
                s12*(S2 -
                   S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
           S1*(S1 + S2))/(Delta*Sigma));
 
   }
 
   //COM G10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM G10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //CLHEP::HepLorentzVector q2=k3+k4;
     CLHEP::HepLorentzVector 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 = 2.*k3.dot(k4);
     s34 = q2.m2();
     Delta = s12*s34 - S1*S2;
 
     return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. +
           4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1,
             k2 + q2, mq) -
           S1*C0DD(k1, q2, mq))*(1./Delta +
           4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2,
             k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta +
           4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) -
           B0DD(k1 + q2, mq))*4.*(s34 +
            S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) -
           B0DD(k2 + q2, mq))*4.*s34/(Delta*S2));
   }
 
   //COM H1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM H1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //return E1(k1,k2,k3,k4,mq)+F1(k1,k2,k3,k4,mq)+G1(k1,k2,k3,k4,mq);
     return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq);
   }
   //COM H4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM H4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //return E4(k1,k2,k3,k4,mq)+F4(k1,k2,k3,k4,mq)+G4(k1,k2,k3,k4,mq);
     return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq);
   }
   //COM H10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM H10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //return E10(k1,k2,k3,k4,mq)+F10(k1,k2,k3,k4,mq)+G10(k1,k2,k3,k4,mq);
     return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq);
   }
   //COM H2(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM H2(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //return -1.*H1(k2,k1,k3,k4,mq);
     return -1.*H1(k2,k1,kh,mq);
   }
   //COM H5(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM H5(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //return -1.*H4(k2,k1,k3,k4,mq);
     return -1.*H4(k2,k1,kh,mq);
   }
   //COM H12(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq)
   COM H12(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq)
   {
     //return -1.*H10(k2,k1,k3,k4,mq);
     return -1.*H10(k2,k1,kh,mq);
   }
 
   // FL and FT functions
   COM FL(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq)
   {
     CLHEP::HepLorentzVector Q = q1 + q2;
     double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
     return -1./(2.*detQ2)*((2.-
            3.*q1.m2()*q2.dot(Q)/detQ2)*(B0DD(q1, mq) -
            B0DD(Q, mq)) + (2. -
            3.*q2.m2()*q1.dot(Q)/detQ2)*(B0DD(q2, mq) -
            B0DD(Q, mq)) - (4.*mq*mq + q1.m2() + q2.m2() +
            Q.m2() - 3.*q1.m2()*q2.m2()*Q.m2()/detQ2)*C0DD(
           q1, q2, mq) - 2.);
   }
   COM FT(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq)
   {
     CLHEP::HepLorentzVector 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);
   }
 
   CLHEP::HepLorentzVector ParityFlip(CLHEP::HepLorentzVector p)
   {
     CLHEP::HepLorentzVector flippedVector;
     flippedVector.setE(p.e());
     flippedVector.setX(-p.x());
     flippedVector.setY(-p.y());
     flippedVector.setZ(-p.z());
     return flippedVector;
   }
   // HC amp for qg->qgH with finite top
   void g_gH_HC(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH, double mq, current &retAns)
   {
     current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1,epsHapart2,conjepsH1part1,conjepsH1part2;
     current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10, ans;
     COM ang1a,sqa1,aH1,Fta,Ft1,h4,h5,h10,h12;
     double prop;
     //cout << mq << endl;
     //double F = 4.*mq*mq/v/16./M_PI/M_PI;
     double F = 4.*mq*mq/v;
     //cout << mq << endl;
     //cout << HVE << endl;
     // 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){
       ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/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());
     }
 
 
     prop = (pa-p1-pH).m2();
 
     //double Acheck = 1./(12.*M_PI*M_PI*v);
 
     cmult(-1./sqrt(2)/ang1a,cura1,conjeps1);
     cmult(1./sqrt(2)/sqa1,cura1,epsa);
 
     //if(mq<1000.){
     Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
     Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
     //}
     //else{
     //  Fta = -Acheck/(pa-pH).m2()/F;
     //  Ft1 = -Acheck/(p1+pH).m2()/F;
     //}
 
     //if(mq<1000.){
     h4 = H4(p1,-pa,pH,mq);
     h5 = H5(p1,-pa,pH,mq);
     h10 = H10(p1,-pa,pH,mq);
     h12 = H12(p1,-pa,pH,mq);
     //}
     //else{
     // h4 = COM(0.,2.)/16./M_PI/M_PI*Acheck/F;
     // h5 = COM(0.,-2.)/16./M_PI/M_PI*Acheck/F;
     // h10 = 0;
     // h12 = 0;
     //}
 
 
     cmult(Fta*pa.dot(pH),epsa,epsHapart1);
     cmult(-1.*Fta*cdot(pHcur,epsa),pacur,epsHapart2);
     cmult(Ft1*cdot(pHcur,conjeps1),p1cur,conjepsH1part1);
     cmult(-Ft1*p1.dot(pH),conjeps1,conjepsH1part2);
     cadd(epsHapart1,epsHapart2,epsHa);
     cadd(conjepsH1part1,conjepsH1part2,conjepsH1);
     aH1 = cdot(pHcur,cura1);
 
     if(gluonforward){
       cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1,conjepsH1,T1);
       cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/ang1a,epsHa,T2);
     }
     else{
       cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/sqa1,conjepsH1,T1);
       cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/ang1a,epsHa,T2);
     }
 
     cmult(sqrt(2.)/ang1a*aH1,epsHa,T3);
     cmult(sqrt(2.)/sqa1*aH1,conjepsH1,T4);
 
     cmult(-sqrt(2.)*Fta*pa.dot(p1)*aH1/sqa1,conjeps1,T5);
     cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/ang1a,epsa,T6);
 
     cmult(-aH1/sqrt(2.)/sqa1*h4*8.*COM(0.,1.)*M_PI*M_PI,conjeps1,T7);
     cmult(aH1/sqrt(2.)/ang1a*h5*8.*COM(0.,1.)*M_PI*M_PI,epsa,T8);
     cmult(aH1*aH1/2./ang1a/sqa1*h10*8.*COM(0.,1.)*M_PI*M_PI,pacur,T9);
     cmult(-aH1*aH1/2./ang1a/sqa1*h12*8.*COM(0.,1.)*M_PI*M_PI,p1cur,T10);
 
     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];
   }
 
   // HNC amp for qg->qgH with finite top
   void g_gH_HNC(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH, double mq, current &retAns)
   {
     double prop;
     //COM F = 4.*mq*mq/v/16./M_PI/M_PI;
     double F = 4.*mq*mq/v;
     COM ang1a,Fta,Ft1,phase,aH1,oneHa,h1,h2,h4,h5,h10,h12,sqa1,Falpha,Fbeta;
     current ans,conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur,p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1,conjepsH1part2,T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a,T11b,T12a,T12b,T13;
     // 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);
     j(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);
     aH1 = cdot(pHcur,cura1);
     oneHa = cdot(pHcur,cur1a);
 
     if(gluonforward){
       ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/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());
     }
 
     prop = (pa-p1-pH).m2();
 
     cmult(1./sqrt(2)/sqa1,cur1a,epsa);
     cmult(-1./sqrt(2)/sqa1,cura1,conjeps1);
     phase = cdot(conjeps1,epsa);
     Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
     Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
     Falpha = FT(p1-pa,pa-p1-pH,mq);
     Fbeta = FL(p1-pa,pa-p1-pH,mq);
 
     h1 = H1(p1,-pa,pH,mq);
     h2 = H2(p1,-pa,pH,mq);
     h4 = H4(p1,-pa,pH,mq);
     h5 = H5(p1,-pa,pH,mq);
     h10 = H10(p1,-pa,pH,mq);
     h12 = H12(p1,-pa,pH,mq);
 
     cmult(Fta*pa.dot(pH),epsa,epsHapart1);
     cmult(-1.*Fta*cdot(pHcur,epsa),pacur,epsHapart2);
     cmult(Ft1*cdot(pHcur,conjeps1),p1cur,conjepsH1part1);
     cmult(-Ft1*p1.dot(pH),conjeps1,conjepsH1part2);
     cadd(epsHapart1,epsHapart2,epsHa);
     cadd(conjepsH1part1,conjepsH1part2,conjepsH1);
 
     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);
 
     }
 
     COM boxdiagFact = 8.*COM(0.,1.)*M_PI*M_PI;
 
     cmult(aH1*sqrt(2.)/sqa1,epsHa,T3);
     cmult(oneHa*sqrt(2.)/sqa1,conjepsH1,T4);
     cmult(-2.*phase*Fta*pa.dot(pH),p1cur,T5a);
     cmult(2.*phase*Ft1*p1.dot(pH),pacur,T5b);
     cmult(-sqrt(2.)*Fta*p1.dot(pa)*oneHa/sqa1,conjeps1,T6);
     cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/sqa1,epsa,T7);
 
     cmult(-boxdiagFact*phase*h2,pacur,T8a);
     cmult(boxdiagFact*phase*h1,p1cur,T8b);
     cmult(boxdiagFact*aH1/sqrt(2.)/sqa1*h5,epsa,T9);
     cmult(-boxdiagFact*oneHa/sqrt(2.)/sqa1*h4,conjeps1,T10);
     cmult(boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h10,pacur,T11a);
     cmult(-boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h12,p1cur,T11b);
 
     cmult(-phase/(pa-p1).m2()*Falpha*(p1-pa).dot(pa-p1-pH),paplusp1cur,T12a);
     cmult(phase/(pa-p1).m2()*Falpha*(pa+p1).dot(pa-p1-pH),p1minuspacur,T12b);
     cmult(-phase*Fbeta*(pa-p1-pH).m2(),paplusp1cur,T13);
 
     //cout << phase << endl;
 
     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(CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double mq, bool includeBottom, double mq2)
 {
 
   current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip;
   current retAns,retAnsb;
   j(p2out,true,p2in,true,cur2bplus);
   j(p2out,false,p2in,false,cur2bminus);
   j(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip);
   j(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))*12./4./(3.*8.);
   //return (abs2(app1))*12./(3.*8.);
   // factor = 12 (colour sum) /4 (hel avg) /(3*8) (col avg)
 }
 #endif // RHEJ_BUILD_WITH_QCDLOOP
 
 double C2gHgm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH)
 {
   static double A=1./(3.*M_PI*v);
   // 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/p1p*conj(p1perp)*p3perp+p1p/p2p*p1perp*conj(p3perp));
   temp=temp*conj(temp);
   return temp.real();
 }
 
 double C2gHgp(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH)
 {
   static double A=1./(3.*M_PI*v);
   // 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 og the opposite
   s12=p1.invariantMass2(-p2);
   if (p2.pz()>0.) { // case considered in hep-ph/0301013
     php=pH.plus();
     phm=pH.minus();
     p1p=p1.plus();
   } else { // opposite case
     php=pH.minus();
     phm=pH.plus();
     p1p=p1.minus();
   }
   p1perp=p1.px()+COM(0,1)*p1.py();
   phperp=pH.px()+COM(0,1)*pH.py();
   p3perp=-(p1perp+phperp);
 
   COM temp=-COM(0,1)*A/(2.*s12)*(conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p)+s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm)-pow(conj(p3perp)+(1.+php/p1p)*conj(p1perp),2)/((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH)))));
   temp=temp*conj(temp);
   return temp.real();
 }
 
 double C2qHqm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH)
 {
   static double A=1./(3.*M_PI*v);
   // 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 og the opposite
   s12=p1.invariantMass2(-p2);
   if (p2.pz()>0.) { // case considered in hep-ph/0301013
     p2p=p2.plus();
     p1p=p1.plus();
   } else { // opposite case
     p2p=p2.minus();
     p1p=p1.minus();
   }
   p1perp=p1.px()+COM(0,1)*p1.py();
   phperp=pH.px()+COM(0,1)*pH.py();
   p3perp=-(p1perp+phperp);
 
   COM temp=A/(2.*s12)*(sqrt(p2p/p1p)*p3perp*conj(p1perp)+sqrt(p1p/p2p)*p1perp*conj(p3perp));
   temp=temp*conj(temp);
   return temp.real();
 }