Page MenuHomeHEPForge
Files F8309902

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/EW/CollinearSudakov.cc b/MatrixElement/EW/CollinearSudakov.cc
--- a/MatrixElement/EW/CollinearSudakov.cc
+++ b/MatrixElement/EW/CollinearSudakov.cc
@@ -1,657 +1,1192 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the CollinearSudakov class.
//
#include "CollinearSudakov.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "GroupInvariants.h"
#include "ElectroWeakReweighter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
+namespace {
+
+void DuplicateColumn0(boost::numeric::ublas::matrix<Complex> &orig) {
+ for (unsigned int i=0; i<orig.size1(); i++) {
+ for (unsigned int j=1; j<orig.size2(); j++) {
+ orig(i,j) = orig(i,0);
+ }
+ }
+}
+
+}
+
CollinearSudakov::CollinearSudakov() : K_ORDER_(3),
B_ORDER_(2)
{}
CollinearSudakov::~CollinearSudakov() {}
IBPtr CollinearSudakov::clone() const {
return new_ptr(*this);
}
IBPtr CollinearSudakov::fullclone() const {
return new_ptr(*this);
}
void CollinearSudakov::persistentOutput(PersistentOStream & os) const {
os << K_ORDER_ << B_ORDER_;
}
void CollinearSudakov::persistentInput(PersistentIStream & is, int) {
is >> K_ORDER_ >> B_ORDER_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<CollinearSudakov,Interfaced>
describeHerwigCollinearSudakov("Herwig::CollinearSudakov", "HwMEEW.so");
void CollinearSudakov::Init() {
static ClassDocumentation<CollinearSudakov> documentation
("The CollinearSudakov class implements the collinear evolution");
}
Complex CollinearSudakov::highScaleIntegral(bool SU3, bool SU2, double Y,
Energy2 s, Energy mu_h, Energy mu_l,
bool fermion, bool longitudinal,
double yukFactor) {
SU3_ = SU3;
SU2_ = SU2;
Y_ = Y;
s_ = s;
fermion_ = fermion;
longitudinal_ = longitudinal;
yukFactor_ = yukFactor;
// perform the integral
Complex result;
high_ = true;
// real piece
real_ = true;
result.real(integrator_.value(*this,mu_h,mu_l));
// imaginary piece
real_ = false;
result.imag(integrator_.value(*this,mu_h,mu_l));
// return the answer
return exp(result);
}
Complex CollinearSudakov::lowScaleIntegral(bool SU3, double Q, Energy2 s,
Energy mu_h, Energy mu_l, bool fermion,
double boostFactor) {
SU3_ = SU3;
Q_ = Q;
s_ = s;
fermion_ = fermion;
boostFactor_ = boostFactor;
// perform the integral
Complex result;
high_ = false;
// real piece
real_ = true;
result.real(integrator_.value(*this,mu_h,mu_l));
// imaginary piece
real_ = false;
result.imag(integrator_.value(*this,mu_h,mu_l));
// return the answer
return exp(result);
}
InvEnergy CollinearSudakov::highScaleIntegrand(Energy mu) const {
using namespace GroupInvariants;
using Constants::pi;
Complex gamma = 0.0;
// Include K-factor Contributions (Cusps):
GaugeContributions cusp = cuspContributions(mu,K_ORDER_,high_);
// Include B-factors (B1 for U1, B2 for SU2, B3 for SU3):
GaugeContributions nonCusp = BContributions(mu,B_ORDER_,fermion_,longitudinal_);
// common log
Complex plog = PlusLog(s_/sqr(mu));
// SU(3)
if(SU3_) {
if(fermion_)
gamma += C_F(3)*cusp.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));
+ }
+ }
+
+ return result;
+}
diff --git a/MatrixElement/EW/CollinearSudakov.h b/MatrixElement/EW/CollinearSudakov.h
--- a/MatrixElement/EW/CollinearSudakov.h
+++ b/MatrixElement/EW/CollinearSudakov.h
@@ -1,513 +1,527 @@
// -*- C++ -*-
#ifndef Herwig_CollinearSudakov_H
#define Herwig_CollinearSudakov_H
//
// This is the declaration of the CollinearSudakov class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/Utilities/GSLIntegrator.h"
+#include <boost/numeric/ublas/matrix.hpp>
+#include "EWProcess.h"
#include "CollinearSudakov.fh"
namespace Herwig {
using namespace ThePEG;
/**
* Struct for the wavefunction corrections
*/
struct WaveFunctionCorrections {
Complex RW;
Complex RA;
Complex RZ;
Complex RAtoZ;
Complex RZtoA;
Complex RPhi;
Complex EW;
Complex EZ;
Complex RPhi3;
Complex RH;
Complex tLuLDiff;
Complex bLdLDiff;
Complex tRuRDiff;
// The following are constants from parameter integrals:
Complex fFW0;
Complex fF0W;
Complex fFZZ;
Complex aHH;
Complex aZZ;
Complex aW0;
Complex a0W;
Complex bHH;
Complex bZZ;
Complex cHH;
Complex cZZ;
Complex cW0;
Complex atHH;
Complex atZZ;
Complex atW0;
Complex at0W;
Complex ctHH;
Complex ctZZ;
Complex ctW0;
Complex btHH;
Complex btZZ;
Complex fs10;
Complex fs1ZW;
Complex fsWZWZ;
Complex fsZW1;
Complex fs01;
Complex fsHW1;
Complex fsHZ1;
Complex fs1HW;
Complex fs1HZ;
};
/**
* Here is the documentation of the CollinearSudakov class.
*
* @see \ref CollinearSudakovInterfaces "The interfaces"
* defined for CollinearSudakov.
*/
class CollinearSudakov: public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
CollinearSudakov();
/**
* The destructor.
*/
virtual ~CollinearSudakov();
//@}
public:
+
+ /**
+ * Evalaute the electroweka matching as a matrix
+ */
+ boost::numeric::ublas::matrix<Complex>
+ electroWeakMatching(Energy EWScale, Energy2 s,
+ Herwig::EWProcess::Process process,
+ bool oneLoop);
+
+public:
+
+ /**
+ * Make plots for tests
+ */
+ void makePlots();
+
+protected:
/**
* Evaluate the high scale contributions
*/
void evaluateHighScale(Energy highScale, Energy EWScale, Energy2 S);
/**
* Evaluate the low scale contributions
*/
void evaluateLowScale(Energy EWScale, Energy lowScale, Energy2 S);
/**
* Evaluate the matching
*/
void evaluateMatching(Energy EWScale,Energy2 S);
- /**
- * Make plots for tests
- */
- void makePlots();
-
public:
/**
* The operator to be integrated
*/
InvEnergy operator ()(Energy mu) const {
if(high_) return highScaleIntegrand(mu);
else return lowScaleIntegrand(mu);
}
/** Argument type for GaussianIntegrator */
typedef Energy ArgType;
/** Return type for GaussianIntegrator */
typedef InvEnergy ValType;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* The integral of the high scale part of the Sudakov
*/
Complex highScaleIntegral(bool SU3, bool SU2, double Y,
Energy2 s, Energy mu_h, Energy mu_l, bool fermion,
bool longitudinal, double yukFactor);
/**
* the integral of the low scale part of the Sudakov
*/
Complex lowScaleIntegral(bool SU3, double Q, Energy2 s,
Energy mu_h, Energy mu_l, bool fermion,
double boostFactor);
protected:
/**
* High-scale integrand
*/
InvEnergy highScaleIntegrand(Energy mu) const;
/**
* Low-scale integrand
*/
InvEnergy lowScaleIntegrand(Energy mu) const;
/**
* Calculate the wavefunction corrections
*/
WaveFunctionCorrections waveFunctionCorrections(Energy EWScale);
/**
* Collinear matiching for W
*/
Complex CollinearDw(Energy2 s, Energy EWScale);
/**
* Collinear matching for Z
*/
Complex CollinearDz(Energy2 s, Energy EWScale);
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
CollinearSudakov & operator=(const CollinearSudakov &);
private:
/**
* Parameters for the integrand
*/
//@{
/**
* Whether high or low scale
*/
bool high_;
/**
* Whether real of imaginary part
*/
bool real_;
/**
* \f$SU(3)\f$
*/
bool SU3_;
/**
*
*/
bool SU2_;
/**
*
*/
double Y_;
/**
*
*/
Energy2 s_;
/**
*
*/
bool fermion_;
/**
*
*/
bool longitudinal_;
/**
*
*/
double yukFactor_;
/**
*
*/
double boostFactor_;
/**
*
*/
double Q_;
//@}
/**
* Parameters
*/
//@{
/**
* Order for the K terms
*/
int K_ORDER_;
/**
* Order for the B terms
*/
int B_ORDER_;
//@}
/**
* Integrator
*/
GSLIntegrator integrator_;
private:
/**
* Storage of the high scale pieces
*/
//@{
/**
*
*/
Complex highColW_;
/**
*
*/
Complex highColB_;
/**
*
*/
Complex highColG_;
/**
*
*/
Complex highColQ_;
/**
*
*/
Complex highColQt_;
/**
*
*/
Complex highColU_;
/**
*
*/
Complex highColtR_;
/**
*
*/
Complex highColD_;
/**
*
*/
Complex highColL_;
/**
*
*/
Complex highColE_;
/**
*
*/
Complex highColPhi_;
//@}
/**
* Storage of the low scale pieces
*/
//@{
/**
*
*/
complex<double> lowColW_;
/**
*
*/
Complex lowColA_;
/**
*
*/
Complex lowColG_;
/**
*
*/
Complex lowColU_;
/**
*
*/
Complex lowColt_;
/**
*
*/
Complex lowColD_;
/**
*
*/
Complex lowColE_;
//@}
/**
* Storage of the matching parameters
*/
//@{
/**
*
*/
Complex ULcollinearCorr_;
/**
*
*/
Complex DLcollinearCorr_;
/**
*
*/
Complex URcollinearCorr_;
/**
*
*/
Complex DRcollinearCorr_;
/**
*
*/
Complex tLcollinearCorr_;
/**
*
*/
Complex tRcollinearCorr_;
/**
*
*/
Complex bLcollinearCorr_;
/**
*
*/
Complex nuLcollinearCorr_;
/**
*
*/
Complex ELcollinearCorr_;
/**
*
*/
Complex ERcollinearCorr_;
/**
*
*/
Complex WtoWcollinearCorr_;
/**
*
*/
Complex WtoZcollinearCorr_;
/**
*
*/
Complex WtoAcollinearCorr_;
/**
*
*/
Complex BtoZcollinearCorr_;
/**
*
*/
Complex BtoAcollinearCorr_;
/**
*
*/
Complex PhitoWcollinearCorr_;
/**
*
*/
Complex PhitoZcollinearCorr_;
/**
*
*/
Complex PhitoHcollinearCorr_;
/**
*
*/
Complex GcollinearCorr_;
//@}
};
}
#endif /* Herwig_CollinearSudakov_H */
diff --git a/MatrixElement/EW/ElectroWeakMatching.cc b/MatrixElement/EW/ElectroWeakMatching.cc
--- a/MatrixElement/EW/ElectroWeakMatching.cc
+++ b/MatrixElement/EW/ElectroWeakMatching.cc
@@ -1,917 +1,874 @@
// -*- C++ -*-
//
// ElectroWeakMatching.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
//
#include "ElectroWeakMatching.h"
#include "ElectroWeakReweighter.h"
#include "GroupInvariants.h"
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/vector.hpp>
using namespace Herwig;
using namespace ElectroWeakMatching;
using namespace GroupInvariants;
using namespace EWProcess;
-namespace {
-
-Complex getT(Energy2 s, Energy2 t) {
- return MinusLog(-t/GeV2) - MinusLog(-s/GeV2);
-}
-
-Complex getU(Energy2 s, Energy2 u) {
- return MinusLog(-u/GeV2) - MinusLog(-s/GeV2);
-}
-
-boost::numeric::ublas::matrix<Complex> Gamma2(Complex U, Complex T) {
- boost::numeric::ublas::matrix<Complex> output(2,2);
- static const Complex I(0,1.0);
- using Constants::pi;
- output(0,0) = (-3.0/2.0)*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;
-}
-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;
-}
-
-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;
-}
-}
-
boost::numeric::ublas::matrix<Complex>
ElectroWeakMatching::electroWeakMatching(Energy mu,
Energy2 s, Energy2 t, Energy2 u,
Herwig::EWProcess::Process process,
bool oneLoop) {
static const Complex I(0,1.0);
using Constants::pi;
Complex T = getT(s,t);
Complex U = getU(s,u);
// Z-Couplings
double g_Lu = ElectroWeakReweighter::coupling()->g_Lu(mu);
double g_Ld = ElectroWeakReweighter::coupling()->g_Ld(mu);
double g_Le = ElectroWeakReweighter::coupling()->g_Le(mu);
double g_Lnu = ElectroWeakReweighter::coupling()->g_Lnu(mu);
double g_Ru = ElectroWeakReweighter::coupling()->g_Ru(mu);
double g_Rd = ElectroWeakReweighter::coupling()->g_Rd(mu);
double g_Re = ElectroWeakReweighter::coupling()->g_Re(mu);
double g_W = ElectroWeakReweighter::coupling()->g_W(mu);
double g_phiPlus = ElectroWeakReweighter::coupling()->g_phiPlus(mu);
// Weinberg Angle:
double cos2 = ElectroWeakReweighter::coupling()->Cos2thW(mu);
double sin2 = 1.0-cos2;
double cos = sqrt(cos2);
double sin = sqrt(sin2);
boost::numeric::ublas::matrix<Complex> R0,G2,Dw,Dz;
switch (process) {
case QQQQ:
case QQQQiden:
case QtQtQQ:
{
unsigned int numGauge = 4, numBrokenGauge = 12;
boost::numeric::ublas::matrix<Complex> R0=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
boost::numeric::ublas::matrix<Complex> G2=boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
boost::numeric::ublas::matrix<Complex> Dw=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
boost::numeric::ublas::matrix<Complex> Dz=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
R0(0,0) = R0(3,0) = 0.25;
R0(1,0) = R0(2,0) = -0.25;
R0(4,0) = R0(5,0) = 0.5;
R0(6,3) = R0(7,3) = R0(8,3) = R0(9,3) = 1.0;
R0(6,2) = R0(9,2) = 0.25;
R0(7,2) = R0(8,2) = -0.25;
R0(10,2) = R0(11,2) = 0.5;
if (oneLoop) {
double g11 = g_Lu;
double g12 = g_Ld;
double g21 = g_Lu;
double g22 = g_Ld;
for(unsigned int ix=0;ix<numBrokenGauge;++ix) {
Dw(ix,ix) = 0.5*I*pi;
}
Complex w1 = -0.5*(T-U);
Complex w2 = -0.5*(T+U);
for(unsigned int ix=0;ix<numBrokenGauge;ix+=6) {
Dw(ix+0,ix+0) += w1;
Dw(ix+3,ix+3) += w1;
Dw(ix+1,ix+1) +=-w1;
Dw(ix+2,ix+2) +=-w1;
Dw(ix+4,ix+4) += w2;
Dw(ix+5,ix+5) += w2;
}
Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
- 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
for(unsigned int ix=0;ix<numBrokenGauge;ix+=6) {
Dz(ix+0,ix+0) = z1;
Dz(ix+1,ix+1) = z2;
Dz(ix+2,ix+2) = z3;
Dz(ix+3,ix+3) = z4;
Dz(ix+4,ix+4) = z5;
Dz(ix+5,ix+5) = z5;
}
boost::numeric::ublas::matrix<Complex> gam2 = Gamma2(U,T);
G2(0,0) += gam2(0,0);
G2(0,1) += gam2(0,1);
G2(1,0) += gam2(1,0);
G2(1,1) += gam2(1,1);
G2(2,2) += gam2(0,0);
G2(2,3) += gam2(0,1);
G2(3,2) += gam2(1,0);
G2(3,3) += gam2(1,1);
}
}
break;
case QQUU:
case QtQtUU:
case QQtRtR:
{
unsigned int numGauge = 2, numBrokenGauge = 4;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,0) = R0(2,1) = R0(3,1) = 1.0;
if (oneLoop) {
double g11 = g_Lu;
double g12 = g_Ld;
//double g21 = g_Ru;
double g22 = g_Ru;
Complex w1 = 0.25*I*pi;
for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
for(unsigned int ix=0;ix<numBrokenGauge;ix+=2) {
Dz(ix+0,ix+0) = z1;
Dz(ix+1,ix+1) = z2;
}
G2 = Gamma2Singlet();
}
}
break;
case QQDD:
case QtQtDD:
{
unsigned int numGauge = 2, numBrokenGauge = 4;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,0) = R0(2,1) = R0(3,1) = 1.0;
if (oneLoop) {
double g11 = g_Lu;
double g12 = g_Ld;
//double g21 = g_Rd;
double g22 = g_Rd;
Complex w1 = 0.25*I*pi;
for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
for(unsigned int ix=0;ix<numBrokenGauge;ix+=2) {
Dz(ix+0,ix+0) = z1;
Dz(ix+1,ix+1) = z2;
}
G2 = Gamma2Singlet();
}
}
break;
case QQLL:
{
unsigned int numGauge = 2, numBrokenGauge = 6;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
R0(0,0) = R0(3,0) = 0.25;
R0(1,0) = R0(2,0) = -0.25;
R0(4,0) = R0(5,0) = 0.5;
if (oneLoop) {
double g11 = g_Lu;
double g12 = g_Ld;
double g21 = g_Lnu;
double g22 = g_Le;
for (unsigned int i=0; i<6; ++i) {
Dw(i,i) = 0.5*I*pi;
}
Complex w1 = (-1.0/2.0)*(T-U);
Complex w2 = (-1.0/2.0)*(T+U);
Dw(0,0) += w1;
Dw(3,3) += w1;
Dw(1,1) += -1.0*w1;
Dw(2,2) += -1.0*w1;
Dw(4,4) += w2;
Dw(5,5) += w2;
Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
- 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
Dz(0,0) = z1;
Dz(1,1) = z2;
Dz(2,2) = z3;
Dz(3,3) = z4;
Dz(4,4) = Dz(5,5) = z5;
G2 = Gamma2(U,T);
}
}
break;
case QQEE:
{
unsigned int numGauge = 1, numBrokenGauge = 2;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,0) = 1.0;
if (oneLoop) {
double g11 = g_Lu;
double g12 = g_Ld;
//double g21 = g_Re;
double g22 = g_Re;
Complex w1 = 0.25*I*pi;
Dw(0,0) = Dw(1,1) = w1;
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
Dz(0,0) = z1;
Dz(1,1) = z2;
G2(0,0) = Gamma2Singlet()(0,0);
}
}
break;
case UUUU:
case UUUUiden:
case tRtRUU:
{
unsigned int numGauge = 2, numBrokenGauge = 2;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,1) = 1.0;
if (oneLoop) {
double g11 = g_Ru;
//double g12 = g_Ru;
//double g21 = g_Ru;
double g22 = g_Ru;
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Dz(0,0) = Dz(1,1) = z1;
}
}
break;
case UUDD:
case tRtRDD:
{
unsigned int numGauge = 2, numBrokenGauge = 2;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,1) = 1.0;
if (oneLoop) {
double g11 = g_Ru;
//double g12 = g_Ru;
//double g21 = g_Rd;
double g22 = g_Rd;
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Dz(0,0) = Dz(1,1) = z1;
}
}
break;
case UULL:
{
unsigned int numGauge = 1, numBrokenGauge = 2;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,0) = 1.0;
if (oneLoop) {
double g11 = g_Lnu;
double g12 = g_Le;
//double g21 = g_Ru;
double g22 = g_Ru;
Complex w1 = 0.25*I*pi;
Dw(0,0) = Dw(1,1) = w1;
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
Dz(0,0) = z1;
Dz(1,1) = z2;
G2(0,0) = Gamma2Singlet()(0,0);
}
}
break;
case UUEE:
{
unsigned int numGauge = 1, numBrokenGauge = 1;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = 1.0;
if (oneLoop) {
double g11 = g_Ru;
//double g12 = g_Ru;
//double g21 = g_Re;
double g22 = g_Re;
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Dz(0,0) = z1;
}
}
break;
case DDDD:
case DDDDiden:
{
unsigned int numGauge = 2, numBrokenGauge = 2;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,1) = 1.0;
if (oneLoop) {
double g11 = g_Rd;
//double g12 = g_Rd;
//double g21 = g_Rd;
double g22 = g_Rd;
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Dz(0,0) = Dz(1,1) = z1;
}
}
break;
case DDLL:
{
unsigned int numGauge = 1, numBrokenGauge = 2;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,0) = 1.0;
if (oneLoop) {
double g11 = g_Lnu;
double g12 = g_Le;
//double g21 = g_Rd;
double g22 = g_Rd;
Complex w1 = 0.25*I*pi;
Dw(0,0) = Dw(1,1) = w1;
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
Dz(0,0) = z1;
Dz(1,1) = z2;
G2(0,0) = Gamma2Singlet()(0,0);
}
}
break;
case DDEE:
{
unsigned int numGauge = 1, numBrokenGauge = 1;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0 *= 0.0; Dw = Dz *= 0.0;
R0(0,0) = 1.0;
if (oneLoop) {
double g11 = g_Rd;
//double g12 = g_Rd;
//double g21 = g_Re;
double g22 = g_Re;
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Dz(0,0) = z1;
}
}
break;
case LLLL:
case LLLLiden:
{
unsigned int numGauge = 2, numBrokenGauge = 6;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
R0(0,0) = R0(3,0) = 0.25;
R0(1,0) = R0(2,0) = -0.25;
R0(4,0) = R0(5,0) = 0.5;
if (oneLoop) {
double g11 = g_Lnu;
double g12 = g_Le;
double g21 = g_Lnu;
double g22 = g_Le;
for (int i=0; i<6; i++) {
Dw(i,i) = 0.5*I*pi;
}
Complex w1 = (-1.0/2.0)*(T-U);
Complex w2 = (-1.0/2.0)*(T+U);
Dw(0,0) += w1;
Dw(3,3) += w1;
Dw(1,1) += -1.0*w1;
Dw(2,2) += -1.0*w1;
Dw(4,4) += w2;
Dw(5,5) += w2;
Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
- 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
Dz(0,0) = z1;
Dz(1,1) = z2;
Dz(2,2) = z3;
Dz(3,3) = z4;
Dz(4,4) = Dz(5,5) = z5;
G2 = Gamma2(U,T);
}
}
break;
case LLEE:
{
unsigned int numGauge = 1, numBrokenGauge = 2;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(1,0) = 1.0;
if (oneLoop) {
double g11 = g_Lnu;
double g12 = g_Le;
//double g21 = g_Re;
double g22 = g_Re;
Complex w1 = 0.25*I*pi;
Dw(0,0) = Dw(1,1) = w1;
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
Dz(0,0) = z1;
Dz(1,1) = z2;
G2(0,0) = Gamma2Singlet()(0,0);
}
}
break;
case EEEE:
case EEEEiden:
{
unsigned int numGauge = 1, numBrokenGauge = 1;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = 1.0;
if (oneLoop) {
double g11 = g_Re;
//double g12 = g_Re;
//double g21 = g_Re;
double g22 = g_Re;
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
Dz(0,0) = z1;
}
}
break;
case QQWW:
case LLWW:
{
unsigned int numGauge = 5, numBrokenGauge = 20;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = 1.0; R0(0,1) = 0.5;
R0(1,0) = 1.0; R0(1,1) = -0.5;
R0(2,0) = cos2; R0(2,2) = -0.5*sin*cos; R0(2,3) = -0.5*sin*cos; R0(2,4) = sin2;
R0(3,0) = sin*cos; R0(3,2) = 0.5*cos2; R0(3,3) = -0.5*sin2; R0(3,4) = -sin*cos;
R0(4,0) = sin*cos; R0(4,2) = -0.5*sin2; R0(4,3) = 0.5*cos2; R0(4,4) = -sin*cos;
R0(5,0) = sin2; R0(5,2) = 0.5*sin*cos; R0(5,3) = 0.5*sin*cos; R0(5,4) = cos2;
R0(6,0) = 1.0; R0(6,1) = -0.5;
R0(7,0) = 1.0; R0(7,1) = 0.5;
R0(8,0) = cos2; R0(8,2) = 0.5*sin*cos; R0(8,3) = 0.5*sin*cos; R0(8,4) = sin2;
R0(9,0) = sin*cos; R0(9,2) = -0.5*cos2; R0(9,3) = 0.5*sin2; R0(9,4) = -sin*cos;
R0(10,0) = sin*cos; R0(10,2) = 0.5*sin2; R0(10,3) = -0.5*cos2; R0(10,4) = -sin*cos;
R0(11,0) = sin2; R0(11,2) = -0.5*sin*cos; R0(11,3) = -0.5*sin*cos; R0(11,4) = cos2;
R0(12,1) = -cos/sqrt(2.0); R0(12,3) = -sin/sqrt(2.0);
R0(13,1) = -sin/sqrt(2.0); R0(13,3) = cos/sqrt(2.0);
R0(14,1) = cos/sqrt(2.0); R0(14,2) = -sin/sqrt(2.0);
R0(15,1) = sin/sqrt(2.0); R0(15,2) = cos/sqrt(2.0);
R0(16,1) = -cos/sqrt(2.0); R0(16,2) = -sin/sqrt(2.0);
R0(17,1) = -sin/sqrt(2.0); R0(17,2) = cos/sqrt(2.0);
R0(18,1) = cos/sqrt(2.0); R0(18,3) = -sin/sqrt(2.0);
R0(19,1) = sin/sqrt(2.0); R0(19,3) = cos/sqrt(2.0);
if (oneLoop) {
double gW = g_W;
double g1(0.),g2(0.);
if (process==QQWW) {
g1 = g_Lu;
g2 = g_Ld;
}
else if (process==LLWW) {
g1 = g_Lnu;
g2 = g_Le;
}
Complex w1 = T-U+5.0/4.0*I*pi;
Complex w2 = -T+U+5.0/4.0*I*pi;
Complex w3 = -0.5*(T+U) + 3.0/4.0*I*pi;
Complex w4 = 0.25*I*pi;
Dw(0,0) = Dw(7,7) = w1;
Dw(1,1) = Dw(6,6) = w2;
for (unsigned int i=12; i<20; i++) {
Dw(i,i) = w3;
}
Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w4;
Dw(8,8) = Dw(9,9) = Dw(10,10) = Dw(11,11) = w4;
Complex z1 = 2.0*g1*gW*(U-T) - I*pi*(g1*g1+gW*gW);
Complex z2 = 2.0*g1*gW*(T-U) - I*pi*(g1*g1+gW*gW);
Complex z3 = 2.0*g2*gW*(U-T) - I*pi*(g2*g2+gW*gW);
Complex z4 = 2.0*g2*gW*(T-U) - I*pi*(g2*g2+gW*gW);
Complex z5 = -(g2*gW)*T + (g1*gW)*U - I*pi*(g1*g2+g1*gW-g2*gW);
Complex z6 = (g1*gW)*T - (g2*gW)*U - I*pi*(g1*g2+g1*gW-g2*gW);
Complex z7 = -I*pi*g1*g1;
Complex z8 = -I*pi*g2*g2;
Dz(0,0) = z1;
Dz(1,1) = z2;
Dz(2,2) = Dz(3,3) = Dz(4,4) = Dz(5,5) = z7;
Dz(6,6) = z3;
Dz(7,7) = z4;
Dz(8,8) = Dz(9,9) = Dz(10,10) = Dz(11,11) = z8;
Dz(12,12) = Dz(13,13) = Dz(16,16) = Dz(17,17) = z5;
Dz(14,14) = Dz(15,15) = Dz(18,18) = Dz(19,19) = z6;
G2 = Gamma2w(U,T);
}
}
break;
case QQPhiPhi:
case LLPhiPhi:
{
unsigned int numGauge = 2, numBrokenGauge = 14;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = 0.25; R0(0,1) = 1.0;
R0(1,0) = -1.0/8.0; R0(1,1) = 0.5;
R0(2,0) = I/8.0; R0(2,1) = -I/2.0;
R0(3,0) = -I/8.0; R0(3,1) = I/2.0;
R0(4,0) = -1.0/8.0; R0(4,1) = 1.0/2.0;
R0(5,0) = -1.0/4.0; R0(5,1) = 1.0;
R0(6,0) = 1.0/8.0; R0(6,1) = 1.0/2.0;
R0(7,0) = -I/8.0; R0(7,1) = -I/2.0;
R0(8,0) = I/8.0; R0(8,1) = I/2.0;
R0(9,0) = 1.0/8.0; R0(9,1) = 1.0/2.0;
R0(10,0) = -1.0/(2.0*sqrt(2.0));
R0(11,0) = I/(2.0*sqrt(2.0));
R0(12,0) = -1.0/(2.0*sqrt(2.0));
R0(13,0) = -I/(2.0*sqrt(2.0));
if (oneLoop) {
double g1(0.),g2(0.);
if (process==QQWW) {
g1 = g_Lu;
g2 = g_Ld;
}
else if (process==LLWW) {
g1 = g_Lnu;
g2 = g_Le;
}
else
assert(false);
double g3 = g_phiPlus;
Complex w0 = 0.25*I*pi;
Complex w1 = 0.5*(T-U) + 0.5*I*pi;
Complex w2 = -0.5*(T-U) + 0.5*I*pi;
Complex w3 = 0.25*I*(T-U);
Complex w4 = -0.25*(T+U) + 0.25*I*pi;
Dw(0,0) = w2;
Dw(1,1) = w0; Dw(1,2) = w3; Dw(1,3) = -w3; Dw(1,4) = w0;
Dw(2,1) = -w3; Dw(2,2) = w0; Dw(2,3) = -w0; Dw(2,4) = -w3;
Dw(3,1) = w3; Dw(3,2) = -w0; Dw(3,3) = w0; Dw(3,4) = w3;
Dw(4,1) = w0; Dw(4,2) = w3; Dw(4,3) = -w3; Dw(4,4) = w0;
Dw(5,5) = w1;
Dw(6,6) = w0; Dw(6,7) = -w3; Dw(6,8) = w3; Dw(6,9) = w0;
Dw(7,6) = w3; Dw(7,7) = w0; Dw(7,8) = -w0; Dw(7,9) = w3;
Dw(8,6) = -w3; Dw(8,7) = -w0; Dw(8,8) = w0; Dw(8,9) = -w3;
Dw(9,6) = w0; Dw(9,7) = -w3; Dw(9,8) = w3; Dw(9,9) = w0;
Dw(10,10) = w4; Dw(10,11) = I*w4;
Dw(11,10) = -I*w4; Dw(11,11) = w4;
Dw(12,12) = w4; Dw(12,13) = -I*w4;
Dw(13,12) = I*w4; Dw(13,13) = w4;
Complex z1 = 2.0*g3*g1*(T-U) - I*pi*(g3*g3+g1*g1);
Complex z2 = 2.0*g3*g2*(T-U) - I*pi*(g3*g3+g2*g2);
Complex z3 = -I*pi*g1*g1;
Complex z4 = 0.5*I*g1*(T-U);
Complex z5 = 0.25*I*pi;
Complex z6 = -I*pi*g2*g2;
Complex z7 = 0.5*I*g2*(T-U);
Complex z8 = g3*g1*T-g3*g2*U-I*pi*(g1*g2-g2*g3+g1*g3);
Complex z9 = 0.5*I*g2*T-0.5*I*g1*U+pi/2.0*g2-pi/2.0*g1+pi/2.0*g3;
Dz(0,0) = z1;
Dz(1,1) = z3; Dz(1,2) = -z4; Dz(1,3) = z4; Dz(1,4) = -z5;
Dz(2,1) = z4; Dz(2,2) = z3; Dz(2,3) = z5; Dz(2,4) = z4;
Dz(3,1) = -z4; Dz(3,2) = z5; Dz(3,3) = z3; Dz(3,4) = -z4;
Dz(4,1) = -z5; Dz(4,2) = -z4; Dz(4,3) = z4; Dz(4,4) = z3;
Dz(5,5) = z2;
Dz(6,6) = z6; Dz(6,7) = -z7; Dz(6,8) = z7; Dz(6,9) = -z5;
Dz(7,6) = z7; Dz(7,7) = z6; Dz(7,8) = z5; Dz(7,9) = z7;
Dz(8,6) = -z7; Dz(8,7) = z5; Dz(8,8) = z6; Dz(8,9) = -z7;
Dz(9,6) = -z5; Dz(9,7) = -z7; Dz(9,8) = z7; Dz(9,9) = z6;
Dz(10,10) = z8; Dz(10,11) = -z9;
Dz(11,10) = z9; Dz(11,11) = z8;
Dz(12,12) = z8; Dz(12,13) = z9;
Dz(13,12) = -z9; Dz(13,13) = z8;
G2 = Gamma2(U,T);
}
}
break;
case QQWG:
{
unsigned int numGauge = 1, numBrokenGauge = 6;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = 1.0/sqrt(2);
R0(1,0) = 1.0/sqrt(2);
R0(2,0) = cos/2.0;
R0(3,0) = sin/2.0;
R0(4,0) = -cos/2.0;
R0(5,0) = -sin/2.0;
if (oneLoop) {
double g1 = g_Lu;
double g2 = g_Ld;
double gW = g_W;
Complex w1 = -0.5*(T+U) + 0.75*I*pi;
Complex w2 = 0.25*I*pi;
Dw(0,0) = Dw(1,1) = w1;
Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;
Complex z1 = gW*g1*T - gW*g2*U - I*pi*(g1*g2+g1*gW-g2*gW);
Complex z2 = gW*g1*U - gW*g2*T - I*pi*(g2*g1+g1*gW-g2*gW);
Complex z3 = -I*pi*g1*g1;
Complex z4 = -I*pi*g2*g2;
Dz(0,0) = z1;
Dz(1,1) = z2;
Dz(2,2) = z3;
Dz(3,3) = z3;
Dz(4,4) = z4;
Dz(5,5) = z4;
G2(0,0) = -7.0/4.0*I*pi + (U+T);
}
}
break;
case QQBG:
{
unsigned int numGauge = 1, numBrokenGauge = 4;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = -sin;
R0(1,0) = cos;
R0(2,0) = -sin;
R0(3,0) = cos;
if (oneLoop) {
double g1 = g_Lu;
double g2 = g_Ld;
Complex w2 = 0.25*I*pi;
Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = w2;
Complex z3 = -I*pi*g1*g1;
Complex z4 = -I*pi*g2*g2;
Dz(0,0) = z3;
Dz(1,1) = z3;
Dz(2,2) = z4;
Dz(3,3) = z4;
G2(0,0) = Gamma2Singlet()(0,0);
}
}
break;
case QQGG:
case QtQtGG:
{
unsigned int numGauge = 3, numBrokenGauge = 6;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = R0(3,0) = 1.0;
R0(1,1) = R0(4,1) = 1.0;
R0(2,2) = R0(5,2) = 1.0;
if (oneLoop) {
double g1 = g_Lu;
double g2 = g_Ld;
Complex w2 = 0.25*I*pi;
Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;
Complex z3 = -I*pi*g1*g1;
Complex z4 = -I*pi*g2*g2;
Dz(0,0) = Dz(1,1) = Dz(2,2) = z3;
Dz(3,3) = Dz(4,4) = Dz(5,5) = z4;
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 *= 0.0;
G2(0,0) = G2(1,1) = G2(2,2) = Gamma2Singlet()(0,0);
}
}
break;
case UUBB:
case DDBB:
case EEBB:
{
unsigned int numGauge = 1, numBrokenGauge = 4;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = sin2;
R0(1,0) = -sin*cos;
R0(2,0) = -sin*cos;
R0(3,0) = cos2;
if (oneLoop) {
double g1(0.);
if (process==UUBB) {
g1 = g_Ru;
}
else if (process==DDBB) {
g1 = g_Rd;
}
else if (process==EEBB) {
g1 = g_Re;
}
else
assert(false);
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = -I*pi*g1*g1;
Dz(0,0) = Dz(1,1) = Dz(2,2) = Dz(3,3) = z1;
}
}
break;
case UUPhiPhi:
case DDPhiPhi:
case EEPhiPhi:
{
unsigned int numGauge = 1, numBrokenGauge = 5;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = 1.0;
R0(1,0) = 0.5;
R0(2,0) = -0.5*I;
R0(3,0) = 0.5*I;
R0(4,0) = 0.5;
if (oneLoop) {
double g1(0.);
if (process==UUPhiPhi) {
g1 = g_Ru;
}
else if (process==DDPhiPhi) {
g1 = g_Rd;
}
else if (process==EEPhiPhi) {
g1 = g_Re;
}
double g3 = g_phiPlus;
Dw(0,0) = Dw(1,4) = Dw(4,1) = 0.25*I*pi;
Dw(2,3) = Dw(3,2) = -0.25*I*pi;
Complex z1 = 2.0*g3*g1*(T-U) - I*pi*(g3*g3+g1*g1);
Complex z2 = 0.5*I*g1*g1;
Complex z3 = -I*pi*g1*g1;
Complex z4 = 0.25*I*pi;
Dz(0,0) = z1;
Dz(1,1) = z3; Dz(1,2) = -z2; Dz(1,3) = z2; Dz(1,4) = -z4;
Dz(2,1) = z2; Dz(2,2) = z3; Dz(2,3) = z4; Dz(2,4) = z2;
Dz(3,1) = -z2; Dz(3,2) = z4; Dz(3,3) = z3; Dz(3,4) = -z2;
Dz(4,1) = -z4; Dz(4,2) = -z2; Dz(4,3) = z2; Dz(4,4) = z3;
G2(0,0) = Gamma2Singlet()(0,0);
}
}
break;
case UUBG:
case DDBG:
{
unsigned int numGauge = 1, numBrokenGauge = 2;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
R0(0,0) = -sin;
R0(1,0) = cos;
if (oneLoop) {
double g1(0.);
if (process==UUBG) {
g1 = g_Ru;
}
else if (process==DDBG) {
g1 = g_Rd;
}
else
assert(false);
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = -I*pi*g1*g1;
Dz(0,0) = Dz(1,1) = z1;
}
}
break;
case UUGG:
case tRtRGG:
case DDGG:
{
unsigned int numGauge = 3, numBrokenGauge = 3;
R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
if (oneLoop) {
double g1(0.);
if ((process==UUGG)||(process==tRtRGG)) {
g1 = g_Ru;
}
else if (process==DDGG) {
g1 = g_Rd;
}
else
assert(false);
// There is no Dw contribution for two SU(2) singlets.
Complex z1 = -I*pi*g1*g1;
Dz(0,0) = Dz(1,1) = Dz(2,2) = z1;
}
}
break;
default:
assert(false);
}
double aW = ElectroWeakReweighter::coupling()->aW(mu);
double aZ = ElectroWeakReweighter::coupling()->aZ(mu);
Energy mZ = ElectroWeakReweighter::coupling()->mZ();
Energy mW = ElectroWeakReweighter::coupling()->mW();
if (!oneLoop) {
return R0;
}
boost::numeric::ublas::matrix<Complex> temp(R0.size1(),R0.size2());
boost::numeric::ublas::axpy_prod(R0,G2,temp);
R0+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
boost::numeric::ublas::axpy_prod(Dw,R0,temp);
R0+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
boost::numeric::ublas::axpy_prod(Dz,R0,temp);
R0+=aZ/(4.0*pi)*4.0*log(mZ/mu)*temp;
return R0;
}
diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
--- a/MatrixElement/EW/GroupInvariants.h
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -1,269 +1,311 @@
// -*- C++ -*-
//
// GroupInvariants.h is a part of Herwig - A multi-purpose Monte Carlo event generator
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
//
#ifndef HERWIG_GroupInvariants_H
#define HERWIG_GroupInvariants_H
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Config/Unitsystem.h"
#include <cassert>
+#include <boost/numeric/ublas/matrix.hpp>
namespace Herwig {
using namespace ThePEG;
namespace GroupInvariants {
/**
* Simple struct for storing the different gauge contributions
*/
struct GaugeContributions {
/**
* Default Constructor
*/
GaugeContributions(double inSU3=0.,
double inSU2=0., double inU1=0.)
: SU3(inSU3),SU2(inSU2),U1(inU1)
{}
/**
* \f$SU(3)\f$
*/
double SU3;
/**
* \f$SU(2)\f$
*/
double SU2;
/**
* \f$U(1)\f$
*/
double U1;
};
/**
* The \f$SU(N)\f$ \f$C_A\f$
*/
inline double C_A(unsigned int N) {
return N !=1 ? double(N) : 0.;
}
/**
* The \f$SU(N)\f$ \f$C_F\f$
*/
inline double C_F(unsigned int N) {
return N !=1 ? 0.5*(double(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;
+ }
+
/**
* 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:51 PM (14 h, 46 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4017665
Default Alt Text
(98 KB)

Event Timeline