Merge branch 'main' of github.com:acts-project/acts into remove-min-f…

…ield-from-param-est

.zenodo.json

@@ -112,11 +112,15 @@
       "affiliation": "CERN / University of Amsterdam",
       "name": "Stephen Nicholas Swatman",
       "orcid": "0000-0002-3747-3229"
-    }
+    },
     {
       "affiliation": "CERN / TU Wien",
       "name": "Felix Russo",
       "orcid": "0009-0005-8975-2245"
+    },
+    {
+      "affiliation": "UC Berkeley",
+      "name": "Carlo Varni"
     }
   ],
   "access_right": "open",

AUTHORS

@@ -21,3 +21,4 @@ The following people have contributed to the project (in alphabetical order):
 - Andreas Stefl, CERN, TU Wien
 - Stephen Nicholas Swatman, CERN, University of Amsterdam
 - Roman Urmanov, Weizmann Institute of Science
+- Carlo Varni, UC Berkeley

CITATION.cff

@@ -89,6 +89,9 @@ authors:
   family-names: Russo
   affiliation: CERN / TU Wien
   orcid: https://orcid.org/0009-0005-8975-2245
+- given-names: Carlo
+  family-names: Varni
+  affiliation: UC Berkeley
 version: 10.0.0
 date-released: 2021-07-28
 repository-code: https://github.com/acts-project/acts

Core/include/Acts/EventData/ProxyAccessor.hpp

@@ -75,7 +75,7 @@ struct ProxyAccessorBase {
   /// Create the accessor from a string key
   /// @param _key the key
   constexpr ProxyAccessorBase(const std::string& _key)
-      : key{hashString(_key)} {}
+      : key{hashStringDynamic(_key)} {}
 
   /// Access the stored key on the proxy given as an argument. Mutable version
   /// @tparam proxy_t the type of the proxy

Core/include/Acts/EventData/TrackContainer.hpp

@@ -252,7 +252,7 @@ class TrackContainer {
   /// Check if this track container has a specific dynamic column
   /// @param key the key to check for
   constexpr bool hasColumn(const std::string& key) const {
-    return m_container->hasColumn_impl(hashString(key));
+    return m_container->hasColumn_impl(hashStringDynamic(key));
   }
 
   /// Check if a this track container has a specific dynamic column

Core/include/Acts/EventData/TrackProxy.hpp

@@ -710,7 +710,7 @@ class TrackProxy {
   constexpr T& component(std::string_view key)
     requires(!ReadOnly)
   {
-    return m_container->template component<T>(hashString(key), m_index);
+    return m_container->template component<T>(hashStringDynamic(key), m_index);
   }
 
   /// Retrieve a const reference to a component
@@ -738,7 +738,7 @@ class TrackProxy {
   /// @return Const reference to the component given by @p key
   template <typename T>
   constexpr const T& component(std::string_view key) const {
-    return m_container->template component<T>(hashString(key), m_index);
+    return m_container->template component<T>(hashStringDynamic(key), m_index);
   }
 
   /// @}

Core/include/Acts/EventData/TrackStateProxy.hpp

@@ -1101,7 +1101,7 @@ class TrackStateProxy {
   /// @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));
+    return has(hashStringDynamic(key));
   }
 
   /// Retrieve a mutable reference to a component
@@ -1135,7 +1135,7 @@ class TrackStateProxy {
   constexpr T& component(std::string_view key)
     requires(!ReadOnly)
   {
-    return m_traj->template component<T>(hashString(key), m_istate);
+    return m_traj->template component<T>(hashStringDynamic(key), m_istate);
   }
 
   /// Retrieve a const reference to a component
@@ -1163,7 +1163,7 @@ class TrackStateProxy {
   /// @return Const reference to the component given by @p key
   template <typename T>
   constexpr const T& component(std::string_view key) const {
-    return m_traj->template component<T>(hashString(key), m_istate);
+    return m_traj->template component<T>(hashStringDynamic(key), m_istate);
   }
 
   /// @}

Core/include/Acts/EventData/VectorMultiTrajectory.hpp

@@ -456,7 +456,7 @@ class VectorMultiTrajectory final
 
   template <typename T>
   void addColumn_impl(std::string_view key) {
-    HashedString hashedKey = hashString(key);
+    HashedString hashedKey = hashStringDynamic(key);
     m_dynamic.insert({hashedKey, std::make_unique<detail::DynamicColumn<T>>()});
   }
 

Core/include/Acts/EventData/VectorTrackContainer.hpp

@@ -225,7 +225,7 @@ class VectorTrackContainer final : public detail_vtc::VectorTrackContainerBase {
 
   template <typename T>
   constexpr void addColumn_impl(const std::string_view& key) {
-    HashedString hashedKey = hashString(key);
+    HashedString hashedKey = hashStringDynamic(key);
     m_dynamic.insert({hashedKey, std::make_unique<detail::DynamicColumn<T>>()});
   }
 

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

@@ -1103,8 +1103,8 @@ class MultiTrajectoryTestsCommon {
     auto test = [&](const std::string& col, auto value) {
       using T = decltype(value);
       std::string col2 = col + "_2";
-      HashedString h{hashString(col)};
-      HashedString h2{hashString(col2)};
+      HashedString h{hashStringDynamic(col)};
+      HashedString h2{hashStringDynamic(col2)};
 
       trajectory_t traj = m_factory.create();
       BOOST_CHECK(!traj.hasColumn(h));
@@ -1188,7 +1188,7 @@ class MultiTrajectoryTestsCommon {
       }
     };
 
-    runTest([](const std::string& c) { return hashString(c.c_str()); });
+    runTest([](const std::string& c) { return hashStringDynamic(c.c_str()); });
     // runTest([](const std::string& c) { return c.c_str(); });
     // runTest([](const std::string& c) { return c; });
     // runTest([](std::string_view c) { return c; });

Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp

@@ -14,14 +14,94 @@
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
+#include "Acts/Utilities/Zip.hpp"
 
 #include <array>
 #include <cmath>
 #include <iterator>
 #include <optional>
+#include <stdexcept>
 
 namespace Acts {
 
+/// Estimate the full track parameters from three space points
+///
+/// This method is based on the conformal map transformation. It estimates the
+/// full free track parameters, i.e. (x, y, z, t, dx, dy, dz, q/p) at the
+/// bottom space point. The bottom space is assumed to be the first element
+/// in the range defined by the iterators. The magnetic field (which might be
+/// along any direction) is also necessary for the momentum estimation.
+///
+/// This is a purely spatial estimation, i.e. the time parameter will be set to
+/// 0.
+///
+/// It resembles the method used in ATLAS for the track parameters
+/// estimated from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane
+/// here:
+/// https://acode-browser.usatlas.bnl.gov/lxr/source/athena/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
+///
+/// @tparam spacepoint_iterator_t  The type of space point iterator
+///
+/// @param sp0 is the bottom space point
+/// @param sp1 is the middle space point
+/// @param sp2 is the top space point
+/// @param bField is the magnetic field vector
+///
+/// @return the free parameters
+FreeVector estimateTrackParamsFromSeed(const Vector3& sp0, const Vector3& sp1,
+                                       const Vector3& sp2,
+                                       const Vector3& bField);
+
+/// Estimate the full track parameters from three space points
+///
+/// This method is based on the conformal map transformation. It estimates the
+/// full free track parameters, i.e. (x, y, z, t, dx, dy, dz, q/p) at the
+/// bottom space point. The bottom space is assumed to be the first element
+/// in the range defined by the iterators. The magnetic field (which might be
+/// along any direction) is also necessary for the momentum estimation.
+///
+/// It resembles the method used in ATLAS for the track parameters
+/// estimated from seed, i.e. the function InDet::SiTrackMaker_xk::getAtaPlane
+/// here:
+/// https://acode-browser.usatlas.bnl.gov/lxr/source/athena/InnerDetector/InDetRecTools/SiTrackMakerTool_xk/src/SiTrackMaker_xk.cxx
+///
+/// @tparam spacepoint_iterator_t  The type of space point iterator
+///
+/// @param spRange is the range of space points
+/// @param bField is the magnetic field vector
+///
+/// @return the free parameters
+template <std::ranges::range spacepoint_range_t>
+FreeVector estimateTrackParamsFromSeed(spacepoint_range_t spRange,
+                                       const Vector3& bField) {
+  // Check the number of provided space points
+  if (spRange.size() != 3) {
+    throw std::invalid_argument(
+        "There should be exactly three space points provided.");
+  }
+
+  // The global positions of the bottom, middle and space points
+  std::array<Vector3, 3> spPositions = {Vector3::Zero(), Vector3::Zero(),
+                                        Vector3::Zero()};
+  std::array<std::optional<double>, 3> spTimes = {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 (auto [sp, spPosition, spTime] :
+       Acts::zip(spRange, spPositions, spTimes)) {
+    if (sp == nullptr) {
+      throw std::invalid_argument("Empty space point found.");
+    }
+    spPosition = Vector3(sp->x(), sp->y(), sp->z());
+    spTime = sp->t();
+  }
+
+  FreeVector params = estimateTrackParamsFromSeed(
+      spPositions[0], spPositions[1], spPositions[2], bField);
+  params[eFreeTime] = spTimes[0].value_or(0);
+  return params;
+}
+
 /// Estimate the full track parameters from three space points
 ///
 /// This method is based on the conformal map transformation. It estimates the
@@ -60,9 +140,6 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
     return std::nullopt;
   }
 
-  // Convert bField to Tesla
-  ActsScalar bFieldInTesla = bField.norm() / UnitConstants::T;
-
   // The global positions of the bottom, middle and space points
   std::array<Vector3, 3> spGlobalPositions = {Vector3::Zero(), Vector3::Zero(),
                                               Vector3::Zero()};
@@ -129,7 +206,7 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   int sign = ia > 0 ? -1 : 1;
   const ActsScalar R = circleCenter.norm();
   ActsScalar invTanTheta =
-      local2.z() / (2.f * R * std::asin(local2.head<2>().norm() / (2.f * R)));
+      local2.z() / (2 * R * std::asin(local2.head<2>().norm() / (2 * R)));
   // The momentum direction in the new frame (the center of the circle has the
   // coordinate (-1.*A/(2*B), 1./(2*B)))
   ActsScalar A = -circleCenter(0) / circleCenter(1);
@@ -159,9 +236,10 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   params[eBoundLoc1] = bottomLocalPos.y();
   params[eBoundTime] = spGlobalTimes[0].value_or(0.);
 
+  ActsScalar bFieldStrength = bField.norm();
   // 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 = sign * (UnitConstants::m) / (0.3 * bFieldInTesla * R);
+  ActsScalar qOverPt = sign / (bFieldStrength * R);
   // The estimated q/p in [GeV/c]^-1
   params[eBoundQOverP] = qOverPt / fastHypot(1., invTanTheta);
 

Core/include/Acts/Utilities/HashedString.hpp

@@ -8,9 +8,11 @@
 
 #pragma once
 
+#include <csignal>
 #include <cstddef>
 #include <cstdint>
 #include <string_view>
+#include <type_traits>
 #include <utility>
 
 namespace Acts {
@@ -35,7 +37,11 @@ constexpr int length(const char* str) {
 }
 }  // namespace detail
 
-constexpr HashedString hashString(std::string_view s) {
+consteval HashedString hashString(std::string_view s) {
+  return detail::fnv1a_32(s);
+}
+
+constexpr HashedString hashStringDynamic(std::string_view s) {
   return detail::fnv1a_32(s);
 }
 

Core/src/Geometry/TrackingVolume.cpp

@@ -96,6 +96,13 @@ const TrackingVolume* TrackingVolume::lowestTrackingVolume(
     }
   }
 
+  // @TODO: Abstract this into an accelerateable structure
+  for (const auto& volume : volumes()) {
+    if (volume.inside(position, tol)) {
+      return volume.lowestTrackingVolume(gctx, position, tol);
+    }
+  }
+
   // there is no lower sub structure
   return this;
 }

Core/src/Seeding/EstimateTrackParamsFromSeed.cpp

@@ -12,6 +12,84 @@
 
 #include <numbers>
 
+Acts::FreeVector Acts::estimateTrackParamsFromSeed(const Vector3& sp0,
+                                                   const Vector3& sp1,
+                                                   const Vector3& sp2,
+                                                   const Vector3& bField) {
+  // Define a new coordinate frame with its origin at the bottom space point, z
+  // axis long the magnetic field direction and y axis perpendicular to vector
+  // from the bottom to middle space point. Hence, the projection of the middle
+  // space point on the transverse plane will be located at the x axis of the
+  // new frame.
+  Vector3 relVec = sp1 - sp0;
+  Vector3 newZAxis = bField.normalized();
+  Vector3 newYAxis = newZAxis.cross(relVec).normalized();
+  Vector3 newXAxis = newYAxis.cross(newZAxis);
+  RotationMatrix3 rotation;
+  rotation.col(0) = newXAxis;
+  rotation.col(1) = newYAxis;
+  rotation.col(2) = newZAxis;
+  // The center of the new frame is at the bottom space point
+  Translation3 trans(sp0);
+  // The transform which constructs the new frame
+  Transform3 transform(trans * rotation);
+
+  // The coordinate of the middle and top space point in the new frame
+  Vector3 local1 = transform.inverse() * sp1;
+  Vector3 local2 = transform.inverse() * sp2;
+
+  // In the new frame the bottom sp is at the origin, while the middle
+  // sp in along the x axis. As such, the x-coordinate of the circle is
+  // at: x-middle / 2.
+  // The y coordinate can be found by using the straight line passing
+  // between the mid point between the middle and top sp and perpendicular to
+  // the line connecting them
+  Vector2 circleCenter;
+  circleCenter(0) = 0.5 * local1(0);
+
+  ActsScalar deltaX21 = local2(0) - local1(0);
+  ActsScalar sumX21 = local2(0) + local1(0);
+  // straight line connecting the two points
+  // y = a * x + c (we don't care about c right now)
+  // we simply need the slope
+  // we compute 1./a since this is what we need for the following computation
+  ActsScalar ia = deltaX21 / local2(1);
+  // Perpendicular line is then y = -1/a *x + b
+  // we can evaluate b given we know a already by imposing
+  // the line passes through P = (0.5 * (x2 + x1), 0.5 * y2)
+  ActsScalar b = 0.5 * (local2(1) + ia * sumX21);
+  circleCenter(1) = -ia * circleCenter(0) + b;
+  // Radius is a signed distance between circleCenter and first sp, which is at
+  // (0, 0) in the new frame. Sign depends on the slope a (positive vs negative)
+  int sign = ia > 0 ? -1 : 1;
+  const ActsScalar R = circleCenter.norm();
+  ActsScalar invTanTheta =
+      local2.z() / (2 * R * std::asin(local2.head<2>().norm() / (2 * R)));
+  // The momentum direction in the new frame (the center of the circle has the
+  // coordinate (-1.*A/(2*B), 1./(2*B)))
+  ActsScalar A = -circleCenter(0) / circleCenter(1);
+  Vector3 transDirection(1., A, fastHypot(1, A) * invTanTheta);
+  // Transform it back to the original frame
+  Vector3 direction = rotation * transDirection.normalized();
+
+  // Initialize the free parameters vector
+  FreeVector params = FreeVector::Zero();
+
+  // The bottom space point position
+  params.segment<3>(eFreePos0) = sp0;
+
+  // The estimated direction
+  params.segment<3>(eFreeDir0) = direction;
+
+  // 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 = sign / (bField.norm() * R);
+  // The estimated q/p in [GeV/c]^-1
+  params[eFreeQOverP] = qOverPt / fastHypot(1., invTanTheta);
+
+  return params;
+}
+
 Acts::BoundMatrix Acts::estimateTrackParamCovariance(
     const EstimateTrackParamCovarianceConfig& config, const BoundVector& params,
     bool hasTime) {

Examples/Io/Root/src/RootMaterialTrackReader.cpp

@@ -89,6 +89,10 @@ RootMaterialTrackReader::RootMaterialTrackReader(const Config& config,
 
   // Sort the entry numbers of the events
   {
+    // necessary to guarantee that m_inputChain->GetV1() is valid for the
+    // entire range
+    m_inputChain->SetEstimate(nentries);
+
     m_entryNumbers.resize(nentries);
     m_inputChain->Draw("event_id", "", "goff");
     RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(),