diff --git a/MatrixElement/EW/EWCouplings.cc b/MatrixElement/EW/EWCouplings.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/EWCouplings.cc
@@ -0,0 +1,614 @@
+// -*- 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;
+  // 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);
+
+    // 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
+}
+
+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;
+  Cd[0] = 1.0/2.0;
+  Cd[1] = 3.0/2.0;
+  Cd[1] = 2.0;
+  Ce[0] = 3.0/2.0;
+  Ce[1] = 1.0/2.0;
+  Ce[1] = 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;
+      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;
+  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) {
+      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;
+         
+      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;
+  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]);
+  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;
+    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;
+    
+    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);
+         
+      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;
+    }
+  }
+}
+
+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.fh b/MatrixElement/EW/EWCouplings.fh
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/EWCouplings.fh
@@ -0,0 +1,18 @@
+// -*- C++ -*-
+//
+// This is the forward declaration of the EWCouplings class.
+//
+#ifndef Herwig_EWCouplings_FH
+#define Herwig_EWCouplings_FH
+
+#include "ThePEG/Config/ThePEG.h"
+
+namespace Herwig {
+
+class EWCouplings;
+
+ThePEG_DECLARE_POINTERS(Herwig::EWCouplings,EWCouplingsPtr);
+
+}
+
+#endif
diff --git a/MatrixElement/EW/EWCouplings.h b/MatrixElement/EW/EWCouplings.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/EWCouplings.h
@@ -0,0 +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_;}
+   
+  /**
+   *  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
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/ElectroWeakReweighter.cc
@@ -0,0 +1,66 @@
+// -*- 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();
+  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);
+}
diff --git a/MatrixElement/EW/ElectroWeakReweighter.h b/MatrixElement/EW/ElectroWeakReweighter.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/ElectroWeakReweighter.h
@@ -0,0 +1,108 @@
+// -*- 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;
+
+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_;
+
+};
+
+}
+
+#endif /* Herwig_ElectroWeakReweighter_H */
diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -0,0 +1,195 @@
+// -*- 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
+
+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).                                                        //
+// //////////////////////////////////////////////////////////////////////////
+
+namespace GroupInvariants {
+
+  /**
+   *  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.;
+  }
+
+  /*
+   *  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;
+  }
+}
+
+
+
+// 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/Makefile.am b/MatrixElement/EW/Makefile.am
copy from MatrixElement/Makefile.am
copy to MatrixElement/EW/Makefile.am
--- a/MatrixElement/Makefile.am
+++ b/MatrixElement/EW/Makefile.am
@@ -1,11 +1,5 @@
-SUBDIRS = General Lepton Hadron DIS Powheg Gamma Matchbox Reweighters
-
-noinst_LTLIBRARIES = libHwME.la 
-
-libHwME_la_SOURCES =  \
-HwMEBase.h HwMEBase.fh HwMEBase.cc \
-MEfftoVH.h MEfftoVH.cc \
-MEfftoffH.h MEfftoffH.cc \
-HardVertex.fh HardVertex.h HardVertex.cc \
-ProductionMatrixElement.h ProductionMatrixElement.cc \
-DrellYanBase.h DrellYanBase.cc
+pkglib_LTLIBRARIES = HwMEEW.la
+HwMEEW_la_SOURCES = GroupInvariants.h \
+ElectroWeakReweigter.h ElectroWeakReweighter.cc\
+EWCouplings.h EWCouplings.fh EWCouplings.cc
+HwMEEW_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
diff --git a/MatrixElement/Makefile.am b/MatrixElement/Makefile.am
--- a/MatrixElement/Makefile.am
+++ b/MatrixElement/Makefile.am
@@ -1,11 +1,11 @@
-SUBDIRS = General Lepton Hadron DIS Powheg Gamma Matchbox Reweighters
+SUBDIRS = General Lepton Hadron DIS Powheg Gamma Matchbox Reweighters EW
 
 noinst_LTLIBRARIES = libHwME.la 
 
 libHwME_la_SOURCES =  \
 HwMEBase.h HwMEBase.fh HwMEBase.cc \
 MEfftoVH.h MEfftoVH.cc \
 MEfftoffH.h MEfftoffH.cc \
 HardVertex.fh HardVertex.h HardVertex.cc \
 ProductionMatrixElement.h ProductionMatrixElement.cc \
 DrellYanBase.h DrellYanBase.cc
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,244 +1,245 @@
 dnl Process this file with autoconf to produce a configure script.
 
 AC_PREREQ([2.63])
 AC_INIT([Herwig],[trunk],[herwig@projects.hepforge.org],[Herwig])
 AC_CONFIG_SRCDIR([Utilities/HerwigStrategy.cc])
 AC_CONFIG_AUX_DIR([Config])
 AC_CONFIG_MACRO_DIR([m4])
 AC_CONFIG_HEADERS([Config/config.h])
 dnl AC_PRESERVE_HELP_ORDER
 AC_CANONICAL_HOST
 
 dnl === disable debug symbols by default =====
 if test "x$CXXFLAGS" = "x"; then
    CXXFLAGS=-O3
 fi
 if test "x$CFLAGS" = "x"; then
    CFLAGS=-O3
 fi
 
 AC_LANG([C++])
 
 AM_INIT_AUTOMAKE([1.11 subdir-objects gnu dist-bzip2 no-dist-gzip -Wall])
 m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
 m4_ifdef([AM_PROG_AR], [AM_PROG_AR])
 
 dnl Checks for C++ compiler. Handle C++11 flags.
 AC_PROG_CXX
 AX_CXX_COMPILE_STDCXX([11],[noext],[mandatory])
 
 dnl check for POSIX
 AC_CHECK_HEADER([unistd.h],[],
       [AC_MSG_ERROR([Herwig needs "unistd.h". Non-POSIX systems are not supported.])])
 
 dnl Checks for programs.
 AC_PROG_INSTALL
 AC_PROG_MAKE_SET
 AC_PROG_LN_S
 
 dnl modified search order
 AC_PROG_FC([gfortran g95 g77]) 
 dnl xlf95 f95 fort ifort ifc efc pgf95 lf95 ftn xlf90 f90 pgf90 pghpf epcf90 xlf f77 frt pgf77 cf77 fort77 fl32 af77])
 AC_LANG_PUSH([Fortran])
 AC_MSG_CHECKING([if the Fortran compiler ($FC) works])
 AC_COMPILE_IFELSE(
    	AC_LANG_PROGRAM([],[      print *[,]"Hello"]),
 	[AC_MSG_RESULT([yes])],
 	[AC_MSG_RESULT([no])
 	 AC_MSG_ERROR([A Fortran compiler is required to build Herwig.])
 	]
 )
 AC_LANG_POP([Fortran])
 
 
 LT_PREREQ([2.2.6])
 LT_INIT([disable-static dlopen pic-only])
 
 dnl ####################################
 dnl ####################################
 
 dnl for Doc/fixinterfaces.pl
 AC_PATH_PROG(PERL, perl)
 
 dnl for Models/Feynrules
 AM_PATH_PYTHON([2.6],, [:])
 AM_CONDITIONAL([HAVE_PYTHON], [test "x$PYTHON" != "x:"])
 
 HERWIG_CHECK_GSL
 
 HERWIG_CHECK_THEPEG
 
 BOOST_REQUIRE([1.41])
 BOOST_FIND_HEADER([boost/array.hpp])
 BOOST_FIND_HEADER([boost/numeric/ublas/io.hpp])
 BOOST_FIND_HEADER([boost/numeric/ublas/matrix.hpp])
 BOOST_FIND_HEADER([boost/numeric/ublas/matrix_proxy.hpp])
 BOOST_FIND_HEADER([boost/numeric/ublas/matrix_sparse.hpp])
 BOOST_FIND_HEADER([boost/numeric/ublas/symmetric.hpp])
 BOOST_FIND_HEADER([boost/numeric/ublas/vector.hpp])
 BOOST_FIND_HEADER([boost/operators.hpp])
 BOOST_FIND_HEADER([boost/progress.hpp])
 BOOST_FIND_HEADER([boost/scoped_array.hpp])
 BOOST_FIND_HEADER([boost/scoped_ptr.hpp])
 BOOST_FIND_HEADER([boost/utility.hpp])
 BOOST_FILESYSTEM([mt])
 BOOST_TEST()
 
 HERWIG_CHECK_VBFNLO
 
 HERWIG_CHECK_NJET
 
 HERWIG_CHECK_GOSAM
 
 HERWIG_CHECK_GOSAM_CONTRIB
 
 HERWIG_CHECK_OPENLOOPS
 
 HERWIG_CHECK_MADGRAPH
 
 HERWIG_CHECK_EVTGEN
 HERWIG_CHECK_PYTHIA
 
 HERWIG_COMPILERFLAGS
 
 HERWIG_LOOPTOOLS
 
 HERWIG_PDF_PATH
 
 FASTJET_CHECK_FASTJET
 
 HERWIG_CHECK_ABS_BUG
 
 HERWIG_ENABLE_MODELS
 
 SHARED_FLAG=-shared
 AM_CONDITIONAL(NEED_APPLE_FIXES,
 		[test "xx${host/darwin/foundit}xx" != "xx${host}xx"])
 if test "xx${host/darwin/foundit}xx" != "xx${host}xx"; then
    APPLE_DSO_FLAGS=-Wl,-undefined,dynamic_lookup
    SHARED_FLAG=-bundle
 fi
 AC_SUBST([APPLE_DSO_FLAGS])
 AC_SUBST([SHARED_FLAG])
 
 AC_CONFIG_FILES([UnderlyingEvent/Makefile
 		Models/Makefile
 		Models/StandardModel/Makefile
 		Models/RSModel/Makefile
 		Models/General/Makefile
 		Models/Susy/Makefile
 		Models/Susy/NMSSM/Makefile
 		Models/Susy/RPV/Makefile
 		Models/UED/Makefile
 		Models/LH/Makefile
 		Models/LHTP/Makefile
 		Models/Transplanckian/Makefile
 		Models/Leptoquarks/Makefile
 		Models/Zprime/Makefile
 		Models/TTbAsymm/Makefile
 		Models/Feynrules/Makefile
 		Models/Feynrules/python/Makefile-FR
 		Models/ADD/Makefile
 		Models/Sextet/Makefile
 		Decay/Makefile
 		Decay/FormFactors/Makefile
 		Decay/Tau/Makefile
 		Decay/Baryon/Makefile
 		Decay/VectorMeson/Makefile
 		Decay/Perturbative/Makefile
                 Decay/ScalarMeson/Makefile
                 Decay/TensorMeson/Makefile
                 Decay/WeakCurrents/Makefile
 		Decay/Partonic/Makefile
 		Decay/General/Makefile	
 		Decay/Radiation/Makefile
 		Decay/EvtGen/Makefile
                 Doc/refman.conf
                 Doc/refman.h
 		PDT/Makefile
 		PDF/Makefile
                 MatrixElement/Makefile
                 MatrixElement/General/Makefile
                 MatrixElement/Lepton/Makefile
                 MatrixElement/Hadron/Makefile
                 MatrixElement/DIS/Makefile
                 MatrixElement/Powheg/Makefile
                 MatrixElement/Gamma/Makefile
                 MatrixElement/Reweighters/Makefile
+                MatrixElement/EW/Makefile
                 MatrixElement/Matchbox/Makefile
                 MatrixElement/Matchbox/Base/Makefile
                 MatrixElement/Matchbox/Utility/Makefile
                 MatrixElement/Matchbox/Phasespace/Makefile
                 MatrixElement/Matchbox/Dipoles/Makefile
                 MatrixElement/Matchbox/InsertionOperators/Makefile
                 MatrixElement/Matchbox/Matching/Makefile
                 MatrixElement/Matchbox/Cuts/Makefile
                 MatrixElement/Matchbox/Scales/Makefile
                 MatrixElement/Matchbox/Scales/MatchboxScale.cc
                 MatrixElement/Matchbox/ColorFull/Makefile
                 MatrixElement/Matchbox/CVolver/Makefile
                 MatrixElement/Matchbox/Builtin/Makefile
                 MatrixElement/Matchbox/Builtin/Amplitudes/Makefile
                 MatrixElement/Matchbox/Tests/Makefile
                 MatrixElement/Matchbox/External/Makefile
                 MatrixElement/Matchbox/External/BLHAGeneric/Makefile
                 MatrixElement/Matchbox/External/VBFNLO/Makefile
                 MatrixElement/Matchbox/External/NJet/Makefile
 		MatrixElement/Matchbox/External/GoSam/Makefile
                 MatrixElement/Matchbox/External/OpenLoops/Makefile
                 MatrixElement/Matchbox/External/MadGraph/Makefile
                 MatrixElement/Matchbox/External/MadGraph/mg2herwig.py
 		Sampling/Makefile
 		Sampling/CellGrids/Makefile
                 Shower/Makefile
 		Shower/Matching/Makefile
 		DipoleShower/Makefile
 		DipoleShower/Base/Makefile
 		DipoleShower/Kernels/Makefile
 		DipoleShower/Kinematics/Makefile
 		DipoleShower/Utility/Makefile
 		DipoleShower/AlphaS/Makefile
                 Utilities/Makefile
                 Utilities/XML/Makefile
                 Utilities/Statistics/Makefile
                 Hadronization/Makefile
 		lib/Makefile
 		include/Makefile
 		src/Makefile
 		src/defaults/Makefile
 		src/snippets/Makefile
 		src/Matchbox/Makefile
 		src/herwig-config
 		Doc/Makefile
 		Doc/HerwigDefaults.in
 		Looptools/Makefile
 		Analysis/Makefile
 		src/Makefile-UserModules
 		src/defaults/Analysis.in
 		src/defaults/MatchboxDefaults.in
 		src/defaults/Decays.in
 		src/defaults/decayers.in
 		src/defaults/setup.gosam.in
 		src/Matchbox/LO-DefaultShower.in
 		src/Matchbox/LO-DipoleShower.in
 		src/Matchbox/MCatLO-DefaultShower.in
 		src/Matchbox/MCatLO-DipoleShower.in
 		src/Matchbox/LO-NoShower.in
 		src/Matchbox/MCatNLO-DefaultShower.in
 		src/Matchbox/MCatNLO-DipoleShower.in
 		src/Matchbox/NLO-NoShower.in
 		src/Matchbox/Powheg-DefaultShower.in
 		src/Matchbox/Powheg-DipoleShower.in
 		Contrib/Makefile
 		Contrib/make_makefiles.sh
 		Tests/Makefile
 		Makefile])
 
 AC_CONFIG_LINKS([Doc/BSMlibs.in:Doc/BSMlibs.in])
 AC_CONFIG_FILES([Doc/fixinterfaces.pl],[chmod +x Doc/fixinterfaces.pl])
 
 HERWIG_OVERVIEW
 
 AC_CONFIG_COMMANDS([summary],[cat config.herwig])
 
 AC_OUTPUT