Merge branch 'main' into feature/gsf-clamp-bha-to-range

Core/CMakeLists.txt

@@ -110,6 +110,7 @@ add_subdirectory(src/MagneticField)
 add_subdirectory(src/Material)
 add_subdirectory(src/Navigation)
 add_subdirectory(src/Propagator)
+add_subdirectory(src/Seeding)
 add_subdirectory(src/Surfaces)
 add_subdirectory(src/TrackFinding)
 add_subdirectory(src/TrackFitting)

Core/include/Acts/EventData/GenericBoundTrackParameters.hpp

@@ -250,6 +250,17 @@ class GenericBoundTrackParameters {
     return m_surface->referenceFrame(geoCtx, position(geoCtx), momentum());
   }
 
+  /// Reflect the parameters in place.
+  void reflectInPlace() { m_params = reflectBoundParameters(m_params); }
+
+  /// Reflect the parameters.
+  /// @return Reflected parameters.
+  GenericBoundTrackParameters<ParticleHypothesis> reflect() const {
+    GenericBoundTrackParameters<ParticleHypothesis> reflected = *this;
+    reflected.reflectInPlace();
+    return reflected;
+  }
+
  private:
   BoundVector m_params;
   std::optional<BoundSquareMatrix> m_cov;

Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp

@@ -111,6 +111,14 @@ class GenericCurvilinearTrackParameters
   Vector3 position() const {
     return GenericBoundTrackParameters<ParticleHypothesis>::position({});
   }
+
+  /// Reflect the parameters.
+  /// @return Reflected parameters.
+  GenericCurvilinearTrackParameters<ParticleHypothesis> reflect() const {
+    GenericCurvilinearTrackParameters<ParticleHypothesis> reflected = *this;
+    reflected.reflectInPlace();
+    return reflected;
+  }
 };
 
 }  // namespace Acts

Core/include/Acts/EventData/GenericFreeTrackParameters.hpp

@@ -12,9 +12,11 @@
 #include "Acts/Definitions/Common.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/EventData/TrackParametersConcept.hpp"
+#include "Acts/EventData/TransformationHelpers.hpp"
 #include "Acts/EventData/detail/PrintParameters.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
 #include "Acts/Utilities/UnitVectors.hpp"
+#include "Acts/Utilities/VectorHelpers.hpp"
 
 #include <cassert>
 #include <cmath>
@@ -55,6 +57,29 @@ class GenericFreeTrackParameters {
         m_cov(std::move(cov)),
         m_particleHypothesis(std::move(particleHypothesis)) {}
 
+  /// Construct from four-position, direction, absolute momentum, and charge.
+  ///
+  /// @param pos4 Track position/time four-vector
+  /// @param dir Track direction three-vector; normalization is ignored.
+  /// @param qOverP Charge over momentum
+  /// @param cov Free parameters covariance matrix
+  /// @param particleHypothesis Particle hypothesis
+  GenericFreeTrackParameters(const Vector4& pos4, const Vector3& dir,
+                             Scalar qOverP, std::optional<CovarianceMatrix> cov,
+                             ParticleHypothesis particleHypothesis)
+      : m_params(FreeVector::Zero()),
+        m_cov(std::move(cov)),
+        m_particleHypothesis(std::move(particleHypothesis)) {
+    m_params[eFreePos0] = pos4[ePos0];
+    m_params[eFreePos1] = pos4[ePos1];
+    m_params[eFreePos2] = pos4[ePos2];
+    m_params[eFreeTime] = pos4[eTime];
+    m_params[eFreeDir0] = dir[eMom0];
+    m_params[eFreeDir1] = dir[eMom1];
+    m_params[eFreeDir2] = dir[eMom2];
+    m_params[eFreeQOverP] = qOverP;
+  }
+
   /// Construct from four-position, angles, absolute momentum, and charge.
   ///
   /// @param pos4 Track position/time four-vector
@@ -135,9 +160,9 @@ class GenericFreeTrackParameters {
   Scalar time() const { return m_params[eFreeTime]; }
 
   /// Phi direction.
-  Scalar phi() const { return phi(direction()); }
+  Scalar phi() const { return VectorHelpers::phi(direction()); }
   /// Theta direction.
-  Scalar theta() const { return theta(direction()); }
+  Scalar theta() const { return VectorHelpers::theta(direction()); }
   /// Charge over momentum.
   Scalar qOverP() const { return m_params[eFreeQOverP]; }
 
@@ -175,6 +200,17 @@ class GenericFreeTrackParameters {
     return m_particleHypothesis;
   }
 
+  /// Reflect the parameters in place.
+  void reflectInPlace() { m_params = reflectFreeParameters(m_params); }
+
+  /// Reflect the parameters.
+  /// @return Reflected parameters.
+  GenericFreeTrackParameters<ParticleHypothesis> reflect() const {
+    GenericFreeTrackParameters<ParticleHypothesis> reflected = *this;
+    reflected.reflectInPlace();
+    return reflected;
+  }
+
  private:
   FreeVector m_params;
   std::optional<FreeSquareMatrix> m_cov;

Core/include/Acts/EventData/TrackProxy.hpp

@@ -324,7 +324,8 @@ class TrackProxy {
     return std::distance(tsRange.begin(), tsRange.end());
   }
 
-  /// Return the number of measurements for the track. Const version
+  /// Return a mutable reference to the number of measurements for the track.
+  /// Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return The number of measurements
   unsigned int& nMeasurements()
@@ -333,8 +334,7 @@ class TrackProxy {
     return component<unsigned int, hashString("nMeasurements")>();
   }
 
-  /// Return a mutable reference to the number of measurements for the track.
-  /// Mutable version
+  /// Return the number of measurements for the track. Const version
   /// @return The number of measurements
   unsigned int nMeasurements() const {
     return component<unsigned int, hashString("nMeasurements")>();

Core/include/Acts/EventData/TransformationHelpers.hpp

@@ -13,11 +13,39 @@
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Utilities/Result.hpp"
+#include "Acts/Utilities/detail/periodic.hpp"
 
 namespace Acts {
 
 class Surface;
 
+/// Reflect bound track parameters.
+///
+/// @param boundParams Bound track parameters vector
+/// @return Reflected bound track parameters vector
+inline BoundVector reflectBoundParameters(const BoundVector& boundParams) {
+  BoundVector reflected = boundParams;
+  auto [phi, theta] = detail::normalizePhiTheta(
+      boundParams[eBoundPhi] - M_PI, M_PI - boundParams[eBoundTheta]);
+  reflected[eBoundPhi] = phi;
+  reflected[eBoundTheta] = theta;
+  reflected[eBoundQOverP] = -boundParams[eBoundQOverP];
+  return reflected;
+}
+
+/// Reflect free track parameters.
+///
+/// @param freeParams Free track parameters vector
+/// @return Reflected free track parameters vector
+inline FreeVector reflectFreeParameters(const FreeVector& freeParams) {
+  FreeVector reflected = freeParams;
+  reflected[eFreeDir0] = -freeParams[eFreeDir0];
+  reflected[eFreeDir1] = -freeParams[eFreeDir1];
+  reflected[eFreeDir2] = -freeParams[eFreeDir2];
+  reflected[eFreeQOverP] = -freeParams[eFreeQOverP];
+  return reflected;
+}
+
 /// Transform bound track parameters into equivalent free track parameters.
 ///
 /// @param surface Surface onto which the input parameters are bound

Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp

@@ -287,4 +287,41 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   return params;
 }
 
+/// Configuration for the estimation of the covariance matrix of the track
+/// parameters with `estimateTrackParamCovariance`.
+struct EstimateTrackParamCovarianceConfig {
+  /// The initial sigmas for the track parameters
+  BoundVector initialSigmas = {1. * UnitConstants::mm,
+                               1. * UnitConstants::mm,
+                               1. * UnitConstants::degree,
+                               1. * UnitConstants::degree,
+                               1. * UnitConstants::e / UnitConstants::GeV,
+                               1. * UnitConstants::ns};
+
+  /// The initial relative uncertainty of the q/pt
+  double initialSigmaPtRel = 0.1;
+
+  /// The inflation factors for the variances of the track parameters
+  BoundVector initialVarInflation = {1., 1., 1., 1., 1., 1.};
+  /// The inflation factor for time uncertainty if the time parameter was not
+  /// estimated
+  double noTimeVarInflation = 100.;
+};
+
+/// Estimate the covariance matrix of the given track parameters based on the
+/// provided configuration. The assumption is that we can model the uncertainty
+/// of the track parameters as a diagonal matrix with the provided initial
+/// sigmas. The inflation factors are used to inflate the initial variances
+/// based on the provided configuration. The uncertainty of q/p is estimated
+/// based on the relative uncertainty of the q/pt and the theta uncertainty.
+///
+/// @param config is the configuration for the estimation
+/// @param params is the track parameters
+/// @param hasTime is true if the track parameters have time
+///
+/// @return the covariance matrix of the track parameters
+BoundMatrix estimateTrackParamCovariance(
+    const EstimateTrackParamCovarianceConfig& config, const BoundVector& params,
+    bool hasTime);
+
 }  // namespace Acts

Core/include/Acts/TrackFinding/TrackSelector.hpp

@@ -142,15 +142,15 @@ class TrackSelector {
     std::vector<Config> cutSets = {};
 
     /// Eta bin edges for varying cuts by eta
-    std::vector<double> absEtaEdges = {};
+    std::vector<double> absEtaEdges = {0, inf};
 
     /// Get the number of eta bins
     /// @return Number of eta bins
     std::size_t nEtaBins() const { return absEtaEdges.size() - 1; }
 
     /// Construct an empty (accepts everything) configuration.
     /// Results in a single cut set and one abs eta bin from 0 to infinity.
-    EtaBinnedConfig() : cutSets{{}}, absEtaEdges{{0, inf}} {};
+    EtaBinnedConfig() : cutSets{{}} {};
 
     /// Constructor to create a config object that is not upper-bounded.
     /// This is useful to use the "fluent" API to populate the configuration.
@@ -163,13 +163,12 @@ class TrackSelector {
     /// @param absEtaEdgesIn is the vector of eta bin edges
     EtaBinnedConfig(std::vector<double> absEtaEdgesIn)
         : absEtaEdges{std::move(absEtaEdgesIn)} {
-      cutSets.resize(absEtaEdges.size() - 1);
+      cutSets.resize(nEtaBins());
     }
 
     /// Auto-converting constructor from a single cut configuration.
     /// Results in a single absolute eta bin from 0 to infinity.
-    EtaBinnedConfig(Config cutSet)
-        : cutSets{std::move(cutSet)}, absEtaEdges{{0, inf}} {}
+    EtaBinnedConfig(Config cutSet) : cutSets{std::move(cutSet)} {}
 
     /// Add a new eta bin with the given upper bound.
     /// @param etaMax Upper bound of the new eta bin
@@ -195,11 +194,17 @@ class TrackSelector {
     /// @return True if the configuration has a bin for the given eta
     bool hasCuts(double eta) const;
 
-    /// Get the index of the eta bin for a given eta
+    /// Get the index of the eta bin for a given eta.
+    /// throws an exception if Eta is outside the abs eta bin edges.
     /// @param eta Eta value
     /// @return Index of the eta bin
     std::size_t binIndex(double eta) const;
 
+    /// Get the index of the eta bin for a given eta
+    /// @param eta Eta value
+    /// @return Index of the eta bin, or >= nEtaBins() if Eta is outside the abs eta bin edges.
+    std::size_t binIndexNoCheck(double eta) const;
+
     /// Get the cuts for a given eta
     /// @param eta Eta value
     /// @return Cuts for the given eta
@@ -237,8 +242,7 @@ class TrackSelector {
 
  private:
   EtaBinnedConfig m_cfg;
-  bool m_isUnbinned;
-  bool m_noEtaCuts;
+  bool m_isUnbinned = false;
 };
 
 inline TrackSelector::Config& TrackSelector::Config::loc0(double min,
@@ -350,14 +354,22 @@ inline bool TrackSelector::EtaBinnedConfig::hasCuts(double eta) const {
 }
 
 inline std::size_t TrackSelector::EtaBinnedConfig::binIndex(double eta) const {
-  if (!hasCuts(eta)) {
+  std::size_t index = binIndexNoCheck(eta);
+  if (!(index < nEtaBins())) {
     throw std::invalid_argument{"Eta is outside the abs eta bin edges"};
   }
+  return index;
+}
 
+inline std::size_t TrackSelector::EtaBinnedConfig::binIndexNoCheck(
+    double eta) const {
   auto binIt =
       std::upper_bound(absEtaEdges.begin(), absEtaEdges.end(), std::abs(eta));
-  std::size_t index = std::distance(absEtaEdges.begin(), binIt) - 1;
-  return index;
+  std::size_t index = std::distance(absEtaEdges.begin(), binIt);
+  if (index == 0) {
+    index = absEtaEdges.size() + 1;  // positive value to check for underflow
+  }
+  return index - 1;
 }
 
 inline const TrackSelector::Config& TrackSelector::EtaBinnedConfig::getCuts(
@@ -428,8 +440,8 @@ bool TrackSelector::isValidTrack(const track_proxy_t& track) const {
 
   return track.hasReferenceSurface() &&
          within(track.transverseMomentum(), cuts.ptMin, cuts.ptMax) &&
-         (m_noEtaCuts || (within(absEta(), cuts.absEtaMin, cuts.absEtaMax) &&
-                          within(_eta, cuts.etaMin, cuts.etaMax))) &&
+         (!m_isUnbinned || (within(absEta(), cuts.absEtaMin, cuts.absEtaMax) &&
+                            within(_eta, cuts.etaMin, cuts.etaMax))) &&
          within(track.phi(), cuts.phiMin, cuts.phiMax) &&
          within(track.loc0(), cuts.loc0Min, cuts.loc0Max) &&
          within(track.loc1(), cuts.loc1Min, cuts.loc1Max) &&
@@ -452,26 +464,19 @@ inline TrackSelector::TrackSelector(
         "TrackSelector cut / eta bin configuration is inconsistent"};
   }
 
-  m_isUnbinned = false;
   if (m_cfg.nEtaBins() == 1) {
     static const std::vector<double> infVec = {0, inf};
-    bool limitEta = m_cfg.absEtaEdges != infVec;
-    m_isUnbinned = !limitEta;  // single bin, no eta edges given
-
-    const Config& cuts = m_cfg.cutSets[0];
-
-    if (limitEta && (cuts.etaMin != -inf || cuts.etaMax != inf ||
-                     cuts.absEtaMin != 0.0 || cuts.absEtaMax != inf)) {
-      throw std::invalid_argument{
-          "Explicit eta cuts are only valid for single eta bin"};
-    }
+    m_isUnbinned =
+        m_cfg.absEtaEdges == infVec;  // single bin, no eta edges given
   }
 
-  m_noEtaCuts = m_isUnbinned;
-  for (const auto& cuts : m_cfg.cutSets) {
-    if (cuts.etaMin != -inf || cuts.etaMax != inf || cuts.absEtaMin != 0.0 ||
-        cuts.absEtaMax != inf) {
-      m_noEtaCuts = false;
+  if (!m_isUnbinned) {
+    for (const auto& cuts : m_cfg.cutSets) {
+      if (cuts.etaMin != -inf || cuts.etaMax != inf || cuts.absEtaMin != 0.0 ||
+          cuts.absEtaMax != inf) {
+        throw std::invalid_argument{
+            "Explicit eta cuts are only valid for single eta bin"};
+      }
     }
   }
 }

Core/src/Seeding/CMakeLists.txt

@@ -0,0 +1 @@
+target_sources(ActsCore PRIVATE EstimateTrackParamsFromSeed.cpp)