Page MenuHomeHEPForge
Files F7877582

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/currents.hh b/include/HEJ/currents.hh
index 6d1802d..8424508 100644
--- a/include/HEJ/currents.hh
+++ b/include/HEJ/currents.hh
@@ -1,1208 +1,1208 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions.
*
* This file contains all the necessary functions to compute the current
* contractions for all valid HEJ processes. PJETS, H+JETS and W+JETS along
* with some unordered counterparts.
*
* @TODO add a namespace
*/
#pragma once
#include <complex>
#include <vector>
#include <ostream>
#include <CLHEP/Vector/LorentzVector.h>
typedef std::complex<double> COM;
typedef COM current[4];
typedef CLHEP::HepLorentzVector HLV;
//! Square of qQ->qenuQ W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qenuQ Scattering
*
* This returns the square of the current contractions in qQ->qenuQ scattering
* with an emission of a W Boson.
*/
double ME_W_qQ (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQ->qbarenuQ W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qbarQ->qbarenuQ Scattering
*
* This returns the square of the current contractions in qbarQ->qbarenuQ
* scattering with an emission of a W Boson.
*/
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qenuQbar W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qenuQbar Scattering
*
* This returns the square of the current contractions in qQbar->qenuQbar
* scattering with an emission of a W Boson.
*/
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarenuQbar W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarenuQbar Scattering
*
* This returns the square of the current contractions in qbarQbar->qbarenuQbar
* scattering with an emission of a W Boson.
*/
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qenug W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qg->qenug Scattering
*
* This returns the square of the current contractions in qg->qenug scattering
* with an emission of a W Boson.
*/
double ME_W_qg (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarenug W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qbarg->qbarenug Scattering
*
* This returns the square of the current contractions in qbarg->qbarenug
* scattering with an emission of a W Boson.
*/
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in);
//! qQg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qQ->qQg Scattering
*
* This returns the square of the current contractions in qQg->qQg scattering
* with an emission of a W Boson.
*/
double ME_W_unob_qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! qbarQg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qbarQ->qbarQg Scattering
*
* This returns the square of the current contractions in qbarQg->qbarQg
* scattering with an emission of a W Boson.
*/
double ME_W_unob_qbarQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! qQbarg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qQbar->qQbarg Scattering
*
* This returns the square of the current contractions in qQbarg->qQbarg
* scattering with an emission of a W Boson.
*/
double ME_W_unob_qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! qbarQbarg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qbarQbar->qbarQbarg Scattering
*
* This returns the square of the current contractions in qbarQbarg->qbarQbarg
* scattering with an emission of a W Boson.
*/
double ME_W_unob_qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! Wjets Unordered forwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qQ->gqQ Scattering
*
/ * This returns the square of the current contractions in qQg->gqQ scattering
* with an emission of a W Boson.
*/
double ME_W_unof_qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! Wjets Unordered forwards opposite leg to W
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qbarQ->gqbarQ Scattering
*
* This returns the square of the current contractions in qbarQg->gqbarQ
* scattering with an emission of a W Boson.
*/
double ME_W_unof_qbarQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! Wjets Unordered forwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qQbar->gqQbar Scattering
*
* This returns the square of the current contractions in qQbarg->gqQbar
* scattering with an emission of a W Boson.
*/
double ME_W_unof_qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! Wjets Unordered forwards opposite leg to W
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qbarQbar->gqbarQbar Scattering
*
* This returns the square of the current contractions in qbarQbarg->gqbarQbar
* scattering with an emission of a W Boson.
*/
double ME_W_unof_qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! W+uno same leg
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qQ->qQg Scattering
*
* This returns the square of the current contractions in gqQ->gqQ scattering
* with an emission of a W Boson.
*/
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! W+uno same leg. quark anti-quark
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qQbar->gqQbar Scattering
*
* This returns the square of the current contractions in gqQbar->gqQbar
* scattering with an emission of a W Boson. (Unordered Same Leg)
*/
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! W+uno same leg. quark gluon
/**
* @param p1out Momentum of final state quark a
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state gluon b
* @param p2in Momentum of intial state gluon b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qg->gqg Scattering
*
* This returns the square of the current contractions in qg->gqg scattering
* with an emission of a W Boson.
*/
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! W+uno same leg. anti-quark quark
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qbarQ->gqbarQ Scattering
*
* This returns the square of the current contractions in qbarQ->gqbarQ
* scattering with an emission of a W Boson.
*/
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! W+uno same leg. anti-quark anti-quark
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for qbarQbar->gqbarQbar Scattering
*
* This returns the square of the current contractions in gqbarQbar->qbarQbar
* scattering with an emission of a W Boson.
*/
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! W+uno same leg. anti-quark gluon
/**
* @param p1out Momentum of final state anti-quark a
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state gluon b
* @param p2in Momentum of intial state gluon b
* @param pg Momentum of final state unordered gluon
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @returns Square of the current contractions for ->gqbarg Scattering
*
* This returns the square of the current contractions in qbarg->gqbarg
* scattering with an emission of a W Boson.
*/
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pg,
HLV plbar, HLV pl);
//! W+Extremal qqx. qxqQ
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state anti-quark b
* @param p2in Momentum of final state gluon b
* @returns Square of the current contractions for ->qbarqQ Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout,HLV plbar,HLV pl, HLV pqbarout, HLV p2out, HLV p2in);
//! W+Extremal qqx. qqxQ
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state anti-quark b
* @param p2in Momentum of final state gluon b
* @returns Square of the current contractions for ->qqbarQ Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqout,HLV plbar,HLV pl, HLV pqbarout, HLV p2out, HLV p2in);
//! W+Extremal qqx. gg->qxqg
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state gluon b
* @param p2in Momentum of final state gluon b
* @returns Square of the current contractions for gg->qbarqg Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout,HLV plbar,HLV pl, HLV pqbarout, HLV p2out, HLV p2in);
//! W+Extremal qqx. gg->qqxg
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state gluon a
* @param p2in Momentum of final state gluon b
* @returns Square of the current contractions for gg->qqbarg Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout,HLV plbar,HLV pl, HLV pqout, HLV p2out, HLV p2in);
//! W+Extremal qqx. gg->qqxg. qqx on forwards leg, W emission backwards leg.
/**
* @param pa Momentum of initial state (anti-)quark
* @param pb Momentum of initial state gluon
* @param p1 Momentum of final state (anti-)quark (after W emission)
* @param p2 Momentum of final state anti-quark
* @param p3 Momentum of final state quark
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param aqlinepa Is opposite extremal leg to qqx a quark or antiquark line
* @returns Square of the current contractions for gq->qqbarqW Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated via current contraction of existing currents.
* Assumes qqx split from forwards leg, W emission from backwards leg.
* Switch input (pa<->pb, p1<->pn) if calculating forwards qqx.
*/
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2, HLV p3,HLV plbar,HLV pl, bool aqlinepa);
//! W+Jets qqxCentral. qqx W emission.
/**
* @param pa Momentum of initial state particle a
* @param pb Momentum of initial state particle b
* @param pl Momentum of final state lepton
* @param plbar Momentum of final state anti-lepton
* @param partons Vector of outgoing parton momenta
* @param aqlinepa True= pa is anti-quark
* @param aqlinepb True= pb is anti-quark
* @param qqxmarker Ordering of the qqbar pair produced (qqx vs qxq)
* @param nabove Number of lipatov vertices "above" qqbar pair
* @param nbelow Number of lipatov vertices "below" qqbar pair
* @returns Square of the current contractions for qq>qQQbarWq Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_WCenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove);
//emission from backwards leg
//! W+Jets qqxCentral. W emission from backwards leg.
/**
* @param ka Momentum of initial state particle a
* @param kb Momentum of initial state particle b
* @param pl Momentum of final state lepton
* @param plbar Momentum of final state anti-lepton
* @param partons outgoing parton momenta
* @param aqlinepa True= pa is anti-quark
* @param aqlinepb True= pb is anti-quark
* @param qqxmarker Ordering of the qqbar pair produced (qqx vs qxq)
* @param nabove Number of lipatov vertices "above" qqbar pair
* @param nbelow Number of lipatov vertices "below" qqbar pair
* @param forwards Swap to emission off front leg TODO:remove so args can be const
* @returns Square of the current contractions for qq>qQQbarWq Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double ME_W_Cenqqx_qq(HLV ka, HLV kb,HLV pl,HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
int nbelow, bool forwards); //Doing
//! Square of qQ->qQ Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qQ Scattering
*/
-double jM2qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+double ME_qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qQbar Scattering
*
* @note this can be used for qbarQ->qbarQ Scattering by inputting arguments
* appropriately.
*/
-double jM2qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+double ME_qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*/
-double jM2qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+double ME_qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qg->qg Scattering
*
* @note this can be used for gq->gq Scattering by inputting arguments
* appropriately.
*/
-double jM2qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+double ME_qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* @note this can be used for gqbar->gqbar Scattering by inputting arguments
* appropriately.
*/
-double jM2qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+double ME_qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of gg->gg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for gg->gg Scattering
*/
-double jM2gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
+double ME_gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of gg->gg Higgs+Jets Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @param q1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for gg->gg Scattering
*
* g~p1 g~p2
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if g is backward, q1 is forward)
*/
-double MH2gg (HLV p1out, HLV p1in,
+double ME_H_gg (HLV p1out, HLV p1in,
HLV p2out, HLV p2in,
HLV q1, HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of gq->gq Higgs+Jets Scattering Current with Higgs before Gluon
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @param pH Momentum of Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contraction
*/
-double MH2gq_outsideH(HLV p1out, HLV p1in,
+double ME_Houtside_gq(HLV p1out, HLV p1in,
HLV p2out, HLV p2in,
HLV pH,
double mt,
bool include_bottom, double mb);
//! Square of qg->qg Higgs+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @param q1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qg->qg Scattering
*
* q~p1 g~p2 (i.e. ALWAYS p1 for quark, p2 for gluon)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if g is backward, q1 is forward)
*/
-double MH2qg (HLV p1out, HLV p1in,
+double ME_H_qg (HLV p1out, HLV p1in,
HLV p2out, HLV p2in,
HLV q1, HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qbarg->qbarg Higgs+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @param q1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* qbar~p1 g~p2 (i.e. ALWAYS p1 for anti-quark, p2 for gluon)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if g is backward, q1 is forward)
*/
-double MH2qbarg (HLV p1out, HLV p1in,
+double ME_H_qbarg (HLV p1out, HLV p1in,
HLV p2out, HLV p2in,
HLV q1, HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qQ->qQ Higgs+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param q1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qQ->qQ Scattering
*
* q~p1 Q~p2 (i.e. ALWAYS p1 for quark, p2 for quark)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if Q is backward, q1 is forward)
*/
-double MH2qQ (HLV p1out, HLV p1in,
+double ME_H_qQ (HLV p1out, HLV p1in,
HLV p2out, HLV p2in,
HLV q1, HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qQbar->qQbar Higgs+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param q1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qQ->qQ Scattering
*
* q~p1 Qbar~p2 (i.e. ALWAYS p1 for quark, p2 for anti-quark)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if Qbar is backward, q1 is forward)
*/
-double MH2qQbar (HLV p1out, HLV p1in,
+double ME_H_qQbar (HLV p1out, HLV p1in,
HLV p2out, HLV p2in,
HLV q1, HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qbarQ->qbarQ Higgs+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param q1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qbarQ->qbarQ Scattering
*
* qbar~p1 Q~p2 (i.e. ALWAYS p1 for anti-quark, p2 for quark)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if Q is backward, q1 is forward)
*/
-double MH2qbarQ (HLV p1out, HLV p1in,
+double ME_H_qbarQ (HLV p1out, HLV p1in,
HLV p2out, HLV p2in,
HLV q1, HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qbarQbar->qbarQbar Higgs+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param q1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*
* qbar~p1 Qbar~p2 (i.e. ALWAYS p1 for anti-quark, p2 for anti-quark)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if Qbar is backward, q1 is forward)
*/
-double MH2qbarQbar (HLV p1out, HLV p1in,
+double ME_H_qbarQbar (HLV p1out, HLV p1in,
HLV p2out, HLV p2in,
HLV q1, HLV qH2,
double mt,
bool include_bottom, double mb);
//! @name Unordered forward
//! @{
//! Square of qQ->gqQ Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered gluon
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qQ->gqQ Scattering
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
-double jM2unogqHQ (HLV pg, HLV p1out,
+double ME_H_unof_qQ (HLV pg, HLV p1out,
HLV p1in, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qQbar->gqQbar Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered gluon
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qQbar->gqQbar Scattering
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
-double jM2unogqHQbar (HLV pg, HLV p1out,
+double ME_H_unof_qQbar (HLV pg, HLV p1out,
HLV p1in, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qbarQ->gqbarQ Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered gluon
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qbarQ->gqbarQ Scattering
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
-double jM2unogqbarHQ (HLV pg, HLV p1out,
+double ME_H_unof_qbarQ (HLV pg, HLV p1out,
HLV p1in, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qbarQbar->gqbarQbar Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered gluon
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qbarQbar->gqbarQbar Scattering
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
-double jM2unogqbarHQbar (HLV pg, HLV p1out,
+double ME_H_unof_qbarQbar (HLV pg, HLV p1out,
HLV p1in, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qg->gqg Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered gluon
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qg->gqg Scattering
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
-double jM2unogqHg (HLV pg, HLV p1out,
+double ME_H_unof_qg (HLV pg, HLV p1out,
HLV p1in, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qbarg->gqbarg Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered gluon
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qbarg->gbarg Scattering
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
-double jM2unogqbarHg (HLV pg, HLV p1out,
+double ME_H_unof_qbarg (HLV pg, HLV p1out,
HLV p1in, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! @}
//! @name Unordered backwards
//! @{
//! Square of qbarQ->qbarQg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered b gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qbarQ->qbarQg Scattering
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
-double jM2unobqbarHQg (HLV p1out, HLV p1in,
+double ME_H_unob_qbarQ (HLV p1out, HLV p1in,
HLV pg, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qQ->qQg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered b gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qQ->qQg Scattering
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
-double jM2unobqHQg (HLV p1out, HLV p1in,
+double ME_H_unob_qQ (HLV p1out, HLV p1in,
HLV pg, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qQbar->qQbarg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered b gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qQbar->qQbarg Scattering
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
-double jM2unobqHQbarg (HLV p1out, HLV p1in,
+double ME_H_unob_qQbar (HLV p1out, HLV p1in,
HLV pg, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of qbarQbar->qbarQbarg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered b gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for qbarQbar->qbarQbarg Scattering
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
-double jM2unobqbarHQbarg (HLV p1out, HLV p1in,
+double ME_H_unob_qbarQbar (HLV p1out, HLV p1in,
HLV pg, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of gQbar->gQbarg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered b gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for gQbar->gQbarg Scattering
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
-double jM2unobgHQbarg (HLV p1out, HLV p1in,
+double ME_H_unob_gQbar (HLV p1out, HLV p1in,
HLV pg, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! Square of gQ->gQg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered b gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @param qH1 Momentum of t-channel propagator before Higgs
* @param qH2 Momentum of t-channel propagator after Higgs
* @param mt Top quark mass
* @param include_bottom Specifies whether bottom corrections are included
* @param mb Bottom quark mass
* @returns Square of the current contractions for gQ->gQg Scattering
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
-double jM2unobgHQg (HLV p1out, HLV p1in,
+double ME_H_unob_gQ (HLV p1out, HLV p1in,
HLV pg, HLV p2out,
HLV p2in, HLV qH1,
HLV qH2,
double mt,
bool include_bottom, double mb);
//! @}
//! @name impact factors for Higgs + jet
//! @{
//! Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
/**
* @param p2 Momentum of Particle 2
* @param p1 Momentum of Particle 1
* @param pH Momentum of Higgs
* @returns Value of Eq. (4.22) in Hep-ph/0301013 with modifications
*
* This gives the impact factor. First it determines first whether this is the
* case p1p\sim php>>p3p or the opposite
*/
double C2gHgm(HLV p2, HLV p1, HLV pH);
//! Implements Eq. (4.23) in hep-ph/0301013 with modifications to incoming plus momenta
/**
* @param p2 Momentum of Particle 2
* @param p1 Momentum of Particle 1
* @param pH Momentum of Higgs
* @returns Value of Eq. (4.23) in Hep-ph/0301013
*
* This gives the impact factor. First it determines first whether this is the
* case p1p\sim php>>p3p or the opposite
*/
double C2gHgp(HLV p2, HLV p1, HLV pH);
//! Implements Eq. (4.22) in hep-ph/0301013
/**
* @param p2 Momentum of Particle 2
* @param p1 Momentum of Particle 1
* @param pH Momentum of Higgs
* @returns Value of Eq. (4.22) in Hep-ph/0301013
*
* This gives the impact factor. First it determines first whether this is the
* case p1p\sim php>>p3p or the opposite
*/
double C2qHqm(HLV p2, HLV p1, HLV pH);
//! @}
/** \class CCurrent currents.hh "include/HEJ/currents.hh"
* \brief This is the a new class structure for currents.
*/
class CCurrent
{
public:
CCurrent(COM sc0, COM sc1, COM sc2, COM sc3)
:c0(sc0),c1(sc1),c2(sc2),c3(sc3)
{};
CCurrent(const HLV p)
{
c0=p.e();
c1=p.px();
c2=p.py();
c3=p.pz();
};
CCurrent()
{};
CCurrent operator+(const CCurrent& other);
CCurrent operator-(const CCurrent& other);
CCurrent operator*(const double x);
CCurrent operator*(const COM x);
CCurrent operator/(const double x);
CCurrent operator/(const COM x);
friend std::ostream& operator<<(std::ostream& os, const CCurrent& cur);
COM dot(HLV p1);
COM dot(CCurrent p1);
COM c0,c1,c2,c3;
};
/* std::ostream& operator <<(std::ostream& os, const CCurrent& cur); */
CCurrent operator * ( double x, CCurrent& m);
CCurrent operator * ( COM x, CCurrent& m);
CCurrent operator / ( double x, CCurrent& m);
CCurrent operator / ( COM x, CCurrent& m);
//! Current <incoming state | mu | outgoing state>
/**
* These functions are a mess. There are many more defined in the source file
* than declared in the header - and the arguments are mislabelled in some
* cases. Need to investigate.
*/
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur);
//! Current <outgoing state | mu | outgoing state>
/**
* These functions are a mess. There are many more defined in the source file
* than declared in the header - and the arguments are mislabelled in some
* cases. Need to investigate.
*/
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* These functions are a mess. There are many more defined in the source file
* than declared in the header - and the arguments are mislabelled in some
* cases. Need to investigate.
*/
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* These functions are a mess. There are many more defined in the source file
* than declared in the header - and the arguments are mislabelled in some
* cases. Need to investigate.
*/
CCurrent joi (HLV pout, bool helout, HLV pin, bool helin);
//! Current <incoming state | mu | outgoing state>
/**
* These functions are a mess. There are many more defined in the source file
* than declared in the header - and the arguments are mislabelled in some
* cases. Need to investigate.
*/
CCurrent jio (HLV pout, bool helout, HLV pin, bool helin);
//! Current <outgoing state | mu | outgoing state>
/**
* These functions are a mess. There are many more defined in the source file
* than declared in the header - and the arguments are mislabelled in some
* cases. Need to investigate.
*/
CCurrent joo (HLV pout, bool helout, HLV pin, bool helin);
inline COM cdot(const current & j1, const current & j2)
{
return j1[0]*j2[0]-j1[1]*j2[1]-j1[2]*j2[2]-j1[3]*j2[3];
}
inline COM cdot(const HLV & p, const current & j1) {
return j1[0]*p.e()-j1[1]*p.x()-j1[2]*p.y()-j1[3]*p.z();
}
inline void cmult(const COM & factor, const current & j1, current &cur)
{
cur[0]=factor*j1[0];
cur[1]=factor*j1[1];
cur[2]=factor*j1[2];
cur[3]=factor*j1[3];
}
// WHY!?!
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, const current & j5, current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0]+j4[0]+j5[0];
sum[1]=j1[1]+j2[1]+j3[1]+j4[1]+j5[1];
sum[2]=j1[2]+j2[2]+j3[2]+j4[2]+j5[2];
sum[3]=j1[3]+j2[3]+j3[3]+j4[3]+j5[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, current &sum) {
sum[0] = j1[0] + j2[0] + j3[0] + j4[0];
sum[1] = j1[1] + j2[1] + j3[1] + j4[1];
sum[2] = j1[2] + j2[2] + j3[2] + j4[2];
sum[3] = j1[3] + j2[3] + j3[3] + j4[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0];
sum[1]=j1[1]+j2[1]+j3[1];
sum[2]=j1[2]+j2[2]+j3[2];
sum[3]=j1[3]+j2[3]+j3[3];
}
inline void cadd(const current & j1, const current & j2, current &sum)
{
sum[0]=j1[0]+j2[0];
sum[1]=j1[1]+j2[1];
sum[2]=j1[2]+j2[2];
sum[3]=j1[3]+j2[3];
}
inline double abs2(const COM & a)
{
return (a*conj(a)).real();
}
inline double vabs2(const CCurrent & cur)
{
return abs2(cur.c0)-abs2(cur.c1)-abs2(cur.c2)-abs2(cur.c3);
}
inline double vre(const CCurrent & a, const CCurrent & b)
{
return real(a.c0*conj(b.c0)-a.c1*conj(b.c1)-a.c2*conj(b.c2)-a.c3*conj(b.c3));
}
// @TODO: These are not currents and should be moved elsewhere.
double K_g(double p1minus, double paminus);
double K_g(HLV const & pout, HLV const & pin);
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 9519b5f..5d35ee0 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1755 +1,1755 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/currents.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/utility.hh"
namespace HEJ{
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
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;
}
Weights MatrixElement::operator()(
Event const & event
) const {
return tree(event)*virtual_corrections(event);
}
Weights MatrixElement::tree(
Event const & event
) const {
return tree_param(event)*tree_kin(event);
}
Weights MatrixElement::tree_param(
Event const & event
) const {
if(! is_HEJ(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
Weights MatrixElement::virtual_corrections(
Event const & event
) const {
if(! is_HEJ(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = virtual_corrections(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = virtual_corrections(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
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(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
bool wc = true;
bool wqq = false;
// With extremal qqx or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::FKL) {
if (in[0].type != partons[0].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::qqxexb) {
q -= partons[1].p;
++first_idx;
if (abs(partons[0].type) != abs(partons[1].type)){
q -= WBoson.p;
wc = false;
}
}
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(abs(backquark->type) != abs((backquark+1)->type)) {
wqq=true;
wc=false;
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
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)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqx) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
if (wc) q -= WBoson.p;
assert(
nearby(q, -1*pb, norm)
|| is_AWZH_boson(partons.back().type)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
double MatrixElement::virtual_corrections(
Event const & event,
double mur
) const{
auto const & in = event.incoming();
auto const & out = event.outgoing();
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
const auto AWZH_boson = std::find_if(
begin(out), end(out),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson != end(out) && abs(AWZH_boson->type) == pid::Wp){
return virtual_corrections_W(event, mur, *AWZH_boson);
}
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 boson, extremal qqx 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)
|| event.type() == event_type::unob
|| event.type() == event_type::qqxexb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return 0;
last_idx = std::distance(begin(out), backquark);
first_idx_qqx = last_idx+1;
}
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)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqx) q -= out[last_idx+1].p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
} // namespace HEJ
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));
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 lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
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 lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
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);
+ return ME_gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
- return jM2qg(pn,pb,p1,pa);
+ return ME_qg(pn,pb,p1,pa);
else
- return jM2qbarg(pn,pb,p1,pa);
+ return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
- return jM2qg(p1,pa,pn,pb);
+ return ME_qg(p1,pa,pn,pb);
else
- return jM2qbarg(p1,pa,pn,pb);
+ return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
- return jM2qQ(pn,pb,p1,pa);
+ return ME_qQ(pn,pb,p1,pa);
else
- return jM2qQbar(pn,pb,p1,pa);
+ return ME_qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
- return jM2qQbar(p1,pa,pn,pb);
+ return ME_qQbar(p1,pa,pn,pb);
else
- return jM2qbarQbar(pn,pb,p1,pa);
+ return ME_qbarQbar(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
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @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.
assert(!(aptype==21 && bptype==21));
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_W_qg(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_W_qg(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qbarg(p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (bptype>0) {
if (aptype>0)
return ME_W_qQ(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa);
}
else {
if (aptype>0)
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa);
else
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return ME_W_qQ(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0)
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb);
else
return ME_W_qbarQbar(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
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @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 ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else{
if (aptype>0)
return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else
return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl);
else //qqbar
return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else //qbarqbar
return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
}
}
/** 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
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @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 ME_Wuno_qg(pn, pb, p1, pa, pg, plbar, pl);
else
return ME_Wuno_qbarg(pn, pb, p1, pa, pg, plbar, pl);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return ME_Wuno_qQ(pn,pb,p1,pa,pg,plbar,pl);
else
return ME_Wuno_qQbar(pn,pb,p1,pa,pg,plbar,pl);
}
else{
if (aptype>0)
return ME_Wuno_qbarQ(pn,pb,p1,pa,pg,plbar,pl);
else
return ME_Wuno_qbarQbar(pn,pb,p1,pa,pg,plbar,pl);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return ME_W_unof_qQ(p1,pa,pn,pb,pg,plbar,pl);
// return ME_W_unof_qQ(pn,pb,p1,pa,pg,plbar,pl);
else //qqbar
return ME_W_unof_qQbar(p1,pa,pn,pb,pg,plbar,pl);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return ME_W_unof_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
else //qbarqbar
return ME_W_unof_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
}
}
}
}
/** \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
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @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 ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn,pb)*CFbackward;}
else {
return ME_WExqqx_qbarqg(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 ME_WExqqx_qqbarQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
else{
return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
}
else { // W Must be emitted from forwards leg.
if(bptype > 0){
if (swapQuarkAntiquark){
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;}
else{
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;}
} else {
if (swapQuarkAntiquark){
return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;}
else{
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, true)*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
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @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 ME_WExqqx_qqbarg(pb, pqbar, plbar, pl, pq, p1,pa)*CFforward;}
else {
return ME_WExqqx_qbarqg(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 ME_WExqqx_qbarqQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
else {
return ME_WExqqx_qqbarQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
}
// W Must be emitted from backwards leg.
if (aptype > 0){
if (swapQuarkAntiquark){
return ME_W_Exqqx_QQq(pa,pb, p1, pqbar, pq, plbar, pl, false)*CFforward;}
else{
return ME_W_Exqqx_QQq(pa,pb, p1, pq, pqbar, plbar, pl, false)*CFforward;}
} else
{
if (swapQuarkAntiquark){
return ME_W_Exqqx_QQq(pa,pb, p1, pqbar, pq, plbar, pl, true)*CFforward;}
else{
return ME_W_Exqqx_QQq(pa,pb, p1, pq, pqbar, plbar, pl, true)*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 wt=1.;
if (aptype==21) wt*=K_g(partons.front(),pa)/HEJ::C_F;
if (bptype==21) wt*=K_g(partons.back(),pb)/HEJ::C_F;
if (aptype <=0 && bptype <=0){ // Both External AntiQuark
if (wqq==1){//emission from central qqbar
return wt*ME_WCenqqx_qq(pa, pb, pl,plbar, partons,true,true,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true,true,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(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*ME_WCenqqx_qq(pa, pb, pl,plbar, partons, false, true,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons,false,true,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(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*ME_WCenqqx_qq(pa, pb, pl,plbar, partons, true, false,
swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, true, false,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(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*ME_WCenqqx_qq(pa, pb, pl, plbar, partons, false, false,
swapQuarkAntiquark, nabove);}
else if (wc==1){//emission from b leg
return wt*ME_W_Cenqqx_qq(pa, pb, pl,plbar, partons, false, false,
swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*ME_W_Cenqqx_qq(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);
+ return ME_H_gg(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.;
+ return ME_H_qg(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.;
+ return ME_H_qbarg(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.;
+ return ME_H_qg(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.;
+ return ME_H_qbarg(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.);
+ return ME_H_qQ(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.);
+ return ME_H_qQbar(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.);
+ return ME_H_qQbar(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.);
+ return ME_H_qbarQbar(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);
+ return ME_H_unof_qg(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);
+ return ME_H_unof_qbarg(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);
+ return ME_H_unof_qQ(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);
+ return ME_H_unof_qQbar(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);
+ return ME_H_unof_qbarQ(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);
+ return ME_H_unof_qbarQbar(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);
+ return ME_H_unob_gQ(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);
+ return ME_H_unob_gQbar(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);
+ return ME_H_unob_qQ(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);
+ return ME_H_unob_qbarQ(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);
+ return ME_H_unob_qQbar(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);
+ return ME_H_unob_qbarQbar(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(HEJ::MatrixElementConfig const & config) {
#ifndef HEJ_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 HEJ{
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
double MatrixElement::tree_kin(
Event const & ev
) const {
if(! is_HEJ(ev.type())) return 0.;
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return tree_kin_jets(ev);
switch(AWZH_boson->type){
case pid::Higgs:
return tree_kin_Higgs(ev);
case pid::Wp:
case pid::Wm:
return tree_kin_W(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("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 lambda
){
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, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
// TODO: avoid reclustering
fastjet::ClusterSequence cs(to_PseudoJet(out_partons), ev.jet_def());
const auto jets = sorted_by_rapidity(cs.inclusive_jets(ev.min_jet_pt()));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission or qqx
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
assert(jets.size() >= 3);
++most_backward;
}
else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
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(HEJ::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(
Event const & ev
) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if(is_uno(ev.type())){
throw not_implemented("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,
param_.regulator_lambda
);
}
namespace{
double tree_kin_W_FKL(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
double tree_kin_W_unob(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto pg = to_HepLorentzVector(partons[0]);
auto p1 = to_HepLorentzVector(partons[1]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - p1 -pg;
if (aptype!=partons[1].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_unob_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_unof(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 2]);
auto pg = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_unof_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxb(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
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 begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - pq - pqbar;
if (partons[1].type!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_qqxb_current(
aptype, bptype, pa, pb,
pq, pqbar, pn, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxf(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
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 begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_qqxf_current(
aptype, bptype, pa, pb,
pq, pqbar, p1, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxmid(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda
) {
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
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 (aptype==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_1 = (backmidquark);
auto begin_ladder_2 = (backmidquark+2);
auto end_ladder = end(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark);
int nbelow = std::distance(begin_ladder_2, 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]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqxt, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const & decays(ev.decays());
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 pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == unordered_backward){
return tree_kin_W_unob(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == unordered_forward){
return tree_kin_W_unof(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqxb(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqxf(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
if(ev.type() == central_qqx){
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda);
}
double MatrixElement::tree_kin_Higgs(
Event const & ev
) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev);
}
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
#ifdef HEJ_BUILD_WITH_QCDLOOP
// 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;
}
#endif
// 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 HEJ_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(
+ return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*ME_Houtside_gq(
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(
Event const & ev
) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
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,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_last(
Event const & ev
) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
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(ev);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
ev
);
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,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_between(
Event const & ev
) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
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);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| 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(ev.type() == unob){
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(ev.type() == unof){
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,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s) {
const auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](auto const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())) return 1.;
switch(AWZH_boson->type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return gw*gw*gw*gw/4.;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
double MatrixElement::tree_param(
Event const & ev,
double mur
) const{
assert(is_HEJ(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*beta0*log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s)*res;
}
} // namespace HEJ
diff --git a/src/currents.cc b/src/currents.cc
index 5ac358b..3112848 100644
--- a/src/currents.cc
+++ b/src/currents.cc
@@ -1,1395 +1,1394 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/currents.hh"
#include <iostream>
#include <limits>
#include <utility>
#include <vector>
#ifdef HEJ_BUILD_WITH_QCDLOOP
#include "qcdloop/qcdloop.h"
#endif
#include "HEJ/Constants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/PDG_codes.hh"
const COM looprwfactor = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4);
constexpr double infinity = std::numeric_limits<double>::infinity();
namespace {
// Loop integrals
#ifdef HEJ_BUILD_WITH_QCDLOOP
COM B0DD(HLV q, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_B0 = [](){
ql::Bubble<std::complex<double>,double,double> ql_B0;
ql_B0.setCacheSize(100);
return ql_B0;
}();
static std::vector<double> masses(2);
static std::vector<double> momenta(1);
for(auto & m: masses) m = mq*mq;
momenta.front() = q.m2();
ql_B0.integral(result, 1, masses, momenta);
return result[0];
}
COM C0DD(HLV q1, HLV q2, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_C0 = [](){
ql::Triangle<std::complex<double>,double,double> ql_C0;
ql_C0.setCacheSize(100);
return ql_C0;
}();
static std::vector<double> masses(3);
static std::vector<double> momenta(3);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = (q1+q2).m2();
ql_C0.integral(result, 1, masses, momenta);
return result[0];
}
COM D0DD(HLV q1,HLV q2, HLV q3, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_D0 = [](){
ql::Box<std::complex<double>,double,double> ql_D0;
ql_D0.setCacheSize(100);
return ql_D0;
}();
static std::vector<double> masses(4);
static std::vector<double> momenta(6);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = q3.m2();
momenta[3] = (q1+q2+q3).m2();
momenta[4] = (q1+q2).m2();
momenta[5] = (q2+q3).m2();
ql_D0.integral(result, 1, masses, momenta);
return result[0];
}
COM A1(HLV q1, HLV q2, double mt)
// As given in Eq. (B.2) of VDD
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
if (mt < 0.)
std::cerr<<"Problem in A1! mt = "<<mt<<std::endl;
mt2=mt*mt;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 4.*mt2/Delta3*(Q2-q12-q22)
-1.-4.*q12*q22/Delta3-12.*q12*q22*Q2/Delta3/Delta3*(q12+q22-Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q2,mt)-B0DD(Q,mt) )
* ( 2.*q22/Delta3+12.*q12*q22/Delta3/Delta3*(q22-q12+Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q1,mt)-B0DD(Q,mt) )
* ( 2.*q12/Delta3+12.*q12*q22/Delta3/Delta3*(q12-q22+Q2) )
- 2./Delta3/16/M_PI/M_PI*(q12+q22-Q2);
return ans;
}
COM A2(HLV q1, HLV q2, double mt)
// As given in Eq. (B.2) of VDD, but with high energy limit
// of invariants taken.
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
if (mt < 0.)
std::cerr<<"Problem in A2! mt = "<<mt<<std::endl;
mt2=mt*mt;
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 2.*mt2+1./2.*(q12+q22-Q2)
+2.*q12*q22*Q2/Delta3 )
+looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))
*q22*(q22-q12-Q2)/Delta3
+looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))
*q12*(q12-q22-Q2)/Delta3+1./16/M_PI/M_PI;
return ans;
}
#else // no QCDloop
COM A1(HLV, HLV, double) {
throw std::logic_error{"A1 called without QCDloop support"};
}
COM A2(HLV, HLV, double) {
throw std::logic_error{"A2 called without QCDloop support"};
}
#endif
void to_current(const HLV & q, current & ret){
ret[0]=q.e();
ret[1]=q.x();
ret[2]=q.y();
ret[3]=q.z();
}
} // namespace anonymous
// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
// @TODO: this is not a current and should be moved somewhere else
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(HEJ::C_A - 1./HEJ::C_A) + 1./HEJ::C_A;
}
double K_g(
HLV const & pout,
HLV const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
CCurrent CCurrent::operator+(const CCurrent& other)
{
COM result_c0=c0 + other.c0;
COM result_c1=c1 + other.c1;
COM result_c2=c2 + other.c2;
COM result_c3=c3 + other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator-(const CCurrent& other)
{
COM result_c0=c0 - other.c0;
COM result_c1=c1 - other.c1;
COM result_c2=c2 - other.c2;
COM result_c3=c3 - other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const double x)
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const double x)
{
COM result_c0=CCurrent::c0/x;
COM result_c1=CCurrent::c1/x;
COM result_c2=CCurrent::c2/x;
COM result_c3=CCurrent::c3/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const COM x)
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const COM x)
{
COM result_c0=(CCurrent::c0)/x;
COM result_c1=(CCurrent::c1)/x;
COM result_c2=(CCurrent::c2)/x;
COM result_c3=(CCurrent::c3)/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
{
os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
return os;
}
CCurrent operator * ( double x, CCurrent& m)
{
return m*x;
}
CCurrent operator * ( COM x, CCurrent& m)
{
return m*x;
}
CCurrent operator / ( double x, CCurrent& m)
{
return m/x;
}
CCurrent operator / ( COM x, CCurrent& m)
{
return m/x;
}
COM CCurrent::dot(HLV p1)
{
// Current goes (E,px,py,pz)
// Vector goes (px,py,pz,E)
return p1[3]*c0-p1[0]*c1-p1[1]*c2-p1[2]*c3;
}
COM CCurrent::dot(CCurrent p1)
{
return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3;
}
//Current Functions
// Current for <outgoing state | mu | incoming state>
/// @TODO always use this instead of "j"
/// @TODO isn't this jio with flipt helicities?
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
const double sqpop = sqrt(pout.plus());
const double sqpom = sqrt(pout.minus());
const COM poperp = pout.x() + COM(0, 1) * pout.y();
if (helout != helin) {
throw std::invalid_argument{"Non-matching helicities"};
} else if (helout == false) { // negative helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(pin.plus());
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * poperp / abs(poperp);
cur[2] = -COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(pin.minus());
cur[0] = -sqpom * sqpim * poperp / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = COM(0,1) * cur[1];
cur[3] = -cur[0];
}
} else { // positive helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(pin.plus());
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
cur[2] = COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(pin.minus());
cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = -COM(0,1) * cur[1];
cur[3] = -cur[0];
}
}
}
CCurrent joi (HLV pout, bool helout, HLV pin, bool helin)
{
current cur;
joi(pout, helout, pin, helin, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
// Current for <incoming state | mu | outgoing state>
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) {
cur[0] = 0.0;
cur[1] = 0.0;
cur[2] = 0.0;
cur[3] = 0.0;
const double sqpop = sqrt(pout.plus());
const double sqpom = sqrt(pout.minus());
const COM poperp = pout.x() + COM(0, 1) * pout.y();
if (helout != helin) {
throw std::invalid_argument{"Non-matching helicities"};
} else if (helout == false) { // negative helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(pin.plus());
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
cur[2] = COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(pin.minus());
cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = -COM(0,1) * cur[1];
cur[3] = -cur[0];
}
} else { // positive helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(pin.plus());
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * poperp / abs(poperp);
cur[2] = -COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(pin.minus());
cur[0] = -sqpom * sqpim * poperp / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = COM(0,1) * cur[1];
cur[3] = -cur[0];
}
}
}
CCurrent jio (HLV pin, bool helin, HLV pout, bool helout)
{
current cur;
jio(pin, helin, pout, helout, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
// Current for <outgoing state | mu | outgoing state>
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {
// Zero our current
cur[0] = 0.0;
cur[1] = 0.0;
cur[2] = 0.0;
cur[3] = 0.0;
if (heli!=helj) {
throw std::invalid_argument{"Non-matching helicities"};
} else if ( heli == true ) { // If positive helicity swap momenta
std::swap(pi,pj);
}
const double sqpjp = sqrt(pj.plus());
const double sqpjm = sqrt(pj.minus());
const double sqpip = sqrt(pi.plus());
const double sqpim = sqrt(pi.minus());
const COM piperp = pi.x() + COM(0,1) * pi.y();
const COM pjperp = pj.x() + COM(0,1) * pj.y();
const COM phasei = piperp / abs(piperp);
const COM phasej = pjperp / abs(pjperp);
cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
}
CCurrent joo (HLV pi, bool heli, HLV pj, bool helj)
{
current cur;
joo(pi, heli, pj, helj, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
//@{
double j_j(HLV p1out, HLV p1in, HLV p2out, HLV p2in, bool aqlineb, bool aqlinef){
HLV q1=p1in-p1out;
HLV q2=-(p2in-p2out);
current mj1m,mj1p,mj2m,mj2p;
if(aqlineb){
jio(p1in,true,p1out,true,mj1p);
jio(p1in,false,p1out,false,mj1m);
} else{
joi(p1out,true,p1in,true,mj1p);
joi(p1out,false,p1in,false,mj1m);
}
if(aqlinef){
jio(p2in,true,p2out,true,mj2p);
jio(p2in,false,p2out,false,mj2m);
} else{
joi(p2out,true,p2in,true,mj2p);
joi(p2out,false,p2in,false,mj2m);
}
COM Mmp=cdot(mj1m,mj2p);
COM Mmm=cdot(mj1m,mj2m);
COM Mpp=cdot(mj1p,mj2p);
COM Mpm=cdot(mj1p,mj2m);
double sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*(sst)/(q1.m2()*q2.m2());
}
-double jM2qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+double ME_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false);
}
-double jM2qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+double ME_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, true);
}
-double jM2qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+double ME_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, true, true);
}
-double jM2qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+double ME_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false)*K_g(p2out, p2in)/HEJ::C_F;
}
-double jM2qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+double ME_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, true, false)*K_g(p2out, p2in)/(HEJ::C_F);
}
-double jM2gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in){
+double ME_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
return j_j(p1out, p1in, p2out, p2in, false, false)*K_g(p1out, p1in)*K_g(p2out, p2in)/(HEJ::C_F*HEJ::C_F);
}
//@}
namespace {
/**
* @brief Higgs vertex contracted with current @param C1 and @param C2
*/
COM cHdot(const current & C1, const current & C2, const current & q1,
const current & q2, double mt, bool incBot, double mb)
{
if (mt == infinity) {
return (cdot(C1,C2)*cdot(q1,q2)-cdot(C1,q2)*cdot(C2,q1))/(6*M_PI*HEJ::vev);
}
else {
HLV vq1,vq2;
vq1.set(q1[1].real(),q1[2].real(),q1[3].real(),q1[0].real());
vq2.set(q2[1].real(),q2[2].real(),q2[3].real(),q2[0].real());
// first minus sign obtained because of q1-difference to VDD
// Factor is because 4 mt^2 g^2/HEJ::vev A1 -> 16 pi mt^2/HEJ::vev alphas,
// and we divide by a factor 4 at the amp sqaured level later
// which I absorb here (i.e. I divide by 2)
/// @TODO move factor 1/2 from S to |ME|^2 => consistent with general notation
if(!(incBot))
return 8.*M_PI*mt*mt/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt));
else
return 8.*M_PI*mt*mt/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt))
+ 8.*M_PI*mb*mb/HEJ::vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)-cdot(C1,C2)*A2(-vq1,vq2,mb));
}
}
} // namespace anonymous
//@{
double j_h_j(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV q1, HLV q2, double mt, bool incBot,
double mb, bool aqlineb, bool aqlinef){
current j1p,j1m,j2p,j2m, q1v, q2v;
-
if(aqlineb){
jio (p1in,true,p1out,true,j1p);
jio (p1in,false,p1out,false,j1m);
} else {
joi (p1out,true,p1in,true,j1p);
joi (p1out,false,p1in,false,j1m);
}
if(aqlinef){
jio (p2in,true,p2out,true,j2p);
jio (p2in,false,p2out,false,j2m);
} else {
joi (p2out,true,p2in,true,j2p);
joi (p2out,false,p2in,false,j2m);
}
to_current(q1, q1v);
to_current(q2, q2v);
COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb);
COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb);
COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb);
COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb);
double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm);
return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
}
-double MH2qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
+double ME_H_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, false);
}
-double MH2qQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
+double ME_H_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, true);
}
-double MH2qbarQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
+double ME_H_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, true, false);
}
-double MH2qbarQbar (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
+double ME_H_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, true, true);
}
-double MH2qg (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
+double ME_H_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, false)*K_g(p2out,p2in)/HEJ::C_A;
}
-double MH2qbarg (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
+double ME_H_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, true, false)*K_g(p2out,p2in)/HEJ::C_A;
}
-double MH2gg (HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
+double ME_H_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV q1, HLV q2,
double mt, bool incBot, double mb){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, false, false)*(K_g(p2out,p2in)/HEJ::C_A)*(K_g(p1out,p1in)/HEJ::C_A);
}
//@}
namespace {
//@{
/// @brief Higgs vertex contracted with one current
- CCurrent jH (HLV pout, bool helout, HLV pin,
+ CCurrent jH(HLV pout, bool helout, HLV pin,
bool helin, HLV q1, HLV q2,
double mt, bool incBot, double mb)
{
CCurrent j2 = joi(pout,helout,pin,helin);
CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());
if(mt == infinity)
return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*HEJ::vev);
else
{
if(incBot)
return (-16.*M_PI*mb*mb/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mb)-16.*M_PI*mb*mb/HEJ::vev*j2*A2(-q1,q2,mb))
+ (-16.*M_PI*mt*mt/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/HEJ::vev*j2*A2(-q1,q2,mt));
else
return (-16.*M_PI*mt*mt/HEJ::vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/HEJ::vev*j2*A2(-q1,q2,mt));
}
}
//@}
} // namespace anonymous
//@{
-double j_h_juno (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1, HLV qH2,
+double j_h_juno(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1, HLV qH2,
double mt, bool incBot, double mb, bool aqlineb
){
// This construction is taking rapidity order: pg > p1out >> p2out
HLV q1=p1in-p1out; // Top End
HLV q2=-(p2in-p2out); // Bottom End
HLV qg=p1in-p1out-pg; // Extra bit post-gluon
CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p, jgam,jgap,j2gm,j2gp;
if(aqlineb){
mj1p=jio(p1in,true,p1out,true);
mj1m=jio(p1in,false,p1out,false);
j2gp=joo(pg,true,p1out,true);
j2gm=joo(pg,false,p1out,false);
jgap=jio(p1in,true,pg,true);
jgam=jio(p1in,false,pg,false);
} else {
mj1p=joi(p1out,true,p1in,true);
mj1m=joi(p1out,false,p1in,false);
j2gp=joo(p1out,true,pg,true);
j2gm=joo(p1out,false,pg,false);
jgap=joi(pg,true,p1in,true);
jgam=joi(pg,false,p1in,false);
}
mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb);
mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb);
// Dot products of these which occur again and again
COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones
COM MHmm=mj1m.dot(mjH2m);
COM MHpp=mj1p.dot(mjH2p);
COM MHpm=mj1p.dot(mjH2m);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm;
CCurrent p2o(p2out), p2i(p2in), p1o(p1out), p1i(p1in), qsum(q1+qg);
Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmm/2.) )/q1.m2();
Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmp/2.) )/q1.m2();
Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpm/2.) )/q1.m2();
Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpp/2.) )/q1.m2();
U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2();
U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2();
U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2();
U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2();
U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2();
const double cf=HEJ::C_F;
double amm,amp,apm,app;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());
// Now add the t-channels for the Higgs
double th=qH1.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
// Factor of (Cf/Ca) for each quark to match ME_H_qQ.
ampsq*=HEJ::C_F*HEJ::C_F/HEJ::C_A/HEJ::C_A;
return ampsq;
}
-double jM2unogqHQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1,
+double ME_H_unof_qQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV qH1,
HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, false);
}
-double jM2unogqbarHQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
+double ME_H_unof_qbarQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, true);
}
-double jM2unogqHQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
+double ME_H_unof_qQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, false);
}
-double jM2unogqbarHQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
+double ME_H_unof_qbarQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, true);
}
-double jM2unogqHg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
+double ME_H_unof_qg(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, false)*K_g(p2out,p2in)/HEJ::C_F;
}
-double jM2unogqbarHg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
+double ME_H_unof_qbarg(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, true)*K_g(p2out,p2in)/HEJ::C_F;
}
-double jM2unobqHQg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
+double ME_H_unob_qQ(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, false);
}
-double jM2unobqbarHQg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
+double ME_H_unob_qbarQ(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, false);
}
-double jM2unobqHQbarg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
+double ME_H_unob_qQbar(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, true);
}
-double jM2unobqbarHQbarg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
+double ME_H_unob_qbarQbar(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, true);
}
-double jM2unobgHQg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
+double ME_H_unob_gQ(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, false)*K_g(p1out,p1in)/HEJ::C_F;
}
-double jM2unobgHQbarg (HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
+double ME_H_unob_gQbar(HLV p1out, HLV p1in, HLV pg, HLV p2out, HLV p2in,
HLV qH1, HLV qH2, double mt, bool incBot, double mb){
return j_h_juno(pg, p2out, p2in, p1out, p1in, qH2, qH1, mt, incBot, mb, true)*K_g(p1out,p1in)/HEJ::C_F;
}
//@}
// Begin finite mass stuff
#ifdef HEJ_BUILD_WITH_QCDLOOP
namespace {
// All the stuff needed for the box functions in qg->qgH now...
COM E1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2=-(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) +
S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) +
2.*(s34 + S1)*(s34 + S1)/Delta +
S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
2.*(s34 + S1)/Delta +
S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) -
C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) -
C0DD(k1 + k2, q2, mq)*2.*s34/
S1 - (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 +
S1)/(S1*Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(2.*s34*(s34 +
S1)*(S1 - S2)/(Delta*Sigma) +
2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S1)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM F1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-S2*D0DD(k1, k2, q2,
mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) -
s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
S2*pow((s34 + S2),2)/Delta/Delta) - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12) -
C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2,
mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 +
S2)*(S2 - S1)/(Delta*Sigma) + (B0DD(
k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S2)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM G1(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor*(S2*D0DD(k1, q2, k2,
mq)*(Delta/s12/s12 - 4.*mq*mq/s12) -
S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) -
C0DD(k1, q2, mq))*4.*mq*mq/
s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./
s12 + (B0DD(k1 + q2, mq) -
B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1)));
}
COM E4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k2, k1, q2,
mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) -
s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S1),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta +
2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) +
s12*(S1 -
S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*((2.*s12*(2.*s12*s34 - S1*(S1 + S2) +
s34*(S2 - S1)))/(Delta*Sigma)));
}
COM F4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k1, k2, q2,
mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) +
s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S2),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 +
S2)/Delta + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 -
S1*(S1 + S2) +
s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(2.*s12*(2.*s12*s34 - S2*(S1 + S2) +
s34*(S1 - S2))/(Delta*Sigma)));
}
COM G4(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2,
mq)*(Delta/s12 + (s12 + S1)/2. -
4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD(
k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD(
k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*2./(s12 + S2));
}
COM E10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta +
12.*mq*mq*S1*(s34 + S1)/Delta/Delta -
4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) +
C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) -
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) +
8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s34*(4.*s12 + 3.*S1 +
S2)/(Delta*Sigma) +
8.*s12*s34*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2,
mq))*(12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*S1*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM F10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (s12*D0DD(k1, k2, q2,
mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta +
12.*mq*mq*s34*(s12 + S1)/Delta/Delta -
4.*s12*pow((s34 + S2),2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*(s34 + S2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) -
C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) +
4.*s12*s34*(S2 - S1)/(Delta*Sigma) +
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) - (B0DD(
k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) +
8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(-12*s34*(2*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s12*s34*s34/(S2*Delta*Delta) +
4.*s34*S1/(Delta*Sigma) -
4.*s34*(s12*s34*(2.*s12 + S2) -
S1*S1*(2.*s12 +
S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*(2.*s34 + S1)/(Delta*Sigma) -
8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) +
s12*(S2 -
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM G10(HLV k1, HLV k2, HLV kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. +
4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2,
k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*4.*(s34 +
S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) -
B0DD(k2 + q2, mq))*4.*s34/(Delta*S2));
}
COM H1(HLV k1, HLV k2, HLV kh, double mq){
return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq);
}
COM H4(HLV k1, HLV k2, HLV kh, double mq){
return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq);
}
COM H10(HLV k1, HLV k2, HLV kh, double mq){
return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq);
}
COM H2(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H1(k2,k1,kh,mq);
}
COM H5(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H4(k2,k1,kh,mq);
}
COM H12(HLV k1, HLV k2, HLV kh, double mq){
return -1.*H10(k2,k1,kh,mq);
}
// FL and FT functions
COM FL(HLV q1, HLV q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*((2.-
3.*q1.m2()*q2.dot(Q)/detQ2)*(B0DD(q1, mq) -
B0DD(Q, mq)) + (2. -
3.*q2.m2()*q1.dot(Q)/detQ2)*(B0DD(q2, mq) -
B0DD(Q, mq)) - (4.*mq*mq + q1.m2() + q2.m2() +
Q.m2() - 3.*q1.m2()*q2.m2()*Q.m2()/detQ2)*C0DD(
q1, q2, mq) - 2.);
}
COM FT(HLV q1, HLV q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*(Q.m2()*(B0DD(q1, mq) + B0DD(q2, mq) - 2.*B0DD(Q, mq) -
2.*q1.dot(q2)*C0DD(q1, q2, mq)) + (q1.m2() -
q2.m2()) *(B0DD(q1, mq) - B0DD(q2, mq))) -
q1.dot(q2)*FL(q1, q2, mq);
}
HLV ParityFlip(HLV p){
HLV flippedVector;
flippedVector.setE(p.e());
flippedVector.setX(-p.x());
flippedVector.setY(-p.y());
flippedVector.setZ(-p.z());
return flippedVector;
}
/// @brief HC amp for qg->qgH with finite top (i.e. j^{++}_H)
void g_gH_HC(HLV pa, HLV p1,
HLV pH, double mq, current &retAns)
{
current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1,
epsHapart2,conjepsH1part1,conjepsH1part2;
COM ang1a,sqa1;
const double F = 4.*mq*mq/HEJ::vev;
// Easier to have the whole thing as current object so I can use cdot functionality.
// Means I need to write pa,p1 as current objects
to_current(pa, pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
//HEJ gauge
jio(pa,false,p1,false,cura1);
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
} else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(-1./sqrt(2)/ang1a,cura1,conjeps1);
cmult(1./sqrt(2)/sqa1,cura1,epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM h4 = H4(p1,-pa,pH,mq);
const COM h5 = H5(p1,-pa,pH,mq);
const COM h10 = H10(p1,-pa,pH,mq);
const COM h12 = H12(p1,-pa,pH,mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
const COM aH1 = cdot(pHcur, cura1);
current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/ang1a, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/ang1a, epsHa, T2);
}
cmult(sqrt(2.)/ang1a*aH1, epsHa, T3);
cmult(sqrt(2.)/sqa1*aH1, conjepsH1, T4);
cmult(-sqrt(2.)*Fta*pa.dot(p1)*aH1/sqa1, conjeps1, T5);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/ang1a, epsa, T6);
cmult(-aH1/sqrt(2.)/sqa1*h4*8.*COM(0.,1.)*M_PI*M_PI, conjeps1, T7);
cmult(aH1/sqrt(2.)/ang1a*h5*8.*COM(0.,1.)*M_PI*M_PI, epsa, T8);
cmult(aH1*aH1/2./ang1a/sqa1*h10*8.*COM(0.,1.)*M_PI*M_PI, pacur, T9);
cmult(-aH1*aH1/2./ang1a/sqa1*h12*8.*COM(0.,1.)*M_PI*M_PI, p1cur, T10);
current ans;
for(int i=0;i<4;i++)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5[i]+T6[i]+T7[i]+T8[i]+T9[i]+T10[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
/// @brief HNC amp for qg->qgH with finite top (i.e. j^{+-}_H)
void g_gH_HNC(HLV pa, HLV p1, HLV pH, double mq, current &retAns)
{
const double F = 4.*mq*mq/HEJ::vev;
COM ang1a,sqa1;
current conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur,
p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1,
conjepsH1part2;
// Find here if pa, meaning the gluon, is forward or backward
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
jio(pa,true,p1,true,cura1);
joi(p1,true,pa,true,cur1a);
to_current(pa,pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
to_current(pa+p1,paplusp1cur);
to_current(p1-pa,p1minuspacur);
const COM aH1 = cdot(pHcur,cura1);
const COM oneHa = std::conj(aH1); // = cdot(pHcur,cur1a)
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
}
else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(1./sqrt(2)/sqa1, cur1a, epsa);
cmult(-1./sqrt(2)/sqa1, cura1, conjeps1);
const COM phase = cdot(conjeps1, epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM Falpha = FT(p1-pa,pa-p1-pH,mq);
const COM Fbeta = FL(p1-pa,pa-p1-pH,mq);
const COM h1 = H1(p1,-pa, pH, mq);
const COM h2 = H2(p1,-pa, pH, mq);
const COM h4 = H4(p1,-pa, pH, mq);
const COM h5 = H5(p1,-pa, pH, mq);
const COM h10 = H10(p1,-pa, pH, mq);
const COM h12 = H12(p1,-pa, pH, mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
current T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a,
T11b,T12a,T12b,T13;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/sqa1, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())*((p1.x()+COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, epsHa, T2);
}
const COM boxdiagFact = 8.*COM(0.,1.)*M_PI*M_PI;
cmult(aH1*sqrt(2.)/sqa1, epsHa, T3);
cmult(oneHa*sqrt(2.)/sqa1, conjepsH1, T4);
cmult(-2.*phase*Fta*pa.dot(pH), p1cur, T5a);
cmult(2.*phase*Ft1*p1.dot(pH), pacur, T5b);
cmult(-sqrt(2.)*Fta*p1.dot(pa)*oneHa/sqa1, conjeps1, T6);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/sqa1, epsa, T7);
cmult(-boxdiagFact*phase*h2, pacur, T8a);
cmult(boxdiagFact*phase*h1, p1cur, T8b);
cmult(boxdiagFact*aH1/sqrt(2.)/sqa1*h5, epsa, T9);
cmult(-boxdiagFact*oneHa/sqrt(2.)/sqa1*h4, conjeps1, T10);
cmult(boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h10, pacur, T11a);
cmult(-boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h12, p1cur, T11b);
cmult(-phase/(pa-p1).m2()*Falpha*(p1-pa).dot(pa-p1-pH), paplusp1cur, T12a);
cmult(phase/(pa-p1).m2()*Falpha*(pa+p1).dot(pa-p1-pH), p1minuspacur, T12b);
cmult(-phase*Fbeta*(pa-p1-pH).m2(), paplusp1cur, T13);
current ans;
for(int i=0;i<4;i++)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]+T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
} // namespace anonymous
// JDC - new amplitude with Higgs emitted close to gluon with full mt effects. Keep usual HEJ-style function call
-double MH2gq_outsideH(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pH, double mq, bool includeBottom, double mq2)
+double ME_Houtside_gq(HLV p1out, HLV p1in, HLV p2out, HLV p2in, HLV pH, double mq, bool includeBottom, double mq2)
{
current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip;
current retAns,retAnsb;
joi(p2out,true,p2in,true,cur2bplus);
joi(p2out,false,p2in,false,cur2bminus);
joi(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip);
joi(ParityFlip(p2out),false,ParityFlip(p2in),false,cur2bminusFlip);
COM app1,app2,apm1,apm2;
COM app3, app4, apm3, apm4;
if(!includeBottom)
{
g_gH_HC(p1in,p1out,pH,mq,retAns);
app1=cdot(retAns,cur2bplus);
app2=cdot(retAns,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
app3=cdot(retAns,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,retAns);
apm1=cdot(retAns,cur2bplus);
apm2=cdot(retAns,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
apm3=cdot(retAns,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip);
} else {
g_gH_HC(p1in,p1out,pH,mq,retAns);
g_gH_HC(p1in,p1out,pH,mq2,retAnsb);
app1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
app2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb);
app3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,retAns);
g_gH_HNC(p1in,p1out,pH,mq2,retAnsb);
apm1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
apm2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb);
apm3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
}
return abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1)
+ abs2(apm2) + abs2(apm3) + abs2(apm4);
}
#endif // HEJ_BUILD_WITH_QCDLOOP
double C2gHgm(HLV p2, HLV p1, HLV pH)
{
static double A=1./(3.*M_PI*HEJ::vev);
// Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
double s12,p1p,p2p;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p og the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p1p=p1.plus();
p2p=p2.plus();
} else { // opposite case
p1p=p1.minus();
p2p=p2.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=COM(0,1)*A/(2.*s12)*(p2p/p1p*conj(p1perp)*p3perp+p1p/p2p*p1perp*conj(p3perp));
temp=temp*conj(temp);
return temp.real();
}
double C2gHgp(HLV p2, HLV p1, HLV pH)
{
static double A=1./(3.*M_PI*HEJ::vev);
// Implements Eq. (4.23) in hep-ph/0301013
double s12,php,p1p,phm;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
php=pH.plus();
phm=pH.minus();
p1p=p1.plus();
} else { // opposite case
php=pH.minus();
phm=pH.plus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=-COM(0,1)*A/(2.*s12)*( conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p)
+s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm)
-pow(conj(p3perp)
+(1.+php/p1p)*conj(p1perp),2)/((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH)))) );
temp=temp*conj(temp);
return temp.real();
}
double C2qHqm(HLV p2, HLV p1, HLV pH)
{
static double A=1./(3.*M_PI*HEJ::vev);
// Implements Eq. (4.22) in hep-ph/0301013
double s12,p2p,p1p;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p2p=p2.plus();
p1p=p1.plus();
} else { // opposite case
p2p=p2.minus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=A/(2.*s12)*( sqrt(p2p/p1p)*p3perp*conj(p1perp)
+sqrt(p1p/p2p)*p1perp*conj(p3perp) );
temp=temp*conj(temp);
return temp.real();
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:50 PM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805026
Default Alt Text
(167 KB)

Event Timeline