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()