diff --git a/Reconstruction/eflowRec/eflowRec/PFCalcRadialEnergyProfiles.h b/Reconstruction/eflowRec/eflowRec/PFCalcRadialEnergyProfiles.h
new file mode 100644
index 0000000000000000000000000000000000000000..ff564ce3d4dcdf775cfedf9537c80138b5e33524
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFCalcRadialEnergyProfiles.h
@@ -0,0 +1,23 @@
+/*
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PFCALCRADIALENERGYPROFILES
+#define PFCALCRADIALENERGYPROFILES
+
+#include "AthenaBaseComps/AthMessaging.h"
+#include "eflowRec/PFData.h"
+#include "GaudiKernel/Bootstrap.h"
+#include "GaudiKernel/IMessageSvc.h"
+#include "GaudiKernel/ISvcLocator.h"
+
+class PFCalcRadialEnergyProfiles : public AthMessaging {
+
+  public:
+    PFCalcRadialEnergyProfiles()  : AthMessaging(Gaudi::svcLocator()->service<IMessageSvc>("MessageSvc"),"PFCalcRadialEnergyProfiles ") {};
+    ~PFCalcRadialEnergyProfiles() {};
+
+    void calculate(const PFData& data);
+
+};
+#endif
\ No newline at end of file
diff --git a/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h b/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h
deleted file mode 100644
index a43eb5f223e5fba3862d6e0068ecb0d16322e656..0000000000000000000000000000000000000000
--- a/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h
+++ /dev/null
@@ -1,104 +0,0 @@
-/*
-  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef PFCELLLEVELSUBTRACTIONTOOL_H
-#define PFCELLLEVELSUBTRACTIONTOOL_H
-
-#include "AthenaBaseComps/AthAlgTool.h"
-#include "eflowRec/IPFSubtractionTool.h"
-#include "GaudiKernel/ToolHandle.h"
-
-#include "xAODCaloEvent/CaloCluster.h"
-#include "xAODTracking/TrackParticle.h"
-
-#include "eflowRec/eflowCaloObject.h"
-#include "eflowRec/eflowCellList.h"
-#include "eflowRec/eflowEEtaBinnedParameters.h"
-#include "eflowRec/PFMatchPositions.h"
-#include "eflowRec/EtaPhiLUT.h"
-
-#include <vector>
-
-class eflowCaloObjectContainer;
-class eflowRecTrackContainer;
-class eflowRecClusterContainer;
-class IEFlowCellEOverPTool;
-class PFTrackClusterMatchingTool;
-class eflowRecTrack;
-
-namespace PFMatch {
-  class TrackEtaPhiInFixedLayersProvider;
-}
-
-namespace eflowRec {
-  class EtaPhiLUT;
-}
-
-class PFCellLevelSubtractionTool : public extends<AthAlgTool, IPFSubtractionTool> {
-public:
-
-  PFCellLevelSubtractionTool(const std::string& type,const std::string& name,const IInterface* parent);
-
-  ~PFCellLevelSubtractionTool();
-
-  StatusCode initialize();
-  void execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer* recTrackContainer, eflowRecClusterContainer* recClusterContainer) const;
-  StatusCode finalize();
-
- private:
-
-  struct eflowData {
-    eflowCaloObjectContainer* caloObjects;
-    eflowRecTrackContainer* tracks;
-    eflowRecClusterContainer* clusters;
-    eflowRec::EtaPhiLUT clusterLUT;
-  };
-
-  void calculateRadialEnergyProfiles(eflowData& data) const;
-  void performSubtraction(eflowData& data) const;
-  static bool isEOverPFail(double expectedEnergy, double sigma, double clusterEnergy, bool consistencySigmaCut, bool useGoldenMode) ;
-  bool canAnnihilated(double expectedEnergy, double sigma, double clusterEnergy) const;
-
-  unsigned int matchAndCreateEflowCaloObj(unsigned int n, eflowData& data) const;
-  static std::string printTrack(const xAOD::TrackParticle* track) ;
-  static std::string printCluster(const xAOD::CaloCluster* cluster) ;
-  void printAllClusters(const eflowRecClusterContainer& recClusterContainer) const;
-
-  // Need a track position provider to preselect clusters
-  std::unique_ptr<PFMatch::TrackEtaPhiInFixedLayersProvider> m_trkpos;
-
-  /** Default track-cluster matching tool */
-  ToolHandle<PFTrackClusterMatchingTool> m_matchingTool{this,"PFTrackClusterMatchingTool","PFTrackClusterMatchingTool/CalObjBldMatchingTool","The track-cluster matching tool"};
-  /* Track-cluster matching tools for calculating the pull */
-  ToolHandle<PFTrackClusterMatchingTool> m_matchingToolForPull_015{this,"PFTrackClusterMatchingTool_015","PFTrackClusterMatchingTool/PFPullMatchingTool_015","The 0.15 track-cluster matching tool to calculate the pull"};
-  ToolHandle<PFTrackClusterMatchingTool> m_matchingToolForPull_02{this,"PFTrackClusterMatchingTool_02","PFTrackClusterMatchingTool/PFPullMatchingTool_02","The 0.2 track-cluster matching tool to calculate the pull"};
-  
-  /* Tools for "shower simulation" */
-  std::unique_ptr<eflowEEtaBinnedParameters> m_binnedParameters;
-  ToolHandle<IEFlowCellEOverPTool> m_theEOverPTool{this,"eflowCellEOverPTool","eflowCellEOverPTool","Energy Flow E/P Values and Shower Paremeters Tool"};
-
-  /** Parameter that controls whether to use retain remaining calorimeter energy in track-cluster system, after charged shower subtraction */
-  Gaudi::Property<double> m_subtractionSigmaCut{this,"SubtractionSigmaCut",1.5,"Parameter that controls whether to use retain remaining calorimeter energy in track-cluster system, after charged shower subtraction"};
-  /** Parameter that controls whether a track, in a track-cluster system, will be processed by the split shower recovery algorithm */
-  Gaudi::Property<double> m_consistencySigmaCut{this,"ConsistencySigmaCut",1.0,"Parameter that controls whether a track, in a track-cluster system, will be processed by the split shower recovery algorithm"};
-
-  /** Toggle EOverP algorithm mode, whereby no charged shower subtraction is performed */
-  Gaudi::Property<bool> m_calcEOverP{this,"CalcEOverP",false,"Toggle EOverP algorithm mode, whereby no charged shower subtraction is performed"};
-
-  /** "Number of clusters to match to each track" */
-  Gaudi::Property<int> m_nMatchesInCellLevelSubtraction{this,"nMatchesInCellLevelSubtraction",1,"Number of clusters to match to each track"};
-
-  /** Toggle whether to use golden mode, whereby we only use idealised track-cluster matches within +- N sigma of expected mean e/p */
-  Gaudi::Property<std::string> m_goldenModeString{this,"goldenModeString","","Toggle whether to use golden mode, whereby we only use idealised track-cluster matches within +- N sigma of expected mean e/p"};
-  bool m_runInGoldenMode{false};
-  
-  /** Toggle whether to use updated 2015 charged shower subtraction, which disables the shower subtraction in high calorimeter energy density regions  */
-  Gaudi::Property<bool> m_useUpdated2015ChargedShowerSubtraction{this,"useUpdated2015ChargedShowerSubtraction",true,"Toggle whether to use updated 2015 charged shower subtraction, which disables the shower subtraction in high calorimeter energy density region"};
-
-  /** Toggle whether we have the HLLHC setup */
-  Gaudi::Property<bool> m_isHLLHC{this,"isHLLHC",false,"Toggle whether we have the HLLHC setup"};
-
-};
-
-#endif 
diff --git a/Reconstruction/eflowRec/eflowRec/PFClusterFiller.h b/Reconstruction/eflowRec/eflowRec/PFClusterFiller.h
new file mode 100644
index 0000000000000000000000000000000000000000..e8f7a4af4719871a294061e7e91b5140f1ba495c
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFClusterFiller.h
@@ -0,0 +1,17 @@
+#ifndef PFCLUSTERFILLER_H
+#define PFCLUSTERFILLER_H
+
+#include "eflowRec/eflowRecCluster.h"
+#include "eflowRec/PFData.h"
+
+class PFClusterFiller {
+
+public:
+  PFClusterFiller(){};
+  ~PFClusterFiller(){};
+
+  void fillClustersToRecover(PFData &data) const;
+  void fillClustersToConsider(PFData &data, eflowRecClusterContainer &recClusterContainer) const;
+
+};
+#endif
\ No newline at end of file
diff --git a/Reconstruction/eflowRec/eflowRec/PFData.h b/Reconstruction/eflowRec/eflowRec/PFData.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ecaf0f24d0406b55a88f74a7c0c60294ecec6ab
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFData.h
@@ -0,0 +1,18 @@
+#ifndef PFDATA_H
+#define PFDATA_H
+
+#include "eflowRec/EtaPhiLUT.h"
+
+class eflowCaloObjectContainer;
+class eflowRecTrack;
+class eflowRecCluster;
+
+struct PFData
+  {
+    eflowCaloObjectContainer *caloObjects;
+    std::vector<eflowRecTrack *> tracks;
+    std::vector<eflowRecCluster *> clusters;
+    eflowRec::EtaPhiLUT clusterLUT;
+  };
+
+#endif
\ No newline at end of file
diff --git a/Reconstruction/eflowRec/eflowRec/PFRecoverSplitShowersTool.h b/Reconstruction/eflowRec/eflowRec/PFRecoverSplitShowersTool.h
deleted file mode 100644
index adc97218ee59239b9343ea0386488046976ed6e3..0000000000000000000000000000000000000000
--- a/Reconstruction/eflowRec/eflowRec/PFRecoverSplitShowersTool.h
+++ /dev/null
@@ -1,80 +0,0 @@
-/*
-  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
-*/
-
-#ifndef PFRECOVERSPLITSHOWERSTOOL_H
-#define PFRECOVERSPLITSHOWERSTOOL_H
-
-#include "AthenaBaseComps/AthAlgTool.h"
-
-#include <unordered_set>
-
-#include "GaudiKernel/ToolHandle.h"
-#include "xAODCaloEvent/CaloCluster.h"
-#include "eflowRec/IPFSubtractionTool.h"
-
-class eflowCaloObjectContainer;
-class eflowRecCluster;
-class eflowRecTrack;
-class eflowTrackCaloPoints;
-class eflowEEtaBinnedParameters;
-class IEFlowCellEOverPTool;
-class PFTrackClusterMatchingTool;
-class eflowRingSubtractionManager;
-class eflowRecTrackContainer;
-class eflowRecClusterContainer;
-class eflowCaloObject;
-
-class PFRecoverSplitShowersTool : public extends<AthAlgTool, IPFSubtractionTool> {
-  
- public:
-
-  PFRecoverSplitShowersTool(const std::string& type,const std::string& name,const IInterface* parent);
-
-  ~PFRecoverSplitShowersTool();
-
-  static const InterfaceID& interfaceID();
-
-  virtual StatusCode initialize();
-  void execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer*, eflowRecClusterContainer*) const;
-  virtual StatusCode finalize();
-
- private:
-
-  struct eflowData {
-    eflowCaloObjectContainer* caloObjects;
-    std::unordered_set<eflowRecCluster*> clustersToConsider;
-    std::vector<eflowRecTrack*> tracksToRecover;
-  };
-
-  static void fillClustersToConsider(eflowData&) ;
-  void fillTracksToRecover(eflowData&) const;
-
-  unsigned int matchAndCreateEflowCaloObj(eflowData&) const;
-  void performRecovery(unsigned int const nOriginalObj, eflowData&) const;
-  void subtractTrackFromClusters(const eflowTrackCaloPoints& trackCalo, eflowRingSubtractionManager& ranking, eflowRecTrack* efRecTrack, std::vector<xAOD::CaloCluster*> clusterSubtractionList) const;
-  static double getSumEnergy(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusters) ;
-
-  void printClusterList(std::vector<xAOD::CaloCluster*>& clusters, std::string prefix) const;
-  void performSubtraction(eflowCaloObject* thisEflowCaloObject) const;
-
-  /* Tool for getting e/p values and hadronic shower cell ordering principle parameters */
-  ToolHandle<IEFlowCellEOverPTool> m_theEOverPTool{this,"eflowCellEOverPTool","eflowCellEOverPTool","Energy Flow E/P Values and Shower Parameters Tool"};
-
-  std::unique_ptr<eflowEEtaBinnedParameters> m_binnedParameters;
-
-  /** Parameter that controls whether to use retain remaining calorimeter energy in track-cluster system, after charged shower subtraction */
-  Gaudi::Property<double> m_subtractionSigmaCut{this,"SubtractionSigmaCut",1.5,"Parameter that controls whether to use retain remaining calorimeter energy in track-cluster system, after charged shower subtraction"};
-
-  /** Toggle whether to recover isolated tracks */
-  Gaudi::Property<bool> m_recoverIsolatedTracks{this,"RecoverIsolatedTracks",false,"Toggle whether to recover isolated tracks"};
-  
-  /** Toggle whether to use updated 2015 charged shower subtraction, which disables the shower subtraction in high calorimeter energy density regions  */
-  Gaudi::Property<bool> m_useUpdated2015ChargedShowerSubtraction{this,"useUpdated2015ChargedShowerSubtraction",true,"Toggle whether to use updated 2015 charged shower subtraction, which disables the shower subtraction in high calorimeter energy density region"};
-  
-  /** Toggle whether we have the HLLHC setup */
-  Gaudi::Property<bool> m_isHLLHC{this,"isHLLHC",false,"Toggle whether we have the HLLHC setup"};
-  
-};
-
-#endif
diff --git a/Reconstruction/eflowRec/eflowRec/PFSubtractionTool.h b/Reconstruction/eflowRec/eflowRec/PFSubtractionTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..12608686878631bf19040f34ce17979e6be3eda6
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFSubtractionTool.h
@@ -0,0 +1,91 @@
+/*
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PFSUBTRACTIONTOOL_H
+#define PFSUBTRACTIONTOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "eflowRec/eflowCaloObject.h"
+#include "eflowRec/EtaPhiLUT.h"
+#include "eflowRec/IPFSubtractionTool.h"
+#include "eflowRec/PFData.h"
+#include "eflowRec/PFMatchPositions.h"
+#include "eflowRec/PFTrackClusterMatchingTool.h"
+
+#include "xAODCaloEvent/CaloCluster.h"
+#include "xAODTracking/TrackParticle.h"
+
+
+
+class eflowCaloObjectContainer;
+class eflowEEtaBinnedParameters;
+class eflowRecClusterContainer;
+class eflowRecTrackContainer;
+class IEFlowCellEOverPTool;
+
+class PFSubtractionTool : public extends<AthAlgTool, IPFSubtractionTool>
+{
+
+public:
+  PFSubtractionTool(const std::string &type, const std::string &name, const IInterface *parent);
+  ~PFSubtractionTool();
+
+  StatusCode initialize();
+  void execute(eflowCaloObjectContainer *theEflowCaloObjectContainer, eflowRecTrackContainer *recTrackContainer, eflowRecClusterContainer *recClusterContainer) const;
+  StatusCode finalize();
+
+private:  
+
+  /** This matches ID tracks and CaloClusters, and then creates eflowCaloObjects */
+  unsigned int matchAndCreateEflowCaloObj(PFData &data) const;
+
+  void performSubtraction(const unsigned int& startingPoint,PFData &data) const;
+  void performSubtraction(eflowCaloObject& thisEflowCaloObject) const;
+
+  bool isEOverPFail(double expectedEnergy, double sigma, double clusterEnergy) const;
+
+  bool canAnnihilate(double expectedEnergy, double sigma, double clusterEnergy) const;
+
+  static std::string printTrack(const xAOD::TrackParticle* track);
+  static std::string printCluster(const xAOD::CaloCluster* cluster);
+  void printAllClusters(const eflowRecClusterContainer& recClusterContainer) const;
+
+  /** Tool for getting e/p values and hadronic shower cell ordering principle parameters */
+  ToolHandle<IEFlowCellEOverPTool> m_theEOverPTool{this, "eflowCellEOverPTool", "eflowCellEOverPTool", "Energy Flow E/P Values and Shower Parameters Tool"};
+
+  std::unique_ptr<eflowEEtaBinnedParameters> m_binnedParameters;
+
+  /** Track position provider to be used to preselect clusters */
+  std::unique_ptr<PFMatch::TrackEtaPhiInFixedLayersProvider> m_trkpos;
+
+  /** Default track-cluster matching tool */
+  ToolHandle<PFTrackClusterMatchingTool> m_theMatchingTool{this, "PFTrackClusterMatchingTool", "PFTrackClusterMatchingTool/CalObjBldMatchingTool", "The track-cluster matching tool"};
+
+  /* Track-cluster matching tools for calculating the pull */
+  ToolHandle<PFTrackClusterMatchingTool> m_theMatchingToolForPull_015{this, "PFTrackClusterMatchingTool_015", "PFTrackClusterMatchingTool/PFPullMatchingTool_015", "The 0.15 track-cluster matching tool to calculate the pull"};
+  ToolHandle<PFTrackClusterMatchingTool> m_theMatchingToolForPull_02{this, "PFTrackClusterMatchingTool_02", "PFTrackClusterMatchingTool/PFPullMatchingTool_02", "The 0.2 track-cluster matching tool to calculate the pull"};
+
+  /** Toggle whether we are recovering split showers or not */
+  Gaudi::Property<bool> m_recoverSplitShowers{this,"RecoverSplitShowers",false,"Toggle whether we are recovering split showers or not"};
+
+  /** Number of clusters to match to each track if not doing recover split shower subtraction */
+  Gaudi::Property<int> m_nClusterMatchesToUse{this, "nClusterMatchesToUse", 1, "Number of clusters to match to each track"};
+
+  /** Toggle whether we have the HLLHC setup */
+  Gaudi::Property<bool> m_isHLLHC{this, "isHLLHC", false, "Toggle whether we have the HLLHC setup"};
+
+  /** Toggle EOverP algorithm mode, whereby no charged shower subtraction is performed */
+  Gaudi::Property<bool> m_calcEOverP{this, "CalcEOverP", false, "Toggle EOverP algorithm mode, whereby no charged shower subtraction is performed"};
+
+  /** Parameter that controls whether a track, in a track-cluster system, will be processed by the split shower recovery algorithm */
+  Gaudi::Property<double> m_consistencySigmaCut{this, "ConsistencySigmaCut", 1.0, "Parameter that controls whether a track, in a track-cluster system, will be processed by the split shower recovery algorithm"};
+
+  /** Parameter that controls whether to use retain remaining calorimeter energy in track-cluster system, after charged shower subtraction */
+  Gaudi::Property<double> m_subtractionSigmaCut{this, "SubtractionSigmaCut", 1.5, "Parameter that controls whether to use retain remaining calorimeter energy in track-cluster system, after charged shower subtraction"};
+
+};
+
+#endif
diff --git a/Reconstruction/eflowRec/eflowRec/PFTrackFiller.h b/Reconstruction/eflowRec/eflowRec/PFTrackFiller.h
new file mode 100644
index 0000000000000000000000000000000000000000..63706e29307155e3c7cf8bc379dd08b5987d802b
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFTrackFiller.h
@@ -0,0 +1,17 @@
+#ifndef PFTRACKFILLER_H
+#define PFTRACKFILLER_H
+
+#include "eflowRec/eflowRecTrack.h"
+#include "eflowRec/PFData.h"
+
+class PFTrackFiller {
+
+public:
+  PFTrackFiller(){};
+  ~PFTrackFiller(){};
+
+  void fillTracksToRecover(PFData &data) const;
+  void fillTracksToConsider(PFData &data, eflowRecTrackContainer &recTrackContainer) const;
+
+};
+#endif
\ No newline at end of file
diff --git a/Reconstruction/eflowRec/python/PFCfg.py b/Reconstruction/eflowRec/python/PFCfg.py
index a10d86ad9747d66d0b9e982ebc5082c60f88c57a..99d899d10cdc3056c4bf77438ecf951b591634e5 100644
--- a/Reconstruction/eflowRec/python/PFCfg.py
+++ b/Reconstruction/eflowRec/python/PFCfg.py
@@ -64,7 +64,7 @@ def getPFTrackClusterMatchingTool(inputFlags,matchCut,distanceType,clusterPositi
 
 
 def getPFCellLevelSubtractionTool(inputFlags,toolName):
-    PFCellLevelSubtractionToolFactory = CompFactory.PFCellLevelSubtractionTool
+    PFCellLevelSubtractionToolFactory = CompFactory.PFSubtractionTool
     PFCellLevelSubtractionTool = PFCellLevelSubtractionToolFactory(toolName)
 
     eflowCellEOverPTool_Run2_mc20_JetETMiss = CompFactory.eflowCellEOverPTool_Run2_mc20_JetETMiss
@@ -72,9 +72,9 @@ def getPFCellLevelSubtractionTool(inputFlags,toolName):
 
     if(inputFlags.PF.EOverPMode):
         PFCellLevelSubtractionTool.CalcEOverP = True
-        PFCellLevelSubtractionTool.nMatchesInCellLevelSubtraction = -1
+        PFCellLevelSubtractionTool.nClusterMatchesToUse = -1
     else:
-        PFCellLevelSubtractionTool.nMatchesInCellLevelSubtraction = 1
+        PFCellLevelSubtractionTool.nClusterMatchesToUse = 1
 
     if(inputFlags.PF.EOverPMode):
         PFCellLevelSubtractionTool.PFTrackClusterMatchingTool = getPFTrackClusterMatchingTool(inputFlags,0.2,"EtaPhiSquareDistance","PlainEtaPhi","CalObjBldMatchingTool")
@@ -87,15 +87,13 @@ def getPFCellLevelSubtractionTool(inputFlags,toolName):
     return PFCellLevelSubtractionTool
 
 def getPFRecoverSplitShowersTool(inputFlags,toolName):
-    PFRecoverSplitShowersToolFactory = CompFactory.PFRecoverSplitShowersTool
+    PFRecoverSplitShowersToolFactory = CompFactory.PFSubtractionTool
     PFRecoverSplitShowersTool = PFRecoverSplitShowersToolFactory(toolName)
 
     eflowCellEOverPTool_Run2_mc20_JetETMiss = CompFactory.eflowCellEOverPTool_Run2_mc20_JetETMiss
     PFRecoverSplitShowersTool.eflowCellEOverPTool = eflowCellEOverPTool_Run2_mc20_JetETMiss("eflowCellEOverPTool_Run2_mc20_JetETMiss_Recover")
 
-    PFRecoverSplitShowersTool.RecoverIsolatedTracks = inputFlags.PF.recoverIsolatedTracks
-
-    PFRecoverSplitShowersTool.useUpdated2015ChargedShowerSubtraction = inputFlags.PF.useUpdated2015ChargedShowerSubtraction
+    PFRecoverSplitShowersTool.RecoverSplitShowers = True
 
     return PFRecoverSplitShowersTool
 
diff --git a/Reconstruction/eflowRec/python/PFHLTSequence.py b/Reconstruction/eflowRec/python/PFHLTSequence.py
index da221408c441e8ead37a307322dc6810ec1d224d..846cb75a1b8286236e0acebcfd2ea72afab5cfe9 100644
--- a/Reconstruction/eflowRec/python/PFHLTSequence.py
+++ b/Reconstruction/eflowRec/python/PFHLTSequence.py
@@ -187,7 +187,7 @@ def getPFAlg(flags, clustersin, tracktype):
 
     # Default energy subtraction where a single cluster satisfies the expected
     # track calo energy
-    PFCellLevelSubtractionTool = eflowRecConf.PFCellLevelSubtractionTool(
+    PFCellLevelSubtractionTool = eflowRecConf.PFSubtractionTool(
         "PFCellLevelSubtractionTool",
         eflowCellEOverPTool=CellEOverPTool,
         # Uses a deltaR' cut (deltaR corrected for cluster width in eta/phi) to
@@ -208,8 +208,9 @@ def getPFAlg(flags, clustersin, tracktype):
     # A second cell-level subtraction tool that handles cases where more than one
     # cluster is needed to recover the full track expected energy
     # Reuse the default E/P subtraction tool
-    PFRecoverSplitShowersTool = eflowRecConf.PFRecoverSplitShowersTool(
-        "PFRecoverSplitShowersTool", eflowCellEOverPTool=CellEOverPTool
+    PFRecoverSplitShowersTool = eflowRecConf.PFSubtractionTool(
+        "PFRecoverSplitShowersTool", eflowCellEOverPTool=CellEOverPTool,
+        RecoverSplitShowers = True
     )
 
     # Configure moment calculation using topocluster moment calculator
diff --git a/Reconstruction/eflowRec/src/PFCalcRadialEnergyProfiles.cxx b/Reconstruction/eflowRec/src/PFCalcRadialEnergyProfiles.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4749f4c119c352c485cf38265f8738f8448b12d4
--- /dev/null
+++ b/Reconstruction/eflowRec/src/PFCalcRadialEnergyProfiles.cxx
@@ -0,0 +1,132 @@
+#include "eflowRec/PFCalcRadialEnergyProfiles.h"
+#include "eflowRec/eflowCaloObject.h"
+#include "eflowRec/eflowCellList.h"
+#include "eflowRec/eflowRecCluster.h"
+#include "eflowRec/eflowRecTrack.h"
+#include "eflowRec/eflowRingThicknesses.h"
+#include "eflowRec/eflowSubtractor.h"
+#include "eflowRec/eflowTrackClusterLink.h"
+#include "xAODCaloEvent/CaloCluster.h"
+
+void PFCalcRadialEnergyProfiles::calculate(const PFData& data){
+
+  ATH_MSG_DEBUG("Accessed radial energy profile function");
+
+  for (auto thisEflowCaloObject : *data.caloObjects){
+
+    //If there are no clusters available then we cannot calculate any calorimeter shower profiles
+    if (thisEflowCaloObject->nClusters() < 1 ) continue;
+
+    const std::vector<std::pair<eflowTrackClusterLink*,std::pair<float,float> > > matchedTrackList = thisEflowCaloObject->efRecLink();
+
+    for( auto track: matchedTrackList){
+
+      eflowRecTrack* efRecTrack = (track.first)->getTrack();
+          
+      std::vector<eflowRecCluster*> matchedClusters;
+      matchedClusters.clear();
+      std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
+      for (auto thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
+
+      std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
+      for (auto thisEFlowRecCluster : matchedClusters) clusterSubtractionList.push_back(std::pair(thisEFlowRecCluster->getCluster(),false));
+
+      eflowCellList calorimeterCellList;
+      eflowSubtract::Subtractor::makeOrderedCellList(efRecTrack->getTrackCaloPoints(),clusterSubtractionList,calorimeterCellList);
+      
+      eflowRingThicknesses ringThicknessGenerator;
+      
+      std::vector<int> layerToStoreVector;
+      std::vector<float> radiusToStoreVector;
+      std::vector<float> avgEdensityToStoreVector;
+
+      //Loop over calorimeter layers and in each layer we will calculate radial energy profiles.
+      for (int i=0; i < eflowCalo::nRegions ;i++){
+	
+      	eflowCaloENUM layer = (eflowCaloENUM)i;
+      	ATH_MSG_DEBUG("layer is "<<layer);
+      	double ringThickness = ringThicknessGenerator.ringThickness((eflowCaloENUM)i);
+      	ATH_MSG_DEBUG("ring thickness is "<<ringThickness);
+	
+      	double eta_extr = calorimeterCellList.etaFF(layer);
+      	ATH_MSG_DEBUG("extrapolated eta ["<<layer<<"] is "<<eta_extr);
+      	double phi_extr = calorimeterCellList.phiFF(layer);
+      	ATH_MSG_DEBUG("extrapolated phi ["<<layer<<"] is "<<phi_extr);
+    
+      	if (eta_extr == -999.0){
+      	  continue;
+      	}
+        		
+        //Loop over rings of calorimeter cells going out in increasing radius from the track impact point in this calorimeter layer
+        for (unsigned  int indexOfRing = 0; indexOfRing < 100; indexOfRing++){
+
+          CellIt beginRing = calorimeterCellList.getLowerBound((eflowCaloENUM)i, ringThickness*(indexOfRing));
+          if(beginRing == calorimeterCellList.end()) break;
+
+          int totalCellsinRing = 0;
+          double totalEnergyPerRing = 0;
+      	  double energyDensityPerRing = 0;
+      	  double averageEnergyDensityPerRing = 0;
+
+          //Get the list of calorimeter cells in this cell ring
+      	  std::vector<std::pair<const CaloCell*,int> > tempVector = (*beginRing).second;
+
+          //Loop over the calorimeter cells in this cell ring and calculate total and average energy densities in this cell ring.
+          for (auto thisPair : tempVector){
+      	    const CaloDetDescrElement* DDE = (thisPair.first)->caloDDE();
+	          CaloCell_ID::CaloSample sampling = DDE->getSampling();
+
+            if(eflowCalo::translateSampl(sampling)!=(eflowCaloENUM)i) break;
+
+            ATH_MSG_DEBUG(" cell eta and phi are " << (thisPair.first)->eta() << " and " << (thisPair.first)->phi() << " with index " << thisPair.second << " and sampling of " << sampling);
+	          ATH_MSG_DEBUG(" cell energy is " << (thisPair.first)->energy());
+            
+      	    totalCellsinRing += 1;
+	    
+      	    totalEnergyPerRing += (thisPair.first)->energy();
+      	    double totalEnergyCell = (thisPair.first)->energy();
+      	    ATH_MSG_DEBUG(" Total E per Cell is " << totalEnergyCell);
+      	    ATH_MSG_DEBUG(" Total E per Ring is " << totalEnergyPerRing);
+	    
+      	    double cellVolume = DDE->volume();
+      	    ATH_MSG_DEBUG(" cell volume is " << cellVolume/1000.);
+	    
+      	    double energyDensityCell = totalEnergyCell/(cellVolume/1000.);
+      	    ATH_MSG_DEBUG(" E density per Cell is " << energyDensityCell);
+      	    ATH_MSG_DEBUG(" Initial added E density per Cell is " << energyDensityPerRing);
+      	    energyDensityPerRing += energyDensityCell;
+      	    ATH_MSG_DEBUG(" Final added E density per Cell is " << energyDensityPerRing);
+      	    averageEnergyDensityPerRing = energyDensityPerRing/((totalCellsinRing)*(efRecTrack->getTrack()->e()/1000.));
+
+          }
+
+          ATH_MSG_DEBUG(" track eta is " << efRecTrack->getTrack()->eta());
+      	  ATH_MSG_DEBUG(" track E is " << efRecTrack->getTrack()->e()/1000.);
+      	  ATH_MSG_DEBUG(" Average E density per Ring is " << averageEnergyDensityPerRing);
+	  
+          //Fill up the vectors with energy density, layer and ring radii
+      	  if (averageEnergyDensityPerRing != 0){
+      	    avgEdensityToStoreVector.push_back(averageEnergyDensityPerRing);
+      	    int layerToStore = (eflowCaloENUM)i;
+      	    ATH_MSG_DEBUG("layerToStore is "<< layerToStore);
+      	    layerToStoreVector.push_back(layerToStore);
+      	    double radiusToStore = (indexOfRing)*ringThickness;
+      	    ATH_MSG_DEBUG("radiusToStore is "<< radiusToStore);
+      	    radiusToStoreVector.push_back(radiusToStore);
+      	  }
+      	  else {ATH_MSG_DEBUG("averageEnergyDensityPerRing = 0");}
+
+        }//loop over rings of calorimeter cells
+
+      }//loop over the calorimeter layers
+
+      //Add the vectors with energy density, layer and ring radii to this eflowRecTrack
+      efRecTrack->setLayerCellOrderVector(layerToStoreVector);
+      efRecTrack->setRadiusCellOrderVector(radiusToStoreVector);
+      efRecTrack->setAvgEDensityCellOrderVector(avgEdensityToStoreVector);
+
+    }//loop over tracks matched to clusters
+
+  }//loop on eflowCaloObjectContainer
+
+}
diff --git a/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx b/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx
deleted file mode 100644
index b1e9b66f81a42d998f1dc148175579d84daa4ffe..0000000000000000000000000000000000000000
--- a/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx
+++ /dev/null
@@ -1,518 +0,0 @@
-/*
-  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
-*/
-
-#include "eflowRec/PFCellLevelSubtractionTool.h"
-
-#include "eflowRec/eflowRecCluster.h"
-#include "eflowRec/eflowRecTrack.h"
-#include "eflowRec/eflowTrackClusterLink.h"
-#include "eflowRec/eflowDatabase.h"
-#include "eflowRec/eflowDepthCalculator.h"
-#include "eflowRec/eflowTrackCaloPoints.h"
-#include "eflowRec/IEFlowCellEOverPTool.h"
-#include "eflowRec/eflowEEtaBinnedParameters.h"
-#include "eflowRec/eflowLayerIntegrator.h"
-#include "eflowRec/PFTrackClusterMatchingTool.h"
-#include "eflowRec/eflowRingSubtractionManager.h"
-#include "eflowRec/eflowCellSubtractionFacilitator.h"
-#include "eflowRec/eflowSubtractor.h"
-#include "eflowRec/eflowRingThicknesses.h"
-#include "eflowRec/PFSubtractionStatusSetter.h"
-#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"
-
-#include "CaloEvent/CaloCluster.h"
-
-#include "xAODCaloEvent/CaloCluster.h"
-#include "xAODCaloEvent/CaloClusterKineHelper.h"
-
-#include "xAODPFlow/PFO.h"
-
-#include "eflowRec/eflowCaloObjectMaker.h"
-#include "GaudiKernel/MsgStream.h"
-
-#include "GaudiKernel/ListItem.h"
-
-#include <algorithm>
-#include <iostream>
-#include <cmath>
-#include <vector>
-#include <map>
-
-
-using namespace eflowSubtract;
-
-PFCellLevelSubtractionTool::PFCellLevelSubtractionTool(const std::string& type,const std::string& name,const IInterface* parent) :
-  base_class( type, name, parent),
-  m_binnedParameters(nullptr)
-{
-}
-
-PFCellLevelSubtractionTool::~PFCellLevelSubtractionTool() {
-}
-
-StatusCode PFCellLevelSubtractionTool::initialize(){
-  StatusCode sc;
-
-  if (m_matchingTool.retrieve().isFailure()) {
-    msg(MSG::WARNING) << "Cannot find PFTrackClusterMatchingTool" << endmsg;
-  }
-
-  if (m_matchingToolForPull_02.retrieve().isFailure()) {
-    msg(MSG::WARNING) << "Cannot find PFTrackClusterMatchingTool_2" << endmsg;
-  }
-
-  m_binnedParameters = std::make_unique<eflowEEtaBinnedParameters>();
-
-  sc = m_theEOverPTool->fillBinnedParameters(m_binnedParameters.get());
-
-  m_runInGoldenMode = ((m_goldenModeString == "golden1") || (m_goldenModeString == "golden2"));
-
-  if (sc.isFailure()) {
-    msg(MSG::WARNING) << "Could not execute eflowCellEOverPTool " << endmsg;
-    return StatusCode::SUCCESS;
-  }
-
-  m_trkpos.reset( dynamic_cast<PFMatch::TrackEtaPhiInFixedLayersProvider*>(PFMatch::TrackPositionFactory::Get("EM2EtaPhi").release()) );
-  if(!m_trkpos.get()) {
-    ATH_MSG_ERROR("Failed to get TrackPositionProvider for cluster preselection!");
-    return StatusCode::FAILURE;
-  }
-
-  return StatusCode::SUCCESS;
-
-}
-
-void PFCellLevelSubtractionTool::execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer* recTrackContainer, eflowRecClusterContainer* recClusterContainer) const {
-
-  ATH_MSG_VERBOSE("Executing PFCellLevelSubtractionTool");
-
-  eflowData data;
-  data.caloObjects = theEflowCaloObjectContainer;
-  data.tracks = recTrackContainer;
-  data.clusters = recClusterContainer;
-  data.clusterLUT.fill(*recClusterContainer);
-
-  /* Add each track to its best matching cluster's eflowCaloObject */
-  matchAndCreateEflowCaloObj(m_nMatchesInCellLevelSubtraction, data);
-
-  if (msgLvl(MSG::DEBUG)) printAllClusters(*recClusterContainer);
-
-  /* Check e/p mode - only perform subtraction if not in this mode */
-  if (!m_calcEOverP) {performSubtraction(data);}
-
-  ATH_MSG_DEBUG("Have executed performSubtraction");
-
-  /* Check e/p mode - only perform radial profiles calculations if in this mode */
-  if (m_calcEOverP) {calculateRadialEnergyProfiles(data);}
-
-  ATH_MSG_DEBUG("Have executed calculateRadialEnergyProfiles");
-
-}
-
-unsigned int PFCellLevelSubtractionTool::matchAndCreateEflowCaloObj(unsigned int n, eflowData& data) const {
-  unsigned int nMatches(0);
-
-  const EventContext& ctx = Gaudi::Hive::currentContext();
-
-  /* Create eflowTrackClusterLink after matching */
-  unsigned int nTrack = data.tracks->size();
-  for (unsigned int iTrack=0; iTrack<nTrack; ++iTrack) {
-    eflowRecTrack *thisEfRecTrack = static_cast<eflowRecTrack*>(data.tracks->at(iTrack));
-
-    // Preselect clusters in a local cone
-    eflowRecMatchTrack matchTrack(thisEfRecTrack);
-    PFMatch::EtaPhi etaphi = m_trkpos->getPosition(&matchTrack);
-    float tracketa = etaphi.getEta();
-    float trackphi = etaphi.getPhiD();
-    float dRpresel = 0.4; // Some safety margin, but still << full detector
-    std::vector<eflowRecCluster*> nearbyClusters = data.clusterLUT.clustersInCone(tracketa,trackphi,dRpresel);
-
-     /* Add cluster matches needed for pull calculation*/
-    std::vector<std::pair<eflowRecCluster*,float> > bestClusters_02 = m_matchingToolForPull_02->doMatches(thisEfRecTrack, data.clusters, -1);
-
-    for (auto& matchpair : bestClusters_02) {
-      eflowRecCluster* theCluster = matchpair.first;
-      float distancesq = matchpair.second;
-      eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, theCluster, ctx);
-      if(distancesq<0.15*0.15) {
-	// Narrower cone is a subset of the selected clusters
-	// Distance returned is deltaR^2
-	thisEfRecTrack->addAlternativeClusterMatch(trackClusterLink,"cone_015");
-      }
-      thisEfRecTrack->addAlternativeClusterMatch(trackClusterLink,"cone_02");
-    }
-
-    std::vector<std::pair<eflowRecCluster*,float> > bestClusters = m_matchingTool->doMatches(thisEfRecTrack, data.clusters, n);
-    if (bestClusters.empty()) { continue; }
-
-    /* Matched cluster: create TrackClusterLink and add it to both the track and the cluster (eflowCaloObject will be created later) */
-    nMatches++;
-    for (auto& matchpair : bestClusters) {
-      eflowRecCluster* theCluster = matchpair.first;
-      eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, theCluster, ctx);
-      thisEfRecTrack->addClusterMatch(trackClusterLink);
-      theCluster->addTrackMatch(trackClusterLink);
-    }
-
-  }
-
-  /* Create 3 types eflowCaloObjects: track-only, cluster-only, track-cluster-link */
-  unsigned int nCaloObjects = eflowCaloObjectMaker::makeTrkCluCaloObjects(data.tracks, data.clusters, data.caloObjects);
-  ATH_MSG_DEBUG("PFCellLevelSubtractionTool created total " << nCaloObjects << " CaloObjects.");
-
-  // Should move to a common header or similar, as shared with PFRecoverSplitShowersTool
-  const double gaussianRadius = 0.032;
-  const double gaussianRadiusError = 1.0e-3;
-  const double maximumRadiusSigma = 3.0;
-
-  eflowLayerIntegrator integrator(gaussianRadius, gaussianRadiusError, maximumRadiusSigma, m_isHLLHC);
-
-  /* integrate cells; determine FLI; eoverp */
-  for (unsigned int iCalo=0; iCalo<data.caloObjects->size(); ++iCalo) {
-    eflowCaloObject *thisEflowCaloObject = static_cast<eflowCaloObject*>(data.caloObjects->at(iCalo));
-
-    thisEflowCaloObject->simulateShower(&integrator, m_binnedParameters.get(), m_useUpdated2015ChargedShowerSubtraction);
-  }
-
-  return nMatches;
-}
-
-void PFCellLevelSubtractionTool::calculateRadialEnergyProfiles(eflowData& data) const {
-
-  ATH_MSG_DEBUG("Accessed radial energy profile function");
-
-  unsigned int nEFCaloObs = data.caloObjects->size();
-
-  for (unsigned int iEFCalOb = 0; iEFCalOb < nEFCaloObs; ++iEFCalOb) {
-
-    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iEFCalOb);
-
-    unsigned int nClusters = thisEflowCaloObject->nClusters();
-
-    if (nClusters < 1) {
-    continue;
-    }
-
-    const std::vector<std::pair<eflowTrackClusterLink*,std::pair<float,float> > >& matchedTrackList = thisEflowCaloObject->efRecLink();
-
-    unsigned int nTrackMatches = thisEflowCaloObject->nTracks();
-
-    for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
-
-      eflowRecTrack* efRecTrack = (matchedTrackList[iTrack].first)->getTrack();
-
-      std::vector<eflowRecCluster*> matchedClusters;
-      matchedClusters.clear();
-      std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
-      for (auto *thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
-
-      std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
-      clusterSubtractionList.reserve(matchedClusters.size());
-
-      for (auto *thisEFlowRecCluster : matchedClusters) clusterSubtractionList.emplace_back(thisEFlowRecCluster->getCluster(),false);
-
-      eflowCellList calorimeterCellList;
-      Subtractor::makeOrderedCellList(efRecTrack->getTrackCaloPoints(),clusterSubtractionList,calorimeterCellList);
-      std::vector<int> layerToStoreVector;
-      std::vector<float> radiusToStoreVector;
-      std::vector<float> avgEdensityToStoreVector;
-
-      for (int i=0; i < eflowCalo::nRegions ;i++){
-
-	eflowCaloENUM layer = (eflowCaloENUM)i;
-	ATH_MSG_DEBUG("layer is "<<layer);
-	double ringThickness = eflowRingThicknesses::ringThickness((eflowCaloENUM)i);
-	ATH_MSG_DEBUG("ring thickness is "<<ringThickness);
-
-	double eta_extr = calorimeterCellList.etaFF(layer);
-	ATH_MSG_DEBUG("extrapolated eta ["<<layer<<"] is "<<eta_extr);
-	double phi_extr = calorimeterCellList.phiFF(layer);
-	ATH_MSG_DEBUG("extrapolated phi ["<<layer<<"] is "<<phi_extr);
-
-	if (eta_extr == -999.0){
-	  continue;
-	}
-
-	int indexofRing = 0;
-	int layerToStore = -999;
-
-	double radiusToStore = 0;
-	double totalEnergyPerCell = 0;
-
-	double energyDensityPerCell = -666;
-	double totalEnergyPerRing = 0;
-
-	double cellVolume = 0;
-
-	/* 100 is chosen as a number higher than the number of cells found in a normal list */
-	bool breakloop = false;
-	for (unsigned int n=1; n<100; n++){
-
-	  CellIt beginRing = calorimeterCellList.getLowerBound((eflowCaloENUM)i, ringThickness*(n-1));
-
-	  if(beginRing == calorimeterCellList.end()){
-	    break;
-	  }
-
-	  indexofRing += 1;
-
-	  int totalCellsinRing = 0;
-	  double energyDensityPerRing = 0;
-	  double averageEnergyDensityPerRing = 0;
-
-	  std::vector<std::pair<const CaloCell*,int> > tempVector = (*beginRing).second;
-
-	  for (auto thisPair : tempVector){
-	    const CaloDetDescrElement* DDE = (thisPair.first)->caloDDE();
-	    CaloCell_ID::CaloSample sampling = DDE->getSampling();
-
-	    if(eflowCalo::translateSampl(sampling)!=(eflowCaloENUM)i){
-	      breakloop = true;
-	      break;
-	    }
-
-	    ATH_MSG_DEBUG(" cell eta and phi are " << (thisPair.first)->eta() << " and " << (thisPair.first)->phi() << " with index " << thisPair.second << " and sampling of " << sampling);
-	    ATH_MSG_DEBUG(" cell energy is " << (thisPair.first)->energy());
-
-	    totalCellsinRing += 1;
-
-	    totalEnergyPerRing += (thisPair.first)->energy();
-	    totalEnergyPerCell = (thisPair.first)->energy();
-	    ATH_MSG_DEBUG(" Total E per Cell is " << totalEnergyPerCell);
-	    ATH_MSG_DEBUG(" Total E per Ring is " << totalEnergyPerRing);
-
-	    cellVolume = DDE->volume();
-	    ATH_MSG_DEBUG(" cell volume is " << cellVolume/1000.);
-
-	    energyDensityPerCell = totalEnergyPerCell/(cellVolume/1000.);
-	    ATH_MSG_DEBUG(" E density per Cell is " << energyDensityPerCell);
-	    ATH_MSG_DEBUG(" Initial added E density per Cell is " << energyDensityPerRing);
-	    energyDensityPerRing += energyDensityPerCell;
-	    ATH_MSG_DEBUG(" Final added E density per Cell is " << energyDensityPerRing);
-	    averageEnergyDensityPerRing = energyDensityPerRing/((totalCellsinRing)*(efRecTrack->getTrack()->e()/1000.));
-	  }
-
-	  ATH_MSG_DEBUG(" track eta is " << efRecTrack->getTrack()->eta());
-	  ATH_MSG_DEBUG(" track E is " << efRecTrack->getTrack()->e()/1000.);
-	  ATH_MSG_DEBUG(" Average E density per Ring is " << averageEnergyDensityPerRing);
-
-	  if (averageEnergyDensityPerRing != 0){
-	    avgEdensityToStoreVector.push_back(averageEnergyDensityPerRing);
-	    layerToStore = (eflowCaloENUM)i;
-	    ATH_MSG_DEBUG("layerToStore is "<< layerToStore);
-	    layerToStoreVector.push_back(layerToStore);
-	    radiusToStore = (indexofRing)*ringThickness;
-	    ATH_MSG_DEBUG("radiusToStore is "<< radiusToStore);
-	    radiusToStoreVector.push_back(radiusToStore);
-	  }
-	  else {ATH_MSG_DEBUG("averageEnergyDensityPerRing = 0");}
-	  if (breakloop) break;
-	}//loop on 100 cells
-      }//loop on calo regions
-
-      efRecTrack->setLayerCellOrderVector(layerToStoreVector);
-      efRecTrack->setRadiusCellOrderVector(radiusToStoreVector);
-      efRecTrack->setAvgEDensityCellOrderVector(avgEdensityToStoreVector);
-
-      layerToStoreVector.clear();
-      radiusToStoreVector.clear();
-      avgEdensityToStoreVector.clear();
-    }//loop on track matches
-  }//loop on eflowCaloObjects
-}
-
-void PFCellLevelSubtractionTool::performSubtraction(eflowData& data) const {
-
-  ATH_MSG_DEBUG("In performSubtraction");
-
-  PFSubtractionStatusSetter pfSubtractionStatusSetter;
-  PFSubtractionEnergyRatioCalculator pfSubtractionEnergyRatioCalculator;
-  if( msgLevel( MSG::DEBUG ) ) {
-    pfSubtractionStatusSetter.setLevel(MSG::DEBUG);
-    pfSubtractionEnergyRatioCalculator.setLevel(MSG::DEBUG);
-  }
-
-  unsigned int nEFCaloObs = data.caloObjects->size();
-  for (unsigned int iEFCalOb = 0; iEFCalOb < nEFCaloObs; ++iEFCalOb) {
-    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iEFCalOb);
-
-    unsigned int nClusters = thisEflowCaloObject->nClusters();
-
-    ATH_MSG_DEBUG("Have got an eflowCaloObject with " << nClusters << " clusters");
-
-    if (nClusters < 1) {
-      continue;
-    }
-
-    /* Get matched cluster via Links */
-
-    double expectedEnergy = thisEflowCaloObject->getExpectedEnergy();
-    double clusterEnergy = thisEflowCaloObject->getClusterEnergy();
-    double expectedSigma = sqrt(thisEflowCaloObject->getExpectedVariance());
-
-    ATH_MSG_DEBUG("Have got reference values for eflowCaloObject");
-
-    /* Check e/p */
-    if (isEOverPFail(expectedEnergy, expectedSigma, clusterEnergy, m_consistencySigmaCut,
-                     m_runInGoldenMode)
-        || (m_runInGoldenMode && thisEflowCaloObject->nTracks() > 1)) {
-      continue;
-    }
-
-    ATH_MSG_DEBUG("eflowCaloObject passed eOverP checks");
-
-    /* Do subtraction */
-
-    unsigned int nTrackMatches = thisEflowCaloObject->nTracks();
-    if (nTrackMatches == 0 || m_runInGoldenMode) {
-      continue;
-    }
-
-    const std::vector<std::pair<eflowTrackClusterLink*, std::pair<float,float> > >& matchedTrackList = thisEflowCaloObject->efRecLink();
-
-    if( msgLevel( MSG::DEBUG ) ){ 
-      for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
-	      const xAOD::TrackParticle* thisTrack = (matchedTrackList[iTrack].first)->getTrack()->getTrack();
-        ATH_MSG_DEBUG("eflowCaloObject has track with E, pt and eta " << thisTrack->e() << ", " << thisTrack->pt() << " and " << thisTrack->eta());
-      }
-    }
-
-    ATH_MSG_DEBUG("About to perform subtraction for this eflowCaloObject");
-
-    /* Check if we can annihilate right away */
-    if (canAnnihilated(expectedEnergy, expectedSigma, clusterEnergy)) {
-
-      std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterList;
-      std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
-      unsigned nCluster = thisEflowCaloObject->nClusters();
-      for (unsigned iCluster = 0; iCluster < nCluster; ++iCluster) clusterList.emplace_back(thisEflowCaloObject->efRecCluster(iCluster)->getCluster(),false);
-
-      ATH_MSG_DEBUG("We are going to annihilate. ExpectedEnergy, expectedSigma and clusterEnergy are " << expectedEnergy << ", " << expectedSigma << " and " << clusterEnergy);
-      if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
-
-      pfSubtractionStatusSetter.markAllTracksAnnihStatus(*thisEflowCaloObject);
-      Subtractor::annihilateClusters(clusterList);
-
-      if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterList) ATH_MSG_DEBUG("Have Annihilated cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
-      
-    } else {
-
-      /* Subtract the track from all matched clusters */
-
-      ATH_MSG_DEBUG("Subtract tracks one by one");
-
-      for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
-	      eflowRecTrack* efRecTrack = (matchedTrackList[iTrack].first)->getTrack();
-        unsigned int trackIndex = efRecTrack->getTrack()->index();
-
-	      ATH_MSG_DEBUG("Have got eflowRecTrack number " << iTrack << " for this eflowCaloObject");
-
-	      /* Can't subtract without e/p */
-	      if (!efRecTrack->hasBin()) {
-	        continue;
-	      }
-
-	      if (efRecTrack->isInDenseEnvironment()) continue;
-
-	      ATH_MSG_DEBUG("Have bin and am not in dense environment for this eflowCaloObject");
-
-	      std::vector<eflowRecCluster*> matchedClusters;
-	      matchedClusters.clear();
-	      std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
-	      for (auto *thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
-
-	      ATH_MSG_DEBUG("Have filled matchedClusters list for this eflowCaloObject");
-
-	      std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
-        std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
-	      for (auto *thisEFlowRecCluster : matchedClusters) {
-          xAOD::CaloCluster* thisCluster = thisEFlowRecCluster->getCluster();
-          clusterSubtractionList.emplace_back(thisCluster,false);
-          clusterEnergyMap[thisCluster] = thisCluster->e();
-        }
-
-	      ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");
-
-        bool debugToggle = false;
-        if( msgLevel( MSG::DEBUG ) ) debugToggle = true;
-	      Subtractor::subtractTracksFromClusters(efRecTrack, clusterSubtractionList,debugToggle);
-
-	      ATH_MSG_DEBUG("Have performed subtraction for this eflowCaloObject");
-
-	      /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */
-	      if (canAnnihilated(0, expectedSigma, clusterEnergy)) {
-          if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
-	        Subtractor::annihilateClusters(clusterSubtractionList);
-	        //Now we should mark all of these clusters as being subtracted
-          std::vector<std::pair<float, float> > clusterSubtractedEnergyRatios;
-          pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);           
-          pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject, trackIndex);
-	      }
-        else{
-          std::vector<std::pair<float, float> > clusterSubtractedEnergyRatios;
-          pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatios(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);          
-	        pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject, trackIndex);
-        }
-
-	      ATH_MSG_DEBUG("Have checked if can annihilate clusters for this eflowCaloOject");
-
-      }
-    }
-
-    /* Flag tracks as subtracted */
-    for (unsigned int iTrack = 0; iTrack < thisEflowCaloObject->nTracks(); ++iTrack) {
-      if (!thisEflowCaloObject->efRecTrack(iTrack)->isInDenseEnvironment()) thisEflowCaloObject->efRecTrack(iTrack)->setSubtracted();
-    }
-  }
-}
-
-StatusCode PFCellLevelSubtractionTool::finalize(){
-
-  return StatusCode::SUCCESS;
-}
-
-
-bool PFCellLevelSubtractionTool::isEOverPFail(double expectedEnergy, double sigma, double clusterEnergy, bool consistencySigmaCut, bool useGoldenMode) {
-  if ((expectedEnergy == 0) && (clusterEnergy > 0)) {
-    return false;
-  }
-
-  bool result = useGoldenMode ? fabs(clusterEnergy - expectedEnergy) > consistencySigmaCut*sigma
-                              : clusterEnergy < expectedEnergy - consistencySigmaCut*sigma;
-  return result;
-}
-
-bool PFCellLevelSubtractionTool::canAnnihilated(double expectedEnergy, double sigma, double clusterEnergy) const {
-   return  clusterEnergy - expectedEnergy < m_subtractionSigmaCut * sigma;
-}
-
-std::string PFCellLevelSubtractionTool::printTrack(const xAOD::TrackParticle* track) {
-  std::stringstream result;
-  result << " track with E, eta and phi "<< track->e() << ", " << track->eta() << " and " << track->phi();
-  return result.str();
-}
-
-std::string PFCellLevelSubtractionTool::printCluster(const xAOD::CaloCluster* cluster) {
-  std::stringstream result;
-  result << " cluster with E, eta and phi of " << cluster->e() << ", " << cluster->eta() << " and " << cluster->phi();
-  return result.str();
-}
-
-void PFCellLevelSubtractionTool::printAllClusters(const eflowRecClusterContainer& recClusterContainer) const {
-  unsigned int nClusters = recClusterContainer.size();
-  for (unsigned int iCluster = 0; iCluster < nClusters; ++iCluster) {
-    /* Create the eflowCaloObject and put it in the container */
-    const eflowRecCluster* thisEFRecCluster = recClusterContainer.at(iCluster);
-    if (thisEFRecCluster->getTrackMatches().empty()) {
-      ATH_MSG_DEBUG("Isolated" << printCluster(thisEFRecCluster->getCluster()));
-    } else {
-      ATH_MSG_DEBUG("Matched" << printCluster(thisEFRecCluster->getCluster()));
-      unsigned int nTrackMatches = thisEFRecCluster->getNTracks();
-      std::vector<eflowTrackClusterLink*> theTrackLinks = thisEFRecCluster->getTrackMatches();
-      for (unsigned int iTrackMatch = 0; iTrackMatch < nTrackMatches; ++iTrackMatch) {
-	ATH_MSG_DEBUG("Matched" << printTrack(theTrackLinks[iTrackMatch]->getTrack()->getTrack()));
-      }
-    }
-  }
-}
diff --git a/Reconstruction/eflowRec/src/PFClusterFiller.cxx b/Reconstruction/eflowRec/src/PFClusterFiller.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..43618857a6646466ec13db116ada26ac4efe76fd
--- /dev/null
+++ b/Reconstruction/eflowRec/src/PFClusterFiller.cxx
@@ -0,0 +1,30 @@
+#include "eflowRec/PFClusterFiller.h"
+
+#include "eflowRec/eflowCaloObject.h"
+
+void PFClusterFiller::fillClustersToRecover(PFData &data) const{
+
+  for (auto thisEflowCaloObject : *data.caloObjects )
+  {
+    if (thisEflowCaloObject->nClusters() == 0)continue;
+
+    for (unsigned i = 0; i < thisEflowCaloObject->nClusters(); ++i)
+    {
+      /* Skip empty clusters (subtraction remnants) */
+      const CaloClusterCellLink *theCellLink = thisEflowCaloObject->efRecCluster(i)->getCluster()->getCellLinks();
+      if (0 == (int)theCellLink->size()) continue;
+
+      thisEflowCaloObject->efRecCluster(i)->clearTrackMatches();
+      data.clusters.push_back(thisEflowCaloObject->efRecCluster(i));
+      thisEflowCaloObject->clearClusters();
+    }
+  }
+
+}
+
+void PFClusterFiller::fillClustersToConsider(PFData &data, eflowRecClusterContainer &recClusterContainer) const{
+
+  for (unsigned int count = 0; count < recClusterContainer.size(); count++)
+    data.clusters.push_back(recClusterContainer[count]);
+
+}
\ No newline at end of file
diff --git a/Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx b/Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx
deleted file mode 100644
index bf520bf952dbc55b44227135eae65f2739b294f9..0000000000000000000000000000000000000000
--- a/Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx
+++ /dev/null
@@ -1,307 +0,0 @@
-/*
-  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
-*/
-
-#include "eflowRec/PFRecoverSplitShowersTool.h"
-#include "eflowRec/eflowUtil.h"
-
-#include "eflowRec/eflowRecCluster.h"
-#include "eflowRec/eflowRecTrack.h"
-#include "eflowRec/eflowTrackClusterLink.h"
-#include "eflowRec/eflowCaloObject.h"
-#include "eflowRec/eflowCaloObjectMaker.h"
-#include "eflowRec/eflowTrackCaloPoints.h"
-#include "eflowRec/eflowCellList.h"
-#include "eflowRec/IEFlowCellEOverPTool.h"
-#include "eflowRec/PFTrackClusterMatchingTool.h"
-#include "eflowRec/eflowEEtaBinnedParameters.h"
-#include "eflowRec/eflowLayerIntegrator.h"
-#include "eflowRec/eflowRingSubtractionManager.h"
-#include "eflowRec/eflowCellSubtractionFacilitator.h"
-#include "eflowRec/eflowSubtractor.h"
-#include "eflowRec/PFSubtractionStatusSetter.h"
-#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"
-
-#include "CaloEvent/CaloClusterContainer.h"
-#include "xAODCaloEvent/CaloClusterKineHelper.h"
-
-using namespace eflowSubtract;
-
-PFRecoverSplitShowersTool::PFRecoverSplitShowersTool(const std::string& type,const std::string& name,const IInterface* parent):
-  base_class(type, name, parent),
-  m_binnedParameters(std::make_unique<eflowEEtaBinnedParameters>())
-{
-}
-
-PFRecoverSplitShowersTool::~PFRecoverSplitShowersTool() {}
-
-StatusCode PFRecoverSplitShowersTool::initialize(){
-
-  if (m_theEOverPTool.retrieve().isFailure()){
-    ATH_MSG_WARNING("Cannot find eflowEOverPTool");
-    return StatusCode::SUCCESS;
-  }
-
-  if (m_theEOverPTool->fillBinnedParameters(m_binnedParameters.get()).isFailure()){
-    ATH_MSG_WARNING("Could not execute eflowCellEOverPTool");
-    return StatusCode::SUCCESS;
-  }
-
-  return StatusCode::SUCCESS;
-}
-
-void PFRecoverSplitShowersTool::execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer* eflowRecTracks, eflowRecClusterContainer* eflowRecClusters) const {
-
-  ATH_MSG_DEBUG("Executing");
-
-  eflowData data;
-  data.caloObjects = theEflowCaloObjectContainer;
-  data.tracksToRecover.reserve(eflowRecTracks->size());
-  data.clustersToConsider.reserve(eflowRecClusters->size());
-
-  fillTracksToRecover(data);
-  fillClustersToConsider(data);
-
-  /* Add each track to its matching cluster's eflowCaloObject */
-  int const nOriginalObj = matchAndCreateEflowCaloObj(data);
-
-  /* Extrapolate tracks, match clusters, do shower simulation, run cell subtraction */
-  performRecovery(nOriginalObj,data);
-
-}
-
-StatusCode PFRecoverSplitShowersTool::finalize(){
-  return StatusCode::SUCCESS;
-
-}
-
-void PFRecoverSplitShowersTool::fillClustersToConsider(eflowData& data) {
-
-  data.clustersToConsider.clear();
-
-  for (auto thisEflowCaloObject : *data.caloObjects){
-
-    if (thisEflowCaloObject->nClusters() == 0) { continue; }
-
-    for(unsigned i=0; i<thisEflowCaloObject->nClusters(); ++i){
-        /* Skip empty clusters (subtraction remnants) */
-        const CaloClusterCellLink* theCellLink = thisEflowCaloObject->efRecCluster(i)->getCluster()->getCellLinks();
-        if (0 == (int)theCellLink->size()){ continue; }
-
-        thisEflowCaloObject->efRecCluster(i)->clearTrackMatches();
-        data.clustersToConsider.insert(thisEflowCaloObject->efRecCluster(i));
-        thisEflowCaloObject->clearClusters();
-    }
-  }
-}
-
-void PFRecoverSplitShowersTool::fillTracksToRecover(eflowData& data) const {
-
-  data.tracksToRecover.clear(); 
-  for (auto thisEflowCaloObject : *data.caloObjects){
-
-    /* Skip isolated tracks if flag set */
-    if (!m_recoverIsolatedTracks && thisEflowCaloObject->nClusters() == 0) {
-      unsigned int nTrk = thisEflowCaloObject->nTracks();
-      // But make sure we get eflowObjects from them
-      for (unsigned int iTrk = 0; iTrk < nTrk; ++iTrk) {
-	      eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrk);
-      	if (!thisEfRecTrack->isSubtracted()) thisEfRecTrack->setSubtracted();
-      }
-      continue;
-    }
-
-    /* Add all tracks on the eflowCaloObject that haven't been subtracted yet*/
-    std::vector<eflowRecTrack*> updatedTracks; updatedTracks.clear();
-    unsigned int nTracks = thisEflowCaloObject->nTracks();
-
-    /* For cluster-only CaloObjects */
-    if(nTracks == 0) continue;
-
-    /* For track-only and track-cluster CaloObjects */
-    for (unsigned int iTrack = 0; iTrack < nTracks; ++iTrack){
-      eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrack);
-      if (thisEfRecTrack->isSubtracted()){
-        updatedTracks.push_back(thisEfRecTrack);
-        continue;
-      }
-      thisEfRecTrack->clearClusterMatches();
-      data.tracksToRecover.push_back(thisEfRecTrack);
-    }
-    thisEflowCaloObject->clearTracks();
-    if (!updatedTracks.empty()) {
-      thisEflowCaloObject->addTracks(updatedTracks);
-      updatedTracks.clear();
-    } else {
-      thisEflowCaloObject->clearLinks();
-    }
-  }
-
-  std::sort(data.tracksToRecover.begin(),data.tracksToRecover.end(),eflowRecTrack::SortDescendingPt());
-
-}
-
-unsigned int PFRecoverSplitShowersTool::matchAndCreateEflowCaloObj(eflowData& data) const {
-
-  const EventContext& ctx = Gaudi::Hive::currentContext();
-
-  /* Cache the original number of eflowCaloObjects */
-  const unsigned int nCaloObj = data.caloObjects->size();
-
-  /* loop tracks in data.tracksToRecover and do matching */
-  for (auto *thisEfRecTrack : data.tracksToRecover){
- 
-    /* No point to do anything if bin does not exist */
-    if (!thisEfRecTrack->hasBin()) {
-      std::unique_ptr<eflowCaloObject> thisEflowCaloObject = std::make_unique<eflowCaloObject>();
-      thisEflowCaloObject->addTrack(thisEfRecTrack);
-      data.caloObjects->push_back(std::move(thisEflowCaloObject));
-      continue;
-    }
-    if (msgLvl(MSG::DEBUG)){
-      const xAOD::TrackParticle* track = thisEfRecTrack->getTrack();
-      ATH_MSG_DEBUG("Recovering charged EFO with e,pt, eta and phi " << track->e() << ", " << track->pt() << ", " << track->eta() << " and " << track->phi());
-    }
-    // Get list of matched clusters in the dR<0.2 cone -- already identified
-    const std::vector<eflowTrackClusterLink*>* matchedClusters_02 = thisEfRecTrack->getAlternativeClusterMatches("cone_02");
-    if (!matchedClusters_02) { continue; }
-
-    if (msgLvl(MSG::DEBUG)){
-      for (auto *trkClusLink : *matchedClusters_02) {
-	const eflowRecCluster* thisEFRecCluster = trkClusLink->getCluster();
-	ATH_MSG_DEBUG("Have matched cluster with e, pt, eta, phi of " << thisEFRecCluster->getCluster()->e() << ", " <<  thisEFRecCluster->getCluster()->pt() << ", " << thisEFRecCluster->getCluster()->eta() << " and " << thisEFRecCluster->getCluster()->phi());
-      }
-    }
-
-    if (matchedClusters_02->empty()) { continue; }
-
-    /* Matched cluster: create TrackClusterLink and add it to both the track and the cluster (eflowCaloObject will be created later) */
-    for (auto *trkClusLink : *matchedClusters_02){
-      eflowRecCluster* thisEFRecCluster = trkClusLink->getCluster();
-      // Look up whether this cluster is intended for recovery
-      if( data.clustersToConsider.find(trkClusLink->getCluster()) == data.clustersToConsider.end() ) {continue;}
-      eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack,thisEFRecCluster, ctx);
-      thisEfRecTrack->addClusterMatch(trackClusterLink);
-      thisEFRecCluster->addTrackMatch(trackClusterLink);
-    }
-  }
-
-  /* Create all eflowCaloObjects that contain eflowRecCluster */
-  eflowCaloObjectMaker makeCaloObject;
-  std::vector<eflowRecCluster*> v_clustersToConsider(data.clustersToConsider.begin(),data.clustersToConsider.end());
-  std::sort(v_clustersToConsider.begin(),v_clustersToConsider.end(),eflowRecCluster::SortDescendingPt());
-  unsigned int nCaloObjects = eflowCaloObjectMaker::makeTrkCluCaloObjects(data.tracksToRecover, v_clustersToConsider,
-								   data.caloObjects);
-  ATH_MSG_DEBUG("PFRecoverSplitShowersTool created " << nCaloObjects << " CaloObjects");
-
-  const double gaussianRadius = 0.032;
-  const double gaussianRadiusError = 1.0e-3;
-  const double maximumRadiusSigma = 3.0;
-  
-  // Should move to a common header or similar, as shared with PFCellLevelSubtractionTool
-  eflowLayerIntegrator integrator(gaussianRadius, gaussianRadiusError, maximumRadiusSigma, m_isHLLHC);
-
-  /* integrate cells; determine FLI; eoverp */
-  for (unsigned int iCalo = nCaloObj; iCalo < data.caloObjects->size(); ++iCalo) {
-    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iCalo);
-    thisEflowCaloObject->simulateShower(&integrator, m_binnedParameters.get(), m_useUpdated2015ChargedShowerSubtraction);
-  }
-  return nCaloObj;
-}
-
-void PFRecoverSplitShowersTool::performSubtraction(eflowCaloObject* thisEflowCaloObject) const {
-
-  PFSubtractionStatusSetter pfSubtractionStatusSetter;
-  PFSubtractionEnergyRatioCalculator pfSubtractionEnergyRatioCalculator;
-  if( msgLevel( MSG::DEBUG ) ) {
-    pfSubtractionStatusSetter.setLevel(MSG::DEBUG);
-    pfSubtractionEnergyRatioCalculator.setLevel(MSG::DEBUG);
-  }
-
-  for (unsigned iTrack = 0; iTrack < thisEflowCaloObject->nTracks(); ++iTrack) {
-    eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrack);
-    unsigned int trackIndex = thisEfRecTrack->getTrack()->index();
-
-    ATH_MSG_DEBUG("About to recover track with e, pt, eta and phi of " << thisEfRecTrack->getTrack()->e() << ", " << thisEfRecTrack->getTrack()->pt() << ", " << thisEfRecTrack->getTrack()->eta() << " and "
-    << thisEfRecTrack->getTrack()->eta());
-    /* Get matched cluster via Links */
-    std::vector<eflowRecCluster*> matchedClusters;
-    matchedClusters.clear();
-    std::vector<eflowTrackClusterLink*> links = thisEfRecTrack->getClusterMatches();
-    for ( auto *thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
-    std::sort(matchedClusters.begin(),matchedClusters.end(),eflowRecCluster::SortDescendingPt());
-
-    if (msgLvl(MSG::DEBUG)){
-      for (auto *thisClus : matchedClusters) ATH_MSG_DEBUG("Cluster " << thisClus->getCluster()->index() << " with e,pt, eta and phi of " << thisClus->getCluster()->e() << ", "<< thisClus->getCluster()->pt() << ", " << thisClus->getCluster()->eta() << " and " << thisClus->getCluster()->phi() << " will be subtracted");
-    }
-    /* Do subtraction */
-    std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
-    clusterSubtractionList.reserve(matchedClusters.size());
-    std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
-    for (auto *thisEFlowRecCluster : matchedClusters) {
-      xAOD::CaloCluster* thisCluster = thisEFlowRecCluster->getCluster();
-      clusterSubtractionList.emplace_back(thisCluster,false);
-      clusterEnergyMap[thisCluster] = thisCluster->e();
-    }
-
-    ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");
-
-    if (getSumEnergy(clusterSubtractionList) - thisEfRecTrack->getEExpect() < m_subtractionSigmaCut
-        * sqrt(thisEfRecTrack->getVarEExpect())) {
-      /* Check if we can annihilate right away */
-      if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
-    
-      Subtractor::annihilateClusters(clusterSubtractionList);
-      //Now we should mark all of these clusters as being subtracted
-      //Now need to mark which clusters were modified in the subtraction procedure
-      std::vector<std::pair<float,float> > clusterSubtractedEnergyRatios;
-      pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);       
-      pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
-    } else {
-      /* Subtract the track from all matched clusters */
-      const bool debugToggle = msgLvl(MSG::DEBUG);
-      Subtractor::subtractTracksFromClusters(thisEfRecTrack, clusterSubtractionList, debugToggle);                
-
-      ATH_MSG_DEBUG("Have performed subtraction for this eflowCaloObject");
-
-      /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */      
-      if (getSumEnergy(clusterSubtractionList) < m_subtractionSigmaCut
-          * sqrt(thisEfRecTrack->getVarEExpect())) {
-        if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating remnant cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
-        Subtractor::annihilateClusters(clusterSubtractionList);
-	      //Now we should mark all of these clusters as being subtracted
-        std::vector<std::pair<float,float> > clusterSubtractedEnergyRatios;
-        pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios); 	      
-        pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
-      }
-      else{
-        std::vector<std::pair<float,float> > clusterSubtractedEnergyRatios;        
-        pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatios(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);    
-	      pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
-      } 
-
-    }
-
-    if (msgLvl(MSG::DEBUG)){
-      for (auto *thisClus : matchedClusters) ATH_MSG_DEBUG("Cluster with e,pt, eta and phi of " << thisClus->getCluster()->e() << ", "<< thisClus->getCluster()->pt() << ", " << thisClus->getCluster()->eta() << " and " << thisClus->getCluster()->phi() << " has been subtracted");
-    } 
-
-    /* Flag tracks as subtracted */
-    if (!thisEfRecTrack->isSubtracted()) thisEfRecTrack->setSubtracted();
-  }
-}
-
-void PFRecoverSplitShowersTool::performRecovery(unsigned int const nOriginalObj, eflowData& data) const {
-  unsigned int nEFCaloObs = data.caloObjects->size();
-  for (unsigned int iCalo = nOriginalObj; iCalo < nEFCaloObs; ++iCalo) {
-    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iCalo);
-    performSubtraction(thisEflowCaloObject);
-  }
-
-}
-
-double PFRecoverSplitShowersTool::getSumEnergy(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusters) {
-  double result = 0.0;
-  for (auto thisPair : clusters) result += (thisPair.first)->e();
-  return result;
-}
diff --git a/Reconstruction/eflowRec/src/PFSubtractionTool.cxx b/Reconstruction/eflowRec/src/PFSubtractionTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c9b04a2f76c48c50def12703fde811244297475d
--- /dev/null
+++ b/Reconstruction/eflowRec/src/PFSubtractionTool.cxx
@@ -0,0 +1,447 @@
+#include "eflowRec/PFSubtractionTool.h"
+
+#include "eflowRec/eflowCaloObject.h"
+#include "eflowRec/eflowCaloObjectMaker.h"
+#include "eflowRec/eflowEEtaBinnedParameters.h"
+#include "eflowRec/eflowLayerIntegrator.h"
+#include "eflowRec/eflowRecTrack.h"
+#include "eflowRec/eflowTrackClusterLink.h"
+#include "eflowRec/IEFlowCellEOverPTool.h"
+#include "eflowRec/PFClusterFiller.h"
+#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"
+#include "eflowRec/PFSubtractionStatusSetter.h"
+#include "eflowRec/eflowSubtractor.h"
+#include "eflowRec/PFTrackFiller.h"
+#include "eflowRec/PFCalcRadialEnergyProfiles.h"
+
+using namespace eflowSubtract;
+
+PFSubtractionTool::PFSubtractionTool(const std::string &type, const std::string &name, const IInterface *parent) : base_class(type, name, parent),
+                                                                                                                   m_binnedParameters(std::make_unique<eflowEEtaBinnedParameters>())
+{
+}
+
+PFSubtractionTool::~PFSubtractionTool()
+{
+}
+
+StatusCode PFSubtractionTool::initialize()
+{
+
+  ATH_CHECK(m_theEOverPTool.retrieve());
+
+  ATH_CHECK(m_theEOverPTool->fillBinnedParameters(m_binnedParameters.get()));
+
+  m_trkpos.reset(dynamic_cast<PFMatch::TrackEtaPhiInFixedLayersProvider *>(PFMatch::TrackPositionFactory::Get("EM2EtaPhi").release()));
+  if (!m_trkpos.get())
+  {
+    ATH_MSG_ERROR("Failed to get TrackPositionProvider for cluster preselection!");
+    return StatusCode::FAILURE;
+  }
+
+  //Retrieve track-cluster matching tools
+  ATH_CHECK(m_theMatchingTool.retrieve());
+  ATH_CHECK(m_theMatchingToolForPull_015.retrieve());
+  ATH_CHECK(m_theMatchingToolForPull_02.retrieve());
+
+  return StatusCode::SUCCESS;
+}
+
+void PFSubtractionTool::execute(eflowCaloObjectContainer *theEflowCaloObjectContainer, eflowRecTrackContainer *recTrackContainer, eflowRecClusterContainer *recClusterContainer) const
+{
+
+  ATH_MSG_DEBUG("Executing");
+
+  PFData data;
+  data.caloObjects = theEflowCaloObjectContainer;
+  const PFTrackFiller pfTrackFiller;
+  if (!m_recoverSplitShowers) pfTrackFiller.fillTracksToConsider(data, *recTrackContainer);
+  else pfTrackFiller.fillTracksToRecover(data);
+
+  const PFClusterFiller pfClusterFiller;
+  if (!m_recoverSplitShowers) pfClusterFiller.fillClustersToConsider(data, *recClusterContainer);
+  else pfClusterFiller.fillClustersToRecover(data);
+
+  ATH_MSG_DEBUG("This event has " << data.tracks.size() << " tracks " << data.clusters.size() << " clusters ");
+
+  unsigned int numMatches = matchAndCreateEflowCaloObj(data);
+
+  if (msgLvl(MSG::DEBUG)) printAllClusters(*recClusterContainer);
+
+  if (!m_calcEOverP){
+    if (!m_recoverSplitShowers) performSubtraction(0,data);
+    else performSubtraction(numMatches,data);
+  }
+  else{
+    PFCalcRadialEnergyProfiles pfCalc;
+    pfCalc.calculate(data);
+  }
+
+}
+
+unsigned int PFSubtractionTool::matchAndCreateEflowCaloObj(PFData &data) const{
+
+  //Counts up how many tracks found at least 1 calorimeter cluster matched to it.
+  unsigned int nMatches(0);
+
+  /* Cache the original number of eflowCaloObjects, if there were any */
+  const unsigned int nCaloObj = data.caloObjects->size();
+  const EventContext &ctx = Gaudi::Hive::currentContext();
+
+  /* loop tracks in data.tracks and do matching */
+  for (auto thisEfRecTrack : data.tracks)
+  {
+    /** No point to do anything if e/p reference bin does not exist */
+    if (!thisEfRecTrack->hasBin()) {
+      std::unique_ptr<eflowCaloObject> thisEflowCaloObject = std::make_unique<eflowCaloObject>();
+      thisEflowCaloObject->addTrack(thisEfRecTrack);
+      data.caloObjects->push_back(std::move(thisEflowCaloObject));
+      continue;
+    }
+
+    if (msgLvl(MSG::DEBUG))
+    {
+      const xAOD::TrackParticle *track = thisEfRecTrack->getTrack();
+      ATH_MSG_DEBUG("Matching track with e,pt, eta and phi " << track->e() << ", " << track->pt() << ", " << track->eta() << " and " << track->phi());
+    }
+
+    std::vector<eflowTrackClusterLink*> bestClusters;
+
+    if (!m_recoverSplitShowers){
+      /** Add cluster matches needed for pull calculation (in eflowCaloObject::simulateShowers) which is used to determine whether to run the charged shower subtraction or not.
+      / Clusters in both a cone of 0.15 and 0.2 are needed for this.
+      / The clusters in a cone of 0.2 are also used as the matched cluster list for recover split showers mode.    
+      **/
+      std::vector<std::pair<eflowRecCluster *, float>> bestClusters_02 = m_theMatchingToolForPull_02->doMatches(thisEfRecTrack, data.clusters, -1);
+      for (auto &matchpair : bestClusters_02)
+      {
+        eflowRecCluster *theCluster = matchpair.first;
+        float distancesq = matchpair.second;
+        eflowTrackClusterLink *trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, theCluster, ctx);
+        if (distancesq < 0.15 * 0.15)
+        {
+          // Narrower cone is a subset of the selected clusters
+          // Distance returned is deltaR^2
+          thisEfRecTrack->addAlternativeClusterMatch(trackClusterLink, "cone_015");
+        }
+        thisEfRecTrack->addAlternativeClusterMatch(trackClusterLink, "cone_02");
+      }//loop over bestClusters_02
+
+      //This matching scheme is used to match the calorimeter cluster(s) to be used in the charged showers subtraction for this track.
+      std::vector<std::pair<eflowRecCluster *, float>> matchedClusters = m_theMatchingTool->doMatches(thisEfRecTrack, data.clusters,m_nClusterMatchesToUse);    
+      for (auto thePair : matchedClusters) bestClusters.push_back(eflowTrackClusterLink::getInstance(thisEfRecTrack, thePair.first, ctx));     
+    }
+    else {
+      const std::vector<eflowTrackClusterLink*>* matchedClusters_02 = thisEfRecTrack->getAlternativeClusterMatches("cone_02");
+      if (!matchedClusters_02) continue;
+      else bestClusters = *matchedClusters_02;
+    }
+
+    if (bestClusters.empty()) continue;
+
+    if (msgLvl(MSG::DEBUG))
+    {
+      for (auto thisClusterLink : bestClusters ) {
+        xAOD::CaloCluster* thisCluster = thisClusterLink->getCluster()->getCluster();
+        ATH_MSG_DEBUG("Matched this track to cluster with e,pt, eta and phi " << thisCluster->e() << ", " << thisCluster->pt() << ", " << thisCluster->eta() << " and " << thisCluster->phi());
+      }
+    }
+
+    nMatches++;
+
+    //loop over the matched calorimeter clusters and associate tracks and clusters to each other as needed.
+    for (auto *trkClusLink : bestClusters){
+
+      eflowRecCluster *thisEFRecCluster = trkClusLink->getCluster();
+
+      if (m_recoverSplitShowers){
+        // Look up whether this cluster is intended for recovery
+        if (std::find(data.clusters.begin(), data.clusters.end(), trkClusLink->getCluster()) == data.clusters.end()) continue;       
+      }
+
+      eflowTrackClusterLink *trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, thisEFRecCluster, ctx);
+      thisEfRecTrack->addClusterMatch(trackClusterLink);
+      thisEFRecCluster->addTrackMatch(trackClusterLink);
+    }
+
+  }
+
+  /* Create 3 types eflowCaloObjects: track-only, cluster-only, track-cluster-link */
+  std::vector<eflowRecCluster *> clusters(data.clusters.begin(), data.clusters.end());
+  if (m_recoverSplitShowers) std::sort(clusters.begin(), clusters.end(), eflowRecCluster::SortDescendingPt());
+  unsigned int nCaloObjects = eflowCaloObjectMaker::makeTrkCluCaloObjects(data.tracks, clusters, data.caloObjects);
+  ATH_MSG_DEBUG("Created  " << nCaloObjects << " eflowCaloObjects.");
+  if (msgLvl(MSG::DEBUG)){
+    for (auto thisEFlowCaloObject : *(data.caloObjects)){
+      ATH_MSG_DEBUG("This eflowCaloObject has " << thisEFlowCaloObject->nTracks() << " tracks and " << thisEFlowCaloObject->nClusters() << " clusters ");
+      for (unsigned int count = 0; count < thisEFlowCaloObject->nTracks(); count++){
+        const xAOD::TrackParticle* thisTrack = thisEFlowCaloObject->efRecTrack(count)->getTrack();
+        ATH_MSG_DEBUG("Have track with e, pt, eta and phi of " << thisTrack->e() << ", " << thisTrack->pt() << ", " << thisTrack->eta() << " and " << thisTrack->phi());
+      }
+      for (unsigned int count = 0; count < thisEFlowCaloObject->nClusters(); count++){
+        const xAOD::CaloCluster* thisCluster = thisEFlowCaloObject->efRecCluster(count)->getCluster();
+        ATH_MSG_DEBUG("Have cluster with e, pt, eta and phi of " << thisCluster->e() << ", " << thisCluster->pt() << ", " << thisCluster->eta() << " and " << thisCluster->phi());
+      }
+    }
+  }
+
+  const double gaussianRadius = 0.032;
+  const double gaussianRadiusError = 1.0e-3;
+  const double maximumRadiusSigma = 3.0;
+
+  eflowLayerIntegrator integrator(gaussianRadius, gaussianRadiusError, maximumRadiusSigma, m_isHLLHC);
+
+  /** Start loop from nCaloObj, which should be zero on a first pass */
+  if (!m_recoverSplitShowers && 0 != nCaloObj) ATH_MSG_WARNING("Not in Split Showers Mode and already have " << nCaloObj << " eflowCaloObjects");
+
+  //For each eflowCaloObject we calculate the expected energy deposit in the calorimeter and cell ordering for subtraction.  
+  for (unsigned int iCalo = nCaloObj; iCalo < data.caloObjects->size(); ++iCalo) {  
+    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iCalo);
+    thisEflowCaloObject->simulateShower(&integrator, m_binnedParameters.get(), true);        
+  }
+
+  if (!m_recoverSplitShowers) return nMatches;
+  else return nCaloObj;
+}
+
+void PFSubtractionTool::performSubtraction(const unsigned int& startingPoint,PFData &data) const{
+  unsigned int nEFCaloObs = data.caloObjects->size();
+  for (unsigned int iCalo = startingPoint; iCalo < nEFCaloObs; ++iCalo) {
+    eflowCaloObject* thisEflowCaloObject = data.caloObjects->at(iCalo);
+    this->performSubtraction(*thisEflowCaloObject);
+  }
+}
+
+void PFSubtractionTool::performSubtraction(eflowCaloObject& thisEflowCaloObject) const{
+
+  ATH_MSG_DEBUG("In performSubtraction");
+
+  //These are used to set whether a track was subtracted and also to calculate how much energy was subtracted
+  PFSubtractionStatusSetter pfSubtractionStatusSetter;
+  PFSubtractionEnergyRatioCalculator pfSubtractionEnergyRatioCalculator;
+  if (msgLevel(MSG::DEBUG))
+  {
+    pfSubtractionStatusSetter.setLevel(MSG::DEBUG);
+    pfSubtractionEnergyRatioCalculator.setLevel(MSG::DEBUG);
+  }
+
+  unsigned int nClusters = thisEflowCaloObject.nClusters();
+  unsigned int nTrackMatches = thisEflowCaloObject.nTracks();
+
+  ATH_MSG_DEBUG("Have got an eflowCaloObject with " << nClusters << " clusters and " << nTrackMatches << " track matches");
+
+  if (msgLevel(MSG::DEBUG)){
+    for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack){
+       eflowRecTrack* thisTrack = thisEflowCaloObject.efRecTrack(iTrack);
+       ATH_MSG_DEBUG("eflowCaloObject has track with E, pt and eta " << thisTrack->getTrack()->e() << ", " << thisTrack->getTrack()->pt() << " and " << thisTrack->getTrack()->eta());
+    }
+  }
+  
+  //To keep historical behaviour when in recover split showers mode allow tracks with no cluster matches to proceed.
+  if (!m_recoverSplitShowers && nClusters < 1) return;  
+
+  //Need at least one track in this eflowCaloObject to continue.
+  if (nTrackMatches < 1) return;
+
+  double expectedEnergy = thisEflowCaloObject.getExpectedEnergy();
+  double clusterEnergy = thisEflowCaloObject.getClusterEnergy();
+  double expectedSigma = sqrt(thisEflowCaloObject.getExpectedVariance());
+
+  /* Check e/p, if on first pass - return if e/p not consistent with expected e/p */
+  if (!m_recoverSplitShowers){
+    if (isEOverPFail(expectedEnergy, expectedSigma, clusterEnergy)) return;
+  }
+  
+  const std::vector<std::pair<eflowTrackClusterLink *, std::pair<float, float>>> &matchedTrackList = thisEflowCaloObject.efRecLink();
+
+  ATH_MSG_DEBUG("Matched Track List has size " << matchedTrackList.size());
+
+  if (msgLevel(MSG::DEBUG))
+  {
+    for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack)
+    {
+      const xAOD::TrackParticle *thisTrack = thisEflowCaloObject.efRecTrack(iTrack)->getTrack();      
+      ATH_MSG_DEBUG("eflowCaloObject has track match with E, pt and eta " << thisTrack->e() << ", " << thisTrack->pt() << " and " << thisTrack->eta());
+    }
+  }
+
+  ATH_MSG_DEBUG("About to perform subtraction for this eflowCaloObject");
+
+  bool wasAnnihilated = false;
+
+  //First deal with non-split showers mode
+  if (!m_recoverSplitShowers){
+    /* Check if we can annihilate right away - true if matched cluster has only the expected energy deposit */
+    if (canAnnihilate(expectedEnergy, expectedSigma, clusterEnergy)){
+
+      wasAnnihilated = true;
+
+      std::vector<std::pair<xAOD::CaloCluster *, bool>> clusterList;
+      std::map<xAOD::CaloCluster *, double> clusterEnergyMap;
+      unsigned nCluster = thisEflowCaloObject.nClusters();
+      for (unsigned iCluster = 0; iCluster < nCluster; ++iCluster){
+        clusterList.emplace_back(thisEflowCaloObject.efRecCluster(iCluster)->getCluster(), false);
+      }
+
+      ATH_MSG_DEBUG("We are going to annihilate. ExpectedEnergy, expectedSigma and clusterEnergy are " << expectedEnergy << ", " << expectedSigma << " and " << clusterEnergy);
+      if (msgLevel(MSG::DEBUG))
+        for (auto thisPair : clusterList)
+          ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
+
+      pfSubtractionStatusSetter.markAllTracksAnnihStatus(thisEflowCaloObject);
+      Subtractor::annihilateClusters(clusterList);
+
+      if (msgLevel(MSG::DEBUG))
+        for (auto thisPair : clusterList)
+          ATH_MSG_DEBUG("Have Annihilated cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
+      
+      /* Flag all tracks in this system as subtracted */
+      for (unsigned iTrack = 0; iTrack < thisEflowCaloObject.nTracks(); ++iTrack){
+        eflowRecTrack *thisEfRecTrack = (matchedTrackList[iTrack].first)->getTrack();
+        if (!thisEfRecTrack->isSubtracted()) thisEfRecTrack->setSubtracted();
+      }
+
+    }//if can annihilate this track-cluster systems matched cluster
+  }//split shower recovery mode or regular mode where above annihilation was not triggered
+  if (m_recoverSplitShowers || !wasAnnihilated){
+
+    for (unsigned iTrack = 0; iTrack < thisEflowCaloObject.nTracks(); ++iTrack){
+
+      eflowRecTrack *thisEfRecTrack = thisEflowCaloObject.efRecTrack(iTrack);
+
+      ATH_MSG_DEBUG("About to subtract track with e, pt, eta and phi of " << thisEfRecTrack->getTrack()->e() << ", " << thisEfRecTrack->getTrack()->pt() << ", " << thisEfRecTrack->getTrack()->eta() << " and "
+                                                                       << thisEfRecTrack->getTrack()->eta());
+      
+      if (!thisEfRecTrack->hasBin()) continue;
+
+      ATH_MSG_DEBUG("Have bin for this eflowCaloObject");
+
+      if (thisEfRecTrack->isInDenseEnvironment() && !m_recoverSplitShowers) continue;
+
+      ATH_MSG_DEBUG("Am not in dense environment for this eflowCaloObject");
+
+      /* Get matched cluster via Links */
+      std::vector<eflowRecCluster *> matchedClusters;
+      std::vector<eflowTrackClusterLink *> links = thisEfRecTrack->getClusterMatches();
+      for (auto *thisEFlowTrackClusterLink : links)
+        matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
+      if (m_recoverSplitShowers) std::sort(matchedClusters.begin(), matchedClusters.end(), eflowRecCluster::SortDescendingPt());
+
+      if (msgLvl(MSG::DEBUG)){
+        for (auto *thisClus : matchedClusters)
+          ATH_MSG_DEBUG("Haved matched cluster " << thisClus->getCluster()->index() << " with e,pt, eta and phi of " << thisClus->getCluster()->e() << ", " << thisClus->getCluster()->pt() << ", " << thisClus->getCluster()->eta() << " and " << thisClus->getCluster()->phi() << " will be subtracted");
+      }
+
+      /* Do subtraction */
+      std::vector<std::pair<xAOD::CaloCluster *, bool>> clusterSubtractionList;
+      clusterSubtractionList.reserve(matchedClusters.size());
+      std::map<xAOD::CaloCluster *, double> clusterEnergyMap;
+      for (auto *thisEFlowRecCluster : matchedClusters){
+        xAOD::CaloCluster *thisCluster = thisEFlowRecCluster->getCluster();
+        clusterSubtractionList.emplace_back(thisCluster, false);
+        clusterEnergyMap[thisCluster] = thisCluster->e();
+      }
+
+      ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");
+
+      unsigned int trackIndex = thisEfRecTrack->getTrack()->index();
+
+      //Previously we only checked this in recover split showers, but makes sense to check it in both passes.
+      auto sumClusEnergy = [](double accumulator, std::pair<xAOD::CaloCluster *, bool> thisPair){ return accumulator += thisPair.first->e();};
+      double totalClusterEnergy = std::accumulate(clusterSubtractionList.begin(),clusterSubtractionList.end(),0.0,sumClusEnergy);      
+
+      /* Check if we can annihilate right away - true if matched cluster has only the expected energy deposit */
+      if(canAnnihilate(thisEfRecTrack->getEExpect(),sqrt(thisEfRecTrack->getVarEExpect()),totalClusterEnergy)){
+        
+        if (msgLevel(MSG::DEBUG))
+          for (auto thisPair : clusterSubtractionList)
+            ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
+
+        Subtractor::annihilateClusters(clusterSubtractionList);
+        //Now we should mark all of these clusters as being subtracted
+        //Now need to mark which clusters were modified in the subtraction procedure
+        std::vector<std::pair<float, float>> clusterSubtractedEnergyRatios;
+        pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList, clusterEnergyMap, clusterSubtractedEnergyRatios);
+        pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, thisEflowCaloObject, trackIndex);
+      }
+      else
+      {
+
+        /* Subtract the track from all matched clusters */
+        const bool debugToggle = msgLvl(MSG::DEBUG);
+        Subtractor::subtractTracksFromClusters(thisEfRecTrack, clusterSubtractionList, debugToggle);
+
+        //recalculate total cluster energy from the clusters afer subtraction
+        totalClusterEnergy = std::accumulate(clusterSubtractionList.begin(),clusterSubtractionList.end(),0.0,sumClusEnergy);        
+
+        /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */
+        if (canAnnihilate(0.0,sqrt(thisEfRecTrack->getVarEExpect()), totalClusterEnergy)){
+
+          if (msgLevel(MSG::DEBUG))
+          for (auto thisPair : clusterSubtractionList)
+            ATH_MSG_DEBUG("Annihilating remnant cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
+          Subtractor::annihilateClusters(clusterSubtractionList);
+          //Now we should mark all of these clusters as being subtracted
+          std::vector<std::pair<float, float>> clusterSubtractedEnergyRatios;
+          pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList, clusterEnergyMap, clusterSubtractedEnergyRatios);
+          pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, thisEflowCaloObject, trackIndex);
+        }//if remove the remnant after cell by cell subtraction
+        else
+        { 
+          std::vector<std::pair<float, float>> clusterSubtractedEnergyRatios;
+          pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatios(clusterSubtractionList, clusterEnergyMap, clusterSubtractedEnergyRatios);
+          pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, thisEflowCaloObject, trackIndex);
+        }//if don't remove the remnant after cell by cell subtraction
+
+      }//if not annihilating, and instead subtracting cell by cell
+
+      ATH_MSG_DEBUG("Have subtracted charged shower for this eflowRecTrack");
+
+      /* Flag tracks as subtracted */
+      if (!thisEfRecTrack->isSubtracted()) thisEfRecTrack->setSubtracted();
+
+    }//loop over tracks in eflowCaloObject
+  }//cell by cell subtraction
+
+}
+
+bool PFSubtractionTool::isEOverPFail(double expectedEnergy, double sigma, double clusterEnergy) const
+{
+  if ((expectedEnergy == 0) && (clusterEnergy > 0)) return false;
+  return clusterEnergy < expectedEnergy - m_consistencySigmaCut * sigma;
+}
+
+bool PFSubtractionTool::canAnnihilate(double expectedEnergy, double sigma, double clusterEnergy) const
+{
+  return clusterEnergy - expectedEnergy < m_subtractionSigmaCut * sigma;
+}
+
+std::string PFSubtractionTool::printTrack(const xAOD::TrackParticle* track) {
+  std::stringstream result;
+  result << " track with E, eta and phi "<< track->e() << ", " << track->eta() << " and " << track->phi();
+  return result.str();
+}
+
+std::string PFSubtractionTool::printCluster(const xAOD::CaloCluster* cluster) {
+  std::stringstream result;
+  result << " cluster with E, eta and phi of " << cluster->e() << ", " << cluster->eta() << " and " << cluster->phi();
+  return result.str();
+}
+
+void PFSubtractionTool::printAllClusters(const eflowRecClusterContainer& recClusterContainer) const {
+
+  for ( auto thisEFRecCluster : recClusterContainer){
+    if (thisEFRecCluster->getTrackMatches().empty()) {
+      ATH_MSG_DEBUG("Isolated" << printCluster(thisEFRecCluster->getCluster()));
+    } else {
+      ATH_MSG_DEBUG("Matched" << printCluster(thisEFRecCluster->getCluster()));
+      std::vector<eflowTrackClusterLink*> theTrackLinks = thisEFRecCluster->getTrackMatches();
+      for ( auto thisTrack : theTrackLinks){
+       ATH_MSG_DEBUG("Matched" << printTrack(thisTrack->getTrack()->getTrack()));
+      }
+    }
+  }
+}
+
+
+StatusCode PFSubtractionTool::finalize() { return StatusCode::SUCCESS; }
diff --git a/Reconstruction/eflowRec/src/PFTrackFiller.cxx b/Reconstruction/eflowRec/src/PFTrackFiller.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7b52c6f4fb2acf0a6ffa144341ef30dca3e7659c
--- /dev/null
+++ b/Reconstruction/eflowRec/src/PFTrackFiller.cxx
@@ -0,0 +1,58 @@
+#include "eflowRec/PFTrackFiller.h"
+
+#include "eflowRec/eflowCaloObject.h"
+
+void PFTrackFiller::fillTracksToRecover(PFData &data) const
+{
+  
+  for (auto thisEflowCaloObject : *data.caloObjects){
+
+    unsigned int nTrk = thisEflowCaloObject->nTracks();
+
+    if (0 == thisEflowCaloObject->nClusters()){
+      //If there are no matched clusters we mark each track in thisEflowCaloObject has one that has been processed ("subtracted")
+      //and continue to the next eflowCaloObject
+      for (unsigned int iTrk = 0; iTrk < nTrk; ++iTrk) {
+        eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrk);
+        if (!thisEfRecTrack->isSubtracted()) thisEfRecTrack->setSubtracted();
+      }//loop on tracks on thisEflowCaloObject
+      continue;
+    }//if thisEflowCaloObject has zero CaloCluster     
+
+    //ignore cluster only systems and go to the next eflowCaloObject
+    if(nTrk == 0) continue;
+
+    //Else process systems with only tracks or combinations of tracks and clusters
+
+    //Create a vector to keep a list of tracks we don't want to consider
+    std::vector<eflowRecTrack*> updatedTracks; 
+      
+    for (unsigned int iTrack = 0; iTrack < nTrk; ++iTrack){
+      eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrack);
+      //If this track was already processed ("subtracted") we continue on to the next track
+      if (thisEfRecTrack->isSubtracted()){
+        updatedTracks.push_back(thisEfRecTrack);
+        continue;
+      }
+      //remove existing matched CaloCluster and add this track to our list that we want to consider
+      //for charged hadron subtraction
+      thisEfRecTrack->clearClusterMatches();
+      data.tracks.push_back(thisEfRecTrack);
+    }//loop on tracks in thisEflowCaloObject
+
+    //Now remove all tracks from thisEflowCaloObject
+    thisEflowCaloObject->clearTracks();
+    //Add back the tracks we won't consider for charged hadron subtraction (if any)
+    if (!updatedTracks.empty()) thisEflowCaloObject->addTracks(updatedTracks);
+    //SHOULD WE ALWAYS DO THIS? TO BE TESTED
+    else thisEflowCaloObject->clearLinks();
+  }
+
+  std::sort(data.tracks.begin(),data.tracks.end(),eflowRecTrack::SortDescendingPt());
+
+}
+
+void PFTrackFiller::fillTracksToConsider(PFData &data, eflowRecTrackContainer &recTrackContainer) const
+{
+  for (unsigned int count = 0; count < recTrackContainer.size(); count++) data.tracks.push_back(recTrackContainer[count]);
+}
\ No newline at end of file
diff --git a/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx b/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
index 9f6af8284896c3f185f11431e44ae0b3edabdcc3..25439bda1d4b807b64d214a4af635ebbc100bde7 100644
--- a/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
+++ b/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
@@ -16,8 +16,7 @@
 #include "eflowRec/PFChargedFlowElementCreatorAlgorithm.h"
 #include "eflowRec/PFNeutralFlowElementCreatorAlgorithm.h"
 #include "eflowRec/PFLCNeutralFlowElementCreatorAlgorithm.h"
-#include "eflowRec/PFCellLevelSubtractionTool.h"
-#include "eflowRec/PFRecoverSplitShowersTool.h"
+#include "eflowRec/PFSubtractionTool.h"
 #include "eflowRec/PFMomentCalculatorTool.h"
 #include "eflowRec/PFClusterCollectionTool.h"
 #include "eflowRec/PFLCCalibTool.h"
@@ -41,8 +40,7 @@ DECLARE_COMPONENT( PFNeutralFlowElementCreatorAlgorithm)
 DECLARE_COMPONENT( PFLCNeutralFlowElementCreatorAlgorithm)
 DECLARE_COMPONENT( PFOChargedCreatorAlgorithm )
 DECLARE_COMPONENT( PFONeutralCreatorAlgorithm )
-DECLARE_COMPONENT( PFCellLevelSubtractionTool )
-DECLARE_COMPONENT( PFRecoverSplitShowersTool )
+DECLARE_COMPONENT( PFSubtractionTool )
 DECLARE_COMPONENT( PFMomentCalculatorTool )
 DECLARE_COMPONENT( PFClusterCollectionTool )
 DECLARE_COMPONENT( PFLCCalibTool )