Merge branch 'main' into ex-digi-particle-measurement-maps

Core/include/Acts/EventData/MultiComponentTrackParameters.hpp

@@ -129,11 +129,11 @@ class MultiComponentBoundTrackParameters {
 
   /// Get the weight and a GenericBoundTrackParameters object for one component
   std::pair<double, Parameters> operator[](std::size_t i) const {
-    return std::make_pair(
+    return {
         std::get<double>(m_components[i]),
         Parameters(m_surface, std::get<BoundVector>(m_components[i]),
                    std::get<std::optional<BoundSquareMatrix>>(m_components[i]),
-                   m_particleHypothesis));
+                   m_particleHypothesis)};
   }
 
   /// Parameters vector.

Core/include/Acts/MagneticField/MultiRangeBField.hpp

@@ -0,0 +1,67 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#pragma once
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/MagneticField/MagneticFieldError.hpp"
+#include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/Utilities/RangeXD.hpp"
+
+namespace Acts {
+
+/// @ingroup MagneticField
+///
+/// @brief Magnetic field provider modelling a magnetic field consisting of
+/// several (potentially overlapping) regions of constant values.
+class MultiRangeBField final : public MagneticFieldProvider {
+ private:
+  struct Cache {
+    explicit Cache(const MagneticFieldContext& /*unused*/);
+
+    std::optional<std::size_t> index = {};
+  };
+
+  using BFieldRange = std::pair<RangeXD<3, double>, Vector3>;
+
+  // The different ranges and their corresponding field vectors. Note that
+  // regions positioned _later_ in this vector take priority over earlier
+  // regions.
+  std::vector<BFieldRange> fieldRanges;
+
+ public:
+  /// @brief Construct a magnetic field from a vector of ranges.
+  ///
+  /// @warning These ranges are listed in increasing order of precedence,
+  /// i.e. ranges further along the vector have higher priority.
+  explicit MultiRangeBField(const std::vector<BFieldRange>& ranges);
+
+  explicit MultiRangeBField(std::vector<BFieldRange>&& ranges);
+
+  /// @brief Construct a cache object.
+  MagneticFieldProvider::Cache makeCache(
+      const MagneticFieldContext& mctx) const override;
+
+  /// @brief Request the value of the magnetic field at a given position.
+  ///
+  /// @param [in] position Global 3D position for the lookup.
+  /// @param [in, out] cache Cache object.
+  /// @returns A successful value containing a field vector if the given
+  /// location is contained inside any of the regions, or a failure value
+  /// otherwise.
+  Result<Vector3> getField(const Vector3& position,
+                           MagneticFieldProvider::Cache& cache) const override;
+
+  /// @brief Get the field gradient at a given position.
+  ///
+  /// @warning This is not currently implemented.
+  Result<Vector3> getFieldGradient(
+      const Vector3& position, ActsMatrix<3, 3>& /*unused*/,
+      MagneticFieldProvider::Cache& cache) const override;
+};
+}  // namespace Acts

Core/include/Acts/Seeding/BinnedGroupIterator.ipp

@@ -64,7 +64,7 @@ Acts::BinnedGroupIterator<grid_t>::operator*() const {
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wstringop-overread"
 #endif
-  return std::make_tuple(std::move(bottoms), global_index, std::move(tops));
+  return {std::move(bottoms), global_index, std::move(tops)};
 #if defined(__GNUC__) && __GNUC__ >= 12 && !defined(__clang__)
 #pragma GCC diagnostic pop
 #endif

Core/include/Acts/Seeding/SeedFinder.ipp

@@ -839,7 +839,7 @@ std::pair<float, float> SeedFinder<external_spacepoint_t, grid_t, platform_t>::
         const external_spacepoint_t& spM,
         const Acts::Range1D<float>& rMiddleSPRange) const {
   if (m_config.useVariableMiddleSPRange) {
-    return std::make_pair(rMiddleSPRange.min(), rMiddleSPRange.max());
+    return {rMiddleSPRange.min(), rMiddleSPRange.max()};
   }
   if (!m_config.rRangeMiddleSP.empty()) {
     /// get zBin position of the middle SP
@@ -848,10 +848,9 @@ std::pair<float, float> SeedFinder<external_spacepoint_t, grid_t, platform_t>::
     int zBin = std::distance(m_config.zBinEdges.begin(), pVal);
     /// protects against zM at the limit of zBinEdges
     zBin == 0 ? zBin : --zBin;
-    return std::make_pair(m_config.rRangeMiddleSP[zBin][0],
-                          m_config.rRangeMiddleSP[zBin][1]);
+    return {m_config.rRangeMiddleSP[zBin][0], m_config.rRangeMiddleSP[zBin][1]};
   }
-  return std::make_pair(m_config.rMinMiddle, m_config.rMaxMiddle);
+  return {m_config.rMinMiddle, m_config.rMaxMiddle};
 }
 
 }  // namespace Acts

Core/include/Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp

@@ -95,7 +95,7 @@ globalTrackParametersCovariance(const traj_t& multiTraj,
     prev_ts = ts;
   });
 
-  return std::make_pair(fullGlobalTrackParamsCov, stateRowIndices);
+  return {fullGlobalTrackParamsCov, stateRowIndices};
 }
 
 }  // namespace Acts::detail

Core/include/Acts/Utilities/AlgebraHelpers.hpp

@@ -178,7 +178,7 @@ inline ActsMatrix<A::RowsAtCompileTime, B::ColsAtCompileTime> blockedMult(
 /// Calculate the inverse of an Eigen matrix after checking if it can be
 /// numerically inverted. This allows to catch potential FPEs before they occur.
 /// For matrices up to 4x4, the inverse is computed directly. For larger
-/// matrices, the FullPivLU is used.
+/// matrices, and dynamic matrices the FullPivLU is used.
 ///
 /// @tparam Derived Eigen derived concrete type
 /// @tparam Result Eigen result type defaulted to input type
@@ -192,12 +192,11 @@ std::optional<ResultType> safeInverse(const MatrixType& m) noexcept {
   constexpr int cols = MatrixType::ColsAtCompileTime;
 
   static_assert(rows == cols);
-  static_assert(rows != -1);
 
   ResultType result;
   bool invertible = false;
 
-  if constexpr (rows > 4) {
+  if constexpr (rows > 4 || rows == -1) {
     Eigen::FullPivLU<MatrixType> mFullPivLU(m);
     if (mFullPivLU.isInvertible()) {
       invertible = true;

Core/include/Acts/Utilities/GridBinFinder.ipp

@@ -45,9 +45,9 @@ std::array<std::pair<int, int>, DIM> Acts::GridBinFinder<DIM>::getSizePerAxis(
           using value_t = typename std::decay_t<decltype(val)>;
           if constexpr (std::is_same_v<int, value_t>) {
             assert(val >= 0);
-            return std::make_pair(-val, val);
+            return {-val, val};
           } else if constexpr (std::is_same_v<std::pair<int, int>, value_t>) {
-            return std::make_pair(-val.first, val.second);
+            return {-val.first, val.second};
           } else {
             assert(locPosition.size() > i);
             assert(locPosition[i] > 0ul);

Core/include/Acts/Utilities/TrackHelpers.hpp

@@ -209,7 +209,7 @@ findTrackStateForExtrapolation(
       }
 
       ACTS_VERBOSE("found intersection at " << intersection.pathLength());
-      return std::make_pair(*first, intersection.pathLength());
+      return std::pair(*first, intersection.pathLength());
     }
 
     case TrackExtrapolationStrategy::last: {
@@ -229,7 +229,7 @@ findTrackStateForExtrapolation(
       }
 
       ACTS_VERBOSE("found intersection at " << intersection.pathLength());
-      return std::make_pair(*last, intersection.pathLength());
+      return std::pair(*last, intersection.pathLength());
     }
 
     case TrackExtrapolationStrategy::firstOrLast: {
@@ -256,13 +256,13 @@ findTrackStateForExtrapolation(
       if (intersectionFirst.isValid() && absDistanceFirst <= absDistanceLast) {
         ACTS_VERBOSE("using first track state with intersection at "
                      << intersectionFirst.pathLength());
-        return std::make_pair(*first, intersectionFirst.pathLength());
+        return std::pair(*first, intersectionFirst.pathLength());
       }
 
       if (intersectionLast.isValid() && absDistanceLast <= absDistanceFirst) {
         ACTS_VERBOSE("using last track state with intersection at "
                      << intersectionLast.pathLength());
-        return std::make_pair(*last, intersectionLast.pathLength());
+        return std::pair(*last, intersectionLast.pathLength());
       }
 
       ACTS_ERROR("no intersection found");
@@ -531,7 +531,7 @@ calculatePredictedResidual(track_state_proxy_t trackState) {
   MeasurementMatrix residualCovariance =
       measurementCovariance + predictedCovariance;
 
-  return std::pair(residual, residualCovariance);
+  return {residual, residualCovariance};
 }
 
 /// Helper function to calculate the filtered residual and its covariance
@@ -568,7 +568,7 @@ calculateFilteredResidual(track_state_proxy_t trackState) {
   MeasurementMatrix residualCovariance =
       measurementCovariance + filteredCovariance;
 
-  return std::pair(residual, residualCovariance);
+  return {residual, residualCovariance};
 }
 
 /// Helper function to calculate the smoothed residual and its covariance
@@ -605,7 +605,7 @@ calculateSmoothedResidual(track_state_proxy_t trackState) {
   MeasurementMatrix residualCovariance =
       measurementCovariance + smoothedCovariance;
 
-  return std::pair(residual, residualCovariance);
+  return {residual, residualCovariance};
 }
 
 /// Helper function to calculate the predicted chi2

Core/src/Detector/detail/IndexedGridFiller.cpp

@@ -47,39 +47,31 @@ std::vector<std::size_t> Acts::Experimental::detail::binSequence(
     rBins.reserve(bmax - bmin + 1u + 2 * expand);
     // handle bmin:/max expand it down (for bound, don't fill underflow)
     if (type == Acts::AxisBoundaryType::Bound) {
-      bmin = (static_cast<int>(bmin) - static_cast<int>(expand) > 0)
-                 ? bmin - expand
-                 : 1u;
+      bmin = bmin > expand ? bmin - expand : 1u;
       bmax = (bmax + expand <= nBins) ? bmax + expand : nBins;
     } else if (type == Acts::AxisBoundaryType::Open) {
-      bmin = (static_cast<int>(bmin) - static_cast<int>(expand) >= 0)
-                 ? bmin - expand
-                 : 0u;
+      bmin = bmin >= expand ? bmin - expand : 0u;
       bmax = (bmax + expand <= nBins + 1u) ? bmax + expand : nBins + 1u;
     }
     fill_linear(bmin, bmax);
   } else {
     // Close case
     std::size_t span = bmax - bmin + 1u + 2 * expand;
     // Safe with respect to the closure point, treat as bound
-    if (2 * span < nBins && (bmax + expand <= nBins) &&
-        (static_cast<int>(bmin) - static_cast<int>(expand) > 0)) {
+    if (2 * span < nBins && (bmax + expand <= nBins) && (bmin > expand)) {
       return binSequence({bmin, bmax}, expand, nBins,
                          Acts::AxisBoundaryType::Bound);
     } else if (2 * span < nBins) {
-      bmin = static_cast<int>(bmin) - static_cast<int>(expand) > 0
-                 ? bmin - expand
-                 : 1u;
+      bmin = bmin > expand ? bmin - expand : 1u;
       bmax = bmax + expand <= nBins ? bmax + expand : nBins;
       fill_linear(bmin, bmax);
       // deal with expansions over the phi boundary
       if (bmax + expand > nBins) {
         std::size_t overstep = (bmax + expand - nBins);
         fill_linear(1u, overstep);
       }
-      if (static_cast<int>(bmin) - static_cast<int>(expand) < 1) {
-        std::size_t understep =
-            abs(static_cast<int>(bmin) - static_cast<int>(expand));
+      if (bmin <= expand) {
+        std::size_t understep = expand - bmin;
         fill_linear(nBins - understep, nBins);
       }
       std::ranges::sort(rBins);

Core/src/MagneticField/CMakeLists.txt

@@ -1,4 +1,8 @@
 target_sources(
     ActsCore
-    PRIVATE BFieldMapUtils.cpp SolenoidBField.cpp MagneticFieldError.cpp
+    PRIVATE
+        BFieldMapUtils.cpp
+        SolenoidBField.cpp
+        MagneticFieldError.cpp
+        MultiRangeBField.cpp
 )

Core/src/MagneticField/MultiRangeBField.cpp

@@ -0,0 +1,78 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#include "Acts/MagneticField/MultiRangeBField.hpp"
+
+namespace Acts {
+
+MultiRangeBField::Cache::Cache(const MagneticFieldContext& /*unused*/) {}
+
+MultiRangeBField::MultiRangeBField(const std::vector<BFieldRange>& ranges)
+    : fieldRanges(ranges) {}
+
+MultiRangeBField::MultiRangeBField(
+    std::vector<MultiRangeBField::BFieldRange>&& ranges)
+    : fieldRanges(std::move(ranges)) {}
+
+MagneticFieldProvider::Cache MultiRangeBField::makeCache(
+    const MagneticFieldContext& mctx) const {
+  return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
+}
+
+Result<Vector3> MultiRangeBField::getField(
+    const Vector3& position, MagneticFieldProvider::Cache& cache) const {
+  // Because we assume that the regions are in increasing order of
+  // precedence, we can iterate over the array, remembering the _last_
+  // region that contained the given point. At the end of the loop, this
+  // region will then be the one we query for its field vector.
+  std::optional<std::size_t> foundRange = {};
+
+  // The cache object for this type contains an optional integer; if the
+  // integer is set, it indicates the index of the region in which the last
+  // access succeeded. This acts as a cache because if the current access is
+  // still within that region, it is guaranteed that none of the regions
+  // with lower priority -- i.e. that are earlier in the region vector --
+  // can be relevant to the current access. Thus, we request the cache index
+  // if it exists and perform a membership check on it; if that succeeds, we
+  // remember the corresponding region as a candidate.
+  if (Cache& lCache = cache.as<Cache>();
+      lCache.index.has_value() &&
+      std::get<0>(fieldRanges[*lCache.index])
+          .contains({position[0], position[1], position[2]})) {
+    foundRange = *lCache.index;
+  }
+
+  // Now, we iterate over the ranges. If we already have a range candidate,
+  // we start iteration just after the existing candidate; otherwise we
+  // start from the beginning.
+  for (std::size_t i = (foundRange.has_value() ? (*foundRange) + 1 : 0);
+       i < fieldRanges.size(); ++i) {
+    if (std::get<0>(fieldRanges[i])
+            .contains({position[0], position[1], position[2]})) {
+      foundRange = i;
+    }
+  }
+
+  // Update the cache using the result of this access.
+  cache.as<Cache>().index = foundRange;
+
+  // If we found a valid range, return the corresponding vector; otherwise
+  // return an out-of-bounds error.
+  if (foundRange.has_value()) {
+    return Result<Vector3>::success(std::get<1>(fieldRanges[*foundRange]));
+  } else {
+    return Result<Vector3>::failure(MagneticFieldError::OutOfBounds);
+  }
+}
+
+Result<Vector3> MultiRangeBField::getFieldGradient(
+    const Vector3& position, ActsMatrix<3, 3>& /*unused*/,
+    MagneticFieldProvider::Cache& cache) const {
+  return getField(position, cache);
+}
+}  // namespace Acts

Core/src/Propagator/detail/CovarianceEngine.cpp

@@ -97,8 +97,7 @@ CurvilinearState detail::curvilinearState(
       pos4, direction, freeParameters[eFreeQOverP], std::move(cov),
       particleHypothesis);
   // Create the curvilinear state
-  return std::make_tuple(std::move(curvilinearParams), fullTransportJacobian,
-                         accumulatedPath);
+  return {std::move(curvilinearParams), fullTransportJacobian, accumulatedPath};
 }
 
 void detail::transportCovarianceToBound(

Core/src/Propagator/detail/SympyCovarianceEngine.cpp

@@ -86,8 +86,7 @@ CurvilinearState sympy::curvilinearState(
       pos4, direction, freeParameters[eFreeQOverP], std::move(cov),
       particleHypothesis);
   // Create the curvilinear state
-  return std::make_tuple(std::move(curvilinearParams), fullTransportJacobian,
-                         accumulatedPath);
+  return {std::move(curvilinearParams), fullTransportJacobian, accumulatedPath};
 }
 
 void sympy::transportCovarianceToBound(

Core/src/Surfaces/detail/AlignmentHelper.cpp

@@ -87,6 +87,6 @@ Acts::detail::RotationToAxes Acts::detail::rotationToLocalAxesDerivative(
                     rotToCompositeLocalYAxis * relRotation(1, 2) +
                     rotToCompositeLocalZAxis * relRotation(2, 2);
 
-  return std::make_tuple(std::move(rotToLocalXAxis), std::move(rotToLocalYAxis),
-                         std::move(rotToLocalZAxis));
+  return {std::move(rotToLocalXAxis), std::move(rotToLocalYAxis),
+          std::move(rotToLocalZAxis)};
 }

Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp

@@ -63,5 +63,5 @@ Acts::detail::AnnulusBoundsHelper::create(const Transform3& transform,
   auto annulusBounds = std::make_shared<AnnulusBounds>(
       rMin, rMax, phiMin, phiMax, originShift, phiShift);
 
-  return std::make_tuple(annulusBounds, boundsTransform);
+  return {annulusBounds, boundsTransform};
 }

Core/src/Utilities/SpacePointUtility.cpp

@@ -99,7 +99,7 @@ SpacePointUtility::globalCoords(
     tcov = std::nullopt;
   }
 
-  return std::make_tuple(globalPos, globalTime, gcov, tcov);
+  return {globalPos, globalTime, gcov, tcov};
 }
 
 Vector2 SpacePointUtility::calcRhoZVars(