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 #include using namespace Herwig; using namespace ElectroWeakMatching; using namespace GroupInvariants; using namespace EWProcess; boost::numeric::ublas::matrix 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 R0,G2,Dw,Dz; - - switch (process) { - case QQQQ: - case QQQQiden: - case QtQtQQ: - { - unsigned int numGauge = 4, numBrokenGauge = 12; - R0=boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); - G2=boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw=boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz=boost::numeric::ublas::zero_matrix(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 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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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 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(numBrokenGauge,numGauge); + G2=boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw=boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz=boost::numeric::ublas::zero_matrix(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 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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - Dz = boost::numeric::ublas::zero_matrix(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 output(R0); - boost::numeric::ublas::matrix 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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numGauge); + G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + Dw = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + Dz = boost::numeric::ublas::zero_matrix(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 output(R0); + boost::numeric::ublas::matrix 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 namespace Herwig { using namespace ThePEG; namespace ElectroWeakMatching { /** * The high energy matching */ boost::numeric::ublas::matrix - electroWeakMatching(Energy mu, - Energy2 s, Energy2 t, Energy2 u, - Herwig::EWProcess::Process process, - bool oneLoop); - - /* /\** */ - /* * Spin\f$\frac12\f$ */ - /* *\/ */ - /* boost::numeric::ublas::matrix > */ - /* SpinHalfMatching(Energy highScale, */ - /* Energy2 s, Energy2 t, Energy2 u, */ - /* EWProcess::Process process, */ - /* bool oneLoop, bool includeAlphaS2); */ - - /* /\** */ - /* * Spin\f$1\f$ */ - /* *\/ */ - /* boost::numeric::ublas::matrix > */ - /* Spin1HighMatching(Energy highScale, */ - /* Energy2 s, Energy2 t, Energy2 u, */ - /* EWProcess::Process process, */ - /* bool oneLoop, bool includeAlphaS2); */ - /* /\** */ - /* * Spin\f$0\f$ */ - /* *\/ */ - /* boost::numeric::ublas::matrix > */ - /* 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 describeHerwigElectroWeakReweighter("Herwig::ElectroWeakReweighter", "HwMEEW.so"); void ElectroWeakReweighter::Init() { static ClassDocumentation documentation ("There is no documentation for the ElectroWeakReweighter class"); static Reference interfaceEWCouplings ("EWCouplings", "The object to calculate the electroweak couplings", &ElectroWeakReweighter::EWCouplings_, false, false, true, false, false); static Reference interfaceCollinearSudakov ("CollinearSudakov", "The collinear Sudakov", &ElectroWeakReweighter::collinearSudakov_, false, false, true, false, false); static Reference interfaceSoftSudakov ("SoftSudakov", "The soft Sudakov", &ElectroWeakReweighter::softSudakov_, false, false, true, false, false); } namespace { void axpy_prod_local(const boost::numeric::ublas::matrix & A, const boost::numeric::ublas::matrix > & B, boost::numeric::ublas::matrix > & C) { assert(A.size2()==B.size1()); C.resize(A.size1(),B.size2()); for(unsigned int ix=0;ix > & A, const boost::numeric::ublas::vector > & B, boost::numeric::ublas::vector & C) { assert(A.size2()==B.size()); C.resize(A.size1()); for(unsigned int ix=0;ix > & A, const boost::numeric::ublas::matrix & B, boost::numeric::ublas::matrix > & C) { assert(A.size2()==B.size1()); C.resize(A.size1(),B.size2()); for(unsigned int ix=0;ixinitialize(); 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 > highMatch_val = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,true,true); boost::numeric::ublas::matrix highRunning_val - = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process); + = softSudakov_->highEnergyRunning(highScale,ewScale,s,t,u,process,0); boost::numeric::ublas::matrix ewMatch_val = - ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true); + ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true,0); boost::numeric::ublas::matrix lowRunning_val = - softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process); + softSudakov_->lowEnergyRunning(ewScale,lowScale,s,t,u,process,0); boost::numeric::ublas::matrix collinearHighRunning_val = collinearSudakov_->highEnergyRunning(highScale,ewScale,s,process,false); boost::numeric::ublas::matrix collinearEWMatch_val = collinearSudakov_->electroWeakMatching(ewScale,s,process,true); boost::numeric::ublas::matrix collinearLowRunning_val = collinearSudakov_->lowEnergyRunning(ewScale,lowScale,s,process); boost::numeric::ublas::matrix lowMatchTemp_val = boost::numeric::ublas::zero_matrix(ewMatch_val.size1(),ewMatch_val.size2()); for (unsigned int ii=0; ii temp(highRunning_val.size1(),collinearHighRunning_val.size2()); boost::numeric::ublas::axpy_prod(highRunning_val,collinearHighRunning_val,temp); boost::numeric::ublas::matrix temp2(collinearLowRunning_val.size1(),lowRunning_val.size2()); boost::numeric::ublas::axpy_prod(collinearLowRunning_val,lowRunning_val,temp2); boost::numeric::ublas::matrix 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 > result(temp2.size1(),highMatch_val.size2()); axpy_prod_local(temp2,highMatch_val,result); for(unsigned int ix=0;ix > & eps3, vector > & 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 eps3Para = (m2+t)/den1*p1 -(m2+u)/den1*p2 +(u-t)/den1*p3; LorentzVector eps3Perp = 2./den2*epsilon(p1,p2,p3); LorentzVector eps4Para = (m2+t)/den1*p2 -(m2+u)/den1*p1 +(u-t)/den1*p4; LorentzVector 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 eps3Para( 1., 0.,0.,0.); LorentzVector eps3Perp( 0.,-1.,0.,0.); LorentzVector eps4Para(-1.,0.,0., 0.); LorentzVector 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 eps3Para( 1., 0.,0.,0.); + LorentzVector eps3Perp( 0.,-1.,0.,0.); + LorentzVector eps4Para(-1.,0.,0., 0.); + LorentzVector 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 eps3Para( 1., 0.,0.,0.); + LorentzVector eps3Perp( 0.,-1.,0.,0.); + LorentzVector eps4Para(-1.,0.,0., 0.); + LorentzVector 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 > - 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 > - 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 > eps3,eps4; SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps3,eps4,0); boost::numeric::ublas::matrix bornME = boost::numeric::ublas::zero_matrix(3,3), EWME = boost::numeric::ublas::zero_matrix(3,3); - // testME = boost::numeric::ublas::zero_matrix(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 > current = iq==0 ? qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) : qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave()); for(unsigned int i1=0;i1<2;++i1) { complex 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 > M(5); Complex d34 = eps3[i1].dot(eps4[i2]); complex 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 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 > 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 > highMatch_val - = HighEnergyMatching::highEnergyMatching(highScale,s,t,u,process,!born,false); + matrix > 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 - ewMatch_val = ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,!born); + ewMatch_val = ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,!born,iswap); matrix collinearEWMatch_val = collinearSudakov_->electroWeakMatching(ewScale,s,process,!born); + // EVOLUTION matrix highRunning_val,lowRunning_val, collinearHighRunning_val,collinearLowRunning_val; + // born process if(born) { highRunning_val = identity_matrix(softSudakov_->numberGauge(process)); lowRunning_val = identity_matrix(softSudakov_->numberBrokenGauge(process)); collinearHighRunning_val = identity_matrix(softSudakov_->numberGauge(process)); collinearLowRunning_val = identity_matrix(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 lowMatchTemp_val = zero_matrix(ewMatch_val.size1(),ewMatch_val.size2()); for (unsigned int ii=0; ii temp(highRunning_val.size1(),collinearHighRunning_val.size2()); axpy_prod(highRunning_val,collinearHighRunning_val,temp); matrix temp2(collinearLowRunning_val.size1(),lowRunning_val.size2()); axpy_prod(collinearLowRunning_val,lowRunning_val,temp2); matrix temp3(temp2.size1(),lowMatchTemp_val.size2()); axpy_prod(temp2,lowMatchTemp_val,temp3); temp2.resize(temp3.size1(),temp.size2()); axpy_prod(temp3,temp,temp2); matrix > result(temp2.size1(),highMatch_val.size2()); axpy_prod_local(temp2,highMatch_val,result); - // for(unsigned int ix=0;ixSU3(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 > - 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 > - 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 > eps1,eps2; SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps1,eps2,1); boost::numeric::ublas::matrix bornME = boost::numeric::ublas::zero_matrix(3,3), EWME = boost::numeric::ublas::zero_matrix(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 > current = iq==0 ? qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) : qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave()); for(unsigned int i1=0;i1<2;++i1) { complex d31 = eps1[i1].dot(p3); for(unsigned int i2=0;i2<2;++i2) { boost::numeric::ublas::vector > M(5); Complex d34 = eps1[i1].dot(eps2[i2]); complex 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 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 > + 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 > eps2,eps3; + SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps2,eps3,2); + // compute the matrix elements + boost::numeric::ublas::matrix + bornME = boost::numeric::ublas::zero_matrix(3,3), + EWME = boost::numeric::ublas::zero_matrix(3,3), + testME = boost::numeric::ublas::zero_matrix(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 > current = iq==0 ? + qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) : + qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave()); + for(unsigned int i1=0;i1<2;++i1) { + complex d31 = eps3[i1].dot(p1); + for(unsigned int i2=0;i2<2;++i2) { + boost::numeric::ublas::vector > M(5); + Complex d34 = eps3[i1].dot(eps2[i2]); + complex 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 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 > + 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 > eps1,eps4; + SackGluonPolarizations(p1,p2,p3,p4,s,t,u,ZERO,eps1,eps4,3); + boost::numeric::ublas::matrix + bornME = boost::numeric::ublas::zero_matrix(3,3), + EWME = boost::numeric::ublas::zero_matrix(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 > current = iq==0 ? + qw.dimensionedWave(). leftCurrent(qbarw.dimensionedWave()) : + qw.dimensionedWave().rightCurrent(qbarw.dimensionedWave()); + for(unsigned int i1=0;i1<2;++i1) { + complex d31 = eps1[i1].dot(p3); + for(unsigned int i2=0;i2<2;++i2) { + boost::numeric::ublas::vector > M(5); + Complex d34 = eps1[i1].dot(eps4[i2]); + complex 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 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 > 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 #include 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 mu0.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 Gamma2(Complex U, Complex T) { boost::numeric::ublas::matrix 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 Gamma2w(Complex U, Complex T) { boost::numeric::ublas::matrix output = boost::numeric::ublas::zero_matrix(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 Gamma2Singlet() { - boost::numeric::ublas::matrix output = boost::numeric::ublas::zero_matrix(2,2); + using namespace boost::numeric::ublas; + matrix output = zero_matrix(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 Gamma2SingletST(Complex T) { + using namespace boost::numeric::ublas; + matrix output = zero_matrix(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 Gamma3(Complex U, Complex T) { boost::numeric::ublas::matrix output = boost::numeric::ublas::zero_matrix(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 Gamma3g(Complex U, Complex T) { boost::numeric::ublas::matrix output = boost::numeric::ublas::zero_matrix(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 Gamma3gST(Complex U, Complex T) { + boost::numeric::ublas::matrix output = boost::numeric::ublas::zero_matrix(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 Gamma3Singlet() { boost::numeric::ublas::matrix output = boost::numeric::ublas::zero_matrix(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 describeHerwigSoftSudakov("Herwig::SoftSudakov", "HwMEEW.so"); void SoftSudakov::Init() { static ClassDocumentation 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 SoftSudakov::evaluateSoft(boost::numeric::ublas::matrix & G3, boost::numeric::ublas::matrix & G2, boost::numeric::ublas::matrix & 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 gamma(NN,NN); for(row_=0;row_ 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 G1, G2, G3; + matrix 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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - boost::numeric::ublas::matrix gam3 = Gamma3(U,T); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); + matrix gam3 = Gamma3(U,T); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - boost::numeric::ublas::matrix gam3 = Gamma3(U,T); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); + matrix gam3 = Gamma3(U,T); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - boost::numeric::ublas::matrix gam3 = Gamma3(U,T); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); + matrix gam3 = Gamma3(U,T); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - boost::numeric::ublas::matrix gam3g = Gamma3g(U,T); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); + matrix 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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam1 = Gamma1(2./3.,0.,T,U); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); Complex gam1 = Gamma1(-1./3.,0.,T,U); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(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(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam1 = Gamma1(-1.,0.,T,U); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); - G2 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); - G3 = boost::numeric::ublas::zero_matrix(numBrokenGauge,numBrokenGauge); + G1 = zero_matrix(numBrokenGauge,numBrokenGauge); + G2 = zero_matrix(numBrokenGauge,numBrokenGauge); + G3 = zero_matrix(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(G1.size1()); + return identity_matrix(G1.size1()); } else { return evaluateSoft(G3,G2,G1,EWScale,lowScale,false); } } boost::numeric::ublas::matrix 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 G1,G2,G3; + matrix 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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - boost::numeric::ublas::matrix gam3 = Gamma3(U,T); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(numGauge,numGauge); + matrix 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 gam2 = Gamma2(U,T); + matrix 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(numGauge,numGauge); + G1 = zero_matrix(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(numGauge,numGauge); + G1 = zero_matrix(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(numGauge,numGauge); + G1 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); + G1 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = Gamma3g(U,T); - Complex gam2s = Gamma2Singlet()(0,0); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(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(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(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(numGauge,numGauge); + G2 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(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(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(numGauge,numGauge); G1(0,0) = Gamma1(-1.0); break; case EEPhiPhi: + assert(iswap==0); numGauge = 1; - G1 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G2 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); - G3 = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + G1 = zero_matrix(numGauge,numGauge); + G2 = zero_matrix(numGauge,numGauge); + G3 = zero_matrix(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 #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 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 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 evaluateSoft(boost::numeric::ublas::matrix & G3, boost::numeric::ublas::matrix & G2, boost::numeric::ublas::matrix & 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 G1_; /** * */ boost::numeric::ublas::matrix G2_; /** * */ boost::numeric::ublas::matrix 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(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(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::initMEQCD2to2; // Definition of the static class description member. void MEQCD2to2::Init() { static ClassDocumentation documentation ("The MEQCD2to2 class implements the QCD 2->2 processes in hadron-hadron" " collisions"); static Parameter interfaceMaximumFlavour ("MaximumFlavour", "The maximum flavour of the quarks in the process", &MEQCD2to2::_maxflavour, 5, 1, 5, false, false, Interface::limited); static Switch 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 MEQCD2to2::diagrams(const DiagramVector & diags) const { // select the diagram, this is easy for us as we have already done it Selector 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 &g1, vector &g2, vector & q, vector & 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 & q, vector & qbar, vector &g1, vector &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 & qin, vector &g2, vector & qout, vector &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 &g1,vector &g2, vector &g3,vector &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 & qin, vector &g2, vector & qout, vector &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 & q1, vector & q2, vector & q3, vector & 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 & q1, vector & q2, vector & q3, vector & 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 & q1, vector & q2, vector & q3, vector & 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 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 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 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 g1,g2; vector q; vector 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 g2,g4; vector qin; vector 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 g1,g2; vector q; vector 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 q1,q4; vector 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 q1,q2; vector 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 g2,g4; vector qin; vector 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 q1,q2; vector 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 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 g1,g2; vector q; vector 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 g2,g4; vector qin; vector 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 g2,g4; vector qin; vector 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 qbar; vector q; vector 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 q1,q2; vector 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 q1,q4; vector 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 q1,q2; vector 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); }