Csv/src/CsvPlanarClusterReader.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019 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/CsvPlanarClusterReader.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Digitization/DigitizationCell.hpp"
0014 #include "Acts/Digitization/DigitizationSourceLink.hpp"
0015 #include "Acts/Digitization/PlanarModuleCluster.hpp"
0016 #include "Acts/EventData/SourceLink.hpp"
0017 #include "Acts/Geometry/GeometryIdentifier.hpp"
0018 #include "Acts/Geometry/TrackingGeometry.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Utilities/Result.hpp"
0021 #include "ActsExamples/EventData/GeometryContainers.hpp"
0022 #include "ActsExamples/EventData/Index.hpp"
0023 #include "ActsExamples/EventData/SimHit.hpp"
0024 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0025 #include "ActsExamples/Utilities/Paths.hpp"
0026 #include "ActsExamples/Utilities/Range.hpp"
0027 #include "ActsFatras/EventData/Barcode.hpp"
0028 #include "ActsFatras/EventData/Hit.hpp"
0029
0030 #include <algorithm>
0031 #include <array>
0032 #include <cstddef>
0033 #include <cstdint>
0034 #include <iterator>
0035 #include <stdexcept>
0036
0037 #include <dfe/dfe_io_dsv.hpp>
0038
0039 #include "CsvOutputData.hpp"
0040
0041 ActsExamples::CsvPlanarClusterReader::CsvPlanarClusterReader(
0042     const ActsExamples::CsvPlanarClusterReader::Config& config,
0043     Acts::Logging::Level level)
0044     : m_cfg(config),
0045       // TODO check that all files (hits,cells,truth) exists
0046       m_eventsRange(determineEventFilesRange(config.inputDir, "hits.csv")),
0047       m_logger(Acts::getDefaultLogger("CsvPlanarClusterReader", level)) {
0048   if (m_cfg.outputClusters.empty()) {
0049     throw std::invalid_argument("Missing cluster output collection");
0050   }
0051   if (m_cfg.outputHitIds.empty()) {
0052     throw std::invalid_argument("Missing hit id output collection");
0053   }
0054   if (m_cfg.outputMeasurementParticlesMap.empty()) {
0055     throw std::invalid_argument("Missing hit-particles map output collection");
0056   }
0057   if (m_cfg.outputSimHits.empty()) {
0058     throw std::invalid_argument("Missing simulated hits output collection");
0059   }
0060   if (!m_cfg.trackingGeometry) {
0061     throw std::invalid_argument("Missing tracking geometry");
0062   }
0063
0064   m_outputClusters.initialize(m_cfg.outputClusters);
0065   m_outputMeasurementParticlesMap.initialize(
0066       m_cfg.outputMeasurementParticlesMap);
0067   m_outputSimHits.initialize(m_cfg.outputSimHits);
0068   m_outputHitIds.initialize(m_cfg.outputHitIds);
0069 }
0070
0071 std::string ActsExamples::CsvPlanarClusterReader::CsvPlanarClusterReader::name()
0072     const {
0073   return "CsvPlanarClusterReader";
0074 }
0075
0076 std::pair<std::size_t, std::size_t>
0077 ActsExamples::CsvPlanarClusterReader::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 /// Convert separate volume/layer/module id into a single geometry identifier.
0100 inline Acts::GeometryIdentifier extractGeometryId(
0101     const ActsExamples::HitData& data) {
0102   return data.geometry_id;
0103 }
0104
0105 struct CompareGeometryId {
0106   bool operator()(const ActsExamples::HitData& left,
0107                   const ActsExamples::HitData& right) const {
0108     auto leftId = extractGeometryId(left).value();
0109     auto rightId = extractGeometryId(right).value();
0110     return leftId < rightId;
0111   }
0112 };
0113
0114 template <typename Data>
0115 inline std::vector<Data> readEverything(
0116     const std::string& inputDir, const std::string& filename,
0117     const std::vector<std::string>& optionalColumns, std::size_t event) {
0118   std::string path = ActsExamples::perEventFilepath(inputDir, filename, event);
0119   dfe::NamedTupleCsvReader<Data> reader(path, optionalColumns);
0120
0121   std::vector<Data> everything;
0122   Data one;
0123   while (reader.read(one)) {
0124     everything.push_back(one);
0125   }
0126
0127   return everything;
0128 }
0129
0130 std::vector<ActsExamples::HitData> readHitsByGeometryId(
0131     const std::string& inputDir, std::size_t event) {
0132   // geometry_id and t are optional columns
0133   auto hits = readEverything<ActsExamples::HitData>(
0134       inputDir, "hits.csv", {"geometry_id", "t"}, event);
0135   // sort same way they will be sorted in the output container
0136   std::sort(hits.begin(), hits.end(), CompareGeometryId{});
0137   return hits;
0138 }
0139
0140 std::vector<ActsExamples::CellDataLegacy> readCellsByHitId(
0141     const std::string& inputDir, std::size_t event) {
0142   // timestamp is an optional element
0143   auto cells = readEverything<ActsExamples::CellDataLegacy>(
0144       inputDir, "cells.csv", {"timestamp"}, event);
0145   // sort for fast hit id look up
0146   std::sort(cells.begin(), cells.end(), CompareHitId{});
0147   return cells;
0148 }
0149
0150 std::vector<ActsExamples::TruthHitData> readTruthHitsByHitId(
0151     const std::string& inputDir, std::size_t event) {
0152   // define all optional columns
0153   std::vector<std::string> optionalColumns = {
0154       "geometry_id", "tt",      "te",     "deltapx",
0155       "deltapy",     "deltapz", "deltae", "index",
0156   };
0157   auto truths = readEverything<ActsExamples::TruthHitData>(
0158       inputDir, "truth.csv", optionalColumns, event);
0159   // sort for fast hit id look up
0160   std::sort(truths.begin(), truths.end(), CompareHitId{});
0161   return truths;
0162 }
0163
0164 }  // namespace
0165
0166 ActsExamples::ProcessCode ActsExamples::CsvPlanarClusterReader::read(
0167     const ActsExamples::AlgorithmContext& ctx) {
0168   // hit_id in the files is not required to be neither continuous nor
0169   // monotonic. internally, we want continuous indices within [0,#hits)
0170   // to simplify data handling. to be able to perform this mapping we first
0171   // read all data into memory before converting to the internal event data
0172   // types.
0173   auto hits = readHitsByGeometryId(m_cfg.inputDir, ctx.eventNumber);
0174   auto cells = readCellsByHitId(m_cfg.inputDir, ctx.eventNumber);
0175   auto truths = readTruthHitsByHitId(m_cfg.inputDir, ctx.eventNumber);
0176
0177   // prepare containers for the hit data using the framework event data types
0178   GeometryIdMultimap<Acts::PlanarModuleCluster> clusters;
0179   std::vector<uint64_t> hitIds;
0180   IndexMultimap<ActsFatras::Barcode> hitParticlesMap;
0181   SimHitContainer simHits;
0182   clusters.reserve(hits.size());
0183   hitIds.reserve(hits.size());
0184   hitParticlesMap.reserve(truths.size());
0185   simHits.reserve(truths.size());
0186
0187   for (const HitData& hit : hits) {
0188     Acts::GeometryIdentifier geoId = extractGeometryId(hit);
0189
0190     // find associated truth/ simulation hits
0191     std::vector<std::size_t> simHitIndices;
0192     {
0193       auto range = makeRange(std::equal_range(truths.begin(), truths.end(),
0194                                               hit.hit_id, CompareHitId{}));
0195       simHitIndices.reserve(range.size());
0196       for (const auto& truth : range) {
0197         const auto simGeometryId = Acts::GeometryIdentifier(truth.geometry_id);
0198         // TODO validate geo id consistency
0199         const auto simParticleId = ActsFatras::Barcode(truth.particle_id);
0200         const auto simIndex = truth.index;
0201         ActsFatras::Hit::Vector4 simPos4{
0202             truth.tx * Acts::UnitConstants::mm,
0203             truth.ty * Acts::UnitConstants::mm,
0204             truth.tz * Acts::UnitConstants::mm,
0205             truth.tt * Acts::UnitConstants::mm,
0206         };
0207         ActsFatras::Hit::Vector4 simMom4{
0208             truth.tpx * Acts::UnitConstants::GeV,
0209             truth.tpy * Acts::UnitConstants::GeV,
0210             truth.tpz * Acts::UnitConstants::GeV,
0211             truth.te * Acts::UnitConstants::GeV,
0212         };
0213         ActsFatras::Hit::Vector4 simDelta4{
0214             truth.deltapx * Acts::UnitConstants::GeV,
0215             truth.deltapy * Acts::UnitConstants::GeV,
0216             truth.deltapz * Acts::UnitConstants::GeV,
0217             truth.deltae * Acts::UnitConstants::GeV,
0218         };
0219
0220         // the cluster stores indices to the underlying simulation hits. thus
0221         // their position in the container must be stable. the preordering of
0222         // hits by geometry id should ensure that new sim hits are always added
0223         // at the end and previously created ones rest at their existing
0224         // locations.
0225         auto inserted = simHits.emplace_hint(simHits.end(), simGeometryId,
0226                                              simParticleId, simPos4, simMom4,
0227                                              simMom4 + simDelta4, simIndex);
0228         if (std::next(inserted) != simHits.end()) {
0229           ACTS_FATAL("Truth hit sorting broke for input hit id " << hit.hit_id);
0230           return ProcessCode::ABORT;
0231         }
0232         simHitIndices.push_back(simHits.index_of(inserted));
0233       }
0234     }
0235
0236     // find matching pixel cell information
0237     std::vector<Acts::DigitizationCell> digitizationCells;
0238     {
0239       auto range = makeRange(std::equal_range(cells.begin(), cells.end(),
0240                                               hit.hit_id, CompareHitId{}));
0241       for (const auto& c : range) {
0242         digitizationCells.emplace_back(c.channel0, c.channel1, c.value);
0243       }
0244     }
0245
0246     // identify hit surface
0247     const Acts::Surface* surface = m_cfg.trackingGeometry->findSurface(geoId);
0248     if (surface == nullptr) {
0249       ACTS_FATAL("Could not retrieve the surface for hit " << hit);
0250       return ProcessCode::ABORT;
0251     }
0252
0253     // transform global hit coordinates into local coordinates on the surface
0254     Acts::Vector3 pos(hit.x * Acts::UnitConstants::mm,
0255                       hit.y * Acts::UnitConstants::mm,
0256                       hit.z * Acts::UnitConstants::mm);
0257     double time = hit.t * Acts::UnitConstants::mm;
0258     Acts::Vector3 mom(1, 1, 1);  // fake momentum
0259     Acts::Vector2 local(0, 0);
0260     auto lpResult = surface->globalToLocal(ctx.geoContext, pos, mom);
0261     if (!lpResult.ok()) {
0262       ACTS_FATAL("Global to local transformation did not succeed.");
0263       return ProcessCode::ABORT;
0264     }
0265     local = lpResult.value();
0266
0267     // TODO what to use as cluster uncertainty?
0268     Acts::ActsSquareMatrix<3> cov = Acts::ActsSquareMatrix<3>::Identity();
0269     // create the planar cluster
0270     Acts::SourceLink sourceLink{
0271         Acts::DigitizationSourceLink(geoId, std::move(simHitIndices))};
0272     Acts::PlanarModuleCluster cluster(
0273         surface->getSharedPtr(), std::move(sourceLink), cov, local[0], local[1],
0274         time, std::move(digitizationCells));
0275
0276     // due to the previous sorting of the raw hit data by geometry id, new
0277     // clusters should always end up at the end of the container. previous
0278     // elements were not touched; cluster indices remain stable and can
0279     // be used to identify the hit.
0280     auto inserted =
0281         clusters.emplace_hint(clusters.end(), geoId, std::move(cluster));
0282     if (std::next(inserted) != clusters.end()) {
0283       ACTS_FATAL("Something went horribly wrong with the hit sorting");
0284       return ProcessCode::ABORT;
0285     }
0286     auto hitIndex = clusters.index_of(inserted);
0287     auto truthRange = makeRange(std::equal_range(truths.begin(), truths.end(),
0288                                                  hit.hit_id, CompareHitId{}));
0289     for (const auto& truth : truthRange) {
0290       hitParticlesMap.emplace_hint(hitParticlesMap.end(), hitIndex,
0291                                    truth.particle_id);
0292     }
0293
0294     // map internal hit/cluster index back to original, non-monotonic hit id
0295     hitIds.push_back(hit.hit_id);
0296   }
0297
0298   // write the data to the EventStore
0299   m_outputClusters(ctx, std::move(clusters));
0300   m_outputHitIds(ctx, std::move(hitIds));
0301   m_outputMeasurementParticlesMap(ctx, std::move(hitParticlesMap));
0302   m_outputSimHits(ctx, std::move(simHits));
0303
0304   return ActsExamples::ProcessCode::SUCCESS;
0305 }