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 +describeHerwigCollinearSudakov("Herwig::CollinearSudakov", "HwMEEW.so"); + +void CollinearSudakov::Init() { + + static ClassDocumentation 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 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 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 W1 = WFC.fFW0-0.5*WFC.aW0-0.5*WFC.cW0; + complex W2 = WFC.fF0W-0.5*WFC.a0W; + complex 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 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 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 HtR = HtL-0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.atW0); + complex HbL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.at0W); + + complex 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 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 F_bL = HbL + aW/(4.0*pi)*0.5*W2; + + Complex Dw = CollinearDw(s,EWScale); + Complex Dz = CollinearDz(s,EWScale); + + complex D_C_UL = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*Dz + 0.5*Dw; + complex D_C_DL = ElectroWeakReweighter::coupling()->g_Ld(EWScale)*ElectroWeakReweighter::coupling()->g_Ld(EWScale)*Dz + 0.5*Dw; + + complex D_C_UR = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*Dz; + complex D_C_DR = ElectroWeakReweighter::coupling()->g_Rd(EWScale)*ElectroWeakReweighter::coupling()->g_Rd(EWScale)*Dz; + + complex D_C_nuL = ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*Dz + 0.5*Dw; + complex D_C_EL = ElectroWeakReweighter::coupling()->g_Le(EWScale)*ElectroWeakReweighter::coupling()->g_Le(EWScale)*Dz + 0.5*Dw; + complex D_C_ER = ElectroWeakReweighter::coupling()->g_Re(EWScale)*ElectroWeakReweighter::coupling()->g_Re(EWScale)*Dz; + + complex 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 D_C_WZ = aW/(4.0*pi)*2.0*(F_W+WFC.fsZW1) + 0.5*WFC.RZ + sqrt(sW2/cW2)*WFC.RAtoZ; + complex D_C_WA = aW/(4.0*pi)*2.0*(F_W+WFC.fs01) + 0.5*WFC.RA + sqrt(cW2/sW2)*WFC.RZtoA; + complex D_C_BZ = 0.5*WFC.RZ - sqrt(cW2/sW2)*WFC.RAtoZ; + complex 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 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 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 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 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 x, vector 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 np; + vector 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 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 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 using namespace Herwig; namespace { Complex trace(boost::numeric::ublas::matrix M) { assert(M.size1()==M.size2()); Complex output(0.); for(unsigned int ix=0;ixmass(); 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 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(highSteps_+1,7); table_(0,0) = logEWScale; for (unsigned int i=1; 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=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 describeHerwigEWCouplings("Herwig::EWCouplings", "HwMEEW.so"); void EWCouplings::Init() { static ClassDocumentation documentation ("The EWCouplings class implements"); static Switch 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 y(N), dydx(N), yout(N); for (unsigned int i=0; i 0) { RK4(y,dydx,x,stepsize,yout); // Advance x and calculate derivatives at new starting point for(unsigned int j=0; j &y, vector &dydx, const double x, const double h, vector &yout) { unsigned int n = y.size(); std::vector 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 & 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 & y, vector &dydx) { // zero the output for (unsigned int i=0; i &y, vector & 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 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 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 b(3), g2(3), gc(3), Cu(3), Cd(3), Ce(3); boost::numeric::ublas::matrix 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 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 &y, vector &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 &dydx) { using Constants::pi; const Complex I(0,1.0); boost::numeric::ublas::identity_matrix Id(3,3); // Yukawa boost::numeric::ublas::matrix 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 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 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 zero3x3(3,3); boost::numeric::ublas::matrix 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 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 &y, vector &dydx) { using Constants::pi; const Complex I(0,1.0); double Ng = 0.5*numberOfFlavours(x); // Yukawa boost::numeric::ublas::matrix 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 gauge(3); for(int l=0; l<3; l++) gauge(l) = y[l]; Complex lambda = y[30]; complex 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 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 (tloghighScale+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)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)> EWCouplings_; + is >> EWCouplings_ >> collinearSudakov_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigElectroWeakReweighter("Herwig::ElectroWeakReweighter", "HwMEEW.so"); void ElectroWeakReweighter::Init() { static ClassDocumentation documentation ("There is no documentation for the ElectroWeakReweighter class"); static Reference interfaceEWCouplings ("EWCouplings", "The object to calculate the electroweak couplings", &ElectroWeakReweighter::EWCouplings_, false, false, true, false, false); + static Reference 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 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 mu0.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