Csv/src/CsvMeasurementReader.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/Csv/CsvMeasurementReader.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "ActsExamples/Digitization/MeasurementCreation.hpp"
0015 #include "ActsExamples/EventData/Cluster.hpp"
0016 #include "ActsExamples/EventData/GeometryContainers.hpp"
0017 #include "ActsExamples/EventData/Index.hpp"
0018 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0019 #include "ActsExamples/EventData/Measurement.hpp"
0020 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0021 #include "ActsExamples/Io/Csv/CsvInputOutput.hpp"
0022 #include "ActsExamples/Utilities/Paths.hpp"
0023
0024 #include <algorithm>
0025 #include <array>
0026 #include <cstdint>
0027 #include <functional>
0028 #include <iterator>
0029 #include <list>
0030 #include <stdexcept>
0031 #include <vector>
0032
0033 #include "CsvOutputData.hpp"
0034
0035 ActsExamples::CsvMeasurementReader::CsvMeasurementReader(
0036     const ActsExamples::CsvMeasurementReader::Config& config,
0037     Acts::Logging::Level level)
0038     : m_cfg(config),
0039       m_eventsRange(
0040           determineEventFilesRange(m_cfg.inputDir, "measurements.csv")),
0041       m_logger(Acts::getDefaultLogger("CsvMeasurementReader", level)) {
0042   if (m_cfg.outputMeasurements.empty()) {
0043     throw std::invalid_argument("Missing measurement output collection");
0044   }
0045
0046   m_outputMeasurements.initialize(m_cfg.outputMeasurements);
0047   m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap);
0048   m_outputClusters.maybeInitialize(m_cfg.outputClusters);
0049   m_outputMeasurementParticlesMap.maybeInitialize(
0050       m_cfg.outputMeasurementParticlesMap);
0051   m_inputHits.maybeInitialize(m_cfg.inputSimHits);
0052
0053   // Check if event ranges match (should also catch missing files)
0054   auto checkRange = [&](const std::string& fileStem) {
0055     const auto hitmapRange = determineEventFilesRange(m_cfg.inputDir, fileStem);
0056     if (hitmapRange.first > m_eventsRange.first ||
0057         hitmapRange.second < m_eventsRange.second) {
0058       throw std::runtime_error("event range mismatch for 'event**-" + fileStem +
0059                                "'");
0060     }
0061   };
0062
0063   checkRange("measurement-simhit-map.csv");
0064   if (!m_cfg.outputClusters.empty()) {
0065     checkRange("cells.csv");
0066   }
0067 }
0068
0069 std::string ActsExamples::CsvMeasurementReader::CsvMeasurementReader::name()
0070     const {
0071   return "CsvMeasurementReader";
0072 }
0073
0074 std::pair<std::size_t, std::size_t>
0075 ActsExamples::CsvMeasurementReader::availableEvents() const {
0076   return m_eventsRange;
0077 }
0078
0079 namespace {
0080 struct CompareHitId {
0081   // support transparent comparison between identifiers and full objects
0082   using is_transparent = void;
0083   template <typename T>
0084   constexpr bool operator()(const T& left, const T& right) const {
0085     return left.hit_id < right.hit_id;
0086   }
0087   template <typename T>
0088   constexpr bool operator()(std::uint64_t left_id, const T& right) const {
0089     return left_id < right.hit_id;
0090   }
0091   template <typename T>
0092   constexpr bool operator()(const T& left, std::uint64_t right_id) const {
0093     return left.hit_id < right_id;
0094   }
0095 };
0096
0097 struct CompareGeometryId {
0098   bool operator()(const ActsExamples::MeasurementData& left,
0099                   const ActsExamples::MeasurementData& right) const {
0100     return left.geometry_id < right.geometry_id;
0101   }
0102 };
0103
0104 template <typename Data>
0105 inline std::vector<Data> readEverything(
0106     const std::string& inputDir, const std::string& filename,
0107     const std::vector<std::string>& optionalColumns, std::size_t event) {
0108   std::string path = ActsExamples::perEventFilepath(inputDir, filename, event);
0109   ActsExamples::NamedTupleCsvReader<Data> reader(path, optionalColumns);
0110
0111   std::vector<Data> everything;
0112   Data one;
0113   while (reader.read(one)) {
0114     everything.push_back(one);
0115   }
0116
0117   return everything;
0118 }
0119
0120 std::vector<ActsExamples::MeasurementData> readMeasurementsByGeometryId(
0121     const std::string& inputDir, std::size_t event) {
0122   // geometry_id and t are optional columns
0123   auto measurements = readEverything<ActsExamples::MeasurementData>(
0124       inputDir, "measurements.csv", {"geometry_id", "t"}, event);
0125   // sort same way they will be sorted in the output container
0126   std::ranges::sort(measurements, CompareGeometryId{});
0127   return measurements;
0128 }
0129
0130 ActsExamples::ClusterContainer makeClusters(
0131     const std::unordered_multimap<std::size_t, ActsExamples::CellData>&
0132         cellDataMap,
0133     std::size_t nMeasurements) {
0134   using namespace ActsExamples;
0135   ClusterContainer clusters;
0136
0137   for (auto index = 0ul; index < nMeasurements; ++index) {
0138     auto [begin, end] = cellDataMap.equal_range(index);
0139
0140     // Fill the channels with the iterators
0141     Cluster cluster;
0142     cluster.channels.reserve(std::distance(begin, end));
0143
0144     for (auto it = begin; it != end; ++it) {
0145       const auto& cellData = it->second;
0146       ActsFatras::Segmentizer::Segment2D dummySegment = {Acts::Vector2::Zero(),
0147                                                          Acts::Vector2::Zero()};
0148
0149       ActsFatras::Segmentizer::Bin2D bin{
0150           static_cast<unsigned int>(cellData.channel0),
0151           static_cast<unsigned int>(cellData.channel1)};
0152
0153       cluster.channels.emplace_back(bin, dummySegment, cellData.value);
0154     }
0155
0156     // update the iterator
0157
0158     // Compute cluster size
0159     if (!cluster.channels.empty()) {
0160       auto compareX = [](const auto& a, const auto& b) {
0161         return a.bin[0] < b.bin[0];
0162       };
0163       auto compareY = [](const auto& a, const auto& b) {
0164         return a.bin[1] < b.bin[1];
0165       };
0166
0167       auto [minX, maxX] = std::minmax_element(cluster.channels.begin(),
0168                                               cluster.channels.end(), compareX);
0169       auto [minY, maxY] = std::minmax_element(cluster.channels.begin(),
0170                                               cluster.channels.end(), compareY);
0171       cluster.sizeLoc0 = 1 + maxX->bin[0] - minX->bin[0];
0172       cluster.sizeLoc1 = 1 + maxY->bin[1] - minY->bin[1];
0173     }
0174
0175     clusters.push_back(cluster);
0176   }
0177   return clusters;
0178 }
0179
0180 }  // namespace
0181
0182 ActsExamples::ProcessCode ActsExamples::CsvMeasurementReader::read(
0183     const ActsExamples::AlgorithmContext& ctx) {
0184   // hit_id in the files is not required to be neither continuous nor
0185   // monotonic. internally, we want continuous indices within [0,#hits)
0186   // to simplify data handling. to be able to perform this mapping we first
0187   // read all data into memory before converting to the internal event data
0188   // types.
0189   //
0190   // Note: the cell data is optional
0191   auto measurementData =
0192       readMeasurementsByGeometryId(m_cfg.inputDir, ctx.eventNumber);
0193
0194   // Prepare containers for the hit data using the framework event data types
0195   MeasurementContainer tmpMeasurements;
0196   GeometryIdMultimap<ConstVariableBoundMeasurementProxy> orderedMeasurements;
0197   IndexMultimap<Index> measurementSimHitsMap;
0198
0199   tmpMeasurements.reserve(measurementData.size());
0200   orderedMeasurements.reserve(measurementData.size());
0201   // Safe long as we have single particle to sim hit association
0202   measurementSimHitsMap.reserve(measurementData.size());
0203
0204   auto measurementSimHitLinkData =
0205       readEverything<ActsExamples::MeasurementSimHitLink>(
0206           m_cfg.inputDir, "measurement-simhit-map.csv", {}, ctx.eventNumber);
0207   for (auto mshLink : measurementSimHitLinkData) {
0208     measurementSimHitsMap.emplace_hint(measurementSimHitsMap.end(),
0209                                        mshLink.measurement_id, mshLink.hit_id);
0210   }
0211
0212   for (const MeasurementData& m : measurementData) {
0213     Acts::GeometryIdentifier geoId = m.geometry_id;
0214
0215     // Create the measurement
0216     DigitizedParameters dParameters;
0217     for (unsigned int ipar = 0;
0218          ipar < static_cast<unsigned int>(Acts::eBoundSize); ++ipar) {
0219       if (((m.local_key) & (1 << (ipar + 1))) != 0) {
0220         dParameters.indices.push_back(static_cast<Acts::BoundIndices>(ipar));
0221         switch (ipar) {
0222           case static_cast<unsigned int>(Acts::eBoundLoc0): {
0223             dParameters.values.push_back(m.local0);
0224             dParameters.variances.push_back(m.var_local0);
0225           }; break;
0226           case static_cast<unsigned int>(Acts::eBoundLoc1): {
0227             dParameters.values.push_back(m.local1);
0228             dParameters.variances.push_back(m.var_local1);
0229           }; break;
0230           case static_cast<unsigned int>(Acts::eBoundPhi): {
0231             dParameters.values.push_back(m.phi);
0232             dParameters.variances.push_back(m.var_phi);
0233           }; break;
0234           case static_cast<unsigned int>(Acts::eBoundTheta): {
0235             dParameters.values.push_back(m.theta);
0236             dParameters.variances.push_back(m.var_theta);
0237           }; break;
0238           case static_cast<unsigned int>(Acts::eBoundTime): {
0239             dParameters.values.push_back(m.time);
0240             dParameters.variances.push_back(m.var_time);
0241           }; break;
0242           default:
0243             break;
0244         }
0245       }
0246     }
0247
0248     // The measurement container is unordered and the index under which
0249     // the measurement will be stored is known before adding it.
0250     auto measurement = createMeasurement(tmpMeasurements, geoId, dParameters);
0251
0252     // Due to the previous sorting of the raw hit data by geometry id, new
0253     // measurements should always end up at the end of the container. previous
0254     // elements were not touched; cluster indices remain stable and can
0255     // be used to identify the m.
0256     auto inserted = orderedMeasurements.emplace_hint(orderedMeasurements.end(),
0257                                                      geoId, measurement);
0258     if (std::next(inserted) != orderedMeasurements.end()) {
0259       ACTS_FATAL("Something went horribly wrong with the hit sorting");
0260       return ProcessCode::ABORT;
0261     }
0262   }
0263
0264   MeasurementContainer measurements;
0265   for (auto& [_, meas] : orderedMeasurements) {
0266     measurements.emplaceMeasurement(meas.size(), meas.geometryId(), meas);
0267   }
0268
0269   // Generate measurement-particles-map
0270   if (m_inputHits.isInitialized() &&
0271       m_outputMeasurementParticlesMap.isInitialized()) {
0272     const auto hits = m_inputHits(ctx);
0273
0274     IndexMultimap<ActsFatras::Barcode> outputMap;
0275
0276     for (const auto& [measIdx, hitIdx] : measurementSimHitsMap) {
0277       const auto& hit = hits.nth(hitIdx);
0278       outputMap.emplace(measIdx, hit->particleId());
0279     }
0280
0281     m_outputMeasurementParticlesMap(ctx, std::move(outputMap));
0282   }
0283
0284   // Write the data to the EventStore
0285   m_outputMeasurements(ctx, std::move(measurements));
0286   m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap));
0287
0288   /////////////////////////
0289   // Cluster information //
0290   /////////////////////////
0291
0292   if (m_cfg.outputClusters.empty()) {
0293     return ActsExamples::ProcessCode::SUCCESS;
0294   }
0295
0296   std::vector<ActsExamples::CellData> cellData;
0297
0298   // This allows seamless import of files created with an older version where
0299   // the measurement_id-column is still named hit_id
0300   try {
0301     cellData = readEverything<ActsExamples::CellData>(
0302         m_cfg.inputDir, "cells.csv", {"timestamp"}, ctx.eventNumber);
0303   } catch (std::runtime_error& e) {
0304     // Rethrow exception if it is not about the measurement_id-column
0305     if (std::string(e.what()).find("Missing header column 'measurement_id'") ==
0306         std::string::npos) {
0307       throw;
0308     }
0309
0310     const auto oldCellData = readEverything<ActsExamples::CellDataLegacy>(
0311         m_cfg.inputDir, "cells.csv", {"timestamp"}, ctx.eventNumber);
0312
0313     auto fromLegacy = [](const CellDataLegacy& old) {
0314       return CellData{old.geometry_id, old.hit_id,    old.channel0,
0315                       old.channel1,    old.timestamp, old.value};
0316     };
0317
0318     cellData.resize(oldCellData.size());
0319     std::transform(oldCellData.begin(), oldCellData.end(), cellData.begin(),
0320                    fromLegacy);
0321   }
0322
0323   std::unordered_multimap<std::size_t, ActsExamples::CellData> cellDataMap;
0324   for (const auto& cd : cellData) {
0325     cellDataMap.emplace(cd.measurement_id, cd);
0326   }
0327
0328   auto clusters = makeClusters(cellDataMap, orderedMeasurements.size());
0329   m_outputClusters(ctx, std::move(clusters));
0330
0331   return ActsExamples::ProcessCode::SUCCESS;
0332 }