Track extrapolation to ECAL calorimeter surface

k4ActsTracking/include/k4ActsTracking/ACTSSeededCKFTrackingAlg.hxx

@@ -106,8 +106,9 @@ public:
                                                             Acts::MagneticFieldProvider::Cache& magCache) const;
 
   StatusCode tracking(const std::vector<Acts::BoundTrackParameters>& paramseeds, const CKF& trackFinder,
-                      const TrackFinderOptions& ckfOptions, Acts::MagneticFieldProvider::Cache& magCache,
-                      edm4hep::TrackCollection& trackCollection) const;
+                      const TrackFinderOptions& ckfOptions, const Propagator& extrapPropagator,
+                      const Acts::PerigeeSurface& perigeeSurface, Propagator::Options<>& extrapOptions,
+                      Acts::MagneticFieldProvider::Cache& magCache, edm4hep::TrackCollection& trackCollection) const;
 
 protected:
   /**

k4ActsTracking/include/k4ActsTracking/Helpers.hxx

@@ -40,15 +40,30 @@
 #include <Acts/EventData/ParticleHypothesis.hpp>
 #include <Acts/EventData/VectorMultiTrajectory.hpp>
 #include <Acts/EventData/VectorTrackContainer.hpp>
+#include <Acts/Geometry/GeometryContext.hpp>
+#include <Acts/MagneticField/MagneticFieldContext.hpp>
 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
+#include <Acts/Propagator/EigenStepper.hpp>
+#include <Acts/Propagator/Propagator.hpp>
+#include <Acts/Propagator/VoidNavigator.hpp>
 #include <Acts/TrackFinding/CombinatorialKalmanFilter.hpp>
 #include <Acts/TrackFitting/KalmanFitter.hpp>
 
 // ACTSTracking
+#include "k4ActsTracking/IActsGeoSvc.h"
 #include "k4ActsTracking/SourceLink.hxx"
 
+// Standard
+#include <optional>
+
 namespace ACTSTracking {
 
+  /// Propagator used to extrapolate fitted tracks out to the calorimeter face.
+  /// It uses a VoidNavigator (no tracking geometry), so it can reach target
+  /// surfaces that lie outside the tracking-geometry world volume. Material
+  /// between the tracker and the calorimeter is not accounted for (field-only).
+  using CaloFacePropagator = Acts::Propagator<Acts::EigenStepper<>, Acts::VoidNavigator>;
+
   using TrackResult =
       Acts::TrackContainer<Acts::VectorTrackContainer, Acts::VectorMultiTrajectory, std::shared_ptr>::TrackProxy;
 
@@ -106,4 +121,28 @@ namespace ACTSTracking {
  */
   Acts::ParticleHypothesis convertParticle(const edm4hep::MCParticle mcParticle);
 
+  //! Extrapolate track parameters to the calorimeter face.
+  /**
+ * Selects, among the calorimeter-face surfaces, the one the track reaches first
+ * (smallest positive path length with a valid, bounds-checked intersection) and
+ * propagates the given parameters to it. The barrel is a ring of planar
+ * surfaces (one per polygon side) and the endcaps are flat discs; intersecting
+ * all candidates naturally handles both the barrel-sector choice and the
+ * barrel/endcap transition.
+ *
+ * \param propagator Field-only propagator (VoidNavigator)
+ * \param start      Track parameters to extrapolate from (e.g. at the last hit)
+ * \param surfaces   Calorimeter-face surfaces from IActsGeoSvc
+ * \param gctx       Geometry context
+ * \param mctx       Magnetic-field context
+ *
+ * \return Bound track parameters at the calorimeter face, or std::nullopt if no
+ *         surface is reached or the propagation fails.
+ */
+  std::optional<Acts::BoundTrackParameters> extrapolateToCaloFace(const CaloFacePropagator&            propagator,
+                                                                  const Acts::BoundTrackParameters&    start,
+                                                                  const IActsGeoSvc::CaloFaceSurfaces& surfaces,
+                                                                  const Acts::GeometryContext&         gctx,
+                                                                  const Acts::MagneticFieldContext&    mctx);
+
 }  // namespace ACTSTracking

k4ActsTracking/include/k4ActsTracking/IActsGeoSvc.h

@@ -27,6 +27,7 @@
 #include <memory>
 #include <string>
 #include <unordered_map>
+#include <vector>
 
 namespace dd4hep {
   namespace rec {
@@ -44,13 +45,30 @@ class GAUDI_API IActsGeoSvc : virtual public IService {
 public:
   using CellIDSurfaceMap = std::unordered_map<dd4hep::CellID, const Acts::Surface*>;
 
+  /// Surfaces approximating the inner face of the electromagnetic calorimeter,
+  /// derived from the DD4hep geometry. These live outside the tracking-geometry
+  /// world volume and are meant as target surfaces for a geometry-free
+  /// extrapolation (see ACTSTracking::extrapolateToCaloFace).
+  ///
+  /// The barrel is a regular polygon, so its inner face is modelled as one
+  /// planar surface per polygon side rather than a cylinder. The endcaps are
+  /// flat discs at constant z.
+  struct CaloFaceSurfaces {
+    std::vector<std::shared_ptr<const Acts::Surface>> barrelFaces;  ///< one plane per polygon side
+    std::shared_ptr<const Acts::Surface>              endcapPos;    ///< disc at +z, may be null
+    std::shared_ptr<const Acts::Surface>              endcapNeg;    ///< disc at -z, may be null
+
+    bool empty() const { return barrelFaces.empty() && !endcapPos && !endcapNeg; }
+  };
+
 public:
   DeclareInterfaceID(IActsGeoSvc, 1, 0);
 
   virtual std::shared_ptr<const Acts::TrackingGeometry>      trackingGeometry() const     = 0;
   virtual std::shared_ptr<const Acts::MagneticFieldProvider> magneticField() const        = 0;
   virtual const CellIDSurfaceMap&                            cellIdToSurfaceMap() const   = 0;
   virtual std::string                                        cellIDEncodingString() const = 0;
+  virtual const CaloFaceSurfaces&                            caloFaceSurfaces() const     = 0;
 
   virtual ~IActsGeoSvc() = default;
 };

k4ActsTracking/src/components/ACTSSeededCKFTrackingAlg.cxx

@@ -43,6 +43,7 @@
 #include <Acts/TrackFitting/GainMatrixUpdater.hpp>
 #include <Acts/Utilities/Logger.hpp>
 #include <Acts/Utilities/RangeXD.hpp>
+#include <Acts/Utilities/TrackHelpers.hpp>
 
 // TBB
 #include <tbb/concurrent_vector.h>
@@ -384,6 +385,12 @@ std::tuple<edm4hep::TrackCollection, edm4hep::TrackCollection> ACTSSeededCKFTrac
   Propagator propagator(std::move(stepper), std::move(navigator));
   CKF        trackFinder(std::move(propagator));
 
+  // For extrapolating the fitted track back to the IP (perigee surface)
+  Stepper    extrapStepper(magneticField());
+  Navigator  extrapNavigator(navigatorCfg);
+  Propagator extrapPropagator(std::move(extrapStepper), std::move(extrapNavigator));
+  auto       perigeeSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
+
   // Set the options
   Acts::MeasurementSelector::Config measurementSelectorCfg = {
       {Acts::GeometryIdentifier(), {{}, {m_CKF_chi2CutOff}, {(std::size_t)(m_CKF_numMeasurementsCutOff)}}}};
@@ -394,9 +401,11 @@ std::tuple<edm4hep::TrackCollection, edm4hep::TrackCollection> ACTSSeededCKFTrac
     pOptions.direction = Acts::Direction::Backward();
   }
 
-  // Construct a perigee surface as the target surface
-  std::shared_ptr<Acts::PerigeeSurface> perigeeSurface =
-      Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
+  Propagator::Options<> extrapOptions{geometryContext(), magneticFieldContext()};
+  extrapOptions.maxSteps = 10000;
+  if (m_propagateBackward) {
+    extrapOptions.direction = Acts::Direction::Backward();
+  }
 
   Acts::GainMatrixUpdater kfUpdater;
 
@@ -471,7 +480,9 @@ std::tuple<edm4hep::TrackCollection, edm4hep::TrackCollection> ACTSSeededCKFTrac
     if (!m_runCKF)
       return;
 
-    if (!tracking(paramseeds, trackFinder, ckfOptions, localMagCache, trackCollection).isSuccess()) {
+    if (!tracking(paramseeds, trackFinder, ckfOptions, extrapPropagator, *perigeeSurface, extrapOptions, localMagCache,
+                  trackCollection)
+             .isSuccess()) {
       warning() << "Tracking failed for this event" << endmsg;
     }
   };  // parallelSeedingAndTracking
@@ -586,6 +597,9 @@ std::vector<Acts::BoundTrackParameters> ACTSSeededCKFTrackingAlg::seedsToParamet
 // CKF tracking,
 StatusCode ACTSSeededCKFTrackingAlg::tracking(const std::vector<Acts::BoundTrackParameters>& paramseeds,
                                               const CKF& trackFinder, const TrackFinderOptions& ckfOptions,
+                                              const Propagator&                   extrapPropagator,
+                                              const Acts::PerigeeSurface&         perigeeSurface,
+                                              Propagator::Options<>&              extrapOptions,
                                               Acts::MagneticFieldProvider::Cache& magCache,
                                               edm4hep::TrackCollection&           trackCollection) const {
   // Initialize track finder
@@ -613,6 +627,16 @@ StatusCode ACTSSeededCKFTrackingAlg::tracking(const std::vector<Acts::BoundTrack
           continue;
         }
 
+        // Extrapolate the fitted track back to the perigee surface at the IP so
+        // that the track parameters (in particular D0 and Z0) are expressed
+        // there. This reproduces the TrackStateAtIP behaviour from older ACTS.
+        auto extrapResult = Acts::extrapolateTrackToReferenceSurface(
+            trackTip, perigeeSurface, extrapPropagator, extrapOptions, Acts::TrackExtrapolationStrategy::firstOrLast);
+        if (!extrapResult.ok()) {
+          warning() << "Track extrapolation to perigee failed: " << extrapResult.error() << endmsg;
+          continue;
+        }
+
         // Helpful debug output
         debug() << "Trajectory Summary" << endmsg;
         debug() << "\tchi2Sum       " << trackTip.chi2() << endmsg;

k4ActsTracking/src/components/ActsGeoSvc.cpp

@@ -37,6 +37,9 @@
 #include <Acts/Geometry/TrackingGeometry.hpp>
 #include <Acts/Geometry/VolumeAttachmentStrategy.hpp>
 #include <Acts/MagneticField/ConstantBField.hpp>
+#include <Acts/Surfaces/DiscSurface.hpp>
+#include <Acts/Surfaces/PlaneSurface.hpp>
+#include <Acts/Surfaces/RectangleBounds.hpp>
 #include <Acts/Surfaces/Surface.hpp>
 #include <Acts/Utilities/AxisDefinitions.hpp>
 #include <Acts/Visualization/ObjVisualization3D.hpp>
@@ -45,7 +48,9 @@
 
 #include <DD4hep/DD4hepUnits.h>
 #include <DD4hep/DetElement.h>
+#include <DD4hep/DetType.h>
 #include <DD4hep/Detector.h>
+#include <DDRec/DetectorData.h>
 
 #include <GaudiKernel/StatusCode.h>
 
@@ -54,6 +59,8 @@
 #include <fmt/ranges.h>
 
 #include <array>
+#include <cmath>
+#include <numbers>
 
 template <> struct fmt::formatter<Acts::GeometryIdentifier> : fmt::ostream_formatter {};
 
@@ -159,5 +166,137 @@ StatusCode ActsGeoSvc::initialize() {
     vis.write(m_objDumpFileName.value());
   }
 
+  if (m_buildCaloSurfaces.value()) {
+    buildCaloFaceSurfaces();
+  }
+
   return StatusCode::SUCCESS;
 }
+
+void ActsGeoSvc::buildCaloFaceSurfaces() {
+  using dd4hep::DetType;
+  using dd4hep::rec::LayeredCalorimeterData;
+  namespace UC = Acts::UnitConstants;
+
+  // LayeredCalorimeterData::extent[] is stored in DD4hep native length units
+  // (cm), so convert to ACTS units with the same scale the blueprint uses.
+  const double lengthScale = Acts::UnitConstants::cm / dd4hep::cm;
+
+  const auto* dd4hepDet = m_geoSvc->getDetector();
+
+  // Locate the electromagnetic calorimeter sub-detectors purely via DetType
+  // flags, so that this works across detector concepts without hard-coding
+  // DetElement names. The dimensions themselves come from the standard DDRec
+  // LayeredCalorimeterData extension (the same source Pandora uses).
+  const auto ecalBarrel =
+      dd4hepDet->detectors(DetType::CALORIMETER | DetType::ELECTROMAGNETIC | DetType::BARREL, DetType::FORWARD);
+  const auto ecalEndcap = dd4hepDet->detectors(DetType::CALORIMETER | DetType::ELECTROMAGNETIC | DetType::ENDCAP, 0);
+
+  // Barrel circumradius (corner radius), needed both for the barrel faces and to
+  // guarantee the endcap discs reach far enough to avoid hermeticity gaps at the
+  // barrel/endcap junction.
+  double barrelCircumradius = 0.0;
+  double barrelHalfZ        = 0.0;
+
+  // --- Barrel: one planar surface per polygon side --------------------------
+  if (ecalBarrel.empty()) {
+    warning() << "No electromagnetic barrel calorimeter found via DetType flags; "
+                 "no barrel calo-face surfaces will be built."
+              << endmsg;
+  } else {
+    const auto* caloData = ecalBarrel.front().extension<LayeredCalorimeterData>(false);
+    if (caloData == nullptr) {
+      warning() << "ECAL barrel DetElement has no LayeredCalorimeterData extension; "
+                   "skipping barrel calo-face surfaces."
+                << endmsg;
+    } else {
+      const int nSides = caloData->inner_symmetry > 0 ? caloData->inner_symmetry : 0;
+      // For a BarrelLayout the calorimeter is centred on z = 0 and extent[] is
+      // {rmin, rmax, zmin=0, zmax=half_length}, so the half-length is extent[3].
+      const double apothem = caloData->extent[0] * lengthScale;  // perpendicular distance to inner face
+      barrelHalfZ          = caloData->extent[3] * lengthScale;  // barrel half-length
+      const double phi0    = caloData->inner_phi0;               // azimuth of first inner-face normal
+
+      if (nSides < 3) {
+        warning() << fmt::format("ECAL barrel has unusable inner_symmetry={}; skipping barrel faces.",
+                                 caloData->inner_symmetry)
+                  << endmsg;
+      } else {
+        const double dPhi      = 2 * std::numbers::pi / nSides;
+        barrelCircumradius     = apothem / std::cos(std::numbers::pi / nSides);
+        const double halfWidth = apothem * std::tan(std::numbers::pi / nSides);
+
+        m_caloFaceSurfaces.barrelFaces.reserve(nSides);
+        for (int i = 0; i < nSides; ++i) {
+          const double  phi  = phi0 + i * dPhi;
+          const double  cphi = std::cos(phi);
+          const double  sphi = std::sin(phi);
+          Acts::Vector3 normal{cphi, sphi, 0};      // local z (surface normal, radial)
+          Acts::Vector3 localX{-sphi, cphi, 0};     // tangential
+          Acts::Vector3 localY{0, 0, 1};            // along global z
+          Acts::Vector3 center = apothem * normal;  // barrel centred on z = 0
+
+          Acts::Transform3 transform = Acts::Transform3::Identity();
+          transform.linear().col(0)  = localX;
+          transform.linear().col(1)  = localY;
+          transform.linear().col(2)  = normal;
+          transform.translation()    = center;
+
+          auto bounds = std::make_shared<Acts::RectangleBounds>(halfWidth, barrelHalfZ);
+          m_caloFaceSurfaces.barrelFaces.push_back(Acts::Surface::makeShared<Acts::PlaneSurface>(transform, bounds));
+        }
+        info() << fmt::format(
+                      "Built ECAL barrel calo face: {} planar faces, apothem={:.1f} mm, circumradius={:.1f} mm, "
+                      "halfZ={:.1f} mm, phi0={:.4f}",
+                      nSides, apothem / UC::mm, barrelCircumradius / UC::mm, barrelHalfZ / UC::mm, phi0)
+               << endmsg;
+      }
+    }
+  }
+
+  // --- Endcaps: flat discs, widened to stay hermetic with the barrel --------
+  if (ecalEndcap.empty()) {
+    warning() << "No electromagnetic endcap calorimeter found via DetType flags; "
+                 "no endcap calo-face surfaces will be built."
+              << endmsg;
+  } else {
+    const auto* caloData = ecalEndcap.front().extension<LayeredCalorimeterData>(false);
+    if (caloData == nullptr) {
+      warning() << "ECAL endcap DetElement has no LayeredCalorimeterData extension; "
+                   "skipping endcap calo-face surfaces."
+                << endmsg;
+    } else {
+      const double rMin    = caloData->extent[0] * lengthScale;
+      double       rMax    = caloData->extent[1] * lengthScale;
+      const double zEndcap = caloData->extent[2] * lengthScale;  // inner-face z (zmin)
+
+      // Guarantee the disc reaches at least the barrel corner radius so there is
+      // no gap at the barrel/endcap junction.
+      if (barrelCircumradius > rMax) {
+        warning() << fmt::format(
+                         "ECAL endcap rMax ({:.1f} mm) is smaller than the barrel circumradius ({:.1f} mm); "
+                         "extending the endcap disc to close the hermeticity gap.",
+                         rMax / UC::mm, barrelCircumradius / UC::mm)
+                  << endmsg;
+        rMax = barrelCircumradius;
+      }
+
+      Acts::Transform3 tPos = Acts::Transform3::Identity();
+      tPos.translation()    = Acts::Vector3{0, 0, zEndcap};
+      Acts::Transform3 tNeg = Acts::Transform3::Identity();
+      tNeg.translation()    = Acts::Vector3{0, 0, -zEndcap};
+
+      m_caloFaceSurfaces.endcapPos = Acts::Surface::makeShared<Acts::DiscSurface>(tPos, rMin, rMax);
+      m_caloFaceSurfaces.endcapNeg = Acts::Surface::makeShared<Acts::DiscSurface>(tNeg, rMin, rMax);
+
+      info() << fmt::format("Built ECAL endcap calo faces: discs at z=+/-{:.1f} mm, rMin={:.1f} mm, rMax={:.1f} mm",
+                            zEndcap / UC::mm, rMin / UC::mm, rMax / UC::mm)
+             << endmsg;
+    }
+  }
+
+  info() << fmt::format("Calo-face surfaces built: {} barrel faces, {} endcap discs",
+                        m_caloFaceSurfaces.barrelFaces.size(),
+                        (m_caloFaceSurfaces.endcapPos ? 1 : 0) + (m_caloFaceSurfaces.endcapNeg ? 1 : 0))
+         << endmsg;
+}

k4ActsTracking/src/components/ActsGeoSvc.h

@@ -53,6 +53,8 @@ class ActsGeoSvc : public extends<Service, IActsGeoSvc> {
 
   std::shared_ptr<const Acts::MagneticFieldProvider> magneticField() const override;
 
+  const CaloFaceSurfaces& caloFaceSurfaces() const override { return m_caloFaceSurfaces; }
+
   ActsGeoSvc(const std::string& name, ISvcLocator* svcLoc);
 
   ~ActsGeoSvc() = default;
@@ -66,6 +68,9 @@ class ActsGeoSvc : public extends<Service, IActsGeoSvc> {
   Gaudi::Property<std::string> m_encodingStringConstant{
       this, "EncodingStringVariable", "GlobalTrackerReadoutID",
       "Name of the DD4hep constant holding the CellID encoding string."};
+  Gaudi::Property<bool> m_buildCaloSurfaces{
+      this, "BuildCaloSurfaces", true,
+      "Whether to build the ECAL inner-face surfaces (for track extrapolation to the calorimeter face)."};
 
   const CellIDSurfaceMap& cellIdToSurfaceMap() const override { return m_cellIDToSurface; }
   std::string             cellIDEncodingString() const override { return m_cellIDEncodingString; }
@@ -75,12 +80,19 @@ class ActsGeoSvc : public extends<Service, IActsGeoSvc> {
 
   using BlueprintPopulationFunc = void(const std::string&, Acts::Experimental::Blueprint&, BlueprintBuilder&);
 
+  /// Build the ECAL inner-face surfaces (m_caloFaceSurfaces) from the DD4hep
+  /// geometry. Surfaces are located via DetType flags and dimensioned from the
+  /// dd4hep::rec::LayeredCalorimeterData extension. Missing ECAL sub-detectors
+  /// or extensions are skipped with a warning rather than treated as an error.
+  void buildCaloFaceSurfaces();
+
   SmartIF<IGeoSvc>                                          m_geoSvc;
   std::shared_ptr<const Acts::TrackingGeometry>             m_trackingGeo{nullptr};
   std::shared_ptr<const Acts::MagneticFieldProvider>        m_magneticField{nullptr};
   std::unordered_map<dd4hep::CellID, const Acts::Surface*>  m_cellIDToSurface{};
   std::unordered_map<std::string, BlueprintPopulationFunc*> m_bluePrintPopulationFuncs{};
   std::string                                               m_cellIDEncodingString{};
+  CaloFaceSurfaces                                          m_caloFaceSurfaces{};
 };
 
 inline std::shared_ptr<const Acts::TrackingGeometry> ActsGeoSvc::trackingGeometry() const { return m_trackingGeo; }