Page MenuHomeHEPForge

No OneTemporary

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

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)

Event Timeline