diff --git a/MatrixElement/EW/EWCouplings.cc b/MatrixElement/EW/EWCouplings.cc
--- a/MatrixElement/EW/EWCouplings.cc
+++ b/MatrixElement/EW/EWCouplings.cc
@@ -1,614 +1,630 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the EWCouplings class.
 //
 
 #include "EWCouplings.h"
 #include "ThePEG/Interface/ClassDocumentation.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/operation.hpp>
 
 using namespace Herwig;
 
 namespace {
 
 Complex trace(boost::numeric::ublas::matrix<Complex> M) {
   assert(M.size1()==M.size2());
   Complex output(0.);
   for(unsigned int ix=0;ix<M.size1();++ix)
     output += M(ix,ix);
   return output;
 }
 }
 
 
 EWCouplings::EWCouplings(unsigned int loops, unsigned int steps, Energy highScale,
 			 Energy lowScale) 
   : ewScale_(91.1876*GeV), highScale_(highScale), lowScale_(lowScale),
     includeSU3_(true), includeEW_(true), initialized_(false), massChoice_(false),
     mZ_(91.1876*GeV), mW_(80.399*GeV), 
     mT_(173.1*GeV), // 179.08045 (should be this?)
     loops_(loops), highSteps_(steps), lowSteps_(steps)
 {}
 
 void EWCouplings::initialize() {
   using Constants::pi;
   if(initialized_) return;
+  initialized_ = true;
   // set the particle masses
   if(massChoice_) {
     mZ_ = getParticleData(ParticleID::Z0   )->mass();
     mW_ = getParticleData(ParticleID::Wplus)->mass();
     mT_ = getParticleData(ParticleID::t    )->mass();
     mH_ = getParticleData(ParticleID::h0   )->mass();
   }
   // logs of scales
   double logEWScale   = log(ewScale_/GeV);
   double logHighScale = log(highScale_/GeV);
   // step size
   double stepsize = (logHighScale - logEWScale)/double(highSteps_);
   // Initialize parameters at the ewScale
   // 32 parameters, mostly zero due massless quarks
   unsigned int N = 32;
   vector<Complex> y(N,0.), dydx(N,0.), yout(N,0.);
   initializeCouplings(y);
   double x = logEWScale;
   derivatives(x,y,dydx);
   // energy scale + 6 parameters: g1,g2,g3,y_t,lambda,vev
   table_ = boost::numeric::ublas::matrix<double>(highSteps_+1,7);
   table_(0,0) = logEWScale;
   for (unsigned int i=1; i<N; i++) {
     if (i <=3)  table_(0,i   ) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
     if (i >=29) table_(0,i-25) = y[i].real();
   }
   int counter = 1;
    
   // Use 4th order runge-kutta to integrate to highScale
   int steps = highSteps_;
   while (steps > 0) {
-//      _RK4(y,dydx,x,stepsize,yout,_loops);
+    RK4(y,dydx,x,stepsize,yout);
 
     // Advance x and calculate derivatives at new starting point
     for(unsigned int j=0; j<N; j++) {
       y[j] = yout[j];  
     }
     x += stepsize;
     derivatives(x,y,dydx);
      
     table_(counter,0) = x;
     for (unsigned int i=1; i<N; i++) {
       if (i<=3 ) table_(counter,i   ) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
       if (i>=29) table_(counter,i-25) = y[i].real();
     }
      
     steps--;
     counter++;
   }
-   
-//    _LowInit(); // Initialize couplings at mu < 91.1876 GeV
+  // Initialize couplings at mu < 91.1876 GeV
+  initializeLow();
 }
 
 EWCouplings::~EWCouplings() {}
 
 IBPtr EWCouplings::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr EWCouplings::fullclone() const {
   return new_ptr(*this);
 }
 
 void EWCouplings::persistentOutput(PersistentOStream & os) const {
   os << ounit(ewScale_,GeV) <<  ounit(highScale_,GeV) <<  ounit(lowScale_,GeV)
      << includeSU3_ << includeEW_ << ounit(mZ_,GeV) << ounit(mW_,GeV)
      << ounit(mT_,GeV) << ounit(mH_,GeV) << massChoice_ << initialized_;
 }
 
 void EWCouplings::persistentInput(PersistentIStream & is, int) {
   is >> iunit(ewScale_,GeV) >>  iunit(highScale_,GeV) >> iunit(lowScale_,GeV)
      >> includeSU3_ >> includeEW_ >> iunit(mZ_,GeV) >> iunit(mW_,GeV) 
      >> iunit(mT_,GeV) >> iunit(mH_,GeV) >> massChoice_ >> initialized_;
 }
 
 
 //The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<EWCouplings,Interfaced>
 describeHerwigEWCouplings("Herwig::EWCouplings", "HwMEEW.so");
 
 void EWCouplings::Init() {
 
   static ClassDocumentation<EWCouplings> documentation
     ("The EWCouplings class implements");
 
   static Switch<EWCouplings,bool> interfaceMassChoice
     ("MassChoice",
      "Where to get the SM particle masses from",
      &EWCouplings::massChoice_, false, false, false);
   static SwitchOption interfaceMassChoiceLocal
     (interfaceMassChoice,
      "Local",
      "Use local values",
      false);
   static SwitchOption interfaceMassChoiceParticleData
     (interfaceMassChoice,
      "ParticleData",
      "Get the values from the ParticleData object",
      true);
 
 }
 
 
 
 void EWCouplings::initializeLow() {
   using Constants::pi;
   // For scales less than ewScale, the only couplings calculated here are those of
   // alpha_EW and alpha3:
    
   double logEWScale  = log(ewScale_ /GeV);
   double loglowScale = log(lowScale_/GeV);
    
   double stepsize = (loglowScale - logEWScale)/double(lowSteps_);
   int steps = lowSteps_;
    
   // Initialize parameters at the ewScale
   unsigned int N=2; // Total # of parameters = 2
   vector<Complex> y(N), dydx(N), yout(N);
    for (unsigned int i=0; i<N; i++) {
       y[i] = 0.0;
       dydx[i] = 0.0;
       yout[i] = 0.0;
    }
    double x = logEWScale;
    
    // Initialize Couplings at the ewScale, including the One-Loop Threshold Matching:
    double a1 = (3.0/5.0)*table_(0,1);
    double a2 = table_(0,2);
    double a3 = table_(0,3);
    double aEM_ewScale = 1./(1./a1+1./a2-1./(6.*pi)*(1.-21.*log(mW()/ewScale_)+16./3.*log(mTatmZ()/ewScale_)));
    double aS_ewScale = 1./(1./a3+1./(3.*pi)*log(ewScale_/mTatmZ()));
    y[0] = sqrt(4.0*pi*aEM_ewScale);
    y[1] = sqrt(4.0*pi*aS_ewScale);
    
    derivatives(x,y,dydx);
    // energy scale + 2 parameters: gEM,g3
    lowTable_.resize(lowSteps_+1,3);
    lowTable_(0,0) = logEWScale;
    for (unsigned int i=1; i<=N; i++) {
       lowTable_(0,i) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
    }
    int counter = 1;
    
    // Use 4th order runge-kutta to integrate to highScale
    while (steps > 0) {
      RK4(y,dydx,x,stepsize,yout);
      // Advance x and calculate derivatives at new starting point
      for(unsigned int j=0; j<N; j++) {
        y[j] = yout[j];  
      }
      x += stepsize;
      derivatives(x,y,dydx);
 		
      lowTable_(counter,0) = x;
      for (unsigned int i=1; i<=N; i++) {
        lowTable_(counter,i) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
      }
 
       steps--;
       counter++;
    }
 }
 
 void EWCouplings::RK4(vector<Complex> &y, vector<Complex> &dydx,
 		      const double x, const double h, vector<Complex> &yout) {
    
   unsigned int n = y.size();
   std::vector<Complex> dym(n), dyt(n), yt(n);
   double hh = h*0.5;
   double h6 = h/6.0;
   double xh = x + hh;
   const Complex I(0,1.0);
 	
   for(unsigned int i=0; i<n; i++) yt[i] = y[i] + hh*dydx[i];
   derivatives(xh, yt, dyt);
   for(unsigned int i=0; i<n; i++) yt[i] = y[i] + hh*dyt[i];
   derivatives(xh, yt, dym);
   for(unsigned int i=0; i<n; i++) {
     yt[i] = y[i] + h*dym[i];
     dym[i] += dyt[i];
   }
   derivatives(x+h, yt, dyt);
   for(unsigned int i=0; i<n; i++) {
     yout[i] = y[i] + h6*(dydx[i] + dyt[i] + 2.0*dym[i]);
   }
 }
 
 
 void EWCouplings::initializeCouplings(vector<Complex> & y) {
   // \todo make these values parameters so they can be changed
   InvEnergy2 gFermi = 1.16637*pow(10.0,-5)/GeV2;
   Energy vev = 1.0/(sqrt(sqrt(2.0)*gFermi)); // vev = 246.221
   
   y[0] = 0.461531463;     // g1 = Sqrt[5/3] * Sqrt[4*pi*a1] with a1(Mz) = 0.01017054
   y[1] = 0.651547066;     // g2 = Sqrt[4*pi*a2] with a2(Mz) = 0.03378168
   y[2] = 1.215650108;     // g3 = Sqrt[4*pi*as] with as(Mz) = 0.1176
   
   // Note lambda_t = sqrt(2.0)*mt/vev only valid for mt(mt); need mt(mZ) here
   // Top Yukawa lambda from Manohar
   //Complex lambda_t = 1.02858; 
   // Top Yukawa lambda from Sascha          
   double lambda_t = 0.991172;  
   // Quartic coupling lambda (need to multiply by a factor of 2 when accessing the quartic coupling)
   double lambda = (mH_/vev)*(mH_/vev);
   y[29] = lambda_t;
   y[30] = lambda;
   y[31] = vev/GeV;
 }
 
 void EWCouplings::derivatives(const double x, vector<Complex> & y, 
 			      vector<Complex> &dydx) {
   // zero the output
   for (unsigned int i=0; i<dydx.size(); i++) dydx[i]=0.0;
   // low scale
   if (y.size()==2) {
     lowBetaGauge(x,y,dydx);
   }
   // high scale
   else {
     betaGauge(x,y,dydx);
     betaYukawa(x,y,dydx);
     betaHiggs(x,y,dydx);
   }
 }
 
 void EWCouplings::betaGauge(const double x, vector<Complex> &y, vector<Complex> & dydx) {
   using Constants::pi;
   using boost::numeric::ublas::axpy_prod;
   using boost::numeric::ublas::herm;
   const Complex I(0,1.0);
   // Yukawa
   boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
   for(unsigned int ix=0;ix<3;++ix) {
     for(unsigned int iy=0;iy<3;++iy) {
       Yuk_u(ix,iy) = y[21+3*ix+iy];
       Yuk_d(ix,iy) = y[12+3*ix+iy];
       Yuk_e(ix,iy) = y[ 3+3*ix+iy];
     }
   }
   // gauge
   boost::numeric::ublas::vector<Complex> gauge(3);
   for(unsigned int l=0; l<3; l++) gauge[l] = y[l];
   // Evaluate beta functions for gauge couplings
   double Ng = 0.5*numberOfFlavours(x);
   boost::numeric::ublas::vector<Complex>  b(3), g2(3), gc(3), Cu(3), Cd(3), Ce(3);
   boost::numeric::ublas::matrix<Complex> B1(3,3),B2(3,3), B3(3,3), B(3,3);
   b[0] = -4.0/3.0*Ng - 1.0/10.0;
   b[1] = 22.0/3.0 - 4.0/3.0*Ng - 1.0/6.0;
   b[2] = 11.0 - 4.0/3.0*Ng;
   for(unsigned int ix=0;ix<3;++ix) {
     for(unsigned int iy=0;iy<3;++iy) {
       B1(ix,iy) = 0.;
       B2(ix,iy) = 0.;
       B3(ix,iy) = 0.;
     }
   }
   B1(1,1) = 136.0/3.0;
   B1(2,2) = 102.0;
   B2(0,0) = 19.0/15.0;
   B2(0,1) = 1.0/5.0;
   B2(0,2) = 11.0/30.0;
   B2(1,0) = 3.0/5.0;
   B2(1,1) = 49.0/3.0;
   B2(1,2) = 3.0/2.0 ;
   B2(2,0) = 44.0/15.0;
   B2(2,1) = 4.0;
   B2(2,2) = 76.0/3.0;
   B3(0,0) = 9.0/50.0;
   B3(0,1) = 3.0/10.0;
   B3(1,0) = 9.0/10.0;
   B3(1,1) = 13.0/6.0;
   B = B1 - Ng*B2 - B3;
   Cu[0] = 17.0/10.0;
   Cu[1] = 3.0/2.0;
-  Cu[1] = 2.0;
+  Cu[2] = 2.0;
   Cd[0] = 1.0/2.0;
   Cd[1] = 3.0/2.0;
-  Cd[1] = 2.0;
+  Cd[2] = 2.0;
   Ce[0] = 3.0/2.0;
   Ce[1] = 1.0/2.0;
-  Ce[1] = 0.0;
+  Ce[2] = 0.0;
   for (int i=0; i<3; i++) {
     g2(i) = pow(gauge(i),2);
   }
   // gc = trans(g2) * B
   axpy_prod(B, g2, gc, true);
   // compute the answer
   if(loops_ >= 1) {
     for(int l=0; l<3; l++) {
       dydx[l] += -b(l)*pow(gauge(l),3)/(16.0*pow(pi,2));
     }
     if (loops_ >= 2) {
-      boost::numeric::ublas::matrix<Complex> temp;
+      boost::numeric::ublas::matrix<Complex> temp(3,3);
       axpy_prod(herm(Yuk_u),Yuk_u,temp);
       Complex tr1 = trace(temp);
       axpy_prod(herm(Yuk_d),Yuk_d,temp);
       Complex tr2 = trace(temp);
       axpy_prod(herm(Yuk_e),Yuk_e,temp);
       Complex tr3 = trace(temp);
       for(int l=0; l<3; l++) {
 	dydx[l] += -pow(gauge(l),3)/pow(16.0*pow(pi,2),2)*
 	  (gc(l) + Cu(l)*tr1 + Cd(l)*tr2 + Ce(l)*tr3 );
       }
     }
   }
 }
 
 void EWCouplings::lowBetaGauge(const double, vector<Complex> &y,
 			       vector<Complex> &dydx) {
   using Constants::pi;
   const Complex I(0,1.0);
   Complex e = y[0],  g3 = y[1];
   // Evaluate beta functions for gauge couplings
   double Nu = 2.0, Nd = 3.0, Nl = 3.0;
   if(loops_ >=1) {
     dydx[0] += (16.0/9.0*Nu + 4.0/9.0*Nd + 4.0/3.0*Nl)*pow(e,3)/pow(4.0*pi,2);
     dydx[1] += (2.0/3.0*(Nu+Nd)-11.0)*pow(g3,3)/pow(4.0*pi,2);
     // Note this also includes the three-loop contribution for alpha_3
     if (loops_ >= 2) { 
       dydx[0] += (64.0/27.0*Nu+4.0/27.0*Nd+4.0*Nl)*pow(e,5)/pow(4.0*pi,4) +
 	(64.0/9.0*Nu+16.0/9.0*Nd)*pow(e,3)*pow(g3,2)/pow(4.0*pi,4);
       dydx[1] += (38.0/3.0*(Nu+Nd)-102.0)*pow(g3,5)/pow(4.0*pi,4) + 
 	(8.0/9.0*Nu+2.0/9.0*Nd)*pow(g3,3)*pow(e,2)/pow(4.0*pi,4) + 
 	(5033.0/18.0*(Nu+Nd)-325.0/54.0*(Nu+Nd)*(Nu+Nd)-2857.0/2.0)*pow(g3,7)/pow(4.0*pi,6);
     }
   }
 }
 
 void EWCouplings::betaYukawa(const double x, vector< Complex > &y, vector<Complex > &dydx) {
   using Constants::pi;
   const Complex I(0,1.0);
   boost::numeric::ublas::identity_matrix<Complex> Id(3,3);
   // Yukawa
   boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
   for(unsigned int ix=0;ix<3;++ix) {
     for(unsigned int iy=0;iy<3;++iy) {
       Yuk_u(ix,iy) = y[21+3*ix+iy];
       Yuk_d(ix,iy) = y[12+3*ix+iy];
       Yuk_e(ix,iy) = y[ 3+3*ix+iy];
     }
   }
   // gauge
   double Ng = 0.5*numberOfFlavours(x);
   boost::numeric::ublas::vector<Complex> gauge(3);
   for(unsigned int l=0; l<3; l++) gauge[l] = y[l];
   Complex lambda = y[30];
   // traces of yukawa matrices
-  boost::numeric::ublas::matrix<Complex> mTemp,MUU,MDD,MLL,MUU2,MDD2,MLL2,MUUDD,MDDUU;
+  boost::numeric::ublas::matrix<Complex> mTemp(3,3),MUU(3,3),MDD(3,3),MLL(3,3),
+    MUU2(3,3),MDD2(3,3),MLL2(3,3),MUUDD(3,3),MDDUU(3,3);
   axpy_prod(herm(Yuk_u),Yuk_u,MUU);
   Complex trU  = trace( MUU);
   axpy_prod(MUU,MUU,MUU2);
   Complex trUU = trace(MUU2);
   axpy_prod(herm(Yuk_d),Yuk_d,MDD);
   Complex trD  = trace( MUU);
   axpy_prod(MDD,MDD,MDD2);
   Complex trDD = trace(MDD2);
   axpy_prod(MUU,MDD,MUUDD);
   Complex trUD = trace(MUUDD);
   axpy_prod(MDD,MUU,MDDUU);
   axpy_prod(herm(Yuk_e),Yuk_e,MLL);
   Complex trL  = trace( MLL);
   axpy_prod(MLL,MLL,MLL2);
   Complex trLL = trace(MLL2);
   Complex g02 = sqr(gauge[0]);
   Complex g12 = sqr(gauge[1]);
   Complex g22 = sqr(gauge[2]);
   // Evaluate beta functions for yukawa couplings
   boost::numeric::ublas::zero_matrix<Complex> zero3x3(3,3);
   boost::numeric::ublas::matrix<Complex> dYuk_udx(zero3x3), dYuk_ddx(zero3x3), dYuk_edx(zero3x3), beta1_u(zero3x3),
     beta1_d(zero3x3), beta1_e(zero3x3), beta2_u(zero3x3), beta2_d(zero3x3), beta2_e(zero3x3);
   Complex Y2 = 3.0*trU+3.0*trD + trL;
   Complex Y4 = (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trU +
     (1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trD + 3.0/4.0*(g02 + g12)*trL;
   Complex chi4 = 27.0/4.0*trUU + 27.0/4.0*trDD + 9.0/4.0*trLL - 6.0/4.0*trUD;
   if(loops_ >= 1) {
     beta1_u = 3.0/2.0*(MUU - MDD) + (Y2 - 17.0/20.0*g02 - 9.0/4.0*g12 - 8.0*g22)*Id;
     beta1_d = 3.0/2.0*(MDD - MUU) + (Y2 - 1.0/4.0*g02 - 9.0/4.0*g12 - 8.0*g22)*Id;
     beta1_e = 3.0/2.0*(MLL) + (Y2 - 9.0/4.0*g02 - 9.0/4.0*g12)*Id;
     
     axpy_prod(Yuk_u,beta1_u,mTemp);
     dYuk_udx += (1.0/(16.0*pow(pi,2)))*mTemp;
     axpy_prod(Yuk_d,beta1_d,mTemp);
     dYuk_ddx += (1.0/(16.0*pow(pi,2)))*mTemp;
     axpy_prod(Yuk_e,beta1_e,mTemp);
     dYuk_edx += (1.0/(16.0*pow(pi,2)))*mTemp;
-      
+    
     if (loops_ >= 2) {
+      Complex l2=sqr(lambda);
       beta2_u = 3.0/2.0*MUU2 - MUUDD - 1.0/4.0*MDDUU + 11.0/4.0*MDD2 + Y2*(5.0/4.0*MDD - 9.0/4.0*MUU) +
-	(-chi4 + 3.0/2.0*pow(lambda,2))*Id - 2.0*lambda*(3.0*MUU + MDD) + (223.0/80.0*g02 + 135.0/16.0*g12 + 16.0*g22)*(MUU) 
-	- (43.0/80.0*g02 - 9.0/16.0*g12 + 16.0*g22)*(MDD) + 5.0/2.0*Y4*Id + 
-	((9.0/200.0 + 29.0/45.0*Ng)*pow(gauge[0],4) - 9.0/20.0*pow(gauge[0]*gauge[1],2) + 19.0/15.0*pow(gauge[0]*gauge[2],2) - (35.0/4.0 - Ng)*pow(gauge[1],4) + 9.0*pow(gauge[1]*gauge[2],2) - (404.0/3.0 - 80.0/9.0*Ng)*pow(gauge[2],4))*Id;
-      beta2_d = 3.0/2.0*MDD2 - MDDUU - 1.0/4.0*MUUDD + 11.0/4.0*MUU2 + Y2*(5.0/4.0*MUU - 9.0/4.0*MDD) + (-chi4 + 3.0/2.0*pow(lambda,2))*Id - 2.0*lambda*(3.0*MDD + MUU) + (187.0/80.0*g02 + 135.0/16.0*g12 + 16.0*g22)*(MDD) - (79.0/80.0*g02 - 9.0/16.0*g12 + 16.0*g22)*(MUU) + 5.0/2.0*Y4*Id - ((29.0/200.0 + 1.0/45.0*Ng)*pow(gauge[0],4) - 27.0/20.0*pow(gauge[0]*gauge[1],2) + 31.0/15.0*pow(gauge[0]*gauge[2],2) - (35.0/4.0 - Ng)*pow(gauge[1],4) + 9.0*pow(gauge[1]*gauge[2],2) - (404.0/3.0 - 80.0/9.0*Ng)*pow(gauge[2],4))*Id;
-      beta2_e = 3.0/2.0*MLL2 - 9.0/4.0*Y2*MLL + (-chi4 + 3.0/2.0*pow(lambda,2))*Id - 6.0*lambda*(MLL) + (387.0/80.0*g02 + 135.0/15.0*g12)*(MLL) + 5.0/2.0*Y4*Id + ((51.0/200.0 + 11.0/5.0*Ng)*pow(gauge[0],4) + 27.0/20.0*pow(gauge[0]*gauge[1],2) - (35.0/4.0 - Ng)*pow(gauge[1],4))*Id;
+	(-chi4 + 3.0/2.0*l2)*Id - 2.0*lambda*(3.0*MUU + MDD) +
+	(223.0/80.0*g02 + 135.0/16.0*g12 + 16.0*g22)*(MUU) -
+	(43.0/80.0*g02 - 9.0/16.0*g12 + 16.0*g22)*(MDD) + 5.0/2.0*Y4*Id + 
+	((9.0/200.0 + 29.0/45.0*Ng)*pow(gauge[0],4) - 9.0/20.0*pow(gauge[0]*gauge[1],2) +
+	 19.0/15.0*pow(gauge[0]*gauge[2],2) - (35.0/4.0 - Ng)*pow(gauge[1],4) + 9.0*pow(gauge[1]*gauge[2],2) - 
+	 (404.0/3.0 - 80.0/9.0*Ng)*pow(gauge[2],4))*Id;
+      beta2_d = 3.0/2.0*MDD2 - MDDUU - 1.0/4.0*MUUDD + 11.0/4.0*MUU2 + Y2*(5.0/4.0*MUU - 9.0/4.0*MDD) + (-chi4 + 3.0/2.0*l2)*Id - 2.0*lambda*(3.0*MDD + MUU) + (187.0/80.0*g02 + 135.0/16.0*g12 + 16.0*g22)*(MDD) - (79.0/80.0*g02 - 9.0/16.0*g12 + 16.0*g22)*(MUU) + 5.0/2.0*Y4*Id - ((29.0/200.0 + 1.0/45.0*Ng)*pow(gauge[0],4) - 27.0/20.0*pow(gauge[0]*gauge[1],2) + 31.0/15.0*pow(gauge[0]*gauge[2],2) - (35.0/4.0 - Ng)*pow(gauge[1],4) + 9.0*pow(gauge[1]*gauge[2],2) - (404.0/3.0 - 80.0/9.0*Ng)*pow(gauge[2],4))*Id;
+      beta2_e = 3.0/2.0*MLL2 - 9.0/4.0*Y2*MLL + (-chi4 + 3.0/2.0*l2)*Id - 6.0*lambda*(MLL) + (387.0/80.0*g02 + 135.0/15.0*g12)*(MLL) + 5.0/2.0*Y4*Id + ((51.0/200.0 + 11.0/5.0*Ng)*pow(gauge[0],4) + 27.0/20.0*pow(gauge[0]*gauge[1],2) - (35.0/4.0 - Ng)*pow(gauge[1],4))*Id;
          
       axpy_prod(Yuk_u,beta2_u,mTemp);
       dYuk_udx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
       axpy_prod(Yuk_d,beta2_d,mTemp);
       dYuk_ddx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
       axpy_prod(Yuk_e,beta2_e,mTemp);
       dYuk_edx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
     }
   }
   
   boost::numeric::ublas::vector<Complex> temp(27);
 
   for(unsigned int ix=0;ix<3;++ix) {
     for(unsigned int iy=0;iy<3;++iy) {
       dydx[ 3+ix+3*iy] = dYuk_edx(ix,iy);
       dydx[12+ix+3*iy] = dYuk_ddx(ix,iy);
       dydx[21+ix+3*iy] = dYuk_udx(ix,iy);
     }
   }
 }
 
 void EWCouplings::betaHiggs(const double x, vector<Complex> &y,
 			    vector<Complex> &dydx) {
   using Constants::pi;
   const Complex I(0,1.0);
   double Ng = 0.5*numberOfFlavours(x);
   // Yukawa
   boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
   for(unsigned int ix=0;ix<3;++ix) {
     for(unsigned int iy=0;iy<3;++iy) {
       Yuk_u(ix,iy) = y[21+3*ix+iy];
       Yuk_d(ix,iy) = y[12+3*ix+iy];
       Yuk_e(ix,iy) = y[ 3+3*ix+iy];
     }
   }
   // gauge
   boost::numeric::ublas::vector<Complex> gauge(3);
   for(int l=0; l<3; l++) gauge(l) = y[l];
   Complex lambda = y[30];
   complex<Energy> vev = y[31]*GeV;
   // Evaluate beta functions for higgs coupling
   Complex beta1_lambda(0.), beta2_lambda(0.), gamma1_vev(0.), gamma2_vev(0.),
     Y2(0.), H(0.), Y4(0.), chi4(0.);
   // traces of yukawa matrices
-  boost::numeric::ublas::matrix<Complex> temp,temp2,MUU,MDD,MLL;
+  boost::numeric::ublas::matrix<Complex> temp(3,3),temp2(3,3),MUU(3,3),MDD(3,3),MLL(3,3);
   axpy_prod(herm(Yuk_u),Yuk_u,MUU);
   Complex trU  = trace( MUU);
   axpy_prod(MUU,MUU,temp);
   Complex trUU = trace(temp);
   axpy_prod(MUU,temp,temp2);
   Complex trUUU = trace(temp2);
   axpy_prod(herm(Yuk_d),Yuk_d,MDD);
   Complex trD  = trace( MUU);
   axpy_prod(MDD,MDD,temp);
   Complex trDD = trace(temp);
   axpy_prod(MDD,temp,temp2);
   Complex trDDD = trace(temp2);
   axpy_prod(MUU,MDD,temp);
   Complex trUD = trace(temp);
   axpy_prod(herm(Yuk_e),Yuk_e,MLL);
   Complex trL  = trace( MLL);
   axpy_prod(MLL,MLL,temp);
   Complex trLL = trace(temp);
   axpy_prod(MLL,temp,temp2);
   Complex trLLL = trace(temp2);
   axpy_prod(MUU+MDD,MDD,temp);
   axpy_prod(MUU,temp,temp2);
   Complex trUUDD = trace(temp2);
 
   Complex g02 = sqr(gauge[0]);
   Complex g12 = sqr(gauge[1]);
   Complex g22 = sqr(gauge[2]);
+  Complex g04 = sqr(g02);
+  Complex g14 = sqr(g12);
+  double pi2 = sqr(pi);
   Y2 = 3.0*trU+3.0*trD + trL;
   Y4 = (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*(trU) + (1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*(trD) + 3.0/4.0*(g02 + g12)*(trL);
   chi4 = 27.0/4.0*trUU + 27.0/4.0*trDD + 9.0/4.0*trLL - 6.0/4.0*trUD;
   H = 3.0*trUU + 3.0*trDD + trLL;
    
   if(loops_ >= 1) {
-    beta1_lambda = 12.0*pow(lambda,2) - (9.0/5.0*g02 + 9.0*g12)*lambda + 9.0/4.0*(3.0/25.0*pow(gauge[0],4)+2.0/5.0*pow(gauge[0]*gauge[1],2) + pow(gauge[1],4)) + 4.0*Y2*lambda - 4.0*H;
+    Complex l2=sqr(lambda);
+    beta1_lambda = 12.0*l2 - (9.0/5.0*g02 + 9.0*g12)*lambda + 9.0/4.0*(3.0/25.0*g04+2.0/5.0*g02*g12 + g14) + 4.0*Y2*lambda - 4.0*H;
     gamma1_vev = 9.0/4.0*(1.0/5.0*g02+g12)-Y2;
       
-    dydx[30] += 1.0/(16.0*pow(pi,2))*beta1_lambda;
-    dydx[31] += vev/(16.0*pow(pi,2))*gamma1_vev/GeV;
+    dydx[30] += 1.0/(16.0*pi2)*beta1_lambda;
+    dydx[31] += vev/(16.0*pi2)*gamma1_vev/GeV;
     
     if (loops_ >= 2) {
-      beta2_lambda = -78.0*pow(lambda,3) + 18.0*(3.0/5.0*g02 + 3.0*g12)*pow(lambda,2) - ( (313.0/8.0 - 10.0*Ng)*pow(gauge[1],4) - 117.0/20.0*pow(gauge[0]*gauge[1],2) + 9.0/25.0*(229.0/4.0+50.0/9.0*Ng)*pow(gauge[0],4) )*lambda + (497.0/8.0 - 8.0*Ng)*pow(gauge[1],6) - 3.0/5.0*(97.0/24.0 + 8.0/3.0*Ng)*g02*pow(gauge[1],4) - 9.0/25.0*(239.0/24.0 + 40.0/9.0*Ng)*pow(gauge[0],4)*g12 - 27.0/125.0*(59.0/24.0 + 40.0/9.0*Ng)*pow(gauge[0],6) - 64.0*g22*(trUU + trDD) - 8.0/5.0*g02*(2.0*trUU - trDD + 3.0*trLL) - 3.0/2.0*pow(gauge[1],4)*Y4 + 10.0*lambda*( (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trU + (1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trD + 3.0/4.0*(g02 + g12)*trL ) + 3.0/5.0*g02*( (-57.0/10.0*g02 + 21.0*g12 )*trU + (3.0/2.0*g02 + 9.0*g12)*trD + (-15.0/2.0*g02 + 11.0*g12)*trL ) - 24.0*pow(lambda,2)*Y2 - lambda*H + 6.0*lambda*trUD + 20.0*(3.0*trUUU + 3.0*trDDD + trLLL) - 12.0*trUUDD;
-      gamma2_vev = -3.0/2.0*pow(lambda,2) - 5.0/2.0*Y4 + chi4 - 27.0/80.0*pow(gauge[0]*gauge[1],2) - (93.0/800.0 + 1.0/2.0*Ng)*pow(gauge[0],4) + (511.0/32.0 - 5.0/2.0*Ng)*pow(gauge[1],4);
+      beta2_lambda = -78.0*lambda*l2 + 18.0*(3.0/5.0*g02 + 3.0*g12)*l2 - 
+	( (313.0/8.0 - 10.0*Ng)*g14 - 117.0/20.0*g02*g12 + 9.0/25.0*(229.0/4.0+50.0/9.0*Ng)*g04 )*lambda + 
+	(497.0/8.0 - 8.0*Ng)*g14*g12 - 3.0/5.0*(97.0/24.0 + 8.0/3.0*Ng)*g02*g14 - 
+	9.0/25.0*(239.0/24.0 + 40.0/9.0*Ng)*g04*g12 - 27.0/125.0*(59.0/24.0 + 40.0/9.0*Ng)*g04*g02 - 
+	64.0*g22*(trUU + trDD) - 8.0/5.0*g02*(2.0*trUU - trDD + 3.0*trLL) - 3.0/2.0*g14*Y4 + 
+	10.0*lambda*( (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trU + (1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trD + 3.0/4.0*(g02 + g12)*trL ) + 
+	3.0/5.0*g02*( (-57.0/10.0*g02 + 21.0*g12 )*trU + (3.0/2.0*g02 + 9.0*g12)*trD + (-15.0/2.0*g02 + 11.0*g12)*trL ) - 24.0*l2*Y2 - lambda*H + 
+	6.0*lambda*trUD + 20.0*(3.0*trUUU + 3.0*trDDD + trLLL) - 12.0*trUUDD;
+      gamma2_vev = -3.0/2.0*l2 - 5.0/2.0*Y4 + chi4 - 27.0/80.0*g02*g12 - 
+	(93.0/800.0 + 1.0/2.0*Ng)*g04 + (511.0/32.0 - 5.0/2.0*Ng)*g14;
          
-      dydx[30] += 1.0/pow(16.0*pow(pi,2),2)*beta2_lambda;
-      dydx[31] += vev/pow(16.0*pow(pi,2),2)*gamma2_vev/GeV;
+      dydx[30] += 1.0/pow(16.0*pi2,2)*beta2_lambda;
+      dydx[31] += vev/pow(16.0*pi2,2)*gamma2_vev/GeV;
     }
   }
 }
 
 double EWCouplings::interpolate(double t, int paramIndex) {
   double stepsize = table_(1,0)-table_(0,0);
   double tol = 0.001*stepsize;
 	
    double logewScale   = log(ewScale_/GeV);
    double loghighScale = log(highScale_/GeV);
    
    if (t<logewScale-tol || t>loghighScale+tol) {
      cerr << "Stepsize: " << stepsize << std::endl;
      cerr << "tol: " << tol << std::endl;
      cerr << "logewScale: " << logewScale << std::endl;
      cerr << "loghighScale: " << loghighScale << std::endl;
      cerr << "paramIndex: " << paramIndex << std::endl;
      cerr << "t: " << t << std::endl;
      
      cerr << "Couplings::_Interp(double t, int parmIndex) trying to obtain parameter ";
      cerr << "value outside the available range. Returning zero." << std::endl;
      assert(false);
    }
    
    // return value at EW scale
    if (abs(t-logewScale)<tol) return table_(0,paramIndex);
    
    // return value at high scale
    if (abs(t-loghighScale)<tol) {
      return table_(table_.size1()-1,paramIndex);
    }
 	
    unsigned int numSteps = int((t-table_(0,0))/stepsize);
    
    // Linear Interpolation:
    double x1 = table_(numSteps,0);
    double y1 = table_(numSteps,paramIndex);
    double x2 = table_(numSteps+1,0);
    double y2 = table_(numSteps+1,paramIndex);
    return y1+((y2-y1)/(x2-x1))*(t-x1);
 }
 
 double EWCouplings::interpolateLow(double t, int paramIndex) {
   double stepsize = lowTable_(0,0)-lowTable_(1,0);
   double tol = 0.00001*stepsize;
 	
   double logewScale  = log(ewScale_ /GeV);
   double loglowScale = log(lowScale_/GeV);
    
   if (t<loglowScale-tol || t>logewScale+tol) {
     cerr<< "Couplings::_LowInterp(double t, int parmIndex) trying to obtain parameter ";
     cerr << "value outside the available range. Returning zero." << std::endl;
     assert(false);
   }
    
   if (abs(t-logewScale)<tol) {
     return lowTable_(0,paramIndex);
   }
    
   if (std::abs(t-loglowScale)<tol) {
     return lowTable_(lowTable_.size1()-1,paramIndex);
   }
 	
   int numSteps = (int)((lowTable_(0,0)-t)/stepsize);
-   
   // Linear Interpolation:
   double x1 = lowTable_(numSteps,0);
   double y1 = lowTable_(numSteps,paramIndex);
   double x2 = lowTable_(numSteps+1,0);
   double y2 = lowTable_(numSteps+1,paramIndex);
-  
   return y1+((y2-y1)/(x2-x1))*(t-x1);
 }
 
 
 
 
 
 
 
 
 
 
 
 
 
 
diff --git a/MatrixElement/EW/EWCouplings.h b/MatrixElement/EW/EWCouplings.h
--- a/MatrixElement/EW/EWCouplings.h
+++ b/MatrixElement/EW/EWCouplings.h
@@ -1,449 +1,449 @@
 // -*- C++ -*-
 #ifndef Herwig_EWCouplings_H
 #define Herwig_EWCouplings_H
 //
 // This is the declaration of the EWCouplings class.
 //
 
 #include "ThePEG/Interface/Interfaced.h"
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/matrix.hpp>
 #include "EWCouplings.fh"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
 /**
  * Here is the documentation of the EWCouplings class.
  *
  * @see \ref EWCouplingsInterfaces "The interfaces"
  * defined for EWCouplings.
  */
 class EWCouplings: public Interfaced {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * The default constructor.
    */
   EWCouplings(unsigned int loops=2, unsigned int steps=200,
 	      Energy highScale=10.*TeV, Energy lowScale=10.*GeV);
 
   /**
    * The destructor.
    */
   virtual ~EWCouplings();
   //@}
 
   /**
    *  Initialize the running couplings
    */
   void initialize();
 
   /**
    * Number of dynamical quarks at $\log\mu = x$ (in GeV)
    * N.B.Integrate out top quark at Mz, not at Mt.
    */
   unsigned int numberOfFlavours(double x) {
     return x >= std::log(ewScale_/GeV) ? 6 : 5;
   }
    
 public:
    
   /**
    *  Set whether or not to include \f$SU(3)\f$ running
    */
   void SU3(bool includeSU3) {includeSU3_ = includeSU3;}
 
   /**
    *  Whether or not to include \f$SU(3)\f$ running
    */
   bool SU3() { return includeSU3_;}
     
   /**
    *  Set whether or not to include EW running
    */  
   void EW(bool includeEW) {includeEW_ = includeEW;}
 
   /**
    *  Whether or not to include EW running
    */
-  bool isEWon() { return includeEW_;}
+  bool EW() { return includeEW_;}
    
   /**
    *  alpha for the U1 gauge coupling at energy mu (in GeV):
    */
   double a1(Energy mu) {
     if (includeEW_) {
       if (mu>=ewScale_) {
 	return (3.0/5.0)*interpolate(log(mu/GeV),1);
       }
       return interpolateLow(log(mu/GeV),1);
     }
     else
       return 0.0;
   }
 
   /**
    *  alpha for the SU2 gauge coupling at energy mu (in GeV):
    */
   double a2(Energy mu) {
     if (includeEW_) {
       if (mu<ewScale_) {
 	return 0.0;
       }
       else
 	return interpolate(log(mu/GeV),2);
     }
     else
       return 0.0;
   }
 
   /**
    *  alpha for the SU3 gauge coupling at energy mu (in GeV):
    */
   double a3(Energy mu) {
     if(includeSU3_) {
       if (mu>=ewScale_) {
 	return interpolate(log(mu/GeV),3);
       }
       else {
 	return interpolateLow(log(mu/GeV),2);
       }
     }
     else
       return 0.0;
   }
 
   /**
    *  alpha for EM
    */
   double aEM(Energy mu) {
     if (includeEW_) {
       if (mu<=ewScale_) {
 	return interpolateLow(log(mu/GeV),1);
       }
       else {
 	double alpha1=a1(mu);
 	double alpha2=a2(mu);
 	return alpha1*alpha2/(alpha1+alpha2);
       }
     }
     return 0.0;
   }
 
   double aS(Energy mu) {
     if(includeSU3_) {
       if (mu<=ewScale_) {
 	return interpolateLow(log(mu/GeV),2);
       }
       else {
 	return interpolate(log(mu/GeV),3);
       }
     }
     else return 0.0;
   }
 
   /**
    *  Top quark Yukawa coupling
    */
   double y_t(Energy mu) {
     if (includeEW_) {
       if(mu<ewScale_)
 	return 0.0;
       else
 	return interpolate(log(mu/GeV),4);
     }
     else
       return 0.0;
   }
 
   /**
    * Quartic scalar coupling lambda (normalization different, hence factor of 2):
    */
   double lambda(Energy mu) {return 2.0*interpolate(log(mu/GeV),5);}
   
   /**
    *  VEV
    */
   Energy vev(Energy mu) {return interpolate(log(mu/GeV),6)*GeV;}
       
   /**
    * \f$\lambda_t\f$
    */
   double lambda_t(Energy mu) {return y_t(mu);}
       
   /**
    *  Z couplings
    */
   //@{
   double g_Lu(Energy mu) {return 0.5-(2.0/3.0)*Sin2thW(mu);}
   double g_Ld(Energy mu) {return -0.5+(1.0/3.0)*Sin2thW(mu);}
   double g_Le(Energy mu) {return (-1.0/2.0)+Sin2thW(mu);}
   double g_Lnu(Energy ) {return (0.5);}
   double g_Ru(Energy mu) {return (-2.0/3.0)*Sin2thW(mu);}
   double g_Rd(Energy mu) {return (1.0/3.0)*Sin2thW(mu);}
   double g_Re(Energy mu) {return Sin2thW(mu);}
   double g_W(Energy mu) {return Cos2thW(mu);}
   double g_phiPlus(Energy mu) {return 0.5-Sin2thW(mu);}
   //@}
 
   /**
    *  \f\cos^2\theta_W\f$
    */
   double Cos2thW(Energy ) {
     //\todo why this value?
     //return mW()*mW()/(mZ()*mZ());
     return (1.-0.2314);
   }
 
   /**
    *  \f\sin^2\theta_W\f$
    */
   double Sin2thW(Energy mu) {return 1.-Cos2thW(mu);}
 
   double aW(Energy mu) {return a2(mu);}
 
   double aZ(Energy mu) {
     //return a2(mu)/Cos2thW(mu); // Same thing, actually
     return a1(mu)+a2(mu);
   }
 
 public: 
 
   /**
    *  Masses of the Standard Model particles
    */
   //@{
   /**
    *  Z mass
    */
   Energy mZ() const { return mZ_;}
  
   /**
    *  W Mass
    */
   Energy mW() const { return mW_;}
 
   /**
    *  Top quark mass
    */
   Energy mT() const {return mT_;}
 
   Energy mTatmZ() { return mT_;}
 
   /**
    *  Higg boson mass
    */
   Energy mH() const {return mH_;}
   //@}
 
 public:
 
   /** @name Functions used by the persistent I/O system. */
   //@{
   /**
    * Function used to write out object persistently.
    * @param os the persistent output stream written to.
    */
   void persistentOutput(PersistentOStream & os) const;
 
   /**
    * Function used to read in object persistently.
    * @param is the persistent input stream read from.
    * @param version the version number of the object when written.
    */
   void persistentInput(PersistentIStream & is, int version);
   //@}
 
   /**
    * The standard Init function used to initialize the interfaces.
    * Called exactly once for each class by the class description system
    * before the main function starts or
    * when this class is dynamically loaded.
    */
   static void Init();
 
 protected:
 
   /** @name Clone Methods. */
   //@{
   /**
    * Make a simple clone of this object.
    * @return a pointer to the new object.
    */
   virtual IBPtr clone() const;
 
   /** Make a clone of this object, possibly modifying the cloned object
    * to make it sane.
    * @return a pointer to the new object.
    */
   virtual IBPtr fullclone() const;
   //@}
 
 private:
 
   /**
    *  Set the initial value of the couplings
    */
   void initializeCouplings(vector<Complex> & y);
 
   /**
    * Assigns numerical values to beta functions 
    * Takes in a point x = log(mu) and the values of y[i] at x and assigns dydx[i] the 
    * value beta[i](mu).  The function Derivs farms out the plugging in to the three
    * functions BetaGauge, BetaYukawa, and BetaHiggs, which evaluates the beta functions
    * for the gauge couplings, yukawa matrices, and higgs quartic coupling/vev, respectively.
    */
   void derivatives(const double x, vector< Complex > &y,
 		   vector< Complex > & dydx);
 
   /**
    * Beta function for the gauge interactions
    */
   void betaGauge(const double x, vector<Complex> &y, vector<Complex> &dydx);
 
   /**
    * Beta function for the gauge interactions at low scales
    */
   void lowBetaGauge(const double x, vector<Complex> &y, vector<Complex> &dydx);
 
   /**
    * Beta function for the Yukawa interactions
    */
   void betaYukawa(const double x, vector<Complex> &y, vector<Complex> &dydx);
 
   /**
    * Beta function for the Higgs interactions
    */
   void betaHiggs(const double x, vector<Complex> &y, vector<Complex> &dydx);
 
   /**
    *  Update the couplings using 4-th order RK
    * Takes in a vector y[i] of function values at a point x and the values of the 
    * first derivatives dydx[i] ( = dy[i]/dx ) alon with a step size stepsize. The 
    * function then defines assigns the value of y[i](x+stepsize) to the variable yout[i].
    * (Adapted from sample code in Numerical Recipes in C++ Press, Teukolsky, et. al.)
    */
   void RK4(vector<Complex> & y, vector<Complex> &dydx, const double x, 
 	   const double stepsize, vector<Complex> &yout);
 
   /**
    *  Initialize the low energy parameters
    */
   void initializeLow();
 
   /**
    *  Interpolate the table,  t = ln(mu)
    */
   double interpolate(double t, int paramIndex);
 
   /**
    *  Interpolate the tabel, t = ln(mu)
    */
   double interpolateLow(double t, int paramIndex);
 
 private:
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   EWCouplings & operator=(const EWCouplings &);
 
 private:
 
   /**
    *  Electoweak Scale
    */
   Energy ewScale_;
 
   /**
    *  High Scale
    */
   Energy highScale_;
 
   /**
    *  Low Scale
    */
   Energy lowScale_;
 
   /**
    *  Whether or not to include SU3
    */
   bool includeSU3_;
 
   /**
    *  Whether or not to include EW
    */
   bool includeEW_;
 
   /**
    *  Whether or not the running couplings have been initialized
    */
   bool initialized_;
 
   /**
    *  Masses of Standard Model particles
    */
   //@{
   /**
    *  Mass Choice
    */
   bool massChoice_;
 
   /**
    *   Z mass 
    */
   Energy mZ_;
 
   /**
    *   W mass 
    */
   Energy mW_;
 
   /**
    *  Top mass
    */
   Energy mT_;
 
   /**
    *  Higgs boson mass
    */
   Energy mH_;
   //@}
 
   /**
    *  Number of loops
    */
   unsigned int loops_;
 
   /**
    *  Number of steps for Runga-Kutta (High scale)
    */
   unsigned int highSteps_;
 
   /**
    *  Number of steps for Runga-Kutta (Low scale)
    */
   unsigned int lowSteps_;
 
   /**
    *  Matrix to store the parameters
    */
   boost::numeric::ublas::matrix<double> table_;
 
   /**
    *  Matrix to store the low energy parameters.
    *  This will hold only aEM and a3 at mu<ewScale
    */
   boost::numeric::ublas::matrix<double> lowTable_;
 
 };
 
 }
 
 #endif /* Herwig_EWCouplings_H */
diff --git a/MatrixElement/EW/ElectroWeakReweighter.cc b/MatrixElement/EW/ElectroWeakReweighter.cc
--- a/MatrixElement/EW/ElectroWeakReweighter.cc
+++ b/MatrixElement/EW/ElectroWeakReweighter.cc
@@ -1,66 +1,95 @@
 // -*- 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/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"
 
 using namespace Herwig;
 
 ElectroWeakReweighter::ElectroWeakReweighter() {}
 
 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_;
 }
 
 void ElectroWeakReweighter::persistentInput(PersistentIStream & is, int) {
   is >> EWCouplings_;
 }
 
 
 // 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);
 
 }
 
 double ElectroWeakReweighter::weight() const {
   EWCouplings_->initialize();
+  staticEWCouplings_ = EWCouplings_;
+  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";
+  }
+
+
   cerr << "testing got here ???\n";
   cerr <<  subProcess() << "\n";
   cerr << *subProcess() << "\n";
   cerr << subProcess()->outgoing()[0] << *subProcess()->outgoing()[0] << "\n";
   cerr << subProcess()->outgoing()[0]->spinInfo() << "\n";
   cerr << subProcess()->outgoing()[0]->spinInfo()->productionVertex() << "\n";
   assert(false);
+  staticEWCouplings_ = tEWCouplingsPtr();
 }
diff --git a/MatrixElement/EW/ElectroWeakReweighter.h b/MatrixElement/EW/ElectroWeakReweighter.h
--- a/MatrixElement/EW/ElectroWeakReweighter.h
+++ b/MatrixElement/EW/ElectroWeakReweighter.h
@@ -1,108 +1,121 @@
 // -*- C++ -*-
 #ifndef Herwig_ElectroWeakReweighter_H
 #define Herwig_ElectroWeakReweighter_H
 //
 // This is the declaration of the ElectroWeakReweighter class.
 //
 
 #include "ThePEG/MatrixElement/ReweightBase.h"
 #include "EWCouplings.h"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
 /**
  * The ElectroWeakReweighter class.
  *
  * @see \ref ElectroWeakReweighterInterfaces "The interfaces"
  * defined for ElectroWeakReweighter.
  */
 class ElectroWeakReweighter: public ReweightBase {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * The default constructor.
    */
   ElectroWeakReweighter();
 
   /**
    * The destructor.
    */
   virtual ~ElectroWeakReweighter();
   //@}  
 
 public:
 
   /**
    * Return the weight for the kinematical configuation provided by
    * the assigned XComb object (in the LastXCombInfo base class).
    */
   virtual double weight() const;
 
+  /**
+   *
+   */
+  static tEWCouplingsPtr coupling() {
+    assert(staticEWCouplings_);
+    return staticEWCouplings_;
+  }
+
 public:
 
   /** @name Functions used by the persistent I/O system. */
   //@{
   /**
    * Function used to write out object persistently.
    * @param os the persistent output stream written to.
    */
   void persistentOutput(PersistentOStream & os) const;
 
   /**
    * Function used to read in object persistently.
    * @param is the persistent input stream read from.
    * @param version the version number of the object when written.
    */
   void persistentInput(PersistentIStream & is, int version);
   //@}
 
   /**
    * The standard Init function used to initialize the interfaces.
    * Called exactly once for each class by the class description system
    * before the main function starts or
    * when this class is dynamically loaded.
    */
   static void Init();
 
 protected:
 
   /** @name Clone Methods. */
   //@{
   /**
    * Make a simple clone of this object.
    * @return a pointer to the new object.
    */
   virtual IBPtr clone() const;
 
   /** Make a clone of this object, possibly modifying the cloned object
    * to make it sane.
    * @return a pointer to the new object.
    */
   virtual IBPtr fullclone() const;
   //@}
 
 private:
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   ElectroWeakReweighter & operator=(const ElectroWeakReweighter &);
 
 private:
 
   /**
    *  The Electroweak Couplings
    */
   EWCouplingsPtr EWCouplings_;
 
+  /**
+   *  The couplings to allow global access
+   */
+  static tEWCouplingsPtr staticEWCouplings_;
+
 };
 
 }
 
 #endif /* Herwig_ElectroWeakReweighter_H */
diff --git a/MatrixElement/EW/GroupInvariants.cc b/MatrixElement/EW/GroupInvariants.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/GroupInvariants.cc
@@ -0,0 +1,268 @@
+// -*- C++ -*-
+//
+// GroupInvariants.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 "GroupInvariants.h"
+#include "ElectroWeakReweighter.h"
+
+using namespace Herwig;
+using namespace GroupInvariants;
+
+double GroupInvariants::K_Factor(unsigned int i,
+				 unsigned int N, bool high) {
+  using Constants::pi;
+  using Constants::zeta3;
+  // Find K_i for SU(N) (or, for U(1) if N==1)
+  // Relevent for finding the cusp anom. dim.
+  if (N!=1 && N!=2 && N!=3) {
+    std::cerr << "Error. In AnomDim, function K_Factor, N!={1,2,3}" 
+	      << std::endl;
+    assert(false);
+  }
+   
+  double CA = C_A(N);
+  double CF = C_F(N);
+  double nF = n_F(N,high);
+  double nS = n_S(N,high);
+  double TF = T_F(N,high);
+  double tS = t_S(N,high);
+  
+  if (i==1) {
+    return (67.0/9.0-1.0/3.0*pi*pi)*CA - 20.0/9.0*nF*TF - 8.0/9.0*tS*nS;
+  }
+  else if (i==2) {
+    return (245.0/6.0-134.0/27.0*pi*pi+22.0/3.0*zeta3+11.0/45.0*pi*pi*pi*pi)*CA*CA + 
+      (-418.0/27.0+40.0/27.0*pi*pi-56.0/3.0*zeta3)*CA*TF*nF + 
+      (-55.0/3.0+16.0*zeta3)*CF*TF*nF - 16.0/27.0*TF*TF*nF*nF;
+  }
+  else {
+    std::cerr << "Error. In AnomDim, function K_Factor, i!={1,2}" << std::endl;
+    assert(false);
+  }
+}
+
+GaugeContributions GroupInvariants::cuspContributions(Energy mu, int K_ORDER,
+						      bool high) {
+  using ThePEG::Constants::pi;
+  if (K_ORDER > 3) {
+    std::cerr << "Cusp anom dim. requested for K_ORDER>3.\n";
+    assert(false);
+  }
+  // couplings (alphas) for U1, SU2, and SU3
+  double a[3];
+  if (high) {
+    a[0] = ElectroWeakReweighter::coupling()->a1(mu)/(4.0*pi);
+    a[1] = ElectroWeakReweighter::coupling()->a2(mu)/(4.0*pi);
+    a[2] = ElectroWeakReweighter::coupling()->a3(mu)/(4.0*pi);
+  }
+  else {
+    a[0] = ElectroWeakReweighter::coupling()->aEM(mu)/(4.0*pi);
+    a[1] = 0.0;
+    a[2] = ElectroWeakReweighter::coupling()->aS(mu)/(4.0*pi);
+  }
+	
+  // gammaN[0] is really 4.0*C_F, but there is a factor of C_F included in the 
+  // G1-G3 matrices passed from HighRunning.h to the soft integrator. 
+  double gamma1[3];
+  gamma1[0] = 4.0;
+  gamma1[1] = K_Factor(1,1,high)*gamma1[0];
+  gamma1[2] = K_Factor(2,1,high)*gamma1[0];
+  double gamma2[3];
+  gamma2[0] = 4.0;
+  gamma2[1] = K_Factor(1,2,high)*gamma2[0];
+  gamma2[2] = K_Factor(2,2,high)*gamma2[0];
+  double gamma3[3];
+  gamma3[0] = 4.0;
+  gamma3[1] = K_Factor(1,3,high)*gamma3[0];
+  gamma3[2] = K_Factor(2,3,high)*gamma3[0];
+   
+  GaugeContributions result;
+  if (K_ORDER==0) {
+    return result;
+  }
+  // LO bit
+  if (K_ORDER>=1) {
+    result.U1 += a[0]*gamma1[0];
+    result.SU2 += a[1]*gamma2[0];
+    result.SU3 += a[2]*gamma3[0];
+  }
+  // NLO bit
+  if (K_ORDER>=2) {
+    result.U1 += a[0]*a[0]*gamma1[1];
+    result.SU2 += a[1]*a[1]*gamma2[1];
+    result.SU3 += a[2]*a[2]*gamma3[1];
+  }
+  // NNLO bit
+  if (K_ORDER>=3) {
+    result.U1 += a[0]*a[0]*a[0]*gamma1[2];
+    result.SU2 += a[1]*a[1]*a[1]*gamma2[2];
+    result.SU3 += a[2]*a[2]*a[2]*gamma3[2];
+  }
+  return result;
+}
+
+double GroupInvariants::B_Factor(int i, int N, bool fermion,
+				 bool longitudinal) {
+  using Constants::pi;
+  using Constants::zeta3;
+  // Find B_i for SU(N) (or, for U(1) if N==1)
+  // Relevent for finding the collinear non-cusp anom. dim.
+  if (N!=1 && N!=2 && N!=3) {
+    std::cerr << "Error. In AnomDim, function B_Factor, N!={1,2,3}\n";
+    assert(false);
+  }
+   
+  double CA = C_A(N);
+  double CF = C_F(N);
+  double nF = n_F(N,true);
+  double nS = n_S(N,true);
+  double TF = T_F(N,true);
+  double tS = t_S(N,true);
+  
+  if (longitudinal) {
+    if (i==1)      return -8.0*CF;
+    // Two loop non-cusp not known for scalars
+    else if (i==2) return 0.0;
+    else assert(false);
+  }
+  else if (fermion) {
+    if (i==1) return -6.0*CF;
+    else if (i==2) {
+      return (4.0*pi*pi - 48.0*zeta3 - 3.0)*CF*CF + 
+	(52.0*zeta3 - 11.0*pi*pi/3.0 - 961.0/27.0)*CA*CF + 
+	(4.0*pi*pi/3.0 + 260.0/27.0)*CF*TF*nF + 
+	(pi*pi/6.0 + 167.0/54.0)*CF*nS*tS*2.0;
+    }
+    else
+      assert(false);
+  }
+  else {
+    if (i==1) {
+      return -2.0*( 11.0*CA/3.0 - 4.0*TF*nF/3.0 - nS*tS/3.0 );
+    }
+    else if (i==2) {
+      return CA*CA*(11.0*pi*pi/9.0+4.0*zeta3-1384.0/27.0) + 
+	2.0*CA*nF*TF*(256.0/27.0-2.0*pi*pi/9.0) + 8.0*CF*nF*TF;
+    }
+    else assert(false);
+  }
+}
+
+double GroupInvariants::B_Factor_Low(int i, int N, bool fermion,
+				     double boostFactor) {
+  using Constants::pi;
+  using Constants::zeta3;
+  // Find B_i for SU(N) (or, for U(1) if N==1)
+  // Relevent for finding the collinear non-cusp anom. dim.
+  if (N!=1 && N!=2 && N!=3) {
+    std::cerr << "Error. In AnomDim, function B_Factor, N!={1,2,3}\n";
+    assert(false);
+  }
+   
+  double CA = C_A(N);
+  double CF = C_F(N);
+  double nF = n_F(N,false);
+  double nS = n_S(N,false);
+  double TF = T_F(N,false);
+  double tS = t_S(N,false);
+   
+  if (abs(boostFactor)>0.001) {
+    if (i==1)      return -4.0*CF;
+    // Two loop non-cusp not known for bHQET top, W_L, and W_T fields.
+    else if (i==2) return 0.0;
+    else assert(false);
+  }
+  else if (fermion) {
+    if (i==1) return -6.0*CF;
+    else if (i==2) {
+      return (4.0*pi*pi - 48.0*zeta3 - 3.0)*CF*CF + 
+	(52.0*zeta3 - 11.0*pi*pi/3.0 - 961.0/27.0)*CA*CF + 
+	(4.0*pi*pi/3.0 + 260.0/27.0)*CF*TF*nF + 
+	(pi*pi/6.0 + 167.0/54.0)*CF*nS*tS*2.0;
+    }
+    else
+      assert(false);
+  }
+  // Gluon and Photon are the only things left... use Gauge Boson NonCusps:
+  else {
+    if (i==1) return -2.0*( 11.0*CA/3.0 - 4.0*TF*nF/3.0 - nS*tS/3.0 );
+    else if (i==2) {
+      return CA*CA*(11.0*pi*pi/9.0+4.0*zeta3-1384.0/27.0) + 
+	2.0*CA*nF*TF*(256.0/27.0-2.0*pi*pi/9.0) + 8.0*CF*nF*TF;
+    }
+    else 
+      assert(false);
+  }
+}
+
+GaugeContributions GroupInvariants::BContributions(Energy mu,
+						   int B_ORDER,
+						   bool fermion,
+						   bool longitudinal) {
+  using Constants::pi;
+  // NOTE! THIS RETURNS 2*Gamma+sigma, AND SHOULD BE MULTIPLIED BY 1/2 IN
+  // THE COLLINEAR ANOMALOUS DIMENSION INTEGRAND
+  if (B_ORDER>2) {
+    std::cerr << "Non-cusp collinear anom dim. requested for B_ORDER>2.\n";
+    assert(false);
+  }
+	
+  double a[3]; // alpha/(4pi) for U1, SU2, and SU3
+	
+  a[0] = ElectroWeakReweighter::coupling()->a1(mu)/(4.0*pi);
+  a[1] = ElectroWeakReweighter::coupling()->a2(mu)/(4.0*pi);
+  a[2] = ElectroWeakReweighter::coupling()->a3(mu)/(4.0*pi);
+
+  GaugeContributions result;
+  if (B_ORDER==0) {
+    return result;
+  }
+  if (B_ORDER>=1) {
+    result.SU3 += a[2]*B_Factor(1,3,fermion,longitudinal);
+    result.SU2 += a[1]*B_Factor(1,2,fermion,longitudinal);
+    result.U1  += a[0]*B_Factor(1,1,fermion,longitudinal);
+  }
+  if (B_ORDER>=2) {
+    result.SU3 += a[2]*a[2]*B_Factor(2,3,fermion,longitudinal);
+    result.SU2 += a[1]*a[1]*B_Factor(2,2,fermion,longitudinal);
+    result.U1  += a[0]*a[0]*B_Factor(2,1,fermion,longitudinal);
+  }
+  return result;
+}
+
+GaugeContributions GroupInvariants::BContributionsLow(Energy mu,
+						      int B_ORDER,
+						      bool fermion,
+						      double boostFactor) {
+  using Constants::pi;
+  // NOTE! THIS RETURNS 2*Gamma+sigma, AND SHOULD BE MULTIPLIED BY 1/2 IN
+  // THE COLLINEAR ANOMALOUS DIMENSION INTEGRAND
+  if (B_ORDER>2) {
+    std::cerr << "Non-cusp collinear anom dim. requested for B_ORDER>2.\n";
+  }
+  // alpha/(4pi) for U1, SU2, and SU3
+  double a[3]; 
+  a[0] = ElectroWeakReweighter::coupling()->aEM(mu)/(4.0*pi);
+  a[1] = 0.0;
+  a[2] = ElectroWeakReweighter::coupling()->aS(mu)/(4.0*pi);
+  GaugeContributions result;
+  if (B_ORDER==0) {
+    return result;
+  }
+  if (B_ORDER>=1) {
+    result.SU3 += a[2]*B_Factor_Low(1,3,fermion,boostFactor);
+    result.SU2 += a[1]*B_Factor_Low(1,2,fermion,boostFactor);
+    result.U1  += a[0]*B_Factor_Low(1,1,fermion,boostFactor);
+  }
+  if (B_ORDER>=2) {
+    result.SU3 += a[2]*a[2]*B_Factor_Low(2,3,fermion,boostFactor);
+    result.SU2 += a[1]*a[1]*B_Factor_Low(2,2,fermion,boostFactor);
+    result.U1  += a[0]*a[0]*B_Factor_Low(2,1,fermion,boostFactor);
+  }
+  return result;
+}
diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
--- a/MatrixElement/EW/GroupInvariants.h
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -1,195 +1,234 @@
 // -*- C++ -*-
 //
 // GroupInvariants.h 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.
 //
 //
 //
 #ifndef HERWIG_GroupInvariants_H
 #define HERWIG_GroupInvariants_H
+#include "ThePEG/Config/ThePEG.h"
+#include "ThePEG/Config/Unitsystem.h"
+#include <cassert>
 
 namespace Herwig {
-// //////////////////////////////////////////////////////////////////////////
-// // Function header file for Group Theory Invariants                     //
-// // Brian Shotwell (bshotwell@ucsd.edu)                                  //
-// //                                                                      //
-// // This file is a collection of group theory invariants for SU(N)/U(1). //
-// // It includes separate namespaces for high and low scales, and for     //
-// // whether the Fermions are to be treated as Dirac or Weyl.             //
-// //                                                                      //
-// // These factors include appropriate "generalizations" to U(1), and     //
-// // they include a factor of 5/3 (high scale) or 20/3 (low scale) for    //
-// // T_F for U(1).                                                        //
-// //////////////////////////////////////////////////////////////////////////
+using namespace ThePEG;
 
 namespace GroupInvariants {
 
   /**
+   *  Simple struct for storing the different gauge contributions
+   */
+  struct GaugeContributions {
+
+    /**
+     *  Default Constructor
+     */
+    GaugeContributions(double inSU3=0.,
+		       double inSU2=0., double inU1=0.) 
+      : SU3(inSU3),SU2(inSU2),U1(inU1)
+    {}
+    /**
+     * \f$SU(3)\f$
+     */
+    double SU3;
+
+    /**
+     *  \f$SU(2)\f$
+     */
+    double SU2;
+
+    /**
+     * \f$U(1)\f$
+     */
+    double U1;
+  };
+
+
+  /**
    *  The \f$SU(N)\f$ \f$C_A\f$ 
    */
   double C_A(unsigned int N) {
     return N !=1 ? double(N) : 0.;
   }
 
   /**
    *  The \f$SU(N)\f$ \f$C_F\f$ 
    */
   double C_F(unsigned int N) {
-    return N !=1 0.5*(double(N*N))-1.)/double(N) ? : 1.;
+    return N !=1 ? 0.5*(double(N*N)-1.)/double(N) : 1.;
   }
 
   /*
    *  The \f$SU(N)\f$ \f$C_d\f$
    */
   double C_d(unsigned int N) {
     return (double(N*N)-4.)/double(N);
   }
 
   /**
    *  The \f$SU(N)\f$\f$C_1\f$ 
    */
   double C_1(unsigned int N) {
     double N2(N*N);
     return 0.25*(N2-1.0)/N2;
   }
+
+  /**
+   *  \f$T_F\f$
+   */
+  double T_F(unsigned int N, bool high) {
+    if(high) {
+      return N !=1 ? 0.5 : 5.0/3.0;
+    }
+    else {
+      return N !=1 ? 0.5 : 20.0/3.0;
+    }
+  }
+
+  /** 
+   *   \f$t_S\f$
+   */
+  double t_S(unsigned int, bool ) {
+    return 0.5;
+  }
+
+  /**
+   * / Number of complex scalars in the fundamental rep. of SU(N)/U(1)
+   */
+  double n_S(unsigned int N, bool high) {
+    if(high) {
+      if(N==2 || N==1) return 1.0;
+      else if(N==3)    return 0.0;
+      else assert(false);
+    }
+    else {
+      if(N>=1&&N<=3) return 0.;
+      else assert(false);
+    }
+  }
+
+  /**
+   * Number of Dirac Fermions in the fund. rep. of SU(N) (or U(1) for N==1)
+   */
+  double n_F(unsigned int N, bool high) {
+    if(high) {
+      if(N==1) return 3.0;
+      else if(N==2 || N==3) return 6.0;
+      else assert(false);
+    }
+    else {
+      if(N==1) return 1.0;
+      else if(N==2) return 0.0;
+      else if(N==3) return 5.0;
+      else assert(false);
+    }
+  } 
+  
+  /**
+   * Find K_i for gauge group N. high=false for running at mu<mZ
+   */
+  double K_Factor(unsigned int i,unsigned int N, bool high);
+
+  /**
+   *   Find B_i for gauge group N, high energy
+   */
+  double B_Factor(int i, int N, bool fermion, bool longitudinal);
+
+  /**
+   *   Find B_i for gauge group N, low  energy
+   */
+  double B_Factor_Low(int i, int N, bool fermion, double boostFactor);
+
+  /**
+   * Contributions to the Cusps
+   */
+  GaugeContributions cuspContributions(Energy mu, int K_ORDER, bool high);
+
+  /**
+   * Contributions to B, high energy
+   */
+  GaugeContributions BContributions(Energy mu, int B_ORDER,
+				    bool fermion, bool longitudinal);
+
+  /**
+   * Contributions to B, low  energy
+   */
+  GaugeContributions BContributionsLow(Energy mu, int B_ORDER,
+				       bool fermion, double boostFactor);
 }
 
 
 
 // namespace DiracHigh {
 
 
-//     double n_S(int N) { // Number of complex scalars in the fundamental rep. of SU(N)/U(1)
-//         if(N==2 || N==1) return 1.0;
-//         if(N==3) return 0.0;
-//         std::cout << "Error! SU(N), N != (1, 2 or 3) used for n_S in ";
-//         std::cout << "GroupInvariants.h but not defined." << std::endl;
-//         return 0.0;
-//     }
-
-//     double n_F(int N) { // Number of Dirac Fermions in the fund. rep. of SU(N) (or U(1) for N==1)
-// 		if(N==1) return 3.0;
-//         //if(N==1) return 6.0;
-//         if(N==2) return 6.0;
-//         if(N==3) return 6.0;
-//         std::cout << "Error! SU(N), N != (1, 2 or 3) used for n_F in ";
-//         std::cout << "GroupInvariants.h but not defined." << std::endl;
-//         return 0.0;
-//     }
-
 //     double n_g() { // Number of fermion generations (only used in gauge boson HighCMatching)
 //         return 3.0;
 //     }
 
-//     double T_F(int N) { 
-//         if (N==1) {
-//             return 5.0/3.0;
-// 			//return 1.0;
-//         }
-        
-//         return 0.5; // I believe T(N) is such that tr(t^a t^b) = T delta(ab)
-//     }
-
-//     double t_S(int N) { return 0.5+0.0*N; } // Analog of T_F but for scalars.
-//                                             // 0.0*N included to stop receiving a stupid warning.
 
 // }
 
 
 // namespace WeylHigh {
 	
 // 	double n_S(int N) { // Number of complex scalars in the fundamental rep. of SU(N)/U(1)
 // 		if(N==2 || N==1) return 1.0;
 // 		if(N==3) return 0.0;
 // 		std::cout << "Error! SU(N), N != (1, 2 or 3) used for n_S in ";
 //         std::cout << "GroupInvariants.h but not defined." << std::endl;
 // 		return 0.0;
 // 	}
 // 	double n_F(int N) { // Number of Weyl Fermions in the fundamental rep. of SU(N)
 // 		if(N==2) return 12.0;
 // 		if(N==3) return 12.0;
 // 		std::cout << "Error! SU(N), N != (2 or 3) used for n_F in ";
 //         std::cout << "GroupInvariants.h but not defined." << std::endl;
 // 		return 0.0;
 // 	}
 //     double n_g() { // Number of fermion generations (only used in gauge boson HighCMatching)
 //         return 3.0;
 //     }
 // 	double T_F(int N) { return 0.5+0.0*N; } // I believe T(N) is such that tr(t^a t^b) = T delta(ab)
 // 	// 0.0*N included to stop receiving a stupid warning.
 	
 // 	double t_S(int N) { return 0.5+0.0*N; } // Analog of T_F but for scalars.
 // 	// 0.0*N included to stop receiving a stupid warning.
 	
 // }
 
-
-// namespace DiracLow {
-    
-    
-//     double n_S(int N) { // Number of complex scalars in the fundamental rep. of SU(3) or U(1) 
-//         if(N==1) return 0.0;
-//         if(N==2) return 0.0;
-//         if(N==3) return 0.0;
-//         std::cout << "Error! SU(N), N != (1, 2 or 3) used for n_S in ";
-//         std::cout << "GroupInvariants.h but not defined." << std::endl;
-//         return 0.0;
-//     }
-    
-//     double n_F(int N) { // Number of Dirac Fermions in the fund. rep. of SU(3) or U(1)
-// 		if(N==3) return 5.0;
-//         if(N==2) return 0.0;
-// 		if(N==1) return 1.0;
-//         std::cout << "Error! SU(N), N != (1, 2 or 3) used for n_F in ";
-//         std::cout << "GroupInvariants.h but not defined." << std::endl;
-//         return 0.0;
-//     }
-    
-//     double T_F(int N) {
-//         if (N==1) {
-//             return 20.0/3.0;
-//         }
-        
-//         return 0.5+0.0*N; // 0.0*N included to stop receiving a stupid warning.
-//     }
-    
-//     double t_S(int N) { return 0.5+0.0*N; } // Analog of T_F but for scalars.
-//     // 0.0*N included to stop receiving a stupid warning.
-    
-// }
-
-
 // namespace WeylLow {
 	
 // 	double n_S(int N) { // Number of complex scalars in the fundamental rep. of SU(N)
 //         if(N==1) return 0.0;
 //         if(N==3) return 0.0;
 //         std::cout << "Error! SU(N), N != (1 or 3) used for n_S in ";
 //         std::cout << "GroupInvariants.h but not defined." << std::endl;
 //         return 0.0;
 // 	}
 // 	double n_F(int N) { // Number of Weyl Fermions in the fundamental rep. of SU(N)
 // 		if(N==3) return 10.0;
 // 		if(N==1) return 2.0;
 //         std::cout << "Error! SU(N), N != (1 or 3) used for n_F in ";
 //         std::cout << "GroupInvariants.h but not defined." << std::endl;
 //         return 0.0;
 // 	}
 // 	double T_F(int N) { return 0.5+0.0*N; } // I believe T(N) is such that tr(t^a t^b) = T delta(ab)
 // 	// 0.0*N included to stop receiving a stupid warning.
 	
 // 	double t_S(int N) { return 0.5+0.0*N; } // Analog of T_F but for scalars.
 // 	// 0.0*N included to stop receiving a stupid warning.
 	
 // }
 
 
 // #endif // GROUP_INVARIANTS_H
 
 
 
 }
 
 #endif // HERWIG_GroupInvariants_H 
diff --git a/MatrixElement/EW/Makefile.am b/MatrixElement/EW/Makefile.am
--- a/MatrixElement/EW/Makefile.am
+++ b/MatrixElement/EW/Makefile.am
@@ -1,5 +1,5 @@
 pkglib_LTLIBRARIES = HwMEEW.la
-HwMEEW_la_SOURCES = GroupInvariants.h \
-ElectroWeakReweigter.h ElectroWeakReweighter.cc\
+HwMEEW_la_SOURCES = GroupInvariants.h GroupInvariants.cc \
+ElectroWeakReweigter.h ElectroWeakReweighter.cc \
 EWCouplings.h EWCouplings.fh EWCouplings.cc
 HwMEEW_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0