diff --git a/MatrixElement/EW/ElectroWeakMatching.cc b/MatrixElement/EW/ElectroWeakMatching.cc --- a/MatrixElement/EW/ElectroWeakMatching.cc +++ b/MatrixElement/EW/ElectroWeakMatching.cc @@ -1,915 +1,1015 @@ // -*- 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,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); 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; + Complex w0(0.),w1(0.),w2(0.); + Complex z1(0.),z2(0.),z3(0.),z4(0.),z5(0.); + boost::numeric::ublas::matrix gam2; + if(iswap==0) { + w0 = 0.5*I*pi; + w1 = -0.5*(T-U); + w2 = -0.5*(T+U); + z1 = 2.0*g11*g21*(T-U) - I*pi*(g11*g11+g21*g21); + z2 = 2.0*g21*g12*(T-U) - I*pi*(g21*g21+g12*g12); + z3 = 2.0*g22*g11*(T-U) - I*pi*(g22*g22+g11*g11); + z4 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22); + z5 = (g11*g21+g12*g22)*T - (g21*g12+g11*g22)*U + - 0.5*I*pi*(g11*g11+g12*g12+g21*g21+g22*g22); + gam2 = Gamma2(U,T); + } + else if(iswap==1) { + w0 = -0.5*(T-I*pi); + w1 = 0.5*U; + w2 = -0.5*(U-2.*T); + z1 = 2.0*g11*g21*(-U) + ( T- I*pi)*(g11*g11+g21*g21); + z2 = 2.0*g21*g12*(-U) + ( T- I*pi)*(g21*g21+g12*g12); + z3 = 2.0*g22*g11*(-U) + ( T- I*pi)*(g22*g22+g11*g11); + z4 = 2.0*g12*g22*(-U) + ( T- I*pi)*(g12*g12+g22*g22); + z5 = -(g11*g21+g12*g22)*T - (g21*g12+g11*g22)*(U-T) + + 0.5*(T-I*pi)*(g11*g11+g12*g12+g21*g21+g22*g22); + gam2 = Gamma2ST(U,T); + } + else if(iswap==2) { + w0 = -0.5*(U-I*pi); + w1 = -0.5*T; + w2 = -0.5*(T-2.*U); + z1 = 2.0*g11*g21*T +(U-I*pi)*(g11*g11+g21*g21); + z2 = 2.0*g21*g12*T +(U-I*pi)*(g21*g21+g12*g12); + z3 = 2.0*g22*g11*T +(U-I*pi)*(g22*g22+g11*g11); + z4 = 2.0*g12*g22*T +(U-I*pi)*(g12*g12+g22*g22); + z5 = (g11*g21+g12*g22)*(T-U) + (g21*g12+g11*g22)*U + + 0.5*(U-I*pi)*(g11*g11+g12*g12+g21*g21+g22*g22); + gam2 = Gamma2SU(U,T); + } + else + assert(false); + // Dw for(unsigned int ix=0;ix gam2 = Gamma2(U,T); + // G2 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: { - 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; + Complex w1(0.),z1(0.),z2(0.); + if(iswap==0) { + w1 = 0.25*I*pi; + z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22); + z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22); + G2 = Gamma2Singlet(); + } + else if(iswap==1) { + w1 = -0.25*(T-I*pi); + z1 = -2.0*g11*g22*U +(T-I*pi)*(g11*g11+g22*g22); + z2 = -2.0*g12*g22*U +(T-I*pi)*(g12*g12+g22*g22); + G2 = Gamma2SingletST(T); + } + else if(iswap==2) { + w1 = -0.25*(U-I*pi); + z1 = 2.0*g11*g22*T +(U-I*pi)*(g11*g11+g22*g22); + z2 = 2.0*g12*g22*T +(U-I*pi)*(g12*g12+g22*g22); + G2 = Gamma2SingletSU(U); + } + else + assert(false); 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; + Complex w1(0.),z1(0.),z2(0.); + if(iswap==0) { + w1 = 0.25*I*pi; + z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22); + z2 = 2.0*g12*g22*(T-U) - I*pi*(g12*g12+g22*g22); + G2 = Gamma2Singlet(); + } + else if(iswap==1) { + w1 =-0.25*(T-I*pi); + z1 =-2.0*g11*g22*U + (T-I*pi)*(g11*g11+g22*g22); + z2 =-2.0*g12*g22*U + (T-I*pi)*(g12*g12+g22*g22); + G2 = Gamma2SingletST(T); + } + else if(iswap==2) { + w1 =-0.25*(U-I*pi); + z1 = 2.0*g11*g22*T + (U-I*pi)*(g11*g11+g22*g22); + z2 = 2.0*g12*g22*T + (U-I*pi)*(g12*g12+g22*g22); + G2 = Gamma2Singlet(); + } + else + assert(false); 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); + Complex z1(0.); + if(iswap==0) { + z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22); + } + else if(iswap==1) { + z1 =-2.0*g11*g22*U + (T-I*pi)*(g11*g11+g22*g22); + } + else if(iswap==2) { + 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); + Complex z1(0.); + if(iswap==0) { + z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22); + } + else if(iswap==1) { + z1 =-2.0*g11*g22*U + (T-I*pi)*(g11*g11+g22*g22); + } + else if(iswap==2) { + z1 = 2.0*g11*g22*T + (U-I*pi)*(g11*g11+g22*g22); + } + else + assert(false); 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); + Complex z1(0.); + if(iswap==0) { + z1 = 2.0*g11*g22*(T-U) - I*pi*(g11*g11+g22*g22); + } + else if(iswap==1) { + z1 =-2.0*g11*g22*U +(T-I*pi)*(g11*g11+g22*g22); + } + else if(iswap==2) { + z1 = 2.0*g11*g22*T +(U-I*pi)*(g11*g11+g22*g22); + } + else + assert(false); 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 if(iswap==2) { + w2 = 0.25*(-U+I*pi); + z3 = (U-I*pi)*g1*g1; + z4 = (U-I*pi)*g2*g2; + G2(0,0) = G2(1,1) = G2(2,2) = Gamma2SingletSU(U)(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 if(iswap==2) { + z1 = (U-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/ElectroWeakReweighter.cc b/MatrixElement/EW/ElectroWeakReweighter.cc --- a/MatrixElement/EW/ElectroWeakReweighter.cc +++ b/MatrixElement/EW/ElectroWeakReweighter.cc @@ -1,907 +1,1924 @@ // -*- 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"; + // 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); + else if((subProcess()->incoming().first ->id() > 0 && + subProcess()->incoming().first ->id()<= 5 && + subProcess()->incoming().second->id() < 0 && + subProcess()->incoming().second->id()>=-5) || + (subProcess()->incoming().second->id() > 0 && + subProcess()->incoming().second->id()<= 5 && + subProcess()->incoming().first ->id() < 0 && + subProcess()->incoming().first ->id()>=-5)) { + // identical flavour q qbar + if(subProcess()->incoming().first ->id() == -subProcess()->incoming().second->id()) { + // q qbar -> gg + if(subProcess()->outgoing()[0]->id()==ParticleID::g && + subProcess()->outgoing()[1]->id()==ParticleID::g) + return reweightqqbargg(); + // q qbar -> q' q'bar + else if(subProcess()->outgoing()[0]->id() == -subProcess()->outgoing()[1]->id() && + abs(subProcess()->outgoing()[0]->id())<=5) + return reweightqqbarqqbarS(); + } + // different flavour q qbar + else { + if((subProcess()->outgoing()[0]->id() > 0 && + subProcess()->outgoing()[0]->id()<= 5 && + subProcess()->outgoing()[1]->id() < 0 && + subProcess()->outgoing()[1]->id()>=-5) || + (subProcess()->outgoing()[1]->id() > 0 && + subProcess()->outgoing()[1]->id()<= 5 && + subProcess()->outgoing()[0]->id() < 0 && + subProcess()->outgoing()[0]->id()>=-5)) { + return reweightqqbarqqbarT(); + } + 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); } + // processes with q q initial-state + else if( subProcess()->incoming().first ->id()> 0 && + subProcess()->incoming().first ->id()<=5 && + subProcess()->incoming().second->id()> 0 && + subProcess()->incoming().second->id()<=5 ) { + if(subProcess()->outgoing()[0]->id()> 0 && + subProcess()->outgoing()[0]->id()<=5 && + subProcess()->outgoing()[1]->id()> 0 && + subProcess()->outgoing()[1]->id()<=5) + return reweightqqqq(); + else + assert(false); + } + // processes with qbar qbar initial-state + else if( subProcess()->incoming().first ->id()< 0 && + subProcess()->incoming().first ->id()>= -5 && + subProcess()->incoming().second->id()< 0 && + subProcess()->incoming().second->id()>= -5 ) { + if(subProcess()->outgoing()[0]->id()< 0 && + subProcess()->outgoing()[0]->id()>= -5 && + subProcess()->outgoing()[1]->id()< 0 && + subProcess()->outgoing()[1]->id()>= -5) + return reweightqbarqbarqbarqbar(); + 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,0); boost::numeric::ublas::matrix ewMatch_val = ElectroWeakMatching::electroWeakMatching(ewScale,s,t,u,process,true,0); boost::numeric::ublas::matrix lowRunning_val = 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; // 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 ,0); EWQQGGweights = evaluateRunning(EWProcess::QQGG,s,t,u,false,0); } else { 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 ,0); EWRRGGweights = evaluateRunning(EWProcess::UUGG,s,t,u,false,0); } else { bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,0); EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false,0); } SpinorWaveFunction qw(p1,q ,incoming); SpinorBarWaveFunction qbarw(p2,qbar,incoming); vector > 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); 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]); // 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)); } } } } } 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; } boost::numeric::ublas::matrix > ElectroWeakReweighter::evaluateRunning(EWProcess::Process process, Energy2 s, 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; 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 if(iswap==2) + highMatch_val = HighEnergyMatching::highEnergyMatching(highScale,u,t,s,process,!born,false); else assert(false); // low energy matching matrix 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,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); // reset the couplings coupling()->SU3(SU3save); coupling()-> EW( EWsave); // return the answer return result; } double ElectroWeakReweighter::reweightggqqbar() const { // momenta and invariants Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum(); Lorentz5Momentum p2 = subProcess()->incoming().second->momentum(); Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum(); Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum(); tcPDPtr qbar = subProcess()->outgoing()[0]->dataPtr(); tcPDPtr q = subProcess()->outgoing()[1]->dataPtr(); if(q->id()<0) { swap(p3,p4 ); swap(q ,qbar); } Energy2 s = (p1+p2).m2(); Energy2 t = (p1-p4).m2(); Energy2 u = (p1-p3).m2(); // boost to partonic rest frame and rescale momenta of outgoing // so zero mass Lorentz5Momentum psum=p1+p2; LorentzRotation boost(-psum.boostVector()); p1 *= boost; p2 *= boost; p3 *= boost; p4 *= boost; p3.setMass(ZERO); p3.rescaleRho(); p4.setMass(ZERO); p4.rescaleRho(); // LO and EW matrix element coefficents boost::numeric::ublas::matrix > bornQQGGweights,bornRRGGweights,EWQQGGweights,EWRRGGweights; // quark left doublet if(q->id()<5) { 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 ,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 ,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 { bornRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,true ,0); EWRRGGweights = evaluateRunning(EWProcess::DDGG,s,t,u,false,0); } 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; } + +double ElectroWeakReweighter::reweightqqbarqqbarS() const { + // momenta and invariants + Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum(); + tcPDPtr q1 = subProcess()->incoming().first ->dataPtr(); + Lorentz5Momentum p2 = subProcess()->incoming().second->momentum(); + tcPDPtr q1bar = subProcess()->incoming().second->dataPtr(); + if(q1->id()<0) { + swap(p1,p2 ); + swap(q1 ,q1bar); + } + Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum(); + tcPDPtr q2bar = subProcess()->outgoing()[0]->dataPtr(); + Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum(); + tcPDPtr q2 = subProcess()->outgoing()[1]->dataPtr(); + if(q2bar->id()>0) { + swap(p3,p4 ); + swap(q2 ,q2bar); + } + 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; + p3.setMass(ZERO); + p3.rescaleRho(); + p4.setMass(ZERO); + p4.rescaleRho(); + // LO and EW corrected matrix element coefficients + boost::numeric::ublas::matrix > + bornLLLLWeights,bornLLRRWeights,bornRRLLWeights,bornRRRRWeights, + EWLLLLWeights,EWLLRRWeights,EWRRLLWeights,EWRRRRWeights; + bool ident = q1->id()==q2->id(); + // LL -> LL + if((q1->id()<=4&& q2->id()<=4)|| (q1->id()==5 && q2->id()==5)) { + if(!ident) { + bornLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,true ,0); + EWLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,false,0); + } + else { + bornLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,true ,0); + EWLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,false,0); + } + } + else if(q1->id()==5 || q2->id()>=5) { + bornLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,true ,0); + EWLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,false,0); + } + else + assert(false); + // RR -> LL + if(q1->id()%2==0) { + if(q2->id()<5) { + bornRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,0); + EWRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,0); + } + else { + bornRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,0); + EWRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,0); + } + } + else { + if(q2->id()<5) { + bornRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,0); + EWRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,0); + } + else { + bornRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,0); + EWRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,0); + } + } + // LL -> RR + if(q1->id()<=4) { + if(q2->id()%2!=0) { + bornLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,0); + EWLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,0); + } + else if (q2->id()==6) { + bornLLRRWeights = evaluateRunning(EWProcess::QQtRtR,s,t,u,true ,0); + EWLLRRWeights = evaluateRunning(EWProcess::QQtRtR,s,t,u,false,0); + } + else { + bornLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,0); + EWLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,0); + } + } + else { + if(q2->id()%2!=0) { + bornLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,0); + EWLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,0); + } + else { + bornLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,0); + EWLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,0); + } + } + // RR -> RR + if(q1->id()%2==0) { + if(q2->id()==6) { + bornRRRRWeights = evaluateRunning(EWProcess::tRtRUU,s,t,u,true ,0); + EWRRRRWeights = evaluateRunning(EWProcess::tRtRUU,s,t,u,false,0); + } + else if(q2->id()%2==0) { + if(ident) { + bornRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,true ,0); + EWRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,false,0); + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,true ,0); + EWRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,false,0); + } + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,0); + EWRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,false,0); + } + } + else { + if(q2->id()==6) { + bornRRRRWeights = evaluateRunning(EWProcess::tRtRDD,s,t,u,true ,0); + EWRRRRWeights = evaluateRunning(EWProcess::tRtRDD,s,t,u,false,0); + } + else if(q2->id()%2==0) { + bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,0); + EWRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,false,0); + } + else { + if(ident) { + bornRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,true ,0); + EWRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,false,0); + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,true ,0); + EWRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,false,0); + } + } + } + // extra terms for identical particles + boost::numeric::ublas::matrix > + borntChannelWeights,EWtChannelWeights; + if(ident) { + if(q1->id()%2==0) { + borntChannelWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,1); + EWtChannelWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,1); + } + else if(q1->id()==5) { + borntChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,1); + EWtChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,1); + } + else { + borntChannelWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,1); + EWtChannelWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,1); + } + } + SpinorWaveFunction q1w(p1,q1 ,incoming); + SpinorBarWaveFunction q1barw(p2,q1bar,incoming); + SpinorWaveFunction q2barw(p3,q2bar,outgoing); + SpinorBarWaveFunction q2w(p4,q2 ,outgoing); + boost::numeric::ublas::matrix + bornME = boost::numeric::ublas::zero_matrix(2,2), + EWME = boost::numeric::ublas::zero_matrix(2,2); + for(unsigned int iq1=0;iq1<2;++iq1) { + if(iq1==0) { + q1w.reset (0); + q1barw.reset(1); + } + else { + q1w.reset (1); + q1barw.reset(0); + } + LorentzVector > current1 = + q1w.dimensionedWave().vectorCurrent(q1barw.dimensionedWave()); + for(unsigned int iq2=0;iq2<2;++iq2) { + if(iq2==0) { + q2w.reset (0); + q2barw.reset(1); + } + else { + q2w.reset (1); + q2barw.reset(0); + } + LorentzVector > current2 = + q2barw.dimensionedWave().vectorCurrent(q2w.dimensionedWave()); + complex amp = current1.dot(current2); + vector Cborn(2),CEW(2); + // amplitudes + if(iq1==0) { + // LL + if(iq2==0) { + unsigned int ioff; + if(q1->id()%2==0) { + ioff = q2->id()%2==0 ? 0 : 2; + } + else { + ioff = q2->id()%2==0 ? 1 : 3; + } + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornLLLLWeights(6*ix+ioff,0); + CEW [ix] = amp* EWLLLLWeights(6*ix+ioff,0); + } + } + // LR + else { + unsigned int ioff = q1->id()%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornLLRRWeights(2*ix+ioff,0); + CEW [ix] = amp* EWLLRRWeights(2*ix+ioff,0); + } + } + } + else { + if(iq2==0) { + unsigned int ioff=q2->id()%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornRRLLWeights(2*ix+ioff,0); + CEW [ix] = amp* EWRRLLWeights(2*ix+ioff,0); + } + } + else { + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornRRRRWeights(ix,0); + CEW [ix] = amp* EWRRRRWeights(ix,0); + } + } + } + // square + for(unsigned int ix=0;ix<2;++ix) { + for(unsigned int iy=0;iy<2;++iy) { + bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]); + EWME (ix,iy) += CEW [ix]*conj(CEW [iy]); + } + } + } + } + // extra t-channel pieces if identical flavours + if(ident) { + for(unsigned int iq1=0;iq1<2;++iq1) { + q1w.reset(iq1); + q2w.reset(iq1); + LorentzVector > current1 = + q1w.dimensionedWave().vectorCurrent(q2w.dimensionedWave()); + q1barw.reset(iq1); + q2barw.reset(iq1); + LorentzVector > current2 = + q2barw.dimensionedWave().vectorCurrent(q1barw.dimensionedWave()); + complex amp = current1.dot(current2); + vector Cborn(2),CEW(2); + unsigned int ioff = q1->id()%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*borntChannelWeights(2*ix+ioff,0); + CEW [ix] = amp* EWtChannelWeights(2*ix+ioff,0); + } + // square + for(unsigned int ix=0;ix<2;++ix) { + for(unsigned int iy=0;iy<2;++iy) { + bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]); + EWME (ix,iy) += CEW [ix]*conj(CEW [iy]); + } + } + } + } + // colour factors + double born = 2.*real(bornME(0,0))+9.*real(bornME(1,1)); + double EW = 2.*real( EWME(0,0))+9.*real( EWME(1,1)); + return EW/born; +} + +double ElectroWeakReweighter::reweightqqbarqqbarT() const { + // momenta and invariants + Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum(); + tcPDPtr q1 = subProcess()->incoming().first ->dataPtr(); + Lorentz5Momentum p2 = subProcess()->incoming().second->momentum(); + tcPDPtr q1bar = subProcess()->incoming().second->dataPtr(); + if(q1->id()<0) { + swap(p1,p2 ); + swap(q1 ,q1bar); + } + Lorentz5Momentum p3 = subProcess()->outgoing()[0]->momentum(); + tcPDPtr q2bar = subProcess()->outgoing()[0]->dataPtr(); + Lorentz5Momentum p4 = subProcess()->outgoing()[1]->momentum(); + tcPDPtr q2 = subProcess()->outgoing()[1]->dataPtr(); + if(q2bar->id()>0) { + swap(p3,p4 ); + swap(q2 ,q2bar); + } + 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; + p3.setMass(ZERO); + p3.rescaleRho(); + p4.setMass(ZERO); + p4.rescaleRho(); + assert(q1==q2 && q1bar==q2bar); + assert( q1->id() != -q1bar->id() && q2->id() != -q2bar->id() ); + // LO and EW corrected matrix element coefficients + boost::numeric::ublas::matrix > + bornLLLLWeights,bornLLRRWeights,bornRRLLWeights,bornRRRRWeights, + EWLLLLWeights,EWLLRRWeights,EWRRLLWeights,EWRRRRWeights; + // LL + if( q1->id() == ParticleID::b || + q1bar->id() == ParticleID::bbar ) { + bornLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,true ,1); + EWLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,false,1); + } + else { + bornLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,true ,1); + EWLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,false,1); + } + // RR -> LL + if(q1->id()%2==0) { + if(q1bar->id()==ParticleID::bbar) { + bornRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,1); + EWRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,1); + } + else { + bornRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,1); + EWRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,1); + } + } + else { + if(q1bar->id()==ParticleID::bbar) { + bornRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,1); + EWRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,1); + } + else { + bornRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,1); + EWRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,1); + } + } + // LL -> RR + if(abs(q1bar->id())%2==0) { + if(q1->id()==ParticleID::b) { + bornLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,1); + EWLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,1); + } + else { + bornLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,1); + EWLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,1); + } + } + else { + if(q1->id()==ParticleID::b) { + bornLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,1); + EWLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,1); + } + else { + bornLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,1); + EWLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,1); + } + } + // RR -> RR + if(q1->id()%2==0) { + if(abs(q1bar->id())%2==0) { + bornRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,true ,1); + EWRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,false,1); + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,1); + EWRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,false,1); + } + } + else { + if(abs(q1bar->id())%2==0) { + bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,1); + EWRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,false,1); + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,true ,1); + EWRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,false,1); + } + } + // calculate the spinors + SpinorWaveFunction q1w(p1,q1 ,incoming); + SpinorBarWaveFunction q1barw(p2,q1bar,incoming); + SpinorWaveFunction q2barw(p3,q2bar,outgoing); + SpinorBarWaveFunction q2w(p4,q2 ,outgoing); + boost::numeric::ublas::matrix + bornME = boost::numeric::ublas::zero_matrix(2,2), + EWME = boost::numeric::ublas::zero_matrix(2,2); + for(unsigned int iq1=0;iq1<2;++iq1) { + q1w.reset(iq1); + q2w.reset(iq1); + LorentzVector > current1 = + q1w.dimensionedWave().vectorCurrent(q2w.dimensionedWave()); + for(unsigned int iq2=0;iq2<2;++iq2) { + q1barw.reset(iq2); + q2barw.reset(iq2); + LorentzVector > current2 = + q2barw.dimensionedWave().vectorCurrent(q1barw.dimensionedWave()); + // calculate the amplitude + complex amp = current1.dot(current2); + vector Cborn(2),CEW(2); + if(iq1==0) { + // LL RR + if(iq2==0) { + unsigned int ioff = q1->id()%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornLLRRWeights(2*ix+ioff,0); + CEW [ix] = amp* EWLLRRWeights(2*ix+ioff,0); + } + } + // LL LL + else { + unsigned int ioff; + if(q1->id()%2==0) { + ioff = abs(q1bar->id())%2==0 ? 0 : 2; + } + else { + ioff = abs(q1bar->id())%2==0 ? 1 : 3; + } + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornLLLLWeights(6*ix+ioff,0); + CEW [ix] = amp* EWLLLLWeights(6*ix+ioff,0); + } + } + } + else { + // RR RR + if(iq2==0) { + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornRRRRWeights(ix,0); + CEW [ix] = amp* EWRRRRWeights(ix,0); + } + } + // RR LL + else { + unsigned int ioff=abs(q1bar->id())%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornRRLLWeights(2*ix+ioff,0); + CEW [ix] = amp* EWRRLLWeights(2*ix+ioff,0); + } + } + } + // square + for(unsigned int ix=0;ix<2;++ix) { + for(unsigned int iy=0;iy<2;++iy) { + bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]); + EWME (ix,iy) += CEW [ix]*conj(CEW [iy]); + } + } + } + } + // colour factors + double born = 2.*real(bornME(0,0))+9.*real(bornME(1,1)); + double EW = 2.*real( EWME(0,0))+9.*real( EWME(1,1)); + return EW/born; +} + +double ElectroWeakReweighter::reweightqqqq() const { + // momenta and invariants + Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum(); + tcPDPtr q1 = subProcess()->incoming().first ->dataPtr(); + Lorentz5Momentum p2 = subProcess()->incoming().second->momentum(); + tcPDPtr q2 = subProcess()->incoming().second->dataPtr(); + Lorentz5Momentum p3 = subProcess()->outgoing()[0] ->momentum(); + tcPDPtr q3 = subProcess()->outgoing()[0] ->dataPtr(); + Lorentz5Momentum p4 = subProcess()->outgoing()[1] ->momentum(); + tcPDPtr q4 = subProcess()->outgoing()[1] ->dataPtr(); + if(q1->id()!=q3->id()) { + swap(q3,q4); + swap(p3,p4); + } + assert(q1->id()==q3->id()); + assert(q2->id()==q4->id()); + 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; + p3.setMass(ZERO); + p3.rescaleRho(); + p4.setMass(ZERO); + p4.rescaleRho(); + // LO and EW corrected matrix element coefficients + boost::numeric::ublas::matrix > + bornLLLLWeights,bornLLRRWeights,bornRRLLWeights,bornRRRRWeights, + EWLLLLWeights,EWLLRRWeights,EWRRLLWeights,EWRRRRWeights; + bool ident = q1->id()==q2->id(); + // LL -> LL + if((q1->id()<=4&& q2->id()<=4)|| (q1->id()==5 && q2->id()==5)) { + if(!ident) { + bornLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,true ,2); + EWLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,false,2); + } + else { + bornLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,true ,2); + EWLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,false,2); + } + } + else if(q1->id()==5 || q2->id()==5) { + bornLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,true ,2); + EWLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,false,2); + } + else + assert(false); + // RR -> LL + if(q1->id()%2==0) { + if(q2->id()<5) { + bornRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,2); + EWRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,2); + } + else { + bornRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,2); + EWRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,2); + } + } + else { + if(q2->id()<5) { + bornRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,2); + EWRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,2); + } + else { + bornRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,2); + EWRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,2); + } + } + // LL -> RR + if(q1->id()<=4) { + if(q2->id()%2!=0) { + bornLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,2); + EWLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,2); + } + else { + bornLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,2); + EWLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,2); + } + } + else { + if(q2->id()%2!=0) { + bornLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,2); + EWLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,2); + } + else { + bornLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,2); + EWLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,2); + } + } + // RR -> RR + if(q1->id()%2==0) { + if(q2->id()%2==0) { + if(ident) { + bornRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,false,2); + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,false,2); + } + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,false,2); + } + } + else { + if(q2->id()%2==0) { + bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,false,2); + } + else { + if(ident) { + bornRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,false,2); + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,false,2); + } + } + } + // extra terms for identical particles + boost::numeric::ublas::matrix > + borntChannelWeights,EWtChannelWeights; + if(ident) { + if(q1->id()%2==0) { + borntChannelWeights = evaluateRunning(EWProcess::QQUU,s,u,t,true ,2); + EWtChannelWeights = evaluateRunning(EWProcess::QQUU,s,u,t,false,2); + } + else if(q1->id()==5) { + borntChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,u,t,true ,2); + EWtChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,u,t,false,2); + } + else { + borntChannelWeights = evaluateRunning(EWProcess::QQDD,s,u,t,true ,2); + EWtChannelWeights = evaluateRunning(EWProcess::QQDD,s,u,t,false,2); + } + } + SpinorWaveFunction q1w(p1,q1,incoming); + SpinorWaveFunction q2w(p2,q2,incoming); + SpinorBarWaveFunction q3w(p3,q3,outgoing); + SpinorBarWaveFunction q4w(p4,q4,outgoing); + boost::numeric::ublas::matrix + bornME = boost::numeric::ublas::zero_matrix(2,2), + EWME = boost::numeric::ublas::zero_matrix(2,2); + for(unsigned int iq1=0;iq1<2;++iq1) { + q1w.reset(iq1); + q3w.reset(iq1); + LorentzVector > current1 = + q1w.dimensionedWave().vectorCurrent(q3w.dimensionedWave()); + for(unsigned int iq2=0;iq2<2;++iq2) { + q2w.reset(iq2); + q4w.reset(iq2); + LorentzVector > current2 = + q2w.dimensionedWave().vectorCurrent(q4w.dimensionedWave()); + complex amp = current1.dot(current2); + vector Cborn(2),CEW(2); + // amplitudes + if(iq1==0) { + // LL + if(iq2==0) { + unsigned int ioff; + if(q1->id()%2==0) { + ioff = q2->id()%2==0 ? 0 : 2; + } + else { + ioff = q2->id()%2==0 ? 1 : 3; + } + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornLLLLWeights(6*ix+ioff,0); + CEW [ix] = amp* EWLLLLWeights(6*ix+ioff,0); + } + } + // LR + else { + unsigned int ioff = q1->id()%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornLLRRWeights(2*ix+ioff,0); + CEW [ix] = amp* EWLLRRWeights(2*ix+ioff,0); + } + } + } + else { + if(iq2==0) { + unsigned int ioff=q2->id()%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornRRLLWeights(2*ix+ioff,0); + CEW [ix] = amp* EWRRLLWeights(2*ix+ioff,0); + } + } + else { + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornRRRRWeights(ix,0); + CEW [ix] = amp* EWRRRRWeights(ix,0); + } + } + } + // square + for(unsigned int ix=0;ix<2;++ix) { + for(unsigned int iy=0;iy<2;++iy) { + bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]); + EWME (ix,iy) += CEW [ix]*conj(CEW [iy]); + } + } + } + } + // extra u-channel pieces if identical flavours + if(ident) { + for(unsigned int iq1=0;iq1<2;++iq1) { + q1w.reset(iq1); + q4w.reset(iq1); + LorentzVector > current1 = + q1w.dimensionedWave().vectorCurrent(q4w.dimensionedWave()); + if(iq1==0) { + q2w.reset(1); + q3w.reset(1); + } + else { + q2w.reset(0); + q3w.reset(0); + } + LorentzVector > current2 = + q2w.dimensionedWave().vectorCurrent(q3w.dimensionedWave()); + complex amp = current1.dot(current2); + vector Cborn(2),CEW(2); + unsigned int ioff = q1->id()%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*borntChannelWeights(2*ix+ioff,0); + CEW [ix] = amp* EWtChannelWeights(2*ix+ioff,0); + } + // square + for(unsigned int ix=0;ix<2;++ix) { + for(unsigned int iy=0;iy<2;++iy) { + bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]); + EWME (ix,iy) += CEW [ix]*conj(CEW [iy]); + } + } + } + } + // colour factors + double born = 2.*real(bornME(0,0))+9.*real(bornME(1,1)); + double EW = 2.*real( EWME(0,0))+9.*real( EWME(1,1)); + return EW/born; +} + +double ElectroWeakReweighter::reweightqbarqbarqbarqbar() const { + // momenta and invariants + Lorentz5Momentum p1 = subProcess()->incoming().first ->momentum(); + tcPDPtr qbar1 = subProcess()->incoming().first ->dataPtr(); + Lorentz5Momentum p2 = subProcess()->incoming().second->momentum(); + tcPDPtr qbar2 = subProcess()->incoming().second->dataPtr(); + Lorentz5Momentum p3 = subProcess()->outgoing()[0] ->momentum(); + tcPDPtr qbar3 = subProcess()->outgoing()[0] ->dataPtr(); + Lorentz5Momentum p4 = subProcess()->outgoing()[1] ->momentum(); + tcPDPtr qbar4 = subProcess()->outgoing()[1] ->dataPtr(); + if(qbar1->id()!=qbar3->id()) { + swap(qbar3,qbar4); + swap(p3,p4); + } + assert(qbar1->id()==qbar3->id()); + assert(qbar2->id()==qbar4->id()); + 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; + p3.setMass(ZERO); + p3.rescaleRho(); + p4.setMass(ZERO); + p4.rescaleRho(); + // LO and EW corrected matrix element coefficients + boost::numeric::ublas::matrix > + bornLLLLWeights,bornLLRRWeights,bornRRLLWeights,bornRRRRWeights, + EWLLLLWeights,EWLLRRWeights,EWRRLLWeights,EWRRRRWeights; + bool ident = qbar1->id()==qbar2->id(); + // LL -> LL + if((abs(qbar1->id())<=4 && abs(qbar2->id())<=4) || + (abs(qbar1->id())==5 && abs(qbar2->id())==5)) { + if(!ident) { + bornLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,true ,2); + EWLLLLWeights = evaluateRunning(EWProcess::QQQQ,s,t,u,false,2); + } + else { + bornLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,true ,2); + EWLLLLWeights = evaluateRunning(EWProcess::QQQQiden,s,t,u,false,2); + } + } + else if(abs(qbar1->id())==5 || abs(qbar2->id())==5) { + bornLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,true ,2); + EWLLLLWeights = evaluateRunning(EWProcess::QtQtQQ,s,t,u,false,2); + } + else + assert(false); + // RR -> LL + if(abs(qbar1->id())%2==0) { + if(abs(qbar2->id())<5) { + bornRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,2); + EWRRLLWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,2); + } + else { + bornRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,2); + EWRRLLWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,2); + } + } + else { + if(abs(qbar2->id())<5) { + bornRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,2); + EWRRLLWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,2); + } + else { + bornRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,2); + EWRRLLWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,2); + } + } + // LL -> RR + if(abs(qbar1->id())<=4) { + if(abs(qbar2->id())%2!=0) { + bornLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,true ,2); + EWLLRRWeights = evaluateRunning(EWProcess::QQDD,s,t,u,false,2); + } + else { + bornLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,true ,2); + EWLLRRWeights = evaluateRunning(EWProcess::QQUU,s,t,u,false,2); + } + } + else { + if(abs(qbar2->id())%2!=0) { + bornLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,true ,2); + EWLLRRWeights = evaluateRunning(EWProcess::QtQtDD,s,t,u,false,2); + } + else { + bornLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,true ,2); + EWLLRRWeights = evaluateRunning(EWProcess::QtQtUU,s,t,u,false,2); + } + } + // RR -> RR + if(abs(qbar1->id())%2==0) { + if(abs(qbar2->id())%2==0) { + if(ident) { + bornRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::UUUUiden,s,t,u,false,2); + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::UUUU,s,t,u,false,2); + } + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,false,2); + } + } + else { + if(abs(qbar2->id())%2==0) { + bornRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::UUDD,s,t,u,false,2); + } + else { + if(ident) { + bornRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::DDDDiden,s,t,u,false,2); + } + else { + bornRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,true ,2); + EWRRRRWeights = evaluateRunning(EWProcess::DDDD,s,t,u,false,2); + } + } + } + // extra terms for identical particles + boost::numeric::ublas::matrix > + borntChannelWeights,EWtChannelWeights; + if(ident) { + if(abs(qbar1->id())%2==0) { + borntChannelWeights = evaluateRunning(EWProcess::QQUU,s,u,t,true ,2); + EWtChannelWeights = evaluateRunning(EWProcess::QQUU,s,u,t,false,2); + } + else if(abs(qbar1->id())==5) { + borntChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,u,t,true ,2); + EWtChannelWeights = evaluateRunning(EWProcess::QtQtDD,s,u,t,false,2); + } + else { + borntChannelWeights = evaluateRunning(EWProcess::QQDD,s,u,t,true ,2); + EWtChannelWeights = evaluateRunning(EWProcess::QQDD,s,u,t,false,2); + } + } + SpinorBarWaveFunction qbar1w(p1,qbar1,incoming); + SpinorBarWaveFunction qbar2w(p2,qbar2,incoming); + SpinorWaveFunction qbar3w(p3,qbar3,outgoing); + SpinorWaveFunction qbar4w(p4,qbar4,outgoing); + boost::numeric::ublas::matrix + bornME = boost::numeric::ublas::zero_matrix(2,2), + EWME = boost::numeric::ublas::zero_matrix(2,2); + for(unsigned int iq1=0;iq1<2;++iq1) { + qbar1w.reset(iq1); + qbar3w.reset(iq1); + LorentzVector > current1 = + qbar3w.dimensionedWave().vectorCurrent(qbar1w.dimensionedWave()); + for(unsigned int iq2=0;iq2<2;++iq2) { + qbar2w.reset(iq2); + qbar4w.reset(iq2); + LorentzVector > current2 = + qbar4w.dimensionedWave().vectorCurrent(qbar2w.dimensionedWave()); + complex amp = current1.dot(current2); + vector Cborn(2),CEW(2); + // amplitudes + if(iq1==1) { + // LL + if(iq2==1) { + unsigned int ioff; + if(abs(qbar1->id())%2==0) { + ioff = abs(qbar2->id())%2==0 ? 0 : 2; + } + else { + ioff = abs(qbar2->id())%2==0 ? 1 : 3; + } + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornLLLLWeights(6*ix+ioff,0); + CEW [ix] = amp* EWLLLLWeights(6*ix+ioff,0); + } + } + // LR + else { + unsigned int ioff = abs(qbar1->id())%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornLLRRWeights(2*ix+ioff,0); + CEW [ix] = amp* EWLLRRWeights(2*ix+ioff,0); + } + } + } + else { + if(iq2==1) { + unsigned int ioff=abs(qbar2->id())%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornRRLLWeights(2*ix+ioff,0); + CEW [ix] = amp* EWRRLLWeights(2*ix+ioff,0); + } + } + else { + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*bornRRRRWeights(ix,0); + CEW [ix] = amp* EWRRRRWeights(ix,0); + } + } + } + // square + for(unsigned int ix=0;ix<2;++ix) { + for(unsigned int iy=0;iy<2;++iy) { + bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]); + EWME (ix,iy) += CEW [ix]*conj(CEW [iy]); + } + } + } + } + // extra u-channel pieces if identical flavours + if(ident) { + for(unsigned int iq1=0;iq1<2;++iq1) { + qbar1w.reset(iq1); + qbar4w.reset(iq1); + LorentzVector > current1 = + qbar4w.dimensionedWave().vectorCurrent(qbar1w.dimensionedWave()); + if(iq1==0) { + qbar2w.reset(1); + qbar3w.reset(1); + } + else { + qbar2w.reset(0); + qbar3w.reset(0); + } + LorentzVector > current2 = + qbar3w.dimensionedWave().vectorCurrent(qbar2w.dimensionedWave()); + complex amp = current1.dot(current2); + vector Cborn(2),CEW(2); + unsigned int ioff = abs(qbar1->id())%2==0 ? 0 : 1; + for(unsigned int ix=0;ix<2;++ix) { + Cborn[ix] = amp*borntChannelWeights(2*ix+ioff,0); + CEW [ix] = amp* EWtChannelWeights(2*ix+ioff,0); + } + // square + for(unsigned int ix=0;ix<2;++ix) { + for(unsigned int iy=0;iy<2;++iy) { + bornME(ix,iy) += Cborn[ix]*conj(Cborn[iy]); + EWME (ix,iy) += CEW [ix]*conj(CEW [iy]); + } + } + } + } + // colour factors + double born = 2.*real(bornME(0,0))+9.*real(bornME(1,1)); + double EW = 2.*real( EWME(0,0))+9.*real( EWME(1,1)); + 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,175 +1,195 @@ // -*- 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; + + /** + * Reweight \f$q\bar{q}\to q'\bar{q'}\f$ (s-channel) + */ + double reweightqqbarqqbarS() const; + + /** + * Reweight \f$q\bar{q}\to q'\bar{q'}\f$ (t-channel) + */ + double reweightqqbarqqbarT() const; + + /** + * Reweight \f$qq \to qq\f$ + */ + double reweightqqqq() const; + + /** + * Reweight \f$\bar{q}\bar{q} \to \bar{q}\bar{q}\f$ + */ + double reweightqbarqbarqbarqbar() 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, 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,407 +1,495 @@ // -*- 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 Gamma2ST(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./2.*(T-I*pi) + U -2.*T; + output(1,1) = 3./2.*(T-I*pi); + output(0,1) = -2.0*U; + output(1,0) = -(3.0/8.0)*U; + return output; + } + + inline boost::numeric::ublas::matrix Gamma2SU(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./2.*(U-I*pi) + T - 2.*U; + output(1,1) = 3./2.*(U-I*pi); + output(0,1) = 2.0*T; + output(1,0) = (3.0/8.0)*T; + 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() { 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*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 boost::numeric::ublas::matrix Gamma2SingletSU(Complex U) { + 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*(U-I*pi); + return output; + } + inline Complex Gamma1(double hypercharge) { Complex I(0,1.0); 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 Gamma1SU(double hypercharge,Complex U) { + Complex I(0,1.0); + return (U-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 Gamma1SU(double y1, double y2, Complex T, Complex U) { + Complex I(0,1.0); + return (U-I*Constants::pi)*(y1*y1+y2*y2) + 2.0*y1*y2*T; + } + 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 Complex Gamma1SU(double y1, double y2, double y3, double y4, + Complex T, Complex U) { + Complex I(0,1.0); + return (U-I*Constants::pi)*(y1*y1+y2*y2+y3*y3+y4*y4)/2.0 + + (y1*y4+y2*y3)*(T-U) + (y1*y3+y2*y4)*U; + } 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 Gamma3ST(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./3.*(T-I*pi); + output(1,1) = 8./3.*(T-I*pi); + output(0,0) += -3.*T + 2./3.*U; + output(0,1) = -2.*U; + output(1,0) = -4./9.*U; + return output; + } + + inline boost::numeric::ublas::matrix Gamma3SU(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./3.*(U-I*pi); + output(1,1) = 8./3.*(U-I*pi); + output(0,0) += 7./3.*T -3.*U; + output(0,1) = 2.*T; + output(1,0) = 4./9.*T; + 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 Gamma3gSU(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) = -T; + output(1,1) = 3./2.*(T-2.*U); + output(1,2) = -3./2.*T; + output(2,0) = -2.0*T; + output(2,1) =-5./6.*T; + output(2,2) = 3./2.*(T-2.*U); + for (unsigned int i=0; i<3; i++) { + output(i,i) += 13./3.*(U-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,1146 +1,1321 @@ // -*- 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, 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); matrix G1, G2, G3; unsigned int numBrokenGauge; switch (process) { case QQQQ: case QQQQiden: case QtQtQQ: { - assert(iswap==0); numBrokenGauge = 12; G1 = zero_matrix(numBrokenGauge,numBrokenGauge); G2 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); - matrix gam3 = Gamma3(U,T); + matrix gam3; + if(iswap==0) { + gam3 = Gamma3(U,T); + G1(0,0) = G1(6,6) = Gamma1(2.0/3.0,2.0/3.0,2.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(7,7) = Gamma1(-1.0/3.0,-1.0/3.0,2.0/3.0,2.0/3.0,T,U); + G1(2,2) = G1(8,8) = Gamma1(2.0/3.0,2.0/3.0,-1.0/3.0,-1.0/3.0,T,U); + G1(3,3) = G1(9,9) = Gamma1(-1.0/3.0,-1.0/3.0,-1.0/3.0,-1.0/3.0,T,U); + G1(4,4) = G1(10,10) = Gamma1(-1.0/3.0,2.0/3.0,2.0/3.0,-1.0/3.0,T,U); + G1(5,5) = G1(11,11) = Gamma1(2.0/3.0,-1.0/3.0,-1.0/3.0,2.0/3.0,T,U); + } + else if(iswap==1) { + gam3 = Gamma3ST(U,T); + G1(0,0) = G1(6,6) = Gamma1ST(2.0/3.0,2.0/3.0,2.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(7,7) = Gamma1ST(-1.0/3.0,-1.0/3.0,2.0/3.0,2.0/3.0,T,U); + G1(2,2) = G1(8,8) = Gamma1ST(2.0/3.0,2.0/3.0,-1.0/3.0,-1.0/3.0,T,U); + G1(3,3) = G1(9,9) = Gamma1ST(-1.0/3.0,-1.0/3.0,-1.0/3.0,-1.0/3.0,T,U); + G1(4,4) = G1(10,10) = Gamma1ST(-1.0/3.0,2.0/3.0,2.0/3.0,-1.0/3.0,T,U); + G1(5,5) = G1(11,11) = Gamma1ST(2.0/3.0,-1.0/3.0,-1.0/3.0,2.0/3.0,T,U); + } + else if(iswap==2) { + gam3 = Gamma3SU(U,T); + G1(0,0) = G1(6,6) = Gamma1SU( 2.0/3.0,2.0/3.0,2.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(7,7) = Gamma1SU(-1.0/3.0,-1.0/3.0,2.0/3.0,2.0/3.0,T,U); + G1(2,2) = G1(8,8) = Gamma1SU( 2.0/3.0,2.0/3.0,-1.0/3.0,-1.0/3.0,T,U); + G1(3,3) = G1(9,9) = Gamma1SU(-1.0/3.0,-1.0/3.0,-1.0/3.0,-1.0/3.0,T,U); + G1(4,4) = G1(10,10) = Gamma1SU(-1.0/3.0,2.0/3.0,2.0/3.0,-1.0/3.0,T,U); + G1(5,5) = G1(11,11) = Gamma1SU( 2.0/3.0,-1.0/3.0,-1.0/3.0,2.0/3.0,T,U); + } + else + assert(false); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); G2 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); - matrix gam3 = Gamma3(U,T); + matrix gam3; + if(iswap==0) { + gam3 = Gamma3(U,T); + G1(0,0) = G1(2,2) = Gamma1(2.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(3,3) = Gamma1(2.0/3.0,-1.0/3.0,T,U); + } + else if(iswap==1) { + gam3 = Gamma3ST(U,T); + G1(0,0) = G1(2,2) = Gamma1ST(2.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(3,3) = Gamma1ST(2.0/3.0,-1.0/3.0,T,U); + } + else if(iswap==2) { + gam3 = Gamma3SU(U,T); + G1(0,0) = G1(2,2) = Gamma1SU(2.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(3,3) = Gamma1SU(2.0/3.0,-1.0/3.0,T,U); + } + else + assert(false); + for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); G2 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); - matrix gam3 = Gamma3(U,T); + matrix gam3; + if(iswap==0) { + gam3 = Gamma3(U,T); + G1(0,0) = G1(2,2) = Gamma1(-1.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(3,3) = Gamma1(-1.0/3.0,-1.0/3.0,T,U); + } + else if(iswap==1) { + gam3 = Gamma3ST(U,T); + G1(0,0) = G1(2,2) = Gamma1ST(-1.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(3,3) = Gamma1ST(-1.0/3.0,-1.0/3.0,T,U); + } + else if(iswap==2) { + gam3 = Gamma3SU(U,T); + G1(0,0) = G1(2,2) = Gamma1SU(-1.0/3.0,2.0/3.0,T,U); + G1(1,1) = G1(3,3) = Gamma1SU(-1.0/3.0,-1.0/3.0,T,U); + } + else + assert(false); for (unsigned int i=0; i(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 = 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 = 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); + if(iswap==0) { + G3 = Gamma3(U,T); + G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,2.0/3.0,T,U); + } + else if(iswap==1) { + G3 = Gamma3ST(U,T); + G1(0,0) = G1(1,1) = Gamma1ST(2.0/3.0,2.0/3.0,T,U); + } + else if(iswap==2) { + G3 = Gamma3SU(U,T); + G1(0,0) = G1(1,1) = Gamma1SU(2.0/3.0,2.0/3.0,T,U); + } + else + assert(false); break; case UUDD: case tRtRDD: - assert(iswap==0); numBrokenGauge = 2; 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); + if(iswap==0) { + G3 = Gamma3(U,T); + G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,2.0/3.0,T,U); + } + else if(iswap==1) { + G3 = Gamma3ST(U,T); + G1(0,0) = G1(1,1) = Gamma1ST(-1.0/3.0,2.0/3.0,T,U); + } + else if(iswap==2) { + G3 = Gamma3SU(U,T); + G1(0,0) = G1(1,1) = Gamma1SU(-1.0/3.0,2.0/3.0,T,U); + } + else + assert(false); break; case UULL: assert(iswap==0); numBrokenGauge = 2; 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 = 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 = 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); + if(iswap==0) { + G3 = Gamma3(U,T); + G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,T,U); + } + else if(iswap==1) { + G3 = Gamma3ST(U,T); + G1(0,0) = G1(1,1) = Gamma1ST(-1.0/3.0,-1.0/3.0,T,U); + } + else if(iswap==2) { + G3 = Gamma3SU(U,T); + G1(0,0) = G1(1,1) = Gamma1SU(-1.0/3.0,-1.0/3.0,T,U); + } + else + assert(false); break; case DDLL: assert(iswap==0); numBrokenGauge = 2; 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 = 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 = 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 = 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 = 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 = 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 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); G2 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); G2 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); for (unsigned int i=0; i(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 if(iswap==2) { + gam3g = Gamma3gSU(U,T); + gam1a = Gamma1SU( 2./3.,0.,T,U); + gam1b = Gamma1SU(-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) = gam1a; G1(3,3) = G1(4,4) = G1(5,5) = gam1b; } break; case LLWW: assert(iswap==0); numBrokenGauge = 20; 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 = 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 = 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 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(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 = 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 if(iswap==2) { + G3 = Gamma3gSU(U,T); + gam1 = Gamma1SU(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 = 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 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); Complex gam3s = Gamma3Singlet()(0,0); for (unsigned int i=0; i(numBrokenGauge,numBrokenGauge); G2 = zero_matrix(numBrokenGauge,numBrokenGauge); G3 = zero_matrix(numBrokenGauge,numBrokenGauge); for (unsigned int i=0; i(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 = 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 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, 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); matrix G1,G2,G3; unsigned int numGauge; switch (process) { case QQQQ: case QQQQiden: case QtQtQQ: { - assert(iswap==0); numGauge = 4; G1 = zero_matrix(numGauge,numGauge); G2 = zero_matrix(numGauge,numGauge); G3 = zero_matrix(numGauge,numGauge); - matrix gam3 = Gamma3(U,T); + matrix gam3,gam2; + if(iswap==0) { + gam3 = Gamma3(U,T); + gam2 = Gamma2(U,T); + G1(0,0) = G1(1,1) = G1(2,2) = G1(3,3) = Gamma1(1.0/6.0,1.0/6.0,T,U); + } + else if(iswap==1) { + gam3 = Gamma3ST(U,T); + gam2 = Gamma2ST(U,T); + G1(0,0) = G1(1,1) = G1(2,2) = G1(3,3) = Gamma1ST(1.0/6.0,1.0/6.0,T,U); + } + else if(iswap==2) { + gam3 = Gamma3SU(U,T); + gam2 = Gamma2SU(U,T); + G1(0,0) = G1(1,1) = G1(2,2) = G1(3,3) = Gamma1SU(1.0/6.0,1.0/6.0,T,U); + } + else + assert(false); 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); - 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 = 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); + if(iswap==0) { + G3 = Gamma3(U,T); + G2 = Gamma2Singlet(); + G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,1.0/6.0,T,U); + } + else if(iswap==1) { + G3 = Gamma3ST(U,T); + G2 = Gamma2SingletST(T); + G1(0,0) = G1(1,1) = Gamma1ST(2.0/3.0,1.0/6.0,T,U); + } + else if(iswap==2) { + G3 = Gamma3SU(U,T); + G2 = Gamma2SingletSU(U); + G1(0,0) = G1(1,1) = Gamma1SU(2.0/3.0,1.0/6.0,T,U); + } + else + assert(false); break; case QQDD: case QtQtDD: - assert(iswap==0); numGauge = 2; 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); + if(iswap==0) { + G3 = Gamma3(U,T); + G2 = Gamma2Singlet(); + G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,1.0/6.0,T,U); + } + else if(iswap==1) { + G3 = Gamma3ST(U,T); + G2 = Gamma2SingletST(T); + G1(0,0) = G1(1,1) = Gamma1ST(-1.0/3.0,1.0/6.0,T,U); + } + else if(iswap==2) { + G3 = Gamma3SU(U,T); + G2 = Gamma2SingletSU(U); + G1(0,0) = G1(1,1) = Gamma1SU(-1.0/3.0,1.0/6.0,T,U); + } + else + assert(false); break; case QQLL: assert(iswap==0); numGauge = 2; 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 = 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 = 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); + if(iswap==0) { + G3 = Gamma3(U,T); + G1(0,0) = G1(1,1) = Gamma1(2.0/3.0,2.0/3.0,T,U); + } + else if(iswap==1) { + G3 = Gamma3ST(U,T); + G1(0,0) = G1(1,1) = Gamma1ST(2.0/3.0,2.0/3.0,T,U); + } + else if(iswap==2) { + G3 = Gamma3SU(U,T); + G1(0,0) = G1(1,1) = Gamma1SU(2.0/3.0,2.0/3.0,T,U); + } + else + assert(false); break; case UUDD: case tRtRDD: - assert(iswap==0); numGauge = 2; 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); + if(iswap==0) { + G3 = Gamma3(U,T); + G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,2.0/3.0,T,U); + } + else if(iswap==1) { + G3 = Gamma3ST(U,T); + G1(0,0) = G1(1,1) = Gamma1ST(-1.0/3.0,2.0/3.0,T,U); + } + else if(iswap==2) { + G3 = Gamma3SU(U,T); + G1(0,0) = G1(1,1) = Gamma1SU(-1.0/3.0,2.0/3.0,T,U); + } + else + assert(false); break; case UULL: assert(iswap==0); numGauge = 1; 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 = 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 = 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); + if(iswap==0) { + G3 = Gamma3(U,T); + G1(0,0) = G1(1,1) = Gamma1(-1.0/3.0,-1.0/3.0,T,U); + } + else if(iswap==1) { + G3 = Gamma3ST(U,T); + G1(0,0) = G1(1,1) = Gamma1ST(-1.0/3.0,-1.0/3.0,T,U); + } + else if(iswap==2) { + G3 = Gamma3SU(U,T); + G1(0,0) = G1(1,1) = Gamma1SU(-1.0/3.0,-1.0/3.0,T,U); + } + else + assert(false); break; case DDLL: assert(iswap==0); numGauge = 1; 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 if(iswap==2) { + G3 = Gamma3gSU(U,T); + gam2s = Gamma2SingletSU(U)(0,0); + gam1 = Gamma1SU(1.0/6.0,U); + } else assert(false); for (unsigned int i=0; i<3; i++) { G2(i,i) = gam2s; G1(i,i) = gam1; } } break; case LLWW: assert(iswap==0); numGauge = 5; 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 = 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 = 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 = 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 = 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 = 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 if(iswap==2) { + G3 = Gamma3gSU(U,T); + gam1 = Gamma1SU(Y,U); + } else assert(false); for (unsigned int i=0; i<3; i++) { G1(i,i) = gam1; } } break; case DDBB: assert(iswap==0); numGauge = 1; 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 = 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 = 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 EEBB: assert(iswap==0); numGauge = 1; 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 = 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); } }