HepMC3/src/HepMC3Reader.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Io/HepMC3/HepMC3Reader.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Utilities/Logger.hpp"
0014 #include "Acts/Utilities/ScopedTimer.hpp"
0015 #include "Acts/Utilities/ThrowAssert.hpp"
0016 #include "ActsExamples/Framework/ProcessCode.hpp"
0017 #include "ActsExamples/Io/HepMC3/HepMC3Metadata.hpp"
0018 #include "ActsExamples/Io/HepMC3/HepMC3Util.hpp"
0019 #include "ActsExamples/Utilities/MultiplicityGenerators.hpp"
0020
0021 #include <filesystem>
0022 #include <memory>
0023
0024 #include <HepMC3/GenEvent.h>
0025 #include <HepMC3/Print.h>
0026 #include <HepMC3/Reader.h>
0027 #include <HepMC3/Units.h>
0028 #include <boost/algorithm/string/join.hpp>
0029
0030 using namespace Acts::UnitLiterals;
0031
0032 namespace ActsExamples {
0033
0034 HepMC3Reader::HepMC3Reader(const HepMC3Reader::Config& cfg,
0035                            Acts::Logging::Level lvl)
0036     : m_cfg(cfg), m_logger(Acts::getDefaultLogger("HepMC3Reader", lvl)) {
0037   if (m_cfg.outputEvent.empty()) {
0038     throw std::invalid_argument("Missing output collection");
0039   }
0040
0041   m_outputEvent.initialize(m_cfg.outputEvent);
0042
0043   // Validate: exactly one of inputPath or inputs must be set
0044   bool hasInputPath = m_cfg.inputPath.has_value();
0045   bool hasInputs = !m_cfg.inputs.empty();
0046
0047   if (!hasInputPath && !hasInputs) {
0048     throw std::invalid_argument(
0049         "HepMC3 reader requires either 'inputPath' or 'inputs' to be set");
0050   }
0051
0052   if (hasInputPath && hasInputs) {
0053     throw std::invalid_argument(
0054         "HepMC3 reader: 'inputPath' and 'inputs' are mutually exclusive. "
0055         "Use 'inputPath' for single file or 'inputs' for multiple files.");
0056   }
0057
0058   // If inputPath is set, create a single input with default multiplicity
0059   // generator
0060   if (hasInputPath) {
0061     Input input;
0062     input.path = m_cfg.inputPath.value();
0063     input.multiplicityGenerator =
0064         std::make_shared<FixedMultiplicityGenerator>(1);
0065
0066     auto reader = HepMC3Util::deduceReader(input.path);
0067     m_inputs.emplace_back(reader, input.path, input.multiplicityGenerator);
0068   } else {
0069     // Use the provided inputs
0070     for (const auto& input : m_cfg.inputs) {
0071       if (!input.multiplicityGenerator) {
0072         throw std::invalid_argument(
0073             "All Input objects must have a multiplicityGenerator set");
0074       }
0075       auto reader = HepMC3Util::deduceReader(input.path);
0076       m_inputs.emplace_back(reader, input.path, input.multiplicityGenerator);
0077     }
0078   }
0079
0080   if (m_cfg.numEvents.has_value()) {
0081     m_eventsRange = {0, m_cfg.numEvents.value()};
0082   } else {
0083     // Need to make a temporary reader to determine the number of events
0084     Acts::ScopedTimer timer("Determining number of events by reading", logger(),
0085                             Acts::Logging::DEBUG);
0086     // This uses the first reader that's configured, with the assumption that
0087     // this is the hard-scatter event
0088     m_eventsRange = {0, determineNumEvents(m_inputs.front().path)};
0089   }
0090
0091   // Check if any input uses a non-Fixed multiplicity generator
0092   m_hasNonFixedMultiplicity =
0093       std::ranges::any_of(m_inputs, [](const auto& input) {
0094         // Check if this is NOT a FixedMultiplicityGenerator
0095         return dynamic_cast<const FixedMultiplicityGenerator*>(
0096                    input.multiplicityGenerator.get()) == nullptr;
0097       });
0098
0099   // Check if randomNumbers is required
0100   // RNG is needed if we have a vertex generator or any non-Fixed multiplicity
0101   // generator
0102   m_needsRng = (m_cfg.vertexGenerator != nullptr) || m_hasNonFixedMultiplicity;
0103
0104   if (m_needsRng && m_cfg.randomNumbers == nullptr) {
0105     throw std::invalid_argument(
0106         "randomNumbers must be set if vertexGenerator or any non-Fixed "
0107         "multiplicityGenerator is used");
0108   }
0109
0110   ACTS_DEBUG("HepMC3Reader: " << m_eventsRange.first << " - "
0111                               << m_eventsRange.second << " events");
0112 }
0113
0114 HepMC3Reader::~HepMC3Reader() = default;
0115
0116 std::string HepMC3Reader::name() const {
0117   return "HepMC3Reader";
0118 }
0119
0120 std::pair<std::size_t, std::size_t> HepMC3Reader::availableEvents() const {
0121   return m_eventsRange;
0122 }
0123
0124 std::size_t HepMC3Reader::determineNumEvents(
0125     const std::filesystem::path& path) const {
0126   std::optional metadata = HepMC3Metadata::readSidecar(path);
0127
0128   if (metadata.has_value()) {
0129     std::size_t numEvents = metadata.value().numEvents;
0130     ACTS_INFO("HepMC3Reader: Found sidecar metadata file for "
0131               << path << " with " << numEvents << " events");
0132     return numEvents;
0133   }
0134
0135   auto reader = HepMC3Util::deduceReader(m_inputs.front().path);
0136   ACTS_INFO(
0137       "HepMC3Reader: Number of events not specified, will read the "
0138       "whole file to determine the number of events. This might take a while "
0139       "and is usually not what you want.");
0140   std::size_t numEvents = 0;
0141   auto event =
0142       std::make_shared<HepMC3::GenEvent>(HepMC3::Units::GEV, HepMC3::Units::MM);
0143   while (!reader->failed()) {
0144     reader->read_event(*event);
0145     if (!reader->failed()) {
0146       ++numEvents;
0147     }
0148   }
0149   return numEvents;
0150 }
0151
0152 std::shared_ptr<HepMC3::GenEvent> HepMC3Reader::makeEvent() {
0153   return std::make_shared<HepMC3::GenEvent>(HepMC3::Units::GEV,
0154                                             HepMC3::Units::MM);
0155 }
0156
0157 ProcessCode HepMC3Reader::skip(std::size_t events) {
0158   using enum ProcessCode;
0159   if (events == 0) {
0160     return SUCCESS;
0161   }
0162
0163   ACTS_DEBUG("Skipping " << events << " logical events");
0164
0165   // Check if all multiplicity generators are Fixed (deterministic)
0166   // Use the precomputed flag
0167   if (m_hasNonFixedMultiplicity) {
0168     ACTS_ERROR(
0169         "Cannot skip events with non-Fixed multiplicityGenerator (e.g., "
0170         "PoissonMultiplicityGenerator). Skipping requires knowing the exact "
0171         "number of physical events to skip from each input file, which is only "
0172         "possible with deterministic (Fixed) multiplicity generators.");
0173     return ABORT;
0174   }
0175
0176   // For each logical event to skip, evaluate the Fixed multiplicity generators
0177   // to determine how many physical events to skip from each input file
0178   for (std::size_t logicalEvent = 0; logicalEvent < events; ++logicalEvent) {
0179     for (const auto& input : m_inputs) {
0180       // Must be FixedMultiplicityGenerator (checked above), so downcast is safe
0181       const auto* fixedGen = static_cast<const FixedMultiplicityGenerator*>(
0182           input.multiplicityGenerator.get());
0183       std::size_t count = fixedGen->n;
0184
0185       ACTS_VERBOSE("Skipping " << count << " events from " << input.path
0186                                << " for logical event "
0187                                << (m_nextEvent + logicalEvent));
0188       if (!input.reader->skip(static_cast<int>(count))) {
0189         ACTS_ERROR("Error skipping " << count << " events from " << input.path);
0190         return ABORT;
0191       }
0192     }
0193   }
0194
0195   m_nextEvent += events;
0196
0197   return SUCCESS;
0198 }
0199
0200 ProcessCode HepMC3Reader::read(const AlgorithmContext& ctx) {
0201   std::shared_ptr<HepMC3::GenEvent> event;
0202
0203   using enum ProcessCode;
0204
0205   if (readSingleFile(ctx, event) != SUCCESS) {
0206     return ABORT;
0207   }
0208
0209   throw_assert(event != nullptr, "Event should not be null");
0210
0211   if (m_cfg.printListing) {
0212     ACTS_VERBOSE("Read event:\n"
0213                  << [&]() {
0214                       std::stringstream ss;
0215                       HepMC3::Print::listing(ss, *event);
0216                       return ss.str();
0217                     }());
0218   }
0219
0220   if (m_cfg.checkEventNumber &&
0221       static_cast<std::size_t>(event->event_number()) != ctx.eventNumber) {
0222     ACTS_ERROR("HepMC3Reader: Event number mismatch. Expected "
0223                << ctx.eventNumber << ", but got " << event->event_number()
0224                << ". You can turn this off in the configuration if your events "
0225                   "were not written in order.");
0226     return ABORT;
0227   }
0228
0229   m_outputEvent(ctx, std::move(event));
0230
0231   return SUCCESS;
0232 }
0233
0234 ProcessCode HepMC3Reader::readSingleFile(
0235     const ActsExamples::AlgorithmContext& ctx,
0236     std::shared_ptr<HepMC3::GenEvent>& outputEvent) {
0237   using enum ProcessCode;
0238   ACTS_VERBOSE("Reading from single file");
0239
0240   std::vector<std::shared_ptr<HepMC3::GenEvent>> events;
0241   {
0242     std::scoped_lock lock(m_queueMutex);
0243
0244     if (m_bufferError) {
0245       ACTS_ERROR("Buffer error (maybe in other thread), aborting");
0246       return ABORT;
0247     }
0248
0249     // Check if we already read this event on another thread
0250     if (!m_events.empty() && ctx.eventNumber < m_nextEvent) {
0251       if (readCached(ctx, events) != SUCCESS) {
0252         ACTS_ERROR("Error reading event " << ctx.eventNumber);
0253         return ABORT;
0254       }
0255     } else {
0256       if (readBuffer(ctx, events) != SUCCESS) {
0257         ACTS_ERROR("Error reading event " << ctx.eventNumber);
0258         return ABORT;
0259       }
0260     }
0261   }
0262
0263   if (m_cfg.vertexGenerator != nullptr) {
0264     Acts::ScopedTimer timer("Shifting events to vertex", logger(),
0265                             Acts::Logging::DEBUG);
0266     auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
0267     for (auto& event : events) {
0268       auto vertexPosition = (*m_cfg.vertexGenerator)(rng);
0269
0270       ACTS_VERBOSE("Shifting event to " << vertexPosition.transpose());
0271       // Our internal time unit is ctau, so is HepMC3's, make sure we convert
0272       // to mm
0273       HepMC3::FourVector vtxPosHepMC(vertexPosition[Acts::eFreePos0] / 1_mm,
0274                                      vertexPosition[Acts::eFreePos1] / 1_mm,
0275                                      vertexPosition[Acts::eFreePos2] / 1_mm,
0276                                      vertexPosition[Acts::eFreeTime] / 1_mm);
0277       event->shift_position_to(vtxPosHepMC);
0278     }
0279   }
0280
0281   outputEvent = makeEvent();
0282
0283   std::vector<const HepMC3::GenEvent*> eventPtrs;
0284   eventPtrs.reserve(events.size());
0285   std::ranges::transform(events, std::back_inserter(eventPtrs),
0286                          [](auto& event) { return event.get(); });
0287   {
0288     Acts::ScopedTimer timer("Merging events", logger(), Acts::Logging::DEBUG);
0289     HepMC3Util::mergeEvents(*outputEvent, eventPtrs, logger());
0290   }
0291
0292   outputEvent->set_event_number(events.front()->event_number());
0293
0294   return SUCCESS;
0295 }
0296
0297 ProcessCode HepMC3Reader::readCached(
0298     const ActsExamples::AlgorithmContext& ctx,
0299     std::vector<std::shared_ptr<HepMC3::GenEvent>>& events) {
0300   ACTS_VERBOSE("Already read event " << ctx.eventNumber);
0301   auto it = std::ranges::find_if(
0302       m_events, [&ctx](auto& elem) { return elem.first == ctx.eventNumber; });
0303
0304   if (it == m_events.end()) {
0305     if (m_cfg.numEvents.has_value()) {
0306       ACTS_ERROR(
0307           "Event "
0308           << ctx.eventNumber
0309           << " was not found in the queue. Most likely the manually "
0310              "configured event count of the input file(s) of "
0311           << m_cfg.numEvents.value()
0312           << " is larger than the number of events found in the file(s)");
0313     } else {
0314       ACTS_ERROR(
0315           "Event " << ctx.eventNumber
0316                    << "  should be in the queue, but is not. This is a bug.");
0317     }
0318     m_bufferError = true;
0319     return ProcessCode::ABORT;
0320   }
0321
0322   auto [num, genEvents] = std::move(*it);
0323   m_events.erase(it);
0324   events = std::move(genEvents);
0325   return ProcessCode::SUCCESS;
0326 }
0327
0328 ProcessCode HepMC3Reader::readBuffer(
0329     const ActsExamples::AlgorithmContext& ctx,
0330     std::vector<std::shared_ptr<HepMC3::GenEvent>>& outputEvents) {
0331   using enum ProcessCode;
0332
0333   ACTS_VERBOSE("Next event to read is: " << m_nextEvent);
0334
0335   std::size_t eventsToRead = ctx.eventNumber - m_nextEvent + 1;
0336   ACTS_VERBOSE("event_number=" << ctx.eventNumber
0337                                << " next_event=" << m_nextEvent
0338                                << " events_to_read=" << eventsToRead);
0339
0340   if (m_events.size() + eventsToRead > m_cfg.maxEventBufferSize) {
0341     ACTS_ERROR("Event buffer size of "
0342                << m_cfg.maxEventBufferSize
0343                << " would be exceeded. This can happen in case there are many "
0344                   "threads and processing happens strongly out-of-order");
0345     m_bufferError = true;
0346     return ABORT;
0347   }
0348
0349   do {
0350     ACTS_VERBOSE("Reading next event as number " << m_nextEvent);
0351     std::size_t thisEvent = m_nextEvent;
0352     m_nextEvent += 1;
0353
0354     std::vector<std::shared_ptr<HepMC3::GenEvent>> events;
0355
0356     if (readLogicalEvent(ctx, events) != SUCCESS) {
0357       ACTS_ERROR("Error reading event " << thisEvent);
0358       return ABORT;
0359     }
0360
0361     if (events.empty()) {
0362       ACTS_ERROR("No events read from file, this is a bug");
0363       return ABORT;
0364     }
0365
0366     m_events.emplace_back(thisEvent, std::move(events));
0367
0368   } while (m_nextEvent <= ctx.eventNumber);
0369
0370   ACTS_VERBOSE("Queue is now: [" << [&]() {
0371     std::vector<std::string> numbers;
0372     std::ranges::transform(
0373         m_events, std::back_inserter(numbers),
0374         [](const auto& v) { return std::to_string(v.first); });
0375     return boost::algorithm::join(numbers, ", ");
0376   }() << "]");
0377
0378   throw_assert(!m_events.empty(), "Event should not be empty, this is a bug.");
0379
0380   m_maxEventBufferSize = std::max(m_maxEventBufferSize, m_events.size());
0381
0382   auto [num, genEvents] = std::move(m_events.back());
0383   m_events.pop_back();
0384
0385   ACTS_VERBOSE("Popping event " << num << " from queue");
0386
0387   outputEvents = std::move(genEvents);
0388
0389   return SUCCESS;
0390 }
0391
0392 ProcessCode HepMC3Reader::readLogicalEvent(
0393     const ActsExamples::AlgorithmContext& ctx,
0394     std::vector<std::shared_ptr<HepMC3::GenEvent>>& events) {
0395   using enum ProcessCode;
0396   ACTS_VERBOSE("Reading logical event " << ctx.eventNumber);
0397   Acts::ScopedTimer timer("Reading logical event", logger(),
0398                           Acts::Logging::DEBUG);
0399
0400   // @TODO: Add the index as an attribute to the event and it's content
0401
0402   // Spawn RNG for multiplicity generators if needed
0403   std::optional<RandomEngine> rng;
0404   if (m_needsRng) {
0405     rng = m_cfg.randomNumbers->spawnGenerator(ctx);
0406   }
0407
0408   for (std::size_t inputIndex = 0; inputIndex < m_inputs.size(); ++inputIndex) {
0409     auto& reader = m_inputs[inputIndex].reader;
0410     auto& path = m_inputs[inputIndex].path;
0411     auto& multiplicityGenerator = m_inputs[inputIndex].multiplicityGenerator;
0412
0413     // Use multiplicityGenerator to determine count
0414     std::size_t count = 0;
0415     if (rng.has_value()) {
0416       count = (*multiplicityGenerator)(*rng);
0417     } else {
0418       // Must be FixedMultiplicityGenerator if no RNG, so downcast is safe
0419       const auto& fixedGen = static_cast<const FixedMultiplicityGenerator&>(
0420           *multiplicityGenerator);
0421       count = fixedGen.n;
0422     }
0423
0424     ACTS_VERBOSE("Reading " << count << " events from " << path);
0425     for (std::size_t i = 0; i < count; ++i) {
0426       auto event = makeEvent();
0427
0428       reader->read_event(*event);
0429       if (reader->failed()) {
0430         ACTS_ERROR("Error reading event " << i << " (input index = "
0431                                           << inputIndex << ") from " << path);
0432         if (inputIndex > 0) {
0433           ACTS_ERROR("-> since this is input file index "
0434                      << inputIndex
0435                      << ", this probably means that the "
0436                         "input file has "
0437                         "fewer events than expected.");
0438         }
0439         return ABORT;
0440       }
0441       events.push_back(std::move(event));
0442       m_inputs[inputIndex].eventsRead++;
0443     }
0444   }
0445
0446   ACTS_VERBOSE("Read " << events.size() << " events in total from all files");
0447
0448   return SUCCESS;
0449 }
0450
0451 ProcessCode HepMC3Reader::finalize() {
0452   ACTS_VERBOSE("Closing " << m_inputs.size() << " input files");
0453
0454   std::size_t totalEventsRead = 0;
0455
0456   // Print summary of events read from each input file
0457   if (m_inputs.size() > 1) {
0458     ACTS_INFO("Events read per input file:");
0459     for (std::size_t i = 0; i < m_inputs.size(); ++i) {
0460       const auto& input = m_inputs[i];
0461       ACTS_INFO("  [" << i << "] " << input.path.filename() << ": "
0462                       << input.eventsRead << " events");
0463       totalEventsRead += input.eventsRead;
0464       input.reader->close();
0465     }
0466     ACTS_INFO("Total physical events read: " << totalEventsRead);
0467   } else {
0468     // Single input file
0469     for (const auto& input : m_inputs) {
0470       ACTS_INFO("Events read from " << input.path.filename() << ": "
0471                                     << input.eventsRead);
0472       input.reader->close();
0473     }
0474   }
0475
0476   ACTS_DEBUG("Maximum event buffer size was: "
0477              << m_maxEventBufferSize << ", limit=" << m_cfg.maxEventBufferSize);
0478   return ProcessCode::SUCCESS;
0479 }
0480
0481 }  // namespace ActsExamples