diff --git a/src/jets.cc b/src/jets.cc
index 4d41515..2ec1465 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,842 +1,830 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019
  *  \copyright GPLv2 or later
  */
 #include "HEJ/jets.hh"
 
 #include "HEJ/Tensor.hh"
 #include "HEJ/Constants.hh"
 
 // Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
 // @TODO: this is not a current and should be moved somewhere else
 double K_g(double p1minus, double paminus) {
     return 1./2.*(p1minus/paminus + paminus/p1minus)*(HEJ::C_A - 1./HEJ::C_A) + 1./HEJ::C_A;
 }
 double K_g(
   HLV const & pout,
   HLV const & pin
   ) {
   if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
   return K_g(pout.minus(), pin.minus());
 }
 
 CCurrent CCurrent::operator+(const CCurrent& other)
 {
     COM result_c0=c0 + other.c0;
     COM result_c1=c1 + other.c1;
     COM result_c2=c2 + other.c2;
     COM result_c3=c3 + other.c3;
 
     return CCurrent(result_c0,result_c1,result_c2,result_c3);
 }
 
 CCurrent CCurrent::operator-(const CCurrent& other)
 {
     COM result_c0=c0 - other.c0;
     COM result_c1=c1 - other.c1;
     COM result_c2=c2 - other.c2;
     COM result_c3=c3 - other.c3;
 
     return CCurrent(result_c0,result_c1,result_c2,result_c3);
 }
 
 CCurrent CCurrent::operator*(const double x)
 {
     COM result_c0=x*CCurrent::c0;
     COM result_c1=x*CCurrent::c1;
     COM result_c2=x*CCurrent::c2;
     COM result_c3=x*CCurrent::c3;
 
     return CCurrent(result_c0,result_c1,result_c2,result_c3);
 }
 
 CCurrent CCurrent::operator/(const double x)
 {
     COM result_c0=CCurrent::c0/x;
     COM result_c1=CCurrent::c1/x;
     COM result_c2=CCurrent::c2/x;
     COM result_c3=CCurrent::c3/x;
 
     return CCurrent(result_c0,result_c1,result_c2,result_c3);
 }
 
 CCurrent CCurrent::operator*(const COM x)
 {
     COM result_c0=x*CCurrent::c0;
     COM result_c1=x*CCurrent::c1;
     COM result_c2=x*CCurrent::c2;
     COM result_c3=x*CCurrent::c3;
 
     return CCurrent(result_c0,result_c1,result_c2,result_c3);
 }
 
 CCurrent CCurrent::operator/(const COM x)
 {
     COM result_c0=(CCurrent::c0)/x;
     COM result_c1=(CCurrent::c1)/x;
     COM result_c2=(CCurrent::c2)/x;
     COM result_c3=(CCurrent::c3)/x;
 
     return CCurrent(result_c0,result_c1,result_c2,result_c3);
 }
 
 std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
 {
     os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
     return os;
 }
 
 CCurrent operator * ( double x, CCurrent& m)
 {
     return m*x;
 }
 
 CCurrent operator * ( COM x, CCurrent& m)
 {
     return m*x;
 }
 
 CCurrent operator / ( double x, CCurrent& m)
 {
     return m/x;
 }
 
 CCurrent operator / ( COM x, CCurrent& m)
 {
     return m/x;
 }
 
 COM CCurrent::dot(HLV p1)
 {
     //  Current goes (E,px,py,pz)
     //  Vector goes (px,py,pz,E)
     return p1[3]*c0-p1[0]*c1-p1[1]*c2-p1[2]*c3;
 }
 
 COM CCurrent::dot(CCurrent p1)
 {
     return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3;
 }
 
 //Current Functions
 void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
   cur[0]=0.;
   cur[1]=0.;
   cur[2]=0.;
   cur[3]=0.;
 
   const double sqpop = sqrt(std::abs(pout.plus()));
   const double sqpom = sqrt(std::abs(pout.minus()));
 
   // Allow for "jii" format
   const COM poperp = (pout.x()==0 && pout.y() ==0) ? -1 : pout.x()+COM(0,1)*pout.y();
 
   if (helout != helin) {
     throw std::invalid_argument{"Non-matching helicities"};
   } else if (helout == false) { // negative helicity
     if (pin.plus() > pin.minus()) { // if forward
       const double sqpip = sqrt(std::abs(pin.plus()));
       cur[0] = sqpop * sqpip;
       cur[1] = sqpom * sqpip * poperp / std::abs(poperp);
       cur[2] = -COM(0,1) * cur[1];
       cur[3] = cur[0];
     } else { // if backward
       const double sqpim = sqrt(std::abs(pin.minus()));
       cur[0] = -sqpom * sqpim * poperp / std::abs(poperp);
       cur[1] = -sqpim * sqpop;
       cur[2] = COM(0,1) * cur[1];
       cur[3] = -cur[0];
     }
   } else { // positive helicity
     if (pin.plus() > pin.minus()) { // if forward
       const double sqpip = sqrt(std::abs(pin.plus()));
       cur[0] = sqpop * sqpip;
       cur[1] = sqpom * sqpip * conj(poperp) / std::abs(poperp);
       cur[2] = COM(0,1) * cur[1];
       cur[3] = cur[0];
     } else { // if backward
       double sqpim = sqrt(std::abs(pin.minus()));
       cur[0] = -sqpom * sqpim * conj(poperp) / std::abs(poperp);
       cur[1] = -sqpim * sqpop;
       cur[2] = -COM(0,1) * cur[1];
       cur[3] = -cur[0];
     }
   }
 }
 
 CCurrent joi (HLV pout, bool helout, HLV pin, bool helin)
 {
   current cur;
   joi(pout, helout, pin, helin, cur);
   return CCurrent(cur[0],cur[1],cur[2],cur[3]);
 }
 
 void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) {
   joi(pout, !helout, pin, !helin, cur);
 }
 
 CCurrent jio (HLV pin, bool helin, HLV pout, bool helout)
 {
   current cur;
   jio(pin, helin, pout, helout, cur);
   return CCurrent(cur[0],cur[1],cur[2],cur[3]);
 }
 
 void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {
 
   // Zero our current
   cur[0] = 0.0;
   cur[1] = 0.0;
   cur[2] = 0.0;
   cur[3] = 0.0;
   if (heli!=helj) {
     throw std::invalid_argument{"Non-matching helicities"};
   } else if ( heli == true ) { // If positive helicity swap momenta
     std::swap(pi,pj);
   }
 
   const double sqpjp = sqrt(std::abs(pj.plus() ));
   const double sqpjm = sqrt(std::abs(pj.minus()));
   const double sqpip = sqrt(std::abs(pi.plus() ));
   const double sqpim = sqrt(std::abs(pi.minus()));
   // Allow for "jii" format
   const COM piperp = (pi.x()==0 && pi.y() ==0) ? -1 : pi.x()+COM(0,1)*pi.y();
   const COM pjperp = (pj.x()==0 && pj.y() ==0) ? -1 : pj.x()+COM(0,1)*pj.y();
 
   const COM phasei = piperp / std::abs(piperp);
   const COM phasej = pjperp / std::abs(pjperp);
 
   cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
   cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
   cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
   cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
 }
 
 CCurrent joo (HLV pi, bool heli, HLV pj, bool helj)
 {
   current cur;
   joo(pi, heli, pj, helj, cur);
   return CCurrent(cur[0],cur[1],cur[2],cur[3]);
 }
 namespace{
 //@{
   /**
    * @brief Pure Jet FKL Contributions, function to handle all incoming types.
    * @param p1out             Outgoing Particle 1.
    * @param p1in              Incoming Particle 1.
    * @param p2out             Outgoing Particle 2
    * @param p2in              Incoming Particle 2
    *
    * Calculates j_\mu    j^\mu.
    * Handles all possible incoming states. Helicity doesn't matter since we sum
    * over all of them.
    */
   double j_j(HLV const & p1out, HLV const & p1in,
              HLV const & p2out, HLV const & p2in
   ){
     HLV const q1=p1in-p1out;
     HLV const q2=-(p2in-p2out);
     current mj1m,mj1p,mj2m,mj2p;
 
     // Note need to flip helicities in anti-quark case.
     joi(p1out, false, p1in, false, mj1p);
     joi(p1out, true,  p1in, true,  mj1m);
     joi(p2out, false, p2in, false, mj2p);
     joi(p2out, true,  p2in, true,  mj2m);
 
     COM const Mmp=cdot(mj1m,mj2p);
     COM const Mmm=cdot(mj1m,mj2m);
     COM const Mpp=cdot(mj1p,mj2p);
     COM const Mpm=cdot(mj1p,mj2m);
 
     double const sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);
 
     // Multiply by Cf^2
     return HEJ::C_F*HEJ::C_F*(sst)/(q1.m2()*q2.m2());
   }
 } //anonymous namespace
 double ME_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return j_j(p1out, p1in, p2out, p2in);
 }
 
 double ME_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return j_j(p1out, p1in, p2out, p2in);
 }
 
 double ME_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return j_j(p1out, p1in, p2out, p2in);
 }
 
 double ME_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/HEJ::C_F;
 }
 
 double ME_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/(HEJ::C_F);
 }
 
 double ME_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return j_j(p1out, p1in, p2out, p2in)*K_g(p1out, p1in)*K_g(p2out, p2in)/(HEJ::C_F*HEJ::C_F);
 }
 //@}
 
 namespace{
   double juno_j(HLV const & pg, HLV const & p1out,
                 HLV const & p1in, HLV const & p2out, HLV const & p2in
   ){
     //  This construction is taking rapidity order: pg > p1out >> p2out
     HLV q1=p1in-p1out;  // Top End
     HLV q2=-(p2in-p2out);   // Bottom End
     HLV qg=p1in-p1out-pg;  // Extra bit post-gluon
 
     // Note <p1|eps|pa> current split into two by gauge choice.
     // See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
     CCurrent mj1p=joi(p1out, false, p1in, false);
     CCurrent mj1m=joi(p1out,  true, p1in,  true);
     CCurrent jgap=joi(pg,    false, p1in, false);
     CCurrent jgam=joi(pg,     true, p1in,  true);
 
     // Note for function joo(): <p1+|pg+> = <pg-|p1->.
     CCurrent j2gp=joo(p1out, false, pg, false);
     CCurrent j2gm=joo(p1out,  true, pg,  true);
 
     CCurrent mj2p=joi(p2out, false, p2in, false);
     CCurrent mj2m=joi(p2out,  true, p2in,  true);
 
     // Dot products of these which occur again and again
     COM Mmp=mj1m.dot(mj2p);
     COM Mmm=mj1m.dot(mj2m);
     COM Mpp=mj1p.dot(mj2p);
     COM Mpm=mj1p.dot(mj2m);
 
     CCurrent p1o(p1out),p2o(p2out),p2i(p2in),qsum(q1+qg),p1i(p1in);
 
     CCurrent Lmm=(qsum*(Mmm)+(-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m
                   +(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
     CCurrent Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p
                   +(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
     CCurrent Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m
                   +(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
     CCurrent Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p
                   +(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
 
     CCurrent U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
     CCurrent U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
     CCurrent U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
     CCurrent U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
     CCurrent U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
     CCurrent U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
     CCurrent U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
     CCurrent U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
 
     constexpr double cf=HEJ::C_F;
 
     double amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
     double amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
     double apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
     double app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
     double ampsq=-(amm+amp+apm+app);
 
     //Divide by t-channels
     ampsq/=q2.m2()*qg.m2();
     ampsq/=16.;
 
     // Factor of (Cf/Ca) for each quark to match j_j.
     ampsq*=(HEJ::C_F*HEJ::C_F)/(HEJ::C_A*HEJ::C_A);
 
     return ampsq;
 
   }
 }
 
 //Unordered bits for pure jet
 double ME_unob_qQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return juno_j(pg, p1out, p1in, p2out, p2in);
 }
 
 double ME_unob_qbarQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return juno_j(pg, p1out, p1in, p2out, p2in);
 }
 
 double ME_unob_qQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return juno_j(pg, p1out, p1in, p2out, p2in);
 }
 
 double ME_unob_qbarQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return juno_j(pg, p1out, p1in, p2out, p2in);
 }
 
 double ME_unob_qg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/HEJ::C_F;
 }
 
 double ME_unob_qbarg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
   return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/HEJ::C_F;
 }
 
 namespace {
   void eps(HLV refmom, HLV kb, bool hel, current &ep){
     //Positive helicity eps has negative helicity choices for spinors and vice versa
     joi(refmom,hel,kb,hel,ep);
     double norm;
     if(kb.z()<0.) norm = sqrt(2.*refmom.plus()*kb.minus());
     else          norm = sqrt(2.*refmom.minus()*kb.plus());
     // Normalise
     std::for_each(begin(ep), end(ep), [&,norm](auto & val){val/=norm;});
   }
 
   COM qggm1(HLV pa, HLV pb, HLV p1, HLV p2, HLV p3, bool helchain,
             bool heltop, bool helb,HLV refmom){
     // Since everything is defined with currents, need to use compeleness relation
     // to expand p slash. i.e. pslash = |p><p|. Only one helicity 'survives' as
     // defined by the helicities of the spinors at the end of the chain.
     current cur33, cur23, curb3, cur2b, cur1a, ep;
     joo(p3, helchain, p3, helchain,cur33);
     joo(p2,helchain,p3,helchain,cur23);
     jio(pb,helchain,p3,helchain,curb3);
     joi(p2,helchain,pb,helchain,cur2b);
     joi(p1, heltop, pa, heltop,cur1a);
 
     const double t2 = (p3-pb)*(p3-pb);
     //Calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
     COM v1[4][4];
     for(int u=0; u<4;++u){
       for(int v=0; v<4;++v){
         v1[u][v]=(cur23[u]*cur33[v]-cur2b[u]*curb3[v])/t2*(-1.);
       }
     }
     //Dot in current and eps
     //Metric tensor
     auto eta = HEJ::metric();
     //eps
     eps(refmom,pb,helb, ep);
     COM M1=0.;
     // Perform Contraction: g^{ik} j_{1a}_k * v1_i^j eps^l g_lj
     for(int i=0;i<4;++i){
       for(int j=0;j<4;++j){
         M1+= eta(i,i) *cur1a[i]*(v1[i][j])*ep[j]*eta(j,j);
       }
     }
     return M1;
   }
 
   COM qggm2(HLV pa, HLV pb, HLV p1, HLV p2, HLV p3, bool helchain, bool heltop,
             bool helb,HLV refmom){
     // Since everything is defined with currents, need to use completeness relation
     // to expand p slash. i.e. pslash = |p><p|. Only one helicity 'survives' as
     // defined by the helicities of the spinors at the end of the chain.
     current cur22, cur23, curb3, cur2b, cur1a, ep;
     joo(p2, helchain, p2, helchain, cur22);
     joo(p2, helchain, p3, helchain, cur23);
     jio(pb, helchain, p3, helchain, curb3);
     joi(p2, helchain, pb, helchain, cur2b);
     joi(p1, heltop,   pa, heltop,   cur1a);
 
     const double t2t = (p2-pb)*(p2-pb);
     //Calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
     COM v2[4][4]={};
     for(int u=0; u<4;++u){
       for(int v=0; v<4; ++v){
         v2[u][v]=(cur22[v]*cur23[u]-cur2b[v]*curb3[u])/t2t;
       }
     }
     //Dot in current and eps
     auto eta=HEJ::metric();
     //eps
     eps(refmom,pb,helb, ep);
     COM M2=0.;
     // Perform Contraction: g^{ik} j_{1a}_k * v2_i^j eps^l g_lj
     for(int i=0;i<4;++i){
       for(int j=0;j<4;++j){
         M2+= eta(i,i)*cur1a[i]*(v2[i][j])*ep[j]*eta(j,j);
       }
     }
     return M2;
   }
 
   COM qggm3(HLV pa, HLV pb, HLV p1, HLV p2, HLV p3, bool helchain, bool heltop,
             bool helb,HLV refmom){
     current qqcur,ep,cur1a;
     const double s23 = (p2+p3)*(p2+p3);
     joo(p2,helchain,p3,helchain,qqcur);
     joi(p1, heltop, pa, heltop,cur1a);
     //Redefine relevant momenta as currents - for ease of calling correct part of vector
     const current kb{pb.e(), pb.x(), pb.y(), pb.z()};
     const current k2{p2.e(), p2.x(), p2.y(), p2.z()};
     const current k3{p3.e(), p3.x(), p3.y(), p3.z()};
     auto eta=HEJ::metric();
     //Calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
     COM V3g[4][4]={};
     const COM kbqq=kb[0]*qqcur[0] -kb[1]*qqcur[1] -kb[2]*qqcur[2] -kb[3]*qqcur[3];
     for(int u=0;u<4;++u){
       for(int v=0;v<4;++v){
         V3g[u][v] += 2.*COM(0.,1.)*(((k2[v]+k3[v])*qqcur[u] - (kb[u])*qqcur[v])+
                                     kbqq*eta(u,v))/s23;
       }
     }
     eps(refmom,pb,helb, ep);
     COM M3=0.;
     // Perform Contraction: g^{ik} j_{1a}_k * (v2_i^j) eps^l g_lj
     for(int i=0;i<4;++i){
       for(int j=0;j<4;++j){
         M3+= eta(i,i)*cur1a[i]*(V3g[i][j])*ep[j]*eta(j,j);
       }
     }
     return M3;
   }
 
   //Now the function to give helicity/colour sum/average
   double MqgtqQQ(HLV pa, HLV pb, HLV p1, HLV p2, HLV p3){
     // 4 indepedent helicity choices (complex conjugation related).
     //Need to evalute each independent hel configuration and store that result somewhere
     const COM Mmmm1 = qggm1(pa,pb,p1,p2,p3,false,false,false, pa);
     const COM Mmmm2 = qggm2(pa,pb,p1,p2,p3,false,false,false, pa);
     const COM Mmmm3 = qggm3(pa,pb,p1,p2,p3,false,false,false, pa);
     const COM Mmmp1 = qggm1(pa,pb,p1,p2,p3,false,true, false, pa);
     const COM Mmmp2 = qggm2(pa,pb,p1,p2,p3,false,true, false, pa);
     const COM Mmmp3 = qggm3(pa,pb,p1,p2,p3,false,true, false, pa);
     const COM Mpmm1 = qggm1(pa,pb,p1,p2,p3,false,false,true,  pa);
     const COM Mpmm2 = qggm2(pa,pb,p1,p2,p3,false,false,true,  pa);
     const COM Mpmm3 = qggm3(pa,pb,p1,p2,p3,false,false,true,  pa);
     const COM Mpmp1 = qggm1(pa,pb,p1,p2,p3,false,true, true,  pa);
     const COM Mpmp2 = qggm2(pa,pb,p1,p2,p3,false,true, true,  pa);
     const COM Mpmp3 = qggm3(pa,pb,p1,p2,p3,false,true, true,  pa);
 
     //Colour factors:
     const COM cm1m1 = 8./3.;
     const COM cm2m2 = 8./3.;
     const COM cm3m3 = 6.;
     const COM cm1m2 = -1./3.;
     const COM cm1m3 = -3.*COM(0.,1.);
     const COM cm2m3 = 3.*COM(0.,1.);
 
     //Sqaure and sum for each helicity config:
     const double Mmmm = real(cm1m1*pow(abs(Mmmm1),2)+cm2m2*pow(abs(Mmmm2),2)+
                              cm3m3*pow(abs(Mmmm3),2)+2.*real(cm1m2*Mmmm1*conj(Mmmm2))+
                              2.*real(cm1m3*Mmmm1*conj(Mmmm3))+2.*real(cm2m3*Mmmm2*conj(Mmmm3)));
     const double Mmmp = real(cm1m1*pow(abs(Mmmp1),2)+cm2m2*pow(abs(Mmmp2),2)+
                              cm3m3*pow(abs(Mmmp3),2)+2.*real(cm1m2*Mmmp1*conj(Mmmp2))+
                              2.*real(cm1m3*Mmmp1*conj(Mmmp3))+2.*real(cm2m3*Mmmp2*conj(Mmmp3)));
     const double Mpmm = real(cm1m1*pow(abs(Mpmm1),2)+cm2m2*pow(abs(Mpmm2),2)+
                              cm3m3*pow(abs(Mpmm3),2)+2.*real(cm1m2*Mpmm1*conj(Mpmm2))+
                              2.*real(cm1m3*Mpmm1*conj(Mpmm3))+2.*real(cm2m3*Mpmm2*conj(Mpmm3)));
     const double Mpmp = real(cm1m1*pow(abs(Mpmp1),2)+cm2m2*pow(abs(Mpmp2),2)+
                              cm3m3*pow(abs(Mpmp3),2)+2.*real(cm1m2*Mpmp1*conj(Mpmp2))+
                              2.*real(cm1m3*Mpmp1*conj(Mpmp3))+2.*real(cm2m3*Mpmp2*conj(Mpmp3)));
 
     // Factor of 2 for the helicity for conjugate configurations
     return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/3.)/(pa-p1).m2()/(p2+p3-pb).m2();
   }
 }
 
 // Extremal qqx
 double ME_Exqqx_qbarqQ(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
   return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout);
 }
 
 double ME_Exqqx_qqbarQ(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
   return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout);
 }
 
 double ME_Exqqx_qbarqg(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
   return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout)*K_g(p2out,p2in)/HEJ::C_F;
 }
 
 double ME_Exqqx_qqbarg(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
   return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout)*K_g(p2out,p2in)/HEJ::C_F;
 }
 
 namespace {
 
   void CurrentMatrix(current j1, current j2, COM array[4][4]){
     for(int i=0;i<4;i++){
       for(int j=0;j<4;j++){
         array[i][j]=j1[i]*j2[j];
       }
     }
   }
 
 
 //qqbar produced in the middle
   COM m1(current jtop, current jbot, bool hel2, HLV ka, HLV kb, HLV k2, HLV k3,
          std::vector<HLV> partons, unsigned int nabove){
     //Define a load of invaraints I need
     const double s23 = 2.*(k2*k3);
     const double sa2 = 2.*(ka*k2);
     const double sa3 = 2.*(ka*k3);
     const double s12 = 2.*(partons.front()*k2);
     const double s13 = 2.*(partons.front()*k3);
     const double sb2 = 2.*(kb*k2);
     const double sb3 = 2.*(kb*k3);
     const double s42 = 2.*(partons.back()*k2);
     const double s43 = 2.*(partons.back()*k3);
     HLV q1=ka-partons.front();
     for(unsigned int i=1;i<nabove+1;i++)
       q1-=partons.at(i);
     const HLV q2=q1-partons.at(nabove+2)-partons.at(nabove+1);
     const double t1 = q1.m2();
     const double t3 = q2.m2();
     //Easier to have everything in terms of 'currents'
     //(E,px,py,pz). To make the distinction between actual currents of
     //the form ubar gamma u and 4-vectors being placed under the
     //'current' class (to make dot products work out), all actual
     //currents will have either a j or 'cur' in their variable name.
     current cur23,curka,curkb,curk1,curk2,curk3,curk4,qc1,qc2;
     //From what I gather, joo is what I use for two outgoing momenta,
     //jio with one outgoing and one incoming (in that order) and j for
     //lines when pa,pb are on the right of the spinor product
     joo(k2,hel2,k3,hel2,cur23);
     curka[0]=ka.e();
     curka[1]=ka.px();
     curka[2]=ka.py();
     curka[3]=ka.pz();
     curkb[0]=kb.e();
     curkb[1]=kb.px();
     curkb[2]=kb.py();
     curkb[3]=kb.pz();
     curk1[0]=partons.front().e();
     curk1[1]=partons.front().px();
     curk1[2]=partons.front().py();
     curk1[3]=partons.front().pz();
     curk2[0]=k2.e();
     curk2[1]=k2.px();
     curk2[2]=k2.py();
     curk2[3]=k2.pz();
     curk3[0]=k3.e();
     curk3[1]=k3.px();
     curk3[2]=k3.py();
     curk3[3]=k3.pz();
     curk4[0]=partons.back().e();
     curk4[1]=partons.back().px();
     curk4[2]=partons.back().py();
     curk4[3]=partons.back().pz();
     qc1[0]=q1.e();
     qc1[1]=q1.px();
     qc1[2]=q1.py();
     qc1[3]=q1.pz();
     qc2[0]=q2.e();
     qc2[1]=q2.px();
     qc2[2]=q2.py();
     qc2[3]=q2.pz();
 
 
     //Metric tensor
     auto eta=HEJ::metric();
     //Create the two bits of this vertex
     COM veik[4][4],v3g[4][4];
     for(int i=0;i<4;i++) {
       for(int j=0;j<4;j++){
         veik[i][j] = (cdot(cur23,curka)*(t1/(sa2+sa3))+cdot(cur23,curk1)*
                       (t1/(s12+s13))-cdot(cur23,curkb)*(t3/(sb2+sb3))-
                       cdot(cur23,curk4)*(t3/(s42+s43)))*eta(i,j);
       }
     }
     for(int i=0;i<4;i++){
       for(int j=0;j<4;j++){
         v3g[i][j] = qc1[j]*cur23[i]+curk2[j]*cur23[i]+curk3[j]*cur23[i]+
           qc2[i]*cur23[j]-curk2[i]*cur23[j]-curk3[i]*cur23[j]-
           (cdot(qc1,cur23)+cdot(qc2,cur23))*eta(i,j);
       }
     }
     //Now dot in the currents - potential problem here with Lorentz
     //indicies, so check this
     COM M1=0;
     for(int i=0;i<4;i++){
       for(int j=0;j<4;j++){
-        for(int k=0; k<4; k++){
-          for(int l=0; l<4;l++){
-            M1+= eta(i,k)*jtop[k]*(veik[i][j]+v3g[i][j])*jbot[l]*eta(l,j);
-          }
-        }
+        M1+= eta(i,i)*jtop[i]*(veik[i][j]+v3g[i][j])*jbot[j]*eta(j,j);
       }
     }
     M1/=s23;
     return M1;
   }
 
   COM m2 (current jtop, current jbot, bool hel2, HLV ka, HLV k2,
           HLV k3, std::vector<HLV> partons, unsigned int nabove){
     //In order to get correct momentum dependence in the vertex, forst
     //have to work with CCurrent objects and then convert to 'current'
     current cur22,cur23,cur2q,curq3;
     COM qarray[4][4]={};
     COM temp[4][4]={};
     joo(k2,hel2,k2,hel2,cur22);
     joo(k2,hel2,k3,hel2,cur23);
     joi(k2,hel2,ka,hel2,cur2q);
     jio(ka,hel2,k3,hel2,curq3);
     CurrentMatrix(cur2q, curq3, qarray);
     for(unsigned int i =0; i<nabove+1; i++){
       joo(k2,hel2,partons.at(i),hel2,cur2q);
       joo(partons.at(i),hel2,k3,hel2,curq3);
       CurrentMatrix(cur2q, curq3, temp);
       for(int ii=0;ii<4;ii++){
         for(int jj=0;jj<4;jj++){
           qarray[ii][jj]=qarray[ii][jj]-temp[ii][jj];
         }
       }
     }
     HLV qt=ka-k2;
     for(unsigned int i=0; i<nabove+1;i++){
       qt-=partons.at(i);
     }
     const double t2=qt*qt;
     //Metric tensor
     auto eta=HEJ::metric();
     COM tempv[4][4];
     for(int i=0; i<4;i++){
       for(int j=0;j<4;j++){
         tempv[i][j] = COM(0.,1.)*(qarray[i][j]-cur22[i]*cur23[j]);
       }
     }
     COM M2=0.;
     for(int i=0;i<4;i++){
       for(int j=0;j<4;j++){
-        for(int k=0; k<4; k++){
-          for(int l=0; l<4;l++){
-            M2+= eta(i,k)*jtop[k]*(tempv[i][j])*jbot[l]*eta(l,j);
-          }
-        }
+        M2+= eta(i,i)*jtop[i]*(tempv[i][j])*jbot[j]*eta(j,j);
       }
     }
     M2/=t2;
     return M2;
   }
 
   COM m3 (current jtop, current jbot, bool hel2, HLV ka, HLV k2,
           HLV k3, std::vector<HLV> partons, unsigned int nabove){
     COM M3=0.;
 
     current cur23,cur33,cur2q,curq3;
     COM qarray[4][4]={};
     COM temp[4][4]={};
     joo(k3,hel2,k3,hel2,cur33);
     joo(k2,hel2,k3,hel2,cur23);
     joi(k2,hel2,ka,hel2,cur2q);
     jio(ka,hel2,k3,hel2,curq3);
     CurrentMatrix(cur2q, curq3, qarray);
     for(unsigned int i =0; i<nabove+1; i++){
       joo(k2,hel2,partons.at(i),hel2,cur2q);
       joo(partons.at(i),hel2,k3,hel2,curq3);
       CurrentMatrix(cur2q, curq3, temp);
       for(int ii=0;ii<4;ii++){
         for(int jj=0;jj<4;jj++){
           qarray[ii][jj]=qarray[ii][jj]-temp[ii][jj];
         }
       }
     }
     HLV qt=ka-k3;
     for(unsigned int i=0; i<nabove+1;i++){
       qt-=partons.at(i);
     }
     const double t2t=qt*qt;
 
     //Metric tensor
     auto eta=HEJ::metric();
     COM tempv[4][4];
 
     for(int i=0; i<4;i++){
       for(int j=0;j<4;j++){
         tempv[i][j] = COM(0.,-1.)*(qarray[j][i]-cur23[j]*cur33[i]);
       }
     }
 
     for(int i=0;i<4;i++){
       for(int j=0;j<4;j++){
-        for(int k=0; k<4; k++){
-          for(int l=0; l<4;l++){
-            M3+= eta(i,k)*jtop[k]*(tempv[i][j])*jbot[l]*eta(l,j);
-          }
-        }
+        M3+= eta(i,i)*jtop[i]*(tempv[i][j])*jbot[j]*eta(j,j);
       }
     }
     M3/= t2t;
     return M3;
   }
 }
 
 double ME_Cenqqx_qq(HLV ka, HLV kb, std::vector<HLV> partons, bool aqlinepa,
                     bool aqlinepb, bool qqxmarker, int nabove){
   //Partons in the 'wrong' ordering, so reverse it. Just put ka =
   //forward and kb = backward into the function call
   std::reverse(partons.begin(),partons.end());
 
   //Get all the possible outer currents
   current j1p,j1m,j4p,j4m;
   if(!(aqlinepa)){
     joi(partons.front(),true,ka,true,j1p);
     joi(partons.front(),false,ka,false,j1m);
   }
   if(aqlinepa){
     jio(ka,true,partons.front(),true,j1p);
     jio(ka,false,partons.front(),false,j1m);
   }
   if(!(aqlinepb)){
     joi(partons.back(),true,kb,true,j4p);
     joi(partons.back(),false,kb,false,j4m);
   }
   if(aqlinepb){
     jio(kb,true,partons.back(),true,j4p);
     jio(kb,false,partons.back(),false,j4m);
   }
 
   HLV k2,k3;
   if(!(qqxmarker)){
     k2=partons.at(nabove+1);
     k3=partons.at(nabove+2);
   }
   else{
     k2=partons.at(nabove+2);
     k3=partons.at(nabove+1);
   }
 
   //8 helicity choices we can make, but only 4 indepedent ones
   //(complex conjugation related).
 
   const COM Mmmm1 = m1(j1m,j4m,false,ka,kb,k2,k3,partons,nabove);
   const COM Mmmm2 = m2(j1m,j4m,false,ka,   k2,k3,partons,nabove);
   const COM Mmmm3 = m3(j1m,j4m,false,ka,   k2,k3,partons,nabove);
   const COM Mmmp1 = m1(j1m,j4m,true, ka,kb,k2,k3,partons,nabove);
   const COM Mmmp2 = m2(j1m,j4m,true, ka,   k2,k3,partons,nabove);
   const COM Mmmp3 = m3(j1m,j4m,true, ka,   k2,k3,partons,nabove);
   const COM Mpmm1 = m1(j1p,j4m,false,ka,kb,k2,k3,partons,nabove);
   const COM Mpmm2 = m2(j1p,j4m,false,ka,   k2,k3,partons,nabove);
   const COM Mpmm3 = m3(j1p,j4m,false,ka,   k2,k3,partons,nabove);
   const COM Mpmp1 = m1(j1p,j4m,true, ka,kb,k2,k3,partons,nabove);
   const COM Mpmp2 = m2(j1p,j4m,true, ka,   k2,k3,partons,nabove);
   const COM Mpmp3 = m3(j1p,j4m,true, ka,   k2,k3,partons,nabove);
 
   //Colour factors:
   const COM cm1m1=3.;
   const COM cm2m2=4./3.;
   const COM cm3m3=4./3.;
   const COM cm1m2 =3./2.*COM(0.,1.);
   const COM cm1m3 = -3./2.*COM(0.,1.);
   const COM cm2m3 = -1./6.;
 
   //Square and sum for each helicity config:
   const double Mmmm = real(cm1m1*pow(abs(Mmmm1),2)+cm2m2*pow(abs(Mmmm2),2)+
               cm3m3*pow(abs(Mmmm3),2)+2.*real(cm1m2*Mmmm1*conj(Mmmm2))+
               2.*real(cm1m3*Mmmm1*conj(Mmmm3))+2.*real(cm2m3*Mmmm2*conj(Mmmm3)));
   const double Mmmp = real(cm1m1*pow(abs(Mmmp1),2)+cm2m2*pow(abs(Mmmp2),2)+
               cm3m3*pow(abs(Mmmp3),2)+2.*real(cm1m2*Mmmp1*conj(Mmmp2))+
               2.*real(cm1m3*Mmmp1*conj(Mmmp3))+2.*real(cm2m3*Mmmp2*conj(Mmmp3)));
   const double Mpmm = real(cm1m1*pow(abs(Mpmm1),2)+cm2m2*pow(abs(Mpmm2),2)+
               cm3m3*pow(abs(Mpmm3),2)+2.*real(cm1m2*Mpmm1*conj(Mpmm2))+
               2.*real(cm1m3*Mpmm1*conj(Mpmm3))+2.*real(cm2m3*Mpmm2*conj(Mpmm3)));
   const double Mpmp = real(cm1m1*pow(abs(Mpmp1),2)+cm2m2*pow(abs(Mpmp2),2)+
               cm3m3*pow(abs(Mpmp3),2)+2.*real(cm1m2*Mpmp1*conj(Mpmp2))+
               2.*real(cm1m3*Mpmp1*conj(Mpmp3))+2.*real(cm2m3*Mpmp2*conj(Mpmp3)));
 
   //Result (averaged, without coupling or t-channel props). Factor of
   //2 for the 4 helicity configurations I didn't work out explicitly
   HLV prop1 = ka;
   for(int i=0; i<=nabove; i++){
     prop1 -= partons[i];
   }
   const HLV prop2 = prop1 - k2 - k3;
 
   return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/9./4.) /
     ((ka-partons.front()).m2()*(kb-partons.back()).m2()*prop1.m2()*prop2.m2());
 }