Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8309909
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
214 KB
Subscribers
None
View Options
diff --git a/MatrixElement/EW/ElectroWeakMatching.cc b/MatrixElement/EW/ElectroWeakMatching.cc
--- a/MatrixElement/EW/ElectroWeakMatching.cc
+++ b/MatrixElement/EW/ElectroWeakMatching.cc
@@ -1,876 +1,915 @@
// -*- C++ -*-
//
// ElectroWeakMatching.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
//
#include "ElectroWeakMatching.h"
#include "ElectroWeakReweighter.h"
#include "GroupInvariants.h"
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/vector.hpp>
using namespace Herwig;
using namespace ElectroWeakMatching;
using namespace GroupInvariants;
using namespace EWProcess;
boost::numeric::ublas::matrix<Complex>
ElectroWeakMatching::electroWeakMatching(Energy mu,
Energy2 s, Energy2 t, Energy2 u,
Herwig::EWProcess::Process process,
- bool oneLoop) {
+ bool oneLoop,unsigned int iswap) {
static const Complex I(0,1.0);
using Constants::pi;
Complex T = getT(s,t);
Complex U = getU(s,u);
- // Z-Couplings
- double g_Lu = ElectroWeakReweighter::coupling()->g_Lu(mu);
- double g_Ld = ElectroWeakReweighter::coupling()->g_Ld(mu);
- double g_Le = ElectroWeakReweighter::coupling()->g_Le(mu);
- double g_Lnu = ElectroWeakReweighter::coupling()->g_Lnu(mu);
- double g_Ru = ElectroWeakReweighter::coupling()->g_Ru(mu);
- double g_Rd = ElectroWeakReweighter::coupling()->g_Rd(mu);
- double g_Re = ElectroWeakReweighter::coupling()->g_Re(mu);
- double g_W = ElectroWeakReweighter::coupling()->g_W(mu);
- double g_phiPlus = ElectroWeakReweighter::coupling()->g_phiPlus(mu);
- // Weinberg Angle:
- double cos2 = ElectroWeakReweighter::coupling()->Cos2thW(mu);
- double sin2 = 1.0-cos2;
- double cos = sqrt(cos2);
- double sin = sqrt(sin2);
+ // Z-Couplings
+ double g_Lu = ElectroWeakReweighter::coupling()->g_Lu(mu);
+ double g_Ld = ElectroWeakReweighter::coupling()->g_Ld(mu);
+ double g_Le = ElectroWeakReweighter::coupling()->g_Le(mu);
+ double g_Lnu = ElectroWeakReweighter::coupling()->g_Lnu(mu);
+ double g_Ru = ElectroWeakReweighter::coupling()->g_Ru(mu);
+ double g_Rd = ElectroWeakReweighter::coupling()->g_Rd(mu);
+ double g_Re = ElectroWeakReweighter::coupling()->g_Re(mu);
+ double g_W = ElectroWeakReweighter::coupling()->g_W(mu);
+ double g_phiPlus = ElectroWeakReweighter::coupling()->g_phiPlus(mu);
+ // Weinberg Angle:
+ double cos2 = ElectroWeakReweighter::coupling()->Cos2thW(mu);
+ double sin2 = 1.0-cos2;
+ double cos = sqrt(cos2);
+ double sin = sqrt(sin2);
- boost::numeric::ublas::matrix<Complex> R0,G2,Dw,Dz;
-
- switch (process) {
- case QQQQ:
- case QQQQiden:
- case QtQtQQ:
- {
- unsigned int numGauge = 4, numBrokenGauge = 12;
- R0=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2=boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
- R0(0,0) = R0(3,0) = 0.25;
- R0(1,0) = R0(2,0) = -0.25;
- R0(4,0) = R0(5,0) = 0.5;
- R0(6,3) = R0(7,3) = R0(8,3) = R0(9,3) = 1.0;
- R0(6,2) = R0(9,2) = 0.25;
- R0(7,2) = R0(8,2) = -0.25;
- R0(10,2) = R0(11,2) = 0.5;
- if (oneLoop) {
- double g11 = g_Lu;
- double g12 = g_Ld;
- double g21 = g_Lu;
- double g22 = g_Ld;
- for(unsigned int ix=0;ix<numBrokenGauge;++ix) {
- Dw(ix,ix) = 0.5*I*pi;
- }
- Complex w1 = -0.5*(T-U);
- Complex w2 = -0.5*(T+U);
- for(unsigned int ix=0;ix<numBrokenGauge;ix+=6) {
- Dw(ix+0,ix+0) += w1;
- Dw(ix+3,ix+3) += w1;
- Dw(ix+1,ix+1) +=-w1;
- Dw(ix+2,ix+2) +=-w1;
- Dw(ix+4,ix+4) += w2;
- Dw(ix+5,ix+5) += w2;
- }
- Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
- Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
- Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
- Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
- - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
- for(unsigned int ix=0;ix<numBrokenGauge;ix+=6) {
- Dz(ix+0,ix+0) = z1;
- Dz(ix+1,ix+1) = z2;
- Dz(ix+2,ix+2) = z3;
- Dz(ix+3,ix+3) = z4;
- Dz(ix+4,ix+4) = z5;
- Dz(ix+5,ix+5) = z5;
- }
- boost::numeric::ublas::matrix<Complex> gam2 = Gamma2(U,T);
- G2(0,0) += gam2(0,0);
- G2(0,1) += gam2(0,1);
- G2(1,0) += gam2(1,0);
- G2(1,1) += gam2(1,1);
- G2(2,2) += gam2(0,0);
- G2(2,3) += gam2(0,1);
- G2(3,2) += gam2(1,0);
- G2(3,3) += gam2(1,1);
- }
- }
- break;
- case QQUU:
- case QtQtUU:
- case QQtRtR:
- {
- unsigned int numGauge = 2, numBrokenGauge = 4;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,0) = R0(2,1) = R0(3,1) = 1.0;
- if (oneLoop) {
- double g11 = g_Lu;
- double g12 = g_Ld;
- //double g21 = g_Ru;
- double g22 = g_Ru;
-
- Complex w1 = 0.25*I*pi;
- for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;
-
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- for(unsigned int ix=0;ix<numBrokenGauge;ix+=2) {
- Dz(ix+0,ix+0) = z1;
- Dz(ix+1,ix+1) = z2;
- }
- G2 = Gamma2Singlet();
- }
- }
- break;
- case QQDD:
- case QtQtDD:
- {
- unsigned int numGauge = 2, numBrokenGauge = 4;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,0) = R0(2,1) = R0(3,1) = 1.0;
- if (oneLoop) {
- double g11 = g_Lu;
- double g12 = g_Ld;
- //double g21 = g_Rd;
- double g22 = g_Rd;
-
- Complex w1 = 0.25*I*pi;
- for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;
-
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- for(unsigned int ix=0;ix<numBrokenGauge;ix+=2) {
- Dz(ix+0,ix+0) = z1;
- Dz(ix+1,ix+1) = z2;
- }
- G2 = Gamma2Singlet();
- }
- }
- break;
- case QQLL:
- {
- unsigned int numGauge = 2, numBrokenGauge = 6;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
- R0(0,0) = R0(3,0) = 0.25;
- R0(1,0) = R0(2,0) = -0.25;
- R0(4,0) = R0(5,0) = 0.5;
- if (oneLoop) {
- double g11 = g_Lu;
- double g12 = g_Ld;
- double g21 = g_Lnu;
- double g22 = g_Le;
-
- for (unsigned int i=0; i<6; ++i) {
- Dw(i,i) = 0.5*I*pi;
- }
- Complex w1 = (-1.0/2.0)*(T-U);
- Complex w2 = (-1.0/2.0)*(T+U);
- Dw(0,0) += w1;
- Dw(3,3) += w1;
- Dw(1,1) += -1.0*w1;
- Dw(2,2) += -1.0*w1;
- Dw(4,4) += w2;
- Dw(5,5) += w2;
-
- Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
- Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
- Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
- Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
- - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
- Dz(0,0) = z1;
- Dz(1,1) = z2;
- Dz(2,2) = z3;
- Dz(3,3) = z4;
- Dz(4,4) = Dz(5,5) = z5;
-
- G2 = Gamma2(U,T);
- }
- }
- break;
- case QQEE:
- {
- unsigned int numGauge = 1, numBrokenGauge = 2;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,0) = 1.0;
- if (oneLoop) {
- double g11 = g_Lu;
- double g12 = g_Ld;
- //double g21 = g_Re;
- double g22 = g_Re;
-
- Complex w1 = 0.25*I*pi;
- Dw(0,0) = Dw(1,1) = w1;
-
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- Dz(0,0) = z1;
- Dz(1,1) = z2;
-
- G2(0,0) = Gamma2Singlet()(0,0);
- }
- }
- break;
- case UUUU:
- case UUUUiden:
- case tRtRUU:
- {
- unsigned int numGauge = 2, numBrokenGauge = 2;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,1) = 1.0;
- if (oneLoop) {
- double g11 = g_Ru;
- //double g12 = g_Ru;
- //double g21 = g_Ru;
- double g22 = g_Ru;
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Dz(0,0) = Dz(1,1) = z1;
- }
- }
- break;
- case UUDD:
- case tRtRDD:
- {
- unsigned int numGauge = 2, numBrokenGauge = 2;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,1) = 1.0;
- if (oneLoop) {
- double g11 = g_Ru;
- //double g12 = g_Ru;
- //double g21 = g_Rd;
- double g22 = g_Rd;
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Dz(0,0) = Dz(1,1) = z1;
- }
- }
- break;
- case UULL:
- {
- unsigned int numGauge = 1, numBrokenGauge = 2;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,0) = 1.0;
- if (oneLoop) {
- double g11 = g_Lnu;
- double g12 = g_Le;
- //double g21 = g_Ru;
- double g22 = g_Ru;
-
- Complex w1 = 0.25*I*pi;
- Dw(0,0) = Dw(1,1) = w1;
-
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- Dz(0,0) = z1;
- Dz(1,1) = z2;
-
- G2(0,0) = Gamma2Singlet()(0,0);
- }
- }
- break;
- case UUEE:
- {
- unsigned int numGauge = 1, numBrokenGauge = 1;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = 1.0;
- if (oneLoop) {
- double g11 = g_Ru;
- //double g12 = g_Ru;
- //double g21 = g_Re;
- double g22 = g_Re;
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Dz(0,0) = z1;
- }
- }
- break;
- case DDDD:
- case DDDDiden:
- {
- unsigned int numGauge = 2, numBrokenGauge = 2;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,1) = 1.0;
- if (oneLoop) {
- double g11 = g_Rd;
- //double g12 = g_Rd;
- //double g21 = g_Rd;
- double g22 = g_Rd;
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Dz(0,0) = Dz(1,1) = z1;
- }
- }
- break;
- case DDLL:
- {
- unsigned int numGauge = 1, numBrokenGauge = 2;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,0) = 1.0;
- if (oneLoop) {
- double g11 = g_Lnu;
- double g12 = g_Le;
- //double g21 = g_Rd;
- double g22 = g_Rd;
-
- Complex w1 = 0.25*I*pi;
- Dw(0,0) = Dw(1,1) = w1;
-
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- Dz(0,0) = z1;
- Dz(1,1) = z2;
-
- G2(0,0) = Gamma2Singlet()(0,0);
- }
- }
- break;
- case DDEE:
- {
- unsigned int numGauge = 1, numBrokenGauge = 1;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0 *= 0.0; Dw = Dz *= 0.0;
- R0(0,0) = 1.0;
- if (oneLoop) {
- double g11 = g_Rd;
- //double g12 = g_Rd;
- //double g21 = g_Re;
- double g22 = g_Re;
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Dz(0,0) = z1;
- }
- }
- break;
- case LLLL:
- case LLLLiden:
- {
- unsigned int numGauge = 2, numBrokenGauge = 6;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
- R0(0,0) = R0(3,0) = 0.25;
- R0(1,0) = R0(2,0) = -0.25;
- R0(4,0) = R0(5,0) = 0.5;
- if (oneLoop) {
- double g11 = g_Lnu;
- double g12 = g_Le;
- double g21 = g_Lnu;
- double g22 = g_Le;
-
- for (int i=0; i<6; i++) {
- Dw(i,i) = 0.5*I*pi;
- }
- Complex w1 = (-1.0/2.0)*(T-U);
- Complex w2 = (-1.0/2.0)*(T+U);
- Dw(0,0) += w1;
- Dw(3,3) += w1;
- Dw(1,1) += -1.0*w1;
- Dw(2,2) += -1.0*w1;
- Dw(4,4) += w2;
- Dw(5,5) += w2;
-
- Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
- Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
- Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
- Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
- - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
- Dz(0,0) = z1;
- Dz(1,1) = z2;
- Dz(2,2) = z3;
- Dz(3,3) = z4;
- Dz(4,4) = Dz(5,5) = z5;
-
- G2 = Gamma2(U,T);
- }
- }
- break;
- case LLEE:
- {
- unsigned int numGauge = 1, numBrokenGauge = 2;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,0) = 1.0;
- if (oneLoop) {
- double g11 = g_Lnu;
- double g12 = g_Le;
- //double g21 = g_Re;
- double g22 = g_Re;
-
- Complex w1 = 0.25*I*pi;
- Dw(0,0) = Dw(1,1) = w1;
-
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
- Dz(0,0) = z1;
- Dz(1,1) = z2;
-
- G2(0,0) = Gamma2Singlet()(0,0);
- }
- }
- break;
- case EEEE:
- case EEEEiden:
- {
- unsigned int numGauge = 1, numBrokenGauge = 1;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = 1.0;
- if (oneLoop) {
- double g11 = g_Re;
- //double g12 = g_Re;
- //double g21 = g_Re;
- double g22 = g_Re;
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
- Dz(0,0) = z1;
- }
- }
- break;
- case QQWW:
- case LLWW:
- {
- unsigned int numGauge = 5, numBrokenGauge = 20;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = 1.0; R0(0,1) = 0.5;
- R0(1,0) = 1.0; R0(1,1) = -0.5;
- R0(2,0) = cos2; R0(2,2) = -0.5*sin*cos; R0(2,3) = -0.5*sin*cos; R0(2,4) = sin2;
- R0(3,0) = sin*cos; R0(3,2) = 0.5*cos2; R0(3,3) = -0.5*sin2; R0(3,4) = -sin*cos;
- R0(4,0) = sin*cos; R0(4,2) = -0.5*sin2; R0(4,3) = 0.5*cos2; R0(4,4) = -sin*cos;
- R0(5,0) = sin2; R0(5,2) = 0.5*sin*cos; R0(5,3) = 0.5*sin*cos; R0(5,4) = cos2;
- R0(6,0) = 1.0; R0(6,1) = -0.5;
- R0(7,0) = 1.0; R0(7,1) = 0.5;
- R0(8,0) = cos2; R0(8,2) = 0.5*sin*cos; R0(8,3) = 0.5*sin*cos; R0(8,4) = sin2;
- R0(9,0) = sin*cos; R0(9,2) = -0.5*cos2; R0(9,3) = 0.5*sin2; R0(9,4) = -sin*cos;
- R0(10,0) = sin*cos; R0(10,2) = 0.5*sin2; R0(10,3) = -0.5*cos2; R0(10,4) = -sin*cos;
- R0(11,0) = sin2; R0(11,2) = -0.5*sin*cos; R0(11,3) = -0.5*sin*cos; R0(11,4) = cos2;
- R0(12,1) = -cos/sqrt(2.0); R0(12,3) = -sin/sqrt(2.0);
- R0(13,1) = -sin/sqrt(2.0); R0(13,3) = cos/sqrt(2.0);
- R0(14,1) = cos/sqrt(2.0); R0(14,2) = -sin/sqrt(2.0);
- R0(15,1) = sin/sqrt(2.0); R0(15,2) = cos/sqrt(2.0);
- R0(16,1) = -cos/sqrt(2.0); R0(16,2) = -sin/sqrt(2.0);
- R0(17,1) = -sin/sqrt(2.0); R0(17,2) = cos/sqrt(2.0);
- R0(18,1) = cos/sqrt(2.0); R0(18,3) = -sin/sqrt(2.0);
- R0(19,1) = sin/sqrt(2.0); R0(19,3) = cos/sqrt(2.0);
- if (oneLoop) {
- double gW = g_W;
- double g1(0.),g2(0.);
- if (process==QQWW) {
- g1 = g_Lu;
- g2 = g_Ld;
- }
- else if (process==LLWW) {
- g1 = g_Lnu;
- g2 = g_Le;
- }
-
- Complex w1 = T-U+5.0/4.0*I*pi;
- Complex w2 = -T+U+5.0/4.0*I*pi;
- Complex w3 = -0.5*(T+U) + 3.0/4.0*I*pi;
- Complex w4 = 0.25*I*pi;
- Dw(0,0) = Dw(7,7) = w1;
- Dw(1,1) = Dw(6,6) = w2;
- for (unsigned int i=12; i<20; i++) {
- Dw(i,i) = w3;
- }
- Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w4;
- Dw(8,8) = Dw(9,9) = Dw(10,10) = Dw(11,11) = w4;
-
- Complex z1 = 2.0*g1*gW*(U-T) - I*pi*(g1*g1+gW*gW);
- Complex z2 = 2.0*g1*gW*(T-U) - I*pi*(g1*g1+gW*gW);
- Complex z3 = 2.0*g2*gW*(U-T) - I*pi*(g2*g2+gW*gW);
- Complex z4 = 2.0*g2*gW*(T-U) - I*pi*(g2*g2+gW*gW);
- Complex z5 = -(g2*gW)*T + (g1*gW)*U - I*pi*(g1*g2+g1*gW-g2*gW);
- Complex z6 = (g1*gW)*T - (g2*gW)*U - I*pi*(g1*g2+g1*gW-g2*gW);
- Complex z7 = -I*pi*g1*g1;
- Complex z8 = -I*pi*g2*g2;
- Dz(0,0) = z1;
- Dz(1,1) = z2;
- Dz(2,2) = Dz(3,3) = Dz(4,4) = Dz(5,5) = z7;
- Dz(6,6) = z3;
- Dz(7,7) = z4;
- Dz(8,8) = Dz(9,9) = Dz(10,10) = Dz(11,11) = z8;
- Dz(12,12) = Dz(13,13) = Dz(16,16) = Dz(17,17) = z5;
- Dz(14,14) = Dz(15,15) = Dz(18,18) = Dz(19,19) = z6;
-
- G2 = Gamma2w(U,T);
- }
- }
- break;
- case QQPhiPhi:
- case LLPhiPhi:
- {
- unsigned int numGauge = 2, numBrokenGauge = 14;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = 0.25; R0(0,1) = 1.0;
- R0(1,0) = -1.0/8.0; R0(1,1) = 0.5;
- R0(2,0) = I/8.0; R0(2,1) = -I/2.0;
- R0(3,0) = -I/8.0; R0(3,1) = I/2.0;
- R0(4,0) = -1.0/8.0; R0(4,1) = 1.0/2.0;
- R0(5,0) = -1.0/4.0; R0(5,1) = 1.0;
- R0(6,0) = 1.0/8.0; R0(6,1) = 1.0/2.0;
- R0(7,0) = -I/8.0; R0(7,1) = -I/2.0;
- R0(8,0) = I/8.0; R0(8,1) = I/2.0;
- R0(9,0) = 1.0/8.0; R0(9,1) = 1.0/2.0;
- R0(10,0) = -1.0/(2.0*sqrt(2.0));
- R0(11,0) = I/(2.0*sqrt(2.0));
- R0(12,0) = -1.0/(2.0*sqrt(2.0));
- R0(13,0) = -I/(2.0*sqrt(2.0));
-
- if (oneLoop) {
- double g1(0.),g2(0.);
- if (process==QQPhiPhi) {
- g1 = g_Lu;
- g2 = g_Ld;
- }
- else if (process==LLPhiPhi) {
- g1 = g_Lnu;
- g2 = g_Le;
- }
- else
- assert(false);
- double g3 = g_phiPlus;
-
- Complex w0 = 0.25*I*pi;
- Complex w1 = 0.5*(T-U) + 0.5*I*pi;
- Complex w2 = -0.5*(T-U) + 0.5*I*pi;
- Complex w3 = 0.25*I*(T-U);
- Complex w4 = -0.25*(T+U) + 0.25*I*pi;
- Dw(0,0) = w2;
- Dw(1,1) = w0; Dw(1,2) = w3; Dw(1,3) = -w3; Dw(1,4) = w0;
- Dw(2,1) = -w3; Dw(2,2) = w0; Dw(2,3) = -w0; Dw(2,4) = -w3;
- Dw(3,1) = w3; Dw(3,2) = -w0; Dw(3,3) = w0; Dw(3,4) = w3;
- Dw(4,1) = w0; Dw(4,2) = w3; Dw(4,3) = -w3; Dw(4,4) = w0;
- Dw(5,5) = w1;
- Dw(6,6) = w0; Dw(6,7) = -w3; Dw(6,8) = w3; Dw(6,9) = w0;
- Dw(7,6) = w3; Dw(7,7) = w0; Dw(7,8) = -w0; Dw(7,9) = w3;
- Dw(8,6) = -w3; Dw(8,7) = -w0; Dw(8,8) = w0; Dw(8,9) = -w3;
- Dw(9,6) = w0; Dw(9,7) = -w3; Dw(9,8) = w3; Dw(9,9) = w0;
- Dw(10,10) = w4; Dw(10,11) = I*w4;
- Dw(11,10) = -I*w4; Dw(11,11) = w4;
- Dw(12,12) = w4; Dw(12,13) = -I*w4;
- Dw(13,12) = I*w4; Dw(13,13) = w4;
-
- Complex z1 = 2.0*g3*g1*(T-U) - I*pi*(g3*g3+g1*g1);
- Complex z2 = 2.0*g3*g2*(T-U) - I*pi*(g3*g3+g2*g2);
- Complex z3 = -I*pi*g1*g1;
- Complex z4 = 0.5*I*g1*(T-U);
- Complex z5 = 0.25*I*pi;
- Complex z6 = -I*pi*g2*g2;
- Complex z7 = 0.5*I*g2*(T-U);
- Complex z8 = g3*g1*T-g3*g2*U-I*pi*(g1*g2-g2*g3+g1*g3);
- Complex z9 = 0.5*I*g2*T-0.5*I*g1*U+pi/2.0*g2-pi/2.0*g1+pi/2.0*g3;
- Dz(0,0) = z1;
- Dz(1,1) = z3; Dz(1,2) = -z4; Dz(1,3) = z4; Dz(1,4) = -z5;
- Dz(2,1) = z4; Dz(2,2) = z3; Dz(2,3) = z5; Dz(2,4) = z4;
- Dz(3,1) = -z4; Dz(3,2) = z5; Dz(3,3) = z3; Dz(3,4) = -z4;
- Dz(4,1) = -z5; Dz(4,2) = -z4; Dz(4,3) = z4; Dz(4,4) = z3;
- Dz(5,5) = z2;
- Dz(6,6) = z6; Dz(6,7) = -z7; Dz(6,8) = z7; Dz(6,9) = -z5;
- Dz(7,6) = z7; Dz(7,7) = z6; Dz(7,8) = z5; Dz(7,9) = z7;
- Dz(8,6) = -z7; Dz(8,7) = z5; Dz(8,8) = z6; Dz(8,9) = -z7;
- Dz(9,6) = -z5; Dz(9,7) = -z7; Dz(9,8) = z7; Dz(9,9) = z6;
- Dz(10,10) = z8; Dz(10,11) = -z9;
- Dz(11,10) = z9; Dz(11,11) = z8;
- Dz(12,12) = z8; Dz(12,13) = z9;
- Dz(13,12) = -z9; Dz(13,13) = z8;
-
- G2 = Gamma2(U,T);
- }
- }
- break;
- case QQWG:
- {
- unsigned int numGauge = 1, numBrokenGauge = 6;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = 1.0/sqrt(2);
- R0(1,0) = 1.0/sqrt(2);
- R0(2,0) = cos/2.0;
- R0(3,0) = sin/2.0;
- R0(4,0) = -cos/2.0;
- R0(5,0) = -sin/2.0;
- if (oneLoop) {
- double g1 = g_Lu;
- double g2 = g_Ld;
- double gW = g_W;
-
- Complex w1 = -0.5*(T+U) + 0.75*I*pi;
- Complex w2 = 0.25*I*pi;
- Dw(0,0) = Dw(1,1) = w1;
- Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;
-
- Complex z1 = gW*g1*T - gW*g2*U - I*pi*(g1*g2+g1*gW-g2*gW);
- Complex z2 = gW*g1*U - gW*g2*T - I*pi*(g2*g1+g1*gW-g2*gW);
- Complex z3 = -I*pi*g1*g1;
- Complex z4 = -I*pi*g2*g2;
- Dz(0,0) = z1;
- Dz(1,1) = z2;
- Dz(2,2) = z3;
- Dz(3,3) = z3;
- Dz(4,4) = z4;
- Dz(5,5) = z4;
-
- G2(0,0) = -7.0/4.0*I*pi + (U+T);
- }
- }
- break;
- case QQBG:
- {
- unsigned int numGauge = 1, numBrokenGauge = 4;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = -sin;
- R0(1,0) = cos;
- R0(2,0) = -sin;
- R0(3,0) = cos;
- if (oneLoop) {
- double g1 = g_Lu;
- double g2 = g_Ld;
- Complex w2 = 0.25*I*pi;
- Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = w2;
- Complex z3 = -I*pi*g1*g1;
- Complex z4 = -I*pi*g2*g2;
- Dz(0,0) = z3;
- Dz(1,1) = z3;
- Dz(2,2) = z4;
- Dz(3,3) = z4;
- G2(0,0) = Gamma2Singlet()(0,0);
+ boost::numeric::ublas::matrix<Complex> R0,G2,Dw,Dz;
+
+ switch (process) {
+ case QQQQ:
+ case QQQQiden:
+ case QtQtQQ:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 4, numBrokenGauge = 12;
+ R0=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2=boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz=boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
+ R0(0,0) = R0(3,0) = 0.25;
+ R0(1,0) = R0(2,0) = -0.25;
+ R0(4,0) = R0(5,0) = 0.5;
+ R0(6,3) = R0(7,3) = R0(8,3) = R0(9,3) = 1.0;
+ R0(6,2) = R0(9,2) = 0.25;
+ R0(7,2) = R0(8,2) = -0.25;
+ R0(10,2) = R0(11,2) = 0.5;
+ if (oneLoop) {
+ double g11 = g_Lu;
+ double g12 = g_Ld;
+ double g21 = g_Lu;
+ double g22 = g_Ld;
+ for(unsigned int ix=0;ix<numBrokenGauge;++ix) {
+ Dw(ix,ix) = 0.5*I*pi;
}
+ Complex w1 = -0.5*(T-U);
+ Complex w2 = -0.5*(T+U);
+ for(unsigned int ix=0;ix<numBrokenGauge;ix+=6) {
+ Dw(ix+0,ix+0) += w1;
+ Dw(ix+3,ix+3) += w1;
+ Dw(ix+1,ix+1) +=-w1;
+ Dw(ix+2,ix+2) +=-w1;
+ Dw(ix+4,ix+4) += w2;
+ Dw(ix+5,ix+5) += w2;
+ }
+ Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
+ Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
+ Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
+ Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
+ - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
+ for(unsigned int ix=0;ix<numBrokenGauge;ix+=6) {
+ Dz(ix+0,ix+0) = z1;
+ Dz(ix+1,ix+1) = z2;
+ Dz(ix+2,ix+2) = z3;
+ Dz(ix+3,ix+3) = z4;
+ Dz(ix+4,ix+4) = z5;
+ Dz(ix+5,ix+5) = z5;
+ }
+ boost::numeric::ublas::matrix<Complex> gam2 = Gamma2(U,T);
+ G2(0,0) += gam2(0,0);
+ G2(0,1) += gam2(0,1);
+ G2(1,0) += gam2(1,0);
+ G2(1,1) += gam2(1,1);
+ G2(2,2) += gam2(0,0);
+ G2(2,3) += gam2(0,1);
+ G2(3,2) += gam2(1,0);
+ G2(3,3) += gam2(1,1);
}
- break;
- case QQGG:
- case QtQtGG:
- {
- unsigned int numGauge = 3, numBrokenGauge = 6;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(3,0) = 1.0;
- R0(1,1) = R0(4,1) = 1.0;
- R0(2,2) = R0(5,2) = 1.0;
- if (oneLoop) {
- double g1 = g_Lu;
- double g2 = g_Ld;
- Complex w2 = 0.25*I*pi;
- Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;
- Complex z3 = -I*pi*g1*g1;
- Complex z4 = -I*pi*g2*g2;
- Dz(0,0) = Dz(1,1) = Dz(2,2) = z3;
- Dz(3,3) = Dz(4,4) = Dz(5,5) = z4;
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 *= 0.0;
- G2(0,0) = G2(1,1) = G2(2,2) = Gamma2Singlet()(0,0);
- }
- }
- break;
- case UUBB:
- case DDBB:
- case EEBB:
- {
- unsigned int numGauge = 1, numBrokenGauge = 4;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = sin2;
- R0(1,0) = -sin*cos;
- R0(2,0) = -sin*cos;
- R0(3,0) = cos2;
- if (oneLoop) {
- double g1(0.);
- if (process==UUBB) {
- g1 = g_Ru;
- }
- else if (process==DDBB) {
- g1 = g_Rd;
- }
- else if (process==EEBB) {
- g1 = g_Re;
- }
- else
- assert(false);
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = -I*pi*g1*g1;
- Dz(0,0) = Dz(1,1) = Dz(2,2) = Dz(3,3) = z1;
- }
- }
- break;
- case UUPhiPhi:
- case DDPhiPhi:
- case EEPhiPhi:
- {
- unsigned int numGauge = 1, numBrokenGauge = 5;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = 1.0;
- R0(1,0) = 0.5;
- R0(2,0) = -0.5*I;
- R0(3,0) = 0.5*I;
- R0(4,0) = 0.5;
- if (oneLoop) {
- double g1(0.);
- if (process==UUPhiPhi) {
- g1 = g_Ru;
- }
- else if (process==DDPhiPhi) {
- g1 = g_Rd;
- }
- else if (process==EEPhiPhi) {
- g1 = g_Re;
- }
- double g3 = g_phiPlus;
- Dw(0,0) = Dw(1,4) = Dw(4,1) = 0.25*I*pi;
- Dw(2,3) = Dw(3,2) = -0.25*I*pi;
- Complex z1 = 2.0*g3*g1*(T-U) - I*pi*(g3*g3+g1*g1);
- Complex z2 = 0.5*I*g1*g1;
- Complex z3 = -I*pi*g1*g1;
- Complex z4 = 0.25*I*pi;
- Dz(0,0) = z1;
- Dz(1,1) = z3; Dz(1,2) = -z2; Dz(1,3) = z2; Dz(1,4) = -z4;
- Dz(2,1) = z2; Dz(2,2) = z3; Dz(2,3) = z4; Dz(2,4) = z2;
- Dz(3,1) = -z2; Dz(3,2) = z4; Dz(3,3) = z3; Dz(3,4) = -z2;
- Dz(4,1) = -z4; Dz(4,2) = -z2; Dz(4,3) = z2; Dz(4,4) = z3;
- G2(0,0) = Gamma2Singlet()(0,0);
+ }
+ break;
+ case QQUU:
+ case QtQtUU:
+ case QQtRtR:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 2, numBrokenGauge = 4;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,0) = R0(2,1) = R0(3,1) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Lu;
+ double g12 = g_Ld;
+ //double g21 = g_Ru;
+ double g22 = g_Ru;
+
+ Complex w1 = 0.25*I*pi;
+ for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;
+
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ for(unsigned int ix=0;ix<numBrokenGauge;ix+=2) {
+ Dz(ix+0,ix+0) = z1;
+ Dz(ix+1,ix+1) = z2;
}
- }
- break;
- case UUBG:
- case DDBG:
- {
- unsigned int numGauge = 1, numBrokenGauge = 2;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = -sin;
- R0(1,0) = cos;
- if (oneLoop) {
- double g1(0.);
- if (process==UUBG) {
- g1 = g_Ru;
- }
- else if (process==DDBG) {
- g1 = g_Rd;
- }
- else
- assert(false);
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = -I*pi*g1*g1;
- Dz(0,0) = Dz(1,1) = z1;
- }
- }
- break;
- case UUGG:
- case tRtRGG:
- case DDGG:
- {
- unsigned int numGauge = 3, numBrokenGauge = 3;
- R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- R0(0,0) = R0(1,1) = R0(2,2) = 1.0;
- if (oneLoop) {
- double g1(0.);
- if ((process==UUGG)||(process==tRtRGG)) {
- g1 = g_Ru;
- }
- else if (process==DDGG) {
- g1 = g_Rd;
- }
- else
- assert(false);
- // There is no Dw contribution for two SU(2) singlets.
- Complex z1 = -I*pi*g1*g1;
- Dz(0,0) = Dz(1,1) = Dz(2,2) = z1;
- }
- }
- break;
- default:
- assert(false);
- }
-
- double aW = ElectroWeakReweighter::coupling()->aW(mu);
- double aZ = ElectroWeakReweighter::coupling()->aZ(mu);
- Energy mZ = ElectroWeakReweighter::coupling()->mZ();
- Energy mW = ElectroWeakReweighter::coupling()->mW();
-
- if (!oneLoop) {
- return R0;
- }
- boost::numeric::ublas::matrix<Complex> output(R0);
- boost::numeric::ublas::matrix<Complex> temp(R0.size1(),R0.size2());
- boost::numeric::ublas::axpy_prod(R0,G2,temp);
- output+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
- boost::numeric::ublas::axpy_prod(Dw,R0,temp);
- output+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
- boost::numeric::ublas::axpy_prod(Dz,R0,temp);
- output+=aZ/(4.0*pi)*4.0*log(mZ/mu)*temp;
- return output;
+ G2 = Gamma2Singlet();
+ }
+ }
+ break;
+ case QQDD:
+ case QtQtDD:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 2, numBrokenGauge = 4;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,0) = R0(2,1) = R0(3,1) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Lu;
+ double g12 = g_Ld;
+ //double g21 = g_Rd;
+ double g22 = g_Rd;
+
+ Complex w1 = 0.25*I*pi;
+ for(unsigned int ix=0;ix<numBrokenGauge;++ix) Dw(ix,ix) = w1;
+
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ for(unsigned int ix=0;ix<numBrokenGauge;ix+=2) {
+ Dz(ix+0,ix+0) = z1;
+ Dz(ix+1,ix+1) = z2;
+ }
+ G2 = Gamma2Singlet();
+ }
+ }
+ break;
+ case QQLL:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 2, numBrokenGauge = 6;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
+ R0(0,0) = R0(3,0) = 0.25;
+ R0(1,0) = R0(2,0) = -0.25;
+ R0(4,0) = R0(5,0) = 0.5;
+ if (oneLoop) {
+ double g11 = g_Lu;
+ double g12 = g_Ld;
+ double g21 = g_Lnu;
+ double g22 = g_Le;
+
+ for (unsigned int i=0; i<6; ++i) {
+ Dw(i,i) = 0.5*I*pi;
+ }
+ Complex w1 = (-1.0/2.0)*(T-U);
+ Complex w2 = (-1.0/2.0)*(T+U);
+ Dw(0,0) += w1;
+ Dw(3,3) += w1;
+ Dw(1,1) += -1.0*w1;
+ Dw(2,2) += -1.0*w1;
+ Dw(4,4) += w2;
+ Dw(5,5) += w2;
+
+ Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
+ Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
+ Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
+ Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
+ - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
+ Dz(0,0) = z1;
+ Dz(1,1) = z2;
+ Dz(2,2) = z3;
+ Dz(3,3) = z4;
+ Dz(4,4) = Dz(5,5) = z5;
+
+ G2 = Gamma2(U,T);
+ }
+ }
+ break;
+ case QQEE:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 2;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,0) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Lu;
+ double g12 = g_Ld;
+ //double g21 = g_Re;
+ double g22 = g_Re;
+
+ Complex w1 = 0.25*I*pi;
+ Dw(0,0) = Dw(1,1) = w1;
+
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ Dz(0,0) = z1;
+ Dz(1,1) = z2;
+
+ G2(0,0) = Gamma2Singlet()(0,0);
+ }
+ }
+ break;
+ case UUUU:
+ case UUUUiden:
+ case tRtRUU:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 2, numBrokenGauge = 2;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,1) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Ru;
+ //double g12 = g_Ru;
+ //double g21 = g_Ru;
+ double g22 = g_Ru;
+ // There is no Dw contribution for two SU(2) singlets.
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Dz(0,0) = Dz(1,1) = z1;
+ }
+ }
+ break;
+ case UUDD:
+ case tRtRDD:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 2, numBrokenGauge = 2;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,1) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Ru;
+ //double g12 = g_Ru;
+ //double g21 = g_Rd;
+ double g22 = g_Rd;
+ // There is no Dw contribution for two SU(2) singlets.
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Dz(0,0) = Dz(1,1) = z1;
+ }
+ }
+ break;
+ case UULL:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 2;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,0) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Lnu;
+ double g12 = g_Le;
+ //double g21 = g_Ru;
+ double g22 = g_Ru;
+
+ Complex w1 = 0.25*I*pi;
+ Dw(0,0) = Dw(1,1) = w1;
+
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ Dz(0,0) = z1;
+ Dz(1,1) = z2;
+
+ G2(0,0) = Gamma2Singlet()(0,0);
+ }
+ }
+ break;
+ case UUEE:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 1;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Ru;
+ //double g12 = g_Ru;
+ //double g21 = g_Re;
+ double g22 = g_Re;
+ // There is no Dw contribution for two SU(2) singlets.
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Dz(0,0) = z1;
+ }
+ }
+ break;
+ case DDDD:
+ case DDDDiden:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 2, numBrokenGauge = 2;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,1) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Rd;
+ //double g12 = g_Rd;
+ //double g21 = g_Rd;
+ double g22 = g_Rd;
+ // There is no Dw contribution for two SU(2) singlets.
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Dz(0,0) = Dz(1,1) = z1;
+ }
+ }
+ break;
+ case DDLL:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 2;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,0) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Lnu;
+ double g12 = g_Le;
+ //double g21 = g_Rd;
+ double g22 = g_Rd;
+
+ Complex w1 = 0.25*I*pi;
+ Dw(0,0) = Dw(1,1) = w1;
+
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ Dz(0,0) = z1;
+ Dz(1,1) = z2;
+
+ G2(0,0) = Gamma2Singlet()(0,0);
+ }
+ }
+ break;
+ case DDEE:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 1;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0 *= 0.0; Dw = Dz *= 0.0;
+ R0(0,0) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Rd;
+ //double g12 = g_Rd;
+ //double g21 = g_Re;
+ double g22 = g_Re;
+ // There is no Dw contribution for two SU(2) singlets.
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Dz(0,0) = z1;
+ }
+ }
+ break;
+ case LLLL:
+ case LLLLiden:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 2, numBrokenGauge = 6;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,1) = R0(1,1) = R0(2,1) = R0(3,1) = 1.0;
+ R0(0,0) = R0(3,0) = 0.25;
+ R0(1,0) = R0(2,0) = -0.25;
+ R0(4,0) = R0(5,0) = 0.5;
+ if (oneLoop) {
+ double g11 = g_Lnu;
+ double g12 = g_Le;
+ double g21 = g_Lnu;
+ double g22 = g_Le;
+
+ for (int i=0; i<6; i++) {
+ Dw(i,i) = 0.5*I*pi;
+ }
+ Complex w1 = (-1.0/2.0)*(T-U);
+ Complex w2 = (-1.0/2.0)*(T+U);
+ Dw(0,0) += w1;
+ Dw(3,3) += w1;
+ Dw(1,1) += -1.0*w1;
+ Dw(2,2) += -1.0*w1;
+ Dw(4,4) += w2;
+ Dw(5,5) += w2;
+
+ Complex z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21);
+ Complex z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12);
+ Complex z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11);
+ Complex z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ Complex z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U
+ - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22);
+ Dz(0,0) = z1;
+ Dz(1,1) = z2;
+ Dz(2,2) = z3;
+ Dz(3,3) = z4;
+ Dz(4,4) = Dz(5,5) = z5;
+
+ G2 = Gamma2(U,T);
+ }
+ }
+ break;
+ case LLEE:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 2;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,0) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Lnu;
+ double g12 = g_Le;
+ //double g21 = g_Re;
+ double g22 = g_Re;
+
+ Complex w1 = 0.25*I*pi;
+ Dw(0,0) = Dw(1,1) = w1;
+
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Complex z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22);
+ Dz(0,0) = z1;
+ Dz(1,1) = z2;
+
+ G2(0,0) = Gamma2Singlet()(0,0);
+ }
+ }
+ break;
+ case EEEE:
+ case EEEEiden:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 1;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = 1.0;
+ if (oneLoop) {
+ double g11 = g_Re;
+ //double g12 = g_Re;
+ //double g21 = g_Re;
+ double g22 = g_Re;
+ // There is no Dw contribution for two SU(2) singlets.
+ Complex z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22);
+ Dz(0,0) = z1;
+ }
+ }
+ break;
+ case QQWW:
+ case LLWW:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 5, numBrokenGauge = 20;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = 1.0; R0(0,1) = 0.5;
+ R0(1,0) = 1.0; R0(1,1) = -0.5;
+ R0(2,0) = cos2; R0(2,2) = -0.5*sin*cos; R0(2,3) = -0.5*sin*cos; R0(2,4) = sin2;
+ R0(3,0) = sin*cos; R0(3,2) = 0.5*cos2; R0(3,3) = -0.5*sin2; R0(3,4) = -sin*cos;
+ R0(4,0) = sin*cos; R0(4,2) = -0.5*sin2; R0(4,3) = 0.5*cos2; R0(4,4) = -sin*cos;
+ R0(5,0) = sin2; R0(5,2) = 0.5*sin*cos; R0(5,3) = 0.5*sin*cos; R0(5,4) = cos2;
+ R0(6,0) = 1.0; R0(6,1) = -0.5;
+ R0(7,0) = 1.0; R0(7,1) = 0.5;
+ R0(8,0) = cos2; R0(8,2) = 0.5*sin*cos; R0(8,3) = 0.5*sin*cos; R0(8,4) = sin2;
+ R0(9,0) = sin*cos; R0(9,2) = -0.5*cos2; R0(9,3) = 0.5*sin2; R0(9,4) = -sin*cos;
+ R0(10,0) = sin*cos; R0(10,2) = 0.5*sin2; R0(10,3) = -0.5*cos2; R0(10,4) = -sin*cos;
+ R0(11,0) = sin2; R0(11,2) = -0.5*sin*cos; R0(11,3) = -0.5*sin*cos; R0(11,4) = cos2;
+ R0(12,1) = -cos/sqrt(2.0); R0(12,3) = -sin/sqrt(2.0);
+ R0(13,1) = -sin/sqrt(2.0); R0(13,3) = cos/sqrt(2.0);
+ R0(14,1) = cos/sqrt(2.0); R0(14,2) = -sin/sqrt(2.0);
+ R0(15,1) = sin/sqrt(2.0); R0(15,2) = cos/sqrt(2.0);
+ R0(16,1) = -cos/sqrt(2.0); R0(16,2) = -sin/sqrt(2.0);
+ R0(17,1) = -sin/sqrt(2.0); R0(17,2) = cos/sqrt(2.0);
+ R0(18,1) = cos/sqrt(2.0); R0(18,3) = -sin/sqrt(2.0);
+ R0(19,1) = sin/sqrt(2.0); R0(19,3) = cos/sqrt(2.0);
+ if (oneLoop) {
+ double gW = g_W;
+ double g1(0.),g2(0.);
+ if (process==QQWW) {
+ g1 = g_Lu;
+ g2 = g_Ld;
+ }
+ else if (process==LLWW) {
+ g1 = g_Lnu;
+ g2 = g_Le;
+ }
+
+ Complex w1 = T-U+5.0/4.0*I*pi;
+ Complex w2 = -T+U+5.0/4.0*I*pi;
+ Complex w3 = -0.5*(T+U) + 3.0/4.0*I*pi;
+ Complex w4 = 0.25*I*pi;
+ Dw(0,0) = Dw(7,7) = w1;
+ Dw(1,1) = Dw(6,6) = w2;
+ for (unsigned int i=12; i<20; i++) {
+ Dw(i,i) = w3;
+ }
+ Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w4;
+ Dw(8,8) = Dw(9,9) = Dw(10,10) = Dw(11,11) = w4;
+
+ Complex z1 = 2.0*g1*gW*(U-T) - I*pi*(g1*g1+gW*gW);
+ Complex z2 = 2.0*g1*gW*(T-U) - I*pi*(g1*g1+gW*gW);
+ Complex z3 = 2.0*g2*gW*(U-T) - I*pi*(g2*g2+gW*gW);
+ Complex z4 = 2.0*g2*gW*(T-U) - I*pi*(g2*g2+gW*gW);
+ Complex z5 = -(g2*gW)*T + (g1*gW)*U - I*pi*(g1*g2+g1*gW-g2*gW);
+ Complex z6 = (g1*gW)*T - (g2*gW)*U - I*pi*(g1*g2+g1*gW-g2*gW);
+ Complex z7 = -I*pi*g1*g1;
+ Complex z8 = -I*pi*g2*g2;
+ Dz(0,0) = z1;
+ Dz(1,1) = z2;
+ Dz(2,2) = Dz(3,3) = Dz(4,4) = Dz(5,5) = z7;
+ Dz(6,6) = z3;
+ Dz(7,7) = z4;
+ Dz(8,8) = Dz(9,9) = Dz(10,10) = Dz(11,11) = z8;
+ Dz(12,12) = Dz(13,13) = Dz(16,16) = Dz(17,17) = z5;
+ Dz(14,14) = Dz(15,15) = Dz(18,18) = Dz(19,19) = z6;
+
+ G2 = Gamma2w(U,T);
+ }
+ }
+ break;
+ case QQPhiPhi:
+ case LLPhiPhi:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 2, numBrokenGauge = 14;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = 0.25; R0(0,1) = 1.0;
+ R0(1,0) = -1.0/8.0; R0(1,1) = 0.5;
+ R0(2,0) = I/8.0; R0(2,1) = -I/2.0;
+ R0(3,0) = -I/8.0; R0(3,1) = I/2.0;
+ R0(4,0) = -1.0/8.0; R0(4,1) = 1.0/2.0;
+ R0(5,0) = -1.0/4.0; R0(5,1) = 1.0;
+ R0(6,0) = 1.0/8.0; R0(6,1) = 1.0/2.0;
+ R0(7,0) = -I/8.0; R0(7,1) = -I/2.0;
+ R0(8,0) = I/8.0; R0(8,1) = I/2.0;
+ R0(9,0) = 1.0/8.0; R0(9,1) = 1.0/2.0;
+ R0(10,0) = -1.0/(2.0*sqrt(2.0));
+ R0(11,0) = I/(2.0*sqrt(2.0));
+ R0(12,0) = -1.0/(2.0*sqrt(2.0));
+ R0(13,0) = -I/(2.0*sqrt(2.0));
+
+ if (oneLoop) {
+ double g1(0.),g2(0.);
+ if (process==QQPhiPhi) {
+ g1 = g_Lu;
+ g2 = g_Ld;
+ }
+ else if (process==LLPhiPhi) {
+ g1 = g_Lnu;
+ g2 = g_Le;
+ }
+ else
+ assert(false);
+ double g3 = g_phiPlus;
+
+ Complex w0 = 0.25*I*pi;
+ Complex w1 = 0.5*(T-U) + 0.5*I*pi;
+ Complex w2 = -0.5*(T-U) + 0.5*I*pi;
+ Complex w3 = 0.25*I*(T-U);
+ Complex w4 = -0.25*(T+U) + 0.25*I*pi;
+ Dw(0,0) = w2;
+ Dw(1,1) = w0; Dw(1,2) = w3; Dw(1,3) = -w3; Dw(1,4) = w0;
+ Dw(2,1) = -w3; Dw(2,2) = w0; Dw(2,3) = -w0; Dw(2,4) = -w3;
+ Dw(3,1) = w3; Dw(3,2) = -w0; Dw(3,3) = w0; Dw(3,4) = w3;
+ Dw(4,1) = w0; Dw(4,2) = w3; Dw(4,3) = -w3; Dw(4,4) = w0;
+ Dw(5,5) = w1;
+ Dw(6,6) = w0; Dw(6,7) = -w3; Dw(6,8) = w3; Dw(6,9) = w0;
+ Dw(7,6) = w3; Dw(7,7) = w0; Dw(7,8) = -w0; Dw(7,9) = w3;
+ Dw(8,6) = -w3; Dw(8,7) = -w0; Dw(8,8) = w0; Dw(8,9) = -w3;
+ Dw(9,6) = w0; Dw(9,7) = -w3; Dw(9,8) = w3; Dw(9,9) = w0;
+ Dw(10,10) = w4; Dw(10,11) = I*w4;
+ Dw(11,10) = -I*w4; Dw(11,11) = w4;
+ Dw(12,12) = w4; Dw(12,13) = -I*w4;
+ Dw(13,12) = I*w4; Dw(13,13) = w4;
+
+ Complex z1 = 2.0*g3*g1*(T-U) - I*pi*(g3*g3+g1*g1);
+ Complex z2 = 2.0*g3*g2*(T-U) - I*pi*(g3*g3+g2*g2);
+ Complex z3 = -I*pi*g1*g1;
+ Complex z4 = 0.5*I*g1*(T-U);
+ Complex z5 = 0.25*I*pi;
+ Complex z6 = -I*pi*g2*g2;
+ Complex z7 = 0.5*I*g2*(T-U);
+ Complex z8 = g3*g1*T-g3*g2*U-I*pi*(g1*g2-g2*g3+g1*g3);
+ Complex z9 = 0.5*I*g2*T-0.5*I*g1*U+pi/2.0*g2-pi/2.0*g1+pi/2.0*g3;
+ Dz(0,0) = z1;
+ Dz(1,1) = z3; Dz(1,2) = -z4; Dz(1,3) = z4; Dz(1,4) = -z5;
+ Dz(2,1) = z4; Dz(2,2) = z3; Dz(2,3) = z5; Dz(2,4) = z4;
+ Dz(3,1) = -z4; Dz(3,2) = z5; Dz(3,3) = z3; Dz(3,4) = -z4;
+ Dz(4,1) = -z5; Dz(4,2) = -z4; Dz(4,3) = z4; Dz(4,4) = z3;
+ Dz(5,5) = z2;
+ Dz(6,6) = z6; Dz(6,7) = -z7; Dz(6,8) = z7; Dz(6,9) = -z5;
+ Dz(7,6) = z7; Dz(7,7) = z6; Dz(7,8) = z5; Dz(7,9) = z7;
+ Dz(8,6) = -z7; Dz(8,7) = z5; Dz(8,8) = z6; Dz(8,9) = -z7;
+ Dz(9,6) = -z5; Dz(9,7) = -z7; Dz(9,8) = z7; Dz(9,9) = z6;
+ Dz(10,10) = z8; Dz(10,11) = -z9;
+ Dz(11,10) = z9; Dz(11,11) = z8;
+ Dz(12,12) = z8; Dz(12,13) = z9;
+ Dz(13,12) = -z9; Dz(13,13) = z8;
+
+ G2 = Gamma2(U,T);
+ }
+ }
+ break;
+ case QQWG:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 6;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = 1.0/sqrt(2);
+ R0(1,0) = 1.0/sqrt(2);
+ R0(2,0) = cos/2.0;
+ R0(3,0) = sin/2.0;
+ R0(4,0) = -cos/2.0;
+ R0(5,0) = -sin/2.0;
+ if (oneLoop) {
+ double g1 = g_Lu;
+ double g2 = g_Ld;
+ double gW = g_W;
+
+ Complex w1 = -0.5*(T+U) + 0.75*I*pi;
+ Complex w2 = 0.25*I*pi;
+ Dw(0,0) = Dw(1,1) = w1;
+ Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;
+
+ Complex z1 = gW*g1*T - gW*g2*U - I*pi*(g1*g2+g1*gW-g2*gW);
+ Complex z2 = gW*g1*U - gW*g2*T - I*pi*(g2*g1+g1*gW-g2*gW);
+ Complex z3 = -I*pi*g1*g1;
+ Complex z4 = -I*pi*g2*g2;
+ Dz(0,0) = z1;
+ Dz(1,1) = z2;
+ Dz(2,2) = z3;
+ Dz(3,3) = z3;
+ Dz(4,4) = z4;
+ Dz(5,5) = z4;
+
+ G2(0,0) = -7.0/4.0*I*pi + (U+T);
+ }
+ }
+ break;
+ case QQBG:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 4;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = -sin;
+ R0(1,0) = cos;
+ R0(2,0) = -sin;
+ R0(3,0) = cos;
+ if (oneLoop) {
+ double g1 = g_Lu;
+ double g2 = g_Ld;
+ Complex w2 = 0.25*I*pi;
+ Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = w2;
+ Complex z3 = -I*pi*g1*g1;
+ Complex z4 = -I*pi*g2*g2;
+ Dz(0,0) = z3;
+ Dz(1,1) = z3;
+ Dz(2,2) = z4;
+ Dz(3,3) = z4;
+ G2(0,0) = Gamma2Singlet()(0,0);
+ }
+ }
+ break;
+ case QQGG:
+ case QtQtGG:
+ {
+ unsigned int numGauge = 3, numBrokenGauge = 6;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(3,0) = 1.0;
+ R0(1,1) = R0(4,1) = 1.0;
+ R0(2,2) = R0(5,2) = 1.0;
+ double g1 = g_Lu;
+ double g2 = g_Ld;
+ Complex w2(0.),z3(0.),z4(0.);
+ if (oneLoop) {
+ if(iswap==0) {
+ w2 = 0.25*I*pi;
+ z3 = -I*pi*g1*g1;
+ z4 = -I*pi*g2*g2;
+ G2(0,0) = G2(1,1) = G2(2,2) = Gamma2Singlet()(0,0);
+ }
+ else if(iswap==1) {
+ w2 = 0.25*(-T+I*pi);
+ z3 = (T-I*pi)*sqr(g1);
+ z4 = (T-I*pi)*sqr(g2);
+ G2(0,0) = G2(1,1) = G2(2,2) = Gamma2SingletST(T)(0,0);
+ }
+ else
+ assert(false);
+ Dw(0,0) = Dw(1,1) = Dw(2,2) = Dw(3,3) = Dw(4,4) = Dw(5,5) = w2;
+ Dz(0,0) = Dz(1,1) = Dz(2,2) = z3;
+ Dz(3,3) = Dz(4,4) = Dz(5,5) = z4;
+ }
+ }
+ break;
+ case UUBB:
+ case DDBB:
+ case EEBB:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 4;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = sin2;
+ R0(1,0) = -sin*cos;
+ R0(2,0) = -sin*cos;
+ R0(3,0) = cos2;
+ if (oneLoop) {
+ double g1(0.);
+ if (process==UUBB) {
+ g1 = g_Ru;
+ }
+ else if (process==DDBB) {
+ g1 = g_Rd;
+ }
+ else if (process==EEBB) {
+ g1 = g_Re;
+ }
+ else
+ assert(false);
+ // There is no Dw contribution for two SU(2) singlets.
+ Complex z1 = -I*pi*g1*g1;
+ Dz(0,0) = Dz(1,1) = Dz(2,2) = Dz(3,3) = z1;
+ }
+ }
+ break;
+ case UUPhiPhi:
+ case DDPhiPhi:
+ case EEPhiPhi:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 5;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = 1.0;
+ R0(1,0) = 0.5;
+ R0(2,0) = -0.5*I;
+ R0(3,0) = 0.5*I;
+ R0(4,0) = 0.5;
+ if (oneLoop) {
+ double g1(0.);
+ if (process==UUPhiPhi) {
+ g1 = g_Ru;
+ }
+ else if (process==DDPhiPhi) {
+ g1 = g_Rd;
+ }
+ else if (process==EEPhiPhi) {
+ g1 = g_Re;
+ }
+ double g3 = g_phiPlus;
+ Dw(0,0) = Dw(1,4) = Dw(4,1) = 0.25*I*pi;
+ Dw(2,3) = Dw(3,2) = -0.25*I*pi;
+ Complex z1 = 2.0*g3*g1*(T-U) - I*pi*(g3*g3+g1*g1);
+ Complex z2 = 0.5*I*g1*g1;
+ Complex z3 = -I*pi*g1*g1;
+ Complex z4 = 0.25*I*pi;
+ Dz(0,0) = z1;
+ Dz(1,1) = z3; Dz(1,2) = -z2; Dz(1,3) = z2; Dz(1,4) = -z4;
+ Dz(2,1) = z2; Dz(2,2) = z3; Dz(2,3) = z4; Dz(2,4) = z2;
+ Dz(3,1) = -z2; Dz(3,2) = z4; Dz(3,3) = z3; Dz(3,4) = -z2;
+ Dz(4,1) = -z4; Dz(4,2) = -z2; Dz(4,3) = z2; Dz(4,4) = z3;
+ G2(0,0) = Gamma2Singlet()(0,0);
+ }
+ }
+ break;
+ case UUBG:
+ case DDBG:
+ {
+ assert(iswap==0);
+ unsigned int numGauge = 1, numBrokenGauge = 2;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = -sin;
+ R0(1,0) = cos;
+ if (oneLoop) {
+ double g1(0.);
+ if (process==UUBG) {
+ g1 = g_Ru;
+ }
+ else if (process==DDBG) {
+ g1 = g_Rd;
+ }
+ else
+ assert(false);
+ // There is no Dw contribution for two SU(2) singlets.
+ Complex z1 = -I*pi*g1*g1;
+ Dz(0,0) = Dz(1,1) = z1;
+ }
+ }
+ break;
+ case UUGG:
+ case tRtRGG:
+ case DDGG:
+ {
+ unsigned int numGauge = 3, numBrokenGauge = 3;
+ R0 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numGauge);
+ G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ Dw = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Dz = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ R0(0,0) = R0(1,1) = R0(2,2) = 1.0;
+ if (oneLoop) {
+ double g1(0.);
+ if ((process==UUGG)||(process==tRtRGG)) {
+ g1 = g_Ru;
+ }
+ else if (process==DDGG) {
+ g1 = g_Rd;
+ }
+ else
+ assert(false);
+ Complex z1(0.);
+ // There is no Dw contribution for two SU(2) singlets.
+ if(iswap==0) {
+ z1 = -I*pi*sqr(g1);
+ }
+ else if(iswap==1) {
+ z1 = (T-I*pi)*sqr(g1);
+ }
+ else
+ assert(false);
+ Dz(0,0) = Dz(1,1) = Dz(2,2) = z1;
+ }
+ }
+ break;
+ default:
+ assert(false);
+ }
+
+ double aW = ElectroWeakReweighter::coupling()->aW(mu);
+ double aZ = ElectroWeakReweighter::coupling()->aZ(mu);
+ Energy mZ = ElectroWeakReweighter::coupling()->mZ();
+ Energy mW = ElectroWeakReweighter::coupling()->mW();
+
+ if (!oneLoop) {
+ return R0;
+ }
+ boost::numeric::ublas::matrix<Complex> output(R0);
+ boost::numeric::ublas::matrix<Complex> temp(R0.size1(),R0.size2());
+ boost::numeric::ublas::axpy_prod(R0,G2,temp);
+ output+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
+ boost::numeric::ublas::axpy_prod(Dw,R0,temp);
+ output+=aW/(4.0*pi)*4.0*log(mW/mu)*temp;
+ boost::numeric::ublas::axpy_prod(Dz,R0,temp);
+ output+=aZ/(4.0*pi)*4.0*log(mZ/mu)*temp;
+ return output;
}
diff --git a/MatrixElement/EW/ElectroWeakMatching.h b/MatrixElement/EW/ElectroWeakMatching.h
--- a/MatrixElement/EW/ElectroWeakMatching.h
+++ b/MatrixElement/EW/ElectroWeakMatching.h
@@ -1,59 +1,32 @@
// -*- C++ -*-
//
// ElectroWeakMatching.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_ElectroWeakMatching_H
#define HERWIG_ElectroWeakMatching_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 ElectroWeakMatching {
/**
* The high energy matching
*/
boost::numeric::ublas::matrix<Complex>
- electroWeakMatching(Energy mu,
- Energy2 s, Energy2 t, Energy2 u,
- Herwig::EWProcess::Process process,
- bool oneLoop);
-
- /* /\** */
- /* * 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); */
+ electroWeakMatching(Energy mu, Energy2 s, Energy2 t, Energy2 u,
+ Herwig::EWProcess::Process process,
+ bool oneLoop,unsigned int iswap);
}
}
#endif // HERWIG_ElectroWeakMatching_H
diff --git a/MatrixElement/EW/ElectroWeakReweighter.cc b/MatrixElement/EW/ElectroWeakReweighter.cc
--- a/MatrixElement/EW/ElectroWeakReweighter.cc
+++ b/MatrixElement/EW/ElectroWeakReweighter.cc
@@ -1,661 +1,907 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ElectroWeakReweighter class.
//
#include "ElectroWeakReweighter.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/operation.hpp"
#include "EWProcess.h"
#include "HighEnergyMatching.h"
#include "ElectroWeakMatching.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/epsilon.h"
using namespace Herwig;
tEWCouplingsPtr ElectroWeakReweighter::staticEWCouplings_ = tEWCouplingsPtr();
ElectroWeakReweighter::ElectroWeakReweighter() {}
ElectroWeakReweighter::~ElectroWeakReweighter() {}
IBPtr ElectroWeakReweighter::clone() const {
return new_ptr(*this);
}
IBPtr ElectroWeakReweighter::fullclone() const {
return new_ptr(*this);
}
void ElectroWeakReweighter::persistentOutput(PersistentOStream & os) const {
os << EWCouplings_ << collinearSudakov_ << softSudakov_;
}
void ElectroWeakReweighter::persistentInput(PersistentIStream & is, int) {
is >> EWCouplings_ >> collinearSudakov_ >> softSudakov_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<ElectroWeakReweighter,ReweightBase>
describeHerwigElectroWeakReweighter("Herwig::ElectroWeakReweighter", "HwMEEW.so");
void ElectroWeakReweighter::Init() {
static ClassDocumentation<ElectroWeakReweighter> documentation
("There is no documentation for the ElectroWeakReweighter class");
static Reference<ElectroWeakReweighter,EWCouplings> interfaceEWCouplings
("EWCouplings",
"The object to calculate the electroweak couplings",
&ElectroWeakReweighter::EWCouplings_, false, false, true, false, false);
static Reference<ElectroWeakReweighter,CollinearSudakov> interfaceCollinearSudakov
("CollinearSudakov",
"The collinear Sudakov",
&ElectroWeakReweighter::collinearSudakov_, false, false, true, false, false);
static Reference<ElectroWeakReweighter,SoftSudakov> interfaceSoftSudakov
("SoftSudakov",
"The soft Sudakov",
&ElectroWeakReweighter::softSudakov_, false, false, true, false, false);
}
namespace {
void axpy_prod_local(const boost::numeric::ublas::matrix<Complex> & A,
const boost::numeric::ublas::matrix<complex<InvEnergy2> > & B,
boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
assert(A.size2()==B.size1());
C.resize(A.size1(),B.size2());
for(unsigned int ix=0;ix<A.size1();++ix) {
for(unsigned int iy=0;iy<B.size2();++iy) {
C(ix,iy) = ZERO;
for(unsigned int iz=0;iz<A.size2();++iz) {
C(ix,iy) += A(ix,iz)*B(iz,iy);
}
}
}
}
void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> > & A,
const boost::numeric::ublas::vector<complex<Energy2> > & B,
boost::numeric::ublas::vector<Complex > & C) {
assert(A.size2()==B.size());
C.resize(A.size1());
for(unsigned int ix=0;ix<A.size1();++ix) {
C(ix) = ZERO;
for(unsigned int iz=0;iz<A.size2();++iz) {
C(ix) += A(ix,iz)*B(iz);
}
}
}
void axpy_prod_local(const boost::numeric::ublas::matrix<complex<InvEnergy2> > & A,
const boost::numeric::ublas::matrix<Complex> & B,
boost::numeric::ublas::matrix<complex<InvEnergy2> > & C) {
assert(A.size2()==B.size1());
C.resize(A.size1(),B.size2());
for(unsigned int ix=0;ix<A.size1();++ix) {
for(unsigned int iy=0;iy<B.size2();++iy) {
C(ix,iy) = ZERO;
for(unsigned int iz=0;iz<A.size2();++iz) {
C(ix,iy) += A(ix,iz)*B(iz,iy);
}
}
}
}
}
double ElectroWeakReweighter::weight() const {
EWCouplings_->initialize();
staticEWCouplings_ = EWCouplings_;
// cerr << "aEM\n";
// for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.1) {
// cerr << scale/GeV << " "
// << EWCouplings_->aEM(scale) << "\n";
// }
// cerr << "aS\n";
// for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
// cerr << scale/GeV << " "
// << EWCouplings_->aS(scale) << "\n";
// }
// cerr << "y_t\n";
// for(Energy scale=10.*GeV; scale<10*TeV; scale *= 1.4) {
// cerr << scale/GeV << " "
// << EWCouplings_->y_t(scale) << "\n";
// }
// cerr << "lambda\n";
// for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
// cerr << scale/GeV << " "
// << EWCouplings_->lambda(scale) << "\n";
// }
// cerr << "vev\n";
// for(Energy scale=91.2*GeV; scale<10*TeV; scale *= 1.4) {
// cerr << scale/GeV << " "
// << EWCouplings_->vev(scale)/GeV << "\n";
// }
// collinearSudakov_->makePlots();
// Energy2 s = sqr(5000.*GeV);
// Energy2 t = -0.25*s;
// Energy2 u = -0.75*s;
// testEvolution(s,t,u);
cerr << subProcess() << "\n";
cerr << *subProcess() << "\n";
cerr << subProcess()->outgoing()[0] << *subProcess()->outgoing()[0] << "\n";
cerr << subProcess()->outgoing()[0]->spinInfo() << "\n";
cerr << subProcess()->outgoing()[0]->spinInfo()->productionVertex() << "\n";
if(subProcess()->outgoing().size()!=2)
return 1.;
+ // processes with gg initial-state
if(subProcess()->incoming().first->id()==ParticleID::g &&
subProcess()->incoming().second->id()==ParticleID::g) {
if(subProcess()->outgoing()[0]->id()==ParticleID::g &&
subProcess()->outgoing()[1]->id()==ParticleID::g)
return 1.;
else if(abs(subProcess()->outgoing()[0]->id())<=5 &&
subProcess()->outgoing()[0]->id()==-subProcess()->outgoing()[1]->id()) {
return reweightggqqbar();
}
else
assert(false);
}
+ // processes with q qbar initial-state
else if(abs(subProcess()->incoming().first->id())<=5 &&
subProcess()->incoming().first->id()==-subProcess()->incoming().second->id()) {
if(subProcess()->outgoing()[0]->id()==ParticleID::g &&
subProcess()->outgoing()[1]->id()==ParticleID::g)
return reweightqqbargg();
else
assert(false);
}
+ // processes with q g initial-state
+ else if((subProcess()->incoming().first ->id()> 0 &&
+ subProcess()->incoming().first ->id()<=5 &&
+ subProcess()->incoming().second->id()==ParticleID::g) ||
+ (subProcess()->incoming().second->id()> 0 &&
+ subProcess()->incoming().second->id()<=5 &&
+ subProcess()->incoming().first ->id()==ParticleID::g)) {
+ // qg -> qg
+ if((subProcess()->outgoing()[0]->id()> 0 &&
+ subProcess()->outgoing()[0]->id()<=5 &&
+ subProcess()->outgoing()[1]->id()==ParticleID::g) ||
+ (subProcess()->outgoing()[1]->id()> 0 &&
+ subProcess()->outgoing()[1]->id()<=5 &&
+ subProcess()->outgoing()[0]->id()==ParticleID::g))
+ return reweightqgqg();
+ // unknown
+ else
+ assert(false);
+ }
+ // processes with qbar g initial-state
+ else if((subProcess()->incoming().first ->id()>=-5 &&
+ subProcess()->incoming().first ->id()< 0 &&
+ subProcess()->incoming().second->id()==ParticleID::g) ||
+ (subProcess()->incoming().second->id()>=-5 &&
+ subProcess()->incoming().second->id()< 0 &&
+ subProcess()->incoming().first ->id()==ParticleID::g)) {
+ if((subProcess()->outgoing()[0]->id()>=-5 &&
+ subProcess()->outgoing()[0]->id()< 0 &&
+ subProcess()->outgoing()[1]->id()==ParticleID::g) ||
+ (subProcess()->outgoing()[1]->id()>=-5 &&
+ subProcess()->outgoing()[1]->id()< 0 &&
+ subProcess()->outgoing()[0]->id()==ParticleID::g))
+ return reweightqbargqbarg();
+ else
+ assert(false);
+ }
+ // unknown initial-state
else
assert(false);
assert(false);
staticEWCouplings_ = tEWCouplingsPtr();
}
void ElectroWeakReweighter::testEvolution(Energy2 s,Energy2 t, Energy2 u) const {
Energy highScale = sqrt(s);
Energy ewScale = coupling()->mZ();
Energy lowScale = 50.0*GeV;
for (unsigned int i=0; i<45;++i) {
EWProcess::Process process = (EWProcess::Process)i;
cerr << "process " << process << "\n";
// result for all EW and QCD SCET contributions:
boost::numeric::ublas::matrix<complex<InvEnergy2> > highMatch_val
= HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,true,true);
boost::numeric::ublas::matrix<Complex> highRunning_val
- = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process);
+ = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process,0);
boost::numeric::ublas::matrix<Complex> ewMatch_val =
- ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true);
+ ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true,0);
boost::numeric::ublas::matrix<Complex> lowRunning_val =
- softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process);
+ softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process,0);
boost::numeric::ublas::matrix<Complex> collinearHighRunning_val =
collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false);
boost::numeric::ublas::matrix<Complex> collinearEWMatch_val =
collinearSudakov_->electroWeakMatching(ewScale,s,process,true);
boost::numeric::ublas::matrix<Complex> collinearLowRunning_val =
collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process);
boost::numeric::ublas::matrix<Complex> lowMatchTemp_val =
boost::numeric::ublas::zero_matrix<Complex>(ewMatch_val.size1(),ewMatch_val.size2());
for (unsigned int ii=0; ii<ewMatch_val.size1(); ++ii) {
for (unsigned int jj=0; jj<ewMatch_val.size2(); ++jj) {
lowMatchTemp_val(ii,jj) = collinearEWMatch_val(ii,jj)*ewMatch_val(ii,jj);
}
}
boost::numeric::ublas::matrix<Complex> temp(highRunning_val.size1(),collinearHighRunning_val.size2());
boost::numeric::ublas::axpy_prod(highRunning_val,collinearHighRunning_val,temp);
boost::numeric::ublas::matrix<Complex> temp2(collinearLowRunning_val.size1(),lowRunning_val.size2());
boost::numeric::ublas::axpy_prod(collinearLowRunning_val,lowRunning_val,temp2);
boost::numeric::ublas::matrix<Complex> temp3(temp2.size1(),lowMatchTemp_val.size2());
boost::numeric::ublas::axpy_prod(temp2,lowMatchTemp_val,temp3);
temp2.resize(temp3.size1(),temp.size2());
boost::numeric::ublas::axpy_prod(temp3,temp,temp2);
boost::numeric::ublas::matrix<complex<InvEnergy2> > result(temp2.size1(),highMatch_val.size2());
axpy_prod_local(temp2,highMatch_val,result);
for(unsigned int ix=0;ix<result.size1();++ix) {
for(unsigned int iy=0;iy<result.size2();++iy) {
cerr << s*result(ix,iy) << " ";
}
cerr << "\n";
}
}
}
namespace {
void SackGluonPolarizations(Lorentz5Momentum &p1,
Lorentz5Momentum &p2,
Lorentz5Momentum &p3,
Lorentz5Momentum &p4,
Energy2 s, Energy2 t, Energy2 u, Energy2 m2,
vector<LorentzVector<Complex> > & eps3,
vector<LorentzVector<Complex> > & eps4,
unsigned int iopt) {
static const Complex I(0.,1.);
// p1 is p-, p2 is p+
// p3 is k-, p4 is k+
// both final-state
if(iopt==0) {
// swap t and u due Aneesh's defn
Energy3 den1 = sqrt((u*t-sqr(m2))*(s-4.*m2));
Energy3 den2 = sqrt(s*(u*t-sqr(m2)));
LorentzVector<Complex> eps3Para = (m2+t)/den1*p1 -(m2+u)/den1*p2 +(u-t)/den1*p3;
LorentzVector<Complex> eps3Perp = 2./den2*epsilon(p1,p2,p3);
LorentzVector<Complex> eps4Para = (m2+t)/den1*p2 -(m2+u)/den1*p1 +(u-t)/den1*p4;
LorentzVector<Complex> eps4Perp = 2./den2*epsilon(p1,p2,p4);
eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
if(m2!=ZERO) assert(false);
}
// both initial-state
else if(iopt==1) {
if(m2!=ZERO) assert(false);
LorentzVector<Complex> eps3Para( 1., 0.,0.,0.);
LorentzVector<Complex> eps3Perp( 0.,-1.,0.,0.);
LorentzVector<Complex> eps4Para(-1.,0.,0., 0.);
LorentzVector<Complex> eps4Perp( 0., 1.,0.,0.);
eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
}
+ else if(iopt==2) {
+ // rotation into the 2,3 Breit frame
+ Lorentz5Momentum pa = p3-p2;
+ Axis axis(pa.vect().unit());
+ LorentzRotation rot;
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ if ( sinth > 1.e-9 )
+ rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ rot.rotateX(Constants::pi);
+ rot.boostZ( pa.e()/pa.vect().mag());
+ Lorentz5Momentum ptemp=rot*p2;
+ Boost trans = -1./ptemp.e()*ptemp.vect();
+ trans.setZ(0.);
+ rot.boost(trans);
+ LorentzVector<Complex> eps3Para( 1., 0.,0.,0.);
+ LorentzVector<Complex> eps3Perp( 0.,-1.,0.,0.);
+ LorentzVector<Complex> eps4Para(-1.,0.,0., 0.);
+ LorentzVector<Complex> eps4Perp( 0., 1.,0.,0.);
+ eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
+ eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
+ eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
+ eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
+ rot = rot.invert();
+ for(unsigned int ix=0;ix<2;++ix) {
+ eps3[ix] *=rot;
+ eps4[ix] *=rot;
+ }
+ }
+ else if(iopt==3) {
+ // rotation into the 1,4 Breit frame
+ Lorentz5Momentum pa = p4-p1;
+ Axis axis(pa.vect().unit());
+ LorentzRotation rot;
+ double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
+ if ( sinth > 1.e-9 )
+ rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ rot.rotateX(Constants::pi);
+ rot.boostZ( pa.e()/pa.vect().mag());
+ Lorentz5Momentum ptemp=rot*p1;
+ Boost trans = -1./ptemp.e()*ptemp.vect();
+ trans.setZ(0.);
+ rot.boost(trans);
+ LorentzVector<Complex> eps3Para( 1., 0.,0.,0.);
+ LorentzVector<Complex> eps3Perp( 0.,-1.,0.,0.);
+ LorentzVector<Complex> eps4Para(-1.,0.,0., 0.);
+ LorentzVector<Complex> eps4Perp( 0., 1.,0.,0.);
+ eps3.push_back(sqrt(0.5)*(eps3Para+I*eps3Perp));
+ eps3.push_back(sqrt(0.5)*(eps3Para-I*eps3Perp));
+ eps4.push_back(sqrt(0.5)*(eps4Para+I*eps4Perp));
+ eps4.push_back(sqrt(0.5)*(eps4Para-I*eps4Perp));
+ rot = rot.invert();
+ for(unsigned int ix=0;ix<2;++ix) {
+ eps3[ix] *=rot;
+ eps4[ix] *=rot;
+ }
+ }
+ else
+ assert(false);
}
}
double ElectroWeakReweighter::reweightqqbargg() const {
// momenta and invariants
Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum();
tcPDPtr q = subProcess()->incoming().first ->dataPtr();
Lorentz5Momentum p2 = subProcess()->incoming().second->momentum();
tcPDPtr qbar = subProcess()->incoming().second->dataPtr();
if(subProcess()->incoming().first->id()<0) {
swap(p1,p2 );
swap(q ,qbar);
}
Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
tcPDPtr g = subProcess()->outgoing()[1]->dataPtr();
Energy2 s = (p1+p2).m2();
Energy2 t = (p1-p4).m2();
Energy2 u = (p1-p3).m2();
- // // boost to partonci rest frame
- // Lorentz5Momentum psum=p1+p2;
- // LorentzRotation boost(-psum.boostVector());
- // p1 *= boost;
- // p2 *= boost;
- // p3 *= boost;
- // p4 *= boost;
- // cerr << "testing momenta in reweight A " << p1/GeV << "\n";
- // cerr << "testing momenta in reweight B " << p2/GeV << "\n";
- // cerr << "testing momenta in reweight C " << p3/GeV << "\n";
- // cerr << "testing momenta in reweight D " << p4/GeV << "\n";
- // LO matrix element coefficents
+ // boost to partonci rest frame
+ Lorentz5Momentum psum=p1+p2;
+ LorentzRotation boost(-psum.boostVector());
+ p1 *= boost;
+ p2 *= boost;
+ p3 *= boost;
+ p4 *= boost;
+ // LO and EW corrected matrix element coefficients
boost::numeric::ublas::matrix<complex<InvEnergy2> >
- bornQQGGweights,bornRRGGweights;
+ bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights;
// quark left doublet
if(q->id()!=5) {
- bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true);
+ bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true ,0);
+ EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false,0);
}
else {
- bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true);
+ bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true ,0);
+ EWQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,false,0);
}
// quark right singlet
- if(abs(subProcess()->incoming().first->id())%2==0)
- bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true);
- else
- bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true);
- // EW corrected matrix element coefficients
- boost::numeric::ublas::matrix<complex<InvEnergy2> >
- EWQQGGweights,EWRRGGweights;
- // quark left doublet
- if(q->id()!=5) {
- EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false);
+ if(abs(subProcess()->incoming().first->id())%2==0) {
+ bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true ,0);
+ EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false,0);
}
else {
- EWQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,false);
+ bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,0);
+ EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false,0);
}
- // quakr right singlet
- if(abs(subProcess()->incoming().first->id())%2==0)
- EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false);
- else
- EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false);
- // cerr << "testing matrices\n";
- // for(unsigned int ix=0;ix<bornRRGGweights.size1();++ix) {
- // for(unsigned int iy=0;iy<bornRRGGweights.size2();++iy) {
- // cerr << bornRRGGweights(ix,iy)*GeV2 << " ";
- // }
- // cerr << "\n";
- // }
- // for(unsigned int ix=0;ix<bornQQGGweights.size1();++ix) {
- // for(unsigned int iy=0;iy<bornQQGGweights.size2();++iy) {
- // cerr << bornQQGGweights(ix,iy)*GeV2 << " ";
- // }
- // cerr << "\n";
- // }
-
SpinorWaveFunction qw(p1,q ,incoming);
SpinorBarWaveFunction qbarw(p2,qbar,incoming);
- // VectorWaveFunction g1w(p3,getParticleData(ParticleID::g),outgoing);
- // VectorWaveFunction g2w(p4,getParticleData(ParticleID::g),outgoing);
vector<LorentzVector<Complex> > eps3,eps4;
SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps3,eps4,0);
boost::numeric::ublas::matrix<Complex>
bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
EWME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
- // testME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
for(unsigned int iq=0;iq<2;++iq) {
if(iq==0) {
qw.reset (0);
qbarw.reset(1);
}
else {
qw.reset (1);
qbarw.reset(0);
}
LorentzVector<complex<Energy> > current = iq==0 ?
qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
for(unsigned int i1=0;i1<2;++i1) {
complex<Energy> d31 = eps3[i1].dot(p1);
for(unsigned int i2=0;i2<2;++i2) {
// g1w.reset(2*i1);
// g2w.reset(2*i2);
boost::numeric::ublas::vector<complex<Energy2> > M(5);
Complex d34 = eps3[i1].dot(eps4[i2]);
complex<Energy> d42 = eps4[i2].dot(p2);
// M0 in paper
M(0) = qw.dimensionedWave().slash(eps3[i1])
.slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps4[i2]);
-
- // // really t channel
- // Complex MEU = -Complex(0.,1.)*qw.dimensionedWave().slash(g1w.wave())
- // .slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(g2w.wave())/u;
- // // really u channel
- // Complex MET = -Complex(0.,1.)*qw.dimensionedWave().slash(g2w.wave())
- // .slash(p3-p2).vectorCurrent(qbarw.dimensionedWave()).dot(g1w.wave())/t;
- // // s channel
- // Complex MES = -Complex(0.,1.)*(current.dot(p3-p4)*g1w.wave().dot(g2w.wave())
- // -2.*current.dot(g1w.wave())*(g2w.wave().dot(p3))
- // -2.*current.dot(g2w.wave())*(g1w.wave().dot(p4)))/s;
-
- // // really t channel
- // Complex MEU = -Complex(0.,1.)*qw.dimensionedWave().slash(eps3[i1])
- // .slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps4[i2])/u;
- // // really u channel
- // Complex MET = -Complex(0.,1.)*qw.dimensionedWave().slash(eps4[i2])
- // .slash(p3-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps3[i1])/t;
- // // s channel
- // Complex MES = -Complex(0.,1.)*(current.dot(p3-p4)*eps3[i1].dot(eps4[i2])
- // -2.*current.dot(eps3[i1])*(eps4[i2].dot(p3))
- // -2.*current.dot(eps4[i2])*(eps3[i1].dot(p4)))/s;
- // cerr << "NEW flows " << MEU+MES << " " << MET-MES << "\n";
-
-
-
-
- // cerr << "testing new U " << MET << "\n";
// M4 in paper
M(2) = current.dot(eps4[i2])*d31;
// M5 in paper
M(3) = -current.dot(eps3[i1])*d42;
// M1 in paper (missing factor)
M(1) = current.dot(p4);
// M6 in paper
M(4) = M(1)*d31*d42/GeV2;
// M1 final factor
M(1) *= d34;
// coefficient of different contributions
boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3),Ctest(3);
// Ctest(0) = 1./6.*( MEU+MET);
// Ctest(1) = 0.5*( MEU+MET);
// Ctest(2) = 0.5*(MEU+MES-MET+MES);
if(iq==0) {
axpy_prod_local(bornQQGGweights,M,Cborn);
axpy_prod_local(EWQQGGweights ,M,CEW );
}
else {
axpy_prod_local(bornRRGGweights,M,Cborn);
axpy_prod_local(EWRRGGweights ,M,CEW );
}
unsigned int ioff = (Cborn.size()==6 && q->id()%2!=0) ? 3 : 0;
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
bornME(ix,iy) += Cborn(ix+ioff)*conj(Cborn(iy+ioff));
EWME (ix,iy) += CEW (ix+ioff)*conj(CEW (iy+ioff));
- // testME (ix,iy) += Ctest (ix)*conj(Ctest (iy));
}
}
}
}
}
double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
double EW = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
- // double test = 24.*real(testME(0,0))+20./3.*real(testME(1,1))+12.*real(testME(2,2));
-
- // double gs2 = 4.*Constants::pi*ElectroWeakReweighter::coupling()->a3(sqrt(s));
-
- // cerr << "testing born A " << 0.125*born/sqr(gs2)/9. << "\n";
- // cerr << "testing born B " << 0.125*test/9. << "\n";
-
return EW/born;
}
boost::numeric::ublas::matrix<complex<InvEnergy2> >
ElectroWeakReweighter::evaluateRunning(EWProcess::Process process, Energy2 s,
- Energy2 t, Energy2 u, bool born) const {
+ Energy2 t, Energy2 u, bool born,
+ unsigned int iswap) const {
using namespace boost::numeric::ublas;
bool SU3save = coupling()->SU3();
bool EWsave = coupling()-> EW();
Energy highScale = sqrt(s);
Energy ewScale = coupling()->mZ();
Energy lowScale = ewScale;
// result for all EW and QCD SCET contributions:
+ // MATCHING CONTRIBUTIONS
// high energy matching
- matrix<complex<InvEnergy2> > highMatch_val
- = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,!born,false);
+ matrix<complex<InvEnergy2> > highMatch_val;
+ if(iswap==0)
+ highMatch_val = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,!born,false);
+ else if(iswap==1)
+ highMatch_val = HighEnergyMatching::highEnergyMatching(highScale,t,s,u,process,!born,false);
+ else
+ assert(false);
// low energy matching
matrix<Complex>
- ewMatch_val = ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,!born);
+ ewMatch_val = ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,!born,iswap);
matrix<Complex> collinearEWMatch_val =
collinearSudakov_->electroWeakMatching(ewScale,s,process,!born);
+ // EVOLUTION
matrix<Complex> highRunning_val,lowRunning_val,
collinearHighRunning_val,collinearLowRunning_val;
+ // born process
if(born) {
highRunning_val = identity_matrix<Complex>(softSudakov_->numberGauge(process));
lowRunning_val = identity_matrix<Complex>(softSudakov_->numberBrokenGauge(process));
collinearHighRunning_val = identity_matrix<Complex>(softSudakov_->numberGauge(process));
collinearLowRunning_val = identity_matrix<Complex>(softSudakov_->numberBrokenGauge(process));
}
+ // EW corrected
else {
coupling()->SU3(false);
coupling()-> EW( true);
- highRunning_val = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process);
- lowRunning_val = softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process);
+ highRunning_val = softSudakov_->highEnergyRunning(highScale, ewScale,s,t,u,process,iswap);
+ lowRunning_val = softSudakov_->lowEnergyRunning ( ewScale,lowScale,s,t,u,process,iswap);
collinearHighRunning_val = collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false);
collinearLowRunning_val = collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process);
};
matrix<Complex> lowMatchTemp_val =
zero_matrix<Complex>(ewMatch_val.size1(),ewMatch_val.size2());
for (unsigned int ii=0; ii<ewMatch_val.size1(); ++ii) {
for (unsigned int jj=0; jj<ewMatch_val.size2(); ++jj) {
lowMatchTemp_val(ii,jj) = collinearEWMatch_val(ii,jj)*ewMatch_val(ii,jj);
}
}
// perform all the multiplications
matrix<Complex> temp(highRunning_val.size1(),collinearHighRunning_val.size2());
axpy_prod(highRunning_val,collinearHighRunning_val,temp);
matrix<Complex> temp2(collinearLowRunning_val.size1(),lowRunning_val.size2());
axpy_prod(collinearLowRunning_val,lowRunning_val,temp2);
matrix<Complex> temp3(temp2.size1(),lowMatchTemp_val.size2());
axpy_prod(temp2,lowMatchTemp_val,temp3);
temp2.resize(temp3.size1(),temp.size2());
axpy_prod(temp3,temp,temp2);
matrix<complex<InvEnergy2> > result(temp2.size1(),highMatch_val.size2());
axpy_prod_local(temp2,highMatch_val,result);
- // for(unsigned int ix=0;ix<result.size1();++ix) {
- // for(unsigned int iy=0;iy<result.size2();++iy) {
- // cerr << s*result(ix,iy) << " ";
- // }
- // cerr << "\n";
- // }
// reset the couplings
coupling()->SU3(SU3save);
coupling()-> EW( EWsave);
// return the answer
return result;
}
double ElectroWeakReweighter::reweightggqqbar() const {
// momenta and invariants
Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum();
Lorentz5Momentum p2 = subProcess()->incoming().second->momentum();
Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
tcPDPtr qbar = subProcess()->outgoing()[0]->dataPtr();
tcPDPtr q = subProcess()->outgoing()[1]->dataPtr();
if(q->id()<0) {
swap(p3,p4 );
swap(q ,qbar);
}
Energy2 s = (p1+p2).m2();
Energy2 t = (p1-p4).m2();
Energy2 u = (p1-p3).m2();
// boost to partonic rest frame and rescale momenta of outgoing
// so zero mass
Lorentz5Momentum psum=p1+p2;
LorentzRotation boost(-psum.boostVector());
p1 *= boost;
p2 *= boost;
p3 *= boost;
p4 *= boost;
p3.setMass(ZERO);
p3.rescaleRho();
p4.setMass(ZERO);
p4.rescaleRho();
- // LO matrix element coefficents
+ // LO and EW matrix element coefficents
boost::numeric::ublas::matrix<complex<InvEnergy2> >
- bornQQGGweights,bornRRGGweights;
+ bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights;
// quark left doublet
if(q->id()<5) {
- bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true);
+ bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true ,0);
+ EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false,0);
}
else {
- bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true);
+ bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true ,0);
+ EWQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,false,0);
}
// quark right singlet
if(q->id()==0) {
- if(q->id()==6)
- bornRRGGweights = evaluateRunning(EWProcess::tRtRGG,s,t,u,true);
- else
- bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true);
- }
- else
- bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true);
- // EW corrected matrix element coefficients
- boost::numeric::ublas::matrix<complex<InvEnergy2> >
- EWQQGGweights,EWRRGGweights;
- // quark left doublet
- if(q->id()<5) {
- EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false);
+ if(q->id()==6) {
+ bornRRGGweights = evaluateRunning(EWProcess::tRtRGG,s,t,u,true ,0);
+ EWRRGGweights = evaluateRunning(EWProcess::tRtRGG,s,t,u,false,0);
+ }
+ else {
+ bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true ,0);
+ EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false,0);
+ }
}
else {
- EWQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,false);
+ bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,0);
+ EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false,0);
}
- // quark right singlet
- if(q->id()%2==0) {
- if(q->id()==6)
- EWRRGGweights = evaluateRunning(EWProcess::tRtRGG,s,t,u,false);
- else
- EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false);
- }
- else
- EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false);
SpinorWaveFunction qw(p4,qbar,incoming);
SpinorBarWaveFunction qbarw(p3,q ,incoming);
vector<LorentzVector<Complex> > eps1,eps2;
SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps1,eps2,1);
boost::numeric::ublas::matrix<Complex>
bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
EWME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
// helicities of outgoing quarks
for(unsigned int iq=0;iq<2;++iq) {
if(iq==0) {
qw.reset (0);
qbarw.reset(1);
}
else {
qw.reset (1);
qbarw.reset(0);
}
LorentzVector<complex<Energy> > current = iq==0 ?
qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
for(unsigned int i1=0;i1<2;++i1) {
complex<Energy> d31 = eps1[i1].dot(p3);
for(unsigned int i2=0;i2<2;++i2) {
boost::numeric::ublas::vector<complex<Energy2> > M(5);
Complex d34 = eps1[i1].dot(eps2[i2]);
complex<Energy> d42 = eps2[i2].dot(p4);
// M0 in paper
M(0) = qw.dimensionedWave().slash(eps1[i1])
.slash(p2-p4).vectorCurrent(qbarw.dimensionedWave()).dot(eps2[i2]);
// M4 in paper
M(2) = current.dot(eps2[i2])*d31;
// M5 in paper
M(3) = -current.dot(eps1[i1])*d42;
// M1 in paper (missing factor)
M(1) = current.dot(p2);
// M6 in paper
M(4) = M(1)*d31*d42/GeV2;
// M1 final factor
M(1) *= d34;
// coefficient of different contributions
boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3);
if(iq==0) {
axpy_prod_local(bornQQGGweights,M,Cborn);
axpy_prod_local(EWQQGGweights ,M,CEW );
}
else {
axpy_prod_local(bornRRGGweights,M,Cborn);
axpy_prod_local(EWRRGGweights ,M,CEW );
}
unsigned int ioff = (Cborn.size()==6 && q->id()%2!=0) ? 3 : 0;
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
bornME(ix,iy) += Cborn(ix+ioff)*conj(Cborn(iy+ioff));
EWME (ix,iy) += CEW (ix+ioff)*conj(CEW (iy+ioff));
}
}
}
}
}
double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
double EW = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
return EW/born;
}
+
+double ElectroWeakReweighter::reweightqgqg() const {
+ // momenta and invariants
+ Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum();
+ Lorentz5Momentum p2 = subProcess()->incoming().second->momentum();
+ tcPDPtr q;
+ if(subProcess()->incoming().first->id()!=ParticleID::g) {
+ q = subProcess()->incoming().first ->dataPtr();
+ }
+ else {
+ q = subProcess()->incoming().second->dataPtr();
+ swap(p1,p2);
+ }
+ Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
+ Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
+ if(subProcess()->outgoing()[0]->id()!=ParticleID::g)
+ swap(p3,p4);
+ Energy2 s = (p1+p2).m2();
+ Energy2 t = (p1-p4).m2();
+ Energy2 u = (p1-p3).m2();
+ // boost to partonic rest frame
+ Lorentz5Momentum psum=p1+p2;
+ LorentzRotation boost(-psum.boostVector());
+ p1 *= boost;
+ p2 *= boost;
+ p3 *= boost;
+ p4 *= boost;
+ // LO and EW corrected matrix element coefficients
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights;
+ // quark left doublet
+ if(q->id()!=5) {
+ bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true ,1);
+ EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false,1);
+ }
+ else {
+ bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true ,1);
+ EWQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,false,1);
+ }
+ // quark right singlet
+ if(abs(q->id())%2==0) {
+ bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true ,1);
+ EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false,1);
+ }
+ else {
+ bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,1);
+ EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false,1);
+ }
+ SpinorWaveFunction qw(p1,q,incoming);
+ SpinorBarWaveFunction qbarw(p4,q,outgoing);
+ vector<LorentzVector<Complex> > eps2,eps3;
+ SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps2,eps3,2);
+ // compute the matrix elements
+ boost::numeric::ublas::matrix<Complex>
+ bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
+ EWME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
+ testME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
+ for(unsigned int iq=0;iq<2;++iq) {
+ if(iq==0) {
+ qw.reset (0);
+ qbarw.reset(0);
+ }
+ else {
+ qw.reset (1);
+ qbarw.reset(1);
+ }
+ LorentzVector<complex<Energy> > current = iq==0 ?
+ qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
+ qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
+ for(unsigned int i1=0;i1<2;++i1) {
+ complex<Energy> d31 = eps3[i1].dot(p1);
+ for(unsigned int i2=0;i2<2;++i2) {
+ boost::numeric::ublas::vector<complex<Energy2> > M(5);
+ Complex d34 = eps3[i1].dot(eps2[i2]);
+ complex<Energy> d42 = eps2[i2].dot(p4);
+ // M0 in paper
+ M(0) = qw.dimensionedWave().slash(eps3[i1])
+ .slash(p2-p4).vectorCurrent(qbarw.dimensionedWave()).dot(eps2[i2]);
+ // M4 in paper
+ M(2) = current.dot(eps2[i2])*d31;
+ // M5 in paper
+ M(3) = -current.dot(eps3[i1])*d42;
+ // M1 in paper (missing factor)
+ M(1) = current.dot(p2);
+ // M6 in paper
+ M(4) = M(1)*d31*d42/GeV2;
+ // M1 final factor
+ M(1) *= d34;
+ // coefficient of different contributions
+ boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3);
+ if(iq==0) {
+ axpy_prod_local(bornQQGGweights,M,Cborn);
+ axpy_prod_local(EWQQGGweights ,M,CEW );
+ }
+ else {
+ axpy_prod_local(bornRRGGweights,M,Cborn);
+ axpy_prod_local(EWRRGGweights ,M,CEW );
+ }
+ unsigned int ioff = (Cborn.size()==6 && q->id()%2!=0) ? 3 : 0;
+ for(unsigned int ix=0;ix<3;++ix) {
+ for(unsigned int iy=0;iy<3;++iy) {
+ bornME(ix,iy) += Cborn(ix+ioff)*conj(Cborn(iy+ioff));
+ EWME (ix,iy) += CEW (ix+ioff)*conj(CEW (iy+ioff));
+ }
+ }
+ }
+ }
+ }
+ double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
+ double EW = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
+ return EW/born;
+}
+
+double ElectroWeakReweighter::reweightqbargqbarg() const {
+ // momenta and invariants
+ Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum();
+ Lorentz5Momentum p2 = subProcess()->incoming().second->momentum();
+ tcPDPtr qbar;
+ if(subProcess()->incoming().first->id()==ParticleID::g) {
+ qbar = subProcess()->incoming().second->dataPtr();
+ }
+ else {
+ qbar = subProcess()->incoming().first ->dataPtr();
+ swap(p1,p2);
+ }
+ Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum();
+ Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum();
+ if(subProcess()->outgoing()[0]->id()==ParticleID::g)
+ swap(p3,p4);
+ Energy2 s = (p1+p2).m2();
+ Energy2 t = (p1-p4).m2();
+ Energy2 u = (p1-p3).m2();
+ // boost to partonci rest frame
+ Lorentz5Momentum psum=p1+p2;
+ LorentzRotation boost(-psum.boostVector());
+ p1 *= boost;
+ p2 *= boost;
+ p3 *= boost;
+ p4 *= boost;
+ // LO and EW corrected matrix element coefficients
+ boost::numeric::ublas::matrix<complex<InvEnergy2> >
+ bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights;
+ // quark left doublet
+ if(qbar->id()!=-5) {
+ bornQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,true ,1);
+ EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false,1);
+ }
+ else {
+ bornQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,true ,1);
+ EWQQGGweights = evaluateRunning(EWProcess::QtQtGG,s,t,u,false,1);
+ }
+ // quark right singlet
+ if(abs(qbar->id())%2==0) {
+ bornRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,true ,1);
+ EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false,1);
+ }
+ else {
+ bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,1);
+ EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false,1);
+ }
+ SpinorWaveFunction qw(p3,qbar,outgoing);
+ SpinorBarWaveFunction qbarw(p2,qbar,incoming);
+ vector<LorentzVector<Complex> > eps1,eps4;
+ SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps1,eps4,3);
+ boost::numeric::ublas::matrix<Complex>
+ bornME = boost::numeric::ublas::zero_matrix<Complex>(3,3),
+ EWME = boost::numeric::ublas::zero_matrix<Complex>(3,3);
+ for(unsigned int iq=0;iq<2;++iq) {
+ if(iq==0) {
+ qw.reset (1);
+ qbarw.reset(1);
+ }
+ else {
+ qw.reset (0);
+ qbarw.reset(0);
+ }
+ LorentzVector<complex<Energy> > current = iq==0 ?
+ qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) :
+ qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave());
+ for(unsigned int i1=0;i1<2;++i1) {
+ complex<Energy> d31 = eps1[i1].dot(p3);
+ for(unsigned int i2=0;i2<2;++i2) {
+ boost::numeric::ublas::vector<complex<Energy2> > M(5);
+ Complex d34 = eps1[i1].dot(eps4[i2]);
+ complex<Energy> d42 = eps4[i2].dot(p2);
+ // M0 in paper
+ M(0) = qw.dimensionedWave().slash(eps1[i1])
+ .slash(p4-p2).vectorCurrent(qbarw.dimensionedWave()).dot(eps4[i2]);
+ // M4 in paper
+ M(2) = current.dot(eps4[i2])*d31;
+ // M5 in paper
+ M(3) = -current.dot(eps1[i1])*d42;
+ // M1 in paper (missing factor)
+ M(1) = current.dot(p4);
+ // M6 in paper
+ M(4) = M(1)*d31*d42/GeV2;
+ // M1 final factor
+ M(1) *= d34;
+ // coefficient of different contributions
+ boost::numeric::ublas::vector<Complex> Cborn(3),CEW(3);
+ if(iq==0) {
+ axpy_prod_local(bornQQGGweights,M,Cborn);
+ axpy_prod_local(EWQQGGweights ,M,CEW );
+ }
+ else {
+ axpy_prod_local(bornRRGGweights,M,Cborn);
+ axpy_prod_local(EWRRGGweights ,M,CEW );
+ }
+ unsigned int ioff = (Cborn.size()==6 && abs(qbar->id())%2!=0) ? 3 : 0;
+ for(unsigned int ix=0;ix<3;++ix) {
+ for(unsigned int iy=0;iy<3;++iy) {
+ bornME(ix,iy) += Cborn(ix+ioff)*conj(Cborn(iy+ioff));
+ EWME (ix,iy) += CEW (ix+ioff)*conj(CEW (iy+ioff));
+ }
+ }
+ }
+ }
+ }
+ double born = 24.*real(bornME(0,0))+20./3.*real(bornME(1,1))+12.*real(bornME(2,2));
+ double EW = 24.*real(EWME(0,0))+20./3.*real(EWME(1,1))+12.*real(EWME(2,2));
+ return EW/born;
+}
diff --git a/MatrixElement/EW/ElectroWeakReweighter.h b/MatrixElement/EW/ElectroWeakReweighter.h
--- a/MatrixElement/EW/ElectroWeakReweighter.h
+++ b/MatrixElement/EW/ElectroWeakReweighter.h
@@ -1,164 +1,175 @@
// -*- C++ -*-
#ifndef Herwig_ElectroWeakReweighter_H
#define Herwig_ElectroWeakReweighter_H
//
// This is the declaration of the ElectroWeakReweighter class.
//
#include "ThePEG/MatrixElement/ReweightBase.h"
#include "EWCouplings.h"
#include "CollinearSudakov.h"
#include "SoftSudakov.h"
namespace Herwig {
using namespace ThePEG;
/**
* The ElectroWeakReweighter class.
*
* @see \ref ElectroWeakReweighterInterfaces "The interfaces"
* defined for ElectroWeakReweighter.
*/
class ElectroWeakReweighter: public ReweightBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ElectroWeakReweighter();
/**
* The destructor.
*/
virtual ~ElectroWeakReweighter();
//@}
public:
/**
* Return the weight for the kinematical configuation provided by
* the assigned XComb object (in the LastXCombInfo base class).
*/
virtual double weight() const;
/**
*
*/
static tEWCouplingsPtr coupling() {
assert(staticEWCouplings_);
return staticEWCouplings_;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Functions to reweight specific processes
*/
//@{
/**
* Reweight \f$g g\to q\bar{q}\f$
*/
double reweightggqqbar() const;
/**
* Reweight \f$q\bar{q}\to g g\f$
*/
double reweightqqbargg() const;
+
+ /**
+ * Reweight \f$q g\to qg\f$
+ */
+ double reweightqgqg() const;
+
+ /**
+ * Reweight \f$q g\to qg\f$
+ */
+ double reweightqbargqbarg() const;
//@}
protected:
/**
* Check the evolution for a fixed s,t,u
*/
void testEvolution(Energy2 s,Energy2 t, Energy2 u) const;
/**
* Evalaute the running
*/
boost::numeric::ublas::matrix<complex<InvEnergy2> >
evaluateRunning(EWProcess::Process process, Energy2 s,
- Energy2 t, Energy2 u, bool born) const;
+ Energy2 t, Energy2 u, bool born,
+ unsigned int iswap) const;
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ElectroWeakReweighter & operator=(const ElectroWeakReweighter &);
private:
/**
* The Electroweak Couplings
*/
EWCouplingsPtr EWCouplings_;
/**
* The Collinear Sudakov
*/
CollinearSudakovPtr collinearSudakov_;
/**
* The Soft Sudakov
*/
SoftSudakovPtr softSudakov_;
/**
* The couplings to allow global access
*/
static tEWCouplingsPtr staticEWCouplings_;
};
}
#endif /* Herwig_ElectroWeakReweighter_H */
diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h
--- a/MatrixElement/EW/GroupInvariants.h
+++ b/MatrixElement/EW/GroupInvariants.h
@@ -1,364 +1,407 @@
// -*- C++ -*-
//
// GroupInvariants.h is a part of Herwig - A multi-purpose Monte Carlo event generator
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
//
#ifndef HERWIG_GroupInvariants_H
#define HERWIG_GroupInvariants_H
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Config/Unitsystem.h"
#include <cassert>
#include <boost/numeric/ublas/matrix.hpp>
namespace Herwig {
using namespace ThePEG;
namespace GroupInvariants {
/**
* Simple struct for storing the different gauge contributions
*/
struct GaugeContributions {
/**
* Default Constructor
*/
GaugeContributions(double inSU3=0.,
double inSU2=0., double inU1=0.)
: SU3(inSU3),SU2(inSU2),U1(inU1)
{}
/**
* \f$SU(3)\f$
*/
double SU3;
/**
* \f$SU(2)\f$
*/
double SU2;
/**
* \f$U(1)\f$
*/
double U1;
};
/**
* The \f$SU(N)\f$ \f$C_A\f$
*/
inline double C_A(unsigned int N) {
return N !=1 ? double(N) : 0.;
}
/**
* The \f$SU(N)\f$ \f$C_F\f$
*/
inline double C_F(unsigned int N) {
return N !=1 ? 0.5*(double(N*N)-1.)/double(N) : 1.;
}
/*
* The \f$SU(N)\f$ \f$C_d\f$
*/
inline double C_d(unsigned int N) {
return (double(N*N)-4.)/double(N);
}
/**
* The \f$SU(N)\f$\f$C_1\f$
*/
inline double C_1(unsigned int N) {
double N2(N*N);
return 0.25*(N2-1.0)/N2;
}
/**
* \f$T_F\f$
*/
inline double T_F(unsigned int N, bool high) {
if(high) {
return N !=1 ? 0.5 : 5.0/3.0;
}
else {
return N !=1 ? 0.5 : 20.0/3.0;
}
}
/**
* \f$t_S\f$
*/
inline double t_S(unsigned int, bool ) {
return 0.5;
}
/**
* / Number of complex scalars in the fundamental rep. of SU(N)/U(1)
*/
inline double n_S(unsigned int N, bool high) {
if(high) {
if(N==2 || N==1) return 1.0;
else if(N==3) return 0.0;
else assert(false);
}
else {
if(N>=1&&N<=3) return 0.;
else assert(false);
}
}
/**
* Number of Dirac Fermions in the fund. rep. of SU(N) (or U(1) for N==1)
*/
inline double n_F(unsigned int N, bool high) {
if(high) {
if(N==1) return 3.0;
else if(N==2 || N==3) return 6.0;
else assert(false);
}
else {
if(N==1) return 1.0;
else if(N==2) return 0.0;
else if(N==3) return 5.0;
else assert(false);
}
}
/**
* Find K_i for gauge group N. high=false for running at mu<mZ
*/
double K_Factor(unsigned int i,unsigned int N, bool high);
/**
* Find B_i for gauge group N, high energy
*/
double B_Factor(int i, int N, bool fermion, bool longitudinal);
/**
* Find B_i for gauge group N, low energy
*/
double B_Factor_Low(int i, int N, bool fermion, double boostFactor);
/**
* Contributions to the Cusps
*/
GaugeContributions cuspContributions(Energy mu, int K_ORDER, bool high);
/**
* Contributions to B, high energy
*/
GaugeContributions BContributions(Energy mu, int B_ORDER,
bool fermion, bool longitudinal);
/**
* Contributions to B, low energy
*/
GaugeContributions BContributionsLow(Energy mu, int B_ORDER,
bool fermion, double boostFactor);
inline Complex PlusLog(double arg) {
static const Complex I(0,1.0);
if (arg>0.0)
return log(arg);
else if (arg<0.0)
return log(-arg)+I*Constants::pi;
else
assert(false);
}
inline Complex MinusLog(double arg) {
static const Complex I(0,1.0);
if (arg>0.0)
return log(arg);
else if (arg<0.0)
return log(-arg)-I*Constants::pi;
else
assert(false);
}
inline Complex getT(Energy2 s, Energy2 t) {
return MinusLog(-t/GeV2) - MinusLog(-s/GeV2);
}
inline Complex getU(Energy2 s, Energy2 u) {
return MinusLog(-u/GeV2) - MinusLog(-s/GeV2);
}
inline boost::numeric::ublas::matrix<Complex> Gamma2(Complex U, Complex T) {
boost::numeric::ublas::matrix<Complex> output(2,2);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) = (-3.0/2.0)*I*pi + (T+U);
output(1,1) = (-3.0/2.0)*I*pi;
output(0,1) = 2.0*(T-U);
output(1,0) = (3.0/8.0)*(T-U);
return output;
}
inline boost::numeric::ublas::matrix<Complex> Gamma2w(Complex U, Complex T) {
boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(5,5);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) += -I*pi*11.0/4.0;
output(0,1) += U-T;
output(1,0) += 2.0*(U-T);
output(1,1) += -I*pi*11.0/4.0 + (T+U);
output(2,2) += -7.0/4.0*I*pi + (U+T);
output(3,3) += -7.0/4.0*I*pi + (U+T);
output(4,4) += -3.0/4.0*I*pi;
return output;
}
inline boost::numeric::ublas::matrix<Complex> Gamma2Singlet() {
- boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(2,2);
+ using namespace boost::numeric::ublas;
+ matrix<Complex> output = zero_matrix<Complex>(2,2);
static const Complex I(0,1.0);
using Constants::pi;
- output(0,0) = output(1,1) = -3.0/4.0*I*pi;
+ output(0,0) = output(1,1) = -0.75*I*pi;
+ return output;
+ }
+
+ inline boost::numeric::ublas::matrix<Complex> Gamma2SingletST(Complex T) {
+ using namespace boost::numeric::ublas;
+ matrix<Complex> output = zero_matrix<Complex>(2,2);
+ static const Complex I(0,1.0);
+ using Constants::pi;
+ output(0,0) = output(1,1) = 0.75*(T-I*pi);
return output;
}
inline Complex Gamma1(double hypercharge) {
Complex I(0,1.0);
- return -I*Constants::pi*(hypercharge*hypercharge);
+ return -I*Constants::pi*sqr(hypercharge);
+ }
+
+ inline Complex Gamma1ST(double hypercharge,Complex T) {
+ Complex I(0,1.0);
+ return (T-I*Constants::pi)*sqr(hypercharge);
}
inline Complex Gamma1(double y1, double y2, Complex T, Complex U) {
Complex I(0,1.0);
return -I*Constants::pi*(y1*y1+y2*y2) + 2.0*y1*y2*(T-U);
}
+ inline Complex Gamma1ST(double y1, double y2, Complex T, Complex U) {
+ Complex I(0,1.0);
+ return (T-I*Constants::pi)*(y1*y1+y2*y2) - 2.0*y1*y2*U;
+ }
+
inline Complex Gamma1(double y1, double y2, double y3, double y4,
Complex T, Complex U) {
Complex I(0,1.0);
return -I*Constants::pi*(y1*y1+y2*y2+y3*y3+y4*y4)/2.0 +
(y1*y4+y2*y3)*T - (y1*y3+y2*y4)*U;
}
+
+ inline Complex Gamma1ST(double y1, double y2, double y3, double y4,
+ Complex T, Complex U) {
+ Complex I(0,1.0);
+ return (T-I*Constants::pi)*(y1*y1+y2*y2+y3*y3+y4*y4)/2.0 -
+ (y1*y4+y2*y3)*T - (y1*y3+y2*y4)*(U-T);
+ }
inline boost::numeric::ublas::matrix<Complex> Gamma3(Complex U, Complex T) {
boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(2,2);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) = -(8.0/3.0)*I*pi;
output(1,1) = -(8.0/3.0)*I*pi;
output(0,0) += (7.0/3.0)*T + (2.0/3.0)*U;
output(0,1) = 2.0*(T-U);
output(1,0) = (4.0/9.0)*(T-U);
return output;
}
inline boost::numeric::ublas::matrix<Complex> Gamma3g(Complex U, Complex T) {
boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(3,3);
static const Complex I(0,1.0);
using Constants::pi;
output(0,2) = U-T;
output(1,1) = 3.0/2.0*(T+U);
output(1,2) = 3.0/2.0*(U-T);
output(2,0) = 2.0*(U-T);
output(2,1) = 5.0/6.0*(U-T);
output(2,2) = 3.0/2.0*(T+U);
for (unsigned int i=0; i<3; i++) {
output(i,i) += -13.0/3.0*I*pi;
}
return output;
}
+ inline boost::numeric::ublas::matrix<Complex> Gamma3gST(Complex U, Complex T) {
+ boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(3,3);
+ static const Complex I(0,1.0);
+ using Constants::pi;
+ output(0,2) = U;
+ output(1,1) = 3.0/2.0*(U-2.*T);
+ output(1,2) = 3.0/2.0*U;
+ output(2,0) = 2.0*U;
+ output(2,1) = 5.0/6.0*U;
+ output(2,2) = 3.0/2.0*(U-2.*T);
+ for (unsigned int i=0; i<3; i++) {
+ output(i,i) += 13.0/3.0*(T-I*pi);
+ }
+ return output;
+ }
+
inline boost::numeric::ublas::matrix<Complex> Gamma3Singlet() {
boost::numeric::ublas::matrix<Complex> output = boost::numeric::ublas::zero_matrix<Complex>(2,2);
static const Complex I(0,1.0);
using Constants::pi;
output(0,0) = output(1,1) = -4.0/3.0*I*pi;
return output;
}
/**
* Number of fermion generations (only used in gauge boson HighCMatching)
*/
inline double n_g() { return 3.0; }
/**
* Number of complex scalars in the fundamental rep. of SU(N)
*/
inline double nSWeyl(unsigned int N, bool high) {
if(high) {
if(N==2 || N==1) return 1.0;
else if (N==3) return 0.0;
else assert(false);
}
else {
if( N==1 || N==3 ) return 0.0;
else assert(false);
}
}
/**
* Number of Weyl Fermions in the fundamental rep. of SU(N)
*/
inline double nFWeyl(unsigned int N, bool high) {
if(high) {
if(N==2 || N==3) return 12.0;
else assert(false);
}
else {
if(N==3) return 10.0;
else if(N==1) return 2.0;
else assert(false);
}
}
inline double TFWeyl(unsigned int) {
return 0.5;
}
inline double tSWeyl(unsigned int) {
return 0.5;
}
inline Complex WFunction(Energy mu, Energy2 s) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ln = MinusLog(-s/(mu*mu));
return (-1.0*ln*ln + 3.0*ln+pi*pi/6.0-8.0);
}
/**
* \fX_N\f% function, v is either t or u
*/
inline Complex XNFunction(unsigned int N, Energy mu, Energy2 s, Energy2 v) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ls = MinusLog(-s/(mu*mu));
return (2.0*C_F(N)*WFunction(mu,s) +
C_A(N)*(2.0*ls*ls - 2.0*MinusLog((s+v)/(mu*mu))*ls -
11.0/3.0*ls + pi*pi + 85.0/9.0) +
(2.0/3.0*ls - 10.0/9.0) * TFWeyl(N) * nFWeyl(N,true) +
(1.0/3.0*ls - 8.0/9.0) * TFWeyl(N) * nSWeyl(N,true));
}
/**
* \f$\Pi_1\f$ function
*/
inline Complex PI1_function(Energy mu, Energy2 s) {
assert(abs(s)>ZERO);
return ((41.0/6.0)*MinusLog(-s/(mu*mu))-104.0/9.0);
}
/**
* \f$\tilde{f}\f$ function, v is either t or u
*/
inline Complex fTildeFunction(Energy mu, Energy2 s, Energy2 v) {
using Constants::pi;
assert(abs(s)>ZERO);
Complex ls = MinusLog(-s/GeV2), lv = MinusLog(-v/GeV2);
Complex lsv = MinusLog((s+v)/GeV2);
return (-2.0*double(s/(s+v))*(lv-ls) +
double(s*(s+2.0*v)/((s+v)*(s+v))) * ((lv-ls)*(lv-ls) + pi*pi) +
4.0*MinusLog(-s/(mu*mu))*(lv-lsv));
}
}
}
#endif // HERWIG_GroupInvariants_H
diff --git a/MatrixElement/EW/SoftSudakov.cc b/MatrixElement/EW/SoftSudakov.cc
--- a/MatrixElement/EW/SoftSudakov.cc
+++ b/MatrixElement/EW/SoftSudakov.cc
@@ -1,1050 +1,1146 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SoftSudakov class.
//
#include "SoftSudakov.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "GroupInvariants.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "expm-1.h"
using namespace Herwig;
using namespace GroupInvariants;
SoftSudakov::SoftSudakov() : K_ORDER_(3), integrator_(0.,1e-5,1000) {}
SoftSudakov::~SoftSudakov() {}
IBPtr SoftSudakov::clone() const {
return new_ptr(*this);
}
IBPtr SoftSudakov::fullclone() const {
return new_ptr(*this);
}
void SoftSudakov::persistentOutput(PersistentOStream & os) const {
os << K_ORDER_;
}
void SoftSudakov::persistentInput(PersistentIStream & is, int) {
is >> K_ORDER_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<SoftSudakov,Interfaced>
describeHerwigSoftSudakov("Herwig::SoftSudakov", "HwMEEW.so");
void SoftSudakov::Init() {
static ClassDocumentation<SoftSudakov> documentation
("The SoftSudakov class implements the soft EW Sudakov");
}
InvEnergy SoftSudakov::operator ()(Energy mu) const {
// Include K-factor Contributions (Cusps):
GaugeContributions cusp = cuspContributions(mu,K_ORDER_,high_);
Complex gamma = cusp.SU3*G3_(row_,col_) + cusp.SU2*G2_(row_,col_) + cusp.U1*G1_(row_,col_);
if (real_) {
return gamma.real()/mu;
}
else {
return gamma.imag()/mu;
}
}
boost::numeric::ublas::matrix<Complex>
SoftSudakov::evaluateSoft(boost::numeric::ublas::matrix<Complex> & G3,
boost::numeric::ublas::matrix<Complex> & G2,
boost::numeric::ublas::matrix<Complex> & G1,
Energy mu_h, Energy mu_l, bool high) {
assert( G3.size1() == G2.size1() && G3.size1() == G1.size1() &&
G3.size2() == G2.size2() && G3.size2() == G1.size2() &&
G3.size1() == G3.size2());
G3_ = G3;
G2_ = G2;
G1_ = G1;
high_ = high;
unsigned int NN = G3_.size1();
// gamma is the matrix to be numerically integrated to run the coefficients.
boost::numeric::ublas::matrix<Complex> gamma(NN,NN);
for(row_=0;row_<NN;++row_) {
for(col_=0;col_<NN;++col_) {
if(G3_(row_,col_) == 0. && G2_(row_,col_) == 0. && G1_(row_,col_) == 0.) {
gamma(row_,col_) = 0.;
}
else {
real_ = true;
gamma(row_,col_).real(integrator_.value(*this,mu_h,mu_l));
real_ = false;
gamma(row_,col_).imag(integrator_.value(*this,mu_h,mu_l));
}
}
}
// Resummed:
return boost::numeric::ublas::expm_pad(gamma,7);
}
boost::numeric::ublas::matrix<Complex>
SoftSudakov::lowEnergyRunning(Energy EWScale, Energy lowScale,
Energy2 s, Energy2 t, Energy2 u,
- Herwig::EWProcess::Process process) {
+ Herwig::EWProcess::Process process,
+ unsigned int iswap) {
using namespace EWProcess;
+ using namespace boost::numeric::ublas;
using Constants::pi;
static const Complex I(0,1.0);
Complex T = getT(s,t), U = getU(s,u);
- boost::numeric::ublas::matrix<Complex> G1, G2, G3;
+ matrix<Complex> G1, G2, G3;
unsigned int numBrokenGauge;
switch (process) {
case QQQQ:
case QQQQiden:
case QtQtQQ:
{
+ assert(iswap==0);
numBrokenGauge = 12;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ matrix<Complex> gam3 = Gamma3(U,T);
for (unsigned int i=0; i<numBrokenGauge/2; i++) {
G3(i,i) += gam3(0,0);
G3(i,i+6) += gam3(0,1);
G3(i+6,i) += gam3(1,0);
G3(i+6,i+6) += gam3(1,1);
}
G1(0,0) = G1(6,6) = Gamma1(2.0/3.0,2.0/3.0,2.0/3.0,2.0/3.0,T,U);
G1(1,1) = G1(7,7) = Gamma1(-1.0/3.0,-1.0/3.0,2.0/3.0,2.0/3.0,T,U);
G1(2,2) = G1(8,8) = Gamma1(2.0/3.0,2.0/3.0,-1.0/3.0,-1.0/3.0,T,U);
G1(3,3) = G1(9,9) = Gamma1(-1.0/3.0,-1.0/3.0,-1.0/3.0,-1.0/3.0,T,U);
G1(4,4) = G1(10,10) = Gamma1(-1.0/3.0,2.0/3.0,2.0/3.0,-1.0/3.0,T,U);
G1(5,5) = G1(11,11) = Gamma1(2.0/3.0,-1.0/3.0,-1.0/3.0,2.0/3.0,T,U);
}
break;
case QQUU:
case QtQtUU:
case QQtRtR:
{
+ assert(iswap==0);
numBrokenGauge = 4;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ matrix<Complex> gam3 = Gamma3(U,T);
for (unsigned int i=0; i<numBrokenGauge/2; i++) {
G3(i,i) += gam3(0,0);
G3(i,i+2) += gam3(0,1);
G3(i+2,i) += gam3(1,0);
G3(i+2,i+2) += gam3(1,1);
}
G1(0,0) = G1(2,2) = Gamma1(2.0/3.0,2.0/3.0,T,U);
G1(1,1) = G1(3,3) = Gamma1(2.0/3.0,-1.0/3.0,T,U);
}
break;
case QQDD:
case QtQtDD:
{
+ assert(iswap==0);
numBrokenGauge = 4;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ matrix<Complex> gam3 = Gamma3(U,T);
for (unsigned int i=0; i<numBrokenGauge/2; i++) {
G3(i,i) += gam3(0,0);
G3(i,i+2) += gam3(0,1);
G3(i+2,i) += gam3(1,0);
G3(i+2,i+2) += gam3(1,1);
}
G1(0,0) = G1(2,2) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
G1(1,1) = G1(3,3) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
}
break;
case QQLL:
{
+ assert(iswap==0);
numBrokenGauge = 6;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(2.0/3.0,2.0/3.0,0.0,0.0,T,U);
G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,0.0,0.0,T,U);
G1(2,2) = Gamma1(2.0/3.0,2.0/3.0,-1.0,-1.0,T,U);
G1(3,3) = Gamma1(-1.0/3.0,-1.0/3.0,-1.0,-1.0,T,U);
G1(4,4) = Gamma1(-1.0/3.0,2.0/3.0,0.0,-1.0,T,U);
G1(5,5) = Gamma1(2.0/3.0,-1.0/3.0,-1.0,0.0,T,U);
}
break;
case QQEE:
{
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<2; i++) {
G3(i,i) += gam3s;
}
G1(0,0) = Gamma1(2.0/3.0,-1.0,T,U);
G1(1,1) = Gamma1(-1.0/3.0,-1.0,T,U);
}
break;
case UUUU:
case UUUUiden:
case tRtRUU:
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,2.0/3.0,T,U);
break;
case UUDD:
case tRtRDD:
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
break;
case UULL:
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3(0,0) = G3(1,1) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(2.0/3.0,0.0,T,U);
G1(1,1) = Gamma1(2.0/3.0,-1.0,T,U);
break;
case UUEE:
+ assert(iswap==0);
numBrokenGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(2.0/3.0,-1.0,T,U);
break;
case DDDD:
case DDDDiden:
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
break;
case DDLL:
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3(0,0) = G3(1,1) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/3.0,0.0,T,U);
G1(1,1) = Gamma1(-1.0/3.0,-1.0,T,U);
break;
case DDEE:
+ assert(iswap==0);
numBrokenGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/3.0,-1.0,T,U);
break;
case LLLL:
case LLLLiden:
+ assert(iswap==0);
numBrokenGauge = 6;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(0.0,0.0,0.0,0.0,T,U);
G1(1,1) = Gamma1(-1.0,-1.0,0.0,0.0,T,U);
G1(2,2) = Gamma1(0.0,0.0,-1.0,-1.0,T,U);
G1(3,3) = Gamma1(-1.0,-1.0,-1.0,-1.0,T,U);
G1(4,4) = Gamma1(-1.0,0.0,0.0,-1.0,T,U);
G1(5,5) = Gamma1(0.0,-1.0,-1.0,0.0,T,U);
break;
case LLEE:
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(0.0,-1.0,T,U);
G1(1,1) = Gamma1(-1.0,-1.0,T,U);
break;
case EEEE:
case EEEEiden:
+ assert(iswap==0);
numBrokenGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(-1.0,-1.0,T,U);
break;
case QQWW:
{
+ assert(iswap==0);
numBrokenGauge = 20;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(2./3.,2./3.,-1.,-1.,T,U);
G1(1,1) = Gamma1(2./3.,2./3.,1.,1.,T,U);
G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(4,4) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(5,5) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(6,6) = Gamma1(-1./3.,-1./3.,-1.,-1.,T,U);
G1(7,7) = Gamma1(-1./3.,-1./3.,1.,1.,T,U);
G1(8,8) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(9,9) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(10,10) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(11,11) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(12,12) = Gamma1(-1./3.,2./3.,0.,-1.,T,U);
G1(13,13) = Gamma1(-1./3.,2./3.,0.,-1.,T,U);
G1(14,14) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(15,15) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(16,16) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
G1(17,17) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
G1(18,18) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
G1(19,19) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
}
break;
case QQPhiPhi:
{
+ assert(iswap==0);
numBrokenGauge = 14;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(2./3.,2./3.,1.,1.,T,U);
G1(1,1) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(4,4) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(5,5) = Gamma1(-1./3.,-1./3.,1.,1.,T,U);
G1(6,6) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(7,7) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(8,8) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(9,9) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(10,10) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(11,11) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(12,12) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
G1(13,13) = Gamma1(2./3.,-1./3.,0.,1.,T,U);
}
break;
case QQWG:
+ assert(iswap==0);
numBrokenGauge = 6;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = -17.0/6.0*I*pi + 3.0/2.0*(U+T);
}
G1(0,0) = Gamma1(-1./3.,2./3.,1.,0.,T,U);
G1(1,1) = Gamma1(2./3.,-1./3.,-1.,0.,T,U);
G1(2,2) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(3,3) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(4,4) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(5,5) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
break;
case QQBG:
+ assert(iswap==0);
numBrokenGauge = 4;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
}
G1(0,0) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(1,1) = Gamma1(2./3.,2./3.,0.,0.,T,U);
G1(2,2) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
G1(3,3) = Gamma1(-1./3.,-1./3.,0.,0.,T,U);
break;
case QQGG:
case QtQtGG:
{
numBrokenGauge = 6;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- boost::numeric::ublas::matrix<Complex> gam3g = Gamma3g(U,T);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ matrix<Complex> gam3g;
+ Complex gam1a(0.),gam1b(0.);
+ if(iswap==0) {
+ gam3g = Gamma3g(U,T);
+ gam1a = Gamma1( 2./3.,0.,T,U);
+ gam1b = Gamma1(-1./3.,0.,T,U);
+ }
+ else if(iswap==1) {
+ gam3g = Gamma3gST(U,T);
+ gam1a = Gamma1ST( 2./3.,0.,T,U);
+ gam1b = Gamma1ST(-1./3.,0.,T,U);
+ }
+ else
+ assert(false);
for(unsigned int ix=0;ix<3;++ix) {
for(unsigned int iy=0;iy<3;++iy) {
G3(ix ,iy ) = gam3g(ix,iy);
G3(ix+3,iy+3) = gam3g(ix,iy);
}
}
- G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(2./3.,0.,T,U);
- G1(3,3) = G1(4,4) = G1(5,5) = Gamma1(-1./3.,0.,T,U);
+ G1(0,0) = G1(1,1) = G1(2,2) = gam1a;
+ G1(3,3) = G1(4,4) = G1(5,5) = gam1b;
}
break;
case LLWW:
+ assert(iswap==0);
numBrokenGauge = 20;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(0.,0.,-1.,-1.,T,U);
G1(1,1) = Gamma1(0.,0.,1.,1.,T,U);
G1(2,2) = Gamma1(0.,0.,0.,0.,T,U);
G1(3,3) = Gamma1(0.,0.,0.,0.,T,U);
G1(4,4) = Gamma1(0.,0.,0.,0.,T,U);
G1(5,5) = Gamma1(0.,0.,0.,0.,T,U);
G1(6,6) = Gamma1(-1.,-1.,-1.,-1.,T,U);
G1(7,7) = Gamma1(-1.,-1.,1.,1.,T,U);
G1(8,8) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(9,9) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(10,10) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(11,11) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(12,12) = Gamma1(-1.,0.,0.,-1.,T,U);
G1(13,13) = Gamma1(-1.,0.,0.,-1.,T,U);
G1(14,14) = Gamma1(-1.,0.,1.,0.,T,U);
G1(15,15) = Gamma1(-1.,0.,1.,0.,T,U);
G1(16,16) = Gamma1(0.,-1.,-1.,0.,T,U);
G1(17,17) = Gamma1(0.,-1.,-1.,0.,T,U);
G1(18,18) = Gamma1(0.,-1.,0.,1.,T,U);
G1(19,19) = Gamma1(0.,-1.,0.,1.,T,U);
break;
case LLPhiPhi:
+ assert(iswap==0);
numBrokenGauge = 14;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(0.,0.,1.,1.,T,U);
G1(1,1) = Gamma1(0.,0.,0.,0.,T,U);
G1(2,2) = Gamma1(0.,0.,0.,0.,T,U);
G1(3,3) = Gamma1(0.,0.,0.,0.,T,U);
G1(4,4) = Gamma1(0.,0.,0.,0.,T,U);
G1(5,5) = Gamma1(-1.,-1.,1.,1.,T,U);
G1(6,6) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(7,7) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(8,8) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(9,9) = Gamma1(-1.,-1.,0.,0.,T,U);
G1(10,10) = Gamma1(-1.,0.,1.,0.,T,U);
G1(11,11) = Gamma1(-1.,0.,1.,0.,T,U);
G1(12,12) = Gamma1(0.,-1.,0.,1.,T,U);
G1(13,13) = Gamma1(0.,-1.,0.,1.,T,U);
break;
case UUBB:
{
+ assert(iswap==0);
numBrokenGauge = 4;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
G1(i,i) = Gamma1(2./3.,0.,T,U);
}
}
break;
case UUPhiPhi:
{
+ assert(iswap==0);
numBrokenGauge = 5;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(2.0/3.0,2.0/3.0,1.,1.,T,U);
G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(2./3.,0.,T,U);
}
break;
case UUBG:
{
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam1 = Gamma1(2./3.,0.,T,U);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
G1(i,i) = gam1;
}
}
break;
+ case DDGG:
case UUGG:
case tRtRGG:
- numBrokenGauge = 3;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = Gamma3g(U,T);
- G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(2./3.,0.,T,U);
+ {
+ numBrokenGauge = 3;
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ Complex gam1(0.);
+ double Y = process==DDGG ? -1./3. : 2./3.;
+ if(iswap==0) {
+ G3 = Gamma3g(U,T);
+ gam1 = Gamma1(Y,0.,T,U);
+ }
+ else if(iswap==1) {
+ G3 = Gamma3gST(U,T);
+ gam1 = Gamma1ST(Y,0.,T,U);
+ }
+ else
+ assert(false);
+ G1(0,0) = G1(1,1) = G1(2,2) = gam1;
+ }
break;
case DDBB:
{
+ assert(iswap==0);
numBrokenGauge = 4;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
Complex gam1 = Gamma1(-1./3.,0.,T,U);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
G1(i,i) = gam1;
}
}
break;
case DDPhiPhi:
{
+ assert(iswap==0);
numBrokenGauge = 5;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = gam3s;
}
G1(0,0) = Gamma1(-1.0/3.0,-1.0/3.0,1.,1.,T,U);
G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(-1./3.,0.,T,U);
}
break;
case DDBG:
+ assert(iswap==0);
numBrokenGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G3(i,i) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
}
G1(0,0) = G1(1,1) = Gamma1(-1./3.,0.,T,U);
break;
- case DDGG:
- numBrokenGauge = 3;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = Gamma3g(U,T);
- G1(0,0) = G1(1,1) = G1(2,2) = Gamma1(-1./3.,0.,T,U);
- break;
case EEBB:
{
+ assert(iswap==0);
numBrokenGauge = 4;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
Complex gam1 = Gamma1(-1.,0.,T,U);
for (unsigned int i=0; i<numBrokenGauge; i++) {
G1(i,i) = gam1;
};
}
break;
case EEPhiPhi:
+ assert(iswap==0);
numBrokenGauge = 5;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G1 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G2 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
+ G3 = zero_matrix<Complex>(numBrokenGauge,numBrokenGauge);
G1(0,0) = Gamma1(-1.,1.,T,U);
G1(1,1) = G1(2,2) = G1(3,3) = G1(4,4) = Gamma1(-1.,0.,T,U);
break;
default:
assert(false);
}
// return the answer
if (EWScale==lowScale) {
- return boost::numeric::ublas::identity_matrix<Complex>(G1.size1());
+ return identity_matrix<Complex>(G1.size1());
}
else {
return evaluateSoft(G3,G2,G1,EWScale,lowScale,false);
}
}
boost::numeric::ublas::matrix<Complex>
SoftSudakov::highEnergyRunning(Energy highScale, Energy EWScale,
Energy2 s, Energy2 t, Energy2 u,
- Herwig::EWProcess::Process process) {
+ Herwig::EWProcess::Process process,
+ unsigned int iswap) {
using namespace EWProcess;
+ using namespace boost::numeric::ublas;
using Constants::pi;
static const Complex I(0,1.0);
Complex T = getT(s,t), U = getU(s,u);
- boost::numeric::ublas::matrix<Complex> G1,G2,G3;
+ matrix<Complex> G1,G2,G3;
unsigned int numGauge;
switch (process) {
case QQQQ:
case QQQQiden:
case QtQtQQ:
{
+ assert(iswap==0);
numGauge = 4;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- boost::numeric::ublas::matrix<Complex> gam3 = Gamma3(U,T);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
+ matrix<Complex> gam3 = Gamma3(U,T);
G3(0,0) += gam3(0,0);
G3(0,2) += gam3(0,1);
G3(2,0) += gam3(1,0);
G3(2,2) += gam3(1,1);
G3(1,1) += gam3(0,0);
G3(1,3) += gam3(0,1);
G3(3,1) += gam3(1,0);
G3(3,3) += gam3(1,1);
- boost::numeric::ublas::matrix<Complex> gam2 = Gamma2(U,T);
+ matrix<Complex> gam2 = Gamma2(U,T);
G2(0,0) += gam2(0,0);
G2(0,1) += gam2(0,1);
G2(1,0) += gam2(1,0);
G2(1,1) += gam2(1,1);
G2(2,2) += gam2(0,0);
G2(2,3) += gam2(0,1);
G2(3,2) += gam2(1,0);
G2(3,3) += gam2(1,1);
G1(0,0) = G1(1,1) = G1(2,2) = G1(3,3) = Gamma1(1.0/6.0,1.0/6.0,T,U);
}
break;
case QQUU:
case QtQtUU:
case QQtRtR:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G2 = Gamma2Singlet();
G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,1.0/6.0,T,U);
break;
case QQDD:
case QtQtDD:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G2 = Gamma2Singlet();
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,1.0/6.0,T,U);
break;
case QQLL:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3Singlet();
G2 = Gamma2(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/2.0,1.0/6.0,T,U);
break;
case QQEE:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(-1.0,1.0/6.0,T,U);
break;
case UUUU:
case UUUUiden:
case tRtRUU:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,2.0/3.0,T,U);
break;
case UUDD:
case tRtRDD:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,2.0/3.0,T,U);
break;
case UULL:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/2.0,2.0/3.0,T,U);
break;
case UUEE:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0,2.0/3.0,T,U);
break;
case DDDD:
case DDDDiden:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,T,U);
break;
case DDLL:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/2.0,-1.0/3.0,T,U);
break;
case DDEE:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0,-1.0/3.0,T,U);
break;
case LLLL:
case LLLLiden:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G2 = Gamma2(U,T);
G1(0,0) = G1(1,1) = Gamma1(-1.0/2.0,-1.0/2.0,T,U);
break;
case LLEE:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(-1.0,-1.0/2.0,T,U);
break;
case EEEE:
case EEEEiden:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G1(0,0) = Gamma1(-1.0,-1.0,T,U);
break;
case QQWW:
{
+ assert(iswap==0);
numGauge = 5;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
Complex gam3s = Gamma3Singlet()(0,0);
for (unsigned int i=0; i<5; i++) {
G3(i,i) = gam3s;
G1(i,i) = Gamma1(1.0/6.0);
}
G2 = Gamma2w(U,T);
}
break;
case QQPhiPhi:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
G3 = Gamma3Singlet();
G2 = Gamma2(U,T);
G1(0,0) = G1(1,1) = Gamma1(1.0/2.0,1.0/6.0,T,U);
break;
case QQWG:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = -17.0/6.0*I*pi + 3.0/2.0*(U+T);
G2(0,0) = -7.0/4.0*I*pi + (U+T);
G1(0,0) = Gamma1(1.0/6.0);
break;
case QQBG:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
G2(0,0) = -3.0/4.0*I*pi;
G1(0,0) = Gamma1(1.0/6.0);
break;
case QQGG:
case QtQtGG:
{
numGauge = 3;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = Gamma3g(U,T);
- Complex gam2s = Gamma2Singlet()(0,0);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ Complex gam2s,gam1;
+ if(iswap==0) {
+ G3 = Gamma3g(U,T);
+ gam2s = Gamma2Singlet()(0,0);
+ gam1 = Gamma1(1.0/6.0);
+ }
+ else if(iswap==1) {
+ G3 = Gamma3gST(U,T);
+ gam2s = Gamma2SingletST(T)(0,0);
+ gam1 = Gamma1ST(1.0/6.0,T);
+ }
+ else
+ assert(false);
for (unsigned int i=0; i<3; i++) {
G2(i,i) = gam2s;
- G1(i,i) = Gamma1(1.0/6.0);
+ G1(i,i) = gam1;
}
}
break;
case LLWW:
+ assert(iswap==0);
numGauge = 5;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
for (unsigned int i=0; i<5; i++) {
G1(i,i) = Gamma1(-1.0/2.0);
}
G2 = Gamma2w(U,T);
break;
case LLPhiPhi:
+ assert(iswap==0);
numGauge = 2;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G2 = Gamma2(U,T);
G1(0,0) = G1(1,1) = Gamma1(1.0/2.0,-1.0/2.0,T,U);
break;
case UUBB:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(2.0/3.0);
break;
case UUPhiPhi:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(1.0/2.0,2.0/3.0,T,U);
break;
case UUBG:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
G1(0,0) = Gamma1(2.0/3.0);
break;
+ case DDGG:
case UUGG:
- case tRtRGG:
- numGauge = 3;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = Gamma3g(U,T);
- for (unsigned int i=0; i<3; i++) {
- G1(i,i) = Gamma1(2.0/3.0);
+ case tRtRGG:
+ {
+ numGauge = 3;
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ double Y = process==DDGG ? -1./3. : 2./3.;
+ Complex gam1(0.);
+ if(iswap==0) {
+ G3 = Gamma3g(U,T);
+ gam1 = Gamma1(Y);
+ }
+ else if(iswap==1) {
+ G3 = Gamma3gST(U,T);
+ gam1 = Gamma1ST(Y,T);
+ }
+ else
+ assert(false);
+ for (unsigned int i=0; i<3; i++) {
+ G1(i,i) = gam1;
+ }
}
break;
case DDBB:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G1(0,0) = Gamma1(-1.0/3.0);
break;
case DDPhiPhi:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = Gamma3Singlet()(0,0);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(1.0/2.0,-1.0/3.0,T,U);
break;
case DDBG:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G3(0,0) = -4.0/3.0*I*pi + 3.0/2.0*(U+T-I*pi);
G1(0,0) = Gamma1(-1.0/3.0);
break;
- case DDGG:
- numGauge = 3;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = Gamma3g(U,T);
- for (unsigned int i=0; i<3; i++) {
- G1(i,i) = Gamma1(-1.0/3.0);
- }
- break;
case EEBB:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G1(0,0) = Gamma1(-1.0);
break;
case EEPhiPhi:
+ assert(iswap==0);
numGauge = 1;
- G1 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G2 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
- G3 = boost::numeric::ublas::zero_matrix<Complex>(numGauge,numGauge);
+ G1 = zero_matrix<Complex>(numGauge,numGauge);
+ G2 = zero_matrix<Complex>(numGauge,numGauge);
+ G3 = zero_matrix<Complex>(numGauge,numGauge);
G2(0,0) = Gamma2Singlet()(0,0);
G1(0,0) = Gamma1(1.0/2.0,-1.0,T,U);
break;
default:
assert(false);
}
return evaluateSoft(G3,G2,G1,highScale,EWScale,true);
}
unsigned int SoftSudakov::numberGauge(Herwig::EWProcess::Process process) {
using namespace EWProcess;
switch (process) {
case QQQQ:
case QQQQiden:
case QtQtQQ:
return 4;
case QQUU:
case QtQtUU:
case QQtRtR:
return 2;
case QQDD:
case QtQtDD:
return 2;
case QQLL:
return 2;
case QQEE:
return 1;
case UUUU:
case UUUUiden:
case tRtRUU:
return 2;
case UUDD:
case tRtRDD:
return 2;
case UULL:
return 1;
case UUEE:
return 1;
case DDDD:
case DDDDiden:
return 2;
case DDLL:
return 1;
case DDEE:
return 1;
case LLLL:
case LLLLiden:
return 2;
case LLEE:
return 1;
case EEEE:
case EEEEiden:
return 1;
case QQWW:
return 5;
case QQPhiPhi:
return 2;
case QQWG:
return 1;
case QQBG:
return 1;
case QQGG:
case QtQtGG:
return 3;
case LLWW:
return 5;
case LLPhiPhi:
return 2;
case UUBB:
return 1;
case UUPhiPhi:
return 1;
case UUBG:
return 1;
case UUGG:
case tRtRGG:
return 3;
case DDBB:
return 1;
case DDPhiPhi:
return 1;
case DDBG:
return 1;
case DDGG:
return 3;
case EEBB:
return 1;
case EEPhiPhi:
return 1;
default:
assert(false);
}
}
unsigned int SoftSudakov::numberBrokenGauge(Herwig::EWProcess::Process process) {
using namespace EWProcess;
switch (process) {
case QQQQ:
case QQQQiden:
case QtQtQQ:
return 12;
case QQUU:
case QtQtUU:
case QQtRtR:
return 4;
case QQDD:
case QtQtDD:
return 4;
case QQLL:
return 6;
case QQEE:
return 2;
case UUUU:
case UUUUiden:
case tRtRUU:
return 2;
case UUDD:
case tRtRDD:
return 2;
case UULL:
return 2;
case UUEE:
return 1;
case DDDD:
case DDDDiden:
return 2;
case DDLL:
return 2;
case DDEE:
return 1;
case LLLL:
case LLLLiden:
return 6;
case EEEE:
case EEEEiden:
return 1;
case QQWW:
return 20;
case QQPhiPhi:
return 14;
case QQWG:
return 6;
case QQBG:
return 4;
case QQGG:
case QtQtGG:
return 6;
case LLWW:
return 20;
case LLPhiPhi:
return 14;
case UUBB:
return 4;
case UUPhiPhi:
return 5;
case UUBG:
return 2;
case UUGG:
case tRtRGG:
return 3;
case DDBB:
return 4;
case DDPhiPhi:
return 5;
case DDBG:
return 2;
case DDGG:
return 3;
case EEBB:
return 4;
case EEPhiPhi:
return 5;
default:
assert(false);
}
}
diff --git a/MatrixElement/EW/SoftSudakov.h b/MatrixElement/EW/SoftSudakov.h
--- a/MatrixElement/EW/SoftSudakov.h
+++ b/MatrixElement/EW/SoftSudakov.h
@@ -1,191 +1,193 @@
// -*- C++ -*-
#ifndef Herwig_SoftSudakov_H
#define Herwig_SoftSudakov_H
//
// This is the declaration of the SoftSudakov class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/Utilities/GSLIntegrator.h"
#include <boost/numeric/ublas/matrix.hpp>
#include "EWProcess.h"
#include "SoftSudakov.fh"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the SoftSudakov class.
*
* @see \ref SoftSudakovInterfaces "The interfaces"
* defined for SoftSudakov.
*/
class SoftSudakov: public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SoftSudakov();
/**
* The destructor.
*/
virtual ~SoftSudakov();
//@}
public:
/**
* Low energy soft evolution
*/
boost::numeric::ublas::matrix<Complex>
lowEnergyRunning(Energy EWScale, Energy lowScale,
Energy2 s, Energy2 t, Energy2 u,
- Herwig::EWProcess::Process process);
-
+ Herwig::EWProcess::Process process,
+ unsigned int iswap);
+
/**
* Evalaute the high energy running as a matrix
*/
boost::numeric::ublas::matrix<Complex>
highEnergyRunning(Energy highScale, Energy EWScale,
Energy2 s, Energy2 t, Energy2 u,
- Herwig::EWProcess::Process process);
+ Herwig::EWProcess::Process process,
+ unsigned int iswap);
/**
* Number of operators for the broken theory at low energy
*/
unsigned int numberBrokenGauge(Herwig::EWProcess::Process process);
/**
* Number of operators for the unbroken theory at high energy
*/
unsigned int numberGauge(Herwig::EWProcess::Process process);
protected:
/**
* Evaluate the soft evolution
*/
boost::numeric::ublas::matrix<Complex> evaluateSoft(boost::numeric::ublas::matrix<Complex> & G3,
boost::numeric::ublas::matrix<Complex> & G2,
boost::numeric::ublas::matrix<Complex> & G1,
Energy mu_h, Energy mu_l, bool high);
public:
/**
* The operator to be integrated
*/
InvEnergy operator ()(Energy mu) const;
/** Argument type for GaussianIntegrator */
typedef Energy ArgType;
/** Return type for GaussianIntegrator */
typedef InvEnergy ValType;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SoftSudakov & operator=(const SoftSudakov &);
private:
/**
* Order for K
*/
unsigned int K_ORDER_;
/**
* Integrator
*/
GSLIntegrator integrator_;
/**
* Whether doing the high or low scale contribution
*/
bool high_;
/**
* Column
*/
unsigned int row_;
/**
* Row
*/
unsigned int col_;
/**
*
*/
bool real_;
/**
*
*/
boost::numeric::ublas::matrix<Complex> G1_;
/**
*
*/
boost::numeric::ublas::matrix<Complex> G2_;
/**
*
*/
boost::numeric::ublas::matrix<Complex> G3_;
};
}
#endif /* Herwig_SoftSudakov_H */
diff --git a/MatrixElement/Hadron/MEQCD2to2.cc b/MatrixElement/Hadron/MEQCD2to2.cc
--- a/MatrixElement/Hadron/MEQCD2to2.cc
+++ b/MatrixElement/Hadron/MEQCD2to2.cc
@@ -1,1179 +1,1181 @@
// -*- C++ -*-
//
// MEQCD2to2.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEQCD2to2 class.
//
#include "MEQCD2to2.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Cuts/Cuts.h"
#include "Herwig/MatrixElement/HardVertex.h"
using namespace Herwig;MEQCD2to2::MEQCD2to2():_maxflavour(5),_process(0) {
massOption(vector<unsigned int>(2,0));
}
void MEQCD2to2::rebind(const TranslationMap & trans)
{
_ggggvertex = trans.translate(_ggggvertex);
_gggvertex = trans.translate( _gggvertex);
_qqgvertex = trans.translate( _qqgvertex);
_gluon = trans.translate( _gluon);
for(unsigned int ix=0;ix<_quark.size();++ix)
{_quark[ix]=trans.translate(_quark[ix]);}
for(unsigned int ix=0;ix<_antiquark.size();++ix)
{_antiquark[ix]=trans.translate(_quark[ix]);}
HwMEBase::rebind(trans);
}
IVector MEQCD2to2::getReferences() {
IVector ret = HwMEBase::getReferences();
ret.push_back(_ggggvertex);
ret.push_back(_gggvertex);
ret.push_back(_qqgvertex);
ret.push_back(_gluon);
for(unsigned int ix=0;ix<_quark.size();++ix)
{ret.push_back(_quark[ix]);}
for(unsigned int ix=0;ix<_antiquark.size();++ix)
{ret.push_back(_antiquark[ix]);}
return ret;
}
void MEQCD2to2::doinit() {
// call the base class
HwMEBase::doinit();
// get the vedrtex pointers from the SM object
tcHwSMPtr hwsm= dynamic_ptr_cast<tcHwSMPtr>(standardModel());
// do the initialisation
if(hwsm) {
_qqgvertex = hwsm->vertexFFG();
_gggvertex = hwsm->vertexGGG();
_ggggvertex = hwsm->vertexGGGG();
}
else throw InitException() << "Wrong type of StandardModel object in "
<< "MEQCD2to2::doinit() the Herwig version must be used"
<< Exception::runerror;
// get the particle data objects
_gluon=getParticleData(ParticleID::g);
for(int ix=1;ix<=int(_maxflavour);++ix) {
_quark.push_back( getParticleData( ix));
_antiquark.push_back(getParticleData(-ix));
}
}
Energy2 MEQCD2to2::scale() const {
Energy2 s(sHat()),u(uHat()),t(tHat());
return 2.*s*t*u/(s*s+t*t+u*u);
}
void MEQCD2to2::persistentOutput(PersistentOStream & os) const {
os << _ggggvertex << _gggvertex << _qqgvertex << _maxflavour
<< _process << _gluon << _quark << _antiquark;
}
void MEQCD2to2::persistentInput(PersistentIStream & is, int) {
is >> _ggggvertex >> _gggvertex >> _qqgvertex >> _maxflavour
>> _process >> _gluon >> _quark >> _antiquark;
}
unsigned int MEQCD2to2::orderInAlphaS() const {
return 2;
}
unsigned int MEQCD2to2::orderInAlphaEW() const {
return 0;
}
ClassDescription<MEQCD2to2> MEQCD2to2::initMEQCD2to2;
// Definition of the static class description member.
void MEQCD2to2::Init() {
static ClassDocumentation<MEQCD2to2> documentation
("The MEQCD2to2 class implements the QCD 2->2 processes in hadron-hadron"
" collisions");
static Parameter<MEQCD2to2,unsigned int> interfaceMaximumFlavour
("MaximumFlavour",
"The maximum flavour of the quarks in the process",
&MEQCD2to2::_maxflavour, 5, 1, 5,
false, false, Interface::limited);
static Switch<MEQCD2to2,unsigned int> interfaceProcess
("Process",
"Which subprocesses to include",
&MEQCD2to2::_process, 0, false, false);
static SwitchOption interfaceProcessAll
(interfaceProcess,
"All",
"Include all subprocesses",
0);
static SwitchOption interfaceProcess1
(interfaceProcess,
"gg2gg",
"Include only gg->gg subprocesses",
1);
static SwitchOption interfaceProcess2
(interfaceProcess,
"gg2qqbar",
"Include only gg -> q qbar processes",
2);
static SwitchOption interfaceProcessqqbargg
(interfaceProcess,
"qqbar2gg",
"Include only q qbar -> gg processes",
3);
static SwitchOption interfaceProcessqgqg
(interfaceProcess,
"qg2qg",
"Include only q g -> q g processes",
4);
static SwitchOption interfaceProcessqbargqbarg
(interfaceProcess,
"qbarg2qbarg",
"Include only qbar g -> qbar g processes",
5);
static SwitchOption interfaceProcessqqqq
(interfaceProcess,
"qq2qq",
"Include only q q -> q q processes",
6);
static SwitchOption interfaceProcessqbarqbarqbarqbar
(interfaceProcess,
"qbarqbar2qbarqbar",
"Include only qbar qbar -> qbar qbar processes",
7);
static SwitchOption interfaceProcessqqbarqqbar
(interfaceProcess,
"qqbar2qqbar",
"Include only q qbar -> q qbar processes",
8);
}
Selector<MEBase::DiagramIndex>
MEQCD2to2::diagrams(const DiagramVector & diags) const {
// select the diagram, this is easy for us as we have already done it
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
if(diags[i]->id()==-int(_diagram)) sel.insert(1.0, i);
else sel.insert(0., i);
}
return sel;
}
double MEQCD2to2::gg2qqbarME(vector<VectorWaveFunction> &g1,
vector<VectorWaveFunction> &g2,
vector<SpinorBarWaveFunction> & q,
vector<SpinorWaveFunction> & qbar,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
interv=_gggvertex->evaluate(mt,5,_gluon,g1[ihel1],g2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
//first t-channel diagram
inters =_qqgvertex->evaluate(mt,5,qbar[ohel2].particle(),
qbar[ohel2],g2[ihel2]);
diag[0]=_qqgvertex->evaluate(mt,inters,q[ohel1],g1[ihel1]);
//second t-channel diagram
inters =_qqgvertex->evaluate(mt,5,qbar[ohel2].particle(),
qbar[ohel2],g1[ihel1]);
diag[1]=_qqgvertex->evaluate(mt,inters,q[ohel1],g2[ihel2]);
// s-channel diagram
diag[2]=_qqgvertex->evaluate(mt,qbar[ohel2],q[ohel1],interv);
// colour flows
flow[0]=diag[0]+diag[2];
flow[1]=diag[1]-diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(2*ihel1,2*ihel2,ohel1,ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 48.*(1./6./u/t-3./8./s/s)*(t*t+u*u)*sqr(alphas)/output << endl;
// select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=4+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output/48.;
}
double MEQCD2to2::qqbar2ggME(vector<SpinorWaveFunction> & q,
vector<SpinorBarWaveFunction> & qbar,
vector<VectorWaveFunction> &g1,
vector<VectorWaveFunction> &g2,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
interv=_qqgvertex->evaluate(mt,5,_gluon,q[ihel1],qbar[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first t-channel diagram
inters=_qqgvertex->evaluate(mt,5,q[ihel1].particle()->CC(),
q[ihel1],g1[ohel1]);
diag[0]=_qqgvertex->evaluate(mt,inters,qbar[ihel2],g2[ohel2]);
// second t-channel diagram
inters=_qqgvertex->evaluate(mt,5,q[ihel1].particle()->CC(),
q[ihel1],g2[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,inters,qbar[ihel2],g1[ohel1]);
// s-channel diagram
diag[2]=_gggvertex->evaluate(mt,g1[ohel1],g2[ohel2],interv);
// colour flows
flow[0]=diag[0]-diag[2];
flow[1]=diag[1]+diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,2*ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 27./2.*0.5*(32./27./u/t-8./3./s/s)*(t*t+u*u)*sqr(alphas)/output << endl;
//select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=7+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return 2.*output/27.;
}
double MEQCD2to2::qg2qgME(vector<SpinorWaveFunction> & qin,
vector<VectorWaveFunction> &g2,
vector<SpinorBarWaveFunction> & qout,
vector<VectorWaveFunction> &g4,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorWaveFunction inters,inters2;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
inters=_qqgvertex->evaluate(mt,5,qin[ihel1].particle()->CC(),
qin[ihel1],g2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// s-channel diagram
diag[0]=_qqgvertex->evaluate(mt,inters,qout[ohel1],g4[ohel2]);
// first t-channel
inters2=_qqgvertex->evaluate(mt,5,qin[ihel1].particle()->CC(),
qin[ihel1],g4[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,inters2,qout[ohel1],g2[ihel2]);
// second t-channel
interv=_qqgvertex->evaluate(mt,5,_gluon,qin[ihel1],qout[ohel1]);
diag[2]=_gggvertex->evaluate(mt,g2[ihel2],g4[ohel2],interv);
// colour flows
flow[0]=diag[0]-diag[2];
flow[1]=diag[1]+diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,2*ihel2,ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
- //double alphas(4.*pi*SM().alphaS(mt));
+ // double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 18./output*(-4./9./s/u+1./t/t)*(s*s+u*u)*sqr(alphas) << endl;
//select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=10+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output/18.;
}
double MEQCD2to2::gg2ggME(vector<VectorWaveFunction> &g1,vector<VectorWaveFunction> &g2,
vector<VectorWaveFunction> &g3,vector<VectorWaveFunction> &g4,
unsigned int iflow) const {
// colour factors for different flows
static const double c1 = 4.*( sqr(9.)/8.-3.*9./8.+1.-0.75/9.);
static const double c2 = 4.*(-0.25*9. +1.-0.75/9.);
// scale
Energy2 mt(scale());
// // matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1,PDT::Spin1,
PDT::Spin1,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[3]={0.,0.,0.};
Complex diag[3],flow[3];
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// s-channel diagram
diag[0]=_ggggvertex->evaluate(mt,1,g3[ohel1],g1[ihel1],
g4[ohel2],g2[ihel2]);
// t-channel
diag[1]=_ggggvertex->evaluate(mt,1,g1[ihel1],g2[ihel2],
g3[ohel1],g4[ohel2]);
// u-channel
diag[2]=_ggggvertex->evaluate(mt,1,g2[ihel2],g1[ihel1],
g3[ohel1],g4[ohel2]);
// colour flows
flow[0] = diag[0]-diag[2];
flow[1] = -diag[0]-diag[1];
flow[2] = diag[1]+diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) {
sumdiag[ix] += norm(diag[ix]);
sumflow[ix] += norm(flow[ix]);
}
// total
output += c1*(norm(flow[0])+norm(flow[1])+norm(flow[2]))
+2.*c2*real(flow[0]*conj(flow[1])+flow[0]*conj(flow[2])+
flow[1]*conj(flow[2]));
// store the me if needed
if(iflow!=0) _me(2*ihel1,2*ihel2,2*ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// spin, colour and identical particle factorsxs
output /= 4.*64.*2.;
// test code vs me from ESW
// Energy2 u(uHat()),t(tHat()),s(sHat());
// using Constants::pi;
// double alphas(4.*pi*SM().alphaS(mt));
// cerr << "testing matrix element "
// << 1./output*9./4.*(3.-t*u/s/s-s*u/t/t-s*t/u/u)*sqr(alphas) << endl;
// select a colour flow
_flow=1+UseRandom::rnd3(sumflow[0],sumflow[1],sumflow[2]);
// and diagram
if(_flow==1) _diagram = 1+2*UseRandom::rnd2(sumdiag[0],sumdiag[2]);
else if(_flow==2) _diagram = 1+ UseRandom::rnd2(sumdiag[0],sumdiag[1]);
else if(_flow==3) _diagram = 2+ UseRandom::rnd2(sumdiag[1],sumdiag[2]);
// final part of colour and spin factors
return output;
}
double MEQCD2to2::qbarg2qbargME(vector<SpinorBarWaveFunction> & qin,
vector<VectorWaveFunction> &g2,
vector<SpinorWaveFunction> & qout,
vector<VectorWaveFunction> &g4,
unsigned int iflow) const {
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1,
PDT::Spin1Half,PDT::Spin1));
// calculate the matrix element
double output(0.),sumdiag[3]={0.,0.,0.},sumflow[2]={0.,0.};
Complex diag[3],flow[2];
VectorWaveFunction interv;
SpinorBarWaveFunction inters,inters2;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
inters=_qqgvertex->evaluate(mt,5,qin[ihel1].particle()->CC(),
qin[ihel1],g2[ihel2]);
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// s-channel diagram
diag[0]=_qqgvertex->evaluate(mt,qout[ohel1],inters,g4[ohel2]);
// first t-channel
inters2=_qqgvertex->evaluate(mt,5,qin[ihel1].particle()->CC(),
qin[ihel1],g4[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,qout[ohel1],inters2,g2[ihel2]);
// second t-channel
interv=_qqgvertex->evaluate(mt,5,_gluon,qout[ohel1],qin[ihel1]);
diag[2]=_gggvertex->evaluate(mt,g2[ihel2],g4[ohel2],interv);
// colour flows
flow[0]=diag[0]+diag[2];
flow[1]=diag[1]-diag[2];
// sums
for(unsigned int ix=0;ix<3;++ix) sumdiag[ix] += norm(diag[ix]);
for(unsigned int ix=0;ix<2;++ix) sumflow[ix] += norm(flow[ix]);
// total
output +=real(flow[0]*conj(flow[0])+flow[1]*conj(flow[1])
-0.25*flow[0]*conj(flow[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,2*ihel2,ohel1,2*ohel2)=flow[iflow-1];
}
}
}
}
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
- //double alphas(4.*pi*SM().alphaS(mt));
+ using Constants::pi;
+ double alphas(4.*pi*SM().alphaS(mt));
//cerr << "testing matrix element "
// << 18./output*(-4./9./s/u+1./t/t)*(s*s+u*u)*sqr(alphas) << endl;
//select a colour flow
_flow=1+UseRandom::rnd2(sumflow[0],sumflow[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=13+UseRandom::rnd3(sumdiag[0],sumdiag[1],sumdiag[2]);
+ cerr << "testing ME " << output/18./sqr(alphas) << "\n";
// final part of colour and spin factors
return output/18.;
}
double MEQCD2to2::qq2qqME(vector<SpinorWaveFunction> & q1,
vector<SpinorWaveFunction> & q2,
vector<SpinorBarWaveFunction> & q3,
vector<SpinorBarWaveFunction> & q4,
unsigned int iflow) const {
// identify special case of identical quarks
bool identical(q1[0].id()==q2[0].id());
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q3[ohel1]);
diag[0] = _qqgvertex->evaluate(mt,q2[ihel2],q4[ohel2],interv);
// second diagram if identical
if(identical) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q4[ohel2]);
diag[1]=_qqgvertex->evaluate(mt,q2[ihel2],q3[ohel1],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];
}
}
}
}
// identical particle symmetry factor if needed
if(identical) output*=0.5;
// test code vs me from ESW
//Energy2 u(uHat()),t(tHat()),s(sHat());
//double alphas(4.*pi*SM().alphaS(mt));
//if(identical)
// {cerr << "testing matrix element A "
// << 18./output*0.5*(4./9.*((s*s+u*u)/t/t+(s*s+t*t)/u/u)
// -8./27.*s*s/u/t)*sqr(alphas) << endl;}
//else
// {cerr << "testing matrix element B "
// << 18./output*(4./9.*(s*s+u*u)/t/t)*sqr(alphas) << endl;}
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=16+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
double MEQCD2to2::qbarqbar2qbarqbarME(vector<SpinorBarWaveFunction> & q1,
vector<SpinorBarWaveFunction> & q2,
vector<SpinorWaveFunction> & q3,
vector<SpinorWaveFunction> & q4,
unsigned int iflow) const {
// identify special case of identical quarks
bool identical(q1[0].id()==q2[0].id());
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0)
{_me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));}
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
interv = _qqgvertex->evaluate(mt,5,_gluon,q3[ohel1],q1[ihel1]);
diag[0] = _qqgvertex->evaluate(mt,q4[ohel2],q2[ihel2],interv);
// second diagram if identical
if(identical) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q4[ohel2],q1[ihel1]);
diag[1]=_qqgvertex->evaluate(mt,q3[ohel1],q2[ihel2],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0) _me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];
}
}
}
}
// identical particle symmetry factor if needed
if(identical){output*=0.5;}
// test code vs me from ESW
// Energy2 u(uHat()),t(tHat()),s(sHat());
// double alphas(4.*pi*SM().alphaS(mt));
// if(identical)
// {cerr << "testing matrix element A "
// << 18./output*0.5*(4./9.*((s*s+u*u)/t/t+(s*s+t*t)/u/u)
// -8./27.*s*s/u/t)*sqr(alphas) << endl;}
// else
// {cerr << "testing matrix element B "
// << 18./output*(4./9.*(s*s+u*u)/t/t)*sqr(alphas) << endl;}
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=18+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
double MEQCD2to2::qqbar2qqbarME(vector<SpinorWaveFunction> & q1,
vector<SpinorBarWaveFunction> & q2,
vector<SpinorBarWaveFunction> & q3,
vector<SpinorWaveFunction> & q4,
unsigned int iflow) const {
// type of process
bool diagon[2]={q1[0].id()== -q2[0].id(),
q1[0].id()== -q3[0].id()};
// scale
Energy2 mt(scale());
// matrix element to be stored
if(iflow!=0) _me.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));
// calculate the matrix element
double output(0.),sumdiag[2]={0.,0.};
Complex diag[2];
VectorWaveFunction interv;
for(unsigned int ihel1=0;ihel1<2;++ihel1) {
for(unsigned int ihel2=0;ihel2<2;++ihel2) {
for(unsigned int ohel1=0;ohel1<2;++ohel1) {
for(unsigned int ohel2=0;ohel2<2;++ohel2) {
// first diagram
if(diagon[0]) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q2[ihel2]);
diag[0] = _qqgvertex->evaluate(mt,q4[ohel2],q3[ohel1],interv);
}
else diag[0]=0.;
// second diagram
if(diagon[1]) {
interv = _qqgvertex->evaluate(mt,5,_gluon,q1[ihel1],q3[ohel1]);
diag[1]=_qqgvertex->evaluate(mt,q4[ohel2],q2[ihel2],interv);
}
else diag[1]=0.;
// sum of diagrams
for(unsigned int ix=0;ix<2;++ix) sumdiag[ix] += norm(diag[ix]);
// total
output +=real(diag[0]*conj(diag[0])+diag[1]*conj(diag[1])
+2./3.*diag[0]*conj(diag[1]));
// store the me if needed
if(iflow!=0){_me(ihel1,ihel2,ohel1,ohel2)=diag[iflow-1];}
}
}
}
}
// test code vs me from ESW
// Energy2 u(uHat()),t(tHat()),s(sHat());
// double alphas(4.*pi*SM().alphaS(mt));
// if(diagon[0]&&diagon[1]) {
// cerr << "testing matrix element A "
// << q1[0].id() << " " << q2[0].id() << " -> "
// << q3[0].id() << " " << q4[0].id() << " "
// << 18./output*0.5*(4./9.*((s*s+u*u)/t/t+(u*u+t*t)/s/s)
// -8./27.*u*u/s/t)*sqr(alphas) << endl;
// }
// else if(diagon[0]) {
// cerr << "testing matrix element B "
// << q1[0].id() << " " << q2[0].id() << " -> "
// << q3[0].id() << " " << q4[0].id() << " "
// << 18./output*(4./9.*(t*t+u*u)/s/s)*sqr(alphas) << endl;
// }
// else if(diagon[1]) {
// cerr << "testing matrix element C "
// << q1[0].id() << " " << q2[0].id() << " -> "
// << q3[0].id() << " " << q4[0].id() << " "
// << 18./output*(4./9.*(s*s+u*u)/t/t)*sqr(alphas) << endl;
// }
//select a colour flow
_flow=1+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// select a diagram ensuring it is one of those in the selected colour flow
sumdiag[_flow%2]=0.;
_diagram=20+UseRandom::rnd2(sumdiag[0],sumdiag[1]);
// final part of colour and spin factors
return output/18.;
}
void MEQCD2to2::getDiagrams() const {
// gg-> gg subprocess
if(_process==0||_process==1) {
// s-channel
add(new_ptr((Tree2toNDiagram(2),_gluon,_gluon, 1, _gluon,
3,_gluon, 3, _gluon, -1)));
// first t-channel
add(new_ptr((Tree2toNDiagram(3),_gluon,_gluon,_gluon,
1,_gluon, 2,_gluon,-2)));
// second t-channel
add(new_ptr((Tree2toNDiagram(3),_gluon,_gluon,_gluon,
2,_gluon, 1,_gluon,-3)));
}
// processes involving one quark line
for(unsigned int ix=0;ix<_maxflavour;++ix) {
// gg -> q qbar subprocesses
if(_process==0||_process==2) {
// first t-channel
add(new_ptr((Tree2toNDiagram(3),_gluon,_antiquark[ix],_gluon,
1,_quark[ix], 2,_antiquark[ix],-4)));
// interchange
add(new_ptr((Tree2toNDiagram(3),_gluon,_antiquark[ix],_gluon,
2,_quark[ix], 1,_antiquark[ix],-5)));
// s-channel
add(new_ptr((Tree2toNDiagram(2),_gluon,_gluon, 1, _gluon,
3,_quark[ix], 3, _antiquark[ix], -6)));
}
// q qbar -> g g subprocesses
if(_process==0||_process==3) {
// first t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_antiquark[ix],_antiquark[ix],
1,_gluon, 2,_gluon,-7)));
// second t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_antiquark[ix],_antiquark[ix],
2,_gluon, 1,_gluon,-8)));
// s-channel
add(new_ptr((Tree2toNDiagram(2),_quark[ix], _antiquark[ix],
1, _gluon, 3, _gluon, 3, _gluon,-9)));
}
// q g -> q g subprocesses
if(_process==0||_process==4) {
// s-channel
add(new_ptr((Tree2toNDiagram(2),_quark[ix], _gluon,
1, _quark[ix], 3, _quark[ix], 3, _gluon,-10)));
// quark t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_quark[ix],_gluon,
2,_quark[ix],1,_gluon,-11)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_gluon,_gluon,
1,_quark[ix],2,_gluon,-12)));
}
// qbar g -> qbar g subprocesses
if(_process==0||_process==5) {
// s-channel
add(new_ptr((Tree2toNDiagram(2),_antiquark[ix], _gluon,
1, _antiquark[ix], 3, _antiquark[ix], 3, _gluon,-13)));
// quark t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[ix],_antiquark[ix],_gluon,
2,_antiquark[ix],1,_gluon,-14)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[ix],_gluon,_gluon,
1,_antiquark[ix],2,_gluon,-15)));
}
// processes involving two quark lines
for(unsigned int iy=ix;iy<_maxflavour;++iy) {
// q q -> q q subprocesses
if(_process==0||_process==6) {
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_gluon,_quark[iy],
1,_quark[ix],2,_quark[iy],-16)));
// exchange for identical quarks
if(ix==iy)
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_gluon,_quark[iy],
2,_quark[ix],1,_quark[iy],-17)));
}
// qbar qbar -> qbar qbar subprocesses
if(_process==0||_process==7) {
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_antiquark[ix],_gluon,_antiquark[iy],
1,_antiquark[ix],2,_antiquark[iy],-18)));
// exchange for identical quarks
if(ix==iy)
add(new_ptr((Tree2toNDiagram(3),_antiquark[ix],_gluon,_antiquark[iy],
2,_antiquark[ix],1,_antiquark[iy],-19)));
}
}
for(unsigned int iy=0;iy<_maxflavour;++iy) {
// q qbar -> q qbar
if(_process==0||_process==8) {
// gluon s-channel
add(new_ptr((Tree2toNDiagram(2),_quark[ix], _antiquark[ix],
1, _gluon, 3, _quark[iy], 3, _antiquark[iy],-20)));
// gluon t-channel
add(new_ptr((Tree2toNDiagram(3),_quark[ix],_gluon,_antiquark[iy],
1,_quark[ix],2,_antiquark[iy],-21)));
}
}
}
}
Selector<const ColourLines *>
MEQCD2to2::colourGeometries(tcDiagPtr diag) const {
// colour lines for gg to gg
static const ColourLines cgggg[12]={ColourLines("1 -2, -1 -3 -5, 5 -4, 2 3 4"),// A_2 s
ColourLines("-1 2, 1 3 5, -5 4, -2 -3 -4"),// A_1 s
ColourLines("1 5, -1 -2 3, -3 -4, -5 2 4"),// A_1 u
ColourLines("-1 -5, 1 2 -3, 3 4, 5 -2 -4"),// A_2 u
ColourLines("1 -2, -1 -3 -4, 4 -5, 2 3 5"),// B_2 s
ColourLines("-1 2, 1 3 4, -4 5, -2 -3 -5"),// B_1 s
ColourLines("1 4, -1 -2 3, -3 -5, -4 2 5"),// B_1 t
ColourLines("-1 -4, 1 2 -3, 3 5, 4 -2 -5"),// B_2 t
ColourLines("1 4, -1 -2 -5, 3 5, -3 2 -4"),// C_1 t
ColourLines("-1 -4, 1 2 5, -3 -5, 3 -2 4"),// C_2 t
ColourLines("1 5, -1 -2 -4, 3 4, -3 2 -5"),// C_1 u
ColourLines("-1 -5, 1 2 4, -3 -4, 3 -2 5") // C_2 u
};
// colour lines for gg to q qbar
static const ColourLines cggqq[4]={ColourLines("1 4, -1 -2 3, -3 -5"),
ColourLines("3 4, -3 -2 1, -1 -5"),
ColourLines("2 -1, 1 3 4, -2 -3 -5"),
ColourLines("1 -2, -1 -3 -5, 2 3 4")};
// colour lines for q qbar to gg
static const ColourLines cqqgg[4]={ColourLines("1 4, -4 -2 5, -3 -5"),
ColourLines("1 5, -3 -4, 4 -2 -5"),
ColourLines("1 3 4, -4 5, -2 -3 -5"),
ColourLines("1 3 5, -5 4, -2 -3 -4")};
// colour lines for q g to q g
static const ColourLines cqgqg[4]={ColourLines("1 -2, 2 3 5, 4 -5"),
ColourLines("1 5, 3 4,-3 2 -5 "),
ColourLines("1 2 -3, 3 5, -5 -2 4"),
ColourLines("1 -2 5,3 2 4,-3 -5")};
// colour lines for qbar g -> qbar g
static const ColourLines cqbgqbg[4]={ColourLines("-1 2, -2 -3 -5, -4 5"),
ColourLines("-1 -5, -3 -4, 3 -2 5"),
ColourLines("-1 -2 3, -3 -5, 5 2 -4"),
ColourLines("-1 2 -5,-3 -2 -4, 3 5")};
// colour lines for q q -> q q
static const ColourLines cqqqq[2]={ColourLines("1 2 5,3 -2 4"),
ColourLines("1 2 4,3 -2 5")};
// colour lines for qbar qbar -> qbar qbar
static const ColourLines cqbqbqbqb[2]={ColourLines("-1 -2 -5,-3 2 -4"),
ColourLines("-1 -2 -4,-3 2 -5")};
// colour lines for q qbar -> q qbar
static const ColourLines cqqbqqb[2]={ColourLines("1 3 4,-2 -3 -5"),
ColourLines("1 2 -3,4 -2 -5")};
// select the colour flow (as all ready picked just insert answer)
Selector<const ColourLines *> sel;
switch(abs(diag->id())) {
// gg -> gg subprocess
case 1:
if(_flow==1) {
sel.insert(0.5, &cgggg[0]);
sel.insert(0.5, &cgggg[1]);
}
else {
sel.insert(0.5, &cgggg[4]);
sel.insert(0.5, &cgggg[5]);
}
break;
case 2:
if(_flow==2) {
sel.insert(0.5, &cgggg[6]);
sel.insert(0.5, &cgggg[7]);
}
else {
sel.insert(0.5, &cgggg[8]);
sel.insert(0.5, &cgggg[9]);
}
break;
case 3:
if(_flow==1) {
sel.insert(0.5, &cgggg[2]);
sel.insert(0.5, &cgggg[3]);
}
else {
sel.insert(0.5, &cgggg[10]);
sel.insert(0.5, &cgggg[11]);
}
break;
// gg -> q qbar subprocess
case 4: case 5:
sel.insert(1.0, &cggqq[abs(diag->id())-4]);
break;
case 6:
sel.insert(1.0, &cggqq[1+_flow]);
break;
// q qbar -> gg subprocess
case 7: case 8:
sel.insert(1.0, &cqqgg[abs(diag->id())-7]);
break;
case 9:
sel.insert(1.0, &cqqgg[1+_flow]);
break;
// q g -> q g subprocess
case 10: case 11:
sel.insert(1.0, &cqgqg[abs(diag->id())-10]);
break;
case 12:
sel.insert(1.0, &cqgqg[1+_flow]);
break;
// q g -> q g subprocess
case 13: case 14:
sel.insert(1.0, &cqbgqbg[abs(diag->id())-13]);
break;
case 15:
sel.insert(1.0, &cqbgqbg[1+_flow]);
break;
// q q -> q q subprocess
case 16: case 17:
sel.insert(1.0, &cqqqq[abs(diag->id())-16]);
break;
// qbar qbar -> qbar qbar subprocess
case 18: case 19:
sel.insert(1.0, &cqbqbqbqb[abs(diag->id())-18]);
break;
// q qbar -> q qbar subprocess
case 20: case 21:
sel.insert(1.0, &cqqbqqb[abs(diag->id())-20]);
break;
}
return sel;
}
double MEQCD2to2::me2() const {
// total matrix element
double me(0.);
// gg initiated processes
if(mePartonData()[0]->id()==ParticleID::g&&mePartonData()[1]->id()==ParticleID::g) {
// gg -> gg
if(mePartonData()[2]->id()==ParticleID::g) {
VectorWaveFunction g1w(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
VectorWaveFunction g3w(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g1,g2,g3,g4;
for(unsigned int ix=0;ix<2;++ix) {
g1w.reset(2*ix);g1.push_back(g1w);
g2w.reset(2*ix);g2.push_back(g2w);
g3w.reset(2*ix);g3.push_back(g3w);
g4w.reset(2*ix);g4.push_back(g4w);
}
// calculate the matrix element
me = gg2ggME(g1,g2,g3,g4,0);
}
// gg -> q qbar
else {
VectorWaveFunction g1w(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qw(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction qbarw(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g1,g2;
vector<SpinorBarWaveFunction> q;
vector<SpinorWaveFunction> qbar;
for(unsigned int ix=0;ix<2;++ix) {
g1w.reset(2*ix);g1.push_back(g1w);
g2w.reset(2*ix);g2.push_back(g2w);
qw.reset(ix);q.push_back(qw);
qbarw.reset(ix);qbar.push_back(qbarw);
}
// calculate the matrix element
me=gg2qqbarME(g1,g2,q,qbar,0);
}
}
// quark first processes
else if(mePartonData()[0]->id()>0) {
// q g -> q g
if(mePartonData()[1]->id()==ParticleID::g) {
SpinorWaveFunction qinw(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction qoutw(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g2,g4;
vector<SpinorWaveFunction> qin;
vector<SpinorBarWaveFunction> qout;
for(unsigned int ix=0;ix<2;++ix) {
qinw.reset(ix);qin.push_back(qinw);
g2w.reset(2*ix);g2.push_back(g2w);
qoutw.reset(ix);qout.push_back(qoutw);
g4w.reset(2*ix);g4.push_back(g4w);
}
// calculate the matrix element
me = qg2qgME(qin,g2,qout,g4,0);
}
else if(mePartonData()[1]->id()<0) {
// q qbar initiated processes( q qbar -> gg)
if(mePartonData()[2]->id()==ParticleID::g) {
SpinorWaveFunction qw(meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction qbarw(meMomenta()[1],mePartonData()[1],incoming);
VectorWaveFunction g1w(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g2w(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g1,g2;
vector<SpinorWaveFunction> q;
vector<SpinorBarWaveFunction> qbar;
for(unsigned int ix=0;ix<2;++ix) {
qw.reset(ix);q.push_back(qw);
qbarw.reset(ix);qbar.push_back(qbarw);
g1w.reset(2*ix);g1.push_back(g1w);
g2w.reset(2*ix);g2.push_back(g2w);
}
// calculate the matrix element
me = qqbar2ggME(q,qbar,g1,g2,0);
}
// q qbar to q qbar
else {
SpinorWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorWaveFunction> q1,q4;
vector<SpinorBarWaveFunction> q2,q3;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qqbar2qqbarME(q1,q2,q3,q4,0);
}
}
// q q -> q q
else if(mePartonData()[1]->id()>0) {
SpinorWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming);
SpinorWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorBarWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing);
SpinorBarWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorWaveFunction> q1,q2;
vector<SpinorBarWaveFunction> q3,q4;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qq2qqME(q1,q2,q3,q4,0);
}
}
// antiquark first processes
else if(mePartonData()[0]->id()<0) {
// qbar g -> qbar g
if(mePartonData()[1]->id()==ParticleID::g) {
SpinorBarWaveFunction qinw(meMomenta()[0],mePartonData()[0],incoming);
VectorWaveFunction g2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorWaveFunction qoutw(meMomenta()[2],mePartonData()[2],outgoing);
VectorWaveFunction g4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<VectorWaveFunction> g2,g4;
vector<SpinorBarWaveFunction> qin;
vector<SpinorWaveFunction> qout;
for(unsigned int ix=0;ix<2;++ix) {
qinw.reset(ix);qin.push_back(qinw);
g2w.reset(2*ix);g2.push_back(g2w);
qoutw.reset(ix);qout.push_back(qoutw);
g4w.reset(2*ix);g4.push_back(g4w);
}
// calculate the matrix element
me = qbarg2qbargME(qin,g2,qout,g4,0);
}
// qbar qbar -> qbar qbar
else if(mePartonData()[1]->id()<0) {
SpinorBarWaveFunction q1w(meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction q2w(meMomenta()[1],mePartonData()[1],incoming);
SpinorWaveFunction q3w(meMomenta()[2],mePartonData()[2],outgoing);
SpinorWaveFunction q4w(meMomenta()[3],mePartonData()[3],outgoing);
vector<SpinorBarWaveFunction> q1,q2;
vector<SpinorWaveFunction> q3,q4;
for(unsigned int ix=0;ix<2;++ix) {
q1w.reset(ix);q1.push_back(q1w);
q2w.reset(ix);q2.push_back(q2w);
q3w.reset(ix);q3.push_back(q3w);
q4w.reset(ix);q4.push_back(q4w);
}
// calculate the matrix element
me = qbarqbar2qbarqbarME(q1,q2,q3,q4,0);
}
}
else throw Exception() << "Unknown process in MEQCD2to2::me2()"
<< Exception::abortnow;
// return the answer
return me;
}
void MEQCD2to2::constructVertex(tSubProPtr sub) {
// extract the particles in the hard process
ParticleVector hard;
hard.push_back(sub->incoming().first);hard.push_back(sub->incoming().second);
hard.push_back(sub->outgoing()[0]);hard.push_back(sub->outgoing()[1]);
// order of particles
unsigned int order[4]={0,1,2,3};
// identify the process and calculate the matrix element
if(hard[0]->id()==ParticleID::g&&hard[1]->id()==ParticleID::g) {
// gg -> gg
if(hard[2]->id()==ParticleID::g) {
vector<VectorWaveFunction> g1,g2,g3,g4;
VectorWaveFunction(g1,hard[0],incoming,false,true,true);
VectorWaveFunction(g2,hard[1],incoming,false,true,true);
VectorWaveFunction(g3,hard[2],outgoing,true ,true,true);
VectorWaveFunction(g4,hard[3],outgoing,true ,true,true);
g1[1]=g1[2];g2[1]=g2[2];g3[1]=g3[2];g4[1]=g4[2];
gg2ggME(g1,g2,g3,g4,_flow);
}
// gg -> q qbar
else {
if(hard[2]->id()<0) swap(order[2],order[3]);
vector<VectorWaveFunction> g1,g2;
vector<SpinorBarWaveFunction> q;
vector<SpinorWaveFunction> qbar;
VectorWaveFunction( g1,hard[ 0 ],incoming,false,true,true);
VectorWaveFunction( g2,hard[ 1 ],incoming,false,true,true);
SpinorBarWaveFunction(q ,hard[order[2]],outgoing,true ,true);
SpinorWaveFunction( qbar,hard[order[3]],outgoing,true ,true);
g1[1]=g1[2];g2[1]=g2[2];
gg2qqbarME(g1,g2,q,qbar,_flow);
}
}
else if(hard[0]->id()==ParticleID::g||hard[1]->id()==ParticleID::g) {
if(hard[0]->id()==ParticleID::g) swap(order[0],order[1]);
if(hard[2]->id()==ParticleID::g) swap(order[2],order[3]);
// q g -> q g
if(hard[order[0]]->id()>0) {
vector<VectorWaveFunction> g2,g4;
vector<SpinorWaveFunction> qin;
vector<SpinorBarWaveFunction> qout;
SpinorWaveFunction( qin,hard[order[0]],incoming,false,true);
VectorWaveFunction( g2,hard[order[1]],incoming,false,true,true);
SpinorBarWaveFunction(qout,hard[order[2]],outgoing,true ,true);
VectorWaveFunction( g4,hard[order[3]],outgoing,true ,true,true);
g2[1]=g2[2];g4[1]=g4[2];
qg2qgME(qin,g2,qout,g4,_flow);
}
// qbar g -> qbar g
else {
vector<VectorWaveFunction> g2,g4;
vector<SpinorBarWaveFunction> qin;
vector<SpinorWaveFunction> qout;
SpinorBarWaveFunction( qin,hard[order[0]],incoming,false,true);
VectorWaveFunction( g2,hard[order[1]],incoming,false,true,true);
SpinorWaveFunction( qout,hard[order[2]],outgoing,true ,true);
VectorWaveFunction( g4,hard[order[3]],outgoing,true ,true,true);
g2[1]=g2[2];g4[1]=g4[2];
qbarg2qbargME(qin,g2,qout,g4,_flow);
}
}
else if(hard[0]->id()>0||hard[1]->id()>0) {
if(hard[2]->id()==ParticleID::g) {
if(hard[0]->id()<0) swap(order[0],order[1]);
vector<SpinorBarWaveFunction> qbar;
vector<SpinorWaveFunction> q;
vector<VectorWaveFunction> g3,g4;
SpinorWaveFunction( q ,hard[order[0]],incoming,false,true);
SpinorBarWaveFunction(qbar,hard[order[1]],incoming,false,true);
VectorWaveFunction( g3,hard[ 2 ],outgoing,true ,true,true);
VectorWaveFunction( g4,hard[ 3 ],outgoing,true ,true,true);
g3[1]=g3[2];g4[1]=g4[2];
qqbar2ggME(q,qbar,g3,g4,_flow);
}
// q q -> q q
else if(hard[0]->id()>0&&hard[1]->id()>0) {
if(hard[2]->id()!=hard[0]->id()) swap(order[2],order[3]);
vector<SpinorWaveFunction> q1,q2;
vector<SpinorBarWaveFunction> q3,q4;
SpinorWaveFunction( q1,hard[order[0]],incoming,false,true);
SpinorWaveFunction( q2,hard[order[1]],incoming,false,true);
SpinorBarWaveFunction(q3,hard[order[2]],outgoing,true ,true);
SpinorBarWaveFunction(q4,hard[order[3]],outgoing,true ,true);
qq2qqME(q1,q2,q3,q4,_flow);
}
// q qbar -> q qbar
else {
if(hard[0]->id()<0) swap(order[0],order[1]);
if(hard[2]->id()<0) swap(order[2],order[3]);
vector<SpinorWaveFunction> q1,q4;
vector<SpinorBarWaveFunction> q2,q3;
SpinorWaveFunction( q1,hard[order[0]],incoming,false,true);
SpinorBarWaveFunction(q2,hard[order[1]],incoming,false,true);
SpinorBarWaveFunction(q3,hard[order[2]],outgoing,true ,true);
SpinorWaveFunction( q4,hard[order[3]],outgoing,true ,true);
qqbar2qqbarME(q1,q2,q3,q4,_flow);
}
}
else if (hard[0]->id()<0&&hard[1]->id()<0) {
if(hard[2]->id()!=hard[0]->id()) swap(order[2],order[3]);
vector<SpinorBarWaveFunction> q1,q2;
vector<SpinorWaveFunction> q3,q4;
SpinorBarWaveFunction(q1,hard[order[0]],incoming,false,true);
SpinorBarWaveFunction(q2,hard[order[1]],incoming,false,true);
SpinorWaveFunction( q3,hard[order[2]],outgoing,true ,true);
SpinorWaveFunction( q4,hard[order[3]],outgoing,true ,true);
qbarqbar2qbarqbarME(q1,q2,q3,q4,_flow);
}
else throw Exception() << "Unknown process in MEQCD2to2::constructVertex()"
<< Exception::runerror;
// construct the vertex
HardVertexPtr hardvertex=new_ptr(HardVertex());
// set the matrix element for the vertex
hardvertex->ME(_me);
// set the pointers and to and from the vertex
for(unsigned int ix=0;ix<4;++ix)
hard[order[ix]]->spinInfo()->productionVertex(hardvertex);
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 4:52 PM (14 h, 52 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4016339
Default Alt Text
(214 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment