Page MenuHomeHEPForge
Files F8309833

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/EW/EWCouplings.cc b/MatrixElement/EW/EWCouplings.cc
--- a/MatrixElement/EW/EWCouplings.cc
+++ b/MatrixElement/EW/EWCouplings.cc
@@ -1,629 +1,734 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the EWCouplings class.
//
#include "EWCouplings.h"
#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <boost/numeric/ublas/operation.hpp>
using namespace Herwig;
namespace {
Complex trace(boost::numeric::ublas::matrix<Complex> M) {
assert(M.size1()==M.size2());
Complex output(0.);
for(unsigned int ix=0;ix<M.size1();++ix)
output += M(ix,ix);
return output;
}
}
-EWCouplings::EWCouplings(unsigned int loops, unsigned int steps, Energy highScale,
+EWCouplings::EWCouplings(unsigned int loops, unsigned int steps,
Energy lowScale)
- : ewScale_(91.1876*GeV), highScale_(highScale), lowScale_(lowScale),
+ : ewScale_(91.1876*GeV), highScale_(14.*TeV), lowScale_(lowScale),
includeSU3_(true), includeEW_(true), initialized_(false), massChoice_(false),
mZ_(91.1876*GeV), mW_(80.399*GeV),
mT_(173.1*GeV), // 179.08045 (should be this?)
mH_(125.0*GeV),
- loops_(loops), highSteps_(steps), lowSteps_(steps)
+ loops_(loops), highSteps_(steps), lowSteps_(steps),
+ a1MZ_(0.01017054),a2MZ_(0.03378168),aSMZ_(0.1176),lambdat_(0.991172)
{}
void EWCouplings::initialize() {
using Constants::pi;
if(initialized_) return;
initialized_ = true;
+ highScale_=generator()->maximumCMEnergy();
// set the particle masses
if(massChoice_) {
mZ_ = getParticleData(ParticleID::Z0 )->mass();
mW_ = getParticleData(ParticleID::Wplus)->mass();
mT_ = getParticleData(ParticleID::t )->mass();
mH_ = getParticleData(ParticleID::h0 )->mass();
}
// logs of scales
double logEWScale = log(ewScale_/GeV);
double logHighScale = log(highScale_/GeV);
// step size
double stepsize = (logHighScale - logEWScale)/double(highSteps_);
// Initialize parameters at the ewScale
// 32 parameters, mostly zero due massless quarks
unsigned int N = 32;
vector<Complex> y(N,0.), dydx(N,0.), yout(N,0.);
initializeCouplings(y);
double x = logEWScale;
derivatives(x,y,dydx);
// energy scale + 6 parameters: g1,g2,g3,y_t,lambda,vev
table_ = boost::numeric::ublas::matrix<double>(highSteps_+1,7);
table_(0,0) = logEWScale;
for (unsigned int i=1; i<N; i++) {
if (i <=3) table_(0,i ) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
if (i >=29) table_(0,i-25) = y[i].real();
}
int counter = 1;
// Use 4th order runge-kutta to integrate to highScale
int steps = highSteps_;
while (steps > 0) {
RK4(y,dydx,x,stepsize,yout);
// Advance x and calculate derivatives at new starting point
for(unsigned int j=0; j<N; j++) {
y[j] = yout[j];
}
x += stepsize;
derivatives(x,y,dydx);
table_(counter,0) = x;
for (unsigned int i=1; i<N; i++) {
if (i<=3 ) table_(counter,i ) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
if (i>=29) table_(counter,i-25) = y[i].real();
}
steps--;
counter++;
}
// Initialize couplings at mu < 91.1876 GeV
initializeLow();
}
EWCouplings::~EWCouplings() {}
IBPtr EWCouplings::clone() const {
return new_ptr(*this);
}
IBPtr EWCouplings::fullclone() const {
return new_ptr(*this);
}
void EWCouplings::persistentOutput(PersistentOStream & os) const {
- os << ounit(ewScale_,GeV) << ounit(highScale_,GeV) << ounit(lowScale_,GeV)
+ os << ounit(ewScale_,GeV) << ounit(highScale_,GeV) << ounit(lowScale_,GeV)
<< includeSU3_ << includeEW_ << ounit(mZ_,GeV) << ounit(mW_,GeV)
- << ounit(mT_,GeV) << ounit(mH_,GeV) << massChoice_ << initialized_;
+ << ounit(mT_,GeV) << ounit(mH_,GeV) << massChoice_ << initialized_
+ << loops_ << highSteps_ << lowSteps_
+ << a1MZ_ << a2MZ_ << aSMZ_ << lambdat_;
+ os << table_.size1() << table_.size2();
+ for(unsigned int ix=0;ix<table_.size1();++ix)
+ for(unsigned int iy=0;iy<table_.size2();++iy)
+ os << table_(ix,iy);
+ os << lowTable_.size1() << lowTable_.size2();
+ for(unsigned int ix=0;ix<lowTable_.size1();++ix)
+ for(unsigned int iy=0;iy<lowTable_.size2();++iy)
+ os << lowTable_(ix,iy);
}
void EWCouplings::persistentInput(PersistentIStream & is, int) {
- is >> iunit(ewScale_,GeV) >> iunit(highScale_,GeV) >> iunit(lowScale_,GeV)
+ is >> iunit(ewScale_,GeV) >> iunit(highScale_,GeV) >> iunit(lowScale_,GeV)
>> includeSU3_ >> includeEW_ >> iunit(mZ_,GeV) >> iunit(mW_,GeV)
- >> iunit(mT_,GeV) >> iunit(mH_,GeV) >> massChoice_ >> initialized_;
+ >> iunit(mT_,GeV) >> iunit(mH_,GeV) >> massChoice_ >> initialized_
+ >> loops_ >> highSteps_ >> lowSteps_
+ >> a1MZ_ >> a2MZ_ >> aSMZ_ >> lambdat_;
+ unsigned int size1,size2;
+ is >> size1 >> size2;
+ table_.resize(size1,size2);
+ for(unsigned int ix=0;ix<table_.size1();++ix)
+ for(unsigned int iy=0;iy<table_.size2();++iy)
+ is >> table_(ix,iy);
+ is >> size1 >> size2;
+ lowTable_.resize(size1,size2);
+ for(unsigned int ix=0;ix<lowTable_.size1();++ix)
+ for(unsigned int iy=0;iy<lowTable_.size2();++iy)
+ is >> lowTable_(ix,iy);
}
//The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<EWCouplings,Interfaced>
describeHerwigEWCouplings("Herwig::EWCouplings", "HwMEEW.so");
void EWCouplings::Init() {
static ClassDocumentation<EWCouplings> documentation
("The EWCouplings class implements");
static Switch<EWCouplings,bool> interfaceMassChoice
("MassChoice",
"Where to get the SM particle masses from",
&EWCouplings::massChoice_, false, false, false);
static SwitchOption interfaceMassChoiceLocal
(interfaceMassChoice,
"Local",
"Use local values",
false);
static SwitchOption interfaceMassChoiceParticleData
(interfaceMassChoice,
"ParticleData",
"Get the values from the ParticleData object",
true);
+ static Parameter<EWCouplings,Energy> interfaceEWScale
+ ("EWScale",
+ "The electroweak scale for matching between high and low energy running",
+ &EWCouplings::ewScale_, GeV, 91.1876*GeV, 10.0*GeV, 10000.0*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,Energy> interfaceLowScale
+ ("LowScale",
+ "The low energy scale at which to stop the running",
+ &EWCouplings::lowScale_, GeV, 10.*GeV, 1.0*GeV, 100.0*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,Energy> interfacemZ
+ ("mZ",
+ "The mass of the Z boson",
+ &EWCouplings::mZ_, 91.1876*GeV, 90.*GeV, 92.*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,Energy> interfacemW
+ ("mW",
+ "The mass of the W boson",
+ &EWCouplings::mW_, 80.399*GeV, 75.*GeV, 85.*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,Energy> interfacemT
+ ("mT",
+ "The mass of the top quark",
+ &EWCouplings::mT_, 173.1*GeV, 100.*GeV, 1000.*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,Energy> interfacemH
+ ("mH",
+ "The mass of the Higgs boson",
+ &EWCouplings::mH_, 125.*GeV, 100.*GeV, 1000.*GeV,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,unsigned int> interfaceLoops
+ ("Loops",
+ "The number of loops",
+ &EWCouplings::loops_, 2, 1, 3,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,unsigned int> interfaceHighSteps
+ ("HighSteps",
+ "The number of steps for the Runga-Kutta and interpolation of the couplings above mZ",
+ &EWCouplings::highSteps_, 200, 10, 1000000,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,unsigned int> interfaceLowSteps
+ ("LowSteps",
+ "The number of steps for the Runga-Kutta and interpolation of the couplings below mZ",
+ &EWCouplings::lowSteps_, 200, 10, 1000000,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,double> interfacea1MZ
+ ("Alpha1MZ",
+ "The value of a1(MZ)",
+ &EWCouplings::a1MZ_, 0.01017054, 0.0, 1.,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,double> interfacea2MZ
+ ("Alpha2MZ",
+ "The value of a2(MZ)",
+ &EWCouplings::a2MZ_, 0.03378168, 0.0, 1.,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,double> interfaceasMZ
+ ("AlphasMZ",
+ "The value of as(MZ)",
+ &EWCouplings::aSMZ_, 0.1176, 0.0, 1.,
+ false, false, Interface::limited);
+
+ static Parameter<EWCouplings,double> interfaceLambdaT
+ ("LambdaT",
+ "The top quark Yukawa at the matching scale",
+ &EWCouplings::lambdat_, 0.991172, 0.5, 2.0,
+ false, false, Interface::limited);
+
}
void EWCouplings::initializeLow() {
using Constants::pi;
// For scales less than ewScale, the only couplings calculated here are those of
// alpha_EW and alpha3:
double logEWScale = log(ewScale_ /GeV);
double loglowScale = log(lowScale_/GeV);
double stepsize = (loglowScale - logEWScale)/double(lowSteps_);
int steps = lowSteps_;
// Initialize parameters at the ewScale
unsigned int N=2; // Total # of parameters = 2
vector<Complex> y(N), dydx(N), yout(N);
for (unsigned int i=0; i<N; i++) {
y[i] = 0.0;
dydx[i] = 0.0;
yout[i] = 0.0;
}
double x = logEWScale;
// Initialize Couplings at the ewScale, including the One-Loop Threshold Matching:
double a1 = (3.0/5.0)*table_(0,1);
double a2 = table_(0,2);
double a3 = table_(0,3);
double aEM_ewScale = 1./(1./a1+1./a2-1./(6.*pi)*(1.-21.*log(mW()/ewScale_)+16./3.*log(mTatmZ()/ewScale_)));
double aS_ewScale = 1./(1./a3+1./(3.*pi)*log(ewScale_/mTatmZ()));
y[0] = sqrt(4.0*pi*aEM_ewScale);
y[1] = sqrt(4.0*pi*aS_ewScale);
derivatives(x,y,dydx);
// energy scale + 2 parameters: gEM,g3
lowTable_.resize(lowSteps_+1,3);
lowTable_(0,0) = logEWScale;
for (unsigned int i=1; i<=N; i++) {
lowTable_(0,i) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
}
int counter = 1;
// Use 4th order runge-kutta to integrate to highScale
while (steps > 0) {
RK4(y,dydx,x,stepsize,yout);
// Advance x and calculate derivatives at new starting point
for(unsigned int j=0; j<N; j++) {
y[j] = yout[j];
}
x += stepsize;
derivatives(x,y,dydx);
lowTable_(counter,0) = x;
for (unsigned int i=1; i<=N; i++) {
lowTable_(counter,i) = (y[i-1].real()*y[i-1].real())/(4.0*pi);
}
steps--;
counter++;
}
}
void EWCouplings::RK4(vector<Complex> &y, vector<Complex> &dydx,
const double x, const double h, vector<Complex> &yout) {
unsigned int n = y.size();
std::vector<Complex> dym(n), dyt(n), yt(n);
double hh = h*0.5;
double h6 = h/6.0;
double xh = x + hh;
const Complex I(0,1.0);
for(unsigned int i=0; i<n; i++) yt[i] = y[i] + hh*dydx[i];
derivatives(xh, yt, dyt);
for(unsigned int i=0; i<n; i++) yt[i] = y[i] + hh*dyt[i];
derivatives(xh, yt, dym);
for(unsigned int i=0; i<n; i++) {
yt[i] = y[i] + h*dym[i];
dym[i] += dyt[i];
}
derivatives(x+h, yt, dyt);
for(unsigned int i=0; i<n; i++) {
yout[i] = y[i] + h6*(dydx[i] + dyt[i] + 2.0*dym[i]);
}
}
void EWCouplings::initializeCouplings(vector<Complex> & y) {
+ using Constants::pi;
// \todo make these values parameters so they can be changed
InvEnergy2 gFermi = 1.16637*pow(10.0,-5)/GeV2;
Energy vev = 1.0/(sqrt(sqrt(2.0)*gFermi)); // vev = 246.221
- y[0] = 0.461531463; // g1 = Sqrt[5/3] * Sqrt[4*pi*a1] with a1(Mz) = 0.01017054
- y[1] = 0.651547066; // g2 = Sqrt[4*pi*a2] with a2(Mz) = 0.03378168
- y[2] = 1.215650108; // g3 = Sqrt[4*pi*as] with as(Mz) = 0.1176
+ y[0] = sqrt(5./3.*4.*pi*a1MZ_); // g1 = Sqrt[5/3] * Sqrt[4*pi*a1]
+ y[1] = sqrt(4.*pi*a2MZ_); // g2 = Sqrt[4*pi*a2]
+ y[2] = sqrt(4.*pi*aSMZ_); // g3 = Sqrt[4*pi*as]
// Note lambda_t = sqrt(2.0)*mt/vev only valid for mt(mt); need mt(mZ) here
// Top Yukawa lambda from Manohar
//Complex lambda_t = 1.02858;
// Top Yukawa lambda from Sascha
- double lambda_t = 0.991172;
+ double lambda_t =lambdat_;
// Quartic coupling lambda (need to multiply by a factor of 2 when accessing the quartic coupling)
double lambda = (mH_/vev)*(mH_/vev);
y[29] = lambda_t;
y[30] = lambda;
y[31] = vev/GeV;
}
void EWCouplings::derivatives(const double x, vector<Complex> & y,
vector<Complex> &dydx) {
// zero the output
for (unsigned int i=0; i<dydx.size(); i++) dydx[i]=0.0;
// low scale
if (y.size()==2) {
lowBetaGauge(x,y,dydx);
}
// high scale
else {
betaGauge(x,y,dydx);
betaYukawa(x,y,dydx);
betaHiggs(x,y,dydx);
}
}
void EWCouplings::betaGauge(const double x, vector<Complex> &y, vector<Complex> & dydx) {
using Constants::pi;
using boost::numeric::ublas::axpy_prod;
using boost::numeric::ublas::herm;
const Complex I(0,1.0);
// Yukawa
boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
Yuk_u(ix,iy) = y[21+3*ix+iy];
Yuk_d(ix,iy) = y[12+3*ix+iy];
Yuk_e(ix,iy) = y[ 3+3*ix+iy];
}
}
// gauge
boost::numeric::ublas::vector<Complex> gauge(3);
for(unsigned int l=0; l<3; l++) gauge[l] = y[l];
// Evaluate beta functions for gauge couplings
double Ng = 0.5*numberOfFlavours(x);
boost::numeric::ublas::vector<Complex> b(3), g2(3), gc(3), Cu(3), Cd(3), Ce(3);
boost::numeric::ublas::matrix<Complex> B1(3,3),B2(3,3), B3(3,3), B(3,3);
b[0] = -4.0/3.0*Ng - 1.0/10.0;
b[1] = 22.0/3.0 - 4.0/3.0*Ng - 1.0/6.0;
b[2] = 11.0 - 4.0/3.0*Ng;
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
B1(ix,iy) = 0.;
B2(ix,iy) = 0.;
B3(ix,iy) = 0.;
}
}
B1(1,1) = 136.0/3.0;
B1(2,2) = 102.0;
B2(0,0) = 19.0/15.0;
B2(0,1) = 1.0/5.0;
B2(0,2) = 11.0/30.0;
B2(1,0) = 3.0/5.0;
B2(1,1) = 49.0/3.0;
B2(1,2) = 3.0/2.0 ;
B2(2,0) = 44.0/15.0;
B2(2,1) = 4.0;
B2(2,2) = 76.0/3.0;
B3(0,0) = 9.0/50.0;
B3(0,1) = 3.0/10.0;
B3(1,0) = 9.0/10.0;
B3(1,1) = 13.0/6.0;
B = B1 - Ng*B2 - B3;
Cu[0] = 17.0/10.0;
Cu[1] = 3.0/2.0;
Cu[2] = 2.0;
Cd[0] = 1.0/2.0;
Cd[1] = 3.0/2.0;
Cd[2] = 2.0;
Ce[0] = 3.0/2.0;
Ce[1] = 1.0/2.0;
Ce[2] = 0.0;
for (int i=0; i<3; i++) {
g2(i) = pow(gauge(i),2);
}
// gc = trans(g2) * B
axpy_prod(g2, B, gc, true);
// compute the answer
if(loops_ >= 1) {
for(int l=0; l<3; l++) {
dydx[l] = -b(l)*pow(gauge(l),3)/(16.0*pow(pi,2));
}
if (loops_ >= 2) {
boost::numeric::ublas::matrix<Complex> temp(3,3);
axpy_prod(herm(Yuk_u),Yuk_u,temp);
Complex tr1 = trace(temp);
axpy_prod(herm(Yuk_d),Yuk_d,temp);
Complex tr2 = trace(temp);
axpy_prod(herm(Yuk_e),Yuk_e,temp);
Complex tr3 = trace(temp);
for(int l=0; l<3; l++) {
dydx[l] += -pow(gauge(l),3)/pow(16.0*pow(pi,2),2)*
(gc(l) + Cu(l)*tr1 + Cd(l)*tr2 + Ce(l)*tr3 );
}
}
}
}
void EWCouplings::lowBetaGauge(const double, vector<Complex> &y,
vector<Complex> &dydx) {
using Constants::pi;
const Complex I(0,1.0);
Complex e = y[0], g3 = y[1];
// Evaluate beta functions for gauge couplings
double Nu = 2.0, Nd = 3.0, Nl = 3.0;
if(loops_ >=1) {
dydx[0] += (16.0/9.0*Nu + 4.0/9.0*Nd + 4.0/3.0*Nl)*pow(e,3)/pow(4.0*pi,2);
dydx[1] += (2.0/3.0*(Nu+Nd)-11.0)*pow(g3,3)/pow(4.0*pi,2);
// Note this also includes the three-loop contribution for alpha_3
if (loops_ >= 2) {
dydx[0] += (64.0/27.0*Nu+4.0/27.0*Nd+4.0*Nl)*pow(e,5)/pow(4.0*pi,4) +
(64.0/9.0*Nu+16.0/9.0*Nd)*pow(e,3)*pow(g3,2)/pow(4.0*pi,4);
dydx[1] += (38.0/3.0*(Nu+Nd)-102.0)*pow(g3,5)/pow(4.0*pi,4) +
(8.0/9.0*Nu+2.0/9.0*Nd)*pow(g3,3)*pow(e,2)/pow(4.0*pi,4) +
(5033.0/18.0*(Nu+Nd)-325.0/54.0*(Nu+Nd)*(Nu+Nd)-2857.0/2.0)*pow(g3,7)/pow(4.0*pi,6);
}
}
}
void EWCouplings::betaYukawa(const double x, vector< Complex > &y, vector<Complex > &dydx) {
using Constants::pi;
const Complex I(0,1.0);
boost::numeric::ublas::identity_matrix<Complex> Id(3,3);
// Yukawa
boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
Yuk_u(ix,iy) = y[21+3*ix+iy];
Yuk_d(ix,iy) = y[12+3*ix+iy];
Yuk_e(ix,iy) = y[ 3+3*ix+iy];
}
}
// gauge
double Ng = 0.5*numberOfFlavours(x);
boost::numeric::ublas::vector<Complex> gauge(3);
for(unsigned int l=0; l<3; l++) gauge[l] = y[l];
Complex lambda = y[30];
// traces of yukawa matrices
boost::numeric::ublas::matrix<Complex> mTemp(3,3),MUU(3,3),MDD(3,3),MLL(3,3),
MUU2(3,3),MDD2(3,3),MLL2(3,3),MUUDD(3,3),MDDUU(3,3);
axpy_prod(herm(Yuk_u),Yuk_u,MUU);
Complex trU = trace( MUU);
axpy_prod(MUU,MUU,MUU2);
Complex trUU = trace(MUU2);
axpy_prod(herm(Yuk_d),Yuk_d,MDD);
Complex trD = trace( MDD);
axpy_prod(MDD,MDD,MDD2);
Complex trDD = trace(MDD2);
axpy_prod(MUU,MDD,MUUDD);
Complex trUD = trace(MUUDD);
axpy_prod(MDD,MUU,MDDUU);
axpy_prod(herm(Yuk_e),Yuk_e,MLL);
Complex trL = trace( MLL);
axpy_prod(MLL,MLL,MLL2);
Complex trLL = trace(MLL2);
Complex g02 = sqr(gauge[0]);
Complex g12 = sqr(gauge[1]);
Complex g22 = sqr(gauge[2]);
// Evaluate beta functions for yukawa couplings
boost::numeric::ublas::zero_matrix<Complex> zero3x3(3,3);
boost::numeric::ublas::matrix<Complex> dYuk_udx(zero3x3), dYuk_ddx(zero3x3), dYuk_edx(zero3x3), beta1_u(zero3x3),
beta1_d(zero3x3), beta1_e(zero3x3), beta2_u(zero3x3), beta2_d(zero3x3), beta2_e(zero3x3);
Complex Y2 = 3.0*trU+3.0*trD + trL;
Complex Y4 = (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trU +
(1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trD + 3.0/4.0*(g02 + g12)*trL;
Complex chi4 = 27.0/4.0*trUU + 27.0/4.0*trDD + 9.0/4.0*trLL - 6.0/4.0*trUD;
if(loops_ >= 1) {
beta1_u = 3.0/2.0*(MUU - MDD) + (Y2 - 17.0/20.0*g02 - 9.0/4.0*g12 - 8.0*g22)*Id;
beta1_d = 3.0/2.0*(MDD - MUU) + (Y2 - 1.0/4.0*g02 - 9.0/4.0*g12 - 8.0*g22)*Id;
beta1_e = 3.0/2.0*(MLL) + (Y2 - 9.0/4.0*g02 - 9.0/4.0*g12)*Id;
axpy_prod(Yuk_u,beta1_u,mTemp);
dYuk_udx += (1.0/(16.0*pow(pi,2)))*mTemp;
axpy_prod(Yuk_d,beta1_d,mTemp);
dYuk_ddx += (1.0/(16.0*pow(pi,2)))*mTemp;
axpy_prod(Yuk_e,beta1_e,mTemp);
dYuk_edx += (1.0/(16.0*pow(pi,2)))*mTemp;
if (loops_ >= 2) {
Complex l2=sqr(lambda);
beta2_u = 3.0/2.0*MUU2 - MUUDD - 1.0/4.0*MDDUU + 11.0/4.0*MDD2 + Y2*(5.0/4.0*MDD - 9.0/4.0*MUU) +
(-chi4 + 3.0/2.0*l2)*Id - 2.0*lambda*(3.0*MUU + MDD) +
(223.0/80.0*g02 + 135.0/16.0*g12 + 16.0*g22)*(MUU) -
(43.0/80.0*g02 - 9.0/16.0*g12 + 16.0*g22)*(MDD) + 5.0/2.0*Y4*Id +
((9.0/200.0 + 29.0/45.0*Ng)*pow(gauge[0],4) - 9.0/20.0*pow(gauge[0]*gauge[1],2) +
19.0/15.0*pow(gauge[0]*gauge[2],2) - (35.0/4.0 - Ng)*pow(gauge[1],4) + 9.0*pow(gauge[1]*gauge[2],2) -
(404.0/3.0 - 80.0/9.0*Ng)*pow(gauge[2],4))*Id;
beta2_d = 3.0/2.0*MDD2 - MDDUU - 1.0/4.0*MUUDD + 11.0/4.0*MUU2 + Y2*(5.0/4.0*MUU - 9.0/4.0*MDD) + (-chi4 + 3.0/2.0*l2)*Id - 2.0*lambda*(3.0*MDD + MUU) + (187.0/80.0*g02 + 135.0/16.0*g12 + 16.0*g22)*(MDD) - (79.0/80.0*g02 - 9.0/16.0*g12 + 16.0*g22)*(MUU) + 5.0/2.0*Y4*Id - ((29.0/200.0 + 1.0/45.0*Ng)*pow(gauge[0],4) - 27.0/20.0*pow(gauge[0]*gauge[1],2) + 31.0/15.0*pow(gauge[0]*gauge[2],2) - (35.0/4.0 - Ng)*pow(gauge[1],4) + 9.0*pow(gauge[1]*gauge[2],2) - (404.0/3.0 - 80.0/9.0*Ng)*pow(gauge[2],4))*Id;
beta2_e = 3.0/2.0*MLL2 - 9.0/4.0*Y2*MLL + (-chi4 + 3.0/2.0*l2)*Id - 6.0*lambda*(MLL) + (387.0/80.0*g02 + 135.0/15.0*g12)*(MLL) + 5.0/2.0*Y4*Id + ((51.0/200.0 + 11.0/5.0*Ng)*pow(gauge[0],4) + 27.0/20.0*pow(gauge[0]*gauge[1],2) - (35.0/4.0 - Ng)*pow(gauge[1],4))*Id;
axpy_prod(Yuk_u,beta2_u,mTemp);
dYuk_udx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
axpy_prod(Yuk_d,beta2_d,mTemp);
dYuk_ddx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
axpy_prod(Yuk_e,beta2_e,mTemp);
dYuk_edx += (1.0/pow(16.0*pow(pi,2),2))*mTemp;
}
}
boost::numeric::ublas::vector<Complex> temp(27);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
dydx[ 3+ix+3*iy] = dYuk_edx(ix,iy);
dydx[12+ix+3*iy] = dYuk_ddx(ix,iy);
dydx[21+ix+3*iy] = dYuk_udx(ix,iy);
}
}
}
void EWCouplings::betaHiggs(const double x, vector<Complex> &y,
vector<Complex> &dydx) {
using Constants::pi;
const Complex I(0,1.0);
double Ng = 0.5*numberOfFlavours(x);
// Yukawa
boost::numeric::ublas::matrix<Complex> Yuk_u(3,3), Yuk_d(3,3), Yuk_e(3,3);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
Yuk_u(ix,iy) = y[21+3*ix+iy];
Yuk_d(ix,iy) = y[12+3*ix+iy];
Yuk_e(ix,iy) = y[ 3+3*ix+iy];
}
}
// gauge
boost::numeric::ublas::vector<Complex> gauge(3);
for(int l=0; l<3; l++) gauge(l) = y[l];
Complex lambda = y[30];
complex<Energy> vev = y[31]*GeV;
// Evaluate beta functions for higgs coupling
Complex beta1_lambda(0.), beta2_lambda(0.), gamma1_vev(0.), gamma2_vev(0.),
Y2(0.), H(0.), Y4(0.), chi4(0.);
// traces of yukawa matrices
boost::numeric::ublas::matrix<Complex> temp(3,3),temp2(3,3),MUU(3,3),MDD(3,3),MLL(3,3);
axpy_prod(herm(Yuk_u),Yuk_u,MUU);
Complex trU = trace( MUU);
axpy_prod(MUU,MUU,temp);
Complex trUU = trace(temp);
axpy_prod(MUU,temp,temp2);
Complex trUUU = trace(temp2);
axpy_prod(herm(Yuk_d),Yuk_d,MDD);
Complex trD = trace( MDD);
axpy_prod(MDD,MDD,temp);
Complex trDD = trace(temp);
axpy_prod(MDD,temp,temp2);
Complex trDDD = trace(temp2);
axpy_prod(MUU,MDD,temp);
Complex trUD = trace(temp);
axpy_prod(herm(Yuk_e),Yuk_e,MLL);
Complex trL = trace( MLL);
axpy_prod(MLL,MLL,temp);
Complex trLL = trace(temp);
axpy_prod(MLL,temp,temp2);
Complex trLLL = trace(temp2);
axpy_prod(MUU+MDD,MDD,temp);
axpy_prod(MUU,temp,temp2);
Complex trUUDD = trace(temp2);
Complex g02 = sqr(gauge[0]);
Complex g12 = sqr(gauge[1]);
Complex g22 = sqr(gauge[2]);
Complex g04 = sqr(g02);
Complex g14 = sqr(g12);
double pi2 = sqr(pi);
Y2 = 3.0*trU+3.0*trD + trL;
Y4 = (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*(trU) + (1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*(trD) + 3.0/4.0*(g02 + g12)*(trL);
chi4 = 27.0/4.0*trUU + 27.0/4.0*trDD + 9.0/4.0*trLL - 6.0/4.0*trUD;
H = 3.0*trUU + 3.0*trDD + trLL;
if(loops_ >= 1) {
Complex l2=sqr(lambda);
beta1_lambda = 12.0*l2 - (9.0/5.0*g02 + 9.0*g12)*lambda + 9.0/4.0*(3.0/25.0*g04+2.0/5.0*g02*g12 + g14) + 4.0*Y2*lambda - 4.0*H;
gamma1_vev = 9.0/4.0*(1.0/5.0*g02+g12)-Y2;
dydx[30] += 1.0/(16.0*pi2)*beta1_lambda;
dydx[31] += vev/(16.0*pi2)*gamma1_vev/GeV;
if (loops_ >= 2) {
beta2_lambda = -78.0*lambda*l2 + 18.0*(3.0/5.0*g02 + 3.0*g12)*l2 -
( (313.0/8.0 - 10.0*Ng)*g14 - 117.0/20.0*g02*g12 + 9.0/25.0*(229.0/4.0+50.0/9.0*Ng)*g04 )*lambda +
(497.0/8.0 - 8.0*Ng)*g14*g12 - 3.0/5.0*(97.0/24.0 + 8.0/3.0*Ng)*g02*g14 -
9.0/25.0*(239.0/24.0 + 40.0/9.0*Ng)*g04*g12 - 27.0/125.0*(59.0/24.0 + 40.0/9.0*Ng)*g04*g02 -
64.0*g22*(trUU + trDD) - 8.0/5.0*g02*(2.0*trUU - trDD + 3.0*trLL) - 3.0/2.0*g14*Y4 +
10.0*lambda*( (17.0/20.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trU + (1.0/4.0*g02 + 9.0/4.0*g12 + 8.0*g22)*trD + 3.0/4.0*(g02 + g12)*trL ) +
3.0/5.0*g02*( (-57.0/10.0*g02 + 21.0*g12 )*trU + (3.0/2.0*g02 + 9.0*g12)*trD + (-15.0/2.0*g02 + 11.0*g12)*trL ) - 24.0*l2*Y2 - lambda*H +
6.0*lambda*trUD + 20.0*(3.0*trUUU + 3.0*trDDD + trLLL) - 12.0*trUUDD;
gamma2_vev = -3.0/2.0*l2 - 5.0/2.0*Y4 + chi4 - 27.0/80.0*g02*g12 -
(93.0/800.0 + 1.0/2.0*Ng)*g04 + (511.0/32.0 - 5.0/2.0*Ng)*g14;
dydx[30] += 1.0/pow(16.0*pi2,2)*beta2_lambda;
dydx[31] += vev/pow(16.0*pi2,2)*gamma2_vev/GeV;
}
}
}
double EWCouplings::interpolate(double t, int paramIndex) {
double stepsize = table_(1,0)-table_(0,0);
double tol = 0.001*stepsize;
double logewScale = log(ewScale_/GeV);
double loghighScale = log(highScale_/GeV);
if (t<logewScale-tol || t>loghighScale+tol) {
cerr << "Stepsize: " << stepsize << std::endl;
cerr << "tol: " << tol << std::endl;
cerr << "logewScale: " << logewScale << std::endl;
cerr << "loghighScale: " << loghighScale << std::endl;
cerr << "paramIndex: " << paramIndex << std::endl;
cerr << "t: " << t << std::endl;
cerr << "Couplings::_Interp(double t, int parmIndex) trying to obtain parameter ";
cerr << "value outside the available range. Returning zero." << std::endl;
assert(false);
}
// return value at EW scale
if (abs(t-logewScale)<tol) return table_(0,paramIndex);
// return value at high scale
if (abs(t-loghighScale)<tol) {
return table_(table_.size1()-1,paramIndex);
}
unsigned int numSteps = int((t-table_(0,0))/stepsize);
// Linear Interpolation:
double x1 = table_(numSteps,0);
double y1 = table_(numSteps,paramIndex);
double x2 = table_(numSteps+1,0);
double y2 = table_(numSteps+1,paramIndex);
return y1+((y2-y1)/(x2-x1))*(t-x1);
}
double EWCouplings::interpolateLow(double t, int paramIndex) {
double stepsize = lowTable_(0,0)-lowTable_(1,0);
double tol = 0.00001*stepsize;
double logewScale = log(ewScale_ /GeV);
double loglowScale = log(lowScale_/GeV);
if (t<loglowScale-tol || t>logewScale+tol) {
cerr<< "Couplings::_LowInterp(double t, int parmIndex) trying to obtain parameter ";
cerr << "value outside the available range. Returning zero." << std::endl;
assert(false);
}
if (abs(t-logewScale)<tol) {
return lowTable_(0,paramIndex);
}
if (std::abs(t-loglowScale)<tol) {
return lowTable_(lowTable_.size1()-1,paramIndex);
}
int numSteps = (int)((lowTable_(0,0)-t)/stepsize);
// Linear Interpolation:
double x1 = lowTable_(numSteps,0);
double y1 = lowTable_(numSteps,paramIndex);
double x2 = lowTable_(numSteps+1,0);
double y2 = lowTable_(numSteps+1,paramIndex);
return y1+((y2-y1)/(x2-x1))*(t-x1);
}
diff --git a/MatrixElement/EW/EWCouplings.h b/MatrixElement/EW/EWCouplings.h
--- a/MatrixElement/EW/EWCouplings.h
+++ b/MatrixElement/EW/EWCouplings.h
@@ -1,449 +1,473 @@
// -*- C++ -*-
#ifndef Herwig_EWCouplings_H
#define Herwig_EWCouplings_H
//
// This is the declaration of the EWCouplings class.
//
#include "ThePEG/Interface/Interfaced.h"
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include "EWCouplings.fh"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the EWCouplings class.
*
* @see \ref EWCouplingsInterfaces "The interfaces"
* defined for EWCouplings.
*/
class EWCouplings: public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
- EWCouplings(unsigned int loops=2, unsigned int steps=200,
- Energy highScale=14.*TeV, Energy lowScale=10.*GeV);
+ EWCouplings(unsigned int loops=2, unsigned int steps=200, Energy lowScale=10.*GeV);
/**
* The destructor.
*/
virtual ~EWCouplings();
//@}
/**
* Initialize the running couplings
*/
void initialize();
/**
* Number of dynamical quarks at $\log\mu = x$ (in GeV)
* N.B.Integrate out top quark at Mz, not at Mt.
*/
unsigned int numberOfFlavours(double x) {
return x >= std::log(ewScale_/GeV) ? 6 : 5;
}
public:
/**
* Set whether or not to include \f$SU(3)\f$ running
*/
void SU3(bool includeSU3) {includeSU3_ = includeSU3;}
/**
* Whether or not to include \f$SU(3)\f$ running
*/
bool SU3() { return includeSU3_;}
/**
* Set whether or not to include EW running
*/
void EW(bool includeEW) {includeEW_ = includeEW;}
/**
* Whether or not to include EW running
*/
bool EW() { return includeEW_;}
/**
* alpha for the U1 gauge coupling at energy mu (in GeV):
*/
double a1(Energy mu) {
if (includeEW_) {
if (mu>=ewScale_) {
return (3.0/5.0)*interpolate(log(mu/GeV),1);
}
return interpolateLow(log(mu/GeV),1);
}
else
return 0.0;
}
/**
* alpha for the SU2 gauge coupling at energy mu (in GeV):
*/
double a2(Energy mu) {
if (includeEW_) {
if (mu<ewScale_) {
return 0.0;
}
else
return interpolate(log(mu/GeV),2);
}
else
return 0.0;
}
/**
* alpha for the SU3 gauge coupling at energy mu (in GeV):
*/
double a3(Energy mu) {
if(includeSU3_) {
if (mu>=ewScale_) {
return interpolate(log(mu/GeV),3);
}
else {
return interpolateLow(log(mu/GeV),2);
}
}
else
return 0.0;
}
/**
* alpha for EM
*/
double aEM(Energy mu) {
if (includeEW_) {
if (mu<=ewScale_) {
return interpolateLow(log(mu/GeV),1);
}
else {
double alpha1=a1(mu);
double alpha2=a2(mu);
return alpha1*alpha2/(alpha1+alpha2);
}
}
return 0.0;
}
double aS(Energy mu) {
if(includeSU3_) {
if (mu<=ewScale_) {
return interpolateLow(log(mu/GeV),2);
}
else {
return interpolate(log(mu/GeV),3);
}
}
else return 0.0;
}
/**
* Top quark Yukawa coupling
*/
double y_t(Energy mu) {
if (includeEW_) {
if(mu<ewScale_)
return 0.0;
else
return interpolate(log(mu/GeV),4);
}
else
return 0.0;
}
/**
* Quartic scalar coupling lambda (normalization different, hence factor of 2):
*/
double lambda(Energy mu) {return 2.0*interpolate(log(mu/GeV),5);}
/**
* VEV
*/
Energy vev(Energy mu) {return interpolate(log(mu/GeV),6)*GeV;}
/**
* \f$\lambda_t\f$
*/
double lambda_t(Energy mu) {return y_t(mu);}
/**
* Z couplings
*/
//@{
double g_Lu(Energy mu) {return 0.5-(2.0/3.0)*Sin2thW(mu);}
double g_Ld(Energy mu) {return -0.5+(1.0/3.0)*Sin2thW(mu);}
double g_Le(Energy mu) {return (-1.0/2.0)+Sin2thW(mu);}
double g_Lnu(Energy ) {return (0.5);}
double g_Ru(Energy mu) {return (-2.0/3.0)*Sin2thW(mu);}
double g_Rd(Energy mu) {return (1.0/3.0)*Sin2thW(mu);}
double g_Re(Energy mu) {return Sin2thW(mu);}
double g_W(Energy mu) {return Cos2thW(mu);}
double g_phiPlus(Energy mu) {return 0.5-Sin2thW(mu);}
//@}
/**
* \f\cos^2\theta_W\f$
*/
double Cos2thW(Energy ) {
//\todo why this value?
//return mW()*mW()/(mZ()*mZ());
return (1.-0.2314);
}
/**
* \f\sin^2\theta_W\f$
*/
double Sin2thW(Energy mu) {return 1.-Cos2thW(mu);}
double aW(Energy mu) {return a2(mu);}
double aZ(Energy mu) {
//return a2(mu)/Cos2thW(mu); // Same thing, actually
return a1(mu)+a2(mu);
}
public:
/**
* Masses of the Standard Model particles
*/
//@{
/**
* Z mass
*/
Energy mZ() const { return mZ_;}
/**
* W Mass
*/
Energy mW() const { return mW_;}
/**
* Top quark mass
*/
Energy mT() const {return mT_;}
Energy mTatmZ() { return mT_;}
/**
* Higg boson mass
*/
Energy mH() const {return mH_;}
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* Set the initial value of the couplings
*/
void initializeCouplings(vector<Complex> & y);
/**
* Assigns numerical values to beta functions
* Takes in a point x = log(mu) and the values of y[i] at x and assigns dydx[i] the
* value beta[i](mu). The function Derivs farms out the plugging in to the three
* functions BetaGauge, BetaYukawa, and BetaHiggs, which evaluates the beta functions
* for the gauge couplings, yukawa matrices, and higgs quartic coupling/vev, respectively.
*/
void derivatives(const double x, vector< Complex > &y,
vector< Complex > & dydx);
/**
* Beta function for the gauge interactions
*/
void betaGauge(const double x, vector<Complex> &y, vector<Complex> &dydx);
/**
* Beta function for the gauge interactions at low scales
*/
void lowBetaGauge(const double x, vector<Complex> &y, vector<Complex> &dydx);
/**
* Beta function for the Yukawa interactions
*/
void betaYukawa(const double x, vector<Complex> &y, vector<Complex> &dydx);
/**
* Beta function for the Higgs interactions
*/
void betaHiggs(const double x, vector<Complex> &y, vector<Complex> &dydx);
/**
* Update the couplings using 4-th order RK
* Takes in a vector y[i] of function values at a point x and the values of the
* first derivatives dydx[i] ( = dy[i]/dx ) alon with a step size stepsize. The
* function then defines assigns the value of y[i](x+stepsize) to the variable yout[i].
* (Adapted from sample code in Numerical Recipes in C++ Press, Teukolsky, et. al.)
*/
void RK4(vector<Complex> & y, vector<Complex> &dydx, const double x,
const double stepsize, vector<Complex> &yout);
/**
* Initialize the low energy parameters
*/
void initializeLow();
/**
* Interpolate the table, t = ln(mu)
*/
double interpolate(double t, int paramIndex);
/**
* Interpolate the tabel, t = ln(mu)
*/
double interpolateLow(double t, int paramIndex);
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
EWCouplings & operator=(const EWCouplings &);
private:
/**
* Electoweak Scale
*/
Energy ewScale_;
/**
* High Scale
*/
Energy highScale_;
/**
* Low Scale
*/
Energy lowScale_;
/**
* Whether or not to include SU3
*/
bool includeSU3_;
/**
* Whether or not to include EW
*/
bool includeEW_;
/**
* Whether or not the running couplings have been initialized
*/
bool initialized_;
/**
* Masses of Standard Model particles
*/
//@{
/**
* Mass Choice
*/
bool massChoice_;
/**
* Z mass
*/
Energy mZ_;
/**
* W mass
*/
Energy mW_;
/**
* Top mass
*/
Energy mT_;
/**
* Higgs boson mass
*/
Energy mH_;
//@}
/**
* Number of loops
*/
unsigned int loops_;
/**
* Number of steps for Runga-Kutta (High scale)
*/
unsigned int highSteps_;
/**
* Number of steps for Runga-Kutta (Low scale)
*/
unsigned int lowSteps_;
/**
* Matrix to store the parameters
*/
boost::numeric::ublas::matrix<double> table_;
/**
* Matrix to store the low energy parameters.
* This will hold only aEM and a3 at mu<ewScale
*/
boost::numeric::ublas::matrix<double> lowTable_;
+ /**
+ * Input values of the couplings
+ */
+ //@{
+ /**
+ * \f%\alpha_1(M_Z)\f$
+ */
+ double a1MZ_;
+
+ /**
+ * \f%\alpha_2(M_Z)\f$
+ */
+ double a2MZ_;
+
+ /**
+ * \f%\alpha_S(M_Z)\f$
+ */
+ double aSMZ_;
+
+ /**
+ * \f$\lambda_t\f$
+ */
+ double lambdat_;
+ //@}
+
};
}
#endif /* Herwig_EWCouplings_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 4:39 PM (21 h, 49 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023483
Default Alt Text
(36 KB)

Event Timeline