Page MenuHomeHEPForge

No OneTemporary

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

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)

Event Timeline