diff --git a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/HoughDataPerSec.h b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/HoughDataPerSec.h
new file mode 100644
index 0000000000000000000000000000000000000000..3dcd46c589b5faa452ccdcbf69e02b7cd83aa5af
--- /dev/null
+++ b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/HoughDataPerSec.h
@@ -0,0 +1,62 @@
+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MUONHOUGHPATTERNTOOLS_HOUGHDATAPERSEC_H
+#define MUONHOUGHPATTERNTOOLS_HOUGHDATAPERSEC_H
+
+#include "AthenaKernel/CLASS_DEF.h"
+#include "MuonLayerHough/MuonLayerHough.h"
+#include "MuonLayerHough/MuonPhiLayerHough.h"
+#include <map>
+#include <set>
+#include <vector>
+
+
+namespace Muon {
+  struct HoughDataPerSec {
+    typedef std::vector<MuonHough::Hit*>    HitVec;
+    typedef std::vector< HitVec >           RegionHitVec;
+    typedef std::vector<MuonHough::PhiHit*> PhiHitVec;
+    typedef std::vector< PhiHitVec >        RegionPhiHitVec;
+    typedef std::vector<MuonHough::MuonLayerHough::Maximum*>    MaximumVec;
+    typedef std::vector<MuonHough::MuonPhiLayerHough::Maximum*> PhiMaximumVec; 
+    typedef std::map<MuonHough::MuonLayerHough::Maximum*, MaximumVec > MaximumAssociationMap;
+    typedef std::vector< MaximumVec >       RegionMaximumVec;
+    typedef std::vector< PhiMaximumVec >    RegionPhiMaximumVec;
+
+    HoughDataPerSec() {
+      sector = -1;
+      hitVec.resize(MuonStationIndex::sectorLayerHashMax());
+      maxVec.resize(MuonStationIndex::sectorLayerHashMax());
+      phiHitVec.resize(MuonStationIndex::DetectorRegionIndexMax);
+      phiMaxVec.resize(MuonStationIndex::DetectorRegionIndexMax);
+      nlayersWithMaxima.resize(MuonStationIndex::DetectorRegionIndexMax);
+      nphilayersWithMaxima.resize(MuonStationIndex::DetectorRegionIndexMax);
+      nmaxHitsInRegion.resize(MuonStationIndex::DetectorRegionIndexMax);
+      nphimaxHitsInRegion.resize(MuonStationIndex::DetectorRegionIndexMax);
+    }
+    void cleanUp();
+    int                   sector;
+    RegionHitVec          hitVec;
+    RegionPhiHitVec       phiHitVec;
+    RegionMaximumVec      maxVec;
+    RegionPhiMaximumVec   phiMaxVec;
+    std::vector<int>      nlayersWithMaxima;
+    std::vector<int>      nphilayersWithMaxima;
+    std::vector<int>      nmaxHitsInRegion;
+    std::vector<int>      nphimaxHitsInRegion;
+    MaximumAssociationMap maxAssociationMap; // stores association of a given maximium with other maxima in neighbouring sectors
+    std::set<MuonHough::MuonLayerHough::Maximum*> associatedToOtherSector; // used to flagged maxima that were already associated to another sector
+    // returns the number of phi and eta hits in the region with most eta hits
+    // regions with phi hits are always prefered over regions without
+    int maxEtaHits() const{
+      return std::max( nmaxHitsInRegion[0], std::max( nmaxHitsInRegion[1], nmaxHitsInRegion[2] ) );
+    }
+  };
+}
+
+CLASS_DEF(Muon::HoughDataPerSec, 163257499, 1)
+CLASS_DEF(std::vector<Muon::HoughDataPerSec>, 118215228, 1)
+
+#endif
\ No newline at end of file
diff --git a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/MuonHoughPatternFinderTool.h b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/MuonHoughPatternFinderTool.h
index 48c97808e849702aeab28522a0b91794f5d80686..4c7bcd39c6491c1afbb60d93a7564112aefd2348 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/MuonHoughPatternFinderTool.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/MuonHoughPatternFinderTool.h
@@ -58,11 +58,12 @@ namespace Muon {
     virtual StatusCode finalize();
 
     /** find patterns for a give set of MuonPrepData collections + optionally CSC segment combinations */
-    MuonPatternCombinationCollection* find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
-					    const std::vector<const CscPrepDataCollection*>& cscCols,  
-					    const std::vector<const TgcPrepDataCollection*>& tgcCols,  
-					    const std::vector<const RpcPrepDataCollection*>& rpcCols,  
-					    const MuonSegmentCombinationCollection* cscSegmentCombis ) const;
+    std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<std::vector<Muon::HoughDataPerSec>>>
+    find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
+          const std::vector<const CscPrepDataCollection*>& cscCols,  
+          const std::vector<const TgcPrepDataCollection*>& tgcCols,  
+          const std::vector<const RpcPrepDataCollection*>& rpcCols,  
+          const MuonSegmentCombinationCollection* cscSegmentCombis ) const;
 
   private:
 
diff --git a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/MuonLayerHoughTool.h b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/MuonLayerHoughTool.h
index 2eeeb29b503e6eb832df6ae05375668345a5363f..d8afe07d870205a7b857fad16b4a0ff70a165734 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/MuonLayerHoughTool.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/MuonHoughPatternTools/MuonLayerHoughTool.h
@@ -26,6 +26,7 @@
 #include "MuonLayerHough/MuonLayerHough.h"
 #include "MuonLayerHough/MuonPhiLayerHough.h"
 #include "MuonLayerHough/MuonLayerHoughSelector.h"
+#include "MuonHoughPatternTools/HoughDataPerSec.h"
 
 #include "TMath.h"
 #include <set>
@@ -52,8 +53,6 @@ namespace MuonHough {
   class HitDebugInfo;
 }
 
-
-
 static const InterfaceID IID_MuonLayerHoughTool("Muon::MuonLayerHoughTool",1,0);
 
 namespace Muon {
@@ -76,50 +75,20 @@ namespace Muon {
     typedef std::vector<CollectionsPerSector>       CollectionsPerSectorVec;
     typedef CollectionsPerSectorVec::const_iterator CollectionsPerSectorCit;
 
-    typedef std::vector<MuonHough::Hit*>    HitVec;
-    typedef std::vector< HitVec >           RegionHitVec;
-    typedef std::vector<MuonHough::PhiHit*> PhiHitVec;
-    typedef std::vector< PhiHitVec >        RegionPhiHitVec;
-    typedef std::vector<MuonHough::MuonLayerHough::Maximum*>    MaximumVec;
-    typedef std::vector<MuonHough::MuonPhiLayerHough::Maximum*> PhiMaximumVec; 
-    typedef std::map<MuonHough::MuonLayerHough::Maximum*, MaximumVec > MaximumAssociationMap;
-
-    typedef std::vector< MaximumVec >       RegionMaximumVec;
-    typedef std::vector< PhiMaximumVec >    RegionPhiMaximumVec;
-
-    struct HoughDataPerSector {
-      HoughDataPerSector() {
-	sector = -1;
-	hitVec.resize(MuonStationIndex::sectorLayerHashMax());
-	maxVec.resize(MuonStationIndex::sectorLayerHashMax());
-	phiHitVec.resize(MuonStationIndex::DetectorRegionIndexMax);
-	phiMaxVec.resize(MuonStationIndex::DetectorRegionIndexMax);
-	nlayersWithMaxima.resize(MuonStationIndex::DetectorRegionIndexMax);
-	nphilayersWithMaxima.resize(MuonStationIndex::DetectorRegionIndexMax);
-	nmaxHitsInRegion.resize(MuonStationIndex::DetectorRegionIndexMax);
-	nphimaxHitsInRegion.resize(MuonStationIndex::DetectorRegionIndexMax);
-      }
-      void cleanUp();
-      int                   sector;
-      RegionHitVec          hitVec;
-      RegionPhiHitVec       phiHitVec;
-      RegionMaximumVec      maxVec;
-      RegionPhiMaximumVec   phiMaxVec;
-      std::vector<int>      nlayersWithMaxima;
-      std::vector<int>      nphilayersWithMaxima;
-      std::vector<int>      nmaxHitsInRegion;
-      std::vector<int>      nphimaxHitsInRegion;
-      MaximumAssociationMap maxAssociationMap; // stores association of a given maximium with other maxima in neighbouring sectors
-      std::set<MuonHough::MuonLayerHough::Maximum*> associatedToOtherSector; // used to flagged maxima that were already associated to another sector
-      // returns the number of phi and eta hits in the region with most eta hits
-      // regions with phi hits are always prefered over regions without
-      int maxEtaHits() const{
-	return std::max( nmaxHitsInRegion[0], std::max( nmaxHitsInRegion[1], nmaxHitsInRegion[2] ) );
-      }
-    };
-
+    typedef HoughDataPerSec::HitVec HitVec;
+    typedef HoughDataPerSec::RegionHitVec RegionHitVec;
+    typedef HoughDataPerSec::PhiHitVec PhiHitVec;
+    typedef HoughDataPerSec::RegionPhiHitVec RegionPhiHitVec;
+    typedef HoughDataPerSec::MaximumVec MaximumVec;
+    typedef HoughDataPerSec::PhiMaximumVec PhiMaximumVec; 
+    typedef HoughDataPerSec::MaximumAssociationMap MaximumAssociationMap;
+    typedef HoughDataPerSec::RegionMaximumVec RegionMaximumVec;
+    typedef HoughDataPerSec::RegionPhiMaximumVec RegionPhiMaximumVec;
+
+    typedef HoughDataPerSec                       HoughDataPerSector;
     typedef std::vector<HoughDataPerSector>       HoughDataPerSectorVec;
     typedef HoughDataPerSectorVec::const_iterator HoughDataPerSectorCit;
+
     
     
     class Road {
@@ -153,11 +122,12 @@ namespace Muon {
     /** @brief access to tool interface */
     static const InterfaceID& interfaceID() { return IID_MuonLayerHoughTool; }
 
-    StatusCode initialize();
+    virtual StatusCode initialize() override;
 
-    StatusCode finalize();
+    virtual StatusCode finalize() override;
     
-    MuonPatternCombinationCollection* analyse( const MdtPrepDataContainer*  mdtCont,  
+    std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> analyse(
+                 const MdtPrepDataContainer*  mdtCont,  
 					       const CscPrepDataContainer*  cscCols,  
 					       const TgcPrepDataContainer*  tgcCont,  
 					       const RpcPrepDataContainer*  rpcCont,
@@ -165,11 +135,12 @@ namespace Muon {
 					       const MMPrepDataContainer*  mmCont ) const;
 
     /** find patterns for a give set of MuonPrepData collections + optionally CSC segment combinations */
-    MuonPatternCombinationCollection* find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
-					    const std::vector<const CscPrepDataCollection*>& cscCols,  
-					    const std::vector<const TgcPrepDataCollection*>& tgcCols,  
-					    const std::vector<const RpcPrepDataCollection*>& rpcCols,  
-					    const MuonSegmentCombinationCollection* ) const;
+    virtual std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>>
+    find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
+          const std::vector<const CscPrepDataCollection*>& cscCols,  
+          const std::vector<const TgcPrepDataCollection*>& tgcCols,  
+          const std::vector<const RpcPrepDataCollection*>& rpcCols,  
+          const MuonSegmentCombinationCollection* ) const override;
     void reset() const;
 
     void getSectors( const Amg::Vector3D& pos, std::vector<int>& sectors ) const;
@@ -182,17 +153,28 @@ namespace Muon {
 
     int sublay( const Identifier& id, float z = 0 ) const; // the z value is only used for the tgcs
 
-    const MuonHough::MuonDetectorHough& houghTransforms() const { return m_detectorHoughTransforms; }
-    const HoughDataPerSectorVec&        houghData() const { return m_houghDataPerSectorVec; }
+    //const MuonHough::MuonDetectorHough& houghTransforms() const { return m_detectorHoughTransforms; }
+    //const HoughDataPerSectorVec&        houghData() const { return m_houghDataPerSectorVec; }
 
     /**  incident service handle for EndEvent */
-    void handle(const Incident& inc);// maybe in the future clear per event
+    virtual void handle(const Incident& inc) override;// maybe in the future clear per event
 
   private:
 
-    MuonPatternCombinationCollection* analyse() const;
+    struct State {
+      MaximumVec seedMaxima;
+      MuonHough::MuonDetectorHough detectorHoughTransforms;
+      std::unique_ptr<HoughDataPerSectorVec> houghDataPerSectorVec { std::make_unique<HoughDataPerSectorVec>() };
+      std::set<Identifier> truthHits;
+      std::set<Identifier> foundTruthHits;
+      std::set<Identifier> outputTruthHits;
+      std::vector<TgcHitClusteringObj*> tgcClusteringObjs;
+    };
+  
+    std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> analyse(State& state) const;
 
-    void fillHitsPerSector(  const CollectionsPerSector& hashes,
+    void fillHitsPerSector(  State& state, 
+           const CollectionsPerSector& hashes,
 			     const MdtPrepDataContainer*  mdtCont,  
 			     const CscPrepDataContainer*  /*cscCont*/,  
 			     const TgcPrepDataContainer*  tgcCont,  
@@ -201,16 +183,18 @@ namespace Muon {
 			     const MMPrepDataContainer*   mmCont,
 			     HoughDataPerSector& houghData ) const;
 
-    void fill( const MdtPrepDataCollection& mdts, HitVec& hits ) const;
-    void fill( const TgcPrepDataCollection& tgcs, HitVec& hits, PhiHitVec& phiHits, int sector ) const;
-    void fill( const RpcPrepDataCollection& rpcs, HitVec& hits, PhiHitVec& phiHits ) const;
-    void fill( const MMPrepDataCollection& mdts, HitVec& hits ) const;
-    void fill( const sTgcPrepDataCollection& stgcs, HitVec& hits, PhiHitVec& phiHits, int sector ) const;
+    void fill( State& state, const MdtPrepDataCollection& mdts, HitVec& hits ) const;
+    void fill( State& state, const TgcPrepDataCollection& tgcs, HitVec& hits, PhiHitVec& phiHits, int sector ) const;
+    void fill( State& state, const RpcPrepDataCollection& rpcs, HitVec& hits, PhiHitVec& phiHits ) const;
+    void fill( State& state, const MMPrepDataCollection& mdts, HitVec& hits ) const;
+    void fill( State& state, const sTgcPrepDataCollection& stgcs, HitVec& hits, PhiHitVec& phiHits, int sector ) const;
 
-    bool findMaxima( MuonHough::MuonLayerHough& hough,
+    bool findMaxima( State& state, 
+         MuonHough::MuonLayerHough& hough,
 		     HitVec& hits, 
 		     MaximumVec& maxima ) const ;
-    bool findMaxima( MuonHough::MuonPhiLayerHough& hough,
+    bool findMaxima( State& state,
+         MuonHough::MuonPhiLayerHough& hough,
 		     PhiHitVec& hits, 
 		     PhiMaximumVec& maxima,
 		     int sector ) const;
@@ -222,7 +206,7 @@ namespace Muon {
 
     void associateMaximaInNeighbouringSectors( HoughDataPerSector& houghData, std::vector<HoughDataPerSector>& houghDataPerSectorVec ) const;
 
-    void extendSeed( Road& road, HoughDataPerSector& sectorData ) const; // const;
+    void extendSeed(State& state, Road& road, HoughDataPerSector& sectorData ) const; // const;
     void associatePhiMaxima( Road& road, PhiMaximumVec& phiMaxima ) const;
 
     double combinedPeakheight( double ph,  double ph1,  double ph2, double phn, double rot, int layer, int /*region*/ ) const;
@@ -231,7 +215,7 @@ namespace Muon {
     void createPatternCombinations( std::vector< MaximumVec >& maxima,
 				    MuonPatternCombinationCollection& patternCombis ) const;
 
-    void createPatternCombinations( std::map< MuonHough::MuonPhiLayerHough::Maximum*, MaximumVec >& phiEtaAssociations,
+    void createPatternCombinations(State& state, std::map< MuonHough::MuonPhiLayerHough::Maximum*, MaximumVec >& phiEtaAssociations,
 				    MuonPatternCombinationCollection& patternCombis ) const;
 
     void fillNtuple( HoughDataPerSectorVec& houghDataPerSectorVec ) const;
@@ -239,23 +223,19 @@ namespace Muon {
     void insertHash( const IdentifierHash& hash, const Identifier& id );
     void insertHash( int sector, const IdentifierHash& hash, const Identifier& id );
 
-    void matchTruth( const PRD_MultiTruthCollection& truthCol, const Identifier& id, MuonHough::HitDebugInfo& debug ) const;
+    void matchTruth( State& state, const PRD_MultiTruthCollection& truthCol, const Identifier& id, MuonHough::HitDebugInfo& debug ) const;
     void initializeSectorMapping();
     void getTruth() const;
     void printTruthSummary( std::set<Identifier>& truth, std::set<Identifier>& found ) const;
 
-    void buildRoads( std::vector<Road>& roads ) const;
+    void buildRoads(State& state, std::vector<Road>& roads ) const;
     void mergePhiMaxima( Road& road ) const;
 
-    mutable MaximumVec m_seedMaxima;
     bool m_useSeeds;
 
-    mutable MuonHough::MuonDetectorHough m_detectorHoughTransforms;
-    mutable HoughDataPerSectorVec        m_houghDataPerSectorVec;
-
     ToolHandle<MuonIdHelperTool> m_idHelper;
     ToolHandle<MuonEDMPrinterTool> m_printer;
-    mutable ToolHandle<Muon::IMuonTruthSummaryTool>         m_truthSummaryTool;
+    mutable ToolHandle<Muon::IMuonTruthSummaryTool> m_truthSummaryTool;
     const MuonGM::MuonDetectorManager* m_detMgr;
 
     std::vector<MuonHough::MuonLayerHoughSelector> m_selectors;
@@ -268,10 +248,6 @@ namespace Muon {
     SG::ReadHandleKeyArray< PRD_MultiTruthCollection >       m_truthNames; 
     SG::ReadHandleKey<xAOD::TruthParticleContainer>       m_MuonTruthParticlesKey;
     SG::ReadHandleKey<xAOD::MuonSegmentContainer>       m_MuonTruthSegmentsKey;
-
-    mutable std::set<Identifier>            m_truthHits;
-    mutable std::set<Identifier>            m_foundTruthHits;
-    mutable std::set<Identifier>            m_outputTruthHits;
     
     bool m_useRpcTimeVeto;
     bool m_requireTriggerConfirmationNSW;
@@ -284,7 +260,6 @@ namespace Muon {
 
     unsigned int m_ntechnologies;
     CollectionsPerSectorVec m_collectionsPerSector;
-    mutable std::vector<TgcHitClusteringObj*> m_tgcClusteringObjs;
 
     ServiceHandle< IIncidentSvc >  m_incidentSvc;
     MuonSectorMapping              m_sectorMapping;
diff --git a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/HoughDataPerSec.cxx b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/HoughDataPerSec.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..fea940d7ebdff6a9d0d182cb9cd3c07bf73ba929
--- /dev/null
+++ b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/HoughDataPerSec.cxx
@@ -0,0 +1,27 @@
+/*
+  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MuonHoughPatternTools/HoughDataPerSec.h"
+
+namespace Muon {
+
+void HoughDataPerSec::cleanUp() {
+  for(RegionHitVec::iterator it=hitVec.begin();it!=hitVec.end();++it)
+    for( HitVec::iterator it2=it->begin();it2!=it->end();++it2 ) delete *it2;
+  hitVec.clear();
+
+  for(RegionPhiHitVec::iterator it=phiHitVec.begin();it!=phiHitVec.end();++it)
+    for( PhiHitVec::iterator it2=it->begin();it2!=it->end();++it2 ) delete *it2;
+  phiHitVec.clear();
+
+  for(RegionMaximumVec::iterator it=maxVec.begin();it!=maxVec.end();++it)
+    for( MaximumVec::iterator it2=it->begin();it2!=it->end();++it2 ) delete *it2;
+  maxVec.clear();
+
+  for(RegionPhiMaximumVec::iterator it=phiMaxVec.begin();it!=phiMaxVec.end();++it)
+    for( PhiMaximumVec::iterator it2=it->begin();it2!=it->end();++it2 ) delete *it2;
+  phiMaxVec.clear();
+}
+
+}
\ No newline at end of file
diff --git a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonHoughPatternFinderTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonHoughPatternFinderTool.cxx
index 7c5d73157c72bd82bd4ce014fbfdc3cfdbbfb9c1..40eba039d280e2533d6e025cf4151f93f028ca9a 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonHoughPatternFinderTool.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonHoughPatternFinderTool.cxx
@@ -156,18 +156,19 @@ namespace Muon {
     return StatusCode::SUCCESS; 
   }
 
-  MuonPatternCombinationCollection* MuonHoughPatternFinderTool::find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
-								      const std::vector<const CscPrepDataCollection*>& cscCols,  
-								      const std::vector<const TgcPrepDataCollection*>& tgcCols,  
-								      const std::vector<const RpcPrepDataCollection*>& rpcCols,  
-								      const MuonSegmentCombinationCollection* cscSegmentCombis ) const {
+  std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<std::vector<HoughDataPerSec>>>
+  MuonHoughPatternFinderTool::find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
+                                    const std::vector<const CscPrepDataCollection*>& cscCols,  
+                                    const std::vector<const TgcPrepDataCollection*>& tgcCols,  
+                                    const std::vector<const RpcPrepDataCollection*>& rpcCols,  
+                                    const MuonSegmentCombinationCollection* cscSegmentCombis ) const {
     // read event_data:
     const MuonHoughHitContainer* hitcontainer = getAllHits( mdtCols, cscCols, tgcCols, rpcCols, cscSegmentCombis );
 
     // analyse data
-    MuonPatternCombinationCollection* patCombiCol = 0;
+    std::unique_ptr<MuonPatternCombinationCollection> patCombiCol;
     if( hitcontainer ) {
-      patCombiCol = analyse( *hitcontainer );
+      patCombiCol.reset(analyse( *hitcontainer ));
     }else{
       ATH_MSG_INFO (" No hit container created! ");
     }
@@ -178,7 +179,7 @@ namespace Muon {
     // ensure we always output a collection
     if( !patCombiCol ){
       ATH_MSG_DEBUG (" NO pattern combinations found, creating empty collection ");
-      patCombiCol = new MuonPatternCombinationCollection();
+      patCombiCol.reset(new MuonPatternCombinationCollection());
     }
 
     // summary
@@ -194,7 +195,7 @@ namespace Muon {
     ATH_MSG_VERBOSE ("execute(end) ");
 
     // return result
-    return patCombiCol;
+    return {std::move(patCombiCol), nullptr};
   } 
 
   void MuonHoughPatternFinderTool::cleanUp() const {
diff --git a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughAlg.cxx b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughAlg.cxx
index a7f11ef20ea499b64722afc09d11903a79893a68..e22b98edbeb074747bf2620f869bf95de294d690 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughAlg.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughAlg.cxx
@@ -54,6 +54,7 @@ StatusCode MuonLayerHoughAlg::initialize()
   ATH_CHECK( m_keysTgc.initialize());
   ATH_CHECK( m_keyMM.initialize()  );
   ATH_CHECK( m_combis.initialize() );
+  ATH_CHECK( m_houghDataPerSectorVecKey.initialize() );
 
   return StatusCode::SUCCESS; 
 }
@@ -69,7 +70,8 @@ StatusCode MuonLayerHoughAlg::execute()
   const Muon::MMPrepDataContainer* mmPrds =GetObject(m_keyMM);
 
   ATH_MSG_VERBOSE("calling layer tool ");
-  std::unique_ptr<MuonPatternCombinationCollection> combis(m_layerTool->analyse(mdtPrds,cscPrds,tgcPrds,rpcPrds,stgcPrds,mmPrds));
+  auto [combis, houghDataPerSectorVec] = m_layerTool->analyse(mdtPrds,cscPrds,tgcPrds,rpcPrds,stgcPrds,mmPrds);
+
   if( combis ){
 
     SG::WriteHandle<MuonPatternCombinationCollection> Handle(m_combis);
@@ -83,6 +85,12 @@ StatusCode MuonLayerHoughAlg::execute()
     }
   }
 
+  // write hough data to SG
+  if (houghDataPerSectorVec) {
+    SG::WriteHandle<Muon::MuonLayerHoughTool::HoughDataPerSectorVec> handle {m_houghDataPerSectorVecKey};
+    ATH_CHECK(handle.record(std::move(houghDataPerSectorVec)));
+  }
+
   return StatusCode::SUCCESS;
 } // execute
 
diff --git a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughAlg.h b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughAlg.h
index 10f176371ee8a9ded6fddcc06d92636e17b78352..77e41d54a15c0baa05d5033f622bd723121a6599 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughAlg.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughAlg.h
@@ -36,6 +36,8 @@ class MuonLayerHoughAlg : public AthAlgorithm
   SG::ReadHandleKey<Muon::MMPrepDataContainer>    m_keyMM;
 
   SG::WriteHandleKey<MuonPatternCombinationCollection> m_combis;
+  SG::WriteHandleKey<Muon::MuonLayerHoughTool::HoughDataPerSectorVec> m_houghDataPerSectorVecKey {this, 
+    "Key_MuonLayerHoughToolHoughDataPerSectorVec", "HoughDataPerSectorVec", "HoughDataPerSectorVec key"};
   ToolHandle<Muon::MuonEDMPrinterTool> m_printer;     
   ToolHandle<Muon::MuonLayerHoughTool> m_layerTool;     
   bool m_printSummary;
diff --git a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughTool.cxx
index 0a61abd3b3ede14f265c45bc0cbbd26a9dd383e0..3070ff15dc2e0801f0822c67a299b311b2a09279 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughTool.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools/src/MuonLayerHoughTool.cxx
@@ -15,11 +15,11 @@
 #include "TrkTruthData/PRD_MultiTruthCollection.h"
 #include "HepMC/GenEvent.h"
 #include "GaudiKernel/IIncidentSvc.h"
+#include "GaudiKernel/ConcurrencyFlags.h"
 #include "CxxUtils/sincos.h"
 #include "xAODTruth/TruthParticle.h"
 #include "xAODTruth/TruthParticleContainer.h"
 #include "xAODMuon/MuonSegmentContainer.h"
-
 namespace Muon {
 
   MuonLayerHoughTool::MuonLayerHoughTool(const std::string& type, const std::string& name, const IInterface* parent):
@@ -74,12 +74,17 @@ namespace Muon {
     ATH_CHECK( detStore()->retrieve( m_detMgr ) );
 
     if( m_doNtuple ){
-      TDirectory* cdir = gDirectory;
-      m_file = new TFile("HitNtuple.root","RECREATE");
-      m_tree = new TTree("data","data");
-      m_ntuple = new MuonHough::HitNtuple();
-      m_ntuple->initForWrite(*m_tree);
-      gDirectory = cdir;
+      if (Gaudi::Concurrency::ConcurrencyFlags::concurrent()) {
+        // Disabled under concurrency due to thread-safety concerns, but we want to keep it as a debug tool
+        ATH_MSG_DEBUG("HitNtuple disabled because of concurrency");
+      } else {
+        TDirectory* cdir = gDirectory;
+        m_file = new TFile("HitNtuple.root","RECREATE");
+        m_tree = new TTree("data","data");
+        m_ntuple = new MuonHough::HitNtuple();
+        m_ntuple->initForWrite(*m_tree);
+        gDirectory = cdir;
+      }
     }else{
       m_file = 0;
       m_tree = 0;
@@ -165,7 +170,6 @@ namespace Muon {
     //   }
     // }
 
-
     return StatusCode::SUCCESS;
   }
 
@@ -235,39 +239,26 @@ namespace Muon {
   }
   
   void MuonLayerHoughTool::reset() const {
-    m_foundTruthHits.clear();
-    m_outputTruthHits.clear();
-    m_truthHits.clear();
-    // ?!? auto
-    for( HoughDataPerSectorVec::iterator hdit=m_houghDataPerSectorVec.begin();hdit!=m_houghDataPerSectorVec.end();++hdit ) hdit->cleanUp();
-    m_houghDataPerSectorVec.clear();
-    m_seedMaxima.clear();
-    std::vector<TgcHitClusteringObj*>::iterator clit = m_tgcClusteringObjs.begin();
-    std::vector<TgcHitClusteringObj*>::iterator clit_end = m_tgcClusteringObjs.end();
-    for( ;clit!=clit_end;++clit ) delete *clit;
-    m_tgcClusteringObjs.clear();
-
-    m_detectorHoughTransforms.reset();
-
     if( m_ntuple )  m_ntuple->reset();
-
   }
 
 
-  MuonPatternCombinationCollection* MuonLayerHoughTool::find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
-							      const std::vector<const CscPrepDataCollection*>& ,  
-							      const std::vector<const TgcPrepDataCollection*>& tgcCols,  
-							      const std::vector<const RpcPrepDataCollection*>& rpcCols,  
-							      const MuonSegmentCombinationCollection* ) const {                                              
-
+  auto MuonLayerHoughTool::find( 
+      const std::vector<const MdtPrepDataCollection*>& mdtCols,  
+      const std::vector<const CscPrepDataCollection*>& ,  
+      const std::vector<const TgcPrepDataCollection*>& tgcCols,  
+      const std::vector<const RpcPrepDataCollection*>& rpcCols,  
+      const MuonSegmentCombinationCollection* ) const 
+      -> std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> {
     reset();
+    State state;
     ATH_MSG_DEBUG("MuonLayerHoughTool::find");
     if( m_doTruth ) getTruth();
 
 
     // create structure to hold data per sector and set the sector indices
-    m_houghDataPerSectorVec.resize(16);
-    for( unsigned int i=0;i<m_houghDataPerSectorVec.size();++i ) m_houghDataPerSectorVec[i].sector=i+1;
+    (*state.houghDataPerSectorVec).resize(16);
+    for( unsigned int i=0;i<(*state.houghDataPerSectorVec).size();++i ) (*state.houghDataPerSectorVec)[i].sector=i+1;
 
     // return DetectorRegionIndex and sectorLayerHash
     auto getHashes = [this]( const Identifier& id ){
@@ -282,7 +273,7 @@ namespace Muon {
       Identifier id = col->identify();
       int sector = m_idHelper->sector(id);
       auto hashes = getHashes(id);
-      fill(*col,m_houghDataPerSectorVec[sector-1].hitVec[hashes.second]);
+      fill(state,*col,(*state.houghDataPerSectorVec)[sector-1].hitVec[hashes.second]);
     }
 
     for( auto col : rpcCols ){
@@ -290,7 +281,7 @@ namespace Muon {
       Identifier id = col->identify();
       int sector = m_idHelper->sector(id);
       auto hashes = getHashes(id);
-      fill(*col,m_houghDataPerSectorVec[sector-1].hitVec[hashes.second],m_houghDataPerSectorVec[sector-1].phiHitVec[hashes.first]);
+      fill(state,*col,(*state.houghDataPerSectorVec)[sector-1].hitVec[hashes.second],(*state.houghDataPerSectorVec)[sector-1].phiHitVec[hashes.first]);
     }
 
     auto hashInSector = [this]( IdentifierHash hash, int sector, unsigned int sectorLayerHash ) {
@@ -304,36 +295,39 @@ namespace Muon {
       int sector = m_idHelper->sector(id);
       auto hashes = getHashes(id);
       // fill current sector
-      fill(*col,m_houghDataPerSectorVec[sector-1].hitVec[hashes.second],m_houghDataPerSectorVec[sector-1].phiHitVec[hashes.first],sector);
+      fill(state, *col,(*state.houghDataPerSectorVec)[sector-1].hitVec[hashes.second],
+          (*state.houghDataPerSectorVec)[sector-1].phiHitVec[hashes.first],sector);
 
       // fill neighbours if in overlap
       int neighbourSectorDown = sector == 1 ? 16 : sector-1;
       if( hashInSector(col->identifyHash(),neighbourSectorDown,hashes.second) ) 
-        fill(*col,m_houghDataPerSectorVec[neighbourSectorDown-1].hitVec[hashes.second],
-             m_houghDataPerSectorVec[neighbourSectorDown-1].phiHitVec[hashes.first],neighbourSectorDown);
+        fill(state, *col,(*state.houghDataPerSectorVec)[neighbourSectorDown-1].hitVec[hashes.second],
+             (*state.houghDataPerSectorVec)[neighbourSectorDown-1].phiHitVec[hashes.first],neighbourSectorDown);
 
       int neighbourSectorUp   = sector == 16 ? 1 : sector+1;
       if( hashInSector(col->identifyHash(),neighbourSectorUp,hashes.second) ) 
-        fill(*col,m_houghDataPerSectorVec[neighbourSectorUp-1].hitVec[hashes.second],
-             m_houghDataPerSectorVec[neighbourSectorUp-1].phiHitVec[hashes.first],neighbourSectorUp);
+        fill(state, *col,(*state.houghDataPerSectorVec)[neighbourSectorUp-1].hitVec[hashes.second],
+             (*state.houghDataPerSectorVec)[neighbourSectorUp-1].phiHitVec[hashes.first],neighbourSectorUp);
       
     }
     
-    return analyse();
+    return analyse(state);
   }
 
-  MuonPatternCombinationCollection* MuonLayerHoughTool::analyse( const MdtPrepDataContainer*  mdtCont,
-                                                                 const CscPrepDataContainer*  cscCont,
-                                                                 const TgcPrepDataContainer*  tgcCont,
-                                                                 const RpcPrepDataContainer*  rpcCont,
-                                                                 const sTgcPrepDataContainer* stgcCont,  
-                                                                 const MMPrepDataContainer*   mmCont ) const {
-
+  auto MuonLayerHoughTool::analyse(
+      const MdtPrepDataContainer*  mdtCont,
+      const CscPrepDataContainer*  cscCont,
+      const TgcPrepDataContainer*  tgcCont,
+      const RpcPrepDataContainer*  rpcCont,
+      const sTgcPrepDataContainer* stgcCont,  
+      const MMPrepDataContainer*   mmCont ) const 
+      -> std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> {
     reset();
+    State state; 
     ATH_MSG_DEBUG("MuonLayerHoughTool::analyse");
     if( m_doTruth ) getTruth();
 
-    m_houghDataPerSectorVec.resize(16);
+    (*state.houghDataPerSectorVec).resize(16);
 
     // loops over all sectors, contains hashes for technology and region and chamber (?)
     CollectionsPerSectorCit sit = m_collectionsPerSector.begin();
@@ -342,22 +336,23 @@ namespace Muon {
             
       ATH_MSG_DEBUG("analyse: Filling hits sector " << sit->sector);
 
-      HoughDataPerSector& houghData = m_houghDataPerSectorVec[sit->sector-1]; 
+      HoughDataPerSector& houghData = (*state.houghDataPerSectorVec)[sit->sector-1]; 
       houghData.sector = sit->sector;
 
       // fill hits for this sector -> hitsVec and PhiHitsVec are known now
-      fillHitsPerSector( *sit, mdtCont,cscCont,tgcCont,rpcCont,stgcCont,mmCont,houghData);
+      fillHitsPerSector( state, *sit, mdtCont,cscCont,tgcCont,rpcCont,stgcCont,mmCont,houghData);
 
     }
-    return analyse();
+    return analyse(state);
   }
   
-  MuonPatternCombinationCollection* MuonLayerHoughTool::analyse() const {    
-
-    MuonPatternCombinationCollection* patternCombis = new MuonPatternCombinationCollection();
+  auto MuonLayerHoughTool::analyse(State& state) const 
+      -> std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<HoughDataPerSectorVec>> {    
 
+    auto patternCombis = std::make_unique<MuonPatternCombinationCollection>();
+    
     // loop over data and fill the hough transform
-    for( auto& houghData : m_houghDataPerSectorVec ){
+    for( auto& houghData : (*state.houghDataPerSectorVec) ){
 
       ATH_MSG_DEBUG("analyse: Filling Hough sector " << houghData.sector);
 
@@ -375,7 +370,7 @@ namespace Muon {
         MuonStationIndex::StIndex             index = MuonStationIndex::toStationIndex(region,layer);
 
         // get Hough transform
-        MuonHough::MuonLayerHough& hough = m_detectorHoughTransforms.hough( houghData.sector,  region, layer );
+        MuonHough::MuonLayerHough& hough = state.detectorHoughTransforms.hough( houghData.sector,  region, layer );
        
         ATH_MSG_DEBUG("analyse: Filling Summary: loc s"   << houghData.sector 
                       << " "           << MuonStationIndex::regionName(region) 
@@ -384,7 +379,7 @@ namespace Muon {
                       << " etaHits: "  << hits.size() );
 
         // look for maxima using hough in eta per layer
-        if( !findMaxima(hough,hits,houghData.maxVec[layerHash]) || houghData.maxVec[layerHash].empty() ) continue;
+        if( !findMaxima(state, hough,hits,houghData.maxVec[layerHash]) || houghData.maxVec[layerHash].empty() ) continue;
 
         ++houghData.nlayersWithMaxima[region];
         houghData.nmaxHitsInRegion[region] += houghData.maxVec[layerHash].front()->max;
@@ -400,7 +395,7 @@ namespace Muon {
 
     if( m_useSeeds ){
       std::vector<Road> roads;
-      buildRoads(roads);
+      buildRoads(state, roads);
       
       // create association map
       ATH_MSG_DEBUG("analyse: Building pattern combinations using roads " << roads.size() );
@@ -423,15 +418,15 @@ namespace Muon {
             phiEtaAssMap[&max] = road.maxima;
           }
         }
-        createPatternCombinations(phiEtaAssMap,*patternCombis);
+        createPatternCombinations(state,phiEtaAssMap,*patternCombis);
         createPatternCombinations(unassociatedEtaMaxima,*patternCombis);
       }
 
     }
     else{
       // now that the full hough transform is filled, order sectors by maxima
-      std::vector<HoughDataPerSector*> sectorData(m_houghDataPerSectorVec.size());
-      for( unsigned int i=0;i<m_houghDataPerSectorVec.size();++i) sectorData[i] = &m_houghDataPerSectorVec[i];
+      std::vector<HoughDataPerSector*> sectorData((*state.houghDataPerSectorVec).size());
+      for( unsigned int i=0;i<(*state.houghDataPerSectorVec).size();++i) sectorData[i] = &(*state.houghDataPerSectorVec)[i];
       std::stable_sort(sectorData.begin(),sectorData.end(),SortHoughDataPerSector());
 
       std::vector<HoughDataPerSector*>::iterator spit = sectorData.begin();
@@ -453,7 +448,7 @@ namespace Muon {
                         << " layers with phi maxima " << houghData.nphilayersWithMaxima[region] << " hits " << houghData.nphimaxHitsInRegion[region] );
   
           // look for maxima in the overlap regions of sectors
-          associateMaximaInNeighbouringSectors(houghData,m_houghDataPerSectorVec);
+          associateMaximaInNeighbouringSectors(houghData,(*state.houghDataPerSectorVec));
 
           // layers in this region
           int nlayers = MuonStationIndex::LayerIndexMax;
@@ -464,7 +459,7 @@ namespace Muon {
           associateMaximaToPhiMaxima( region, houghData,  phiEtaAssociations, unassociatedEtaMaxima );
 
           // create pattern combinations for combined patterns 
-          createPatternCombinations(phiEtaAssociations, *patternCombis );
+          createPatternCombinations(state, phiEtaAssociations, *patternCombis );
   
           // create pattern combinations for unassociated patterns 
           createPatternCombinations(unassociatedEtaMaxima,*patternCombis);
@@ -473,7 +468,7 @@ namespace Muon {
     }
     
     if( m_ntuple ) {
-      fillNtuple(m_houghDataPerSectorVec);
+      fillNtuple((*state.houghDataPerSectorVec));
       m_tree->Fill();
     }
 
@@ -481,21 +476,21 @@ namespace Muon {
 
     if( m_doTruth && msgLvl(MSG::DEBUG) ){
       ATH_MSG_DEBUG("Hough performance ");
-      printTruthSummary(m_truthHits,m_foundTruthHits);
+      printTruthSummary(state.truthHits,state.foundTruthHits);
       ATH_MSG_DEBUG("Association performance ");
-      printTruthSummary(m_foundTruthHits,m_outputTruthHits);
+      printTruthSummary(state.foundTruthHits,state.outputTruthHits);
     }
 
-    return patternCombis;
+    return {std::move(patternCombis), std::move(state.houghDataPerSectorVec)};
   }
 
-  void MuonLayerHoughTool::buildRoads( std::vector<MuonLayerHoughTool::Road>& roads ) const {
+  void MuonLayerHoughTool::buildRoads(State& state, std::vector<MuonLayerHoughTool::Road>& roads ) const {
     // sort maxima according to hits
-    std::stable_sort( m_seedMaxima.begin(),m_seedMaxima.end(), [](const MuonHough::MuonLayerHough::Maximum* m1,
+    std::stable_sort( state.seedMaxima.begin(),state.seedMaxima.end(), [](const MuonHough::MuonLayerHough::Maximum* m1,
                                                                   const MuonHough::MuonLayerHough::Maximum* m2 ){ return m1->max > m2->max; } );
     // loop over seed maxima (which are maxima) that pass certain thresholds detailed in cut_values
     std::set<const MuonHough::MuonLayerHough::Maximum*> associatedMaxima;
-    for( auto maxit = m_seedMaxima.begin();maxit!=m_seedMaxima.end(); ++maxit ){
+    for( auto maxit = state.seedMaxima.begin();maxit!=state.seedMaxima.end(); ++maxit ){
       // if this maximum is already in the set of associated maxima, do not do anything
       if( associatedMaxima.count(*maxit) ) continue;
 
@@ -518,7 +513,7 @@ namespace Muon {
       bool isNSW=m_idHelper->issTgc(seed.hits[0]->prd->identify()) || m_idHelper->isMM(seed.hits[0]->prd->identify());
       // extend seed within the current sector
       // sector indices have an offset of -1 because the numbering of the sectors are from 1 to 16 but the indices in the vertices are of course 0 to 15
-      extendSeed( road, m_houghDataPerSectorVec[sector-1] );
+      extendSeed( state, road, (*state.houghDataPerSectorVec)[sector-1] );
 
       // look for maxima in the overlap regions of sectors
       int sectorN            = sector-1;
@@ -527,28 +522,28 @@ namespace Muon {
       if(sectorP>16) sectorP = 1;
 
       // associate the road with phi maxima
-      associatePhiMaxima( road, m_houghDataPerSectorVec[sector-1].phiMaxVec[region] );
+      associatePhiMaxima( road, (*state.houghDataPerSectorVec)[sector-1].phiMaxVec[region] );
       //
       if(m_addSectors && isNSW) {
-        extendSeed( road, m_houghDataPerSectorVec[sectorN-1] );
-        associatePhiMaxima( road, m_houghDataPerSectorVec[sectorN-1].phiMaxVec[region] );
-        extendSeed( road, m_houghDataPerSectorVec[sectorP-1] );
-        associatePhiMaxima( road, m_houghDataPerSectorVec[sectorP-1].phiMaxVec[region] );
+        extendSeed(state, road, (*state.houghDataPerSectorVec)[sectorN-1] );
+        associatePhiMaxima( road, (*state.houghDataPerSectorVec)[sectorN-1].phiMaxVec[region] );
+        extendSeed(state, road, (*state.houghDataPerSectorVec)[sectorP-1] );
+        associatePhiMaxima( road, (*state.houghDataPerSectorVec)[sectorP-1].phiMaxVec[region] );
       }
 
       if( road.neighbouringRegion != MuonStationIndex::DetectorRegionUnknown ) {
-        associatePhiMaxima( road, m_houghDataPerSectorVec[sector-1].phiMaxVec[road.neighbouringRegion] );
+        associatePhiMaxima( road, (*state.houghDataPerSectorVec)[sector-1].phiMaxVec[road.neighbouringRegion] );
       }
       // if close to a sector boundary, try adding maxima in that sector as well
       if( road.neighbouringSector != -1 ) {
         ATH_MSG_DEBUG("  Adding neighbouring sector " << road.neighbouringSector );
-        extendSeed( road, m_houghDataPerSectorVec[road.neighbouringSector-1] );
-        associatePhiMaxima( road, m_houghDataPerSectorVec[road.neighbouringSector-1].phiMaxVec[region] );
+        extendSeed(state, road, (*state.houghDataPerSectorVec)[road.neighbouringSector-1] );
+        associatePhiMaxima( road, (*state.houghDataPerSectorVec)[road.neighbouringSector-1].phiMaxVec[region] );
       }
 
       // finally deal with the case that we have both neighbouring region and sector
       if( road.neighbouringRegion != MuonStationIndex::DetectorRegionUnknown && road.neighbouringSector != -1) {
-        associatePhiMaxima( road, m_houghDataPerSectorVec[ road.neighbouringSector-1].phiMaxVec[road.neighbouringRegion] );
+        associatePhiMaxima( road, (*state.houghDataPerSectorVec)[ road.neighbouringSector-1].phiMaxVec[road.neighbouringRegion] );
       }
 
       // merge phi maxima
@@ -689,7 +684,7 @@ namespace Muon {
   // roads are combinations of maxima
   
   
-  void MuonLayerHoughTool::extendSeed( MuonLayerHoughTool::Road& road, MuonLayerHoughTool::HoughDataPerSector& sectorData ) const { //const {
+  void MuonLayerHoughTool::extendSeed(State& state, MuonLayerHoughTool::Road& road, MuonLayerHoughTool::HoughDataPerSector& sectorData ) const { //const {
     if( !road.seed ) return;
     
     RegionMaximumVec& maxVec = sectorData.maxVec;
@@ -808,7 +803,7 @@ namespace Muon {
       }
     }
       
-    MuonHough::MuonPhiLayerHough& phiHough = m_detectorHoughTransforms.phiHough( region ); // get phi transform in the same region as the seed
+    MuonHough::MuonPhiLayerHough& phiHough = state.detectorHoughTransforms.phiHough( region ); // get phi transform in the same region as the seed
             
     // gather phiHits in sector that match the etaHits of the maximum
     PhiHitVec phiHitsInMaximum;
@@ -829,7 +824,7 @@ namespace Muon {
                                << " phiHitsInMaxima " << phiHitsInMaximum.size()
                                << " phi hits:  "      << phiHits.size() );
 
-    if( !findMaxima(phiHough,phiHitsInMaximum,sectorData.phiMaxVec[region], sectorData.sector) || sectorData.phiMaxVec[region].empty() ) {
+    if( !findMaxima(state, phiHough,phiHitsInMaximum,sectorData.phiMaxVec[region], sectorData.sector) || sectorData.phiMaxVec[region].empty() ) {
       ATH_MSG_DEBUG("extendSeed: No phi maxima found in  s" << sectorData.sector << " " << MuonStationIndex::regionName(region) );
       return;
     }
@@ -1417,7 +1412,7 @@ namespace Muon {
     patternCombis.push_back(combi);
   }
 
-  void MuonLayerHoughTool::createPatternCombinations( std::map< MuonHough::MuonPhiLayerHough::Maximum*, MuonLayerHoughTool::MaximumVec >& phiEtaAssociations,
+  void MuonLayerHoughTool::createPatternCombinations(State& state, std::map< MuonHough::MuonPhiLayerHough::Maximum*, MuonLayerHoughTool::MaximumVec >& phiEtaAssociations,
                   MuonPatternCombinationCollection& patternCombis ) const {
     
     ATH_MSG_DEBUG("Creating pattern combinations from eta/phi combinations " << phiEtaAssociations.size() );
@@ -1585,7 +1580,7 @@ namespace Muon {
   
   if( m_doTruth ){
     for( std::vector<const Trk::PrepRawData*>::iterator it=prds.begin();it!=prds.end();++it ){
-      if( m_truthHits.count((*it)->identify()) ) m_outputTruthHits.insert((*it)->identify());
+      if( state.truthHits.count((*it)->identify()) ) state.outputTruthHits.insert((*it)->identify());
       if( !m_truthSummaryTool.empty() ) m_truthSummaryTool->add((*it)->identify(),2);
     }
   }
@@ -1602,7 +1597,8 @@ namespace Muon {
     }
   }
 
-  bool MuonLayerHoughTool::findMaxima( MuonHough::MuonLayerHough& hough,
+  bool MuonLayerHoughTool::findMaxima( State& state,
+                                       MuonHough::MuonLayerHough& hough,
                                        MuonLayerHoughTool::HitVec& hits, 
                                        MuonLayerHoughTool::MaximumVec& maxima ) const {
     if( hits.empty() ) return false;
@@ -1657,7 +1653,7 @@ namespace Muon {
 
           if( m_doTruth ){
             if( !m_truthSummaryTool.empty() ) m_truthSummaryTool->add(id,1);
-            if( m_truthHits.count(id) )       m_foundTruthHits.insert(id);
+            if( state.truthHits.count(id) )       state.foundTruthHits.insert(id);
           }
 
           ATH_MSG_VERBOSE("findMaxima: hit " << hit.layer << "  " << m_idHelper->toString(id) << " hits " << nhits );
@@ -1667,7 +1663,7 @@ namespace Muon {
         if( nmdt > 0 || (nmm + nstgc) > 0) {
           maxima.push_back( new MuonHough::MuonLayerHough::Maximum(maximum) );
           // add to seed list if 
-          if( maximum.max > m_selectors[hough.m_descriptor.chIndex].getCutValue(maximum.pos) ) m_seedMaxima.push_back(maxima.back());        
+          if( maximum.max > m_selectors[hough.m_descriptor.chIndex].getCutValue(maximum.pos) ) state.seedMaxima.push_back(maxima.back());        
           ++nmaxima;
         }
         hough.fillLayer2(maximum.hits,true);
@@ -1683,7 +1679,8 @@ namespace Muon {
     return true;
   }
 
-  bool MuonLayerHoughTool::findMaxima( MuonHough::MuonPhiLayerHough& hough,
+  bool MuonLayerHoughTool::findMaxima( State& state,
+                                       MuonHough::MuonPhiLayerHough& hough,
                                        MuonLayerHoughTool::PhiHitVec& hits, 
                                        MuonLayerHoughTool::PhiMaximumVec& maxima,
                                        int sector ) const{
@@ -1717,7 +1714,7 @@ namespace Muon {
         
           if( m_doTruth ){
             if( !m_truthSummaryTool.empty() ) m_truthSummaryTool->add(id,1);
-            if( m_truthHits.count(id) )       m_foundTruthHits.insert(id);
+            if( state.truthHits.count(id) )       state.foundTruthHits.insert(id);
           }
           
           int nhits = hit.tgc ? hit.tgc->phiCluster.hitList.size() : 1;
@@ -1771,7 +1768,8 @@ namespace Muon {
     return true;
   }
 
-  void MuonLayerHoughTool::fillHitsPerSector(  const MuonLayerHoughTool::CollectionsPerSector& collectionsPerSector,
+  void MuonLayerHoughTool::fillHitsPerSector(  State& state,
+                                               const MuonLayerHoughTool::CollectionsPerSector& collectionsPerSector,
                                                const MdtPrepDataContainer*  mdtCont,  
                                                const CscPrepDataContainer*  /*cscCont*/,  
                                                const TgcPrepDataContainer*  tgcCont,  
@@ -1793,23 +1791,23 @@ namespace Muon {
           // !?! else if made by Felix
           if( mdtCont && mdtCont->size()>0 && tech == MuonStationIndex::MDT ) {
             MdtPrepDataContainer::const_iterator pos = mdtCont->indexFind(*iit);
-            if( pos != mdtCont->end() ) fill(**pos,houghData.hitVec[layerHash]);
+            if( pos != mdtCont->end() ) fill(state,**pos,houghData.hitVec[layerHash]);
           }
           else if( rpcCont && rpcCont->size()>0 && tech == MuonStationIndex::RPC ) {
             RpcPrepDataContainer::const_iterator pos = rpcCont->indexFind(*iit);
-            if( pos != rpcCont->end() ) fill(**pos,houghData.hitVec[layerHash],houghData.phiHitVec[regionLayer.first]);
+            if( pos != rpcCont->end() ) fill(state,**pos,houghData.hitVec[layerHash],houghData.phiHitVec[regionLayer.first]);
           }
           else if( tgcCont && tgcCont->size()>0 && tech == MuonStationIndex::TGC ) {
             TgcPrepDataContainer::const_iterator pos = tgcCont->indexFind(*iit);
-            if( pos != tgcCont->end() ) fill(**pos,houghData.hitVec[layerHash],houghData.phiHitVec[regionLayer.first],collectionsPerSector.sector);
+            if( pos != tgcCont->end() ) fill(state,**pos,houghData.hitVec[layerHash],houghData.phiHitVec[regionLayer.first],collectionsPerSector.sector);
           }
           else if( stgcCont && stgcCont->size()>0 && tech == MuonStationIndex::STGC ) {
             sTgcPrepDataContainer::const_iterator pos = stgcCont->indexFind(*iit);
-            if( pos != stgcCont->end() ) fill(**pos,houghData.hitVec[layerHash],houghData.phiHitVec[regionLayer.first],collectionsPerSector.sector);
+            if( pos != stgcCont->end() ) fill(state,**pos,houghData.hitVec[layerHash],houghData.phiHitVec[regionLayer.first],collectionsPerSector.sector);
           }
           else if( mmCont && mmCont->size()>0 && tech == MuonStationIndex::MM ) {
             MMPrepDataContainer::const_iterator pos = mmCont->indexFind(*iit);
-            if( pos != mmCont->end() ) fill(**pos,houghData.hitVec[layerHash]);
+            if( pos != mmCont->end() ) fill(state,**pos,houghData.hitVec[layerHash]);
           }
         }
       }
@@ -1847,7 +1845,7 @@ namespace Muon {
     }
   }
 
-  void MuonLayerHoughTool::matchTruth( const PRD_MultiTruthCollection& truthCol, const Identifier& id, MuonHough::HitDebugInfo& debug ) const {
+  void MuonLayerHoughTool::matchTruth(State& state, const PRD_MultiTruthCollection& truthCol, const Identifier& id, MuonHough::HitDebugInfo& debug ) const {
     typedef PRD_MultiTruthCollection::const_iterator iprdt;
     std::pair<iprdt, iprdt> range = truthCol.equal_range(id);
     
@@ -1860,13 +1858,13 @@ namespace Muon {
         if( link.cptr() && abs(link.cptr()->pdg_id()) == 13 ){
           debug.barcode = link.barcode();
           debug.pdgId = link.cptr()->pdg_id();
-          m_truthHits.insert(id);
+          state.truthHits.insert(id);
         }
       }
     }
   }
 
-  void MuonLayerHoughTool::fill( const MdtPrepDataCollection& mdts, MuonLayerHoughTool::HitVec& hits ) const {
+  void MuonLayerHoughTool::fill( State& state, const MdtPrepDataCollection& mdts, MuonLayerHoughTool::HitVec& hits ) const {
     
     if( mdts.empty() ) return;
     auto truthCollections = m_truthNames.makeHandles();
@@ -1900,7 +1898,7 @@ namespace Muon {
       debug->time = prd.tdc();
       debug->r = prd.localPosition()[Trk::locR];
       
-      if( technology < truthCollections.size() ) matchTruth(*truthCollections[technology],id,*debug);
+      if( technology < truthCollections.size() ) matchTruth(state, *truthCollections[technology],id,*debug);
       MuonHough::Hit* hit = new MuonHough::Hit(sublayer,x,ymin,ymax,1.,debug,&prd);
       hits.push_back(hit);
     }
@@ -1913,7 +1911,7 @@ namespace Muon {
 
   }
 
-  void MuonLayerHoughTool::fill( const RpcPrepDataCollection& rpcs, MuonLayerHoughTool::HitVec& hits, MuonLayerHoughTool::PhiHitVec& phiHits ) const {
+  void MuonLayerHoughTool::fill( State& state, const RpcPrepDataCollection& rpcs, MuonLayerHoughTool::HitVec& hits, MuonLayerHoughTool::PhiHitVec& phiHits ) const {
     
     if( rpcs.empty() ) return;
     auto truthCollections = m_truthNames.makeHandles();
@@ -1947,7 +1945,7 @@ namespace Muon {
       debug->isEtaPhi = (neta && nphi);
       debug->trigConfirm = 1;
       debug->time = prd.time();
-      if( technology < truthCollections.size() ) matchTruth(*truthCollections[technology],id,*debug);
+      if( technology < truthCollections.size() ) matchTruth(state, *truthCollections[technology],id,*debug);
 
       float weight = (neta && nphi) ? 2 : 1;
       if( m_idHelper->rpcIdHelper().measuresPhi(id) ) {
@@ -1971,7 +1969,7 @@ namespace Muon {
     }
   }
 
-  void MuonLayerHoughTool::fill( const MMPrepDataCollection& mms, MuonLayerHoughTool::HitVec& hits ) const {
+  void MuonLayerHoughTool::fill( State& state, const MMPrepDataCollection& mms, MuonLayerHoughTool::HitVec& hits ) const {
     
     if( mms.empty() ) return;
     auto truthCollections = m_truthNames.makeHandles();
@@ -1997,14 +1995,14 @@ namespace Muon {
       float ymax = y + stripCor;
       MuonHough::HitDebugInfo* debug = new MuonHough::HitDebugInfo(technology,sector,region,layer,sublayer);
       debug->r = stripCor;
-      if( technology < truthCollections.size() ) matchTruth(*truthCollections[technology],id,*debug);
+      if( technology < truthCollections.size() ) matchTruth(state, *truthCollections[technology],id,*debug);
 
       MuonHough::Hit* hit = new MuonHough::Hit(sublayer,x,ymin,ymax,1.,debug,&prd);
       hits.push_back(hit);
     }
   }
 
-  void MuonLayerHoughTool::fill( const sTgcPrepDataCollection& stgcs, MuonLayerHoughTool::HitVec& hits, MuonLayerHoughTool::PhiHitVec& phiHits, int selectedSector ) const {
+  void MuonLayerHoughTool::fill( State& state, const sTgcPrepDataCollection& stgcs, MuonLayerHoughTool::HitVec& hits, MuonLayerHoughTool::PhiHitVec& phiHits, int selectedSector ) const {
     
     if( stgcs.empty() ) return;
     auto truthCollections = m_truthNames.makeHandles();
@@ -2034,7 +2032,7 @@ namespace Muon {
       debug->isEtaPhi = 1;
       debug->trigConfirm = (prd.getBcBitMap() & sTgcPrepData::BCBIT_CURRENT) == sTgcPrepData::BCBIT_CURRENT;
       debug->time = prd.getBcBitMap();
-      if( technology < truthCollections.size() ) matchTruth(*truthCollections[technology],id,*debug);
+      if( technology < truthCollections.size() ) matchTruth(state, *truthCollections[technology],id,*debug);
 
       if( m_idHelper->stgcIdHelper().channelType(id) == 1 ) {
         // eta strips
@@ -2122,12 +2120,12 @@ namespace Muon {
     }
   }
 
-  void MuonLayerHoughTool::fill( const TgcPrepDataCollection& tgcs, MuonLayerHoughTool::HitVec& hits,
+  void MuonLayerHoughTool::fill( State& state, const TgcPrepDataCollection& tgcs, MuonLayerHoughTool::HitVec& hits,
          MuonLayerHoughTool::PhiHitVec& phiHits, int sector ) const {
     
     if( tgcs.empty() ) return;
-    m_tgcClusteringObjs.push_back( new TgcHitClusteringObj(m_idHelper->tgcIdHelper()) );
-    TgcHitClusteringObj& clustering = *m_tgcClusteringObjs.back();
+    state.tgcClusteringObjs.push_back( new TgcHitClusteringObj(m_idHelper->tgcIdHelper()) );
+    TgcHitClusteringObj& clustering = *state.tgcClusteringObjs.back();
     std::vector<const TgcPrepData*> prds;
     prds.insert(prds.begin(),tgcs.begin(),tgcs.end());
     clustering.cluster(prds);
@@ -2197,7 +2195,7 @@ namespace Muon {
   debug->clusterLayers = cl.etaCluster.layers();
   debug->isEtaPhi = cl.phiCluster.layers();
   debug->time = cl.etaCluster.hitList.front()->getBcBitMap();
-        if( technology < truthCollections.size() ) matchTruth(*truthCollections[technology],id,*debug);
+        if( technology < truthCollections.size() ) matchTruth(state, *truthCollections[technology],id,*debug);
 
   MuonHough::HitDebugInfo* phiDebug = new MuonHough::HitDebugInfo(*debug);
   phiDebug->clusterSize = cl.phiCluster.hitList.size();
@@ -2365,24 +2363,6 @@ namespace Muon {
     if( msgLvl(MSG::DEBUG) ) msg(MSG::DEBUG) << endmsg;
   }
 
-  void MuonLayerHoughTool::HoughDataPerSector::cleanUp() {
-    for(RegionHitVec::iterator it=hitVec.begin();it!=hitVec.end();++it)
-      for( HitVec::iterator it2=it->begin();it2!=it->end();++it2 ) delete *it2;
-    hitVec.clear();
-
-    for(RegionPhiHitVec::iterator it=phiHitVec.begin();it!=phiHitVec.end();++it)
-      for( PhiHitVec::iterator it2=it->begin();it2!=it->end();++it2 ) delete *it2;
-    phiHitVec.clear();
-
-    for(RegionMaximumVec::iterator it=maxVec.begin();it!=maxVec.end();++it)
-      for( MaximumVec::iterator it2=it->begin();it2!=it->end();++it2 ) delete *it2;
-    maxVec.clear();
-
-    for(RegionPhiMaximumVec::iterator it=phiMaxVec.begin();it!=phiMaxVec.end();++it)
-      for( PhiMaximumVec::iterator it2=it->begin();it2!=it->end();++it2 ) delete *it2;
-    phiMaxVec.clear();
-  }
-
   void MuonLayerHoughTool::printTruthSummary( std::set<Identifier>& truth, std::set<Identifier>& found ) const {
     if( truth.size() == found.size() ){
       ATH_MSG_DEBUG(" All hits found: truth " << truth.size() << " found " << found.size() );
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/CMakeLists.txt b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/CMakeLists.txt
index b7b63bf22183e6748c4e3bcf3a8d5db435ce910a..adeda837853846d212aca0cc344512de1229a079 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/CMakeLists.txt
+++ b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/CMakeLists.txt
@@ -14,6 +14,7 @@ atlas_depends_on_subdirs( PUBLIC
 			  Generators/GeneratorObjects
                           GaudiKernel
                           MuonSpectrometer/MuonIdHelpers
+                          MuonSpectrometer/MuonReconstruction/MuonPatternFinders/MuonPatternFinderTools/MuonHoughPatternTools
                           MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonLayerEvent
                           MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonPattern
                           MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonPrepRawData
diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonHoughPatternFinderTool.h b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonHoughPatternFinderTool.h
index 19f7b023ac504c9ca1bd46842683bf6cd10385e8..5eb0e9821c65629b0437e1b28566a900894e5d38 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonHoughPatternFinderTool.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonRecTools/MuonRecToolInterfaces/MuonRecToolInterfaces/IMuonHoughPatternFinderTool.h
@@ -15,6 +15,8 @@
 #include "MuonSegment/MuonSegmentCombinationCollection.h"
 #include "MuonPattern/MuonPatternCombinationCollection.h"
 
+#include "MuonHoughPatternTools/HoughDataPerSec.h"
+
 #include <vector>
 
 static const InterfaceID IID_IMuonHoughPatternFinderTool("Muon::IMuonHoughPatternFinderTool",1,0);
@@ -33,11 +35,12 @@ namespace Muon {
     static const InterfaceID& interfaceID();
 
     /** find patterns for a give set of MuonPrepData collections + optionally CSC segment combinations */
-    virtual MuonPatternCombinationCollection* find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
-						    const std::vector<const CscPrepDataCollection*>& cscCols,  
-						    const std::vector<const TgcPrepDataCollection*>& tgcCols,  
-						    const std::vector<const RpcPrepDataCollection*>& rpcCols,  
-						    const MuonSegmentCombinationCollection* cscSegmentCombis ) const = 0;
+    virtual std::pair<std::unique_ptr<MuonPatternCombinationCollection>, std::unique_ptr<std::vector<Muon::HoughDataPerSec>>>
+    find( const std::vector<const MdtPrepDataCollection*>& mdtCols,  
+          const std::vector<const CscPrepDataCollection*>& cscCols,  
+          const std::vector<const TgcPrepDataCollection*>& tgcCols,  
+          const std::vector<const RpcPrepDataCollection*>& rpcCols,  
+          const MuonSegmentCombinationCollection* cscSegmentCombis ) const = 0;
 
   };
   
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentCombinationFinder.cxx b/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentCombinationFinder.cxx
index 01ef6472e8142ab1ce02ca15d2e83a24341975b3..2f69d38e365cc0c6b9fc234e49ce99ae99ca9aa0 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentCombinationFinder.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentCombinationFinder.cxx
@@ -224,7 +224,9 @@ void Muon::MooSegmentCombinationFinder::findSegments( const std::vector<const Md
     if( m_doMdtSegments ){
       // search for global patterns 
       auditorBefore( m_houghPatternFinder );
-      output.patternCombinations = m_houghPatternFinder->find( mdtCols, cscCols, tgcCols, rpcCols, csc4dSegmentCombinations.get() );
+      auto [combis, houghData] = m_houghPatternFinder->find( mdtCols, cscCols, tgcCols, rpcCols, csc4dSegmentCombinations.get() );
+      output.patternCombinations = combis.release();
+      output.houghDataPerSectorVec = std::move(houghData);
       auditorAfter( m_houghPatternFinder, output.patternCombinations );
       printSummary( "Pattern finding", output.patternCombinations );
 
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentFinderAlg.cxx b/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentFinderAlg.cxx
index 0c8d0f5636cd9e1f2cf5ebdf2edaf63125a226be..c1e0b81a2c5f17cd4d95b73b8df3dbf4ca2993ed 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentFinderAlg.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentFinderAlg.cxx
@@ -77,6 +77,7 @@ StatusCode MooSegmentFinderAlg::initialize()
 
   ATH_CHECK( m_patternCombiLocation.initialize() );
   ATH_CHECK( m_segmentLocation.initialize() );
+  ATH_CHECK( m_houghDataPerSectorVecKey.initialize() );
   
   return StatusCode::SUCCESS; 
 }
@@ -126,6 +127,14 @@ StatusCode MooSegmentFinderAlg::execute()
     }
   }
 
+  // write hough data to SG
+  if (output.houghDataPerSectorVec) {
+    SG::WriteHandle<std::vector<Muon::HoughDataPerSec>> handle {m_houghDataPerSectorVecKey};
+    ATH_CHECK(handle.record(std::move(output.houghDataPerSectorVec)));
+  } else {
+    ATH_MSG_VERBOSE("HoughDataPerSectorVec was empty, key: " << m_houghDataPerSectorVecKey.key());
+  }
+
   //do cluster based segment finding
   if (m_doTGCClust || m_doRPCClust){
     const PRD_MultiTruthCollection* tgcTruthColl=0;
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentFinderAlg.h b/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentFinderAlg.h
index 781784d425d570660786ca6e59609b547af34437..4aaf2a3a4f10657941892493fb4e35e2d1ef3e90 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentFinderAlg.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MooSegmentCombinationFinder/src/MooSegmentFinderAlg.h
@@ -19,6 +19,8 @@
 #include "TrkSegment/SegmentCollection.h"
 #include "MuonPrepRawData/MuonPrepDataContainer.h"
 #include "MuonPattern/MuonPatternCombinationCollection.h"
+#include "MuonHoughPatternTools/HoughDataPerSec.h"
+
 class MsgStream;
 class StoreGateSvc;
 
@@ -72,6 +74,8 @@ class MooSegmentFinderAlg : public AthAlgorithm
   
   SG::WriteHandleKey<MuonPatternCombinationCollection>   m_patternCombiLocation;
   SG::WriteHandleKey<Trk::SegmentCollection>                   m_segmentLocation;
+  SG::WriteHandleKey<std::vector<Muon::HoughDataPerSec>> m_houghDataPerSectorVecKey {this, 
+    "Key_MuonLayerHoughToolHoughDataPerSectorVec", "HoughDataPerSectorVec", "HoughDataPerSectorVec key"};
 
   ToolHandle<Muon::IMooSegmentCombinationFinder> m_segmentFinder;     //<! pointer to the segment finder
   ToolHandle<Muon::IMuonClusterSegmentFinder> m_clusterSegMaker;
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MuonSegmentCombinerToolInterfaces/MuonSegmentCombinerToolInterfaces/IMooSegmentCombinationFinder.h b/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MuonSegmentCombinerToolInterfaces/MuonSegmentCombinerToolInterfaces/IMooSegmentCombinationFinder.h
index e25899196ac3cf8058e1ecbc0c26ea45828f7369..e2e59b91e51aca5697e09337c4d6a5cedda9b6c8 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MuonSegmentCombinerToolInterfaces/MuonSegmentCombinerToolInterfaces/IMooSegmentCombinationFinder.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentCombiners/MuonSegmentCombinerTools/MuonSegmentCombinerToolInterfaces/MuonSegmentCombinerToolInterfaces/IMooSegmentCombinationFinder.h
@@ -17,6 +17,7 @@
 #include "MuonPrepRawData/CscPrepDataCollection.h"
 #include "MuonPrepRawData/RpcPrepDataCollection.h"
 #include "MuonPrepRawData/TgcPrepDataCollection.h"
+#include "MuonHoughPatternTools/HoughDataPerSec.h"
 
 namespace Muon 
 {
@@ -30,6 +31,7 @@ namespace Muon
     struct Output {
       MuonPatternCombinationCollection* patternCombinations;
       Trk::SegmentCollection*           segmentCollection;
+      std::unique_ptr<std::vector<HoughDataPerSec>> houghDataPerSectorVec;
 
     Output() : patternCombinations(0),segmentCollection(0) {}
     };
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MuonLayerSegmentMakerTools/src/MuonLayerSegmentFinderTool.cxx b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MuonLayerSegmentMakerTools/src/MuonLayerSegmentFinderTool.cxx
index f4b29b8d460682ec1b3697c907a341b5773a4320..ea588d6cf0e082c494fb2bfe0ca38b32a988de9b 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MuonLayerSegmentMakerTools/src/MuonLayerSegmentFinderTool.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MuonLayerSegmentMakerTools/src/MuonLayerSegmentFinderTool.cxx
@@ -67,7 +67,7 @@ namespace Muon {
     ATH_CHECK(m_clusterSegMakerNSW.retrieve());
     ATH_CHECK(m_layerHoughTool.retrieve());
     if( !m_recoValidationTool.empty() ) ATH_CHECK(m_recoValidationTool.retrieve());
-
+    ATH_CHECK(m_houghDataPerSectorVecKey.initialize());
     return StatusCode::SUCCESS;
   }
 
@@ -101,20 +101,26 @@ namespace Muon {
     MuonStationIndex::DetectorRegionIndex regionIndex = intersection.layerSurface.regionIndex;
     MuonStationIndex::LayerIndex  layerIndex  = intersection.layerSurface.layerIndex;
 
+    // get hough data
+    SG::ReadHandle<MuonLayerHoughTool::HoughDataPerSectorVec> houghDataPerSectorVec {m_houghDataPerSectorVecKey};
+    if (!houghDataPerSectorVec.isValid()) {
+      ATH_MSG_ERROR("Hough data per sector vector not found");
+      return;
+    }
 
     // sanity check
-    if( static_cast<int>(m_layerHoughTool->houghData().size()) <= sector-1 ){
-      ATH_MSG_WARNING(" m_layerHoughTool->houghData().size() smaller than sector " << m_layerHoughTool->houghData().size()
+    if( static_cast<int>(houghDataPerSectorVec->size()) <= sector-1 ){
+      ATH_MSG_WARNING(" m_layerHoughTool->houghData().size() smaller than sector " << houghDataPerSectorVec->size()
                       << " sector " << sector );
       return;
     }
 
     // get hough maxima in the layer
     unsigned int sectorLayerHash = MuonStationIndex::sectorLayerHash( regionIndex,layerIndex );
-    const MuonLayerHoughTool::HoughDataPerSector& houghDataPerSector = m_layerHoughTool->houghData()[sector-1];
+    const MuonLayerHoughTool::HoughDataPerSector& houghDataPerSector = (*houghDataPerSectorVec)[sector-1];
     ATH_MSG_DEBUG(" findMdtSegmentsFromHough: sector " << sector << " " << MuonStationIndex::regionName(regionIndex)
                   << " " << MuonStationIndex::layerName(layerIndex) << " sector hash " << sectorLayerHash
-                  << " houghData " << m_layerHoughTool->houghData().size() << " " << houghDataPerSector.maxVec.size());
+                  << " houghData " << houghDataPerSectorVec->size() << " " << houghDataPerSector.maxVec.size());
 
     // sanity check
     if( houghDataPerSector.maxVec.size() <= sectorLayerHash ){
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MuonLayerSegmentMakerTools/src/MuonLayerSegmentFinderTool.h b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MuonLayerSegmentMakerTools/src/MuonLayerSegmentFinderTool.h
index 43a860fd113a87176891d46a5344ddb2c464a0dd..8ca36cf34206957079eeabae875e2f4b5f7cb968 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MuonLayerSegmentMakerTools/src/MuonLayerSegmentFinderTool.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MuonLayerSegmentMakerTools/src/MuonLayerSegmentFinderTool.h
@@ -16,6 +16,8 @@
 #include "MuonLayerEvent/MuonSystemExtension.h"
 #include "MuonDetDescrUtils/MuonSectorMapping.h"
 
+#include "MuonHoughPatternTools/MuonLayerHoughTool.h"
+
 class ICscSegmentFinder;
 namespace Muon {
   
@@ -29,7 +31,6 @@ namespace Muon {
   class IMuonClusterSegmentFinder;
   class IMuonClusterSegmentFinderTool;
   class IMuonRecoValidationTool;
-  class MuonLayerHoughTool;
   class MdtDriftCircleOnTrack;
   class MuonClusterOnTrack;
   class MuonLayerSegmentFinderTool : virtual public IMuonLayerSegmentFinderTool,  public AthAlgTool {
@@ -66,6 +67,10 @@ namespace Muon {
                           const std::vector<const MuonClusterOnTrack*>&    clusters,
                           std::vector< std::shared_ptr<const Muon::MuonSegment> >& segments ) const;
 
+    /** storegate */
+    SG::ReadHandleKey<MuonLayerHoughTool::HoughDataPerSectorVec> m_houghDataPerSectorVecKey {this, 
+        "Key_MuonLayerHoughToolHoughDataPerSectorVec", "HoughDataPerSectorVec", "HoughDataPerSectorVec key"};
+
     /** tool handles */
     ToolHandle<MuonIdHelperTool>                      m_idHelper; 
     ToolHandle<MuonEDMPrinterTool>                    m_printer; 
diff --git a/Reconstruction/MuonIdentification/MuGirl/src/MuGirlRecoTool.cxx b/Reconstruction/MuonIdentification/MuGirl/src/MuGirlRecoTool.cxx
index ddfe49110a00c9f09a1860c058187758c21201f1..b7ab1b422d55ffd6e22003b9c125f3dbe7d011d2 100644
--- a/Reconstruction/MuonIdentification/MuGirl/src/MuGirlRecoTool.cxx
+++ b/Reconstruction/MuonIdentification/MuGirl/src/MuGirlRecoTool.cxx
@@ -230,6 +230,8 @@ StatusCode MuGirlRecoTool::initialize() {
 
   if (!m_caloExtensionTool.empty()) ATH_CHECK(m_caloExtensionTool.retrieve());
 
+  ATH_CHECK(m_houghDataPerSectorVecKey.initialize());
+
   ATH_MSG_INFO("MuGirlrecoTool initialized successfully");
 
   return StatusCode::SUCCESS;
@@ -704,12 +706,16 @@ StatusCode MuGirlRecoTool::processHoughData() {
   int NumMaxima = 0;
 
   if (m_pMuonLayerHoughTool) {
-    Muon::MuonLayerHoughTool::HoughDataPerSectorVec data = m_pMuonLayerHoughTool->houghData();
-    ATH_MSG_DEBUG(data.size() << " sector are present in the HoughData");
-    for (unsigned int sector = 0; sector < data.size(); sector++) {
+    SG::ReadHandle<Muon::MuonLayerHoughTool::HoughDataPerSectorVec> data {m_houghDataPerSectorVecKey};
+    if (!data.isValid()) {
+      ATH_MSG_ERROR("Hough data per sector vector not found");
+      return StatusCode::FAILURE;
+    }
+    ATH_MSG_DEBUG(data->size() << " sector are present in the HoughData");
+    for (unsigned int sector = 0; sector < data->size(); sector++) {
       int sector_id = -99;
       bool hits_in_same_sector = true;
-      Muon::MuonLayerHoughTool::HoughDataPerSector sector_data = data.at(sector);
+      Muon::MuonLayerHoughTool::HoughDataPerSector sector_data = data->at(sector);
       ATH_MSG_DEBUG("----------------------------- Sector " << sector
                                                             << " -----------------------------");
       for (unsigned int region = 0; region < sector_data.maxVec.size(); region++) {
diff --git a/Reconstruction/MuonIdentification/MuGirl/src/MuGirlRecoTool.h b/Reconstruction/MuonIdentification/MuGirl/src/MuGirlRecoTool.h
index 372c889f95aa99059f4ca0b494b5f95d44c68a5d..46b5a98e9f5704c0068c27db18f68d7cda5f2cec 100644
--- a/Reconstruction/MuonIdentification/MuGirl/src/MuGirlRecoTool.h
+++ b/Reconstruction/MuonIdentification/MuGirl/src/MuGirlRecoTool.h
@@ -23,6 +23,7 @@
 
 // Hough transform
 #include "MuonLayerHough/MuonLayerHough.h"
+#include "MuonHoughPatternTools/MuonLayerHoughTool.h"
 
 
 // Forward declarations
@@ -239,6 +240,10 @@ namespace MuGirlNS {
 
     ServiceHandle<MagField::IMagFieldSvc>               m_magFieldSvc;             /**< Magnetic Field Service */
 
+    // Storegate
+    SG::ReadHandleKey<Muon::MuonLayerHoughTool::HoughDataPerSectorVec> m_houghDataPerSectorVecKey {this, 
+        "Key_MuonLayerHoughToolHoughDataPerSectorVec", "HoughDataPerSectorVec", "HoughDataPerSectorVec key"};
+
   }; // class MuGirlRecoTool
 
 
diff --git a/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.cxx b/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.cxx
index 799aa146666be15936af043516c608ad96ffa9ab..0c9ea9f84811d2c786407a89535e218d865b4125 100644
--- a/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.cxx
+++ b/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.cxx
@@ -140,6 +140,7 @@ namespace MuonCombined {
     ATH_CHECK(m_insideOutRecoTool.retrieve());
     ATH_CHECK(m_updator.retrieve());
     ATH_CHECK(m_calibrationDbTool.retrieve());
+    ATH_CHECK(m_houghDataPerSectorVecKey.initialize());
     
     if( m_doTruth ){
       // add pdgs from jobO to set
@@ -1463,15 +1464,22 @@ namespace MuonCombined {
     Muon::MuonStationIndex::DetectorRegionIndex regionIndex = intersection.layerSurface.regionIndex;
     Muon::MuonStationIndex::LayerIndex  layerIndex  = intersection.layerSurface.layerIndex;
 
+    // get hough data
+    SG::ReadHandle<Muon::MuonLayerHoughTool::HoughDataPerSectorVec> houghDataPerSectorVec {m_houghDataPerSectorVecKey};
+    if (!houghDataPerSectorVec.isValid()) {
+      ATH_MSG_ERROR("Hough data per sector vector not found");
+      return;
+    }
+
     // sanity check
-    if( static_cast<int>(m_layerHoughTool->houghData().size()) <= sector-1 ){
-      ATH_MSG_WARNING( " sector " << sector << " larger than the available sectors in the Hough tool: " << m_layerHoughTool->houghData().size() );
+    if( static_cast<int>(houghDataPerSectorVec->size()) <= sector-1 ){
+      ATH_MSG_WARNING( " sector " << sector << " larger than the available sectors in the Hough tool: " << houghDataPerSectorVec->size() );
       return;
     }
 
     // get hough maxima in the layer
     unsigned int sectorLayerHash = Muon::MuonStationIndex::sectorLayerHash( regionIndex,layerIndex );
-    const Muon::MuonLayerHoughTool::HoughDataPerSector& houghDataPerSector = m_layerHoughTool->houghData()[sector-1];
+    const Muon::MuonLayerHoughTool::HoughDataPerSector& houghDataPerSector = (*houghDataPerSectorVec)[sector-1];
 
     // sanity check
     if( houghDataPerSector.maxVec.size() <= sectorLayerHash ){
diff --git a/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.h b/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.h
index 76e7b67ed6b2ec8756e110475560d643a05d1334..68d028fcec635916c982a6f626022a57f0fe1836 100644
--- a/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.h
+++ b/Reconstruction/MuonIdentification/MuonCombinedTrackFindingTools/src/MuonStauRecoTool.h
@@ -231,6 +231,9 @@ namespace MuonCombined {
     void rpcTimeCalibration( const Identifier& id, float& time, float& error ) const;
     void segmentTimeCalibration( const Identifier& id, float& time, float& error ) const;
 
+    /** storegate */
+    SG::ReadHandleKey<Muon::MuonLayerHoughTool::HoughDataPerSectorVec> m_houghDataPerSectorVecKey {this, 
+        "Key_MuonLayerHoughToolHoughDataPerSectorVec", "HoughDataPerSectorVec", "HoughDataPerSectorVec key"};
 
     /** tool handles */
     ToolHandle<Muon::MuonIdHelperTool>               m_idHelper;