Root/src/RootSimHitReader.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/Root/RootSimHitReader.hpp"
0010
0011 #include "Acts/Definitions/PdgParticle.hpp"
0012 #include "Acts/Utilities/Logger.hpp"
0013 #include "ActsExamples/EventData/SimParticle.hpp"
0014 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0015 #include "ActsFatras/EventData/ProcessType.hpp"
0016
0017 #include <algorithm>
0018 #include <cstdint>
0019 #include <stdexcept>
0020
0021 #include <TChain.h>
0022 #include <TMathBase.h>
0023
0024 namespace ActsExamples {
0025
0026 RootSimHitReader::RootSimHitReader(const RootSimHitReader::Config& config,
0027                                    Acts::Logging::Level level)
0028     : IReader(),
0029       m_cfg(config),
0030       m_logger(Acts::getDefaultLogger(name(), level)) {
0031   m_inputChain = std::make_unique<TChain>(m_cfg.treeName.c_str());
0032
0033   if (m_cfg.filePath.empty()) {
0034     throw std::invalid_argument("Missing input filename");
0035   }
0036   if (m_cfg.treeName.empty()) {
0037     throw std::invalid_argument("Missing tree name");
0038   }
0039
0040   m_outputSimHits.initialize(m_cfg.outputSimHits);
0041
0042   // Set the branches
0043   auto setBranches = [&]<class T>(const auto& keys, T& columns) {
0044     using MappedType = typename std::remove_reference_t<T>::mapped_type;
0045     for (auto key : keys) {
0046       MappedType a{};  // 0 or nullptr
0047       columns.emplace(key, a);
0048     }
0049     for (auto key : keys) {
0050       m_inputChain->SetBranchAddress(key, &columns.at(key));
0051     }
0052   };
0053
0054   setBranches(m_floatKeys, m_floatColumns);
0055   setBranches(m_uint32Keys, m_uint32Columns);
0056   setBranches(m_uint64Keys, m_uint64Columns);
0057   setBranches(m_int32Keys, m_int32Columns);
0058   setBranches(m_vecUint32Keys, m_vecUint32Columns);
0059
0060   // add file to the input chain
0061   m_inputChain->Add(m_cfg.filePath.c_str());
0062   m_inputChain->LoadTree(0);
0063   ACTS_DEBUG("Adding File " << m_cfg.filePath << " to tree '" << m_cfg.treeName
0064                             << "'.");
0065
0066   // Because each hit is stored in a single entry in the root file, we need to
0067   // scan the file first for the positions of the events in the file in order to
0068   // efficiently read the events later on.
0069   // TODO change the file format to store one event per entry
0070
0071   // Disable all branches and only enable event-id for a first scan of the file
0072   m_inputChain->SetBranchStatus("*", false);
0073   m_inputChain->SetBranchStatus("event_id", true);
0074
0075   auto nEntries = static_cast<std::size_t>(m_inputChain->GetEntriesFast());
0076   if (nEntries == 0) {
0077     throw std::runtime_error("Did not find any entries in input file");
0078   }
0079
0080   // Add the first entry
0081   m_inputChain->GetEntry(0);
0082   m_eventMap.push_back({m_uint32Columns.at("event_id"), 0ul, 0ul});
0083
0084   // Go through all entries and store the position of new events
0085   for (auto i = 1ul; i < nEntries; ++i) {
0086     m_inputChain->GetEntry(i);
0087     const auto evtId = m_uint32Columns.at("event_id");
0088
0089     if (evtId != std::get<0>(m_eventMap.back())) {
0090       std::get<2>(m_eventMap.back()) = i;
0091       m_eventMap.push_back({evtId, i, i});
0092     }
0093   }
0094
0095   std::get<2>(m_eventMap.back()) = nEntries;
0096
0097   // Sort by event id
0098   std::ranges::sort(m_eventMap, {},
0099                     [](const auto& m) { return std::get<0>(m); });
0100
0101   // Re-Enable all branches
0102   m_inputChain->SetBranchStatus("*", true);
0103   ACTS_DEBUG("Event range: " << availableEvents().first << " - "
0104                              << availableEvents().second);
0105 }
0106
0107 RootSimHitReader::~RootSimHitReader() = default;
0108
0109 std::pair<std::size_t, std::size_t> RootSimHitReader::availableEvents() const {
0110   return {std::get<0>(m_eventMap.front()), std::get<0>(m_eventMap.back()) + 1};
0111 }
0112
0113 ProcessCode RootSimHitReader::read(const AlgorithmContext& context) {
0114   auto it = std::ranges::find_if(m_eventMap, [&](const auto& a) {
0115     return std::get<0>(a) == context.eventNumber;
0116   });
0117
0118   if (it == m_eventMap.end()) {
0119     // explicitly warn if it happens for the first or last event as that might
0120     // indicate a human error
0121     if ((context.eventNumber == availableEvents().first) &&
0122         (context.eventNumber == availableEvents().second - 1)) {
0123       ACTS_WARNING("Reading empty event: " << context.eventNumber);
0124     } else {
0125       ACTS_DEBUG("Reading empty event: " << context.eventNumber);
0126     }
0127
0128     m_outputSimHits(context, {});
0129
0130     // Return success flag
0131     return ProcessCode::SUCCESS;
0132   }
0133
0134   // lock the mutex
0135   std::lock_guard<std::mutex> lock(m_read_mutex);
0136
0137   ACTS_DEBUG("Reading event: " << std::get<0>(*it)
0138                                << " stored in entries: " << std::get<1>(*it)
0139                                << " - " << std::get<2>(*it));
0140
0141   SimHitContainer hits;
0142   for (auto entry = std::get<1>(*it); entry < std::get<2>(*it); ++entry) {
0143     m_inputChain->GetEntry(entry);
0144
0145     auto eventId = m_uint32Columns.at("event_id");
0146     if (eventId != context.eventNumber) {
0147       break;
0148     }
0149
0150     const Acts::GeometryIdentifier geoid{m_uint64Columns.at("geometry_id")};
0151     const SimBarcode pid =
0152         SimBarcode().withData(*m_vecUint32Columns.at("barcode"));
0153     const auto index = m_int32Columns.at("index");
0154
0155     const Acts::Vector4 pos4 = {
0156         m_floatColumns.at("tx") * Acts::UnitConstants::mm,
0157         m_floatColumns.at("ty") * Acts::UnitConstants::mm,
0158         m_floatColumns.at("tz") * Acts::UnitConstants::mm,
0159         m_floatColumns.at("tt") * Acts::UnitConstants::mm,
0160     };
0161
0162     const Acts::Vector4 before4 = {
0163         m_floatColumns.at("tpx") * Acts::UnitConstants::GeV,
0164         m_floatColumns.at("tpy") * Acts::UnitConstants::GeV,
0165         m_floatColumns.at("tpz") * Acts::UnitConstants::GeV,
0166         m_floatColumns.at("te") * Acts::UnitConstants::GeV,
0167     };
0168
0169     const Acts::Vector4 delta = {
0170         m_floatColumns.at("deltapx") * Acts::UnitConstants::GeV,
0171         m_floatColumns.at("deltapy") * Acts::UnitConstants::GeV,
0172         m_floatColumns.at("deltapz") * Acts::UnitConstants::GeV,
0173         m_floatColumns.at("deltae") * Acts::UnitConstants::GeV,
0174     };
0175
0176     SimHit hit(geoid, pid, pos4, before4, before4 + delta, index);
0177
0178     hits.insert(hit);
0179   }
0180
0181   ACTS_DEBUG("Read " << hits.size() << " hits for event "
0182                      << context.eventNumber);
0183
0184   m_outputSimHits(context, std::move(hits));
0185
0186   // Return success flag
0187   return ProcessCode::SUCCESS;
0188 }
0189
0190 }  // namespace ActsExamples