Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8725641
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
66 KB
Subscribers
None
View Options
diff --git a/analysis-plugins/include/ROOTReader_decl.hh b/analysis-plugins/include/ROOTReader_decl.hh
index b754699..3052693 100644
--- a/analysis-plugins/include/ROOTReader_decl.hh
+++ b/analysis-plugins/include/ROOTReader_decl.hh
@@ -1,101 +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 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::Event to_Event(BranchEntry const & event_entry) const;
+ RHEJ::UnclusteredEvent to_UnclusteredEvent(
+ BranchEntry const & event_entry
+ ) const;
long long idx_;
TTree* events_;
mutable ClusteredEvent 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 16c73b2..ad9de9a 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::Event const & ev){
+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_);
- RHEJ::Event ev = to_Event(entry_);
+ 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};
return ev_;
}
inline
- RHEJ::Event EventIterator::to_Event(
+ RHEJ::UnclusteredEvent EventIterator::to_UnclusteredEvent(
EventIterator::BranchEntry const & event_entry
) const{
- RHEJ::Event ev;
+ 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/cmake/Modules/FindHepMC.cmake b/cmake/Modules/FindHepMC.cmake
index 68043d2..28f5f10 100644
--- a/cmake/Modules/FindHepMC.cmake
+++ b/cmake/Modules/FindHepMC.cmake
@@ -1,106 +1,107 @@
#.rst:
# FindHepMC
# -------
#
# Finds the HepMC library.
#
# This will define the following variables::
#
# HepMC_FOUND - True if the system has the HepMC library
# HepMC_VERSION - The version of the HepMC library which was found
# HepMC_VERSION_MAJOR - Major version (first number in version string)
# HepMC_VERSION_MINOR - Minor version (second number in version string)
# HepMC_VERSION_PATCH - Patch version (third number in version string)
#
# and the following imported targets::
#
# HepMC::HepMC - The HepMC library
#
# To provide a hint for the preferred HepMC installation, set
# ``HepMC_ROOT`` to the directory containing said installation.
#
#=============================================================================
# Copyright 2016 Andreas Maier (andreas.maier@durham.ac.uk)
#
# Distributed under the OSI-approved BSD License (the "License");
# see accompanying file Copyright.txt for details.
#
# This software is distributed WITHOUT ANY WARRANTY; without even the
# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the License for more information.
#=============================================================================
# (To distribute this file outside of CMake, substitute the full
# License text for the above reference.)
find_package(PkgConfig)
pkg_check_modules(PC_HepMC QUIET HepMC)
+unset(HepMC_INCLUDE_DIR)
find_path(HepMC_INCLUDE_DIR
NAMES HepMCDefs.h
PATHS ${HepMC_ROOT} ${PC_HepMC_INCLUDE_DIRS}
PATH_SUFFIXES HepMC
)
find_library(HepMC_LIBRARY
NAMES HepMC
PATHS ${HepMC_ROOT} ${PC_HepMC_LIBRARY_DIRS}
)
if(DEFINED HepMC_INCLUDE_DIR)
file(STRINGS "${HepMC_INCLUDE_DIR}/HepMCDefs.h" _HepMC_DEFS)
file(STRINGS "${HepMC_INCLUDE_DIR}/Version.h" _HepMC_VERS)
string(
REGEX MATCH
"#define *HEPMC_VERSION *\"[^\"]*\""
_HepMC_VERSION_STR "${_HepMC_DEFS} ${_HepMC_VERS}"
)
string(
REGEX MATCH
"[0-9]*\\.[0-9]*\\.[0-9]*"
HepMC_VERSION ${_HepMC_VERSION_STR}
)
string(
REGEX MATCH
"^[0-9]*"
HepMC_VERSION_MAJOR ${HepMC_VERSION}
)
string(
REGEX REPLACE
"^[0-9]*\\.([0-9]*)" "\\1"
HepMC_VERSION_MINOR ${HepMC_VERSION}
)
string(
REGEX MATCH
"[0-9]*$"
HepMC_VERSION_PATCH ${HepMC_VERSION}
)
message(STATUS "HepMC version: ${HepMC_VERSION}")
endif()
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(HepMC
FOUND_VAR HepMC_FOUND
REQUIRED_VARS
HepMC_LIBRARY
HepMC_INCLUDE_DIR
VERSION_VAR HepMC_VERSION
)
if(HepMC_FOUND)
set(HepMC_LIBRARIES ${HepMC_LIBRARY})
set(HepMC_INCLUDE_DIRS ${HepMC_INCLUDE_DIR})
set(HepMC_DEFINITIONS ${PC_HepMC_CFLAGS_OTHER})
endif()
if(HepMC_FOUND AND NOT TARGET HepMC::HepMC)
add_library(HepMC::HepMC UNKNOWN IMPORTED)
set_target_properties(HepMC::HepMC PROPERTIES
IMPORTED_LOCATION "${HepMC_LIBRARY}"
INTERFACE_COMPILE_OPTIONS "${PC_HepMC_CFLAGS_OTHER}"
INTERFACE_INCLUDE_DIRECTORIES "${HepMC_INCLUDE_DIR}"
)
endif()
mark_as_advanced(
HepMC_INCLUDE_DIR
HepMC_LIBRARY
)
diff --git a/include/RHEJ/Event.hh b/include/RHEJ/Event.hh
index b022edb..91b031b 100644
--- a/include/RHEJ/Event.hh
+++ b/include/RHEJ/Event.hh
@@ -1,121 +1,122 @@
/** \file Event.hh
* \brief Details the parameters of an Event.
*
* Contains the EventParameters struct and the Event struct. Also
- * contains the ClusteredEvent class.
+ * 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
+ * \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 Event Event.hh "include/RHEJ/Event.hh"
+ /** \struct UnclusteredEvent Event.hh "include/RHEJ/Event.hh"
* \brief Event Struct which contains Particle Content.
*/
- struct Event{
- //! Default Event Constructor
- Event() = default;
-
- //! Specified Event Constructor
- Event(LHEF::HEPEUP const & hepeup);
- std::array<Sparticle, 2> incoming; /**< Incoming Particles */
- std::vector<Sparticle> outgoing; /**< Outgoing Particles */
- EventParameters central; /**< Central Value? */
- std::vector<EventParameters> variations; /**< Variations? */
+ 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
*
- * This class contains information about a clustered events jets.
+ * This class contains information about a clustered events jets.
*/
class ClusteredEvent{
public:
//! Default ClusteredEvent Constructor
ClusteredEvent() = default;
-
+
//! Specified ClusteredEvent Constructor
ClusteredEvent(
- Event ev,
+ UnclusteredEvent ev,
fastjet::JetDefinition const & jet_def, double min_jet_pt
);
std::vector<fastjet::PseudoJet> jets() const;
- Event const & unclustered() 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:
- Event ev_;
+ UnclusteredEvent ev_;
fastjet::ClusterSequence cs_;
double min_jet_pt_;
};
double shat(ClusteredEvent const & ev);
LHEF::HEPEUP to_HEPEUP(ClusteredEvent const & event, LHEF::HEPRUP *);
}
diff --git a/src/Event.cc b/src/Event.cc
index a8a12f5..0749a72 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;
}
- Event::Event(LHEF::HEPEUP const & hepeup):
+ 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 ev,
+ 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{
return cs_.inclusive_jets(min_jet_pt_);
}
double shat(ClusteredEvent 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 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 bba3878..2180589 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();
- Event to_Event(PhaseSpacePoint const & psp){
- Event result;
+ 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
){
return reweight(input_ev, num_events, identity);
}
std::vector<ClusteredEvent> EventReweighter::reweight(
ClusteredEvent 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,
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;
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_Event(std::move(psp)), jet_def_, jetptmin_
+ 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
) 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
) 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{
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{
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{
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{
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/event_classes.cc b/src/event_classes.cc
index 4151713..cea54dd 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){
return event.jets().size() >= 2;
}
bool is_unordered_backward(ClusteredEvent 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 = orig_ev.unclustered();
+ 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){
// check whether the mirrored event is FKL with unordered backward emission
return is_unordered_backward(mirrored(ev));
}
}
EventClass classify(ClusteredEvent 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 2e02c5e..370587c 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{reader.hepeup},
+ 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 cbbb490..dfc5956 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{reader.hepeup},
+ 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 a636578..ca42703 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::Event tmp;
+ 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};
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 52805e7..6e40ad9 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{
- RHEJ::Event{reader.hepeup},
+ 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 5a34bdd..e763806 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{
- RHEJ::Event{reader.hepeup},
+ 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::Event tmp_ev;
+ RHEJ::UnclusteredEvent tmp_ev;
tmp_ev.incoming = psp.incoming();
tmp_ev.outgoing = psp.outgoing();
tmp_ev.central = {0,0,0};
RHEJ::ClusteredEvent 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
Details
Attached
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
4243609
Default Alt Text
(66 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment