diff --git a/MatrixElement/EW/CollinearSudakov.cc b/MatrixElement/EW/CollinearSudakov.cc
--- a/MatrixElement/EW/CollinearSudakov.cc
+++ b/MatrixElement/EW/CollinearSudakov.cc
@@ -1,2084 +1,2075 @@
 // -*- C++ -*-
 //
 // This is the implementation of the non-inlined, non-templated member
 // functions of the CollinearSudakov class.
 //
 
 #include "CollinearSudakov.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 "ElectroWeakReweighter.h"
 #include "ThePEG/Persistency/PersistentOStream.h"
 #include "ThePEG/Persistency/PersistentIStream.h"
 
 using namespace Herwig;
 
 namespace {
 
 void DuplicateColumn0(boost::numeric::ublas::matrix<Complex> &orig) {
   for (unsigned int i=0; i<orig.size1(); i++) {
     for (unsigned int j=1; j<orig.size2(); j++) {
       orig(i,j) = orig(i,0);
     }
   }
 }
 
 }
 
 CollinearSudakov::CollinearSudakov() : K_ORDER_(3),
 				       B_ORDER_(2)
 {}
 
 CollinearSudakov::~CollinearSudakov() {}
 
 IBPtr CollinearSudakov::clone() const {
   return new_ptr(*this);
 }
 
 IBPtr CollinearSudakov::fullclone() const {
   return new_ptr(*this);
 }
 
 void CollinearSudakov::persistentOutput(PersistentOStream & os) const {
   os << K_ORDER_ << B_ORDER_;
 }
 
 void CollinearSudakov::persistentInput(PersistentIStream & is, int) {
   is >> K_ORDER_ >> B_ORDER_;
 }
 
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<CollinearSudakov,Interfaced>
 describeHerwigCollinearSudakov("Herwig::CollinearSudakov", "HwMEEW.so");
 
 void CollinearSudakov::Init() {
 
   static ClassDocumentation<CollinearSudakov> documentation
     ("The CollinearSudakov class implements the collinear evolution");
 
 
 }
 
 
 Complex CollinearSudakov::highScaleIntegral(bool SU3, bool SU2, double Y,
 					    Energy2 s, Energy mu_h, Energy mu_l,
 					    bool fermion, bool longitudinal,
 					    double yukFactor) {
   SU3_          = SU3;
   SU2_          = SU2;
   Y_            = Y;
   s_            = s;
   fermion_      = fermion;
   longitudinal_ = longitudinal;
   yukFactor_    = yukFactor;
   // perform the integral
   Complex result;
   high_ = true;
   // real piece
   real_ = true;
   result.real(integrator_.value(*this,mu_h,mu_l));
   // imaginary piece
   real_ = false;
   result.imag(integrator_.value(*this,mu_h,mu_l));
   // return the answer
   return exp(result);
 }
 
 Complex  CollinearSudakov::lowScaleIntegral(bool SU3, double Q, Energy2 s, 
 					    Energy mu_h, Energy mu_l, bool fermion, 
 					    double boostFactor) {
   SU3_          = SU3;
   Q_            = Q;
   s_            = s;
   fermion_      = fermion;
   boostFactor_  = boostFactor;
   // perform the integral
   Complex result;
   high_ = false;
   // real piece
   real_ = true;
   result.real(integrator_.value(*this,mu_h,mu_l));
   // imaginary piece
   real_ = false;
   result.imag(integrator_.value(*this,mu_h,mu_l));
   // return the answer
   return exp(result);
 }
 
 InvEnergy CollinearSudakov::highScaleIntegrand(Energy mu) const {
   using namespace GroupInvariants;
   using Constants::pi;
   Complex gamma = 0.0;
   // Include K-factor Contributions (Cusps):
   GaugeContributions cusp = cuspContributions(mu,K_ORDER_,high_);
   // Include B-factors (B1 for U1, B2 for SU2, B3 for SU3):
   GaugeContributions nonCusp = BContributions(mu,B_ORDER_,fermion_,longitudinal_);
   // common log
   Complex plog = PlusLog(s_/sqr(mu));
   // SU(3)
   if(SU3_) {
     if(fermion_)
       gamma += C_F(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
     else
       gamma += (C_A(3))*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
   }
   // SU(2)
   if(SU2_) {
     if (fermion_ || longitudinal_ )
       gamma += (C_F(2))*cusp.SU2*0.5*plog + 0.5*nonCusp.SU2;
     else
       gamma += (C_A(2))*cusp.SU2*0.5*plog + 0.5*nonCusp.SU2;
   }
   
   if (fermion_ || longitudinal_ ) {
     gamma += sqr(Y_)*(cusp.U1*0.5*plog + 0.5*nonCusp.U1);
   }
   else {
     // U(1) Gauge boson
     if (!SU3_ && !SU2_ && abs(Y_)<0.001) {
       gamma += 0.5*nonCusp.U1;
     }
   }
   // top Yukawa piece
   double y_t = ElectroWeakReweighter::coupling()->y_t(mu);
   gamma += yukFactor_*sqr(y_t)/(16.0*sqr(pi));
   // return the answer
   return real_ ? gamma.real()/mu : gamma.imag()/mu;
 }
 
 InvEnergy CollinearSudakov::lowScaleIntegrand(Energy mu) const {
   using namespace GroupInvariants;
   using Constants::pi;
   Complex gamma = 0.0;
   // Include K-factor Contributions (Cusps):
   GaugeContributions cusp = cuspContributions(mu,K_ORDER_,false);
   // Include B-factors (B1 for U1, B2 for SU2, B3 for SU3):
   GaugeContributions nonCusp = BContributionsLow(mu,B_ORDER_,fermion_,boostFactor_);
   // common log
   Complex plog = PlusLog(s_/sqr(mu));
   Complex blog(0.);
   if(boostFactor_ >0.001) blog = PlusLog(4.0*boostFactor_);
   // SU(3)
   if (SU3_) {
     if (fermion_) {
       // not a bHQET top quark field
       if (abs(boostFactor_)<0.001)
  	gamma += C_F(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
       else
   	gamma += C_F(3)*cusp.SU3*0.5*blog + 0.5*nonCusp.SU3;
     }
     else {
       gamma += C_A(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
     }
   }
   // fermions
   if (fermion_) {
     // not a bHQET top quark field
     if (abs(boostFactor_)<0.001)
       gamma += sqr(Q_)*(cusp.U1*0.5*plog + 0.5*nonCusp.U1);
     else
       gamma += sqr(Q_)*(cusp.U1*0.5*blog + 0.5*nonCusp.U1);
   }
   else {
     // i.e., not a fermion, not a bHQ, not a gluon => photon
     if (abs(boostFactor_)<0.001 && !SU3_)
       gamma += 0.5*nonCusp.U1;
     // i.e., W treated as a bHQET field
     else if (abs(boostFactor_)>0.001) 
       gamma += sqr(Q_)*(cusp.U1*0.5*blog + 0.5*nonCusp.U1);
   }
   // return the answer
   return real_ ? gamma.real()/mu : gamma.imag()/mu;
 }
 
 void CollinearSudakov::evaluateHighScale(Energy highScale, Energy EWScale, Energy2 S) {
   double yCoeffQt(0.), yCoefftR(0.), yCoeffPhi(0.);
   if (K_ORDER_ >= 1) {
     yCoeffQt = 0.5;
     yCoefftR = 1.0;
     yCoeffPhi = 3.0;
   }
   highColW_   = highScaleIntegral(false,true ,0.0   ,S,highScale,EWScale,false,false,0.0);
   highColB_   = highScaleIntegral(false,false,0.0   ,S,highScale,EWScale,false,false,0.0);
   highColG_   = highScaleIntegral(true ,false,0.0   ,S,highScale,EWScale,false,false,0.0);
   highColQ_   = highScaleIntegral(true ,true ,1./6. ,S,highScale,EWScale,true,false,0.0);
   highColQt_  = highScaleIntegral(true ,true ,1./6. ,S,highScale,EWScale,true,false,yCoeffQt);
   highColU_   = highScaleIntegral(true ,false,2./3. ,S,highScale,EWScale,true,false,0.0);
   highColtR_  = highScaleIntegral(true ,false,2./3. ,S,highScale,EWScale,true,false,yCoefftR);
   highColD_   = highScaleIntegral(true ,false,-1./3.,S,highScale,EWScale,true,false,0.0);
   highColL_   = highScaleIntegral(false,true ,-1./2.,S,highScale,EWScale,true,false,0.0);
   highColE_   = highScaleIntegral(false,false,-1.   ,S,highScale,EWScale,true,false,0.0);
   highColPhi_ = highScaleIntegral(false,true ,1./2. ,S,highScale,EWScale,false,true,yCoeffPhi);
 }
 
 void CollinearSudakov::evaluateLowScale(Energy EWScale, Energy lowScale, Energy2 S) {
   lowColW_ = lowScaleIntegral(false,1.    ,S,EWScale,lowScale,false,
 			      S/sqr(2.*ElectroWeakReweighter::coupling()->mW()));
   lowColA_ = lowScaleIntegral(false,0.    ,S,EWScale,lowScale,false,0.0);
   lowColG_ = lowScaleIntegral(true ,0.    ,S,EWScale,lowScale,false,0.0);
   lowColU_ = lowScaleIntegral(true ,2./3. ,S,EWScale,lowScale,true,0.0);
   lowColt_ = lowScaleIntegral(true ,2./3. ,S,EWScale,lowScale,true,
 			      S/sqr(2.*ElectroWeakReweighter::coupling()->mT()));
   lowColD_ = lowScaleIntegral(true ,-1./3.,S,EWScale,lowScale,true,0.0);
   lowColE_ = lowScaleIntegral(false,-1.   ,S,EWScale,lowScale,true,0.0);
 }
 
 void CollinearSudakov::evaluateMatching(Energy EWScale,Energy2 s) {
   using Constants::pi;
   using GroupInvariants::PlusLog;
   static const Complex I(0,1.0);
   // wave function corrections  
   WaveFunctionCorrections WFC = waveFunctionCorrections(EWScale);
    
   double aS   = ElectroWeakReweighter::coupling()->a3(EWScale);
   double aEM  = ElectroWeakReweighter::coupling()->aEM(EWScale);
   double aW   = ElectroWeakReweighter::coupling()->aW(EWScale);
   double aZ   = ElectroWeakReweighter::coupling()->aZ(EWScale);
   double cW2  = ElectroWeakReweighter::coupling()->Cos2thW(EWScale);
   double sW2  = ElectroWeakReweighter::coupling()->Sin2thW(EWScale);
   Energy mW   = ElectroWeakReweighter::coupling()->mW();
   Energy mZ   = ElectroWeakReweighter::coupling()->mZ();
   Energy mH   = ElectroWeakReweighter::coupling()->mH();
   Energy mT   = ElectroWeakReweighter::coupling()->mT();
   double gPHI = ElectroWeakReweighter::coupling()->g_phiPlus(EWScale);
   double y_t  = ElectroWeakReweighter::coupling()->y_t(EWScale);
    
   double lz = log(mZ/EWScale);
   double lw = log(mW/EWScale);
 
    complex<double> F_W = 4.0*lw*0.5*PlusLog(s/EWScale/EWScale)-2.0*lw*lw-2.0*lw-5.0*pi*pi/12.0+1.0;
    complex<double> F_Z = 4.0*lz*0.5*PlusLog(s/EWScale/EWScale)-2.0*lz*lz-2.0*lz-5.0*pi*pi/12.0+1.0;
    
    // Taken from Manohar... along with his formulae for F_tL, F_tR, and F_bL (for 3rd/1st Gen. Wavefunction Differences)
    complex<double> W1 = WFC.fFW0-0.5*WFC.aW0-0.5*WFC.cW0;
    complex<double> W2 = WFC.fF0W-0.5*WFC.a0W;
    complex<double> U1 = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*(WFC.fFZZ-0.5*WFC.aZZ) - 
                         0.5*WFC.cZZ*(ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale) + 
                                      ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)) + 
                         ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*WFC.bZZ;
    complex<double> U2 = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*(WFC.fFZZ-0.5*WFC.aZZ) -
                         0.5*WFC.cZZ*(ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale) + 
                                      ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)) +
                         ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*WFC.bZZ;
    complex<double> HtL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-0.5*log(mH/EWScale)-0.5*lz+
                                                     0.5*WFC.atHH+0.5*WFC.atZZ+WFC.ctHH+WFC.ctZZ+
                                                     WFC.ctW0-WFC.btHH+WFC.btZZ);
    complex<double> HtR = HtL-0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.atW0);
    complex<double> HbL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.at0W);
    
    complex<double> F_tL = (4.0/3.0*aS+4.0/9.0*aEM)/(4.0*pi)*(2.0*log(mT/EWScale)*log(mT/EWScale)-log(mT/EWScale)+
                                                              pi*pi/12.0+2.0) + HtL + 
                            aW/(4.0*pi)*0.5*W1 + aZ/(4.0*pi)*U1;
    
    complex<double> F_tR = (4.0/3.0*aS+4.0/9.0*aEM)/(4.0*pi)*(2.0*log(mT/EWScale)*log(mT/EWScale)-log(mT/EWScale)+
                                                              pi*pi/12.0+2.0) + HtR - 
                            aW/(4.0*pi)*0.25*WFC.cW0 + aZ/(4.0*pi)*U2;
    
    complex<double> F_bL = HbL + aW/(4.0*pi)*0.5*W2;
    
    Complex Dw = CollinearDw(s,EWScale);
    Complex Dz = CollinearDz(s,EWScale);
 
    complex<double> D_C_UL = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*Dz + 0.5*Dw;
    complex<double> D_C_DL = ElectroWeakReweighter::coupling()->g_Ld(EWScale)*ElectroWeakReweighter::coupling()->g_Ld(EWScale)*Dz + 0.5*Dw;
    
    complex<double> D_C_UR = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*Dz;
    complex<double> D_C_DR = ElectroWeakReweighter::coupling()->g_Rd(EWScale)*ElectroWeakReweighter::coupling()->g_Rd(EWScale)*Dz;
    
    complex<double> D_C_nuL = ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*Dz + 0.5*Dw;
    complex<double> D_C_EL = ElectroWeakReweighter::coupling()->g_Le(EWScale)*ElectroWeakReweighter::coupling()->g_Le(EWScale)*Dz + 0.5*Dw;
    complex<double> D_C_ER = ElectroWeakReweighter::coupling()->g_Re(EWScale)*ElectroWeakReweighter::coupling()->g_Re(EWScale)*Dz;
    
    complex<double> D_C_WW = aW/(4.0*pi)*cW2*(F_Z+WFC.fsWZWZ) + 
      aW/(4.0*pi)*cW2*(F_W+WFC.fs1ZW) + aW/(4.0*pi)*sW2*(F_W+WFC.fs10) + 
      aW/(4.0*pi)*sW2*(2.0*lw*lw-
 		      2.0*lw+pi*pi/12.0+2.0) + 0.5*WFC.RW;
    complex<double> D_C_WZ = aW/(4.0*pi)*2.0*(F_W+WFC.fsZW1) + 0.5*WFC.RZ + sqrt(sW2/cW2)*WFC.RAtoZ;
    complex<double> D_C_WA = aW/(4.0*pi)*2.0*(F_W+WFC.fs01) + 0.5*WFC.RA + sqrt(cW2/sW2)*WFC.RZtoA;
    complex<double> D_C_BZ = 0.5*WFC.RZ - sqrt(cW2/sW2)*WFC.RAtoZ;
    complex<double> D_C_BA = 0.5*WFC.RA - sqrt(sW2/cW2)*WFC.RZtoA;
    
    // The WFC.RW and WFC.RZ are used on purpose (instead of WFC.RPhi and WFC.RPhi3, respectively):
    complex<double> D_C_PhiW = aW/(4.0*pi)*0.25*(F_W+WFC.fs1HW) + 
      aW/(4.0*pi)*0.25*(F_W+WFC.fs1ZW) + aZ/(4.0*pi)*gPHI*gPHI*(F_Z+WFC.fsWZWZ) + 
      aW/(4.0*pi)*sW2*(2.0*lw*lw-2.0*lw+pi*pi/12.0+2.0) + 
      0.5*WFC.RW + WFC.EW;
    complex<double> D_C_PhiZ = aW/(4.0*pi)*0.5*(F_W+WFC.fsZW1) + aZ/(4.0*pi)*0.25*(F_Z+WFC.fs1HZ) + 0.5*WFC.RZ + WFC.EZ;
    complex<double> D_C_PhiH = aW/(4.0*pi)*0.5*(F_W+WFC.fsHW1) + aZ/(4.0*pi)*0.25*(F_Z+WFC.fsHZ1) + 0.5*WFC.RH;
    
    complex<double> D_C_GG = aS/(4.0*pi)*2.0/3.0*log(mT/EWScale);
    
    ULcollinearCorr_     = exp(D_C_UL);
    DLcollinearCorr_     = exp(D_C_DL);
    URcollinearCorr_     = exp(D_C_UR);
    DRcollinearCorr_     = exp(D_C_DR);
    
    tLcollinearCorr_     = exp(D_C_UL+F_tL);
    tRcollinearCorr_     = exp(D_C_UR+F_tR);
    bLcollinearCorr_     = exp(D_C_DL+F_bL);
    
    nuLcollinearCorr_    = exp(D_C_nuL);
    ELcollinearCorr_     = exp(D_C_EL);
    ERcollinearCorr_     = exp(D_C_ER);
    
    WtoWcollinearCorr_   = exp(D_C_WW);
    WtoZcollinearCorr_   = exp(D_C_WZ);
    WtoAcollinearCorr_   = exp(D_C_WA);
    BtoZcollinearCorr_   = exp(D_C_BZ);
    BtoAcollinearCorr_   = exp(D_C_BA);
    
    PhitoWcollinearCorr_ = exp(D_C_PhiW);
    PhitoZcollinearCorr_ = exp(D_C_PhiZ);
    PhitoHcollinearCorr_ = exp(D_C_PhiH);
    
    GcollinearCorr_      = exp(D_C_GG);
 }
 
 WaveFunctionCorrections CollinearSudakov::waveFunctionCorrections(Energy EWScale) {
   static const Complex I(0.,1.);
   using Constants::pi;
   double lZ = 2.0*log(ElectroWeakReweighter::coupling()->mZ()/EWScale);
   WaveFunctionCorrections WFC;
   // From Manohar, 2/12/12: (these assume mH=125, mZ=91.1876, mW=80.399)
   WFC.RW    = (0.8410283377963967 - 9.424777961271568*I) + 1.2785863646210789*lZ;
   WFC.RA    = (1.4835982362022198 + 1.855775680704845*pow(10.,-9)*I) - 0.27209907467584415*lZ;
   WFC.RZ    = (1.5114724841549798 - 9.926944419863688*I) + 1.0834802397165764*lZ;
   WFC.RAtoZ = (0.3667485032811715 - 2.2770907130064835*I) - 1.2994544609942593*lZ;
   WFC.RZtoA = -0.2095310079712942 + 0.8320191107808546*lZ;
   WFC.RH    = (12.229832449946716 - 1.7643103462419842*10.0*pow(10.,-12)*I) + 5.309998583664737*lZ;
   WFC.RPhi  = (5.569012418081201 + 1.5439133581417356*0.10*pow(10.,-9)*I) + 5.309998583664737*lZ;
   WFC.RPhi3 = (8.945333042265943 + 5.499309445612249*pow(10.,-12)*I) + 5.309998583664737*lZ;
   WFC.EW    = (3.967645734304811 + 4.712388980384717*I) + 2.238332625165702*lZ;
   WFC.EZ    = (5.916079892937651 + 4.96347220970469*I) + 2.1132591719740788*lZ;
   
   WFC.RW    *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.RA    *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.RZ    *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.RAtoZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.RZtoA *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.RPhi  *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.EW    *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.EZ    *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.RPhi3 *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   WFC.RH    *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
   
   
   WFC.RW    = WFC.RW.real();
   WFC.RA    = WFC.RA.real();
   WFC.RZ    = WFC.RZ.real();
   WFC.RAtoZ = WFC.RAtoZ.real();
   WFC.RZtoA = WFC.RZtoA.real();
   WFC.RPhi  = WFC.RPhi.real();
   WFC.RPhi3 = WFC.RPhi3.real();
   WFC.RH    = WFC.RH.real();
   
   // Also from Manohar, 2/10/12
   WFC.fFW0 = -3.7946842553189453 - 4.709019671388589*I;
   WFC.fF0W = 3.8181790485161176;
   WFC.fFZZ = 1.364503250377989 + 0.*I;
   WFC.aHH  = -0.474396977740686 + 0.*I;
   WFC.aZZ  = -0.6806563210877914 + 0.*I;
   WFC.aW0  = 0.49036102506811907 + 1.9323351450971642*I;
   WFC.a0W  = -1.2184776671065072;
   WFC.bHH  = 1.234775745474991 + 0.*I;
   WFC.bZZ  = 1.7303526426747613 + 0.*I;
   WFC.cHH  = 0.33140608274387473 + 0.*I;
   WFC.cZZ  = 0.4961363208382017 + 0.*I;
   WFC.cW0  = -1.005299829777395 + 1.063048500002757*I;
   WFC.atHH = -0.237198488870343 + 0.*I;
   WFC.atZZ = -0.3403281605438957 + 0.*I;
   WFC.atW0 = 0.24518051253405954 + 0.9661675725485821*I;
   WFC.at0W = -0.6092388335532536;
   WFC.ctHH = 0.16570304137193737 + 0.*I;
   WFC.ctZZ = 0.24806816041910085 + 0.*I;
   WFC.ctW0 = -0.5026499148886975 + 0.5315242500013785*I;
   WFC.btHH = -0.30869393636874776 + 0.*I;
   WFC.btZZ = -0.4325881606686903 + 0.*I;
   
   WFC.fs10   = -2.289868133696459;
   WFC.fs1ZW  = 1.8627978596240622;
   WFC.fsWZWZ = 1.1866922529667008;
   WFC.fsZW1  = 1.0840307611156266;
   WFC.fs01   = 2.2898681336964595;
   WFC.fsHW1  = -0.32306745562682404;
   WFC.fsHZ1  = 0.4042992116255279;
   WFC.fs1HW  = 3.330954543719127;
   WFC.fs1HZ  = 2.69768201932412;
   return WFC;
 }
 
 Complex CollinearSudakov::CollinearDw(Energy2 s, Energy EWScale) {
   using Constants::pi;
   using GroupInvariants::PlusLog;
   double lw = log(ElectroWeakReweighter::coupling()->mW()/EWScale);
   //s = s/2.; // This should not be here... I think this is a discrepency with Sascha
   return ElectroWeakReweighter::coupling()->aW(EWScale)/(4.0*pi)*
     (4.0*lw*0.5*PlusLog(s/EWScale/EWScale) - 2.0*lw*lw -
      3.0*lw - 5.0/12.0*pi*pi + 9.0/4.0);
 }
 
 Complex CollinearSudakov::CollinearDz(Energy2 s, Energy EWScale) {
   using Constants::pi;
   using GroupInvariants::PlusLog;
   double lz = log(ElectroWeakReweighter::coupling()->mZ()/EWScale);
   return ElectroWeakReweighter::coupling()->aZ(EWScale)/(4.0*pi)*
     (4.0*lz*0.5*PlusLog(s/EWScale/EWScale) - 2.0*lz*lz -
      3.0*lz - 5.0/12.0*pi*pi + 9.0/4.0);
 }
 
 namespace {
 
 void writeLine(ofstream & file, vector<Energy> x, vector<double> y,
 	       string name,string title,
 	       string colour, string style) {
   file << "# BEGIN HISTO1D "+name +"\n";
   file << "SmoothLine=1\n";
   file << "ErrorBars=0\n";
   file << "Path="+name+"\n";
   file << "Title="+title+"\n";
   file << "LineColor="+colour+"\n";
   file << "LineStyle="+style +"\n";
   file << "# xlow	 xhigh	 val	 errminus	 errplus\n";
   for(unsigned int ix=0;ix<x.size();++ix)
     file << (x[ix]-MeV)/GeV << "\t" << (x[ix]+MeV)/GeV << "\t"
 	 << y[ix] << "\t" << 0. << "\t" << 0. << "\n";
   file << "# END HISTO1D\n";
 }
 
 }
 
 void CollinearSudakov::makePlots() {
   vector<Energy> np;
   vector<double> uL,uR,tL,tR,dL,dR,bL,bR,nuL,eL,eR,Wgamma,Bgamma,g,WT,WL,WZT,BZT,ZL,H;
   Energy mZ = ElectroWeakReweighter::coupling()->mZ();
   for(Energy x=200.*GeV;x<5000.*GeV;x*=1.02) {
     Energy2 s(sqr(x));
     np.push_back(x);
     evaluateHighScale(x,mZ,s);
     evaluateMatching(mZ,s);
     uL    .push_back(real(highColQ_   * ULcollinearCorr_    ));
     uR    .push_back(real(highColU_   * URcollinearCorr_    ));
     tL    .push_back(real(highColQt_  * tLcollinearCorr_    ));
     tR    .push_back(real(highColtR_  * tRcollinearCorr_    ));
     dL    .push_back(real(highColQ_   * DLcollinearCorr_    ));
     dR    .push_back(real(highColD_   * DRcollinearCorr_    ));
     bL    .push_back(real(highColQt_  * bLcollinearCorr_    ));
     bR    .push_back(real(highColD_   * DRcollinearCorr_    ));
     nuL   .push_back(real(highColL_   * nuLcollinearCorr_   ));
     eL    .push_back(real(highColL_   * ELcollinearCorr_    ));
     eR    .push_back(real(highColE_   * ERcollinearCorr_    ));
     Wgamma.push_back(real(highColW_   * WtoAcollinearCorr_  ));
     Bgamma.push_back(real(highColB_   * BtoAcollinearCorr_  ));
     g     .push_back(real(highColG_   * GcollinearCorr_     ));
     WT    .push_back(real(highColW_   * WtoWcollinearCorr_  ));
     WL    .push_back(real(highColPhi_ * PhitoWcollinearCorr_));
     WZT   .push_back(real(highColW_   * WtoZcollinearCorr_  ));
     BZT   .push_back(real(highColB_   * BtoZcollinearCorr_  ));
     ZL    .push_back(real(highColPhi_ * PhitoZcollinearCorr_));
     H     .push_back(real(highColPhi_ * PhitoHcollinearCorr_));
   }
   ofstream fig1a("fig1a.dat");
   fig1a << "#BEGIN PLOT\n";
   fig1a << "LegendOnly=/uL /uR /tL /tR\n";
   fig1a << "DrawOnly=/uL /uR /tL /tR\n";
   fig1a << "Title=Figure 1a\n";
   fig1a << "RatioPlot=0\n";
   fig1a << "LogX=1\n";
   fig1a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig1a << "YLabel=u, t\n";
   fig1a << "Legend=1\n";
   fig1a << "XMin=250.\n";
   fig1a << "YMin=0.7\n";
   fig1a << "YMax=1.05\n";
   fig1a << "# END PLOT\n";
   writeLine(fig1a,np,uL,"/uL","$u_L$","green","dotted"   );
   writeLine(fig1a,np,uR,"/uR","$u_R$","cyan" ,"solid"    );
   writeLine(fig1a,np,tL,"/tL","$t_L$","red"  ,"dashed"   );
   writeLine(fig1a,np,tR,"/tR","$t_R$","blue" ,"dotdashed");
   fig1a.close();
   ofstream fig1b("fig1b.dat");
   fig1b << "#BEGIN PLOT\n";
   fig1b << "LegendOnly=/dL /dR /bL /bR\n";
   fig1b << "DrawOnly=/dL /dR /bL /bR\n";
   fig1b << "Title=Figure 1b\n";
   fig1b << "RatioPlot=0\n";
   fig1b << "LogX=1\n";
   fig1b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig1b << "YLabel=d, b\n";
   fig1b << "Legend=1\n";
   fig1b << "XMin=250.\n";
   fig1b << "YMin=0.7\n";
   fig1b << "YMax=1.05\n";
   fig1b << "# END PLOT\n";
   writeLine(fig1b,np,dL,"/dL","$d_L$","green","dotted"   );
   writeLine(fig1b,np,dR,"/dR","$d_R$","cyan" ,"solid"    );
   writeLine(fig1b,np,bL,"/bL","$b_L$","red"  ,"dashed"   );
   writeLine(fig1b,np,bR,"/bR","$b_R$","blue" ,"dotdashed");
   fig1b.close();
   ofstream fig2("fig2.dat");
   fig2 << "#BEGIN PLOT\n";
   fig2 << "LegendOnly=/uL /uR /dL /dR\n";
   fig2 << "DrawOnly=/uL /uR /dL /dR\n";
   fig2 << "Title=Figure 2\n";
   fig2 << "RatioPlot=0\n";
   fig2 << "LogX=1\n";
   fig2 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig2 << "YLabel=u, d\n";
   fig2 << "Legend=1\n";
   fig2 << "XMin=250.\n";
   fig2 << "YMin=0.7\n";
   fig2 << "YMax=1.05\n";
   fig2 << "# END PLOT\n";
   writeLine(fig2,np,uL,"/uL","$u_L$","green","dotted"   );
   writeLine(fig2,np,uR,"/uR","$u_R$","cyan" ,"solid"    );
   writeLine(fig2,np,dL,"/dL","$d_L$","red"  ,"dashed"   );
   writeLine(fig2,np,dR,"/dR","$d_R$","blue" ,"dotdashed");
   fig2.close();
   ofstream fig3("fig3.dat");
   fig3 << "#BEGIN PLOT\n";
   fig3 << "LegendOnly=/nuL /eL /eR\n";
   fig3 << "DrawOnly=/nuL /eL /eR\n";
   fig3 << "Title=Figure 3\n";
   fig3 << "RatioPlot=0\n";
   fig3 << "LogX=1\n";
   fig3 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig3 << "YLabel=$\\nu$, $e$\n";
   fig3 << "Legend=1\n";
   fig3 << "XMin=250.\n";
   fig3 << "YMin=0.9\n";
   fig3 << "YMax=1.05\n";
   fig3 << "# END PLOT\n";
   writeLine(fig3,np,nuL,"/nuL","$\\nu_L$","blue","dashed");
   writeLine(fig3,np, eL,"/eL" ,"$e_L$"   ,"red" ,"dotted");
   writeLine(fig3,np, eR,"/eR" ,"$e_R$"   ,"red" ,"solid" );
   fig3.close();
   ofstream fig5("fig5.dat");
   fig5 << "#BEGIN PLOT\n";
   fig5 << "LegendOnly=/g /Wgamma /Bgamma\n";
   fig5 << "DrawOnly=/g /Wgamma /Bgamma\n";
   fig5 << "Title=Figure 5\n";
   fig5 << "RatioPlot=0\n";
   fig5 << "LogX=1\n";
   fig5 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig5 << "YLabel=$\\gamma$, g\n";
   fig5 << "Legend=1\n";
   fig5 << "XMin=250.\n";
   fig5 << "YMin=0.7\n";
   fig5 << "YMax=1.05\n";
   fig5 << "# END PLOT\n";
   writeLine(fig5,np,g     ,"/g"     ,"$g$"           ,"blue","dashed");
   writeLine(fig5,np,Wgamma,"/Wgamma","$W\\to\\gamma$","red" ,"solid" );
   writeLine(fig5,np,Bgamma,"/Bgamma","$B\\to\\gamma$","red" ,"dotted");
   fig5.close();
   ofstream fig6a("fig6a.dat");
   fig6a << "#BEGIN PLOT\n";
   fig6a << "LegendOnly=/WT /WL\n";
   fig6a << "DrawOnly=/WT /WL\n";
   fig6a << "Title=Figure 6a\n";
   fig6a << "RatioPlot=0\n";
   fig6a << "LogX=1\n";
   fig6a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig6a << "YLabel=$Z_L$, $Z_T$, $H$\n";
   fig6a << "Legend=1\n";
   fig6a << "XMin=250.\n";
   fig6a << "YMin=0.7\n";
   fig6a << "YMax=1.05\n";
   fig6a << "# END PLOT\n";
   writeLine(fig6a,np,WT,"/WT","$W_T$","red"  ,"solid");
   writeLine(fig6a,np,WL,"/WL","$W_L$","blue" ,"dashed" );
   fig6a.close();
-
-
   ofstream fig6b("fig6b.dat");
   fig6b << "#BEGIN PLOT\n";
   fig6b << "LegendOnly=/WZT /BZT /ZL /H\n";
   fig6b << "DrawOnly=/WZT /BZT /ZL /H\n";
   fig6b << "Title=Figure 6b\n";
   fig6b << "RatioPlot=0\n";
   fig6b << "LogX=1\n";
   fig6b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig6b << "YLabel=d, b\n";
   fig6b << "Legend=1\n";
   fig6b << "XMin=250.\n";
   fig6b << "YMin=0.7\n";
   fig6b << "YMax=1.05\n";
   fig6b << "# END PLOT\n";
   writeLine(fig6b,np,WZT,"/WZT","$W\\to Z_T$","red"  ,"solid"    );
   writeLine(fig6b,np,BZT,"/BZT","$B\\to Z_T$","red"  ,"dotted"   );
   writeLine(fig6b,np,ZL ,"/ZL ","$Z_L$"      ,"blue" ,"dashed"   );
   writeLine(fig6b,np,H  ,"/H  ","$H$"        ,"green","dotdashed");
   fig6b.close();
 
 
   np.clear();
   vector<double> e30,e50,q30,q50,g30,g50;
   for(Energy x=200.*GeV;x<5000.*GeV;x*=1.02) {
     Energy2 s(sqr(x));
     np.push_back(x);
     evaluateLowScale(mZ,30.*GeV,s);
     e30.push_back(real(lowColE_));
     q30.push_back(real(lowColU_));
     g30.push_back(real(lowColG_));
     evaluateLowScale(mZ,50.*GeV,s);
     e50.push_back(real(lowColE_));
     q50.push_back(real(lowColU_));
     g50.push_back(real(lowColG_));
   }
   ofstream fig4a("fig4a.dat");
   fig4a << "#BEGIN PLOT\n";
   fig4a << "LegendOnly=/e30 /e50\n";
   fig4a << "DrawOnly=/e30 /e50\n";
   fig4a << "Title=Figure 4a\n";
   fig4a << "RatioPlot=0\n";
   fig4a << "LogX=1\n";
   fig4a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig4a << "YLabel=e\n";
   fig4a << "Legend=1\n";
   fig4a << "XMin=250.\n";
   fig4a << "YMin=0.7\n";
   fig4a << "YMax=1.05\n";
   fig4a << "# END PLOT\n";
   writeLine(fig4a,np,e30,"/e30","e $(\\mu_f=30\\,\\mathrm{GeV})$","red" ,"solid" );
   writeLine(fig4a,np,e50,"/e50","e $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dashed");
   fig4a.close();
   ofstream fig4b("fig4b.dat");
   fig4b << "#BEGIN PLOT\n";
   fig4b << "LegendOnly=/q30 /q50 /g30 /g50\n";
   fig4b << "DrawOnly=/q30 /q50 /g30 /g50\n";
   fig4b << "Title=Figure 4a\n";
   fig4b << "RatioPlot=0\n";
   fig4b << "LogX=1\n";
   fig4b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
   fig4b << "YLabel=e\n";
   fig4b << "Legend=1\n";
   fig4b << "XMin=250.\n";
   fig4b << "YMin=0.5\n";
   fig4b << "YMax=1.05\n";
   fig4b << "# END PLOT\n";
   writeLine(fig4b,np,q30,"/q30","q $(\\mu_f=30\\,\\mathrm{GeV})$","red" ,"solid" );
   writeLine(fig4b,np,q50,"/q50","q $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dashed");
   writeLine(fig4b,np,g30,"/g30","g $(\\mu_f=30\\,\\mathrm{GeV})$","green" ,"dotted" );
   writeLine(fig4b,np,g50,"/g50","g $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dotdashed");
   fig4b.close();
 }
 
 boost::numeric::ublas::matrix<Complex>
 CollinearSudakov::electroWeakMatching(Energy EWScale, Energy2 s,
 				      Herwig::EWProcess::Process process,
 				      bool oneLoop) {
   using namespace EWProcess;
   // calculate the matching coefficients
   evaluateMatching(EWScale,s);
   // fill the matrix
   boost::numeric::ublas::matrix<Complex> result(1,1);
   switch (process) {         
   case QQQQ:
   case QQQQiden: 
     {
       unsigned int numGauge = 4, numBrokenGauge = 12;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
       result(0,0) = result(6,0) = ULcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_;
       result(3,0) = result(9,0) = DLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_;
       for (int i=0; i<12; i++) {
 	if (i!=0 && i!=3 && i!=6 && i!=9) {
 	  result(i,0) = ULcollinearCorr_*ULcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_;
 	}
       }
       DuplicateColumn0(result);
     }
     break;
   case QtQtQQ:
     {
       unsigned int numGauge = 4, numBrokenGauge = 12;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
       result(0,0) = result(6,0) = ULcollinearCorr_*ULcollinearCorr_*tLcollinearCorr_*tLcollinearCorr_;
       result(3,0) = result(9,0) = DLcollinearCorr_*DLcollinearCorr_*bLcollinearCorr_*bLcollinearCorr_;
       for (int i=0; i<12; i++) {
 	if (i==4 || i==5 || i==10 || i==11) {
 	  result(i,0) = ULcollinearCorr_*tLcollinearCorr_*DLcollinearCorr_*bLcollinearCorr_;
 	}
 	else if (i==1 || i==7) {
 	  result(i,0) = DLcollinearCorr_*DLcollinearCorr_*tLcollinearCorr_*tLcollinearCorr_;
 	}
 	else if (i==2 || i==8) {
 	  result(i,0) = ULcollinearCorr_*ULcollinearCorr_*bLcollinearCorr_*bLcollinearCorr_;
 	}
       }
       DuplicateColumn0(result);
     }
     break;
   case QQUU:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
       result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QtQtUU:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
       result(1,0) = result(3,0) = bLcollinearCorr_*bLcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QQtRtR:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*tRcollinearCorr_*tRcollinearCorr_;
       result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*tRcollinearCorr_*tRcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QQDD:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QtQtDD:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       result(1,0) = result(3,0) = bLcollinearCorr_*bLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QQLL:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 6;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_;
       result(1,0) = nuLcollinearCorr_*nuLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_;
       result(2,0) = ELcollinearCorr_*ELcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_;
       result(3,0) = ELcollinearCorr_*ELcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_;
       result(4,0) = result(5,0) = nuLcollinearCorr_*ELcollinearCorr_*ULcollinearCorr_*DLcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QQEE:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = ULcollinearCorr_*ULcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
       result(1,0) = DLcollinearCorr_*DLcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case UUUU:
   case UUUUiden:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = URcollinearCorr_*URcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case tRtRUU:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = tRcollinearCorr_*tRcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case UUDD:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = URcollinearCorr_*URcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case tRtRDD:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = tRcollinearCorr_*tRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case UULL:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
       result(1,0) = ELcollinearCorr_*ELcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case UUEE:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 1;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = URcollinearCorr_*URcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case DDDD:
   case DDDDiden:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = DRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case DDLL:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       result(1,0) = ELcollinearCorr_*ELcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case DDEE:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 1;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = DRcollinearCorr_*DRcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case LLLL:
   case LLLLiden:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 6;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_;
       result(1,0) = nuLcollinearCorr_*nuLcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_;
       result(2,0) = ELcollinearCorr_*ELcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_;
       result(3,0) = ELcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_;
       result(4,0) = result(5,0) = nuLcollinearCorr_*ELcollinearCorr_*nuLcollinearCorr_*ELcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case LLEE:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
       result(1,0) = ELcollinearCorr_*ELcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case EEEE:
   case EEEEiden:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 1;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = ERcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QQWW:
   case LLWW:
     {
       unsigned int numGauge = 5;
       unsigned int numBrokenGauge = 20;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       for (unsigned int row = 0; row < result.size1(); row++) {
 	for (unsigned int col = 0; col < result.size2(); col++) {
 	  
 	  // Boson Collinear Corr_ections:
 	  if (col==0 || col==1) {
 	    if (row==0 || row==1 || row==6 || row==7) result(row,col) = (WtoWcollinearCorr_*WtoWcollinearCorr_);
 	    if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*WtoZcollinearCorr_);
 	    if (row==3 || row==4 || row==9 || row==10) result(row,col) = (WtoZcollinearCorr_*WtoAcollinearCorr_);
 	    if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*WtoAcollinearCorr_);
 	    if (row==12 || row==14) result(row,col) = (WtoWcollinearCorr_*WtoZcollinearCorr_);
 	    if (row==13 || row==15) result(row,col) = (WtoWcollinearCorr_*WtoAcollinearCorr_);
 	    if (row==16 || row==18) result(row,col) = (WtoWcollinearCorr_*WtoZcollinearCorr_);
 	    if (row==17 || row==19) result(row,col) = (WtoWcollinearCorr_*WtoAcollinearCorr_);
 	  }
 	  if (col==2) {
 	    if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==3 || row==9) result(row,col) = (WtoZcollinearCorr_*BtoAcollinearCorr_);
 	    if (row==4 || row==10) result(row,col) = (WtoAcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*BtoAcollinearCorr_);
 	    if (row==14) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==15) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_);
 	    if (row==16) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==17) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_);
 	  }
 	  if (col==3) {
 	    if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==3 || row==9) result(row,col) = (WtoAcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==4 || row==10) result(row,col) = (WtoZcollinearCorr_*BtoAcollinearCorr_);
 	    if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*BtoAcollinearCorr_);
 	    if (row==12) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==13) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_);
 	    if (row==18) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==19) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_);
 	  }
 	  if (col==4) {
 	    if (row==2 || row==8) result(row,col) = (BtoZcollinearCorr_*BtoZcollinearCorr_);
 	    if (row==3 || row==4 || row==9 || row==10) result(row,col) = (BtoZcollinearCorr_*BtoAcollinearCorr_);
 	    if (row==5 || row==11) result(row,col) = (BtoAcollinearCorr_*BtoAcollinearCorr_);
 	  }
           
 	  // Particle Collinear Corr_ections:
 	  if (process==QQWW) {
 	    if (row<6) result(row,col) *= (ULcollinearCorr_*ULcollinearCorr_);
 	    if ((row>=6)&&(row<12)) result(row,col) *= (DLcollinearCorr_*DLcollinearCorr_);
 	    if (row>=12) result(row,col) *= (ULcollinearCorr_*DLcollinearCorr_);
 	  }
 	  else if (process==LLWW) {
 	    if (row<6) result(row,col) *= (nuLcollinearCorr_*nuLcollinearCorr_);
 	    if ((row>=6)&&(row<12)) result(row,col) *= (ELcollinearCorr_*ELcollinearCorr_);
 	    if (row>=12) result(row,col) *= (nuLcollinearCorr_*ELcollinearCorr_);
 	  }
 	}
       }
     }
     break;
   case QQPhiPhi:
   case LLPhiPhi:
     {
       unsigned int numGauge = 2;
       unsigned int numBrokenGauge = 14;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       for (unsigned int row = 0; row < result.size1(); row++) {
 	
 	// Boson Colinear Corr_ections:
 	if (row==0 || row==5) result(row,0) = (PhitoWcollinearCorr_*PhitoWcollinearCorr_);
 	if (row==1 || row==6) result(row,0) = (PhitoZcollinearCorr_*PhitoZcollinearCorr_);
 	if (row==2 || row==3 || row==7 || row==8) result(row,0) = (PhitoZcollinearCorr_*PhitoHcollinearCorr_);
 	if (row==4 || row==9) result(row,0) = (PhitoHcollinearCorr_*PhitoHcollinearCorr_);
 	if (row==10) result(row,0) = (PhitoWcollinearCorr_*PhitoZcollinearCorr_);
 	if (row==11) result(row,0) = (PhitoWcollinearCorr_*PhitoHcollinearCorr_);
 	if (row==12) result(row,0) = (PhitoWcollinearCorr_*PhitoZcollinearCorr_);
 	if (row==13) result(row,0) = (PhitoWcollinearCorr_*PhitoHcollinearCorr_);
         
 	// Particle Colinear Corr_ections:
 	if (process==QQPhiPhi) {
 	  if (row<5) result(row,0) *= (ULcollinearCorr_*ULcollinearCorr_);
 	  if ((row>=5)&&(row<10)) result(row,0) *= (DLcollinearCorr_*DLcollinearCorr_);
 	  if (row>=10) result(row,0) *= (ULcollinearCorr_*DLcollinearCorr_);
 	}
 	else if (process==LLPhiPhi) {
 	  if (row<5) result(row,0) *= (nuLcollinearCorr_*nuLcollinearCorr_);
 	  if ((row>=5)&&(row<10)) result(row,0) *= (ELcollinearCorr_*ELcollinearCorr_);
 	  if (row>=10) result(row,0) *= (nuLcollinearCorr_*ELcollinearCorr_);
 	}
       }
       DuplicateColumn0(result);
     }
     break;
   case QQWG:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 6;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = ULcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoWcollinearCorr_;
       result(2,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*WtoZcollinearCorr_;
       result(3,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*WtoAcollinearCorr_;
       result(4,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoZcollinearCorr_;
       result(5,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoAcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QQBG:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_;
       result(1,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_;
       result(2,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_;
       result(3,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QQGG:
     {
       unsigned int numGauge = 3;
       unsigned int numBrokenGauge = 6;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
       result(3,0) = result(4,0) = result(5,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case QtQtGG:
     {
       unsigned int numGauge = 3;
       unsigned int numBrokenGauge = 6;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
       result(3,0) = result(4,0) = result(5,0) = bLcollinearCorr_*bLcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case UUBB:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = URcollinearCorr_*URcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_;
       result(1,0) = URcollinearCorr_*URcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_;
       result(2,0) = URcollinearCorr_*URcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_;
       result(3,0) = URcollinearCorr_*URcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case UUPhiPhi:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 5;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = URcollinearCorr_*URcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_;
       result(1,0) = URcollinearCorr_*URcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_;
       result(2,0) = URcollinearCorr_*URcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_;
       result(3,0) = URcollinearCorr_*URcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_;
       result(4,0) = URcollinearCorr_*URcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case UUBG:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_;
       result(1,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case UUGG:
     {
       unsigned int numGauge = 3;
       unsigned int numBrokenGauge = 3;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = result(2,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case tRtRGG:
     {
       unsigned int numGauge = 3;
       unsigned int numBrokenGauge = 3;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = result(2,0) = tRcollinearCorr_*tRcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case DDBB:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = DRcollinearCorr_*DRcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_;
       result(1,0) = DRcollinearCorr_*DRcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_;
       result(2,0) = DRcollinearCorr_*DRcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_;
       result(3,0) = DRcollinearCorr_*DRcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case DDPhiPhi:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 5;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_;
       result(1,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_;
       result(2,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_;
       result(3,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_;
       result(4,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;      
   case DDBG:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 2;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_;
       result(1,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case DDGG:
     {
       unsigned int numGauge = 3;
       unsigned int numBrokenGauge = 3;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = result(1,0) = result(2,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case EEBB:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 4;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = ERcollinearCorr_*ERcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_;
       result(1,0) = ERcollinearCorr_*ERcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_;
       result(2,0) = ERcollinearCorr_*ERcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_;
       result(3,0) = ERcollinearCorr_*ERcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   case EEPhiPhi:
     {
       unsigned int numGauge = 1;
       unsigned int numBrokenGauge = 5;
       result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
       result(0,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_;
       result(1,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_;
       result(2,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_;
       result(3,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_;
       result(4,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_;
       DuplicateColumn0(result);
     }
     break;
   default:
     assert(false);
   }
 
   // This is done at the end instead of the beginning for result.size1() and cols()
   if (!oneLoop) { 
     boost::numeric::ublas::matrix<Complex> OnesMatrix(result.size1(),result.size2());
     for (unsigned int i=0; i<OnesMatrix.size1(); i++) {
       for (unsigned int j=0; j<OnesMatrix.size2(); j++) {
 	OnesMatrix(i,j) = 1.0;
       }
     }
     return OnesMatrix;
   }
-   
-  // Only include the following for the FO calculation:
-  for (unsigned int i=0; i<result.size1(); i++) {
-    for (unsigned int j=0; j<result.size2(); j++) {
-      result(i,j) = 1.0 + log(result(i,j));
-    }
-  }
   
   return result;
 }
 
 boost::numeric::ublas::matrix<Complex>
 CollinearSudakov::highEnergyRunning(Energy highScale, Energy EWScale, Energy2 s,
 				    Herwig::EWProcess::Process process,
 				    bool fixedOrder) {
   using namespace EWProcess;
   // perform the calculation
   evaluateHighScale(highScale,EWScale,s);
   Complex colW(highColW_);
   Complex colB(highColB_); 
   Complex colG(highColG_);
   Complex colQ(highColQ_);
   Complex colQt(highColQt_);
   Complex colU(highColU_);
   Complex coltR(highColtR_);
   Complex colD(highColD_);
   Complex colL(highColL_);
   Complex colE(highColE_);
   Complex colPhi(highColPhi_);
   if (fixedOrder) {
     /* colX not necessarily positive for s = (1000TeV)^2 for the following:
        colW = log(colW.real());
        colB = log(colB.real());
        colPhi = log(colPhi.real());
     */
     colG = log(colG.real());
     colQ = log(colQ.real());
     colQt = log(colQt.real());
     colU = log(colU.real());
     coltR = log(coltR.real());
     colD = log(colD.real());
     colL = log(colL.real());
     colE = log(colE.real());
   }
   // set up the matrix
   boost::numeric::ublas::matrix<Complex> result;
   unsigned int numGauge(0);
   switch (process) {
          
       case QQQQ:
       case QQQQiden:
       case QtQtQQ:
          numGauge = 4;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (process!=QtQtQQ) {
             for (unsigned int i=0; i<numGauge; i++) {
                if (fixedOrder) {
                   result(i,i) = 1.0+colQ+colQ+colQ+colQ;
                }
                else {
                   result(i,i) = colQ*colQ*colQ*colQ;
                }
             }
          }
          else {
             for (unsigned int i=0; i<numGauge; i++) {
                if (fixedOrder) {
                   result(i,i) = 1.0+colQt+colQt+colQ+colQ;
                }
                else {
                   result(i,i) = colQt*colQt*colQ*colQ;
                }
             }
          }
          break;
          
       case QQUU:
       case QtQtUU:
       case QQtRtR:
          numGauge = 2;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (process==QQUU) {
             for (unsigned int i=0; i<numGauge; i++) {
                if (fixedOrder) {
                   result(i,i) = 1.0+colQ+colQ+colU+colU;
                }
                else {
                   result(i,i) = colQ*colQ*colU*colU;
                }
             }
          }
          else if (process==QtQtUU) {
             for (unsigned int i=0; i<numGauge; i++) {
                if (fixedOrder) {
                   result(i,i) = 1.0+colQt+colQt+colU+colU;
                }
                else {
                   result(i,i) = colQt*colQt*colU*colU;
                }
             }
          }
          else if (process==QQtRtR) {
             for (unsigned int i=0; i<numGauge; i++) {
                if (fixedOrder) {
                   result(i,i) = 1.0+colQ+colQ+coltR+coltR;
                }
                else {
                   result(i,i) = colQ*colQ*coltR*coltR;
                }
             }
          }
          break;
          
       case QQDD:
       case QtQtDD:
          numGauge = 2;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (process==QQDD) {
             for (unsigned int i=0; i<numGauge; i++) {
                if (fixedOrder) {
                   result(i,i) = 1.0+colQ+colQ+colD+colD;
                }
                else {
                   result(i,i) = colQ*colQ*colD*colD;
                }
             }
          }
          else {
             for (unsigned int i=0; i<numGauge; i++) {
                if (fixedOrder) {
                   result(i,i) = 1.0+colQt+colQt+colD+colD;
                }
                else {
                   result(i,i) = colQt*colQt*colD*colD;
                }
             }
          }
          break;
          
       case QQLL:
          numGauge = 2;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          for (unsigned int i=0; i<numGauge; i++) {
    	   if (fixedOrder) {
    	     result(i,i) = 1.0+colQ+colQ+colL+colL;
    	   }
    	   else {
    	     result(i,i) = colQ*colQ*colL*colL;
    	   }
          }
          break;
          
    case QQEE:
      numGauge = 1;
      result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
      for (unsigned int i=0; i<numGauge; i++) {
        if (fixedOrder) {
    	 result(i,i) = 1.0+colQ+colQ+colE+colE;
        }
        else {
    	 result(i,i) = colQ*colQ*colE*colE;
        }
      }
      break;
      
    case UUUU:
    case UUUUiden:
    case tRtRUU:
      numGauge = 2;
      result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
      if (process!=tRtRUU) {
        for (unsigned int i=0; i<numGauge; i++) {
    	 if (fixedOrder) {
    	   result(i,i) = 1.0+colU+colU+colU+colU;
    	 }
    	 else {
    	   result(i,i) = colU*colU*colU*colU;
    	 }
        }
      }
      else {
        for (unsigned int i=0; i<numGauge; i++) {
    	 if (fixedOrder) {
    	   result(i,i) = 1.0+coltR+coltR+colU+colU;
    	 }
    	 else {
    	   result(i,i) = coltR*coltR*colU*colU;
    	 }
        }
      }
      break;
      
    case UUDD:
    case tRtRDD:
      numGauge = 2;
      result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
      if (process==UUDD) {
        for (unsigned int i=0; i<numGauge; i++) {
    	 if (fixedOrder) {
    	   result(i,i) = 1.0+colU+colU+colD+colD;
    	 }
    	 else {
    	   result(i,i) = colU*colU*colD*colD;
    	 }
        }
      }
      else {
        for (unsigned int i=0; i<numGauge; i++) {
    	 if (fixedOrder) {
    	   result(i,i) = 1.0+coltR+coltR+colD+colD;
    	 }
    	 else {
    	   result(i,i) = coltR*coltR*colD*colD;
    	 }
        }
      }
      break;
      
    case UULL:
      numGauge = 1;
      result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
      for (unsigned int i=0; i<numGauge; i++) {
        if (fixedOrder) {
    	 result(i,i) = 1.0+colU+colU+colL+colL;
        }
        else {
    	 result(i,i) = colU*colU*colL*colL;
        }
      }
      break;
      
    case UUEE:
      numGauge = 1;
      result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
      for (unsigned int i=0; i<numGauge; i++) {
        if (fixedOrder) {
    	 result(i,i) = 1.0+colU+colU+colE+colE;
        }
        else {
    	 result(i,i) = colU*colU*colE*colE;
        }
      }
      break;
      
    case DDDD:
    case DDDDiden:
      numGauge = 2;
      result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
      for (unsigned int i=0; i<numGauge; i++) {
        if (fixedOrder) {
    	 result(i,i) = 1.0+colD+colD+colD+colD;
        }
        else {
    	 result(i,i) = colD*colD*colD*colD;
        }
      }
      break;
      
    case DDLL:
      numGauge = 1;
      result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
      for (unsigned int i=0; i<numGauge; i++) {
        if (fixedOrder) {
    	 result(i,i) = 1.0+colD+colD+colL+colL;
        }
        else {
    	 result(i,i) = colD*colD*colL*colL;
             }
          }
          break;
          
       case DDEE:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          for (unsigned int i=0; i<numGauge; i++) {
             if (fixedOrder) {
                result(i,i) = 1.0+colD+colD+colE+colE;
             }
             else {
                result(i,i) = colD*colD*colE*colE;
             }
          }
          break;
          
       case LLLL:
       case LLLLiden:
          numGauge = 2;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          for (unsigned int i=0; i<numGauge; i++) {
             if (fixedOrder) {
                result(i,i) = 1.0+colL+colL+colL+colL;
             }
             else {
                result(i,i) = colL*colL*colL*colL;
             }
          }
          break;
          
       case LLEE:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          for (unsigned int i=0; i<numGauge; i++) {
             if (fixedOrder) {
                result(i,i) = 1.0+colL+colL+colE+colE;
             }
             else {
                result(i,i) = colL*colL*colE*colE;
             }
          }
          break;
          
       case EEEE:
       case EEEEiden:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          for (unsigned int i=0; i<numGauge; i++) {
             if (fixedOrder) {
                result(i,i) = 1.0+colE+colE+colE+colE;
             }
             else {
                result(i,i) = colE*colE*colE*colE;
             }
          }
          break;
          
       case QQWW:
          numGauge = 5;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = result(1,1) = 1.0+colQ+colQ+colW+colW;
             result(2,2) = result(3,3) = 1.0+colQ+colQ+colW+colB;
             result(4,4) = 1.0+colQ+colQ+colB+colB;
          }
          else {
             result(0,0) = result(1,1) = colQ*colQ*colW*colW;
             result(2,2) = result(3,3) = colQ*colQ*colW*colB;
             result(4,4) = colQ*colQ*colB*colB;
          }
          break;
          
       case QQPhiPhi:
          numGauge = 2;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = result(1,1) = 1.0+colQ+colQ+colPhi+colPhi;
          }
          else {
             result(0,0) = result(1,1) = colQ*colQ*colPhi*colPhi;
          }
          break;
          
       case QQWG:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colQ+colQ+colW+colG;
          }
          else {
             result(0,0) = colQ*colQ*colW*colG;
          }
          break;
          
       case QQBG:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colQ+colQ+colB+colG;
          }
          else {
             result(0,0) = colQ*colQ*colB*colG;
          }
          break;
          
       case QQGG:
       case QtQtGG:
          numGauge = 3;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (process==QQGG) {
             if (fixedOrder) {
                result(0,0) = result(1,1) = result(2,2) = 1.0+colQ+colQ+colG+colG;
             }
             else {
                result(0,0) = result(1,1) = result(2,2) = colQ*colQ*colG*colG;
             }
          }
          else {
             if (fixedOrder) {
                result(0,0) = result(1,1) = result(2,2) = 1.0+colQt+colQt+colG+colG;
             }
             else {
                result(0,0) = result(1,1) = result(2,2) = colQt*colQt*colG*colG;
             }
          }
          break;
          
       case LLWW:
          numGauge = 5;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = result(1,1) = 1.0+colL+colL+colW+colW;
             result(2,2) = result(3,3) = 1.0+colL+colL+colW+colB;
             result(4,4) = 1.0+colL+colL+colB+colB;
          }
          else {
             result(0,0) = result(1,1) = colL*colL*colW*colW;
             result(2,2) = result(3,3) = colL*colL*colW*colB;
             result(4,4) = colL*colL*colB*colB;
          }
          break;
          
       case LLPhiPhi:
          numGauge = 2;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = result(1,1) = 1.0+colL+colL+colPhi+colPhi;
          }
          else {
             result(0,0) = result(1,1) = colL*colL*colPhi*colPhi;
          }
          break;
          
       case UUBB:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colU+colU+colB+colB;
          }
          else {
             result(0,0) = colU*colU*colB*colB;
          }
          break;
          
       case UUPhiPhi:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colU+colU+colPhi+colPhi;
          }
          else {
             result(0,0) = colU*colU*colPhi*colPhi;
          }
          break;
          
       case UUBG:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colU+colU+colB+colG;
          }
          else {
             result(0,0) = colU*colU*colB*colG;
          }
          break;
          
       case UUGG:
       case tRtRGG:
          numGauge = 3;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (process==UUGG) {
             if (fixedOrder) {
                result(0,0) = result(1,1) = result(2,2) = 1.0+colU+colU+colG+colG;
             }
             else {
                result(0,0) = result(1,1) = result(2,2) = colU*colU*colG*colG;
             }
          }
          else {
             if (fixedOrder) {
                result(0,0) = result(1,1) = result(2,2) = 1.0+coltR+coltR+colG+colG;
             }
             else {
                result(0,0) = result(1,1) = result(2,2) = coltR*coltR*colG*colG;
             }
          }
          break;
          
       case DDBB:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colD+colD+colB+colB;
          }
          else {
             result(0,0) = colD*colD*colB*colB;
          }
          break;
          
       case DDPhiPhi:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colD+colD+colPhi+colPhi;
          }
          else {
             result(0,0) = colD*colD*colPhi*colPhi;
          }
          break;
          
       case DDBG:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colD+colD+colB+colG;
          }
          else {
             result(0,0) = colD*colD*colB*colG;
          }
          break;
          
       case DDGG:
          numGauge = 3;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = result(1,1) = result(2,2) = 1.0+colD+colD+colG+colG;
          }
          else {
             result(0,0) = result(1,1) = result(2,2) = colD*colD*colG*colG;
          }
          break;
          
       case EEBB:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colE+colE+colB+colB;
          }
          else {
             result(0,0) = colE*colE*colB*colB;
          }
          break;
          
       case EEPhiPhi:
          numGauge = 1;
          result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
          if (fixedOrder) {
             result(0,0) = 1.0+colE+colE+colPhi+colPhi;
          }
          else {
             result(0,0) = colE*colE*colPhi*colPhi;
          }
          break;
          
   default:
     assert(false);
   }
    
    return result;
 }
 
 boost::numeric::ublas::matrix<Complex>
 CollinearSudakov::lowEnergyRunning(Energy EWScale, Energy lowScale, Energy2 s,
 				   Herwig::EWProcess::Process process) {
   using namespace EWProcess;
   // evaluate the running
   evaluateLowScale(EWScale,lowScale,s);
 
   Complex colUL = lowColU_;
   Complex colDL = lowColD_;
   Complex colUR = lowColU_;
   Complex colDR = lowColD_;
   
   Complex coltL = lowColt_;
   Complex coltR = lowColt_;
   Complex colbL = lowColD_;
   
   Complex colnuL = 1.0;
   Complex colEL = lowColE_;
   Complex colER = lowColE_;
   
   Complex colW = lowColW_;
   Complex colZ = 1.0;
   Complex colA = lowColA_;
   
   Complex colPhi = lowColW_;
   Complex colPhi3 = 1.0;
   Complex colH = 1.0;
   
   Complex colG = lowColG_;
   
   // calculate the matrix
   boost::numeric::ublas::matrix<Complex> result;
   unsigned int numBrokenGauge;
   switch (process) {
   case QQQQ:
   case QQQQiden:
     numBrokenGauge = 12;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(6,6) = colUL*colUL*colUL*colUL;
     result(3,3) = result(9,9) = colDL*colDL*colDL*colDL;
     for (unsigned int i=0; i<12; i++) {
       if (i!=0 && i!=3 && i!=6 && i!=9) {
 	result(i,i) = colUL*colUL*colDL*colDL;
       }
     }
     break;
   case QtQtQQ:
     numBrokenGauge = 12;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(6,6) = colUL*colUL*coltL*coltL;
     result(3,3) = result(9,9) = colDL*colDL*colbL*colbL;
     for (unsigned int i=0; i<12; i++) {
       if (i==4 || i==5 || i==10 || i==11) {
 	result(i,i) = colUL*coltL*colDL*colbL;
       }
       else if (i==1 || i==7) {
 	result(i,i) = colDL*colDL*coltL*coltL;
       }
       else if (i==2 || i==8) {
 	result(i,i) = colUL*colUL*colbL*colbL;
       }
     }
     break;
   case QQUU:
     numBrokenGauge = 4;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(2,2) = colUL*colUL*colUR*colUR;
     result(1,1) = result(3,3) = colDL*colDL*colUR*colUR;
     break;
   case QtQtUU:
     numBrokenGauge = 4;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(2,2) = coltL*coltL*colUR*colUR;
     result(1,1) = result(3,3) = colbL*colbL*colUR*colUR;
     break;
   case QQtRtR:
     numBrokenGauge = 4;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(2,2) = colUL*colUL*coltR*coltR;
     result(1,1) = result(3,3) = colDL*colDL*coltR*coltR;
     break;
   case QQDD:
     numBrokenGauge = 4;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(2,2) = colUL*colUL*colDR*colDR;
     result(1,1) = result(3,3) = colDL*colDL*colDR*colDR;
     break;
   case QtQtDD:
     numBrokenGauge = 4;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(2,2) = coltL*coltL*colDR*colDR;
     result(1,1) = result(3,3) = colbL*colbL*colDR*colDR;
     break;
   case QQLL:
     numBrokenGauge = 6;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colnuL*colnuL*colUL*colUL;
     result(1,1) = colnuL*colnuL*colDL*colDL;
     result(2,2) = colEL*colEL*colUL*colUL;
     result(3,3) = colEL*colEL*colDL*colDL;
     result(4,4) = result(5,5) = colnuL*colEL*colUL*colDL;
     break;
   case QQEE:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colUL*colUL*colER*colER;
     result(1,1) = colDL*colDL*colER*colER;
     break;
   case UUUU:
   case UUUUiden:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = colUR*colUR*colUR*colUR;
     break;
   case tRtRUU:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = coltR*coltR*colUR*colUR;
     break;
   case UUDD:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = colUR*colUR*colDR*colDR;
     break;
   case tRtRDD:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = coltR*coltR*colDR*colDR;
     break;
   case UULL:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colnuL*colnuL*colUR*colUR;
     result(1,1) = colEL*colEL*colUR*colUR;
     break;
   case UUEE:
     numBrokenGauge = 1;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colUR*colUR*colER*colER;
     break;
   case DDDD:
   case DDDDiden:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = colDR*colDR*colDR*colDR;
     break;
   case DDLL:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colnuL*colnuL*colDR*colDR;
     result(1,1) = colEL*colEL*colDR*colDR;
     break;
   case DDEE:
     numBrokenGauge = 1;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colDR*colDR*colER*colER;
     break;
   case LLLL:
   case LLLLiden:
     numBrokenGauge = 6;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colnuL*colnuL*colnuL*colnuL;
     result(1,1) = colnuL*colnuL*colEL*colEL;
     result(2,2) = colEL*colEL*colnuL*colnuL;
     result(3,3) = colEL*colEL*colEL*colEL;
     result(4,4) = result(5,5) = colnuL*colEL*colnuL*colEL;
     break;
   case LLEE:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colnuL*colnuL*colER*colER;
     result(1,1) = colEL*colEL*colER*colER;
     break;
   case EEEE:
   case EEEEiden:
     numBrokenGauge = 1;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colER*colER*colER*colER;
     break;
   case QQWW:
   case LLWW:
     numBrokenGauge = 20;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     for (unsigned int row = 0; row < result.size1(); row++) {
       // Boson Collinear Corrections:
       if (row==0 || row==1 || row==6 || row==7) result(row,row) = (colW*colW);
       else if (row==2 || row==8) result(row,row) = (colZ*colZ);
       else if (row==3 || row==4 || row==9 || row==10) result(row,row) = (colZ*colA);
       else if (row==5 || row==11) result(row,row) = (colA*colA);
       else if (row==12 || row==14) result(row,row) = (colW*colZ);
       else if (row==13 || row==15) result(row,row) = (colW*colA);
       else if (row==16 || row==18) result(row,row) = (colW*colZ);
       else if (row==17 || row==19) result(row,row) = (colW*colA);
       
       // Particle Collinear Corrections:
       if (process==QQWW) {
 	if (row<6) result(row,row) *= (colUL*colUL);
 	if ((row>=6)&&(row<12)) result(row,row) *= (colDL*colDL);
 	if (row>=12) result(row,row) *= (colUL*colDL);
       }
       else if (process==LLWW) {
 	if (row<6) result(row,row) *= (colnuL*colnuL);
 	if ((row>=6)&&(row<12)) result(row,row) *= (colEL*colEL);
 	if (row>=12) result(row,row) *= (colnuL*colEL);
       }
     }
     break;
   case QQPhiPhi:
   case LLPhiPhi:
     numBrokenGauge = 14;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     for (unsigned int row = 0; row < result.size1(); row++) {
       
       // Boson Colinear Corrections:
       if (row==0 || row==5) result(row,row) = (colPhi*colPhi);
       else if (row==1 || row==6) result(row,row) = (colPhi3*colPhi3);
       else if (row==2 || row==3 || row==7 || row==8) result(row,row) = (colPhi3*colH);
       else if (row==4 || row==9) result(row,row) = (colH*colH);
       else if (row==10) result(row,row) = (colPhi*colPhi3);
       else if (row==11) result(row,row) = (colPhi*colH);
       else if (row==12) result(row,row) = (colPhi*colPhi3);
       else if (row==13) result(row,row) = (colPhi*colH);
       
       // Particle Colinear Corrections:
       if (process==QQPhiPhi) {
 	if (row<5) result(row,row) *= (colUL*colUL);
 	if ((row>=5)&&(row<10)) result(row,row) *= (colDL*colDL);
 	if (row>=10) result(row,row) *= (colUL*colDL);
       }
       else if (process==LLPhiPhi) {
 	if (row<5) result(row,row) *= (colnuL*colnuL);
 	if ((row>=5)&&(row<10)) result(row,row) *= (colEL*colEL);
 	if (row>=10) result(row,row) *= (colnuL*colEL);
       }
     }
     break;
   case QQWG:
     numBrokenGauge = 6;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = colUL*colDL*colG*colW;
     result(2,2) = colUL*colUL*colG*colZ;
     result(3,3) = colUL*colUL*colG*colA;
     result(4,4) = colDL*colDL*colG*colZ;
     result(5,5) = colDL*colDL*colG*colA;
     break;
   case QQBG:
     numBrokenGauge = 4;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colUL*colUL*colG*colZ;
     result(1,1) = colUL*colUL*colG*colA;
     result(2,2) = colDL*colDL*colG*colZ;
     result(3,3) = colDL*colDL*colG*colA;
     break;
   case QQGG:
     numBrokenGauge = 6;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = result(2,2) = colUL*colUL*colG*colG;
     result(3,3) = result(4,4) = result(5,5) = colDL*colDL*colG*colG;
     break;
   case QtQtGG:
     numBrokenGauge = 6;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = result(2,2) = coltL*coltL*colG*colG;
     result(3,3) = result(4,4) = result(5,5) = colbL*colbL*colG*colG;
     break;
       case UUBB:
 	numBrokenGauge = 4;
 	result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
 	result(0,0) = colUR*colUR*colZ*colZ;
 	result(1,1) = colUR*colUR*colZ*colA;
 	result(2,2) = colUR*colUR*colA*colZ;
 	result(3,3) = colUR*colUR*colA*colA;
 	break;
   case UUPhiPhi:
     numBrokenGauge = 5;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colUR*colUR*colPhi*colPhi;
     result(1,1) = colUR*colUR*colPhi3*colPhi3;
     result(2,2) = colUR*colUR*colH*colPhi3;
     result(3,3) = colUR*colUR*colPhi3*colH;
     result(4,4) = colUR*colUR*colH*colH;
     break;
   case UUBG:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colUR*colUR*colG*colZ;
     result(1,1) = colUR*colUR*colG*colA;
     break;
   case UUGG:
     numBrokenGauge = 3;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = result(2,2) = colUR*colUR*colG*colG;
     break;
   case tRtRGG:
     numBrokenGauge = 3;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = result(2,2) = coltR*coltR*colG*colG;
     break;
   case DDBB:
     numBrokenGauge = 4;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colDR*colDR*colZ*colZ;
     result(1,1) = colDR*colDR*colZ*colA;
     result(2,2) = colDR*colDR*colA*colZ;
     result(3,3) = colDR*colDR*colA*colA;
     break;
   case DDPhiPhi:
     numBrokenGauge = 5;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colDR*colDR*colPhi*colPhi;
     result(1,1) = colDR*colDR*colPhi3*colPhi3;
     result(2,2) = colDR*colDR*colH*colPhi3;
     result(3,3) = colDR*colDR*colPhi3*colH;
     result(4,4) = colDR*colDR*colH*colH;
     break;
   case DDBG:
     numBrokenGauge = 2;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colDR*colDR*colG*colZ;
     result(1,1) = colDR*colDR*colG*colA;
     break;
     
   case DDGG:
     numBrokenGauge = 3;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = result(1,1) = result(2,2) = colDR*colDR*colG*colG;
     break;
   case EEBB:
     numBrokenGauge = 4;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colER*colER*colZ*colZ;
     result(1,1) = colER*colER*colZ*colA;
     result(2,2) = colER*colER*colA*colZ;
     result(3,3) = colER*colER*colA*colA;
     break;
   case EEPhiPhi:
     numBrokenGauge = 5;
     result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
     result(0,0) = colER*colER*colPhi*colPhi;
     result(1,1) = colER*colER*colPhi3*colPhi3;
     result(2,2) = colER*colER*colH*colPhi3;
     result(3,3) = colER*colER*colPhi3*colH;
     result(4,4) = colER*colER*colH*colH;
     break;
   default:
     assert(false);
   }
   return result;
 }
 
 
 
 
diff --git a/MatrixElement/EW/ElectroWeakMatching.cc b/MatrixElement/EW/ElectroWeakMatching.cc
--- a/MatrixElement/EW/ElectroWeakMatching.cc
+++ b/MatrixElement/EW/ElectroWeakMatching.cc
@@ -1,874 +1,876 @@
 // -*- C++ -*-
 //
 // ElectroWeakMatching.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
 //
 // Herwig is licenced under version 2 of the GPL, see COPYING for details.
 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
 //
 //
 //
 #include "ElectroWeakMatching.h"
 #include "ElectroWeakReweighter.h"
 #include "GroupInvariants.h"
 #include <boost/numeric/ublas/operation.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 
 using namespace Herwig;
 using namespace ElectroWeakMatching;
 using namespace GroupInvariants;
 using namespace EWProcess;
 
 boost::numeric::ublas::matrix<Complex>
 ElectroWeakMatching::electroWeakMatching(Energy mu,
 					 Energy2 s, Energy2 t, Energy2 u,
 					 Herwig::EWProcess::Process process,
 					 bool oneLoop) {
   static const Complex I(0,1.0);
   using Constants::pi;
   Complex T = getT(s,t);
   Complex U = getU(s,u);
    // Z-Couplings
    double g_Lu = ElectroWeakReweighter::coupling()->g_Lu(mu);
    double g_Ld = ElectroWeakReweighter::coupling()->g_Ld(mu);
    double g_Le = ElectroWeakReweighter::coupling()->g_Le(mu);
    double g_Lnu = ElectroWeakReweighter::coupling()->g_Lnu(mu);
    double g_Ru = ElectroWeakReweighter::coupling()->g_Ru(mu);
    double g_Rd = ElectroWeakReweighter::coupling()->g_Rd(mu);
    double g_Re = ElectroWeakReweighter::coupling()->g_Re(mu);
    double g_W = ElectroWeakReweighter::coupling()->g_W(mu);
    double g_phiPlus = ElectroWeakReweighter::coupling()->g_phiPlus(mu);
    // Weinberg Angle:
    double cos2 = ElectroWeakReweighter::coupling()->Cos2thW(mu);
    double sin2 = 1.0-cos2;
    double cos = sqrt(cos2);
    double sin = sqrt(sin2);
 
    boost::numeric::ublas::matrix<Complex> R0,G2,Dw,Dz;
    
    switch (process) {
    case QQQQ:
    case QQQQiden:
    case QtQtQQ:
      {
        unsigned int numGauge = 4, numBrokenGauge = 12;
-       boost::numeric::ublas::matrix<Complex> R0=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
-       boost::numeric::ublas::matrix<Complex> G2=boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
-       boost::numeric::ublas::matrix<Complex> Dw=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
-       boost::numeric::ublas::matrix<Complex> Dz=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+       R0=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+       G2=boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+       Dw=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+       Dz=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
        R0(0,0) = R0(3,0) = 0.25;
        R0(1,0) = R0(2,0) = -0.25;
        R0(4,0) = R0(5,0) = 0.5;
        R0(6,3) = R0(7,3) = R0(8,3) = R0(9,3) = 1.0;
        R0(6,2) = R0(9,2) = 0.25;
        R0(7,2) = R0(8,2) = -0.25;
        R0(10,2) = R0(11,2) = 0.5;
        if (oneLoop) {
 	 double g11 = g_Lu;
 	 double g12 = g_Ld;
 	 double g21 = g_Lu;
 	 double g22 = g_Ld;
 	 for(unsigned int ix=0;ix<numBrokenGauge;++ix) {
 	   Dw(ix,ix) = 0.5*I*pi;
 	 }
 	 Complex w1 = -0.5*(T-U);
 	 Complex w2 = -0.5*(T+U);
 	 for(unsigned int ix=0;ix<numBrokenGauge;ix+=6) {
 	   Dw(ix+0,ix+0) += w1;
 	   Dw(ix+3,ix+3) += w1;
 	   Dw(ix+1,ix+1) +=-w1;
 	   Dw(ix+2,ix+2) +=-w1;
 	   Dw(ix+4,ix+4) += w2; 
 	   Dw(ix+5,ix+5) += w2;
 	 }
 	 Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
 	 Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
 	 Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
 	 Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U 
 	   - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
 	 for(unsigned int ix=0;ix<numBrokenGauge;ix+=6) {
 	   Dz(ix+0,ix+0) = z1;
 	   Dz(ix+1,ix+1) = z2;
 	   Dz(ix+2,ix+2) = z3;
 	   Dz(ix+3,ix+3) = z4;
 	   Dz(ix+4,ix+4) = z5;
 	   Dz(ix+5,ix+5) = z5;
 	 }  
 	 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);
        }
      }
      break;
    case QQUU:
    case QtQtUU:
    case QQtRtR:
      {
        unsigned int numGauge = 2, numBrokenGauge = 4;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,0) = R0(2,1) = R0(3,1) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Lu;
 	 double g12 = g_Ld;
 	 //double g21 = g_Ru;
 	 double g22 = g_Ru;
          
 	 Complex w1 = 0.25*I*pi;
 	 for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;
 	 
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 for(unsigned int ix=0;ix<numBrokenGauge;ix+=2) {
 	   Dz(ix+0,ix+0) = z1;
 	   Dz(ix+1,ix+1) = z2;
 	 }
 	 G2 = Gamma2Singlet();
        }
      }
      break;
    case QQDD:
    case QtQtDD:
      {
        unsigned int numGauge = 2, numBrokenGauge = 4;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,0) = R0(2,1) = R0(3,1) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Lu;
 	 double g12 = g_Ld;
 	 //double g21 = g_Rd;
 	 double g22 = g_Rd;
          
 	 Complex w1 = 0.25*I*pi;
 	 for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;
          
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 for(unsigned int ix=0;ix<numBrokenGauge;ix+=2) {
 	   Dz(ix+0,ix+0) = z1;
 	   Dz(ix+1,ix+1) = z2;
 	 }
 	 G2 = Gamma2Singlet();
        }
      }
      break;
    case QQLL:
      {
        unsigned int numGauge = 2, numBrokenGauge = 6;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
        R0(0,0) = R0(3,0) = 0.25;
        R0(1,0) = R0(2,0) = -0.25;
        R0(4,0) = R0(5,0) = 0.5;
        if (oneLoop) {
 	 double g11 = g_Lu;
 	 double g12 = g_Ld;
 	 double g21 = g_Lnu;
 	 double g22 = g_Le;
 	 
 	 for (unsigned int i=0; i<6; ++i) {
 	   Dw(i,i) = 0.5*I*pi;
 	 }
 	 Complex w1 = (-1.0/2.0)*(T-U);
 	 Complex w2 = (-1.0/2.0)*(T+U);
 	 Dw(0,0) += w1;
 	 Dw(3,3) += w1;
 	 Dw(1,1) += -1.0*w1;
 	 Dw(2,2) += -1.0*w1;
 	 Dw(4,4) += w2; 
 	 Dw(5,5) += w2;
          
 	 Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
 	 Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
 	 Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
 	 Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U 
 	   - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z2;
 	 Dz(2,2) = z3;
 	 Dz(3,3) = z4;
 	 Dz(4,4) = Dz(5,5) = z5;
          
 	 G2 = Gamma2(U,T);
        }
      }
      break;
    case QQEE:
      {
        unsigned int numGauge = 1, numBrokenGauge = 2;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,0) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Lu;
 	 double g12 = g_Ld;
 	 //double g21 = g_Re;
 	 double g22 = g_Re;
             
 	 Complex w1 = 0.25*I*pi;
 	 Dw(0,0) = Dw(1,1) = w1;
             
          Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z2;
             
 	 G2(0,0) = Gamma2Singlet()(0,0);
        }
      }
      break;
    case UUUU:
    case UUUUiden:
    case tRtRUU:
      {
        unsigned int numGauge = 2, numBrokenGauge = 2;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,1) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Ru;
 	 //double g12 = g_Ru;
 	 //double g21 = g_Ru;
 	 double g22 = g_Ru;
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Dz(0,0) = Dz(1,1) = z1;
        }
      }
      break;
    case UUDD:
    case tRtRDD:
      {
        unsigned int numGauge = 2, numBrokenGauge = 2;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,1) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Ru;
 	 //double g12 = g_Ru;
 	 //double g21 = g_Rd;
 	 double g22 = g_Rd;
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Dz(0,0) = Dz(1,1) = z1;
        }
      }
      break;
    case UULL:
      {
        unsigned int numGauge = 1, numBrokenGauge = 2;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,0) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Lnu;
 	 double g12 = g_Le;
 	 //double g21 = g_Ru;
 	 double g22 = g_Ru;
          
 	 Complex w1 = 0.25*I*pi;
 	 Dw(0,0) = Dw(1,1) = w1;
             
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z2;
             
 	 G2(0,0) = Gamma2Singlet()(0,0);
        }
      }
      break;
    case UUEE:
      {
        unsigned int numGauge = 1, numBrokenGauge = 1;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Ru;
 	 //double g12 = g_Ru;
 	 //double g21 = g_Re;
 	 double g22 = g_Re;
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Dz(0,0) = z1;
        }
      }
      break;
    case DDDD:
    case DDDDiden:
      {
        unsigned int numGauge = 2, numBrokenGauge = 2;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,1) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Rd;
 	 //double g12 = g_Rd;
 	 //double g21 = g_Rd;
 	 double g22 = g_Rd;
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Dz(0,0) = Dz(1,1) = z1;
        }
      }
      break;
    case DDLL:
      {
        unsigned int numGauge = 1, numBrokenGauge = 2;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,0) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Lnu;
 	 double g12 = g_Le;
 	 //double g21 = g_Rd;
 	 double g22 = g_Rd;
          
 	 Complex w1 = 0.25*I*pi;
 	 Dw(0,0) = Dw(1,1) = w1;
          
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z2;
          
 	 G2(0,0) = Gamma2Singlet()(0,0);
        }
      }
      break;
    case DDEE:
      {
        unsigned int numGauge = 1, numBrokenGauge = 1;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0 *= 0.0; Dw = Dz *= 0.0;
        R0(0,0) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Rd;
 	 //double g12 = g_Rd;
 	 //double g21 = g_Re;
 	 double g22 = g_Re;
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Dz(0,0) = z1;
        }
      }
      break;
    case LLLL:
    case LLLLiden:
      {
        unsigned int numGauge = 2, numBrokenGauge = 6;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
        R0(0,0) = R0(3,0) = 0.25;
        R0(1,0) = R0(2,0) = -0.25;
        R0(4,0) = R0(5,0) = 0.5;
        if (oneLoop) {
 	 double g11 = g_Lnu;
 	 double g12 = g_Le;
 	 double g21 = g_Lnu;
 	 double g22 = g_Le;
          
 	 for (int i=0; i<6; i++) {
 	   Dw(i,i) = 0.5*I*pi;
 	 }
 	 Complex w1 = (-1.0/2.0)*(T-U);
 	 Complex w2 = (-1.0/2.0)*(T+U);
 	 Dw(0,0) += w1;
 	 Dw(3,3) += w1;
 	 Dw(1,1) += -1.0*w1;
 	 Dw(2,2) += -1.0*w1;
 	 Dw(4,4) += w2; 
 	 Dw(5,5) += w2;
          
 	 Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
 	 Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
 	 Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
 	 Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U 
 	   - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z2;
 	 Dz(2,2) = z3;
 	 Dz(3,3) = z4;
 	 Dz(4,4) = Dz(5,5) = z5;
          
 	 G2 = Gamma2(U,T);
        }
      }
      break;
    case LLEE:
      {
        unsigned int numGauge = 1, numBrokenGauge = 2;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(1,0) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Lnu;
 	 double g12 = g_Le;
 	 //double g21 = g_Re;
 	 double g22 = g_Re;
 	 
 	 Complex w1 = 0.25*I*pi;
 	 Dw(0,0) = Dw(1,1) = w1;
          
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z2;
          
 	 G2(0,0) = Gamma2Singlet()(0,0);
        }
      }
      break;
    case EEEE:
    case EEEEiden:
      {
        unsigned int numGauge = 1, numBrokenGauge = 1;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = 1.0;
        if (oneLoop) {
 	 double g11 = g_Re;
 	 //double g12 = g_Re;
 	 //double g21 = g_Re;
 	 double g22 = g_Re;
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
 	 Dz(0,0) = z1;
        }
      }
      break;
    case QQWW:
    case LLWW:
      {
        unsigned int numGauge = 5, numBrokenGauge = 20;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = 1.0; R0(0,1) = 0.5;
        R0(1,0) = 1.0; R0(1,1) = -0.5;
        R0(2,0) = cos2; R0(2,2) = -0.5*sin*cos; R0(2,3) = -0.5*sin*cos; R0(2,4) = sin2;
        R0(3,0) = sin*cos; R0(3,2) = 0.5*cos2; R0(3,3) = -0.5*sin2; R0(3,4) = -sin*cos;
        R0(4,0) = sin*cos; R0(4,2) = -0.5*sin2; R0(4,3) = 0.5*cos2; R0(4,4) = -sin*cos;
        R0(5,0) = sin2; R0(5,2) = 0.5*sin*cos; R0(5,3) = 0.5*sin*cos; R0(5,4) = cos2;
        R0(6,0) = 1.0; R0(6,1) = -0.5;
        R0(7,0) = 1.0; R0(7,1) = 0.5;
        R0(8,0) = cos2; R0(8,2) = 0.5*sin*cos; R0(8,3) = 0.5*sin*cos; R0(8,4) = sin2;
        R0(9,0) = sin*cos; R0(9,2) = -0.5*cos2; R0(9,3) = 0.5*sin2; R0(9,4) = -sin*cos;
        R0(10,0) = sin*cos; R0(10,2) = 0.5*sin2; R0(10,3) = -0.5*cos2; R0(10,4) = -sin*cos;
        R0(11,0) = sin2; R0(11,2) = -0.5*sin*cos; R0(11,3) = -0.5*sin*cos; R0(11,4) = cos2;
        R0(12,1) = -cos/sqrt(2.0); R0(12,3) = -sin/sqrt(2.0);
        R0(13,1) = -sin/sqrt(2.0); R0(13,3) = cos/sqrt(2.0);
        R0(14,1) = cos/sqrt(2.0); R0(14,2) = -sin/sqrt(2.0);
        R0(15,1) = sin/sqrt(2.0); R0(15,2) = cos/sqrt(2.0);
        R0(16,1) = -cos/sqrt(2.0); R0(16,2) = -sin/sqrt(2.0);
        R0(17,1) = -sin/sqrt(2.0); R0(17,2) = cos/sqrt(2.0);
        R0(18,1) = cos/sqrt(2.0); R0(18,3) = -sin/sqrt(2.0);
        R0(19,1) = sin/sqrt(2.0); R0(19,3) = cos/sqrt(2.0);
        if (oneLoop) {
 	 double gW = g_W;
 	 double g1(0.),g2(0.);
 	 if (process==QQWW) {
 	   g1 = g_Lu;
 	   g2 = g_Ld;
 	 }
 	 else if (process==LLWW) {
 	   g1 = g_Lnu;
 	   g2 = g_Le;
 	 }
             
 	 Complex w1 = T-U+5.0/4.0*I*pi;
 	 Complex w2 = -T+U+5.0/4.0*I*pi;
 	 Complex w3 = -0.5*(T+U) + 3.0/4.0*I*pi;
 	 Complex w4 = 0.25*I*pi;
 	 Dw(0,0) = Dw(7,7) = w1;
 	 Dw(1,1) = Dw(6,6) = w2;
 	 for (unsigned int i=12; i<20; i++) {
 	   Dw(i,i) = w3;
 	 }
 	 Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w4;
 	 Dw(8,8) = Dw(9,9) = Dw(10,10) = Dw(11,11) = w4;
          
 	 Complex z1 = 2.0*g1*gW*(U-T) - I*pi*(g1*g1+gW*gW);
 	 Complex z2 = 2.0*g1*gW*(T-U) - I*pi*(g1*g1+gW*gW);
 	 Complex z3 = 2.0*g2*gW*(U-T) - I*pi*(g2*g2+gW*gW);
 	 Complex z4 = 2.0*g2*gW*(T-U) - I*pi*(g2*g2+gW*gW);
 	 Complex z5 = -(g2*gW)*T + (g1*gW)*U - I*pi*(g1*g2+g1*gW-g2*gW);
 	 Complex z6 = (g1*gW)*T - (g2*gW)*U - I*pi*(g1*g2+g1*gW-g2*gW);
 	 Complex z7 = -I*pi*g1*g1;
 	 Complex z8 = -I*pi*g2*g2;
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z2;
 	 Dz(2,2) = Dz(3,3) = Dz(4,4) = Dz(5,5) = z7;
 	 Dz(6,6) = z3;
 	 Dz(7,7) = z4;
 	 Dz(8,8) = Dz(9,9) = Dz(10,10) = Dz(11,11) = z8;
 	 Dz(12,12) = Dz(13,13) = Dz(16,16) = Dz(17,17) = z5;
 	 Dz(14,14) = Dz(15,15) = Dz(18,18) = Dz(19,19) = z6;
          
 	 G2 = Gamma2w(U,T);
        }
      }
      break;
    case QQPhiPhi:
    case LLPhiPhi:
      {
        unsigned int numGauge = 2, numBrokenGauge = 14;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = 0.25;		R0(0,1) = 1.0;
        R0(1,0) = -1.0/8.0;	R0(1,1) = 0.5;
        R0(2,0) = I/8.0;		R0(2,1) = -I/2.0;
        R0(3,0) = -I/8.0;    R0(3,1) = I/2.0;
        R0(4,0) = -1.0/8.0;	R0(4,1) = 1.0/2.0;
        R0(5,0) = -1.0/4.0;	R0(5,1) = 1.0;
        R0(6,0) = 1.0/8.0;	R0(6,1) = 1.0/2.0;
        R0(7,0) = -I/8.0;    R0(7,1) = -I/2.0;
        R0(8,0) = I/8.0;		R0(8,1) = I/2.0;
        R0(9,0) = 1.0/8.0; 	R0(9,1) = 1.0/2.0;
        R0(10,0) = -1.0/(2.0*sqrt(2.0));
        R0(11,0) = I/(2.0*sqrt(2.0));
        R0(12,0) = -1.0/(2.0*sqrt(2.0));
        R0(13,0) = -I/(2.0*sqrt(2.0));
        
        if (oneLoop) {
 	 double g1(0.),g2(0.);
-	 if (process==QQWW) {
+	 if (process==QQPhiPhi) {
 	   g1 = g_Lu;
 	   g2 = g_Ld;
 	 }
-	 else if (process==LLWW) {
+	 else if (process==LLPhiPhi) {
 	   g1 = g_Lnu;
 	   g2 = g_Le;
 	 }
 	 else
 	   assert(false);
 	 double g3 = g_phiPlus;
          
 	 Complex w0 = 0.25*I*pi;
 	 Complex w1 = 0.5*(T-U) + 0.5*I*pi;
 	 Complex w2 = -0.5*(T-U) + 0.5*I*pi;
 	 Complex w3 = 0.25*I*(T-U);
 	 Complex w4 = -0.25*(T+U) + 0.25*I*pi;
 	 Dw(0,0) = w2;
 	 Dw(1,1) = w0; Dw(1,2) = w3; Dw(1,3) = -w3; Dw(1,4) = w0;
 	 Dw(2,1) = -w3; Dw(2,2) = w0; Dw(2,3) = -w0; Dw(2,4) = -w3;
 	 Dw(3,1) = w3; Dw(3,2) = -w0; Dw(3,3) = w0; Dw(3,4) = w3;
 	 Dw(4,1) = w0; Dw(4,2) = w3; Dw(4,3) = -w3; Dw(4,4) = w0;
 	 Dw(5,5) = w1;
 	 Dw(6,6) = w0; Dw(6,7) = -w3; Dw(6,8) = w3; Dw(6,9) = w0;
 	 Dw(7,6) = w3; Dw(7,7) = w0; Dw(7,8) = -w0; Dw(7,9) = w3;
 	 Dw(8,6) = -w3; Dw(8,7) = -w0; Dw(8,8) = w0; Dw(8,9) = -w3;
 	 Dw(9,6) = w0; Dw(9,7) = -w3; Dw(9,8) = w3; Dw(9,9) = w0;
 	 Dw(10,10) = w4; Dw(10,11) = I*w4;
 	 Dw(11,10) = -I*w4; Dw(11,11) = w4;
 	 Dw(12,12) = w4; Dw(12,13) = -I*w4;
 	 Dw(13,12) = I*w4; Dw(13,13) = w4;
          
 	 Complex z1 = 2.0*g3*g1*(T-U) - I*pi*(g3*g3+g1*g1);
 	 Complex z2 = 2.0*g3*g2*(T-U) - I*pi*(g3*g3+g2*g2);
 	 Complex z3 = -I*pi*g1*g1;
 	 Complex z4 = 0.5*I*g1*(T-U);
 	 Complex z5 = 0.25*I*pi;
 	 Complex z6 = -I*pi*g2*g2;
 	 Complex z7 = 0.5*I*g2*(T-U);
 	 Complex z8 = g3*g1*T-g3*g2*U-I*pi*(g1*g2-g2*g3+g1*g3);
 	 Complex z9 = 0.5*I*g2*T-0.5*I*g1*U+pi/2.0*g2-pi/2.0*g1+pi/2.0*g3;
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z3; Dz(1,2) = -z4; Dz(1,3) = z4; Dz(1,4) = -z5;
 	 Dz(2,1) = z4; Dz(2,2) = z3; Dz(2,3) = z5; Dz(2,4) = z4;
 	 Dz(3,1) = -z4; Dz(3,2) = z5; Dz(3,3) = z3; Dz(3,4) = -z4;
 	 Dz(4,1) = -z5; Dz(4,2) = -z4; Dz(4,3) = z4; Dz(4,4) = z3;
 	 Dz(5,5) = z2;
 	 Dz(6,6) = z6; Dz(6,7) = -z7; Dz(6,8) = z7; Dz(6,9) = -z5;
 	 Dz(7,6) = z7; Dz(7,7) = z6; Dz(7,8) = z5; Dz(7,9) = z7;
 	 Dz(8,6) = -z7; Dz(8,7) = z5; Dz(8,8) = z6; Dz(8,9) = -z7;
 	 Dz(9,6) = -z5; Dz(9,7) = -z7; Dz(9,8) = z7; Dz(9,9) = z6;
 	 Dz(10,10) = z8; Dz(10,11) = -z9;
 	 Dz(11,10) = z9; Dz(11,11) = z8;
 	 Dz(12,12) = z8; Dz(12,13) = z9;
 	 Dz(13,12) = -z9; Dz(13,13) = z8;
 
 	 G2 = Gamma2(U,T);
        }
      }
      break;
    case QQWG:
      {
        unsigned int numGauge = 1, numBrokenGauge = 6;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = 1.0/sqrt(2);
        R0(1,0) = 1.0/sqrt(2);
        R0(2,0) = cos/2.0;
        R0(3,0) = sin/2.0;
        R0(4,0) = -cos/2.0;
        R0(5,0) = -sin/2.0;
        if (oneLoop) {
 	 double g1 = g_Lu;
 	 double g2 = g_Ld;
 	 double gW = g_W;
 	 
 	 Complex w1 = -0.5*(T+U) + 0.75*I*pi;
 	 Complex w2 = 0.25*I*pi;
 	 Dw(0,0) = Dw(1,1) = w1;
 	 Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;
          
 	 Complex z1 = gW*g1*T - gW*g2*U - I*pi*(g1*g2+g1*gW-g2*gW);
 	 Complex z2 = gW*g1*U - gW*g2*T - I*pi*(g2*g1+g1*gW-g2*gW);
 	 Complex z3 = -I*pi*g1*g1;
 	 Complex z4 = -I*pi*g2*g2;
 	 Dz(0,0) = z1;
 	 Dz(1,1) = z2;
 	 Dz(2,2) = z3;
 	 Dz(3,3) = z3;
 	 Dz(4,4) = z4;
 	 Dz(5,5) = z4;
          
 	 G2(0,0) = -7.0/4.0*I*pi + (U+T);
        }
      }
      break;
     case QQBG:
       {
 	unsigned int numGauge = 1, numBrokenGauge = 4;
 	R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
 	G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
 	Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
 	Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
 	R0(0,0) = -sin;
 	R0(1,0) = cos;
 	R0(2,0) = -sin;
 	R0(3,0) = cos;
 	if (oneLoop) {
 	  double g1 = g_Lu;
 	  double g2 = g_Ld;
 	  Complex w2 = 0.25*I*pi;
 	  Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = w2;
 	  Complex z3 = -I*pi*g1*g1;
 	  Complex z4 = -I*pi*g2*g2;
 	  Dz(0,0) = z3;
 	  Dz(1,1) = z3;
 	  Dz(2,2) = z4;
 	  Dz(3,3) = z4;
 	  G2(0,0) = Gamma2Singlet()(0,0);
 	}
       }
       break;
    case QQGG:
    case QtQtGG:
      {
        unsigned int numGauge = 3, numBrokenGauge = 6;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge); 
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = R0(3,0) = 1.0;
        R0(1,1) = R0(4,1) = 1.0;
        R0(2,2) = R0(5,2) = 1.0;
        if (oneLoop) {
 	 double g1 = g_Lu;
 	 double g2 = g_Ld;
 	 Complex w2 = 0.25*I*pi;
 	 Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;            
 	 Complex z3 = -I*pi*g1*g1;
 	 Complex z4 = -I*pi*g2*g2;
 	 Dz(0,0) = Dz(1,1) = Dz(2,2) = z3;
 	 Dz(3,3) = Dz(4,4) = Dz(5,5) = z4;
 	 G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
 	 G2 *= 0.0;
 	 G2(0,0) = G2(1,1) = G2(2,2) = Gamma2Singlet()(0,0);
        }
      }
      break;
    case UUBB:
    case DDBB:
    case EEBB:
      {
        unsigned int numGauge = 1, numBrokenGauge = 4;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = sin2;
        R0(1,0) = -sin*cos;
        R0(2,0) = -sin*cos;
        R0(3,0) = cos2;
        if (oneLoop) {
 	 double g1(0.);
 	 if (process==UUBB) {
 	   g1 = g_Ru;
 	 }
 	 else if (process==DDBB) {
 	   g1 = g_Rd;
 	 }
 	 else if (process==EEBB) {
 	   g1 = g_Re;
 	 }
 	 else
 	   assert(false);
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = -I*pi*g1*g1;
 	 Dz(0,0) = Dz(1,1) = Dz(2,2) = Dz(3,3) = z1;
        }
      }
      break;
    case UUPhiPhi:
    case DDPhiPhi:
    case EEPhiPhi:
      {
        unsigned int numGauge = 1, numBrokenGauge = 5;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = 1.0;
         R0(1,0) = 0.5;
         R0(2,0) = -0.5*I;
         R0(3,0) = 0.5*I;
         R0(4,0) = 0.5;
         if (oneLoop) {
 	  double g1(0.);
             if (process==UUPhiPhi) {
                g1 = g_Ru;
             }
             else if (process==DDPhiPhi) {
                g1 = g_Rd;
             }
             else if (process==EEPhiPhi) {
                g1 = g_Re;
             }
 	    double g3 = g_phiPlus;
             Dw(0,0) = Dw(1,4) = Dw(4,1) = 0.25*I*pi;
             Dw(2,3) = Dw(3,2) = -0.25*I*pi;
 	    Complex z1 = 2.0*g3*g1*(T-U) - I*pi*(g3*g3+g1*g1);
 	    Complex z2 = 0.5*I*g1*g1;
 	    Complex z3 = -I*pi*g1*g1;
 	    Complex z4 = 0.25*I*pi;
             Dz(0,0) = z1;
             Dz(1,1) = z3;		Dz(1,2) = -z2;		Dz(1,3) = z2;		Dz(1,4) = -z4;
             Dz(2,1) = z2;		Dz(2,2) = z3;		Dz(2,3) = z4;		Dz(2,4) = z2;
             Dz(3,1) = -z2;		Dz(3,2) = z4;		Dz(3,3) = z3;		Dz(3,4) = -z2;
             Dz(4,1) = -z4;		Dz(4,2) = -z2;		Dz(4,3) = z2;		Dz(4,4) = z3;
             G2(0,0) = Gamma2Singlet()(0,0);
 	}
      }
      break;
    case UUBG:
    case DDBG:
      {
        unsigned int numGauge = 1, numBrokenGauge = 2;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        R0(0,0) = -sin;
        R0(1,0) = cos;
        if (oneLoop) {
 	 double g1(0.);
 	   if (process==UUBG) {
 	     g1 = g_Ru;
 	   }
 	   else if (process==DDBG) {
 	     g1 = g_Rd;
 	   }
 	   else
 	     assert(false);
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = -I*pi*g1*g1;
 	 Dz(0,0) = Dz(1,1) = z1;
        }
      }
      break;
    case UUGG:
    case tRtRGG:
    case DDGG:
      {
        unsigned int numGauge = 3, numBrokenGauge = 3;
        R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
        G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
        Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
        Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+       R0(0,0) = R0(1,1) = R0(2,2) = 1.0;
        if (oneLoop) {
 	 double g1(0.);
 	 if ((process==UUGG)||(process==tRtRGG)) {
 	   g1 = g_Ru;
 	 }
 	 else if (process==DDGG) {
 	   g1 = g_Rd;
 	 }
 	 else
 	   assert(false);
 	 // There is no Dw contribution for two SU(2) singlets.
 	 Complex z1 = -I*pi*g1*g1;
 	 Dz(0,0) = Dz(1,1) = Dz(2,2) = z1;
        }
      }
      break;
    default:
      assert(false);
    }
    
    double aW = ElectroWeakReweighter::coupling()->aW(mu);
    double aZ = ElectroWeakReweighter::coupling()->aZ(mu);
    Energy mZ = ElectroWeakReweighter::coupling()->mZ();
    Energy mW = ElectroWeakReweighter::coupling()->mW();
    
    if (!oneLoop) {
       return R0;
-   } 
+   }
+   boost::numeric::ublas::matrix<Complex> output(R0);
    boost::numeric::ublas::matrix<Complex> temp(R0.size1(),R0.size2());
    boost::numeric::ublas::axpy_prod(R0,G2,temp);
-   R0+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
+   output+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
    boost::numeric::ublas::axpy_prod(Dw,R0,temp);
-   R0+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
+   output+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
    boost::numeric::ublas::axpy_prod(Dz,R0,temp);
-   R0+=aZ/(4.0*pi)*4.0*log(mZ/mu)*temp;
-   return R0;
+   output+=aZ/(4.0*pi)*4.0*log(mZ/mu)*temp;
+   return output;
 }
diff --git a/MatrixElement/EW/ElectroWeakReweighter.cc b/MatrixElement/EW/ElectroWeakReweighter.cc
--- a/MatrixElement/EW/ElectroWeakReweighter.cc
+++ b/MatrixElement/EW/ElectroWeakReweighter.cc
@@ -1,103 +1,200 @@
 // -*- 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"
+#include "boost/numeric/ublas/matrix.hpp"
+#include "boost/numeric/ublas/operation.hpp"
+#include "EWProcess.h"
+#include "HighEnergyMatching.h"
+#include "ElectroWeakMatching.h"
 
 using namespace Herwig;
 
 tEWCouplingsPtr ElectroWeakReweighter::staticEWCouplings_ = tEWCouplingsPtr();
 
 
 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_ << collinearSudakov_;
+  os << EWCouplings_ << collinearSudakov_ << softSudakov_;
 }
 
 void ElectroWeakReweighter::persistentInput(PersistentIStream & is, int) {
-  is >> EWCouplings_ >> collinearSudakov_;
+  is >> EWCouplings_ >> collinearSudakov_ >> softSudakov_;
 }
 
 
 // The following static variable is needed for the type
 // description system in ThePEG.
 DescribeClass<ElectroWeakReweighter,ReweightBase>
   describeHerwigElectroWeakReweighter("Herwig::ElectroWeakReweighter", "HwMEEW.so");
 
 void ElectroWeakReweighter::Init() {
 
   static ClassDocumentation<ElectroWeakReweighter> documentation
     ("There is no documentation for the ElectroWeakReweighter class");
 
   static Reference<ElectroWeakReweighter,EWCouplings> interfaceEWCouplings
     ("EWCouplings",
      "The object to calculate the electroweak couplings",
      &ElectroWeakReweighter::EWCouplings_, false, false, true, false, false);
 
   static Reference<ElectroWeakReweighter,CollinearSudakov> interfaceCollinearSudakov
     ("CollinearSudakov",
      "The collinear Sudakov",
      &ElectroWeakReweighter::collinearSudakov_, false, false, true, false, false);
 
+  static Reference<ElectroWeakReweighter,SoftSudakov> interfaceSoftSudakov
+    ("SoftSudakov",
+     "The soft Sudakov",
+     &ElectroWeakReweighter::softSudakov_, false, false, true, false, false);
+
 }
 
+namespace {
+
+void axpy_prod_local(const boost::numeric::ublas::matrix<Complex>  & A,
+		     const boost::numeric::ublas::matrix<complex<InvEnergy2> > & B,
+		     boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
+  assert(A.size2()==B.size1());
+  C.resize(A.size1(),B.size2());
+  for(unsigned int ix=0;ix<A.size1();++ix) {
+    for(unsigned int iy=0;iy<B.size2();++iy) {
+      C(ix,iy) = ZERO;
+      for(unsigned int iz=0;iz<A.size2();++iz) {
+	C(ix,iy) += A(ix,iz)*B(iz,iy);
+      }
+    }
+  }
+}
+
+void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> > & A,
+		     const boost::numeric::ublas::matrix<Complex> & B,
+		     boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
+  assert(A.size2()==B.size1());
+  C.resize(A.size1(),B.size2());
+  for(unsigned int ix=0;ix<A.size1();++ix) {
+    for(unsigned int iy=0;iy<B.size2();++iy) {
+      C(ix,iy) = ZERO;
+      for(unsigned int iz=0;iz<A.size2();++iz) {
+	C(ix,iy) += A(ix,iz)*B(iz,iy);
+      }
+    }
+  }
+}
+
+}
 double ElectroWeakReweighter::weight() const {
   EWCouplings_->initialize();
   staticEWCouplings_ = EWCouplings_;
-  cerr << "aEM\n";
-  for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.1) {
-    cerr << scale/GeV << " " 
-  	 << EWCouplings_->aEM(scale) << "\n";
-  }
-  cerr << "aS\n";
-  for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
-    cerr << scale/GeV << " " 
-  	 << EWCouplings_->aS(scale) << "\n";
-  }
-  cerr << "y_t\n";
-  for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
-    cerr << scale/GeV << " " 
-  	 << EWCouplings_->y_t(scale) << "\n";
-  }
-  cerr << "lambda\n";
-  for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
-    cerr << scale/GeV << " " 
-  	 << EWCouplings_->lambda(scale) << "\n";
-  }
-  cerr << "vev\n";
-  for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
-    cerr << scale/GeV << " " 
-  	 << EWCouplings_->vev(scale)/GeV << "\n";
-  }
+  // cerr << "aEM\n";
+  // for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.1) {
+  //   cerr << scale/GeV << " " 
+  // 	 << EWCouplings_->aEM(scale) << "\n";
+  // }
+  // cerr << "aS\n";
+  // for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
+  //   cerr << scale/GeV << " " 
+  // 	 << EWCouplings_->aS(scale) << "\n";
+  // }
+  // cerr << "y_t\n";
+  // for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
+  //   cerr << scale/GeV << " " 
+  // 	 << EWCouplings_->y_t(scale) << "\n";
+  // }
+  // cerr << "lambda\n";
+  // for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
+  //   cerr << scale/GeV << " " 
+  // 	 << EWCouplings_->lambda(scale) << "\n";
+  // }
+  // cerr << "vev\n";
+  // for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
+  //   cerr << scale/GeV << " " 
+  // 	 << EWCouplings_->vev(scale)/GeV << "\n";
+  // }
   collinearSudakov_->makePlots();
+  Energy2 s = sqr(5000.*GeV);
+  Energy2 t = -0.25*s;
+  Energy2 u = -0.75*s;
+  testEvolution(s,t,u);
 
-  // cerr << "testing got here ???\n";
+
+
+
   // cerr <<  subProcess() << "\n";
   // cerr << *subProcess() << "\n";
   // cerr << subProcess()->outgoing()[0] << *subProcess()->outgoing()[0] << "\n";
   // cerr << subProcess()->outgoing()[0]->spinInfo() << "\n";
   // cerr << subProcess()->outgoing()[0]->spinInfo()->productionVertex() << "\n";
   assert(false);
   staticEWCouplings_ = tEWCouplingsPtr();
 }
+
+void ElectroWeakReweighter::testEvolution(Energy2 s,Energy2 t, Energy2 u) const {
+  Energy highScale = sqrt(s);
+  Energy ewScale = coupling()->mZ();
+  Energy lowScale = 50.0*GeV;
+  for (unsigned int i=0; i<45;++i) {
+    EWProcess::Process process = (EWProcess::Process)i;
+    cerr << "process " << process << "\n";
+    // result for all EW and QCD SCET contributions:
+    boost::numeric::ublas::matrix<complex<InvEnergy2> > highMatch_val 
+      = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,true,true);
+    boost::numeric::ublas::matrix<Complex> highRunning_val
+      = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process);
+    boost::numeric::ublas::matrix<Complex> ewMatch_val = 
+      ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true);
+    boost::numeric::ublas::matrix<Complex> lowRunning_val = 
+      softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process);
+    boost::numeric::ublas::matrix<Complex> collinearHighRunning_val =
+      collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false);
+    boost::numeric::ublas::matrix<Complex> collinearEWMatch_val =
+      collinearSudakov_->electroWeakMatching(ewScale,s,process,true);
+    boost::numeric::ublas::matrix<Complex> collinearLowRunning_val =
+      collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process);
+    boost::numeric::ublas::matrix<Complex> lowMatchTemp_val = 
+      boost::numeric::ublas::zero_matrix<Complex>(ewMatch_val.size1(),ewMatch_val.size2());
+    for (unsigned int ii=0; ii<ewMatch_val.size1(); ++ii) {
+      for (unsigned int jj=0; jj<ewMatch_val.size2(); ++jj) {
+	lowMatchTemp_val(ii,jj) = collinearEWMatch_val(ii,jj)*ewMatch_val(ii,jj);
+      }
+    }
+    boost::numeric::ublas::matrix<Complex> temp(highRunning_val.size1(),collinearHighRunning_val.size2());
+    boost::numeric::ublas::axpy_prod(highRunning_val,collinearHighRunning_val,temp);
+    boost::numeric::ublas::matrix<Complex> temp2(collinearLowRunning_val.size1(),lowRunning_val.size2());
+    boost::numeric::ublas::axpy_prod(collinearLowRunning_val,lowRunning_val,temp2);
+    boost::numeric::ublas::matrix<Complex> temp3(temp2.size1(),lowMatchTemp_val.size2());
+    boost::numeric::ublas::axpy_prod(temp2,lowMatchTemp_val,temp3);
+    temp2.resize(temp3.size1(),temp.size2());
+    boost::numeric::ublas::axpy_prod(temp3,temp,temp2);
+    boost::numeric::ublas::matrix<complex<InvEnergy2> >  result(temp2.size1(),highMatch_val.size2());
+    axpy_prod_local(temp2,highMatch_val,result);
+    for(unsigned int ix=0;ix<result.size1();++ix) {
+      for(unsigned int iy=0;iy<result.size2();++iy) {
+	cerr << s*result(ix,iy) << " ";
+      }
+      cerr << "\n";
+    }
+  }
+}
diff --git a/MatrixElement/EW/ElectroWeakReweighter.h b/MatrixElement/EW/ElectroWeakReweighter.h
--- a/MatrixElement/EW/ElectroWeakReweighter.h
+++ b/MatrixElement/EW/ElectroWeakReweighter.h
@@ -1,127 +1,140 @@
 // -*- 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"
 #include "CollinearSudakov.h"
+#include "SoftSudakov.h"
 
 namespace Herwig {
 
 using namespace ThePEG;
 
 /**
  * The ElectroWeakReweighter class.
  *
  * @see \ref ElectroWeakReweighterInterfaces "The interfaces"
  * defined for ElectroWeakReweighter.
  */
 class ElectroWeakReweighter: public ReweightBase {
 
 public:
 
   /** @name Standard constructors and destructors. */
   //@{
   /**
    * The default constructor.
    */
   ElectroWeakReweighter();
 
   /**
    * The destructor.
    */
   virtual ~ElectroWeakReweighter();
   //@}  
 
 public:
 
   /**
    * Return the weight for the kinematical configuation provided by
    * the assigned XComb object (in the LastXCombInfo base class).
    */
   virtual double weight() const;
 
   /**
    *
    */
   static tEWCouplingsPtr coupling() {
     assert(staticEWCouplings_);
     return staticEWCouplings_;
   }
 
 public:
 
   /** @name Functions used by the persistent I/O system. */
   //@{
   /**
    * Function used to write out object persistently.
    * @param os the persistent output stream written to.
    */
   void persistentOutput(PersistentOStream & os) const;
 
   /**
    * Function used to read in object persistently.
    * @param is the persistent input stream read from.
    * @param version the version number of the object when written.
    */
   void persistentInput(PersistentIStream & is, int version);
   //@}
 
   /**
    * The standard Init function used to initialize the interfaces.
    * Called exactly once for each class by the class description system
    * before the main function starts or
    * when this class is dynamically loaded.
    */
   static void Init();
 
 protected:
 
+  /**
+   *  Check the evolution for a fixed s,t,u
+   */
+  void testEvolution(Energy2 s,Energy2 t, Energy2 u) const;
+
+protected:
+
   /** @name Clone Methods. */
   //@{
   /**
    * Make a simple clone of this object.
    * @return a pointer to the new object.
    */
   virtual IBPtr clone() const;
 
   /** Make a clone of this object, possibly modifying the cloned object
    * to make it sane.
    * @return a pointer to the new object.
    */
   virtual IBPtr fullclone() const;
   //@}
 
 private:
 
   /**
    * The assignment operator is private and must never be called.
    * In fact, it should not even be implemented.
    */
   ElectroWeakReweighter & operator=(const ElectroWeakReweighter &);
 
 private:
 
   /**
    *  The Electroweak Couplings
    */
   EWCouplingsPtr EWCouplings_;
 
   /**
    *  The Collinear Sudakov
    */
   CollinearSudakovPtr collinearSudakov_;
 
   /**
+   *  The Soft Sudakov
+   */
+  SoftSudakovPtr softSudakov_;
+
+  /**
    *  The couplings to allow global access
    */
   static tEWCouplingsPtr staticEWCouplings_;
 
 };
 
 }
 
 #endif /* Herwig_ElectroWeakReweighter_H */
diff --git a/MatrixElement/EW/Makefile.am b/MatrixElement/EW/Makefile.am
--- a/MatrixElement/EW/Makefile.am
+++ b/MatrixElement/EW/Makefile.am
@@ -1,9 +1,10 @@
 pkglib_LTLIBRARIES = HwMEEW.la
 HwMEEW_la_SOURCES = EWProcess.h GroupInvariants.h GroupInvariants.cc \
 ElectroWeakReweigter.h ElectroWeakReweighter.cc \
 SoftSudakov.h SoftSudakov.cc \
 CollinearSudakov.h CollinearSudakov.cc \
 HighEnergyMatching.h HighEnergyMatching.cc \
 ElectroWeakMatching.h ElectroWeakMatching.cc \
-EWCouplings.h EWCouplings.fh EWCouplings.cc
+EWCouplings.h EWCouplings.fh EWCouplings.cc \
+expm-1.hpp
 HwMEEW_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
diff --git a/MatrixElement/EW/SoftSudakov.cc b/MatrixElement/EW/SoftSudakov.cc
--- a/MatrixElement/EW/SoftSudakov.cc
+++ b/MatrixElement/EW/SoftSudakov.cc
@@ -1,881 +1,885 @@
 // -*- 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"
+#include "expm-1.h"
 
 using namespace Herwig;
 using namespace GroupInvariants;
 
-SoftSudakov::SoftSudakov() : K_ORDER_(3) {}
+SoftSudakov::SoftSudakov() : K_ORDER_(3), integrator_(0.,1e-5,1000) {}
 
 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));
+      if(G3_(row_,col_) == 0. && G2_(row_,col_) == 0. && G1_(row_,col_) == 0.) {
+	gamma(row_,col_) = 0.;
+      }
+      else {
+	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;
+  return boost::numeric::ublas::expm_pad(gamma,7);
 }
 
 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);
   }
-   
+  // return the answer
   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);
+      G3 = Gamma3g(U,T);
       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/expm-1.h b/MatrixElement/EW/expm-1.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/expm-1.h
@@ -0,0 +1,147 @@
+//
+//  Copyright (c) 2007
+//  Tsai, Dung-Bang	
+//  National Taiwan University, Department of Physics
+// 
+//  E-Mail : dbtsai (at) gmail.com
+//  Begine : 2007/11/20
+//  Last modify : 2007/11/22
+//  Version : v0.1
+//
+//  EXPGM_PAD computes the matrix exponential exp(H) for general matrixs,
+//  including complex and real matrixs using the irreducible (p,p) degree
+//  rational Pade approximation to the exponential 
+//  exp(z) = r(z)=(+/-)( I+2*(Q(z)/P(z))).
+//
+//  Usage : 
+//
+//    U = expm_pad(H)
+//    U = expm_pad(H, p)
+//    
+//    where p is internally set to 6 (recommended and gererally satisfactory).
+//
+//  See also MATLAB supplied functions, EXPM and EXPM1.
+//
+//  Reference :
+//  EXPOKIT, Software Package for Computing Matrix Exponentials.
+//  ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
+//
+//  Permission to use, copy, modify, distribute and sell this software
+//  and its documentation for any purpose is hereby granted without fee,
+//  provided that the above copyright notice appear in all copies and
+//  that both that copyright notice and this permission notice appear
+//  in supporting documentation.  The authors make no representations
+//  about the suitability of this software for any purpose.
+//  It is provided "as is" without express or implied warranty.
+//
+
+#ifndef _BOOST_UBLAS_EXPM_
+#define _BOOST_UBLAS_EXPM_
+#include <complex>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+
+namespace boost { namespace numeric { namespace ublas {
+
+template<typename MATRIX> MATRIX expm_pad(const MATRIX &H, const int p = 6) {
+  typedef typename MATRIX::value_type value_type;
+  typedef typename MATRIX::size_type size_type;
+  typedef double real_value_type;	// Correct me. Need to modify.
+  assert(H.size1() == H.size2());	
+  const size_type n = H.size1();
+  const identity_matrix<value_type> I(n);
+  matrix<value_type> U(n,n),H2(n,n),P(n,n),Q(n,n);
+  real_value_type norm = 0.0;
+
+  // Calcuate Pade coefficients  (1-based instead of 0-based as in the c vector)
+  vector<real_value_type> c(p+2);
+  c(1)=1;  
+  for(size_type i = 1; i <= p; ++i) 
+    c(i+1) = c(i) * ((p + 1.0 - i)/(i * (2.0 * p + 1 - i)));
+  // Calcuate the infinty norm of H, which is defined as the largest row sum of a matrix
+  for(size_type i=0; i<n; ++i)
+    {
+      real_value_type temp = 0.0;
+      for(size_type j=0;j<n;j++)
+	temp += std::abs<real_value_type>(H(j,i)); // Correct me, if H is complex, can I use that abs?
+      norm = std::max<real_value_type>(norm, temp);
+    }
+  if (norm == 0.0) 
+    {
+      boost::throw_exception(boost::numeric::ublas::bad_argument());
+      std::cerr<<"Error! Null input in the routine EXPM_PAD.\n";
+      exit(0);
+    }
+  // Scaling, seek s such that || H*2^(-s) || < 1/2, and set scale = 2^(-s)
+  int s = 0;
+  real_value_type scale = 1.0;
+  if(norm > 0.5) {
+    s = std::max<int>(0, static_cast<int>((log(norm) / log(2.0) + 2.0)));
+    scale /= static_cast<real_value_type>(std::pow(2.0, s));
+    U.assign(scale * H); // Here U is used as temp value due to that H is const
+  }
+  else
+    U.assign(H);
+  // Horner evaluation of the irreducible fraction, see the following ref above.
+  // Initialise P (numerator) and Q (denominator) 
+  H2.assign( prod(U, U) );
+  Q.assign( c(p+1)*I );
+  P.assign( c(p)*I );
+  size_type odd = 1;
+  for( size_type k = p - 1; k > 0; --k)
+    {
+      if( odd == 1)
+	{
+	  Q = ( prod(Q, H2) + c(k) * I ); 
+	}
+      else
+	{
+	  P = ( prod(P, H2) + c(k) * I );
+	}
+      odd = 1 - odd;
+    }
+  if( odd == 1)
+    {
+      Q = ( prod(Q, U) );	
+      Q -= P ;
+      //U.assign( -(I + 2*(Q\P)));
+    }
+  else
+    {
+      P = (prod(P, U));
+      Q -= P;
+      //U.assign( I + 2*(Q\P));
+    }
+  // In origine expokit package, they use lapack ZGESV to obtain inverse matrix,
+  // and in that ZGESV routine, it uses LU decomposition for obtaing inverse matrix.
+  // Since in ublas, there is no matrix inversion template, I simply use the build-in
+  // LU decompostion package in ublas, and back substitute by myself.
+  //
+  //////////////// Implement Matrix Inversion ///////////////////////
+  permutation_matrix<size_type> pm(n); 
+  int res = lu_factorize(Q, pm);
+  if( res != 0)
+    {
+      std::cerr << "Error in the matrix inversion in template expm_pad.\n";
+      exit(0);
+    }
+  H2 = I;  // H2 is not needed anymore, so it is temporary used as identity matrix for substituting.
+  
+  lu_substitute(Q, pm, H2); 
+  if( odd == 1) 
+    U.assign( -(I + 2.0 * prod(H2, P)));
+  else
+    U.assign( I + 2.0 * prod(H2, P));
+  // Squaring 
+  for(size_t i = 0; i < s; ++i)
+    {
+      U = (prod(U,U));
+    }
+  return U;
+ }
+ 
+}}}
+
+
+#endif