Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8310048
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
162 KB
Subscribers
None
View Options
diff --git a/MatrixElement/EW/CollinearSudakov.cc b/MatrixElement/EW/CollinearSudakov.cc
--- a/MatrixElement/EW/CollinearSudakov.cc
+++ b/MatrixElement/EW/CollinearSudakov.cc
@@ -1,2084 +1,2075 @@
// -*- 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 + 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;
}
boost::numeric::ublas::matrix<Complex>
CollinearSudakov::lowEnergyRunning(Energy EWScale, Energy lowScale, Energy2 s,
Herwig::EWProcess::Process process) {
using namespace EWProcess;
// evaluate the running
evaluateLowScale(EWScale,lowScale,s);
Complex colUL = lowColU_;
Complex colDL = lowColD_;
Complex colUR = lowColU_;
Complex colDR = lowColD_;
Complex coltL = lowColt_;
Complex coltR = lowColt_;
Complex colbL = lowColD_;
Complex colnuL = 1.0;
Complex colEL = lowColE_;
Complex colER = lowColE_;
Complex colW = lowColW_;
Complex colZ = 1.0;
Complex colA = lowColA_;
Complex colPhi = lowColW_;
Complex colPhi3 = 1.0;
Complex colH = 1.0;
Complex colG = lowColG_;
// calculate the matrix
boost::numeric::ublas::matrix<Complex> result;
unsigned int numBrokenGauge;
switch (process) {
case QQQQ:
case QQQQiden:
numBrokenGauge = 12;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(6,6) = colUL*colUL*colUL*colUL;
result(3,3) = result(9,9) = colDL*colDL*colDL*colDL;
for (unsigned int i=0; i<12; i++) {
if (i!=0 && i!=3 && i!=6 && i!=9) {
result(i,i) = colUL*colUL*colDL*colDL;
}
}
break;
case QtQtQQ:
numBrokenGauge = 12;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(6,6) = colUL*colUL*coltL*coltL;
result(3,3) = result(9,9) = colDL*colDL*colbL*colbL;
for (unsigned int i=0; i<12; i++) {
if (i==4 || i==5 || i==10 || i==11) {
result(i,i) = colUL*coltL*colDL*colbL;
}
else if (i==1 || i==7) {
result(i,i) = colDL*colDL*coltL*coltL;
}
else if (i==2 || i==8) {
result(i,i) = colUL*colUL*colbL*colbL;
}
}
break;
case QQUU:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(2,2) = colUL*colUL*colUR*colUR;
result(1,1) = result(3,3) = colDL*colDL*colUR*colUR;
break;
case QtQtUU:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(2,2) = coltL*coltL*colUR*colUR;
result(1,1) = result(3,3) = colbL*colbL*colUR*colUR;
break;
case QQtRtR:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(2,2) = colUL*colUL*coltR*coltR;
result(1,1) = result(3,3) = colDL*colDL*coltR*coltR;
break;
case QQDD:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(2,2) = colUL*colUL*colDR*colDR;
result(1,1) = result(3,3) = colDL*colDL*colDR*colDR;
break;
case QtQtDD:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(2,2) = coltL*coltL*colDR*colDR;
result(1,1) = result(3,3) = colbL*colbL*colDR*colDR;
break;
case QQLL:
numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colnuL*colnuL*colUL*colUL;
result(1,1) = colnuL*colnuL*colDL*colDL;
result(2,2) = colEL*colEL*colUL*colUL;
result(3,3) = colEL*colEL*colDL*colDL;
result(4,4) = result(5,5) = colnuL*colEL*colUL*colDL;
break;
case QQEE:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colUL*colUL*colER*colER;
result(1,1) = colDL*colDL*colER*colER;
break;
case UUUU:
case UUUUiden:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = colUR*colUR*colUR*colUR;
break;
case tRtRUU:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = coltR*coltR*colUR*colUR;
break;
case UUDD:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = colUR*colUR*colDR*colDR;
break;
case tRtRDD:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = coltR*coltR*colDR*colDR;
break;
case UULL:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colnuL*colnuL*colUR*colUR;
result(1,1) = colEL*colEL*colUR*colUR;
break;
case UUEE:
numBrokenGauge = 1;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colUR*colUR*colER*colER;
break;
case DDDD:
case DDDDiden:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = colDR*colDR*colDR*colDR;
break;
case DDLL:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colnuL*colnuL*colDR*colDR;
result(1,1) = colEL*colEL*colDR*colDR;
break;
case DDEE:
numBrokenGauge = 1;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colDR*colDR*colER*colER;
break;
case LLLL:
case LLLLiden:
numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colnuL*colnuL*colnuL*colnuL;
result(1,1) = colnuL*colnuL*colEL*colEL;
result(2,2) = colEL*colEL*colnuL*colnuL;
result(3,3) = colEL*colEL*colEL*colEL;
result(4,4) = result(5,5) = colnuL*colEL*colnuL*colEL;
break;
case LLEE:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colnuL*colnuL*colER*colER;
result(1,1) = colEL*colEL*colER*colER;
break;
case EEEE:
case EEEEiden:
numBrokenGauge = 1;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colER*colER*colER*colER;
break;
case QQWW:
case LLWW:
numBrokenGauge = 20;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
for (unsigned int row = 0; row < result.size1(); row++) {
// Boson Collinear Corrections:
if (row==0 || row==1 || row==6 || row==7) result(row,row) = (colW*colW);
else if (row==2 || row==8) result(row,row) = (colZ*colZ);
else if (row==3 || row==4 || row==9 || row==10) result(row,row) = (colZ*colA);
else if (row==5 || row==11) result(row,row) = (colA*colA);
else if (row==12 || row==14) result(row,row) = (colW*colZ);
else if (row==13 || row==15) result(row,row) = (colW*colA);
else if (row==16 || row==18) result(row,row) = (colW*colZ);
else if (row==17 || row==19) result(row,row) = (colW*colA);
// Particle Collinear Corrections:
if (process==QQWW) {
if (row<6) result(row,row) *= (colUL*colUL);
if ((row>=6)&&(row<12)) result(row,row) *= (colDL*colDL);
if (row>=12) result(row,row) *= (colUL*colDL);
}
else if (process==LLWW) {
if (row<6) result(row,row) *= (colnuL*colnuL);
if ((row>=6)&&(row<12)) result(row,row) *= (colEL*colEL);
if (row>=12) result(row,row) *= (colnuL*colEL);
}
}
break;
case QQPhiPhi:
case LLPhiPhi:
numBrokenGauge = 14;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
for (unsigned int row = 0; row < result.size1(); row++) {
// Boson Colinear Corrections:
if (row==0 || row==5) result(row,row) = (colPhi*colPhi);
else if (row==1 || row==6) result(row,row) = (colPhi3*colPhi3);
else if (row==2 || row==3 || row==7 || row==8) result(row,row) = (colPhi3*colH);
else if (row==4 || row==9) result(row,row) = (colH*colH);
else if (row==10) result(row,row) = (colPhi*colPhi3);
else if (row==11) result(row,row) = (colPhi*colH);
else if (row==12) result(row,row) = (colPhi*colPhi3);
else if (row==13) result(row,row) = (colPhi*colH);
// Particle Colinear Corrections:
if (process==QQPhiPhi) {
if (row<5) result(row,row) *= (colUL*colUL);
if ((row>=5)&&(row<10)) result(row,row) *= (colDL*colDL);
if (row>=10) result(row,row) *= (colUL*colDL);
}
else if (process==LLPhiPhi) {
if (row<5) result(row,row) *= (colnuL*colnuL);
if ((row>=5)&&(row<10)) result(row,row) *= (colEL*colEL);
if (row>=10) result(row,row) *= (colnuL*colEL);
}
}
break;
case QQWG:
numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = colUL*colDL*colG*colW;
result(2,2) = colUL*colUL*colG*colZ;
result(3,3) = colUL*colUL*colG*colA;
result(4,4) = colDL*colDL*colG*colZ;
result(5,5) = colDL*colDL*colG*colA;
break;
case QQBG:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colUL*colUL*colG*colZ;
result(1,1) = colUL*colUL*colG*colA;
result(2,2) = colDL*colDL*colG*colZ;
result(3,3) = colDL*colDL*colG*colA;
break;
case QQGG:
numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = result(2,2) = colUL*colUL*colG*colG;
result(3,3) = result(4,4) = result(5,5) = colDL*colDL*colG*colG;
break;
case QtQtGG:
numBrokenGauge = 6;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = result(2,2) = coltL*coltL*colG*colG;
result(3,3) = result(4,4) = result(5,5) = colbL*colbL*colG*colG;
break;
case UUBB:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colUR*colUR*colZ*colZ;
result(1,1) = colUR*colUR*colZ*colA;
result(2,2) = colUR*colUR*colA*colZ;
result(3,3) = colUR*colUR*colA*colA;
break;
case UUPhiPhi:
numBrokenGauge = 5;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colUR*colUR*colPhi*colPhi;
result(1,1) = colUR*colUR*colPhi3*colPhi3;
result(2,2) = colUR*colUR*colH*colPhi3;
result(3,3) = colUR*colUR*colPhi3*colH;
result(4,4) = colUR*colUR*colH*colH;
break;
case UUBG:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colUR*colUR*colG*colZ;
result(1,1) = colUR*colUR*colG*colA;
break;
case UUGG:
numBrokenGauge = 3;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = result(2,2) = colUR*colUR*colG*colG;
break;
case tRtRGG:
numBrokenGauge = 3;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = result(2,2) = coltR*coltR*colG*colG;
break;
case DDBB:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colDR*colDR*colZ*colZ;
result(1,1) = colDR*colDR*colZ*colA;
result(2,2) = colDR*colDR*colA*colZ;
result(3,3) = colDR*colDR*colA*colA;
break;
case DDPhiPhi:
numBrokenGauge = 5;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colDR*colDR*colPhi*colPhi;
result(1,1) = colDR*colDR*colPhi3*colPhi3;
result(2,2) = colDR*colDR*colH*colPhi3;
result(3,3) = colDR*colDR*colPhi3*colH;
result(4,4) = colDR*colDR*colH*colH;
break;
case DDBG:
numBrokenGauge = 2;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colDR*colDR*colG*colZ;
result(1,1) = colDR*colDR*colG*colA;
break;
case DDGG:
numBrokenGauge = 3;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = result(1,1) = result(2,2) = colDR*colDR*colG*colG;
break;
case EEBB:
numBrokenGauge = 4;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colER*colER*colZ*colZ;
result(1,1) = colER*colER*colZ*colA;
result(2,2) = colER*colER*colA*colZ;
result(3,3) = colER*colER*colA*colA;
break;
case EEPhiPhi:
numBrokenGauge = 5;
result = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
result(0,0) = colER*colER*colPhi*colPhi;
result(1,1) = colER*colER*colPhi3*colPhi3;
result(2,2) = colER*colER*colH*colPhi3;
result(3,3) = colER*colER*colPhi3*colH;
result(4,4) = colER*colER*colH*colH;
break;
default:
assert(false);
}
return result;
}
diff --git a/MatrixElement/EW/ElectroWeakMatching.cc b/MatrixElement/EW/ElectroWeakMatching.cc
--- a/MatrixElement/EW/ElectroWeakMatching.cc
+++ b/MatrixElement/EW/ElectroWeakMatching.cc
@@ -1,874 +1,876 @@
// -*- 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;
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=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;
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) {
+ if (process==QQPhiPhi) {
g1 = g_Lu;
g2 = g_Ld;
}
- else if (process==LLWW) {
+ else if (process==LLPhiPhi) {
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);
+ R0(0,0) = R0(1,1) = R0(2,2) = 1.0;
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> output(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;
+ output+=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;
+ output+=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;
+ output+=aZ/(4.0*pi)*4.0*log(mZ/mu)*temp;
+ return output;
}
diff --git a/MatrixElement/EW/ElectroWeakReweighter.cc b/MatrixElement/EW/ElectroWeakReweighter.cc
--- a/MatrixElement/EW/ElectroWeakReweighter.cc
+++ b/MatrixElement/EW/ElectroWeakReweighter.cc
@@ -1,103 +1,200 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ElectroWeakReweighter class.
//
#include "ElectroWeakReweighter.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
+#include "boost/numeric/ublas/matrix.hpp"
+#include "boost/numeric/ublas/operation.hpp"
+#include "EWProcess.h"
+#include "HighEnergyMatching.h"
+#include "ElectroWeakMatching.h"
using namespace Herwig;
tEWCouplingsPtr ElectroWeakReweighter::staticEWCouplings_ = tEWCouplingsPtr();
ElectroWeakReweighter::ElectroWeakReweighter() {}
ElectroWeakReweighter::~ElectroWeakReweighter() {}
IBPtr ElectroWeakReweighter::clone() const {
return new_ptr(*this);
}
IBPtr ElectroWeakReweighter::fullclone() const {
return new_ptr(*this);
}
void ElectroWeakReweighter::persistentOutput(PersistentOStream & os) const {
- os << EWCouplings_ << collinearSudakov_;
+ os << EWCouplings_ << collinearSudakov_ << softSudakov_;
}
void ElectroWeakReweighter::persistentInput(PersistentIStream & is, int) {
- is >> EWCouplings_ >> collinearSudakov_;
+ is >> EWCouplings_ >> collinearSudakov_ >> softSudakov_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<ElectroWeakReweighter,ReweightBase>
describeHerwigElectroWeakReweighter("Herwig::ElectroWeakReweighter", "HwMEEW.so");
void ElectroWeakReweighter::Init() {
static ClassDocumentation<ElectroWeakReweighter> documentation
("There is no documentation for the ElectroWeakReweighter class");
static Reference<ElectroWeakReweighter,EWCouplings> interfaceEWCouplings
("EWCouplings",
"The object to calculate the electroweak couplings",
&ElectroWeakReweighter::EWCouplings_, false, false, true, false, false);
static Reference<ElectroWeakReweighter,CollinearSudakov> interfaceCollinearSudakov
("CollinearSudakov",
"The collinear Sudakov",
&ElectroWeakReweighter::collinearSudakov_, false, false, true, false, false);
+ static Reference<ElectroWeakReweighter,SoftSudakov> interfaceSoftSudakov
+ ("SoftSudakov",
+ "The soft Sudakov",
+ &ElectroWeakReweighter::softSudakov_, false, false, true, false, false);
+
}
+namespace {
+
+void axpy_prod_local(const boost::numeric::ublas::matrix<Complex> & A,
+ const boost::numeric::ublas::matrix<complex<InvEnergy2> > & B,
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
+ assert(A.size2()==B.size1());
+ C.resize(A.size1(),B.size2());
+ for(unsigned int ix=0;ix<A.size1();++ix) {
+ for(unsigned int iy=0;iy<B.size2();++iy) {
+ C(ix,iy) = ZERO;
+ for(unsigned int iz=0;iz<A.size2();++iz) {
+ C(ix,iy) += A(ix,iz)*B(iz,iy);
+ }
+ }
+ }
+}
+
+void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> > & A,
+ const boost::numeric::ublas::matrix<Complex> & B,
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
+ assert(A.size2()==B.size1());
+ C.resize(A.size1(),B.size2());
+ for(unsigned int ix=0;ix<A.size1();++ix) {
+ for(unsigned int iy=0;iy<B.size2();++iy) {
+ C(ix,iy) = ZERO;
+ for(unsigned int iz=0;iz<A.size2();++iz) {
+ C(ix,iy) += A(ix,iz)*B(iz,iy);
+ }
+ }
+ }
+}
+
+}
double ElectroWeakReweighter::weight() const {
EWCouplings_->initialize();
staticEWCouplings_ = EWCouplings_;
- cerr << "aEM\n";
- for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.1) {
- cerr << scale/GeV << " "
- << EWCouplings_->aEM(scale) << "\n";
- }
- cerr << "aS\n";
- for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
- cerr << scale/GeV << " "
- << EWCouplings_->aS(scale) << "\n";
- }
- cerr << "y_t\n";
- for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
- cerr << scale/GeV << " "
- << EWCouplings_->y_t(scale) << "\n";
- }
- cerr << "lambda\n";
- for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
- cerr << scale/GeV << " "
- << EWCouplings_->lambda(scale) << "\n";
- }
- cerr << "vev\n";
- for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
- cerr << scale/GeV << " "
- << EWCouplings_->vev(scale)/GeV << "\n";
- }
+ // cerr << "aEM\n";
+ // for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.1) {
+ // cerr << scale/GeV << " "
+ // << EWCouplings_->aEM(scale) << "\n";
+ // }
+ // cerr << "aS\n";
+ // for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
+ // cerr << scale/GeV << " "
+ // << EWCouplings_->aS(scale) << "\n";
+ // }
+ // cerr << "y_t\n";
+ // for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
+ // cerr << scale/GeV << " "
+ // << EWCouplings_->y_t(scale) << "\n";
+ // }
+ // cerr << "lambda\n";
+ // for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
+ // cerr << scale/GeV << " "
+ // << EWCouplings_->lambda(scale) << "\n";
+ // }
+ // cerr << "vev\n";
+ // for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
+ // cerr << scale/GeV << " "
+ // << EWCouplings_->vev(scale)/GeV << "\n";
+ // }
collinearSudakov_->makePlots();
+ Energy2 s = sqr(5000.*GeV);
+ Energy2 t = -0.25*s;
+ Energy2 u = -0.75*s;
+ testEvolution(s,t,u);
- // cerr << "testing got here ???\n";
+
+
+
// cerr << subProcess() << "\n";
// cerr << *subProcess() << "\n";
// cerr << subProcess()->outgoing()[0] << *subProcess()->outgoing()[0] << "\n";
// cerr << subProcess()->outgoing()[0]->spinInfo() << "\n";
// cerr << subProcess()->outgoing()[0]->spinInfo()->productionVertex() << "\n";
assert(false);
staticEWCouplings_ = tEWCouplingsPtr();
}
+
+void ElectroWeakReweighter::testEvolution(Energy2 s,Energy2 t, Energy2 u) const {
+ Energy highScale = sqrt(s);
+ Energy ewScale = coupling()->mZ();
+ Energy lowScale = 50.0*GeV;
+ for (unsigned int i=0; i<45;++i) {
+ EWProcess::Process process = (EWProcess::Process)i;
+ cerr << "process " << process << "\n";
+ // result for all EW and QCD SCET contributions:
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highMatch_val
+ = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,true,true);
+ boost::numeric::ublas::matrix<Complex> highRunning_val
+ = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process);
+ boost::numeric::ublas::matrix<Complex> ewMatch_val =
+ ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true);
+ boost::numeric::ublas::matrix<Complex> lowRunning_val =
+ softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process);
+ boost::numeric::ublas::matrix<Complex> collinearHighRunning_val =
+ collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false);
+ boost::numeric::ublas::matrix<Complex> collinearEWMatch_val =
+ collinearSudakov_->electroWeakMatching(ewScale,s,process,true);
+ boost::numeric::ublas::matrix<Complex> collinearLowRunning_val =
+ collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process);
+ boost::numeric::ublas::matrix<Complex> lowMatchTemp_val =
+ boost::numeric::ublas::zero_matrix<Complex>(ewMatch_val.size1(),ewMatch_val.size2());
+ for (unsigned int ii=0; ii<ewMatch_val.size1(); ++ii) {
+ for (unsigned int jj=0; jj<ewMatch_val.size2(); ++jj) {
+ lowMatchTemp_val(ii,jj) = collinearEWMatch_val(ii,jj)*ewMatch_val(ii,jj);
+ }
+ }
+ boost::numeric::ublas::matrix<Complex> temp(highRunning_val.size1(),collinearHighRunning_val.size2());
+ boost::numeric::ublas::axpy_prod(highRunning_val,collinearHighRunning_val,temp);
+ boost::numeric::ublas::matrix<Complex> temp2(collinearLowRunning_val.size1(),lowRunning_val.size2());
+ boost::numeric::ublas::axpy_prod(collinearLowRunning_val,lowRunning_val,temp2);
+ boost::numeric::ublas::matrix<Complex> temp3(temp2.size1(),lowMatchTemp_val.size2());
+ boost::numeric::ublas::axpy_prod(temp2,lowMatchTemp_val,temp3);
+ temp2.resize(temp3.size1(),temp.size2());
+ boost::numeric::ublas::axpy_prod(temp3,temp,temp2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > result(temp2.size1(),highMatch_val.size2());
+ axpy_prod_local(temp2,highMatch_val,result);
+ for(unsigned int ix=0;ix<result.size1();++ix) {
+ for(unsigned int iy=0;iy<result.size2();++iy) {
+ cerr << s*result(ix,iy) << " ";
+ }
+ cerr << "\n";
+ }
+ }
+}
diff --git a/MatrixElement/EW/ElectroWeakReweighter.h b/MatrixElement/EW/ElectroWeakReweighter.h
--- a/MatrixElement/EW/ElectroWeakReweighter.h
+++ b/MatrixElement/EW/ElectroWeakReweighter.h
@@ -1,127 +1,140 @@
// -*- C++ -*-
#ifndef Herwig_ElectroWeakReweighter_H
#define Herwig_ElectroWeakReweighter_H
//
// This is the declaration of the ElectroWeakReweighter class.
//
#include "ThePEG/MatrixElement/ReweightBase.h"
#include "EWCouplings.h"
#include "CollinearSudakov.h"
+#include "SoftSudakov.h"
namespace Herwig {
using namespace ThePEG;
/**
* The ElectroWeakReweighter class.
*
* @see \ref ElectroWeakReweighterInterfaces "The interfaces"
* defined for ElectroWeakReweighter.
*/
class ElectroWeakReweighter: public ReweightBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ElectroWeakReweighter();
/**
* The destructor.
*/
virtual ~ElectroWeakReweighter();
//@}
public:
/**
* Return the weight for the kinematical configuation provided by
* the assigned XComb object (in the LastXCombInfo base class).
*/
virtual double weight() const;
/**
*
*/
static tEWCouplingsPtr coupling() {
assert(staticEWCouplings_);
return staticEWCouplings_;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
+ /**
+ * Check the evolution for a fixed s,t,u
+ */
+ void testEvolution(Energy2 s,Energy2 t, Energy2 u) const;
+
+protected:
+
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ElectroWeakReweighter & operator=(const ElectroWeakReweighter &);
private:
/**
* The Electroweak Couplings
*/
EWCouplingsPtr EWCouplings_;
/**
* The Collinear Sudakov
*/
CollinearSudakovPtr collinearSudakov_;
/**
+ * The Soft Sudakov
+ */
+ SoftSudakovPtr softSudakov_;
+
+ /**
* The couplings to allow global access
*/
static tEWCouplingsPtr staticEWCouplings_;
};
}
#endif /* Herwig_ElectroWeakReweighter_H */
diff --git a/MatrixElement/EW/Makefile.am b/MatrixElement/EW/Makefile.am
--- a/MatrixElement/EW/Makefile.am
+++ b/MatrixElement/EW/Makefile.am
@@ -1,9 +1,10 @@
pkglib_LTLIBRARIES = HwMEEW.la
HwMEEW_la_SOURCES = EWProcess.h GroupInvariants.h GroupInvariants.cc \
ElectroWeakReweigter.h ElectroWeakReweighter.cc \
SoftSudakov.h SoftSudakov.cc \
CollinearSudakov.h CollinearSudakov.cc \
HighEnergyMatching.h HighEnergyMatching.cc \
ElectroWeakMatching.h ElectroWeakMatching.cc \
-EWCouplings.h EWCouplings.fh EWCouplings.cc
+EWCouplings.h EWCouplings.fh EWCouplings.cc \
+expm-1.hpp
HwMEEW_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
diff --git a/MatrixElement/EW/SoftSudakov.cc b/MatrixElement/EW/SoftSudakov.cc
--- a/MatrixElement/EW/SoftSudakov.cc
+++ b/MatrixElement/EW/SoftSudakov.cc
@@ -1,881 +1,885 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SoftSudakov class.
//
#include "SoftSudakov.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 "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
+#include "expm-1.h"
using namespace Herwig;
using namespace GroupInvariants;
-SoftSudakov::SoftSudakov() : K_ORDER_(3) {}
+SoftSudakov::SoftSudakov() : K_ORDER_(3), integrator_(0.,1e-5,1000) {}
SoftSudakov::~SoftSudakov() {}
IBPtr SoftSudakov::clone() const {
return new_ptr(*this);
}
IBPtr SoftSudakov::fullclone() const {
return new_ptr(*this);
}
void SoftSudakov::persistentOutput(PersistentOStream & os) const {
os << K_ORDER_;
}
void SoftSudakov::persistentInput(PersistentIStream & is, int) {
is >> K_ORDER_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SoftSudakov,Interfaced>
describeHerwigSoftSudakov("Herwig::SoftSudakov", "HwMEEW.so");
void SoftSudakov::Init() {
static ClassDocumentation<SoftSudakov> documentation
("The SoftSudakov class implements the soft EW Sudakov");
}
InvEnergy SoftSudakov::operator ()(Energy mu) const {
// Include K-factor Contributions (Cusps):
GaugeContributions cusp = cuspContributions(mu,K_ORDER_,high_);
Complex gamma = cusp.SU3*G3_(row_,col_) + cusp.SU2*G2_(row_,col_) + cusp.U1*G1_(row_,col_);
if (real_) {
return gamma.real()/mu;
}
else {
return gamma.imag()/mu;
}
}
boost::numeric::ublas::matrix<Complex>
SoftSudakov::evaluateSoft(boost::numeric::ublas::matrix<Complex> & G3,
boost::numeric::ublas::matrix<Complex> & G2,
boost::numeric::ublas::matrix<Complex> & G1,
Energy mu_h, Energy mu_l, bool high) {
assert( G3.size1() == G2.size1() && G3.size1() == G1.size1() &&
G3.size2() == G2.size2() && G3.size2() == G1.size2() &&
G3.size1() == G3.size2());
G3_ = G3;
G2_ = G2;
G1_ = G1;
high_ = high;
unsigned int NN = G3_.size1();
// gamma is the matrix to be numerically integrated to run the coefficients.
boost::numeric::ublas::matrix<Complex> gamma(NN,NN);
for(row_=0;row_<NN;++row_) {
for(col_=0;col_<NN;++col_) {
- real_ = true;
- gamma(row_,col_).real(integrator_.value(*this,mu_h,mu_l));
- real_ = false;
- gamma(row_,col_).imag(integrator_.value(*this,mu_h,mu_l));
+ if(G3_(row_,col_) == 0. && G2_(row_,col_) == 0. && G1_(row_,col_) == 0.) {
+ gamma(row_,col_) = 0.;
+ }
+ else {
+ real_ = true;
+ gamma(row_,col_).real(integrator_.value(*this,mu_h,mu_l));
+ real_ = false;
+ gamma(row_,col_).imag(integrator_.value(*this,mu_h,mu_l));
+ }
}
}
// Resummed:
- //return gamma.exp();
- // Fixed Order:
- return boost::numeric::ublas::identity_matrix<Complex>(NN,NN) + gamma;
+ return boost::numeric::ublas::expm_pad(gamma,7);
}
boost::numeric::ublas::matrix<Complex>
SoftSudakov::lowEnergyRunning(Energy EWScale, Energy lowScale,
Energy2 s, Energy2 t, Energy2 u,
Herwig::EWProcess::Process process) {
using namespace EWProcess;
using Constants::pi;
static const Complex I(0,1.0);
Complex T = getT(s,t), U = getU(s,u);
boost::numeric::ublas::matrix<Complex> G1, G2, G3;
unsigned int numBrokenGauge;
switch (process) {
case QQQQ:
case QQQQiden:
case QtQtQQ:
{
numBrokenGauge = 12;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
for (unsigned int i=0; i<numBrokenGauge/2; i++) {
G3(i,i) += gam3(0,0);
G3(i,i+6) += gam3(0,1);
G3(i+6,i) += gam3(1,0);
G3(i+6,i+6) += gam3(1,1);
}
G1(0,0) = G1(6,6) = Gamma1(2.0/3.0,2.0/3.0,2.0/3.0,2.0/3.0,T,U);
G1(1,1) = G1(7,7) = Gamma1(-1.0/3.0,-1.0/3.0,2.0/3.0,2.0/3.0,T,U);
G1(2,2) = G1(8,8) = Gamma1(2.0/3.0,2.0/3.0,-1.0/3.0,-1.0/3.0,T,U);
G1(3,3) = G1(9,9) = Gamma1(-1.0/3.0,-1.0/3.0,-1.0/3.0,-1.0/3.0,T,U);
G1(4,4) = G1(10,10) = Gamma1(-1.0/3.0,2.0/3.0,2.0/3.0,-1.0/3.0,T,U);
G1(5,5) = G1(11,11) = Gamma1(2.0/3.0,-1.0/3.0,-1.0/3.0,2.0/3.0,T,U);
}
break;
case QQUU:
case QtQtUU:
case QQtRtR:
{
numBrokenGauge = 4;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
for (unsigned int i=0; i<numBrokenGauge/2; i++) {
G3(i,i) += gam3(0,0);
G3(i,i+2) += gam3(0,1);
G3(i+2,i) += gam3(1,0);
G3(i+2,i+2) += gam3(1,1);
}
G1(0,0) = G1(2,2) = Gamma1(2.0/3.0,2.0/3.0,T,U);
G1(1,1) = G1(3,3) = Gamma1(2.0/3.0,-1.0/3.0,T,U);
}
break;
case QQDD:
case QtQtDD:
{
numBrokenGauge = 4;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
for (unsigned int i=0; i<numBrokenGauge/2; i++) {
G3(i,i) += gam3(0,0);
G3(i,i+2) += gam3(0,1);
G3(i+2,i) += gam3(1,0);
G3(i+2,i+2) += gam3(1,1);
}
G1(0,0) = G1(2,2) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
G1(1,1) = G1(3,3) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
}
break;
case QQLL:
{
numBrokenGauge = 6;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(2.0/3.0,2.0/3.0,0.0,0.0,T,U);
G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,0.0,0.0,T,U);
G1(2,2) = Gamma1(2.0/3.0,2.0/3.0,-1.0,-1.0,T,U);
G1(3,3) = Gamma1(-1.0/3.0,-1.0/3.0,-1.0,-1.0,T,U);
G1(4,4) = Gamma1(-1.0/3.0,2.0/3.0,0.0,-1.0,T,U);
G1(5,5) = Gamma1(2.0/3.0,-1.0/3.0,-1.0,0.0,T,U);
}
break;
case QQEE:
{
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1 *= 0.0; G2 *= 0.0; G3 *= 0.0;
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<2; i++) {
G3(i,i) += gam3s;
}
G1(0,0) = Gamma1(2.0/3.0,-1.0,T,U);
G1(1,1) = Gamma1(-1.0/3.0,-1.0,T,U);
}
break;
case UUUU:
case UUUUiden:
case tRtRUU:
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1 *= 0.0; G2 *= 0.0; G3 *= 0.0;
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,2.0/3.0,T,U);
break;
case UUDD:
case tRtRDD:
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
break;
case UULL:
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3(0,0) = G3(1,1) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(2.0/3.0,0.0,T,U);
G1(1,1) = Gamma1(2.0/3.0,-1.0,T,U);
break;
case UUEE:
numBrokenGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(2.0/3.0,-1.0,T,U);
break;
case DDDD:
case DDDDiden:
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
break;
case DDLL:
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3(0,0) = G3(1,1) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/3.0,0.0,T,U);
G1(1,1) = Gamma1(-1.0/3.0,-1.0,T,U);
break;
case DDEE:
numBrokenGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/3.0,-1.0,T,U);
break;
case LLLL:
case LLLLiden:
numBrokenGauge = 6;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(0.0,0.0,0.0,0.0,T,U);
G1(1,1) = Gamma1(-1.0,-1.0,0.0,0.0,T,U);
G1(2,2) = Gamma1(0.0,0.0,-1.0,-1.0,T,U);
G1(3,3) = Gamma1(-1.0,-1.0,-1.0,-1.0,T,U);
G1(4,4) = Gamma1(-1.0,0.0,0.0,-1.0,T,U);
G1(5,5) = Gamma1(0.0,-1.0,-1.0,0.0,T,U);
break;
case LLEE:
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(0.0,-1.0,T,U);
G1(1,1) = Gamma1(-1.0,-1.0,T,U);
break;
case EEEE:
case EEEEiden:
numBrokenGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(-1.0,-1.0,T,U);
break;
case QQWW:
{
numBrokenGauge = 20;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(2./3.,2./3.,-1.,-1.,T,U);
G1(1,1) = Gamma1(2./3.,2./3.,1.,1.,T,U);
G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(4,4) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(5,5) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(6,6) = Gamma1(-1./3.,-1./3.,-1.,-1.,T,U);
G1(7,7) = Gamma1(-1./3.,-1./3.,1.,1.,T,U);
G1(8,8) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(9,9) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(10,10) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(11,11) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(12,12) = Gamma1(-1./3.,2./3.,0.,-1.,T,U);
G1(13,13) = Gamma1(-1./3.,2./3.,0.,-1.,T,U);
G1(14,14) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(15,15) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(16,16) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
G1(17,17) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
G1(18,18) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
G1(19,19) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
}
break;
case QQPhiPhi:
{
numBrokenGauge = 14;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(2./3.,2./3.,1.,1.,T,U);
G1(1,1) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(4,4) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(5,5) = Gamma1(-1./3.,-1./3.,1.,1.,T,U);
G1(6,6) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(7,7) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(8,8) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(9,9) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(10,10) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(11,11) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(12,12) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
G1(13,13) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
}
break;
case QQWG:
numBrokenGauge = 6;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = -17.0/6.0*I*pi + 3.0/2.0*(U+T);
}
G1(0,0) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(1,1) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(4,4) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(5,5) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
break;
case QQBG:
numBrokenGauge = 4;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
}
G1(0,0) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(1,1) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(2,2) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(3,3) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
break;
case QQGG:
case QtQtGG:
{
numBrokenGauge = 6;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
boost::numeric::ublas::matrix<Complex> gam3g = Gamma3g(U,T);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
G3(ix ,iy ) = gam3g(ix,iy);
G3(ix+3,iy+3) = gam3g(ix,iy);
}
}
G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(2./3.,0.,T,U);
G1(3,3) = G1(4,4) = G1(5,5) = Gamma1(-1./3.,0.,T,U);
}
break;
case LLWW:
numBrokenGauge = 20;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(0.,0.,-1.,-1.,T,U);
G1(1,1) = Gamma1(0.,0.,1.,1.,T,U);
G1(2,2) = Gamma1(0.,0.,0.,0.,T,U);
G1(3,3) = Gamma1(0.,0.,0.,0.,T,U);
G1(4,4) = Gamma1(0.,0.,0.,0.,T,U);
G1(5,5) = Gamma1(0.,0.,0.,0.,T,U);
G1(6,6) = Gamma1(-1.,-1.,-1.,-1.,T,U);
G1(7,7) = Gamma1(-1.,-1.,1.,1.,T,U);
G1(8,8) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(9,9) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(10,10) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(11,11) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(12,12) = Gamma1(-1.,0.,0.,-1.,T,U);
G1(13,13) = Gamma1(-1.,0.,0.,-1.,T,U);
G1(14,14) = Gamma1(-1.,0.,1.,0.,T,U);
G1(15,15) = Gamma1(-1.,0.,1.,0.,T,U);
G1(16,16) = Gamma1(0.,-1.,-1.,0.,T,U);
G1(17,17) = Gamma1(0.,-1.,-1.,0.,T,U);
G1(18,18) = Gamma1(0.,-1.,0.,1.,T,U);
G1(19,19) = Gamma1(0.,-1.,0.,1.,T,U);
break;
case LLPhiPhi:
numBrokenGauge = 14;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1 *= 0.0; G2 *= 0.0; G3 *= 0.0;
G1(0,0) = Gamma1(0.,0.,1.,1.,T,U);
G1(1,1) = Gamma1(0.,0.,0.,0.,T,U);
G1(2,2) = Gamma1(0.,0.,0.,0.,T,U);
G1(3,3) = Gamma1(0.,0.,0.,0.,T,U);
G1(4,4) = Gamma1(0.,0.,0.,0.,T,U);
G1(5,5) = Gamma1(-1.,-1.,1.,1.,T,U);
G1(6,6) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(7,7) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(8,8) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(9,9) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(10,10) = Gamma1(-1.,0.,1.,0.,T,U);
G1(11,11) = Gamma1(-1.,0.,1.,0.,T,U);
G1(12,12) = Gamma1(0.,-1.,0.,1.,T,U);
G1(13,13) = Gamma1(0.,-1.,0.,1.,T,U);
break;
case UUBB:
{
numBrokenGauge = 4;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
G1(i,i) = Gamma1(2./3.,0.,T,U);
}
}
break;
case UUPhiPhi:
{
numBrokenGauge = 5;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(2.0/3.0,2.0/3.0,1.,1.,T,U);
G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(2./3.,0.,T,U);
}
break;
case UUBG:
{
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam1 = Gamma1(2./3.,0.,T,U);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
G1(i,i) = gam1;
}
}
break;
case UUGG:
case tRtRGG:
numBrokenGauge = 3;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = Gamma3g(U,T);
G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(2./3.,0.,T,U);
break;
case DDBB:
{
numBrokenGauge = 4;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
Complex gam1 = Gamma1(-1./3.,0.,T,U);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
G1(i,i) = gam1;
}
}
break;
case DDPhiPhi:
{
numBrokenGauge = 5;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(-1.0/3.0,-1.0/3.0,1.,1.,T,U);
G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(-1./3.,0.,T,U);
}
break;
case DDBG:
numBrokenGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
}
G1(0,0) = G1(1,1) = Gamma1(-1./3.,0.,T,U);
break;
case DDGG:
numBrokenGauge = 3;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = Gamma3g(U,T);
G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(-1./3.,0.,T,U);
break;
case EEBB:
{
numBrokenGauge = 4;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam1 = Gamma1(-1.,0.,T,U);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G1(i,i) = gam1;
};
}
break;
case EEPhiPhi:
numBrokenGauge = 5;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(-1.,1.,T,U);
G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(-1.,0.,T,U);
break;
default:
assert(false);
}
-
+ // return the answer
if (EWScale==lowScale) {
return boost::numeric::ublas::identity_matrix<Complex>(G1.size1());
}
else {
return evaluateSoft(G3,G2,G1,EWScale,lowScale,false);
}
}
boost::numeric::ublas::matrix<Complex>
SoftSudakov::highEnergyRunning(Energy highScale, Energy EWScale,
Energy2 s, Energy2 t, Energy2 u,
Herwig::EWProcess::Process process) {
using namespace EWProcess;
using Constants::pi;
static const Complex I(0,1.0);
Complex T = getT(s,t), U = getU(s,u);
boost::numeric::ublas::matrix<Complex> G1,G2,G3;
unsigned int numGauge;
switch (process) {
case QQQQ:
case QQQQiden:
case QtQtQQ:
{
numGauge = 4;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
G3(0,0) += gam3(0,0);
G3(0,2) += gam3(0,1);
G3(2,0) += gam3(1,0);
G3(2,2) += gam3(1,1);
G3(1,1) += gam3(0,0);
G3(1,3) += gam3(0,1);
G3(3,1) += gam3(1,0);
G3(3,3) += gam3(1,1);
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);
G1(0,0) = G1(1,1) = G1(2,2) = G1(3,3) = Gamma1(1.0/6.0,1.0/6.0,T,U);
}
break;
case QQUU:
case QtQtUU:
case QQtRtR:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G2 = Gamma2Singlet();
G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,1.0/6.0,T,U);
break;
case QQDD:
case QtQtDD:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G2 = Gamma2Singlet();
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,1.0/6.0,T,U);
break;
case QQLL:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3Singlet();
G2 = Gamma2(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/2.0,1.0/6.0,T,U);
break;
case QQEE:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(-1.0,1.0/6.0,T,U);
break;
case UUUU:
case UUUUiden:
case tRtRUU:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,2.0/3.0,T,U);
break;
case UUDD:
case tRtRDD:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
break;
case UULL:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/2.0,2.0/3.0,T,U);
break;
case UUEE:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0,2.0/3.0,T,U);
break;
case DDDD:
case DDDDiden:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
break;
case DDLL:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/2.0,-1.0/3.0,T,U);
break;
case DDEE:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0,-1.0/3.0,T,U);
break;
case LLLL:
case LLLLiden:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = Gamma2(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/2.0,-1.0/2.0,T,U);
break;
case LLEE:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(-1.0,-1.0/2.0,T,U);
break;
case EEEE:
case EEEEiden:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G1(0,0) = Gamma1(-1.0,-1.0,T,U);
break;
case QQWW:
{
numGauge = 5;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<5; i++) {
G3(i,i) = gam3s;
G1(i,i) = Gamma1(1.0/6.0);
}
G2 = Gamma2w(U,T);
}
break;
case QQPhiPhi:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3Singlet();
G2 = Gamma2(U,T);
G1(0,0) = G1(1,1) = Gamma1(1.0/2.0,1.0/6.0,T,U);
break;
case QQWG:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = -17.0/6.0*I*pi + 3.0/2.0*(U+T);
G2(0,0) = -7.0/4.0*I*pi + (U+T);
G1(0,0) = Gamma1(1.0/6.0);
break;
case QQBG:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
G2(0,0) = -3.0/4.0*I*pi;
G1(0,0) = Gamma1(1.0/6.0);
break;
case QQGG:
case QtQtGG:
{
numGauge = 3;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = Gamma3g(T,U);
+ G3 = Gamma3g(U,T);
Complex gam2s = Gamma2Singlet()(0,0);
for (unsigned int i=0; i<3; i++) {
G2(i,i) = gam2s;
G1(i,i) = Gamma1(1.0/6.0);
}
}
break;
case LLWW:
numGauge = 5;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
for (unsigned int i=0; i<5; i++) {
G1(i,i) = Gamma1(-1.0/2.0);
}
G2 = Gamma2w(U,T);
break;
case LLPhiPhi:
numGauge = 2;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = Gamma2(U,T);
G1(0,0) = G1(1,1) = Gamma1(1.0/2.0,-1.0/2.0,T,U);
break;
case UUBB:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(2.0/3.0);
break;
case UUPhiPhi:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(1.0/2.0,2.0/3.0,T,U);
break;
case UUBG:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
G1(0,0) = Gamma1(2.0/3.0);
break;
case UUGG:
case tRtRGG:
numGauge = 3;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3g(U,T);
for (unsigned int i=0; i<3; i++) {
G1(i,i) = Gamma1(2.0/3.0);
}
break;
case DDBB:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/3.0);
break;
case DDPhiPhi:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(1.0/2.0,-1.0/3.0,T,U);
break;
case DDBG:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
G1(0,0) = Gamma1(-1.0/3.0);
break;
case DDGG:
numGauge = 3;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3g(U,T);
for (unsigned int i=0; i<3; i++) {
G1(i,i) = Gamma1(-1.0/3.0);
}
break;
case EEBB:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G1(0,0) = Gamma1(-1.0);
break;
case EEPhiPhi:
numGauge = 1;
G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(1.0/2.0,-1.0,T,U);
break;
default:
assert(false);
}
return evaluateSoft(G3,G2,G1,highScale,EWScale,true);
}
diff --git a/MatrixElement/EW/expm-1.h b/MatrixElement/EW/expm-1.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/expm-1.h
@@ -0,0 +1,147 @@
+//
+// Copyright (c) 2007
+// Tsai, Dung-Bang
+// National Taiwan University, Department of Physics
+//
+// E-Mail : dbtsai (at) gmail.com
+// Begine : 2007/11/20
+// Last modify : 2007/11/22
+// Version : v0.1
+//
+// EXPGM_PAD computes the matrix exponential exp(H) for general matrixs,
+// including complex and real matrixs using the irreducible (p,p) degree
+// rational Pade approximation to the exponential
+// exp(z) = r(z)=(+/-)( I+2*(Q(z)/P(z))).
+//
+// Usage :
+//
+// U = expm_pad(H)
+// U = expm_pad(H, p)
+//
+// where p is internally set to 6 (recommended and gererally satisfactory).
+//
+// See also MATLAB supplied functions, EXPM and EXPM1.
+//
+// Reference :
+// EXPOKIT, Software Package for Computing Matrix Exponentials.
+// ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
+//
+// Permission to use, copy, modify, distribute and sell this software
+// and its documentation for any purpose is hereby granted without fee,
+// provided that the above copyright notice appear in all copies and
+// that both that copyright notice and this permission notice appear
+// in supporting documentation. The authors make no representations
+// about the suitability of this software for any purpose.
+// It is provided "as is" without express or implied warranty.
+//
+
+#ifndef _BOOST_UBLAS_EXPM_
+#define _BOOST_UBLAS_EXPM_
+#include <complex>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+
+namespace boost { namespace numeric { namespace ublas {
+
+template<typename MATRIX> MATRIX expm_pad(const MATRIX &H, const int p = 6) {
+ typedef typename MATRIX::value_type value_type;
+ typedef typename MATRIX::size_type size_type;
+ typedef double real_value_type; // Correct me. Need to modify.
+ assert(H.size1() == H.size2());
+ const size_type n = H.size1();
+ const identity_matrix<value_type> I(n);
+ matrix<value_type> U(n,n),H2(n,n),P(n,n),Q(n,n);
+ real_value_type norm = 0.0;
+
+ // Calcuate Pade coefficients (1-based instead of 0-based as in the c vector)
+ vector<real_value_type> c(p+2);
+ c(1)=1;
+ for(size_type i = 1; i <= p; ++i)
+ c(i+1) = c(i) * ((p + 1.0 - i)/(i * (2.0 * p + 1 - i)));
+ // Calcuate the infinty norm of H, which is defined as the largest row sum of a matrix
+ for(size_type i=0; i<n; ++i)
+ {
+ real_value_type temp = 0.0;
+ for(size_type j=0;j<n;j++)
+ temp += std::abs<real_value_type>(H(j,i)); // Correct me, if H is complex, can I use that abs?
+ norm = std::max<real_value_type>(norm, temp);
+ }
+ if (norm == 0.0)
+ {
+ boost::throw_exception(boost::numeric::ublas::bad_argument());
+ std::cerr<<"Error! Null input in the routine EXPM_PAD.\n";
+ exit(0);
+ }
+ // Scaling, seek s such that || H*2^(-s) || < 1/2, and set scale = 2^(-s)
+ int s = 0;
+ real_value_type scale = 1.0;
+ if(norm > 0.5) {
+ s = std::max<int>(0, static_cast<int>((log(norm) / log(2.0) + 2.0)));
+ scale /= static_cast<real_value_type>(std::pow(2.0, s));
+ U.assign(scale * H); // Here U is used as temp value due to that H is const
+ }
+ else
+ U.assign(H);
+ // Horner evaluation of the irreducible fraction, see the following ref above.
+ // Initialise P (numerator) and Q (denominator)
+ H2.assign( prod(U, U) );
+ Q.assign( c(p+1)*I );
+ P.assign( c(p)*I );
+ size_type odd = 1;
+ for( size_type k = p - 1; k > 0; --k)
+ {
+ if( odd == 1)
+ {
+ Q = ( prod(Q, H2) + c(k) * I );
+ }
+ else
+ {
+ P = ( prod(P, H2) + c(k) * I );
+ }
+ odd = 1 - odd;
+ }
+ if( odd == 1)
+ {
+ Q = ( prod(Q, U) );
+ Q -= P ;
+ //U.assign( -(I + 2*(Q\P)));
+ }
+ else
+ {
+ P = (prod(P, U));
+ Q -= P;
+ //U.assign( I + 2*(Q\P));
+ }
+ // In origine expokit package, they use lapack ZGESV to obtain inverse matrix,
+ // and in that ZGESV routine, it uses LU decomposition for obtaing inverse matrix.
+ // Since in ublas, there is no matrix inversion template, I simply use the build-in
+ // LU decompostion package in ublas, and back substitute by myself.
+ //
+ //////////////// Implement Matrix Inversion ///////////////////////
+ permutation_matrix<size_type> pm(n);
+ int res = lu_factorize(Q, pm);
+ if( res != 0)
+ {
+ std::cerr << "Error in the matrix inversion in template expm_pad.\n";
+ exit(0);
+ }
+ H2 = I; // H2 is not needed anymore, so it is temporary used as identity matrix for substituting.
+
+ lu_substitute(Q, pm, H2);
+ if( odd == 1)
+ U.assign( -(I + 2.0 * prod(H2, P)));
+ else
+ U.assign( I + 2.0 * prod(H2, P));
+ // Squaring
+ for(size_t i = 0; i < s; ++i)
+ {
+ U = (prod(U,U));
+ }
+ return U;
+ }
+
+}}}
+
+
+#endif
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:13 PM (15 h, 4 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023599
Default Alt Text
(162 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment