Csv/src/CsvMeasurementReader.cpp

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