diff --git a/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root b/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root index 3452dfba57a..37d7662eadc 100644 Binary files a/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root and b/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root differ diff --git a/Core/include/Acts/Detector/GeometryIdGenerator.hpp b/Core/include/Acts/Detector/GeometryIdGenerator.hpp index 013ffb6c210..925a34ac85e 100644 --- a/Core/include/Acts/Detector/GeometryIdGenerator.hpp +++ b/Core/include/Acts/Detector/GeometryIdGenerator.hpp @@ -77,8 +77,8 @@ class GeometryIdGenerator final : public IGeometryIdGenerator { ~GeometryIdGenerator() override = default; - /// @brief Interface method to generata a geometry id cache - /// @return a geometry id cache decorated in a std::any object + /// @brief Interface method to generate a geometry id cache + /// @return a geometry id cache wrapped in a std::any object IGeometryIdGenerator::GeoIdCache generateCache() const final; /// @brief Method for assigning a geometry id to a detector volume @@ -121,5 +121,102 @@ class GeometryIdGenerator final : public IGeometryIdGenerator { std::unique_ptr m_logger; }; +/// This is a chained tgeometry id generator that will be in seuqnce +/// @tparam generators_t the generators that will be called in sequence +/// +/// @note the generators are expected to be of pointer type +template +class ChainedGeometryIdGenerator : public IGeometryIdGenerator { + public: + struct Cache { + /// The caches + std::array + storage; + }; + + /// The stored generators + std::tuple generators; + + /// Constructor for chained generators_t in a tuple, this will unroll + /// the tuple and call them in sequence + /// + /// @param gens the updators to be called in chain + /// @param mlogger is the logging instance + ChainedGeometryIdGenerator(const std::tuple&& gens, + std::unique_ptr mlogger = + getDefaultLogger("ChainedGeometryIdGenerator", + Logging::INFO)) + : generators(std::move(gens)), m_logger(std::move(mlogger)) {} + + /// @brief Interface method to generate a geometry id cache + /// @return a geometry id cache wrapped in a std::any object + IGeometryIdGenerator::GeoIdCache generateCache() const final { + // Unfold the tuple and add the attachers + Cache cache; + std::size_t it = 0; + std::apply( + [&](auto&&... generator) { + ((cache.storage[it++] = generator->generateCache()), ...); + }, + generators); + return cache; + } + + /// @brief Method for assigning a geometry id to a detector volume + /// + /// @param cache is the cache object for e.g. object counting + /// @param dVolume the detector volume to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& cache, + DetectorVolume& dVolume) const final { + ACTS_VERBOSE("Assigning chained geometry id to volume."); + assign(cache, dVolume); + } + + /// @brief Method for assigning a geometry id to a portal + /// + /// @param cache is the cache object for e.g. object counting + /// @param portal the portal to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& cache, + Portal& portal) const final { + ACTS_VERBOSE("Assigning chained geometry id to portal."); + assign(cache, portal); + } + + /// @brief Method for assigning a geometry id to a surface + /// + /// @param cache is the cache object for e.g. object counting + /// @param surface the surface to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& cache, + Surface& surface) const final { + ACTS_VERBOSE("Assigning chained geometry id to surface."); + assign(cache, surface); + } + + private: + /// @brief Helper to run through the chain of generators + /// + /// @tparam gometry_object_t the geometry object type + /// + /// @param the cache object with the array of sub caches + /// @param object the object to assign the geometry id to + template + void assign(IGeometryIdGenerator::GeoIdCache& cache, + gometry_object_t& object) const { + std::size_t it = 0; + auto& sCache = std::any_cast(cache); + std::apply( + [&](auto&&... generator) { + (generator->assignGeometryId(sCache.storage[it++], object), ...); + }, + generators); + } + + /// Private access method to the logger + const Logger& logger() const { return *m_logger; } + + /// logging instance + std::unique_ptr m_logger; +}; + } // namespace Experimental } // namespace Acts diff --git a/Core/include/Acts/Detector/GeometryIdMapper.hpp b/Core/include/Acts/Detector/GeometryIdMapper.hpp new file mode 100644 index 00000000000..00037f26497 --- /dev/null +++ b/Core/include/Acts/Detector/GeometryIdMapper.hpp @@ -0,0 +1,148 @@ +// 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/. + +#pragma once + +#include "Acts/Detector/DetectorVolume.hpp" +#include "Acts/Detector/Portal.hpp" +#include "Acts/Detector/interface/IGeometryIdGenerator.hpp" +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Surfaces/Surface.hpp" + +#include + +namespace Acts { +namespace Experimental { + +/// @brief This is a mapper of geometry ids, which can be used to +/// assign predefined geometry ids to objects +/// +/// @tparam SourceIdentifier is the type of the source identifier +/// @tparam SourceCapture is the type of the source capture function/struct +/// +/// The source capture function/struct is a callable object/function that +/// can navigate from the provided surface to the source identifier. Usually +/// this would happen via the associated detector element. +/// +/// The only requirement is that the source identifier can be established +/// from the object that receives the target geometry id itself. +template +class GeometryIdMapper final : public IGeometryIdGenerator { + public: + /// @brief Nested config struct + struct Config { + /// The source identifier to target geometry Id map + std::map sourceTargetMap; + /// The source capture function + SourceCapture sourceCapture = SourceCapture(); + }; + + /// @brief Cache object + struct Cache { + unsigned int volumeCounter = 0u; + unsigned int portalCounter = 0u; + unsigned int surfaceCounter = 0u; + }; + + /// @brief Constructor with config + /// + /// @param cfg is the geometry configuration object + /// @param mlogger is the logging instance + GeometryIdMapper(const Config& cfg, + std::unique_ptr mlogger = + getDefaultLogger("GeometryIdMapper", Logging::INFO)) + : m_cfg(cfg), m_logger(std::move(mlogger)) {} + + ~GeometryIdMapper() override = default; + + /// @brief Interface method to generate a geometry id cache + /// @return a geometry id cache wrapped in a std::any object + IGeometryIdGenerator::GeoIdCache generateCache() const final { + return Cache{0u, 0u, 0u}; + } + + /// @brief Method for assigning a geometry id to a detector volume + /// + /// @param cache is the cache object for e.g. object counting + /// @param dVolume the detector volume to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& cache, + DetectorVolume& dVolume) const final { + auto& sCache = std::any_cast(cache); + /// Retrieve the source id for the detector volume + SourceIdentifier vID = m_cfg.sourceCapture(dVolume); + auto source = m_cfg.sourceTargetMap.find(vID); + if (source != m_cfg.sourceTargetMap.end()) { + dVolume.assignGeometryId(source->second); + ACTS_VERBOSE("Assigning geometry id " << source->second << " to volume " + << dVolume.name() << " with id " + << vID); + sCache.volumeCounter++; + } + + // Portals + std::for_each(dVolume.portalPtrs().begin(), dVolume.portalPtrs().end(), + [&](auto& portal) { assignGeometryId(cache, *portal); }); + + // Surfaces + std::for_each(dVolume.surfacePtrs().begin(), dVolume.surfacePtrs().end(), + [&](auto& surface) { assignGeometryId(cache, *surface); }); + + // Sub volumes + std::for_each(dVolume.volumePtrs().begin(), dVolume.volumePtrs().end(), + [&](auto& volume) { assignGeometryId(cache, *volume); }); + } + + /// @brief Method for assigning a geometry id to a portal + /// + /// @param cache is the cache object for e.g. object counting + /// @param portal the portal to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& cache, + Portal& portal) const final { + auto& sCache = std::any_cast(cache); + /// Retrieve the source id for the portal + SourceIdentifier pID = m_cfg.sourceCapture(portal); + auto source = m_cfg.sourceTargetMap.find(pID); + if (source != m_cfg.sourceTargetMap.end()) { + portal.surface().assignGeometryId(source->second); + ACTS_VERBOSE("Assigning geometry id " << source->second << " to portal " + << " with id " << pID); + sCache.portalCounter++; + } + } + + /// @brief Method for assigning a geometry id to a surface + /// + /// @param cache is the cache object for e.g. object counting + /// @param surface the surface to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& cache, + Surface& surface) const final { + auto& sCache = std::any_cast(cache); + /// Retrieve the source id for the surface + SourceIdentifier sID = m_cfg.sourceCapture(surface); + auto source = m_cfg.sourceTargetMap.find(sID); + if (source != m_cfg.sourceTargetMap.end()) { + ACTS_VERBOSE("Assigning geometry id " << source->second << " to surface " + << " with id " << sID); + surface.assignGeometryId(source->second); + sCache.surfaceCounter++; + } + } + + private: + /// Configuration object + Config m_cfg; + + /// Private access method to the logger + const Logger& logger() const { return *m_logger; } + + /// logging instance + std::unique_ptr m_logger; +}; + +} // namespace Experimental +} // namespace Acts diff --git a/Core/include/Acts/Detector/interface/IGeometryIdGenerator.hpp b/Core/include/Acts/Detector/interface/IGeometryIdGenerator.hpp index 3f6c0bd4552..bb161cb9d58 100644 --- a/Core/include/Acts/Detector/interface/IGeometryIdGenerator.hpp +++ b/Core/include/Acts/Detector/interface/IGeometryIdGenerator.hpp @@ -28,8 +28,8 @@ class IGeometryIdGenerator { virtual ~IGeometryIdGenerator() = default; - /// @brief Virtual interface method to generata a geometry id cache - /// @return a geometry id cache decorated in a std::any object + /// @brief Virtual interface method to generate a geometry id cache + /// @return a geometry id cache wrapped in a std::any object virtual GeoIdCache generateCache() const = 0; /// The virtual interface definition for assigning a geometry id to diff --git a/Core/include/Acts/EventData/MultiTrajectory.hpp b/Core/include/Acts/EventData/MultiTrajectory.hpp index babe672e2be..17107013b36 100644 --- a/Core/include/Acts/EventData/MultiTrajectory.hpp +++ b/Core/include/Acts/EventData/MultiTrajectory.hpp @@ -148,29 +148,34 @@ struct IsReadOnlyMultiTrajectory; /// iterating over specific sub-components. template class MultiTrajectory { - public: using Derived = derived_t; static constexpr bool ReadOnly = IsReadOnlyMultiTrajectory::value; // Pull out type alias and re-expose them for ease of use. - // static constexpr unsigned int MeasurementSizeMax = MultiTrajectoryTraits::MeasurementSizeMax; - private: friend class TrackStateProxy; friend class TrackStateProxy; template friend class MultiTrajectory; public: + /// Alias for the const version of a track state proxy, with the same + /// backends as this container using ConstTrackStateProxy = Acts::TrackStateProxy; + + /// Alias for the mutable version of a track state proxy, with the same + /// backends as this container using TrackStateProxy = Acts::TrackStateProxy; + /// The index type of the track state container using IndexType = typename TrackStateProxy::IndexType; + + /// Sentinel value that indicates an invalid index static constexpr IndexType kInvalid = TrackStateProxy::kInvalid; protected: @@ -202,7 +207,16 @@ class MultiTrajectory { } public: + /// @anchor track_state_container_track_access + /// @name MultiTrajectory track state (proxy) access and manipulation + /// + /// These methods allow accessing track states, i.e. adding or retrieving a + /// track state proxy that points at a specific track state in the container. + /// + /// @{ + /// Access a read-only point on the trajectory by index. + /// @note Only available if the MultiTrajectory is not read-only /// @param istate The index to access /// @return Read only proxy to the stored track state ConstTrackStateProxy getTrackState(IndexType istate) const { @@ -210,6 +224,7 @@ class MultiTrajectory { } /// Access a writable point on the trajectory by index. + /// @note Only available if the MultiTrajectory is not read-only /// @param istate The index to access /// @return Read-write proxy to the stored track state template > @@ -217,6 +232,20 @@ class MultiTrajectory { return {*this, istate}; } + /// Add a track state without providing explicit information. Which components + /// of the track state are initialized/allocated can be controlled via @p mask + /// @note Only available if the MultiTrajectory is not read-only + /// @param mask The bitmask that instructs which components to allocate and + /// which to leave invalid + /// @param iprevious index of the previous state, kInvalid if first + /// @return Index of the newly added track state + template > + constexpr IndexType addTrackState( + TrackStatePropMask mask = TrackStatePropMask::All, + IndexType iprevious = kInvalid) { + return self().addTrackState_impl(mask, iprevious); + } + /// Add a track state to the container and return a track state proxy to it /// This effectively calls @c addTrackState and @c getTrackState /// @note Only available if the track state container is not read-only @@ -228,6 +257,12 @@ class MultiTrajectory { return getTrackState(addTrackState(mask, iprevious)); } + /// @} + + /// @anchor track_state_container_iteration + /// @name MultiTrajectory track state iteration + /// @{ + /// Visit all previous states starting at a given endpoint. /// /// @param iendpoint index of the last state @@ -235,6 +270,45 @@ class MultiTrajectory { template void visitBackwards(IndexType iendpoint, F&& callable) const; + /// Apply a function to all previous states starting at a given endpoint. + /// + /// @param iendpoint index of the last state + /// @param callable modifying functor to be called with each point + /// + /// @warning If the trajectory contains multiple components with common + /// points, this can have an impact on the other components. + /// @note Only available if the MultiTrajectory is not read-only + template > + void applyBackwards(IndexType iendpoint, F&& callable) { + static_assert(detail_lt::VisitorConcept, + "Callable needs to satisfy VisitorConcept"); + + if (iendpoint == MultiTrajectoryTraits::kInvalid) { + throw std::runtime_error( + "Cannot apply backwards with kInvalid as endpoint"); + } + + while (true) { + auto ts = getTrackState(iendpoint); + if constexpr (std::is_same_v, + bool>) { + bool proceed = callable(ts); + // this point has no parent and ends the trajectory, or a break was + // requested + if (!proceed || !ts.hasPrevious()) { + break; + } + } else { + callable(ts); + // this point has no parent and ends the trajectory + if (!ts.hasPrevious()) { + break; + } + } + iendpoint = ts.previous(); + } + } + /// Range for the track states from @p iendpoint to the trajectory start /// @param iendpoint Trajectory entry point to start from /// @return Iterator pair to iterate over @@ -251,6 +325,7 @@ class MultiTrajectory { /// Range for the track states from @p iendpoint to the trajectory start, /// i.e from the outside in. + /// @note Only available if the MultiTrajectory is not read-only /// @param iendpoint Trajectory entry point to start from /// @return Iterator pair to iterate over /// @note Mutable version @@ -283,6 +358,7 @@ class MultiTrajectory { /// Range for the track states from @p istartpoint to the trajectory end, /// i.e from inside out + /// @note Only available if the MultiTrajectory is not read-only /// @param istartpoint Trajectory state index for the innermost track /// state to start from /// @return Iterator pair to iterate over @@ -297,79 +373,20 @@ class MultiTrajectory { return range_t{getTrackState(istartpoint)}; } - /// Apply a function to all previous states starting at a given endpoint. - /// - /// @param iendpoint index of the last state - /// @param callable modifying functor to be called with each point - /// - /// @warning If the trajectory contains multiple components with common - /// points, this can have an impact on the other components. - template > - void applyBackwards(IndexType iendpoint, F&& callable) { - static_assert(detail_lt::VisitorConcept, - "Callable needs to satisfy VisitorConcept"); - - if (iendpoint == MultiTrajectoryTraits::kInvalid) { - throw std::runtime_error( - "Cannot apply backwards with kInvalid as endpoint"); - } - - while (true) { - auto ts = getTrackState(iendpoint); - if constexpr (std::is_same_v, - bool>) { - bool proceed = callable(ts); - // this point has no parent and ends the trajectory, or a break was - // requested - if (!proceed || !ts.hasPrevious()) { - break; - } - } else { - callable(ts); - // this point has no parent and ends the trajectory - if (!ts.hasPrevious()) { - break; - } - } - iendpoint = ts.previous(); - } - } - - auto&& convertToReadOnly() const { - auto&& cv = self().convertToReadOnly_impl(); - static_assert( - IsReadOnlyMultiTrajectory::value, - "convertToReadOnly_impl does not return something that reports " - "being ReadOnly."); - return cv; - } - - /// Clear the @c MultiTrajectory. Leaves the underlying storage untouched - template > - constexpr void clear() { - self().clear_impl(); - } - - /// Returns the number of track states contained - constexpr IndexType size() const { return self().size_impl(); } + /// @} - /// Add a track state without providing explicit information. Which components - /// of the track state are initialized/allocated can be controlled via @p mask - /// @param mask The bitmask that instructs which components to allocate and - /// which to leave invalid - /// @param iprevious index of the previous state, kInvalid if first - /// @return Index of the newly added track state - template > - constexpr IndexType addTrackState( - TrackStatePropMask mask = TrackStatePropMask::All, - IndexType iprevious = kInvalid) { - return self().addTrackState_impl(mask, iprevious); - } + /// @anchor track_state_container_columns + /// @name MultiTrajectory column management + /// MultiTrajectory can manage a set of common static columns, and dynamic + /// columns that can be added at runtime. This set of methods allows you to + /// manage the dynamic columns. + /// @{ /// Add a column to the @c MultiTrajectory /// @tparam T Type of the column values to add /// @note This takes a string argument rather than a hashed string to maintain /// compatibility with backends. + /// @note Only available if the MultiTrajectory is not read-only template > constexpr void addColumn(const std::string& key) { self().template addColumn_impl(key); @@ -382,6 +399,19 @@ class MultiTrajectory { return self().hasColumn_impl(key); } + /// @} + + /// Clear the @c MultiTrajectory. Leaves the underlying storage untouched + /// @note Only available if the MultiTrajectory is not read-only + template > + constexpr void clear() { + self().clear_impl(); + } + + /// Returns the number of track states contained + /// @return The number of track states + constexpr IndexType size() const { return self().size_impl(); } + protected: // These are internal helper functions which the @c TrackStateProxy class talks to diff --git a/Core/include/Acts/EventData/ProxyAccessor.hpp b/Core/include/Acts/EventData/ProxyAccessor.hpp new file mode 100644 index 00000000000..a24aa757813 --- /dev/null +++ b/Core/include/Acts/EventData/ProxyAccessor.hpp @@ -0,0 +1,117 @@ +// 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/. + +#pragma once + +#include "Acts/EventData/TrackProxy.hpp" +#include "Acts/EventData/TrackStateProxy.hpp" +#include "Acts/Utilities/Concepts.hpp" +#include "Acts/Utilities/HashedString.hpp" + +#include + +#if defined(__cpp_concepts) +#include +namespace Acts::detail { + +template +concept MutableProxyType = requires(T t, HashedString key) { + requires !T::ReadOnly; + + { + t.template component(key) + } -> std::same_as>; +}; + +template +concept ConstProxyType = requires(T t, HashedString key) { + requires T::ReadOnly; + { t.template component(key) } -> std::same_as; +}; + +template +concept ProxyType = (MutableProxyType || ConstProxyType)&&requires { + typename T::ConstProxyType; + + requires ConstProxyType; +}; + +} // namespace Acts::detail +#endif + +namespace Acts { + +namespace detail { +template +struct associatedConstProxy; + +template +struct associatedConstProxy> { + using type = TrackStateProxy; +}; + +template class holder_t, bool read_only> +struct associatedConstProxy< + TrackProxy> { + using type = TrackProxy; +}; + +} // namespace detail + +/// Utility class that eases accessing dynamic columns in track and track state +/// containers +/// @tparam T the type of the value to access +/// @tparam ReadOnly true if this is a const accessor +template +struct ProxyAccessorBase { + HashedString key; + + /// Create the accessor from an already-hashed string key + /// @param _key the key + constexpr ProxyAccessorBase(HashedString _key) : key{_key} {} + + /// Create the accessor from a string key + /// @param _key the key + ProxyAccessorBase(const std::string& _key) : key{hashString(_key)} {} + + /// Access the stored key on the proxy given as an argument. Mutable version + /// @tparam proxy_t the type of the proxy + /// @param proxy the proxy object to access + /// @return mutable reference to the column behind the key + template > + T& operator()(proxy_t proxy) const { + static_assert(!proxy_t::ReadOnly, + "Cannot get mutable ref for const track proxy"); + return proxy.template component(key); + } + + /// Access the stored key on the proxy given as an argument. Const version + /// @tparam proxy_t the type of the track proxy + /// @param proxy the proxy to access + /// @return const reference to the column behind the key + template > + const T& operator()(proxy_t proxy) const { + if constexpr (proxy_t::ReadOnly) { + return proxy.template component(key); + + } else { + using const_proxy_t = typename proxy_t::ConstProxyType; + const_proxy_t cproxy{proxy}; + return cproxy.template component(key); + } + } +}; + +template +using ProxyAccessor = ProxyAccessorBase; +template +using ConstProxyAccessor = ProxyAccessorBase; +} // namespace Acts diff --git a/Core/include/Acts/EventData/TrackContainer.hpp b/Core/include/Acts/EventData/TrackContainer.hpp index 82d5e3a3ff5..1940c0cc344 100644 --- a/Core/include/Acts/EventData/TrackContainer.hpp +++ b/Core/include/Acts/EventData/TrackContainer.hpp @@ -40,8 +40,12 @@ template class holder_t = detail::RefHolder> class TrackContainer { public: + /// Indicates if this track container is read-only, or if it can be modified static constexpr bool ReadOnly = IsReadOnlyTrackContainer::value; + + /// Indicates if the track state container is read-only, or if it can be + /// modified static constexpr bool TrackStateReadOnly = IsReadOnlyMultiTrajectory::value; @@ -49,11 +53,19 @@ class TrackContainer { "Either both track container and track state container need to " "be readonly or both have to be readwrite"); + /// The index type of the track container, taken from the container backend using IndexType = TrackIndexType; + + /// Sentinel value that indicates an invalid index static constexpr IndexType kInvalid = kTrackIndexInvalid; + /// Alias for the mutable version of a track proxy, with the same backends as + /// this container using TrackProxy = Acts::TrackProxy; + + /// Alias for the const version of a track proxy, with the same backends as + /// this container using ConstTrackProxy = Acts::TrackProxy; @@ -62,6 +74,17 @@ class TrackContainer { friend ConstTrackProxy; #endif + /// @anchor track_container_construction + /// @name TrackContainer construction + /// + /// Constructors for the track container by using a set of backends + /// (track + track state). The container can either take ownership of the + /// backends or just hold references to them. This is driven by the @c + /// holder_t template parameter, where you can also supply a custom holder + /// type. Template deduction is used to try to guess the correct holder type. + /// + /// @{ + /// Constructor from a track container backend and a track state container /// backend /// @param container the track container backend @@ -96,6 +119,16 @@ class TrackContainer { TrackContainer(const track_container_t& container, const traj_t& traj) : m_container{&container}, m_traj{&traj} {} + /// @} + + /// @anchor track_container_track_access + /// @name TrackContainer track (proxy) access and manipulation + /// + /// These methods allow accessing tracks, i.e. adding or retrieving a track + /// proxy that points at a specific track in the container. + /// + /// @{ + /// Get a const track proxy for a track index /// @param itrack the track index in the container /// @return A const track proxy for the index @@ -104,6 +137,7 @@ class TrackContainer { } /// Get a mutable track proxy for a track index + /// @note Only available if the track container is not read-only /// @param itrack the track index in the container /// @return A mutable track proxy for the index template > @@ -111,14 +145,9 @@ class TrackContainer { return {*this, itrack}; } - /// Get the size of the track container - /// @return the sixe - constexpr IndexType size() const { - return m_container->size_impl(); - } - /// Add a track to the container. Note this only creates the logical track and /// allocates memory. You can combine this with @c getTrack to obtain a track proxy + /// @note Only available if the track container is not read-only /// @return the index to the newly added track template > IndexType addTrack() { @@ -137,14 +166,58 @@ class TrackContainer { } /// Remove a track at index @p itrack from the container - /// @note This invalidates all track proxies! + /// @note Only available if the track container is not read-only + /// @note This invalidates track proxies that point to tracks with larger + /// indices than @p itrack! /// @param itrack The index of the track to remove template > void removeTrack(IndexType itrack) { m_container->removeTrack_impl(itrack); } + /// Get a mutable iterator to the first track in the container + /// @note Only available if the track container is not read-only + /// @return a mutable iterator to the first track + template > + auto begin() { + return detail_tc::TrackProxyIterator, + TrackProxy, false>{*this, 0}; + } + + /// Get a past-the-end iterator for this container + /// @note Only available if the track container is not read-only + /// @return a past-the-end iterator + template > + auto end() { + return detail_tc::TrackProxyIterator, + TrackProxy, false>{*this, size()}; + } + + /// Get an const iterator to the first track in the container + /// @return a const iterator to the first track + auto begin() const { + return detail_tc::TrackProxyIterator, + ConstTrackProxy, true>{*this, 0}; + } + + /// Get a past-the-end iterator for this container + /// @return a past-the-end iterator + auto end() const { + return detail_tc::TrackProxyIterator, + ConstTrackProxy, true>{*this, size()}; + } + + /// @} + + /// @anchor track_container_columns + /// @name TrackContainer column management + /// TrackContainer can manage a set of common static columns, and dynamic + /// columns that can be added at runtime. This set of methods allows you to + /// manage the dynamic columns. + /// @{ + /// Add a dymanic column to the track container + /// @note Only available if the track container is not read-only /// @param key the name of the column to be added template > constexpr void addColumn(const std::string& key) { @@ -163,7 +236,29 @@ class TrackContainer { return m_container->hasColumn_impl(key); } + /// Helper function to make this track container match the dynamic columns of + /// another one. This will only work if the track container supports this + /// source, and depends on the implementation details of the dynamic columns + /// of the container + /// @note Only available if the track container is not read-only + /// @tparam other_track_container_t Type of the other track container + /// @param other The other track container + template > + void ensureDynamicColumns(const other_track_container_t& other) { + container().ensureDynamicColumns_impl(other.container()); + } + + /// @} + + /// @anchor track_congtainer_backend_access + /// @name TrackContainer backend access + /// These methods allow accessing the backend of the track container. In most + /// cases, this is not necessary for interacting with the container. + /// @{ + /// Get a mutable reference to the track container backend + /// @note Only available if the track container is not read-only /// @return a mutable reference to the backend template > auto& container() { @@ -177,6 +272,7 @@ class TrackContainer { } /// Get a mutable reference to the track state container backend + /// @note Only available if the track container is not read-only /// @return a mutable reference to the backend template > auto& trackStateContainer() { @@ -185,6 +281,7 @@ class TrackContainer { /// Retrieve the holder of the track state container /// @return The track state container including it's holder + /// @note Only available if the track container is not read-only template > auto& trackStateContainerHolder() { return m_traj; @@ -202,49 +299,16 @@ class TrackContainer { return m_traj; } - /// Get a mutable iterator to the first track in the container - /// @return a mutable iterator to the first track - template > - auto begin() { - return detail_tc::TrackProxyIterator, - TrackProxy, false>{*this, 0}; - } - - /// Get a past-the-end iterator for this container - /// @return a past-the-end iterator - template > - auto end() { - return detail_tc::TrackProxyIterator, - TrackProxy, false>{*this, size()}; - } - - /// Get an const iterator to the first track in the container - /// @return a const iterator to the first track - auto begin() const { - return detail_tc::TrackProxyIterator, - ConstTrackProxy, true>{*this, 0}; - } + /// @} - /// Get a past-the-end iterator for this container - /// @return a past-the-end iterator - auto end() const { - return detail_tc::TrackProxyIterator, - ConstTrackProxy, true>{*this, size()}; - } - - /// Helper function to make this track container match the dynamic columns of - /// another one. This will only work if the track container supports this - /// source, and depends on the implementation details of the dynamic columns - /// of the container - /// @tparam other_track_container_t Type of the other track container - /// @param other The other track container - template > - void ensureDynamicColumns(const other_track_container_t& other) { - container().ensureDynamicColumns_impl(other.container()); + /// Get the size (number of tracks) of the track container + /// @return the sixe + constexpr IndexType size() const { + return m_container->size_impl(); } /// Clear the content of the track container + /// @note Only available if the track container is not read-only template > void clear() { m_container->clear(); @@ -350,51 +414,4 @@ template TrackContainer; -/// Utility class that eases accessing dynamic columns in track containers -/// @tparam T the type of the value to access -/// @tparam ReadOnly true if this is a const accessor -template -struct TrackAccessorBase { - HashedString key; - - /// Create the accessor from an already-hashed string key - /// @param _key the key - TrackAccessorBase(HashedString _key) : key{_key} {} - /// Create the accessor from a string key - /// @param _key the key - TrackAccessorBase(const std::string& _key) : key{hashString(_key)} {} - - /// Access the stored key on the track given as an argument. Mutable version - /// @tparam track_proxy_t the type of the track proxy - /// @param track the track to access - /// @return mutable reference to the column behind the key - template > - T& operator()(track_proxy_t track) const { - static_assert(!track_proxy_t::ReadOnly, - "Cannot get mutable ref for const track proxy"); - return track.template component(key); - } - - /// Access the stored key on the track given as an argument. COnst version - /// @tparam track_proxy_t the type of the track proxy - /// @param track the track to access - /// @return const reference to the column behind the key - template > - const T& operator()(track_proxy_t track) const { - if constexpr (track_proxy_t::ReadOnly) { - return track.template component(key); - } else { - typename track_proxy_t::ConstTrackProxy ctrack{track}; - return ctrack.template component(key); - } - } -}; - -template -using TrackAccessor = TrackAccessorBase; -template -using ConstTrackAccessor = TrackAccessorBase; - } // namespace Acts diff --git a/Core/include/Acts/EventData/TrackProxy.hpp b/Core/include/Acts/EventData/TrackProxy.hpp index 5a36a38afee..43cfd70832f 100644 --- a/Core/include/Acts/EventData/TrackProxy.hpp +++ b/Core/include/Acts/EventData/TrackProxy.hpp @@ -11,6 +11,7 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/ParticleHypothesis.hpp" +#include "Acts/EventData/TrackContainerBackendConcept.hpp" #include "Acts/EventData/TrackParameters.hpp" #include "Acts/EventData/TrackStatePropMask.hpp" #include "Acts/Utilities/Concepts.hpp" @@ -143,7 +144,12 @@ class TrackProxyIterator { } // namespace detail_tc -/// Proxy class representing a single track +/// Proxy class representing a single track. +/// This class provides a **view** into an associated @ref TrackContainer, and +/// has **reference semantics**. You can think of it as a pointer to a vector +/// of tracks, which exposes an object-oriented interface for accessing the +/// track properties. +/// /// @tparam track_container_t the container backend /// @tparam trajectory_t the track state container backend /// @tparam holder_t ownership management class for the backend @@ -152,37 +158,74 @@ template class holder_t, bool read_only = true> class TrackProxy { public: + /// Indicates whether this track proxy is read-only or if it can be modified static constexpr bool ReadOnly = read_only; + + /// The track container backend given as a template parameter using Container = track_container_t; + + /// The track state container backend given as a template parameter using Trajectory = trajectory_t; + /// The index type of the track container + using IndexType = typename Container::IndexType; + + /// Sentinel value that indicates an invalid index + static constexpr IndexType kInvalid = Container::kInvalid; + + /// Alias for the mutable version of this track proxy, with the same backends using MutableTrackProxy = TrackProxy; + + /// Alias for the const version of this track proxy, with the same backends using ConstTrackProxy = TrackProxy; + /// Alias for an associated const track proxy, with the same backends + using ConstProxyType = ConstTrackProxy; + + /// Alias for an associated mutable track state proxy, with the same backends using TrackStateProxy = typename Trajectory::TrackStateProxy; + + /// Alias for an associated const track state proxy, with the same backends using ConstTrackStateProxy = typename Trajectory::ConstTrackStateProxy; + /// Map-type for a bound parameter vector. This has reference semantics, i.e. + /// points at a matrix by an internal pointer. using Parameters = typename detail_lt::Types::CoefficientsMap; + + /// Same as @ref Parameters, but with const semantics using ConstParameters = typename detail_lt::Types::CoefficientsMap; + /// Map-type for a bound covariance. This has reference semantics, i.e. + /// points at a matrix by an internal pointer. using Covariance = typename detail_lt::Types::CovarianceMap; + + /// Same as @ref Covariance, but with const semantics using ConstCovariance = typename detail_lt::Types::CovarianceMap; - using IndexType = typename Container::IndexType; - static constexpr IndexType kInvalid = Container::kInvalid; - #ifndef DOXYGEN friend TrackContainer; friend MutableTrackProxy; friend ConstTrackProxy; + // Track proxies are friends, not food! + template class H, bool R> + friend class TrackProxy; #endif + /// @anchor track_proxy_construct + /// @name Constructors and assignment operator + /// + /// Public constructors and assignment operators for @c TrackProxy only + /// allow construction from another @c TrackProxy. You should generally + /// not have to construct @c TrackProxy manually. + /// + /// @{ + /// Copy constructor from a mutable track proxy. This is always valid, either /// mutable to mutable or mutable to const /// @param other the other track state proxy @@ -198,62 +241,26 @@ class TrackProxy { return *this; } - /// Retrieve a mutable reference to a component - /// @tparam T The type of the component to access - /// @tparam key String key for the component to access - /// @return Mutable reference to the component given by @p key - template > - constexpr T& component() { - return m_container->template component(m_index); - } - - /// Retrieve a mutable reference to a component - /// @tparam T The type of the component to access - /// @param key String key for the component to access - /// @return Mutable reference to the component given by @p key - template > - constexpr T& component(HashedString key) { - return m_container->template component(key, m_index); - } - - /// Retrieve a mutable reference to a component - /// @tparam T The type of the component to access - /// @param key String key for the component to access - /// @note This might hash the @p key at runtime instead of compile-time - /// @return Mutable reference to the component given by @p key - template > - constexpr T& component(std::string_view key) { - return m_container->template component(hashString(key), m_index); - } - - /// Retrieve a const reference to a component - /// @tparam T The type of the component to access - /// @tparam key String key for the component to access - /// @return Const reference to the component given by @p key - template - constexpr const T& component() const { - return m_container->template component(m_index); - } + /// @} - /// Retrieve a const reference to a component - /// @tparam T The type of the component to access - /// @param key String key for the component to access - /// @return Const reference to the component given by @p key - template - constexpr const T& component(HashedString key) const { - return m_container->template component(key, m_index); + /// Equality operator with another track proxy + /// Checks the container identity and the track index + /// @return True if the track proxies refer to the same track + bool operator==(const TrackProxy& other) const { + return &(*m_container) == &(*other.m_container) && m_index == other.m_index; } - /// Retrieve a const reference to a component - /// @tparam T The type of the component to access - /// @param key String key for the component to access - /// @note This might hash the @p key at runtime instead of compile-time - /// @return Const reference to the component given by @p key - template - constexpr const T& component(std::string_view key) const { - return m_container->template component(hashString(key), m_index); - } + /// @anchor track_proxy_props + /// @name TrackProxy properties + /// Methods that give access to the properties of a track represented by + /// @c TrackProxy. + /// + /// Many of these methods come in a @c const and a non-@c const version. The + /// non-@c const version is only available if you have an instance of + /// @c TrackProxy that does not have the @c read_only template parameter set to + /// @c true, even if you hold it as an lvalue. + /// + /// @{ /// Get the tip index, i.e. the entry point into the track state container /// @return the tip index by value @@ -271,6 +278,7 @@ class TrackProxy { /// Get a mutable reference to the tip index, i.e. the entry point into the /// track container + /// @note Only available if the track proxy is not read-only /// @return mutable reference to the tip index template > IndexType& tipIndex() { @@ -280,45 +288,13 @@ class TrackProxy { /// Index of the stem, i.e. the innermost track state of the track. /// This might be invalid, signifying that the track state is not /// forward-linked. + /// @note Only available if the track proxy is not read-only /// @return mutable reference to the stem index template > IndexType& stemIndex() { return component(hashString("stemIndex")); } - /// Return a const track state proxy to the innermost track state - /// @note This is only available, if the track is forward linked - /// @return The innermost track state proxy - auto innermostTrackState() const { - using proxy_t = decltype(m_container->trackStateContainer().getTrackState( - std::declval())); - - IndexType stem = component(hashString("stemIndex")); - if (stem == kInvalid) { - return std::optional{}; - } else { - return std::optional{ - m_container->trackStateContainer().getTrackState(stem)}; - } - } - - /// Return a mutable track state proxy to the innermost track state - /// @note This is only available, if the track is forward linked - /// @return The innermost track state proxy - template > - auto innermostTrackState() { - using proxy_t = decltype(m_container->trackStateContainer().getTrackState( - std::declval())); - - IndexType stem = component(hashString("stemIndex")); - if (stem == kInvalid) { - return std::optional{}; - } else { - return std::optional{ - m_container->trackStateContainer().getTrackState(stem)}; - } - } - /// Get the reference surface of the track (e.g. the perigee) /// @return the reference surface const Surface& referenceSurface() const { @@ -335,7 +311,7 @@ class TrackProxy { } // NOLINTEND(performance-unnecessary-value-param) - /// Return whether a reference surface is associated to this track + /// Returns whether the track has a reference surface or not /// @return whether a surface exists or not bool hasReferenceSurface() const { // @TODO: This could be more efficient @@ -358,6 +334,7 @@ class TrackProxy { /// Get the parameters of the track at the reference surface (e.g. perigee). /// Mutable version + /// @note Only available if the track proxy is not read-only /// @return Proxy vector for the parameters template > Parameters parameters() { @@ -366,6 +343,7 @@ class TrackProxy { /// Get the covariance of the track at the reference surface (e.g. perigee). /// Mutable version + /// @note Only available if the track proxy is not read-only /// @return Proxy matrix for the covariance template > Covariance covariance() { @@ -415,6 +393,7 @@ class TrackProxy { } /// Set a new particle hypothesis for this track + /// @note Only available if the track proxy is not read-only /// @param particleHypothesis The particle hypothesis to set template > void setParticleHypothesis(const ParticleHypothesis& particleHypothesis) { @@ -453,65 +432,6 @@ class TrackProxy { return absoluteMomentum() * direction(); } - /// Get a range over the track states of this track. Return value is - /// compatible with range based for loop. Const version - /// @note This range is from the outside inwards! - /// @return Track state range to iterate over - auto trackStatesReversed() const { - return m_container->reverseTrackStateRange(m_index); - } - - /// Get a range over the track states of this track. Return value is - /// compatible with range based for loop. Mutable version - /// @note This range is from the outside inwards! - /// @return Track state range to iterate over - template > - auto trackStatesReversed() { - return m_container->reverseTrackStateRange(m_index); - } - - /// Get a range over the track states of this track. Return value is - /// compatible with range based for loop. Const version - /// @note This range is from the inside out! - /// @return Track state range to iterate over - auto trackStates() const { - return m_container->forwardTrackStateRange(m_index); - } - - /// Get a range over the track states of this track. Return value is - /// compatible with range based for loop. Mutable version - /// @note This range is from the inside out! - /// @return Track state range to iterate over - template > - auto trackStates() { - return m_container->forwardTrackStateRange(m_index); - } - - /// Forward connect a track, i.e. set indices from the inside out - /// on all track states. - template > - void linkForward() { - IndexType last = kInvalid; - for (auto ts : trackStatesReversed()) { - ts.template component(hashString("next")) = last; - last = ts.index(); - } - stemIndex() = last; - } - - /// Append a track state to this track. This will modify the tip index to - /// point at the newly created track state, which will be directly after the - /// previous track state at tip index. - /// @param mask The allocation prop mask for the new track state - /// @return The newly added track state - template > - auto appendTrackState(TrackStatePropMask mask = TrackStatePropMask::All) { - auto& tsc = m_container->trackStateContainer(); - auto ts = tsc.getTrackState(tsc.addTrackState(mask, tipIndex())); - tipIndex() = ts.index(); - return ts; - } - /// Return the number of track states associated to this track /// @note This is calculated by iterating over the track states which is /// somewhat expensive. Consider caching this value if you need It @@ -529,87 +449,93 @@ class TrackProxy { } /// Return the number of measurements for the track. Const version + /// @note Only available if the track proxy is not read-only /// @return The number of measurements template > unsigned int& nMeasurements() { - return component(hashString("nMeasurements")); + return component(); } /// Return a mutable reference to the number of measurements for the track. /// Mutable version /// @return The number of measurements unsigned int nMeasurements() const { - return component(hashString("nMeasurements")); + return component(); } /// Return a mutable reference to the number of holes for the track. /// Mutable version + /// @note Only available if the track proxy is not read-only /// @return The number of holes template > unsigned int& nHoles() { - return component(hashString("nHoles")); + return component(); } /// Return the number of measurements for the track. Const version /// @return The number of measurements unsigned int nHoles() const { - return component(hashString("nHoles")); + return component(); } /// Return a mutable reference to the number of outliers for the track. /// Mutable version + /// @note Only available if the track proxy is not read-only /// @return The number of outliers template > unsigned int& nOutliers() { - return component(hashString("nOutliers")); + return component(); } /// Return the number of outliers for the track. Const version /// @return The number of outliers unsigned int nOutliers() const { - return component(hashString("nOutliers")); + return component(); } /// Return a mutable reference to the number of shared hits for the track. /// Mutable version + /// @note Only available if the track proxy is not read-only /// @return The number of shared hits template > unsigned int& nSharedHits() { - return component(hashString("nSharedHits")); + return component(); } /// Return the number of shared hits for the track. Const version /// @return The number of shared hits unsigned int nSharedHits() const { - return component(hashString("nSharedHits")); + return component(); } /// Return a mutable reference to the chi squared /// Mutable version + /// @note Only available if the track proxy is not read-only /// @return The chi squared template > float& chi2() { - return component(hashString("chi2")); + return component(); } /// Return the chi squared for the track. Const version /// @return The chi squared float chi2() const { - return component(hashString("chi2")); + return component(); } /// Return a mutable reference to the number of degrees of freedom for the /// track. Mutable version + /// @note Only available if the track proxy is not read-only /// @return The number of degrees of freedom template > unsigned int& nDoF() { - return component(hashString("ndf")); + return component(); } /// Return the number of degrees of freedom for the track. Const version /// @return The number of degrees of freedom unsigned int nDoF() const { - return component(hashString("ndf")); + return component(); } /// Return the index of this track in the track container @@ -619,22 +545,127 @@ class TrackProxy { return m_index; } - /// Return the track parameters at the reference surface - /// @note The parameters are created on the fly - /// @return the track parameters - BoundTrackParameters createParametersAtReference() const { - return BoundTrackParameters(referenceSurface().getSharedPtr(), parameters(), - covariance(), particleHypothesis()); + /// @} + + /// @anchor track_proxy_track_states + /// @name TrackProxy track state access + /// Methods that give access to the track states of a track represented by @c TrackProxy. + /// @{ + + /// Return a const track state proxy to the innermost track state + /// @note This is only available, if the track is forward linked + /// @return The innermost track state proxy + auto innermostTrackState() const { + using proxy_t = decltype(m_container->trackStateContainer().getTrackState( + std::declval())); + + IndexType stem = component(); + if (stem == kInvalid) { + return std::optional{}; + } else { + return std::optional{ + m_container->trackStateContainer().getTrackState(stem)}; + } } - /// Return a reference to the track container backend, mutable version. - /// @return reference to the track container backend + /// Return a mutable track state proxy to the innermost track state + /// @note This is only available, if the track is forward linked + /// @note Only available if the track proxy is not read-only + /// @return The innermost track state proxy template > - auto& container() { - return *m_container; + auto innermostTrackState() { + using proxy_t = decltype(m_container->trackStateContainer().getTrackState( + std::declval())); + + IndexType stem = component(hashString("stemIndex")); + if (stem == kInvalid) { + return std::optional{}; + } else { + return std::optional{ + m_container->trackStateContainer().getTrackState(stem)}; + } + } + + /// Get a range over the track states of this track. Return value is + /// compatible with range based for loop. Const version + /// @note This range is from the outside inwards! + /// @return Track state range to iterate over + auto trackStatesReversed() const { + return m_container->reverseTrackStateRange(m_index); + } + + /// Get a range over the track states of this track. Return value is + /// compatible with range based for loop. Mutable version + /// @note Only available if the track proxy is not read-only + /// @note This range is from the outside inwards! + /// @return Track state range to iterate over + template > + auto trackStatesReversed() { + return m_container->reverseTrackStateRange(m_index); + } + + /// Get a range over the track states of this track. Return value is + /// compatible with range based for loop. This overload returns a const-only + /// track state range, which means you cannot modify the track states obtained + /// in the iteration. + /// @note This range is from the inside out! + /// @warning This access direction is only possible if the track states are + /// **forward-linked**. + /// @return Track state range to iterate over + auto trackStates() const { + return m_container->forwardTrackStateRange(m_index); + } + + /// Get a range over the track states of this track. + /// Return value is compatible with range based for loop. + /// This overload returns a mutable track state range, which means you + /// can modify the track states obtained in the iteration. + /// @note Only available if the track proxy is not read-only + /// @note This range is from the inside out! + /// @warning This access direction is only possible if the track states are + /// **forward-linked**. + /// @return Track state range to iterate over + template > + auto trackStates() { + return m_container->forwardTrackStateRange(m_index); + } + + /// @} + + /// @anchor track_proxy_track_state_manipulation + /// @name TrackProxy track state manipulation + /// Methods that manipulate the track states of a track represented by @c TrackProxy. + /// @{ + + /// Forward connect a track. + /// This means setting indices from the inside out on all track states. + /// @note Only available if the track proxy is not read-only + template > + void linkForward() { + IndexType last = kInvalid; + for (auto ts : trackStatesReversed()) { + ts.template component(hashString("next")) = last; + last = ts.index(); + } + stemIndex() = last; + } + + /// Append a track state to this track. + /// This will modify the tip index to point at the newly created track state, + /// which will be directly after the previous track state at tip index. + /// @note Only available if the track proxy is not read-only + /// @param mask The allocation prop mask for the new track state + /// @return The newly added track state + template > + auto appendTrackState(TrackStatePropMask mask = TrackStatePropMask::All) { + auto& tsc = m_container->trackStateContainer(); + auto ts = tsc.getTrackState(tsc.addTrackState(mask, tipIndex())); + tipIndex() = ts.index(); + return ts; } /// Copy the content of another track proxy into this one + /// @note Only available if the track proxy is not read-only /// @tparam track_proxy_t the other track proxy's type /// @param other The track proxy /// @param copyTrackStates Copy the track state sequence from @p other @@ -679,6 +710,7 @@ class TrackProxy { /// Reverse the ordering of track states for this track /// Afterwards, the previous endpoint of the track state sequence will be the /// "innermost" track state + /// @note Only available if the track proxy is not read-only /// @note This is dangerous with branching track state sequences, as it will break them /// @note This also automatically forward-links the track! /// @param invertJacobians Whether to invert the Jacobians of the track states @@ -714,23 +746,104 @@ class TrackProxy { } } - /// Return a reference to the track container backend, const version. + /// @} + + /// @anchor track_proxy_generic_component + /// @name TrackProxy generic component access + /// Methods that give access to generic components of a track represented by + /// @c TrackProxy. Internally, a compile-time hash of the component name is + /// used to identify which component is being requested. Most of the named + /// methods in @ref track_proxy_props "TrackProxy properties" use these + /// methods to retrieve the actual data. + /// + /// A number of overloads exist, where you can either supply the + /// @ref HashedString @c key as a template parameter or a runtime argument. The + /// former has the advantage of being guaranteed to be evaluated at + /// compile-time. + /// + /// @{ + + /// Retrieve a mutable reference to a component + /// @tparam T The type of the component to access + /// @tparam key String key for the component to access + /// @return Mutable reference to the component given by @p key + template > + constexpr T& component() { + return m_container->template component(m_index); + } + + /// Retrieve a mutable reference to a component + /// @tparam T The type of the component to access + /// @param key String key for the component to access + /// @return Mutable reference to the component given by @p key + template > + constexpr T& component(HashedString key) { + return m_container->template component(key, m_index); + } + + /// Retrieve a mutable reference to a component + /// @tparam T The type of the component to access + /// @param key String key for the component to access + /// @note This might hash the @p key at runtime instead of compile-time + /// @return Mutable reference to the component given by @p key + template > + constexpr T& component(std::string_view key) { + return m_container->template component(hashString(key), m_index); + } + + /// Retrieve a const reference to a component + /// @tparam T The type of the component to access + /// @tparam key String key for the component to access + /// @return Const reference to the component given by @p key + template + constexpr const T& component() const { + return m_container->template component(m_index); + } + + /// Retrieve a const reference to a component + /// @tparam T The type of the component to access + /// @param key String key for the component to access + /// @return Const reference to the component given by @p key + template + constexpr const T& component(HashedString key) const { + return m_container->template component(key, m_index); + } + + /// Retrieve a const reference to a component + /// @tparam T The type of the component to access + /// @param key String key for the component to access + /// @note This might hash the @p key at runtime instead of compile-time + /// @return Const reference to the component given by @p key + template + constexpr const T& component(std::string_view key) const { + return m_container->template component(hashString(key), m_index); + } + + /// @} + + /// Return the track parameters at the reference surface + /// @note The parameters are created on the fly + /// @return the track parameters + BoundTrackParameters createParametersAtReference() const { + return BoundTrackParameters(referenceSurface().getSharedPtr(), parameters(), + covariance(), particleHypothesis()); + } + + /// Return a reference to the track container backend, mutable version. + /// @note Only available if the track proxy is not read-only /// @return reference to the track container backend - const auto& container() const { + template > + auto& container() { return *m_container; } - /// Equality operator with another track proxy - /// Checks the container identity and the track index - /// @return True if the track proxies refer to the same track - bool operator==(const TrackProxy& other) const { - return &(*m_container) == &(*other.m_container) && m_index == other.m_index; + /// Return a reference to the track container backend, const version. + /// @return reference to the track container backend + const auto& container() const { + return *m_container; } - // Track proxies are friends, not food! - template class H, bool R> - friend class TrackProxy; - private: TrackProxy(detail_tc::ConstIf, ReadOnly>& container, diff --git a/Core/include/Acts/EventData/TrackStateProxy.hpp b/Core/include/Acts/EventData/TrackStateProxy.hpp index b82ad298a15..e52ac0ecf64 100644 --- a/Core/include/Acts/EventData/TrackStateProxy.hpp +++ b/Core/include/Acts/EventData/TrackStateProxy.hpp @@ -120,46 +120,91 @@ struct TrackStateTraits { /// /// @tparam SourceLink Type to link back to an original measurement /// @tparam M Maximum number of measurement dimensions -/// @tparam ReadOnly true for read-only access to underlying storage -template +/// @tparam read_only true for read-only access to underlying storage +template class TrackStateProxy { public: + /// Indicates whether this track state proxy is read-only or if it can be + /// modified + static constexpr bool ReadOnly = read_only; + + /// Alias for an associated const track state proxy, with the same backends + using ConstProxyType = TrackStateProxy; + + /// Map-type for a bound parameter vector. This has reference semantics, i.e. + /// points at a matrix by an internal pointer. using Parameters = typename TrackStateTraits::Parameters; - using Covariance = typename TrackStateTraits::Covariance; + + /// Same as @ref Parameters, but with const semantics using ConstParameters = typename TrackStateTraits::Parameters; + + /// Map-type for a bound covariance. This has reference semantics, i.e. + /// points at a matrix by an internal pointer. + using Covariance = typename TrackStateTraits::Covariance; + + /// Same as @ref Covariance, but with const semantics using ConstCovariance = typename TrackStateTraits::Covariance; + /// Map-type for a measurement vector, where the local measurement dimension + /// is variable. template using Measurement = typename TrackStateTraits::Measurement; + + /// Same as @c Measurement, but with const semantics + template + using ConstMeasurement = typename TrackStateTraits::Measurement; + + /// Map-type for a measurement covariance matrix, where the local measurement + /// dimension is variable. template using MeasurementCovariance = typename TrackStateTraits::MeasurementCovariance; - template - using ConstMeasurement = typename TrackStateTraits::Measurement; + + /// Same as @ref MeasurementCovariance, but with const semantics template using ConstMeasurementCovariance = typename TrackStateTraits::MeasurementCovariance; + /// The index type of the track state container using IndexType = TrackIndexType; + + /// Sentinel value that indicates an invalid index static constexpr IndexType kInvalid = kTrackIndexInvalid; - // as opposed to the types above, this is an actual Matrix (rather than a - // map) - // @TODO: Does not copy flags, because this fails: can't have col major row - // vector, but that's required for 1xN projection matrices below. + /// Matrix representing the projector (measurement mapping function) for a + /// measurement. This is not a map type, but an actual matrix. This matrix + /// is always \f$M \times M\f$, even if the local measurement dimension is lower. + /// The actual \f$N\times M\f$ projector is given by the top \f$N\f$ rows. using Projector = typename TrackStateTraits::Projector; + + /// Dynamic variant of the projector matrix + /// @warning Using this type is discouraged, as it has a runtime overhead using EffectiveProjector = Eigen::Matrix::ProjectorFlags, M, eBoundSize>; + /// The track state container backend given as a template parameter using Trajectory = trajectory_t; - // Constructor and assignment operator to construct ReadOnly TrackStateProxy - // from ReadWrite (mutable -> const) + /// @anchor track_state_proxy_construct + /// @name Constructors and assignment operator + /// + /// Public constructors and assignment operators for @c TrackStateProxy only + /// allow construction from another @c TrackStateProxy. You should generally + /// not have to construct @c TrackStateProxy manually. + /// + /// @{ + + /// Constructor and assignment operator to construct TrackStateProxy + /// from mutable + /// @param other The other TrackStateProxy to construct from TrackStateProxy(const TrackStateProxy& other) : m_traj{other.m_traj}, m_istate{other.m_istate} {} + /// Assignment operator to from mutable @c TrackStateProxy + /// @param other The other TrackStateProxy to construct from + /// @return Reference to this TrackStateProxy TrackStateProxy& operator=( const TrackStateProxy& other) { m_traj = other.m_traj; @@ -168,11 +213,58 @@ class TrackStateProxy { return *this; } + /// @} + + /// @anchor track_state_proxy_props + /// @name Track state properties + /// + /// Properties of the track state represented by @c TrackStateProxy. + /// + /// Many of these methods come in a @c const and a non-@c const version. The + /// non-@c const version is only available if you have an instance of + /// @c TrackStateProxy that does not have the @c read_only template parameter set to + /// @c true, even if you hold it as an lvalue. + /// + /// The track states each have an index in the track state container. The + /// sequence of track states is implemented as a one or two-way linked list, + /// which uses indices into the same container. + /// + /// Each track state has a @c previous index, which points at the track state + /// immediately preceding. A track state with a @c previous index of @c + /// kInvalid is the first (innermost) track state in a track or track + /// candidate. This is also referred to as a *stem* at the track level. + /// + /// During track finding and fitting, track states are usually appended to the + /// sequence, populating the @c previous index of the new track state. Combinatorial + /// track finding can produce track states which fork in this way, by having + /// more than one track state with the same @c previous index. + /// + /// The track states have static, optional and dynamic properties. Static + /// properties are always present, and can always be retrieved. Optional + /// components use an extra indirection mechanism that coordinates with the + /// backend to allow both not having the component set, or sharing it with + /// other track states. An example is a branching trajectory from track + /// finding which shares the same predicted parameter vector and associated + /// covariance. + /// + /// Optional components are + /// - predicted parameters and covariance + /// - filtered parameters and covariance + /// - smoothed parameters and covariance + /// - jacobian + /// - calibrated measurement info including projector + /// + /// They can be unset via @ref unset, @ref getMask can be used to check which + /// components are present. The first four are shareable between track + /// states via @ref shareFrom. + /// + /// @{ + /// Index within the trajectory. /// @return the index IndexType index() const { return m_istate; } - /// Return the index of the track state 'previous' in the track sequence + /// Return the index of the track state `previous` in the track sequence /// @return The index of the previous track state. IndexType previous() const { return component(); @@ -180,6 +272,7 @@ class TrackStateProxy { /// Return a mutable reference to the index of the track state 'previous' in /// the track sequence + /// @note Only available if the track state proxy is not read-only /// @return The index of the previous track state. template > IndexType& previous() { @@ -197,191 +290,8 @@ class TrackStateProxy { /// @return The generated mask TrackStatePropMask getMask() const; - /// Share a shareable component within this track state - /// @param shareSource Which component to share from - /// @param shareTarget Which component to share as. This should be different from - /// as @p shareSource, e.g. predicted can be shared as filtered. - /// @note Shareable components are predicted, filtered, smoothed, calibrated, jacobian, - /// or projector. See @c TrackStatePropMask. - template > - void shareFrom(TrackStatePropMask shareSource, - TrackStatePropMask shareTarget) { - shareFrom(*this, shareSource, shareTarget); - } - - /// Share a shareable component from another track state. - /// @param other Track state proxy to share component from - /// @param component Which component to share. - /// @note Shareable components are predicted, filtered, smoothed, calibrated, jacobian, - /// or projector. See @c TrackStatePropMask. - /// @note The track states both need to be stored in the - /// same @c MultiTrajectory instance - template > - void shareFrom(const TrackStateProxy& other, - TrackStatePropMask component) { - shareFrom(other, component, component); - } - - /// Share a shareable component from another track state - /// @param shareSource Which component to share from - /// @param shareTarget Which component to share as. This can be be different from - /// as @p shareSource, e.g. predicted can be shared as filtered. - /// @note Shareable components are predicted, filtered, smoothed, calibrated, jacobian, - /// or projector. See @c TrackStatePropMask. - template > - void shareFrom(const TrackStateProxy& other, - TrackStatePropMask shareSource, - TrackStatePropMask shareTarget) { - assert(m_traj == other.m_traj && - "Cannot share components across MultiTrajectories"); - - assert(ACTS_CHECK_BIT(other.getMask(), shareSource) && - "Source has incompatible allocation"); - - m_traj->self().shareFrom(m_istate, other.m_istate, shareSource, - shareTarget); - } - - /// Copy the contents of another track state proxy into this one - /// @param other The other track state to copy from - /// @param mask An optional mask to determine what to copy from - /// @param onlyAllocated Whether to only copy allocated components - /// @note If the this track state proxy does not have compatible allocations - /// with the source track state proxy, and @p onlyAllocated is false, - /// an exception is thrown. - /// @note The mask parameter will not cause a copy of components that are - /// not allocated in the source track state proxy. - template > - void copyFrom(const track_state_proxy_t& other, - TrackStatePropMask mask = TrackStatePropMask::All, - bool onlyAllocated = true) { - using PM = TrackStatePropMask; - - // @TODO: How to support arbitrary columns here? - - if (onlyAllocated) { - auto dest = getMask(); - auto src = other.getMask() & - mask; // combine what we have with what we want to copy - - if (ACTS_CHECK_BIT(src, PM::Calibrated)) { - // on-demand allocate calibrated - dest |= PM::Calibrated; - } - - if ((static_cast>( - (src ^ dest) & src) != 0 || - dest == TrackStatePropMask::None || - src == TrackStatePropMask::None) && - mask != TrackStatePropMask::None) { - throw std::runtime_error( - "Attempt track state copy with incompatible allocations"); - } - - // we're sure now this has correct allocations, so just copy - if (ACTS_CHECK_BIT(src, PM::Predicted)) { - predicted() = other.predicted(); - predictedCovariance() = other.predictedCovariance(); - } - - if (ACTS_CHECK_BIT(src, PM::Filtered)) { - filtered() = other.filtered(); - filteredCovariance() = other.filteredCovariance(); - } - - if (ACTS_CHECK_BIT(src, PM::Smoothed)) { - smoothed() = other.smoothed(); - smoothedCovariance() = other.smoothedCovariance(); - } - - if (other.hasUncalibratedSourceLink()) { - setUncalibratedSourceLink(other.getUncalibratedSourceLink()); - } - - if (ACTS_CHECK_BIT(src, PM::Jacobian)) { - jacobian() = other.jacobian(); - } - - if (ACTS_CHECK_BIT(src, PM::Calibrated)) { - allocateCalibrated(other.calibratedSize()); - - // workaround for gcc8 bug: - // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86594 - auto* self = this; - visit_measurement(other.calibratedSize(), [&](auto N) { - constexpr int measdim = decltype(N)::value; - self->template calibrated() = - other.template calibrated(); - self->template calibratedCovariance() = - other.template calibratedCovariance(); - }); - - setProjectorBitset(other.projectorBitset()); - } - } else { - if (ACTS_CHECK_BIT(mask, PM::Predicted) && - has() && - other.template has()) { - predicted() = other.predicted(); - predictedCovariance() = other.predictedCovariance(); - } - - if (ACTS_CHECK_BIT(mask, PM::Filtered) && has() && - other.template has()) { - filtered() = other.filtered(); - filteredCovariance() = other.filteredCovariance(); - } - - if (ACTS_CHECK_BIT(mask, PM::Smoothed) && has() && - other.template has()) { - smoothed() = other.smoothed(); - smoothedCovariance() = other.smoothedCovariance(); - } - - if (other.hasUncalibratedSourceLink()) { - setUncalibratedSourceLink(other.getUncalibratedSourceLink()); - } - - if (ACTS_CHECK_BIT(mask, PM::Jacobian) && has() && - other.template has()) { - jacobian() = other.jacobian(); - } - - if (ACTS_CHECK_BIT(mask, PM::Calibrated) && - has() && - other.template has()) { - allocateCalibrated(other.calibratedSize()); - - // workaround for gcc8 bug: - // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86594 - auto* self = this; - visit_measurement(other.calibratedSize(), [&](auto N) { - constexpr int measdim = decltype(N)::value; - self->template calibrated() = - other.template calibrated(); - self->template calibratedCovariance() = - other.template calibratedCovariance(); - }); - - setProjectorBitset(other.projectorBitset()); - } - } - - chi2() = other.chi2(); - pathLength() = other.pathLength(); - typeFlags() = other.typeFlags(); - - if (other.hasReferenceSurface()) { - setReferenceSurface(other.referenceSurface().getSharedPtr()); - } - - m_traj->copyDynamicFrom(m_istate, other.container(), other.index()); - } - /// Unset an optional track state component + /// @note Only available if the track state proxy is not read-only /// @param target The component to unset template > void unset(TrackStatePropMask target) { @@ -414,85 +324,59 @@ class TrackStateProxy { } // NOLINTEND(performance-unnecessary-value-param) - /// Check if a component is set - /// @tparam key Hashed string key to check for - /// @return true if the component exists, false if not - template - constexpr bool has() const { - return m_traj->template has(m_istate); + /// Getter/setter for chi2 value associated with the track state + /// This overload returns a mutable reference, which allows setting a new + /// value directly into the backing store. + /// @note this overload is only enabled in case the proxy is not read-only + /// @return Mutable reference to the chi2 value + template > + float& chi2() { + return component(); } - /// Check if a component is set - /// @param key Hashed string key to check for - /// @return true if the component exists, false if not - constexpr bool has(HashedString key) const { - return m_traj->has(key, m_istate); + /// Getter for the chi2 value associated with the track state. + /// This overload returns a copy of the chi2 value, and thus does not allow + /// modification of the value in the backing storage. + /// @return the chi2 value of the track state + float chi2() const { return component(); } + + /// Getter for the path length associated with the track state. + /// This overloaded is only enabled if not read-only, and returns a mutable + /// reference. + /// @return Mutable reference to the pathlength. + template > + double& pathLength() { + return component(); } - /// Check if a component is set - /// @param key String key to check for - /// @note This might hash the @p key at runtime instead of compile-time - /// @return true if the component exists, false if not - constexpr bool has(std::string_view key) const { - return has(hashString(key)); + /// Getter for the path length. Returns a copy of the path length value. + /// @return The path length of this track state + double pathLength() const { + return component(); } - /// Retrieve a mutable reference to a component - /// @tparam T The type of the component to access - /// @tparam key String key for the component to access - /// @return Mutable reference to the component given by @p key - template > - constexpr T& component() { - return m_traj->template component(m_istate); + /// Getter for the type flags associated with the track state. + /// This overloaded is only enabled if not read-only, and returns a mutable + /// reference. + /// @return reference to the type flags. + template > + TrackStateType typeFlags() { + return TrackStateType{ + component()}; } - /// Retrieve a mutable reference to a component - /// @tparam T The type of the component to access - /// @param key String key for the component to access - /// @return Mutable reference to the component given by @p key - template > - constexpr T& component(HashedString key) { - return m_traj->template component(key, m_istate); + /// Getter for the type flags. Returns a copy of the type flags value. + /// @return The type flags of this track state + ConstTrackStateType typeFlags() const { + return ConstTrackStateType{ + component()}; } - /// Retrieve a mutable reference to a component - /// @tparam T The type of the component to access - /// @param key String key for the component to access - /// @note This might hash the @p key at runtime instead of compile-time - /// @return Mutable reference to the component given by @p key - template > - constexpr T& component(std::string_view key) { - return m_traj->template component(hashString(key), m_istate); - } + /// @} - /// Retrieve a const reference to a component - /// @tparam T The type of the component to access - /// @tparam key String key for the component to access - /// @return Const reference to the component given by @p key - template - constexpr const T& component() const { - return m_traj->template component(m_istate); - } - - /// Retrieve a const reference to a component - /// @tparam T The type of the component to access - /// @param key String key for the component to access - /// @return Const reference to the component given by @p key - template - constexpr const T& component(HashedString key) const { - return m_traj->template component(key, m_istate); - } - - /// Retrieve a const reference to a component - /// @tparam T The type of the component to access - /// @param key String key for the component to access - /// @note This might hash the @p key at runtime instead of compile-time - /// @return Const reference to the component given by @p key - template - constexpr const T& component(std::string_view key) const { - return m_traj->template component(hashString(key), m_istate); - } + /// @anchor track_state_proxy_params + /// @name Track state parameters + /// @{ /// Track parameters vector. This tries to be somewhat smart and return the /// first parameters that are set in this order: predicted -> filtered -> @@ -646,6 +530,38 @@ class TrackStateProxy { /// @return Whether it is set bool hasJacobian() const { return has(); } + /// @} + + /// @anchor track_state_proxy_meas + /// @name Track state measurement properties + /// + /// Properties of the measurement associated with the track state represented. + /// This consists of a vector and an associated square matrix of a measurement + /// dimension which is between one and the size of the track parametrization. + /// The measurement coordinate frame is required to be a strict subset of the + /// bound track parametrization on the local geometry coordinate frame, i.e. + /// using a pure projector matrix to convert from the bound parametrization to + /// the measurement frame is possible. + /// + /// The track state stores the parameter vector and covariance, and the + /// backend is given the possibility to do so in a jagged way, i.e. only + /// storing the number of values needed. This requires calling + /// @ref allocateCalibrated before storing the measurements + /// (even if it might be a no-op). + /// + /// The projector matrix is packed as a bitset, which is converted to a matrix + /// on-demand (and therefore returned by value). + /// + /// A convenience function to assign this from the @ref Measurement class + /// is provided, although it's use is discouraged. + /// + /// The track state also includes a @ref SourceLink which acts as a proxy + /// to the original uncalibrated measurement that the calibrated measurement + /// was derived from. It is set and returned by value, to allow unpacking / + /// repacking by the backend, if needed. + /// + /// @{ + /// Returns the projector (measurement mapping function) for this track /// state. It is derived from the uncalibrated measurement /// @note This function returns the overallocated projector. This means it @@ -661,8 +577,8 @@ class TrackStateProxy { /// Returns the projector (measurement mapping function) for this track /// state. It is derived from the uncalibrated measurement - /// @note This function returns the effective projector. This means it - /// is of dimension NxM, where N is the actual dimension of the + /// @warning This function returns the effective projector. This means it + /// is of dimension \f$N\times M\f$, where \f$N\f$ is the actual dimension of the /// measurement. /// @return The effective projector EffectiveProjector effectiveProjector() const { @@ -702,10 +618,11 @@ class TrackStateProxy { projectorBitset.to_ullong(); } - /// Get the projector bitset, a compressed form of a projection matrix + /// Directly get the projector bitset, a compressed form of a projection + /// matrix /// @note This is mainly to copy explicitly a projector from one state - /// to another. Use the `projector` or `effectiveProjector` method if - /// you want to access the matrix. + /// to another. Use the `projector` or `effectiveProjector` method if + /// you want to access the matrix. /// @return The projector bitset ProjectorBitset projectorBitset() const { assert(has()); @@ -716,7 +633,8 @@ class TrackStateProxy { /// @param proj The projector bitset /// /// @note This is mainly to copy explicitly a projector from one state - /// to another. If you have a projection matrix, set it with `setProjector`. + /// to another. If you have a projection matrix, set it with + /// `setProjector`. template > void setProjectorBitset(ProjectorBitset proj) { assert(has()); @@ -765,20 +683,18 @@ class TrackStateProxy { return m_traj->self().template measurement(m_istate); } - /// Full calibrated measurement covariance matrix. The effective covariance - /// is located in the top left corner, everything else is zeroed. + /// Const full calibrated measurement covariance matrix. The effective + /// covariance is located in the top left corner, everything else is zeroed. /// @return The measurement covariance matrix - /// @note Const version template ConstMeasurementCovariance calibratedCovariance() const { assert(has()); return m_traj->self().template measurementCovariance(m_istate); } - /// Full calibrated measurement covariance matrix. The effective covariance - /// is located in the top left corner, everything else is zeroed. + /// Mutable full calibrated measurement covariance matrix. The effective + /// covariance is located in the top left corner, everything else is zeroed. /// @return The measurement covariance matrix - /// @note Mutable version template > MeasurementCovariance calibratedCovariance() { @@ -786,9 +702,9 @@ class TrackStateProxy { return m_traj->self().template measurementCovariance(m_istate); } - /// Dynamic measurement vector with only the valid dimensions. + /// Mutable dynamic measurement vector with only the valid dimensions. + /// @warning The dynamic vector has a runtime overhead! /// @return The effective calibrated measurement vector - /// @note Mutable version template > auto effectiveCalibrated() { // repackage the data pointer to a dynamic map type @@ -803,9 +719,9 @@ class TrackStateProxy { }); } - /// Dynamic measurement vector with only the valid dimensions. + /// Const dynamic measurement vector with only the valid dimensions. + /// @warning The dynamic matrix has a runtime overhead! /// @return The effective calibrated measurement vector - /// @note Const version auto effectiveCalibrated() const { // repackage the data pointer to a dynamic map type // workaround for gcc8 bug: @@ -819,9 +735,10 @@ class TrackStateProxy { }); } - /// Dynamic measurement covariance matrix with only the valid dimensions. + /// Mutable dynamic measurement covariance matrix with only the valid + /// dimensions. + /// @warning The dynamic matrix has a runtime overhead! /// @return The effective calibrated covariance matrix - /// @note Mutable version template > auto effectiveCalibratedCovariance() { // repackage the data pointer to a dynamic map type @@ -836,9 +753,10 @@ class TrackStateProxy { }); } - /// Dynamic measurement covariance matrix with only the valid dimensions. + /// Const dynamic measurement covariance matrix with only the valid + /// dimensions. + /// @warning The dynamic matrix has a runtime overhead! /// @return The effective calibrated covariance matrix - /// @note Const version auto effectiveCalibratedCovariance() const { // repackage the data pointer to a dynamic map type // workaround for gcc8 bug: @@ -853,20 +771,24 @@ class TrackStateProxy { } /// Return the (dynamic) number of dimensions stored for this measurement. - /// @note The underlying storage is overallocated to MeasurementSizeMax - /// regardless of this value + /// @note Depending on the backend, this size is used to determine the + /// memory range of the measurement vector and covariance. /// @return The number of dimensions IndexType calibratedSize() const { return m_traj->calibratedSize(m_istate); } - /// Overwrite existing measurement data. + /// Set the calibrated measurement of this track state from a measurement + /// object. This is a convenience function to set the calibrated measurement + /// from the @c Acts::Measurement object. In general, you should implement this + /// functionality specifically for an (experiment-specific) uncalibrated + /// measurement EDM. /// /// @tparam kMeasurementSize Size of the calibrated measurement /// @param meas The measurement object to set /// - /// @note This assumes this TrackState stores it's own calibrated - /// measurement. **If storage is shared with another TrackState, both will - /// be overwritten!**. Also assumes none of the calibrated components is - /// *invalid* (i.e. unset) for this TrackState. + /// @warning This assumes this TrackState stores it's own calibrated + /// measurement. **If storage is shared with another TrackState, both + /// will be overwritten!**. Also assumes none of the calibrated + /// components is *invalid* (i.e. unset) for this TrackState. /// @note This does not set the reference surface. template > @@ -888,67 +810,306 @@ class TrackStateProxy { setProjector(meas.projector()); } + /// Allocate storage to be able to store a measurement of size @p measdim. + /// This must be called **before** setting the measurement content. void allocateCalibrated(std::size_t measdim) { m_traj->allocateCalibrated(m_istate, measdim); } - /// Getter/setter for chi2 value associated with the track state - /// This overload returns a mutable reference, which allows setting a new - /// value directly into the backing store. - /// @note this overload is only enabled in case the proxy is not read-only - /// @return Mutable reference to the chi2 value + /// @} + + /// @anchor track_state_share_copy + /// @name Sharing and copying + /// + /// Methods to share and copy track state components. Sharing means setting up + /// more than one track state to point to the same component. + /// + /// Shareable components are + /// - predicted parameters and covariance + /// - filtered parameters and covariance + /// - smoothed parameters and covariance + /// - jacobian + /// + /// See @ref TrackStatePropMask. + /// + /// @{ + + /// Share a shareable component **within** this track state + /// @param shareSource Which component to share from + /// @param shareTarget Which component to share as. This should be different from + /// as @p shareSource, e.g. predicted can be shared as filtered. template > - float& chi2() { - return component(); + void shareFrom(TrackStatePropMask shareSource, + TrackStatePropMask shareTarget) { + shareFrom(*this, shareSource, shareTarget); } - /// Getter for the chi2 value associated with the track state. - /// This overload returns a copy of the chi2 value, and thus does not allow - /// modification of the value in the backing storage. - /// @return the chi2 value of the track state - float chi2() const { return component(); } + /// Share a shareable component from another track state. + /// @param other Track state proxy to share component from + /// @param component Which component to share. + /// @note The track states both need to be stored in the + /// same @c MultiTrajectory instance + template > + void shareFrom(const TrackStateProxy& other, + TrackStatePropMask component) { + shareFrom(other, component, component); + } - /// Getter for the path length associated with the track state. - /// This overloaded is only enabled if not read-only, and returns a mutable - /// reference. - /// @return Mutable reference to the pathlength. - template > - double& pathLength() { - return component(); + /// Share a shareable component from another track state + /// @param other Track state proxy to share component(s) from + /// @param shareSource Which component to share from + /// @param shareTarget Which component to share as. This can be be different from + /// as @p shareSource, e.g. predicted can be shared as filtered. + /// @note Shareable components are predicted, filtered, smoothed, calibrated, jacobian, + /// or projector. See @c TrackStatePropMask. + template > + void shareFrom(const TrackStateProxy& other, + TrackStatePropMask shareSource, + TrackStatePropMask shareTarget) { + assert(m_traj == other.m_traj && + "Cannot share components across MultiTrajectories"); + + assert(ACTS_CHECK_BIT(other.getMask(), shareSource) && + "Source has incompatible allocation"); + + m_traj->self().shareFrom(m_istate, other.m_istate, shareSource, + shareTarget); } - /// Getter for the path length. Returns a copy of the path length value. - /// @return The path length of this track state - double pathLength() const { - return component(); + /// Copy the contents of another track state proxy into this one + /// @param other The other track state to copy from + /// @param mask An optional mask to determine what to copy from + /// @param onlyAllocated Whether to only copy allocated components + /// @note If the this track state proxy does not have compatible allocations + /// with the source track state proxy, and @p onlyAllocated is false, + /// an exception is thrown. + /// @note The mask parameter will not cause a copy of components that are + /// not allocated in the source track state proxy. + template > + void copyFrom(const track_state_proxy_t& other, + TrackStatePropMask mask = TrackStatePropMask::All, + bool onlyAllocated = true) { + using PM = TrackStatePropMask; + + if (onlyAllocated) { + auto dest = getMask(); + auto src = other.getMask() & + mask; // combine what we have with what we want to copy + + if (ACTS_CHECK_BIT(src, PM::Calibrated)) { + // on-demand allocate calibrated + dest |= PM::Calibrated; + } + + if ((static_cast>( + (src ^ dest) & src) != 0 || + dest == TrackStatePropMask::None || + src == TrackStatePropMask::None) && + mask != TrackStatePropMask::None) { + throw std::runtime_error( + "Attempt track state copy with incompatible allocations"); + } + + // we're sure now this has correct allocations, so just copy + if (ACTS_CHECK_BIT(src, PM::Predicted)) { + predicted() = other.predicted(); + predictedCovariance() = other.predictedCovariance(); + } + + if (ACTS_CHECK_BIT(src, PM::Filtered)) { + filtered() = other.filtered(); + filteredCovariance() = other.filteredCovariance(); + } + + if (ACTS_CHECK_BIT(src, PM::Smoothed)) { + smoothed() = other.smoothed(); + smoothedCovariance() = other.smoothedCovariance(); + } + + if (other.hasUncalibratedSourceLink()) { + setUncalibratedSourceLink(other.getUncalibratedSourceLink()); + } + + if (ACTS_CHECK_BIT(src, PM::Jacobian)) { + jacobian() = other.jacobian(); + } + + if (ACTS_CHECK_BIT(src, PM::Calibrated)) { + allocateCalibrated(other.calibratedSize()); + + // workaround for gcc8 bug: + // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86594 + auto* self = this; + visit_measurement(other.calibratedSize(), [&](auto N) { + constexpr int measdim = decltype(N)::value; + self->template calibrated() = + other.template calibrated(); + self->template calibratedCovariance() = + other.template calibratedCovariance(); + }); + + setProjectorBitset(other.projectorBitset()); + } + } else { + if (ACTS_CHECK_BIT(mask, PM::Predicted) && + has() && + other.template has()) { + predicted() = other.predicted(); + predictedCovariance() = other.predictedCovariance(); + } + + if (ACTS_CHECK_BIT(mask, PM::Filtered) && has() && + other.template has()) { + filtered() = other.filtered(); + filteredCovariance() = other.filteredCovariance(); + } + + if (ACTS_CHECK_BIT(mask, PM::Smoothed) && has() && + other.template has()) { + smoothed() = other.smoothed(); + smoothedCovariance() = other.smoothedCovariance(); + } + + if (other.hasUncalibratedSourceLink()) { + setUncalibratedSourceLink(other.getUncalibratedSourceLink()); + } + + if (ACTS_CHECK_BIT(mask, PM::Jacobian) && has() && + other.template has()) { + jacobian() = other.jacobian(); + } + + if (ACTS_CHECK_BIT(mask, PM::Calibrated) && + has() && + other.template has()) { + allocateCalibrated(other.calibratedSize()); + + // workaround for gcc8 bug: + // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=86594 + auto* self = this; + visit_measurement(other.calibratedSize(), [&](auto N) { + constexpr int measdim = decltype(N)::value; + self->template calibrated() = + other.template calibrated(); + self->template calibratedCovariance() = + other.template calibratedCovariance(); + }); + + setProjectorBitset(other.projectorBitset()); + } + } + + chi2() = other.chi2(); + pathLength() = other.pathLength(); + typeFlags() = other.typeFlags(); + + if (other.hasReferenceSurface()) { + setReferenceSurface(other.referenceSurface().getSharedPtr()); + } + + m_traj->copyDynamicFrom(m_istate, other.container(), other.index()); } - /// Getter for the type flags associated with the track state. - /// This overloaded is only enabled if not read-only, and returns a mutable - /// reference. - /// @return reference to the type flags. - template > - TrackStateType typeFlags() { - return TrackStateType{ - component()}; + /// @} + + /// @anchor track_state_proxy_generic_component + /// @name Track state proxy Generic component access + /// @{ + + /// Check if a component is set + /// @tparam key Hashed string key to check for + /// @return true if the component exists, false if not + template + constexpr bool has() const { + return m_traj->template has(m_istate); } - /// Getter for the type flags. Returns a copy of the type flags value. - /// @return The type flags of this track state - ConstTrackStateType typeFlags() const { - return ConstTrackStateType{ - component()}; + /// Check if a component is set + /// @param key Hashed string key to check for + /// @return true if the component exists, false if not + constexpr bool has(HashedString key) const { + return m_traj->has(key, m_istate); } - /// Get a mutable reference to the track state container backend - /// @return a mutable reference to the backend + /// Check if a component is set + /// @param key String key to check for + /// @note This might hash the @p key at runtime instead of compile-time + /// @return true if the component exists, false if not + constexpr bool has(std::string_view key) const { + return has(hashString(key)); + } + + /// Retrieve a mutable reference to a component + /// @tparam T The type of the component to access + /// @tparam key String key for the component to access + /// @return Mutable reference to the component given by @p key + template > + constexpr T& component() { + return m_traj->template component(m_istate); + } + + /// Retrieve a mutable reference to a component + /// @tparam T The type of the component to access + /// @param key String key for the component to access + /// @return Mutable reference to the component given by @p key + template > + constexpr T& component(HashedString key) { + return m_traj->template component(key, m_istate); + } + + /// Retrieve a mutable reference to a component + /// @tparam T The type of the component to access + /// @param key String key for the component to access + /// @note This might hash the @p key at runtime instead of compile-time + /// @return Mutable reference to the component given by @p key + template > + constexpr T& component(std::string_view key) { + return m_traj->template component(hashString(key), m_istate); + } + + /// Retrieve a const reference to a component + /// @tparam T The type of the component to access + /// @tparam key String key for the component to access + /// @return Const reference to the component given by @p key + template + constexpr const T& component() const { + return m_traj->template component(m_istate); + } + + /// Retrieve a const reference to a component + /// @tparam T The type of the component to access + /// @param key String key for the component to access + /// @return Const reference to the component given by @p key + template + constexpr const T& component(HashedString key) const { + return m_traj->template component(key, m_istate); + } + + /// Retrieve a const reference to a component + /// @tparam T The type of the component to access + /// @param key String key for the component to access + /// @note This might hash the @p key at runtime instead of compile-time + /// @return Const reference to the component given by @p key + template + constexpr const T& component(std::string_view key) const { + return m_traj->template component(hashString(key), m_istate); + } + + /// @} + + /// Return a mutable reference to the underlying backend container + /// @return A reference to the backend container template > MultiTrajectory& trajectory() { return *m_traj; } - /// Get a const reference to the track state container backend - /// @return a const reference to the backend + /// Return a const reference to the underlying backend container + /// @return A const reference to the backend container const MultiTrajectory& trajectory() const { return *m_traj; } /// Get a mutable reference to the track state container backend diff --git a/Core/include/Acts/EventData/VectorMultiTrajectory.hpp b/Core/include/Acts/EventData/VectorMultiTrajectory.hpp index d73e6d9d223..5a6aa87e62e 100644 --- a/Core/include/Acts/EventData/VectorMultiTrajectory.hpp +++ b/Core/include/Acts/EventData/VectorMultiTrajectory.hpp @@ -182,6 +182,7 @@ class VectorMultiTrajectoryBase { for (const auto& [key, value] : other.m_dynamic) { m_dynamic.insert({key, value->clone()}); } + m_dynamicKeys = other.m_dynamicKeys; }; VectorMultiTrajectoryBase(VectorMultiTrajectoryBase&& other) = default; @@ -333,6 +334,7 @@ class VectorMultiTrajectoryBase { // be handled in a smart way by moving but not sure. std::vector> m_referenceSurfaces; + std::vector m_dynamicKeys; std::unordered_map> m_dynamic; }; diff --git a/Core/include/Acts/EventData/VectorTrackContainer.hpp b/Core/include/Acts/EventData/VectorTrackContainer.hpp index 0183ce8d801..50bdbc2cd88 100644 --- a/Core/include/Acts/EventData/VectorTrackContainer.hpp +++ b/Core/include/Acts/EventData/VectorTrackContainer.hpp @@ -188,6 +188,7 @@ class VectorTrackContainerBase { std::unordered_map> m_dynamic; + std::vector m_dynamicKeys; }; } // namespace detail_vtc diff --git a/Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp b/Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp index 8b97da0f8b7..8693a396104 100644 --- a/Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp +++ b/Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp @@ -463,7 +463,7 @@ class MultiTrajectoryTestsCommon { } BOOST_CHECK(ts.hasProjector()); - ActsMatrix fullProj; + ActsMatrix fullProj; fullProj.setZero(); { Acts::GeometryContext gctx; diff --git a/Core/include/Acts/Geometry/ApproachDescriptor.hpp b/Core/include/Acts/Geometry/ApproachDescriptor.hpp index 662cbc0faac..5f05a941672 100644 --- a/Core/include/Acts/Geometry/ApproachDescriptor.hpp +++ b/Core/include/Acts/Geometry/ApproachDescriptor.hpp @@ -43,17 +43,16 @@ class ApproachDescriptor { /// @param position is the position from start of the search /// @param direction is the direction at the start of the search /// @param bcheck is the boundary check directive - /// @param pLimit The path limit - /// @param oLimit The overstep limit - /// @param tolerance The surface tolerance + /// @param nearLimit The minimum distance for an intersection to be considered + /// @param farLimit The maximum distance for an intersection to be considered /// /// @return is a surface intersection virtual SurfaceIntersection approachSurface(const GeometryContext& gctx, const Vector3& position, const Vector3& direction, const BoundaryCheck& bcheck, - double pLimit, double oLimit, - double tolerance) const = 0; + double nearLimit, + double farLimit) const = 0; /// Get all the contained surfaces /// @return all contained surfaces of this approach descriptor diff --git a/Core/include/Acts/Geometry/BoundarySurfaceT.hpp b/Core/include/Acts/Geometry/BoundarySurfaceT.hpp index 09578def199..9179aec59f7 100644 --- a/Core/include/Acts/Geometry/BoundarySurfaceT.hpp +++ b/Core/include/Acts/Geometry/BoundarySurfaceT.hpp @@ -103,13 +103,12 @@ class BoundarySurfaceT { /// /// @param gctx The current geometry context object, e.g. alignment /// @param pos The global position on surface - /// @param mom The direction on the surface - /// @param dir is an additional direction corrective + /// @param dir The direction on the surface /// /// @return The attached volume at that position virtual const volume_t* attachedVolume(const GeometryContext& gctx, - const Vector3& pos, const Vector3& mom, - Direction dir) const; + const Vector3& pos, + const Vector3& dir) const; /// templated onBoundary method /// @@ -181,11 +180,10 @@ void BoundarySurfaceT::attachVolumeArray( template const volume_t* BoundarySurfaceT::attachedVolume( - const GeometryContext& gctx, const Vector3& pos, const Vector3& mom, - Direction dir) const { + const GeometryContext& gctx, const Vector3& pos, const Vector3& dir) const { const volume_t* attVolume = nullptr; // dot product with normal vector to distinguish inside/outside - if ((surfaceRepresentation().normal(gctx, pos)).dot(dir * mom) > 0.) { + if ((surfaceRepresentation().normal(gctx, pos)).dot(dir) > 0.) { attVolume = m_alongVolumeArray ? m_alongVolumeArray->object(pos).get() : m_alongVolume; } else { diff --git a/Core/include/Acts/Geometry/GenericApproachDescriptor.hpp b/Core/include/Acts/Geometry/GenericApproachDescriptor.hpp index 2c69d92d068..53dffaa1eaa 100644 --- a/Core/include/Acts/Geometry/GenericApproachDescriptor.hpp +++ b/Core/include/Acts/Geometry/GenericApproachDescriptor.hpp @@ -57,17 +57,16 @@ class GenericApproachDescriptor : public ApproachDescriptor { /// @param position The global position to start the approach from /// @param direction The momentum vector /// @param bcheck The boundary check prescription - /// @param pLimit The path limit - /// @param oLimit The overstep limit - /// @param tolerance The surface tolerance + /// @param nearLimit The minimum distance for an intersection to be considered + /// @param farLimit The maximum distance for an intersection to be considered /// /// @return : a @c SurfaceIntersection SurfaceIntersection approachSurface(const GeometryContext& gctx, const Vector3& position, const Vector3& direction, const BoundaryCheck& bcheck, - double pLimit, double oLimit, - double tolerance) const override; + double nearLimit, + double farLimit) const override; /// return all contained surfaces of this approach descriptor const std::vector& containedSurfaces() const override; diff --git a/Core/include/Acts/Geometry/NavigationLayer.hpp b/Core/include/Acts/Geometry/NavigationLayer.hpp index 57d020a4665..24f54653006 100644 --- a/Core/include/Acts/Geometry/NavigationLayer.hpp +++ b/Core/include/Acts/Geometry/NavigationLayer.hpp @@ -26,7 +26,6 @@ namespace Acts { /// Class to be used for gaps in Volumes as a navigational link. /// Navigation Layers have a surface representation, but should usually never be /// propagated to. - class NavigationLayer : public Layer { public: /// Factory Constructor - the surface representation is given by pointer diff --git a/Core/include/Acts/Navigation/DetectorNavigator.hpp b/Core/include/Acts/Navigation/DetectorNavigator.hpp index 1301017f8bc..1b0ce0322f6 100644 --- a/Core/include/Acts/Navigation/DetectorNavigator.hpp +++ b/Core/include/Acts/Navigation/DetectorNavigator.hpp @@ -83,16 +83,6 @@ class DetectorNavigator { return result; } - void resetState(State& state, const GeometryContext& /*geoContext*/, - const Vector3& /*pos*/, const Vector3& /*dir*/, - const Surface* /*ssurface*/, - const Surface* /*tsurface*/) const { - // Reset everything first - state = State(); - - // TODO fill state - } - const Surface* currentSurface(const State& state) const { return state.currentSurface; } @@ -202,24 +192,24 @@ class DetectorNavigator { << posInfo(state, stepper) << "stepping through portal"); nState.surfaceCandidates.clear(); - nState.surfaceCandidate = nState.surfaceCandidates.cend(); + nState.surfaceCandidateIndex = 0; nState.currentPortal->updateDetectorVolume(state.geoContext, nState); initializeTarget(state, stepper); } - for (; nState.surfaceCandidate != nState.surfaceCandidates.cend(); - ++nState.surfaceCandidate) { + for (; nState.surfaceCandidateIndex != nState.surfaceCandidates.size(); + ++nState.surfaceCandidateIndex) { // Screen output how much is left to try ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper) - << std::distance(nState.surfaceCandidate, - nState.surfaceCandidates.cend()) + << (nState.surfaceCandidates.size() - + nState.surfaceCandidateIndex) << " out of " << nState.surfaceCandidates.size() << " surfaces remain to try."); // Take the surface - const auto& c = *(nState.surfaceCandidate); + const auto& c = nState.surfaceCandidate(); const auto& surface = (c.surface != nullptr) ? (*c.surface) : (c.portal->surface()); // Screen output which surface you are on @@ -271,7 +261,7 @@ class DetectorNavigator { return; } - if (nState.surfaceCandidate == nState.surfaceCandidates.cend()) { + if (nState.surfaceCandidateIndex == nState.surfaceCandidates.size()) { ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper) << "no surface candidates - waiting for target call"); @@ -281,12 +271,12 @@ class DetectorNavigator { const Portal* nextPortal = nullptr; const Surface* nextSurface = nullptr; bool isPortal = false; - bool boundaryCheck = nState.surfaceCandidate->boundaryCheck; + bool boundaryCheck = nState.surfaceCandidate().boundaryCheck; - if (nState.surfaceCandidate->surface != nullptr) { - nextSurface = nState.surfaceCandidate->surface; - } else if (nState.surfaceCandidate->portal != nullptr) { - nextPortal = nState.surfaceCandidate->portal; + if (nState.surfaceCandidate().surface != nullptr) { + nextSurface = nState.surfaceCandidate().surface; + } else if (nState.surfaceCandidate().portal != nullptr) { + nextPortal = nState.surfaceCandidate().portal; nextSurface = &nextPortal->surface(); isPortal = true; } else { @@ -299,7 +289,7 @@ class DetectorNavigator { // TODO not sure about the boundary check auto surfaceStatus = stepper.updateSurfaceStatus( state.stepping, *nextSurface, - nState.surfaceCandidate->objectIntersection.index(), + nState.surfaceCandidate().objectIntersection.index(), state.options.direction, BoundaryCheck(boundaryCheck), state.options.surfaceTolerance, logger()); @@ -327,7 +317,7 @@ class DetectorNavigator { ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper) << "current surface set to " << nState.currentSurface->geometryId()); - ++nState.surfaceCandidate; + ++nState.surfaceCandidateIndex; } } } @@ -426,7 +416,7 @@ class DetectorNavigator { return pathToA < pathToB; }); // Set the surface candidate - nState.surfaceCandidate = nCandidates.begin(); + nState.surfaceCandidateIndex = 0; } template diff --git a/Core/include/Acts/Navigation/NavigationState.hpp b/Core/include/Acts/Navigation/NavigationState.hpp index ae72ce67657..1cd752d940a 100644 --- a/Core/include/Acts/Navigation/NavigationState.hpp +++ b/Core/include/Acts/Navigation/NavigationState.hpp @@ -9,12 +9,14 @@ #pragma once #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Units.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Surfaces/BoundaryCheck.hpp" #include "Acts/Utilities/Delegate.hpp" #include "Acts/Utilities/Intersection.hpp" #include +#include #include namespace Acts { @@ -80,16 +82,20 @@ struct NavigationState { /// That are the candidate surfaces to process SurfaceCandidates surfaceCandidates = {}; - SurfaceCandidates::const_iterator surfaceCandidate = surfaceCandidates.cend(); + std::size_t surfaceCandidateIndex = 0; /// Boundary directives for surfaces BoundaryCheck surfaceBoundaryCheck = BoundaryCheck(true); /// An overstep tolerance - ActsScalar overstepTolerance = -0.1; + ActsScalar overstepTolerance = -100 * UnitConstants::um; /// Auxiliary attached information std::any auxiliary; + + const SurfaceCandidate& surfaceCandidate() const { + return surfaceCandidates.at(surfaceCandidateIndex); + } }; } // namespace Experimental diff --git a/Core/include/Acts/Navigation/NavigationStateUpdaters.hpp b/Core/include/Acts/Navigation/NavigationStateUpdaters.hpp index 0658ccfe8de..3edfefb8c0a 100644 --- a/Core/include/Acts/Navigation/NavigationStateUpdaters.hpp +++ b/Core/include/Acts/Navigation/NavigationStateUpdaters.hpp @@ -12,6 +12,7 @@ #include "Acts/Definitions/Common.hpp" #include "Acts/Detector/Portal.hpp" #include "Acts/Navigation/NavigationDelegates.hpp" +#include "Acts/Navigation/NavigationState.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/BinningType.hpp" #include "Acts/Utilities/Enumerate.hpp" @@ -38,9 +39,10 @@ inline void updateCandidates(const GeometryContext& gctx, NavigationState& nState) { const auto& position = nState.position; const auto& direction = nState.direction; - auto& nCandidates = nState.surfaceCandidates; - for (auto& c : nCandidates) { + NavigationState::SurfaceCandidates nextSurfaceCandidates; + + for (NavigationState::SurfaceCandidate c : nState.surfaceCandidates) { // Get the surface representation: either native surfcae of portal const Surface& sRep = c.surface != nullptr ? *c.surface : c.portal->surface(); @@ -50,7 +52,14 @@ inline void updateCandidates(const GeometryContext& gctx, auto sIntersection = sRep.intersect(gctx, position, direction, c.boundaryCheck, s_onSurfaceTolerance); c.objectIntersection = sIntersection[c.objectIntersection.index()]; + + if (c.objectIntersection && + c.objectIntersection.pathLength() > nState.overstepTolerance) { + nextSurfaceCandidates.emplace_back(std::move(c)); + } } + + nState.surfaceCandidates = std::move(nextSurfaceCandidates); } /// @brief This sets a single object, e.g. single surface or single volume diff --git a/Core/include/Acts/Propagator/DirectNavigator.hpp b/Core/include/Acts/Propagator/DirectNavigator.hpp index 9bd1e500120..d3b95b392cb 100644 --- a/Core/include/Acts/Propagator/DirectNavigator.hpp +++ b/Core/include/Acts/Propagator/DirectNavigator.hpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-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 @@ -24,13 +24,11 @@ namespace Acts { -/// DirectNavigator class +/// This is a fully guided navigator that progresses through a pre-given +/// sequence of surfaces. /// -/// This is a fully guided navigator that progresses through -/// a pre-given sequence of surfaces. -/// -/// This can either be used as a validation tool, for truth -/// tracking, or track refitting +/// This can either be used as a validation tool, for truth tracking, or track +/// refitting. class DirectNavigator { public: /// The sequentially crossed surfaces @@ -41,10 +39,9 @@ class DirectNavigator { getDefaultLogger("DirectNavigator", Logging::INFO)) : m_logger{std::move(_logger)} {} - /// Nested Actor struct, called Initializer + /// @brief Nested Actor struct, called Initializer /// - /// This is needed for the initialization of the - /// surface sequence + /// This is needed for the initialization of the surface sequence. struct Initializer { /// The Surface sequence SurfaceSequence navSurfaces = {}; @@ -94,14 +91,13 @@ class DirectNavigator { } }; - /// Nested State struct + /// @brief Nested State struct /// - /// It acts as an internal state which is - /// created for every propagation/extrapolation step - /// and keep thread-local navigation information + /// It acts as an internal state which is created for every + /// propagation/extrapolation step and keep thread-local navigation + /// information struct State { - /// Externally provided surfaces - expected to be ordered - /// along the path + /// Externally provided surfaces - expected to be ordered along the path SurfaceSequence navSurfaces = {}; /// Iterator the next surface @@ -194,8 +190,7 @@ class DirectNavigator { void initialize(propagator_state_t& state, const stepper_t& stepper) const { (void)stepper; - // Call the navigation helper prior to actual navigation - ACTS_VERBOSE(volInfo(state) << "Initialization."); + ACTS_VERBOSE(volInfo(state) << "initialize"); // We set the current surface to the start surface state.navigation.currentSurface = state.navigation.startSurface; @@ -215,8 +210,7 @@ class DirectNavigator { /// @param [in] stepper Stepper in use template void preStep(propagator_state_t& state, const stepper_t& stepper) const { - // Screen output - ACTS_VERBOSE("Entering navigator::target."); + ACTS_VERBOSE(volInfo(state) << "pre step"); // Navigator target always resets the current surface state.navigation.currentSurface = nullptr; @@ -230,12 +224,12 @@ class DirectNavigator { // Establish & update the surface status // TODO we do not know the intersection index - passing the closer one const auto& surface = **state.navigation.navSurfaceIter; + const double farLimit = std::numeric_limits::max(); const auto index = chooseIntersection( state.geoContext, surface, stepper.position(state.stepping), state.options.direction * stepper.direction(state.stepping), - BoundaryCheck(false), std::numeric_limits::max(), - stepper.overstepLimit(state.stepping), + BoundaryCheck(false), m_nearLimit, farLimit, state.options.surfaceTolerance) .index(); auto surfaceStatus = stepper.updateSurfaceStatus( @@ -272,8 +266,7 @@ class DirectNavigator { /// @param [in] stepper Stepper in use template void postStep(propagator_state_t& state, const stepper_t& stepper) const { - // Screen output - ACTS_VERBOSE("Entering navigator::postStep."); + ACTS_VERBOSE(volInfo(state) << "post step"); // Navigator post step always resets the current surface state.navigation.currentSurface = nullptr; @@ -288,12 +281,12 @@ class DirectNavigator { // Establish the surface status // TODO we do not know the intersection index - passing the closer one const auto& surface = **state.navigation.navSurfaceIter; + const double farLimit = std::numeric_limits::max(); const auto index = chooseIntersection( state.geoContext, surface, stepper.position(state.stepping), state.options.direction * stepper.direction(state.stepping), - BoundaryCheck(false), std::numeric_limits::max(), - stepper.overstepLimit(state.stepping), + BoundaryCheck(false), m_nearLimit, farLimit, state.options.surfaceTolerance) .index(); auto surfaceStatus = stepper.updateSurfaceStatus( @@ -321,24 +314,22 @@ class DirectNavigator { private: template std::string volInfo(const propagator_state_t& state) const { - return (state.navigation.currentVolume + return (state.navigation.currentVolume != nullptr ? state.navigation.currentVolume->volumeName() : "No Volume") + " | "; } - ObjectIntersection chooseIntersection(const GeometryContext& gctx, - const Surface& surface, - const Vector3& position, - const Vector3& direction, - const BoundaryCheck& bcheck, - double pLimit, double oLimit, - double tolerance) const { + ObjectIntersection chooseIntersection( + const GeometryContext& gctx, const Surface& surface, + const Vector3& position, const Vector3& direction, + const BoundaryCheck& bcheck, double nearLimit, double farLimit, + double tolerance) const { auto intersections = surface.intersect(gctx, position, direction, bcheck, tolerance); for (auto& intersection : intersections.split()) { - if (detail::checkIntersection(intersection, pLimit, oLimit, tolerance, + if (detail::checkIntersection(intersection, nearLimit, farLimit, logger())) { return intersection; } @@ -350,6 +341,12 @@ class DirectNavigator { const Logger& logger() const { return *m_logger; } std::unique_ptr m_logger; + + // TODO https://github.com/acts-project/acts/issues/2738 + /// Distance limit to discard intersections "behind us" + /// @note this is only necessary because some surfaces have more than one + /// intersection + double m_nearLimit = -100 * UnitConstants::um; }; } // namespace Acts diff --git a/Core/include/Acts/Propagator/Navigator.hpp b/Core/include/Acts/Propagator/Navigator.hpp index dc9a61a436a..2babb8d5e3c 100644 --- a/Core/include/Acts/Propagator/Navigator.hpp +++ b/Core/include/Acts/Propagator/Navigator.hpp @@ -51,42 +51,34 @@ struct NavigationOptions { /// object to check against: at end const object_t* endObject = nullptr; - /// Target surface to exclude - const Surface* targetSurface = nullptr; /// External surface identifier for which the boundary check is ignored std::vector externalSurfaces = {}; - /// The maximum path limit for this navigation step - double pathLimit = std::numeric_limits::max(); - - /// The overstep tolerance for this navigation step - double overstepLimit = 0; + /// The minimum distance for a surface to be considered + double nearLimit = 0; + /// The maximum distance for a surface to be considered + double farLimit = std::numeric_limits::max(); /// Force intersection with boundaries bool forceIntersectBoundaries = false; }; -/// Navigator class +/// @brief Steers the propagation through the geometry by adjusting the step +/// size and providing the next surface to be targeted. /// -/// This is an Actor to be added to the ActorList in order to navigate -/// through the static tracking geometry setup. +/// The Navigator is part of the propagation and responsible for steering +/// the step size in order to encounter all the relevant surfaces which are +/// intersected by the trajectory. /// /// The current navigation stage is cached in the state struct and updated -/// when necessary. If any surface in the extrapolation flow is hit, it is -/// set to the propagation satate, such that other actors can deal with it. -/// This navigation actor thus always needs to run first! -/// It does two things: it figures out the order of volumes, layers and -/// surfaces. +/// when necessary. If any surface in the extrapolation flow is hit, it is +/// set to the navigation state, such that other actors can deal with it. /// -/// The current target surface is the surface pointed to by of the index -/// for the surfaces, layers or volume boundaries. -/// If a surface is found, the state.navigation.currentSurface -/// pointer is set. This enables subsequent actors to react. Secondly, this -/// actor uses the ordered indices to figure out which surface, layer or -/// volume boundary is _supposed_ to be hit next. It then sets the maximum -/// step size to the path length found out by straight line intersection. If -/// the state is not on surface, it also re-computes the step size, to make -/// sure we end up at the desired surface. +/// The current target surface is referenced by an index which points into +/// the navigation candidates. The navigation candidates are ordered by the +/// path length to the surface. If a surface is hit, the +/// `state.navigation.currentSurface` pointer is set. This actors to observe +/// that we are on a surface. /// class Navigator { public: @@ -115,7 +107,6 @@ class Navigator { /// Tracking Geometry for this Navigator std::shared_ptr trackingGeometry{nullptr}; - /// Configuration for this Navigator /// stop at every sensitive surface (whether it has material or not) bool resolveSensitive = true; /// stop at every material surface (whether it is passive or not) @@ -128,11 +119,10 @@ class Navigator { BoundaryCheck boundaryCheckLayerResolving = BoundaryCheck(true); }; - /// Nested State struct + /// @brief Nested State struct /// - /// It acts as an internal state which is - /// created for every propagation/extrapolation step - /// and keep thread-local navigation information + /// It acts as an internal state which is created for every propagation and + /// meant to keep thread-local navigation information. struct State { // Navigation on surface level /// the vector of navigation surfaces to work through @@ -261,15 +251,13 @@ class Navigator { std::pair(geoid.layer(), geoid)); } - /// @brief Initialize call - start of propagation + /// @brief Initialize call - start of navigation /// /// @tparam propagator_state_t The state type of the propagator /// @tparam stepper_t The type of stepper used for the propagation /// /// @param [in,out] state is the propagation state object /// @param [in] stepper Stepper in use - /// - /// @return boolean return triggers exit to stepper template void initialize(propagator_state_t& state, const stepper_t& stepper) const { // Call the navigation helper prior to actual navigation @@ -389,7 +377,7 @@ class Navigator { state.navigation.currentSurface = nullptr; } - /// @brief Navigator post step call, will be called in two modes + /// @brief Navigator post step call /// /// (a) It initializes the Navigation stream if start volume is /// not yet defined: @@ -490,7 +478,7 @@ class Navigator { auto boundary = state.navigation.navBoundary().second; state.navigation.currentVolume = boundary->attachedVolume( state.geoContext, stepper.position(state.stepping), - stepper.direction(state.stepping), state.options.direction); + state.options.direction * stepper.direction(state.stepping)); // No volume anymore : end of known world if (!state.navigation.currentVolume) { ACTS_VERBOSE( @@ -743,7 +731,7 @@ class Navigator { navOpts.resolveMaterial = m_cfg.resolveMaterial; navOpts.resolvePassive = m_cfg.resolvePassive; navOpts.endObject = state.navigation.targetSurface; - navOpts.overstepLimit = stepper.overstepLimit(state.stepping); + navOpts.nearLimit = stepper.overstepLimit(state.stepping); double opening_angle = 0; // Preliminary version of the frustum opening angle estimation. @@ -928,9 +916,9 @@ class Navigator { NavigationOptions navOpts; // Exclude the current surface in case it's a boundary navOpts.startObject = state.navigation.currentSurface; - navOpts.pathLimit = + navOpts.nearLimit = stepper.overstepLimit(state.stepping); + navOpts.farLimit = stepper.getStepSize(state.stepping, ConstrainedStep::aborter); - navOpts.overstepLimit = stepper.overstepLimit(state.stepping); navOpts.forceIntersectBoundaries = state.navigation.forceIntersectBoundaries; @@ -1134,13 +1122,13 @@ class Navigator { navOpts.externalSurfaces.push_back(itSurface->second); } } + // No overstepping on start layer, otherwise ask the stepper + navOpts.nearLimit = (cLayer != nullptr) + ? state.options.surfaceTolerance + : stepper.overstepLimit(state.stepping); // Check the limit - navOpts.pathLimit = + navOpts.farLimit = stepper.getStepSize(state.stepping, ConstrainedStep::aborter); - // No overstepping on start layer, otherwise ask the stepper - navOpts.overstepLimit = (cLayer != nullptr) - ? state.options.surfaceTolerance - : stepper.overstepLimit(state.stepping); // get the surfaces state.navigation.navSurfaces = navLayer->compatibleSurfaces( @@ -1204,11 +1192,9 @@ class Navigator { navOpts.resolveMaterial = m_cfg.resolveMaterial; navOpts.resolvePassive = m_cfg.resolvePassive; navOpts.startObject = startLayer; - // Set also the target surface - navOpts.targetSurface = state.navigation.targetSurface; - navOpts.pathLimit = + navOpts.nearLimit = stepper.overstepLimit(state.stepping); + navOpts.farLimit = stepper.getStepSize(state.stepping, ConstrainedStep::aborter); - navOpts.overstepLimit = stepper.overstepLimit(state.stepping); // Request the compatible layers state.navigation.navLayers = state.navigation.currentVolume->compatibleLayers( diff --git a/Core/include/Acts/Propagator/Propagator.hpp b/Core/include/Acts/Propagator/Propagator.hpp index a3c914a8cc0..8eee3fcf2db 100644 --- a/Core/include/Acts/Propagator/Propagator.hpp +++ b/Core/include/Acts/Propagator/Propagator.hpp @@ -23,8 +23,8 @@ #include "Acts/Propagator/ActionList.hpp" #include "Acts/Propagator/StandardAborters.hpp" #include "Acts/Propagator/StepperConcept.hpp" +#include "Acts/Propagator/VoidNavigator.hpp" #include "Acts/Propagator/detail/ParameterTraits.hpp" -#include "Acts/Propagator/detail/VoidPropagatorComponents.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/Result.hpp" @@ -203,7 +203,7 @@ struct PropagatorOptions : public PropagatorPlainOptions { /// - a type mapping for: (initial track parameter type and destination /// surface type) -> type of internal state object /// -template +template class Propagator final { /// Re-define bound track parameters dependent on the stepper using StepperBoundTrackParameters = diff --git a/Core/include/Acts/Propagator/RiddersPropagator.hpp b/Core/include/Acts/Propagator/RiddersPropagator.hpp index 7934bc266c5..01c050ddfec 100644 --- a/Core/include/Acts/Propagator/RiddersPropagator.hpp +++ b/Core/include/Acts/Propagator/RiddersPropagator.hpp @@ -85,8 +85,8 @@ class RiddersPropagator { /// /// @param [in] propagator Underlying propagator that will be used /// @param [in] config Config for the Ridders propagation - RiddersPropagator(propagator_t& propagator, Config config = {}) - : m_propagator(propagator), m_config(std::move(config)) {} + RiddersPropagator(propagator_t propagator, Config config = {}) + : m_propagator(std::move(propagator)), m_config(std::move(config)) {} /// @brief Constructor building a propagator /// @@ -96,10 +96,10 @@ class RiddersPropagator { /// @param [in] stepper Stepper that will be used /// @param [in] navigator Navigator that will be used /// @param [in] config Config for the Ridders propagation - template + template RiddersPropagator(stepper_t stepper, navigator_t navigator = navigator_t(), Config config = {}) - : m_propagator(Propagator(stepper, navigator)), + : m_propagator(std::move(stepper), std::move(navigator)), m_config(std::move(config)) {} /// @brief Propagation method targeting curvilinear parameters diff --git a/Core/include/Acts/Propagator/StandardAborters.hpp b/Core/include/Acts/Propagator/StandardAborters.hpp index 7c93b79198e..c1bafd7603d 100644 --- a/Core/include/Acts/Propagator/StandardAborters.hpp +++ b/Core/include/Acts/Propagator/StandardAborters.hpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2017-2018 CERN for the benefit of the Acts project +// Copyright (C) 2017-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 @@ -67,15 +67,20 @@ struct PathLimitReached { } }; -/// This is the condition that the Surface has been reached -/// it then triggers an propagation abort of the propagation +/// This is the condition that the Surface has been reached it then triggers a +/// propagation abort struct SurfaceReached { const Surface* surface = nullptr; BoundaryCheck boundaryCheck = BoundaryCheck(true); - std::optional overstepLimit; + + // TODO https://github.com/acts-project/acts/issues/2738 + /// Distance limit to discard intersections "behind us" + /// @note this is only necessary because some surfaces have more than one + /// intersection + std::optional overrideNearLimit; SurfaceReached() = default; - SurfaceReached(double oLimit) : overstepLimit(oLimit) {} + SurfaceReached(double nearLimit) : overrideNearLimit(nearLimit) {} /// boolean operator for abort condition without using the result /// @@ -101,10 +106,14 @@ struct SurfaceReached { return true; } - const double pLimit = + // not blindly using the stepper overstep limit here because it does not + // always work for perigee surfaces. + // note: the near limit is necessary for surfaces with more than one + // intersections in order to discard the ones which are behind us. + const double nearLimit = + overrideNearLimit.value_or(stepper.overstepLimit(state.stepping)); + const double farLimit = state.stepping.stepSize.value(ConstrainedStep::aborter); - const double oLimit = - overstepLimit.value_or(stepper.overstepLimit(state.stepping)); const double tolerance = state.options.surfaceTolerance; const auto sIntersection = surface->intersect( @@ -113,33 +122,40 @@ struct SurfaceReached { boundaryCheck, tolerance); const auto closest = sIntersection.closest(); + bool reached = false; + if (closest.status() == Intersection3D::Status::onSurface) { const double distance = closest.pathLength(); ACTS_VERBOSE( "SurfaceReached aborter | " "Target surface reached at distance (tolerance) " << distance << " (" << tolerance << ")"); - return true; + reached = true; } + bool intersectionFound = false; + for (const auto& intersection : sIntersection.split()) { if (intersection && - detail::checkIntersection(intersection.intersection(), pLimit, oLimit, - tolerance, logger)) { + detail::checkIntersection(intersection.intersection(), nearLimit, + farLimit, logger)) { stepper.updateStepSize(state.stepping, intersection.pathLength(), ConstrainedStep::aborter, false); ACTS_VERBOSE( "SurfaceReached aborter | " "Target stepSize (surface) updated to " << stepper.outputStepSize(state.stepping)); - return false; + intersectionFound = true; + break; } } - ACTS_VERBOSE( - "SurfaceReached aborter | " - "Target intersection not found. Maybe next time?"); - return false; + if (!intersectionFound) { + ACTS_VERBOSE( + "SurfaceReached aborter | " + "Target intersection not found. Maybe next time?"); + } + return reached; } }; diff --git a/Core/include/Acts/Propagator/SurfaceCollector.hpp b/Core/include/Acts/Propagator/SurfaceCollector.hpp index 130fad7ebef..4f996019233 100644 --- a/Core/include/Acts/Propagator/SurfaceCollector.hpp +++ b/Core/include/Acts/Propagator/SurfaceCollector.hpp @@ -8,6 +8,7 @@ #pragma once +#include "Acts/Propagator/Propagator.hpp" #include "Acts/Surfaces/Surface.hpp" #include @@ -94,6 +95,10 @@ struct SurfaceCollector { void operator()(propagator_state_t& state, const stepper_t& stepper, const navigator_t& navigator, result_type& result, const Logger& logger) const { + if (state.stage == PropagatorStage::postPropagation) { + return; + } + auto currentSurface = navigator.currentSurface(state.navigation); // The current surface has been assigned by the navigator diff --git a/Core/include/Acts/Propagator/detail/VoidPropagatorComponents.hpp b/Core/include/Acts/Propagator/VoidNavigator.hpp similarity index 96% rename from Core/include/Acts/Propagator/detail/VoidPropagatorComponents.hpp rename to Core/include/Acts/Propagator/VoidNavigator.hpp index f8f44f275a6..fd8d7fc50e8 100644 --- a/Core/include/Acts/Propagator/detail/VoidPropagatorComponents.hpp +++ b/Core/include/Acts/Propagator/VoidNavigator.hpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2018 CERN for the benefit of the Acts project +// Copyright (C) 2018-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 @@ -9,7 +9,8 @@ #pragma once namespace Acts { -namespace detail { + +class Surface; /// @brief The void navigator struct as a default navigator /// @@ -107,5 +108,4 @@ struct VoidNavigator { const stepper_t& /*stepper*/) const {} }; -} // namespace detail } // namespace Acts diff --git a/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp b/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp index e5b61043844..12a6db38f7c 100644 --- a/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp +++ b/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp @@ -138,5 +138,18 @@ void transportCovarianceToCurvilinear(BoundSquareMatrix& boundCovariance, BoundToFreeMatrix& boundToFreeJacobian, const Vector3& direction); +/// Convert bound track parameters to another bound surface. +/// @pre The @p targetSurface must intersect with the surface attached to +/// @p boundParameters, and the parameters must be on-surface on the +/// target surface. +/// @param gctx The geometry context. +/// @param boundParameters The bound track parameters to convert. +/// @param targetSurface The target surface. +/// @param bField The magnetic field at the target surface. +/// @return The converted bound track parameters. +Result boundToBoundConversion( + const GeometryContext& gctx, const BoundTrackParameters& boundParameters, + const Surface& targetSurface, const Vector3& bField = Vector3::Zero()); + } // namespace detail } // namespace Acts diff --git a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp index 359233ad033..d2f3a7ef2f1 100644 --- a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp +++ b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp @@ -36,8 +36,8 @@ Acts::Intersection3D::Status updateSingleSurfaceStatus( const Surface& surface, std::uint8_t index, Direction navDir, const BoundaryCheck& bcheck, ActsScalar surfaceTolerance, const Logger& logger) { - ACTS_VERBOSE( - "Update single surface status for surface: " << surface.geometryId()); + ACTS_VERBOSE("Update single surface status for surface: " + << surface.geometryId() << " index " << (int)index); auto sIntersection = surface.intersect( state.geoContext, stepper.position(state), @@ -52,12 +52,11 @@ Acts::Intersection3D::Status updateSingleSurfaceStatus( } // Path and overstep limit checking - const double pLimit = state.stepSize.value(ConstrainedStep::aborter); - const double oLimit = stepper.overstepLimit(state); + const double nearLimit = stepper.overstepLimit(state); + const double farLimit = state.stepSize.value(ConstrainedStep::aborter); - if (sIntersection && - detail::checkIntersection(sIntersection.intersection(), pLimit, oLimit, - surfaceTolerance, logger)) { + if (sIntersection && detail::checkIntersection(sIntersection.intersection(), + nearLimit, farLimit, logger)) { ACTS_VERBOSE("Surface is reachable"); stepper.updateStepSize(state, sIntersection.pathLength(), ConstrainedStep::actor); diff --git a/Core/include/Acts/Seeding/SpacePointGrid.ipp b/Core/include/Acts/Seeding/SpacePointGrid.ipp index d1e4ccf86f1..5577df32acf 100644 --- a/Core/include/Acts/Seeding/SpacePointGrid.ipp +++ b/Core/include/Acts/Seeding/SpacePointGrid.ipp @@ -115,24 +115,24 @@ Acts::SpacePointGridCreator::createGrid( // seeds // FIXME: zBinSize must include scattering float zBinSize = config.cotThetaMax * config.deltaRMax; - int zBins = - std::max(1, (int)std::floor((config.zMax - config.zMin) / zBinSize)); + float zBins = + std::max(1.f, std::floor((config.zMax - config.zMin) / zBinSize)); - for (int bin = 0; bin <= zBins; bin++) { + for (int bin = 0; bin <= static_cast(zBins); bin++) { AxisScalar edge = - config.zMin + bin * ((config.zMax - config.zMin) / (float)zBins); + config.zMin + bin * ((config.zMax - config.zMin) / zBins); zValues.push_back(edge); } } else { // Use the zBinEdges defined in the config - for (auto& bin : config.zBinEdges) { + for (float bin : config.zBinEdges) { zValues.push_back(bin); } } detail::Axis - zAxis(zValues); + zAxis(std::move(zValues)); return std::make_unique>( - std::make_tuple(phiAxis, zAxis)); + std::make_tuple(std::move(phiAxis), std::move(zAxis))); } diff --git a/Core/include/Acts/Utilities/Intersection.hpp b/Core/include/Acts/Utilities/Intersection.hpp index a197ddf206c..7d276b3375d 100644 --- a/Core/include/Acts/Utilities/Intersection.hpp +++ b/Core/include/Acts/Utilities/Intersection.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Tolerance.hpp" #include "Acts/Utilities/Logger.hpp" #include @@ -230,26 +231,25 @@ namespace detail { /// path-limit and overstep-limit /// /// @tparam intersection_t Type of the intersection object -/// @tparam logger_t The logger type, which defaults to std::false_type to -/// prevent the generation of logging code /// /// @param intersection The intersection to check -/// @param pLimit The path-limit -/// @param oLimit The overstep-limit -/// @param tolerance The tolerance that is applied to the path-limit criterion +/// @param nearLimit The minimum distance for an intersection to be considered +/// @param farLimit The maximum distance for an intersection to be considered /// @param logger A optionally supplied logger which prints out a lot of infos -/// at VERBOSE level -template -bool checkIntersection(const intersection_t& intersection, double pLimit, - double oLimit, double tolerance, +/// at VERBOSE level +template +bool checkIntersection(const intersection_t& intersection, double nearLimit, + double farLimit, const Logger& logger = getDummyLogger()) { - const double cLimit = intersection.pathLength(); + const double distance = intersection.pathLength(); + // TODO why? + const double tolerance = s_onSurfaceTolerance; - ACTS_VERBOSE(" -> pLimit, oLimit, cLimit: " << pLimit << ", " << oLimit - << ", " << cLimit); + ACTS_VERBOSE(" -> near limit, far limit, distance: " + << nearLimit << ", " << farLimit << ", " << distance); - const bool coCriterion = cLimit > oLimit; - const bool cpCriterion = std::abs(cLimit) < std::abs(pLimit) + tolerance; + const bool coCriterion = distance > nearLimit; + const bool cpCriterion = std::abs(distance) < std::abs(farLimit) + tolerance; const bool accept = coCriterion && cpCriterion; @@ -259,12 +259,12 @@ bool checkIntersection(const intersection_t& intersection, double pLimit, ACTS_VERBOSE("Intersection is OUTSIDE limit because: "); if (!coCriterion) { ACTS_VERBOSE("- intersection path length " - << cLimit << " <= overstep limit " << oLimit); + << distance << " <= near limit " << nearLimit); } if (!cpCriterion) { ACTS_VERBOSE("- intersection path length " - << std::abs(cLimit) << " is over the path limit " - << (std::abs(pLimit) + tolerance) + << std::abs(distance) << " is over the far limit " + << (std::abs(farLimit) + tolerance) << " (including tolerance of " << tolerance << ")"); } } diff --git a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp index fd47b7ac85c..bf7d52d41d6 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp @@ -8,7 +8,6 @@ #pragma once -#include "Acts/Definitions/Algebra.hpp" #include "Acts/EventData/TrackParameters.hpp" #include "Acts/Utilities/Result.hpp" #include "Acts/Vertexing/AdaptiveGridTrackDensity.hpp" @@ -75,9 +74,6 @@ class AdaptiveGridDensityVertexFinder { // Map from input track to corresponding track density map std::unordered_map trackDensities; - // Map to store bool if track has passed track selection or not - std::unordered_map trackSelectionMap; - // Store tracks that have been removed from track collection. These // tracks will be removed from the main grid std::vector tracksToRemove; diff --git a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp index bfa6c0ae6eb..a6a3ba15240 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveGridDensityVertexFinder.ipp @@ -14,72 +14,65 @@ auto Acts::AdaptiveGridDensityVertexFinder::find( // Remove density contributions from tracks removed from track collection if (m_cfg.cacheGridStateForTrackRemoval && state.isInitialized && !state.tracksToRemove.empty()) { - // Bool to check if removable tracks, that pass selection, still exist - bool couldRemoveTracks = false; for (auto trk : state.tracksToRemove) { - if (!state.trackSelectionMap.at(trk)) { + auto it = state.trackDensities.find(trk); + if (it == state.trackDensities.end()) { // Track was never added to grid, so cannot remove it continue; } - couldRemoveTracks = true; - auto trackDensityMap = state.trackDensities.at(trk); - m_cfg.gridDensity.subtractTrack(trackDensityMap, state.mainDensityMap); - } - if (!couldRemoveTracks) { - // No tracks were removed anymore - // Return empty seed, i.e. vertex at constraint position - // (Note: Upstream finder should check for this break condition) - std::vector> seedVec{vertexingOptions.constraint}; - return seedVec; + m_cfg.gridDensity.subtractTrack(it->second, state.mainDensityMap); } } else { - state.mainDensityMap.clear(); + state.mainDensityMap = DensityMap(); // Fill with track densities for (auto trk : trackVector) { const BoundTrackParameters& trkParams = m_extractParameters(*trk); // Take only tracks that fulfill selection criteria if (!doesPassTrackSelection(trkParams)) { - if (m_cfg.cacheGridStateForTrackRemoval) { - state.trackSelectionMap[trk] = false; - } continue; } auto trackDensityMap = m_cfg.gridDensity.addTrack(trkParams, state.mainDensityMap); // Cache track density contribution to main grid if enabled if (m_cfg.cacheGridStateForTrackRemoval) { - state.trackDensities[trk] = trackDensityMap; - state.trackSelectionMap[trk] = true; + state.trackDensities[trk] = std::move(trackDensityMap); } } state.isInitialized = true; } + if (state.mainDensityMap.empty()) { + // No tracks passed selection + // Return empty seed, i.e. vertex at constraint position + // (Note: Upstream finder should check for this break condition) + std::vector> seedVec{vertexingOptions.constraint}; + return seedVec; + } + double z = 0; double t = 0; double zWidth = 0; - if (!state.mainDensityMap.empty()) { - if (!m_cfg.estimateSeedWidth) { - // Get z value of highest density bin - auto maxZTRes = m_cfg.gridDensity.getMaxZTPosition(state.mainDensityMap); - if (!maxZTRes.ok()) { - return maxZTRes.error(); - } - z = (*maxZTRes).first; - t = (*maxZTRes).second; - } else { - // Get z value of highest density bin and width - auto maxZTResAndWidth = - m_cfg.gridDensity.getMaxZTPositionAndWidth(state.mainDensityMap); - - if (!maxZTResAndWidth.ok()) { - return maxZTResAndWidth.error(); - } - z = (*maxZTResAndWidth).first.first; - t = (*maxZTResAndWidth).first.second; - zWidth = (*maxZTResAndWidth).second; + if (!m_cfg.estimateSeedWidth) { + // Get z value of highest density bin + auto maxZTRes = m_cfg.gridDensity.getMaxZTPosition(state.mainDensityMap); + + if (!maxZTRes.ok()) { + return maxZTRes.error(); + } + z = (*maxZTRes).first; + t = (*maxZTRes).second; + } else { + // Get z value of highest density bin and width + auto maxZTResAndWidth = + m_cfg.gridDensity.getMaxZTPositionAndWidth(state.mainDensityMap); + + if (!maxZTResAndWidth.ok()) { + return maxZTResAndWidth.error(); } + z = (*maxZTResAndWidth).first.first; + t = (*maxZTResAndWidth).first.second; + zWidth = (*maxZTResAndWidth).second; } // Construct output vertex, t will be 0 if temporalTrkGridSize == 1 diff --git a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp index 85366761bbe..bfd563479f1 100644 --- a/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp @@ -11,8 +11,7 @@ #include "Acts/EventData/TrackParameters.hpp" #include "Acts/Utilities/Result.hpp" -#include - +#include // TODO use flat unordered map #include namespace Acts { @@ -29,81 +28,65 @@ namespace Acts { /// Unlike in the GaussianGridTrackDensity, the overall density map /// grows adaptively when tracks densities are added to the grid. class AdaptiveGridTrackDensity { - // Assert odd spatial and temporal track grid size - public: - // The first (second) integer indicates the bin's z (t) position - using Bin = std::pair; - // Mapping between bins and track densities - using DensityMap = std::unordered_map>; - // Coordinates in the z-t plane; the t value will be set to 0 if time - // vertex seeding is disabled - using ZTPosition = std::pair; - // z-t position of a maximum and its width - using ZTPositionAndWidth = std::pair; + /// The first (second) integer indicates the bin's z (t) position + using Bin = std::pair; + /// Mapping between bins and track densities + using DensityMap = boost::container::flat_map; + /// Coordinates in the z-t plane; the t value will be set to 0 if time + /// vertex seeding is disabled + using ZTPosition = std::pair; + /// z-t position of a maximum and its width + using ZTPositionAndWidth = std::pair; + /// Optional grid size range + using GridSizeRange = + std::pair, std::optional>; /// The configuration struct struct Config { - /// Default constructor - Config() = default; - - // In total, a track is represented by a grid of size - // spatialTrkGridSize * temporalTrkGridSize - // Number of bins per track in z direction - unsigned int spatialTrkGridSize = 15; - - // Spatial extent of a bin in d0 and z0 direction, should always be set to a - // positive value - float spatialBinExtent = 0.; // mm - - // Number of bins per track in t direction - unsigned int temporalTrkGridSize = 1; - - // Temporal extent of a bin, should be set to 0 if time vertex seeding is - // disabled (i.e., if temporalTrkGridSize = 1) - float temporalBinExtent = 0.; // mm - - // Do NOT use just the z-bin with the highest - // track density, but instead check (up to) - // first three density maxima (only those that have - // a maximum relative deviation of 'relativeDensityDev' - // from the main maximum) and take the z-bin of the - // maximum with the highest surrounding density sum + /// Spatial extent of a bin in d0 and z0 direction, should always be set to + /// a positive value + double spatialBinExtent = 15 * UnitConstants::um; + + /// Number of standard deviations that the grid covers in z direction + double nSpatialTrkSigmas = 3.0; + + /// Temporal extent of a bin, not used if useTime == true + double temporalBinExtent = 19 * UnitConstants::mm; + + /// Number of standard deviations that the grid covers in t direction, not + /// used if useTime == true + double nTemporalTrkSigmas = 3.0; + + /// Spatial window for filling the density map + std::pair spatialWindow = {-250 * UnitConstants::mm, + 250 * UnitConstants::mm}; + /// Temporal window for filling the density map + std::pair temporalWindow = {-10 * UnitConstants::ns, + 10 * UnitConstants::ns}; + + /// Optional minimal and maximal number of bins in z direction + GridSizeRange spatialTrkGridSizeRange = {std::nullopt, std::nullopt}; + /// Optional minimal and maximal number of bins in time direction + GridSizeRange temporalTrkGridSizeRange = {std::nullopt, std::nullopt}; + + bool useTime = false; + + /// Do NOT use just the z-bin with the highest + /// track density, but instead check (up to) + /// first three density maxima (only those that have + /// a maximum relative deviation of 'relativeDensityDev' + /// from the main maximum) and take the z-bin of the + /// maximum with the highest surrounding density sum bool useHighestSumZPosition = false; - // The maximum relative density deviation from the main - // maximum to consider the second and third maximum for - // the highest-sum approach from above - float maxRelativeDensityDev = 0.01; + /// The maximum relative density deviation from the main + /// maximum to consider the second and third maximum for + /// the highest-sum approach from above + double maxRelativeDensityDev = 0.01; }; - AdaptiveGridTrackDensity(const Config& cfg) : m_cfg(cfg) { - // Check that spatial and temporal track grid size are odd - if (!bool(m_cfg.spatialTrkGridSize % 2) || - !bool(m_cfg.temporalTrkGridSize % 2)) { - throw std::invalid_argument( - "Please choose an odd spatialTrkGridSize and temporalTrkGridSize."); - } - } - - /// @brief Calculates the bin center from the bin number - /// @param bin Bin number - /// @param binExtent Bin extent - /// @return Bin center - static float getBinCenter(int bin, float binExtent); - - /// @brief Calculates the bin number corresponding to a d, z, or time value - /// @param value d, z, or time value - /// @param binExtent Bin extent - /// @return Bin number - static int getBin(float value, float binExtent); - - /// @brief Finds the maximum density of a DensityMap - /// @param densityMap Map between bins and corresponding density - /// values - /// @return Iterator of the map entry with the highest density - DensityMap::const_iterator highestDensityEntry( - const DensityMap& densityMap) const; + AdaptiveGridTrackDensity(const Config& cfg); /// @brief Returns the z and t coordinate of maximum (surrounding) /// track density @@ -145,7 +128,68 @@ class AdaptiveGridTrackDensity { void subtractTrack(const DensityMap& trackDensityMap, DensityMap& mainDensityMap) const; + // TODO this should not be public + /// @brief Calculates the bin center from the bin number + /// @param bin Bin number + /// @param binExtent Bin extent + /// @return Bin center + static double getBinCenter(std::int32_t bin, double binExtent); + private: + Config m_cfg; + + /// @brief Calculates the bin number corresponding to a d, z, or time value + /// @param value d, z, or time value + /// @param binExtent Bin extent + /// @return Bin number + static std::int32_t getBin(double value, double binExtent); + + /// @brief Calculates the grid size in z or time direction + /// @param sigma Standard deviation of the track density + /// @param trkSigmas Number of standard deviations that the grid + /// covers in z or time direction + /// @param binExtent Bin extent + /// @param trkGridSizeRange Optional minimal and maximal number of bins + /// in z or time direction + /// @return Grid size + static std::uint32_t getTrkGridSize(double sigma, double trkSigmas, + double binExtent, + const GridSizeRange& trkGridSizeRange); + + /// @brief Calculates the bin number corresponding to a z value + /// @param value z value + /// @return Bin number + std::int32_t getSpatialBin(double value) const; + /// @brief Calculates the bin number corresponding to a time value + /// @param value Time value + /// @return Bin number + std::int32_t getTemporalBin(double value) const; + + /// @brief Calculates the spatial bin center corresponding to a bin number + /// @param bin Bin number + /// @return Bin center + double getSpatialBinCenter(std::int32_t bin) const; + /// @brief Calculates the temporal bin center corresponding to a bin number + /// @param bin Bin number + /// @return Bin center + double getTemporalBinCenter(std::int32_t bin) const; + + /// @brief Calculates the grid size in z direction + /// @param sigma Standard deviation of the track density + /// @return Grid size + std::uint32_t getSpatialTrkGridSize(double sigma) const; + /// @brief Calculates the grid size in time direction + /// @param sigma Standard deviation of the track density + /// @return Grid size + std::uint32_t getTemporalTrkGridSize(double sigma) const; + + /// @brief Finds the maximum density of a DensityMap + /// @param densityMap Map between bins and corresponding density + /// values + /// @return Iterator of the map entry with the highest density + DensityMap::const_iterator highestDensityEntry( + const DensityMap& densityMap) const; + /// @brief Function that creates a track density map, i.e., a map from bins /// to the corresponding density values for a single track. /// @@ -155,7 +199,9 @@ class AdaptiveGridTrackDensity { /// @param cov 3x3 impact parameter covariance matrix DensityMap createTrackGrid(const Acts::Vector3& impactParams, const Bin& centralBin, - const Acts::SquareMatrix3& cov) const; + const Acts::SquareMatrix3& cov, + std::uint32_t spatialTrkGridSize, + std::uint32_t temporalTrkGridSize) const; /// @brief Function that estimates the seed width in z direction based /// on the full width at half maximum (FWHM) of the maximum density peak @@ -167,22 +213,8 @@ class AdaptiveGridTrackDensity { /// @param maxZT z-t position of the maximum density value /// /// @return The width - Result estimateSeedWidth(const DensityMap& densityMap, - const ZTPosition& maxZT) const; - - /// @brief Helper to retrieve values of an nDim-dimensional normal - /// distribution - /// @note The constant prefactor (2 * pi)^(- nDim / 2) is discarded - /// - /// @param args Coordinates where the Gaussian should be evaluated - /// @note args must be in a coordinate system with origin at the mean - /// values of the Gaussian - /// @param cov Covariance matrix - /// - /// @return Multivariate Gaussian evaluated at args - template - static float multivariateGaussian(const Acts::ActsVector& args, - const Acts::ActsSquareMatrix& cov); + Result estimateSeedWidth(const DensityMap& densityMap, + const ZTPosition& maxZT) const; /// @brief Checks (up to) first three density maxima that have a /// maximum relative deviation of 'relativeDensityDev' from the @@ -201,9 +233,7 @@ class AdaptiveGridTrackDensity { /// @param bin Bin whose neighbors in z we want to sum up /// /// @return The density sum - float getDensitySum(const DensityMap& densityMap, const Bin& bin) const; - - Config m_cfg; + double getDensitySum(const DensityMap& densityMap, const Bin& bin) const; }; } // namespace Acts diff --git a/Core/src/Detector/DetectorVolume.cpp b/Core/src/Detector/DetectorVolume.cpp index 2b49b3476a0..c233f36ae3b 100644 --- a/Core/src/Detector/DetectorVolume.cpp +++ b/Core/src/Detector/DetectorVolume.cpp @@ -240,7 +240,7 @@ void Acts::Experimental::DetectorVolume::updateNavigationState( const GeometryContext& gctx, NavigationState& nState) const { nState.currentVolume = this; m_surfaceCandidatesUpdater(gctx, nState); - nState.surfaceCandidate = nState.surfaceCandidates.begin(); + nState.surfaceCandidateIndex = 0; } void Acts::Experimental::DetectorVolume::assignSurfaceCandidatesUpdater( diff --git a/Core/src/EventData/VectorTrackContainer.cpp b/Core/src/EventData/VectorTrackContainer.cpp index 1b75a5eb987..27f7425182a 100644 --- a/Core/src/EventData/VectorTrackContainer.cpp +++ b/Core/src/EventData/VectorTrackContainer.cpp @@ -34,6 +34,7 @@ VectorTrackContainerBase::VectorTrackContainerBase( for (const auto& [key, value] : other.m_dynamic) { m_dynamic.insert({key, value->clone()}); } + m_dynamicKeys = other.m_dynamicKeys; assert(checkConsistency()); } } // namespace detail_vtc diff --git a/Core/src/Geometry/Extent.cpp b/Core/src/Geometry/Extent.cpp index 175607c7abb..c60b4fb80bb 100644 --- a/Core/src/Geometry/Extent.cpp +++ b/Core/src/Geometry/Extent.cpp @@ -84,7 +84,7 @@ void Acts::Extent::addConstrain(const Acts::Extent& rhs, for (const auto& bValue : s_binningValues) { if (rhs.constrains(bValue) && !constrains(bValue)) { const auto& cRange = rhs.range(bValue); - m_range[bValue].setMin(cRange.min() + envelope[bValue][0u]); + m_range[bValue].setMin(cRange.min() - envelope[bValue][0u]); m_range[bValue].setMax(cRange.max() + envelope[bValue][1u]); m_constrains.set(bValue); } diff --git a/Core/src/Geometry/GenericApproachDescriptor.cpp b/Core/src/Geometry/GenericApproachDescriptor.cpp index 9307cde34d6..a6f5eca3ade 100644 --- a/Core/src/Geometry/GenericApproachDescriptor.cpp +++ b/Core/src/Geometry/GenericApproachDescriptor.cpp @@ -25,16 +25,16 @@ void Acts::GenericApproachDescriptor::registerLayer(const Layer& lay) { Acts::SurfaceIntersection Acts::GenericApproachDescriptor::approachSurface( const GeometryContext& gctx, const Vector3& position, - const Vector3& direction, const BoundaryCheck& bcheck, double pLimit, - double oLimit, double tolerance) const { + const Vector3& direction, const BoundaryCheck& bcheck, double nearLimit, + double farLimit) const { // almost always 2 - boost::container::small_vector sIntersections; + boost::container::small_vector sIntersections; sIntersections.reserve(m_surfaceCache.size()); for (const auto& sf : m_surfaceCache) { auto sfIntersection = sf->intersect(gctx, position, direction, bcheck); for (const auto& intersection : sfIntersection.split()) { if (intersection && - detail::checkIntersection(intersection, pLimit, oLimit, tolerance)) { + detail::checkIntersection(intersection, nearLimit, farLimit)) { sIntersections.push_back(intersection); } } diff --git a/Core/src/Geometry/Layer.cpp b/Core/src/Geometry/Layer.cpp index 705188d3dae..bbb2210ade1 100644 --- a/Core/src/Geometry/Layer.cpp +++ b/Core/src/Geometry/Layer.cpp @@ -122,8 +122,8 @@ Acts::Layer::compatibleSurfaces( // (0) End surface check // @todo: - we might be able to skip this by use of options.pathLimit // check if you have to stop at the endSurface - double pathLimit = options.pathLimit; - double overstepLimit = options.overstepLimit; + double nearLimit = options.nearLimit; + double farLimit = options.farLimit; if (options.endObject != nullptr) { // intersect the end surface // - it is the final one don't use the boundary check at all @@ -136,7 +136,7 @@ Acts::Layer::compatibleSurfaces( // -> do not return compatible surfaces since they may lead you on a wrong // navigation path if (endInter) { - pathLimit = endInter.pathLength(); + farLimit = endInter.pathLength(); } else { return sIntersections; } @@ -147,7 +147,7 @@ Acts::Layer::compatibleSurfaces( // -> this avoids punch through for cylinders double pCorrection = surfaceRepresentation().pathCorrection(gctx, position, direction); - pathLimit = 1.5 * thickness() * pCorrection; + farLimit = 1.5 * thickness() * pCorrection; } // lemma 0 : accept the surface @@ -188,8 +188,7 @@ Acts::Layer::compatibleSurfaces( for (const auto& sfi : sfmi.split()) { // check if intersection is valid and pathLimit has not been exceeded if (sfi && - detail::checkIntersection(sfi.intersection(), pathLimit, - overstepLimit, s_onSurfaceTolerance)) { + detail::checkIntersection(sfi.intersection(), nearLimit, farLimit)) { sIntersections.push_back(sfi); } } @@ -272,17 +271,15 @@ Acts::SurfaceIntersection Acts::Layer::surfaceOnApproach( (surfaceRepresentation().surfaceMaterial() != nullptr)); // The Limits: current path & overstepping - double pLimit = options.pathLimit; - double oLimit = options.overstepLimit; - // TODO this should be configurable - double tolerance = s_onSurfaceTolerance; + double nearLimit = options.nearLimit; + double farLimit = options.farLimit; // Helper function to find valid intersection auto findValidIntersection = [&](const SurfaceMultiIntersection& sfmi) -> SurfaceIntersection { for (const auto& sfi : sfmi.split()) { - if (sfi && detail::checkIntersection(sfi.intersection(), pLimit, oLimit, - tolerance)) { + if (sfi && + detail::checkIntersection(sfi.intersection(), nearLimit, farLimit)) { return sfi; } } @@ -294,8 +291,7 @@ Acts::SurfaceIntersection Acts::Layer::surfaceOnApproach( // Approach descriptor present and resolving is necessary if (m_approachDescriptor && (resolvePS || resolveMS)) { SurfaceIntersection aSurface = m_approachDescriptor->approachSurface( - gctx, position, direction, options.boundaryCheck, pLimit, oLimit, - tolerance); + gctx, position, direction, options.boundaryCheck, nearLimit, farLimit); return aSurface; } diff --git a/Core/src/Geometry/TrackingVolume.cpp b/Core/src/Geometry/TrackingVolume.cpp index bc282f551b2..1fd2187c159 100644 --- a/Core/src/Geometry/TrackingVolume.cpp +++ b/Core/src/Geometry/TrackingVolume.cpp @@ -482,8 +482,8 @@ Acts::TrackingVolume::compatibleBoundaries( boost::container::small_vector bIntersections; // The Limits: current, path & overstepping - double pLimit = options.pathLimit; - double oLimit = 0; + double nearLimit = 0; + double farLimit = options.farLimit; // Helper function to test intersection auto checkIntersection = @@ -496,28 +496,26 @@ Acts::TrackingVolume::compatibleBoundaries( if (options.forceIntersectBoundaries) { const bool coCriterion = - std::abs(sIntersection.pathLength()) < std::abs(oLimit); + std::abs(sIntersection.pathLength()) < std::abs(nearLimit); ACTS_VERBOSE("Forcing intersection with surface " << bSurface->surfaceRepresentation().geometryId()); if (coCriterion) { ACTS_VERBOSE("Intersection forced successfully "); ACTS_VERBOSE("- intersection path length " << std::abs(sIntersection.pathLength()) - << " < overstep limit " << std::abs(oLimit)); + << " < overstep limit " << std::abs(nearLimit)); return BoundaryIntersection(sIntersection, bSurface); } ACTS_VERBOSE("Can't force intersection: "); ACTS_VERBOSE("- intersection path length " << std::abs(sIntersection.pathLength()) - << " > overstep limit " << std::abs(oLimit)); + << " > overstep limit " << std::abs(nearLimit)); } ACTS_VERBOSE("Check intersection with surface " << bSurface->surfaceRepresentation().geometryId()); - if (detail::checkIntersection( - sIntersection.intersection(), pLimit, oLimit, - s_onSurfaceTolerance, logger)) { + if (detail::checkIntersection(sIntersection.intersection(), nearLimit, + farLimit, logger)) { return BoundaryIntersection(sIntersection, bSurface); } } @@ -614,10 +612,9 @@ Acts::TrackingVolume::compatibleLayers( auto atIntersection = tLayer->surfaceOnApproach(gctx, position, direction, options); auto path = atIntersection.pathLength(); - bool withinLimit = std::abs(path) <= std::abs(options.pathLimit); + bool withinLimit = std::abs(path) <= std::abs(options.farLimit); // Intersection is ok - take it (move to surface on approach) - if (atIntersection && - (atIntersection.object() != options.targetSurface) && withinLimit) { + if (atIntersection && withinLimit) { // create a layer intersection lIntersections.push_back(LayerIntersection(atIntersection, tLayer)); } @@ -676,8 +673,8 @@ Acts::TrackingVolume::compatibleSurfacesFromHierarchy( sIntersections.reserve(20); // arbitrary // The limits for this navigation step - double pLimit = options.pathLimit; - double oLimit = options.overstepLimit; + double nearLimit = options.nearLimit; + double farLimit = options.farLimit; if (m_bvhTop == nullptr) { return sIntersections; @@ -703,7 +700,7 @@ Acts::TrackingVolume::compatibleSurfacesFromHierarchy( auto sfmi = srf.intersect(gctx, position, direction, BoundaryCheck(false)); for (const auto& sfi : sfmi.split()) { - if (sfi && sfi.pathLength() > oLimit && sfi.pathLength() <= pLimit) { + if (sfi && detail::checkIntersection(sfi, nearLimit, farLimit)) { sIntersections.push_back(sfi); } } diff --git a/Core/src/Propagator/detail/CovarianceEngine.cpp b/Core/src/Propagator/detail/CovarianceEngine.cpp index df6788d54f3..e196f24a97f 100644 --- a/Core/src/Propagator/detail/CovarianceEngine.cpp +++ b/Core/src/Propagator/detail/CovarianceEngine.cpp @@ -16,6 +16,7 @@ #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp" #include "Acts/EventData/detail/TransformationBoundToFree.hpp" #include "Acts/EventData/detail/TransformationFreeToBound.hpp" +#include "Acts/Propagator/detail/JacobianEngine.hpp" #include "Acts/Utilities/AlgebraHelpers.hpp" #include "Acts/Utilities/JacobianHelpers.hpp" #include "Acts/Utilities/Result.hpp" @@ -408,5 +409,48 @@ void transportCovarianceToCurvilinear(BoundSquareMatrix& boundCovariance, boundToFreeJacobian, direction); } +Acts::Result boundToBoundConversion( + const GeometryContext& gctx, const BoundTrackParameters& boundParameters, + const Surface& targetSurface, const Vector3& bField) { + const auto& sourceSurface = boundParameters.referenceSurface(); + + Acts::FreeVector freePars = Acts::detail::transformBoundToFreeParameters( + sourceSurface, gctx, boundParameters.parameters()); + + auto res = Acts::detail::transformFreeToBoundParameters(freePars, + targetSurface, gctx); + + if (!res.ok()) { + return res.error(); + } + Acts::BoundVector parOut = *res; + + std::optional covOut = std::nullopt; + + if (boundParameters.covariance().has_value()) { + Acts::BoundToFreeMatrix boundToFreeJacobian = + sourceSurface.boundToFreeJacobian(gctx, boundParameters.parameters()); + + Acts::FreeMatrix freeTransportJacobian = FreeMatrix::Identity(); + + FreeVector freeToPathDerivatives = FreeVector::Zero(); + freeToPathDerivatives.head<3>() = freePars.segment<3>(eFreeDir0); + + freeToPathDerivatives.segment<3>(eFreeDir0) = + bField.cross(freePars.segment<3>(eFreeDir0)); + + BoundMatrix boundToBoundJac = detail::boundToBoundTransportJacobian( + gctx, freePars, boundToFreeJacobian, freeTransportJacobian, + freeToPathDerivatives, targetSurface); + + covOut = boundToBoundJac * (*boundParameters.covariance()) * + boundToBoundJac.transpose(); + } + + return Acts::BoundTrackParameters{targetSurface.getSharedPtr(), parOut, + covOut, + boundParameters.particleHypothesis()}; +} + } // namespace detail } // namespace Acts diff --git a/Core/src/Vertexing/AdaptiveGridTrackDensity.cpp b/Core/src/Vertexing/AdaptiveGridTrackDensity.cpp index 1491e188870..bff02fa8910 100644 --- a/Core/src/Vertexing/AdaptiveGridTrackDensity.cpp +++ b/Core/src/Vertexing/AdaptiveGridTrackDensity.cpp @@ -5,6 +5,7 @@ // 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 "Acts/Vertexing/AdaptiveGridTrackDensity.hpp" #include "Acts/Utilities/AlgebraHelpers.hpp" @@ -12,27 +13,132 @@ #include -float Acts::AdaptiveGridTrackDensity::getBinCenter(int bin, float binExtent) { +namespace Acts { + +namespace { + +/// @brief Helper to retrieve values of an nDim-dimensional normal +/// distribution +/// @note The constant prefactor (2 * pi)^(- nDim / 2) is discarded +/// +/// @param args Coordinates where the Gaussian should be evaluated +/// @note args must be in a coordinate system with origin at the mean +/// values of the Gaussian +/// @param cov Covariance matrix +/// +/// @return Multivariate Gaussian evaluated at args +template +double multivariateGaussian(const ActsVector& args, + const ActsSquareMatrix& cov) { + double exponent = -0.5 * args.transpose().dot(cov.inverse() * args); + double gaussianDensity = safeExp(exponent) / std::sqrt(cov.determinant()); + return gaussianDensity; +} + +} // namespace + +double AdaptiveGridTrackDensity::getBinCenter(std::int32_t bin, + double binExtent) { return bin * binExtent; } -int Acts::AdaptiveGridTrackDensity::getBin(float value, float binExtent) { +std::int32_t AdaptiveGridTrackDensity::getBin(double value, double binExtent) { return static_cast(std::floor(value / binExtent - 0.5) + 1); } -typename Acts::AdaptiveGridTrackDensity::DensityMap::const_iterator -Acts::AdaptiveGridTrackDensity::highestDensityEntry( +std::uint32_t AdaptiveGridTrackDensity::getTrkGridSize( + double sigma, double nTrkSigmas, double binExtent, + const GridSizeRange& trkGridSizeRange) { + std::uint32_t size = + static_cast(std::ceil(2 * nTrkSigmas * sigma / binExtent)); + // Make sure the grid size is odd + if (size % 2 == 0) { + ++size; + } + // Make sure the grid size is within the allowed range + if (trkGridSizeRange.first && size < trkGridSizeRange.first.value()) { + size = trkGridSizeRange.first.value(); + } + if (trkGridSizeRange.second && size > trkGridSizeRange.second.value()) { + size = trkGridSizeRange.second.value(); + } + return size; +} + +std::int32_t AdaptiveGridTrackDensity::getSpatialBin(double value) const { + return getBin(value, m_cfg.spatialBinExtent); +} + +std::int32_t AdaptiveGridTrackDensity::getTemporalBin(double value) const { + if (!m_cfg.useTime) { + return 0; + } + return getBin(value, m_cfg.temporalBinExtent); +} + +double AdaptiveGridTrackDensity::getSpatialBinCenter(std::int32_t bin) const { + return getBinCenter(bin, m_cfg.spatialBinExtent); +} + +double AdaptiveGridTrackDensity::getTemporalBinCenter(std::int32_t bin) const { + if (!m_cfg.useTime) { + return 0.; + } + return getBinCenter(bin, m_cfg.temporalBinExtent); +} + +std::uint32_t AdaptiveGridTrackDensity::getSpatialTrkGridSize( + double sigma) const { + return getTrkGridSize(sigma, m_cfg.nSpatialTrkSigmas, m_cfg.spatialBinExtent, + m_cfg.spatialTrkGridSizeRange); +} + +std::uint32_t AdaptiveGridTrackDensity::getTemporalTrkGridSize( + double sigma) const { + if (!m_cfg.useTime) { + return 1; + } + return getTrkGridSize(sigma, m_cfg.nTemporalTrkSigmas, + m_cfg.temporalBinExtent, + m_cfg.temporalTrkGridSizeRange); +} + +AdaptiveGridTrackDensity::AdaptiveGridTrackDensity(const Config& cfg) + : m_cfg(cfg) { + if (m_cfg.spatialTrkGridSizeRange.first && + m_cfg.spatialTrkGridSizeRange.first.value() % 2 == 0) { + throw std::invalid_argument( + "AdaptiveGridTrackDensity: spatialTrkGridSizeRange.first must be odd"); + } + if (m_cfg.spatialTrkGridSizeRange.second && + m_cfg.spatialTrkGridSizeRange.second.value() % 2 == 0) { + throw std::invalid_argument( + "AdaptiveGridTrackDensity: spatialTrkGridSizeRange.second must be odd"); + } + if (m_cfg.temporalTrkGridSizeRange.first && + m_cfg.temporalTrkGridSizeRange.first.value() % 2 == 0) { + throw std::invalid_argument( + "AdaptiveGridTrackDensity: temporalTrkGridSizeRange.first must be odd"); + } + if (m_cfg.temporalTrkGridSizeRange.second && + m_cfg.temporalTrkGridSizeRange.second.value() % 2 == 0) { + throw std::invalid_argument( + "AdaptiveGridTrackDensity: temporalTrkGridSizeRange.second must be " + "odd"); + } +} + +AdaptiveGridTrackDensity::DensityMap::const_iterator +AdaptiveGridTrackDensity::highestDensityEntry( const DensityMap& densityMap) const { auto maxEntry = std::max_element( std::begin(densityMap), std::end(densityMap), - [](const auto& densityEntry1, const auto& densityEntry2) { - return densityEntry1.second < densityEntry2.second; - }); + [](const auto& a, const auto& b) { return a.second < b.second; }); return maxEntry; } -Acts::Result -Acts::AdaptiveGridTrackDensity::getMaxZTPosition(DensityMap& densityMap) const { +Result +AdaptiveGridTrackDensity::getMaxZTPosition(DensityMap& densityMap) const { if (densityMap.empty()) { return VertexingError::EmptyInput; } @@ -46,22 +152,15 @@ Acts::AdaptiveGridTrackDensity::getMaxZTPosition(DensityMap& densityMap) const { bin = highestDensitySumBin(densityMap); } - // Derive corresponding z value - float maxZ = getBinCenter(bin.first, m_cfg.spatialBinExtent); + // Derive corresponding z and t value + double maxZ = getSpatialBinCenter(bin.first); + double maxT = getTemporalBinCenter(bin.second); - ZTPosition maxValues = std::make_pair(maxZ, 0.); - - // Get t value of the maximum if we do time vertex seeding - if (m_cfg.temporalTrkGridSize > 1) { - float maxT = getBinCenter(bin.second, m_cfg.temporalBinExtent); - maxValues.second = maxT; - } - - return maxValues; + return std::make_pair(maxZ, maxT); } -Acts::Result -Acts::AdaptiveGridTrackDensity::getMaxZTPositionAndWidth( +Result +AdaptiveGridTrackDensity::getMaxZTPositionAndWidth( DensityMap& densityMap) const { // Get z value where the density is the highest auto maxZTRes = getMaxZTPosition(densityMap); @@ -75,219 +174,205 @@ Acts::AdaptiveGridTrackDensity::getMaxZTPositionAndWidth( if (!widthRes.ok()) { return widthRes.error(); } - float width = *widthRes; + double width = *widthRes; ZTPositionAndWidth maxZTAndWidth{maxZT, width}; return maxZTAndWidth; } -typename Acts::AdaptiveGridTrackDensity::DensityMap -Acts::AdaptiveGridTrackDensity::addTrack(const Acts::BoundTrackParameters& trk, - DensityMap& mainDensityMap) const { +AdaptiveGridTrackDensity::DensityMap AdaptiveGridTrackDensity::addTrack( + const BoundTrackParameters& trk, DensityMap& mainDensityMap) const { ActsVector<3> impactParams = trk.impactParameters(); ActsSquareMatrix<3> cov = trk.impactParameterCovariance().value(); + std::uint32_t spatialTrkGridSize = + getSpatialTrkGridSize(std::sqrt(cov(1, 1))); + std::uint32_t temporalTrkGridSize = + getTemporalTrkGridSize(std::sqrt(cov(2, 2))); + // Calculate bin in d direction - int centralDBin = getBin(impactParams(0), m_cfg.spatialBinExtent); + std::int32_t centralDBin = getBin(impactParams(0), m_cfg.spatialBinExtent); // Check if current track affects grid density - if (std::abs(centralDBin) > (m_cfg.spatialTrkGridSize - 1) / 2.) { - DensityMap emptyTrackDensityMap; - return emptyTrackDensityMap; + if (std::abs(centralDBin) > (spatialTrkGridSize - 1) / 2.) { + // Return empty map + return {}; } - // Calculate bin in z direction - int centralZBin = getBin(impactParams(1), m_cfg.spatialBinExtent); - // If we don't do time vertex seeding, the time index is set to 0 - Bin centralBin = std::make_pair(centralZBin, 0.); + // Calculate bin in z and t direction + std::int32_t centralZBin = getSpatialBin(impactParams(1)); + std::int32_t centralTBin = getTemporalBin(impactParams(2)); - // Calculate bin in t direction if we do time vertex seeding - if (m_cfg.temporalTrkGridSize > 1) { - int centralTBin = getBin(impactParams(2), m_cfg.temporalBinExtent); - centralBin.second = centralTBin; - } + Bin centralBin = {centralZBin, centralTBin}; - DensityMap trackDensityMap = createTrackGrid(impactParams, centralBin, cov); + DensityMap trackDensityMap = createTrackGrid( + impactParams, centralBin, cov, spatialTrkGridSize, temporalTrkGridSize); - for (const auto& densityEntry : trackDensityMap) { - Bin bin = densityEntry.first; - float trackDensity = densityEntry.second; - // Check if z bin is already part of the main grid - if (mainDensityMap.count(bin) == 1) { - mainDensityMap.at(bin) += trackDensity; - } else { - mainDensityMap[bin] = trackDensity; - } + for (const auto& [bin, density] : trackDensityMap) { + mainDensityMap[bin] += density; } return trackDensityMap; } -void Acts::AdaptiveGridTrackDensity::subtractTrack( - const DensityMap& trackDensityMap, DensityMap& mainDensityMap) const { - for (auto it = trackDensityMap.begin(); it != trackDensityMap.end(); it++) { - mainDensityMap.at(it->first) -= it->second; +void AdaptiveGridTrackDensity::subtractTrack(const DensityMap& trackDensityMap, + DensityMap& mainDensityMap) const { + for (const auto& [bin, density] : trackDensityMap) { + mainDensityMap[bin] -= density; } } -typename Acts::AdaptiveGridTrackDensity::DensityMap -Acts::AdaptiveGridTrackDensity::createTrackGrid( - const Acts::Vector3& impactParams, const Bin& centralBin, - const Acts::SquareMatrix3& cov) const { +AdaptiveGridTrackDensity::DensityMap AdaptiveGridTrackDensity::createTrackGrid( + const Vector3& impactParams, const Bin& centralBin, + const SquareMatrix3& cov, std::uint32_t spatialTrkGridSize, + std::uint32_t temporalTrkGridSize) const { DensityMap trackDensityMap; - int halfSpatialTrkGridSize = (m_cfg.spatialTrkGridSize - 1) / 2; - int firstZBin = centralBin.first - halfSpatialTrkGridSize; + std::uint32_t halfSpatialTrkGridSize = (spatialTrkGridSize - 1) / 2; + std::int32_t firstZBin = centralBin.first - halfSpatialTrkGridSize; // If we don't do time vertex seeding, firstTBin will be 0. - int halfTemporalTrkGridSize = (m_cfg.temporalTrkGridSize - 1) / 2; - int firstTBin = centralBin.second - halfTemporalTrkGridSize; + std::uint32_t halfTemporalTrkGridSize = (temporalTrkGridSize - 1) / 2; + std::int32_t firstTBin = centralBin.second - halfTemporalTrkGridSize; // Loop over bins - for (unsigned int i = 0; i < m_cfg.temporalTrkGridSize; i++) { - int tBin = firstTBin + i; - // If we don't do vertex time seeding, we set the time to 0 since it will be - // discarded in the for loop below anyways - float t = 0; - if (m_cfg.temporalTrkGridSize > 1) { - t = getBinCenter(tBin, m_cfg.temporalBinExtent); + for (std::uint32_t i = 0; i < temporalTrkGridSize; i++) { + std::int32_t tBin = firstTBin + i; + double t = getTemporalBinCenter(tBin); + if (t < m_cfg.temporalWindow.first || t > m_cfg.temporalWindow.second) { + continue; } - for (unsigned int j = 0; j < m_cfg.spatialTrkGridSize; j++) { - int zBin = firstZBin + j; - float z = getBinCenter(zBin, m_cfg.spatialBinExtent); + for (std::uint32_t j = 0; j < spatialTrkGridSize; j++) { + std::int32_t zBin = firstZBin + j; + double z = getSpatialBinCenter(zBin); + if (z < m_cfg.spatialWindow.first || z > m_cfg.spatialWindow.second) { + continue; + } // Bin coordinates in the d-z-t plane - Acts::Vector3 binCoords(0., z, t); + Vector3 binCoords(0., z, t); // Transformation to coordinate system with origin at the track center binCoords -= impactParams; - Bin bin = std::make_pair(zBin, tBin); - if (m_cfg.temporalTrkGridSize == 1) { - trackDensityMap[bin] = multivariateGaussian<2>( - binCoords.head<2>(), cov.topLeftCorner<2, 2>()); + Bin bin = {zBin, tBin}; + double density = 0; + if (m_cfg.useTime) { + density = multivariateGaussian<3>(binCoords, cov); } else { - trackDensityMap[bin] = multivariateGaussian<3>(binCoords, cov); + density = multivariateGaussian<2>(binCoords.head<2>(), + cov.topLeftCorner<2, 2>()); + } + // Only add density if it is positive (otherwise it is 0) + if (density > 0) { + trackDensityMap[bin] = density; } } } + return trackDensityMap; } -Acts::Result Acts::AdaptiveGridTrackDensity::estimateSeedWidth( +Result AdaptiveGridTrackDensity::estimateSeedWidth( const DensityMap& densityMap, const ZTPosition& maxZT) const { if (densityMap.empty()) { return VertexingError::EmptyInput; } - // Get z bin of max density - int zMaxBin = getBin(maxZT.first, m_cfg.spatialBinExtent); - int tMaxBin = 0; - // Fill the time bin with a non-zero value if we do time vertex seeding - if (m_cfg.temporalTrkGridSize > 1) { - tMaxBin = getBin(maxZT.second, m_cfg.temporalBinExtent); - } - const float maxValue = densityMap.at(std::make_pair(zMaxBin, tMaxBin)); + // Get z and t bin of max density + std::int32_t zMaxBin = getBin(maxZT.first, m_cfg.spatialBinExtent); + std::int32_t tMaxBin = getBin(maxZT.second, m_cfg.temporalBinExtent); + double maxValue = densityMap.at({zMaxBin, tMaxBin}); - int rhmBin = zMaxBin; - float gridValue = maxValue; + std::int32_t rhmBin = zMaxBin; + double gridValue = maxValue; // Boolean indicating whether we find a filled bin that has a densityValue <= // maxValue/2 bool binFilled = true; while (gridValue > maxValue / 2) { // Check if we are still operating on continuous z values - if (densityMap.count(std::make_pair(rhmBin + 1, tMaxBin)) == 0) { + if (densityMap.count({rhmBin + 1, tMaxBin}) == 0) { binFilled = false; break; } rhmBin += 1; - gridValue = densityMap.at(std::make_pair(rhmBin, tMaxBin)); + gridValue = densityMap.at({rhmBin, tMaxBin}); } // Use linear approximation to find better z value for FWHM between bins - float rightDensity = 0; + double rightDensity = 0; if (binFilled) { - rightDensity = densityMap.at(std::make_pair(rhmBin, tMaxBin)); + rightDensity = densityMap.at({rhmBin, tMaxBin}); } - float leftDensity = densityMap.at(std::make_pair(rhmBin - 1, tMaxBin)); - float deltaZ1 = m_cfg.spatialBinExtent * (maxValue / 2 - leftDensity) / - (rightDensity - leftDensity); + double leftDensity = densityMap.at({rhmBin - 1, tMaxBin}); + double deltaZ1 = m_cfg.spatialBinExtent * (maxValue / 2 - leftDensity) / + (rightDensity - leftDensity); - int lhmBin = zMaxBin; + std::int32_t lhmBin = zMaxBin; gridValue = maxValue; binFilled = true; while (gridValue > maxValue / 2) { // Check if we are still operating on continuous z values - if (densityMap.count(std::make_pair(lhmBin - 1, tMaxBin)) == 0) { + if (densityMap.count({lhmBin - 1, tMaxBin}) == 0) { binFilled = false; break; } lhmBin -= 1; - gridValue = densityMap.at(std::make_pair(lhmBin, tMaxBin)); + gridValue = densityMap.at({lhmBin, tMaxBin}); } // Use linear approximation to find better z value for FWHM between bins - rightDensity = densityMap.at(std::make_pair(lhmBin + 1, tMaxBin)); + rightDensity = densityMap.at({lhmBin + 1, tMaxBin}); if (binFilled) { - leftDensity = densityMap.at(std::make_pair(lhmBin, tMaxBin)); + leftDensity = densityMap.at({lhmBin, tMaxBin}); } else { leftDensity = 0; } - float deltaZ2 = m_cfg.spatialBinExtent * (rightDensity - maxValue / 2) / - (rightDensity - leftDensity); + double deltaZ2 = m_cfg.spatialBinExtent * (rightDensity - maxValue / 2) / + (rightDensity - leftDensity); - float fwhm = (rhmBin - lhmBin) * m_cfg.spatialBinExtent - deltaZ1 - deltaZ2; + double fwhm = (rhmBin - lhmBin) * m_cfg.spatialBinExtent - deltaZ1 - deltaZ2; // FWHM = 2.355 * sigma - float width = fwhm / 2.355f; + double width = fwhm / 2.355; - return std::isnormal(width) ? width : 0.0f; -} - -template -float Acts::AdaptiveGridTrackDensity::multivariateGaussian( - const Acts::ActsVector& args, - const Acts::ActsSquareMatrix& cov) { - float expo = -0.5 * args.transpose().dot(cov.inverse() * args); - float gaussianDensity = safeExp(expo) / std::sqrt(cov.determinant()); - return gaussianDensity; + return std::isnormal(width) ? width : 0.0; } -typename Acts::AdaptiveGridTrackDensity::Bin -Acts::AdaptiveGridTrackDensity::highestDensitySumBin( +AdaptiveGridTrackDensity::Bin AdaptiveGridTrackDensity::highestDensitySumBin( DensityMap& densityMap) const { // The global maximum auto firstMax = highestDensityEntry(densityMap); Bin binFirstMax = firstMax->first; - float valueFirstMax = firstMax->second; - float firstSum = getDensitySum(densityMap, binFirstMax); + double valueFirstMax = firstMax->second; + double firstSum = getDensitySum(densityMap, binFirstMax); // Smaller maxima must have a density of at least: // valueFirstMax - densityDeviation - float densityDeviation = valueFirstMax * m_cfg.maxRelativeDensityDev; + double densityDeviation = valueFirstMax * m_cfg.maxRelativeDensityDev; // Get the second highest maximum - densityMap.at(binFirstMax) = 0; + densityMap[binFirstMax] = 0; auto secondMax = highestDensityEntry(densityMap); Bin binSecondMax = secondMax->first; - float valueSecondMax = secondMax->second; - float secondSum = 0; + double valueSecondMax = secondMax->second; + double secondSum = 0; if (valueFirstMax - valueSecondMax < densityDeviation) { secondSum = getDensitySum(densityMap, binSecondMax); } else { // If the second maximum is not sufficiently large the third maximum won't // be either - densityMap.at(binFirstMax) = valueFirstMax; + densityMap[binFirstMax] = valueFirstMax; return binFirstMax; } // Get the third highest maximum - densityMap.at(binSecondMax) = 0; + densityMap[binSecondMax] = 0; auto thirdMax = highestDensityEntry(densityMap); Bin binThirdMax = thirdMax->first; - float valueThirdMax = thirdMax->second; - float thirdSum = 0; + double valueThirdMax = thirdMax->second; + double thirdSum = 0; if (valueFirstMax - valueThirdMax < densityDeviation) { thirdSum = getDensitySum(densityMap, binThirdMax); } // Revert back to original values - densityMap.at(binFirstMax) = valueFirstMax; - densityMap.at(binSecondMax) = valueSecondMax; + densityMap[binFirstMax] = valueFirstMax; + densityMap[binSecondMax] = valueSecondMax; // Return the z bin position of the highest density sum if (secondSum > firstSum && secondSum > thirdSum) { @@ -299,25 +384,20 @@ Acts::AdaptiveGridTrackDensity::highestDensitySumBin( return binFirstMax; } -float Acts::AdaptiveGridTrackDensity::getDensitySum( - const DensityMap& densityMap, const Bin& bin) const { +double AdaptiveGridTrackDensity::getDensitySum(const DensityMap& densityMap, + const Bin& bin) const { + auto valueOrZero = [&densityMap](const Bin& b) { + if (auto it = densityMap.find(b); it != densityMap.end()) { + return it->second; + } + return 0.0f; + }; + // Add density from the bin. - float sum = densityMap.at(bin); // Check if neighboring bins are part of the densityMap and add them (if they // are not part of the map, we assume them to be 0). - // Note that each key in a map is unique; the .count() function therefore - // returns either 0 or 1. - Bin binShifted = bin; - // Add density from the neighboring bin in -z direction. - binShifted.first -= 1; - if (densityMap.count(binShifted) == 1) { - sum += densityMap.at(binShifted); - } - - // Add density from the neighboring bin in +z direction. - binShifted.first += 2; - if (densityMap.count(binShifted) == 1) { - sum += densityMap.at(binShifted); - } - return sum; + return valueOrZero(bin) + valueOrZero({bin.first - 1, bin.second}) + + valueOrZero({bin.first + 1, bin.second}); } + +} // namespace Acts diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp index aedede24d2a..2b06e55799d 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp @@ -11,6 +11,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Direction.hpp" #include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/ProxyAccessor.hpp" #include "Acts/EventData/TrackContainer.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" #include "Acts/EventData/VectorTrackContainer.hpp" @@ -135,7 +136,7 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( tracks.addColumn("trackGroup"); tracksTemp.addColumn("trackGroup"); - Acts::TrackAccessor seedNumber("trackGroup"); + Acts::ProxyAccessor seedNumber("trackGroup"); unsigned int nSeed = 0; diff --git a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp index 3e5a7f1bdfb..493b5799d14 100644 --- a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp +++ b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp @@ -8,6 +8,7 @@ #include "ActsExamples/TrackFindingExaTrkX/TrackFindingFromPrototrackAlgorithm.hpp" +#include "Acts/EventData/ProxyAccessor.hpp" #include "ActsExamples/EventData/IndexSourceLink.hpp" #include "ActsExamples/EventData/MeasurementCalibration.hpp" @@ -125,7 +126,7 @@ ActsExamples::ProcessCode TrackFindingFromPrototrackAlgorithm::execute( TrackContainer tracks(trackContainer, trackStateContainer); tracks.addColumn("trackGroup"); - Acts::TrackAccessor seedNumber("trackGroup"); + Acts::ProxyAccessor seedNumber("trackGroup"); std::size_t nSeed = 0; std::size_t nFailed = 0; diff --git a/Examples/Algorithms/Utilities/src/TracksToTrajectories.cpp b/Examples/Algorithms/Utilities/src/TracksToTrajectories.cpp index 40a6f60cbba..99835202091 100644 --- a/Examples/Algorithms/Utilities/src/TracksToTrajectories.cpp +++ b/Examples/Algorithms/Utilities/src/TracksToTrajectories.cpp @@ -10,6 +10,7 @@ #include "Acts/EventData/GenericBoundTrackParameters.hpp" #include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/ProxyAccessor.hpp" #include "Acts/EventData/TrackContainer.hpp" #include "Acts/EventData/TrackProxy.hpp" #include "Acts/Surfaces/Surface.hpp" @@ -37,7 +38,7 @@ ProcessCode TracksToTrajectories::execute(const AlgorithmContext& ctx) const { TrajectoriesContainer trajectories; trajectories.reserve(tracks.size()); - static const Acts::ConstTrackAccessor seedNumber("trackGroup"); + static const Acts::ConstProxyAccessor seedNumber("trackGroup"); if (tracks.hasColumn(Acts::hashString("trackGroup"))) { // track group by seed is available, produce grouped trajectories diff --git a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp index c4e828a534c..1e1f736832d 100644 --- a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp +++ b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp @@ -112,4 +112,5 @@ class AdaptiveMultiVertexFinderAlgorithm final : public IAlgorithm { WriteDataHandle m_outputVertices{this, "OutputVertices"}; }; + } // namespace ActsExamples diff --git a/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp b/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp index a21c39c3f92..cf0a11d4ca2 100644 --- a/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp +++ b/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp @@ -64,14 +64,11 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute( } else if (m_cfg.seedFinder == SeedFinder::AdaptiveGridSeeder) { // Set up track density used during vertex seeding Acts::AdaptiveGridTrackDensity::Config trkDensityCfg; - // Number of bins in z-direction - trkDensityCfg.spatialTrkGridSize = 9; // Bin extent in z-direction - trkDensityCfg.spatialBinExtent = 0.015; - // Number of bins in t-direction - trkDensityCfg.temporalTrkGridSize = 5; + trkDensityCfg.spatialBinExtent = 15. * Acts::UnitConstants::um; // Bin extent in t-direction - trkDensityCfg.temporalBinExtent = 19.; + trkDensityCfg.temporalBinExtent = 19. * Acts::UnitConstants::mm; + trkDensityCfg.useTime = m_cfg.useTime; Acts::AdaptiveGridTrackDensity trkDensity(trkDensityCfg); // Set up vertex seeder and finder @@ -128,7 +125,7 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::executeAfterSeederChoice( // Set the initial variance of the 4D vertex position. Since time is on a // numerical scale, we have to provide a greater value in the corresponding // dimension. - finderConfig.initialVariances << 1e+2, 1e+2, 1e+2, 1e+7; + finderConfig.initialVariances << 1e+2, 1e+2, 1e+2, 1e+8; finderConfig.tracksMaxZinterval = 1. * Acts::UnitConstants::mm; finderConfig.maxIterations = 200; finderConfig.useTime = m_cfg.useTime; diff --git a/Examples/Algorithms/Vertexing/src/IterativeVertexFinderAlgorithm.cpp b/Examples/Algorithms/Vertexing/src/IterativeVertexFinderAlgorithm.cpp index 0052381886e..c432da57910 100644 --- a/Examples/Algorithms/Vertexing/src/IterativeVertexFinderAlgorithm.cpp +++ b/Examples/Algorithms/Vertexing/src/IterativeVertexFinderAlgorithm.cpp @@ -10,7 +10,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Propagator/EigenStepper.hpp" -#include "Acts/Propagator/detail/VoidPropagatorComponents.hpp" +#include "Acts/Propagator/VoidNavigator.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/Result.hpp" #include "Acts/Vertexing/IterativeVertexFinder.hpp" @@ -71,9 +71,8 @@ ActsExamples::ProcessCode ActsExamples::IterativeVertexFinderAlgorithm::execute( Acts::EigenStepper<> stepper(m_cfg.bField); // Set up propagator with void navigator - auto propagator = - std::make_shared(stepper, Acts::detail::VoidNavigator{}, - logger().cloneWithSuffix("Propagator")); + auto propagator = std::make_shared( + stepper, Acts::VoidNavigator{}, logger().cloneWithSuffix("Propagator")); // Setup the vertex fitter Fitter::Config vertexFitterCfg; Fitter vertexFitter(vertexFitterCfg, diff --git a/Examples/Algorithms/Vertexing/src/VertexFitterAlgorithm.cpp b/Examples/Algorithms/Vertexing/src/VertexFitterAlgorithm.cpp index 7889a8f0ef6..05ea2f092f2 100644 --- a/Examples/Algorithms/Vertexing/src/VertexFitterAlgorithm.cpp +++ b/Examples/Algorithms/Vertexing/src/VertexFitterAlgorithm.cpp @@ -11,7 +11,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/MagneticField/MagneticFieldProvider.hpp" #include "Acts/Propagator/EigenStepper.hpp" -#include "Acts/Propagator/detail/VoidPropagatorComponents.hpp" +#include "Acts/Propagator/VoidNavigator.hpp" #include "Acts/Utilities/Result.hpp" #include "Acts/Vertexing/Vertex.hpp" #include "ActsExamples/EventData/ProtoVertex.hpp" @@ -45,7 +45,7 @@ ActsExamples::ProcessCode ActsExamples::VertexFitterAlgorithm::execute( // Setup the propagator with void navigator auto propagator = std::make_shared( - stepper, Acts::detail::VoidNavigator{}, logger().cloneWithSuffix("Prop")); + stepper, Acts::VoidNavigator{}, logger().cloneWithSuffix("Prop")); PropagatorOptions propagatorOpts(ctx.geoContext, ctx.magFieldContext); // Setup the vertex fitter VertexFitter::Config vertexFitterCfg; diff --git a/Examples/Io/Csv/src/CsvTrackWriter.cpp b/Examples/Io/Csv/src/CsvTrackWriter.cpp index d6b3b336ce4..f56f9f8761f 100644 --- a/Examples/Io/Csv/src/CsvTrackWriter.cpp +++ b/Examples/Io/Csv/src/CsvTrackWriter.cpp @@ -10,6 +10,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/ProxyAccessor.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" #include "Acts/Utilities/Helpers.hpp" #include "Acts/Utilities/MultiIndex.hpp" @@ -109,7 +110,7 @@ ProcessCode CsvTrackWriter::writeT(const AlgorithmContext& context, nMajorityHits = particleHitCount.front().hitCount; } - static const Acts::ConstTrackAccessor seedNumber( + static const Acts::ConstProxyAccessor seedNumber( "trackGroup"); // track info diff --git a/Examples/Python/src/DD4hepComponent.cpp b/Examples/Python/src/DD4hepComponent.cpp index 3da8ba9799f..5199de27c73 100644 --- a/Examples/Python/src/DD4hepComponent.cpp +++ b/Examples/Python/src/DD4hepComponent.cpp @@ -6,10 +6,13 @@ // 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 "Acts/Detector/GeometryIdGenerator.hpp" #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp" #include "Acts/Plugins/DD4hep/DD4hepDetectorStructure.hpp" #include "Acts/Plugins/DD4hep/DD4hepFieldAdapter.hpp" +#include "Acts/Plugins/DD4hep/DD4hepIdentifierMapper.hpp" #include "Acts/Plugins/Python/Utilities.hpp" +#include "Acts/Utilities/Logger.hpp" #include "ActsExamples/DD4hepDetector/DD4hepDetector.hpp" #include "ActsExamples/DD4hepDetector/DD4hepGeometryService.hpp" #include "ActsExamples/Framework/IContextDecorator.hpp" @@ -69,6 +72,39 @@ PYBIND11_MODULE(ActsPythonBindingsDD4hep, m) { m, "DD4hepDetectorElement"); } + { + m.def("createDD4hepIdGeoIdMap", + [](const Acts::TrackingGeometry& tGeometry) + -> std::map { + // The surface visitor + struct DD4hepIdGrabber { + std::map + dd4hepIdGeoIdMap; + + void operator()(const Acts::Surface* surface) { + const auto* dde = surface->associatedDetectorElement(); + const auto* dd4hepDetElement = + dynamic_cast(dde); + // Check if it is valid + if (dd4hepDetElement) { + dd4hep::DDSegmentation::VolumeID dd4hepID = + dd4hepDetElement->sourceElement().volumeID(); + auto geoID = surface->geometryId(); + dd4hepIdGeoIdMap[dd4hepID] = geoID; + } + } + }; + + // Create an instance + DD4hepIdGrabber dd4hepIdGrabber; + // Visit the surfaces & return what you have + tGeometry.visitSurfaces(dd4hepIdGrabber); + return dd4hepIdGrabber.dd4hepIdGeoIdMap; + }); + } + { using Options = Acts::Experimental::DD4hepDetectorStructure::Options; auto o = py::class_(m, "DD4hepDetectorOptions").def(py::init<>()); @@ -78,6 +114,41 @@ PYBIND11_MODULE(ActsPythonBindingsDD4hep, m) { ACTS_PYTHON_STRUCT_END(); patchKwargsConstructor(o); + + m.def( + "attachDD4hepGeoIdMapper", + [](Acts::Experimental::DD4hepDetectorStructure::Options& options, + const std::map& dd4hepIdGeoIdMap) { + // The Geo mapper + auto geoIdMapper = + std::make_shared( + Acts::DD4hepIdentifierMapper::Config{dd4hepIdGeoIdMap}, + Acts::getDefaultLogger("GeometryIdMapper", options.logLevel)); + + // A remaining recursive logger + auto geoIdGenerator = + std::make_shared( + Acts::Experimental::GeometryIdGenerator::Config{}, + Acts::getDefaultLogger("GeometryIdGenerator", + options.logLevel)); + + std::tuple< + std::shared_ptr, + std::shared_ptr> + chainedGenerators = {geoIdGenerator, geoIdMapper}; + + auto chainedGeoIdGenerator = std::make_shared< + const Acts::Experimental::ChainedGeometryIdGenerator< + std::shared_ptr< + const Acts::Experimental::GeometryIdGenerator>, + std::shared_ptr>>( + std::move(chainedGenerators), + Acts::getDefaultLogger("ChainedGeometryIdGenerator", + options.logLevel)); + + options.geoIdGenerator = chainedGeoIdGenerator; + }); } { diff --git a/Examples/Python/src/Geometry.cpp b/Examples/Python/src/Geometry.cpp index f5d37c93996..e02768736f4 100644 --- a/Examples/Python/src/Geometry.cpp +++ b/Examples/Python/src/Geometry.cpp @@ -79,7 +79,8 @@ void addGeometry(Context& ctx) { .def("boundary", &Acts::GeometryIdentifier::boundary) .def("approach", &Acts::GeometryIdentifier::approach) .def("sensitive", &Acts::GeometryIdentifier::sensitive) - .def("extra", &Acts::GeometryIdentifier::extra); + .def("extra", &Acts::GeometryIdentifier::extra) + .def("value", &Acts::GeometryIdentifier::value); } { diff --git a/Examples/Scripts/Python/detector_creation.py b/Examples/Scripts/Python/detector_creation.py index 0cea0942b08..c10674d4703 100644 --- a/Examples/Scripts/Python/detector_creation.py +++ b/Examples/Scripts/Python/detector_creation.py @@ -8,6 +8,7 @@ DD4hepGeometryService, ) +import json from common import getOpenDataDetectorDirectory @@ -24,7 +25,28 @@ dd4hepGeometryService = DD4hepGeometryService(dd4hepConfig) dd4hepDetector = DD4hepDetector(dd4hepGeometryService) + cOptions = DD4hepDetectorOptions(logLevel=acts.logging.VERBOSE, emulateToGraph="") + + # Uncomment if you want to use the geometry id mapping + # This map can be produced with the 'geometry.py' script + geoIdMappingFile = "odd-dd4hep-geoid-mapping-wo-extra.json" + if geoIdMappingFile is not None: + # Load the geometry id mapping json file + with open(geoIdMappingFile) as f: + # load the file as is + geometry_id_mapping = json.load(f) + # create a dictionary with GeometryIdentifier as value + geometry_id_mapping_patched = { + int(k): acts.GeometryIdentifier(int(v)) + for k, v in geometry_id_mapping.items() + } + # patch the options struct + acts.examples.dd4hep.attachDD4hepGeoIdMapper( + cOptions, geometry_id_mapping_patched + ) + # Context and options geoContext = acts.GeometryContext() - cOptions = DD4hepDetectorOptions(logLevel=acts.logging.INFO, emulateToGraph="") [detector, contextors, store] = dd4hepDetector.finalize(geoContext, cOptions) + + acts.examples.writeDetectorToJsonDetray(geoContext, detector, "odd-detray") diff --git a/Examples/Scripts/Python/geometry.py b/Examples/Scripts/Python/geometry.py index a8b2c349c30..2fe61641fcc 100755 --- a/Examples/Scripts/Python/geometry.py +++ b/Examples/Scripts/Python/geometry.py @@ -16,6 +16,7 @@ ) import acts +import json from acts import MaterialMapJsonConverter @@ -30,7 +31,6 @@ def runGeometry( outputJson=True, outputRoot=True, ): - for ievt in range(events): eventStore = WhiteBoard(name=f"EventStore#{ievt}", level=acts.logging.INFO) ialg = 0 @@ -43,6 +43,8 @@ def runGeometry( raise RuntimeError("Failed to decorate event context") if outputCsv: + # if not os.path.isdir(outputDir + "/csv"): + # os.makedirs(outputDir + "/csv") writer = CsvTrackingGeometryWriter( level=acts.logging.INFO, trackingGeometry=trackingGeometry, @@ -58,6 +60,8 @@ def runGeometry( writer.write(context, trackingGeometry) if outputJson: + # if not os.path.isdir(outputDir + "/json"): + # os.makedirs(outputDir + "/json") writer = JsonSurfacesWriter( level=acts.logging.INFO, trackingGeometry=trackingGeometry, @@ -90,6 +94,17 @@ def runGeometry( if "__main__" == __name__: detector, trackingGeometry, decorators = AlignedDetector.create() # detector, trackingGeometry, decorators = GenericDetector.create() - # detector, trackingGeometry, decorators = getOpenDataDetector(getOpenDataDetectorDirectory()) + # detector, trackingGeometry, decorators = getOpenDataDetector( + # getOpenDataDetectorDirectory() + # ) runGeometry(trackingGeometry, decorators, outputDir=os.getcwd()) + + # Uncomment if you want to create the geometry id mapping for DD4hep + # dd4hepIdGeoIdMap = acts.examples.dd4hep.createDD4hepIdGeoIdMap(trackingGeometry) + # dd4hepIdGeoIdValueMap = {} + # for key, value in dd4hepIdGeoIdMap.items(): + # dd4hepIdGeoIdValueMap[key] = value.value() + + # with open('odd-dd4hep-geoid-mapping.json', 'w') as outfile: + # json.dump(dd4hepIdGeoIdValueMap, outfile) diff --git a/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp b/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp index 97c78cbd62a..583111fe2c0 100644 --- a/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp +++ b/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp @@ -38,6 +38,9 @@ class ISurfaceMaterial; /// class DD4hepDetectorElement : public TGeoDetectorElement { public: + // Define the context type + using DD4hepVolumeID = dd4hep::DDSegmentation::VolumeID; + /// Broadcast the context type using ContextType = GeometryContext; @@ -81,6 +84,9 @@ class DD4hepDetectorElement : public TGeoDetectorElement { ~DD4hepDetectorElement() override = default; + // Give access to the DD4hep detector element + const dd4hep::DetElement& sourceElement() const { return m_detElement; } + private: /// DD4hep detector element dd4hep::DetElement m_detElement; diff --git a/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepDetectorStructure.hpp b/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepDetectorStructure.hpp index d8eb088b4d8..abb482cc7fc 100644 --- a/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepDetectorStructure.hpp +++ b/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepDetectorStructure.hpp @@ -42,6 +42,8 @@ class DD4hepDetectorStructure { /// If this string is not empty, the detector is not built, but only /// a building graph is generated as `.dot` file with the given name std::string emulateToGraph = ""; + /// A Top level geometry id generator + std::shared_ptr geoIdGenerator = nullptr; }; /// Constructor with from file name diff --git a/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepIdentifierMapper.hpp b/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepIdentifierMapper.hpp new file mode 100644 index 00000000000..fc9c8ee3453 --- /dev/null +++ b/Plugins/DD4hep/include/Acts/Plugins/DD4hep/DD4hepIdentifierMapper.hpp @@ -0,0 +1,59 @@ +// 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/. + +#pragma once + +#include "Acts/Detector/GeometryIdMapper.hpp" +#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp" +#include "Acts/Surfaces/Surface.hpp" + +#include + +namespace Acts { + +using DD4hepIdentifier = DD4hepDetectorElement::DD4hepVolumeID; + +/// @brief This struct helps to capture the DD4hep identifier +/// from the DD4hepDetectorElement from a Acts::Surface +/// +/// It can be used for mapping/assigning geometry ids to surfaces +struct DD4hepIdentifierCapture { + /// @brief Return an invalid identifier for volumes as they are not directly translated + /// @return maximum value + DD4hepIdentifier operator()( + const Acts::Experimental::DetectorVolume& /*unused*/) const { + return std::numeric_limits::max(); + } + + /// @brief Return an invalid identifier for portal objects as they are not directly translated + /// @return maximum value + DD4hepIdentifier operator()( + const Acts::Experimental::Portal& /*unused*/) const { + return std::numeric_limits::max(); + } + + /// @brief Return the DD4hep identifier from the DD4hepDetectorElement associated to the surface + /// @param surface the Acts::Surface object + /// @return the DD4hep identifier (or maximum value if now DD4hepDetectorElement is associated) + DD4hepIdentifier operator()(const Acts::Surface& surface) const { + // Get the DD4hepDetectorElement + const auto* dde = surface.associatedDetectorElement(); + const auto* dd4hepDetElement = + dynamic_cast(dde); + // Check if it is valid + if (dd4hepDetElement) { + return dd4hepDetElement->sourceElement().volumeID(); + } // Return the maximum value + return std::numeric_limits::max(); + } +}; + +using DD4hepIdentifierMapper = + Experimental::GeometryIdMapper; + +} // namespace Acts diff --git a/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp b/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp index 82d1a3cd4ec..31194af124f 100644 --- a/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp +++ b/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp @@ -161,11 +161,11 @@ std::shared_ptr volumeBuilder_dd4hep( auto volumeHelper = cylinderVolumeHelper_dd4hep(logger); // create local logger for conversion ACTS_VERBOSE("Processing detector element: " << subDetector.name()); - dd4hep::DetType subDetType{subDetector.typeFlag()}; ACTS_VERBOSE("SubDetector type is: [" << subDetType << "], compound: " << (subDetector.type() == "compound" ? "yes" : "no")); + if (subDetector.type() == "compound") { ACTS_VERBOSE("Subdetector : '" << subDetector.name() << "' has type compound "); @@ -591,8 +591,13 @@ void collectSubDetectors_dd4hep(dd4hep::DetElement& detElement, dd4hep::DetElement childDetElement = child.second; dd4hep::DetType type{childDetElement.typeFlag()}; if (childDetElement.type() == "compound") { - subdetectors.push_back(childDetElement); - continue; + // Check if the compound is excluded from assembly + // This is needed to eventually exclude compounds of pixel, strip, etc. + // from barrel / endcap parsing + if (getParamOr("acts_legacy_assembly", childDetElement, true)) { + subdetectors.push_back(childDetElement); + continue; + } } if (type.is(dd4hep::DetType::TRACKER)) { diff --git a/Plugins/DD4hep/src/DD4hepBlueprintFactory.cpp b/Plugins/DD4hep/src/DD4hepBlueprintFactory.cpp index e97d83e7c8b..281a827aedc 100644 --- a/Plugins/DD4hep/src/DD4hepBlueprintFactory.cpp +++ b/Plugins/DD4hep/src/DD4hepBlueprintFactory.cpp @@ -228,8 +228,10 @@ Acts::Experimental::DD4hepBlueprintFactory::extractInternals( // Check if the extent should be measured auto interenalsMeasure = Acts::getParamOr( baseName + "_internals_measure", dd4hepElement, ""); - auto internalsClearance = Acts::getParamOr( - baseName + "_internals_clearance", dd4hepElement, 0.); + auto internalsClearance = + unitLength * + Acts::getParamOr(baseName + "_internals_clearance", + dd4hepElement, 0.); auto internalBinningValues = stringToBinningValues(interenalsMeasure); if (!internalBinningValues.empty()) { ACTS_VERBOSE(" - internals extent measurement requested"); @@ -238,6 +240,7 @@ Acts::Experimental::DD4hepBlueprintFactory::extractInternals( for (const auto& bv : internalBinningValues) { ACTS_VERBOSE(" -> measuring extent for " << binningValueNames()[bv]); + ACTS_VERBOSE(" -> with clearance :" << internalsClearance); clearance[bv] = {internalsClearance, internalsClearance}; } internalsExtent.setEnvelope(clearance); diff --git a/Plugins/DD4hep/src/DD4hepDetectorStructure.cpp b/Plugins/DD4hep/src/DD4hepDetectorStructure.cpp index d7dbc6c55ff..ae3683d5ca3 100644 --- a/Plugins/DD4hep/src/DD4hepDetectorStructure.cpp +++ b/Plugins/DD4hep/src/DD4hepDetectorStructure.cpp @@ -10,6 +10,7 @@ #include "Acts/Detector/CylindricalContainerBuilder.hpp" #include "Acts/Detector/DetectorBuilder.hpp" +#include "Acts/Detector/GeometryIdGenerator.hpp" #include "Acts/Detector/detail/BlueprintDrawer.hpp" #include "Acts/Detector/detail/BlueprintHelper.hpp" #include "Acts/Plugins/DD4hep/DD4hepBlueprintFactory.hpp" @@ -87,7 +88,9 @@ Acts::Experimental::DD4hepDetectorStructure::construct( "*** DD4hep : auto generated cylindrical detector builder ***"; dCfg.name = "Cylindrical detector from DD4hep blueprint"; dCfg.builder = detectorBuilder; - dCfg.geoIdGenerator = dd4hepBlueprint->geoIdGenerator; + dCfg.geoIdGenerator = options.geoIdGenerator != nullptr + ? options.geoIdGenerator + : dd4hepBlueprint->geoIdGenerator; detector = DetectorBuilder(dCfg, getDefaultLogger("DD4hepDetectorBuilder", options.logLevel)) .construct(gctx); diff --git a/Plugins/Json/src/DetectorVolumeJsonConverter.cpp b/Plugins/Json/src/DetectorVolumeJsonConverter.cpp index 8a5077768c8..8397c7a690d 100644 --- a/Plugins/Json/src/DetectorVolumeJsonConverter.cpp +++ b/Plugins/Json/src/DetectorVolumeJsonConverter.cpp @@ -120,12 +120,14 @@ nlohmann::json Acts::DetectorVolumeJsonConverter::toJsonDetray( int vIndex = findVolume(&volume, detectorVolumes); jVolume["index"] = vIndex; + std::size_t sIndex = 0; // Write the surfaces - patch bounds & augment with self links nlohmann::json jSurfaces; for (const auto& s : volume.surfaces()) { auto jSurface = SurfaceJsonConverter::toJsonDetray(gctx, *s, options.surfaceOptions); DetrayJsonHelper::addVolumeLink(jSurface["mask"], vIndex); + jSurface["index_in_coll"] = sIndex++; jSurfaces.push_back(jSurface); } @@ -140,7 +142,10 @@ nlohmann::json Acts::DetectorVolumeJsonConverter::toJsonDetray( (toJsonDetray(gctx, *p, ip, volume, orientedSurfaces, detectorVolumes, options.portalOptions)); std::for_each(jPortalSurfaces.begin(), jPortalSurfaces.end(), - [&](auto& jSurface) { jSurfaces.push_back(jSurface); }); + [&](auto& jSurface) { + jSurface["index_in_coll"] = sIndex++; + jSurfaces.push_back(jSurface); + }); } jVolume["surfaces"] = jSurfaces; diff --git a/Plugins/Json/src/PortalJsonConverter.cpp b/Plugins/Json/src/PortalJsonConverter.cpp index ae06f2ef26f..7f4af5aea17 100644 --- a/Plugins/Json/src/PortalJsonConverter.cpp +++ b/Plugins/Json/src/PortalJsonConverter.cpp @@ -92,7 +92,7 @@ std::vector Acts::PortalJsonConverter::toJsonDetray( const auto surfaceNormal = surface.normal(gctx, surfaceCenter); // Get the direction from the volume to the surface, correct link const auto volumeToSurface = surfaceCenter - volumeCenter; - if (volumeToSurface.dot(surfaceNormal) < 0) { + if (volumeToSurface.dot(surfaceNormal) < 0.) { outside = 0u; } } else { @@ -176,7 +176,8 @@ std::vector Acts::PortalJsonConverter::toJsonDetray( if (boundaries[ib - 1] >= clipRange[1u]) { break; } - if (highEdge <= clipRange[0u]) { + if (highEdge <= clipRange[0u] || + std::abs(highEdge - clipRange[0u]) < 1e-5) { continue; } if (highEdge > clipRange[1u]) { @@ -186,7 +187,6 @@ std::vector Acts::PortalJsonConverter::toJsonDetray( clippedBoundaries.push_back(highEdge); } } - // Interpret the clipped information // // Clipped cylinder case diff --git a/Plugins/Podio/include/Acts/Plugins/Podio/PodioTrackContainer.hpp b/Plugins/Podio/include/Acts/Plugins/Podio/PodioTrackContainer.hpp index 4ae77b59051..89cc7c06f46 100644 --- a/Plugins/Podio/include/Acts/Plugins/Podio/PodioTrackContainer.hpp +++ b/Plugins/Podio/include/Acts/Plugins/Podio/PodioTrackContainer.hpp @@ -294,6 +294,7 @@ class MutablePodioTrackContainer : public PodioTrackContainerBase { friend PodioTrackContainerBase; std::unique_ptr m_collection; + std::vector m_dynamicKeys; std::unordered_map> m_dynamic; @@ -383,6 +384,7 @@ class ConstPodioTrackContainer : public PodioTrackContainerBase { std::unordered_map> m_dynamic; + std::vector m_dynamicKeys; }; ACTS_STATIC_CHECK_CONCEPT(ConstTrackContainerBackend, ConstPodioTrackContainer); diff --git a/Plugins/Podio/include/Acts/Plugins/Podio/PodioTrackStateContainer.hpp b/Plugins/Podio/include/Acts/Plugins/Podio/PodioTrackStateContainer.hpp index 0ba4a0540eb..442b6d21590 100644 --- a/Plugins/Podio/include/Acts/Plugins/Podio/PodioTrackStateContainer.hpp +++ b/Plugins/Podio/include/Acts/Plugins/Podio/PodioTrackStateContainer.hpp @@ -339,6 +339,7 @@ class ConstPodioTrackStateContainer final std::unordered_map> m_dynamic; + std::vector m_dynamicKeys; }; static_assert(IsReadOnlyMultiTrajectory::value, @@ -655,15 +656,13 @@ class MutablePodioTrackStateContainer final std::unordered_map> m_dynamic; + std::vector m_dynamicKeys; }; static_assert( !IsReadOnlyMultiTrajectory::value, "MutablePodioTrackStateContainer should not be read-only"); -static_assert(!MutablePodioTrackStateContainer::ReadOnly, - "MutablePodioTrackStateContainer should not be read-only"); - ACTS_STATIC_CHECK_CONCEPT(MutableMultiTrajectoryBackend, MutablePodioTrackStateContainer); diff --git a/Tests/IntegrationTests/PropagationDenseConstant.cpp b/Tests/IntegrationTests/PropagationDenseConstant.cpp index db901b4f801..b5a9eee554d 100644 --- a/Tests/IntegrationTests/PropagationDenseConstant.cpp +++ b/Tests/IntegrationTests/PropagationDenseConstant.cpp @@ -80,7 +80,7 @@ inline Propagator makePropagator(double bz) { auto magField = std::make_shared(Acts::Vector3(0.0, 0.0, bz)); Stepper stepper(std::move(magField)); - auto logger = getDefaultLogger("Dense", Logging::INFO); + auto logger = getDefaultLogger("Nominal", Logging::INFO); return Propagator( std::move(stepper), Acts::Navigator{{makeDetector()}, logger->cloneWithSuffix("Nav")}, @@ -92,8 +92,12 @@ inline RiddersPropagator makeRiddersPropagator(double bz) { auto magField = std::make_shared(Acts::Vector3(0.0, 0.0, bz)); Stepper stepper(std::move(magField)); - return RiddersPropagator(std::move(stepper), - Acts::Navigator{{makeDetector()}}); + auto logger = getDefaultLogger("Ridder", Logging::INFO); + auto propagator = Propagator( + std::move(stepper), + Acts::Navigator{{makeDetector()}, logger->cloneWithSuffix("Nav")}, + logger->cloneWithSuffix("Prop")); + return RiddersPropagator(std::move(propagator)); } } // namespace diff --git a/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp b/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp index a0c875ec6e7..9aa704bd21f 100644 --- a/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp +++ b/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp @@ -180,8 +180,9 @@ BOOST_AUTO_TEST_CASE(CuboidWithCuboid) { outerBox->updateNavigationState(tContext, nState); - // We should have 12 candidates, 6 inner, 6 outer portals - BOOST_CHECK_EQUAL(nState.surfaceCandidates.size(), 12u); + // We should have 12 candidates, 6 inner, 6 outer portals but only 3 are + // reachable + BOOST_CHECK_EQUAL(nState.surfaceCandidates.size(), 3u); } BOOST_AUTO_TEST_CASE(CylinderWithSurfacesTestExtractors) { diff --git a/Tests/UnitTests/Core/Detector/GeometryIdGeneratorTests.cpp b/Tests/UnitTests/Core/Detector/GeometryIdGeneratorTests.cpp index 3f1ae0c59b1..3bb4ccf137a 100644 --- a/Tests/UnitTests/Core/Detector/GeometryIdGeneratorTests.cpp +++ b/Tests/UnitTests/Core/Detector/GeometryIdGeneratorTests.cpp @@ -77,6 +77,55 @@ std::vector> createVolumes( } } // namespace +/// @brief Test struct to increment the layer id by one +struct GeoIdIncrementer : public IGeometryIdGenerator { + struct Cache {}; + + /// @brief Interface method to generate a geometry id cache + /// @return a geometry id cache wrapped in a std::any object + IGeometryIdGenerator::GeoIdCache generateCache() const final { + // Unfold the tuple and add the attachers + return Cache{}; + } + + /// @brief Method for assigning a geometry id to a detector volume + /// + /// @param cache is the cache object for e.g. object counting + /// @param dVolume the detector volume to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& /*unused*/, + DetectorVolume& dVolume) const final { + auto vgid = dVolume.geometryId(); + vgid.setVolume(vgid.volume() + 1); + dVolume.assignGeometryId(vgid); + } + + /// @brief Method for assigning a geometry id to a portal + /// + /// @param cache is the cache object for e.g. object counting + /// @param portal the portal to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& /*unused*/, + Portal& portal) const final { + auto pgid = portal.surface().geometryId(); + pgid.setBoundary(pgid.boundary() + 1); + portal.surface().assignGeometryId(pgid); + } + + /// @brief Method for assigning a geometry id to a surface + /// + /// @param cache is the cache object for e.g. object counting + /// @param surface the surface to assign the geometry id to + void assignGeometryId(IGeometryIdGenerator::GeoIdCache& /*unused*/, + Surface& surface) const final { + auto sgid = surface.geometryId(); + if (sgid.sensitive() != 0u) { + sgid.setSensitive(sgid.sensitive() + 1); + } else { + sgid.setPassive(sgid.passive() + 1); + } + surface.assignGeometryId(sgid); + } +}; + BOOST_AUTO_TEST_SUITE(Detector) BOOST_AUTO_TEST_CASE(SequentialGeoIdGeneratorReset) { @@ -174,4 +223,39 @@ BOOST_AUTO_TEST_CASE(ContainerGeoIdGenerator) { BOOST_CHECK_EQUAL(volumes[2]->geometryId().layer(), 3); } +BOOST_AUTO_TEST_CASE(ChainedGeoIdGenerator) { + std::vector> detectorStore; + + auto volumes = createVolumes(detectorStore); + + GeometryIdGenerator::Config cfg; + + cfg.containerMode = true; + cfg.containerId = 15; + auto cgenerator = std::make_shared( + cfg, getDefaultLogger("ContainerIdGenerator", Logging::VERBOSE)); + + auto igenerator = std::make_shared(); + + std::tuple, + std::shared_ptr> + geoGenerators = {cgenerator, igenerator}; + + ChainedGeometryIdGenerator, + std::shared_ptr> + generator(std::move(geoGenerators)); + + auto cache = generator.generateCache(); + for (auto& volume : volumes) { + generator.assignGeometryId(cache, *volume); + } + + BOOST_CHECK_EQUAL(volumes[0]->geometryId().volume(), 16); + BOOST_CHECK_EQUAL(volumes[0]->geometryId().layer(), 1); + BOOST_CHECK_EQUAL(volumes[1]->geometryId().volume(), 16); + BOOST_CHECK_EQUAL(volumes[1]->geometryId().layer(), 2); + BOOST_CHECK_EQUAL(volumes[2]->geometryId().volume(), 16); + BOOST_CHECK_EQUAL(volumes[2]->geometryId().layer(), 3); +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/Tests/UnitTests/Core/EventData/MultiTrajectoryTests.cpp b/Tests/UnitTests/Core/EventData/MultiTrajectoryTests.cpp index 422d31c3879..7ecf486080a 100644 --- a/Tests/UnitTests/Core/EventData/MultiTrajectoryTests.cpp +++ b/Tests/UnitTests/Core/EventData/MultiTrajectoryTests.cpp @@ -14,6 +14,7 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/ProxyAccessor.hpp" #include "Acts/EventData/TrackParameters.hpp" #include "Acts/EventData/TrackStatePropMask.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" @@ -41,6 +42,7 @@ using namespace Acts; using namespace Acts::UnitLiterals; using namespace Acts::Test; using namespace Acts::detail::Test; +using namespace Acts::HashedStringLiteral; namespace bd = boost::unit_test::data; using ParametersVector = BoundTrackParameters::ParametersVector; @@ -75,7 +77,7 @@ BOOST_AUTO_TEST_CASE(ConstCorrectness) { VectorMultiTrajectory t; auto i0 = t.addTrackState(); - BOOST_CHECK(!t.ReadOnly); + BOOST_CHECK(!IsReadOnlyMultiTrajectory::value); { VectorMultiTrajectory::TrackStateProxy tsp = t.getTrackState(i0); @@ -244,4 +246,31 @@ BOOST_AUTO_TEST_CASE(MemoryStats) { } } +BOOST_AUTO_TEST_CASE(Accessors) { + VectorMultiTrajectory mtj; + mtj.addColumn("ndof"); + mtj.addColumn("super_chi2"); + + auto ts = mtj.getTrackState(mtj.addTrackState()); + + ProxyAccessor ndof("ndof"); + ConstProxyAccessor ndofConst("ndof"); + ProxyAccessor superChi2("super_chi2"); + ConstProxyAccessor superChi2Const("super_chi2"); + + ndof(ts) = 65; + BOOST_CHECK_EQUAL((ts.component()), 65); + BOOST_CHECK_EQUAL(ndofConst(ts), 65); + + // should not compile + // ndofConst(ts) = 66; + + superChi2(ts) = 123.45; + BOOST_CHECK_EQUAL((ts.component()), 123.45); + BOOST_CHECK_EQUAL(superChi2Const(ts), 123.45); + + // should not compile + // superChi2Const(ts) = 66.66; +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/Tests/UnitTests/Core/EventData/TrackContainerComplianceTests.cpp b/Tests/UnitTests/Core/EventData/TrackContainerComplianceTests.cpp index a4e18501ef5..cc6bbd782cc 100644 --- a/Tests/UnitTests/Core/EventData/TrackContainerComplianceTests.cpp +++ b/Tests/UnitTests/Core/EventData/TrackContainerComplianceTests.cpp @@ -7,6 +7,7 @@ // file, You can obtain one at http://mozilla.org/MPL/2.0/. #include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/ProxyAccessor.hpp" #include "Acts/EventData/TrackContainer.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" #include "Acts/EventData/VectorTrackContainer.hpp" @@ -73,7 +74,7 @@ ACTS_DOES_NOT_COMPILE_SUITE_BEGIN(BuildFromConstRef) auto t = tc.getTrack(tc.addTrack()); (void)t; - ConstTrackAccessor caccNMeasuements("nMeasurements"); + ConstProxyAccessor caccNMeasuements("nMeasurements"); ACTS_DOES_NOT_COMPILE_BEGIN(ConstAccessorMutate) caccNMeasuements(t) = 66; ACTS_DOES_NOT_COMPILE_END() diff --git a/Tests/UnitTests/Core/EventData/TrackTests.cpp b/Tests/UnitTests/Core/EventData/TrackTests.cpp index bab9f803586..f2c1afe8679 100644 --- a/Tests/UnitTests/Core/EventData/TrackTests.cpp +++ b/Tests/UnitTests/Core/EventData/TrackTests.cpp @@ -13,6 +13,7 @@ #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/ProxyAccessor.hpp" #include "Acts/EventData/TrackContainer.hpp" #include "Acts/EventData/TrackHelpers.hpp" #include "Acts/EventData/TrackProxy.hpp" @@ -260,8 +261,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Build, factory_t, holder_types) { t.setReferenceSurface(surface); BOOST_CHECK_EQUAL(surface.get(), &t.referenceSurface()); - TrackAccessor accNMeasuements("nMeasurements"); - ConstTrackAccessor caccNMeasuements("nMeasurements"); + ProxyAccessor accNMeasuements("nMeasurements"); + ConstProxyAccessor caccNMeasuements("nMeasurements"); t.nMeasurements() = 42; BOOST_CHECK_EQUAL(t2.nMeasurements(), 42); diff --git a/Tests/UnitTests/Core/Geometry/GenericApproachDescriptorTests.cpp b/Tests/UnitTests/Core/Geometry/GenericApproachDescriptorTests.cpp index 154831b2e12..daa7a2944c3 100644 --- a/Tests/UnitTests/Core/Geometry/GenericApproachDescriptorTests.cpp +++ b/Tests/UnitTests/Core/Geometry/GenericApproachDescriptorTests.cpp @@ -66,9 +66,8 @@ BOOST_AUTO_TEST_CASE(GenericApproachDescriptorProperties) { }; Vector3 zDir{0., 0., 1.}; BoundaryCheck bcheck{true}; - double pLimit = std::numeric_limits::max(); - double oLimit = -100 * UnitConstants::um; - double tolerance = s_onSurfaceTolerance; + double nearLimit = -100 * UnitConstants::um; + double farLimit = std::numeric_limits::max(); // std::vector> someSurfaces{ Surface::makeShared(), Surface::makeShared()}; @@ -78,7 +77,7 @@ BOOST_AUTO_TEST_CASE(GenericApproachDescriptorProperties) { BOOST_CHECK_NO_THROW(approachDescriptor.registerLayer(aLayer)); // approachSurface SurfaceIntersection surfIntersection = approachDescriptor.approachSurface( - tgContext, origin, zDir, bcheck, pLimit, oLimit, tolerance); + tgContext, origin, zDir, bcheck, nearLimit, farLimit); double expectedIntersection = 20.0; // property of SurfaceStub CHECK_CLOSE_REL(surfIntersection.pathLength(), expectedIntersection, 1e-6); // containedSurfaces() @@ -98,9 +97,8 @@ BOOST_AUTO_TEST_CASE(GenericApproachNoOverstepping) { Vector3 origin{0., -0.5, 1.}; Vector3 direction{0., 1., 0.}; BoundaryCheck bcheck{true}; - double pLimit = std::numeric_limits::max(); - double oLimit = -100 * UnitConstants::um; - double tolerance = s_onSurfaceTolerance; + double nearLimit = -100 * UnitConstants::um; + double farLimit = std::numeric_limits::max(); auto conCyl = Surface::makeShared(Transform3::Identity(), 10., 20.); @@ -110,7 +108,7 @@ BOOST_AUTO_TEST_CASE(GenericApproachNoOverstepping) { GenericApproachDescriptor gad(approachSurface); auto sfIntersection = gad.approachSurface( - GeometryContext(), origin, direction, bcheck, pLimit, oLimit, tolerance); + GeometryContext(), origin, direction, bcheck, nearLimit, farLimit); // No overstepping allowed, the preferred solution should be the forward one CHECK_CLOSE_ABS(sfIntersection.pathLength(), 10.5, s_epsilon); diff --git a/Tests/UnitTests/Core/Navigation/MultiWireNavigationTests.cpp b/Tests/UnitTests/Core/Navigation/MultiWireNavigationTests.cpp index 9b2b8e048ef..84712cfb47b 100644 --- a/Tests/UnitTests/Core/Navigation/MultiWireNavigationTests.cpp +++ b/Tests/UnitTests/Core/Navigation/MultiWireNavigationTests.cpp @@ -92,8 +92,9 @@ BOOST_AUTO_TEST_CASE(Navigation_in_Indexed_Surfaces) { nState.currentVolume = volumes.front().get(); nState.currentVolume->updateNavigationState(tContext, nState); - // check the surface candidates after update (12 surfaces + 6 portals) - BOOST_CHECK_EQUAL(nState.surfaceCandidates.size(), 18u); + // check the surface candidates after update (12 surfaces + 6 portals but only + // 5 are reachable) + BOOST_CHECK_EQUAL(nState.surfaceCandidates.size(), 5u); } BOOST_AUTO_TEST_SUITE_END() diff --git a/Tests/UnitTests/Core/Propagator/CovarianceEngineTests.cpp b/Tests/UnitTests/Core/Propagator/CovarianceEngineTests.cpp index 8a2d1fdc3cf..4dcee9f3d07 100644 --- a/Tests/UnitTests/Core/Propagator/CovarianceEngineTests.cpp +++ b/Tests/UnitTests/Core/Propagator/CovarianceEngineTests.cpp @@ -6,27 +6,45 @@ // 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 #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/ParticleHypothesis.hpp" #include "Acts/EventData/TrackParameters.hpp" #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp" +#include "Acts/EventData/detail/TransformationBoundToFree.hpp" #include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/MagneticField/ConstantBField.hpp" +#include "Acts/MagneticField/MagneticFieldContext.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/VoidNavigator.hpp" #include "Acts/Propagator/detail/CovarianceEngine.hpp" +#include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Surfaces/PlaneSurface.hpp" #include "Acts/Surfaces/Surface.hpp" +#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" #include "Acts/Utilities/Result.hpp" #include #include #include +#include #include #include +namespace bdata = boost::unit_test::data; + namespace Acts { namespace Test { +Acts::GeometryContext gctx; +Acts::MagneticFieldContext mctx; + +using namespace Acts::UnitLiterals; + using Covariance = BoundSquareMatrix; using Jacobian = BoundMatrix; @@ -164,5 +182,328 @@ BOOST_AUTO_TEST_CASE(covariance_engine_test) { BOOST_CHECK_NE(*(std::get<0>(boundResult).covariance()), Covariance::Identity()); } + +std::pair boundToBound(const BoundVector& parIn, + const BoundMatrix& covIn, + const Surface& srfA, + const Surface& srfB, + const Vector3& bField) { + Acts::BoundTrackParameters boundParamIn{srfA.getSharedPtr(), parIn, covIn, + ParticleHypothesis::pion()}; + + auto converted = + detail::boundToBoundConversion(gctx, boundParamIn, srfB, bField).value(); + + return {converted.parameters(), converted.covariance().value()}; +} + +using propagator_t = Propagator, VoidNavigator>; + +BoundVector localToLocal(const propagator_t& prop, const BoundVector& local, + const Surface& src, const Surface& dst) { + PropagatorOptions<> options{gctx, mctx}; + options.stepTolerance = 1e-10; + options.surfaceTolerance = 1e-10; + + BoundTrackParameters start{src.getSharedPtr(), local, std::nullopt, + ParticleHypothesis::pion()}; + + auto res = prop.propagate(start, dst, options).value(); + auto endParameters = res.endParameters.value(); + + BOOST_CHECK_EQUAL(&endParameters.referenceSurface(), &dst); + + BoundVector out = endParameters.parameters(); + out[eBoundTime] = local[eBoundTime]; + return out; +} + +propagator_t makePropagator(const Vector3& bField) { + return propagator_t{EigenStepper<>{std::make_shared(bField)}, + VoidNavigator{}}; +} + +BoundMatrix numericalBoundToBoundJacobian(const propagator_t& prop, + const BoundVector& parA, + const Surface& srfA, + const Surface& srfB) { + double h = 1e-4; + BoundMatrix J; + for (std::size_t i = 0; i < 6; i++) { + for (std::size_t j = 0; j < 6; j++) { + BoundVector parInitial1 = parA; + BoundVector parInitial2 = parA; + parInitial1[i] -= h; + parInitial2[i] += h; + BoundVector parFinal1 = localToLocal(prop, parInitial1, srfA, srfB); + BoundVector parFinal2 = localToLocal(prop, parInitial2, srfA, srfB); + + J(j, i) = (parFinal2[j] - parFinal1[j]) / (2 * h); + } + } + return J; +} + +unsigned int getNextSeed() { + static unsigned int seed = 10; + return ++seed; +} + +auto makeDist(double a, double b) { + return bdata::random( + (bdata::engine = std::mt19937{}, bdata::seed = getNextSeed(), + bdata::distribution = std::uniform_real_distribution(a, b))); +} + +const auto locDist = makeDist(-5_mm, 5_mm); +const auto bFieldDist = makeDist(0, 3_T); +const auto angleDist = makeDist(-2 * M_PI, 2 * M_PI); +const auto posDist = makeDist(-50_mm, 50_mm); + +#define MAKE_SURFACE() \ + [&]() { \ + Transform3 transformA = Transform3::Identity(); \ + transformA = AngleAxis3(Rx, Vector3::UnitX()); \ + transformA = AngleAxis3(Ry, Vector3::UnitY()); \ + transformA = AngleAxis3(Rz, Vector3::UnitZ()); \ + transformA.translation() << gx, gy, gz; \ + return Surface::makeShared(transformA); \ + }() + +BOOST_DATA_TEST_CASE(CovarianceConversionSamePlane, + (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^ + angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^ + locDist ^ locDist) ^ + bdata::xrange(100), + Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l0, l1, index) { + (void)index; + const Vector3 bField{Bx, By, Bz}; + + auto planeSurfaceA = MAKE_SURFACE(); + auto planeSurfaceB = + Surface::makeShared(planeSurfaceA->transform(gctx)); + + BoundMatrix covA; + covA.setZero(); + covA.diagonal() << 1, 2, 3, 4, 5, 6; + + BoundVector parA; + parA << l0, l1, M_PI / 4., M_PI_2 * 0.9, -1 / 1_GeV, 5_ns; + + // identical surface, this should work + auto [parB, covB] = + boundToBound(parA, covA, *planeSurfaceA, *planeSurfaceB, bField); + + // these should be the same because the plane surface are the same + CHECK_CLOSE_ABS(parA, parB, 1e-9); + CHECK_CLOSE_COVARIANCE(covA, covB, 1e-9); + + // now go back + auto [parA2, covA2] = + boundToBound(parB, covB, *planeSurfaceB, *planeSurfaceA, bField); + CHECK_CLOSE_ABS(parA, parA2, 1e-9); + CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-9); + + auto prop = makePropagator(bField); + + BoundMatrix J = + numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *planeSurfaceB); + BoundMatrix covC = J * covA * J.transpose(); + CHECK_CLOSE_COVARIANCE(covB, covC, 1e-6); +} + +BOOST_DATA_TEST_CASE(CovarianceConversionRotatedPlane, + (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^ + angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^ + locDist ^ locDist ^ angleDist) ^ + bdata::xrange(100), + Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l0, l1, angle, index) { + (void)index; + const Vector3 bField{Bx, By, Bz}; + + auto planeSurfaceA = MAKE_SURFACE(); + + Transform3 transform; + transform = planeSurfaceA->transform(gctx).rotation(); + transform = AngleAxis3(angle, planeSurfaceA->normal(gctx)) * transform; + transform.translation() = planeSurfaceA->transform(gctx).translation(); + auto planeSurfaceB = Surface::makeShared(transform); + + // sanity check that the normal didn't change + CHECK_CLOSE_ABS(planeSurfaceA->normal(gctx), planeSurfaceB->normal(gctx), + 1e-9); + + BoundMatrix covA; + covA.setZero(); + covA.diagonal() << 1, 2, 3, 4, 5, 6; + + BoundVector parA; + parA << l0, l1, M_PI / 4., M_PI_2 * 0.9, -1 / 1_GeV, 5_ns; + + auto [parB, covB] = + boundToBound(parA, covA, *planeSurfaceA, *planeSurfaceB, bField); + BoundVector exp = parA; + // loc0 and loc1 are rotated + exp.head<2>() = Eigen::Rotation2D(-angle) * parA.head<2>(); + + CHECK_CLOSE_ABS(exp, parB, 1e-9); + + // now go back + auto [parA2, covA2] = + boundToBound(parB, covB, *planeSurfaceB, *planeSurfaceA, bField); + CHECK_CLOSE_ABS(parA, parA2, 1e-9); + CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-9); + + auto prop = makePropagator(bField); + BoundMatrix J = + numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *planeSurfaceB); + BoundMatrix covC = J * covA * J.transpose(); + CHECK_CLOSE_COVARIANCE(covB, covC, 1e-6); +} + +BOOST_DATA_TEST_CASE(CovarianceConversionL0TiltedPlane, + (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^ + angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^ + locDist ^ angleDist) ^ + bdata::xrange(100), + Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l1, angle, index) { + (void)index; + const Vector3 bField{Bx, By, Bz}; + + auto planeSurfaceA = MAKE_SURFACE(); + + // make plane that is slightly rotated + Transform3 transform; + transform = planeSurfaceA->transform(gctx).rotation(); + + // figure out rotation axis along local x + Vector3 axis = planeSurfaceA->transform(gctx).rotation() * Vector3::UnitY(); + transform = AngleAxis3(angle, axis) * transform; + + transform.translation() = planeSurfaceA->transform(gctx).translation(); + + auto planeSurfaceB = Surface::makeShared(transform); + + BoundVector parA; + // loc 0 must be zero so we're on the intersection of both surfaces. + parA << 0, l1, M_PI / 4., M_PI_2 * 0.9, -1 / 1_GeV, 5_ns; + + BoundMatrix covA; + covA.setZero(); + covA.diagonal() << 1, 2, 3, 4, 5, 6; + + auto [parB, covB] = + boundToBound(parA, covA, *planeSurfaceA, *planeSurfaceB, bField); + + // now go back + auto [parA2, covA2] = + boundToBound(parB, covB, *planeSurfaceB, *planeSurfaceA, bField); + CHECK_CLOSE_ABS(parA, parA2, 1e-9); + CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-7); + + auto prop = makePropagator(bField); + BoundMatrix J = + numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *planeSurfaceB); + BoundMatrix covC = J * covA * J.transpose(); + CHECK_CLOSE_OR_SMALL((covB.template topLeftCorner<2, 2>()), + (covC.template topLeftCorner<2, 2>()), 1e-7, 1e-9); + CHECK_CLOSE_OR_SMALL(covB.diagonal(), covC.diagonal(), 1e-7, 1e-9); +} + +BOOST_DATA_TEST_CASE(CovarianceConversionL1TiltedPlane, + (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^ + angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^ + locDist ^ angleDist) ^ + bdata::xrange(100), + Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l0, angle, index) { + (void)index; + const Vector3 bField{Bx, By, Bz}; + + auto planeSurfaceA = MAKE_SURFACE(); + + // make plane that is slightly rotated + Transform3 transform; + transform = planeSurfaceA->transform(gctx).rotation(); + + Vector3 axis = planeSurfaceA->transform(gctx).rotation() * Vector3::UnitX(); + transform = AngleAxis3(angle, axis) * transform; + + transform.translation() = planeSurfaceA->transform(gctx).translation(); + + auto planeSurfaceB = Surface::makeShared(transform); + + BoundVector parA; + // loc 1 must be zero so we're on the intersection of both surfaces. + parA << l0, 0, M_PI / 4., M_PI_2 * 0.9, -1 / 1_GeV, 5_ns; + + BoundMatrix covA; + covA.setZero(); + covA.diagonal() << 1, 2, 3, 4, 5, 6; + + auto [parB, covB] = + boundToBound(parA, covA, *planeSurfaceA, *planeSurfaceB, bField); + + // now go back + auto [parA2, covA2] = + boundToBound(parB, covB, *planeSurfaceB, *planeSurfaceA, bField); + CHECK_CLOSE_ABS(parA, parA2, 1e-9); + // tolerance is a bit higher here + CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-6); + + auto prop = makePropagator(bField); + BoundMatrix J = + numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *planeSurfaceB); + BoundMatrix covC = J * covA * J.transpose(); + CHECK_CLOSE_OR_SMALL((covB.template topLeftCorner<2, 2>()), + (covC.template topLeftCorner<2, 2>()), 1e-6, 1e-9); + CHECK_CLOSE_OR_SMALL(covB.diagonal(), covC.diagonal(), 1e-6, 1e-9); +} + +BOOST_DATA_TEST_CASE(CovarianceConversionPerigee, + (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^ + angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^ + locDist ^ locDist ^ angleDist ^ angleDist ^ angleDist) ^ + bdata::xrange(100), + Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l0, l1, pRx, pRy, pRz, + index) { + (void)index; + const Vector3 bField{Bx, By, Bz}; + + auto planeSurfaceA = MAKE_SURFACE(); + + BoundVector parA; + parA << l0, l1, M_PI / 4., M_PI_2 * 0.9, -1 / 1_GeV, 5_ns; + + BoundMatrix covA; + covA.setZero(); + covA.diagonal() << 1, 2, 3, 4, 5, 6; + + Vector3 global = planeSurfaceA->localToGlobal(gctx, parA.head<2>()); + + Transform3 transform; + transform.setIdentity(); + transform.rotate(AngleAxis3(pRx, Vector3::UnitX())); + transform.rotate(AngleAxis3(pRy, Vector3::UnitY())); + transform.rotate(AngleAxis3(pRz, Vector3::UnitZ())); + transform.translation() = global; + + auto perigee = Surface::makeShared(transform); + + auto [parB, covB] = + boundToBound(parA, covA, *planeSurfaceA, *perigee, bField); + + // now go back + auto [parA2, covA2] = + boundToBound(parB, covB, *perigee, *planeSurfaceA, bField); + CHECK_CLOSE_ABS(parA, parA2, 1e-9); + CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-9); + + auto prop = makePropagator(bField); + BoundMatrix J = + numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *perigee); + BoundMatrix covC = J * covA * J.transpose(); + CHECK_CLOSE_ABS(covB, covC, 1e-7); +} + } // namespace Test } // namespace Acts diff --git a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp index 9d507d5dbc5..5af7d925dcd 100644 --- a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp +++ b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp @@ -458,15 +458,6 @@ BOOST_AUTO_TEST_CASE(Navigator_target_methods) { navCfg.resolvePassive = false; Navigator navigator{navCfg}; - // create a navigator for the Bounding Volume Hierarchy test - CubicBVHTrackingGeometry grid(20, 1000, 5); - Navigator::Config bvhNavCfg; - bvhNavCfg.trackingGeometry = grid.trackingGeometry; - bvhNavCfg.resolveSensitive = true; - bvhNavCfg.resolveMaterial = true; - bvhNavCfg.resolvePassive = false; - Navigator BVHNavigator{bvhNavCfg}; - // position and direction vector Vector4 position4(0., 0., 0, 0); Vector3 momentum(1., 1., 0); @@ -701,6 +692,17 @@ BOOST_AUTO_TEST_CASE(Navigator_target_methods) { std::cout << state.options.debugString << std::endl; state.options.debugString = ""; } +} + +BOOST_AUTO_TEST_CASE(Navigator_target_methods_BVH) { + // create a navigator for the Bounding Volume Hierarchy test + CubicBVHTrackingGeometry grid(20, 1000, 5); + Navigator::Config bvhNavCfg; + bvhNavCfg.trackingGeometry = grid.trackingGeometry; + bvhNavCfg.resolveSensitive = true; + bvhNavCfg.resolveMaterial = true; + bvhNavCfg.resolvePassive = false; + Navigator bvhNavigator{bvhNavCfg}; // test the navigation in a bounding volume hierarchy // ---------------------------------------------- @@ -710,37 +712,37 @@ BOOST_AUTO_TEST_CASE(Navigator_target_methods) { } // position and direction vector - Vector4 BVHPosition4(0., 0., 0, 0); - Vector3 BVHMomentum(1., 1., 0.); + Vector4 bvhPosition4(0., 0., 0, 0); + Vector3 bvhMomentum(1., 1., 0.); // the propagator cache - PropagatorState BVHState; - BVHState.options.debug = debug; + PropagatorState bvhState; + bvhState.options.debug = debug; // the stepper cache - BVHState.stepping.pos4 = BVHPosition4; - BVHState.stepping.dir = BVHMomentum.normalized(); + bvhState.stepping.pos4 = bvhPosition4; + bvhState.stepping.dir = bvhMomentum.normalized(); // Stepper - PropagatorState::Stepper BVHStepper; + PropagatorState::Stepper bvhStepper; - BVHNavigator.initialize(BVHState, BVHStepper); + bvhNavigator.initialize(bvhState, bvhStepper); // Check that the currentVolume is set - BOOST_CHECK_NE(BVHState.navigation.currentVolume, nullptr); + BOOST_CHECK_NE(bvhState.navigation.currentVolume, nullptr); // Check that the currentVolume is the startVolume - BOOST_CHECK_EQUAL(BVHState.navigation.currentVolume, - BVHState.navigation.startVolume); + BOOST_CHECK_EQUAL(bvhState.navigation.currentVolume, + bvhState.navigation.startVolume); // Check that the currentSurface is reset to: - BOOST_CHECK_EQUAL(BVHState.navigation.currentSurface, nullptr); + BOOST_CHECK_EQUAL(bvhState.navigation.currentSurface, nullptr); // No layer has been found - BOOST_CHECK_EQUAL(BVHState.navigation.navLayers.size(), 0u); + BOOST_CHECK_EQUAL(bvhState.navigation.navLayers.size(), 0u); // ACTORS-ABORTERS-TARGET - navigator.preStep(BVHState, BVHStepper); + bvhNavigator.preStep(bvhState, bvhStepper); // Still no layer - BOOST_CHECK_EQUAL(BVHState.navigation.navLayers.size(), 0u); + BOOST_CHECK_EQUAL(bvhState.navigation.navLayers.size(), 0u); // Surfaces have been found - BOOST_CHECK_EQUAL(BVHState.navigation.navSurfaces.size(), 42u); + BOOST_CHECK_EQUAL(bvhState.navigation.navSurfaces.size(), 42u); } } // namespace Test diff --git a/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp b/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp index 4e17b47b82e..048f29c99b2 100644 --- a/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp +++ b/Tests/UnitTests/Core/TrackFitting/FitterTestsCommon.hpp @@ -12,6 +12,7 @@ #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/MultiTrajectory.hpp" +#include "Acts/EventData/ProxyAccessor.hpp" #include "Acts/EventData/SourceLink.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" #include "Acts/EventData/VectorTrackContainer.hpp" @@ -176,8 +177,8 @@ struct FitterTester { // this is the default option. set anyway for consistency options.referenceSurface = nullptr; - Acts::ConstTrackAccessor reversed{"reversed"}; - Acts::ConstTrackAccessor smoothed{"smoothed"}; + Acts::ConstProxyAccessor reversed{"reversed"}; + Acts::ConstProxyAccessor smoothed{"smoothed"}; auto doTest = [&](bool diag) { Acts::TrackContainer tracks{Acts::VectorTrackContainer{}, @@ -249,8 +250,8 @@ struct FitterTester { BOOST_CHECK(tracks.hasColumn("reversed")); BOOST_CHECK(tracks.hasColumn("smoothed")); - Acts::ConstTrackAccessor reversed{"reversed"}; - Acts::ConstTrackAccessor smoothed{"smoothed"}; + Acts::ConstProxyAccessor reversed{"reversed"}; + Acts::ConstProxyAccessor smoothed{"smoothed"}; // check the output status flags if (doDiag) { @@ -305,8 +306,8 @@ struct FitterTester { BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size()); BOOST_CHECK_EQUAL(track.nHoles(), 0u); - Acts::ConstTrackAccessor reversed{"reversed"}; - Acts::ConstTrackAccessor smoothed{"smoothed"}; + Acts::ConstProxyAccessor reversed{"reversed"}; + Acts::ConstProxyAccessor smoothed{"smoothed"}; // check the output status flags if (doDiag) { BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed); @@ -358,8 +359,8 @@ struct FitterTester { BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size()); BOOST_CHECK_EQUAL(track.nHoles(), 0u); - Acts::ConstTrackAccessor reversed{"reversed"}; - Acts::ConstTrackAccessor smoothed{"smoothed"}; + Acts::ConstProxyAccessor reversed{"reversed"}; + Acts::ConstProxyAccessor smoothed{"smoothed"}; // check the output status flags if (doDiag) { @@ -388,8 +389,8 @@ struct FitterTester { tracks.addColumn("reversed"); tracks.addColumn("smoothed"); - Acts::ConstTrackAccessor reversed{"reversed"}; - Acts::ConstTrackAccessor smoothed{"smoothed"}; + Acts::ConstProxyAccessor reversed{"reversed"}; + Acts::ConstProxyAccessor smoothed{"smoothed"}; // fit w/ all hits in order { @@ -448,8 +449,8 @@ struct FitterTester { tracks.addColumn("reversed"); tracks.addColumn("smoothed"); - Acts::ConstTrackAccessor reversed{"reversed"}; - Acts::ConstTrackAccessor smoothed{"smoothed"}; + Acts::ConstProxyAccessor reversed{"reversed"}; + Acts::ConstProxyAccessor smoothed{"smoothed"}; // always keep the first and last measurement. leaving those in seems to not // count the respective surfaces as holes. @@ -498,8 +499,8 @@ struct FitterTester { tracks.addColumn("reversed"); tracks.addColumn("smoothed"); - Acts::ConstTrackAccessor reversed{"reversed"}; - Acts::ConstTrackAccessor smoothed{"smoothed"}; + Acts::ConstProxyAccessor reversed{"reversed"}; + Acts::ConstProxyAccessor smoothed{"smoothed"}; for (std::size_t i = 0; i < sourceLinks.size(); ++i) { // replace the i-th measurement with an outlier @@ -547,8 +548,8 @@ struct FitterTester { tracks.addColumn("reversed"); tracks.addColumn("smoothed"); - Acts::ConstTrackAccessor reversed{"reversed"}; - Acts::ConstTrackAccessor smoothed{"smoothed"}; + Acts::ConstProxyAccessor reversed{"reversed"}; + Acts::ConstProxyAccessor smoothed{"smoothed"}; auto sourceLinks = prepareSourceLinks(measurements.sourceLinks); diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp index 17fd0376c03..2277197ce1a 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveGridTrackDensityTests.cpp @@ -12,8 +12,6 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/GenericBoundTrackParameters.hpp" #include "Acts/EventData/ParticleHypothesis.hpp" #include "Acts/EventData/TrackParameters.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" @@ -22,12 +20,9 @@ #include "Acts/Utilities/Result.hpp" #include "Acts/Vertexing/AdaptiveGridTrackDensity.hpp" -#include -#include #include #include #include -#include namespace bdata = boost::unit_test::data; using namespace Acts::UnitLiterals; @@ -49,19 +44,17 @@ Covariance makeRandomCovariance(int seed = 31415) { } BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { - using Vector2 = Eigen::Matrix; - using Matrix2 = Eigen::Matrix; // Using a large track grid so we can choose a small bin size - const unsigned int spatialTrkGridSize = 4001; + const std::uint32_t spatialTrkGridSize = 4001; // Arbitrary (but small) bin size - const float binExtent = 3.1e-4; + const double binExtent = 3.1e-4; // Arbitrary impact parameters - const float d0 = 0.4; - const float z0 = -0.2; + const double d0 = 0.4; + const double z0 = -0.2; Vector2 impactParameters{d0, z0}; Covariance covMat = makeRandomCovariance(); - Matrix2 subCovMat = covMat.block<2, 2>(0, 0).cast(); + SquareMatrix2 subCovMat = covMat.block<2, 2>(0, 0); BoundVector paramVec; paramVec << d0, z0, 0, 0, 0, 0; @@ -73,7 +66,9 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { ParticleHypothesis::pion()); AdaptiveGridTrackDensity::Config cfg; - cfg.spatialTrkGridSize = spatialTrkGridSize; + // force track to have exactly spatialTrkGridSize spatial bins for testing + // purposes + cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg.spatialBinExtent = binExtent; AdaptiveGridTrackDensity grid(cfg); @@ -83,66 +78,67 @@ BOOST_AUTO_TEST_CASE(compare_to_analytical_solution_for_single_track) { // Add track auto trackDensityMap = grid.addTrack(params1, mainDensityMap); - float relTol = 1e-5; - float small = 1e-5; + double relTol = 1e-5; + double small = 1e-5; auto gaussian2D = [&](const Vector2& args, const Vector2& mus, - const Matrix2& sigmas) { + const SquareMatrix2& sigmas) { Vector2 diffs = args - mus; - float coef = 1 / std::sqrt(sigmas.determinant()); - float expo = -0.5 * diffs.transpose().dot(sigmas.inverse() * diffs); + double coef = 1 / std::sqrt(sigmas.determinant()); + double expo = -0.5 * diffs.transpose().dot(sigmas.inverse() * diffs); return coef * std::exp(expo); }; - for (auto const& it : mainDensityMap) { + for (const auto& [bin, _] : mainDensityMap) { // Extract variables for better readability - int zBin = it.first.first; - float density = it.second; + std::int32_t zBin = bin.first; + float density = mainDensityMap.at(bin); // Argument for 2D gaussian Vector2 dzVec{0., grid.getBinCenter(zBin, binExtent)}; // Compute correct density... float correctDensity = gaussian2D(dzVec, impactParameters, subCovMat); // ... and check if our result is equivalent - CHECK_CLOSE_OR_SMALL(density, correctDensity, relTol, small); + CHECK_CLOSE_OR_SMALL(correctDensity, density, relTol, small); } // Analytical maximum of the Gaussian (can be obtained by expressing the // exponent as (az - b)^2 + c and noting correctMaxZ = b/a) - float correctMaxZ = + double correctMaxZ = -0.5 * (subCovMat(0, 1) + subCovMat(1, 0)) / subCovMat(0, 0) * d0 + z0; // Analytical FWHM of the Gaussian (result similar to // https://en.wikipedia.org/wiki/Full_width_at_half_maximum#Normal_distribution // but the calculation needs to be slightly modified in our case) - float correctFWHM = 2. * std::sqrt(2 * std::log(2.) * - subCovMat.determinant() / subCovMat(0, 0)); + double correctFWHM = + 2. * + std::sqrt(2 * std::log(2.) * subCovMat.determinant() / subCovMat(0, 0)); // Estimate maximum z position and seed width auto res = grid.getMaxZTPositionAndWidth(mainDensityMap); BOOST_CHECK(res.ok()); // Extract variables for better readability... - float maxZ = res.value().first.first; - float fwhm = res.value().second * 2.355f; + double maxZ = res.value().first.first; + double fwhm = res.value().second * 2.355; // ... and check if they are correct (note: the optimization is not as exact // as the density values). - float relTolOptimization = 1e-3; - CHECK_CLOSE_REL(maxZ, correctMaxZ, relTolOptimization); - CHECK_CLOSE_REL(fwhm, correctFWHM, relTolOptimization); + double relTolOptimization = 1e-3; + CHECK_CLOSE_REL(correctMaxZ, maxZ, relTolOptimization); + CHECK_CLOSE_REL(correctFWHM, fwhm, relTolOptimization); } BOOST_AUTO_TEST_CASE( compare_to_analytical_solution_for_single_track_with_time) { // Number of bins in z- and t-direction - const unsigned int spatialTrkGridSize = 401; - const unsigned int temporalTrkGridSize = 401; + const std::uint32_t spatialTrkGridSize = 401; + const std::uint32_t temporalTrkGridSize = 401; // Bin extents - const float spatialBinExtent = 3.1e-3; - const float temporalBinExtent = 3.1e-3; + const double spatialBinExtent = 3.1e-3; + const double temporalBinExtent = 3.1e-3; // Arbitrary impact parameters - const float d0 = -0.1; - const float z0 = -0.2; - const float t0 = 0.1; + const double d0 = -0.1; + const double z0 = -0.2; + const double t0 = 0.1; Vector3 impactParameters{d0, z0, t0}; // symmetric covariance matrix @@ -160,10 +156,15 @@ BOOST_AUTO_TEST_CASE( ActsSquareMatrix<3> ipCov = params.impactParameterCovariance().value(); AdaptiveGridTrackDensity::Config cfg; - cfg.spatialTrkGridSize = spatialTrkGridSize; + // force track to have exactly spatialTrkGridSize spatial bins for testing + // purposes + cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg.spatialBinExtent = spatialBinExtent; - cfg.temporalTrkGridSize = temporalTrkGridSize; + // force track to have exactly temporalTrkGridSize temporal bins for testing + // purposes + cfg.temporalTrkGridSizeRange = {temporalTrkGridSize, temporalTrkGridSize}; cfg.temporalBinExtent = temporalBinExtent; + cfg.useTime = true; AdaptiveGridTrackDensity grid(cfg); // Empty map @@ -172,22 +173,21 @@ BOOST_AUTO_TEST_CASE( // Add track auto trackDensityMap = grid.addTrack(params, mainDensityMap); - float relTol = 1e-5; - float small = 1e-5; + double relTol = 1e-5; + double small = 1e-5; auto gaussian3D = [&](const Vector3& args, const Vector3& mus, const SquareMatrix3& sigmas) { Vector3 diffs = args - mus; - float coef = 1 / std::sqrt(sigmas.determinant()); - float expo = -0.5 * diffs.transpose().dot(sigmas.inverse() * diffs); + double coef = 1 / std::sqrt(sigmas.determinant()); + double expo = -0.5 * diffs.transpose().dot(sigmas.inverse() * diffs); return coef * std::exp(expo); }; - for (auto const& it : mainDensityMap) { + for (const auto& [bin, density] : mainDensityMap) { // Extract variables for better readability - float z = grid.getBinCenter(it.first.first, spatialBinExtent); - float t = grid.getBinCenter(it.first.second, temporalBinExtent); - float density = it.second; + double z = grid.getBinCenter(bin.first, spatialBinExtent); + double t = grid.getBinCenter(bin.second, temporalBinExtent); // Argument for 3D gaussian Vector3 dztVec{0., z, t}; @@ -195,7 +195,7 @@ BOOST_AUTO_TEST_CASE( float correctDensity = gaussian3D(dztVec, impactParameters, ipCov); // ... and check if our result is equivalent - CHECK_CLOSE_OR_SMALL(density, correctDensity, relTol, small); + CHECK_CLOSE_OR_SMALL(correctDensity, density, relTol, small); } // The analytical calculations of the following can be found here: @@ -221,24 +221,24 @@ BOOST_AUTO_TEST_CASE( BOOST_CHECK(res.ok()); // Extract variables for better readability... - float maxZ = res.value().first.first; - float maxT = res.value().first.second; - float fwhm = res.value().second * 2.355f; + double maxZ = res.value().first.first; + double maxT = res.value().first.second; + double fwhm = res.value().second * 2.355; // ... and check if they are correct (note: the optimization is not as exact // as the density values). - float relTolOptimization = 1e-1; - CHECK_CLOSE_REL(maxZ, correctMaxZ, relTolOptimization); - CHECK_CLOSE_REL(maxT, correctMaxT, relTolOptimization); - CHECK_CLOSE_REL(fwhm, correctFWHM, relTolOptimization); + double relTolOptimization = 1e-1; + CHECK_CLOSE_REL(correctMaxZ, maxZ, relTolOptimization); + CHECK_CLOSE_REL(correctMaxT, maxT, relTolOptimization); + CHECK_CLOSE_REL(correctFWHM, fwhm, relTolOptimization); } BOOST_AUTO_TEST_CASE(seed_width_estimation) { // Dummy track grid size (not needed for this unit test) - const unsigned int spatialTrkGridSize = 1; - float binExtent = 2.; + const std::uint32_t spatialTrkGridSize = 1; + double binExtent = 2.; AdaptiveGridTrackDensity::Config cfg; - cfg.spatialTrkGridSize = spatialTrkGridSize; + cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg.spatialBinExtent = binExtent; AdaptiveGridTrackDensity grid(cfg); @@ -246,14 +246,14 @@ BOOST_AUTO_TEST_CASE(seed_width_estimation) { AdaptiveGridTrackDensity::DensityMap mainDensityMap; // z-position of the maximum density - float correctMaxZ = -2.; + double correctMaxZ = -2.; // Fill map with a triangular track density. // We use an isoscele triangle with a maximum density value of 1 and a width // of 20 mm. The linear approximation we use during the seed width estimation // should be exact in this case. - for (int i = -6; i <= 4; i++) { - mainDensityMap[std::make_pair(i, 0)] = + for (std::int32_t i = -6; i <= 4; i++) { + mainDensityMap[{i, 0}] = 1.0 - 0.1 * std::abs(correctMaxZ - grid.getBinCenter(i, binExtent)); } @@ -262,22 +262,22 @@ BOOST_AUTO_TEST_CASE(seed_width_estimation) { BOOST_CHECK(res.ok()); // Check if we found the correct maximum - float maxZ = res.value().first.first; - BOOST_CHECK_EQUAL(maxZ, correctMaxZ); + double maxZ = res.value().first.first; + BOOST_CHECK_EQUAL(correctMaxZ, maxZ); // Calculate full width at half maximum from seed width and check if it's // correct - float fwhm = res.value().second * 2.355f; - BOOST_CHECK_EQUAL(fwhm, 10.); + double fwhm = res.value().second * 2.355; + CHECK_CLOSE_REL(10., fwhm, 1e-5); } BOOST_AUTO_TEST_CASE(track_adding) { - const unsigned int spatialTrkGridSize = 15; + const std::uint32_t spatialTrkGridSize = 15; double binExtent = 0.1; // mm AdaptiveGridTrackDensity::Config cfg; - cfg.spatialTrkGridSize = spatialTrkGridSize; + cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg.spatialBinExtent = binExtent; AdaptiveGridTrackDensity grid(cfg); @@ -317,49 +317,50 @@ BOOST_AUTO_TEST_CASE(track_adding) { // Track should have been entirely added to both grids trackDensityMap = grid.addTrack(params1, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), spatialTrkGridSize); + BOOST_CHECK_EQUAL(spatialTrkGridSize, mainDensityMap.size()); // Track should have been entirely added to both grids trackDensityMap = grid.addTrack(params2, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 2 * spatialTrkGridSize); + BOOST_CHECK_EQUAL(2 * spatialTrkGridSize, mainDensityMap.size()); // Track 3 has overlap of 2 bins with track 1 trackDensityMap = grid.addTrack(params3, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * spatialTrkGridSize - 2); + BOOST_CHECK_EQUAL(3 * spatialTrkGridSize - 2, mainDensityMap.size()); // Add first track again, should *not* introduce new z entries trackDensityMap = grid.addTrack(params1, mainDensityMap); - BOOST_CHECK_EQUAL(mainDensityMap.size(), 3 * spatialTrkGridSize - 2); + BOOST_CHECK_EQUAL(3 * spatialTrkGridSize - 2, mainDensityMap.size()); } BOOST_AUTO_TEST_CASE(max_z_t_and_width) { - const unsigned int spatialTrkGridSize = 29; - const unsigned int temporalTrkGridSize = 29; + const std::uint32_t spatialTrkGridSize = 29; + const std::uint32_t temporalTrkGridSize = 29; // spatial and temporal bin extent double binExtent = 0.05; // 1D grid of z values AdaptiveGridTrackDensity::Config cfg1D; - cfg1D.spatialTrkGridSize = spatialTrkGridSize; + cfg1D.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg1D.spatialBinExtent = binExtent; AdaptiveGridTrackDensity grid1D(cfg1D); // 2D grid of z and t values AdaptiveGridTrackDensity::Config cfg2D; - cfg2D.spatialTrkGridSize = spatialTrkGridSize; + cfg2D.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg2D.spatialBinExtent = binExtent; - cfg2D.temporalTrkGridSize = temporalTrkGridSize; + cfg2D.temporalTrkGridSizeRange = {temporalTrkGridSize, temporalTrkGridSize}; cfg2D.temporalBinExtent = binExtent; + cfg2D.useTime = true; AdaptiveGridTrackDensity grid2D(cfg2D); // Create some test tracks Covariance covMat(Covariance::Identity() * 0.005); - float z0Trk1 = 0.25; - float t0Trk1 = 0.05; - float z0Trk2 = -10.95; - float t0Trk2 = 0.1; + double z0Trk1 = 0.25; + double t0Trk1 = 0.05; + double z0Trk2 = -10.95; + double t0Trk2 = 0.1; BoundVector paramVec1; paramVec1 << 0.02, z0Trk1, 0, 0, 0, t0Trk1; BoundVector paramVec2; @@ -383,18 +384,18 @@ BOOST_AUTO_TEST_CASE(max_z_t_and_width) { auto firstRes = grid1D.getMaxZTPosition(mainDensityMap1D); BOOST_CHECK(firstRes.ok()); // Maximum should be at z0Trk1 position ... - BOOST_CHECK_EQUAL((*firstRes).first, z0Trk1); + BOOST_CHECK_EQUAL(z0Trk1, (*firstRes).first); // ... and the corresponding time should be set to 0 - BOOST_CHECK_EQUAL((*firstRes).second, 0.); + BOOST_CHECK_EQUAL(0., (*firstRes).second); // Add first track to 2D grid auto trackDensityMap2D = grid2D.addTrack(params1, mainDensityMap2D); auto firstRes2D = grid2D.getMaxZTPosition(mainDensityMap2D); BOOST_CHECK(firstRes2D.ok()); // Maximum should be at z0Trk1 position ... - BOOST_CHECK_EQUAL((*firstRes2D).first, z0Trk1); + BOOST_CHECK_EQUAL(z0Trk1, (*firstRes2D).first); // ... and the corresponding time should be at t0Trk1 - BOOST_CHECK_EQUAL((*firstRes2D).second, t0Trk1); + BOOST_CHECK_EQUAL(t0Trk1, (*firstRes2D).second); // Add second track to spatial grid trackDensityMap = grid1D.addTrack(params2, mainDensityMap1D); @@ -403,11 +404,11 @@ BOOST_AUTO_TEST_CASE(max_z_t_and_width) { BOOST_CHECK(secondRes.ok()); // Trk 2 is closer to z-axis and should thus yield higher density values. // Therefore, the new maximum is at z0Trk2 ... - BOOST_CHECK_EQUAL((*secondRes).first.first, z0Trk2); + CHECK_CLOSE_OR_SMALL(z0Trk2, (*secondRes).first.first, 1e-5, 1e-5); // ... the corresponding time should be set to 0... - BOOST_CHECK_EQUAL((*secondRes).first.second, 0.); + BOOST_CHECK_EQUAL(0., (*secondRes).first.second); // ... and it should have a positive width - BOOST_CHECK_GT((*secondRes).second, 0); + BOOST_CHECK_LT(0, (*secondRes).second); // Add second track to 2D grid trackDensityMap = grid2D.addTrack(params2, mainDensityMap2D); @@ -416,20 +417,20 @@ BOOST_AUTO_TEST_CASE(max_z_t_and_width) { BOOST_CHECK(secondRes2D.ok()); // Trk 2 is closer to z-axis and should thus yield higher density values. // Therefore, the new maximum is at z0Trk2 ... - BOOST_CHECK_EQUAL((*secondRes2D).first.first, z0Trk2); + CHECK_CLOSE_OR_SMALL(z0Trk2, (*secondRes2D).first.first, 1e-5, 1e-5); // ... the corresponding time should be at t0Trk2 ... - BOOST_CHECK_EQUAL((*secondRes2D).first.second, t0Trk2); + BOOST_CHECK_EQUAL(t0Trk2, (*secondRes2D).first.second); // ... and it should have approximately the same width in z direction CHECK_CLOSE_OR_SMALL((*secondRes2D).second, (*secondRes).second, 1e-5, 1e-5); } BOOST_AUTO_TEST_CASE(highest_density_sum) { - const unsigned int spatialTrkGridSize = 29; + const std::uint32_t spatialTrkGridSize = 29; double binExtent = 0.05; // mm AdaptiveGridTrackDensity::Config cfg; - cfg.spatialTrkGridSize = spatialTrkGridSize; + cfg.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg.spatialBinExtent = binExtent; cfg.useHighestSumZPosition = true; @@ -438,8 +439,8 @@ BOOST_AUTO_TEST_CASE(highest_density_sum) { // Create some test tracks Covariance covMat(Covariance::Identity() * 0.005); - float z0Trk1 = 0.25; - float z0Trk2 = -10.95; + double z0Trk1 = 0.25; + double z0Trk2 = -10.95; BoundVector paramVec1; paramVec1 << 0.01, z0Trk1, 0, 0, 0, 0; BoundVector paramVec2; @@ -463,7 +464,7 @@ BOOST_AUTO_TEST_CASE(highest_density_sum) { auto res1 = grid.getMaxZTPosition(mainDensityMap); BOOST_CHECK(res1.ok()); // Maximum should be at z0Trk1 position - BOOST_CHECK_EQUAL((*res1).first, z0Trk1); + BOOST_CHECK_EQUAL(z0Trk1, (*res1).first); // Add second track trackDensityMap = grid.addTrack(params2, mainDensityMap); @@ -471,52 +472,52 @@ BOOST_AUTO_TEST_CASE(highest_density_sum) { BOOST_CHECK(res2.ok()); // Trk 2 is closer to z-axis and should yield higher density values // New maximum is therefore at z0Trk2 - BOOST_CHECK_EQUAL((*res2).first, z0Trk2); + CHECK_CLOSE_REL(z0Trk2, (*res2).first, 1e-5); // Add small density values around the maximum of track 1 - const float densityToAdd = 0.5; - mainDensityMap.at(std::make_pair(4, 0)) += densityToAdd; - mainDensityMap.at(std::make_pair(6, 0)) += densityToAdd; + mainDensityMap[{4, 0}] += 0.5; + mainDensityMap[{6, 0}] += 0.5; auto res3 = grid.getMaxZTPosition(mainDensityMap); BOOST_CHECK(res3.ok()); // Trk 2 still has the highest peak density value, however, the small // added densities for track 1 around its maximum should now lead to // a prediction at z0Trk1 again - BOOST_CHECK_EQUAL((*res3).first, z0Trk1); + BOOST_CHECK_EQUAL(z0Trk1, (*res3).first); } BOOST_AUTO_TEST_CASE(track_removing) { - const unsigned int spatialTrkGridSize = 29; - const unsigned int temporalTrkGridSize = 29; + const std::uint32_t spatialTrkGridSize = 29; + const std::uint32_t temporalTrkGridSize = 29; - int trkGridSize = spatialTrkGridSize * temporalTrkGridSize; + std::uint32_t trkGridSize = spatialTrkGridSize * temporalTrkGridSize; // bin extent in z and t direction double binExtent = 0.05; // 1D grid AdaptiveGridTrackDensity::Config cfg1D; - cfg1D.spatialTrkGridSize = spatialTrkGridSize; + cfg1D.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg1D.spatialBinExtent = binExtent; AdaptiveGridTrackDensity grid1D(cfg1D); // 2D grid AdaptiveGridTrackDensity::Config cfg2D; - cfg2D.spatialTrkGridSize = spatialTrkGridSize; + cfg2D.spatialTrkGridSizeRange = {spatialTrkGridSize, spatialTrkGridSize}; cfg2D.spatialBinExtent = binExtent; - cfg2D.temporalTrkGridSize = temporalTrkGridSize; + cfg2D.temporalTrkGridSizeRange = {temporalTrkGridSize, temporalTrkGridSize}; cfg2D.temporalBinExtent = binExtent; + cfg2D.useTime = true; AdaptiveGridTrackDensity grid2D(cfg2D); // Create some test tracks Covariance covMat = makeRandomCovariance(); // Define z0 values for test tracks - float z0Trk1 = -0.45; - float t0Trk1 = -0.15; - float z0Trk2 = -0.25; - float t0Trk2 = t0Trk1; + double z0Trk1 = -0.45; + double t0Trk1 = -0.15; + double z0Trk2 = -0.25; + double t0Trk2 = t0Trk1; BoundVector paramVec0; paramVec0 << 0.1, z0Trk1, 0, 0, 0, t0Trk1; @@ -538,9 +539,9 @@ BOOST_AUTO_TEST_CASE(track_removing) { // Lambda for calculating total density auto densitySum = [](const auto& densityMap) { - float sum = 0.; - for (auto it = densityMap.begin(); it != densityMap.end(); it++) { - sum += it->second; + double sum = 0.; + for (const auto& [_, density] : densityMap) { + sum += density; } return sum; }; @@ -549,23 +550,23 @@ BOOST_AUTO_TEST_CASE(track_removing) { auto firstTrackDensityMap1D = grid1D.addTrack(params0, mainDensityMap1D); BOOST_CHECK(!mainDensityMap1D.empty()); // Grid size should match spatialTrkGridSize - BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); - float firstDensitySum1D = densitySum(mainDensityMap1D); + BOOST_CHECK_EQUAL(spatialTrkGridSize, mainDensityMap1D.size()); + double firstDensitySum1D = densitySum(mainDensityMap1D); // Add track 0 to 2D grid auto firstTrackDensityMap2D = grid2D.addTrack(params0, mainDensityMap2D); BOOST_CHECK(!mainDensityMap2D.empty()); // Grid size should match spatialTrkGridSize - BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); - float firstDensitySum2D = densitySum(mainDensityMap2D); + BOOST_CHECK_EQUAL(trkGridSize, mainDensityMap2D.size()); + double firstDensitySum2D = densitySum(mainDensityMap2D); // Add track 0 again to 1D grid firstTrackDensityMap1D = grid1D.addTrack(params0, mainDensityMap1D); BOOST_CHECK(!mainDensityMap1D.empty()); // Grid size should still match spatialTrkGridSize - BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); + BOOST_CHECK_EQUAL(spatialTrkGridSize, mainDensityMap1D.size()); // Calculate new total density ... - float secondDensitySum1D = densitySum(mainDensityMap1D); + double secondDensitySum1D = densitySum(mainDensityMap1D); // ... and check that it's twice as large as before BOOST_CHECK_EQUAL(2 * firstDensitySum1D, secondDensitySum1D); @@ -573,56 +574,58 @@ BOOST_AUTO_TEST_CASE(track_removing) { firstTrackDensityMap2D = grid2D.addTrack(params0, mainDensityMap2D); BOOST_CHECK(!mainDensityMap2D.empty()); // Grid size should still match trkGridSize - BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); + BOOST_CHECK_EQUAL(trkGridSize, mainDensityMap2D.size()); // Calculate new total density ... - float secondDensitySum2D = densitySum(mainDensityMap2D); + double secondDensitySum2D = densitySum(mainDensityMap2D); // ... and check that it's twice as large as before BOOST_CHECK_EQUAL(2 * firstDensitySum2D, secondDensitySum2D); // Remove track 0 from 1D grid grid1D.subtractTrack(firstTrackDensityMap1D, mainDensityMap1D); // Calculate new total density - float thirdDensitySum1D = densitySum(mainDensityMap1D); + double thirdDensitySum1D = densitySum(mainDensityMap1D); // Density should be old one again BOOST_CHECK_EQUAL(firstDensitySum1D, thirdDensitySum1D); // Grid size should still match spatialTrkGridSize (removal does not change // grid size) - BOOST_CHECK_EQUAL(mainDensityMap1D.size(), spatialTrkGridSize); + BOOST_CHECK_EQUAL(spatialTrkGridSize, mainDensityMap1D.size()); // Remove track 0 from 2D grid grid2D.subtractTrack(firstTrackDensityMap2D, mainDensityMap2D); // Calculate new total density - float thirdDensitySum2D = densitySum(mainDensityMap2D); + double thirdDensitySum2D = densitySum(mainDensityMap2D); // Density should be old one again BOOST_CHECK_EQUAL(firstDensitySum2D, thirdDensitySum2D); // Grid size should still match trkGridSize (removal does not change grid // size) - BOOST_CHECK_EQUAL(mainDensityMap2D.size(), trkGridSize); + BOOST_CHECK_EQUAL(trkGridSize, mainDensityMap2D.size()); // Add track 1 to 1D grid (overlaps with track 0!) auto secondTrackDensityMap1D = grid1D.addTrack(params1, mainDensityMap1D); - int nNonOverlappingBins1D = int(std::abs(z0Trk1 - z0Trk2) / binExtent + 1); - BOOST_CHECK_EQUAL(mainDensityMap1D.size(), - spatialTrkGridSize + nNonOverlappingBins1D); - float fourthDensitySum1D = densitySum(mainDensityMap1D); + std::uint32_t nNonOverlappingBins1D = + (std::uint32_t)(std::abs(z0Trk1 - z0Trk2) / binExtent); + BOOST_CHECK_EQUAL(spatialTrkGridSize + nNonOverlappingBins1D, + mainDensityMap1D.size()); + double fourthDensitySum1D = densitySum(mainDensityMap1D); // Add track 1 to 2D grid (overlaps with track 0!) auto secondTrackDensityMap2D = grid2D.addTrack(params1, mainDensityMap2D); - int nNonOverlappingBins2D = nNonOverlappingBins1D * temporalTrkGridSize; - BOOST_CHECK_EQUAL(mainDensityMap2D.size(), - trkGridSize + nNonOverlappingBins2D); - float fourthDensitySum2D = densitySum(mainDensityMap2D); + std::uint32_t nNonOverlappingBins2D = + nNonOverlappingBins1D * temporalTrkGridSize; + BOOST_CHECK_EQUAL(trkGridSize + nNonOverlappingBins2D, + mainDensityMap2D.size()); + double fourthDensitySum2D = densitySum(mainDensityMap2D); // Remove second track 0 from 1D grid grid1D.subtractTrack(firstTrackDensityMap1D, mainDensityMap1D); - float fifthDensitySum1D = densitySum(mainDensityMap1D); + double fifthDensitySum1D = densitySum(mainDensityMap1D); // Density should match differences of removed tracks CHECK_CLOSE_REL(fifthDensitySum1D, fourthDensitySum1D - firstDensitySum1D, 1e-5); // Remove second track 0 from 2D grid grid2D.subtractTrack(firstTrackDensityMap2D, mainDensityMap2D); - float fifthDensitySum2D = densitySum(mainDensityMap2D); + double fifthDensitySum2D = densitySum(mainDensityMap2D); // Density should match differences of removed tracks CHECK_CLOSE_REL(fifthDensitySum2D, fourthDensitySum2D - firstDensitySum2D, 1e-5); @@ -630,20 +633,20 @@ BOOST_AUTO_TEST_CASE(track_removing) { // Remove track 1 from 1D grid grid1D.subtractTrack(secondTrackDensityMap1D, mainDensityMap1D); // Size should not have changed - BOOST_CHECK_EQUAL(mainDensityMap1D.size(), - spatialTrkGridSize + nNonOverlappingBins1D); - float sixthDensitySum1D = densitySum(mainDensityMap1D); + BOOST_CHECK_EQUAL(spatialTrkGridSize + nNonOverlappingBins1D, + mainDensityMap1D.size()); + double sixthDensitySum1D = densitySum(mainDensityMap1D); // 1D grid is now empty after all tracks were removed - CHECK_CLOSE_ABS(sixthDensitySum1D, 0., 1e-4); + CHECK_CLOSE_ABS(0., sixthDensitySum1D, 1e-4); // Remove track 1 from 2D grid grid2D.subtractTrack(secondTrackDensityMap2D, mainDensityMap2D); // Size should not have changed - BOOST_CHECK_EQUAL(mainDensityMap2D.size(), - trkGridSize + nNonOverlappingBins2D); - float sixthDensitySum2D = densitySum(mainDensityMap2D); + BOOST_CHECK_EQUAL(trkGridSize + nNonOverlappingBins2D, + mainDensityMap2D.size()); + double sixthDensitySum2D = densitySum(mainDensityMap2D); // 2D grid is now empty after all tracks were removed - CHECK_CLOSE_ABS(sixthDensitySum2D, 0., 1e-4); + CHECK_CLOSE_ABS(0., sixthDensitySum2D, 1e-4); } } // namespace Test diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp index 51e24623943..4d823a6fc1a 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp @@ -560,7 +560,9 @@ BOOST_AUTO_TEST_CASE( // Grid density used during vertex seed finding AdaptiveGridTrackDensity::Config gridDensityCfg; - gridDensityCfg.spatialTrkGridSize = 55; + // force track to have exactly spatialTrkGridSize spatial bins for testing + // purposes + gridDensityCfg.spatialTrkGridSizeRange = {55, 55}; gridDensityCfg.spatialBinExtent = 0.05; AdaptiveGridTrackDensity gridDensity(gridDensityCfg); diff --git a/Tests/UnitTests/Core/Vertexing/CMakeLists.txt b/Tests/UnitTests/Core/Vertexing/CMakeLists.txt index 2cf6309574f..c0a7ac4fa0c 100644 --- a/Tests/UnitTests/Core/Vertexing/CMakeLists.txt +++ b/Tests/UnitTests/Core/Vertexing/CMakeLists.txt @@ -11,4 +11,4 @@ add_unittest(TrackDensityVertexFinder TrackDensityVertexFinderTests.cpp) add_unittest(GaussianGridTrackDensity GaussianGridTrackDensityTests.cpp) add_unittest(GridDensityVertexFinder GridDensityVertexFinderTests.cpp) add_unittest(AdaptiveGridTrackDensity AdaptiveGridTrackDensityTests.cpp) -add_unittest(SingleSeedVertexFinder SingleSeedVertexFinderTests.cpp) \ No newline at end of file +add_unittest(SingleSeedVertexFinder SingleSeedVertexFinderTests.cpp) diff --git a/Tests/UnitTests/Core/Vertexing/GridDensityVertexFinderTests.cpp b/Tests/UnitTests/Core/Vertexing/GridDensityVertexFinderTests.cpp index a73169c7b80..55f092032d1 100644 --- a/Tests/UnitTests/Core/Vertexing/GridDensityVertexFinderTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/GridDensityVertexFinderTests.cpp @@ -14,7 +14,6 @@ #include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" -#include "Acts/EventData/Charge.hpp" #include "Acts/EventData/GenericBoundTrackParameters.hpp" #include "Acts/EventData/TrackParameters.hpp" #include "Acts/Geometry/GeometryContext.hpp" @@ -32,12 +31,9 @@ #include "Acts/Vertexing/Vertex.hpp" #include "Acts/Vertexing/VertexingOptions.hpp" -#include -#include #include #include #include -#include #include #include @@ -76,7 +72,7 @@ std::uniform_real_distribution etaDist(-4., 4.); /// of track density values /// BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_test) { - bool debugMode = true; + bool debugMode = false; // Note that the AdaptiveGridTrackDensity and the GaussianGridTrackDensity // only furnish exactly the same results for uneven mainGridSize, where the @@ -109,7 +105,7 @@ BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_test) { // Use custom grid density here with same bin size as Finder1 AdaptiveGridTrackDensity::Config adaptiveDensityConfig; - adaptiveDensityConfig.spatialTrkGridSize = trkGridSize; + adaptiveDensityConfig.spatialTrkGridSizeRange = {trkGridSize, trkGridSize}; adaptiveDensityConfig.spatialBinExtent = 2. / 30.01 * 1_mm; AdaptiveGridTrackDensity adaptiveDensity(adaptiveDensityConfig); @@ -174,7 +170,8 @@ BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_test) { BOOST_CHECK(!(*res1).empty()); Vector3 result1 = (*res1).back().position(); if (debugMode) { - std::cout << "Vertex position result 1: " << result1 << std::endl; + std::cout << "Vertex position result 1: " << result1.transpose() + << std::endl; } CHECK_CLOSE_ABS(result1[eZ], zVertexPos1, 1_mm); zResult1 = result1[eZ]; @@ -185,18 +182,19 @@ BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_test) { BOOST_CHECK(!(*res2).empty()); Vector3 result2 = (*res2).back().position(); if (debugMode) { - std::cout << "Vertex position result 2: " << result2 << std::endl; + std::cout << "Vertex position result 2: " << result2.transpose() + << std::endl; } CHECK_CLOSE_ABS(result2[eZ], zVertexPos1, 1_mm); zResult2 = result2[eZ]; } // Both finders should give same results - BOOST_CHECK_EQUAL(zResult1, zResult2); + CHECK_CLOSE_REL(zResult1, zResult2, 1e-5); } BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_track_caching_test) { - bool debugMode = true; + bool debugMode = false; const int mainGridSize = 3001; const int trkGridSize = 35; @@ -225,7 +223,7 @@ BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_track_caching_test) { // Use custom grid density here with same bin size as Finder1 AdaptiveGridTrackDensity::Config adaptiveDensityConfig; - adaptiveDensityConfig.spatialTrkGridSize = trkGridSize; + adaptiveDensityConfig.spatialTrkGridSizeRange = {trkGridSize, trkGridSize}; adaptiveDensityConfig.spatialBinExtent = 2. / 30.01 * 1_mm; adaptiveDensityConfig.useHighestSumZPosition = true; AdaptiveGridTrackDensity adaptiveDensity(adaptiveDensityConfig); @@ -289,7 +287,7 @@ BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_track_caching_test) { BOOST_CHECK(!(*res1).empty()); Vector3 result = (*res1).back().position(); if (debugMode) { - std::cout << "Vertex position after first fill 1: " << result + std::cout << "Vertex position after first fill 1: " << result.transpose() << std::endl; } CHECK_CLOSE_ABS(result[eZ], zVertexPos1, 1_mm); @@ -304,14 +302,14 @@ BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_track_caching_test) { BOOST_CHECK(!(*res2).empty()); Vector3 result = (*res2).back().position(); if (debugMode) { - std::cout << "Vertex position after first fill 2: " << result + std::cout << "Vertex position after first fill 2: " << result.transpose() << std::endl; } CHECK_CLOSE_ABS(result[eZ], zVertexPos1, 1_mm); zResult2 = result[eZ]; } - BOOST_CHECK_EQUAL(zResult1, zResult2); + CHECK_CLOSE_REL(zResult1, zResult2, 1e-5); int trkCount = 0; std::vector removedTracks; @@ -357,14 +355,14 @@ BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_track_caching_test) { zResult2 = result[eZ]; } - BOOST_CHECK_EQUAL(zResult1, zResult2); + CHECK_CLOSE_REL(zResult1, zResult2, 1e-5); } /// /// @brief Unit test for GridDensityVertexFinder with seed with estimation /// BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_seed_width_test) { - bool debugMode = true; + bool debugMode = false; const int mainGridSize = 3001; const int trkGridSize = 35; @@ -391,7 +389,7 @@ BOOST_AUTO_TEST_CASE(grid_density_vertex_finder_seed_width_test) { // Use custom grid density here with same bin size as Finder1 AdaptiveGridTrackDensity::Config adaptiveDensityConfig; - adaptiveDensityConfig.spatialTrkGridSize = trkGridSize; + adaptiveDensityConfig.spatialTrkGridSizeRange = {trkGridSize, trkGridSize}; adaptiveDensityConfig.spatialBinExtent = 2. / 30.01 * 1_mm; AdaptiveGridTrackDensity adaptiveDensity(adaptiveDensityConfig); diff --git a/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp b/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp index 4020d767a5e..cd983d348e2 100644 --- a/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/ImpactPointEstimatorTests.cpp @@ -26,7 +26,7 @@ #include "Acts/Propagator/EigenStepper.hpp" #include "Acts/Propagator/Propagator.hpp" #include "Acts/Propagator/StraightLineStepper.hpp" -#include "Acts/Propagator/detail/VoidPropagatorComponents.hpp" +#include "Acts/Propagator/VoidNavigator.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" @@ -98,7 +98,7 @@ Estimator makeEstimator(double bZ) { Stepper stepper(field); Estimator::Config cfg(field, std::make_shared( - std::move(stepper), detail::VoidNavigator(), + std::move(stepper), VoidNavigator(), getDefaultLogger("Prop", Logging::Level::WARNING))); return Estimator(cfg); } diff --git a/docs/core/geometry/layerless/layerless.md b/docs/core/geometry/layerless/layerless.md index 0debd65e2a7..7f18561640e 100644 --- a/docs/core/geometry/layerless/layerless.md +++ b/docs/core/geometry/layerless/layerless.md @@ -111,7 +111,7 @@ The detector object is the holder class of all geometry objects, it has to conta - a volume finder delegate (as `Acts::Experimental::DetectorVolumeFinder`) that allows to uniquely associate a point in space with a contained volume of the detector. :::{note} -When the detector is constructed, name duplicates are checked for and if found a `std::exception` is thrown. Similarly, when sensitive surfaces are provided and duplicate `Acts::GeometryIdentifier` objects are found during detector construction a `std::exception` is thrown. The latter can be avoided by using an appropriate (set of) `Acts::GeometyIdGenerator` tool(s) which will guarantee some level of uniqueness. +When the detector is constructed, name duplicates are checked for and if found a `std::exception` is thrown. Similarly, when sensitive surfaces are provided and duplicate `Acts::GeometryIdentifier` objects are found during detector construction a `std::exception` is thrown. The latter can be avoided by using an appropriate (set of) `Acts::geometryIdGenerator` tool(s) which will guarantee some level of uniqueness. ::: :::{figure} ../figures/ODD_Detector.png diff --git a/docs/core/reconstruction/track_fitting.md b/docs/core/reconstruction/track_fitting.md index 659ae7d87f1..cf617e798da 100644 --- a/docs/core/reconstruction/track_fitting.md +++ b/docs/core/reconstruction/track_fitting.md @@ -6,9 +6,9 @@ We can run the track fitting algorithms, after we allocated all hits to single t It is not necessary, that all points of a track are present. Currently, we have implementations for three different fitters: -* Kalman Filter -* GSF -* Global Chi-Square Fitter (GX2F) [wip] +* Kalman Filter (KF) +* Gaussian Sum Filter (GSF) +* Global Chi-Square Fitter (GX2F) [in development] Even though all of them are least-squares fits, the concepts are quite different. Therefore, we should not expect identical results from all of them. @@ -126,10 +126,255 @@ The {class}`Acts::AtlasBetheHeitlerApprox` is constructed with two parameterizat * *R Frühwirth*, A Gaussian-mixture approximation of the Bethe–Heitler model of electron energy loss by bremsstrahlung, 2003, see [here](https://doi.org/10.1016/S0010-4655(03)00292-3) (gx2f_core)= -## Global Chi-Square Fitter (GX2F) [wip] +## Global Chi-Square Fitter (GX2F) + +In general the *GX2F* is a weighted least squares fit, minimising the + +$$ +\chi^2 = \sum_i \frac{r_i^2}{\sigma_i^2} +$$ + +of a track. +Here, $r_i$ are our residuals that we weight with $\sigma_i^2$, the covariance of the measurement (a detector property). +Unlike the *KF* and the *GSF*, the *GX2F* looks at all measurements at the same time and iteratively minimises the starting parameters. + +With the *GX2F* we can obtain the final parameters $\vec\alpha_n$ from starting parameters $\vec\alpha_0$. +We set the $\chi^2 = \chi^2(\vec\alpha)$ as a function of the track parameters, but the $\chi^2$-minimisation could be used for many other problems. +Even in the context of track fitting, we are quite free on how to use the *GX2F*. +Especially the residuals $r_i$ can have many interpretations. +Most of the time we will see them as the distance between a measurement and our prediction. +But we can also use scattering angles, energy loss, ... as residuals. +Therefore, the subscript $i$ stands most of the time for a measurement surface, since we want to go over all of them. + +This chapter on the *GX2F* guides through: +- Mathematical description of the base algorithm +- Mathematical description of the multiple scattering +- (coming soon) Mathematical description of the energy loss +- Implementation in ACTS +- Pros/Cons + +### Mathematical description of the base algorithm + +:::{note} +The mathematical derivation is shortened at some places. +There will be a publication including the full derivation coming soon. +::: + +To begin with, there will be a short overview on the algorithm. +Later in this section, each step is described in more detail. +1. Minimise the $\chi^2$ function +2. Update the initial parameters (iteratively) +3. Calculate the covariance for the final parameters + +But before going into detail, we need to introduce a few symbols. +As already mentioned, we have our track parameters $\vec\alpha$ that we want to fit. +To fit them we, we need to calculate our residuals as + +$$ +r_i = m_i - f_i^m(\vec\alpha) +$$ + +where $f^m(\vec\alpha)$ is the projection of our propagation function $f(\vec\alpha)$ into the measurement dimension. +Basically, if we have a pixel measurement we project onto the surface, discarding all angular information. +This projection could be different for each measurement surface. + +#### 1. Minimise the $\chi^2$ function + +We expect the minimum of the $\chi^2$ function at + +$$ +\frac{\partial\chi^2(\vec\alpha)}{\partial\vec\alpha} = 0. +$$ + +To find the zero(s) of this function we could use any method, but we will stick to a modified [Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method), +since it requires just another derivative of the $\chi^2$ function. + +#### 2. Update the initial parameters (iteratively) + +Since we are using the Newton-Raphson method to find the minimum of the $\chi^2$ function, we need to iterate. +Each iteration (should) give as improved parameters $\vec\alpha$. +While iterating we update a system, therefore we want to bring it in this form: + +$$ +\vec\alpha_{n+i} = \vec\alpha_n + \vec{\delta\alpha}_n. +$$ + +After some derivations of the $\chi^2$ function and the Newton-Raphson method, we find matrix equation to calculate $\vec{\delta\alpha}_n$: + +$$ +[a_{kl}] \vec{\delta\alpha}_n = \vec b +$$ + +with + +$$ +a_{kl} = \sum_{i=1}^N \frac{1}{\sigma_i^2} \frac{\partial f_i^m(\vec\alpha)}{\partial \alpha_k}\frac{\partial f_i^m(\vec\alpha)}{\partial \alpha_l}\\ +$$ + +(where we omitted second order derivatives) and + +$$ +b_k = \sum_{i=1}^N \frac{r_i}{\sigma_i^2} \frac{\partial f_i^m(\vec\alpha)}{\partial \alpha_k}. +$$ + +At first sight, these expression might seem intimidating and hard to compute. +But having a closer look, we see, that those derivatives already exist in our framework. +All derivatives are elements of the Jacobian + +$$ +\mathbf{J} = \begin{pmatrix} + \cdot & \dots & \cdot\\ + \vdots & \frac{\partial f^m(\vec\alpha)}{\partial \alpha_k} & \vdots\\ + \cdot & \dots & \cdot + \end{pmatrix}. +$$ + +At this point we got all information to perform a parameter update and repeat until the parameters $\vec\alpha$ converge. + +#### 3. Calculate the covariance for the final parameters + +The calculation of the covariance of the final parameters is quite simple compared to the steps before: + +$$ +cov_{\vec\alpha} = [a_{kl}]^{-1} +$$ + +Since it only depends on the $[a_{kl}]$ of the last iteration, the *GX2F* does not need an initial estimate for the covariance. + +### Mathematical description of the multiple scattering + +To describe multiple scattering, the *GX2F* can fit the scattering angles as they were normal parameters. +Of course, fitting more parameters increases the dimensions of all matrices. +This makes it computationally more expensive to. + +But first shortly recap on multiple scattering. +To describe scattering, on a surface, only the two angles $\theta$ and $\phi$ are needed, where: +- $\theta$ is the angle between the extrapolation of the incoming trajectory and the scattered trajectory +- $\phi$ is the rotation around the extrapolation of the incoming trajectory + +This description is only valid for thin materials. +To model thicker materials, one could in theory add multiple thin materials together. +It can be shown, that it is enough to two sets of $\theta$ and $\phi$ on both sides of the material. +We could name them $\theta_{in}$, $\theta_{out}$, $\phi_{in}$, and $\phi_{out}$. +But in the end they are just multiple parameters in our fit. +That's why we will only look at $\theta$ and $\phi$ (like for thin materials). + +By defining residuals and covariances for the scattering angles, we can put them into our $\chi^2$ function. +For the residuals we choose (since the expected angle is 0) + +$$ +r_s = \begin{cases} + 0 - \theta_s(\vec\alpha) \\ + 0 - \sin(\theta_{loc})\phi_s(\vec\alpha) + \end{cases} +$$ + +with $\theta_{loc}$ the angle between incoming trajectory and normal of the surface. +(We cannot have angle information $\phi$ if we are perpendicular.) +For the covariances we use the Highland form [^cornelissen] as + +$$ +\sigma_{scat} = \frac{13.6 \text{MeV}}{\beta c p} Z\prime \sqrt{\frac{x}{X_0}} \left( 1 + 0.038 \ln{\frac{x}{X_0}} \right) +$$ + +with +- $x$ ... material layer with thickness +- $X_0$ ... radiation length +- $p$ ... particle momentum +- $Z\prime$ ... charge number +- $\beta c$ ... velocity + +Combining those terms we can write our $\chi^2$ function including multiple scattering as + +$$ +\chi^2 = \sum_{i=1}^N \frac{r_i^2}{\sigma_i^2} + \sum_{s}^S \left(\frac{\theta_s^2}{\sigma_s^2} + \frac{\sin^2{(\theta_{loc})}\phi_s^2}{\sigma_s^2}\right) +$$ + +Note, that both scattering angles have the same covariance. + +### (coming soon) Mathematical description of the energy loss [wip] :::{todo} -Write GX2F documentation +Write *GX2F*: Mathematical description of the energy loss + +The development work on the energy loss has not finished yet. ::: -[^billoir]: https://twiki.cern.ch/twiki/pub/LHCb/ParametrizedKalman/paramKalmanV01.pdf +### Implementation in ACTS + +The implementation is in some points similar to the KF, since the KF interface was chosen as a starting point. +This makes it easier to replace both fitters with each other. +The structure of the *GX2F* implementation follows coarsely the mathematical outline given above. +It is best to start reading the implementation from `fit()`: +1. Set up the fitter: + - Actor + - Aborter + - Propagator + - Variables we need longer than one iteration +2. Iterate + 1. Update parameters + 2. Propagate through geometry + 3. Collect: + - Residuals + - Covariances + - Projected Jacobians + 4. Loop over collection and calculate and sum over: + - $\chi^2$ + - $[a_{kl}]$ + - $\vec b$ + 5. Solve $[a_{kl}] \vec{\delta\alpha}_n = \vec b$ + 6. Check for convergence +3. Calculate covariance of the final parameters +4. Prepare and return the final track + +#### Configuration + +Here is a simplified example of the configuration of the fitter. + +```cpp +template +struct Gx2FitterOptions { +Gx2FitterOptions( ... ) : ... {} + +Gx2FitterOptions() = delete; + +... +//common options: +// geoContext, magFieldContext, calibrationContext, extensions, +// propagatorPlainOptions, referenceSurface, multipleScattering, +// energyLoss, freeToBoundCorrection + +/// Max number of iterations during the fit (abort condition) +size_t nUpdateMax = 5; + +/// Disables the QoP fit in case of missing B-field +bool zeroField = false; + +/// Check for convergence (abort condition). Set to 0 to skip. +double relChi2changeCutOff = 1e-7; +}; +``` + +Common options like the geometry context or toggling of the energy loss are similar to the other fitters. +For now there are three *GX2F* specific options: +1. `nUpdateMax` sets an abort condition for the parameter update as a maximum number of iterations allowed. +We do not really want to use this condition, but it stops the fit in case of poor convergence. +2. `zeroField` toggles the q/p-fit. +If there is no magnetic field, we get no q/p-information. +This would crash the fitter when needing matrix inverses. +When this option is set to `true`, most of the matrices will omit the q/p-rows and -columns. +3. `relChi2changeCutOff` is the desired convergence criterion. +We compare at each step of the iteration the current to the previous $\chi^2$. +If the relative change is small enough, we finish the fit. + +### Pros/Cons + +There are some reasons for and against the *GX2F*. +The biggest issue of the *GX2F* is its performance. +Currently, the most expensive part is the propagation. +Since we need to do a full propagation each iteration, we end up with at least 4-5 full propagation. +This is a lot compared to the 2 propagations of the *KF*. +However, since the *GX2F* is a global fitter, it can easier resolve left-right-ambiguous measurements, like in the TRT (Transition Radiation Tracker – straw tubes). + +[^billoir]: [https://twiki.cern.ch/twiki/pub/LHCb/ParametrizedKalman/paramKalmanV01.pdf](https://twiki.cern.ch/twiki/pub/LHCb/ParametrizedKalman/paramKalmanV01.pdf) +[^cornelissen]: [https://cds.cern.ch/record/1005181/files/thesis-2006-072.pdf#page=80](https://cds.cern.ch/record/1005181/files/thesis-2006-072.pdf)