diff --git a/MatrixElement/EW/ElectroWeakReweighter.cc b/MatrixElement/EW/ElectroWeakReweighter.cc
--- a/MatrixElement/EW/ElectroWeakReweighter.cc
+++ b/MatrixElement/EW/ElectroWeakReweighter.cc
@@ -1,1985 +1,1989 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the ElectroWeakReweighter class.
 //
 
 #include "ElectroWeakReweighter.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/Interface/Reference.h"
 #include "ThePEG/Interface/Switch.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/Repository/UseRandom.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 #include "boost/numeric/ublas/matrix.hpp"
 #include "boost/numeric/ublas/operation.hpp"
 #include "EWProcess.h"
 #include "HighEnergyMatching.h"
 #include "ElectroWeakMatching.h"
 #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
 #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
 #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
 #include "ThePEG/Helicity/epsilon.h"
 #include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
 #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
 #include "ThePEG/Handlers/StandardXComb.h"
 
 using namespace Herwig;
 
 tEWCouplingsPtr ElectroWeakReweighter::staticEWCouplings_ = tEWCouplingsPtr();
 
 
 ElectroWeakReweighter::ElectroWeakReweighter() : testing_(false)
 {}
 
 ElectroWeakReweighter::~ElectroWeakReweighter() {}
 
 IBPtr ElectroWeakReweighter::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr ElectroWeakReweighter::fullclone() const {
   return new_ptr(*this);
 }
 
 void ElectroWeakReweighter::persistentOutput(PersistentOStream & os) const {
   os << EWCouplings_ << collinearSudakov_ << softSudakov_ << testing_;
 }
 
 void ElectroWeakReweighter::persistentInput(PersistentIStream & is, int) {
   is >> EWCouplings_ >> collinearSudakov_ >> softSudakov_ >> testing_;
 }
 
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<ElectroWeakReweighter,ReweightBase>
 describeHerwigElectroWeakReweighter("Herwig::ElectroWeakReweighter", "HwMEEW.so");
 
 void ElectroWeakReweighter::Init() {
 
   static ClassDocumentation<ElectroWeakReweighter> documentation
     ("There is no documentation for the ElectroWeakReweighter class");
 
   static Reference<ElectroWeakReweighter,EWCouplings> interfaceEWCouplings
     ("EWCouplings",
      "The object to calculate the electroweak couplings",
      &ElectroWeakReweighter::EWCouplings_, false, false, true, false, false);
 
   static Reference<ElectroWeakReweighter,CollinearSudakov> interfaceCollinearSudakov
     ("CollinearSudakov",
      "The collinear Sudakov",
      &ElectroWeakReweighter::collinearSudakov_, false, false, true, false, false);
 
   static Reference<ElectroWeakReweighter,SoftSudakov> interfaceSoftSudakov
     ("SoftSudakov",
      "The soft Sudakov",
      &ElectroWeakReweighter::softSudakov_, false, false, true, false, false);
 
   static Switch<ElectroWeakReweighter,bool> interfaceTesting
     ("Testing",
      "Whether or not to output testing information",
      &ElectroWeakReweighter::testing_, false, false, false);
   static SwitchOption interfaceTestingYes
     (interfaceTesting,
      "Yes",
      "Output the information",
      true);
   static SwitchOption interfaceTestingNo
     (interfaceTesting,
      "No",
      "Don't output the information",
      false);
 
 }
 
 void ElectroWeakReweighter::doinit() {
   ReweightBase::doinit();
   if(!testing_) return;
   // testing output
   cerr << "aEM\n";
   for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.1) {
     cerr << scale/GeV << " " 
   	 << EWCouplings_->aEM(scale) << "\n";
   }
   cerr << "aS\n";
   for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
     cerr << scale/GeV << " " 
   	 << EWCouplings_->aS(scale) << "\n";
   }
   cerr << "y_t\n";
   for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
     cerr << scale/GeV << " " 
   	 << EWCouplings_->y_t(scale) << "\n";
   }
   cerr << "lambda\n";
   for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
     cerr << scale/GeV << " " 
   	 << EWCouplings_->lambda(scale) << "\n";
   }
   cerr << "vev\n";
   for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
     cerr << scale/GeV << " " 
   	 << EWCouplings_->vev(scale)/GeV << "\n";
   }
   collinearSudakov_->makePlots();
   Energy2 s = sqr(5000.*GeV);
   Energy2 t = -0.25*s;
   Energy2 u = -0.75*s;
   testEvolution(s,t,u);
 }
 
 namespace {
 #ifdef ThePEG_HAS_UNITS_CHECKING
 void axpy_prod_local(const boost::numeric::ublas::matrix<Complex>  & A,
 		     const boost::numeric::ublas::matrix<complex<InvEnergy2> > & B,
 		     boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
   assert(A.size2()==B.size1());
   C.resize(A.size1(),B.size2());
   for(unsigned int ix=0;ix<A.size1();++ix) {
     for(unsigned int iy=0;iy<B.size2();++iy) {
       C(ix,iy) = ZERO;
       for(unsigned int iz=0;iz<A.size2();++iz) {
 	C(ix,iy) += A(ix,iz)*B(iz,iy);
       }
     }
   }
 }
 
 void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> >  & A,
 		     const boost::numeric::ublas::vector<complex<Energy2> >      & B,
 		     boost::numeric::ublas::vector<Complex > & C) {
   assert(A.size2()==B.size());
   C.resize(A.size1());
   for(unsigned int ix=0;ix<A.size1();++ix) {
     C(ix) = ZERO;
       for(unsigned int iz=0;iz<A.size2();++iz) {
 	C(ix) += A(ix,iz)*B(iz);
       }
   }
 }
 
 void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> > & A,
 		     const boost::numeric::ublas::matrix<Complex> & B,
 		     boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
   assert(A.size2()==B.size1());
   C.resize(A.size1(),B.size2());
   for(unsigned int ix=0;ix<A.size1();++ix) {
     for(unsigned int iy=0;iy<B.size2();++iy) {
       C(ix,iy) = ZERO;
       for(unsigned int iz=0;iz<A.size2();++iz) {
 	C(ix,iy) += A(ix,iz)*B(iz,iy);
       }
     }
   }
 }
 
 #else
 void axpy_prod_local(const boost::numeric::ublas::matrix<Complex> & A,
 		     const boost::numeric::ublas::matrix<Complex> & B,
 		     boost::numeric::ublas::matrix<Complex> & C) {
+  assert(A.size2()==B.size1());
+  C.resize(A.size1(),B.size2());
   axpy_prod(A,B,C);
 }
 
 void axpy_prod_local(const boost::numeric::ublas::matrix<Complex> & A,
 		     const boost::numeric::ublas::vector<Complex> & B,
 		     boost::numeric::ublas::vector<Complex> & C) {
+  assert(A.size2()==B.size());
+  C.resize(A.size1());
   axpy_prod(A,B,C);
 }
 #endif
 
 }
 
 double ElectroWeakReweighter::weight() const {
   EWCouplings_->initialize();
   staticEWCouplings_ = EWCouplings_;
   // cast the XComb
   Ptr<StandardXComb>::ptr sxc = dynamic_ptr_cast<Ptr<StandardXComb>::ptr>(lastXCombPtr());
   // if the Herwig XComb
   if(sxc) {
     // get information about the type of event
     Ptr<SubtractedME>::tptr      subme = dynamic_ptr_cast<Ptr<SubtractedME>::tptr>(sxc->matrixElement());
     Ptr<MatchboxMEBase>::tptr       me = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(sxc->matrixElement());
     Ptr<SubtractionDipole>::tptr dipme = dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(sxc->matrixElement());
     bool isHEvent(false),isSEvent(false);
     if(subme) {
       if ( subme->realShowerSubtraction() )
 	isHEvent = true;
       else if ( subme->virtualShowerSubtraction() ||  subme->loopSimSubtraction() )
 	isSEvent = true;
     }
     // H or S event of virtual return 1.
     if(isHEvent || isSEvent || (me && me->oneLoopNoBorn()))
       return 1.;
     // cerr << "testing after type check\n";
     // cerr << "testing pointers " << subme << " " << me << " " << dipme << "\n";
     // cerr << "testing event type " << isHEvent << " " << isSEvent << " " << "\n";
     // if(subme) cerr << subme->fullName() << "\n";
     // if(   me) {
     //   cerr <<    me->fullName() << "\n";
     //   cerr << me->oneLoopNoBorn() << " " << me->oneLoopNoLoops() << "\n";
     // }
     // if(dipme) cerr << dipme->fullName() << "\n";
   }
   // cerr <<  subProcess() << "\n";
   // cerr << *subProcess() << "\n";
   // only 2->2 processes
   if(subProcess()->outgoing().size()!=2) return 1.;
   // processes with gg initial-state
   if(subProcess()->incoming().first->id()==ParticleID::g &&
      subProcess()->incoming().second->id()==ParticleID::g) {
     if(subProcess()->outgoing()[0]->id()==ParticleID::g &&
        subProcess()->outgoing()[1]->id()==ParticleID::g)
       return 1.;
     else if(abs(subProcess()->outgoing()[0]->id())<=6 &&
 	    subProcess()->outgoing()[0]->id()==-subProcess()->outgoing()[1]->id()) {
       return reweightggqqbar();
     }
     else
       assert(false);
   }
   // processes with q qbar initial-state
   else if((subProcess()->incoming().first ->id() > 0 &&
 	   subProcess()->incoming().first ->id()<= 5 &&
 	   subProcess()->incoming().second->id() < 0 &&
 	   subProcess()->incoming().second->id()>=-5) ||
 	  (subProcess()->incoming().second->id() > 0 &&
 	   subProcess()->incoming().second->id()<= 5 &&
 	   subProcess()->incoming().first ->id() < 0 &&
 	   subProcess()->incoming().first ->id()>=-5)) {
     // identical flavour q qbar
     if(subProcess()->incoming().first ->id() == -subProcess()->incoming().second->id()) {
       // q qbar -> gg
       if(subProcess()->outgoing()[0]->id()==ParticleID::g &&
 	 subProcess()->outgoing()[1]->id()==ParticleID::g)
 	return reweightqqbargg();
       // q qbar -> q' q'bar
       else if(subProcess()->outgoing()[0]->id() == -subProcess()->outgoing()[1]->id() &&
 	      abs(subProcess()->outgoing()[0]->id())<=6)
 	return reweightqqbarqqbarS();
     }
     // different flavour q qbar
     else {
       if((subProcess()->outgoing()[0]->id() > 0 &&
 	  subProcess()->outgoing()[0]->id()<= 5 &&
 	  subProcess()->outgoing()[1]->id() < 0 &&
 	  subProcess()->outgoing()[1]->id()>=-5) ||
 	 (subProcess()->outgoing()[1]->id() > 0 &&
 	  subProcess()->outgoing()[1]->id()<= 5 &&
 	  subProcess()->outgoing()[0]->id() < 0 &&
 	  subProcess()->outgoing()[0]->id()>=-5)) {
 	return reweightqqbarqqbarT();
       }
       else
 	assert(false);
     }
   }
   // processes with q g initial-state
   else if((subProcess()->incoming().first ->id()> 0 &&
 	   subProcess()->incoming().first ->id()<=5 &&
 	   subProcess()->incoming().second->id()==ParticleID::g) ||
 	  (subProcess()->incoming().second->id()> 0 &&
 	   subProcess()->incoming().second->id()<=5 &&
 	   subProcess()->incoming().first ->id()==ParticleID::g)) {
     // qg -> qg
     if((subProcess()->outgoing()[0]->id()> 0 &&
 	subProcess()->outgoing()[0]->id()<=5 &&
 	subProcess()->outgoing()[1]->id()==ParticleID::g) ||
        (subProcess()->outgoing()[1]->id()> 0 &&
 	subProcess()->outgoing()[1]->id()<=5 &&
 	subProcess()->outgoing()[0]->id()==ParticleID::g))
       return reweightqgqg();
     // unknown
     else
       assert(false);
   }
   // processes with qbar g initial-state
   else if((subProcess()->incoming().first ->id()>=-5 &&
 	   subProcess()->incoming().first ->id()<  0 &&
 	   subProcess()->incoming().second->id()==ParticleID::g) ||
 	  (subProcess()->incoming().second->id()>=-5 &&
 	   subProcess()->incoming().second->id()<  0 &&
 	   subProcess()->incoming().first ->id()==ParticleID::g)) {
     if((subProcess()->outgoing()[0]->id()>=-5 &&
 	subProcess()->outgoing()[0]->id()<  0 &&
 	subProcess()->outgoing()[1]->id()==ParticleID::g) ||
        (subProcess()->outgoing()[1]->id()>=-5 &&
 	subProcess()->outgoing()[1]->id()<  0 &&
 	subProcess()->outgoing()[0]->id()==ParticleID::g))
       return reweightqbargqbarg();
     else
       assert(false);
   }
   // processes with q q initial-state
   else if( subProcess()->incoming().first ->id()> 0 &&
 	   subProcess()->incoming().first ->id()<=5 &&
 	   subProcess()->incoming().second->id()> 0 &&
 	   subProcess()->incoming().second->id()<=5 ) {
     if(subProcess()->outgoing()[0]->id()> 0 &&
        subProcess()->outgoing()[0]->id()<=5 &&
        subProcess()->outgoing()[1]->id()> 0 &&
        subProcess()->outgoing()[1]->id()<=5)
       return reweightqqqq();
     else
       assert(false);
   }
   // processes with qbar qbar initial-state
   else if( subProcess()->incoming().first ->id()<   0 &&
 	   subProcess()->incoming().first ->id()>= -5 &&
 	   subProcess()->incoming().second->id()<   0 &&
 	   subProcess()->incoming().second->id()>= -5 ) {
     if(subProcess()->outgoing()[0]->id()<   0 &&
        subProcess()->outgoing()[0]->id()>= -5 &&
        subProcess()->outgoing()[1]->id()<   0 &&
        subProcess()->outgoing()[1]->id()>= -5)
       return reweightqbarqbarqbarqbar();
     else
       assert(false);
   }
   // unknown initial-state
   else
     assert(false);
   assert(false);
   staticEWCouplings_ = tEWCouplingsPtr();
 }
 
 void ElectroWeakReweighter::testEvolution(Energy2 s,Energy2 t, Energy2 u) const {
   Energy highScale = sqrt(s);
   Energy ewScale = coupling()->mZ();
   Energy lowScale = 50.0*GeV;
   for (unsigned int i=0; i<45;++i) {
     EWProcess::Process process = (EWProcess::Process)i;
     cerr << "process " << process << "\n";
     // result for all EW and QCD SCET contributions:
     boost::numeric::ublas::matrix<complex<InvEnergy2> > highMatch_val 
       = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,true,true);
     boost::numeric::ublas::matrix<Complex> highRunning_val
       = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process,0);
     boost::numeric::ublas::matrix<Complex> ewMatch_val = 
       ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true,0);
     boost::numeric::ublas::matrix<Complex> lowRunning_val = 
       softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process,0);
     boost::numeric::ublas::matrix<Complex> collinearHighRunning_val =
       collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false);
     boost::numeric::ublas::matrix<Complex> collinearEWMatch_val =
       collinearSudakov_->electroWeakMatching(ewScale,s,process,true);
     boost::numeric::ublas::matrix<Complex> collinearLowRunning_val =
       collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process);
     boost::numeric::ublas::matrix<Complex> lowMatchTemp_val = 
       boost::numeric::ublas::zero_matrix<Complex>(ewMatch_val.size1(),ewMatch_val.size2());
     for (unsigned int ii=0; ii<ewMatch_val.size1(); ++ii) {
       for (unsigned int jj=0; jj<ewMatch_val.size2(); ++jj) {
 	lowMatchTemp_val(ii,jj) = collinearEWMatch_val(ii,jj)*ewMatch_val(ii,jj);
       }
     }
     boost::numeric::ublas::matrix<Complex> temp(highRunning_val.size1(),collinearHighRunning_val.size2());
     boost::numeric::ublas::axpy_prod(highRunning_val,collinearHighRunning_val,temp);
     boost::numeric::ublas::matrix<Complex> temp2(collinearLowRunning_val.size1(),lowRunning_val.size2());
     boost::numeric::ublas::axpy_prod(collinearLowRunning_val,lowRunning_val,temp2);
     boost::numeric::ublas::matrix<Complex> temp3(temp2.size1(),lowMatchTemp_val.size2());
     boost::numeric::ublas::axpy_prod(temp2,lowMatchTemp_val,temp3);
     temp2.resize(temp3.size1(),temp.size2());
     boost::numeric::ublas::axpy_prod(temp3,temp,temp2);
     boost::numeric::ublas::matrix<complex<InvEnergy2> >  result(temp2.size1(),highMatch_val.size2());
     axpy_prod_local(temp2,highMatch_val,result);
     for(unsigned int ix=0;ix<result.size1();++ix) {
       for(unsigned int iy=0;iy<result.size2();++iy) {
 	cerr << s*result(ix,iy) << " ";
       }
       cerr << "\n";
     }
   }
 }
 
 namespace {
 
 void SackGluonPolarizations(Lorentz5Momentum &p1,
 			    Lorentz5Momentum &p2,
 			    Lorentz5Momentum &p3,
 			    Lorentz5Momentum &p4,
 			    Energy2 s, Energy2 t, Energy2 u, Energy2 m2,
 			    vector<LorentzVector<Complex> > & eps3,
 			    vector<LorentzVector<Complex> > & eps4,
 			    unsigned int iopt) {
   static const Complex I(0.,1.);
   // p1 is p-, p2 is p+
   // p3 is k-, p4 is k+
   // both final-state
   if(iopt==0) {
     // swap t and u due Aneesh's defn
     Energy3 den1 = sqrt((u*t-sqr(m2))*(s-4.*m2));
     Energy3 den2 = sqrt(s*(u*t-sqr(m2)));
     LorentzVector<Complex> eps3Para = (m2+t)/den1*p1 -(m2+u)/den1*p2 +(u-t)/den1*p3;
     LorentzVector<Complex> eps3Perp = 2./den2*epsilon(p1,p2,p3);
     LorentzVector<Complex> eps4Para = (m2+t)/den1*p2 -(m2+u)/den1*p1 +(u-t)/den1*p4;
     LorentzVector<Complex> eps4Perp = 2./den2*epsilon(p1,p2,p4);
     eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
     eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
     eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
     eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
     if(m2!=ZERO) assert(false);
   }
   // both initial-state
   else if(iopt==1) {
     if(m2!=ZERO) assert(false);
     LorentzVector<Complex> eps3Para( 1., 0.,0.,0.);
     LorentzVector<Complex> eps3Perp( 0.,-1.,0.,0.);
     LorentzVector<Complex> eps4Para(-1.,0.,0., 0.);
     LorentzVector<Complex> eps4Perp( 0., 1.,0.,0.);
     eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
     eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
     eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
     eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
   }
   else if(iopt==2) {
     // rotation into the 2,3 Breit frame
     Lorentz5Momentum pa = p3-p2;
     Axis axis(pa.vect().unit());
     LorentzRotation rot;
     double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
     if ( sinth > 1.e-9 )
       rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
     rot.rotateX(Constants::pi);
     rot.boostZ( pa.e()/pa.vect().mag());
     Lorentz5Momentum ptemp=rot*p2;
     Boost trans = -1./ptemp.e()*ptemp.vect();
     trans.setZ(0.);
     rot.boost(trans);
     LorentzVector<Complex> eps3Para( 1., 0.,0.,0.);
     LorentzVector<Complex> eps3Perp( 0.,-1.,0.,0.);
     LorentzVector<Complex> eps4Para(-1.,0.,0., 0.);
     LorentzVector<Complex> eps4Perp( 0., 1.,0.,0.);
     eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
     eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
     eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
     eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
     rot = rot.invert();
     for(unsigned int ix=0;ix<2;++ix) {
       eps3[ix] *=rot;
       eps4[ix] *=rot;
     }
   }
   else if(iopt==3) {
     // rotation into the 1,4 Breit frame
     Lorentz5Momentum pa = p4-p1;
     Axis axis(pa.vect().unit());
     LorentzRotation rot;
     double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
     if ( sinth > 1.e-9 )
       rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
     rot.rotateX(Constants::pi);
     rot.boostZ( pa.e()/pa.vect().mag());
     Lorentz5Momentum ptemp=rot*p1;
     Boost trans = -1./ptemp.e()*ptemp.vect();
     trans.setZ(0.);
     rot.boost(trans);
     LorentzVector<Complex> eps3Para( 1., 0.,0.,0.);
     LorentzVector<Complex> eps3Perp( 0.,-1.,0.,0.);
     LorentzVector<Complex> eps4Para(-1.,0.,0., 0.);
     LorentzVector<Complex> eps4Perp( 0., 1.,0.,0.);
     eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
     eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
     eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
     eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
     rot = rot.invert();
     for(unsigned int ix=0;ix<2;++ix) {
       eps3[ix] *=rot;
       eps4[ix] *=rot;
     }
   }
   else
     assert(false);
 }
 
 }
 
 double ElectroWeakReweighter::reweightqqbargg() const {
   // momenta and invariants
   Lorentz5Momentum p1   = subProcess()->incoming().first ->momentum();
   tcPDPtr          q    = subProcess()->incoming().first ->dataPtr();
   Lorentz5Momentum p2   = subProcess()->incoming().second->momentum();
   tcPDPtr          qbar = subProcess()->incoming().second->dataPtr();
   if(subProcess()->incoming().first->id()<0) {
     swap(p1,p2  );
     swap(q ,qbar);
   }
   Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
   Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
   tcPDPtr           g = subProcess()->outgoing()[1]->dataPtr();
   Energy2 s = (p1+p2).m2();
   Energy2 t = (p1-p4).m2();
   Energy2 u = (p1-p3).m2();
   // boost to partonci rest frame
   Lorentz5Momentum psum=p1+p2;
   LorentzRotation boost(-psum.boostVector());
   p1 *= boost;
   p2 *= boost;
   p3 *= boost;
   p4 *= boost;
   // LO and EW corrected matrix element coefficients
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights;
   // quark left doublet
   if(q->id()!=5) {
     bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true ,0);
     EWQQGGweights   = evaluateRunning(EWProcess::QQGG,s,t,u,false,0);
   }
   else {
     bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true ,0);
     EWQQGGweights   = evaluateRunning(EWProcess::QtQtGG,s,t,u,false,0);
   }
   // quark right singlet
   if(abs(subProcess()->incoming().first->id())%2==0) {
     bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true ,0);
     EWRRGGweights   = evaluateRunning(EWProcess::UUGG,s,t,u,false,0);
   }
   else {
     bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,0);
     EWRRGGweights   = evaluateRunning(EWProcess::DDGG,s,t,u,false,0);
   }
   SpinorWaveFunction       qw(p1,q   ,incoming);
   SpinorBarWaveFunction qbarw(p2,qbar,incoming);
   vector<LorentzVector<Complex> > eps3,eps4;
   SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps3,eps4,0);
   boost::numeric::ublas::matrix<Complex>
     bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
     EWME   = boost::numeric::ublas::zero_matrix<Complex>(3,3);
   for(unsigned int iq=0;iq<2;++iq) {
     if(iq==0) {
       qw.reset   (0);
       qbarw.reset(1);
     }
     else {
       qw.reset   (1);
       qbarw.reset(0);
     }
     LorentzVector<complex<Energy> > current  = iq==0 ?
       qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
       qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
     for(unsigned int i1=0;i1<2;++i1) {
       complex<Energy> d31 = eps3[i1].dot(p1);
       for(unsigned int i2=0;i2<2;++i2) {
 	// g1w.reset(2*i1);
 	// g2w.reset(2*i2);
 	boost::numeric::ublas::vector<complex<Energy2> > M(5);
 	Complex         d34 = eps3[i1].dot(eps4[i2]);
 	complex<Energy> d42 = eps4[i2].dot(p2);
 	// M0 in paper
 	M(0) = qw.dimensionedWave().slash(eps3[i1])
 	  .slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps4[i2]);
 	// M4 in paper
 	M(2) =  current.dot(eps4[i2])*d31;
 	// M5 in paper
 	M(3) = -current.dot(eps3[i1])*d42;
 	// M1 in paper (missing factor)
 	M(1) =  current.dot(p4);
 	// M6 in paper
 	M(4) = M(1)*d31*d42/GeV2;
 	// M1 final factor
 	M(1) *= d34;
 	// coefficient of different contributions
 	boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3),Ctest(3);
 
 
 	// Ctest(0) = 1./6.*( MEU+MET);
 	// Ctest(1) = 0.5*( MEU+MET);
 	// Ctest(2) = 0.5*(MEU+MES-MET+MES);
 	if(iq==0) {
 	  axpy_prod_local(bornQQGGweights,M,Cborn);
 	  axpy_prod_local(EWQQGGweights  ,M,CEW  );
 	}
 	else {
 	  axpy_prod_local(bornRRGGweights,M,Cborn);
 	  axpy_prod_local(EWRRGGweights  ,M,CEW  );
 	}
 	unsigned int ioff = (Cborn.size()==6 && q->id()%2!=0) ? 3 : 0;
 	for(unsigned int ix=0;ix<3;++ix) {
 	  for(unsigned int iy=0;iy<3;++iy) {
 	    bornME(ix,iy) += Cborn(ix+ioff)*conj(Cborn(iy+ioff));
 	    EWME  (ix,iy) += CEW  (ix+ioff)*conj(CEW  (iy+ioff));
 	  }
 	}
       }
     }
   }
   double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
   double EW   = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
   return EW/born;
 }
 
 boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 ElectroWeakReweighter::evaluateRunning(EWProcess::Process process, Energy2 s,
 				       Energy2 t, Energy2 u, bool born,
 				       unsigned int iswap) const {
   using namespace boost::numeric::ublas;
   bool SU3save = coupling()->SU3();
   bool EWsave  = coupling()-> EW();
   Energy highScale = sqrt(s);
   Energy   ewScale = coupling()->mZ();
   Energy  lowScale = ewScale;
   // result for all EW and QCD SCET contributions:
   // MATCHING CONTRIBUTIONS
   // high energy matching
   matrix<complex<InvEnergy2> > highMatch_val;
   if(iswap==0) 
     highMatch_val = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,!born,false);
   else if(iswap==1)
     highMatch_val = HighEnergyMatching::highEnergyMatching(highScale,t,s,u,process,!born,false);
   else if(iswap==2)
     highMatch_val = HighEnergyMatching::highEnergyMatching(highScale,u,t,s,process,!born,false);
   else
     assert(false);
   // low energy matching
   matrix<Complex> 
     ewMatch_val = ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,!born,iswap);
   matrix<Complex> collinearEWMatch_val =
     collinearSudakov_->electroWeakMatching(ewScale,s,process,!born);
   // EVOLUTION
   matrix<Complex> highRunning_val,lowRunning_val,
     collinearHighRunning_val,collinearLowRunning_val;
   // born process
   if(born) {
     highRunning_val = identity_matrix<Complex>(softSudakov_->numberGauge(process));
     lowRunning_val  = identity_matrix<Complex>(softSudakov_->numberBrokenGauge(process));
     collinearHighRunning_val = identity_matrix<Complex>(softSudakov_->numberGauge(process));
     collinearLowRunning_val  = identity_matrix<Complex>(softSudakov_->numberBrokenGauge(process));
   }
   // EW corrected
   else {
     coupling()->SU3(false);
     coupling()-> EW( true);
     highRunning_val = softSudakov_->highEnergyRunning(highScale, ewScale,s,t,u,process,iswap);
     lowRunning_val  = softSudakov_->lowEnergyRunning (  ewScale,lowScale,s,t,u,process,iswap);
     collinearHighRunning_val = collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false);
     collinearLowRunning_val = collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process);
   };
   matrix<Complex> lowMatchTemp_val = 
     zero_matrix<Complex>(ewMatch_val.size1(),ewMatch_val.size2());
   for (unsigned int ii=0; ii<ewMatch_val.size1(); ++ii) {
     for (unsigned int jj=0; jj<ewMatch_val.size2(); ++jj) {
       lowMatchTemp_val(ii,jj) = collinearEWMatch_val(ii,jj)*ewMatch_val(ii,jj);
     }
   }
   // perform all the multiplications
   matrix<Complex> temp(highRunning_val.size1(),collinearHighRunning_val.size2());
   axpy_prod(highRunning_val,collinearHighRunning_val,temp);
   matrix<Complex> temp2(collinearLowRunning_val.size1(),lowRunning_val.size2());
   axpy_prod(collinearLowRunning_val,lowRunning_val,temp2);
   matrix<Complex> temp3(temp2.size1(),lowMatchTemp_val.size2());
   axpy_prod(temp2,lowMatchTemp_val,temp3);
   temp2.resize(temp3.size1(),temp.size2());
   axpy_prod(temp3,temp,temp2);
   matrix<complex<InvEnergy2> >  result(temp2.size1(),highMatch_val.size2());
   axpy_prod_local(temp2,highMatch_val,result);
   // reset the couplings
   coupling()->SU3(SU3save);
   coupling()-> EW( EWsave);
   // return the answer
   return result;
 }
 
 
 double ElectroWeakReweighter::reweightggqqbar() const {
   // momenta and invariants
   Lorentz5Momentum p1   = subProcess()->incoming().first ->momentum();
   Lorentz5Momentum p2   = subProcess()->incoming().second->momentum();
   Lorentz5Momentum p3   = subProcess()->outgoing()[0]->momentum();
   Lorentz5Momentum p4   = subProcess()->outgoing()[1]->momentum();
   tcPDPtr          qbar = subProcess()->outgoing()[0]->dataPtr();
   tcPDPtr          q    = subProcess()->outgoing()[1]->dataPtr();
   if(q->id()<0) {
     swap(p3,p4  );
     swap(q ,qbar);
   }
   Energy2 s = (p1+p2).m2();
   Energy2 t = (p1-p4).m2();
   Energy2 u = (p1-p3).m2();
   // boost to partonic rest frame and rescale momenta of outgoing
   // so zero mass
   Lorentz5Momentum psum=p1+p2;
   LorentzRotation boost(-psum.boostVector());
   p1 *= boost;
   p2 *= boost;
   p3 *= boost;
   p4 *= boost;
   p3.setMass(ZERO);
   p3.rescaleRho();
   p4.setMass(ZERO);
   p4.rescaleRho();
   // LO and EW matrix element coefficents
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights;
   // quark left doublet
   if(q->id()<5) {
     bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true ,0);
     EWQQGGweights   = evaluateRunning(EWProcess::QQGG,s,t,u,false,0);
   }
   else {
     bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true ,0);
     EWQQGGweights   = evaluateRunning(EWProcess::QtQtGG,s,t,u,false,0);
   }
   // quark right singlet
   if(q->id()==0) {
     if(q->id()==6) {
       bornRRGGweights = evaluateRunning(EWProcess::tRtRGG,s,t,u,true ,0);
       EWRRGGweights   = evaluateRunning(EWProcess::tRtRGG,s,t,u,false,0);
     }
     else {
       bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true ,0);
       EWRRGGweights   = evaluateRunning(EWProcess::UUGG,s,t,u,false,0);
     }
   }
   else {
     bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,0);
     EWRRGGweights   = evaluateRunning(EWProcess::DDGG,s,t,u,false,0);
   }
   SpinorWaveFunction       qw(p4,qbar,incoming);
   SpinorBarWaveFunction qbarw(p3,q   ,incoming);
   vector<LorentzVector<Complex> > eps1,eps2;
   SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps1,eps2,1);
   boost::numeric::ublas::matrix<Complex>
     bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
     EWME   = boost::numeric::ublas::zero_matrix<Complex>(3,3);
   // helicities of outgoing quarks
   for(unsigned int iq=0;iq<2;++iq) {
     if(iq==0) {
       qw.reset   (0);
       qbarw.reset(1);
     }
     else {
       qw.reset   (1);
       qbarw.reset(0);
     }
     LorentzVector<complex<Energy> > current  = iq==0 ?
       qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
       qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
     for(unsigned int i1=0;i1<2;++i1) {
       complex<Energy> d31 = eps1[i1].dot(p3);
       for(unsigned int i2=0;i2<2;++i2) {
   	boost::numeric::ublas::vector<complex<Energy2> > M(5);
   	Complex         d34 = eps1[i1].dot(eps2[i2]);
   	complex<Energy> d42 = eps2[i2].dot(p4);
   	// M0 in paper
   	M(0) = qw.dimensionedWave().slash(eps1[i1])
   	  .slash(p2-p4).vectorCurrent(qbarw.dimensionedWave()).dot(eps2[i2]);
   	// M4 in paper
   	M(2) =  current.dot(eps2[i2])*d31;
   	// M5 in paper
   	M(3) = -current.dot(eps1[i1])*d42;
   	// M1 in paper (missing factor)
   	M(1) =  current.dot(p2);
   	// M6 in paper
   	M(4) = M(1)*d31*d42/GeV2;
   	// M1 final factor
   	M(1) *= d34;
   	// coefficient of different contributions
  	boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3);
   	if(iq==0) {
   	  axpy_prod_local(bornQQGGweights,M,Cborn);
   	  axpy_prod_local(EWQQGGweights  ,M,CEW  );
   	}
   	else {
   	  axpy_prod_local(bornRRGGweights,M,Cborn);
   	  axpy_prod_local(EWRRGGweights  ,M,CEW  );
   	}
 	unsigned int ioff = (Cborn.size()==6 && q->id()%2!=0) ? 3 : 0;
   	for(unsigned int ix=0;ix<3;++ix) {
   	  for(unsigned int iy=0;iy<3;++iy) {
   	    bornME(ix,iy) += Cborn(ix+ioff)*conj(Cborn(iy+ioff));
   	    EWME  (ix,iy) += CEW  (ix+ioff)*conj(CEW  (iy+ioff));
   	  }
   	}
       }
     }
   }
   double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
   double EW   = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
   return EW/born;
 }
 
 double ElectroWeakReweighter::reweightqgqg() const {
   // momenta and invariants
   Lorentz5Momentum p1   = subProcess()->incoming().first ->momentum();
   Lorentz5Momentum p2   = subProcess()->incoming().second->momentum();
   tcPDPtr q;
   if(subProcess()->incoming().first->id()!=ParticleID::g) {
     q = subProcess()->incoming().first ->dataPtr();
   }
   else {
     q = subProcess()->incoming().second->dataPtr();
     swap(p1,p2);
   }
   Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
   Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
   if(subProcess()->outgoing()[0]->id()!=ParticleID::g)
     swap(p3,p4);
   Energy2 s = (p1+p2).m2();
   Energy2 t = (p1-p4).m2();
   Energy2 u = (p1-p3).m2();
   // boost to partonic rest frame
   Lorentz5Momentum psum=p1+p2;
   LorentzRotation boost(-psum.boostVector());
   p1 *= boost;
   p2 *= boost;
   p3 *= boost;
   p4 *= boost;
   // LO and EW corrected matrix element coefficients
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights;
   // quark left doublet
   if(q->id()!=5) {
     bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true ,1);
     EWQQGGweights   = evaluateRunning(EWProcess::QQGG,s,t,u,false,1);
   }
   else {
     bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true ,1);
     EWQQGGweights   = evaluateRunning(EWProcess::QtQtGG,s,t,u,false,1);
   }
   // quark right singlet
   if(abs(q->id())%2==0) {
     bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true ,1);
     EWRRGGweights   = evaluateRunning(EWProcess::UUGG,s,t,u,false,1);
   }
   else {
     bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,1);
     EWRRGGweights   = evaluateRunning(EWProcess::DDGG,s,t,u,false,1);
   }
   SpinorWaveFunction       qw(p1,q,incoming);
   SpinorBarWaveFunction qbarw(p4,q,outgoing);
   vector<LorentzVector<Complex> > eps2,eps3;
   SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps2,eps3,2);
   // compute the matrix elements
   boost::numeric::ublas::matrix<Complex>
     bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
     EWME   = boost::numeric::ublas::zero_matrix<Complex>(3,3),
     testME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
   for(unsigned int iq=0;iq<2;++iq) {
     if(iq==0) {
       qw.reset   (0);
       qbarw.reset(0);
     }
     else {
       qw.reset   (1);
       qbarw.reset(1);
     }
     LorentzVector<complex<Energy> > current  = iq==0 ?
       qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
       qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
     for(unsigned int i1=0;i1<2;++i1) {
       complex<Energy> d31 = eps3[i1].dot(p1);
       for(unsigned int i2=0;i2<2;++i2) {
   	boost::numeric::ublas::vector<complex<Energy2> > M(5);
    	Complex         d34 = eps3[i1].dot(eps2[i2]);
    	complex<Energy> d42 = eps2[i2].dot(p4);
    	// M0 in paper
   	M(0) = qw.dimensionedWave().slash(eps3[i1])
    	  .slash(p2-p4).vectorCurrent(qbarw.dimensionedWave()).dot(eps2[i2]);
   	// M4 in paper
   	M(2) =  current.dot(eps2[i2])*d31;
   	// M5 in paper
   	M(3) = -current.dot(eps3[i1])*d42;
   	// M1 in paper (missing factor)
   	M(1) =  current.dot(p2);
  	// M6 in paper
  	M(4) = M(1)*d31*d42/GeV2;
  	// M1 final factor
  	M(1) *= d34;
  	// coefficient of different contributions
  	boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3);
   	if(iq==0) {
   	  axpy_prod_local(bornQQGGweights,M,Cborn);
   	  axpy_prod_local(EWQQGGweights  ,M,CEW  );
   	}
   	else {
   	  axpy_prod_local(bornRRGGweights,M,Cborn);
   	  axpy_prod_local(EWRRGGweights  ,M,CEW  );
   	}
    	unsigned int ioff = (Cborn.size()==6 && q->id()%2!=0) ? 3 : 0;
   	for(unsigned int ix=0;ix<3;++ix) {
   	  for(unsigned int iy=0;iy<3;++iy) {
   	    bornME(ix,iy) += Cborn(ix+ioff)*conj(Cborn(iy+ioff));
   	    EWME  (ix,iy) += CEW  (ix+ioff)*conj(CEW  (iy+ioff));
   	  }
   	}
       }
     }
   }
   double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
   double EW   = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
   return EW/born;
 }
 
 double ElectroWeakReweighter::reweightqbargqbarg() const {
   // momenta and invariants
   Lorentz5Momentum p1   = subProcess()->incoming().first ->momentum();
   Lorentz5Momentum p2   = subProcess()->incoming().second->momentum();
   tcPDPtr          qbar;
   if(subProcess()->incoming().first->id()==ParticleID::g) {
     qbar = subProcess()->incoming().second->dataPtr();
   }
   else {
     qbar = subProcess()->incoming().first ->dataPtr();
     swap(p1,p2);
   }
   Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
   Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
   if(subProcess()->outgoing()[0]->id()==ParticleID::g)
     swap(p3,p4);
   Energy2 s = (p1+p2).m2();
   Energy2 t = (p1-p4).m2();
   Energy2 u = (p1-p3).m2();
   // boost to partonci rest frame
   Lorentz5Momentum psum=p1+p2;
   LorentzRotation boost(-psum.boostVector());
   p1 *= boost;
   p2 *= boost;
   p3 *= boost;
   p4 *= boost;
   // LO and EW corrected matrix element coefficients
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights;
   // quark left doublet
   if(qbar->id()!=-5) {
     bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true ,1);
     EWQQGGweights   = evaluateRunning(EWProcess::QQGG,s,t,u,false,1);
   }
   else {
     bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true ,1);
     EWQQGGweights   = evaluateRunning(EWProcess::QtQtGG,s,t,u,false,1);
   }
   // quark right singlet
   if(abs(qbar->id())%2==0) {
     bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true ,1);
     EWRRGGweights   = evaluateRunning(EWProcess::UUGG,s,t,u,false,1);
   }
   else {
     bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,1);
     EWRRGGweights   = evaluateRunning(EWProcess::DDGG,s,t,u,false,1);
   }
   SpinorWaveFunction       qw(p3,qbar,outgoing);
   SpinorBarWaveFunction qbarw(p2,qbar,incoming);
   vector<LorentzVector<Complex> > eps1,eps4;
   SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps1,eps4,3);
   boost::numeric::ublas::matrix<Complex>
     bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
     EWME   = boost::numeric::ublas::zero_matrix<Complex>(3,3);
   for(unsigned int iq=0;iq<2;++iq) {
     if(iq==0) {
       qw.reset   (1);
       qbarw.reset(1);
     }
     else {
       qw.reset   (0);
       qbarw.reset(0);
     }
     LorentzVector<complex<Energy> > current  = iq==0 ?
       qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
       qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
     for(unsigned int i1=0;i1<2;++i1) {
       complex<Energy> d31 = eps1[i1].dot(p3);
       for(unsigned int i2=0;i2<2;++i2) {
    	boost::numeric::ublas::vector<complex<Energy2> > M(5);
    	Complex         d34 = eps1[i1].dot(eps4[i2]);
    	complex<Energy> d42 = eps4[i2].dot(p2);
   	// M0 in paper
   	M(0) = qw.dimensionedWave().slash(eps1[i1])
   	  .slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps4[i2]);
    	// M4 in paper
    	M(2) =  current.dot(eps4[i2])*d31;
   	// M5 in paper
   	M(3) = -current.dot(eps1[i1])*d42;
   	// M1 in paper (missing factor)
   	M(1) =  current.dot(p4);
   	// M6 in paper
   	M(4) = M(1)*d31*d42/GeV2;
   	// M1 final factor
   	M(1) *= d34;
   	// coefficient of different contributions
   	boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3);
   	if(iq==0) {
   	  axpy_prod_local(bornQQGGweights,M,Cborn);
   	  axpy_prod_local(EWQQGGweights  ,M,CEW  );
   	}
   	else {
   	  axpy_prod_local(bornRRGGweights,M,Cborn);
   	  axpy_prod_local(EWRRGGweights  ,M,CEW  );
   	}
   	unsigned int ioff = (Cborn.size()==6 && abs(qbar->id())%2!=0) ? 3 : 0;
   	for(unsigned int ix=0;ix<3;++ix) {
   	  for(unsigned int iy=0;iy<3;++iy) {
   	    bornME(ix,iy) += Cborn(ix+ioff)*conj(Cborn(iy+ioff));
   	    EWME  (ix,iy) += CEW  (ix+ioff)*conj(CEW  (iy+ioff));
  	  }
  	}
       }
     }
   }
   double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
   double EW   = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
   return EW/born;
 }
 
 double ElectroWeakReweighter::reweightqqbarqqbarS() const {
   // momenta and invariants
   Lorentz5Momentum p1    = subProcess()->incoming().first ->momentum();
   tcPDPtr          q1    = subProcess()->incoming().first ->dataPtr();
   Lorentz5Momentum p2    = subProcess()->incoming().second->momentum();
   tcPDPtr          q1bar = subProcess()->incoming().second->dataPtr();
   if(q1->id()<0) {
     swap(p1,p2    );
     swap(q1 ,q1bar);
   }
   Lorentz5Momentum p3    = subProcess()->outgoing()[0]->momentum();
   tcPDPtr          q2bar = subProcess()->outgoing()[0]->dataPtr();
   Lorentz5Momentum p4    = subProcess()->outgoing()[1]->momentum();
   tcPDPtr          q2    = subProcess()->outgoing()[1]->dataPtr();
   if(q2bar->id()>0) {
     swap(p3,p4    );
     swap(q2 ,q2bar);
   }
   Energy2 s = (p1+p2).m2();
   Energy2 t = (p1-p4).m2();
   Energy2 u = (p1-p3).m2();
   // boost to partonci rest frame
   Lorentz5Momentum psum=p1+p2;
   LorentzRotation boost(-psum.boostVector());
   p1 *= boost;
   p2 *= boost;
   p3 *= boost;
   p4 *= boost;
   p3.setMass(ZERO);
   p3.rescaleRho();
   p4.setMass(ZERO);
   p4.rescaleRho();
   // LO and EW corrected matrix element coefficients
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     bornLLLLWeights,bornLLRRWeights,bornRRLLWeights,bornRRRRWeights,
     EWLLLLWeights,EWLLRRWeights,EWRRLLWeights,EWRRRRWeights;
   bool ident = q1->id()==q2->id();
   // LL -> LL
   if((q1->id()<=4&& q2->id()<=4)|| (q1->id()==5 && q2->id()==5)) {
     if(!ident) {
       bornLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,true ,0);
       EWLLLLWeights   = evaluateRunning(EWProcess::QQQQ,s,t,u,false,0);
     }
     else {
       bornLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,true ,0);
       EWLLLLWeights   = evaluateRunning(EWProcess::QQQQiden,s,t,u,false,0);
     }
   }
   else if(q1->id()==5 || q2->id()>=5) {
     bornLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,true ,0);
     EWLLLLWeights   = evaluateRunning(EWProcess::QtQtQQ,s,t,u,false,0);
   }
   else
     assert(false);
   // RR -> LL
   if(q1->id()%2==0) {
     if(q2->id()<5) {
       bornRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,0);
       EWRRLLWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,0);
     }
     else {
       bornRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,0);
       EWRRLLWeights   = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,0);
     }
   }
   else {
     if(q2->id()<5) {
       bornRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,0);
       EWRRLLWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,0);
     }
     else {
       bornRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,0);
       EWRRLLWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,0);
     }
   }
   // LL -> RR
   if(q1->id()<=4) {
     if(q2->id()%2!=0) {
       bornLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,0);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,0);
     }
     else if (q2->id()==6) {
       bornLLRRWeights = evaluateRunning(EWProcess::QQtRtR,s,t,u,true ,0);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQtRtR,s,t,u,false,0);
     }
     else {
       bornLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,0);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,0);
     }
   }
   else {
     if(q2->id()%2!=0) {
       bornLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,0);
       EWLLRRWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,0);
     }
     else {
       bornLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,0);
       EWLLRRWeights   = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,0);
     }
   }
   // RR -> RR
   if(q1->id()%2==0) {
     if(q2->id()==6) {
       bornRRRRWeights = evaluateRunning(EWProcess::tRtRUU,s,t,u,true ,0);
       EWRRRRWeights   = evaluateRunning(EWProcess::tRtRUU,s,t,u,false,0);
     }
     else if(q2->id()%2==0) {
       if(ident) {
 	bornRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,true ,0);
 	EWRRRRWeights   = evaluateRunning(EWProcess::UUUUiden,s,t,u,false,0);
       }
       else {
 	bornRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,true ,0);
 	EWRRRRWeights   = evaluateRunning(EWProcess::UUUU,s,t,u,false,0);
       }
     }
     else {
       bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,0);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUDD,s,t,u,false,0);
     }
   }
   else {
     if(q2->id()==6) {
       bornRRRRWeights = evaluateRunning(EWProcess::tRtRDD,s,t,u,true ,0);
       EWRRRRWeights   = evaluateRunning(EWProcess::tRtRDD,s,t,u,false,0);
     }
     else if(q2->id()%2==0) {
       bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,0);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUDD,s,t,u,false,0);
     }
     else {
       if(ident) {
 	bornRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,true ,0);
 	EWRRRRWeights   = evaluateRunning(EWProcess::DDDDiden,s,t,u,false,0);
       }
       else {
 	bornRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,true ,0);
 	EWRRRRWeights   = evaluateRunning(EWProcess::DDDD,s,t,u,false,0);
       }
     }
   }
   // extra terms for identical particles
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     borntChannelWeights,EWtChannelWeights;
   if(ident) {
     if(q1->id()%2==0) {
       borntChannelWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,1);
       EWtChannelWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,1);
     }
     else if(q1->id()==5) {
       borntChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,1);
       EWtChannelWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,1);
     }
     else {
       borntChannelWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,1);
       EWtChannelWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,1);
     }
   }
   SpinorWaveFunction       q1w(p1,q1   ,incoming);
   SpinorBarWaveFunction q1barw(p2,q1bar,incoming);
   SpinorWaveFunction    q2barw(p3,q2bar,outgoing);
   SpinorBarWaveFunction    q2w(p4,q2   ,outgoing);
   boost::numeric::ublas::matrix<Complex>
     bornME = boost::numeric::ublas::zero_matrix<Complex>(2,2),
     EWME   = boost::numeric::ublas::zero_matrix<Complex>(2,2);
   for(unsigned int iq1=0;iq1<2;++iq1) {
     if(iq1==0) {
       q1w.reset   (0);
       q1barw.reset(1);
     }
     else {
       q1w.reset   (1);
       q1barw.reset(0);
     }
     LorentzVector<complex<Energy> > current1 =
       q1w.dimensionedWave().vectorCurrent(q1barw.dimensionedWave());
     for(unsigned int iq2=0;iq2<2;++iq2) {
       if(iq2==0) {
 	q2w.reset   (0);
 	q2barw.reset(1);
       }
       else {
 	q2w.reset   (1);
 	q2barw.reset(0);
       }
       LorentzVector<complex<Energy> > current2 =
 	q2barw.dimensionedWave().vectorCurrent(q2w.dimensionedWave());
       complex<Energy2> amp = current1.dot(current2);
       vector<Complex> Cborn(2),CEW(2);
       // amplitudes
       if(iq1==0) {
 	// LL
 	if(iq2==0) {
 	  unsigned int ioff;
 	  if(q1->id()%2==0) {
 	    ioff = q2->id()%2==0 ? 0 : 2;
 	  }
 	  else {
 	    ioff = q2->id()%2==0 ? 1 : 3;
 	  }
 	  for(unsigned int ix=0;ix<2;++ix) {
 	    Cborn[ix] = amp*bornLLLLWeights(6*ix+ioff,0);
 	    CEW  [ix] = amp*  EWLLLLWeights(6*ix+ioff,0);
 	  }
 	}
 	// LR
 	else {
 	  unsigned int ioff =  q1->id()%2==0 ? 0 : 1;
 	  for(unsigned int ix=0;ix<2;++ix) {
 	    Cborn[ix] = amp*bornLLRRWeights(2*ix+ioff,0);
 	    CEW  [ix] = amp*  EWLLRRWeights(2*ix+ioff,0);
 	  }
 	}
       }
       else {
 	if(iq2==0) {
 	  unsigned int ioff=q2->id()%2==0 ? 0 : 1;
 	  for(unsigned int ix=0;ix<2;++ix) {
 	    Cborn[ix] = amp*bornRRLLWeights(2*ix+ioff,0);
 	    CEW  [ix] = amp*  EWRRLLWeights(2*ix+ioff,0);
 	  }
 	}
 	else {
 	  for(unsigned int ix=0;ix<2;++ix) {
 	    Cborn[ix] = amp*bornRRRRWeights(ix,0);
 	    CEW  [ix] = amp*  EWRRRRWeights(ix,0);
 	  }
 	}
       }
       // square
       for(unsigned int ix=0;ix<2;++ix) {
 	for(unsigned int iy=0;iy<2;++iy) {
 	  bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]);
 	  EWME  (ix,iy) += CEW  [ix]*conj(CEW  [iy]);
 	}
       }
     }
   }
   // extra t-channel pieces if identical flavours
   if(ident) {
     for(unsigned int iq1=0;iq1<2;++iq1) {
       q1w.reset(iq1);
       q2w.reset(iq1);
       LorentzVector<complex<Energy> > current1 =
 	q1w.dimensionedWave().vectorCurrent(q2w.dimensionedWave());
       q1barw.reset(iq1);
       q2barw.reset(iq1);
       LorentzVector<complex<Energy> > current2 =
 	q2barw.dimensionedWave().vectorCurrent(q1barw.dimensionedWave());
       complex<Energy2> amp = current1.dot(current2);
       vector<Complex> Cborn(2),CEW(2);
       unsigned int ioff =  q1->id()%2==0 ? 0 : 1;
       for(unsigned int ix=0;ix<2;++ix) {
 	Cborn[ix] = amp*borntChannelWeights(2*ix+ioff,0);
 	CEW  [ix] = amp*  EWtChannelWeights(2*ix+ioff,0);
       }
       // square
       for(unsigned int ix=0;ix<2;++ix) {
 	for(unsigned int iy=0;iy<2;++iy) {
 	  bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]);
 	  EWME  (ix,iy) += CEW  [ix]*conj(CEW  [iy]);
 	}
       }
     }
   }
   // colour factors
   double born = 2.*real(bornME(0,0))+9.*real(bornME(1,1));
   double EW   = 2.*real(  EWME(0,0))+9.*real(  EWME(1,1));
   return EW/born;
 }
 
 double ElectroWeakReweighter::reweightqqbarqqbarT() const {
   // momenta and invariants
   Lorentz5Momentum p1    = subProcess()->incoming().first ->momentum();
   tcPDPtr          q1    = subProcess()->incoming().first ->dataPtr();
   Lorentz5Momentum p2    = subProcess()->incoming().second->momentum();
   tcPDPtr          q1bar = subProcess()->incoming().second->dataPtr();
   if(q1->id()<0) {
     swap(p1,p2    );
     swap(q1 ,q1bar);
   }
   Lorentz5Momentum p3    = subProcess()->outgoing()[0]->momentum();
   tcPDPtr          q2bar = subProcess()->outgoing()[0]->dataPtr();
   Lorentz5Momentum p4    = subProcess()->outgoing()[1]->momentum();
   tcPDPtr          q2    = subProcess()->outgoing()[1]->dataPtr();
   if(q2bar->id()>0) {
     swap(p3,p4    );
     swap(q2 ,q2bar);
   }
   Energy2 s = (p1+p2).m2();
   Energy2 t = (p1-p4).m2();
   Energy2 u = (p1-p3).m2();
   // boost to partonci rest frame
   Lorentz5Momentum psum=p1+p2;
   LorentzRotation boost(-psum.boostVector());
   p1 *= boost;
   p2 *= boost;
   p3 *= boost;
   p4 *= boost;
   p3.setMass(ZERO);
   p3.rescaleRho();
   p4.setMass(ZERO);
   p4.rescaleRho();
   assert(q1==q2 && q1bar==q2bar);
   assert( q1->id() != -q1bar->id() && q2->id() != -q2bar->id() );
   // LO and EW corrected matrix element coefficients
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     bornLLLLWeights,bornLLRRWeights,bornRRLLWeights,bornRRRRWeights,
     EWLLLLWeights,EWLLRRWeights,EWRRLLWeights,EWRRRRWeights;
   // LL
   if( q1->id()    == ParticleID::b || 
       q1bar->id() == ParticleID::bbar ) {
     bornLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,true ,1);
     EWLLLLWeights   = evaluateRunning(EWProcess::QtQtQQ,s,t,u,false,1);
   }
   else {
     bornLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,true ,1);
     EWLLLLWeights   = evaluateRunning(EWProcess::QQQQ,s,t,u,false,1);
   }
   // RR -> LL
   if(q1->id()%2==0) {
     if(q1bar->id()==ParticleID::bbar) {
       bornRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,1);
       EWRRLLWeights   = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,1);
     }
     else {
       bornRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,1);
       EWRRLLWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,1);
     }
   }
   else {
     if(q1bar->id()==ParticleID::bbar) {
       bornRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,1);
       EWRRLLWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,1);
     }
     else {
       bornRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,1);
       EWRRLLWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,1);
     }
   }
   // LL -> RR
   if(abs(q1bar->id())%2==0) {
     if(q1->id()==ParticleID::b) {
       bornLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,1);
       EWLLRRWeights   = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,1);
     }
     else {
       bornLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,1);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,1);
     }
   }
   else {
     if(q1->id()==ParticleID::b) {
       bornLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,1);
       EWLLRRWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,1);
     }
     else {
       bornLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,1);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,1);
     }
   }
   // RR -> RR
   if(q1->id()%2==0) {
     if(abs(q1bar->id())%2==0) {
       bornRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,true ,1);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUUU,s,t,u,false,1);
     }
     else {
       bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,1);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUDD,s,t,u,false,1);
     }
   }
   else {
     if(abs(q1bar->id())%2==0) {
       bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,1);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUDD,s,t,u,false,1);
     }
     else {
       bornRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,true ,1);
       EWRRRRWeights   = evaluateRunning(EWProcess::DDDD,s,t,u,false,1);
     }
   }
   // calculate the spinors
   SpinorWaveFunction       q1w(p1,q1   ,incoming);
   SpinorBarWaveFunction q1barw(p2,q1bar,incoming);
   SpinorWaveFunction    q2barw(p3,q2bar,outgoing);
   SpinorBarWaveFunction    q2w(p4,q2   ,outgoing);
   boost::numeric::ublas::matrix<Complex>
     bornME = boost::numeric::ublas::zero_matrix<Complex>(2,2),
     EWME   = boost::numeric::ublas::zero_matrix<Complex>(2,2);
   for(unsigned int iq1=0;iq1<2;++iq1) {
     q1w.reset(iq1);
     q2w.reset(iq1);
     LorentzVector<complex<Energy> > current1 =
       q1w.dimensionedWave().vectorCurrent(q2w.dimensionedWave());
     for(unsigned int iq2=0;iq2<2;++iq2) {
       q1barw.reset(iq2);
       q2barw.reset(iq2);
       LorentzVector<complex<Energy> > current2 =
 	q2barw.dimensionedWave().vectorCurrent(q1barw.dimensionedWave());
       // calculate the amplitude
       complex<Energy2> amp = current1.dot(current2);
       vector<Complex> Cborn(2),CEW(2);
       if(iq1==0) {
 	// LL RR
 	if(iq2==0) {
 	  unsigned int ioff =  q1->id()%2==0 ? 0 : 1;
 	  for(unsigned int ix=0;ix<2;++ix) {
 	    Cborn[ix] = amp*bornLLRRWeights(2*ix+ioff,0);
 	    CEW  [ix] = amp*  EWLLRRWeights(2*ix+ioff,0);
 	  }
 	}
 	// LL LL
 	else {
 	  unsigned int ioff;
 	  if(q1->id()%2==0) {
 	    ioff = abs(q1bar->id())%2==0 ? 0 : 2;
 	  }
 	  else {
 	    ioff = abs(q1bar->id())%2==0 ? 1 : 3;
 	  }
 	  for(unsigned int ix=0;ix<2;++ix) {
 	    Cborn[ix] = amp*bornLLLLWeights(6*ix+ioff,0);
 	    CEW  [ix] = amp*  EWLLLLWeights(6*ix+ioff,0);
  	  }
 	}
       }
       else {
 	// RR RR
 	if(iq2==0) {
 	  for(unsigned int ix=0;ix<2;++ix) {
 	    Cborn[ix] = amp*bornRRRRWeights(ix,0);
 	    CEW  [ix] = amp*  EWRRRRWeights(ix,0);
 	  }
 	}
 	// RR LL
 	else {
 	  unsigned int ioff=abs(q1bar->id())%2==0 ? 0 : 1;
 	  for(unsigned int ix=0;ix<2;++ix) {
 	    Cborn[ix] = amp*bornRRLLWeights(2*ix+ioff,0);
 	    CEW  [ix] = amp*  EWRRLLWeights(2*ix+ioff,0);
 	  }
 	}
       }
       // square
       for(unsigned int ix=0;ix<2;++ix) {
 	for(unsigned int iy=0;iy<2;++iy) {
 	  bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]);
 	  EWME  (ix,iy) += CEW  [ix]*conj(CEW  [iy]);
 	}
       }
     }
   }
   // colour factors
   double born = 2.*real(bornME(0,0))+9.*real(bornME(1,1));
   double EW   = 2.*real(  EWME(0,0))+9.*real(  EWME(1,1));
   return EW/born;
 }
 
 double ElectroWeakReweighter::reweightqqqq() const {
   // momenta and invariants
   Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum();
   tcPDPtr          q1 = subProcess()->incoming().first ->dataPtr();
   Lorentz5Momentum p2 = subProcess()->incoming().second->momentum();
   tcPDPtr          q2 = subProcess()->incoming().second->dataPtr();
   Lorentz5Momentum p3 = subProcess()->outgoing()[0]    ->momentum();
   tcPDPtr          q3 = subProcess()->outgoing()[0]    ->dataPtr();
   Lorentz5Momentum p4 = subProcess()->outgoing()[1]    ->momentum();
   tcPDPtr          q4 = subProcess()->outgoing()[1]    ->dataPtr();
   if(q1->id()!=q3->id()) {
     swap(q3,q4);
     swap(p3,p4);
   }
   assert(q1->id()==q3->id());
   assert(q2->id()==q4->id());
   Energy2 s = (p1+p2).m2();
   Energy2 t = (p1-p4).m2();
   Energy2 u = (p1-p3).m2();
   // boost to partonci rest frame
   Lorentz5Momentum psum=p1+p2;
   LorentzRotation boost(-psum.boostVector());
   p1 *= boost;
   p2 *= boost;
   p3 *= boost;
   p4 *= boost;
   p3.setMass(ZERO);
   p3.rescaleRho();
   p4.setMass(ZERO);
   p4.rescaleRho();
   // LO and EW corrected matrix element coefficients
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     bornLLLLWeights,bornLLRRWeights,bornRRLLWeights,bornRRRRWeights,
     EWLLLLWeights,EWLLRRWeights,EWRRLLWeights,EWRRRRWeights;
   bool ident = q1->id()==q2->id();
   // LL -> LL
   if((q1->id()<=4&& q2->id()<=4)|| (q1->id()==5 && q2->id()==5)) {
     if(!ident) {
       bornLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,true ,2);
       EWLLLLWeights   = evaluateRunning(EWProcess::QQQQ,s,t,u,false,2);
     }
     else {
       bornLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,true ,2);
       EWLLLLWeights   = evaluateRunning(EWProcess::QQQQiden,s,t,u,false,2);
     }
   }
   else if(q1->id()==5 || q2->id()==5) {
     bornLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,true ,2);
     EWLLLLWeights   = evaluateRunning(EWProcess::QtQtQQ,s,t,u,false,2);
   }
   else
     assert(false);
   // RR -> LL
   if(q1->id()%2==0) {
     if(q2->id()<5) {
       bornRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,2);
       EWRRLLWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,2);
     }
     else {
       bornRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,2);
       EWRRLLWeights   = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,2);
     }
   }
   else {
     if(q2->id()<5) {
       bornRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,2);
       EWRRLLWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,2);
     }
     else {
       bornRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,2);
       EWRRLLWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,2);
     }
   }
   // LL -> RR
   if(q1->id()<=4) {
     if(q2->id()%2!=0) {
       bornLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,2);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,2);
     }
     else {
       bornLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,2);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,2);
     }
   }
   else {
     if(q2->id()%2!=0) {
       bornLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,2);
       EWLLRRWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,2);
     }
     else {
       bornLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,2);
       EWLLRRWeights   = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,2);
     }
   }
   // RR -> RR
   if(q1->id()%2==0) {
     if(q2->id()%2==0) {
       if(ident) {
   	bornRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,true ,2);
   	EWRRRRWeights   = evaluateRunning(EWProcess::UUUUiden,s,t,u,false,2);
       }
       else {
   	bornRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,true ,2);
   	EWRRRRWeights   = evaluateRunning(EWProcess::UUUU,s,t,u,false,2);
       }
     }
     else {
       bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,2);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUDD,s,t,u,false,2);
     }
   }
   else {
     if(q2->id()%2==0) {
       bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,2);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUDD,s,t,u,false,2);
     }
     else {
       if(ident) {
   	bornRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,true ,2);
   	EWRRRRWeights   = evaluateRunning(EWProcess::DDDDiden,s,t,u,false,2);
       }
       else {
   	bornRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,true ,2);
   	EWRRRRWeights   = evaluateRunning(EWProcess::DDDD,s,t,u,false,2);
       }
     }
   }
   // extra terms for identical particles
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     borntChannelWeights,EWtChannelWeights;
   if(ident) {
     if(q1->id()%2==0) {
       borntChannelWeights = evaluateRunning(EWProcess::QQUU,s,u,t,true ,2);
       EWtChannelWeights   = evaluateRunning(EWProcess::QQUU,s,u,t,false,2);
     }
     else if(q1->id()==5) {
       borntChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,u,t,true ,2);
       EWtChannelWeights   = evaluateRunning(EWProcess::QtQtDD,s,u,t,false,2);
     }
     else {
       borntChannelWeights = evaluateRunning(EWProcess::QQDD,s,u,t,true ,2);
       EWtChannelWeights   = evaluateRunning(EWProcess::QQDD,s,u,t,false,2);
     }
   }
   SpinorWaveFunction    q1w(p1,q1,incoming);
   SpinorWaveFunction    q2w(p2,q2,incoming);
   SpinorBarWaveFunction q3w(p3,q3,outgoing);
   SpinorBarWaveFunction q4w(p4,q4,outgoing);
   boost::numeric::ublas::matrix<Complex>
     bornME = boost::numeric::ublas::zero_matrix<Complex>(2,2),
     EWME   = boost::numeric::ublas::zero_matrix<Complex>(2,2);
   for(unsigned int iq1=0;iq1<2;++iq1) {
     q1w.reset(iq1);
     q3w.reset(iq1);
     LorentzVector<complex<Energy> > current1 =
       q1w.dimensionedWave().vectorCurrent(q3w.dimensionedWave());
     for(unsigned int iq2=0;iq2<2;++iq2) {
       q2w.reset(iq2);
       q4w.reset(iq2);
       LorentzVector<complex<Energy> > current2 =
   	q2w.dimensionedWave().vectorCurrent(q4w.dimensionedWave());
       complex<Energy2> amp = current1.dot(current2);
       vector<Complex> Cborn(2),CEW(2);
       // amplitudes
       if(iq1==0) {
  	// LL
  	if(iq2==0) {
  	  unsigned int ioff;
  	  if(q1->id()%2==0) {
  	    ioff = q2->id()%2==0 ? 0 : 2;
  	  }
  	  else {
  	    ioff = q2->id()%2==0 ? 1 : 3;
  	  }
  	  for(unsigned int ix=0;ix<2;++ix) {
  	    Cborn[ix] = amp*bornLLLLWeights(6*ix+ioff,0);
  	    CEW  [ix] = amp*  EWLLLLWeights(6*ix+ioff,0);
  	  }
  	}
  	// LR
  	else {
  	  unsigned int ioff =  q1->id()%2==0 ? 0 : 1;
  	  for(unsigned int ix=0;ix<2;++ix) {
  	    Cborn[ix] = amp*bornLLRRWeights(2*ix+ioff,0);
  	    CEW  [ix] = amp*  EWLLRRWeights(2*ix+ioff,0);
  	  }
  	}
       }
       else {
   	if(iq2==0) {
   	  unsigned int ioff=q2->id()%2==0 ? 0 : 1;
   	  for(unsigned int ix=0;ix<2;++ix) {
   	    Cborn[ix] = amp*bornRRLLWeights(2*ix+ioff,0);
   	    CEW  [ix] = amp*  EWRRLLWeights(2*ix+ioff,0);
   	  }
   	}
   	else {
   	  for(unsigned int ix=0;ix<2;++ix) {
   	    Cborn[ix] = amp*bornRRRRWeights(ix,0);
   	    CEW  [ix] = amp*  EWRRRRWeights(ix,0);
   	  }
   	}
       }
       // square
       for(unsigned int ix=0;ix<2;++ix) {
   	for(unsigned int iy=0;iy<2;++iy) {
   	  bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]);
   	  EWME  (ix,iy) += CEW  [ix]*conj(CEW  [iy]);
   	}
       }
     }
   }
   // extra u-channel pieces if identical flavours
   if(ident) {
     for(unsigned int iq1=0;iq1<2;++iq1) {
       q1w.reset(iq1);
       q4w.reset(iq1);
       LorentzVector<complex<Energy> > current1 =
 	q1w.dimensionedWave().vectorCurrent(q4w.dimensionedWave());
       if(iq1==0) {
 	q2w.reset(1);
 	q3w.reset(1);
       }
       else {
 	q2w.reset(0);
 	q3w.reset(0);
       }
       LorentzVector<complex<Energy> > current2 =
 	q2w.dimensionedWave().vectorCurrent(q3w.dimensionedWave());
       complex<Energy2> amp = current1.dot(current2);
       vector<Complex> Cborn(2),CEW(2);
       unsigned int ioff =  q1->id()%2==0 ? 0 : 1;
       for(unsigned int ix=0;ix<2;++ix) {
   	Cborn[ix] = amp*borntChannelWeights(2*ix+ioff,0);
    	CEW  [ix] = amp*  EWtChannelWeights(2*ix+ioff,0);
       }
       // square
       for(unsigned int ix=0;ix<2;++ix) {
   	for(unsigned int iy=0;iy<2;++iy) {
   	  bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]);
   	  EWME  (ix,iy) += CEW  [ix]*conj(CEW  [iy]);
  	}
       }
     }
   }
   // colour factors
   double born = 2.*real(bornME(0,0))+9.*real(bornME(1,1));
   double EW   = 2.*real(  EWME(0,0))+9.*real(  EWME(1,1));
   return EW/born;
 }
 
 double ElectroWeakReweighter::reweightqbarqbarqbarqbar() const {
   // momenta and invariants
   Lorentz5Momentum p1    = subProcess()->incoming().first ->momentum();
   tcPDPtr          qbar1 = subProcess()->incoming().first ->dataPtr();
   Lorentz5Momentum p2    = subProcess()->incoming().second->momentum();
   tcPDPtr          qbar2 = subProcess()->incoming().second->dataPtr();
   Lorentz5Momentum p3    = subProcess()->outgoing()[0]    ->momentum();
   tcPDPtr          qbar3 = subProcess()->outgoing()[0]    ->dataPtr();
   Lorentz5Momentum p4    = subProcess()->outgoing()[1]    ->momentum();
   tcPDPtr          qbar4 = subProcess()->outgoing()[1]    ->dataPtr();
   if(qbar1->id()!=qbar3->id()) {
     swap(qbar3,qbar4);
     swap(p3,p4);
   }
   assert(qbar1->id()==qbar3->id());
   assert(qbar2->id()==qbar4->id());
   Energy2 s = (p1+p2).m2();
   Energy2 t = (p1-p4).m2();
   Energy2 u = (p1-p3).m2();
   // boost to partonic rest frame
   Lorentz5Momentum psum=p1+p2;
   LorentzRotation boost(-psum.boostVector());
   p1 *= boost;
   p2 *= boost;
   p3 *= boost;
   p4 *= boost;
   p3.setMass(ZERO);
   p3.rescaleRho();
   p4.setMass(ZERO);
   p4.rescaleRho();
   // LO and EW corrected matrix element coefficients
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     bornLLLLWeights,bornLLRRWeights,bornRRLLWeights,bornRRRRWeights,
     EWLLLLWeights,EWLLRRWeights,EWRRLLWeights,EWRRRRWeights;
   bool ident = qbar1->id()==qbar2->id();
   // LL -> LL
   if((abs(qbar1->id())<=4 && abs(qbar2->id())<=4) ||
      (abs(qbar1->id())==5 && abs(qbar2->id())==5)) {
     if(!ident) {
       bornLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,true ,2);
       EWLLLLWeights   = evaluateRunning(EWProcess::QQQQ,s,t,u,false,2);
     }
     else {
       bornLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,true ,2);
       EWLLLLWeights   = evaluateRunning(EWProcess::QQQQiden,s,t,u,false,2);
     }
   }
   else if(abs(qbar1->id())==5 || abs(qbar2->id())==5) {
     bornLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,true ,2);
     EWLLLLWeights   = evaluateRunning(EWProcess::QtQtQQ,s,t,u,false,2);
   }
   else
     assert(false);
   // RR -> LL
   if(abs(qbar1->id())%2==0) {
     if(abs(qbar2->id())<5) {
       bornRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,2);
       EWRRLLWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,2);
     }
     else {
       bornRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,2);
       EWRRLLWeights   = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,2);
     }
   }
   else {
     if(abs(qbar2->id())<5) {
       bornRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,2);
       EWRRLLWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,2);
     }
     else {
       bornRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,2);
       EWRRLLWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,2);
     }
   }
   // LL -> RR
   if(abs(qbar1->id())<=4) {
     if(abs(qbar2->id())%2!=0) {
       bornLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,2);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQDD,s,t,u,false,2);
     }
     else {
       bornLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,2);
       EWLLRRWeights   = evaluateRunning(EWProcess::QQUU,s,t,u,false,2);
     }
   }
   else {
     if(abs(qbar2->id())%2!=0) {
       bornLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,2);
       EWLLRRWeights   = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,2);
     }
     else {
       bornLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,2);
       EWLLRRWeights   = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,2);
     }
   }
   // RR -> RR
   if(abs(qbar1->id())%2==0) {
     if(abs(qbar2->id())%2==0) {
       if(ident) {
   	bornRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,true ,2);
   	EWRRRRWeights   = evaluateRunning(EWProcess::UUUUiden,s,t,u,false,2);
       }
       else {
   	bornRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,true ,2);
   	EWRRRRWeights   = evaluateRunning(EWProcess::UUUU,s,t,u,false,2);
       }
     }
     else {
       bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,2);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUDD,s,t,u,false,2);
     }
   }
   else {
     if(abs(qbar2->id())%2==0) {
       bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,2);
       EWRRRRWeights   = evaluateRunning(EWProcess::UUDD,s,t,u,false,2);
     }
     else {
       if(ident) {
   	bornRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,true ,2);
   	EWRRRRWeights   = evaluateRunning(EWProcess::DDDDiden,s,t,u,false,2);
       }
       else {
   	bornRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,true ,2);
   	EWRRRRWeights   = evaluateRunning(EWProcess::DDDD,s,t,u,false,2);
       }
     }
   }
   // extra terms for identical particles
   boost::numeric::ublas::matrix<complex<InvEnergy2> >
     borntChannelWeights,EWtChannelWeights;
   if(ident) {
     if(abs(qbar1->id())%2==0) {
       borntChannelWeights = evaluateRunning(EWProcess::QQUU,s,u,t,true ,2);
       EWtChannelWeights   = evaluateRunning(EWProcess::QQUU,s,u,t,false,2);
     }
     else if(abs(qbar1->id())==5) {
       borntChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,u,t,true ,2);
       EWtChannelWeights   = evaluateRunning(EWProcess::QtQtDD,s,u,t,false,2);
     }
     else {
       borntChannelWeights = evaluateRunning(EWProcess::QQDD,s,u,t,true ,2);
       EWtChannelWeights   = evaluateRunning(EWProcess::QQDD,s,u,t,false,2);
     }
   }
   SpinorBarWaveFunction qbar1w(p1,qbar1,incoming);
   SpinorBarWaveFunction qbar2w(p2,qbar2,incoming);
   SpinorWaveFunction    qbar3w(p3,qbar3,outgoing);
   SpinorWaveFunction    qbar4w(p4,qbar4,outgoing);
   boost::numeric::ublas::matrix<Complex>
     bornME = boost::numeric::ublas::zero_matrix<Complex>(2,2),
     EWME   = boost::numeric::ublas::zero_matrix<Complex>(2,2);
   for(unsigned int iq1=0;iq1<2;++iq1) {
     qbar1w.reset(iq1);
     qbar3w.reset(iq1);
     LorentzVector<complex<Energy> > current1 =
       qbar3w.dimensionedWave().vectorCurrent(qbar1w.dimensionedWave());
     for(unsigned int iq2=0;iq2<2;++iq2) {
       qbar2w.reset(iq2);
       qbar4w.reset(iq2);
       LorentzVector<complex<Energy> > current2 =
   	qbar4w.dimensionedWave().vectorCurrent(qbar2w.dimensionedWave());
       complex<Energy2> amp = current1.dot(current2);
       vector<Complex> Cborn(2),CEW(2);
       // amplitudes
       if(iq1==1) {
  	// LL
  	if(iq2==1) {
  	  unsigned int ioff;
  	  if(abs(qbar1->id())%2==0) {
  	    ioff = abs(qbar2->id())%2==0 ? 0 : 2;
  	  }
  	  else {
  	    ioff = abs(qbar2->id())%2==0 ? 1 : 3;
  	  }
  	  for(unsigned int ix=0;ix<2;++ix) {
  	    Cborn[ix] = amp*bornLLLLWeights(6*ix+ioff,0);
  	    CEW  [ix] = amp*  EWLLLLWeights(6*ix+ioff,0);
  	  }
  	}
  	// LR
  	else {
  	  unsigned int ioff =  abs(qbar1->id())%2==0 ? 0 : 1;
  	  for(unsigned int ix=0;ix<2;++ix) {
  	    Cborn[ix] = amp*bornLLRRWeights(2*ix+ioff,0);
  	    CEW  [ix] = amp*  EWLLRRWeights(2*ix+ioff,0);
  	  }
  	}
       }
       else {
   	if(iq2==1) {
   	  unsigned int ioff=abs(qbar2->id())%2==0 ? 0 : 1;
   	  for(unsigned int ix=0;ix<2;++ix) {
   	    Cborn[ix] = amp*bornRRLLWeights(2*ix+ioff,0);
   	    CEW  [ix] = amp*  EWRRLLWeights(2*ix+ioff,0);
   	  }
   	}
   	else {
   	  for(unsigned int ix=0;ix<2;++ix) {
   	    Cborn[ix] = amp*bornRRRRWeights(ix,0);
   	    CEW  [ix] = amp*  EWRRRRWeights(ix,0);
   	  }
   	}
       }
       // square
       for(unsigned int ix=0;ix<2;++ix) {
  	for(unsigned int iy=0;iy<2;++iy) {
  	  bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]);
  	  EWME  (ix,iy) += CEW  [ix]*conj(CEW  [iy]);
  	}
       }
     }
   }
   // extra u-channel pieces if identical flavours
   if(ident) {
     for(unsigned int iq1=0;iq1<2;++iq1) {
       qbar1w.reset(iq1);
       qbar4w.reset(iq1);
       LorentzVector<complex<Energy> > current1 =
 	qbar4w.dimensionedWave().vectorCurrent(qbar1w.dimensionedWave());
       if(iq1==0) {
 	qbar2w.reset(1);
 	qbar3w.reset(1);
       }
       else {
 	qbar2w.reset(0);
 	qbar3w.reset(0);
       }
       LorentzVector<complex<Energy> > current2 =
 	qbar3w.dimensionedWave().vectorCurrent(qbar2w.dimensionedWave());
       complex<Energy2> amp = current1.dot(current2);
       vector<Complex> Cborn(2),CEW(2);
       unsigned int ioff =  abs(qbar1->id())%2==0 ? 0 : 1;
       for(unsigned int ix=0;ix<2;++ix) {
   	Cborn[ix] = amp*borntChannelWeights(2*ix+ioff,0);
    	CEW  [ix] = amp*  EWtChannelWeights(2*ix+ioff,0);
       }
       // square
       for(unsigned int ix=0;ix<2;++ix) {
   	for(unsigned int iy=0;iy<2;++iy) {
   	  bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]);
   	  EWME  (ix,iy) += CEW  [ix]*conj(CEW  [iy]);
  	}
       }
     }
   }
   // colour factors
   double born = 2.*real(bornME(0,0))+9.*real(bornME(1,1));
   double EW   = 2.*real(  EWME(0,0))+9.*real(  EWME(1,1));
   return EW/born;
 }
diff --git a/MatrixElement/EW/HighEnergyMatching.cc b/MatrixElement/EW/HighEnergyMatching.cc
--- a/MatrixElement/EW/HighEnergyMatching.cc
+++ b/MatrixElement/EW/HighEnergyMatching.cc
@@ -1,1179 +1,1181 @@
 // -*- C++ -*-
 //
 // HighEnergyMatching.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 //
 // Herwig is licenced under version 2 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 //
 //
 #include "HighEnergyMatching.h"
 #include "ElectroWeakReweighter.h"
 #include "GroupInvariants.h"
 #include <boost/numeric/ublas/operation.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 
 using namespace Herwig;
 using namespace HighEnergyMatching;
 using namespace GroupInvariants;
 using namespace EWProcess;
 
 namespace {
 
 /**
  *   \f$M_N\f$, this matrix is used for identical particle scattering
  */
 boost::numeric::ublas::matrix<Complex> M_N(unsigned int suN) {
   double N(suN);
   boost::numeric::ublas::matrix<Complex> M(2,2);
   M(0,0) = -1.0/N;
   M(0,1) = 2.0;
   M(1,0) = 0.5 - 1.0/(2.0*N*N);
   M(1,1) = 1.0/N;
   return M;
 }
 
 #ifdef ThePEG_HAS_UNITS_CHECKING
 void axpy_prod_local(const boost::numeric::ublas::matrix<Complex>  & A,
 		     const boost::numeric::ublas::matrix<complex<InvEnergy2> > & B,
 		     boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
   assert(A.size2()==B.size1());
   C.resize(A.size1(),B.size2());
   for(unsigned int ix=0;ix<A.size1();++ix) {
     for(unsigned int iy=0;iy<B.size2();++iy) {
       C(ix,iy) = ZERO;
       for(unsigned int iz=0;iz<A.size2();++iz) {
 	C(ix,iy) += A(ix,iz)*B(iz,iy);
       }
     }
   }
 }
 
 void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> > & A,
 		       const boost::numeric::ublas::matrix<Complex> & B,
 		       boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
   assert(A.size2()==B.size1());
   C.resize(A.size1(),B.size2());
   for(unsigned int ix=0;ix<A.size1();++ix) {
     for(unsigned int iy=0;iy<B.size2();++iy) {
       C(ix,iy) = ZERO;
       for(unsigned int iz=0;iz<A.size2();++iz) {
 	C(ix,iy) += A(ix,iz)*B(iz,iy);
       }
     }
   }
 }
 #else
 void axpy_prod_local(const boost::numeric::ublas::matrix<Complex> & A,
 		     const boost::numeric::ublas::matrix<Complex> & B,
 		     boost::numeric::ublas::matrix<Complex> & C) {
+  assert(A.size2()==B.size1());
+  C.resize(A.size1(),B.size2());
   axpy_prod(A,B,C);
 }
 #endif
 
 }
 
 boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 HighEnergyMatching::highEnergyMatching(Energy highScale, 
 				       Energy2 s, Energy2 t, Energy2 u,
 				       EWProcess::Process process,
 				       bool oneLoop, bool includeAlphaS2) {
   switch (process) {
   case QQQQ:     case QQQQiden:
   case QtQtQQ:   case QQUU:
   case QtQtUU:   case QQtRtR:
   case QQDD:     case QtQtDD:
   case QQLL:     case QQEE:
   case UUUU:     case UUUUiden:
   case tRtRUU:   case UUDD:
   case tRtRDD:   case UULL:
   case UUEE:     case DDDD:
   case DDDDiden: case DDLL:
   case DDEE:     case LLLL:
   case LLLLiden: case LLEE:
   case EEEE:     case EEEEiden:
     return SpinHalfMatching(highScale,s,t,u,process,oneLoop,includeAlphaS2);
     break;         
   case QQWW:
   case QQWG:
   case QQBG:
   case QQGG:
   case QtQtGG:
   case LLWW:
   case UUBB:
   case UUBG:
   case UUGG:
   case tRtRGG:
   case DDBB:
   case DDBG:
   case DDGG:
   case EEBB:
     return Spin1HighMatching(highScale,s,t,u,process,oneLoop,includeAlphaS2);
     break;
   case QQPhiPhi:
   case LLPhiPhi:
   case UUPhiPhi:   
   case DDPhiPhi:
   case EEPhiPhi:
     return Spin0HighMatching(highScale,s,t,u,process,oneLoop,includeAlphaS2);
     break;
   default:
     assert(false);
     break;
   }
 }
 
 boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 HighEnergyMatching::SpinHalfMatching(Energy highScale, 
 				     Energy2 s, Energy2 t, Energy2 u,
 				     EWProcess::Process process, 
 				     bool oneLoop, bool includeAlphaS2) {
   using Constants::pi;
   boost::numeric::ublas::matrix<complex<InvEnergy2> > highC;
   Energy Q = highScale;
   double a1 = ElectroWeakReweighter::coupling()->a1(Q);
   double a2 = ElectroWeakReweighter::coupling()->a2(Q);
   double a3 = ElectroWeakReweighter::coupling()->a3(Q);
   double y_t = ElectroWeakReweighter::coupling()->y_t(Q);
   unsigned int order = !oneLoop ? 0 : 1;
   double Yi(0.),Yf(0.);
    
   Complex Ls = MinusLog(-s/(Q*Q));
    
   double C_A_2 = C_A(2);
   double C_A_3 = C_A(3);
   double C_F_2 = C_F(2);
   double C_F_3 = C_F(3);
   double C_d_2 = C_d(2);
   double C_d_3 = C_d(3);
   double C_1_2 = C_1(2);
   double C_1_3 = C_1(3);
   Complex W = WFunction(Q,s);
   Complex X_2_st    = XNFunction(2,Q,s,t);
   //Complex X_2_su    = XNFunction(2,Q,s,u);
   Complex X_3_st    = XNFunction(3,Q,s,t);
   Complex X_3_su    = XNFunction(3,Q,s,u);
   Complex PI1       = PI1_function(Q,s);
   Complex fTilde_st = fTildeFunction(Q,s,t);
   Complex fTilde_su = fTildeFunction(Q,s,u);
    
   switch (process) {
      
   case QQQQ:
     // NOTE this 4x1 column vector highC is given by (C_11,C_21,C_12,C_22), 
     // where C_12 is the coeff. for the term that transforms under SU(2) but not SU(3)
     Yi = Yf = 1./6.;
     highC.resize(4,1);
     highC(0,0) = ZERO;
     highC(2,0) = 4.0*pi*a2 / s;
     highC(1,0) = 4.0*pi*a3 / s;
     highC(3,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(-2.0*a2*a3*fTilde_st);
       highC(2,0) += (1.0/s)*(a2*a2*(X_2_st-(C_d_2+C_A_2)/4.0*fTilde_st) + 
 			     2.0*(a1*a2*Yi*Yf+a2*a3*C_F_3)*W - 
 			     2.0*a1*a2*Yi*Yf*fTilde_st);
       highC(1,0) += (1.0/s)*(2.0*(a1*a3*Yi*Yf+a2*a3*C_F_2)*W - 
 			     2.0*a1*a3*Yi*Yf*fTilde_st);
       highC(3,0) += (1.0/s)*(-1.0*(a2*a2*C_1_2+a1*a1*Yi*Yi*Yf*Yf)*fTilde_st + 
 			     a1*a1*Yi*Yf*PI1 + 
 			     2.0*(a1*a3*Yi*Yf*C_F_3+a1*a2*Yi*Yf*C_F_2+a1*a1*Yi*Yi*Yf*Yf)*W);
       if (includeAlphaS2) {
 	highC(1,0) += (1.0/s)*(a3*a3*(X_3_st-(C_d_3+C_A_3)/4.0*fTilde_st));
 	
 	highC(3,0) += (1.0/s)*(-1.0*(a3*a3*C_1_3)*fTilde_st);
       }
     }
     break;
      
   case QQQQiden:
     {    
       Process parentCase = QQQQ;
       boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 	highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 	highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st(4,1);         
       boost::numeric::ublas::matrix<complex<InvEnergy2> > highCs_ts_2x2(2,2);
       highCs_ts_2x2(0,0) = highCs_ts(0,0);
       highCs_ts_2x2(1,0) = highCs_ts(1,0);
       highCs_ts_2x2(0,1) = highCs_ts(2,0);
       highCs_ts_2x2(1,1) = highCs_ts(3,0);
       boost::numeric::ublas::matrix<Complex> temp(2,2);
       temp = boost::numeric::ublas::trans(M_N(3));
       boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st_2x2(2,2),temp2(2,2);
       axpy_prod_local(highCs_ts_2x2,temp,temp2);
       axpy_prod_local(M_N(2),temp2,highCt_st_2x2);
       highCt_st(0,0) = highCt_st_2x2(0,0);
       highCt_st(1,0) = highCt_st_2x2(1,0);
       highCt_st(2,0) = highCt_st_2x2(0,1);
       highCt_st(3,0) = highCt_st_2x2(1,1);
       highC = highCs_st + highCt_st;
     }
     break;
   case QtQtQQ:
     {
       highC.resize(4,1);
       Process parentCase = QQQQ;
       highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       double Y = 1.0/6.0; // Hypercharge of the non-3rd-gen doublet (still a quark doublet).
       if (order >= 1) {
 	highC(2,0) += y_t*y_t*a2/(4.0*pi*s)*(3.0/2.0-0.5*Ls);
 	highC(1,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0/2.0-0.5*Ls);
 	highC(3,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(-5.0/12.0-1.0/12.0*Ls);
       }
     }
     break;
   case QQUU:
     Yi = 1./6.; Yf = 2./3.;
     highC.resize(2,1);
     highC(0,0) = 4.0*pi*a3 / s;
     highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*((a1*a3*(Yi*Yi+Yf*Yf)+a2*a3*C_F_2)*W + 
 			     2.0*a1*a3*Yi*Yf*fTilde_su);
       highC(1,0) += (1.0/s)*(a1*a1*Yi*Yf*PI1 + 
 			     (a1*a2*Yi*Yf*C_F_2+2.0*a1*a3*Yi*Yf*C_F_3+
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
       if (includeAlphaS2) {
       	highC(0,0) += (1.0/s)*(a3*a3*(X_3_su+(C_d_3-C_A_3)/4.0*fTilde_su));
       	highC(1,0) += (1.0/s)*((a3*a3*C_1_3+a1*a1*Yi*Yi*Yf*Yf)*fTilde_su);
       }
     }
     break;
   case QtQtUU:
     {
       highC.resize(2,1);
       Process parentCase = QQUU;
       highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       double Y = 2./3.;
       if (order >= 1) {
 	highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0/2.0-0.5*Ls);
 	highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(-5.0/12.0-1.0/12.0*Ls);
       }
     }
     break;         
   case QQtRtR:
     {
       highC.resize(2,1);
       Process parentCase = QQUU;
       highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       double Y = 1.0/6.0;
       if (order >= 1) {
 	highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0-Ls);
 	highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(5.0/3.0-2.0/3.0*Ls);
       }
     }
     break;
   case QQDD:
     Yi = 1./6.; Yf = -1./3.;
     highC.resize(2,1);
     highC(0,0) = 4.0*pi*a3 / s;
     highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*((a1*a3*(Yi*Yi+Yf*Yf)+a2*a3*C_F_2)*W + 
 			     2.0*a1*a3*Yi*Yf*fTilde_su);
       highC(1,0) += (1.0/s)*(a1*a1*Yi*Yf*PI1 + 
 			     (a1*a2*Yi*Yf*C_F_2+2.0*a1*a3*Yi*Yf*C_F_3+
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
       if (includeAlphaS2) {
       	highC(0,0) += (1.0/s)*(a3*a3*(X_3_su+(C_d_3-C_A_3)/4.0*fTilde_su));
       	highC(1,0) += (1.0/s)*((a3*a3*C_1_3+a1*a1*Yi*Yi*Yf*Yf)*fTilde_su);
       }
     }
     break;
   case QtQtDD:
     {
       highC.resize(2,1);
       Process parentCase = QQDD;
       highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       double Y = -1./3.;
       if (order >= 1) {
 	highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0/2.0-0.5*Ls);
 	highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(-5.0/12.0-1.0/12.0*Ls);
       }
     }
     break;
   case QQLL:
     Yi = 1./6.; Yf = -1./2.;
     highC.resize(2,1);
     highC(0,0) = 4.0*pi*a2 / s;
     highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(a2*a2*(X_2_st-(C_d_2+C_A_2)/4.0*fTilde_st) + 
 			     (a2*a3*C_F_3 + a1*a2*(Yi*Yi+Yf*Yf))*W - 
 			     2.0*a1*a2*Yi*Yf*fTilde_st);
       highC(1,0) += (1.0/s)*(-1.0*(a2*a2*C_1_2+a1*a1*Yi*Yi*Yf*Yf)*fTilde_st + 
 			     a1*a1*Yi*Yf*PI1 + 
 			     (a1*a3*Yi*Yf*C_F_3+2.0*a1*a2*Yi*Yf*C_F_2+
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
     }
     break;
   case QQEE:
     Yi = 1./6.; Yf = -1.;
     highC.resize(1,1);
     highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(a1*a1*Yi*Yi*Yf*Yf*fTilde_su + a1*a1*Yi*Yf*PI1 + 
 			     (a1*a3*Yi*Yf*C_F_3 + a1*a2*Yi*Yf*C_F_2 + 
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
     }
     break;
   case UUUU:
     Yi = Yf = 2./3.;
     highC.resize(2,1);
     highC(0,0) = 4.0*pi*a3 / s;
     highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(-2.0*a1*a3*Yi*Yf*fTilde_st +
 			     a1*a3*(Yi*Yi+Yf*Yf)*W);
       highC(1,0) += (1.0/s)*(-1.0*(a1*a1*Yi*Yi*Yf*Yf)*fTilde_st + 
 			     a1*a1*Yi*Yf*PI1 + 
 			     (2.0*a1*a3*Yi*Yf*C_F_3+a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
       if (includeAlphaS2) {
          highC(0,0) += (1.0/s)*(a3*a3*(X_3_st-(C_d_3+C_A_3)/4.0*fTilde_st));
          highC(1,0) += (1.0/s)*(-1.0*(a3*a3*C_1_3)*fTilde_st);
       }
     }
     break;
   case UUUUiden:
     {
       Process parentCase = UUUU;
       boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 	highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> >
 	highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st;
 	axpy_prod_local(M_N(3),highCs_ts,highCt_st);
       highC = highCs_st + highCt_st;
     }
     break;
   case tRtRUU:
     {
       highC.resize(2,1);
       Process parentCase = UUUU;
       highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       double Y = 2./3.;
       if (order >= 1) {
 	highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0-Ls);
 	highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(5.0/3.0-2.0/3.0*Ls);
       }
     }
     break;
   case UUDD:
     Yi = 2./3.; Yf = -1./3.;
     highC.resize(2,1);
     highC(0,0) = 4.0*pi*a3 / s;
     highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(-2.0*a1*a3*Yi*Yf*fTilde_st +
 			     a1*a3*(Yi*Yi+Yf*Yf)*W);
       highC(1,0) += (1.0/s)*(-1.0*(a1*a1*Yi*Yi*Yf*Yf)*fTilde_st + 
 			     a1*a1*Yi*Yf*PI1 + 
 			     (2.0*a1*a3*Yi*Yf*C_F_3+a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
       if (includeAlphaS2) {
 	highC(0,0) += (1.0/s)*(a3*a3*(X_3_st-(C_d_3+C_A_3)/4.0*fTilde_st));
 	highC(1,0) += (1.0/s)*(-1.0*(a3*a3*C_1_3)*fTilde_st);
       }
     }
     break;
   case tRtRDD:
     {
       highC.resize(2,1);
       Process parentCase = UUDD;
       highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       double Y = -1./3.;
       if (order >= 1) {
 	highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0-Ls);
 	highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(5.0/3.0-2.0/3.0*Ls);
       }
     }
     break;
   case UULL:
     Yi = 2./3.; Yf = -0.5;
     highC.resize(1,1);
     highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(a1*a1*Yi*Yi*Yf*Yf*fTilde_su + a1*a1*Yi*Yf*PI1 + 
 			     (a1*a3*Yi*Yf*C_F_3 + a1*a2*Yi*Yf*C_F_2 + 
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
     }
     break;
   case UUEE:
     Yi = 2./3.; Yf = -1.;
     highC.resize(1,1);
     highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(-1.0*a1*a1*Yi*Yi*Yf*Yf*fTilde_st + a1*a1*Yi*Yf*PI1 + 
 			     (a1*a3*Yi*Yf*C_F_3 + 
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
     }
     break;
   case DDDD:
     Yi = Yf = -1./3.;
     highC.resize(2,1);
     highC(0,0) = 4.0*pi*a3 / s;
     highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(-2.0*a1*a3*Yi*Yf*fTilde_st +
 			     a1*a3*(Yi*Yi+Yf*Yf)*W);
       highC(1,0) += (1.0/s)*(-1.0*(a1*a1*Yi*Yi*Yf*Yf)*fTilde_st + 
 			     a1*a1*Yi*Yf*PI1 + 
 			     (2.0*a1*a3*Yi*Yf*C_F_3+a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
       if (includeAlphaS2) {
          highC(0,0) += (1.0/s)*(a3*a3*(X_3_st-(C_d_3+C_A_3)/4.0*fTilde_st));
          highC(1,0) += (1.0/s)*(-1.0*(a3*a3*C_1_3)*fTilde_st);
       }
     }
     break;
   case DDDDiden:
     {
       Process parentCase = DDDD;
       boost::numeric::ublas::matrix<complex<InvEnergy2> >
 	highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> >
 	highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st;
       axpy_prod_local(M_N(3),highCs_ts,highCt_st);
       highC = highCs_st + highCt_st;
     }
     break;
   case DDLL:
     Yi = -1./3.; Yf = -0.5;
     highC.resize(1,1);
     highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(a1*a1*Yi*Yi*Yf*Yf*fTilde_su + a1*a1*Yi*Yf*PI1 + 
 			     (a1*a3*Yi*Yf*C_F_3 + a1*a2*Yi*Yf*C_F_2 + 
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
     }
     break;         
   case DDEE:
     Yi = -1./3.; Yf = -1.;
     highC.resize(1,1);
     highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(-1.0*a1*a1*Yi*Yi*Yf*Yf*fTilde_st + a1*a1*Yi*Yf*PI1 + 
 			     (a1*a3*Yi*Yf*C_F_3 + 
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
     }
     break;         
   case LLLL:
     Yi = Yf = -0.5;
     highC.resize(2,1);
     highC(0,0) = 4.0*pi*a2 / s;
     highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(a2*a2*(X_2_st-(C_d_2+C_A_2)/4.0*fTilde_st) + 
 			     2.0*a1*a2*Yi*Yf*W -
 			     2.0*a1*a2*Yi*Yf*fTilde_st);
       highC(1,0) += (1.0/s)*(-1.0*(a2*a2*C_1_2+a1*a1*Yi*Yi*Yf*Yf)*fTilde_st + 
 			     a1*a1*Yi*Yf*PI1 + 
 			     2.0*(a1*a2*Yi*Yf*C_F_2+a1*a1*Yi*Yi*Yf*Yf)*W);
     }
     break;
   case LLLLiden:
     {
       Process parentCase = LLLL;
       boost::numeric::ublas::matrix<complex<InvEnergy2> >
   	highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> >
   	highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st;
       axpy_prod_local(M_N(2), highCs_ts, highCt_st);
       highC = highCs_st + highCt_st;
     }
     break;
   case LLEE:
     Yi = -0.5; Yf = -1.;
     highC.resize(1,1);
     highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(a1*a1*Yi*Yi*Yf*Yf*fTilde_su + a1*a1*Yi*Yf*PI1 + 
 			     (a1*a2*Yi*Yf*C_F_2 +  
 			      a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
     }
     break;
   case EEEE:
     Yi = Yf = -1.;
     highC.resize(1,1);
     highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
     if (order >= 1) {
       highC(0,0) += (1.0/s)*(-1.0*a1*a1*Yi*Yi*Yf*Yf*fTilde_st + a1*a1*Yi*Yf*PI1 + 
 			     2.0*a1*a1*Yi*Yi*Yf*Yf*W);
     }
     break;
   case EEEEiden: 
     {
       Process parentCase = EEEE;
       boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 	highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 	highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
       boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 	highCt_st = highCs_ts;
       highC = highCs_st + highCt_st;
     }
     break;
   default:
     assert(false);
   }
   return highC;
 }
 
 boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 HighEnergyMatching::Spin1HighMatching(Energy highScale, 
 				      Energy2 s, Energy2 t, Energy2 u,
 				      EWProcess::Process process,
 				      bool oneLoop, bool includeAlphaS2) {
   using Constants::pi;
    
   unsigned int order = !oneLoop ? 0 : 1;
   // (If crossed graphs, swap s and t here)
   Complex L_s = MinusLog(-s/(highScale*highScale));
   Complex L_t = MinusLog(-t/(highScale*highScale));
   Complex L_u = MinusLog(-u/(highScale*highScale));
   Complex L_s2 = L_s*L_s;
   Complex L_t2 = L_t*L_t;
   Complex L_u2 = L_u*L_u;
    
   // Tree-Level:
   // Topology types defined. Note each of these is a vector of 5 entries. They are the coefficients 
   // for the dirac structures M_0, M_1, M_4, M_5, and M_6 for vector boson production.
   boost::numeric::ublas::vector<complex<InvEnergy2> > R1(5);
   for(unsigned int ix=0;ix<5;++ix) R1[ix] = ZERO;
   R1[0] = -1.0/t;
   R1[1] = -2.0/t;
   R1[2] = ZERO;
   R1[3] = ZERO;
   R1[4] = ZERO;
   boost::numeric::ublas::vector<complex<InvEnergy2> > R1_bar(5);
   for(unsigned int ix=0;ix<5;++ix) R1_bar[ix] = ZERO;
   R1_bar[0] = -1.0/u;
   boost::numeric::ublas::vector<complex<InvEnergy2> > R2(5);
   for(unsigned int ix=0;ix<5;++ix) R2[ix] = ZERO;
   R2[1] = -1.0/s*2.0;
   // Topologies T1:
   boost::numeric::ublas::vector<complex<InvEnergy2> > T1a(5);
   for(unsigned int ix=0;ix<5;++ix) T1a[ix] = ZERO;
   T1a[0] = 1.0/(t*u)*(-3.0*t*L_s2-(s+4.0*t)*L_t2+2.0*(s+4.0*t)*L_s*L_t+2.0*u*L_t-
                       pi*pi*(7.0/6.0*s+25.0/6.0*t)-4.0*u);
   T1a[1] = 1.0/(u*u*t*s)*(0.5*t*(9.0*s*s+14.0*s*t+7.0*t*t)*L_s2+s*(2.0*s+t)*(s+2.0*t)*L_t2-
                           2.0*(2.0*s*s*s+9.0*s*s*t+10.0*s*t*t+4.0*t*t*t)*L_s*L_t-
                           2.0*t*t*u*L_s-2.0*u*s*(2.0*s+3.0*t)*L_t+
                           pi*pi*(7.0/3.0*s*s*s+125.0/12.0*s*s*t+71.0/6.0*s*t*t+
                                  19.0/4.0*t*t*t)-
                           8.0*s*s*s-20.0*s*s*t-16.0*s*t*t-4.0*t*t*t);
   T1a[2] = 1.0/(t*u*u)*(-t*(3.0*s+4.0*t)*L_s2-(s*s+5.0*s*t+5.0*t*t)*L_t2+
                         2.0*t*(3.0*s+4.0*t)*L_s*L_t+2.0*u*t*(2.0*s+t)*L_s/s-
                         2.0*u*t*L_t+pi*pi*(s*s/6.0-8.0/3.0*s*t-23.0/6.0*t*t)+
                         4.0*t*t*t/s+4.0*s*t+8.0*t*t);
   T1a[3] = T1a[2];
   T1a[4] = GeV2/(t*u*u*u)*(-4.0*t*(s+2.0*t)*(L_s-L_t)*(L_s-L_t)+
 			   4.0*u*(3.0*s+5.0*t)*(L_s-L_t)-4.0*pi*pi*t*(s+2.0*t)-4.0*u*u);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T1b(5);
   for(unsigned int ix=0;ix<5;++ix) T1b[ix] = ZERO;
   T1b[0] = 1.0/(t*u*s*s)*(-s*t*(2.0*s+3.0*t)*L_u2+s*u*(s+3.0*t)*L_t2+
 			  2.0*s*(s*s+3.0*s*t+3.0*t*t)*L_u*L_t+s*s*t*L_u+s*s*u*L_t-
 			  pi*pi*(7.0/6.0*s*s*s+3.0*s*s*t+3.0*s*t*t)+2.0*s*s*s);
   T1b[1] = 1.0/(t*s*s*u)*(3.0*s*t*u*L_u2+s*u*(2.0*s+3.0*t)*L_t2-
 			  2.0*s*u*(2.0*s+3.0*t)*L_u*L_t+2.0*s*s*u*L_t-
 			  pi*pi*(7.0/3.0*s*s*s+16.0/3.0*s*s*t+3.0*s*t*t)+
 			  4.0*s*s*s+4.0*s*s*t);
   T1b[2] = 1.0/(t*u*s*s)*(-3.0*s*t*u*(L_u-L_t)*(L_u-L_t)+4.0*s*s*t*L_u+4.0*s*s*u*L_t+
 			  pi*pi*(3.0*s*s*t+3.0*s*t*t)+8.0*s*s*s);
   T1b[3] = 1.0/(t*u*s*s)*(s*t*(2.0*s+3.0*t)*L_u2-s*u*(s+3.0*t)*L_t2+6.0*s*t*u*L_u*L_t+
 			  pi*pi*(-1.0/6.0*s*s*s+3.0*s*s*t+3.0*s*t*t));
   T1b[4] = 12.0*GeV2/(t*u)*(L_t-L_u);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T1c(5);
   for(unsigned int ix=0;ix<5;++ix) T1c[ix] = ZERO;
   T1c[0] = 1.0/t*(2.0*L_s*L_t-7.0*pi*pi/6.0-L_t2);
   T1c[1] = 1.0/(t*u*u)*(s*t*L_s2-(2.0*s*s+3.0*s*t+2.0*t*t)*L_t2+
 			2.0*(2.0*s*s+3.0*s*t+2.0*t*t)*L_s*L_t+2.0*t*u*(L_s-L_t)-
 			pi*pi*(7.0/3.0*s*s+11.0/3.0*s*t+7.0/3.0*t*t));
   T1c[2] = 1.0/(t*u*u)*(t*(3.0*s+2.0*t)*(L_s-L_t)*(L_s-L_t)+2.0*u*t*L_s+
 			2.0*u*(2.0*s+t)*L_t+pi*pi*t*(3.0*s+2.0*t)+8.0*u*u);
   T1c[3] = T1c[2];
   T1c[4] = GeV2/(t*u*u*u)*(4.0*t*(2.0*s+t)*(L_s-L_t)*(L_s-L_t)-4.0*u*(3.0*s+t)*(L_s-L_t)+
 			   4.0*pi*pi*t*(2.0*s+t)-4.0*u*u);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T1d(5);
   for(unsigned int ix=0;ix<5;++ix) T1d[ix] = ZERO;
   T1d[2] = 1.0/s*(-2.0*L_s+4.0);
   T1d[3] = T1d[2];
   boost::numeric::ublas::vector<complex<InvEnergy2> > T1a_bar(5);
   for(unsigned int ix=0;ix<5;++ix) T1a_bar[ix] = ZERO;
   T1a_bar[0] = 1.0/(u*t)*(-3.0*u*L_s2-(s+4.0*u)*L_u2+2.0*(s+4.0*u)*L_s*L_u+2.0*t*L_u-
 			  pi*pi*(7.0/6.0*s+25.0/6.0*u)-4.0*t);
   T1a_bar[1] = 2.0*T1a_bar[0] - 
     1.0/(t*t*u*s)*(0.5*u*(9.0*s*s+14.0*s*u+7.0*u*u)*L_s2+s*(2.0*s+u)*(s+2.0*u)*L_u2-
 		   2.0*(2.0*s*s*s+9.0*s*s*u+10.0*s*u*u+4.0*u*u*u)*L_s*L_u-
 		   2.0*u*u*t*L_s-2.0*t*s*(2.0*s+3.0*u)*L_u+
 		   pi*pi*(7.0/3.0*s*s*s+125.0/12.0*s*s*u+71.0/6.0*s*u*u+
 			  19.0/4.0*u*u*u)-
 		   8.0*s*s*s-20.0*s*s*u-16.0*s*u*u-4.0*u*u*u);
   T1a_bar[2] = 1.0/(u*t*t)*(-u*(3.0*s+4.0*u)*L_s2-(s*s+5.0*s*u+5.0*u*u)*L_u2+
 			    2.0*u*(3.0*s+4.0*u)*L_s*L_u+2.0*t*u*(2.0*s+u)*L_s/s-
 			    2.0*t*u*L_u+pi*pi*(s*s/6.0-8.0/3.0*s*u-23.0/6.0*u*u)+
 			    4.0*u*u*u/s+4.0*s*u+8.0*u*u);
   T1a_bar[3] = T1a_bar[2];
   T1a_bar[4] = -GeV2/(u*t*t*t)*(-4.0*u*(s+2.0*u)*(L_s-L_u)*(L_s-L_u)+
 				4.0*t*(3.0*s+5.0*u)*(L_s-L_u)-4.0*pi*pi*u*(s+2.0*u)-4.0*t*t);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T1b_bar(5);
   for(unsigned int ix=0;ix<5;++ix) T1b_bar[ix] = ZERO;
   T1b_bar[0] = 1.0/(u*t*s*s)*(-s*u*(2.0*s+3.0*u)*L_t2+s*t*(s+3.0*u)*L_u2+
 			      2.0*s*(s*s+3.0*s*u+3.0*u*u)*L_t*L_u+
 			      s*s*u*L_t+s*s*t*L_u-
 			      pi*pi*(7.0/6.0*s*s*s+3.0*s*s*u+3.0*s*u*u)+2.0*s*s*s);
   T1b_bar[1] = 2.0*T1b_bar[0] - 
     1.0/(u*s*s*t)*(3.0*s*u*t*L_t2+s*t*(2.0*s+3.0*u)*L_u2-
 		   2.0*s*t*(2.0*s+3.0*u)*L_t*L_u+2.0*s*s*t*L_u-
 		   pi*pi*(7.0/3.0*s*s*s+16.0/3.0*s*s*u+3.0*s*u*u)+
 		   4.0*s*s*s+4.0*s*s*u);
   T1b_bar[3] = 1.0/(u*t*s*s)*(-3.0*s*u*t*(L_t-L_u)*(L_t-L_u)+4.0*s*s*u*L_t+4.0*s*s*t*L_u+
 			      pi*pi*(3.0*s*s*u+3.0*s*u*u)+8.0*s*s*s);
   T1b_bar[2] = 1.0/(u*t*s*s)*(s*u*(2.0*s+3.0*u)*L_t2-s*t*(s+3.0*u)*L_u2+6.0*s*u*t*L_t*L_u+
 			      pi*pi*(-1.0/6.0*s*s*s+3.0*s*s*u+3.0*s*u*u));
   T1b_bar[4] = -12.0*GeV2/(u*t)*(L_u-L_t);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T1c_bar(5);
   for(unsigned int ix=0;ix<5;++ix) T1c_bar[ix] = ZERO;
   T1c_bar[0] = 1.0/u*(2.0*L_s*L_u-7.0*pi*pi/6.0-L_u2);
   T1c_bar[1] = 2.0*T1c_bar[0] - 
     1.0/(u*t*t)*(s*u*L_s2-(2.0*s*s+3.0*s*u+2.0*u*u)*L_u2+
 		 2.0*(2.0*s*s+3.0*s*u+2.0*u*u)*L_s*L_u+2.0*u*t*(L_s-L_u)-
 		 pi*pi*(7.0/3.0*s*s+11.0/3.0*s*u+7.0/3.0*u*u));
   T1c_bar[2] = 1.0/(u*t*t)*(u*(3.0*s+2.0*u)*(L_s-L_u)*(L_s-L_u)+2.0*t*u*L_s+
                             2.0*t*(2.0*s+u)*L_u+pi*pi*u*(3.0*s+2.0*u)+8.0*t*t);
   T1c_bar[3] = T1c_bar[2];
   T1c_bar[4] = -GeV2/(u*t*t*t)*(4.0*u*(2.0*s+u)*(L_s-L_u)*(L_s-L_u)-4.0*t*(3.0*s+u)*(L_s-L_u)+
                                4.0*pi*pi*u*(2.0*s+u)-4.0*t*t);
   // Topologies T2:
   boost::numeric::ublas::vector<complex<InvEnergy2> > T2a(5);
   for(unsigned int ix=0;ix<5;++ix) T2a[ix] = ZERO;
   T2a[1] = 1.0/s*(L_s/6.0-11.0/18.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T2b(5);
   for(unsigned int ix=0;ix<5;++ix) T2b[ix] = ZERO;
   T2b[1] = 1.0/s*(-2.0/3.0*L_s+22.0/9.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T2c(5);
   for(unsigned int ix=0;ix<5;++ix) T2c[ix] = ZERO;
   T2c[1] = 1.0/s*(3.0/2.0*L_s2-17.0/2.0*L_s-pi*pi/4.0+95.0/6.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T2d(5);
   for(unsigned int ix=0;ix<5;++ix) T2d[ix] = ZERO;
   T2d[1] = 1.0/s*(-4.0/3.0*L_s+14.0/9.0);
   // Topologies T3:
   boost::numeric::ublas::vector<complex<InvEnergy2> > T3a(5);
   for(unsigned int ix=0;ix<5;++ix) T3a[ix] = ZERO;
   T3a[0] = 1.0/t*(L_t2-pi*pi/6.0);
   T3a[1] = 2.0*T3a[0];
   T3a[3] = 1.0/t*(-L_t2+pi*pi/6.0+2.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T3b(5);
   for(unsigned int ix=0;ix<5;++ix) T3b[ix] = ZERO;
   T3b[0] = 1.0/t*(-L_t+4.0);
   T3b[1] = 2.0*T3b[0];
   T3b[3] = 1.0/t*(4.0*L_t-10.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T3a_bar(5);
   for(unsigned int ix=0;ix<5;++ix) T3a_bar[ix] = ZERO;
   T3a_bar[0] = 1.0/u*(L_u2-pi*pi/6.0);
   T3a_bar[2] = 1.0/u*(-L_u2+pi*pi/6.0+2.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T3b_bar(5);
   for(unsigned int ix=0;ix<5;++ix) T3b_bar[ix] = ZERO;
   T3b_bar[0] = 1.0/u*(-L_u+4.0);
   T3b_bar[2] = 1.0/u*(4.0*L_u-10.0);
   // Topologies T4:
   boost::numeric::ublas::vector<complex<InvEnergy2> > T4a(5);
   for(unsigned int ix=0;ix<5;++ix) T4a[ix] = ZERO;
   T4a[0] = 1.0/t*(L_t2-pi*pi/6.0);
   T4a[1] = 2.0*T4a[0];
   T4a[2] = 1.0/t*(-L_t2+pi*pi/6.0+2.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T4b(5);
   for(unsigned int ix=0;ix<5;++ix) T4b[ix] = ZERO;
   T4b[0] = 1.0/t*(-L_t+4.0);
   T4b[1] = 2.0*T4b[0];
   T4b[2] = 1.0/t*(4.0*L_t-10.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T4a_bar(5);
   for(unsigned int ix=0;ix<5;++ix) T4a_bar[ix] = ZERO;
   T4a_bar[0] = 1.0/u*(L_u2-pi*pi/6.0);
   T4a_bar[3] = 1.0/u*(-L_u2+pi*pi/6.0+2.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T4b_bar(5);
   for(unsigned int ix=0;ix<5;++ix) T4b_bar[ix] = ZERO;
   T4b_bar[0] = 1.0/u*(-L_u+4.0);
   T4b_bar[3] = 1.0/u*(4.0*L_u-10.0);
   // Topologies T5:
   boost::numeric::ublas::vector<complex<InvEnergy2> > T5a(5);
   for(unsigned int ix=0;ix<5;++ix) T5a[ix] = ZERO;
   T5a[1] = 1.0/s*(2.0*L_s-4.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T5b(5);
   for(unsigned int ix=0;ix<5;++ix) T5b[ix] = ZERO;
   T5b[1] = 1.0/s*(-2.0*L_s2+6.0*L_s-16.0+pi*pi/3.0);
   // Topologies T6:
   boost::numeric::ublas::vector<complex<InvEnergy2> > T6a(5);
   for(unsigned int ix=0;ix<5;++ix) T6a[ix] = ZERO;
   T6a[1] = 1.0/s*(-19.0/6.0*L_s+58.0/9.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T6b(5);
   for(unsigned int ix=0;ix<5;++ix) T6b[ix] = ZERO;
   T6b[1] = 1.0/s*(-1.0/6.0*L_s+4.0/9.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T6c(5);
   for(unsigned int ix=0;ix<5;++ix) T6c[ix] = ZERO;
   T6c[1] = 1.0/s*(2.0/3.0*L_s-16.0/9.0);
   boost::numeric::ublas::vector<complex<InvEnergy2> > T6d(5);
   for(unsigned int ix=0;ix<5;++ix) T6d[ix] = ZERO;
   T6d[1] = 1.0/s*(4.0/3.0*L_s-20.0/9.0);
   // Topology T7:
   boost::numeric::ublas::vector<complex<InvEnergy2> > T7(5);
   for(unsigned int ix=0;ix<5;++ix) T7[ix] = ZERO;
   T7[0] = 1.0/t*(-L_t+1.0);
   T7[1] = 2.0*T7[0];
   boost::numeric::ublas::vector<complex<InvEnergy2> > T7_bar(5);
   for(unsigned int ix=0;ix<5;++ix) T7_bar[ix] = ZERO;
   T7_bar[0] = 1.0/u*(-L_u+1.0);
   // Group Theory Factors / SM parameters needed for matrix elements:
   double a1 = ElectroWeakReweighter::coupling()->a1(highScale);
   double a2 = ElectroWeakReweighter::coupling()->a2(highScale);
   double a3 = ElectroWeakReweighter::coupling()->a3(highScale);
   double y_t = ElectroWeakReweighter::coupling()->y_t(highScale);
   // Traces over complex scalars and weyl fermions.
   double T_CS_3 = 0.0;
   double T_CS_2 = 0.5;
   //double T_CS_1 = 0.5;
   double T_WF_3 = 2.0*3.0;
   double T_WF_2 = 2.0*3.0;
   //double T_WF_1 = 10.0/3.0*3.0;
   double C_A_3 = 3.0;
   double C_A_2 = 2.0;
   double C_A_1 = 0.0;
   double C_F_3 = 4.0/3.0;
   double C_F_2 = 3.0/4.0;
   double C_F_1 = 1.0;
   // This is the coefficient of the delta term in G_TT
   double G_TT = 0.5;
   // This is the coeffidient of d^ABC in G_TT (non-zero for SU(3))
   double G_TT_3_D = 0.25*C_A_3; 
   double G_f = 1.0;
   // Factors TBD after fermion helicity is specified:
   double Y_Q(0.), G_Plus_U1(0.);
   double G_Plus_SU2 = 0.25;
   double G_Plus_SU3 = 1./6.;
   double G_Plus_SU3_D = 0.5;
   double Lambda_Q(0.);
   // the _s and _ew are the alpha3 and alpha1/2 parts of Lambda_Q
   double Lambda_Q_s(0.); 
   double Lambda_Q_ew(0.);
   double rho12(0.), rho13(0.);
   double rho23 = sqrt(a2*a3);
   double tRorQ = 1.0;
   boost::numeric::ublas::matrix<complex<InvEnergy2> > highC(1,1);
   switch (process) {
   case QQWW: case LLWW:
     {
       // Finish Group Theory Factors:
       if (process==QQWW) {
 	Y_Q = 1.0/6.0;
 	G_Plus_U1 = Y_Q*Y_Q;
 	Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
 	rho12 = Y_Q*sqrt(a1*a2);
       }
       else if (process==LLWW) {
 	Y_Q = -1.0/2.0;
 	G_Plus_U1 = Y_Q*Y_Q;
 	Lambda_Q = C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
 	rho12 = Y_Q*sqrt(a1*a2);
       }
       highC.resize(5,5);
       for (int i=0; i<5; i++) {
 	highC(0,i) = G_Plus_SU2*(4.0*pi*a2)*(R1[i]+R1_bar[i]);
 	highC(1,i) = G_f*4.0*pi*a2*(-0.5*R1[i]+0.5*R1_bar[i]-R2[i]);
 	highC(2,i) = rho12*4.0*pi*(R1[i]+R1_bar[i]);
 	highC(3,i) = rho12*4.0*pi*(R1[i]+R1_bar[i]);
 	highC(4,i) = G_Plus_U1*(4.0*pi*a1)*(R1[i]+R1_bar[i]);
 	if (order>=1) {
 	  highC(0,i) += G_Plus_SU2*((-0.5*a2*a2*C_A_2)*(T1b[i]+T1b_bar[i])+
 				    a2*(Lambda_Q-a2*C_A_2)*(T1c[i]+T1c_bar[i])+
 				    0.5*a2*a2*C_A_2*(T3a[i]+T3a_bar[i])+
 				    a2*(Lambda_Q-0.5*a2*C_A_2)*(T3b[i]+T3b_bar[i])+
 				    0.5*a2*a2*C_A_2*(T4a[i]+T4a_bar[i])+
 				    a2*(Lambda_Q-0.5*a2*C_A_2)*(T4b[i]+T4b_bar[i])+
 				    a2*Lambda_Q*(T7[i]+T7_bar[i])) + 
 	    G_TT*(-a2*a2*(T1a[i]+T1a_bar[i])+a2*a2*(T1b[i]+T1b_bar[i])+
 		  a2*a2*(T1c[i]+T1c_bar[i])+2.0*a2*a2*T1d[i]);
 	  highC(1,i) += G_f*(0.25*a2*a2*C_A_2*(T1a[i]-T1a_bar[i])+
 			     a2*(0.25*a2*C_A_2-0.5*Lambda_Q)*(T1c[i]-T1c_bar[i])+
 			     0.5*a2*a2*C_A_2*T2a[i]+a2*a2*T_CS_2*T2b[i]-
 			     0.5*a2*a2*C_A_2*T2c[i]+a2*a2*T_WF_2*T2d[i]-
 			     0.25*a2*a2*C_A_2*(T3a[i]-T3a_bar[i])-
 			     0.5*a2*(Lambda_Q-0.5*a2*C_A_2)*(T3b[i]-T3b_bar[i])-
 			     0.25*a2*a2*C_A_2*(T4a[i]-T4a_bar[i])-
 			     0.5*a2*(Lambda_Q-0.5*a2*C_A_2)*(T4b[i]-T4b_bar[i])+
 			     0.5*a2*a2*C_A_2*T5a[i]+a2*(Lambda_Q-0.5*a2*C_A_2)*T5b[i]+
 			     a2*a2*C_A_2*T6a[i]+a2*a2*C_A_2*T6b[i]+
 			     a2*a2*T_CS_2*T6c[i]+a2*a2*T_WF_2*T6d[i]-
 			     0.5*a2*Lambda_Q*(T7[i]-T7_bar[i]));
 	  highC(2,i) += rho12*(-0.5*a1*C_A_1*T1b[i]-0.5*a2*C_A_2*T1b_bar[i]+
 			       (Lambda_Q-0.5*a1*C_A_1-0.5*a2*C_A_2)*(T1c[i]+T1c_bar[i])+
 			       0.5*a1*C_A_1*T3a[i]+0.5*a2*C_A_2*T3a_bar[i]+
 			       (Lambda_Q-0.5*a1*C_A_1)*T3b[i]+(Lambda_Q-0.5*a2*C_A_2)*T3b_bar[i]+
 			       0.5*a2*C_A_2*T4a[i]+0.5*a1*C_A_1*T4a_bar[i]+
 			       (Lambda_Q-0.5*a2*C_A_2)*T4b[i]+(Lambda_Q-0.5*a1*C_A_1)*T4b_bar[i]+
 			       Lambda_Q*(T7[i]+T7_bar[i]));
 	  highC(3,i) += rho12*(-0.5*a2*C_A_2*T1b[i]-0.5*a1*C_A_1*T1b_bar[i]+
 			       (Lambda_Q-0.5*a2*C_A_2-0.5*a1*C_A_1)*(T1c[i]+T1c_bar[i])+
 			       0.5*a2*C_A_2*T3a[i]+0.5*a1*C_A_1*T3a_bar[i]+
 			       (Lambda_Q-0.5*a2*C_A_2)*T3b[i]+(Lambda_Q-0.5*a1*C_A_1)*T3b_bar[i]+
 			       0.5*a1*C_A_1*T4a[i]+0.5*a2*C_A_2*T4a_bar[i]+
 			       (Lambda_Q-0.5*a1*C_A_1)*T4b[i]+(Lambda_Q-0.5*a2*C_A_2)*T4b_bar[i]+
 			       Lambda_Q*(T7[i]+T7_bar[i]));
 	  highC(4,i) += G_Plus_U1*((-0.5*a1*a1*C_A_1)*(T1b[i]+T1b_bar[i])+
 				   a1*(Lambda_Q-a1*C_A_1)*(T1c[i]+T1c_bar[i])+
 				   0.5*a1*a1*C_A_1*(T3a[i]+T3a_bar[i])+
 				   a1*(Lambda_Q-0.5*a1*C_A_1)*(T3b[i]+T3b_bar[i])+
 				   0.5*a1*a1*C_A_1*(T4a[i]+T4a_bar[i])+
 				   a1*(Lambda_Q-0.5*a1*C_A_1)*(T4b[i]+T4b_bar[i])+
 				   a1*Lambda_Q*(T7[i]+T7_bar[i]));
 	}
       }
     }
     break;
   case UUBB: case DDBB: case EEBB: 
     {
       // Finish Group Theory Factors:
       if (process==UUBB) {
 	Y_Q = 2.0/3.0;
 	G_Plus_U1 = Y_Q*Y_Q;
 	Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
       }
       else if (process==DDBB) {
 	Y_Q = -1.0/3.0;
 	G_Plus_U1 = Y_Q*Y_Q;
 	Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
       }
       else if (process==EEBB) {
 	Y_Q = -1.0;
 	G_Plus_U1 = Y_Q*Y_Q;
 	Lambda_Q = Y_Q*Y_Q*C_F_1*a1;
       }
       highC.resize(1,5);
       for (int i=0; i<5; i++) {
 	highC(0,i) = G_Plus_U1*(4.0*pi*a1)*(R1[i]+R1_bar[i]);
 	if (order>=1) {
 	  highC(0,i) += G_Plus_U1*((-0.5*a1*a1*C_A_1)*(T1b[i]+T1b_bar[i])+
 				   a1*(Lambda_Q-a1*C_A_1)*(T1c[i]+T1c_bar[i])+
 				   0.5*a1*a1*C_A_1*(T3a[i]+T3a_bar[i])+
 				   a1*(Lambda_Q-0.5*a1*C_A_1)*(T3b[i]+T3b_bar[i])+
 				   0.5*a1*a1*C_A_1*(T4a[i]+T4a_bar[i])+
 				   a1*(Lambda_Q-0.5*a1*C_A_1)*(T4b[i]+T4b_bar[i])+
 				   a1*Lambda_Q*(T7[i]+T7_bar[i]));
 	}
       }
     }
     break;
   case QQWG:
     {
       // Finish Group Theory Factors:
       Y_Q = 1./6.;
       Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
       
       highC.resize(1,5);
       
       for (int i=0; i<5; i++) {
 	highC(0,i) = rho23*4.0*pi*(R1[i]+R1_bar[i]);
         
 	if (order>=1) {
 	  
 	  highC(0,i) += rho23*(-0.5*a3*C_A_3*T1b[i]-0.5*a2*C_A_2*T1b_bar[i]+
 			       (Lambda_Q-0.5*a3*C_A_3-0.5*a2*C_A_2)*(T1c[i]+T1c_bar[i])+
 			       0.5*a3*C_A_3*T3a[i]+0.5*a2*C_A_2*T3a_bar[i]+
 			       (Lambda_Q-0.5*a3*C_A_3)*T3b[i]+(Lambda_Q-0.5*a2*C_A_2)*T3b_bar[i]+
 			       0.5*a2*C_A_2*T4a[i]+0.5*a3*C_A_3*T4a_bar[i]+
 			       (Lambda_Q-0.5*a2*C_A_2)*T4b[i]+(Lambda_Q-0.5*a3*C_A_3)*T4b_bar[i]+
 			       Lambda_Q*(T7[i]+T7_bar[i]));
 	}
       }
     }
     break;
   case QQBG:
     {
       // Finish Group Theory Factors:
       Y_Q = 1.0/6.0;
       Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
       rho13 = Y_Q*sqrt(a1*a3);
       
       highC.resize(1,5);
       
       for (int i=0; i<5; i++) {
 	highC(0,i) = rho13*4.0*pi*(R1[i]+R1_bar[i]);
         
 	if (order>=1) {
 	  
 	  highC(0,i) += rho13*(-0.5*a3*C_A_3*T1b[i]-0.5*a1*C_A_1*T1b_bar[i]+
 			       (Lambda_Q-0.5*a3*C_A_3-0.5*a1*C_A_1)*(T1c[i]+T1c_bar[i])+
 			       0.5*a3*C_A_3*T3a[i]+0.5*a1*C_A_1*T3a_bar[i]+
 			       (Lambda_Q-0.5*a3*C_A_3)*T3b[i]+(Lambda_Q-0.5*a1*C_A_1)*T3b_bar[i]+
 			       0.5*a1*C_A_1*T4a[i]+0.5*a3*C_A_3*T4a_bar[i]+
 			       (Lambda_Q-0.5*a1*C_A_1)*T4b[i]+(Lambda_Q-0.5*a3*C_A_3)*T4b_bar[i]+
 			       Lambda_Q*(T7[i]+T7_bar[i]));
 	}
       }
     }
     break;
   case UUBG: case DDBG:
      {
         // Finish Group Theory Factors:
         if (process==UUBG) {
            Y_Q = 2.0/3.0;
            Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
            rho13 = Y_Q*sqrt(a1*a3);
         }
         else if (process==DDBG) {
            Y_Q = -1.0/3.0;
            Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
            rho13 = Y_Q*sqrt(a1*a3);
         }
      
         highC.resize(1,5);
      
         for (int i=0; i<5; i++) {
            highC(0,i) = rho13*4.0*pi*(R1[i]+R1_bar[i]);
         
            if (order>=1) {
            
               highC(0,i) += rho13*(-0.5*a3*C_A_3*T1b[i]-0.5*a1*C_A_1*T1b_bar[i]+
                                     (Lambda_Q-0.5*a3*C_A_3-0.5*a1*C_A_1)*(T1c[i]+T1c_bar[i])+
                                     0.5*a3*C_A_3*T3a[i]+0.5*a1*C_A_1*T3a_bar[i]+
                                     (Lambda_Q-0.5*a3*C_A_3)*T3b[i]+(Lambda_Q-0.5*a1*C_A_1)*T3b_bar[i]+
                                     0.5*a1*C_A_1*T4a[i]+0.5*a3*C_A_3*T4a_bar[i]+
                                     (Lambda_Q-0.5*a1*C_A_1)*T4b[i]+(Lambda_Q-0.5*a3*C_A_3)*T4b_bar[i]+
                                     Lambda_Q*(T7[i]+T7_bar[i]));
            }
         }
      }
     break;
   case QQGG: case QtQtGG: case UUGG:
   case tRtRGG: case DDGG:
     {
       // Finish Group Theory Factors:
       if (process==QQGG || process==QtQtGG) {
 	Y_Q = 1.0/6.0;
 	Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
 	Lambda_Q_s = C_F_3*a3;
 	Lambda_Q_ew = C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
       }
       else if (process==UUGG || process==tRtRGG) {
 	Y_Q = 2.0/3.0;
 	Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
 	Lambda_Q_s = C_F_3*a3;
 	Lambda_Q_ew = Y_Q*Y_Q*C_F_1*a1;
       }
       else if (process==DDGG || process==tRtRGG) {
 	Y_Q = -1.0/3.0;
 	Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
 	Lambda_Q_s = C_F_3*a3;
 	Lambda_Q_ew = Y_Q*Y_Q*C_F_1*a1;
       }
       
       highC.resize(3,5);
          
       for (int i=0; i<5; i++) {
 	highC(0,i) = G_Plus_SU3*(4.0*pi*a3)*(R1[i]+R1_bar[i]);
 	highC(1,i) = G_Plus_SU3_D*(4.0*pi*a3)*(R1[i]+R1_bar[i]);
 	highC(2,i) = G_f*4.0*pi*a3*(-0.5*R1[i]+0.5*R1_bar[i]-R2[i]);
         
 	if (order>=1) {
 	  highC(0,i) += G_Plus_SU3*(a3*(Lambda_Q_ew)*(T1c[i]+T1c_bar[i])+
 				    a3*(Lambda_Q_ew)*(T3b[i]+T3b_bar[i])+
 				    a3*(Lambda_Q_ew)*(T4b[i]+T4b_bar[i])+
 				    a3*Lambda_Q_ew*(T7[i]+T7_bar[i]));
 	  highC(1,i) += G_Plus_SU3_D*(a3*(Lambda_Q_ew)*(T1c[i]+T1c_bar[i])+
 				      a3*(Lambda_Q_ew)*(T3b[i]+T3b_bar[i])+
 				      a3*(Lambda_Q_ew)*(T4b[i]+T4b_bar[i])+
 				      a3*Lambda_Q_ew*(T7[i]+T7_bar[i]));
 	  highC(2,i) += G_f*(a3*(-0.5*Lambda_Q_ew)*(T1c[i]-T1c_bar[i])-
 			     0.5*a3*(Lambda_Q_ew)*(T3b[i]-T3b_bar[i])-
 			     0.5*a3*(Lambda_Q_ew)*(T4b[i]-T4b_bar[i])+
 			     a3*(Lambda_Q_ew)*T5b[i]-
 			     0.5*a3*Lambda_Q_ew*(T7[i]-T7_bar[i]));
 	  if (includeAlphaS2) {
 	    highC(0,i) += G_Plus_SU3*((-0.5*a3*a3*C_A_3)*(T1b[i]+T1b_bar[i])+
 				      a3*(Lambda_Q_s-a3*C_A_3)*(T1c[i]+T1c_bar[i])+
 				      0.5*a3*a3*C_A_3*(T3a[i]+T3a_bar[i])+
 				      a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T3b[i]+T3b_bar[i])+
 				      0.5*a3*a3*C_A_3*(T4a[i]+T4a_bar[i])+
 				      a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T4b[i]+T4b_bar[i])+
 				      a3*Lambda_Q_s*(T7[i]+T7_bar[i])) + 
 	      G_TT*(-a3*a3*(T1a[i]+T1a_bar[i])+a3*a3*(T1b[i]+T1b_bar[i])+
 		    a3*a3*(T1c[i]+T1c_bar[i])+2.0*a3*a3*T1d[i]);
 	    highC(1,i) += G_Plus_SU3_D*((-0.5*a3*a3*C_A_3)*(T1b[i]+T1b_bar[i])+
 					a3*(Lambda_Q_s-a3*C_A_3)*(T1c[i]+T1c_bar[i])+
 					0.5*a3*a3*C_A_3*(T3a[i]+T3a_bar[i])+
 					a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T3b[i]+T3b_bar[i])+
 					0.5*a3*a3*C_A_3*(T4a[i]+T4a_bar[i])+
 					a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T4b[i]+T4b_bar[i])+
 					a3*Lambda_Q_s*(T7[i]+T7_bar[i])) + 
 	      G_TT_3_D*(-a3*a3*(T1a[i]+T1a_bar[i])+a3*a3*(T1b[i]+T1b_bar[i])+
 			a3*a3*(T1c[i]+T1c_bar[i])+2.0*a3*a3*T1d[i]);
 	    highC(2,i) += G_f*(0.25*a3*a3*C_A_3*(T1a[i]-T1a_bar[i])+
 			       a3*(0.25*a3*C_A_3-0.5*Lambda_Q_s)*(T1c[i]-T1c_bar[i])+
 			       0.5*a3*a3*C_A_3*T2a[i]+a3*a3*T_CS_3*T2b[i]-
 			       0.5*a3*a3*C_A_3*T2c[i]+a3*a3*T_WF_3*T2d[i]-
 			       0.25*a3*a3*C_A_3*(T3a[i]-T3a_bar[i])-
 			       0.5*a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T3b[i]-T3b_bar[i])-
 			       0.25*a3*a3*C_A_3*(T4a[i]-T4a_bar[i])-
 			       0.5*a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T4b[i]-T4b_bar[i])+
 			       0.5*a3*a3*C_A_3*T5a[i]+a3*(Lambda_Q_s-0.5*a3*C_A_3)*T5b[i]+
 			       a3*a3*C_A_3*T6a[i]+a3*a3*C_A_3*T6b[i]+
 			       a3*a3*T_CS_3*T6c[i]+a3*a3*T_WF_3*T6d[i]-
 			       0.5*a3*Lambda_Q_s*(T7[i]-T7_bar[i]));
 	  }
 	}
       }
       
       if ( (process==QtQtGG||process==tRtRGG) && order>=1) {
             
 	if (process==tRtRGG) {
 	  tRorQ = 2.0;
 	}
 	else {
 	  tRorQ = 1.0;
 	}
 	highC(0,0) += tRorQ*(-1.0*(s*((s+t)*L_t - t*L_u)*y_t*y_t*a3)/(48.*pi*t*u*s));
 	highC(0,1) += tRorQ*((s*L_t*y_t*y_t*a3)/(24.*pi*t*s));
 	highC(0,2) += tRorQ*(-(s*s*y_t*y_t*a3)/((24.*pi*s*t+24.*pi*t*t)*s));
 	highC(0,3) += tRorQ*(-(s*s*y_t*y_t*a3)/((24.*pi*s*t+24.*pi*t*t)*s));
 	highC(1,0) += tRorQ*(-1.0*(s*((s+t)*L_t - t*L_u)*y_t*y_t*a3)/(16.*pi*t*u*s));
 	highC(1,1) += tRorQ*((s*L_t*y_t*y_t*a3)/(8.*pi*t*s));
 	highC(1,2) += tRorQ*(-(s*s*y_t*y_t*a3)/((8.*pi*s*t+8.*pi*t*t)*s));
 	highC(1,3) += tRorQ*(-(s*s*y_t*y_t*a3)/((8.*pi*s*t+8.*pi*t*t)*s));
 	highC(2,0) += tRorQ*((s*((s+t)*L_t + t*L_u)*y_t*y_t*a3)/(16.*pi*t*u*s));
 	highC(2,1) += tRorQ*(((2.*t-2.*t*L_s-s*L_t)*y_t*y_t*a3)/(8.*pi*t*s));
 	highC(2,2) += tRorQ*(-1.0*(s*(s+2.*t)*y_t*y_t*a3)/(8.*pi*t*u*s));
 	highC(2,3) += tRorQ*(-1.0*(s*(s+2.*t)*y_t*y_t*a3)/(8.*pi*t*u*s));
       }
     }
     break;
   default:
     assert(false);
   }  
   return highC;
 }
 
 boost::numeric::ublas::matrix<complex<InvEnergy2> > 
 HighEnergyMatching::Spin0HighMatching(Energy highScale, 
 				      Energy2 s, Energy2 t, Energy2 u,
 				      EWProcess::Process process, 
 				      bool oneLoop, bool ) {
   using Constants::pi;
   unsigned int order = !oneLoop? 0 : 1;
   // (If crossed graphs, swap s and t here)
   Complex L_s = MinusLog(-s/(highScale*highScale));
   Complex L_t = MinusLog(-t/(highScale*highScale));
   Complex L_u = MinusLog(-u/(highScale*highScale));
   Complex L_s2 = L_s*L_s;
   Complex L_t2 = L_t*L_t;
   Complex L_u2 = L_u*L_u;
 
   // Tree-Level:
   complex<InvEnergy2> S1 = 2.0/s;
    
   // Topology T1:
   complex<InvEnergy2> T1b = (-L_s2/(2.0*u)*(7.0*t/s+3.0)+2.0/u*L_t2+L_s*L_t*4.0/u*(t-u)/s+
 			     L_s*2.0/s-4.0/s-pi*pi/(4.0*u)*(11.0+19.0*t/s));
   complex<InvEnergy2> T1b_bar = -1.0*(-L_s2/(2.0*t)*(7.0*u/s+3.0)+2.0/t*L_u2+L_s*L_u*4.0/t*(u-t)/s+
 				      L_s*2.0/s-4.0/s-pi*pi/(4.0*t)*(11.0+19.0*u/s));
    
   // Topologies T2:
   complex<InvEnergy2> T2a = 1.0/s*(-2.0*L_s2+8.0*L_s-16.0+pi*pi/3.0);
   complex<InvEnergy2> T2b = 1.0/s*(0.5*L_s2+2.0*L_s-4.0-pi*pi/12.0);
   
   // Topologies T5:
   complex<InvEnergy2> T5a = 1.0/s*(-2.0*L_s2+6.0*L_s+pi*pi/3.0-16.0);
   complex<InvEnergy2> T5b = 1.0/s*(2.0*L_s-4.0);
   
   // Topologies T6:
   complex<InvEnergy2> T6a = 1.0/s*(-19.0/6.0*L_s+58.0/9.0);
   complex<InvEnergy2> T6b = 1.0/s*(-1.0/6.0*L_s+4.0/9.0);
   complex<InvEnergy2> T6c = 1.0/s*(2.0/3.0*L_s-16.0/9.0);
   complex<InvEnergy2> T6d = 1.0/s*(4.0/3.0*L_s-20.0/9.0);
    
   // Group Theory Factors / SM parameters needed for matrix elements:
   double a1 = ElectroWeakReweighter::coupling()->a1(highScale);
   double a2 = ElectroWeakReweighter::coupling()->a2(highScale);
   double a3 = ElectroWeakReweighter::coupling()->a3(highScale);
   double y_t = ElectroWeakReweighter::coupling()->y_t(highScale);
   double Y_phi = 1.0/2.0;
   double C_F_3 = 4.0/3.0;
   double C_F_2 = 3.0/4.0;
   double C_F_1 = 1.0;
   double n_g = 3.0;
   double n_S = 1.0;
   // Factors TBD after fermion helicity is specified:
   double Y_Q(0.), Lambda_Q(0.), Lambda_phi(0.);
   boost::numeric::ublas::matrix<complex<InvEnergy2> > highC(1,1);
   switch (process) {
     
   case QQPhiPhi: case LLPhiPhi:
     // Finish Group Theory Factors:
     if (process==QQPhiPhi) {
       Y_Q = 1.0/6.0;
       Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
       Lambda_phi = C_F_2*a2+Y_phi*Y_phi*a1;
     }
     else if (process==LLPhiPhi) {
       Y_Q = -1.0/2.0;
       Lambda_Q = C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
       Lambda_phi = C_F_2*a2+Y_phi*Y_phi*a1;
     }
     highC.resize(2,1);
     highC(0,0) = S1*(4.0*pi*a2);
     highC(1,0) = S1*(4.0*pi*a1*Y_Q*Y_phi);
     if (order>=1) {
       highC(0,0) += T1b*(0.5*a2*a2+2.0*a1*a2*Y_Q*Y_phi) + 
 	T1b_bar*(-0.5*a2*a2+2.0*a1*a2*Y_Q*Y_phi) + 
 	T2a*(-a2*a2+Lambda_phi*a2) + T2b*a2*a2 + 
 	T5a*(-a2*a2+Lambda_Q*a2) + T5b*a2*a2 + 
 	T6a*2.0*a2*a2 + T6b*2.0*a2*a2 + 
 	T6c*0.5*a2*a2*n_S + T6d*2.0*a2*a2*n_g;
       highC(1,0) += T1b*(3.0/16.0*a2*a2+a1*a1*Y_Q*Y_Q*Y_phi*Y_phi) + 
 	T1b_bar*(3.0/16.0*a2*a2+a1*a1*Y_Q*Y_Q*Y_phi*Y_phi) + 
 	T2a*(Lambda_phi*a1*Y_Q*Y_phi) + T5a*(Lambda_Q*a1*Y_Q*Y_phi) + 
 	T6c*(2.0*a1*a1*n_S*Y_Q*Y_phi*Y_phi*Y_phi) + 
 	T6d*(10.0/3.0*a1*a1*n_g*Y_Q*Y_phi);
       // Top Quark contributions:
       highC(0,0) += -3.0*y_t*y_t*a2/(4.0*pi)/s*(2.0*L_s-4.0);
       highC(1,0) += -3.0*y_t*y_t*a1/(4.0*pi)*(Y_Q*Y_phi)/s*(2.0*L_s-4.0);
     }
     break;
   case UUPhiPhi: case DDPhiPhi: case EEPhiPhi:
     // Finish Group Theory Factors:
     if (process==UUPhiPhi) {
       Y_Q = 2.0/3.0;
       Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
       Lambda_phi = C_F_2*a2 + Y_phi*Y_phi*a1;
     }
     else if (process==DDPhiPhi) {
       Y_Q = -1.0/3.0;
       Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
       Lambda_phi = C_F_2*a2 + Y_phi*Y_phi*a1;
     }
     else if (process==EEPhiPhi) {
       Y_Q = -1.0;
       Lambda_Q = Y_Q*Y_Q*C_F_1*a1;
       Lambda_phi = C_F_2*a2 + Y_phi*Y_phi*a1;
     }
     
     highC.resize(1,1);
     highC(0,0) = ZERO;
     
     highC(0,0) = S1*(4.0*pi*a1*Y_Q*Y_phi);
     
     if (order>=1) {
       highC(0,0) += T1b*(a1*a1*Y_Q*Y_Q*Y_phi*Y_phi) + 
 	T1b_bar*(a1*a1*Y_Q*Y_Q*Y_phi*Y_phi) + 
 	T2a*(Lambda_phi*a1*Y_Q*Y_phi) + T5a*(Lambda_Q*a1*Y_Q*Y_phi) + 
 	T6c*(2.0*a1*a1*n_S*Y_Q*Y_phi*Y_phi*Y_phi) + 
 	T6d*(10.0/3.0*a1*a1*n_g*Y_Q*Y_phi);
       // Top Quark Contribution:
       highC(0,0) += -3.0*y_t*y_t*a1/(4.0*pi)*(Y_Q*Y_phi)/s*(2.0*L_s-4.0);
     }
     break;
   default:
     assert(false);
   }
   return highC;
 }