From 6ff7d492dc0b8cbc65487edbc7c0a6fe7e3bedf3 Mon Sep 17 00:00:00 2001 From: ssdetlab <113530373+ssdetlab@users.noreply.github.com> Date: Fri, 22 Nov 2024 00:59:53 +0200 Subject: [PATCH] feat: Track parameters lookup estimation examples (#3823) The PR adding the simulation based track parameter estimation. The algorithm runs the `Fatras` simulation of the detector setup. A set of grids is imposed onto the user-defined reference layers of the tracking detector (e.g. first tracking layers sensitive surfaces). `CurvilinearTrackParameters` at the vertex and the reference tracking layers are recorded for the simulated particles and stored into the bins of the grids. After the simulation is finished, the contents of the grids' bins are averaged and the grids containing the correspondence between the particles' intersection points at the reference layers, track parameters at the vertex and track parameters at the reference layers are constructed. The constructed grids are stored into the json files. The grids are then readout and used for track parameters estimation in seeding algorithms, e.g. connected to the `PathSeeder`. This PR contains the lookup generation part of the algorithm. When the interfaces, realisation and the general idea are agreed upon, the second part with the validation of the estimated lookups is going to be put up. --- .../EventData/detail/TrackParametersUtils.hpp | 49 ++++ .../TrackParamsLookupAccumulator.hpp | 179 ++++++++++++++ .../Algorithms/TrackFinding/CMakeLists.txt | 1 + .../TrackFinding/ITrackParamsLookupReader.hpp | 27 +++ .../TrackFinding/ITrackParamsLookupWriter.hpp | 27 +++ .../TrackParamsLookupEstimation.hpp | 78 +++++++ .../TrackFinding/TrackParamsLookupTable.hpp | 44 ++++ .../src/TrackParamsLookupEstimation.cpp | 104 +++++++++ .../Io/Json/JsonTrackParamsLookupReader.hpp | 101 ++++++++ .../Io/Json/JsonTrackParamsLookupWriter.hpp | 69 ++++++ Examples/Python/src/Json.cpp | 51 ++++ Examples/Python/src/Output.cpp | 10 + Examples/Python/src/TrackFinding.cpp | 6 + ...elescope_track_params_lookup_generation.py | 111 +++++++++ .../Json/TrackParametersJsonConverter.hpp | 31 +-- .../Core/TrackFinding/CMakeLists.txt | 1 + .../TrackParamsLookupAccumulatorTests.cpp | 221 ++++++++++++++++++ 17 files changed, 1087 insertions(+), 23 deletions(-) create mode 100644 Core/include/Acts/EventData/detail/TrackParametersUtils.hpp create mode 100644 Core/include/Acts/TrackFinding/TrackParamsLookupAccumulator.hpp create mode 100644 Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/ITrackParamsLookupReader.hpp create mode 100644 Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp create mode 100644 Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsLookupEstimation.hpp create mode 100644 Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsLookupTable.hpp create mode 100644 Examples/Algorithms/TrackFinding/src/TrackParamsLookupEstimation.cpp create mode 100644 Examples/Io/Json/include/ActsExamples/Io/Json/JsonTrackParamsLookupReader.hpp create mode 100644 Examples/Io/Json/include/ActsExamples/Io/Json/JsonTrackParamsLookupWriter.hpp create mode 100644 Examples/Scripts/Python/telescope_track_params_lookup_generation.py create mode 100644 Tests/UnitTests/Core/TrackFinding/TrackParamsLookupAccumulatorTests.cpp diff --git a/Core/include/Acts/EventData/detail/TrackParametersUtils.hpp b/Core/include/Acts/EventData/detail/TrackParametersUtils.hpp new file mode 100644 index 00000000000..2efcef0e763 --- /dev/null +++ b/Core/include/Acts/EventData/detail/TrackParametersUtils.hpp @@ -0,0 +1,49 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/EventData/GenericBoundTrackParameters.hpp" +#include "Acts/EventData/TrackParametersConcept.hpp" + +namespace Acts::detail { + +/// @brief Shorthand for Bound or Free track parameters +template +concept isBoundOrFreeTrackParams = + Acts::FreeTrackParametersConcept || + Acts::BoundTrackParametersConcept; + +/// @brief Shorthand for GenericBoundTrackParameters +template +concept isGenericBoundTrackParams = + std::same_as>; + +/// @brief Concept that restricts the type of the +/// accumulation grid cell +template +concept TrackParamsGrid = requires { + typename grid_t::value_type::first_type; + typename grid_t::value_type::second_type; + + requires isBoundOrFreeTrackParams< + typename grid_t::value_type::first_type::element_type>; + requires isBoundOrFreeTrackParams< + typename grid_t::value_type::second_type::element_type>; + + requires requires(typename grid_t::value_type val) { + { + val.first + } -> std::same_as< + std::shared_ptr&>; + { val.second } -> std::same_as; + }; +}; + +} // namespace Acts::detail diff --git a/Core/include/Acts/TrackFinding/TrackParamsLookupAccumulator.hpp b/Core/include/Acts/TrackFinding/TrackParamsLookupAccumulator.hpp new file mode 100644 index 00000000000..a1e1173335d --- /dev/null +++ b/Core/include/Acts/TrackFinding/TrackParamsLookupAccumulator.hpp @@ -0,0 +1,179 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/EventData/detail/TrackParametersUtils.hpp" +#include "Acts/Geometry/GeometryContext.hpp" + +#include +#include +#include +#include +#include + +namespace Acts { + +/// @brief Class to accumulate and average track lookup tables +/// +/// @tparam Grid type for track parameters accumulation +/// +/// This class is used to accumulate track parameters in +/// reference layer grids and average them to create a lookup +/// table for track parameter estimation in seeding +/// +/// @note Geometry context is left to be handled by the user +/// outside of accumulation +template +class TrackParamsLookupAccumulator { + public: + using LookupGrid = grid_t; + using TrackParameters = typename std::pointer_traits< + typename grid_t::value_type::first_type>::element_type; + + /// @brief Constructor + explicit TrackParamsLookupAccumulator(grid_t grid) + : m_grid(std::move(grid)) {} + + /// @brief Add track parameters to the accumulator + /// + /// @param ipTrackParameters Track parameters at the IP + /// @param refTrackParameters Track parameters at the reference layer + /// @param position Local position of the track hit on the reference layer + void addTrack(const TrackParameters& ipTrackParameters, + const TrackParameters& refTrackParameters, + const Vector2& position) { + std::lock_guard lock(m_gridMutex); + + auto bin = m_grid.localBinsFromPosition(position); + + if (m_countGrid[bin] == 0) { + m_grid.atLocalBins(bin).first = + std::make_shared(ipTrackParameters); + m_grid.atLocalBins(bin).second = + std::make_shared(refTrackParameters); + + m_countGrid.at(bin)++; + return; + } + + *m_grid.atLocalBins(bin).first = + addTrackParameters(*m_grid.atLocalBins(bin).first, ipTrackParameters); + *m_grid.atLocalBins(bin).second = + addTrackParameters(*m_grid.atLocalBins(bin).second, refTrackParameters); + m_countGrid.at(bin)++; + } + + /// @brief Finalize the lookup table + /// + /// @return Grid with the bin track parameters averaged + LookupGrid finalizeLookup() { + auto meanTrack = [&](const TrackParameters& track, std::size_t count) { + if constexpr (detail::isGenericBoundTrackParams) { + Acts::GeometryContext gctx; + + auto res = TrackParameters::create( + track.referenceSurface().getSharedPtr(), gctx, + track.fourPosition(gctx) / count, track.momentum().normalized(), + count * track.charge() / track.momentum().norm(), + track.covariance(), track.particleHypothesis()); + + if (!res.ok()) { + throw std::invalid_argument("Bound track grid finalization failed"); + } + return res.value(); + } else { + return TrackParameters(track.fourPosition() / count, + track.momentum().normalized(), + count * track.charge() / track.momentum().norm(), + track.covariance(), track.particleHypothesis()); + } + }; + + for (auto [bin, count] : m_countGrid) { + if (count == 0) { + continue; + } + *m_grid.atLocalBins(bin).first = + meanTrack(*m_grid.atLocalBins(bin).first, count); + *m_grid.atLocalBins(bin).second = + meanTrack(*m_grid.atLocalBins(bin).second, count); + } + + return m_grid; + } + + private: + /// @brief Add two track parameters + /// + /// @param a First track parameter in the sum + /// @param b Second track parameter in the sum + /// + /// @return Sum of track parameters a + b + /// + /// @note Covariances of the track parameters + /// are not added and instead assumed to be + /// generated by the same random process for + /// both a and b, making its averaging redundant + TrackParameters addTrackParameters(const TrackParameters& a, + const TrackParameters& b) { + if (a.particleHypothesis() != b.particleHypothesis()) { + throw std::invalid_argument( + "Cannot accumulate track parameters with different particle " + "hypotheses"); + } + if (a.charge() != b.charge()) { + throw std::invalid_argument( + "Cannot accumulate track parameters with different charges"); + } + if constexpr (detail::isGenericBoundTrackParams) { + if (a.referenceSurface() != b.referenceSurface()) { + throw std::invalid_argument( + "Cannot accumulate bound track parameters with different reference " + "surfaces"); + } + } + + Acts::Vector3 momentum = a.momentum() + b.momentum(); + + // Assume track parameters being i.i.d. + if constexpr (detail::isGenericBoundTrackParams) { + Acts::GeometryContext gctx; + + Acts::Vector4 fourPosition = a.fourPosition(gctx) + b.fourPosition(gctx); + + auto res = TrackParameters::create( + a.referenceSurface().getSharedPtr(), gctx, fourPosition, + momentum.normalized(), a.charge() / momentum.norm(), a.covariance(), + a.particleHypothesis()); + + if (!res.ok()) { + throw std::runtime_error("Invalid bound track parameters"); + } + return res.value(); + } else { + Acts::Vector4 fourPosition = a.fourPosition() + b.fourPosition(); + return TrackParameters(fourPosition, momentum.normalized(), + a.charge() / momentum.norm(), a.covariance(), + a.particleHypothesis()); + } + } + + /// Grids to accumulate IP and reference + /// layer track parameters + LookupGrid m_grid; + + /// Mutex for protecting grid access + std::mutex m_gridMutex; + + /// Map to keep the accumulation count + /// in the occupied grid bins + std::map, std::size_t> m_countGrid; +}; + +} // namespace Acts diff --git a/Examples/Algorithms/TrackFinding/CMakeLists.txt b/Examples/Algorithms/TrackFinding/CMakeLists.txt index 38781494068..33b924e35aa 100644 --- a/Examples/Algorithms/TrackFinding/CMakeLists.txt +++ b/Examples/Algorithms/TrackFinding/CMakeLists.txt @@ -10,6 +10,7 @@ add_library( src/TrackParamsEstimationAlgorithm.cpp src/MuonHoughSeeder.cpp src/GbtsSeedingAlgorithm.cpp + src/TrackParamsLookupEstimation.cpp ) target_include_directories( diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/ITrackParamsLookupReader.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/ITrackParamsLookupReader.hpp new file mode 100644 index 00000000000..5777f90ef18 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/ITrackParamsLookupReader.hpp @@ -0,0 +1,27 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "ActsExamples/TrackFinding/TrackParamsLookupTable.hpp" + +namespace ActsExamples { + +/// @brief Interface for reading track parameter lookup tables +class ITrackParamsLookupReader { + public: + /// Virtual Destructor + virtual ~ITrackParamsLookupReader() = default; + + /// Reader method + /// + /// @param path the path to the file to read + virtual TrackParamsLookup readLookup(const std::string& path) const = 0; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp new file mode 100644 index 00000000000..9cddadae062 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp @@ -0,0 +1,27 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "ActsExamples/TrackFinding/TrackParamsLookupTable.hpp" + +namespace ActsExamples { + +/// @brief Interface for writing track parameter lookup tables +class ITrackParamsLookupWriter { + public: + /// Virtual Destructor + virtual ~ITrackParamsLookupWriter() = default; + + /// Writer method + /// + /// @param lookup track lookup to write + virtual void writeLookup(const TrackParamsLookup& lookup) const = 0; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsLookupEstimation.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsLookupEstimation.hpp new file mode 100644 index 00000000000..e50991cd944 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsLookupEstimation.hpp @@ -0,0 +1,78 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/TrackFinding/TrackParamsLookupAccumulator.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/SimParticle.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" +#include "ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp" + +#include + +namespace ActsExamples { + +/// @brief Algorithm to estimate track parameters lookup tables +/// +/// This algorithm is used to estimate track parameters lookup tables +/// for track parameter estimation in seeding. The algorithm imposes +/// grids onto the reference tracking layers and accumulates track +/// parameters in the grid bins. The track parameters are then averaged +/// to create a lookup table for track parameter estimation in seeding. +class TrackParamsLookupEstimation : public IAlgorithm { + public: + using TrackParamsLookupAccumulator = + Acts::TrackParamsLookupAccumulator; + + /// @brief Nested configuration struct + struct Config { + /// Reference tracking layers + std::unordered_map + refLayers; + /// Binning of the grid to be emposed + /// onto the reference layers + std::pair bins; + /// Input SimHit container + std::string inputHits = "InputHits"; + /// Input SimParticle container + std::string inputParticles = "InputParticles"; + /// Track lookup writers + std::vector> + trackLookupGridWriters{}; + }; + + /// @brief Constructor + TrackParamsLookupEstimation(const Config& config, Acts::Logging::Level level); + + /// @brief The execute method + ProcessCode execute(const AlgorithmContext& ctx) const override; + + ProcessCode finalize() override; + + /// Get readonly access to the config parameters + const Config& config() const { return m_cfg; } + + private: + /// Configuration + Config m_cfg; + + /// Input data handles + ReadDataHandle m_inputParticles{this, + "InputSimParticles"}; + + ReadDataHandle m_inputSimHits{this, "InputSimHits"}; + + /// Accumulators for the track parameters + std::unordered_map> + m_accumulators; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsLookupTable.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsLookupTable.hpp new file mode 100644 index 00000000000..5b9c1fabd83 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsLookupTable.hpp @@ -0,0 +1,44 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/Utilities/Grid.hpp" +#include "Acts/Utilities/GridAxisGenerators.hpp" + +#include +#include + +namespace ActsExamples { + +using TrackParamsLookupPair = + std::pair, + std::shared_ptr>; + +/// @brief Track parameters lookup table axis used +/// in the track estimation algorithm +using TrackParamsLookupAxis = + Acts::Axis; + +/// @brief Track parameters lookup table axis generator +/// used in the track estimation algorithm +using TrackParamsLookupAxisGen = Acts::GridAxisGenerators::EqOpenEqOpen; + +/// @brief Lookup grid for track parameters estimation +/// in a given layer +using TrackParamsLookupGrid = + Acts::Grid; + +/// @brief Lookup table for track parameters estimation +/// in the track estimation algorithm +using TrackParamsLookup = + std::unordered_map; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/src/TrackParamsLookupEstimation.cpp b/Examples/Algorithms/TrackFinding/src/TrackParamsLookupEstimation.cpp new file mode 100644 index 00000000000..4e4b283ddf9 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/src/TrackParamsLookupEstimation.cpp @@ -0,0 +1,104 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include "ActsExamples/TrackFinding/TrackParamsLookupEstimation.hpp" + +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "ActsExamples/Framework/ProcessCode.hpp" + +ActsExamples::TrackParamsLookupEstimation::TrackParamsLookupEstimation( + const Config& config, Acts::Logging::Level level) + : IAlgorithm("TrackParamsLookupEstimation", level), m_cfg(config) { + // Iterate over the reference layers and create + // track parameter accumulators + for (const auto& [geoId, refSurface] : m_cfg.refLayers) { + // Get bounds to construct the accumulator grid + auto bounds = + dynamic_cast(&refSurface->bounds()); + + if (bounds == nullptr) { + throw std::invalid_argument("Only rectangle bounds supported"); + } + if (refSurface->type() != Acts::Surface::SurfaceType::Plane) { + throw std::invalid_argument("Only plane surfaces supported"); + } + + // Initialize the accumulator grid + auto halfX = bounds->halfLengthX(); + auto halfY = bounds->halfLengthY(); + + TrackParamsLookupAxisGen axisGen{ + {-halfX, halfX}, m_cfg.bins.first, {-halfY, halfY}, m_cfg.bins.second}; + + // Each reference layer has its own accumulator + m_accumulators[geoId] = std::make_unique( + TrackParamsLookupGrid(axisGen())); + } + + m_inputParticles.initialize(m_cfg.inputParticles); + m_inputSimHits.initialize(m_cfg.inputHits); +} + +ActsExamples::ProcessCode +ActsExamples::TrackParamsLookupEstimation::finalize() { + // Finiliaze the lookup tables and write them + ActsExamples::TrackParamsLookup lookup; + for (auto& [id, acc] : m_accumulators) { + lookup.insert({id, acc->finalizeLookup()}); + } + for (const auto& writer : m_cfg.trackLookupGridWriters) { + writer->writeLookup(lookup); + } + + return ActsExamples::ProcessCode::SUCCESS; +}; + +ActsExamples::ProcessCode ActsExamples::TrackParamsLookupEstimation::execute( + const ActsExamples::AlgorithmContext& ctx) const { + // Get the particles and hits + const auto& particles = m_inputParticles(ctx); + const auto& hits = m_inputSimHits(ctx); + + // Iterate over the reference layer hits and + // accumulate the track parameters + for (const auto& [geoId, refSurface] : m_cfg.refLayers) { + // Get reference layer hits + auto refLayerHits = hits.equal_range(geoId); + + for (auto hit = refLayerHits.first; hit != refLayerHits.second; ++hit) { + // Get the corresponding particle + const auto& id = hit->particleId(); + const auto& particle = particles.find(id); + + if (particle == particles.end()) { + throw std::invalid_argument("Particle not found"); + } + + // Hit stores the reference layer parameters + auto refLayerPars = Acts::CurvilinearTrackParameters( + hit->fourPosition(), hit->direction(), particle->qOverP(), + std::nullopt, particle->hypothesis()); + + // Particle stores the IP parameters + auto ipPars = Acts::CurvilinearTrackParameters( + particle->fourPosition(), particle->direction(), particle->qOverP(), + std::nullopt, particle->hypothesis()); + + // Get the local position of the hit + auto localPos = refSurface + ->globalToLocal(ctx.geoContext, hit->position(), + Acts::Vector3{0, 1, 0}) + .value(); + + // Add the track parameters to the accumulator grid + m_accumulators.at(geoId)->addTrack(ipPars, refLayerPars, localPos); + } + } + + return ActsExamples::ProcessCode::SUCCESS; +} diff --git a/Examples/Io/Json/include/ActsExamples/Io/Json/JsonTrackParamsLookupReader.hpp b/Examples/Io/Json/include/ActsExamples/Io/Json/JsonTrackParamsLookupReader.hpp new file mode 100644 index 00000000000..94abea8b903 --- /dev/null +++ b/Examples/Io/Json/include/ActsExamples/Io/Json/JsonTrackParamsLookupReader.hpp @@ -0,0 +1,101 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Plugins/Json/GridJsonConverter.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "ActsExamples/TrackFinding/ITrackParamsLookupReader.hpp" + +#include + +#include + +namespace ActsExamples { + +/// @brief Json reader for track parameter lookup tables +/// +/// This reader is used to read track parameter lookup tables +/// from a json file to be later used in track parameter estimation +/// for seeding +class JsonTrackParamsLookupReader final : public ITrackParamsLookupReader { + public: + /// @brief Nested configuration struct + struct Config { + /// Reference tracking layers + std::unordered_map + refLayers; + /// Binning of the grid to be emposed + /// onto the reference layers + std::pair bins; + }; + + explicit JsonTrackParamsLookupReader(const Config& config) : m_cfg(config) {}; + + ~JsonTrackParamsLookupReader() override = default; + + /// @brief Read the lookup from a json file + /// + /// @param path path to the json file + /// + /// @return lookup table for track parameter estimation + TrackParamsLookup readLookup(const std::string& path) const override { + // Read the json file + std::ifstream ifj(path); + nlohmann::json jLookup; + ifj >> jLookup; + + TrackParamsLookup lookup; + // Iterate over the json and deserialize the grids + for (const auto& jGrid : jLookup) { + Acts::GeometryIdentifier id(jGrid["geo_id"]); + + if (!m_cfg.refLayers.contains(id)) { + throw std::invalid_argument("Geometry identifier not found"); + } + + const auto* refSurface = m_cfg.refLayers.at(id); + + // Get bounds to construct the lookup grid + auto bounds = + dynamic_cast(&refSurface->bounds()); + + if (bounds == nullptr) { + throw std::invalid_argument("Only rectangle bounds supported"); + } + + // Axis is not deserilizable, so we need to recreate it + auto halfX = bounds->halfLengthX(); + auto halfY = bounds->halfLengthY(); + + TrackParamsLookupAxisGen axisGen{{-halfX, halfX}, + m_cfg.bins.first, + {-halfY, halfY}, + m_cfg.bins.second}; + + // Deserialize the grid + TrackParamsLookupGrid grid = + Acts::GridJsonConverter::fromJson( + jGrid["grid"], axisGen); + + lookup.try_emplace(id, std::move(grid)); + } + + return lookup; + }; + + /// Readonly access to the config + const Config& config() const { return m_cfg; } + + private: + /// The config of the writer + Config m_cfg; +}; + +} // namespace ActsExamples diff --git a/Examples/Io/Json/include/ActsExamples/Io/Json/JsonTrackParamsLookupWriter.hpp b/Examples/Io/Json/include/ActsExamples/Io/Json/JsonTrackParamsLookupWriter.hpp new file mode 100644 index 00000000000..63a2c083618 --- /dev/null +++ b/Examples/Io/Json/include/ActsExamples/Io/Json/JsonTrackParamsLookupWriter.hpp @@ -0,0 +1,69 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Plugins/Json/GridJsonConverter.hpp" +#include "ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp" + +#include + +#include + +namespace ActsExamples { + +/// @brief Json writer for track parameter lookup tables +/// +/// This writer is used to write track parameter lookup tables +/// to a json file to be later used in track parameter estimation +/// for seeding +class JsonTrackParamsLookupWriter final : public ITrackParamsLookupWriter { + public: + /// @brief Nested configuration struct + struct Config { + /// Output file name + std::string path; + }; + + /// Constructor + /// + /// @param config The configuration struct of the writer + explicit JsonTrackParamsLookupWriter(const Config& config) : m_cfg(config) {}; + + /// Virtual destructor + ~JsonTrackParamsLookupWriter() override = default; + + /// Write out track parameters lookup table + /// + /// @param lookup The lookup to write + void writeLookup(const TrackParamsLookup& lookup) const override { + nlohmann::json jLookup; + + // Iterate over the lookup and serialize the grids + for (const auto& [id, grid] : lookup) { + nlohmann::json jGrid; + jGrid["geo_id"] = id.value(); + jGrid["grid"] = Acts::GridJsonConverter::toJson(grid); + + jLookup.push_back(jGrid); + } + + // Write the json file + std::ofstream ofj(m_cfg.path, std::ios::out); + ofj << std::setw(4) << jLookup << std::endl; + }; + + /// Readonly access to the config + const Config& config() const { return m_cfg; } + + private: + /// The config of the writer + Config m_cfg; +}; + +} // namespace ActsExamples diff --git a/Examples/Python/src/Json.cpp b/Examples/Python/src/Json.cpp index be367612ac2..12baaa9f6d3 100644 --- a/Examples/Python/src/Json.cpp +++ b/Examples/Python/src/Json.cpp @@ -19,6 +19,8 @@ #include "ActsExamples/Io/Json/JsonMaterialWriter.hpp" #include "ActsExamples/Io/Json/JsonSurfacesReader.hpp" #include "ActsExamples/Io/Json/JsonSurfacesWriter.hpp" +#include "ActsExamples/Io/Json/JsonTrackParamsLookupReader.hpp" +#include "ActsExamples/Io/Json/JsonTrackParamsLookupWriter.hpp" #include #include @@ -37,6 +39,11 @@ class IMaterialDecorator; namespace ActsExamples { class IMaterialWriter; class IWriter; + +namespace Experimental { +class ITrackParamsLookupWriter; +} // namespace Experimental + } // namespace ActsExamples namespace py = pybind11; @@ -111,6 +118,50 @@ void addJson(Context& ctx) { ACTS_PYTHON_STRUCT_END(); } + { + using IWriter = ActsExamples::ITrackParamsLookupWriter; + using Writer = ActsExamples::JsonTrackParamsLookupWriter; + using Config = Writer::Config; + + auto cls = py::class_>( + mex, "JsonTrackParamsLookupWriter") + .def(py::init(), py::arg("config")) + .def("writeLookup", &Writer::writeLookup) + .def_property_readonly("config", &Writer::config); + + auto c = py::class_(cls, "Config") + .def(py::init<>()) + .def(py::init(), py::arg("path")); + + ACTS_PYTHON_STRUCT_BEGIN(c, Config); + ACTS_PYTHON_MEMBER(path); + ACTS_PYTHON_STRUCT_END(); + } + + { + using IReader = ActsExamples::ITrackParamsLookupReader; + using Reader = ActsExamples::JsonTrackParamsLookupReader; + using Config = Reader::Config; + + auto cls = py::class_>( + mex, "JsonTrackParamsLookupReader") + .def(py::init(), py::arg("config")) + .def("readLookup", &Reader::readLookup) + .def_property_readonly("config", &Reader::config); + + auto c = py::class_(cls, "Config") + .def(py::init<>()) + .def(py::init, + std::pair>(), + py::arg("refLayers"), py::arg("bins")); + + ACTS_PYTHON_STRUCT_BEGIN(c, Config); + ACTS_PYTHON_MEMBER(refLayers); + ACTS_PYTHON_MEMBER(bins); + ACTS_PYTHON_STRUCT_END(); + } + { auto cls = py::class_ #include @@ -279,6 +281,14 @@ void addOutput(Context& ctx) { py::class_>( mex, "IMaterialWriter"); + py::class_>( + mex, "ITrackParamsLookupWriter"); + + py::class_>( + mex, "ITrackParamsLookupReader"); + { using Writer = ActsExamples::RootMaterialWriter; auto w = py::class_>( diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index ad3ad364ce7..54eb7e76641 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -28,6 +28,7 @@ #include "ActsExamples/TrackFinding/SpacePointMaker.hpp" #include "ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp" #include "ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp" +#include "ActsExamples/TrackFinding/TrackParamsLookupEstimation.hpp" #include #include @@ -293,6 +294,11 @@ void addTrackFinding(Context& ctx) { magneticField, bFieldMin, initialSigmas, initialSigmaPtRel, initialVarInflation, noTimeVarInflation, particleHypothesis); + ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::TrackParamsLookupEstimation, mex, + "TrackParamsLookupEstimation", refLayers, bins, + inputHits, inputParticles, + trackLookupGridWriters); + { using Alg = ActsExamples::TrackFindingAlgorithm; using Config = Alg::Config; diff --git a/Examples/Scripts/Python/telescope_track_params_lookup_generation.py b/Examples/Scripts/Python/telescope_track_params_lookup_generation.py new file mode 100644 index 00000000000..ecdffc20ec3 --- /dev/null +++ b/Examples/Scripts/Python/telescope_track_params_lookup_generation.py @@ -0,0 +1,111 @@ +#!/usr/bin/env python3 + +import argparse + +import acts +import acts.examples +from acts.examples.simulation import ( + addParticleGun, + addFatras, + MomentumConfig, + EtaConfig, + PhiConfig, + ParticleConfig, +) + +u = acts.UnitConstants + + +def estimateLookup(trackingGeometry, numEvents, outputPath): + + # Set up the dipole magnetic field + field = acts.ConstantBField(acts.Vector3(50 * u.T, 0, 0)) + + # Fatras simulation of muons + rnd = acts.examples.RandomNumbers(seed=42) + + s = acts.examples.Sequencer( + events=numEvents, numThreads=1, logLevel=acts.logging.INFO + ) + + vertexGen = acts.examples.GaussianVertexGenerator( + stddev=acts.Vector4(0, 0, 0, 0), mean=acts.Vector4(0, 9, 0, 0) + ) + + addParticleGun( + s=s, + etaConfig=EtaConfig(10.0, 10.0), + phiConfig=PhiConfig(0, 0), + momentumConfig=MomentumConfig(0.5 * u.GeV, 10 * u.GeV), + particleConfig=ParticleConfig(1, acts.PdgParticle.eMuon, False), + multiplicity=1, + rnd=rnd, + vtxGen=vertexGen, + ) + + addFatras( + s, + trackingGeometry, + field, + inputParticles="particles_input", + outputSimHits="sim_hits", + rnd=rnd, + preSelectParticles=None, + ) + + # Set up the track lookup grid writer + jsonWriterConfig = acts.examples.JsonTrackParamsLookupWriter.Config(path=outputPath) + jsonWriter = acts.examples.JsonTrackParamsLookupWriter(jsonWriterConfig) + + # Set up the track estimation algorithm + surfaces = list(trackingGeometry.geoIdSurfaceMap().values()) + refSurface = surfaces[0] + refGeometryId = refSurface.geometryId() + + trackEstConfig = acts.examples.TrackParamsLookupEstimation.Config( + refLayers={refGeometryId: refSurface}, + bins=(1, 1000), + inputHits="sim_hits", + inputParticles="particles_input", + trackLookupGridWriters=[jsonWriter], + ) + trackEstAlg = acts.examples.TrackParamsLookupEstimation( + trackEstConfig, acts.logging.INFO + ) + + s.addAlgorithm(trackEstAlg) + + s.run() + + +if __name__ == "__main__": + p = argparse.ArgumentParser() + + p.add_argument( + "-n", + "--events", + type=int, + default=100000, + help="Number of events for lookup estimation", + ) + p.add_argument( + "-o", + "--output", + type=str, + default="lookup.json", + help="Output lookup file name", + ) + + args = p.parse_args() + + # Initialize the geometry + detector, trackingGeometry, decorators = acts.examples.TelescopeDetector.create( + bounds=[4, 10], + positions=[30, 60, 90], + stereos=[0, 0, 0], + binValue=2, + surfaceType=0, + ) + + # Estimate the lookup + estimateLookup(trackingGeometry, args.events, args.output) diff --git a/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp index ebf7d5c6054..d9670858566 100644 --- a/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp +++ b/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp @@ -8,28 +8,13 @@ #pragma once -#include "Acts/EventData/TrackParameters.hpp" -#include "Acts/Plugins/Json/ActsJson.hpp" +#include "Acts/Definitions/PdgParticle.hpp" +#include "Acts/EventData/GenericBoundTrackParameters.hpp" +#include "Acts/EventData/detail/TrackParametersUtils.hpp" #include "Acts/Plugins/Json/SurfaceJsonConverter.hpp" #include -namespace { - -// Alias to bound adl_serializer specialization -// only to track parameters -template -concept TrackParameters = Acts::FreeTrackParametersConcept || - Acts::BoundTrackParametersConcept; - -// Shorthand for bound track parameters -template -concept IsGenericBound = - std::same_as>; - -} // namespace - namespace Acts { NLOHMANN_JSON_SERIALIZE_ENUM(Acts::PdgParticle, @@ -67,7 +52,7 @@ namespace nlohmann { /// convention is followed. /// /// @tparam parameters_t The track parameters type -template +template struct adl_serializer { /// Covariance matrix type attached to the parameters using CovarianceMatrix = typename parameters_t::CovarianceMatrix; @@ -101,7 +86,7 @@ struct adl_serializer { // Bound track parameters have // reference surface attached // and position takes a geometry context - if constexpr (IsGenericBound) { + if constexpr (Acts::detail::isGenericBoundTrackParams) { Acts::GeometryContext gctx; j["position"] = t.fourPosition(gctx); @@ -152,7 +137,7 @@ struct adl_serializer { // reference surface attached // and constructor is hidden // behind a factory method - if constexpr (IsGenericBound) { + if constexpr (Acts::detail::isGenericBoundTrackParams) { Acts::GeometryContext gctx; auto referenceSurface = Acts::SurfaceJsonConverter::fromJson(j.at("referenceSurface")); @@ -178,7 +163,7 @@ struct adl_serializer { /// convention is followed. /// /// @tparam parameters_t The track parameters type -template +template struct adl_serializer> { using CovarianceMatrix = typename parameters_t::CovarianceMatrix; static void to_json(nlohmann::json& j, @@ -202,7 +187,7 @@ struct adl_serializer> { /// convention is followed. /// /// @tparam parameters_t The track parameters type -template +template struct adl_serializer> { using CovarianceMatrix = typename parameters_t::CovarianceMatrix; static void to_json(nlohmann::json& j, diff --git a/Tests/UnitTests/Core/TrackFinding/CMakeLists.txt b/Tests/UnitTests/Core/TrackFinding/CMakeLists.txt index 523200299ac..bf0c520d672 100644 --- a/Tests/UnitTests/Core/TrackFinding/CMakeLists.txt +++ b/Tests/UnitTests/Core/TrackFinding/CMakeLists.txt @@ -1,2 +1,3 @@ add_unittest(CombinatorialKalmanFilter CombinatorialKalmanFilterTests.cpp) add_unittest(TrackSelector TrackSelectorTests.cpp) +add_unittest(TrackParamsLookupAccumulator TrackParamsLookupAccumulatorTests.cpp) diff --git a/Tests/UnitTests/Core/TrackFinding/TrackParamsLookupAccumulatorTests.cpp b/Tests/UnitTests/Core/TrackFinding/TrackParamsLookupAccumulatorTests.cpp new file mode 100644 index 00000000000..10e12747fce --- /dev/null +++ b/Tests/UnitTests/Core/TrackFinding/TrackParamsLookupAccumulatorTests.cpp @@ -0,0 +1,221 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include +#include +#include + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/EventData/ParticleHypothesis.hpp" +#include "Acts/EventData/TrackParameters.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" +#include "Acts/TrackFinding/TrackParamsLookupAccumulator.hpp" +#include "Acts/Utilities/AxisFwd.hpp" +#include "Acts/Utilities/Grid.hpp" +#include "Acts/Utilities/GridAxisGenerators.hpp" + +#include +#include +#include +#include +#include +#include + +BOOST_AUTO_TEST_SUITE(TrackParamsLookupAccumulator) + +Acts::GeometryContext gctx; + +using Axis = + Acts::Axis; +using AxisGen = Acts::GridAxisGenerators::EqOpenEqOpen; + +using CellBound = std::pair, + std::shared_ptr>; + +using GridBound = Acts::Grid; +using AccBound = Acts::TrackParamsLookupAccumulator; + +using CellCurvilinear = + std::pair, + std::shared_ptr>; + +using GridCurvilinear = Acts::Grid; +using AccCurvilinear = Acts::TrackParamsLookupAccumulator; + +using CellFree = std::pair, + std::shared_ptr>; + +using GridFree = Acts::Grid; +using AccFree = Acts::TrackParamsLookupAccumulator; + +AxisGen axisGen{{-1, 1}, 2, {-1, 1}, 2}; + +BOOST_AUTO_TEST_CASE(Exceptions) { + // Instantiate grid + GridBound grid(axisGen()); + AccBound acc(grid); + + // Create a reference surface for bound parameters + auto transform = Acts::Transform3::Identity(); + auto bounds1 = std::make_shared(1, 1); + auto bounds2 = std::make_shared(2, 2); + + auto surf1 = + Acts::Surface::makeShared(transform, bounds1); + + auto surf2 = + Acts::Surface::makeShared(transform, bounds2); + + // Create parameters to accumulate + Acts::Vector4 pos{1, 2, 0, 4}; + Acts::Vector3 dir{1, 0, 0}; + Acts::ActsScalar P = 1; + + auto hypothesis1 = Acts::ParticleHypothesis::electron(); + auto hypothesis2 = Acts::ParticleHypothesis::muon(); + + auto pars1 = Acts::BoundTrackParameters::create(surf1, gctx, pos, dir, 1. / P, + std::nullopt, hypothesis1) + .value(); + + auto pars2 = Acts::BoundTrackParameters::create(surf2, gctx, pos, dir, 1. / P, + std::nullopt, hypothesis1) + .value(); + + auto pars3 = Acts::BoundTrackParameters::create(surf1, gctx, pos, dir, 1. / P, + std::nullopt, hypothesis2) + .value(); + + auto pars4 = Acts::BoundTrackParameters::create( + surf1, gctx, pos, dir, -1. / P, std::nullopt, hypothesis2) + .value(); + + // Get the point of the grid + auto bin = grid.localBinsFromGlobalBin(2); + auto center = grid.binCenter(bin); + Acts::Vector2 loc{center.at(0), center.at(1)}; + + // Fill in grid + acc.addTrack(pars1, pars1, loc); + + // Different reference surfaces + BOOST_CHECK_THROW(acc.addTrack(pars2, pars2, loc), std::invalid_argument); + + // Different particle hypotheses + BOOST_CHECK_THROW(acc.addTrack(pars3, pars3, loc), std::invalid_argument); + + // Different charges + BOOST_CHECK_THROW(acc.addTrack(pars4, pars4, loc), std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE(Accumulation) { + // Instantiate grids + GridBound gridBound(axisGen()); + AccBound accBound(gridBound); + + GridCurvilinear gridCurvilinear(axisGen()); + AccCurvilinear accCurvilinear(gridCurvilinear); + + GridFree gridFree(axisGen()); + AccFree accFree(gridFree); + + // Create a reference surface for bound parameters + auto transform = Acts::Transform3::Identity(); + auto bounds = std::make_shared(1, 1); + auto surf = Acts::Surface::makeShared(transform, bounds); + + auto hypothesis = Acts::ParticleHypothesis::electron(); + + std::vector avgPoss; + std::vector avgMoms; + Acts::Vector4 pos{1, 2, 0, 4}; + for (std::size_t i = 0; i < gridBound.size(); i++) { + // Create parameters to accumulate + std::array fourPositions = {pos * (i + 1), pos * (i + 2), + pos * (i + 3), pos * (i + 4)}; + + std::array thetas = { + std::numbers::pi / (i + 1), std::numbers::pi / (i + 2), + std::numbers::pi / (i + 3), std::numbers::pi / (i + 4)}; + + std::array phis = { + 2 * std::numbers::pi / (i + 1), 2 * std::numbers::pi / (i + 2), + 2 * std::numbers::pi / (i + 3), 2 * std::numbers::pi / (i + 4)}; + + Acts::ActsScalar P = 1.5 * (i + 1); + + // Get the point of the grid + auto bin = gridBound.localBinsFromGlobalBin(i); + auto center = gridBound.binCenter(bin); + Acts::Vector2 loc{center.at(0), center.at(1)}; + + // Accumulate + Acts::Vector4 avgPos{0, 0, 0, 0}; + Acts::Vector3 avgMom{0, 0, 0}; + for (std::size_t j = 0; j < 4; j++) { + Acts::Vector3 direction{std::sin(thetas.at(j)) * std::cos(phis.at(j)), + std::sin(thetas.at(j)) * std::sin(phis.at(j)), + std::cos(thetas.at(j))}; + + avgPos += fourPositions.at(j); + avgMom += P * direction; + + // Fill in each grid + auto parsBound = Acts::BoundTrackParameters::create( + surf, gctx, fourPositions.at(j), direction, 1. / P, + std::nullopt, hypothesis) + .value(); + + auto parsCurvilinear = Acts::CurvilinearTrackParameters( + fourPositions.at(j), direction, 1. / P, std::nullopt, hypothesis); + + auto parsFree = Acts::FreeTrackParameters( + fourPositions.at(j), direction, 1. / P, std::nullopt, hypothesis); + + accBound.addTrack(parsBound, parsBound, loc); + accCurvilinear.addTrack(parsCurvilinear, parsCurvilinear, loc); + accFree.addTrack(parsFree, parsFree, loc); + } + avgPoss.push_back(avgPos / fourPositions.size()); + avgMoms.push_back(avgMom / fourPositions.size()); + } + + // Finalize and compare + GridBound avgGridBound = accBound.finalizeLookup(); + GridCurvilinear avgGridCurvilinear = accCurvilinear.finalizeLookup(); + GridFree avgGridFree = accFree.finalizeLookup(); + for (std::size_t i = 0; i < avgGridBound.size(); i++) { + auto [ipBound, refBound] = avgGridBound.at(i); + auto [ipCurvilinear, refCurvilinear] = avgGridCurvilinear.at(i); + auto [ipFree, refFree] = avgGridFree.at(i); + + Acts::Vector4 avgPos = avgPoss.at(i); + + Acts::Vector3 avgMom = avgMoms.at(i); + Acts::Vector3 avgDir = avgMom.normalized(); + Acts::ActsScalar avgP = avgMom.norm(); + + CHECK_CLOSE_ABS(ipBound->fourPosition(gctx), avgPos, 1e-3); + CHECK_CLOSE_ABS(ipBound->direction(), avgDir, 1e-3); + CHECK_CLOSE_ABS(ipBound->absoluteMomentum(), avgP, 1e-3); + + CHECK_CLOSE_ABS(ipCurvilinear->fourPosition(), avgPos, 1e-3); + CHECK_CLOSE_ABS(ipCurvilinear->direction(), avgDir, 1e-3); + CHECK_CLOSE_ABS(ipCurvilinear->absoluteMomentum(), avgP, 1e-3); + + CHECK_CLOSE_ABS(ipFree->fourPosition(), avgPos, 1e-3); + CHECK_CLOSE_ABS(ipFree->direction(), avgDir, 1e-3); + CHECK_CLOSE_ABS(ipFree->absoluteMomentum(), avgP, 1e-3); + } +} + +BOOST_AUTO_TEST_SUITE_END()