From cb43097df02ba8232b91f21e836ec9f64923666e Mon Sep 17 00:00:00 2001
From: Xiaocong Ai <xiaocong.ai@cern.ch>
Date: Thu, 27 Jun 2024 22:52:48 +0200
Subject: [PATCH] include states and track writer in CKF2

---
 .../Acts/FaserActsKalmanFilter/CMakeLists.txt |    8 +-
 .../FaserActsKalmanFilter/IndexSourceLink.h   |    3 -
 .../Acts/FaserActsKalmanFilter/src/CKF2.cxx   |   64 +-
 .../Acts/FaserActsKalmanFilter/src/CKF2.h     |    2 +-
 .../src/CreateTrkTrackTool.cxx                |    6 +-
 .../FaserActsKalmanFilter/src/EffPlotTool.cxx |    2 +-
 .../FaserActsKalmanFilter/src/EffPlotTool.h   |   16 +-
 .../src/FaserActsKalmanFilterAlg.cxx          |    8 +-
 .../src/FaserActsKalmanFilterAlg.h            |    4 +-
 .../src/KalmanFitterTool.cxx                  |   17 +-
 .../src/KalmanFitterTool.h                    |    4 +-
 .../FaserActsKalmanFilter/src/ResPlotTool.cxx |    2 +-
 .../FaserActsKalmanFilter/src/ResPlotTool.h   |    2 +-
 .../src/RootTrajectoryStatesWriterTool.cxx    | 1138 +++++++++--------
 .../src/RootTrajectoryStatesWriterTool.h      |   76 +-
 15 files changed, 674 insertions(+), 678 deletions(-)

diff --git a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
index 3e6a66d5b..f52d8c1ea 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
+++ b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
@@ -55,14 +55,14 @@ atlas_add_component(FaserActsKalmanFilter
     PerformanceWriterTool.h
     PlotHelpers.h
     ResPlotTool.h
-#todo    RootTrajectoryStatesWriterTool.h
+    RootTrajectoryStatesWriterTool.h
     RootTrajectorySummaryWriterTool.h
     SeedingAlg.h
 #    SegmentFitClusterTrackFinderTool.h
 #    SegmentFitTrackFinderTool.h
     SimWriterTool.h
-#todo    SPSeedBasedInitialParameterTool.h
-#todo    SPSimpleInitialParameterTool.h
+    SPSeedBasedInitialParameterTool.h
+    SPSimpleInitialParameterTool.h
     SummaryPlotTool.h
     TrackClassification.h
     TrackSeedWriterTool.h
@@ -95,7 +95,7 @@ atlas_add_component(FaserActsKalmanFilter
     src/PlotHelpers.cxx
 #todo    src/ResPlotTool.cxx
     src/SeedingAlg.cxx
-#todo    src/RootTrajectoryStatesWriterTool.cxx
+    src/RootTrajectoryStatesWriterTool.cxx
     src/RootTrajectorySummaryWriterTool.cxx
 #    src/SegmentFitClusterTrackFinderTool.cxx
 #    src/SegmentFitTrackFinderTool.cxx
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IndexSourceLink.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IndexSourceLink.h
index a9b59c5f5..63e872c77 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IndexSourceLink.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IndexSourceLink.h
@@ -15,9 +15,6 @@
 
 #include <cassert>
 
-//toberemoved
-#include <boost/container/flat_map.hpp>
-
 #include "FaserActsGeometryContainers.h"
 // #include "TrackerPrepRawData/FaserSCT_Cluster.h"
 namespace Tracker
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index 58bc75199..bf405a669 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -53,9 +53,9 @@ StatusCode CKF2::initialize() {
 //todo  if (m_performanceWriter && !m_noDiagnostics) {
 //todo    ATH_CHECK(m_performanceWriterTool.retrieve());
 //todo  }
-//todo  if (m_statesWriter && !m_noDiagnostics) {
-//todo    ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
-//todo  }
+  if (m_statesWriter && !m_noDiagnostics) {
+    ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
+  }
   if (m_summaryWriter && !m_noDiagnostics) {
     ATH_CHECK(m_trajectorySummaryWriterTool.retrieve());
   }
@@ -139,14 +139,6 @@ StatusCode CKF2::execute() {
   Acts::ProxyAccessor<unsigned int> seedNumber("trackGroup");
   
 
-//  // Set the CombinatorialKalmanFilter options
-//  CKF2::TrackFinderOptions options(
-//      gctx, magFieldContext, calibContext,
-//      IndexSourceLinkAccessor(), MeasurementCalibrator(*measurements),
-//      Acts::MeasurementSelector(measurementSelectorCfg),
-//      Acts::LoggerWrapper{*m_logger}, pOptions, &(*initialSurface));
-
-
   //1) calibrator
   MeasurementCalibratorAdapter calibrator(MeasurementCalibrator(), *measurements);
 
@@ -237,38 +229,36 @@ StatusCode CKF2::execute() {
     });
   }
 
-/*
-  for (const FaserActsRecMultiTrajectory &traj : selectedTrajectories) {
-    const auto params = traj.trackParameters(traj.tips().front());
-    ATH_MSG_DEBUG("Fitted parameters");
-    ATH_MSG_DEBUG("  params:   " << params.parameters().transpose());
-    ATH_MSG_DEBUG("  position: " << params.position(gctx).transpose());
-    ATH_MSG_DEBUG("  momentum: " << params.momentum().transpose());
-    ATH_MSG_DEBUG("  charge:   " << params.charge());
-    std::unique_ptr<Trk::Track> track = m_createTrkTrackTool->createTrack(gctx, traj, m_backwardPropagation);
-    if (track != nullptr) {
-      m_numberOfSelectedTracks++;
-      std::unique_ptr<Trk::Track> track2 = m_kalmanFitterTool1->fit(
-        ctx, gctx, track.get(), Acts::BoundVector::Zero(), m_isMC);
-      if (track2) {
-        outputTracks->push_back(std::move(track2));
-      } else {
-        outputTracks->push_back(std::move(track));
-        ATH_MSG_WARNING("Re-Fit failed.");
-      }
+
+  for (const auto & track : selectedTracks) {
+    //if(track.hasReferenceSurface){ 
+    //  const auto& parameter = track.parameters();
+    //  ATH_MSG_DEBUG("Fitted parameters");
+    //  ATH_MSG_DEBUG("  params:   " << params.parameters().transpose());
+    //  ATH_MSG_DEBUG("  position: " << params.position(gctx).transpose());
+    //  ATH_MSG_DEBUG("  momentum: " << params.momentum().transpose());
+    //  ATH_MSG_DEBUG("  charge:   " << params.charge());
+    //}
+    std::unique_ptr<Trk::Track> trk = m_createTrkTrackTool->createTrack(gctx, track, m_backwardPropagation);
+    m_numberOfSelectedTracks++;
+    std::unique_ptr<Trk::Track> trk2 = m_kalmanFitterTool1->fit(
+      ctx, gctx, trk.get(), Acts::BoundVector::Zero(), m_isMC);
+    if (trk2) {
+      outputTracks->push_back(std::move(trk2));
     } else {
-      ATH_MSG_WARNING("CKF failed.");
+      outputTracks->push_back(std::move(trk));
+      ATH_MSG_WARNING("Re-Fit failed.");
     }
   }
 
-  // run the performance writer
+//  // @todo run the performance writer
+//  if  (m_performanceWriter && !m_noDiagnostics) {
+//    ATH_CHECK(m_performanceWriterTool->write(gctx, selectedTrajectories));
+//  }
+
   if (m_statesWriter && !m_noDiagnostics) {
-    ATH_CHECK(m_trajectoryStatesWriterTool->write(gctx, selectedTrajectories, m_isMC));
-  }
-  if  (m_performanceWriter && !m_noDiagnostics) {
-    ATH_CHECK(m_performanceWriterTool->write(gctx, selectedTrajectories));
+    ATH_CHECK(m_trajectoryStatesWriterTool->write(gctx, selectedTracks, m_isMC));
   }
-*/
   if (m_summaryWriter && !m_noDiagnostics) {
     ATH_CHECK(m_trajectorySummaryWriterTool->write(gctx, selectedTracks, m_isMC));
   }
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
index 9dff400da..a7d46feaf 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.h
@@ -104,7 +104,7 @@ private:
   ToolHandle<ITrackSeedTool> m_trackSeedTool {this, "TrackSeed", "ClusterTrackSeedTool"};
   ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
 //todo  ToolHandle<PerformanceWriterTool> m_performanceWriterTool {this, "PerformanceWriterTool", "PerformanceWriterTool"};
-//todo  ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
+  ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
   ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"};
   ToolHandle<KalmanFitterTool> m_kalmanFitterTool1 {this, "KalmanFitterTool1", "KalmanFitterTool"};
   ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
index 8a7454564..f3d772b0c 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CreateTrkTrackTool.cxx
@@ -31,19 +31,19 @@ CreateTrkTrackTool::createTrack(const Acts::GeometryContext &gctx, const CreateT
       if (flag.test(Acts::TrackStateFlag::HoleFlag)) {
         //todo: make the particle hypothesis configurable
         const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
-                                                   state.predicted(), state.predictedCovariance(), Acts::ParticleHypothesis::pion());
+                                                   state.predicted(), state.predictedCovariance(), Acts::ParticleHypothesis::muon());
         parm = std::move(ConvertActsTrackParameterToATLAS(actsParam, gctx));
         typePattern.set(Trk::TrackStateOnSurface::Hole);
       }
       else if (flag.test(Acts::TrackStateFlag::OutlierFlag)) {
         const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
-                                                   state.filtered(), state.filteredCovariance(), Acts::ParticleHypothesis::pion());
+                                                   state.filtered(), state.filteredCovariance(), Acts::ParticleHypothesis::muon());
         parm = std::move(ConvertActsTrackParameterToATLAS(actsParam, gctx));
         typePattern.set(Trk::TrackStateOnSurface::Outlier);
       }
       else {
         const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
-                                                   state.smoothed(), state.smoothedCovariance(), Acts::ParticleHypothesis::pion());
+                                                   state.smoothed(), state.smoothedCovariance(), Acts::ParticleHypothesis::muon());
         parm = std::move(ConvertActsTrackParameterToATLAS(actsParam, gctx));
         typePattern.set(Trk::TrackStateOnSurface::Measurement);
       }
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/EffPlotTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/EffPlotTool.cxx
index e757c9c8d..75be6fa30 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/EffPlotTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/EffPlotTool.cxx
@@ -17,7 +17,7 @@ void EffPlotTool::book(EffPlotTool::EffPlotCache &effPlotCache) const {
 }
 
 void EffPlotTool::fill(EffPlotTool::EffPlotCache &effPlotCache,
-                       const HepMC::GenParticle* truthParticle, bool status) const {
+                       const HepMC3::GenParticle* truthParticle, bool status) const {
   const auto t_phi = truthParticle->momentum().phi();
   const auto t_eta = truthParticle->momentum().eta();
   const auto t_pT = truthParticle->momentum().perp() * m_MeV2GeV;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/EffPlotTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/EffPlotTool.h
index 63b92ad1a..74df93ce2 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/EffPlotTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/EffPlotTool.h
@@ -3,7 +3,7 @@
 
 #include "PlotHelpers.h"
 #include "Acts/EventData/TrackParameters.hpp"
-#include "HepMC/GenParticle.h"
+#include "HepMC3/GenParticle.h"
 #include "TEfficiency.h"
 #include "TProfile.h"
 #include <map>
@@ -14,14 +14,18 @@ public:
   std::map<std::string, PlotHelpers::Binning> m_varBinning = {
       {"Eta", PlotHelpers::Binning("#eta", 40, 4, 12)},
       {"Phi", PlotHelpers::Binning("#phi", 100, -3.15, 3.15)},
-      {"Pt", PlotHelpers::Binning("pT [GeV/c]", 40, 0, 20)}
+      {"Pt", PlotHelpers::Binning("pT [GeV/c]", 40, 0, 20)},
+      {"DeltaR", PlotHelpers::Binning("#Delta R", 100, 0, 0.3)}
   };
 
   /// @brief Nested Cache struct
   struct EffPlotCache {
-    TEfficiency* trackEff_vs_pT;   ///< Tracking efficiency vs pT
-    TEfficiency* trackEff_vs_eta;  ///< Tracking efficiency vs eta
-    TEfficiency* trackEff_vs_phi;  ///< Tracking efficiency vs phi
+    TEfficiency* trackEff_vs_pT{nullptr};   ///< Tracking efficiency vs pT
+    TEfficiency* trackEff_vs_eta{nullptr};  ///< Tracking efficiency vs eta
+    TEfficiency* trackEff_vs_phi{nullptr};  ///< Tracking efficiency vs phi
+    TEfficiency* trackEff_vs_DeltaR{
+        nullptr};  ///< Tracking efficiency vs distance to the closest truth
+                   ///< particle
   };
 
   /// Constructor
@@ -38,7 +42,7 @@ public:
   /// @param effPlotCache cache object for efficiency plots
   /// @param truthParticle the truth Particle
   /// @param status the reconstruction status
-  void fill(EffPlotCache& effPlotCache, const HepMC::GenParticle* truthParticle, bool status) const;
+  void fill(EffPlotCache& effPlotCache, const HepMC3::GenParticle* truthParticle, bool status) const;
 
   /// @brief write the efficiency plots to file
   ///
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx
index 7a8d541d8..41579b41a 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.cxx
@@ -81,7 +81,7 @@ StatusCode FaserActsKalmanFilterAlg::initialize() {
   ATH_CHECK(m_trackingGeometryTool.retrieve());
   ATH_CHECK(m_trackFinderTool.retrieve());
 //todo  ATH_CHECK(m_trajectoryWriterTool.retrieve());
-//todo  ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
+  ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
 //  ATH_CHECK(m_protoTrackWriterTool.retrieve());
   ATH_CHECK(m_trackCollection.initialize());
   ATH_CHECK(detStore()->retrieve(m_idHelper,"FaserSCT_ID"));
@@ -242,7 +242,7 @@ FaserActsKalmanFilterAlg::makeTrack(Acts::GeometryContext& geoCtx, TrackFitterRe
           //@todo: ParticleHypothesis? 
           const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
                                                      state.predicted(),
-                                                     state.predictedCovariance(), Acts::ParticleHypothesis::pion());
+                                                     state.predictedCovariance(), Acts::ParticleHypothesis::muon());
           parm = std::move(ConvertActsTrackParameterToATLAS(actsParam, geoCtx));
           // auto boundaryCheck = m_boundaryCheckTool->boundaryCheck(*p arm);
           typePattern.set(Trk::TrackStateOnSurface::Hole);
@@ -251,7 +251,7 @@ FaserActsKalmanFilterAlg::makeTrack(Acts::GeometryContext& geoCtx, TrackFitterRe
         else if (flag.test(Acts::TrackStateFlag::OutlierFlag) == true) {
           //@todo: ParticleHypothesis? 
           const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
-                                                     state.filtered(), state.filteredCovariance(), Acts::ParticleHypothesis::pion());
+                                                     state.filtered(), state.filteredCovariance(), Acts::ParticleHypothesis::muon());
           parm = std::move(ConvertActsTrackParameterToATLAS(actsParam, geoCtx));
           typePattern.set(Trk::TrackStateOnSurface::Outlier);
         }
@@ -259,7 +259,7 @@ FaserActsKalmanFilterAlg::makeTrack(Acts::GeometryContext& geoCtx, TrackFitterRe
         else {
           //@todo: ParticleHypothesis? 
           const Acts::BoundTrackParameters actsParam(state.referenceSurface().getSharedPtr(),
-                                                     state.smoothed(), state.smoothedCovariance(), Acts::ParticleHypothesis::pion());
+                                                     state.smoothed(), state.smoothedCovariance(), Acts::ParticleHypothesis::muon());
           actsSmoothedParam.push_back(std::make_unique<const Acts::BoundTrackParameters>(Acts::BoundTrackParameters(actsParam)));
           //  const auto& psurface=actsParam.referenceSurface();
           Acts::Vector2 local(actsParam.parameters()[Acts::eBoundLoc0], actsParam.parameters()[Acts::eBoundLoc1]);
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h
index 9ef73935d..de7cb8b63 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/FaserActsKalmanFilterAlg.h
@@ -56,7 +56,7 @@
 #include "TrackFitterFunction.h"
 
 //#include "ProtoTrackWriterTool.h"
-//todo #include "RootTrajectoryStatesWriterTool.h"
+#include "RootTrajectoryStatesWriterTool.h"
 
 // STL
 #include <memory>
@@ -110,7 +110,7 @@ private:
   ToolHandle<ITrackFinderTool> m_trackFinderTool {this, "TrackFinderTool", "TruthTrackFinderTool"};
   ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
 //todo  ToolHandle<TrajectoryWriterTool> m_trajectoryWriterTool {this, "OutputTool", "TrajectoryWriterTool"};
-//todo  ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
+  ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
 //  ToolHandle<ProtoTrackWriterTool> m_protoTrackWriterTool {this, "ProtoTrackWriterTool", "ProtoTrackWriterTool"};
 
   SG::ReadHandleKey<McEventCollection> m_mcEventKey { this, "McEventCollection", "BeamTruthEvent" };
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
index f5f64c970..73365b135 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.cxx
@@ -22,8 +22,8 @@ StatusCode KalmanFitterTool::initialize() {
   ATH_CHECK(m_trackingGeometryTool.retrieve());
   ATH_CHECK(m_createTrkTrackTool.retrieve());
   ATH_CHECK(detStore()->retrieve(m_idHelper,"FaserSCT_ID"));
-//todo  if (m_statesWriter && !m_noDiagnostics) ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
-//todo  if (m_summaryWriter && !m_noDiagnostics) ATH_CHECK(m_trajectorySummaryWriterTool.retrieve());
+  if (m_statesWriter && !m_noDiagnostics) ATH_CHECK(m_trajectoryStatesWriterTool.retrieve());
+  if (m_summaryWriter && !m_noDiagnostics) ATH_CHECK(m_trajectorySummaryWriterTool.retrieve());
   if (m_actsLogging == "VERBOSE" && !m_noDiagnostics) {
     m_logger = Acts::getDefaultLogger("KalmanFitter", Acts::Logging::VERBOSE);
   } else if (m_actsLogging == "DEBUG" && !m_noDiagnostics) {
@@ -417,13 +417,12 @@ KalmanFitterTool::fit(const EventContext &ctx, const Acts::GeometryContext &gctx
     newTrack = m_createTrkTrackTool->createTrack(gctx, track);
   }
 
-//todo
-//  if (m_statesWriter && !m_noDiagnostics) {
-//    StatusCode statusStatesWriterTool = m_trajectoryStatesWriterTool->write(gctx, myTrajectories, isMC);
-//  }
-//  if (m_summaryWriter && !m_noDiagnostics) {
-//    StatusCode statusSummaryWriterTool = m_trajectorySummaryWriterTool->write(gctx, myTrajectories, isMC);
-//  }
+  if (m_statesWriter && !m_noDiagnostics) {
+    StatusCode statusStatesWriterTool = m_trajectoryStatesWriterTool->write(gctx, tracks, isMC);
+  }
+  if (m_summaryWriter && !m_noDiagnostics) {
+    StatusCode statusSummaryWriterTool = m_trajectorySummaryWriterTool->write(gctx, tracks, isMC);
+  }
 
   return newTrack;
 }
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
index 8484d37d8..36a1afa07 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/KalmanFitterTool.h
@@ -89,8 +89,8 @@ private:
 
   SG::ReadCondHandleKey<FaserFieldCacheCondObj> m_fieldCondObjInputKey {this, "FaserFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};
   ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
-//todo  ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
-//todo  ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"};
+  ToolHandle<RootTrajectoryStatesWriterTool> m_trajectoryStatesWriterTool {this, "RootTrajectoryStatesWriterTool", "RootTrajectoryStatesWriterTool"};
+  ToolHandle<RootTrajectorySummaryWriterTool> m_trajectorySummaryWriterTool {this, "RootTrajectorySummaryWriterTool", "RootTrajectorySummaryWriterTool"};
   ToolHandle<CreateTrkTrackTool> m_createTrkTrackTool {this, "CreateTrkTrackTool", "CreateTrkTrackTool"};
 
 //  Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> getExtensions();
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.cxx
index 78c29b9ea..0e8d08e4e 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.cxx
@@ -1,6 +1,6 @@
 #include "ResPlotTool.h"
 #include "Acts/Utilities/Helpers.hpp"
-#include "HepMC/GenVertex.h"
+#include "HepMC3/GenVertex.h"
 
 void ResPlotTool::book(ResPlotTool::ResPlotCache& resPlotCache) const {
   PlotHelpers::Binning bEta = m_varBinning.at("Eta");
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.h
index 9b27be2b0..2f3aea284 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/ResPlotTool.h
@@ -4,7 +4,7 @@
 #include "PlotHelpers.h"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
-#include "HepMC/GenParticle.h"
+#include "HepMC3/GenParticle.h"
 #include "TH1F.h"
 #include "TH2F.h"
 #include <map>
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
index c46c1a0b9..9982aa0db 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.cxx
@@ -1,3 +1,9 @@
+// Amg
+
+// Ensure Eigen plugin comes first
+#include "EventPrimitives/EventPrimitives.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+
 #include "RootTrajectoryStatesWriterTool.h"
 
 #include "TrackerPrepRawData/FaserSCT_Cluster.h"
@@ -14,8 +20,7 @@
 #include <TFile.h>
 #include <TTree.h>
 
-constexpr float NaNfloat = std::numeric_limits<float>::quiet_NaN();
-constexpr float NaNint = std::numeric_limits<int>::quiet_NaN();
+//constexpr float nan = std::numeric_limits<float>::quiet_NaN();
 
 using Acts::VectorHelpers::eta;
 using Acts::VectorHelpers::perp;
@@ -51,8 +56,7 @@ StatusCode RootTrajectoryStatesWriterTool::initialize() {
     m_outputTree = new TTree("tree", "tree");
 
     m_outputTree->Branch("event_nr", &m_eventNr);
-    m_outputTree->Branch("multiTraj_nr", &m_multiTrajNr);
-    m_outputTree->Branch("subTraj_nr", &m_subTrajNr);
+    m_outputTree->Branch("track_nr", &m_trackNr);
 
     m_outputTree->Branch("t_x", &m_t_x);
     m_outputTree->Branch("t_y", &m_t_y);
@@ -91,110 +95,148 @@ StatusCode RootTrajectoryStatesWriterTool::initialize() {
     m_outputTree->Branch("pull_y_hit", &m_pull_y_hit);
     m_outputTree->Branch("dim_hit", &m_dim_hit);
 
-    m_outputTree->Branch("nPredicted", &m_nParams[0]);
-    m_outputTree->Branch("predicted", &m_hasParams[0]);
-    m_outputTree->Branch("eLOC0_prt", &m_eLOC0[0]);
-    m_outputTree->Branch("eLOC1_prt", &m_eLOC1[0]);
-    m_outputTree->Branch("ePHI_prt", &m_ePHI[0]);
-    m_outputTree->Branch("eTHETA_prt", &m_eTHETA[0]);
-    m_outputTree->Branch("eQOP_prt", &m_eQOP[0]);
-    m_outputTree->Branch("eT_prt", &m_eT[0]);
-    m_outputTree->Branch("res_eLOC0_prt", &m_res_eLOC0[0]);
-    m_outputTree->Branch("res_eLOC1_prt", &m_res_eLOC1[0]);
-    m_outputTree->Branch("res_ePHI_prt", &m_res_ePHI[0]);
-    m_outputTree->Branch("res_eTHETA_prt", &m_res_eTHETA[0]);
-    m_outputTree->Branch("res_eQOP_prt", &m_res_eQOP[0]);
-    m_outputTree->Branch("res_eT_prt", &m_res_eT[0]);
-    m_outputTree->Branch("err_eLOC0_prt", &m_err_eLOC0[0]);
-    m_outputTree->Branch("err_eLOC1_prt", &m_err_eLOC1[0]);
-    m_outputTree->Branch("err_ePHI_prt", &m_err_ePHI[0]);
-    m_outputTree->Branch("err_eTHETA_prt", &m_err_eTHETA[0]);
-    m_outputTree->Branch("err_eQOP_prt", &m_err_eQOP[0]);
-    m_outputTree->Branch("err_eT_prt", &m_err_eT[0]);
-    m_outputTree->Branch("pull_eLOC0_prt", &m_pull_eLOC0[0]);
-    m_outputTree->Branch("pull_eLOC1_prt", &m_pull_eLOC1[0]);
-    m_outputTree->Branch("pull_ePHI_prt", &m_pull_ePHI[0]);
-    m_outputTree->Branch("pull_eTHETA_prt", &m_pull_eTHETA[0]);
-    m_outputTree->Branch("pull_eQOP_prt", &m_pull_eQOP[0]);
-    m_outputTree->Branch("pull_eT_prt", &m_pull_eT[0]);
-    m_outputTree->Branch("g_x_prt", &m_x[0]);
-    m_outputTree->Branch("g_y_prt", &m_y[0]);
-    m_outputTree->Branch("g_z_prt", &m_z[0]);
-    m_outputTree->Branch("px_prt", &m_px[0]);
-    m_outputTree->Branch("py_prt", &m_py[0]);
-    m_outputTree->Branch("pz_prt", &m_pz[0]);
-    m_outputTree->Branch("eta_prt", &m_eta[0]);
-    m_outputTree->Branch("pT_prt", &m_pT[0]);
-
-    m_outputTree->Branch("nFiltered", &m_nParams[1]);
-    m_outputTree->Branch("filtered", &m_hasParams[1]);
-    m_outputTree->Branch("eLOC0_flt", &m_eLOC0[1]);
-    m_outputTree->Branch("eLOC1_flt", &m_eLOC1[1]);
-    m_outputTree->Branch("ePHI_flt", &m_ePHI[1]);
-    m_outputTree->Branch("eTHETA_flt", &m_eTHETA[1]);
-    m_outputTree->Branch("eQOP_flt", &m_eQOP[1]);
-    m_outputTree->Branch("eT_flt", &m_eT[1]);
-    m_outputTree->Branch("res_eLOC0_flt", &m_res_eLOC0[1]);
-    m_outputTree->Branch("res_eLOC1_flt", &m_res_eLOC1[1]);
-    m_outputTree->Branch("res_ePHI_flt", &m_res_ePHI[1]);
-    m_outputTree->Branch("res_eTHETA_flt", &m_res_eTHETA[1]);
-    m_outputTree->Branch("res_eQOP_flt", &m_res_eQOP[1]);
-    m_outputTree->Branch("res_eT_flt", &m_res_eT[1]);
-    m_outputTree->Branch("err_eLOC0_flt", &m_err_eLOC0[1]);
-    m_outputTree->Branch("err_eLOC1_flt", &m_err_eLOC1[1]);
-    m_outputTree->Branch("err_ePHI_flt", &m_err_ePHI[1]);
-    m_outputTree->Branch("err_eTHETA_flt", &m_err_eTHETA[1]);
-    m_outputTree->Branch("err_eQOP_flt", &m_err_eQOP[1]);
-    m_outputTree->Branch("err_eT_flt", &m_err_eT[1]);
-    m_outputTree->Branch("pull_eLOC0_flt", &m_pull_eLOC0[1]);
-    m_outputTree->Branch("pull_eLOC1_flt", &m_pull_eLOC1[1]);
-    m_outputTree->Branch("pull_ePHI_flt", &m_pull_ePHI[1]);
-    m_outputTree->Branch("pull_eTHETA_flt", &m_pull_eTHETA[1]);
-    m_outputTree->Branch("pull_eQOP_flt", &m_pull_eQOP[1]);
-    m_outputTree->Branch("pull_eT_flt", &m_pull_eT[1]);
-    m_outputTree->Branch("g_x_flt", &m_x[1]);
-    m_outputTree->Branch("g_y_flt", &m_y[1]);
-    m_outputTree->Branch("g_z_flt", &m_z[1]);
-    m_outputTree->Branch("px_flt", &m_px[1]);
-    m_outputTree->Branch("py_flt", &m_py[1]);
-    m_outputTree->Branch("pz_flt", &m_pz[1]);
-    m_outputTree->Branch("eta_flt", &m_eta[1]);
-    m_outputTree->Branch("pT_flt", &m_pT[1]);
-
-    m_outputTree->Branch("nSmoothed", &m_nParams[2]);
-    m_outputTree->Branch("smoothed", &m_hasParams[2]);
-    m_outputTree->Branch("eLOC0_smt", &m_eLOC0[2]);
-    m_outputTree->Branch("eLOC1_smt", &m_eLOC1[2]);
-    m_outputTree->Branch("ePHI_smt", &m_ePHI[2]);
-    m_outputTree->Branch("eTHETA_smt", &m_eTHETA[2]);
-    m_outputTree->Branch("eQOP_smt", &m_eQOP[2]);
-    m_outputTree->Branch("eT_smt", &m_eT[2]);
-    m_outputTree->Branch("res_eLOC0_smt", &m_res_eLOC0[2]);
-    m_outputTree->Branch("res_eLOC1_smt", &m_res_eLOC1[2]);
-    m_outputTree->Branch("res_ePHI_smt", &m_res_ePHI[2]);
-    m_outputTree->Branch("res_eTHETA_smt", &m_res_eTHETA[2]);
-    m_outputTree->Branch("res_eQOP_smt", &m_res_eQOP[2]);
-    m_outputTree->Branch("res_eT_smt", &m_res_eT[2]);
-    m_outputTree->Branch("err_eLOC0_smt", &m_err_eLOC0[2]);
-    m_outputTree->Branch("err_eLOC1_smt", &m_err_eLOC1[2]);
-    m_outputTree->Branch("err_ePHI_smt", &m_err_ePHI[2]);
-    m_outputTree->Branch("err_eTHETA_smt", &m_err_eTHETA[2]);
-    m_outputTree->Branch("err_eQOP_smt", &m_err_eQOP[2]);
-    m_outputTree->Branch("err_eT_smt", &m_err_eT[2]);
-    m_outputTree->Branch("pull_eLOC0_smt", &m_pull_eLOC0[2]);
-    m_outputTree->Branch("pull_eLOC1_smt", &m_pull_eLOC1[2]);
-    m_outputTree->Branch("pull_ePHI_smt", &m_pull_ePHI[2]);
-    m_outputTree->Branch("pull_eTHETA_smt", &m_pull_eTHETA[2]);
-    m_outputTree->Branch("pull_eQOP_smt", &m_pull_eQOP[2]);
-    m_outputTree->Branch("pull_eT_smt", &m_pull_eT[2]);
-    m_outputTree->Branch("g_x_smt", &m_x[2]);
-    m_outputTree->Branch("g_y_smt", &m_y[2]);
-    m_outputTree->Branch("g_z_smt", &m_z[2]);
-    m_outputTree->Branch("px_smt", &m_px[2]);
-    m_outputTree->Branch("py_smt", &m_py[2]);
-    m_outputTree->Branch("pz_smt", &m_pz[2]);
-    m_outputTree->Branch("eta_smt", &m_eta[2]);
-    m_outputTree->Branch("pT_smt", &m_pT[2]);
+    m_outputTree->Branch("nPredicted", &m_nParams[ePredicted]);
+    m_outputTree->Branch("predicted", &m_hasParams[ePredicted]);
+    m_outputTree->Branch("eLOC0_prt", &m_eLOC0[ePredicted]);
+    m_outputTree->Branch("eLOC1_prt", &m_eLOC1[ePredicted]);
+    m_outputTree->Branch("ePHI_prt", &m_ePHI[ePredicted]);
+    m_outputTree->Branch("eTHETA_prt", &m_eTHETA[ePredicted]);
+    m_outputTree->Branch("eQOP_prt", &m_eQOP[ePredicted]);
+    m_outputTree->Branch("eT_prt", &m_eT[ePredicted]);
+    m_outputTree->Branch("res_eLOC0_prt", &m_res_eLOC0[ePredicted]);
+    m_outputTree->Branch("res_eLOC1_prt", &m_res_eLOC1[ePredicted]);
+    m_outputTree->Branch("res_ePHI_prt", &m_res_ePHI[ePredicted]);
+    m_outputTree->Branch("res_eTHETA_prt", &m_res_eTHETA[ePredicted]);
+    m_outputTree->Branch("res_eQOP_prt", &m_res_eQOP[ePredicted]);
+    m_outputTree->Branch("res_eT_prt", &m_res_eT[ePredicted]);
+    m_outputTree->Branch("err_eLOC0_prt", &m_err_eLOC0[ePredicted]);
+    m_outputTree->Branch("err_eLOC1_prt", &m_err_eLOC1[ePredicted]);
+    m_outputTree->Branch("err_ePHI_prt", &m_err_ePHI[ePredicted]);
+    m_outputTree->Branch("err_eTHETA_prt", &m_err_eTHETA[ePredicted]);
+    m_outputTree->Branch("err_eQOP_prt", &m_err_eQOP[ePredicted]);
+    m_outputTree->Branch("err_eT_prt", &m_err_eT[ePredicted]);
+    m_outputTree->Branch("pull_eLOC0_prt", &m_pull_eLOC0[ePredicted]);
+    m_outputTree->Branch("pull_eLOC1_prt", &m_pull_eLOC1[ePredicted]);
+    m_outputTree->Branch("pull_ePHI_prt", &m_pull_ePHI[ePredicted]);
+    m_outputTree->Branch("pull_eTHETA_prt", &m_pull_eTHETA[ePredicted]);
+    m_outputTree->Branch("pull_eQOP_prt", &m_pull_eQOP[ePredicted]);
+    m_outputTree->Branch("pull_eT_prt", &m_pull_eT[ePredicted]);
+    m_outputTree->Branch("g_x_prt", &m_x[ePredicted]);
+    m_outputTree->Branch("g_y_prt", &m_y[ePredicted]);
+    m_outputTree->Branch("g_z_prt", &m_z[ePredicted]);
+    m_outputTree->Branch("px_prt", &m_px[ePredicted]);
+    m_outputTree->Branch("py_prt", &m_py[ePredicted]);
+    m_outputTree->Branch("pz_prt", &m_pz[ePredicted]);
+    m_outputTree->Branch("eta_prt", &m_eta[ePredicted]);
+    m_outputTree->Branch("pT_prt", &m_pT[ePredicted]);
+
+    m_outputTree->Branch("nFiltered", &m_nParams[eFiltered]);
+    m_outputTree->Branch("filtered", &m_hasParams[eFiltered]);
+    m_outputTree->Branch("eLOC0_flt", &m_eLOC0[eFiltered]);
+    m_outputTree->Branch("eLOC1_flt", &m_eLOC1[eFiltered]);
+    m_outputTree->Branch("ePHI_flt", &m_ePHI[eFiltered]);
+    m_outputTree->Branch("eTHETA_flt", &m_eTHETA[eFiltered]);
+    m_outputTree->Branch("eQOP_flt", &m_eQOP[eFiltered]);
+    m_outputTree->Branch("eT_flt", &m_eT[eFiltered]);
+    m_outputTree->Branch("res_eLOC0_flt", &m_res_eLOC0[eFiltered]);
+    m_outputTree->Branch("res_eLOC1_flt", &m_res_eLOC1[eFiltered]);
+    m_outputTree->Branch("res_ePHI_flt", &m_res_ePHI[eFiltered]);
+    m_outputTree->Branch("res_eTHETA_flt", &m_res_eTHETA[eFiltered]);
+    m_outputTree->Branch("res_eQOP_flt", &m_res_eQOP[eFiltered]);
+    m_outputTree->Branch("res_eT_flt", &m_res_eT[eFiltered]);
+    m_outputTree->Branch("err_eLOC0_flt", &m_err_eLOC0[eFiltered]);
+    m_outputTree->Branch("err_eLOC1_flt", &m_err_eLOC1[eFiltered]);
+    m_outputTree->Branch("err_ePHI_flt", &m_err_ePHI[eFiltered]);
+    m_outputTree->Branch("err_eTHETA_flt", &m_err_eTHETA[eFiltered]);
+    m_outputTree->Branch("err_eQOP_flt", &m_err_eQOP[eFiltered]);
+    m_outputTree->Branch("err_eT_flt", &m_err_eT[eFiltered]);
+    m_outputTree->Branch("pull_eLOC0_flt", &m_pull_eLOC0[eFiltered]);
+    m_outputTree->Branch("pull_eLOC1_flt", &m_pull_eLOC1[eFiltered]);
+    m_outputTree->Branch("pull_ePHI_flt", &m_pull_ePHI[eFiltered]);
+    m_outputTree->Branch("pull_eTHETA_flt", &m_pull_eTHETA[eFiltered]);
+    m_outputTree->Branch("pull_eQOP_flt", &m_pull_eQOP[eFiltered]);
+    m_outputTree->Branch("pull_eT_flt", &m_pull_eT[eFiltered]);
+    m_outputTree->Branch("g_x_flt", &m_x[eFiltered]);
+    m_outputTree->Branch("g_y_flt", &m_y[eFiltered]);
+    m_outputTree->Branch("g_z_flt", &m_z[eFiltered]);
+    m_outputTree->Branch("px_flt", &m_px[eFiltered]);
+    m_outputTree->Branch("py_flt", &m_py[eFiltered]);
+    m_outputTree->Branch("pz_flt", &m_pz[eFiltered]);
+    m_outputTree->Branch("eta_flt", &m_eta[eFiltered]);
+    m_outputTree->Branch("pT_flt", &m_pT[eFiltered]);
+
+    m_outputTree->Branch("nSmoothed", &m_nParams[eSmoothed]);
+    m_outputTree->Branch("smoothed", &m_hasParams[eSmoothed]);
+    m_outputTree->Branch("eLOC0_smt", &m_eLOC0[eSmoothed]);
+    m_outputTree->Branch("eLOC1_smt", &m_eLOC1[eSmoothed]);
+    m_outputTree->Branch("ePHI_smt", &m_ePHI[eSmoothed]);
+    m_outputTree->Branch("eTHETA_smt", &m_eTHETA[eSmoothed]);
+    m_outputTree->Branch("eQOP_smt", &m_eQOP[eSmoothed]);
+    m_outputTree->Branch("eT_smt", &m_eT[eSmoothed]);
+    m_outputTree->Branch("res_eLOC0_smt", &m_res_eLOC0[eSmoothed]);
+    m_outputTree->Branch("res_eLOC1_smt", &m_res_eLOC1[eSmoothed]);
+    m_outputTree->Branch("res_ePHI_smt", &m_res_ePHI[eSmoothed]);
+    m_outputTree->Branch("res_eTHETA_smt", &m_res_eTHETA[eSmoothed]);
+    m_outputTree->Branch("res_eQOP_smt", &m_res_eQOP[eSmoothed]);
+    m_outputTree->Branch("res_eT_smt", &m_res_eT[eSmoothed]);
+    m_outputTree->Branch("err_eLOC0_smt", &m_err_eLOC0[eSmoothed]);
+    m_outputTree->Branch("err_eLOC1_smt", &m_err_eLOC1[eSmoothed]);
+    m_outputTree->Branch("err_ePHI_smt", &m_err_ePHI[eSmoothed]);
+    m_outputTree->Branch("err_eTHETA_smt", &m_err_eTHETA[eSmoothed]);
+    m_outputTree->Branch("err_eQOP_smt", &m_err_eQOP[eSmoothed]);
+    m_outputTree->Branch("err_eT_smt", &m_err_eT[eSmoothed]);
+    m_outputTree->Branch("pull_eLOC0_smt", &m_pull_eLOC0[eSmoothed]);
+    m_outputTree->Branch("pull_eLOC1_smt", &m_pull_eLOC1[eSmoothed]);
+    m_outputTree->Branch("pull_ePHI_smt", &m_pull_ePHI[eSmoothed]);
+    m_outputTree->Branch("pull_eTHETA_smt", &m_pull_eTHETA[eSmoothed]);
+    m_outputTree->Branch("pull_eQOP_smt", &m_pull_eQOP[eSmoothed]);
+    m_outputTree->Branch("pull_eT_smt", &m_pull_eT[eSmoothed]);
+    m_outputTree->Branch("g_x_smt", &m_x[eSmoothed]);
+    m_outputTree->Branch("g_y_smt", &m_y[eSmoothed]);
+    m_outputTree->Branch("g_z_smt", &m_z[eSmoothed]);
+    m_outputTree->Branch("px_smt", &m_px[eSmoothed]);
+    m_outputTree->Branch("py_smt", &m_py[eSmoothed]);
+    m_outputTree->Branch("pz_smt", &m_pz[eSmoothed]);
+    m_outputTree->Branch("eta_smt", &m_eta[eSmoothed]);
+    m_outputTree->Branch("pT_smt", &m_pT[eSmoothed]);
+
+
+    m_outputTree->Branch("nUnbiased", &m_nParams[eUnbiased]);
+    m_outputTree->Branch("unbiased", &m_hasParams[eUnbiased]);
+    m_outputTree->Branch("eLOC0_ubs", &m_eLOC0[eUnbiased]);
+    m_outputTree->Branch("eLOC1_ubs", &m_eLOC1[eUnbiased]);
+    m_outputTree->Branch("ePHI_ubs", &m_ePHI[eUnbiased]);
+    m_outputTree->Branch("eTHETA_ubs", &m_eTHETA[eUnbiased]);
+    m_outputTree->Branch("eQOP_ubs", &m_eQOP[eUnbiased]);
+    m_outputTree->Branch("eT_ubs", &m_eT[eUnbiased]);
+    m_outputTree->Branch("res_eLOC0_ubs", &m_res_eLOC0[eUnbiased]);
+    m_outputTree->Branch("res_eLOC1_ubs", &m_res_eLOC1[eUnbiased]);
+    m_outputTree->Branch("res_ePHI_ubs", &m_res_ePHI[eUnbiased]);
+    m_outputTree->Branch("res_eTHETA_ubs", &m_res_eTHETA[eUnbiased]);
+    m_outputTree->Branch("res_eQOP_ubs", &m_res_eQOP[eUnbiased]);
+    m_outputTree->Branch("res_eT_ubs", &m_res_eT[eUnbiased]);
+    m_outputTree->Branch("err_eLOC0_ubs", &m_err_eLOC0[eUnbiased]);
+    m_outputTree->Branch("err_eLOC1_ubs", &m_err_eLOC1[eUnbiased]);
+    m_outputTree->Branch("err_ePHI_ubs", &m_err_ePHI[eUnbiased]);
+    m_outputTree->Branch("err_eTHETA_ubs", &m_err_eTHETA[eUnbiased]);
+    m_outputTree->Branch("err_eQOP_ubs", &m_err_eQOP[eUnbiased]);
+    m_outputTree->Branch("err_eT_ubs", &m_err_eT[eUnbiased]);
+    m_outputTree->Branch("pull_eLOC0_ubs", &m_pull_eLOC0[eUnbiased]);
+    m_outputTree->Branch("pull_eLOC1_ubs", &m_pull_eLOC1[eUnbiased]);
+    m_outputTree->Branch("pull_ePHI_ubs", &m_pull_ePHI[eUnbiased]);
+    m_outputTree->Branch("pull_eTHETA_ubs", &m_pull_eTHETA[eUnbiased]);
+    m_outputTree->Branch("pull_eQOP_ubs", &m_pull_eQOP[eUnbiased]);
+    m_outputTree->Branch("pull_eT_ubs", &m_pull_eT[eUnbiased]);
+    m_outputTree->Branch("g_x_ubs", &m_x[eUnbiased]);
+    m_outputTree->Branch("g_y_ubs", &m_y[eUnbiased]);
+    m_outputTree->Branch("g_z_ubs", &m_z[eUnbiased]);
+    m_outputTree->Branch("px_ubs", &m_px[eUnbiased]);
+    m_outputTree->Branch("py_ubs", &m_py[eUnbiased]);
+    m_outputTree->Branch("pz_ubs", &m_pz[eUnbiased]);
+    m_outputTree->Branch("eta_ubs", &m_eta[eUnbiased]);
+    m_outputTree->Branch("pT_ubs", &m_pT[eUnbiased]);
+
+
 
     m_outputTree->Branch("chi2", &m_chi2);
   }
@@ -212,7 +254,8 @@ StatusCode RootTrajectoryStatesWriterTool::finalize() {
   return StatusCode::SUCCESS;
 }
 
-StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories, bool isMC) const {
+StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gctx, const FaserActsTrackContainer& tracks, bool isMC) const {
+  float nan = std::numeric_limits<float>::quiet_NaN();
 
   if (m_outputFile == nullptr)
     return StatusCode::SUCCESS;
@@ -235,8 +278,8 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
       return StatusCode::FAILURE;
     }
     ATH_MSG_VERBOSE("Found " << mcEvents->front()->particles_size() << " particles.");
-    for (const HepMC::GenParticle* particle : mcEvents->front()->particle_range()) {
-      particles[particle->barcode()] = particle;
+    for (const auto& particle : mcEvents->front()->particles()) {
+      particles[particle->id()] = &(*particle);
     }
 
     SG::ReadHandle<TrackerSimDataCollection> simDataHandle {m_simDataCollectionKey, ctx};
@@ -252,475 +295,436 @@ StatusCode RootTrajectoryStatesWriterTool::write(const Acts::GeometryContext& gc
     }
   }
 
-
-  // Loop over the trajectories
-  for (size_t itraj = 0; itraj < trajectories.size(); ++itraj) {
-    const auto& traj = trajectories[itraj];
-
-    if (traj.empty()) {
-      ATH_MSG_WARNING("Empty trajectories object " << itraj);
-      continue;
+  for (const auto& track : tracks) {
+    m_trackNr = track.index();
+
+    // Collect the track summary info
+    m_nMeasurements = track.nMeasurements();
+    m_nStates = track.nTrackStates();
+
+
+    // Get the majority truth particle to this track
+    int barcode = std::numeric_limits<int>::quiet_NaN();
+    // int truthQ = std::numeric_limits<int>::quiet_NaN();
+    // float truthMomentum = nan;
+    float truthLOC0 = nan;
+    float truthLOC1 = nan;
+    float truthPHI = nan;
+    float truthTHETA = nan;
+    float truthQOP = nan;
+    float truthTIME = nan;
+
+    if (isMC) {
+      // truthQ = 1;
+      // truthMomentum = 1;
+      identifyContributingParticles(*simData, track, particleHitCounts);
+      if (not particleHitCounts.empty()) {
+        // Get the barcode of the majority truth particle
+        barcode = particleHitCounts.front().particleId;
+        // Find the truth particle via the barcode
+        auto ip = particles.find(barcode);
+        if (ip != particles.end()) {
+          // const auto& particle = ip->second;
+          ATH_MSG_DEBUG("Find the truth particle with barcode = " << barcode);
+          // Get the truth particle charge
+          // FIXME find better way to access charge of simulated particle, this does not work for
+          // pions which have a positive pdg code (211) and positive charge
+          // truthQ = particle->pdg_id() > 0 ? -1 : 1;
+          // truthMomentum = particle->momentum().rho() * m_MeV2GeV;
+        } else {
+          ATH_MSG_WARNING("Truth particle with barcode = " << barcode << " not found!");
+        }
+      }
     }
 
-    // The trajectory index
-    m_multiTrajNr = itraj;
-
-    // The trajectory entry indices and the multiTrajectory
-    const auto& mj = traj.multiTrajectory();
-    const auto& trackTips = traj.tips();
-
-    // Loop over the entry indices for the subtrajectories
-    for (unsigned int isubtraj = 0; isubtraj < trackTips.size(); ++isubtraj) {
-      // The subtrajectory index
-      m_subTrajNr = isubtraj;
-      // The entry index for this subtrajectory
-      const auto& trackTip = trackTips[isubtraj];
-      // Collect the trajectory summary info
-      auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
-      m_nMeasurements = trajState.nMeasurements;
-      m_nStates = trajState.nStates;
-
-      // Get the majority truth particle to this track
-      int barcode = NaNint;
-      // int truthQ = NaNint;
-      // float truthMomentum = NaNfloat;
-      float truthLOC0 = NaNfloat;
-      float truthLOC1 = NaNfloat;
-      float truthPHI = NaNfloat;
-      float truthTHETA = NaNfloat;
-      float truthQOP = NaNfloat;
-      float truthTIME = NaNfloat;
-
-      if (isMC) {
-        // truthQ = 1;
-        // truthMomentum = 1;
-        identifyContributingParticles(*simData, traj, trackTip, particleHitCounts);
-        if (not particleHitCounts.empty()) {
-          // Get the barcode of the majority truth particle
-          barcode = particleHitCounts.front().particleId;
-          // Find the truth particle via the barcode
-          auto ip = particles.find(barcode);
-          if (ip != particles.end()) {
-            // const auto& particle = ip->second;
-            ATH_MSG_DEBUG("Find the truth particle with barcode = " << barcode);
-            // Get the truth particle charge
-            // FIXME find better way to access charge of simulated particle, this does not work for
-            // pions which have a positive pdg code (211) and positive charge
-            // truthQ = particle->pdg_id() > 0 ? -1 : 1;
-            // truthMomentum = particle->momentum().rho() * m_MeV2GeV;
-          } else {
-            ATH_MSG_WARNING("Truth particle with barcode = " << barcode << " not found!");
-          }
-        }
+    // Get the trackStates on the trajectory
+    m_nParams = {0, 0, 0, 0};
+    
+    for (const auto& state : track.trackStatesReversed()) {
+      const auto& surface = state.referenceSurface();
+ 
+      // we only fill the track states with non-outlier measurement
+      auto typeFlags = state.typeFlags();
+      if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
+        continue;
       }
 
-      // Get the trackStates on the trajectory
-      m_nParams = {0, 0, 0};
-      using ConstTrackStateProxy =
-          Acts::detail_lt::TrackStateProxy<IndexSourceLink, 6, true>;
-      mj.visitBackwards(trackTip, [&](const ConstTrackStateProxy& state) {
-        // we only fill the track states with non-outlier measurement
-        auto typeFlags = state.typeFlags();
-        if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
-          return true;
-        }
+      // get the truth hits corresponding to this trackState
+      IndexSourceLink sl = state.getUncalibratedSourceLink().template get<IndexSourceLink>();
+      const Tracker::FaserSCT_Cluster* cluster = sl.hit(); 
+      Identifier id = cluster->identify();
+      Identifier waferId = m_idHelper->wafer_id(id);
+      // if (siHitMap.count(std::make_pair(barcode, waferId)) == 0) {
+      //   ATH_MSG_WARNING("no FaserSiHit for hit with id " << id << " from particle " << barcode);
+      //   return true;
+      // }
+
+      Acts::Vector3 truthUnitDir(1,1,1);
+      if (isMC && siHitMap.count(std::make_pair(barcode, waferId)) != 0) {
+        const FaserSiHit* siHit = siHitMap.find(std::make_pair(barcode, waferId))->second;
+        HepGeom::Point3D localStartPos = siHit->localStartPosition();
+        HepGeom::Point3D localEndPos = siHit->localEndPosition();
+        HepGeom::Point3D<double> localPos = 0.5 * (localEndPos + localStartPos);
+        auto truthLocal = Acts::Vector2(localPos.y(), localPos.z());
+        const TrackerDD::SiDetectorElement* element = m_detMgr->getDetectorElement(id);
+        const HepGeom::Point3D<double> globalStartPosition =
+            Amg::EigenTransformToCLHEP(element->transformHit()) * localStartPos;
+        const HepGeom::Point3D<double> globalEndPosition =
+            Amg::EigenTransformToCLHEP(element->transformHit()) * localEndPos;
+        auto globalPosition = 0.5 * (globalStartPosition + globalEndPosition);
+        auto globalDirection = globalEndPosition - globalStartPosition;
+        truthUnitDir = Acts::Vector3(globalDirection.x(), globalDirection.y(), globalDirection.z()).normalized();
+        auto truthPos = Acts::Vector3(globalPosition.x() , globalPosition.y(), globalPosition.z());
+        // FIXME get truthQOP for each state
+
+        // fill the truth hit info
+        m_t_x.push_back(truthPos[Acts::ePos0]);
+        m_t_y.push_back(truthPos[Acts::ePos1]);
+        m_t_z.push_back(truthPos[Acts::ePos2]);
+        m_t_dx.push_back(truthUnitDir[Acts::eMom0]);
+        m_t_dy.push_back(truthUnitDir[Acts::eMom1]);
+        m_t_dz.push_back(truthUnitDir[Acts::eMom2]);
+
+        // get the truth track parameter at this track State
+        float truthLOC0 = truthLocal[Acts::ePos0];
+        float truthLOC1 = truthLocal[Acts::ePos1];
+        float truthTIME = siHit->meanTime();
+        float truthPHI = phi(truthUnitDir);
+        float truthTHETA = theta(truthUnitDir);
+
+        // fill the truth track parameter at this track State
+        m_t_eLOC0.push_back(truthLOC0);
+        m_t_eLOC1.push_back(truthLOC1);
+        m_t_ePHI.push_back(truthPHI);
+        m_t_eTHETA.push_back(truthTHETA);
+        m_t_eQOP.push_back(truthQOP);
+        m_t_eT.push_back(truthTIME);
+      } else {
+        m_t_x.push_back(nan);
+        m_t_y.push_back(nan);
+        m_t_z.push_back(nan);
+        m_t_dx.push_back(nan);
+        m_t_dy.push_back(nan);
+        m_t_dz.push_back(nan);
+        m_t_eLOC0.push_back(nan);
+        m_t_eLOC1.push_back(nan);
+        m_t_ePHI.push_back(nan);
+        m_t_eTHETA.push_back(nan);
+        m_t_eQOP.push_back(nan);
+        m_t_eT.push_back(nan);
+      }
 
-        const auto& surface = state.referenceSurface();
-
-        // get the truth hits corresponding to this trackState
-        Identifier id = state.uncalibrated().hit()->identify();
-        Identifier waferId = m_idHelper->wafer_id(id);
-        // if (siHitMap.count(std::make_pair(barcode, waferId)) == 0) {
-        //   ATH_MSG_WARNING("no FaserSiHit for hit with id " << id << " from particle " << barcode);
-        //   return true;
-        // }
-
-        if (isMC && siHitMap.count(std::make_pair(barcode, waferId)) != 0) {
-          const FaserSiHit* siHit = siHitMap.find(std::make_pair(barcode, waferId))->second;
-          HepGeom::Point3D localStartPos = siHit->localStartPosition();
-          HepGeom::Point3D localEndPos = siHit->localEndPosition();
-          HepGeom::Point3D<double> localPos = 0.5 * (localEndPos + localStartPos);
-          auto truthLocal = Acts::Vector2(localPos.y(), localPos.z());
-          const TrackerDD::SiDetectorElement* element = m_detMgr->getDetectorElement(id);
-          const HepGeom::Point3D<double> globalStartPosition =
-              Amg::EigenTransformToCLHEP(element->transformHit()) * localStartPos;
-          const HepGeom::Point3D<double> globalEndPosition =
-              Amg::EigenTransformToCLHEP(element->transformHit()) * localEndPos;
-          auto globalPosition = 0.5 * (globalStartPosition + globalEndPosition);
-          auto globalDirection = globalEndPosition - globalStartPosition;
-          auto truthUnitDir = Acts::Vector3(globalDirection.x(), globalDirection.y(), globalDirection.z()).normalized();
-          auto truthPos = Acts::Vector3(globalPosition.x() , globalPosition.y(), globalPosition.z());
-          // FIXME get truthQOP for each state
-
-          // fill the truth hit info
-          m_t_x.push_back(truthPos[Acts::ePos0]);
-          m_t_y.push_back(truthPos[Acts::ePos1]);
-          m_t_z.push_back(truthPos[Acts::ePos2]);
-          m_t_dx.push_back(truthUnitDir[Acts::eMom0]);
-          m_t_dy.push_back(truthUnitDir[Acts::eMom1]);
-          m_t_dz.push_back(truthUnitDir[Acts::eMom2]);
-
-          // get the truth track parameter at this track State
-          float truthLOC0 = truthLocal[Acts::ePos0];
-          float truthLOC1 = truthLocal[Acts::ePos1];
-          float truthTIME = siHit->meanTime();
-          float truthPHI = phi(truthUnitDir);
-          float truthTHETA = theta(truthUnitDir);
-
-          // fill the truth track parameter at this track State
-          m_t_eLOC0.push_back(truthLOC0);
-          m_t_eLOC1.push_back(truthLOC1);
-          m_t_ePHI.push_back(truthPHI);
-          m_t_eTHETA.push_back(truthTHETA);
-          m_t_eQOP.push_back(truthQOP);
-          m_t_eT.push_back(truthTIME);
-        } else {
-          m_t_x.push_back(NaNfloat);
-          m_t_y.push_back(NaNfloat);
-          m_t_z.push_back(NaNfloat);
-          m_t_dx.push_back(NaNfloat);
-          m_t_dy.push_back(NaNfloat);
-          m_t_dz.push_back(NaNfloat);
-          m_t_eLOC0.push_back(NaNfloat);
-          m_t_eLOC1.push_back(NaNfloat);
-          m_t_ePHI.push_back(NaNfloat);
-          m_t_eTHETA.push_back(NaNfloat);
-          m_t_eT.push_back(NaNfloat);
+      // get the geometry ID
+      auto geoID = surface.geometryId();
+      m_volumeID.push_back(geoID.volume());
+      m_layerID.push_back(geoID.layer());
+      m_moduleID.push_back(geoID.sensitive());
+
+      // get wafer information
+      m_station.push_back(m_idHelper->station(id));
+      m_layer.push_back(m_idHelper->layer(id));
+      m_phi_module.push_back(m_idHelper->phi_module(id));
+      m_eta_module.push_back(m_idHelper->eta_module(id));
+      m_side.push_back(m_idHelper->side(id));
+
+      // get the path length
+      m_pathLength.push_back(state.pathLength());
+
+      // fill the chi2
+      m_chi2.push_back(state.chi2());
+
+      // expand the local measurements into the full bound space
+      Acts::BoundVector meas = state.effectiveProjector().transpose() *
+                               state.effectiveCalibrated();
+      // extract local and global position
+      Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
+      Acts::Vector3 global =
+          surface.localToGlobal(gctx, local, truthUnitDir);
+
+      // fill the measurement info
+      m_lx_hit.push_back(local[Acts::ePos0]);
+      m_ly_hit.push_back(local[Acts::ePos1]);
+      m_x_hit.push_back(global[Acts::ePos0]);
+      m_y_hit.push_back(global[Acts::ePos1]);
+      m_z_hit.push_back(global[Acts::ePos2]);
+     
+
+      // lambda to get the fitted track parameters
+      auto getTrackParams = [&](unsigned int ipar)
+          -> std::optional<std::pair<Acts::BoundVector, Acts::BoundMatrix>> {
+        if (ipar == ePredicted && state.hasPredicted()) {
+          return std::make_pair(state.predicted(), state.predictedCovariance());
+        }
+        if (ipar == eFiltered && state.hasFiltered()) {
+          return std::make_pair(state.filtered(), state.filteredCovariance());
+        }
+        if (ipar == eSmoothed && state.hasSmoothed()) {
+          return std::make_pair(state.smoothed(), state.smoothedCovariance());
+        }
+        if (ipar == eUnbiased && state.hasSmoothed() && state.hasProjector()) {
+          // calculate the unbiased track parameters (i.e. fitted track
+          // parameters with this measurement removed) using Eq.(12a)-Eq.(12c)
+          // of NIMA 262, 444 (1987)
+          auto m = state.effectiveCalibrated();
+          auto H = state.effectiveProjector();
+          auto V = state.effectiveCalibratedCovariance();
+          auto K =
+              (state.smoothedCovariance() * H.transpose() *
+               (H * state.smoothedCovariance() * H.transpose() - V).inverse())
+                  .eval();
+          auto unbiasedParamsVec =
+              state.smoothed() + K * (m - H * state.smoothed());
+          auto unbiasedParamsCov =
+              state.smoothedCovariance() - K * H * state.smoothedCovariance();
+          return std::make_pair(unbiasedParamsVec, unbiasedParamsCov);
         }
+        return std::nullopt;
+      };
 
-        // get the geometry ID
-        auto geoID = surface.geometryId();
-        m_volumeID.push_back(geoID.volume());
-        m_layerID.push_back(geoID.layer());
-        m_moduleID.push_back(geoID.sensitive());
-
-        // get wafer information
-        m_station.push_back(m_idHelper->station(id));
-        m_layer.push_back(m_idHelper->layer(id));
-        m_phi_module.push_back(m_idHelper->phi_module(id));
-        m_eta_module.push_back(m_idHelper->eta_module(id));
-        m_side.push_back(m_idHelper->side(id));
-
-        // get the path length
-        m_pathLength.push_back(state.pathLength());
-
-        // expand the local measurements into the full bound space
-        Acts::BoundVector meas = state.projector().transpose() * state.calibrated();
-        // extract local and global position
-        Acts::Vector2 local(meas[Acts::eBoundLoc0], meas[Acts::eBoundLoc1]);
-        Acts::Vector3 mom(1, 1, 1);
-        Acts::Vector3 global = surface.localToGlobal(gctx, local, mom);
-
-        // fill the measurement info
-        m_lx_hit.push_back(local[Acts::ePos0]);
-        m_ly_hit.push_back(local[Acts::ePos1]);
-        m_x_hit.push_back(global[Acts::ePos0]);
-        m_y_hit.push_back(global[Acts::ePos1]);
-        m_z_hit.push_back(global[Acts::ePos2]);
-
-        // status of the fitted track parameters
-        std::array<bool, 3> hasParams = {false, false, false};
-        // optional fitted track parameters
-        std::optional<std::pair<Acts::BoundVector, Acts::BoundMatrix>>
-            trackParamsOpt = std::nullopt;
-        // lambda to get the fitted track parameters
-        auto getTrackParams = [&](unsigned int ipar) {
-          if (ipar == 0 && state.hasPredicted()) {
-            hasParams[0] = true;
-            m_nParams[0]++;
-            trackParamsOpt =
-                std::make_pair(state.predicted(), state.predictedCovariance());
-          } else if (ipar == 1 && state.hasFiltered()) {
-            hasParams[1] = true;
-            m_nParams[1]++;
-            trackParamsOpt =
-                std::make_pair(state.filtered(), state.filteredCovariance());
-          } else if (ipar == 2 && state.hasSmoothed()) {
-            hasParams[2] = true;
-            m_nParams[2]++;
-            trackParamsOpt =
-                std::make_pair(state.smoothed(), state.smoothedCovariance());
+
+      // fill the fitted track parameters
+      for (unsigned int ipar = 0; ipar < eSize; ++ipar) {
+        // get the fitted track parameters
+        auto trackParamsOpt = getTrackParams(ipar); 
+        // fill the track parameters status
+        m_hasParams[ipar].push_back(trackParamsOpt.has_value());
+
+        if (!trackParamsOpt) {
+          if (ipar == ePredicted) {
+            // push default values if no track parameters
+            m_res_x_hit.push_back(nan);
+            m_res_y_hit.push_back(nan);
+            m_err_x_hit.push_back(nan);
+            m_err_y_hit.push_back(nan);
+            m_pull_x_hit.push_back(nan);
+            m_pull_y_hit.push_back(nan);
+            m_dim_hit.push_back(0);
           }
-        };
-
-        // fill the fitted track parameters
-        for (unsigned int ipar = 0; ipar < 3; ++ipar) {
-          // get the fitted track parameters
-          getTrackParams(ipar);
-          if (trackParamsOpt) {
-            const auto& [parameters, covariance] = *trackParamsOpt;
-            if (ipar == 0) {
-              //
-              // local hit residual info
-              auto H = state.effectiveProjector();
-              auto resCov = state.effectiveCalibratedCovariance() +
-                            H * covariance * H.transpose();
-              auto res = state.effectiveCalibrated() - H * parameters;
-              m_res_x_hit.push_back(res[Acts::eBoundLoc0]);
-//              m_res_y_hit.push_back(res[Acts::eBoundLoc1]);
-              m_res_y_hit.push_back(NaNfloat);
-              m_err_x_hit.push_back(
-                  sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
-//              m_err_y_hit.push_back(
-//                  sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
-              m_err_x_hit.push_back(NaNfloat);
-              m_res_y_hit.push_back(NaNfloat);
-              m_pull_x_hit.push_back(
-                  res[Acts::eBoundLoc0] /
-                  sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
-//              m_pull_y_hit.push_back(
-//                  res[Acts::eBoundLoc1] /
-//                  sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
-              m_pull_y_hit.push_back(NaNfloat);
-              m_dim_hit.push_back(state.calibratedSize());
-            }
-
-            // track parameters
-            m_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0]);
-            m_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1]);
-            m_ePHI[ipar].push_back(parameters[Acts::eBoundPhi]);
-            m_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta]);
-            m_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP]);
-            m_eT[ipar].push_back(parameters[Acts::eBoundTime]);
-
-            // track parameters residual
-            float resPhi;
-            if (isMC) {
-              m_res_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0] - truthLOC0);
-              m_res_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1] - truthLOC1);
-              resPhi = Acts::detail::difference_periodic<float>( parameters[Acts::eBoundPhi], truthPHI,
-                                                                 static_cast<float>(2 * M_PI));
-              m_res_ePHI[ipar].push_back(resPhi);
-              m_res_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta] - truthTHETA);
-              m_res_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP] - truthQOP);
-              m_res_eT[ipar].push_back(parameters[Acts::eBoundTime] - truthTIME);
-
-              // track parameters error
-              m_err_eLOC0[ipar].push_back(
-                  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
-              m_err_eLOC1[ipar].push_back(
-                  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
-              m_err_ePHI[ipar].push_back(
-                  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
-              m_err_eTHETA[ipar].push_back(
-                  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
-              m_err_eQOP[ipar].push_back(
-                  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
-              m_err_eT[ipar].push_back(
-                  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
-
-              // track parameters pull
-              m_pull_eLOC0[ipar].push_back(
-                  (parameters[Acts::eBoundLoc0] - truthLOC0) /
-                  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
-              m_pull_eLOC1[ipar].push_back(
-                  (parameters[Acts::eBoundLoc1] - truthLOC1) /
-                  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
-              m_pull_ePHI[ipar].push_back(
-                  resPhi / sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
-              m_pull_eTHETA[ipar].push_back(
-                  (parameters[Acts::eBoundTheta] - truthTHETA) /
-                  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
-              m_pull_eQOP[ipar].push_back(
-                  (parameters[Acts::eBoundQOverP] - truthQOP) /
-                  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
-              m_pull_eT[ipar].push_back(
-                  (parameters[Acts::eBoundTime] - truthTIME) /
-                  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
-            } else {
-              if (ipar == 0) {
-                // push default values if no track parameters
-                m_res_x_hit.push_back(NaNfloat);
-                m_res_y_hit.push_back(NaNfloat);
-                m_err_x_hit.push_back(NaNfloat);
-                m_err_y_hit.push_back(NaNfloat);
-                m_pull_x_hit.push_back(NaNfloat);
-                m_pull_y_hit.push_back(NaNfloat);
-                m_dim_hit.push_back(NaNint);
-              }
-              // push default values if no track parameters
-              // m_eLOC0[ipar].push_back(NaNfloat);
-              // m_eLOC1[ipar].push_back(NaNfloat);
-              // m_ePHI[ipar].push_back(NaNfloat);
-              // m_eTHETA[ipar].push_back(NaNfloat);
-              // m_eQOP[ipar].push_back(NaNfloat);
-              // m_eT[ipar].push_back(NaNfloat);
-              m_res_eLOC0[ipar].push_back(NaNfloat);
-              m_res_eLOC1[ipar].push_back(NaNfloat);
-              m_res_ePHI[ipar].push_back(NaNfloat);
-              m_res_eTHETA[ipar].push_back(NaNfloat);
-              m_res_eQOP[ipar].push_back(NaNfloat);
-              m_res_eT[ipar].push_back(NaNfloat);
-              m_err_eLOC0[ipar].push_back(NaNfloat);
-              m_err_eLOC1[ipar].push_back(NaNfloat);
-              m_err_ePHI[ipar].push_back(NaNfloat);
-              m_err_eTHETA[ipar].push_back(NaNfloat);
-              m_err_eQOP[ipar].push_back(NaNfloat);
-              m_err_eT[ipar].push_back(NaNfloat);
-              m_pull_eLOC0[ipar].push_back(NaNfloat);
-              m_pull_eLOC1[ipar].push_back(NaNfloat);
-              m_pull_ePHI[ipar].push_back(NaNfloat);
-              m_pull_eTHETA[ipar].push_back(NaNfloat);
-              m_pull_eQOP[ipar].push_back(NaNfloat);
-              m_pull_eT[ipar].push_back(NaNfloat);
-              // m_x[ipar].push_back(NaNfloat);
-              // m_y[ipar].push_back(NaNfloat);
-              // m_z[ipar].push_back(NaNfloat);
-              // m_px[ipar].push_back(NaNfloat);
-              // m_py[ipar].push_back(NaNfloat);
-              // m_pz[ipar].push_back(NaNfloat);
-              // m_pT[ipar].push_back(NaNfloat);
-              // m_eta[ipar].push_back(NaNfloat);
-            }
-
-            // further track parameter info
-            Acts::FreeVector freeParams =
-                Acts::detail::transformBoundToFreeParameters(surface, gctx,
-                                                             parameters);
-            m_x[ipar].push_back(freeParams[Acts::eFreePos0]);
-            m_y[ipar].push_back(freeParams[Acts::eFreePos1]);
-            m_z[ipar].push_back(freeParams[Acts::eFreePos2]);
-            auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
-            m_px[ipar].push_back(p * freeParams[Acts::eFreeDir0]);
-            m_py[ipar].push_back(p * freeParams[Acts::eFreeDir1]);
-            m_pz[ipar].push_back(p * freeParams[Acts::eFreeDir2]);
-            m_pT[ipar].push_back(p * std::hypot(freeParams[Acts::eFreeDir0],
-                                                freeParams[Acts::eFreeDir1]));
-            m_eta[ipar].push_back(Acts::VectorHelpers::eta(
-                freeParams.segment<3>(Acts::eFreeDir0)));
+        // push default values if no track parameters
+          m_eLOC0[ipar].push_back(nan);
+          m_eLOC1[ipar].push_back(nan);
+          m_ePHI[ipar].push_back(nan);
+          m_eTHETA[ipar].push_back(nan);
+          m_eQOP[ipar].push_back(nan);
+          m_eT[ipar].push_back(nan);
+          m_res_eLOC0[ipar].push_back(nan);
+          m_res_eLOC1[ipar].push_back(nan);
+          m_res_ePHI[ipar].push_back(nan);
+          m_res_eTHETA[ipar].push_back(nan);
+          m_res_eQOP[ipar].push_back(nan);
+          m_res_eT[ipar].push_back(nan);
+          m_err_eLOC0[ipar].push_back(nan);
+          m_err_eLOC1[ipar].push_back(nan);
+          m_err_ePHI[ipar].push_back(nan);
+          m_err_eTHETA[ipar].push_back(nan);
+          m_err_eQOP[ipar].push_back(nan);
+          m_err_eT[ipar].push_back(nan);
+          m_pull_eLOC0[ipar].push_back(nan);
+          m_pull_eLOC1[ipar].push_back(nan);
+          m_pull_ePHI[ipar].push_back(nan);
+          m_pull_eTHETA[ipar].push_back(nan);
+          m_pull_eQOP[ipar].push_back(nan);
+          m_pull_eT[ipar].push_back(nan);
+          m_x[ipar].push_back(nan);
+          m_y[ipar].push_back(nan);
+          m_z[ipar].push_back(nan);
+          m_px[ipar].push_back(nan);
+          m_py[ipar].push_back(nan);
+          m_pz[ipar].push_back(nan);
+          m_pT[ipar].push_back(nan);
+          m_eta[ipar].push_back(nan);
+
+          continue;
+        }
+
+        ++m_nParams[ipar];
+        const auto& [parameters, covariance] = *trackParamsOpt;
+
+        // track parameters
+        m_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0]);
+        m_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1]);
+        m_ePHI[ipar].push_back(parameters[Acts::eBoundPhi]);
+        m_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta]);
+        m_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP]);
+        m_eT[ipar].push_back(parameters[Acts::eBoundTime]);
+
+        // track parameters residual
+        float resPhi;
+        if (isMC) {
+          m_res_eLOC0[ipar].push_back(parameters[Acts::eBoundLoc0] - truthLOC0);
+          m_res_eLOC1[ipar].push_back(parameters[Acts::eBoundLoc1] - truthLOC1);
+          resPhi = Acts::detail::difference_periodic<float>( parameters[Acts::eBoundPhi], truthPHI,
+                                                             static_cast<float>(2 * M_PI));
+          m_res_ePHI[ipar].push_back(resPhi);
+          m_res_eTHETA[ipar].push_back(parameters[Acts::eBoundTheta] - truthTHETA);
+          m_res_eQOP[ipar].push_back(parameters[Acts::eBoundQOverP] - truthQOP);
+          m_res_eT[ipar].push_back(parameters[Acts::eBoundTime] - truthTIME);
+
+          // track parameters error
+          m_err_eLOC0[ipar].push_back(
+              sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+          m_err_eLOC1[ipar].push_back(
+              sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
+          m_err_ePHI[ipar].push_back(
+              sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
+          m_err_eTHETA[ipar].push_back(
+              sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
+          m_err_eQOP[ipar].push_back(
+              sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
+          m_err_eT[ipar].push_back(
+              sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
+
+          // track parameters pull
+          m_pull_eLOC0[ipar].push_back(
+              (parameters[Acts::eBoundLoc0] - truthLOC0) /
+              sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+          m_pull_eLOC1[ipar].push_back(
+              (parameters[Acts::eBoundLoc1] - truthLOC1) /
+              sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
+          m_pull_ePHI[ipar].push_back(
+              resPhi / sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
+          m_pull_eTHETA[ipar].push_back(
+              (parameters[Acts::eBoundTheta] - truthTHETA) /
+              sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
+          m_pull_eQOP[ipar].push_back(
+              (parameters[Acts::eBoundQOverP] - truthQOP) /
+              sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
+          m_pull_eT[ipar].push_back(
+              (parameters[Acts::eBoundTime] - truthTIME) /
+              sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
+        } 
+
+        // further track parameter info
+        Acts::FreeVector freeParams =
+            Acts::detail::transformBoundToFreeParameters(surface, gctx,
+                                                         parameters);
+        m_x[ipar].push_back(freeParams[Acts::eFreePos0]);
+        m_y[ipar].push_back(freeParams[Acts::eFreePos1]);
+        m_z[ipar].push_back(freeParams[Acts::eFreePos2]);
+        auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
+        m_px[ipar].push_back(p * freeParams[Acts::eFreeDir0]);
+        m_py[ipar].push_back(p * freeParams[Acts::eFreeDir1]);
+        m_pz[ipar].push_back(p * freeParams[Acts::eFreeDir2]);
+        m_pT[ipar].push_back(p * std::hypot(freeParams[Acts::eFreeDir0],
+                                            freeParams[Acts::eFreeDir1]));
+        m_eta[ipar].push_back(Acts::VectorHelpers::eta(
+            freeParams.segment<3>(Acts::eFreeDir0)));
+      
+
+        if (ipar == ePredicted) {
+          // local hit residual info
+          auto H = state.effectiveProjector();
+          auto V = state.effectiveCalibratedCovariance();
+          auto resCov = V + H * covariance * H.transpose();
+          Acts::ActsDynamicVector res(state.calibratedSize());
+          res.setZero();
+
+          res = state.effectiveCalibrated() - H * parameters;
+
+          m_res_x_hit.push_back(res[Acts::eBoundLoc0]);
+          m_err_x_hit.push_back(
+              std::sqrt(V(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+          m_pull_x_hit.push_back(
+              res[Acts::eBoundLoc0] /
+              std::sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
+
+          if (state.calibratedSize() >= 2) {
+            m_res_y_hit.push_back(res[Acts::eBoundLoc1]);
+            m_err_y_hit.push_back(
+                std::sqrt(V(Acts::eBoundLoc1, Acts::eBoundLoc1)));
+            m_pull_y_hit.push_back(
+                res[Acts::eBoundLoc1] /
+                std::sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
           } else {
-            if (ipar == 0) {
-              // push default values if no track parameters
-              m_res_x_hit.push_back(NaNfloat);
-              m_res_y_hit.push_back(NaNfloat);
-              m_err_x_hit.push_back(NaNfloat);
-              m_err_y_hit.push_back(NaNfloat);
-              m_pull_x_hit.push_back(NaNfloat);
-              m_pull_y_hit.push_back(NaNfloat);
-              m_dim_hit.push_back(NaNint);
-            }
-            // push default values if no track parameters
-            m_eLOC0[ipar].push_back(NaNfloat);
-            m_eLOC1[ipar].push_back(NaNfloat);
-            m_ePHI[ipar].push_back(NaNfloat);
-            m_eTHETA[ipar].push_back(NaNfloat);
-            m_eQOP[ipar].push_back(NaNfloat);
-            m_eT[ipar].push_back(NaNfloat);
-            m_res_eLOC0[ipar].push_back(NaNfloat);
-            m_res_eLOC1[ipar].push_back(NaNfloat);
-            m_res_ePHI[ipar].push_back(NaNfloat);
-            m_res_eTHETA[ipar].push_back(NaNfloat);
-            m_res_eQOP[ipar].push_back(NaNfloat);
-            m_res_eT[ipar].push_back(NaNfloat);
-            m_err_eLOC0[ipar].push_back(NaNfloat);
-            m_err_eLOC1[ipar].push_back(NaNfloat);
-            m_err_ePHI[ipar].push_back(NaNfloat);
-            m_err_eTHETA[ipar].push_back(NaNfloat);
-            m_err_eQOP[ipar].push_back(NaNfloat);
-            m_err_eT[ipar].push_back(NaNfloat);
-            m_pull_eLOC0[ipar].push_back(NaNfloat);
-            m_pull_eLOC1[ipar].push_back(NaNfloat);
-            m_pull_ePHI[ipar].push_back(NaNfloat);
-            m_pull_eTHETA[ipar].push_back(NaNfloat);
-            m_pull_eQOP[ipar].push_back(NaNfloat);
-            m_pull_eT[ipar].push_back(NaNfloat);
-            m_x[ipar].push_back(NaNfloat);
-            m_y[ipar].push_back(NaNfloat);
-            m_z[ipar].push_back(NaNfloat);
-            m_px[ipar].push_back(NaNfloat);
-            m_py[ipar].push_back(NaNfloat);
-            m_pz[ipar].push_back(NaNfloat);
-            m_pT[ipar].push_back(NaNfloat);
-            m_eta[ipar].push_back(NaNfloat);
+            m_res_y_hit.push_back(nan);
+            m_err_y_hit.push_back(nan);
+            m_pull_y_hit.push_back(nan);
           }
-          // fill the track parameters status
-          m_hasParams[ipar].push_back(hasParams[ipar]);
+
+          m_dim_hit.push_back(state.calibratedSize());
         }
 
-        // fill the chi2
-        m_chi2.push_back(state.chi2());
-
-        return true;
-      });  // all states
-
-      // fill the variables for one track to tree
-      m_outputTree->Fill();
-
-      // now reset
-      m_t_x.clear();
-      m_t_y.clear();
-      m_t_z.clear();
-      m_t_dx.clear();
-      m_t_dy.clear();
-      m_t_dz.clear();
-      m_t_eLOC0.clear();
-      m_t_eLOC1.clear();
-      m_t_ePHI.clear();
-      m_t_eTHETA.clear();
-      m_t_eQOP.clear();
-      m_t_eT.clear();
-
-      m_volumeID.clear();
-      m_layerID.clear();
-      m_moduleID.clear();
-      m_station.clear();
-      m_layer.clear();
-      m_phi_module.clear();
-      m_eta_module.clear();
-      m_side.clear();
-      m_pathLength.clear();
-      m_lx_hit.clear();
-      m_ly_hit.clear();
-      m_x_hit.clear();
-      m_y_hit.clear();
-      m_z_hit.clear();
-      m_res_x_hit.clear();
-      m_res_y_hit.clear();
-      m_err_x_hit.clear();
-      m_err_y_hit.clear();
-      m_pull_x_hit.clear();
-      m_pull_y_hit.clear();
-      m_dim_hit.clear();
-
-      for (unsigned int ipar = 0; ipar < 3; ++ipar) {
-        m_hasParams[ipar].clear();
-        m_eLOC0[ipar].clear();
-        m_eLOC1[ipar].clear();
-        m_ePHI[ipar].clear();
-        m_eTHETA[ipar].clear();
-        m_eQOP[ipar].clear();
-        m_eT[ipar].clear();
-        m_res_eLOC0[ipar].clear();
-        m_res_eLOC1[ipar].clear();
-        m_res_ePHI[ipar].clear();
-        m_res_eTHETA[ipar].clear();
-        m_res_eQOP[ipar].clear();
-        m_res_eT[ipar].clear();
-        m_err_eLOC0[ipar].clear();
-        m_err_eLOC1[ipar].clear();
-        m_err_ePHI[ipar].clear();
-        m_err_eTHETA[ipar].clear();
-        m_err_eQOP[ipar].clear();
-        m_err_eT[ipar].clear();
-        m_pull_eLOC0[ipar].clear();
-        m_pull_eLOC1[ipar].clear();
-        m_pull_ePHI[ipar].clear();
-        m_pull_eTHETA[ipar].clear();
-        m_pull_eQOP[ipar].clear();
-        m_pull_eT[ipar].clear();
-        m_x[ipar].clear();
-        m_y[ipar].clear();
-        m_z[ipar].clear();
-        m_px[ipar].clear();
-        m_py[ipar].clear();
-        m_pz[ipar].clear();
-        m_eta[ipar].clear();
-        m_pT[ipar].clear();
-      }
+      } // loop over parameters
+
+    }// loop over states 
+    
+    // fill the variables for one track to tree
+    m_outputTree->Fill();
+
+    // now reset
+    m_t_x.clear();
+    m_t_y.clear();
+    m_t_z.clear();
+    m_t_dx.clear();
+    m_t_dy.clear();
+    m_t_dz.clear();
+    m_t_eLOC0.clear();
+    m_t_eLOC1.clear();
+    m_t_ePHI.clear();
+    m_t_eTHETA.clear();
+    m_t_eQOP.clear();
+    m_t_eT.clear();
+
+    m_volumeID.clear();
+    m_layerID.clear();
+    m_moduleID.clear();
+    m_station.clear();
+    m_layer.clear();
+    m_phi_module.clear();
+    m_eta_module.clear();
+    m_side.clear();
+    m_pathLength.clear();
+    m_lx_hit.clear();
+    m_ly_hit.clear();
+    m_x_hit.clear();
+    m_y_hit.clear();
+    m_z_hit.clear();
+    m_res_x_hit.clear();
+    m_res_y_hit.clear();
+    m_err_x_hit.clear();
+    m_err_y_hit.clear();
+    m_pull_x_hit.clear();
+    m_pull_y_hit.clear();
+    m_dim_hit.clear();
+
+    for (unsigned int ipar = 0; ipar < eSize; ++ipar) {
+      m_hasParams[ipar].clear();
+      m_eLOC0[ipar].clear();
+      m_eLOC1[ipar].clear();
+      m_ePHI[ipar].clear();
+      m_eTHETA[ipar].clear();
+      m_eQOP[ipar].clear();
+      m_eT[ipar].clear();
+      m_res_eLOC0[ipar].clear();
+      m_res_eLOC1[ipar].clear();
+      m_res_ePHI[ipar].clear();
+      m_res_eTHETA[ipar].clear();
+      m_res_eQOP[ipar].clear();
+      m_res_eT[ipar].clear();
+      m_err_eLOC0[ipar].clear();
+      m_err_eLOC1[ipar].clear();
+      m_err_ePHI[ipar].clear();
+      m_err_eTHETA[ipar].clear();
+      m_err_eQOP[ipar].clear();
+      m_err_eT[ipar].clear();
+      m_pull_eLOC0[ipar].clear();
+      m_pull_eLOC1[ipar].clear();
+      m_pull_ePHI[ipar].clear();
+      m_pull_eTHETA[ipar].clear();
+      m_pull_eQOP[ipar].clear();
+      m_pull_eT[ipar].clear();
+      m_x[ipar].clear();
+      m_y[ipar].clear();
+      m_z[ipar].clear();
+      m_px[ipar].clear();
+      m_py[ipar].clear();
+      m_pz[ipar].clear();
+      m_eta[ipar].clear();
+      m_pT[ipar].clear();
+    }
 
-      m_chi2.clear();
-    }  // all subtrajectories
-  }    // all trajectories
+    m_chi2.clear();
+  }  // all tracks 
 
   return StatusCode::SUCCESS;
 }
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.h
index 266de4bd0..dfd6890b0 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/RootTrajectoryStatesWriterTool.h
@@ -1,6 +1,7 @@
 #ifndef FASERACTSKALMANFILTER_ROOTTRAJECTORYSTATESWRITERTOOL_H
 #define FASERACTSKALMANFILTER_ROOTTRAJECTORYSTATESWRITERTOOL_H
 
+#include "FaserActsTrack.h"
 #include "GeoPrimitives/GeoPrimitives.h"
 #include "AthenaBaseComps/AthAlgTool.h"
 #include "Acts/EventData/TrackParameters.hpp"
@@ -30,7 +31,7 @@ public:
   StatusCode initialize() override;
   StatusCode finalize() override;
 
-  StatusCode write(const Acts::GeometryContext& gctx, const TrajectoriesContainer& trajectories, bool isMC) const;
+  StatusCode write(const Acts::GeometryContext& gctx, const FaserActsTrackContainer& tracks, bool isMC) const;
 
 private:
   SG::ReadHandleKey<McEventCollection> m_mcEventCollectionKey {this, "McEventCollection", "TruthEvent"};
@@ -48,9 +49,10 @@ private:
   TFile* m_outputFile;
   TTree* m_outputTree;
 
+  enum ParameterType { ePredicted, eFiltered, eSmoothed, eUnbiased, eSize };
+
   mutable uint32_t m_eventNr{0};         ///< the event number
-  mutable uint32_t m_multiTrajNr{0};     ///< the multi-trajectory number
-  mutable unsigned int m_subTrajNr{0};   ///< the multi-trajectory sub-trajectory number
+  mutable uint32_t m_trackNr{0};  ///< The track number in event
 
   mutable std::vector<float> m_t_x;  ///< Global truth hit position x
   mutable std::vector<float> m_t_y;  ///< Global truth hit position y
@@ -90,40 +92,40 @@ private:
   mutable std::vector<float> m_pull_y_hit;  ///< hit pull y
   mutable std::vector<int> m_dim_hit;       ///< dimension of measurement
 
-  mutable std::array<int, 3> m_nParams;  ///< number of states which have filtered/predicted/smoothed parameters
-  mutable std::array<std::vector<bool>, 3> m_hasParams;  ///< status of the filtered/predicted/smoothed parameters
-  mutable std::array<std::vector<float>, 3> m_eLOC0;  ///< predicted/filtered/smoothed parameter eLOC0
-  mutable std::array<std::vector<float>, 3> m_eLOC1;  ///< predicted/filtered/smoothed parameter eLOC1
-  mutable std::array<std::vector<float>, 3> m_ePHI;  ///< predicted/filtered/smoothed parameter ePHI
-  mutable std::array<std::vector<float>, 3> m_eTHETA;  ///< predicted/filtered/smoothed parameter eTHETA
-  mutable std::array<std::vector<float>, 3> m_eQOP;  ///< predicted/filtered/smoothed parameter eQOP
-  mutable std::array<std::vector<float>, 3> m_eT;  ///< predicted/filtered/smoothed parameter eT
-  mutable std::array<std::vector<float>, 3> m_res_eLOC0;  ///< predicted/filtered/smoothed parameter eLOC0 residual
-  mutable std::array<std::vector<float>, 3> m_res_eLOC1;  ///< predicted/filtered/smoothed parameter eLOC1 residual
-  mutable std::array<std::vector<float>, 3> m_res_ePHI;  ///< predicted/filtered/smoothed parameter ePHI residual
-  mutable std::array<std::vector<float>, 3> m_res_eTHETA;  ///< predicted/filtered/smoothed parameter eTHETA residual
-  mutable std::array<std::vector<float>, 3> m_res_eQOP;  ///< predicted/filtered/smoothed parameter eQOP residual
-  mutable std::array<std::vector<float>, 3> m_res_eT;  ///< predicted/filtered/smoothed parameter eT residual
-  mutable std::array<std::vector<float>, 3> m_err_eLOC0;  ///< predicted/filtered/smoothed parameter eLOC0 error
-  mutable std::array<std::vector<float>, 3> m_err_eLOC1;  ///< predicted/filtered/smoothed parameter eLOC1 error
-  mutable std::array<std::vector<float>, 3> m_err_ePHI;  ///< predicted/filtered/smoothed parameter ePHI error
-  mutable std::array<std::vector<float>, 3> m_err_eTHETA;  ///< predicted/filtered/smoothed parameter eTHETA error
-  mutable std::array<std::vector<float>, 3> m_err_eQOP;  ///< predicted/filtered/smoothed parameter eQOP error
-  mutable std::array<std::vector<float>, 3> m_err_eT;  ///< predicted/filtered/smoothed parameter eT error
-  mutable std::array<std::vector<float>, 3> m_pull_eLOC0;  ///< predicted/filtered/smoothed parameter eLOC0 pull
-  mutable std::array<std::vector<float>, 3> m_pull_eLOC1;  ///< predicted/filtered/smoothed parameter eLOC1 pull
-  mutable std::array<std::vector<float>, 3> m_pull_ePHI;  ///< predicted/filtered/smoothed parameter ePHI pull
-  mutable std::array<std::vector<float>, 3> m_pull_eTHETA;  ///< predicted/filtered/smoothed parameter eTHETA pull
-  mutable std::array<std::vector<float>, 3> m_pull_eQOP;  ///< predicted/filtered/smoothed parameter eQOP pull
-  mutable std::array<std::vector<float>, 3> m_pull_eT;  ///< predicted/filtered/smoothed parameter eT pull
-  mutable std::array<std::vector<float>, 3> m_x;  ///< predicted/filtered/smoothed parameter global x
-  mutable std::array<std::vector<float>, 3> m_y;  ///< predicted/filtered/smoothed parameter global y
-  mutable std::array<std::vector<float>, 3> m_z;  ///< predicted/filtered/smoothed parameter global z
-  mutable std::array<std::vector<float>, 3> m_px;  ///< predicted/filtered/smoothed parameter px
-  mutable std::array<std::vector<float>, 3> m_py;  ///< predicted/filtered/smoothed parameter py
-  mutable std::array<std::vector<float>, 3> m_pz;  ///< predicted/filtered/smoothed parameter pz
-  mutable std::array<std::vector<float>, 3> m_eta;  ///< predicted/filtered/smoothed parameter eta
-  mutable std::array<std::vector<float>, 3> m_pT;  ///< predicted/filtered/smoothed parameter pT
+  mutable std::array<int, eSize> m_nParams;  ///< number of states which have filtered/predicted/smoothed parameters
+  mutable std::array<std::vector<bool>, eSize> m_hasParams;  ///< status of the filtered/predicted/smoothed parameters
+  mutable std::array<std::vector<float>, eSize> m_eLOC0;  ///< predicted/filtered/smoothed parameter eLOC0
+  mutable std::array<std::vector<float>, eSize> m_eLOC1;  ///< predicted/filtered/smoothed parameter eLOC1
+  mutable std::array<std::vector<float>, eSize> m_ePHI;  ///< predicted/filtered/smoothed parameter ePHI
+  mutable std::array<std::vector<float>, eSize> m_eTHETA;  ///< predicted/filtered/smoothed parameter eTHETA
+  mutable std::array<std::vector<float>, eSize> m_eQOP;  ///< predicted/filtered/smoothed parameter eQOP
+  mutable std::array<std::vector<float>, eSize> m_eT;  ///< predicted/filtered/smoothed parameter eT
+  mutable std::array<std::vector<float>, eSize> m_res_eLOC0;  ///< predicted/filtered/smoothed parameter eLOC0 residual
+  mutable std::array<std::vector<float>, eSize> m_res_eLOC1;  ///< predicted/filtered/smoothed parameter eLOC1 residual
+  mutable std::array<std::vector<float>, eSize> m_res_ePHI;  ///< predicted/filtered/smoothed parameter ePHI residual
+  mutable std::array<std::vector<float>, eSize> m_res_eTHETA;  ///< predicted/filtered/smoothed parameter eTHETA residual
+  mutable std::array<std::vector<float>, eSize> m_res_eQOP;  ///< predicted/filtered/smoothed parameter eQOP residual
+  mutable std::array<std::vector<float>, eSize> m_res_eT;  ///< predicted/filtered/smoothed parameter eT residual
+  mutable std::array<std::vector<float>, eSize> m_err_eLOC0;  ///< predicted/filtered/smoothed parameter eLOC0 error
+  mutable std::array<std::vector<float>, eSize> m_err_eLOC1;  ///< predicted/filtered/smoothed parameter eLOC1 error
+  mutable std::array<std::vector<float>, eSize> m_err_ePHI;  ///< predicted/filtered/smoothed parameter ePHI error
+  mutable std::array<std::vector<float>, eSize> m_err_eTHETA;  ///< predicted/filtered/smoothed parameter eTHETA error
+  mutable std::array<std::vector<float>, eSize> m_err_eQOP;  ///< predicted/filtered/smoothed parameter eQOP error
+  mutable std::array<std::vector<float>, eSize> m_err_eT;  ///< predicted/filtered/smoothed parameter eT error
+  mutable std::array<std::vector<float>, eSize> m_pull_eLOC0;  ///< predicted/filtered/smoothed parameter eLOC0 pull
+  mutable std::array<std::vector<float>, eSize> m_pull_eLOC1;  ///< predicted/filtered/smoothed parameter eLOC1 pull
+  mutable std::array<std::vector<float>, eSize> m_pull_ePHI;  ///< predicted/filtered/smoothed parameter ePHI pull
+  mutable std::array<std::vector<float>, eSize> m_pull_eTHETA;  ///< predicted/filtered/smoothed parameter eTHETA pull
+  mutable std::array<std::vector<float>, eSize> m_pull_eQOP;  ///< predicted/filtered/smoothed parameter eQOP pull
+  mutable std::array<std::vector<float>, eSize> m_pull_eT;  ///< predicted/filtered/smoothed parameter eT pull
+  mutable std::array<std::vector<float>, eSize> m_x;  ///< predicted/filtered/smoothed parameter global x
+  mutable std::array<std::vector<float>, eSize> m_y;  ///< predicted/filtered/smoothed parameter global y
+  mutable std::array<std::vector<float>, eSize> m_z;  ///< predicted/filtered/smoothed parameter global z
+  mutable std::array<std::vector<float>, eSize> m_px;  ///< predicted/filtered/smoothed parameter px
+  mutable std::array<std::vector<float>, eSize> m_py;  ///< predicted/filtered/smoothed parameter py
+  mutable std::array<std::vector<float>, eSize> m_pz;  ///< predicted/filtered/smoothed parameter pz
+  mutable std::array<std::vector<float>, eSize> m_eta;  ///< predicted/filtered/smoothed parameter eta
+  mutable std::array<std::vector<float>, eSize> m_pT;  ///< predicted/filtered/smoothed parameter pT
 
   mutable  std::vector<float> m_chi2;  ///< chisq from filtering
 };
-- 
GitLab