diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
--- a/MatrixElement/EW/GroupInvariants.h
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -1,328 +1,364 @@
 // -*- 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>
 #include <boost/numeric/ublas/matrix.hpp>
 
 namespace Herwig {
 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$ 
    */
   inline double C_A(unsigned int N) {
     return N !=1 ? double(N) : 0.;
   }
 
   /**
    *  The \f$SU(N)\f$ \f$C_F\f$ 
    */
   inline 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$
    */
   inline double C_d(unsigned int N) {
     return (double(N*N)-4.)/double(N);
   }
 
   /**
    *  The \f$SU(N)\f$\f$C_1\f$ 
    */
   inline double C_1(unsigned int N) {
     double N2(N*N);
     return 0.25*(N2-1.0)/N2;
   }
 
   /**
    *  \f$T_F\f$
    */
   inline 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$
    */
   inline double t_S(unsigned int, bool ) {
     return 0.5;
   }
 
   /**
    * / Number of complex scalars in the fundamental rep. of SU(N)/U(1)
    */
   inline 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)
    */
   inline 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);
 
   inline Complex PlusLog(double arg) {
     static const Complex I(0,1.0);
     if (arg>0.0) 
       return log(arg);
     else if (arg<0.0)
       return log(-arg)+I*Constants::pi;
     else 
       assert(false);
   }
 
   inline Complex MinusLog(double arg) {
     static const Complex I(0,1.0);
     if (arg>0.0) 
       return log(arg);
     else if (arg<0.0)
       return log(-arg)-I*Constants::pi;
     else 
       assert(false);
   }
 
   inline Complex getT(Energy2 s, Energy2 t) {
     return MinusLog(-t/GeV2) - MinusLog(-s/GeV2);
   }
   
   inline Complex getU(Energy2 s, Energy2 u) {
     return MinusLog(-u/GeV2) - MinusLog(-s/GeV2);
   }
 
   inline boost::numeric::ublas::matrix<Complex> Gamma2(Complex U, Complex T) {
     boost::numeric::ublas::matrix<Complex> output(2,2);
     static const Complex I(0,1.0);
     using Constants::pi;
     output(0,0) = (-3.0/2.0)*I*pi + (T+U);
     output(1,1) = (-3.0/2.0)*I*pi;
     output(0,1) = 2.0*(T-U);
     output(1,0) = (3.0/8.0)*(T-U);
     return output;
   }
 
   inline boost::numeric::ublas::matrix<Complex> Gamma2w(Complex U, Complex T) {
     boost::numeric::ublas::matrix<Complex> output =  boost::numeric::ublas::zero_matrix<Complex>(5,5);
     static const Complex I(0,1.0);
     using Constants::pi;
     output(0,0) += -I*pi*11.0/4.0;
     output(0,1) += U-T;
     output(1,0) += 2.0*(U-T);
     output(1,1) += -I*pi*11.0/4.0 + (T+U);
     output(2,2) += -7.0/4.0*I*pi + (U+T);
     output(3,3) += -7.0/4.0*I*pi + (U+T);
     output(4,4) += -3.0/4.0*I*pi;
     return output;
   }
 
   inline boost::numeric::ublas::matrix<Complex> Gamma2Singlet() {
     boost::numeric::ublas::matrix<Complex> output =  boost::numeric::ublas::zero_matrix<Complex>(2,2);
     static const Complex I(0,1.0);
     using Constants::pi;
     output(0,0) = output(1,1) = -3.0/4.0*I*pi;
     return output;
   }
 
   inline Complex Gamma1(double hypercharge) {
     Complex I(0,1.0);
     return -I*Constants::pi*(hypercharge*hypercharge);
   }
   
   inline Complex Gamma1(double y1, double y2, Complex T, Complex U) {
     Complex I(0,1.0);
     return -I*Constants::pi*(y1*y1+y2*y2) + 2.0*y1*y2*(T-U);
   }
   
   inline Complex Gamma1(double y1, double y2, double y3, double y4,
 			Complex T, Complex U) {
     Complex I(0,1.0);
     return -I*Constants::pi*(y1*y1+y2*y2+y3*y3+y4*y4)/2.0 + 
       (y1*y4+y2*y3)*T - (y1*y3+y2*y4)*U;
   }
 
+  inline boost::numeric::ublas::matrix<Complex> Gamma3(Complex U, Complex T) {
+    boost::numeric::ublas::matrix<Complex> output =  boost::numeric::ublas::zero_matrix<Complex>(2,2);
+    static const Complex I(0,1.0);
+    using Constants::pi;
+    output(0,0) = -(8.0/3.0)*I*pi;
+    output(1,1) = -(8.0/3.0)*I*pi;
+    output(0,0) += (7.0/3.0)*T + (2.0/3.0)*U;
+    output(0,1) = 2.0*(T-U);
+    output(1,0) = (4.0/9.0)*(T-U);
+    return output;
+  }
+
+  inline boost::numeric::ublas::matrix<Complex> Gamma3g(Complex U, Complex T) {
+    boost::numeric::ublas::matrix<Complex> output =  boost::numeric::ublas::zero_matrix<Complex>(3,3);
+    static const Complex I(0,1.0);
+    using Constants::pi;
+    output(0,2) = U-T;
+    output(1,1) = 3.0/2.0*(T+U);
+    output(1,2) = 3.0/2.0*(U-T);
+    output(2,0) = 2.0*(U-T);
+    output(2,1) = 5.0/6.0*(U-T);
+    output(2,2) = 3.0/2.0*(T+U);
+    for (unsigned int i=0; i<3; i++) {
+      output(i,i) += -13.0/3.0*I*pi;
+    }
+    return output;
+  }
+
+  inline boost::numeric::ublas::matrix<Complex> Gamma3Singlet() {
+    boost::numeric::ublas::matrix<Complex> output =  boost::numeric::ublas::zero_matrix<Complex>(2,2);
+    static const Complex I(0,1.0);
+    using Constants::pi;
+    output(0,0) = output(1,1) = -4.0/3.0*I*pi;
+    return output;
+  }
+
   /**
    * Number of fermion generations (only used in gauge boson HighCMatching)
    */
   inline double n_g() { return 3.0; }
 
   /**
    * Number of complex scalars in the fundamental rep. of SU(N)
    */
   inline double nSWeyl(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.0;
       else assert(false);
     }
   }
 
   /**
    * Number of Weyl Fermions in the fundamental rep. of SU(N)
    */
   inline double nFWeyl(unsigned int N, bool high) {
     if(high) {
       if(N==2 || N==3) return 12.0;
       else assert(false);
     }
     else {
       if(N==3)      return 10.0;
       else if(N==1) return  2.0;
       else assert(false);
     }
   }
 
   inline double TFWeyl(unsigned int) {
     return 0.5;
   }
 
   inline double tSWeyl(unsigned int) {
     return 0.5;
   }
 
   inline Complex WFunction(Energy mu, Energy2 s) {
     using Constants::pi;
     assert(abs(s)>ZERO);
     Complex ln = MinusLog(-s/(mu*mu));
     return (-1.0*ln*ln + 3.0*ln+pi*pi/6.0-8.0);
   }
 
   /**
    * \fX_N\f% function, v is either t or u
    */
   inline Complex XNFunction(unsigned int N, Energy mu, Energy2 s, Energy2 v) {
     using Constants::pi;
     assert(abs(s)>ZERO);
     Complex ls = MinusLog(-s/(mu*mu));
     return (2.0*C_F(N)*WFunction(mu,s) + 
   	    C_A(N)*(2.0*ls*ls - 2.0*MinusLog((s+v)/(mu*mu))*ls - 
   		   11.0/3.0*ls + pi*pi + 85.0/9.0) + 
   	    (2.0/3.0*ls - 10.0/9.0) * TFWeyl(N) * nFWeyl(N,true) + 
   	    (1.0/3.0*ls - 8.0/9.0) * TFWeyl(N) * nSWeyl(N,true));
   }
 
   /**
    *  \f$\Pi_1\f$ function
    */
   inline Complex PI1_function(Energy mu, Energy2 s) {
     assert(abs(s)>ZERO);
     return ((41.0/6.0)*MinusLog(-s/(mu*mu))-104.0/9.0);
   }
 
   /**
    *  \f$\tilde{f}\f$ function, v is either t or u
    */
   inline Complex fTildeFunction(Energy mu, Energy2 s, Energy2 v) {
     using Constants::pi;
     assert(abs(s)>ZERO);
     Complex ls  = MinusLog(-s/GeV2), lv = MinusLog(-v/GeV2);
     Complex lsv = MinusLog((s+v)/GeV2);
     return (-2.0*double(s/(s+v))*(lv-ls) + 
    	    double(s*(s+2.0*v)/((s+v)*(s+v))) * ((lv-ls)*(lv-ls) + pi*pi) + 
   	    4.0*MinusLog(-s/(mu*mu))*(lv-lsv));
   }
 }
 }
 
 #endif // HERWIG_GroupInvariants_H 
diff --git a/MatrixElement/EW/SoftSudakov.cc b/MatrixElement/EW/SoftSudakov.cc
--- a/MatrixElement/EW/SoftSudakov.cc
+++ b/MatrixElement/EW/SoftSudakov.cc
@@ -1,92 +1,881 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the SoftSudakov class.
 //
 
 #include "SoftSudakov.h"
 #include "ThePEG/Interface/ClassDocumentation.h"
 #include "ThePEG/EventRecord/Particle.h"
 #include "ThePEG/Repository/UseRandom.h"
 #include "ThePEG/Repository/EventGenerator.h"
 #include "ThePEG/Utilities/DescribeClass.h"
 #include "GroupInvariants.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 using namespace GroupInvariants;
 
 SoftSudakov::SoftSudakov() : K_ORDER_(3) {}
 
 SoftSudakov::~SoftSudakov() {}
 
 IBPtr SoftSudakov::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr SoftSudakov::fullclone() const {
   return new_ptr(*this);
 }
 
 void SoftSudakov::persistentOutput(PersistentOStream & os) const {
   os << K_ORDER_;
 }
 
 void SoftSudakov::persistentInput(PersistentIStream & is, int) {
   is >> K_ORDER_;
 }
 
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<SoftSudakov,Interfaced>
 describeHerwigSoftSudakov("Herwig::SoftSudakov", "HwMEEW.so");
 
 void SoftSudakov::Init() {
 
   static ClassDocumentation<SoftSudakov> documentation
     ("The SoftSudakov class implements the soft EW Sudakov");
 
 }
 
 InvEnergy SoftSudakov::operator ()(Energy mu) const {
   // Include K-factor Contributions (Cusps):
   GaugeContributions cusp = cuspContributions(mu,K_ORDER_,high_);
   Complex gamma = cusp.SU3*G3_(row_,col_) + cusp.SU2*G2_(row_,col_) + cusp.U1*G1_(row_,col_);
   if (real_) {
     return gamma.real()/mu;
   }
   else {
     return gamma.imag()/mu;
   }
 }
 
 boost::numeric::ublas::matrix<Complex> 
 SoftSudakov::evaluateSoft(boost::numeric::ublas::matrix<Complex> & G3,
 			  boost::numeric::ublas::matrix<Complex> & G2,
 			  boost::numeric::ublas::matrix<Complex> & G1,
 			  Energy mu_h, Energy mu_l, bool high) {
   assert( G3.size1() == G2.size1() && G3.size1() == G1.size1() &&
 	  G3.size2() == G2.size2() && G3.size2() == G1.size2() && 
 	  G3.size1() == G3.size2());
   G3_ = G3;
   G2_ = G2;
   G1_ = G1;
   high_ = high;
   unsigned int NN = G3_.size1();
   // gamma is the matrix to be numerically integrated to run the coefficients.
   boost::numeric::ublas::matrix<Complex> gamma(NN,NN);
   for(row_=0;row_<NN;++row_) {
     for(col_=0;col_<NN;++col_) {
       real_ = true;
       gamma(row_,col_).real(integrator_.value(*this,mu_h,mu_l));
       real_ = false;
       gamma(row_,col_).imag(integrator_.value(*this,mu_h,mu_l));
     }
   }
   // Resummed:
   //return gamma.exp();
   // Fixed Order:
   return boost::numeric::ublas::identity_matrix<Complex>(NN,NN) + gamma;
 }
+
+boost::numeric::ublas::matrix<Complex>
+SoftSudakov::lowEnergyRunning(Energy EWScale, Energy lowScale, 
+			      Energy2 s, Energy2 t, Energy2 u, 
+			      Herwig::EWProcess::Process process) {
+  using namespace EWProcess;
+  using Constants::pi;
+  static const Complex I(0,1.0);
+  Complex T = getT(s,t), U = getU(s,u);
+  boost::numeric::ublas::matrix<Complex> G1, G2, G3;
+  unsigned int numBrokenGauge;
+  switch (process) {
+  case QQQQ:
+  case QQQQiden:
+  case QtQtQQ:
+    {
+      numBrokenGauge = 12;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
+      for (unsigned int i=0; i<numBrokenGauge/2; i++) {
+	G3(i,i)     += gam3(0,0);
+	G3(i,i+6)   += gam3(0,1);
+	G3(i+6,i)   += gam3(1,0);
+	G3(i+6,i+6) += gam3(1,1);
+      }
+      G1(0,0) = G1(6,6)   = Gamma1(2.0/3.0,2.0/3.0,2.0/3.0,2.0/3.0,T,U);
+      G1(1,1) = G1(7,7)   = Gamma1(-1.0/3.0,-1.0/3.0,2.0/3.0,2.0/3.0,T,U);
+      G1(2,2) = G1(8,8)   = Gamma1(2.0/3.0,2.0/3.0,-1.0/3.0,-1.0/3.0,T,U);
+      G1(3,3) = G1(9,9)   = Gamma1(-1.0/3.0,-1.0/3.0,-1.0/3.0,-1.0/3.0,T,U);
+      G1(4,4) = G1(10,10) = Gamma1(-1.0/3.0,2.0/3.0,2.0/3.0,-1.0/3.0,T,U);
+      G1(5,5) = G1(11,11) = Gamma1(2.0/3.0,-1.0/3.0,-1.0/3.0,2.0/3.0,T,U);
+    }
+    break;
+  case QQUU:
+  case QtQtUU:
+  case QQtRtR:
+    {
+      numBrokenGauge = 4;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
+      for (unsigned int i=0; i<numBrokenGauge/2; i++) {
+	G3(i,i)     += gam3(0,0);
+	G3(i,i+2)   += gam3(0,1);
+	G3(i+2,i)   += gam3(1,0);
+	G3(i+2,i+2) += gam3(1,1);
+      }
+      G1(0,0) = G1(2,2) = Gamma1(2.0/3.0,2.0/3.0,T,U);
+      G1(1,1) = G1(3,3) = Gamma1(2.0/3.0,-1.0/3.0,T,U);
+    }
+    break;
+  case QQDD:
+  case QtQtDD:
+    {
+      numBrokenGauge = 4;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
+      for (unsigned int i=0; i<numBrokenGauge/2; i++) {
+	G3(i,i) += gam3(0,0);
+	G3(i,i+2) += gam3(0,1);
+	G3(i+2,i) += gam3(1,0);
+	G3(i+2,i+2) += gam3(1,1);
+      }
+      G1(0,0) = G1(2,2) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
+      G1(1,1) = G1(3,3) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
+    }
+    break;
+  case QQLL:
+    {
+      numBrokenGauge = 6;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam3s = Gamma3Singlet()(0,0);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+            G3(i,i) = gam3s;
+      }
+      G1(0,0) = Gamma1(2.0/3.0,2.0/3.0,0.0,0.0,T,U);
+      G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,0.0,0.0,T,U);
+      G1(2,2) = Gamma1(2.0/3.0,2.0/3.0,-1.0,-1.0,T,U);
+      G1(3,3) = Gamma1(-1.0/3.0,-1.0/3.0,-1.0,-1.0,T,U);
+      G1(4,4) = Gamma1(-1.0/3.0,2.0/3.0,0.0,-1.0,T,U);
+      G1(5,5) = Gamma1(2.0/3.0,-1.0/3.0,-1.0,0.0,T,U);
+    }
+    break;
+  case QQEE:
+    {
+      numBrokenGauge = 2;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G1 *= 0.0; G2 *= 0.0; G3 *= 0.0;
+      Complex gam3s = Gamma3Singlet()(0,0);
+      for (unsigned int i=0; i<2; i++) {
+	G3(i,i) += gam3s;
+      }
+      G1(0,0) = Gamma1(2.0/3.0,-1.0,T,U);
+      G1(1,1) = Gamma1(-1.0/3.0,-1.0,T,U);
+    }
+    break;
+  case UUUU:
+  case UUUUiden:
+  case tRtRUU:
+    numBrokenGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G1 *= 0.0; G2 *= 0.0; G3 *= 0.0;
+    G3 = Gamma3(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,2.0/3.0,T,U);
+    break;
+  case UUDD:
+  case tRtRDD:
+    numBrokenGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = Gamma3(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
+    break;
+  case UULL:
+    numBrokenGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3(0,0) = G3(1,1) = Gamma3Singlet()(0,0);
+    G1(0,0) = Gamma1(2.0/3.0,0.0,T,U);
+    G1(1,1) = Gamma1(2.0/3.0,-1.0,T,U);
+    break;
+  case UUEE:
+    numBrokenGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G1(0,0) = Gamma1(2.0/3.0,-1.0,T,U);
+    break;
+  case DDDD:
+  case DDDDiden:
+    numBrokenGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = Gamma3(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
+    break;
+  case DDLL:
+    numBrokenGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3(0,0) = G3(1,1) = Gamma3Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0/3.0,0.0,T,U);
+    G1(1,1) = Gamma1(-1.0/3.0,-1.0,T,U);
+    break;
+  case DDEE:
+    numBrokenGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0/3.0,-1.0,T,U);
+    break;
+  case LLLL:
+  case LLLLiden:
+    numBrokenGauge = 6;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G1(0,0) = Gamma1(0.0,0.0,0.0,0.0,T,U);
+    G1(1,1) = Gamma1(-1.0,-1.0,0.0,0.0,T,U);
+    G1(2,2) = Gamma1(0.0,0.0,-1.0,-1.0,T,U);
+    G1(3,3) = Gamma1(-1.0,-1.0,-1.0,-1.0,T,U);
+    G1(4,4) = Gamma1(-1.0,0.0,0.0,-1.0,T,U);
+    G1(5,5) = Gamma1(0.0,-1.0,-1.0,0.0,T,U);
+    break;
+    
+  case LLEE:
+    numBrokenGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G1(0,0) = Gamma1(0.0,-1.0,T,U);
+    G1(1,1) = Gamma1(-1.0,-1.0,T,U);
+    break;
+  case EEEE:
+  case EEEEiden:
+    numBrokenGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G1(0,0) = Gamma1(-1.0,-1.0,T,U);
+    break;
+  case QQWW:
+    {
+      numBrokenGauge = 20;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam3s = Gamma3Singlet()(0,0);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+	G3(i,i) = gam3s;
+      }
+      G1(0,0) = Gamma1(2./3.,2./3.,-1.,-1.,T,U);
+      G1(1,1) = Gamma1(2./3.,2./3.,1.,1.,T,U);
+      G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+      G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+      G1(4,4) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+      G1(5,5) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+      G1(6,6) = Gamma1(-1./3.,-1./3.,-1.,-1.,T,U);
+      G1(7,7) = Gamma1(-1./3.,-1./3.,1.,1.,T,U);
+      G1(8,8) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+      G1(9,9) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+      G1(10,10) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+      G1(11,11) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+      G1(12,12) = Gamma1(-1./3.,2./3.,0.,-1.,T,U);
+      G1(13,13) = Gamma1(-1./3.,2./3.,0.,-1.,T,U);
+      G1(14,14) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
+      G1(15,15) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
+      G1(16,16) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
+      G1(17,17) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
+      G1(18,18) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
+      G1(19,19) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
+    }
+    break;
+  case QQPhiPhi:
+    {
+      numBrokenGauge = 14;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam3s = Gamma3Singlet()(0,0);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+	G3(i,i) = gam3s;
+      }
+      G1(0,0) = Gamma1(2./3.,2./3.,1.,1.,T,U);
+      G1(1,1) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+      G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+      G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+      G1(4,4) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+      G1(5,5) = Gamma1(-1./3.,-1./3.,1.,1.,T,U);
+      G1(6,6) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+      G1(7,7) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+      G1(8,8) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+      G1(9,9) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+      G1(10,10) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
+      G1(11,11) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
+      G1(12,12) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
+      G1(13,13) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
+    }
+    break;
+  case QQWG:
+    numBrokenGauge = 6;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    for (unsigned int i=0; i<numBrokenGauge; i++) {
+      G3(i,i) = -17.0/6.0*I*pi + 3.0/2.0*(U+T);
+    }
+    G1(0,0) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
+    G1(1,1) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
+    G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+    G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+    G1(4,4) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+    G1(5,5) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+    break;
+  case QQBG:
+    numBrokenGauge = 4;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    for (unsigned int i=0; i<numBrokenGauge; i++) {
+      G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
+    }
+    G1(0,0) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+    G1(1,1) = Gamma1(2./3.,2./3.,0.,0.,T,U);
+    G1(2,2) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+    G1(3,3) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
+    break;
+  case QQGG:
+  case QtQtGG:
+    {
+      numBrokenGauge = 6;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      boost::numeric::ublas::matrix<Complex> gam3g = Gamma3g(U,T);
+      for(unsigned int ix=0;ix<3;++ix) {
+	for(unsigned int iy=0;iy<3;++iy) {
+	  G3(ix  ,iy  ) = gam3g(ix,iy);
+	  G3(ix+3,iy+3) = gam3g(ix,iy);
+	}
+      }
+      G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(2./3.,0.,T,U);
+      G1(3,3) = G1(4,4) = G1(5,5) = Gamma1(-1./3.,0.,T,U);
+    }
+    break;
+  case LLWW:
+    numBrokenGauge = 20;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G1(0,0) = Gamma1(0.,0.,-1.,-1.,T,U);
+    G1(1,1) = Gamma1(0.,0.,1.,1.,T,U);
+    G1(2,2) = Gamma1(0.,0.,0.,0.,T,U);
+    G1(3,3) = Gamma1(0.,0.,0.,0.,T,U);
+    G1(4,4) = Gamma1(0.,0.,0.,0.,T,U);
+    G1(5,5) = Gamma1(0.,0.,0.,0.,T,U);
+    G1(6,6) = Gamma1(-1.,-1.,-1.,-1.,T,U);
+    G1(7,7) = Gamma1(-1.,-1.,1.,1.,T,U);
+    G1(8,8) = Gamma1(-1.,-1.,0.,0.,T,U);
+    G1(9,9) = Gamma1(-1.,-1.,0.,0.,T,U);
+    G1(10,10) = Gamma1(-1.,-1.,0.,0.,T,U);
+    G1(11,11) = Gamma1(-1.,-1.,0.,0.,T,U);
+    G1(12,12) = Gamma1(-1.,0.,0.,-1.,T,U);
+    G1(13,13) = Gamma1(-1.,0.,0.,-1.,T,U);
+    G1(14,14) = Gamma1(-1.,0.,1.,0.,T,U);
+    G1(15,15) = Gamma1(-1.,0.,1.,0.,T,U);
+    G1(16,16) = Gamma1(0.,-1.,-1.,0.,T,U);
+    G1(17,17) = Gamma1(0.,-1.,-1.,0.,T,U);
+    G1(18,18) = Gamma1(0.,-1.,0.,1.,T,U);
+    G1(19,19) = Gamma1(0.,-1.,0.,1.,T,U);
+    break;
+  case LLPhiPhi:
+    numBrokenGauge = 14;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G1 *= 0.0; G2 *= 0.0; G3 *= 0.0;
+    G1(0,0) = Gamma1(0.,0.,1.,1.,T,U);
+    G1(1,1) = Gamma1(0.,0.,0.,0.,T,U);
+    G1(2,2) = Gamma1(0.,0.,0.,0.,T,U);
+    G1(3,3) = Gamma1(0.,0.,0.,0.,T,U);
+    G1(4,4) = Gamma1(0.,0.,0.,0.,T,U);
+    G1(5,5) = Gamma1(-1.,-1.,1.,1.,T,U);
+    G1(6,6) = Gamma1(-1.,-1.,0.,0.,T,U);
+    G1(7,7) = Gamma1(-1.,-1.,0.,0.,T,U);
+    G1(8,8) = Gamma1(-1.,-1.,0.,0.,T,U);
+    G1(9,9) = Gamma1(-1.,-1.,0.,0.,T,U);
+    G1(10,10) = Gamma1(-1.,0.,1.,0.,T,U);
+    G1(11,11) = Gamma1(-1.,0.,1.,0.,T,U);
+    G1(12,12) = Gamma1(0.,-1.,0.,1.,T,U);
+    G1(13,13) = Gamma1(0.,-1.,0.,1.,T,U);
+    break;
+  case UUBB:
+    {
+      numBrokenGauge = 4;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam3s = Gamma3Singlet()(0,0);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+	G3(i,i) = gam3s;
+	G1(i,i) = Gamma1(2./3.,0.,T,U);
+      }
+    }
+    break;
+  case UUPhiPhi:
+    {
+      numBrokenGauge = 5;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam3s = Gamma3Singlet()(0,0);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+	G3(i,i) = gam3s;
+      }
+      G1(0,0) = Gamma1(2.0/3.0,2.0/3.0,1.,1.,T,U);
+      G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(2./3.,0.,T,U);
+    }
+    break;
+  case UUBG:
+    {
+      numBrokenGauge = 2;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam1 = Gamma1(2./3.,0.,T,U);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+	G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
+	G1(i,i) = gam1;
+      }
+    }
+    break;
+  case UUGG:
+  case tRtRGG:
+    numBrokenGauge = 3;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = Gamma3g(U,T);
+    G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(2./3.,0.,T,U);
+    break;
+  case DDBB:
+    {
+      numBrokenGauge = 4;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam3s = Gamma3Singlet()(0,0);
+      Complex gam1  = Gamma1(-1./3.,0.,T,U);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+	G3(i,i) = gam3s;
+	G1(i,i) = gam1;
+      }
+    }
+    break;
+  case DDPhiPhi:
+    {
+      numBrokenGauge = 5;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam3s = Gamma3Singlet()(0,0);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+	G3(i,i) = gam3s;
+      }
+      G1(0,0) = Gamma1(-1.0/3.0,-1.0/3.0,1.,1.,T,U);
+      G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(-1./3.,0.,T,U);
+    }
+    break;
+  case DDBG:
+    numBrokenGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    for (unsigned int i=0; i<numBrokenGauge; i++) {
+      G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
+    }
+    G1(0,0) = G1(1,1) = Gamma1(-1./3.,0.,T,U);
+    break;
+  case DDGG:
+    numBrokenGauge = 3;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = Gamma3g(U,T);
+    G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(-1./3.,0.,T,U);
+    break;
+  case EEBB:
+    {
+      numBrokenGauge = 4;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+      Complex gam1 = Gamma1(-1.,0.,T,U);
+      for (unsigned int i=0; i<numBrokenGauge; i++) {
+	G1(i,i) = gam1;
+      };
+    }
+    break;
+  case EEPhiPhi:
+    numBrokenGauge = 5;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+    G1(0,0) = Gamma1(-1.,1.,T,U);
+    G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(-1.,0.,T,U);
+    break;         
+  default:
+    assert(false);
+  }
+   
+  if (EWScale==lowScale) {
+    return boost::numeric::ublas::identity_matrix<Complex>(G1.size1());
+  }
+  else {
+    return evaluateSoft(G3,G2,G1,EWScale,lowScale,false);
+  }
+}
+
+boost::numeric::ublas::matrix<Complex>
+SoftSudakov::highEnergyRunning(Energy highScale, Energy EWScale,
+			       Energy2 s, Energy2 t, Energy2 u,
+			       Herwig::EWProcess::Process process) {
+  using namespace EWProcess;
+  using Constants::pi;
+  static const Complex I(0,1.0);
+  Complex T = getT(s,t), U = getU(s,u);
+  boost::numeric::ublas::matrix<Complex> G1,G2,G3;
+  unsigned int numGauge;
+  switch (process) {
+  case QQQQ:
+  case QQQQiden:
+  case QtQtQQ:
+    {
+      numGauge = 4;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+      boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
+      G3(0,0) += gam3(0,0);
+      G3(0,2) += gam3(0,1);
+      G3(2,0) += gam3(1,0);
+      G3(2,2) += gam3(1,1);
+      G3(1,1) += gam3(0,0);
+      G3(1,3) += gam3(0,1);
+      G3(3,1) += gam3(1,0);
+      G3(3,3) += gam3(1,1);
+      boost::numeric::ublas::matrix<Complex> gam2 = Gamma2(U,T);
+      G2(0,0) += gam2(0,0);
+      G2(0,1) += gam2(0,1);
+      G2(1,0) += gam2(1,0);
+      G2(1,1) += gam2(1,1);
+      G2(2,2) += gam2(0,0);
+      G2(2,3) += gam2(0,1);
+      G2(3,2) += gam2(1,0);
+      G2(3,3) += gam2(1,1);
+      G1(0,0) = G1(1,1) = G1(2,2) = G1(3,3) = Gamma1(1.0/6.0,1.0/6.0,T,U);
+    }
+    break;
+  case QQUU:
+  case QtQtUU:
+  case QQtRtR:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = Gamma3(U,T);
+    G2 = Gamma2Singlet();
+    G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,1.0/6.0,T,U);
+    break;
+  case QQDD:
+  case QtQtDD:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = Gamma3(U,T);
+    G2 = Gamma2Singlet();
+    G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,1.0/6.0,T,U);
+    break;
+  case QQLL:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = Gamma3Singlet();
+    G2 = Gamma2(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(-1.0/2.0,1.0/6.0,T,U);
+    break;
+  case QQEE:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge); 
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G2(0,0) = Gamma2Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0,1.0/6.0,T,U);
+    break;
+  case UUUU:
+  case UUUUiden:
+  case tRtRUU:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = Gamma3(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,2.0/3.0,T,U);
+    break;
+  case UUDD:
+  case tRtRDD:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = Gamma3(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
+    break;
+  case UULL:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G2(0,0) = Gamma2Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0/2.0,2.0/3.0,T,U);
+    break;
+  case UUEE:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0,2.0/3.0,T,U);
+    break;
+  case DDDD:
+  case DDDDiden:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = Gamma3(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
+    break;
+  case DDLL:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G2(0,0) = Gamma2Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0/2.0,-1.0/3.0,T,U);
+    break;
+  case DDEE:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0,-1.0/3.0,T,U);
+    break;
+  case LLLL:
+  case LLLLiden:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = Gamma2(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(-1.0/2.0,-1.0/2.0,T,U);
+    break;
+  case LLEE:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2(0,0) = Gamma2Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0,-1.0/2.0,T,U);
+    break;
+  case EEEE:
+  case EEEEiden:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G1(0,0) = Gamma1(-1.0,-1.0,T,U);
+    break;
+  case QQWW:
+    {
+      numGauge = 5;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+      Complex gam3s = Gamma3Singlet()(0,0);
+      for (unsigned int i=0; i<5; i++) {
+	G3(i,i) = gam3s;
+	G1(i,i) = Gamma1(1.0/6.0);
+      }
+      G2 = Gamma2w(U,T);
+    }
+    break;
+  case QQPhiPhi:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = Gamma3Singlet();
+    G2 = Gamma2(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(1.0/2.0,1.0/6.0,T,U);
+    break;
+  case QQWG:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = -17.0/6.0*I*pi + 3.0/2.0*(U+T);
+    G2(0,0) = -7.0/4.0*I*pi + (U+T);
+    G1(0,0) = Gamma1(1.0/6.0);
+    break;
+  case QQBG:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
+    G2(0,0) = -3.0/4.0*I*pi;
+    G1(0,0) = Gamma1(1.0/6.0);
+    break;
+  case QQGG:
+  case QtQtGG:
+    {
+      numGauge = 3;
+      G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+      G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+      G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+      G3 = Gamma3g(T,U);
+      Complex gam2s = Gamma2Singlet()(0,0);
+      for (unsigned int i=0; i<3; i++) {
+	G2(i,i) = gam2s;
+	G1(i,i) = Gamma1(1.0/6.0);
+      }
+    }
+    break;
+  case LLWW:
+    numGauge = 5;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    for (unsigned int i=0; i<5; i++) {
+      G1(i,i) = Gamma1(-1.0/2.0);
+    }
+    G2 = Gamma2w(U,T);
+    break;
+  case LLPhiPhi:
+    numGauge = 2;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = Gamma2(U,T);
+    G1(0,0) = G1(1,1) = Gamma1(1.0/2.0,-1.0/2.0,T,U);
+    break;
+  case UUBB:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G1(0,0) = Gamma1(2.0/3.0);
+    break;
+  case UUPhiPhi:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G2(0,0) = Gamma2Singlet()(0,0);
+    G1(0,0) = Gamma1(1.0/2.0,2.0/3.0,T,U);
+    break;
+  case UUBG:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
+    G1(0,0) = Gamma1(2.0/3.0);
+    break;
+  case UUGG:
+  case tRtRGG:
+    numGauge = 3;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = Gamma3g(U,T);
+    for (unsigned int i=0; i<3; i++) {
+      G1(i,i) = Gamma1(2.0/3.0);
+    }
+    break;
+  case DDBB:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G1(0,0) = Gamma1(-1.0/3.0);
+    break;
+  case DDPhiPhi:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = Gamma3Singlet()(0,0);
+    G2(0,0) = Gamma2Singlet()(0,0);
+    G1(0,0) = Gamma1(1.0/2.0,-1.0/3.0,T,U);
+    break;
+  case DDBG:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
+    G1(0,0) = Gamma1(-1.0/3.0);
+    break;
+  case DDGG:
+    numGauge = 3;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge); 
+    G3 = Gamma3g(U,T);
+    for (unsigned int i=0; i<3; i++) {
+      G1(i,i) = Gamma1(-1.0/3.0);
+    }
+    break;
+  case EEBB:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G1(0,0) = Gamma1(-1.0);
+    break;
+  case EEPhiPhi:
+    numGauge = 1;
+    G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+    G2(0,0) = Gamma2Singlet()(0,0);
+    G1(0,0) = Gamma1(1.0/2.0,-1.0,T,U);
+    break;
+  default:
+    assert(false);
+  }
+  return evaluateSoft(G3,G2,G1,highScale,EWScale,true);
+}
diff --git a/MatrixElement/EW/SoftSudakov.h b/MatrixElement/EW/SoftSudakov.h
--- a/MatrixElement/EW/SoftSudakov.h
+++ b/MatrixElement/EW/SoftSudakov.h
@@ -1,162 +1,181 @@
 // -*- C++ -*-
 #ifndef Herwig_SoftSudakov_H
 #define Herwig_SoftSudakov_H
 //
 // This is the declaration of the SoftSudakov class.
 //
 
 #include "ThePEG/Interface/Interfaced.h"
 #include "Herwig/Utilities/GSLIntegrator.h"
 #include <boost/numeric/ublas/matrix.hpp>
+#include "EWProcess.h"
 #include "SoftSudakov.fh"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
 /**
  * Here is the documentation of the SoftSudakov class.
  *
  * @see \ref SoftSudakovInterfaces "The interfaces"
  * defined for SoftSudakov.
  */
 class SoftSudakov: public Interfaced {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * The default constructor.
    */
   SoftSudakov();
 
   /**
    * The destructor.
    */
   virtual ~SoftSudakov();
   //@}
 
 public:
+
+  /**
+   *  Low energy soft evolution
+   */
+  boost::numeric::ublas::matrix<Complex>
+  lowEnergyRunning(Energy EWScale, Energy lowScale, 
+		   Energy2 s, Energy2 t, Energy2 u,
+		   Herwig::EWProcess::Process process);
+
+  /**
+   *  Evalaute the high energy running as a matrix
+   */
+  boost::numeric::ublas::matrix<Complex>
+  highEnergyRunning(Energy highScale, Energy EWScale,
+		    Energy2 s, Energy2 t, Energy2 u,
+		    Herwig::EWProcess::Process process);
+
+protected:
       
   /**
-   *  Evaluate bthe soft evolution 
+   *  Evaluate the soft evolution 
    */
   boost::numeric::ublas::matrix<Complex> evaluateSoft(boost::numeric::ublas::matrix<Complex> & G3,
 						      boost::numeric::ublas::matrix<Complex> & G2,
 						      boost::numeric::ublas::matrix<Complex> & G1,
 						      Energy mu_h, Energy mu_l, bool high);
 
 public:
 
   /**
    *  The operator to be integrated
    */
   InvEnergy operator ()(Energy mu) const;
   /** Argument type for GaussianIntegrator */
   typedef Energy ArgType;
   /** Return type for GaussianIntegrator */
   typedef InvEnergy ValType;
 
 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.
    */
   SoftSudakov & operator=(const SoftSudakov &);
 
 private:
 
   /**
    *  Order for K
    */
   unsigned int K_ORDER_;
 
   /**
    *  Integrator
    */
   GSLIntegrator integrator_;
 
   /**
    *  Whether doing the high or low scale contribution
    */
   bool high_;
 
   /**
    *  Column
    */
   unsigned int row_; 
 
   /**
    *  Row
    */
   unsigned int col_; 
 
   /**
    *
    */
   bool real_; 
 
   /**
    *
    */
   boost::numeric::ublas::matrix<Complex> G1_;
 
   /**
    *
    */
   boost::numeric::ublas::matrix<Complex> G2_;
 
   /**
    *
    */
   boost::numeric::ublas::matrix<Complex> G3_;
 };
 
 }
 
 #endif /* Herwig_SoftSudakov_H */