Page MenuHomeHEPForge
Files F19251914

FKS.cc
No OneTemporary
Actions

Size
31 KB
Referenced Files
None
Subscribers
None

FKS.cc
View Options

/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2022
* \copyright GPLv2 or later
*/
#include "HEJ/FKS.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <iterator>
#include <limits>
#include <unordered_map>
#include <utility>
#include <gsl/gsl_sf_dilog.h> //Dilog for FKS implementation
#include "fastjet/PseudoJet.hh"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/LorentzVector.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "HEJ/event_types.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/utility.hh"
#include "HEJ/jets.hh"
#include "HEJ/Config.hh"
#include "HEJ/MatrixElement.hh"
namespace HEJ {
namespace FKS {
//////// FKS functions /////////
//Shorthand for paper references:
//FKS: arXiv:hep-ph/9512328
//ES: DOI 10.1016/0550-3213(86)90232-4
//MadFKS: 0908.4272
double ME2_NLO_FKS(Event const & ev, FKSConfig const & FKSParams) {
//Read all partons from event
auto const & incoming = ev.incoming();
auto const out_partons = filter_partons(ev.outgoing());
//HEJ labelling of incoming partons: pa has negative z component
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
//Calculate boost to partonic COM frame
const double Ea= pa.e();
const double Eb= pb.e();
double beta_to_COM = (1-(Eb/Ea))/(1+(Eb/Ea));
//Switch to FKS labelling of partons
CLHEP::HepLorentzVector kp=pb;
CLHEP::HepLorentzVector km=pa;
//Perform boost for incoming partons
kp.boostZ(beta_to_COM);
km.boostZ(beta_to_COM);
//For now, use central renormalisation scale
const double mur2 = pow(ev.central().mur,2);
//Temporarily require u(kp) d(km) incoming
if(!(incoming[0].type==1&&incoming[1].type==2)){
throw std::logic_error("FKS treatment currently requires u(kp) d(km) incoming.");
}
if(out_partons.size()==2){
//HEJ labelling of outgoing partons: y1<y2, f1=fa, f2=fb.
//Caution: .front() refers to position in list, ordered from lowest to highest rapidity.
const auto p1 = to_HepLorentzVector(out_partons.front());
const auto p2 = to_HepLorentzVector(out_partons.back());
//Temporarily require u(p2) d(p1) outgoing, i.e. FKL initial event.
if(!(out_partons.front().type==1&&out_partons.back().type==2)){
throw std::logic_error("FKS treatment currently requires u(kp) d(km) incoming.");
}
/* // Print event
std::cout<<"pa="<<pa<<std::endl;
std::cout<<"fa="<<incoming[0].type<<std::endl;
std::cout<<"pb="<<pb<<std::endl;
std::cout<<"fb="<<incoming[1].type<<std::endl;
std::cout<<"p1="<<p1<<std::endl;
std::cout<<"f1="<<out_partons.front().type<<std::endl;
std::cout<<"p2="<<p2<<std::endl;
std::cout<<"f2="<<out_partons.back().type<<std::endl;
*/
//Change to FKS labelling of partons
//Temporarily q=u=2, Q=d=1
CLHEP::HepLorentzVector kq=p2;
CLHEP::HepLorentzVector kQ=p1;
//Boost along z direction to partonic COM frame
kq.boostZ(beta_to_COM);
kQ.boostZ(beta_to_COM);
//Return C+S+V in (4.3) of MadFKS
//But factor out (as/2/p2)gs^4
//Q: does d\sigma^V include running coupling?
//return M_q_Q_to_q_Q(kp,km,kq,kQ);
return Regulated_2parton_q_Q_to_q_Q(mur2,FKSParams,kp,km,kq,kQ);
} else if(out_partons.size()==3){
return 0.;
} else {
throw std::logic_error("FKS treatment currently supports only 2 or 3 partons in final state.");
}
}
///// Exact squared amplitudes - required for validating FKS NLO implementation. /////
//Square of 0->qb Qb q Q amplitude. No flux or averaging factors. C.f. table 2 of ES.
double M2_qb_Qb_q_Q(
CLHEP::HepLorentzVector const & pqb,
CLHEP::HepLorentzVector const & pQb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pQ
){
const double sqQ= 2.*pq.dot(pQ);
const double sqqb= -2.*pq.dot(pqb);
const double sqQb= -2.*pq.dot(pQb);
const double kin = 2.*(sqQ*sqQ+sqQb*sqQb)/(sqqb*sqqb);
const double col = N_C*N_C-1;
return col*kin;
}
//Implements (3.4) of ES.
//Converted to all-outgoing convention, eps->0 limit taken,
//and reorganised kinematics into IR stable form.
double M2_qb_Qb_q_Q_g(
CLHEP::HepLorentzVector const & pqb,
CLHEP::HepLorentzVector const & pQb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pQ,
CLHEP::HepLorentzVector const & pg
){
//Overcomplete set of invariants, but allows compact numerator:
const double sqbQb= 2.*pqb.dot(pQb); //s
const double sqQ= 2.*pq.dot(pQ); //s'
const double sqbq= 2.*pqb.dot(pq); //t
const double sQbQ= 2.*pQb.dot(pQ); //t'
const double sqbQ= 2.*pqb.dot(pQ); //u
const double sQbq= 2.*pQb.dot(pq); //u'
const double kin0 = (sqbQ*sqbQ + sqbQb*sqbQb + sqQ*sqQ + sQbq*sQbq)/2.;
/*
//Old form of kinematics, numerically unstable in soft limit.
double kin1 = ((s14 + s23)*(-s14*s23 + s13*s24 + s12*s34) + s23*(s12*s24 + s13*s34) + s14*(s12*s13 + s24*s34))/4.;
double kin2 = (2.*s13*(s14 + s23)*s24 + 2.*s14*s23*(s13 + s24) + (s12 + s34)*(-s14*s23 - s13*s24 + s12*s34))/4.;
*/
//We need pairs of the following invariants to obtain IR-convergent numerator:
const double sqbg= 2.*pqb.dot(pg);
const double sQbg= 2.*pQb.dot(pg);
const double sqg= 2.*pq.dot(pg);
const double sQg= 2.*pQ.dot(pg);
const double kin1s = (1/4.)*(sqbQ*(sQbg*sqg + sqbg*sQg)-sqbg*sQg*(sQbg + sqg));
const double kin2s = (1/4.)*(sqbQ*(sQbg*sqg + sqbg*sQg) + sqbg*(sQbg*sqg + sQbg*sQg - 2.*sqg*sQg)
-2.*sqbQb*(sqbg*sQbg + sqg*sQg) + (sqbg*sqg + sQbg*sQg)*sqbq);
const double col1 = (N_C*N_C-1.)*(N_C*N_C-1.)/N_C;
const double col2 = (N_C*N_C-1.)/N_C;
return 16.*kin0*(col1*kin1s-col2*kin2s)/(sqbq*sQbQ*sqbg*sQbg*sqg*sQg);
}
//Matrix element for q Q -> q Q, including flux and averaging factors
//See (3.23) of MadFKS
double M_q_Q_to_q_Q(
CLHEP::HepLorentzVector const & pqin,
CLHEP::HepLorentzVector const & pQin,
CLHEP::HepLorentzVector const & pqout,
CLHEP::HepLorentzVector const & pQout
){
const double flux = 2.*(2.*pqin.dot(pQin));
constexpr double avg = omega_q*omega_q;
return M2_qb_Qb_q_Q(-pqin,-pQin,pqout,pQout)/(flux*avg);
}
//Matrix element for q Q -> q Q g, including flux and averaging factors
//See (3.23) of MadFKS
double M_q_Q_to_q_Q_g(
CLHEP::HepLorentzVector const & pqin,
CLHEP::HepLorentzVector const & pQin,
CLHEP::HepLorentzVector const & pqout,
CLHEP::HepLorentzVector const & pQout,
CLHEP::HepLorentzVector const & pgout
){
const double flux = 2.*(2.*pqin.dot(pQin));
constexpr double avg = omega_q*omega_q;
return M2_qb_Qb_q_Q_g(-pqin,-pQin,pqout,pQout,pgout)/(flux*avg);
}
///// 2-parton functions /////
//Implements C+S+V terms in (4.3) of MadFKS
double Regulated_2parton_q_Q_to_q_Q(
double mur2,
FKSConfig const & FKSParams,
CLHEP::HepLorentzVector const & pqin,
CLHEP::HepLorentzVector const & pQin,
CLHEP::HepLorentzVector const & pqout,
CLHEP::HepLorentzVector const & pQout
){
const double Shat = 2.*pqin.dot(pQin);
const double sqrtS = sqrt(Shat);
//We don't include Born contribution, (4.4) of MadFKS
//But we do include running coupling
const double tA= 0.;
//Collinear contribution, (4.5) of MadFKS
const double Eq = pqout.e();
const double EQ = pQout.e();
const double tC= Q_q_Q_to_q_Q(mur2, FKSParams, sqrtS, Eq, EQ);
//Soft contribution, (4.12) of MadFKS
//For eikonals we use all-outgoing flavour labels:
//i.e. kp=qbar, km=Qbar, kq=q, kQ=Q.
//We have factored out C_0 = T_F*T_F*(N_C*N_C-1.).
constexpr double C_qQbar = T_F*(N_C*N_C-2.)/N_C;
constexpr double C_qqbar = T_F*(-1.)/N_C;
constexpr double C_qQ = T_F*(-2.)/N_C;
double EikonalSum = -C_qQ*( Regulated_Integrated_Eikonal_mn(FKSParams, Shat, pqin, pQin)
+Regulated_Integrated_Eikonal_mn(FKSParams, Shat, pqout,pQout))
+C_qQbar*(Regulated_Integrated_Eikonal_mn(FKSParams, Shat, pQin, pqout)
+Regulated_Integrated_Eikonal_mn(FKSParams, Shat, pqin, pQout))
+C_qqbar*(Regulated_Integrated_Eikonal_mn(FKSParams, Shat, pqin, pqout)
+Regulated_Integrated_Eikonal_mn(FKSParams, Shat, pQin, pQout));
//Factor of 2 from (2-delta) in (3.24) of MadFKS
const double tS = 2.*EikonalSum;
//Virtual contribution (CDR Scheme)
const double tV = Regulated_Virtuals_q_Q_to_q_Q(mur2, FKSParams,pqin,pQin,pqout,pQout);
return (tA+tC+tS+tV)*M_q_Q_to_q_Q(pqin,pQin,pqout,pQout);
}
//Implements (A.6) of MadFKS (~(4.22) of FKS)
double Regulated_Integrated_Eikonal_mn(
FKSConfig const & FKSParams,
double Shat,
CLHEP::HepLorentzVector const & km,
CLHEP::HepLorentzVector const & kn
){
const double Q2 = pow(FKSParams.Q_ES,2);
const double xi_cut = FKSParams.xi_cut;
const double arg1 = xi_cut*xi_cut*Shat/Q2;
const double arg2 = km.dot(kn)/2./km.e()/kn.e();
//Using: log(1.-z)*log(z)= M_PI*M_PI/6.-Li2(z)-Li2(1-z)
return (1./2.)*pow(log(arg1),2)+(1./2.)*pow(log(arg2),2)
+log(arg1)*log(arg2)-M_PI*M_PI/6.
+gsl_sf_dilog(1-arg2);
}
//Nonsingular part of the virtual corrections, square bracket of (2.9) of ES
//Real part has been explicitly taken below, which explains the differeing pi^2 coeffs.
//Physical process: pqin + pQin -> pqout + pQout
//ES labelling: p1 + p2 -> p3 + p4
double Regulated_Virtuals_q_Q_to_q_Q(
double mu2,
FKSConfig const & FKSParams,
CLHEP::HepLorentzVector const & pqin,
CLHEP::HepLorentzVector const & pQin,
CLHEP::HepLorentzVector const & pqout,
CLHEP::HepLorentzVector const & pQout
){
const double Q2 = pow(FKSParams.Q_ES,2);
const double s= 2.*pqin.dot(pQin);
const double t= -2.*pqin.dot(pqout);
const double u= -2.*pqin.dot(pQout);
constexpr double col1 = (N_C*N_C-1)/2./N_C;
const double kin1 = -16. +log(-t/Q2)*(8.*log(s/Q2)-8.*log(-u/Q2)-2*log(-t/Q2)+6.)
-2.*(s*s-u*u)/(s*s+u*u)*(M_PI*M_PI +pow(log(-t/Q2)-log(s/Q2),2) + pow(log(-t/Q2)-log(-u/Q2),2) )
+2.*(s+u)/(s*s+u*u)*((s+u)*(log(-u/Q2)-log(s/Q2))+(u-s)*(2.*log(-t/Q2)-log(s/Q2)-log(-u/Q2)));
constexpr double col2= N_C;
const double kin2= 85./9.+M_PI*M_PI +2.*log(-t/Q2)*(log(-t/Q2)+log(-u/Q2)-2*log(s/Q2))
+(s*s-u*u)/(s*s+u*u)/2.*(1.0*M_PI*M_PI+2*pow(log(-t/Q2)-log(s/Q2),2) + pow(log(-t/Q2)-log(-u/Q2),2) )
-s*t/(s*s+u*u)*(log(-t/Q2)-log(-u/Q2))
+2.*u*t/(s*s+u*u)*(log(-t/Q2)-log(s/Q2))
+11./3.*(log(mu2/Q2)-log(-t/Q2));
constexpr double col3= 2.*N_F/3.;
const double kin3= (log(-t/Q2)-log(mu2/Q2)) -5./3.;
return col1*kin1 +col2*kin2 +col3*kin3;
}
//Helper functions for (5.6) of FKS/(4.6) of MadFKS
double Q_q(
FKSConfig const & FKSParams,
double sqrtS,
double Eq
){
const double Q2 = pow(FKSParams.Q_ES,2);
const double xi_cut = FKSParams.xi_cut;
const double delta_o = FKSParams.delta_o;
return gamma_prime_q
-log(sqrtS*sqrtS*delta_o/2./Q2)*(gamma_q-2.*C_F*log(2.*Eq/xi_cut/sqrtS))
+2.*C_F*( pow(log(2.*Eq/sqrtS),2)- pow(log(xi_cut),2) )
-2.*gamma_q*log(2.*Eq/sqrtS);
}
double Q_g(
FKSConfig const & FKSParams,
double sqrtS,
double Eg
){
const double Q2 = pow(FKSParams.Q_ES,2);
const double xi_cut = FKSParams.xi_cut;
const double delta_o = FKSParams.delta_o;
return gamma_prime_g
-log(sqrtS*sqrtS*delta_o/2./Q2)*(gamma_g-2.*C_A*log(2.*Eg/xi_cut/sqrtS))
+2.*C_A*( pow(log(2.*Eg/sqrtS),2)- pow(log(xi_cut),2) )
-2.*gamma_g*log(2.*Eg/sqrtS);
}
//Implement (5.6) of FKS
double Q_q_Q_to_q_Q(
double mur2,
FKSConfig const & FKSParams,
double sqrtS,
double Eq,
double EQ
){
const double Q2 = pow(FKSParams.Q_ES,2);
const double xi_cut = FKSParams.xi_cut;
return -log(mur2/Q2)*2.*(gamma_q +2.*C_F*log(xi_cut))
+ Q_q(FKSParams,sqrtS,Eq)
+ Q_q(FKSParams,sqrtS,EQ);
}
#if 0
//Implements Lorentz invariant form of (2.78) of 0709.2092
double diplus(
CLHEP::HepLorentzVector const & ka,
CLHEP::HepLorentzVector const & kb,
CLHEP::HepLorentzVector const & ki
){
double A = param_.fks.A_FNO;
double B = param_.fks.B_FNO;
const double sai = 2.*ka.dot(ki);
const double sbi = 2.*kb.dot(ki);
double fA = (sai+sbi)/4.;
double fB = 4.*sai/(pow(sai+sbi,2));
return pow(fA,A)*pow(fB,B);
}
//Implements Lorentz invariant form of (2.78) of 0709.2092
double diminus(
CLHEP::HepLorentzVector const & ka,
CLHEP::HepLorentzVector const & kb,
CLHEP::HepLorentzVector const & ki
){
double A = param_.fks.A_FNO;
double B = param_.fks.B_FNO;
const double sai = 2.*ka.dot(ki);
const double sbi = 2.*kb.dot(ki);
double fA = (sai+sbi)/4.;
double fB = 4.*sbi/(pow(sai+sbi,2));
return pow(fA,A)*pow(fB,B);
}
//Implements Lorentz invariant form of (2.75) of 0709.2092 with h=1
double dij(
CLHEP::HepLorentzVector const & ka,
CLHEP::HepLorentzVector const & kb,
CLHEP::HepLorentzVector const & ki,
CLHEP::HepLorentzVector const & kj
){
double A = param_.fks.A_FNO;
double B = param_.fks.B_FNO;
const double sab = 2.*ka.dot(kb);
const double sij = 2.*ki.dot(kj);
const double sai = 2.*ka.dot(ki);
const double sbi = 2.*kb.dot(ki);
const double saj = 2.*ka.dot(kj);
const double sbj = 2.*kb.dot(kj);
double fA = (sai+sbi)*(saj+sbj)/(4.*sab);
double fB = 2.*sab*sij/((sai+sbi)*(saj+sbj));
return pow(fA,A)*pow(fB,B);
}
void counter_event_final(
CLHEP::HepLorentzVector & k1c,
CLHEP::HepLorentzVector & k2c,
CLHEP::HepLorentzVector & k3c,
CLHEP::HepLorentzVector & k4c,
CLHEP::HepLorentzVector & k5c,
CLHEP::HepLorentzVector const & k1,
CLHEP::HepLorentzVector const & k2,
CLHEP::HepLorentzVector const & k3,
CLHEP::HepLorentzVector const & k4,
CLHEP::HepLorentzVector const & k5,
CLHEP::HepLorentzVector const & ki,
CLHEP::HepLorentzVector const & kj,
CLHEP::HepLorentzVector const & kk,
std::string type,
double zeta_col,
double xi_soft
) {
double sqrts_2 = k1.e();
double xi_i = ki.e()/sqrts_2;
double kihat_x = ki.px()/sqrts_2/xi_i;
double kihat_y = ki.py()/sqrts_2/xi_i;
double kihat_z = ki.pz()/sqrts_2/xi_i;
// We get the R matrix from the relation ki = p R where p = (0, 0, 1)
// | cos(theta) sin(phi)sin(theta) cos(phi)sin(theta) |
// R = | 0 cos(phi) -sin(phi) |
// | -sin(theta) sin(phi)cos(theta) cos(phi)cos(theta) |
// ki = p R = ( -sin(theta) sin(phi)cos(theta) cos(phi)cos(theta) )
double sintheta = -kihat_x;
double tanphi = kihat_y/kihat_z;
double phi = atan(tanphi);
double costheta = kihat_y/sin(phi);
double cosphi = cos(phi), sinphi = sin(phi);
double xi_j = kj.e()/sqrts_2;
double zeta_j = kj.pz()/(sqrts_2/xi_j);
double cosphi_j = kj.px()/(sqrts_2*xi_j*pow(1-zeta_j*zeta_j, 0.5));
double sinphi_j = kj.py()/(sqrts_2*xi_j*pow(1-zeta_j*zeta_j, 0.5));
double zeta = 1.0 - zeta_col;
// collinear p
double pj_x = pow(1.0-zeta*zeta, 0.5)*cosphi_j;
double pj_y = pow(1.0-zeta*zeta, 0.5)*sinphi_j;
double pj_z = zeta;
// deduce kj from kj = pR
double kj_e = sqrts_2*xi_j;
double kj_x = sqrts_2*xi_j*(pj_x*costheta - pj_z*sintheta);
double kj_y = sqrts_2*xi_j*(pj_x*sinphi*sintheta + pj_y*cosphi + pj_z*sinphi*costheta);
double kj_z = sqrts_2*xi_j*(pj_x*cosphi*sintheta - pj_y*sinphi + pj_z*cosphi*costheta);
// counter event: ka = pa, kb = pb
k1c = k1;
k2c = k2;
double kjt = pow(kj_x*kj_x + kj_y*kj_y, 0.5);
std::vector <double> mom_kj_unit = {kj_e/kjt, kj_x/kjt, kj_y/kjt, kj_z/kjt};
double kt1 = pow(ki.px()*ki.px() + ki.py()*ki.py(), 0.5);
double Y1 = (ki.e()+ki.pz())/kt1;
double Y2 = mom_kj_unit[0]+mom_kj_unit[3];
double cosphi1 = ki.px()/kt1, sinphi1 = ki.py()/kt1;
double cosphi2 = mom_kj_unit[1], sinphi2 = mom_kj_unit[2];
double A = 2.0*sqrts_2 - kt1*Y1, B = 2.0*sqrts_2 - kt1/Y1, C = -kt1*cosphi1, D = -kt1*sinphi1;
double kt2 = (C*C+D*D-B*A)/(2.0*C*cosphi2+2.0*D*sinphi2-B*Y2-A/Y2);
std::vector <double> mom_kj = {kt2*mom_kj_unit[0], kt2*mom_kj_unit[1], kt2*mom_kj_unit[2], kt2*mom_kj_unit[3]};
k3c = CLHEP::HepLorentzVector(ki.px(), ki.py(), ki.pz(), ki.e());
k4c = CLHEP::HepLorentzVector(kj.px(), kj.py(), kj.pz(), kj.e());
k5c = k1c+k2c-k3c-k4c;
}
//Implements (todo: refer to sum of integrals in Jeremy's notes)
double FSC_contribution(
double xi_s,
double zeta_cr,
CLHEP::HepLorentzVector const & k1,
CLHEP::HepLorentzVector const & k2,
CLHEP::HepLorentzVector const & k3,
CLHEP::HepLorentzVector const & k4,
CLHEP::HepLorentzVector const & kg,
CLHEP::HepLorentzVector const & ki,
CLHEP::HepLorentzVector const & kj,
CLHEP::HepLorentzVector const & kk
){
const double xi_cut= param_.fks.xi_cut;
const double delta_o = param_.fks.delta_o;
//Original event parameters:
const double sqrtS = sqrt(2.*k1.dot(k2));
const double xi_i = 2.*ki.e()/sqrtS;
const double zeta_ij = 1.-ki.dot(kj)/(ki.e()*kj.e());
// Counter events - soft, collinear +, collinear -, soft and collinear +, soft and collinear -
CLHEP::HepLorentzVector k1s, k2s, k3s, k4s, k5s;
CLHEP::HepLorentzVector k1c, k2c, k3c, k4c, k5c;
CLHEP::HepLorentzVector k1sc, k2sc, k3sc, k4sc, k5sc;
counter_event_final(k1s, k2s, k3s, k4s, k5s, k1, k2, k3, k4, kg, ki, kj, kk, "soft", 1e-9, 1e-9);
counter_event_final(k1c, k2c, k3c, k4c, k5c, k1, k2, k3, k4, kg, ki, kj, kk, "col", 1e-9, 1e-9);
counter_event_final(k1sc, k2sc, k3sc, k4sc, k5sc, k1, k2, k3, k4, kg, ki, kj, kk, "soft_col", 1e-9, 1e-9);
//Damped counterterms:
const double ME2dampedFS = xi_i*xi_i*(1-zeta_ij)*MatrixElement::M2_qb_Qb_q_Q_g(-k1, -k2, k4, k3, kg);
const double ME2dampedFS_s = xi_s*xi_s*(1-zeta_ij)*MatrixElement::M2_qb_Qb_q_Q_g(-k1s, -k2s, k4s, k3s, k5s);
const double ME2dampedFS_c = xi_i*xi_i*(1-zeta_cr)*MatrixElement::M2_qb_Qb_q_Q_g(-k1c, -k2c, k4c, k3c, k5c);
const double ME2dampedFS_sc = xi_s*xi_s*(1-zeta_cr)*MatrixElement::M2_qb_Qb_q_Q_g(-k1sc, -k2sc, k4sc, k3sc, k5sc);
//Add the contributions
if(xi_i>xi_cut){
double sum = 0.;
if(zeta_ij<1.-delta_o){
sum += ME2dampedFS;
}
if(zeta_ij>1.-delta_o){
sum += ME2dampedFS-ME2dampedFS_c;
}
return sum/(1-zeta_ij)/xi_i/xi_i;
} else{ //xig<xicut
double sum = 0.;
if(zeta_ij<1.-delta_o){
sum += ME2dampedFS-ME2dampedFS_s;
}
if(zeta_ij>1.-delta_o){
sum += ME2dampedFS-ME2dampedFS_s-ME2dampedFS_c+ME2dampedFS_sc;
}
return sum/(1-zeta_ij)/xi_i/xi_i;
}
}
void counter_event_initial(
CLHEP::HepLorentzVector & k1c,
CLHEP::HepLorentzVector & k2c,
CLHEP::HepLorentzVector & k3c,
CLHEP::HepLorentzVector & k4c,
CLHEP::HepLorentzVector & k5c,
CLHEP::HepLorentzVector const & k1,
CLHEP::HepLorentzVector const & k2,
CLHEP::HepLorentzVector const & k3,
CLHEP::HepLorentzVector const & k4,
CLHEP::HepLorentzVector const & k5,
std::string type,
double zeta_col,
double xi_soft
){
k1c = k1;
k2c = k2;
double sqrts_2 = k1.e();
double xi_i = k5.e()/sqrts_2;
double zeta_i = k5.pz()/(sqrts_2*xi_i);
double zeta_p = 1.-zeta_col;
double zeta_m = -1.+zeta_col;
double r_xi = (xi_soft/xi_i);
double r_z_p = pow(1-zeta_p*zeta_p, 0.5)/pow(1-zeta_i*zeta_i, 0.5);
double r_z_m = pow(1-zeta_m*zeta_m, 0.5)/pow(1-zeta_i*zeta_i, 0.5);
std::vector <double> momentum_gluon = {0.0, 0.0, 0.0, 0.0};
if (type == "soft") {
momentum_gluon = {r_xi*k5.e(), r_xi*k5.px(), r_xi*k5.py(), r_xi*k5.pz()};
}
if (type == "soft_col_p1") {
momentum_gluon = {sqrts_2*xi_soft, r_xi*r_z_p*k5.px(), r_xi*r_z_p*k5.py(), sqrts_2*xi_soft*zeta_p};
}
if (type == "soft_col_m1") {
momentum_gluon = {sqrts_2*xi_soft, r_xi*r_z_m*k5.px(), r_xi*r_z_m*k5.py(), sqrts_2*xi_soft*zeta_m};
}
if (type == "col_p1") {
momentum_gluon = {sqrts_2*xi_i, r_z_p*k5.px(), r_z_p*k5.py(), sqrts_2*xi_i*zeta_p};
}
if (type == "col_m1") {
momentum_gluon = {sqrts_2*xi_i, r_z_m*k5.px(), r_z_m*k5.py(), sqrts_2*xi_i*zeta_m};
}
k5c = CLHEP::HepLorentzVector(momentum_gluon[1], momentum_gluon[2], momentum_gluon[3], momentum_gluon[0]);
CLHEP::HepLorentzVector P = k1c+k2c-k5c;
double PX = P.px();
double PY = P.py();
double kt1 = pow(k3.px()*k3.px()+k3.py()*k3.py(), 0.5);
double phi1 = acos(k3.px()/kt1);
double phi2 = atan((PY-kt1*sin(phi1))/(PX-kt1*cos(phi1)));
double kt2 = (PX-kt1*cos(phi1))/cos(phi2);
double Pp = P.e()+P.pz();
double Pm = P.e()-P.pz();
double a = kt2*Pm;
double b = kt1*kt1-kt2*kt2-Pp*Pm;
double c = Pp*kt2;
double Y2 = (-b+sqrt(b*b-4.0*a*c))/(2.0*a);
double y2 = log(Y2);
double Y1 = kt1/(Pm-kt2/Y2);
double y1 = log(Y1);
k3c = CLHEP::HepLorentzVector(kt1*cos(phi1), kt1*sin(phi1), kt1*sinh(y1), kt1*cosh(y1));
k4c = CLHEP::HepLorentzVector(kt2*cos(phi2), kt2*sin(phi2), kt2*sinh(y2), kt2*cosh(y2));
}
//Implements (todo: refer to sum of integrals in Jeremy's notes)
double ISC_contribution(
double xi_s,
double zeta_c,
CLHEP::HepLorentzVector const & k1,
CLHEP::HepLorentzVector const & k2,
CLHEP::HepLorentzVector const & k3,
CLHEP::HepLorentzVector const & k4,
CLHEP::HepLorentzVector const & k5,
int j
){
const double xi_cut= param_.fks.xi_cut;
const double delta_i = param_.fks.delta_i;
const double sqrtS = sqrt(2.*k1.dot(k2));
//Assign the momentum kj of region Sj
CLHEP::HepLorentzVector kj;
if (j==3){
kj = k3;
} else if (j==4){
kj = k4;
} else if (j==5){
kj = k5;
} else {
throw std::logic_error("int j not in {3,4,5}");
}
//Parton j original parameters:
const double xi_j = 2.*kj.e()/sqrtS;
const double zeta_j = kj.z()/kj.e();
// Counter events - soft, collinear +, collinear -, soft and collinear +, soft and collinear -
CLHEP::HepLorentzVector k1s, k2s, k3s, k4s, k5s;
CLHEP::HepLorentzVector k1cp, k2cp, k3cp, k4cp, k5cp;
CLHEP::HepLorentzVector k1cm, k2cm, k3cm, k4cm, k5cm;
CLHEP::HepLorentzVector k1scp, k2scp, k3scp, k4scp, k5scp;
CLHEP::HepLorentzVector k1scm, k2scm, k3scm, k4scm, k5scm;
counter_event_initial(k1s, k2s, k3s, k4s, k5s, k1, k2, k3, k4, k5, "soft", 1e-9, 1e-9);
counter_event_initial(k1cp, k2cp, k3cp, k4cp, k5cp, k1, k2, k3, k4, k5, "col_p1", 1e-9, 1e-9);
counter_event_initial(k1cm, k2cm, k3cm, k4cm, k5cm, k1, k2, k3, k4, k5, "col_m1", 1e-9, 1e-9);
counter_event_initial(k1scp, k2scp, k3scp, k4scp, k5scp, k1, k2, k3, k4, k5, "soft_col_p1", 1e-9, 1e-9);
counter_event_initial(k1scm, k2scm, k3scm, k4scm, k5scm, k1, k2, k3, k4, k5, "soft_col_m1", 1e-9, 1e-9);
//Initialise counterterms:
double ME2dampedIS_f=0;
double ME2dampedIS_s=0;
double ME2dampedIS_cp=0;
double ME2dampedIS_cm=0;
double ME2dampedIS_scp=0;
double ME2dampedIS_scm=0;
//Construct counterterms, placing kjbar in appropriate slot:
if (j==3){
ME2dampedIS_f = xi_j*xi_j*(1.-zeta_j*zeta_j)*MatrixElement::M2_qb_Qb_q_Q_g(-k1,-k2,k5,k3,k4);
ME2dampedIS_s = xi_s*xi_s*(1.-zeta_j*zeta_j)*MatrixElement::M2_qb_Qb_q_Q_g(-k1s,k2s,k5s,k3s,k4s);
ME2dampedIS_cp = xi_j*xi_j*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1cp,-k2cp,k5cp,k3cp,k4cp);
ME2dampedIS_cm = xi_j*xi_j*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1cm,-k2cm,k5cm,k3cm,k4cm);
ME2dampedIS_scp = xi_s*xi_s*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1scp,-k2scp,k5scp,k3scp,k4scp);
ME2dampedIS_scm = xi_s*xi_s*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1scm,-k2scm,k5scm,k3scm,k4scm);
} else if (j==4){
ME2dampedIS_f = xi_j*xi_j*(1.-zeta_j*zeta_j)*MatrixElement::M2_qb_Qb_q_Q_g(-k1,-k2,k4,k5,k3);
ME2dampedIS_s = xi_s*xi_s*(1.-zeta_j*zeta_j)*MatrixElement::M2_qb_Qb_q_Q_g(-k1s,-k2s,k4s,k5s,k3s);
ME2dampedIS_cp = xi_j*xi_j*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1cp,-k2cp,k4cp,k5cp,k3cp);
ME2dampedIS_cm = xi_j*xi_j*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1cm,-k2cm,k4cm,k5cm,k3cm);
ME2dampedIS_scp = xi_s*xi_s*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1scp,-k2scp,k4scp,k5scp,k3scp);
ME2dampedIS_scm = xi_s*xi_s*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1scm,-k2scm,k4scm,k5scm,k3scm);
} else if (j==5){
ME2dampedIS_f = xi_j*xi_j*(1.-zeta_j*zeta_j)*MatrixElement::M2_qb_Qb_q_Q_g(-k1,-k2,k4,k3,k5);
ME2dampedIS_s = xi_s*xi_s*(1.-zeta_j*zeta_j)*MatrixElement::M2_qb_Qb_q_Q_g(-k1s,-k2s,k4s,k3s,k5s);
ME2dampedIS_cp = xi_j*xi_j*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1cp,-k2cp,k4cp,k3cp,k5cp);
ME2dampedIS_cm = xi_j*xi_j*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1cm,-k2cm,k4cm,k3cm,k5cm);
ME2dampedIS_scp = xi_s*xi_s*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1scp,-k2scp,k4scp,k3scp,k5scp);
ME2dampedIS_scm = xi_s*xi_s*(1.-zeta_c*zeta_c)*MatrixElement::M2_qb_Qb_q_Q_g(-k1scm,-k2scm,k4scm,k3scm,k5scm);
} else {
throw std::logic_error("int j not in {3,4,5}");
}
//Add the contributions
if(xi_j>xi_cut){
double sum = 0.;
if(zeta_j<1.-delta_i){
sum += ME2dampedIS_f/(1-zeta_j);
}
if(zeta_j>-1.+delta_i){
sum += ME2dampedIS_f/(1+zeta_j);
}
if(zeta_j>1.-delta_i){
sum += (ME2dampedIS_f - ME2dampedIS_cp)/(1-zeta_j) ;//ME2dampedIS - Mtilde(1,xij)
}
if(zeta_j<-1.+delta_i){
sum += (ME2dampedIS_f - ME2dampedIS_cm)/(1+zeta_j); //ME2dampedIS - Mtilde(-1,xij)
}
return sum/xi_j/xi_j;
} else{ //xij<xicut
double sum = 0.;
if(zeta_j<1.-delta_i){
sum += (ME2dampedIS_f - ME2dampedIS_s)/(1-zeta_j) ; //ME2dampedIS - Mtilde(zetaj,0)
}
if(zeta_j>-1.+delta_i){
sum += (ME2dampedIS_f - ME2dampedIS_s)/(1+zeta_j) ; //ME2dampedIS - Mtilde(zetaj,0)
}
if(zeta_j>1.-delta_i){
sum += (ME2dampedIS_f - ME2dampedIS_cp- ME2dampedIS_s + ME2dampedIS_scp)/(1-zeta_j) ; //ME2dampedIS - Mtilde(1,xij)- Mtilde(zetaj,0) + Mtilde(1,0)
}
if(zeta_j<-1.+delta_i){
sum += (ME2dampedIS_f - ME2dampedIS_cm- ME2dampedIS_s + ME2dampedIS_scm)/(1+zeta_j) ; //ME2dampedIS - Mtilde(-1,xij)- Mtilde(zetaj,0) + Mtilde(-1,0)
}
return sum/xi_j/xi_j;
}
}
//FKS treatment
double ME2_NLO_FKS(Event const & ev) {
std::cout<<"ME2_NLO_FKS"<<std::endl;
auto const & incoming = ev.incoming();
//Avoid tag_extremal_jet_partons.
auto out_partons = filter_partons(ev.outgoing());
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
//Calculate boost required to partonic COM frame
const double Ea= pa.e();
const double Eb= pb.e();
double beta_to_COM = (1-(Eb/Ea))/(1+(Eb/Ea));
//Perform boost for incoming partons
CLHEP::HepLorentzVector k1=pb;
CLHEP::HepLorentzVector k2=pa;
k1.boostZ(beta_to_COM);
k2.boostZ(beta_to_COM);
//Strong coupling, use central renormalisation scale
const double alpha_s = alpha_s_(ev.central().mur);
const double mur2 = pow(ev.central().mur,2);
if(out_partons.size()==2){
const auto p1 = to_HepLorentzVector(out_partons.front());
const auto pn = to_HepLorentzVector(out_partons.back());
//Change to FKS labelling of partons
CLHEP::HepLorentzVector k3=p1;
CLHEP::HepLorentzVector k4=pn;
//Boost along z direction to partonic COM frame
k3.boostZ(beta_to_COM);
k4.boostZ(beta_to_COM);
return alpha_s*Regulated_2parton_ME2(mur2,k1,k2,k3,k4);
} else if(out_partons.size()==3){
//Define the momenta based on event classification.
CLHEP::HepLorentzVector pq;
CLHEP::HepLorentzVector pQ;
CLHEP::HepLorentzVector pg;
if (ev.type()==event_type::FKL){
pq = to_HepLorentzVector(out_partons.at(0));
pg = to_HepLorentzVector(out_partons.at(1));
pQ = to_HepLorentzVector(out_partons.at(2));
} else if (ev.type()==event_type::unordered_backward){
//std::cout<<"unob"<<std::endl;
pg = to_HepLorentzVector(out_partons.at(0));
pq = to_HepLorentzVector(out_partons.at(1));
pQ = to_HepLorentzVector(out_partons.at(2));
} else if (ev.type()==event_type::unordered_forward){
//std::cout<<"unof"<<std::endl;
pq = to_HepLorentzVector(out_partons.at(0));
pQ = to_HepLorentzVector(out_partons.at(1));
pg = to_HepLorentzVector(out_partons.at(2));
} else {
return 0.;
}
//ki = pi in partonic COM frame.
CLHEP::HepLorentzVector kq=pq;
CLHEP::HepLorentzVector kQ=pQ;
CLHEP::HepLorentzVector kg=pg;
kq.boostZ(beta_to_COM);
kQ.boostZ(beta_to_COM);
kg.boostZ(beta_to_COM);
// //Print event info:
// std::cout << "After boost" << std::endl;
// std::cout<<"p1= "<< k1<<std::endl;
// std::cout<<"p2= "<< k2<<std::endl;
// std::cout<<"pq= "<< kq<<std::endl;
// std::cout<<"pQ= "<< kQ<<std::endl;
// std::cout<<"pg= "<< kg<<std::endl;
//Measurement functions for (pb + pa -> pq + pQ + pg )
// (2.67) and (2.68) of 0709.2092, expressed in terms of invariants.
const double dgplus = diplus(pa, pb, pg);
const double dqplus = diplus(pa, pb, pq);
const double dQplus = diplus(pa, pb, pQ);
const double dgminus = diminus(pa, pb, pg);
const double dqminus = diminus(pa, pb, pq);
const double dQminus = diminus(pa, pb, pQ);
const double dgq = dij(pa, pb, pg, pq);
const double dgQ = dij(pa, pb, pg, pQ);
const double Denom = (1/dgplus+1/dqplus+1/dQplus+1/dgminus+1/dqminus+1/dQminus+1/dgq+1/dgQ);
const double Sgplus = (1/Denom)*(1/dgplus);
const double Sqplus = (1/Denom)*(1/dqplus);
const double SQplus = (1/Denom)*(1/dQplus);
const double Sgminus = (1/Denom)*(1/dgminus);
const double Sqminus = (1/Denom)*(1/dqminus);
const double SQminus = (1/Denom)*(1/dQminus);
const double Sgq= (1/Denom)*(1/dgq);
const double SgQ= (1/Denom)*(1/dgQ);
//Counterevents, numerical 0 and 1:
//const double xi_s = param_.fks.xi_s;
double xi_s = 1e-9;
const double zeta_c = 1.-xi_s;
double ME2_3parton= 0.;
//Sum FSC contributions (kg||kq), (kg||kQ).
ME2_3parton += Sgq*FSC_contribution(xi_s, zeta_c, k1, k2, kq, kQ, kg, kg, kq, kQ);
ME2_3parton += SgQ*FSC_contribution(xi_s, zeta_c, k1, k2, kq, kQ, kg, kg, kQ, kq);
//Sum ISC contributions:
ME2_3parton += (Sgplus+Sgminus)*ISC_contribution(xi_s, zeta_c, k1, k2, kq, kQ, kg, 5);
ME2_3parton += (Sqplus+Sqminus)*ISC_contribution(xi_s, zeta_c, k1, k2, kq, kQ, kg, 3);
ME2_3parton += (SQplus+SQminus)*ISC_contribution(xi_s, zeta_c, k1, k2, kq, kQ, kg, 4);
return alpha_s*ME2_3parton/omega_q/omega_q; //Divide by 2s here?
//End 3-parton treatment
} else{
throw std::logic_error("FKS treatment currently supports only 2 or 3 partons in final state.");
}
}
#endif
} }

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:12 AM (15 h, 23 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6549950
Default Alt Text
FKS.cc (31 KB)

Event Timeline

HEPForge.org · .