Merge branch 'main' into ML-to-Core

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/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/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/include/Acts/Visualization/Interpolation3D.hpp

@@ -0,0 +1,96 @@
+// 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 <unsupported/Eigen/Splines>
+
+namespace Acts::Interpolation3D {
+
+/// @brief Helper function to interpolate points using a spline
+/// from Eigen
+///
+/// The only requirement is that the input trajectory type has
+/// a method empty() and size() and that the elements can be
+/// accessed with operator[] and have themselves a operator[] to
+/// access the coordinates.
+///
+/// @tparam input_trajectory_type input trajectory type
+///
+/// @param inputsRaw input vector points
+/// @param nPoints number of interpolation points
+/// @param keepOriginalHits keep the original hits in the trajectory
+///
+/// @return std::vector<Acts::Vector3> interpolated points
+template <typename trajectory_type>
+trajectory_type spline(const trajectory_type& inputsRaw, std::size_t nPoints,
+                       bool keepOriginalHits = false) {
+  trajectory_type output;
+  if (inputsRaw.empty()) {
+    return output;
+  }
+
+  using InputVectorType = typename trajectory_type::value_type;
+
+  std::vector<Vector3> inputs;
+  // If input type is a vector of Vector3 we can use it directly
+  if constexpr (std::is_same_v<trajectory_type, std::vector<Vector3>>) {
+    inputs = inputsRaw;
+  } else {
+    inputs.reserve(inputsRaw.size());
+    for (const auto& input : inputsRaw) {
+      inputs.push_back(Vector3(input[0], input[1], input[2]));
+    }
+  }
+
+  // Don't do anything if we have less than 3 points or less interpolation
+  // points than input points
+  if (inputsRaw.size() < 3 || nPoints <= inputsRaw.size()) {
+    return inputsRaw;
+  } else {
+    Eigen::MatrixXd points(3, inputs.size());
+    for (std::size_t i = 0; i < inputs.size(); ++i) {
+      points.col(i) = inputs[i].transpose();
+    }
+    Eigen::Spline<double, 3> spline3D =
+        Eigen::SplineFitting<Eigen::Spline<double, 3>>::Interpolate(points, 2);
+
+    double step = 1. / (nPoints - 1);
+    for (std::size_t i = 0; i < nPoints; ++i) {
+      double t = i * step;
+      InputVectorType point;
+      point[0] = spline3D(t)[0];
+      point[1] = spline3D(t)[1];
+      point[2] = spline3D(t)[2];
+      output.push_back(point);
+    }
+  }
+  // If we want to keep the original hits, we add them to the output
+  // (first and last are there anyway)
+  if (keepOriginalHits) {
+    output.insert(output.begin(), inputsRaw.begin() + 1, inputsRaw.end() - 1);
+    // We need to sort the output in distance to first
+    std::sort(output.begin(), output.end(),
+              [&inputs](const auto& a, const auto& b) {
+                const auto ifront = inputs.front();
+                double da2 = (a[0] - ifront[0]) * (a[0] - ifront[0]) +
+                             (a[1] - ifront[1]) * (a[1] - ifront[1]) +
+                             (a[2] - ifront[2]) * (a[2] - ifront[2]);
+                double db2 = (b[0] - ifront[0]) * (b[0] - ifront[0]) +
+                             (b[1] - ifront[1]) * (b[1] - ifront[1]) +
+                             (b[2] - ifront[2]) * (b[2] - ifront[2]);
+                return da2 < db2;
+              });
+  }
+
+  return output;
+}
+
+}  // namespace Acts::Interpolation3D

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

Core/src/Vertexing/AdaptiveGridTrackDensity.cpp

@@ -156,7 +156,7 @@ AdaptiveGridTrackDensity::getMaxZTPosition(DensityMap& densityMap) const {
   double maxZ = getSpatialBinCenter(bin.first);
   double maxT = getTemporalBinCenter(bin.second);
 
-  return std::make_pair(maxZ, maxT);
+  return std::pair(maxZ, maxT);
 }
 
 Result<AdaptiveGridTrackDensity::ZTPositionAndWidth>

Core/src/Vertexing/ImpactPointEstimator.cpp

@@ -217,7 +217,7 @@ Result<std::pair<Vector4, Vector3>> getDistanceAndMomentumImpl(
     Vector4 deltaRStraightTrack{Vector4::Zero()};
     deltaRStraightTrack.head<nDim>() = pcaStraightTrack - vtxPos;
 
-    return std::make_pair(deltaRStraightTrack, momDirStraightTrack);
+    return std::pair(deltaRStraightTrack, momDirStraightTrack);
   }
 
   // Charged particles in a constant B field follow a helical trajectory. In
@@ -291,7 +291,7 @@ Result<std::pair<Vector4, Vector3>> getDistanceAndMomentumImpl(
   Vector4 deltaR{Vector4::Zero()};
   deltaR.head<nDim>() = pca - vtxPos;
 
-  return std::make_pair(deltaR, momDir);
+  return std::pair(deltaR, momDir);
 }
 
 }  // namespace

Core/src/Vertexing/Vertex.cpp

@@ -74,7 +74,7 @@ const std::vector<TrackAtVertex>& Vertex::tracks() const {
 }
 
 std::pair<double, double> Vertex::fitQuality() const {
-  return std::pair<double, double>(m_chiSquared, m_numberDoF);
+  return {m_chiSquared, m_numberDoF};
 }
 
 void Vertex::setPosition(const Vector3& position) {

Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp

@@ -48,6 +48,10 @@ class DigitizationAlgorithm final : public IAlgorithm {
     std::string outputMeasurementParticlesMap = "measurement_particles_map";
     /// Output collection to map measured hits to simulated hits.
     std::string outputMeasurementSimHitsMap = "measurement_simhits_map";
+    /// Output collection to map particles to measurements.
+    std::string outputParticleMeasurementsMap = "particle_measurements_map";
+    /// Output collection to map particles to simulated hits.
+    std::string outputSimHitMeasurementsMap = "simhit_measurements_map";
 
     /// Map of surface by identifier to allow local - to global
     std::unordered_map<Acts::GeometryIdentifier, const Acts::Surface*>
@@ -122,7 +126,7 @@ class DigitizationAlgorithm final : public IAlgorithm {
   Config m_cfg;
   /// Digitizers within geometry hierarchy
   Acts::GeometryHierarchyMap<Digitizer> m_digitizers;
-  /// Geometric digtizer
+  /// Geometric digitizer
   ActsFatras::Channelizer m_channelizer;
 
   using CellsMap =
@@ -140,6 +144,11 @@ class DigitizationAlgorithm final : public IAlgorithm {
   WriteDataHandle<IndexMultimap<Index>> m_outputMeasurementSimHitsMap{
       this, "OutputMeasurementSimHitsMap"};
 
+  WriteDataHandle<InverseMultimap<SimBarcode>> m_outputParticleMeasurementsMap{
+      this, "OutputParticleMeasurementsMap"};
+  WriteDataHandle<InverseMultimap<Index>> m_outputSimHitMeasurementsMap{
+      this, "OutputSimHitMeasurementsMap"};
+
   /// Construct a fixed-size smearer from a configuration.
   ///
   /// It's templated on the smearing dimension given by @tparam kSmearDIM
@@ -153,7 +162,7 @@ class DigitizationAlgorithm final : public IAlgorithm {
     // Copy the geometric configuration
     impl.geometric = cfg.geometricDigiConfig;
     // Prepare the smearing configuration
-    for (int i = 0; i < static_cast<int>(kSmearDIM); ++i) {
+    for (std::size_t i = 0; i < kSmearDIM; ++i) {
       impl.smearing.indices[i] = cfg.smearingDigiConfig.at(i).index;
       impl.smearing.smearFunctions[i] =
           cfg.smearingDigiConfig.at(i).smearFunction;