Combine ExaTrkX CKF (this is just to see the diff)

.gitignore

@@ -35,3 +35,9 @@ build*
 !.gitlab-ci.yml
 !.clang-tidy
 !.oclint
+
+Examples/Scripts/Python/csv
+Examples/Scripts/Python/onnx_models
+Examples/Scripts/Python/torchscript_models/
+*.old
+*.old2

CMakeLists.txt

@@ -42,6 +42,8 @@ option(ACTS_BUILD_PLUGIN_GEANT4 "Build Geant4 plugin" OFF)
 option(ACTS_BUILD_PLUGIN_EXATRKX "Build the Exa.TrkX plugin" OFF)
 option(ACTS_EXATRKX_ENABLE_ONNX "Build the Onnx backend for the exatrkx plugin" OFF)
 option(ACTS_EXATRKX_ENABLE_TORCH "Build the torchscript backend for the exatrkx plugin" ON)
+option(ACTS_EXATRKX_ENABLE_STUB "Only build a stub backend for the exatrkx plugin" OFF)
+option(ACTS_USE_SYSTEM_ACTSDD4HEP "Use the ActsDD4hep glue library provided by the system instead of building it" OFF)
 option(ACTS_BUILD_PLUGIN_IDENTIFICATION "Build Identification plugin" OFF)
 option(ACTS_BUILD_PLUGIN_JSON "Build json plugin" OFF)
 option(ACTS_USE_SYSTEM_NLOHMANN_JSON "Use nlohmann::json provided by the system instead of the bundled version" ${ACTS_USE_SYSTEM_LIBS})

Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp

@@ -312,7 +312,7 @@ class CombinatorialKalmanFilter {
         m_updaterLogger{m_logger->cloneWithSuffix("Updater")},
         m_smootherLogger{m_logger->cloneWithSuffix("Smoother")} {}
 
- private:
+ public:
   using KalmanNavigator = typename propagator_t::Navigator;
 
   /// The propagator for the transport and material update

Core/include/Acts/Utilities/ContainerPrinter.hpp

@@ -0,0 +1,53 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2023 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 http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <limits>
+#include <ostream>
+
+namespace Acts {
+
+template <typename iterator_t>
+struct ContainerPrinter {
+  iterator_t begin;
+  iterator_t end;
+
+  template <typename container_t>
+  ContainerPrinter(const container_t &c, std::size_t max) {
+    begin = c.cbegin();
+    const auto n = std::min(c.size(), max);
+    end = begin;
+    std::advance(end, n);
+    assert(std::distance(begin, end) == static_cast<std::ptrdiff_t>(n));
+  }
+
+  template <typename container_t>
+  ContainerPrinter(const container_t &c) : ContainerPrinter(c, c.size()) {}
+
+  ContainerPrinter(iterator_t a, iterator_t b) : begin(a), end(b) {}
+};
+
+template <typename iterator_t>
+std::ostream &operator<<(std::ostream &os,
+                         const ContainerPrinter<iterator_t> &p) {
+  for (auto it = p.begin; it != p.end; ++it) {
+    os << *it << " ";
+  }
+  return os;
+}
+
+template <class container_t>
+ContainerPrinter(const container_t &, std::size_t)
+    -> ContainerPrinter<typename std::decay_t<container_t>::const_iterator>;
+
+template <class container_t>
+ContainerPrinter(const container_t &)
+    -> ContainerPrinter<typename std::decay_t<container_t>::const_iterator>;
+
+}  // namespace Acts

Core/include/Acts/Utilities/KDTree.hpp

@@ -12,6 +12,7 @@
 
 #include <algorithm>
 #include <array>
+#include <cassert>
 #include <cmath>
 #include <functional>
 #include <memory>

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

@@ -26,8 +26,7 @@
 #include "ActsExamples/Framework/ProcessCode.hpp"
 #include "ActsExamples/Framework/RandomNumbers.hpp"
 #include "ActsFatras/Digitization/Channelizer.hpp"
-#include "ActsFatras/Digitization/PlanarSurfaceDrift.hpp"
-#include "ActsFatras/Digitization/PlanarSurfaceMask.hpp"
+#include "ActsFatras/Digitization/Segmentizer.hpp"
 #include "ActsFatras/Digitization/UncorrelatedHitSmearer.hpp"
 
 #include <cstddef>
@@ -69,20 +68,6 @@ class DigitizationAlgorithm final : public IAlgorithm {
   const DigitizationConfig& config() const { return m_cfg; }
 
  private:
-  /// Helper method for the geometric channelizing part
-  ///
-  /// @param geoCfg is the geometric digitization configuration
-  /// @param hit the Simultated hit
-  /// @param surface the Surface on which this is supposed to happen
-  /// @param gctx the Geometry context
-  /// @param rng the Random number engine for the drift smearing
-  ///
-  /// @return the list of channels
-  std::vector<ActsFatras::Channelizer::ChannelSegment> channelizing(
-      const GeometricConfig& geoCfg, const SimHit& hit,
-      const Acts::Surface& surface, const Acts::GeometryContext& gctx,
-      RandomEngine& rng) const;
-
   /// Helper method for creating digitized parameters from clusters
   ///
   /// @todo ADD random smearing
@@ -93,7 +78,7 @@ class DigitizationAlgorithm final : public IAlgorithm {
   /// @return the list of digitized parameters
   DigitizedParameters localParameters(
       const GeometricConfig& geoCfg,
-      const std::vector<ActsFatras::Channelizer::ChannelSegment>& channels,
+      const std::vector<ActsFatras::Segmentizer::ChannelSegment>& channels,
       RandomEngine& rng) const;
 
   /// Nested smearer struct that holds geometric digitizer and smearing
@@ -113,9 +98,7 @@ class DigitizationAlgorithm final : public IAlgorithm {
   DigitizationConfig m_cfg;
   /// Digitizers within geometry hierarchy
   Acts::GeometryHierarchyMap<Digitizer> m_digitizers;
-  /// Geometric digtizers
-  ActsFatras::PlanarSurfaceDrift m_surfaceDrift;
-  ActsFatras::PlanarSurfaceMask m_surfaceMask;
+  /// Geometric digtizer
   ActsFatras::Channelizer m_channelizer;
 
   ReadDataHandle<SimHitContainer> m_simContainerReadHandle{this,

Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp

@@ -205,16 +205,20 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
               ACTS_VERBOSE("Configured to geometric digitize "
                            << digitizer.geometric.indices.size()
                            << " parameters.");
-              auto channels = channelizing(digitizer.geometric, simHit,
-                                           *surfacePtr, ctx.geoContext, rng);
-              if (channels.empty()) {
+              const auto& cfg = digitizer.geometric;
+              Acts::Vector3 driftDir = cfg.drift(simHit.position(), rng);
+              auto channelsRes = m_channelizer.channelize(
+                  simHit, *surfacePtr, ctx.geoContext, driftDir,
+                  cfg.segmentation, cfg.thickness);
+              if (!channelsRes.ok() || channelsRes->empty()) {
                 ACTS_DEBUG(
                     "Geometric channelization did not work, skipping this hit.")
                 continue;
               }
-              ACTS_VERBOSE("Activated " << channels.size()
+              ACTS_VERBOSE("Activated " << channelsRes->size()
                                         << " channels for this hit.");
-              dParameters = localParameters(digitizer.geometric, channels, rng);
+              dParameters =
+                  localParameters(digitizer.geometric, *channelsRes, rng);
             }
 
             // Smearing part - (optionally) rest
@@ -292,30 +296,10 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
   return ProcessCode::SUCCESS;
 }
 
-std::vector<ActsFatras::Channelizer::ChannelSegment>
-ActsExamples::DigitizationAlgorithm::channelizing(
-    const GeometricConfig& geoCfg, const SimHit& hit,
-    const Acts::Surface& surface, const Acts::GeometryContext& gctx,
-    RandomEngine& rng) const {
-  Acts::Vector3 driftDir = geoCfg.drift(hit.position(), rng);
-
-  auto driftedSegment =
-      m_surfaceDrift.toReadout(gctx, surface, geoCfg.thickness, hit.position(),
-                               hit.direction(), driftDir);
-  auto maskedSegmentRes = m_surfaceMask.apply(surface, driftedSegment);
-  if (maskedSegmentRes.ok()) {
-    auto maskedSegment = maskedSegmentRes.value();
-    // Now Channelize
-    return m_channelizer.segments(gctx, surface, geoCfg.segmentation,
-                                  maskedSegment);
-  }
-  return {};
-}
-
 ActsExamples::DigitizedParameters
 ActsExamples::DigitizationAlgorithm::localParameters(
     const GeometricConfig& geoCfg,
-    const std::vector<ActsFatras::Channelizer::ChannelSegment>& channels,
+    const std::vector<ActsFatras::Segmentizer::ChannelSegment>& channels,
     RandomEngine& rng) const {
   DigitizedParameters dParameters;
 

Examples/Algorithms/Geant4/src/SensitiveSteppingAction.cpp

@@ -28,6 +28,7 @@
 #include <G4Track.hh>
 #include <G4UnitsTable.hh>
 #include <G4VPhysicalVolume.hh>
+#include <boost/describe.hpp>
 
 class G4PrimaryParticle;
 

Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp

@@ -141,6 +141,10 @@ class TrackFindingAlgorithm final : public IAlgorithm {
   mutable std::atomic<size_t> m_nTotalSeeds{0};
   mutable std::atomic<size_t> m_nFailedSeeds{0};
 
+  mutable std::mutex m_mutex;
+  mutable std::vector<unsigned> m_nTracksPerSeeds;
+  mutable std::vector<unsigned> m_nSelTracksPerSeeds;
+
   mutable tbb::combinable<Acts::VectorMultiTrajectory::Statistics>
       m_memoryStatistics{[]() {
         auto mtj = std::make_shared<Acts::VectorMultiTrajectory>();

Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/CombinedKfAndCkf.hpp

@@ -0,0 +1,137 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2022 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 http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
+#include "Acts/TrackFitting/KalmanFitter.hpp"
+
+namespace Acts {
+
+template <typename propagator_t, typename traj_t>
+struct CombinedKfAndCkf {
+  propagator_t m_propagator;
+  std::unique_ptr<const Logger> m_logger;
+  KalmanFitter<propagator_t, traj_t> m_kalmanFitter;
+
+  CombinedKfAndCkf(propagator_t pPropagator,
+                   std::unique_ptr<const Logger> _logger =
+                       getDefaultLogger("KfCkfComb", Logging::INFO))
+      : m_propagator(std::move(pPropagator)),
+        m_logger(std::move(_logger)),
+        m_kalmanFitter(pPropagator, m_logger->cloneWithSuffix(":KF")) {}
+
+  const Logger& logger() const { return *m_logger; }
+
+  template <typename sli_kf_t, typename sli_ckf_t, typename start_parameters_t,
+            typename track_container_t>
+  auto runKalmanFitter(
+      sli_kf_t it, sli_kf_t end, const start_parameters_t& sParameters,
+      const CombinatorialKalmanFilterOptions<sli_ckf_t, traj_t>& ckfOptions,
+      track_container_t& trackContainer) const {
+    KalmanFitterExtensions<traj_t> extensions;
+    extensions.calibrator = ckfOptions.extensions.calibrator;
+    extensions.updater = ckfOptions.extensions.updater;
+    extensions.smoother = ckfOptions.extensions.smoother;
+    //     extensions.outlierFinder = ckfOptions.extensions.outlierFinder;
+
+    KalmanFitterOptions<traj_t> kfOptions(
+        ckfOptions.geoContext, ckfOptions.magFieldContext,
+        ckfOptions.calibrationContext, extensions,
+        ckfOptions.propagatorPlainOptions, ckfOptions.referenceSurface,
+        ckfOptions.multipleScattering, ckfOptions.energyLoss, false);
+
+    return m_kalmanFitter.fit(it, end, sParameters, kfOptions, trackContainer);
+  }
+
+  template <typename sli_kf_t, typename sli_ckf_t, typename start_parameters_t,
+            typename track_container_t, template <typename> class holder_t>
+  auto findTracks(
+      sli_kf_t it, sli_kf_t end, const start_parameters_t& kfStartParameters,
+      const CombinatorialKalmanFilterOptions<sli_ckf_t, traj_t>& ckfOptions,
+      TrackContainer<track_container_t, traj_t, holder_t>& trackContainer) const
+      -> Result<CombinatorialKalmanFilterResult<traj_t>> {
+    // The KF run
+    auto kfResult =
+        runKalmanFitter(it, end, kfStartParameters, ckfOptions, trackContainer);
+
+    if (!kfResult.ok()) {
+      return kfResult.error();
+    }
+
+    ACTS_INFO("Done KF fitting");
+
+    // The CKF run
+    using SourceLinkAccessor = SourceLinkAccessorDelegate<sli_ckf_t>;
+
+    using ThisCkf = CombinatorialKalmanFilter<propagator_t, traj_t>;
+    using Aborter = typename ThisCkf::template Aborter<SourceLinkAccessor,
+                                                       BoundTrackParameters>;
+    using Actor = typename ThisCkf::template Actor<SourceLinkAccessor,
+                                                   BoundTrackParameters>;
+    using Actors = ActionList<Actor>;
+    using Aborters = AbortList<Aborter>;
+
+    PropagatorOptions<Actors, Aborters> propOptions(ckfOptions.geoContext,
+                                                    ckfOptions.magFieldContext);
+
+    propOptions.setPlainOptions(ckfOptions.propagatorPlainOptions);
+
+    auto& combKalmanActor = propOptions.actionList.template get<Actor>();
+    combKalmanActor.targetSurface = ckfOptions.referenceSurface;
+    combKalmanActor.multipleScattering = ckfOptions.multipleScattering;
+    combKalmanActor.energyLoss = ckfOptions.energyLoss;
+    combKalmanActor.smoothing = ckfOptions.smoothing;
+    combKalmanActor.m_sourcelinkAccessor = ckfOptions.sourcelinkAccessor;
+    combKalmanActor.m_extensions = ckfOptions.extensions;
+
+    // Prepare the start parameters
+    const auto state = *kfResult->trackStates().end();
+
+    BoundTrackParameters ckfStartParameters(
+        state.referenceSurface().getSharedPtr(), state.filtered(),
+        state.filteredCovariance());
+
+    // Prepare the result
+    using CkfResultType =
+        typename propagator_t::template action_list_t_result_t<
+            CurvilinearTrackParameters, Actors>;
+    CkfResultType inputResult;
+
+    auto& ckfResult =
+        inputResult.template get<CombinatorialKalmanFilterResult<traj_t>>();
+    ckfResult.fittedStates = &kfResult->container().trackStateContainer();
+    ckfResult.lastMeasurementIndices.push_back(kfResult->tipIndex());
+    ckfResult.lastTrackIndices.push_back(kfResult->tipIndex());
+
+    const CombinatorialKalmanFilterTipState tipState{
+        kfResult->nMeasurements() + kfResult->nHoles(),
+        kfResult->nTrackStates(), kfResult->nMeasurements(), 0ul,
+        kfResult->nHoles()};
+
+    ckfResult.activeTips.push_back({kfResult->tipIndex(), tipState});
+
+    // Run the CombinatorialKalmanFilter.
+    auto result = m_propagator.template propagate(
+        ckfStartParameters, propOptions, std::move(inputResult));
+
+    if (!result.ok()) {
+      ACTS_ERROR("Propapation failed: " << result.error() << " "
+                                        << result.error().message()
+                                        << " with the initial parameters:\n"
+                                        << ckfStartParameters.parameters());
+      return result.error();
+    }
+
+    ACTS_INFO("Done CKF fitting");
+
+    return result->template get<CombinatorialKalmanFilterResult<traj_t>>();
+  }
+};  // namespace Acts
+
+}  // namespace Acts