Page MenuHomeHEPForge
Files F8309915

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
--- a/MatrixElement/EW/GroupInvariants.h
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -1,328 +1,364 @@
// -*- C++ -*-
//
// GroupInvariants.h is a part of Herwig - A multi-purpose Monte Carlo event generator
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
//
#ifndef HERWIG_GroupInvariants_H
#define HERWIG_GroupInvariants_H
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Config/Unitsystem.h"
#include <cassert>
#include <boost/numeric/ublas/matrix.hpp>
namespace Herwig {
using namespace ThePEG;
namespace GroupInvariants {
/**
* Simple struct for storing the different gauge contributions
*/
struct GaugeContributions {
/**
* Default Constructor
*/
GaugeContributions(double inSU3=0.,
double inSU2=0., double inU1=0.)
: SU3(inSU3),SU2(inSU2),U1(inU1)
{}
/**
* \f$SU(3)\f$
*/
double SU3;
/**
* \f$SU(2)\f$
*/
double SU2;
/**
* \f$U(1)\f$
*/
double U1;
};
/**
* The \f$SU(N)\f$ \f$C_A\f$
*/
inline double C_A(unsigned int N) {
return N !=1 ? double(N) : 0.;
}
/**
* The \f$SU(N)\f$ \f$C_F\f$
*/
inline double C_F(unsigned int N) {
return N !=1 ? 0.5*(double(N*N)-1.)/double(N) : 1.;
}
/*
* The \f$SU(N)\f$ \f$C_d\f$
*/
inline double C_d(unsigned int N) {
return (double(N*N)-4.)/double(N);
}
/**
* The \f$SU(N)\f$\f$C_1\f$
*/
inline double C_1(unsigned int N) {
double N2(N*N);
return 0.25*(N2-1.0)/N2;
}
/**
* \f$T_F\f$
*/
inline double T_F(unsigned int N, bool high) {
if(high) {
return N !=1 ? 0.5 : 5.0/3.0;
}
else {
return N !=1 ? 0.5 : 20.0/3.0;
}
}
/**
* \f$t_S\f$
*/
inline double t_S(unsigned int, bool ) {
return 0.5;
}
/**
* / Number of complex scalars in the fundamental rep. of SU(N)/U(1)
*/
inline double n_S(unsigned int N, bool high) {
if(high) {
if(N==2 || N==1) return 1.0;
else if(N==3) return 0.0;
else assert(false);
}
else {
if(N>=1&&N<=3) return 0.;
else assert(false);
}
}
/**
* Number of Dirac Fermions in the fund. rep. of SU(N) (or U(1) for N==1)
*/
inline double n_F(unsigned int N, bool high) {
if(high) {
if(N==1) return 3.0;
else if(N==2 || N==3) return 6.0;
else assert(false);
}
else {
if(N==1) return 1.0;
else if(N==2) return 0.0;
else if(N==3) return 5.0;
else assert(false);
}
}
/**
* Find K_i for gauge group N. high=false for running at mu<mZ
*/
double K_Factor(unsigned int i,unsigned int N, bool high);
/**
* Find B_i for gauge group N, high energy
*/
double B_Factor(int i, int N, bool fermion, bool longitudinal);
/**
* Find B_i for gauge group N, low energy
*/
double B_Factor_Low(int i, int N, bool fermion, double boostFactor);
/**
* Contributions to the Cusps
*/
GaugeContributions cuspContributions(Energy mu, int K_ORDER, bool high);
/**
* Contributions to B, high energy
*/
GaugeContributions BContributions(Energy mu, int B_ORDER,
bool fermion, bool longitudinal);
/**
* Contributions to B, low energy
*/
GaugeContributions BContributionsLow(Energy mu, int B_ORDER,
bool fermion, double boostFactor);
inline Complex PlusLog(double arg) {
static const Complex I(0,1.0);
if (arg>0.0)
return log(arg);
else if (arg<0.0)
return log(-arg)+I*Constants::pi;
else
assert(false);
}
inline Complex MinusLog(double arg) {
static const Complex I(0,1.0);
if (arg>0.0)
return log(arg);
else if (arg<0.0)
return log(-arg)-I*Constants::pi;
else
assert(false);
}
inline Complex getT(Energy2 s, Energy2 t) {
return MinusLog(-t/GeV2) - MinusLog(-s/GeV2);
}
inline Complex getU(Energy2 s, Energy2 u) {
return MinusLog(-u/GeV2) - MinusLog(-s/GeV2);
}
inline boost::numeric::ublas::matrix<Complex> Gamma2(Complex U, Complex T) {
boost::numeric::ublas::matrix<Complex> output(2,2);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) = (-3.0/2.0)*I*pi + (T+U);
output(1,1) = (-3.0/2.0)*I*pi;
output(0,1) = 2.0*(T-U);
output(1,0) = (3.0/8.0)*(T-U);
return output;
}
inline boost::numeric::ublas::matrix<Complex> Gamma2w(Complex U, Complex T) {
boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(5,5);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) += -I*pi*11.0/4.0;
output(0,1) += U-T;
output(1,0) += 2.0*(U-T);
output(1,1) += -I*pi*11.0/4.0 + (T+U);
output(2,2) += -7.0/4.0*I*pi + (U+T);
output(3,3) += -7.0/4.0*I*pi + (U+T);
output(4,4) += -3.0/4.0*I*pi;
return output;
}
inline boost::numeric::ublas::matrix<Complex> Gamma2Singlet() {
boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(2,2);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) = output(1,1) = -3.0/4.0*I*pi;
return output;
}
inline Complex Gamma1(double hypercharge) {
Complex I(0,1.0);
return -I*Constants::pi*(hypercharge*hypercharge);
}
inline Complex Gamma1(double y1, double y2, Complex T, Complex U) {
Complex I(0,1.0);
return -I*Constants::pi*(y1*y1+y2*y2) + 2.0*y1*y2*(T-U);
}
inline Complex Gamma1(double y1, double y2, double y3, double y4,
Complex T, Complex U) {
Complex I(0,1.0);
return -I*Constants::pi*(y1*y1+y2*y2+y3*y3+y4*y4)/2.0 +
(y1*y4+y2*y3)*T - (y1*y3+y2*y4)*U;
}
+ inline boost::numeric::ublas::matrix<Complex> Gamma3(Complex U, Complex T) {
+ boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(2,2);
+ static const Complex I(0,1.0);
+ using Constants::pi;
+ output(0,0) = -(8.0/3.0)*I*pi;
+ output(1,1) = -(8.0/3.0)*I*pi;
+ output(0,0) += (7.0/3.0)*T + (2.0/3.0)*U;
+ output(0,1) = 2.0*(T-U);
+ output(1,0) = (4.0/9.0)*(T-U);
+ return output;
+ }
+
+ inline boost::numeric::ublas::matrix<Complex> Gamma3g(Complex U, Complex T) {
+ boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(3,3);
+ static const Complex I(0,1.0);
+ using Constants::pi;
+ output(0,2) = U-T;
+ output(1,1) = 3.0/2.0*(T+U);
+ output(1,2) = 3.0/2.0*(U-T);
+ output(2,0) = 2.0*(U-T);
+ output(2,1) = 5.0/6.0*(U-T);
+ output(2,2) = 3.0/2.0*(T+U);
+ for (unsigned int i=0; i<3; i++) {
+ output(i,i) += -13.0/3.0*I*pi;
+ }
+ return output;
+ }
+
+ inline boost::numeric::ublas::matrix<Complex> Gamma3Singlet() {
+ boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(2,2);
+ static const Complex I(0,1.0);
+ using Constants::pi;
+ output(0,0) = output(1,1) = -4.0/3.0*I*pi;
+ return output;
+ }
+
/**
* Number of fermion generations (only used in gauge boson HighCMatching)
*/
inline double n_g() { return 3.0; }
/**
* Number of complex scalars in the fundamental rep. of SU(N)
*/
inline double nSWeyl(unsigned int N, bool high) {
if(high) {
if(N==2 || N==1) return 1.0;
else if (N==3) return 0.0;
else assert(false);
}
else {
if( N==1 || N==3 ) return 0.0;
else assert(false);
}
}
/**
* Number of Weyl Fermions in the fundamental rep. of SU(N)
*/
inline double nFWeyl(unsigned int N, bool high) {
if(high) {
if(N==2 || N==3) return 12.0;
else assert(false);
}
else {
if(N==3) return 10.0;
else if(N==1) return 2.0;
else assert(false);
}
}
inline double TFWeyl(unsigned int) {
return 0.5;
}
inline double tSWeyl(unsigned int) {
return 0.5;
}
inline Complex WFunction(Energy mu, Energy2 s) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ln = MinusLog(-s/(mu*mu));
return (-1.0*ln*ln + 3.0*ln+pi*pi/6.0-8.0);
}
/**
* \fX_N\f% function, v is either t or u
*/
inline Complex XNFunction(unsigned int N, Energy mu, Energy2 s, Energy2 v) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ls = MinusLog(-s/(mu*mu));
return (2.0*C_F(N)*WFunction(mu,s) +
C_A(N)*(2.0*ls*ls - 2.0*MinusLog((s+v)/(mu*mu))*ls -
11.0/3.0*ls + pi*pi + 85.0/9.0) +
(2.0/3.0*ls - 10.0/9.0) * TFWeyl(N) * nFWeyl(N,true) +
(1.0/3.0*ls - 8.0/9.0) * TFWeyl(N) * nSWeyl(N,true));
}
/**
* \f$\Pi_1\f$ function
*/
inline Complex PI1_function(Energy mu, Energy2 s) {
assert(abs(s)>ZERO);
return ((41.0/6.0)*MinusLog(-s/(mu*mu))-104.0/9.0);
}
/**
* \f$\tilde{f}\f$ function, v is either t or u
*/
inline Complex fTildeFunction(Energy mu, Energy2 s, Energy2 v) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ls = MinusLog(-s/GeV2), lv = MinusLog(-v/GeV2);
Complex lsv = MinusLog((s+v)/GeV2);
return (-2.0*double(s/(s+v))*(lv-ls) +
double(s*(s+2.0*v)/((s+v)*(s+v))) * ((lv-ls)*(lv-ls) + pi*pi) +
4.0*MinusLog(-s/(mu*mu))*(lv-lsv));
}
}
}
#endif // HERWIG_GroupInvariants_H
diff --git a/MatrixElement/EW/SoftSudakov.cc b/MatrixElement/EW/SoftSudakov.cc
--- a/MatrixElement/EW/SoftSudakov.cc
+++ b/MatrixElement/EW/SoftSudakov.cc
@@ -1,92 +1,881 @@
// -*- 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"
using namespace Herwig;
using namespace GroupInvariants;
SoftSudakov::SoftSudakov() : K_ORDER_(3) {}
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));
}
}
// Resummed:
//return gamma.exp();
// Fixed Order:
return boost::numeric::ublas::identity_matrix<Complex>(NN,NN) + gamma;
}
+
+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);
+ }
+
+ 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);
+ 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/SoftSudakov.h b/MatrixElement/EW/SoftSudakov.h
--- a/MatrixElement/EW/SoftSudakov.h
+++ b/MatrixElement/EW/SoftSudakov.h
@@ -1,162 +1,181 @@
// -*- C++ -*-
#ifndef Herwig_SoftSudakov_H
#define Herwig_SoftSudakov_H
//
// This is the declaration of the SoftSudakov class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/Utilities/GSLIntegrator.h"
#include <boost/numeric/ublas/matrix.hpp>
+#include "EWProcess.h"
#include "SoftSudakov.fh"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the SoftSudakov class.
*
* @see \ref SoftSudakovInterfaces "The interfaces"
* defined for SoftSudakov.
*/
class SoftSudakov: public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SoftSudakov();
/**
* The destructor.
*/
virtual ~SoftSudakov();
//@}
public:
+
+ /**
+ * Low energy soft evolution
+ */
+ boost::numeric::ublas::matrix<Complex>
+ lowEnergyRunning(Energy EWScale, Energy lowScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ Herwig::EWProcess::Process process);
+
+ /**
+ * Evalaute the high energy running as a matrix
+ */
+ boost::numeric::ublas::matrix<Complex>
+ highEnergyRunning(Energy highScale, Energy EWScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ Herwig::EWProcess::Process process);
+
+protected:
/**
- * Evaluate bthe soft evolution
+ * Evaluate the soft evolution
*/
boost::numeric::ublas::matrix<Complex> 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);
public:
/**
* The operator to be integrated
*/
InvEnergy operator ()(Energy mu) const;
/** Argument type for GaussianIntegrator */
typedef Energy ArgType;
/** Return type for GaussianIntegrator */
typedef InvEnergy ValType;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @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.
*/
SoftSudakov & operator=(const SoftSudakov &);
private:
/**
* Order for K
*/
unsigned int K_ORDER_;
/**
* Integrator
*/
GSLIntegrator integrator_;
/**
* Whether doing the high or low scale contribution
*/
bool high_;
/**
* Column
*/
unsigned int row_;
/**
* Row
*/
unsigned int col_;
/**
*
*/
bool real_;
/**
*
*/
boost::numeric::ublas::matrix<Complex> G1_;
/**
*
*/
boost::numeric::ublas::matrix<Complex> G2_;
/**
*
*/
boost::numeric::ublas::matrix<Complex> G3_;
};
}
#endif /* Herwig_SoftSudakov_H */

File Metadata

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

Event Timeline