diff --git a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvMeasurementReader.hpp b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvMeasurementReader.hpp index 1386ff1575c..a6cbcaa99f9 100644 --- a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvMeasurementReader.hpp +++ b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvMeasurementReader.hpp @@ -15,6 +15,8 @@ #include "ActsExamples/EventData/Index.hpp" #include "ActsExamples/EventData/IndexSourceLink.hpp" #include "ActsExamples/EventData/Measurement.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/Framework/DataHandle.hpp" #include "ActsExamples/Framework/IReader.hpp" #include "ActsExamples/Framework/ProcessCode.hpp" @@ -62,6 +64,12 @@ class CsvMeasurementReader final : public IReader { std::string outputSourceLinks; /// Output cluster collection (optional). std::string outputClusters; + + /// Input SimHits for measurment-particle map (optional) + std::string inputSimHits; + /// Output measurement to particle collection (optional) + /// @note Only filled if inputSimHits is given + std::string outputMeasurementParticlesMap; }; /// Construct the cluster reader. @@ -98,6 +106,11 @@ class CsvMeasurementReader final : public IReader { this, "OutputSourceLinks"}; WriteDataHandle m_outputClusters{this, "OutputClusters"}; + + WriteDataHandle> + m_outputMeasurementParticlesMap{this, "OutputMeasurementParticlesMap"}; + + ReadDataHandle m_inputHits{this, "InputHits"}; }; } // namespace ActsExamples diff --git a/Examples/Io/Csv/src/CsvMeasurementReader.cpp b/Examples/Io/Csv/src/CsvMeasurementReader.cpp index 88be4531219..5ab071cb001 100644 --- a/Examples/Io/Csv/src/CsvMeasurementReader.cpp +++ b/Examples/Io/Csv/src/CsvMeasurementReader.cpp @@ -37,7 +37,6 @@ ActsExamples::CsvMeasurementReader::CsvMeasurementReader( const ActsExamples::CsvMeasurementReader::Config& config, Acts::Logging::Level level) : m_cfg(config), - // TODO check that all files (hits,cells,truth) exists m_eventsRange( determineEventFilesRange(m_cfg.inputDir, "measurements.csv")), m_logger(Acts::getDefaultLogger("CsvMeasurementReader", level)) { @@ -49,6 +48,24 @@ ActsExamples::CsvMeasurementReader::CsvMeasurementReader( m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap); m_outputSourceLinks.initialize(m_cfg.outputSourceLinks); m_outputClusters.maybeInitialize(m_cfg.outputClusters); + m_outputMeasurementParticlesMap.maybeInitialize( + m_cfg.outputMeasurementParticlesMap); + m_inputHits.maybeInitialize(m_cfg.inputSimHits); + + // Check if event ranges match (should also catch missing files) + auto checkRange = [&](const std::string& fileStem) { + const auto hitmapRange = determineEventFilesRange(m_cfg.inputDir, fileStem); + if (hitmapRange.first > m_eventsRange.first or + hitmapRange.second < m_eventsRange.second) { + throw std::runtime_error("event range mismatch for 'event**-" + fileStem + + "'"); + } + }; + + checkRange("measurement-simhit-map.csv"); + if (not m_cfg.outputClusters.empty()) { + checkRange("cells.csv"); + } } std::string ActsExamples::CsvMeasurementReader::CsvMeasurementReader::name() @@ -112,14 +129,54 @@ std::vector readMeasurementsByGeometryId( return measurements; } -std::vector readCellsByHitId( - const std::string& inputDir, size_t event) { - // timestamp is an optional element - auto cells = readEverything(inputDir, "cells.csv", - {"timestamp"}, event); - // sort for fast hit id look up - std::sort(cells.begin(), cells.end(), CompareHitId{}); - return cells; +ActsExamples::ClusterContainer makeClusters( + const std::unordered_multimap& + cellDataMap, + std::size_t nMeasurments) { + using namespace ActsExamples; + ClusterContainer clusters; + + for (auto index = 0ul; index < nMeasurments; ++index) { + auto [begin, end] = cellDataMap.equal_range(index); + + // Fill the channels with the iterators + Cluster cluster; + cluster.channels.reserve(std::distance(begin, end)); + + for (auto it = begin; it != end; ++it) { + const auto& cellData = it->second; + ActsFatras::Channelizer::Segment2D dummySegment = {Acts::Vector2::Zero(), + Acts::Vector2::Zero()}; + + ActsFatras::Channelizer::Bin2D bin{ + static_cast(cellData.channel0), + static_cast(cellData.channel1)}; + + cluster.channels.emplace_back(bin, dummySegment, cellData.value); + } + + // update the iterator + + // Compute cluster size + if (not cluster.channels.empty()) { + auto compareX = [](const auto& a, const auto& b) { + return a.bin[0] < b.bin[0]; + }; + auto compareY = [](const auto& a, const auto& b) { + return a.bin[1] < b.bin[1]; + }; + + auto [minX, maxX] = std::minmax_element(cluster.channels.begin(), + cluster.channels.end(), compareX); + auto [minY, maxY] = std::minmax_element(cluster.channels.begin(), + cluster.channels.end(), compareY); + cluster.sizeLoc0 = 1 + maxX->bin[0] - minX->bin[0]; + cluster.sizeLoc1 = 1 + maxY->bin[1] - minY->bin[1]; + } + + clusters.push_back(cluster); + } + return clusters; } } // namespace @@ -136,14 +193,8 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementReader::read( auto measurementData = readMeasurementsByGeometryId(m_cfg.inputDir, ctx.eventNumber); - std::vector cellData = {}; - if (not m_cfg.outputClusters.empty()) { - cellData = readCellsByHitId(m_cfg.inputDir, ctx.eventNumber); - } - // Prepare containers for the hit data using the framework event data types GeometryIdMultimap orderedMeasurements; - ClusterContainer clusters; IndexMultimap measurementSimHitsMap; IndexSourceLinkContainer sourceLinks; // need list here for stable addresses @@ -199,8 +250,8 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementReader::read( // The measurement container is unordered and the index under which // the measurement will be stored is known before adding it. - Index hitIdx = orderedMeasurements.size(); - IndexSourceLink& sourceLink = sourceLinkStorage.emplace_back(geoId, hitIdx); + const Index index = orderedMeasurements.size(); + IndexSourceLink& sourceLink = sourceLinkStorage.emplace_back(geoId, index); auto measurement = createMeasurement(dParameters, sourceLink); // Due to the previous sorting of the raw hit data by geometry id, new @@ -222,13 +273,68 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementReader::read( measurements.emplace_back(std::move(meas)); } + // Generate measurment-particles-map + if (m_inputHits.isInitialized() && + m_outputMeasurementParticlesMap.isInitialized()) { + const auto hits = m_inputHits(ctx); + + IndexMultimap outputMap; + + for (const auto& [measIdx, hitIdx] : measurementSimHitsMap) { + const auto& hit = hits.nth(hitIdx); + outputMap.emplace(measIdx, hit->particleId()); + } + + m_outputMeasurementParticlesMap(ctx, std::move(outputMap)); + } + // Write the data to the EventStore m_outputMeasurements(ctx, std::move(measurements)); m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap)); m_outputSourceLinks(ctx, std::move(sourceLinks)); - if (not clusters.empty()) { - m_outputClusters(ctx, std::move(clusters)); + + ///////////////////////// + // Cluster information // + ///////////////////////// + + if (m_cfg.outputClusters.empty()) { + return ActsExamples::ProcessCode::SUCCESS; } + std::vector cellData; + + // This allows seamless import of files created with a older version where + // the measurment_id-column is still named hit_id + try { + cellData = readEverything( + m_cfg.inputDir, "cells.csv", {"timestamp"}, ctx.eventNumber); + } catch (std::runtime_error& e) { + // Rethrow exception if it is not about the measurement_id-column + if (std::string(e.what()).find("Missing header column 'measurement_id'") == + std::string::npos) { + throw; + } + + const auto oldCellData = readEverything( + m_cfg.inputDir, "cells.csv", {"timestamp"}, ctx.eventNumber); + + auto fromLegacy = [](const CellDataLegacy& old) { + return CellData{old.geometry_id, old.hit_id, old.channel0, + old.channel1, old.timestamp, old.value}; + }; + + cellData.resize(oldCellData.size()); + std::transform(oldCellData.begin(), oldCellData.end(), cellData.begin(), + fromLegacy); + } + + std::unordered_multimap cellDataMap; + for (const auto& cd : cellData) { + cellDataMap.emplace(cd.measurement_id, cd); + } + + auto clusters = makeClusters(cellDataMap, orderedMeasurements.size()); + m_outputClusters(ctx, std::move(clusters)); + return ActsExamples::ProcessCode::SUCCESS; } diff --git a/Examples/Io/Csv/src/CsvMeasurementWriter.cpp b/Examples/Io/Csv/src/CsvMeasurementWriter.cpp index 691a4c908e3..7c89fad5bd9 100644 --- a/Examples/Io/Csv/src/CsvMeasurementWriter.cpp +++ b/Examples/Io/Csv/src/CsvMeasurementWriter.cpp @@ -83,18 +83,18 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementWriter::writeT( MeasurementData meas; CellData cell; - // Will be reused as hit counter + // Will be reused as measurement counter meas.measurement_id = 0; ACTS_VERBOSE("Writing " << measurements.size() << " measurements in this event."); - for (Index hitIdx = 0u; hitIdx < measurements.size(); ++hitIdx) { - const auto& measurement = measurements[hitIdx]; + for (Index measIdx = 0u; measIdx < measurements.size(); ++measIdx) { + const auto& measurement = measurements[measIdx]; - auto simHitIndices = makeRange(measurementSimHitsMap.equal_range(hitIdx)); + auto simHitIndices = makeRange(measurementSimHitsMap.equal_range(measIdx)); for (auto [_, simHitIdx] : simHitIndices) { - writerMeasurementSimHitMap.append({hitIdx, simHitIdx}); + writerMeasurementSimHitMap.append({measIdx, simHitIdx}); } std::visit( @@ -131,9 +131,9 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementWriter::writeT( // CLUSTER / channel information ------------------------------ if (not clusters.empty() && writerCells) { - auto cluster = clusters[hitIdx]; + auto cluster = clusters[measIdx]; cell.geometry_id = meas.geometry_id; - cell.hit_id = meas.measurement_id; + cell.measurement_id = meas.measurement_id; for (auto& c : cluster.channels) { cell.channel0 = c.bin[0]; cell.channel1 = c.bin[1]; diff --git a/Examples/Io/Csv/src/CsvOutputData.hpp b/Examples/Io/Csv/src/CsvOutputData.hpp index cd58b22b1a6..1903306fb0c 100644 --- a/Examples/Io/Csv/src/CsvOutputData.hpp +++ b/Examples/Io/Csv/src/CsvOutputData.hpp @@ -139,9 +139,28 @@ struct MeasurementData { struct CellData { /// Hit surface identifier. uint64_t geometry_id = 0u; - /// Event-unique hit identifier. As defined for the simulated hit above and - /// used to link back to it; same value can appear multiple times for clusters - /// with more than one active cell. + /// Event-unique measurement identifier. As defined for the measurement above + /// and used to link back to it; same value can appear multiple times for + /// clusters with more than one active cell. + uint64_t measurement_id = 0; + /// Digital cell address/ channel + int32_t channel0 = 0, channel1 = 0; + /// Digital cell timestamp. Not available in the TrackML datasets. + float timestamp = 0; + /// (Digital) measured cell value, e.g. amplitude or time-over-threshold. + float value = 0; + + DFE_NAMEDTUPLE(CellData, geometry_id, measurement_id, channel0, channel1, + timestamp, value); +}; + +// uses hit id +struct CellDataLegacy { + /// Hit surface identifier. + uint64_t geometry_id = 0u; + /// Event-unique measurement identifier. As defined for the measurement above + /// and used to link back to it; same value can appear multiple times for + /// clusters with more than one active cell. uint64_t hit_id = 0; /// Digital cell address/ channel int32_t channel0 = 0, channel1 = 0; @@ -150,8 +169,8 @@ struct CellData { /// (Digital) measured cell value, e.g. amplitude or time-over-threshold. float value = 0; - DFE_NAMEDTUPLE(CellData, geometry_id, hit_id, channel0, channel1, timestamp, - value); + DFE_NAMEDTUPLE(CellDataLegacy, geometry_id, hit_id, channel0, channel1, + timestamp, value); }; struct SurfaceData { diff --git a/Examples/Io/Csv/src/CsvPlanarClusterReader.cpp b/Examples/Io/Csv/src/CsvPlanarClusterReader.cpp index 024580e8e6e..f6068fc5f45 100644 --- a/Examples/Io/Csv/src/CsvPlanarClusterReader.cpp +++ b/Examples/Io/Csv/src/CsvPlanarClusterReader.cpp @@ -137,11 +137,11 @@ std::vector readHitsByGeometryId( return hits; } -std::vector readCellsByHitId( +std::vector readCellsByHitId( const std::string& inputDir, size_t event) { // timestamp is an optional element - auto cells = readEverything(inputDir, "cells.csv", - {"timestamp"}, event); + auto cells = readEverything( + inputDir, "cells.csv", {"timestamp"}, event); // sort for fast hit id look up std::sort(cells.begin(), cells.end(), CompareHitId{}); return cells; diff --git a/Examples/Io/Csv/src/CsvPlanarClusterWriter.cpp b/Examples/Io/Csv/src/CsvPlanarClusterWriter.cpp index 19b00eb070f..04ebc841540 100644 --- a/Examples/Io/Csv/src/CsvPlanarClusterWriter.cpp +++ b/Examples/Io/Csv/src/CsvPlanarClusterWriter.cpp @@ -65,13 +65,13 @@ ActsExamples::ProcessCode ActsExamples::CsvPlanarClusterWriter::writeT( perEventFilepath(m_cfg.outputDir, "truth.csv", ctx.eventNumber); dfe::NamedTupleCsvWriter writerHits(pathHits, m_cfg.outputPrecision); - dfe::NamedTupleCsvWriter writerCells(pathCells, - m_cfg.outputPrecision); + dfe::NamedTupleCsvWriter writerCells(pathCells, + m_cfg.outputPrecision); dfe::NamedTupleCsvWriter writerTruth(pathTruth, m_cfg.outputPrecision); HitData hit; - CellData cell; + CellDataLegacy cell; TruthHitData truth; // will be reused as hit counter hit.hit_id = 0; diff --git a/Examples/Python/src/Input.cpp b/Examples/Python/src/Input.cpp index c73bed7f28b..4d365e9cc3b 100644 --- a/Examples/Python/src/Input.cpp +++ b/Examples/Python/src/Input.cpp @@ -55,10 +55,10 @@ void addInput(Context& ctx) { "CsvParticleReader", inputDir, inputStem, outputParticles); - ACTS_PYTHON_DECLARE_READER(ActsExamples::CsvMeasurementReader, mex, - "CsvMeasurementReader", inputDir, - outputMeasurements, outputMeasurementSimHitsMap, - outputSourceLinks, outputClusters); + ACTS_PYTHON_DECLARE_READER( + ActsExamples::CsvMeasurementReader, mex, "CsvMeasurementReader", inputDir, + outputMeasurements, outputMeasurementSimHitsMap, outputSourceLinks, + outputClusters, outputMeasurementParticlesMap, inputSimHits); ACTS_PYTHON_DECLARE_READER(ActsExamples::CsvPlanarClusterReader, mex, "CsvPlanarClusterReader", inputDir, outputClusters, diff --git a/Examples/Python/tests/test_reader.py b/Examples/Python/tests/test_reader.py index e032c3679dd..f745fd75320 100644 --- a/Examples/Python/tests/test_reader.py +++ b/Examples/Python/tests/test_reader.py @@ -173,18 +173,40 @@ def test_csv_meas_reader(tmp_path, fatras, trk_geo, conf_const): out = tmp_path / "csv" out.mkdir() - config = CsvMeasurementWriter.Config( - inputMeasurements=digiAlg.config.outputMeasurements, - inputClusters=digiAlg.config.outputClusters, - inputMeasurementSimHitsMap=digiAlg.config.outputMeasurementSimHitsMap, - outputDir=str(out), + s.addWriter( + CsvMeasurementWriter( + level=acts.logging.INFO, + inputMeasurements=digiAlg.config.outputMeasurements, + inputClusters=digiAlg.config.outputClusters, + inputMeasurementSimHitsMap=digiAlg.config.outputMeasurementSimHitsMap, + outputDir=str(out), + ) ) - s.addWriter(CsvMeasurementWriter(level=acts.logging.INFO, config=config)) + + # Write hits, so we can later construct the measurment-particles-map + s.addWriter( + CsvSimHitWriter( + level=acts.logging.INFO, + inputSimHits=simAlg.config.outputSimHits, + outputDir=str(out), + outputStem="hits", + ) + ) + s.run() # read back in s = Sequencer(numThreads=1) + s.addReader( + CsvSimHitReader( + level=acts.logging.INFO, + outputSimHits=simAlg.config.outputSimHits, + inputDir=str(out), + inputStem="hits", + ) + ) + s.addReader( conf_const( CsvMeasurementReader, @@ -192,13 +214,15 @@ def test_csv_meas_reader(tmp_path, fatras, trk_geo, conf_const): outputMeasurements="measurements", outputMeasurementSimHitsMap="simhitsmap", outputSourceLinks="sourcelinks", + outputMeasurementParticlesMap="meas_ptcl_map", + inputSimHits=simAlg.config.outputSimHits, inputDir=str(out), ) ) algs = [ AssertCollectionExistsAlg(k, f"check_alg_{k}", acts.logging.WARNING) - for k in ("measurements", "simhitsmap", "sourcelinks") + for k in ("measurements", "simhitsmap", "sourcelinks", "meas_ptcl_map") ] for alg in algs: s.addAlgorithm(alg) diff --git a/Tests/UnitTests/Examples/Io/CMakeLists.txt b/Tests/UnitTests/Examples/Io/CMakeLists.txt index f4a992fe62c..2c5d082a6c9 100644 --- a/Tests/UnitTests/Examples/Io/CMakeLists.txt +++ b/Tests/UnitTests/Examples/Io/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory_if(Json ACTS_BUILD_PLUGIN_JSON) add_subdirectory(Root) +add_subdirectory(Csv) diff --git a/Tests/UnitTests/Examples/Io/Csv/CMakeLists.txt b/Tests/UnitTests/Examples/Io/Csv/CMakeLists.txt new file mode 100644 index 00000000000..f6d76c22b41 --- /dev/null +++ b/Tests/UnitTests/Examples/Io/Csv/CMakeLists.txt @@ -0,0 +1,3 @@ +set(unittest_extra_libraries ActsExamplesDigitization ActsExamplesIoCsv) + +add_unittest(MeasurementReaderWriter MeasurementReaderWriterTests.cpp) diff --git a/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp b/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp new file mode 100644 index 00000000000..1c8013130f1 --- /dev/null +++ b/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp @@ -0,0 +1,181 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include + +#include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" +#include "Acts/Tests/CommonHelpers/WhiteBoardUtilities.hpp" +#include "Acts/Utilities/BinUtility.hpp" +#include "Acts/Utilities/Zip.hpp" +#include "ActsExamples/Digitization/DigitizationConfig.hpp" +#include "ActsExamples/Digitization/Smearers.hpp" +#include "ActsExamples/EventData/Cluster.hpp" +#include "ActsExamples/EventData/Measurement.hpp" +#include "ActsExamples/Io/Csv/CsvMeasurementReader.hpp" +#include "ActsExamples/Io/Csv/CsvMeasurementWriter.hpp" + +#include +#include +#include + +using namespace ActsExamples; +using namespace Acts::Test; + +BOOST_AUTO_TEST_CASE(CsvMeasurmentRoundTrip) { + IndexSourceLinkContainer sourceLinksOriginal; + MeasurementContainer measOriginal; + ClusterContainer clusterOriginal; + IndexMultimap mapOriginal; + + //////////////////////////// + // Create some dummy data // + //////////////////////////// + const std::size_t nMeasurements = 3; + Acts::GeometryIdentifier someGeoId{298453}; + + std::mt19937 gen(23); + std::uniform_int_distribution disti(1, 10); + std::uniform_real_distribution distf(0.0, 1.0); + + for (auto i = 0ul; i < nMeasurements; ++i) { + IndexSourceLink sl(someGeoId, static_cast(i)); + sourceLinksOriginal.insert(sl); + + Acts::Vector2 p = Acts::Vector2::Random(); + Acts::SymMatrix2 c = Acts::SymMatrix2::Random(); + + // NOTE this fails: + // auto m = Acts::makeMeasurement(sl, p, c, eBoundLoc0, eBoundTime) + // because we dont support non-consecutive parameters here for now + auto m = Acts::makeMeasurement(Acts::SourceLink{sl}, p, c, Acts::eBoundLoc0, + Acts::eBoundLoc1); + + measOriginal.push_back(m); + + ActsExamples::Cluster cl; + + using Bin2D = ActsFatras::Channelizer::Bin2D; + using Seg2D = ActsFatras::Channelizer::Segment2D; + + // We have two cluster shapes which are displaced randomly + const auto o = disti(gen); + cl.channels.emplace_back( + Bin2D{o + 0, o + 0}, + Seg2D{Acts::Vector2::Random(), Acts::Vector2::Random()}, distf(gen)); + cl.channels.emplace_back( + Bin2D{o + 0, o + 1}, + Seg2D{Acts::Vector2::Random(), Acts::Vector2::Random()}, distf(gen)); + if (distf(gen) < 0.5) { + cl.channels.emplace_back( + Bin2D{o + 0, o + 2}, + Seg2D{Acts::Vector2::Random(), Acts::Vector2::Random()}, distf(gen)); + cl.sizeLoc0 = 1; + cl.sizeLoc1 = 3; + } else { + cl.channels.emplace_back( + Bin2D{o + 1, o + 0}, + Seg2D{Acts::Vector2::Random(), Acts::Vector2::Random()}, distf(gen)); + cl.sizeLoc0 = 2; + cl.sizeLoc1 = 2; + } + + clusterOriginal.push_back(cl); + + // Just generate some random hitid + mapOriginal.insert(std::pair{i, disti(gen)}); + } + + /////////// + // Write // + /////////// + CsvMeasurementWriter::Config writerConfig; + writerConfig.inputClusters = "clusters"; + writerConfig.inputMeasurements = "meas"; + writerConfig.inputMeasurementSimHitsMap = "map"; + writerConfig.outputDir = ""; + + CsvMeasurementWriter writer(writerConfig, Acts::Logging::WARNING); + + auto writeTool = + GenericReadWriteTool<>() + .add(writerConfig.inputMeasurements, measOriginal) + .add(writerConfig.inputClusters, clusterOriginal) + .add(writerConfig.inputMeasurementSimHitsMap, mapOriginal); + + writeTool.write(writer); + + ////////////////// + // Write & Read // + ////////////////// + CsvMeasurementReader::Config readerConfig; + readerConfig.inputDir = writerConfig.outputDir; + readerConfig.outputMeasurements = writerConfig.inputMeasurements; + readerConfig.outputMeasurementSimHitsMap = + writerConfig.inputMeasurementSimHitsMap; + readerConfig.outputClusters = writerConfig.inputClusters; + readerConfig.outputSourceLinks = "sourcelinks"; + + CsvMeasurementReader reader(readerConfig, Acts::Logging::WARNING); + + auto readTool = + writeTool.add(readerConfig.outputSourceLinks, sourceLinksOriginal); + + const auto [measRead, clusterRead, mapRead, sourceLinksRead] = + readTool.read(reader); + + /////////// + // Check // + /////////// + auto checkMeasurementClose = [](const auto &m1, const auto &m2) { + constexpr auto SizeA = std::decay_t::size(); + constexpr auto SizeB = std::decay_t::size(); + if constexpr (SizeA == SizeB) { + CHECK_CLOSE_REL(m1.parameters(), m2.parameters(), 1e-4); + } + }; + + static_assert( + std::is_same_v, decltype(measOriginal)>); + BOOST_REQUIRE(measRead.size() == measOriginal.size()); + for (const auto &[a, b] : Acts::zip(measRead, measOriginal)) { + std::visit(checkMeasurementClose, a, b); + } + + static_assert(std::is_same_v, + decltype(clusterOriginal)>); + BOOST_REQUIRE(clusterRead.size() == clusterOriginal.size()); + for (auto [a, b] : Acts::zip(clusterRead, clusterOriginal)) { + BOOST_REQUIRE(a.sizeLoc0 == b.sizeLoc0); + BOOST_REQUIRE(a.sizeLoc1 == b.sizeLoc1); + + for (const auto &ca : a.channels) { + auto match = [&](const auto &cb) { + return ca.bin == cb.bin && + std::abs(ca.activation - cb.activation) < 1.e-4; + }; + + BOOST_CHECK(std::any_of(b.channels.begin(), b.channels.end(), match)); + } + } + + static_assert( + std::is_same_v, decltype(mapOriginal)>); + BOOST_REQUIRE(mapRead.size() == mapOriginal.size()); + for (const auto &[a, b] : Acts::zip(mapRead, mapOriginal)) { + BOOST_REQUIRE(a == b); + } + + static_assert(std::is_same_v, + decltype(sourceLinksOriginal)>); + BOOST_REQUIRE(sourceLinksRead.size() == sourceLinksOriginal.size()); + for (const auto &[a, b] : Acts::zip(sourceLinksRead, sourceLinksOriginal)) { + BOOST_REQUIRE(a.geometryId() == b.geometryId()); + BOOST_REQUIRE(a.index() == b.index()); + } +}