Page MenuHomeHEPForge
Files F10664206

No OneTemporary
Actions

View Options

diff --git a/.clang-tidy b/.clang-tidy
index 02674ff..0f06cf4 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -1,39 +1,40 @@
Checks: '
abseil*,
android*,
boost*,
bugprone*,
cert*,
+ cppcoreguidelines*,
linuxkernel*,
llvm*,
-clang-analyzer-optin.cplusplus.VirtualCall,
-cppcoreguidelines-avoid-magic-numbers,
-cppcoreguidelines-pro-bounds-array-to-pointer-decay,
-cppcoreguidelines-pro-bounds-constant-array-index,
-fuchsia*,
-google-build-using-namespace,
-google-readability*,
-google-runtime-int,
-google-runtime-references,
-hicpp*,
-llvm-header-guard,
-modernize-use-trailing-return-type,
-readability-braces-around-statements,
'
CheckOptions:
- key: readability-identifier-naming.PrivateMemberSuffix
value: '_'
- key: readability-identifier-naming.GlobalConstantCase
value: UPPER_CASE
- key: readability-identifier-naming.GlobalVariableCase
value: UPPER_CASE
- key: misc-non-private-member-variables-in-classes.IgnoreClassesWithAllMemberVariablesBeingPublic
value: '1'
- key: llvm-namespace-comment.ShortNamespaceLines
value: '10'
- key: readability-magic-numbers.IgnoreAllFloatingPointValues
value: '1'
- key: cppcoreguidelines-special-member-functions.AllowSoleDefaultDtor
value: '1'
- key: performance-move-const-arg.CheckTriviallyCopyableMove
value: '0'
diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh
index 33f151e..b08f39d 100644
--- a/include/HEJ/Config.hh
+++ b/include/HEJ/Config.hh
@@ -1,234 +1,234 @@
/** \file
* \brief HEJ 2 configuration parameters
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <map>
#include <string>
#include <utility>
#include <vector>
#include "fastjet/JetDefinition.hh"
#include "yaml-cpp/yaml.h"
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/Fraction.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/ScaleFunction.hh"
#include "HEJ/event_types.hh"
#include "HEJ/optional.hh"
#include "HEJ/output_formats.hh"
namespace HEJ{
//! Jet parameters
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
- double min_pt; /**< Minimum Jet Transverse Momentum */
+ 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;
+ double max_ratio{};
};
//! Settings for random number generator
struct RNGConfig {
//! Random number generator name
std::string name;
//! Optional initial seed
optional<std::string> seed;
};
//! Settings for partial unweighting
struct PartialUnweightConfig {
//! Number of trials for training
size_t trials;
//! Maximum distance in standard deviations from mean logarithmic weight
double max_dev;
};
/**! 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>;
//! Possible setting for the event weight
enum class WeightType{
weighted, //!< weighted events
unweighted_resum, //!< unweighted only resummation part
partially_unweighted //!< mixed weighted and unweighted
};
/**! 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
//! \deprecated This will be removed in future versions.
//! Use \ref max_ext_soft_pt_fraction instead.
double min_extparton_pt = 0;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
- Fraction<double> max_ext_soft_pt_fraction;
+ Fraction<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
- size_t trials;
+ size_t trials{};
//! Maximal number of events
optional<size_t> max_events;
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
- bool log_correction;
+ 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 analysis
//! @deprecated use analyses_parameters instead
YAML::Node analysis_parameters;
//! %Parameters for custom analyses
std::vector<YAML::Node> analyses_parameters;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
//! elector weak parameters
EWConstants ew_parameters;
//! Type of event weight e.g. (un)weighted
WeightType weight_type;
//! Settings for partial unweighting
optional<PartialUnweightConfig> unweight_config;
};
//! Configuration options for the PhaseSpacePoint class
struct PhaseSpacePointConfig {
//! Properties of resummation jets
JetParameters jet_param;
//! Minimum transverse momentum for extremal partons
//! \deprecated This will be removed in future versions.
//! Use \ref max_ext_soft_pt_fraction instead.
double min_extparton_pt = 0;
//! Maximum transverse momentum fraction from soft radiation in extremal jets
Fraction<double> max_ext_soft_pt_fraction;
};
//! Configuration options for the MatrixElement class
struct MatrixElementConfig {
MatrixElementConfig() = default;
MatrixElementConfig(
bool log_correction,
HiggsCouplingSettings Higgs_coupling,
EWConstants ew_parameters,
double regulator_lambda = CLAMBDA
):
log_correction{log_correction},
Higgs_coupling{std::move(Higgs_coupling)},
ew_parameters{std::move(ew_parameters)},
regulator_lambda{regulator_lambda}
{}
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
- bool log_correction;
+ bool log_correction{};
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
//! elector weak parameters
EWConstants ew_parameters;
//! 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;
//! Access properties of resummation jets
JetParameters & jet_param() {
return psp_config.jet_param;}
//! Access properties of resummation jets (const version)
JetParameters const & jet_param() const {
return psp_config.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,
conf.ew_parameters, 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.treat
};
}
} // namespace HEJ
diff --git a/include/HEJ/EWConstants.hh b/include/HEJ/EWConstants.hh
index ff27185..0a6cb08 100644
--- a/include/HEJ/EWConstants.hh
+++ b/include/HEJ/EWConstants.hh
@@ -1,90 +1,90 @@
/** \file
* \brief Defines the electro weak parameters
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <cmath>
#include "HEJ/PDG_codes.hh"
#include "HEJ/exceptions.hh"
namespace HEJ {
//! collection of basic particle properties
struct ParticleProperties {
double mass; //!< Mass
double width; //!< Decay width
};
//! Collection of electro-weak constants
class EWConstants {
public:
EWConstants() = default;
//! initialise by Vacuum expectation value & boson properties
EWConstants(
double vev, //!< vacuum expectation value
ParticleProperties const & Wprop, //!< W boson mass & width
ParticleProperties const & Zprop, //!< Z boson mass & width
ParticleProperties const & Hprop //!< Higgs boson mass & width
): set{true}, vev_{vev}, Wprop_{Wprop}, Zprop_{Zprop}, Hprop_{Hprop} {}
//! set constants by Vacuum expectation value & boson properties
void set_vevWZH(
double vev, //!< vacuum expectation value
ParticleProperties const & Wprop, //!< W boson mass & width
ParticleProperties const & Zprop, //!< Z boson mass & width
ParticleProperties const & Hprop //!< Higgs boson mass & width
){
set = true; vev_= vev; Wprop_= Wprop; Zprop_= Zprop; Hprop_= Hprop;
}
//! vacuum expectation value
double vev() const {check_set(); return vev_;}
//! Properties of the W boson
ParticleProperties const & Wprop() const {check_set(); return Wprop_;}
//! Properties of the Z boson
ParticleProperties const & Zprop() const {check_set(); return Zprop_;}
//! Properties of the Higgs boson
ParticleProperties const & Hprop() const {check_set(); return Hprop_;}
//! access Properties by boson id
ParticleProperties const & prop(ParticleID const id) const {
using namespace pid;
switch(id){
case Wp:
case Wm:
return Wprop();
case Z:
return Zprop();
case h:
return Hprop();
default:
throw std::invalid_argument("No properties available for particle "+name(id));
}
}
//! cosine of Weinberg angle
double cos_tw() const {return Wprop().mass/Zprop().mass;}
//! cosine square of Weinberg angle
double cos2_tw() const {return cos_tw()*cos_tw();}
//! sinus Weinberg angle
double sin_tw() const {return std::sqrt(sin2_tw());}
//! sinus square of Weinberg angle
double sin2_tw() const {return 1. - cos2_tw();}
//! elector magnetic coupling
double alpha_em() const {return e2()/4./M_PI;}
//! weak coupling
double alpha_w() const {return gw2()/2.;}
private:
double gw2() const {return 4*Wprop().mass/vev()*Wprop().mass/vev();}
double e2() const {return gw2()*sin2_tw();}
void check_set() const {
if(!set) throw std::invalid_argument("EW constants not specified");
}
bool set{false};
- double vev_;
- ParticleProperties Wprop_;
- ParticleProperties Zprop_;
- ParticleProperties Hprop_;
+ double vev_{};
+ ParticleProperties Wprop_{};
+ ParticleProperties Zprop_{};
+ ParticleProperties Hprop_{};
};
} // namespace HEJ
diff --git a/include/HEJ/HepMC2Interface.hh b/include/HEJ/HepMC2Interface.hh
index e678c40..c89e91f 100644
--- a/include/HEJ/HepMC2Interface.hh
+++ b/include/HEJ/HepMC2Interface.hh
@@ -1,78 +1,78 @@
/** \file
* \brief Header file for the HepMC2Interface
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <cstddef>
#include "HEJ/PDG_codes.hh"
namespace HepMC {
class GenCrossSection;
class GenEvent;
}
namespace LHEF {
class HEPRUP;
}
namespace HEJ {
class Event;
//! This class converts the Events into HepMC::GenEvents
/**
* \details The output is depended on the HepMC version HEJ is compiled with,
* both HepMC 2 and HepMC 3 are supported. If HEJ 2 is compiled
* without HepMC calling this interface will throw an error.
*
* This interface will also keep track of the cross section of all the events that
* being fed into it.
*/
class HepMC2Interface{
public:
HepMC2Interface(LHEF::HEPRUP const &);
~HepMC2Interface();
/**
* \brief main function to convert an event into HepMC::GenEvent
*
* \param event Event to convert
* \param weight_index optional selection of specific weight
* (negative value gives central weight)
*/
HepMC::GenEvent operator()(Event const & event,
int weight_index = -1);
/**
* \brief Sets the central value from \p event to \p out_ev
*
* \param out_ev HepMC::GenEvent to write to
* \param event Event to convert
* \param weight_index optional selection of specific weight
* (negative value gives "central")
*
* This overwrites the central value of \p out_ev.
*/
void set_central(HepMC::GenEvent & out_ev, Event const & event,
int weight_index = -1);
protected:
/**
* \brief initialise generic event infomrations (not central weights)
*
* \param event Event to convert
*/
HepMC::GenEvent init_event(Event const & event) const;
private:
- std::array<ParticleID,2> beam_particle_;
- std::array<double,2> beam_energy_;
+ std::array<ParticleID,2> beam_particle_{};
+ std::array<double,2> beam_energy_{};
std::size_t event_count_;
double tot_weight_;
double tot_weight2_;
HepMC::GenCrossSection cross_section() const;
};
} // namespace HEJ
diff --git a/include/HEJ/HepMC3Interface.hh b/include/HEJ/HepMC3Interface.hh
index 9ecdca5..fad936a 100644
--- a/include/HEJ/HepMC3Interface.hh
+++ b/include/HEJ/HepMC3Interface.hh
@@ -1,82 +1,82 @@
/** \file
* \brief Header file for the HepMC3Interface
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <cstddef>
#include <memory>
#include "HEJ/PDG_codes.hh"
namespace HepMC3 {
class GenCrossSection;
class GenEvent;
class GenRunInfo;
}
namespace LHEF {
class HEPRUP;
}
namespace HEJ {
class Event;
//! This class converts the Events into HepMC3::GenEvents
/**
* \details The output is depended on the HepMC3 version HEJ is compiled with,
* both HepMC3 2 and HepMC3 3 are supported. If HEJ 2 is compiled
* without HepMC3 calling this interface will throw an error.
*
* This interface will also keep track of the cross section of all the events that
* being fed into it.
*/
class HepMC3Interface{
public:
HepMC3Interface(LHEF::HEPRUP const &);
~HepMC3Interface();
/**
* \brief main function to convert an event into HepMC3::GenEvent
*
* \param event Event to convert
* \param weight_index optional selection of specific weight
* (negative value gives central weight)
*/
HepMC3::GenEvent operator()(Event const & event, int weight_index = -1);
/**
* \brief Sets the central value from \p event to \p out_ev
*
* \param out_ev HepMC3::GenEvent to write to
* \param event Event to convert
* \param weight_index optional selection of specific weight
* (negative value gives "central")
*
* This overwrites the central value of \p out_ev.
*/
void set_central(HepMC3::GenEvent & out_ev, Event const & event,
int weight_index = -1);
//! Pointer to generic run informations
std::shared_ptr<HepMC3::GenRunInfo> run_info;
protected:
/**
* \brief initialise generic event infomrations (not central weights)
*
* \param event Event to convert
*/
HepMC3::GenEvent init_event(Event const & event) const;
private:
- std::array<ParticleID,2> beam_particle_;
- std::array<double,2> beam_energy_;
+ std::array<ParticleID,2> beam_particle_{};
+ std::array<double,2> beam_energy_{};
std::size_t event_count_;
double tot_weight_;
double tot_weight2_;
std::shared_ptr<HepMC3::GenCrossSection> xs_;
};
} // namespace HEJ
diff --git a/include/HEJ/Parameters.hh b/include/HEJ/Parameters.hh
index f9a91f1..a2febdc 100644
--- a/include/HEJ/Parameters.hh
+++ b/include/HEJ/Parameters.hh
@@ -1,216 +1,216 @@
/** \file
* \brief Containers for Parameter variations, e.g. different Weights
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#pragma once
#include <cstddef>
#include <memory>
#include <string>
#include <utility>
#include <vector>
#include "HEJ/exceptions.hh"
namespace HEJ{
//! Collection of parameters, e.g. Weights, assigned to a single event
/**
* A number of member functions of the MatrixElement class return Parameters
* objects containing the squares of the matrix elements for the various
* scale choices.
*/
template<class T>
struct Parameters {
- T central;
+ T central{};
std::vector<T> variations;
template<class T_ext>
Parameters<T>& operator*=(Parameters<T_ext> const & other);
Parameters<T>& operator*=(double factor);
template<class T_ext>
Parameters<T>& operator/=(Parameters<T_ext> const & other);
Parameters<T>& operator/=(double factor);
template<class T_ext>
Parameters<T>& operator+=(Parameters<T_ext> const & other);
template<class T_ext>
Parameters<T>& operator-=(Parameters<T_ext> const & other);
};
template<class T1, class T2> inline
Parameters<T1> operator*(Parameters<T1> a, Parameters<T2> const & b) {
a*=b;
return a;
}
template<class T> inline
Parameters<T> operator*(Parameters<T> a, double b) {
a*=b;
return a;
}
template<class T> inline
Parameters<T> operator*(double b, Parameters<T> a) {
a*=b;
return a;
}
template<class T1, class T2> inline
Parameters<T1> operator/(Parameters<T1> a, Parameters<T2> const & b) {
a/=b;
return a;
}
template<class T> inline
Parameters<T> operator/(Parameters<T> a, double b) {
a/=b;
return a;
}
template<class T1, class T2> inline
Parameters<T1> operator+(Parameters<T1> a, Parameters<T2> const & b) {
a+=b;
return a;
}
template<class T1, class T2> inline
Parameters<T1> operator-(Parameters<T1> a, Parameters<T2> const & b) {
a-=b;
return a;
}
//! Alias for weight container, e.g. used by the MatrixElement
using Weights = Parameters<double>;
//! Description of event parameters, see also EventParameters
struct ParameterDescription {
//! Name of central scale choice (e.g. "H_T/2")
std::string scale_name;
//! Actual renormalisation scale divided by central scale
- double mur_factor;
+ double mur_factor{};
//! Actual factorisation scale divided by central scale
- double muf_factor;
+ double muf_factor{};
ParameterDescription() = default;
ParameterDescription(
std::string scale_name, double mur_factor, double muf_factor
):
scale_name{std::move(scale_name)},
mur_factor{mur_factor}, muf_factor{muf_factor}
{}
};
//! generate human readable string name
std::string to_string(ParameterDescription const & p);
//! generate simplified string, intended for easy parsing
//! Format: Scale_SCALENAME_MuRxx_MuFyy
std::string to_simple_string(ParameterDescription const & p);
//! Event parameters
struct EventParameters{
- double mur; /**< Value of the Renormalisation Scale */
- double muf; /**< Value of the Factorisation Scale */
- double weight; /**< Event Weight */
+ double mur{}; /**< Value of the Renormalisation Scale */
+ double muf{}; /**< Value of the Factorisation Scale */
+ double weight{}; /**< Event Weight */
//! Optional description
std::shared_ptr<ParameterDescription> description = nullptr;
//! multiply weight by factor
EventParameters& operator*=(double factor){
weight*=factor;
return *this;
}
//! divide weight by factor
EventParameters& operator/=(double factor){
weight/=factor;
return *this;
}
//! add number to weight
EventParameters& operator+=(double factor){
weight+=factor;
return *this;
}
//! add number to weight
EventParameters& operator-=(double factor){
weight-=factor;
return *this;
}
};
inline EventParameters operator*(EventParameters a, double b){
a*=b;
return a;
}
inline EventParameters operator*(double b, EventParameters a){
a*=b;
return a;
}
inline EventParameters operator/(EventParameters a, double b){
a/=b;
return a;
}
//! @{
//! @internal Implementation of template functions
template<class T>
template<class T_ext>
Parameters<T>& Parameters<T>::operator*=(Parameters<T_ext> const & other) {
if(other.variations.size() != variations.size()) {
throw std::invalid_argument{"Wrong number of Parameters"};
}
central *= other.central;
for(std::size_t i = 0; i < variations.size(); ++i) {
variations[i] *= other.variations[i];
}
return *this;
}
template<class T>
Parameters<T>& Parameters<T>::operator*=(double factor) {
central *= factor;
for(auto & wt: variations) wt *= factor;
return *this;
}
template<class T>
template<class T_ext>
Parameters<T>& Parameters<T>::operator/=(Parameters<T_ext> const & other) {
if(other.variations.size() != variations.size()) {
throw std::invalid_argument{"Wrong number of Parameters"};
}
central /= other.central;
for(std::size_t i = 0; i < variations.size(); ++i) {
variations[i] /= other.variations[i];
}
return *this;
}
template<class T>
Parameters<T>& Parameters<T>::operator/=(double factor) {
central /= factor;
for(auto & wt: variations) wt /= factor;
return *this;
}
template<class T>
template<class T_ext>
Parameters<T>& Parameters<T>::operator+=(Parameters<T_ext> const & other) {
if(other.variations.size() != variations.size()) {
throw std::invalid_argument{"Wrong number of Parameters"};
}
central += other.central;
for(std::size_t i = 0; i < variations.size(); ++i) {
variations[i] += other.variations[i];
}
return *this;
}
template<class T>
template<class T_ext>
Parameters<T>& Parameters<T>::operator-=(Parameters<T_ext> const & other) {
if(other.variations.size() != variations.size()) {
throw std::invalid_argument{"Wrong number of Parameters"};
}
central -= other.central;
for(std::size_t i = 0; i < variations.size(); ++i) {
variations[i] -= other.variations[i];
}
return *this;
}
//! @}
} // namespace HEJ
diff --git a/include/HEJ/ScaleFunction.hh b/include/HEJ/ScaleFunction.hh
index db49349..29aa9d9 100644
--- a/include/HEJ/ScaleFunction.hh
+++ b/include/HEJ/ScaleFunction.hh
@@ -1,174 +1,174 @@
/** \file
* \brief Functions to calculate the (renormalisation and factorisation) scales for an event
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <functional>
#include <memory>
#include <string>
#include <utility>
#include <vector>
namespace HEJ{
class Event;
//! Class to calculate the scale associated with an event
class ScaleFunction {
public:
//! Constructor
/**
* @param name Name of the scale choice (e.g. H_T)
* @param fun Function used to calculate the scale
*/
ScaleFunction(std::string name, std::function<double(Event const &)> fun):
name_{std::move(name)},
fun_{std::move(fun)}
{}
//! Name of the scale choice
std::string const & name() const {
return name_;
}
//! Calculate the scale associated with an event
double operator()(Event const & ev) const {
return fun_(ev);
}
private:
friend ScaleFunction operator*(double factor, ScaleFunction base_scale);
friend ScaleFunction operator/(ScaleFunction base_scale, double denom);
friend ScaleFunction operator*(ScaleFunction func1, ScaleFunction func2);
friend ScaleFunction operator/(ScaleFunction func1, ScaleFunction func2);
std::string name_;
std::function<double(Event const &)> fun_;
};
//! Multiply a scale choice by a constant factor
/**
* For example, multiplying 0.5 and a scale function for H_T
* will result in a scale function for H_T/2.
*/
ScaleFunction operator*(double factor, ScaleFunction base_scale);
//! Multiply a scale choice by a second one
/**
* For example, multiplying H_T and m_j1j2
* will result in a scale function for H_T*m_j1j2.
*/
ScaleFunction operator*(ScaleFunction factor, ScaleFunction base_scale);
//! Divide a scale choice by a constant factor
/**
* For example, dividing a scale function for H_T by 2
* will result in a scale function for H_T/2.
*/
ScaleFunction operator/(ScaleFunction base_scale, double denom);
//! Divide a scale choice by a second one
/**
* For example, dividing a scale function for H_T by m_j1j2
* will result in a scale function for H_T/m_j1j2.
*/
ScaleFunction operator/(ScaleFunction base_scale, ScaleFunction denom);
//! Calculate \f$H_T\f$ for the input event
/**
* \f$H_T\f$ is the sum of the (scalar) transverse momenta of all
* final-state particles
*/
double H_T(Event const &);
//! The maximum of all (scalar) jet transverse momentum
double max_jet_pt(Event const &);
//! The invariant mass of the sum of all jet momenta
double jet_invariant_mass(Event const &);
//! Invariant mass of the two hardest jets
double m_j1j2(Event const &);
//! Functor that returns a fixed scale regardless of the input event
class FixedScale {
public:
explicit FixedScale(double mu): mu_{mu} {}
double operator()(Event const &) const {
return mu_;
}
private:
double mu_;
};
struct ParameterDescription;
//! Generate combinations of renormalisation and factorisation scales
class ScaleGenerator{
public:
ScaleGenerator() = default;
/** \brief Constructor
* @param scale_functions_begin Iterator to first base scale
* @param scale_functions_end Iterator past last base scale
* @param scale_factors_begin Iterator to first scale factor
* @param scale_factors_end Iterator past last scale factor
* @param max_scale_ratio Maximum ratio between renormalisation
* and factorisation scale
*/
template<class ScaleFunIterator, class FactorIterator>
ScaleGenerator(
ScaleFunIterator scale_functions_begin,
ScaleFunIterator scale_functions_end,
FactorIterator scale_factors_begin,
FactorIterator scale_factors_end,
double max_scale_ratio
):
scales_(scale_functions_begin, scale_functions_end),
scale_factors_(scale_factors_begin, scale_factors_end),
max_scale_ratio_{max_scale_ratio}
{
gen_descriptions();
}
/** \brief Constructor
* @param scales Base scales
* @param scale_factors Factors to multiply the base scales
* @param max_scale_ratio Maximum ratio between renormalisation
* and factorisation scale
*/
ScaleGenerator(
std::vector<ScaleFunction> scales,
std::vector<double> scale_factors,
double max_scale_ratio
):
scales_(std::move(scales)),
scale_factors_(std::move(scale_factors)),
max_scale_ratio_{max_scale_ratio}
{
gen_descriptions();
}
/** \brief Adjust event parameters, adding scale variation
*
* The central renormalisation and factorisation scale of the returned
* event is given be the first base scale passed to the constructor.
* The scale variation (stored in event.variation()) is constructed as
* follows: For each base scale according to the arguments of the
* constructor we add one variation where both renormalisation and
* factorisation scale are set according to the current base scale.
* Then, all combinations where the base renormalisation and factorisation
* scales are multiplied by one of the scale factors are added.
* The case were both scales are multiplied by one is skipped.
* Scale combinations where the ratio is larger than the maximum scale ratio
* set in the constructor are likewise discarded.
*/
Event operator()(Event event) const;
private:
void gen_descriptions();
std::vector<ScaleFunction> scales_;
std::vector<double> scale_factors_;
std::vector<std::shared_ptr<ParameterDescription>> descriptions_;
- double max_scale_ratio_;
+ double max_scale_ratio_{};
};
} // namespace HEJ
diff --git a/src/HDF5Writer.cc b/src/HDF5Writer.cc
index 59251c7..d776808 100644
--- a/src/HDF5Writer.cc
+++ b/src/HDF5Writer.cc
@@ -1,426 +1,424 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/HDF5Writer.hh"
#include <cassert>
#include "LHEF/LHEF.h"
#include "HEJ/ConfigFlags.hh"
#ifdef HEJ_BUILD_WITH_HDF5
#include <cmath>
#include <cstddef>
#include <iterator>
#include <type_traits>
#include <utility>
#include <vector>
#include "highfive/H5File.hpp"
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#else
#include "HEJ/exceptions.hh"
#endif
#ifdef HEJ_BUILD_WITH_HDF5
namespace HEJ{
using HighFive::File;
using HighFive::DataSpace;
namespace {
using std::size_t;
constexpr size_t chunk_size = 1000;
constexpr unsigned compression_level = 3;
size_t to_index(event_type::EventType const type){
return type==0?0:std::floor(std::log2(static_cast<size_t>(type)))+1;
}
template<typename T>
void write_dataset(HighFive::Group & group, std::string const & name, T val) {
using data_t = std::decay_t<T>;
group.createDataSet<data_t>(name, DataSpace::From(val)).write(val);
}
template<typename T>
void write_dataset(
HighFive::Group & group, std::string const & name,
std::vector<T> const & val
) {
using data_t = std::decay_t<T>;
group.createDataSet<data_t>(name, DataSpace::From(val)).write(val);
}
struct Cache {
explicit Cache(size_t capacity):
capacity{capacity}
{
nparticles.reserve(capacity);
start.reserve(capacity);
pid.reserve(capacity);
weight.reserve(capacity);
scale.reserve(capacity);
fscale.reserve(capacity);
rscale.reserve(capacity);
aqed.reserve(capacity);
aqcd.reserve(capacity);
trials.reserve(capacity);
npLO.reserve(capacity);
npNLO.reserve(capacity);
}
void fill(Event const & ev) {
const auto hepeup = to_HEPEUP(ev, nullptr);
// HEJ event to get nice wrapper
const auto num_partons = std::distance(ev.cbegin_partons(),
ev.cend_partons());
assert(num_partons>0);
// Number of partons for LO matching, HEJ requires at least 2 partons
npLO.emplace_back(num_partons>1?num_partons-2:num_partons);
// Number of real emissions in NLO, HEJ is LO -> -1
npNLO.emplace_back(-1);
fill_event_params(hepeup);
fill_event_particles(hepeup);
}
void fill_event_params(LHEF::HEPEUP const & ev) {
nparticles.emplace_back(ev.NUP);
start.emplace_back(particle_pos);
pid.emplace_back(ev.IDPRUP);
weight.emplace_back(ev.XWGTUP);
scale.emplace_back(ev.SCALUP);
fscale.emplace_back(ev.scales.muf);
rscale.emplace_back(ev.scales.mur);
aqed.emplace_back(ev.AQEDUP);
aqcd.emplace_back(ev.AQCDUP);
// set first trial=1 for first event
// -> sum(trials) = 1 -> xs=sum(weights)/sum(trials) as in Sherpa
if(particle_pos == 0){
trials.emplace_back(1.);
} else {
trials.emplace_back(0.);
}
particle_pos += ev.NUP;
}
void fill_event_particles(LHEF::HEPEUP const & ev) {
id.insert(end(id), begin(ev.IDUP), end(ev.IDUP));
status.insert(end(status), begin(ev.ISTUP), end(ev.ISTUP));
lifetime.insert(end(lifetime), begin(ev.VTIMUP), end(ev.VTIMUP));
spin.insert(end(spin), begin(ev.SPINUP), end(ev.SPINUP));
for(int i = 0; i < ev.NUP; ++i) {
mother1.emplace_back(ev.MOTHUP[i].first);
mother2.emplace_back(ev.MOTHUP[i].second);
color1.emplace_back(ev.ICOLUP[i].first);
color2.emplace_back(ev.ICOLUP[i].second);
px.emplace_back(ev.PUP[i][0]);
py.emplace_back(ev.PUP[i][1]);
pz.emplace_back(ev.PUP[i][2]);
e.emplace_back(ev.PUP[i][3]);
m.emplace_back(ev.PUP[i][4]);
}
}
bool is_full() const {
return nparticles.size() >= capacity;
}
void clear() {
nparticles.clear();
start.clear();
pid.clear();
id.clear();
status.clear();
mother1.clear();
mother2.clear();
color1.clear();
color2.clear();
weight.clear();
scale.clear();
fscale.clear();
rscale.clear();
aqed.clear();
aqcd.clear();
trials.clear();
npLO.clear();
npNLO.clear();
px.clear();
py.clear();
pz.clear();
e.clear();
m.clear();
lifetime.clear();
spin.clear();
}
size_t capacity;
std::vector<int> nparticles, start, pid, id, status,
mother1, mother2, color1, color2, npLO, npNLO;
std::vector<double> weight, scale, fscale, rscale, aqed,
aqcd, trials, px, py, pz, e, m, lifetime, spin;
-
- private:
size_t particle_pos = 0;
};
} // namespace
struct HDF5Writer::HDF5WriterImpl{
File file;
LHEF::HEPRUP heprup;
Cache cache{chunk_size};
size_t event_idx = 0;
size_t particle_idx = 0;
HDF5WriterImpl(std::string const & filename, LHEF::HEPRUP && hepr):
file{filename, File::ReadWrite | File::Create | File::Truncate},
heprup{std::move(hepr)}
{
// TODO: code duplication with Les Houches Writer
const size_t max_number_types = to_index(event_type::last_type)+1;
heprup.NPRUP = max_number_types;
// ids of event types
heprup.LPRUP.clear();
heprup.LPRUP.reserve(max_number_types);
heprup.LPRUP.emplace_back(0);
for(size_t i=event_type::first_type+1; i<=event_type::last_type; i*=2) {
heprup.LPRUP.emplace_back(i);
}
heprup.XSECUP = std::vector<double>(max_number_types, 0.);
heprup.XERRUP = std::vector<double>(max_number_types, 0.);
heprup.XMAXUP = std::vector<double>(max_number_types, 0.);
write_init();
create_event_group();
create_particle_group();
}
void write_init() {
auto init = file.createGroup("init");
write_dataset(init, "PDFgroupA" , heprup.PDFGUP.first);
write_dataset(init, "PDFgroupB" , heprup.PDFGUP.second);
write_dataset(init, "PDFsetA" , heprup.PDFSUP.first);
write_dataset(init, "PDFsetB" , heprup.PDFSUP.second);
write_dataset(init, "beamA" , heprup.IDBMUP.first);
write_dataset(init, "beamB" , heprup.IDBMUP.second);
write_dataset(init, "energyA" , heprup.EBMUP.first);
write_dataset(init, "energyB" , heprup.EBMUP.second);
write_dataset(init, "numProcesses" , heprup.NPRUP);
write_dataset(init, "weightingStrategy", heprup.IDWTUP);
auto proc_info = file.createGroup("procInfo");
write_dataset(proc_info, "procId", heprup.LPRUP);
}
static HighFive::DataSetCreateProps const & hdf5_chunk() {
static const auto props = [](){
HighFive::DataSetCreateProps props;
props.add(HighFive::Chunking({chunk_size}));
props.add(HighFive::Deflate(compression_level));
return props;
}();
return props;
}
void create_event_group() {
static const auto dim = DataSpace({0}, {DataSpace::UNLIMITED});
auto events = file.createGroup("event");
events.createDataSet<int>("nparticles", dim, hdf5_chunk());
events.createDataSet<int>("start", dim, hdf5_chunk());
events.createDataSet<int>("pid", dim, hdf5_chunk());
events.createDataSet<double>("weight", dim, hdf5_chunk());
events.createDataSet<double>("scale", dim, hdf5_chunk());
events.createDataSet<double>("fscale", dim, hdf5_chunk());
events.createDataSet<double>("rscale", dim, hdf5_chunk());
events.createDataSet<double>("aqed", dim, hdf5_chunk());
events.createDataSet<double>("aqcd", dim, hdf5_chunk());
events.createDataSet<double>("trials", dim, hdf5_chunk());
events.createDataSet<int>("npLO", dim, hdf5_chunk());
events.createDataSet<int>("npNLO", dim, hdf5_chunk());
}
void resize_event_group(size_t new_size) {
auto events = file.getGroup("event");
events.getDataSet("nparticles").resize({new_size});
events.getDataSet("start").resize({new_size});
events.getDataSet("pid").resize({new_size});
events.getDataSet("weight").resize({new_size});
events.getDataSet("scale").resize({new_size});
events.getDataSet("fscale").resize({new_size});
events.getDataSet("rscale").resize({new_size});
events.getDataSet("aqed").resize({new_size});
events.getDataSet("aqcd").resize({new_size});
events.getDataSet("trials").resize({new_size});
events.getDataSet("npLO").resize({new_size});
events.getDataSet("npNLO").resize({new_size});
}
void create_particle_group() {
static const auto dim = DataSpace({0}, {DataSpace::UNLIMITED});
auto particles = file.createGroup("particle");
particles.createDataSet<int>("id", dim, hdf5_chunk());
particles.createDataSet<int>("status", dim, hdf5_chunk());
particles.createDataSet<int>("mother1", dim, hdf5_chunk());
particles.createDataSet<int>("mother2", dim, hdf5_chunk());
particles.createDataSet<int>("color1", dim, hdf5_chunk());
particles.createDataSet<int>("color2", dim, hdf5_chunk());
particles.createDataSet<double>("px", dim, hdf5_chunk());
particles.createDataSet<double>("py", dim, hdf5_chunk());
particles.createDataSet<double>("pz", dim, hdf5_chunk());
particles.createDataSet<double>("e", dim, hdf5_chunk());
particles.createDataSet<double>("m", dim, hdf5_chunk());
particles.createDataSet<double>("lifetime", dim, hdf5_chunk());
particles.createDataSet<double>("spin", dim, hdf5_chunk());
}
void resize_particle_group(size_t new_size) {
auto particles = file.getGroup("particle");
particles.getDataSet("id").resize({new_size});
particles.getDataSet("status").resize({new_size});
particles.getDataSet("mother1").resize({new_size});
particles.getDataSet("mother2").resize({new_size});
particles.getDataSet("color1").resize({new_size});
particles.getDataSet("color2").resize({new_size});
particles.getDataSet("px").resize({new_size});
particles.getDataSet("py").resize({new_size});
particles.getDataSet("pz").resize({new_size});
particles.getDataSet("e").resize({new_size});
particles.getDataSet("m").resize({new_size});
particles.getDataSet("lifetime").resize({new_size});
particles.getDataSet("spin").resize({new_size});
}
void write(Event const & ev){
cache.fill(ev);
if(cache.is_full()) {
dump_cache();
}
const double wt = ev.central().weight;
const size_t idx = to_index(ev.type());
heprup.XSECUP[idx] += wt;
heprup.XERRUP[idx] += wt*wt;
if(wt > heprup.XMAXUP[idx]){
heprup.XMAXUP[idx] = wt;
}
}
void dump_cache() {
write_event_params();
write_event_particles();
cache.clear();
}
void write_event_params() {
auto events = file.getGroup("event");
// choose arbitrary dataset to find size
const auto dataset = events.getDataSet("nparticles");
const size_t size = dataset.getSpace().getDimensions().front();
resize_event_group(size + cache.nparticles.size());
#define WRITE_FROM_CACHE(GROUP, PROPERTY) \
GROUP.getDataSet(#PROPERTY).select({size}, {cache.PROPERTY.size()}).write(cache.PROPERTY)
WRITE_FROM_CACHE(events, nparticles);
WRITE_FROM_CACHE(events, start);
WRITE_FROM_CACHE(events, pid);
WRITE_FROM_CACHE(events, weight);
WRITE_FROM_CACHE(events, scale);
WRITE_FROM_CACHE(events, fscale);
WRITE_FROM_CACHE(events, rscale);
WRITE_FROM_CACHE(events, aqed);
WRITE_FROM_CACHE(events, aqcd);
WRITE_FROM_CACHE(events, trials);
WRITE_FROM_CACHE(events, npLO);
WRITE_FROM_CACHE(events, npNLO);
}
void write_event_particles() {
auto particles = file.getGroup("particle");
// choose arbitrary dataset to find size
const auto dataset = particles.getDataSet("id");
const size_t size = dataset.getSpace().getDimensions().front();
resize_particle_group(size + cache.id.size());
WRITE_FROM_CACHE(particles, id);
WRITE_FROM_CACHE(particles, status);
WRITE_FROM_CACHE(particles, lifetime);
WRITE_FROM_CACHE(particles, spin);
WRITE_FROM_CACHE(particles, mother1);
WRITE_FROM_CACHE(particles, mother2);
WRITE_FROM_CACHE(particles, color1);
WRITE_FROM_CACHE(particles, color2);
WRITE_FROM_CACHE(particles, px);
WRITE_FROM_CACHE(particles, py);
WRITE_FROM_CACHE(particles, pz);
WRITE_FROM_CACHE(particles, e);
WRITE_FROM_CACHE(particles, m);
}
#undef WRITE_FROM_CACHE
void finish(){
dump_cache();
auto proc_info = file.getGroup("procInfo");
write_dataset(proc_info, "xSection", heprup.XSECUP);
write_dataset(proc_info, "error", heprup.XERRUP);
write_dataset(proc_info, "unitWeight", heprup.XMAXUP);
file.flush();
}
};
HDF5Writer::HDF5Writer(std::string const & filename, LHEF::HEPRUP heprup):
impl_{std::make_unique<HDF5Writer::HDF5WriterImpl>(
filename, std::move(heprup))}
{}
void HDF5Writer::write(Event const & ev){
impl_->write(ev);
}
void HDF5Writer::finish(){
impl_->finish();
}
} // namespace HEJ
#else // no HDF5 support
namespace HEJ{
class HDF5Writer::HDF5WriterImpl{};
HDF5Writer::HDF5Writer(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument{
"Failed to create HDF5 writer: "
"HEJ 2 was built without HDF5 support"
};
}
void HDF5Writer::write(Event const &){
assert(false);
}
void HDF5Writer::finish(){
assert(false);
}
}
#endif
namespace HEJ {
HDF5Writer::~HDF5Writer() = default;
}
diff --git a/src/Hjets.cc b/src/Hjets.cc
index 83d794a..a324cb0 100644
--- a/src/Hjets.cc
+++ b/src/Hjets.cc
@@ -1,1133 +1,1133 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Hjets.hh"
#include <cmath>
#include <complex>
#include <limits>
#include "CLHEP/Vector/LorentzVector.h"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/Constants.hh"
#ifdef HEJ_BUILD_WITH_QCDLOOP
#include <cassert>
#include "qcdloop/qcdloop.h"
#else
#include "HEJ/exceptions.hh"
#endif
namespace HEJ {
namespace currents {
namespace {
// short hand for math functions
using std::abs;
using std::conj;
using std::pow;
using std::sqrt;
const COM looprwfactor = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4);
constexpr double infinity = std::numeric_limits<double>::infinity();
// Loop integrals
#ifdef HEJ_BUILD_WITH_QCDLOOP
COM B0DD(HLV const & q, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_B0 = [](){
ql::Bubble<std::complex<double>,double,double> ql_B0;
ql_B0.setCacheSize(100);
return ql_B0;
}();
static std::vector<double> masses(2);
static std::vector<double> momenta(1);
for(auto & m: masses) m = mq*mq;
momenta.front() = q.m2();
ql_B0.integral(result, 1, masses, momenta);
return result[0];
}
COM C0DD(HLV const & q1, HLV const & q2, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_C0 = [](){
ql::Triangle<std::complex<double>,double,double> ql_C0;
ql_C0.setCacheSize(100);
return ql_C0;
}();
static std::vector<double> masses(3);
static std::vector<double> momenta(3);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = (q1+q2).m2();
ql_C0.integral(result, 1, masses, momenta);
return result[0];
}
COM D0DD(HLV const & q1, HLV const & q2, HLV q3, double mq)
{
static std::vector<std::complex<double>> result(3);
static auto ql_D0 = [](){
ql::Box<std::complex<double>,double,double> ql_D0;
ql_D0.setCacheSize(100);
return ql_D0;
}();
static std::vector<double> masses(4);
static std::vector<double> momenta(6);
for(auto & m: masses) m = mq*mq;
momenta[0] = q1.m2();
momenta[1] = q2.m2();
momenta[2] = q3.m2();
momenta[3] = (q1+q2+q3).m2();
momenta[4] = (q1+q2).m2();
momenta[5] = (q2+q3).m2();
ql_D0.integral(result, 1, masses, momenta);
return result[0];
}
COM A1(HLV const & q1, HLV const & q2, double mt)
// As given in Eq. (B.2) of VDD
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
assert(mt > 0.);
mt2=mt*mt;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 4.*mt2/Delta3*(Q2-q12-q22)
-1.-4.*q12*q22/Delta3-12.*q12*q22*Q2/Delta3/Delta3*(q12+q22-Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q2,mt)-B0DD(Q,mt) )
* ( 2.*q22/Delta3+12.*q12*q22/Delta3/Delta3*(q22-q12+Q2) )
- looprwfactor*COM(0,-1)*( B0DD(q1,mt)-B0DD(Q,mt) )
* ( 2.*q12/Delta3+12.*q12*q22/Delta3/Delta3*(q12-q22+Q2) )
- 2./Delta3/16/M_PI/M_PI*(q12+q22-Q2);
return ans;
}
COM A2(HLV const & q1, HLV const & q2, double mt)
// As given in Eq. (B.2) of VDD, but with high energy limit
// of invariants taken.
{
double q12,q22,Q2;
HLV Q;
double Delta3,mt2;
COM ans(COM(0.,0.));
assert(mt > 0.);
mt2=mt*mt;
q12=q1.m2();
q22=q2.m2();
Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD
Q2=Q.m2();
Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2;
ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*( 2.*mt2+1./2.*(q12+q22-Q2)
+2.*q12*q22*Q2/Delta3 )
+looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))
*q22*(q22-q12-Q2)/Delta3
+looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))
*q12*(q12-q22-Q2)/Delta3+1./16/M_PI/M_PI;
return ans;
}
#else // no QCDloop
COM A1(HLV const &, HLV const &, double) {
throw std::logic_error{"A1 called without QCDloop support"};
}
COM A2(HLV const & , HLV const & , double) {
throw std::logic_error{"A2 called without QCDloop support"};
}
#endif
void to_current(const HLV & q, current & ret){
ret[0]=q.e();
ret[1]=q.x();
ret[2]=q.y();
ret[3]=q.z();
}
/**
* @brief Higgs vertex contracted with current @param C1 and @param C2
*/
COM cHdot(const current & C1, const current & C2, const current & q1,
const current & q2, double mt, bool incBot, double mb, double vev)
{
if (mt == infinity) {
return (cdot(C1,C2)*cdot(q1,q2)-cdot(C1,q2)*cdot(C2,q1))/(3*M_PI*vev);
}
else {
HLV vq1,vq2;
vq1.set(q1[1].real(),q1[2].real(),q1[3].real(),q1[0].real());
vq2.set(q2[1].real(),q2[2].real(),q2[3].real(),q2[0].real());
// first minus sign obtained because of q1-difference to VDD
// Factor is because 4 mt^2 g^2/vev A1 -> 16 pi mt^2/vev alphas,
if(!(incBot))
return 16.*M_PI*mt*mt/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)
-cdot(C1,C2)*A2(-vq1,vq2,mt));
else
return 16.*M_PI*mt*mt/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)
-cdot(C1,C2)*A2(-vq1,vq2,mt))
+ 16.*M_PI*mb*mb/vev*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)
-cdot(C1,C2)*A2(-vq1,vq2,mb));
}
}
//@{
/**
* @brief Higgs+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @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 q1 t-channel momenta into higgs vertex
* @param q2 t-channel momenta out of higgs vertex
* @param mt top mass (inf or value)
* @param incBot Bool, to include bottom mass (true) or not (false)?
* @param mb bottom mass (value)
* @param pg Unordered Gluon momenta
*
* Calculates j^\mu H j_\mu. FKL with higgs vertex somewhere in the FKL chain.
* Handles all possible incoming states.
*/
double j_h_j(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
current j1p,j1m,j2p,j2m, q1v, q2v;
// Note need to flip helicities in anti-quark case.
joi(p1out, false, p1in, false, j1p);
joi(p1out, true, p1in, true, j1m);
joi(p2out, false, p2in, false, j2p);
joi(p2out, true, p2in, true, j2m);
to_current(q1, q1v);
to_current(q2, q2v);
COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb, vev);
COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb, vev);
COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb, vev);
COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb, vev);
// average over helicities
const double sst=(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm))/4.;
return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2());
}
} // namespace
double ME_H_qQ(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2, double mt,
bool incBot, double mb, double vev
){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qQbar(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qbarQ(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qbarQbar(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev);
}
double ME_H_qg(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/C_A;
}
double ME_H_qbarg(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/C_A;
}
double ME_H_gg(HLV const & p1out, HLV const & p1in, HLV const & p2out,
HLV const & p2in, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
return j_h_j(p1out, p1in, p2out, p2in, q1, q2, mt, incBot, mb, vev)
* K_g(p2out,p2in)/C_A * K_g(p1out,p1in)/C_A;
}
//@}
namespace {
//@{
/// @brief Higgs vertex contracted with one current
CCurrent jH(HLV const & pout, bool helout, HLV const & pin,
bool helin, HLV const & q1, HLV const & q2,
double mt, bool incBot, double mb, double vev
){
CCurrent j2 = joi(pout,helout,pin,helin);
CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz());
if(mt == infinity)
return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*vev);
else
{
if(incBot)
return (-16.*M_PI*mb*mb/vev*j2.dot(q1)*jq2*A1(-q1,q2,mb)
-16.*M_PI*mb*mb/vev*j2*A2(-q1,q2,mb))
+ (-16.*M_PI*mt*mt/vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)
-16.*M_PI*mt*mt/vev*j2*A2(-q1,q2,mt));
else
return (-16.*M_PI*mt*mt/vev*j2.dot(q1)*jq2*A1(-q1,q2,mt)
-16.*M_PI*mt*mt/vev*j2*A2(-q1,q2,mt));
}
}
//@}
//@{
/**
* @brief Higgs+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (W emission)
* @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 q1 t-channel momenta into higgs vertex
* @param q2 t-channel momenta out of higgs vertex
* @param mt top mass (inf or value)
* @param incBot Bool, to include bottom mass (true) or not (false)?
* @param mb bottom mass (value)
*
* Calculates j_{uno}^\mu H j_\mu. Unordered with higgs vertex
* somewhere in the FKL chain. Handles all possible incoming states.
*/
double juno_h_j(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in,
HLV const & qH1, HLV const & qH2,
double mt, bool incBot, double mb, double vev
){
// This construction is taking rapidity order: pg > p1out >> p2out
HLV q1=p1in-p1out; // Top End
HLV q2=-(p2in-p2out); // Bottom End
HLV qg=p1in-p1out-pg; // Extra bit post-gluon
// Note <p1|eps|pa> current split into two by gauge choice.
// See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
CCurrent mj1p=joi(p1out,false, p1in,false);
CCurrent mj1m=joi(p1out, true, p1in, true);
CCurrent jgap=joi(pg, false, p1in,false);
CCurrent jgam=joi(pg, true, p1in, true);
// Note for function joo(): <p1+|pg+> = <pg-|p1->.
CCurrent j2gp=joo(p1out, false, pg, false);
CCurrent j2gm=joo(p1out, true, pg, true);
CCurrent mjH2p=jH(p2out, true,p2in, true,qH1,qH2,mt,incBot,mb, vev);
CCurrent mjH2m=jH(p2out,false,p2in,false,qH1,qH2,mt,incBot,mb, vev);
// Dot products of these which occur again and again
COM MHmp=mj1m.dot(mjH2p);
COM MHmm=mj1m.dot(mjH2m);
COM MHpp=mj1p.dot(mjH2p);
COM MHpm=mj1p.dot(mjH2m);
CCurrent p2o(p2out), p2i(p2in), p1o(p1out), p1i(p1in), qsum(q1+qg);
CCurrent Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmm/2.) )/q1.m2();
CCurrent Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m + 2.*mj1m.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHmp/2.) )/q1.m2();
CCurrent Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2m
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpm/2.) )/q1.m2();
CCurrent Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p + 2.*mj1p.dot(pg)*mjH2p
+ ( p2o/pg.dot(p2out) + p2i/pg.dot(p2in) )*( qg.m2()*MHpp/2.) )/q1.m2();
CCurrent U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2();
CCurrent U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2();
CCurrent U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2();
CCurrent U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2();
CCurrent U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2();
CCurrent U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2();
CCurrent U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2();
CCurrent U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2();
constexpr double cf=C_F;
double amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
double amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2());
// Now add the t-channels for the Higgs
ampsq/=qH1.m2()*qg.m2();
ampsq/=16.;
// Factor of (Cf/Ca) for each quark to match ME_H_qQ.
ampsq*=C_F*C_F/C_A/C_A;
return ampsq;
}
} // namespace
double ME_H_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & qH1,
HLV const & qH2, double mt, bool incBot, double mb,
double vev
){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & qH1,
HLV const & qH2, double mt, bool incBot, double mb,
double vev
){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & qH1,
HLV const & qH2, double mt, bool incBot, double mb,
double vev
){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & qH1,
HLV const & qH2, double mt, bool incBot, double mb,
double vev
){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev);
}
double ME_H_unob_gQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & qH1,
HLV const & qH2, double mt, bool incBot, double mb,
double vev
){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev)
*K_g(p2out,p2in)/C_F;
}
double ME_H_unob_gQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & qH1,
HLV const & qH2, double mt, bool incBot, double mb,
double vev
){
return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, incBot, mb, vev)
*K_g(p2out,p2in)/C_F;
}
//@}
// Begin finite mass stuff
#ifdef HEJ_BUILD_WITH_QCDLOOP
namespace {
// All the stuff needed for the box functions in qg->qgH now...
COM E1(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2=-(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) +
S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) +
2.*(s34 + S1)*(s34 + S1)/Delta +
S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
2.*(s34 + S1)/Delta +
S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) -
C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) -
C0DD(k1 + k2, q2, mq)*2.*s34/
S1 - (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 +
S1)/(S1*Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(2.*s34*(s34 +
S1)*(S1 - S2)/(Delta*Sigma) +
2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S1)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM F1(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-S2*D0DD(k1, k2, q2,
mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) -
s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2,
mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) +
S2*pow((s34 + S2),2)/Delta/Delta)
- (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12)
- C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2,
mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 +
S2)*(S2 - S1)/(Delta*Sigma) + (B0DD(
k1 + k2, mq) -
B0DD(k1 + k2 + q2,
mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 +
S2)*(2.*s12*s34 -
S2*(S1 + S2))/(Delta*Sigma));
}
COM G1(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor*(S2*D0DD(k1, q2, k2,
mq)*(Delta/s12/s12 - 4.*mq*mq/s12) -
S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) -
C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) -
C0DD(k1, q2, mq))*4.*mq*mq/
s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./
s12 + (B0DD(k1 + q2, mq) -
B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1)));
}
COM E4(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k2, k1, q2,
mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) -
s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2,
k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S1),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta +
2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD(
q2, mq) - B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) +
s12*(S1 - S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))
*((2.*s12*(2.*s12*s34 - S1*(S1 + S2) + s34*(S2 - S1)))/(Delta*Sigma)));
}
COM F4(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (-s12*D0DD(k1, k2, q2,
mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) +
s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) +
s12*pow((s34 + S2),2)/Delta/Delta) -
C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*2.*(s34 +
S2)/Delta + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 -
S1*(S1 + S2) +
s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))
*(2.*s12*(2.*s12*s34 - S2*(S1 + S2) + s34*(S1 - S2))/(Delta*Sigma)));
}
COM G4(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2,
mq)*(Delta/s12 + (s12 + S1)/2. -
4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./
s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD(
k1 + q2, k2, mq) -
S2*C0DD(k2, q2, mq))*(1./
s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD(
k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*2./(s12 + S2));
}
COM E10(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta +
12.*mq*mq*S1*(s34 + S1)/Delta/Delta -
4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) -
s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) +
C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) -
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) +
8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s34*(4.*s12 + 3.*S1 +
S2)/(Delta*Sigma) +
8.*s12*s34*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2,
mq))*(12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*S1*(s34*(s12 + S2) -
S1*(s34 +
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM F10(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, Sigma, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
Sigma = 4.*s12*s34 - pow(S1+S2,2);
return looprwfactor* (s12*D0DD(k1, k2, q2,
mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta +
12.*mq*mq*s34*(s12 + S1)/Delta/Delta -
4.*s12*pow((s34 + S2),2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) -
s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) -
S2*C0DD(k2, q2, mq))*(1./Delta +
4.*mq*mq*S1/Delta/Delta -
4.*s12*(s34 + S2)/Delta/Delta -
4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) -
C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) +
4.*s12*s34*(S2 - S1)/(Delta*Sigma) +
4.*(s12 -
2.*mq*mq)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma)) - (B0DD(
k2 + q2, mq) -
B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) +
8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) -
B0DD(k1 + k2 + q2, mq) +
s12*C0DD(k1 + k2, q2,
mq))*(-12*s34*(2*s12 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) -
4.*s12*s34*s34/(S2*Delta*Delta) +
4.*s34*S1/(Delta*Sigma) -
4.*s34*(s12*s34*(2.*s12 + S2) -
S1*S1*(2.*s12 +
S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) -
B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 +
S2)*(2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma*Sigma) +
8.*s12*(2.*s34 + S1)/(Delta*Sigma) -
8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) +
s12*(S2 -
S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 -
S1*(S1 + S2))/(Delta*Sigma));
}
COM G10(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
HLV q2 = -(k1+k2+kh);
double Delta, S1, S2, s12, s34;
S1 = 2.*k1.dot(q2);
S2 = 2.*k2.dot(q2);
s12 = 2.*k1.dot(k2);
s34 = q2.m2();
Delta = s12*s34 - S1*S2;
return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. +
4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1,
k2 + q2, mq) -
S1*C0DD(k1, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2,
k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta +
4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) -
B0DD(k1 + q2, mq))*4.*(s34 +
S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) -
B0DD(k2 + q2, mq))*4.*s34/(Delta*S2));
}
COM H1(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq);
}
COM H4(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq);
}
COM H10(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq);
}
COM H2(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
return -1.*H1(k2,k1,kh,mq);
}
COM H5(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
return -1.*H4(k2,k1,kh,mq);
}
COM H12(HLV const & k1, HLV const & k2, HLV const & kh, double mq){
return -1.*H10(k2,k1,kh,mq);
}
// FL and FT functions
COM FL(HLV const & q1, HLV const & q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*((2.-
3.*q1.m2()*q2.dot(Q)/detQ2)*(B0DD(q1, mq) -
B0DD(Q, mq)) + (2. -
3.*q2.m2()*q1.dot(Q)/detQ2)*(B0DD(q2, mq) -
B0DD(Q, mq)) - (4.*mq*mq + q1.m2() + q2.m2() +
Q.m2() - 3.*q1.m2()*q2.m2()*Q.m2()/detQ2)*C0DD(
q1, q2, mq) - 2.);
}
COM FT(HLV const & q1, HLV const & q2, double mq){
HLV Q = q1 + q2;
double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2);
return -1./(2.*detQ2)*(Q.m2()*(B0DD(q1, mq) + B0DD(q2, mq) - 2.*B0DD(Q, mq) -
2.*q1.dot(q2)*C0DD(q1, q2, mq)) + (q1.m2() -
q2.m2()) *(B0DD(q1, mq) - B0DD(q2, mq))) -
q1.dot(q2)*FL(q1, q2, mq);
}
HLV ParityFlip(HLV const & p){
HLV flippedVector;
flippedVector.setE(p.e());
flippedVector.setX(-p.x());
flippedVector.setY(-p.y());
flippedVector.setZ(-p.z());
return flippedVector;
}
/// @brief HC amp for qg->qgH with finite top (i.e. j^{++}_H)
void g_gH_HC(HLV const & pa, HLV const & p1,
HLV pH, double mq, double vev, current &retAns)
{
current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1,
epsHapart2,conjepsH1part1,conjepsH1part2;
COM ang1a,sqa1;
const double F = 4.*mq*mq/vev;
// Easier to have the whole thing as current object so I can use cdot functionality.
// Means I need to write pa,p1 as current objects
to_current(pa, pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
//HEJ gauge
jio(pa,false,p1,false,cura1);
if(gluonforward){
// sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp)
ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp();
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp();
} else {
ang1a = sqrt(pa.minus()*p1.plus());
sqa1 = sqrt(pa.minus()*p1.plus());
}
const double prop = (pa-p1-pH).m2();
cmult(-1./sqrt(2.)/ang1a,cura1,conjeps1);
cmult(1./sqrt(2.)/sqa1,cura1,epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM h4 = H4(p1,-pa,pH,mq);
const COM h5 = H5(p1,-pa,pH,mq);
const COM h10 = H10(p1,-pa,pH,mq);
const COM h12 = H12(p1,-pa,pH,mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
const COM aH1 = cdot(pHcur, cura1);
current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/ang1a, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())
*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/ang1a, epsHa, T2);
}
cmult(sqrt(2.)/ang1a*aH1, epsHa, T3);
cmult(sqrt(2.)/sqa1*aH1, conjepsH1, T4);
cmult(-sqrt(2.)*Fta*pa.dot(p1)*aH1/sqa1, conjeps1, T5);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/ang1a, epsa, T6);
cmult(-aH1/sqrt(2.)/sqa1*h4*8.*COM(0.,1.)*M_PI*M_PI, conjeps1, T7);
cmult(aH1/sqrt(2.)/ang1a*h5*8.*COM(0.,1.)*M_PI*M_PI, epsa, T8);
cmult(aH1*aH1/2./ang1a/sqa1*h10*8.*COM(0.,1.)*M_PI*M_PI, pacur, T9);
cmult(-aH1*aH1/2./ang1a/sqa1*h12*8.*COM(0.,1.)*M_PI*M_PI, p1cur, T10);
current ans;
for(int i=0;i<4;++i)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5[i]+T6[i]+T7[i]+T8[i]+T9[i]+T10[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
/// @brief HNC amp for qg->qgH with finite top (i.e. j^{+-}_H)
void g_gH_HNC(HLV const & pa, HLV const & p1, HLV const & pH,
double mq, double vev, current &retAns
){
const double F = 4.*mq*mq/vev;
current conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur,
p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1,
conjepsH1part2;
// Find here if pa, meaning the gluon, is forward or backward
bool gluonforward = true;
if(pa.z() < 0)
gluonforward = false;
jio(pa,true,p1,true,cura1);
joi(p1,true,pa,true,cur1a);
to_current(pa,pacur);
to_current(p1,p1cur);
to_current(pH,pHcur);
to_current(pa+p1,paplusp1cur);
to_current(p1-pa,p1minuspacur);
const COM aH1 = cdot(pHcur,cura1);
const COM oneHa = conj(aH1); // = cdot(pHcur,cur1a)
COM const sqa1 = gluonforward
// sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp)
?sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp()
:sqrt(pa.minus()*p1.plus());
const double prop = (pa-p1-pH).m2();
cmult(1./sqrt(2.)/sqa1, cur1a, epsa);
cmult(-1./sqrt(2.)/sqa1, cura1, conjeps1);
const COM phase = cdot(conjeps1, epsa);
const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2();
const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2();
const COM Falpha = FT(p1-pa,pa-p1-pH,mq);
const COM Fbeta = FL(p1-pa,pa-p1-pH,mq);
const COM h1 = H1(p1,-pa, pH, mq);
const COM h2 = H2(p1,-pa, pH, mq);
const COM h4 = H4(p1,-pa, pH, mq);
const COM h5 = H5(p1,-pa, pH, mq);
const COM h10 = H10(p1,-pa, pH, mq);
const COM h12 = H12(p1,-pa, pH, mq);
cmult(Fta*pa.dot(pH), epsa, epsHapart1);
cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2);
cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1);
cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2);
cadd(epsHapart1, epsHapart2, epsHa);
cadd(conjepsH1part1, conjepsH1part2, conjepsH1);
current T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a,
T11b,T12a,T12b,T13;
if(gluonforward){
cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1);
cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/sqa1, epsHa, T2);
}
else{
cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, conjepsH1, T1);
cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())*((p1.x()+COM(0.,1.)*p1.y())/p1.perp())
*prop/sqa1, epsHa, T2);
}
const COM boxdiagFact = 8.*COM(0.,1.)*M_PI*M_PI;
cmult(aH1*sqrt(2.)/sqa1, epsHa, T3);
cmult(oneHa*sqrt(2.)/sqa1, conjepsH1, T4);
cmult(-2.*phase*Fta*pa.dot(pH), p1cur, T5a);
cmult(2.*phase*Ft1*p1.dot(pH), pacur, T5b);
cmult(-sqrt(2.)*Fta*p1.dot(pa)*oneHa/sqa1, conjeps1, T6);
cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/sqa1, epsa, T7);
cmult(-boxdiagFact*phase*h2, pacur, T8a);
cmult(boxdiagFact*phase*h1, p1cur, T8b);
cmult(boxdiagFact*aH1/sqrt(2.)/sqa1*h5, epsa, T9);
cmult(-boxdiagFact*oneHa/sqrt(2.)/sqa1*h4, conjeps1, T10);
cmult(boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h10, pacur, T11a);
cmult(-boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h12, p1cur, T11b);
cmult(-phase/(pa-p1).m2()*Falpha*(p1-pa).dot(pa-p1-pH), paplusp1cur, T12a);
cmult(phase/(pa-p1).m2()*Falpha*(pa+p1).dot(pa-p1-pH), p1minuspacur, T12b);
cmult(-phase*Fbeta*(pa-p1-pH).m2(), paplusp1cur, T13);
current ans;
for(int i=0;i<4;++i)
{
ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]
+T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i];
}
retAns[0] = F/prop*ans[0];
retAns[1] = F/prop*ans[1];
retAns[2] = F/prop*ans[2];
retAns[3] = F/prop*ans[3];
}
} // namespace anonymous
// JDC - new amplitude with Higgs emitted close to gluon with full mt effects.
// Keep usual HEJ-style function call
double ME_Houtside_gq(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in, HLV const & pH,
double mq, bool includeBottom, double mq2, double vev
){
current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip;
current retAns,retAnsb;
joi(p2out,true,p2in,true,cur2bplus);
joi(p2out,false,p2in,false,cur2bminus);
joi(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip);
joi(ParityFlip(p2out),false,ParityFlip(p2in),false,cur2bminusFlip);
COM app1,app2,apm1,apm2;
COM app3, app4, apm3, apm4;
if(!includeBottom)
{
g_gH_HC(p1in,p1out,pH,mq,vev,retAns);
app1=cdot(retAns,cur2bplus);
app2=cdot(retAns,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
app3=cdot(retAns,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,vev,retAns);
apm1=cdot(retAns,cur2bplus);
apm2=cdot(retAns,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
apm3=cdot(retAns,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip);
} else {
g_gH_HC(p1in,p1out,pH,mq,vev,retAns);
g_gH_HC(p1in,p1out,pH,mq2,vev,retAnsb);
app1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
app2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,vev,retAnsb);
app3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
app4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
// And non-conserving bits
g_gH_HNC(p1in,p1out,pH,mq,vev,retAns);
g_gH_HNC(p1in,p1out,pH,mq2,vev,retAnsb);
apm1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus);
apm2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,vev,retAns);
g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,vev,retAnsb);
apm3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip);
apm4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip);
}
return abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1)
+ abs2(apm2) + abs2(apm3) + abs2(apm4);
}
#endif // HEJ_BUILD_WITH_QCDLOOP
double C2gHgm(HLV const & p2, HLV const & p1, HLV const & pH, double vev){
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
- double s12,p1p,p2p;
+ double s12 = NAN,p1p = NAN,p2p = NAN;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p1p=p1.plus();
p2p=p2.plus();
} else { // opposite case
p1p=p1.minus();
p2p=p2.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=COM(0,1)*A/(2.*s12)*(p2p/p1p*conj(p1perp)*p3perp
+p1p/p2p*p1perp*conj(p3perp));
temp=temp*conj(temp);
return temp.real();
}
double C2gHgp(HLV const & p2, HLV const & p1, HLV const & pH, double vev){
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.23) in hep-ph/0301013
- double s12,php,p1p,phm;
+ double s12 = NAN,php = NAN,p1p = NAN,phm = NAN;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
php=pH.plus();
phm=pH.minus();
p1p=p1.plus();
} else { // opposite case
php=pH.minus();
phm=pH.plus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=-COM(0,1)*A/(2.*s12)*( conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p)
+s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm)
-pow(conj(p3perp)+(1.+php/p1p)*conj(p1perp),2)
/((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH)))) );
temp=temp*conj(temp);
return temp.real();
}
double C2qHqm(HLV const & p2, HLV const & p1, HLV const & pH, double vev){
const double A=1./(3.*M_PI*vev);
// Implements Eq. (4.21) in hep-ph/0301013
- double s12,p2p,p1p;
+ double s12 = NAN,p2p = NAN,p1p = NAN;
COM p1perp,p3perp,phperp;
// Determine first whether this is the case p1p\sim php>>p3p or the opposite
s12=p1.invariantMass2(-p2);
if (p2.pz()>0.) { // case considered in hep-ph/0301013
p2p=p2.plus();
p1p=p1.plus();
} else { // opposite case
p2p=p2.minus();
p1p=p1.minus();
}
p1perp=p1.px()+COM(0,1)*p1.py();
phperp=pH.px()+COM(0,1)*pH.py();
p3perp=-(p1perp+phperp);
COM temp=A/(2.*s12)*( sqrt(p2p/p1p)*p3perp*conj(p1perp)
+sqrt(p1p/p2p)*p1perp*conj(p3perp) );
temp=temp*conj(temp);
return temp.real();
}
} // namespace currents
} // namespace HEJ
diff --git a/src/JetSplitter.cc b/src/JetSplitter.cc
index c2ab802..4a1d034 100644
--- a/src/JetSplitter.cc
+++ b/src/JetSplitter.cc
@@ -1,191 +1,191 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/JetSplitter.hh"
#include <array>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <numeric>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/PseudoJet.hh"
#include "HEJ/Constants.hh"
#include "HEJ/RNG.hh"
#include "HEJ/exceptions.hh"
namespace HEJ {
namespace {
template<class Iterator>
bool same_pt_and_rapidity(
Iterator begin, Iterator end,
fastjet::PseudoJet const & jet
){
constexpr double ep = 1e-2;
const fastjet::PseudoJet reconstructed_jet = std::accumulate(
begin, end, fastjet::PseudoJet{}
);
return
(std::abs(reconstructed_jet.pt() - jet.pt()) < ep)
&& (std::abs(reconstructed_jet.rapidity() - jet.rapidity()) < ep)
;
}
bool all_in_one_jet(
std::vector<fastjet::PseudoJet> const & partons,
fastjet::JetDefinition jet_def, double min_jet_pt
){
fastjet::ClusterSequence ev(partons, jet_def);
const std::vector<fastjet::PseudoJet> testjet = ev.inclusive_jets(min_jet_pt);
return testjet.size() == 1u
&& testjet[0].constituents().size() == partons.size();
}
} // namespace
JetSplitter::JetSplitter(
fastjet::JetDefinition jet_def, double min_pt
):
min_jet_pt_{min_pt},
jet_def_{std::move(jet_def)}
{}
using SplitResult = JetSplitter::SplitResult;
SplitResult JetSplitter::split(
fastjet::PseudoJet const & j2split, int ncons, RNG & ran
) const{
if(ncons <= 0) {
throw std::invalid_argument{
"number of requested jet constituents less than 1"
};
}
double swt = 1.;
std::vector<fastjet::PseudoJet> jcons;
if(ncons == 1){
jcons.emplace_back(j2split);
jcons.back().set_user_index(0);
return {jcons, swt};
}
if(ncons == 2){
return Split2(j2split, ran);
}
const double R_max = R_factor*jet_def_.R();
assert(R_max < M_PI);
double pt_remaining = j2split.pt();
const double phi_jet = j2split.phi();
const double y_jet = j2split.rapidity();
for(int i = 0; i < ncons - 1; ++i){
/**
* Generate rapidity and azimuthal angle with a distance
* R = sqrt(delta_y^2 + delta_phi^2) < R_max
* from the jet centre
*/
const double R = R_max*ran.flat();
const double theta = 2*M_PI*ran.flat();
const double delta_phi = R*std::cos(theta);
const double delta_y = R*std::sin(theta);
/**
* Generate pt such that the total contribution of all partons
* along the jet pt axis does not exceed the jet pt
*/
const double pt_max = pt_remaining/std::cos(delta_phi);
assert(pt_max > 0);
if(pt_max < CMINPT) return {}; // no pt remaining for this parton
const double pt = (pt_max - CMINPT)*ran.flat() + CMINPT;
pt_remaining -= pt*std::cos(delta_phi);
jcons.emplace_back(
pt*std::cos(phi_jet + delta_phi), pt*std::sin(phi_jet + delta_phi),
pt*std::sinh(y_jet + delta_y), pt*std::cosh(y_jet + delta_y)
);
jcons.back().set_user_index(i);
swt *= 2*M_PI*R*R_max*pt*(pt_max - CMINPT);
}
const fastjet::PseudoJet p_total = std::accumulate(
jcons.cbegin(), jcons.cend(), fastjet::PseudoJet{}
);
// Calculate the pt of the last parton
const double last_px = j2split.px() - p_total.px();
const double last_py = j2split.py() - p_total.py();
const double last_pt = std::sqrt(last_px*last_px + last_py*last_py);
if(last_pt < CMINPT) return {};
// Calculate the rapidity of the last parton using the requirement that the
// new jet must have the same rapidity as the LO jet.
const double exp_2y_jet = (j2split.e() + j2split.pz())/(j2split.e() - j2split.pz());
const double bb = (p_total.e()+p_total.pz()) - exp_2y_jet*(p_total.e()-p_total.pz());
const double lasty = std::log(
( -bb+std::sqrt(bb*bb+4.*exp_2y_jet*last_pt*last_pt) ) / ( 2.*last_pt ));
jcons.emplace_back(
last_px, last_py, last_pt*std::sinh(lasty), last_pt*std::cosh(lasty)
);
jcons.back().set_user_index(ncons-1);
assert(same_pt_and_rapidity(jcons.cbegin(), jcons.cend(), j2split));
// Test that the last parton is not too far away from the jet centre.
if (jcons.back().delta_R(j2split) > R_max) return {};
if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
return {jcons, swt};
}
double JetSplitter::sample_distance_2p(double & wt, RNG & ran) const{
static constexpr double x_small = 0.1;
static constexpr double p_small = 0.4;
const double pR = ran.flat();
if(pR < p_small){
wt *= x_small/p_small;
return x_small/p_small*pR;
}
wt *= (1-x_small)/(1-p_small);
return (1-x_small)/(1-p_small)*(pR-p_small) + x_small;
}
SplitResult JetSplitter::Split2(
fastjet::PseudoJet const & j2split, RNG & ran
) const {
static constexpr std::size_t ncons = 2;
std::vector<fastjet::PseudoJet> jcons(ncons);
- std::array<double, ncons> R, phi, y, pt;
+ std::array<double, ncons> R{}, phi{}, y{}, pt{};
double wt = 1;
const double theta = 2*M_PI*ran.flat(); // angle in y-phi plane
// empiric observation: we are always within the jet radius
R[0] = sample_distance_2p(wt, ran)*jet_def_.R();
R[1] = -sample_distance_2p(wt, ran)*jet_def_.R();
for(std::size_t i = 0; i <= 1; ++i){
phi[i] = j2split.phi() + R[i]*std::cos(theta);
y[i] = j2split.rapidity() + R[i]*std::sin(theta);
}
for(std::size_t i = 0; i <= 1; ++i){
pt[i] = (j2split.py() - std::tan(phi[1-i])*j2split.px())/
(std::sin(phi[i]) - std::tan(phi[1-i])*std::cos(phi[i]));
if(pt[i] < CMINPT) return {};
jcons[i].reset_PtYPhiM(pt[i], y[i], phi[i]);
jcons[i].set_user_index(i);
}
assert(same_pt_and_rapidity(jcons.cbegin(), jcons.cend(), j2split));
if(! all_in_one_jet(jcons, jet_def_, min_jet_pt_)) return {};
wt *= 2*M_PI*pt[0]*R[0]*jet_def_.R()*jet_def_.R();
// from transformation of delta(R[1] - ...) to delta(pt[0] - ...)
const double dphi0 = phi[0] - j2split.phi();
const double ptJ = j2split.pt();
const double jacobian = std::cos(theta)*pt[1]*pt[1]/(ptJ*std::sin(dphi0));
return {jcons, jacobian*wt};
}
} // namespace HEJ
diff --git a/src/LorentzVector.cc b/src/LorentzVector.cc
index a4d9bd0..4b9727c 100644
--- a/src/LorentzVector.cc
+++ b/src/LorentzVector.cc
@@ -1,87 +1,87 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2020
* \copyright GPLv2 or later
*/
#include "HEJ/LorentzVector.hh"
#include <cmath>
#include <utility>
namespace HEJ {
namespace {
using COM = std::complex<double>;
COM perp(CLHEP::HepLorentzVector const & p) {
return COM{p.x(), p.y()};
}
void angle_aux(CLHEP::HepLorentzVector const & p, COM & perp_hat, double & p_plus, double & p_minus) {
const COM p_perp = perp(p);
if(p_perp == COM{0., 0.}){
perp_hat = -1.;
if(std::abs(std::abs(p.z())/p.t() - 1.) > 1e-9){
throw std::invalid_argument("zero pT particle with |pz| != E");
}
if(p.z() > 0){
p_plus = p.plus();
p_minus = 0.;
}else{
p_plus = 0.;
p_minus = p.minus();
}
}else{
perp_hat = p_perp/std::abs(p_perp);
p_plus = p.plus();
p_minus = p.minus();
}
}
} // anonymous namespace
// "angle" product angle(pi, pj) = <i j>
// see eq. (53) (\eqref{eq:angle_product}) in developer manual
COM angle(
CLHEP::HepLorentzVector const & pi,
CLHEP::HepLorentzVector const & pj
) {
COM pi_perp_hat, pj_perp_hat;
- double pi_plus, pi_minus, pj_plus, pj_minus;
+ double pi_plus = NAN, pi_minus = NAN, pj_plus = NAN, pj_minus = NAN;
angle_aux(pi, pi_perp_hat, pi_plus, pi_minus);
angle_aux(pj, pj_perp_hat, pj_plus, pj_minus);
return
+ std::sqrt(pi_minus*pj_plus)*pi_perp_hat
- std::sqrt(pi_plus*pj_minus)*pj_perp_hat;
}
// "square" product square(pi, pj) = [i j]
// see eq. (54) (\eqref{eq:square_product}) in developer manual
COM square(
CLHEP::HepLorentzVector const & pi,
CLHEP::HepLorentzVector const & pj
) {
return -std::conj(angle(pi, pj));
}
// see eq:P_massive, eq:P_massive_p, eq:P_massive_q, , eq:P_massive_plus,
// eq:P_massive_minus in the developer manual
std::pair<CLHEP::HepLorentzVector, CLHEP::HepLorentzVector>
split_into_lightlike(CLHEP::HepLorentzVector const & P) {
const double pt = P.perp();
if(P.plus() > 0){
const double y = std::log(P.plus()/pt);
const double E = P.m2()/(2.*P.plus());
return std::make_pair(
CLHEP::HepLorentzVector{P.x(), P.y(), pt*std::sinh(y), pt*std::cosh(y)},
CLHEP::HepLorentzVector{0., 0., -E, E}
);
} else {
const double y = std::log(pt/P.minus());
const double E = P.m2()/(2.*P.minus());
return std::make_pair(
CLHEP::HepLorentzVector{P.x(), P.y(), pt*std::sinh(y), pt*std::cosh(y)},
CLHEP::HepLorentzVector{0., 0., +E, E}
);
}
}
} // namespace HEJ
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index 739bcbb..6717179 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,2216 +1,2216 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <iterator>
#include <limits>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "fastjet/PseudoJet.hh"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/Constants.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/Event.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/Particle.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Zjets.hh"
#include "HEJ/event_types.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/jets.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*std::log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
return (
1. + alpha_s/(4.*M_PI)*beta0*std::log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
std::vector <Weights> virtual_part=virtual_corrections(event);
if(tree_kin_part.size() != virtual_part.size()) {
throw std::logic_error("tree and virtuals have different sizes");
}
Weights sum = Weights{0., std::vector<double>(event.variations().size(), 0.)};
for(size_t i=0; i<tree_kin_part.size(); ++i) {
sum += tree_kin_part.at(i)*virtual_part.at(i);
}
return tree_param(event)*sum;
}
Weights MatrixElement::tree(Event const & event) const {
std::vector <double> tree_kin_part=tree_kin(event);
double sum = 0.;
for(size_t i=0; i<tree_kin_part.size(); ++i) {
sum += tree_kin_part.at(i);
}
return tree_param(event)*sum;
}
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;
}
std::vector<Weights> MatrixElement::virtual_corrections(Event const & event) const {
if(! is_resummable(event.type())) {
return {Weights{0., std::vector<double>(event.variations().size(), 0.)}};
}
// only compute once for each renormalisation scale
std::unordered_map<double, std::vector<double> > known_vec;
std::vector<double> central_vec=virtual_corrections(event, event.central().mur);
known_vec.emplace(event.central().mur, central_vec);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
if(ME_it == end(known_vec)) {
known_vec.emplace(var.mur, virtual_corrections(event, var.mur));
}
}
// At this stage known_vec contains one vector of virtual corrections for each mur value
// Now put this into a vector of Weights
std::vector<Weights> result_vec;
for(size_t i=0; i<central_vec.size(); ++i) {
Weights result;
result.central = central_vec.at(i);
for(auto const & var: event.variations()) {
const auto ME_it = known_vec.find(var.mur);
result.variations.emplace_back(ME_it->second.at(i));
}
result_vec.emplace_back(result);
}
return result_vec;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
const 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;
std::size_t first_idx = 0;
std::size_t last_idx = partons.size() - 1;
#ifndef NDEBUG
bool wc = true;
#endif
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::unob) {
q -= partons[1].p;
++first_idx;
if (in[0].type != partons[1].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else if (event.type() == event_type::qqxexb) {
q -= partons[1].p;
++first_idx;
if (std::abs(partons[0].type) != std::abs(partons[1].type)){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
else {
if(event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
if (in[0].type != partons[0].type ){
q -= WBoson.p;
#ifndef NDEBUG
wc=false;
#endif
}
}
std::size_t first_idx_qqx = last_idx;
std::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(std::abs(backquark->type) != std::abs((backquark+1)->type)) {
wqq=true;
#ifndef NDEBUG
wc=false;
#endif
}
last_idx = std::distance(begin(partons), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(std::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(std::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;
}
#ifndef NDEBUG
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
);
#endif
return std::exp(exponent);
}
std::vector <double> MatrixElement::virtual_corrections_Z_qq(
Event const & event,
const double mur,
Particle const & ZBoson
) 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;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p;
fastjet::PseudoJet q_b = pa - partons[0].p;
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q_t -= partons[1].p;
q_b -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum_top=0., sum_bot=0., sum_mix=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
const double dy = partons[j+1].rapidity() - partons[j].rapidity();
const double tmp_top = omega0(alpha_s, mur, q_t)*dy;
const double tmp_bot = omega0(alpha_s, mur, q_b)*dy;
sum_top += tmp_top;
sum_bot += tmp_bot;
sum_mix += (tmp_top + tmp_bot) / 2.;
q_t -= partons[j+1].p;
q_b -= partons[j+1].p;
}
return {exp(sum_top), exp(sum_bot), exp(sum_mix)};
}
double MatrixElement::virtual_corrections_Z_qg(
Event const & event,
const double mur,
Particle const & ZBoson,
const bool is_gq_event
) 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;
#endif
assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{}));
assert(partons.size() >= 2);
assert(pa.pz() < pb.pz());
// If this is a gq event, don't subtract the Z momentum from first q
fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p);
size_t first_idx = 0;
size_t last_idx = partons.size() - 1;
// Unordered gluon does not contribute to the virtual corrections
if (event.type() == event_type::unob) {
// Gluon is partons[0] and is already subtracted
// partons[1] is the backward quark
q -= partons[1].p;
++first_idx;
} else if (event.type() == event_type::unof) {
// End sum at forward quark
--last_idx;
}
double sum=0.;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity()
- partons[j].rapidity());
q -= partons[j+1].p;
}
return exp(sum);
}
std::vector<double> MatrixElement::virtual_corrections(
Event const & event,
const 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) && std::abs(AWZH_boson->type) == pid::Wp){
return {virtual_corrections_W(event, mur, *AWZH_boson)};
}
if(AWZH_boson != end(out) && AWZH_boson->type == pid::Z_photon_mix){
if(is_gluon(in.back().type)){
// This is a qg event
return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, false)};
}else if(is_gluon(in.front().type)){
// This is a gq event
return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, true)};
}else{
// This is a qq event
return virtual_corrections_Z_qq(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;
std::size_t first_idx = 0;
std::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;
}
std::size_t first_idx_qqx = last_idx;
std::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(std::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(std::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 {std::exp(exponent)};
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ 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;
}
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, p1, 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,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, p1, p2);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, p1, p2);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda) {
return Cls;
} else {
return Cls - 4.0 / (kperp * kperp);
}
}
CLHEP::HepLorentzVector CLipatov(
CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
){
const CLHEP::HepLorentzVector p5 = qav-qbv;
const CLHEP::HepLorentzVector CL = -(qav+qbv)
+ 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;
}
//! 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
){
const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, pim, pip, pom, pop);
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,
const double lambda
) {
const double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()));
const double kperp=(qav-qbv).perp();
if (kperp>lambda)
return Cls;
return Cls-4./(kperp*kperp);
}
double C2Lipatov_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip,
CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop
) {
const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, pim, pip, pom, pop);
const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, pim, pip, pom, pop);
return -CL_t.dot(CL_b);
}
double C2Lipatovots_Mix(
CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t,
CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2,
const double lambda
) {
const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, pa, pb, p1, p2)
/ sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2());
const double kperp = (qav_t - qbv_t).perp();
if (kperp > lambda) {
return Cls;
} else {
return Cls - 4.0 / (kperp * kperp);
}
}
/** 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(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pg,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
using namespace currents;
assert(aptype!=pid::gluon); // aptype cannot be gluon
if (bptype==pid::gluon) {
if (is_quark(aptype))
return ME_unob_qg(pg,p1,pa,pn,pb);
else
return ME_unob_qbarg(pg,p1,pa,pn,pb);
}
else if (is_antiquark(bptype)) {
if (is_quark(aptype))
return ME_unob_qQbar(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQbar(pg,p1,pa,pn,pb);
}
else { //bptype == quark
if (is_quark(aptype))
return ME_unob_qQ(pg,p1,pa,pn,pb);
else
return ME_unob_qbarQ(pg,p1,pa,pn,pb);
}
throw std::logic_error("unreachable");
}
/** Matrix element squared for tree-level current-current scattering
* @param bptype Particle b PDG ID
* @param pgin Incoming gluon momentum
* @param pq Quark from splitting Momentum
* @param pqbar Anti-quark from splitting Momentum
* @param pn Particle n Momentum
* @param pb Particle b Momentum
* @param swap_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal.
* @returns ME Squared for Tree-Level Current-Current Scattering
*
* @note The qqxf contribution can be calculated by reversing the argument ordering.
*/
double ME_qqx_current(
ParticleID bptype,
CLHEP::HepLorentzVector const & pgin,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
bool const swap_q_qx
){
using namespace currents;
if (bptype==pid::gluon) {
if (swap_q_qx) // pq extremal
return ME_Exqqx_qqbarg(pgin,pq,pqbar,pn,pb);
else // pqbar extremal
return ME_Exqqx_qbarqg(pgin,pq,pqbar,pn,pb);
}
else { // b leg quark line
if (swap_q_qx) //extremal pq
return ME_Exqqx_qqbarQ(pgin,pq,pqbar,pn,pb);
else
return ME_Exqqx_qbarqQ(pgin,pq,pqbar,pn,pb);
}
throw std::logic_error("unreachable");
}
/* \brief Matrix element squared for central qqx tree-level current-current
* scattering
*
* @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
* @returns ME Squared for qqxmid Tree-Level Current-Current Scattering
*/
double ME_qqxmid_current(
ParticleID aptype, ParticleID bptype, int nabove,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons
){
using namespace currents;
// 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==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;
return wt*ME_Cenqqx_qq(pa, pb, partons, is_antiquark(bptype),
is_antiquark(aptype), swap_q_qx, nabove);
}
/** 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(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector const & pn,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & pa
){
using namespace currents;
if (aptype==pid::gluon && bptype==pid::gluon) {
return ME_gg(pn,pb,p1,pa);
} else if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_qg(pn,pb,p1,pa);
else
return ME_qbarg(pn,pb,p1,pa);
}
else if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_qg(p1,pa,pn,pb);
else
return ME_qbarg(p1,pa,pn,pb);
}
else { // they are both quark
if (is_quark(bptype)) {
if (is_quark(aptype))
return ME_qQ(pn,pb,p1,pa);
else
return ME_qQbar(pn,pb,p1,pa);
}
else {
if (is_quark(aptype))
return ME_qQbar(p1,pa,pn,pb);
else
return ME_qbarQbar(pn,pb,p1,pa);
}
}
throw std::logic_error("unreachable");
}
/** 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(
ParticleID aptype, ParticleID 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, ParticleProperties const & Wprop
){
using namespace currents;
// We know it cannot be gg incoming.
assert(!(aptype==pid::gluon && bptype==pid::gluon));
if (aptype==pid::gluon && bptype!=pid::gluon) {
if (is_quark(bptype))
return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop);
}
else if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop);
else
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 (is_quark(bptype)) {
if (is_quark(aptype))
return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
else {
if (is_quark(aptype))
return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop);
else
return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop);
}
}
else{ // emission off a, (first argument paout)
if (is_quark(aptype)) {
if (is_quark(bptype))
return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
else { // a is anti-quark
if (is_quark(bptype))
return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop);
else
return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop);
}
}
}
throw std::logic_error("unreachable");
}
/** 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(
ParticleID aptype, ParticleID 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, ParticleProperties const & Wprop
){
using namespace currents;
// we know they are not both gluons
assert(bptype != pid::gluon || aptype != pid::gluon);
if (bptype == pid::gluon && aptype != pid::gluon) {
// b gluon => W emission off a
if (is_quark(aptype))
return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop);
else
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 (is_quark(bptype)){
if (is_quark(aptype))
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,Wprop);
}
else{
if (is_quark(aptype))
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,Wprop);
}
}
else {// wc == false, emission off a, i.e. a is first current
if (is_quark(aptype)) {
if (is_quark(bptype)) //qq
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,Wprop);
}
else { // a is anti-quark
if (is_quark(bptype)) //qbarq
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,Wprop);
}
}
}
throw std::logic_error("unreachable");
}
/** \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 swap_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal.
* @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(
ParticleID aptype, ParticleID 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,
ParticleProperties const & Wprop
){
using namespace currents;
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
const double CFbackward = K_g( (swap_q_qx)?pq:pqbar ,pa)/C_F;
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==pid::gluon && bptype==pid::gluon) {
//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, Wprop)
* CFbackward;
else
return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop)
* CFbackward;
}
else {
assert(aptype==pid::gluon && bptype!=pid::gluon );
//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, Wprop)
* CFbackward;
else
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, is_antiquark(bptype), Wprop)
* CFbackward;
else
return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop)
* CFbackward;
}
}
throw std::logic_error("unreachable");
}
/* \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(
ParticleID aptype, ParticleID bptype,
int nabove, int nbelow,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & pq,
CLHEP::HepLorentzVector const & pqbar,
std::vector<CLHEP::HepLorentzVector> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
using namespace currents;
// 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==pid::gluon) wt*=K_g(partons.front(),pa)/C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F;
if(wqq)
return wt*ME_WCenqqx_qq(pa, pb, pl, plbar, partons,
is_antiquark(bptype),is_antiquark(aptype),
swap_q_qx, nabove, Wprop);
return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons,
is_antiquark(bptype), is_antiquark(aptype),
swap_q_qx, nabove, nbelow, wc, Wprop);
}
/** Matrix element squared for tree-level current-current scattering With Z+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 plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for Tree-Level Current-Current Scattering
*/
std::vector<double> ME_Z_current(
const ParticleID aptype, const ParticleID 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,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if(is_anyquark(aptype) && is_gluon(bptype)){
// This is a qg event
const double current_factor=ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
return { current_factor };
}else if(is_gluon(aptype) && is_anyquark(bptype)){
// This is a gq event
const double current_factor=ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw);
return { current_factor };
}else if(is_anyquark(aptype) && is_anyquark(bptype)){
// This is a qq event
return ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}else{
// It cannot be gg incoming
throw std::logic_error("Bad particle types");
}
}
/** Matrix element squared for backwards uno tree-level current-current
* scattering With Z+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 plbar Final state positron momentum
* @param pl Final state electron momentum
* @param Zprop Z properties
* @param stw2 Value of sin(theta_w)^2
* @param ctw Value of cos(theta_w)
* @returns ME Squared for unob Tree-Level Current-Current Scattering
*
* @note The unof contribution can be calculated by reversing the argument ordering.
*/
std::vector<double> ME_Z_uno_current(
const ParticleID aptype, const ParticleID 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,
ParticleProperties const & Zprop,
const double stw2, const double ctw
){
using namespace currents;
// we know they are not both gluons
assert(!is_gluon(aptype) || !is_gluon(bptype));
if (is_anyquark(aptype) && is_gluon(bptype)) {
// This is a qg event
const double current_factor=ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
return { current_factor };
} else if (is_gluon(aptype) && is_anyquark(bptype)) {
// This is a gq event
const double current_factor=ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw);
return { current_factor };
} else if (is_anyquark(aptype) && is_anyquark(bptype)) {
// This is a qq event
return ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw);
}
throw std::logic_error("unreachable");
}
/** \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(
ParticleID aptype, ParticleID 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
){
using namespace currents;
if (aptype==pid::gluon && bptype==pid::gluon)
// gg initial state
return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev);
else if (aptype==pid::gluon&&bptype!=pid::gluon) {
if (is_quark(bptype))
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==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
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 (is_quark(bptype)) {
if (is_quark(aptype))
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 (is_quark(aptype))
return ME_H_qbarQ(pn,pb,p1,pa,-qHp1,-qH,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("unreachable");
}
/** \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(
ParticleID aptype, ParticleID 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
){
using namespace currents;
if (bptype==pid::gluon && aptype!=pid::gluon) {
if (is_quark(aptype))
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 (is_quark(aptype)) {
if (is_quark(bptype))
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 (is_quark(bptype))
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("unreachable");
}
CLHEP::HepLorentzVector to_HepLorentzVector(Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(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
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
std::vector<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)};
case pid::Z_photon_mix:
return tree_kin_Z(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;
}
template<class InputIterator>
std::vector <double> FKL_ladder_weight_mix(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0_t, CLHEP::HepLorentzVector const & q0_b,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
const double lambda
){
double wt_top = 1;
double wt_bot = 1;
double wt_mix = 1;
auto qi_t = q0_t;
auto qi_b = q0_b;
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_t = qi_t - g;
const auto qip1_b = qi_b - g;
if(treat_as_extremal(*gluon_it)){
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, lambda)*C_A;
} else{
wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, p1, pn, lambda)*C_A;
wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A;
}
qi_t = qip1_t;
qi_b = qip1_b;
}
return {wt_top, wt_bot, wt_mix};
}
} // namespace
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;
}
auto const & 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(std::size_t i = 0; i < out_partons.size(); ++i){
assert(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 {
double tree_kin_jets_qqxmid(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
std::vector<Particle> const & partons,
double lambda
){
CLHEP::HepLorentzVector 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;
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);
}
const int nabove = std::distance(begin_ladder, backmidquark);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_qqxmid_current(
aptype, bptype, nabove, pa, pb,
pq, pqbar, partonsHLV
);
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*ladder_factor;
}
template<class InIter, class partIter>
double tree_kin_jets_qqx(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, double lambda){
const bool swap_q_qx = is_quark(*BeginPart);
const auto pgin = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-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));
assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon.
const double current_factor = ME_qqx_current(
(EndIn-1)->type, pgin, pq, pqbar, pn, pb, swap_q_qx
)/(4.*(N_C*N_C - 1.));
const double ladder_factor = FKL_ladder_weight(
(BeginPart+2), (EndPart-1),
pgin-pq-pqbar, pgin, pb, p1, pn, lambda
);
return current_factor*ladder_factor;
}
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;
}
} // namespace
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()==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()==event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
else if (ev.type()==event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
else if (ev.type()==event_type::extremal_qqxb){
return tree_kin_jets_qqx(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
else if (ev.type()==event_type::extremal_qqxf){
return tree_kin_jets_qqx(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
else if (ev.type()==event_type::central_qqx){
return tree_kin_jets_qqxmid(incoming[0].type, incoming[1].type,
to_HepLorentzVector(incoming[0]),
to_HepLorentzVector(incoming[1]),
partons, param_.regulator_lambda);
}
else {
throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets");
}
}
namespace{
double tree_kin_W_FKL(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector plbar, CLHEP::HepLorentzVector pl,
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, 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 CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
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, 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 CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
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, 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(
ParticleID aptype, ParticleID bptype,
CLHEP::HepLorentzVector pa,
CLHEP::HepLorentzVector pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector plbar, CLHEP::HepLorentzVector pl,
double lambda, ParticleProperties const & Wprop
){
CLHEP::HepLorentzVector 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.front());
auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqx emit W
bool wc = !wqq && (aptype==partons.front().type); //leg b emits w
assert(!wqq || !wc);
if(wqq){ // emission from qqx
qqxt -= pl + plbar;
} else if(!wc) { // emission from leg a
q0 -= pl + plbar;
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);
}
const int nabove = std::distance(begin_ladder, backmidquark);
const int nbelow = std::distance(begin_ladder_2, end_ladder);
std::vector<CLHEP::HepLorentzVector> partonsHLV;
partonsHLV.reserve(partons.size());
for (std::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, 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
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
#ifndef NDEBUG
// assert that there is exactly one decay corresponding to the W
assert(ev.decays().size() == 1);
auto const & w_boson{
std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(),
[] (Particle const & p) -> bool {
return std::abs(p.type) == ParticleID::Wp;
}) };
assert(w_boson != ev.outgoing().cend());
assert( static_cast<long int>(ev.decays().cbegin()->first)
== std::distance(ev.outgoing().cbegin(), w_boson) );
#endif
// find decay products of W
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) )
|| ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) );
// get lepton & neutrino
CLHEP::HepLorentzVector plbar, pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.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_.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_.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_.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_.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_.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_.ew_parameters.Wprop());
}
namespace{
std::vector <double> tree_kin_Z_FKL(
const ParticleID aptype, const ParticleID bptype,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
std::vector<Particle> const & partons,
CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw
){
const auto p1 = to_HepLorentzVector(partons[0]);
const auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
const auto begin_ladder = cbegin(partons) + 1;
const auto end_ladder = cend(partons) - 1;
const std::vector <double> current_factor = ME_Z_current(
aptype, bptype, pn, pb, p1, pa,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
}else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-p1;
ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q0 = pa-p1-plbar-pl;
const auto q1 = pa-p1;
ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
template<class InIter, class partIter>
std::vector <double> tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart,
const CLHEP::HepLorentzVector & plbar,
const CLHEP::HepLorentzVector & pl,
const double lambda, ParticleProperties const & Zprop,
const double stw2, const double ctw){
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));
const ParticleID aptype = (BeginIn)->type;
const ParticleID bptype = (BeginIn+1)->type;
const std::vector <double> current_factor = ME_Z_uno_current(
aptype, bptype, pn, pb, p1, pa, pg,
plbar, pl, Zprop, stw2, ctw
);
std::vector <double> ladder_factor;
if(is_gluon(bptype)){
// This is a qg event
const auto q0 = pa-pg-p1-plbar-pl;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else if(is_gluon(aptype)){
// This is a gq event
const auto q0 = pa-pg-p1;
ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1,
q0, pa, pb, p1, pn, lambda));
}else{
// This is a qq event
const auto q0 = pa-pg-p1-plbar-pl;
const auto q1 = pa-pg-p1;
ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1,
q0, q1, pa, pb, p1, pn, lambda);
}
std::vector <double> result;
for(size_t i=0; i<current_factor.size(); ++i){
result.push_back(current_factor.at(i)*ladder_factor.at(i));
}
return result;
}
} // namespace
std::vector<double> MatrixElement::tree_kin_Z(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
// find decay products of Z
auto const & decay{ ev.decays().cbegin()->second };
assert(decay.size() == 2);
assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0))
&& decay.at(0).type==-decay.at(1).type);
// get leptons
CLHEP::HepLorentzVector plbar, pl;
if (decay.at(0).type < 0){
plbar = to_HepLorentzVector(decay.at(0));
pl = to_HepLorentzVector(decay.at(1));
}
else{
pl = to_HepLorentzVector(decay.at(0));
plbar = to_HepLorentzVector(decay.at(1));
}
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto partons = tag_extremal_jet_partons(ev);
const double stw2 = param_.ew_parameters.sin2_tw();
const double ctw = param_.ew_parameters.cos_tw();
if(ev.type() == FKL){
return tree_kin_Z_FKL(incoming[0].type, incoming[1].type,
pa, pb, partons, plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_backward){
return tree_kin_Z_uno(cbegin(incoming), cbegin(partons),
cend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
if(ev.type() == unordered_forward){
return tree_kin_Z_uno(crbegin(incoming), crbegin(partons),
crend(partons), plbar, pl,
param_.regulator_lambda,
param_.ew_parameters.Zprop(),
stw2, ctw);
}
throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets");
}
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
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == pid::gluon) return currents::K_g(pout, pin);
return C_F;
}
#endif
// Colour factor in strict MRK limit
double K_MRK(ParticleID type) {
return (type == pid::gluon)?C_A:C_F;
}
} // namespace
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{
using namespace currents;
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,
CLHEP::HepLorentzVector const & qH,
CLHEP::HepLorentzVector const & 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
);
}
} // namespace
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;
+ double current_factor = NAN;
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 = C_A*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 = C_A*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:
case pid::Z_photon_mix:
return alpha_w*alpha_w;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
} // namespace
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*std::log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ
diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc
index d2615b6..d6aba12 100644
--- a/src/YAMLreader.cc
+++ b/src/YAMLreader.cc
@@ -1,565 +1,565 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \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/ConfigFlags.hh"
#include "HEJ/Constants.hh"
#include "HEJ/ScaleFunction.hh"
#include "HEJ/event_types.hh"
#include "HEJ/output_formats.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",
"scales", "scale factors", "max scale ratio", "import scales",
"log correction", "event output", "analysis", "analyses", "vev",
"regulator parameter", "max events"
};
// 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] = "";
}
for(auto && opt: {"FKL", "unordered", "extremal qqx", "central qqx", "non-resummable"}){
supported["event treatment"][opt] = "";
}
for(auto && particle_type: {"Higgs", "W", "Z"}){
for(auto && particle_opt: {"mass", "width"}){
supported["particle properties"][particle_type][particle_opt] = "";
}
}
for(auto && opt: {"type", "trials", "max deviation"}){
supported["unweight"][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;
}
WeightType to_weight_type(std::string const & setting){
if(setting == "weighted")
return WeightType::weighted;
if(setting =="resummation")
return WeightType::unweighted_resum;
if(setting =="partial")
return WeightType::partially_unweighted;
throw std::invalid_argument{"Unknown weight type \"" + setting + "\""};
}
} // namespace
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>());
}
void set_from_yaml(WeightType & setting, YAML::Node const & yaml){
setting = to_weight_type(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;
+ double R = NAN;
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;
}
ParticleProperties get_particle_properties(
YAML::Node const & node, std::string const & entry,
std::string const & boson
){
- ParticleProperties result;
+ ParticleProperties result{};
set_from_yaml(result.mass, node, entry, boson, "mass");
set_from_yaml(result.width, node, entry, boson, "width");
return result;
}
EWConstants get_ew_parameters(YAML::Node const & node){
EWConstants result;
- double vev;
+ double vev = NAN;
set_from_yaml(vev, node, "vev");
result.set_vevWZH(vev,
get_particle_properties(node, "particle properties", "W"),
get_particle_properties(node, "particle properties", "Z"),
get_particle_properties(node, "particle properties", "Higgs")
);
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},
{"HepMC2", FileFormat::HepMC2},
{"HepMC3", FileFormat::HepMC3},
{"HDF5", FileFormat::HDF5}
};
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;
if(suffix == "hepmc3") return FileFormat::HepMC3;
if(suffix == "hepmc2") return FileFormat::HepMC2;
if(suffix == "hdf5") return FileFormat::HDF5;
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 analyses sub-options
* those depend on the analysis being used and should be checked there
* similar for "import scales"
*/
if(name != "analyses" && 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;
+ std::size_t pos = 0;
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
) {
void * 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));
+ known.emplace(scale, reinterpret_cast<EventScale>(sym)); // NOLINT
}
}
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 & node, std::string const & entry
){
using namespace event_type;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{FKL, EventTreatment::discard},
{unob, EventTreatment::discard},
{unof, EventTreatment::discard},
{qqxexb, EventTreatment::discard},
{qqxexf, EventTreatment::discard},
{qqxmid, EventTreatment::discard},
{non_resummable, EventTreatment::discard}
};
set_from_yaml(treat.at(FKL), node, entry, "FKL");
set_from_yaml(treat.at(unob), node, entry, "unordered");
treat.at(unof) = treat.at(unob);
set_from_yaml(treat.at(qqxexb), node, entry, "extremal qqx");
treat.at(qqxexf) = treat.at(qqxexb);
set_from_yaml(treat.at(qqxmid), node, entry, "central qqx");
set_from_yaml(treat.at(non_resummable), node, entry, "non-resummable");
if(treat[non_resummable] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight non-resummable 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_if_defined(config.min_extparton_pt, yaml, "min extparton pt");
if(config.min_extparton_pt!=0)
std::cerr << "WARNING: \"min extparton pt\" is deprecated."
<< " Please use \"max ext soft pt fraction\" instead.\n";
set_from_yaml(
config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction"
);
// Sets the standard value, then changes this if defined
config.regulator_lambda=CLAMBDA;
set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter");
set_from_yaml_if_defined(config.max_events, yaml, "max events");
set_from_yaml(config.trials, yaml, "trials");
config.weight_type = WeightType::weighted;
set_from_yaml_if_defined(config.weight_type, yaml, "unweight", "type");
if(config.weight_type == WeightType::partially_unweighted) {
config.unweight_config = PartialUnweightConfig{};
set_from_yaml(
config.unweight_config->trials, yaml,
"unweight", "trials"
);
set_from_yaml(
config.unweight_config->max_dev, yaml,
"unweight", "max deviation"
);
}
else if(yaml["unweight"]) {
for(auto && opt: {"trials", "max deviation"}) {
if(yaml["unweight"][opt]) {
throw std::invalid_argument{
"'unweight: " + std::string{opt} + "' "
"is only supported if 'unweight: type' is set to 'partial'"
};
}
}
}
set_from_yaml(config.log_correction, yaml, "log correction");
config.treat = get_event_treatment(yaml, "event treatment");
set_from_yaml_if_defined(config.output, yaml, "event output");
config.rng = to_RNGConfig(yaml, "random generator");
set_from_yaml_if_defined(config.analyses_parameters, yaml, "analyses");
if(yaml["analysis"]){
std::cerr <<
"WARNING: Configuration entry 'analysis' is deprecated. "
" Use 'analyses' instead.\n";
set_from_yaml(config.analysis_parameters, yaml, "analysis");
if(!config.analysis_parameters.IsNull()){
config.analyses_parameters.push_back(config.analysis_parameters);
}
}
config.scales = to_ScaleConfig(yaml);
config.ew_parameters = get_ew_parameters(yaml);
config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling");
return config;
}
} // namespace
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
diff --git a/src/get_analysis.cc b/src/get_analysis.cc
index 8020d39..7b34046 100644
--- a/src/get_analysis.cc
+++ b/src/get_analysis.cc
@@ -1,51 +1,51 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/get_analysis.hh"
#include <dlfcn.h>
#include <string>
#include "yaml-cpp/yaml.h"
#include "HEJ/EmptyAnalysis.hh"
#include "HEJ/RivetAnalysis.hh"
namespace HEJ{
std::unique_ptr<Analysis> get_analysis(
YAML::Node const & parameters, LHEF::HEPRUP const & heprup
){
if(!parameters["plugin"]){
if(parameters["rivet"])
return RivetAnalysis::create(parameters, heprup);
return EmptyAnalysis::create(parameters, heprup);
}
using AnalysisMaker = std::unique_ptr<Analysis> (*)(
YAML::Node const &, LHEF::HEPRUP const &);
const auto plugin_name = parameters["plugin"].as<std::string>();
void * handle = dlopen(plugin_name.c_str(), RTLD_NOW);
char * error = dlerror();
if(error != nullptr) throw std::runtime_error(error);
void * sym = dlsym(handle, "make_analysis");
error = dlerror();
if(error != nullptr) throw std::runtime_error(error);
- auto make_analysis = reinterpret_cast<AnalysisMaker>(sym);
+ auto make_analysis = reinterpret_cast<AnalysisMaker>(sym); // NOLINT
return make_analysis(parameters, heprup);
}
std::vector<std::unique_ptr<Analysis>> get_analyses(
std::vector<YAML::Node> const & parameters, LHEF::HEPRUP const & heprup
){
std::vector<std::unique_ptr<Analysis>> anas;
for(auto const & param: parameters){
anas.emplace_back(get_analysis(param, heprup));
}
return anas;
}
} // namespace HEJ
diff --git a/src/jets.cc b/src/jets.cc
index bc40e40..05b2670 100644
--- a/src/jets.cc
+++ b/src/jets.cc
@@ -1,855 +1,855 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include <algorithm>
#include <cmath>
#include "HEJ/Constants.hh"
namespace {
// short hand for math functions
using std::abs;
using std::conj;
using std::pow;
using std::sqrt;
double metric(const size_t mu, const size_t nu) {
if(mu != nu) return 0.;
return (mu == 0)?1.:-1.;
}
} // Anonymous Namespace
namespace HEJ {
namespace currents {
// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
// @TODO: this is not a current and should be moved somewhere else
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
HLV const & pout,
HLV const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
CCurrent CCurrent::operator+(const CCurrent& other) const
{
COM result_c0=c0 + other.c0;
COM result_c1=c1 + other.c1;
COM result_c2=c2 + other.c2;
COM result_c3=c3 + other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator-(const CCurrent& other) const
{
COM result_c0=c0 - other.c0;
COM result_c1=c1 - other.c1;
COM result_c2=c2 - other.c2;
COM result_c3=c3 - other.c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const double x) const
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const double x) const
{
COM result_c0=CCurrent::c0/x;
COM result_c1=CCurrent::c1/x;
COM result_c2=CCurrent::c2/x;
COM result_c3=CCurrent::c3/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator*(const COM x) const
{
COM result_c0=x*CCurrent::c0;
COM result_c1=x*CCurrent::c1;
COM result_c2=x*CCurrent::c2;
COM result_c3=x*CCurrent::c3;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
CCurrent CCurrent::operator/(const COM x) const
{
COM result_c0=(CCurrent::c0)/x;
COM result_c1=(CCurrent::c1)/x;
COM result_c2=(CCurrent::c2)/x;
COM result_c3=(CCurrent::c3)/x;
return CCurrent(result_c0,result_c1,result_c2,result_c3);
}
std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
{
os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
return os;
}
CCurrent operator * ( double x, CCurrent const & m)
{
return m*x;
}
CCurrent operator * ( COM const & x, CCurrent const & m)
{
return m*x;
}
CCurrent operator / ( double x, CCurrent const & m)
{
return m/x;
}
CCurrent operator / ( COM const & x, CCurrent const & m)
{
return m/x;
}
COM CCurrent::dot(HLV const & p1) const
{
// Current goes (E,px,py,pz)
// Vector goes (px,py,pz,E)
return p1[3]*c0-p1[0]*c1-p1[1]*c2-p1[2]*c3;
}
COM CCurrent::dot(CCurrent const & p1) const
{
return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3;
}
//Current Functions
void joi(HLV const & pout, bool helout,
HLV const & pin, bool helin, current &cur
){
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
const double sqpop = sqrt(abs(pout.plus()));
const double sqpom = sqrt(abs(pout.minus()));
// Allow for "jii" format
const COM poperp = (pout.x()==0 && pout.y() ==0) ? -1 : pout.x()+COM(0,1)*pout.y();
if (helout != helin) {
throw std::invalid_argument{"Non-matching helicities"};
} else if (helout == false) { // negative helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(abs(pin.plus()));
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * poperp / abs(poperp);
cur[2] = -COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
const double sqpim = sqrt(abs(pin.minus()));
cur[0] = -sqpom * sqpim * poperp / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = COM(0,1) * cur[1];
cur[3] = -cur[0];
}
} else { // positive helicity
if (pin.plus() > pin.minus()) { // if forward
const double sqpip = sqrt(abs(pin.plus()));
cur[0] = sqpop * sqpip;
cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
cur[2] = COM(0,1) * cur[1];
cur[3] = cur[0];
} else { // if backward
double sqpim = sqrt(abs(pin.minus()));
cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
cur[1] = -sqpim * sqpop;
cur[2] = -COM(0,1) * cur[1];
cur[3] = -cur[0];
}
}
}
CCurrent joi(HLV const & pout, bool helout, HLV const & pin, bool helin)
{
current cur;
joi(pout, helout, pin, helin, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void jio(HLV const & pin, bool helin,
HLV const & pout, bool helout, current &cur
) {
joi(pout, !helout, pin, !helin, cur);
}
CCurrent jio(HLV const & pin, bool helin, HLV const & pout, bool helout){
current cur;
jio(pin, helin, pout, helout, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {
// Zero our current
cur[0] = 0.0;
cur[1] = 0.0;
cur[2] = 0.0;
cur[3] = 0.0;
if (heli!=helj) {
throw std::invalid_argument{"Non-matching helicities"};
} else if ( heli == true ) { // If positive helicity swap momenta
std::swap(pi,pj);
}
const double sqpjp = sqrt(abs(pj.plus() ));
const double sqpjm = sqrt(abs(pj.minus()));
const double sqpip = sqrt(abs(pi.plus() ));
const double sqpim = sqrt(abs(pi.minus()));
// Allow for "jii" format
const COM piperp = (pi.x()==0 && pi.y() ==0) ? -1 : pi.x()+COM(0,1)*pi.y();
const COM pjperp = (pj.x()==0 && pj.y() ==0) ? -1 : pj.x()+COM(0,1)*pj.y();
const COM phasei = piperp / abs(piperp);
const COM phasej = pjperp / abs(pjperp);
cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
}
CCurrent joo(HLV const & pi, bool heli, HLV const & pj, bool helj)
{
current cur;
joo(pi, heli, pj, helj, cur);
return CCurrent(cur[0],cur[1],cur[2],cur[3]);
}
namespace {
//@{
/**
* @brief Pure Jet FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1.
* @param p1in Incoming Particle 1.
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
*
* Calculates j_\mu j^\mu.
* Handles all possible incoming states. Helicity doesn't matter since we sum
* over all of them.
*/
double j_j(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
HLV const q1=p1in-p1out;
HLV const q2=-(p2in-p2out);
current mj1m,mj1p,mj2m,mj2p;
// Note need to flip helicities in anti-quark case.
joi(p1out, false, p1in, false, mj1p);
joi(p1out, true, p1in, true, mj1m);
joi(p2out, false, p2in, false, mj2p);
joi(p2out, true, p2in, true, mj2m);
COM const Mmp=cdot(mj1m,mj2p);
COM const Mmm=cdot(mj1m,mj2m);
COM const Mpp=cdot(mj1p,mj2p);
COM const Mpm=cdot(mj1p,mj2m);
double const sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);
// Multiply by Cf^2
return C_F*C_F*(sst)/(q1.m2()*q2.m2());
}
} // Anonymous Namespace
double ME_qQ(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qbarQbar(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in);
}
double ME_qg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/C_F;
}
double ME_qbarg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/(C_F);
}
double ME_gg(HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return j_j(p1out, p1in, p2out, p2in)*K_g(p1out, p1in)*K_g(p2out, p2in)/(C_F*C_F);
}
//@}
namespace {
double juno_j(HLV const & pg, HLV const & p1out,
HLV const & p1in, HLV const & p2out, HLV const & p2in
){
// This construction is taking rapidity order: pg > p1out >> p2out
HLV q1=p1in-p1out; // Top End
HLV q2=-(p2in-p2out); // Bottom End
HLV qg=p1in-p1out-pg; // Extra bit post-gluon
// Note <p1|eps|pa> current split into two by gauge choice.
// See James C's Thesis (p72). <p1|eps|pa> -> <p1|pg><pg|pa>
CCurrent mj1p=joi(p1out, false, p1in, false);
CCurrent mj1m=joi(p1out, true, p1in, true);
CCurrent jgap=joi(pg, false, p1in, false);
CCurrent jgam=joi(pg, true, p1in, true);
// Note for function joo(): <p1+|pg+> = <pg-|p1->.
CCurrent j2gp=joo(p1out, false, pg, false);
CCurrent j2gm=joo(p1out, true, pg, true);
CCurrent mj2p=joi(p2out, false, p2in, false);
CCurrent mj2m=joi(p2out, true, p2in, true);
// Dot products of these which occur again and again
COM Mmp=mj1m.dot(mj2p);
COM Mmm=mj1m.dot(mj2m);
COM Mpp=mj1p.dot(mj2p);
COM Mpm=mj1p.dot(mj2m);
CCurrent p1o(p1out),p2o(p2out),p2i(p2in),qsum(q1+qg),p1i(p1in);
CCurrent Lmm=(qsum*(Mmm)+(-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmm/2.))/q1.m2();
CCurrent Lmp=(qsum*(Mmp) + (-2.*mj2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2p
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mmp/2.))/q1.m2();
CCurrent Lpm=(qsum*(Mpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpm/2.))/q1.m2();
CCurrent Lpp=(qsum*(Mpp) + (-2.*mj2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2p
+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*Mpp/2.))/q1.m2();
CCurrent U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*Mmm)/(p1out+pg).m2();
CCurrent U1mp=(jgam.dot(mj2p)*j2gm+2.*p1o*Mmp)/(p1out+pg).m2();
CCurrent U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*Mpm)/(p1out+pg).m2();
CCurrent U1pp=(jgap.dot(mj2p)*j2gp+2.*p1o*Mpp)/(p1out+pg).m2();
CCurrent U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*Mmm)/(p1in-pg).m2();
CCurrent U2mp=((-1.)*j2gm.dot(mj2p)*jgam+2.*p1i*Mmp)/(p1in-pg).m2();
CCurrent U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*Mpm)/(p1in-pg).m2();
CCurrent U2pp=((-1.)*j2gp.dot(mj2p)*jgap+2.*p1i*Mpp)/(p1in-pg).m2();
constexpr double cf=C_F;
double amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
double amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp);
double ampsq=-(amm+amp+apm+app);
//Divide by t-channels
ampsq/=q2.m2()*qg.m2();
ampsq/=16.;
// Factor of (Cf/Ca) for each quark to match j_j.
ampsq*=(C_F*C_F)/(C_A*C_A);
return ampsq;
}
} // Anonymous Namespace
//Unordered bits for pure jet
double ME_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in);
}
double ME_unob_qg( HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F;
}
double ME_unob_qbarg(HLV const & pg, HLV const & p1out, HLV const & p1in,
HLV const & p2out, HLV const & p2in
){
return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/C_F;
}
namespace {
void eps(HLV const & refmom, HLV const & kb, bool hel, current &ep){
//Positive helicity eps has negative helicity choices for spinors and vice versa
joi(refmom,hel,kb,hel,ep);
- double norm;
+ double norm = NAN;
if(kb.z()<0.) norm = sqrt(2.*refmom.plus()*kb.minus());
else norm = sqrt(2.*refmom.minus()*kb.plus());
// Normalise
std::for_each(begin(ep), end(ep), [&,norm](auto & val){val/=norm;});
}
COM qggm1(HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2,
HLV const & p3, bool helchain, bool heltop, bool helb, HLV const & refmom
){
// Since everything is defined with currents, need to use compeleness relation
// to expand p slash. i.e. pslash = |p><p|. Only one helicity 'survives' as
// defined by the helicities of the spinors at the end of the chain.
current cur33, cur23, curb3, cur2b, cur1a, ep;
joo(p3, helchain, p3, helchain,cur33);
joo(p2,helchain,p3,helchain,cur23);
jio(pb,helchain,p3,helchain,curb3);
joi(p2,helchain,pb,helchain,cur2b);
joi(p1, heltop, pa, heltop,cur1a);
const double t2 = (p3-pb)*(p3-pb);
//Calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
COM v1[4][4];
for(int u=0; u<4;++u){
for(int v=0; v<4;++v){
v1[u][v]=(cur23[u]*cur33[v]-cur2b[u]*curb3[v])/t2*(-1.);
}
}
//Dot in current and eps
//eps
eps(refmom,pb,helb, ep);
COM M1=0.;
// Perform Contraction: g^{ik} j_{1a}_k * v1_i^j eps^l g_lj
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M1+= metric(i,i) *cur1a[i]*(v1[i][j])*ep[j]*metric(j,j);
}
}
return M1;
}
COM qggm2(HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2,
HLV const & p3, bool helchain, bool heltop, bool helb, HLV const & refmom
){
// Since everything is defined with currents, need to use completeness relation
// to expand p slash. i.e. pslash = |p><p|. Only one helicity 'survives' as
// defined by the helicities of the spinors at the end of the chain.
current cur22, cur23, curb3, cur2b, cur1a, ep;
joo(p2, helchain, p2, helchain, cur22);
joo(p2, helchain, p3, helchain, cur23);
jio(pb, helchain, p3, helchain, curb3);
joi(p2, helchain, pb, helchain, cur2b);
joi(p1, heltop, pa, heltop, cur1a);
const double t2t = (p2-pb)*(p2-pb);
//Calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
COM v2[4][4]={};
for(int u=0; u<4;++u){
for(int v=0; v<4; ++v){
v2[u][v]=(cur22[v]*cur23[u]-cur2b[v]*curb3[u])/t2t;
}
}
//Dot in current and eps
//eps
eps(refmom,pb,helb, ep);
COM M2=0.;
// Perform Contraction: g^{ik} j_{1a}_k * v2_i^j eps^l g_lj
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M2+= metric(i,i)*cur1a[i]*(v2[i][j])*ep[j]*metric(j,j);
}
}
return M2;
}
COM qggm3(HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2,
HLV const & p3, bool helchain, bool heltop, bool helb, HLV const & refmom
){
current qqcur,ep,cur1a;
const double s23 = (p2+p3)*(p2+p3);
joo(p2,helchain,p3,helchain,qqcur);
joi(p1, heltop, pa, heltop,cur1a);
//Redefine relevant momenta as currents - for ease of calling correct part of vector
const current kb{pb.e(), pb.x(), pb.y(), pb.z()};
const current k2{p2.e(), p2.x(), p2.y(), p2.z()};
const current k3{p3.e(), p3.x(), p3.y(), p3.z()};
//Calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
COM V3g[4][4]={};
const COM kbqq=kb[0]*qqcur[0] -kb[1]*qqcur[1] -kb[2]*qqcur[2] -kb[3]*qqcur[3];
for(int u=0;u<4;++u){
for(int v=0;v<4;++v){
V3g[u][v] += 2.*COM(0.,1.)*(((k2[v]+k3[v])*qqcur[u] - (kb[u])*qqcur[v])+
kbqq*metric(u,v))/s23;
}
}
eps(refmom,pb,helb, ep);
COM M3=0.;
// Perform Contraction: g^{ik} j_{1a}_k * (v2_i^j) eps^l g_lj
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M3+= metric(i,i)*cur1a[i]*(V3g[i][j])*ep[j]*metric(j,j);
}
}
return M3;
}
//Now the function to give helicity/colour sum/average
double MqgtqQQ(HLV const & pa, HLV const & pb,
HLV const & p1, HLV const & p2, HLV const & p3
){
// 4 indepedent helicity choices (complex conjugation related).
//Need to evalute each independent hel configuration and store that result somewhere
const COM Mmmm1 = qggm1(pa,pb,p1,p2,p3,false,false,false, pa);
const COM Mmmm2 = qggm2(pa,pb,p1,p2,p3,false,false,false, pa);
const COM Mmmm3 = qggm3(pa,pb,p1,p2,p3,false,false,false, pa);
const COM Mmmp1 = qggm1(pa,pb,p1,p2,p3,false,true, false, pa);
const COM Mmmp2 = qggm2(pa,pb,p1,p2,p3,false,true, false, pa);
const COM Mmmp3 = qggm3(pa,pb,p1,p2,p3,false,true, false, pa);
const COM Mpmm1 = qggm1(pa,pb,p1,p2,p3,false,false,true, pa);
const COM Mpmm2 = qggm2(pa,pb,p1,p2,p3,false,false,true, pa);
const COM Mpmm3 = qggm3(pa,pb,p1,p2,p3,false,false,true, pa);
const COM Mpmp1 = qggm1(pa,pb,p1,p2,p3,false,true, true, pa);
const COM Mpmp2 = qggm2(pa,pb,p1,p2,p3,false,true, true, pa);
const COM Mpmp3 = qggm3(pa,pb,p1,p2,p3,false,true, true, pa);
//Colour factors:
const COM cm1m1 = 8./3.;
const COM cm2m2 = 8./3.;
const COM cm3m3 = 6.;
const COM cm1m2 = -1./3.;
const COM cm1m3 = -3.*COM(0.,1.);
const COM cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
const 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)));
const double Mmmp = real(cm1m1*pow(abs(Mmmp1),2)+cm2m2*pow(abs(Mmmp2),2)+
cm3m3*pow(abs(Mmmp3),2)+2.*real(cm1m2*Mmmp1*conj(Mmmp2))+
2.*real(cm1m3*Mmmp1*conj(Mmmp3))+2.*real(cm2m3*Mmmp2*conj(Mmmp3)));
const 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)));
const double Mpmp = real(cm1m1*pow(abs(Mpmp1),2)+cm2m2*pow(abs(Mpmp2),2)+
cm3m3*pow(abs(Mpmp3),2)+2.*real(cm1m2*Mpmp1*conj(Mpmp2))+
2.*real(cm1m3*Mpmp1*conj(Mpmp3))+2.*real(cm2m3*Mpmp2*conj(Mpmp3)));
// Factor of 2 for the helicity for conjugate configurations
return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/3.)/(pa-p1).m2()/(p2+p3-pb).m2();
}
} // Anonymous Namespace
// Extremal qqx
double ME_Exqqx_qbarqQ(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout);
}
double ME_Exqqx_qqbarQ(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout);
}
double ME_Exqqx_qbarqg(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqout, pqbarout)*K_g(p2out,p2in)/C_F;
}
double ME_Exqqx_qqbarg(HLV const & pgin, HLV const & pqout,
HLV const & pqbarout, HLV const & p2out, HLV const & p2in
){
return MqgtqQQ(p2in, pgin, p2out, pqbarout, pqout)*K_g(p2out,p2in)/C_F;
}
namespace {
void CurrentMatrix(current const & j1, current const & j2,
COM array[4][4]
){
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
array[i][j]=j1[i]*j2[j];
}
}
}
//qqbar produced in the middle
COM m1(current const & jtop, current const & jbot, bool hel2,
HLV const & ka, HLV const & kb, HLV const & k2, HLV const & k3,
std::vector<HLV> const & partons, unsigned int nabove
){
const double s23 = 2.*(k2*k3);
const double sa2 = 2.*(ka*k2);
const double sa3 = 2.*(ka*k3);
const double s12 = 2.*(partons.front()*k2);
const double s13 = 2.*(partons.front()*k3);
const double sb2 = 2.*(kb*k2);
const double sb3 = 2.*(kb*k3);
const double s42 = 2.*(partons.back()*k2);
const double s43 = 2.*(partons.back()*k3);
HLV q1=ka-partons.front();
for(unsigned int i=1;i<nabove+1;++i)
q1-=partons.at(i);
const HLV q2=q1-partons.at(nabove+2)-partons.at(nabove+1);
const double t1 = q1.m2();
const double t3 = q2.m2();
current cur23;
joo(k2,hel2,k3,hel2,cur23);
const current curka{ka.e(),ka.px(),ka.py(),ka.pz()};
const current curkb{kb.e(),kb.px(),kb.py(),kb.pz()};
const current curk1{partons.front().e(),partons.front().px(),
partons.front().py(),partons.front().pz()};
const current curk2{k2.e(),k2.px(),k2.py(),k2.pz()};
const current curk3{k3.e(),k3.px(),k3.py(),k3.pz()};
const current curk4{partons.back().e(),partons.back().px(),
partons.back().py(),partons.back().pz()};
const current qc1{q1.e(),q1.px(),q1.py(),q1.pz()};
const current qc2{q2.e(),q2.px(),q2.py(),q2.pz()};
//Create the two bits of this vertex
COM veik[4][4],v3g[4][4];
for(int i=0;i<4;++i) {
for(int j=0;j<4;++j){
veik[i][j] = (cdot(cur23,curka)*(t1/(sa2+sa3))+cdot(cur23,curk1)*
(t1/(s12+s13))-cdot(cur23,curkb)*(t3/(sb2+sb3))-
cdot(cur23,curk4)*(t3/(s42+s43)))*metric(i,j);
}
}
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
v3g[i][j] = qc1[j]*cur23[i]+curk2[j]*cur23[i]+curk3[j]*cur23[i]+
qc2[i]*cur23[j]-curk2[i]*cur23[j]-curk3[i]*cur23[j]-
(cdot(qc1,cur23)+cdot(qc2,cur23))*metric(i,j);
}
}
//Now dot in the currents - potential problem here with Lorentz
//indicies, so check this
COM M1=0;
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M1+= metric(i,i)*jtop[i]*(veik[i][j]+v3g[i][j])*jbot[j]*metric(j,j);
}
}
M1/=s23;
return M1;
}
COM m2(current const & jtop, current const & jbot, bool hel2,
HLV const & ka, HLV k2, HLV const & k3,
std::vector<HLV> const & partons, unsigned int nabove
){
//In order to get correct momentum dependence in the vertex, forst
//have to work with CCurrent objects and then convert to 'current'
current cur22,cur23,cur2q,curq3;
COM qarray[4][4]={};
COM temp[4][4]={};
joo(k2,hel2,k2,hel2,cur22);
joo(k2,hel2,k3,hel2,cur23);
joi(k2,hel2,ka,hel2,cur2q);
jio(ka,hel2,k3,hel2,curq3);
CurrentMatrix(cur2q, curq3, qarray);
for(unsigned int i =0; i<nabove+1; ++i){
joo(k2,hel2,partons.at(i),hel2,cur2q);
joo(partons.at(i),hel2,k3,hel2,curq3);
CurrentMatrix(cur2q, curq3, temp);
for(int ii=0;ii<4;++ii){
for(int jj=0;jj<4;++jj){
qarray[ii][jj]=qarray[ii][jj]-temp[ii][jj];
}
}
}
HLV qt=ka-k2;
for(unsigned int i=0; i<nabove+1;++i){
qt-=partons.at(i);
}
const double t2=qt*qt;
COM tempv[4][4];
for(int i=0; i<4;++i){
for(int j=0;j<4;++j){
tempv[i][j] = COM(0.,1.)*(qarray[i][j]-cur22[i]*cur23[j]);
}
}
COM M2=0.;
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M2+= metric(i,i)*jtop[i]*(tempv[i][j])*jbot[j]*metric(j,j);
}
}
M2/=t2;
return M2;
}
COM m3(current const & jtop, current const & jbot, bool hel2,
HLV const & ka, HLV const & k2, HLV const & k3,
std::vector<HLV> const & partons, unsigned int nabove
){
COM M3=0.;
current cur23,cur33,cur2q,curq3;
COM qarray[4][4]={};
COM temp[4][4]={};
joo(k3,hel2,k3,hel2,cur33);
joo(k2,hel2,k3,hel2,cur23);
joi(k2,hel2,ka,hel2,cur2q);
jio(ka,hel2,k3,hel2,curq3);
CurrentMatrix(cur2q, curq3, qarray);
for(unsigned int i =0; i<nabove+1; ++i){
joo(k2,hel2,partons.at(i),hel2,cur2q);
joo(partons.at(i),hel2,k3,hel2,curq3);
CurrentMatrix(cur2q, curq3, temp);
for(int ii=0;ii<4;++ii){
for(int jj=0;jj<4;++jj){
qarray[ii][jj]=qarray[ii][jj]-temp[ii][jj];
}
}
}
HLV qt=ka-k3;
for(unsigned int i=0; i<nabove+1;++i){
qt-=partons.at(i);
}
const double t2t=qt*qt;
COM tempv[4][4];
for(int i=0; i<4;++i){
for(int j=0;j<4;++j){
tempv[i][j] = COM(0.,-1.)*(qarray[j][i]-cur23[j]*cur33[i]);
}
}
for(int i=0;i<4;++i){
for(int j=0;j<4;++j){
M3+= metric(i,i)*jtop[i]*(tempv[i][j])*jbot[j]*metric(j,j);
}
}
M3/= t2t;
return M3;
}
} // Anonymous Namespace
double ME_Cenqqx_qq(HLV const & ka, HLV const & kb,
std::vector<HLV> const & partons, bool aqlinepa,
bool aqlinepb, bool qqxmarker, int nabove
){
//Get all the possible outer currents
current j1p,j1m,j4p,j4m;
if(!(aqlinepa)){
joi(partons.front(),true,ka,true,j1p);
joi(partons.front(),false,ka,false,j1m);
}
if(aqlinepa){
jio(ka,true,partons.front(),true,j1p);
jio(ka,false,partons.front(),false,j1m);
}
if(!(aqlinepb)){
joi(partons.back(),true,kb,true,j4p);
joi(partons.back(),false,kb,false,j4m);
}
if(aqlinepb){
jio(kb,true,partons.back(),true,j4p);
jio(kb,false,partons.back(),false,j4m);
}
HLV k2,k3;
if(!(qqxmarker)){
k2=partons.at(nabove+1);
k3=partons.at(nabove+2);
}
else{
k2=partons.at(nabove+2);
k3=partons.at(nabove+1);
}
//8 helicity choices we can make, but only 4 indepedent ones
//(complex conjugation related).
const COM Mmmm1 = m1(j1m,j4m,false,ka,kb,k2,k3,partons,nabove);
const COM Mmmm2 = m2(j1m,j4m,false,ka, k2,k3,partons,nabove);
const COM Mmmm3 = m3(j1m,j4m,false,ka, k2,k3,partons,nabove);
const COM Mmmp1 = m1(j1m,j4m,true, ka,kb,k2,k3,partons,nabove);
const COM Mmmp2 = m2(j1m,j4m,true, ka, k2,k3,partons,nabove);
const COM Mmmp3 = m3(j1m,j4m,true, ka, k2,k3,partons,nabove);
const COM Mpmm1 = m1(j1p,j4m,false,ka,kb,k2,k3,partons,nabove);
const COM Mpmm2 = m2(j1p,j4m,false,ka, k2,k3,partons,nabove);
const COM Mpmm3 = m3(j1p,j4m,false,ka, k2,k3,partons,nabove);
const COM Mpmp1 = m1(j1p,j4m,true, ka,kb,k2,k3,partons,nabove);
const COM Mpmp2 = m2(j1p,j4m,true, ka, k2,k3,partons,nabove);
const COM Mpmp3 = m3(j1p,j4m,true, ka, k2,k3,partons,nabove);
//Colour factors:
const COM cm1m1=3.;
const COM cm2m2=4./3.;
const COM cm3m3=4./3.;
const COM cm1m2 =3./2.*COM(0.,1.);
const COM cm1m3 = -3./2.*COM(0.,1.);
const COM cm2m3 = -1./6.;
//Square and sum for each helicity config:
const 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)));
const double Mmmp = real(cm1m1*pow(abs(Mmmp1),2)+cm2m2*pow(abs(Mmmp2),2)+
cm3m3*pow(abs(Mmmp3),2)+2.*real(cm1m2*Mmmp1*conj(Mmmp2))+
2.*real(cm1m3*Mmmp1*conj(Mmmp3))+2.*real(cm2m3*Mmmp2*conj(Mmmp3)));
const 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)));
const double Mpmp = real(cm1m1*pow(abs(Mpmp1),2)+cm2m2*pow(abs(Mpmp2),2)+
cm3m3*pow(abs(Mpmp3),2)+2.*real(cm1m2*Mpmp1*conj(Mpmp2))+
2.*real(cm1m3*Mpmp1*conj(Mpmp3))+2.*real(cm2m3*Mpmp2*conj(Mpmp3)));
//Result (averaged, without coupling or t-channel props). Factor of
//2 for the 4 helicity configurations I didn't work out explicitly
HLV prop1 = ka;
for(int i=0; i<=nabove; ++i){
prop1 -= partons[i];
}
const HLV prop2 = prop1 - k2 - k3;
return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/9./4.) /
((ka-partons.front()).m2()*(kb-partons.back()).m2()*prop1.m2()*prop2.m2());
}
} // namespace currents
} // namespace HEJ

File Metadata

Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:32 AM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4856825
Default Alt Text
(230 KB)

Event Timeline