diff --git a/MatrixElement/EW/EWCouplings.cc b/MatrixElement/EW/EWCouplings.cc
--- a/MatrixElement/EW/EWCouplings.cc
+++ b/MatrixElement/EW/EWCouplings.cc
@@ -1,629 +1,734 @@
 // -*- 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/Parameter.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,
+EWCouplings::EWCouplings(unsigned int loops, unsigned int steps,
 			 Energy lowScale) 
-  : ewScale_(91.1876*GeV), highScale_(highScale), lowScale_(lowScale),
+  : ewScale_(91.1876*GeV), highScale_(14.*TeV), 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?)
     mH_(125.0*GeV),
-    loops_(loops), highSteps_(steps), lowSteps_(steps)
+    loops_(loops), highSteps_(steps), lowSteps_(steps),
+    a1MZ_(0.01017054),a2MZ_(0.03378168),aSMZ_(0.1176),lambdat_(0.991172)
 {}
 
 void EWCouplings::initialize() {
   using Constants::pi;
   if(initialized_) return;
   initialized_ = true;
+  highScale_=generator()->maximumCMEnergy();
   // 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);
 
     // 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++;
   }
   // 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)
+  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_;
+     << ounit(mT_,GeV) << ounit(mH_,GeV) << massChoice_ << initialized_
+     << loops_ << highSteps_ << lowSteps_
+     << a1MZ_ << a2MZ_ << aSMZ_ << lambdat_;
+  os << table_.size1() << table_.size2();
+  for(unsigned int ix=0;ix<table_.size1();++ix)
+    for(unsigned int iy=0;iy<table_.size2();++iy)
+      os << table_(ix,iy);
+  os << lowTable_.size1() << lowTable_.size2();
+  for(unsigned int ix=0;ix<lowTable_.size1();++ix)
+    for(unsigned int iy=0;iy<lowTable_.size2();++iy)
+      os << lowTable_(ix,iy);
 }
 
 void EWCouplings::persistentInput(PersistentIStream & is, int) {
-  is >> iunit(ewScale_,GeV) >>  iunit(highScale_,GeV) >> iunit(lowScale_,GeV)
+  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_;
+     >> iunit(mT_,GeV) >> iunit(mH_,GeV) >> massChoice_ >> initialized_
+     >> loops_ >> highSteps_ >> lowSteps_
+     >> a1MZ_ >> a2MZ_ >> aSMZ_ >> lambdat_;
+  unsigned int size1,size2;
+  is >> size1 >> size2;
+  table_.resize(size1,size2);
+  for(unsigned int ix=0;ix<table_.size1();++ix)
+    for(unsigned int iy=0;iy<table_.size2();++iy)
+      is >> table_(ix,iy);
+  is >> size1 >> size2;
+  lowTable_.resize(size1,size2);
+  for(unsigned int ix=0;ix<lowTable_.size1();++ix)
+    for(unsigned int iy=0;iy<lowTable_.size2();++iy)
+      is >> lowTable_(ix,iy);
 }
 
 
 //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);
 
+  static Parameter<EWCouplings,Energy> interfaceEWScale
+    ("EWScale",
+     "The electroweak scale for matching between high and low energy running",
+     &EWCouplings::ewScale_, GeV, 91.1876*GeV, 10.0*GeV, 10000.0*GeV,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,Energy> interfaceLowScale
+    ("LowScale",
+     "The low energy scale at which to stop the running",
+     &EWCouplings::lowScale_, GeV, 10.*GeV, 1.0*GeV, 100.0*GeV,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,Energy> interfacemZ
+    ("mZ",
+     "The mass of the Z boson",
+     &EWCouplings::mZ_, 91.1876*GeV, 90.*GeV, 92.*GeV,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,Energy> interfacemW
+    ("mW",
+     "The mass of the W boson",
+     &EWCouplings::mW_, 80.399*GeV, 75.*GeV, 85.*GeV,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,Energy> interfacemT
+    ("mT",
+     "The mass of the top quark",
+     &EWCouplings::mT_, 173.1*GeV, 100.*GeV, 1000.*GeV,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,Energy> interfacemH
+    ("mH",
+     "The mass of the Higgs boson",
+     &EWCouplings::mH_, 125.*GeV, 100.*GeV, 1000.*GeV,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,unsigned int> interfaceLoops
+    ("Loops",
+     "The number of loops",
+     &EWCouplings::loops_, 2, 1, 3,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,unsigned int> interfaceHighSteps
+    ("HighSteps",
+     "The number of steps for the Runga-Kutta and interpolation of the couplings above mZ",
+     &EWCouplings::highSteps_, 200, 10, 1000000,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,unsigned int> interfaceLowSteps
+    ("LowSteps",
+     "The number of steps for the Runga-Kutta and interpolation of the couplings below mZ",
+     &EWCouplings::lowSteps_, 200, 10, 1000000,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,double> interfacea1MZ
+    ("Alpha1MZ",
+     "The value of a1(MZ)",
+     &EWCouplings::a1MZ_, 0.01017054, 0.0, 1.,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,double> interfacea2MZ
+    ("Alpha2MZ",
+     "The value of a2(MZ)",
+     &EWCouplings::a2MZ_, 0.03378168, 0.0, 1.,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,double> interfaceasMZ
+    ("AlphasMZ",
+     "The value of as(MZ)",
+     &EWCouplings::aSMZ_, 0.1176, 0.0, 1.,
+     false, false, Interface::limited);
+
+  static Parameter<EWCouplings,double> interfaceLambdaT
+    ("LambdaT",
+     "The top quark Yukawa at the matching scale",
+     &EWCouplings::lambdat_, 0.991172, 0.5, 2.0,
+     false, false, Interface::limited);
+
 }
 
 
 
 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) {
+  using Constants::pi;
   // \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
+  y[0] = sqrt(5./3.*4.*pi*a1MZ_); // g1 = Sqrt[5/3] * Sqrt[4*pi*a1]
+  y[1] = sqrt(4.*pi*a2MZ_);     // g2 = Sqrt[4*pi*a2]
+  y[2] = sqrt(4.*pi*aSMZ_);     // g3 = Sqrt[4*pi*as]
   // 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;  
+  double lambda_t =lambdat_;
   // 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[2] = 2.0;
   Cd[0] = 1.0/2.0;
   Cd[1] = 3.0/2.0;
   Cd[2] = 2.0;
   Ce[0] = 3.0/2.0;
   Ce[1] = 1.0/2.0;
   Ce[2] = 0.0;
   for (int i=0; i<3; i++) {
     g2(i) = pow(gauge(i),2);
   }
   // gc = trans(g2) * B
   axpy_prod(g2, B, 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(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(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( MDD);
   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*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(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( MDD);
   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) {
     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*pi2)*beta1_lambda;
     dydx[31] += vev/(16.0*pi2)*gamma1_vev/GeV;
     
     if (loops_ >= 2) {
       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*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,473 @@
 // -*- 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=14.*TeV, Energy lowScale=10.*GeV);
+  EWCouplings(unsigned int loops=2, unsigned int steps=200, 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 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_;
 
+  /**
+   *  Input values of the couplings
+   */
+  //@{
+  /**
+   *   \f%\alpha_1(M_Z)\f$
+   */
+  double a1MZ_;
+
+  /**
+   *   \f%\alpha_2(M_Z)\f$
+   */
+  double a2MZ_;
+
+  /**
+   *   \f%\alpha_S(M_Z)\f$
+   */
+  double aSMZ_;
+
+  /**
+   *  \f$\lambda_t\f$
+   */
+  double lambdat_;
+  //@}
+
 };
 
 }
 
 #endif /* Herwig_EWCouplings_H */