diff --git a/include/RHEJ/MatrixElement.hh b/include/RHEJ/MatrixElement.hh
index 012fc6f..bd8b92b 100644
--- a/include/RHEJ/MatrixElement.hh
+++ b/include/RHEJ/MatrixElement.hh
@@ -1,240 +1,254 @@
/** \file
* \brief Contains the MatrixElement Class
*/
#pragma once
#include <functional>
#include "RHEJ/config.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/HiggsCouplingSettings.hh"
#include "CLHEP/Vector/LorentzVector.h"
namespace RHEJ{
//! Class to calculate the squares of matrix elements
class MatrixElement{
public:
/** \brief MatrixElement Constructor
* @param alpha_s Function taking the renormalisation scale
* and returning the strong coupling constant
* @param conf General matrix element settings
*/
MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
);
/**
* \brief regulated HEJ matrix element
* @param mur Value of the renormalisation scale
* @param incoming Incoming particles
* @param outgoing Outgoing particles
* @param check_momenta Special treatment for partons inside extremal jets
* @returns The HEJ matrix element including virtual corrections
*
* cf. eq. (22) in \cite Andersen:2011hs
* Incoming particles should be ordered by ascending z momentum.
* Outgoing particles should be ordered by ascending rapidity.
*
* \internal Relation to standard HEJ Met2: MatrixElement = Met2*shat^2/(pdfta*pdftb)
*/
double 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;
//! HEJ tree-level matrix element
/**
* @param mur Value of the renormalisation scale
* @param incoming Incoming particles
* @param outgoing Outgoing particles
* @param check_momenta Special treatment for partons inside extremal jets
* @returns The HEJ matrix element without virtual corrections
*
* cf. eq. (22) in \cite Andersen:2011hs
* Incoming particles should be ordered by ascending z momentum.
* Outgoing particles should be ordered by ascending rapidity.
*/
double 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;
//! HEJ tree-level matrix element - parametric part
/**
* @param mur Value of the renormalisation scale
* @param incoming Incoming particles
* @param outgoing Outgoing particles
* @returns The parametric part of the tree matrix element
*
* cf. eq. (22) in \cite Andersen:2011hs
*
* The tree level matrix element factorises into a parametric part
* which depends on the theory parameters (alpha_s and scale)
* and a kinematic part comprising the dependence on the particle momenta
* and colour factors. This function returns the former.
*/
double tree_param(
double mur,
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing
) const;
//! HEJ tree-level matrix element - kinematic part
/**
* @param incoming Incoming particles
* @param outgoing Outgoing particles
* @param check_momenta Special treatment for partons inside extremal jets
* @returns The kinematic part of the tree matrix element
*
* cf. eq. (22) in \cite Andersen:2011hs
* Incoming particles should be ordered by ascending z momentum.
* Outgoing particles should be ordered by ascending rapidity.
*
* The tree level matrix element factorises into a parametric part
* which depends on the theory parameters (alpha_s and scale)
* and a kinematic part comprising the dependence on the particle momenta
* and colour factors. This function returns the latter.
*/
double 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;
/**
* \brief Calculates the Virtual Corrections
* @param mur Value of the renormalisation scale
* @param in Incoming particles
* @param out Outgoing particles
* @returns The Virtual Corrections of the Matrix Element
*
* Incoming particles should be ordered by ascending z momentum.
* Outgoing particles should be ordered by ascending rapidity.
*
* The all order virtual corrections to LL in the MRK limit is
* given by replacing 1/t in the scattering amplitude according to the
* lipatov ansatz.
*
* cf. second-to-last line of eq. (22) in \cite Andersen:2011hs
* note that indices are off by one, i.e. out[0].p corresponds to p_1
*/
double virtual_corrections(
double mur,
std::array<Particle, 2> const & in,
std::vector<Particle> const & out
) const;
private:
//! \internal cf. last line of eq. (22) in \cite Andersen:2011hs
double omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j, double lambda
) const;
double tree_kin_jets(
std::array<Particle, 2> const & incoming,
std::vector<Particle> partons,
bool check_momenta
) const;
double 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;
double tree_kin_W_FKL(
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;
double tree_kin_W_unob(
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;
double tree_kin_W_unof(
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;
+ double tree_kin_W_qqxb(
+ 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;
+ double tree_kin_W_qqxf(
+ 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;
double tree_kin_Higgs(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const;
double tree_kin_Higgs_first(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const;
double tree_kin_Higgs_last(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const;
/**
* \internal
* \brief Higgs inbetween extremal partons.
*
* Note that in the case of unordered emission, the Higgs is *always*
* treated as if in between the extremal (FKL) partons, even if its
* rapidity is outside the extremal parton rapidities
*/
double tree_kin_Higgs_between(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing,
bool check_momenta
) const;
double tree_param_partons(
double alpha_s, double mur,
std::vector<Particle> const & partons
) const;
std::vector<int> in_extremal_jet_indices(
std::vector<fastjet::PseudoJet> const & partons
) const;
std::vector<Particle> tag_extremal_jet_partons(
std::array<Particle, 2> const & incoming,
std::vector<Particle> out_partons, bool check_momenta
) const;
double MH2_forwardH(
CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
pid::ParticleID type2,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const;
std::function<double (double)> alpha_s_;
MatrixElementConfig param_;
};
}
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 716b3bf..593de0a 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1268 +1,1542 @@
#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 ) {//b gluon => W emission off a 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;}
+ }
+ }
+ }
+
+ /* \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) {// a gluon => W emission off b 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;}
+ }
+ }
+
/** \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 {
if(has_unob_gluon(incoming, outgoing)){
throw std::logic_error("unordered emission not yet implemented for W+jets");
//return tree_kin_W_unob(incoming, outgoing, check_momenta);
}
else if(has_unof_gluon(incoming, outgoing)){
throw std::logic_error("unordered emission not yet implemented for W+jets");
// return tree_kin_W_unof(incoming, outgoing, check_momenta);
}
else if(has_Ex_qqx(incoming, outgoing)){
throw std::logic_error("Extremal qqx not yet implemented for W+jets");
// return tree_kin_W_Exqqx(incoming, outgoing, check_momenta);
}
else if(has_mid_qqx(outgoing)){
throw std::logic_error("Central qqx not yet implemented for W+jets");
// return tree_kin_W_qqxCentral(incoming, outgoing, check_momenta);
}
else{
return tree_kin_W_FKL(incoming, outgoing, decays, WPlus, check_momenta);
}
}
double MatrixElement::tree_kin_W_FKL(
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 {
const auto the_W = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
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 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::unordered_map<int, std::vector<Particle>> const & decays,
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; }
);
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 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::unordered_map<int, std::vector<Particle>> const & decays,
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; }
);
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 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::unordered_map<int, std::vector<Particle>> const & decays,
+ 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; }
+ );
+
+ 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 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::unordered_map<int, std::vector<Particle>> const & decays,
+ 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; }
+ );
+
+ 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 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_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