diff --git a/MuonSpectrometer/MuonCnv/MuonPrepRawDataProviderTools/src/MuonPRDSelectionTool.cxx b/MuonSpectrometer/MuonCnv/MuonPrepRawDataProviderTools/src/MuonPRDSelectionTool.cxx
index 1157a8049399e8bfe16a1193d42c19f859eb705c..51106b60ec58a0f4248a8f5a59b65bf4794ca5c1 100644
--- a/MuonSpectrometer/MuonCnv/MuonPrepRawDataProviderTools/src/MuonPRDSelectionTool.cxx
+++ b/MuonSpectrometer/MuonCnv/MuonPrepRawDataProviderTools/src/MuonPRDSelectionTool.cxx
@@ -75,14 +75,12 @@ namespace Muon {
   }
 
 
-  const MdtDriftCircleOnTrack* MuonPRDSelectionTool::calibrateAndSelect( const MuonSystemExtension::Intersection& intersection, const MdtPrepData& mdt ) const {
+  const MdtDriftCircleOnTrack* MuonPRDSelectionTool::calibrateAndSelect( const MuonSystemExtension::Intersection& intersection, const MdtPrepData& mdt, double beta) const {
 
     // calculate intersection with tube
     const Amg::Vector3D& direction = intersection.trackParameters->momentum();
     Amg::Vector3D intersect = intersectMDT(mdt,intersection.trackParameters->position(),direction,true);
 
-    // get the error in the precision plane
-    double err_precision = Amg::error(*intersection.trackParameters->covariance(),Trk::locX);
 
     // calculate local position of the intersection in tube frame
     const Identifier& id = mdt.identify();
@@ -99,7 +97,7 @@ namespace Muon {
     double distanceAlongTube = localPosition[Trk::locZ];
 
     if( msgLvl(MSG::VERBOSE) ) msg(MSG::VERBOSE) << " Intersected " << m_idHelperSvc->toString(id) << " distance to wire " << localPosition[Trk::locR] 
-                                                 << " error " << err_precision << " along tube (%) " << distanceAlongTube/tubeHalfLen;
+                                                 << " error " << Amg::error(*intersection.trackParameters->covariance(),Trk::locX) << " along tube (%) " << distanceAlongTube/tubeHalfLen;
 
     if( std::abs(distanceAlongTube) > tubeHalfLen + m_secondCoordinateCut ) {
       if( msgLvl(MSG::VERBOSE) ) msg(MSG::VERBOSE) << " outside tube second coordinate range, dropping " << endmsg;
@@ -120,7 +118,7 @@ namespace Muon {
     if( msgLvl(MSG::VERBOSE) ) msg(MSG::VERBOSE) << endmsg;
 
     // calibrate hit
-    const MdtDriftCircleOnTrack* mdtROT = m_mdtCreator->createRIO_OnTrack( mdt, intersect,  &direction ); 
+    const MdtDriftCircleOnTrack* mdtROT = m_mdtCreator->createRIO_OnTrack( mdt, intersect,  &direction, 0., nullptr, beta, 0.); 
     if( !mdtROT ) ATH_MSG_VERBOSE(" Failed to calibrate " << m_idHelperSvc->toString(id));
     return mdtROT;    
   }
diff --git a/MuonSpectrometer/MuonCnv/MuonPrepRawDataProviderTools/src/MuonPRDSelectionTool.h b/MuonSpectrometer/MuonCnv/MuonPrepRawDataProviderTools/src/MuonPRDSelectionTool.h
index 1bf85f5454ca62697ef9060ccbc602a80ec62de0..b963ea1813e2d5df5799aaa5acacd0715a0e44c9 100644
--- a/MuonSpectrometer/MuonCnv/MuonPrepRawDataProviderTools/src/MuonPRDSelectionTool.h
+++ b/MuonSpectrometer/MuonCnv/MuonPrepRawDataProviderTools/src/MuonPRDSelectionTool.h
@@ -36,7 +36,7 @@ namespace Muon {
     bool calibrateAndSelect( const MuonSystemExtension::Intersection& intersection, const MuonLayerPrepRawData& layerPrepRawData, MuonLayerROTs& layerROTs ) const;
 
     /** calibrate and select MDTs in a collection */
-    bool calibrateAndSelectMdt( const MuonSystemExtension::Intersection& intersection, const MdtPrepDataCollection& prds, std::vector<const MdtDriftCircleOnTrack*>& rots ) const {
+    bool calibrateAndSelectMdt( const MuonSystemExtension::Intersection& intersection, const MdtPrepDataCollection& prds, std::vector<const MdtDriftCircleOnTrack*>& rots) const {
       for( MdtPrepDataCollection::const_iterator it = prds.begin(); it != prds.end();++it ){
         const MdtDriftCircleOnTrack* rot = calibrateAndSelect(intersection,**it);
         if( rot ) rots.push_back(rot);
@@ -55,7 +55,7 @@ namespace Muon {
     }
     
     /** IMuonPRDSelectionTool interface: calibrate and select single MDT */
-    const MdtDriftCircleOnTrack* calibrateAndSelect( const MuonSystemExtension::Intersection& intersection, const MdtPrepData& mdt ) const;
+    const MdtDriftCircleOnTrack* calibrateAndSelect( const MuonSystemExtension::Intersection& intersection, const MdtPrepData& mdt, double beta = 1.) const;
 
     /** IMuonPRDSelectionTool interface: calibrate and select single cluster */
     const MuonClusterOnTrack* calibrateAndSelect( const Trk::TrackParameters& pars, const MuonCluster& clus ) const;
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonPRDSelectionTool.h b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonPRDSelectionTool.h
index 69b314acdf29ae881a6a168451e2c38865fba4d2..a7b9584ec2ac52bf07612b86c073ac768fdaa781 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonPRDSelectionTool.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonPRDSelectionTool.h
@@ -37,7 +37,7 @@ namespace Muon {
 
         /** IMuonPRDSelectionTool interface:  calibrate and select single MDT */
         virtual const MdtDriftCircleOnTrack* calibrateAndSelect(const MuonSystemExtension::Intersection& intersection,
-                                                                const MdtPrepData& mdt) const = 0;
+                                                                const MdtPrepData& mdt, double beta = 1.) const = 0;
 
         /** IMuonPRDSelectionTool interface:  calibrate and select single cluster */
         virtual const MuonClusterOnTrack* calibrateAndSelect(const Trk::TrackParameters& pars, const MuonCluster& clus) const = 0;
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonSegmentMaker.h b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonSegmentMaker.h
index 900c51f4b2883d04bd8da1088160df3d9d4f1d3d..992ff3dc7d95e225e9731ed6ed924cc811023385 100755
--- a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonSegmentMaker.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonSegmentMaker.h
@@ -109,7 +109,7 @@ namespace Muon {
         */
         virtual void find(const Amg::Vector3D& gpos, const Amg::Vector3D& gdir, const std::vector<const MdtDriftCircleOnTrack*>& mdts,
                           const std::vector<const MuonClusterOnTrack*>& clusters, bool updatePhi = false,
-                          Trk::SegmentCollection* segColl = nullptr, double momentum = 1e9, double sinAngleCut = 0.) const = 0;
+                          Trk::SegmentCollection* segColl = nullptr, double momentum = 1e9, double sinAngleCut = 0., double beta = 1.) const = 0;
 
         /** @brief seeded segment search starting from a list of MdtDriftCircleOnTrack objects and a list of MuonClusterOnTrack objects
             @param road an estimate of the position and direction of the muon in the chamber
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonTGRecTools/src/MuonSystemExtensionTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonTGRecTools/src/MuonSystemExtensionTool.cxx
index 639fe1d07020294a2136ef1ebb993a2b54e54e8f..e34f4c3db0b06d11f3539763c63e13e436e401a6 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonTGRecTools/src/MuonSystemExtensionTool.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonTGRecTools/src/MuonSystemExtensionTool.cxx
@@ -1,10 +1,11 @@
 /*
-  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "MuonSystemExtensionTool.h"
-
 #include "EventPrimitives/EventPrimitivesHelpers.h"
+#include "GeoPrimitives/GeoPrimitivesHelpers.h"
+#include "GeoPrimitives/GeoPrimitivesToStringConverter.h"
 #include "EventPrimitives/EventPrimitivesCovarianceHelpers.h"
 #include "MuonCombinedEvent/InDetCandidate.h"
 #include "MuonCombinedEvent/TagBase.h"
@@ -15,20 +16,26 @@
 #include "TrkSurfaces/DiscSurface.h"
 #include "TrkSurfaces/PlaneSurface.h"
 
-namespace {
-    std::string to_string(const Trk::TrackParameters& pars) {
-        std::stringstream sstr;
-        sstr<<"position: "<<pars.position().x()<<", "<<pars.position().y()<<", "<<pars.position().z()
-            <<" -- momentum: "<<pars.momentum().x()<<", "<<pars.momentum().y()<<", "<<pars.momentum().z()<<" ";
-        if (pars.covariance())
-            sstr<<" -- uncertainties Lx:"<<Amg::error(*pars.covariance(), Trk::locX)<<
-                                 ", Ly: "<< Amg::error(*pars.covariance(), Trk::locY);
-        return sstr.str();
+namespace{
+    template <class T>
+    std::string to_string(const std::vector<T>& v) {
+        std::stringstream ostr{};
+        ostr<<"{";
+        for (auto x: v)ostr<<x<<", ";
+        ostr<<"}";
+        return ostr.str();
+    }
+    /** @brief Calculate the path from the origin along the track parameters */
+    double pathAlongPars(const Trk::TrackParameters& pars, const Amg::Vector3D& pos) {
+        const Amg::Vector3D parDir{pars.momentum().unit()};
+        return (pars.position().dot(parDir) > 0 ? 1 : -1) * (parDir.dot(pos));
     }
 }
-
 namespace Muon {
 
+
+
+
     MuonSystemExtensionTool::MuonSystemExtensionTool(const std::string& type, const std::string& name, const IInterface* parent) :
         AthAlgTool(type, name, parent) {
         declareInterface<IMuonSystemExtensionTool>(this);
@@ -45,17 +52,14 @@ namespace Muon {
 
     bool MuonSystemExtensionTool::initializeGeometry() {
         // initialize reference surfaces hash vectors per region and sector
-        m_referenceSurfaces.resize(MuonStationIndex::DetectorRegionIndexMax);
-        for (auto& vec : m_referenceSurfaces) vec.resize(MuonStationIndex::numberOfSectors());
 
         // first build the barrel, we need different reference planes for all sectors
         ATH_MSG_DEBUG("Initializing barrel ");
         for (unsigned int sector = 1; sector <= 16; ++sector) {
             // get rotation into sector frame
             double sectorPhi = m_sectorMapping.sectorPhi(sector);
-            Amg::AngleAxis3D sectorRotation(sectorPhi, Amg::Vector3D(0., 0., 1.));
+            const Amg::Transform3D sectorRotation(Amg::getRotateZ3D(sectorPhi));
             if (!initializeGeometryBarrel(sector, sectorRotation)) return false;
-
             if (!initializeGeometryEndcap(sector, MuonStationIndex::EndcapA, sectorRotation)) return false;
             if (!initializeGeometryEndcap(sector, MuonStationIndex::EndcapC, sectorRotation)) return false;
         }
@@ -64,7 +68,7 @@ namespace Muon {
     }
 
     bool MuonSystemExtensionTool::initializeGeometryEndcap(int sector, MuonStationIndex::DetectorRegionIndex regionIndex,
-                                                           const Amg::AngleAxis3D& sectorRotation) {
+                                                           const Amg::Transform3D& sectorRotation) {
         ATH_MSG_DEBUG("Initializing endcap: sector " << sector << " " << MuonStationIndex::regionName(regionIndex));
 
         SurfaceVec& surfaces = m_referenceSurfaces[regionIndex][sector - 1];
@@ -77,38 +81,29 @@ namespace Muon {
             // calculate reference position from radial position of the layer
             MuonStationIndex::LayerIndex layer = MuonStationIndex::toLayerIndex(stLayer);
             MuonChamberLayerDescriptor layerDescriptor = chamberLayerDescription.getDescriptor(sector, regionIndex, layer);
-            Amg::Vector3D globalPosition(0., 0., layerDescriptor.referencePosition);
-
             // reference transform + surface
-            Amg::Transform3D trans(sectorRotation);  //*Amg::AngleAxis3D(xToZRotation,Amg::Vector3D(0.,1.,0.)
-            trans.pretranslate(globalPosition);
-
-            Trk::PlaneSurface* surface = new Trk::PlaneSurface(trans);
-            MuonLayerSurface data(MuonLayerSurface::SurfacePtr(surface), sector, regionIndex, layer);
-            surfaces.push_back(std::move(data));
-
+            Amg::Transform3D trans{Amg::getTranslateZ3D(layerDescriptor.referencePosition) *sectorRotation};  //*Amg::AngleAxis3D(xToZRotation,Amg::Vector3D(0.,1.,0.)
+            std::unique_ptr<Trk::PlaneSurface> surface = std::make_unique<Trk::PlaneSurface>(trans);
             // sanity checks
             if (msgLvl(MSG::VERBOSE)) {
-                for (unsigned int test = 0; test < 2; ++test) {
-                    Amg::Vector3D lpos3d = surface->transform().inverse() * globalPosition;
-                    Amg::Vector2D lpos;
-                    surface->globalToLocal(globalPosition, globalPosition, lpos);
-
-                    ATH_MSG_VERBOSE(" sector " << sector << " layer " << MuonStationIndex::layerName(layer) << " hit x "
-                                               << globalPosition.x() << " hit z " << globalPosition.z() << " lpos3d " << lpos3d.x() << " "
-                                               << lpos3d.y() << " " << lpos3d.z() << " lpos " << lpos[Trk::loc1] << " " << lpos[Trk::loc2]
-                                               << " center " << surface->center().z() << " normal: phi " << surface->normal().phi()
-                                               << " theta " << surface->normal().theta() << " normal: x " << surface->normal().x() << " y "
-                                               << surface->normal().y() << " z " << surface->normal().z());
-                    globalPosition[0] += 100;
-                }
+                ATH_MSG_VERBOSE("initializeGeometryEndcap() --  sector " << sector << " layer " << MuonStationIndex::layerName(layer) << " surface "
+                                            <<Amg::toString(surface->transform())
+                                            << " center " << Amg::toString(surface->center()) << " theta " << surface->normal().theta() 
+                                            << " normal:  " <<Amg::toString(surface->normal()));
+                
+                
             }
+
+            MuonLayerSurface data(std::move(surface), sector, regionIndex, layer);
+            surfaces.push_back(std::move(data));
+
+
         }
         ATH_MSG_VERBOSE(" Total number of surfaces " << surfaces.size());
         return true;
     }
 
-    bool MuonSystemExtensionTool::initializeGeometryBarrel(int sector, const Amg::AngleAxis3D& sectorRotation) {
+    bool MuonSystemExtensionTool::initializeGeometryBarrel(int sector, const Amg::Transform3D& sectorRotation) {
         MuonChamberLayerDescription chamberLayerDescription;
 
         SurfaceVec& surfaces = m_referenceSurfaces[MuonStationIndex::Barrel][sector - 1];
@@ -125,28 +120,22 @@ namespace Muon {
             Amg::Vector3D globalPosition = sectorRotation * positionInSector;
 
             // reference transform + surface
-            Amg::Transform3D trans(sectorRotation * Amg::AngleAxis3D(xToZRotation, Amg::Vector3D(0., 1., 0.)));
+            Amg::Transform3D trans(sectorRotation * Amg::getRotateY3D(xToZRotation));
             trans.pretranslate(globalPosition);
-            Trk::PlaneSurface* surface = new Trk::PlaneSurface(trans);
-
-            MuonLayerSurface data(MuonLayerSurface::SurfacePtr(surface), sector, MuonStationIndex::Barrel, layer);
-            surfaces.push_back(std::move(data));
-
+            std::unique_ptr<Trk::PlaneSurface> surface = std::make_unique<Trk::PlaneSurface>(trans);
             // sanity checks
             if (msgLvl(MSG::VERBOSE)) {
-                Amg::Vector3D lpos3d = surface->transform().inverse() * globalPosition;
-                // Amg::Vector3D lpos3d2 = surface->transform().inverse()*globalPosition2;
-                Amg::Vector2D lpos;
-                surface->globalToLocal(globalPosition, globalPosition, lpos);
                 double sectorPhi = m_sectorMapping.sectorPhi(sector);
-
                 ATH_MSG_VERBOSE(" sector " << sector << " layer " << MuonStationIndex::layerName(layer) << " phi " << sectorPhi
                                            << " ref theta " << globalPosition.theta() << " phi " << globalPosition.phi() << " r "
-                                           << globalPosition.perp() << " pos " << globalPosition.x() << " " << globalPosition.y() << " "
-                                           << globalPosition.z() << " lpos3d " << lpos3d.x() << " " << lpos3d.y() << " " << lpos3d.z()
-                                           << " normal: x " << surface->normal().x() << " y " << surface->normal().y() << " z "
-                                           << surface->normal().z());
+                                           << globalPosition.perp() << " pos " << Amg::toString(globalPosition)
+                                           << " lpos3d " << Amg::toString(surface->transform().inverse() * globalPosition) 
+                                           << " normal: " << Amg::toString(surface->normal()));
             }
+            MuonLayerSurface data(std::move(surface), sector, MuonStationIndex::Barrel, layer);
+            surfaces.push_back(std::move(data));
+
+
         }
         return true;
     }
@@ -188,49 +177,46 @@ namespace Muon {
         std::vector<MuonSystemExtension::Intersection> intersections;
 
         // garbage collection
-        std::vector<std::shared_ptr<const Trk::TrackParameters> > trackParametersVec;
+        std::vector<std::shared_ptr<Trk::TrackParameters> > trackParametersVec;
 
         // loop over reference surfaces
         for (const Muon::MuonLayerSurface& it : surfaces) {
             // extrapolate to next layer
             const Trk::Surface& surface = *it.surfacePtr;
-            ATH_MSG_VERBOSE(" startPars: "<<to_string(*currentPars));
+            ATH_MSG_VERBOSE(" startPars: "<<m_printer->print(*currentPars));
             ATH_MSG_VERBOSE(" destination: sector " << it.sector << "  " << MuonStationIndex::regionName(it.regionIndex) << "  "
                                                     << MuonStationIndex::layerName(it.layerIndex) << " phi  " << surface.center().phi()
                                                     << " r " << surface.center().perp() << " z " << surface.center().z());
 
-            std::unique_ptr<const Trk::TrackParameters> exPars{
-                m_extrapolator->extrapolate(ctx, *currentPars, surface, Trk::alongMomentum, false, Trk::muon)};
+            std::unique_ptr<Trk::TrackParameters> exPars{m_extrapolator->extrapolate(ctx, *currentPars, surface, 
+                                                                                     Trk::alongMomentum, false, Trk::muon)};
             if (!exPars) {
                 ATH_MSG_VERBOSE("extrapolation failed, trying next layer ");
                 continue;
             }
-            ATH_MSG_DEBUG("Extrapolated in event "<<ctx.eventID().event_number()<<" track @ "<<to_string(*exPars));
+            ATH_MSG_DEBUG("Extrapolated in event "<<ctx.eventID().event_number()<<" track @ "<<m_printer->print(*exPars));
 
             //    reject intersections with very big uncertainties (parallel to surface)
-            constexpr double TenMeter = 10. * Gaudi::Units::meter;
             const double sigma_lx = Amg::error(*exPars->covariance(), Trk::locX);
             const double sigma_ly = Amg::error(*exPars->covariance(), Trk::locY);
-            if (!Amg::hasPositiveDiagElems(*exPars->covariance()) || std::max(sigma_lx, sigma_ly) > TenMeter){
-                ATH_MSG_DEBUG("Reject bad extrapolation "<<to_string(*exPars));
+            if (!Amg::hasPositiveDiagElems(*exPars->covariance()) || std::max(sigma_lx, sigma_ly) > 1. * Gaudi::Units::meter){
+                ATH_MSG_DEBUG("Reject bad extrapolation "<<m_printer->print(*exPars));
                 continue;
             }
 
             // create shared pointer and add to garbage collection
-            std::shared_ptr<const Trk::TrackParameters> sharedPtr{std::move(exPars)};
+            std::shared_ptr<Trk::TrackParameters> sharedPtr{std::move(exPars)};
             trackParametersVec.emplace_back(sharedPtr);
 
-            MuonSystemExtension::Intersection intersection = makeInterSection(sharedPtr, it);
-            if (intersection.trackParameters) intersections.push_back(std::move(intersection));
-
+            intersections.emplace_back(sharedPtr, it);
             constexpr double TenCm = 10. * Gaudi::Units::cm;
             if(std::hypot(sigma_lx, sigma_ly) < std::max(TenCm, 10. * std::hypot(Amg::error(*currentPars->covariance(), Trk::locX),
                                                                                  Amg::error(*currentPars->covariance(), Trk::locY)))) {
                 // only update the parameters if errors don't blow up
                 currentPars = sharedPtr.get();
             } else {
-                ATH_MSG_DEBUG("Extrapolation reached at "<<to_string(*sharedPtr)<<" has larger uncertainties than "<<
-                                 to_string(*currentPars));
+                ATH_MSG_DEBUG("Extrapolation reached at "<<m_printer->print(*sharedPtr)<<" has larger uncertainties than "<<
+                              m_printer->print(*currentPars));
             }
         }
         ATH_MSG_DEBUG("Completed extrapolation: destinations " << surfaces.size() << " intersections " << intersections.size());
@@ -239,38 +225,10 @@ namespace Muon {
                          <<" will accept the candidate: "<< (!cache.requireSystemExtension ? "si": "no"));
             return !cache.requireSystemExtension;
         }
-        cache.candidate->setExtension(
-            std::make_unique<MuonSystemExtension>(cache.candidate->getCaloExtension()->muonEntryLayerIntersection(), std::move(intersections)));
+        cache.candidate->setExtension(std::make_unique<MuonSystemExtension>(cache.candidate->getCaloExtension()->muonEntryLayerIntersection(), 
+                                                                            std::move(intersections)));
         return true;
     }
-    MuonSystemExtension::Intersection MuonSystemExtensionTool::makeInterSection(const std::shared_ptr<const Trk::TrackParameters>& sharedPtr, const MuonLayerSurface& it) const{
-        if (sharedPtr && it.regionIndex != MuonStationIndex::Barrel) {
-            // in the endcaps we need to update the sector and check that we are in the correct sector
-            double phi = sharedPtr->position().phi();
-            std::vector<int> sectors;
-            m_sectorMapping.getSectors(phi, sectors);
-
-            // loop over sectors, if the surface sector and the sector in the loop are either both large or small, pick
-            // the sector
-            int sector = -1;
-            for (const int& sec : sectors) {
-                if (it.sector % 2 == sec % 2) {
-                    sector = sec;
-                    break;
-                }
-            }
-            ATH_MSG_DEBUG(" sector " << it.sector << " selected " << sector << " nsectors " << sectors);
-            // only select if we found a matching sector in the set
-            if (sector != -1) {
-                MuonSystemExtension::Intersection intersection{sharedPtr, it};
-                intersection.layerSurface.sector = sector;
-                return intersection;
-            }
-            return MuonSystemExtension::Intersection{nullptr, it};
-        }
-        return  MuonSystemExtension::Intersection{sharedPtr, it};
-    }
-
     MuonSystemExtensionTool::SurfaceVec MuonSystemExtensionTool::getSurfacesForIntersection(const Trk::TrackParameters& muonEntryPars,
                                                                                             const SystemExtensionCache& cache) const {
         // if in endcaps pick endcap surfaces
@@ -304,90 +262,123 @@ namespace Muon {
             surfaces.insert(surfaces.end(), m_referenceSurfaces[regionIndex][sector - 1].begin(),
                             m_referenceSurfaces[regionIndex][sector - 1].end());
         }
-        auto sortFunction = (regionIndex == MuonStationIndex::Barrel)
-                                ? [](const MuonLayerSurface& s1,
-                                     const MuonLayerSurface& s2) { return s1.surfacePtr->center().perp() < s2.surfacePtr->center().perp(); }
-                                : [](const MuonLayerSurface& s1, const MuonLayerSurface& s2) {
-                                      return std::abs(s1.surfacePtr->center().z()) < std::abs(s2.surfacePtr->center().z());
-                                  };
-
-        std::stable_sort(surfaces.begin(), surfaces.end(), sortFunction);
+        std::stable_sort(surfaces.begin(), surfaces.end(), 
+                         [&muonEntryPars](const MuonLayerSurface& s1, const MuonLayerSurface& s2) {
+                            return std::abs(pathAlongPars(muonEntryPars,s1.surfacePtr->center())) <
+                                   std::abs(pathAlongPars(muonEntryPars,s2.surfacePtr->center()));
+                            
+                         });
+        if (msgLvl(MSG::VERBOSE)) {
+            for (auto& s1 : surfaces) {
+                ATH_MSG_VERBOSE("Surface "<<Muon::MuonStationIndex::layerName(s1.layerIndex)
+                              <<", center: "<<Amg::toString(s1.surfacePtr->center())
+                              <<", pathAlongPars "<<pathAlongPars(muonEntryPars,s1.surfacePtr->center())
+                              <<std::endl<<(*s1.surfacePtr));
+
+        	}
+        }
         return surfaces;
     }
     bool MuonSystemExtensionTool::muonLayerInterSections(const EventContext& ctx,
-                                            const MuonCombined::TagBase& cmb_tag,
-                                            SystemExtensionCache& cache) const{
+                                                         const MuonCombined::TagBase& cmbTag,
+                                                         SystemExtensionCache& cache) const{
 
-        const Trk::Track* cmb_trk = cmb_tag.primaryTrack();
-        if (!cmb_trk){
-            ATH_MSG_WARNING("A combined tag without any track? Please check "<<cmb_tag.toString());
+        const Trk::Track* cmbTrk = cmbTag.primaryTrack();
+        if (!cmbTrk){
+            ATH_MSG_WARNING("A combined tag without any track? Please check "<<cmbTag.toString());
             return false;
         }
         /// Get the proper surfaces for the intersections
-        const Trk::TrackParameters* entry_pars = cache.candidate->getCaloExtension()->muonEntryLayerIntersection();
-        const Trk::TrackStates* const cmb_tsos = cmb_trk->trackStateOnSurfaces();
-        /// Copy all track parameters in the MS into a separate vector
-        std::vector<const Trk::TrackStateOnSurface*> ms_tsos;
-        ms_tsos.reserve(cmb_tsos->size());
-        std::copy_if(cmb_tsos->begin(), cmb_tsos->end(), ms_tsos.end(),[this, entry_pars](const Trk::TrackStateOnSurface* tsos){
-            if  (!(tsos->type(Trk::TrackStateOnSurface::Perigee) || tsos->type(Trk::TrackStateOnSurface::Parameter)))
-             return false;
-            const Trk::TrackParameters* params = tsos->trackParameters();
-            if (!params) return false;
-            const Trk::Surface& surf = params->associatedSurface();
-            if (!surf.associatedDetectorElement()) {
-                return entry_pars->associatedSurface().normal().dot(params->position() - entry_pars->associatedSurface().center())>0;
+        const Trk::TrackParameters* entryPars = cache.candidate->getCaloExtension()->muonEntryLayerIntersection();
+        std::vector<const Trk::TrackStateOnSurface*> cmbParVec{};
+
+        ATH_MSG_VERBOSE("Calo entry parameters "<<m_printer->print(*entryPars));
+        std::vector<Muon::MuonSystemExtension::Intersection> intersections{};
+        Muon::MuonLayerSurface lastSurf{};
+        
+        const Trk::TrackStates::const_iterator end_itr = cmbTrk->trackStateOnSurfaces()->end();
+        for (Trk::TrackStates::const_iterator itr = cmbTrk->trackStateOnSurfaces()->begin(); itr != end_itr; ++itr) {
+            const Trk::TrackStateOnSurface* msTSOS{*itr};
+            ATH_MSG_VERBOSE("Track state "<<m_printer->print(*msTSOS));
+
+            if (!msTSOS->measurementOnTrack()) {
+                continue;
+            } 
+            const Identifier measId = m_edmHelperSvc->getIdentifier(*msTSOS->measurementOnTrack());
+            const Trk::TrackParameters& msPars{*msTSOS->trackParameters()};
+            /** The parameters are not muon paramaters or the parameters are before the entrance parameters */
+            if (!m_idHelperSvc->isMuon(measId) || pathAlongPars(msPars, msPars.position() - entryPars->position())< 0.){
+                continue;
+            }
+               
+            Muon::MuonStationIndex::DetectorRegionIndex regionIdx = m_idHelperSvc->regionIndex(measId);
+            Muon::MuonStationIndex::LayerIndex layerIdx = m_idHelperSvc->layerIndex(measId);
+            if (layerIdx == Muon::MuonStationIndex::BarrelExtended) {
+                regionIdx = Muon::MuonStationIndex::DetectorRegionIndex::Barrel;
+                layerIdx = Muon::MuonStationIndex::Inner;
             }
-            return m_idHelperSvc->isMuon(surf.associatedDetectorElement()->identify());
-        });
-
-        if (ms_tsos.empty()){
-            ATH_MSG_VERBOSE("A combined track without track parameters in the MS. THat should not happen");
-            return false;
-        }
-
-        SurfaceVec surfaces = getSurfacesForIntersection(*entry_pars, cache);
-        if (surfaces.empty()) {
-            ATH_MSG_DEBUG("MS extenion without surfaces.. "<<cache.candidate->toString());
-            return false;
-        }
-
-        std::vector<Muon::MuonSystemExtension::Intersection> intersections;
-
-        std::vector<const Trk::TrackStateOnSurface*>::const_iterator itr = ms_tsos.begin();
-        std::vector<const Trk::TrackStateOnSurface*>::const_iterator end  = ms_tsos.end();
-        for (const Muon::MuonLayerSurface& it : surfaces) {
-            const Trk::Surface& surf = *(it.surfacePtr);
-            /// Calculate the minimum distance to the plane
-            itr = std::min_element(itr, end, [&surf](const Trk::TrackStateOnSurface* a, const Trk::TrackStateOnSurface* b) {
-                return std::abs(surf.normal().dot(a->trackParameters()->position() - surf.center())) <
-                       std::abs(surf.normal().dot(b->trackParameters()->position() - surf.center()));
-            });
-            const Trk::TrackParameters* trk_pars = (*itr)->trackParameters();
-
-            const Amg::Vector3D pos = trk_pars->position();
-            const bool isBarrel = it.regionIndex == MuonStationIndex::Barrel;
-            ATH_MSG_VERBOSE("Distance of closest parameters to the surface: "<< std::sqrt(std::abs(surf.normal().dot(pos - surf.center()))));
-
-
 
-            const Trk::PropDirection dir = (isBarrel && pos.perp2() < surf.center().perp2()) ||
-                                           (!isBarrel && std::abs(pos.z()) < std::abs(surf.center().z()) ) ? Trk::alongMomentum : Trk::oppositeMomentum;
-            std::shared_ptr<const Trk::TrackParameters> exPars{m_extrapolator->extrapolate(ctx, *trk_pars, surf, dir, false, Trk::muon)};
+            const int sector = m_sectorMapping.getSector(msTSOS->measurementOnTrack()->associatedSurface().center().phi());
+       
+            /// The track parameters belong the same surface. Do nothing
+            if (lastSurf.layerIndex == layerIdx && lastSurf.regionIndex == regionIdx && lastSurf.sector == sector) {
+                continue;
+            }
+            const SurfaceVec& refSurfaces = m_referenceSurfaces[regionIdx][sector - 1];
+            SurfaceVec::const_iterator surfItr = std::find_if(refSurfaces.begin(), refSurfaces.end(),
+                                                    [&layerIdx](const MuonLayerSurface& surf){
+                                                    return surf.layerIndex == layerIdx;
+                                                    });
+            if (surfItr == refSurfaces.end()) {
+                ATH_MSG_WARNING("Failed to find a reference surface matching to combined parameters "<<m_printer->print(*msTSOS)
+                                <<", sector: "<<sector <<", layer: "<<Muon::MuonStationIndex::layerName(layerIdx)
+                                <<", region index: "<<Muon::MuonStationIndex::regionName(regionIdx));
+                continue;
+            }
+            lastSurf = (*surfItr);
+            const Trk::Surface& target{*lastSurf.surfacePtr};
+                
+            /// Check whether there're measurement parameters that are closer to the target surface
+            for (Trk::TrackStates::const_iterator closePars = itr+1; closePars != end_itr; ++closePars) {
+                const Trk::TrackStateOnSurface * tsosInChamb{*closePars};
+                if (!tsosInChamb->measurementOnTrack()) continue; 
+                
+                const Trk::TrackParameters& chPars{*tsosInChamb->trackParameters()};
+                if (std::abs(pathAlongPars(chPars, chPars.position() - target.center())) >
+                    std::abs(pathAlongPars(*msTSOS->trackParameters(), 
+                                            msTSOS->trackParameters()->position() - target.center()))) {
+                    break;
+                }
+                itr =  closePars -1;
+                msTSOS = *itr;
+            }
+                
+            std::unique_ptr<Trk::TrackParameters> exPars{m_extrapolator->extrapolate(ctx, *msTSOS->trackParameters(), 
+                                                                                     target, Trk::anyDirection, false, Trk::muon)};
             if (!exPars) {
-                ATH_MSG_WARNING("Failed to extrapolate track @ "<<to_string(*trk_pars)<<" to surface "<<std::endl<<surf);
+                ATH_MSG_VERBOSE(__LINE__<<" - Failed to extrapolate track @ "<<m_printer->print(*msTSOS)
+                                <<", sector: "<<sector
+                                <<", layer: "<<Muon::MuonStationIndex::layerName(layerIdx)
+                                <<", region index: "<<Muon::MuonStationIndex::regionName(regionIdx)
+                                <<" to surface "<<std::endl<<target<<std::endl
+                                <<", sector: "<<lastSurf.sector
+                                <<", layer: "<<Muon::MuonStationIndex::layerName(lastSurf.layerIndex)
+                                <<", region index: "<<Muon::MuonStationIndex::regionName(lastSurf.regionIndex)
+                                <<" pathAlongPars "<<pathAlongPars(*msTSOS->trackParameters(), target.center())
+                                <<", dir: "<<Amg::toString(msTSOS->trackParameters()->momentum().unit()));
                 continue;
             }
-            MuonSystemExtension::Intersection intersect = makeInterSection(exPars,it);
-            if (intersect.trackParameters) intersections.push_back(std::move(intersect));
+            intersections.emplace_back(std::move(exPars), lastSurf);
+            
         }
-        /// We are trying to reconstruct the layer extensions of a combined track.
-        /// There is no way that this can fail
+        
+        ATH_MSG_VERBOSE("Selected "<<intersections.size()<<" intersections");
+        
         if (intersections.empty()) {
             ATH_MSG_DEBUG("Failed to find the intersections for the combined track");
             return false;
         }
-        cache.candidate->setExtension(std::make_unique<MuonSystemExtension>(entry_pars, std::move(intersections)));
+        cache.candidate->setExtension(std::make_unique<MuonSystemExtension>(entryPars, std::move(intersections)));
 
         return true;
     }
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonTGRecTools/src/MuonSystemExtensionTool.h b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonTGRecTools/src/MuonSystemExtensionTool.h
index 67c486205e865b087ce5a0356e29085ee4b66661..e89b49035a54b1e04a3f49694d86fc92b81a4025 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonTGRecTools/src/MuonSystemExtensionTool.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonTGRecTools/src/MuonSystemExtensionTool.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef MUON_MUONSYSTEMEXTENSIONTOOL_H
@@ -12,6 +12,7 @@
 #include "MuonDetDescrUtils/MuonSectorMapping.h"
 #include "MuonLayerEvent/MuonLayerSurface.h"
 #include "MuonRecToolInterfaces/IMuonSystemExtensionTool.h"
+#include "MuonRecHelperTools/MuonEDMPrinterTool.h"
 #include "MuonStationIndex/MuonStationIndex.h"
 #include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
 #include "TrkExInterfaces/IExtrapolator.h"
@@ -47,32 +48,30 @@ namespace Muon {
     private:
         /** initialize geometry */
         bool initializeGeometry();
-        bool initializeGeometryBarrel(int sector, const Amg::AngleAxis3D& sectorRotation);
+        bool initializeGeometryBarrel(int sector, const Amg::Transform3D& sectorRotation);
         bool initializeGeometryEndcap(int sector, MuonStationIndex::DetectorRegionIndex regionIndex,
-                                      const Amg::AngleAxis3D& sectorRotation);
+                                      const Amg::Transform3D& sectorRotation);
 
         /** get surfaces to be intersected for a given start parameters */
         SurfaceVec getSurfacesForIntersection(const Trk::TrackParameters& muonEntryPars, const SystemExtensionCache& cache) const;
-
-        MuonSystemExtension::Intersection makeInterSection(const std::shared_ptr<const Trk::TrackParameters>& pars, const MuonLayerSurface& surf) const;
-
         
         ToolHandle<Trk::IParticleCaloExtensionTool> m_caloExtensionTool{
             this,
             "ParticleCaloExtensionTool",
             "Trk::ParticleCaloExtensionTool/ParticleCaloExtensionTool",
         };
-        ToolHandle<Trk::IExtrapolator> m_extrapolator{
-            this,
-            "Extrapolator",
-            "Trk::Extrapolator/AtlasExtrapolator",
-        };
+        ToolHandle<Trk::IExtrapolator> m_extrapolator{this,"Extrapolator",""};
+        
+        PublicToolHandle<MuonEDMPrinterTool> m_printer{this, "Printer", "Muon::MuonEDMPrinterTool/MuonEDMPrinterTool"};
         
         ServiceHandle<Muon::IMuonIdHelperSvc> m_idHelperSvc{this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"};
 
-      
+        ServiceHandle<IMuonEDMHelperSvc> m_edmHelperSvc{this,"edmHelper","Muon::MuonEDMHelperSvc/MuonEDMHelperSvc",
+                                                        "Handle to the service providing the IMuonEDMHelperSvc interface"};
+
         /** reference surfaces per region and sector */
-        std::vector<std::vector<SurfaceVec> > m_referenceSurfaces;
+        std::array<std::array<SurfaceVec, 16> , 
+                   MuonStationIndex::DetectorRegionIndexMax > m_referenceSurfaces{};
 
         /** sector mapping helper */
         MuonSectorMapping m_sectorMapping;
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/DCMathSegmentMaker.cxx b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/DCMathSegmentMaker.cxx
index faee75ebe1e7ca4065dd93dfd56b4bbfbedc3e92..c4125a0097ada53726cc0244c0691a31a74da406 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/DCMathSegmentMaker.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/DCMathSegmentMaker.cxx
@@ -90,7 +90,7 @@ namespace Muon {
     void DCMathSegmentMaker::find(const Amg::Vector3D& roadpos, const Amg::Vector3D& roaddir,
                                   const std::vector<const MdtDriftCircleOnTrack*>& mdts,
                                   const std::vector<const MuonClusterOnTrack*>& clusters, bool hasPhiMeasurements,
-                                  Trk::SegmentCollection* segColl, double momentum, double sinAngleCut) const {
+                                  Trk::SegmentCollection* segColl, double momentum, double sinAngleCut, double beta) const {
         const EventContext& ctx = Gaudi::Hive::currentContext();
         if (m_doTimeOutChecks && Athena::Timeout::instance().reached()) {
             ATH_MSG_DEBUG("Timeout reached. Aborting sequence.");
@@ -224,7 +224,7 @@ namespace Muon {
         // loop over segments
         segmentCreationInfo sInfo(spVecs, multiGeo.get(), gToStation, amdbToGlobal, phimin, phimax);
         for (TrkDriftCircleMath::Segment& seg : segs) {
-            std::unique_ptr<MuonSegment> segment = createSegment(ctx, seg, chid, roadpos, roaddir2, mdts, hasPhiMeasurements, sInfo);
+            std::unique_ptr<MuonSegment> segment = createSegment(ctx, seg, chid, roadpos, roaddir2, mdts, hasPhiMeasurements, sInfo, beta);
             if (segment) segColl->push_back(segment.release());
         }
         ATH_MSG_DEBUG(" Done ");
@@ -233,7 +233,7 @@ namespace Muon {
     std::unique_ptr<MuonSegment> DCMathSegmentMaker::createSegment(const EventContext& ctx, TrkDriftCircleMath::Segment& segment, const Identifier& chid,
                                                    const Amg::Vector3D& roadpos, const Amg::Vector3D& roaddir2,
                                                    const std::vector<const MdtDriftCircleOnTrack*>& mdts, bool hasPhiMeasurements,
-                                                   segmentCreationInfo& sInfo) const {
+                                                   segmentCreationInfo& sInfo, double beta) const {
         bool isEndcap = m_idHelperSvc->isEndcap(chid);
         // find all curved segments
         MuonStationIndex::ChIndex chIndex = m_idHelperSvc->chamberIndex(chid);
@@ -323,7 +323,7 @@ namespace Muon {
         std::set<Identifier> deltaVec;
         std::set<Identifier> outoftimeVec;
 
-        associateMDTsToSegment(gdir, segment, mdts, sInfo.geom, sInfo.globalTrans, sInfo.amdbTrans, deltaVec, outoftimeVec, rioDistVec);
+        associateMDTsToSegment(gdir, segment, mdts, sInfo.geom, sInfo.globalTrans, sInfo.amdbTrans, deltaVec, outoftimeVec, rioDistVec, beta);
         std::vector<std::pair<double, std::unique_ptr<const Trk::MeasurementBase>>> garbage_collector;
 
         TrkDriftCircleMath::DCSLHitSelector hitSelector;
@@ -1192,7 +1192,7 @@ namespace Muon {
         const Amg::Vector3D& gdir, TrkDriftCircleMath::Segment& segment, const std::vector<const MdtDriftCircleOnTrack*>& mdts,
         const TrkDriftCircleMath::ChamberGeometry* multiGeo, const Amg::Transform3D& gToStation, const Amg::Transform3D& amdbToGlobal,
         std::set<Identifier>& deltaVec, std::set<Identifier>& outoftimeVec,
-        std::vector<std::pair<double,  std::unique_ptr<const Trk::MeasurementBase>> >& rioDistVec) const {
+        std::vector<std::pair<double,  std::unique_ptr<const Trk::MeasurementBase>> >& rioDistVec, double beta) const {
         // clear result vectors
 
         // convert segment parameters + x position from road
@@ -1264,7 +1264,7 @@ namespace Muon {
             bool hasT0 = segment.hasT0Shift();
             if (!hasT0 || m_mdtCreatorT0.empty()) {
                 // ATH_MSG_VERBOSE(" recalibrate MDT hit");
-                nonconstDC.reset(m_mdtCreator->createRIO_OnTrack(*riodc->prepRawData(), mdtGP, &gdir));
+                nonconstDC.reset(m_mdtCreator->createRIO_OnTrack(*riodc->prepRawData(), mdtGP, &gdir, 0., nullptr, beta, 0.));
                 if (hasT0) ATH_MSG_WARNING("Attempted to change t0 without a properly configured MDT creator tool. ");
             } else {
                 ATH_MSG_VERBOSE(" recalibrate MDT hit with shift " << segment.t0Shift()<<" "<<m_printer->print(*riodc));
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/DCMathSegmentMaker.h b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/DCMathSegmentMaker.h
index df6c4877612d86d550cf21ee79dbe44a7814da5c..230c12cbe5648e42ba2b339aeacdcbc92e22be2a 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/DCMathSegmentMaker.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/DCMathSegmentMaker/src/DCMathSegmentMaker.h
@@ -246,7 +246,7 @@ namespace Muon {
         */
         void find(const Amg::Vector3D& gpos, const Amg::Vector3D& gdir, const std::vector<const MdtDriftCircleOnTrack*>& mdts,
                   const std::vector<const MuonClusterOnTrack*>& clusters, bool hasPhiMeasurements = false,
-                  Trk::SegmentCollection* segColl = nullptr, double momentum = 1e9, double sinAngleCut = 0) const;
+                  Trk::SegmentCollection* segColl = nullptr, double momentum = 1e9, double sinAngleCut = 0, double beta = 1. ) const;
 
         /** find segments starting from:
             - a track prediction
@@ -296,7 +296,7 @@ namespace Muon {
             const Amg::Vector3D& gdir, TrkDriftCircleMath::Segment& segment, const std::vector<const MdtDriftCircleOnTrack*>& mdts,
             const TrkDriftCircleMath::ChamberGeometry* multiGeo, const Amg::Transform3D& gToStation, const Amg::Transform3D& amdbToGlobal,
             std::set<Identifier>& deltaVec, std::set<Identifier>& outoftimeVec,
-            std::vector<std::pair<double,  std::unique_ptr<const Trk::MeasurementBase>> >& rioDistVec) const;
+            std::vector<std::pair<double,  std::unique_ptr<const Trk::MeasurementBase>> >& rioDistVec, double beta = 1. ) const;
         std::pair<std::pair<int, int>, bool> associateClustersToSegment(
             const TrkDriftCircleMath::Segment& segment, const Identifier& chid, const Amg::Transform3D& gToStation, ClusterVecPair& spVecs,
             double phimin, double phimax, std::vector<std::pair<double, std::unique_ptr<const Trk::MeasurementBase>> >& rioDistVec) const;
@@ -336,7 +336,7 @@ namespace Muon {
 
         std::unique_ptr<MuonSegment> createSegment(const EventContext& ctx, TrkDriftCircleMath::Segment& segment, const Identifier& chid, const Amg::Vector3D& roadpos,
                                    const Amg::Vector3D& roaddir2, const std::vector<const MdtDriftCircleOnTrack*>& mdts,
-                                   bool hasPhiMeasurements, segmentCreationInfo& sInfo) const;
+                                   bool hasPhiMeasurements, segmentCreationInfo& sInfo, double beta = 1. ) const;
 
         const MdtPrepData* findMdt(const EventContext& ctx, const Identifier& id) const;
 
diff --git a/Reconstruction/MuonIdentification/MuonCombinedAlgs/src/MuonInDetToMuonSystemExtensionAlg.cxx b/Reconstruction/MuonIdentification/MuonCombinedAlgs/src/MuonInDetToMuonSystemExtensionAlg.cxx
index 6305c1390cd78e13ef53d9f20f17feba8aa79569..b496195260553c3d4e2a7b8ea4ede442c55cd7e6 100644
--- a/Reconstruction/MuonIdentification/MuonCombinedAlgs/src/MuonInDetToMuonSystemExtensionAlg.cxx
+++ b/Reconstruction/MuonIdentification/MuonCombinedAlgs/src/MuonInDetToMuonSystemExtensionAlg.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "MuonInDetToMuonSystemExtensionAlg.h"
@@ -368,6 +368,9 @@ StatusCode MuonInDetToMuonSystemExtensionAlg::createStaus(const EventContext& ct
     }
 
     SG::WriteHandle<InDetCandidateCollection> indetCandidateCollection(m_stauInDetCandKey, ctx);
+    // sort candidates beofre storing, otherwise ordering in container of stau segments can be inconsistent
+    std::sort(stau_cache.outputContainer->begin(),stau_cache.outputContainer->end(),[](const MuonCombined::InDetCandidate*a, const MuonCombined::InDetCandidate*b){
+        return a->indetTrackParticle().pt() < b->indetTrackParticle().pt();});
     ATH_CHECK(indetCandidateCollection.record(std::move(stau_cache.outputContainer)));
     return StatusCode::SUCCESS;
 }
diff --git a/Reconstruction/MuonIdentification/MuonCombinedConfig/python/MuonCombinedRecToolsConfig.py b/Reconstruction/MuonIdentification/MuonCombinedConfig/python/MuonCombinedRecToolsConfig.py
index d693d33fa04522561e8e0dd178f3ccada7ffb188..fe05a2d4bcd799ad9ee64da3ba2952a5bc65d06e 100644
--- a/Reconstruction/MuonIdentification/MuonCombinedConfig/python/MuonCombinedRecToolsConfig.py
+++ b/Reconstruction/MuonIdentification/MuonCombinedConfig/python/MuonCombinedRecToolsConfig.py
@@ -106,11 +106,10 @@ def MuonMaterialProviderToolCfg(flags,  name="MuonTrkMaterialProviderTool", **kw
 
 
 def MuonSegmentHitSummaryToolCfg(flags, name="MuonSegmentHitSummaryTool", **kwargs):
-    from MuonConfig.MuonGeometryConfig import MuonDetectorCondAlgCfg
-
-    result = MuonEDMPrinterToolCfg(flags)
-    kwargs.setdefault("Printer", result.getPrimary())
-    result.merge(MuonDetectorCondAlgCfg(flags))
+    from MuonConfig.MuonGeometryConfig import MuonGeoModelCfg
+    result = ComponentAccumulator()
+    kwargs.setdefault("Printer", result.getPrimaryAndMerge(MuonEDMPrinterToolCfg(flags)))
+    result.merge(MuonGeoModelCfg(flags))
     kwargs.setdefault("DetectorManagerKey", "MuonDetectorManager")
     tool = CompFactory.Muon.MuonSegmentHitSummaryTool(name, **kwargs)
     result.setPrivateTools(tool)
@@ -1165,14 +1164,15 @@ def MuonSystemExtensionToolCfg(flags, **kwargs):
     result = ComponentAccumulator()
 
     from TrackToCalo.TrackToCaloConfig import ParticleCaloExtensionToolCfg
-    particle_calo_extension_tool = result.popToolsAndMerge(
-        ParticleCaloExtensionToolCfg(flags, name='MuonParticleCaloExtensionTool'))
+    from MuonConfig.MuonRecToolsConfig import MuonEDMPrinterToolCfg
 
-    atlas_extrapolator = result.popToolsAndMerge(AtlasExtrapolatorCfg(flags))
+    kwargs.setdefault("Extrapolator", result.popToolsAndMerge(AtlasExtrapolatorCfg(flags)))
+    kwargs.setdefault("Printer", result.addPublicTool(result.popToolsAndMerge(MuonEDMPrinterToolCfg(flags))) )
+    kwargs.setdefault("ParticleCaloExtensionTool",
+                        result.popToolsAndMerge(ParticleCaloExtensionToolCfg(flags, 
+                                                                             name='MuonParticleCaloExtensionTool')))
 
-    muon_ext_tool = CompFactory.Muon.MuonSystemExtensionTool("MuonSystemExtensionTool",
-                                                             ParticleCaloExtensionTool=particle_calo_extension_tool,
-                                                             Extrapolator=atlas_extrapolator)
+    muon_ext_tool = CompFactory.Muon.MuonSystemExtensionTool("MuonSystemExtensionTool", **kwargs)
     result.setPrivateTools(muon_ext_tool)
     return result
 
@@ -1253,4 +1253,4 @@ def MuonSegmentTagToolCfg(flags, name="MuonSegmentTagTool", **kwargs):
 
     result.setPrivateTools(
         CompFactory.MuonCombined.MuonSegmentTagTool(name, **kwargs))
-    return result
+    return result
\ No newline at end of file
diff --git a/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.cxx b/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.cxx
index 841f8d44f622608f4d041aaa4029498670149677..f2a83736736f773cf2add1b5a04165c9b8e35e2e 100644
--- a/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.cxx
+++ b/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.cxx
@@ -193,6 +193,8 @@ namespace MuonCombined {
                                                         << " layerDataVec size" << candidate->layerDataVec.size() << " hits size"
                                                         << candidate->hits.size());
 
+            float beta = candidate->betaFitResult.beta;
+        
             // loop over layers and perform segment finding, collect segments per layer
             for (const auto& layerData : candidate->layerDataVec) {
                 // store segments in layer
@@ -201,7 +203,7 @@ namespace MuonCombined {
                 // loop over maxima
                 for (const auto& maximumData : layerData.maximumDataVec) {
                     // find segments for intersection
-                    findSegments(layerData.intersection, *maximumData, segments, m_muonPRDSelectionToolStau, m_segmentMaker);
+                    findSegments(layerData.intersection, *maximumData, segments, m_muonPRDSelectionToolStau, m_segmentMaker, beta);
                 }
 
                 // skip if no segment were found
@@ -1136,14 +1138,17 @@ namespace MuonCombined {
     void MuonStauRecoTool::findSegments(const Muon::MuonSystemExtension::Intersection& intersection, MaximumData& maximumData,
                                         std::vector<std::shared_ptr<const Muon::MuonSegment>>& segments,
                                         const ToolHandle<Muon::IMuonPRDSelectionTool>& muonPRDSelectionTool,
-                                        const ToolHandle<Muon::IMuonSegmentMaker>& segmentMaker) const {
+                                        const ToolHandle<Muon::IMuonSegmentMaker>& segmentMaker,
+                                        float beta) const {
+
         const MuonHough::MuonLayerHough::Maximum& maximum = *maximumData.maximum;
         const std::vector<std::shared_ptr<const Muon::MuonClusterOnTrack>>& phiClusterOnTracks = maximumData.phiClusterOnTracks;
 
         // lambda to handle calibration and selection of MDTs
         auto handleMdt = [intersection, muonPRDSelectionTool](const Muon::MdtPrepData& prd,
-                                                              std::vector<const Muon::MdtDriftCircleOnTrack*>& mdts) {
-            const Muon::MdtDriftCircleOnTrack* mdt = muonPRDSelectionTool->calibrateAndSelect(intersection, prd);
+                                                             std::vector<const Muon::MdtDriftCircleOnTrack*>& mdts,
+                                                              float beta) {
+            const Muon::MdtDriftCircleOnTrack* mdt = muonPRDSelectionTool->calibrateAndSelect(intersection, prd, beta);
             if (mdt) mdts.push_back(mdt);
         };
 
@@ -1175,7 +1180,7 @@ namespace MuonCombined {
             } else if ((*hit)->prd) {
                 Identifier id = (*hit)->prd->identify();
                 if (m_idHelperSvc->isMdt(id))
-                    handleMdt(static_cast<const Muon::MdtPrepData&>(*(*hit)->prd), mdts);
+                    handleMdt(static_cast<const Muon::MdtPrepData&>(*(*hit)->prd), mdts, beta);
                 else
                     handleCluster(static_cast<const Muon::MuonCluster&>(*(*hit)->prd), clusters);
             }
@@ -1194,7 +1199,7 @@ namespace MuonCombined {
             // run segment finder
             std::unique_ptr<Trk::SegmentCollection> segColl(new Trk::SegmentCollection(SG::VIEW_ELEMENTS));
             segmentMaker->find(intersection.trackParameters->position(), intersection.trackParameters->momentum(), mdts, clusters,
-                               !clusters.empty(), segColl.get(), intersection.trackParameters->momentum().mag());
+                               !clusters.empty(), segColl.get(), intersection.trackParameters->momentum().mag(), 0, beta);
             if (segColl) {
                 Trk::SegmentCollection::iterator sit = segColl->begin();
                 Trk::SegmentCollection::iterator sit_end = segColl->end();
@@ -1344,7 +1349,9 @@ namespace MuonCombined {
                              : intersection.trackParameters->parameters()[Trk::locX];
 
         float z = intersection.trackParameters->position().z();
-        float errx = Amg::error(*intersection.trackParameters->covariance(), Trk::locX);
+        float errx = intersection.trackParameters->covariance() ?
+                     Amg::error(*intersection.trackParameters->covariance(), Trk::locX) : 0.;
+
         float x = barrelLike ? z : r;
         float y = barrelLike ? r : z;
         float theta = std::atan2(y, x);
@@ -1356,8 +1363,8 @@ namespace MuonCombined {
         // lambda to handle calibration and selection of clusters
         auto handleCluster = [intersection, this](const Muon::MuonCluster& prd,
                                                   std::vector<std::shared_ptr<const Muon::MuonClusterOnTrack>>& clusters) {
-            const Muon::MuonClusterOnTrack* cluster = m_muonPRDSelectionTool->calibrateAndSelect(intersection, prd);
-            if (cluster) clusters.push_back(std::shared_ptr<const Muon::MuonClusterOnTrack>(cluster));
+            std::unique_ptr<const Muon::MuonClusterOnTrack> cluster{m_muonPRDSelectionTool->calibrateAndSelect(intersection, prd)};
+            if (cluster) clusters.push_back(std::move(cluster));
         };
 
         // loop over maxima and associate phi hits with the extrapolation, should optimize this but calculating the residual with the phi
@@ -1373,7 +1380,8 @@ namespace MuonCombined {
                     Identifier id = hit->tgc->phiCluster.hitList.front()->identify();
                     if (m_idHelperSvc->layerIndex(id) != intersection.layerSurface.layerIndex) continue;
                     for (const Muon::MuonCluster* prd : hit->tgc->phiCluster.hitList) handleCluster(*prd, phiClusterOnTracks);
-                } else if (hit->prd && !(hit->prd->type(Trk::PrepRawDataType::sTgcPrepData) || hit->prd->type(Trk::PrepRawDataType::MMPrepData))) {
+                } else if (hit->prd && !(hit->prd->type(Trk::PrepRawDataType::sTgcPrepData) || 
+                                         hit->prd->type(Trk::PrepRawDataType::MMPrepData))) {
                     const Identifier id = hit->prd->identify();
                     if (m_idHelperSvc->layerIndex(id) != intersection.layerSurface.layerIndex) continue;
                     handleCluster(static_cast<const Muon::MuonCluster&>(*hit->prd), phiClusterOnTracks);
@@ -1386,20 +1394,22 @@ namespace MuonCombined {
                                                 << " angle " << theta);
 
         // loop over maxima and associate them to the extrapolation
-        Muon::MuonLayerHoughTool::MaximumVec::const_iterator mit = maxVec.begin();
-        Muon::MuonLayerHoughTool::MaximumVec::const_iterator mit_end = maxVec.end();
-        for (; mit != mit_end; ++mit) {
-            const MuonHough::MuonLayerHough::Maximum& maximum = **mit;
-            if (std::find_if(maximum.hits.begin(),maximum.hits.end(),[](const std::shared_ptr<MuonHough::Hit>& hit){
-                return hit->prd && (hit->prd->type(Trk::PrepRawDataType::sTgcPrepData) || hit->prd->type(Trk::PrepRawDataType::MMPrepData));
-            }) != maximum.hits.end()) continue;
+        for (const auto& mit : maxVec) {
+            const MuonHough::MuonLayerHough::Maximum& maximum = *mit;
+            if (std::find_if(maximum.hits.begin(),maximum.hits.end(),
+                             [](const std::shared_ptr<MuonHough::Hit>& hit){
+                                return hit->prd && (hit->prd->type(Trk::PrepRawDataType::sTgcPrepData) || 
+                                                    hit->prd->type(Trk::PrepRawDataType::MMPrepData));
+                            }) != maximum.hits.end()) continue;
             float residual = maximum.pos - x;
             float residualTheta = maximum.theta - theta;
             float refPos = (maximum.hough != nullptr) ? maximum.hough->m_descriptor.referencePosition : 0;
             float maxwidth = (maximum.binposmax - maximum.binposmin);
-            if (maximum.hough) maxwidth *= maximum.hough->m_binsize;
-
-            float pull = residual / std::sqrt(errx * errx + maxwidth * maxwidth / 12.);
+            if (maximum.hough){
+                maxwidth *= maximum.hough->m_binsize;
+            }
+            const float pullUncert = std::sqrt(errx * errx + maxwidth * maxwidth / 12.);  
+            float pull = residual / (pullUncert > std::numeric_limits<float>::epsilon() ? pullUncert : 1.) ;
 
             ATH_MSG_DEBUG("   Hough maximum " << maximum.max << " position (" << refPos << "," << maximum.pos << ") residual " << residual
                                               << " pull " << pull << " angle " << maximum.theta << " residual " << residualTheta);
diff --git a/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.h b/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.h
index 8db6d66ca524af5f81fea3a6cef9f6f9391fd199..e7ca65fdd5166e2a929ffe688dd5f255b1d03148 100644
--- a/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.h
+++ b/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.h
@@ -164,7 +164,9 @@ namespace MuonCombined {
         void findSegments(const Muon::MuonSystemExtension::Intersection& intersection, MaximumData& maximumData,
                           std::vector<std::shared_ptr<const Muon::MuonSegment>>& t0fittedSegments,
                           const ToolHandle<Muon::IMuonPRDSelectionTool>& muonPRDSelectionTool,
-                          const ToolHandle<Muon::IMuonSegmentMaker>& segmentMaker) const;
+                          const ToolHandle<Muon::IMuonSegmentMaker>& segmentMaker,
+                          float beta = 1.) const;
+
 
         /** associate Hough maxima and associate time measurements */
         bool extractTimeMeasurements(const EventContext& ctx, const Muon::MuonSystemExtension& muonSystemExtension,
diff --git a/Reconstruction/MuonIdentification/MuonSegmentTaggers/MuonSegmentTaggerTools/src/MuTagMatchingTool.cxx b/Reconstruction/MuonIdentification/MuonSegmentTaggers/MuonSegmentTaggerTools/src/MuTagMatchingTool.cxx
index 0bb3ce6ff8db93d729f48c057abd3329edfcd947..cc924e332a11a2d278ecbb896e79c684affbca64 100644
--- a/Reconstruction/MuonIdentification/MuonSegmentTaggers/MuonSegmentTaggerTools/src/MuTagMatchingTool.cxx
+++ b/Reconstruction/MuonIdentification/MuonSegmentTaggers/MuonSegmentTaggerTools/src/MuTagMatchingTool.cxx
@@ -773,6 +773,8 @@ void MuTagMatchingTool::calculateLocalAngleErrors(const Muon::MuonSegment& segme
 
 void MuTagMatchingTool::calculateLocalAngleErrors(const Trk::AtaPlane& exTrack, double& angleXZerror, double& angleYZerror,
                                                   double& covLocYYZ) const {
+                                                    
+    if (!exTrack.covariance()) return;
     // Parameters are described as Trk::LocX, Trk::locY, Trk::phi, Trk::theta
     // So the errormatrix of the track 'localErrorMatrix' still holds global angle representation!!!!
     // retrieve Jabcobian to transform the global errors err_phi,err_theta to local errors err_alphaXZ, err_alphaYZ
diff --git a/Tools/PROCTools/data/q442_AOD_content.ref b/Tools/PROCTools/data/q442_AOD_content.ref
index 72a7b5060db53d1b1f507203e48fd803b301caab..002dd38806d9e2591fc203c2bca8b0982bff36f2 100644
--- a/Tools/PROCTools/data/q442_AOD_content.ref
+++ b/Tools/PROCTools/data/q442_AOD_content.ref
@@ -79,6 +79,31 @@ CombinedMuonsLRTTrackParticles
 CombinedMuonsLRTTrackParticlesAux.
 CombinedStauTrackParticles
 CombinedStauTrackParticlesAux.
+CombinedStauTrackParticlesAuxDyn.TRTdEdx
+CombinedStauTrackParticlesAuxDyn.TRTdEdxUsedHits
+CombinedStauTrackParticlesAuxDyn.alignEffectChId
+CombinedStauTrackParticlesAuxDyn.alignEffectDeltaAngle
+CombinedStauTrackParticlesAuxDyn.alignEffectDeltaTrans
+CombinedStauTrackParticlesAuxDyn.alignEffectSigmaDeltaAngle
+CombinedStauTrackParticlesAuxDyn.alignEffectSigmaDeltaTrans
+CombinedStauTrackParticlesAuxDyn.deltaphi_0
+CombinedStauTrackParticlesAuxDyn.deltaphi_1
+CombinedStauTrackParticlesAuxDyn.deltatheta_0
+CombinedStauTrackParticlesAuxDyn.deltatheta_1
+CombinedStauTrackParticlesAuxDyn.eProbabilityNN
+CombinedStauTrackParticlesAuxDyn.parameterPX
+CombinedStauTrackParticlesAuxDyn.parameterPY
+CombinedStauTrackParticlesAuxDyn.parameterPZ
+CombinedStauTrackParticlesAuxDyn.parameterPosition
+CombinedStauTrackParticlesAuxDyn.parameterX
+CombinedStauTrackParticlesAuxDyn.parameterY
+CombinedStauTrackParticlesAuxDyn.parameterZ
+CombinedStauTrackParticlesAuxDyn.sigmadeltaphi_0
+CombinedStauTrackParticlesAuxDyn.sigmadeltaphi_1
+CombinedStauTrackParticlesAuxDyn.sigmadeltatheta_0
+CombinedStauTrackParticlesAuxDyn.sigmadeltatheta_1
+CombinedStauTrackParticlesAuxDyn.trackLink
+CombinedStauTrackParticlesAuxDyn.trackParameterCovarianceMatrices
 ConditionsRun
 DataHeader
 DataHeaderForm
@@ -980,10 +1005,51 @@ PrimaryVerticesAuxDyn.sumPt2
 RunNumber
 SlowMuons
 SlowMuonsAux.
+SlowMuonsAuxDyn.hitEnergy
+SlowMuonsAuxDyn.hitError
+SlowMuonsAuxDyn.hitIdentifier
+SlowMuonsAuxDyn.hitPositionX
+SlowMuonsAuxDyn.hitPositionY
+SlowMuonsAuxDyn.hitPositionZ
+SlowMuonsAuxDyn.hitPropagationTime
+SlowMuonsAuxDyn.hitShift
+SlowMuonsAuxDyn.hitTOF
+SlowMuonsAuxDyn.hitTechnology
 StauSegments
 StauSegmentsAux.
+StauSegmentsAuxDyn.clusterTime
+StauSegmentsAuxDyn.clusterTimeError
+StauSegmentsAuxDyn.clusterTimeValid
+StauSegmentsAuxDyn.muonSegment
 Staus
 StausAux.
+StausAuxDyn.CT_EL_Type
+StausAuxDyn.CT_ET_FSRCandidateEnergy
+StausAuxDyn.ET_Core
+StausAuxDyn.ET_EMCore
+StausAuxDyn.ET_HECCore
+StausAuxDyn.ET_TileCore
+StausAuxDyn.FSR_CandidateEnergy
+StausAuxDyn.combinedTrackOutBoundsPrecisionHits
+StausAuxDyn.d0_staco
+StausAuxDyn.extendedClosePrecisionHits
+StausAuxDyn.extendedOutBoundsPrecisionHits
+StausAuxDyn.innerClosePrecisionHits
+StausAuxDyn.innerOutBoundsPrecisionHits
+StausAuxDyn.isEndcapGoodLayers
+StausAuxDyn.isSmallGoodSectors
+StausAuxDyn.middleClosePrecisionHits
+StausAuxDyn.middleOutBoundsPrecisionHits
+StausAuxDyn.nUnspoiledCscHits
+StausAuxDyn.numEnergyLossPerTrack
+StausAuxDyn.numberOfGoodPrecisionLayers
+StausAuxDyn.outerClosePrecisionHits
+StausAuxDyn.outerOutBoundsPrecisionHits
+StausAuxDyn.phi0_staco
+StausAuxDyn.qOverPErr_staco
+StausAuxDyn.qOverP_staco
+StausAuxDyn.theta_staco
+StausAuxDyn.z0_staco
 StreamAOD
 TauFinalPi0s
 TauFinalPi0sAux.
diff --git a/Tools/PROCTools/data/q449_AOD_digest.ref b/Tools/PROCTools/data/q449_AOD_digest.ref
index b1440b2836b6288f5b232d137acdffe3df3c7f46..01b831400863d26a72f363d3fee93a33987bd1a0 100644
--- a/Tools/PROCTools/data/q449_AOD_digest.ref
+++ b/Tools/PROCTools/data/q449_AOD_digest.ref
@@ -7,7 +7,7 @@
       431493  1096113771         294         409          77           6           0          10           0          10          12           0          12
       431493  1096114154         374         572          87           6           0           6           0           6          18           0          18
       431493  1096114168         132         151          17           4           0           2           0           2           6           0           6
-      431493  1096115980         390         519         130          10           1          20           0          20          22           0          22
+      431493  1096115980         390         519         130          10           0          20           0          20          22           0          22
       431493  1096116200         309         410          57           6           1           4           0           4          12           0          12
       431493  1096118782         294         357          95           7           1          12           0          12          23           0          23
       431493  1096119418         307         337          81           6           0           5           0           5           9           0           9
@@ -47,7 +47,7 @@
       431493  1096144036         268         343          65           5           0          12           0          12          14           0          14
       431493  1096145110         241         257          18           2           0           5           0           5           8           0           8
       431493  1096146639         283         422          77           6           0           3           0           3          18           0          18
-      431493  1096146735         321         378         105          10           2          21           0          21          26           0          26
+      431493  1096146735         321         378         105          10           1          21           0          21          26           0          26
       431493  1096147568         392         600          41           4           2           4           0           4          14           0          14
       431493  1096150983         272         320          22           5           0           2           0           2          13           0          13
       431493  1096152978         365         498          87           7           2           8           0           8          11           0          11
diff --git a/Tools/WorkflowTestRunner/python/References.py b/Tools/WorkflowTestRunner/python/References.py
index 8001884a75f44cd8944998716878926baec06819..61e92ad401c0d7eb278d2eff240ff83bcc2fbe74 100644
--- a/Tools/WorkflowTestRunner/python/References.py
+++ b/Tools/WorkflowTestRunner/python/References.py
@@ -22,10 +22,10 @@ references_map = {
     "d1759": "v4",
     "d1912": "v5",
     # Reco
-    "q442": "v9",
-    "q449": "v26",
-    "q452": "v10",
-    "q454": "v23",
+    "q442": "v10",
+    "q449": "v27",
+    "q452": "v11",
+    "q454": "v24",
     # Derivations
     "data_PHYS_Run2": "v3",
     "data_PHYS_Run3": "v3",