diff --git a/RecoParticleFlow/PFClusterProducer/BuildFile.xml b/RecoParticleFlow/PFClusterProducer/BuildFile.xml
index fcee687e8e160..30ca19dee8af4 100644
--- a/RecoParticleFlow/PFClusterProducer/BuildFile.xml
+++ b/RecoParticleFlow/PFClusterProducer/BuildFile.xml
@@ -12,6 +12,7 @@
 <use   name="Calibration/Tools"/>
 <use   name="CalibCalorimetry/EcalTPGTools"/>
 <use   name="DataFormats/HGCRecHit"/>
+<use   name="RecoLocalCalo/HGCalRecAlgos"/>
 <use   name="vdt_headers"/>
 <use   name="rootmath"/>
 <use   name="root"/>
diff --git a/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticCluster.h b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticCluster.h
new file mode 100644
index 0000000000000..d0e2a78c86cf9
--- /dev/null
+++ b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticCluster.h
@@ -0,0 +1,140 @@
+#ifndef __RecoParticleFlow_PFClusterProducer_RealisticCluster_H__
+#define __RecoParticleFlow_PFClusterProducer_RealisticCluster_H__
+
+#include <array>
+#include <vector>
+#include <algorithm>
+
+class RealisticCluster
+{
+        using Hit3DPosition = std::array<float,3>;
+
+    public:
+
+        // for each SimCluster and for each layer, we store the position of the most energetic hit of the simcluster in the layer
+
+        struct LayerInfo
+        {
+            Hit3DPosition centerOfGravityAtLayer_;
+            Hit3DPosition maxHitPosAtLayer_;
+            float maxEnergyHitAtLayer_;
+        };
+
+        RealisticCluster():
+            totalEnergy(0.f),
+            exclusiveEnergy(0.f),
+            visible(true)
+        {
+
+        }
+
+        void increaseEnergy(float value)
+        {
+            totalEnergy+=value;
+        }
+        void increaseExclusiveEnergy(float value)
+        {
+            exclusiveEnergy+=value;
+        }
+
+        float getExclusiveEnergyFraction() const
+        {
+            float fraction = 0.f;
+            if(totalEnergy>0.f){
+                fraction = exclusiveEnergy/totalEnergy;
+            }
+            return fraction;
+        }
+
+        float getEnergy() const
+        {
+            return totalEnergy;
+        }
+
+        float getExclusiveEnergy() const
+        {
+            return exclusiveEnergy;
+        }
+
+        bool isVisible() const
+        {
+            return visible;
+        }
+
+        void setVisible(bool vis)
+        {
+            visible = vis;
+        }
+
+        void setCenterOfGravity(unsigned int layerId, const Hit3DPosition& position)
+        {
+            layerInfo_[layerId].centerOfGravityAtLayer_ = position;
+        }
+
+        Hit3DPosition getCenterOfGravity(unsigned int layerId) const
+        {
+            return layerInfo_[layerId].centerOfGravityAtLayer_ ;
+        }
+
+        bool setMaxEnergyHit(unsigned int layerId, float newEnergy, const Hit3DPosition position)
+        {
+            if (newEnergy > layerInfo_[layerId].maxEnergyHitAtLayer_)
+            {
+                layerInfo_[layerId].maxEnergyHitAtLayer_ = newEnergy;
+                layerInfo_[layerId].maxHitPosAtLayer_ = position;
+                return true;
+            }
+            else
+                return false;
+        }
+
+        Hit3DPosition getMaxEnergyPosition (unsigned int layerId) const
+        {
+            return layerInfo_[layerId].maxHitPosAtLayer_;
+        }
+
+        float getMaxEnergy(unsigned int layerId) const
+        {
+            return layerInfo_[layerId].maxEnergyHitAtLayer_;
+        }
+
+        void setLayersNum(unsigned int numberOfLayers)
+        {
+            layerInfo_.resize(numberOfLayers);
+        }
+
+        unsigned int getLayersNum() const
+        {
+            return layerInfo_.size();
+        }
+
+        void addHitAndFraction(unsigned int hit, float fraction)
+        {
+            hitIdsAndFractions_.emplace_back(hit,fraction);
+        }
+
+        void modifyFractionForHitId(float fraction, unsigned int hitId)
+        {
+            auto it = std::find_if( hitIdsAndFractions_.begin(), hitIdsAndFractions_.end(),
+                [&hitId](const std::pair<unsigned int, float>& element){ return element.first == hitId;} );
+
+            it->second = fraction;
+        }
+
+        void modifyFractionByIndex(float fraction, unsigned int index)
+        {
+            hitIdsAndFractions_[index].second = fraction;
+        }
+
+        const std::vector< std::pair<unsigned int, float> > & hitsIdsAndFractions() const { return hitIdsAndFractions_; }
+
+    private:
+        std::vector<std::pair<unsigned int, float> > hitIdsAndFractions_;
+        std::vector<LayerInfo> layerInfo_;
+
+        float totalEnergy;
+        float exclusiveEnergy;
+        bool visible;
+};
+
+#endif
diff --git a/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h
new file mode 100644
index 0000000000000..b6d9944ed76d1
--- /dev/null
+++ b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h
@@ -0,0 +1,337 @@
+#ifndef __RecoParticleFlow_PFClusterProducer_RealisticHitToClusterAssociator_H__
+#define __RecoParticleFlow_PFClusterProducer_RealisticHitToClusterAssociator_H__
+/////////////////////////
+// Author: Felice Pantaleo
+// Date:   30/06/2017
+// Email: felice@cern.ch
+/////////////////////////
+#include <vector>
+#include <unordered_map>
+#include "RealisticCluster.h"
+
+namespace
+{
+
+float getDecayLength(unsigned int layer, unsigned int fhOffset, unsigned int bhOffset)
+{
+    constexpr float eeDecayLengthInLayer = 2.f;
+    constexpr float fhDecayLengthInLayer = 1.5f;
+    constexpr float bhDecayLengthInLayer = 1.f;
+
+    if (layer <= fhOffset)
+        return eeDecayLengthInLayer;
+    else if (layer > fhOffset && layer <= bhOffset)
+        return fhDecayLengthInLayer;
+    else
+        return bhDecayLengthInLayer;
+}
+}
+
+
+class RealisticHitToClusterAssociator
+{
+        using Hit3DPosition = std::array<float,3>;
+
+    public:
+
+        struct RealisticHit
+        {
+                struct HitToCluster
+                {
+                        unsigned int simClusterId_;
+                        float mcEnergyFraction_;
+                        float distanceFromMaxHit_;
+                        float realisticEnergyFraction_;
+                };
+
+                Hit3DPosition hitPosition_;
+                float totalEnergy_;
+                unsigned int layerId_;
+                std::vector<HitToCluster> hitToCluster_;
+        };
+
+        void init(std::size_t numberOfHits, std::size_t numberOfSimClusters,
+                std::size_t numberOfLayers)
+        {
+            realisticHits_.resize(numberOfHits);
+            realisticSimClusters_.resize(numberOfSimClusters);
+            for(auto& sc: realisticSimClusters_)
+                sc.setLayersNum(numberOfLayers);
+        }
+
+        void insertHitPosition(float x, float y, float z, unsigned int hitIndex)
+        {
+            realisticHits_[hitIndex].hitPosition_ = {{x,y,z}};
+        }
+
+        void insertLayerId(unsigned int layerId, unsigned int hitIndex)
+        {
+            realisticHits_[hitIndex].layerId_ = layerId;
+        }
+
+        void insertHitEnergy(float energy, unsigned int hitIndex)
+        {
+            realisticHits_[hitIndex].totalEnergy_ = energy;
+        }
+
+        void insertSimClusterIdAndFraction(unsigned int scIdx, float fraction,
+                unsigned int hitIndex, float associatedEnergy)
+        {
+            realisticHits_[hitIndex].hitToCluster_.emplace_back(RealisticHit::HitToCluster{scIdx, fraction, 0.f, 0.f});
+            realisticSimClusters_[scIdx].setMaxEnergyHit(realisticHits_[hitIndex].layerId_, associatedEnergy, realisticHits_[hitIndex].hitPosition_);
+        }
+
+        float XYdistanceFromMaxHit(unsigned int hitId, unsigned int clusterId)
+        {
+            auto l = realisticHits_[hitId].layerId_;
+            const auto& maxHitPosition = realisticSimClusters_[clusterId].getMaxEnergyPosition(l);
+            float distanceSquared = std::pow((realisticHits_[hitId].hitPosition_[0] - maxHitPosition[0]),2) + std::pow((realisticHits_[hitId].hitPosition_[1] - maxHitPosition[1]),2);
+            return std::sqrt(distanceSquared);
+        }
+
+        float XYdistanceFromPointOnSameLayer(unsigned int hitId, const Hit3DPosition& point)
+        {
+            float distanceSquared = std::pow((realisticHits_[hitId].hitPosition_[0] - point[0]),2) + std::pow((realisticHits_[hitId].hitPosition_[1] - point[1]),2);
+            return std::sqrt(distanceSquared);
+        }
+
+        void computeAssociation( float exclusiveFraction, bool useMCFractionsForExclEnergy, unsigned int fhOffset, unsigned int bhOffset)
+        {
+            //if more than exclusiveFraction of a hit's energy belongs to a cluster, that rechit is not counted as shared
+            unsigned int numberOfHits = realisticHits_.size();
+            std::vector<float> partialEnergies;
+            for(unsigned int hitId = 0; hitId < numberOfHits; ++hitId)
+            {
+                partialEnergies.clear();
+                std::vector<unsigned int> removeAssociation;
+                auto& realisticHit = realisticHits_[hitId];
+                unsigned int numberOfClusters = realisticHit.hitToCluster_.size();
+                if(numberOfClusters == 1)
+                {
+                    unsigned int simClusterId = realisticHit.hitToCluster_[0].simClusterId_;
+                    float assignedFraction = 1.f;
+                    realisticHit.hitToCluster_[0].realisticEnergyFraction_ = assignedFraction;
+                    float assignedEnergy = realisticHit.totalEnergy_;
+                    realisticSimClusters_[simClusterId].increaseEnergy(assignedEnergy);
+                    realisticSimClusters_[simClusterId].addHitAndFraction(hitId, assignedFraction);
+                    realisticSimClusters_[simClusterId].increaseExclusiveEnergy(assignedEnergy);
+                }
+                else
+                {
+                    partialEnergies.resize(numberOfClusters,0.f);
+                    unsigned int layer = realisticHit.layerId_;
+                    float sumE = 0.f;
+                    float energyDecayLength = getDecayLength(layer, fhOffset, bhOffset);
+                    for(unsigned int clId = 0; clId < numberOfClusters; ++clId )
+                    {
+                        auto simClusterId = realisticHit.hitToCluster_[clId].simClusterId_;
+                        realisticHit.hitToCluster_[clId].distanceFromMaxHit_ = XYdistanceFromMaxHit(hitId,simClusterId);
+                        // partial energy is computed based on the distance from the maximum energy hit and its energy
+                        // partial energy is only needed to compute a fraction and it's not the energy assigned to the cluster
+                        auto maxEnergyOnLayer = realisticSimClusters_[simClusterId].getMaxEnergy(layer);
+                        if(maxEnergyOnLayer>0.f)
+                        {
+                            partialEnergies[clId] = maxEnergyOnLayer * std::exp(-realisticHit.hitToCluster_[clId].distanceFromMaxHit_/energyDecayLength);
+                        }
+                        sumE += partialEnergies[clId];
+                    }
+                    if(sumE > 0.f)
+                    {
+                        float invSumE = 1.f/sumE;
+                        for(unsigned int clId = 0; clId < numberOfClusters; ++clId )
+                        {
+                            unsigned int simClusterIndex = realisticHit.hitToCluster_[clId].simClusterId_;
+                            float assignedFraction = partialEnergies[clId]*invSumE;
+                            if(assignedFraction > 1e-3)
+                            {
+                                realisticHit.hitToCluster_[clId].realisticEnergyFraction_ = assignedFraction;
+                                float assignedEnergy = assignedFraction *realisticHit.totalEnergy_;
+                                realisticSimClusters_[simClusterIndex].increaseEnergy(assignedEnergy);
+                                realisticSimClusters_[simClusterIndex].addHitAndFraction(hitId, assignedFraction);
+                                // if the hits energy belongs for more than exclusiveFraction to a cluster, also the cluster's
+                                // exclusive energy is increased. The exclusive energy will be needed to evaluate if
+                                // a realistic cluster will be invisible, i.e. absorbed by other clusters
+                                if( (useMCFractionsForExclEnergy and realisticHit.hitToCluster_[clId].mcEnergyFraction_ > exclusiveFraction) or
+                                        (!useMCFractionsForExclEnergy and assignedFraction > exclusiveFraction) )
+                                {
+                                    realisticSimClusters_[simClusterIndex].increaseExclusiveEnergy(assignedEnergy);
+                                }
+                            }
+                            else
+                            {
+                                removeAssociation.push_back(simClusterIndex);
+                            }
+                        }
+                    }
+                }
+
+                while(!removeAssociation.empty())
+                {
+                    auto clusterToRemove = removeAssociation.back();
+                    removeAssociation.pop_back();
+
+                    realisticHit.hitToCluster_.erase(std::remove_if(realisticHit.hitToCluster_.begin(), realisticHit.hitToCluster_.end(), [clusterToRemove](const RealisticHit::HitToCluster& x)
+                    {
+                        return x.simClusterId_ == clusterToRemove;
+                    }), realisticHit.hitToCluster_.end());
+                }
+            }
+        }
+
+        void findAndMergeInvisibleClusters(float invisibleFraction, float exclusiveFraction)
+        {
+            unsigned int numberOfRealSimClusters = realisticSimClusters_.size();
+
+            for(unsigned int clId= 0; clId < numberOfRealSimClusters; ++clId)
+            {
+                if(realisticSimClusters_[clId].getExclusiveEnergyFraction() < invisibleFraction)
+                {
+                    realisticSimClusters_[clId].setVisible(false);
+                    auto& hAndF = realisticSimClusters_[clId].hitsIdsAndFractions();
+                    std::unordered_map < unsigned int, float> energyInNeighbors;
+                    float totalSharedEnergy=0.f;
+
+                    for(auto& elt : hAndF)
+                    {
+                        unsigned int hitId = elt.first;
+                        float fraction = elt.second;
+                        auto& realisticHit = realisticHits_[hitId];
+
+                        if(realisticHit.hitToCluster_.size() >1 && fraction < 1.f)
+                        {
+                            float correction = 1.f - fraction;
+                            unsigned int numberOfClusters = realisticHit.hitToCluster_.size();
+                            int clusterToRemove = -1;
+                            for(unsigned int i = 0; i< numberOfClusters; ++i)
+                            {
+                                auto simClusterIndex = realisticHit.hitToCluster_[i].simClusterId_;
+                                if(simClusterIndex == clId)
+                                {
+                                    clusterToRemove = i;
+                                }
+                                else
+                                if(realisticSimClusters_[simClusterIndex].isVisible())
+                                {
+                                    float oldFraction = realisticHit.hitToCluster_[i].realisticEnergyFraction_;
+                                    float newFraction = oldFraction/correction;
+                                    float oldEnergy = oldFraction*realisticHit.totalEnergy_;
+                                    float newEnergy= newFraction*realisticHit.totalEnergy_;
+                                    float sharedEnergy = newEnergy-oldEnergy;
+                                    energyInNeighbors[simClusterIndex] +=sharedEnergy;
+                                    totalSharedEnergy +=sharedEnergy;
+                                    realisticSimClusters_[simClusterIndex].increaseEnergy(sharedEnergy);
+                                    realisticSimClusters_[simClusterIndex].modifyFractionForHitId(newFraction, hitId);
+                                    realisticHit.hitToCluster_[i].realisticEnergyFraction_ = newFraction;
+                                    if(newFraction > exclusiveFraction)
+                                    {
+                                        realisticSimClusters_[simClusterIndex].increaseExclusiveEnergy(sharedEnergy);
+                                        if(oldFraction <=exclusiveFraction)
+                                        {
+                                            realisticSimClusters_[simClusterIndex].increaseExclusiveEnergy(oldEnergy);
+                                        }
+                                    }
+                                }
+                            }
+                            realisticSimClusters_[realisticHit.hitToCluster_[clusterToRemove].simClusterId_].modifyFractionForHitId(0.f, hitId);
+                            realisticHit.hitToCluster_.erase(realisticHit.hitToCluster_.begin()+clusterToRemove);
+                        }
+                    }
+
+                    for(auto& elt : hAndF)
+                    {
+                        unsigned int hitId = elt.first;
+                        auto& realisticHit = realisticHits_[hitId];
+                        if(realisticHit.hitToCluster_.size()==1 and realisticHit.hitToCluster_[0].simClusterId_ ==  clId and totalSharedEnergy > 0.f)
+                        {
+                            for (auto& pair: energyInNeighbors)
+                            {
+                                // hits that belonged completely to the absorbed cluster are redistributed
+                                // based on the fraction of energy shared in the shared hits
+                                float sharedFraction = pair.second/totalSharedEnergy;
+                                float assignedEnergy = realisticHit.totalEnergy_*sharedFraction;
+                                realisticSimClusters_[pair.first].increaseEnergy(assignedEnergy);
+                                realisticSimClusters_[pair.first].addHitAndFraction(hitId, sharedFraction);
+                                realisticHit.hitToCluster_.emplace_back(RealisticHit::HitToCluster{pair.first, 0.f, -1.f, sharedFraction});
+                                if(sharedFraction > exclusiveFraction)
+                                    realisticSimClusters_[pair.first].increaseExclusiveEnergy(assignedEnergy);
+                            }
+                        }
+                    }
+
+                }
+            }
+        }
+
+        void findCentersOfGravity()
+        {
+            for(auto& cluster : realisticSimClusters_)
+            {
+                if(cluster.isVisible())
+                {
+                    unsigned int layersNum = cluster.getLayersNum();
+                    std::vector<float> totalEnergyPerLayer(layersNum, 0.f);
+                    std::vector<float> xEnergyPerLayer(layersNum, 0.f);
+                    std::vector<float> yEnergyPerLayer(layersNum, 0.f);
+                    std::vector<float> zPositionPerLayer(layersNum, 0.f);
+                    const auto& hAndF = cluster.hitsIdsAndFractions();
+                    for(auto& elt : hAndF)
+                    {
+                        auto hitId = elt.first;
+                        auto fraction = elt.second;
+                        const auto & hit = realisticHits_[hitId];
+                        const auto & hitPos = hit.hitPosition_;
+                        auto layerId = hit.layerId_;
+                        auto hitEinCluster = hit.totalEnergy_*fraction;
+                        totalEnergyPerLayer[layerId]+= hitEinCluster;
+                        xEnergyPerLayer[layerId] += hitPos[0]*hitEinCluster;
+                        yEnergyPerLayer[layerId] += hitPos[1]*hitEinCluster;
+                        zPositionPerLayer[layerId] = hitPos[2];
+                    }
+                    Hit3DPosition centerOfGravity;
+                    for(unsigned int layerId=0; layerId<layersNum; layerId++)
+                    {
+                        auto energyOnLayer = totalEnergyPerLayer[layerId];
+                        if(energyOnLayer > 0.f)
+                        {
+                            centerOfGravity = {{xEnergyPerLayer[layerId]/energyOnLayer,yEnergyPerLayer[layerId]/energyOnLayer, zPositionPerLayer[layerId] }};
+                            cluster.setCenterOfGravity(layerId,centerOfGravity );
+                        }
+                    }
+                }
+            }
+        }
+
+        void filterHitsByDistance(float maxDistance)
+        {
+            for(auto& cluster : realisticSimClusters_)
+            {
+                if(cluster.isVisible())
+                {
+                    auto& hAndF = cluster.hitsIdsAndFractions();
+                    for(unsigned int i = 0; i<hAndF.size(); ++i)
+                    {
+                        auto hitId = hAndF[i].first;
+                        const auto & hit = realisticHits_[hitId];
+                        auto layerId = hit.layerId_;
+                        if(XYdistanceFromPointOnSameLayer(hitId, cluster.getCenterOfGravity(layerId)) > maxDistance)
+                        {
+                            cluster.modifyFractionByIndex(0.f, i);
+                        }
+                    }
+                }
+            }
+        }
+
+        const std::vector< RealisticCluster > & realisticClusters() const
+        {   return realisticSimClusters_;}
+
+    private:
+        // the vector of the Realistic SimClusters
+        std::vector< RealisticCluster > realisticSimClusters_;
+        // the vector of the Realistic SimClusters
+        std::vector< RealisticHit > realisticHits_;
+    };
+
+
+
+#endif
diff --git a/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticSimClusterMapper.cc b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticSimClusterMapper.cc
new file mode 100644
index 0000000000000..0ddb4f8ec4315
--- /dev/null
+++ b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticSimClusterMapper.cc
@@ -0,0 +1,128 @@
+/////////////////////////
+// Author: Felice Pantaleo
+// Date:   30/06/2017
+// Email: felice@cern.ch
+/////////////////////////
+#include <unordered_map>
+#include "RealisticSimClusterMapper.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "FWCore/Framework/interface/Event.h"
+
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
+#include "RealisticHitToClusterAssociator.h"
+
+#ifdef PFLOW_DEBUG
+#define LOGVERB(x) edm::LogVerbatim(x)
+#define LOGWARN(x) edm::LogWarning(x)
+#define LOGERR(x) edm::LogError(x)
+#define LOGDRESSED(x) edm::LogInfo(x)
+#else
+#define LOGVERB(x) LogTrace(x)
+#define LOGWARN(x) edm::LogWarning(x)
+#define LOGERR(x) edm::LogError(x)
+#define LOGDRESSED(x) LogDebug(x)
+#endif
+
+void RealisticSimClusterMapper::updateEvent(const edm::Event& ev)
+{
+    ev.getByToken(simClusterToken_, simClusterH_);
+}
+
+void RealisticSimClusterMapper::update(const edm::EventSetup& es)
+{
+    rhtools_.getEventSetup(es);
+
+}
+
+void RealisticSimClusterMapper::buildClusters(const edm::Handle<reco::PFRecHitCollection>& input,
+        const std::vector<bool>& rechitMask, const std::vector<bool>& seedable,
+        reco::PFClusterCollection& output)
+{
+    const SimClusterCollection& simClusters = *simClusterH_;
+    auto const& hits = *input;
+    RealisticHitToClusterAssociator realisticAssociator;
+    const int numberOfLayers = rhtools_.getLayer(ForwardSubdetector::ForwardEmpty);
+    realisticAssociator.init(hits.size(), simClusters.size(), numberOfLayers + 1);
+    // for quick indexing back to hit energy
+    std::unordered_map < uint32_t, size_t > detIdToIndex(hits.size());
+    for (uint32_t i = 0; i < hits.size(); ++i)
+    {
+        detIdToIndex[hits[i].detId()] = i;
+        auto ref = makeRefhit(input, i);
+        const auto& hitPos = rhtools_.getPosition(ref->detId());
+
+        realisticAssociator.insertHitPosition(hitPos.x(), hitPos.y(), hitPos.z(), i);
+        realisticAssociator.insertHitEnergy(ref->energy(), i);
+        realisticAssociator.insertLayerId(rhtools_.getLayerWithOffset(ref->detId()), i);
+
+    }
+
+    for (unsigned int ic = 0; ic < simClusters.size(); ++ic)
+    {
+        const auto & sc = simClusters[ic];
+        const auto& hitsAndFractions = sc.hits_and_fractions();
+        for (const auto& hAndF : hitsAndFractions)
+        {
+            auto itr = detIdToIndex.find(hAndF.first);
+            if (itr == detIdToIndex.end())
+            {
+                continue; // hit wasn't saved in reco or did not pass the SNR threshold
+            }
+            auto hitId = itr->second;
+            auto ref = makeRefhit(input, hitId);
+            float fraction = hAndF.second;
+            float associatedEnergy = fraction * ref->energy();
+            realisticAssociator.insertSimClusterIdAndFraction(ic, fraction, hitId,
+                    associatedEnergy);
+        }
+    }
+    realisticAssociator.computeAssociation(exclusiveFraction_, useMCFractionsForExclEnergy_,
+            rhtools_.lastLayerEE(), rhtools_.lastLayerFH());
+    realisticAssociator.findAndMergeInvisibleClusters(invisibleFraction_, exclusiveFraction_);
+    realisticAssociator.findCentersOfGravity();
+    if(maxDistanceFilter_)
+        realisticAssociator.filterHitsByDistance(maxDistance_);
+
+    const auto& realisticClusters = realisticAssociator.realisticClusters();
+    unsigned int nClusters = realisticClusters.size();
+    for (unsigned ic = 0; ic < nClusters; ++ic)
+    {
+        float highest_energy = 0.0f;
+        output.emplace_back();
+        reco::PFCluster& back = output.back();
+        edm::Ref < std::vector<reco::PFRecHit> > seed;
+        if (realisticClusters[ic].isVisible())
+        {
+            const auto& hitsIdsAndFractions = realisticClusters[ic].hitsIdsAndFractions();
+            for (const auto& idAndF : hitsIdsAndFractions)
+            {
+                auto fraction = idAndF.second;
+                if (fraction > 0.f)
+                {
+                    auto ref = makeRefhit(input, idAndF.first);
+                    back.addRecHitFraction(reco::PFRecHitFraction(ref, fraction));
+                    const float hit_energy = fraction * ref->energy();
+                    if (hit_energy > highest_energy || highest_energy == 0.0)
+                    {
+                        highest_energy = hit_energy;
+                        seed = ref;
+                    }
+                }
+            }
+        }
+        if (back.hitsAndFractions().size() != 0)
+        {
+            back.setSeed(seed->detId());
+            back.setEnergy(realisticClusters[ic].getEnergy());
+            back.setCorrectedEnergy(realisticClusters[ic].getEnergy());
+        }
+        else
+        {
+            back.setSeed(-1);
+            back.setEnergy(0.f);
+        }
+    }
+}
+
diff --git a/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticSimClusterMapper.h b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticSimClusterMapper.h
new file mode 100644
index 0000000000000..ce8024ae8b232
--- /dev/null
+++ b/RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticSimClusterMapper.h
@@ -0,0 +1,55 @@
+#ifndef __RecoParticleFlow_PFClusterProducer_RealisticSimClusterMapper_H__
+#define __RecoParticleFlow_PFClusterProducer_RealisticSimClusterMapper_H__
+/////////////////////////
+// Author: Felice Pantaleo
+// Date:   30/06/2017
+// Email: felice@cern.ch
+/////////////////////////
+#include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
+#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+
+class RealisticSimClusterMapper : public InitialClusteringStepBase {
+ public:
+ RealisticSimClusterMapper(const edm::ParameterSet& conf,
+			 edm::ConsumesCollector& sumes) :
+    InitialClusteringStepBase(conf,sumes),
+    invisibleFraction_(conf.getParameter<double>("invisibleFraction")),
+    exclusiveFraction_(conf.getParameter<double>("exclusiveFraction")),
+    maxDistanceFilter_(conf.getParameter<bool>("maxDistanceFilter")),
+    maxDistance_(conf.getParameter<double>("maxDistance")),
+    useMCFractionsForExclEnergy_(conf.getParameter<bool>("useMCFractionsForExclEnergy"))
+    {
+      simClusterToken_ = sumes.consumes<SimClusterCollection>(conf.getParameter<edm::InputTag>("simClusterSrc"));
+    }
+  virtual ~RealisticSimClusterMapper() {}
+  RealisticSimClusterMapper(const RealisticSimClusterMapper&) = delete;
+  RealisticSimClusterMapper& operator=(const RealisticSimClusterMapper&) = delete;
+
+  virtual void updateEvent(const edm::Event&) override final;
+  virtual void update(const edm::EventSetup&) override final;
+
+  void buildClusters(const edm::Handle<reco::PFRecHitCollection>&,
+		     const std::vector<bool>&,
+		     const std::vector<bool>&, 
+		     reco::PFClusterCollection&) override;
+  
+ private:  
+  hgcal::RecHitTools rhtools_;
+  const float invisibleFraction_ = 0.3f;
+  const float exclusiveFraction_ = 0.7f;
+  const bool maxDistanceFilter_ = false;
+  const float maxDistance_ = 10.f;
+  const bool useMCFractionsForExclEnergy_ = false;
+  edm::EDGetTokenT<SimClusterCollection> simClusterToken_;
+  edm::Handle<SimClusterCollection> simClusterH_;
+  
+};
+
+DEFINE_EDM_PLUGIN(InitialClusteringStepFactory,
+		  RealisticSimClusterMapper,
+		  "RealisticSimClusterMapper");
+
+#endif
diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGC_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGC_cfi.py
index 784e58e02775e..781af07c71167 100644
--- a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGC_cfi.py
+++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGC_cfi.py
@@ -14,12 +14,19 @@
 
 # initial step clusterizer
 _simClusterMapper_HGCal = cms.PSet(
-    algoName = cms.string("GenericSimClusterMapper"),
+    algoName = cms.string("RealisticSimClusterMapper"),
+    exclusiveFraction = cms.double(0.6),
+    invisibleFraction = cms.double(0.6),
+    maxDistanceFilter = cms.bool(True),
+    #filtering out hits outside a cylinder of 10cm radius, built around the center of gravity per each layer
+    maxDistance =  cms.double(10.0),
+    useMCFractionsForExclEnergy = cms.bool(False),    
     thresholdsByDetector = cms.VPSet(
     ),
     simClusterSrc = cms.InputTag("mix:MergedCaloTruth")
 )
 
+
 #position calculations
 _positionCalcPCA_HGCal = cms.PSet(
         algoName = cms.string("Cluster3DPCACalculator"),
diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGC_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGC_cfi.py
index d6f0dc2142f83..680f1213357dd 100644
--- a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGC_cfi.py
+++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGC_cfi.py
@@ -24,10 +24,10 @@
              qualityTests = cms.VPSet(
 # Enabling PFRecHitQTestHGCalThresholdSNR will filter out of the PFRecHits, all the HGCRecHits with energy not exceeding
 # 5 sigma noise                 
-#                 cms.PSet(
-#                     name = cms.string("PFRecHitQTestHGCalThresholdSNR"),
-#                     thresholdSNR = cms.double(5.0),
-#                    ),
+                cms.PSet(
+                    name = cms.string("PFRecHitQTestHGCalThresholdSNR"),
+                    thresholdSNR = cms.double(5.0),
+                   ),
                 )              
            ),
            cms.PSet(
@@ -37,10 +37,10 @@
              qualityTests = cms.VPSet(
 # Enabling PFRecHitQTestHGCalThresholdSNR will filter out of the PFRecHits, all the HGCRecHits with energy not exceeding
 # 5 sigma noise                     
-#                 cms.PSet(
-#                     name = cms.string("PFRecHitQTestHGCalThresholdSNR"),
-#                     thresholdSNR = cms.double(5.0),
-#                    ),                 
+                cms.PSet(
+                    name = cms.string("PFRecHitQTestHGCalThresholdSNR"),
+                    thresholdSNR = cms.double(5.0),
+                   ),                 
                 )
            ),
            cms.PSet(
@@ -50,10 +50,10 @@
              qualityTests = cms.VPSet(
 # Enabling PFRecHitQTestHGCalThresholdSNR will filter out of the PFRecHits, all the HGCRecHits with energy not exceeding
 # 5 sigma noise                     
-#                 cms.PSet(
-#                     name = cms.string("PFRecHitQTestHGCalThresholdSNR"),
-#                     thresholdSNR = cms.double(5.0),
-#                    ),                 
+                cms.PSet(
+                    name = cms.string("PFRecHitQTestHGCalThresholdSNR"),
+                    thresholdSNR = cms.double(5.0),
+                   ),                 
                 )
            )
     )