Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8309840
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
71 KB
Subscribers
None
View Options
diff --git a/MatrixElement/EW/CollinearSudakov.cc b/MatrixElement/EW/CollinearSudakov.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/CollinearSudakov.cc
@@ -0,0 +1,657 @@
+// -*- 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;
+
+CollinearSudakov::CollinearSudakov() : K_ORDER_(3),
+ B_ORDER_(2)
+{}
+
+CollinearSudakov::~CollinearSudakov() {}
+
+IBPtr CollinearSudakov::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr CollinearSudakov::fullclone() const {
+ return new_ptr(*this);
+}
+
+void CollinearSudakov::persistentOutput(PersistentOStream & os) const {
+ os << K_ORDER_ << B_ORDER_;
+}
+
+void CollinearSudakov::persistentInput(PersistentIStream & is, int) {
+ is >> K_ORDER_ >> B_ORDER_;
+}
+
+
+// The following static variable is needed for the type
+// description system in ThePEG.
+DescribeClass<CollinearSudakov,Interfaced>
+describeHerwigCollinearSudakov("Herwig::CollinearSudakov", "HwMEEW.so");
+
+void CollinearSudakov::Init() {
+
+ static ClassDocumentation<CollinearSudakov> documentation
+ ("The CollinearSudakov class implements the collinear evolution");
+
+
+}
+
+
+Complex CollinearSudakov::highScaleIntegral(bool SU3, bool SU2, double Y,
+ Energy2 s, Energy mu_h, Energy mu_l,
+ bool fermion, bool longitudinal,
+ double yukFactor) {
+ SU3_ = SU3;
+ SU2_ = SU2;
+ Y_ = Y;
+ s_ = s;
+ fermion_ = fermion;
+ longitudinal_ = longitudinal;
+ yukFactor_ = yukFactor;
+ // perform the integral
+ Complex result;
+ high_ = true;
+ // real piece
+ real_ = true;
+ result.real(integrator_.value(*this,mu_h,mu_l));
+ // imaginary piece
+ real_ = false;
+ result.imag(integrator_.value(*this,mu_h,mu_l));
+ // return the answer
+ return exp(result);
+}
+
+Complex CollinearSudakov::lowScaleIntegral(bool SU3, double Q, Energy2 s,
+ Energy mu_h, Energy mu_l, bool fermion,
+ double boostFactor) {
+ SU3_ = SU3;
+ Q_ = Q;
+ s_ = s;
+ fermion_ = fermion;
+ boostFactor_ = boostFactor;
+ // perform the integral
+ Complex result;
+ high_ = false;
+ // real piece
+ real_ = true;
+ result.real(integrator_.value(*this,mu_h,mu_l));
+ // imaginary piece
+ real_ = false;
+ result.imag(integrator_.value(*this,mu_h,mu_l));
+ // return the answer
+ return exp(result);
+}
+
+InvEnergy CollinearSudakov::highScaleIntegrand(Energy mu) const {
+ using namespace GroupInvariants;
+ using Constants::pi;
+ Complex gamma = 0.0;
+ // Include K-factor Contributions (Cusps):
+ GaugeContributions cusp = cuspContributions(mu,K_ORDER_,high_);
+ // Include B-factors (B1 for U1, B2 for SU2, B3 for SU3):
+ GaugeContributions nonCusp = BContributions(mu,B_ORDER_,fermion_,longitudinal_);
+ // common log
+ Complex plog = PlusLog(s_/sqr(mu));
+ // SU(3)
+ if(SU3_) {
+ if(fermion_)
+ gamma += C_F(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
+ else
+ gamma += (C_A(3))*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
+ }
+ // SU(2)
+ if(SU2_) {
+ if (fermion_ || longitudinal_ )
+ gamma += (C_F(2))*cusp.SU2*0.5*plog + 0.5*nonCusp.SU2;
+ else
+ gamma += (C_A(2))*cusp.SU2*0.5*plog + 0.5*nonCusp.SU2;
+ }
+
+ if (fermion_ || longitudinal_ ) {
+ gamma += sqr(Y_)*(cusp.U1*0.5*plog + 0.5*nonCusp.U1);
+ }
+ else {
+ // U(1) Gauge boson
+ if (!SU3_ && !SU2_ && abs(Y_)<0.001) {
+ gamma += 0.5*nonCusp.U1;
+ }
+ }
+ // top Yukawa piece
+ double y_t = ElectroWeakReweighter::coupling()->y_t(mu);
+ gamma += yukFactor_*sqr(y_t)/(16.0*sqr(pi));
+ // return the answer
+ return real_ ? gamma.real()/mu : gamma.imag()/mu;
+}
+
+InvEnergy CollinearSudakov::lowScaleIntegrand(Energy mu) const {
+ using namespace GroupInvariants;
+ using Constants::pi;
+ Complex gamma = 0.0;
+ // Include K-factor Contributions (Cusps):
+ GaugeContributions cusp = cuspContributions(mu,K_ORDER_,false);
+ // Include B-factors (B1 for U1, B2 for SU2, B3 for SU3):
+ GaugeContributions nonCusp = BContributionsLow(mu,B_ORDER_,fermion_,boostFactor_);
+ // common log
+ Complex plog = PlusLog(s_/sqr(mu));
+ Complex blog(0.);
+ if(boostFactor_ >0.001) blog = PlusLog(4.0*boostFactor_);
+ // SU(3)
+ if (SU3_) {
+ if (fermion_) {
+ // not a bHQET top quark field
+ if (abs(boostFactor_)<0.001)
+ gamma += C_F(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
+ else
+ gamma += C_F(3)*cusp.SU3*0.5*blog + 0.5*nonCusp.SU3;
+ }
+ else {
+ gamma += C_A(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
+ }
+ }
+ // fermions
+ if (fermion_) {
+ // not a bHQET top quark field
+ if (abs(boostFactor_)<0.001)
+ gamma += sqr(Q_)*(cusp.U1*0.5*plog + 0.5*nonCusp.U1);
+ else
+ gamma += sqr(Q_)*(cusp.U1*0.5*blog + 0.5*nonCusp.U1);
+ }
+ else {
+ // i.e., not a fermion, not a bHQ, not a gluon => photon
+ if (abs(boostFactor_)<0.001 && !SU3_)
+ gamma += 0.5*nonCusp.U1;
+ // i.e., W treated as a bHQET field
+ else if (abs(boostFactor_)>0.001)
+ gamma += sqr(Q_)*(cusp.U1*0.5*blog + 0.5*nonCusp.U1);
+ }
+ // return the answer
+ return real_ ? gamma.real()/mu : gamma.imag()/mu;
+}
+
+void CollinearSudakov::evaluateHighScale(Energy highScale, Energy EWScale, Energy2 S) {
+ double yCoeffQt(0.), yCoefftR(0.), yCoeffPhi(0.);
+ if (K_ORDER_ >= 1) {
+ yCoeffQt = 0.5;
+ yCoefftR = 1.0;
+ yCoeffPhi = 3.0;
+ }
+ highColW_ = highScaleIntegral(false,true ,0.0 ,S,highScale,EWScale,false,false,0.0);
+ highColB_ = highScaleIntegral(false,false,0.0 ,S,highScale,EWScale,false,false,0.0);
+ highColG_ = highScaleIntegral(true ,false,0.0 ,S,highScale,EWScale,false,false,0.0);
+ highColQ_ = highScaleIntegral(true ,true ,1./6. ,S,highScale,EWScale,true,false,0.0);
+ highColQt_ = highScaleIntegral(true ,true ,1./6. ,S,highScale,EWScale,true,false,yCoeffQt);
+ highColU_ = highScaleIntegral(true ,false,2./3. ,S,highScale,EWScale,true,false,0.0);
+ highColtR_ = highScaleIntegral(true ,false,2./3. ,S,highScale,EWScale,true,false,yCoefftR);
+ highColD_ = highScaleIntegral(true ,false,-1./3.,S,highScale,EWScale,true,false,0.0);
+ highColL_ = highScaleIntegral(false,true ,-1./2.,S,highScale,EWScale,true,false,0.0);
+ highColE_ = highScaleIntegral(false,false,-1. ,S,highScale,EWScale,true,false,0.0);
+ highColPhi_ = highScaleIntegral(false,true ,1./2. ,S,highScale,EWScale,false,true,yCoeffPhi);
+}
+
+void CollinearSudakov::evaluateLowScale(Energy EWScale, Energy lowScale, Energy2 S) {
+ lowColW_ = lowScaleIntegral(false,1. ,S,EWScale,lowScale,false,
+ S/sqr(2.*ElectroWeakReweighter::coupling()->mW()));
+ lowColA_ = lowScaleIntegral(false,0. ,S,EWScale,lowScale,false,0.0);
+ lowColG_ = lowScaleIntegral(true ,0. ,S,EWScale,lowScale,false,0.0);
+ lowColU_ = lowScaleIntegral(true ,2./3. ,S,EWScale,lowScale,true,0.0);
+ lowColt_ = lowScaleIntegral(true ,2./3. ,S,EWScale,lowScale,true,
+ S/sqr(2.*ElectroWeakReweighter::coupling()->mT()));
+ lowColD_ = lowScaleIntegral(true ,-1./3.,S,EWScale,lowScale,true,0.0);
+ lowColE_ = lowScaleIntegral(false,-1. ,S,EWScale,lowScale,true,0.0);
+}
+
+void CollinearSudakov::evaluateMatching(Energy EWScale,Energy2 s) {
+ using Constants::pi;
+ using GroupInvariants::PlusLog;
+ static const Complex I(0,1.0);
+ // wave function corrections
+ WaveFunctionCorrections WFC = waveFunctionCorrections(EWScale);
+
+ double aS = ElectroWeakReweighter::coupling()->a3(EWScale);
+ double aEM = ElectroWeakReweighter::coupling()->aEM(EWScale);
+ double aW = ElectroWeakReweighter::coupling()->aW(EWScale);
+ double aZ = ElectroWeakReweighter::coupling()->aZ(EWScale);
+ double cW2 = ElectroWeakReweighter::coupling()->Cos2thW(EWScale);
+ double sW2 = ElectroWeakReweighter::coupling()->Sin2thW(EWScale);
+ Energy mW = ElectroWeakReweighter::coupling()->mW();
+ Energy mZ = ElectroWeakReweighter::coupling()->mZ();
+ Energy mH = ElectroWeakReweighter::coupling()->mH();
+ Energy mT = ElectroWeakReweighter::coupling()->mT();
+ double gPHI = ElectroWeakReweighter::coupling()->g_phiPlus(EWScale);
+ double y_t = ElectroWeakReweighter::coupling()->y_t(EWScale);
+
+ double lz = log(mZ/EWScale);
+ double lw = log(mW/EWScale);
+
+ complex<double> F_W = 4.0*lw*0.5*PlusLog(s/EWScale/EWScale)-2.0*lw*lw-2.0*lw-5.0*pi*pi/12.0+1.0;
+ complex<double> F_Z = 4.0*lz*0.5*PlusLog(s/EWScale/EWScale)-2.0*lz*lz-2.0*lz-5.0*pi*pi/12.0+1.0;
+
+ // Taken from Manohar... along with his formulae for F_tL, F_tR, and F_bL (for 3rd/1st Gen. Wavefunction Differences)
+ complex<double> W1 = WFC.fFW0-0.5*WFC.aW0-0.5*WFC.cW0;
+ complex<double> W2 = WFC.fF0W-0.5*WFC.a0W;
+ complex<double> U1 = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*(WFC.fFZZ-0.5*WFC.aZZ) -
+ 0.5*WFC.cZZ*(ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale) +
+ ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)) +
+ ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*WFC.bZZ;
+ complex<double> U2 = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*(WFC.fFZZ-0.5*WFC.aZZ) -
+ 0.5*WFC.cZZ*(ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale) +
+ ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)) +
+ ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*WFC.bZZ;
+ complex<double> HtL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-0.5*log(mH/EWScale)-0.5*lz+
+ 0.5*WFC.atHH+0.5*WFC.atZZ+WFC.ctHH+WFC.ctZZ+
+ WFC.ctW0-WFC.btHH+WFC.btZZ);
+ complex<double> HtR = HtL-0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.atW0);
+ complex<double> HbL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.at0W);
+
+ complex<double> F_tL = (4.0/3.0*aS+4.0/9.0*aEM)/(4.0*pi)*(2.0*log(mT/EWScale)*log(mT/EWScale)-log(mT/EWScale)+
+ pi*pi/12.0+2.0) + HtL +
+ aW/(4.0*pi)*0.5*W1 + aZ/(4.0*pi)*U1;
+
+ complex<double> F_tR = (4.0/3.0*aS+4.0/9.0*aEM)/(4.0*pi)*(2.0*log(mT/EWScale)*log(mT/EWScale)-log(mT/EWScale)+
+ pi*pi/12.0+2.0) + HtR -
+ aW/(4.0*pi)*0.25*WFC.cW0 + aZ/(4.0*pi)*U2;
+
+ complex<double> F_bL = HbL + aW/(4.0*pi)*0.5*W2;
+
+ Complex Dw = CollinearDw(s,EWScale);
+ Complex Dz = CollinearDz(s,EWScale);
+
+ complex<double> D_C_UL = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*Dz + 0.5*Dw;
+ complex<double> D_C_DL = ElectroWeakReweighter::coupling()->g_Ld(EWScale)*ElectroWeakReweighter::coupling()->g_Ld(EWScale)*Dz + 0.5*Dw;
+
+ complex<double> D_C_UR = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*Dz;
+ complex<double> D_C_DR = ElectroWeakReweighter::coupling()->g_Rd(EWScale)*ElectroWeakReweighter::coupling()->g_Rd(EWScale)*Dz;
+
+ complex<double> D_C_nuL = ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*Dz + 0.5*Dw;
+ complex<double> D_C_EL = ElectroWeakReweighter::coupling()->g_Le(EWScale)*ElectroWeakReweighter::coupling()->g_Le(EWScale)*Dz + 0.5*Dw;
+ complex<double> D_C_ER = ElectroWeakReweighter::coupling()->g_Re(EWScale)*ElectroWeakReweighter::coupling()->g_Re(EWScale)*Dz;
+
+ complex<double> D_C_WW = aW/(4.0*pi)*cW2*(F_Z+WFC.fsWZWZ) +
+ aW/(4.0*pi)*cW2*(F_W+WFC.fs1ZW) + aW/(4.0*pi)*sW2*(F_W+WFC.fs10) +
+ aW/(4.0*pi)*sW2*(2.0*lw*lw-
+ 2.0*lw+pi*pi/12.0+2.0) + 0.5*WFC.RW;
+ complex<double> D_C_WZ = aW/(4.0*pi)*2.0*(F_W+WFC.fsZW1) + 0.5*WFC.RZ + sqrt(sW2/cW2)*WFC.RAtoZ;
+ complex<double> D_C_WA = aW/(4.0*pi)*2.0*(F_W+WFC.fs01) + 0.5*WFC.RA + sqrt(cW2/sW2)*WFC.RZtoA;
+ complex<double> D_C_BZ = 0.5*WFC.RZ - sqrt(cW2/sW2)*WFC.RAtoZ;
+ complex<double> D_C_BA = 0.5*WFC.RA - sqrt(sW2/cW2)*WFC.RZtoA;
+
+ // The WFC.RW and WFC.RZ are used on purpose (instead of WFC.RPhi and WFC.RPhi3, respectively):
+ complex<double> D_C_PhiW = aW/(4.0*pi)*0.25*(F_W+WFC.fs1HW) +
+ aW/(4.0*pi)*0.25*(F_W+WFC.fs1ZW) + aZ/(4.0*pi)*gPHI*gPHI*(F_Z+WFC.fsWZWZ) +
+ aW/(4.0*pi)*sW2*(2.0*lw*lw-2.0*lw+pi*pi/12.0+2.0) +
+ 0.5*WFC.RW + WFC.EW;
+ complex<double> D_C_PhiZ = aW/(4.0*pi)*0.5*(F_W+WFC.fsZW1) + aZ/(4.0*pi)*0.25*(F_Z+WFC.fs1HZ) + 0.5*WFC.RZ + WFC.EZ;
+ complex<double> D_C_PhiH = aW/(4.0*pi)*0.5*(F_W+WFC.fsHW1) + aZ/(4.0*pi)*0.25*(F_Z+WFC.fsHZ1) + 0.5*WFC.RH;
+
+ complex<double> D_C_GG = aS/(4.0*pi)*2.0/3.0*log(mT/EWScale);
+
+ ULcollinearCorr_ = exp(D_C_UL);
+ DLcollinearCorr_ = exp(D_C_DL);
+ URcollinearCorr_ = exp(D_C_UR);
+ DRcollinearCorr_ = exp(D_C_DR);
+
+ tLcollinearCorr_ = exp(D_C_UL+F_tL);
+ tRcollinearCorr_ = exp(D_C_UR+F_tR);
+ bLcollinearCorr_ = exp(D_C_DL+F_bL);
+
+ nuLcollinearCorr_ = exp(D_C_nuL);
+ ELcollinearCorr_ = exp(D_C_EL);
+ ERcollinearCorr_ = exp(D_C_ER);
+
+ WtoWcollinearCorr_ = exp(D_C_WW);
+ WtoZcollinearCorr_ = exp(D_C_WZ);
+ WtoAcollinearCorr_ = exp(D_C_WA);
+ BtoZcollinearCorr_ = exp(D_C_BZ);
+ BtoAcollinearCorr_ = exp(D_C_BA);
+
+ PhitoWcollinearCorr_ = exp(D_C_PhiW);
+ PhitoZcollinearCorr_ = exp(D_C_PhiZ);
+ PhitoHcollinearCorr_ = exp(D_C_PhiH);
+
+ GcollinearCorr_ = exp(D_C_GG);
+}
+
+WaveFunctionCorrections CollinearSudakov::waveFunctionCorrections(Energy EWScale) {
+ static const Complex I(0.,1.);
+ using Constants::pi;
+ double lZ = 2.0*log(ElectroWeakReweighter::coupling()->mZ()/EWScale);
+ WaveFunctionCorrections WFC;
+ // From Manohar, 2/12/12: (these assume mH=125, mZ=91.1876, mW=80.399)
+ WFC.RW = (0.8410283377963967 - 9.424777961271568*I) + 1.2785863646210789*lZ;
+ WFC.RA = (1.4835982362022198 + 1.855775680704845*pow(10.,-9)*I) - 0.27209907467584415*lZ;
+ WFC.RZ = (1.5114724841549798 - 9.926944419863688*I) + 1.0834802397165764*lZ;
+ WFC.RAtoZ = (0.3667485032811715 - 2.2770907130064835*I) - 1.2994544609942593*lZ;
+ WFC.RZtoA = -0.2095310079712942 + 0.8320191107808546*lZ;
+ WFC.RH = (12.229832449946716 - 1.7643103462419842*10.0*pow(10.,-12)*I) + 5.309998583664737*lZ;
+ WFC.RPhi = (5.569012418081201 + 1.5439133581417356*0.10*pow(10.,-9)*I) + 5.309998583664737*lZ;
+ WFC.RPhi3 = (8.945333042265943 + 5.499309445612249*pow(10.,-12)*I) + 5.309998583664737*lZ;
+ WFC.EW = (3.967645734304811 + 4.712388980384717*I) + 2.238332625165702*lZ;
+ WFC.EZ = (5.916079892937651 + 4.96347220970469*I) + 2.1132591719740788*lZ;
+
+ WFC.RW *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.RA *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.RZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.RAtoZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.RZtoA *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.RPhi *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.EW *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.EZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.RPhi3 *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+ WFC.RH *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
+
+
+ WFC.RW = WFC.RW.real();
+ WFC.RA = WFC.RA.real();
+ WFC.RZ = WFC.RZ.real();
+ WFC.RAtoZ = WFC.RAtoZ.real();
+ WFC.RZtoA = WFC.RZtoA.real();
+ WFC.RPhi = WFC.RPhi.real();
+ WFC.RPhi3 = WFC.RPhi3.real();
+ WFC.RH = WFC.RH.real();
+
+ // Also from Manohar, 2/10/12
+ WFC.fFW0 = -3.7946842553189453 - 4.709019671388589*I;
+ WFC.fF0W = 3.8181790485161176;
+ WFC.fFZZ = 1.364503250377989 + 0.*I;
+ WFC.aHH = -0.474396977740686 + 0.*I;
+ WFC.aZZ = -0.6806563210877914 + 0.*I;
+ WFC.aW0 = 0.49036102506811907 + 1.9323351450971642*I;
+ WFC.a0W = -1.2184776671065072;
+ WFC.bHH = 1.234775745474991 + 0.*I;
+ WFC.bZZ = 1.7303526426747613 + 0.*I;
+ WFC.cHH = 0.33140608274387473 + 0.*I;
+ WFC.cZZ = 0.4961363208382017 + 0.*I;
+ WFC.cW0 = -1.005299829777395 + 1.063048500002757*I;
+ WFC.atHH = -0.237198488870343 + 0.*I;
+ WFC.atZZ = -0.3403281605438957 + 0.*I;
+ WFC.atW0 = 0.24518051253405954 + 0.9661675725485821*I;
+ WFC.at0W = -0.6092388335532536;
+ WFC.ctHH = 0.16570304137193737 + 0.*I;
+ WFC.ctZZ = 0.24806816041910085 + 0.*I;
+ WFC.ctW0 = -0.5026499148886975 + 0.5315242500013785*I;
+ WFC.btHH = -0.30869393636874776 + 0.*I;
+ WFC.btZZ = -0.4325881606686903 + 0.*I;
+
+ WFC.fs10 = -2.289868133696459;
+ WFC.fs1ZW = 1.8627978596240622;
+ WFC.fsWZWZ = 1.1866922529667008;
+ WFC.fsZW1 = 1.0840307611156266;
+ WFC.fs01 = 2.2898681336964595;
+ WFC.fsHW1 = -0.32306745562682404;
+ WFC.fsHZ1 = 0.4042992116255279;
+ WFC.fs1HW = 3.330954543719127;
+ WFC.fs1HZ = 2.69768201932412;
+ return WFC;
+}
+
+Complex CollinearSudakov::CollinearDw(Energy2 s, Energy EWScale) {
+ using Constants::pi;
+ using GroupInvariants::PlusLog;
+ double lw = log(ElectroWeakReweighter::coupling()->mW()/EWScale);
+ //s = s/2.; // This should not be here... I think this is a discrepency with Sascha
+ return ElectroWeakReweighter::coupling()->aW(EWScale)/(4.0*pi)*
+ (4.0*lw*0.5*PlusLog(s/EWScale/EWScale) - 2.0*lw*lw -
+ 3.0*lw - 5.0/12.0*pi*pi + 9.0/4.0);
+}
+
+Complex CollinearSudakov::CollinearDz(Energy2 s, Energy EWScale) {
+ using Constants::pi;
+ using GroupInvariants::PlusLog;
+ double lz = log(ElectroWeakReweighter::coupling()->mZ()/EWScale);
+ return ElectroWeakReweighter::coupling()->aZ(EWScale)/(4.0*pi)*
+ (4.0*lz*0.5*PlusLog(s/EWScale/EWScale) - 2.0*lz*lz -
+ 3.0*lz - 5.0/12.0*pi*pi + 9.0/4.0);
+}
+
+namespace {
+
+void writeLine(ofstream & file, vector<Energy> x, vector<double> y,
+ string name,string title,
+ string colour, string style) {
+ file << "# BEGIN HISTO1D "+name +"\n";
+ file << "SmoothLine=1\n";
+ file << "ErrorBars=0\n";
+ file << "Path="+name+"\n";
+ file << "Title="+title+"\n";
+ file << "LineColor="+colour+"\n";
+ file << "LineStyle="+style +"\n";
+ file << "# xlow xhigh val errminus errplus\n";
+ for(unsigned int ix=0;ix<x.size();++ix)
+ file << (x[ix]-MeV)/GeV << "\t" << (x[ix]+MeV)/GeV << "\t"
+ << y[ix] << "\t" << 0. << "\t" << 0. << "\n";
+ file << "# END HISTO1D\n";
+}
+
+}
+
+void CollinearSudakov::makePlots() {
+ vector<Energy> np;
+ vector<double> uL,uR,tL,tR,dL,dR,bL,bR,nuL,eL,eR,Wgamma,Bgamma,g,WT,WL,WZT,BZT,ZL,H;
+ Energy mZ = ElectroWeakReweighter::coupling()->mZ();
+ for(Energy x=200.*GeV;x<5000.*GeV;x*=1.02) {
+ Energy2 s(sqr(x));
+ np.push_back(x);
+ evaluateHighScale(x,mZ,s);
+ evaluateMatching(mZ,s);
+ uL .push_back(real(highColQ_ * ULcollinearCorr_ ));
+ uR .push_back(real(highColU_ * URcollinearCorr_ ));
+ tL .push_back(real(highColQt_ * tLcollinearCorr_ ));
+ tR .push_back(real(highColtR_ * tRcollinearCorr_ ));
+ dL .push_back(real(highColQ_ * DLcollinearCorr_ ));
+ dR .push_back(real(highColD_ * DRcollinearCorr_ ));
+ bL .push_back(real(highColQt_ * bLcollinearCorr_ ));
+ bR .push_back(real(highColD_ * DRcollinearCorr_ ));
+ nuL .push_back(real(highColL_ * nuLcollinearCorr_ ));
+ eL .push_back(real(highColL_ * ELcollinearCorr_ ));
+ eR .push_back(real(highColE_ * ERcollinearCorr_ ));
+ Wgamma.push_back(real(highColW_ * WtoAcollinearCorr_ ));
+ Bgamma.push_back(real(highColB_ * BtoAcollinearCorr_ ));
+ g .push_back(real(highColG_ * GcollinearCorr_ ));
+ WT .push_back(real(highColW_ * WtoWcollinearCorr_ ));
+ WL .push_back(real(highColPhi_ * PhitoWcollinearCorr_));
+ WZT .push_back(real(highColW_ * WtoZcollinearCorr_ ));
+ BZT .push_back(real(highColB_ * BtoZcollinearCorr_ ));
+ ZL .push_back(real(highColPhi_ * PhitoZcollinearCorr_));
+ H .push_back(real(highColPhi_ * PhitoHcollinearCorr_));
+ }
+ ofstream fig1a("fig1a.dat");
+ fig1a << "#BEGIN PLOT\n";
+ fig1a << "LegendOnly=/uL /uR /tL /tR\n";
+ fig1a << "DrawOnly=/uL /uR /tL /tR\n";
+ fig1a << "Title=Figure 1a\n";
+ fig1a << "RatioPlot=0\n";
+ fig1a << "LogX=1\n";
+ fig1a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig1a << "YLabel=u, t\n";
+ fig1a << "Legend=1\n";
+ fig1a << "XMin=250.\n";
+ fig1a << "YMin=0.7\n";
+ fig1a << "YMax=1.05\n";
+ fig1a << "# END PLOT\n";
+ writeLine(fig1a,np,uL,"/uL","$u_L$","green","dotted" );
+ writeLine(fig1a,np,uR,"/uR","$u_R$","cyan" ,"solid" );
+ writeLine(fig1a,np,tL,"/tL","$t_L$","red" ,"dashed" );
+ writeLine(fig1a,np,tR,"/tR","$t_R$","blue" ,"dotdashed");
+ fig1a.close();
+ ofstream fig1b("fig1b.dat");
+ fig1b << "#BEGIN PLOT\n";
+ fig1b << "LegendOnly=/dL /dR /bL /bR\n";
+ fig1b << "DrawOnly=/dL /dR /bL /bR\n";
+ fig1b << "Title=Figure 1b\n";
+ fig1b << "RatioPlot=0\n";
+ fig1b << "LogX=1\n";
+ fig1b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig1b << "YLabel=d, b\n";
+ fig1b << "Legend=1\n";
+ fig1b << "XMin=250.\n";
+ fig1b << "YMin=0.7\n";
+ fig1b << "YMax=1.05\n";
+ fig1b << "# END PLOT\n";
+ writeLine(fig1b,np,dL,"/dL","$d_L$","green","dotted" );
+ writeLine(fig1b,np,dR,"/dR","$d_R$","cyan" ,"solid" );
+ writeLine(fig1b,np,bL,"/bL","$b_L$","red" ,"dashed" );
+ writeLine(fig1b,np,bR,"/bR","$b_R$","blue" ,"dotdashed");
+ fig1b.close();
+ ofstream fig2("fig2.dat");
+ fig2 << "#BEGIN PLOT\n";
+ fig2 << "LegendOnly=/uL /uR /dL /dR\n";
+ fig2 << "DrawOnly=/uL /uR /dL /dR\n";
+ fig2 << "Title=Figure 2\n";
+ fig2 << "RatioPlot=0\n";
+ fig2 << "LogX=1\n";
+ fig2 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig2 << "YLabel=u, d\n";
+ fig2 << "Legend=1\n";
+ fig2 << "XMin=250.\n";
+ fig2 << "YMin=0.7\n";
+ fig2 << "YMax=1.05\n";
+ fig2 << "# END PLOT\n";
+ writeLine(fig2,np,uL,"/uL","$u_L$","green","dotted" );
+ writeLine(fig2,np,uR,"/uR","$u_R$","cyan" ,"solid" );
+ writeLine(fig2,np,dL,"/dL","$d_L$","red" ,"dashed" );
+ writeLine(fig2,np,dR,"/dR","$d_R$","blue" ,"dotdashed");
+ fig2.close();
+ ofstream fig3("fig3.dat");
+ fig3 << "#BEGIN PLOT\n";
+ fig3 << "LegendOnly=/nuL /eL /eR\n";
+ fig3 << "DrawOnly=/nuL /eL /eR\n";
+ fig3 << "Title=Figure 3\n";
+ fig3 << "RatioPlot=0\n";
+ fig3 << "LogX=1\n";
+ fig3 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig3 << "YLabel=$\\nu$, $e$\n";
+ fig3 << "Legend=1\n";
+ fig3 << "XMin=250.\n";
+ fig3 << "YMin=0.9\n";
+ fig3 << "YMax=1.05\n";
+ fig3 << "# END PLOT\n";
+ writeLine(fig3,np,nuL,"/nuL","$\\nu_L$","blue","dashed");
+ writeLine(fig3,np, eL,"/eL" ,"$e_L$" ,"red" ,"dotted");
+ writeLine(fig3,np, eR,"/eR" ,"$e_R$" ,"red" ,"solid" );
+ fig3.close();
+ ofstream fig5("fig5.dat");
+ fig5 << "#BEGIN PLOT\n";
+ fig5 << "LegendOnly=/g /Wgamma /Bgamma\n";
+ fig5 << "DrawOnly=/g /Wgamma /Bgamma\n";
+ fig5 << "Title=Figure 5\n";
+ fig5 << "RatioPlot=0\n";
+ fig5 << "LogX=1\n";
+ fig5 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig5 << "YLabel=$\\gamma$, g\n";
+ fig5 << "Legend=1\n";
+ fig5 << "XMin=250.\n";
+ fig5 << "YMin=0.7\n";
+ fig5 << "YMax=1.05\n";
+ fig5 << "# END PLOT\n";
+ writeLine(fig5,np,g ,"/g" ,"$g$" ,"blue","dashed");
+ writeLine(fig5,np,Wgamma,"/Wgamma","$W\\to\\gamma$","red" ,"solid" );
+ writeLine(fig5,np,Bgamma,"/Bgamma","$B\\to\\gamma$","red" ,"dotted");
+ fig5.close();
+ ofstream fig6a("fig6a.dat");
+ fig6a << "#BEGIN PLOT\n";
+ fig6a << "LegendOnly=/WT /WL\n";
+ fig6a << "DrawOnly=/WT /WL\n";
+ fig6a << "Title=Figure 6a\n";
+ fig6a << "RatioPlot=0\n";
+ fig6a << "LogX=1\n";
+ fig6a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig6a << "YLabel=$Z_L$, $Z_T$, $H$\n";
+ fig6a << "Legend=1\n";
+ fig6a << "XMin=250.\n";
+ fig6a << "YMin=0.7\n";
+ fig6a << "YMax=1.05\n";
+ fig6a << "# END PLOT\n";
+ writeLine(fig6a,np,WT,"/WT","$W_T$","red" ,"solid");
+ writeLine(fig6a,np,WL,"/WL","$W_L$","blue" ,"dashed" );
+ fig6a.close();
+
+
+ ofstream fig6b("fig6b.dat");
+ fig6b << "#BEGIN PLOT\n";
+ fig6b << "LegendOnly=/WZT /BZT /ZL /H\n";
+ fig6b << "DrawOnly=/WZT /BZT /ZL /H\n";
+ fig6b << "Title=Figure 6b\n";
+ fig6b << "RatioPlot=0\n";
+ fig6b << "LogX=1\n";
+ fig6b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig6b << "YLabel=d, b\n";
+ fig6b << "Legend=1\n";
+ fig6b << "XMin=250.\n";
+ fig6b << "YMin=0.7\n";
+ fig6b << "YMax=1.05\n";
+ fig6b << "# END PLOT\n";
+ writeLine(fig6b,np,WZT,"/WZT","$W\\to Z_T$","red" ,"solid" );
+ writeLine(fig6b,np,BZT,"/BZT","$B\\to Z_T$","red" ,"dotted" );
+ writeLine(fig6b,np,ZL ,"/ZL ","$Z_L$" ,"blue" ,"dashed" );
+ writeLine(fig6b,np,H ,"/H ","$H$" ,"green","dotdashed");
+ fig6b.close();
+
+
+ np.clear();
+ vector<double> e30,e50,q30,q50,g30,g50;
+ for(Energy x=200.*GeV;x<5000.*GeV;x*=1.02) {
+ Energy2 s(sqr(x));
+ np.push_back(x);
+ evaluateLowScale(mZ,30.*GeV,s);
+ e30.push_back(real(lowColE_));
+ q30.push_back(real(lowColU_));
+ g30.push_back(real(lowColG_));
+ evaluateLowScale(mZ,50.*GeV,s);
+ e50.push_back(real(lowColE_));
+ q50.push_back(real(lowColU_));
+ g50.push_back(real(lowColG_));
+ }
+ ofstream fig4a("fig4a.dat");
+ fig4a << "#BEGIN PLOT\n";
+ fig4a << "LegendOnly=/e30 /e50\n";
+ fig4a << "DrawOnly=/e30 /e50\n";
+ fig4a << "Title=Figure 4a\n";
+ fig4a << "RatioPlot=0\n";
+ fig4a << "LogX=1\n";
+ fig4a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig4a << "YLabel=e\n";
+ fig4a << "Legend=1\n";
+ fig4a << "XMin=250.\n";
+ fig4a << "YMin=0.7\n";
+ fig4a << "YMax=1.05\n";
+ fig4a << "# END PLOT\n";
+ writeLine(fig4a,np,e30,"/e30","e $(\\mu_f=30\\,\\mathrm{GeV})$","red" ,"solid" );
+ writeLine(fig4a,np,e50,"/e50","e $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dashed");
+ fig4a.close();
+ ofstream fig4b("fig4b.dat");
+ fig4b << "#BEGIN PLOT\n";
+ fig4b << "LegendOnly=/q30 /q50 /g30 /g50\n";
+ fig4b << "DrawOnly=/q30 /q50 /g30 /g50\n";
+ fig4b << "Title=Figure 4a\n";
+ fig4b << "RatioPlot=0\n";
+ fig4b << "LogX=1\n";
+ fig4b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
+ fig4b << "YLabel=e\n";
+ fig4b << "Legend=1\n";
+ fig4b << "XMin=250.\n";
+ fig4b << "YMin=0.5\n";
+ fig4b << "YMax=1.05\n";
+ fig4b << "# END PLOT\n";
+ writeLine(fig4b,np,q30,"/q30","q $(\\mu_f=30\\,\\mathrm{GeV})$","red" ,"solid" );
+ writeLine(fig4b,np,q50,"/q50","q $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dashed");
+ writeLine(fig4b,np,g30,"/g30","g $(\\mu_f=30\\,\\mathrm{GeV})$","green" ,"dotted" );
+ writeLine(fig4b,np,g50,"/g50","g $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dotdashed");
+ fig4b.close();
+}
diff --git a/MatrixElement/EW/CollinearSudakov.fh b/MatrixElement/EW/CollinearSudakov.fh
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/CollinearSudakov.fh
@@ -0,0 +1,18 @@
+// -*- C++ -*-
+//
+// This is the forward declaration of the CollinearSudakov class.
+//
+#ifndef Herwig_CollinearSudakov_FH
+#define Herwig_CollinearSudakov_FH
+
+#include "ThePEG/Config/ThePEG.h"
+
+namespace Herwig {
+
+class CollinearSudakov;
+
+ThePEG_DECLARE_POINTERS(Herwig::CollinearSudakov,CollinearSudakovPtr);
+
+}
+
+#endif
diff --git a/MatrixElement/EW/CollinearSudakov.h b/MatrixElement/EW/CollinearSudakov.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/CollinearSudakov.h
@@ -0,0 +1,513 @@
+// -*- 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 "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:
+
+ /**
+ * 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/EWCouplings.cc b/MatrixElement/EW/EWCouplings.cc
--- a/MatrixElement/EW/EWCouplings.cc
+++ b/MatrixElement/EW/EWCouplings.cc
@@ -1,630 +1,631 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the EWCouplings class.
//
#include "EWCouplings.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <boost/numeric/ublas/operation.hpp>
using namespace Herwig;
namespace {
Complex trace(boost::numeric::ublas::matrix<Complex> M) {
assert(M.size1()==M.size2());
Complex output(0.);
for(unsigned int ix=0;ix<M.size1();++ix)
output += M(ix,ix);
return output;
}
}
EWCouplings::EWCouplings(unsigned int loops, unsigned int steps, Energy highScale,
Energy lowScale)
: ewScale_(91.1876*GeV), highScale_(highScale), lowScale_(lowScale),
includeSU3_(true), includeEW_(true), initialized_(false), massChoice_(false),
mZ_(91.1876*GeV), mW_(80.399*GeV),
mT_(173.1*GeV), // 179.08045 (should be this?)
+ mH_(125.0*GeV),
loops_(loops), highSteps_(steps), lowSteps_(steps)
{}
void EWCouplings::initialize() {
using Constants::pi;
if(initialized_) return;
initialized_ = true;
// set the particle masses
if(massChoice_) {
mZ_ = getParticleData(ParticleID::Z0 )->mass();
mW_ = getParticleData(ParticleID::Wplus)->mass();
mT_ = getParticleData(ParticleID::t )->mass();
mH_ = getParticleData(ParticleID::h0 )->mass();
}
// logs of scales
double logEWScale = log(ewScale_/GeV);
double logHighScale = log(highScale_/GeV);
// step size
double stepsize = (logHighScale - logEWScale)/double(highSteps_);
// Initialize parameters at the ewScale
// 32 parameters, mostly zero due massless quarks
unsigned int N = 32;
vector<Complex> y(N,0.), dydx(N,0.), yout(N,0.);
initializeCouplings(y);
double x = logEWScale;
derivatives(x,y,dydx);
// energy scale + 6 parameters: g1,g2,g3,y_t,lambda,vev
table_ = boost::numeric::ublas::matrix<double>(highSteps_+1,7);
table_(0,0) = logEWScale;
for (unsigned int i=1; i<N; i++) {
if (i <=3) table_(0,i ) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
if (i >=29) table_(0,i-25) = y[i].real();
}
int counter = 1;
// Use 4th order runge-kutta to integrate to highScale
int steps = highSteps_;
while (steps > 0) {
RK4(y,dydx,x,stepsize,yout);
// Advance x and calculate derivatives at new starting point
for(unsigned int j=0; j<N; j++) {
y[j] = yout[j];
}
x += stepsize;
derivatives(x,y,dydx);
table_(counter,0) = x;
for (unsigned int i=1; i<N; i++) {
if (i<=3 ) table_(counter,i ) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
if (i>=29) table_(counter,i-25) = y[i].real();
}
steps--;
counter++;
}
// Initialize couplings at mu < 91.1876 GeV
initializeLow();
}
EWCouplings::~EWCouplings() {}
IBPtr EWCouplings::clone() const {
return new_ptr(*this);
}
IBPtr EWCouplings::fullclone() const {
return new_ptr(*this);
}
void EWCouplings::persistentOutput(PersistentOStream & os) const {
os << ounit(ewScale_,GeV) << ounit(highScale_,GeV) << ounit(lowScale_,GeV)
<< includeSU3_ << includeEW_ << ounit(mZ_,GeV) << ounit(mW_,GeV)
<< ounit(mT_,GeV) << ounit(mH_,GeV) << massChoice_ << initialized_;
}
void EWCouplings::persistentInput(PersistentIStream & is, int) {
is >> iunit(ewScale_,GeV) >> iunit(highScale_,GeV) >> iunit(lowScale_,GeV)
>> includeSU3_ >> includeEW_ >> iunit(mZ_,GeV) >> iunit(mW_,GeV)
>> iunit(mT_,GeV) >> iunit(mH_,GeV) >> massChoice_ >> initialized_;
}
//The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<EWCouplings,Interfaced>
describeHerwigEWCouplings("Herwig::EWCouplings", "HwMEEW.so");
void EWCouplings::Init() {
static ClassDocumentation<EWCouplings> documentation
("The EWCouplings class implements");
static Switch<EWCouplings,bool> interfaceMassChoice
("MassChoice",
"Where to get the SM particle masses from",
&EWCouplings::massChoice_, false, false, false);
static SwitchOption interfaceMassChoiceLocal
(interfaceMassChoice,
"Local",
"Use local values",
false);
static SwitchOption interfaceMassChoiceParticleData
(interfaceMassChoice,
"ParticleData",
"Get the values from the ParticleData object",
true);
}
void EWCouplings::initializeLow() {
using Constants::pi;
// For scales less than ewScale, the only couplings calculated here are those of
// alpha_EW and alpha3:
double logEWScale = log(ewScale_ /GeV);
double loglowScale = log(lowScale_/GeV);
double stepsize = (loglowScale - logEWScale)/double(lowSteps_);
int steps = lowSteps_;
// Initialize parameters at the ewScale
unsigned int N=2; // Total # of parameters = 2
vector<Complex> y(N), dydx(N), yout(N);
for (unsigned int i=0; i<N; i++) {
y[i] = 0.0;
dydx[i] = 0.0;
yout[i] = 0.0;
}
double x = logEWScale;
// Initialize Couplings at the ewScale, including the One-Loop Threshold Matching:
double a1 = (3.0/5.0)*table_(0,1);
double a2 = table_(0,2);
double a3 = table_(0,3);
double aEM_ewScale = 1./(1./a1+1./a2-1./(6.*pi)*(1.-21.*log(mW()/ewScale_)+16./3.*log(mTatmZ()/ewScale_)));
double aS_ewScale = 1./(1./a3+1./(3.*pi)*log(ewScale_/mTatmZ()));
y[0] = sqrt(4.0*pi*aEM_ewScale);
y[1] = sqrt(4.0*pi*aS_ewScale);
derivatives(x,y,dydx);
// energy scale + 2 parameters: gEM,g3
lowTable_.resize(lowSteps_+1,3);
lowTable_(0,0) = logEWScale;
for (unsigned int i=1; i<=N; i++) {
lowTable_(0,i) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
}
int counter = 1;
// Use 4th order runge-kutta to integrate to highScale
while (steps > 0) {
RK4(y,dydx,x,stepsize,yout);
// Advance x and calculate derivatives at new starting point
for(unsigned int j=0; j<N; j++) {
y[j] = yout[j];
}
x += stepsize;
derivatives(x,y,dydx);
lowTable_(counter,0) = x;
for (unsigned int i=1; i<=N; i++) {
lowTable_(counter,i) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
}
steps--;
counter++;
}
}
void EWCouplings::RK4(vector<Complex> &y, vector<Complex> &dydx,
const double x, const double h, vector<Complex> &yout) {
unsigned int n = y.size();
std::vector<Complex> dym(n), dyt(n), yt(n);
double hh = h*0.5;
double h6 = h/6.0;
double xh = x + hh;
const Complex I(0,1.0);
for(unsigned int i=0; i<n; i++) yt[i] = y[i] + hh*dydx[i];
derivatives(xh, yt, dyt);
for(unsigned int i=0; i<n; i++) yt[i] = y[i] + hh*dyt[i];
derivatives(xh, yt, dym);
for(unsigned int i=0; i<n; i++) {
yt[i] = y[i] + h*dym[i];
dym[i] += dyt[i];
}
derivatives(x+h, yt, dyt);
for(unsigned int i=0; i<n; i++) {
yout[i] = y[i] + h6*(dydx[i] + dyt[i] + 2.0*dym[i]);
}
}
void EWCouplings::initializeCouplings(vector<Complex> & y) {
// \todo make these values parameters so they can be changed
InvEnergy2 gFermi = 1.16637*pow(10.0,-5)/GeV2;
Energy vev = 1.0/(sqrt(sqrt(2.0)*gFermi)); // vev = 246.221
y[0] = 0.461531463; // g1 = Sqrt[5/3] * Sqrt[4*pi*a1] with a1(Mz) = 0.01017054
y[1] = 0.651547066; // g2 = Sqrt[4*pi*a2] with a2(Mz) = 0.03378168
y[2] = 1.215650108; // g3 = Sqrt[4*pi*as] with as(Mz) = 0.1176
// Note lambda_t = sqrt(2.0)*mt/vev only valid for mt(mt); need mt(mZ) here
// Top Yukawa lambda from Manohar
//Complex lambda_t = 1.02858;
// Top Yukawa lambda from Sascha
double lambda_t = 0.991172;
// Quartic coupling lambda (need to multiply by a factor of 2 when accessing the quartic coupling)
double lambda = (mH_/vev)*(mH_/vev);
y[29] = lambda_t;
y[30] = lambda;
y[31] = vev/GeV;
}
void EWCouplings::derivatives(const double x, vector<Complex> & y,
vector<Complex> &dydx) {
// zero the output
for (unsigned int i=0; i<dydx.size(); i++) dydx[i]=0.0;
// low scale
if (y.size()==2) {
lowBetaGauge(x,y,dydx);
}
// high scale
else {
betaGauge(x,y,dydx);
betaYukawa(x,y,dydx);
betaHiggs(x,y,dydx);
}
}
void EWCouplings::betaGauge(const double x, vector<Complex> &y, vector<Complex> & dydx) {
using Constants::pi;
using boost::numeric::ublas::axpy_prod;
using boost::numeric::ublas::herm;
const Complex I(0,1.0);
// Yukawa
boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
Yuk_u(ix,iy) = y[21+3*ix+iy];
Yuk_d(ix,iy) = y[12+3*ix+iy];
Yuk_e(ix,iy) = y[ 3+3*ix+iy];
}
}
// gauge
boost::numeric::ublas::vector<Complex> gauge(3);
for(unsigned int l=0; l<3; l++) gauge[l] = y[l];
// Evaluate beta functions for gauge couplings
double Ng = 0.5*numberOfFlavours(x);
boost::numeric::ublas::vector<Complex> b(3), g2(3), gc(3), Cu(3), Cd(3), Ce(3);
boost::numeric::ublas::matrix<Complex> B1(3,3),B2(3,3), B3(3,3), B(3,3);
b[0] = -4.0/3.0*Ng - 1.0/10.0;
b[1] = 22.0/3.0 - 4.0/3.0*Ng - 1.0/6.0;
b[2] = 11.0 - 4.0/3.0*Ng;
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
B1(ix,iy) = 0.;
B2(ix,iy) = 0.;
B3(ix,iy) = 0.;
}
}
B1(1,1) = 136.0/3.0;
B1(2,2) = 102.0;
B2(0,0) = 19.0/15.0;
B2(0,1) = 1.0/5.0;
B2(0,2) = 11.0/30.0;
B2(1,0) = 3.0/5.0;
B2(1,1) = 49.0/3.0;
B2(1,2) = 3.0/2.0 ;
B2(2,0) = 44.0/15.0;
B2(2,1) = 4.0;
B2(2,2) = 76.0/3.0;
B3(0,0) = 9.0/50.0;
B3(0,1) = 3.0/10.0;
B3(1,0) = 9.0/10.0;
B3(1,1) = 13.0/6.0;
B = B1 - Ng*B2 - B3;
Cu[0] = 17.0/10.0;
Cu[1] = 3.0/2.0;
Cu[2] = 2.0;
Cd[0] = 1.0/2.0;
Cd[1] = 3.0/2.0;
Cd[2] = 2.0;
Ce[0] = 3.0/2.0;
Ce[1] = 1.0/2.0;
Ce[2] = 0.0;
for (int i=0; i<3; i++) {
g2(i) = pow(gauge(i),2);
}
// gc = trans(g2) * B
axpy_prod(B, g2, gc, true);
// compute the answer
if(loops_ >= 1) {
for(int l=0; l<3; l++) {
dydx[l] += -b(l)*pow(gauge(l),3)/(16.0*pow(pi,2));
}
if (loops_ >= 2) {
boost::numeric::ublas::matrix<Complex> temp(3,3);
axpy_prod(herm(Yuk_u),Yuk_u,temp);
Complex tr1 = trace(temp);
axpy_prod(herm(Yuk_d),Yuk_d,temp);
Complex tr2 = trace(temp);
axpy_prod(herm(Yuk_e),Yuk_e,temp);
Complex tr3 = trace(temp);
for(int l=0; l<3; l++) {
dydx[l] += -pow(gauge(l),3)/pow(16.0*pow(pi,2),2)*
(gc(l) + Cu(l)*tr1 + Cd(l)*tr2 + Ce(l)*tr3 );
}
}
}
}
void EWCouplings::lowBetaGauge(const double, vector<Complex> &y,
vector<Complex> &dydx) {
using Constants::pi;
const Complex I(0,1.0);
Complex e = y[0], g3 = y[1];
// Evaluate beta functions for gauge couplings
double Nu = 2.0, Nd = 3.0, Nl = 3.0;
if(loops_ >=1) {
dydx[0] += (16.0/9.0*Nu + 4.0/9.0*Nd + 4.0/3.0*Nl)*pow(e,3)/pow(4.0*pi,2);
dydx[1] += (2.0/3.0*(Nu+Nd)-11.0)*pow(g3,3)/pow(4.0*pi,2);
// Note this also includes the three-loop contribution for alpha_3
if (loops_ >= 2) {
dydx[0] += (64.0/27.0*Nu+4.0/27.0*Nd+4.0*Nl)*pow(e,5)/pow(4.0*pi,4) +
(64.0/9.0*Nu+16.0/9.0*Nd)*pow(e,3)*pow(g3,2)/pow(4.0*pi,4);
dydx[1] += (38.0/3.0*(Nu+Nd)-102.0)*pow(g3,5)/pow(4.0*pi,4) +
(8.0/9.0*Nu+2.0/9.0*Nd)*pow(g3,3)*pow(e,2)/pow(4.0*pi,4) +
(5033.0/18.0*(Nu+Nd)-325.0/54.0*(Nu+Nd)*(Nu+Nd)-2857.0/2.0)*pow(g3,7)/pow(4.0*pi,6);
}
}
}
void EWCouplings::betaYukawa(const double x, vector< Complex > &y, vector<Complex > &dydx) {
using Constants::pi;
const Complex I(0,1.0);
boost::numeric::ublas::identity_matrix<Complex> Id(3,3);
// Yukawa
boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
Yuk_u(ix,iy) = y[21+3*ix+iy];
Yuk_d(ix,iy) = y[12+3*ix+iy];
Yuk_e(ix,iy) = y[ 3+3*ix+iy];
}
}
// gauge
double Ng = 0.5*numberOfFlavours(x);
boost::numeric::ublas::vector<Complex> gauge(3);
for(unsigned int l=0; l<3; l++) gauge[l] = y[l];
Complex lambda = y[30];
// traces of yukawa matrices
boost::numeric::ublas::matrix<Complex> mTemp(3,3),MUU(3,3),MDD(3,3),MLL(3,3),
MUU2(3,3),MDD2(3,3),MLL2(3,3),MUUDD(3,3),MDDUU(3,3);
axpy_prod(herm(Yuk_u),Yuk_u,MUU);
Complex trU = trace( MUU);
axpy_prod(MUU,MUU,MUU2);
Complex trUU = trace(MUU2);
axpy_prod(herm(Yuk_d),Yuk_d,MDD);
Complex trD = trace( MUU);
axpy_prod(MDD,MDD,MDD2);
Complex trDD = trace(MDD2);
axpy_prod(MUU,MDD,MUUDD);
Complex trUD = trace(MUUDD);
axpy_prod(MDD,MUU,MDDUU);
axpy_prod(herm(Yuk_e),Yuk_e,MLL);
Complex trL = trace( MLL);
axpy_prod(MLL,MLL,MLL2);
Complex trLL = trace(MLL2);
Complex g02 = sqr(gauge[0]);
Complex g12 = sqr(gauge[1]);
Complex g22 = sqr(gauge[2]);
// Evaluate beta functions for yukawa couplings
boost::numeric::ublas::zero_matrix<Complex> zero3x3(3,3);
boost::numeric::ublas::matrix<Complex> dYuk_udx(zero3x3), dYuk_ddx(zero3x3), dYuk_edx(zero3x3), beta1_u(zero3x3),
beta1_d(zero3x3), beta1_e(zero3x3), beta2_u(zero3x3), beta2_d(zero3x3), beta2_e(zero3x3);
Complex Y2 = 3.0*trU+3.0*trD + trL;
Complex Y4 = (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trU +
(1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trD + 3.0/4.0*(g02 + g12)*trL;
Complex chi4 = 27.0/4.0*trUU + 27.0/4.0*trDD + 9.0/4.0*trLL - 6.0/4.0*trUD;
if(loops_ >= 1) {
beta1_u = 3.0/2.0*(MUU - MDD) + (Y2 - 17.0/20.0*g02 - 9.0/4.0*g12 - 8.0*g22)*Id;
beta1_d = 3.0/2.0*(MDD - MUU) + (Y2 - 1.0/4.0*g02 - 9.0/4.0*g12 - 8.0*g22)*Id;
beta1_e = 3.0/2.0*(MLL) + (Y2 - 9.0/4.0*g02 - 9.0/4.0*g12)*Id;
axpy_prod(Yuk_u,beta1_u,mTemp);
dYuk_udx += (1.0/(16.0*pow(pi,2)))*mTemp;
axpy_prod(Yuk_d,beta1_d,mTemp);
dYuk_ddx += (1.0/(16.0*pow(pi,2)))*mTemp;
axpy_prod(Yuk_e,beta1_e,mTemp);
dYuk_edx += (1.0/(16.0*pow(pi,2)))*mTemp;
if (loops_ >= 2) {
Complex l2=sqr(lambda);
beta2_u = 3.0/2.0*MUU2 - MUUDD - 1.0/4.0*MDDUU + 11.0/4.0*MDD2 + Y2*(5.0/4.0*MDD - 9.0/4.0*MUU) +
(-chi4 + 3.0/2.0*l2)*Id - 2.0*lambda*(3.0*MUU + MDD) +
(223.0/80.0*g02 + 135.0/16.0*g12 + 16.0*g22)*(MUU) -
(43.0/80.0*g02 - 9.0/16.0*g12 + 16.0*g22)*(MDD) + 5.0/2.0*Y4*Id +
((9.0/200.0 + 29.0/45.0*Ng)*pow(gauge[0],4) - 9.0/20.0*pow(gauge[0]*gauge[1],2) +
19.0/15.0*pow(gauge[0]*gauge[2],2) - (35.0/4.0 - Ng)*pow(gauge[1],4) + 9.0*pow(gauge[1]*gauge[2],2) -
(404.0/3.0 - 80.0/9.0*Ng)*pow(gauge[2],4))*Id;
beta2_d = 3.0/2.0*MDD2 - MDDUU - 1.0/4.0*MUUDD + 11.0/4.0*MUU2 + Y2*(5.0/4.0*MUU - 9.0/4.0*MDD) + (-chi4 + 3.0/2.0*l2)*Id - 2.0*lambda*(3.0*MDD + MUU) + (187.0/80.0*g02 + 135.0/16.0*g12 + 16.0*g22)*(MDD) - (79.0/80.0*g02 - 9.0/16.0*g12 + 16.0*g22)*(MUU) + 5.0/2.0*Y4*Id - ((29.0/200.0 + 1.0/45.0*Ng)*pow(gauge[0],4) - 27.0/20.0*pow(gauge[0]*gauge[1],2) + 31.0/15.0*pow(gauge[0]*gauge[2],2) - (35.0/4.0 - Ng)*pow(gauge[1],4) + 9.0*pow(gauge[1]*gauge[2],2) - (404.0/3.0 - 80.0/9.0*Ng)*pow(gauge[2],4))*Id;
beta2_e = 3.0/2.0*MLL2 - 9.0/4.0*Y2*MLL + (-chi4 + 3.0/2.0*l2)*Id - 6.0*lambda*(MLL) + (387.0/80.0*g02 + 135.0/15.0*g12)*(MLL) + 5.0/2.0*Y4*Id + ((51.0/200.0 + 11.0/5.0*Ng)*pow(gauge[0],4) + 27.0/20.0*pow(gauge[0]*gauge[1],2) - (35.0/4.0 - Ng)*pow(gauge[1],4))*Id;
axpy_prod(Yuk_u,beta2_u,mTemp);
dYuk_udx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
axpy_prod(Yuk_d,beta2_d,mTemp);
dYuk_ddx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
axpy_prod(Yuk_e,beta2_e,mTemp);
dYuk_edx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
}
}
boost::numeric::ublas::vector<Complex> temp(27);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
dydx[ 3+ix+3*iy] = dYuk_edx(ix,iy);
dydx[12+ix+3*iy] = dYuk_ddx(ix,iy);
dydx[21+ix+3*iy] = dYuk_udx(ix,iy);
}
}
}
void EWCouplings::betaHiggs(const double x, vector<Complex> &y,
vector<Complex> &dydx) {
using Constants::pi;
const Complex I(0,1.0);
double Ng = 0.5*numberOfFlavours(x);
// Yukawa
boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
Yuk_u(ix,iy) = y[21+3*ix+iy];
Yuk_d(ix,iy) = y[12+3*ix+iy];
Yuk_e(ix,iy) = y[ 3+3*ix+iy];
}
}
// gauge
boost::numeric::ublas::vector<Complex> gauge(3);
for(int l=0; l<3; l++) gauge(l) = y[l];
Complex lambda = y[30];
complex<Energy> vev = y[31]*GeV;
// Evaluate beta functions for higgs coupling
Complex beta1_lambda(0.), beta2_lambda(0.), gamma1_vev(0.), gamma2_vev(0.),
Y2(0.), H(0.), Y4(0.), chi4(0.);
// traces of yukawa matrices
boost::numeric::ublas::matrix<Complex> temp(3,3),temp2(3,3),MUU(3,3),MDD(3,3),MLL(3,3);
axpy_prod(herm(Yuk_u),Yuk_u,MUU);
Complex trU = trace( MUU);
axpy_prod(MUU,MUU,temp);
Complex trUU = trace(temp);
axpy_prod(MUU,temp,temp2);
Complex trUUU = trace(temp2);
axpy_prod(herm(Yuk_d),Yuk_d,MDD);
Complex trD = trace( MUU);
axpy_prod(MDD,MDD,temp);
Complex trDD = trace(temp);
axpy_prod(MDD,temp,temp2);
Complex trDDD = trace(temp2);
axpy_prod(MUU,MDD,temp);
Complex trUD = trace(temp);
axpy_prod(herm(Yuk_e),Yuk_e,MLL);
Complex trL = trace( MLL);
axpy_prod(MLL,MLL,temp);
Complex trLL = trace(temp);
axpy_prod(MLL,temp,temp2);
Complex trLLL = trace(temp2);
axpy_prod(MUU+MDD,MDD,temp);
axpy_prod(MUU,temp,temp2);
Complex trUUDD = trace(temp2);
Complex g02 = sqr(gauge[0]);
Complex g12 = sqr(gauge[1]);
Complex g22 = sqr(gauge[2]);
Complex g04 = sqr(g02);
Complex g14 = sqr(g12);
double pi2 = sqr(pi);
Y2 = 3.0*trU+3.0*trD + trL;
Y4 = (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*(trU) + (1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*(trD) + 3.0/4.0*(g02 + g12)*(trL);
chi4 = 27.0/4.0*trUU + 27.0/4.0*trDD + 9.0/4.0*trLL - 6.0/4.0*trUD;
H = 3.0*trUU + 3.0*trDD + trLL;
if(loops_ >= 1) {
Complex l2=sqr(lambda);
beta1_lambda = 12.0*l2 - (9.0/5.0*g02 + 9.0*g12)*lambda + 9.0/4.0*(3.0/25.0*g04+2.0/5.0*g02*g12 + g14) + 4.0*Y2*lambda - 4.0*H;
gamma1_vev = 9.0/4.0*(1.0/5.0*g02+g12)-Y2;
dydx[30] += 1.0/(16.0*pi2)*beta1_lambda;
dydx[31] += vev/(16.0*pi2)*gamma1_vev/GeV;
if (loops_ >= 2) {
beta2_lambda = -78.0*lambda*l2 + 18.0*(3.0/5.0*g02 + 3.0*g12)*l2 -
( (313.0/8.0 - 10.0*Ng)*g14 - 117.0/20.0*g02*g12 + 9.0/25.0*(229.0/4.0+50.0/9.0*Ng)*g04 )*lambda +
(497.0/8.0 - 8.0*Ng)*g14*g12 - 3.0/5.0*(97.0/24.0 + 8.0/3.0*Ng)*g02*g14 -
9.0/25.0*(239.0/24.0 + 40.0/9.0*Ng)*g04*g12 - 27.0/125.0*(59.0/24.0 + 40.0/9.0*Ng)*g04*g02 -
64.0*g22*(trUU + trDD) - 8.0/5.0*g02*(2.0*trUU - trDD + 3.0*trLL) - 3.0/2.0*g14*Y4 +
10.0*lambda*( (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trU + (1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trD + 3.0/4.0*(g02 + g12)*trL ) +
3.0/5.0*g02*( (-57.0/10.0*g02 + 21.0*g12 )*trU + (3.0/2.0*g02 + 9.0*g12)*trD + (-15.0/2.0*g02 + 11.0*g12)*trL ) - 24.0*l2*Y2 - lambda*H +
6.0*lambda*trUD + 20.0*(3.0*trUUU + 3.0*trDDD + trLLL) - 12.0*trUUDD;
gamma2_vev = -3.0/2.0*l2 - 5.0/2.0*Y4 + chi4 - 27.0/80.0*g02*g12 -
(93.0/800.0 + 1.0/2.0*Ng)*g04 + (511.0/32.0 - 5.0/2.0*Ng)*g14;
dydx[30] += 1.0/pow(16.0*pi2,2)*beta2_lambda;
dydx[31] += vev/pow(16.0*pi2,2)*gamma2_vev/GeV;
}
}
}
double EWCouplings::interpolate(double t, int paramIndex) {
double stepsize = table_(1,0)-table_(0,0);
double tol = 0.001*stepsize;
double logewScale = log(ewScale_/GeV);
double loghighScale = log(highScale_/GeV);
if (t<logewScale-tol || t>loghighScale+tol) {
cerr << "Stepsize: " << stepsize << std::endl;
cerr << "tol: " << tol << std::endl;
cerr << "logewScale: " << logewScale << std::endl;
cerr << "loghighScale: " << loghighScale << std::endl;
cerr << "paramIndex: " << paramIndex << std::endl;
cerr << "t: " << t << std::endl;
cerr << "Couplings::_Interp(double t, int parmIndex) trying to obtain parameter ";
cerr << "value outside the available range. Returning zero." << std::endl;
assert(false);
}
// return value at EW scale
if (abs(t-logewScale)<tol) return table_(0,paramIndex);
// return value at high scale
if (abs(t-loghighScale)<tol) {
return table_(table_.size1()-1,paramIndex);
}
unsigned int numSteps = int((t-table_(0,0))/stepsize);
// Linear Interpolation:
double x1 = table_(numSteps,0);
double y1 = table_(numSteps,paramIndex);
double x2 = table_(numSteps+1,0);
double y2 = table_(numSteps+1,paramIndex);
return y1+((y2-y1)/(x2-x1))*(t-x1);
}
double EWCouplings::interpolateLow(double t, int paramIndex) {
double stepsize = lowTable_(0,0)-lowTable_(1,0);
double tol = 0.00001*stepsize;
double logewScale = log(ewScale_ /GeV);
double loglowScale = log(lowScale_/GeV);
if (t<loglowScale-tol || t>logewScale+tol) {
cerr<< "Couplings::_LowInterp(double t, int parmIndex) trying to obtain parameter ";
cerr << "value outside the available range. Returning zero." << std::endl;
assert(false);
}
if (abs(t-logewScale)<tol) {
return lowTable_(0,paramIndex);
}
if (std::abs(t-loglowScale)<tol) {
return lowTable_(lowTable_.size1()-1,paramIndex);
}
int numSteps = (int)((lowTable_(0,0)-t)/stepsize);
// Linear Interpolation:
double x1 = lowTable_(numSteps,0);
double y1 = lowTable_(numSteps,paramIndex);
double x2 = lowTable_(numSteps+1,0);
double y2 = lowTable_(numSteps+1,paramIndex);
return y1+((y2-y1)/(x2-x1))*(t-x1);
}
diff --git a/MatrixElement/EW/ElectroWeakReweighter.cc b/MatrixElement/EW/ElectroWeakReweighter.cc
--- a/MatrixElement/EW/ElectroWeakReweighter.cc
+++ b/MatrixElement/EW/ElectroWeakReweighter.cc
@@ -1,95 +1,103 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ElectroWeakReweighter class.
//
#include "ElectroWeakReweighter.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
+tEWCouplingsPtr ElectroWeakReweighter::staticEWCouplings_ = tEWCouplingsPtr();
+
+
ElectroWeakReweighter::ElectroWeakReweighter() {}
ElectroWeakReweighter::~ElectroWeakReweighter() {}
IBPtr ElectroWeakReweighter::clone() const {
return new_ptr(*this);
}
IBPtr ElectroWeakReweighter::fullclone() const {
return new_ptr(*this);
}
void ElectroWeakReweighter::persistentOutput(PersistentOStream & os) const {
- os << EWCouplings_;
+ os << EWCouplings_ << collinearSudakov_;
}
void ElectroWeakReweighter::persistentInput(PersistentIStream & is, int) {
- is >> EWCouplings_;
+ is >> EWCouplings_ >> collinearSudakov_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<ElectroWeakReweighter,ReweightBase>
describeHerwigElectroWeakReweighter("Herwig::ElectroWeakReweighter", "HwMEEW.so");
void ElectroWeakReweighter::Init() {
static ClassDocumentation<ElectroWeakReweighter> documentation
("There is no documentation for the ElectroWeakReweighter class");
static Reference<ElectroWeakReweighter,EWCouplings> interfaceEWCouplings
("EWCouplings",
"The object to calculate the electroweak couplings",
&ElectroWeakReweighter::EWCouplings_, false, false, true, false, false);
+ static Reference<ElectroWeakReweighter,CollinearSudakov> interfaceCollinearSudakov
+ ("CollinearSudakov",
+ "The collinear Sudakov",
+ &ElectroWeakReweighter::collinearSudakov_, false, false, true, false, false);
+
}
double ElectroWeakReweighter::weight() const {
EWCouplings_->initialize();
staticEWCouplings_ = EWCouplings_;
cerr << "aEM\n";
for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.1) {
cerr << scale/GeV << " "
<< EWCouplings_->aEM(scale) << "\n";
}
cerr << "aS\n";
for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
cerr << scale/GeV << " "
<< EWCouplings_->aS(scale) << "\n";
}
cerr << "y_t\n";
for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
cerr << scale/GeV << " "
<< EWCouplings_->y_t(scale) << "\n";
}
cerr << "lambda\n";
for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
cerr << scale/GeV << " "
<< EWCouplings_->lambda(scale) << "\n";
}
cerr << "vev\n";
for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
cerr << scale/GeV << " "
<< EWCouplings_->vev(scale)/GeV << "\n";
}
+ collinearSudakov_->makePlots();
-
- cerr << "testing got here ???\n";
- cerr << subProcess() << "\n";
- cerr << *subProcess() << "\n";
- cerr << subProcess()->outgoing()[0] << *subProcess()->outgoing()[0] << "\n";
- cerr << subProcess()->outgoing()[0]->spinInfo() << "\n";
- cerr << subProcess()->outgoing()[0]->spinInfo()->productionVertex() << "\n";
+ // cerr << "testing got here ???\n";
+ // cerr << subProcess() << "\n";
+ // cerr << *subProcess() << "\n";
+ // cerr << subProcess()->outgoing()[0] << *subProcess()->outgoing()[0] << "\n";
+ // cerr << subProcess()->outgoing()[0]->spinInfo() << "\n";
+ // cerr << subProcess()->outgoing()[0]->spinInfo()->productionVertex() << "\n";
assert(false);
staticEWCouplings_ = tEWCouplingsPtr();
}
diff --git a/MatrixElement/EW/ElectroWeakReweighter.h b/MatrixElement/EW/ElectroWeakReweighter.h
--- a/MatrixElement/EW/ElectroWeakReweighter.h
+++ b/MatrixElement/EW/ElectroWeakReweighter.h
@@ -1,121 +1,127 @@
// -*- C++ -*-
#ifndef Herwig_ElectroWeakReweighter_H
#define Herwig_ElectroWeakReweighter_H
//
// This is the declaration of the ElectroWeakReweighter class.
//
#include "ThePEG/MatrixElement/ReweightBase.h"
#include "EWCouplings.h"
+#include "CollinearSudakov.h"
namespace Herwig {
using namespace ThePEG;
/**
* The ElectroWeakReweighter class.
*
* @see \ref ElectroWeakReweighterInterfaces "The interfaces"
* defined for ElectroWeakReweighter.
*/
class ElectroWeakReweighter: public ReweightBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ElectroWeakReweighter();
/**
* The destructor.
*/
virtual ~ElectroWeakReweighter();
//@}
public:
/**
* Return the weight for the kinematical configuation provided by
* the assigned XComb object (in the LastXCombInfo base class).
*/
virtual double weight() const;
/**
*
*/
static tEWCouplingsPtr coupling() {
assert(staticEWCouplings_);
return staticEWCouplings_;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ElectroWeakReweighter & operator=(const ElectroWeakReweighter &);
private:
/**
* The Electroweak Couplings
*/
EWCouplingsPtr EWCouplings_;
/**
+ * The Collinear Sudakov
+ */
+ CollinearSudakovPtr collinearSudakov_;
+
+ /**
* The couplings to allow global access
*/
static tEWCouplingsPtr staticEWCouplings_;
};
}
#endif /* Herwig_ElectroWeakReweighter_H */
diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
--- a/MatrixElement/EW/GroupInvariants.h
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -1,234 +1,254 @@
// -*- 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>
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$
*/
- double C_A(unsigned int N) {
+ inline double C_A(unsigned int N) {
return N !=1 ? double(N) : 0.;
}
/**
* The \f$SU(N)\f$ \f$C_F\f$
*/
- double C_F(unsigned int N) {
+ inline double C_F(unsigned int N) {
return N !=1 ? 0.5*(double(N*N)-1.)/double(N) : 1.;
}
/*
* The \f$SU(N)\f$ \f$C_d\f$
*/
- double C_d(unsigned int N) {
+ inline double C_d(unsigned int N) {
return (double(N*N)-4.)/double(N);
}
/**
* The \f$SU(N)\f$\f$C_1\f$
*/
- double C_1(unsigned int N) {
+ inline double C_1(unsigned int N) {
double N2(N*N);
return 0.25*(N2-1.0)/N2;
}
/**
* \f$T_F\f$
*/
- double T_F(unsigned int N, bool high) {
+ 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$
*/
- double t_S(unsigned int, bool ) {
+ inline double t_S(unsigned int, bool ) {
return 0.5;
}
/**
* / Number of complex scalars in the fundamental rep. of SU(N)/U(1)
*/
- double n_S(unsigned int N, bool high) {
+ 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)
*/
- double n_F(unsigned int N, bool high) {
+ 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);
+ }
}
// namespace DiracHigh {
// double n_g() { // Number of fermion generations (only used in gauge boson HighCMatching)
// return 3.0;
// }
// }
// namespace WeylHigh {
// double n_S(int N) { // Number of complex scalars in the fundamental rep. of SU(N)/U(1)
// if(N==2 || N==1) return 1.0;
// if(N==3) return 0.0;
// std::cout << "Error! SU(N), N != (1, 2 or 3) used for n_S in ";
// std::cout << "GroupInvariants.h but not defined." << std::endl;
// return 0.0;
// }
// double n_F(int N) { // Number of Weyl Fermions in the fundamental rep. of SU(N)
// if(N==2) return 12.0;
// if(N==3) return 12.0;
// std::cout << "Error! SU(N), N != (2 or 3) used for n_F in ";
// std::cout << "GroupInvariants.h but not defined." << std::endl;
// return 0.0;
// }
// double n_g() { // Number of fermion generations (only used in gauge boson HighCMatching)
// return 3.0;
// }
// double T_F(int N) { return 0.5+0.0*N; } // I believe T(N) is such that tr(t^a t^b) = T delta(ab)
// // 0.0*N included to stop receiving a stupid warning.
// double t_S(int N) { return 0.5+0.0*N; } // Analog of T_F but for scalars.
// // 0.0*N included to stop receiving a stupid warning.
// }
// namespace WeylLow {
// double n_S(int N) { // Number of complex scalars in the fundamental rep. of SU(N)
// if(N==1) return 0.0;
// if(N==3) return 0.0;
// std::cout << "Error! SU(N), N != (1 or 3) used for n_S in ";
// std::cout << "GroupInvariants.h but not defined." << std::endl;
// return 0.0;
// }
// double n_F(int N) { // Number of Weyl Fermions in the fundamental rep. of SU(N)
// if(N==3) return 10.0;
// if(N==1) return 2.0;
// std::cout << "Error! SU(N), N != (1 or 3) used for n_F in ";
// std::cout << "GroupInvariants.h but not defined." << std::endl;
// return 0.0;
// }
// double T_F(int N) { return 0.5+0.0*N; } // I believe T(N) is such that tr(t^a t^b) = T delta(ab)
// // 0.0*N included to stop receiving a stupid warning.
// double t_S(int N) { return 0.5+0.0*N; } // Analog of T_F but for scalars.
// // 0.0*N included to stop receiving a stupid warning.
// }
// #endif // GROUP_INVARIANTS_H
}
#endif // HERWIG_GroupInvariants_H
diff --git a/MatrixElement/EW/Makefile.am b/MatrixElement/EW/Makefile.am
--- a/MatrixElement/EW/Makefile.am
+++ b/MatrixElement/EW/Makefile.am
@@ -1,5 +1,6 @@
pkglib_LTLIBRARIES = HwMEEW.la
HwMEEW_la_SOURCES = GroupInvariants.h GroupInvariants.cc \
ElectroWeakReweigter.h ElectroWeakReweighter.cc \
+CollinearSudakov.h CollinearSudakov.cc \
EWCouplings.h EWCouplings.fh EWCouplings.cc
HwMEEW_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 4:40 PM (21 h, 41 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023487
Default Alt Text
(71 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment