Merge branch 'fix-g4-detector-with-fatras' of github.com:andiwand/act…

…s into fix-g4-detector-with-fatras

CI/physmon/physmon.py

@@ -200,7 +200,6 @@ def run_ckf_tracking(truthSmearedSeeded, truthEstimatedSeeded, label):
             field,
             CKFPerformanceConfig(ptMin=400.0 * u.MeV, nMeasurementsMin=6),
             TrackSelectorRanges(
-                removeNeutral=True,
                 loc0=(-4.0 * u.mm, 4.0 * u.mm),
                 pt=(500 * u.MeV, None),
             ),
@@ -316,7 +315,6 @@ def run_vertexing(fitter, mu, events):
                 pt=(500 * u.MeV, None),
                 loc0=(-4.0 * u.mm, 4.0 * u.mm),
                 absEta=(None, 3.0),
-                removeNeutral=True,
             ),
         )
 

Core/include/Acts/EventData/Track.hpp

@@ -8,8 +8,13 @@
 
 #pragma once
 
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/EventData/Charge.hpp"
 #include "Acts/EventData/MultiTrajectory.hpp"
 #include "Acts/Utilities/HashedString.hpp"
+#include "Acts/Utilities/UnitVectors.hpp"
 
 #include <any>
 #include <cstddef>
@@ -206,6 +211,58 @@ class TrackProxy {
     return m_container->covariance(m_index);
   }
 
+  ActsScalar charge() const {
+    // Currently, neutral tracks are not supported here
+    // @TODO: Evaluate if/how neutral 'tracks' should be accounted for
+    return SinglyCharged{}.extractCharge(parameters()[eBoundQOverP]);
+  }
+
+  /// Access the theta parameter of the track at the reference surface
+  /// @return The theta parameter
+  ActsScalar theta() const { return parameters()[eBoundTheta]; }
+
+  /// Access the phi parameter of the track at the reference surface
+  /// @return The phi parameter
+  ActsScalar phi() const { return parameters()[eBoundPhi]; }
+
+  /// Access the loc0 parameter of the track at the reference surface
+  /// @return The loc0 parameter
+  ActsScalar loc0() const { return parameters()[eBoundLoc0]; }
+
+  /// Access the loc1 parameter of the track at the reference surface
+  /// @return The loc1 parameter
+  ActsScalar loc1() const { return parameters()[eBoundLoc1]; }
+
+  /// Access the time parameter of the track at the reference surface
+  /// @return The time parameter
+  ActsScalar time() const { return parameters()[eBoundTime]; }
+
+  /// Access the q/p (curvature) parameter of the track at the reference surface
+  /// @return The q/p parameter
+  ActsScalar qOverP() const { return parameters()[eBoundQOverP]; }
+
+  /// Get the absolute momentum of the tack
+  /// @return The absolute track momentum
+  ActsScalar absoluteMomentum() const {
+    return SinglyCharged{}.extractMomentum(qOverP());
+  }
+
+  /// Get the transverse momentum of the track
+  /// @return The track transverse momentum value
+  ActsScalar transverseMomentum() const {
+    return std::sin(theta()) * absoluteMomentum();
+  }
+
+  /// Get a unit vector along the track direction at the reference surface
+  /// @return The direction unit vector
+  Vector3 unitDirection() const {
+    return makeDirectionUnitFromPhiTheta(phi(), theta());
+  }
+
+  /// Get the global momentum vector
+  /// @return the global momentum vector
+  Vector3 momentum() const { return absoluteMomentum() * unitDirection(); }
+
   /// Get a range over the track states of this track. Return value is
   /// compatible with range based for loop. Const version
   /// @return Track state range to iterate over
@@ -522,6 +579,14 @@ class TrackContainer {
     return m_container->addTrack_impl();
   }
 
+  /// Remove a track at index @p itrack from the container
+  /// @note This invalidates all track proxies!
+  /// @param itrack The index of the track to remmove
+  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
+  void removeTrack(IndexType itrack) {
+    m_container->removeTrack_impl(itrack);
+  }
+
   /// Add a dymanic column to the track container
   /// @param key the name of the column to be added
   template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
@@ -559,10 +624,21 @@ class TrackContainer {
     return *m_traj;
   }
 
+  /// Retrieve the holder of the track state container
+  /// @return The track state container including it's holder
+  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
+  auto& trackStateContainerHolder() {
+    return m_traj;
+  }
+
   /// Get a const reference to the track state container backend
   /// @return a const reference to the backend
   const auto& trackStateContainer() const { return *m_traj; }
 
+  /// Retrieve the holder of the track state container
+  /// @return The track state container including it's holder
+  const auto& trackStateContainerHolder() const { return m_traj; }
+
   /// Get a mutable iterator to the first track in the container
   /// @return a mutable iterator to the first track
   template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>

Core/include/Acts/EventData/VectorTrackContainer.hpp

@@ -37,15 +37,7 @@ class VectorTrackContainerBase {
  protected:
   VectorTrackContainerBase() = default;
 
-  VectorTrackContainerBase(const VectorTrackContainerBase& other)
-      : m_tipIndex{other.m_tipIndex},
-        m_params{other.m_params},
-        m_cov{other.m_cov},
-        m_referenceSurfaces{other.m_referenceSurfaces} {
-    for (const auto& [key, value] : other.m_dynamic) {
-      m_dynamic.insert({key, value->clone()});
-    }
-  };
+  VectorTrackContainerBase(const VectorTrackContainerBase& other);
 
   VectorTrackContainerBase(VectorTrackContainerBase&& other) = default;
 
@@ -88,6 +80,31 @@ class VectorTrackContainerBase {
     }
   }
 
+  bool checkConsistency() const {
+    size_t size = m_tipIndex.size();
+    (void)size;
+
+    bool result = true;
+    result = result && m_tipIndex.size() == size;
+    assert(result);
+    result = result && m_params.size() == size;
+    assert(result);
+    result = result && m_cov.size() == size;
+    assert(result);
+    result = result && m_referenceSurfaces.size() == size;
+    assert(result);
+    result = result && m_nMeasurements.size() == size;
+    assert(result);
+    result = result && m_nHoles.size() == size;
+
+    for (const auto& [key, col] : m_dynamic) {
+      (void)key;
+      result = result && col->size() == size;
+    }
+
+    return result;
+  }
+
  public:
   constexpr bool hasColumn_impl(HashedString key) const {
     using namespace Acts::HashedStringLiteral;
@@ -97,7 +114,10 @@ class VectorTrackContainerBase {
     }
   }
 
-  std::size_t size_impl() const { return m_tipIndex.size(); }
+  std::size_t size_impl() const {
+    assert(checkConsistency());
+    return m_tipIndex.size();
+  }
   // END INTERFACE HELPER
 
   std::vector<IndexType> m_tipIndex;
@@ -115,6 +135,8 @@ class VectorTrackContainerBase {
 }  // namespace detail_vtc
 
 class VectorTrackContainer;
+class ConstVectorTrackContainer;
+
 template <>
 struct IsReadOnlyTrackContainer<VectorTrackContainer> : std::false_type {};
 
@@ -124,6 +146,8 @@ class VectorTrackContainer final : public detail_vtc::VectorTrackContainerBase {
   VectorTrackContainer(const VectorTrackContainer& other) = default;
   VectorTrackContainer(VectorTrackContainer&&) = default;
 
+  VectorTrackContainer(const ConstVectorTrackContainer& other);
+
  public:
   // BEGIN INTERFACE
 
@@ -137,23 +161,9 @@ class VectorTrackContainer final : public detail_vtc::VectorTrackContainerBase {
         *this, key, itrack);
   }
 
-  IndexType addTrack_impl() {
-    m_tipIndex.emplace_back();
+  IndexType addTrack_impl();
 
-    m_params.emplace_back();
-    m_cov.emplace_back();
-    m_referenceSurfaces.emplace_back();
-
-    m_nMeasurements.emplace_back();
-    m_nHoles.emplace_back();
-
-    // dynamic columns
-    for (auto& [key, vec] : m_dynamic) {
-      vec->add();
-    }
-
-    return m_tipIndex.size() - 1;
-  }
+  void removeTrack_impl(IndexType itrack);
 
   template <typename T>
   constexpr void addColumn_impl(const std::string& key) {
@@ -191,11 +201,15 @@ class ConstVectorTrackContainer final
 
   ConstVectorTrackContainer(const ConstVectorTrackContainer& other) = default;
   ConstVectorTrackContainer(const VectorTrackContainer& other)
-      : VectorTrackContainerBase{other} {}
+      : VectorTrackContainerBase{other} {
+    assert(checkConsistency());
+  }
 
   ConstVectorTrackContainer(ConstVectorTrackContainer&&) = default;
   ConstVectorTrackContainer(VectorTrackContainer&& other)
-      : VectorTrackContainerBase{std::move(other)} {}
+      : VectorTrackContainerBase{std::move(other)} {
+    assert(checkConsistency());
+  }
 
  public:
   // BEGIN INTERFACE
@@ -216,4 +230,10 @@ class ConstVectorTrackContainer final
   // END INTERFACE
 };
 
+inline VectorTrackContainer::VectorTrackContainer(
+    const ConstVectorTrackContainer& other)
+    : VectorTrackContainerBase{other} {
+  assert(checkConsistency());
+}
+
 }  // namespace Acts

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

@@ -21,6 +21,8 @@ struct DynamicColumnBase {
 
   virtual void add() = 0;
   virtual void clear() = 0;
+  virtual void erase(size_t i) = 0;
+  virtual size_t size() const = 0;
 
   virtual std::unique_ptr<DynamicColumnBase> clone() const = 0;
 };
@@ -39,6 +41,8 @@ struct DynamicColumn : public DynamicColumnBase {
 
   void add() override { m_vector.emplace_back(); }
   void clear() override { m_vector.clear(); }
+  void erase(size_t i) override { m_vector.erase(m_vector.begin() + i); }
+  size_t size() const override { return m_vector.size(); }
 
   std::unique_ptr<DynamicColumnBase> clone() const override {
     return std::make_unique<DynamicColumn<T>>(*this);
@@ -65,6 +69,8 @@ struct DynamicColumn<bool> : public DynamicColumnBase {
 
   void add() override { m_vector.emplace_back(); }
   void clear() override { m_vector.clear(); }
+  void erase(size_t i) override { m_vector.erase(m_vector.begin() + i); }
+  size_t size() const override { return m_vector.size(); }
 
   std::unique_ptr<DynamicColumnBase> clone() const override {
     return std::make_unique<DynamicColumn<bool>>(*this);

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -164,7 +164,8 @@ class EigenStepper {
   };
 
   /// Constructor requires knowledge of the detector's magnetic field
-  EigenStepper(std::shared_ptr<const MagneticFieldProvider> bField);
+  EigenStepper(std::shared_ptr<const MagneticFieldProvider> bField,
+               double overstepLimit = 100 * UnitConstants::um);
 
   template <typename charge_t>
   State makeState(std::reference_wrapper<const GeometryContext> gctx,
@@ -402,8 +403,8 @@ class EigenStepper {
   /// Magnetic field inside of the detector
   std::shared_ptr<const MagneticFieldProvider> m_bField;
 
-  /// Overstep limit: could/should be dynamic
-  double m_overstepLimit = 100 * UnitConstants::um;
+  /// Overstep limit
+  double m_overstepLimit;
 };
 }  // namespace Acts
 

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -11,8 +11,9 @@
 
 template <typename E, typename A>
 Acts::EigenStepper<E, A>::EigenStepper(
-    std::shared_ptr<const MagneticFieldProvider> bField)
-    : m_bField(std::move(bField)) {}
+
+    std::shared_ptr<const MagneticFieldProvider> bField, double overstepLimit)
+    : m_bField(std::move(bField)), m_overstepLimit(overstepLimit) {}
 
 template <typename E, typename A>
 template <typename charge_t>

Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp

@@ -28,7 +28,7 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
     return SingleStepper::boundState(cmpIt->state, surface, transportCov,
                                      freeToBoundCorrection);
   } else {
-    SmallVector<std::pair<double, BoundTrackParameters>> states;
+    SmallVector<std::tuple<double, BoundVector, BoundSymMatrix>> states;
     double accumulatedPathLength = 0.0;
     int failedBoundTransforms = 0;
 
@@ -37,8 +37,10 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
                                           transportCov, freeToBoundCorrection);
 
       if (bs.ok()) {
-        states.push_back(
-            {state.components[i].weight, std::get<BoundTrackParameters>(*bs)});
+        const auto& btp = std::get<BoundTrackParameters>(*bs);
+        states.emplace_back(
+            state.components[i].weight, btp.parameters(),
+            btp.covariance().value_or(Acts::BoundSymMatrix::Zero()));
         accumulatedPathLength +=
             std::get<double>(*bs) * state.components[i].weight;
       } else {
@@ -54,16 +56,16 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
       return MultiStepperError::SomeComponentsConversionToBoundFailed;
     }
 
-    const auto proj = [&](const auto& wbs) {
-      return std::tie(wbs.first, wbs.second.parameters(),
-                      wbs.second.covariance());
-    };
-
     auto [params, cov] =
         detail::angleDescriptionSwitch(surface, [&](const auto& desc) {
-          return detail::combineGaussianMixture(states, proj, desc);
+          return detail::combineGaussianMixture(states, Acts::Identity{}, desc);
         });
 
+    std::optional<BoundSymMatrix> finalCov = std::nullopt;
+    if (cov != BoundSymMatrix::Zero()) {
+      finalCov = cov;
+    }
+
     return BoundState{BoundTrackParameters(surface.getSharedPtr(), params, cov),
                       Jacobian::Zero(), accumulatedPathLength};
   }
@@ -90,6 +92,7 @@ auto MultiEigenStepperLoop<E, R, A>::curvilinearState(State& state,
     ActsScalar qop = 0.0;
     BoundSymMatrix cov = BoundSymMatrix::Zero();
     ActsScalar pathLenth = 0.0;
+    ActsScalar sumOfWeights = 0.0;
 
     for (auto i = 0ul; i < numberComponents(state); ++i) {
       const auto [cp, jac, pl] = SingleStepper::curvilinearState(
@@ -102,10 +105,23 @@ auto MultiEigenStepperLoop<E, R, A>::curvilinearState(State& state,
         cov += state.components[i].weight * *cp.covariance();
       }
       pathLenth += state.components[i].weight * pathLenth;
+      sumOfWeights += state.components[i].weight;
+    }
+
+    pos4 /= sumOfWeights;
+    dir /= sumOfWeights;
+    qop /= sumOfWeights;
+    pathLenth /= sumOfWeights;
+    cov /= sumOfWeights;
+
+    std::optional<BoundSymMatrix> finalCov = std::nullopt;
+    if (cov != BoundSymMatrix::Zero()) {
+      finalCov = cov;
     }
 
-    return CurvilinearState{CurvilinearTrackParameters(pos4, dir, qop, cov),
-                            Jacobian::Zero(), pathLenth};
+    return CurvilinearState{
+        CurvilinearTrackParameters(pos4, dir, qop, finalCov), Jacobian::Zero(),
+        pathLenth};
   }
 }