Extrapolation to calo surface

ACTSTracking/ACTSSeededCKFTrackingProc.hxx

@@ -60,6 +60,10 @@ class ACTSSeededCKFTrackingProc : public ACTSProcBase {
   bool _runCKF = true;
   bool _propagateBackward = false;
 
+  // Extrapolation to calo settings
+  float _caloFaceR = 1857; //mm
+  float _caloFaceZ = 2307; //mm
+
   // Seed finding configuration
   float _seedFinding_rMax = 150;
   float _seedFinding_deltaRMin = 5;

ACTSTracking/Helpers.hxx

@@ -83,6 +83,11 @@ EVENT::Track* ACTS2Marlin_track(
     std::shared_ptr<Acts::MagneticFieldProvider> magneticField,
     Acts::MagneticFieldProvider::Cache& magCache);
 
+EVENT::Track* ACTS2Marlin_track(
+    const TrackResult& fitter_res,
+    std::shared_ptr<Acts::MagneticFieldProvider> magneticField,
+    Acts::MagneticFieldProvider::Cache& magCache, double caloFaceR, double caloFaceZ, Acts::GeometryContext geoContext, Acts::MagneticFieldContext magFieldContext, std::shared_ptr<const Acts::TrackingGeometry> trackingGeo);
+
 //! Convert ACTS track state class to Marlin class
 /**
  * \param location Location where the track state is defined (ie: `AtIP`)

src/ACTSSeededCKFTrackingProc.cxx

@@ -72,6 +72,15 @@ ACTSSeededCKFTrackingProc::ACTSSeededCKFTrackingProc()
                              "Track error estimate, local position (mm).",
                              _initialTrackError_pos, 10_um);
 
+  // Extrapolation to calo surface
+  registerProcessorParameter("CaloFace_Radius",
+                             "ECAL Inner Radius (mm).",
+                             _caloFaceR, _caloFaceR);
+
+  registerProcessorParameter("CaloFace_Z",
+                             "ECAL half length (mm).",
+                             _caloFaceZ, _caloFaceZ);
+
   // Seeding configurations
   registerProcessorParameter(
       "SeedingLayers",
@@ -682,7 +691,7 @@ void ACTSSeededCKFTrackingProc::processEvent(LCEvent *evt) {
 
           // Make track object
           EVENT::Track *track = ACTSTracking::ACTS2Marlin_track(
-              trackTip, magneticField(), magCache);
+              trackTip, magneticField(), magCache, _caloFaceR, _caloFaceZ, geometryContext(), magneticFieldContext(), trackingGeometry());
 
           // Save results
           trackCollection->addElement(track);

src/Helpers.cxx

@@ -6,10 +6,24 @@
 #include <UTIL/CellIDDecoder.h>
 #include <UTIL/LCTrackerConf.h>
 
+#include <Acts/Surfaces/CylinderSurface.hpp>
+#include <Acts/Surfaces/DiscSurface.hpp>
+#include <Acts/Definitions/Algebra.hpp>
+#include <Acts/EventData/ParticleHypothesis.hpp>
+
 #include <filesystem>
 
 #include "config.h"
 
+#include <Acts/Propagator/EigenStepper.hpp>
+#include <Acts/Propagator/Navigator.hpp>
+#include <Acts/Propagator/Propagator.hpp>
+
+
+using Stepper = Acts::EigenStepper<>;
+using Navigator = Acts::Navigator;
+using Propagator = Acts::Propagator<Stepper, Navigator>;
+
 namespace ACTSTracking {
 
 std::string findFile(const std::string& inpath) {
@@ -301,6 +315,161 @@ EVENT::Track* ACTS2Marlin_track(
   return track;
 }
 
+// Conversion with extrapolation to calorimeter face
+EVENT::Track* ACTS2Marlin_track(
+    const TrackResult& fitter_res,
+    std::shared_ptr<Acts::MagneticFieldProvider> magneticField,
+    Acts::MagneticFieldProvider::Cache& magCache, double caloFaceR, double caloFaceZ, 
+    Acts::GeometryContext geoContext, Acts::MagneticFieldContext magFieldContext, 
+    std::shared_ptr<const Acts::TrackingGeometry> trackingGeo)
+{
+  IMPL::TrackImpl* track = new IMPL::TrackImpl;
+
+  track->setChi2(fitter_res.chi2());
+  track->setNdf(fitter_res.nDoF());
+
+  const Acts::Vector3 zeroPos(0, 0, 0);
+  Acts::Result<Acts::Vector3> fieldRes = magneticField->getField(zeroPos, magCache);
+  if (!fieldRes.ok()) {
+    throw std::runtime_error("Field lookup error: " + fieldRes.error().value());
+  }
+  Acts::Vector3 field = *fieldRes;
+
+  const Acts::BoundVector& params = fitter_res.parameters();
+  const Acts::BoundMatrix& covariance = fitter_res.covariance();
+  EVENT::TrackState* trackStateAtIP = ACTSTracking::ACTS2Marlin_trackState(
+      EVENT::TrackState::AtIP, params, covariance, field[2] / Acts::UnitConstants::T);
+  track->trackStates().push_back(trackStateAtIP);
+
+  EVENT::TrackerHitVec hitsOnTrack;
+  EVENT::TrackStateVec statesOnTrack;
+
+  for (const auto& trk_state : fitter_res.trackStatesReversed())
+  {
+    if (!trk_state.hasUncalibratedSourceLink()) continue;
+
+    auto sl = trk_state.getUncalibratedSourceLink()
+                       .get<ACTSTracking::SourceLink>();
+    EVENT::TrackerHit* curr_hit = sl.lciohit();
+    hitsOnTrack.push_back(curr_hit);
+
+    const Acts::Vector3 hitPos(curr_hit->getPosition()[0],
+                               curr_hit->getPosition()[1],
+                               curr_hit->getPosition()[2]);
+
+    Acts::Result<Acts::Vector3> fieldRes =
+        magneticField->getField(hitPos, magCache);
+    if (!fieldRes.ok()) {
+      throw std::runtime_error("Field lookup error: " +
+                               fieldRes.error().value());
+    }
+    Acts::Vector3 field = *fieldRes;
+
+    EVENT::TrackState* trackState = ACTSTracking::ACTS2Marlin_trackState(
+            EVENT::TrackState::AtOther, trk_state.smoothed(),
+            trk_state.smoothedCovariance(), field[2] / Acts::UnitConstants::T);
+    statesOnTrack.push_back(trackState);
+  }
+
+  // Create the CaloSurface
+  auto caloCylinder = std::make_shared<Acts::CylinderBounds>(caloFaceR, caloFaceZ);
+  auto caloSurface = Acts::Surface::makeShared<Acts::CylinderSurface>(Acts::Transform3::Identity(), caloCylinder);
+
+  // Define the circle dimensions (circles at both ends of the cylinder)
+  Acts::Translation3 circlePosition1(0, 0, -caloFaceZ); // circle at -z end
+  Acts::Translation3 circlePosition2(0, 0, caloFaceZ);  // circle at +z end
+
+  // Create the circle surfaces
+  auto circleSurface1 = Acts::Surface::makeShared<Acts::DiscSurface>(Acts::Transform3(circlePosition1), 0. ,caloFaceR);
+  auto circleSurface2 = Acts::Surface::makeShared<Acts::DiscSurface>(Acts::Transform3(circlePosition2), 0., caloFaceR);
+
+  // define start parameters - swap this out with some smart call
+  double d0 = params[Acts::eBoundLoc0];
+  double z0 = params[Acts::eBoundLoc1];
+  double phi = params[Acts::eBoundPhi];
+  double theta = params[Acts::eBoundTheta];
+  double qoverp = params[Acts::eBoundQOverP];
+  double time = params[Acts::eBoundTime];
+
+  Acts::Vector3 pos(d0 * cos(phi), d0 * sin(phi), z0);
+  Acts::CurvilinearTrackParameters start(Acts::VectorHelpers::makeVector4(pos, time), phi, theta, qoverp, covariance, Acts::ParticleHypothesis::pion());
+
+  // Set propagator options
+  Acts::PropagatorOptions myCaloPropOptions(geoContext, magFieldContext);
+  myCaloPropOptions.pathLimit = 20 * Acts::UnitConstants::m;
+
+  //Let's try with our private propagator
+  // Configurations
+  Navigator::Config navigatorCfg{trackingGeo};
+  navigatorCfg.resolvePassive = false;
+  navigatorCfg.resolveMaterial = true;
+  navigatorCfg.resolveSensitive = true;
+
+  // construct all components for the fitter
+  Stepper stepper(magneticField);
+  Navigator navigator(navigatorCfg);
+  Propagator mypropagator(std::move(stepper), std::move(navigator));
+  auto resultProp = mypropagator.propagate(start, *caloSurface, myCaloPropOptions);
+  if (resultProp.ok()) {
+    auto end_parameters = resultProp.value().endParameters;
+    const Acts::BoundMatrix& atCaloCovariance = *(end_parameters->covariance());
+
+    EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState(
+        EVENT::TrackState::AtCalorimeter, end_parameters->parameters(), atCaloCovariance, field[2] / Acts::UnitConstants::T);
+    track->trackStates().push_back(trackStateAtCalo);
+  }
+  else {
+    if (theta > M_PI/2){
+        auto resultPropP = mypropagator.propagate(start, *circleSurface1, myCaloPropOptions);
+        if (resultPropP.ok()) {
+          auto end_parametersP = resultPropP.value().endParameters;
+          const Acts::BoundMatrix& atCaloCovariance = *(end_parametersP->covariance());
+
+          EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState(
+            EVENT::TrackState::AtCalorimeter, end_parametersP->parameters(), atCaloCovariance, field[2] / Acts::UnitConstants::T);
+          track->trackStates().push_back(trackStateAtCalo);
+        }
+        else{
+          std::cout << "Failed propagation! " << std::endl;
+        }
+    }
+    else{
+        auto resultPropM = mypropagator.propagate(start, *circleSurface2, myCaloPropOptions);
+        if (resultPropM.ok()) {
+          auto end_parametersM = resultPropM.value().endParameters;
+          const Acts::BoundMatrix& atCaloCovariance = *(end_parametersM->covariance());
+
+          EVENT::TrackState* trackStateAtCalo = ACTSTracking::ACTS2Marlin_trackState(
+            EVENT::TrackState::AtCalorimeter, end_parametersM->parameters(), atCaloCovariance, field[2] / Acts::UnitConstants::T);
+          track->trackStates().push_back(trackStateAtCalo);
+        }
+        else{
+          std::cout << "Failed propagation!" << std::endl;
+        }      
+    }
+  }
+
+  std::reverse(hitsOnTrack.begin(), hitsOnTrack.end());
+  std::reverse(statesOnTrack.begin(), statesOnTrack.end());
+
+  for (EVENT::TrackerHit* hit : hitsOnTrack) {
+    track->addHit(hit);
+  }
+
+  if (statesOnTrack.size() > 0) {
+    dynamic_cast<IMPL::TrackStateImpl*>(statesOnTrack.back())
+        ->setLocation(EVENT::TrackState::AtLastHit);
+    dynamic_cast<IMPL::TrackStateImpl*>(statesOnTrack.front())
+        ->setLocation(EVENT::TrackState::AtFirstHit);
+  }
+
+  EVENT::TrackStateVec& myTrackStates = track->trackStates();
+  myTrackStates.insert(myTrackStates.end(), statesOnTrack.begin(),
+                       statesOnTrack.end());
+
+  return track;
+}
+
 
 
 EVENT::TrackState* ACTS2Marlin_trackState(