diff --git a/include/HEJ/Constants.hh b/include/HEJ/Constants.hh
index be79bae..0bcd2ac 100644
--- a/include/HEJ/Constants.hh
+++ b/include/HEJ/Constants.hh
@@ -1,39 +1,33 @@
/** \file
* \brief Header file defining all global constants used for HEJ
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
namespace HEJ{
/// @name QCD parameters
//@{
constexpr double N_C = 3.; //!< number of Colours
constexpr double C_A = N_C; //!< \f$C_A\f$
constexpr double C_F = (N_C*N_C - 1.)/(2.*N_C); //!< \f$C_F\f$
constexpr double t_f = 0.5; //!< \f$t_f\f$
constexpr double n_f = 5.; //!< number light flavours
constexpr double beta0 = 11./3.*C_A - 4./3.*t_f*n_f; //!< \f$\beta_0\f$
//@}
-/// @name QFT parameters
-//@{
- constexpr double MW = 80.419; //!< The W mass in GeV/c^2
- constexpr double GammaW = 2.0476; //!< the W width in GeV/c^2
-
- //@}
/// @name Generation Parameters
//@{
//! @brief Default scale for virtual correction
//! \f$\lambda\f$ cf. eq. (20) in \cite Andersen:2011hs
constexpr double CLAMBDA = 0.2;
constexpr double CMINPT = 0.2; //!< minimal \f$p_t\f$ of all partons
//@}
/// @name Conventional Parameters
//@{
//! Value of first colour for colour dressing, according to LHE convention
//! \cite Boos:2001cv
constexpr int COLOUR_OFFSET = 501;
//@}
}
diff --git a/include/HEJ/MatrixElement.hh b/include/HEJ/MatrixElement.hh
index 7f9e372..e62ae7e 100644
--- a/include/HEJ/MatrixElement.hh
+++ b/include/HEJ/MatrixElement.hh
@@ -1,187 +1,187 @@
/** \file
* \brief Contains the MatrixElement Class
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <functional>
#include <vector>
#include "fastjet/PseudoJet.hh"
-#include "HEJ/PDG_codes.hh"
-#include "HEJ/Parameters.hh"
#include "HEJ/config.hh"
+#include "HEJ/Parameters.hh"
+#include "HEJ/PDG_codes.hh"
namespace CLHEP {
class HepLorentzVector;
}
namespace HEJ{
class Event;
class Particle;
//! Class to calculate the squares of matrix elements
class MatrixElement{
public:
/** \brief MatrixElement Constructor
* @param alpha_s Function taking the renormalisation scale
* and returning the strong coupling constant
* @param conf General matrix element settings
*/
MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
);
/**
* \brief squares of regulated HEJ matrix elements
* @param event The event for which to calculate matrix elements
* @returns The squares of HEJ matrix elements including virtual corrections
*
* This function returns one value for the central parameter choice
* and one additional value for each entry in \ref Event.variations().
* See eq. (22) in \cite Andersen:2011hs for the definition of the squared
* matrix element.
*
* \internal Relation to standard HEJ Met2: MatrixElement = Met2*shat^2/(pdfta*pdftb)
*/
Weights operator()(Event const & event) const;
//! Squares of HEJ tree-level matrix elements
/**
* @param event The event for which to calculate matrix elements
* @returns The squares of HEJ matrix elements without virtual corrections
*
* cf. eq. (22) in \cite Andersen:2011hs
*/
Weights tree(Event const & event) const;
/**
* \brief Virtual corrections to matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The virtual corrections to the squares of the matrix elements
*
* The all order virtual corrections to LL in the MRK limit is
* given by replacing 1/t in the scattering amplitude according to the
* lipatov ansatz.
*
* cf. second-to-last line of eq. (22) in \cite Andersen:2011hs
* note that indices are off by one, i.e. out[0].p corresponds to p_1
*/
Weights virtual_corrections(Event const & event) const;
/**
* \brief Scale-dependent part of tree-level matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The scale-dependent part of the squares of the
* tree-level matrix elements
*
* The tree-level matrix elements factorises into a renormalisation-scale
* dependent part, given by the strong coupling to some power, and a
* scale-independent remainder. This function only returns the former parts
* for the central scale choice and all \ref Event.variations().
*
* @see tree, tree_kin
*/
Weights tree_param(Event const & event) const;
/**
* \brief Kinematic part of tree-level matrix element squares
* @param event The event for which to calculate matrix elements
* @returns The kinematic part of the squares of the
* tree-level matrix elements
*
* The tree-level matrix elements factorises into a renormalisation-scale
* dependent part, given by the strong coupling to some power, and a
* scale-independent remainder. This function only returns the latter part.
* Since it does not depend on the parameter variations, only a single value
* is returned.
*
* @see tree, tree_param
*/
double tree_kin(Event const & event) const;
private:
double tree_param(
Event const & event,
double mur
) const;
double virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const;
double virtual_corrections(
Event const & event,
double mur
) const;
//! \internal cf. last line of eq. (22) in \cite Andersen:2011hs
double omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const;
double tree_kin_jets(
Event const & ev
) const;
double tree_kin_W(
Event const & ev
) const;
double tree_kin_Higgs(
Event const & ev
) const;
double tree_kin_Higgs_first(
Event const & ev
) const;
double tree_kin_Higgs_last(
Event const & ev
) const;
/**
* \internal
* \brief Higgs inbetween extremal partons.
*
* Note that in the case of unordered emission, the Higgs is *always*
* treated as if in between the extremal (FKL) partons, even if its
* rapidity is outside the extremal parton rapidities
*/
double tree_kin_Higgs_between(
Event const & ev
) const;
double tree_param_partons(
double alpha_s, double mur,
std::vector<Particle> const & partons
) const;
std::vector<int> in_extremal_jet_indices(
std::vector<fastjet::PseudoJet> const & partons
) const;
std::vector<Particle> tag_extremal_jet_partons(
Event const & ev
) const;
double MH2_forwardH(
CLHEP::HepLorentzVector const & p1out,
CLHEP::HepLorentzVector const & p1in,
pid::ParticleID type2,
CLHEP::HepLorentzVector const & p2out,
CLHEP::HepLorentzVector const & p2in,
CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const;
std::function<double (double)> alpha_s_;
MatrixElementConfig param_;
};
}
diff --git a/include/HEJ/Wjets.hh b/include/HEJ/Wjets.hh
index 9ad9d78..fedf84a 100644
--- a/include/HEJ/Wjets.hh
+++ b/include/HEJ/Wjets.hh
@@ -1,409 +1,458 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in W+Jets.
*
* This file contains all the W+Jet specific components to compute
* the current contractions for valid HEJ processes, to form a full
* W+Jets ME, currently one would have to use functions from the
* jets.hh header also. We can calculate all subleading processes for
* W+Jets.
*
* @TODO add a namespace
*/
#pragma once
#include <vector>
#include <CLHEP/Vector/LorentzVector.h>
typedef CLHEP::HepLorentzVector HLV;
+namespace HEJ{
+ class ParticleProperties;
+}
+
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_W_qQ (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_W_qg (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_W_unob_qQ (HLV p1out, HLV p1in, HLV p2out, HLV p2in,HLV pg,
+ HLV plbar, HLV pl, HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ HLV plbar, HLV pl,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ HLV plbar, HLV pl,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ HLV plbar, HLV pl,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ HLV plbar, HLV pl, HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ HLV plbar, HLV pl, HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ HLV plbar, HLV pl,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ HLV plbar, HLV pl, HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,HLV pg,
+ HLV plbar, HLV pl,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ HLV plbar, HLV pl, HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout,HLV plbar,HLV pl, HLV pqbarout,
+ HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_WExqqx_qqbarQ(HLV pgin, HLV pqout,HLV plbar,HLV pl, HLV pqbarout,
+ HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_WExqqx_qbarqg(HLV pgin, HLV pqout,HLV plbar,HLV pl, HLV pqbarout,
+ HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout,HLV plbar,HLV pl, HLV pqout,
+ HLV p2out, HLV p2in,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2, HLV p3,HLV plbar, HLV pl,
+ bool aqlinepa, HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
+ HEJ::ParticleProperties const & wprop);
//! 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
+ * @param wpro Mass and width of the W boson
* @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);
+ int nbelow, bool forwards,
+ HEJ::ParticleProperties const & wprop);
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index addb071..dad1f93 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1506 +1,1512 @@
/**
* \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 "HEJ/Constants.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/jets.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_resummable(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_resummable(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 const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & 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 const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
//! Lipatov vertex
double C2Lipatov( // B
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim,
CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom,
CLHEP::HepLorentzVector const & pop
){
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 const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
double lambda
) {
double kperp=(qav-qbv).perp();
if (kperp>lambda)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
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 pg Unordered gluon momentum
* @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
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_uno_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
assert(aptype!=21); // aptype cannot be gluon
if (bptype==21) {
if (aptype > 0)
return ME_unob_qg(pg,p1,pa,pn,pb);
else
return ME_unob_qbarg(pg,p1,pa,pn,pb);
}
else if (bptype<0) { // ----- || -----
if (aptype > 0)
return ME_unob_qQbar(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQbar(pg,p1,pa,pn,pb);
}
else { //bptype == quark
if (aptype > 0)
return ME_unob_qQ(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQ(pg,p1,pa,pn,pb);
}
}
/** 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 ME_gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_qg(pn,pb,p1,pa);
else
return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return ME_qg(p1,pa,pn,pb);
else
return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_qQ(pn,pb,p1,pa);
else
return ME_qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return ME_qQbar(p1,pa,pn,pb);
else
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
+ bool const wc, ParticleProperties const & Wprop
){
// 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);
+ return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
else
- return ME_W_qbarg(pn,plbar,pl,pb,p1,pa);
+ return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
- return ME_W_qg(p1,plbar,pl,pa,pn,pb);
+ return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
else
- return ME_W_qbarg(p1,plbar,pl,pa,pn,pb);
+ return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop);
}
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);
+ return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
- return ME_W_qQbar(pn,plbar,pl,pb,p1,pa);
+ return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
else {
if (aptype>0)
- return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa);
+ return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
- return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa);
+ return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
- return ME_W_qQ(p1,plbar,pl,pa,pn,pb);
+ return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
- return ME_W_qQbar(p1,plbar,pl,pa,pn,pb);
+ return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
else { // a is anti-quark
if (bptype > 0)
- return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb);
+ return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
- return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb);
+ return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
}
}
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
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
double ME_W_uno_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
+ bool const wc, ParticleProperties const & Wprop
){
// 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);
+ return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
- return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else { // they are both quark
if (wc) {// 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);
+ return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
- return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else{
if (aptype>0)
- return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
- return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
}
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);
+ return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else //qqbar
- return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
- return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else //qbarqbar
- return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl);
+ return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop);
}
}
}
throw std::logic_error("unknown particle types");
}
/** \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
*
* @note calculate forwards qqx contribution by reversing argument ordering.
*/
double ME_W_qqx_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 swap_q_qx, bool const wc
+ bool const swap_q_qx, bool const wc,
+ ParticleProperties const & Wprop
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
// const bool swap_q_qx= pqbar.rapidity() > pq.rapidity();
const double CFbackward = K_g( (swap_q_qx)?pq:pqbar ,pa)/HEJ::C_F;
// 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 (swap_q_qx)
- return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb)*CFbackward;
+ return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)*CFbackward;
else
- return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;
+ return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)*CFbackward;
}
else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
if (!wc){ // W Emitted from backwards qqx
if (swap_q_qx)
- return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb)*CFbackward;
+ return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop)*CFbackward;
else
- return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;
+ return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)*CFbackward;
}
else { // W Must be emitted from forwards leg.
if (swap_q_qx)
- return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, bptype<0)*CFbackward;
+ return ME_W_Exqqx_QQq(pb, pa, pn, pqbar, pq, plbar, pl, bptype<0, Wprop)*CFbackward;
else
- return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, bptype<0)*CFbackward;
+ return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, bptype<0, Wprop)*CFbackward;
}
}
throw std::logic_error("Incompatible incoming particle types with qqxb");
}
/* \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> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
- bool const wqq, bool const wc
+ bool const wqq, bool const wc,
+ ParticleProperties const & Wprop
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
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(wqq)
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons,(bptype<0),(aptype<0),
- swap_q_qx, nabove);
+ swap_q_qx, nabove, Wprop);
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, (bptype<0), (aptype<0),
- swap_q_qx, nabove, nbelow, wc);
+ swap_q_qx, nabove, nbelow, wc, Wprop);
}
/** \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, double vev
){
if (aptype==21&&bptype==21) // gg initial state
return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
else
return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
else
return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return ME_H_qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
}
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 pg 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
*
* @note This function assumes unordered gluon backwards from pa-p1 current.
* For unof, reverse call order
*/
double ME_Higgs_current_uno(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pg,
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, double vev
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
else {
if (bptype>0)
return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
else
return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev);
}
}
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_resummable(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;
}
const auto & jets = ev.jets();
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 = ev.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;
}
namespace {
template<class InIter, class partIter>
double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const double current_factor = ME_uno_current(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pa-p1-pg, pa, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
}
double MatrixElement::tree_kin_jets(Event const & ev) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if (ev.type()==HEJ::event_type::FKL){
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
);
}
else if (ev.type()==HEJ::event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
else {
throw std::logic_error("Can only reweight FKL or uno processes in Pure Jets");
}
}
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
+ double lambda, ParticleProperties const & Wprop
){
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
bool wc = aptype==partons[0].type; //leg b emits w
auto q0 = pa - p1;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
- p1, pa, plbar, pl, wc
+ p1, pa, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_uno(InIter BeginIn, partIter BeginPart,
partIter EndPart, const HLV & plbar, const HLV & pl,
- double lambda){
+ double lambda, ParticleProperties const & Wprop){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w
auto q0 = pa - p1 - pg;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_uno_current(
(BeginIn)->type, (BeginIn+1)->type, pn, pb,
- p1, pa, pg, plbar, pl, wc
+ p1, pa, pg, plbar, pl, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_W_qqx(InIter BeginIn, partIter BeginPart,
partIter EndPart, const HLV & plbar, const HLV & pl,
- double lambda){
+ double lambda, ParticleProperties const & Wprop){
const bool swap_q_qx=is_quark(*BeginPart);
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(BeginIn+1));
const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1)));
const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?1:0)));
const auto p1 = to_HepLorentzVector(*(BeginPart));
const auto pn = to_HepLorentzVector(*(EndPart-1));
const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w
auto q0 = pa - pq - pqbar;
if(!wc)
q0 -= pl + plbar;
const double current_factor = ME_W_qqx_current(
(BeginIn)->type, (BeginIn+1)->type, pa, pb,
- pq, pqbar, pn, plbar, pl, swap_q_qx, wc
+ pq, pqbar, pn, plbar, pl, swap_q_qx, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
BeginPart+2, EndPart-1,
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
+ double lambda, ParticleProperties const & Wprop
){
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;
}
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder_1 = (backmidquark);
const auto begin_ladder_2 = (backmidquark+2);
const auto end_ladder = cend(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
+ pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop
);
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() == FKL){
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
- param_.regulator_lambda);
+ param_.regulator_lambda,
+ param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_backward){
return tree_kin_W_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
- param_.regulator_lambda);
+ param_.regulator_lambda,
+ param_.ew_parameters.Wprop());
}
if(ev.type() == unordered_forward){
return tree_kin_W_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
- param_.regulator_lambda);
+ param_.regulator_lambda,
+ param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqx(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
- param_.regulator_lambda);
+ param_.regulator_lambda,
+ param_.ew_parameters.Wprop());
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqx(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
- param_.regulator_lambda);
+ param_.regulator_lambda,
+ param_.ew_parameters.Wprop());
}
assert(ev.type() == central_qqx);
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
- param_.regulator_lambda);
+ param_.regulator_lambda,
+ param_.ew_parameters.Wprop());
}
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 jets.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 const & p1out,
CLHEP::HepLorentzVector const & p1in,
ParticleID type2,
CLHEP::HepLorentzVector const & p2out,
CLHEP::HepLorentzVector const & p2in,
CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
const double vev = param_.ew_parameters.vev();
// 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*ME_Houtside_gq(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH,vev) + C2gHgm(p1in,p1out,pH,vev)
)/(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
);
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, const HLV & qH, const HLV & qHp1,
double mt, bool inc_bot, double mb, double vev){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb, vev
);
}
}
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() == FKL){
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,
param_.ew_parameters.vev()
);
}
else if(ev.type() == unob){
current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
begin(incoming), end(incoming), begin(partons),
end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
rbegin(incoming), rend(incoming), rbegin(partons),
rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
throw std::logic_error("Can only reweight FKL or uno processes in H+Jets");
}
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, double alpha_w) {
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 alpha_w*alpha_w;
// 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_resummable(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, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ
diff --git a/src/Wjets.cc b/src/Wjets.cc
index ce810f8..5490f4d 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1087 +1,1136 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
-#include "HEJ/jets.hh"
#include "HEJ/Wjets.hh"
-#include "HEJ/Tensor.hh"
-#include "HEJ/Constants.hh"
#include <array>
-
#include <iostream>
+#include "HEJ/Constants.hh"
+#include "HEJ/EWConstants.hh"
+#include "HEJ/jets.hh"
+#include "HEJ/Tensor.hh"
+
using HEJ::Tensor;
using HEJ::init_sigma_index;
using HEJ::metric;
using HEJ::rank3_current;
using HEJ::rank5_current;
using HEJ::eps;
using HEJ::to_tensor;
using HEJ::Helicity;
using HEJ::angle;
using HEJ::square;
using HEJ::flip;
+using HEJ::ParticleProperties;
namespace helicity = HEJ::helicity;
namespace { // Helper Functions
// FKL W Helper Functions
- double WProp (const HLV & plbar, const HLV & pl){
- COM propW = COM(0.,-1.)/( (pl+plbar).m2() - HEJ::MW*HEJ::MW + COM(0.,1.)*HEJ::MW*HEJ::GammaW);
+ double WProp (const HLV & plbar, const HLV & pl, ParticleProperties const & wprop){
+ COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass
+ + COM(0.,1.)*wprop.mass*wprop.width);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
namespace {
// FKL current including W emission off negative helicities
// See eq. (87) {eq:jW-} in developer manual
// Note that the terms are rearranged
Tensor<1> jW_minus(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl
){
using HEJ::helicity::minus;
const double tWin = (pa-pl-plbar).m2();
const double tWout = (p1+pl+plbar).m2();
// C++ arithmetic operators are evaluated left-to-right,
// so the following first computes complex scalar coefficients,
// which then multiply a current, reducing the number
// of multiplications
return 2.*(
+ angle(p1, pl)*square(p1, plbar)/tWout
+ square(pa, plbar)*angle(pa, pl)/tWin
)*HEJ::current(p1, pa, helicity::minus)
+ 2.*angle(p1, pl)*square(pl, plbar)/tWout
*HEJ::current(pl, pa, helicity::minus)
+ 2.*square(pa, plbar)*angle(pl, plbar)/tWin
*HEJ::current(p1, plbar, helicity::minus);
}
}
// FKL current including W emission
// see eqs. (87), (88) {eq:jW-}, {eq:jW+} in developer manual
Tensor<1> jW(
HLV const & pa, HLV const & p1,
HLV const & plbar, HLV const & pl,
Helicity h
){
if(h == helicity::minus) {
return jW_minus(pa, p1, plbar, pl);
}
return jW_minus(pa, p1, pl, plbar).complex_conj();
}
/**
* @brief W+Jets Unordered Contribution Helper Functions
* @returns result of equation (4.1.28) in Helen's Thesis (p.100)
*/
double jM2Wuno(HLV pg, HLV p1,HLV plbar, HLV pl, HLV pa, Helicity h1,
- HLV p2, HLV pb, Helicity h2, Helicity pol
- ){
+ HLV p2, HLV pb, Helicity h2, Helicity pol,
+ ParticleProperties const & wprop
+ ){
//@TODO Simplify the below (less Tensor class?)
init_sigma_index();
HLV pW = pl+plbar;
HLV q1g=pa-pW-p1-pg;
HLV q1 = pa-p1-pW;
HLV q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double tb2 = (pb-p2).m2();
const double tb2g = (pb-p2-pg).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
//use p1 as ref vec in pol tensor
Tensor<1> epsg = eps(pg,p2,pol);
Tensor<1> epsW = HEJ::current(pl,plbar,helicity::minus);
Tensor<1> j2b = HEJ::current(p2,pb,h2);
Tensor<1> Tq1q2 = to_tensor((q1+q2)/taW1 + (pb/pb.dot(pg)
+p2/p2.dot(pg)) * tb2/(2*tb2g));
Tensor<3> J31a = rank3_current(p1, pa, h1);
Tensor<2> J2_qaW =J31a.contract((pa-pW)/taW, 2);
Tensor<2> J2_p1W =J31a.contract((p1+pW)/s1W, 2);
Tensor<3> L1a = outer(Tq1q2, J2_qaW);
Tensor<3> L1b = outer(Tq1q2, J2_p1W);
Tensor<3> L2a = outer(-pg-q1,J2_qaW)/taW1;
Tensor<3> L2b = outer(-pg-q1, J2_p1W)/taW1;
Tensor<3> L3 = outer(metric(), J2_qaW.contract(pg-q2,1)+J2_p1W.contract(pg-q2,2))/taW1;
Tensor<3> L(0.);
Tensor<5> J51a = rank5_current(p1, pa, h1);
Tensor<4> J_qaW = J51a.contract((pa-pW),4);
Tensor<4> J_qag = J51a.contract(pa-pg,4);
Tensor<4> J_p1gW = J51a.contract(p1+pg+pW,4);
Tensor<3> U1a = J_qaW.contract(p1+pg,2);
Tensor<3> U1b = J_p1gW.contract(p1+pg,2);
Tensor<3> U1c = J_p1gW.contract(p1+pW,2);
Tensor<3> U1(0.);
Tensor<3> U2a = J_qaW.contract(pa-pg-pW,2);
Tensor<3> U2b = J_qag.contract(pa-pg-pW,2);
Tensor<3> U2c = J_qag.contract(p1+pW,2);
Tensor<3> U2(0.);
for(int nu=0; nu<4;nu++){
for(int mu=0;mu<4;mu++){
for(int rho=0;rho<4;rho++){
L(nu, mu, rho) = L1a(nu,mu,rho) + L1b(nu,rho,mu)
+ L2a(mu,nu,rho) + L2b(mu,rho,nu) + L3(mu,nu,rho);
U1(nu, mu, rho) = U1a(nu, mu, rho) / (s1g*taW)
+ U1b(nu,rho,mu) / (s1g*s1gW) + U1c(rho,nu,mu) / (s1W*s1gW);
U2(nu,mu,rho) = U2a(mu,nu,rho) / (taWg*taW)
+ U2b(mu,rho,nu) / (taWg*tag) + U2c(rho,mu,nu) / (s1W*tag);
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
double amp = HEJ::C_A*HEJ::C_F*HEJ::C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
double t1 = q1g.m2();
double t2 = q2.m2();
- double WPropfact = WProp(plbar, pl);
-
//Divide by WProp
- amp*=WPropfact;
+ amp*=WProp(plbar, pl, wprop);
//Divide by t-channels
amp/=(t1*t2);
return amp;
}
// Relevant Wqqx Helper Functions.
//g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
Tensor <1> gtqqxW(HLV pq,HLV pqbar,HLV pl,HLV plbar){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//blank 3 Gamma Current
Tensor<3> JV23 = rank3_current(pq,pqbar,helicity::minus);
// Components of g->qqW before W Contraction
Tensor<2> JV1 = JV23.contract((pq + pl + plbar),2)/(s2AB);
Tensor<2> JV2 = JV23.contract((pqbar + pl + plbar),2)/(s3AB);
// g->qqW Current. Note Minus between terms due to momentum flow.
// Also note: (-I)^2 from W vert. (I) from Quark prop.
Tensor<1> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
return JVCur;
}
// Helper Functions Calculate the Crossed Contribution
Tensor <2> MCrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar, HLV pl,
HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCross?
// Useful propagator factors
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double tcro1=(q3+pq).m2();
double tcro2=(q1-pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar,helicity::minus);
//Blank 5 gamma Current
Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_q3q = J523.contract((q3 + pq),2);
Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
// Components of Crossed Vertex Contribution
Tensor<3> Xcro1 = J_q3q.contract((pqbar + pl + plbar),3);
Tensor<3> Xcro2 = J_q3q.contract((q1 - pqbar),3);
Tensor<3> Xcro3 = J_2AB.contract((q1 - pqbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
Tensor<2> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
Tensor<2> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
//Initialise the Crossed Vertex Object
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro(mu,nu) = -(-Xcro1Cont(nu,mu)+Xcro2Cont(nu,mu)+Xcro3Cont(nu,mu));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar,
HLV pl, HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncross?
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pl - plbar - pq - pqbar;
double tunc1 = (q1-pq).m2();
double tunc2 = (q3+pqbar).m2();
// Define llx current.
Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//Blank 5 gamma Current
Tensor<5> J523 = rank5_current(pq,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_2AB = J523.contract((pq + pl + plbar),2);
Tensor<4> J_q1q = J523.contract((q1 - pq),2);
// 2 Contractions taken care of.
Tensor<3> Xunc1 = J_2AB.contract((q3 + pqbar),3);
Tensor<3> Xunc2 = J_q1q.contract((q3 + pqbar),3);
Tensor<3> Xunc3 = J_q1q.contract((pqbar + pl + plbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
Tensor<2> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
Tensor<2> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
//Initialise the Uncrossed Vertex Object
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc(mu,nu) = -(- Xunc1Cont(mu,nu)+Xunc2Cont(mu,nu) +Xunc3Cont(mu,nu));
}
}
return Xunc;
}
// Helper Functions Calculate the g->qqxW (Eikonal) Contributions
Tensor <2> MSymW(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
HLV pl,HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MSym?
double sa2=(pa+pq).m2();
double s12=(p1+pq).m2();
double sa3=(pa+pqbar).m2();
double s13=(p1+pqbar).m2();
double saA=(pa+pl).m2();
double s1A=(p1+pl).m2();
double saB=(pa+plbar).m2();
double s1B=(p1+plbar).m2();
double sb2=(pb+pq).m2();
double s42=(p4+pq).m2();
double sb3=(pb+pqbar).m2();
double s43=(p4+pqbar).m2();
double sbA=(pb+pl).m2();
double s4A=(p4+pl).m2();
double sbB=(pb+plbar).m2();
double s4B=(p4+plbar).m2();
double s23AB=(pl+plbar+pq+pqbar).m2();
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3=q1-pq-pqbar-plbar-pl;
double t1 = (q1).m2();
double t3 = (q3).m2();
// g->qqW Current (Factors of sqrt2 dealt with in this function.)
Tensor<1> JV = gtqqxW(pq,pqbar,pl,plbar);
// 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13+s1A+s1B))
+ pa*(t1/(sa2+sa3+saA+saB)) );
Tensor<2> X1aCont = X1a.contract(JV,3);
//4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43+s4A+s4B))
+ pb*(t3/(sb2+sb3+sbA+sbB)) );
Tensor<2> X4bCont = X4b.contract(JV,3);
//Set up each term of 3G diagram.
Tensor<3> X3g1 = outer(q1+pq+pqbar+pl+plbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar-pl-plbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(JV,3);
Tensor<2> X3g2Cont = X3g2.contract(JV,2);
Tensor<2> X3g3Cont = X3g3.contract(JV,1);
// XSym is an amalgamation of x1a, X4b and X3g.
// Makes sense from a colour factor point of view.
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym(mu,nu) = (X3g1Cont(nu,mu) + X3g2Cont(mu,nu) - X3g3Cont(nu,mu))
+ (X1aCont(mu,nu) - X4bCont(mu,nu));
}
}
return Xsym/s23AB;
}
Tensor <2> MCross(HLV pa, HLV pq, HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2> XCroCont = J323.contract((q1-pqbar),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro(mu,nu) = XCroCont(nu,mu);
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncross(HLV pa, HLV pq,HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
//Blank 3 gamma Current
Tensor<3> J323 = rank3_current(pq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2> XUncCont = J323.contract((q1-pq),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc(mu,nu) = -XUncCont(mu,nu);
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2> MSym(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
std::vector<HLV> partons, Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MsymW?
HLV q1, q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1-pq-pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
double s23 = (pq+pqbar).m2();
double sa2 = (pa+pq).m2();
double sa3 = (pa+pqbar).m2();
double s12 = (p1+pq).m2();
double s13 = (p1+pqbar).m2();
double sb2 = (pb+pq).m2();
double sb3 = (pb+pqbar).m2();
double s42 = (p4+pq).m2();
double s43 = (p4+pqbar).m2();
Tensor<1> qqxCur = HEJ::current(pq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3> X1a = outer(metric(), p1*(t1/(s12+s13))+ pa*(t1/(sa2+sa3)));
Tensor<2> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3> X4b = outer(metric(), p4*(t3/(s42+s43)) + pb*(t3/(sb2+sb3)));
Tensor<2> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3> X3g1 = outer(q1+pq+pqbar, metric());
Tensor<3> X3g2 = outer(q3-pq-pqbar, metric());
Tensor<3> X3g3 = outer(q1+q3, metric());
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym(mu, nu) = COM(0,1) * ( (X3g1Cont(nu,mu) + X3g2Cont(mu,nu)
- X3g3Cont(nu,mu)) + (X1aCont(mu,nu) - X4bCont(mu,nu)) );
}
}
return Xsym/s23;
}
//! W+Jets FKL Contributions
/**
* @brief W+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_\mu.
* Handles all possible incoming states.
*/
double jW_j( HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
- bool aqlineb, bool /* aqlinef */
+ bool aqlineb, bool /* aqlinef */,
+ ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out);
- const double WPropfact = WProp(plbar, pl);
+ const double WPropfact = WProp(plbar, pl, wprop);
const auto j_W = COM{0,-1}*jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
double Msqr = 0.;
for(const auto h: {plus, minus}) {
const auto j = HEJ::current(p2out, p2in, h);
Msqr += abs2(j_W.contract(j, 1));
}
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*Msqr/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
} // Anonymous Namespace
-double ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
- return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false);
+double ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop);
}
-double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
- return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true);
+double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true, wprop);
}
-double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
- return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false);
+double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop);
}
-double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
- return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true);
+double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true, wprop);
}
-double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
- return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false)*K_g(p2out, p2in)/HEJ::C_F;
+double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false, wprop)
+ *K_g(p2out, p2in)/HEJ::C_F;
}
-double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
- return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false)*K_g(p2out, p2in)/HEJ::C_F;
+double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false, wprop)
+ *K_g(p2out, p2in)/HEJ::C_F;
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param pg Unordered Gluon momenta
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_{uno}_\mu. Ie, unordered with W emission opposite side.
* Handles all possible incoming states.
*/
double jW_juno(HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out,
- HLV p2in, HLV pg, bool aqlineb, bool aqlinef){
+ HLV p2in, HLV pg, bool aqlineb, bool aqlinef,
+ ParticleProperties const & wprop
+ ){
using helicity::minus;
using helicity::plus;
const HLV q1=p1in-p1out-plbar-pl;
const HLV q2=-(p2in-p2out-pg);
const HLV q3=-(p2in-p2out);
const Helicity fhel = aqlinef?plus:minus;
const auto j_W = jW(p1in, p1out, plbar, pl, aqlineb?plus:minus);
const auto mj2p = HEJ::current(p2out, p2in, flip(fhel));
const auto mj2m = HEJ::current(p2out, p2in, fhel);
const auto jgbp = HEJ::current(pg, p2in, flip(fhel));
const auto jgbm = HEJ::current(pg, p2in, fhel);
const auto j2gp = HEJ::current(p2out, pg, flip(fhel));
const auto j2gm = HEJ::current(p2out, pg, fhel);
// Dot products of these which occur again and again
COM MWmp=j_W.dot(mj2p); // And now for the Higgs ones
COM MWmm=j_W.dot(mj2m);
const auto qsum = to_tensor(q2+q3);
const auto p1o = to_tensor(p1out);
const auto p1i = to_tensor(p1in);
const auto p2o = to_tensor(p2out);
const auto p2i = to_tensor(p2in);
const auto Lmm=( (-1.)*qsum*(MWmm) + (-2.*COM{j_W.dot(pg)})*mj2m + 2.*COM{mj2m.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmm/2. ) )/q3.m2();
const auto Lmp=( (-1.)*qsum*(MWmp) + (-2.*COM{j_W.dot(pg)})*mj2p + 2.*COM{mj2p.dot(pg)}*j_W
+ ( p1o/pg.dot(p1out) + p1i/pg.dot(p1in) )*( q2.m2()*MWmp/2. ) )/q3.m2();
const auto U1mm=(COM{jgbm.dot(j_W)}*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
const auto U1mp=(COM{jgbp.dot(j_W)}*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
const auto U2mm=((-1.)*COM{j2gm.dot(j_W)}*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
const auto U2mp=((-1.)*COM{j2gp.dot(j_W)}*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
double amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*abs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
- ampsq*=WProp(plbar, pl);
+ ampsq*=WProp(plbar, pl, wprop);
return ampsq/((16)*(q2.m2()*q1.m2()));
}
}
double ME_W_unob_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl,
+ ParticleProperties const & wprop
){
- return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false);
+ return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false, wprop);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl,
+ ParticleProperties const & wprop
){
- return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true);
+ return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true, wprop);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl,
+ ParticleProperties const & wprop
){
- return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false);
+ return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false, wprop);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl,
+ ParticleProperties const & wprop
){
- return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true);
+ return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true, wprop);
}
namespace{
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same side.
* Handles all possible incoming states. Note this handles both forward and back-
* -ward Wuno emission. For forward, ensure p1out is the uno and W emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean up
* ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(HLV pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
- HLV p2out, HLV p2in, bool aqlineb
+ HLV p2out, HLV p2in, bool aqlineb,
+ ParticleProperties const & wprop
){
//Calculate different Helicity choices
const Helicity h = aqlineb?helicity::plus:helicity::minus;
- double ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::plus,helicity::plus);
- double ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::plus,helicity::minus);
- double ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::minus,helicity::plus);
- double ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::minus,helicity::minus);
+ double ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
+ helicity::plus,helicity::plus, wprop);
+ double ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
+ helicity::plus,helicity::minus, wprop);
+ double ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
+ helicity::minus,helicity::plus, wprop);
+ double ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,
+ helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states
return (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
}
}
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
- return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false);
+ return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl,
+ ParticleProperties const & wprop
){
- return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false);
+ return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl,
+ ParticleProperties const & wprop
){
- return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true);
+ return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl,
+ ParticleProperties const & wprop
){
- return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true);
+ return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl, ParticleProperties const & wprop
){
- return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false)*K_g(p2out, p2in)/HEJ::C_F;
+ return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false, wprop)
+ *K_g(p2out, p2in)/HEJ::C_F;
}
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
- HLV pg, HLV plbar, HLV pl
+ HLV pg, HLV plbar, HLV pl,
+ ParticleProperties const & wprop
){
- return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true)*K_g(p2out, p2in)/HEJ::C_F;
+ return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true, wprop)
+ *K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Extremal qqx Contributions, function to handle all incoming types.
* @param pgin Incoming gluon which will split into qqx.
* @param pqout Quark of extremal qqx outgoing (W-Emission).
* @param plbar Outgoing anti-lepton momenta
* @param pl Outgoing lepton momenta
* @param pqbarout Anti-quark of extremal qqx pair. (W-Emission)
* @param pout Outgoing Particle 2 (end of FKL chain)
* @param p2in Incoming Particle 2
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W_{qqx} ^\mu j_\mu. Ie, Ex-QQX with W emission same side.
* Handles all possible incoming states. Calculated via crossing symmetry from jWuno_j.
*/
double jWqqx_j(HLV pgin, HLV pqout, HLV plbar, HLV pl,
- HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef){
+ HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef,
+ ParticleProperties const & wprop
+){
//Calculate Different Helicity Configurations.
const Helicity h = aqlinef?helicity::plus:helicity::minus;
- double ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::plus,helicity::plus);
- double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::plus,helicity::minus);
- double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::minus,helicity::plus);
- double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::minus,helicity::minus);
+ double ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
+ helicity::plus,helicity::plus, wprop);
+ double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
+ helicity::plus,helicity::minus, wprop);
+ double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
+ helicity::minus,helicity::plus, wprop);
+ double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,
+ helicity::minus,helicity::minus, wprop);
//Helicity sum and average over initial states.
double ME2 = (ME2mpp + ME2mpm + ME2mmp + ME2mmm)/(4.*HEJ::C_A*HEJ::C_A);
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout, HLV plbar, HLV pl,
- HLV pqbarout, HLV p2out, HLV p2in){
- return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false);
+ HLV pqbarout, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop);
}
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
- HLV pqout, HLV p2out, HLV p2in){
- return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true);
+ HLV pqout, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop);
}
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout, HLV plbar, HLV pl,
- HLV pqbarout, HLV p2out, HLV p2in){
- return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false)*K_g(p2out,p2in)/HEJ::C_F;
+ HLV pqbarout, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false, wprop)
+ *K_g(p2out,p2in)/HEJ::C_F;
}
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
- HLV pqout, HLV p2out, HLV p2in){
- return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true)*K_g(p2out,p2in)/HEJ::C_F;
+ HLV pqout, HLV p2out, HLV p2in,
+ ParticleProperties const & wprop
+){
+ return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true, wprop)
+ *K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm1(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below. (Less Tensor class use?)
double t1 = (p3-pb)*(p3-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p3-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,2)/t1;
return gqqCur*(-1);
}
//Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm2(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double t1 = (p2-pb)*(p2-pb);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb,refmom, helg);
Tensor<3> qqCurBlank = rank3_current(p2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(p2-pb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm3(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s23 = (p2+p3)*(p2+p3);
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<3> qqCurBlank1 = outer(p2+p3, metric())/s23;
Tensor<3> qqCurBlank2 = outer(pb, metric())/s23;
Tensor<1> Cur23 = HEJ::current(p2, p3,hel2);
Tensor<2> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2,
- HLV p3,HLV plbar, HLV pl, bool aqlinepa
+ HLV p3,HLV plbar, HLV pl, bool aqlinepa,
+ ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
// 2 independent helicity choices (complex conjugation related).
Tensor<1> TMmmm1 = qggm1(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,minus,minus, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,minus,minus, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,minus,plus, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,minus,plus, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa?plus:minus);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1);
COM Mmmm2 = TMmmm2.contract(cur1a,1);
COM Mmmm3 = TMmmm3.contract(cur1a,1);
COM Mpmm1 = TMpmm1.contract(cur1a,1);
COM Mpmm2 = TMpmm2.contract(cur1a,1);
COM Mpmm3 = TMpmm3.contract(cur1a,1);
//Colour factors:
COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
cm1m1=8./3.;
cm2m2=8./3.;
cm3m3=6.;
cm1m2 =-1./3.;
cm1m3 = -3.*COM(0.,1.);
cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
double Mmmm = real( cm1m1*pow(abs(Mmmm1),2) + cm2m2*pow(abs(Mmmm2),2)
+ cm3m3*pow(abs(Mmmm3),2) + 2.*real(cm1m2*Mmmm1*conj(Mmmm2))
+ 2.*real(cm1m3*Mmmm1*conj(Mmmm3))
+ 2.*real(cm2m3*Mmmm2*conj(Mmmm3)) );
double Mpmm = real( cm1m1*pow(abs(Mpmm1),2) + cm2m2*pow(abs(Mpmm2),2)
+ cm3m3*pow(abs(Mpmm3),2) + 2.*real(cm1m2*Mpmm1*conj(Mpmm2))
+ 2.*real(cm1m3*Mpmm1*conj(Mpmm3))
+ 2.*real(cm2m3*Mpmm2*conj(Mpmm3)) );
// Divide by WProp
- double WPropfact = WProp(plbar, pl);
-
+ const double WPropfact = WProp(plbar, pl, wprop);
return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
}
// W+Jets qqxCentral
double ME_WCenqqx_qq(HLV pa, HLV pb,HLV pl, HLV plbar, std::vector<HLV> partons,
- bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove
+ bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
+ ParticleProperties const & wprop
){
init_sigma_index();
HLV pq, pqbar, p1, p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am, T4bm, T1ap, T4bp;
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, helicity::plus);
T1am = HEJ::current(p1, pa, helicity::minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, helicity::plus);
T1am = HEJ::current(pa, p1, helicity::minus);}
if(!(aqlinepb)){
T4bp = HEJ::current(p4, pb, helicity::plus);
T4bm = HEJ::current(p4, pb, helicity::minus);}
else if(aqlinepb){
T4bp = HEJ::current(pb, p4, helicity::plus);
T4bm = HEJ::current(pb, p4, helicity::minus);}
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xsym = MSymW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is
// also the choice in qqbar helicity.
// (- - hel choice)
COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1));
COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1));
COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1));
COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1));
COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
- double WPropfact = WProp(plbar, pl);
- amp*=WPropfact;
+ amp*=WProp(plbar, pl, wprop);
return amp;
}
// no wqq emission
double ME_W_Cenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons,
bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove,
- int nbelow, bool forwards
+ int nbelow, bool forwards, ParticleProperties const & wprop
){
using helicity::minus;
using helicity::plus;
init_sigma_index();
if (!forwards){ //If Emission from Leg a instead, flip process.
std::swap(pa, pb);
std::reverse(partons.begin(),partons.end());
std::swap(aqlinepa, aqlinepb);
qqxmarker = !qqxmarker;
std::swap(nabove,nbelow);
}
HLV pq, pqbar, p1,p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
T1ap = HEJ::current(p1, pa, plus);
T1am = HEJ::current(p1, pa, minus);}
else if(aqlinepa){
T1ap = HEJ::current(pa, p1, plus);
T1am = HEJ::current(pa, p1, minus);}
Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb?plus:minus);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, minus, nabove);
Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, minus, nabove);
Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, plus, nabove);
Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, plus, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1));
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
HLV q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
- double WPropfact = WProp(plbar, pl);
- amp*=WPropfact;
+ amp*=WProp(plbar, pl, wprop);
return amp;
}