Merge branch 'main' into return

.zenodo.json

@@ -120,7 +120,8 @@
     },
     {
       "affiliation": "UC Berkeley",
-      "name": "Carlo Varni"
+      "name": "Carlo Varni",
+      "orcid": "0000-0001-6733-4310"
     }
   ],
   "access_right": "open",

Core/include/Acts/EventData/MultiTrajectory.hpp

@@ -696,7 +696,7 @@ class MultiTrajectory {
 
     visit_measurement(measdim, [this, istate]<std::size_t DIM>(
                                    std::integral_constant<std::size_t, DIM>) {
-      self().template allocateCalibrated_impl(
+      self().allocateCalibrated_impl(
           istate, ActsVector<DIM>{ActsVector<DIM>::Zero()},
           ActsSquareMatrix<DIM>{ActsSquareMatrix<DIM>::Zero()});
     });
@@ -705,7 +705,7 @@ class MultiTrajectory {
   template <std::size_t measdim, typename val_t, typename cov_t>
   void allocateCalibrated(IndexType istate, const Eigen::DenseBase<val_t>& val,
                           const Eigen::DenseBase<cov_t>& cov) {
-    self().template allocateCalibrated_impl(istate, val, cov);
+    self().allocateCalibrated_impl(istate, val, cov);
   }
 
   void setUncalibratedSourceLink(IndexType istate, SourceLink&& sourceLink)

Core/include/Acts/EventData/MultiTrajectoryBackendConcept.hpp

@@ -131,13 +131,13 @@ concept MutableMultiTrajectoryBackend =
       { v.template addColumn_impl<double>(col) };
 
       {
-        v.template allocateCalibrated_impl(istate, ActsVector<1>{},
-                                           ActsSquareMatrix<1>{})
+        v.allocateCalibrated_impl(istate, ActsVector<1>{},
+                                  ActsSquareMatrix<1>{})
       };
       // Assuming intermediate values also work
       {
-        v.template allocateCalibrated_impl(istate, ActsVector<eBoundSize>{},
-                                           ActsSquareMatrix<eBoundSize>{})
+        v.allocateCalibrated_impl(istate, ActsVector<eBoundSize>{},
+                                  ActsSquareMatrix<eBoundSize>{})
       };
 
       { v.setUncalibratedSourceLink_impl(istate, std::move(sl)) };

Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp

@@ -1225,7 +1225,7 @@ class MultiTrajectoryTestsCommon {
 
     auto [par, cov] = generateBoundParametersCovariance(rng, {});
 
-    ts.template allocateCalibrated(par.head<3>(), cov.topLeftCorner<3, 3>());
+    ts.allocateCalibrated(par.head<3>(), cov.topLeftCorner<3, 3>());
 
     BOOST_CHECK_EQUAL(ts.calibratedSize(), 3);
     BOOST_CHECK_EQUAL(ts.template calibrated<3>(), par.head<3>());
@@ -1239,17 +1239,17 @@ class MultiTrajectoryTestsCommon {
     BOOST_CHECK_EQUAL(ts.template calibratedCovariance<3>(),
                       ActsSquareMatrix<3>::Zero());
 
-    ts.template allocateCalibrated(par2.head<3>(), cov2.topLeftCorner<3, 3>());
+    ts.allocateCalibrated(par2.head<3>(), cov2.topLeftCorner<3, 3>());
     BOOST_CHECK_EQUAL(ts.calibratedSize(), 3);
     // The values are re-assigned
     BOOST_CHECK_EQUAL(ts.template calibrated<3>(), par2.head<3>());
     BOOST_CHECK_EQUAL(ts.template calibratedCovariance<3>(),
                       (cov2.topLeftCorner<3, 3>()));
 
     // Re-allocation with a different measurement dimension is an error
-    BOOST_CHECK_THROW(ts.template allocateCalibrated(
-                          par2.head<4>(), cov2.topLeftCorner<4, 4>()),
-                      std::invalid_argument);
+    BOOST_CHECK_THROW(
+        ts.allocateCalibrated(par2.head<4>(), cov2.topLeftCorner<4, 4>()),
+        std::invalid_argument);
   }
 };
 }  // namespace Acts::detail::Test

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/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/GridAccessHelpers.hpp

@@ -212,7 +212,7 @@ class LocalSubspace : public IBoundToGridLocal {
   }
 };
 
-class BoundCylinderToZPhi final : public IBoundToGridLocal {
+class BoundCylinderToZPhi : public IBoundToGridLocal {
  public:
   double radius = 1.;
   double shift = 0.;

Core/include/Acts/Utilities/TrackHelpers.hpp

@@ -517,10 +517,11 @@ calculatePredictedResidual(track_state_proxy_t trackState) {
   }
 
   auto subspaceHelper =
-      trackState.template fixedBoundSubspaceHelper<nMeasurementDim>();
+      trackState.template projectorSubspaceHelper<nMeasurementDim>();
 
-  auto measurement = trackState.calibrated();
-  auto measurementCovariance = trackState.calibratedCovariance();
+  auto measurement = trackState.template calibrated<nMeasurementDim>();
+  auto measurementCovariance =
+      trackState.template calibratedCovariance<nMeasurementDim>();
   MeasurementVector predicted =
       subspaceHelper.projectVector(trackState.predicted());
   MeasurementMatrix predictedCovariance =
@@ -553,10 +554,11 @@ calculateFilteredResidual(track_state_proxy_t trackState) {
   }
 
   auto subspaceHelper =
-      trackState.template fixedBoundSubspaceHelper<nMeasurementDim>();
+      trackState.template projectorSubspaceHelper<nMeasurementDim>();
 
-  auto measurement = trackState.calibrated();
-  auto measurementCovariance = trackState.calibratedCovariance();
+  auto measurement = trackState.template calibrated<nMeasurementDim>();
+  auto measurementCovariance =
+      trackState.template calibratedCovariance<nMeasurementDim>();
   MeasurementVector filtered =
       subspaceHelper.projectVector(trackState.filtered());
   MeasurementMatrix filteredCovariance =
@@ -589,10 +591,11 @@ calculateSmoothedResidual(track_state_proxy_t trackState) {
   }
 
   auto subspaceHelper =
-      trackState.template fixedBoundSubspaceHelper<nMeasurementDim>();
+      trackState.template projectorSubspaceHelper<nMeasurementDim>();
 
-  auto measurement = trackState.calibrated();
-  auto measurementCovariance = trackState.calibratedCovariance();
+  auto measurement = trackState.template calibrated<nMeasurementDim>();
+  auto measurementCovariance =
+      trackState.template calibratedCovariance<nMeasurementDim>();
   MeasurementVector smoothed =
       subspaceHelper.projectVector(trackState.smoothed());
   MeasurementMatrix smoothedCovariance =
@@ -620,11 +623,13 @@ double calculatePredictedChi2(track_state_proxy_t trackState) {
 
   return visit_measurement(
       trackState.calibratedSize(),
-      [&]<std::size_t measdim>(std::integral_constant<std::size_t, measdim>) {
+      [&]<std::size_t measdim>(
+          std::integral_constant<std::size_t, measdim>) -> double {
         auto [residual, residualCovariance] =
             calculatePredictedResidual<measdim>(trackState);
 
-        return residual.transpose() * residualCovariance.inverse() * residual;
+        return (residual.transpose() * residualCovariance.inverse() * residual)
+            .eval()(0, 0);
       });
 }
 
@@ -643,11 +648,13 @@ double calculateFilteredChi2(track_state_proxy_t trackState) {
 
   return visit_measurement(
       trackState.calibratedSize(),
-      [&]<std::size_t measdim>(std::integral_constant<std::size_t, measdim>) {
+      [&]<std::size_t measdim>(
+          std::integral_constant<std::size_t, measdim>) -> double {
         auto [residual, residualCovariance] =
             calculateFilteredResidual<measdim>(trackState);
 
-        return residual.transpose() * residualCovariance.inverse() * residual;
+        return (residual.transpose() * residualCovariance.inverse() * residual)
+            .eval()(0, 0);
       });
 }
 
@@ -666,11 +673,13 @@ double calculateSmoothedChi2(track_state_proxy_t trackState) {
 
   return visit_measurement(
       trackState.calibratedSize(),
-      [&]<std::size_t measdim>(std::integral_constant<std::size_t, measdim>) {
+      [&]<std::size_t measdim>(
+          std::integral_constant<std::size_t, measdim>) -> double {
         auto [residual, residualCovariance] =
             calculateSmoothedResidual<measdim>(trackState);
 
-        return residual.transpose() * residualCovariance.inverse() * residual;
+        return (residual.transpose() * residualCovariance.inverse() * residual)
+            .eval()(0, 0);
       });
 }
 

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
 )