Page MenuHomeHEPForge
Files F8723456

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/Event.hh b/include/RHEJ/Event.hh
index 0a5eb17..8fc6aa6 100644
--- a/include/RHEJ/Event.hh
+++ b/include/RHEJ/Event.hh
@@ -1,163 +1,166 @@
/** \file
* \brief Declares the Event class and helpers
*
*/
#pragma once
#include <string>
+#include <memory>
#include <unordered_map>
#include "RHEJ/utility.hh"
#include "RHEJ/event_types.hh"
#include "LHEF/LHEF.h"
#include "fastjet/JetDefinition.hh"
#include "fastjet/ClusterSequence.hh"
namespace RHEJ{
//! Event parameters
struct EventParameters{
double mur; /**< Value of the Renormalisation Scale */
double muf; /**< Value of the Factorisation Scale */
double weight; /**< Event Weight */
+ //! Optional description (e.g. scale names)
+ std::shared_ptr<std::string> description = nullptr;
};
//! An event before jet clustering
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<int, std::vector<Particle>> decays;
//! Central parameter (e.g. scale) choice
EventParameters central;
std::vector<EventParameters> variations; /**< For parameter variation */
};
/** An event with clustered jets
*
* This is the main reversed HEJ event class.
* It contains kinematic information including jet clustering,
* parameter (e.g. scale) settings and the event weight.
*/
class Event{
public:
//! Default Event Constructor
Event() = default;
//! Event Constructor adding jet clustering to an unclustered event
Event(
UnclusteredEvent 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
UnclusteredEvent const & unclustered() const {
return ev_;
}
//! Central parameter choice
EventParameters const & central() const{
return ev_.central;
}
//! Central parameter choice
EventParameters & central(){
return ev_.central;
}
//! Incoming particles
std::array<Particle, 2> const & incoming() const{
return ev_.incoming;
}
//! Outgoing particles
std::vector<Particle> const & outgoing() const{
return ev_.outgoing;
}
//! Particle decays
/**
* The key in the returned map corresponds to the index in the
* vector returned by outgoing()
*/
std::unordered_map<int, std::vector<Particle>> const & decays() const{
return ev_.decays;
}
//! Parameter (scale) variations
std::vector<EventParameters> const & variations() const{
return ev_.variations;
}
//! Parameter (scale) variations
std::vector<EventParameters> & variations(){
return ev_.variations;
}
//! Parameter (scale) variation
/**
* @param i Index of the requested variation
*/
EventParameters const & variations(size_t i) const{
return ev_.variations[i];
}
//! Parameter (scale) variation
/**
* @param i Index of the requested variation
*/
EventParameters & variations(size_t i){
return ev_.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:
UnclusteredEvent ev_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
event_type::EventType type_;
};
//! Square of the partonic centre-of-mass energy
double shat(Event const & ev);
//! Convert an event to a LHEF::HEPEUP
LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *);
}
diff --git a/include/RHEJ/ScaleFunction.hh b/include/RHEJ/ScaleFunction.hh
index 2f0a836..c7de1b9 100644
--- a/include/RHEJ/ScaleFunction.hh
+++ b/include/RHEJ/ScaleFunction.hh
@@ -1,145 +1,153 @@
/** \file
* \brief Functions to calculate the (renormalisation and factorisation) scales for an event
*/
#pragma once
#include <functional>
#include <vector>
#include <string>
+#include <memory>
namespace RHEJ{
class Event;
//! Class to calculate the scale associated with an event
class ScaleFunction {
public:
//! Constructor
/**
* @param name Name of the scale choice (e.g. H_T)
* @param fun Function used to calculate the scale
*/
ScaleFunction(std::string name, std::function<double(Event const &)> fun):
name_{std::move(name)},
fun_{std::move(fun)}
{}
//! Name of the scale choice
std::string const & name() const {
return name_;
}
//! Calculate the scale associated with an event
double operator()(Event const & ev) const {
return fun_(ev);
}
private:
friend ScaleFunction operator*(double factor, ScaleFunction base_scale);
friend ScaleFunction operator/(ScaleFunction base_scale, double denom);
std::string name_;
std::function<double(Event const &)> fun_;
};
//! Multiply a scale choice by a constant factor
/**
* For example, multiplying 0.5 and a scale function for H_T
* will result in a scale function for H_T/2.
*/
ScaleFunction operator*(double factor, ScaleFunction base_scale);
//! Divide a scale choice by a constant factor
/**
* For example, dividing a scale function for H_T by 2
* will result in a scale function for H_T/2.
*/
ScaleFunction operator/(ScaleFunction base_scale, double denom);
//! Calculate \f$H_T\f$ for the input event
/**
* \f$H_T\f$ is the sum of the (scalar) transverse momenta of all
* final-state particles
*/
double H_T(Event const &);
//! The maximum of all (scalar) jet transverse momentum
double max_jet_pt(Event const &);
//! The invariant mass of the sum of all jet momenta
double jet_invariant_mass(Event const &);
//! Invariant mass of the two hardest jets
double m_j1j2(Event const &);
//! Functor that returns a fixed scale regardless of the input event
class FixedScale {
public:
explicit FixedScale(double mu): mu_{mu} {}
double operator()(Event const &) const {
return mu_;
}
private:
double mu_;
};
//! Generate combinations of renormalisation and factorisation scales
class ScaleGenerator{
public:
ScaleGenerator() = default;
/** \brief Constructor
* @param scale_functions_begin Iterator to first base scale
* @param scale_functions_end Iterator past last base scale
* @param scale_factors_begin Iterator to first scale factor
* @param scale_factors_end Iterator past last scale factor
* @param max_scale_ratio Maximum ratio between renormalisation
* and factorisation scale
*/
template<class ScaleFunIterator, class FactorIterator>
ScaleGenerator(
ScaleFunIterator scale_functions_begin,
ScaleFunIterator scale_functions_end,
FactorIterator scale_factors_begin,
FactorIterator scale_factors_end,
double max_scale_ratio
):
scales_(scale_functions_begin, scale_functions_end),
scale_factors_(scale_factors_begin, scale_factors_end),
max_scale_ratio_{max_scale_ratio}
- {}
+ {
+ gen_scale_names();
+ }
/** \brief Constructor
* @param scales Base scales
* @param scale_factors Factors to multiply the base scales
* @param max_scale_ratio Maximum ratio between renormalisation
* and factorisation scale
*/
ScaleGenerator(
std::vector<ScaleFunction> scales,
std::vector<double> scale_factors,
double max_scale_ratio
):
scales_(std::move(scales)),
scale_factors_(std::move(scale_factors)),
max_scale_ratio_{max_scale_ratio}
- {}
+ {
+ gen_scale_names();
+ }
/** \brief Adjust event parameters, adding scale variation
*
* The central renormalisation and factorisation scale of the returned
* event is given be the first base scale passed to the constructor.
* The scale variation (stored in event.variation()) is constructed as
* follows: For each base scale according to the arguments of the
* constructor we add one variation where both renormalisation and
* factorisation scale are set according to the current base scale.
* Then, all combinations where the base renormalisation and factorisation
* scales are multiplied by one of the scale factors are added.
* The case were both scales are multiplied by one is skipped.
* Scale combinations where the ratio is larger than the maximum scale ratio
* set in the constructor are likewise discarded.
*/
Event operator()(Event event) const;
private:
+ void gen_scale_names();
+
std::vector<ScaleFunction> scales_;
std::vector<double> scale_factors_;
+ std::vector<std::shared_ptr<std::string>> scale_names_;
double max_scale_ratio_;
};
}
diff --git a/src/ScaleFunction.cc b/src/ScaleFunction.cc
index 7390c54..3004c5f 100644
--- a/src/ScaleFunction.cc
+++ b/src/ScaleFunction.cc
@@ -1,101 +1,150 @@
#include "RHEJ/ScaleFunction.hh"
#include <numeric>
#include <cassert>
#include "RHEJ/Event.hh"
namespace RHEJ{
double H_T(Event const & ev){
double result = 0.;
for(size_t i = 0; i < ev.outgoing().size(); ++i){
auto const decay_products = ev.decays().find(i);
if(decay_products == end(ev.decays())){
result += ev.outgoing()[i].perp();
}
else{
for(auto const & particle: decay_products->second){
result += particle.perp();
}
}
}
return result;
}
double max_jet_pt(Event const & ev) {
return sorted_by_pt(ev.jets()).front().pt();
}
double jet_invariant_mass(Event const & ev) {
fastjet::PseudoJet sum;
for(const auto & jet: ev.jets()) sum+=jet;
return sum.m();
}
double m_j1j2(Event const & ev) {
const auto jets = sorted_by_pt(ev.jets());
assert(jets.size() >= 2);
return (jets[0] + jets[1]).m();
}
ScaleFunction operator*(double factor, ScaleFunction base_scale) {
base_scale.name_.insert(0, std::to_string(factor) + '*');
auto new_fun =
[factor,fun{std::move(base_scale.fun_)}](RHEJ::Event const & ev) {
return factor*fun(ev);
};
base_scale.fun_ = std::move(new_fun);
return base_scale;
}
ScaleFunction operator/(ScaleFunction base_scale, double denom) {
base_scale.name_.append('/' + std::to_string(denom));
auto new_fun =
[denom,fun{std::move(base_scale.fun_)}](RHEJ::Event const & ev) {
return fun(ev)/denom;
};
base_scale.fun_ = std::move(new_fun);
return base_scale;
}
- Event ScaleGenerator::operator()(Event ev) const {
- if(! ev.variations().empty()) {
- throw std::invalid_argument{"Input event already has scale variation"};
+ namespace {
+ std::shared_ptr<std::string> description(
+ std::string const & mur_name, std::string const & muf_name
+ ) {
+ return std::make_shared<std::string>(
+ "mur=" + mur_name + ", muf=" + muf_name
+ );
}
+ }
+
+ // TODO: significant logic duplication with operator()
+ void ScaleGenerator::gen_scale_names() {
if(scales_.empty()) {
throw std::logic_error{"Need at least one scale"};
}
+ scale_names_.emplace_back(
+ description(scales_.front().name(), scales_.front().name())
+ );
+
+ for(auto & base_scale: scales_){
+ const auto base_name = base_scale.name();
+ scale_names_.emplace_back(description(base_name, base_name));
+ //multiplicative scale variation
+ for(double mur_factor: scale_factors_){
+ const auto mur_name = std::to_string(mur_factor) + '*' + base_name;
+ for(double muf_factor: scale_factors_){
+ if(muf_factor == 1. && mur_factor == 1.) continue;
+ const auto muf_name = std::to_string(muf_factor) + '*' + base_name;
+ if(
+ mur_factor/muf_factor < 1/max_scale_ratio_
+ || mur_factor/muf_factor > max_scale_ratio_
+ ) continue;
+ scale_names_.emplace_back(description(mur_name, muf_name));
+ }
+ }
+ }
+ }
+
+ Event ScaleGenerator::operator()(Event ev) const {
+ if(! ev.variations().empty()) {
+ throw std::invalid_argument{"Input event already has scale variation"};
+ }
+ assert(!scales_.empty());
+ assert(!scale_names_.empty());
+ size_t name_idx = 0;
const double mu_central = (scales_.front())(ev);
ev.central().mur = mu_central;
ev.central().muf = mu_central;
+ assert(name_idx < scale_names_.size());
+ assert(scale_names_[name_idx]);
+ ev.central().description = scale_names_[name_idx++];
// check if we are doing scale variation at all
if(scales_.size() == 1 && scale_factors_.empty()) return ev;
for(auto & base_scale: scales_){
const double mu_base = base_scale(ev);
+ assert(name_idx < scale_names_.size());
+ assert(scale_names_[name_idx]);
ev.variations().emplace_back(
- EventParameters{mu_base, mu_base, ev.central().weight}
+ EventParameters{
+ mu_base, mu_base, ev.central().weight, scale_names_[name_idx++]
+ }
);
//multiplicative scale variation
for(double mur_factor: scale_factors_){
const double mur = mu_base*mur_factor;
for(double muf_factor: scale_factors_){
if(muf_factor == 1. && mur_factor == 1.) continue;
const double muf = mu_base*muf_factor;
if(
mur/muf < 1/max_scale_ratio_
|| mur/muf > max_scale_ratio_
) continue;
+ assert(name_idx < scale_names_.size());
+ assert(scale_names_[name_idx]);
ev.variations().emplace_back(
- EventParameters{mur, muf, ev.central().weight}
+ EventParameters{
+ mur, muf, ev.central().weight, scale_names_[name_idx++]
+ }
);
}
}
};
return ev;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 8:36 PM (16 h, 38 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242257
Default Alt Text
(15 KB)

Event Timeline