No OneTemporary
Actions

Size

98 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,657 +1,1192 @@
	// -- 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));
	+ }
	+ }
	+
	+ 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,513 +1,527 @@
	// -- 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
	+ */
	+ boost::numeric::ublas::matrix<Complex>
	+ electroWeakMatching(Energy EWScale, Energy2 s,
	+ Herwig::EWProcess::Process process,
	+ bool oneLoop);
	+
	+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);

	- /**
	- * Make plots for tests
	- */
	- void makePlots();
	-
	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/ElectroWeakMatching.cc b/MatrixElement/EW/ElectroWeakMatching.cc
	--- a/MatrixElement/EW/ElectroWeakMatching.cc
	+++ b/MatrixElement/EW/ElectroWeakMatching.cc
	@@ -1,917 +1,874 @@
	// -- 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;

	-namespace {
	-
	-Complex getT(Energy2 s, Energy2 t) {
	- return MinusLog(-t/GeV2) - MinusLog(-s/GeV2);
	-}
	-
	-Complex getU(Energy2 s, Energy2 u) {
	- return MinusLog(-u/GeV2) - MinusLog(-s/GeV2);
	-}
	-
	-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;
	-}
	-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;
	-}
	-
	-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;
	-}
	-}
	-
	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(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.5Ipi;
	}
	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.0g11g21(T-U) - Ipi(g11g11+g21*g21);
	Complex z2 = 2.0g21g12(T-U) - Ipi(g21g21+g12*g12);
	Complex z3 = 2.0g22g11(T-U) - Ipi(g22g22+g11*g11);
	Complex z4 = 2.0g12g22(T-U) - Ipi(g12g12+g22*g22);
	Complex z5 = (g11g21+g12g22)T - (g21g12+g11g22)U
	- 0.5Ipi(g11g11+g12g12+g21g21+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.25Ipi;
	for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;

	Complex z1 = 2.0g11g22(T-U) - Ipi(g11g11+g22*g22);
	Complex z2 = 2.0g12g22(T-U) - Ipi(g12g12+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.25Ipi;
	for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;

	Complex z1 = 2.0g11g22(T-U) - Ipi(g11g11+g22*g22);
	Complex z2 = 2.0g12g22(T-U) - Ipi(g12g12+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.5Ipi;
	}
	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.0g11g21(T-U) - Ipi(g11g11+g21*g21);
	Complex z2 = 2.0g21g12(T-U) - Ipi(g21g21+g12*g12);
	Complex z3 = 2.0g22g11(T-U) - Ipi(g22g22+g11*g11);
	Complex z4 = 2.0g12g22(T-U) - Ipi(g12g12+g22*g22);
	Complex z5 = (g11g21+g12g22)T - (g21g12+g11g22)U
	- 0.5Ipi(g11g11+g12g12+g21g21+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.25Ipi;
	Dw(0,0) = Dw(1,1) = w1;

	Complex z1 = 2.0g11g22(T-U) - Ipi(g11g11+g22*g22);
	Complex z2 = 2.0g12g22(T-U) - Ipi(g12g12+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.0g11g22(T-U) - Ipi(g11g11+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.0g11g22(T-U) - Ipi(g11g11+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.25Ipi;
	Dw(0,0) = Dw(1,1) = w1;

	Complex z1 = 2.0g11g22(T-U) - Ipi(g11g11+g22*g22);
	Complex z2 = 2.0g12g22(T-U) - Ipi(g12g12+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.0g11g22(T-U) - Ipi(g11g11+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.0g11g22(T-U) - Ipi(g11g11+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.25Ipi;
	Dw(0,0) = Dw(1,1) = w1;

	Complex z1 = 2.0g11g22(T-U) - Ipi(g11g11+g22*g22);
	Complex z2 = 2.0g12g22(T-U) - Ipi(g12g12+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.0g11g22(T-U) - Ipi(g11g11+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.5Ipi;
	}
	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.0g11g21(T-U) - Ipi(g11g11+g21*g21);
	Complex z2 = 2.0g21g12(T-U) - Ipi(g21g21+g12*g12);
	Complex z3 = 2.0g22g11(T-U) - Ipi(g22g22+g11*g11);
	Complex z4 = 2.0g12g22(T-U) - Ipi(g12g12+g22*g22);
	Complex z5 = (g11g21+g12g22)T - (g21g12+g11g22)U
	- 0.5Ipi(g11g11+g12g12+g21g21+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.25Ipi;
	Dw(0,0) = Dw(1,1) = w1;

	Complex z1 = 2.0g11g22(T-U) - Ipi(g11g11+g22*g22);
	Complex z2 = 2.0g12g22(T-U) - Ipi(g12g12+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.0g11g22(T-U) - Ipi(g11g11+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.5sincos; R0(2,3) = -0.5sincos; R0(2,4) = sin2;
	R0(3,0) = sincos; R0(3,2) = 0.5cos2; R0(3,3) = -0.5sin2; R0(3,4) = -sincos;
	R0(4,0) = sincos; R0(4,2) = -0.5sin2; R0(4,3) = 0.5cos2; R0(4,4) = -sincos;
	R0(5,0) = sin2; R0(5,2) = 0.5sincos; R0(5,3) = 0.5sincos; 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.5sincos; R0(8,3) = 0.5sincos; R0(8,4) = sin2;
	R0(9,0) = sincos; R0(9,2) = -0.5cos2; R0(9,3) = 0.5sin2; R0(9,4) = -sincos;
	R0(10,0) = sincos; R0(10,2) = 0.5sin2; R0(10,3) = -0.5cos2; R0(10,4) = -sincos;
	R0(11,0) = sin2; R0(11,2) = -0.5sincos; R0(11,3) = -0.5sincos; 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.0Ipi;
	Complex w2 = -T+U+5.0/4.0Ipi;
	Complex w3 = -0.5(T+U) + 3.0/4.0I*pi;
	Complex w4 = 0.25Ipi;
	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.0g1gW(U-T) - Ipi(g1g1+gW*gW);
	Complex z2 = 2.0g1gW(T-U) - Ipi(g1g1+gW*gW);
	Complex z3 = 2.0g2gW(U-T) - Ipi(g2g2+gW*gW);
	Complex z4 = 2.0g2gW(T-U) - Ipi(g2g2+gW*gW);
	Complex z5 = -(g2gW)T + (g1gW)U - Ipi(g1g2+g1gW-g2*gW);
	Complex z6 = (g1gW)T - (g2gW)U - Ipi(g1g2+g1gW-g2*gW);
	Complex z7 = -Ipig1*g1;
	Complex z8 = -Ipig2*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) {
	g1 = g_Lu;
	g2 = g_Ld;
	}
	else if (process==LLWW) {
	g1 = g_Lnu;
	g2 = g_Le;
	}
	else
	assert(false);
	double g3 = g_phiPlus;

	Complex w0 = 0.25Ipi;
	Complex w1 = 0.5(T-U) + 0.5I*pi;
	Complex w2 = -0.5(T-U) + 0.5I*pi;
	Complex w3 = 0.25I(T-U);
	Complex w4 = -0.25(T+U) + 0.25I*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.0g3g1(T-U) - Ipi(g3g3+g1*g1);
	Complex z2 = 2.0g3g2(T-U) - Ipi(g3g3+g2*g2);
	Complex z3 = -Ipig1*g1;
	Complex z4 = 0.5Ig1*(T-U);
	Complex z5 = 0.25Ipi;
	Complex z6 = -Ipig2*g2;
	Complex z7 = 0.5Ig2*(T-U);
	Complex z8 = g3g1T-g3g2U-Ipi(g1g2-g2g3+g1*g3);
	Complex z9 = 0.5Ig2T-0.5Ig1U+pi/2.0g2-pi/2.0g1+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.75I*pi;
	Complex w2 = 0.25Ipi;
	Dw(0,0) = Dw(1,1) = w1;
	Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;

	Complex z1 = gWg1T - gWg2U - Ipi(g1g2+g1gW-g2*gW);
	Complex z2 = gWg1U - gWg2T - Ipi(g2g1+g1gW-g2*gW);
	Complex z3 = -Ipig1*g1;
	Complex z4 = -Ipig2*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.0Ipi + (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.25Ipi;
	Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = w2;
	Complex z3 = -Ipig1*g1;
	Complex z4 = -Ipig2*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.25Ipi;
	Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;
	Complex z3 = -Ipig1*g1;
	Complex z4 = -Ipig2*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 = -Ipig1*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.25Ipi;
	Dw(2,3) = Dw(3,2) = -0.25Ipi;
	Complex z1 = 2.0g3g1(T-U) - Ipi(g3g3+g1*g1);
	Complex z2 = 0.5Ig1*g1;
	Complex z3 = -Ipig1*g1;
	Complex z4 = 0.25Ipi;
	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 = -Ipig1*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);
	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 = -Ipig1*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> temp(R0.size1(),R0.size2());
	boost::numeric::ublas::axpy_prod(R0,G2,temp);
	R0+=aW/(4.0pi)4.0log(mW/mu)temp;
	boost::numeric::ublas::axpy_prod(Dw,R0,temp);
	R0+=aW/(4.0pi)4.0log(mW/mu)temp;
	boost::numeric::ublas::axpy_prod(Dz,R0,temp);
	R0+=aZ/(4.0pi)4.0log(mZ/mu)temp;
	return R0;
	}
	diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
	--- a/MatrixElement/EW/GroupInvariants.h
	+++ b/MatrixElement/EW/GroupInvariants.h
	@@ -1,269 +1,311 @@
	// -- 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;
	+ }
	+
	/**
	* 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:51 PM (18 h, 12 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4017665
Default Alt Text: (98 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions