Page MenuHomeHEPForge
Files F8723873

No OneTemporary
Actions

View Options

diff --git a/Changes-API.md b/Changes-API.md
index 8df426b..ab0bf51 100644
--- a/Changes-API.md
+++ b/Changes-API.md
@@ -1,56 +1,59 @@
# Changelog for HEJ API
This log lists only changes on the HEJ API. These are primarily code changes
relevant for calling HEJ as an API. This file should only be read as an addition
to `Changes.md`, where the main features are documented.
## Version 2.X
### 2.X.0
* Made `MatrixElement.tree_kin(...)` and `MatrixElement.tree_param(...)` public
* New class `CrossSectionAccumulator` to keep track of Cross Section of the
different subproccess
* New template struct `Parameters` similar to old `Weights`
- `Weights` are now an alias for `Parameters<double>`. Calling `Weights` did
not change
- `Weights.hh` was replaced by `Parameters.hh`. The old `Weights.hh` header
will be removed in HEJ Version 2.3.0
* Function to multiplication and division of `EventParameters.weight` by double
- This can be combined with `Parameters`, e.g.
`Parameters<EventParameters>*Weights`, see also `Events.parameters()`
- Moved `EventParameters` to `Parameters.hh` header
* Restructured `Event` class
- `Event` can now only be build from a (new) `Event::EventData` class
- Removed default constructor for `Event`
- `Event::EventData` replaces the old `UnclusteredEvent` struct.
- `UnclusteredEvent` is now deprecated, and will be removed in HEJ Version
2.3.0
- Removed `Event.unclustered()` function
- Added new member function `Events.parameters()`, to directly access
(underlying) `Parameters<EventParameters>`
* Added optional Colour charges to particles (`Particle.colour`)
- Colour connection in the HEJ limit can be generated via
`Event.generate_colours` (automatically done in the resummation)
+* New abstact `EventReader` class, as base for reading events from files
+ - Moved LHE file reader to `HEJ::LesHouchesReader`
+ - New `HEJ::HDF5Reader` to read `hdf5` files
## Version 2.0
### 2.0.5
* no further changes to API
### 2.0.4
* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
### 2.0.3
* no further changes to API
### 2.0.2
* no further changes to API
### 2.0.1
* no further changes to API
diff --git a/Changes.md b/Changes.md
index f1d36fb..22ecbcd 100644
--- a/Changes.md
+++ b/Changes.md
@@ -1,43 +1,46 @@
# Changelog
This is the log for changes to the HEJ program. Further changes to the HEJ API
are documented in `Changes-API.md`. If you are using HEJ as a library, please
also read the changes there.
## Version 2.X
### 2.X.0
* Allow multiplication and division of multiple scale functions e.g.
`H_T/2*m_j1j2`
* Print cross sections at end of run
* Follow HepMC convention for particle Status codes: incoming = 11,
decaying = 2, outgoing = 1 (unchanged)
* Partons now have a Colour charge
- Colours are read from and written to LHE files
- For reweighted events the colours are created according to leading colour in
the FKL limit
* Allow changing the regulator lambda in input (`regulator parameter`, only for
advanced users)
+* Added support to read `hdf5` event files suggested in
+ [arXiv:1905.05120](https://arxiv.org/abs/1905.05120) (needs
+ [HighFive](https://github.com/BlueBrain/HighFive))
## 2.0.5
* Fixed event classification for input not ordered in rapidity
### 2.0.4
* Fixed wrong path of `HEJ_INCLUDE_DIR` in `hej-config.cmake`
### 2.0.3
* Fixed parsing of (numerical factor) * (base scale) in configuration
* Don't change scale names, but sanitise Rivet output file names instead
### 2.0.2
* Changed scale names to `"_over_"` and `"_times_"` for proper file names (was
`"/"` and `"*"` before)
### 2.0.1
* Fixed name of fixed-order generator in error message.
diff --git a/include/HEJ/EventReader.hh b/include/HEJ/EventReader.hh
index 205e3ff..ff920b6 100644
--- a/include/HEJ/EventReader.hh
+++ b/include/HEJ/EventReader.hh
@@ -1,44 +1,44 @@
/** \file
* \brief Header file for event reader interface
*
* This header defines an abstract base class for reading events from files.
*
* \authors Jeppe Andersen, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
-#include<string>
-#include<memory>
+#include <memory>
+#include <string>
#include "LHEF/LHEF.h"
namespace HEJ{
class EventData;
// Abstract base class for reading events from files
struct EventReader {
//! Read an event
virtual bool read_event() = 0;
//! Access header text
virtual std::string const & header() const = 0;
//! Access run information
virtual LHEF::HEPRUP const & heprup() const = 0;
//! Access last read event
virtual LHEF::HEPEUP const & hepeup() const = 0;
virtual ~EventReader() = default;
};
//! Factory function for event readers
/**
* @param infile The name of the input file
* @returns A pointer to an instance of an EventReader
* for the input file
*/
std::unique_ptr<EventReader> make_reader(std::string const & filename);
}
diff --git a/include/HEJ/LesHouchesReader.hh b/include/HEJ/LesHouchesReader.hh
index e79cc08..56656f6 100644
--- a/include/HEJ/LesHouchesReader.hh
+++ b/include/HEJ/LesHouchesReader.hh
@@ -1,52 +1,53 @@
/** \file
* \brief Header file for reading events in the Les Houches Event File format.
*
* \authors Jeppe Andersen, Marian Heil, Andreas Maier, Jennifer Smillie
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
-#include "HEJ/EventReader.hh"
+#include "LHEF/LHEF.h"
+
#include "HEJ/Event.hh"
+#include "HEJ/EventReader.hh"
#include "HEJ/stream.hh"
-#include "LHEF/LHEF.h"
namespace HEJ{
//! Class for writing events to a file in the Les Houches Event File format
class LesHouchesReader : public EventReader{
public:
//! Contruct object reading from the given file
explicit LesHouchesReader(std::string const & filename):
stream_{filename},
reader_{stream_}
{}
//! Read an event
bool read_event() override {
return reader_.readEvent();
}
//! Access header text
std::string const & header() const override {
return reader_.headerBlock;
}
//! Access run information
LHEF::HEPRUP const & heprup() const override {
return reader_.heprup;
}
//! Access last read event
LHEF::HEPEUP const & hepeup() const override {
return reader_.hepeup;
}
private:
HEJ::istream stream_;
LHEF::Reader reader_;
};
}
diff --git a/src/EventReader.cc b/src/EventReader.cc
index a141c6c..0067f31 100644
--- a/src/EventReader.cc
+++ b/src/EventReader.cc
@@ -1,19 +1,19 @@
#include "HEJ/EventReader.hh"
-#include "HEJ/LesHouchesReader.hh"
#include "HEJ/HDF5Reader.hh"
+#include "HEJ/LesHouchesReader.hh"
#include "HEJ/utility.hh"
namespace HEJ {
std::unique_ptr<EventReader> make_reader(std::string const & filename) {
try {
return std::make_unique<LesHouchesReader>(filename);
}
catch(std::runtime_error&) {
#if HEJ_BUILD_WITH_HDF5
return std::make_unique<HDF5Reader>(filename);
#else
throw;
#endif
}
}
}
diff --git a/src/HDF5Reader.cc b/src/HDF5Reader.cc
index 346ba72..75fd647 100644
--- a/src/HDF5Reader.cc
+++ b/src/HDF5Reader.cc
@@ -1,289 +1,290 @@
#include "HEJ/HDF5Reader.hh"
#ifdef HEJ_BUILD_WITH_HDF5
+
#include <numeric>
#include "highfive/H5File.hpp"
namespace HEJ {
namespace {
constexpr std::size_t chunk_size = 10000;
struct ParticleData {
std::vector<int> id;
std::vector<int> status;
std::vector<int> mother1;
std::vector<int> mother2;
std::vector<int> color1;
std::vector<int> color2;
std::vector<double> px;
std::vector<double> py;
std::vector<double> pz;
std::vector<double> e;
std::vector<double> m;
std::vector<double> lifetime;
std::vector<double> spin;
};
struct EventRecords {
std::vector<int> particle_start;
std::vector<int> nparticles;
std::vector<int> pid;
std::vector<int> weight;
std::vector<double> scale;
std::vector<double> fscale;
std::vector<double> rscale;
std::vector<double> aqed;
std::vector<double> aqcd;
ParticleData particles;
};
class ConstEventIterator: public std::iterator<
std::bidirectional_iterator_tag,
LHEF::HEPEUP,
std::ptrdiff_t,
const LHEF::HEPEUP*,
LHEF::HEPEUP const &
>
{
public:
using iterator = ConstEventIterator;
friend iterator cbegin(EventRecords const & records) noexcept;
friend iterator cend(EventRecords const & records) noexcept;
iterator& operator++() {
particle_offset_ += records_.get().nparticles[idx_];
++idx_;
return *this;
}
iterator& operator--() {
--idx_;
particle_offset_ -= records_.get().nparticles[idx_];
return *this;
}
iterator operator--(int) {
auto res = *this;
--(*this);
return res;
}
bool operator==(iterator const & other) const {
return idx_ == other.idx_;
}
bool operator!=(iterator other) const {
return !(*this == other);
}
value_type operator*() const {
value_type hepeup{};
auto const & r = records_.get();
hepeup.NUP = r.nparticles[idx_];
hepeup.IDPRUP = r.pid[idx_];
hepeup.XWGTUP = r.weight[idx_];
hepeup.weights.emplace_back(hepeup.XWGTUP, nullptr);
hepeup.SCALUP = r.scale[idx_];
hepeup.scales.muf = r.fscale[idx_];
hepeup.scales.mur = r.rscale[idx_];
hepeup.AQEDUP = r.aqed[idx_];
hepeup.AQCDUP = r.aqcd[idx_];
const size_t start = particle_offset_;
const size_t end = start + hepeup.NUP;
auto const & p = r.particles;
hepeup.IDUP = std::vector<long>(begin(p.id)+start, begin(p.id)+end);
hepeup.ISTUP = std::vector<int>(begin(p.status)+start, begin(p.status)+end);
hepeup.VTIMUP = std::vector<double>(begin(p.lifetime)+start, begin(p.lifetime)+end);
hepeup.SPINUP = std::vector<double>(begin(p.spin)+start, begin(p.spin)+end);
hepeup.MOTHUP.resize(hepeup.NUP);
hepeup.ICOLUP.resize(hepeup.NUP);
hepeup.PUP.resize(hepeup.NUP);
for(size_t i = 0; i < hepeup.MOTHUP.size(); ++i) {
const size_t idx = start + i;
assert(idx < end);
hepeup.MOTHUP[i] = std::make_pair(p.mother1[idx], p.mother2[idx]);
hepeup.ICOLUP[i] = std::make_pair(p.color1[idx], p.color2[idx]);
hepeup.PUP[i] = std::vector<double>{
p.px[idx], p.py[idx], p.pz[idx], p.e[idx], p.m[idx]
};
}
return hepeup;
}
private:
explicit ConstEventIterator(EventRecords const & records):
records_{records} {}
std::reference_wrapper<const EventRecords> records_;
size_t idx_ = 0;
size_t particle_offset_ = 0;
};
ConstEventIterator cbegin(EventRecords const & records) noexcept {
return ConstEventIterator{records};
}
ConstEventIterator cend(EventRecords const & records) noexcept {
auto it =ConstEventIterator{records};
it.idx_ = records.aqcd.size(); // or size of any other records member
return it;
}
}
struct HDF5Reader::HDF5ReaderImpl{
HighFive::File file;
std::size_t event_idx;
std::size_t particle_idx;
std::size_t nevents;
EventRecords records;
ConstEventIterator cur_event;
LHEF::HEPRUP heprup;
LHEF::HEPEUP hepeup;
explicit HDF5ReaderImpl(std::string const & filename):
file{filename},
event_idx{0},
particle_idx{0},
nevents{
file.getGroup("event")
.getDataSet("nparticles") // or any other dataset
.getSpace().getDimensions().front()
},
records{},
cur_event{cbegin(records)},
heprup{},
hepeup{}
{
read_heprup();
read_event_records(chunk_size);
}
void read_heprup() {
const auto init = file.getGroup("init");
init.getDataSet("PDFgroupA").read(heprup.PDFGUP.first);
init.getDataSet("PDFgroupB").read(heprup.PDFGUP.second);
init.getDataSet("PDFsetA").read(heprup.PDFSUP.first);
init.getDataSet("PDFsetB").read(heprup.PDFSUP.second);
init.getDataSet("beamA").read(heprup.IDBMUP.first);
init.getDataSet("beamB").read(heprup.IDBMUP.second);
init.getDataSet("energyA").read(heprup.EBMUP.first);
init.getDataSet("energyB").read(heprup.EBMUP.second);
init.getDataSet("numProcesses").read(heprup.NPRUP);
init.getDataSet("weightingStrategy").read(heprup.IDWTUP);
const auto proc_info = file.getGroup("procInfo");
proc_info.getDataSet("procId").read(heprup.LPRUP);
proc_info.getDataSet("xSection").read(heprup.XSECUP);
proc_info.getDataSet("error").read(heprup.XERRUP);
// TODO: is this identification correct?
proc_info.getDataSet("unitWeight").read(heprup.XMAXUP);
}
std::size_t read_event_records(std::size_t count) {
count = std::min(count, nevents-event_idx);
auto events = file.getGroup("event");
events.getDataSet("nparticles").select({event_idx}, {count}).read(records.nparticles);
assert(records.nparticles.size() == count);
events.getDataSet("pid").select({event_idx}, {count}).read(records.pid);
events.getDataSet("weight").select({event_idx}, {count}).read(records.weight);
events.getDataSet("scale").select({event_idx}, {count}).read(records.scale);
events.getDataSet("fscale").select({event_idx}, {count}).read(records.fscale);
events.getDataSet("rscale").select({event_idx}, {count}).read(records.rscale);
events.getDataSet("aqed").select({event_idx}, {count}).read(records.aqed);
events.getDataSet("aqcd").select({event_idx}, {count}).read(records.aqcd);
const std::size_t particle_count = std::accumulate(
begin(records.nparticles), end(records.nparticles), 0
);
auto pdata = file.getGroup("particle");
auto & particles = records.particles;
pdata.getDataSet("id").select({particle_idx}, {particle_count}).read(particles.id);
pdata.getDataSet("status").select({particle_idx}, {particle_count}).read(particles.status);
pdata.getDataSet("mother1").select({particle_idx}, {particle_count}).read(particles.mother1);
pdata.getDataSet("mother2").select({particle_idx}, {particle_count}).read(particles.mother2);
pdata.getDataSet("color1").select({particle_idx}, {particle_count}).read(particles.color1);
pdata.getDataSet("color2").select({particle_idx}, {particle_count}).read(particles.color2);
pdata.getDataSet("px").select({particle_idx}, {particle_count}).read(particles.px);
pdata.getDataSet("px").select({particle_idx}, {particle_count}).read(particles.py);
pdata.getDataSet("pz").select({particle_idx}, {particle_count}).read(particles.pz);
pdata.getDataSet("e").select({particle_idx}, {particle_count}).read(particles.e);
pdata.getDataSet("m").select({particle_idx}, {particle_count}).read(particles.m);
pdata.getDataSet("lifetime").select({particle_idx}, {particle_count}).read(particles.lifetime);
pdata.getDataSet("spin").select({particle_idx}, {particle_count}).read(particles.spin);
event_idx += count;
particle_idx += particle_count;
return count;
}
};
HDF5Reader::HDF5Reader(std::string const & filename):
impl_{
new HDF5Reader::HDF5ReaderImpl{filename},
HDF5Reader::HDF5ReaderImplDeleter{}
}
{}
bool HDF5Reader::read_event() {
if(impl_->cur_event == cend(impl_->records)) {
// end of active chunk, read new events from file
const auto events_read = impl_->read_event_records(chunk_size);
impl_->cur_event = cbegin(impl_->records);
if(events_read == 0) return false;
}
impl_->hepeup = *impl_->cur_event;
++impl_->cur_event;
return true;
}
namespace {
static const std::string nothing = "";
}
std::string const & HDF5Reader::header() const {
return nothing;
}
LHEF::HEPRUP const & HDF5Reader::heprup() const {
return impl_->heprup;
}
LHEF::HEPEUP const & HDF5Reader::hepeup() const {
return impl_->hepeup;
}
}
#else // no HDF5 support
namespace HEJ {
class HDF5Reader::HDF5ReaderImpl{};
HDF5Reader::HDF5Reader(std::string const &){
throw std::invalid_argument{
"Failed to create HDF5 reader: "
"HEJ 2 was built without HDF5 support"
};
}
bool HDF5Reader::read_event() {
throw std::logic_error{"unreachable"};
}
std::string const & HDF5Reader::header() const {
throw std::logic_error{"unreachable"};
}
LHEF::HEPRUP const & HDF5Reader::heprup() const {
throw std::logic_error{"unreachable"};
}
LHEF::HEPEUP const & HDF5Reader::hepeup() const {
throw std::logic_error{"unreachable"};
}
}
#endif
namespace HEJ {
void HDF5Reader::HDF5ReaderImplDeleter::operator()(HDF5ReaderImpl* p) {
delete p;
}
}
diff --git a/t/test_hdf5.cc b/t/test_hdf5.cc
index b37a564..8b14d31 100644
--- a/t/test_hdf5.cc
+++ b/t/test_hdf5.cc
@@ -1,40 +1,40 @@
#include <iostream>
#include "HEJ/EventReader.hh"
int main(int argc, char** argv) {
if(argc != 2) {
- std::cerr << "Usage: test_hdf5 file.hdf5\n";
+ std::cerr << "Usage: " << argv[0] << " file.hdf5\n";
return EXIT_FAILURE;
}
auto reader = HEJ::make_reader(argv[1]);
if(
reader->heprup().EBMUP != std::make_pair(7000., 7000.)
|| reader->heprup().PDFSUP != std::make_pair(13000, 13000)
) {
std::cerr << "Read incorrect init parameters\n";
return EXIT_FAILURE;
}
int nevent = 0;
while(reader->read_event()) {
++nevent;
if(reader->hepeup().NUP != 13) {
std::cerr << "Read wrong number of particles: "
<< reader->hepeup().NUP << " != 13 in event " << nevent;
return EXIT_FAILURE;
}
for(size_t i = 0; i < 2; ++i) {
for(size_t j = 0; j < 2; ++j) {
if(reader->hepeup().PUP[i][j] != 0) {
std::cerr << "Non-vanishing transverse momentum in incoming particle"
" in event " << nevent;
return EXIT_FAILURE;
}
}
}
}
if(nevent != 51200) {
std::cerr << "Wrong number of events " << nevent << " != 51200\n";
return EXIT_FAILURE;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:49 PM (1 d, 4 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242526
Default Alt Text
(18 KB)

Event Timeline