Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8310006
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
63 KB
Subscribers
None
View Options
diff --git a/MatrixElement/EW/EWProcess.h b/MatrixElement/EW/EWProcess.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/EWProcess.h
@@ -0,0 +1,69 @@
+// -*- C++ -*-
+//
+// EWProcess.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_EWProcess_H
+#define HERWIG_EWProcess_H
+
+namespace Herwig {
+
+namespace EWProcess {
+
+ /**
+ * Enumerates the processes for which we have SCET Wilson Coefficients
+ */
+ enum Process { QQQQ,
+ QQQQiden,
+ QtQtQQ,
+ QQUU,
+ QtQtUU,
+ QQtRtR,
+ QQDD,
+ QtQtDD,
+ QQLL,
+ QQEE,
+ UUUU,
+ UUUUiden,
+ tRtRUU,
+ UUDD,
+ tRtRDD,
+ UULL,
+ UUEE,
+ DDDD,
+ DDDDiden,
+ DDLL,
+ DDEE,
+ LLLL,
+ LLLLiden,
+ LLEE,
+ EEEE,
+ EEEEiden,
+ QQWW,
+ QQPhiPhi,
+ QQWG,
+ QQBG,
+ QQGG,
+ QtQtGG,
+ LLWW,
+ LLPhiPhi,
+ UUBB,
+ UUPhiPhi,
+ UUBG,
+ UUGG,
+ tRtRGG,
+ DDBB,
+ DDPhiPhi,
+ DDBG,
+ DDGG,
+ EEBB,
+ EEPhiPhi,
+ LAST};
+}
+}
+
+#endif // HERWIG_EWProcess_H
diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
--- a/MatrixElement/EW/GroupInvariants.h
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -1,254 +1,269 @@
// -*- 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>
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);
}
+
+ /**
+ * 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));
+ }
}
-
-
-
-// namespace DiracHigh {
-
-
-// double n_g() { // Number of fermion generations (only used in gauge boson HighCMatching)
-// return 3.0;
-// }
-
-
-// }
-
-
-// namespace WeylHigh {
-
-// double n_S(int N) { // Number of complex scalars in the fundamental rep. of SU(N)/U(1)
-// if(N==2 || N==1) return 1.0;
-// if(N==3) return 0.0;
-// std::cout << "Error! SU(N), N != (1, 2 or 3) used for n_S in ";
-// std::cout << "GroupInvariants.h but not defined." << std::endl;
-// return 0.0;
-// }
-// double n_F(int N) { // Number of Weyl Fermions in the fundamental rep. of SU(N)
-// if(N==2) return 12.0;
-// if(N==3) return 12.0;
-// std::cout << "Error! SU(N), N != (2 or 3) used for n_F in ";
-// std::cout << "GroupInvariants.h but not defined." << std::endl;
-// return 0.0;
-// }
-// double n_g() { // Number of fermion generations (only used in gauge boson HighCMatching)
-// return 3.0;
-// }
-// double T_F(int N) { return 0.5+0.0*N; } // I believe T(N) is such that tr(t^a t^b) = T delta(ab)
-// // 0.0*N included to stop receiving a stupid warning.
-
-// double t_S(int N) { return 0.5+0.0*N; } // Analog of T_F but for scalars.
-// // 0.0*N included to stop receiving a stupid warning.
-
-// }
-
-// namespace WeylLow {
-
-// double n_S(int N) { // Number of complex scalars in the fundamental rep. of SU(N)
-// if(N==1) return 0.0;
-// if(N==3) return 0.0;
-// std::cout << "Error! SU(N), N != (1 or 3) used for n_S in ";
-// std::cout << "GroupInvariants.h but not defined." << std::endl;
-// return 0.0;
-// }
-// double n_F(int N) { // Number of Weyl Fermions in the fundamental rep. of SU(N)
-// if(N==3) return 10.0;
-// if(N==1) return 2.0;
-// std::cout << "Error! SU(N), N != (1 or 3) used for n_F in ";
-// std::cout << "GroupInvariants.h but not defined." << std::endl;
-// return 0.0;
-// }
-// double T_F(int N) { return 0.5+0.0*N; } // I believe T(N) is such that tr(t^a t^b) = T delta(ab)
-// // 0.0*N included to stop receiving a stupid warning.
-
-// double t_S(int N) { return 0.5+0.0*N; } // Analog of T_F but for scalars.
-// // 0.0*N included to stop receiving a stupid warning.
-
-// }
-
-
-// #endif // GROUP_INVARIANTS_H
-
-
-
}
#endif // HERWIG_GroupInvariants_H
diff --git a/MatrixElement/EW/HighEnergyMatching.cc b/MatrixElement/EW/HighEnergyMatching.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/HighEnergyMatching.cc
@@ -0,0 +1,1171 @@
+// -*- C++ -*-
+//
+// HighEnergyMatching.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 "HighEnergyMatching.h"
+#include "ElectroWeakReweighter.h"
+#include "GroupInvariants.h"
+#include <boost/numeric/ublas/operation.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
+using namespace Herwig;
+using namespace HighEnergyMatching;
+using namespace GroupInvariants;
+using namespace EWProcess;
+
+namespace {
+
+/**
+ * \f$M_N\f$, this matrix is used for identical particle scattering
+ */
+boost::numeric::ublas::matrix<Complex> M_N(unsigned int suN) {
+ double N(suN);
+ boost::numeric::ublas::matrix<Complex> M(2,2);
+ M(0,0) = -1.0/N;
+ M(0,1) = 2.0;
+ M(1,0) = 0.5 - 1.0/(2.0*N*N);
+ M(1,1) = 1.0/N;
+ return M;
+}
+
+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);
+ }
+ }
+ }
+}
+
+}
+
+boost::numeric::ublas::matrix<complex<InvEnergy2> >
+HighEnergyMatching::highEnergyMatching(Energy highScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ EWProcess::Process process,
+ bool oneLoop, bool includeAlphaS2) {
+ switch (process) {
+ case QQQQ: case QQQQiden:
+ case QtQtQQ: case QQUU:
+ case QtQtUU: case QQtRtR:
+ case QQDD: case QtQtDD:
+ case QQLL: case QQEE:
+ case UUUU: case UUUUiden:
+ case tRtRUU: case UUDD:
+ case tRtRDD: case UULL:
+ case UUEE: case DDDD:
+ case DDDDiden: case DDLL:
+ case DDEE: case LLLL:
+ case LLLLiden: case LLEE:
+ case EEEE: case EEEEiden:
+ return SpinHalfMatching(highScale,s,t,u,process,oneLoop,includeAlphaS2);
+ break;
+ case QQWW:
+ case QQWG:
+ case QQBG:
+ case QQGG:
+ case QtQtGG:
+ case LLWW:
+ case UUBB:
+ case UUBG:
+ case UUGG:
+ case tRtRGG:
+ case DDBB:
+ case DDBG:
+ case DDGG:
+ case EEBB:
+ return Spin1HighMatching(highScale,s,t,u,process,oneLoop,includeAlphaS2);
+ break;
+ case QQPhiPhi:
+ case LLPhiPhi:
+ case UUPhiPhi:
+ case DDPhiPhi:
+ case EEPhiPhi:
+ return Spin0HighMatching(highScale,s,t,u,process,oneLoop,includeAlphaS2);
+ break;
+ default:
+ assert(false);
+ break;
+ }
+}
+
+boost::numeric::ublas::matrix<complex<InvEnergy2> >
+HighEnergyMatching::SpinHalfMatching(Energy highScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ EWProcess::Process process,
+ bool oneLoop, bool includeAlphaS2) {
+ using Constants::pi;
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highC;
+ Energy Q = highScale;
+ double a1 = ElectroWeakReweighter::coupling()->a1(Q);
+ double a2 = ElectroWeakReweighter::coupling()->a2(Q);
+ double a3 = ElectroWeakReweighter::coupling()->a3(Q);
+ double y_t = ElectroWeakReweighter::coupling()->y_t(Q);
+ unsigned int order = !oneLoop ? 0 : 1;
+ double Yi(0.),Yf(0.);
+
+ Complex Ls = MinusLog(-s/(Q*Q));
+
+ double C_A_2 = C_A(2);
+ double C_A_3 = C_A(3);
+ double C_F_2 = C_F(2);
+ double C_F_3 = C_F(3);
+ double C_d_2 = C_d(2);
+ double C_d_3 = C_d(3);
+ double C_1_2 = C_1(2);
+ double C_1_3 = C_1(3);
+ Complex W = WFunction(Q,s);
+ Complex X_2_st = XNFunction(2,Q,s,t);
+ //Complex X_2_su = XNFunction(2,Q,s,u);
+ Complex X_3_st = XNFunction(3,Q,s,t);
+ Complex X_3_su = XNFunction(3,Q,s,u);
+ Complex PI1 = PI1_function(Q,s);
+ Complex fTilde_st = fTildeFunction(Q,s,t);
+ Complex fTilde_su = fTildeFunction(Q,s,u);
+
+ switch (process) {
+
+ case QQQQ:
+ // NOTE this 4x1 column vector highC is given by (C_11,C_21,C_12,C_22),
+ // where C_12 is the coeff. for the term that transforms under SU(2) but not SU(3)
+ Yi = Yf = 1./6.;
+ highC.resize(4,1);
+ highC(0,0) = ZERO;
+ highC(2,0) = 4.0*pi*a2 / s;
+ highC(1,0) = 4.0*pi*a3 / s;
+ highC(3,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(-2.0*a2*a3*fTilde_st);
+ highC(2,0) += (1.0/s)*(a2*a2*(X_2_st-(C_d_2+C_A_2)/4.0*fTilde_st) +
+ 2.0*(a1*a2*Yi*Yf+a2*a3*C_F_3)*W -
+ 2.0*a1*a2*Yi*Yf*fTilde_st);
+ highC(1,0) += (1.0/s)*(2.0*(a1*a3*Yi*Yf+a2*a3*C_F_2)*W -
+ 2.0*a1*a3*Yi*Yf*fTilde_st);
+ highC(3,0) += (1.0/s)*(-1.0*(a2*a2*C_1_2+a1*a1*Yi*Yi*Yf*Yf)*fTilde_st +
+ a1*a1*Yi*Yf*PI1 +
+ 2.0*(a1*a3*Yi*Yf*C_F_3+a1*a2*Yi*Yf*C_F_2+a1*a1*Yi*Yi*Yf*Yf)*W);
+ if (includeAlphaS2) {
+ highC(1,0) += (1.0/s)*(a3*a3*(X_3_st-(C_d_3+C_A_3)/4.0*fTilde_st));
+
+ highC(3,0) += (1.0/s)*(-1.0*(a3*a3*C_1_3)*fTilde_st);
+ }
+ }
+ break;
+
+ case QQQQiden:
+ {
+ Process parentCase = QQQQ;
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st(4,1);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highCs_ts_2x2(2,2);
+ highCs_ts_2x2(0,0) = highCs_ts(0,0);
+ highCs_ts_2x2(1,0) = highCs_ts(1,0);
+ highCs_ts_2x2(0,1) = highCs_ts(2,0);
+ highCs_ts_2x2(1,1) = highCs_ts(3,0);
+ boost::numeric::ublas::matrix<Complex> temp(2,2);
+ temp = boost::numeric::ublas::trans(M_N(3));
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st_2x2(2,2),temp2(2,2);
+ axpy_prod_local(highCs_ts_2x2,temp,temp2);
+ axpy_prod_local(M_N(2),temp2,highCt_st_2x2);
+ highCt_st(0,0) = highCt_st_2x2(0,0);
+ highCt_st(1,0) = highCt_st_2x2(1,0);
+ highCt_st(2,0) = highCt_st_2x2(0,1);
+ highCt_st(3,0) = highCt_st_2x2(1,1);
+ highC = highCs_st + highCt_st;
+ }
+ break;
+ case QtQtQQ:
+ {
+ highC.resize(4,1);
+ Process parentCase = QQQQ;
+ highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ double Y = 1.0/6.0; // Hypercharge of the non-3rd-gen doublet (still a quark doublet).
+ if (order >= 1) {
+ highC(2,0) += y_t*y_t*a2/(4.0*pi*s)*(3.0/2.0-0.5*Ls);
+ highC(1,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0/2.0-0.5*Ls);
+ highC(3,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(-5.0/12.0-1.0/12.0*Ls);
+ }
+ }
+ break;
+ case QQUU:
+ Yi = 1./6.; Yf = 2./3.;
+ highC.resize(2,1);
+ highC(0,0) = 4.0*pi*a3 / s;
+ highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*((a1*a3*(Yi*Yi+Yf*Yf)+a2*a3*C_F_2)*W +
+ 2.0*a1*a3*Yi*Yf*fTilde_su);
+ highC(1,0) += (1.0/s)*(a1*a1*Yi*Yf*PI1 +
+ (a1*a2*Yi*Yf*C_F_2+2.0*a1*a3*Yi*Yf*C_F_3+
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ if (includeAlphaS2) {
+ highC(0,0) += (1.0/s)*(a3*a3*(X_3_su+(C_d_3-C_A_3)/4.0*fTilde_su));
+ highC(1,0) += (1.0/s)*((a3*a3*C_1_3+a1*a1*Yi*Yi*Yf*Yf)*fTilde_su);
+ }
+ }
+ break;
+ case QtQtUU:
+ {
+ highC.resize(2,1);
+ Process parentCase = QQUU;
+ highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ double Y = 2./3.;
+ if (order >= 1) {
+ highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0/2.0-0.5*Ls);
+ highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(-5.0/12.0-1.0/12.0*Ls);
+ }
+ }
+ break;
+ case QQtRtR:
+ {
+ highC.resize(2,1);
+ Process parentCase = QQUU;
+ highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ double Y = 1.0/6.0;
+ if (order >= 1) {
+ highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0-Ls);
+ highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(5.0/3.0-2.0/3.0*Ls);
+ }
+ }
+ break;
+ case QQDD:
+ Yi = 1./6.; Yf = -1./3.;
+ highC.resize(2,1);
+ highC(0,0) = 4.0*pi*a3 / s;
+ highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*((a1*a3*(Yi*Yi+Yf*Yf)+a2*a3*C_F_2)*W +
+ 2.0*a1*a3*Yi*Yf*fTilde_su);
+ highC(1,0) += (1.0/s)*(a1*a1*Yi*Yf*PI1 +
+ (a1*a2*Yi*Yf*C_F_2+2.0*a1*a3*Yi*Yf*C_F_3+
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ if (includeAlphaS2) {
+ highC(0,0) += (1.0/s)*(a3*a3*(X_3_su+(C_d_3-C_A_3)/4.0*fTilde_su));
+ highC(1,0) += (1.0/s)*((a3*a3*C_1_3+a1*a1*Yi*Yi*Yf*Yf)*fTilde_su);
+ }
+ }
+ break;
+ case QtQtDD:
+ {
+ highC.resize(2,1);
+ Process parentCase = QQDD;
+ highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ double Y = -1./3.;
+ if (order >= 1) {
+ highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0/2.0-0.5*Ls);
+ highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(-5.0/12.0-1.0/12.0*Ls);
+ }
+ }
+ break;
+ case QQLL:
+ Yi = 1./6.; Yf = -1./2.;
+ highC.resize(2,1);
+ highC(0,0) = 4.0*pi*a2 / s;
+ highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(a2*a2*(X_2_st-(C_d_2+C_A_2)/4.0*fTilde_st) +
+ (a2*a3*C_F_3 + a1*a2*(Yi*Yi+Yf*Yf))*W -
+ 2.0*a1*a2*Yi*Yf*fTilde_st);
+ highC(1,0) += (1.0/s)*(-1.0*(a2*a2*C_1_2+a1*a1*Yi*Yi*Yf*Yf)*fTilde_st +
+ a1*a1*Yi*Yf*PI1 +
+ (a1*a3*Yi*Yf*C_F_3+2.0*a1*a2*Yi*Yf*C_F_2+
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ }
+ break;
+ case QQEE:
+ Yi = 1./6.; Yf = -1.;
+ highC.resize(1,1);
+ highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(a1*a1*Yi*Yi*Yf*Yf*fTilde_su + a1*a1*Yi*Yf*PI1 +
+ (a1*a3*Yi*Yf*C_F_3 + a1*a2*Yi*Yf*C_F_2 +
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ }
+ break;
+ case UUUU:
+ Yi = Yf = 2./3.;
+ highC.resize(2,1);
+ highC(0,0) = 4.0*pi*a3 / s;
+ highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(-2.0*a1*a3*Yi*Yf*fTilde_st +
+ a1*a3*(Yi*Yi+Yf*Yf)*W);
+ highC(1,0) += (1.0/s)*(-1.0*(a1*a1*Yi*Yi*Yf*Yf)*fTilde_st +
+ a1*a1*Yi*Yf*PI1 +
+ (2.0*a1*a3*Yi*Yf*C_F_3+a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ if (includeAlphaS2) {
+ highC(0,0) += (1.0/s)*(a3*a3*(X_3_st-(C_d_3+C_A_3)/4.0*fTilde_st));
+ highC(1,0) += (1.0/s)*(-1.0*(a3*a3*C_1_3)*fTilde_st);
+ }
+ }
+ break;
+ case UUUUiden:
+ {
+ Process parentCase = UUUU;
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st;
+ axpy_prod_local(M_N(3),highCs_ts,highCt_st);
+ highC = highCs_st + highCt_st;
+ }
+ break;
+ case tRtRUU:
+ {
+ highC.resize(2,1);
+ Process parentCase = UUUU;
+ highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ double Y = 2./3.;
+ if (order >= 1) {
+ highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0-Ls);
+ highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(5.0/3.0-2.0/3.0*Ls);
+ }
+ }
+ break;
+ case UUDD:
+ Yi = 2./3.; Yf = -1./3.;
+ highC.resize(2,1);
+ highC(0,0) = 4.0*pi*a3 / s;
+ highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(-2.0*a1*a3*Yi*Yf*fTilde_st +
+ a1*a3*(Yi*Yi+Yf*Yf)*W);
+ highC(1,0) += (1.0/s)*(-1.0*(a1*a1*Yi*Yi*Yf*Yf)*fTilde_st +
+ a1*a1*Yi*Yf*PI1 +
+ (2.0*a1*a3*Yi*Yf*C_F_3+a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ if (includeAlphaS2) {
+ highC(0,0) += (1.0/s)*(a3*a3*(X_3_st-(C_d_3+C_A_3)/4.0*fTilde_st));
+ highC(1,0) += (1.0/s)*(-1.0*(a3*a3*C_1_3)*fTilde_st);
+ }
+ }
+ break;
+ case tRtRDD:
+ {
+ highC.resize(2,1);
+ Process parentCase = UUDD;
+ highC = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ double Y = -1./3.;
+ if (order >= 1) {
+ highC(0,0) += y_t*y_t*a3/(4.0*pi*s)*(1.0-Ls);
+ highC(1,0) += y_t*y_t*a1*Y/(4.0*pi*s)*(5.0/3.0-2.0/3.0*Ls);
+ }
+ }
+ break;
+ case UULL:
+ Yi = 2./3.; Yf = -0.5;
+ highC.resize(1,1);
+ highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(a1*a1*Yi*Yi*Yf*Yf*fTilde_su + a1*a1*Yi*Yf*PI1 +
+ (a1*a3*Yi*Yf*C_F_3 + a1*a2*Yi*Yf*C_F_2 +
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ }
+ break;
+ case UUEE:
+ Yi = 2./3.; Yf = -1.;
+ highC.resize(1,1);
+ highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(-1.0*a1*a1*Yi*Yi*Yf*Yf*fTilde_st + a1*a1*Yi*Yf*PI1 +
+ (a1*a3*Yi*Yf*C_F_3 +
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ }
+ break;
+ case DDDD:
+ Yi = Yf = -1./3.;
+ highC.resize(2,1);
+ highC(0,0) = 4.0*pi*a3 / s;
+ highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(-2.0*a1*a3*Yi*Yf*fTilde_st +
+ a1*a3*(Yi*Yi+Yf*Yf)*W);
+ highC(1,0) += (1.0/s)*(-1.0*(a1*a1*Yi*Yi*Yf*Yf)*fTilde_st +
+ a1*a1*Yi*Yf*PI1 +
+ (2.0*a1*a3*Yi*Yf*C_F_3+a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ if (includeAlphaS2) {
+ highC(0,0) += (1.0/s)*(a3*a3*(X_3_st-(C_d_3+C_A_3)/4.0*fTilde_st));
+ highC(1,0) += (1.0/s)*(-1.0*(a3*a3*C_1_3)*fTilde_st);
+ }
+ }
+ break;
+ case DDDDiden:
+ {
+ Process parentCase = DDDD;
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st;
+ axpy_prod_local(M_N(3),highCs_ts,highCt_st);
+ highC = highCs_st + highCt_st;
+ }
+ break;
+ case DDLL:
+ Yi = -1./3.; Yf = -0.5;
+ highC.resize(1,1);
+ highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(a1*a1*Yi*Yi*Yf*Yf*fTilde_su + a1*a1*Yi*Yf*PI1 +
+ (a1*a3*Yi*Yf*C_F_3 + a1*a2*Yi*Yf*C_F_2 +
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ }
+ break;
+ case DDEE:
+ Yi = -1./3.; Yf = -1.;
+ highC.resize(1,1);
+ highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(-1.0*a1*a1*Yi*Yi*Yf*Yf*fTilde_st + a1*a1*Yi*Yf*PI1 +
+ (a1*a3*Yi*Yf*C_F_3 +
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ }
+ break;
+ case LLLL:
+ Yi = Yf = -0.5;
+ highC.resize(2,1);
+ highC(0,0) = 4.0*pi*a2 / s;
+ highC(1,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(a2*a2*(X_2_st-(C_d_2+C_A_2)/4.0*fTilde_st) +
+ 2.0*a1*a2*Yi*Yf*W -
+ 2.0*a1*a2*Yi*Yf*fTilde_st);
+ highC(1,0) += (1.0/s)*(-1.0*(a2*a2*C_1_2+a1*a1*Yi*Yi*Yf*Yf)*fTilde_st +
+ a1*a1*Yi*Yf*PI1 +
+ 2.0*(a1*a2*Yi*Yf*C_F_2+a1*a1*Yi*Yi*Yf*Yf)*W);
+ }
+ break;
+ case LLLLiden:
+ {
+ Process parentCase = LLLL;
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highCt_st;
+ axpy_prod_local(M_N(2), highCs_ts, highCt_st);
+ highC = highCs_st + highCt_st;
+ }
+ break;
+ case LLEE:
+ Yi = -0.5; Yf = -1.;
+ highC.resize(1,1);
+ highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(a1*a1*Yi*Yi*Yf*Yf*fTilde_su + a1*a1*Yi*Yf*PI1 +
+ (a1*a2*Yi*Yf*C_F_2 +
+ a1*a1*(Yi*Yi*Yi*Yf+Yf*Yf*Yf*Yi))*W);
+ }
+ break;
+ case EEEE:
+ Yi = Yf = -1.;
+ highC.resize(1,1);
+ highC(0,0) = 4.0*pi*a1*Yi*Yf / s;
+ if (order >= 1) {
+ highC(0,0) += (1.0/s)*(-1.0*a1*a1*Yi*Yi*Yf*Yf*fTilde_st + a1*a1*Yi*Yf*PI1 +
+ 2.0*a1*a1*Yi*Yi*Yf*Yf*W);
+ }
+ break;
+ case EEEEiden:
+ {
+ Process parentCase = EEEE;
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_st = highEnergyMatching(Q,s,t,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCs_ts = highEnergyMatching(Q,t,s,u,parentCase,oneLoop,includeAlphaS2);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highCt_st = highCs_ts;
+ highC = highCs_st + highCt_st;
+ }
+ break;
+ default:
+ assert(false);
+ }
+ return highC;
+}
+
+boost::numeric::ublas::matrix<complex<InvEnergy2> >
+HighEnergyMatching::Spin1HighMatching(Energy highScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ EWProcess::Process process,
+ bool oneLoop, bool includeAlphaS2) {
+ using Constants::pi;
+
+ unsigned int order = !oneLoop ? 0 : 1;
+ // (If crossed graphs, swap s and t here)
+ Complex L_s = MinusLog(-s/(highScale*highScale));
+ Complex L_t = MinusLog(-t/(highScale*highScale));
+ Complex L_u = MinusLog(-u/(highScale*highScale));
+ Complex L_s2 = L_s*L_s;
+ Complex L_t2 = L_t*L_t;
+ Complex L_u2 = L_u*L_u;
+
+ // Tree-Level:
+ // Topology types defined. Note each of these is a vector of 5 entries. They are the coefficients
+ // for the dirac structures M_0, M_1, M_4, M_5, and M_6 for vector boson production.
+ boost::numeric::ublas::vector<complex<InvEnergy2> > R1(5);
+ for(unsigned int ix=0;ix<5;++ix) R1[ix] = ZERO;
+ R1[0] = -1.0/t;
+ R1[1] = -2.0/t;
+ R1[2] = ZERO;
+ R1[3] = ZERO;
+ R1[4] = ZERO;
+ boost::numeric::ublas::vector<complex<InvEnergy2> > R1_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) R1_bar[ix] = ZERO;
+ R1_bar[0] = -1.0/u;
+ boost::numeric::ublas::vector<complex<InvEnergy2> > R2(5);
+ for(unsigned int ix=0;ix<5;++ix) R2[ix] = ZERO;
+ R2[1] = -1.0/s*2.0;
+ // Topologies T1:
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T1a(5);
+ for(unsigned int ix=0;ix<5;++ix) T1a[ix] = ZERO;
+ T1a[0] = 1.0/(t*u)*(-3.0*t*L_s2-(s+4.0*t)*L_t2+2.0*(s+4.0*t)*L_s*L_t+2.0*u*L_t-
+ pi*pi*(7.0/6.0*s+25.0/6.0*t)-4.0*u);
+ T1a[1] = 1.0/(u*u*t*s)*(0.5*t*(9.0*s*s+14.0*s*t+7.0*t*t)*L_s2+s*(2.0*s+t)*(s+2.0*t)*L_t2-
+ 2.0*(2.0*s*s*s+9.0*s*s*t+10.0*s*t*t+4.0*t*t*t)*L_s*L_t-
+ 2.0*t*t*u*L_s-2.0*u*s*(2.0*s+3.0*t)*L_t+
+ pi*pi*(7.0/3.0*s*s*s+125.0/12.0*s*s*t+71.0/6.0*s*t*t+
+ 19.0/4.0*t*t*t)-
+ 8.0*s*s*s-20.0*s*s*t-16.0*s*t*t-4.0*t*t*t);
+ T1a[2] = 1.0/(t*u*u)*(-t*(3.0*s+4.0*t)*L_s2-(s*s+5.0*s*t+5.0*t*t)*L_t2+
+ 2.0*t*(3.0*s+4.0*t)*L_s*L_t+2.0*u*t*(2.0*s+t)*L_s/s-
+ 2.0*u*t*L_t+pi*pi*(s*s/6.0-8.0/3.0*s*t-23.0/6.0*t*t)+
+ 4.0*t*t*t/s+4.0*s*t+8.0*t*t);
+ T1a[3] = T1a[2];
+ T1a[4] = GeV2/(t*u*u*u)*(-4.0*t*(s+2.0*t)*(L_s-L_t)*(L_s-L_t)+
+ 4.0*u*(3.0*s+5.0*t)*(L_s-L_t)-4.0*pi*pi*t*(s+2.0*t)-4.0*u*u);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T1b(5);
+ for(unsigned int ix=0;ix<5;++ix) T1b[ix] = ZERO;
+ T1b[0] = 1.0/(t*u*s*s)*(-s*t*(2.0*s+3.0*t)*L_u2+s*u*(s+3.0*t)*L_t2+
+ 2.0*s*(s*s+3.0*s*t+3.0*t*t)*L_u*L_t+s*s*t*L_u+s*s*u*L_t-
+ pi*pi*(7.0/6.0*s*s*s+3.0*s*s*t+3.0*s*t*t)+2.0*s*s*s);
+ T1b[1] = 1.0/(t*s*s*u)*(3.0*s*t*u*L_u2+s*u*(2.0*s+3.0*t)*L_t2-
+ 2.0*s*u*(2.0*s+3.0*t)*L_u*L_t+2.0*s*s*u*L_t-
+ pi*pi*(7.0/3.0*s*s*s+16.0/3.0*s*s*t+3.0*s*t*t)+
+ 4.0*s*s*s+4.0*s*s*t);
+ T1b[2] = 1.0/(t*u*s*s)*(-3.0*s*t*u*(L_u-L_t)*(L_u-L_t)+4.0*s*s*t*L_u+4.0*s*s*u*L_t+
+ pi*pi*(3.0*s*s*t+3.0*s*t*t)+8.0*s*s*s);
+ T1b[3] = 1.0/(t*u*s*s)*(s*t*(2.0*s+3.0*t)*L_u2-s*u*(s+3.0*t)*L_t2+6.0*s*t*u*L_u*L_t+
+ pi*pi*(-1.0/6.0*s*s*s+3.0*s*s*t+3.0*s*t*t));
+ T1b[4] = 12.0*GeV2/(t*u)*(L_t-L_u);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T1c(5);
+ for(unsigned int ix=0;ix<5;++ix) T1c[ix] = ZERO;
+ T1c[0] = 1.0/t*(2.0*L_s*L_t-7.0*pi*pi/6.0-L_t2);
+ T1c[1] = 1.0/(t*u*u)*(s*t*L_s2-(2.0*s*s+3.0*s*t+2.0*t*t)*L_t2+
+ 2.0*(2.0*s*s+3.0*s*t+2.0*t*t)*L_s*L_t+2.0*t*u*(L_s-L_t)-
+ pi*pi*(7.0/3.0*s*s+11.0/3.0*s*t+7.0/3.0*t*t));
+ T1c[2] = 1.0/(t*u*u)*(t*(3.0*s+2.0*t)*(L_s-L_t)*(L_s-L_t)+2.0*u*t*L_s+
+ 2.0*u*(2.0*s+t)*L_t+pi*pi*t*(3.0*s+2.0*t)+8.0*u*u);
+ T1c[3] = T1c[2];
+ T1c[4] = GeV2/(t*u*u*u)*(4.0*t*(2.0*s+t)*(L_s-L_t)*(L_s-L_t)-4.0*u*(3.0*s+t)*(L_s-L_t)+
+ 4.0*pi*pi*t*(2.0*s+t)-4.0*u*u);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T1d(5);
+ for(unsigned int ix=0;ix<5;++ix) T1d[ix] = ZERO;
+ T1d[2] = 1.0/s*(-2.0*L_s+4.0);
+ T1d[3] = T1d[2];
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T1a_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) T1a_bar[ix] = ZERO;
+ T1a_bar[0] = 1.0/(u*t)*(-3.0*u*L_s2-(s+4.0*u)*L_u2+2.0*(s+4.0*u)*L_s*L_u+2.0*t*L_u-
+ pi*pi*(7.0/6.0*s+25.0/6.0*u)-4.0*t);
+ T1a_bar[1] = 2.0*T1a_bar[0] -
+ 1.0/(t*t*u*s)*(0.5*u*(9.0*s*s+14.0*s*u+7.0*u*u)*L_s2+s*(2.0*s+u)*(s+2.0*u)*L_u2-
+ 2.0*(2.0*s*s*s+9.0*s*s*u+10.0*s*u*u+4.0*u*u*u)*L_s*L_u-
+ 2.0*u*u*t*L_s-2.0*t*s*(2.0*s+3.0*u)*L_u+
+ pi*pi*(7.0/3.0*s*s*s+125.0/12.0*s*s*u+71.0/6.0*s*u*u+
+ 19.0/4.0*u*u*u)-
+ 8.0*s*s*s-20.0*s*s*u-16.0*s*u*u-4.0*u*u*u);
+ T1a_bar[2] = 1.0/(u*t*t)*(-u*(3.0*s+4.0*u)*L_s2-(s*s+5.0*s*u+5.0*u*u)*L_u2+
+ 2.0*u*(3.0*s+4.0*u)*L_s*L_u+2.0*t*u*(2.0*s+u)*L_s/s-
+ 2.0*t*u*L_u+pi*pi*(s*s/6.0-8.0/3.0*s*u-23.0/6.0*u*u)+
+ 4.0*u*u*u/s+4.0*s*u+8.0*u*u);
+ T1a_bar[3] = T1a_bar[2];
+ T1a_bar[4] = -GeV2/(u*t*t*t)*(-4.0*u*(s+2.0*u)*(L_s-L_u)*(L_s-L_u)+
+ 4.0*t*(3.0*s+5.0*u)*(L_s-L_u)-4.0*pi*pi*u*(s+2.0*u)-4.0*t*t);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T1b_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) T1b_bar[ix] = ZERO;
+ T1b_bar[0] = 1.0/(u*t*s*s)*(-s*u*(2.0*s+3.0*u)*L_t2+s*t*(s+3.0*u)*L_u2+
+ 2.0*s*(s*s+3.0*s*u+3.0*u*u)*L_t*L_u+
+ s*s*u*L_t+s*s*t*L_u-
+ pi*pi*(7.0/6.0*s*s*s+3.0*s*s*u+3.0*s*u*u)+2.0*s*s*s);
+ T1b_bar[1] = 2.0*T1b_bar[0] -
+ 1.0/(u*s*s*t)*(3.0*s*u*t*L_t2+s*t*(2.0*s+3.0*u)*L_u2-
+ 2.0*s*t*(2.0*s+3.0*u)*L_t*L_u+2.0*s*s*t*L_u-
+ pi*pi*(7.0/3.0*s*s*s+16.0/3.0*s*s*u+3.0*s*u*u)+
+ 4.0*s*s*s+4.0*s*s*u);
+ T1b_bar[3] = 1.0/(u*t*s*s)*(-3.0*s*u*t*(L_t-L_u)*(L_t-L_u)+4.0*s*s*u*L_t+4.0*s*s*t*L_u+
+ pi*pi*(3.0*s*s*u+3.0*s*u*u)+8.0*s*s*s);
+ T1b_bar[2] = 1.0/(u*t*s*s)*(s*u*(2.0*s+3.0*u)*L_t2-s*t*(s+3.0*u)*L_u2+6.0*s*u*t*L_t*L_u+
+ pi*pi*(-1.0/6.0*s*s*s+3.0*s*s*u+3.0*s*u*u));
+ T1b_bar[4] = -12.0*GeV2/(u*t)*(L_u-L_t);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T1c_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) T1c_bar[ix] = ZERO;
+ T1c_bar[0] = 1.0/u*(2.0*L_s*L_u-7.0*pi*pi/6.0-L_u2);
+ T1c_bar[1] = 2.0*T1c_bar[0] -
+ 1.0/(u*t*t)*(s*u*L_s2-(2.0*s*s+3.0*s*u+2.0*u*u)*L_u2+
+ 2.0*(2.0*s*s+3.0*s*u+2.0*u*u)*L_s*L_u+2.0*u*t*(L_s-L_u)-
+ pi*pi*(7.0/3.0*s*s+11.0/3.0*s*u+7.0/3.0*u*u));
+ T1c_bar[2] = 1.0/(u*t*t)*(u*(3.0*s+2.0*u)*(L_s-L_u)*(L_s-L_u)+2.0*t*u*L_s+
+ 2.0*t*(2.0*s+u)*L_u+pi*pi*u*(3.0*s+2.0*u)+8.0*t*t);
+ T1c_bar[3] = T1c_bar[2];
+ T1c_bar[4] = -GeV2/(u*t*t*t)*(4.0*u*(2.0*s+u)*(L_s-L_u)*(L_s-L_u)-4.0*t*(3.0*s+u)*(L_s-L_u)+
+ 4.0*pi*pi*u*(2.0*s+u)-4.0*t*t);
+ // Topologies T2:
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T2a(5);
+ for(unsigned int ix=0;ix<5;++ix) T2a[ix] = ZERO;
+ T2a[1] = 1.0/s*(L_s/6.0-11.0/18.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T2b(5);
+ for(unsigned int ix=0;ix<5;++ix) T2b[ix] = ZERO;
+ T2b[1] = 1.0/s*(-2.0/3.0*L_s+22.0/9.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T2c(5);
+ for(unsigned int ix=0;ix<5;++ix) T2c[ix] = ZERO;
+ T2c[1] = 1.0/s*(3.0/2.0*L_s2-17.0/2.0*L_s-pi*pi/4.0+95.0/6.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T2d(5);
+ for(unsigned int ix=0;ix<5;++ix) T2d[ix] = ZERO;
+ T2d[1] = 1.0/s*(-4.0/3.0*L_s+14.0/9.0);
+ // Topologies T3:
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T3a(5);
+ for(unsigned int ix=0;ix<5;++ix) T3a[ix] = ZERO;
+ T3a[0] = 1.0/t*(L_t2-pi*pi/6.0);
+ T3a[1] = 2.0*T3a[0];
+ T3a[3] = 1.0/t*(-L_t2+pi*pi/6.0+2.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T3b(5);
+ for(unsigned int ix=0;ix<5;++ix) T3b[ix] = ZERO;
+ T3b[0] = 1.0/t*(-L_t+4.0);
+ T3b[1] = 2.0*T3b[0];
+ T3b[3] = 1.0/t*(4.0*L_t-10.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T3a_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) T3a_bar[ix] = ZERO;
+ T3a_bar[0] = 1.0/u*(L_u2-pi*pi/6.0);
+ T3a_bar[2] = 1.0/u*(-L_u2+pi*pi/6.0+2.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T3b_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) T3b_bar[ix] = ZERO;
+ T3b_bar[0] = 1.0/u*(-L_u+4.0);
+ T3b_bar[2] = 1.0/u*(4.0*L_u-10.0);
+ // Topologies T4:
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T4a(5);
+ for(unsigned int ix=0;ix<5;++ix) T4a[ix] = ZERO;
+ T4a[0] = 1.0/t*(L_t2-pi*pi/6.0);
+ T4a[1] = 2.0*T4a[0];
+ T4a[2] = 1.0/t*(-L_t2+pi*pi/6.0+2.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T4b(5);
+ for(unsigned int ix=0;ix<5;++ix) T4b[ix] = ZERO;
+ T4b[0] = 1.0/t*(-L_t+4.0);
+ T4b[1] = 2.0*T4b[0];
+ T4b[2] = 1.0/t*(4.0*L_t-10.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T4a_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) T4a_bar[ix] = ZERO;
+ T4a_bar[0] = 1.0/u*(L_u2-pi*pi/6.0);
+ T4a_bar[3] = 1.0/u*(-L_u2+pi*pi/6.0+2.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T4b_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) T4b_bar[ix] = ZERO;
+ T4b_bar[0] = 1.0/u*(-L_u+4.0);
+ T4b_bar[3] = 1.0/u*(4.0*L_u-10.0);
+ // Topologies T5:
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T5a(5);
+ for(unsigned int ix=0;ix<5;++ix) T5a[ix] = ZERO;
+ T5a[1] = 1.0/s*(2.0*L_s-4.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T5b(5);
+ for(unsigned int ix=0;ix<5;++ix) T5b[ix] = ZERO;
+ T5b[1] = 1.0/s*(-2.0*L_s2+6.0*L_s-16.0+pi*pi/3.0);
+ // Topologies T6:
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T6a(5);
+ for(unsigned int ix=0;ix<5;++ix) T6a[ix] = ZERO;
+ T6a[1] = 1.0/s*(-19.0/6.0*L_s+58.0/9.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T6b(5);
+ for(unsigned int ix=0;ix<5;++ix) T6b[ix] = ZERO;
+ T6b[1] = 1.0/s*(-1.0/6.0*L_s+4.0/9.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T6c(5);
+ for(unsigned int ix=0;ix<5;++ix) T6c[ix] = ZERO;
+ T6c[1] = 1.0/s*(2.0/3.0*L_s-16.0/9.0);
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T6d(5);
+ for(unsigned int ix=0;ix<5;++ix) T6d[ix] = ZERO;
+ T6d[1] = 1.0/s*(4.0/3.0*L_s-20.0/9.0);
+ // Topology T7:
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T7(5);
+ for(unsigned int ix=0;ix<5;++ix) T7[ix] = ZERO;
+ T7[0] = 1.0/t*(-L_t+1.0);
+ T7[1] = 2.0*T7[0];
+ boost::numeric::ublas::vector<complex<InvEnergy2> > T7_bar(5);
+ for(unsigned int ix=0;ix<5;++ix) T7_bar[ix] = ZERO;
+ T7_bar[0] = 1.0/u*(-L_u+1.0);
+ // Group Theory Factors / SM parameters needed for matrix elements:
+ double a1 = ElectroWeakReweighter::coupling()->a1(highScale);
+ double a2 = ElectroWeakReweighter::coupling()->a2(highScale);
+ double a3 = ElectroWeakReweighter::coupling()->a3(highScale);
+ double y_t = ElectroWeakReweighter::coupling()->y_t(highScale);
+ // Traces over complex scalars and weyl fermions.
+ double T_CS_3 = 0.0;
+ double T_CS_2 = 0.5;
+ //double T_CS_1 = 0.5;
+ double T_WF_3 = 2.0*3.0;
+ double T_WF_2 = 2.0*3.0;
+ //double T_WF_1 = 10.0/3.0*3.0;
+ double C_A_3 = 3.0;
+ double C_A_2 = 2.0;
+ double C_A_1 = 0.0;
+ double C_F_3 = 4.0/3.0;
+ double C_F_2 = 3.0/4.0;
+ double C_F_1 = 1.0;
+ // This is the coefficient of the delta term in G_TT
+ double G_TT = 0.5;
+ // This is the coeffidient of d^ABC in G_TT (non-zero for SU(3))
+ double G_TT_3_D = 0.25*C_A_3;
+ double G_f = 1.0;
+ // Factors TBD after fermion helicity is specified:
+ double Y_Q(0.), G_Plus_U1(0.);
+ double G_Plus_SU2 = 0.25;
+ double G_Plus_SU3 = 1./6.;
+ double G_Plus_SU3_D = 0.5;
+ double Lambda_Q(0.);
+ // the _s and _ew are the alpha3 and alpha1/2 parts of Lambda_Q
+ double Lambda_Q_s(0.);
+ double Lambda_Q_ew(0.);
+ double rho12(0.), rho13(0.);
+ double rho23 = sqrt(a2*a3);
+ double tRorQ = 1.0;
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highC(1,1);
+ switch (process) {
+ case QQWW: case LLWW:
+ {
+ // Finish Group Theory Factors:
+ if (process==QQWW) {
+ Y_Q = 1.0/6.0;
+ G_Plus_U1 = Y_Q*Y_Q;
+ Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
+ rho12 = Y_Q*sqrt(a1*a2);
+ }
+ else if (process==LLWW) {
+ Y_Q = -1.0/2.0;
+ G_Plus_U1 = Y_Q*Y_Q;
+ Lambda_Q = C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
+ rho12 = Y_Q*sqrt(a1*a2);
+ }
+ highC.resize(5,5);
+ for (int i=0; i<5; i++) {
+ highC(0,i) = G_Plus_SU2*(4.0*pi*a2)*(R1[i]+R1_bar[i]);
+ highC(1,i) = G_f*4.0*pi*a2*(-0.5*R1[i]+0.5*R1_bar[i]-R2[i]);
+ highC(2,i) = rho12*4.0*pi*(R1[i]+R1_bar[i]);
+ highC(3,i) = rho12*4.0*pi*(R1[i]+R1_bar[i]);
+ highC(4,i) = G_Plus_U1*(4.0*pi*a1)*(R1[i]+R1_bar[i]);
+ if (order>=1) {
+ highC(0,i) += G_Plus_SU2*((-0.5*a2*a2*C_A_2)*(T1b[i]+T1b_bar[i])+
+ a2*(Lambda_Q-a2*C_A_2)*(T1c[i]+T1c_bar[i])+
+ 0.5*a2*a2*C_A_2*(T3a[i]+T3a_bar[i])+
+ a2*(Lambda_Q-0.5*a2*C_A_2)*(T3b[i]+T3b_bar[i])+
+ 0.5*a2*a2*C_A_2*(T4a[i]+T4a_bar[i])+
+ a2*(Lambda_Q-0.5*a2*C_A_2)*(T4b[i]+T4b_bar[i])+
+ a2*Lambda_Q*(T7[i]+T7_bar[i])) +
+ G_TT*(-a2*a2*(T1a[i]+T1a_bar[i])+a2*a2*(T1b[i]+T1b_bar[i])+
+ a2*a2*(T1c[i]+T1c_bar[i])+2.0*a2*a2*T1d[i]);
+ highC(1,i) += G_f*(0.25*a2*a2*C_A_2*(T1a[i]-T1a_bar[i])+
+ a2*(0.25*a2*C_A_2-0.5*Lambda_Q)*(T1c[i]-T1c_bar[i])+
+ 0.5*a2*a2*C_A_2*T2a[i]+a2*a2*T_CS_2*T2b[i]-
+ 0.5*a2*a2*C_A_2*T2c[i]+a2*a2*T_WF_2*T2d[i]-
+ 0.25*a2*a2*C_A_2*(T3a[i]-T3a_bar[i])-
+ 0.5*a2*(Lambda_Q-0.5*a2*C_A_2)*(T3b[i]-T3b_bar[i])-
+ 0.25*a2*a2*C_A_2*(T4a[i]-T4a_bar[i])-
+ 0.5*a2*(Lambda_Q-0.5*a2*C_A_2)*(T4b[i]-T4b_bar[i])+
+ 0.5*a2*a2*C_A_2*T5a[i]+a2*(Lambda_Q-0.5*a2*C_A_2)*T5b[i]+
+ a2*a2*C_A_2*T6a[i]+a2*a2*C_A_2*T6b[i]+
+ a2*a2*T_CS_2*T6c[i]+a2*a2*T_WF_2*T6d[i]-
+ 0.5*a2*Lambda_Q*(T7[i]-T7_bar[i]));
+ highC(2,i) += rho12*(-0.5*a1*C_A_1*T1b[i]-0.5*a2*C_A_2*T1b_bar[i]+
+ (Lambda_Q-0.5*a1*C_A_1-0.5*a2*C_A_2)*(T1c[i]+T1c_bar[i])+
+ 0.5*a1*C_A_1*T3a[i]+0.5*a2*C_A_2*T3a_bar[i]+
+ (Lambda_Q-0.5*a1*C_A_1)*T3b[i]+(Lambda_Q-0.5*a2*C_A_2)*T3b_bar[i]+
+ 0.5*a2*C_A_2*T4a[i]+0.5*a1*C_A_1*T4a_bar[i]+
+ (Lambda_Q-0.5*a2*C_A_2)*T4b[i]+(Lambda_Q-0.5*a1*C_A_1)*T4b_bar[i]+
+ Lambda_Q*(T7[i]+T7_bar[i]));
+ highC(3,i) += rho12*(-0.5*a2*C_A_2*T1b[i]-0.5*a1*C_A_1*T1b_bar[i]+
+ (Lambda_Q-0.5*a2*C_A_2-0.5*a1*C_A_1)*(T1c[i]+T1c_bar[i])+
+ 0.5*a2*C_A_2*T3a[i]+0.5*a1*C_A_1*T3a_bar[i]+
+ (Lambda_Q-0.5*a2*C_A_2)*T3b[i]+(Lambda_Q-0.5*a1*C_A_1)*T3b_bar[i]+
+ 0.5*a1*C_A_1*T4a[i]+0.5*a2*C_A_2*T4a_bar[i]+
+ (Lambda_Q-0.5*a1*C_A_1)*T4b[i]+(Lambda_Q-0.5*a2*C_A_2)*T4b_bar[i]+
+ Lambda_Q*(T7[i]+T7_bar[i]));
+ highC(4,i) += G_Plus_U1*((-0.5*a1*a1*C_A_1)*(T1b[i]+T1b_bar[i])+
+ a1*(Lambda_Q-a1*C_A_1)*(T1c[i]+T1c_bar[i])+
+ 0.5*a1*a1*C_A_1*(T3a[i]+T3a_bar[i])+
+ a1*(Lambda_Q-0.5*a1*C_A_1)*(T3b[i]+T3b_bar[i])+
+ 0.5*a1*a1*C_A_1*(T4a[i]+T4a_bar[i])+
+ a1*(Lambda_Q-0.5*a1*C_A_1)*(T4b[i]+T4b_bar[i])+
+ a1*Lambda_Q*(T7[i]+T7_bar[i]));
+ }
+ }
+ }
+ break;
+ case UUBB: case DDBB: case EEBB:
+ {
+ // Finish Group Theory Factors:
+ if (process==UUBB) {
+ Y_Q = 2.0/3.0;
+ G_Plus_U1 = Y_Q*Y_Q;
+ Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
+ }
+ else if (process==DDBB) {
+ Y_Q = -1.0/3.0;
+ G_Plus_U1 = Y_Q*Y_Q;
+ Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
+ }
+ else if (process==EEBB) {
+ Y_Q = -1.0;
+ G_Plus_U1 = Y_Q*Y_Q;
+ Lambda_Q = Y_Q*Y_Q*C_F_1*a1;
+ }
+ highC.resize(1,5);
+ for (int i=0; i<5; i++) {
+ highC(0,i) = G_Plus_U1*(4.0*pi*a1)*(R1[i]+R1_bar[i]);
+ if (order>=1) {
+ highC(0,i) += G_Plus_U1*((-0.5*a1*a1*C_A_1)*(T1b[i]+T1b_bar[i])+
+ a1*(Lambda_Q-a1*C_A_1)*(T1c[i]+T1c_bar[i])+
+ 0.5*a1*a1*C_A_1*(T3a[i]+T3a_bar[i])+
+ a1*(Lambda_Q-0.5*a1*C_A_1)*(T3b[i]+T3b_bar[i])+
+ 0.5*a1*a1*C_A_1*(T4a[i]+T4a_bar[i])+
+ a1*(Lambda_Q-0.5*a1*C_A_1)*(T4b[i]+T4b_bar[i])+
+ a1*Lambda_Q*(T7[i]+T7_bar[i]));
+ }
+ }
+ }
+ break;
+ case QQWG:
+ {
+ // Finish Group Theory Factors:
+ Y_Q = 1./6.;
+ Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
+
+ highC.resize(1,5);
+
+ for (int i=0; i<5; i++) {
+ highC(0,i) = rho23*4.0*pi*(R1[i]+R1_bar[i]);
+
+ if (order>=1) {
+
+ highC(0,i) += rho23*(-0.5*a3*C_A_3*T1b[i]-0.5*a2*C_A_2*T1b_bar[i]+
+ (Lambda_Q-0.5*a3*C_A_3-0.5*a2*C_A_2)*(T1c[i]+T1c_bar[i])+
+ 0.5*a3*C_A_3*T3a[i]+0.5*a2*C_A_2*T3a_bar[i]+
+ (Lambda_Q-0.5*a3*C_A_3)*T3b[i]+(Lambda_Q-0.5*a2*C_A_2)*T3b_bar[i]+
+ 0.5*a2*C_A_2*T4a[i]+0.5*a3*C_A_3*T4a_bar[i]+
+ (Lambda_Q-0.5*a2*C_A_2)*T4b[i]+(Lambda_Q-0.5*a3*C_A_3)*T4b_bar[i]+
+ Lambda_Q*(T7[i]+T7_bar[i]));
+ }
+ }
+ }
+ break;
+ case QQBG:
+ {
+ // Finish Group Theory Factors:
+ Y_Q = 1.0/6.0;
+ Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
+ rho13 = Y_Q*sqrt(a1*a3);
+
+ highC.resize(1,5);
+
+ for (int i=0; i<5; i++) {
+ highC(0,i) = rho13*4.0*pi*(R1[i]+R1_bar[i]);
+
+ if (order>=1) {
+
+ highC(0,i) += rho13*(-0.5*a3*C_A_3*T1b[i]-0.5*a1*C_A_1*T1b_bar[i]+
+ (Lambda_Q-0.5*a3*C_A_3-0.5*a1*C_A_1)*(T1c[i]+T1c_bar[i])+
+ 0.5*a3*C_A_3*T3a[i]+0.5*a1*C_A_1*T3a_bar[i]+
+ (Lambda_Q-0.5*a3*C_A_3)*T3b[i]+(Lambda_Q-0.5*a1*C_A_1)*T3b_bar[i]+
+ 0.5*a1*C_A_1*T4a[i]+0.5*a3*C_A_3*T4a_bar[i]+
+ (Lambda_Q-0.5*a1*C_A_1)*T4b[i]+(Lambda_Q-0.5*a3*C_A_3)*T4b_bar[i]+
+ Lambda_Q*(T7[i]+T7_bar[i]));
+ }
+ }
+ }
+ break;
+ case UUBG: case DDBG:
+ {
+ // Finish Group Theory Factors:
+ if (process==UUBG) {
+ Y_Q = 2.0/3.0;
+ Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
+ rho13 = Y_Q*sqrt(a1*a3);
+ }
+ else if (process==DDBG) {
+ Y_Q = -1.0/3.0;
+ Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
+ rho13 = Y_Q*sqrt(a1*a3);
+ }
+
+ highC.resize(1,5);
+
+ for (int i=0; i<5; i++) {
+ highC(0,i) = rho13*4.0*pi*(R1[i]+R1_bar[i]);
+
+ if (order>=1) {
+
+ highC(0,i) += rho13*(-0.5*a3*C_A_3*T1b[i]-0.5*a1*C_A_1*T1b_bar[i]+
+ (Lambda_Q-0.5*a3*C_A_3-0.5*a1*C_A_1)*(T1c[i]+T1c_bar[i])+
+ 0.5*a3*C_A_3*T3a[i]+0.5*a1*C_A_1*T3a_bar[i]+
+ (Lambda_Q-0.5*a3*C_A_3)*T3b[i]+(Lambda_Q-0.5*a1*C_A_1)*T3b_bar[i]+
+ 0.5*a1*C_A_1*T4a[i]+0.5*a3*C_A_3*T4a_bar[i]+
+ (Lambda_Q-0.5*a1*C_A_1)*T4b[i]+(Lambda_Q-0.5*a3*C_A_3)*T4b_bar[i]+
+ Lambda_Q*(T7[i]+T7_bar[i]));
+ }
+ }
+ }
+ break;
+ case QQGG: case QtQtGG: case UUGG:
+ case tRtRGG: case DDGG:
+ {
+ // Finish Group Theory Factors:
+ if (process==QQGG || process==QtQtGG) {
+ Y_Q = 1.0/6.0;
+ Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
+ Lambda_Q_s = C_F_3*a3;
+ Lambda_Q_ew = C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
+ }
+ else if (process==UUGG || process==tRtRGG) {
+ Y_Q = 2.0/3.0;
+ Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
+ Lambda_Q_s = C_F_3*a3;
+ Lambda_Q_ew = Y_Q*Y_Q*C_F_1*a1;
+ }
+ else if (process==DDGG || process==tRtRGG) {
+ Y_Q = -1.0/3.0;
+ Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
+ Lambda_Q_s = C_F_3*a3;
+ Lambda_Q_ew = Y_Q*Y_Q*C_F_1*a1;
+ }
+
+ highC.resize(3,5);
+
+ for (int i=0; i<5; i++) {
+ highC(0,i) = G_Plus_SU3*(4.0*pi*a3)*(R1[i]+R1_bar[i]);
+ highC(1,i) = G_Plus_SU3_D*(4.0*pi*a3)*(R1[i]+R1_bar[i]);
+ highC(2,i) = G_f*4.0*pi*a3*(-0.5*R1[i]+0.5*R1_bar[i]-R2[i]);
+
+ if (order>=1) {
+ highC(0,i) += G_Plus_SU3*(a3*(Lambda_Q_ew)*(T1c[i]+T1c_bar[i])+
+ a3*(Lambda_Q_ew)*(T3b[i]+T3b_bar[i])+
+ a3*(Lambda_Q_ew)*(T4b[i]+T4b_bar[i])+
+ a3*Lambda_Q_ew*(T7[i]+T7_bar[i]));
+ highC(1,i) += G_Plus_SU3_D*(a3*(Lambda_Q_ew)*(T1c[i]+T1c_bar[i])+
+ a3*(Lambda_Q_ew)*(T3b[i]+T3b_bar[i])+
+ a3*(Lambda_Q_ew)*(T4b[i]+T4b_bar[i])+
+ a3*Lambda_Q_ew*(T7[i]+T7_bar[i]));
+ highC(2,i) += G_f*(a3*(-0.5*Lambda_Q_ew)*(T1c[i]-T1c_bar[i])-
+ 0.5*a3*(Lambda_Q_ew)*(T3b[i]-T3b_bar[i])-
+ 0.5*a3*(Lambda_Q_ew)*(T4b[i]-T4b_bar[i])+
+ a3*(Lambda_Q_ew)*T5b[i]-
+ 0.5*a3*Lambda_Q_ew*(T7[i]-T7_bar[i]));
+ if (includeAlphaS2) {
+ highC(0,i) += G_Plus_SU3*((-0.5*a3*a3*C_A_3)*(T1b[i]+T1b_bar[i])+
+ a3*(Lambda_Q_s-a3*C_A_3)*(T1c[i]+T1c_bar[i])+
+ 0.5*a3*a3*C_A_3*(T3a[i]+T3a_bar[i])+
+ a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T3b[i]+T3b_bar[i])+
+ 0.5*a3*a3*C_A_3*(T4a[i]+T4a_bar[i])+
+ a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T4b[i]+T4b_bar[i])+
+ a3*Lambda_Q_s*(T7[i]+T7_bar[i])) +
+ G_TT*(-a3*a3*(T1a[i]+T1a_bar[i])+a3*a3*(T1b[i]+T1b_bar[i])+
+ a3*a3*(T1c[i]+T1c_bar[i])+2.0*a3*a3*T1d[i]);
+ highC(1,i) += G_Plus_SU3_D*((-0.5*a3*a3*C_A_3)*(T1b[i]+T1b_bar[i])+
+ a3*(Lambda_Q_s-a3*C_A_3)*(T1c[i]+T1c_bar[i])+
+ 0.5*a3*a3*C_A_3*(T3a[i]+T3a_bar[i])+
+ a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T3b[i]+T3b_bar[i])+
+ 0.5*a3*a3*C_A_3*(T4a[i]+T4a_bar[i])+
+ a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T4b[i]+T4b_bar[i])+
+ a3*Lambda_Q_s*(T7[i]+T7_bar[i])) +
+ G_TT_3_D*(-a3*a3*(T1a[i]+T1a_bar[i])+a3*a3*(T1b[i]+T1b_bar[i])+
+ a3*a3*(T1c[i]+T1c_bar[i])+2.0*a3*a3*T1d[i]);
+ highC(2,i) += G_f*(0.25*a3*a3*C_A_3*(T1a[i]-T1a_bar[i])+
+ a3*(0.25*a3*C_A_3-0.5*Lambda_Q_s)*(T1c[i]-T1c_bar[i])+
+ 0.5*a3*a3*C_A_3*T2a[i]+a3*a3*T_CS_3*T2b[i]-
+ 0.5*a3*a3*C_A_3*T2c[i]+a3*a3*T_WF_3*T2d[i]-
+ 0.25*a3*a3*C_A_3*(T3a[i]-T3a_bar[i])-
+ 0.5*a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T3b[i]-T3b_bar[i])-
+ 0.25*a3*a3*C_A_3*(T4a[i]-T4a_bar[i])-
+ 0.5*a3*(Lambda_Q_s-0.5*a3*C_A_3)*(T4b[i]-T4b_bar[i])+
+ 0.5*a3*a3*C_A_3*T5a[i]+a3*(Lambda_Q_s-0.5*a3*C_A_3)*T5b[i]+
+ a3*a3*C_A_3*T6a[i]+a3*a3*C_A_3*T6b[i]+
+ a3*a3*T_CS_3*T6c[i]+a3*a3*T_WF_3*T6d[i]-
+ 0.5*a3*Lambda_Q_s*(T7[i]-T7_bar[i]));
+ }
+ }
+ }
+
+ if ( (process==QtQtGG||process==tRtRGG) && order>=1) {
+
+ if (process==tRtRGG) {
+ tRorQ = 2.0;
+ }
+ else {
+ tRorQ = 1.0;
+ }
+ highC(0,0) += tRorQ*(-1.0*(s*((s+t)*L_t - t*L_u)*y_t*y_t*a3)/(48.*pi*t*u*s));
+ highC(0,1) += tRorQ*((s*L_t*y_t*y_t*a3)/(24.*pi*t*s));
+ highC(0,2) += tRorQ*(-(s*s*y_t*y_t*a3)/((24.*pi*s*t+24.*pi*t*t)*s));
+ highC(0,3) += tRorQ*(-(s*s*y_t*y_t*a3)/((24.*pi*s*t+24.*pi*t*t)*s));
+ highC(1,0) += tRorQ*(-1.0*(s*((s+t)*L_t - t*L_u)*y_t*y_t*a3)/(16.*pi*t*u*s));
+ highC(1,1) += tRorQ*((s*L_t*y_t*y_t*a3)/(8.*pi*t*s));
+ highC(1,2) += tRorQ*(-(s*s*y_t*y_t*a3)/((8.*pi*s*t+8.*pi*t*t)*s));
+ highC(1,3) += tRorQ*(-(s*s*y_t*y_t*a3)/((8.*pi*s*t+8.*pi*t*t)*s));
+ highC(2,0) += tRorQ*((s*((s+t)*L_t + t*L_u)*y_t*y_t*a3)/(16.*pi*t*u*s));
+ highC(2,1) += tRorQ*(((2.*t-2.*t*L_s-s*L_t)*y_t*y_t*a3)/(8.*pi*t*s));
+ highC(2,2) += tRorQ*(-1.0*(s*(s+2.*t)*y_t*y_t*a3)/(8.*pi*t*u*s));
+ highC(2,3) += tRorQ*(-1.0*(s*(s+2.*t)*y_t*y_t*a3)/(8.*pi*t*u*s));
+ }
+ }
+ break;
+ default:
+ assert(false);
+ }
+ return highC;
+}
+
+boost::numeric::ublas::matrix<complex<InvEnergy2> >
+HighEnergyMatching::Spin0HighMatching(Energy highScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ EWProcess::Process process,
+ bool oneLoop, bool ) {
+ using Constants::pi;
+ unsigned int order = !oneLoop? 0 : 1;
+ // (If crossed graphs, swap s and t here)
+ Complex L_s = MinusLog(-s/(highScale*highScale));
+ Complex L_t = MinusLog(-t/(highScale*highScale));
+ Complex L_u = MinusLog(-u/(highScale*highScale));
+ Complex L_s2 = L_s*L_s;
+ Complex L_t2 = L_t*L_t;
+ Complex L_u2 = L_u*L_u;
+
+ // Tree-Level:
+ complex<InvEnergy2> S1 = 2.0/s;
+
+ // Topology T1:
+ complex<InvEnergy2> T1b = (-L_s2/(2.0*u)*(7.0*t/s+3.0)+2.0/u*L_t2+L_s*L_t*4.0/u*(t-u)/s+
+ L_s*2.0/s-4.0/s-pi*pi/(4.0*u)*(11.0+19.0*t/s));
+ complex<InvEnergy2> T1b_bar = -1.0*(-L_s2/(2.0*t)*(7.0*u/s+3.0)+2.0/t*L_u2+L_s*L_u*4.0/t*(u-t)/s+
+ L_s*2.0/s-4.0/s-pi*pi/(4.0*t)*(11.0+19.0*u/s));
+
+ // Topologies T2:
+ complex<InvEnergy2> T2a = 1.0/s*(-2.0*L_s2+8.0*L_s-16.0+pi*pi/3.0);
+ complex<InvEnergy2> T2b = 1.0/s*(0.5*L_s2+2.0*L_s-4.0-pi*pi/12.0);
+
+ // Topologies T5:
+ complex<InvEnergy2> T5a = 1.0/s*(-2.0*L_s2+6.0*L_s+pi*pi/3.0-16.0);
+ complex<InvEnergy2> T5b = 1.0/s*(2.0*L_s-4.0);
+
+ // Topologies T6:
+ complex<InvEnergy2> T6a = 1.0/s*(-19.0/6.0*L_s+58.0/9.0);
+ complex<InvEnergy2> T6b = 1.0/s*(-1.0/6.0*L_s+4.0/9.0);
+ complex<InvEnergy2> T6c = 1.0/s*(2.0/3.0*L_s-16.0/9.0);
+ complex<InvEnergy2> T6d = 1.0/s*(4.0/3.0*L_s-20.0/9.0);
+
+ // Group Theory Factors / SM parameters needed for matrix elements:
+ double a1 = ElectroWeakReweighter::coupling()->a1(highScale);
+ double a2 = ElectroWeakReweighter::coupling()->a2(highScale);
+ double a3 = ElectroWeakReweighter::coupling()->a3(highScale);
+ double y_t = ElectroWeakReweighter::coupling()->y_t(highScale);
+ double Y_phi = 1.0/2.0;
+ double C_F_3 = 4.0/3.0;
+ double C_F_2 = 3.0/4.0;
+ double C_F_1 = 1.0;
+ double n_g = 3.0;
+ double n_S = 1.0;
+ // Factors TBD after fermion helicity is specified:
+ double Y_Q(0.), Lambda_Q(0.), Lambda_phi(0.);
+ boost::numeric::ublas::matrix<complex<InvEnergy2> > highC(1,1);
+ switch (process) {
+
+ case QQPhiPhi: case LLPhiPhi:
+ // Finish Group Theory Factors:
+ if (process==QQPhiPhi) {
+ Y_Q = 1.0/6.0;
+ Lambda_Q = C_F_3*a3 + C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
+ Lambda_phi = C_F_2*a2+Y_phi*Y_phi*a1;
+ }
+ else if (process==LLPhiPhi) {
+ Y_Q = -1.0/2.0;
+ Lambda_Q = C_F_2*a2 + Y_Q*Y_Q*C_F_1*a1;
+ Lambda_phi = C_F_2*a2+Y_phi*Y_phi*a1;
+ }
+ highC.resize(2,1);
+ highC(0,0) = S1*(4.0*pi*a2);
+ highC(1,0) = S1*(4.0*pi*a1*Y_Q*Y_phi);
+ if (order>=1) {
+ highC(0,0) += T1b*(0.5*a2*a2+2.0*a1*a2*Y_Q*Y_phi) +
+ T1b_bar*(-0.5*a2*a2+2.0*a1*a2*Y_Q*Y_phi) +
+ T2a*(-a2*a2+Lambda_phi*a2) + T2b*a2*a2 +
+ T5a*(-a2*a2+Lambda_Q*a2) + T5b*a2*a2 +
+ T6a*2.0*a2*a2 + T6b*2.0*a2*a2 +
+ T6c*0.5*a2*a2*n_S + T6d*2.0*a2*a2*n_g;
+ highC(1,0) += T1b*(3.0/16.0*a2*a2+a1*a1*Y_Q*Y_Q*Y_phi*Y_phi) +
+ T1b_bar*(3.0/16.0*a2*a2+a1*a1*Y_Q*Y_Q*Y_phi*Y_phi) +
+ T2a*(Lambda_phi*a1*Y_Q*Y_phi) + T5a*(Lambda_Q*a1*Y_Q*Y_phi) +
+ T6c*(2.0*a1*a1*n_S*Y_Q*Y_phi*Y_phi*Y_phi) +
+ T6d*(10.0/3.0*a1*a1*n_g*Y_Q*Y_phi);
+ // Top Quark contributions:
+ highC(0,0) += -3.0*y_t*y_t*a2/(4.0*pi)/s*(2.0*L_s-4.0);
+ highC(1,0) += -3.0*y_t*y_t*a1/(4.0*pi)*(Y_Q*Y_phi)/s*(2.0*L_s-4.0);
+ }
+ break;
+ case UUPhiPhi: case DDPhiPhi: case EEPhiPhi:
+ // Finish Group Theory Factors:
+ if (process==UUPhiPhi) {
+ Y_Q = 2.0/3.0;
+ Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
+ Lambda_phi = C_F_2*a2 + Y_phi*Y_phi*a1;
+ }
+ else if (process==DDPhiPhi) {
+ Y_Q = -1.0/3.0;
+ Lambda_Q = C_F_3*a3 + Y_Q*Y_Q*C_F_1*a1;
+ Lambda_phi = C_F_2*a2 + Y_phi*Y_phi*a1;
+ }
+ else if (process==EEPhiPhi) {
+ Y_Q = -1.0;
+ Lambda_Q = Y_Q*Y_Q*C_F_1*a1;
+ Lambda_phi = C_F_2*a2 + Y_phi*Y_phi*a1;
+ }
+
+ highC.resize(1,1);
+ highC(0,0) = ZERO;
+
+ highC(0,0) = S1*(4.0*pi*a1*Y_Q*Y_phi);
+
+ if (order>=1) {
+ highC(0,0) += T1b*(a1*a1*Y_Q*Y_Q*Y_phi*Y_phi) +
+ T1b_bar*(a1*a1*Y_Q*Y_Q*Y_phi*Y_phi) +
+ T2a*(Lambda_phi*a1*Y_Q*Y_phi) + T5a*(Lambda_Q*a1*Y_Q*Y_phi) +
+ T6c*(2.0*a1*a1*n_S*Y_Q*Y_phi*Y_phi*Y_phi) +
+ T6d*(10.0/3.0*a1*a1*n_g*Y_Q*Y_phi);
+ // Top Quark Contribution:
+ highC(0,0) += -3.0*y_t*y_t*a1/(4.0*pi)*(Y_Q*Y_phi)/s*(2.0*L_s-4.0);
+ }
+ break;
+ default:
+ assert(false);
+ }
+ return highC;
+}
diff --git a/MatrixElement/EW/HighEnergyMatching.h b/MatrixElement/EW/HighEnergyMatching.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/HighEnergyMatching.h
@@ -0,0 +1,59 @@
+// -*- C++ -*-
+//
+// HighEnergyMatching.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_HighEnergyMatching_H
+#define HERWIG_HighEnergyMatching_H
+#include "ThePEG/Config/ThePEG.h"
+#include "ThePEG/Config/Unitsystem.h"
+#include "EWProcess.h"
+#include <boost/numeric/ublas/matrix.hpp>
+
+namespace Herwig {
+using namespace ThePEG;
+
+namespace HighEnergyMatching {
+
+ /**
+ * The high energy matching
+ */
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ highEnergyMatching(Energy highScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ Herwig::EWProcess::Process process,
+ bool oneLoop, bool includeAlphaS2);
+
+ /**
+ * Spin\f$\frac12\f$
+ */
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ SpinHalfMatching(Energy highScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ EWProcess::Process process,
+ bool oneLoop, bool includeAlphaS2);
+
+ /**
+ * Spin\f$1\f$
+ */
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ Spin1HighMatching(Energy highScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ EWProcess::Process process,
+ bool oneLoop, bool includeAlphaS2);
+ /**
+ * Spin\f$0\f$
+ */
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ Spin0HighMatching(Energy highScale,
+ Energy2 s, Energy2 t, Energy2 u,
+ EWProcess::Process process,
+ bool oneLoop, bool includeAlphaS2);
+}
+}
+
+#endif // HERWIG_HighEnergyMatching_H
diff --git a/MatrixElement/EW/Makefile.am b/MatrixElement/EW/Makefile.am
--- a/MatrixElement/EW/Makefile.am
+++ b/MatrixElement/EW/Makefile.am
@@ -1,6 +1,8 @@
pkglib_LTLIBRARIES = HwMEEW.la
-HwMEEW_la_SOURCES = GroupInvariants.h GroupInvariants.cc \
+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 \
EWCouplings.h EWCouplings.fh EWCouplings.cc
HwMEEW_la_LDFLAGS = $(AM_LDFLAGS) -module -version-info 1:0:0
diff --git a/MatrixElement/EW/SoftSudakov.cc b/MatrixElement/EW/SoftSudakov.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/SoftSudakov.cc
@@ -0,0 +1,92 @@
+// -*- 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;
+}
diff --git a/MatrixElement/EW/SoftSudakov.fh b/MatrixElement/EW/SoftSudakov.fh
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/SoftSudakov.fh
@@ -0,0 +1,18 @@
+// -*- C++ -*-
+//
+// This is the forward declaration of the SoftSudakov class.
+//
+#ifndef Herwig_SoftSudakov_FH
+#define Herwig_SoftSudakov_FH
+
+#include "ThePEG/Config/ThePEG.h"
+
+namespace Herwig {
+
+class SoftSudakov;
+
+ThePEG_DECLARE_POINTERS(Herwig::SoftSudakov,SoftSudakovPtr);
+
+}
+
+#endif
diff --git a/MatrixElement/EW/SoftSudakov.h b/MatrixElement/EW/SoftSudakov.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/EW/SoftSudakov.h
@@ -0,0 +1,162 @@
+// -*- 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 "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:
+
+ /**
+ * Evaluate bthe 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
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:07 PM (18 h, 57 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023576
Default Alt Text
(63 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment