diff --git a/src/Wjets.cc b/src/Wjets.cc
index a59afab..3e92a2e 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,2073 +1,2067 @@
 #include "HEJ/currents.hh"
 #include "HEJ/utility.hh"
 #include "HEJ/Tensor.hh"
 #include "HEJ/Constants.hh"
 
 #include <array>
 
 #include <iostream>
 
 
 
 namespace { // Helper Functions
   // FKL W Helper Functions
   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);
       joi(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);
       joi(pnu,helnu,pin,helin,temp3);
       joi(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];
     }
   }
 
   double WProp (const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl){
     COM propW = COM(0.,-1.)/((pl+plbar).m2() -HEJ::MW*HEJ::MW + COM(0.,1.)*HEJ::MW*HEJ::GammaW);
     double PropFactor=(propW*conj(propW)).real();
     return PropFactor;
   }
 
 CCurrent jW (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin)
 {
 
   COM cur[4];
 
   cur[0]=0.;
   cur[1]=0.;
   cur[2]=0.;
   cur[3]=0.;
   CCurrent sum(0.,0.,0.,0.);
 
   // 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());
 
     CCurrent temp2,temp3,temp5;
     CCurrent t65 = joo(pnu,helnu,pe,hele);
     CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
     CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
 
     COM brac615=t65.dot(vout);
     COM brac645=t65.dot(vin);
 
 
     // prod1565 and prod6465 are zero for Ws (not Zs)!!
     temp2 = joo(pout,helout,pnu,helout);
     COM prod1665=temp2.dot(t65);
     temp3 = j(pe,helin,pin,helin);
     COM prod5465=temp3.dot(t65);
 
     temp2=joo(pout,helout,pe,helout);
     temp3=j(pnu,helnu,pin,helin);
     temp5=j(pout,helout,pin,helin);
 
     CCurrent term1,term2,term3;
     term1=(2.*brac615/ta+2.*brac645/tb)*temp5;
     term2=(prod1665/ta)*temp3;
     term3=(-prod5465/tb)*temp2;
 
     sum=term1+term2+term3;
   }
 
   return sum;
 }
 
 CCurrent jWbar (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin)
 {
 
     COM cur[4];
 
     cur[0]=0.;
     cur[1]=0.;
     cur[2]=0.;
     cur[3]=0.;
     CCurrent sum(0.,0.,0.,0.);
 
   // 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());
 
     CCurrent temp2,temp3,temp5;
     CCurrent t65 = joo(pnu,helnu,pe,hele);
     CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
     CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
 
     COM brac615=t65.dot(vout);
     COM brac645=t65.dot(vin);
 
     // prod1565 and prod6465 are zero for Ws (not Zs)!!
     temp2 = joo(pe,helout,pout,helout);  //  temp2 is <5|alpha|1>
     COM prod5165=temp2.dot(t65);
     temp3 = jio(pin,helin,pnu,helin);      // temp3 is <4|alpha|6>
     COM prod4665=temp3.dot(t65);
 
     temp2=joo(pnu,helout,pout,helout);  // temp2 is now <6|mu|1>
     temp3=jio(pin,helin,pe,helin);        // temp3 is now <4|mu|5>
     temp5=jio(pin,helin,pout,helout);     //  temp5 is <4|mu|1>
 
     CCurrent term1,term2,term3;
     term1 =(-2.*brac615/ta-2.*brac645/tb)*temp5;
     term2 =(-prod5165/ta)*temp3;
     term3 =(prod4665/tb)*temp2;
 
     sum = term1 + term2 + term3;
   }
 
   return sum;
 }
 
 
   // Extremal quark current with W emission. Using Tensor class rather than CCurrent
   Tensor <1,4> jW(HLV pin, HLV pout, HLV plbar, HLV pl, bool aqline){
     // Build the external quark line W Emmision
     Tensor<1,4> ABCurr = TCurrent(pl, false, plbar, false);
     Tensor<1,4> Tp4W = Construct1Tensor((pout+pl+plbar));//p4+pw
     Tensor<1,4> TpbW = Construct1Tensor((pin-pl-plbar));//pb-pw
 
     Tensor<3,4> J4bBlank;
     if (aqline){
       J4bBlank = T3Current(pin,false,pout,false);
     }
     else{
       J4bBlank = T3Current(pout,false,pin,false);
     }
     double t4AB = (pout+pl+plbar).m2();
     double tbAB = (pin-pl-plbar).m2();
 
     Tensor<2,4> J4b1 = (J4bBlank.contract(Tp4W,2))/t4AB;
     Tensor<2,4> J4b2 = (J4bBlank.contract(TpbW,2))/tbAB;
 
     Tensor<2,4> T4bmMom(0.);
 
     if (aqline){
       for(int mu=0; mu<4;mu++){
         for(int nu=0;nu<4;nu++){
           T4bmMom.Set(mu,nu, (J4b1.at(nu,mu) + J4b2.at(mu,nu))*(COM(0,-1)));
         }
       }
     }
     else{
       for(int mu=0; mu<4;mu++){
         for(int nu=0;nu<4;nu++){
           T4bmMom.Set(nu,mu, (J4b1.at(nu,mu) + J4b2.at(mu,nu))*(COM(0,1)));
         }
       }
     }
     Tensor<1,4> T4bm = T4bmMom.contract(ABCurr,1);
 
     return T4bm;
   }
 
 
 
   // Relevant W+Jets Unordered Contribution Helper Functions
   // W+Jets Uno
   double jM2Wuno(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector plbar, CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pa, bool h1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector pb, bool h2, bool pol)
   {
     static bool is_sigma_index_set(false);
 
     if(!is_sigma_index_set){
       //std::cout<<"Setting sigma_index...." << std::endl;
       if(init_sigma_index())
   is_sigma_index_set = true;
       else
   return 0.;
     }
 
     CLHEP::HepLorentzVector pW = pl+plbar;
     CLHEP::HepLorentzVector q1g=pa-pW-p1-pg;
     CLHEP::HepLorentzVector q1 = pa-p1-pW;
     CLHEP::HepLorentzVector q2 = p2-pb;
 
     const double taW  = (pa-pW).m2();
     const double taW1 = (pa-pW-p1).m2();
     const double tb2  = (pb-p2).m2();
     const double tb2g = (pb-p2-pg).m2();
     const double s1W  = (p1+pW).m2();
     const double s1gW = (p1+pW+pg).m2();
     const double s1g  = (p1+pg).m2();
     const double tag  = (pa-pg).m2();
     const double taWg = (pa-pW-pg).m2();
 
     //use p1 as ref vec in pol tensor
     Tensor<1,4> epsg = eps(pg,p2,pol);
     Tensor<1,4> epsW = TCurrent(pl,false,plbar,false);
     Tensor<1,4> j2b = TCurrent(p2,h2,pb,h2);
 
     Tensor<1,4> Tq1q2 = Construct1Tensor((q1+q2)/taW1 + (pb/pb.dot(pg)
       + p2/p2.dot(pg)) * tb2/(2*tb2g));
     Tensor<1,4> Tq1g = Construct1Tensor((-pg-q1))/taW1;
     Tensor<1,4> Tq2g = Construct1Tensor((pg-q2));
     Tensor<1,4> TqaW = Construct1Tensor((pa-pW));//pa-pw
     Tensor<1,4> Tqag = Construct1Tensor((pa-pg));
     Tensor<1,4> TqaWg = Construct1Tensor((pa-pg-pW));
     Tensor<1,4> Tp1g = Construct1Tensor((p1+pg));
     Tensor<1,4> Tp1W = Construct1Tensor((p1+pW));//p1+pw
     Tensor<1,4> Tp1gW = Construct1Tensor((p1+pg+pW));//p1+pw+pg
 
     Tensor<2,4> g=Metric();
 
     Tensor<3,4> J31a = T3Current(p1, h1, pa, h1);
     Tensor<2,4> J2_qaW =J31a.contract(TqaW/taW, 2);
     Tensor<2,4> J2_p1W =J31a.contract(Tp1W/s1W, 2);
     Tensor<3,4> L1a =J2_qaW.leftprod(Tq1q2);
     Tensor<3,4> L1b =J2_p1W.leftprod(Tq1q2);
     Tensor<3,4> L2a = J2_qaW.leftprod(Tq1g);
     Tensor<3,4> L2b = J2_p1W.leftprod(Tq1g);
     Tensor<3,4> L3 = (g.rightprod(J2_qaW.contract(Tq2g,1)+J2_p1W.contract(Tq2g,2)))/taW1;
     Tensor<3,4> L(0.);
 
     Tensor<5,4> J51a = T5Current(p1, h1, pa, h1);
 
     Tensor<4,4> J_qaW = J51a.contract(TqaW,4);
     Tensor<4,4> J_qag = J51a.contract(Tqag,4);
     Tensor<4,4> J_p1gW = J51a.contract(Tp1gW,4);
 
     Tensor<3,4> U1a = J_qaW.contract(Tp1g,2);
     Tensor<3,4> U1b = J_p1gW.contract(Tp1g,2);
     Tensor<3,4> U1c = J_p1gW.contract(Tp1W,2);
     Tensor<3,4> U1(0.);
 
     Tensor<3,4> U2a = J_qaW.contract(TqaWg,2);
     Tensor<3,4> U2b = J_qag.contract(TqaWg,2);
     Tensor<3,4> U2c = J_qag.contract(Tp1W,2);
     Tensor<3,4> U2(0.);
 
     for(int nu=0; nu<4;nu++){
       for(int mu=0;mu<4;mu++){
         for(int rho=0;rho<4;rho++){
           L.Set(nu, mu, rho, L1a.at(nu,mu,rho) + L1b.at(nu,rho,mu)
             + L2a.at(mu,nu,rho) + L2b.at(mu,rho,nu) + L3.at(mu,nu,rho));
           U1.Set(nu, mu, rho, U1a.at(nu, mu, rho) / (s1g*taW)
             + U1b.at(nu,rho,mu) / (s1g*s1gW) + U1c.at(rho,nu,mu) / (s1W*s1gW));
           U2.Set(nu,mu,rho,U2a.at(mu,nu,rho) / (taWg*taW)
             + U2b.at(mu,rho,nu) / (taWg*tag) + U2c.at(rho,mu,nu) / (s1W*tag));
         }
       }
     }
 
     COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1)).at(0);
     COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1)).at(0);
 
     double amp = HEJ::C_A*HEJ::C_F*HEJ::C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
 
     double t1 = q1g.m2();
     double t2 = q2.m2();
 
     double WPropfact = WProp(plbar, pl);
 
     //Divide by WProp
     amp*=WPropfact;
 
     //Divide by t-channels
     amp/=(t1*t2);
 
     //Average over initial states
     amp/=(4.*HEJ::C_A*HEJ::C_A);
 
     return amp;
   }
 
 
 
   // Relevant Wqqx Helper Functions.
   //g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
   Tensor <1,4> gtqqxW(CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar){
 
     double s2AB=(pl+plbar+pq).m2();
     double s3AB=(pl+plbar+pqbar).m2();
 
     Tensor<1,4> Tpq = Construct1Tensor(pq);
     Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
     Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
 
     // Define llx current.
     Tensor<1,4> ABCur = TCurrent(pl, false, plbar, false);
 
     //blank 3 Gamma Current
     Tensor<3,4> JV23 = T3Current(pq,false,pqbar,false);
 
     // Components of g->qqW before W Contraction
     Tensor<2,4> JV1 = JV23.contract((Tpq + TAB),2)/(s2AB);
     Tensor<2,4> JV2 = JV23.contract((Tpqbar + TAB),2)/(s3AB);
 
     // g->qqW Current. Note Minus between terms due to momentum flow.
     // Also note: (-I)^2 from W vert. (I) from Quark prop.
     Tensor<1,4> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
 
     return JVCur;
   }
 
   // Helper Functions Calculate the Crossed Contribution
   Tensor <2,4> MCrossW(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
 
     // Useful propagator factors
     double s2AB=(pl+plbar+pq).m2();
     double s3AB=(pl+plbar+pqbar).m2();
 
     CLHEP::HepLorentzVector q1, q3;
     q1=pa;
     for(int i=0; i<nabove+1;i++){
       q1=q1-partons.at(i);
     }
     q3 = q1 - pq - pqbar - pl - plbar;
 
     double tcro1=(q3+pq).m2();
     double tcro2=(q1-pqbar).m2();
 
     Tensor<1,4> Tp1 = Construct1Tensor(p1);
     Tensor<1,4> Tp4 = Construct1Tensor(p4);
     Tensor<1,4> Tpa = Construct1Tensor(pa);
     Tensor<1,4> Tpb = Construct1Tensor(pb);
     Tensor<1,4> Tpq = Construct1Tensor(pq);
     Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
     Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
     Tensor<1,4> Tq1 = Construct1Tensor(q1);
     Tensor<1,4> Tq3 = Construct1Tensor(q3);
     Tensor<2,4> g=Metric();
 
     // Define llx current.
     Tensor<1,4> ABCur = TCurrent(pl, false, plbar,false);
 
     //Blank 5 gamma Current
     Tensor<5,4> J523 = T5Current(pq,false,pqbar,false);
 
     // 4 gamma currents (with 1 contraction already).
     Tensor<4,4> J_q3q = J523.contract((Tq3+Tpq),2);
     Tensor<4,4> J_2AB = J523.contract((Tpq+TAB),2);
 
     // Components of Crossed Vertex Contribution
     Tensor<3,4> Xcro1 = J_q3q.contract((Tpqbar + TAB),3);
     Tensor<3,4> Xcro2 = J_q3q.contract((Tq1-Tpqbar),3);
     Tensor<3,4> Xcro3 = J_2AB.contract((Tq1-Tpqbar),3);
 
     // Term Denominators Taken Care of at this stage
     Tensor<2,4> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
     Tensor<2,4> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
     Tensor<2,4> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
 
     //Initialise the Crossed Vertex Object
     Tensor<2,4> Xcro(0.);
 
     for(int mu=0; mu<4;mu++){
       for(int nu=0;nu<4;nu++){
         Xcro.Set(mu,nu, -(-Xcro1Cont.at(nu,mu)+Xcro2Cont.at(nu,mu)+Xcro3Cont.at(nu,mu)));
       }
     }
 
 
     return Xcro;
   }
 
   // Helper Functions Calculate the Uncrossed Contribution
   Tensor <2,4> MUncrossW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
 
     double s2AB=(pl+plbar+pq).m2();
     double s3AB=(pl+plbar+pqbar).m2();
 
     CLHEP::HepLorentzVector q1, q3;
     q1=pa;
     for(int i=0; i<nabove+1;i++){
       q1=q1-partons.at(i);
     }
     q3 = q1 - pl - plbar - pq - pqbar;
     double tunc1 = (q1-pq).m2();
     double tunc2 = (q3+pqbar).m2();
 
     Tensor<1,4> Tp1 = Construct1Tensor(p1);
     Tensor<1,4> Tp4 = Construct1Tensor(p4);
     Tensor<1,4> Tpa = Construct1Tensor(pa);
     Tensor<1,4> Tpb = Construct1Tensor(pb);
     Tensor<1,4> Tpq = Construct1Tensor(pq);
     Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
     Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
     Tensor<1,4> Tq1 = Construct1Tensor(q1);
     Tensor<1,4> Tq3 = Construct1Tensor(q3);
     Tensor<2,4> g=Metric();
 
 
     // Define llx current.
     Tensor<1,4> ABCur = TCurrent(pl, false, plbar, false);
 
     //Blank 5 gamma Current
     Tensor<5,4> J523 = T5Current(pq,false,pqbar,false);
 
     // 4 gamma currents (with 1 contraction already).
     Tensor<4,4> J_2AB = J523.contract((Tpq+TAB),2);
     Tensor<4,4> J_q1q = J523.contract((Tq1-Tpq),2);
 
     // 2 Contractions taken care of.
     Tensor<3,4> Xunc1 = J_2AB.contract((Tq3+Tpqbar),3);
     Tensor<3,4> Xunc2 = J_q1q.contract((Tq3+Tpqbar),3);
     Tensor<3,4> Xunc3 = J_q1q.contract((Tpqbar+TAB),3);
 
     // Term Denominators Taken Care of at this stage
     Tensor<2,4> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
     Tensor<2,4> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
     Tensor<2,4> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
 
     //Initialise the Uncrossed Vertex Object
     Tensor<2,4> Xunc(0.);
 
     for(int mu=0; mu<4;mu++){
       for(int nu=0;nu<4;nu++){
         Xunc.Set(mu,nu,-(- Xunc1Cont.at(mu,nu)+Xunc2Cont.at(mu,nu) +Xunc3Cont.at(mu,nu)));
       }
     }
 
     return Xunc;
   }
 
 
   // Helper Functions Calculate the g->qqxW (Eikonal) Contributions
   Tensor <2,4> MSymW(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
 
     double sa2=(pa+pq).m2();
     double s12=(p1+pq).m2();
     double sa3=(pa+pqbar).m2();
     double s13=(p1+pqbar).m2();
     double saA=(pa+pl).m2();
     double s1A=(p1+pl).m2();
     double saB=(pa+plbar).m2();
     double s1B=(p1+plbar).m2();
     double sb2=(pb+pq).m2();
     double s42=(p4+pq).m2();
     double sb3=(pb+pqbar).m2();
     double s43=(p4+pqbar).m2();
     double sbA=(pb+pl).m2();
     double s4A=(p4+pl).m2();
     double sbB=(pb+plbar).m2();
     double s4B=(p4+plbar).m2();
     double s23AB=(pl+plbar+pq+pqbar).m2();
 
     CLHEP::HepLorentzVector q1,q3;
     q1=pa;
     for(int i=0;i<nabove+1;i++){
       q1-=partons.at(i);
     }
     q3=q1-pq-pqbar-plbar-pl;
     double t1 = (q1).m2();
     double t3 = (q3).m2();
 
     //Define Tensors to be used
     Tensor<1,4> Tp1 = Construct1Tensor(p1);
     Tensor<1,4> Tp4 = Construct1Tensor(p4);
     Tensor<1,4> Tpa = Construct1Tensor(pa);
     Tensor<1,4> Tpb = Construct1Tensor(pb);
     Tensor<1,4> Tpq = Construct1Tensor(pq);
     Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
     Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
     Tensor<1,4> Tq1 = Construct1Tensor(q1);
     Tensor<1,4> Tq3 = Construct1Tensor(q3);
     Tensor<2,4> g=Metric();
 
     // g->qqW Current (Factors of sqrt2 dealt with in this function.)
     Tensor<1,4> JV = gtqqxW(pq,pqbar,pl,plbar);
 
     // 1a gluon emisson Contribution
     Tensor<3,4> X1a = g.rightprod(Tp1*(t1/(s12+s13+s1A+s1B)) + Tpa*(t1/(sa2+sa3+saA+saB)));
     Tensor<2,4> X1aCont = X1a.contract(JV,3);
 
     //4b gluon emission Contribution
     Tensor<3,4> X4b = g.rightprod(Tp4*(t3/(s42+s43+s4A+s4B)) + Tpb*(t3/(sb2+sb3+sbA+sbB)));
     Tensor<2,4> X4bCont = X4b.contract(JV,3);
 
     //Set up each term of 3G diagram.
     Tensor<3,4> X3g1 = g.leftprod(Tq1+Tpq+Tpqbar+TAB);
     Tensor<3,4> X3g2 = g.leftprod(Tq3-Tpq-Tpqbar-TAB);
     Tensor<3,4> X3g3 = g.leftprod((Tq1+Tq3));
 
     // Note the contraction of indices changes term by term
     Tensor<2,4> X3g1Cont = X3g1.contract(JV,3);
     Tensor<2,4> X3g2Cont = X3g2.contract(JV,2);
     Tensor<2,4> X3g3Cont = X3g3.contract(JV,1);
 
     // XSym is an amalgamation of x1a, X4b and X3g. Makes sense from a colour factor point of view.
     Tensor<2,4>Xsym(0.);
 
     for(int mu=0; mu<4;mu++){
       for(int nu=0;nu<4;nu++){
   Xsym.Set(mu,nu, (X3g1Cont.at(nu,mu) + X3g2Cont.at(mu,nu) - X3g3Cont.at(nu,mu))
      + (X1aCont.at(mu,nu) - X4bCont.at(mu,nu)) );
       }
     }
     return Xsym/s23AB;
   }
 
   Tensor <2,4> MCross(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
 
     CLHEP::HepLorentzVector q1;
     q1=pa;
     for(int i=0;i<nabove+1;i++){
       q1-=partons.at(i);
     }
 
     double t2=(q1-pqbar).m2();
 
     Tensor<1,4> Tq1 = Construct1Tensor(q1-pqbar);
 
     //Blank 3 gamma Current
     Tensor<3,4> J323 = T3Current(pq,hq,pqbar,hq);
 
     // 2 gamma current (with 1 contraction already).
     Tensor<2,4> XCroCont = J323.contract((Tq1),2)/(t2);
 
     //Initialise the Crossed Vertex
     Tensor<2,4> Xcro(0.);
 
     for(int mu=0; mu<4;mu++){
       for(int nu=0;nu<4;nu++){
         Xcro.Set(mu,nu, (XCroCont.at(nu,mu)));
       }
     }
 
     return Xcro;
   }
 
 
   // Helper Functions Calculate the Uncrossed Contribution
   Tensor <2,4> MUncross(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
 
     CLHEP::HepLorentzVector q1;
     q1=pa;
     for(int i=0;i<nabove+1;i++){
       q1-=partons.at(i);
     }
     double t2 = (q1-pq).m2();
 
     Tensor<1,4> Tq1 = Construct1Tensor(q1-pq);
 
     //Blank 3 gamma Current
     Tensor<3,4> J323 = T3Current(pq,hq,pqbar,hq);
 
     // 2 gamma currents (with 1 contraction already).
     Tensor<2,4> XUncCont = J323.contract((Tq1),2)/t2;
 
     //Initialise the Uncrossed Vertex
     Tensor<2,4> Xunc(0.);
 
     for(int mu=0; mu<4;mu++){
       for(int nu=0;nu<4;nu++){
   Xunc.Set(mu,nu,-(XUncCont.at(mu,nu)));
       }
     }
 
     return Xunc;
   }
 
 
   // Helper Functions Calculate the Eikonal Contributions
   Tensor <2,4> MSym(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
 
     CLHEP::HepLorentzVector q1, q3;
     q1=pa;
     for(int i=0;i<nabove+1;i++){
       q1-=partons.at(i);
     }
     q3 = q1-pq-pqbar;
     double t1 = (q1).m2();
     double t3 = (q3).m2();
 
     double s23 = (pq+pqbar).m2();
     double sa2 = (pa+pq).m2();
     double sa3 = (pa+pqbar).m2();
     double s12 = (p1+pq).m2();
     double s13 = (p1+pqbar).m2();
     double sb2 = (pb+pq).m2();
     double sb3 = (pb+pqbar).m2();
     double s42 = (p4+pq).m2();
     double s43 = (p4+pqbar).m2();
 
     //Define Tensors to be used
     Tensor<1,4> Tp1 = Construct1Tensor(p1);
     Tensor<1,4> Tp4 = Construct1Tensor(p4);
     Tensor<1,4> Tpa = Construct1Tensor(pa);
     Tensor<1,4> Tpb = Construct1Tensor(pb);
     Tensor<1,4> Tpq = Construct1Tensor(pq);
     Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
     Tensor<1,4> Tq1 = Construct1Tensor(q1);
     Tensor<1,4> Tq3 = Construct1Tensor(q3);
     Tensor<2,4> g=Metric();
 
     Tensor<1,4> qqxCur = TCurrent(pq, hq, pqbar, hq);
 
     // // 1a gluon emisson Contribution
     Tensor<3,4> X1a = g.rightprod(Tp1*(t1/(s12+s13))+Tpa*(t1/(sa2+sa3)));
     Tensor<2,4> X1aCont = X1a.contract(qqxCur,3);
 
     // //4b gluon emission Contribution
     Tensor<3,4> X4b = g.rightprod(Tp4*(t3/(s42+s43)) + Tpb*(t3/(sb2+sb3)));
     Tensor<2,4> X4bCont = X4b.contract(qqxCur,3);
 
     // New Formulation Corresponding to New Analytics
     Tensor<3,4> X3g1 = g.leftprod(Tq1+Tpq+Tpqbar);
     Tensor<3,4> X3g2 = g.leftprod(Tq3-Tpq-Tpqbar);
     Tensor<3,4> X3g3 = g.leftprod((Tq1+Tq3));
 
     // Note the contraction of indices changes term by term
     Tensor<2,4> X3g1Cont = X3g1.contract(qqxCur,3);
     Tensor<2,4> X3g2Cont = X3g2.contract(qqxCur,2);
     Tensor<2,4> X3g3Cont = X3g3.contract(qqxCur,1);
 
     Tensor<2,4>Xsym(0.);
 
     for(int mu=0; mu<4;mu++){
       for(int nu=0;nu<4;nu++){
         Xsym.Set(mu, nu, COM(0,1) * ( (X3g1Cont.at(nu,mu) + X3g2Cont.at(mu,nu)
           - X3g3Cont.at(nu,mu)) + (X1aCont.at(mu,nu) - X4bCont.at(mu,nu)) ) );
       }
     }
     return Xsym/s23;
   }
 } // Anonymous Namespace helper functions
 
 // W+Jets FKL Contributions
 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);
   joi(p2out,true,p2in,true,mj2p);
   joi(p2out,false,p2in,false,mj2m);
 
   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);
 
   double WPropfact = WProp(pe, pnu);
 
-  // Leave division by colour and Helicity avg until Tree files
-  // Leave multi. of couplings to later
+  // Division by colour and Helicity average (Nc2-1)(4)
   // Multiply by Cf^2
-  return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());
+  return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
 
 }
 
 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);
 
   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);
 
   double WPropfact = WProp(pe, pnu);
 
-  // Leave division by colour and Helicity avg until Tree files
-  // Leave multi. of couplings to later
+  // Division by colour and Helicity average (Nc2-1)(4)
   // Multiply by Cf^2
-  return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());
+  return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
 
 }
 
 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);
   joi(p2out,true,p2in,true,mj2p);
   joi(p2out,false,p2in,false,mj2m);
 
   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);
 
   double WPropfact = WProp(pe, pnu);
 
-  // Leave division by colour and Helicity avg until Tree files
-  // Leave multi. of couplings to later
+  // Division by colour and Helicity average (Nc2-1)(4)
   // Multiply by Cf^2
-  return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());
+  return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
 
 }
 
 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);
 
   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);
 
   double WPropfact = WProp(pe, pnu);
 
-  // Leave division by colour and Helicity avg until Tree files
-  // Leave multi. of couplings to later
+  // Division by colour and Helicity average (Nc2-1)(4)
   // Multiply by Cf^2
-  return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2());
+  return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
 
 }
 
 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);
 
   joi(p2out,true,p2in,true,mj2p);
   joi(p2out,false,p2in,false,mj2m);
 
   // mj1m.mj2p
   COM Mmp=cdot(mj1m,mj2p);
 
   // mj1m.mj2m
   COM Mmm=cdot(mj1m,mj2m);
 
   const double K = K_g(p2out, p2in);
 
   // sum of spinor strings ||^2
   double a2Mmp=abs2(Mmp);
   double a2Mmm=abs2(Mmm);
   double sst = K/HEJ::C_A*(a2Mmp+a2Mmm);
 
   double WPropfact = WProp(pe, pnu);
 
-  // Leave division by colour and Helicity avg until Tree files
-  // Leave multi. of couplings to later
+  // Division by colour and Helicity average (Nc2-1)(4)
   // Multiply by Cf*Ca=4
-  return HEJ::C_F*HEJ::C_A*WPropfact*sst/(q1.m2()*q2.m2());
+  return HEJ::C_F*HEJ::C_A*WPropfact*sst/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
 
 }
 
 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);
 
   joi(p2out,true,p2in,true,mj2p);
   joi(p2out,false,p2in,false,mj2m);
 
   // mj1m.mj2p
   COM Mmp=cdot(mj1m,mj2p);
 
   // mj1m.mj2m
   COM Mmm=cdot(mj1m,mj2m);
 
   const double K = K_g(p2out, p2in);
 
   // sum of spinor strings ||^2
   double a2Mmp=abs2(Mmp);
   double a2Mmm=abs2(Mmm);
   double sst = K/HEJ::C_A*(a2Mmp+a2Mmm);
 
   double WPropfact = WProp(pe, pnu);
 
-//   // Leave division by colour and Helicity avg until Tree files
-  // Leave multi. of couplings to later
+  // Division by colour and Helicity average (Nc2-1)(4)
   // Multiply by Cf*Ca=4
-  return HEJ::C_F*HEJ::C_A*WPropfact*sst/(q1.m2()*q2.m2());
+  return HEJ::C_F*HEJ::C_A*WPropfact*sst/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
 
 }
 
 
 // W+Jets Unordered Contributions
 //qQ->qQWg_unob
 double junobMWqQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   CCurrent mj1m,mj2p,mj2m;
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out-pg);
   CLHEP::HepLorentzVector q3=-(p2in-p2out);
 
   mj1m=jW(p1out,false,pe,false,pnu,false,p1in,false);
   mj2p=j(p2out,true,p2in,true);
   mj2m=j(p2out,false,p2in,false);
 
 
   // Dot products of these which occur again and again
   COM MWmp=mj1m.dot(mj2p);  // And now for the Higgs ones
   COM MWmm=mj1m.dot(mj2m);
 
   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);
 
 
   Lmm=((-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmm/2.))/q3.m2();
   Lmp=((-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmp/2.))/q3.m2();
 
   U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
   U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
 
   U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
   U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
 
 
   double amm,amp;
 
   amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
   amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
 
   double ampsq=-(amm+amp);
   //Divide by WProp
   double WPropfact = WProp(pe, pnu);
   ampsq*=WPropfact;
 
 
   // Now add the t-channels
   double th=q2.m2()*q1.m2();
   ampsq/=th;
   ampsq/=16.;
 
 
   return ampsq;
 
 
 }
 
 //qQbar->qQbarWg_unob
 double junobMWqQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   CCurrent mj1m,mj2p,mj2m;
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out-pg);
   CLHEP::HepLorentzVector q3=-(p2in-p2out);
 
   mj1m=jW(p1out,false,pe,false,pnu,false,p1in,false);
   mj2p=jio(p2in,true,p2out,true);
   mj2m=jio(p2in,false,p2out,false);
 
 
   // Dot products of these which occur again and again
   COM MWmp=mj1m.dot(mj2p);  // And now for the Higgs ones
   COM MWmm=mj1m.dot(mj2m);
 
   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);
 
 
   Lmm=((-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmm/2.))/q3.m2();
   Lmp=((-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmp/2.))/q3.m2();
 
   U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
   U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
 
   U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
   U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
 
   double amm,amp;
 
   amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
   amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
 
   double ampsq=-(amm+amp);
   //Divide by WProp
   double WPropfact = WProp(pe, pnu);
   ampsq*=WPropfact;
 
   // Now add the t-channels
   double th=q2.m2()*q1.m2();
   ampsq/=th;
   ampsq/=16.;
 
 
   return ampsq;
 
 
 }
 
 //qbarQ->qbarQWg_unob
 double junobMWqbarQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   CCurrent mj1m,mj2p,mj2m;
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out-pg);
   CLHEP::HepLorentzVector q3=-(p2in-p2out);
 
   mj1m=jWbar(p1out,false,pe,false,pnu,false,p1in,false);
   mj2p=j(p2out,true,p2in,true);
   mj2m=j(p2out,false,p2in,false);
 
 
   // Dot products of these which occur again and again
   COM MWmp=mj1m.dot(mj2p);  // And now for the Higgs ones
   COM MWmm=mj1m.dot(mj2m);
 
   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);
 
 
   Lmm=((-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmm/2.))/q3.m2();
   Lmp=((-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmp/2.))/q3.m2();
 
   U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
   U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
 
   U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
   U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
 
 
   double amm,amp;
 
   amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
   amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
 
   double ampsq=-(amm+amp);
   //Divide by WProp
   double WPropfact = WProp(pe, pnu);
   ampsq*=WPropfact;
 
   // Now add the t-channels
   double th=q2.m2()*q1.m2();
   ampsq/=th;
   ampsq/=16.;
 
   return ampsq;
 
 
 }
 
 //qbarQbar->qbarQbarWg_unob
 double junobMWqbarQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   CCurrent mj1m,mj2p,mj2m;
   CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
   CLHEP::HepLorentzVector q2=-(p2in-p2out-pg);
   CLHEP::HepLorentzVector q3=-(p2in-p2out);
 
   mj1m=jWbar(p1out,false,pe,false,pnu,false,p1in,false);
   mj2p=jio(p2in,true,p2out,true);
   mj2m=jio(p2in,false,p2out,false);
 
 
   // Dot products of these which occur again and again
   COM MWmp=mj1m.dot(mj2p);  // And now for the Higgs ones
   COM MWmm=mj1m.dot(mj2m);
 
   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);
 
 
   Lmm=((-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmm/2.))/q3.m2();
   Lmp=((-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmp/2.))/q3.m2();
 
   U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
   U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
 
   U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
   U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
 
 
   double amm,amp;
 
   amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
   amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
 
   double ampsq=-(amm+amp);
   //Divide by WProp
   double WPropfact = WProp(pe, pnu);
   ampsq*=WPropfact;
 
   // Now add the t-channels
   double th=q2.m2()*q1.m2();
   ampsq/=th;
   ampsq/=16.;
 
   return ampsq;
 
 
 }
 
 ////////////////////////////////////////////////////////////////////
 //qQ->qQWg_unof
 double junofMWgqQ (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   CCurrent mj2m,mj1p,mj1m;
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector qg=p1in-p1out-pg;
   CLHEP::HepLorentzVector q2=-(p2in-p2out-pe-pnu);
 
   mj2m=jW(p2out,false,pe,false,pnu,false,p2in,false);
   mj1p=j(p1out,true,p1in,true);
   mj1m=j(p1out,false,p1in,false);
 
 
   // Dot products of these which occur again and again
   COM MWpm=mj1p.dot(mj2m);  // And now for the Higgs ones
   COM MWmm=mj1m.dot(mj2m);
 
   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*(MWmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWmm/2.))/q1.m2();
   Lpm=(qsum*(MWpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWpm/2.))/q1.m2();
 
 
   U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*MWmm)/(p1out+pg).m2();
   U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*MWpm)/(p1out+pg).m2();
 
   U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*MWmm)/(p1in-pg).m2();
   U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*MWpm)/(p1in-pg).m2();
 
   double amm,apm;
   amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
   apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
 
   double ampsq=-(apm+amm);
   //Divide by WProp
   double WPropfact = WProp(pe, pnu);
   ampsq*=WPropfact;
 
   // Now add the t-channels
   double th=q2.m2()*qg.m2();
   ampsq/=th;
   ampsq/=16.;
 
   return ampsq;
 
 
 }
 
 //qQbar->qQbarWg_unof
 double junofMWgqQbar (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   CCurrent mj2m,mj1p,mj1m;
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector qg=p1in-p1out-pg;
   CLHEP::HepLorentzVector q2=-(p2in-p2out-pe-pnu);
 
   mj2m=jWbar(p2out,false,pe,false,pnu,false,p2in,false);
   mj1p=j(p1out,true,p1in,true);
   mj1m=j(p1out,false,p1in,false);
 
 
   // Dot products of these which occur again and again
   COM MWpm=mj1p.dot(mj2m);  // And now for the Higgs ones
   COM MWmm=mj1m.dot(mj2m);
 
   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*(MWmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWmm/2.))/q1.m2();
   Lpm=(qsum*(MWpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWpm/2.))/q1.m2();
 
 
   U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*MWmm)/(p1out+pg).m2();
   U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*MWpm)/(p1out+pg).m2();
 
   U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*MWmm)/(p1in-pg).m2();
   U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*MWpm)/(p1in-pg).m2();
 
   double amm,apm;
 
   amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
   apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
 
   double ampsq=-(apm+amm);
   //Divide by WProp
   double WPropfact = WProp(pe, pnu);
   ampsq*=WPropfact;
 
   // Now add the t-channels
   double th=q2.m2()*qg.m2();
   ampsq/=th;
   ampsq/=16.;
 
   return ampsq;
 }
 
 //qbarQ->qbarQWg_unof
 double junofMWgqbarQ (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   CCurrent mj2m,mj1p,mj1m;
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector qg=p1in-p1out-pg;
   CLHEP::HepLorentzVector q2=-(p2in-p2out-pe-pnu);
 
   mj2m=jW(p2out,false,pe,false,pnu,false,p2in,false);
   mj1p=jio(p1in,true,p1out,true);
   mj1m=jio(p1in,false,p1out,false);
 
 
   // Dot products of these which occur again and again
   COM MWpm=mj1p.dot(mj2m);  // And now for the Higgs ones
   COM MWmm=mj1m.dot(mj2m);
 
   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*(MWmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWmm/2.))/q1.m2();
   Lpm=(qsum*(MWpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWpm/2.))/q1.m2();
 
 
   U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*MWmm)/(p1out+pg).m2();
   U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*MWpm)/(p1out+pg).m2();
 
   U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*MWmm)/(p1in-pg).m2();
   U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*MWpm)/(p1in-pg).m2();
 
   double amm,apm;
   amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
   apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
 
   double ampsq=-(apm+amm);
   //Divide by WProp
   double WPropfact = WProp(pe, pnu);
   ampsq*=WPropfact;
 
   // Now add the t-channels
   double th=q2.m2()*qg.m2();
   ampsq/=th;
   ampsq/=16.;
 
   return ampsq;
 
 
 }
 
 //qbarQbar->qbarQbarWg_unof
 double junofMWgqbarQbar (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in)
 // Calculates the square of the current contractions for qQ->qenuQ scattering
 // p1: quark (with W emittance)
 // p2: Quark
 {
   CCurrent mj2m,mj1p,mj1m;
   CLHEP::HepLorentzVector q1=p1in-p1out;
   CLHEP::HepLorentzVector qg=p1in-p1out-pg;
   CLHEP::HepLorentzVector q2=-(p2in-p2out-pe-pnu);
 
   mj2m=jWbar(p2out,false,pe,false,pnu,false,p2in,false);
   mj1p=jio(p1in,true,p1out,true);
   mj1m=jio(p1in,false,p1out,false);
 
 
   // Dot products of these which occur again and again
   COM MWpm=mj1p.dot(mj2m);  // And now for the Higgs ones
   COM MWmm=mj1m.dot(mj2m);
 
   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*(MWmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWmm/2.))/q1.m2();
   Lpm=(qsum*(MWpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWpm/2.))/q1.m2();
 
   U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*MWmm)/(p1out+pg).m2();
   U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*MWpm)/(p1out+pg).m2();
 
   U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*MWmm)/(p1in-pg).m2();
   U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*MWpm)/(p1in-pg).m2();
 
   double amm,apm;
 
   amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
   apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
 
 
 
   double ampsq=-(apm+amm);
   //Divide by WProp
   double WPropfact = WProp(pe, pnu);
   ampsq*=WPropfact;
 
 
   // Now add the t-channels
   double th=q2.m2()*qg.m2();
   ampsq/=th;
   ampsq/=16.;
   return ampsq;
 
 
 }
 
 
 
 ///TODO make this comment more visible
 /// Naming scheme jM2-Wuno-g-({q/qbar}{Q/Qbar/g})
 ///TODO Spit naming for more complicated functions?
 /// e.g. jM2WqqtoqQQq -> jM2_Wqq_to_qQQq
 double jM2WunogqQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   return ME2;
 
 }
 
 //same as function above but actually obtaining the antiquark line by crossing symmetry, where p1out and p1in are expected to be negative.
 //should give same result as jM2WunogqbarQ below (verified)
 double jM2WunogqQ_crossqQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   return ME2;
 
 }
 
 double jM2WunogqQbar(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   return ME2;
 
 }
 
 double jM2Wunogqg(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   double ratio; // p2-/pb- in the notes
   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 cam = ( (HEJ::C_A - 1/HEJ::C_A)*(ratio + 1./ratio)/2. + 1/HEJ::C_A)/HEJ::C_F;
   ME2*=cam;
 
   return ME2;
 
 }
 
 double jM2WunogqbarQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   return ME2;
 
 }
 
 double jM2WunogqbarQbar(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   return ME2;
 
 }
 
 double jM2Wunogqbarg(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   double ratio; // p2-/pb- in the notes
   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 cam = ( (HEJ::C_A - 1/HEJ::C_A)*(ratio + 1./ratio)/2. + 1/HEJ::C_A)/HEJ::C_F;
   ME2*=cam;
 
   return ME2;
 
 }
 
 
 // W+Jets qqxExtremal
 // W+Jets qqxExtremal Currents - wqq emission
 double jM2WgQtoqbarqQ(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqbarout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   //Correct colour averaging
   ME2*=(3.0/8.0);
 
   return ME2;
 
 }
 
 double jM2WgQtoqqbarQ(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqbarout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in){
 
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   //Correct colour averaging
   ME2*=(3.0/8.0);
 
   return ME2;
 
 }
 
 
 double jM2Wggtoqbarqg(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqbarout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
 {
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   double ratio; // p2-/pb- in the notes
   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 cam = ( (HEJ::C_A - 1/HEJ::C_A)*(ratio + 1./ratio)/2. + 1/HEJ::C_A)/HEJ::C_F;
   ME2*=cam;
 
   //Correct colour averaging
   ME2*=(3.0/8.0);
 
   return ME2;
 }
 
 double jM2Wggtoqqbarg(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqbarout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in){
 
   //COM temp;
   double ME2mpp=0.;
   double ME2mpm=0.;
   double ME2mmp=0.;
   double ME2mmm=0.;
   double ME2;
 
   ME2mpp = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,true,true);
   ME2mpm = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,true,false);
   ME2mmp = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,false,true);
   ME2mmm = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,false,false);
 
   //Helicity sum
   ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
 
   double ratio; // p2-/pb- in the notes
   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 cam = ( (HEJ::C_A - 1/HEJ::C_A)*(ratio + 1./ratio)/2. + 1/HEJ::C_A)/HEJ::C_F;
   ME2*=cam;
 
   //Correct colour averaging
   ME2*=(3.0/8.0);
 
   return ME2;
 
 }
 
 
 namespace {
 //Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
   Tensor<1,4> qggm1(CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool hel2, bool helg, CLHEP::HepLorentzVector refmom){
 
     double t1 = (p3-pb)*(p3-pb);
     Tensor<1,4> Tp3 = Construct1Tensor((p3));//p3
     Tensor<1,4> Tpb = Construct1Tensor((pb));//pb
     // Gauge choice in polarisation tensor. (see JC's Thesis)
     Tensor<1,4> epsg = eps(pb, refmom, helg);
     Tensor<3,4> qqCurBlank = T3Current(p2,hel2,p3,hel2);
     Tensor<2,4> qqCur = qqCurBlank.contract(Tp3-Tpb,2);
     Tensor<1,4> gqqCur = qqCur.contract(epsg,2)/t1;
 
     return gqqCur*(-1);
 }
 
 //Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
   Tensor<1,4> qggm2(CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool hel2, bool helg, CLHEP::HepLorentzVector refmom){
 
     double t1 = (p2-pb)*(p2-pb);
     Tensor<1,4> Tp2 = Construct1Tensor((p2));//p2
     Tensor<1,4> Tpb = Construct1Tensor((pb));//pb
     // Gauge choice in polarisation tensor. (see JC's Thesis)
     Tensor<1,4> epsg = eps(pb,refmom, helg);
     Tensor<3,4> qqCurBlank = T3Current(p2,hel2,p3,hel2);
     Tensor<2,4> qqCur = qqCurBlank.contract(Tp2-Tpb,2);
     Tensor<1,4> gqqCur = qqCur.contract(epsg,1)/t1;
 
     return gqqCur;
 }
 
 //Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
   Tensor<1,4> qggm3(CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool hel2, bool helg, CLHEP::HepLorentzVector refmom){
 
     double s23 = (p2+p3)*(p2+p3);
     Tensor<1,4> Tp2 = Construct1Tensor((p2));//p2
     Tensor<1,4> Tp3 = Construct1Tensor((p3));//p3
     Tensor<1,4> Tpb = Construct1Tensor((pb));//pb
     // Gauge choice in polarisation tensor. (see JC's Thesis)
     Tensor<1,4> epsg = eps(pb, refmom, helg);
     Tensor<2,4> g=Metric();
     Tensor<3,4> qqCurBlank1 = g.leftprod(Tp2+Tp3)/s23;
     Tensor<3,4> qqCurBlank2 = g.leftprod(Tpb)/s23;
     Tensor<1,4> Cur23 = TCurrent(p2,hel2, p3,hel2);
 
     Tensor<2,4> qqCur1 = qqCurBlank1.contract(Cur23,3);
     Tensor<2,4> qqCur2 = qqCurBlank2.contract(Cur23,3);
     Tensor<2,4> qqCur3 = qqCurBlank2.contract(Cur23,1);
 
     Tensor<1,4> gqqCur = (qqCur1.contract(epsg,1)
                           - qqCur2.contract(epsg,2)
                           + qqCur3.contract(epsg,1))*2*COM(0,1);
     return gqqCur;
 }
 }
 
 // no wqq emission
 double jM2WgqtoQQqW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1,  CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, bool aqlinepa){
 
   static bool is_sigma_index_set(false);
   if(!is_sigma_index_set){
     if(init_sigma_index())
       is_sigma_index_set = true;
     else
       return 0.;}
 
   // 2 independent helicity choices (complex conjugation related).
   Tensor<1,4> TMmmm1 = qggm1(pb,p2,p3,false,false, pa);
   Tensor<1,4> TMmmm2 = qggm2(pb,p2,p3,false,false, pa);
   Tensor<1,4> TMmmm3 = qggm3(pb,p2,p3,false,false, pa);
   Tensor<1,4> TMpmm1 = qggm1(pb,p2,p3,false,true, pa);
   Tensor<1,4> TMpmm2 = qggm2(pb,p2,p3,false,true, pa);
   Tensor<1,4> TMpmm3 = qggm3(pb,p2,p3,false,true, pa);
 
   // Build the external quark line W Emmision
   Tensor<1,4> cur1a = jW(pa,p1,plbar,pl, aqlinepa);
 
   //Contract with the qqxCurrent.
   COM Mmmm1 = TMmmm1.contract(cur1a,1).at(0);
   COM Mmmm2 = TMmmm2.contract(cur1a,1).at(0);
   COM Mmmm3 = TMmmm3.contract(cur1a,1).at(0);
   COM Mpmm1 = TMpmm1.contract(cur1a,1).at(0);
   COM Mpmm2 = TMpmm2.contract(cur1a,1).at(0);
   COM Mpmm3 = TMpmm3.contract(cur1a,1).at(0);
 
   //Colour factors:
   COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
   cm1m1=8./3.;
   cm2m2=8./3.;
   cm3m3=6.;
   cm1m2 =-1./3.;
   cm1m3 = -3.*COM(0.,1.);
   cm2m3 = 3.*COM(0.,1.);
 
   //Sqaure and sum for each helicity config:
   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)));
   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)));
 
   // Divide by WProp
   double WPropfact = WProp(plbar, pl);
 
   return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
 }
 
 
 // W+Jets qqxCentral
 double jM2WqqtoqQQq(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove)
 {
 
   static bool is_sigma_index_set(false);
   if(!is_sigma_index_set){
     if(init_sigma_index())
       is_sigma_index_set = true;
     else
       return 0.;}
 
   HLV pq, pqbar, p1, p4;
   if (qqxmarker){
     pqbar = partons[nabove+1];
     pq = partons[nabove+2];}
   else{
     pq = partons[nabove+1];
     pqbar = partons[nabove+2];}
 
   p1 = partons.front();
   p4 = partons.back();
 
   Tensor<1,4> T1am, T4bm, T1ap, T4bp;
   if(!(aqlinepa)){
     T1ap = TCurrent(p1, true, pa, true);
     T1am = TCurrent(p1, false, pa, false);}
   else if(aqlinepa){
     T1ap = TCurrent(pa, true, p1, true);
     T1am = TCurrent(pa, false, p1, false);}
   if(!(aqlinepb)){
     T4bp = TCurrent(p4, true, pb, true);
     T4bm = TCurrent(p4, false, pb, false);}
   else if(aqlinepb){
     T4bp = TCurrent(pb, true, p4, true);
     T4bm = TCurrent(pb, false, p4, false);}
 
   // Calculate the 3 separate contributions to the effective vertex
   Tensor<2,4> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
   Tensor<2,4> Xcro = MCrossW(  pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
   Tensor<2,4> Xsym = MSymW(    pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
 
   // 4 Different Helicity Choices (Differs from Pure Jet Case, where there is also the choice in qqbar helicity.
   // (- - hel choice)
   COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1)).at(0);
   COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1)).at(0);
   COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1)).at(0);
   // (- + hel choice)
   COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1)).at(0);
   COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1)).at(0);
   COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1)).at(0);
   // (+ - hel choice)
   COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1)).at(0);
   COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1)).at(0);
   COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1)).at(0);
   // (+ + hel choice)
   COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1)).at(0);
   COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1)).at(0);
   COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1)).at(0);
 
   //Colour factors:
   COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
   cmsms=3.;
   cmumu=4./3.;
   cmcmc=4./3.;
   cmsmu =3./2.*COM(0.,1.);
   cmsmc = -3./2.*COM(0.,1.);
   cmumc = -1./6.;
 
   // Work Out Interference in each case of helicity:
   double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
              +cmumu*pow(abs(M_mmUnc),2)
              +cmcmc*pow(abs(M_mmCro),2)
              +2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
              +2.*real(cmsmc*M_mmSym*conj(M_mmCro))
              +2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
 
   double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
              +cmumu*pow(abs(M_mpUnc),2)
              +cmcmc*pow(abs(M_mpCro),2)
              +2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
              +2.*real(cmsmc*M_mpSym*conj(M_mpCro))
              +2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
 
   double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
              +cmumu*pow(abs(M_pmUnc),2)
              +cmcmc*pow(abs(M_pmCro),2)
              +2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
              +2.*real(cmsmc*M_pmSym*conj(M_pmCro))
              +2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
 
   double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
              +cmumu*pow(abs(M_ppUnc),2)
              +cmcmc*pow(abs(M_ppCro),2)
              +2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
              +2.*real(cmsmc*M_ppSym*conj(M_ppCro))
              +2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
 
   double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
 
 
 
   CLHEP::HepLorentzVector q1,q3;
   q1=pa;
   for(int i=0;i<nabove+1;i++){
     q1-=partons.at(i);
   }
   q3 = q1 - pq - pqbar - pl - plbar;
 
   double t1 = (q1).m2();
   double t3 = (q3).m2();
 
   //Divide by t-channels
   amp/=(t1*t1*t3*t3);
 
   //Divide by WProp
   double WPropfact = WProp(plbar, pl);
   amp*=WPropfact;
 
 
   return amp;
 }
 
 // no wqq emission
 double jM2WqqtoqQQqW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<CLHEP::HepLorentzVector> partons, bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove, int nbelow, bool forwards){
 
   static bool is_sigma_index_set(false);
   if(!is_sigma_index_set){
     if(init_sigma_index())
       is_sigma_index_set = true;
     else
       return 0.;
   }
 
   if (!forwards){ //If Emission from Leg a instead, flip process.
     HLV dummymom = pa;
     bool dummybool= aqlinepa;
     int dummyint = nabove;
     pa = pb;
     pb = dummymom;
     std::reverse(partons.begin(),partons.end());
     qqxmarker = !(qqxmarker);
     aqlinepa = aqlinepb;
     aqlinepb = dummybool;
     nabove = nbelow;
     nbelow = dummyint;
   }
 
   HLV pq, pqbar, p1,p4;
   if (qqxmarker){
     pqbar = partons[nabove+1];
     pq = partons[nabove+2];}
   else{
     pq = partons[nabove+1];
     pqbar = partons[nabove+2];}
 
   p1 = partons.front();
   p4 = partons.back();
 
   Tensor<1,4> T1am(0.), T1ap(0.);
   if(!(aqlinepa)){
     T1ap = TCurrent(p1, true, pa, true);
     T1am = TCurrent(p1, false, pa, false);}
   else if(aqlinepa){
     T1ap = TCurrent(pa, true, p1, true);
     T1am = TCurrent(pa, false, p1, false);}
 
   Tensor <1,4> T4bm = jW(pb, p4, plbar, pl, aqlinepb);
 
   // Calculate the 3 separate contributions to the effective vertex
   Tensor<2,4> Xunc_m = MUncross(pa, pq, pqbar,partons, false, nabove);
   Tensor<2,4> Xcro_m = MCross(  pa, pq, pqbar,partons, false, nabove);
   Tensor<2,4> Xsym_m = MSym(    pa, p1, pb, p4, pq, pqbar, partons, false, nabove);
 
   Tensor<2,4> Xunc_p = MUncross(pa, pq, pqbar,partons, true, nabove);
   Tensor<2,4> Xcro_p = MCross(  pa, pq, pqbar,partons, true, nabove);
   Tensor<2,4> Xsym_p = MSym(    pa, p1, pb, p4, pq, pqbar, partons, true, nabove);
 
 
   // (- - hel choice)
   COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1)).at(0);
   COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1)).at(0);
   COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1)).at(0);
   // (- + hel choice)
   COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1)).at(0);
   COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1)).at(0);
   COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1)).at(0);
   // (+ - hel choice)
   COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
   COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
   COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
   // (+ + hel choice)
   COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
   COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
   COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
 
   //Colour factors:
   COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
   cmsms=3.;
   cmumu=4./3.;
   cmcmc=4./3.;
   cmsmu =3./2.*COM(0.,1.);
   cmsmc = -3./2.*COM(0.,1.);
   cmumc = -1./6.;
 
   // Work Out Interference in each case of helicity:
   double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
              +cmumu*pow(abs(M_mmUnc),2)
              +cmcmc*pow(abs(M_mmCro),2)
              +2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
              +2.*real(cmsmc*M_mmSym*conj(M_mmCro))
              +2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
 
   double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
              +cmumu*pow(abs(M_mpUnc),2)
              +cmcmc*pow(abs(M_mpCro),2)
              +2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
              +2.*real(cmsmc*M_mpSym*conj(M_mpCro))
              +2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
 
   double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
              +cmumu*pow(abs(M_pmUnc),2)
              +cmcmc*pow(abs(M_pmCro),2)
              +2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
              +2.*real(cmsmc*M_pmSym*conj(M_pmCro))
              +2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
 
   double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
              +cmumu*pow(abs(M_ppUnc),2)
              +cmcmc*pow(abs(M_ppCro),2)
              +2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
              +2.*real(cmsmc*M_ppSym*conj(M_ppCro))
              +2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
 
   double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
 
   CLHEP::HepLorentzVector q1,q3;
   q1=pa;
   for(int i=0;i<nabove+1;i++){
     q1-=partons.at(i);
   }
   q3 = q1 - pq - pqbar;
 
   double t1 = (q1).m2();
   double t3 = (q3).m2();
 
   //Divide by t-channels
   amp/=(t1*t1*t3*t3);
 
   //Divide by WProp
   double WPropfact = WProp(plbar, pl);
   amp*=WPropfact;
 
   return amp;
 }