Page MenuHomeHEPForge

No OneTemporary

diff --git a/MatrixElement/EW/CollinearSudakov.cc b/MatrixElement/EW/CollinearSudakov.cc
--- a/MatrixElement/EW/CollinearSudakov.cc
+++ b/MatrixElement/EW/CollinearSudakov.cc
@@ -1,1192 +1,1730 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the CollinearSudakov class.
//
#include "CollinearSudakov.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "GroupInvariants.h"
#include "ElectroWeakReweighter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
namespace {
void DuplicateColumn0(boost::numeric::ublas::matrix<Complex> &orig) {
for (unsigned int i=0; i<orig.size1(); i++) {
for (unsigned int j=1; j<orig.size2(); j++) {
orig(i,j) = orig(i,0);
}
}
}
}
CollinearSudakov::CollinearSudakov() : K_ORDER_(3),
B_ORDER_(2)
{}
CollinearSudakov::~CollinearSudakov() {}
IBPtr CollinearSudakov::clone() const {
return new_ptr(*this);
}
IBPtr CollinearSudakov::fullclone() const {
return new_ptr(*this);
}
void CollinearSudakov::persistentOutput(PersistentOStream & os) const {
os << K_ORDER_ << B_ORDER_;
}
void CollinearSudakov::persistentInput(PersistentIStream & is, int) {
is >> K_ORDER_ >> B_ORDER_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<CollinearSudakov,Interfaced>
describeHerwigCollinearSudakov("Herwig::CollinearSudakov", "HwMEEW.so");
void CollinearSudakov::Init() {
static ClassDocumentation<CollinearSudakov> documentation
("The CollinearSudakov class implements the collinear evolution");
}
Complex CollinearSudakov::highScaleIntegral(bool SU3, bool SU2, double Y,
Energy2 s, Energy mu_h, Energy mu_l,
bool fermion, bool longitudinal,
double yukFactor) {
SU3_ = SU3;
SU2_ = SU2;
Y_ = Y;
s_ = s;
fermion_ = fermion;
longitudinal_ = longitudinal;
yukFactor_ = yukFactor;
// perform the integral
Complex result;
high_ = true;
// real piece
real_ = true;
result.real(integrator_.value(*this,mu_h,mu_l));
// imaginary piece
real_ = false;
result.imag(integrator_.value(*this,mu_h,mu_l));
// return the answer
return exp(result);
}
Complex CollinearSudakov::lowScaleIntegral(bool SU3, double Q, Energy2 s,
Energy mu_h, Energy mu_l, bool fermion,
double boostFactor) {
SU3_ = SU3;
Q_ = Q;
s_ = s;
fermion_ = fermion;
boostFactor_ = boostFactor;
// perform the integral
Complex result;
high_ = false;
// real piece
real_ = true;
result.real(integrator_.value(*this,mu_h,mu_l));
// imaginary piece
real_ = false;
result.imag(integrator_.value(*this,mu_h,mu_l));
// return the answer
return exp(result);
}
InvEnergy CollinearSudakov::highScaleIntegrand(Energy mu) const {
using namespace GroupInvariants;
using Constants::pi;
Complex gamma = 0.0;
// Include K-factor Contributions (Cusps):
GaugeContributions cusp = cuspContributions(mu,K_ORDER_,high_);
// Include B-factors (B1 for U1, B2 for SU2, B3 for SU3):
GaugeContributions nonCusp = BContributions(mu,B_ORDER_,fermion_,longitudinal_);
// common log
Complex plog = PlusLog(s_/sqr(mu));
// SU(3)
if(SU3_) {
if(fermion_)
gamma += C_F(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
else
gamma += (C_A(3))*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
}
// SU(2)
if(SU2_) {
if (fermion_ || longitudinal_ )
gamma += (C_F(2))*cusp.SU2*0.5*plog + 0.5*nonCusp.SU2;
else
gamma += (C_A(2))*cusp.SU2*0.5*plog + 0.5*nonCusp.SU2;
}
if (fermion_ || longitudinal_ ) {
gamma += sqr(Y_)*(cusp.U1*0.5*plog + 0.5*nonCusp.U1);
}
else {
// U(1) Gauge boson
if (!SU3_ && !SU2_ && abs(Y_)<0.001) {
gamma += 0.5*nonCusp.U1;
}
}
// top Yukawa piece
double y_t = ElectroWeakReweighter::coupling()->y_t(mu);
gamma += yukFactor_*sqr(y_t)/(16.0*sqr(pi));
// return the answer
return real_ ? gamma.real()/mu : gamma.imag()/mu;
}
InvEnergy CollinearSudakov::lowScaleIntegrand(Energy mu) const {
using namespace GroupInvariants;
using Constants::pi;
Complex gamma = 0.0;
// Include K-factor Contributions (Cusps):
GaugeContributions cusp = cuspContributions(mu,K_ORDER_,false);
// Include B-factors (B1 for U1, B2 for SU2, B3 for SU3):
GaugeContributions nonCusp = BContributionsLow(mu,B_ORDER_,fermion_,boostFactor_);
// common log
Complex plog = PlusLog(s_/sqr(mu));
Complex blog(0.);
if(boostFactor_ >0.001) blog = PlusLog(4.0*boostFactor_);
// SU(3)
if (SU3_) {
if (fermion_) {
// not a bHQET top quark field
if (abs(boostFactor_)<0.001)
gamma += C_F(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
else
gamma += C_F(3)*cusp.SU3*0.5*blog + 0.5*nonCusp.SU3;
}
else {
gamma += C_A(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3;
}
}
// fermions
if (fermion_) {
// not a bHQET top quark field
if (abs(boostFactor_)<0.001)
gamma += sqr(Q_)*(cusp.U1*0.5*plog + 0.5*nonCusp.U1);
else
gamma += sqr(Q_)*(cusp.U1*0.5*blog + 0.5*nonCusp.U1);
}
else {
// i.e., not a fermion, not a bHQ, not a gluon => photon
if (abs(boostFactor_)<0.001 && !SU3_)
gamma += 0.5*nonCusp.U1;
// i.e., W treated as a bHQET field
else if (abs(boostFactor_)>0.001)
gamma += sqr(Q_)*(cusp.U1*0.5*blog + 0.5*nonCusp.U1);
}
// return the answer
return real_ ? gamma.real()/mu : gamma.imag()/mu;
}
void CollinearSudakov::evaluateHighScale(Energy highScale, Energy EWScale, Energy2 S) {
double yCoeffQt(0.), yCoefftR(0.), yCoeffPhi(0.);
if (K_ORDER_ >= 1) {
yCoeffQt = 0.5;
yCoefftR = 1.0;
yCoeffPhi = 3.0;
}
highColW_ = highScaleIntegral(false,true ,0.0 ,S,highScale,EWScale,false,false,0.0);
highColB_ = highScaleIntegral(false,false,0.0 ,S,highScale,EWScale,false,false,0.0);
highColG_ = highScaleIntegral(true ,false,0.0 ,S,highScale,EWScale,false,false,0.0);
highColQ_ = highScaleIntegral(true ,true ,1./6. ,S,highScale,EWScale,true,false,0.0);
highColQt_ = highScaleIntegral(true ,true ,1./6. ,S,highScale,EWScale,true,false,yCoeffQt);
highColU_ = highScaleIntegral(true ,false,2./3. ,S,highScale,EWScale,true,false,0.0);
highColtR_ = highScaleIntegral(true ,false,2./3. ,S,highScale,EWScale,true,false,yCoefftR);
highColD_ = highScaleIntegral(true ,false,-1./3.,S,highScale,EWScale,true,false,0.0);
highColL_ = highScaleIntegral(false,true ,-1./2.,S,highScale,EWScale,true,false,0.0);
highColE_ = highScaleIntegral(false,false,-1. ,S,highScale,EWScale,true,false,0.0);
highColPhi_ = highScaleIntegral(false,true ,1./2. ,S,highScale,EWScale,false,true,yCoeffPhi);
}
void CollinearSudakov::evaluateLowScale(Energy EWScale, Energy lowScale, Energy2 S) {
lowColW_ = lowScaleIntegral(false,1. ,S,EWScale,lowScale,false,
S/sqr(2.*ElectroWeakReweighter::coupling()->mW()));
lowColA_ = lowScaleIntegral(false,0. ,S,EWScale,lowScale,false,0.0);
lowColG_ = lowScaleIntegral(true ,0. ,S,EWScale,lowScale,false,0.0);
lowColU_ = lowScaleIntegral(true ,2./3. ,S,EWScale,lowScale,true,0.0);
lowColt_ = lowScaleIntegral(true ,2./3. ,S,EWScale,lowScale,true,
S/sqr(2.*ElectroWeakReweighter::coupling()->mT()));
lowColD_ = lowScaleIntegral(true ,-1./3.,S,EWScale,lowScale,true,0.0);
lowColE_ = lowScaleIntegral(false,-1. ,S,EWScale,lowScale,true,0.0);
}
void CollinearSudakov::evaluateMatching(Energy EWScale,Energy2 s) {
using Constants::pi;
using GroupInvariants::PlusLog;
static const Complex I(0,1.0);
// wave function corrections
WaveFunctionCorrections WFC = waveFunctionCorrections(EWScale);
double aS = ElectroWeakReweighter::coupling()->a3(EWScale);
double aEM = ElectroWeakReweighter::coupling()->aEM(EWScale);
double aW = ElectroWeakReweighter::coupling()->aW(EWScale);
double aZ = ElectroWeakReweighter::coupling()->aZ(EWScale);
double cW2 = ElectroWeakReweighter::coupling()->Cos2thW(EWScale);
double sW2 = ElectroWeakReweighter::coupling()->Sin2thW(EWScale);
Energy mW = ElectroWeakReweighter::coupling()->mW();
Energy mZ = ElectroWeakReweighter::coupling()->mZ();
Energy mH = ElectroWeakReweighter::coupling()->mH();
Energy mT = ElectroWeakReweighter::coupling()->mT();
double gPHI = ElectroWeakReweighter::coupling()->g_phiPlus(EWScale);
double y_t = ElectroWeakReweighter::coupling()->y_t(EWScale);
double lz = log(mZ/EWScale);
double lw = log(mW/EWScale);
complex<double> F_W = 4.0*lw*0.5*PlusLog(s/EWScale/EWScale)-2.0*lw*lw-2.0*lw-5.0*pi*pi/12.0+1.0;
complex<double> F_Z = 4.0*lz*0.5*PlusLog(s/EWScale/EWScale)-2.0*lz*lz-2.0*lz-5.0*pi*pi/12.0+1.0;
// Taken from Manohar... along with his formulae for F_tL, F_tR, and F_bL (for 3rd/1st Gen. Wavefunction Differences)
complex<double> W1 = WFC.fFW0-0.5*WFC.aW0-0.5*WFC.cW0;
complex<double> W2 = WFC.fF0W-0.5*WFC.a0W;
complex<double> U1 = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*(WFC.fFZZ-0.5*WFC.aZZ) -
0.5*WFC.cZZ*(ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale) +
ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)) +
ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*WFC.bZZ;
complex<double> U2 = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*(WFC.fFZZ-0.5*WFC.aZZ) -
0.5*WFC.cZZ*(ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale) +
ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)) +
ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*WFC.bZZ;
complex<double> HtL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-0.5*log(mH/EWScale)-0.5*lz+
0.5*WFC.atHH+0.5*WFC.atZZ+WFC.ctHH+WFC.ctZZ+
WFC.ctW0-WFC.btHH+WFC.btZZ);
complex<double> HtR = HtL-0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.atW0);
complex<double> HbL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.at0W);
complex<double> F_tL = (4.0/3.0*aS+4.0/9.0*aEM)/(4.0*pi)*(2.0*log(mT/EWScale)*log(mT/EWScale)-log(mT/EWScale)+
pi*pi/12.0+2.0) + HtL +
aW/(4.0*pi)*0.5*W1 + aZ/(4.0*pi)*U1;
complex<double> F_tR = (4.0/3.0*aS+4.0/9.0*aEM)/(4.0*pi)*(2.0*log(mT/EWScale)*log(mT/EWScale)-log(mT/EWScale)+
pi*pi/12.0+2.0) + HtR -
aW/(4.0*pi)*0.25*WFC.cW0 + aZ/(4.0*pi)*U2;
complex<double> F_bL = HbL + aW/(4.0*pi)*0.5*W2;
Complex Dw = CollinearDw(s,EWScale);
Complex Dz = CollinearDz(s,EWScale);
complex<double> D_C_UL = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*Dz + 0.5*Dw;
complex<double> D_C_DL = ElectroWeakReweighter::coupling()->g_Ld(EWScale)*ElectroWeakReweighter::coupling()->g_Ld(EWScale)*Dz + 0.5*Dw;
complex<double> D_C_UR = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*Dz;
complex<double> D_C_DR = ElectroWeakReweighter::coupling()->g_Rd(EWScale)*ElectroWeakReweighter::coupling()->g_Rd(EWScale)*Dz;
complex<double> D_C_nuL = ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*Dz + 0.5*Dw;
complex<double> D_C_EL = ElectroWeakReweighter::coupling()->g_Le(EWScale)*ElectroWeakReweighter::coupling()->g_Le(EWScale)*Dz + 0.5*Dw;
complex<double> D_C_ER = ElectroWeakReweighter::coupling()->g_Re(EWScale)*ElectroWeakReweighter::coupling()->g_Re(EWScale)*Dz;
complex<double> D_C_WW = aW/(4.0*pi)*cW2*(F_Z+WFC.fsWZWZ) +
aW/(4.0*pi)*cW2*(F_W+WFC.fs1ZW) + aW/(4.0*pi)*sW2*(F_W+WFC.fs10) +
aW/(4.0*pi)*sW2*(2.0*lw*lw-
2.0*lw+pi*pi/12.0+2.0) + 0.5*WFC.RW;
complex<double> D_C_WZ = aW/(4.0*pi)*2.0*(F_W+WFC.fsZW1) + 0.5*WFC.RZ + sqrt(sW2/cW2)*WFC.RAtoZ;
complex<double> D_C_WA = aW/(4.0*pi)*2.0*(F_W+WFC.fs01) + 0.5*WFC.RA + sqrt(cW2/sW2)*WFC.RZtoA;
complex<double> D_C_BZ = 0.5*WFC.RZ - sqrt(cW2/sW2)*WFC.RAtoZ;
complex<double> D_C_BA = 0.5*WFC.RA - sqrt(sW2/cW2)*WFC.RZtoA;
// The WFC.RW and WFC.RZ are used on purpose (instead of WFC.RPhi and WFC.RPhi3, respectively):
complex<double> D_C_PhiW = aW/(4.0*pi)*0.25*(F_W+WFC.fs1HW) +
aW/(4.0*pi)*0.25*(F_W+WFC.fs1ZW) + aZ/(4.0*pi)*gPHI*gPHI*(F_Z+WFC.fsWZWZ) +
aW/(4.0*pi)*sW2*(2.0*lw*lw-2.0*lw+pi*pi/12.0+2.0) +
0.5*WFC.RW + WFC.EW;
complex<double> D_C_PhiZ = aW/(4.0*pi)*0.5*(F_W+WFC.fsZW1) + aZ/(4.0*pi)*0.25*(F_Z+WFC.fs1HZ) + 0.5*WFC.RZ + WFC.EZ;
complex<double> D_C_PhiH = aW/(4.0*pi)*0.5*(F_W+WFC.fsHW1) + aZ/(4.0*pi)*0.25*(F_Z+WFC.fsHZ1) + 0.5*WFC.RH;
complex<double> D_C_GG = aS/(4.0*pi)*2.0/3.0*log(mT/EWScale);
ULcollinearCorr_ = exp(D_C_UL);
DLcollinearCorr_ = exp(D_C_DL);
URcollinearCorr_ = exp(D_C_UR);
DRcollinearCorr_ = exp(D_C_DR);
tLcollinearCorr_ = exp(D_C_UL+F_tL);
tRcollinearCorr_ = exp(D_C_UR+F_tR);
bLcollinearCorr_ = exp(D_C_DL+F_bL);
nuLcollinearCorr_ = exp(D_C_nuL);
ELcollinearCorr_ = exp(D_C_EL);
ERcollinearCorr_ = exp(D_C_ER);
WtoWcollinearCorr_ = exp(D_C_WW);
WtoZcollinearCorr_ = exp(D_C_WZ);
WtoAcollinearCorr_ = exp(D_C_WA);
BtoZcollinearCorr_ = exp(D_C_BZ);
BtoAcollinearCorr_ = exp(D_C_BA);
PhitoWcollinearCorr_ = exp(D_C_PhiW);
PhitoZcollinearCorr_ = exp(D_C_PhiZ);
PhitoHcollinearCorr_ = exp(D_C_PhiH);
GcollinearCorr_ = exp(D_C_GG);
}
WaveFunctionCorrections CollinearSudakov::waveFunctionCorrections(Energy EWScale) {
static const Complex I(0.,1.);
using Constants::pi;
double lZ = 2.0*log(ElectroWeakReweighter::coupling()->mZ()/EWScale);
WaveFunctionCorrections WFC;
// From Manohar, 2/12/12: (these assume mH=125, mZ=91.1876, mW=80.399)
WFC.RW = (0.8410283377963967 - 9.424777961271568*I) + 1.2785863646210789*lZ;
WFC.RA = (1.4835982362022198 + 1.855775680704845*pow(10.,-9)*I) - 0.27209907467584415*lZ;
WFC.RZ = (1.5114724841549798 - 9.926944419863688*I) + 1.0834802397165764*lZ;
WFC.RAtoZ = (0.3667485032811715 - 2.2770907130064835*I) - 1.2994544609942593*lZ;
WFC.RZtoA = -0.2095310079712942 + 0.8320191107808546*lZ;
WFC.RH = (12.229832449946716 - 1.7643103462419842*10.0*pow(10.,-12)*I) + 5.309998583664737*lZ;
WFC.RPhi = (5.569012418081201 + 1.5439133581417356*0.10*pow(10.,-9)*I) + 5.309998583664737*lZ;
WFC.RPhi3 = (8.945333042265943 + 5.499309445612249*pow(10.,-12)*I) + 5.309998583664737*lZ;
WFC.EW = (3.967645734304811 + 4.712388980384717*I) + 2.238332625165702*lZ;
WFC.EZ = (5.916079892937651 + 4.96347220970469*I) + 2.1132591719740788*lZ;
WFC.RW *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.RA *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.RZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.RAtoZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.RZtoA *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.RPhi *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.EW *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.EZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.RPhi3 *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.RH *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi);
WFC.RW = WFC.RW.real();
WFC.RA = WFC.RA.real();
WFC.RZ = WFC.RZ.real();
WFC.RAtoZ = WFC.RAtoZ.real();
WFC.RZtoA = WFC.RZtoA.real();
WFC.RPhi = WFC.RPhi.real();
WFC.RPhi3 = WFC.RPhi3.real();
WFC.RH = WFC.RH.real();
// Also from Manohar, 2/10/12
WFC.fFW0 = -3.7946842553189453 - 4.709019671388589*I;
WFC.fF0W = 3.8181790485161176;
WFC.fFZZ = 1.364503250377989 + 0.*I;
WFC.aHH = -0.474396977740686 + 0.*I;
WFC.aZZ = -0.6806563210877914 + 0.*I;
WFC.aW0 = 0.49036102506811907 + 1.9323351450971642*I;
WFC.a0W = -1.2184776671065072;
WFC.bHH = 1.234775745474991 + 0.*I;
WFC.bZZ = 1.7303526426747613 + 0.*I;
WFC.cHH = 0.33140608274387473 + 0.*I;
WFC.cZZ = 0.4961363208382017 + 0.*I;
WFC.cW0 = -1.005299829777395 + 1.063048500002757*I;
WFC.atHH = -0.237198488870343 + 0.*I;
WFC.atZZ = -0.3403281605438957 + 0.*I;
WFC.atW0 = 0.24518051253405954 + 0.9661675725485821*I;
WFC.at0W = -0.6092388335532536;
WFC.ctHH = 0.16570304137193737 + 0.*I;
WFC.ctZZ = 0.24806816041910085 + 0.*I;
WFC.ctW0 = -0.5026499148886975 + 0.5315242500013785*I;
WFC.btHH = -0.30869393636874776 + 0.*I;
WFC.btZZ = -0.4325881606686903 + 0.*I;
WFC.fs10 = -2.289868133696459;
WFC.fs1ZW = 1.8627978596240622;
WFC.fsWZWZ = 1.1866922529667008;
WFC.fsZW1 = 1.0840307611156266;
WFC.fs01 = 2.2898681336964595;
WFC.fsHW1 = -0.32306745562682404;
WFC.fsHZ1 = 0.4042992116255279;
WFC.fs1HW = 3.330954543719127;
WFC.fs1HZ = 2.69768201932412;
return WFC;
}
Complex CollinearSudakov::CollinearDw(Energy2 s, Energy EWScale) {
using Constants::pi;
using GroupInvariants::PlusLog;
double lw = log(ElectroWeakReweighter::coupling()->mW()/EWScale);
//s = s/2.; // This should not be here... I think this is a discrepency with Sascha
return ElectroWeakReweighter::coupling()->aW(EWScale)/(4.0*pi)*
(4.0*lw*0.5*PlusLog(s/EWScale/EWScale) - 2.0*lw*lw -
3.0*lw - 5.0/12.0*pi*pi + 9.0/4.0);
}
Complex CollinearSudakov::CollinearDz(Energy2 s, Energy EWScale) {
using Constants::pi;
using GroupInvariants::PlusLog;
double lz = log(ElectroWeakReweighter::coupling()->mZ()/EWScale);
return ElectroWeakReweighter::coupling()->aZ(EWScale)/(4.0*pi)*
(4.0*lz*0.5*PlusLog(s/EWScale/EWScale) - 2.0*lz*lz -
3.0*lz - 5.0/12.0*pi*pi + 9.0/4.0);
}
namespace {
void writeLine(ofstream & file, vector<Energy> x, vector<double> y,
string name,string title,
string colour, string style) {
file << "# BEGIN HISTO1D "+name +"\n";
file << "SmoothLine=1\n";
file << "ErrorBars=0\n";
file << "Path="+name+"\n";
file << "Title="+title+"\n";
file << "LineColor="+colour+"\n";
file << "LineStyle="+style +"\n";
file << "# xlow xhigh val errminus errplus\n";
for(unsigned int ix=0;ix<x.size();++ix)
file << (x[ix]-MeV)/GeV << "\t" << (x[ix]+MeV)/GeV << "\t"
<< y[ix] << "\t" << 0. << "\t" << 0. << "\n";
file << "# END HISTO1D\n";
}
}
void CollinearSudakov::makePlots() {
vector<Energy> np;
vector<double> uL,uR,tL,tR,dL,dR,bL,bR,nuL,eL,eR,Wgamma,Bgamma,g,WT,WL,WZT,BZT,ZL,H;
Energy mZ = ElectroWeakReweighter::coupling()->mZ();
for(Energy x=200.*GeV;x<5000.*GeV;x*=1.02) {
Energy2 s(sqr(x));
np.push_back(x);
evaluateHighScale(x,mZ,s);
evaluateMatching(mZ,s);
uL .push_back(real(highColQ_ * ULcollinearCorr_ ));
uR .push_back(real(highColU_ * URcollinearCorr_ ));
tL .push_back(real(highColQt_ * tLcollinearCorr_ ));
tR .push_back(real(highColtR_ * tRcollinearCorr_ ));
dL .push_back(real(highColQ_ * DLcollinearCorr_ ));
dR .push_back(real(highColD_ * DRcollinearCorr_ ));
bL .push_back(real(highColQt_ * bLcollinearCorr_ ));
bR .push_back(real(highColD_ * DRcollinearCorr_ ));
nuL .push_back(real(highColL_ * nuLcollinearCorr_ ));
eL .push_back(real(highColL_ * ELcollinearCorr_ ));
eR .push_back(real(highColE_ * ERcollinearCorr_ ));
Wgamma.push_back(real(highColW_ * WtoAcollinearCorr_ ));
Bgamma.push_back(real(highColB_ * BtoAcollinearCorr_ ));
g .push_back(real(highColG_ * GcollinearCorr_ ));
WT .push_back(real(highColW_ * WtoWcollinearCorr_ ));
WL .push_back(real(highColPhi_ * PhitoWcollinearCorr_));
WZT .push_back(real(highColW_ * WtoZcollinearCorr_ ));
BZT .push_back(real(highColB_ * BtoZcollinearCorr_ ));
ZL .push_back(real(highColPhi_ * PhitoZcollinearCorr_));
H .push_back(real(highColPhi_ * PhitoHcollinearCorr_));
}
ofstream fig1a("fig1a.dat");
fig1a << "#BEGIN PLOT\n";
fig1a << "LegendOnly=/uL /uR /tL /tR\n";
fig1a << "DrawOnly=/uL /uR /tL /tR\n";
fig1a << "Title=Figure 1a\n";
fig1a << "RatioPlot=0\n";
fig1a << "LogX=1\n";
fig1a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig1a << "YLabel=u, t\n";
fig1a << "Legend=1\n";
fig1a << "XMin=250.\n";
fig1a << "YMin=0.7\n";
fig1a << "YMax=1.05\n";
fig1a << "# END PLOT\n";
writeLine(fig1a,np,uL,"/uL","$u_L$","green","dotted" );
writeLine(fig1a,np,uR,"/uR","$u_R$","cyan" ,"solid" );
writeLine(fig1a,np,tL,"/tL","$t_L$","red" ,"dashed" );
writeLine(fig1a,np,tR,"/tR","$t_R$","blue" ,"dotdashed");
fig1a.close();
ofstream fig1b("fig1b.dat");
fig1b << "#BEGIN PLOT\n";
fig1b << "LegendOnly=/dL /dR /bL /bR\n";
fig1b << "DrawOnly=/dL /dR /bL /bR\n";
fig1b << "Title=Figure 1b\n";
fig1b << "RatioPlot=0\n";
fig1b << "LogX=1\n";
fig1b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig1b << "YLabel=d, b\n";
fig1b << "Legend=1\n";
fig1b << "XMin=250.\n";
fig1b << "YMin=0.7\n";
fig1b << "YMax=1.05\n";
fig1b << "# END PLOT\n";
writeLine(fig1b,np,dL,"/dL","$d_L$","green","dotted" );
writeLine(fig1b,np,dR,"/dR","$d_R$","cyan" ,"solid" );
writeLine(fig1b,np,bL,"/bL","$b_L$","red" ,"dashed" );
writeLine(fig1b,np,bR,"/bR","$b_R$","blue" ,"dotdashed");
fig1b.close();
ofstream fig2("fig2.dat");
fig2 << "#BEGIN PLOT\n";
fig2 << "LegendOnly=/uL /uR /dL /dR\n";
fig2 << "DrawOnly=/uL /uR /dL /dR\n";
fig2 << "Title=Figure 2\n";
fig2 << "RatioPlot=0\n";
fig2 << "LogX=1\n";
fig2 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig2 << "YLabel=u, d\n";
fig2 << "Legend=1\n";
fig2 << "XMin=250.\n";
fig2 << "YMin=0.7\n";
fig2 << "YMax=1.05\n";
fig2 << "# END PLOT\n";
writeLine(fig2,np,uL,"/uL","$u_L$","green","dotted" );
writeLine(fig2,np,uR,"/uR","$u_R$","cyan" ,"solid" );
writeLine(fig2,np,dL,"/dL","$d_L$","red" ,"dashed" );
writeLine(fig2,np,dR,"/dR","$d_R$","blue" ,"dotdashed");
fig2.close();
ofstream fig3("fig3.dat");
fig3 << "#BEGIN PLOT\n";
fig3 << "LegendOnly=/nuL /eL /eR\n";
fig3 << "DrawOnly=/nuL /eL /eR\n";
fig3 << "Title=Figure 3\n";
fig3 << "RatioPlot=0\n";
fig3 << "LogX=1\n";
fig3 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig3 << "YLabel=$\\nu$, $e$\n";
fig3 << "Legend=1\n";
fig3 << "XMin=250.\n";
fig3 << "YMin=0.9\n";
fig3 << "YMax=1.05\n";
fig3 << "# END PLOT\n";
writeLine(fig3,np,nuL,"/nuL","$\\nu_L$","blue","dashed");
writeLine(fig3,np, eL,"/eL" ,"$e_L$" ,"red" ,"dotted");
writeLine(fig3,np, eR,"/eR" ,"$e_R$" ,"red" ,"solid" );
fig3.close();
ofstream fig5("fig5.dat");
fig5 << "#BEGIN PLOT\n";
fig5 << "LegendOnly=/g /Wgamma /Bgamma\n";
fig5 << "DrawOnly=/g /Wgamma /Bgamma\n";
fig5 << "Title=Figure 5\n";
fig5 << "RatioPlot=0\n";
fig5 << "LogX=1\n";
fig5 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig5 << "YLabel=$\\gamma$, g\n";
fig5 << "Legend=1\n";
fig5 << "XMin=250.\n";
fig5 << "YMin=0.7\n";
fig5 << "YMax=1.05\n";
fig5 << "# END PLOT\n";
writeLine(fig5,np,g ,"/g" ,"$g$" ,"blue","dashed");
writeLine(fig5,np,Wgamma,"/Wgamma","$W\\to\\gamma$","red" ,"solid" );
writeLine(fig5,np,Bgamma,"/Bgamma","$B\\to\\gamma$","red" ,"dotted");
fig5.close();
ofstream fig6a("fig6a.dat");
fig6a << "#BEGIN PLOT\n";
fig6a << "LegendOnly=/WT /WL\n";
fig6a << "DrawOnly=/WT /WL\n";
fig6a << "Title=Figure 6a\n";
fig6a << "RatioPlot=0\n";
fig6a << "LogX=1\n";
fig6a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig6a << "YLabel=$Z_L$, $Z_T$, $H$\n";
fig6a << "Legend=1\n";
fig6a << "XMin=250.\n";
fig6a << "YMin=0.7\n";
fig6a << "YMax=1.05\n";
fig6a << "# END PLOT\n";
writeLine(fig6a,np,WT,"/WT","$W_T$","red" ,"solid");
writeLine(fig6a,np,WL,"/WL","$W_L$","blue" ,"dashed" );
fig6a.close();
ofstream fig6b("fig6b.dat");
fig6b << "#BEGIN PLOT\n";
fig6b << "LegendOnly=/WZT /BZT /ZL /H\n";
fig6b << "DrawOnly=/WZT /BZT /ZL /H\n";
fig6b << "Title=Figure 6b\n";
fig6b << "RatioPlot=0\n";
fig6b << "LogX=1\n";
fig6b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig6b << "YLabel=d, b\n";
fig6b << "Legend=1\n";
fig6b << "XMin=250.\n";
fig6b << "YMin=0.7\n";
fig6b << "YMax=1.05\n";
fig6b << "# END PLOT\n";
writeLine(fig6b,np,WZT,"/WZT","$W\\to Z_T$","red" ,"solid" );
writeLine(fig6b,np,BZT,"/BZT","$B\\to Z_T$","red" ,"dotted" );
writeLine(fig6b,np,ZL ,"/ZL ","$Z_L$" ,"blue" ,"dashed" );
writeLine(fig6b,np,H ,"/H ","$H$" ,"green","dotdashed");
fig6b.close();
np.clear();
vector<double> e30,e50,q30,q50,g30,g50;
for(Energy x=200.*GeV;x<5000.*GeV;x*=1.02) {
Energy2 s(sqr(x));
np.push_back(x);
evaluateLowScale(mZ,30.*GeV,s);
e30.push_back(real(lowColE_));
q30.push_back(real(lowColU_));
g30.push_back(real(lowColG_));
evaluateLowScale(mZ,50.*GeV,s);
e50.push_back(real(lowColE_));
q50.push_back(real(lowColU_));
g50.push_back(real(lowColG_));
}
ofstream fig4a("fig4a.dat");
fig4a << "#BEGIN PLOT\n";
fig4a << "LegendOnly=/e30 /e50\n";
fig4a << "DrawOnly=/e30 /e50\n";
fig4a << "Title=Figure 4a\n";
fig4a << "RatioPlot=0\n";
fig4a << "LogX=1\n";
fig4a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig4a << "YLabel=e\n";
fig4a << "Legend=1\n";
fig4a << "XMin=250.\n";
fig4a << "YMin=0.7\n";
fig4a << "YMax=1.05\n";
fig4a << "# END PLOT\n";
writeLine(fig4a,np,e30,"/e30","e $(\\mu_f=30\\,\\mathrm{GeV})$","red" ,"solid" );
writeLine(fig4a,np,e50,"/e50","e $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dashed");
fig4a.close();
ofstream fig4b("fig4b.dat");
fig4b << "#BEGIN PLOT\n";
fig4b << "LegendOnly=/q30 /q50 /g30 /g50\n";
fig4b << "DrawOnly=/q30 /q50 /g30 /g50\n";
fig4b << "Title=Figure 4a\n";
fig4b << "RatioPlot=0\n";
fig4b << "LogX=1\n";
fig4b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n";
fig4b << "YLabel=e\n";
fig4b << "Legend=1\n";
fig4b << "XMin=250.\n";
fig4b << "YMin=0.5\n";
fig4b << "YMax=1.05\n";
fig4b << "# END PLOT\n";
writeLine(fig4b,np,q30,"/q30","q $(\\mu_f=30\\,\\mathrm{GeV})$","red" ,"solid" );
writeLine(fig4b,np,q50,"/q50","q $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dashed");
writeLine(fig4b,np,g30,"/g30","g $(\\mu_f=30\\,\\mathrm{GeV})$","green" ,"dotted" );
writeLine(fig4b,np,g50,"/g50","g $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dotdashed");
fig4b.close();
}
boost::numeric::ublas::matrix<Complex>
CollinearSudakov::electroWeakMatching(Energy EWScale, Energy2 s,
Herwig::EWProcess::Process process,
bool oneLoop) {
using namespace EWProcess;
// calculate the matching coefficients
evaluateMatching(EWScale,s);
// fill the matrix
boost::numeric::ublas::matrix<Complex> result(1,1);
switch (process) {
case QQQQ:
case QQQQiden:
{
unsigned int numGauge = 4, numBrokenGauge = 12;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
result(0,0) = result(6,0) = ULcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_;
result(3,0) = result(9,0) = DLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_;
for (int i=0; i<12; i++) {
if (i!=0 && i!=3 && i!=6 && i!=9) {
result(i,0) = ULcollinearCorr_*ULcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_;
}
}
DuplicateColumn0(result);
}
break;
case QtQtQQ:
{
unsigned int numGauge = 4, numBrokenGauge = 12;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
result(0,0) = result(6,0) = ULcollinearCorr_*ULcollinearCorr_*tLcollinearCorr_*tLcollinearCorr_;
result(3,0) = result(9,0) = DLcollinearCorr_*DLcollinearCorr_*bLcollinearCorr_*bLcollinearCorr_;
for (int i=0; i<12; i++) {
if (i==4 || i==5 || i==10 || i==11) {
result(i,0) = ULcollinearCorr_*tLcollinearCorr_*DLcollinearCorr_*bLcollinearCorr_;
}
else if (i==1 || i==7) {
result(i,0) = DLcollinearCorr_*DLcollinearCorr_*tLcollinearCorr_*tLcollinearCorr_;
}
else if (i==2 || i==8) {
result(i,0) = ULcollinearCorr_*ULcollinearCorr_*bLcollinearCorr_*bLcollinearCorr_;
}
}
DuplicateColumn0(result);
}
break;
case QQUU:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QtQtUU:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
result(1,0) = result(3,0) = bLcollinearCorr_*bLcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QQtRtR:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*tRcollinearCorr_*tRcollinearCorr_;
result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*tRcollinearCorr_*tRcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QQDD:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QtQtDD:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
result(1,0) = result(3,0) = bLcollinearCorr_*bLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QQLL:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_;
result(1,0) = nuLcollinearCorr_*nuLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_;
result(2,0) = ELcollinearCorr_*ELcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_;
result(3,0) = ELcollinearCorr_*ELcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_;
result(4,0) = result(5,0) = nuLcollinearCorr_*ELcollinearCorr_*ULcollinearCorr_*DLcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QQEE:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = ULcollinearCorr_*ULcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
result(1,0) = DLcollinearCorr_*DLcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
DuplicateColumn0(result);
}
break;
case UUUU:
case UUUUiden:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = URcollinearCorr_*URcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
DuplicateColumn0(result);
}
break;
case tRtRUU:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = tRcollinearCorr_*tRcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
DuplicateColumn0(result);
}
break;
case UUDD:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = URcollinearCorr_*URcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
DuplicateColumn0(result);
}
break;
case tRtRDD:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = tRcollinearCorr_*tRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
DuplicateColumn0(result);
}
break;
case UULL:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
result(1,0) = ELcollinearCorr_*ELcollinearCorr_*URcollinearCorr_*URcollinearCorr_;
DuplicateColumn0(result);
}
break;
case UUEE:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 1;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = URcollinearCorr_*URcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
DuplicateColumn0(result);
}
break;
case DDDD:
case DDDDiden:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = DRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
DuplicateColumn0(result);
}
break;
case DDLL:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
result(1,0) = ELcollinearCorr_*ELcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_;
DuplicateColumn0(result);
}
break;
case DDEE:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 1;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = DRcollinearCorr_*DRcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
DuplicateColumn0(result);
}
break;
case LLLL:
case LLLLiden:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_;
result(1,0) = nuLcollinearCorr_*nuLcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_;
result(2,0) = ELcollinearCorr_*ELcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_;
result(3,0) = ELcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_;
result(4,0) = result(5,0) = nuLcollinearCorr_*ELcollinearCorr_*nuLcollinearCorr_*ELcollinearCorr_;
DuplicateColumn0(result);
}
break;
case LLEE:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
result(1,0) = ELcollinearCorr_*ELcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
DuplicateColumn0(result);
}
break;
case EEEE:
case EEEEiden:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 1;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = ERcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QQWW:
case LLWW:
{
unsigned int numGauge = 5;
unsigned int numBrokenGauge = 20;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
for (unsigned int row = 0; row < result.size1(); row++) {
for (unsigned int col = 0; col < result.size2(); col++) {
// Boson Collinear Corr_ections:
if (col==0 || col==1) {
if (row==0 || row==1 || row==6 || row==7) result(row,col) = (WtoWcollinearCorr_*WtoWcollinearCorr_);
if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*WtoZcollinearCorr_);
if (row==3 || row==4 || row==9 || row==10) result(row,col) = (WtoZcollinearCorr_*WtoAcollinearCorr_);
if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*WtoAcollinearCorr_);
if (row==12 || row==14) result(row,col) = (WtoWcollinearCorr_*WtoZcollinearCorr_);
if (row==13 || row==15) result(row,col) = (WtoWcollinearCorr_*WtoAcollinearCorr_);
if (row==16 || row==18) result(row,col) = (WtoWcollinearCorr_*WtoZcollinearCorr_);
if (row==17 || row==19) result(row,col) = (WtoWcollinearCorr_*WtoAcollinearCorr_);
}
if (col==2) {
if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*BtoZcollinearCorr_);
if (row==3 || row==9) result(row,col) = (WtoZcollinearCorr_*BtoAcollinearCorr_);
if (row==4 || row==10) result(row,col) = (WtoAcollinearCorr_*BtoZcollinearCorr_);
if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*BtoAcollinearCorr_);
if (row==14) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_);
if (row==15) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_);
if (row==16) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_);
if (row==17) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_);
}
if (col==3) {
if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*BtoZcollinearCorr_);
if (row==3 || row==9) result(row,col) = (WtoAcollinearCorr_*BtoZcollinearCorr_);
if (row==4 || row==10) result(row,col) = (WtoZcollinearCorr_*BtoAcollinearCorr_);
if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*BtoAcollinearCorr_);
if (row==12) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_);
if (row==13) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_);
if (row==18) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_);
if (row==19) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_);
}
if (col==4) {
if (row==2 || row==8) result(row,col) = (BtoZcollinearCorr_*BtoZcollinearCorr_);
if (row==3 || row==4 || row==9 || row==10) result(row,col) = (BtoZcollinearCorr_*BtoAcollinearCorr_);
if (row==5 || row==11) result(row,col) = (BtoAcollinearCorr_*BtoAcollinearCorr_);
}
// Particle Collinear Corr_ections:
if (process==QQWW) {
if (row<6) result(row,col) *= (ULcollinearCorr_*ULcollinearCorr_);
if ((row>=6)&&(row<12)) result(row,col) *= (DLcollinearCorr_*DLcollinearCorr_);
if (row>=12) result(row,col) *= (ULcollinearCorr_*DLcollinearCorr_);
}
else if (process==LLWW) {
if (row<6) result(row,col) *= (nuLcollinearCorr_*nuLcollinearCorr_);
if ((row>=6)&&(row<12)) result(row,col) *= (ELcollinearCorr_*ELcollinearCorr_);
if (row>=12) result(row,col) *= (nuLcollinearCorr_*ELcollinearCorr_);
}
}
}
}
break;
case QQPhiPhi:
case LLPhiPhi:
{
unsigned int numGauge = 2;
unsigned int numBrokenGauge = 14;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
for (unsigned int row = 0; row < result.size1(); row++) {
// Boson Colinear Corr_ections:
if (row==0 || row==5) result(row,0) = (PhitoWcollinearCorr_*PhitoWcollinearCorr_);
if (row==1 || row==6) result(row,0) = (PhitoZcollinearCorr_*PhitoZcollinearCorr_);
if (row==2 || row==3 || row==7 || row==8) result(row,0) = (PhitoZcollinearCorr_*PhitoHcollinearCorr_);
if (row==4 || row==9) result(row,0) = (PhitoHcollinearCorr_*PhitoHcollinearCorr_);
if (row==10) result(row,0) = (PhitoWcollinearCorr_*PhitoZcollinearCorr_);
if (row==11) result(row,0) = (PhitoWcollinearCorr_*PhitoHcollinearCorr_);
if (row==12) result(row,0) = (PhitoWcollinearCorr_*PhitoZcollinearCorr_);
if (row==13) result(row,0) = (PhitoWcollinearCorr_*PhitoHcollinearCorr_);
// Particle Colinear Corr_ections:
if (process==QQPhiPhi) {
if (row<5) result(row,0) *= (ULcollinearCorr_*ULcollinearCorr_);
if ((row>=5)&&(row<10)) result(row,0) *= (DLcollinearCorr_*DLcollinearCorr_);
if (row>=10) result(row,0) *= (ULcollinearCorr_*DLcollinearCorr_);
}
else if (process==LLPhiPhi) {
if (row<5) result(row,0) *= (nuLcollinearCorr_*nuLcollinearCorr_);
if ((row>=5)&&(row<10)) result(row,0) *= (ELcollinearCorr_*ELcollinearCorr_);
if (row>=10) result(row,0) *= (nuLcollinearCorr_*ELcollinearCorr_);
}
}
DuplicateColumn0(result);
}
break;
case QQWG:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = ULcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoWcollinearCorr_;
result(2,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*WtoZcollinearCorr_;
result(3,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*WtoAcollinearCorr_;
result(4,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoZcollinearCorr_;
result(5,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoAcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QQBG:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_;
result(1,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_;
result(2,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_;
result(3,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QQGG:
{
unsigned int numGauge = 3;
unsigned int numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
result(3,0) = result(4,0) = result(5,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
DuplicateColumn0(result);
}
break;
case QtQtGG:
{
unsigned int numGauge = 3;
unsigned int numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
result(3,0) = result(4,0) = result(5,0) = bLcollinearCorr_*bLcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
DuplicateColumn0(result);
}
break;
case UUBB:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = URcollinearCorr_*URcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_;
result(1,0) = URcollinearCorr_*URcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_;
result(2,0) = URcollinearCorr_*URcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_;
result(3,0) = URcollinearCorr_*URcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_;
DuplicateColumn0(result);
}
break;
case UUPhiPhi:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 5;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = URcollinearCorr_*URcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_;
result(1,0) = URcollinearCorr_*URcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_;
result(2,0) = URcollinearCorr_*URcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_;
result(3,0) = URcollinearCorr_*URcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_;
result(4,0) = URcollinearCorr_*URcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_;
DuplicateColumn0(result);
}
break;
case UUBG:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_;
result(1,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_;
DuplicateColumn0(result);
}
break;
case UUGG:
{
unsigned int numGauge = 3;
unsigned int numBrokenGauge = 3;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = result(2,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
DuplicateColumn0(result);
}
break;
case tRtRGG:
{
unsigned int numGauge = 3;
unsigned int numBrokenGauge = 3;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = result(2,0) = tRcollinearCorr_*tRcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
DuplicateColumn0(result);
}
break;
case DDBB:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = DRcollinearCorr_*DRcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_;
result(1,0) = DRcollinearCorr_*DRcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_;
result(2,0) = DRcollinearCorr_*DRcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_;
result(3,0) = DRcollinearCorr_*DRcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_;
DuplicateColumn0(result);
}
break;
case DDPhiPhi:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 5;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_;
result(1,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_;
result(2,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_;
result(3,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_;
result(4,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_;
DuplicateColumn0(result);
}
break;
case DDBG:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_;
result(1,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_;
DuplicateColumn0(result);
}
break;
case DDGG:
{
unsigned int numGauge = 3;
unsigned int numBrokenGauge = 3;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = result(1,0) = result(2,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*GcollinearCorr_;
DuplicateColumn0(result);
}
break;
case EEBB:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = ERcollinearCorr_*ERcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_;
result(1,0) = ERcollinearCorr_*ERcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_;
result(2,0) = ERcollinearCorr_*ERcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_;
result(3,0) = ERcollinearCorr_*ERcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_;
DuplicateColumn0(result);
}
break;
case EEPhiPhi:
{
unsigned int numGauge = 1;
unsigned int numBrokenGauge = 5;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge); result *= 0.0;
result(0,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_;
result(1,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_;
result(2,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_;
result(3,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_;
result(4,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_;
DuplicateColumn0(result);
}
break;
default:
assert(false);
}
// This is done at the end instead of the beginning for result.size1() and cols()
if (!oneLoop) {
boost::numeric::ublas::matrix<Complex> OnesMatrix(result.size1(),result.size2());
for (unsigned int i=0; i<OnesMatrix.size1(); i++) {
for (unsigned int j=0; j<OnesMatrix.size2(); j++) {
OnesMatrix(i,j) = 1.0;
}
}
return OnesMatrix;
}
// Only include the following for the FO calculation:
for (unsigned int i=0; i<result.size1(); i++) {
for (unsigned int j=0; j<result.size2(); j++) {
- result(i,j) = 1.0 + std::log(result(i,j));
+ result(i,j) = 1.0 + log(result(i,j));
}
}
return result;
}
+
+boost::numeric::ublas::matrix<Complex>
+CollinearSudakov::highEnergyRunning(Energy highScale, Energy EWScale, Energy2 s,
+ Herwig::EWProcess::Process process,
+ bool fixedOrder) {
+ using namespace EWProcess;
+ // perform the calculation
+ evaluateHighScale(highScale,EWScale,s);
+ Complex colW(highColW_);
+ Complex colB(highColB_);
+ Complex colG(highColG_);
+ Complex colQ(highColQ_);
+ Complex colQt(highColQt_);
+ Complex colU(highColU_);
+ Complex coltR(highColtR_);
+ Complex colD(highColD_);
+ Complex colL(highColL_);
+ Complex colE(highColE_);
+ Complex colPhi(highColPhi_);
+ if (fixedOrder) {
+ /* colX not necessarily positive for s = (1000TeV)^2 for the following:
+ colW = log(colW.real());
+ colB = log(colB.real());
+ colPhi = log(colPhi.real());
+ */
+ colG = log(colG.real());
+ colQ = log(colQ.real());
+ colQt = log(colQt.real());
+ colU = log(colU.real());
+ coltR = log(coltR.real());
+ colD = log(colD.real());
+ colL = log(colL.real());
+ colE = log(colE.real());
+ }
+ // set up the matrix
+ boost::numeric::ublas::matrix<Complex> result;
+ unsigned int numGauge(0);
+ switch (process) {
+
+ case QQQQ:
+ case QQQQiden:
+ case QtQtQQ:
+ numGauge = 4;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (process!=QtQtQQ) {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQ+colQ+colQ+colQ;
+ }
+ else {
+ result(i,i) = colQ*colQ*colQ*colQ;
+ }
+ }
+ }
+ else {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQt+colQt+colQ+colQ;
+ }
+ else {
+ result(i,i) = colQt*colQt*colQ*colQ;
+ }
+ }
+ }
+ break;
+
+ case QQUU:
+ case QtQtUU:
+ case QQtRtR:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (process==QQUU) {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQ+colQ+colU+colU;
+ }
+ else {
+ result(i,i) = colQ*colQ*colU*colU;
+ }
+ }
+ }
+ else if (process==QtQtUU) {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQt+colQt+colU+colU;
+ }
+ else {
+ result(i,i) = colQt*colQt*colU*colU;
+ }
+ }
+ }
+ else if (process==QQtRtR) {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQ+colQ+coltR+coltR;
+ }
+ else {
+ result(i,i) = colQ*colQ*coltR*coltR;
+ }
+ }
+ }
+ break;
+
+ case QQDD:
+ case QtQtDD:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (process==QQDD) {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQ+colQ+colD+colD;
+ }
+ else {
+ result(i,i) = colQ*colQ*colD*colD;
+ }
+ }
+ }
+ else {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQt+colQt+colD+colD;
+ }
+ else {
+ result(i,i) = colQt*colQt*colD*colD;
+ }
+ }
+ }
+ break;
+
+ case QQLL:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQ+colQ+colL+colL;
+ }
+ else {
+ result(i,i) = colQ*colQ*colL*colL;
+ }
+ }
+ break;
+
+ case QQEE:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colQ+colQ+colE+colE;
+ }
+ else {
+ result(i,i) = colQ*colQ*colE*colE;
+ }
+ }
+ break;
+
+ case UUUU:
+ case UUUUiden:
+ case tRtRUU:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (process!=tRtRUU) {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colU+colU+colU+colU;
+ }
+ else {
+ result(i,i) = colU*colU*colU*colU;
+ }
+ }
+ }
+ else {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+coltR+coltR+colU+colU;
+ }
+ else {
+ result(i,i) = coltR*coltR*colU*colU;
+ }
+ }
+ }
+ break;
+
+ case UUDD:
+ case tRtRDD:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (process==UUDD) {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colU+colU+colD+colD;
+ }
+ else {
+ result(i,i) = colU*colU*colD*colD;
+ }
+ }
+ }
+ else {
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+coltR+coltR+colD+colD;
+ }
+ else {
+ result(i,i) = coltR*coltR*colD*colD;
+ }
+ }
+ }
+ break;
+
+ case UULL:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colU+colU+colL+colL;
+ }
+ else {
+ result(i,i) = colU*colU*colL*colL;
+ }
+ }
+ break;
+
+ case UUEE:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colU+colU+colE+colE;
+ }
+ else {
+ result(i,i) = colU*colU*colE*colE;
+ }
+ }
+ break;
+
+ case DDDD:
+ case DDDDiden:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colD+colD+colD+colD;
+ }
+ else {
+ result(i,i) = colD*colD*colD*colD;
+ }
+ }
+ break;
+
+ case DDLL:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colD+colD+colL+colL;
+ }
+ else {
+ result(i,i) = colD*colD*colL*colL;
+ }
+ }
+ break;
+
+ case DDEE:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colD+colD+colE+colE;
+ }
+ else {
+ result(i,i) = colD*colD*colE*colE;
+ }
+ }
+ break;
+
+ case LLLL:
+ case LLLLiden:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colL+colL+colL+colL;
+ }
+ else {
+ result(i,i) = colL*colL*colL*colL;
+ }
+ }
+ break;
+
+ case LLEE:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colL+colL+colE+colE;
+ }
+ else {
+ result(i,i) = colL*colL*colE*colE;
+ }
+ }
+ break;
+
+ case EEEE:
+ case EEEEiden:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ for (unsigned int i=0; i<numGauge; i++) {
+ if (fixedOrder) {
+ result(i,i) = 1.0+colE+colE+colE+colE;
+ }
+ else {
+ result(i,i) = colE*colE*colE*colE;
+ }
+ }
+ break;
+
+ case QQWW:
+ numGauge = 5;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = 1.0+colQ+colQ+colW+colW;
+ result(2,2) = result(3,3) = 1.0+colQ+colQ+colW+colB;
+ result(4,4) = 1.0+colQ+colQ+colB+colB;
+ }
+ else {
+ result(0,0) = result(1,1) = colQ*colQ*colW*colW;
+ result(2,2) = result(3,3) = colQ*colQ*colW*colB;
+ result(4,4) = colQ*colQ*colB*colB;
+ }
+ break;
+
+ case QQPhiPhi:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = 1.0+colQ+colQ+colPhi+colPhi;
+ }
+ else {
+ result(0,0) = result(1,1) = colQ*colQ*colPhi*colPhi;
+ }
+ break;
+
+ case QQWG:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colQ+colQ+colW+colG;
+ }
+ else {
+ result(0,0) = colQ*colQ*colW*colG;
+ }
+ break;
+
+ case QQBG:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colQ+colQ+colB+colG;
+ }
+ else {
+ result(0,0) = colQ*colQ*colB*colG;
+ }
+ break;
+
+ case QQGG:
+ case QtQtGG:
+ numGauge = 3;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (process==QQGG) {
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = result(2,2) = 1.0+colQ+colQ+colG+colG;
+ }
+ else {
+ result(0,0) = result(1,1) = result(2,2) = colQ*colQ*colG*colG;
+ }
+ }
+ else {
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = result(2,2) = 1.0+colQt+colQt+colG+colG;
+ }
+ else {
+ result(0,0) = result(1,1) = result(2,2) = colQt*colQt*colG*colG;
+ }
+ }
+ break;
+
+ case LLWW:
+ numGauge = 5;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = 1.0+colL+colL+colW+colW;
+ result(2,2) = result(3,3) = 1.0+colL+colL+colW+colB;
+ result(4,4) = 1.0+colL+colL+colB+colB;
+ }
+ else {
+ result(0,0) = result(1,1) = colL*colL*colW*colW;
+ result(2,2) = result(3,3) = colL*colL*colW*colB;
+ result(4,4) = colL*colL*colB*colB;
+ }
+ break;
+
+ case LLPhiPhi:
+ numGauge = 2;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = 1.0+colL+colL+colPhi+colPhi;
+ }
+ else {
+ result(0,0) = result(1,1) = colL*colL*colPhi*colPhi;
+ }
+ break;
+
+ case UUBB:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colU+colU+colB+colB;
+ }
+ else {
+ result(0,0) = colU*colU*colB*colB;
+ }
+ break;
+
+ case UUPhiPhi:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colU+colU+colPhi+colPhi;
+ }
+ else {
+ result(0,0) = colU*colU*colPhi*colPhi;
+ }
+ break;
+
+ case UUBG:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colU+colU+colB+colG;
+ }
+ else {
+ result(0,0) = colU*colU*colB*colG;
+ }
+ break;
+
+ case UUGG:
+ case tRtRGG:
+ numGauge = 3;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (process==UUGG) {
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = result(2,2) = 1.0+colU+colU+colG+colG;
+ }
+ else {
+ result(0,0) = result(1,1) = result(2,2) = colU*colU*colG*colG;
+ }
+ }
+ else {
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = result(2,2) = 1.0+coltR+coltR+colG+colG;
+ }
+ else {
+ result(0,0) = result(1,1) = result(2,2) = coltR*coltR*colG*colG;
+ }
+ }
+ break;
+
+ case DDBB:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colD+colD+colB+colB;
+ }
+ else {
+ result(0,0) = colD*colD*colB*colB;
+ }
+ break;
+
+ case DDPhiPhi:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colD+colD+colPhi+colPhi;
+ }
+ else {
+ result(0,0) = colD*colD*colPhi*colPhi;
+ }
+ break;
+
+ case DDBG:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colD+colD+colB+colG;
+ }
+ else {
+ result(0,0) = colD*colD*colB*colG;
+ }
+ break;
+
+ case DDGG:
+ numGauge = 3;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = result(1,1) = result(2,2) = 1.0+colD+colD+colG+colG;
+ }
+ else {
+ result(0,0) = result(1,1) = result(2,2) = colD*colD*colG*colG;
+ }
+ break;
+
+ case EEBB:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colE+colE+colB+colB;
+ }
+ else {
+ result(0,0) = colE*colE*colB*colB;
+ }
+ break;
+
+ case EEPhiPhi:
+ numGauge = 1;
+ result = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ if (fixedOrder) {
+ result(0,0) = 1.0+colE+colE+colPhi+colPhi;
+ }
+ else {
+ result(0,0) = colE*colE*colPhi*colPhi;
+ }
+ break;
+
+ default:
+ assert(false);
+ }
+
+ return result;
+}
diff --git a/MatrixElement/EW/CollinearSudakov.h b/MatrixElement/EW/CollinearSudakov.h
--- a/MatrixElement/EW/CollinearSudakov.h
+++ b/MatrixElement/EW/CollinearSudakov.h
@@ -1,527 +1,535 @@
// -*- C++ -*-
#ifndef Herwig_CollinearSudakov_H
#define Herwig_CollinearSudakov_H
//
// This is the declaration of the CollinearSudakov class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/Utilities/GSLIntegrator.h"
#include <boost/numeric/ublas/matrix.hpp>
#include "EWProcess.h"
#include "CollinearSudakov.fh"
namespace Herwig {
using namespace ThePEG;
/**
* Struct for the wavefunction corrections
*/
struct WaveFunctionCorrections {
Complex RW;
Complex RA;
Complex RZ;
Complex RAtoZ;
Complex RZtoA;
Complex RPhi;
Complex EW;
Complex EZ;
Complex RPhi3;
Complex RH;
Complex tLuLDiff;
Complex bLdLDiff;
Complex tRuRDiff;
// The following are constants from parameter integrals:
Complex fFW0;
Complex fF0W;
Complex fFZZ;
Complex aHH;
Complex aZZ;
Complex aW0;
Complex a0W;
Complex bHH;
Complex bZZ;
Complex cHH;
Complex cZZ;
Complex cW0;
Complex atHH;
Complex atZZ;
Complex atW0;
Complex at0W;
Complex ctHH;
Complex ctZZ;
Complex ctW0;
Complex btHH;
Complex btZZ;
Complex fs10;
Complex fs1ZW;
Complex fsWZWZ;
Complex fsZW1;
Complex fs01;
Complex fsHW1;
Complex fsHZ1;
Complex fs1HW;
Complex fs1HZ;
};
/**
* Here is the documentation of the CollinearSudakov class.
*
* @see \ref CollinearSudakovInterfaces "The interfaces"
* defined for CollinearSudakov.
*/
class CollinearSudakov: public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
CollinearSudakov();
/**
* The destructor.
*/
virtual ~CollinearSudakov();
//@}
public:
/**
- * Evalaute the electroweka matching as a matrix
+ * Evalaute the electroweak matching as a matrix
*/
boost::numeric::ublas::matrix<Complex>
electroWeakMatching(Energy EWScale, Energy2 s,
Herwig::EWProcess::Process process,
bool oneLoop);
+ /**
+ * Evalaute the high energy running as a matrix
+ */
+ boost::numeric::ublas::matrix<Complex>
+ highEnergyRunning(Energy highScale, Energy EWScale, Energy2 s,
+ Herwig::EWProcess::Process process,
+ bool fixedOrder);
+
public:
/**
* Make plots for tests
*/
void makePlots();
protected:
/**
* Evaluate the high scale contributions
*/
void evaluateHighScale(Energy highScale, Energy EWScale, Energy2 S);
/**
* Evaluate the low scale contributions
*/
void evaluateLowScale(Energy EWScale, Energy lowScale, Energy2 S);
/**
* Evaluate the matching
*/
void evaluateMatching(Energy EWScale,Energy2 S);
public:
/**
* The operator to be integrated
*/
InvEnergy operator ()(Energy mu) const {
if(high_) return highScaleIntegrand(mu);
else return lowScaleIntegrand(mu);
}
/** Argument type for GaussianIntegrator */
typedef Energy ArgType;
/** Return type for GaussianIntegrator */
typedef InvEnergy ValType;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* The integral of the high scale part of the Sudakov
*/
Complex highScaleIntegral(bool SU3, bool SU2, double Y,
Energy2 s, Energy mu_h, Energy mu_l, bool fermion,
bool longitudinal, double yukFactor);
/**
* the integral of the low scale part of the Sudakov
*/
Complex lowScaleIntegral(bool SU3, double Q, Energy2 s,
Energy mu_h, Energy mu_l, bool fermion,
double boostFactor);
protected:
/**
* High-scale integrand
*/
InvEnergy highScaleIntegrand(Energy mu) const;
/**
* Low-scale integrand
*/
InvEnergy lowScaleIntegrand(Energy mu) const;
/**
* Calculate the wavefunction corrections
*/
WaveFunctionCorrections waveFunctionCorrections(Energy EWScale);
/**
* Collinear matiching for W
*/
Complex CollinearDw(Energy2 s, Energy EWScale);
/**
* Collinear matching for Z
*/
Complex CollinearDz(Energy2 s, Energy EWScale);
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
CollinearSudakov & operator=(const CollinearSudakov &);
private:
/**
* Parameters for the integrand
*/
//@{
/**
* Whether high or low scale
*/
bool high_;
/**
* Whether real of imaginary part
*/
bool real_;
/**
* \f$SU(3)\f$
*/
bool SU3_;
/**
*
*/
bool SU2_;
/**
*
*/
double Y_;
/**
*
*/
Energy2 s_;
/**
*
*/
bool fermion_;
/**
*
*/
bool longitudinal_;
/**
*
*/
double yukFactor_;
/**
*
*/
double boostFactor_;
/**
*
*/
double Q_;
//@}
/**
* Parameters
*/
//@{
/**
* Order for the K terms
*/
int K_ORDER_;
/**
* Order for the B terms
*/
int B_ORDER_;
//@}
/**
* Integrator
*/
GSLIntegrator integrator_;
private:
/**
* Storage of the high scale pieces
*/
//@{
/**
*
*/
Complex highColW_;
/**
*
*/
Complex highColB_;
/**
*
*/
Complex highColG_;
/**
*
*/
Complex highColQ_;
/**
*
*/
Complex highColQt_;
/**
*
*/
Complex highColU_;
/**
*
*/
Complex highColtR_;
/**
*
*/
Complex highColD_;
/**
*
*/
Complex highColL_;
/**
*
*/
Complex highColE_;
/**
*
*/
Complex highColPhi_;
//@}
/**
* Storage of the low scale pieces
*/
//@{
/**
*
*/
complex<double> lowColW_;
/**
*
*/
Complex lowColA_;
/**
*
*/
Complex lowColG_;
/**
*
*/
Complex lowColU_;
/**
*
*/
Complex lowColt_;
/**
*
*/
Complex lowColD_;
/**
*
*/
Complex lowColE_;
//@}
/**
* Storage of the matching parameters
*/
//@{
/**
*
*/
Complex ULcollinearCorr_;
/**
*
*/
Complex DLcollinearCorr_;
/**
*
*/
Complex URcollinearCorr_;
/**
*
*/
Complex DRcollinearCorr_;
/**
*
*/
Complex tLcollinearCorr_;
/**
*
*/
Complex tRcollinearCorr_;
/**
*
*/
Complex bLcollinearCorr_;
/**
*
*/
Complex nuLcollinearCorr_;
/**
*
*/
Complex ELcollinearCorr_;
/**
*
*/
Complex ERcollinearCorr_;
/**
*
*/
Complex WtoWcollinearCorr_;
/**
*
*/
Complex WtoZcollinearCorr_;
/**
*
*/
Complex WtoAcollinearCorr_;
/**
*
*/
Complex BtoZcollinearCorr_;
/**
*
*/
Complex BtoAcollinearCorr_;
/**
*
*/
Complex PhitoWcollinearCorr_;
/**
*
*/
Complex PhitoZcollinearCorr_;
/**
*
*/
Complex PhitoHcollinearCorr_;
/**
*
*/
Complex GcollinearCorr_;
//@}
};
}
#endif /* Herwig_CollinearSudakov_H */
diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
--- a/MatrixElement/EW/GroupInvariants.h
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -1,311 +1,328 @@
// -*- C++ -*-
//
// GroupInvariants.h is a part of Herwig - A multi-purpose Monte Carlo event generator
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
//
#ifndef HERWIG_GroupInvariants_H
#define HERWIG_GroupInvariants_H
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Config/Unitsystem.h"
#include <cassert>
#include <boost/numeric/ublas/matrix.hpp>
namespace Herwig {
using namespace ThePEG;
namespace GroupInvariants {
/**
* Simple struct for storing the different gauge contributions
*/
struct GaugeContributions {
/**
* Default Constructor
*/
GaugeContributions(double inSU3=0.,
double inSU2=0., double inU1=0.)
: SU3(inSU3),SU2(inSU2),U1(inU1)
{}
/**
* \f$SU(3)\f$
*/
double SU3;
/**
* \f$SU(2)\f$
*/
double SU2;
/**
* \f$U(1)\f$
*/
double U1;
};
/**
* The \f$SU(N)\f$ \f$C_A\f$
*/
inline double C_A(unsigned int N) {
return N !=1 ? double(N) : 0.;
}
/**
* The \f$SU(N)\f$ \f$C_F\f$
*/
inline double C_F(unsigned int N) {
return N !=1 ? 0.5*(double(N*N)-1.)/double(N) : 1.;
}
/*
* The \f$SU(N)\f$ \f$C_d\f$
*/
inline double C_d(unsigned int N) {
return (double(N*N)-4.)/double(N);
}
/**
* The \f$SU(N)\f$\f$C_1\f$
*/
inline double C_1(unsigned int N) {
double N2(N*N);
return 0.25*(N2-1.0)/N2;
}
/**
* \f$T_F\f$
*/
inline double T_F(unsigned int N, bool high) {
if(high) {
return N !=1 ? 0.5 : 5.0/3.0;
}
else {
return N !=1 ? 0.5 : 20.0/3.0;
}
}
/**
* \f$t_S\f$
*/
inline double t_S(unsigned int, bool ) {
return 0.5;
}
/**
* / Number of complex scalars in the fundamental rep. of SU(N)/U(1)
*/
inline double n_S(unsigned int N, bool high) {
if(high) {
if(N==2 || N==1) return 1.0;
else if(N==3) return 0.0;
else assert(false);
}
else {
if(N>=1&&N<=3) return 0.;
else assert(false);
}
}
/**
* Number of Dirac Fermions in the fund. rep. of SU(N) (or U(1) for N==1)
*/
inline double n_F(unsigned int N, bool high) {
if(high) {
if(N==1) return 3.0;
else if(N==2 || N==3) return 6.0;
else assert(false);
}
else {
if(N==1) return 1.0;
else if(N==2) return 0.0;
else if(N==3) return 5.0;
else assert(false);
}
}
/**
* Find K_i for gauge group N. high=false for running at mu<mZ
*/
double K_Factor(unsigned int i,unsigned int N, bool high);
/**
* Find B_i for gauge group N, high energy
*/
double B_Factor(int i, int N, bool fermion, bool longitudinal);
/**
* Find B_i for gauge group N, low energy
*/
double B_Factor_Low(int i, int N, bool fermion, double boostFactor);
/**
* Contributions to the Cusps
*/
GaugeContributions cuspContributions(Energy mu, int K_ORDER, bool high);
/**
* Contributions to B, high energy
*/
GaugeContributions BContributions(Energy mu, int B_ORDER,
bool fermion, bool longitudinal);
/**
* Contributions to B, low energy
*/
GaugeContributions BContributionsLow(Energy mu, int B_ORDER,
bool fermion, double boostFactor);
inline Complex PlusLog(double arg) {
static const Complex I(0,1.0);
if (arg>0.0)
return log(arg);
else if (arg<0.0)
return log(-arg)+I*Constants::pi;
else
assert(false);
}
inline Complex MinusLog(double arg) {
static const Complex I(0,1.0);
if (arg>0.0)
return log(arg);
else if (arg<0.0)
return log(-arg)-I*Constants::pi;
else
assert(false);
}
inline Complex getT(Energy2 s, Energy2 t) {
return MinusLog(-t/GeV2) - MinusLog(-s/GeV2);
}
inline Complex getU(Energy2 s, Energy2 u) {
return MinusLog(-u/GeV2) - MinusLog(-s/GeV2);
}
inline boost::numeric::ublas::matrix<Complex> Gamma2(Complex U, Complex T) {
boost::numeric::ublas::matrix<Complex> output(2,2);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) = (-3.0/2.0)*I*pi + (T+U);
output(1,1) = (-3.0/2.0)*I*pi;
output(0,1) = 2.0*(T-U);
output(1,0) = (3.0/8.0)*(T-U);
return output;
}
inline boost::numeric::ublas::matrix<Complex> Gamma2w(Complex U, Complex T) {
boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(5,5);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) += -I*pi*11.0/4.0;
output(0,1) += U-T;
output(1,0) += 2.0*(U-T);
output(1,1) += -I*pi*11.0/4.0 + (T+U);
output(2,2) += -7.0/4.0*I*pi + (U+T);
output(3,3) += -7.0/4.0*I*pi + (U+T);
output(4,4) += -3.0/4.0*I*pi;
return output;
}
inline boost::numeric::ublas::matrix<Complex> Gamma2Singlet() {
boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(2,2);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) = output(1,1) = -3.0/4.0*I*pi;
return output;
}
+ inline Complex Gamma1(double hypercharge) {
+ Complex I(0,1.0);
+ return -I*Constants::pi*(hypercharge*hypercharge);
+ }
+
+ inline Complex Gamma1(double y1, double y2, Complex T, Complex U) {
+ Complex I(0,1.0);
+ return -I*Constants::pi*(y1*y1+y2*y2) + 2.0*y1*y2*(T-U);
+ }
+
+ inline Complex Gamma1(double y1, double y2, double y3, double y4,
+ Complex T, Complex U) {
+ Complex I(0,1.0);
+ return -I*Constants::pi*(y1*y1+y2*y2+y3*y3+y4*y4)/2.0 +
+ (y1*y4+y2*y3)*T - (y1*y3+y2*y4)*U;
+ }
+
/**
* Number of fermion generations (only used in gauge boson HighCMatching)
*/
inline double n_g() { return 3.0; }
/**
* Number of complex scalars in the fundamental rep. of SU(N)
*/
inline double nSWeyl(unsigned int N, bool high) {
if(high) {
if(N==2 || N==1) return 1.0;
else if (N==3) return 0.0;
else assert(false);
}
else {
if( N==1 || N==3 ) return 0.0;
else assert(false);
}
}
/**
* Number of Weyl Fermions in the fundamental rep. of SU(N)
*/
inline double nFWeyl(unsigned int N, bool high) {
if(high) {
if(N==2 || N==3) return 12.0;
else assert(false);
}
else {
if(N==3) return 10.0;
else if(N==1) return 2.0;
else assert(false);
}
}
inline double TFWeyl(unsigned int) {
return 0.5;
}
inline double tSWeyl(unsigned int) {
return 0.5;
}
inline Complex WFunction(Energy mu, Energy2 s) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ln = MinusLog(-s/(mu*mu));
return (-1.0*ln*ln + 3.0*ln+pi*pi/6.0-8.0);
}
/**
* \fX_N\f% function, v is either t or u
*/
inline Complex XNFunction(unsigned int N, Energy mu, Energy2 s, Energy2 v) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ls = MinusLog(-s/(mu*mu));
return (2.0*C_F(N)*WFunction(mu,s) +
C_A(N)*(2.0*ls*ls - 2.0*MinusLog((s+v)/(mu*mu))*ls -
11.0/3.0*ls + pi*pi + 85.0/9.0) +
(2.0/3.0*ls - 10.0/9.0) * TFWeyl(N) * nFWeyl(N,true) +
(1.0/3.0*ls - 8.0/9.0) * TFWeyl(N) * nSWeyl(N,true));
}
/**
* \f$\Pi_1\f$ function
*/
inline Complex PI1_function(Energy mu, Energy2 s) {
assert(abs(s)>ZERO);
return ((41.0/6.0)*MinusLog(-s/(mu*mu))-104.0/9.0);
}
/**
* \f$\tilde{f}\f$ function, v is either t or u
*/
inline Complex fTildeFunction(Energy mu, Energy2 s, Energy2 v) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ls = MinusLog(-s/GeV2), lv = MinusLog(-v/GeV2);
Complex lsv = MinusLog((s+v)/GeV2);
return (-2.0*double(s/(s+v))*(lv-ls) +
double(s*(s+2.0*v)/((s+v)*(s+v))) * ((lv-ls)*(lv-ls) + pi*pi) +
4.0*MinusLog(-s/(mu*mu))*(lv-lsv));
}
}
}
#endif // HERWIG_GroupInvariants_H

File Metadata

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

Event Timeline