Page MenuHomeHEPForge
Files F8725642

No OneTemporary
Actions

View Options

diff --git a/analysis-plugins/include/ROOTReader_decl.hh b/analysis-plugins/include/ROOTReader_decl.hh
index 3052693..e75a033 100644
--- a/analysis-plugins/include/ROOTReader_decl.hh
+++ b/analysis-plugins/include/ROOTReader_decl.hh
@@ -1,103 +1,103 @@
#pragma once
#include <iterator>
#include "RHEJ/Event.hh"
class TTree;
namespace RHEJ_analyses{
class ROOTReader;
class EventIterator{
public:
EventIterator() = default;
EventIterator(EventIterator const & other);
EventIterator(EventIterator && other);
EventIterator & operator=(EventIterator const & other);
EventIterator & operator=(EventIterator && other);
~EventIterator() = default;
using difference_type = long long;
- using value_type = RHEJ::ClusteredEvent;
- using pointer = RHEJ::ClusteredEvent *;
- using reference = RHEJ::ClusteredEvent const &;
+ using value_type = RHEJ::Event;
+ using pointer = RHEJ::Event *;
+ using reference = RHEJ::Event const &;
using iterator_category = std::random_access_iterator_tag;
reference operator*() const;
reference operator[](difference_type n) const;
EventIterator & operator++();
EventIterator operator++(int);
EventIterator & operator--();
EventIterator operator--(int);
EventIterator & operator+=(difference_type);
EventIterator & operator-=(difference_type);
private:
friend class ROOTReader;
friend bool operator==(EventIterator const &, EventIterator const &);
friend bool operator<(EventIterator const &, EventIterator const &);
friend difference_type operator-(
EventIterator const &, EventIterator const &
);
struct BranchEntry{
int id;
int n_out;
std::vector<float> px, py, pz, E;
std::vector<int> PDG_id;
double weight;
int PDG_id_in1, PDG_id_in2;
double muf, mur;
int n_weights;
std::vector<double> weights;
int jet_algo;
double R;
double min_jet_pt;
};
explicit EventIterator(TTree* events);
void reset_branches() const;
RHEJ::UnclusteredEvent to_UnclusteredEvent(
BranchEntry const & event_entry
) const;
long long idx_;
TTree* events_;
- mutable ClusteredEvent ev_;
+ mutable Event ev_;
mutable BranchEntry entry_;
};
bool operator==(EventIterator const &, EventIterator const &);
bool operator!=(EventIterator const &, EventIterator const &);
bool operator<(EventIterator const &, EventIterator const &);
bool operator>(EventIterator const &, EventIterator const &);
bool operator<=(EventIterator const &, EventIterator const &);
bool operator>=(EventIterator const &, EventIterator const &);
EventIterator operator+(EventIterator const &, EventIterator::difference_type);
EventIterator operator+(EventIterator::difference_type, EventIterator const &);
EventIterator operator-(EventIterator const &, EventIterator::difference_type);
EventIterator::difference_type operator-(
EventIterator const &, EventIterator const &
);
class ROOTReader{
public:
using const_iterator = EventIterator;
explicit ROOTReader(TTree * events);
ROOTReader(ROOTReader const & other) = delete;
ROOTReader(ROOTReader && other) = delete;
ROOTReader & operator=(ROOTReader const & other) = delete;
ROOTReader & operator=(ROOTReader && other) = delete;
~ROOTReader() = default;
const_iterator begin() const;
const_iterator end() const;
private:
TTree * events_;
};
}
diff --git a/analysis-plugins/include/ROOTReader_impl.hh b/analysis-plugins/include/ROOTReader_impl.hh
index ad9de9a..5ac1fc6 100644
--- a/analysis-plugins/include/ROOTReader_impl.hh
+++ b/analysis-plugins/include/ROOTReader_impl.hh
@@ -1,270 +1,270 @@
#pragma once
#include "ROOTReader.hh"
#include "TTree.h"
#include "RHEJ/kinematics.hh"
#include "RHEJ/debug.hh"
inline
std::ostream & operator<<(std::ostream & o, RHEJ::Sparticle const & s){
return o << '{' << s.type << ", " << s.p << '}';
}
inline
std::ostream & operator<<(std::ostream & o, RHEJ::EventParameters const & p){
return o << '{' << p.mur << ", " << p.muf << ", " << p.weight << '}';
}
inline
std::ostream & operator<<(std::ostream & o, RHEJ::UnclusteredEvent const & ev){
o << "incoming:\n";
for(auto && p: ev.incoming) o << " " << p << '\n';
o << "outgoing:\n";
for(auto && p: ev.outgoing) o << " " << p << '\n';
o << "central parameters: ";
return o << ev.central << '\n';
if(! ev.variations.empty()){
o << "variation parameters: {" << ev.variations.front();
for(auto it = begin(ev.variations) + 1; it != end(ev.variations); ++it){
o << ", " << (*it);
}
o << "}\n";
}
return o;
}
namespace RHEJ_analyses{
using const_iterator = ROOTReader::const_iterator;
inline
ROOTReader::ROOTReader(TTree * events):
events_{events}
{
}
inline
const_iterator ROOTReader::begin() const{
const_iterator it{events_};
it.idx_ = 0;
return it;
}
inline
const_iterator ROOTReader::end() const{
const_iterator it{events_};
it.idx_ = events_->GetEntries();
return it;
}
using difference_type = EventIterator::difference_type;
using value_type = EventIterator::value_type;
using pointer = EventIterator::pointer;
using reference = EventIterator::reference;
inline
EventIterator::EventIterator(EventIterator const & other):
idx_{other.idx_},
events_{other.events_}
{
reset_branches();
}
inline
EventIterator::EventIterator(EventIterator && other):
idx_{std::move(other.idx_)},
events_{std::move(other.events_)}
{
reset_branches();
}
inline
EventIterator & EventIterator::operator=(EventIterator const & other){
idx_ = other.idx_;
events_ = other.events_;
reset_branches();
return *this;
}
inline
EventIterator & EventIterator::operator=(EventIterator && other){
idx_ = std::move(other.idx_);
events_ = std::move(other.events_);
reset_branches();
return *this;
}
inline
EventIterator::EventIterator(TTree * events):
events_{events}
{
reset_branches();
}
inline
void EventIterator::reset_branches() const{
events_->SetBranchAddress("id", &entry_.id);
events_->SetBranchAddress("nparticle", &entry_.n_out);
events_->SetBranchAddress("weight", &entry_.weight);
events_->SetBranchAddress("id1", &entry_.PDG_id_in1);
events_->SetBranchAddress("id2", &entry_.PDG_id_in2);
events_->SetBranchAddress("fac_scale", &entry_.muf);
events_->SetBranchAddress("ren_scale", &entry_.mur);
events_->SetBranchAddress("nuwgt", &entry_.n_weights);
events_->SetBranchAddress("jet_algo", &entry_.jet_algo);
events_->SetBranchAddress("R", &entry_.R);
events_->SetBranchAddress("min_jet_pt", &entry_.min_jet_pt);
}
inline
reference EventIterator::operator*() const{
// first read number of particles and weights
events_->FindBranch("nparticle")->GetEntry(idx_);
events_->FindBranch("nuwgt")->GetEntry(idx_);
entry_.px.resize(entry_.n_out);
entry_.py.resize(entry_.n_out);
entry_.pz.resize(entry_.n_out);
entry_.E.resize(entry_.n_out);
entry_.PDG_id.resize(entry_.n_out);
entry_.weights.resize(entry_.n_weights);
events_->SetBranchAddress("px", entry_.px.data());
events_->SetBranchAddress("py", entry_.py.data());
events_->SetBranchAddress("pz", entry_.pz.data());
events_->SetBranchAddress("E", entry_.E.data());
events_->SetBranchAddress("kf", entry_.PDG_id.data());
events_->SetBranchAddress("usr_wgts", entry_.weights.data());
// second read to get the whole event
events_->GetEntry(idx_);
auto ev = to_UnclusteredEvent(entry_);
fastjet::JetDefinition jet_def{
static_cast<fastjet::JetAlgorithm>(entry_.jet_algo), entry_.R
};
- ev_ = ClusteredEvent{std::move(ev), std::move(jet_def), entry_.min_jet_pt};
+ ev_ = Event{std::move(ev), std::move(jet_def), entry_.min_jet_pt};
return ev_;
}
inline
RHEJ::UnclusteredEvent EventIterator::to_UnclusteredEvent(
EventIterator::BranchEntry const & event_entry
) const{
RHEJ::UnclusteredEvent ev;
ev.outgoing.resize(event_entry.n_out);
for(size_t i = 0; i < (size_t) event_entry.n_out; ++i){
ev.outgoing[i].p = {
event_entry.px[i], event_entry.py[i],
event_entry.pz[i], event_entry.E[i]
};
ev.outgoing[i].type = static_cast<RHEJ::ParticleID>(event_entry.PDG_id[i]);
}
std::tie(ev.incoming[0].p, ev.incoming[1].p) =
RHEJ::incoming_momenta(ev.outgoing);
ev.incoming[0].type = static_cast<RHEJ::ParticleID>(event_entry.PDG_id_in1);
ev.incoming[1].type = static_cast<RHEJ::ParticleID>(event_entry.PDG_id_in2);
ev.central = {event_entry.mur, event_entry.muf, event_entry.weight};
ev.variations.resize(event_entry.n_weights);
for(size_t i = 0; i < (size_t) event_entry.n_weights; ++i){
ev.variations[i].weight = event_entry.weights[i];
}
return ev;
}
inline
EventIterator & EventIterator::operator++(){
++idx_;
return *this;
}
inline
EventIterator EventIterator::operator++(int){
EventIterator res = *this;
++idx_;
return res;
}
inline
EventIterator & EventIterator::operator--(){
--idx_;
return *this;
}
inline
EventIterator EventIterator::operator--(int){
EventIterator res = *this;
--idx_;
return res;
}
inline
EventIterator & EventIterator::operator+=(difference_type n){
idx_ += n;
return *this;
}
inline
EventIterator & EventIterator::operator-=(difference_type n){
idx_ -= n;
return *this;
}
inline
reference EventIterator::operator[](difference_type n) const{
return *(*this + n);
}
inline
bool operator==(EventIterator const & a, EventIterator const & b){
return a.events_==b.events_ && a.idx_==b.idx_;
}
inline
bool operator!=(EventIterator const & a, EventIterator const & b){
return !(a==b);
}
inline
bool operator<(EventIterator const & a, EventIterator const & b){
return a.idx_ < b.idx_;
}
inline
bool operator<=(EventIterator const & a, EventIterator const & b){
return a==b || a<b;
}
inline
bool operator>(EventIterator const & a, EventIterator const & b){
return !(a<=b);
}
inline
bool operator>=(EventIterator const & a, EventIterator const & b){
return !(a<b);
}
inline
EventIterator operator+(EventIterator const & a, difference_type n){
EventIterator res = a;
return res += n;
}
inline
EventIterator operator+(difference_type n, EventIterator const & a){
return a + n;
}
inline
EventIterator operator-(EventIterator const & a, difference_type n){
return a + (-n);
}
inline
difference_type operator-(
EventIterator const & a, EventIterator const & b
){
return a.idx_ - b.idx_;
}
}
diff --git a/analysis-plugins/include/ROOTWriter_decl.hh b/analysis-plugins/include/ROOTWriter_decl.hh
index 828f7a7..48a973c 100644
--- a/analysis-plugins/include/ROOTWriter_decl.hh
+++ b/analysis-plugins/include/ROOTWriter_decl.hh
@@ -1,43 +1,43 @@
#pragma once
#include <vector>
#include "TTree.h"
namespace RHEJ{
- class ClusteredEvent;
+ class Event;
}
namespace RHEJ_analyses{
- using RHEJ::ClusteredEvent;
+ using RHEJ::Event;
class ROOTWriter{
public:
explicit ROOTWriter(char const * treename);
explicit ROOTWriter(std::string const & treename);
ROOTWriter(ROOTWriter const & other) = delete;
ROOTWriter(ROOTWriter && other) = delete;
ROOTWriter & operator=(ROOTWriter const & other) = delete;
ROOTWriter & operator=(ROOTWriter && other) = delete;
~ROOTWriter();
- void write(ClusteredEvent const & ev);
+ void write(Event const & ev);
private:
TTree events_;
int counter_;
int n_out_;
std::vector<float> px_, py_, pz_, E_;
std::vector<int> PDG_id_;
double weight_;
int PDG_id_in1_, PDG_id_in2_;
double muf_, mur_;
int n_weights_;
std::vector<double> weights_;
int jet_algo_;
double R_;
double min_jet_pt_;
};
}
diff --git a/analysis-plugins/include/ROOTWriter_impl.hh b/analysis-plugins/include/ROOTWriter_impl.hh
index 327f44a..7875d44 100644
--- a/analysis-plugins/include/ROOTWriter_impl.hh
+++ b/analysis-plugins/include/ROOTWriter_impl.hh
@@ -1,84 +1,84 @@
#pragma once
#include "RHEJ/Analysis.hh"
#include "RHEJ/Event.hh"
#include "fastjet/JetDefinition.hh"
namespace RHEJ_analyses{
inline
ROOTWriter::ROOTWriter(char const * treename):
events_{treename, treename},
counter_{-1}
{
events_.Branch("id", &counter_, "id/I");
events_.Branch("nparticle", &n_out_, "nparticle/I");
events_.Branch("px", px_.data(), "px[nparticle]/F");
events_.Branch("py", py_.data(), "py[nparticle]/F");
events_.Branch("pz", pz_.data(), "pz[nparticle]/F");
events_.Branch("E", E_.data(), "E[nparticle]/F");
events_.Branch("kf", PDG_id_.data(), "kf[nparticle]/I");
events_.Branch("weight", &weight_, "weight/D");
events_.Branch("id1", &PDG_id_in1_, "id1/I");
events_.Branch("id2", &PDG_id_in2_, "id2/I");
events_.Branch("fac_scale", &muf_, "fac_scale/D");
events_.Branch("ren_scale", &mur_, "ren_scale/D");
events_.Branch("nuwgt", &n_weights_, "nuwgt/I");
events_.Branch("usr_wgts", weights_.data(), "usr_wgts[nuwgt]/D");
events_.Branch("jet_algo", &jet_algo_, "jet_algo/I");
events_.Branch("R", &R_, "R/D");
events_.Branch("min_jet_pt", &min_jet_pt_, "min_jet_pt/D");
}
inline
ROOTWriter::ROOTWriter(std::string const & filename):
ROOTWriter{filename.c_str()}
{}
inline
ROOTWriter::~ROOTWriter(){
events_.Write();
}
inline
- void ROOTWriter::write(ClusteredEvent const & ev){
+ void ROOTWriter::write(Event const & ev){
using RHEJ::EventParameters;
++counter_;
n_out_ = ev.outgoing().size();
px_.resize(ev.outgoing().size());
py_.resize(ev.outgoing().size());
pz_.resize(ev.outgoing().size());
E_.resize(ev.outgoing().size());
PDG_id_.resize(ev.outgoing().size());
for(size_t i = 0; i < ev.outgoing().size(); ++i){
px_[i] = ev.outgoing()[i].px();
py_[i] = ev.outgoing()[i].py();
pz_[i] = ev.outgoing()[i].pz();
E_[i] = ev.outgoing()[i].E();
PDG_id_[i] = static_cast<int>(ev.outgoing()[i].type);
}
weight_ = ev.central().weight;
PDG_id_in1_ = ev.incoming().front().type;
PDG_id_in2_ = ev.incoming().back().type;
muf_ = ev.central().muf;
mur_ = ev.central().mur;
n_weights_ = ev.variations().size();
weights_.resize(ev.variations().size());
std::transform(
begin(ev.variations()), end(ev.variations()), begin(weights_),
[](EventParameters const & p){ return p.weight; }
);
jet_algo_ = static_cast<int>(ev.jet_def().jet_algorithm());
R_ = ev.jet_def().R();
min_jet_pt_ = ev.min_jet_pt();
events_.SetBranchAddress("px", px_.data());
events_.SetBranchAddress("py", py_.data());
events_.SetBranchAddress("pz", pz_.data());
events_.SetBranchAddress("E", E_.data());
events_.SetBranchAddress("kf", PDG_id_.data());
events_.SetBranchAddress("usr_wgts", weights_.data());
events_.Fill();
}
}
diff --git a/analysis-plugins/include/StandardAnalysis_decl.hh b/analysis-plugins/include/StandardAnalysis_decl.hh
index aac6cb2..9b4bf5f 100644
--- a/analysis-plugins/include/StandardAnalysis_decl.hh
+++ b/analysis-plugins/include/StandardAnalysis_decl.hh
@@ -1,143 +1,143 @@
#pragma once
#include <array>
#include <map>
#include "fastjet/JetDefinition.hh"
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "RHEJ/Analysis.hh"
#include "RHEJ/Event.hh"
#include "Histogram.hh"
class TFile;
namespace YAML{
class Node;
}
namespace RHEJ_analyses{
- using RHEJ::ClusteredEvent;
+ using RHEJ::Event;
namespace detail{
struct AbsWtLess{
- bool operator()(ClusteredEvent const & a, ClusteredEvent const & b){
+ bool operator()(Event const & a, Event const & b){
return std::abs(a.central().weight) < std::abs(b.central().weight);
}
};
}
class StandardAnalysis: public RHEJ::DynamicAnalysis{
public:
using histogram_type = Histogram;
// macro definition to avoid repetition of histogram names
#define RHEJ_HISTOGRAMS(F) F(total), F(np), F(njets),F(pt1), F(pt2), F(pth1), F(pth2), F(qt), F(ya), F(yb), F(ydif), F(yEW), F(ptEW), F(wt), F(wtwt)
#define RHEJ_AS_ENUM(VAR) VAR
#define RHEJ_AS_STRING(VAR) #VAR
enum Id{
RHEJ_HISTOGRAMS(RHEJ_AS_ENUM),
first_histogram = np,
last_histogram = wtwt,
first_histpack = np,
last_histpack = ptEW
};
static constexpr auto hist_names = RHEJ::make_array(
RHEJ_HISTOGRAMS(RHEJ_AS_STRING)
);
#undef RHEJ_HISTOGRAMS
#undef RHEJ_AS_ENUM
#undef RHEJ_AS_STRING
StandardAnalysis(YAML::Node const &);
explicit StandardAnalysis(std::string const & filename);
StandardAnalysis(StandardAnalysis const & other);
StandardAnalysis(StandardAnalysis && other);
StandardAnalysis & operator=(StandardAnalysis const & other);
StandardAnalysis & operator=(StandardAnalysis && other);
- void fill(ClusteredEvent const & e) override;
+ void fill(Event const & e) override;
std::vector<std::string> histogram_names() const override;
std::vector<Parameter> parameters() const override;
void draw_histogram(
std::string const & name, std::string const & filename
) const override;
bool update_parameter(std::string const & name, double value) override;
~StandardAnalysis() override;
static std::unique_ptr<Analysis> create(YAML::Node const &);
static std::unique_ptr<DynamicAnalysis> create(
std::string const & filename
);
private:
void store();
void restore(std::string const & filename);
void write_settings(TFile & out);
void write_histograms(TFile & out);
void write_events(TFile & out);
void read_settings(TFile & in);
void read_histograms(TFile & in);
void read_events(TFile & in);
struct hist_pack{
histogram_type h, hwt, hN;
hist_pack() = default;
hist_pack(
const char * title, const std::string & id,
int numbins, double min, double max
);
void fill(double val, double wt);
void write() const;
int nbins() const;
};
struct KinInfo{
std::vector<RHEJ::Sparticle> out;
std::vector<fastjet::PseudoJet> jets_y;
std::vector<fastjet::PseudoJet> jets_pt;
fastjet::PseudoJet q;
};
struct EventHists{
EventHists();
explicit EventHists(std::string const & id);
void fill(KinInfo const & i, double weight);
void write() const;
std::array<hist_pack, last_histpack + 1> hist_;
histogram_type hist_wt_, hist_wtwt_;
};
void restore_histograms(
std::map<std::string, histogram_type> && hists
);
void update_hist_by_name();
void add_to_hist_by_name(EventHists const & hist);
void add_to_hist_by_name(hist_pack const & pack);
void reset();
EventHists central_;
std::vector<EventHists> hist_;
double log_wt_cut_;
std::string output_;
- std::set<ClusteredEvent, detail::AbsWtLess> bad_events_;
+ std::set<Event, detail::AbsWtLess> bad_events_;
using Hist_ref = std::reference_wrapper<const histogram_type>;
std::map<std::string, Hist_ref> hist_by_name_;
EventHists orig_central_;
std::vector<EventHists> orig_hist_;
double orig_log_wt_cut_;
};
}
diff --git a/analysis-plugins/include/StandardAnalysis_impl.hh b/analysis-plugins/include/StandardAnalysis_impl.hh
index 1e06d04..3a0b19b 100644
--- a/analysis-plugins/include/StandardAnalysis_impl.hh
+++ b/analysis-plugins/include/StandardAnalysis_impl.hh
@@ -1,525 +1,525 @@
#pragma once
#include <iostream>
#include <numeric>
#include "StandardAnalysis.hh"
#include "ROOTWriter.hh"
#include "ROOTReader.hh"
#include "TDirectory.h"
#include "TKey.h"
#include "TCanvas.h"
#include "TImage.h"
#include "TFile.h"
#include "yaml-cpp/yaml.h"
namespace RHEJ_analyses{
namespace detail{
inline
bool disable_TH1_reference(){
TH1::AddDirectory(kFALSE);
std::cout << "Disabled central ROOT histogram references\n";
return true;
};
}
using histogram = StandardAnalysis::histogram_type;
inline
StandardAnalysis::hist_pack::hist_pack(
const char * title, const std::string & id,
int numbins, double min, double max
):
h((title + id).c_str(), title, numbins, min, max),
hwt(
(title + ("wt2" + id)).c_str(),
(title + std::string("wt2")).c_str(),
numbins, min, max
),
hN(
(title + ("N" + id)).c_str(),
(title + std::string("N")).c_str(),
numbins, min, max
)
{
static const bool TH1_reference_disabled = detail::disable_TH1_reference();
(void) TH1_reference_disabled; // avoid compiler warnings
h.sumw2();
}
inline
void StandardAnalysis::hist_pack::write() const{
h.write();
hwt.write();
hN.write();
}
inline
int StandardAnalysis::hist_pack::nbins() const{
return h.nbins();
}
inline
void StandardAnalysis::hist_pack::fill(double val, double wt){
const int binnumber = h.x_axis()->FindBin(val);
const double binwidth = h.x_axis()->GetBinWidth(binnumber);
h.fill(val, wt/binwidth);
hwt.fill(val, wt*wt/binwidth);
hN.fill(val, 1);
}
inline
StandardAnalysis::EventHists::EventHists():
EventHists{""} {}
StandardAnalysis::EventHists::EventHists(std::string const & id){
#define NEWHIST(X,UUID,BINS,MIN,MAX) hist_[X] = hist_pack{#X, UUID, BINS, MIN, MAX};
NEWHIST(total,id,1,-.5,.5);
NEWHIST(np,id,35,-.5,34.5);
NEWHIST(njets,id,7,.5,7.5);
NEWHIST(pt1,id,100,0.,1000.);
NEWHIST(pt2,id,100,0.,1000.);
NEWHIST(pth1,id,100,0.,1000.);
NEWHIST(pth2,id,100,0.,1000.);
NEWHIST(qt,id,100,0.,99.);
NEWHIST(ya,id,40,-5.,5.);
NEWHIST(yb,id,40,-5.,5.);
NEWHIST(ydif,id,32,0.,8.);
NEWHIST(yEW,id,40,-5.,5.);
NEWHIST(ptEW,id,20,0.,200.);
#undef NEWHIST
hist_wt_ = histogram{("wt" + id).c_str(),"wt",800,-30.,50.};
hist_wtwt_ = histogram{("wtwt" + id).c_str(),"wtwt",800,-30.,50.};
}
inline
void StandardAnalysis::EventHists::fill(KinInfo const & i, double weight){
const double ymin = i.jets_y.front().rapidity();
const double ymax = i.jets_y.back().rapidity();
assert(ymax > ymin);
const auto EW_boson = std::find_if(
begin(i.out), end(i.out),
[](RHEJ::Sparticle const & s){ return is_AWZH_boson(s); }
);
const size_t npartons = (EW_boson==end(i.out))?
i.out.size():
i.out.size()-1
;
hist_[total].fill(0., weight);
hist_[np].fill(npartons, weight);
for(size_t jets = 2; jets <= i.jets_y.size(); ++jets){
hist_[njets].fill(jets, weight);
}
hist_[pt1].fill(i.jets_y.front().perp(), weight);
hist_[pt2].fill(i.jets_y.back().perp(), weight);
hist_[pth1].fill(i.jets_pt.front().perp(), weight);
hist_[pth2].fill(i.jets_pt[1].perp(), weight);
hist_[qt].fill(i.q.perp(), weight);
hist_[ya].fill(ymin, weight);
hist_[yb].fill(ymax, weight);
hist_[ydif].fill(ymax - ymin, weight);
if(EW_boson!=end(i.out)){
hist_[yEW].fill(EW_boson->rapidity(), weight);
hist_[ptEW].fill(EW_boson->perp(), weight);
}
const double awt = std::abs(weight);
hist_wt_.fill(std::log(awt), 1);
hist_wtwt_.fill(std::log(awt), awt);
}
inline
void StandardAnalysis::EventHists::write() const{
for(auto hist: hist_){
hist.write();
}
hist_wt_.write();
hist_wtwt_.write();
}
namespace detail{
inline
double get_log_wt_cut(YAML::Node const & cut){
return cut?
cut.as<double>():
std::numeric_limits<double>::infinity();
}
}
inline
StandardAnalysis::StandardAnalysis(YAML::Node const & parameters):
log_wt_cut_{detail::get_log_wt_cut(parameters["log wt cut"])},
output_{parameters["output"].as<std::string>()}
{}
inline
StandardAnalysis::StandardAnalysis(
std::string const & filename
){
restore(filename);
}
inline
StandardAnalysis::~StandardAnalysis(){
if(!output_.empty()) store();
}
inline
StandardAnalysis::StandardAnalysis(StandardAnalysis const & other){
*this = other;
}
inline
StandardAnalysis::StandardAnalysis(StandardAnalysis && other){
*this = std::move(other);
}
inline
StandardAnalysis & StandardAnalysis::operator=(StandardAnalysis const & other){
if(this == &other) return *this;
central_ = other.central_;
hist_ = other.hist_;
log_wt_cut_ = other.log_wt_cut_;
output_ = other.output_;
bad_events_ = other.bad_events_;
orig_central_ = other.orig_central_;
orig_hist_ = other.orig_hist_;
orig_log_wt_cut_ = other.orig_log_wt_cut_;
update_hist_by_name();
return *this;
}
StandardAnalysis & StandardAnalysis::operator=(StandardAnalysis && other){
if(this == &other) return *this;
central_ = std::move(other.central_);
hist_ = std::move(other.hist_);
log_wt_cut_ = std::move(other.log_wt_cut_);
output_ = std::move(other.output_);
bad_events_ = std::move(other.bad_events_);
orig_central_ = std::move(other.orig_central_);
orig_hist_ = std::move(other.orig_hist_);
orig_log_wt_cut_ = std::move(other.orig_log_wt_cut_);
update_hist_by_name();
return *this;
}
inline
std::vector<std::string> StandardAnalysis::histogram_names() const{
std::vector<std::string> result;
for(auto const & entry: hist_by_name_) result.emplace_back(entry.first);
return result;
};
inline
std::vector<RHEJ::DynamicAnalysis::Parameter> StandardAnalysis::parameters() const{
double min = log_wt_cut_;
double max = log_wt_cut_;
if(! bad_events_.empty()){
min = std::min(min, std::log(std::abs(begin(bad_events_)->central().weight)));
max = std::max(max, std::log(std::abs(rbegin(bad_events_)->central().weight)));
}
return {RHEJ::DynamicAnalysis::Parameter{"log wt cut", log_wt_cut_, min, max}};
}
inline
void StandardAnalysis::draw_histogram(
std::string const & name, std::string const & filename
) const{
const auto hist = hist_by_name_.find(name);
if(hist == end(hist_by_name_)){
throw std::invalid_argument{
"histogram with title "+ name + " does not exist"
};
}
TCanvas c;
hist->second.get().draw();
const auto img = std::unique_ptr<TImage>{TImage::Create()};
img->FromPad(&c);
img->WriteImage(filename.c_str(), TImage::EImageFileTypes::kPng);
}
inline
bool StandardAnalysis::update_parameter(std::string const & name, double value){
if(name != "log wt cut"){
throw std::invalid_argument{"invalid parameter " + name};
}
- ClusteredEvent search_event;
+ Event search_event;
search_event.central().weight = std::exp(log_wt_cut_);
const auto old_end = bad_events_.lower_bound(search_event);
search_event.central().weight = std::exp(value);
const auto new_end = bad_events_.lower_bound(search_event);
if(new_end == old_end){
log_wt_cut_ = value;
return false; // no update required
}
if(value < log_wt_cut_){
// we would have to remove events from histograms
// instead, reset and add all events again
reset();
update_parameter(name, value);
return true;
}
log_wt_cut_ = value;
for(auto it = old_end; it != new_end; ++it) fill(*it);
return true;
}
inline
- void StandardAnalysis::fill(ClusteredEvent const & ev){
+ void StandardAnalysis::fill(Event const & ev){
if(std::log(std::abs(ev.central().weight)) > log_wt_cut_){
bad_events_.emplace(ev);
return;
}
const std::vector<fastjet::PseudoJet> partons =
to_PseudoJet(filter_partons(ev.outgoing()));
auto jets_y = sorted_by_rapidity(ev.jets());
auto q =
std::accumulate(partons.begin(), partons.end(), fastjet::PseudoJet{})
- std::accumulate(jets_y.begin(), jets_y.end(), fastjet::PseudoJet{})
;
const KinInfo info = {
ev.outgoing(),
std::move(jets_y),
sorted_by_pt(ev.jets()),
std::move(q)
};
assert(info.jets_y.size() == info.jets_pt.size());
assert(info.jets_y.size() >= 2);
central_.fill(info, ev.central().weight);
if(hist_.empty()){
hist_.reserve(ev.variations().size());
for(size_t i = 0; i < ev.variations().size(); ++i){
hist_.emplace_back('_' + std::to_string(i));
}
}
if(hist_.size() != ev.variations().size()){
throw std::invalid_argument(
"number of weights " + std::to_string(ev.variations().size())
+ " does not match number of histograms "
+ std::to_string(hist_.size())
);
}
for(size_t i = 0; i < ev.variations().size(); ++i){
hist_[i].fill(info, ev.variations(i).weight);
}
}
inline
void StandardAnalysis::write_settings(TFile &){
TTree settings{"settings", "settings"};
settings.Branch("log_wt_cut", &log_wt_cut_, "wtwt_cut/D");
settings.Fill();
settings.Write();
}
inline
void StandardAnalysis::read_settings(TFile & in){
auto tree = static_cast<TTree*>(in.Get("settings"));
if(tree == nullptr){
throw std::runtime_error("no TTree with name \"settings\"");
}
tree->SetBranchAddress("log_wt_cut", &log_wt_cut_);
if(tree->GetEntries() != 1){
throw std::runtime_error("wrong number of entries in \"settings\"");
}
tree->GetEntry(0);
}
inline
void StandardAnalysis::write_histograms(TFile & out){
const auto dir = out.mkdir("histograms");
assert(dir != nullptr);
dir->cd();
central_.write();
for(auto const & hist: hist_) hist.write();
out.cd();
}
inline
void StandardAnalysis::write_events(TFile &){
ROOTWriter writer{"events"};
for(auto && event: bad_events_) writer.write(event);
}
namespace detail{
struct weight_less{
- bool operator()(ClusteredEvent const & a, ClusteredEvent const & b){
+ bool operator()(Event const & a, Event const & b){
return a.central().weight < b.central().weight;
}
- bool operator()(ClusteredEvent const & ev, double wt){
+ bool operator()(Event const & ev, double wt){
return ev.central().weight < wt;
}
- bool operator()(double wt, RHEJ::ClusteredEvent const & ev){
+ bool operator()(double wt, RHEJ::Event const & ev){
return wt < ev.central().weight;
}
};
}
inline
void StandardAnalysis::read_events(TFile & in){
assert(in.IsOpen());
auto event_tree = static_cast<TTree*>(in.Get("events"));
if(event_tree == nullptr){
throw std::runtime_error("no TTree with name \"events\"");
}
ROOTReader event_reader{event_tree};
- bad_events_ = std::set<ClusteredEvent, detail::AbsWtLess>(
+ bad_events_ = std::set<Event, detail::AbsWtLess>(
std::begin(event_reader), std::end(event_reader)
);
}
inline
void StandardAnalysis::store(){
TFile out{output_.c_str(), "RECREATE"};
write_settings(out);
write_histograms(out);
write_events(out);
}
inline
void StandardAnalysis::restore(std::string const & filename){
TFile in{filename.c_str()};
read_settings(in);
read_histograms(in);
read_events(in);
orig_central_ = central_;
orig_hist_ = hist_;
orig_log_wt_cut_ = log_wt_cut_;
}
namespace detail{
using histogram_type = StandardAnalysis::histogram_type;
inline
size_t num_hist_indices(
std::map<std::string, Histogram> const & hists,
std::string const & name
){
size_t result = 0;
while(hists.find(name + "_" + std::to_string(result)) != hists.end()){
++result;
}
return result;
}
inline
std::map<std::string, histogram_type> read_hist_map(TDirectory & in){
TIter keys = in.GetListOfKeys();
TKey* key;
std::map<std::string, histogram_type> hists;
while( (key = static_cast<TKey*>(keys())) ){
auto new_hist = dynamic_cast<TH1D*>(key->ReadObj());
if(! new_hist){
throw std::runtime_error("Encountered unexpected object in ROOT file");
}
hists.emplace(new_hist->GetName(), histogram_type{new_hist});
}
return hists;
}
}
inline
void StandardAnalysis::read_histograms(TFile & in){
const auto dir = in.GetDirectory("histograms");
if(dir == nullptr){
throw std::runtime_error{"no directory \"histograms\""};
}
restore_histograms(detail::read_hist_map(*dir));
update_hist_by_name();
}
inline
void StandardAnalysis::restore_histograms(
std::map<std::string, histogram_type> && hists
){
hist_.resize(detail::num_hist_indices(hists, hist_names[first_histogram]));
for(size_t id = first_histogram; id <= last_histpack; ++id){
central_.hist_[id].h = std::move(hists[hist_names[id]]);
for(size_t i = 0; i < hist_.size(); ++i){
const std::string name = hist_names[id] + ("_" + std::to_string(i));
hist_[i].hist_[id].h = std::move(hists[name]);
}
std::string name = hist_names[id] + std::string{"wt2"};
central_.hist_[id].hwt = std::move(hists[name]);
for(size_t i = 0; i < hist_.size(); ++i){
const std::string name = hist_names[id] + ("wt2_" + std::to_string(i));
hist_[i].hist_[id].hwt = std::move(hists[name]);
}
name = hist_names[id] + std::string{"N"};
central_.hist_[id].hN = std::move(hists[name]);
for(size_t i = 0; i < hist_.size(); ++i){
const std::string name = hist_names[id] + ("N_" + std::to_string(i));
hist_[i].hist_[id].hN = std::move(hists[name]);
}
}
central_.hist_wt_ = std::move(hists["wt"]);
central_.hist_wtwt_ = std::move(hists["wtwt"]);
for(size_t i = 0; i < hist_.size(); ++i){
std::string name = "wt_" + std::to_string(i);
hist_[i].hist_wt_ = std::move(hists[name]);
name = "wtwt_" + std::to_string(i);
hist_[i].hist_wtwt_ = std::move(hists[name]);
}
}
inline
void StandardAnalysis::update_hist_by_name(){
hist_by_name_.clear();
add_to_hist_by_name(central_);
for(auto const & hist: hist_) add_to_hist_by_name(hist);
}
inline
void StandardAnalysis::add_to_hist_by_name(EventHists const & hist){
for(auto const & h: hist.hist_) add_to_hist_by_name(h);
hist_by_name_.emplace(hist.hist_wt_.title(), Hist_ref{hist.hist_wt_});
hist_by_name_.emplace(hist.hist_wtwt_.title(), Hist_ref{hist.hist_wtwt_});
}
inline
void StandardAnalysis::add_to_hist_by_name(hist_pack const & pack){
hist_by_name_.emplace(pack.h.title(), Hist_ref{pack.h});
hist_by_name_.emplace(pack.hwt.title(), Hist_ref{pack.hwt});
hist_by_name_.emplace(pack.hN.title(), Hist_ref{pack.hN});
}
inline
std::unique_ptr<RHEJ::Analysis> StandardAnalysis::create(
YAML::Node const & parameters
){
return std::unique_ptr<Analysis>{new StandardAnalysis{parameters}};
}
inline
std::unique_ptr<RHEJ::DynamicAnalysis> StandardAnalysis::create(
std::string const & filename
){
return std::unique_ptr<RHEJ::DynamicAnalysis>{
new StandardAnalysis{filename}
};
}
inline
void StandardAnalysis::reset(){
central_ = orig_central_;
hist_ = orig_hist_;
log_wt_cut_ = orig_log_wt_cut_;
update_hist_by_name();
}
constexpr
decltype(StandardAnalysis::hist_names) StandardAnalysis::hist_names;
}
diff --git a/analysis-plugins/src/VBF.cc b/analysis-plugins/src/VBF.cc
index b472bb2..95e869f 100644
--- a/analysis-plugins/src/VBF.cc
+++ b/analysis-plugins/src/VBF.cc
@@ -1,93 +1,93 @@
#include "StandardAnalysis.hh"
#include <iostream>
namespace{
using namespace RHEJ_analyses;
class VBF: RHEJ::DynamicAnalysis{
public:
VBF(YAML::Node const & parameters):
analysis_{parameters},
min_dy12_{parameters["min dy12"].as<double>()},
min_m12_{parameters["min m12"].as<double>()}
{}
VBF(std::string const & filename):
analysis_{filename},
/**
* The following lines are only safe if "filename" was actually
* generated using this analysis.
*
* In this case all stored event already fulfil the VBF cuts
* and we can ignore them
*/
min_dy12_{0.},
min_m12_{0.}
{}
- void fill(ClusteredEvent const & ev) override{
+ void fill(Event const & ev) override{
const auto jets_pt = sorted_by_pt(ev.jets());
// tight cuts a la arXiv:1311.6738
if(jets_pt.size() < 2) return;
const double y0 = jets_pt[0].rapidity();
const double y1 = jets_pt[1].rapidity();
if(y0*y1 >= 0) return;
if(std::abs(y0 - y1) <= min_dy12_) return;
if((jets_pt[0] + jets_pt[1]).m2() <= min_m12_*min_m12_) return;
analysis_.fill(ev);
}
std::vector<std::string> histogram_names() const override{
return analysis_.histogram_names();
}
std::vector<Parameter> parameters() const override{
return analysis_.parameters();
}
void draw_histogram(
std::string const & name, std::string const & filename
) const override{
analysis_.draw_histogram(name, filename);
}
bool update_parameter(std::string const & name, double value) override{
return analysis_.update_parameter(name, value);
}
static std::unique_ptr<RHEJ::Analysis> create(
YAML::Node const & parameters
){
return std::unique_ptr<RHEJ::Analysis>{new VBF{parameters}};
}
static std::unique_ptr<DynamicAnalysis> create(
std::string const & filename
){
return std::unique_ptr<RHEJ::DynamicAnalysis>{new VBF{filename}};
}
private:
StandardAnalysis analysis_;
double min_dy12_;
double min_m12_;
};
}
extern "C"
__attribute__((visibility("default")))
std::unique_ptr<RHEJ::Analysis> make_analysis(
YAML::Node const & parameters
){
return VBF::create(parameters);
}
extern "C"
__attribute__((visibility("default")))
std::unique_ptr<RHEJ::Analysis> load_analysis(
std::string const & filename
){
return VBF::create(filename);
}
diff --git a/include/RHEJ/Analysis.hh b/include/RHEJ/Analysis.hh
index f2d1177..04d7eef 100644
--- a/include/RHEJ/Analysis.hh
+++ b/include/RHEJ/Analysis.hh
@@ -1,86 +1,86 @@
/** \file Analysis.hh
* \brief Header file for the Analysis Interface (ABC)
*
* Contains the Analysis struct and the Dynamic Analysis
* struct. Also contains some histogram centric functions
* such as histogram_names(), parameters(), draw_histogram()
* update_parameter()
*/
#pragma once
#include <vector>
#include <string>
//! Reversed HEJ Namespace
namespace RHEJ{
- class ClusteredEvent;
+ class Event;
/** \struct Analysis Analysis.hh "include/RHEJ/Analysis.hh"
* \brief Analysis Struct
*
* Contains a fill function which needs to be overridden by any
* inheriting class before it can be utilized.
*/
struct Analysis{
- virtual void fill(ClusteredEvent const &) = 0;
+ virtual void fill(Event const &) = 0;
virtual ~Analysis() = default;
};
/** \struct DynamicAnalysis Analysis.hh "include/RHEJ/Analysis.hh"
* \brief DynamicAnalysis struct to be used with HEJview
*
* Inherits from the standard Analysis struct.
*/
struct DynamicAnalysis : Analysis {
/** \struct Parameter Analysis.hh "include/RHEJ/Analysis.hh"
* \brief Parameter Struct which contains Parameters of Analysis
*
* Contains all info necessary about a parameter of the analysis.
*/
struct Parameter{
std::string name; /**< Name of parameter */
double value; /**< Value of Paramter */
double min; /**< Min Value allowed for parameter */
double max; /**< Max Value allowed for parameter */
};
//! Virtual Histogram Names Function
/**
* Returns a vector of the names of the histograms.
*/
virtual std::vector<std::string> histogram_names() const = 0;
//! Virtual Parameters Function
/**
* Returns a vector which contains the parameters to be placed into the
* histograms
*/
virtual std::vector<Parameter> parameters() const = 0;
//! Virtual Draw Histogram Function
/**
* @param name Name of Histogram to be drawn
* @param filename Name of Output file containing Histogram
*
* This function will update the appropriate histogram?
*/
virtual void draw_histogram(
std::string const & name, std::string const & filename
) const = 0;
//! Virtual Update Parameter Function
/**
* @param name Name of the Parameter which is to be updated
* @param value Value to be updated to
* @returns 0 or 1 depending on if value is updated?
*
* How is its decided whether the value will be updated?
*/
virtual bool update_parameter(std::string const & name, double value) = 0;
//! Dynamic Analysis Destructor
virtual ~DynamicAnalysis() = default;
};
}
diff --git a/include/RHEJ/CombinedEventWriter.hh b/include/RHEJ/CombinedEventWriter.hh
index 4a4f263..163d133 100644
--- a/include/RHEJ/CombinedEventWriter.hh
+++ b/include/RHEJ/CombinedEventWriter.hh
@@ -1,50 +1,50 @@
/** \file CombinedEventWriter.hh
* \brief Details the CombinedEventWriter class
*
* CombinedEventWriter Class which handles all of the EventWriters being used by a single use
* of RHEJ. It calls all other EventWriters in its list to use their write function and output.
*/
#pragma once
#include <memory>
#include <vector>
#include "RHEJ/EventWriter.hh"
#include "RHEJ/make_writer.hh"
namespace LHEF{
/** \struct HEPRUP CombinedEventWriter.hh "include/RHEJ/CombinedEventWriter.hh"
* \brief HEPRUP struct which contains input event information.
*
* This struct contains important information from the LHEF header such as the beam energy.
*/
struct HEPRUP;
}
namespace RHEJ{
/** \class CombinedEventWriter CombinedEventWriter.hh "include/RHEJ/CombinedEventWriter.hh"
* \brief Intended as an EventWriter to Multiple Output types.
*
* Inherits from EventWriter, and then uses its write function which calls all of the event
* writers write function to output to their specific output files in their specific format.
* This handles this complexity of having multiple EventWriters.
*/
class CombinedEventWriter: public EventWriter{
public:
//! CombinedEventWriter Constructor
CombinedEventWriter(
std::vector<OutputFile> const & outfiles, /**< Vector of Output file Types */
LHEF::HEPRUP const & heprup /**< Contains important information from the Input */
);
- void write(ClusteredEvent const &) override;
+ void write(Event const &) override;
private:
//! A vector of the EventWriters
std::vector<std::unique_ptr<EventWriter>> writers_;
};
}
diff --git a/include/RHEJ/EmptyAnalysis.hh b/include/RHEJ/EmptyAnalysis.hh
index 1da3222..1733892 100644
--- a/include/RHEJ/EmptyAnalysis.hh
+++ b/include/RHEJ/EmptyAnalysis.hh
@@ -1,31 +1,31 @@
/** \file EmptyAnalysis.hh
* \brief Details the EmptyAnalysis Class. (The trivial Analysis case.)
*/
#pragma once
#include <memory>
#include "RHEJ/Analysis.hh"
//! YAML Namespace
namespace YAML{
class Node;
}
namespace RHEJ{
/** \struct EmptyAnalysis EmptyAnalysis.hh "include/RHEJ/EmptyAnalysis.hh"
* \brief EmptyAnalysis struct. Trivial Analysis case.
*
* Inherits from the Analysis Struct and defines its own fill function (as necessary)
* This function does nothing (which is the trivial case.
*/
struct EmptyAnalysis: Analysis{
static std::unique_ptr<Analysis> create(YAML::Node const & parameters);
- virtual void fill(ClusteredEvent const &) override;
+ virtual void fill(Event const &) override;
virtual ~EmptyAnalysis() override = default;
};
}
diff --git a/include/RHEJ/Event.hh b/include/RHEJ/Event.hh
index 91b031b..1a048be 100644
--- a/include/RHEJ/Event.hh
+++ b/include/RHEJ/Event.hh
@@ -1,122 +1,123 @@
/** \file Event.hh
* \brief Details the parameters of an Event.
*
* Contains the EventParameters struct and the Event struct. Also
* contains the ClusteredEvent class.
*/
#pragma once
#include <string>
#include "utility.hh"
#include "LHEF/LHEF.h"
#include "fastjet/JetDefinition.hh"
#include "fastjet/ClusterSequence.hh"
namespace RHEJ{
/** \struct EventParameters Event.hh "include/RHEJ/Event.hh"
* \brief A struct containing the parameters of an event
*/
struct EventParameters{
double mur; /**< Value of the Renormalisation Scale */
double muf; /**< Value of the Factorisation Scale */
double weight; /**< Event Weight */
};
/** \struct UnclusteredEvent Event.hh "include/RHEJ/Event.hh"
* \brief Event Struct which contains Particle Content.
*/
struct UnclusteredEvent{
//! Default Constructor
UnclusteredEvent() = default;
//! Constructor from LesHouches event information
UnclusteredEvent(LHEF::HEPEUP const & hepeup);
std::array<Sparticle, 2> incoming; /**< Incoming Particles */
std::vector<Sparticle> outgoing; /**< Outgoing Particles */
//! Central parameter (e.g. scale) choice
EventParameters central;
std::vector<EventParameters> variations; /**< For parameter variation */
};
- /** \class ClusteredEvent Event.hh "include/RHEJ/Event.hh"
- * \brief Class for ClusteredEvents
+ /** \class Event Event.hh "include/RHEJ/Event.hh"
+ * \brief Class for Events
*
- * This class contains information about a clustered events jets.
+ * This is the default reversed HEJ event class.
+ * It contains kinematic information including jet clustering,
+ * parameter (e.g. scale) settings and the event weight.
*/
- class ClusteredEvent{
+ class Event{
public:
- //! Default ClusteredEvent Constructor
- ClusteredEvent() = default;
-
- //! Specified ClusteredEvent Constructor
- ClusteredEvent(
+ //! 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
);
std::vector<fastjet::PseudoJet> jets() const;
UnclusteredEvent const & unclustered() const {
return ev_;
}
EventParameters const & central() const{
return ev_.central;
}
EventParameters & central(){
return ev_.central;
}
std::array<Sparticle, 2> const & incoming() const{
return ev_.incoming;
}
std::vector<Sparticle> const & outgoing() const{
return ev_.outgoing;
}
std::vector<EventParameters> const & variations() const{
return ev_.variations;
}
std::vector<EventParameters> & variations(){
return ev_.variations;
}
EventParameters const & variations(size_t i) const{
return ev_.variations[i];
}
EventParameters & variations(size_t i){
return ev_.variations[i];
}
std::vector<int> particle_jet_indices(
std::vector<fastjet::PseudoJet> const & jets
) const{
return cs_.particle_jet_indices(jets);
}
fastjet::JetDefinition const & jet_def() const{
return cs_.jet_def();
}
double min_jet_pt() const{
return min_jet_pt_;
}
private:
UnclusteredEvent ev_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
};
- double shat(ClusteredEvent const & ev);
+ double shat(Event const & ev);
- LHEF::HEPEUP to_HEPEUP(ClusteredEvent const & event, LHEF::HEPRUP *);
+ LHEF::HEPEUP to_HEPEUP(Event const & event, LHEF::HEPRUP *);
}
diff --git a/include/RHEJ/EventReweighter.hh b/include/RHEJ/EventReweighter.hh
index 105f68a..4de8cf2 100644
--- a/include/RHEJ/EventReweighter.hh
+++ b/include/RHEJ/EventReweighter.hh
@@ -1,193 +1,192 @@
/** \file EventReweighter.hh
* \brief Main class for reweighting events according to the reversed HEJ method
*
* EventReweighter Class is the main class used within RHEJ. It reweights the
* resummation events.
*/
#pragma once
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/PDF.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/event_classes.hh"
#include "RHEJ/config.hh"
#include "RHEJ/PhaseSpacePoint.hh"
#include "RHEJ/MatrixElement.hh"
namespace RHEJ{
class ScaleGenerator;
/** \struct Beam EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief Beam Struct. Contains beam energy and incoming particle content.
*
* Contains the value of the energy of the beam and the identity of the two incoming Particles
*/
struct Beam{
double E;
std::array<ParticleID, 2> type;
};
/** \class EventReweighter EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief EventReweighter Main class for reweighting events in RHEJ.
*
* This is the class which reweights the resummation phase space points once they have been
* generated by PhaseSpacePoint. This class also handles how many phase space points are generated
* for every fixed order event which is processed.
*/
class EventReweighter{
using EventClass = event_class::EventClass;
public:
EventReweighter(
Beam beam, /**< Beam Energy */
int pdf_id, /**< PDF ID */
fastjet::JetDefinition jet_def, /**< Jet Definition */
double jetptmin, /**< Jet Pt minimum */
EventTreatMap treat, /**< Treat??? */
double extpartonptmin, /**< External Parton Pt Minimum */
double max_ext_soft_pt_fraction, /**< External Parton Pt Maximum */
bool log_corr /**< Log Correct On or Off */
);
EventReweighter(
LHEF::HEPRUP const & heprup, /**< HEPRUP???? */
fastjet::JetDefinition jet_def, /**< Jet Definition */
double jetptmin, /**< Jet Pt minimum */
EventTreatMap treat, /**< Treat??? */
double extpartonptmin, /**< External Parton Pt Minimum */
double max_ext_soft_pt_fraction, /**< External Parton Pt Maximum */
bool log_corr /**< Log Correct On or Off */
);
PDF const & pdf() const; /**< PDF Used */
- std::vector<ClusteredEvent> reweight(
- ClusteredEvent const & ev, /**< Clustered Event For Reweighing */
+ std::vector<Event> reweight(
+ Event const & ev, /**< Event For Reweighing */
int num_events /**< Number of Events */
);
- std::vector<ClusteredEvent> reweight(
- ClusteredEvent const & ev, /**< Clustered Event For Reweighing */
+ std::vector<Event> reweight(
+ Event const & ev, /**< Event For Reweighing */
int num_events, /**< Number of Events */
ScaleGenerator const & gen_scales /**< Generated Scales see config.hh */
);
//! Constructor for last_event_class
EventClass last_event_class() const;
private:
/** \struct EventFactors EventReweighter.hh "include/RHEJ/EventReweighter.hh"
* \brief EventFactors Struct containing central value and variations
*
* Contains the Central Value and Variation around
*/
struct EventFactors{
double central; /**< Central Value for the Factor */
std::vector<double> variations; /**< Vector of Variations from central value */
};
template<typename... T>
PDF const & pdf(T&& ...);
- std::vector<ClusteredEvent> gen_res_events(
- ClusteredEvent const & ev, int num_events
+ std::vector<Event> gen_res_events(
+ Event const & ev, int num_events
);
- std::vector<ClusteredEvent> rescale(
- ClusteredEvent const & Born_ev, std::vector<ClusteredEvent> events
+ std::vector<Event> rescale(
+ Event const & Born_ev, std::vector<Event> events
) const;
-
/**
* \brief Do the Jets pass the resummation Cuts?
*
- * @param ev ClusteredEvent in Question
+ * @param ev Event in Question
* @returns 0 or 1 depending on if ev passes Jet Cuts
*/
- bool jets_pass_resummation_cuts(ClusteredEvent const & ev) const;
+ bool jets_pass_resummation_cuts(Event const & ev) const;
/**
* \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.
*/
- EventFactors pdf_factors(ClusteredEvent const & ev) const;
+ EventFactors pdf_factors(Event const & ev) const;
/**
* \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.
*/
- EventFactors matrix_elements(ClusteredEvent const & ev) const;
+ EventFactors matrix_elements(Event const & ev) const;
/**
* \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.
*/
- EventFactors fixed_order_scale_ME(ClusteredEvent const & ev) const;
+ EventFactors fixed_order_scale_ME(Event const & ev) const;
/**
* \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(ClusteredEvent const & ev) const;
+ double tree_matrix_element(Event const & ev) const;
/**
* \brief Generation of Resummation Phase Space Points
*
* @param ev Event In Question
* @returns Resummation PhaseSpacePoint
*
* Where is this even called? I can't find it with a grep -r in top directory???
*/
PhaseSpacePoint gen_PhaseSpacePoint(Event const & ev) const;
EventClass last_event_class_; /**< Class of the last event */
double extpartonptmin_; /**< External Parton Minimum Pt */
double max_ext_soft_pt_fraction_; /**< External Parton Maximum Pt fraction */
double jetptmin_; /**< Jet Pt Minimum threshold */
double E_beam_; /**< Energy of Beam */
fastjet::JetDefinition jet_def_; /**< Jet Definition being used */
PDF pdf_; /**< Relevant PDF */
MatrixElement MEt2_; /**< Matrix Element Object */
EventTreatMap treat_;
};
template<typename... T>
PDF const & EventReweighter::pdf(T&&... t){
return pdf_ = PDF{std::forward<T>(t)...};
}
}
diff --git a/include/RHEJ/EventWriter.hh b/include/RHEJ/EventWriter.hh
index 71103fd..c121049 100644
--- a/include/RHEJ/EventWriter.hh
+++ b/include/RHEJ/EventWriter.hh
@@ -1,24 +1,24 @@
/** \file EventWriter.hh
* \brief Contains the EventWriter Class which is inherited by all writers.
*/
#pragma once
#include <string>
#include <exception>
namespace RHEJ{
- class ClusteredEvent;
+ class Event;
/** \struct EventWriter EventWriter.hh "include/RHEJ/EventWriter.hh"
* \brief The basis of every other writer in RHEJ
*
* All other writers inherit from this, regardless of final
* output. Includes the virtual function write and a destructor
*/
struct EventWriter{
- virtual void write(ClusteredEvent const &) = 0;
+ virtual void write(Event const &) = 0;
virtual ~EventWriter() = default;
};
}
diff --git a/include/RHEJ/HepMCWriter.hh b/include/RHEJ/HepMCWriter.hh
index 04f1108..042056b 100644
--- a/include/RHEJ/HepMCWriter.hh
+++ b/include/RHEJ/HepMCWriter.hh
@@ -1,43 +1,43 @@
/** \file HepMCWriter.hh
* \brief Contains the EventWriter which is necessary for HepMC Output.
*/
#pragma once
#ifdef RHEJ_BUILD_WITH_HepMC
#include "HepMC/WriterAscii.h"
#endif
#include <string>
#include "RHEJ/EventWriter.hh"
#include "LHEF/LHEF.h"
namespace RHEJ{
class Event;
/** \class HepMCWriter HepMCWriter.hh "include/RHEJ/HepMCWriter.hh"
* \brief This is an event writer specifically for HepMC output.
* as such it inherits everything from the EventWriter class
* and also includes HepMCWriter constructors
*/
class HepMCWriter: public EventWriter{
public:
HepMCWriter(std::string const & file, LHEF::HEPRUP heprup);
HepMCWriter & operator=(HepMCWriter const & other) = delete;
HepMCWriter(HepMCWriter const & other) = delete;
HepMCWriter(HepMCWriter && other) = delete;
HepMCWriter & operator=(HepMCWriter && other) = delete;
~HepMCWriter() override;
- void write(ClusteredEvent const & ev) override;
+ void write(Event const & ev) override;
#ifdef RHEJ_BUILD_WITH_HepMC
private:
HepMC::shared_ptr<HepMC::GenRunInfo> runinfo_;
HepMC::WriterAscii writer_;
#endif
};
}
diff --git a/include/RHEJ/LesHouchesWriter.hh b/include/RHEJ/LesHouchesWriter.hh
index 53e21a7..a312644 100644
--- a/include/RHEJ/LesHouchesWriter.hh
+++ b/include/RHEJ/LesHouchesWriter.hh
@@ -1,55 +1,55 @@
/** \file LesHouchesWriter.hh
* \brief Contains the writer for LesHouches Output
*/
#pragma once
#include <fstream>
#include "RHEJ/EventWriter.hh"
#include "LHEF/LHEF.h"
namespace RHEJ{
class Event;
/** \class LesHouchesWriter LesHouchesWriter.hh "include/RHEJ/LesHouchesWriter.hh"
* \brief EventWriter Class specifically for LesHouches Output
*
* This is an event writer specifically for Les Houches Output. This
* inherits from the EventWriter Class and contains its own
* contructors. TODO: in principle, this class should be movable but that
* somehow(?) breaks the write member function.
**/
class LesHouchesWriter : public EventWriter{
public:
LesHouchesWriter(std::string const & file, LHEF::HEPRUP heprup);
LesHouchesWriter(LesHouchesWriter const & other) = delete;
LesHouchesWriter & operator=(LesHouchesWriter const & other) = delete;
// TODO: in principle, this class should be movable
// but that somehow(?) breaks the write member function
LesHouchesWriter(LesHouchesWriter && other) = delete;
LesHouchesWriter & operator=(LesHouchesWriter && other) = delete;
~LesHouchesWriter() override;
- void write(ClusteredEvent const & ev) override;
+ void write(Event const & ev) override;
private:
void write_init(){
writer_->init();
}
void rewrite_init();
LHEF::HEPRUP & heprup(){
return writer_->heprup;
}
LHEF::HEPEUP & hepeup(){
return writer_->hepeup;
}
std::fstream out_;
std::unique_ptr<LHEF::Writer> writer_;
};
}
diff --git a/include/RHEJ/PhaseSpacePoint.hh b/include/RHEJ/PhaseSpacePoint.hh
index dfb075e..261a52d 100644
--- a/include/RHEJ/PhaseSpacePoint.hh
+++ b/include/RHEJ/PhaseSpacePoint.hh
@@ -1,156 +1,156 @@
/** \file PhaseSpacePoint.hh
* \brief Contains the PhaseSpacePoint Class
*/
#pragma once
#include <vector>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RanluxEngine.h>
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/Splitter.hh"
namespace RHEJ{
//! A point in resummation phase space
class PhaseSpacePoint{
public:
//! Default PhaseSpacePoint Constructor
PhaseSpacePoint() = default;
//! PhaseSpacePoint Constructor
/**
* @param ev Clustered Jet Event
* @param jet_def The Jet Definition Used
* @param jetptmin Minimum Jet Transverse Momentum
* @param extpartonptmin Minimum Parton Transverse Momentum
* @param max_ext_soft_pt Maximum Parton Transverse Momentum?
*/
PhaseSpacePoint(
- ClusteredEvent const & ev,
+ Event const & ev,
fastjet::JetDefinition jet_def, double jetptmin,
double extpartonptmin, double max_ext_soft_pt
);
//! Get Weight Function
/**
* @returns Weight of Event
*/
double weight() const{
return weight_;
}
//! Get Incoming Function
/**
* @returns Incoming Particles
*/
std::array<Sparticle, 2> const & incoming() const{
return incoming_;
}
//! Get Outgoing Function
/**
* @returns Outgoing Particles
*/
std::vector<Sparticle> const & outgoing() const{
return outgoing_;
}
static constexpr double ccut = 0.2; // universal min pt cut
static constexpr int ng_max = 1000; // maximum number of extra gluons
static void reset_ranlux(std::string const & init_file);
static void reset_ranlux(char const * init_file);
private:
std::vector<fastjet::PseudoJet> cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const;
bool pass_resummation_cuts(
std::vector<fastjet::PseudoJet> const & jets
) const;
bool pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const;
int sample_ng(std::vector<fastjet::PseudoJet> const & Born_jets);
int sample_ng_jets(int ng, std::vector<fastjet::PseudoJet> const & Born_jets);
double probability_in_jet(
std::vector<fastjet::PseudoJet> const & Born_jets
) const;
std::vector<fastjet::PseudoJet> gen_non_jet(
int ng_non_jet,
double ptmin, double ptmax
);
void rescale_rapidities(
std::vector<fastjet::PseudoJet> & partons,
double ymin, double ymax
);
std::vector<fastjet::PseudoJet> reshuffle(
std::vector<fastjet::PseudoJet> const & Born_jets,
fastjet::PseudoJet const & q
);
bool jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const;
void reconstruct_incoming(std::array<Sparticle, 2> const & Born_incoming);
double phase_space_normalisation(
int num_Born_jets,
int num_res_partons
) const;
std::vector<fastjet::PseudoJet> split(
std::vector<fastjet::PseudoJet> const & jets, int ng_jets
);
std::vector<int> distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets
);
std::vector<fastjet::PseudoJet> split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np_in_jet
);
bool split_preserved_jets(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<fastjet::PseudoJet> const & jet_partons
) const;
template<class Particle>
Particle const & most_backward_FKL(
std::vector<Particle> const & partons
) const;
template<class Particle>
Particle const & most_forward_FKL(
std::vector<Particle> const & partons
) const;
template<class Particle>
Particle & most_backward_FKL(std::vector<Particle> & partons) const;
template<class Particle>
Particle & most_forward_FKL(std::vector<Particle> & partons) const;
bool extremal_FKL_ok(
std::vector<fastjet::PseudoJet> const & partons
) const;
bool momentum_conserved(double ep) const;
bool unob_, unof_;
double weight_;
double extpartonptmin_;
double max_ext_soft_pt_fraction_;
double jetptmin_;
fastjet::JetDefinition jet_def_;
std::array<Sparticle, 2> incoming_;
std::vector<Sparticle> outgoing_;
static CLHEP::Ranlux64Engine ran_;
HejSplit splitter_;
};
}
diff --git a/include/RHEJ/ScaleFunction.hh b/include/RHEJ/ScaleFunction.hh
index 6a7ce0f..c2a598a 100644
--- a/include/RHEJ/ScaleFunction.hh
+++ b/include/RHEJ/ScaleFunction.hh
@@ -1,99 +1,99 @@
/** \file ScaleFunction.hh
* \brief ScaleFunction Header file. Contains the different scale choice structs.
*/
#pragma once
#include <memory>
#include "fastjet/JetDefinition.hh"
#include "RHEJ/Event.hh"
namespace RHEJ{
/** \struct ScaleFunction ScaleFunction.hh "include/RHEJ/ScaleFunction.hh"
* \brief Contains EventParameters and ClusteredEvent struct
*/
struct ScaleFunction{
- virtual EventParameters operator()(ClusteredEvent const & ev) const = 0;
+ virtual EventParameters operator()(Event const & ev) const = 0;
virtual ~ScaleFunction(){};
};
/** \struct InputScales ScaleFunction.hh "include/RHEJ/ScaleFunction.hh"
* \brief InputScales Struct which inherits from ScaleFunction Struct.
*
* This is the struct used when the user inputs a scale as the scale
* choice used in calculation
*/
struct InputScales: ScaleFunction{
- EventParameters operator()(ClusteredEvent const & ev) const override{
+ EventParameters operator()(Event const & ev) const override{
return ev.central();
}
};
/** \class FixedScale ScaleFunction.hh "include/RHEJ/ScaleFunction.hh"
* \brief FixedScale Class which inherits from ScaleFuncton
*
* This is the class used when the user selects fixed scale choice
* in the input card.
*/
class FixedScale: public ScaleFunction{
public:
FixedScale(double scale): scale_{scale} {}
- EventParameters operator()(ClusteredEvent const & ev) const override{
+ EventParameters operator()(Event const & ev) const override{
return {scale_, scale_, ev.central().weight};
}
private:
double scale_;
};
/** \struct Ht ScaleFunction.hh "include/RHEJ/ScaleFunction.hh"
* \brief Ht Struct which inherits from ScaleFunction
*
* This is the struct used when the user selects Ht as the scale choice
* in the input card.
*/
struct Ht: public ScaleFunction{
public:
- EventParameters operator()(ClusteredEvent const & ev) const override;
+ EventParameters operator()(Event const & ev) const override;
};
/** \struct MaxJetPperp ScaleFunction.hh "include/RHEJ/ScaleFunction.hh"
* \brief MaxJetPperp Struct which inherits from ScaleFunction.
*
* This is the struct used when the user selects MaxJetPperp as the scale choice
* in the input card.
*/
struct MaxJetPperp: public ScaleFunction{
- EventParameters operator()(ClusteredEvent const & ev) const override;
+ EventParameters operator()(Event const & ev) const override;
};
/** \struct JetInvariantMass ScaleFunction.hh "include/RHEJ/ScaleFunction.hh"
* \brief JetInvariantMass Struct which inherits from ScaleFunction.
*
* This is the struct used when the user selects JetInvarientMass as the scale choice
* in the input card.
*/
struct JetInvariantMass: public ScaleFunction{
- EventParameters operator()(ClusteredEvent const & ev) const override;
+ EventParameters operator()(Event const & ev) const override;
};
/** \class Product ScaleFunction.hh "include/RHEJ/ScaleFunction.hh"
* \brief Product Class which inherits from ScaleFunction.
*
* This is the class used when the user selects Product as the scale choice
* in the input card. THIS DOESN'T SEEM RIGHT.
*/
class Product: public ScaleFunction{
public:
Product(double factor, std::unique_ptr<ScaleFunction> base_scale):
factor_{factor},
base_scale_{std::move(base_scale)}
{}
- EventParameters operator()(ClusteredEvent const & ev) const override;
+ EventParameters operator()(Event const & ev) const override;
private:
double factor_;
std::unique_ptr<ScaleFunction> base_scale_;
};
}
diff --git a/include/RHEJ/config.hh b/include/RHEJ/config.hh
index f6296ea..bf6a3ba 100644
--- a/include/RHEJ/config.hh
+++ b/include/RHEJ/config.hh
@@ -1,106 +1,106 @@
/** \file config.hh
* \brief The file which handles the configuration file parameters
*
* Contains the JetParameters Struct and ScaleConfig Struct. Also
* contains the ScaleGenerator Class and the EventTreatment class.
* Config struct is also defined here. The configuration files
* parameters are read and then stored within this objects.
*/
#pragma once
#include <string>
#include <memory>
#include "fastjet/JetDefinition.hh"
#include "RHEJ/ScaleFunction.hh"
#include "RHEJ/optional.hh"
#include "RHEJ/output_formats.hh"
#include "RHEJ/event_classes.hh"
#include "yaml-cpp/yaml.h"
namespace RHEJ{
/**! \struct JetParameters config.hh "include/RHEJ/config.hh"
* \brief Contains Jet Definition and Minimum Pt to be a Jet.
*
* Contains the Fastjet definition of a Jet and a threshold (minimum) value for pt
*/
struct JetParameters{
fastjet::JetDefinition def; /**< Jet Definition */
double min_pt; /**< Minimum Jet Pt */
};
/**! \struct ScaleConfig config.hh "include/RHEJ/config.hh"
* \brief Sets up the possible choices for scale variation?
*
* There is a vector which contains the possible scale choices and
* also vector of doubles with the scale factors along with a max
* scale ratio.
*/
struct ScaleConfig{
std::vector<std::unique_ptr<ScaleFunction>> scales; /**< Vector of possible Scale choices */
std::vector<double> scale_factors; /**< Vector of possible Scale Factors */
double max_scale_ratio; /**< Maximum ratio for the scales */
};
/**! \class ScaleGenerator config.hh "include/RHEJ/config.hh"
* \brief A Class used to generate scales
*
* This class has a clustered event class and scale config struct
* within it.
*/
class ScaleGenerator{
public:
ScaleGenerator() = default;
explicit ScaleGenerator(ScaleConfig && settings):
settings_{std::move(settings)}
{}
- ClusteredEvent operator()(ClusteredEvent ev) const;
+ Event operator()(Event ev) const;
ScaleConfig const & settings() const;
private:
ScaleConfig settings_;
};
/**
* \brief Enumeration which deals with how to treat events.
*
* The program will decide on how to treat an event based on
* the value of this enumeration.
*/
enum class EventTreatment{
reweight, /**< Reweigh the event */
keep, /**< Keep the event */
discard, /**< Discard the event */
};
using EventTreatMap = std::map<event_class::EventClass, EventTreatment>;
/**! \struct Config config.hh "include/RHEJ/config.hh"
* \brief Config Struct for user input parameters.
*
* The struct which handles the input card given by the user.
*/
struct Config {
ScaleGenerator scale_gen; /**< Scale */
JetParameters resummation_jets; /**< Resummation Jets */
JetParameters fixed_order_jets; /**< Fixed-Order Jets */
double min_extparton_pt; /**< Min External Parton Pt*/
double max_ext_soft_pt_fraction; /**< Min External Parton Pt Fraction */
int trials; /**< Number of resummation points to generate per FO */
bool log_correction; /**< Log Correct On or Off */
bool unweight; /**< Unweight On or Off */
std::vector<OutputFile> output; /**< Output File Type */
optional<std::string> RanLux_init; /**< Ranlux File Choice*/
EventTreatMap treat; /**< Event Treat Map? */
YAML::Node analysis_parameters; /**< Analysis Parameters */
};
Config load_config(std::string const & config_file);
}
diff --git a/include/RHEJ/event_classes.hh b/include/RHEJ/event_classes.hh
index 8d6ae57..30adb8d 100644
--- a/include/RHEJ/event_classes.hh
+++ b/include/RHEJ/event_classes.hh
@@ -1,49 +1,49 @@
/** \file event_classes.hh
* \brief This file details the classification of events.
*
* This file makes use of a macro in order to avoid repetition of EventClass names.
* To this end, the Documentation of the EventClass enumeration below is incomplete.
* The possible event classifications are stored within that enumeration.
*/
#pragma once
#include "RHEJ/utility.hh"
namespace RHEJ{
// macro definition to avoid repetition of EventClass names
#define RHEJ_EVENT_CLASSES(F) F(FKL), F(unordered_backward), F(unordered_forward), F(nonFKL), F(no_2_jets), F(bad_final_state)
#define RHEJ_AS_ENUM(VAR) VAR
#define RHEJ_AS_STRING(VAR) #VAR
//! Event_Class NameSpace
namespace event_class{
/** \enum EventClass
* \brief EventClass enumeration gives different possible event classes
- *
+ *
* This enumeration is used to distinguish between different event types.
*/
enum EventClass: size_t{
RHEJ_EVENT_CLASSES(RHEJ_AS_ENUM), /**< Macro of Possible States */
unob = unordered_backward, /**< Unordered Emission by backwards Jet */
unof = unordered_forward, /**< Unordered Emission by forwards Jet */
first_class = FKL, /**< FKL Event */
last_class = bad_final_state /**< Bad Final State Event */
};
-
+
static constexpr auto names = make_array(
RHEJ_EVENT_CLASSES(RHEJ_AS_STRING)
);
}
#undef RHEJ_HISTOGRAMS
#undef RHEJ_AS_ENUM
#undef RHEJ_AS_STRING
- class ClusteredEvent;
+ class Event;
//! A Function used to classify a clustered Event
- event_class::EventClass classify(ClusteredEvent const & ev);
+ event_class::EventClass classify(Event const & ev);
}
diff --git a/src/CombinedEventWriter.cc b/src/CombinedEventWriter.cc
index f7a852b..5b1cf38 100644
--- a/src/CombinedEventWriter.cc
+++ b/src/CombinedEventWriter.cc
@@ -1,21 +1,21 @@
#include "RHEJ/CombinedEventWriter.hh"
namespace RHEJ{
CombinedEventWriter::CombinedEventWriter(
std::vector<OutputFile> const & outfiles,
LHEF::HEPRUP const & heprup
){
writers_.reserve(outfiles.size());
for(OutputFile const & outfile: outfiles){
writers_.emplace_back(
make_format_writer(outfile.format, outfile.name, heprup)
);
}
}
- void CombinedEventWriter::write(ClusteredEvent const & ev){
+ void CombinedEventWriter::write(Event const & ev){
for(auto & writer: writers_) writer->write(ev);
}
}
diff --git a/src/EmptyAnalysis.cc b/src/EmptyAnalysis.cc
index 63963a8..76067f5 100644
--- a/src/EmptyAnalysis.cc
+++ b/src/EmptyAnalysis.cc
@@ -1,42 +1,42 @@
#include "RHEJ/EmptyAnalysis.hh"
#include <iostream>
#include <string>
#include "yaml-cpp/yaml.h"
namespace RHEJ{
namespace{
void warn_ignore(std::string const & what){
std::cerr << "WARNING: ignoring unknown analysis parameter "
<< what << '\n';
}
void warn_unless_empty(YAML::Node const & parameters){
using YAML::NodeType;
switch(parameters.Type()){
case NodeType::Null:
case NodeType::Undefined:
return;
case NodeType::Scalar:
return warn_ignore(parameters.as<std::string>());
case NodeType::Sequence:
for(auto && param: parameters) warn_ignore(param.as<std::string>());
return;
case NodeType::Map:
for(auto && param: parameters) warn_ignore(param.first.as<std::string>());
return;
}
}
}
std::unique_ptr<Analysis> EmptyAnalysis::create(
YAML::Node const & parameters
){
warn_unless_empty(parameters);
return std::unique_ptr<Analysis>{new EmptyAnalysis{}};
}
- void EmptyAnalysis::fill(ClusteredEvent const &){
+ void EmptyAnalysis::fill(Event const &){
// do nothing
}
}
diff --git a/src/Event.cc b/src/Event.cc
index 0749a72..40f9a13 100644
--- a/src/Event.cc
+++ b/src/Event.cc
@@ -1,137 +1,137 @@
#include "RHEJ/Event.hh"
#include "RHEJ/debug.hh"
namespace RHEJ{
namespace{
constexpr int status_in = -1;
constexpr int status_out = 1;
}
UnclusteredEvent::UnclusteredEvent(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) {
// Only consider final state particles
if (abs(hepeup.ISTUP[i]) != status_out) continue;
// Save particle information
Sparticle 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]
}
};
if (hepeup.ISTUP[i]>0) outgoing.emplace_back(std::move(particle));
else{
if(in_idx > incoming.size()) {
throw std::invalid_argument{
"Event has too many incoming particles"
};
}
incoming[in_idx++] = std::move(particle);
}
}
std::sort(
begin(incoming), end(incoming),
[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
);
std::sort(begin(outgoing), end(outgoing), rapidity_less{});
}
- ClusteredEvent::ClusteredEvent(
+ Event::Event(
UnclusteredEvent ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
):
ev_{std::move(ev)},
cs_{to_PseudoJet(filter_partons(ev_.outgoing)), jet_def},
min_jet_pt_{min_jet_pt}
{}
- std::vector<fastjet::PseudoJet> ClusteredEvent::jets() const{
+ std::vector<fastjet::PseudoJet> Event::jets() const{
return cs_.inclusive_jets(min_jet_pt_);
}
- double shat(ClusteredEvent const & ev){
+ 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<Sparticle, 2> const & incoming,
std::vector<Sparticle> 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(ClusteredEvent const & event, LHEF::HEPRUP * heprup){
+ 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);
}
const size_t num_particles =
event.incoming().size() + event.outgoing().size();
result.NUP = num_particles;
// the following entries are pretty much meaningless
result.IDPRUP = 1;
result.AQEDUP = 1./128.;
// 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(Sparticle 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(Sparticle const & out: event.outgoing()){
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()});
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()),
[](Sparticle 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)
);
}
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 2180589..a02567d 100644
--- a/src/EventReweighter.cc
+++ b/src/EventReweighter.cc
@@ -1,316 +1,316 @@
#include "RHEJ/EventReweighter.hh"
#include <string>
#include <unordered_map>
#include "RHEJ/PhaseSpacePoint.hh"
#include "RHEJ/PDG_codes.hh"
#include "RHEJ/config.hh" // ScaleGenerator definition
#include "RHEJ/debug.hh"
namespace RHEJ{
using EventClass = event_class::EventClass;
namespace {
static_assert(
std::numeric_limits<double>::has_quiet_NaN,
"no quiet NaN for double"
);
constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
UnclusteredEvent to_UnclusteredEvent(PhaseSpacePoint const & psp){
UnclusteredEvent result;
result.incoming = psp.incoming();
std::sort(
begin(result.incoming), end(result.incoming),
[](Sparticle o1, Sparticle o2){return o1.p.pz()<o2.p.pz();}
);
assert(result.incoming.size() == 2);
result.outgoing = psp.outgoing();
assert(
std::is_sorted(
begin(result.outgoing), end(result.outgoing),
rapidity_less{}
)
);
assert(result.outgoing.size() >= 2);
result.central.mur = NaN;
result.central.muf = NaN;
result.central.weight = psp.weight();
return result;
}
} // anonymous namespace
EventReweighter::EventReweighter(
LHEF::HEPRUP const & heprup,
fastjet::JetDefinition jet_def, double jetptmin,
EventTreatMap treat,
double extpartonptmin, double max_ext_soft_pt_fraction,
bool log_corr
):
EventReweighter{
RHEJ::Beam{
heprup.EBMUP.first,
{{
static_cast<RHEJ::ParticleID>(heprup.IDBMUP.first),
static_cast<RHEJ::ParticleID>(heprup.IDBMUP.second)
}}
},
heprup.PDFSUP.first,
jet_def, jetptmin,
std::move(treat),
extpartonptmin, max_ext_soft_pt_fraction,
log_corr
}
{
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,
fastjet::JetDefinition jet_def, double jetptmin,
EventTreatMap treat,
double extpartonptmin, double max_ext_soft_pt_fraction,
bool log_corr
):
extpartonptmin_{extpartonptmin},
max_ext_soft_pt_fraction_{max_ext_soft_pt_fraction},
jetptmin_{jetptmin},
E_beam_{beam.E},
jet_def_{jet_def},
pdf_{pdf_id, beam.type.front(), beam.type.back()},
MEt2_{jet_def, jetptmin, log_corr},
treat_{std::move(treat)}
{}
PDF const & EventReweighter::pdf() const{
return pdf_;
}
EventClass EventReweighter::last_event_class() const{
return last_event_class_;
}
namespace{
static_assert(
std::numeric_limits<double>::has_infinity, "infinity not supported"
);
// there is no simple way to create a vector of move-only objects
// this is a clumsy work-around
ScaleGenerator make_identity(){
std::vector<std::unique_ptr<ScaleFunction>> scale_fun;
scale_fun.emplace_back(new InputScales);
return ScaleGenerator{
ScaleConfig{
std::move(scale_fun), {}, std::numeric_limits<double>::infinity()
}
};
}
const ScaleGenerator identity{make_identity()};
}
- std::vector<ClusteredEvent> EventReweighter::reweight(
- ClusteredEvent const & input_ev, int num_events
+ std::vector<Event> EventReweighter::reweight(
+ Event const & input_ev, int num_events
){
return reweight(input_ev, num_events, identity);
}
- std::vector<ClusteredEvent> EventReweighter::reweight(
- ClusteredEvent const & input_ev, int num_events,
+ std::vector<Event> EventReweighter::reweight(
+ Event const & input_ev, int num_events,
ScaleGenerator const & gen_scales
){
auto res_events = gen_res_events(input_ev, num_events);
if(res_events.empty()) return {};
for(auto & event: res_events) event = gen_scales(event);
return rescale(input_ev, std::move(res_events));
}
// main generation/reweighting function
// generate phase space points and divide out Born factors
- std::vector<ClusteredEvent> EventReweighter::gen_res_events(
- ClusteredEvent const & ev,
+ std::vector<Event> EventReweighter::gen_res_events(
+ Event const & ev,
int phase_space_points
){
assert(ev.variations().empty());
last_event_class_ = classify(ev);
switch(treat_.at(last_event_class_)){
case EventTreatment::discard: return {};
case EventTreatment::keep:
if(! jets_pass_resummation_cuts(ev)) return {};
else return {ev};
default:;
}
const double Born_shat = shat(ev);
- std::vector<ClusteredEvent> resummation_events;
+ std::vector<Event> resummation_events;
for(int psp_number = 0; psp_number < phase_space_points; ++psp_number){
PhaseSpacePoint psp{
ev,
jet_def_, jetptmin_,
extpartonptmin_, max_ext_soft_pt_fraction_
};
if(psp.weight() == 0.) continue;
resummation_events.emplace_back(
to_UnclusteredEvent(std::move(psp)), jet_def_, jetptmin_
);
auto & new_event = resummation_events.back();
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<ClusteredEvent> EventReweighter::rescale(
- ClusteredEvent const & Born_ev,
- std::vector<ClusteredEvent> 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.central().weight *= pdf.central*ME.central/(Born_pdf*Born_ME);
for(size_t i = 0; i < cur_event.variations().size(); ++i){
cur_event.variations(i).weight *=
pdf.variations[i]*ME.variations[i]/(Born_pdf*Born_ME);
}
}
return events;
};
bool EventReweighter::jets_pass_resummation_cuts(
- ClusteredEvent const & ev
+ Event const & ev
) const{
const auto out_as_PseudoJet = to_PseudoJet(filter_partons(ev.outgoing()));
fastjet::ClusterSequence cs{out_as_PseudoJet, jet_def_};
return cs.inclusive_jets(jetptmin_).size() == ev.jets().size();
}
EventReweighter::EventFactors
- EventReweighter::pdf_factors(ClusteredEvent const & ev) const{
+ 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_;
EventFactors 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 & param: ev.variations()){
const double muf = 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;
}
EventReweighter::EventFactors
- EventReweighter::matrix_elements(ClusteredEvent const & ev) const{
+ EventReweighter::matrix_elements(Event const & ev) const{
assert(treat_.count(last_event_class_) > 0);
if(treat_.find(last_event_class_)->second == EventTreatment::keep){
return fixed_order_scale_ME(ev);
}
// precompute overall kinematic factor
const double ME_kin = MEt2_.tree_kin(ev.incoming(), ev.outgoing(), true);
EventFactors result;
std::unordered_map<double, double> known_ME;
result.central = MEt2_(
pdf_.Halphas(ev.central().mur), ev.central().mur,
ev.incoming(), ev.outgoing(),
true
);
known_ME.emplace(ev.central().mur, result.central);
result.variations.reserve(ev.variations().size());
for(auto const & param: ev.variations()){
const double mur = param.mur;
auto cur_ME = known_ME.find(mur);
if(cur_ME == known_ME.end()){
const double alpha_s = pdf_.Halphas(mur);
const double ME = MEt2_.tree_param(
alpha_s, mur, ev.incoming(), ev.outgoing()
)*ME_kin*MEt2_.virtual_corrections(
alpha_s, mur, ev.incoming(), ev.outgoing()
);
cur_ME = known_ME.emplace(mur, ME).first;
}
result.variations.emplace_back(cur_ME->second);
}
assert(result.variations.size() == ev.variations().size());
return result;
}
- double EventReweighter::tree_matrix_element(ClusteredEvent const & ev) const{
+ double EventReweighter::tree_matrix_element(Event const & ev) const{
assert(ev.variations().empty());
assert(treat_.count(last_event_class_) > 0);
if(treat_.find(last_event_class_)->second == EventTreatment::keep){
return fixed_order_scale_ME(ev).central;
}
return MEt2_.tree(
pdf_.Halphas(ev.central().mur), ev.central().mur,
ev.incoming(), ev.outgoing(),
false
);
}
EventReweighter::EventFactors
- EventReweighter::fixed_order_scale_ME(ClusteredEvent const & ev) const{
+ EventReweighter::fixed_order_scale_ME(Event const & ev) const{
const int alpha_s_power = std::count_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](Sparticle const & p){ return is_parton(p); }
);
EventFactors 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;
}
}
diff --git a/src/HepMCWriter.cc b/src/HepMCWriter.cc
index 8ae2989..25e848f 100644
--- a/src/HepMCWriter.cc
+++ b/src/HepMCWriter.cc
@@ -1,108 +1,108 @@
#include "RHEJ/HepMCWriter.hh"
#include <cassert>
#ifdef RHEJ_BUILD_WITH_HepMC
#include "HepMC/LHEFAttributes.h"
#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"
#include "RHEJ/utility.hh"
#include "RHEJ/Event.hh"
namespace RHEJ{
namespace {
void add_generator_tag(LHEF::HEPRUP & heprup){
heprup.generators.emplace_back(LHEF::XMLTag{});
heprup.generators.back().name = "HEJ";
heprup.generators.back().version = "0.0.1";
}
void reset_weight_info(LHEF::HEPRUP & heprup){
heprup.IDWTUP = -4;
// use placeholders for unknown init block values
// we can overwrite them after processing all events
heprup.XSECUP = {0.};
heprup.XERRUP = {0.};
heprup.XMAXUP = {0.};
}
HepMC::FourVector to_FourVector(Sparticle const & sp){
return {sp.px(), sp.py(), sp.pz(), sp.E()};
}
HepMC::shared_ptr<HepMC::GenRunInfo> & add_HEPRUP(
HepMC::shared_ptr<HepMC::GenRunInfo> & runinfo,
LHEF::HEPRUP heprup
){
add_generator_tag(heprup);
reset_weight_info(heprup);
auto hepr = HepMC::make_shared<HepMC::HEPRUPAttribute>();
hepr->heprup = std::move(heprup);
runinfo->add_attribute("HEPRUP", hepr);
for (int i = 0, N = hepr->heprup.generators.size(); i < N; ++i ){
HepMC::GenRunInfo::ToolInfo tool;
tool.name = hepr->heprup.generators[i].name;
tool.version = hepr->heprup.generators[i].version;
tool.description = hepr->heprup.generators[i].contents;
runinfo->tools().push_back(tool);
}
return runinfo;
}
}
HepMCWriter::HepMCWriter(std::string const & file, LHEF::HEPRUP heprup):
runinfo_{HepMC::make_shared<HepMC::GenRunInfo>()},
writer_{file, add_HEPRUP(runinfo_, std::move(heprup))}
{}
HepMCWriter::~HepMCWriter(){
writer_.close();
}
- void HepMCWriter::write(ClusteredEvent const & ev){
+ void HepMCWriter::write(Event const & ev){
auto vx = HepMC::make_shared<HepMC::GenVertex>();
for(auto const & in: ev.incoming()){
vx->add_particle_in(
HepMC::make_shared<HepMC::GenParticle>(
to_FourVector(in), static_cast<int>(in.type), -1
)
);
}
for(auto const & out: ev.outgoing()){
vx->add_particle_out(
HepMC::make_shared<HepMC::GenParticle>(
to_FourVector(out), static_cast<int>(out.type), +1
)
);
}
HepMC::GenEvent out_ev{runinfo_, HepMC::Units::GEV, HepMC::Units::MM};
out_ev.add_vertex(vx);
out_ev.weights().reserve(ev.variations().size() + 1u);
out_ev.weights().emplace_back(ev.central().weight);
for(auto const & var: ev.variations()){
out_ev.weights().emplace_back(var.weight);
}
writer_.write_event(out_ev);
}
}
#else
namespace RHEJ{
HepMCWriter::HepMCWriter(std::string const &, LHEF::HEPRUP){
throw std::invalid_argument(
"Failed to create HepMC writer: "
"Reversed HEJ was built without HepMC support"
);
}
- void HepMCWriter::write(ClusteredEvent const &){
+ void HepMCWriter::write(Event const &){
assert(false);
}
HepMCWriter::~HepMCWriter() = default;
}
#endif
diff --git a/src/LesHouchesWriter.cc b/src/LesHouchesWriter.cc
index be7198a..2e321d9 100644
--- a/src/LesHouchesWriter.cc
+++ b/src/LesHouchesWriter.cc
@@ -1,82 +1,82 @@
#include <stdexcept>
#include <memory>
#include <cassert>
#include "RHEJ/LesHouchesWriter.hh"
#include "RHEJ/Event.hh"
namespace RHEJ{
namespace{
template<class T, class... Args>
std::unique_ptr<T> make_unique(Args&&... a){
return std::unique_ptr<T>{new T{std::forward<Args>(a)...}};
}
}
LesHouchesWriter::LesHouchesWriter(
std::string const & file, LHEF::HEPRUP heprup
):
out_{file, std::fstream::in | std::fstream::out | std::fstream::trunc},
writer_{RHEJ::make_unique<LHEF::Writer>(out_)}
{
if(! out_.is_open()){
throw std::ios_base::failure("Failed to open " + file);
};
writer_->heprup = std::move(heprup);
writer_->heprup.IDWTUP = -4;
writer_->heprup.generators.emplace_back(LHEF::XMLTag{});
writer_->heprup.generators.back().name = "HEJ";
writer_->heprup.generators.back().version = "0.0.1";
// use placeholders for unknown init block values
// we can overwrite them after processing all events
writer_->heprup.XSECUP = {0.};
writer_->heprup.XERRUP = {0.};
writer_->heprup.XMAXUP = {0.};
write_init();
}
- void LesHouchesWriter::write(ClusteredEvent const & ev){
+ void LesHouchesWriter::write(Event const & ev){
assert(writer_ && out_.is_open());
const double wt = ev.central().weight;
writer_->hepeup = RHEJ::to_HEPEUP(std::move(ev), &heprup());
writer_->writeEvent();
heprup().XSECUP.front() += wt;
heprup().XERRUP.front() += wt*wt;
if(wt > heprup().XMAXUP.front()){
heprup().XMAXUP.front() = wt;
}
}
// this function is called after overwritting the Les Houches init block
// assert that we have overwritten *exactly* the init block,
// i.e. an intact event block is next
void assert_next_event_intact(std::istream & out){
(void) out; // suppress compiler warnings if not in debug mode
#ifndef NDEBUG
std::string line;
getline(out, line);
assert(line == "<event>");
#endif
}
void LesHouchesWriter::rewrite_init(){
assert(writer_ && out_.is_open());
// replace placeholder entries
const auto pos = out_.tellp();
out_.seekp(0);
writer_->init();
assert_next_event_intact(out_);
out_.seekp(pos);
}
LesHouchesWriter::~LesHouchesWriter(){
assert(writer_ && out_.is_open());
heprup().XERRUP.front() = sqrt(heprup().XERRUP.front());
rewrite_init();
}
}
diff --git a/src/PhaseSpacePoint.cc b/src/PhaseSpacePoint.cc
index 5f75ec6..5d17c5d 100644
--- a/src/PhaseSpacePoint.cc
+++ b/src/PhaseSpacePoint.cc
@@ -1,563 +1,563 @@
#include "RHEJ/PhaseSpacePoint.hh"
#include <random>
#include "RHEJ/resummation_jet_momenta.hh"
#include "RHEJ/Jacobian.hh"
#include "RHEJ/RNGWrapper.hh"
#include "RHEJ/uno.hh"
#include "RHEJ/debug.hh"
#include "RHEJ/kinematics.hh"
namespace RHEJ{
namespace{
//generate Ranlux64Engine with fixed, predefined state
/*
* some (all?) of the Ranlux64Engine constructors leave fields
* uninitialised, invoking undefined behaviour. This can be
* circumvented by restoring the state from a file
*/
CLHEP::Ranlux64Engine gen_Ranlux64Engine(){
static const std::string state =
"9876\n"
"0.91280703978419097666\n"
"0.41606065829518357191\n"
"0.99156342622341142601\n"
"0.030922955274050423213\n"
"0.16206278421638486975\n"
"0.76151768001958330956\n"
"0.43765760066092695979\n"
"0.42904698253748563275\n"
"0.11476317525663759511\n"
"0.026620053590963976831\n"
"0.65953715764414511114\n"
"0.30136722624439826745\n"
"3.5527136788005009294e-15 4\n"
"1 202\n";
const std::string file = std::tmpnam(nullptr);
{
std::ofstream out{file};
out << state;
}
CLHEP::Ranlux64Engine result;
result.restoreStatus(file.c_str());
return result;
}
}
CLHEP::Ranlux64Engine PhaseSpacePoint::ran_{gen_Ranlux64Engine()};
void PhaseSpacePoint::reset_ranlux(std::string const & init_file){
reset_ranlux(init_file.c_str());
}
void PhaseSpacePoint::reset_ranlux(char const * init_file){
ran_.restoreStatus(init_file);
}
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 y_min_user_idx = -3;
constexpr int y_max_user_idx = -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 to 2 jet data
// for 3 jets: -0.0131111 + (1.39385 + 0.050085*delta_y)*delta_y
return -0.0213569 + (1.39765 + 0.0498387*delta_y)*delta_y;
}
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::cluster_jets(
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, jet_def_);
return cs.inclusive_jets(jetptmin_);
}
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){
const double ng_mean = estimate_ng_mean(Born_jets);
std::poisson_distribution<int> dist(ng_mean);
RNGWrapper<CLHEP::Ranlux64Engine> rng{ran_};
const int ng = dist(rng);
assert(ng >= 0);
assert(ng < ng_max);
weight_ *= std::tgamma(ng + 1)*std::exp(ng_mean)*std::pow(ng_mean, -ng);
return ng;
}
namespace {
void copy_AWZH_boson(
std::vector<Sparticle> const & from,
std::vector<Sparticle> & to
){
const auto AWZH_boson = std::find_if(
begin(from), end(from),
[](Sparticle const & p){ return is_AWZH_boson(p); }
);
if(AWZH_boson == end(from)) return;
auto insertion_point = std::lower_bound(
begin(to), end(to), *AWZH_boson, rapidity_less{}
);
to.insert(insertion_point, *AWZH_boson);
assert(std::is_sorted(begin(to), end(to), rapidity_less{}));
}
}
PhaseSpacePoint::PhaseSpacePoint(
- ClusteredEvent const & ev,
+ Event const & ev,
fastjet::JetDefinition jet_def, double jetptmin,
double extpartonptmin, double max_ext_soft_pt_fraction
):
unob_{has_unob_gluon(ev.incoming(), ev.outgoing())},
unof_{!unob_ && has_unof_gluon(ev.incoming(), ev.outgoing())},
extpartonptmin_{extpartonptmin},
max_ext_soft_pt_fraction_{max_ext_soft_pt_fraction},
jetptmin_{jetptmin},
jet_def_{jet_def},
splitter_{jet_def.R(), jet_def, jetptmin, ran_}
{
weight_ = 1;
const auto Born_jets = sorted_by_rapidity(ev.jets());
const int ng = sample_ng(Born_jets);
weight_ /= std::tgamma(ng + 1);
const int ng_jets = sample_ng_jets(ng, Born_jets);
std::vector<fastjet::PseudoJet> out_partons = gen_non_jet(
ng - ng_jets, ccut, jetptmin_
);
{
const auto qperp = std::accumulate(
begin(out_partons), end(out_partons),
fastjet::PseudoJet{}
);
const auto jets = reshuffle(Born_jets, qperp);
if(weight_ == 0.) return;
if(! pass_resummation_cuts(jets)){
weight_ = 0.;
return;
}
std::vector<fastjet::PseudoJet> jet_partons = split(jets, ng_jets);
if(weight_ == 0.) return;
rescale_rapidities(
out_partons,
most_backward_FKL(jet_partons).rapidity(),
most_forward_FKL(jet_partons).rapidity()
);
if(! cluster_jets(out_partons).empty()){
weight_ = 0.;
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{}
);
assert(extremal_FKL_ok(out_partons));
}
if(! jets_ok(Born_jets, out_partons)){
weight_ = 0.;
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(Sparticle{pid::gluon, std::move(p)});
}
most_backward_FKL(outgoing_).type = ev.incoming().front().type;
most_forward_FKL(outgoing_).type = ev.incoming().back().type;
copy_AWZH_boson(ev.outgoing(), outgoing_);
assert(!outgoing_.empty());
reconstruct_incoming(ev.incoming());
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::gen_non_jet(
int count, double ptmin, double ptmax
){
// 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){
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 + RHEJ::PhaseSpacePoint::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 partons;
}
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);
}
}
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_eff = splitter_.R_factor*jet_def_.R();
const int njets = Born_jets.size();
const double p_J = (dy > 2.*R_eff*(njets-1))?
(njets-1)*R_eff*R_eff/(2.*dy):
njets*R_eff*R_eff/(2.*(dy + 2.*R_eff));
return std::min(p_J, 1.);
}
int PhaseSpacePoint::sample_ng_jets(
int ng, std::vector<fastjet::PseudoJet> const & Born_jets
){
RNGWrapper<CLHEP::Ranlux64Engine> 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(rng);
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;
std::vector<fastjet::PseudoJet> jets = resummation_jet_momenta(Born_jets, q);
if(jets.empty()){
weight_ = 0;
return {};
}
// transform delta functions to integration over resummation momenta
weight_ /= Jacobian(jets, q);
return jets;
}
std::vector<int> PhaseSpacePoint::distribute_jet_partons(
int ng_jets, std::vector<fastjet::PseudoJet> const & jets
){
size_t first_valid_jet = 0;
size_t num_valid_jets = jets.size();
const double R_eff = 5./3.*jet_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_ && jets[0].delta_R(jets[1]) > R_eff){
++first_valid_jet;
--num_valid_jets;
}
else if(unof_ && 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() == y_min_user_idx;
}
) != 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() == y_min_user_idx;
}
) != 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);
}
} // end anonymous namespace
#endif
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
int ng_jets
){
return split(jets, distribute_jet_partons(ng_jets, jets));
}
bool PhaseSpacePoint::pass_extremal_cuts(
fastjet::PseudoJet const & ext_parton,
fastjet::PseudoJet const & jet
) const{
if(ext_parton.pt() < extpartonptmin_) return false;
return (ext_parton - jet).pt()/ext_parton.pt() < max_ext_soft_pt_fraction_;
}
std::vector<fastjet::PseudoJet> PhaseSpacePoint::split(
std::vector<fastjet::PseudoJet> const & jets,
std::vector<int> const & np
){
assert(! jets.empty());
assert(jets.size() == np.size());
assert(pass_resummation_cuts(jets));
const size_t most_backward_FKL_idx = 0 + unob_;
const size_t most_forward_FKL_idx = jets.size() - 1 - unof_;
std::vector<fastjet::PseudoJet> jet_partons;
// randomly distribute jet gluons among jets
for(size_t i = 0; i < jets.size(); ++i){
weight_ *= splitter_.Split(jets[i], np[i]);
if(weight_ == 0) return {};
assert(
std::all_of(
begin(splitter_.get_jcons()), end(splitter_.get_jcons()),
is_jet_parton
)
);
const auto first_new_parton = jet_partons.insert(
end(jet_partons),
begin(splitter_.get_jcons()), end(splitter_.get_jcons())
);
auto extremal = end(jet_partons);
if((unob_ && i == 0) || i == most_backward_FKL_idx){
// unordered or FKL backward emission
extremal = std::min_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
if(i == most_backward_FKL_idx) extremal->set_user_index(y_min_user_idx);
}
else if((unof_ && i == jets.size() - 1) || i == most_forward_FKL_idx){
// unordered or FKL forward emission
extremal = std::max_element(
first_new_parton, end(jet_partons), rapidity_less{}
);
if(i == most_forward_FKL_idx) extremal->set_user_index(y_max_user_idx);
}
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_FKL_ok(jet_partons)
|| !split_preserved_jets(jets, jet_partons)
){
weight_ = 0.;
return {};
}
return jet_partons;
}
bool PhaseSpacePoint::extremal_FKL_ok(
std::vector<fastjet::PseudoJet> const & partons
) const{
assert(std::is_sorted(begin(partons), end(partons), rapidity_less{}));
return
most_backward_FKL(partons).user_index() == y_min_user_idx
&& most_forward_FKL(partons).user_index() == y_max_user_idx;
}
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 = sorted_by_rapidity(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_];
}
template<class Particle>
Particle const & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> const & partons
) const{
const size_t idx = partons.size() - 1 - unof_;
assert(idx < partons.size());
return partons[idx];
}
template<class Particle>
Particle & PhaseSpacePoint::most_backward_FKL(
std::vector<Particle> & partons
) const{
return partons[0 + unob_];
}
template<class Particle>
Particle & PhaseSpacePoint::most_forward_FKL(
std::vector<Particle> & partons
) const{
const size_t idx = partons.size() - 1 - unof_;
assert(idx < partons.size());
return partons[idx];
}
namespace{
bool contains_idx(
fastjet::PseudoJet const & jet, fastjet::PseudoJet const & parton
){
auto const & constituents = jet.constituents();
const int idx = parton.user_index();
return std::find_if(
begin(constituents), end(constituents),
[idx](fastjet::PseudoJet const & con){return con.user_index() == idx;}
) != end(constituents);
}
}
/**
* final jet test:
* - number of jets must match Born kinematics
* - no partons designated as nonjet may end up inside jets
* - all other outgoing partons *must* end up inside jets
* - the extremal (in rapidity) partons must be inside the extremal jets
* - rapidities must be the same (by construction)
*/
bool PhaseSpacePoint::jets_ok(
std::vector<fastjet::PseudoJet> const & Born_jets,
std::vector<fastjet::PseudoJet> const & partons
) const{
fastjet::ClusterSequence cs(partons, jet_def_);
const auto jets = sorted_by_rapidity(cs.inclusive_jets(jetptmin_));
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;
for(size_t i = 0; i < jets.size(); ++i){
assert(nearby_ep(jets[i].rapidity(), Born_jets[i].rapidity(), 1e-2));
}
return true;
}
namespace{
}
void PhaseSpacePoint::reconstruct_incoming(
std::array<Sparticle, 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(1e-7));
}
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(double ep) const{
fastjet::PseudoJet diff;
for(auto const & in: incoming()) diff += in.p;
for(auto const & out: outgoing()) diff -= out.p;
return nearby_ep(diff, fastjet::PseudoJet{}, ep);
}
}
diff --git a/src/ScaleFunction.cc b/src/ScaleFunction.cc
index dfe80bc..7db75b3 100644
--- a/src/ScaleFunction.cc
+++ b/src/ScaleFunction.cc
@@ -1,37 +1,37 @@
#include "RHEJ/ScaleFunction.hh"
#include <numeric>
namespace RHEJ{
namespace{
- double ht(ClusteredEvent const & ev){
+ double ht(Event const & ev){
double result = 0.;
for(auto const & parton: ev.outgoing()) result += parton.p.perp();
return result;
}
}
- EventParameters Ht::operator()(ClusteredEvent const & ev) const{
+ EventParameters Ht::operator()(Event const & ev) const{
const double mu = ht(ev);
return {mu, mu, ev.central().weight};
}
- EventParameters MaxJetPperp::operator()(ClusteredEvent const & ev) const{
+ EventParameters MaxJetPperp::operator()(Event const & ev) const{
const double mu = sorted_by_pt(ev.jets()).front().pt();
return {mu, mu, ev.central().weight};
}
- EventParameters JetInvariantMass::operator()(ClusteredEvent const & ev) const{
+ EventParameters JetInvariantMass::operator()(Event const & ev) const{
const double mu = std::accumulate(
begin(ev.jets()), end(ev.jets()), fastjet::PseudoJet{}
).m();
return {mu, mu, ev.central().weight};
}
- EventParameters Product::operator()(ClusteredEvent const &ev) const{
+ EventParameters Product::operator()(Event const &ev) const{
EventParameters s = (*base_scale_)(ev);
s.muf *= factor_;
s.mur *= factor_;
return s;
}
}
diff --git a/src/config.cc b/src/config.cc
index 4f2302f..9428585 100644
--- a/src/config.cc
+++ b/src/config.cc
@@ -1,556 +1,556 @@
#include "RHEJ/config.hh"
#include <set>
#include <string>
#include <vector>
#include <iostream>
#include <stdexcept>
namespace RHEJ{
FileFormat to_FileFormat(std::string const & name){
static const std::map<std::string, FileFormat> known = {
{"Les Houches", FileFormat::Les_Houches},
{"HepMC", FileFormat::HepMC}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown file format " + name);
}
return res->second;
}
std::string extract_suffix(std::string const & filename){
size_t separator = filename.rfind('.');
if(separator == filename.npos) return {};
return filename.substr(separator + 1);
}
FileFormat format_from_suffix(std::string const & filename){
const std::string suffix = extract_suffix(filename);
if(suffix == "lhe") return FileFormat::Les_Houches;
if(suffix == "hepmc3") return FileFormat::HepMC;
throw std::invalid_argument{
"Can't determine format for output file " + filename
};
}
}
namespace YAML {
template<>
struct convert<RHEJ::OutputFile> {
static Node encode(RHEJ::OutputFile const & outfile) {
Node node;
node[to_string(outfile.format)] = outfile.name;
return node;
}
static bool decode(Node const & node, RHEJ::OutputFile & out) {
switch(node.Type()){
case NodeType::Map: {
YAML::const_iterator it = node.begin();
out.format = RHEJ::to_FileFormat(it->first.as<std::string>());
out.name = it->second.as<std::string>();
return true;
}
case NodeType::Scalar:
out.name = node.as<std::string>();
out.format = RHEJ::format_from_suffix(out.name);
return true;
default:
return false;
}
}
};
}
namespace RHEJ{
namespace{
struct invalid_type: std::invalid_argument {
explicit invalid_type(std::string const & what):
std::invalid_argument{what} {};
explicit invalid_type(char const * what):
std::invalid_argument{what} {};
};
struct missing_option: std::logic_error {
explicit missing_option(std::string const & what):
std::logic_error{what} {};
explicit missing_option(char const * what):
std::logic_error{what} {};
};
static const std::set<std::string> known = {
{"trials"},
{"min extparton pt"}, {"max ext soft pt fraction"},
{"resummation jets"}, {"fixed order jets"},
{"FKL"}, {"non-FKL"}, {"unordered"},
{"log correction"},
{"event output"},
{"analysis"},
{"unweight"}, {"RanLux init"},
{"scales"}, {"scale factors"}, {"max scale ratio"}
};
static const std::set<std::string> known_jet = {
{"min pt"}, {"algorithm"}, {"R"}
};
template<typename T>
std::string type_string();
template<>
std::string type_string<std::vector<double>>(){
return "array of doubles";
}
template<>
std::string type_string<std::vector<std::string>>(){
return "array of strings";
}
template<>
std::string type_string<int>(){
return "integer";
}
template<>
std::string type_string<double>(){
return "double";
}
template<>
std::string type_string<bool>(){
return "bool";
}
template<>
std::string type_string<std::string>(){
return "string";
}
template<>
std::string type_string<std::vector<OutputFile>>(){
return "array of output files";
}
fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){
using namespace fastjet;
static const std::map<std::string, fastjet::JetAlgorithm> known = {
{"kt", kt_algorithm},
{"cambridge", cambridge_algorithm},
{"antikt", antikt_algorithm},
{"genkt", genkt_algorithm},
{"cambridge for passive", cambridge_for_passive_algorithm},
{"genkt for passive", genkt_for_passive_algorithm},
{"ee kt", ee_kt_algorithm},
{"ee genkt", ee_genkt_algorithm},
{"plugin", plugin_algorithm}
};
const auto res = known.find(algo);
if(res == known.end()){
throw std::invalid_argument("Unknown jet algorithm " + algo);
}
return res->second;
}
EventTreatment to_EventTreatment(std::string const & name){
static const std::map<std::string, EventTreatment> known = {
{"reweight", EventTreatment::reweight},
{"keep", EventTreatment::keep},
{"discard", EventTreatment::discard}
};
const auto res = known.find(name);
if(res == known.end()){
throw std::invalid_argument("Unknown event treatment " + name);
}
return res->second;
}
// helper struct to allow partial template specialisation
template<typename T>
struct existing_as_helper{
static T as(YAML::Node const & config, std::string const & entry){
assert(config[entry]);
try{
return config[entry].as<T>();
}
catch(std::exception const &){
throw invalid_type{
"Entry " + entry + " is not of type " + type_string<T>()
};
}
}
};
template<>
struct existing_as_helper<fastjet::JetAlgorithm>{
static fastjet::JetAlgorithm as(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
const std::string algo_name =
existing_as_helper<std::string>::as(config, entry);
return to_JetAlgorithm(algo_name);
}
};
template<>
struct existing_as_helper<EventTreatment>{
static EventTreatment as(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
const std::string name =
existing_as_helper<std::string>::as(config, entry);
return to_EventTreatment(name);
}
};
template<>
struct existing_as_helper<OutputFile>{
static OutputFile as(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
YAML::convert<RHEJ::OutputFile> converter{};
OutputFile out;
if(converter.decode(config[entry], out)) return out;
throw std::invalid_argument{
"Bad output file setting: " + config[entry].as<std::string>()
};
}
};
template<typename T>
struct existing_as_helper<std::vector<T>>{
static std::vector<T> as(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
// special case: treat a single value like a vector with one element
if(config[entry].IsScalar()){
return {existing_as_helper<T>::as(config, entry)};
}
try{
return config[entry].as<std::vector<T>>();
}
catch(std::exception const & p){
std::cerr << "Error: " << p.what() << '\n';
throw invalid_type{
"Entry " + entry + " is not of type " + type_string<std::vector<T>>()
};
}
}
};
template<typename T>
T existing_as(YAML::Node const & config, std::string const & entry){
return existing_as_helper<T>::as(config, entry);
}
template<typename T>
T as(YAML::Node const & config, std::string const & entry);
template<>
JetParameters existing_as<JetParameters>(
YAML::Node const & config, std::string const & entry
){
assert(config[entry]);
YAML::Node const & jet_config = config[entry];
JetParameters result;
for(auto const & opt: jet_config){
auto const & opt_name = opt.first.as<std::string>();
if(! known_jet.count(opt_name)){
std::cerr << "In entry " + entry + ": "
+ "Unknown option " + opt_name + "\n";
}
}
try{
result.def = fastjet::JetDefinition{
as<fastjet::JetAlgorithm>(jet_config, "algorithm"),
as<double>(jet_config, "R")
};
result.min_pt = as<double>(jet_config, "min pt");
return result;
}
catch(missing_option const & exc){
throw missing_option{"In entry " + entry + ":\n" + exc.what()};
}
catch(invalid_type const & exc){
throw invalid_type{"In entry " + entry + ":\n" + exc.what()};
}
}
template<typename T>
T as(YAML::Node const & config, std::string const & entry){
if(!config[entry]){
throw missing_option{"No entry for option " + entry};
}
return existing_as<T>(config, entry);
}
template<typename T>
optional<T> as_optional(YAML::Node const & config, std::string const & entry){
if(!config[entry]) return {};
return {existing_as<T>(config, entry)};
}
void update_fixed_order_jet_parameters(
JetParameters & fixed_order_jets, YAML::Node const & yaml
){
YAML::Node const & yaml_jets = yaml["fixed order jets"];
assert(yaml_jets);
for(auto const & opt: yaml_jets){
auto const & opt_name = opt.first.as<std::string>();
if(! known_jet.count(opt_name)){
std::cerr << "In entry fixed order jets: "
"Unknown option " + opt_name + "\n";
}
}
try{
auto algo = as_optional<fastjet::JetAlgorithm>(yaml_jets, "algorithm");
if(algo){
fixed_order_jets.def = fastjet::JetDefinition{
*algo, fixed_order_jets.def.R()
};
}
auto R = as_optional<double>(yaml_jets, "R");
if(R){
fixed_order_jets.def = fastjet::JetDefinition{
fixed_order_jets.def.jet_algorithm(), *R
};
}
auto min_pt = as_optional<double>(yaml_jets, "min pt");
if(min_pt) fixed_order_jets.min_pt = *min_pt;
}
catch(missing_option const & exc){
throw missing_option{
std::string{"In entry fixed order jets:\n"} + exc.what()
};
}
catch(invalid_type const & exc){
throw invalid_type{
std::string{"In entry fixed order jets:\n"} + exc.what()
};
}
}
// like std::stod, but throw if not the whole string can be converted
double to_double(std::string const & str){
std::size_t pos;
const double result = std::stod(str, &pos);
if(pos < str.size()){
throw std::invalid_argument(str + " is not a valid double value");
}
return result;
}
// simple (as in non-composite) scale functions
/**
* An example for a simple scale function would be Ht,
* Ht/2 is then composite (take Ht and then divide by 2)
*/
std::unique_ptr<ScaleFunction> make_simple_ScaleFunction_ptr(
std::string const & scale_fun
){
using ret_type = std::unique_ptr<ScaleFunction>;
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
if(scale_fun == "input") return ret_type{new InputScales{}};
if(scale_fun == "Ht") return ret_type{new Ht{}};
if(scale_fun == "max jet pperp") return ret_type{new MaxJetPperp{}};
if(scale_fun == "jet invariant mass"){
return ret_type{new JetInvariantMass{}};
}
try{
const double scale = to_double(scale_fun);
return ret_type{new FixedScale{scale}};
} catch(std::invalid_argument const &){}
throw std::invalid_argument{"Unknown scale choice: " + scale_fun};
}
std::string trim_front(std::string const & str){
const auto new_begin = std::find_if(
begin(str), end(str), [](char c){ return ! std::isspace(c); }
);
return std::string(new_begin, end(str));
}
std::string trim_back(std::string str){
size_t pos = str.size() - 1;
// use guaranteed wrap-around behaviour to check whether we have
// traversed the whole string
for(; pos < str.size() && std::isspace(str[pos]); --pos) {}
str.resize(pos + 1); // note that pos + 1 can be 0
return str;
}
std::unique_ptr<ScaleFunction> make_ScaleFunction_ptr(
std::string const & scale_fun
){
assert(
scale_fun.empty() ||
(!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back()))
);
const size_t delim = scale_fun.find_first_of("*/");
if(delim == scale_fun.npos){
return make_simple_ScaleFunction_ptr(scale_fun);
}
const std::string first = trim_back(std::string{scale_fun, 0, delim});
const std::string second = trim_front(std::string{scale_fun, delim+1});
double factor;
std::unique_ptr<ScaleFunction> fun;
if(scale_fun[delim] == '/'){
factor = 1/to_double(second);
fun = make_simple_ScaleFunction_ptr(first);
}
else{
assert(scale_fun[delim] == '*');
try{
factor = to_double(second);
fun = make_simple_ScaleFunction_ptr(first);
}
catch(std::invalid_argument const &){
factor = to_double(first);
fun = make_simple_ScaleFunction_ptr(second);
}
}
assert(fun != nullptr);
return std::unique_ptr<ScaleFunction>{
new Product{factor, std::move(fun)}
};
}
ScaleConfig to_ScaleConfig(YAML::Node const & yaml){
ScaleConfig config;
const auto scales = as<std::vector<std::string>>(yaml, "scales");
config.scales.reserve(scales.size());
std::transform(
begin(scales), end(scales), std::back_inserter(config.scales),
make_ScaleFunction_ptr
);
auto scale_factors = as_optional<std::vector<double>>(
yaml, "scale factors"
);
if(scale_factors) config.scale_factors = std::move(*scale_factors);
const auto max_scale_ratio = as_optional<double>(yaml, "max scale ratio");
static_assert(
std::numeric_limits<double>::has_infinity, "infinity not supported"
);
config.max_scale_ratio = max_scale_ratio?
*max_scale_ratio:
std::numeric_limits<double>::infinity()
;
return config;
}
EventTreatMap get_event_treatment(
YAML::Node const & yaml
){
using namespace event_class;
EventTreatMap treat {
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard}
};
treat.emplace(FKL, as<EventTreatment>(yaml, "FKL"));
const auto unordered_treatment = as<EventTreatment>(yaml, "unordered");
treat.emplace(unordered_forward, unordered_treatment);
treat.emplace(unordered_backward, unordered_treatment);
treat.emplace(nonFKL, as<EventTreatment>(yaml, "non-FKL"));
if(treat[nonFKL] == EventTreatment::reweight){
throw std::invalid_argument{"Cannot reweight non-FKL events"};
}
return treat;
}
Config to_Config(YAML::Node const & yaml){
for(auto const & opt: yaml){
auto const & opt_name = opt.first.as<std::string>();
if(! known.count(opt_name)){
std::cerr << "WARNING: unknown option " + opt_name + "\n";
}
}
Config config;
config.resummation_jets = as<JetParameters>(yaml, "resummation jets");
config.fixed_order_jets = config.resummation_jets;
if(yaml["fixed order jets"]){
update_fixed_order_jet_parameters(config.fixed_order_jets, yaml);
}
config.min_extparton_pt = as<double>(yaml, "min extparton pt");
config.max_ext_soft_pt_fraction = yaml["max ext soft pt fraction"]?
as<double>(yaml, "max ext soft pt fraction"):
std::numeric_limits<double>::infinity()
;
config.trials = as<int>(yaml, "trials");
config.log_correction = as<bool>(yaml, "log correction");
config.unweight = as<bool>(yaml, "unweight");
config.treat = get_event_treatment(yaml);
if(yaml["event output"]){
config.output = as<std::vector<OutputFile>>(yaml, "event output");
}
config.RanLux_init = as_optional<std::string>(yaml, "RanLux init");
config.analysis_parameters = yaml["analysis"]?yaml["analysis"]:YAML::Node{};
config.scale_gen = ScaleGenerator{to_ScaleConfig(yaml)};
return config;
}
} // anonymous namespace
Config load_config(std::string const & config_file){
try{
return to_Config(YAML::LoadFile(config_file));
}
catch(...){
std::cerr << "Error reading " << config_file << ":\n ";
throw;
}
}
- ClusteredEvent ScaleGenerator::operator()(ClusteredEvent ev) const{
+ Event ScaleGenerator::operator()(Event ev) const{
assert(ev.variations().empty());
assert(! settings_.scales.empty());
auto const & scales = settings_.scales;
auto const & scale_factors = settings_.scale_factors;
const double max_scale_ratio = settings_.max_scale_ratio;
// check if we are doing scale variation at all
if(scales.size() == 1 && scale_factors.empty()){
ev.central() = (*scales.front())(ev);
return ev;
}
for(auto && base_scale: scales){
const EventParameters cur_base = (*base_scale)(ev);
ev.variations().emplace_back(cur_base);
//multiplicative scale variation
for(double mur_factor: scale_factors){
const double mur = cur_base.mur*mur_factor;
for(double muf_factor: scale_factors){
if(muf_factor == 1. && mur_factor == 1.) continue;
const double muf = cur_base.muf*muf_factor;
if(mur/muf < 1/max_scale_ratio || mur/muf > max_scale_ratio) continue;
ev.variations().emplace_back(
EventParameters{mur, muf, cur_base.weight}
);
}
}
};
ev.central() = (*scales.front())(ev);
return ev;
}
}
diff --git a/src/event_classes.cc b/src/event_classes.cc
index cea54dd..fd5920f 100644
--- a/src/event_classes.cc
+++ b/src/event_classes.cc
@@ -1,140 +1,140 @@
#include "RHEJ/event_classes.hh"
#include <array>
#include <vector>
#include <algorithm>
#include "RHEJ/Event.hh"
namespace RHEJ{
using EventClass = event_class::EventClass;
namespace{
// 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<Sparticle> 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, [](Sparticle 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, [](Sparticle 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,
[](Sparticle const & p){ return p.type == pid::gluon; }
);
}
bool is_FKL(
std::array<Sparticle, 2> const & incoming,
std::vector<Sparticle> 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(ClusteredEvent const & event){
+ bool has_2_jets(Event const & event){
return event.jets().size() >= 2;
}
- bool is_unordered_backward(ClusteredEvent const & ev){
+ 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.
if(is_AWZH_boson(out[1])) return false;
if(!is_FKL(in, {next(begin(out)), 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;
}
// return event with all pz -> -pz
- ClusteredEvent mirrored(ClusteredEvent const & orig_ev){
+ Event mirrored(Event const & orig_ev){
UnclusteredEvent mirrored = orig_ev.unclustered();
for(auto & in: mirrored.incoming){
in.p.reset_momentum(in.p.px(), in.p.py(), -in.p.pz(), in.p.E());
}
std::swap(mirrored.incoming[0], mirrored.incoming[1]);
for(auto & out: mirrored.outgoing){
out.p.reset_momentum(out.p.px(), out.p.py(), -out.p.pz(), out.p.E());
}
std::reverse(begin(mirrored.outgoing), end(mirrored.outgoing));
return {mirrored, orig_ev.jet_def(), orig_ev.min_jet_pt()};
}
- bool is_unordered_forward(ClusteredEvent const & ev){
+ bool is_unordered_forward(Event const & ev){
// check whether the mirrored event is FKL with unordered backward emission
return is_unordered_backward(mirrored(ev));
}
}
- EventClass classify(ClusteredEvent const & ev){
+ EventClass classify(Event const & ev){
if(! final_state_ok(ev.outgoing())) return EventClass::bad_final_state;
if(! has_2_jets(ev)) return EventClass::no_2_jets;
if(is_FKL(ev.incoming(), ev.outgoing())) return EventClass::FKL;
if(is_unordered_backward(ev)){
return EventClass::unordered_backward;
}
if(is_unordered_forward(ev)){
return EventClass::unordered_forward;
}
return EventClass::nonFKL;
}
}
diff --git a/src/main.cc b/src/main.cc
index 370587c..9f8e867 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -1,143 +1,143 @@
/**
* Name: main.cc
* Authors: Tuomas Hapola, Andreas Maier <andreas.maier@durham.ac.uk>
*
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include "yaml-cpp/yaml.h"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/EventReweighter.hh"
#include "RHEJ/config.hh"
#include "RHEJ/stream.hh"
int event_number(std::string const & record){
size_t start = record.rfind("Number of Events");
start = record.find_first_of("123456789", start);
if(start == std::string::npos) {
throw std::invalid_argument("no event number record found");
}
const size_t end = record.find_first_not_of("0123456789", start);
return std::stoi(record.substr(start, end - start));
}
RHEJ::Config load_config(char const * filename){
try{
return RHEJ::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
int main(int argn, char** argv) {
using clock = std::chrono::system_clock;
if (argn < 3) {
std::cerr << "\n# Usage:\n./HEJ config_file input_file\n\n";
return EXIT_FAILURE;
}
const auto start_time = clock::now();
// read configuration
const RHEJ::Config config = load_config(argv[1]);
RHEJ::istream in{argv[2]};
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
LHEF::Reader reader{in};
RHEJ::CombinedEventWriter writer{config.output, reader.heprup};
double global_reweight = 1.;
int max_events = std::numeric_limits<int>::max();
if(argn > 3){
max_events = std::stoi(argv[3]);
const int input_events = event_number(reader.headerBlock);
global_reweight = input_events/static_cast<double>(max_events);
std::cout << "processing " << max_events
<< " out of " << input_events << " events\n";
}
RHEJ::EventReweighter rhej{
reader.heprup,
config.resummation_jets.def, config.resummation_jets.min_pt,
config.treat,
config.min_extparton_pt, config.max_ext_soft_pt_fraction,
config.log_correction
};
if(config.RanLux_init){
RHEJ::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
}
int nevent = 0;
std::array<int, RHEJ::event_class::last_class + 1>
nevent_class{0}, nfailed_class{0};
// Loop over the events in the inputfile
while(reader.readEvent()){
// reweight events so that the total cross section is conserved
reader.hepeup.setWeight(0, global_reweight * reader.hepeup.weight());
if(nevent == max_events) break;
nevent++;
if (nevent % 10000 == 0)
std::cout << "event number " << nevent << std::endl;
// calculate rHEJ weight
- RHEJ::ClusteredEvent FO_event{
+ RHEJ::Event FO_event{
RHEJ::UnclusteredEvent{reader.hepeup},
config.fixed_order_jets.def, config.fixed_order_jets.min_pt,
};
auto resummed_events = rhej.reweight(
FO_event,
config.trials, config.scale_gen
);
const auto ev_class = rhej.last_event_class();
++nevent_class[ev_class];
if(resummed_events.empty()) ++nfailed_class[ev_class];
for(auto const & ev: resummed_events){
analysis->fill(ev);
writer.write(ev);
}
} // main event loop
using namespace RHEJ::event_class;
std::cout<< "Events processed: " << nevent << '\n';
for(size_t ev_class = first_class; ev_class <= last_class; ++ev_class){
std::cout << '\t' << names[ev_class] << ": " << nevent_class[ev_class]
<< ", failed to reconstruct " << nfailed_class[ev_class]
<< '\n';
}
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
}
diff --git a/t/check_res.cc b/t/check_res.cc
index dfc5956..85fa950 100644
--- a/t/check_res.cc
+++ b/t/check_res.cc
@@ -1,70 +1,70 @@
#include <iostream>
#include "LHEF/LHEF.h"
#include "RHEJ/stream.hh"
#include "RHEJ/EventReweighter.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 extpartonptmin = 30;
constexpr double max_ext_soft_pt_fraction =
std::numeric_limits<double>::infinity();
constexpr double jetptmin = 35;
constexpr bool log_corr = false;
using EventTreatment = RHEJ::EventTreatment;
using namespace RHEJ::event_class;
RHEJ::EventTreatMap treat{
{no_2_jets, EventTreatment::discard},
{bad_final_state, EventTreatment::discard},
{nonFKL, EventTreatment::discard},
{unof, EventTreatment::discard},
{unob, EventTreatment::discard},
{FKL, EventTreatment::reweight}
};
};
int main(int argn, char** argv) {
if(argn == 5 && std::string(argv[4]) == "uno"){
--argn;
treat[unof] = EventTreatment::reweight;
treat[unob] = 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]);
RHEJ::istream in{argv[1]};
LHEF::Reader reader{in};
RHEJ::EventReweighter rhej{
reader.heprup,
jet_def, jetptmin,
treat,
extpartonptmin, max_ext_soft_pt_fraction,
log_corr
};
double xsec = 0.;
while(reader.readEvent()){
- RHEJ::ClusteredEvent ev{
+ RHEJ::Event ev{
RHEJ::UnclusteredEvent{reader.hepeup},
Born_jet_def, Born_jetptmin
};
auto resummed_events = rhej.reweight(ev, 10);
for(auto const & ev: resummed_events) xsec += ev.central().weight;
}
if(std::abs(xsec - xsec_ref) > tolerance){
std::cerr << "cross section is off: "
<< xsec << " != " << xsec_ref
<< " +- " << tolerance << '\n';
return EXIT_FAILURE;
}
}
diff --git a/t/test_ME_h_3j.cc b/t/test_ME_h_3j.cc
index ca42703..874bfbc 100644
--- a/t/test_ME_h_3j.cc
+++ b/t/test_ME_h_3j.cc
@@ -1,84 +1,84 @@
#include "LHEF/LHEF.h"
#include "RHEJ/MatrixElement.hh"
#include "RHEJ/Event.hh"
struct Event{
std::array<RHEJ::Sparticle, 2> incoming;
std::vector<RHEJ::Sparticle> outgoing;
Event(
RHEJ::Sparticle in1, RHEJ::Sparticle in2,
std::initializer_list<RHEJ::Sparticle> out
):
incoming{std::move(in1), std::move(in2)},
outgoing{out}
{}
};
constexpr
RHEJ::pid::ParticleID PDG(int id){
return static_cast<RHEJ::pid::ParticleID>(id);
}
constexpr double alpha_s = 0.113559;
constexpr double mu = 125.;
constexpr double R = 0.4;
constexpr double min_jet_pt = 50.;
constexpr auto jet_def = fastjet::antikt_algorithm;
constexpr double ep = 1e-4;
void dump(Event const & ev){
LHEF::Writer writer{std::cerr};
std::cerr << std::setprecision(15);
RHEJ::UnclusteredEvent tmp;
tmp.incoming = ev.incoming;
tmp.outgoing = ev.outgoing;
tmp.central = {mu, mu, 0.};
- RHEJ::ClusteredEvent out_ev{std::move(tmp), jet_def, min_jet_pt};
+ RHEJ::Event out_ev{std::move(tmp), jet_def, min_jet_pt};
writer.hepeup = to_HEPEUP(std::move(out_ev), nullptr);
writer.writeEvent();
}
int main(){
std::vector<Event> events = {
#include "ME_test_events.dat"
};
const std::vector<double> ref_weights{
1.20840366448943e-08,3.41482815561959e-18,4.09113550770465e-12,4.62704699562129e-15,4.48439126356004e-15,1.87972263060254e-12,1.44306589954128e-14,9.10660553817556e-13,6.78185181417658e-13,2.2313219279205e-14,5.04225136522914e-10,1.15445795400542e-13,4.25305992011228e-15,6.34464495221583e-11,8.1504882635896e-09,8.13050604628519e-10,1.16718768515638e-12,2.10641257506331e-08,6.62448031524987e-23,4.3319311109536e-07,1.20944494921349e-15,5.92539274973405e-05,4.4796169513938e-17,1.49273221119393e-15,2.43311569643869e-18,2.402863950606e-12,1.0620402696538e-11,5.35907540952987e-14,1.23031328706173e-12,2.79984692016006e-15,9.15681008462036e-19,6.82192590709917e-19,7.89237067022333e-13,1.86125787460036e-14,6.46979689455398e-07,3.04911830916725e-15,8.13163686285418e-15,1.53766124466997e-14,8.32276518653951e-14,8.88400099785712e-11,8.91102728936822e-06,1.80107644544759e-15,2.10582866365414e-09,3.10947796493188e-17,3.91243813845236e-15,3.26654480787734e-17,9.5447871679842e-14,2.59793669360488e-15,4.96134344054012e-13,5.81162371703576e-14,1.111167877034e-09,2.79109572058797e-16,4.46160513513727e-16,6.75218397576417e-15,2.68260561114305e-11,4.28989454505788e-16,5.8329868340234e-16,2.2024897389957e-09,2.57397955688743e-15,2.44654345522128e-14,2.30866752912176e-14,1.41827747056589e-07,1.08893147410797e-14,4.03382910318587e-11,3.1389869623674e-10,3.90505557430021e-13,1.08295237971538e-12,7.454591488958e-15,7.54483306433453e-14,1.02908259302562e-15,5.63531574394432e-14,6.30249429292698e-11,4.53762691808614e-13,1.89024041595359e-16,6.8158270448436e-15,1.14399844221906e-12,8.90256332512462e-11,6.06548092223796e-17,9.98362270779565e-12,4.59859716748231e-15,2.15099414771066e-10,2.11591109618363e-13,4.9742178019335e-09,1.33845681137089e-10,2.71186491627684e-16,8.86739836718398e-13,1.8249669781697e-14,1.92215499935624e-13,5.63538105112205e-14,4.14171568963499e-12,2.2717979957039e-15,5.55611156772101e-15,1.28364738716082e-16,5.3851134765245e-12,1.24211103260246e-08,4.83232530709888e-17,0.000521569751168655,1.6806219324588e-15,6.59977348355555e-11,9.1627970683281e-14,2.76820118228537e-12,4.0198982918336e-11,1.48426883756103e-18,1.40114232815435e-18,1.7965991343524e-11,7.6977628175923e-10,4.62196685302415e-14,2.69246069544232e-08,1.51516950568948e-13,6.99496684443045e-13,5.18573546306991e-16,0.0434630460883225,1.9837330600509e-23,2.67017601462925e-12,9.8878856049861e-16,8.72580364190897e-11,9.78094461016429e-06,8.71781003862465e-17,1.1559862791241e-08,4.84083438123941e-11,0.0139742563571754,1.99847601444777e-13,3.04915596353023e-17,5.93151845056274e-16,2.17435618150856e-14,5.33225221812805e-08,3.87773807827851e-14,4.50810202343348e-08,1.21607462589437e-08,3.1594598104258e-13,3.28656273761685e-16,2.50210924939855e-21,2.38126296013289e-13,1.25433376653286e-11,2.14322725557091e-14,1.49146107213853e-09,1.242585653879e-15,3.55954759815053e-17,3.58323190858796e-10,4.51037144634267e-14,7.55931651143187e-09,5.37698824199559e-14,7.38868781004529e-14,1.69357650230295e-11,5.80639629823124e-10,7.4837350747944e-17,6.96956083393151e-16,1.81834207579059e-10,2.82540028744915e-12,4.0395726561383e-13,1.79694539029015e-12,1.31041890677814e-16,1.30438189219159e-07,1.43221371664219e-10,9.114047729275e-16,7.34577355123687e-09,3.22593068239571e-15,4.15681922353836e-15,3.36408379142654e-16,1.44480854726148e-11,1.01363611458694e-13,1.00219358430676e-12,0.00338137229726986,4.29548017717404e-15,2.60395354258443e-14,4.1476296627427e-14,1.06855911096978e-13,8.22439039784712e-11,2.73629652512136e-10,5.15659349602634e-11,1.01584718933351e-13,6.94380398977803e-16,2.18931000347143e-18,3.07933496617785e-09,1.66390942837878e-16,3.52525907142114e-13,3.78080647432571e-14,1.17556440172277e-12,6.59751630435329e-17,1.21755902521401e-08,2.55450721249666e-11,1.53256550428156e-14,1.08380893991568e-12,4.16098986326189e-13,2.90205083341855e-14,5.511460594056e-18,1.57699398911427e-16,7.00104401760928e-13,3.58046798793545e-12,2.01229361268661e-12,3.73614518690553e-07,2.14710426845508e-14,3.70105600162546e-13,7.98071394035332e-15,7.37834963791502e-16,1.60529982053429e-09,1.08256973923745e-12,1.0334353499219e-13,9.49543099778791e-11,1.9115591364306e-12,4.0448657562838e-16,8.71422438179998e-15,1.25631778664749e-16,3.14479096095934e-19,2.64428535684163e-16,1.12107615349188e-11,1.23625298812219e-13,4.01428912903254e-13,2.70797308155832e-11,6.56279497820654e-09,7.13734059806893e-16,0.000782868291989126,3.76614026148548e-16,1.26257000484765e-13,9.58631950869197e-13,8.31708875806124e-11,2.34133301418612e-13,9.56863619261736e-14,4.64712032868663e-12,2.86185199150703e-12,1.33576465609689e-11,1.10079981208014e-15,1.0127254462106e-12,1.1462054990501e-12,6.6945137657857e-11,5.58956247731427e-17,9.55921087729154e-12,2.34217482968982e-13,9.06915368555449e-19,5.71005344841003e-15,7.30755001164554e-09,6.21879071656432e-07,3.18226338142145e-17,7.85187718829052e-15,1.12784019414768e-13,1.26620458266236e-12,3.69177699852588e-14,5.15075992359208e-08,4.92072565553936e-15,1.85781482577606e-18,4.45444190278706e-16,1.26266755594178e-13,2.4288203970545e-16,3.05176309462854e-13,1.10490541122231e-16,8.47749453725764e-16,5.15877353528789e-11,5.24254149207253e-17,1.91355494316452e-10,5.00556344079427e-17,2.17154626918891e-15,8.84358943960976e-15,3.61414423657153e-11,4.82735747825927e-09,1.29064865560832e-13,8.42127739585314e-10,3.92906216459766e-16,6.04446459041361e-17,5.89524177334148e-12,5.46194519536606e-19,1.34088171433487e-11,5.52860041827642e-15,1.75477788501097e-15,1.73667838135065e-16,1.65397026290053e-13,1.12185521855334e-09,6.08386662697057e-17,4.59141650497559e-07,1.5209325219688e-11,2.23037744763897e-13,2.01952877105987e-13,6.8597135743949e-17,1.21998302502671e-11,3.21171910404898e-17,5.26021980482437e-15,7.19705632175135e-19,1.20963419308456e-09,6.88528237531901e-13,1.98848279147693e-11,4.01262000581642e-12,5.20247739003617e-13,1.30428205622828e-14,4.11298291694491e-13,5.52043003142038e-10,1.23093415613032e-16,7.36616975809533e-17,3.45364902177034e-16,4.34985181424009e-08,5.38804276057675e-14,9.48300760869968e-13,1.22337104642822e-09,8.35787040119371e-14,1.55423397182344e-15,1.42203125672685e-13,3.18903277676011e-13,3.88652741986258e-15,7.04896809011702e-18,1.32646740227354e-21,7.99231856556843e-13,5.95802234870006e-09,1.65149306543534e-12,1.33064116256099e-11,1.56832702886549e-12,1.10299934631458e-08,3.67898989067657e-16,7.84682358054423e-17,2.57971628718788e-09,8.33182695406546e-17,1.18568992167241e-10,1.84355449893326e-12,2.55066898328171e-13,2.83673294248314e-14,1.78606190612042e-12,3.58348766614931e-12,4.59943073492537e-16,9.51522624162865e-13,8.0286195911246e-14,1.36695977212813e-16,2.66996900743536e-13,1.18874419459462e-15,4.20697505402443e-07,1.44157427471053e-12,3.88352207261609e-12,1.10352682706254e-15,9.16382540998255e-16,1.6918487768822e-14,6.27001020277801e-13,1.19982339734628e-15,2.08524585760556e-15,3.98635216491517e-11,2.29164410091261e-16,6.59993812570474e-14,9.09720299660866e-19,1.00791143935241e-13,3.48282488442434e-15,3.82462701684145e-11,1.82912920125976e-14,2.35934504263123e-12,3.07876467756239e-11,3.01958860729077e-05,1.34442938440306e-12,1.98172170171169e-11,2.54803366934403e-14,4.63525528427102e-13,1.98623090357032e-11,1.17532063246858e-12,1.62645374532443e-16,2.23093836031326e-11,5.23718977936426e-13,6.30254161980384e-15,5.64457070654369e-16,5.8598582297477e-15,1.54121478101166e-11,2.54757492800786e-11,7.24784320481052e-13,7.10590291356503e-11,1.90012029508776e-12,5.38421849395334e-05,9.2027530054952e-13,0.00802523961864568,1.69971085315301e-12,4.17128853468654e-15,2.58531014487319e-15,6.72871249369096e-12,3.68392430578061e-16,1.56192266889718e-14,8.09758960836169e-16,1.33670176619002e-07,3.40402526854793e-10,7.14029639915786e-15,1.15907463202726e-14
};
const std::vector<double> ref_virt{
0.000407437434887122,0.0616303680115081,0.0619338272462471,0.00113600815259533,0.0611805205486735,0.000565698082480416,0.516976611025519,0.0361140019339901,0.00145895360905874,0.134071980986426,0.000248975936203567,0.199294625870235,0.24358305091454,4.74361172738652e-05,5.89666196640687e-05,0.00531027816302915,0.000200243383490462,0.0805261711344379,0.00120027603203917,0.00201705178179895,0.339181495419129,0.000229863490455131,0.0523929260778993,0.000263057002724373,0.00244241108460007,0.0261467924657004,5.9657372131903e-05,0.200512468910312,0.00810099141722299,0.013973033333956,0.0108154079325531,0.158574821336736,0.00070998197888614,0.000498608284914822,3.73054997210496e-05,0.0589374030166269,0.00230622721548646,0.0339477696166811,0.000912463632412847,0.000736710687806361,0.00409015383244036,0.245951477936874,0.0837392661445892,0.222078993106222,0.00532590926455614,0.625630979237125,0.0496003388518398,0.10116468200925,0.0661786669164834,0.0210544045965245,0.0244625178289293,0.0146510511459362,0.101257018174142,0.00316763508016474,0.00143243489499069,0.333592621341247,0.00370907687545401,0.00213872056006003,0.01900030889465,0.167144065343861,2.6100961817368e-05,4.91550277446228e-05,0.0166499896061438,0.00011398434602661,0.000359768751250008,0.00225147558163185,0.0317147645964709,0.0842595514695342,0.35407612826878,0.00010526694382132,0.0430608610020923,3.5715013356257e-05,0.000838215519554185,0.00192018265054966,0.0713712284182796,0.000485232329561933,0.00152471141515176,0.000311862603200811,0.166619420201183,0.0199021050517686,1.41779057397457e-05,0.0017318628596011,0.00053760812759845,0.0686125677803084,0.0532620439395034,0.00166666380089786,0.0113533421860016,0.000109349820296255,0.0115667016173654,0.00173679164546385,0.0256607191277514,0.00120155627964714,0.0111815333349668,0.00349052851306921,0.00531735228893806,0.00152975285273546,4.66037248418355e-05,0.0425850611539657,0.00117697776465851,1.23387751381228e-05,0.0038237470220905,0.0905644410585002,0.02820127090728,0.00024367136851271,0.147636957219027,0.000593136877076053,0.0618471045067722,9.07687985864349e-06,0.000158960279301719,0.0887358461047651,0.0332353914977947,7.4341716367417e-06,0.000892008177580303,0.000741022113249675,0.00237741099974918,0.00354813689304541,0.00074173961384706,0.00467141577324516,0.00842296482521057,0.00821661392007147,4.83792120557372e-05,0.00425252516494632,0.000432204419963046,0.00586736518573752,0.000956530707046858,0.000236851602228008,0.0892427215701044,1.38830962250345e-05,0.0080267826500828,0.284615916615109,0.000228095910528196,0.000135342214864163,7.55341699407512e-06,0.000571884401074523,0.0101768016413147,0.0171951107578717,0.00189669876455143,0.00394930354845142,0.00328294862282636,0.00123193340562429,0.00331783483858544,0.00319051646342221,0.00155129659005506,8.87358953340508e-05,0.27623170363077,0.000223827470760779,0.0568015505265356,0.000286151801691539,0.0193226057191209,0.00123755083758374,0.000465248541258071,0.0760142560624272,0.00112490607994365,0.000144842352216183,0.197405681090406,0.00129883979140284,0.0930965339402953,0.0010308437485118,0.670013359364299,0.00432667098624779,0.0454915871138651,0.00251703872973342,0.00423441033466031,0.000345157767746845,0.00161415356368555,0.00416069958683314,0.00164939718639514,0.0111406419016027,0.00655421252542946,0.00851532238518754,0.00348709627275328,0.496059905302338,0.0371917367183897,0.000211531938236732,0.000258942974591912,0.17759337052515,0.00418622722335427,0.166120697713642,0.00360502919950296,0.000407414278464457,0.0360738279541873,0.0288893174873495,0.00047259015965458,1.56547248986137e-05,0.834956154090444,0.00732703650642937,0.00474657800944094,0.00273096554923287,0.257207154892208,0.00469808491585188,0.000270576144489275,0.000136075839471612,0.097205919234904,0.000607126031374257,0.000136619984636834,0.000160087925098193,0.0073858227644151,0.538982397271589,0.00336284303209108,0.00152656630103797,0.367113302445716,7.48038626150403e-05,2.60037367946894e-05,4.68016275421979e-05,0.000158257654442151,0.000139074963089943,0.000665698879408135,0.000343142682314476,0.0686327476602143,0.000242720484154181,0.0343708294206813,0.00114889968298977,0.0586454267607763,0.0681515891243611,0.00825429790468652,0.023992046337532,0.0176721744671367,0.000118812692493978,0.152442537631382,0.017153140472422,0.0064172004743378,3.15311961394686e-05,0.0176719977165195,0.00523267370956132,0.00696013339412051,0.0283129688131134,0.000387153677372906,0.00691199037138631,0.000146275451696295,0.701230156661611,0.000192846626056869,0.000164329434653947,0.00161262788918379,0.138649850457724,0.0298984433917698,0.0318122273046375,3.04827880239274e-05,0.000447390422746664,0.00113026494701542,0.00113285595566244,2.89098339981399e-05,0.00113856357833898,0.0539144571449436,0.0708152617112735,0.000146823107858066,0.0226524231026193,0.00337461891944858,0.000951319767163979,0.00173528759879105,0.027553837556441,0.000644790987241654,0.00202675597720372,0.0904812026317863,0.000121859990107361,0.0778263959001868,0.0078676495423855,0.00106041104696191,0.0460181415526085,0.111418229601818,0.00485578450515022,0.000301615731939099,0.0991074788208953,0.000545043001929007,0.00699096247420379,0.203553031240112,0.00129391506623192,0.0163036481873179,9.84784360850182e-06,0.0011930992876943,0.000233138513444116,0.00327134762224621,0.00934510369593952,0.0257690022271133,0.0586620958023085,0.0044269686967688,0.00458204111385776,0.000575173870997244,0.0196091953126155,0.0172446254918132,0.0023874923763825,0.00203444177442346,0.00590189106361814,0.160798676736245,0.000218465854877618,0.0449130791723244,0.0221907455403962,0.00577930006864582,4.68129848124974e-05,4.39088532565896e-05,0.0957291068496414,0.000222196152405895,0.00244492299475677,0.0649283175892222,0.00246575641284538,0.00395895850035338,0.452609135498808,2.18247701431579e-05,0.000103691495491759,0.0300254483476907,0.000997825085990438,0.00120168803041095,0.000488310402876861,5.02375648302248e-05,0.0889218155230037,0.886196353445691,0.537472442381394,0.000989158503397272,0.00585452117968649,0.000166646041974894,0.0559913557856292,0.271624614149917,1.9604193016588e-05,0.00400455993655262,0.000719953658022318,0.00769221521719457,0.0122222590022294,0.00869539856191144,0.00665334452394116,0.00269456106141176,0.0203144580365219,1.06269256412071e-05,0.00266305904353431,0.000866750257077464,0.0118678506617259,0.00774140571086026,0.472668086097713,0.00449314655439557,0.000363054262255041,0.000123958362353204,0.0593128481039698,6.49758497533422e-05,0.0147733538718899,0.000181422709422017,0.124358363237549,0.0149418567595413,0.000496790567694211,0.00319609174497007,0.171371757877329,0.000132314177013318,0.000118421517677613,0.0384411149022143,0.0560023977555578,0.000545951922843868,0.00923849055867683,7.42750022847054e-05,0.0034657141222543,0.0814951536826408,0.000247436489650107,0.00279135741337402,0.128217103654628,0.000239656955359697,0.0286081920788208,0.0264861644616376,0.000700223039958121,0.000224261722502375,0.206999869456071,0.0245477336638465,0.00020080171595025,0.000534370708560753,0.00114853261836608,0.0140588914482069,0.0491695416950957,0.0219392429988939,0.000187573794241001,2.76073824918907e-05,0.0177905351398601,0.0227102533101537,0.0156338659659345,0.0045130866637909,0.0118789100672146,0.0684121072806861
};
RHEJ::MatrixElement ME{{jet_def, R}, min_jet_pt, false};
assert(events.size() == ref_weights.size());
for(size_t i = 0; i < events.size(); ++i){
std::sort(
begin(events[i].outgoing), end(events[i].outgoing),
RHEJ::rapidity_less{}
);
const double our_ME = ME.tree(
alpha_s, mu, events[i].incoming, events[i].outgoing, true
);
const double deviation = our_ME/ref_weights[i] - 1.;
if(std::abs(deviation) > ep){
std::cerr
<< "Matrix element deviates by factor of " << 1+deviation << ":\n"
"Is " << our_ME << ", should be " << ref_weights[i] << "\n"
"Event:\n";
dump(events[i]);
return EXIT_FAILURE;
}
const double our_virt = ME.virtual_corrections(
alpha_s, mu, events[i].incoming, events[i].outgoing
);
const double virt_deviation = our_virt/ref_virt[i] - 1.;
if(std::abs(virt_deviation) > ep){
std::cerr
<< "Virtual corrections deviate by factor of " << 1+virt_deviation << ":\n"
"Is " << our_virt << ", should be " << ref_virt[i] << "\n"
"Event:\n";
dump(events[i]);
return EXIT_FAILURE;
}
}
}
diff --git a/t/test_classify.cc b/t/test_classify.cc
index 6e40ad9..c888d61 100644
--- a/t/test_classify.cc
+++ b/t/test_classify.cc
@@ -1,51 +1,51 @@
#include "LHEF/LHEF.h"
#include "RHEJ/stream.hh"
#include "RHEJ/event_classes.hh"
#include "RHEJ/Event.hh"
namespace{
constexpr double min_jet_pt = 30.;
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
using namespace RHEJ::event_class;
static const std::vector<EventClass> results{
unob,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,unob,FKL,FKL,FKL,unof,FKL,unob,FKL,
FKL,unob,unob,FKL,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unof,
FKL,FKL,unof,FKL,FKL,FKL,FKL,FKL,unof,FKL,FKL,FKL,unof,FKL,FKL,unob,unof,
FKL,unof,FKL,unob,FKL,FKL,unob,FKL,unob,unof,unob,unof,FKL,FKL,FKL,FKL,FKL,
FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,
FKL,FKL,FKL,unof,FKL,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unob,FKL,
FKL,FKL,FKL,FKL,unob,FKL,unob,unob,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,FKL,unof,unob,FKL
};
}
int main(int argn, char** argv) {
if(argn != 2){
std::cerr << "Usage: test_classify eventfile";
return EXIT_FAILURE;
}
RHEJ::istream in{argv[1]};
LHEF::Reader reader{in};
LHEF::Writer writer{std::cerr};
writer.heprup = reader.heprup;
for(auto const & expected: results){
reader.readEvent();
- const RHEJ::ClusteredEvent ev{
+ const RHEJ::Event ev{
RHEJ::UnclusteredEvent{reader.hepeup},
jet_def, min_jet_pt
};
const auto ev_class = classify(ev);
if(ev_class != expected){
using RHEJ::event_class::names;
writer.hepeup = reader.hepeup;
std::cerr << "wrong classification of event:\n";
writer.writeEvent();
std::cerr << "classified as " << names[ev_class]
<< ", is " << names[expected] << '\n';
return EXIT_FAILURE;
}
}
}
diff --git a/t/test_psp.cc b/t/test_psp.cc
index e763806..4fa423c 100644
--- a/t/test_psp.cc
+++ b/t/test_psp.cc
@@ -1,64 +1,64 @@
#include "LHEF/LHEF.h"
#include "RHEJ/stream.hh"
#include "RHEJ/event_classes.hh"
#include "RHEJ/Event.hh"
#include "RHEJ/PhaseSpacePoint.hh"
namespace{
constexpr int max_trials = 100;
constexpr double extpartonptmin = 45.;
constexpr double max_ext_soft_pt_fraction =
std::numeric_limits<double>::infinity();
const fastjet::JetDefinition jet_def{fastjet::kt_algorithm, 0.4};
constexpr double min_jet_pt = 50;
};
int main(int argn, char** argv) {
if(argn != 2){
std::cerr << "Usage: test_psp eventfile";
return EXIT_FAILURE;
}
RHEJ::istream in{argv[1]};
LHEF::Reader reader{in};
LHEF::Writer writer{std::cerr};
writer.heprup = reader.heprup;
while(reader.readEvent()){
- const RHEJ::ClusteredEvent ev{
+ const RHEJ::Event ev{
RHEJ::UnclusteredEvent{reader.hepeup},
jet_def, min_jet_pt
};
const auto in_class = classify(ev);
for(int trial = 0; trial < max_trials; ++trial){
RHEJ::PhaseSpacePoint psp{
ev, jet_def, min_jet_pt,
extpartonptmin, max_ext_soft_pt_fraction
};
if(psp.weight() != 0){
RHEJ::UnclusteredEvent tmp_ev;
tmp_ev.incoming = psp.incoming();
tmp_ev.outgoing = psp.outgoing();
tmp_ev.central = {0,0,0};
- RHEJ::ClusteredEvent out_ev{
+ RHEJ::Event out_ev{
std::move(tmp_ev),
jet_def, min_jet_pt
};
const auto out_class = classify(out_ev);
if(out_class != in_class){
using RHEJ::event_class::names;
std::cerr << "Wrong class of phase space point:\n"
"original event of class " << names[in_class] << ":\n";
writer.hepeup = reader.hepeup;
writer.writeEvent();
std::cerr << "Phase space point of class " << names[out_class] << ":\n";
writer.hepeup = to_HEPEUP(out_ev, &writer.heprup);
writer.writeEvent();
return EXIT_FAILURE;
}
}
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:15 AM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4222510
Default Alt Text
(169 KB)

Event Timeline