diff --git a/Reconstruction/eflowRec/eflowRec/EtaPhiLUT.h b/Reconstruction/eflowRec/eflowRec/EtaPhiLUT.h
new file mode 100644
index 0000000000000000000000000000000000000000..13d2ffd26832a88b317c8ef07db0addbece460fd
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/EtaPhiLUT.h
@@ -0,0 +1,45 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+// EtaPhi Look-up table for eflowRecClusters
+//     author: Teng Jian Khoo
+//     date:   Jun 2020
+//
+// Inspired by IParticlesLookUpTable from ParticlesInConeTools
+// By binning in phi and sorting in eta, extraction of objects in
+// a small eta/phi region can be accelerated.
+// Useful to reduce the number of distance computations needed
+// for e.g. track-cluster matching
+
+#ifndef EFLOWREC_ETAPHILOOKUPTABLE_H
+#define EFLOWREC_ETAPHILOOKUPTABLE_H
+
+#include <vector>
+
+#include "xAODCaloEvent/CaloCluster.h"
+#include "eflowRec/eflowRecCluster.h"
+
+namespace eflowRec {
+
+  /** 2D look up table for eflowRecClusters */
+  class EtaPhiLUT {
+  public:
+    /** constructor taking the desired binsize */
+  EtaPhiLUT( unsigned int nbins = 50 );
+
+  void fill(eflowRecClusterContainer& clustersin);
+
+  /** collect eflowRecClusters in a given cone */
+  std::vector<eflowRecCluster*> clustersInCone( float eta, float phi, float dr ) const;
+
+  private:
+
+    unsigned int m_nphiBins;    /// number of bins
+    float m_phiBinSize;  /// bin size
+    std::vector< std::vector<eflowRecCluster*> > m_phiBinnedLookUpTable; /// the look-up table
+  };
+
+}	// end of namespace
+
+#endif
diff --git a/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h b/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h
index 622bb23af016a153f02cb089ad34d393424c601c..dbbe152db7693fd25eabfab3d664b5389654aea1 100644
--- a/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h
+++ b/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h
@@ -9,11 +9,14 @@
 #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 "xAODCaloEvent/CaloCluster.h"
-#include "xAODTracking/TrackParticle.h"
+#include "eflowRec/PFMatchPositions.h"
+#include "eflowRec/EtaPhiLUT.h"
 
 #include <vector>
 
@@ -24,6 +27,14 @@ class IEFlowCellEOverPTool;
 class PFTrackClusterMatchingTool;
 class eflowRecTrack;
 
+namespace PFMatch {
+  class TrackEtaPhiInFixedLayersProvider;
+}
+
+namespace eflowRec {
+  class EtaPhiLUT;
+}
+
 class PFCellLevelSubtractionTool : public extends<AthAlgTool, IPFSubtractionTool> {
 public:
 
@@ -41,6 +52,7 @@ public:
     eflowCaloObjectContainer* caloObjects;
     eflowRecTrackContainer* tracks;
     eflowRecClusterContainer* clusters;
+    eflowRec::EtaPhiLUT clusterLUT;
   };
 
   void calculateRadialEnergyProfiles(eflowData& data) const;
@@ -53,6 +65,9 @@ public:
   std::string printCluster(const xAOD::CaloCluster* cluster) const;
   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 */
diff --git a/Reconstruction/eflowRec/src/EtaPhiLUT.cxx b/Reconstruction/eflowRec/src/EtaPhiLUT.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..778fb86f788c8f44ea18e1c5fd99f5f7c9fa98b5
--- /dev/null
+++ b/Reconstruction/eflowRec/src/EtaPhiLUT.cxx
@@ -0,0 +1,85 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <cmath>
+#include <algorithm>
+
+#include "eflowRec/EtaPhiLUT.h"
+
+namespace eflowRec {
+
+  /// define 2*Pi
+  static constexpr float TWOPI = 2*M_PI;
+
+  /// hepler function to ensure phi is within +-Pi
+  float phiInRange(float phi) {
+    while (phi >= M_PI) phi -= TWOPI;
+    while (phi < -M_PI) phi += TWOPI;
+    return phi;
+  }
+
+  /// calculate phi index for a given phi
+  unsigned int phiIndex(float phi, float binsize) { return (phi + M_PI)/binsize; }
+
+  EtaPhiLUT::EtaPhiLUT( unsigned int nbins ) :
+    m_nphiBins(nbins),
+    m_phiBinSize(TWOPI/m_nphiBins),
+    m_phiBinnedLookUpTable(m_nphiBins)
+    {}
+
+  void EtaPhiLUT::fill(eflowRecClusterContainer& clustersin)
+  {
+    for(eflowRecCluster* cluster : clustersin) {
+      int index = phiIndex(cluster->getCluster()->phi(), m_phiBinSize);
+      m_phiBinnedLookUpTable[index].push_back(cluster);
+    }
+
+    for( auto& vec : m_phiBinnedLookUpTable ) {
+      std::stable_sort(vec.begin(),vec.end(),[](const eflowRecCluster* c1, const eflowRecCluster* c2) { return c1->getCluster()->eta() < c2->getCluster()->eta(); });
+    }
+  }
+
+  std::vector<eflowRecCluster*> EtaPhiLUT::clustersInCone( float eta, float phi, float dr ) const {
+    std::vector<eflowRecCluster*> result;
+
+    // Indices corresponding to phi range
+    unsigned int iPhiMin = phiIndex( phiInRange(phi-dr), m_phiBinSize );
+    unsigned int iPhiMax = phiIndex( phiInRange(phi+dr), m_phiBinSize );
+
+    // Extract index ranges to iterate over
+    std::vector< std::pair<int,int> > iPhiRanges;
+    if( iPhiMin < iPhiMax ) {
+      iPhiRanges.push_back( std::make_pair(iPhiMin,iPhiMax) );
+    } else { // special treatment for phi-wrapping
+      iPhiRanges.push_back( std::make_pair(0,iPhiMax) );
+      iPhiRanges.push_back( std::make_pair(iPhiMin,m_nphiBins-1) );
+    }
+
+    float dr2Cut = dr*dr;
+
+    // loop over ranges
+    for( auto& range : iPhiRanges ){
+      unsigned int indexMin = range.first;
+      unsigned int indexMax = range.second;
+      for( ; indexMin <= indexMax; ++indexMin ){
+        const std::vector<eflowRecCluster*>& phiClusters = m_phiBinnedLookUpTable[indexMin];
+
+        // Get iterators for clusters in relevant eta range
+        auto it_min = std::lower_bound (phiClusters.begin(), phiClusters.end(), eta-dr, [] (const eflowRecCluster* cl, float etaval) {return cl->getCluster()->eta()<etaval;} );
+        auto it_max = std::upper_bound (it_min, phiClusters.end(), eta+dr, [] (float etaval, const eflowRecCluster* cl) {return etaval<cl->getCluster()->eta();} );
+
+        // Apply deltaR cut
+        for( ;it_min!=it_max;++it_min ){
+          const xAOD::CaloCluster& xcluster = *(*it_min)->getCluster();
+          float deta = eta - xcluster.eta();
+          float dphi = phiInRange(phi - xcluster.phi());
+          float dr2 = deta*deta + dphi*dphi;
+          if( dr2 < dr2Cut ) result.push_back(*it_min);
+        }
+      }
+    }
+    return result;
+  }
+
+} // namespace
diff --git a/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx b/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx
index e43b11d0f6cc75b3fb00435c747365a182d9869c..6461f12d191356abc2edb0bbcdc9ad74e464c396 100644
--- a/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx
+++ b/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx
@@ -44,7 +44,7 @@ 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() {
@@ -72,6 +72,12 @@ StatusCode PFCellLevelSubtractionTool::initialize(){
     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;
 
 }
@@ -84,6 +90,7 @@ void PFCellLevelSubtractionTool::execute(eflowCaloObjectContainer* theEflowCaloO
   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);
@@ -110,6 +117,14 @@ unsigned int PFCellLevelSubtractionTool::matchAndCreateEflowCaloObj(unsigned int
   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);