diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 704acd8..ea3ce00 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1680 +1,1630 @@
#include "RHEJ/MatrixElement.hh"
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "RHEJ/Constants.hh"
#include "RHEJ/currents.hh"
#include "RHEJ/PDG_codes.hh"
#include "RHEJ/uno.hh"
#include "RHEJ/qqx.hh"
#include "RHEJ/utility.hh"
namespace RHEJ{
//cf. last line of eq. (22) in \ref Andersen:2011hs
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j, double lambda
) const {
const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
// use alpha_s(sqrt(q_j*lambda)), evolved to mur
return (
1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda))
)*result;
}
double MatrixElement::virtual_corrections(
double mur,
std::array<Particle, 2> const & in,
std::vector<Particle> const & out
) const{
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
size_t first_idx = 0;
size_t last_idx = out.size() - 1;
// if there is a Higgs or unordered gluon outside the extremal partons
// then it is not part of the FKL ladder and does not contribute
// to the virtual corrections
if(out.front().type == pid::Higgs || has_unob_gluon(in, out)){
q -= out[1].p;
++first_idx;
}
if(out.back().type == pid::Higgs || has_unof_gluon(in, out)){
--last_idx;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q, CLAMBDA)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| has_unof_gluon(in, out)
);
return exp(exponent);
}
} // namespace RHEJ
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
// cout << "#Fadin qa : "<<qav<<endl;
// cout << "#Fadin qb : "<<qbv<<endl;
// cout << "#Fadin p1 : "<<p1<<endl;
// cout << "#Fadin p2 : "<<p2<<endl;
// cout << "#Fadin p5 : "<<p5<<endl;
// cout << "#Fadin Gauge Check : "<< CL.dot(p5)<<endl;
// cout << "#Fadin C2L : "<< -CL.dot(CL)<<" "<<-CL.dot(CL)/(qav.m2()*qbv.m2())/(4./p5.perp2())<<endl;
// TODO can this dead test go?
// if (-CL.dot(CL)<0.)
// if (fabs(CL.dot(p5))>fabs(CL.dot(CL))) // not sufficient!
// return 0.;
// else
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
double kperp=(qav-qbv).perp();
if (kperp>RHEJ::CLAMBDA)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
}
//! Lipatov vertex
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,
CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
double kperp=(qav-qbv).perp();
if (kperp>RHEJ::CLAMBDA)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
return temp;
}
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
if (aptype==21&&bptype==21) {
return jM2gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2qg(pn,pb,p1,pa);
else
return jM2qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jM2qg(p1,pa,pn,pb);
else
return jM2qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2qQ(pn,pb,p1,pa);
else
return jM2qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return jM2qQbar(p1,pa,pn,pb);
else
return jM2qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// We know it cannot be gg incoming.
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jMWqg(pn,pl,plbar,pb,p1,pa);
else
return jMWqbarg(pn,pl,plbar,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jMWqg(p1,pl,plbar,pa,pn,pb);
else
return jMWqbarg(p1,pl,plbar,pa,pn,pb);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (bptype>0) {
if (aptype>0)
return jMWqQ(pn,pl,plbar,pb,p1,pa);
else
return jMWqQbar(pn,pl,plbar,pb,p1,pa);
}
else {
if (aptype>0)
return jMWqbarQ(pn,pl,plbar,pb,p1,pa);
else
return jMWqbarQbar(pn,pl,plbar,pb,p1,pa);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return jMWqQ(p1,plbar,pl,pa,pn,pb);
else
return jMWqQbar(p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0)
return jMWqbarQ(p1,plbar,pl,pa,pn,pb);
else
return jMWqbarQbar(p1,plbar,pl,pa,pn,pb);
}
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for backwards uno tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*/
double ME_W_unob_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (bptype == 21 && aptype != 21) { // b gluon => W emission off a
if (aptype > 0)
return jM2Wunogqg(pg,p1,plbar,pl,pa,pn,pb);
else
return jM2Wunogqbarg(pg,p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return junobMWqQg(pn,plbar,pl,pb,p1,pa,pg);
else
return junobMWqQbarg(pn,plbar,pl,pb,p1,pa,pg);
}
else{
if (aptype>0)
return junobMWqbarQg(pn,plbar,pl,pb,p1,pa,pg);
else
return junobMWqbarQbarg(pn,plbar,pl,pb,p1,pa,pg);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return jM2WunogqQ(pg,p1,plbar,pl,pa,pn,pb);
else //qqbar
return jM2WunogqQbar(pg,p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return jM2WunogqbarQ(pg,p1,plbar,pl,pa,pn,pb);
else //qbarqbar
return jM2WunogqbarQbar(pg,p1,plbar,pl,pa,pn,pb);
}
}
}
}
/** Matrix element squared for uno forward tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @returns ME Squared for unof Tree-Level Current-Current Scattering
*/
double ME_W_unof_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (aptype==21 && bptype!=21) {//a gluon => W emission off b
if (bptype > 0)
return jM2Wunogqg(pg, pn,plbar, pl, pb, p1, pa);
else
return jM2Wunogqbarg(pg, pn,plbar, pl, pb, p1, pa);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return jM2WunogqQ(pg,pn,plbar,pl,pb,p1,pa);
else
return jM2WunogqQbar(pg,pn,plbar,pl,pb,p1,pa);
}
else{
if (aptype>0)
return jM2WunogqbarQ(pg,pn,plbar,pl,pb,p1,pa);
else
return jM2WunogqbarQbar(pg,pn,plbar,pl,pb,p1,pa);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return junofMWgqQ(pg,pn,pb,p1,plbar,pl,pa);
else //qqbar
return junofMWgqQbar(pg,pn,pb,p1,plbar,pl,pa);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return junofMWgqbarQ(pg,pn,pb,p1,plbar,pl,pa);
else //qbarqbar
return junofMWgqbarQbar(pg,pn,pb,p1,plbar,pl,pa);
}
}
}
}
/** \brief Matrix element squared for backward qqx tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*/
double ME_W_qqxb_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFbackward;
if (pqbar.rapidity() > pq.rapidity()){
swapQuarkAntiquark=true;
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pq.minus())+(pq.minus())/pa.minus())+1./3.)*3./4.;
}
else{
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pqbar.minus())+(pqbar.minus())/pa.minus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark){
return jM2Wggtoqqbarg(pa, pqbar, plbar, pl, pq, pn,pb)*CFbackward;}
else {
return jM2Wggtoqbarqg(pa, pq, plbar, pl, pqbar, pn,pb)*CFbackward;}
}
else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
if (wc!=1){ // W Emitted from backwards qqx
if (swapQuarkAntiquark){
return jM2WgQtoqqbarQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
else{
return jM2WgQtoqbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
}
else { // W Must be emitted from forwards leg.
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl)*CFbackward;}
else{
return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl)*CFbackward;}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxb");
}
}
/* \brief Matrix element squared for forward qqx tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param p1 Final state 1 Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @returns ME Squared for qqxf Tree-Level Current-Current Scattering
*/
double ME_W_qqxf_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFforward;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pq.plus())+(pq.plus())/pb.plus())+1./3.)*3./4.;
}
else{
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pqbar.plus())+(pqbar.plus())/pb.plus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark){
return jM2Wggtoqqbarg(pb, pqbar, plbar, pl, pq, p1,pa)*CFforward;}
else {
return jM2Wggtoqbarqg(pb, pq, plbar, pl, pqbar, p1,pa)*CFforward;}
}
else if (bptype==21&&aptype!=21) {// b gluon => W emission off a or qqx
if (wc==1){ // W Emitted from forwards qqx
if (swapQuarkAntiquark){
return jM2WgQtoqbarqQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
else {
return jM2WgQtoqqbarQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
}
// W Must be emitted from backwards leg.
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl)*CFforward;}
else{
return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl)*CFforward;}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxf");
}
}
/* \brief Matrix element squared for central qqx tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqxpair
* @param nbelow Number of gluons emitted after central qqxpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqx
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_W_qqxmid_current(
int aptype, int bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<HLV> partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
bool swapQuarkAntiquark=false;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
}
double CFforward = (0.5*(3.-1./3.)*(pb.plus()/(partons[partons.size()-1].plus())+(partons[partons.size()-1].plus())/pb.plus())+1./3.)*3./4.;
double CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(partons[0].minus())+(partons[0].minus())/pa.minus())+1./3.)*3./4.;
double wt=1.;
if (aptype==21) wt*=CFbackward;
if (bptype==21) wt*=CFforward;
if (aptype <=0 && bptype <=0){ // Both External AntiQuark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons,true,true, swapQuarkAntiquark, nabove, nbelow);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, false);
}
} // end both antiquark
else if (aptype<=0){ // a is antiquark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove, nbelow);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons,false,true, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove, nbelow, false);
}
} // end a is antiquark
else if (bptype<=0){ // b is antiquark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, false);
}
} //end b is antiquark
else{ //Both Quark or gluon
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl, plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow);}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, false);
}
}
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype==21) // gg initial state
return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered forward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param punof Unordered Particle Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered forward emission
*/
double ME_Higgs_current_unof(
int aptype, int bptype,
CLHEP::HepLorentzVector const & punof,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (aptype>0)
return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param punob Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*/
double ME_Higgs_current_unob(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & punob,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (bptype>0)
return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(RHEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(RHEJ::MatrixElementConfig const & config) {
#ifndef RHEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace anonymous
namespace RHEJ{
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
double MatrixElement::operator()(
double mur,
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
std::unordered_map<int, std::vector<Particle>> const & decays,
bool check_momenta
) const {
return tree(
mur,
incoming, outgoing, decays,
check_momenta
)*virtual_corrections(
mur,
incoming, outgoing
);
}
double MatrixElement::tree_kin(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
std::unordered_map<int, std::vector<Particle>> const & decays,
bool check_momenta
) const {
assert(
std::is_sorted(
incoming.begin(), incoming.end(),
[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
)
);
assert(std::is_sorted(outgoing.begin(), outgoing.end(), rapidity_less{}));
auto AWZH_boson = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(outgoing)){
return tree_kin_jets(incoming, outgoing, check_momenta);
}
switch(AWZH_boson->type){
case pid::Higgs: {
return tree_kin_Higgs(incoming, outgoing, check_momenta);
}
// TODO
case pid::Wp: {
return tree_kin_W(incoming, outgoing, decays, true, check_momenta);
}
case pid::Wm: {
return tree_kin_W(incoming, outgoing, decays, false, check_momenta);
}
case pid::photon:
case pid::Z:
default:
throw std::logic_error("Emission of boson of unsupported type.");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
std::array<Particle, 2> const & incoming,
std::vector<Particle> out_partons, bool check_momenta
) const{
if(!check_momenta){
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param_.jet_param.def);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission
if(has_unob_gluon(incoming, out_partons)){
assert(jets.size() >= 3);
++most_backward;
}
else if(has_unof_gluon(incoming, out_partons)){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = cs.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(RHEJ::is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
extremal_jet_idx:
no_extremal_jet_idx;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
double MatrixElement::tree_kin_jets(
std::array<Particle, 2> const & incoming,
std::vector<Particle> partons,
bool check_momenta
) const {
partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
if(has_unob_gluon(incoming, partons) || has_unof_gluon(incoming, partons)){
throw std::logic_error("unordered emission not implemented for pure jets");
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4*(N_C*N_C - 1))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn
);
}
double MatrixElement::tree_kin_W(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
std::unordered_map<int, std::vector<Particle>> const & decays,
bool WPlus,
bool check_momenta
) const {
HLV plbar, pl;
for (auto& x: decays) {
if (x.second.at(0).type < 0){
plbar = to_HepLorentzVector(x.second.at(0));
pl = to_HepLorentzVector(x.second.at(1));
}
else{
pl = to_HepLorentzVector(x.second.at(0));
plbar = to_HepLorentzVector(x.second.at(1));
}
}
+
+ const auto the_W = std::find_if(
+ begin(outgoing), end(outgoing),
+ [](Particle const & s){ return abs(s.type) == pid::Wp; }
+ );
+
+ std::vector<Particle> partons(begin(outgoing), the_W);
+ partons.insert(end(partons), the_W + 1, end(outgoing));
+ partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
+
if(has_unob_gluon(incoming, outgoing)){
- return tree_kin_W_unob(incoming, outgoing, plbar, pl, WPlus, check_momenta);
+ return tree_kin_W_unob(incoming, partons, plbar, pl, WPlus, check_momenta);
}
else if(has_unof_gluon(incoming, outgoing)){
- return tree_kin_W_unof(incoming, outgoing, plbar, pl, WPlus, check_momenta);
+ return tree_kin_W_unof(incoming, partons, plbar, pl, WPlus, check_momenta);
}
else if(has_Ex_qqxb(incoming, outgoing)){
- return tree_kin_W_qqxb(incoming, outgoing, plbar, pl, WPlus, check_momenta);
+ return tree_kin_W_qqxb(incoming, partons, plbar, pl, WPlus, check_momenta);
}
else if(has_Ex_qqxf(incoming, outgoing)){
- return tree_kin_W_qqxf(incoming, outgoing, plbar, pl, WPlus, check_momenta);
+ return tree_kin_W_qqxf(incoming, partons, plbar, pl, WPlus, check_momenta);
}
else if(has_mid_qqx(outgoing)){
- return tree_kin_W_qqxmid(incoming, outgoing, plbar, pl, WPlus, check_momenta);
+ return tree_kin_W_qqxmid(incoming, partons, plbar, pl, WPlus, check_momenta);
}
else{
- return tree_kin_W_FKL(incoming, outgoing, plbar, pl, WPlus, check_momenta);
+ return tree_kin_W_FKL(incoming, partons, plbar, pl, WPlus, check_momenta);
}
}
double MatrixElement::tree_kin_W_FKL(
std::array<Particle, 2> const & incoming,
- std::vector<Particle> const & outgoing,
+ std::vector<Particle> const & partons,
HLV plbar, HLV pl,
bool WPlus, bool check_momenta
) const {
- const auto the_W = std::find_if(
- begin(outgoing), end(outgoing),
- [](Particle const & s){ return abs(s.type) == pid::Wp; }
- );
-
- const auto pW = to_HepLorentzVector(*the_W);
- std::vector<Particle> partons(begin(outgoing), the_W);
- partons.insert(end(partons), the_W + 1, end(outgoing));
- partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
-
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
bool wc;
if (incoming[0].type==partons[0].type) { //leg b emits w
wc = true;}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, pl, plbar, wc
);
}
else{
current_factor = ME_W_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*9./8.*ladder_factor;
}
double MatrixElement::tree_kin_W_unob(
std::array<Particle, 2> const & incoming,
- std::vector<Particle> const & outgoing,
+ std::vector<Particle> const & partons,
HLV plbar, HLV pl,
bool WPlus, bool check_momenta
) const {
- const auto the_W = std::find_if(
- begin(outgoing), end(outgoing),
- [](Particle const & s){ return abs(s.type) == pid::Wp; }
- );
-
- const auto pW = to_HepLorentzVector(*the_W);
- std::vector<Particle> partons(begin(outgoing), the_W);
- partons.insert(end(partons), the_W + 1, end(outgoing));
- partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
-
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto pg = to_HepLorentzVector(partons[0]);
auto p1 = to_HepLorentzVector(partons[1]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1- pg;
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc;
if (incoming[0].type==partons[1].type) { //leg b emits w
wc = true;}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_unob_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, pg, pl, plbar, wc
);
}
else{
current_factor = ME_W_unob_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, pg, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*9./8.*ladder_factor;
}
double MatrixElement::tree_kin_W_unof(
std::array<Particle, 2> const & incoming,
- std::vector<Particle> const & outgoing,
+ std::vector<Particle> const & partons,
HLV plbar, HLV pl,
bool WPlus, bool check_momenta
) const {
- const auto the_W = std::find_if(
- begin(outgoing), end(outgoing),
- [](Particle const & s){ return abs(s.type) == pid::Wp; }
- );
-
- const auto pW = to_HepLorentzVector(*the_W);
- std::vector<Particle> partons(begin(outgoing), the_W);
- partons.insert(end(partons), the_W + 1, end(outgoing));
- partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
-
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 2]);
auto pg = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc;
if (incoming[0].type==partons[0].type) { //leg b emits w
wc = true;}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_unof_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, pg, pl, plbar, wc
);
}
else{
current_factor = ME_W_unof_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, pg, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*9./8.*ladder_factor;
}
double MatrixElement::tree_kin_W_qqxb(
std::array<Particle, 2> const & incoming,
- std::vector<Particle> const & outgoing,
+ std::vector<Particle> const & partons,
HLV plbar, HLV pl,
bool WPlus, bool check_momenta
) const {
- const auto the_W = std::find_if(
- begin(outgoing), end(outgoing),
- [](Particle const & s){ return abs(s.type) == pid::Wp; }
- );
-
- const auto pW = to_HepLorentzVector(*the_W);
- std::vector<Particle> partons(begin(outgoing), the_W);
- partons.insert(end(partons), the_W + 1, end(outgoing));
- partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
-
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
HLV pq,pqbar;
if(is_quark(partons[0])){
pq = to_HepLorentzVector(partons[0]);
pqbar = to_HepLorentzVector(partons[1]);
}
else{
pq = to_HepLorentzVector(partons[1]);
pqbar = to_HepLorentzVector(partons[0]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - pq - pqbar;
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc;
if (partons[0].type==-partons[1].type) { //leg b emits w
wc = true;}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_qqxb_current(
incoming[0].type, incoming[1].type,
pa, pb, pq, pqbar, pn, pl, plbar, wc
);
}
else{
current_factor = ME_W_qqxb_current(
incoming[0].type, incoming[1].type,
pa, pb, pq, pqbar, pn, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*9./8.*ladder_factor;
}
double MatrixElement::tree_kin_W_qqxf(
std::array<Particle, 2> const & incoming,
- std::vector<Particle> const & outgoing,
+ std::vector<Particle> const & partons,
HLV plbar, HLV pl,
bool WPlus, bool check_momenta
) const {
- const auto the_W = std::find_if(
- begin(outgoing), end(outgoing),
- [](Particle const & s){ return abs(s.type) == pid::Wp; }
- );
-
- const auto pW = to_HepLorentzVector(*the_W);
- std::vector<Particle> partons(begin(outgoing), the_W);
- partons.insert(end(partons), the_W + 1, end(outgoing));
- partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
-
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
HLV pq,pqbar;
if(is_quark(partons[partons.size() - 1])){
pq = to_HepLorentzVector(partons[partons.size() - 1]);
pqbar = to_HepLorentzVector(partons[partons.size() - 2]);
}
else{
pq = to_HepLorentzVector(partons[partons.size() - 2]);
pqbar = to_HepLorentzVector(partons[partons.size() - 1]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc;
if (incoming[0].type==partons[0].type) { //leg b emits w
wc = true;}
else{
wc = false;
q0 -=pl + plbar;
}
double current_factor;
if (WPlus){
current_factor = ME_W_qqxf_current(
incoming[0].type, incoming[1].type,
pa, pb, pq, pqbar, p1, pl, plbar, wc
);
}
else{
current_factor = ME_W_qqxf_current(
incoming[0].type, incoming[1].type,
pa, pb, pq, pqbar, p1, plbar, pl, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn
);
return current_factor*9./8.*ladder_factor;
}
double MatrixElement::tree_kin_W_qqxmid(
std::array<Particle, 2> const & incoming,
- std::vector<Particle> const & outgoing,
+ std::vector<Particle> const & partons,
HLV plbar, HLV pl,
bool WPlus, bool check_momenta
) const {
- const auto the_W = std::find_if(
- begin(outgoing), end(outgoing),
- [](Particle const & s){ return abs(s.type) == pid::Wp; }
- );
-
- const auto pW = to_HepLorentzVector(*the_W);
- std::vector<Particle> partons(begin(outgoing), the_W);
- partons.insert(end(partons), the_W + 1, end(outgoing));
- partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
-
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wc, wqq;
if (backmidquark->type == -(backmidquark+1)->type){ // Central qqx does not emit
wqq=false;
if (incoming[0].type==partons[0].type) {
wc = true;
}
else{
wc = false;
q0-=pl+plbar;
}
}
else{
wqq = true;
wc = false;
qqxt-=pl+plbar;
}
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
auto first_after_qqx = (backmidquark+2);
for(auto parton_it = begin_ladder; parton_it != first_after_qqx; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark-1);
int nbelow = std::distance(first_after_qqx, end_ladder);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
double current_factor;
if (WPlus){
current_factor = ME_W_qqxmid_current(
incoming[0].type, incoming[1].type, nabove, nbelow,
pa, pb, pq, pqbar, partonsHLV, pl, plbar, wqq, wc
);
}
else{
current_factor = ME_W_qqxmid_current(
incoming[0].type, incoming[1].type, nabove, nbelow,
pa, pb, pq, pqbar, partonsHLV, plbar, pl, wqq, wc
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, backmidquark-1,
q0, pa, pb, p1, pn
)*FKL_ladder_weight(
first_after_qqx, end_ladder,
qqxt, pa, pb, p1, pn
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double MatrixElement::tree_kin_Higgs(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
if(has_uno_gluon(incoming, outgoing)){
return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
}
if(outgoing.front().type == pid::Higgs){
return tree_kin_Higgs_first(incoming, outgoing, check_momenta);
}
if(outgoing.back().type == pid::Higgs){
return tree_kin_Higgs_last(incoming, outgoing, check_momenta);
}
return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
// TODO: code duplication with currents.cc
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == ParticleID::gluon) return K_g(pout, pin);
return C_F;
}
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == ParticleID::gluon)?C_A:C_F;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
ParticleID type2,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
// gluon case
#ifdef RHEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
incoming,
std::vector<Particle>(begin(outgoing) + 1, end(outgoing)),
check_momenta
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn
);
}
double MatrixElement::tree_kin_Higgs_last(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
assert(outgoing.back().type == pid::Higgs);
if(outgoing[outgoing.size()-2].type != pid::gluon) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(incoming, outgoing, check_momenta);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
incoming,
std::vector<Particle>(begin(outgoing), end(outgoing) - 1),
check_momenta
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
return MH2_forwardH(
pn, pb, incoming[0].type, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn
);
}
double MatrixElement::tree_kin_Higgs_between(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const {
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
std::vector<Particle> partons(begin(outgoing), the_Higgs);
partons.insert(end(partons), the_Higgs + 1, end(outgoing));
partons = tag_extremal_jet_partons(incoming, partons, check_momenta);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[has_unob_gluon(incoming, outgoing)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - (has_unof_gluon(incoming, outgoing)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
has_unob_gluon(incoming, outgoing)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
has_unof_gluon(incoming, outgoing)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor;
if(has_unob_gluon(incoming, outgoing)){
current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(has_unof_gluon(incoming, outgoing)){
current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double MatrixElement::tree_param_partons(
double alpha_s, double mur,
std::vector<Particle> const & partons
) const{
const double gs2 = 4.*M_PI*alpha_s;
double wt = std::pow(gs2, partons.size());
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(partons.size() >= 2);
for(size_t i = 1; i < partons.size()-1; ++i){
wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp());
}
}
return wt;
}
double MatrixElement::tree_param(
double mur,
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing
) const{
const double alpha_s = alpha_s_(mur);
auto AWZH_boson = std::find_if(
begin(outgoing), end(outgoing),
[](auto const & p){return is_AWZH_boson(p);}
);
double AWZH_coupling = 1.;
if(AWZH_boson != end(outgoing)){
switch(AWZH_boson->type){
case pid::Higgs: {
AWZH_coupling = alpha_s*alpha_s;
break;
}
// TODO
case pid::Wp:{
AWZH_coupling = alpha_w*alpha_w/2;
break;
}
case pid::Wm:{
AWZH_coupling = alpha_w*alpha_w/2;
break;
}
case pid::photon:
case pid::Z:
default:
throw std::logic_error("Emission of boson of unsupported type");
}
}
if(has_unob_gluon(incoming, outgoing)){
return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(outgoing) + 1, end(outgoing)})
);
}
if(has_unof_gluon(incoming, outgoing)){
return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(outgoing), end(outgoing) - 1})
);
}
return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(outgoing));
}
double MatrixElement::tree(
double mur,
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
std::unordered_map<int, std::vector<Particle>> const & decays,
bool check_momenta
) const {
return tree_param(mur, incoming, outgoing)*tree_kin(
incoming, outgoing, decays, check_momenta
);
}
} // namespace RHEJ