Page MenuHomeHEPForge
Files F8723888

No OneTemporary
Actions

View Options

diff --git a/src/HDF5Reader.cc b/src/HDF5Reader.cc
index c146cc7..a5b5e62 100644
--- a/src/HDF5Reader.cc
+++ b/src/HDF5Reader.cc
@@ -1,151 +1,259 @@
#include "HEJ/HDF5Reader.hh"
#ifdef HEJ_BUILD_WITH_HDF5
#include "highfive/H5File.hpp"
namespace HEJ {
+ namespace {
+ 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> 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;
LHEF::HEPRUP heprup;
LHEF::HEPEUP hepeup;
- std::size_t nevent;
- std::size_t offset;
+ EventRecords records;
+ ConstEventIterator cur_event;
+
+ HDF5ReaderImpl():
+ heprup{},
+ hepeup{},
+ records{},
+ cur_event{cbegin(records)}
+ {}
+
+ template<typename T>
+ void read_heprup(HighFive::NodeTraits<T> const & from) {
+ const auto init = from.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 = from.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);
+ }
+
+ template<typename T>
+ void read_event_records(HighFive::NodeTraits<T> const & from) {
+ auto events = from.getGroup("event");
+ events.getDataSet("nparticles").read(records.nparticles);
+ events.getDataSet("pid").read(records.pid);
+ events.getDataSet("weight").read(records.weight);
+ events.getDataSet("scale").read(records.scale);
+ events.getDataSet("fscale").read(records.fscale);
+ events.getDataSet("rscale").read(records.rscale);
+ events.getDataSet("aqed").read(records.aqed);
+ events.getDataSet("aqcd").read(records.aqcd);
+ auto pdata = from.getGroup("particle");
+ auto & particles = records.particles;
+ pdata.getDataSet("id").read(particles.id);
+ pdata.getDataSet("status").read(particles.status);
+ pdata.getDataSet("mother1").read(particles.mother1);
+ pdata.getDataSet("mother2").read(particles.mother2);
+ pdata.getDataSet("color1").read(particles.color1);
+ pdata.getDataSet("color2").read(particles.color2);
+ pdata.getDataSet("px").read(particles.px);
+ pdata.getDataSet("px").read(particles.py);
+ pdata.getDataSet("pz").read(particles.pz);
+ pdata.getDataSet("e").read(particles.e);
+ pdata.getDataSet("m").read(particles.m);
+ pdata.getDataSet("lifetime").read(particles.lifetime);
+ pdata.getDataSet("spin").read(particles.spin);
+ }
};
HDF5Reader::HDF5Reader(std::string const & filename):
impl_{
- new HDF5Reader::HDF5ReaderImpl{HighFive::File{filename}, {}, {}, 0, 0},
+ new HDF5Reader::HDF5ReaderImpl{},
HDF5Reader::HDF5ReaderImplDeleter{}
}
{
- const auto init = impl_->file.getGroup("init");
- auto & heprup = impl_->heprup;
- 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 = impl_->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);
+ HighFive::File file{filename};
+ impl_->read_heprup(file);
+ impl_->read_event_records(file);
}
bool HDF5Reader::read_event() {
- auto events = impl_->file.getGroup("event");
- const std::size_t n = impl_->nevent;
- const std::size_t offset = impl_->offset;
- auto & hepeup = impl_->hepeup;
- events.getDataSet("nparticles").select({n}).read(hepeup.NUP);
- events.getDataSet("pid").select({n}).read(hepeup.IDPRUP);
- events.getDataSet("weight").select({n}).read(hepeup.XWGTUP);
- hepeup.weights.emplace_back(hepeup.XWGTUP, nullptr);
- events.getDataSet("scale").select({n}).read(hepeup.SCALUP);
- events.getDataSet("fscale").select({n}).read(hepeup.scales.muf);
- events.getDataSet("rscale").select({n}).read(hepeup.scales.mur);
- events.getDataSet("aqed").select({n}).read(hepeup.AQEDUP);
- events.getDataSet("aqcd").select({n}).read(hepeup.AQCDUP);
- //events.getDataSet("trials").select({n}).read(hepeup.);
- auto particles = impl_->file.getGroup("particle");
- const size_t npart = hepeup.NUP;
- particles.getDataSet("id").select({offset}, {npart}).read(hepeup.IDUP);
- particles.getDataSet("status").select({offset}, {npart}).read(hepeup.ISTUP);
- std::vector<int> mothers1;
- particles.getDataSet("mother1").select({offset}, {npart}).read(mothers1);
- std::vector<int> mothers2;
- particles.getDataSet("mother2").select({offset}, {npart}).read(mothers2);
- std::vector<int> colors1;
- particles.getDataSet("color1").select({offset}, {npart}).read(colors1);
- std::vector<int> colors2;
- particles.getDataSet("color2").select({offset}, {npart}).read(colors2);
- std::vector<double> px;
- particles.getDataSet("px").select({offset}, {npart}).read(px);
- std::vector<double> py;
- particles.getDataSet("py").select({offset}, {npart}).read(py);
- std::vector<double> pz;
- particles.getDataSet("pz").select({offset}, {npart}).read(pz);
- std::vector<double> e;
- particles.getDataSet("e").select({offset}, {npart}).read(e);
- std::vector<double> m;
- particles.getDataSet("m").select({offset}, {npart}).read(m);
- std::vector<double> lifetimes;
- particles.getDataSet("lifetime").select({offset}, {npart}).read(hepeup.VTIMUP);
- std::vector<double> spins;
- particles.getDataSet("spin").select({offset}, {npart}).read(hepeup.SPINUP);
- hepeup.MOTHUP.resize(npart);
- hepeup.ICOLUP.resize(npart);
- hepeup.PUP.resize(npart);
- for(size_t i = 0; i < npart; ++i) {
- hepeup.MOTHUP[i] = std::make_pair(mothers1[i], mothers2[i]);
- hepeup.ICOLUP[i] = std::make_pair(colors1[i], colors2[i]);
- hepeup.PUP[i] = std::vector<double>{px[i], py[i], pz[i], e[i], m[i]};
- }
- ++impl_->nevent;
- impl_->offset += hepeup.NUP;
+ if(impl_->cur_event == cend(impl_->records)) return false;
+ ++impl_->cur_event;
+ if(impl_->cur_event == cend(impl_->records)) return false;
+ impl_->hepeup = *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{
-
- };
+ 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;
}
}

File Metadata

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

Event Timeline