No OneTemporary
Actions

Size

83 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/EW/CollinearSudakov.cc b/MatrixElement/EW/CollinearSudakov.cc
	--- a/MatrixElement/EW/CollinearSudakov.cc
	+++ b/MatrixElement/EW/CollinearSudakov.cc
	@@ -1,1192 +1,1730 @@
	// -- 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.SU30.5plog + 0.5nonCusp.SU3;
	else
	gamma += (C_A(3))cusp.SU30.5plog + 0.5nonCusp.SU3;
	}
	// SU(2)
	if(SU2_) {
	if (fermion_ \|\| longitudinal_ )
	gamma += (C_F(2))cusp.SU20.5plog + 0.5nonCusp.SU2;
	else
	gamma += (C_A(2))cusp.SU20.5plog + 0.5nonCusp.SU2;
	}

	if (fermion_ \|\| longitudinal_ ) {
	gamma += sqr(Y_)(cusp.U10.5plog + 0.5nonCusp.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.0sqr(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.SU30.5plog + 0.5nonCusp.SU3;
	else
	gamma += C_F(3)cusp.SU30.5blog + 0.5nonCusp.SU3;
	}
	else {
	gamma += C_A(3)cusp.SU30.5plog + 0.5nonCusp.SU3;
	}
	}
	// fermions
	if (fermion_) {
	// not a bHQET top quark field
	if (abs(boostFactor_)<0.001)
	gamma += sqr(Q_)(cusp.U10.5plog + 0.5nonCusp.U1);
	else
	gamma += sqr(Q_)(cusp.U10.5blog + 0.5nonCusp.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.U10.5blog + 0.5nonCusp.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.0lw0.5PlusLog(s/EWScale/EWScale)-2.0lwlw-2.0lw-5.0pipi/12.0+1.0;
	complex<double> F_Z = 4.0lz0.5PlusLog(s/EWScale/EWScale)-2.0lzlz-2.0lz-5.0pipi/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.5WFC.aW0-0.5WFC.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.5WFC.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.5WFC.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.5y_ty_t/(16.0pipi)(0.25-0.5log(mH/EWScale)-0.5*lz+
	0.5WFC.atHH+0.5WFC.atZZ+WFC.ctHH+WFC.ctZZ+
	WFC.ctW0-WFC.btHH+WFC.btZZ);
	complex<double> HtR = HtL-0.5y_ty_t/(16.0pipi)*(0.25-lw+WFC.atW0);
	complex<double> HbL = -0.5y_ty_t/(16.0pipi)*(0.25-lw+WFC.at0W);

	complex<double> F_tL = (4.0/3.0aS+4.0/9.0aEM)/(4.0pi)(2.0log(mT/EWScale)log(mT/EWScale)-log(mT/EWScale)+
	pi*pi/12.0+2.0) + HtL +
	aW/(4.0pi)0.5W1 + aZ/(4.0pi)*U1;

	complex<double> F_tR = (4.0/3.0aS+4.0/9.0aEM)/(4.0pi)(2.0log(mT/EWScale)log(mT/EWScale)-log(mT/EWScale)+
	pi*pi/12.0+2.0) + HtR -
	aW/(4.0pi)0.25WFC.cW0 + aZ/(4.0pi)*U2;

	complex<double> F_bL = HbL + aW/(4.0pi)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.0pi)cW2*(F_Z+WFC.fsWZWZ) +
	aW/(4.0pi)cW2(F_W+WFC.fs1ZW) + aW/(4.0pi)sW2(F_W+WFC.fs10) +
	aW/(4.0pi)sW2(2.0lw*lw-
	2.0lw+pipi/12.0+2.0) + 0.5*WFC.RW;
	complex<double> D_C_WZ = aW/(4.0pi)2.0(F_W+WFC.fsZW1) + 0.5WFC.RZ + sqrt(sW2/cW2)*WFC.RAtoZ;
	complex<double> D_C_WA = aW/(4.0pi)2.0(F_W+WFC.fs01) + 0.5WFC.RA + sqrt(cW2/sW2)*WFC.RZtoA;
	complex<double> D_C_BZ = 0.5WFC.RZ - sqrt(cW2/sW2)WFC.RAtoZ;
	complex<double> D_C_BA = 0.5WFC.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.0pi)0.25*(F_W+WFC.fs1HW) +
	aW/(4.0pi)0.25(F_W+WFC.fs1ZW) + aZ/(4.0pi)gPHIgPHI*(F_Z+WFC.fsWZWZ) +
	aW/(4.0pi)sW2(2.0lwlw-2.0lw+pi*pi/12.0+2.0) +
	0.5*WFC.RW + WFC.EW;
	complex<double> D_C_PhiZ = aW/(4.0pi)0.5(F_W+WFC.fsZW1) + aZ/(4.0pi)0.25(F_Z+WFC.fs1HZ) + 0.5*WFC.RZ + WFC.EZ;
	complex<double> D_C_PhiH = aW/(4.0pi)0.5(F_W+WFC.fsHW1) + aZ/(4.0pi)0.25(F_Z+WFC.fsHZ1) + 0.5*WFC.RH;

	complex<double> D_C_GG = aS/(4.0pi)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.424777961271568I) + 1.2785863646210789lZ;
	WFC.RA = (1.4835982362022198 + 1.855775680704845pow(10.,-9)I) - 0.27209907467584415*lZ;
	WFC.RZ = (1.5114724841549798 - 9.926944419863688I) + 1.0834802397165764lZ;
	WFC.RAtoZ = (0.3667485032811715 - 2.2770907130064835I) - 1.2994544609942593lZ;
	WFC.RZtoA = -0.2095310079712942 + 0.8320191107808546*lZ;
	WFC.RH = (12.229832449946716 - 1.764310346241984210.0pow(10.,-12)I) + 5.309998583664737lZ;
	WFC.RPhi = (5.569012418081201 + 1.54391335814173560.10pow(10.,-9)I) + 5.309998583664737lZ;
	WFC.RPhi3 = (8.945333042265943 + 5.499309445612249pow(10.,-12)I) + 5.309998583664737*lZ;
	WFC.EW = (3.967645734304811 + 4.712388980384717I) + 2.238332625165702lZ;
	WFC.EZ = (5.916079892937651 + 4.96347220970469I) + 2.1132591719740788lZ;

	WFC.RW = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.RA = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.RZ = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.RAtoZ = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.RZtoA = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.RPhi = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.EW = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.EZ = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.RPhi3 = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);
	WFC.RH = ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0pi);


	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.0pi)
	(4.0lw0.5PlusLog(s/EWScale/EWScale) - 2.0lw*lw -
	3.0lw - 5.0/12.0pi*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.0pi)
	(4.0lz0.5PlusLog(s/EWScale/EWScale) - 2.0lz*lz -
	3.0lz - 5.0/12.0pi*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 + std::log(result(i,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) = colQcolQcolQ*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) = colQtcolQtcolQ*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) = colQcolQcolU*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) = colQtcolQtcolU*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) = colQcolQcoltR*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) = colQcolQcolD*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) = colQtcolQtcolD*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) = colQcolQcolL*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) = colQcolQcolE*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) = colUcolUcolU*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) = coltRcoltRcolU*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) = colUcolUcolD*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) = coltRcoltRcolD*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) = colUcolUcolL*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) = colUcolUcolE*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) = colDcolDcolD*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) = colDcolDcolL*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) = colDcolDcolE*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) = colLcolLcolL*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) = colLcolLcolE*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) = colEcolEcolE*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) = colQcolQcolW*colW;
	+ result(2,2) = result(3,3) = colQcolQcolW*colB;
	+ result(4,4) = colQcolQcolB*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) = colQcolQcolPhi*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) = colQcolQcolW*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) = colQcolQcolB*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) = colQcolQcolG*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) = colQtcolQtcolG*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) = colLcolLcolW*colW;
	+ result(2,2) = result(3,3) = colLcolLcolW*colB;
	+ result(4,4) = colLcolLcolB*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) = colLcolLcolPhi*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) = colUcolUcolB*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) = colUcolUcolPhi*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) = colUcolUcolB*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) = colUcolUcolG*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) = coltRcoltRcolG*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) = colDcolDcolB*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) = colDcolDcolPhi*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) = colDcolDcolB*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) = colDcolDcolG*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) = colEcolEcolB*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) = colEcolEcolPhi*colPhi;
	+ }
	+ break;
	+
	+ default:
	+ assert(false);
	+ }
	+
	+ return result;
	+}
	diff --git a/MatrixElement/EW/CollinearSudakov.h b/MatrixElement/EW/CollinearSudakov.h
	--- a/MatrixElement/EW/CollinearSudakov.h
	+++ b/MatrixElement/EW/CollinearSudakov.h
	@@ -1,527 +1,535 @@
	// -- C++ --
	#ifndef Herwig_CollinearSudakov_H
	#define Herwig_CollinearSudakov_H
	//
	// This is the declaration of the CollinearSudakov class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "Herwig/Utilities/GSLIntegrator.h"
	#include <boost/numeric/ublas/matrix.hpp>
	#include "EWProcess.h"
	#include "CollinearSudakov.fh"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* Struct for the wavefunction corrections
	*/
	struct WaveFunctionCorrections {
	Complex RW;
	Complex RA;
	Complex RZ;
	Complex RAtoZ;
	Complex RZtoA;
	Complex RPhi;
	Complex EW;
	Complex EZ;
	Complex RPhi3;
	Complex RH;
	Complex tLuLDiff;
	Complex bLdLDiff;
	Complex tRuRDiff;

	// The following are constants from parameter integrals:

	Complex fFW0;
	Complex fF0W;
	Complex fFZZ;
	Complex aHH;
	Complex aZZ;
	Complex aW0;
	Complex a0W;
	Complex bHH;
	Complex bZZ;
	Complex cHH;
	Complex cZZ;
	Complex cW0;
	Complex atHH;
	Complex atZZ;
	Complex atW0;
	Complex at0W;
	Complex ctHH;
	Complex ctZZ;
	Complex ctW0;
	Complex btHH;
	Complex btZZ;

	Complex fs10;
	Complex fs1ZW;
	Complex fsWZWZ;
	Complex fsZW1;
	Complex fs01;
	Complex fsHW1;
	Complex fsHZ1;
	Complex fs1HW;
	Complex fs1HZ;
	};

	/**
	* Here is the documentation of the CollinearSudakov class.
	*
	* @see \ref CollinearSudakovInterfaces "The interfaces"
	* defined for CollinearSudakov.
	*/
	class CollinearSudakov: public Interfaced {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	CollinearSudakov();

	/**
	* The destructor.
	*/
	virtual ~CollinearSudakov();
	//@}

	public:

	/**
	- * Evalaute the electroweka matching as a matrix
	+ * Evalaute the electroweak matching as a matrix
	*/
	boost::numeric::ublas::matrix<Complex>
	electroWeakMatching(Energy EWScale, Energy2 s,
	Herwig::EWProcess::Process process,
	bool oneLoop);

	+ /**
	+ * Evalaute the high energy running as a matrix
	+ */
	+ boost::numeric::ublas::matrix<Complex>
	+ highEnergyRunning(Energy highScale, Energy EWScale, Energy2 s,
	+ Herwig::EWProcess::Process process,
	+ bool fixedOrder);
	+
	public:

	/**
	* Make plots for tests
	*/
	void makePlots();

	protected:

	/**
	* Evaluate the high scale contributions
	*/
	void evaluateHighScale(Energy highScale, Energy EWScale, Energy2 S);

	/**
	* Evaluate the low scale contributions
	*/
	void evaluateLowScale(Energy EWScale, Energy lowScale, Energy2 S);

	/**
	* Evaluate the matching
	*/
	void evaluateMatching(Energy EWScale,Energy2 S);

	public:

	/**
	* The operator to be integrated
	*/
	InvEnergy operator ()(Energy mu) const {
	if(high_) return highScaleIntegrand(mu);
	else return lowScaleIntegrand(mu);
	}
	/** Argument type for GaussianIntegrator */
	typedef Energy ArgType;
	/** Return type for GaussianIntegrator */
	typedef InvEnergy ValType;

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/**
	* The integral of the high scale part of the Sudakov
	*/
	Complex highScaleIntegral(bool SU3, bool SU2, double Y,
	Energy2 s, Energy mu_h, Energy mu_l, bool fermion,
	bool longitudinal, double yukFactor);

	/**
	* the integral of the low scale part of the Sudakov
	*/
	Complex lowScaleIntegral(bool SU3, double Q, Energy2 s,
	Energy mu_h, Energy mu_l, bool fermion,
	double boostFactor);

	protected:

	/**
	* High-scale integrand
	*/
	InvEnergy highScaleIntegrand(Energy mu) const;

	/**
	* Low-scale integrand
	*/
	InvEnergy lowScaleIntegrand(Energy mu) const;

	/**
	* Calculate the wavefunction corrections
	*/
	WaveFunctionCorrections waveFunctionCorrections(Energy EWScale);

	/**
	* Collinear matiching for W
	*/
	Complex CollinearDw(Energy2 s, Energy EWScale);

	/**
	* Collinear matching for Z
	*/
	Complex CollinearDz(Energy2 s, Energy EWScale);

	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;
	//@}


	// If needed, insert declarations of virtual function defined in the
	// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).


	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	CollinearSudakov & operator=(const CollinearSudakov &);
	private:

	/**
	* Parameters for the integrand
	*/
	//@{
	/**
	* Whether high or low scale
	*/
	bool high_;

	/**
	* Whether real of imaginary part
	*/
	bool real_;

	/**
	* \f$SU(3)\f$
	*/
	bool SU3_;

	/**
	*
	*/
	bool SU2_;

	/**
	*
	*/
	double Y_;

	/**
	*
	*/
	Energy2 s_;

	/**
	*
	*/
	bool fermion_;

	/**
	*
	*/
	bool longitudinal_;

	/**
	*
	*/
	double yukFactor_;

	/**
	*
	*/
	double boostFactor_;

	/**
	*
	*/
	double Q_;
	//@}

	/**
	* Parameters
	*/
	//@{
	/**
	* Order for the K terms
	*/
	int K_ORDER_;

	/**
	* Order for the B terms
	*/
	int B_ORDER_;
	//@}

	/**
	* Integrator
	*/
	GSLIntegrator integrator_;

	private:

	/**
	* Storage of the high scale pieces
	*/
	//@{
	/**
	*
	*/
	Complex highColW_;

	/**
	*
	*/
	Complex highColB_;

	/**
	*
	*/
	Complex highColG_;

	/**
	*
	*/
	Complex highColQ_;

	/**
	*
	*/
	Complex highColQt_;

	/**
	*
	*/
	Complex highColU_;

	/**
	*
	*/
	Complex highColtR_;

	/**
	*
	*/
	Complex highColD_;

	/**
	*
	*/
	Complex highColL_;

	/**
	*
	*/
	Complex highColE_;

	/**
	*
	*/
	Complex highColPhi_;
	//@}

	/**
	* Storage of the low scale pieces
	*/
	//@{
	/**
	*
	*/
	complex<double> lowColW_;

	/**
	*
	*/
	Complex lowColA_;

	/**
	*
	*/
	Complex lowColG_;

	/**
	*
	*/
	Complex lowColU_;

	/**
	*
	*/
	Complex lowColt_;

	/**
	*
	*/
	Complex lowColD_;

	/**
	*
	*/
	Complex lowColE_;
	//@}

	/**
	* Storage of the matching parameters
	*/
	//@{
	/**
	*
	*/
	Complex ULcollinearCorr_;

	/**
	*
	*/
	Complex DLcollinearCorr_;

	/**
	*
	*/
	Complex URcollinearCorr_;

	/**
	*
	*/
	Complex DRcollinearCorr_;

	/**
	*
	*/
	Complex tLcollinearCorr_;

	/**
	*
	*/
	Complex tRcollinearCorr_;

	/**
	*
	*/
	Complex bLcollinearCorr_;

	/**
	*
	*/
	Complex nuLcollinearCorr_;

	/**
	*
	*/
	Complex ELcollinearCorr_;

	/**
	*
	*/
	Complex ERcollinearCorr_;

	/**
	*
	*/
	Complex WtoWcollinearCorr_;

	/**
	*
	*/
	Complex WtoZcollinearCorr_;

	/**
	*
	*/
	Complex WtoAcollinearCorr_;

	/**
	*
	*/
	Complex BtoZcollinearCorr_;

	/**
	*
	*/
	Complex BtoAcollinearCorr_;

	/**
	*
	*/
	Complex PhitoWcollinearCorr_;

	/**
	*
	*/
	Complex PhitoZcollinearCorr_;

	/**
	*
	*/
	Complex PhitoHcollinearCorr_;

	/**
	*
	*/
	Complex GcollinearCorr_;
	//@}
	};

	}

	#endif /* Herwig_CollinearSudakov_H */
	diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
	--- a/MatrixElement/EW/GroupInvariants.h
	+++ b/MatrixElement/EW/GroupInvariants.h
	@@ -1,311 +1,328 @@
	// -- C++ --
	//
	// GroupInvariants.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	//
	#ifndef HERWIG_GroupInvariants_H
	#define HERWIG_GroupInvariants_H
	#include "ThePEG/Config/ThePEG.h"
	#include "ThePEG/Config/Unitsystem.h"
	#include <cassert>
	#include <boost/numeric/ublas/matrix.hpp>

	namespace Herwig {
	using namespace ThePEG;

	namespace GroupInvariants {

	/**
	* Simple struct for storing the different gauge contributions
	*/
	struct GaugeContributions {

	/**
	* Default Constructor
	*/
	GaugeContributions(double inSU3=0.,
	double inSU2=0., double inU1=0.)
	: SU3(inSU3),SU2(inSU2),U1(inU1)
	{}
	/**
	* \f$SU(3)\f$
	*/
	double SU3;

	/**
	* \f$SU(2)\f$
	*/
	double SU2;

	/**
	* \f$U(1)\f$
	*/
	double U1;
	};


	/**
	* The \f$SU(N)\f$ \f$C_A\f$
	*/
	inline double C_A(unsigned int N) {
	return N !=1 ? double(N) : 0.;
	}

	/**
	* The \f$SU(N)\f$ \f$C_F\f$
	*/
	inline double C_F(unsigned int N) {
	return N !=1 ? 0.5(double(NN)-1.)/double(N) : 1.;
	}

	/*
	* The \f$SU(N)\f$ \f$C_d\f$
	*/
	inline double C_d(unsigned int N) {
	return (double(N*N)-4.)/double(N);
	}

	/**
	* The \f$SU(N)\f$\f$C_1\f$
	*/
	inline double C_1(unsigned int N) {
	double N2(N*N);
	return 0.25*(N2-1.0)/N2;
	}

	/**
	* \f$T_F\f$
	*/
	inline double T_F(unsigned int N, bool high) {
	if(high) {
	return N !=1 ? 0.5 : 5.0/3.0;
	}
	else {
	return N !=1 ? 0.5 : 20.0/3.0;
	}
	}

	/**
	* \f$t_S\f$
	*/
	inline double t_S(unsigned int, bool ) {
	return 0.5;
	}

	/**
	* / Number of complex scalars in the fundamental rep. of SU(N)/U(1)
	*/
	inline double n_S(unsigned int N, bool high) {
	if(high) {
	if(N==2 \|\| N==1) return 1.0;
	else if(N==3) return 0.0;
	else assert(false);
	}
	else {
	if(N>=1&&N<=3) return 0.;
	else assert(false);
	}
	}

	/**
	* Number of Dirac Fermions in the fund. rep. of SU(N) (or U(1) for N==1)
	*/
	inline double n_F(unsigned int N, bool high) {
	if(high) {
	if(N==1) return 3.0;
	else if(N==2 \|\| N==3) return 6.0;
	else assert(false);
	}
	else {
	if(N==1) return 1.0;
	else if(N==2) return 0.0;
	else if(N==3) return 5.0;
	else assert(false);
	}
	}

	/**
	* Find K_i for gauge group N. high=false for running at mu<mZ
	*/
	double K_Factor(unsigned int i,unsigned int N, bool high);

	/**
	* Find B_i for gauge group N, high energy
	*/
	double B_Factor(int i, int N, bool fermion, bool longitudinal);

	/**
	* Find B_i for gauge group N, low energy
	*/
	double B_Factor_Low(int i, int N, bool fermion, double boostFactor);

	/**
	* Contributions to the Cusps
	*/
	GaugeContributions cuspContributions(Energy mu, int K_ORDER, bool high);

	/**
	* Contributions to B, high energy
	*/
	GaugeContributions BContributions(Energy mu, int B_ORDER,
	bool fermion, bool longitudinal);

	/**
	* Contributions to B, low energy
	*/
	GaugeContributions BContributionsLow(Energy mu, int B_ORDER,
	bool fermion, double boostFactor);

	inline Complex PlusLog(double arg) {
	static const Complex I(0,1.0);
	if (arg>0.0)
	return log(arg);
	else if (arg<0.0)
	return log(-arg)+I*Constants::pi;
	else
	assert(false);
	}

	inline Complex MinusLog(double arg) {
	static const Complex I(0,1.0);
	if (arg>0.0)
	return log(arg);
	else if (arg<0.0)
	return log(-arg)-I*Constants::pi;
	else
	assert(false);
	}

	inline Complex getT(Energy2 s, Energy2 t) {
	return MinusLog(-t/GeV2) - MinusLog(-s/GeV2);
	}

	inline Complex getU(Energy2 s, Energy2 u) {
	return MinusLog(-u/GeV2) - MinusLog(-s/GeV2);
	}

	inline boost::numeric::ublas::matrix<Complex> Gamma2(Complex U, Complex T) {
	boost::numeric::ublas::matrix<Complex> output(2,2);
	static const Complex I(0,1.0);
	using Constants::pi;
	output(0,0) = (-3.0/2.0)Ipi + (T+U);
	output(1,1) = (-3.0/2.0)Ipi;
	output(0,1) = 2.0*(T-U);
	output(1,0) = (3.0/8.0)*(T-U);
	return output;
	}

	inline boost::numeric::ublas::matrix<Complex> Gamma2w(Complex U, Complex T) {
	boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(5,5);
	static const Complex I(0,1.0);
	using Constants::pi;
	output(0,0) += -Ipi11.0/4.0;
	output(0,1) += U-T;
	output(1,0) += 2.0*(U-T);
	output(1,1) += -Ipi11.0/4.0 + (T+U);
	output(2,2) += -7.0/4.0Ipi + (U+T);
	output(3,3) += -7.0/4.0Ipi + (U+T);
	output(4,4) += -3.0/4.0Ipi;
	return output;
	}

	inline boost::numeric::ublas::matrix<Complex> Gamma2Singlet() {
	boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(2,2);
	static const Complex I(0,1.0);
	using Constants::pi;
	output(0,0) = output(1,1) = -3.0/4.0Ipi;
	return output;
	}

	+ inline Complex Gamma1(double hypercharge) {
	+ Complex I(0,1.0);
	+ return -IConstants::pi(hypercharge*hypercharge);
	+ }
	+
	+ inline Complex Gamma1(double y1, double y2, Complex T, Complex U) {
	+ Complex I(0,1.0);
	+ return -IConstants::pi(y1y1+y2y2) + 2.0y1y2*(T-U);
	+ }
	+
	+ inline Complex Gamma1(double y1, double y2, double y3, double y4,
	+ Complex T, Complex U) {
	+ Complex I(0,1.0);
	+ return -IConstants::pi(y1y1+y2y2+y3y3+y4y4)/2.0 +
	+ (y1y4+y2y3)T - (y1y3+y2y4)U;
	+ }
	+
	/**
	* Number of fermion generations (only used in gauge boson HighCMatching)
	*/
	inline double n_g() { return 3.0; }

	/**
	* Number of complex scalars in the fundamental rep. of SU(N)
	*/
	inline double nSWeyl(unsigned int N, bool high) {
	if(high) {
	if(N==2 \|\| N==1) return 1.0;
	else if (N==3) return 0.0;
	else assert(false);
	}
	else {
	if( N==1 \|\| N==3 ) return 0.0;
	else assert(false);
	}
	}

	/**
	* Number of Weyl Fermions in the fundamental rep. of SU(N)
	*/
	inline double nFWeyl(unsigned int N, bool high) {
	if(high) {
	if(N==2 \|\| N==3) return 12.0;
	else assert(false);
	}
	else {
	if(N==3) return 10.0;
	else if(N==1) return 2.0;
	else assert(false);
	}
	}

	inline double TFWeyl(unsigned int) {
	return 0.5;
	}

	inline double tSWeyl(unsigned int) {
	return 0.5;
	}

	inline Complex WFunction(Energy mu, Energy2 s) {
	using Constants::pi;
	assert(abs(s)>ZERO);
	Complex ln = MinusLog(-s/(mu*mu));
	return (-1.0lnln + 3.0ln+pipi/6.0-8.0);
	}

	/**
	* \fX_N\f% function, v is either t or u
	*/
	inline Complex XNFunction(unsigned int N, Energy mu, Energy2 s, Energy2 v) {
	using Constants::pi;
	assert(abs(s)>ZERO);
	Complex ls = MinusLog(-s/(mu*mu));
	return (2.0C_F(N)WFunction(mu,s) +
	C_A(N)(2.0lsls - 2.0MinusLog((s+v)/(mumu))ls -
	11.0/3.0ls + pipi + 85.0/9.0) +
	(2.0/3.0ls - 10.0/9.0) TFWeyl(N) * nFWeyl(N,true) +
	(1.0/3.0ls - 8.0/9.0) TFWeyl(N) * nSWeyl(N,true));
	}

	/**
	* \f$\Pi_1\f$ function
	*/
	inline Complex PI1_function(Energy mu, Energy2 s) {
	assert(abs(s)>ZERO);
	return ((41.0/6.0)MinusLog(-s/(mumu))-104.0/9.0);
	}

	/**
	* \f$\tilde{f}\f$ function, v is either t or u
	*/
	inline Complex fTildeFunction(Energy mu, Energy2 s, Energy2 v) {
	using Constants::pi;
	assert(abs(s)>ZERO);
	Complex ls = MinusLog(-s/GeV2), lv = MinusLog(-v/GeV2);
	Complex lsv = MinusLog((s+v)/GeV2);
	return (-2.0double(s/(s+v))(lv-ls) +
	double(s(s+2.0v)/((s+v)(s+v))) ((lv-ls)(lv-ls) + pipi) +
	4.0MinusLog(-s/(mumu))*(lv-lsv));
	}
	}
	}

	#endif // HERWIG_GroupInvariants_H

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 4:47 PM (22 h, 3 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023513
Default Alt Text: (83 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions