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