diff --git a/MatrixElement/EW/CollinearSudakov.cc b/MatrixElement/EW/CollinearSudakov.cc --- a/MatrixElement/EW/CollinearSudakov.cc +++ b/MatrixElement/EW/CollinearSudakov.cc @@ -1,1192 +1,1730 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the CollinearSudakov class. // #include "CollinearSudakov.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 "ElectroWeakReweighter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; namespace { void DuplicateColumn0(boost::numeric::ublas::matrix &orig) { for (unsigned int i=0; i> K_ORDER_ >> B_ORDER_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigCollinearSudakov("Herwig::CollinearSudakov", "HwMEEW.so"); void CollinearSudakov::Init() { static ClassDocumentation documentation ("The CollinearSudakov class implements the collinear evolution"); } Complex CollinearSudakov::highScaleIntegral(bool SU3, bool SU2, double Y, Energy2 s, Energy mu_h, Energy mu_l, bool fermion, bool longitudinal, double yukFactor) { SU3_ = SU3; SU2_ = SU2; Y_ = Y; s_ = s; fermion_ = fermion; longitudinal_ = longitudinal; yukFactor_ = yukFactor; // perform the integral Complex result; high_ = true; // real piece real_ = true; result.real(integrator_.value(*this,mu_h,mu_l)); // imaginary piece real_ = false; result.imag(integrator_.value(*this,mu_h,mu_l)); // return the answer return exp(result); } Complex CollinearSudakov::lowScaleIntegral(bool SU3, double Q, Energy2 s, Energy mu_h, Energy mu_l, bool fermion, double boostFactor) { SU3_ = SU3; Q_ = Q; s_ = s; fermion_ = fermion; boostFactor_ = boostFactor; // perform the integral Complex result; high_ = false; // real piece real_ = true; result.real(integrator_.value(*this,mu_h,mu_l)); // imaginary piece real_ = false; result.imag(integrator_.value(*this,mu_h,mu_l)); // return the answer return exp(result); } InvEnergy CollinearSudakov::highScaleIntegrand(Energy mu) const { using namespace GroupInvariants; using Constants::pi; Complex gamma = 0.0; // Include K-factor Contributions (Cusps): GaugeContributions cusp = cuspContributions(mu,K_ORDER_,high_); // Include B-factors (B1 for U1, B2 for SU2, B3 for SU3): GaugeContributions nonCusp = BContributions(mu,B_ORDER_,fermion_,longitudinal_); // common log Complex plog = PlusLog(s_/sqr(mu)); // SU(3) if(SU3_) { if(fermion_) gamma += C_F(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3; else gamma += (C_A(3))*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3; } // SU(2) if(SU2_) { if (fermion_ || longitudinal_ ) gamma += (C_F(2))*cusp.SU2*0.5*plog + 0.5*nonCusp.SU2; else gamma += (C_A(2))*cusp.SU2*0.5*plog + 0.5*nonCusp.SU2; } if (fermion_ || longitudinal_ ) { gamma += sqr(Y_)*(cusp.U1*0.5*plog + 0.5*nonCusp.U1); } else { // U(1) Gauge boson if (!SU3_ && !SU2_ && abs(Y_)<0.001) { gamma += 0.5*nonCusp.U1; } } // top Yukawa piece double y_t = ElectroWeakReweighter::coupling()->y_t(mu); gamma += yukFactor_*sqr(y_t)/(16.0*sqr(pi)); // return the answer return real_ ? gamma.real()/mu : gamma.imag()/mu; } InvEnergy CollinearSudakov::lowScaleIntegrand(Energy mu) const { using namespace GroupInvariants; using Constants::pi; Complex gamma = 0.0; // Include K-factor Contributions (Cusps): GaugeContributions cusp = cuspContributions(mu,K_ORDER_,false); // Include B-factors (B1 for U1, B2 for SU2, B3 for SU3): GaugeContributions nonCusp = BContributionsLow(mu,B_ORDER_,fermion_,boostFactor_); // common log Complex plog = PlusLog(s_/sqr(mu)); Complex blog(0.); if(boostFactor_ >0.001) blog = PlusLog(4.0*boostFactor_); // SU(3) if (SU3_) { if (fermion_) { // not a bHQET top quark field if (abs(boostFactor_)<0.001) gamma += C_F(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3; else gamma += C_F(3)*cusp.SU3*0.5*blog + 0.5*nonCusp.SU3; } else { gamma += C_A(3)*cusp.SU3*0.5*plog + 0.5*nonCusp.SU3; } } // fermions if (fermion_) { // not a bHQET top quark field if (abs(boostFactor_)<0.001) gamma += sqr(Q_)*(cusp.U1*0.5*plog + 0.5*nonCusp.U1); else gamma += sqr(Q_)*(cusp.U1*0.5*blog + 0.5*nonCusp.U1); } else { // i.e., not a fermion, not a bHQ, not a gluon => photon if (abs(boostFactor_)<0.001 && !SU3_) gamma += 0.5*nonCusp.U1; // i.e., W treated as a bHQET field else if (abs(boostFactor_)>0.001) gamma += sqr(Q_)*(cusp.U1*0.5*blog + 0.5*nonCusp.U1); } // return the answer return real_ ? gamma.real()/mu : gamma.imag()/mu; } void CollinearSudakov::evaluateHighScale(Energy highScale, Energy EWScale, Energy2 S) { double yCoeffQt(0.), yCoefftR(0.), yCoeffPhi(0.); if (K_ORDER_ >= 1) { yCoeffQt = 0.5; yCoefftR = 1.0; yCoeffPhi = 3.0; } highColW_ = highScaleIntegral(false,true ,0.0 ,S,highScale,EWScale,false,false,0.0); highColB_ = highScaleIntegral(false,false,0.0 ,S,highScale,EWScale,false,false,0.0); highColG_ = highScaleIntegral(true ,false,0.0 ,S,highScale,EWScale,false,false,0.0); highColQ_ = highScaleIntegral(true ,true ,1./6. ,S,highScale,EWScale,true,false,0.0); highColQt_ = highScaleIntegral(true ,true ,1./6. ,S,highScale,EWScale,true,false,yCoeffQt); highColU_ = highScaleIntegral(true ,false,2./3. ,S,highScale,EWScale,true,false,0.0); highColtR_ = highScaleIntegral(true ,false,2./3. ,S,highScale,EWScale,true,false,yCoefftR); highColD_ = highScaleIntegral(true ,false,-1./3.,S,highScale,EWScale,true,false,0.0); highColL_ = highScaleIntegral(false,true ,-1./2.,S,highScale,EWScale,true,false,0.0); highColE_ = highScaleIntegral(false,false,-1. ,S,highScale,EWScale,true,false,0.0); highColPhi_ = highScaleIntegral(false,true ,1./2. ,S,highScale,EWScale,false,true,yCoeffPhi); } void CollinearSudakov::evaluateLowScale(Energy EWScale, Energy lowScale, Energy2 S) { lowColW_ = lowScaleIntegral(false,1. ,S,EWScale,lowScale,false, S/sqr(2.*ElectroWeakReweighter::coupling()->mW())); lowColA_ = lowScaleIntegral(false,0. ,S,EWScale,lowScale,false,0.0); lowColG_ = lowScaleIntegral(true ,0. ,S,EWScale,lowScale,false,0.0); lowColU_ = lowScaleIntegral(true ,2./3. ,S,EWScale,lowScale,true,0.0); lowColt_ = lowScaleIntegral(true ,2./3. ,S,EWScale,lowScale,true, S/sqr(2.*ElectroWeakReweighter::coupling()->mT())); lowColD_ = lowScaleIntegral(true ,-1./3.,S,EWScale,lowScale,true,0.0); lowColE_ = lowScaleIntegral(false,-1. ,S,EWScale,lowScale,true,0.0); } void CollinearSudakov::evaluateMatching(Energy EWScale,Energy2 s) { using Constants::pi; using GroupInvariants::PlusLog; static const Complex I(0,1.0); // wave function corrections WaveFunctionCorrections WFC = waveFunctionCorrections(EWScale); double aS = ElectroWeakReweighter::coupling()->a3(EWScale); double aEM = ElectroWeakReweighter::coupling()->aEM(EWScale); double aW = ElectroWeakReweighter::coupling()->aW(EWScale); double aZ = ElectroWeakReweighter::coupling()->aZ(EWScale); double cW2 = ElectroWeakReweighter::coupling()->Cos2thW(EWScale); double sW2 = ElectroWeakReweighter::coupling()->Sin2thW(EWScale); Energy mW = ElectroWeakReweighter::coupling()->mW(); Energy mZ = ElectroWeakReweighter::coupling()->mZ(); Energy mH = ElectroWeakReweighter::coupling()->mH(); Energy mT = ElectroWeakReweighter::coupling()->mT(); double gPHI = ElectroWeakReweighter::coupling()->g_phiPlus(EWScale); double y_t = ElectroWeakReweighter::coupling()->y_t(EWScale); double lz = log(mZ/EWScale); double lw = log(mW/EWScale); complex F_W = 4.0*lw*0.5*PlusLog(s/EWScale/EWScale)-2.0*lw*lw-2.0*lw-5.0*pi*pi/12.0+1.0; complex F_Z = 4.0*lz*0.5*PlusLog(s/EWScale/EWScale)-2.0*lz*lz-2.0*lz-5.0*pi*pi/12.0+1.0; // Taken from Manohar... along with his formulae for F_tL, F_tR, and F_bL (for 3rd/1st Gen. Wavefunction Differences) complex W1 = WFC.fFW0-0.5*WFC.aW0-0.5*WFC.cW0; complex W2 = WFC.fF0W-0.5*WFC.a0W; complex U1 = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*(WFC.fFZZ-0.5*WFC.aZZ) - 0.5*WFC.cZZ*(ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale) + ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)) + ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*WFC.bZZ; complex U2 = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*(WFC.fFZZ-0.5*WFC.aZZ) - 0.5*WFC.cZZ*(ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale) + ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)) + ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*WFC.bZZ; complex HtL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-0.5*log(mH/EWScale)-0.5*lz+ 0.5*WFC.atHH+0.5*WFC.atZZ+WFC.ctHH+WFC.ctZZ+ WFC.ctW0-WFC.btHH+WFC.btZZ); complex HtR = HtL-0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.atW0); complex HbL = -0.5*y_t*y_t/(16.0*pi*pi)*(0.25-lw+WFC.at0W); complex F_tL = (4.0/3.0*aS+4.0/9.0*aEM)/(4.0*pi)*(2.0*log(mT/EWScale)*log(mT/EWScale)-log(mT/EWScale)+ pi*pi/12.0+2.0) + HtL + aW/(4.0*pi)*0.5*W1 + aZ/(4.0*pi)*U1; complex F_tR = (4.0/3.0*aS+4.0/9.0*aEM)/(4.0*pi)*(2.0*log(mT/EWScale)*log(mT/EWScale)-log(mT/EWScale)+ pi*pi/12.0+2.0) + HtR - aW/(4.0*pi)*0.25*WFC.cW0 + aZ/(4.0*pi)*U2; complex F_bL = HbL + aW/(4.0*pi)*0.5*W2; Complex Dw = CollinearDw(s,EWScale); Complex Dz = CollinearDz(s,EWScale); complex D_C_UL = ElectroWeakReweighter::coupling()->g_Lu(EWScale)*ElectroWeakReweighter::coupling()->g_Lu(EWScale)*Dz + 0.5*Dw; complex D_C_DL = ElectroWeakReweighter::coupling()->g_Ld(EWScale)*ElectroWeakReweighter::coupling()->g_Ld(EWScale)*Dz + 0.5*Dw; complex D_C_UR = ElectroWeakReweighter::coupling()->g_Ru(EWScale)*ElectroWeakReweighter::coupling()->g_Ru(EWScale)*Dz; complex D_C_DR = ElectroWeakReweighter::coupling()->g_Rd(EWScale)*ElectroWeakReweighter::coupling()->g_Rd(EWScale)*Dz; complex D_C_nuL = ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*ElectroWeakReweighter::coupling()->g_Lnu(EWScale)*Dz + 0.5*Dw; complex D_C_EL = ElectroWeakReweighter::coupling()->g_Le(EWScale)*ElectroWeakReweighter::coupling()->g_Le(EWScale)*Dz + 0.5*Dw; complex D_C_ER = ElectroWeakReweighter::coupling()->g_Re(EWScale)*ElectroWeakReweighter::coupling()->g_Re(EWScale)*Dz; complex D_C_WW = aW/(4.0*pi)*cW2*(F_Z+WFC.fsWZWZ) + aW/(4.0*pi)*cW2*(F_W+WFC.fs1ZW) + aW/(4.0*pi)*sW2*(F_W+WFC.fs10) + aW/(4.0*pi)*sW2*(2.0*lw*lw- 2.0*lw+pi*pi/12.0+2.0) + 0.5*WFC.RW; complex D_C_WZ = aW/(4.0*pi)*2.0*(F_W+WFC.fsZW1) + 0.5*WFC.RZ + sqrt(sW2/cW2)*WFC.RAtoZ; complex D_C_WA = aW/(4.0*pi)*2.0*(F_W+WFC.fs01) + 0.5*WFC.RA + sqrt(cW2/sW2)*WFC.RZtoA; complex D_C_BZ = 0.5*WFC.RZ - sqrt(cW2/sW2)*WFC.RAtoZ; complex D_C_BA = 0.5*WFC.RA - sqrt(sW2/cW2)*WFC.RZtoA; // The WFC.RW and WFC.RZ are used on purpose (instead of WFC.RPhi and WFC.RPhi3, respectively): complex D_C_PhiW = aW/(4.0*pi)*0.25*(F_W+WFC.fs1HW) + aW/(4.0*pi)*0.25*(F_W+WFC.fs1ZW) + aZ/(4.0*pi)*gPHI*gPHI*(F_Z+WFC.fsWZWZ) + aW/(4.0*pi)*sW2*(2.0*lw*lw-2.0*lw+pi*pi/12.0+2.0) + 0.5*WFC.RW + WFC.EW; complex D_C_PhiZ = aW/(4.0*pi)*0.5*(F_W+WFC.fsZW1) + aZ/(4.0*pi)*0.25*(F_Z+WFC.fs1HZ) + 0.5*WFC.RZ + WFC.EZ; complex D_C_PhiH = aW/(4.0*pi)*0.5*(F_W+WFC.fsHW1) + aZ/(4.0*pi)*0.25*(F_Z+WFC.fsHZ1) + 0.5*WFC.RH; complex D_C_GG = aS/(4.0*pi)*2.0/3.0*log(mT/EWScale); ULcollinearCorr_ = exp(D_C_UL); DLcollinearCorr_ = exp(D_C_DL); URcollinearCorr_ = exp(D_C_UR); DRcollinearCorr_ = exp(D_C_DR); tLcollinearCorr_ = exp(D_C_UL+F_tL); tRcollinearCorr_ = exp(D_C_UR+F_tR); bLcollinearCorr_ = exp(D_C_DL+F_bL); nuLcollinearCorr_ = exp(D_C_nuL); ELcollinearCorr_ = exp(D_C_EL); ERcollinearCorr_ = exp(D_C_ER); WtoWcollinearCorr_ = exp(D_C_WW); WtoZcollinearCorr_ = exp(D_C_WZ); WtoAcollinearCorr_ = exp(D_C_WA); BtoZcollinearCorr_ = exp(D_C_BZ); BtoAcollinearCorr_ = exp(D_C_BA); PhitoWcollinearCorr_ = exp(D_C_PhiW); PhitoZcollinearCorr_ = exp(D_C_PhiZ); PhitoHcollinearCorr_ = exp(D_C_PhiH); GcollinearCorr_ = exp(D_C_GG); } WaveFunctionCorrections CollinearSudakov::waveFunctionCorrections(Energy EWScale) { static const Complex I(0.,1.); using Constants::pi; double lZ = 2.0*log(ElectroWeakReweighter::coupling()->mZ()/EWScale); WaveFunctionCorrections WFC; // From Manohar, 2/12/12: (these assume mH=125, mZ=91.1876, mW=80.399) WFC.RW = (0.8410283377963967 - 9.424777961271568*I) + 1.2785863646210789*lZ; WFC.RA = (1.4835982362022198 + 1.855775680704845*pow(10.,-9)*I) - 0.27209907467584415*lZ; WFC.RZ = (1.5114724841549798 - 9.926944419863688*I) + 1.0834802397165764*lZ; WFC.RAtoZ = (0.3667485032811715 - 2.2770907130064835*I) - 1.2994544609942593*lZ; WFC.RZtoA = -0.2095310079712942 + 0.8320191107808546*lZ; WFC.RH = (12.229832449946716 - 1.7643103462419842*10.0*pow(10.,-12)*I) + 5.309998583664737*lZ; WFC.RPhi = (5.569012418081201 + 1.5439133581417356*0.10*pow(10.,-9)*I) + 5.309998583664737*lZ; WFC.RPhi3 = (8.945333042265943 + 5.499309445612249*pow(10.,-12)*I) + 5.309998583664737*lZ; WFC.EW = (3.967645734304811 + 4.712388980384717*I) + 2.238332625165702*lZ; WFC.EZ = (5.916079892937651 + 4.96347220970469*I) + 2.1132591719740788*lZ; WFC.RW *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.RA *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.RZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.RAtoZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.RZtoA *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.RPhi *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.EW *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.EZ *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.RPhi3 *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.RH *= ElectroWeakReweighter::coupling()->a2(EWScale)/(4.0*pi); WFC.RW = WFC.RW.real(); WFC.RA = WFC.RA.real(); WFC.RZ = WFC.RZ.real(); WFC.RAtoZ = WFC.RAtoZ.real(); WFC.RZtoA = WFC.RZtoA.real(); WFC.RPhi = WFC.RPhi.real(); WFC.RPhi3 = WFC.RPhi3.real(); WFC.RH = WFC.RH.real(); // Also from Manohar, 2/10/12 WFC.fFW0 = -3.7946842553189453 - 4.709019671388589*I; WFC.fF0W = 3.8181790485161176; WFC.fFZZ = 1.364503250377989 + 0.*I; WFC.aHH = -0.474396977740686 + 0.*I; WFC.aZZ = -0.6806563210877914 + 0.*I; WFC.aW0 = 0.49036102506811907 + 1.9323351450971642*I; WFC.a0W = -1.2184776671065072; WFC.bHH = 1.234775745474991 + 0.*I; WFC.bZZ = 1.7303526426747613 + 0.*I; WFC.cHH = 0.33140608274387473 + 0.*I; WFC.cZZ = 0.4961363208382017 + 0.*I; WFC.cW0 = -1.005299829777395 + 1.063048500002757*I; WFC.atHH = -0.237198488870343 + 0.*I; WFC.atZZ = -0.3403281605438957 + 0.*I; WFC.atW0 = 0.24518051253405954 + 0.9661675725485821*I; WFC.at0W = -0.6092388335532536; WFC.ctHH = 0.16570304137193737 + 0.*I; WFC.ctZZ = 0.24806816041910085 + 0.*I; WFC.ctW0 = -0.5026499148886975 + 0.5315242500013785*I; WFC.btHH = -0.30869393636874776 + 0.*I; WFC.btZZ = -0.4325881606686903 + 0.*I; WFC.fs10 = -2.289868133696459; WFC.fs1ZW = 1.8627978596240622; WFC.fsWZWZ = 1.1866922529667008; WFC.fsZW1 = 1.0840307611156266; WFC.fs01 = 2.2898681336964595; WFC.fsHW1 = -0.32306745562682404; WFC.fsHZ1 = 0.4042992116255279; WFC.fs1HW = 3.330954543719127; WFC.fs1HZ = 2.69768201932412; return WFC; } Complex CollinearSudakov::CollinearDw(Energy2 s, Energy EWScale) { using Constants::pi; using GroupInvariants::PlusLog; double lw = log(ElectroWeakReweighter::coupling()->mW()/EWScale); //s = s/2.; // This should not be here... I think this is a discrepency with Sascha return ElectroWeakReweighter::coupling()->aW(EWScale)/(4.0*pi)* (4.0*lw*0.5*PlusLog(s/EWScale/EWScale) - 2.0*lw*lw - 3.0*lw - 5.0/12.0*pi*pi + 9.0/4.0); } Complex CollinearSudakov::CollinearDz(Energy2 s, Energy EWScale) { using Constants::pi; using GroupInvariants::PlusLog; double lz = log(ElectroWeakReweighter::coupling()->mZ()/EWScale); return ElectroWeakReweighter::coupling()->aZ(EWScale)/(4.0*pi)* (4.0*lz*0.5*PlusLog(s/EWScale/EWScale) - 2.0*lz*lz - 3.0*lz - 5.0/12.0*pi*pi + 9.0/4.0); } namespace { void writeLine(ofstream & file, vector x, vector y, string name,string title, string colour, string style) { file << "# BEGIN HISTO1D "+name +"\n"; file << "SmoothLine=1\n"; file << "ErrorBars=0\n"; file << "Path="+name+"\n"; file << "Title="+title+"\n"; file << "LineColor="+colour+"\n"; file << "LineStyle="+style +"\n"; file << "# xlow xhigh val errminus errplus\n"; for(unsigned int ix=0;ix np; vector uL,uR,tL,tR,dL,dR,bL,bR,nuL,eL,eR,Wgamma,Bgamma,g,WT,WL,WZT,BZT,ZL,H; Energy mZ = ElectroWeakReweighter::coupling()->mZ(); for(Energy x=200.*GeV;x<5000.*GeV;x*=1.02) { Energy2 s(sqr(x)); np.push_back(x); evaluateHighScale(x,mZ,s); evaluateMatching(mZ,s); uL .push_back(real(highColQ_ * ULcollinearCorr_ )); uR .push_back(real(highColU_ * URcollinearCorr_ )); tL .push_back(real(highColQt_ * tLcollinearCorr_ )); tR .push_back(real(highColtR_ * tRcollinearCorr_ )); dL .push_back(real(highColQ_ * DLcollinearCorr_ )); dR .push_back(real(highColD_ * DRcollinearCorr_ )); bL .push_back(real(highColQt_ * bLcollinearCorr_ )); bR .push_back(real(highColD_ * DRcollinearCorr_ )); nuL .push_back(real(highColL_ * nuLcollinearCorr_ )); eL .push_back(real(highColL_ * ELcollinearCorr_ )); eR .push_back(real(highColE_ * ERcollinearCorr_ )); Wgamma.push_back(real(highColW_ * WtoAcollinearCorr_ )); Bgamma.push_back(real(highColB_ * BtoAcollinearCorr_ )); g .push_back(real(highColG_ * GcollinearCorr_ )); WT .push_back(real(highColW_ * WtoWcollinearCorr_ )); WL .push_back(real(highColPhi_ * PhitoWcollinearCorr_)); WZT .push_back(real(highColW_ * WtoZcollinearCorr_ )); BZT .push_back(real(highColB_ * BtoZcollinearCorr_ )); ZL .push_back(real(highColPhi_ * PhitoZcollinearCorr_)); H .push_back(real(highColPhi_ * PhitoHcollinearCorr_)); } ofstream fig1a("fig1a.dat"); fig1a << "#BEGIN PLOT\n"; fig1a << "LegendOnly=/uL /uR /tL /tR\n"; fig1a << "DrawOnly=/uL /uR /tL /tR\n"; fig1a << "Title=Figure 1a\n"; fig1a << "RatioPlot=0\n"; fig1a << "LogX=1\n"; fig1a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig1a << "YLabel=u, t\n"; fig1a << "Legend=1\n"; fig1a << "XMin=250.\n"; fig1a << "YMin=0.7\n"; fig1a << "YMax=1.05\n"; fig1a << "# END PLOT\n"; writeLine(fig1a,np,uL,"/uL","$u_L$","green","dotted" ); writeLine(fig1a,np,uR,"/uR","$u_R$","cyan" ,"solid" ); writeLine(fig1a,np,tL,"/tL","$t_L$","red" ,"dashed" ); writeLine(fig1a,np,tR,"/tR","$t_R$","blue" ,"dotdashed"); fig1a.close(); ofstream fig1b("fig1b.dat"); fig1b << "#BEGIN PLOT\n"; fig1b << "LegendOnly=/dL /dR /bL /bR\n"; fig1b << "DrawOnly=/dL /dR /bL /bR\n"; fig1b << "Title=Figure 1b\n"; fig1b << "RatioPlot=0\n"; fig1b << "LogX=1\n"; fig1b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig1b << "YLabel=d, b\n"; fig1b << "Legend=1\n"; fig1b << "XMin=250.\n"; fig1b << "YMin=0.7\n"; fig1b << "YMax=1.05\n"; fig1b << "# END PLOT\n"; writeLine(fig1b,np,dL,"/dL","$d_L$","green","dotted" ); writeLine(fig1b,np,dR,"/dR","$d_R$","cyan" ,"solid" ); writeLine(fig1b,np,bL,"/bL","$b_L$","red" ,"dashed" ); writeLine(fig1b,np,bR,"/bR","$b_R$","blue" ,"dotdashed"); fig1b.close(); ofstream fig2("fig2.dat"); fig2 << "#BEGIN PLOT\n"; fig2 << "LegendOnly=/uL /uR /dL /dR\n"; fig2 << "DrawOnly=/uL /uR /dL /dR\n"; fig2 << "Title=Figure 2\n"; fig2 << "RatioPlot=0\n"; fig2 << "LogX=1\n"; fig2 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig2 << "YLabel=u, d\n"; fig2 << "Legend=1\n"; fig2 << "XMin=250.\n"; fig2 << "YMin=0.7\n"; fig2 << "YMax=1.05\n"; fig2 << "# END PLOT\n"; writeLine(fig2,np,uL,"/uL","$u_L$","green","dotted" ); writeLine(fig2,np,uR,"/uR","$u_R$","cyan" ,"solid" ); writeLine(fig2,np,dL,"/dL","$d_L$","red" ,"dashed" ); writeLine(fig2,np,dR,"/dR","$d_R$","blue" ,"dotdashed"); fig2.close(); ofstream fig3("fig3.dat"); fig3 << "#BEGIN PLOT\n"; fig3 << "LegendOnly=/nuL /eL /eR\n"; fig3 << "DrawOnly=/nuL /eL /eR\n"; fig3 << "Title=Figure 3\n"; fig3 << "RatioPlot=0\n"; fig3 << "LogX=1\n"; fig3 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig3 << "YLabel=$\\nu$, $e$\n"; fig3 << "Legend=1\n"; fig3 << "XMin=250.\n"; fig3 << "YMin=0.9\n"; fig3 << "YMax=1.05\n"; fig3 << "# END PLOT\n"; writeLine(fig3,np,nuL,"/nuL","$\\nu_L$","blue","dashed"); writeLine(fig3,np, eL,"/eL" ,"$e_L$" ,"red" ,"dotted"); writeLine(fig3,np, eR,"/eR" ,"$e_R$" ,"red" ,"solid" ); fig3.close(); ofstream fig5("fig5.dat"); fig5 << "#BEGIN PLOT\n"; fig5 << "LegendOnly=/g /Wgamma /Bgamma\n"; fig5 << "DrawOnly=/g /Wgamma /Bgamma\n"; fig5 << "Title=Figure 5\n"; fig5 << "RatioPlot=0\n"; fig5 << "LogX=1\n"; fig5 << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig5 << "YLabel=$\\gamma$, g\n"; fig5 << "Legend=1\n"; fig5 << "XMin=250.\n"; fig5 << "YMin=0.7\n"; fig5 << "YMax=1.05\n"; fig5 << "# END PLOT\n"; writeLine(fig5,np,g ,"/g" ,"$g$" ,"blue","dashed"); writeLine(fig5,np,Wgamma,"/Wgamma","$W\\to\\gamma$","red" ,"solid" ); writeLine(fig5,np,Bgamma,"/Bgamma","$B\\to\\gamma$","red" ,"dotted"); fig5.close(); ofstream fig6a("fig6a.dat"); fig6a << "#BEGIN PLOT\n"; fig6a << "LegendOnly=/WT /WL\n"; fig6a << "DrawOnly=/WT /WL\n"; fig6a << "Title=Figure 6a\n"; fig6a << "RatioPlot=0\n"; fig6a << "LogX=1\n"; fig6a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig6a << "YLabel=$Z_L$, $Z_T$, $H$\n"; fig6a << "Legend=1\n"; fig6a << "XMin=250.\n"; fig6a << "YMin=0.7\n"; fig6a << "YMax=1.05\n"; fig6a << "# END PLOT\n"; writeLine(fig6a,np,WT,"/WT","$W_T$","red" ,"solid"); writeLine(fig6a,np,WL,"/WL","$W_L$","blue" ,"dashed" ); fig6a.close(); ofstream fig6b("fig6b.dat"); fig6b << "#BEGIN PLOT\n"; fig6b << "LegendOnly=/WZT /BZT /ZL /H\n"; fig6b << "DrawOnly=/WZT /BZT /ZL /H\n"; fig6b << "Title=Figure 6b\n"; fig6b << "RatioPlot=0\n"; fig6b << "LogX=1\n"; fig6b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig6b << "YLabel=d, b\n"; fig6b << "Legend=1\n"; fig6b << "XMin=250.\n"; fig6b << "YMin=0.7\n"; fig6b << "YMax=1.05\n"; fig6b << "# END PLOT\n"; writeLine(fig6b,np,WZT,"/WZT","$W\\to Z_T$","red" ,"solid" ); writeLine(fig6b,np,BZT,"/BZT","$B\\to Z_T$","red" ,"dotted" ); writeLine(fig6b,np,ZL ,"/ZL ","$Z_L$" ,"blue" ,"dashed" ); writeLine(fig6b,np,H ,"/H ","$H$" ,"green","dotdashed"); fig6b.close(); np.clear(); vector e30,e50,q30,q50,g30,g50; for(Energy x=200.*GeV;x<5000.*GeV;x*=1.02) { Energy2 s(sqr(x)); np.push_back(x); evaluateLowScale(mZ,30.*GeV,s); e30.push_back(real(lowColE_)); q30.push_back(real(lowColU_)); g30.push_back(real(lowColG_)); evaluateLowScale(mZ,50.*GeV,s); e50.push_back(real(lowColE_)); q50.push_back(real(lowColU_)); g50.push_back(real(lowColG_)); } ofstream fig4a("fig4a.dat"); fig4a << "#BEGIN PLOT\n"; fig4a << "LegendOnly=/e30 /e50\n"; fig4a << "DrawOnly=/e30 /e50\n"; fig4a << "Title=Figure 4a\n"; fig4a << "RatioPlot=0\n"; fig4a << "LogX=1\n"; fig4a << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig4a << "YLabel=e\n"; fig4a << "Legend=1\n"; fig4a << "XMin=250.\n"; fig4a << "YMin=0.7\n"; fig4a << "YMax=1.05\n"; fig4a << "# END PLOT\n"; writeLine(fig4a,np,e30,"/e30","e $(\\mu_f=30\\,\\mathrm{GeV})$","red" ,"solid" ); writeLine(fig4a,np,e50,"/e50","e $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dashed"); fig4a.close(); ofstream fig4b("fig4b.dat"); fig4b << "#BEGIN PLOT\n"; fig4b << "LegendOnly=/q30 /q50 /g30 /g50\n"; fig4b << "DrawOnly=/q30 /q50 /g30 /g50\n"; fig4b << "Title=Figure 4a\n"; fig4b << "RatioPlot=0\n"; fig4b << "LogX=1\n"; fig4b << "XLabel=$\\bar{n}\\cdot p$ [GeV]\n"; fig4b << "YLabel=e\n"; fig4b << "Legend=1\n"; fig4b << "XMin=250.\n"; fig4b << "YMin=0.5\n"; fig4b << "YMax=1.05\n"; fig4b << "# END PLOT\n"; writeLine(fig4b,np,q30,"/q30","q $(\\mu_f=30\\,\\mathrm{GeV})$","red" ,"solid" ); writeLine(fig4b,np,q50,"/q50","q $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dashed"); writeLine(fig4b,np,g30,"/g30","g $(\\mu_f=30\\,\\mathrm{GeV})$","green" ,"dotted" ); writeLine(fig4b,np,g50,"/g50","g $(\\mu_f=50\\,\\mathrm{GeV})$","blue","dotdashed"); fig4b.close(); } boost::numeric::ublas::matrix CollinearSudakov::electroWeakMatching(Energy EWScale, Energy2 s, Herwig::EWProcess::Process process, bool oneLoop) { using namespace EWProcess; // calculate the matching coefficients evaluateMatching(EWScale,s); // fill the matrix boost::numeric::ublas::matrix result(1,1); switch (process) { case QQQQ: case QQQQiden: { unsigned int numGauge = 4, numBrokenGauge = 12; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result(0,0) = result(6,0) = ULcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_; result(3,0) = result(9,0) = DLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_; for (int i=0; i<12; i++) { if (i!=0 && i!=3 && i!=6 && i!=9) { result(i,0) = ULcollinearCorr_*ULcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_; } } DuplicateColumn0(result); } break; case QtQtQQ: { unsigned int numGauge = 4, numBrokenGauge = 12; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result(0,0) = result(6,0) = ULcollinearCorr_*ULcollinearCorr_*tLcollinearCorr_*tLcollinearCorr_; result(3,0) = result(9,0) = DLcollinearCorr_*DLcollinearCorr_*bLcollinearCorr_*bLcollinearCorr_; for (int i=0; i<12; i++) { if (i==4 || i==5 || i==10 || i==11) { result(i,0) = ULcollinearCorr_*tLcollinearCorr_*DLcollinearCorr_*bLcollinearCorr_; } else if (i==1 || i==7) { result(i,0) = DLcollinearCorr_*DLcollinearCorr_*tLcollinearCorr_*tLcollinearCorr_; } else if (i==2 || i==8) { result(i,0) = ULcollinearCorr_*ULcollinearCorr_*bLcollinearCorr_*bLcollinearCorr_; } } DuplicateColumn0(result); } break; case QQUU: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*URcollinearCorr_*URcollinearCorr_; result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*URcollinearCorr_*URcollinearCorr_; DuplicateColumn0(result); } break; case QtQtUU: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*URcollinearCorr_*URcollinearCorr_; result(1,0) = result(3,0) = bLcollinearCorr_*bLcollinearCorr_*URcollinearCorr_*URcollinearCorr_; DuplicateColumn0(result); } break; case QQtRtR: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*tRcollinearCorr_*tRcollinearCorr_; result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*tRcollinearCorr_*tRcollinearCorr_; DuplicateColumn0(result); } break; case QQDD: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; result(1,0) = result(3,0) = DLcollinearCorr_*DLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; DuplicateColumn0(result); } break; case QtQtDD: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; result(1,0) = result(3,0) = bLcollinearCorr_*bLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; DuplicateColumn0(result); } break; case QQLL: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 6; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_; result(1,0) = nuLcollinearCorr_*nuLcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_; result(2,0) = ELcollinearCorr_*ELcollinearCorr_*ULcollinearCorr_*ULcollinearCorr_; result(3,0) = ELcollinearCorr_*ELcollinearCorr_*DLcollinearCorr_*DLcollinearCorr_; result(4,0) = result(5,0) = nuLcollinearCorr_*ELcollinearCorr_*ULcollinearCorr_*DLcollinearCorr_; DuplicateColumn0(result); } break; case QQEE: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = ULcollinearCorr_*ULcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_; result(1,0) = DLcollinearCorr_*DLcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_; DuplicateColumn0(result); } break; case UUUU: case UUUUiden: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = URcollinearCorr_*URcollinearCorr_*URcollinearCorr_*URcollinearCorr_; DuplicateColumn0(result); } break; case tRtRUU: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = tRcollinearCorr_*tRcollinearCorr_*URcollinearCorr_*URcollinearCorr_; DuplicateColumn0(result); } break; case UUDD: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = URcollinearCorr_*URcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; DuplicateColumn0(result); } break; case tRtRDD: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = tRcollinearCorr_*tRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; DuplicateColumn0(result); } break; case UULL: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*URcollinearCorr_*URcollinearCorr_; result(1,0) = ELcollinearCorr_*ELcollinearCorr_*URcollinearCorr_*URcollinearCorr_; DuplicateColumn0(result); } break; case UUEE: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 1; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = URcollinearCorr_*URcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_; DuplicateColumn0(result); } break; case DDDD: case DDDDiden: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = DRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; DuplicateColumn0(result); } break; case DDLL: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; result(1,0) = ELcollinearCorr_*ELcollinearCorr_*DRcollinearCorr_*DRcollinearCorr_; DuplicateColumn0(result); } break; case DDEE: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 1; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = DRcollinearCorr_*DRcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_; DuplicateColumn0(result); } break; case LLLL: case LLLLiden: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 6; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_; result(1,0) = nuLcollinearCorr_*nuLcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_; result(2,0) = ELcollinearCorr_*ELcollinearCorr_*nuLcollinearCorr_*nuLcollinearCorr_; result(3,0) = ELcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_*ELcollinearCorr_; result(4,0) = result(5,0) = nuLcollinearCorr_*ELcollinearCorr_*nuLcollinearCorr_*ELcollinearCorr_; DuplicateColumn0(result); } break; case LLEE: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = nuLcollinearCorr_*nuLcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_; result(1,0) = ELcollinearCorr_*ELcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_; DuplicateColumn0(result); } break; case EEEE: case EEEEiden: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 1; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = ERcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_*ERcollinearCorr_; DuplicateColumn0(result); } break; case QQWW: case LLWW: { unsigned int numGauge = 5; unsigned int numBrokenGauge = 20; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; for (unsigned int row = 0; row < result.size1(); row++) { for (unsigned int col = 0; col < result.size2(); col++) { // Boson Collinear Corr_ections: if (col==0 || col==1) { if (row==0 || row==1 || row==6 || row==7) result(row,col) = (WtoWcollinearCorr_*WtoWcollinearCorr_); if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*WtoZcollinearCorr_); if (row==3 || row==4 || row==9 || row==10) result(row,col) = (WtoZcollinearCorr_*WtoAcollinearCorr_); if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*WtoAcollinearCorr_); if (row==12 || row==14) result(row,col) = (WtoWcollinearCorr_*WtoZcollinearCorr_); if (row==13 || row==15) result(row,col) = (WtoWcollinearCorr_*WtoAcollinearCorr_); if (row==16 || row==18) result(row,col) = (WtoWcollinearCorr_*WtoZcollinearCorr_); if (row==17 || row==19) result(row,col) = (WtoWcollinearCorr_*WtoAcollinearCorr_); } if (col==2) { if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*BtoZcollinearCorr_); if (row==3 || row==9) result(row,col) = (WtoZcollinearCorr_*BtoAcollinearCorr_); if (row==4 || row==10) result(row,col) = (WtoAcollinearCorr_*BtoZcollinearCorr_); if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*BtoAcollinearCorr_); if (row==14) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_); if (row==15) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_); if (row==16) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_); if (row==17) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_); } if (col==3) { if (row==2 || row==8) result(row,col) = (WtoZcollinearCorr_*BtoZcollinearCorr_); if (row==3 || row==9) result(row,col) = (WtoAcollinearCorr_*BtoZcollinearCorr_); if (row==4 || row==10) result(row,col) = (WtoZcollinearCorr_*BtoAcollinearCorr_); if (row==5 || row==11) result(row,col) = (WtoAcollinearCorr_*BtoAcollinearCorr_); if (row==12) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_); if (row==13) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_); if (row==18) result(row,col) = (WtoWcollinearCorr_*BtoZcollinearCorr_); if (row==19) result(row,col) = (WtoWcollinearCorr_*BtoAcollinearCorr_); } if (col==4) { if (row==2 || row==8) result(row,col) = (BtoZcollinearCorr_*BtoZcollinearCorr_); if (row==3 || row==4 || row==9 || row==10) result(row,col) = (BtoZcollinearCorr_*BtoAcollinearCorr_); if (row==5 || row==11) result(row,col) = (BtoAcollinearCorr_*BtoAcollinearCorr_); } // Particle Collinear Corr_ections: if (process==QQWW) { if (row<6) result(row,col) *= (ULcollinearCorr_*ULcollinearCorr_); if ((row>=6)&&(row<12)) result(row,col) *= (DLcollinearCorr_*DLcollinearCorr_); if (row>=12) result(row,col) *= (ULcollinearCorr_*DLcollinearCorr_); } else if (process==LLWW) { if (row<6) result(row,col) *= (nuLcollinearCorr_*nuLcollinearCorr_); if ((row>=6)&&(row<12)) result(row,col) *= (ELcollinearCorr_*ELcollinearCorr_); if (row>=12) result(row,col) *= (nuLcollinearCorr_*ELcollinearCorr_); } } } } break; case QQPhiPhi: case LLPhiPhi: { unsigned int numGauge = 2; unsigned int numBrokenGauge = 14; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; for (unsigned int row = 0; row < result.size1(); row++) { // Boson Colinear Corr_ections: if (row==0 || row==5) result(row,0) = (PhitoWcollinearCorr_*PhitoWcollinearCorr_); if (row==1 || row==6) result(row,0) = (PhitoZcollinearCorr_*PhitoZcollinearCorr_); if (row==2 || row==3 || row==7 || row==8) result(row,0) = (PhitoZcollinearCorr_*PhitoHcollinearCorr_); if (row==4 || row==9) result(row,0) = (PhitoHcollinearCorr_*PhitoHcollinearCorr_); if (row==10) result(row,0) = (PhitoWcollinearCorr_*PhitoZcollinearCorr_); if (row==11) result(row,0) = (PhitoWcollinearCorr_*PhitoHcollinearCorr_); if (row==12) result(row,0) = (PhitoWcollinearCorr_*PhitoZcollinearCorr_); if (row==13) result(row,0) = (PhitoWcollinearCorr_*PhitoHcollinearCorr_); // Particle Colinear Corr_ections: if (process==QQPhiPhi) { if (row<5) result(row,0) *= (ULcollinearCorr_*ULcollinearCorr_); if ((row>=5)&&(row<10)) result(row,0) *= (DLcollinearCorr_*DLcollinearCorr_); if (row>=10) result(row,0) *= (ULcollinearCorr_*DLcollinearCorr_); } else if (process==LLPhiPhi) { if (row<5) result(row,0) *= (nuLcollinearCorr_*nuLcollinearCorr_); if ((row>=5)&&(row<10)) result(row,0) *= (ELcollinearCorr_*ELcollinearCorr_); if (row>=10) result(row,0) *= (nuLcollinearCorr_*ELcollinearCorr_); } } DuplicateColumn0(result); } break; case QQWG: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 6; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = ULcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoWcollinearCorr_; result(2,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*WtoZcollinearCorr_; result(3,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*WtoAcollinearCorr_; result(4,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoZcollinearCorr_; result(5,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*WtoAcollinearCorr_; DuplicateColumn0(result); } break; case QQBG: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_; result(1,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_; result(2,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_; result(3,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_; DuplicateColumn0(result); } break; case QQGG: { unsigned int numGauge = 3; unsigned int numBrokenGauge = 6; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = result(2,0) = ULcollinearCorr_*ULcollinearCorr_*GcollinearCorr_*GcollinearCorr_; result(3,0) = result(4,0) = result(5,0) = DLcollinearCorr_*DLcollinearCorr_*GcollinearCorr_*GcollinearCorr_; DuplicateColumn0(result); } break; case QtQtGG: { unsigned int numGauge = 3; unsigned int numBrokenGauge = 6; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = result(2,0) = tLcollinearCorr_*tLcollinearCorr_*GcollinearCorr_*GcollinearCorr_; result(3,0) = result(4,0) = result(5,0) = bLcollinearCorr_*bLcollinearCorr_*GcollinearCorr_*GcollinearCorr_; DuplicateColumn0(result); } break; case UUBB: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = URcollinearCorr_*URcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_; result(1,0) = URcollinearCorr_*URcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_; result(2,0) = URcollinearCorr_*URcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_; result(3,0) = URcollinearCorr_*URcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_; DuplicateColumn0(result); } break; case UUPhiPhi: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 5; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = URcollinearCorr_*URcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_; result(1,0) = URcollinearCorr_*URcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_; result(2,0) = URcollinearCorr_*URcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_; result(3,0) = URcollinearCorr_*URcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_; result(4,0) = URcollinearCorr_*URcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_; DuplicateColumn0(result); } break; case UUBG: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_; result(1,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_; DuplicateColumn0(result); } break; case UUGG: { unsigned int numGauge = 3; unsigned int numBrokenGauge = 3; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = result(2,0) = URcollinearCorr_*URcollinearCorr_*GcollinearCorr_*GcollinearCorr_; DuplicateColumn0(result); } break; case tRtRGG: { unsigned int numGauge = 3; unsigned int numBrokenGauge = 3; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = result(2,0) = tRcollinearCorr_*tRcollinearCorr_*GcollinearCorr_*GcollinearCorr_; DuplicateColumn0(result); } break; case DDBB: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = DRcollinearCorr_*DRcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_; result(1,0) = DRcollinearCorr_*DRcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_; result(2,0) = DRcollinearCorr_*DRcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_; result(3,0) = DRcollinearCorr_*DRcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_; DuplicateColumn0(result); } break; case DDPhiPhi: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 5; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_; result(1,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_; result(2,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_; result(3,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_; result(4,0) = DRcollinearCorr_*DRcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_; DuplicateColumn0(result); } break; case DDBG: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 2; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*BtoZcollinearCorr_; result(1,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*BtoAcollinearCorr_; DuplicateColumn0(result); } break; case DDGG: { unsigned int numGauge = 3; unsigned int numBrokenGauge = 3; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = result(1,0) = result(2,0) = DRcollinearCorr_*DRcollinearCorr_*GcollinearCorr_*GcollinearCorr_; DuplicateColumn0(result); } break; case EEBB: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 4; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = ERcollinearCorr_*ERcollinearCorr_*BtoZcollinearCorr_*BtoZcollinearCorr_; result(1,0) = ERcollinearCorr_*ERcollinearCorr_*BtoZcollinearCorr_*BtoAcollinearCorr_; result(2,0) = ERcollinearCorr_*ERcollinearCorr_*BtoAcollinearCorr_*BtoZcollinearCorr_; result(3,0) = ERcollinearCorr_*ERcollinearCorr_*BtoAcollinearCorr_*BtoAcollinearCorr_; DuplicateColumn0(result); } break; case EEPhiPhi: { unsigned int numGauge = 1; unsigned int numBrokenGauge = 5; result = boost::numeric::ublas::zero_matrix(numBrokenGauge,numGauge); result *= 0.0; result(0,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoWcollinearCorr_*PhitoWcollinearCorr_; result(1,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoZcollinearCorr_*PhitoZcollinearCorr_; result(2,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoHcollinearCorr_*PhitoZcollinearCorr_; result(3,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoZcollinearCorr_*PhitoHcollinearCorr_; result(4,0) = ERcollinearCorr_*ERcollinearCorr_*PhitoHcollinearCorr_*PhitoHcollinearCorr_; DuplicateColumn0(result); } break; default: assert(false); } // This is done at the end instead of the beginning for result.size1() and cols() if (!oneLoop) { boost::numeric::ublas::matrix OnesMatrix(result.size1(),result.size2()); for (unsigned int i=0; i +CollinearSudakov::highEnergyRunning(Energy highScale, Energy EWScale, Energy2 s, + Herwig::EWProcess::Process process, + bool fixedOrder) { + using namespace EWProcess; + // perform the calculation + evaluateHighScale(highScale,EWScale,s); + Complex colW(highColW_); + Complex colB(highColB_); + Complex colG(highColG_); + Complex colQ(highColQ_); + Complex colQt(highColQt_); + Complex colU(highColU_); + Complex coltR(highColtR_); + Complex colD(highColD_); + Complex colL(highColL_); + Complex colE(highColE_); + Complex colPhi(highColPhi_); + if (fixedOrder) { + /* colX not necessarily positive for s = (1000TeV)^2 for the following: + colW = log(colW.real()); + colB = log(colB.real()); + colPhi = log(colPhi.real()); + */ + colG = log(colG.real()); + colQ = log(colQ.real()); + colQt = log(colQt.real()); + colU = log(colU.real()); + coltR = log(coltR.real()); + colD = log(colD.real()); + colL = log(colL.real()); + colE = log(colE.real()); + } + // set up the matrix + boost::numeric::ublas::matrix result; + unsigned int numGauge(0); + switch (process) { + + case QQQQ: + case QQQQiden: + case QtQtQQ: + numGauge = 4; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (process!=QtQtQQ) { + for (unsigned int i=0; i(numGauge,numGauge); + if (process==QQUU) { + for (unsigned int i=0; i(numGauge,numGauge); + if (process==QQDD) { + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + if (process!=tRtRUU) { + for (unsigned int i=0; i(numGauge,numGauge); + if (process==UUDD) { + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + for (unsigned int i=0; i(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = result(1,1) = 1.0+colQ+colQ+colW+colW; + result(2,2) = result(3,3) = 1.0+colQ+colQ+colW+colB; + result(4,4) = 1.0+colQ+colQ+colB+colB; + } + else { + result(0,0) = result(1,1) = colQ*colQ*colW*colW; + result(2,2) = result(3,3) = colQ*colQ*colW*colB; + result(4,4) = colQ*colQ*colB*colB; + } + break; + + case QQPhiPhi: + numGauge = 2; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = result(1,1) = 1.0+colQ+colQ+colPhi+colPhi; + } + else { + result(0,0) = result(1,1) = colQ*colQ*colPhi*colPhi; + } + break; + + case QQWG: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colQ+colQ+colW+colG; + } + else { + result(0,0) = colQ*colQ*colW*colG; + } + break; + + case QQBG: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colQ+colQ+colB+colG; + } + else { + result(0,0) = colQ*colQ*colB*colG; + } + break; + + case QQGG: + case QtQtGG: + numGauge = 3; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (process==QQGG) { + if (fixedOrder) { + result(0,0) = result(1,1) = result(2,2) = 1.0+colQ+colQ+colG+colG; + } + else { + result(0,0) = result(1,1) = result(2,2) = colQ*colQ*colG*colG; + } + } + else { + if (fixedOrder) { + result(0,0) = result(1,1) = result(2,2) = 1.0+colQt+colQt+colG+colG; + } + else { + result(0,0) = result(1,1) = result(2,2) = colQt*colQt*colG*colG; + } + } + break; + + case LLWW: + numGauge = 5; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = result(1,1) = 1.0+colL+colL+colW+colW; + result(2,2) = result(3,3) = 1.0+colL+colL+colW+colB; + result(4,4) = 1.0+colL+colL+colB+colB; + } + else { + result(0,0) = result(1,1) = colL*colL*colW*colW; + result(2,2) = result(3,3) = colL*colL*colW*colB; + result(4,4) = colL*colL*colB*colB; + } + break; + + case LLPhiPhi: + numGauge = 2; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = result(1,1) = 1.0+colL+colL+colPhi+colPhi; + } + else { + result(0,0) = result(1,1) = colL*colL*colPhi*colPhi; + } + break; + + case UUBB: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colU+colU+colB+colB; + } + else { + result(0,0) = colU*colU*colB*colB; + } + break; + + case UUPhiPhi: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colU+colU+colPhi+colPhi; + } + else { + result(0,0) = colU*colU*colPhi*colPhi; + } + break; + + case UUBG: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colU+colU+colB+colG; + } + else { + result(0,0) = colU*colU*colB*colG; + } + break; + + case UUGG: + case tRtRGG: + numGauge = 3; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (process==UUGG) { + if (fixedOrder) { + result(0,0) = result(1,1) = result(2,2) = 1.0+colU+colU+colG+colG; + } + else { + result(0,0) = result(1,1) = result(2,2) = colU*colU*colG*colG; + } + } + else { + if (fixedOrder) { + result(0,0) = result(1,1) = result(2,2) = 1.0+coltR+coltR+colG+colG; + } + else { + result(0,0) = result(1,1) = result(2,2) = coltR*coltR*colG*colG; + } + } + break; + + case DDBB: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colD+colD+colB+colB; + } + else { + result(0,0) = colD*colD*colB*colB; + } + break; + + case DDPhiPhi: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colD+colD+colPhi+colPhi; + } + else { + result(0,0) = colD*colD*colPhi*colPhi; + } + break; + + case DDBG: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colD+colD+colB+colG; + } + else { + result(0,0) = colD*colD*colB*colG; + } + break; + + case DDGG: + numGauge = 3; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = result(1,1) = result(2,2) = 1.0+colD+colD+colG+colG; + } + else { + result(0,0) = result(1,1) = result(2,2) = colD*colD*colG*colG; + } + break; + + case EEBB: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colE+colE+colB+colB; + } + else { + result(0,0) = colE*colE*colB*colB; + } + break; + + case EEPhiPhi: + numGauge = 1; + result = boost::numeric::ublas::zero_matrix(numGauge,numGauge); + if (fixedOrder) { + result(0,0) = 1.0+colE+colE+colPhi+colPhi; + } + else { + result(0,0) = colE*colE*colPhi*colPhi; + } + break; + + default: + assert(false); + } + + return result; +} diff --git a/MatrixElement/EW/CollinearSudakov.h b/MatrixElement/EW/CollinearSudakov.h --- a/MatrixElement/EW/CollinearSudakov.h +++ b/MatrixElement/EW/CollinearSudakov.h @@ -1,527 +1,535 @@ // -*- C++ -*- #ifndef Herwig_CollinearSudakov_H #define Herwig_CollinearSudakov_H // // This is the declaration of the CollinearSudakov class. // #include "ThePEG/Interface/Interfaced.h" #include "Herwig/Utilities/GSLIntegrator.h" #include #include "EWProcess.h" #include "CollinearSudakov.fh" namespace Herwig { using namespace ThePEG; /** * Struct for the wavefunction corrections */ struct WaveFunctionCorrections { Complex RW; Complex RA; Complex RZ; Complex RAtoZ; Complex RZtoA; Complex RPhi; Complex EW; Complex EZ; Complex RPhi3; Complex RH; Complex tLuLDiff; Complex bLdLDiff; Complex tRuRDiff; // The following are constants from parameter integrals: Complex fFW0; Complex fF0W; Complex fFZZ; Complex aHH; Complex aZZ; Complex aW0; Complex a0W; Complex bHH; Complex bZZ; Complex cHH; Complex cZZ; Complex cW0; Complex atHH; Complex atZZ; Complex atW0; Complex at0W; Complex ctHH; Complex ctZZ; Complex ctW0; Complex btHH; Complex btZZ; Complex fs10; Complex fs1ZW; Complex fsWZWZ; Complex fsZW1; Complex fs01; Complex fsHW1; Complex fsHZ1; Complex fs1HW; Complex fs1HZ; }; /** * Here is the documentation of the CollinearSudakov class. * * @see \ref CollinearSudakovInterfaces "The interfaces" * defined for CollinearSudakov. */ class CollinearSudakov: public Interfaced { public: /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ CollinearSudakov(); /** * The destructor. */ virtual ~CollinearSudakov(); //@} public: /** - * Evalaute the electroweka matching as a matrix + * Evalaute the electroweak matching as a matrix */ boost::numeric::ublas::matrix electroWeakMatching(Energy EWScale, Energy2 s, Herwig::EWProcess::Process process, bool oneLoop); + /** + * Evalaute the high energy running as a matrix + */ + boost::numeric::ublas::matrix + highEnergyRunning(Energy highScale, Energy EWScale, Energy2 s, + Herwig::EWProcess::Process process, + bool fixedOrder); + public: /** * Make plots for tests */ void makePlots(); protected: /** * Evaluate the high scale contributions */ void evaluateHighScale(Energy highScale, Energy EWScale, Energy2 S); /** * Evaluate the low scale contributions */ void evaluateLowScale(Energy EWScale, Energy lowScale, Energy2 S); /** * Evaluate the matching */ void evaluateMatching(Energy EWScale,Energy2 S); public: /** * The operator to be integrated */ InvEnergy operator ()(Energy mu) const { if(high_) return highScaleIntegrand(mu); else return lowScaleIntegrand(mu); } /** Argument type for GaussianIntegrator */ typedef Energy ArgType; /** Return type for GaussianIntegrator */ typedef InvEnergy ValType; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** * The integral of the high scale part of the Sudakov */ Complex highScaleIntegral(bool SU3, bool SU2, double Y, Energy2 s, Energy mu_h, Energy mu_l, bool fermion, bool longitudinal, double yukFactor); /** * the integral of the low scale part of the Sudakov */ Complex lowScaleIntegral(bool SU3, double Q, Energy2 s, Energy mu_h, Energy mu_l, bool fermion, double boostFactor); protected: /** * High-scale integrand */ InvEnergy highScaleIntegrand(Energy mu) const; /** * Low-scale integrand */ InvEnergy lowScaleIntegrand(Energy mu) const; /** * Calculate the wavefunction corrections */ WaveFunctionCorrections waveFunctionCorrections(Energy EWScale); /** * Collinear matiching for W */ Complex CollinearDw(Energy2 s, Energy EWScale); /** * Collinear matching for Z */ Complex CollinearDz(Energy2 s, Energy EWScale); 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; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ CollinearSudakov & operator=(const CollinearSudakov &); private: /** * Parameters for the integrand */ //@{ /** * Whether high or low scale */ bool high_; /** * Whether real of imaginary part */ bool real_; /** * \f$SU(3)\f$ */ bool SU3_; /** * */ bool SU2_; /** * */ double Y_; /** * */ Energy2 s_; /** * */ bool fermion_; /** * */ bool longitudinal_; /** * */ double yukFactor_; /** * */ double boostFactor_; /** * */ double Q_; //@} /** * Parameters */ //@{ /** * Order for the K terms */ int K_ORDER_; /** * Order for the B terms */ int B_ORDER_; //@} /** * Integrator */ GSLIntegrator integrator_; private: /** * Storage of the high scale pieces */ //@{ /** * */ Complex highColW_; /** * */ Complex highColB_; /** * */ Complex highColG_; /** * */ Complex highColQ_; /** * */ Complex highColQt_; /** * */ Complex highColU_; /** * */ Complex highColtR_; /** * */ Complex highColD_; /** * */ Complex highColL_; /** * */ Complex highColE_; /** * */ Complex highColPhi_; //@} /** * Storage of the low scale pieces */ //@{ /** * */ complex lowColW_; /** * */ Complex lowColA_; /** * */ Complex lowColG_; /** * */ Complex lowColU_; /** * */ Complex lowColt_; /** * */ Complex lowColD_; /** * */ Complex lowColE_; //@} /** * Storage of the matching parameters */ //@{ /** * */ Complex ULcollinearCorr_; /** * */ Complex DLcollinearCorr_; /** * */ Complex URcollinearCorr_; /** * */ Complex DRcollinearCorr_; /** * */ Complex tLcollinearCorr_; /** * */ Complex tRcollinearCorr_; /** * */ Complex bLcollinearCorr_; /** * */ Complex nuLcollinearCorr_; /** * */ Complex ELcollinearCorr_; /** * */ Complex ERcollinearCorr_; /** * */ Complex WtoWcollinearCorr_; /** * */ Complex WtoZcollinearCorr_; /** * */ Complex WtoAcollinearCorr_; /** * */ Complex BtoZcollinearCorr_; /** * */ Complex BtoAcollinearCorr_; /** * */ Complex PhitoWcollinearCorr_; /** * */ Complex PhitoZcollinearCorr_; /** * */ Complex PhitoHcollinearCorr_; /** * */ Complex GcollinearCorr_; //@} }; } #endif /* Herwig_CollinearSudakov_H */ diff --git a/MatrixElement/EW/GroupInvariants.h b/MatrixElement/EW/GroupInvariants.h --- a/MatrixElement/EW/GroupInvariants.h +++ b/MatrixElement/EW/GroupInvariants.h @@ -1,311 +1,328 @@ // -*- C++ -*- // // GroupInvariants.h is a part of Herwig - A multi-purpose Monte Carlo event generator // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // #ifndef HERWIG_GroupInvariants_H #define HERWIG_GroupInvariants_H #include "ThePEG/Config/ThePEG.h" #include "ThePEG/Config/Unitsystem.h" #include #include namespace Herwig { using namespace ThePEG; namespace GroupInvariants { /** * Simple struct for storing the different gauge contributions */ struct GaugeContributions { /** * Default Constructor */ GaugeContributions(double inSU3=0., double inSU2=0., double inU1=0.) : SU3(inSU3),SU2(inSU2),U1(inU1) {} /** * \f$SU(3)\f$ */ double SU3; /** * \f$SU(2)\f$ */ double SU2; /** * \f$U(1)\f$ */ double U1; }; /** * The \f$SU(N)\f$ \f$C_A\f$ */ inline double C_A(unsigned int N) { return N !=1 ? double(N) : 0.; } /** * The \f$SU(N)\f$ \f$C_F\f$ */ inline double C_F(unsigned int N) { return N !=1 ? 0.5*(double(N*N)-1.)/double(N) : 1.; } /* * The \f$SU(N)\f$ \f$C_d\f$ */ inline double C_d(unsigned int N) { return (double(N*N)-4.)/double(N); } /** * The \f$SU(N)\f$\f$C_1\f$ */ inline double C_1(unsigned int N) { double N2(N*N); return 0.25*(N2-1.0)/N2; } /** * \f$T_F\f$ */ inline double T_F(unsigned int N, bool high) { if(high) { return N !=1 ? 0.5 : 5.0/3.0; } else { return N !=1 ? 0.5 : 20.0/3.0; } } /** * \f$t_S\f$ */ inline double t_S(unsigned int, bool ) { return 0.5; } /** * / Number of complex scalars in the fundamental rep. of SU(N)/U(1) */ inline double n_S(unsigned int N, bool high) { if(high) { if(N==2 || N==1) return 1.0; else if(N==3) return 0.0; else assert(false); } else { if(N>=1&&N<=3) return 0.; else assert(false); } } /** * Number of Dirac Fermions in the fund. rep. of SU(N) (or U(1) for N==1) */ inline double n_F(unsigned int N, bool high) { if(high) { if(N==1) return 3.0; else if(N==2 || N==3) return 6.0; else assert(false); } else { if(N==1) return 1.0; else if(N==2) return 0.0; else if(N==3) return 5.0; else assert(false); } } /** * Find K_i for gauge group N. high=false for running at mu0.0) return log(arg); else if (arg<0.0) return log(-arg)+I*Constants::pi; else assert(false); } inline Complex MinusLog(double arg) { static const Complex I(0,1.0); if (arg>0.0) return log(arg); else if (arg<0.0) return log(-arg)-I*Constants::pi; else assert(false); } inline Complex getT(Energy2 s, Energy2 t) { return MinusLog(-t/GeV2) - MinusLog(-s/GeV2); } inline Complex getU(Energy2 s, Energy2 u) { return MinusLog(-u/GeV2) - MinusLog(-s/GeV2); } inline boost::numeric::ublas::matrix Gamma2(Complex U, Complex T) { boost::numeric::ublas::matrix output(2,2); static const Complex I(0,1.0); using Constants::pi; output(0,0) = (-3.0/2.0)*I*pi + (T+U); output(1,1) = (-3.0/2.0)*I*pi; output(0,1) = 2.0*(T-U); output(1,0) = (3.0/8.0)*(T-U); return output; } inline boost::numeric::ublas::matrix Gamma2w(Complex U, Complex T) { boost::numeric::ublas::matrix output = boost::numeric::ublas::zero_matrix(5,5); static const Complex I(0,1.0); using Constants::pi; output(0,0) += -I*pi*11.0/4.0; output(0,1) += U-T; output(1,0) += 2.0*(U-T); output(1,1) += -I*pi*11.0/4.0 + (T+U); output(2,2) += -7.0/4.0*I*pi + (U+T); output(3,3) += -7.0/4.0*I*pi + (U+T); output(4,4) += -3.0/4.0*I*pi; return output; } inline boost::numeric::ublas::matrix Gamma2Singlet() { boost::numeric::ublas::matrix output = boost::numeric::ublas::zero_matrix(2,2); static const Complex I(0,1.0); using Constants::pi; output(0,0) = output(1,1) = -3.0/4.0*I*pi; return output; } + inline Complex Gamma1(double hypercharge) { + Complex I(0,1.0); + return -I*Constants::pi*(hypercharge*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 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; + } + /** * 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