Page MenuHomeHEPForge
Files F9501729

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh
index 76d813d..31070cd 100644
--- a/include/HEJ/Config.hh
+++ b/include/HEJ/Config.hh
@@ -1,232 +1,232 @@
/** \file
* \brief HEJ 2 configuration parameters
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <map>
#include <string>
#include <vector>
#include "fastjet/JetDefinition.hh"
#include "yaml-cpp/yaml.h"
#include "HEJ/Constants.hh"
#include "HEJ/event_types.hh"
#include "HEJ/EWConstants.hh"
#include "HEJ/HiggsCouplingSettings.hh"
#include "HEJ/optional.hh"
#include "HEJ/output_formats.hh"
#include "HEJ/ScaleFunction.hh"
namespace HEJ{
//! Jet parameters
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
double min_pt; /**< Minimum Jet Transverse Momentum */
};
//! Settings for scale variation
struct ScaleConfig{
//! Base scale choices
std::vector<ScaleFunction> base;
//! Factors for multiplicative scale variation
std::vector<double> factors;
//! Maximum ratio between renormalisation and factorisation scale
double max_ratio;
};
//! Settings for random number generator
struct RNGConfig {
//! Random number generator name
std::string name;
//! Optional initial seed
optional<std::string> seed;
};
//! Settings for partial unweighting
struct PartialUnweightConfig {
//! Number of trials for training
- int trials;
+ 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
double max_ext_soft_pt_fraction;
//! The regulator lambda for the subtraction terms
double regulator_lambda = CLAMBDA;
//! Number of resummation configurations to generate per fixed-order event
- int trials;
+ 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;
//! 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
HEJ::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
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{Higgs_coupling},
ew_parameters{ew_parameters},
regulator_lambda{regulator_lambda}
{}
//! Whether to include the logarithmic correction from \f$\alpha_s\f$ running
bool log_correction;
//! Settings for effective Higgs-gluon coupling
HiggsCouplingSettings Higgs_coupling;
//! 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/Event.hh b/include/HEJ/Event.hh
index 927c49c..94363f5 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,355 +1,355 @@
/** \file
* \brief Declares the Event class and helpers
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <memory>
#include <string>
#include <unordered_map>
#include <vector>
#include "boost/iterator/filter_iterator.hpp"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/Particle.hh"
#include "HEJ/RNG.hh"
namespace LHEF {
class HEPEUP;
class HEPRUP;
}
namespace fastjet {
class JetDefinition;
}
namespace HEJ {
struct UnclusteredEvent;
/** @brief An event with clustered jets
*
* This is the main HEJ 2 event class.
* It contains kinematic information including jet clustering,
* parameter (e.g. scale) settings and the event weight.
*/
class Event {
public:
class EventData;
using ConstPartonIterator = boost::filter_iterator<
bool (*)(Particle const &),
std::vector<Particle>::const_iterator
>;
using ConstReversePartonIterator = std::reverse_iterator<
ConstPartonIterator>;
//! No default Constructor
Event() = delete;
//! Event Constructor adding jet clustering to an unclustered event
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
[[deprecated("UnclusteredEvent will be replaced by EventData")]]
Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
);
//! @name Particle Access
//! @{
//! Incoming particles
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Outgoing particles
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
//! Iterator to the first outgoing parton
ConstPartonIterator begin_partons() const;
//! Iterator to the first outgoing parton
ConstPartonIterator cbegin_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator end_partons() const;
//! Iterator to the end of the outgoing partons
ConstPartonIterator cend_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator rbegin_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator crbegin_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator rend_partons() const;
//! Reverse Iterator to the first outgoing parton
ConstReversePartonIterator crend_partons() const;
//! Particle decays
/**
* The key in the returned map corresponds to the index in the
* vector returned by outgoing()
*/
std::unordered_map<size_t, std::vector<Particle>> const & decays() const{
return decays_;
}
//! The jets formed by the outgoing partons, sorted in rapidity
std::vector<fastjet::PseudoJet> const & jets() const{
return jets_;
}
//! @}
//! @name Weight variations
//! @{
//! All chosen parameter, i.e. scale choices (const version)
Parameters<EventParameters> const & parameters() const{
return parameters_;
}
//! All chosen parameter, i.e. scale choices
Parameters<EventParameters> & parameters(){
return parameters_;
}
//! Central parameter choice (const version)
EventParameters const & central() const{
return parameters_.central;
}
//! Central parameter choice
EventParameters & central(){
return parameters_.central;
}
//! Parameter (scale) variations (const version)
std::vector<EventParameters> const & variations() const{
return parameters_.variations;
}
//! Parameter (scale) variations
std::vector<EventParameters> & variations(){
return parameters_.variations;
}
//! Parameter (scale) variation (const version)
/**
* @param i Index of the requested variation
*/
EventParameters const & variations(size_t i) const{
return parameters_.variations.at(i);
}
//! Parameter (scale) variation
/**
* @param i Index of the requested variation
*/
EventParameters & variations(size_t i){
return parameters_.variations.at(i);
}
//! @}
//! Indices of the jets the outgoing partons belong to
/**
* @param jets Jets to be tested
* @returns A vector containing, for each outgoing parton,
* the index in the vector of jets the considered parton
* belongs to. If the parton is not inside any of the
* passed jets, the corresponding index is set to -1.
*/
std::vector<int> particle_jet_indices(
std::vector<fastjet::PseudoJet> const & jets
) const {
return cs_.particle_jet_indices(jets);
}
//! particle_jet_indices() of the Event jets()
std::vector<int> particle_jet_indices() const {
return particle_jet_indices(jets());
}
//! Jet definition used for clustering
fastjet::JetDefinition const & jet_def() const{
return cs_.jet_def();
}
//! Minimum jet transverse momentum
double min_jet_pt() const{
return min_jet_pt_;
}
//! Event type
event_type::EventType type() const{
return type_;
}
//! Give colours to each particle
/**
* @returns true if new colours are generated, i.e. same as is_resummable()
* @details Colour ordering is done according to leading colour in the MRK
* limit, see \cite Andersen:2011zd. This only affects \ref
* is_resummable() "HEJ" configurations, all other \ref event_type
* "EventTypes" will be ignored.
* @note This overwrites all previously set colours.
*/
bool generate_colours(HEJ::RNG &);
//! Check that current colours are leading in the high energy limit
/**
* @details Checks that the colour configuration can be split up in
* multiple, rapidity ordered, non-overlapping ladders. Such
* configurations are leading in the MRK limit, see
* \cite Andersen:2011zd
*
* @note This is _not_ to be confused with \ref is_resummable(), however
* for all resummable states it is possible to create a leading colour
* configuration, see generate_colours()
*/
bool is_leading_colour() const;
/**
* @brief Check if given event could have been produced by HEJ
* @details A HEJ state has to fulfil:
* 1. type() has to be \ref is_resummable() "resummable"
* 2. Soft radiation in the tagging jets contributes at most to
* `max_ext_soft_pt_fraction` of the total jet \f$ p_\perp \f$
*
* @note This is true for any resummed stated produced by the
* EventReweighter or any \ref is_resummable() "resummable" Leading
* Order state.
*
* @param max_ext_soft_pt_fraction Maximum transverse momentum fraction from
* soft radiation in extremal jets
* @param min_extparton_pt Absolute minimal \f$ p_\perp \f$,
* \b deprecated use max_ext_soft_pt_fraction
* instead
* @return True if this state could have been produced by HEJ
*/
bool valid_hej_state(
double max_ext_soft_pt_fraction, double min_extparton_pt = 0.) const;
private:
//! \internal
//! @brief Construct Event explicitly from input.
/** This is only intended to be called from EventData.
*
* \warning The input is taken _as is_, sorting and classification has to be
* done externally, i.e. by EventData
*/
Event(
std::array<Particle, 2> && incoming,
std::vector<Particle> && outgoing,
std::unordered_map<size_t, std::vector<Particle>> && decays,
Parameters<EventParameters> && parameters,
fastjet::JetDefinition const & jet_def,
double const min_jet_pt
);
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
std::unordered_map<size_t, std::vector<Particle>> decays_;
std::vector<fastjet::PseudoJet> jets_;
Parameters<EventParameters> parameters_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
event_type::EventType type_;
}; // end class Event
//! Class to store general Event setup, i.e. Phase space and weights
class Event::EventData {
public:
//! Default Constructor
EventData() = default;
//! Constructor from LesHouches event information
EventData(LHEF::HEPEUP const & hepeup);
//! Constructor with all values given
EventData(
std::array<Particle, 2> const & incoming_,
std::vector<Particle> const & outgoing_,
std::unordered_map<size_t, std::vector<Particle>> const & decays_,
Parameters<EventParameters> const & parameters_
):
incoming(incoming_), outgoing(outgoing_),
decays(decays_), parameters(parameters_)
- {};
+ {}
//! Move Constructor with all values given
EventData(
std::array<Particle, 2> && incoming_,
std::vector<Particle> && outgoing_,
std::unordered_map<size_t, std::vector<Particle>> && decays_,
Parameters<EventParameters> && parameters_
):
incoming(std::move(incoming_)), outgoing(std::move(outgoing_)),
decays(std::move(decays_)), parameters(std::move(parameters_))
{}
//! Generate an Event from the stored EventData.
/**
* @details Do jet clustering and classification.
* Use this to generate an Event.
*
* @note Calling this function destroys EventData
*
* @param jet_def Jet definition
* @param min_jet_pt minimal \f$p_T\f$ for each jet
*
* @returns Full clustered and classified event.
*/
Event cluster(
fastjet::JetDefinition const & jet_def, double const min_jet_pt);
//! Alias for cluster()
Event operator()(
fastjet::JetDefinition const & jet_def, double const min_jet_pt){
return cluster(jet_def, min_jet_pt);
}
//! Sort particles in rapidity
void sort();
//! Reconstruct intermediate particles from final-state leptons
/**
* Final-state leptons are created from virtual photons, W, or Z bosons.
* This function tries to reconstruct such intermediate bosons if they
* are not part of the event record.
*/
void reconstruct_intermediate();
std::array<Particle, 2> incoming;
std::vector<Particle> outgoing;
std::unordered_map<size_t, std::vector<Particle>> decays;
Parameters<EventParameters> parameters;
}; // end class EventData
//! Print Event
std::ostream& operator<<(std::ostream & os, Event const & ev);
//! Square of the partonic centre-of-mass energy \f$\hat{s}\f$
double shat(Event const & ev);
//! Convert an event to a LHEF::HEPEUP
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *);
// put deprecated warning at the end, so don't get the warning inside Event.hh,
// additionally doxygen can not identify [[deprecated]] correctly
struct [[deprecated("UnclusteredEvent will be replaced by EventData")]]
UnclusteredEvent;
//! An event before jet clustering
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.2.0
struct UnclusteredEvent{
//! Default Constructor
UnclusteredEvent() = default;
//! Constructor from LesHouches event information
UnclusteredEvent(LHEF::HEPEUP const & hepeup);
std::array<Particle, 2> incoming; /**< Incoming Particles */
std::vector<Particle> outgoing; /**< Outgoing Particles */
//! Particle decays in the format {outgoing index, decay products}
std::unordered_map<size_t, std::vector<Particle>> decays;
//! Central parameter (e.g. scale) choice
EventParameters central;
std::vector<EventParameters> variations; /**< For parameter variation */
};
} // namespace HEJ
diff --git a/include/HEJ/EventReweighter.hh b/include/HEJ/EventReweighter.hh
index 87a0279..3796a3d 100644
--- a/include/HEJ/EventReweighter.hh
+++ b/include/HEJ/EventReweighter.hh
@@ -1,197 +1,197 @@
/** \file
* \brief Declares the EventReweighter class
*
* EventReweighter is the main class used within HEJ 2. It reweights the
* resummation events.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <memory>
#include <vector>
#include "HEJ/Config.hh"
#include "HEJ/event_types.hh"
#include "HEJ/MatrixElement.hh"
#include "HEJ/Parameters.hh"
#include "HEJ/PDF.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/ScaleFunction.hh"
#include "HEJ/StatusCode.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ {
class Event;
class RNG;
//! Beam parameters
/**
* Currently, only symmetric beams are supported,
* so there is a single beam energy.
*/
struct Beam{
double E; /**< Beam energy */
std::array<ParticleID, 2> type; /**< Beam particles */
};
//! Main class for reweighting events in HEJ.
class EventReweighter{
using EventType = event_type::EventType;
public:
EventReweighter(
Beam beam, /**< Beam Energy */
int pdf_id, /**< PDF ID */
ScaleGenerator scale_gen, /**< Scale settings */
EventReweighterConfig conf, /**< Configuration parameters */
std::shared_ptr<RNG> ran /**< Random number generator */
);
EventReweighter(
LHEF::HEPRUP const & heprup, /**< LHEF event header */
ScaleGenerator scale_gen, /**< Scale settings */
EventReweighterConfig conf, /**< Configuration parameters */
std::shared_ptr<RNG> ran /**< Random number generator */
);
//! Get the used pdf
PDF const & pdf() const;
//! Get event treatment
EventTreatment treatment(EventType type) const;
//! Generate resummation events for a given fixed-order event
/**
* @param ev Fixed-order event corresponding
* to the resummation events
* @param num_events Number of trial resummation configurations.
* @returns A vector of resummation events.
*
* The result vector depends on the type of the input event and the
* treatment of different types as specified in the constructor:
*
* \ref reweight The result vector contains between
* 0 and num_events resummation events.
*
* \ref keep If the input event passes the resummation jet cuts
* the result vector contains one event. Otherwise it is empty.
*
* \ref discard The result vector is empty
*/
std::vector<Event> reweight(
Event const & ev,
- int num_events
+ size_t num_events
);
//! Gives all StatusCodes of the last reweight()
/**
* Each StatusCode corresponds to one tried generation. Only good
* StatusCodes generated an event.
*/
std::vector<StatusCode> const & status() const {
return status_;
}
private:
template<typename... T>
PDF const & pdf(T&& ...);
/** \internal
* \brief main generation/reweighting function:
* generate phase space points and divide out Born factors
*/
std::vector<Event> gen_res_events(
Event const & ev, size_t num_events
);
std::vector<Event> rescale(
Event const & Born_ev, std::vector<Event> events
) const;
/** \internal
* \brief Do the Jets pass the resummation Cuts?
*
* @param ev Event in Question
* @returns 0 or 1 depending on if ev passes Jet Cuts
*/
bool jets_pass_resummation_cuts(Event const & ev) const;
/** \internal
* \brief pdf_factors Function
*
* @param ev Event in Question
* @returns EventFactor due to PDFs
*
* Calculates the Central value and the variation due
* to the PDF choice made.
*/
Weights pdf_factors(Event const & ev) const;
/** \internal
* \brief matrix_elements Function
*
* @param ev Event in question
* @returns EventFactor due to MatrixElements
*
* Calculates the Central value and the variation due
* to the Matrix Element.
*/
Weights matrix_elements(Event const & ev) const;
/** \internal
* \brief Scale-dependent part of fixed-order matrix element
*
* @param ev Event in question
* @returns EventFactor scale variation due to FO-ME.
*
* This is only called to compute the scale variation for events where
* we don't do resummation (e.g. non-FKL).
* Since at tree level the scale dependence is just due to alpha_s,
* it is enough to return the alpha_s(mur) factors in the matrix element.
* The rest drops out in the ratio of (output event ME)/(input event ME),
* so we never have to compute it.
*/
Weights fixed_order_scale_ME(Event const & ev) const;
/** \internal
* \brief Computes the tree level matrix element
*
* @param ev Event in Question
* @returns HEJ approximation to Tree level Matrix Element
*
* This computes the HEJ approximation to the tree level FO
* Matrix element which is used within the LO weighting process.
*/
double tree_matrix_element(Event const & ev) const;
//! \internal General parameters
EventReweighterConfig param_;
//! \internal Beam energy
double E_beam_;
//! \internal PDF
PDF pdf_;
//! \internal Object to calculate the square of the matrix element
MatrixElement MEt2_;
//! \internal Object to calculate event renormalisation and factorisation scales
ScaleGenerator scale_gen_;
//! \internal random number generator
std::shared_ptr<RNG> ran_;
//! \internal StatusCode of each attempt
std::vector<StatusCode> status_;
};
template<typename... T>
PDF const & EventReweighter::pdf(T&&... t){
return pdf_ = PDF{std::forward<T>(t)...};
}
} // namespace HEJ
diff --git a/include/HEJ/HDF5Reader.hh b/include/HEJ/HDF5Reader.hh
index 4497d3f..66efc75 100644
--- a/include/HEJ/HDF5Reader.hh
+++ b/include/HEJ/HDF5Reader.hh
@@ -1,50 +1,50 @@
/** \file
* \brief Header file for reading events in the HDF5 event format.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <string>
#include "HEJ/EventReader.hh"
namespace HEJ{
//! Class for reading events from a file in the HDF5 file format
/**
* @details This format is specified in \cite Hoeche:2019rti.
*/
class HDF5Reader : public EventReader{
public:
HDF5Reader() = delete;
//! Contruct object reading from the given file
explicit HDF5Reader(std::string const & filename);
//! Read an event
bool read_event() override;
//! Access header text
std::string const & header() const override;
//! Access run information
LHEF::HEPRUP const & heprup() const override;
//! Access last read event
LHEF::HEPEUP const & hepeup() const override;
//! Get number of events
HEJ::optional<size_t> number_events() const override;
- ~HDF5Reader();
+ ~HDF5Reader() override;
private:
struct HDF5ReaderImpl;
std::unique_ptr<HDF5ReaderImpl> impl_;
};
} // namespace HEJ
diff --git a/include/HEJ/HDF5Writer.hh b/include/HEJ/HDF5Writer.hh
index d2542af..f74464c 100644
--- a/include/HEJ/HDF5Writer.hh
+++ b/include/HEJ/HDF5Writer.hh
@@ -1,54 +1,54 @@
/** \file
* \brief Contains the EventWriter for HDF5 Output.
*
* The output format is specified in arXiv:1905.05120.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include "HEJ/EventWriter.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ {
class Event;
//! This is an event writer specifically for HDF5 output.
/**
* \internal Implementation note: This uses the pimpl ("pointer to
* implementation") idiom. HDF5 support is optional. Without pimpl,
* we would have to specify whether HDF5 is available via the
* preprocessor whenever this header is included. We don't want to
* burden users of the HEJ library (for example the HEJ fixed-order
* generator) with those details
*/
class HDF5Writer: public EventWriter{
public:
//! Constructor
/**
* @param file name of the output file
* @param heprup general process information
*/
HDF5Writer(std::string const & file, LHEF::HEPRUP heprup);
HDF5Writer() = delete;
//! Write an event to the output file
void write(Event const & ev) override;
- ~HDF5Writer();
+ ~HDF5Writer() override;
private:
struct HDF5WriterImpl;
std::unique_ptr<HDF5WriterImpl> impl_;
};
} // namespace HEJ
diff --git a/include/HEJ/HepMC2Writer.hh b/include/HEJ/HepMC2Writer.hh
index 2afe462..7f2deae 100644
--- a/include/HEJ/HepMC2Writer.hh
+++ b/include/HEJ/HepMC2Writer.hh
@@ -1,53 +1,53 @@
/** \file
* \brief Contains the EventWriter for HepMC Output.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include "HEJ/EventWriter.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ {
class Event;
//! This is an event writer specifically for HepMC output.
/**
* \internal Implementation note:
* This uses the pimpl ("pointer to implementation") idiom.
* HepMC support is optional and the implementation depends on the
* HepMC version. Without pimpl, we would have to specify the HepMC version
* via the preprocessor whenever this header is included. We don't want to
* burden users of the HEJ library (for example the HEJ fixed-order generator)
* with those details
*/
class HepMC2Writer: public EventWriter{
public:
//! Constructor
/**
* @param file name of the output file
* @param heprup general process information
*/
HepMC2Writer(std::string const & file, LHEF::HEPRUP heprup);
HepMC2Writer() = delete;
//! Write an event to the output file
void write(Event const & ev) override;
- ~HepMC2Writer();
+ ~HepMC2Writer() override;
private:
struct HepMC2WriterImpl;
std::unique_ptr<HepMC2WriterImpl> impl_;
};
} // namespace HEJ
diff --git a/include/HEJ/HepMC3Writer.hh b/include/HEJ/HepMC3Writer.hh
index c02c8b4..e1d15bc 100644
--- a/include/HEJ/HepMC3Writer.hh
+++ b/include/HEJ/HepMC3Writer.hh
@@ -1,52 +1,52 @@
/** \file
* \brief Contains the EventWriter for HepMC3 Output.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include "HEJ/EventWriter.hh"
namespace LHEF {
class HEPRUP;
}
namespace HEJ {
class Event;
//! This is an event writer specifically for HepMC3 output.
/**
* \internal Implementation note:
* This uses the pimpl ("pointer to implementation") idiom.
* HepMC3 support is optional and the implementation depends on the
* HepMC3 version. Without pimpl, we would have to specify the HepMC3 version
* via the preprocessor whenever this header is included. We don't want to
* burden users of the HEJ library (for example the HEJ fixed-order generator)
* with those details
*/
class HepMC3Writer: public EventWriter{
public:
//! Constructor
/**
* @param file name of the output file
* @param heprup general process information
*/
HepMC3Writer(std::string const & file, LHEF::HEPRUP heprup);
HepMC3Writer() = delete;
//! Write an event to the output file
void write(Event const & ev) override;
- ~HepMC3Writer();
+ ~HepMC3Writer() override;
private:
struct HepMC3WriterImpl;
std::unique_ptr<HepMC3WriterImpl> impl_;
};
} // namespace HEJ
diff --git a/include/HEJ/Mixmax.hh b/include/HEJ/Mixmax.hh
index ff3fe2f..61f2c52 100644
--- a/include/HEJ/Mixmax.hh
+++ b/include/HEJ/Mixmax.hh
@@ -1,35 +1,35 @@
/** \file
* \brief The Mixmax random number generator
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include "CLHEP/Random/MixMaxRng.h"
#include "CLHEP/Random/Randomize.h"
#include "HEJ/RNG.hh"
namespace HEJ {
//! MIXMAX random number generator
/**
* For details on MIXMAX, see \cite Savvidy:2014ana
*/
class Mixmax : public DefaultRNG {
public:
Mixmax() = default;
- Mixmax(long seed): ran_{seed} {};
+ Mixmax(long seed): ran_{seed} {}
//! Generate pseudorandom number between 0 and 1
double flat() override {
return ran_.flat();
- };
+ }
private:
CLHEP::MixMaxRng ran_;
};
}
diff --git a/include/HEJ/Parameters.hh b/include/HEJ/Parameters.hh
index 0973698..84a640a 100644
--- a/include/HEJ/Parameters.hh
+++ b/include/HEJ/Parameters.hh
@@ -1,163 +1,164 @@
/** \file
* \brief Containers for Parameter variations, e.g. different Weights
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#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;
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 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;
}
//! 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;
//! Actual factorisation scale divided by central scale
double muf_factor;
ParameterDescription() = default;
ParameterDescription(
std::string scale_name, double mur_factor, double muf_factor
):
- scale_name{scale_name}, mur_factor{mur_factor}, muf_factor{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 */
//! 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;
}
};
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;
}
//! @}
} // namespace HEJ
diff --git a/include/HEJ/RivetAnalysis.hh b/include/HEJ/RivetAnalysis.hh
index 7c4d1bd..79a7598 100644
--- a/include/HEJ/RivetAnalysis.hh
+++ b/include/HEJ/RivetAnalysis.hh
@@ -1,72 +1,72 @@
/** \file
* \brief HEJ 2 interface to rivet analyses
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <memory>
#include <string>
#include <vector>
#include "LHEF/LHEF.h"
#include "HEJ/Analysis.hh"
#include "HEJ/optional.hh"
namespace Rivet {
class AnalysisHandler;
}
namespace YAML {
class Node;
}
namespace HEJ {
/**
* @brief Class representing a Rivet analysis
*
* This class inherits from Analysis and can therefore be used
* like any other HEJ 2 analysis.
*/
class RivetAnalysis: public HEJ::Analysis {
public:
static std::unique_ptr<Analysis> create(
YAML::Node const & config, LHEF::HEPRUP const & heprup);
//! Constructor
/**
* @param config Configuration parameters
* @param heprup General run informations
*
* config["rivet"] should be the name of a single Rivet analysis or
* a list of Rivet analyses. config["output"] is the prefix for
* the .yoda output files.
*/
RivetAnalysis(YAML::Node const & config, LHEF::HEPRUP const & heprup);
- ~RivetAnalysis();
+ ~RivetAnalysis() override;
//! Pass an event to the underlying Rivet analysis
void fill(HEJ::Event const & event, HEJ::Event const &) override;
bool pass_cuts(HEJ::Event const &, HEJ::Event const &) override
{return true;} //!< no additional cuts are applied
void finalise() override;
private:
std::vector<std::string> analyses_names_;
std::string output_name_;
LHEF::HEPRUP heprup_;
/// struct to organise the infos per rivet run/scale setting
struct rivet_info;
std::vector<rivet_info> rivet_runs_;
/**
* \internal
* @brief Calculates the scale variation from the first event for the output
* file
*/
void init(HEJ::Event const & event);
bool first_event_;
};
} // namespace HEJ
diff --git a/include/HEJ/detail/HepMCInterface_common.hh b/include/HEJ/detail/HepMCInterface_common.hh
index 44a1a16..3de7956 100644
--- a/include/HEJ/detail/HepMCInterface_common.hh
+++ b/include/HEJ/detail/HepMCInterface_common.hh
@@ -1,86 +1,86 @@
/** \file
* \brief Template functions shared between HepMC2 and HepMC3
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include "HEJ/Event.hh"
#include "HEJ/Particle.hh"
namespace HEJ{
namespace detail_HepMC{
template<class FourVector>
FourVector to_FourVector(Particle const & sp){
return {sp.px(), sp.py(), sp.pz(), sp.E()};
}
constexpr int status_beam = 4;
constexpr int status_in = 11;
constexpr int status_decayed = 2;
constexpr int status_out = 1;
// Helper classes depending on HepMC version V
template<int V> struct HepMCVersion;
template<int V>
using GenEvent = typename HepMCVersion<V>::GenEvent;
template<int V>
using Beam = typename HepMCVersion<V>::Beam;
template<int V>
auto make_particle_ptr(
Particle const & sp, int status
);
template<int V>
auto make_vx_ptr();
/**
* @brief Template class to initialise the kinematics for HepMC
*
* @tparam V HepMC major version
* @param event HEJ event
* @param beam beam particles
* @param out_ev new HepMC event
*/
template<int V>
GenEvent<V> HepMC_init_kinematics(
Event const & event, Beam<V> const & beam, GenEvent<V> && out_ev
){
out_ev.set_beam_particles(beam[0], beam[1]);
auto vx = make_vx_ptr<V>();
for(size_t i=0; i<event.incoming().size(); ++i){
auto particle = make_particle_ptr<V>(event.incoming()[i], status_in);
auto vx_beam = make_vx_ptr<V>();
vx_beam->add_particle_in(beam[i]);
vx_beam->add_particle_out(particle);
out_ev.add_vertex(vx_beam);
vx->add_particle_in(particle);
}
for(size_t i=0; i < event.outgoing().size(); ++i){
auto const & out = event.outgoing()[i];
auto particle = make_particle_ptr<V>(out, status_out);
const int status = event.decays().count(i)?status_decayed:status_out;
particle->set_status(status);
if( status == status_decayed ){
auto vx_decay = make_vx_ptr<V>();
vx_decay->add_particle_in(particle);
- for( auto const & out: event.decays().at(i) ){
+ for( auto const & decay: event.decays().at(i) ){
vx_decay->add_particle_out(
- make_particle_ptr<V>(out, status_out)
+ make_particle_ptr<V>(decay, status_out)
);
}
out_ev.add_vertex(vx_decay);
}
vx->add_particle_out(particle);
}
out_ev.add_vertex(vx);
return out_ev;
}
} // namespace detail_HepMC
} // namespace HEJ
diff --git a/include/HEJ/exceptions.hh b/include/HEJ/exceptions.hh
index 765d5a3..29dab25 100644
--- a/include/HEJ/exceptions.hh
+++ b/include/HEJ/exceptions.hh
@@ -1,58 +1,58 @@
/** \file
* \brief Custom exception classes
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <stdexcept>
#include <string>
namespace HEJ {
//! Exception indicating wrong option type
/**
* This exception is thrown if a configuration option has
* the wrong type (e.g. 'trials' is not set to a number)
*/
struct invalid_type: std::invalid_argument {
explicit invalid_type(std::string const & what):
- std::invalid_argument{what} {};
+ std::invalid_argument{what} {}
explicit invalid_type(char const * what):
- std::invalid_argument{what} {};
+ std::invalid_argument{what} {}
};
//! Exception indicating unknown option
/**
* This exception is thrown if an unknown configuration option
* is set (e.g. the 'trials' setting is misspelt as 'trails')
*/
struct unknown_option: std::invalid_argument {
explicit unknown_option(std::string const & what):
- std::invalid_argument{what} {};
+ std::invalid_argument{what} {}
explicit unknown_option(char const * what):
- std::invalid_argument{what} {};
+ std::invalid_argument{what} {}
};
//! Exception indicating missing option setting
/**
* This exception is thrown if a mandatory configuration option
* (e.g. 'trials') is not set.
*/
struct missing_option: std::logic_error {
explicit missing_option(std::string const & what):
- std::logic_error{what} {};
+ std::logic_error{what} {}
explicit missing_option(char const * what):
- std::logic_error{what} {};
+ std::logic_error{what} {}
};
//! Exception indicating functionality that has not been implemented yet
struct not_implemented: std::logic_error {
explicit not_implemented(std::string const & what):
- std::logic_error{what} {};
+ std::logic_error{what} {}
explicit not_implemented(char const * what):
- std::logic_error{what} {};
+ std::logic_error{what} {}
};
} // namespace HEJ
diff --git a/include/HEJ/jets.hh b/include/HEJ/jets.hh
index 43a81ef..ab361c0 100644
--- a/include/HEJ/jets.hh
+++ b/include/HEJ/jets.hh
@@ -1,403 +1,403 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
/** \file
* \brief Functions computing the square of current contractions in pure jets.
*
* This file contains all the necessary functions to compute the
* current contractions for all valid pure jet HEJ processes, which
* so far is FKL and unordered processes. It will also contain some
* pure jet ME components used in other process ME calculations
*
* @TODO add a namespace
*/
#pragma once
#include <complex>
#include <ostream>
#include <vector>
#include "CLHEP/Vector/LorentzVector.h"
typedef std::complex<double> COM;
typedef COM current[4];
typedef CLHEP::HepLorentzVector HLV;
//! Square of qQ->qQ Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qQ Scattering
*/
double ME_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qQbar Scattering
*
* @note this can be used for qbarQ->qbarQ Scattering by inputting arguments
* appropriately.
*/
double ME_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarQbar Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*/
double ME_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qg->qg Scattering
*
* @note this can be used for gq->gq Scattering by inputting arguments
* appropriately.
*/
double ME_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* @note this can be used for gqbar->gqbar Scattering by inputting arguments
* appropriately.
*/
double ME_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of gg->gg Pure Jets Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for gg->gg Scattering
*/
double ME_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in);
// Unordered Backwards contributions:
//! Square of qQ->qQ Pure Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qQ Scattering
*/
double ME_unob_qQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQ->qbarQ Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qbarQ->qbarQ Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qbarQ(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qQbar->qQbar Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qQbar Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarQbar->qbarQbar Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qbarQbar(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qg->qg Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qg->qg Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qg(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of qbarg->qbarg Pure Jets Unordered backwards Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered gluon
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* @note this can be used for unof contributions by inputting
* arguments appropriately.
*/
double ME_unob_qbarg(HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in);
//! Square of gQ->qbarqQ Pure Jets Extremal qqx backwards Scattering Current
/**
* @param pgin Momentum of incoming gluon
* @param pqout Momentum of Quark from split
* @param pqbarout Momentum of Anti-quark from split
* @param p2out Momentum of Outgoing forwards leg
* @param p2in Momentum of Incoming forwards leg
* @returns Square of the current contractions for gQ->qbarqg Scattering
*
* @note this can be used for Exqqxf contributions by inputting
* arguments appropriately.
*/
double ME_Exqqx_qbarqQ(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in);
//! Square of gQ->qqbarQ Pure Jets Extremal qqx backwards Scattering Current
/**
* @param pgin Momentum of incoming gluon
* @param pqout Momentum of Quark from split
* @param pqbarout Momentum of Anti-quark from split
* @param p2out Momentum of Outgoing forwards leg
* @param p2in Momentum of Incoming forwards leg
* @returns Square of the current contractions for gQ->qqbarg Scattering
*
* @note this can be used for Exqqxf contributions by inputting
* arguments appropriately.
*/
double ME_Exqqx_qqbarQ(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in);
//! Square of gg->qbarqg Pure Jets Extremal qqx backwards Scattering Current
/**
* @param pgin Momentum of incoming gluon
* @param pqout Momentum of Quark from split
* @param pqbarout Momentum of Anti-quark from split
* @param p2out Momentum of Outgoing forwards leg
* @param p2in Momentum of Incoming forwards leg
* @returns Square of the current contractions for gg->qbarqg Scattering
*
* @note this can be used for Exqqxf contributions by inputting
* arguments appropriately.
*/
double ME_Exqqx_qbarqg(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in);
//! Square of gg->qqbarg Pure Jets Extremal qqx backwards Scattering Current
/**
* @param pgin Momentum of incoming gluon
* @param pqout Momentum of Quark from split
* @param pqbarout Momentum of Anti-quark from split
* @param p2out Momentum of Outgoing forwards leg
* @param p2in Momentum of Incoming forwards leg
* @returns Square of the current contractions for gg->qqbarg Scattering
*
* @note this can be used for Exqqxf contributions by inputting
* arguments appropriately.
*/
double ME_Exqqx_qqbarg(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in);
//! Square of qq->qQQbarq Pure Jets Central qqx Scattering Current
/**
* @param ka Momentum of incoming leg a
* @param kb Momentum of incoming leg b
* @param partons std::vector<HLV> outgoing partons
* @param aqlinepa Is leg a an anti-quark?
* @param aqlinepb Is leg b an anti-quark?
* @param qqxmarker Is anti-quark further back in rapidity than quark (qqx pair)
* @param nabove Number of gluons emitted above qqx pair (back in rap)
* @returns Square of the current contractions for qq->qQQxq Scattering
*/
double ME_Cenqqx_qq(HLV ka, HLV kb, std::vector<HLV> partons, bool aqlinepa,
bool aqlinepb, bool qqxmarker, int nabove);
/** \class CCurrent jets.hh "include/HEJ/jets.hh"
* \brief This is the a new class structure for currents.
*/
class CCurrent {
public:
CCurrent(COM sc0, COM sc1, COM sc2, COM sc3)
:c0(sc0), c1(sc1), c2(sc2), c3(sc3)
- {};
+ {}
CCurrent(const HLV p) {
c0=p.e();
c1=p.px();
c2=p.py();
c3=p.pz();
}
CCurrent() {}
CCurrent operator+(const CCurrent& other);
CCurrent operator-(const CCurrent& other);
CCurrent operator*(const double x);
CCurrent operator*(const COM x);
CCurrent operator/(const double x);
CCurrent operator/(const COM x);
friend std::ostream& operator<<(std::ostream& os, const CCurrent& cur);
COM dot(HLV p1);
COM dot(CCurrent p1);
COM c0, c1, c2, c3;
};
/* std::ostream& operator <<(std::ostream& os, const CCurrent& cur); */
CCurrent operator*(double x, CCurrent& m);
CCurrent operator*(COM x, CCurrent& m);
CCurrent operator/(double x, CCurrent& m);
CCurrent operator/(COM x, CCurrent& m);
//! Current <incoming state | mu | outgoing state>
/**
* This is a wrapper function around \see joi() note helicity flip to
* give same answer.
*/
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur);
//! Current <outgoing state | mu | outgoing state>
/**
* @param pi bra state momentum
* @param heli helicity of pi
* @param pj ket state momentum
* @param helj helicity of pj. (must be same as heli)
* @param cur reference to current which is saved.
*
* This function is for building <i (out)| mu |j (out)> currents. It
* must be called with pi as the bra, and pj as the ket.
*
* @TODO Remove heli/helj and just have helicity of current as argument.
*/
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* @param pout bra state momentum
* @param helout helicity of pout
* @param pin ket state momentum
* @param helin helicity of pin. (must be same as helout)
* @param cur reference to current which is saved.
*
* This function is for building <out| mu |in> currents. It must be
* called with pout as the bra, and pin as the ket. jio calls this
* with flipped helicity
*
* @TODO Remove helout/helin and just have helicity of current as argument.
*/
void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur);
//! Current <outgoing state | mu | incoming state>
/**
* This is a wrapper function around the void function of the same name. \see joi
*
* @TODO This is never used
*/
CCurrent joi(HLV pout, bool helout, HLV pin, bool helin);
//! Current <incoming state | mu | outgoing state>
/**
* This is a wrapper function around the void function of the same name. \see jio
*/
CCurrent jio(HLV pout, bool helout, HLV pin, bool helin);
//! Current <outgoing state | mu | outgoing state>
/**
* This is a wrapper function around the void function of the same name. \see joo
*/
CCurrent joo(HLV pout, bool helout, HLV pin, bool helin);
inline COM cdot(const current & j1, const current & j2) {
return j1[0]*j2[0]-j1[1]*j2[1]-j1[2]*j2[2]-j1[3]*j2[3];
}
inline COM cdot(const HLV & p, const current & j1) {
return j1[0]*p.e()-j1[1]*p.x()-j1[2]*p.y()-j1[3]*p.z();
}
inline void cmult(const COM & factor, const current & j1, current &cur) {
cur[0]=factor*j1[0];
cur[1]=factor*j1[1];
cur[2]=factor*j1[2];
cur[3]=factor*j1[3];
}
// WHY!?!
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, const current & j5, current &sum
) {
sum[0]=j1[0]+j2[0]+j3[0]+j4[0]+j5[0];
sum[1]=j1[1]+j2[1]+j3[1]+j4[1]+j5[1];
sum[2]=j1[2]+j2[2]+j3[2]+j4[2]+j5[2];
sum[3]=j1[3]+j2[3]+j3[3]+j4[3]+j5[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, current &sum
) {
sum[0] = j1[0] + j2[0] + j3[0] + j4[0];
sum[1] = j1[1] + j2[1] + j3[1] + j4[1];
sum[2] = j1[2] + j2[2] + j3[2] + j4[2];
sum[3] = j1[3] + j2[3] + j3[3] + j4[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
current &sum
) {
sum[0]=j1[0]+j2[0]+j3[0];
sum[1]=j1[1]+j2[1]+j3[1];
sum[2]=j1[2]+j2[2]+j3[2];
sum[3]=j1[3]+j2[3]+j3[3];
}
inline void cadd(const current & j1, const current & j2, current &sum) {
sum[0]=j1[0]+j2[0];
sum[1]=j1[1]+j2[1];
sum[2]=j1[2]+j2[2];
sum[3]=j1[3]+j2[3];
}
inline double abs2(const COM & a) {
return (a*conj(a)).real();
}
inline double vabs2(const CCurrent & cur) {
return abs2(cur.c0)-abs2(cur.c1)-abs2(cur.c2)-abs2(cur.c3);
}
inline double vre(const CCurrent & a, const CCurrent & b) {
return real(a.c0*conj(b.c0)-a.c1*conj(b.c1)-a.c2*conj(b.c2)-a.c3*conj(b.c3));
}
//! @TODO These are not currents and should be moved elsewhere.
double K_g(double p1minus, double paminus);
double K_g(HLV const & pout, HLV const & pin);
diff --git a/src/Event.cc b/src/Event.cc
index b33c3a1..e54bd8a 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,987 +1,987 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Event.hh"
#include <algorithm>
#include <assert.h>
#include <iterator>
#include <numeric>
#include <unordered_set>
#include <utility>
#include "LHEF/LHEF.h"
#include "fastjet/JetDefinition.hh"
#include "HEJ/Constants.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/PDG_codes.hh"
namespace HEJ{
namespace {
constexpr int status_in = -1;
constexpr int status_decayed = 2;
constexpr int status_out = 1;
//! true if leptonic W decay
bool valid_W_decay( int const w_type, // sign of W
std::vector<Particle> const & decays
){
if(decays.size() != 2) // no 1->2 decay
return false;
const int pidsum = decays[0].type + decays[1].type;
if( std::abs(pidsum) != 1 || pidsum != w_type ) // correct charge
return false;
// leptonic decay (only check first, second follows from pidsum)
if( w_type == 1 ) // W+
return is_antilepton(decays[0]) || is_neutrino(decays[0]);
// W-
return is_lepton(decays[0]) || is_antineutrino(decays[0]);
}
/// @name helper functions to determine event type
//@{
/**
* \brief check if final state valid for HEJ
*
* check if there is at most one photon, W, H, Z in the final state
* and all the rest are quarks or gluons
*/
bool final_state_ok(Event const & ev){
std::vector<Particle> const & outgoing = ev.outgoing();
if(ev.decays().size() > 1) // at most one decay
return false;
bool has_AWZH_boson = false;
for( size_t i=0; i<outgoing.size(); ++i ){
auto const & out{ outgoing[i] };
if(is_AWZH_boson(out.type)){
// at most one boson
if(has_AWZH_boson) return false;
has_AWZH_boson = true;
// valid decay for W
if(std::abs(out.type) == ParticleID::Wp){
// exactly 1 decay of W
if( ev.decays().size() != 1 || ev.decays().cbegin()->first != i )
return false;
if( !valid_W_decay(out.type>0?+1:-1, ev.decays().cbegin()->second) )
return false;
}
}
else if(! is_parton(out.type)) return false;
}
return true;
}
/**
* returns all EventTypes implemented in HEJ
*/
size_t implemented_types(std::vector<Particle> const & bosons){
using namespace event_type;
if(bosons.empty()) return FKL | unob | unof | qqxexb | qqxexf | qqxmid;
if(bosons.size()>1) return non_resummable; // multi boson
switch (bosons[0].type) {
case ParticleID::Wp:
case ParticleID::Wm:
return FKL | unob | unof | qqxexb | qqxexf | qqxmid;
case ParticleID::h:
return FKL | unob | unof;
default:
return non_resummable;
}
}
/**
* \brief function which determines if type change is consistent with Wp emission.
* @param in incoming Particle id
* @param out outgoing Particle id
* @param qqx Current both incoming/both outgoing?
*
* \see is_Wm_Change
*/
bool is_Wp_Change(ParticleID in, ParticleID out, bool qqx){
if(!qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out== (in-1);
if( qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out==-(in+1);
return false;
}
/**
* \brief function which determines if type change is consistent with Wm emission.
* @param in incoming Particle id
* @param out outgoing Particle id
* @param qqx Current both incoming/both outgoing?
*
* Ensures that change type of quark line is possible by a flavour changing
* Wm emission. Allows checking of qqx currents also.
*/
bool is_Wm_Change(ParticleID in, ParticleID out, bool qqx){
if(!qqx && (in== 1 || in==-2 || in== 3 || in==-4)) return out== (in+1);
if( qqx && (in==-1 || in== 2 || in==-3 || in== 4)) return out==-(in-1);
return false;
}
/**
* \brief checks if particle type remains same from incoming to outgoing
* @param in incoming Particle
* @param out outgoing Particle
* @param qqx Current both incoming/outgoing?
*/
bool no_flavour_change(ParticleID in, ParticleID out, bool qqx){
const int qqxCurrent = qqx?-1:1;
if(abs(in)<=6 || in==pid::gluon) return (in==out*qqxCurrent);
else return false;
}
bool has_2_jets(Event const & event){
return event.jets().size() >= 2;
}
/**
* \brief check if we have a valid Impact factor
* @param in incoming Particle
* @param out outgoing Particle
* @param qqx Current both incoming/outgoing?
- * @param qqx returns +1 if Wp, -1 if Wm, else 0
+ * @param W_change returns +1 if Wp, -1 if Wm, else 0
*/
bool is_valid_impact_factor(
ParticleID in, ParticleID out, bool qqx, int & W_change
){
if( no_flavour_change(in, out, qqx) ){
return true;
}
if( is_Wp_Change(in, out, qqx) ) {
W_change+=1;
return true;
}
if( is_Wm_Change(in, out, qqx) ) {
W_change-=1;
return true;
}
return false;
}
//! Returns all possible classifications from the impact factors
// the beginning points are changed s.t. after the the classification they
// point to the beginning of the (potential) FKL chain
// sets W_change: + if Wp change
// 0 if no change
// - if Wm change
// This function can be used with forward & backwards iterators
template<class OutIterator>
size_t possible_impact_factors(
ParticleID incoming_id, // incoming
OutIterator & begin_out, OutIterator const & end_out, // outgoing
int & W_change, std::vector<Particle> const & boson,
bool const backward // backward?
){
using namespace event_type;
assert(boson.size() < 2);
// keep track of all states that we don't test
size_t not_tested = qqxmid;
if(backward)
not_tested |= unof | qqxexf;
else
not_tested |= unob | qqxexb;
// Is this LL current?
if( is_valid_impact_factor(incoming_id, begin_out->type, false, W_change) ){
++begin_out;
return not_tested | FKL;
}
// or NLL current?
// -> needs two partons in two different jets
if( std::distance(begin_out, end_out)>=2
){
// Is this unordered emisson?
if( incoming_id!=pid::gluon && begin_out->type==pid::gluon ){
if( is_valid_impact_factor(
incoming_id, (begin_out+1)->type, false, W_change )
){
// veto Higgs inside uno
assert((begin_out+1)<end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < (begin_out+1)->rapidity())
||(!backward && boson.front().rapidity() > (begin_out+1)->rapidity()))
return non_resummable;
}
begin_out+=2;
return not_tested | (backward?unob:unof);
}
}
// Is this QQbar?
else if( incoming_id==pid::gluon ){
if( is_valid_impact_factor(
begin_out->type, (begin_out+1)->type, true, W_change )
){
// veto Higgs inside qqx
assert((begin_out+1)<end_out);
if( !boson.empty() && boson.front().type == ParticleID::h
){
if( (backward && boson.front().rapidity() < (begin_out+1)->rapidity())
||(!backward && boson.front().rapidity() > (begin_out+1)->rapidity()))
return non_resummable;
}
begin_out+=2;
return not_tested | (backward?qqxexb:qqxexf);
}
}
}
return non_resummable;
}
//! Returns all possible classifications from central emissions
// the beginning points are changed s.t. after the the classification they
// point to the end of the emission chain
// sets W_change: + if Wp change
// 0 if no change
// - if Wm change
template<class OutIterator>
size_t possible_central(
OutIterator & begin_out, OutIterator const & end_out,
int & W_change, std::vector<Particle> const & boson
){
using namespace event_type;
assert(boson.size() < 2);
// if we already passed the central chain,
// then it is not a valid all-order state
if(std::distance(begin_out, end_out) < 0) return non_resummable;
// keep track of all states that we don't test
size_t possible = unob | unof
| qqxexb | qqxexf;
// Find the first non-gluon/non-FKL
while( (begin_out->type==pid::gluon) && (begin_out<end_out) ){
++begin_out;
}
// end of chain -> FKL
if( begin_out==end_out ){
return possible | FKL;
}
// is this a qqbar-pair?
// needs two partons in two separate jets
if( is_valid_impact_factor(
begin_out->type, (begin_out+1)->type, true, W_change )
){
// veto Higgs inside qqx
if( !boson.empty() && boson.front().type == ParticleID::h
&& boson.front().rapidity() > begin_out->rapidity()
&& boson.front().rapidity() < (begin_out+1)->rapidity()
){
return non_resummable;
}
begin_out+=2;
// remaining chain should be pure gluon/FKL
for(; begin_out<end_out; ++begin_out){
if(begin_out->type != pid::gluon) return non_resummable;
}
return possible | qqxmid;
}
return non_resummable;
}
/**
* \brief Checks for all event types
* @param ev Event
* @returns Event Type
*
*/
event_type::EventType classify(Event const & ev){
using namespace event_type;
if(! has_2_jets(ev))
return no_2_jets;
// currently we can't handle multiple boson states in the ME. So they are
// considered "bad_final_state" even though the "classify" could work with
// them.
if(! final_state_ok(ev))
return bad_final_state;
// initialise variables
auto const & in = ev.incoming();
auto const & out = filter_partons(ev.outgoing());
assert(std::distance(begin(in), end(in)) == 2);
assert(out.size() >= 2);
assert(std::distance(begin(out), end(out)) >= 2);
assert(std::is_sorted(begin(out), end(out), rapidity_less{}));
auto const boson{ filter_AWZH_bosons(ev.outgoing()) };
// we only allow one boson through final_state_ok
assert(boson.size()<=1);
// keep track of potential W couplings, at the end the sum should be 0
int remaining_Wp = 0;
int remaining_Wm = 0;
if(!boson.empty() && abs(boson.front().type) == ParticleID::Wp ){
if(boson.front().type>0) ++remaining_Wp;
else ++remaining_Wm;
}
int W_change = 0;
// range for current checks
auto begin_out{out.cbegin()};
auto end_out{out.crbegin()};
size_t final_type = ~(no_2_jets | bad_final_state);
// check forward impact factor
final_type &= possible_impact_factors(
in.front().type,
begin_out, end_out.base(),
W_change, boson, true );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
W_change = 0;
// check backward impact factor
final_type &= possible_impact_factors(
in.back().type,
end_out, std::make_reverse_iterator(begin_out),
W_change, boson, false );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
W_change = 0;
// check central emissions
final_type &= possible_central(
begin_out, end_out.base(), W_change, boson );
if( final_type == non_resummable )
return non_resummable;
if(W_change>0) remaining_Wp-=W_change;
else if(W_change<0) remaining_Wm+=W_change;
// Check whether the right number of Ws are present
if( remaining_Wp != 0 || remaining_Wm != 0 ) return non_resummable;
// result has to be unique
if( (final_type & (final_type-1)) != 0) return non_resummable;
// check that each sub processes is implemented
// (has to be done at the end)
if( (final_type & ~implemented_types(boson)) != 0 )
return non_resummable;
return static_cast<EventType>(final_type);
}
//@}
- Particle extract_particle(LHEF::HEPEUP const & hepeup, int i){
+ Particle extract_particle(LHEF::HEPEUP const & hepeup, size_t i){
const ParticleID id = static_cast<ParticleID>(hepeup.IDUP[i]);
const fastjet::PseudoJet momentum{
hepeup.PUP[i][0], hepeup.PUP[i][1],
hepeup.PUP[i][2], hepeup.PUP[i][3]
};
if(is_parton(id))
return Particle{ id, std::move(momentum), hepeup.ICOLUP[i] };
return Particle{ id, std::move(momentum), {} };
}
bool is_decay_product(std::pair<int, int> const & mothers){
if(mothers.first == 0) return false;
return mothers.second == 0 || mothers.first == mothers.second;
}
} // namespace anonymous
Event::EventData::EventData(LHEF::HEPEUP const & hepeup){
parameters.central = EventParameters{
hepeup.scales.mur, hepeup.scales.muf, hepeup.XWGTUP
};
size_t in_idx = 0;
for (int i = 0; i < hepeup.NUP; ++i) {
// skip decay products
// we will add them later on, but we have to ensure that
// the decayed particle is added before
if(is_decay_product(hepeup.MOTHUP[i])) continue;
auto particle = extract_particle(hepeup, i);
// needed to identify mother particles for decay products
particle.p.set_user_index(i+1);
if(hepeup.ISTUP[i] == status_in){
if(in_idx > incoming.size()) {
throw std::invalid_argument{
"Event has too many incoming particles"
};
}
incoming[in_idx++] = std::move(particle);
}
else outgoing.emplace_back(std::move(particle));
}
// add decay products
for (int i = 0; i < hepeup.NUP; ++i) {
if(!is_decay_product(hepeup.MOTHUP[i])) continue;
const int mother_id = hepeup.MOTHUP[i].first;
const auto mother = std::find_if(
begin(outgoing), end(outgoing),
[mother_id](Particle const & particle){
return particle.p.user_index() == mother_id;
}
);
if(mother == end(outgoing)){
throw std::invalid_argument{"invalid decay product parent"};
}
const int mother_idx = std::distance(begin(outgoing), mother);
assert(mother_idx >= 0);
decays[mother_idx].emplace_back(extract_particle(hepeup, i));
}
}
Event::Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double const min_jet_pt
):
Event( Event::EventData{
ev.incoming, ev.outgoing, ev.decays,
Parameters<EventParameters>{ev.central, ev.variations}
}.cluster(jet_def, min_jet_pt) )
{}
//! @TODO remove in HEJ 2.2.0
UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){
Event::EventData const evData{hepeup};
incoming = evData.incoming;
outgoing = evData.outgoing;
decays = evData.decays;
central = evData.parameters.central;
variations = evData.parameters.variations;
}
void Event::EventData::sort(){
// sort particles
std::sort(
begin(incoming), end(incoming),
[](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();}
);
auto old_outgoing = std::move(outgoing);
std::vector<size_t> idx(old_outgoing.size());
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(), [&old_outgoing](size_t i, size_t j){
return old_outgoing[i].rapidity() < old_outgoing[j].rapidity();
});
outgoing.clear();
outgoing.reserve(old_outgoing.size());
for(size_t i: idx) {
outgoing.emplace_back(std::move(old_outgoing[i]));
}
// find decays again
if(!decays.empty()){
auto old_decays = std::move(decays);
decays.clear();
for(size_t i=0; i<idx.size(); ++i) {
auto decay = old_decays.find(idx[i]);
if(decay != old_decays.end())
decays.emplace(i, std::move(decay->second));
}
assert(old_decays.size() == decays.size());
}
}
namespace {
Particle reconstruct_boson(std::vector<Particle> const & leptons) {
Particle decayed_boson;
decayed_boson.p = leptons[0].p + leptons[1].p;
const int pidsum = leptons[0].type + leptons[1].type;
if(pidsum == +1) {
assert(is_antilepton(leptons[0]));
if(is_antineutrino(leptons[0])) {
throw not_implemented{"lepton-flavour violating final state"};
}
assert(is_neutrino(leptons[1]));
// charged antilepton + neutrino means we had a W+
decayed_boson.type = pid::Wp;
}
else if(pidsum == -1) {
assert(is_antilepton(leptons[0]));
if(is_neutrino(leptons[1])) {
throw not_implemented{"lepton-flavour violating final state"};
}
assert(is_antineutrino(leptons[0]));
// charged lepton + antineutrino means we had a W-
decayed_boson.type = pid::Wm;
}
else {
throw not_implemented{
"final state with leptons "
+ name(leptons[0].type)
+ " and "
+ name(leptons[1].type)
};
}
return decayed_boson;
}
}
void Event::EventData::reconstruct_intermediate() {
const auto begin_leptons = std::partition(
begin(outgoing), end(outgoing),
[](Particle const & p) {return !is_anylepton(p);}
);
if(begin_leptons == end(outgoing)) return;
assert(is_anylepton(*begin_leptons));
std::vector<Particle> leptons(begin_leptons, end(outgoing));
outgoing.erase(begin_leptons, end(outgoing));
if(leptons.size() != 2) {
throw not_implemented{"Final states with one or more than two leptons"};
}
std::sort(
begin(leptons), end(leptons),
[](Particle const & p0, Particle const & p1) {
return p0.type < p1.type;
}
);
outgoing.emplace_back(reconstruct_boson(leptons));
decays.emplace(outgoing.size()-1, std::move(leptons));
}
Event Event::EventData::cluster(
fastjet::JetDefinition const & jet_def, double const min_jet_pt
){
sort();
Event ev{ std::move(incoming), std::move(outgoing), std::move(decays),
std::move(parameters),
jet_def, min_jet_pt
};
assert(std::is_sorted(begin(ev.outgoing_), end(ev.outgoing_),
rapidity_less{}));
ev.type_ = classify(ev);
return ev;
}
Event::Event(
std::array<Particle, 2> && incoming,
std::vector<Particle> && outgoing,
std::unordered_map<size_t, std::vector<Particle>> && decays,
Parameters<EventParameters> && parameters,
fastjet::JetDefinition const & jet_def,
double const min_jet_pt
): incoming_{std::move(incoming)},
outgoing_{std::move(outgoing)},
decays_{std::move(decays)},
parameters_{std::move(parameters)},
cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
min_jet_pt_{min_jet_pt}
{
jets_ = sorted_by_rapidity(cs_.inclusive_jets(min_jet_pt_));
}
namespace {
// check that Particles have a reasonable colour
bool correct_colour(Particle const & part){
ParticleID id{ part.type };
if(!is_parton(id))
return !part.colour;
if(!part.colour)
return false;
Colour const & col{ *part.colour };
if(is_quark(id))
return col.first != 0 && col.second == 0;
if(is_antiquark(id))
return col.first == 0 && col.second != 0;
assert(id==ParticleID::gluon);
return col.first != 0 && col.second != 0 && col.first != col.second;
}
}
bool Event::is_leading_colour() const {
if( !correct_colour(incoming()[0]) || !correct_colour(incoming()[1]) )
return false;
Colour line_colour = *incoming()[0].colour;
std::swap(line_colour.first, line_colour.second);
for(auto const & part: outgoing()){
// reasonable colour
if(!correct_colour(part))
return false;
if(!is_parton(part)) // skip colour neutral particles
continue;
// if possible connect to line
if( line_colour.first == part.colour->second )
line_colour.first = part.colour->first;
else if( line_colour.second == part.colour->first )
line_colour.second = part.colour->second;
else
return false;
// no colour singlet exchange/disconnected diagram
if(line_colour.first == line_colour.second)
return false;
}
return (incoming()[1].colour->first == line_colour.first)
&& (incoming()[1].colour->second == line_colour.second);
}
namespace {
void connect_incoming(Particle & in, int & colour, int & anti_colour){
in.colour = std::make_pair(anti_colour, colour);
// gluon
if(in.type == pid::gluon)
return;
if(in.type > 0){
// quark
assert(is_quark(in));
in.colour->second = 0;
colour*=-1;
return;
}
// anti-quark
assert(is_antiquark(in));
in.colour->first = 0;
anti_colour*=-1;
return;
}
}
bool Event::generate_colours(RNG & ran){
// generate only for HEJ events
if(!event_type::is_resummable(type()))
return false;
assert(std::is_sorted(
begin(outgoing()), end(outgoing()), rapidity_less{}));
assert(incoming()[0].pz() < incoming()[1].pz());
// positive (anti-)colour -> can connect
// negative (anti-)colour -> not available/used up by (anti-)quark
int colour = COLOUR_OFFSET;
int anti_colour = colour+1;
// initialise first
connect_incoming(incoming_[0], colour, anti_colour);
for(auto & part: outgoing_){
assert(colour>0 || anti_colour>0);
if(part.type == ParticleID::gluon){
// gluon
if(colour>0 && anti_colour>0){
// on g line => connect to colour OR anti-colour (random)
if(ran.flat() < 0.5){
part.colour = std::make_pair(colour+2,colour);
colour+=2;
} else {
part.colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(colour > 0){
// on q line => connect to available colour
part.colour = std::make_pair(colour+2, colour);
colour+=2;
} else {
assert(colour<0 && anti_colour>0);
// on qx line => connect to available anti-colour
part.colour = std::make_pair(anti_colour, anti_colour+2);
anti_colour+=2;
}
} else if(is_quark(part)) {
// quark
assert(anti_colour>0);
if(colour>0){
// on g line => connect and remove anti-colour
part.colour = std::make_pair(anti_colour, 0);
anti_colour+=2;
anti_colour*=-1;
} else {
// on qx line => new colour
colour*=-1;
part.colour = std::make_pair(colour, 0);
}
} else if(is_antiquark(part)) {
// anti-quark
assert(colour>0);
if(anti_colour>0){
// on g line => connect and remove colour
part.colour = std::make_pair(0, colour);
colour+=2;
colour*=-1;
} else {
// on q line => new anti-colour
anti_colour*=-1;
part.colour = std::make_pair(0, anti_colour);
}
} else { // not a parton
assert(!is_parton(part));
part.colour = {};
}
}
// Connect last
connect_incoming(incoming_[1], anti_colour, colour);
assert(is_leading_colour());
return true;
} // generate_colours
namespace {
bool valid_parton(
std::vector<fastjet::PseudoJet> const & jets,
Particle const & parton, int const idx,
double const max_ext_soft_pt_fraction, double const min_extparton_pt
){
// TODO code overlap with PhaseSpacePoint::pass_extremal_cuts
if(min_extparton_pt > parton.pt()) return false;
if(idx<0) return false;
- assert((int) jets.size()>=idx);
+ assert(static_cast<int>(jets.size())>=idx);
auto const & jet{ jets[idx] };
if( (parton.p - jet).pt()/jet.pt() > max_ext_soft_pt_fraction)
return false;
return true;
}
}
// this should work with multiple types
bool Event::valid_hej_state(double const max_frac,
double const min_pt
) const {
using namespace event_type;
if(!is_resummable(type()))
return false;
auto const & jet_idx{ particle_jet_indices() };
auto idx_begin{ jet_idx.cbegin() };
auto idx_end{ jet_idx.crbegin() };
auto part_begin{ cbegin_partons() };
auto part_end{ crbegin_partons() };
// always seperate extremal jets
if( !valid_parton(jets(), *part_begin, *idx_begin, max_frac, min_pt) )
return false;
++part_begin;
++idx_begin;
if( !valid_parton(jets(), *part_end, *idx_end, max_frac, min_pt) )
return false;
++part_end;
++idx_end;
// unob -> second parton in own jet
if( type() & (unob | qqxexb) ){
if( !valid_parton(jets(), *part_begin, *idx_begin, max_frac, min_pt) )
return false;
++part_begin;
++idx_begin;
}
if( type() & (unof | qqxexf) ){
if( !valid_parton(jets(), *part_end, *idx_end, max_frac, min_pt) )
return false;
++part_end;
++idx_end;
}
if( type() & qqxmid ){
// find qqx pair
auto begin_qqx{ std::find_if( part_begin, part_end.base(),
[](Particle const & part) -> bool {
return part.type != ParticleID::gluon;
}
)};
assert(begin_qqx != part_end.base());
long int qqx_pos{ std::distance(part_begin, begin_qqx) };
assert(qqx_pos >= 0);
idx_begin+=qqx_pos;
if( !( valid_parton(jets(),*begin_qqx, *idx_begin, max_frac,min_pt)
&& valid_parton(jets(),*(++begin_qqx),*(++idx_begin),max_frac,min_pt)
))
return false;
}
return true;
}
Event::ConstPartonIterator Event::begin_partons() const {
return cbegin_partons();
}
Event::ConstPartonIterator Event::cbegin_partons() const {
return boost::make_filter_iterator(
static_cast<bool (*)(Particle const &)>(is_parton),
cbegin(outgoing()),
cend(outgoing())
);
}
Event::ConstPartonIterator Event::end_partons() const {
return cend_partons();
}
Event::ConstPartonIterator Event::cend_partons() const {
return boost::make_filter_iterator(
static_cast<bool (*)(Particle const &)>(is_parton),
cend(outgoing()),
cend(outgoing())
);
}
Event::ConstReversePartonIterator Event::rbegin_partons() const {
return crbegin_partons();
}
Event::ConstReversePartonIterator Event::crbegin_partons() const {
return std::reverse_iterator<ConstPartonIterator>( cend_partons() );
}
Event::ConstReversePartonIterator Event::rend_partons() const {
return crend_partons();
}
Event::ConstReversePartonIterator Event::crend_partons() const {
return std::reverse_iterator<ConstPartonIterator>( cbegin_partons() );
}
namespace {
void print_momentum(std::ostream & os, fastjet::PseudoJet const & part){
const std::streamsize orig_prec = os.precision();
os <<std::scientific<<std::setprecision(6) << "["
<<std::setw(13)<<std::right<< part.px() << ", "
<<std::setw(13)<<std::right<< part.py() << ", "
<<std::setw(13)<<std::right<< part.pz() << ", "
<<std::setw(13)<<std::right<< part.E() << "]"<< std::fixed;
os.precision(orig_prec);
}
void print_colour(std::ostream & os, optional<Colour> const & col){
if(!col)
os << "(no color)"; // American spelling for better alignment
else
os << "(" <<std::setw(3)<<std::right<< col->first
<< ", " <<std::setw(3)<<std::right<< col->second << ")";
}
}
std::ostream& operator<<(std::ostream & os, Event const & ev){
const std::streamsize orig_prec = os.precision();
os <<std::setprecision(4)<<std::fixed;
os << "########## " << event_type::name(ev.type()) << " ##########" << std::endl;
os << "Incoming particles:\n";
for(auto const & in: ev.incoming()){
os <<std::setw(3)<< in.type << ": ";
print_colour(os, in.colour);
os << " ";
print_momentum(os, in.p);
os << std::endl;
}
os << "\nOutgoing particles: " << ev.outgoing().size() << "\n";
for(auto const & out: ev.outgoing()){
os <<std::setw(3)<< out.type << ": ";
print_colour(os, out.colour);
os << " ";
print_momentum(os, out.p);
os << " => rapidity="
<<std::setw(7)<<std::right<< out.rapidity() << std::endl;
}
os << "\nForming Jets: " << ev.jets().size() << "\n";
for(auto const & jet: ev.jets()){
print_momentum(os, jet);
os << " => rapidity="
<<std::setw(7)<<std::right<< jet.rapidity() << std::endl;
}
if(ev.decays().size() > 0 ){
os << "\nDecays: " << ev.decays().size() << "\n";
for(auto const & decay: ev.decays()){
os <<std::setw(3)<< ev.outgoing()[decay.first].type
<< " (" << decay.first << ") to:\n";
for(auto const & out: decay.second){
os <<" "<<std::setw(3)<< out.type << ": ";
print_momentum(os, out.p);
os << " => rapidity="
<<std::setw(7)<<std::right<< out.rapidity() << std::endl;
}
}
}
os << std::defaultfloat;
os.precision(orig_prec);
return os;
}
double shat(Event const & ev){
return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
}
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP * heprup){
LHEF::HEPEUP result;
result.heprup = heprup;
result.weights = {{event.central().weight, nullptr}};
for(auto const & var: event.variations()){
result.weights.emplace_back(var.weight, nullptr);
}
size_t num_particles = event.incoming().size() + event.outgoing().size();
for(auto const & decay: event.decays()) num_particles += decay.second.size();
result.NUP = num_particles;
// the following entries are pretty much meaningless
result.IDPRUP = event.type(); // event type
result.AQEDUP = 1./128.; // alpha_EW
//result.AQCDUP = 0.118 // alpha_QCD
// end meaningless part
result.XWGTUP = event.central().weight;
result.SCALUP = event.central().muf;
result.scales.muf = event.central().muf;
result.scales.mur = event.central().mur;
result.scales.SCALUP = event.central().muf;
result.pdfinfo.p1 = event.incoming().front().type;
result.pdfinfo.p2 = event.incoming().back().type;
result.pdfinfo.scale = event.central().muf;
result.IDUP.reserve(num_particles); // PID
result.ISTUP.reserve(num_particles); // status (in, out, decay)
result.PUP.reserve(num_particles); // momentum
result.MOTHUP.reserve(num_particles); // index mother particle
result.ICOLUP.reserve(num_particles); // colour
// incoming
std::array<Particle, 2> incoming{ event.incoming() };
// First incoming should be positive pz according to LHE standard
// (or at least most (everyone?) do it this way, and Pythia assumes it)
if(incoming[0].pz() < incoming[1].pz())
std::swap(incoming[0], incoming[1]);
for(Particle const & in: incoming){
result.IDUP.emplace_back(in.type);
result.ISTUP.emplace_back(status_in);
result.PUP.push_back({in.p[0], in.p[1], in.p[2], in.p[3], in.p.m()});
result.MOTHUP.emplace_back(0, 0);
assert(in.colour);
result.ICOLUP.emplace_back(*in.colour);
}
// outgoing
for(size_t i = 0; i < event.outgoing().size(); ++i){
Particle const & out = event.outgoing()[i];
result.IDUP.emplace_back(out.type);
const int status = event.decays().count(i)?status_decayed:status_out;
result.ISTUP.emplace_back(status);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
result.MOTHUP.emplace_back(1, 2);
if(out.colour)
result.ICOLUP.emplace_back(*out.colour);
else{
result.ICOLUP.emplace_back(std::make_pair(0,0));
}
}
// decays
for(auto const & decay: event.decays()){
- for(auto const out: decay.second){
+ for(auto const & out: decay.second){
result.IDUP.emplace_back(out.type);
result.ISTUP.emplace_back(status_out);
result.PUP.push_back({out.p[0], out.p[1], out.p[2], out.p[3], out.p.m()});
const size_t mother_idx = 1 + event.incoming().size() + decay.first;
result.MOTHUP.emplace_back(mother_idx, mother_idx);
result.ICOLUP.emplace_back(0,0);
}
}
assert(result.ICOLUP.size() == num_particles);
static constexpr double unknown_spin = 9.; //per Les Houches accord
result.VTIMUP = std::vector<double>(num_particles, unknown_spin);
result.SPINUP = result.VTIMUP;
return result;
}
}
diff --git a/src/EventReweighter.cc b/src/EventReweighter.cc
index bf28520..8911f98 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,284 +1,284 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/EventReweighter.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <string>
#include <unordered_map>
#include <utility>
#include "fastjet/ClusterSequence.hh"
#include "LHEF/LHEF.h"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/PhaseSpacePoint.hh"
#include "HEJ/RNG.hh"
namespace HEJ{
namespace {
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
Event::EventData to_EventData(PhaseSpacePoint const & psp){
Event::EventData result;
result.incoming=psp.incoming();
assert(result.incoming.size() == 2);
result.outgoing=psp.outgoing();
// technically Event::EventData doesn't have to be sorted,
// but PhaseSpacePoint should be anyway
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.decays = psp.decays();
result.parameters.central = {NaN, NaN, psp.weight()};
return result;
}
} // namespace anonymous
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
std::shared_ptr<RNG> ran
):
EventReweighter{
HEJ::Beam{
heprup.EBMUP.first,
{{
static_cast<HEJ::ParticleID>(heprup.IDBMUP.first),
static_cast<HEJ::ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
std::move(scale_gen),
std::move(conf),
std::move(ran)
}
{
if(heprup.EBMUP.second != E_beam_){
throw std::invalid_argument(
"asymmetric beam: " + std::to_string(E_beam_)
+ " ---> <--- " + std::to_string(heprup.EBMUP.second)
);
- };
+ }
if(heprup.PDFSUP.second != pdf_.id()){
throw std::invalid_argument(
"conflicting PDF ids: " + std::to_string(pdf_.id())
+ " vs. " + std::to_string(heprup.PDFSUP.second)
);
}
}
EventReweighter::EventReweighter(
Beam beam,
int pdf_id,
ScaleGenerator scale_gen,
EventReweighterConfig conf,
std::shared_ptr<RNG> ran
):
param_{std::move(conf)},
E_beam_{beam.E},
pdf_{pdf_id, beam.type.front(), beam.type.back()},
MEt2_{
[this](double mu){ return pdf_.Halphas(mu); },
param_.ME_config
},
scale_gen_{std::move(scale_gen)},
ran_{std::move(ran)}
{
assert(ran_);
}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
std::vector<Event> EventReweighter::reweight(
- Event const & input_ev, int num_events
+ Event const & input_ev, size_t num_events
){
auto res_events{ gen_res_events(input_ev, num_events) };
if(res_events.empty()) return {};
for(auto & event: res_events) event = scale_gen_(std::move(event));
return rescale(input_ev, std::move(res_events));
}
EventTreatment EventReweighter::treatment(EventType type) const {
return param_.treat.at(type);
}
std::vector<Event> EventReweighter::gen_res_events(
Event const & ev,
size_t phase_space_points
){
assert(ev.variations().empty());
status_.clear();
switch(treatment(ev.type())){
case EventTreatment::discard: {
status_.emplace_back(StatusCode::discard);
return {};
}
case EventTreatment::keep:
if(! jets_pass_resummation_cuts(ev)) {
status_.emplace_back(StatusCode::failed_resummation_cuts);
return {};
}
else {
status_.emplace_back(StatusCode::good);
return {ev};
}
default:;
}
const double Born_shat = shat(ev);
std::vector<Event> resummation_events;
status_.reserve(phase_space_points);
for(size_t psp_number = 0; psp_number < phase_space_points; ++psp_number){
PhaseSpacePoint psp{ev, param_.psp_config, *ran_};
status_.emplace_back(psp.status());
assert(psp.status() != StatusCode::unspecified);
if(psp.status() != StatusCode::good) continue;
assert(psp.weight() != 0.);
if(psp.incoming()[0].E() > E_beam_ || psp.incoming()[1].E() > E_beam_) {
status_.back() = StatusCode::too_much_energy;
continue;
}
resummation_events.emplace_back(
to_EventData( std::move(psp) ).cluster(
param_.jet_param().def, param_.jet_param().min_pt
)
);
auto & new_event = resummation_events.back();
assert( new_event.valid_hej_state(
param_.psp_config.max_ext_soft_pt_fraction,
param_.psp_config.min_extparton_pt ) );
if( new_event.type() != ev.type() )
throw std::logic_error{"Resummation Event does not match Born event"};
new_event.generate_colours(*ran_);
assert(new_event.variations().empty());
new_event.central().mur = ev.central().mur;
new_event.central().muf = ev.central().muf;
const double resum_shat = shat(new_event);
new_event.central().weight *= ev.central().weight*Born_shat*Born_shat/
(phase_space_points*resum_shat*resum_shat);
}
return resummation_events;
}
std::vector<Event> EventReweighter::rescale(
Event const & Born_ev,
std::vector<Event> events
) const{
const double Born_pdf = pdf_factors(Born_ev).central;
const double Born_ME = tree_matrix_element(Born_ev);
for(auto & cur_event: events){
const auto pdf = pdf_factors(cur_event);
assert(pdf.variations.size() == cur_event.variations().size());
const auto ME = matrix_elements(cur_event);
assert(ME.variations.size() == cur_event.variations().size());
cur_event.parameters() *= pdf*ME/(Born_pdf*Born_ME);
}
return events;
}
bool EventReweighter::jets_pass_resummation_cuts(
Event const & ev
) const{
const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
fastjet::ClusterSequence cs{out_as_PseudoJet, param_.jet_param().def};
return cs.inclusive_jets(param_.jet_param().min_pt).size() == ev.jets().size();
}
Weights EventReweighter::pdf_factors(Event const & ev) const{
auto const & a = ev.incoming().front();
auto const & b = ev.incoming().back();
const double xa = a.p.e()/E_beam_;
const double xb = b.p.e()/E_beam_;
Weights result;
std::unordered_map<double, double> known_pdf;
result.central =
pdf_.pdfpt(0,xa,ev.central().muf,a.type)*
pdf_.pdfpt(1,xb,ev.central().muf,b.type);
known_pdf.emplace(ev.central().muf, result.central);
result.variations.reserve(ev.variations().size());
for(auto const & ev_param: ev.variations()){
const double muf = ev_param.muf;
auto cur_pdf = known_pdf.find(muf);
if(cur_pdf == known_pdf.end()){
cur_pdf = known_pdf.emplace(
muf,
pdf_.pdfpt(0,xa,muf,a.type)*pdf_.pdfpt(1,xb,muf,b.type)
).first;
}
result.variations.emplace_back(cur_pdf->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
Weights
EventReweighter::matrix_elements(Event const & ev) const{
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev);
}
return MEt2_(ev);
}
double EventReweighter::tree_matrix_element(Event const & ev) const{
assert(ev.variations().empty());
assert(param_.treat.count(ev.type()) > 0);
if(param_.treat.find(ev.type())->second == EventTreatment::keep){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(ev).central;
}
Weights
EventReweighter::fixed_order_scale_ME(Event const & ev) const{
int alpha_s_power = 0;
for(auto const & part: ev.outgoing()){
if(is_parton(part))
++alpha_s_power;
else if(part.type == pid::Higgs) {
alpha_s_power += 2;
}
// nothing to do for other uncoloured particles
}
Weights result;
result.central = pow(pdf_.Halphas(ev.central().mur), alpha_s_power);
for(auto const & var: ev.variations()){
result.variations.emplace_back(
pow(pdf_.Halphas(var.mur), alpha_s_power)
);
}
return result;
}
} // namespace HEJ
diff --git a/src/HepMC2Interface.cc b/src/HepMC2Interface.cc
index 68e41c3..c1e619d 100644
--- a/src/HepMC2Interface.cc
+++ b/src/HepMC2Interface.cc
@@ -1,155 +1,155 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HepMC2Interface.hh"
#include "HEJ/exceptions.hh"
#ifdef HEJ_BUILD_WITH_HepMC2
#include <math.h>
#include <utility>
#include "HEJ/detail/HepMCInterface_common.hh"
#include "HEJ/Event.hh"
#include "HEJ/Particle.hh"
#include "LHEF/LHEF.h"
#include "HepMC/GenCrossSection.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
namespace HEJ{
namespace detail_HepMC {
template<>
struct HepMCVersion<2> {
using GenEvent = HepMC::GenEvent;
using Beam = std::array<HepMC::GenParticle*,2>;
};
template<>
auto make_particle_ptr<2> (
Particle const & sp, int status
) {
return new HepMC::GenParticle(
to_FourVector<HepMC::FourVector>(sp),
static_cast<int> (sp.type),
status
);
}
template<>
auto make_vx_ptr<2>() {
return new HepMC::GenVertex();
}
}
HepMC2Interface::HepMC2Interface(LHEF::HEPRUP const & heprup):
event_count_(0.), tot_weight_(0.), tot_weight2_(0.)
{
beam_particle_[0] = static_cast<ParticleID>(heprup.IDBMUP.first);
beam_particle_[1] = static_cast<ParticleID>(heprup.IDBMUP.second);
beam_energy_[0] = heprup.EBMUP.first;
beam_energy_[1] = heprup.EBMUP.second;
}
HepMC::GenCrossSection HepMC2Interface::cross_section() const {
HepMC::GenCrossSection xs;
xs.set_cross_section(tot_weight_, sqrt(tot_weight2_));
return xs;
}
HepMC::GenEvent HepMC2Interface::init_event(Event const & event) const {
const std::array<HepMC::GenParticle*,2> beam {
new HepMC::GenParticle(
HepMC::FourVector(0,0,-beam_energy_[0],beam_energy_[0]),
beam_particle_[0], detail_HepMC::status_beam ),
new HepMC::GenParticle(
HepMC::FourVector(0,0, beam_energy_[1],beam_energy_[1]),
beam_particle_[1], detail_HepMC::status_beam )
};
auto hepmc_ev{ detail_HepMC::HepMC_init_kinematics<2>(
event, beam, HepMC::GenEvent{ HepMC::Units::GEV, HepMC::Units::MM }
) };
hepmc_ev.weights().push_back( event.central().weight );
for(auto const & var: event.variations()){
hepmc_ev.weights().push_back( var.weight );
// no weight name for HepMC2 since rivet3 seem to mix them up
// (could be added via hepmc_ev.weights()[name]=weight)
}
return hepmc_ev;
}
void HepMC2Interface::set_central(
HepMC::GenEvent & out_ev, Event const & event, ssize_t const weight_index
){
EventParameters event_param;
if(weight_index < 0)
event_param = event.central();
- else if ( (size_t) weight_index < event.variations().size())
+ else if ( static_cast<size_t>(weight_index) < event.variations().size())
event_param = event.variations(weight_index);
else
throw std::invalid_argument{
"HepMC2Interface tried to access a weight outside of the variation range."
};
const double wt = event_param.weight;
tot_weight_ += wt;
tot_weight2_ += wt * wt;
++event_count_;
// central always on first
assert(out_ev.weights().size() == event.variations().size()+1);
out_ev.weights()[0] = wt;
out_ev.set_cross_section( cross_section() );
out_ev.set_signal_process_id(event.type());
out_ev.set_event_scale(event_param.mur);
out_ev.set_event_number(event_count_);
/// @TODO add alphaQCD (need function) and alphaQED
/// @TODO output pdf (currently not avaiable from event alone)
}
HepMC::GenEvent HepMC2Interface::operator()(Event const & event,
ssize_t const weight_index
){
HepMC::GenEvent out_ev(init_event(event));
set_central(out_ev, event, weight_index);
return out_ev;
}
}
#else // no HepMC2 => empty class
namespace HepMC {
class GenEvent {};
class GenCrossSection {};
}
namespace HEJ{
HepMC2Interface::HepMC2Interface(LHEF::HEPRUP const &){
throw std::invalid_argument(
"Failed to create HepMC2Interface: "
"HEJ 2 was built without HepMC2 support"
);
}
HepMC::GenEvent HepMC2Interface::operator()(Event const &, ssize_t)
{return HepMC::GenEvent();}
HepMC::GenEvent HepMC2Interface::init_event(Event const &) const
{return HepMC::GenEvent();}
void HepMC2Interface::set_central(HepMC::GenEvent &, Event const &, ssize_t){}
HepMC::GenCrossSection HepMC2Interface::cross_section() const
{return HepMC::GenCrossSection();}
}
#endif
namespace HEJ{
HepMC2Interface::~HepMC2Interface() = default;
}
diff --git a/src/HepMC3Interface.cc b/src/HepMC3Interface.cc
index 06b20e0..a13f6c7 100644
--- a/src/HepMC3Interface.cc
+++ b/src/HepMC3Interface.cc
@@ -1,202 +1,202 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/HepMC3Interface.hh"
#include "HEJ/exceptions.hh"
#ifdef HEJ_BUILD_WITH_HepMC3
#include <math.h>
#include <utility>
#include "HEJ/detail/HepMCInterface_common.hh"
#include "HEJ/Event.hh"
#include "HEJ/Particle.hh"
#include "LHEF/LHEF.h"
#include "HepMC3/GenCrossSection.h"
#include "HepMC3/GenEvent.h"
#include "HepMC3/GenParticle.h"
#include "HepMC3/GenRunInfo.h"
#include "HepMC3/GenVertex.h"
#include "HepMC3/LHEFAttributes.h"
namespace HEJ{
namespace detail_HepMC {
template<>
struct HepMCVersion<3> {
using GenEvent = HepMC3::GenEvent;
using Beam = std::array<HepMC3::GenParticlePtr,2>;
};
template<>
auto make_particle_ptr<3> (
Particle const & sp, int status
) {
return HepMC3::make_shared<HepMC3::GenParticle>(
to_FourVector<HepMC3::FourVector>(sp),
static_cast<int> (sp.type),
status
);
}
template<>
auto make_vx_ptr<3>() {
return HepMC3::make_shared<HepMC3::GenVertex>();
}
}
namespace {
void reset_weight_info(LHEF::HEPRUP & heprup){
heprup.IDWTUP = 2;
// use placeholders for unknown init block values
// we can overwrite them after processing all events
heprup.XSECUP = {0.};
heprup.XERRUP = {0.};
heprup.XMAXUP = {0.};
}
HepMC3::shared_ptr<HepMC3::GenRunInfo> init_runinfo(LHEF::HEPRUP heprup){
reset_weight_info(heprup);
auto runinfo{ HepMC3::make_shared<HepMC3::GenRunInfo>() };
auto hepr{ HepMC3::make_shared<HepMC3::HEPRUPAttribute>() };
hepr->heprup = heprup;
runinfo->add_attribute(std::string("HEPRUP"), hepr);
for(auto const & gen: heprup.generators){
runinfo->tools().emplace_back(
HepMC3::GenRunInfo::ToolInfo{gen.name, gen.version, gen.contents} );
}
return runinfo;
}
std::vector<std::string> get_weight_names(Event const & ev){
std::vector<std::string> names;
names.reserve(ev.variations().size()+1); // +1 from central
names.emplace_back(""); // rivet assumes central band to have no name
for( size_t i=0; i<ev.variations().size(); ++i ){
auto const & var{ ev.variations()[i] };
if(var.description){
names.emplace_back( to_simple_string(*var.description) );
} else {
names.emplace_back( "" );
}
}
assert(names.size() == ev.variations().size()+1);
return names;
}
} // namespace anonymous
HepMC3Interface::HepMC3Interface(LHEF::HEPRUP const & heprup):
run_info{ init_runinfo(heprup) },
event_count_(0.), tot_weight_(0.), tot_weight2_(0.),
xs_{std::make_shared<HepMC3::GenCrossSection>()}
{
beam_particle_[0] = static_cast<ParticleID>(heprup.IDBMUP.first);
beam_particle_[1] = static_cast<ParticleID>(heprup.IDBMUP.second);
beam_energy_[0] = heprup.EBMUP.first;
beam_energy_[1] = heprup.EBMUP.second;
}
HepMC3::GenEvent HepMC3Interface::init_event(Event const & event) const {
const std::array<HepMC3::GenParticlePtr,2> beam {
HepMC3::make_shared<HepMC3::GenParticle>(
HepMC3::FourVector(0,0,-beam_energy_[0],beam_energy_[0]),
beam_particle_[0], detail_HepMC::status_beam ),
HepMC3::make_shared<HepMC3::GenParticle>(
HepMC3::FourVector(0,0, beam_energy_[1],beam_energy_[1]),
beam_particle_[1], detail_HepMC::status_beam )
};
auto hepmc_ev{ detail_HepMC::HepMC_init_kinematics<3>(
event, beam, HepMC3::GenEvent{ HepMC3::Units::GEV, HepMC3::Units::MM }
) };
// set up run specific informations
if( run_info->weight_names().size() != event.variations().size()+1 ){
run_info->set_weight_names( get_weight_names(event) );
}
// order matters: weights in hepmc_ev initialised when registering run_info
hepmc_ev.set_run_info(run_info);
assert(hepmc_ev.weights().size() == event.variations().size()+1);
for(size_t i=0; i<event.variations().size(); ++i){
hepmc_ev.weights()[i+1] = event.variations()[i].weight;
//! @TODO set variation specific cross section
//! the problem is that set_cross_section overwrites everything
}
return hepmc_ev;
}
void HepMC3Interface::set_central(HepMC3::GenEvent & out_ev, Event const & event,
ssize_t const weight_index
){
EventParameters event_param;
if(weight_index < 0)
event_param = event.central();
- else if ( (size_t) weight_index < event.variations().size())
+ else if ( static_cast<size_t>(weight_index) < event.variations().size())
event_param = event.variations(weight_index);
else
throw std::invalid_argument{
"HepMC3Interface tried to access a weight outside of the variation range."
};
const double wt = event_param.weight;
tot_weight_ += wt;
tot_weight2_ += wt * wt;
++event_count_;
// central always on first
assert(out_ev.weights().size() == event.variations().size()+1);
out_ev.weights()[0] = wt;
// out_ev can be setup with a different central scale -> save xs manually
out_ev.set_cross_section(xs_);
assert(out_ev.cross_section() && out_ev.cross_section() == xs_);
// overwrites all previously set xs ...
xs_->set_cross_section(tot_weight_,sqrt(tot_weight2_));
out_ev.set_event_number(event_count_);
/// @TODO add number of attempted events
xs_->set_accepted_events(event_count_);
/// @TODO add alphaQCD (need function) and alphaQED
/// @TODO output pdf (currently not avaiable from event alone)
}
HepMC3::GenEvent HepMC3Interface::operator()(Event const & event,
ssize_t const weight_index
){
HepMC3::GenEvent out_ev(init_event(event));
set_central(out_ev, event, weight_index);
return out_ev;
}
}
#else // no HepMC3 => empty class
namespace HepMC3 {
class GenEvent {};
class GenCrossSection {};
class GenRunInfo {};
}
namespace HEJ{
HepMC3Interface::HepMC3Interface(LHEF::HEPRUP const &){
throw std::invalid_argument(
"Failed to create HepMC3Interface: "
"HEJ 2 was built without HepMC3 support"
);
}
HepMC3::GenEvent HepMC3Interface::operator()(Event const &, ssize_t)
{return HepMC3::GenEvent();}
HepMC3::GenEvent HepMC3Interface::init_event(Event const &) const
{return HepMC3::GenEvent();}
void HepMC3Interface::set_central(HepMC3::GenEvent &, Event const &, ssize_t){}
}
#endif
namespace HEJ{
HepMC3Interface::~HepMC3Interface() = default;
}
diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc
index b95b3d1..0e8d2af 100644
--- a/src/MatrixElement.cc
+++ b/src/MatrixElement.cc
@@ -1,1711 +1,1711 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/MatrixElement.hh"
#include <algorithm>
#include <assert.h>
#include <limits>
#include <math.h>
#include <stddef.h>
#include <unordered_map>
#include <utility>
#include "CLHEP/Vector/LorentzVector.h"
#include "HEJ/Constants.hh"
#include "HEJ/Wjets.hh"
#include "HEJ/Hjets.hh"
#include "HEJ/jets.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#include "HEJ/Particle.hh"
#include "HEJ/utility.hh"
namespace HEJ{
double MatrixElement::omega0(
double alpha_s, double mur,
fastjet::PseudoJet const & q_j
) const {
const double lambda = param_.regulator_lambda;
const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda));
if(! param_.log_correction) return result;
return (
1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda))
)*result;
}
Weights MatrixElement::operator()(Event const & event) const {
return tree(event)*virtual_corrections(event);
}
Weights MatrixElement::tree(Event const & event) const {
return tree_param(event)*tree_kin(event);
}
Weights MatrixElement::tree_param(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = tree_param(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = tree_param(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
Weights MatrixElement::virtual_corrections(Event const & event) const {
if(! is_resummable(event.type())) {
return Weights{0., std::vector<double>(event.variations().size(), 0.)};
}
Weights result;
// only compute once for each renormalisation scale
std::unordered_map<double, double> known;
result.central = virtual_corrections(event, event.central().mur);
known.emplace(event.central().mur, result.central);
for(auto const & var: event.variations()) {
const auto ME_it = known.find(var.mur);
if(ME_it == end(known)) {
const double wt = virtual_corrections(event, var.mur);
result.variations.emplace_back(wt);
known.emplace(var.mur, wt);
}
else {
result.variations.emplace_back(ME_it->second);
}
}
return result;
}
double MatrixElement::virtual_corrections_W(
Event const & event,
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;
size_t first_idx = 0;
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 (abs(partons[0].type) != 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
}
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(partons) + 1, end(partons) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon); }
);
if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0;
if(abs(backquark->type) != abs((backquark+1)->type)) {
wqq=true;
#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(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -=partons[j+1].p;
} // End Loop one
if (last_idx != first_idx_qqx) q -= partons[last_idx+1].p;
if (wqq) q -= WBoson.p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
partons[j+1].rapidity() - partons[j].rapidity()
);
q -= partons[j+1].p;
}
#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 exp(exponent);
}
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) && abs(AWZH_boson->type) == pid::Wp){
return virtual_corrections_W(event, mur, *AWZH_boson);
}
assert(std::is_sorted(out.begin(), out.end(), rapidity_less{}));
assert(out.size() >= 2);
assert(pa.pz() < pb.pz());
fastjet::PseudoJet q = pa - out[0].p;
size_t first_idx = 0;
size_t last_idx = out.size() - 1;
// if there is a Higgs boson, extremal qqx or unordered gluon
// outside the extremal partons then it is not part of the FKL
// ladder and does not contribute to the virtual corrections
if((out.front().type == pid::Higgs)
|| event.type() == event_type::unob
|| event.type() == event_type::qqxexb){
q -= out[1].p;
++first_idx;
}
if((out.back().type == pid::Higgs)
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf){
--last_idx;
}
size_t first_idx_qqx = last_idx;
size_t last_idx_qqx = last_idx;
//if qqxMid event, virtual correction do not occur between
//qqx pair.
if(event.type() == event_type::qqxmid){
const auto backquark = std::find_if(
begin(out) + 1, end(out) - 1 ,
[](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); }
);
if(backquark == end(out) || (backquark+1)->type==pid::gluon) return 0;
last_idx = std::distance(begin(out), backquark);
first_idx_qqx = last_idx+1;
}
double exponent = 0;
const double alpha_s = alpha_s_(mur);
for(size_t j = first_idx; j < last_idx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
if (last_idx != first_idx_qqx) q -= out[last_idx+1].p;
for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){
exponent += omega0(alpha_s, mur, q)*(
out[j+1].rapidity() - out[j].rapidity()
);
q -= out[j+1].p;
}
assert(
nearby(q, -1*pb, norm)
|| out.back().type == pid::Higgs
|| event.type() == event_type::unof
|| event.type() == event_type::qqxexf
);
return exp(exponent);
}
namespace {
//! Lipatov vertex for partons emitted into extremal jets
double C2Lipatov(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2
){
const CLHEP::HepLorentzVector temptrans=-(qav+qbv);
const CLHEP::HepLorentzVector p5=qav-qbv;
const CLHEP::HepLorentzVector CL=temptrans
+ p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2))
- p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2));
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction for partons emitted into extremal jets
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
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);
}
//! 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 temptrans=-(qav+qbv);
const CLHEP::HepLorentzVector p5=qav-qbv;
const CLHEP::HepLorentzVector CL=temptrans
+ qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2.
- qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2.
+ ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom))
+ pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom))
- pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim))
- pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.;
return -CL.dot(CL);
}
//! Lipatov vertex with soft subtraction
double C2Lipatovots(
CLHEP::HepLorentzVector const & qav,
CLHEP::HepLorentzVector const & qbv,
CLHEP::HepLorentzVector const & pa,
CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
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);
}
/** 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
){
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
){
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<HLV> const & partons){
// 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)/HEJ::C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/HEJ::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
){
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
){
// 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
){
// 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
){
// CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal)
const double CFbackward = K_g( (swap_q_qx)?pq:pqbar ,pa)/HEJ::C_F;
// With qqbar we could have 2 incoming gluons and W Emission
if (aptype==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<HLV> const & partons,
CLHEP::HepLorentzVector const & plbar,
CLHEP::HepLorentzVector const & pl,
bool const wqq, bool const wc,
ParticleProperties const & Wprop
){
// CAM factors for the qqx amps, and qqbar ordering (default, pq backwards)
const bool swap_q_qx=pqbar.rapidity() < pq.rapidity();
double wt=1.;
if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/HEJ::C_F;
if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/HEJ::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);
}
/** \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
){
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_qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
else
return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.);
}
}
throw std::logic_error("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
){
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(HEJ::Particle const & particle){
return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()};
}
void validate(HEJ::MatrixElementConfig const & config) {
#ifndef HEJ_BUILD_WITH_QCDLOOP
if(!config.Higgs_coupling.use_impact_factors) {
throw std::invalid_argument{
"Invalid Higgs coupling settings.\n"
"HEJ without QCDloop support can only use impact factors.\n"
"Set use_impact_factors to true or recompile HEJ.\n"
};
}
#endif
if(config.Higgs_coupling.use_impact_factors
&& config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) {
throw std::invalid_argument{
"Conflicting settings: "
"impact factors may only be used in the infinite top mass limit"
};
}
}
} // namespace anonymous
MatrixElement::MatrixElement(
std::function<double (double)> alpha_s,
MatrixElementConfig conf
):
alpha_s_{std::move(alpha_s)},
param_{std::move(conf)}
{
validate(param_);
}
double MatrixElement::tree_kin(
Event const & ev
) const {
if(! is_resummable(ev.type())) return 0.;
auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing()))
return tree_kin_jets(ev);
switch(AWZH_boson->type){
case pid::Higgs:
return tree_kin_Higgs(ev);
case pid::Wp:
case pid::Wm:
return tree_kin_W(ev);
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
namespace{
constexpr int extremal_jet_idx = 1;
constexpr int no_extremal_jet_idx = 0;
bool treat_as_extremal(Particle const & parton){
return parton.p.user_index() == extremal_jet_idx;
}
template<class InputIterator>
double FKL_ladder_weight(
InputIterator begin_gluon, InputIterator end_gluon,
CLHEP::HepLorentzVector const & q0,
CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb,
CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn,
double lambda
){
double wt = 1;
auto qi = q0;
for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){
assert(gluon_it->type == pid::gluon);
const auto g = to_HepLorentzVector(*gluon_it);
const auto qip1 = qi - g;
if(treat_as_extremal(*gluon_it)){
wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A;
} else{
wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A;
}
qi = qip1;
}
return wt;
}
} // namespace anonymous
std::vector<Particle> MatrixElement::tag_extremal_jet_partons(
Event const & ev
) const{
auto out_partons = filter_partons(ev.outgoing());
if(out_partons.size() == ev.jets().size()){
// no additional emissions in extremal jets, don't need to tag anything
for(auto & parton: out_partons){
parton.p.set_user_index(no_extremal_jet_idx);
}
return out_partons;
}
const auto & jets = ev.jets();
assert(jets.size() >= 2);
auto most_backward = begin(jets);
auto most_forward = end(jets) - 1;
// skip jets caused by unordered emission or qqx
if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){
assert(jets.size() >= 3);
++most_backward;
}
else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){
assert(jets.size() >= 3);
--most_forward;
}
const auto extremal_jet_indices = ev.particle_jet_indices(
{*most_backward, *most_forward}
);
assert(extremal_jet_indices.size() == out_partons.size());
for(size_t i = 0; i < out_partons.size(); ++i){
assert(HEJ::is_parton(out_partons[i]));
const int idx = (extremal_jet_indices[i]>=0)?
extremal_jet_idx:
no_extremal_jet_idx;
out_partons[i].p.set_user_index(idx);
}
return out_partons;
}
namespace {
double tree_kin_jets_qqxmid(
ParticleID aptype, ParticleID bptype, HLV pa, HLV pb,
std::vector<Particle> const & partons,
double lambda
){
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons[0]);
auto pn = to_HepLorentzVector(partons[partons.size() - 1]);
auto q0 = pa - p1;
// t-channel momentum after qqx
auto qqxt = q0;
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<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_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;
}
}
double MatrixElement::tree_kin_jets(Event const & ev) const {
auto const & incoming = ev.incoming();
const auto partons = tag_extremal_jet_partons(ev);
if (ev.type()==HEJ::event_type::FKL){
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
return ME_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa
)/(4.*(N_C*N_C - 1.))*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
pa - p1, pa, pb, p1, pn,
param_.regulator_lambda
);
}
else if (ev.type()==HEJ::event_type::unordered_backward){
return tree_kin_jets_uno(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::unordered_forward){
return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::extremal_qqxb){
return tree_kin_jets_qqx(incoming.begin(), incoming.end(),
partons.begin(), partons.end(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::event_type::extremal_qqxf){
return tree_kin_jets_qqx(incoming.rbegin(), incoming.rend(),
partons.rbegin(), partons.rend(),
param_.regulator_lambda);
}
else if (ev.type()==HEJ::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, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV 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 HLV & plbar, const HLV & 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 HLV & plbar, const HLV & 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, HLV pa, HLV pb,
std::vector<Particle> const & partons,
HLV plbar, HLV pl,
double lambda, ParticleProperties const & Wprop
){
HLV pq,pqbar;
const auto backmidquark = std::find_if(
begin(partons)+1, end(partons)-1,
[](Particle const & s){ return s.type != pid::gluon; }
);
assert(backmidquark!=end(partons)-1);
if (is_quark(backmidquark->type)){
pq = to_HepLorentzVector(*backmidquark);
pqbar = to_HepLorentzVector(*(backmidquark+1));
}
else {
pqbar = to_HepLorentzVector(*backmidquark);
pq = to_HepLorentzVector(*(backmidquark+1));
}
auto p1 = to_HepLorentzVector(partons.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 || (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<HLV> partonsHLV;
partonsHLV.reserve(partons.size());
for (size_t i = 0; i != partons.size(); ++i) {
partonsHLV.push_back(to_HepLorentzVector(partons[i]));
}
const double current_factor = ME_W_qqxmid_current(
aptype, bptype, nabove, nbelow, pa, pb,
pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop
);
const double ladder_factor = FKL_ladder_weight(
begin_ladder, end_ladder_1,
q0, pa, pb, p1, pn,
lambda
)*FKL_ladder_weight(
begin_ladder_2, end_ladder,
qqxt, pa, pb, p1, pn,
lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
} // namespace anonymous
double MatrixElement::tree_kin_W(Event const & ev) const {
using namespace event_type;
auto const & incoming(ev.incoming());
#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( (long int) ev.decays().cbegin()->first
+ 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
HLV 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());
}
double MatrixElement::tree_kin_Higgs(Event const & ev) const {
if(is_uno(ev.type())){
return tree_kin_Higgs_between(ev);
}
if(ev.outgoing().front().type == pid::Higgs){
return tree_kin_Higgs_first(ev);
}
if(ev.outgoing().back().type == pid::Higgs){
return tree_kin_Higgs_last(ev);
}
return tree_kin_Higgs_between(ev);
}
namespace {
// Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113
#ifdef HEJ_BUILD_WITH_QCDLOOP
// TODO: code duplication with jets.cc
double K_g(double p1minus, double paminus) {
return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A;
}
double K_g(
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
return K_g(pout.minus(), pin.minus());
}
double K(
ParticleID type,
CLHEP::HepLorentzVector const & pout,
CLHEP::HepLorentzVector const & pin
) {
if(type == pid::gluon) return 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;
}
}
double MatrixElement::MH2_forwardH(
CLHEP::HepLorentzVector const & p1out,
CLHEP::HepLorentzVector const & p1in,
ParticleID type2,
CLHEP::HepLorentzVector const & p2out,
CLHEP::HepLorentzVector const & p2in,
CLHEP::HepLorentzVector const & pH,
double t1, double t2
) const{
ignore(p2out, p2in);
const double shat = p1in.invariantMass2(p2in);
const double vev = param_.ew_parameters.vev();
// gluon case
#ifdef HEJ_BUILD_WITH_QCDLOOP
if(!param_.Higgs_coupling.use_impact_factors){
return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*ME_Houtside_gq(
p1out, p1in, p2out, p2in, pH,
param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom,
param_.Higgs_coupling.mb, vev
)/(4*(N_C*N_C - 1));
}
#endif
return K_MRK(type2)/C_A*9./2.*shat*shat*(
C2gHgp(p1in,p1out,pH,vev) + C2gHgm(p1in,p1out,pH,vev)
)/(t1*t2);
}
double MatrixElement::tree_kin_Higgs_first(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.front().type == pid::Higgs);
if(outgoing[1].type != pid::gluon) {
assert(incoming.front().type == outgoing[1].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.front());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
const auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
const auto q0 = pa - p1 - pH;
const double t1 = q0.m2();
const double t2 = (pn - pb).m2();
return MH2_forwardH(
p1, pa, incoming[1].type, pn, pb, pH,
t1, t2
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
double MatrixElement::tree_kin_Higgs_last(Event const & ev) const {
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
assert(outgoing.back().type == pid::Higgs);
if(outgoing[outgoing.size()-2].type != pid::gluon) {
assert(incoming.back().type == outgoing[outgoing.size()-2].type);
return tree_kin_Higgs_between(ev);
}
const auto pH = to_HepLorentzVector(outgoing.back());
const auto partons = tag_extremal_jet_partons(
ev
);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(partons.front());
const auto pn = to_HepLorentzVector(partons.back());
auto q0 = pa - p1;
const double t1 = q0.m2();
const double t2 = (pn + pH - pb).m2();
return MH2_forwardH(
pn, pb, incoming[0].type, p1, pa, pH,
t2, t1
)*FKL_ladder_weight(
begin(partons) + 1, end(partons) - 1,
q0, pa, pb, p1, pn,
param_.regulator_lambda
);
}
namespace {
template<class InIter, class partIter>
double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart,
partIter EndPart, const HLV & qH, const HLV & qHp1,
double mt, bool inc_bot, double mb, double vev){
const auto pa = to_HepLorentzVector(*BeginIn);
const auto pb = to_HepLorentzVector(*(EndIn-1));
const auto pg = to_HepLorentzVector(*BeginPart);
const auto p1 = to_HepLorentzVector(*(BeginPart+1));
const auto pn = to_HepLorentzVector(*(EndPart-1));
return ME_Higgs_current_uno(
(BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa,
qH, qHp1, mt, inc_bot, mb, vev
);
}
}
double MatrixElement::tree_kin_Higgs_between(Event const & ev) const {
using namespace event_type;
auto const & incoming = ev.incoming();
auto const & outgoing = ev.outgoing();
const auto the_Higgs = std::find_if(
begin(outgoing), end(outgoing),
[](Particle const & s){ return s.type == pid::Higgs; }
);
assert(the_Higgs != end(outgoing));
const auto pH = to_HepLorentzVector(*the_Higgs);
const auto partons = tag_extremal_jet_partons(ev);
const auto pa = to_HepLorentzVector(incoming[0]);
const auto pb = to_HepLorentzVector(incoming[1]);
auto p1 = to_HepLorentzVector(
partons[(ev.type() == unob)?1:0]
);
auto pn = to_HepLorentzVector(
partons[partons.size() - ((ev.type() == unof)?2:1)]
);
auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing));
assert(
(first_after_Higgs == end(partons) && (
(ev.type() == unob)
|| partons.back().type != pid::gluon
))
|| first_after_Higgs->rapidity() >= the_Higgs->rapidity()
);
assert(
(first_after_Higgs == begin(partons) && (
(ev.type() == unof)
|| partons.front().type != pid::gluon
))
|| (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity()
);
// always treat the Higgs as if it were in between the extremal FKL partons
if(first_after_Higgs == begin(partons)) ++first_after_Higgs;
else if(first_after_Higgs == end(partons)) --first_after_Higgs;
// t-channel momentum before Higgs
auto qH = pa;
for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){
qH -= to_HepLorentzVector(*parton_it);
}
auto q0 = pa - p1;
auto begin_ladder = begin(partons) + 1;
auto end_ladder = end(partons) - 1;
double current_factor;
if(ev.type() == FKL){
current_factor = ME_Higgs_current(
incoming[0].type, incoming[1].type,
pn, pb, p1, pa, qH, qH - pH,
param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
}
else if(ev.type() == unob){
current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
begin(incoming), end(incoming), begin(partons),
end(partons), qH, qH-pH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
const auto p_unob = to_HepLorentzVector(partons.front());
q0 -= p_unob;
p1 += p_unob;
++begin_ladder;
}
else if(ev.type() == unof){
current_factor = HEJ::C_A*HEJ::C_A/2*tree_kin_Higgs_uno(
rbegin(incoming), rend(incoming), rbegin(partons),
rend(partons), qH-pH, qH, param_.Higgs_coupling.mt,
param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb,
param_.ew_parameters.vev()
);
pn += to_HepLorentzVector(partons.back());
--end_ladder;
}
else{
throw std::logic_error("Can only reweight FKL or uno processes in H+Jets");
}
const double ladder_factor = FKL_ladder_weight(
begin_ladder, first_after_Higgs,
q0, pa, pb, p1, pn,
param_.regulator_lambda
)*FKL_ladder_weight(
first_after_Higgs, end_ladder,
qH - pH, pa, pb, p1, pn,
param_.regulator_lambda
);
return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor;
}
namespace {
double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) {
const auto AWZH_boson = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](auto const & p){return is_AWZH_boson(p);}
);
if(AWZH_boson == end(ev.outgoing())) return 1.;
switch(AWZH_boson->type){
case pid::Higgs:
return alpha_s*alpha_s;
case pid::Wp:
case pid::Wm:
return alpha_w*alpha_w;
// TODO
case pid::photon:
case pid::Z:
default:
throw not_implemented("Emission of boson of unsupported type");
}
}
}
double MatrixElement::tree_param(Event const & ev, double mur) const {
assert(is_resummable(ev.type()));
const auto begin_partons = ev.begin_partons();
const auto end_partons = ev.end_partons();
const auto num_partons = std::distance(begin_partons, end_partons);
const double alpha_s = alpha_s_(mur);
const double gs2 = 4.*M_PI*alpha_s;
double res = std::pow(gs2, num_partons);
if(param_.log_correction){
// use alpha_s(q_perp), evolved to mur
assert(num_partons >= 2);
const auto first_emission = std::next(begin_partons);
const auto last_emission = std::prev(end_partons);
for(auto parton = first_emission; parton != last_emission; ++parton){
res *= 1. + alpha_s/(2.*M_PI)*beta0*log(mur/parton->perp());
}
}
return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res;
}
} // namespace HEJ
diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
index 9c5b669..73c11d0 100644
--- a/src/PhaseSpacePoint.cc
+++ b/src/PhaseSpacePoint.cc
@@ -1,814 +1,814 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/PhaseSpacePoint.hh"
#include <algorithm>
#include <assert.h>
#include <numeric>
#include <random>
#include "fastjet/ClusterSequence.hh"
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/JetSplitter.hh"
#include "HEJ/kinematics.hh"
#include "HEJ/resummation_jet.hh"
#include "HEJ/utility.hh"
#include "HEJ/PDG_codes.hh"
#include "HEJ/event_types.hh"
namespace HEJ{
namespace {
constexpr int max_jet_user_idx = PhaseSpacePoint::ng_max;
bool is_nonjet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() > max_jet_user_idx;
}
bool is_jet_parton(fastjet::PseudoJet const & parton){
assert(parton.user_index() != -1);
return parton.user_index() <= max_jet_user_idx;
}
// user indices for partons with extremal rapidity
- constexpr int qqxmid1_idx = -9;
- constexpr int qqxmid2_idx = -8;
- constexpr int qqxb_idx = -7;
- constexpr int qqxf_idx = -6;
- constexpr int unob_idx = -5;
- constexpr int unof_idx = -4;
- constexpr int backward_FKL_idx = -3;
- constexpr int forward_FKL_idx = -2;
+ constexpr int qqxmid1_uid = -9;
+ constexpr int qqxmid2_uid = -8;
+ constexpr int qqxb_uid = -7;
+ constexpr int qqxf_uid = -6;
+ constexpr int unob_uid = -5;
+ constexpr int unof_uid = -4;
+ constexpr int backward_FKL_uid = -3;
+ constexpr int forward_FKL_uid = -2;
}
namespace {
double estimate_ng_mean(std::vector<fastjet::PseudoJet> const & Born_jets){
const double delta_y =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
assert(delta_y > 0);
// Formula derived from fit in arXiv:1805.04446 (see Fig. 2)
return 0.975052*delta_y;
}
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
return sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
}
bool PhaseSpacePoint::pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cluster_jets(jets).size() == jets.size();
}
int PhaseSpacePoint::sample_ng(
std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran
){
const double ng_mean = estimate_ng_mean(Born_jets);
std::poisson_distribution<int> dist(ng_mean);
const int ng = dist(ran);
assert(ng >= 0);
assert(ng < ng_max);
weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
return ng;
}
void PhaseSpacePoint::copy_AWZH_boson_from(Event const & event){
auto const & from = event.outgoing();
const auto AWZH_boson = std::find_if(
begin(from), end(from),
[](Particle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson == end(from)) return;
auto insertion_point = std::lower_bound(
begin(outgoing_), end(outgoing_), *AWZH_boson, rapidity_less{}
);
outgoing_.insert(insertion_point, *AWZH_boson);
// copy decay products
const int idx = std::distance(begin(from), AWZH_boson);
assert(idx >= 0);
const auto decay_it = event.decays().find(idx);
if(decay_it != end(event.decays())){
const int new_idx = std::distance(begin(outgoing_), insertion_point);
assert(new_idx >= 0);
assert(outgoing_[new_idx].type == AWZH_boson->type);
decays_.emplace(new_idx, decay_it->second);
}
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
}
namespace {
auto get_first_anyquark_emission(Event const & ev) {
// find born quarks (ignore extremal partons)
auto const firstquark = std::find_if(
std::next(ev.begin_partons()), std::prev(ev.end_partons(), 2),
[](Particle const & s){ return (is_anyquark(s)); }
);
// assert that it is a q-q_bar pair.
assert(std::distance(firstquark, ev.end_partons()) != 2);
assert(
( is_quark(*firstquark) && is_antiquark(*std::next(firstquark)) )
|| ( is_antiquark(*firstquark) && is_quark(*std::next(firstquark)) )
);
return firstquark;
}
//! returns index of most backward q-qbar jet
template<class Iterator>
int get_back_quark_jet(Event const & ev, Iterator firstquark){
// find jets at FO corresponding to the quarks
// technically this isn't necessary for LO
std::vector<int> const born_indices{ ev.particle_jet_indices() };
const auto firstquark_idx = std::distance(ev.begin_partons(), firstquark);
int const firstjet_idx = born_indices[firstquark_idx];
assert(firstjet_idx>0);
assert( born_indices[firstquark_idx+1] == firstjet_idx+1 );
return firstjet_idx;
}
//! returns index of most backward q-qbar jet
int getBackQuarkJet(Event const & ev){
const auto firstquark = get_first_anyquark_emission(ev);
return get_back_quark_jet(ev, firstquark);
}
template<class ConstIterator, class Iterator>
void label_extremal_qqx(
ConstIterator born_begin, ConstIterator born_end,
Iterator first_out
){
// find born quarks
const auto firstquark = std::find_if(
born_begin, born_end-1,
[](Particle const & s){ return (is_anyquark(s)); }
);
assert(firstquark != born_end-1);
const auto secondquark = std::find_if(
firstquark+1, born_end,
[](Particle const & s){ return (is_anyquark(s)); }
);
assert(secondquark != born_end);
assert( ( is_quark(*firstquark) && is_antiquark(*secondquark) )
|| ( is_antiquark(*firstquark) && is_quark(*secondquark) ));
assert(first_out->type == ParticleID::gluon);
assert((first_out+1)->type == ParticleID::gluon);
// copy type from born
first_out->type = firstquark->type;
(first_out+1)->type = secondquark->type;
}
}
void PhaseSpacePoint::label_qqx(Event const & event){
assert(std::is_sorted(begin(outgoing_), end(outgoing_), rapidity_less{}));
assert(filter_partons(outgoing_).size() == outgoing_.size());
if(qqxb_){
label_extremal_qqx( event.outgoing().cbegin(), event.outgoing().cend(),
outgoing_.begin()
);
return;
}
if(qqxf_){ // same as qqxb with reversed order
label_extremal_qqx( event.outgoing().crbegin(), event.outgoing().crend(),
outgoing_.rbegin()
);
return;
}
// central qqx
const auto firstquark = get_first_anyquark_emission(event);
// find jets at FO corresponding to the quarks
// technically this isn't necessary for LO
const auto firstjet_idx = get_back_quark_jet(event, firstquark);
// find corresponding jets after resummation
fastjet::ClusterSequence cs{to_PseudoJet(outgoing_), param_.jet_param.def};
auto const jets = fastjet::sorted_by_rapidity(
cs.inclusive_jets( param_.jet_param.min_pt ));
std::vector<int> const resum_indices{ cs.particle_jet_indices({jets}) };
// assert that jets didn't move
assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx )->rapidity(),
jets[ firstjet_idx ].rapidity(), 1e-2) );
assert(nearby_ep( ( event.jets().cbegin()+firstjet_idx+1 )->rapidity(),
jets[ firstjet_idx+1 ].rapidity(), 1e-2) );
// find last partons in first (central) jet
size_t idx_out = 0;
for(size_t i=resum_indices.size()-2; i>0; --i)
if(resum_indices[i] == firstjet_idx){
idx_out = i;
break;
}
assert(idx_out != 0);
// check that there is sufficient pt in jets from the quarks
const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction;
if (outgoing_[idx_out].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx )->pt()){
weight_=0.;
status_ = StatusCode::wrong_jets;
return;
}
if (outgoing_[idx_out+1].p.pt()<minpartonjetpt*( event.jets().cbegin()+firstjet_idx+1 )->pt()){
weight_=0.;
status_ = StatusCode::wrong_jets;
return;
}
// check that no additional emission between jets
// such configurations are possible if we have an gluon gets generated
// inside the rapidities of the qqx chain, but clusted to a
// differnet/outside jet. Changing this is non trivial
if(resum_indices[idx_out+1] != resum_indices[idx_out]+1){
weight_=0.;
status_ = StatusCode::gluon_in_qqx;
return;
}
outgoing_[idx_out].type = firstquark->type;
outgoing_[idx_out+1].type = std::next(firstquark)->type;
}
void PhaseSpacePoint::label_quarks(Event const & ev){
const auto WEmit = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Particle const & s){ return abs(s.type) == pid::Wp; }
);
if (WEmit != end(ev.outgoing())){
if(!qqxb_) {
const size_t backward_FKL_idx = unob_?1:0;
const auto backward_FKL = std::next(ev.begin_partons(), backward_FKL_idx);
outgoing_[backward_FKL_idx].type = backward_FKL->type;
}
if(!qqxf_) {
const size_t forward_FKL_idx = unof_?1:0;
const auto forward_FKL = std::prev(ev.end_partons(), 1+forward_FKL_idx);
outgoing_.rbegin()[unof_].type = forward_FKL->type;
}
} else {
most_backward_FKL(outgoing_).type = ev.incoming().front().type;
most_forward_FKL(outgoing_).type = ev.incoming().back().type;
}
if(qqxmid_||qqxb_||qqxf_){
label_qqx(ev);
}
}
PhaseSpacePoint::PhaseSpacePoint(
Event const & ev, PhaseSpacePointConfig conf, RNG & ran
):
unob_{ev.type() == event_type::unob},
unof_{ev.type() == event_type::unof},
qqxb_{ev.type() == event_type::qqxexb},
qqxf_{ev.type() == event_type::qqxexf},
qqxmid_{ev.type() == event_type::qqxmid},
param_{std::move(conf)},
status_{unspecified}
{
weight_ = 1;
const auto & Born_jets = ev.jets();
const int ng = sample_ng(Born_jets, ran);
weight_ /= std::tgamma(ng + 1);
const int ng_jets = sample_ng_jets(ng, Born_jets, ran);
std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
ng - ng_jets, CMINPT, param_.jet_param.min_pt, ran
);
int qqxbackjet(-1);
if(qqxmid_){
qqxbackjet = getBackQuarkJet(ev);
}
const auto qperp = std::accumulate(
begin(out_partons), end(out_partons),
fastjet::PseudoJet{}
);
const auto jets = reshuffle(Born_jets, qperp);
if(weight_ == 0.) {
status_ = failed_reshuffle;
return;
}
if(! pass_resummation_cuts(jets)){
status_ = failed_resummation_cuts;
weight_ = 0.;
return;
}
std::vector<fastjet::PseudoJet> jet_partons = split(
jets, ng_jets, qqxbackjet, ran
);
if(weight_ == 0.) {
status_ = StatusCode::failed_split;
return;
}
if(qqxmid_){
rescale_qqx_rapidities(
out_partons, jets,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity(),
qqxbackjet
);
}
else{
rescale_rapidities(
out_partons,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity()
);
}
if(! cluster_jets(out_partons).empty()){
weight_ = 0.;
status_ = StatusCode::empty_jets;
return;
}
std::sort(begin(out_partons), end(out_partons), rapidity_less{});
assert(
std::is_sorted(begin(jet_partons), end(jet_partons), rapidity_less{})
);
const auto first_jet_parton = out_partons.insert(
end(out_partons), begin(jet_partons), end(jet_partons)
);
std::inplace_merge(
begin(out_partons), first_jet_parton, end(out_partons), rapidity_less{}
);
if(! jets_ok(Born_jets, out_partons)){
weight_ = 0.;
status_ = StatusCode::wrong_jets;
return;
}
weight_ *= phase_space_normalisation(Born_jets.size(), out_partons.size());
outgoing_.reserve(out_partons.size() + 1); // one slot for possible A, W, Z, H
for(auto & p: out_partons){
outgoing_.emplace_back(Particle{pid::gluon, std::move(p), {}});
}
assert(!outgoing_.empty());
label_quarks(ev);
if(weight_ == 0.) {
//! @TODO optimise s.t. this is not possible
// status is handled internally
return;
}
copy_AWZH_boson_from(ev);
reconstruct_incoming(ev.incoming());
status_ = StatusCode::good;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
int count, double ptmin, double ptmax, RNG & ran
){
// heuristic parameters for pt sampling
const double ptpar = 1.3 + count/5.;
const double temp1 = atan((ptmax - ptmin)/ptpar);
std::vector<fastjet::PseudoJet> partons(count);
- for(size_t i = 0; i < (size_t) count; ++i){
+ for(size_t i = 0; i < static_cast<size_t>(count); ++i){
const double r1 = ran.flat();
const double pt = ptmin + ptpar*tan(r1*temp1);
const double temp2 = cos(r1*temp1);
const double phi = 2*M_PI*ran.flat();
weight_ *= 2.0*M_PI*pt*ptpar*temp1/(temp2*temp2);
// we don't know the allowed rapidity span yet,
// set a random value to be rescaled later on
const double y = ran.flat();
partons[i].reset_PtYPhiM(pt, y, phi);
// Set user index higher than any jet-parton index
// in order to assert that these are not inside jets
partons[i].set_user_index(i + 1 + ng_max);
assert(ptmin-1e-5 <= partons[i].pt() && partons[i].pt() <= ptmax+1e-5);
}
assert(std::all_of(partons.cbegin(), partons.cend(), is_nonjet_parton));
return sorted_by_rapidity(partons);
}
void PhaseSpacePoint::rescale_qqx_rapidities(
std::vector<fastjet::PseudoJet> & out_partons,
std::vector<fastjet::PseudoJet> const & jets,
const double ymin1, const double ymax2,
const int qqxbackjet
){
const double ymax1 = jets[qqxbackjet].rapidity();
const double ymin2 = jets[qqxbackjet+1].rapidity();
constexpr double ep = 1e-7;
const double tot_y = ymax1 - ymin1 + ymax2 - ymin2;
std::vector<std::reference_wrapper<fastjet::PseudoJet>> refpart(
out_partons.begin(), out_partons.end());
double ratio = (ymax1 - ymin1)/tot_y;
const auto gap{ std::find_if(refpart.begin(), refpart.end(),
[ratio](fastjet::PseudoJet p){
return (p.rapidity()>=ratio);} ) };
double ymin = ymin1;
double ymax = ymax1;
double dy = ymax - ymin - 2*ep;
double offset = 0.;
for(auto it_part=refpart.begin(); it_part<refpart.end(); ++it_part){
if(it_part == gap){
ymin = ymin2;
ymax = ymax2;
dy = ymax - ymin - 2*ep;
offset = ratio;
ratio = 1-ratio;
}
fastjet::PseudoJet & part = *it_part;
assert(offset <= part.rapidity() && part.rapidity() < ratio+offset);
const double y = ymin + ep + dy*((part.rapidity()-offset)/ratio);
part.reset_momentum_PtYPhiM(part.pt(), y, part.phi());
weight_ *= tot_y-4.*ep;
assert(ymin <= part.rapidity() && part.rapidity() <= ymax);
}
assert(is_sorted(begin(out_partons), end(out_partons), rapidity_less{}));
}
void PhaseSpacePoint::rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
){
constexpr double ep = 1e-7;
for(auto & parton: partons){
assert(0 <= parton.rapidity() && parton.rapidity() <= 1);
const double dy = ymax - ymin - 2*ep;
const double y = ymin + ep + dy*parton.rapidity();
parton.reset_momentum_PtYPhiM(parton.pt(), y, parton.phi());
weight_ *= dy;
assert(ymin <= parton.rapidity() && parton.rapidity() <= ymax);
}
}
namespace {
template<typename T, typename... Rest>
auto min(T const & a, T const & b, Rest&&... r) {
using std::min;
return min(a, min(b, std::forward<Rest>(r)...));
}
}
double PhaseSpacePoint::probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const{
assert(std::is_sorted(begin(Born_jets), end(Born_jets), rapidity_less{}));
assert(Born_jets.size() >= 2);
const double dy =
Born_jets.back().rapidity() - Born_jets.front().rapidity();
const double R = param_.jet_param.def.R();
const int njets = Born_jets.size();
const double p_J_y_large = (njets-1)*R*R/(2.*dy);
const double p_J_y0 = njets*R/M_PI;
return min(p_J_y_large, p_J_y0, 1.);
}
int PhaseSpacePoint::sample_ng_jets(
int ng, std::vector<fastjet::PseudoJet> const & Born_jets, RNG & ran
){
const double p_J = probability_in_jet(Born_jets);
std::binomial_distribution<> bin_dist(ng, p_J);
const int ng_J = bin_dist(ran);
weight_ *= std::pow(p_J, -ng_J)*std::pow(1 - p_J, ng_J - ng);
return ng_J;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
){
if(q == fastjet::PseudoJet{0, 0, 0, 0}) return Born_jets;
const auto jets = resummation_jet_momenta(Born_jets, q);
if(jets.empty()){
weight_ = 0;
return {};
}
// additional Jacobian to ensure Born integration over delta gives 1
weight_ *= resummation_jet_weight(Born_jets, q);
return jets;
}
std::vector<int> PhaseSpacePoint::distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets, RNG & ran
){
size_t first_valid_jet = 0;
size_t num_valid_jets = jets.size();
const double R_eff = 5./3.*param_.jet_param.def.R();
// if there is an unordered jet too far away from the FKL jets
// then extra gluon constituents of the unordered jet would
// violate the FKL rapidity ordering
if((unob_||qqxb_) && jets[0].delta_R(jets[1]) > R_eff){
++first_valid_jet;
--num_valid_jets;
}
else if((unof_||qqxf_) && jets[jets.size()-1].delta_R(jets[jets.size()-2]) > R_eff){
--num_valid_jets;
}
std::vector<int> np(jets.size(), 1);
for(int i = 0; i < ng_jets; ++i){
++np[first_valid_jet + ran.flat() * num_valid_jets];
}
weight_ *= std::pow(num_valid_jets, ng_jets);
return np;
}
#ifndef NDEBUG
namespace{
bool tagged_FKL_backward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return std::find_if(
begin(jet_partons), end(jet_partons),
[](fastjet::PseudoJet const & p){
- return p.user_index() == backward_FKL_idx;
+ return p.user_index() == backward_FKL_uid;
}
) != end(jet_partons);
}
bool tagged_FKL_forward(
std::vector<fastjet::PseudoJet> const & jet_partons
){
// the most forward FKL parton is most likely near the end of jet_partons;
// start search from there
return std::find_if(
jet_partons.rbegin(), jet_partons.rend(),
[](fastjet::PseudoJet const & p){
- return p.user_index() == forward_FKL_idx;
+ return p.user_index() == forward_FKL_uid;
}
) != jet_partons.rend();
}
bool tagged_FKL_extremal(
std::vector<fastjet::PseudoJet> const & jet_partons
){
return tagged_FKL_backward(jet_partons) && tagged_FKL_forward(jet_partons);
}
} // namespace anonymous
#endif
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
int ng_jets, size_t qqxbackjet, RNG & ran
){
return split(
jets, distribute_jet_partons(ng_jets, jets, ran), qqxbackjet, ran);
}
bool PhaseSpacePoint::pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const{
if(ext_parton.pt() < param_.min_extparton_pt) return false;
return (ext_parton - jet).pt()/jet.pt() < param_.max_ext_soft_pt_fraction;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np,
size_t qqxbackjet,
RNG & ran
){
assert(! jets.empty());
assert(jets.size() == np.size());
assert(pass_resummation_cuts(jets));
const size_t most_backward_FKL_idx = 0 + unob_ + qqxb_;
const size_t most_forward_FKL_idx = jets.size() - 1 - unof_ - qqxf_;
const auto & jet = param_.jet_param;
const JetSplitter jet_splitter{jet.def, jet.min_pt};
std::vector<fastjet::PseudoJet> jet_partons;
// randomly distribute jet gluons among jets
for(size_t i = 0; i < jets.size(); ++i){
auto split_res = jet_splitter.split(jets[i], np[i], ran);
weight_ *= split_res.weight;
if(weight_ == 0) return {};
assert(
std::all_of(
begin(split_res.constituents), end(split_res.constituents),
is_jet_parton
)
);
const auto first_new_parton = jet_partons.insert(
end(jet_partons),
begin(split_res.constituents), end(split_res.constituents)
);
// mark uno and extremal FKL emissions here so we can check
// their position once all emissions are generated
// also mark qqxmid partons, and apply appropriate pt cut.
auto extremal = end(jet_partons);
if (i == most_backward_FKL_idx){ //FKL backward emission
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
- extremal->set_user_index(backward_FKL_idx);
+ extremal->set_user_index(backward_FKL_uid);
}
else if(((unob_ || qqxb_) && i == 0)){
// unordered/qqxb
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
- extremal->set_user_index((unob_)?unob_idx:qqxb_idx);
+ extremal->set_user_index((unob_)?unob_uid:qqxb_uid);
}
else if (i == most_forward_FKL_idx){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
- extremal->set_user_index(forward_FKL_idx);
+ extremal->set_user_index(forward_FKL_uid);
}
else if(((unof_ || qqxf_) && i == jets.size() - 1)){
// unordered/qqxf
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
- extremal->set_user_index((unof_)?unof_idx:qqxf_idx);
+ extremal->set_user_index((unof_)?unof_uid:qqxf_uid);
}
else if((qqxmid_ && i == qqxbackjet)){
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
- extremal->set_user_index(qqxmid1_idx);
+ extremal->set_user_index(qqxmid1_uid);
}
else if((qqxmid_ && i == qqxbackjet+1)){
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
- extremal->set_user_index(qqxmid2_idx);
+ extremal->set_user_index(qqxmid2_uid);
}
if(
extremal != end(jet_partons)
&& !pass_extremal_cuts(*extremal, jets[i])
){
weight_ = 0;
return {};
}
}
assert(tagged_FKL_extremal(jet_partons));
std::sort(begin(jet_partons), end(jet_partons), rapidity_less{});
if(
!extremal_ok(jet_partons)
|| !split_preserved_jets(jets, jet_partons)
){
weight_ = 0.;
return {};
}
return jet_partons;
}
bool PhaseSpacePoint::extremal_ok(
std::vector<fastjet::PseudoJet> const & partons
) const{
assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
- if(unob_ && partons.front().user_index() != unob_idx) return false;
- if(unof_ && partons.back().user_index() != unof_idx) return false;
- if(qqxb_ && partons.front().user_index() != qqxb_idx) return false;
- if(qqxf_ && partons.back().user_index() != qqxf_idx) return false;
+ if(unob_ && partons.front().user_index() != unob_uid) return false;
+ if(unof_ && partons.back().user_index() != unof_uid) return false;
+ if(qqxb_ && partons.front().user_index() != qqxb_uid) return false;
+ if(qqxf_ && partons.back().user_index() != qqxf_uid) return false;
return
- most_backward_FKL(partons).user_index() == backward_FKL_idx
- && most_forward_FKL(partons).user_index() == forward_FKL_idx;
+ most_backward_FKL(partons).user_index() == backward_FKL_uid
+ && most_forward_FKL(partons).user_index() == forward_FKL_uid;
}
bool PhaseSpacePoint::split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const{
assert(std::is_sorted(begin(jets), end(jets), rapidity_less{}));
const auto split_jets = cluster_jets(jet_partons);
// this can happen if two overlapping jets
// are both split into more than one parton
if(split_jets.size() != jets.size()) return false;
for(size_t i = 0; i < split_jets.size(); ++i){
// this can happen if there are two overlapping jets
// and a parton is assigned to the "wrong" jet
if(!nearby_ep(jets[i].rapidity(), split_jets[i].rapidity(), 1e-2)){
return false;
}
}
return true;
}
template<class Particle>
Particle const & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> const & partons
) const{
return partons[0 + unob_ + qqxb_];
}
template<class Particle>
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t idx = partons.size() - 1 - unof_ - qqxf_;
assert(idx < partons.size());
return partons[idx];
}
template<class Particle>
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
return partons[0 + unob_ + qqxb_];
}
template<class Particle>
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t idx = partons.size() - 1 - unof_ - qqxf_;
assert(idx < partons.size());
return partons[idx];
}
bool PhaseSpacePoint::contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
) const {
auto const & constituents = jet.constituents();
const int idx = parton.user_index();
const bool injet = std::find_if(
begin(constituents), end(constituents),
[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
) != end(constituents);
const double minpartonjetpt = 1. - param_.max_ext_soft_pt_fraction;
return ((parton.pt()>minpartonjetpt*jet.pt())&&injet);
}
bool PhaseSpacePoint::jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, param_.jet_param.def);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt));
if(jets.size() != Born_jets.size()) return false;
int in_jet = 0;
for(size_t i = 0; i < jets.size(); ++i){
assert(jets[i].has_constituents());
for(auto && parton: jets[i].constituents()){
if(is_nonjet_parton(parton)) return false;
}
in_jet += jets[i].constituents().size();
}
const int expect_in_jet = std::count_if(
partons.cbegin(), partons.cend(), is_jet_parton
);
if(in_jet != expect_in_jet) return false;
// note that PseudoJet::contains does not work here
if(! (
contains_idx(most_backward_FKL(jets), most_backward_FKL(partons))
&& contains_idx(most_forward_FKL(jets), most_forward_FKL(partons))
)) return false;
if(unob_ && !contains_idx(jets.front(), partons.front())) return false;
if(qqxb_ && !contains_idx(jets.front(), partons.front())) return false;
if(unof_ && !contains_idx(jets.back(), partons.back())) return false;
if(qqxf_ && !contains_idx(jets.back(), partons.back())) return false;
for(size_t i = 0; i < jets.size(); ++i){
assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
}
return true;
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Particle, 2> const & Born_incoming
){
std::tie(incoming_[0].p, incoming_[1].p) = incoming_momenta(outgoing_);
for(size_t i = 0; i < incoming_.size(); ++i){
incoming_[i].type = Born_incoming[i].type;
}
assert(momentum_conserved());
}
double PhaseSpacePoint::phase_space_normalisation(
int num_Born_jets, int num_out_partons
) const{
return pow(16*pow(M_PI,3), num_Born_jets - num_out_partons);
}
bool PhaseSpacePoint::momentum_conserved() const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
const double norm = diff.E();
for(auto const & out: outgoing()) diff -= out.p;
return nearby(diff, fastjet::PseudoJet{}, norm);
}
} //namespace HEJ
diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
index f5c9ded..165c5eb 100644
--- a/src/bin/HEJ.cc
+++ b/src/bin/HEJ.cc
@@ -1,382 +1,382 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <array>
#include <chrono>
#include <iostream>
#include <limits>
#include <memory>
#include <numeric>
#include "yaml-cpp/yaml.h"
#include "fastjet/ClusterSequence.hh"
#include "HEJ/CombinedEventWriter.hh"
#include "HEJ/Config.hh"
#include "HEJ/CrossSectionAccumulator.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
#include "HEJ/BufferedEventReader.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/get_analysis.hh"
#include "HEJ/make_RNG.hh"
#include "HEJ/optional.hh"
#include "HEJ/ProgressBar.hh"
#include "HEJ/stream.hh"
#include "HEJ/Unweighter.hh"
#include "HEJ/Version.hh"
#include "HEJ/YAMLreader.hh"
HEJ::Config load_config(char const * filename){
try{
return HEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::vector<std::unique_ptr<HEJ::Analysis>> get_analyses(
std::vector<YAML::Node> const & parameters, LHEF::HEPRUP const & heprup
){
try{
return HEJ::get_analyses(parameters, heprup);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
// unique_ptr is a workaround:
// HEJ::optional is a better fit, but gives spurious errors with g++ 7.3.0
std::unique_ptr<HEJ::ProgressBar<double>> make_progress_bar(
std::vector<double> const & xs
) {
if(xs.empty()) return {};
const double Born_xs = std::accumulate(begin(xs), end(xs), 0.);
return std::make_unique<HEJ::ProgressBar<double>>(std::cout, Born_xs);
}
std::string time_to_string(const time_t time){
char s[30];
struct tm * p = localtime(&time);
strftime(s, 30, "%a %b %d %Y %H:%M:%S", p);
return s;
}
HEJ::Event to_event(
LHEF::HEPEUP const & hepeup,
HEJ::JetParameters const & fixed_order_jets
) {
HEJ::Event::EventData event_data{hepeup};
event_data.reconstruct_intermediate();
return HEJ::Event{
std::move(event_data).cluster(
fixed_order_jets.def, fixed_order_jets.min_pt
)
};
}
void unweight(
HEJ::Unweighter & unweighter,
HEJ::WeightType weight_type,
std::vector<HEJ::Event> & events,
HEJ::RNG & ran
) {
if(weight_type == HEJ::WeightType::unweighted_resum){
unweighter.set_cut_to_maxwt(events);
}
events.erase(
unweighter.unweight(begin(events), end(events), ran),
end(events)
);
}
// peek up to nevents events from reader
std::vector<LHEF::HEPEUP> peek_events(
HEJ::BufferedEventReader & reader,
const int nevents
) {
std::vector<LHEF::HEPEUP> events;
while(
static_cast<int>(events.size()) < nevents
&& reader.read_event()
) {
events.emplace_back(reader.hepeup());
}
// put everything back into the reader
for(auto it = rbegin(events); it != rend(events); ++it) {
reader.emplace(*it);
}
return events;
}
void append_resummed_events(
std::vector<HEJ::Event> & resummation_events,
HEJ::EventReweighter & reweighter,
LHEF::HEPEUP const & hepeup,
- const int trials,
+ const size_t trials,
HEJ::JetParameters const & fixed_order_jets
) {
const HEJ::Event FO_event = to_event(hepeup, fixed_order_jets);
if(reweighter.treatment(FO_event.type()) != HEJ::EventTreatment::reweight) {
return;
}
const auto resummed = reweighter.reweight(FO_event, trials);
resummation_events.insert(
end(resummation_events),
begin(resummed), end(resummed)
);
}
void train(
HEJ::Unweighter & unweighter,
HEJ::BufferedEventReader & reader,
HEJ::EventReweighter & reweighter,
- const int total_trials,
+ const size_t total_trials,
const double max_dev,
double reweight_factor,
HEJ::JetParameters const & fixed_order_jets
) {
std::cout << "Reading up to " << total_trials << " training events...\n";
auto FO_events = peek_events(reader, total_trials);
if(FO_events.empty()) {
throw std::runtime_error{
"No events generated to calibrate the unweighting weight!"
"Please increase the number \"trials\" or deactivate the unweighting."
};
}
- const int trials = total_trials/FO_events.size();
+ const size_t trials = total_trials/FO_events.size();
// adjust reweight factor so that the overall normalisation
// is the same as in the full run
reweight_factor *= trials;
for(auto & hepeup: FO_events) {
hepeup.XWGTUP *= reweight_factor;
}
std::cout << "Training unweighter with "
<< trials << '*' << FO_events.size() << " events\n";
auto progress = HEJ::ProgressBar<int>{
std::cout, static_cast<int>(FO_events.size())
};
std::vector<HEJ::Event> resummation_events;
for(auto const & hepeup: FO_events) {
append_resummed_events(
resummation_events,
reweighter, hepeup, trials, fixed_order_jets
);
++progress;
}
unweighter.set_cut_to_peakwt(resummation_events, max_dev);
std::cout << "\nUnweighting events with weight up to "
<< unweighter.get_cut() << '\n';
}
int main(int argn, char** argv) {
using clock = std::chrono::system_clock;
if (argn != 3) {
std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
{
std::cout << "Starting " << HEJ::Version::package_name_full()
<< ", revision " << HEJ::Version::revision() << " ("
<< time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
}
fastjet::ClusterSequence::print_banner();
// read configuration
const HEJ::Config config = load_config(argv[1]);
auto reader = HEJ::make_reader(argv[2]);
assert(reader);
auto heprup{ reader->heprup() };
heprup.generators.emplace_back(LHEF::XMLTag{});
heprup.generators.back().name = HEJ::Version::package_name();
heprup.generators.back().version = HEJ::Version::String();
auto analyses = get_analyses( config.analyses_parameters, heprup );
assert(analyses.empty() || analyses.front() != nullptr);
HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};
double global_reweight = 1.;
const auto & max_events = config.max_events;
// if we need the event number:
if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1 || max_events){
// try to read from LHE head
auto input_events{reader->number_events()};
if(!input_events) {
// else count manually
auto t_reader = HEJ::make_reader(argv[2]);
input_events = 0;
while(t_reader->read_event()) ++(*input_events);
}
if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1){
// IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs
std::cout << "Found IDWTUP " << heprup.IDWTUP << ": "
<< "assuming \"cross section = average weight\".\n"
<< "converting to \"cross section = sum of weights\" ";
global_reweight /= *input_events;
}
if(max_events && (*input_events > *max_events)){
// maximal number of events given
global_reweight *= *input_events/static_cast<double>(*max_events);
std::cout << "Processing " << *max_events
<< " out of " << *input_events << " events\n";
}
}
HEJ::ScaleGenerator scale_gen{
config.scales.base,
config.scales.factors,
config.scales.max_ratio
};
std::shared_ptr<HEJ::RNG> ran{
HEJ::make_RNG(config.rng.name, config.rng.seed)};
assert(ran != nullptr);
HEJ::EventReweighter hej{
reader->heprup(),
std::move(scale_gen),
to_EventReweighterConfig(config),
ran
};
HEJ::optional<HEJ::Unweighter> unweighter{};
if(config.weight_type != HEJ::WeightType::weighted) {
unweighter = HEJ::Unweighter{};
}
if(config.weight_type == HEJ::WeightType::partially_unweighted) {
HEJ::BufferedEventReader buffered_reader{std::move(reader)};
assert(config.unweight_config);
train(
*unweighter,
buffered_reader,
hej,
config.unweight_config->trials,
config.unweight_config->max_dev,
global_reweight/config.trials,
config.fixed_order_jets
);
reader = std::make_unique<HEJ::BufferedEventReader>(
std::move(buffered_reader)
);
}
// status infos & eye candy
size_t nevent = 0;
std::array<int, HEJ::event_type::last_type + 1>
nevent_type{0}, nfailed_type{0};
auto progress = make_progress_bar(reader->heprup().XSECUP);
HEJ::CrossSectionAccumulator xs;
std::map<HEJ::StatusCode, int> status_counter;
size_t total_trials = 0;
size_t total_resum = 0;
// Loop over the events in the input file
while(reader->read_event() && (!max_events || nevent < *max_events) ){
++nevent;
// reweight events so that the total cross section is conserved
auto hepeup = reader->hepeup();
hepeup.XWGTUP *= global_reweight;
const auto FO_event = to_event(hepeup, config.fixed_order_jets);
if(FO_event.central().weight == 0) {
static const bool warned_once = [argv,nevent](){
std::cerr
<< "WARNING: event number " << nevent
<< " in " << argv[2] << " has zero weight. "
"Ignoring this and all further events with vanishing weight.\n";
return true;
}();
(void) warned_once; // shut up compiler warnings
continue;
}
auto resummed_events{ hej.reweight(FO_event, config.trials) };
// some bookkeeping
for(auto const & s: hej.status())
++status_counter[s];
total_trials+=hej.status().size();
++nevent_type[FO_event.type()];
if(resummed_events.empty()) ++nfailed_type[FO_event.type()];
if(unweighter) {
unweight(*unweighter, config.weight_type, resummed_events, *ran);
}
// analysis
for(auto & ev: resummed_events){
//TODO: move pass_cuts to after phase space point generation
bool passed = analyses.empty();
for(auto const & analysis: analyses){
if(analysis->pass_cuts(ev, FO_event)){
passed = true;
analysis->fill(ev, FO_event);
};
}
if(passed){
writer.write(ev);
} else {
ev.parameters()*=0; // do not use discarded events afterwards
}
}
xs.fill_correlated(resummed_events);
total_resum += resummed_events.size();
if(progress) progress->increment(FO_event.central().weight);
} // main event loop
std::cout << '\n';
for(auto const & analysis: analyses){
analysis->finalise();
}
using namespace HEJ::event_type;
std::cout<< "Events processed: " << nevent << " (" << total_resum << " resummed)"<< '\n';
std::cout << '\t' << name(EventType::first_type) << ": "
<< nevent_type[EventType::first_type]
<< ", failed to reconstruct " << nfailed_type[EventType::first_type]
<< '\n';
for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){
std::cout << '\t' << name(static_cast<EventType>(i)) << ": "
<< nevent_type[i]
<< ", failed to reconstruct " << nfailed_type[i]
<< '\n';
}
std::cout << '\n' << xs << '\n';
std::cout << "Generation statistic: "
<< status_counter[HEJ::StatusCode::good] << "/" << total_trials
<< " trials successful.\n";
for(auto && entry: status_counter){
const double fraction = static_cast<double>(entry.second)/total_trials;
const int percent = std::round(100*fraction);
std::cout << std::left << std::setw(17) << (to_string(entry.first) + ":")
<< " [";
for(int i = 0; i < percent/2; ++i) std::cout << '#';
for(int i = percent/2; i < 50; ++i) std::cout << ' ';
std::cout << "] " <<std::setw(2)<<std::right<< percent << "%\n";
}
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nFinished " << HEJ::Version::package_name() << " at "
<< time_to_string(clock::to_time_t(clock::now()))
<< "\n=> Runtime: " << run_time.count() << " sec ("
<< nevent/run_time.count() << " Events/sec).\n";
return EXIT_SUCCESS;
}
diff --git a/t/check_res.cc b/t/check_res.cc
index 3e74e19..7f1ed62 100644
--- a/t/check_res.cc
+++ b/t/check_res.cc
@@ -1,168 +1,168 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <cmath>
#include <iostream>
#include "LHEF/LHEF.h"
#include "HEJ/CrossSectionAccumulator.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReweighter.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/stream.hh"
#include "hej_test.hh"
namespace{
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
const fastjet::JetDefinition Born_jet_def{jet_def};
constexpr double Born_jetptmin = 30;
constexpr double max_ext_soft_pt_fraction = 0.1;
constexpr double jetptmin = 35;
constexpr bool log_corr = false;
const HEJ::ParticleProperties Wprop{80.385, 2.085};
const HEJ::ParticleProperties Zprop{91.187, 2.495};
const HEJ::ParticleProperties Hprop{125, 0.004165};
constexpr double vev = 246.2196508;
using EventTreatment = HEJ::EventTreatment;
using namespace HEJ::event_type;
HEJ::EventTreatMap treat{
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{non_resummable, EventTreatment::discard},
{unof, EventTreatment::discard},
{unob, EventTreatment::discard},
{qqxexb, EventTreatment::discard},
{qqxexf, EventTreatment::discard},
{qqxmid, EventTreatment::discard},
{FKL, EventTreatment::reweight}
};
/// true if colour is allowed for particle
bool correct_colour(HEJ::Particle const & part){
if(HEJ::is_AWZH_boson(part) && !part.colour) return true;
if(!part.colour) return false;
int const colour = part.colour->first;
int const anti_colour = part.colour->second;
if(part.type == HEJ::ParticleID::gluon)
return colour != anti_colour && colour > 0 && anti_colour > 0;
if(HEJ::is_quark(part))
return anti_colour == 0 && colour > 0;
return colour == 0 && anti_colour > 0;
}
bool correct_colour(HEJ::Event const & ev){
if(!HEJ::event_type::is_resummable(ev.type()))
return true;
for(auto const & part: ev.incoming()){
if(!correct_colour(part))
return false;
}
for(auto const & part: ev.outgoing()){
if(!correct_colour(part))
return false;
}
return true;
}
}
int main(int argn, char** argv) {
if(argn == 5 && std::string(argv[4]) == "unof"){
--argn;
treat[unof] = EventTreatment::reweight;
treat[unob] = EventTreatment::discard;
treat[FKL] = EventTreatment::discard;
}
if(argn == 5 && std::string(argv[4]) == "unob"){
--argn;
treat[unof] = EventTreatment::discard;
treat[unob] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "splitf"){
--argn;
treat[qqxexb] = EventTreatment::discard;
treat[qqxexf] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "splitb"){
--argn;
treat[qqxexb] = EventTreatment::reweight;
treat[qqxexf] = EventTreatment::discard;
treat[FKL] = EventTreatment::discard;
}
else if(argn == 5 && std::string(argv[4]) == "qqxmid"){
--argn;
treat[qqxmid] = EventTreatment::reweight;
treat[FKL] = EventTreatment::discard;
}
if(argn != 4){
std::cerr << "Usage: check_res eventfile xsection tolerance [uno]";
return EXIT_FAILURE;
}
const double xsec_ref = std::stod(argv[2]);
const double tolerance = std::stod(argv[3]);
HEJ::istream in{argv[1]};
LHEF::Reader reader{in};
HEJ::PhaseSpacePointConfig psp_conf;
psp_conf.jet_param = HEJ::JetParameters{jet_def, jetptmin};
psp_conf.max_ext_soft_pt_fraction = max_ext_soft_pt_fraction;
HEJ::MatrixElementConfig ME_conf;
ME_conf.log_correction = log_corr;
ME_conf.Higgs_coupling = HEJ::HiggsCouplingSettings{};
ME_conf.ew_parameters.set_vevWZH(vev, Wprop, Zprop, Hprop);
HEJ::EventReweighterConfig conf;
conf.psp_config = std::move(psp_conf);
conf.ME_config = std::move(ME_conf);
conf.treat = treat;
reader.readEvent();
const bool has_Higgs = std::find(
begin(reader.hepeup.IDUP),
end(reader.hepeup.IDUP),
25
) != end(reader.hepeup.IDUP);
const double mu = has_Higgs?125.:91.188;
HEJ::ScaleGenerator scale_gen{
{{std::to_string(mu), HEJ::FixedScale{mu}}}, {}, 1.
};
std::shared_ptr<HEJ::RNG> ran{std::make_shared<HEJ::Mixmax>()};
HEJ::EventReweighter hej{reader.heprup, std::move(scale_gen), conf, ran};
HEJ::CrossSectionAccumulator xs;
do{
auto ev_data = HEJ::Event::EventData{reader.hepeup};
shuffle_particles(ev_data);
ev_data.reconstruct_intermediate();
HEJ::Event ev{
ev_data.cluster(
Born_jet_def, Born_jetptmin
)
};
auto resummed_events = hej.reweight(ev, 100);
- for(auto const & ev: resummed_events) {
- ASSERT(correct_colour(ev));
- ASSERT(std::isfinite(ev.central().weight));
+ for(auto const & res_ev: resummed_events) {
+ ASSERT(correct_colour(res_ev));
+ ASSERT(std::isfinite(res_ev.central().weight));
// we fill the xs uncorrelated since we only want to test the uncertainty
// of the resummation
- xs.fill(ev);
+ xs.fill(res_ev);
}
} while(reader.readEvent());
const double xsec = xs.total().value;
const double xsec_err = std::sqrt(xs.total().error);
const double significance =
std::abs(xsec - xsec_ref) / std::sqrt( xsec_err*xsec_err + tolerance*tolerance );
std::cout << xsec_ref << " +/- " << tolerance << " ~ "
<< xsec << " +- " << xsec_err << " => " << significance << " sigma\n";
if(significance > 3.){
std::cerr << "Cross section is off by over 3 sigma!\n";
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
diff --git a/t/test_colours.cc b/t/test_colours.cc
index 78540c9..5435f44 100644
--- a/t/test_colours.cc
+++ b/t/test_colours.cc
@@ -1,352 +1,352 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <stdexcept>
#include <utility>
#include "HEJ/Constants.hh"
#include "HEJ/Event.hh"
#include "HEJ/RNG.hh"
#include "hej_test.hh"
/// biased RNG to connect always to colour
class dum_rnd: public HEJ::DefaultRNG {
public:
dum_rnd() = default;
double flat() override {
return 0.;
- };
+ }
};
HEJ::Event::EventData decay_boson( HEJ::Event::EventData ev ){
for( size_t i=0; i<ev.outgoing.size(); ++i ){
if( std::abs(ev.outgoing[i].type) == HEJ::ParticleID::Wp){
ev.decays[i] = decay_W(ev.outgoing[i]);
}
}
return ev;
}
void dump_event(HEJ::Event const & ev){
for(auto const & in: ev.incoming()){
std::cerr << "in type=" << in.type
<< ", colour={" << (*in.colour).first
<< ", " << (*in.colour).second << "}\n";
}
for(auto const & out: ev.outgoing()){
std::cerr << "out type=" << out.type << ", colour={";
if(out.colour)
std::cerr << (*out.colour).first << ", " << (*out.colour).second;
else
std::cerr << "non, non";
std::cerr << "}\n";
}
}
/// true if colour is allowed for particle
bool correct_colour(HEJ::Particle const & part){
if(!HEJ::is_parton(part) && !part.colour) return true;
if(!part.colour) return false;
int const colour = part.colour->first;
int const anti_colour = part.colour->second;
if(part.type == HEJ::ParticleID::gluon)
return colour != anti_colour
&& colour >= HEJ::COLOUR_OFFSET
&& anti_colour >= HEJ::COLOUR_OFFSET;
if(HEJ::is_quark(part))
return anti_colour == 0 && colour >= HEJ::COLOUR_OFFSET;
return colour == 0 && anti_colour >= HEJ::COLOUR_OFFSET;
}
bool correct_colour(HEJ::Event const & ev){
if(!ev.is_leading_colour())
return false;
// some of these additional checks are also in ev.is_leading_colour()
for(auto const & part: ev.incoming()){
if(!correct_colour(part))
return false;
}
for(auto const & part: ev.outgoing()){
if(!correct_colour(part))
return false;
}
return true;
}
bool match_expected(
HEJ::Event const & ev,
std::vector<HEJ::Colour> const & expected
){
ASSERT(ev.outgoing().size()+2==expected.size());
for(size_t i=0; i<ev.incoming().size(); ++i){
ASSERT(ev.incoming()[i].colour);
if( *ev.incoming()[i].colour != expected[i])
return false;
}
for(size_t i=2; i<ev.outgoing().size()+2; ++i){
if( ev.outgoing()[i-2].colour ){
if( *ev.outgoing()[i-2].colour != expected[i] )
return false;
} else if( expected[i].first != 0 || expected[i].second != 0)
return false;
}
return true;
}
void check_event(
HEJ::Event::EventData unc_ev,
std::vector<HEJ::Colour> const & expected_colours
){
unc_ev = decay_boson(std::move(unc_ev));
shuffle_particles(unc_ev); // make sure incoming order doesn't matter
HEJ::Event ev{unc_ev.cluster(
fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 30.)
};
ASSERT(HEJ::event_type::is_resummable(ev.type()));
dum_rnd rng;
ASSERT(!ev.is_leading_colour());
ASSERT(ev.generate_colours(rng));
if(!correct_colour(ev)){
std::cerr << "Found illegal colours for event\n";
dump_event(ev);
throw std::invalid_argument("Illegal colour set");
}
if(!match_expected(ev, expected_colours)){
std::cerr << "Colours didn't match expectation. Found\n";
dump_event(ev);
std::cerr << "but expected\n";
for(auto const & col: expected_colours){
std::cerr << "colour={" << col.first << ", " << col.second << "}\n";
}
throw std::logic_error("Colours did not match expectation");
}
}
HEJ::Event::EventData reset_colour(
HEJ::Event::EventData ev, std::vector<HEJ::Colour> const & goal
){
for(size_t i=0; i<2; ++i){
ev.incoming[i].colour = goal[i];
}
for(size_t i=0; i<ev.outgoing.size(); ++i){
auto const & col_goal{ goal[i+2] };
if(col_goal.first == 0 && col_goal.second == 0)
ev.outgoing[i].colour = HEJ::optional<HEJ::Colour>{};
else
ev.outgoing[i].colour = col_goal;
}
return ev;
}
int main() {
HEJ::Event::EventData ev;
std::vector<HEJ::Colour> expected_colours(7);
/// pure gluon (they all have a mass of 4 GeV to allow decays)
ev.incoming[0] = { HEJ::ParticleID::gluon, { 0, 0, -205, 205}, {}};
ev.incoming[1] = { HEJ::ParticleID::gluon, { 0, 0, 279, 279}, {}};
ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-15, -82, -82, 117}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, { 68, 93, 20, 117}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::higgs, {-30, -65, 22, 75}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-12, 92, 76, 120}, {}});
ev.outgoing.push_back({ HEJ::ParticleID::gluon, {-11, -38, 38, 55}, {}});
expected_colours[0] = {502, 501};
expected_colours[1] = {509, 502};
expected_colours[2] = {503, 501};
expected_colours[3] = {505, 503};
expected_colours[4] = { 0, 0};
expected_colours[5] = {507, 505};
expected_colours[6] = {509, 507};
// default colours is always forbidden!
// default: swap last two (anti-)colour -> crossing
ev=reset_colour(ev, expected_colours);
std::swap(ev.outgoing[4].colour, ev.outgoing[3].colour);
check_event(ev, expected_colours);
/// last g to Qx (=> gQx -> g ... Qx)
ev.incoming[1].type = HEJ::ParticleID::d_bar;
ev.outgoing[4].type = HEJ::ParticleID::d_bar;
// => only end changes
expected_colours[1].first = 0;
expected_colours[6].first = 0;
// default: swap last two anti-colours -> last gluon colour singlet
ev=reset_colour(ev, expected_colours);
std::swap(ev.outgoing[4].colour->second, ev.outgoing[3].colour->second);
check_event(ev, expected_colours);
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno forward (=> gQx -> g ... Qx g)
std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
// => uno quarks eats colour and gluon connects to anti-colour
new_expected[5] = {0, expected_colours[3].first};
new_expected[6] = {expected_colours[0].first, expected_colours[0].first+2};
new_expected[1].second += 2; // one more anti-colour in line
// default: swap last two anti-colours -> crossing
new_ev=reset_colour(new_ev, new_expected);
std::swap(new_ev.outgoing[4].colour->second, new_ev.outgoing[3].colour->second);
check_event(new_ev, new_expected);
}
/// swap Qx <-> Q (=> gQ -> g ... Q)
ev.incoming[1].type = HEJ::ParticleID::d;
ev.outgoing[4].type = HEJ::ParticleID::d;
// => swap: colour<->anti && initial<->final
std::swap(expected_colours[1], expected_colours[6]);
std::swap(expected_colours[1].first, expected_colours[1].second);
std::swap(expected_colours[6].first, expected_colours[6].second);
// default: swap incoming <-> outgoing
ev=reset_colour(ev, expected_colours);
std::swap(ev.incoming[0].colour, ev.outgoing[0].colour);
check_event(ev, expected_colours);
/// first g to qx (=> qxQ -> qx ... Q)
ev.incoming[0].type = HEJ::ParticleID::u_bar;
ev.outgoing[0].type = HEJ::ParticleID::u_bar;
expected_colours[0] = { 0, 501};
// => shift anti-colour index one up
expected_colours[1].first -= 2;
expected_colours[5] = expected_colours[3];
expected_colours[3] = expected_colours[2];
expected_colours[2] = { 0, 502};
// default: closed qx->qx g
ev=reset_colour(ev, expected_colours);
ev.outgoing[1].colour->first = ev.outgoing[0].colour->second;
ev.outgoing[1].colour->second = ev.incoming[0].colour->second;
ev.outgoing[4].colour->first = ev.outgoing[3].colour->second;
check_event(ev, expected_colours);
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno backward (=> qxQ -> g qx ... Q)
std::swap(new_ev.outgoing[0].type, new_ev.outgoing[1].type);
// => uno gluon connects to quark colour
new_expected[3] = expected_colours[2];
new_expected[2] = {expected_colours[0].second+2, expected_colours[0].second};
// default: Colourful Higgs
new_ev=reset_colour(new_ev, new_expected);
new_ev.outgoing[2].colour = std::make_pair(1,1);
check_event(new_ev, new_expected);
/// swap qx <-> q (=> qQ -> g q ... Q)
new_ev.incoming[0].type = HEJ::ParticleID::u;
new_ev.outgoing[1].type = HEJ::ParticleID::u;
// => swap: colour<->anti && inital<->final
std::swap(new_expected[0], new_expected[3]);
std::swap(new_expected[0].first, new_expected[0].second);
std::swap(new_expected[3].first, new_expected[3].second);
// => & connect first gluon with remaining anti-colour
new_expected[2] = {new_expected[0].first, new_expected[0].first+2};
// shift colour line one down
new_expected[1].first-=2;
new_expected[5].first-=2;
new_expected[5].second-=2;
// shift anti-colour line one up
new_expected[6].first+=2;
// default: swap 2 quarks -> disconnected
new_ev=reset_colour(new_ev, new_expected);
std::swap(new_ev.outgoing[1].colour, new_ev.outgoing[4].colour);
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// uno forward (=> qxQ -> qx ... Q g)
std::swap(new_ev.outgoing[3].type, new_ev.outgoing[4].type);
// => uno gluon connects to remaining colour
new_expected[5] = expected_colours[6];
new_expected[6] = {expected_colours[3].first+2, expected_colours[3].first};
// default: no colour on last gluon
new_ev=reset_colour(new_ev, new_expected);
new_ev.incoming[1].colour->first = new_ev.outgoing[4].colour->second;
new_ev.outgoing[4].colour = {};
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// qqx backward (=> gQ -> qx q ... Q) with Wp
// => swap: incoming q <-> outgoing gluon
std::swap(new_ev.incoming[0].type, new_ev.outgoing[1].type);
new_ev.outgoing[1].type=static_cast<HEJ::ParticleID>(
-(new_ev.outgoing[1].type+1) );
new_ev.outgoing[2].type = HEJ::ParticleID::Wp;
// incoming q -> outgoing q (colour<->anti)
std::swap(new_expected[0], new_expected[3]);
std::swap(new_expected[3].first, new_expected[3].second);
new_expected[3].first+=2;
new_expected[0].first-=1; // skip one index
// couple first in to first out
new_expected[2].second=new_expected[0].second;
// default: swap qqx <-> first g
new_ev=reset_colour(new_ev, new_expected);
std::swap(new_ev.outgoing[0].colour->second, new_ev.outgoing[3].colour->second);
std::swap(new_ev.outgoing[1].colour->first, new_ev.outgoing[3].colour->first);
check_event(new_ev, new_expected);
}
{
// don't overwrite
auto new_expected = expected_colours;
auto new_ev = ev;
/// qqx forward (=> qx g -> qx ... Qx Q) with Wp
// => swap: incoming Q <-> outgoing gluon
std::swap(new_ev.incoming[1].type, new_ev.outgoing[3].type);
new_ev.outgoing[3].type=static_cast<HEJ::ParticleID>(
-(new_ev.outgoing[3].type+1));
new_ev.outgoing[2].type = HEJ::ParticleID::Wp;
// incoming q -> outgoing q (colour<->anti)
std::swap(new_expected[1], new_expected[5]);
std::swap(new_expected[5].first, new_expected[5].second);
new_expected[5].second-=2;
new_expected[1].second-=1; // skip one index
// couple last in to last out
new_expected[6].first=new_expected[1].first;
// default: uncoloured quark
new_ev=reset_colour(new_ev, new_expected);
new_ev.outgoing[0].colour = {};
check_event(new_ev, new_expected);
// move Higgs to position 1 (=> qx g -> qx h g Qx Q)
std::swap(new_ev.outgoing[1].type, new_ev.outgoing[2].type);
std::swap(new_expected[3], new_expected[4]); // trivial
// default: incoming qx wrong colour
new_ev=reset_colour(new_ev, new_expected);
new_ev.incoming[0].colour->first = 1;
check_event(new_ev, new_expected);
// central qqx (=> qx g -> qx h Q Qx g)
// => swap: Q <-> g
std::swap(new_ev.outgoing[2].type, new_ev.outgoing[4].type);
std::swap(new_expected[4], new_expected[6]);
// gluon was connected on left side, i.e. doesn't matter for QQx
// => couple Q to out qx
new_expected[4].first = new_expected[2].second;
// Qx next in line
new_expected[5].second = new_expected[4].first+2;
// incoming g shifted by one position in line
new_expected[1].first-=2;
new_expected[1].second+=2;
// default: wrong colour in last incoming
new_ev=reset_colour(new_ev, new_expected);
std::swap(new_ev.incoming[1].colour->first,
new_ev.incoming[1].colour->second);
check_event(new_ev, new_expected);
}
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 3:00 PM (5 h, 43 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486787
Default Alt Text
(237 KB)

Event Timeline