fix: avoid using ActsExamples::Trajectories in IterativeVertexFinder

src/algorithms/tracking/IterativeVertexFinder.cc

@@ -4,11 +4,19 @@
 
 #include "IterativeVertexFinder.h"
 
+#include <Acts/Definitions/Algebra.hpp>
+#include <Acts/Definitions/TrackParametrization.hpp>
 #include <Acts/Definitions/Units.hpp>
+#include <Acts/EventData/GenericBoundTrackParameters.hpp>
+#include <Acts/EventData/ParticleHypothesis.hpp>
+#include <Acts/EventData/TrackParameters.hpp>
+#include <Acts/Geometry/TrackingGeometry.hpp>
 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
 #include <Acts/Propagator/EigenStepper.hpp>
 #include <Acts/Propagator/Propagator.hpp>
 #include <Acts/Propagator/VoidNavigator.hpp>
+#include <Acts/Surfaces/PerigeeSurface.hpp>
+#include <Acts/Surfaces/Surface.hpp>
 #include <Acts/Utilities/Logger.hpp>
 #include <Acts/Utilities/Result.hpp>
 #include <Acts/Utilities/detail/ContextType.hpp>
@@ -22,26 +30,43 @@
 #include <Acts/Vertexing/Vertex.hpp>
 #include <Acts/Vertexing/VertexingOptions.hpp>
 #include <Acts/Vertexing/ZScanVertexFinder.hpp>
-#include <ActsExamples/EventData/Track.hpp>
 #include <ActsExamples/EventData/Trajectories.hpp>
+#include <boost/container/detail/std_fwd.hpp>
 #include <edm4eic/Cov4f.h>
+#include <edm4eic/Cov6f.h>
 #include <edm4eic/ReconstructedParticleCollection.h>
 #include <edm4eic/Track.h>
 #include <edm4eic/TrackParameters.h>
 #include <edm4eic/Trajectory.h>
-#include <edm4eic/unit_system.h>
 #include <edm4hep/Vector2f.h>
 #include <edm4hep/Vector4f.h>
 #include <fmt/core.h>
 #include <fmt/format.h>
 #include <podio/RelationRange.h>
 #include <Eigen/Core>
-#include <cmath>
+#include <array>
+#include <cstddef>
+#include <iterator>
+#include <ranges>
 #include <string>
 #include <utility>
 
 #include "extensions/spdlog/SpdlogToActs.h"
 
+// This array relates the Acts and EDM4eic covariance matrices, including
+// the unit conversion to get from Acts units into EDM4eic units.
+//
+// Note: std::map is not constexpr, so we use a constexpr std::array
+// std::array initialization need double braces since arrays are aggregates
+// ref: https://en.cppreference.com/w/cpp/language/aggregate_initialization
+static constexpr std::array<std::pair<Acts::BoundIndices, double>, 6> edm4eic_indexed_units{
+    {{Acts::eBoundLoc0, Acts::UnitConstants::mm},
+     {Acts::eBoundLoc1, Acts::UnitConstants::mm},
+     {Acts::eBoundPhi, 1.},
+     {Acts::eBoundTheta, 1.},
+     {Acts::eBoundQOverP, 1. / Acts::UnitConstants::GeV},
+     {Acts::eBoundTime, Acts::UnitConstants::ns}}};
+
 void eicrecon::IterativeVertexFinder::init(std::shared_ptr<const ActsGeometryProvider> geo_svc,
                                            std::shared_ptr<spdlog::logger> log) {
 
@@ -55,7 +80,7 @@ void eicrecon::IterativeVertexFinder::init(std::shared_ptr<const ActsGeometryPro
 }
 
 std::unique_ptr<edm4eic::VertexCollection> eicrecon::IterativeVertexFinder::produce(
-    std::vector<const ActsExamples::Trajectories*> trajectories,
+    std::vector<const ActsExamples::Trajectories*> /* trajectories */,
     const edm4eic::ReconstructedParticleCollection* reconParticles) {
 
   auto outputVertices = std::make_unique<edm4eic::VertexCollection>();
@@ -112,20 +137,58 @@ std::unique_ptr<edm4eic::VertexCollection> eicrecon::IterativeVertexFinder::prod
   VertexFinderOptions finderOpts(m_geoctx, m_fieldctx);
 
   std::vector<Acts::InputTrack> inputTracks;
-
-  for (const auto& trajectory : trajectories) {
-    auto tips = trajectory->tips();
-    if (tips.empty()) {
-      continue;
-    }
-    /// CKF can provide multiple track trajectories for a single input seed
-    for (auto& tip : tips) {
-      ActsExamples::TrackParameters par = trajectory->trackParameters(tip);
-
-      inputTracks.emplace_back(&(trajectory->trackParameters(tip)));
-      m_log->trace("Track local position at input = {} mm, {} mm",
-                   par.localPosition().x() / Acts::UnitConstants::mm,
-                   par.localPosition().y() / Acts::UnitConstants::mm);
+  std::vector<Acts::BoundTrackParameters> inputTrackParameters;
+  std::vector<const edm4eic::ReconstructedParticle*> inputParticles;
+
+  for (const auto& particle : *reconParticles) {
+    for (const auto& track : particle.getTracks()) {
+      const auto& trajectory = track.getTrajectory();
+      for (const auto& track_parameter : trajectory.getTrackParameters()) {
+
+        std::shared_ptr<const Acts::Surface> surface;
+        // Get reference surface by geometryId
+        if (const auto* surface_ptr =
+                m_geoSvc->trackingGeometry()->findSurface(track_parameter.getSurface())) {
+          surface = surface_ptr->getSharedPtr();
+        } else {
+          // Perigee surface is the default
+          surface = std::dynamic_pointer_cast<const Acts::Surface>(
+              Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.}));
+        }
+
+        // Parameters
+        Acts::BoundVector params;
+        params(Acts::eBoundLoc0) =
+            track_parameter.getLoc().a * Acts::UnitConstants::mm; // cylinder radius
+        params(Acts::eBoundLoc1) =
+            track_parameter.getLoc().b * Acts::UnitConstants::mm; // cylinder length
+        params(Acts::eBoundPhi)    = track_parameter.getPhi();
+        params(Acts::eBoundTheta)  = track_parameter.getTheta();
+        params(Acts::eBoundQOverP) = track_parameter.getQOverP() / Acts::UnitConstants::GeV;
+        params(Acts::eBoundTime)   = track_parameter.getTime() * Acts::UnitConstants::ns;
+
+        // Covariance FIXME stick edm4eic_indexed_units elsewhere
+        Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero();
+        for (std::size_t i = 0; const auto& [a, x] : edm4eic_indexed_units) {
+          for (std::size_t j = 0; const auto& [b, y] : edm4eic_indexed_units) {
+            cov(a, b) = track_parameter.getCovariance()(i, j) * x * y;
+            ++j;
+          }
+          ++i;
+        }
+
+        // Finally create BoundTrackParameters
+        auto& par = inputTrackParameters.emplace_back(surface, params, cov,
+                                                      Acts::ParticleHypothesis::pion());
+
+        inputTracks.emplace_back(&par);
+        m_log->trace("Track local position at input = {} mm, {} mm",
+                     par.localPosition().x() / Acts::UnitConstants::mm,
+                     par.localPosition().y() / Acts::UnitConstants::mm);
+
+        // Store reference to particles
+        inputParticles.push_back(&particle);
+      }
     }
   }
 
@@ -158,29 +221,14 @@ std::unique_ptr<edm4eic::VertexCollection> eicrecon::IterativeVertexFinder::prod
       m_log->trace("Track local position from vertex = {} mm, {} mm",
                    par.localPosition().x() / Acts::UnitConstants::mm,
                    par.localPosition().y() / Acts::UnitConstants::mm);
-      float loc_a = par.localPosition().x();
-      float loc_b = par.localPosition().y();
-
-      for (const auto& part : *reconParticles) {
-        const auto& tracks = part.getTracks();
-        for (const auto& trk : tracks) {
-          const auto& traj    = trk.getTrajectory();
-          const auto& trkPars = traj.getTrackParameters();
-          for (const auto& par : trkPars) {
-            const double EPSILON = 1.0e-4; // mm
-            if (std::abs((par.getLoc().a / edm4eic::unit::mm) - (loc_a / Acts::UnitConstants::mm)) <
-                    EPSILON &&
-                std::abs((par.getLoc().b / edm4eic::unit::mm) - (loc_b / Acts::UnitConstants::mm)) <
-                    EPSILON) {
-              m_log->trace(
-                  "From ReconParticles, track local position [Loc a, Loc b] = {} mm, {} mm",
-                  par.getLoc().a / edm4eic::unit::mm, par.getLoc().b / edm4eic::unit::mm);
-              eicvertex.addToAssociatedParticles(part);
-            } // endif
-          } // end for par
-        } // end for trk
-      } // end for part
+
+      auto inputTrackIter = std::ranges::find(inputTracks, t.originalParams);
+      auto i              = std::distance(inputTracks.begin(), inputTrackIter);
+
+      eicvertex.addToAssociatedParticles(*inputParticles.at(i));
+
     } // end for t
+
     m_log->debug("One vertex found at (x,y,z) = ({}, {}, {}) mm.",
                  vtx.position().x() / Acts::UnitConstants::mm,
                  vtx.position().y() / Acts::UnitConstants::mm,