diff --git a/include/HEJ/Event.hh b/include/HEJ/Event.hh
index 2d11440..2e89f8c 100644
--- a/include/HEJ/Event.hh
+++ b/include/HEJ/Event.hh
@@ -1,342 +1,353 @@
/** \file
* \brief Declares the Event class and helpers
*
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <memory>
#include <string>
#include <unordered_map>
#include <vector>
#include "HEJ/event_types.hh"
#include "HEJ/Particle.hh"
#include "fastjet/ClusterSequence.hh"
namespace LHEF{
class HEPEUP;
class HEPRUP;
}
namespace fastjet{
class JetDefinition;
}
namespace HEJ{
struct ParameterDescription;
//! 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;
};
//! Description of event parameters
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}
{};
};
struct UnclusteredEvent;
- /** An event with clustered jets
+ /** @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.
+ *
+ * \note Use EventData to build this class.
+ * There is no other constructor available.
*/
class Event{
public:
class EventData;
- //! Default Event Constructor
+ //! No default Constructor
Event() = delete;
//! Event Constructor adding jet clustering to an unclustered event
//! @deprecated UnclusteredEvent will be replaced by EventData in HEJ 2.3.0
[[deprecated("UnclusteredEvent will be replaced by EventData")]]
Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
);
//! The jets formed by the outgoing partons
std::vector<fastjet::PseudoJet> jets() const;
//! The corresponding event before jet clustering
// [[deprecated]]
// UnclusteredEvent unclustered() const {
// return UnclusteredEvent(ev_);
// TODO what to do with this?
// }
//! Incoming particles
std::array<Particle, 2> const & incoming() const{
return incoming_;
}
//! Outgoing particles
std::vector<Particle> const & outgoing() const{
return outgoing_;
}
//! 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_;
}
- //! Central parameter choice
+ //! Central parameter choice (const version)
EventParameters const & central() const{
return central_;
}
//! Central parameter choice
EventParameters & central(){
return central_;
}
- //! Parameter (scale) variations
+ //! Parameter (scale) variations (const version)
std::vector<EventParameters> const & variations() const{
return variations_;
}
//! Parameter (scale) variations
std::vector<EventParameters> & variations(){
return variations_;
}
- //! Parameter (scale) variation
+ //! Parameter (scale) variation (const version)
/**
* @param i Index of the requested variation
*/
EventParameters const & variations(size_t i) const{
return variations_[i];
}
//! Parameter (scale) variation
/**
* @param i Index of the requested variation
*/
EventParameters & variations(size_t i){
return variations_[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);
}
//! 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_;
}
private:
- //! \internal Event Constructor adding jet clustering to an bare, unclustered event
+ //! \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> const & incoming,
- std::vector<Particle> const & outgoing,
- std::unordered_map<size_t, std::vector<Particle>> const & decays,
- EventParameters const & central,
- std::vector<EventParameters> const & variations,
+ std::array<Particle, 2> && incoming,
+ std::vector<Particle> && outgoing,
+ std::unordered_map<size_t, std::vector<Particle>> && decays,
+ EventParameters && central,
+ std::vector<EventParameters> && variations,
fastjet::JetDefinition const & jet_def,
double const min_jet_pt
): incoming_{incoming},
outgoing_{outgoing},
decays_{decays},
central_{central},
variations_{variations},
cs_{ to_PseudoJet( filter_partons(outgoing_) ), jet_def },
min_jet_pt_{min_jet_pt}
{};
- //! \internal sort particles
- void sort();
-
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
std::unordered_map<size_t, std::vector<Particle>> decays_;
//! @TODO replace this by "ParameterVariations"
EventParameters central_;
//! @TODO replace this by "ParameterVariations"
std::vector<EventParameters> variations_;
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,
EventParameters const & central,
std::vector<EventParameters> const & variations
):
incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
decays_(std::move(decays)),
central_(std::move(central)), variations_(std::move(variations))
{};
//! 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,
EventParameters && central,
std::vector<EventParameters> && variations
):
incoming_(std::move(incoming)), outgoing_(std::move(outgoing)),
decays_(std::move(decays)),
central_(std::move(central)), variations_(std::move(variations))
{};
//! 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) const;
+ 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) const{
+ fastjet::JetDefinition const & jet_def, double const min_jet_pt){
return cluster(jet_def, min_jet_pt);
};
+ //! Sort particles in rapidity
+ void sort();
+
//! Get Incoming Particles
std::array<Particle, 2> const & get_incoming() const{
return incoming_;
}
//! Set Incoming Particles
void set_incoming(std::array<Particle, 2> const & in){
incoming_ = in;
}
//! Get Outgoing Particles
std::vector<Particle> const & get_outgoing() const{
return outgoing_;
}
//! Set Outgoing Particles
void set_outgoing(std::vector<Particle> const & out){
outgoing_ = out;
}
//! Get Particle decays in the format {outgoing index, decay products}
std::unordered_map<size_t, std::vector<Particle>> const & get_decays() const{
return decays_;
}
//! Set Particle decays in the format {outgoing index, decay products}
void set_decays(
std::unordered_map<size_t, std::vector<Particle>> const & decays
){
decays_ = decays;
}
//! Get Central parameter (e.g. scale) choice
EventParameters const & get_central() const{
return central_;
}
//! Set Central parameter (e.g. scale) choice
void set_central(EventParameters const & central){
central_ = central;
}
//! Get parameter variation
std::vector<EventParameters> const & get_variations() const{
return variations_;
}
//! Set parameter variation
void set_variations(std::vector<EventParameters> const & variations){
variations_ = variations;
}
private:
std::array<Particle, 2> incoming_;
std::vector<Particle> outgoing_;
std::unordered_map<size_t, std::vector<Particle>> decays_;
//! @TODO replace this by "ParameterVariations"
EventParameters central_;
//! @TODO replace this by "ParameterVariations"
std::vector<EventParameters> variations_;
}; // end class EventData
//! 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.3.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 */
};
}
diff --git a/src/Event.cc b/src/Event.cc
index 7d8e184..dd97a9f 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,405 +1,406 @@
/**
* \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Event.hh"
#include <algorithm>
#include <assert.h>
#include <numeric>
#include <utility>
#include "LHEF/LHEF.h"
#include "fastjet/JetDefinition.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;
/// @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(std::vector<Particle> const & outgoing){
bool has_AWZH_boson = false;
for(auto const & out: outgoing){
if(is_AWZH_boson(out.type)){
if(has_AWZH_boson) return false;
has_AWZH_boson = true;
}
else if(! is_parton(out.type)) return false;
}
return true;
}
template<class Iterator>
Iterator remove_AWZH(Iterator begin, Iterator end){
return std::remove_if(
begin, end, [](Particle const & p){return is_AWZH_boson(p);}
);
}
template<class Iterator>
bool valid_outgoing(Iterator begin, Iterator end){
return std::distance(begin, end) >= 2
&& std::is_sorted(begin, end, rapidity_less{})
&& std::count_if(
begin, end, [](Particle const & s){return is_AWZH_boson(s);}
) < 2;
}
/// @note that this changes the outgoing range!
template<class ConstIterator, class Iterator>
bool is_FKL(
ConstIterator begin_incoming, ConstIterator end_incoming,
Iterator begin_outgoing, Iterator end_outgoing
){
assert(std::distance(begin_incoming, end_incoming) == 2);
assert(std::distance(begin_outgoing, end_outgoing) >= 2);
// One photon, W, H, Z in the final state is allowed.
// Remove it for remaining tests,
end_outgoing = remove_AWZH(begin_outgoing, end_outgoing);
// Test if this is a standard FKL configuration.
return
(begin_incoming->type == begin_outgoing->type)
&& ((end_incoming-1)->type == (end_outgoing-1)->type)
&& std::all_of(
begin_outgoing + 1, end_outgoing - 1,
[](Particle const & p){ return p.type == pid::gluon; }
);
}
bool is_FKL(
std::array<Particle, 2> const & incoming,
std::vector<Particle> outgoing
){
assert(std::is_sorted(begin(incoming), end(incoming), pz_less{}));
assert(valid_outgoing(begin(outgoing), end(outgoing)));
return is_FKL(
begin(incoming), end(incoming),
begin(outgoing), end(outgoing)
);
}
bool has_2_jets(Event const & event){
return event.jets().size() >= 2;
}
/**
* \brief Checks whether event is unordered backwards
* @param ev Event
* @returns Is Event Unordered Backwards
*
* - Checks there is more than 3 constuents in the final state
* - Checks there is more than 3 jets
* - Checks the most backwards parton is a gluon
* - Checks the most forwards jet is not a gluon
* - Checks the rest of the event is FKL
* - Checks the second most backwards is not a different boson
* - Checks the unordered gluon actually forms a jet
*/
bool is_unordered_backward(Event const & ev){
auto const & in = ev.incoming();
auto const & out = ev.outgoing();
assert(std::is_sorted(begin(in), end(in), pz_less{}));
assert(valid_outgoing(begin(out), end(out)));
if(out.size() < 3) return false;
if(ev.jets().size() < 3) return false;
if(in.front().type == pid::gluon) return false;
if(out.front().type != pid::gluon) return false;
// When skipping the unordered emission
// the remainder should be a regular FKL event,
// except that the (new) first outgoing particle must not be a A,W,Z,H.
const auto FKL_begin = next(begin(out));
if(is_AWZH_boson(*FKL_begin)) return false;
if(!is_FKL(in, {FKL_begin, end(out)})) return false;
// check that the unordered gluon forms an extra jet
const auto jets = sorted_by_rapidity(ev.jets());
const auto indices = ev.particle_jet_indices({jets.front()});
return indices[0] >= 0 && indices[1] == -1;
}
/**
* \brief Checks for a forward unordered gluon emission
* @param ev Event
* @returns Is the event a forward unordered emission
*
* \see is_unordered_backward
*/
bool is_unordered_forward(Event const & ev){
auto const & in = ev.incoming();
auto const & out = ev.outgoing();
assert(std::is_sorted(begin(in), end(in), pz_less{}));
assert(valid_outgoing(begin(out), end(out)));
if(out.size() < 3) return false;
if(ev.jets().size() < 3) return false;
if(in.back().type == pid::gluon) return false;
if(out.back().type != pid::gluon) return false;
// When skipping the unordered emission
// the remainder should be a regular FKL event,
// except that the (new) last outgoing particle must not be a A,W,Z,H.
const auto FKL_end = prev(end(out));
if(is_AWZH_boson(*prev(FKL_end))) return false;
if(!is_FKL(in, {begin(out), FKL_end})) return false;
// check that the unordered gluon forms an extra jet
const auto jets = sorted_by_rapidity(ev.jets());
const auto indices = ev.particle_jet_indices({jets.back()});
return indices.back() >= 0 && indices[indices.size()-2] == -1;
}
using event_type::EventType;
EventType classify(Event const & ev){
if(! final_state_ok(ev.outgoing())) return EventType::bad_final_state;
if(! has_2_jets(ev)) return EventType::no_2_jets;
if(is_FKL(ev.incoming(), ev.outgoing())) return EventType::FKL;
if(is_unordered_backward(ev)){
return EventType::unordered_backward;
}
if(is_unordered_forward(ev)){
return EventType::unordered_forward;
}
return EventType::nonHEJ;
}
//@}
Particle extract_particle(LHEF::HEPEUP const & hepeup, int i){
return Particle{
static_cast<ParticleID>(hepeup.IDUP[i]),
fastjet::PseudoJet{
hepeup.PUP[i][0], hepeup.PUP[i][1],
hepeup.PUP[i][2], hepeup.PUP[i][3]
}
};
}
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):
central_(EventParameters{
hepeup.scales.mur, hepeup.scales.muf, hepeup.weight()
})
{
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));
}
}
//! @TODO remove in HEJ 2.3.0
Event::Event(
UnclusteredEvent const & ev,
fastjet::JetDefinition const & jet_def, double const min_jet_pt
):
Event( Event::EventData{
ev.incoming, ev.outgoing, ev.decays, ev.central, ev.variations
}.cluster(jet_def, min_jet_pt) )
{}
//! @TODO remove in HEJ 2.3.0
UnclusteredEvent::UnclusteredEvent(LHEF::HEPEUP const & hepeup){
Event::EventData const evData{hepeup};
incoming = evData.get_incoming();
outgoing = evData.get_outgoing();
decays = evData.get_decays();
central = evData.get_central();
variations = evData.get_variations();
}
- void Event::sort(){
+ 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());
}
}
Event Event::EventData::cluster(
fastjet::JetDefinition const & jet_def, double const min_jet_pt
- ) const{
- Event ev{ incoming_, outgoing_, decays_, central_, variations_,
+ ){
+ sort();
+ Event ev{ std::move(incoming_), std::move(outgoing_), std::move(decays_),
+ std::move(central_), std::move(variations_),
jet_def, min_jet_pt
};
-
- ev.sort();
- assert(std::is_sorted(begin(ev.outgoing_), end(ev.outgoing_), rapidity_less{}));
+ assert(std::is_sorted(begin(ev.outgoing_), end(ev.outgoing_),
+ rapidity_less{}));
ev.type_ = classify(ev);
return ev;
}
std::vector<fastjet::PseudoJet> Event::jets() const{
return cs_.inclusive_jets(min_jet_pt_);
}
double shat(Event const & ev){
return (ev.incoming()[0].p + ev.incoming()[1].p).m2();
}
namespace{
// colour flow according to Les Houches standard
// TODO: stub
std::vector<std::pair<int, int>> colour_flow(
std::array<Particle, 2> const & incoming,
std::vector<Particle> const & outgoing
){
std::vector<std::pair<int, int>> result(
incoming.size() + outgoing.size()
);
for(auto & col: result){
col = std::make_pair(-1, -1);
}
return result;
}
}
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()+1; // event ID
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;
for(Particle const & in: event.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);
}
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);
}
result.ICOLUP = colour_flow(
event.incoming(), filter_partons(event.outgoing())
);
if(result.ICOLUP.size() < num_particles){
const size_t AWZH_boson_idx = std::find_if(
begin(event.outgoing()), end(event.outgoing()),
[](Particle const & s){ return is_AWZH_boson(s); }
) - begin(event.outgoing()) + event.incoming().size();
assert(AWZH_boson_idx <= result.ICOLUP.size());
result.ICOLUP.insert(
begin(result.ICOLUP) + AWZH_boson_idx,
std::make_pair(0,0)
);
}
for(auto const & decay: event.decays()){
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;
}
}