Draft: Start to migrate the R4-MS segment fit

Core/include/Acts/EventData/StationSpacePoint.hpp

@@ -0,0 +1,64 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2025 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 https://mozilla.org/MPL/2.0/.
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Utilities/PointerConcept.hpp"
+
+#include <type_traits>
+
+
+namespace Acts{
+    /** @brief Concept definition of the station space points. They're primarly used in composite detectors,
+     *         like the Muon stations in side the ATLAS experiment. The stations usually consist of few layers
+     *         of straw tubes which maybe sandwiched by other strip detector layers. The straws are used to measure
+     *         the passage of the particle in the bending plane, while the strip may supplement the track measurement
+     *         by providing measurements along the straw.
+     *  
+     */
+    template <typename SpacePointType>
+        concept StationSpacePoint = requires(const SpacePointType sp) {
+            /** @brief Local position of the space point measurement. It'either
+             *         the position of the wire or the position of the fired strip in the station */
+            { sp.localPosition() } -> std::same_as<const Acts::Vector3&>;
+            /** @brief Orientation of the sensor, which is either the wire orientation or 
+             *         the strip orientation. Travelling along the direction does not alter the residual */
+            { sp.sensorDirection()} -> std::same_as<const Acts::Vector3&>;
+            /** @brief Normal vector on the strip-plane. */
+            { sp.stripPlaneNormal()} -> std::same_as<const Acts::Vector3&>;
+            /** @brief Radius of the  tube drift-measurement. The returned value is zero for strip measurements */
+            { sp.driftRadius() } -> std::same_as<double>;
+            /** @brief Time when the measurement was taken */
+            { sp.time()} -> std::same_as<double>;
+            /** @brief Return the space point covariance. The upper left 2x2 block
+             *         describes the spatial meaurement uncertainty. The remaining block
+             *         the correlation between the time <-> spatial measurement or the pure time resolution */
+            { sp.covariance() } -> std::same_as<const ActsSquareMatrix<3>&>;
+        };
+
+
+        /** @brief Define the Space Point pointer concept as an ordinary / smart pointer 
+         *         over space points */
+        template <typename SpacePoint_t>
+             concept StationSpacePointPtr = PointerConcept<SpacePoint_t> &&
+                     StationSpacePoint<typename std::remove_pointer<SpacePoint_t>::type>;
+
+        /** @brief A station space point container is any std::container over space points */
+        template <typename ContType_t>
+            concept StationSpacePointContainer = requires(ContType_t cont, const ContType_t const_cont){
+                { cont.begin() } -> std::same_as<typename ContType_t::iterator>;
+                { cont.end() } -> std::same_as<typename ContType_t::iterator>;
+                { const_cont.begin() } -> std::same_as<typename ContType_t::const_iterator>;
+                { const_cont.end() } -> std::same_as<typename ContType_t::const_iterator>;                
+                { const_cont.size() } -> std::same_as<typename ContType_t::size_type>;
+                { const_cont.empty() } -> std::same_as<bool>;
+                requires StationSpacePointPtr<typename ContType_t::value_type>;
+            };
+
+
+    }

Core/include/Acts/Surfaces/detail/LineHelper.hpp

@@ -0,0 +1,104 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2025 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 https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Utilities/Intersection.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
+
+namespace Acts{
+    namespace LineHelper{
+        /** @brief Intersect two straight N-dimensional lines with each other or more generally calculate the point of 
+         *         closest approach of the second line to the first line.
+         *  @param linePosA: Arbitrary point on the first line
+         *  @param lineDirA: Direction of the first line (Unit-length)
+         *  @param linePosB: Arbitrary point on the second line
+         *  @param lineDirB: Direction of the second line (Unit-length) */
+        template <unsigned N>
+        inline Intersection<N> lineIntersect(const ActsVector<N>& linePosA, const ActsVector<N>& lineDirA,
+                                             const ActsVector<N>& linePosB, const ActsVector<N>& lineDirB) {
+
+            static_assert(N>=2, "One dimensional intersect not sensible");
+            /**  Use the formula
+             **    A + lambda dirA  = B + mu dirB
+             **    (A-B) + lambda dirA = mu dirB
+             **    <A-B, dirB> + lambda <dirA,dirB> = mu
+             **     A + lambda dirA = B + (<A-B, dirB> + lambda <dirA,dirB>)dirB
+             **     <A-B,dirA> + lambda <dirA, dirA> = <A-B, dirB><dirA,dirB> + lamda<dirA,dirB><dirA,dirB>
+             **   -> lambda = -(<A-B, dirA> - <A-B, dirB> * <dirA, dirB>) / (1- <dirA,dirB>^2)
+             **   --> mu    =  (<A-B, dirB> - <A-B, dirA> * <dirA, dirB>) / (1- <dirA,dirB>^2) */
+            const double dirDots = lineDirA.dot(lineDirB);
+            const double divisor = (1. - square(dirDots));
+            /// If the two directions are parallel to each other there's no way of intersection
+            if (std::abs(divisor) < std::numeric_limits<double>::epsilon()) {
+                return Intersection<N>::invalid();
+            }
+            const ActsVector<N> AminusB = linePosA - linePosB;
+            const double pathLength =  (AminusB.dot(lineDirB) - AminusB.dot(lineDirA) * dirDots) / divisor;
+
+            return Intersection<N>{linePosB + pathLength * lineDirB, pathLength, IntersectionStatus::onSurface};
+        }
+        /** @brief Intersect the lines of two line surfaces using their respective transforms.
+         *  @param lineSurfTrf1: local -> global transform of the first surface
+         *  @param lineSurfTrf2: local -> global transform of the second surface */
+        inline Intersection3D lineIntersect(const Acts::Transform3& lineSurfTrf1, 
+                                           const Acts::Transform3& lineSurfTrf2) {
+            return lineIntersect<3>(lineSurfTrf1.translation(), lineSurfTrf1.linear().col(2),
+                                    lineSurfTrf2.translation(), lineSurfTrf2.linear().col(2));
+        }
+        /** @brief Intersect a line in 3D space with a plane represented by the Hesse-Normal form
+         *  @param linePos: Arbitrary point on the line to intersect
+         *  @param lineDir: Direction of the line to intersect
+         *  @param planeNorm: Normal vector of the plane
+         *  @param offSet: Offset to move the plane along the normal vector */
+        inline Intersection3D planeIntersect(const Vector3& linePos,
+                                             const Vector3& lineDir,
+                                             const Vector3& planeNorm,
+                                             const double offset) {
+            /** Use the formula: <P, N> - C = 0
+             *   --> insert line equation: <A + lambda * B, N> - C = 0
+             *   --> lambda = (C - <A,N>)/ <N, B> */
+            const double normDot = planeNorm.dot(lineDir); 
+            if (std::abs(normDot) < std::numeric_limits<double>::epsilon()) {
+                return Intersection3D::invalid();
+            }
+            const double path = (offset - linePos.dot(planeNorm)) / normDot;
+            return Intersection3D{linePos + path * lineDir, path, IntersectionStatus::onSurface};
+        }
+        /** @brief Intersect a line in 3D space with a plane represented by the Hesse-Normal form
+         *  @param linePos: Arbitrary point on the line to intersect
+         *  @param lineDir: Direction of the line to intersect
+         *  @param planeNorm: Normal vector of the plane
+         *  @param planePoint: Point on the plane */
+        inline Intersection3D planeIntersect(const Vector3& linePos,
+                                             const Vector3& lineDir,
+                                             const Vector3& planeNorm,
+                                             const Vector3& planePoint) {
+            return planeIntersect(linePos, lineDir, planeNorm, planePoint.dot(planeNorm));
+        }
+        /** @brief Calculates the signed distance between two lines in 3D space 
+            *  @param linePosA: Arbitrary point on the first line
+         *  @param lineDirA: Direction of the first line (Unit-length)
+         *  @param linePosB: Arbitrary point on the second line
+         *  @param lineDirB: Direction of the second line (Unit-length) */
+        inline double signedDistance(const Vector3& linePosA,
+                                     const Vector3& lineDirA,
+                                     const Vector3& linePosB,
+                                     const Vector3& lineDirB) {
+            /** Project the first direction onto the second & renormalize to a unit vector */
+            const double dirDots = lineDirA.dot(lineDirB);
+            const Vector3 AminusB = linePosA - linePosB;
+            if (std::abs(dirDots -1.) <  std::numeric_limits<double>::epsilon()){
+            return (AminusB - lineDirA.dot(AminusB)*lineDirA).norm();
+            }
+            const Vector3 projDir = (lineDirA - dirDots*lineDirB).normalized();
+            return AminusB.cross(lineDirB).dot(projDir);
+        }
+    }
+}

Core/include/Acts/Utilities/MathHelpers.hpp

@@ -27,4 +27,14 @@ constexpr auto fastHypot(T... args) {
   return std::sqrt(hypotSquare(args...));
 }
 
+/// @brief Helper struct to calculate sin/cos in one go
+template <typename T>
+struct sincos{
+  sincos(const T alpha):
+      cs{std::cos(alpha)},
+      sn{std::sin(alpha)}{}
+  T cs{0.};
+  T sn{0.};
+};
+
 }  // namespace Acts

Core/include/Acts/Utilities/PointerConcept.hpp

@@ -0,0 +1,41 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2025 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 https://mozilla.org/MPL/2.0/.
+#pragma once
+#include <concepts>
+#include <type_traits>
+namespace Acts {
+    /** @brief The Pointer concept is an extension of the usual std::is_pointer_v type trait to 
+     *         also include the smart pointers like 
+     *     std::shared_ptr<T>,
+     *     std::unique_ptr<T>
+     *  The smart pointer is required to have an element_type typedef indicating over which data
+     *  type the pointer is constructed, the arrow operator
+     * 
+     *      T* operator->() const;
+     *  and also the dereference operator
+     * 
+     *      T& operator*() const  */
+    template <typename Pointer_t>
+        concept SmartPointerConcept = requires(Pointer_t ptr) {
+        /** @brief arrow operator element_type* operator->() const; */
+        {ptr.operator->()} -> std::same_as<std::add_pointer_t<typename Pointer_t::element_type>>;
+        /** @brief dereference operator element_type& operator->() const; */
+        {ptr.operator*()} -> std::same_as<typename Pointer_t::element_type&>;
+
+    }; 
+    template <typename Pointer_t> 
+        concept PointerConcept = (std::is_pointer_v<Pointer_t> || SmartPointerConcept<Pointer_t>); 
+}
+
+
+namespace std{
+template <Acts::SmartPointerConcept T>
+    struct remove_pointer<T>{
+        typedef typename T::element_type type;
+};
+}

Core/include/Acts/Utilities/UnitVectors.hpp

@@ -9,6 +9,7 @@
 #pragma once
 
 #include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Utilities/MathHelpers.hpp"
 
 #include <cmath>
 #include <limits>
@@ -46,23 +47,33 @@ inline Eigen::Matrix<T, 3, 1> makeDirectionFromPhiEta(T phi, T eta) {
 ///       explicitly cast mismatched input types.
 template <typename T>
 inline Eigen::Matrix<T, 3, 1> makeDirectionFromPhiTheta(T phi, T theta) {
-  const auto cosTheta = std::cos(theta);
-  const auto sinTheta = std::sin(theta);
+  sincos csTheta{theta};
+  sincos csPhi{phi};
   return {
-      std::cos(phi) * sinTheta,
-      std::sin(phi) * sinTheta,
-      cosTheta,
+      csPhi.cs * csTheta.sn,
+      csPhi.sn * csTheta.sn,
+      csTheta.cs,
   };
 }
+/// @brief Construct a normalized direction vector from the tangents of the
+///        x-axis to the z-axis and of the y-axis to the z-axis           
+///
+/// @param tanAlpha: Tangent of the x-axis to the z-axis
+/// @param tanBeta: Tangent of the y-axis to the z-axis 
+template <typename T>
+  inline Eigen::Matrix<T,3, 1> makeDirectionFromAxisTangents(T tanalpha, T tanBeta) {
+      return Eigen::Matrix<T,3, 1>{tanalpha, tanBeta, 1}.normalized();
+  }
+
+
 
 /// Construct a phi and theta angle from a direction vector.
 ///
 /// @param unitDir 3D vector indicating a direction
 ///
 template <typename T>
 inline Eigen::Matrix<T, 2, 1> makePhiThetaFromDirection(
-    Eigen::Matrix<T, 3, 1> unitDir) {
-  unitDir.normalize();
+    const Eigen::Matrix<T, 3, 1>& unitDir) {
   T phi = std::atan2(unitDir[1], unitDir[0]);
   T theta = std::acos(unitDir[2]);
   return {

Examples/Algorithms/TrackFinding/CMakeLists.txt

@@ -9,6 +9,7 @@ add_library(
     src/HoughTransformSeeder.cpp
     src/TrackParamsEstimationAlgorithm.cpp
     src/MuonHoughSeeder.cpp
+    src/MuonSegmentFinder.cpp
     src/GbtsSeedingAlgorithm.cpp
     src/TrackParamsLookupEstimation.cpp
 )

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

@@ -8,31 +8,31 @@
 
 #pragma once
 
-#include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Seeding/HoughTransformUtils.hpp"
 #include "Acts/Utilities/Delegate.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
-#include "ActsExamples/EventData/DriftCircle.hpp"
-#include "ActsExamples/EventData/SimHit.hpp"
+
+#include "ActsExamples/EventData/MuonSegment.hpp"
+#include "ActsExamples/EventData/MuonSpacePoint.hpp"
+#include "ActsExamples/EventData/MuonHoughMaximum.hpp"
+
 #include "ActsExamples/Framework/DataHandle.hpp"
 #include "ActsExamples/Framework/IAlgorithm.hpp"
 #include "ActsExamples/Framework/ProcessCode.hpp"
+#include "Acts/Definitions/Units.hpp"
 
 #include <cstddef>
 #include <memory>
 #include <string>
-#include <unordered_set>
 #include <utility>
 #include <vector>
 
 #include "TCanvas.h"
-#include "TH2D.h"
-#include "TMarker.h"
-#include "TStyle.h"
+
 
 namespace ActsExamples {
-struct AlgorithmContext;
+  struct AlgorithmContext;
 }
 
 namespace ActsExamples {
@@ -48,8 +48,15 @@ class MuonHoughSeeder final : public IAlgorithm {
  public:
   /// config
   struct Config {
-    std::string inSimHits;
-    std::string inDriftCircles;
+    std::string inTruthSegments{};
+    std::string inSpacePoints{};
+    std::string outHoughMax{};
+
+    /** @brief Extra margin added to both y-sides of the eta-hough accumulator plane */
+    double etaPlaneMarginIcept{10.*Acts::UnitConstants::cm};
+    /** @brief Extra margin added to both y-sides of the phi-hough accumulator plane */
+    double phiPlaneMarginIcept{10.*Acts::UnitConstants::cm};
+
   };
 
   MuonHoughSeeder(Config cfg, Acts::Logging::Level lvl);
@@ -70,9 +77,9 @@ class MuonHoughSeeder final : public IAlgorithm {
   std::unique_ptr<const Acts::Logger> m_logger;
   const Acts::Logger& logger() const { return *m_logger; }
 
-  ReadDataHandle<SimHitContainer> m_inputSimHits{this, "InputSimHits"};
-  ReadDataHandle<DriftCircleContainer> m_inputDriftCircles{this,
-                                                           "InputDriftCircles"};
+  ReadDataHandle<MuonSegmentContainer> m_inputTruthSegs{this, "InputTruthSegments"};
+  ReadDataHandle<MuonSpacePointContainer> m_inputSpacePoints{this, "InputSpacePoints"};
+  WriteDataHandle<MuonHoughMaxContainer> m_outputMaxima{this, "OutputHoughMax"};
   /// use ROOT for visualisation
   std::unique_ptr<TCanvas> m_outCanvas;
 };