feat!: Wire time to spacepoints and seeds (#2829)

Currently our spacepoints and seeds do not carry any time information from a potential measurement to the initial track parameters.

In this PR I try to wire optional time information originating from the measurement, carried through spacepoints into seeds and finally into initial track params.

One caveat is that we use bound parameters in between which cannot tell if time is actually measured or not. My intermediate solution was to rely on the projected covariance entry to be 0 in this case.

CI/codespell_ignore.txt

@@ -20,3 +20,4 @@ boxs
 ans
 dthe
 dthe
+vart

Core/include/Acts/Seeding/BinnedSPGroup.ipp

@@ -138,7 +138,7 @@ Acts::BinnedSPGroup<external_spacepoint_t>::BinnedSPGroup(
       continue;
     }
     const external_spacepoint_t& sp = **it;
-    const auto& [spPosition, variance] =
+    const auto& [spPosition, variance, spTime] =
         toGlobal(sp, config.zAlign, config.rAlign, config.sigmaError);
 
     float spX = spPosition[0];
@@ -158,7 +158,7 @@ Acts::BinnedSPGroup<external_spacepoint_t>::BinnedSPGroup(
     }
 
     auto isp = std::make_unique<InternalSpacePoint<external_spacepoint_t>>(
-        counter, sp, spPosition, options.beamPos, variance);
+        counter, sp, spPosition, options.beamPos, variance, spTime);
     // calculate r-Bin index and protect against overflow (underflow not
     // possible)
     std::size_t rIndex =

Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp

@@ -148,15 +148,13 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
 /// @param bFieldMin is the minimum magnetic field required to trigger the
 /// estimation of q/pt
 /// @param logger A logger instance
-/// @param mass is the estimated particle mass
 ///
 /// @return optional bound parameters
 template <typename spacepoint_iterator_t>
 std::optional<BoundVector> estimateTrackParamsFromSeed(
     const GeometryContext& gctx, spacepoint_iterator_t spBegin,
     spacepoint_iterator_t spEnd, const Surface& surface, const Vector3& bField,
-    ActsScalar bFieldMin, const Acts::Logger& logger = getDummyLogger(),
-    ActsScalar mass = 139.57018 * UnitConstants::MeV) {
+    ActsScalar bFieldMin, const Acts::Logger& logger = getDummyLogger()) {
   // Check the number of provided space points
   std::size_t numSP = std::distance(spBegin, spEnd);
   if (numSP != 3) {
@@ -180,6 +178,8 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   // The global positions of the bottom, middle and space points
   std::array<Vector3, 3> spGlobalPositions = {Vector3::Zero(), Vector3::Zero(),
                                               Vector3::Zero()};
+  std::array<std::optional<float>, 3> spGlobalTimes = {
+      std::nullopt, std::nullopt, std::nullopt};
   // The first, second and third space point are assumed to be bottom, middle
   // and top space point, respectively
   for (std::size_t isp = 0; isp < 3; ++isp) {
@@ -190,6 +190,7 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
     }
     const auto& sp = *it;
     spGlobalPositions[isp] = Vector3(sp->x(), sp->y(), sp->z());
+    spGlobalTimes[isp] = sp->t();
   }
 
   // Define a new coordinate frame with its origin at the bottom space point, z
@@ -274,33 +275,14 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   // The estimated loc0 and loc1
   params[eBoundLoc0] = bottomLocalPos.x();
   params[eBoundLoc1] = bottomLocalPos.y();
+  params[eBoundTime] = spGlobalTimes[0].value_or(0.);
 
   // The estimated q/pt in [GeV/c]^-1 (note that the pt is the projection of
   // momentum on the transverse plane of the new frame)
   ActsScalar qOverPt = rho * (UnitConstants::m) / (0.3 * bFieldInTesla);
   // The estimated q/p in [GeV/c]^-1
   params[eBoundQOverP] = qOverPt / std::hypot(1., invTanTheta);
 
-  // The estimated momentum, and its projection along the magnetic field
-  // diretion
-  ActsScalar pInGeV = std::abs(1.0 / params[eBoundQOverP]);
-  ActsScalar pzInGeV = 1.0 / std::abs(qOverPt) * invTanTheta;
-  ActsScalar massInGeV = mass / UnitConstants::GeV;
-  // The estimated velocity, and its projection along the magnetic field
-  // diretion
-  ActsScalar v = pInGeV / std::hypot(pInGeV, massInGeV);
-  ActsScalar vz = pzInGeV / std::hypot(pInGeV, massInGeV);
-  // The z coordinate of the bottom space point along the magnetic field
-  // direction
-  ActsScalar pathz = spGlobalPositions[0].dot(bField) / bField.norm();
-  // The estimated time (use path length along magnetic field only if it's not
-  // zero)
-  if (pathz != 0 && vz != 0) {
-    params[eBoundTime] = pathz / vz;
-  } else {
-    params[eBoundTime] = spGlobalPositions[0].norm() / v;
-  }
-
   if (params.hasNaN()) {
     ACTS_ERROR(
         "The NaN value exists at the estimated track parameters from seed with"

Core/include/Acts/Seeding/InternalSpacePoint.hpp

@@ -10,12 +10,12 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 
-#include <array>
 #include <cmath>
 #include <functional>
-#include <limits>
+#include <optional>
 
 namespace Acts {
+
 template <typename SpacePoint>
 class InternalSpacePoint {
   /////////////////////////////////////////////////////////////////////////////////
@@ -27,7 +27,8 @@ class InternalSpacePoint {
   InternalSpacePoint(std::size_t index, const SpacePoint& sp,
                      const Acts::Vector3& globalPos,
                      const Acts::Vector2& offsetXY,
-                     const Acts::Vector2& variance);
+                     const Acts::Vector2& variance,
+                     std::optional<float> globalTime);
 
   InternalSpacePoint(const InternalSpacePoint<SpacePoint>& sp);
   ~InternalSpacePoint() = default;
@@ -39,6 +40,7 @@ class InternalSpacePoint {
   float x() const { return m_x; }
   float y() const { return m_y; }
   float z() const { return m_z; }
+  std::optional<float> t() const { return m_t; }
   float radius() const { return m_r; }
   float phi() const { return m_phi; }
   float varianceR() const { return m_varianceR; }
@@ -47,13 +49,14 @@ class InternalSpacePoint {
 
  protected:
   std::size_t m_index;
-  float m_x;          // x-coordinate in beam system coordinates
-  float m_y;          // y-coordinate in beam system coordinates
-  float m_z;          // z-coordinate in beam system coordinetes
-  float m_r;          // radius       in beam system coordinates
-  float m_phi;        //
-  float m_varianceR;  //
-  float m_varianceZ;  //
+  float m_x;                 // x-coordinate in beam system coordinates
+  float m_y;                 // y-coordinate in beam system coordinates
+  float m_z;                 // z-coordinate in beam system coordinetes
+  float m_r;                 // radius       in beam system coordinates
+  float m_phi;               //
+  float m_varianceR;         //
+  float m_varianceZ;         //
+  std::optional<float> m_t;  // time
   std::reference_wrapper<const SpacePoint> m_sp;  // external space point
 };
 
@@ -64,7 +67,8 @@ class InternalSpacePoint {
 template <typename SpacePoint>
 inline InternalSpacePoint<SpacePoint>::InternalSpacePoint(
     std::size_t index, const SpacePoint& sp, const Acts::Vector3& globalPos,
-    const Acts::Vector2& offsetXY, const Acts::Vector2& variance)
+    const Acts::Vector2& offsetXY, const Acts::Vector2& variance,
+    std::optional<float> globalTime)
     : m_index(index),
       m_x(globalPos.x() - offsetXY.x()),
       m_y(globalPos.y() - offsetXY.y()),
@@ -73,6 +77,7 @@ inline InternalSpacePoint<SpacePoint>::InternalSpacePoint(
       m_phi(std::atan2(m_y, m_x)),
       m_varianceR(variance.x()),
       m_varianceZ(variance.y()),
+      m_t(globalTime),
       m_sp(sp) {}
 
 /////////////////////////////////////////////////////////////////////////////////

Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp

@@ -742,9 +742,9 @@ void SeedFinderOrthogonal<external_spacepoint_t>::createSeeds(
   spacePointData.resize(spacePoints.size());
 
   for (const external_spacepoint_t *p : spacePoints) {
-    auto [position, variance] = extract_coordinates(p);
+    auto [position, variance, time] = extract_coordinates(p);
     internalSpacePoints.push_back(new InternalSpacePoint<external_spacepoint_t>(
-        counter++, *p, position, options.beamPos, variance));
+        counter++, *p, position, options.beamPos, variance, time));
     // store x,y,z values in extent
     rRangeSPExtent.extend(position);
   }

Core/include/Acts/SpacePointFormation/SpacePointBuilder.hpp

@@ -34,17 +34,17 @@ namespace Acts {
 template <typename spacepoint_t>
 class SpacePointBuilder {
  public:
+  using BuilderFunction = std::function<spacepoint_t(
+      Acts::Vector3, std::optional<double>, Acts::Vector2,
+      std::optional<double>, boost::container::static_vector<SourceLink, 2>)>;
+
   // Constructor
   /// @param cfg The configuration for the space point builder
   /// @param func The function that provides user's SP constructor with global pos, global cov, and sourceLinks.
   /// @param logger The logging instance
-  SpacePointBuilder(const SpacePointBuilderConfig& cfg,
-                    std::function<spacepoint_t(
-                        Acts::Vector3, Acts::Vector2,
-                        boost::container::static_vector<SourceLink, 2>)>
-                        func,
+  SpacePointBuilder(const SpacePointBuilderConfig& cfg, BuilderFunction func,
                     std::unique_ptr<const Logger> logger =
-                        getDefaultLogger("SpamcePointBuilder", Logging::INFO));
+                        getDefaultLogger("SpacePointBuilder", Logging::INFO));
 
   // Default constructor
   SpacePointBuilder() = default;
@@ -83,9 +83,7 @@ class SpacePointBuilder {
   /// @brief Function to create external space point
   /// The constructor of spacepoint_t with Vector3 global pos, Vector2 global
   /// cov, and vector of source link pointers.
-  std::function<spacepoint_t(Acts::Vector3, Acts::Vector2,
-                             boost::container::static_vector<SourceLink, 2>)>
-      m_spConstructor;
+  BuilderFunction m_spConstructor;
 
   /// the logging instance
   std::unique_ptr<const Acts::Logger> m_logger;
@@ -96,4 +94,5 @@ class SpacePointBuilder {
 };
 
 }  // namespace Acts
+
 #include "Acts/SpacePointFormation/detail/SpacePointBuilder.ipp"

Core/include/Acts/SpacePointFormation/SpacePointBuilderOptions.hpp

@@ -15,15 +15,17 @@
 
 namespace Acts {
 
+using ParamCovAccessor =
+    std::function<std::pair<const BoundVector, const BoundSquareMatrix>(
+        const SourceLink&)>;
+
 struct SpacePointBuilderOptions {
   // ends of strip pairs
   std::pair<const std::pair<Vector3, Vector3>,
             const std::pair<Vector3, Vector3>>
       stripEndsPair;
   // accessor of local position and covariance from source link
-  std::function<std::pair<const BoundVector, const BoundSquareMatrix>(
-      const SourceLink&)>
-      paramCovAccessor;
+  ParamCovAccessor paramCovAccessor;
   /// vertex position
   Vector3 vertex = {0., 0., 0.};
   /// Allowed increase of strip length
@@ -34,9 +36,7 @@ struct SpacePointBuilderOptions {
 
 struct StripPairOptions {
   // accessor of local position and covariance from source link
-  std::function<std::pair<const BoundVector, const BoundSquareMatrix>(
-      const SourceLink&)>
-      paramCovAccessor;
+  ParamCovAccessor paramCovAccessor;
   /// vertex position
   Vector3 vertex = {0., 0., 0.};
   /// Accepted squared difference in theta for two clusters

Core/include/Acts/SpacePointFormation/detail/SpacePointBuilder.ipp

@@ -5,14 +5,14 @@
 // 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/Definitions/Algebra.hpp"
+
 namespace Acts {
 
 template <typename spacepoint_t>
 SpacePointBuilder<spacepoint_t>::SpacePointBuilder(
-    const SpacePointBuilderConfig& cfg,
-    std::function<spacepoint_t(Acts::Vector3, Acts::Vector2,
-                               boost::container::static_vector<SourceLink, 2>)>
-        func,
+    const SpacePointBuilderConfig& cfg, BuilderFunction func,
     std::unique_ptr<const Logger> logger)
     : m_config(cfg), m_spConstructor(func), m_logger(std::move(logger)) {
   m_spUtility = std::make_shared<SpacePointUtility>(cfg);
@@ -27,15 +27,15 @@ void SpacePointBuilder<spacepoint_t>::buildSpacePoint(
   const unsigned int num_slinks = sourceLinks.size();
 
   Acts::Vector3 gPos = Acts::Vector3::Zero();
+  std::optional<Acts::ActsScalar> gTime = std::nullopt;
   Acts::Vector2 gCov = Acts::Vector2::Zero();
+  std::optional<Acts::ActsScalar> gCovT = std::nullopt;
 
   if (num_slinks == 1) {  // pixel SP formation
     auto slink = sourceLinks.at(0);
-    auto [param, cov] = opt.paramCovAccessor(sourceLinks.at(0));
-    auto gPosCov = m_spUtility->globalCoords(
+    auto [param, cov] = opt.paramCovAccessor(slink);
+    std::tie(gPos, gTime, gCov, gCovT) = m_spUtility->globalCoords(
         gctx, slink, m_config.slSurfaceAccessor, param, cov);
-    gPos = gPosCov.first;
-    gCov = gPosCov.second;
   } else if (num_slinks == 2) {  // strip SP formation
 
     const auto& ends1 = opt.stripEndsPair.first;
@@ -71,9 +71,9 @@ void SpacePointBuilder<spacepoint_t>::buildSpacePoint(
       gPos = ends1.first + resultPerpProj.value() * spParams.firstBtmToTop;
     }
 
-    double theta =
-        acos(spParams.firstBtmToTop.dot(spParams.secondBtmToTop) /
-             (spParams.firstBtmToTop.norm() * spParams.secondBtmToTop.norm()));
+    double theta = std::acos(
+        spParams.firstBtmToTop.dot(spParams.secondBtmToTop) /
+        (spParams.firstBtmToTop.norm() * spParams.secondBtmToTop.norm()));
 
     gCov = m_spUtility->calcRhoZVars(gctx, sourceLinks.at(0), sourceLinks.at(1),
                                      m_config.slSurfaceAccessor,
@@ -85,7 +85,7 @@ void SpacePointBuilder<spacepoint_t>::buildSpacePoint(
   boost::container::static_vector<SourceLink, 2> slinks(sourceLinks.begin(),
                                                         sourceLinks.end());
 
-  spacePointIt = m_spConstructor(gPos, gCov, std::move(slinks));
+  spacePointIt = m_spConstructor(gPos, gTime, gCov, gCovT, std::move(slinks));
 }
 
 template <typename spacepoint_t>
@@ -108,12 +108,15 @@ void SpacePointBuilder<spacepoint_t>::makeSourceLinkPairs(
       const auto& slinkBack = slinksBack[j];
 
       const auto [paramFront, covFront] = pairOpt.paramCovAccessor(slinkFront);
-      const auto [gposFront, gcovFront] = m_spUtility->globalCoords(
-          gctx, slinkFront, m_config.slSurfaceAccessor, paramFront, covFront);
+      const auto [gposFront, gtimeFront, gcovFront, gcovtFront] =
+          m_spUtility->globalCoords(gctx, slinkFront,
+                                    m_config.slSurfaceAccessor, paramFront,
+                                    covFront);
 
       const auto [paramBack, covBack] = pairOpt.paramCovAccessor(slinkBack);
-      const auto [gposBack, gcovBack] = m_spUtility->globalCoords(
-          gctx, slinkBack, m_config.slSurfaceAccessor, paramBack, covBack);
+      const auto [gposBack, gtimeBack, gcovBack, gcovtBack] =
+          m_spUtility->globalCoords(gctx, slinkBack, m_config.slSurfaceAccessor,
+                                    paramBack, covBack);
 
       auto res = m_spUtility->differenceOfMeasurementsChecked(
           gposFront, gposBack, pairOpt.vertex, pairOpt.diffDist,

Core/include/Acts/Utilities/SpacePointUtility.hpp

@@ -15,6 +15,7 @@
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/SpacePointFormation/SpacePointBuilderConfig.hpp"
+#include "Acts/SpacePointFormation/SpacePointBuilderOptions.hpp"
 #include "Acts/Utilities/Result.hpp"
 
 #include <array>
@@ -74,10 +75,11 @@ class SpacePointUtility {
   /// @param par local position
   /// @param cov local covariance
   /// @return vectors of the global coordinates and covariance of the SourceLink
-  std::pair<Vector3, Vector2> globalCoords(
-      const GeometryContext& gctx, const SourceLink& slink,
-      const SourceLinkSurfaceAccessor& surfaceAccessor, const BoundVector& par,
-      const BoundSquareMatrix& cov) const;
+  std::tuple<Vector3, std::optional<ActsScalar>, Vector2,
+             std::optional<ActsScalar>>
+  globalCoords(const GeometryContext& gctx, const SourceLink& slink,
+               const SourceLinkSurfaceAccessor& surfaceAccessor,
+               const BoundVector& par, const BoundSquareMatrix& cov) const;
 
   /// @brief Get rho and z covariance from the local position and covariance
   /// @param gctx The current geometry context object, e.g. alignment
@@ -98,13 +100,12 @@ class SpacePointUtility {
   /// @param globalPos global position
   /// @param theta The angle between the two strips
   /// @return (rho, z) components of the global covariance
-  Vector2 calcRhoZVars(
-      const GeometryContext& gctx, const SourceLink& slinkFront,
-      const SourceLink& slinkBack,
-      const SourceLinkSurfaceAccessor& surfaceAccessor,
-      const std::function<std::pair<const BoundVector, const BoundSquareMatrix>(
-          SourceLink)>& paramCovAccessor,
-      const Vector3& globalPos, const double theta) const;
+  Vector2 calcRhoZVars(const GeometryContext& gctx,
+                       const SourceLink& slinkFront,
+                       const SourceLink& slinkBack,
+                       const SourceLinkSurfaceAccessor& surfaceAccessor,
+                       const ParamCovAccessor& paramCovAccessor,
+                       const Vector3& globalPos, const double theta) const;
 
   /// @brief This function performs a straight forward calculation of a space
   /// point and returns whether it was successful or not.