diff --git a/include/HEJ/Constants.hh b/include/HEJ/Constants.hh
index 4e2c052..f73bb41 100644
--- a/include/HEJ/Constants.hh
+++ b/include/HEJ/Constants.hh
@@ -1,33 +1,33 @@
/** \file
* \brief Header file defining all global constants used for HEJ
*
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \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 vev = 246.2196508; //!< Higgs vacuum expectation value in GeV
constexpr double gw = 0.653233;
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
//@{
constexpr double CLAMBDA = 0.2; //!< Scale for virtual correction, \f$\lambda\f$ cf. eq. (20) in \cite Andersen:2011hs
- constexpr double CMINPT = CLAMBDA; //!< minimal \f$p_t\f$ of all partons
+ constexpr double CMINPT = 0.2; //!< minimal \f$p_t\f$ of all partons
//@}
}
diff --git a/include/HEJ/MatrixElement.hh b/include/HEJ/MatrixElement.hh
index ca5a1ee..8df28b9 100644
--- a/include/HEJ/MatrixElement.hh
+++ b/include/HEJ/MatrixElement.hh
@@ -1,191 +1,191 @@
/** \file
* \brief Contains the MatrixElement Class
*
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <functional>
#include <vector>
#include "fastjet/PseudoJet.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Weights.hh"
#include "HEJ/config.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, double lambda
+ 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 p1out, CLHEP::HepLorentzVector p1in,
pid::ParticleID type2,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const;
std::function<double (double)> alpha_s_;
MatrixElementConfig param_;
};
}
diff --git a/include/HEJ/config.hh b/include/HEJ/config.hh
index 2b9d9f6..51f973c 100644
--- a/include/HEJ/config.hh
+++ b/include/HEJ/config.hh
@@ -1,175 +1,191 @@
/** \file
* \brief HEJ 2 configuration parameters
*
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <string>
#include "fastjet/JetDefinition.hh"
#include "yaml-cpp/yaml.h"
+#include "HEJ/Constants.hh"
#include "HEJ/event_types.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/optional.hh"
#include "HEJ/output_formats.hh"
#include "HEJ/ScaleFunction.hh"
namespace HEJ{
//! Jet parameters
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
double min_pt; /**< Minimum Jet Transverse Momentum */
};
//! Settings for scale variation
struct ScaleConfig{
//! Base scale choices
std::vector<ScaleFunction> base;
//! Factors for multiplicative scale variation
std::vector<double> factors;
//! Maximum ratio between renormalisation and factorisation scale
double max_ratio;
};
//! Settings for random number generator
struct RNGConfig {
//! Random number generator name
std::string name;
//! Optional initial seed
optional<std::string> seed;
};
/**! Possible treatments for fixed-order input events.
*
* The program will decide on how to treat an event based on
* the value of this enumeration.
*/
enum class EventTreatment{
reweight, /**< Perform resummation */
keep, /**< Keep the event */
discard, /**< Discard the event */
};
//! Container to store the treatments for various event types
using EventTreatMap = std::map<event_type::EventType, EventTreatment>;
/**! Input parameters.
*
* This struct handles stores all configuration parameters
* needed in a HEJ 2 run.
*
* \internal To add a new option:
* 1. Add a member to the Config struct.
* 2. Inside "src/YAMLreader.cc":
* - Add the option name to the "supported" Node in
* get_supported_options.
* - Initialise the new Config member in to_Config.
* The functions set_from_yaml (for mandatory options) and
* set_from_yaml_if_defined (non-mandatory) may be helpful.
* 3. Add a new entry (with short description) to config.yaml
* 4. Update the user documentation in "doc/Sphinx/"
*/
struct Config {
//! Parameters for scale variation
ScaleConfig scales;
//! Resummation jet properties
JetParameters resummation_jets;
//! Fixed-order jet properties
JetParameters fixed_order_jets;
//! Minimum transverse momentum for extremal partons
double min_extparton_pt;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
double max_ext_soft_pt_fraction;
+ //! The regulator lambda for the subtraction terms
+ double regulator_lambda = CLAMBDA;
//! Number of resummation configurations to generate per fixed-order event
int trials;
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
bool log_correction;
//! Event output files names and formats
std::vector<OutputFile> output;
//! Parameters for random number generation
RNGConfig rng;
//! Map to decide what to do for different event types
EventTreatMap treat;
//! Parameters for custom analyses
YAML::Node analysis_parameters;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
};
//! Configuration options for the PhaseSpacePoint class
struct PhaseSpacePointConfig {
//! Properties of resummation jets
JetParameters jet_param;
//! Minimum transverse momentum for extremal partons
double min_extparton_pt;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
double max_ext_soft_pt_fraction;
};
//! Configuration options for the MatrixElement class
struct MatrixElementConfig {
+ MatrixElementConfig() = default;
+ MatrixElementConfig(
+ bool log_correction,
+ HiggsCouplingSettings Higgs_coupling,
+ double regulator_lambda = CLAMBDA
+ ):
+ log_correction(log_correction),
+ Higgs_coupling(Higgs_coupling),
+ regulator_lambda(regulator_lambda)
+ {}
+
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
bool log_correction;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
+ //! The regulator lambda for the subtraction terms
+ double regulator_lambda = CLAMBDA;
};
//! Configuration options for the EventReweighter class
struct EventReweighterConfig {
//! Settings for phase space point generation
PhaseSpacePointConfig psp_config;
//! Settings for matrix element calculation
MatrixElementConfig ME_config;
//! Properties of resummation jets
JetParameters jet_param;
//! Treatment of the various event types
EventTreatMap treat;
};
/**! Extract PhaseSpacePointConfig from Config
*
* \internal We do not provide a PhaseSpacePointConfig constructor from Config
* so that PhaseSpacePointConfig remains an aggregate.
* This faciliates writing client code (e.g. the HEJ fixed-order generator)
* that creates a PhaseSpacePointConfig *without* a Config object.
*
* @see to_MatrixElementConfig, to_EventReweighterConfig
*/
inline
PhaseSpacePointConfig to_PhaseSpacePointConfig(Config const & conf) {
return {
conf.resummation_jets,
conf.min_extparton_pt,
conf.max_ext_soft_pt_fraction
};
}
/**! Extract MatrixElementConfig from Config
*
* @see to_PhaseSpacePointConfig, to_EventReweighterConfig
*/
inline
MatrixElementConfig to_MatrixElementConfig(Config const & conf) {
- return {conf.log_correction, conf.Higgs_coupling};
+ return {conf.log_correction, conf.Higgs_coupling, conf.regulator_lambda};
}
/**! Extract EventReweighterConfig from Config
*
* @see to_PhaseSpacePointConfig, to_MatrixElementConfig
*/
inline
EventReweighterConfig to_EventReweighterConfig(Config const & conf) {
return {
to_PhaseSpacePointConfig(conf),
to_MatrixElementConfig(conf),
conf.resummation_jets, conf.treat
};
}
} // namespace HEJ
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index e153997..17322b6 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1722 +1,1757 @@
/**
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/currents.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/utility.hh"
namespace HEJ{
//cf. last line of eq. (22) in \ref Andersen:2011hs
double MatrixElement::omega0(
double alpha_s, double mur,
- fastjet::PseudoJet const & q_j, double lambda
+ fastjet::PseudoJet const & q_j
) const {
+ const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
// use alpha_s(sqrt(q_j*lambda)), evolved to mur
return (
1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(
Event const & event
) const {
return tree(event)*virtual_corrections(event);
}
Weights MatrixElement::tree(
Event const & event
) const {
return tree_param(event)*tree_kin(event);
}
Weights MatrixElement::tree_param(
Event const & event
) const {
if(! is_HEJ(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
Weights MatrixElement::virtual_corrections(
Event const & event
) const {
if(! is_HEJ(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = virtual_corrections(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = virtual_corrections(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
double mur,
Particle const & WBoson
) const{
auto const & in = event.incoming();
const auto partons = filter_partons(event.outgoing());
fastjet::PseudoJet const & pa = in.front().p;
#ifndef NDEBUG
fastjet::PseudoJet const & pb = in.back().p;
double const norm = (in.front().p + in.back().p).E();
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
bool wc = true;
bool wqq = false;
// With extremal qqx or unordered gluon outside the extremal
// partons then it is not part of the FKL ladder and does not
// contribute to the virtual corrections. W emitted from the
// most backward leg must be taken into account in t-channel
if (event.type() == event_type::FKL) {
if (in[0].type != partons[0].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
wc = false;
}
}
else if (event.type() == event_type::qqxexb) {
q -= partons[1].p;
++first_idx;
if (abs(partons[0].type) != abs(partons[1].type)){
q -= WBoson.p;
wc = false;
}
}
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(abs(backquark->type) != abs((backquark+1)->type)) {
wqq=true;
wc=false;
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
- exponent += omega0(alpha_s, mur, q, CLAMBDA)*(
+ 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, CLAMBDA)*(
+ 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, CLAMBDA)*(
+ 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, CLAMBDA)*(
+ exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
} // namespace HEJ
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
- double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
- CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
- {
+ double C2Lipatovots(
+ CLHEP::HepLorentzVector qav,
+ CLHEP::HepLorentzVector qbv,
+ CLHEP::HepLorentzVector p1,
+ CLHEP::HepLorentzVector p2,
+ double lambda
+ ) {
double kperp=(qav-qbv).perp();
- if (kperp>HEJ::CLAMBDA)
+ if (kperp>lambda)
return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
return Cls-4./(kperp*kperp);
}
}
//! Lipatov vertex
double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip,
CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B
{
CLHEP::HepLorentzVector temptrans=-(qav+qbv);
CLHEP::HepLorentzVector p5=qav-qbv;
CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
- double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv,
- CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
- CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2)
- {
+ double C2Lipatovots(
+ CLHEP::HepLorentzVector qav,
+ CLHEP::HepLorentzVector qbv,
+ CLHEP::HepLorentzVector pa,
+ CLHEP::HepLorentzVector pb,
+ CLHEP::HepLorentzVector p1,
+ CLHEP::HepLorentzVector p2,
+ double lambda
+ ) {
double kperp=(qav-qbv).perp();
- if (kperp>HEJ::CLAMBDA)
+ if (kperp>lambda)
return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2());
else {
double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
double temp=Cls-4./(kperp*kperp);
return temp;
}
}
/** Matrix element squared for tree-level current-current scattering
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
if (aptype==21&&bptype==21) {
return jM2gg(pn,pb,p1,pa);
} else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2qg(pn,pb,p1,pa);
else
return jM2qbarg(pn,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jM2qg(p1,pa,pn,pb);
else
return jM2qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2qQ(pn,pb,p1,pa);
else
return jM2qQbar(pn,pb,p1,pa);
}
else {
if (aptype>0)
return jM2qQbar(p1,pa,pn,pb);
else
return jM2qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
double ME_W_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// We know it cannot be gg incoming.
assert(!(aptype==21 && bptype==21));
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jMWqg(pn,plbar,pl,pb,p1,pa);
else
return jMWqbarg(pn,plbar,pl,pb,p1,pa);
}
else if (bptype==21&&aptype!=21) { // ----- || -----
if (aptype > 0)
return jMWqg(p1,plbar,pl,pa,pn,pb);
else
return jMWqbarg(p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true){ // emission off b, (first argument pbout)
if (bptype>0) {
if (aptype>0)
return jMWqQ(pn,plbar,pl,pb,p1,pa);
else
return jMWqQbar(pn,plbar,pl,pb,p1,pa);
}
else {
if (aptype>0)
return jMWqbarQ(pn,plbar,pl,pb,p1,pa);
else
return jMWqbarQbar(pn,plbar,pl,pb,p1,pa);
}
}
else{ // emission off a, (first argument paout)
if (aptype > 0) {
if (bptype > 0)
return jMWqQ(p1,plbar,pl,pa,pn,pb);
else
return jMWqQbar(p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0)
return jMWqbarQ(p1,plbar,pl,pa,pn,pb);
else
return jMWqbarQbar(p1,plbar,pl,pa,pn,pb);
}
}
}
throw std::logic_error("unknown particle types");
}
/** Matrix element squared for backwards uno tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*/
double ME_W_unob_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (bptype == 21 && aptype != 21) { // b gluon => W emission off a
if (aptype > 0)
return jM2Wunogqg(pg,p1,plbar,pl,pa,pn,pb);
else
return jM2Wunogqbarg(pg,p1,plbar,pl,pa,pn,pb);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return junobMWqQg(pn,plbar,pl,pb,p1,pa,pg);
else
return junobMWqQbarg(pn,plbar,pl,pb,p1,pa,pg);
}
else{
if (aptype>0)
return junobMWqbarQg(pn,plbar,pl,pb,p1,pa,pg);
else
return junobMWqbarQbarg(pn,plbar,pl,pb,p1,pa,pg);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return jM2WunogqQ(pg,p1,plbar,pl,pa,pn,pb);
else //qqbar
return jM2WunogqQbar(pg,p1,plbar,pl,pa,pn,pb);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return jM2WunogqbarQ(pg,p1,plbar,pl,pa,pn,pb);
else //qbarqbar
return jM2WunogqbarQbar(pg,p1,plbar,pl,pa,pn,pb);
}
}
}
}
/** Matrix element squared for uno forward tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param pg Unordered gluon momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for unof Tree-Level Current-Current Scattering
*/
double ME_W_unof_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// we know they are not both gluons
if (aptype==21 && bptype!=21) {//a gluon => W emission off b
if (bptype > 0)
return jM2Wunogqg(pg, pn,plbar, pl, pb, p1, pa);
else
return jM2Wunogqbarg(pg, pn,plbar, pl, pb, p1, pa);
}
else { // they are both quark
if (wc==true) {// emission off b, i.e. b is first current
if (bptype>0){
if (aptype>0)
return jM2WunogqQ(pg,pn,plbar,pl,pb,p1,pa);
else
return jM2WunogqQbar(pg,pn,plbar,pl,pb,p1,pa);
}
else{
if (aptype>0)
return jM2WunogqbarQ(pg,pn,plbar,pl,pb,p1,pa);
else
return jM2WunogqbarQbar(pg,pn,plbar,pl,pb,p1,pa);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (aptype > 0) {
if (bptype > 0) //qq
return junofMWgqQ(pg,pn,pb,p1,plbar,pl,pa);
else //qqbar
return junofMWgqQbar(pg,pn,pb,p1,plbar,pl,pa);
}
else { // a is anti-quark
if (bptype > 0) //qbarq
return junofMWgqbarQ(pg,pn,pb,p1,plbar,pl,pa);
else //qbarqbar
return junofMWgqbarQbar(pg,pn,pb,p1,plbar,pl,pa);
}
}
}
}
/** \brief Matrix element squared for backward qqx tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param pn Final state n Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxb Tree-Level Current-Current Scattering
*/
double ME_W_qqxb_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFbackward;
if (pqbar.rapidity() > pq.rapidity()){
swapQuarkAntiquark=true;
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pq.minus())+(pq.minus())/pa.minus())+1./3.)*3./4.;
}
else{
CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pqbar.minus())+(pqbar.minus())/pa.minus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark){
return jM2Wggtoqqbarg(pa, pqbar, plbar, pl, pq, pn,pb)*CFbackward;}
else {
return jM2Wggtoqbarqg(pa, pq, plbar, pl, pqbar, pn,pb)*CFbackward;}
}
else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx
if (wc!=1){ // W Emitted from backwards qqx
if (swapQuarkAntiquark){
return jM2WgQtoqqbarQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
else{
return jM2WgQtoqbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;}
}
else { // W Must be emitted from forwards leg.
if(bptype > 0){
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;}
else{
return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;}
} else {
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;}
else{
return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl, true)*CFbackward;}
}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxb");
}
}
/* \brief Matrix element squared for forward qqx tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum
* @param pq Final state q Momentum
* @param pqbar Final state qbar Momentum
* @param p1 Final state 1 Momentum
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wc Boolean. True->W Emitted from b. Else; emitted from leg a
* @returns ME Squared for qqxf Tree-Level Current-Current Scattering
*/
double ME_W_qqxf_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
bool swapQuarkAntiquark=false;
double CFforward;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pq.plus())+(pq.plus())/pb.plus())+1./3.)*3./4.;
}
else{
CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pqbar.plus())+(pqbar.plus())/pb.plus())+1./3.)*3./4.;
}
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg
// This will be a wqqx emission as there is no other possible W Emission Site.
if (swapQuarkAntiquark){
return jM2Wggtoqqbarg(pb, pqbar, plbar, pl, pq, p1,pa)*CFforward;}
else {
return jM2Wggtoqbarqg(pb, pq, plbar, pl, pqbar, p1,pa)*CFforward;}
}
else if (bptype==21&&aptype!=21) {// b gluon => W emission off a or qqx
if (wc==1){ // W Emitted from forwards qqx
if (swapQuarkAntiquark){
return jM2WgQtoqbarqQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
else {
return jM2WgQtoqqbarQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;}
}
// W Must be emitted from backwards leg.
if (aptype > 0){
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl, false)*CFforward;}
else{
return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl, false)*CFforward;}
} else
{
if (swapQuarkAntiquark){
return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl, true)*CFforward;}
else{
return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl, true)*CFforward;}
}
}
else{
throw std::logic_error("Incompatible incoming particle types with qqxf");
}
}
/* \brief Matrix element squared for central qqx tree-level current-current scattering With W+Jets
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param nabove Number of gluons emitted before central qqxpair
* @param nbelow Number of gluons emitted after central qqxpair
* @param pa Initial state a Momentum
* @param pb Initial state b Momentum\
* @param pq Final state qbar Momentum
* @param pqbar Final state q Momentum
* @param partons Vector of all outgoing partons
* @param plbar Final state anti-lepton momentum
* @param pl Final state lepton momentum
* @param wqq Boolean. True siginfies W boson is emitted from Central qqx
* @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false.
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_W_qqxmid_current(
int aptype, int bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<HLV> partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
bool swapQuarkAntiquark=false;
if (pqbar.rapidity() < pq.rapidity()){
swapQuarkAntiquark=true;
}
double CFforward = (0.5*(3.-1./3.)*(pb.plus()/(partons[partons.size()-1].plus())+(partons[partons.size()-1].plus())/pb.plus())+1./3.)*3./4.;
double CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(partons[0].minus())+(partons[0].minus())/pa.minus())+1./3.)*3./4.;
double wt=1.;
if (aptype==21) wt*=CFbackward;
if (bptype==21) wt*=CFforward;
if (aptype <=0 && bptype <=0){ // Both External AntiQuark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons,true,true, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, false);
}
} // end both antiquark
else if (aptype<=0){ // a is antiquark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons,false,true, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove, nbelow, false);
}
} // end a is antiquark
else if (bptype<=0){ // b is antiquark
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove);
}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, false);
}
} //end b is antiquark
else{ //Both Quark or gluon
if (wqq==1){//emission from central qqbar
return wt*jM2WqqtoqQQq(pa, pb, pl, plbar, partons, false, false, swapQuarkAntiquark, nabove);}
else if (wc==1){//emission from b leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, true);
}
else { // emission from a leg
return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, false);
}
}
}
/** \brief Matrix element squared for tree-level current-current scattering with Higgs
* @param aptype Particle a PDG ID
* @param bptype Particle b PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared for Tree-Level Current-Current Scattering with Higgs
*/
double ME_Higgs_current(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype==21) // gg initial state
return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else if (aptype==21&&bptype!=21) {
if (bptype > 0)
return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.;
}
else if (bptype==21&&aptype!=21) {
if (aptype > 0)
return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
else
return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.;
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
else {
if (aptype>0)
return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.);
else
return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered forward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param punof Unordered Particle Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered forward emission
*/
double ME_Higgs_current_unof(
int aptype, int bptype,
CLHEP::HepLorentzVector const & punof,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (aptype==21&&bptype!=21) {
if (bptype > 0)
return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (bptype>0) {
if (aptype>0)
return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (aptype>0)
return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
/** \brief Current matrix element squared with Higgs and unordered backward emission
* @param aptype Particle A PDG ID
* @param bptype Particle B PDG ID
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param punob Unordered back Particle Momentum
* @param p1 Particle 1 Momentum
* @param pa Particle a Momentum
* @param qH t-channel momentum before Higgs
* @param qHp1 t-channel momentum after Higgs
* @returns ME Squared with Higgs and unordered backward emission
*/
double ME_Higgs_current_unob(
int aptype, int bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & punob,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs
CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs
double mt, bool include_bottom, double mb
){
if (bptype==21&&aptype!=21) {
if (aptype > 0)
return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else { // they are both quark
if (aptype>0) {
if (bptype>0)
return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
else {
if (bptype>0)
return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
else
return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb);
}
}
throw std::logic_error("unknown particle types");
}
CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(HEJ::MatrixElementConfig const & config) {
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace anonymous
namespace HEJ{
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
double MatrixElement::tree_kin(
Event const & ev
) const {
if(! is_HEJ(ev.type())) return 0.;
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())){
return tree_kin_jets(ev);
}
switch(AWZH_boson->type){
case pid::Higgs: {
return tree_kin_Higgs(ev);
}
case pid::Wp:
case pid::Wm:
return tree_kin_W(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
- CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn
+ 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)*C_A;
+ wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
- wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A;
+ wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
// TODO: avoid reclustering
fastjet::ClusterSequence cs(to_PseudoJet(out_partons), ev.jet_def());
const auto jets = sorted_by_rapidity(cs.inclusive_jets(ev.min_jet_pt()));
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission or qqx
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
assert(jets.size() >= 3);
++most_backward;
}
else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = cs.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(HEJ::is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
extremal_jet_idx:
no_extremal_jet_idx;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
double MatrixElement::tree_kin_jets(
Event const & ev
) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if(is_uno(ev.type())){
throw not_implemented("unordered emission not implemented for pure jets");
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4*(N_C*N_C - 1))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
- pa - p1, pa, pb, p1, pn
+ pa - p1, pa, pb, p1, pn,
+ param_.regulator_lambda
);
}
namespace{
double tree_kin_W_FKL(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
- HLV plbar, HLV pl
+ HLV plbar, HLV pl,
+ double lambda
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_current(
aptype, bptype, pn, pb,
p1, pa, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
- q0, pa, pb, p1, pn
+ q0, pa, pb, p1, pn,
+ lambda
);
return current_factor*ladder_factor;
}
double tree_kin_W_unob(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
- HLV plbar, HLV pl
+ HLV plbar, HLV pl,
+ double lambda
) {
auto pg = to_HepLorentzVector(partons[0]);
auto p1 = to_HepLorentzVector(partons[1]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - p1 -pg;
if (aptype!=partons[1].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_unob_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
- q0, pa, pb, p1, pn
+ q0, pa, pb, p1, pn,
+ lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_unof(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
- HLV plbar, HLV pl
+ HLV plbar, HLV pl,
+ double lambda
) {
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 2]);
auto pg = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_unof_current(
aptype, bptype, pn, pb,
p1, pa, pg, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
- q0, pa, pb, p1, pn
+ q0, pa, pb, p1, pn,
+ lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxb(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
- HLV plbar, HLV pl
+ HLV plbar, HLV pl,
+ double lambda
) {
HLV pq,pqbar;
if(is_quark(partons[0])){
pq = to_HepLorentzVector(partons[0]);
pqbar = to_HepLorentzVector(partons[1]);
}
else{
pq = to_HepLorentzVector(partons[1]);
pqbar = to_HepLorentzVector(partons[0]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 2;
auto end_ladder = end(partons) - 1;
bool wc = true;
auto q0 = pa - pq - pqbar;
if (partons[1].type!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_qqxb_current(
aptype, bptype, pa, pb,
pq, pqbar, pn, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
- q0, pa, pb, p1, pn
+ q0, pa, pb, p1, pn,
+ lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double tree_kin_W_qqxf(
int aptype, int bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
- HLV plbar, HLV pl
+ HLV plbar, HLV pl,
+ double lambda
) {
HLV pq,pqbar;
if(is_quark(partons[partons.size() - 1])){
pq = to_HepLorentzVector(partons[partons.size() - 1]);
pqbar = to_HepLorentzVector(partons[partons.size() - 2]);
}
else{
pq = to_HepLorentzVector(partons[partons.size() - 2]);
pqbar = to_HepLorentzVector(partons[partons.size() - 1]);
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 2;
bool wc = true;
auto q0 = pa - p1;
if (aptype!=partons[0].type) { //leg a emits w
wc = false;
q0 -=pl + plbar;
}
const double current_factor = ME_W_qqxf_current(
aptype, bptype, pa, pb,
pq, pqbar, p1, plbar, pl, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder,
- q0, pa, pb, p1, pn
+ 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
+ HLV plbar, HLV pl,
+ double lambda
) {
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wc, wqq;
if (backmidquark->type == -(backmidquark+1)->type){ // Central qqx does not emit
wqq=false;
if (aptype==partons[0].type) {
wc = true;
}
else{
wc = false;
q0-=pl+plbar;
}
}
else{
wqq = true;
wc = false;
qqxt-=pl+plbar;
}
auto begin_ladder = begin(partons) + 1;
auto end_ladder_1 = (backmidquark);
auto begin_ladder_2 = (backmidquark+2);
auto end_ladder = end(partons) - 1;
for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){
qqxt -= to_HepLorentzVector(*parton_it);
}
int nabove = std::distance(begin_ladder, backmidquark);
int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
- q0, pa, pb, p1, pn
+ q0, pa, pb, p1, pn,
+ lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
- qqxt, pa, pb, p1, pn
+ qqxt, pa, pb, p1, pn,
+ lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
auto const & decays(ev.decays());
HLV plbar, pl;
for (auto& x: decays) {
if (x.second.at(0).type < 0){
plbar = to_HepLorentzVector(x.second.at(0));
pl = to_HepLorentzVector(x.second.at(1));
}
else{
pl = to_HepLorentzVector(x.second.at(0));
plbar = to_HepLorentzVector(x.second.at(1));
}
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
if(ev.type() == unordered_backward){
return tree_kin_W_unob(incoming[0].type, incoming[1].type,
- pa, pb, partons, plbar, pl);
+ pa, pb, partons, plbar, pl,
+ param_.regulator_lambda);
}
if(ev.type() == unordered_forward){
return tree_kin_W_unof(incoming[0].type, incoming[1].type,
- pa, pb, partons, plbar, pl);
+ pa, pb, partons, plbar, pl,
+ param_.regulator_lambda);
}
if(ev.type() == extremal_qqxb){
return tree_kin_W_qqxb(incoming[0].type, incoming[1].type,
- pa, pb, partons, plbar, pl);
+ pa, pb, partons, plbar, pl,
+ param_.regulator_lambda);
}
if(ev.type() == extremal_qqxf){
return tree_kin_W_qqxf(incoming[0].type, incoming[1].type,
- pa, pb, partons, plbar, pl);
+ pa, pb, partons, plbar, pl,
+ param_.regulator_lambda);
}
if(ev.type() == central_qqx){
return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type,
- pa, pb, partons, plbar, pl);
+ pa, pb, partons, plbar, pl,
+ param_.regulator_lambda);
}
return tree_kin_W_FKL(incoming[0].type, incoming[1].type,
- pa, pb, partons, plbar, pl);
+ pa, pb, partons, plbar, pl,
+ param_.regulator_lambda);
}
double MatrixElement::tree_kin_Higgs(
Event const & ev
) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev);
}
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
#ifdef HEJ_BUILD_WITH_QCDLOOP
// TODO: code duplication with currents.cc
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == ParticleID::gluon) return K_g(pout, pin);
return C_F;
}
#endif
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == ParticleID::gluon)?C_A:C_F;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
ParticleID type2,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
// gluon case
#ifdef HEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(
Event const & ev
) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
- q0, pa, pb, p1, pn
+ 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
+ q0, pa, pb, p1, pn,
+ param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_between(
Event const & ev
) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor;
if(ev.type() == unob){
current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno"
incoming[0].type, incoming[1].type,
to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb
);
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
- q0, pa, pb, p1, pn
+ q0, pa, pb, p1, pn,
+ param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
- qH - pH, pa, pb, p1, pn
+ qH - pH, pa, pb, p1, pn,
+ param_.regulator_lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
double MatrixElement::tree_param_partons(
double alpha_s, double mur,
std::vector<Particle> const & partons
) const{
const double gs2 = 4.*M_PI*alpha_s;
double wt = std::pow(gs2, partons.size());
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(partons.size() >= 2);
for(size_t i = 1; i < partons.size()-1; ++i){
wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp());
}
}
return wt;
}
double MatrixElement::tree_param(
Event const & ev,
double mur
) const{
assert(is_HEJ(ev.type()));
const auto & out = ev.outgoing();
const double alpha_s = alpha_s_(mur);
auto AWZH_boson = std::find_if(
begin(out), end(out),
[](auto const & p){return is_AWZH_boson(p);}
);
double AWZH_coupling = 1.;
if(AWZH_boson != end(out)){
switch(AWZH_boson->type){
case pid::Higgs: {
AWZH_coupling = alpha_s*alpha_s;
break;
}
// TODO
case pid::Wp:{
AWZH_coupling = gw*gw*gw*gw/4;
break;
}
case pid::Wm:{
AWZH_coupling = gw*gw*gw*gw/4;
break;
}
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(out) + 1, end(out)})
);
}
if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons(
alpha_s, mur, filter_partons({begin(out), end(out) - 1})
);
}
return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(out));
}
} // namespace HEJ
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index 72919e4..1d51bd4 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,471 +1,475 @@
/**
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/YAMLreader.hh"
#include <algorithm>
#include <iostream>
#include <limits>
#include <map>
#include <string>
#include <unordered_map>
#include <vector>
#include <dlfcn.h>
#include "HEJ/ScaleFunction.hh"
#include "HEJ/event_types.hh"
#include "HEJ/output_formats.hh"
+#include "HEJ/Constants.hh"
namespace HEJ{
class Event;
namespace{
//! Get YAML tree of supported options
/**
* The configuration file is checked against this tree of options
* in assert_all_options_known.
*/
YAML::Node const & get_supported_options(){
const static YAML::Node supported = [](){
YAML::Node supported;
static const auto opts = {
"trials", "min extparton pt", "max ext soft pt fraction",
"FKL", "unordered", "qqx", "non-HEJ",
"scales", "scale factors", "max scale ratio", "import scales",
- "log correction", "event output", "analysis"
+ "log correction", "event output", "analysis", "regulator parameter"
};
// add subnodes to "supported" - the assigned value is irrelevant
for(auto && opt: opts) supported[opt] = "";
for(auto && jet_opt: {"min pt", "algorithm", "R"}){
supported["resummation jets"][jet_opt] = "";
supported["fixed order jets"][jet_opt] = "";
}
for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){
supported["Higgs coupling"][opt] = "";
}
for(auto && opt: {"name", "seed"}){
supported["random generator"][opt] = "";
}
return supported;
}();
return supported;
}
fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){
using namespace fastjet;
static const std::map<std::string, fastjet::JetAlgorithm> known = {
{"kt", kt_algorithm},
{"cambridge", cambridge_algorithm},
{"antikt", antikt_algorithm},
{"cambridge for passive", cambridge_for_passive_algorithm},
{"plugin", plugin_algorithm}
};
const auto res = known.find(algo);
if(res == known.end()){
throw std::invalid_argument("Unknown jet algorithm " + algo);
}
return res->second;
}
EventTreatment to_EventTreatment(std::string const & name){
static const std::map<std::string, EventTreatment> known = {
{"reweight", EventTreatment::reweight},
{"keep", EventTreatment::keep},
{"discard", EventTreatment::discard}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown event treatment " + name);
}
return res->second;
}
} // namespace anonymous
namespace detail{
void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){
setting = to_JetAlgorithm(yaml.as<std::string>());
}
void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){
setting = to_EventTreatment(yaml.as<std::string>());
}
void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){
setting = to_ParticleID(yaml.as<std::string>());
}
} // namespace detail
JetParameters get_jet_parameters(
YAML::Node const & node,
std::string const & entry
){
assert(node);
JetParameters result;
fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm;
double R;
set_from_yaml_if_defined(jet_algo, node, entry, "algorithm");
set_from_yaml(R, node, entry, "R");
result.def = fastjet::JetDefinition{jet_algo, R};
set_from_yaml(result.min_pt, node, entry, "min pt");
return result;
}
RNGConfig to_RNGConfig(
YAML::Node const & node,
std::string const & entry
){
assert(node);
RNGConfig result;
set_from_yaml(result.name, node, entry, "name");
set_from_yaml_if_defined(result.seed, node, entry, "seed");
return result;
}
HiggsCouplingSettings get_Higgs_coupling(
YAML::Node const & node,
std::string const & entry
){
assert(node);
static constexpr double mt_max = 2e4;
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(node[entry]){
throw std::invalid_argument{
"Higgs coupling settings require building HEJ 2 "
"with QCDloop support"
};
}
#endif
HiggsCouplingSettings settings;
set_from_yaml_if_defined(settings.mt, node, entry, "mt");
set_from_yaml_if_defined(settings.mb, node, entry, "mb");
set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom");
set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors");
if(settings.use_impact_factors){
if(settings.mt != std::numeric_limits<double>::infinity()){
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
else{
// huge values of the top mass are numerically unstable
settings.mt = std::min(settings.mt, mt_max);
}
return settings;
}
FileFormat to_FileFormat(std::string const & name){
static const std::map<std::string, FileFormat> known = {
{"Les Houches", FileFormat::Les_Houches},
{"HepMC", FileFormat::HepMC}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown file format " + name);
}
return res->second;
}
std::string extract_suffix(std::string const & filename){
size_t separator = filename.rfind('.');
if(separator == filename.npos) return {};
return filename.substr(separator + 1);
}
FileFormat format_from_suffix(std::string const & filename){
const std::string suffix = extract_suffix(filename);
if(suffix == "lhe") return FileFormat::Les_Houches;
if(suffix == "hepmc") return FileFormat::HepMC;
throw std::invalid_argument{
"Can't determine format for output file " + filename
};
}
void assert_all_options_known(
YAML::Node const & conf, YAML::Node const & supported
){
if(!conf.IsMap()) return;
if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"};
for(auto const & entry: conf){
const auto name = entry.first.as<std::string>();
if(! supported[name]) throw unknown_option{name};
/* check sub-options, e.g. 'resummation jets: min pt'
* we don't check analysis sub-options
* those depend on the analysis being used and should be checked there
* similar for "import scales"
*/
if(name != "analysis" && name != "import scales"){
try{
assert_all_options_known(conf[name], supported[name]);
}
catch(unknown_option const & ex){
throw unknown_option{name + ": " + ex.what()};
}
catch(invalid_type const & ex){
throw invalid_type{name + ": " + ex.what()};
}
}
}
}
} // namespace HEJ
namespace YAML {
Node convert<HEJ::OutputFile>::encode(HEJ::OutputFile const & outfile) {
Node node;
node[to_string(outfile.format)] = outfile.name;
return node;
};
bool convert<HEJ::OutputFile>::decode(Node const & node, HEJ::OutputFile & out) {
switch(node.Type()){
case NodeType::Map: {
YAML::const_iterator it = node.begin();
out.format = HEJ::to_FileFormat(it->first.as<std::string>());
out.name = it->second.as<std::string>();
return true;
}
case NodeType::Scalar:
out.name = node.as<std::string>();
out.format = HEJ::format_from_suffix(out.name);
return true;
default:
return false;
}
}
} // namespace YAML
namespace HEJ{
namespace detail{
void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){
setting = yaml.as<OutputFile>();
}
}
namespace{
void update_fixed_order_jet_parameters(
JetParameters & fixed_order_jets, YAML::Node const & yaml
){
if(!yaml["fixed order jets"]) return;
set_from_yaml_if_defined(
fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt"
);
fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm();
set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm");
double R = fixed_order_jets.def.R();
set_from_yaml_if_defined(R, yaml, "fixed order jets", "R");
fixed_order_jets.def = fastjet::JetDefinition{algo, R};
}
// like std::stod, but throw if not the whole string can be converted
double to_double(std::string const & str){
std::size_t pos;
const double result = std::stod(str, &pos);
if(pos < str.size()){
throw std::invalid_argument(str + " is not a valid double value");
}
return result;
}
using EventScale = double (*)(Event const &);
void import_scale_functions(
std::string const & file,
std::vector<std::string> const & scale_names,
std::unordered_map<std::string, EventScale> & known
) {
auto handle = dlopen(file.c_str(), RTLD_NOW);
char * error = dlerror();
if(error != nullptr) throw std::runtime_error{error};
for(auto const & scale: scale_names) {
void * sym = dlsym(handle, scale.c_str());
error = dlerror();
if(error != nullptr) throw std::runtime_error{error};
known.emplace(scale, reinterpret_cast<EventScale>(sym));
}
}
auto get_scale_map(
YAML::Node const & yaml
) {
std::unordered_map<std::string, EventScale> scale_map;
scale_map.emplace("H_T", H_T);
scale_map.emplace("max jet pperp", max_jet_pt);
scale_map.emplace("jet invariant mass", jet_invariant_mass);
scale_map.emplace("m_j1j2", m_j1j2);
if(yaml["import scales"]) {
if(! yaml["import scales"].IsMap()) {
throw invalid_type{"Entry 'import scales' is not a map"};
}
for(auto const & import: yaml["import scales"]) {
const auto file = import.first.as<std::string>();
const auto scale_names =
import.second.IsSequence()
?import.second.as<std::vector<std::string>>()
:std::vector<std::string>{import.second.as<std::string>()};
import_scale_functions(file, scale_names, scale_map);
}
}
return scale_map;
}
// simple (as in non-composite) scale functions
/**
* An example for a simple scale function would be H_T,
* H_T/2 is then composite (take H_T and then divide by 2)
*/
ScaleFunction parse_simple_ScaleFunction(
std::string const & scale_fun,
std::unordered_map<std::string, EventScale> const & known
) {
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
const auto it = known.find(scale_fun);
if(it != end(known)) return {it->first, it->second};
try{
const double scale = to_double(scale_fun);
return {scale_fun, FixedScale{scale}};
} catch(std::invalid_argument const &){}
throw std::invalid_argument{"Unknown scale choice: " + scale_fun};
}
std::string trim_front(std::string const & str){
const auto new_begin = std::find_if(
begin(str), end(str), [](char c){ return ! std::isspace(c); }
);
return std::string(new_begin, end(str));
}
std::string trim_back(std::string str){
size_t pos = str.size() - 1;
// use guaranteed wrap-around behaviour to check whether we have
// traversed the whole string
for(; pos < str.size() && std::isspace(str[pos]); --pos) {}
str.resize(pos + 1); // note that pos + 1 can be 0
return str;
}
ScaleFunction parse_ScaleFunction(
std::string const & scale_fun,
std::unordered_map<std::string, EventScale> const & known
){
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
// parse from right to left => a/b/c gives (a/b)/c
const size_t delim = scale_fun.find_last_of("*/");
if(delim == scale_fun.npos){
return parse_simple_ScaleFunction(scale_fun, known);
}
const std::string first = trim_back(std::string{scale_fun, 0, delim});
const std::string second = trim_front(std::string{scale_fun, delim+1});
if(scale_fun[delim] == '/'){
return parse_ScaleFunction(first, known)
/ parse_ScaleFunction(second, known);
}
else{
assert(scale_fun[delim] == '*');
return parse_ScaleFunction(first, known)
* parse_ScaleFunction(second, known);
}
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
using namespace event_type;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{FKL, EventTreatment::reweight},
{unob, EventTreatment::keep},
{unof, EventTreatment::keep},
{qqxexb, EventTreatment::keep},
{qqxexf, EventTreatment::keep},
{qqxmid, EventTreatment::keep},
{nonHEJ, EventTreatment::keep}
};
set_from_yaml(treat.at(FKL), yaml, "FKL");
set_from_yaml(treat.at(unob), yaml, "unordered");
treat.at(unof) = treat.at(unob);
set_from_yaml(treat.at(qqxexb), yaml, "qqx");
set_from_yaml(treat.at(qqxexf), yaml, "qqx");
set_from_yaml(treat.at(qqxmid), yaml, "qqx");
set_from_yaml(treat.at(nonHEJ), yaml, "non-HEJ");
if(treat[nonHEJ] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight non-HEJ events"};
}
return treat;
}
Config to_Config(YAML::Node const & yaml){
try{
assert_all_options_known(yaml, get_supported_options());
}
catch(unknown_option const & ex){
throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"};
}
Config config;
config.resummation_jets = get_jet_parameters(yaml, "resummation jets");
config.fixed_order_jets = config.resummation_jets;
update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
set_from_yaml(config.min_extparton_pt, yaml, "min extparton pt");
+ // Sets the standard value, then changes this if defined
+ config.regulator_lambda=CLAMBDA;
+ set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter");
config.max_ext_soft_pt_fraction = std::numeric_limits<double>::infinity();
set_from_yaml_if_defined(
config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
);
set_from_yaml(config.trials, yaml, "trials");
set_from_yaml(config.log_correction, yaml, "log correction");
config.treat = get_event_treatment(yaml);
set_from_yaml_if_defined(config.output, yaml, "event output");
config.rng = to_RNGConfig(yaml, "random generator");
set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis");
config.scales = to_ScaleConfig(yaml);
config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
} // namespace anonymous
ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
ScaleConfig config;
auto scale_funs = get_scale_map(yaml);
std::vector<std::string> scales;
set_from_yaml(scales, yaml, "scales");
config.base.reserve(scales.size());
std::transform(
begin(scales), end(scales), std::back_inserter(config.base),
[scale_funs](auto const & entry){
return parse_ScaleFunction(entry, scale_funs);
}
);
set_from_yaml_if_defined(config.factors, yaml, "scale factors");
config.max_ratio = std::numeric_limits<double>::infinity();
set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio");
return config;
}
Config load_config(std::string const & config_file){
try{
return to_Config(YAML::LoadFile(config_file));
}
catch(...){
std::cerr << "Error reading " << config_file << ":\n ";
throw;
}
}
} // namespace HEJ