From ccedb061dba2533c939a6be2b4f0471fba49bf69 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Wed, 13 Apr 2016 16:33:39 +0200
Subject: [PATCH 01/16] speed up avoiding eta phi computations

---
 .../src/Basic2DGenericPFlowPositionCalc.cc    | 104 +++++++++++-------
 .../src/Basic2DGenericPFlowPositionCalc.h     |  17 ++-
 .../src/ECAL2DPositionCalcWithDepthCorr.cc    |   8 +-
 .../src/PFMultiDepthClusterizer.cc            |  50 ++++-----
 .../src/PFMultiDepthClusterizer.h             |   2 +-
 5 files changed, 112 insertions(+), 69 deletions(-)

diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
index 94e3fd2c436d9..ef168d790da7e 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
@@ -4,7 +4,9 @@
 #include "FWCore/Utilities/interface/isFinite.h"
 
 #include <cmath>
-#include <unordered_map>
+#include "CommonTools/Utils/interface/DynArray.h"
+#include<iterator>
+#include <boost/function_output_iterator.hpp>
 
 #include "vdt/vdtMath.h"
 
@@ -31,47 +33,63 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
   double cl_timeweight=0.0;
   double max_e = 0.0;  
   PFLayer::Layer max_e_layer = PFLayer::NONE;
-  reco::PFRecHitRef refseed;  
   // find the seed and max layer and also calculate time
   //Michalis : Even if we dont use timing in clustering here we should fill
   //the time information for the cluster. This should use the timing resolution(1/E)
   //so the weight should be fraction*E^2
   //calculate a simplistic depth now. The log weighted will be done
   //in different stage  
-  for( const reco::PFRecHitFraction& rhf : cluster.recHitFractions() ) {
-    const reco::PFRecHitRef& refhit = rhf.recHitRef();
-    if( refhit->detId() == cluster.seed() ) refseed = refhit;
-    const double rh_fraction = rhf.fraction();
-    const double rh_rawenergy = refhit->energy();
+
+
+  auto const recHitCollection = &(*cluster.recHitFractions()[0].recHitRef()) - cluster.recHitFractions()[0].recHitRef().key();
+  auto nhits = cluster.recHitFractions().size();
+  struct LHit{ reco::PFRecHit const * hit; double energy; double fraction;};
+  declareDynArray(LHit,nhits,hits);
+  for(auto i=0U; i<nhits; ++i) {
+    auto const & hf = cluster.recHitFractions()[i];
+    auto k = hf.recHitRef().key();
+    auto p = recHitCollection+k;
+    hits[i]= {p,(*p).energy(), hf.fraction()}; 
+  }
+
+  if(_posCalcNCrystals != -1) // sorted to make neighbour search faster
+    std::sort(hits.begin(),hits.end(),[](LHit const& a, LHit const& b) { return a.hit<b.hit;});
+
+  LHit mySeed={nullptr}; 
+  for( auto const & rhf : hits ) {
+    const reco::PFRecHit & refhit = *rhf.hit;
+    if( refhit.detId() == cluster.seed() ) mySeed = rhf;
+    const double rh_fraction = rhf.fraction;
+    const double rh_rawenergy = rhf.energy;
     const double rh_energy = rh_rawenergy * rh_fraction;   
     if( edm::isNotFinite(rh_energy) ) {
       throw cms::Exception("PFClusterAlgo")
-	<<"rechit " << refhit->detId() << " has a NaN energy... " 
+	<<"rechit " << refhit.detId() << " has a NaN energy... " 
 	<< "The input of the particle flow clustering seems to be corrupted.";
     }
     cl_energy += rh_energy;
     // If time resolution is given, calculated weighted average
-    if (_timeResolutionCalcBarrel && _timeResolutionCalcEndcap) {
-      double res2 = 10000.;
-      int cell_layer = (int)refhit->layer();
+    if ( bool(_timeResolutionCalcBarrel) & bool(_timeResolutionCalcEndcap) ) {
+      double res2 = 1.e-4;
+      int cell_layer = (int)refhit.layer();
       if (cell_layer == PFLayer::HCAL_BARREL1 ||
           cell_layer == PFLayer::HCAL_BARREL2 ||
           cell_layer == PFLayer::ECAL_BARREL)
-        res2 = _timeResolutionCalcBarrel->timeResolution2(rh_rawenergy);
+        res2 = 1./_timeResolutionCalcBarrel->timeResolution2(rh_rawenergy);
       else
-        res2 = _timeResolutionCalcEndcap->timeResolution2(rh_rawenergy);
-      cl_time += rh_fraction*refhit->time()/res2;
-      cl_timeweight += rh_fraction/res2;
+        res2 = 1./_timeResolutionCalcEndcap->timeResolution2(rh_rawenergy);
+      cl_time += rh_fraction*refhit.time()*res2;
+      cl_timeweight += rh_fraction*res2;
     }
     else { // assume resolution = 1/E**2
       const double rh_rawenergy2 = rh_rawenergy*rh_rawenergy;
       cl_timeweight+=rh_rawenergy2*rh_fraction;
-      cl_time += rh_rawenergy2*rh_fraction*refhit->time();
+      cl_time += rh_rawenergy2*rh_fraction*refhit.time();
     }
 
     if( rh_energy > max_e ) {
       max_e = rh_energy;
-      max_e_layer = rhf.recHitRef()->layer();
+      max_e_layer = refhit.layer();
     }    
   }
   cluster.setEnergy(cl_energy);
@@ -85,38 +103,50 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
   const reco::PFRecHitRefVector* seedNeighbours = nullptr;
   switch( _posCalcNCrystals ) {
   case 5:
-    seedNeighbours = &refseed->neighbours4();
+    seedNeighbours = &mySeed.hit->neighbours4();
     break;
   case 9:
-    seedNeighbours = &refseed->neighbours8();
-    break;
-  case -1:
+    seedNeighbours = &mySeed.hit->neighbours8();
     break;
   default:
-    assert(0); //bug
+    break;
   }
 
-  for( const reco::PFRecHitFraction& rhf : cluster.recHitFractions() ) {
-    const reco::PFRecHitRef& refhit = rhf.recHitRef();
-    
-    if( refhit != refseed && _posCalcNCrystals != -1 ) {
-      auto pos = std::find(seedNeighbours->begin(),seedNeighbours->end(),
-			   refhit);
-      if( pos == seedNeighbours->end() ) continue;
-    }
-    
-    const double rh_energy = refhit->energy() * ((float)rhf.fraction());
-    const double norm = ( rhf.fraction() < _minFractionInCalc ? 
+  auto compute = [&](LHit const& rhf) {
+    const reco::PFRecHit & refhit = *rhf.hit;  
+    const double rh_energy = rhf.energy * rhf.fraction;
+    const double norm = ( rhf.fraction < _minFractionInCalc ? 
 			  0.0 : 
-			  std::max(0.0,vdt::fast_log(rh_energy/_logWeightDenom)) );
-    const math::XYZPoint& rhpos_xyz = refhit->position();
+			  std::max(0.0,vdt::fast_log(rh_energy*_logWeightDenom)) );
+    const math::XYZPoint& rhpos_xyz = refhit.position();
     x += rhpos_xyz.X() * norm;
     y += rhpos_xyz.Y() * norm;
     z += rhpos_xyz.Z() * norm;
-    depth += refhit->depth()*norm;
-    
+    depth += refhit.depth()*norm;
     position_norm += norm;
+  };
+
+  if(_posCalcNCrystals == -1)
+    for( auto const & rhf : hits ) compute(rhf);
+  else {  // only seed and its neighbours
+     compute(mySeed);
+     // search seedNeighbours to find energy fraction in cluster (sic)
+     unInitDynArray(reco::PFRecHit const *,seedNeighbours->size(),nei);	  
+     for(auto k :seedNeighbours->refVector().keys()){ 
+      nei.push_back(&recHitCollection[k]);
+     }
+     std::sort(nei.begin(),nei.end());
+     struct LHitLess {
+       auto operator()(LHit const &a, reco::PFRecHit const * b) const {return a.hit<b;}
+       auto operator()(reco::PFRecHit const * b, LHit const &a) const {return b<a.hit;}
+     };
+     std::set_intersection(hits.begin(),hits.end(),nei.begin(),nei.end(), 
+       boost::make_function_output_iterator(compute),
+       LHitLess()
+     );
   }
+
+
   if( position_norm < _minAllowedNorm ) {
     edm::LogError("WeirdClusterNormalization") 
       << "PFCluster too far from seeding cell: set position to (0,0,0).";
diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.h b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.h
index b59d06796b5ae..d7d55c633da23 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.h
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.h
@@ -5,6 +5,8 @@
 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
 
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
 #include "RecoParticleFlow/PFClusterProducer/interface/CaloRecHitResolutionProvider.h"
 
 class Basic2DGenericPFlowPositionCalc : public PFCPositionCalculatorBase {
@@ -12,7 +14,7 @@ class Basic2DGenericPFlowPositionCalc : public PFCPositionCalculatorBase {
   Basic2DGenericPFlowPositionCalc(const edm::ParameterSet& conf) :
     PFCPositionCalculatorBase(conf),    
     _posCalcNCrystals(conf.getParameter<int>("posCalcNCrystals")),
-    _logWeightDenom(conf.getParameter<double>("logWeightDenominator")),
+    _logWeightDenom(1./conf.getParameter<double>("logWeightDenominator")),
     _minAllowedNorm(conf.getParameter<double>("minAllowedNormalization"))
 
   {  
@@ -28,7 +30,20 @@ class Basic2DGenericPFlowPositionCalc : public PFCPositionCalculatorBase {
         conf.getParameterSet("timeResolutionCalcEndcap");
         _timeResolutionCalcEndcap.reset(new CaloRecHitResolutionProvider(timeResConf));
     }
+
+   switch( _posCalcNCrystals ) {
+    case 5:
+    case 9:
+    case -1:
+      break;
+    default:
+      edm::LogError("Basic2DGenericPFlowPositionCalc") << "posCalcNCrystals not valid";
+      assert(0); // bug
+   }
+
+
   }
+
   Basic2DGenericPFlowPositionCalc(const Basic2DGenericPFlowPositionCalc&) = delete;
   Basic2DGenericPFlowPositionCalc& operator=(const Basic2DGenericPFlowPositionCalc&) = delete;
 
diff --git a/RecoParticleFlow/PFClusterProducer/src/ECAL2DPositionCalcWithDepthCorr.cc b/RecoParticleFlow/PFClusterProducer/src/ECAL2DPositionCalcWithDepthCorr.cc
index 979913f488515..50fee4becd80c 100644
--- a/RecoParticleFlow/PFClusterProducer/src/ECAL2DPositionCalcWithDepthCorr.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/ECAL2DPositionCalcWithDepthCorr.cc
@@ -79,9 +79,9 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
     cl_energy_float += rh_energyf;
     // If time resolution is given, calculate weighted average
     if (_timeResolutionCalc) {
-      const double res2 = _timeResolutionCalc->timeResolution2(rh_rawenergy);
-      cl_time += rh_fraction*refhit->time()/res2;
-      cl_timeweight += rh_fraction/res2;
+      const double res2 = 1./_timeResolutionCalc->timeResolution2(rh_rawenergy);
+      cl_time += rh_fraction*refhit->time()*res2;
+      cl_timeweight += rh_fraction*res2;
     }
     else { // assume resolution ~ 1/E**2
       const double rh_rawenergy2 = rh_rawenergy*rh_rawenergy;
@@ -115,7 +115,7 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
 
   const CaloCellGeometry* center_cell = 
     ecal_geom->getGeometry(refmax->detId());
-  const double ctreta = center_cell->getPosition().eta();
+  const double ctreta = center_cell->etaPos();
   const double actreta = std::abs(ctreta);
   // need to change T0 if in ES
   if( actreta > preshowerStartEta && actreta < preshowerEndEta ) { 
diff --git a/RecoParticleFlow/PFClusterProducer/src/PFMultiDepthClusterizer.cc b/RecoParticleFlow/PFClusterProducer/src/PFMultiDepthClusterizer.cc
index aec63e5145b48..01a2d54860dc6 100644
--- a/RecoParticleFlow/PFClusterProducer/src/PFMultiDepthClusterizer.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/PFMultiDepthClusterizer.cc
@@ -37,8 +37,8 @@ buildClusters(const reco::PFClusterCollection& input,
 	      const std::vector<bool>& seedable,
 	      reco::PFClusterCollection& output) {
 
-  std::vector<double> etaRMS(input.size(),0.0);
-  std::vector<double> phiRMS(input.size(),0.0);
+  std::vector<double> etaRMS2(input.size(),0.0);
+  std::vector<double> phiRMS2(input.size(),0.0);
 
   //We need to sort the clusters for smaller to larger depth
   //  for (unsigned int i=0;i<input.size();++i)
@@ -47,10 +47,10 @@ buildClusters(const reco::PFClusterCollection& input,
 
 
   //calculate cluster shapes
-  calculateShowerShapes(input,etaRMS,phiRMS);
+  calculateShowerShapes(input,etaRMS2,phiRMS2);
 
   //link
-  std::vector<ClusterLink> links = link(input,etaRMS,phiRMS); 
+  auto && links = link(input,etaRMS2,phiRMS2); 
   //  for (const auto& link: links)
   //    printf("link %d %d %f %f\n",link.from(),link.to(),link.dR(),link.dZ());
 
@@ -60,9 +60,7 @@ buildClusters(const reco::PFClusterCollection& input,
   std::vector<bool> linked(input.size(),false);
 
   //prune
-  std::vector<ClusterLink> prunedLinks;
-  if (links.size())
-    prunedLinks =  prune(links,linked);
+  auto && prunedLinks =  prune(links,linked);
 
   //printf("Pruned links\n")
   //  for (const auto& link: prunedLinks)
@@ -97,15 +95,13 @@ buildClusters(const reco::PFClusterCollection& input,
 
 
 void PFMultiDepthClusterizer::
-calculateShowerShapes(const reco::PFClusterCollection& clusters,std::vector<double>& etaRMS,std::vector<double>& phiRMS) {
-  float etaSum;
-  float phiSum;
+calculateShowerShapes(const reco::PFClusterCollection& clusters,std::vector<double>& etaRMS2,std::vector<double>& phiRMS2) {
 
   //shower shapes. here do not use the fractions 
 
   for( unsigned int i=0;i<clusters.size();++i ) {
     const reco::PFCluster& cluster = clusters[i]; 
-    etaSum=0.0;phiSum=0.0;
+    double etaSum=0.0; double phiSum=0.0;
     auto const & crep = cluster.positionREP();
     for (const auto& frac : cluster.recHitFractions()) {
       auto const & h = *frac.recHitRef();
@@ -114,42 +110,44 @@ calculateShowerShapes(const reco::PFClusterCollection& clusters,std::vector<doub
       phiSum += (frac.fraction()*h.energy())*std::abs(deltaPhi(rep.phi(),crep.phi()));
     }
     //protection for single line : assign ~ tower
-    etaRMS[i] = std::max(etaSum/cluster.energy(),0.1);
-    phiRMS[i] = std::max(phiSum/cluster.energy(),0.1);
-    
+    etaRMS2[i] = std::max(etaSum/cluster.energy(),0.1);
+    etaRMS2[i]*= etaRMS2[i];
+    phiRMS2[i] = std::max(phiSum/cluster.energy(),0.1);
+    phiRMS2[i]*=phiRMS2[i];
   }
 
 }
 
 
 std::vector<PFMultiDepthClusterizer::ClusterLink>   PFMultiDepthClusterizer::
-link(const reco::PFClusterCollection& clusters ,const std::vector<double>& etaRMS,const std::vector<double>& phiRMS) {
+link(const reco::PFClusterCollection& clusters ,const std::vector<double>& etaRMS2,const std::vector<double>& phiRMS2) {
 
   std::vector<ClusterLink> links;
   //loop on all pairs
   for (unsigned int i=0;i<clusters.size();++i)
     for (unsigned int j=0;j<clusters.size();++j) {
-      if (i==j )
-	continue;
+      if (i==j )  continue;
 
       const reco::PFCluster& cluster1 = clusters[i]; 
       const reco::PFCluster& cluster2 = clusters[j]; 
 
-      float dz = (cluster2.depth() - cluster1.depth());
+      auto dz = (cluster2.depth() - cluster1.depth());
 
       //Do not link at the same layer and only link inside out!
-      if (dz<0.0 || fabs(dz)<0.2)
-	continue;
+      if (dz<0.0f || std::abs(dz)<0.2f)  continue;
+ 
       auto const & crep1 = cluster1.positionREP();
       auto const & crep2 = cluster2.positionREP();
-
-      float deta =(crep1.eta()-crep2.eta())*(crep1.eta()-crep2.eta())/(etaRMS[i]*etaRMS[i]+etaRMS[j]*etaRMS[j]);
-      float dphi = deltaPhi(crep1.phi(),crep2.phi())*deltaPhi(crep1.phi(),crep2.phi())/(phiRMS[i]*phiRMS[i]+phiRMS[j]*phiRMS[j]);
+ 
+      auto deta = crep1.eta()-crep2.eta();
+      deta = deta*deta/(etaRMS2[i]+etaRMS2[j]);
+      auto dphi = deltaPhi(crep1.phi(),crep2.phi());
+      dphi = dphi*dphi/(phiRMS2[i]+phiRMS2[j]);
 
       //      printf("Testing Link %d -> %d (%f %f %f %f ) \n",i,j,deta,dphi,cluster1.position().Eta()-cluster2.position().Eta(),deltaPhi(cluster1.position().Phi(),cluster2.position().Phi()));
 
-      if (deta<nSigmaEta_ && dphi<nSigmaPhi_ ) 
-	links.push_back(ClusterLink(i,j,deta+dphi,fabs(dz),cluster1.energy()+cluster2.energy()));
+      if ( (deta<nSigmaEta_) & (dphi<nSigmaPhi_ ) )
+	links.push_back(ClusterLink(i,j,deta+dphi,std::abs(dz),cluster1.energy()+cluster2.energy()));
     }
 
       return links;
@@ -159,12 +157,12 @@ std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::
 prune(std::vector<ClusterLink>& links,std::vector<bool>& linkedClusters) {
       std::vector<ClusterLink> goodLinks ;
       std::vector<bool> mask(links.size(),false);
+      if (links.empty()) return goodLinks;
 
       for (unsigned int i=0;i<links.size()-1;++i) {
 	if (mask[i])
 	  continue;
 	for (unsigned int j=i+1;j<links.size();++j) {
-
 	  if (mask[j])
 	    continue;
 	  
diff --git a/RecoParticleFlow/PFClusterProducer/src/PFMultiDepthClusterizer.h b/RecoParticleFlow/PFClusterProducer/src/PFMultiDepthClusterizer.h
index 4f0aa7bf3f71d..eb3104acb1041 100644
--- a/RecoParticleFlow/PFClusterProducer/src/PFMultiDepthClusterizer.h
+++ b/RecoParticleFlow/PFClusterProducer/src/PFMultiDepthClusterizer.h
@@ -6,7 +6,7 @@
 
 #include <unordered_map>
 
-class PFMultiDepthClusterizer : public PFClusterBuilderBase {
+class PFMultiDepthClusterizer final : public PFClusterBuilderBase {
   typedef PFMultiDepthClusterizer B2DGPF;
  public:
   PFMultiDepthClusterizer(const edm::ParameterSet& conf);

From 90ef330c19dc495a882d3e07bc18305021f0b4a9 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Fri, 15 Apr 2016 18:35:37 +0200
Subject: [PATCH 02/16] add rhoetaphi to calocell (for corners as well)

---
 DataFormats/Math/interface/PtEtaPhiMass.h     | 19 ++++++++++++++++
 .../ParticleFlowReco/interface/PFRecHit.h     |  8 ++-----
 .../CaloGeometry/interface/CaloCellGeometry.h | 22 ++++++++++++++-----
 Geometry/CaloGeometry/src/CaloCellGeometry.cc |  6 ++---
 4 files changed, 40 insertions(+), 15 deletions(-)

diff --git a/DataFormats/Math/interface/PtEtaPhiMass.h b/DataFormats/Math/interface/PtEtaPhiMass.h
index 44fa6deb59022..30134da5f2bc7 100644
--- a/DataFormats/Math/interface/PtEtaPhiMass.h
+++ b/DataFormats/Math/interface/PtEtaPhiMass.h
@@ -34,5 +34,24 @@ class PtEtaPhiMass {
 };
 
 
+class RhoEtaPhi {
+private:
+  float rho_, eta_, phi_;
+public:
+  // default constructor (unitialized)
+  RhoEtaPhi() {}
+
+  //positional constructor (still compatible with Root, c++03)
+  RhoEtaPhi(float irho, float ieta, float iphi):
+    rho_(irho), eta_(ieta), phi_(iphi) {}
+
+  /// transverse momentum
+  float rho() const { return rho_;}
+  /// momentum pseudorapidity
+  float eta() const { return eta_; }
+  /// momentum azimuthal angle
+  float phi() const { return phi_; }
+};
+
 
 #endif
diff --git a/DataFormats/ParticleFlowReco/interface/PFRecHit.h b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
index 51ca21e9db8c6..eddf25565c8c7 100644
--- a/DataFormats/ParticleFlowReco/interface/PFRecHit.h
+++ b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
@@ -7,20 +7,16 @@
 #include <iostream>
 
 #include "DataFormats/Math/interface/Point3D.h"
-#include "Rtypes.h" 
 #include "DataFormats/Math/interface/Vector3D.h"
 #include "Math/GenVector/PositionVector3D.h"
 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
 
-//C decide what is the default rechit index. 
-//C maybe 0 ? -> compression 
-//C then the position is index-1. 
-//C provide a helper class to access the rechit. 
-
 #include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
 #include "DataFormats/Common/interface/RefToBase.h"
 
+#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
+
 
 namespace reco {
 
diff --git a/Geometry/CaloGeometry/interface/CaloCellGeometry.h b/Geometry/CaloGeometry/interface/CaloCellGeometry.h
index a660e39655e35..0dafeb092d2fd 100644
--- a/Geometry/CaloGeometry/interface/CaloCellGeometry.h
+++ b/Geometry/CaloGeometry/interface/CaloCellGeometry.h
@@ -6,11 +6,13 @@
 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
 #include <CLHEP/Geometry/Point3D.h>
 #include <CLHEP/Geometry/Transform3D.h>
+#include "DataFormats/Math/interface/PtEtaPhiMass.h"
+
 #include <vector>
+#include <array>
 #include <string>
 #include <cassert>
 
-#include "FWCore/Utilities/interface/GCC11Compatibility.h"
 
 
 
@@ -77,8 +79,10 @@ class CaloCellGeometry
   const GlobalPoint& getPosition() const {return m_refPoint;}
   const GlobalPoint& getBackPoint() const {return m_backPoint;} 
 
-  float etaPos() const { return m_eta;}
-  float phiPos() const { return m_phi;}
+  RhoEtaPhi const & repPos() const { return m_rep;}
+  float rhoPos() const { return m_rep.rho();}
+  float etaPos() const { return m_rep.eta();}
+  float phiPos() const { return m_rep.phi();}
 
   float etaSpan() const { return m_dEta;}
   float	phiSpan() const	{ return m_dPhi;}
@@ -126,6 +130,7 @@ class CaloCellGeometry
      m_dPhi = std::abs(getCorners()[0].phi() -
                       getCorners()[2].phi());
      initBack();
+     initReps();
   }
 
   virtual void initCorners(CornersVec&) = 0;
@@ -137,18 +142,23 @@ class CaloCellGeometry
                               0.25 * (cv[4].y() + cv[5].y() + cv[6].y() + cv[7].y()),
                               0.25 * (cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z()));   
   }
-
+  void initReps() {
+    for (int i=0;i<k_cornerSize; ++i) m_repCorners[i]= {getCorners()[i].perp(), getCorners()[i].eta(), getCorners()[i].barePhi()}; 
+      
+  }
+  
 
   GlobalPoint         m_refPoint ;
   GlobalPoint         m_backPoint ;
   CornersVec  m_corners  ;
   const CCGFloat*     m_parms    ;
-  float m_eta, m_phi;
+  RhoEtaPhi m_rep;
   float m_dEta;
   float m_dPhi;
-
+  std::array<RhoEtaPhi,k_cornerSize> m_repCorners;
 };
 
 std::ostream& operator<<( std::ostream& s, const CaloCellGeometry& cell ) ;
 
+
 #endif
diff --git a/Geometry/CaloGeometry/src/CaloCellGeometry.cc b/Geometry/CaloGeometry/src/CaloCellGeometry.cc
index e9898a00c0e27..e094bc56fceae 100644
--- a/Geometry/CaloGeometry/src/CaloCellGeometry.cc
+++ b/Geometry/CaloGeometry/src/CaloCellGeometry.cc
@@ -19,7 +19,7 @@ const float CaloCellGeometry::k_ScaleFromDDDtoGeant ( 0.1 ) ;
 CaloCellGeometry::CaloCellGeometry() :
    m_refPoint ( 0., 0., 0. ),
    m_corners  (  ) ,
-   m_parms    ( (CCGFloat*) 0 ), m_eta(0), m_phi(0) 
+   m_parms    ( (CCGFloat*) 0 )
 {}
 
 
@@ -34,7 +34,7 @@ CaloCellGeometry::CaloCellGeometry( CornersVec::const_reference gp ,
    m_refPoint ( gp  ),
    m_corners  ( mgr ),
    m_parms    ( par ),
-   m_eta(gp.eta()), m_phi(gp.phi())
+   m_rep(gp.perp(),gp.eta(),gp.barePhi())
 {}
 
 CaloCellGeometry::CaloCellGeometry( const CornersVec& cv,
@@ -44,7 +44,7 @@ CaloCellGeometry::CaloCellGeometry( const CornersVec& cv,
 		0.25*( cv[0].z() + cv[1].z() + cv[2].z() + cv[3].z() )  ), 
    m_corners  ( cv ),
    m_parms    ( par ),
-   m_eta( m_refPoint.eta()), m_phi(m_refPoint.phi())
+   m_rep( m_refPoint.perp(), m_refPoint.eta(), m_refPoint.barePhi())
 {}
 
 std::ostream& operator<<( std::ostream& s, const CaloCellGeometry& cell ) 

From e82219b711661c2a4f74d9fcefb095fe28cc002f Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sat, 16 Apr 2016 10:55:02 +0200
Subject: [PATCH 03/16] first version of Sli PFRecHit

---
 .../ParticleFlowReco/interface/PFRecHit.h     | 115 +++--------
 DataFormats/ParticleFlowReco/src/PFCluster.cc |  42 ++++
 DataFormats/ParticleFlowReco/src/PFRecHit.cc  | 184 +-----------------
 .../ParticleFlowReco/src/classes_def_2.xml    |   9 +-
 .../CaloGeometry/interface/CaloCellGeometry.h |  11 +-
 5 files changed, 88 insertions(+), 273 deletions(-)

diff --git a/DataFormats/ParticleFlowReco/interface/PFRecHit.h b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
index eddf25565c8c7..e88503e824b63 100644
--- a/DataFormats/ParticleFlowReco/interface/PFRecHit.h
+++ b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
@@ -32,34 +32,25 @@ namespace reco {
 
   public:
 
-    // Next typedef uses double in ROOT 6 rather than Double32_t due to a bug in ROOT 5,
-    // which otherwise would make ROOT5 files unreadable in ROOT6.  This does not increase
-    // the size on disk, because due to the bug, double was actually stored on disk in ROOT 5.
-
-    typedef ROOT::Math::PositionVector3D<ROOT::Math::CylindricalEta3D<double> > REPPoint;
-
-    typedef std::vector<REPPoint> REPPointVector;
-
+    using REPPoint = RhoEtaPhi;
+    using RepCorners = CaloCellGeometry::RepCorners;
+    using REPPointVector = RepCorners;
+    using CornersVec = CaloCellGeometry::CornersVec;
+    
     enum {
       NONE=0
     };
     /// default constructor. Sets energy and position to zero
-    PFRecHit();
+    PFRecHit(){}
 
-    /// constructor from values
-    PFRecHit(unsigned detId,
+    PFRecHit(CaloCellGeometry const * caloCell, unsigned detId,
              PFLayer::Layer layer,
-             double energy, 
-             const math::XYZPoint& posxyz, 
-             const math::XYZVector& axisxyz, 
-             const std::vector< math::XYZPoint >& cornersxyz);
+             double energy) :
+        caloCell_(caloCell),  detId_(detId),
+        layer_(layer), energy_(energy){}
 
-    PFRecHit(unsigned detId,
-             PFLayer::Layer layer,
-             double energy, 
-             double posx, double posy, double posz, 
-             double axisx, double axisy, double axisz);    
 
+    
     /// copy
     PFRecHit(const PFRecHit& other) = default;
     PFRecHit(PFRecHit&& other) = default;
@@ -72,7 +63,6 @@ namespace reco {
 
     void setEnergy( double energy) { energy_ = energy; }
 
-    void calculatePositionREP();
 
     void addNeighbour(short x,short y, short z,const PFRecHitRef&);
     const PFRecHitRef getNeighbour(short x,short y, short z);
@@ -98,11 +88,10 @@ namespace reco {
     }
 
 
-    void      setNWCorner( double posx, double posy, double posz );
-    void      setSWCorner( double posx, double posy, double posz );
-    void      setSECorner( double posx, double posy, double posz );
-    void      setNECorner( double posx, double posy, double posz );
-
+    /// calo cell
+    CaloCellGeometry const & callCell() const { return  *caloCell_; }
+    bool hasCaloCel() const { return caloCell_; }
+    
     /// rechit detId
     unsigned detId() const {return detId_;}
 
@@ -121,30 +110,20 @@ namespace reco {
 
     /// rechit momentum transverse to the beam, squared.
     double pt2() const { return energy_ * energy_ *
-			   ( position_.X()*position_.X() + 
-			     position_.Y()*position_.Y() ) / 
-			   ( position_.X()*position_.X() +
-			     position_.Y()*position_.Y() + 
-			     position_.Z()*position_.Z()) ; }
+	( position().basicVector().perp2()/ position().basicVector().mag2());}
 
 
     /// rechit cell centre x, y, z
-    const math::XYZPoint& position() const { return position_; }
-
-    const REPPoint& positionREP() const { return positionrep_; }
-
-
-    /// rechit cell axis x, y, z
-    const math::XYZVector& getAxisXYZ() const { return axisxyz_; }    
+    GlobalPoint const & position() const { return callCell().getPosition(); }
+    
+    RhoEtaPhi const &  positionREP() const { return callCell().repPos(); }
 
     /// rechit corners
-    const std::vector< math::XYZPoint >& getCornersXYZ() const 
-      { return cornersxyz_; }    
-
-    const std::vector<REPPoint>& getCornersREP() const { return cornersrep_; }
-
-    void size(double& deta, double& dphi) const;
+    CornersVec const & getCornersXYZ() const { return callCell().getCorners(); }    
 
+    RepCorners const & getCornersREP() const { return callCell().getCornersREP();}
+ 
+ 
     /// comparison >= operator
     bool operator>=(const PFRecHit& rhs) const { return (energy_>=rhs.energy_); }
 
@@ -160,50 +139,25 @@ namespace reco {
     friend std::ostream& operator<<(std::ostream& out, 
                                     const reco::PFRecHit& hit);
 
-    const edm::RefToBase<CaloRecHit>& originalRecHit() const {
-      return originalRecHit_;
-    }
-
-    template<typename T> 
-    void setOriginalRecHit(const T& rh) {
-      originalRecHit_ = edm::RefToBase<CaloRecHit>(rh);
-    }
-
   private:
 
-    // original rechit
-    edm::RefToBase<CaloRecHit> originalRecHit_;
-
+    CaloCellGeometry const * caloCell_=nullptr;
+ 
     ///C cell detid - should be detid or index in collection ?
-    unsigned            detId_;             
+    unsigned  int        detId_=0;             
 
     /// rechit layer
-    PFLayer::Layer      layer_;
+    PFLayer::Layer      layer_=PFLayer::NONE;
 
     /// rechit energy 
-    double              energy_;
+    double              energy_=0;
 
     /// time
-    double              time_;
-
+    double              time_=-1;
 
     /// depth
-    int      depth_;
-
-    /// rechit cell centre: x, y, z
-    math::XYZPoint      position_;
-
-    /// rechit cell centre: rho, eta, phi (transient)
-    REPPoint positionrep_;
-
-    /// rechit cell axisxyz
-    math::XYZVector     axisxyz_;
+    int      depth_=0;
 
-    /// rechit cell corners
-    std::vector< math::XYZPoint > cornersxyz_;
-
-    /// rechit cell corners rho/eta/phi
-    std::vector< REPPoint > cornersrep_;
   
     /// indices to existing neighbours (1 common side)
     PFRecHitRefVector   neighbours_;
@@ -212,15 +166,8 @@ namespace reco {
     //Caching the neighbours4/8 per request of Lindsey
     PFRecHitRefVector   neighbours4_;
     PFRecHitRefVector   neighbours8_;
+  };
 
 
-    /// number of corners
-    static const unsigned    nCorners_;
-
-    /// set position of one of the corners
-    void      setCorner( unsigned i, double posx, double posy, double posz );
-  };
-  
 }
-
 #endif
diff --git a/DataFormats/ParticleFlowReco/src/PFCluster.cc b/DataFormats/ParticleFlowReco/src/PFCluster.cc
index 9f4b6f9b30eff..8cf51f955bdf0 100644
--- a/DataFormats/ParticleFlowReco/src/PFCluster.cc
+++ b/DataFormats/ParticleFlowReco/src/PFCluster.cc
@@ -1,5 +1,47 @@
 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
 
+#include "vdt/vdtMath.h"
+#include "Math/GenVector/etaMax.h"
+
+namespace {
+  
+  // an implementation of Eta_FromRhoZ of root libraries using vdt
+  template<typename Scalar>
+  inline Scalar Eta_FromRhoZ_fast(Scalar rho, Scalar z) {    
+    using namespace ROOT::Math;
+    // value to control Taylor expansion of sqrt
+    const Scalar big_z_scaled =
+      std::pow(std::numeric_limits<Scalar>::epsilon(),static_cast<Scalar>(-.25));
+    if (rho > 0) {      
+      Scalar z_scaled = z/rho;
+      if (std::fabs(z_scaled) < big_z_scaled) {
+        return vdt::fast_log(z_scaled+std::sqrt(z_scaled*z_scaled+1.0));
+      } else {
+        // apply correction using first order Taylor expansion of sqrt
+        return  z>0 ? vdt::fast_log(2.0*z_scaled + 0.5/z_scaled) : -vdt::fast_log(-2.0*z_scaled);
+      }
+    }
+    // case vector has rho = 0
+    else if (z==0) {
+      return 0;
+    }
+    else if (z>0) {
+      return z + etaMax<Scalar>();
+    }
+    else {
+      return z - etaMax<Scalar>();
+    }    
+  }
+
+  inline void calculateREP(const math::XYZPoint& pos, double& rho, double& eta, double& phi) {
+    const double z = pos.z();
+    rho = pos.Rho();    
+    eta = Eta_FromRhoZ_fast<double>(rho,z);
+    phi = (pos.x()==0 && pos.y()==0) ? 0 : vdt::fast_atan2(pos.y(), pos.x());
+  }
+
+}
+
 using namespace std;
 using namespace reco;
 
diff --git a/DataFormats/ParticleFlowReco/src/PFRecHit.cc b/DataFormats/ParticleFlowReco/src/PFRecHit.cc
index d458d6e309f1f..11ad16a3f367c 100644
--- a/DataFormats/ParticleFlowReco/src/PFRecHit.cc
+++ b/DataFormats/ParticleFlowReco/src/PFRecHit.cc
@@ -1,162 +1,8 @@
 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
+
 using namespace reco;
 using namespace std;
 
-#include "vdt/vdtMath.h"
-#include "Math/GenVector/etaMax.h"
-
-const unsigned    PFRecHit::nCorners_ = 4;
-
-namespace {
-  
-  // an implementation of Eta_FromRhoZ of root libraries using vdt
-  template<typename Scalar>
-  inline Scalar Eta_FromRhoZ_fast(Scalar rho, Scalar z) {    
-    using namespace ROOT::Math;
-    // value to control Taylor expansion of sqrt
-    const Scalar big_z_scaled =
-      std::pow(std::numeric_limits<Scalar>::epsilon(),static_cast<Scalar>(-.25));
-    if (rho > 0) {      
-      Scalar z_scaled = z/rho;
-      if (std::fabs(z_scaled) < big_z_scaled) {
-        return vdt::fast_log(z_scaled+std::sqrt(z_scaled*z_scaled+1.0));
-      } else {
-        // apply correction using first order Taylor expansion of sqrt
-        return  z>0 ? vdt::fast_log(2.0*z_scaled + 0.5/z_scaled) : -vdt::fast_log(-2.0*z_scaled);
-      }
-    }
-    // case vector has rho = 0
-    else if (z==0) {
-      return 0;
-    }
-    else if (z>0) {
-      return z + etaMax<Scalar>();
-    }
-    else {
-      return z - etaMax<Scalar>();
-    }    
-  }
-
-  inline void calculateREP(const math::XYZPoint& pos, double& rho, double& eta, double& phi) {
-    const double z = pos.z();
-    rho = pos.Rho();    
-    eta = Eta_FromRhoZ_fast<double>(rho,z);
-    phi = (pos.x()==0 && pos.y()==0) ? 0 : vdt::fast_atan2(pos.y(), pos.x());
-  }
-
-}
-
-PFRecHit::PFRecHit() : 
-  detId_(0),
-  layer_(PFLayer::NONE),
-  energy_(0.), 
-  time_(-1.),
-  depth_(0),
-  position_(math::XYZPoint(0.,0.,0.))
-{
-  
-  cornersxyz_.reserve( nCorners_ );
-  for(unsigned i=0; i<nCorners_; i++) { 
-    cornersxyz_.push_back( position_ );    
-  }
-  calculatePositionREP();
-}
-
-
-PFRecHit::PFRecHit(unsigned detId,
-                   PFLayer::Layer layer, 
-                   double energy, 
-                   const math::XYZPoint& position,
-                   const math::XYZVector& axisxyz,
-                   const vector< math::XYZPoint >& cornersxyz) : 
-  detId_(detId),
-  layer_(layer),
-  energy_(energy), 
-  time_(-1.),
-  depth_(0),
-  position_(position),
-  axisxyz_(axisxyz),
-  cornersxyz_(cornersxyz) 
-{
-  calculatePositionREP();
-}
-
-PFRecHit::PFRecHit(unsigned detId,
-                   PFLayer::Layer layer,
-                   double energy, 
-                   double posx, double posy, double posz, 
-                   double axisx, double axisy, double axisz) :
-
-  detId_(detId),
-  layer_(layer),
-  energy_(energy), 
-  time_(-1.),
-  depth_(0),
-  position_(posx, posy, posz),
-  axisxyz_(axisx, axisy, axisz) {
-
-  cornersxyz_.reserve( nCorners_ );
-  for(unsigned i=0; i<nCorners_; i++) { 
-    cornersxyz_.push_back( position_ );    
-  } 
-
-  calculatePositionREP();
-
-}    
-
-
-
-
-void PFRecHit::setNWCorner( double posx, double posy, double posz ) {
-  setCorner(0, posx, posy, posz);
-}
-
-
-void PFRecHit::setSWCorner( double posx, double posy, double posz ) {
-  setCorner(1, posx, posy, posz);
-}
-
-
-void PFRecHit::setSECorner( double posx, double posy, double posz ) {
-  setCorner(2, posx, posy, posz);
-}
-
-
-void PFRecHit::setNECorner( double posx, double posy, double posz ) {
-  setCorner(3, posx, posy, posz);
-}
-
-
-void PFRecHit::setCorner( unsigned i, double posx, double posy, double posz ) {
-  assert( cornersxyz_.size() == nCorners_);
-  assert( i<cornersxyz_.size() );
-  double rho(0), eta(0), phi(0);
-
-  cornersxyz_[i] = math::XYZPoint( posx, posy, posz);
-  const auto& corner = cornersxyz_[i];
-  calculateREP(corner, rho, eta, phi);  
-  cornersrep_[i] = REPPoint( rho,
-			     eta,
-			     phi );
-}
-
-void
-PFRecHit::calculatePositionREP() {
-  double rho(0), eta(0), phi(0);
-  calculateREP(position_,rho,eta,phi);
-
-  positionrep_.SetCoordinates( rho, 
-			       eta, 
-			       phi );
-  cornersrep_.resize(cornersxyz_.size());
-  for( unsigned i = 0; i < cornersxyz_.size(); ++i ) {
-    const auto& corner = cornersxyz_[i];
-    calculateREP(corner,rho,eta,phi);
-    cornersrep_[i].SetCoordinates( rho,
-                                   eta,
-                                   phi );
-  }
-}
 
 void PFRecHit::addNeighbour(short x,short y,short z,const PFRecHitRef& ref) {
   //bitmask interface  to accomodate more advanced naighbour finding [i.e in z as well]
@@ -221,40 +67,16 @@ const PFRecHitRef PFRecHit::getNeighbour(short x,short y,short z) {
 }
 
 
-void PFRecHit::size(double& deta, double& dphi) const {
-
-  double minphi=9999;
-  double maxphi=-9999;
-  double mineta=9999;
-  double maxeta=-9999;
-  for ( unsigned ic=0; ic<cornersxyz_.size(); ++ic ) { 
-    double eta = cornersxyz_[ic].Eta();
-    double phi = cornersxyz_[ic].Phi();
-    
-    if(phi>maxphi) maxphi=phi;
-    if(phi<minphi) minphi=phi;
-    if(eta>maxeta) maxeta=eta;
-    if(eta<mineta) mineta=eta;    
-  }
-
-  deta = maxeta - mineta;
-  dphi = maxphi - minphi;
-}
-
-
 ostream& reco::operator<<(ostream& out, const reco::PFRecHit& hit) {
 
   if(!out) return out;
 
-  //   reco::PFRecHit& nshit = const_cast<reco::PFRecHit& >(hit);
-  //   const reco::PFRecHit::REPPoint& posrep = nshit.positionREP();
-  
-  const  math::XYZPoint& posxyz = hit.position();
+  auto const & pos = hit.positionREP();
 
   out<<"hit id:"<<hit.detId()
      <<" l:"<<hit.layer()
      <<" E:"<<hit.energy()
      <<" t:"<<hit.time()
-     <<" rep:"<<posxyz.Rho()<<","<<posxyz.Eta()<<","<<posxyz.Phi()<<"| N:";
+     <<" rep:"<<pos.rho()<<","<<pos.eta()<<","<<pos.phi()<<"| N:";
   return out;
 }
diff --git a/DataFormats/ParticleFlowReco/src/classes_def_2.xml b/DataFormats/ParticleFlowReco/src/classes_def_2.xml
index 9e93f60b5917a..fd03a39045a4d 100644
--- a/DataFormats/ParticleFlowReco/src/classes_def_2.xml
+++ b/DataFormats/ParticleFlowReco/src/classes_def_2.xml
@@ -36,7 +36,8 @@
   </class>
 
 
-  <class name="reco::PFRecHit" ClassVersion="17">
+  <class name="reco::PFRecHit" ClassVersion="18">
+   <version ClassVersion="18" checksum="  1826948750"/>
    <version ClassVersion="17" checksum="2606786437"/>
    <version ClassVersion="10" checksum="1807687081"/>
    <version ClassVersion="11" checksum="3727193351"/>
@@ -49,10 +50,8 @@
     <![CDATA[]]>
    </ioread>
 
-   <field name="positionrep_" transient="true"/>
-   <field name="nCorners_" transient="true"/>
-   <field name="cornersrep_" transient="true"/>
-
+   <field name="caloCell_" transient="true"/>
+   
   </class>
 
 
diff --git a/Geometry/CaloGeometry/interface/CaloCellGeometry.h b/Geometry/CaloGeometry/interface/CaloCellGeometry.h
index 0dafeb092d2fd..17013f79e6bb2 100644
--- a/Geometry/CaloGeometry/interface/CaloCellGeometry.h
+++ b/Geometry/CaloGeometry/interface/CaloCellGeometry.h
@@ -66,15 +66,20 @@ class CaloCellGeometry
   typedef std::vector<ParVec> ParVecVec ;
   typedef EZMgrFL< CCGFloat > ParMgr ;
 
-  enum CornersSize { k_cornerSize = 8 };
+  static constexpr unsigned int k_cornerSize = 8;
+
+  using RepCorners = std::array<RhoEtaPhi,k_cornerSize>;
+  
 
   static const CCGFloat k_ScaleFromDDDtoGeant ;
 
   virtual ~CaloCellGeometry() ;
       
   /// Returns the corner points of this cell's volume.
-  const CornersVec& getCorners() const { assert(not m_corners.uninitialized()); return m_corners; }
+  CornersVec const & getCorners() const { assert(not m_corners.uninitialized()); return m_corners; }
+  RepCorners const & getCornersREP() const {  return m_repCorners;}
 
+  
   /// Returns the position of reference for this cell 
   const GlobalPoint& getPosition() const {return m_refPoint;}
   const GlobalPoint& getBackPoint() const {return m_backPoint;} 
@@ -143,7 +148,7 @@ class CaloCellGeometry
                               0.25 * (cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z()));   
   }
   void initReps() {
-    for (int i=0;i<k_cornerSize; ++i) m_repCorners[i]= {getCorners()[i].perp(), getCorners()[i].eta(), getCorners()[i].barePhi()}; 
+    for (auto i=0U;i<k_cornerSize; ++i) m_repCorners[i]= {getCorners()[i].perp(), getCorners()[i].eta(), getCorners()[i].barePhi()}; 
       
   }
   

From 5f2fef0eb7210f4997bbc24ec8b20b452cadc096 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sat, 16 Apr 2016 15:13:47 +0200
Subject: [PATCH 04/16] RecoParticleFlow compiles, need to fix HF z corr

---
 .../ParticleFlowReco/interface/PFCluster.h    |  2 +-
 .../ParticleFlowReco/interface/PFRecHit.h     | 26 +++---
 DataFormats/ParticleFlowReco/src/PFRecHit.cc  | 16 ++--
 .../interface/PFEcalRecHitCreator.h           | 53 ++---------
 .../interface/PFEcalRecHitCreatorMaxSample.h  | 58 ++-----------
 .../interface/PFHBHERecHitCreator.h           | 30 ++-----
 .../interface/PFHBHERecHitCreatorMaxSample.h  | 20 +----
 .../interface/PFHFRecHitCreator.h             | 61 +++++--------
 .../interface/PFHGCalRecHitCreator.h          | 17 +---
 .../interface/PFHcalRecHitCreator.h           | 40 +++------
 .../interface/PFPSRecHitCreator.h             | 28 +-----
 .../interface/PFRecHitCaloNavigatorWithTime.h |  9 +-
 .../plugins/PFCTRecHitProducer.cc             | 25 ++----
 .../src/Basic2DGenericPFlowClusterizer.cc     |  2 +-
 .../src/Basic2DGenericPFlowPositionCalc.cc    |  2 +-
 .../src/PFlow2DClusterizerWithTime.cc         |  2 +-
 RecoParticleFlow/PFClusterTools/BuildFile.xml |  1 +
 .../PFClusterTools/src/LinkByRecHit.cc        | 87 +++++++++----------
 .../PFClusterTools/src/PFClusterWidthAlgo.cc  | 26 +++---
 .../PFClusterTools/src/PFPhotonClusters.cc    | 12 +--
 .../plugins/kdtrees/KDTreeLinkerPSEcal.cc     | 23 +++--
 .../plugins/kdtrees/KDTreeLinkerTrackEcal.cc  | 49 +++++------
 .../plugins/kdtrees/KDTreeLinkerTrackHcal.cc  | 51 ++++++-----
 RecoParticleFlow/PFProducer/src/PFAlgo.cc     |  2 +-
 24 files changed, 215 insertions(+), 427 deletions(-)

diff --git a/DataFormats/ParticleFlowReco/interface/PFCluster.h b/DataFormats/ParticleFlowReco/interface/PFCluster.h
index 97b9fe6c5d66c..81f01565ce75c 100644
--- a/DataFormats/ParticleFlowReco/interface/PFCluster.h
+++ b/DataFormats/ParticleFlowReco/interface/PFCluster.h
@@ -90,7 +90,7 @@ namespace reco {
     void         setDepth(double depth) {depth_ = depth;}
     
     /// cluster position: rho, eta, phi
-    const REPPoint&       positionREP() const {return posrep_;}
+    const REPPoint& positionREP() const {return posrep_;}
     
     /// computes posrep_ once and for all
     void calculatePositionREP() {
diff --git a/DataFormats/ParticleFlowReco/interface/PFRecHit.h b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
index e88503e824b63..fcf53b8caaa9b 100644
--- a/DataFormats/ParticleFlowReco/interface/PFRecHit.h
+++ b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
@@ -31,7 +31,7 @@ namespace reco {
   class PFRecHit {
 
   public:
-
+    using PositionType = GlobalPoint::BasicVectorType;
     using REPPoint = RhoEtaPhi;
     using RepCorners = CaloCellGeometry::RepCorners;
     using REPPointVector = RepCorners;
@@ -89,8 +89,8 @@ namespace reco {
 
 
     /// calo cell
-    CaloCellGeometry const & callCell() const { return  *caloCell_; }
-    bool hasCaloCel() const { return caloCell_; }
+    CaloCellGeometry const & caloCell() const { return  *caloCell_; }
+    bool hasCaloCell() const { return caloCell_; }
     
     /// rechit detId
     unsigned detId() const {return detId_;}
@@ -110,18 +110,18 @@ namespace reco {
 
     /// rechit momentum transverse to the beam, squared.
     double pt2() const { return energy_ * energy_ *
-	( position().basicVector().perp2()/ position().basicVector().mag2());}
+	( position().perp2()/ position().mag2());}
 
 
     /// rechit cell centre x, y, z
-    GlobalPoint const & position() const { return callCell().getPosition(); }
+    PositionType const & position() const { return caloCell().getPosition().basicVector(); }
     
-    RhoEtaPhi const &  positionREP() const { return callCell().repPos(); }
+    RhoEtaPhi const &  positionREP() const { return caloCell().repPos(); }
 
     /// rechit corners
-    CornersVec const & getCornersXYZ() const { return callCell().getCorners(); }    
+    CornersVec const & getCornersXYZ() const { return caloCell().getCorners(); }    
 
-    RepCorners const & getCornersREP() const { return callCell().getCornersREP();}
+    RepCorners const & getCornersREP() const { return caloCell().getCornersREP();}
  
  
     /// comparison >= operator
@@ -136,14 +136,13 @@ namespace reco {
     /// comparison < operator
     bool operator< (const PFRecHit& rhs) const { return (energy_< rhs.energy_); }
 
-    friend std::ostream& operator<<(std::ostream& out, 
-                                    const reco::PFRecHit& hit);
-
+ 
   private:
 
+    /// cell geometry
     CaloCellGeometry const * caloCell_=nullptr;
  
-    ///C cell detid - should be detid or index in collection ?
+    ///cell detid
     unsigned  int        detId_=0;             
 
     /// rechit layer
@@ -168,6 +167,7 @@ namespace reco {
     PFRecHitRefVector   neighbours8_;
   };
 
-
 }
+std::ostream& operator<<(std::ostream& out, const reco::PFRecHit& hit);
+
 #endif
diff --git a/DataFormats/ParticleFlowReco/src/PFRecHit.cc b/DataFormats/ParticleFlowReco/src/PFRecHit.cc
index 11ad16a3f367c..f3a80948548a6 100644
--- a/DataFormats/ParticleFlowReco/src/PFRecHit.cc
+++ b/DataFormats/ParticleFlowReco/src/PFRecHit.cc
@@ -1,8 +1,6 @@
 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
 
-using namespace reco;
-using namespace std;
-
+namespace reco {
 
 void PFRecHit::addNeighbour(short x,short y,short z,const PFRecHitRef& ref) {
   //bitmask interface  to accomodate more advanced naighbour finding [i.e in z as well]
@@ -66,17 +64,19 @@ const PFRecHitRef PFRecHit::getNeighbour(short x,short y,short z) {
   return PFRecHitRef();
 }
 
+}
 
-ostream& reco::operator<<(ostream& out, const reco::PFRecHit& hit) {
+std::ostream& operator<<(std::ostream& out, const reco::PFRecHit& hit) {
 
   if(!out) return out;
 
-  auto const & pos = hit.positionREP();
-
   out<<"hit id:"<<hit.detId()
      <<" l:"<<hit.layer()
      <<" E:"<<hit.energy()
-     <<" t:"<<hit.time()
-     <<" rep:"<<pos.rho()<<","<<pos.eta()<<","<<pos.phi()<<"| N:";
+     <<" t:"<<hit.time();
+  if (hit.hasCaloCell()) {
+    auto const & pos = hit.positionREP();
+    out  <<" rep:"<<pos.rho()<<","<<pos.eta()<<","<<pos.phi()<<"|";
+  }
   return out;
 }
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFEcalRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFEcalRecHitCreator.h
index 720d3d0b591f2..60dbbe430b5cf 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFEcalRecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFEcalRecHitCreator.h
@@ -51,14 +51,10 @@ template <typename Geometry,PFLayer::Layer Layer,int Detector>
       iEvent.getByToken(recHitToken_,recHitHandle);
       for(const auto& erh : *recHitHandle ) {      
 	const DetId& detid = erh.detid();
-	double energy = erh.energy();
-	double time = erh.time();
+	auto energy = erh.energy();
+	auto time = erh.time();
 
-	math::XYZVector position;
-	math::XYZVector axis;
-	
-	const CaloCellGeometry *thisCell;
-	thisCell= ecalGeo->getGeometry(detid);
+	const CaloCellGeometry * thisCell= ecalGeo->getGeometry(detid);
   
 	// find rechit geometry
 	if(!thisCell) {
@@ -68,50 +64,11 @@ template <typename Geometry,PFLayer::Layer Layer,int Detector>
 	  continue;
 	}
 
-
-	const auto point  = thisCell->getPosition();
-	position.SetCoordinates ( point.x(),
-				  point.y(),
-				  point.z() );
-  
-	// the axis vector is the difference 
-	const TruncatedPyramid* pyr 
-	  = dynamic_cast< const TruncatedPyramid* > (thisCell);    
-
-
-	if( pyr ) {
-	  auto const point1 = pyr->getPosition(1); 
-	  axis.SetCoordinates( point1.x(), 
-			       point1.y(), 
-			       point1.z() ); 
-
-	  auto const point0 = pyr->getPosition(0); 
-    
-	  math::XYZVector axis0( point0.x(), 
-				 point0.y(), 
-				 point0.z() );
-    
-	  axis -= axis0;    
-	}
-	else continue;
-
-	out->emplace_back( detid.rawId(),Layer,
-			   energy, 
-			   position.x(), position.y(), position.z(), 
-			   axis.x(), axis.y(), axis.z() ); 
+	out->emplace_back(thisCell, detid.rawId(),Layer,
+			   energy); 
 
         auto & rh = out->back();
-	//ECAL has no segmentation so put 1
-	
-	const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
-	assert( corners.size() == 8 );
-
-	rh.setNECorner( corners[0].x(), corners[0].y(),  corners[0].z() );
-	rh.setSECorner( corners[1].x(), corners[1].y(),  corners[1].z() );
-	rh.setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z() );
-	rh.setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z() );
 	
-
 	bool rcleaned = false;
 	bool keep=true;
 
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFEcalRecHitCreatorMaxSample.h b/RecoParticleFlow/PFClusterProducer/interface/PFEcalRecHitCreatorMaxSample.h
index 20215754ed03e..04513fe62802a 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFEcalRecHitCreatorMaxSample.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFEcalRecHitCreatorMaxSample.h
@@ -24,7 +24,7 @@
 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
 
 template <typename Geometry,PFLayer::Layer Layer,int Detector>
-  class PFEcalRecHitCreatorMaxSample :  public  PFRecHitCreatorBase {
+  class PFEcalRecHitCreatorMaxSample final :  public  PFRecHitCreatorBase {
 
  public:  
   PFEcalRecHitCreatorMaxSample(const edm::ParameterSet& iConfig,edm::ConsumesCollector& iC):
@@ -51,14 +51,10 @@ template <typename Geometry,PFLayer::Layer Layer,int Detector>
       iEvent.getByToken(recHitToken_,recHitHandle);
       for(const auto& erh : *recHitHandle ) {      
 	const DetId& detid = erh.detid();
-	double energy = erh.energy();
-	double time = erh.time();
-
-	math::XYZVector position;
-	math::XYZVector axis;
-	
-	const CaloCellGeometry *thisCell;
-	thisCell= ecalGeo->getGeometry(detid);
+	auto energy = erh.energy();
+	auto time = erh.time();
+
+	const CaloCellGeometry *thisCell= ecalGeo->getGeometry(detid);
   
 	// find rechit geometry
 	if(!thisCell) {
@@ -69,49 +65,9 @@ template <typename Geometry,PFLayer::Layer Layer,int Detector>
 	}
 
 
-	auto const point  = thisCell->getPosition();
-
-	position.SetCoordinates (point.x(),
-				 point.y(),
-				 point.z() );
-  
-	// the axis vector is the difference 
-	const TruncatedPyramid* pyr 
-	  = dynamic_cast< const TruncatedPyramid* > (thisCell);    
-
-	if( pyr ) {
-	  auto const point1 = pyr->getPosition(1); 
-	  axis.SetCoordinates( point1.x(), 
-			       point1.y(), 
-			       point1.z() ); 
-
-	  auto const point0 = pyr->getPosition(0); 
-    
-	  math::XYZVector axis0( point0.x(), 
-				 point0.y(), 
-				 point0.z() );
-    
-	  axis -= axis0;    
-	}
-	else continue;
-
-
-	reco::PFRecHit rh( detid.rawId(),Layer,
-			   energy, 
-			   position.x(), position.y(), position.z(), 
-			   axis.x(), axis.y(), axis.z() ); 
-
-
-	//ECAL has no segmentation so put 1
-	
-	const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
-	assert( corners.size() == 8 );
+	reco::PFRecHit rh(thisCell, detid.rawId(),Layer,
+			   energy); 
 
-	rh.setNECorner( corners[0].x(), corners[0].y(),  corners[0].z() );
-	rh.setSECorner( corners[1].x(), corners[1].y(),  corners[1].z() );
-	rh.setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z() );
-	rh.setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z() );
-	
 
 	bool rcleaned = false;
 	bool keep=true;
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHBHERecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFHBHERecHitCreator.h
index 889619d87c682..a032ae9e2651f 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFHBHERecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFHBHERecHitCreator.h
@@ -45,14 +45,11 @@ class PFHBHERecHitCreator :  public  PFRecHitCreatorBase {
 	const HcalDetId& detid = (HcalDetId)erh.detid();
 	HcalSubdetector esd=(HcalSubdetector)detid.subdetId();
 	
-	double energy = erh.energy();
-	double time = erh.time();
-	int depth = detid.depth();
-
-	math::XYZVector position;
-	math::XYZVector axis;
+	auto energy = erh.energy();
+	auto time = erh.time();
+	auto depth = detid.depth();
 	
-	const CaloCellGeometry *thisCell=0;
+	const CaloCellGeometry *thisCell=nullptr;
 	PFLayer::Layer layer = PFLayer::HCAL_BARREL1;
 	switch(esd) {
 	case HcalBarrel:
@@ -76,25 +73,10 @@ class PFHBHERecHitCreator :  public  PFRecHitCreatorBase {
 	  continue;
 	}
 
-	auto const point = thisCell->getPosition();
-	position.SetCoordinates ( point.x(),
-				  point.y(),
-				  point.z() );
-  
-	reco::PFRecHit rh( detid.rawId(),layer,
-			   energy, 
-			   position.x(), position.y(), position.z(), 
-			   0,0,0);
+	reco::PFRecHit rh(thisCell, detid.rawId(),layer,
+			   energy);
 	rh.setTime(time); //Mike: This we will use later
 	rh.setDepth(depth);
-	const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
-	assert( corners.size() == 8 );
-
-	rh.setNECorner( corners[0].x(), corners[0].y(),  corners[0].z());
-	rh.setSECorner( corners[1].x(), corners[1].y(),  corners[1].z());
-	rh.setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z());
-	rh.setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z());
-	
 
 	bool rcleaned = false;
 	bool keep=true;
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHBHERecHitCreatorMaxSample.h b/RecoParticleFlow/PFClusterProducer/interface/PFHBHERecHitCreatorMaxSample.h
index f56c9cfbb8d0f..9977ac9c3cfc5 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFHBHERecHitCreatorMaxSample.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFHBHERecHitCreatorMaxSample.h
@@ -164,8 +164,6 @@ class PFHBHERecHitCreatorMaxSample :  public  PFRecHitCreatorBase {
 	  continue;
 
 	int depth = detid.depth();
-	math::XYZVector position;
-	math::XYZVector axis;
 	
 	const CaloCellGeometry *thisCell=0;
 	PFLayer::Layer layer = PFLayer::HCAL_BARREL1;
@@ -192,27 +190,13 @@ class PFHBHERecHitCreatorMaxSample :  public  PFRecHitCreatorBase {
 	}
 
 
-	auto const point = thisCell->getPosition();
-	position.SetCoordinates ( point.x(),
-				  point.y(),
-				  point.z() );
 
-
-	reco::PFRecHit rh( detid.rawId(),layer,
-			   energy, 
-			   position.x(), position.y(), position.z(), 
-			   0,0,0);
+	reco::PFRecHit rh(thisCell, detid.rawId(),layer,
+			   energy);
 
 	rh.setDepth(depth);
 
 
-	const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
-	assert( corners.size() == 8 );
-
-	rh.setNECorner( corners[0].x(), corners[0].y(),  corners[0].z());
-	rh.setSECorner( corners[1].x(), corners[1].y(),  corners[1].z());
-	rh.setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z());
-	rh.setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z());
 	
 	//	for (unsigned int i=0;i<hitEnergies.size();++i)
 	//	  printf(" %f / %f ,",hitEnergies[i],hitTimes[i]);
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
index 7ee155e669372..ae63dd227c16f 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
@@ -17,7 +17,7 @@
 
 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
 
-class PFHFRecHitCreator :  public  PFRecHitCreatorBase {
+class PFHFRecHitCreator final :  public  PFRecHitCreatorBase {
 
  public:  
   PFHFRecHitCreator(const edm::ParameterSet& iConfig,edm::ConsumesCollector& iC):
@@ -56,16 +56,12 @@ class PFHFRecHitCreator :  public  PFRecHitCreatorBase {
       iEvent.getByToken(recHitToken_,recHitHandle);
       for( const auto& erh : *recHitHandle ) {      
 	const HcalDetId& detid = (HcalDetId)erh.detid();
-	int depth = detid.depth();
+	auto depth = detid.depth();
 
-	double energy = erh.energy();
-	double time = erh.time();
+	auto energy = erh.energy();
+	auto time = erh.time();
 
-	math::XYZVector position;
-	math::XYZVector axis;
-	
-	const CaloCellGeometry *thisCell;
-	thisCell= hcalGeo->getGeometry(detid);
+	const CaloCellGeometry * thisCell= hcalGeo->getGeometry(detid);
 
 	// find rechit geometry
 	if(!thisCell) {
@@ -75,42 +71,31 @@ class PFHFRecHitCreator :  public  PFRecHitCreatorBase {
 	  continue;
 	}
 
-	auto const point =  thisCell->getPosition();
+	// auto const point =  thisCell->getPosition();
 
 
 	PFLayer::Layer layer;
-	double depth_correction;
+	//	double depth_correction;
 	if (depth==1) {
 	  layer = PFLayer::HF_EM;
-          depth_correction = point.z() > 0. ? EM_Depth_ : -EM_Depth_;
+          // depth_correction = point.z() > 0. ? EM_Depth_ : -EM_Depth_;
 	}
 	else {
 	  layer = PFLayer::HF_HAD;
-	  depth_correction = point.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
+	  // depth_correction = point.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
 	}
 
-  
+        /*
 	position.SetCoordinates ( point.x(),
 				  point.y(),
 				  point.z()+depth_correction );
+	*/
 
-
-	reco::PFRecHit rh( detid.rawId(),layer,
-			   energy, 
-			   position.x(), position.y(), position.z(), 
-			   0,0,0);
+	reco::PFRecHit rh(thisCell, detid.rawId(),layer,
+			   energy);
 	rh.setTime(time); 
 	rh.setDepth(depth);
 
-	const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
-	assert( corners.size() == 8 );
-
-	rh.setNECorner( corners[0].x(), corners[0].y(),  corners[0].z()+depth_correction);
-	rh.setSECorner( corners[1].x(), corners[1].y(),  corners[1].z()+depth_correction);
-	rh.setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z()+depth_correction);
-	rh.setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z()+depth_correction);
-	
-
 	bool rcleaned = false;
 	bool keep=true;
 
@@ -149,15 +134,14 @@ class PFHFRecHitCreator :  public  PFRecHitCreatorBase {
 	  lONG=hit.energy();
 	  //find the short hit
 	  HcalDetId shortID (HcalForward, detid.ieta(), detid.iphi(), 2);
-	  const reco::PFRecHit temp(shortID,PFLayer::NONE,0.0,math::XYZPoint(0,0,0),math::XYZVector(0,0,0),std::vector<math::XYZPoint>());
 	  auto found_hit = std::lower_bound(tmpOut.begin(),tmpOut.end(),
-					    temp,
+					    shortID,
 					    [](const reco::PFRecHit& a, 
-					       const reco::PFRecHit& b){
-					      return (HcalDetId)(a.detId()) < (HcalDetId)(b.detId());
+					       HcalDetId b){
+					     return  a.detId() < b.rawId();
 					    });
-	  if( found_hit != tmpOut.end() && (HcalDetId)(found_hit->detId()) == (HcalDetId)(shortID.rawId()) ) {
-	    sHORT = found_hit->energy();
+	if( found_hit != tmpOut.end() && found_hit->detId() == shortID.rawId() ) {
+	  sHORT = found_hit->energy();
 	    //Ask for fraction
 	    double energy = lONG-sHORT;
 
@@ -186,15 +170,14 @@ class PFHFRecHitCreator :  public  PFRecHitCreatorBase {
 	else {
 	  sHORT=hit.energy();
 	  HcalDetId longID (HcalForward, detid.ieta(), detid.iphi(), 1);
-	  const reco::PFRecHit temp(longID,PFLayer::NONE,0.0,math::XYZPoint(0,0,0),math::XYZVector(0,0,0),std::vector<math::XYZPoint>());
 	  auto found_hit = std::lower_bound(tmpOut.begin(),tmpOut.end(),
-					    temp,
+					    longID,
 					    [](const reco::PFRecHit& a, 
-					       const reco::PFRecHit& b){
-					      return (HcalDetId)(a.detId()) < (HcalDetId)(b.detId());
+					       HcalDetId b){
+					      return a.detId() < b.rawId();
 					    });
 	  double energy = 2*sHORT;
-	  if( found_hit != tmpOut.end() && (HcalDetId)(found_hit->detId()) == (HcalDetId)(longID.rawId()) ) {
+	  if( found_hit != tmpOut.end() && found_hit->detId() == longID.rawId() ) {
 	    lONG = found_hit->energy();
 	    //Ask for fraction
 
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHGCalRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFHGCalRecHitCreator.h
index 0ad48a4288c9b..89b3c4b0bc92a 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFHGCalRecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFHGCalRecHitCreator.h
@@ -70,23 +70,12 @@ template <typename DET,PFLayer::Layer Layer,ForwardSubdetector subdet>
 	  continue;
 	}
   
-	const GlobalPoint position( std::move( hgcGeo.getPosition( detid ) ) );
-	//std::cout << "geometry cell position: " << position << std::endl;
 
-	reco::PFRecHit rh( detid.rawId(),Layer,
-			   energy, 
-			   position.x(), position.y(), position.z(), 
-			   0, 0, 0 ); 
+	reco::PFRecHit rh(thisCell, detid.rawId(),Layer,
+			   energy); 
 	
-	rh.setOriginalRecHit(edm::Ref<HGCRecHitCollection>(recHitHandle,i));
+	//  rh.setOriginalRecHit(edm::Ref<HGCRecHitCollection>(recHitHandle,i));
 
-	const HGCalGeometry::CornersVec corners( std::move( hgcGeo.getCorners( detid ) ) );
-	assert( corners.size() == 8 );
-
-	rh.setNECorner( corners[0].x(), corners[0].y(),  corners[0].z() );
-	rh.setSECorner( corners[1].x(), corners[1].y(),  corners[1].z() );
-	rh.setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z() );
-	rh.setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z() );
 	
 	bool rcleaned = false;
 	bool keep=true;
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHcalRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFHcalRecHitCreator.h
index d058f7fdc63e8..5881c4f11f265 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFHcalRecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFHcalRecHitCreator.h
@@ -18,7 +18,7 @@
 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
 
 template <typename Digi, typename Geometry,PFLayer::Layer Layer,int Detector>
-  class PFHcalRecHitCreator :  public  PFRecHitCreatorBase {
+  class PFHcalRecHitCreator final :  public  PFRecHitCreatorBase {
 
  public:  
   PFHcalRecHitCreator(const edm::ParameterSet& iConfig,edm::ConsumesCollector& iC):
@@ -53,15 +53,12 @@ template <typename Digi, typename Geometry,PFLayer::Layer Layer,int Detector>
 	  continue;
 
 
-	double energy = erh.energy();
-	double time = erh.time();
-	int depth =detid.depth();
+	auto energy = erh.energy();
+	auto time = erh.time();
+	auto depth =detid.depth();
 	  
-	math::XYZVector position;
-	math::XYZVector axis;
 	
-	const CaloCellGeometry *thisCell;
-	thisCell= hcalGeo->getGeometry(detid);
+	const CaloCellGeometry * thisCell= hcalGeo->getGeometry(detid);
   
 	// find rechit geometry
 	if(!thisCell) {
@@ -71,30 +68,13 @@ template <typename Digi, typename Geometry,PFLayer::Layer Layer,int Detector>
 	  continue;
 	}
 
-	auto const point  = thisCell->getPosition();
-	position.SetCoordinates ( point.x(),
-				  point.y(),
-				  point.z() );
-  
-
-
 
-	reco::PFRecHit rh( detid.rawId(),Layer,
-			   energy, 
-			   position.x(), position.y(), position.z(), 
-			   0,0,0);
+	reco::PFRecHit rh(thisCell, detid.rawId(),Layer,
+			   energy);
 	rh.setTime(time); //Mike: This we will use later
 	rh.setDepth(depth);
 
-	const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
-	assert( corners.size() == 8 );
-
-	rh.setNECorner( corners[0].x(), corners[0].y(),  corners[0].z());
-	rh.setSECorner( corners[1].x(), corners[1].y(),  corners[1].z());
-	rh.setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z());
-	rh.setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z());
-	
-
+ 
 	bool rcleaned = false;
 	bool keep=true;
 
@@ -107,10 +87,10 @@ template <typename Digi, typename Geometry,PFLayer::Layer Layer,int Detector>
 	}
 	  
 	if(keep) {
-	  out->push_back(rh);
+	  out->push_back(std::move(rh));
 	}
 	else if (rcleaned) 
-	  cleaned->push_back(rh);
+	  cleaned->push_back(std::move(rh));
       }
     }
 
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFPSRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFPSRecHitCreator.h
index 7ac1677794222..ad990e3582e93 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFPSRecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFPSRecHitCreator.h
@@ -24,7 +24,7 @@
 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
 
-class PFPSRecHitCreator :  public  PFRecHitCreatorBase {
+class PFPSRecHitCreator final :  public  PFRecHitCreatorBase {
 
  public:  
   PFPSRecHitCreator(const edm::ParameterSet& iConfig,edm::ConsumesCollector& iC):
@@ -50,9 +50,6 @@ class PFPSRecHitCreator :  public  PFRecHitCreatorBase {
 	ESDetId detid(erh.detid());
 	double energy = erh.energy();
 
-
-	math::XYZVector position;
-
 	PFLayer::Layer layer = PFLayer::NONE;
 	
 	switch( detid.plane() ) {
@@ -70,8 +67,7 @@ class PFPSRecHitCreator :  public  PFRecHitCreatorBase {
  
 
 	
-	const CaloCellGeometry *thisCell;
-	thisCell= psGeometry->getGeometry(detid);
+	const CaloCellGeometry * thisCell= psGeometry->getGeometry(detid);
   
 	// find rechit geometry
 	if(!thisCell) {
@@ -81,28 +77,12 @@ class PFPSRecHitCreator :  public  PFRecHitCreatorBase {
 	  continue;
 	}
 
-	auto const point  = thisCell->getPosition();
-	position.SetCoordinates ( point.x(),
-				  point.y(),
-				  point.z() );
-
-        out->emplace_back(  detid.rawId(),layer,
-			   energy, 
-			   position.x(), position.y(), position.z(), 
-			   0.0,0.0,0.0);
+        out->emplace_back(thisCell, detid.rawId(),layer,
+			   energy);
         auto & rh = out->back();
 	rh.setDepth(detid.plane());
 	rh.setTime(erh.time());
 	
-	const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
-	assert( corners.size() == 8 );
-
-	rh.setNECorner( corners[0].x(), corners[0].y(),  corners[0].z() );
-	rh.setSECorner( corners[1].x(), corners[1].y(),  corners[1].z() );
-	rh.setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z() );
-	rh.setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z() );
-	
-
 	bool rcleaned = false;
 	bool keep=true;
 
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h
index d71da9c889f3a..52dd72d66b381 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h
@@ -119,15 +119,14 @@ class PFRecHitCaloNavigatorWithTime : public PFRecHitNavigatorBase {
   void associateNeighbour(const DetId& id, reco::PFRecHit& hit,std::auto_ptr<reco::PFRecHitCollection>& hits,edm::RefProd<reco::PFRecHitCollection>& refProd,short eta, short phi) {
     double sigma2=10000.0;
     
-    const reco::PFRecHit temp(id,PFLayer::NONE,0.0,math::XYZPoint(0,0,0),math::XYZVector(0,0,0),std::vector<math::XYZPoint>());
 
     auto found_hit = std::lower_bound(hits->begin(),hits->end(),
-				      temp,
+				      id,
 				      [](const reco::PFRecHit& a, 
-					 const reco::PFRecHit& b){
-					if (DetId(a.detId()).det() == DetId::Hcal || DetId(b.detId()).det() == DetId::Hcal) return (HcalDetId)(a.detId()) < (HcalDetId)(b.detId());
+					  DetId b){
+					if (DetId(a.detId()).det() == DetId::Hcal || b.det() == DetId::Hcal) return (HcalDetId)(a.detId()) < (HcalDetId)(b);
 
-					else return a.detId() < b.detId();
+					else return a.detId() < b.rawId();
 				      });
 
 
diff --git a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
index 682c48cdb1b4a..83bf89b7155e7 100644
--- a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
+++ b/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
@@ -856,36 +856,21 @@ PFCTRecHitProducer::createHcalRecHit( const DetId& detid,
       <<layer<<endl;
     return 0;
   }
-  
-  const GlobalPoint& position = thisCell->getPosition();
-  
-  double depth_correction = 0.;
+    
+  // double depth_correction = 0.;
   switch ( layer ) { 
   case PFLayer::HF_EM:
-    depth_correction = position.z() > 0. ? EM_Depth_ : -EM_Depth_;
+    // depth_correction = position.z() > 0. ? EM_Depth_ : -EM_Depth_;
     break;
   case PFLayer::HF_HAD:
-    depth_correction = position.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
+    //  depth_correction = position.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
     break;
   default:
     break;
   }
 
   reco::PFRecHit *rh = 
-    new reco::PFRecHit( newDetId.rawId(),  layer, energy, 
-			position.x(), position.y(), position.z()+depth_correction, 
-			0,0,0 );
- 
-  
-  
-  
-  // set the corners
-  const CaloCellGeometry::CornersVec& corners = thisCell->getCorners();
-
-  rh->setNECorner( corners[0].x(), corners[0].y(),  corners[0].z()+depth_correction );
-  rh->setSECorner( corners[1].x(), corners[1].y(),  corners[1].z()+depth_correction );
-  rh->setSWCorner( corners[2].x(), corners[2].y(),  corners[2].z()+depth_correction );
-  rh->setNWCorner( corners[3].x(), corners[3].y(),  corners[3].z()+depth_correction );
+    new reco::PFRecHit(thisCell, newDetId.rawId(),  layer, energy);
  
   return rh;
 }
diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowClusterizer.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowClusterizer.cc
index fbf33a0eac83f..8e670674e08df 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowClusterizer.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowClusterizer.cc
@@ -168,7 +168,7 @@ growPFClusters(const reco::PFCluster& topo,
     }  
     const double recHitEnergyNorm = 
       _recHitEnergyNorms.find(cell_layer)->second; 
-    const math::XYZPoint& topocellpos_xyz = refhit->position();
+    math::XYZPoint topocellpos_xyz(refhit->position());
     dist2.clear(); frac.clear(); fractot = 0;
     // add rechits to clusters, calculating fraction based on distance
     for( auto& cluster : clusters ) {      
diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
index ef168d790da7e..a12a20882757f 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
@@ -118,7 +118,7 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
     const double norm = ( rhf.fraction < _minFractionInCalc ? 
 			  0.0 : 
 			  std::max(0.0,vdt::fast_log(rh_energy*_logWeightDenom)) );
-    const math::XYZPoint& rhpos_xyz = refhit.position();
+    const math::XYZPoint rhpos_xyz(refhit.position());
     x += rhpos_xyz.X() * norm;
     y += rhpos_xyz.Y() * norm;
     z += rhpos_xyz.Z() * norm;
diff --git a/RecoParticleFlow/PFClusterProducer/src/PFlow2DClusterizerWithTime.cc b/RecoParticleFlow/PFClusterProducer/src/PFlow2DClusterizerWithTime.cc
index 6201c09c4ceb6..6748d4ea8dd37 100644
--- a/RecoParticleFlow/PFClusterProducer/src/PFlow2DClusterizerWithTime.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/PFlow2DClusterizerWithTime.cc
@@ -206,7 +206,7 @@ growPFClusters(const reco::PFCluster& topo,
     }  
     const double recHitEnergyNorm = 
       _recHitEnergyNorms.find(cell_layer)->second; 
-    const math::XYZPoint& topocellpos_xyz = refhit->position();
+    math::XYZPoint topocellpos_xyz(refhit->position());
     dist2.clear(); frac.clear(); fractot = 0;
 
     // add rechits to clusters, calculating fraction based on distance
diff --git a/RecoParticleFlow/PFClusterTools/BuildFile.xml b/RecoParticleFlow/PFClusterTools/BuildFile.xml
index c21137bea2760..bd1e922d39828 100644
--- a/RecoParticleFlow/PFClusterTools/BuildFile.xml
+++ b/RecoParticleFlow/PFClusterTools/BuildFile.xml
@@ -2,6 +2,7 @@
 <use   name="DataFormats/ParticleFlowReco"/>
 <use   name="DataFormats/Math"/>
 <use   name="DataFormats/EcalDetId"/>
+<use   name="Geometry/CaloGeometry"/>
 <use   name="boost"/>
 <use   name="rootphysics"/>
 
diff --git a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
index e94578f2db78a..2a08acc5647ff 100644
--- a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
+++ b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
@@ -18,11 +18,9 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
 #endif
   
   //cluster position
-  double clustereta  = cluster.positionREP().Eta();
-  double clusterphi  = cluster.positionREP().Phi();
-  //double clusterX    = cluster.position().X();
-  //double clusterY    = cluster.position().Y();
-  double clusterZ    = cluster.position().Z();
+  auto clustereta  = cluster.positionREP().Eta();
+  auto clusterphi  = cluster.positionREP().Phi();
+  auto clusterZ    = cluster.position().Z();
 
   bool barrel = false;
   bool hcal = false;
@@ -53,8 +51,7 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
   case PFLayer::ECAL_ENDCAP:
 #ifdef PFLOW_DEBUG
     if( debug )
-      std::cout << "Fetching Ecal Resolution Maps"
-	   << std::endl;
+      std::cout << "Fetching Ecal Resolution Maps"<< std::endl;
 #endif
     // did not reach ecal, cannot be associated with a cluster.
     if( ! atECAL.isValid() ) return -1.;   
@@ -77,8 +74,7 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
   case PFLayer::HCAL_ENDCAP:  
 #ifdef PFLOW_DEBUG
     if( debug )
-      std::cout << "Fetching Hcal Resolution Maps"
-	   << std::endl;
+      std::cout << "Fetching Hcal Resolution Maps" << std::endl;
 #endif
     if( isBrem ) {  
       return  -1.;
@@ -113,8 +109,7 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
   case PFLayer::HCAL_BARREL2: barrel = true; 
 #ifdef PFLOW_DEBUG
     if( debug )
-     std::cout << "Fetching HO Resolution Maps"
-	       << std::endl;
+     std::cout << "Fetching HO Resolution Maps" << std::endl;
 #endif
     if( isBrem ) {  
       return  -1.;
@@ -217,18 +212,17 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
     if(rh.isNull()) continue;
     
     //getting rechit center position
-    const reco::PFRecHit& rechit_cluster = *rh;
-    const math::XYZPoint& posxyz 
+    const auto & rechit_cluster = *rh;
+    const auto & posxyz 
       = rechit_cluster.position();
-    const reco::PFRecHit::REPPoint& posrep 
+    const auto & posrep 
       = rechit_cluster.positionREP();
     
     //getting rechit corners
-    const std::vector< math::XYZPoint >& 
+    const auto & 
       cornersxyz = rechit_cluster.getCornersXYZ();
-    const std::vector<reco::PFRecHit::REPPoint>& corners = 
-      rechit_cluster.getCornersREP();
-    assert(corners.size() == 4);
+    const auto & corners = rechit_cluster.getCornersREP();
+
     
     if( barrel || hcal ){ // barrel case matching in eta/phi 
                           // (and HCAL endcap too!)
@@ -237,13 +231,13 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
       // blown up by 50% (HCAL) to 100% (ECAL) to include cracks & gaps
       // also blown up to account for multiple scattering at low pt.
       double rhsizeEta 
-	= fabs(corners[0].Eta() - corners[2].Eta());
+	= std::abs(corners[0].eta() - corners[2].eta());
       double rhsizePhi 
-	= fabs(corners[0].Phi() - corners[2].Phi());
+	= std::abs(corners[0].phi() - corners[2].phi());
       if ( rhsizePhi > M_PI ) rhsizePhi = 2.*M_PI - rhsizePhi;
       if ( hcal ) { 
 	const double mult = horesolscale * (1.50 + 0.5/fracs.size());
-	rhsizeEta = rhsizeEta * mult + 0.2*fabs(dHEta);
+	rhsizeEta = rhsizeEta * mult + 0.2*std::abs(dHEta);
 	rhsizePhi = rhsizePhi * mult + 0.2*fabs(dHPhi); 
 	
       } else { 
@@ -260,8 +254,8 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
 	     << rechit_cluster.energy() 
 	     << std::endl; 
 	for ( unsigned jc=0; jc<4; ++jc ) 
-	  std::cout<<"corners "<<jc<<" "<<corners[jc].Eta()
-	      <<" "<<corners[jc].Phi()<<std::endl;
+	  std::cout<<"corners "<<jc<<" "<<corners[jc].eta()
+	      <<" "<<corners[jc].phi()<<std::endl;
 	
 	std::cout << "RecHit SizeEta=" << rhsizeEta
 	     << " SizePhi=" << rhsizePhi << std::endl;
@@ -271,8 +265,8 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
       //distance track-rechit center
       // const math::XYZPoint& posxyz 
       // = rechit_cluster.position();
-      double deta = fabs(posrep.Eta() - tracketa);
-      double dphi = fabs(posrep.Phi() - trackphi);
+      double deta = fabs(posrep.eta() - tracketa);
+      double dphi = fabs(posrep.phi() - trackphi);
       if ( dphi > M_PI ) dphi = 2.*M_PI - dphi;
       
 #ifdef PFLOW_DEBUG
@@ -295,11 +289,11 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
       
 #ifdef PFLOW_DEBUG
       if( debug ){
-	const math::XYZPoint& posxyz 
+	const auto & posxyz 
 	  = rechit_cluster.position();
 	
-	std::cout << "RH " << posxyz.X()
-	     << " "   << posxyz.Y()
+	std::cout << "RH " << posxyz.x()
+	     << " "   << posxyz.y()
 	     << std::endl;
 	
 	std::cout << "TRACK " << track_X
@@ -312,16 +306,16 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
       double y[5];
       
       for ( unsigned jc=0; jc<4; ++jc ) {
-	const math::XYZPoint& cornerposxyz = cornersxyz[jc];
+	const auto & cornerposxyz = cornersxyz[jc];
 	const double mult = (1.00+0.50/(fracs.size()*std::min(1.,0.5*trackPt)));
-	x[jc] = cornerposxyz.X() + (cornerposxyz.X()-posxyz.X()) * mult;
-	y[jc] = cornerposxyz.Y() + (cornerposxyz.Y()-posxyz.Y()) * mult;
+	x[jc] = cornerposxyz.x() + (cornerposxyz.y()-posxyz.x()) * mult;
+	y[jc] = cornerposxyz.y() + (cornerposxyz.y()-posxyz.y()) * mult;
 	
 #ifdef PFLOW_DEBUG
 	if( debug ){
 	  std::cout<<"corners "<<jc
-	      << " " << cornerposxyz.X()
-	      << " " << cornerposxyz.Y()
+	      << " " << cornerposxyz.x()
+	      << " " << cornerposxyz.y()
 	      << std::endl;
 	}
 #endif
@@ -433,7 +427,7 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
   //loop rechits
   for(unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
 
-    const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
+    const auto & rh = fracs[rhit].recHitRef();
     double fraction = fracs[rhit].fraction();
     if(fraction < 1E-4) continue;
     if(rh.isNull()) continue;
@@ -442,10 +436,9 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
     const reco::PFRecHit& rechit_cluster = *rh;
     
     //getting rechit corners
-    const std::vector< math::XYZPoint >&  corners = rechit_cluster.getCornersXYZ();
-    assert(corners.size() == 4);
+    const auto &  corners = rechit_cluster.getCornersXYZ();
     
-    const math::XYZPoint& posxyz = rechit_cluster.position() * zPS/zECAL;
+    auto posxyz = rechit_cluster.position() * zPS/zECAL;
 #ifdef PFLOW_DEBUG
     if( debug ){
       std::cout << "Ecal rechit " << posxyz.X() << " "   << posxyz.Y() << std::endl;
@@ -457,16 +450,16 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
     double y[5];
     for ( unsigned jc=0; jc<4; ++jc ) {
       // corner position projected onto the preshower
-      math::XYZPoint cornerpos = corners[jc] * zPS/zECAL;
+      auto cornerpos = corners[jc].basicVector() * zPS/zECAL;
       // Inflate the size by the size of the PS strips, and by 5% to include ECAL cracks.
-      x[jc] = cornerpos.X() + (cornerpos.X()-posxyz.X()) * (0.05 +1.0/fabs((cornerpos.X()-posxyz.X()))*0.5*deltaX);
-      y[jc] = cornerpos.Y() + (cornerpos.Y()-posxyz.Y()) * (0.05 +1.0/fabs((cornerpos.Y()-posxyz.Y()))*0.5*deltaY);
+      x[jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/fabs((cornerpos.x()-posxyz.x()))*0.5*deltaX);
+      y[jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/fabs((cornerpos.y()-posxyz.y()))*0.5*deltaY);
       
 #ifdef PFLOW_DEBUG
       if( debug ){
 	std::cout<<"corners "<<jc
-	    << " " << cornerpos.X() << " " << x[jc] 
-	    << " " << cornerpos.Y() << " " << y[jc]
+	    << " " << cornerpos.x() << " " << x[jc] 
+	    << " " << cornerpos.y() << " " << y[jc]
 	    << std::endl;
       }
 #endif
@@ -489,7 +482,9 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
   }//loop rechits
   
   if( linkedbyrechit ) {
+#ifdef PFLOW_DEBUG
     if( debug ) std::cout << "Cluster PS and Cluster ECAL LINKED BY RECHIT" << std::endl;
+#endif
     double dist = computeDist( xECAL/1000.,yECAL/1000.,
 			       xPS/1000.  ,yPS/1000, 
 			       false);    
@@ -507,8 +502,8 @@ LinkByRecHit::testHFEMAndHFHADByRecHit(const reco::PFCluster& clusterHFEM,
 				      const reco::PFCluster& clusterHFHAD,
 				      bool debug) {
   
-  math::XYZPoint posxyzEM = clusterHFEM.position();
-  math::XYZPoint posxyzHAD = clusterHFHAD.position();
+  auto posxyzEM = clusterHFEM.position();
+  auto posxyzHAD = clusterHFHAD.position();
 
   double dX = posxyzEM.X()-posxyzHAD.X();
   double dY = posxyzEM.Y()-posxyzHAD.Y();
@@ -533,8 +528,8 @@ LinkByRecHit::computeDist( double eta1, double phi1,
 			   double eta2, double phi2,
 			   bool etaPhi )  {
   
-  const double phicor = etaPhi ? normalizedPhi(phi1 - phi2) : phi1 - phi2;
-  const double etadiff = eta1 - eta2;
+  auto phicor = etaPhi ? normalizedPhi(phi1 - phi2) : phi1 - phi2;
+  auto etadiff = eta1 - eta2;
   
   // double chi2 =  
   //  (eta1 - eta2)*(eta1 - eta2) / ( reta1*reta1+ reta2*reta2 ) +
diff --git a/RecoParticleFlow/PFClusterTools/src/PFClusterWidthAlgo.cc b/RecoParticleFlow/PFClusterTools/src/PFClusterWidthAlgo.cc
index ad19da5669b9e..0671d50afb757 100644
--- a/RecoParticleFlow/PFClusterTools/src/PFClusterWidthAlgo.cc
+++ b/RecoParticleFlow/PFClusterTools/src/PFClusterWidthAlgo.cc
@@ -42,9 +42,9 @@ PFClusterWidthAlgo::PFClusterWidthAlgo(const std::vector<const reco::PFCluster *
         double energyHit = RefPFRecHit->energy()*it->fraction();
 
         sclusterE += energyHit;
-        posX += energyHit*RefPFRecHit->position().X();
-        posY += energyHit*RefPFRecHit->position().Y();
-        posZ += energyHit*RefPFRecHit->position().Z();
+        posX += energyHit*RefPFRecHit->position().x();
+        posY += energyHit*RefPFRecHit->position().y();
+        posZ += energyHit*RefPFRecHit->position().z();
       
       }
     } // end for ncluster    
@@ -66,9 +66,9 @@ PFClusterWidthAlgo::PFClusterWidthAlgo(const std::vector<const reco::PFCluster *
 
     //second loop, compute variances
     for(unsigned int icl=0; icl<nclust; ++icl) {
-      const std::vector< reco::PFRecHitFraction >& PFRecHits =  pfclust[icl]->recHitFractions();  
+      const auto & PFRecHits =  pfclust[icl]->recHitFractions();  
       
-      for ( std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); 
+      for ( auto it = PFRecHits.begin(); 
 	    it != PFRecHits.end(); ++it) {
 	const PFRecHitRef& RefPFRecHit = it->recHitRef(); 
         //compute rechit energy taking into account fractions
@@ -83,30 +83,30 @@ PFClusterWidthAlgo::PFClusterWidthAlgo(const std::vector<const reco::PFCluster *
 	  }
 	}
 
-	double dPhi = reco::deltaPhi(RefPFRecHit->position().phi(),scPhi);
-	double dEta = RefPFRecHit->position().eta() - scEta;
+	double dPhi = reco::deltaPhi(RefPFRecHit->positionREP().phi(),scPhi);
+	double dEta = RefPFRecHit->positionREP().eta() - scEta;
 	numeratorEtaWidth += energyHit * dEta * dEta;
 	numeratorPhiWidth += energyHit * dPhi * dPhi;
       }
     } // end for ncluster
 
     //for the first cluster (from GSF) computed sigmaEtaEta
-    const std::vector< reco::PFRecHitFraction >& PFRecHits =  pfclust[0]->recHitFractions();
-    for ( std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); 
+    const auto & PFRecHits =  pfclust[0]->recHitFractions();
+    for ( auto it = PFRecHits.begin(); 
 	  it != PFRecHits.end(); ++it) {
-      const PFRecHitRef& RefPFRecHit = it->recHitRef(); 
+      const auto & RefPFRecHit = it->recHitRef(); 
       if(!RefPFRecHit.isAvailable()) 
 	return;
       double energyHit = RefPFRecHit->energy();
       if (RefPFRecHit->detId() != SeedDetID) {
-	float diffEta =  RefPFRecHit->position().eta() - SeedEta;
+	float diffEta =  RefPFRecHit->positionREP().eta() - SeedEta;
 	sigmaEtaEta_ += (diffEta*diffEta) * (energyHit/SeedClusEnergy);
       }
     }
     if (sigmaEtaEta_ == 0.) sigmaEtaEta_ = 0.00000001;
 
-    etaWidth_ = sqrt(numeratorEtaWidth / denominator);
-    phiWidth_ = sqrt(numeratorPhiWidth / denominator);
+    etaWidth_ = std::sqrt(numeratorEtaWidth / denominator);
+    phiWidth_ = std::sqrt(numeratorPhiWidth / denominator);
     
 
   } // endif ncluster > 0
diff --git a/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc b/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc
index 30bcb3d8bfec3..d2e83672a7f34 100644
--- a/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc
+++ b/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc
@@ -2,7 +2,7 @@
 #include "DataFormats/EcalDetId/interface/EEDetId.h"
 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
 #include "RecoParticleFlow/PFClusterTools/interface/PFPhotonClusters.h"
-
+#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
 #include <TMath.h>
 #include <TVector2.h>
 using namespace reco;
@@ -24,17 +24,17 @@ void PFPhotonClusters::SetSeed(){
   math::XYZVector axis;
   math::XYZVector position;
   DetId idseed;
-  const std::vector< reco::PFRecHitFraction >& PFRecHits=
+  const auto & PFRecHits=
     PFClusterRef_->recHitFractions();
   
-  for(std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin();
+  for(auto it = PFRecHits.begin();
       it != PFRecHits.end(); ++it){
-    const PFRecHitRef& RefPFRecHit = it->recHitRef();
-    double frac=it->fraction();
+    const auto & RefPFRecHit = it->recHitRef();
+    auto frac=it->fraction();
     float E= RefPFRecHit->energy()* frac;
     if(E>PFSeedE){
       PFSeedE=E;  
-      axis=RefPFRecHit->getAxisXYZ();
+      axis=  dynamic_cast<TruncatedPyramid const &>(RefPFRecHit->caloCell()).axis();
       position=RefPFRecHit->position();
       idseed = RefPFRecHit->detId();
     }
diff --git a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc
index fed08dfdc98d4..22423d051b274 100644
--- a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc
+++ b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc
@@ -89,9 +89,9 @@ KDTreeLinkerPSEcal::buildTree(const RecHitSet	&rechitsSet,
       it != rechitsSet.end(); it++) {
 
     const reco::PFRecHit* rh = *it;
-    const math::XYZPoint& posxyz = rh->position();
+    const auto & posxyz = rh->position();
         
-    KDTreeNodeInfo rhinfo (rh, posxyz.X(), posxyz.Y());
+    KDTreeNodeInfo rhinfo (rh, posxyz.x(), posxyz.y());
     eltList.push_back(rhinfo);
   }
 
@@ -171,8 +171,7 @@ KDTreeLinkerPSEcal::searchLinks()
     for(std::vector<KDTreeNodeInfo>::const_iterator rhit = recHits.begin(); 
 	rhit != recHits.end(); ++rhit) {
            
-      const std::vector< math::XYZPoint >& corners = rhit->ptr->getCornersXYZ();
-      if(corners.size() != 4) continue;
+      const auto & corners = rhit->ptr->getCornersXYZ();
 
       // Find all clusters associated to given rechit
       RecHit2BlockEltMap::iterator ret = rechit2ClusterLinks_.find(rhit->ptr);
@@ -181,16 +180,16 @@ KDTreeLinkerPSEcal::searchLinks()
 	  clusterIt != ret->second.end(); clusterIt++) {
 	
 	reco::PFClusterRef clusterref = (*clusterIt)->clusterRef();
-	double clusterz = clusterref->position().Z();
+	double clusterz = clusterref->position().z();
 
-	const math::XYZPoint& posxyz = rhit->ptr->position() * zPS / clusterz;
+	const auto & posxyz = rhit->ptr->position() * zPS / clusterz;
 
 	double x[5];
 	double y[5];
 	for ( unsigned jc=0; jc<4; ++jc ) {
-	  math::XYZPoint cornerpos = corners[jc] * zPS / clusterz;
-	  x[jc] = cornerpos.X() + (cornerpos.X()-posxyz.X()) * (0.05 +1.0/fabs((cornerpos.X()-posxyz.X()))*deltaX/2.);
-	  y[jc] = cornerpos.Y() + (cornerpos.Y()-posxyz.Y()) * (0.05 +1.0/fabs((cornerpos.Y()-posxyz.Y()))*deltaY/2.);
+	  auto cornerpos = corners[jc].basicVector() * zPS / clusterz;
+	  x[jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/fabs((cornerpos.x()-posxyz.x()))*deltaX/2.);
+	  y[jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/fabs((cornerpos.y()-posxyz.y()))*deltaY/2.);
 	}
 
 	x[4] = x[0];
@@ -222,10 +221,10 @@ KDTreeLinkerPSEcal::updatePFBlockEltWithLinks()
     for (BlockEltSet::iterator jt = it->second.begin();
 	 jt != it->second.end(); ++jt) {
 
-      double clusterPhi = (*jt)->clusterRef()->positionREP().Phi();
-      double clusterEta = (*jt)->clusterRef()->positionREP().Eta();
+      double clusterphi = (*jt)->clusterRef()->positionREP().phi();
+      double clustereta = (*jt)->clusterRef()->positionREP().eta();
 
-      multitracks.linkedClusters.push_back(std::make_pair(clusterPhi, clusterEta));
+      multitracks.linkedClusters.push_back(std::make_pair(clusterphi, clustereta));
     }
 
     it->first->setMultilinks(multitracks);
diff --git a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc
index 86af2180c2e21..3dfaaeba31155 100644
--- a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc
+++ b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc
@@ -71,20 +71,20 @@ KDTreeLinkerTrackEcal::buildTree()
     
     const reco::PFRecHit::REPPoint &posrep = (*it)->positionREP();
     
-    KDTreeNodeInfo rh1 (*it, posrep.Eta(), posrep.Phi());
+    KDTreeNodeInfo rh1 (*it, posrep.eta(), posrep.phi());
     eltList.push_back(rh1);
     
     // Here we solve the problem of phi circular set by duplicating some rechits
     // too close to -Pi (or to Pi) and adding (substracting) to them 2 * Pi.
     if (rh1.dim2 > (M_PI - getPhiOffset())) {
       double phi = rh1.dim2 - 2 * M_PI;
-      KDTreeNodeInfo rh2(*it, posrep.Eta(), phi); 
+      KDTreeNodeInfo rh2(*it, posrep.eta(), phi); 
       eltList.push_back(rh2);
     }
 
     if (rh1.dim2 < (M_PI * -1.0 + getPhiOffset())) {
       double phi = rh1.dim2 + 2 * M_PI;
-      KDTreeNodeInfo rh3(*it, posrep.Eta(), phi); 
+      KDTreeNodeInfo rh3(*it, posrep.eta(), phi); 
       eltList.push_back(rh3);
     }
   }
@@ -125,8 +125,8 @@ KDTreeLinkerTrackEcal::searchLinks()
       trackref->extrapolatedPoint( reco::PFTrajectoryPoint::ClosestApproach );
     
     double trackPt = sqrt(atVertex.momentum().Vect().Perp2());
-    double tracketa = atECAL.positionREP().Eta();
-    double trackphi = atECAL.positionREP().Phi();
+    double tracketa = atECAL.positionREP().eta();
+    double trackphi = atECAL.positionREP().phi();
     double trackx = atECAL.position().X();
     double tracky = atECAL.position().Y();
     double trackz = atECAL.position().Z();
@@ -144,18 +144,17 @@ KDTreeLinkerTrackEcal::searchLinks()
     for(std::vector<KDTreeNodeInfo>::const_iterator rhit = recHits.begin(); 
 	rhit != recHits.end(); ++rhit) {
            
-      const std::vector< math::XYZPoint >& cornersxyz      = rhit->ptr->getCornersXYZ();
-      const math::XYZPoint& posxyz			   = rhit->ptr->position();
-      const reco::PFRecHit::REPPoint &rhrep		   = rhit->ptr->positionREP();
-      const std::vector<reco::PFRecHit::REPPoint>& corners = rhit->ptr->getCornersREP();
-      if(corners.size() != 4) continue;
+      const auto & cornersxyz      = rhit->ptr->getCornersXYZ();
+      const auto & posxyz			   = rhit->ptr->position();
+      const auto &rhrep		   = rhit->ptr->positionREP();
+      const auto & corners = rhit->ptr->getCornersREP();
       
-      double rhsizeEta = fabs(corners[0].Eta() - corners[2].Eta());
-      double rhsizePhi = fabs(corners[0].Phi() - corners[2].Phi());
-      if ( rhsizePhi > M_PI ) rhsizePhi = 2.*M_PI - rhsizePhi;
+      double rhsizeeta = fabs(corners[0].eta() - corners[2].eta());
+      double rhsizephi = fabs(corners[0].phi() - corners[2].phi());
+      if ( rhsizephi > M_PI ) rhsizephi = 2.*M_PI - rhsizephi;
       
-      double deta = fabs(rhrep.Eta() - tracketa);
-      double dphi = fabs(rhrep.Phi() - trackphi);
+      double deta = fabs(rhrep.eta() - tracketa);
+      double dphi = fabs(rhrep.phi() - trackphi);
       if ( dphi > M_PI ) dphi = 2.*M_PI - dphi;
       
       // Find all clusters associated to given rechit
@@ -165,18 +164,18 @@ KDTreeLinkerTrackEcal::searchLinks()
 	  clusterIt != ret->second.end(); clusterIt++) {
 	
 	reco::PFClusterRef clusterref = (*clusterIt)->clusterRef();
-	double clusterz = clusterref->position().Z();
+	double clusterz = clusterref->position().z();
 	int fracsNbr = clusterref->recHitFractions().size();
 
 	if (clusterref->layer() == PFLayer::ECAL_BARREL){ // BARREL
 	  // Check if the track is in the barrel
 	  if (fabs(trackz) > 300.) continue;
 
-	  double _rhsizeEta = rhsizeEta * (2.00 + 1.0 / (fracsNbr * std::min(1.,trackPt/2.)));
-	  double _rhsizePhi = rhsizePhi * (2.00 + 1.0 / (fracsNbr * std::min(1.,trackPt/2.)));
+	  double _rhsizeeta = rhsizeeta * (2.00 + 1.0 / (fracsNbr * std::min(1.,trackPt/2.)));
+	  double _rhsizephi = rhsizephi * (2.00 + 1.0 / (fracsNbr * std::min(1.,trackPt/2.)));
 	  
 	  // Check if the track and the cluster are linked
-	  if(deta < (_rhsizeEta / 2.) && dphi < (_rhsizePhi / 2.))
+	  if(deta < (_rhsizeeta / 2.) && dphi < (_rhsizephi / 2.))
 	    target2ClusterLinks_[*it].insert(*clusterIt);
 
 	  
@@ -189,10 +188,10 @@ KDTreeLinkerTrackEcal::searchLinks()
 	  double x[5];
 	  double y[5];
 	  for ( unsigned jc=0; jc<4; ++jc ) {
-	    math::XYZPoint cornerposxyz = cornersxyz[jc];
-	    x[jc] = cornerposxyz.X() + (cornerposxyz.X()-posxyz.X())
+	    auto cornerposxyz = cornersxyz[jc];
+	    x[jc] = cornerposxyz.x() + (cornerposxyz.x()-posxyz.x())
 	      * (1.00+0.50/fracsNbr /std::min(1.,trackPt/2.));
-	    y[jc] = cornerposxyz.Y() + (cornerposxyz.Y()-posxyz.Y())
+	    y[jc] = cornerposxyz.y() + (cornerposxyz.y()-posxyz.y())
 	      * (1.00+0.50/fracsNbr /std::min(1.,trackPt/2.));
 	  }
 	  
@@ -225,10 +224,10 @@ KDTreeLinkerTrackEcal::updatePFBlockEltWithLinks()
     for (BlockEltSet::iterator jt = it->second.begin();
 	 jt != it->second.end(); ++jt) {
 
-      double clusterPhi = (*jt)->clusterRef()->positionREP().Phi();
-      double clusterEta = (*jt)->clusterRef()->positionREP().Eta();
+      double clusterphi = (*jt)->clusterRef()->positionREP().phi();
+      double clustereta = (*jt)->clusterRef()->positionREP().eta();
 
-      multitracks.linkedClusters.push_back(std::make_pair(clusterPhi, clusterEta));
+      multitracks.linkedClusters.push_back(std::make_pair(clusterphi, clustereta));
     }
 
     it->first->setMultilinks(multitracks);
diff --git a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackHcal.cc b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackHcal.cc
index 66e686baafda6..6aa48b7b6fe90 100644
--- a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackHcal.cc
+++ b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackHcal.cc
@@ -74,20 +74,20 @@ KDTreeLinkerTrackHcal::buildTree()
     
     const reco::PFRecHit::REPPoint &posrep = (*it)->positionREP();
     
-    KDTreeNodeInfo rh1 (*it, posrep.Eta(), posrep.Phi());
+    KDTreeNodeInfo rh1 (*it, posrep.eta(), posrep.phi());
     eltList.push_back(rh1);
     
     // Here we solve the problem of phi circular set by duplicating some rechits
     // too close to -Pi (or to Pi) and adding (substracting) to them 2 * Pi.
     if (rh1.dim2 > (M_PI - getPhiOffset())) {
       double phi = rh1.dim2 - 2 * M_PI;
-      KDTreeNodeInfo rh2(*it, posrep.Eta(), phi); 
+      KDTreeNodeInfo rh2(*it, posrep.eta(), phi); 
       eltList.push_back(rh2);
     }
 
     if (rh1.dim2 < (M_PI * -1.0 + getPhiOffset())) {
       double phi = rh1.dim2 + 2 * M_PI;
-      KDTreeNodeInfo rh3(*it, posrep.Eta(), phi); 
+      KDTreeNodeInfo rh3(*it, posrep.eta(), phi); 
       eltList.push_back(rh3);
     }
   }
@@ -122,13 +122,13 @@ KDTreeLinkerTrackHcal::searchLinks()
     // The track didn't reach hcal
     if( ! atHCAL.isValid()) continue;
     
-    double dHEta = atHCALExit.positionREP().Eta() - atHCAL.positionREP().Eta();
-    double dHPhi = atHCALExit.positionREP().Phi() - atHCAL.positionREP().Phi(); 
-    if ( dHPhi > M_PI ) dHPhi = dHPhi - 2. * M_PI;
-    else if ( dHPhi < -M_PI ) dHPhi = dHPhi + 2. * M_PI; 
+    double dHeta = atHCALExit.positionREP().eta() - atHCAL.positionREP().eta();
+    double dHphi = atHCALExit.positionREP().phi() - atHCAL.positionREP().phi(); 
+    if ( dHphi > M_PI ) dHphi = dHphi - 2. * M_PI;
+    else if ( dHphi < -M_PI ) dHphi = dHphi + 2. * M_PI; 
     
-    double tracketa = atHCAL.positionREP().Eta() + 0.1 * dHEta;
-    double trackphi = atHCAL.positionREP().Phi() + 0.1 * dHPhi;
+    double tracketa = atHCAL.positionREP().eta() + 0.1 * dHeta;
+    double trackphi = atHCAL.positionREP().phi() + 0.1 * dHphi;
     
     if (trackphi > M_PI) trackphi -= 2 * M_PI;
     else if (trackphi < -M_PI) trackphi += 2 * M_PI;
@@ -136,29 +136,28 @@ KDTreeLinkerTrackHcal::searchLinks()
     // Estimate the maximal envelope in phi/eta that will be used to find rechit candidates.
     // Same envelope for cap et barrel rechits.
     double inflation = 1.;
-    double rangeEta = (getCristalPhiEtaMaxSize() * (1.5 + 0.5) + 0.2 * fabs(dHEta)) * inflation; 
-    double rangePhi = (getCristalPhiEtaMaxSize() * (1.5 + 0.5) + 0.2 * fabs(dHPhi)) * inflation; 
+    double rangeeta = (getCristalPhiEtaMaxSize() * (1.5 + 0.5) + 0.2 * fabs(dHeta)) * inflation; 
+    double rangephi = (getCristalPhiEtaMaxSize() * (1.5 + 0.5) + 0.2 * fabs(dHphi)) * inflation; 
 
     // We search for all candidate recHits, ie all recHits contained in the maximal size envelope.
     std::vector<KDTreeNodeInfo> recHits;
-    KDTreeBox trackBox(tracketa - rangeEta, tracketa + rangeEta, 
-		       trackphi - rangePhi, trackphi + rangePhi);
+    KDTreeBox trackBox(tracketa - rangeeta, tracketa + rangeeta, 
+		       trackphi - rangephi, trackphi + rangephi);
     tree_.search(trackBox, recHits);
     
     // Here we check all rechit candidates using the non-approximated method.
     for(std::vector<KDTreeNodeInfo>::const_iterator rhit = recHits.begin(); 
 	rhit != recHits.end(); ++rhit) {
 
-      const reco::PFRecHit::REPPoint &rhrep		   = rhit->ptr->positionREP();
-      const std::vector<reco::PFRecHit::REPPoint>& corners = rhit->ptr->getCornersREP();
-      if(corners.size() != 4) continue;
+      const auto &rhrep		   = rhit->ptr->positionREP();
+      const auto & corners = rhit->ptr->getCornersREP();
       
-      double rhsizeEta = fabs(corners[0].Eta() - corners[2].Eta());
-      double rhsizePhi = fabs(corners[0].Phi() - corners[2].Phi());
-      if ( rhsizePhi > M_PI ) rhsizePhi = 2.*M_PI - rhsizePhi;
+      double rhsizeeta = fabs(corners[0].eta() - corners[2].eta());
+      double rhsizephi = fabs(corners[0].phi() - corners[2].phi());
+      if ( rhsizephi > M_PI ) rhsizephi = 2.*M_PI - rhsizephi;
       
-      double deta = fabs(rhrep.Eta() - tracketa);
-      double dphi = fabs(rhrep.Phi() - trackphi);
+      double deta = fabs(rhrep.eta() - tracketa);
+      double dphi = fabs(rhrep.phi() - trackphi);
       if ( dphi > M_PI ) dphi = 2.*M_PI - dphi;
       
       // Find all clusters associated to given rechit
@@ -170,11 +169,11 @@ KDTreeLinkerTrackHcal::searchLinks()
 	const reco::PFClusterRef clusterref = (*clusterIt)->clusterRef();
 	int fracsNbr = clusterref->recHitFractions().size();
 	
-	double _rhsizeEta = rhsizeEta * (1.5 + 0.5 / fracsNbr) + 0.2 * fabs(dHEta);
-	double _rhsizePhi = rhsizePhi * (1.5 + 0.5 / fracsNbr) + 0.2 * fabs(dHPhi);
+	double _rhsizeeta = rhsizeeta * (1.5 + 0.5 / fracsNbr) + 0.2 * fabs(dHeta);
+	double _rhsizephi = rhsizephi * (1.5 + 0.5 / fracsNbr) + 0.2 * fabs(dHphi);
 	
 	// Check if the track and the cluster are linked
-	if(deta < (_rhsizeEta / 2.) && dphi < (_rhsizePhi / 2.))
+	if(deta < (_rhsizeeta / 2.) && dphi < (_rhsizephi / 2.))
 	  cluster2TargetLinks_[*clusterIt].insert(*it);
       }
     }
@@ -197,8 +196,8 @@ KDTreeLinkerTrackHcal::updatePFBlockEltWithLinks()
       reco::PFRecTrackRef trackref = (*jt)->trackRefPF();
       const reco::PFTrajectoryPoint& atHCAL = 
 	trackref->extrapolatedPoint(reco::PFTrajectoryPoint::HCALEntrance);
-      double tracketa = atHCAL.positionREP().Eta();
-      double trackphi = atHCAL.positionREP().Phi();
+      double tracketa = atHCAL.positionREP().eta();
+      double trackphi = atHCAL.positionREP().phi();
       
       multitracks.linkedClusters.push_back(std::make_pair(trackphi, tracketa));
     }
diff --git a/RecoParticleFlow/PFProducer/src/PFAlgo.cc b/RecoParticleFlow/PFProducer/src/PFAlgo.cc
index 6cbc3b38cf3e7..f5aa3d539e6a0 100644
--- a/RecoParticleFlow/PFProducer/src/PFAlgo.cc
+++ b/RecoParticleFlow/PFProducer/src/PFAlgo.cc
@@ -3592,7 +3592,7 @@ PFAlgo::checkCleaning( const reco::PFRecHitCollection& cleanedHits ) {
     // Loop on the candidates
     for(unsigned i=0; i<cleanedHits.size(); ++i) {
       const PFRecHit& hit = cleanedHits[i];
-      double length = std::sqrt(hit.position().Mag2()); 
+      double length = std::sqrt(hit.position().mag2()); 
       double px = hit.energy() * hit.position().x()/length;
       double py = hit.energy() * hit.position().y()/length;
       double pt = std::sqrt(px*px + py*py);

From bec604b01c8123e5b7a0c21a69894dffb318c1fe Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sat, 16 Apr 2016 15:23:15 +0200
Subject: [PATCH 05/16] Fireworks compiles, need to fix corners

---
 Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc b/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
index 13075dfaf0849..08b6349e50624 100644
--- a/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
+++ b/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
@@ -382,7 +382,8 @@ void FWPFCandidateDetailView::addHits( const std::vector<reco::PFRecHit> *hits)
 
    for (std::vector<reco::PFRecHit>::const_iterator it = hits->begin(); it != hits->end(); ++it)
    {
-      const std::vector< math::XYZPoint >& corners = it->getCornersXYZ();
+      // FIXME reading from reco file we will need to get the corners from the geometry...
+      const auto & corners = it->getCornersXYZ();
       if (!isPntInRng(corners[0].eta(), corners[0].phi()))
          continue;
      

From 7f29f130b4e434531d0982ff825806a1c361ffb7 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sat, 16 Apr 2016 16:57:26 +0200
Subject: [PATCH 06/16] add PF geometry for HF

---
 .../plugins/FWPFCandidateDetailView.cc        |  3 +--
 Geometry/CaloGeometry/interface/IdealZPrism.h | 26 ++++++++++++++----
 Geometry/CaloGeometry/src/IdealZPrism.cc      | 27 ++++++++++++++++---
 Geometry/HcalEventSetup/src/moduleDB.cc       |  2 +-
 Geometry/HcalTowerAlgo/src/HcalDDDGeometry.cc |  2 +-
 .../src/HcalDDDGeometryLoader.cc              |  4 +--
 .../src/HcalFlexiHardcodeGeometryLoader.cc    |  2 +-
 Geometry/HcalTowerAlgo/src/HcalGeometry.cc    |  2 +-
 .../src/HcalHardcodeGeometryLoader.cc         |  2 +-
 .../interface/PFHFRecHitCreator.h             |  6 +++++
 .../plugins/PFCTRecHitProducer.cc             |  7 ++++-
 11 files changed, 64 insertions(+), 19 deletions(-)

diff --git a/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc b/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
index 08b6349e50624..13075dfaf0849 100644
--- a/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
+++ b/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
@@ -382,8 +382,7 @@ void FWPFCandidateDetailView::addHits( const std::vector<reco::PFRecHit> *hits)
 
    for (std::vector<reco::PFRecHit>::const_iterator it = hits->begin(); it != hits->end(); ++it)
    {
-      // FIXME reading from reco file we will need to get the corners from the geometry...
-      const auto & corners = it->getCornersXYZ();
+      const std::vector< math::XYZPoint >& corners = it->getCornersXYZ();
       if (!isPntInRng(corners[0].eta(), corners[0].phi()))
          continue;
      
diff --git a/Geometry/CaloGeometry/interface/IdealZPrism.h b/Geometry/CaloGeometry/interface/IdealZPrism.h
index 1ec922691f049..f3c3b2099eded 100644
--- a/Geometry/CaloGeometry/interface/IdealZPrism.h
+++ b/Geometry/CaloGeometry/interface/IdealZPrism.h
@@ -2,6 +2,7 @@
 #define GEOMETRY_CALOGEOMETRY_IDEALZPRISM_H 1
 
 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
+#include <memory>
 
 /** \class IdealZPrism
     
@@ -22,10 +23,12 @@ parameters are eta and phi HALF-widths and the tower z thickness.
 
 \author J. Mans - Minnesota
 */
-class IdealZPrism : public CaloCellGeometry 
+class IdealZPrism final : public CaloCellGeometry 
 {
-   public:
-      
+ public:
+
+  enum DEPTH {None, EM, HADR};
+  
       typedef CaloCellGeometry::CCGFloat CCGFloat ;
       typedef CaloCellGeometry::Pt3D     Pt3D     ;
       typedef CaloCellGeometry::Pt3DVec  Pt3DVec  ;
@@ -38,7 +41,8 @@ class IdealZPrism : public CaloCellGeometry
       
       IdealZPrism( const GlobalPoint& faceCenter , 
 		   CornersMgr*        mgr        ,
-		   const CCGFloat*    parm         ) ;
+		   const CCGFloat*    parm       ,
+			  IdealZPrism::DEPTH depth) ;
       
       virtual ~IdealZPrism() ;
       
@@ -55,7 +59,13 @@ class IdealZPrism : public CaloCellGeometry
       virtual void vocalCorners( Pt3DVec&        vec ,
 				 const CCGFloat* pv  ,
 				 Pt3D&           ref   ) const override;
-      
+
+
+  
+  
+      // corrected geom for PF
+      IdealZPrism const *  forPF() const  { return  m_geoForPF.get();}
+  
    private:
 
       virtual void initCorners(CornersVec& ) override;
@@ -71,6 +81,12 @@ class IdealZPrism : public CaloCellGeometry
       static GlobalPoint etaPhiZ( float eta , 
 				  float phi ,
 				  float z    ) ;
+
+
+private:
+      // corrected geom for PF
+      std::unique_ptr<IdealZPrism> m_geoForPF;
+
 };
 
 std::ostream& operator<<( std::ostream& s , const IdealZPrism& cell ) ;
diff --git a/Geometry/CaloGeometry/src/IdealZPrism.cc b/Geometry/CaloGeometry/src/IdealZPrism.cc
index 6cce524d477d6..fbd5303ea1638 100644
--- a/Geometry/CaloGeometry/src/IdealZPrism.cc
+++ b/Geometry/CaloGeometry/src/IdealZPrism.cc
@@ -9,23 +9,42 @@ IdealZPrism::IdealZPrism()
   : CaloCellGeometry()
 {}
 
+namespace {
+
+  // magic numbers determined by ParticleFlow
+  constexpr float  EMDepthCorrection  = 22.;
+  constexpr float  HADDepthCorrection = 25.;
+
+  GlobalPoint correct(GlobalPoint const & ori, IdealZPrism::DEPTH depth) {
+    if (depth==IdealZPrism::None) return ori;
+    float zcorr = depth==IdealZPrism::EM ?EMDepthCorrection :  HADDepthCorrection;
+    if (ori.z()<0) zcorr = -zcorr;
+    return ori + GlobalVector(0.,0.,zcorr);
+  }
+}
+
 IdealZPrism::IdealZPrism( const IdealZPrism& idzp ) 
   : CaloCellGeometry( idzp )
 {
-  *this = idzp ;
+  if (idzp.forPF()) m_geoForPF.reset(new IdealZPrism(*idzp.forPF()));
 }
 
 IdealZPrism& 
 IdealZPrism::operator=( const IdealZPrism& idzp ) 
 {
-  if( &idzp != this ) CaloCellGeometry::operator=( idzp ) ;
+  if( &idzp != this ) {
+     CaloCellGeometry::operator=( idzp ) ;
+     if (idzp.forPF()) m_geoForPF.reset(new IdealZPrism(*idzp.forPF()));
+  }
   return *this ;
 }
 
 IdealZPrism::IdealZPrism( const GlobalPoint& faceCenter , 
 			  CornersMgr*  mgr              ,
-			  const CCGFloat*    parm         )
-  : CaloCellGeometry ( faceCenter, mgr, parm )   
+			  const CCGFloat*    parm       ,
+			  IdealZPrism::DEPTH depth)
+  : CaloCellGeometry ( faceCenter, mgr, parm ),
+    m_geoForPF(depth==None ? nullptr : new IdealZPrism(correct(faceCenter,depth), mgr, parm, None ))
 {initSpan();}
 
 IdealZPrism::~IdealZPrism() 
diff --git a/Geometry/HcalEventSetup/src/moduleDB.cc b/Geometry/HcalEventSetup/src/moduleDB.cc
index 79ba4bbdd3af8..b0f87182ed12a 100644
--- a/Geometry/HcalEventSetup/src/moduleDB.cc
+++ b/Geometry/HcalEventSetup/src/moduleDB.cc
@@ -68,7 +68,7 @@ CaloGeometryDBEP<HcalGeometry, CaloGeometryDBReader>::produceAligned( const type
    
   ptr->fillDefaultNamedParameters() ;
   
-  ptr->allocateCorners( hcalTopology->ncells() );
+  ptr->allocateCorners( hcalTopology->ncells()+hcalTopology->getHFSize() );
 
   ptr->allocatePar(    dvec.size() ,
 		       HcalGeometry::k_NumberOfParametersPerShape ) ;
diff --git a/Geometry/HcalTowerAlgo/src/HcalDDDGeometry.cc b/Geometry/HcalTowerAlgo/src/HcalDDDGeometry.cc
index 342b825400ad4..506c11add14cf 100644
--- a/Geometry/HcalTowerAlgo/src/HcalDDDGeometry.cc
+++ b/Geometry/HcalTowerAlgo/src/HcalDDDGeometry.cc
@@ -192,7 +192,7 @@ HcalDDDGeometry::newCell( const GlobalPoint& f1 ,
 				   - m_hbCellVec.size() 
 				   - m_heCellVec.size() 
 				   - m_hoCellVec.size() ) ;
-	m_hfCellVec[ index ] = IdealZPrism( f1, cornersMgr(), parm ) ;
+	m_hfCellVec[ index ] = IdealZPrism( f1, cornersMgr(), parm, hId.depth()==1 ? IdealZPrism::EM : IdealZPrism::HADR ) ;
       }
     }
   }
diff --git a/Geometry/HcalTowerAlgo/src/HcalDDDGeometryLoader.cc b/Geometry/HcalTowerAlgo/src/HcalDDDGeometryLoader.cc
index 03bc98fd0bea5..7fecb85a51744 100644
--- a/Geometry/HcalTowerAlgo/src/HcalDDDGeometryLoader.cc
+++ b/Geometry/HcalTowerAlgo/src/HcalDDDGeometryLoader.cc
@@ -31,7 +31,7 @@ HcalDDDGeometryLoader::load(const HcalTopology& topo, DetId::Detector det, int s
   if ( geom->cornersMgr() == 0 ) {
      const unsigned int count (hcalConstants->numberOfCells(HcalBarrel ) +
 			       hcalConstants->numberOfCells(HcalEndcap ) +
-			       hcalConstants->numberOfCells(HcalForward) +
+			       2*hcalConstants->numberOfCells(HcalForward) +
 			       hcalConstants->numberOfCells(HcalOuter  ) );
      geom->allocateCorners( count ) ;
   }
@@ -53,7 +53,7 @@ HcalDDDGeometryLoader::load(const HcalTopology& topo) {
   if( geom->cornersMgr() == 0 ) {
     const unsigned int count (hcalConstants->numberOfCells(HcalBarrel ) +
 			      hcalConstants->numberOfCells(HcalEndcap ) +
-			      hcalConstants->numberOfCells(HcalForward) +
+			      2*hcalConstants->numberOfCells(HcalForward) +
 			      hcalConstants->numberOfCells(HcalOuter  ) );
     geom->allocateCorners( count ) ;
   }
diff --git a/Geometry/HcalTowerAlgo/src/HcalFlexiHardcodeGeometryLoader.cc b/Geometry/HcalTowerAlgo/src/HcalFlexiHardcodeGeometryLoader.cc
index 970dd215a8e88..70cb7a4590230 100644
--- a/Geometry/HcalTowerAlgo/src/HcalFlexiHardcodeGeometryLoader.cc
+++ b/Geometry/HcalTowerAlgo/src/HcalFlexiHardcodeGeometryLoader.cc
@@ -23,7 +23,7 @@ HcalFlexiHardcodeGeometryLoader::HcalFlexiHardcodeGeometryLoader(const edm::Para
 
 CaloSubdetectorGeometry* HcalFlexiHardcodeGeometryLoader::load(const HcalTopology& fTopology, const HcalDDDRecConstants& hcons) {
   CaloSubdetectorGeometry* hcalGeometry = new HcalGeometry (fTopology);
-  if( 0 == hcalGeometry->cornersMgr() ) hcalGeometry->allocateCorners ( fTopology.ncells() );
+  if( 0 == hcalGeometry->cornersMgr() ) hcalGeometry->allocateCorners ( fTopology.ncells()+fTopology.getHFSize() );
   if( 0 == hcalGeometry->parMgr() ) hcalGeometry->allocatePar (hcalGeometry->numberOfShapes(),
 							       HcalGeometry::k_NumberOfParametersPerShape ) ;
 #ifdef DebugLog
diff --git a/Geometry/HcalTowerAlgo/src/HcalGeometry.cc b/Geometry/HcalTowerAlgo/src/HcalGeometry.cc
index 4d9a59f2f835e..71b35995f2927 100644
--- a/Geometry/HcalTowerAlgo/src/HcalGeometry.cc
+++ b/Geometry/HcalTowerAlgo/src/HcalGeometry.cc
@@ -366,7 +366,7 @@ void HcalGeometry::newCell(const GlobalPoint& f1 ,
 			       - m_hbCellVec.size()
 			       - m_heCellVec.size()
 			       - m_hoCellVec.size() ) ;
-    m_hfCellVec[ index ] = IdealZPrism( f1, cornersMgr(), parm ) ;
+    m_hfCellVec[ index ] = IdealZPrism( f1, cornersMgr(), parm,  hid.depth()==1 ? IdealZPrism::EM : IdealZPrism::HADR ) ;
   }
 
   addValidID( detId ) ;
diff --git a/Geometry/HcalTowerAlgo/src/HcalHardcodeGeometryLoader.cc b/Geometry/HcalTowerAlgo/src/HcalHardcodeGeometryLoader.cc
index 7062baacf0505..1fe5b729f1823 100644
--- a/Geometry/HcalTowerAlgo/src/HcalHardcodeGeometryLoader.cc
+++ b/Geometry/HcalTowerAlgo/src/HcalHardcodeGeometryLoader.cc
@@ -35,7 +35,7 @@ CaloSubdetectorGeometry* HcalHardcodeGeometryLoader::load(const HcalTopology& fT
 #endif
   }
   CaloSubdetectorGeometry* hcalGeometry = new HcalGeometry (fTopology);
-  if( 0 == hcalGeometry->cornersMgr() ) hcalGeometry->allocateCorners ( fTopology.ncells() );
+  if( 0 == hcalGeometry->cornersMgr() ) hcalGeometry->allocateCorners ( fTopology.ncells()+fTopology.getHFSize() );
   if( 0 == hcalGeometry->parMgr() ) hcalGeometry->allocatePar (hcalGeometry->numberOfShapes(),
 							       HcalGeometry::k_NumberOfParametersPerShape ) ;
   if (fTopology.mode() == HcalTopologyMode::H2) {  // TB geometry
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
index ae63dd227c16f..fdd8318f9ba23 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
@@ -14,6 +14,8 @@
 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
 #include "Geometry/Records/interface/CaloGeometryRecord.h"
 #include "Geometry/CaloTopology/interface/HcalTopology.h"
+#include "Geometry/CaloGeometry/interface/IdealZPrism.h"
+
 
 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
 
@@ -62,7 +64,11 @@ class PFHFRecHitCreator final :  public  PFRecHitCreatorBase {
 	auto time = erh.time();
 
 	const CaloCellGeometry * thisCell= hcalGeo->getGeometry(detid);
+	auto zp = dynamic_cast<IdealZPrism const*>(thisCell);
+	assert(zp);
+	thisCell = zp->forPF();
 
+	
 	// find rechit geometry
 	if(!thisCell) {
 	  edm::LogError("PFHFRecHitCreator")
diff --git a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
index 83bf89b7155e7..315012ba8db8f 100644
--- a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
+++ b/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
@@ -23,6 +23,7 @@
 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
+#include "Geometry/CaloGeometry/interface/IdealZPrism.h"
 #include "Geometry/Records/interface/CaloGeometryRecord.h"
 
 #include "Geometry/CaloTopology/interface/HcalTopology.h"
@@ -861,9 +862,13 @@ PFCTRecHitProducer::createHcalRecHit( const DetId& detid,
   switch ( layer ) { 
   case PFLayer::HF_EM:
     // depth_correction = position.z() > 0. ? EM_Depth_ : -EM_Depth_;
-    break;
   case PFLayer::HF_HAD:
+  {
     //  depth_correction = position.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
+    auto zp = dynamic_cast<IdealZPrism const*>(thisCell);
+    assert(zp);
+    thisCell = zp->forPF();
+  }
     break;
   default:
     break;

From 588f02b6df2c041380ee3d05019ac700142627cb Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sun, 17 Apr 2016 09:47:34 +0200
Subject: [PATCH 07/16] fix order of corners

---
 .../PFClusterTools/src/LinkByRecHit.cc            | 15 +++++++--------
 .../plugins/kdtrees/KDTreeLinkerPSEcal.cc         |  4 ++--
 .../plugins/kdtrees/KDTreeLinkerTrackEcal.cc      |  8 ++++----
 .../plugins/kdtrees/KDTreeLinkerTrackHcal.cc      |  4 ++--
 4 files changed, 15 insertions(+), 16 deletions(-)

diff --git a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
index 2a08acc5647ff..39805a8592b6b 100644
--- a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
+++ b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
@@ -219,8 +219,7 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
       = rechit_cluster.positionREP();
     
     //getting rechit corners
-    const auto & 
-      cornersxyz = rechit_cluster.getCornersXYZ();
+    const auto & cornersxyz = rechit_cluster.getCornersXYZ();
     const auto & corners = rechit_cluster.getCornersREP();
 
     
@@ -231,9 +230,9 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
       // blown up by 50% (HCAL) to 100% (ECAL) to include cracks & gaps
       // also blown up to account for multiple scattering at low pt.
       double rhsizeEta 
-	= std::abs(corners[0].eta() - corners[2].eta());
+	= std::abs(corners[3].eta() - corners[1].eta());
       double rhsizePhi 
-	= std::abs(corners[0].phi() - corners[2].phi());
+	= std::abs(corners[3].phi() - corners[1].phi());
       if ( rhsizePhi > M_PI ) rhsizePhi = 2.*M_PI - rhsizePhi;
       if ( hcal ) { 
 	const double mult = horesolscale * (1.50 + 0.5/fracs.size());
@@ -308,8 +307,8 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
       for ( unsigned jc=0; jc<4; ++jc ) {
 	const auto & cornerposxyz = cornersxyz[jc];
 	const double mult = (1.00+0.50/(fracs.size()*std::min(1.,0.5*trackPt)));
-	x[jc] = cornerposxyz.x() + (cornerposxyz.y()-posxyz.x()) * mult;
-	y[jc] = cornerposxyz.y() + (cornerposxyz.y()-posxyz.y()) * mult;
+	x[3-jc] = cornerposxyz.x() + (cornerposxyz.y()-posxyz.x()) * mult;
+	y[3-jc] = cornerposxyz.y() + (cornerposxyz.y()-posxyz.y()) * mult;
 	
 #ifdef PFLOW_DEBUG
 	if( debug ){
@@ -452,8 +451,8 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
       // corner position projected onto the preshower
       auto cornerpos = corners[jc].basicVector() * zPS/zECAL;
       // Inflate the size by the size of the PS strips, and by 5% to include ECAL cracks.
-      x[jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/fabs((cornerpos.x()-posxyz.x()))*0.5*deltaX);
-      y[jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/fabs((cornerpos.y()-posxyz.y()))*0.5*deltaY);
+      x[3-jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/fabs((cornerpos.x()-posxyz.x()))*0.5*deltaX);
+      y[3-jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/fabs((cornerpos.y()-posxyz.y()))*0.5*deltaY);
       
 #ifdef PFLOW_DEBUG
       if( debug ){
diff --git a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc
index 22423d051b274..da333abe3b63e 100644
--- a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc
+++ b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc
@@ -188,8 +188,8 @@ KDTreeLinkerPSEcal::searchLinks()
 	double y[5];
 	for ( unsigned jc=0; jc<4; ++jc ) {
 	  auto cornerpos = corners[jc].basicVector() * zPS / clusterz;
-	  x[jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/fabs((cornerpos.x()-posxyz.x()))*deltaX/2.);
-	  y[jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/fabs((cornerpos.y()-posxyz.y()))*deltaY/2.);
+	  x[3-jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/fabs((cornerpos.x()-posxyz.x()))*deltaX/2.);
+	  y[3-jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/fabs((cornerpos.y()-posxyz.y()))*deltaY/2.);
 	}
 
 	x[4] = x[0];
diff --git a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc
index 3dfaaeba31155..ac84786d53ce1 100644
--- a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc
+++ b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc
@@ -149,8 +149,8 @@ KDTreeLinkerTrackEcal::searchLinks()
       const auto &rhrep		   = rhit->ptr->positionREP();
       const auto & corners = rhit->ptr->getCornersREP();
       
-      double rhsizeeta = fabs(corners[0].eta() - corners[2].eta());
-      double rhsizephi = fabs(corners[0].phi() - corners[2].phi());
+      double rhsizeeta = fabs(corners[3].eta() - corners[1].eta());
+      double rhsizephi = fabs(corners[3].phi() - corners[1].phi());
       if ( rhsizephi > M_PI ) rhsizephi = 2.*M_PI - rhsizephi;
       
       double deta = fabs(rhrep.eta() - tracketa);
@@ -189,9 +189,9 @@ KDTreeLinkerTrackEcal::searchLinks()
 	  double y[5];
 	  for ( unsigned jc=0; jc<4; ++jc ) {
 	    auto cornerposxyz = cornersxyz[jc];
-	    x[jc] = cornerposxyz.x() + (cornerposxyz.x()-posxyz.x())
+	    x[3-jc] = cornerposxyz.x() + (cornerposxyz.x()-posxyz.x())
 	      * (1.00+0.50/fracsNbr /std::min(1.,trackPt/2.));
-	    y[jc] = cornerposxyz.y() + (cornerposxyz.y()-posxyz.y())
+	    y[3-jc] = cornerposxyz.y() + (cornerposxyz.y()-posxyz.y())
 	      * (1.00+0.50/fracsNbr /std::min(1.,trackPt/2.));
 	  }
 	  
diff --git a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackHcal.cc b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackHcal.cc
index 6aa48b7b6fe90..394f41c3be0ad 100644
--- a/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackHcal.cc
+++ b/RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackHcal.cc
@@ -152,8 +152,8 @@ KDTreeLinkerTrackHcal::searchLinks()
       const auto &rhrep		   = rhit->ptr->positionREP();
       const auto & corners = rhit->ptr->getCornersREP();
       
-      double rhsizeeta = fabs(corners[0].eta() - corners[2].eta());
-      double rhsizephi = fabs(corners[0].phi() - corners[2].phi());
+      double rhsizeeta = fabs(corners[3].eta() - corners[1].eta());
+      double rhsizephi = fabs(corners[3].phi() - corners[1].phi());
       if ( rhsizephi > M_PI ) rhsizephi = 2.*M_PI - rhsizephi;
       
       double deta = fabs(rhrep.eta() - tracketa);

From 12452df42729984db2447f67f236ca7136e7a429 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sun, 17 Apr 2016 11:20:19 +0200
Subject: [PATCH 08/16] fix mistype, correct axis definition

---
 RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc     | 2 +-
 RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc | 6 ++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
index 39805a8592b6b..6f319d1977509 100644
--- a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
+++ b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
@@ -307,7 +307,7 @@ LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
       for ( unsigned jc=0; jc<4; ++jc ) {
 	const auto & cornerposxyz = cornersxyz[jc];
 	const double mult = (1.00+0.50/(fracs.size()*std::min(1.,0.5*trackPt)));
-	x[3-jc] = cornerposxyz.x() + (cornerposxyz.y()-posxyz.x()) * mult;
+	x[3-jc] = cornerposxyz.x() + (cornerposxyz.x()-posxyz.x()) * mult;
 	y[3-jc] = cornerposxyz.y() + (cornerposxyz.y()-posxyz.y()) * mult;
 	
 #ifdef PFLOW_DEBUG
diff --git a/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc b/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc
index d2e83672a7f34..077b1d9a60407 100644
--- a/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc
+++ b/RecoParticleFlow/PFClusterTools/src/PFPhotonClusters.cc
@@ -33,8 +33,10 @@ void PFPhotonClusters::SetSeed(){
     auto frac=it->fraction();
     float E= RefPFRecHit->energy()* frac;
     if(E>PFSeedE){
-      PFSeedE=E;  
-      axis=  dynamic_cast<TruncatedPyramid const &>(RefPFRecHit->caloCell()).axis();
+      PFSeedE=E;
+      // FIXME will optimize later...
+      auto const & pyr = dynamic_cast<TruncatedPyramid const &>(RefPFRecHit->caloCell());
+      axis = pyr.getPosition(1) - pyr.getPosition(0);
       position=RefPFRecHit->position();
       idseed = RefPFRecHit->detId();
     }

From 1c23df40f54038109d2f2ab45fff3bd11be4e013 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sun, 17 Apr 2016 14:12:56 +0200
Subject: [PATCH 09/16] mugrate to float, protect FireWork

---
 .../ParticleFlowReco/interface/PFRecHit.h     | 14 ++++-----
 .../ParticleFlowReco/src/classes_def_2.xml    |  2 +-
 .../plugins/FWPFCandidateDetailView.cc        |  6 +++-
 .../interface/PFHFRecHitCreator.h             | 29 ++++---------------
 .../interface/PFPSRecHitCreator.h             |  5 ++--
 5 files changed, 20 insertions(+), 36 deletions(-)

diff --git a/DataFormats/ParticleFlowReco/interface/PFRecHit.h b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
index fcf53b8caaa9b..366d9a0af3f52 100644
--- a/DataFormats/ParticleFlowReco/interface/PFRecHit.h
+++ b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
@@ -43,9 +43,9 @@ namespace reco {
     /// default constructor. Sets energy and position to zero
     PFRecHit(){}
 
-    PFRecHit(CaloCellGeometry const * caloCell, unsigned detId,
+    PFRecHit(CaloCellGeometry const * caloCell, unsigned int detId,
              PFLayer::Layer layer,
-             double energy) :
+             float energy) :
         caloCell_(caloCell),  detId_(detId),
         layer_(layer), energy_(energy){}
 
@@ -61,7 +61,7 @@ namespace reco {
     /// destructor
     ~PFRecHit()=default;
 
-    void setEnergy( double energy) { energy_ = energy; }
+    void setEnergy( float energy) { energy_ = energy; }
 
 
     void addNeighbour(short x,short y, short z,const PFRecHitRef&);
@@ -99,11 +99,11 @@ namespace reco {
     PFLayer::Layer layer() const { return layer_; }
 
     /// rechit energy
-    double energy() const { return energy_; }
+    float energy() const { return energy_; }
 
 
     /// timing for cleaned hits
-    double time() const { return time_; }
+    float time() const { return time_; }
 
     /// depth for segemntation
     int  depth() const { return depth_; }
@@ -149,10 +149,10 @@ namespace reco {
     PFLayer::Layer      layer_=PFLayer::NONE;
 
     /// rechit energy 
-    double              energy_=0;
+    float              energy_=0;
 
     /// time
-    double              time_=-1;
+    float              time_=-1;
 
     /// depth
     int      depth_=0;
diff --git a/DataFormats/ParticleFlowReco/src/classes_def_2.xml b/DataFormats/ParticleFlowReco/src/classes_def_2.xml
index fd03a39045a4d..e325ca5f965c5 100644
--- a/DataFormats/ParticleFlowReco/src/classes_def_2.xml
+++ b/DataFormats/ParticleFlowReco/src/classes_def_2.xml
@@ -37,7 +37,7 @@
 
 
   <class name="reco::PFRecHit" ClassVersion="18">
-   <version ClassVersion="18" checksum="  1826948750"/>
+   <version ClassVersion="18" checksum=" 3589980878"/>
    <version ClassVersion="17" checksum="2606786437"/>
    <version ClassVersion="10" checksum="1807687081"/>
    <version ClassVersion="11" checksum="3727193351"/>
diff --git a/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc b/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
index 13075dfaf0849..6cee89cdd160c 100644
--- a/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
+++ b/Fireworks/ParticleFlow/plugins/FWPFCandidateDetailView.cc
@@ -261,6 +261,8 @@ void FWPFCandidateDetailView::voteMaxEtEVal( const std::vector<reco::PFRecHit> *
 {
    if (!hits) return;
 
+   // FIXME: require access to geometry while reading from reco file
+   if ( (!hits->empty()) && hits->front().hasCaloCell())
    for (std::vector<reco::PFRecHit>::const_iterator it = hits->begin(); it != hits->end(); ++it)
    {
       TEveVector centre(it->position().x(), it->position().y(), it->position().z());
@@ -380,9 +382,11 @@ void FWPFCandidateDetailView::addHits( const std::vector<reco::PFRecHit> *hits)
    m_eventList->AddElement(ps);
    ps->SetMainColor(kOrange);
 
+      // FIXME, requires access to geometry
+   if ( (!hits->empty()) && hits->front().hasCaloCell()) 
    for (std::vector<reco::PFRecHit>::const_iterator it = hits->begin(); it != hits->end(); ++it)
    {
-      const std::vector< math::XYZPoint >& corners = it->getCornersXYZ();
+      const auto & corners = it->getCornersXYZ();
       if (!isPntInRng(corners[0].eta(), corners[0].phi()))
          continue;
      
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
index fdd8318f9ba23..2cb287ed09be8 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFHFRecHitCreator.h
@@ -67,7 +67,6 @@ class PFHFRecHitCreator final :  public  PFRecHitCreatorBase {
 	auto zp = dynamic_cast<IdealZPrism const*>(thisCell);
 	assert(zp);
 	thisCell = zp->forPF();
-
 	
 	// find rechit geometry
 	if(!thisCell) {
@@ -77,28 +76,10 @@ class PFHFRecHitCreator final :  public  PFRecHitCreatorBase {
 	  continue;
 	}
 
-	// auto const point =  thisCell->getPosition();
-
-
-	PFLayer::Layer layer;
-	//	double depth_correction;
-	if (depth==1) {
-	  layer = PFLayer::HF_EM;
-          // depth_correction = point.z() > 0. ? EM_Depth_ : -EM_Depth_;
-	}
-	else {
-	  layer = PFLayer::HF_HAD;
-	  // depth_correction = point.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
-	}
-
-        /*
-	position.SetCoordinates ( point.x(),
-				  point.y(),
-				  point.z()+depth_correction );
-	*/
+	PFLayer::Layer layer  =  depth==1 ? PFLayer::HF_EM : PFLayer::HF_HAD;
+       
 
-	reco::PFRecHit rh(thisCell, detid.rawId(),layer,
-			   energy);
+	reco::PFRecHit rh(thisCell, detid.rawId(),layer,energy);
 	rh.setTime(time); 
 	rh.setDepth(depth);
 
@@ -114,10 +95,10 @@ class PFHFRecHitCreator final :  public  PFRecHitCreatorBase {
 	}
 	  
 	if(keep) {
-	  tmpOut.push_back(rh);
+	  tmpOut.push_back(std::move(rh));
 	}
 	else if (rcleaned) 
-	  cleaned->push_back(rh);
+	  cleaned->push_back(std::move(rh));
       }
       //Sort by DetID the collection
       DetIDSorter sorter;
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFPSRecHitCreator.h b/RecoParticleFlow/PFClusterProducer/interface/PFPSRecHitCreator.h
index ad990e3582e93..19724dcae95da 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFPSRecHitCreator.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFPSRecHitCreator.h
@@ -48,7 +48,7 @@ class PFPSRecHitCreator final :  public  PFRecHitCreatorBase {
       iEvent.getByToken(recHitToken_,recHitHandle);
       for( const auto& erh : *recHitHandle ) {      
 	ESDetId detid(erh.detid());
-	double energy = erh.energy();
+	auto energy = erh.energy();
 
 	PFLayer::Layer layer = PFLayer::NONE;
 	
@@ -77,8 +77,7 @@ class PFPSRecHitCreator final :  public  PFRecHitCreatorBase {
 	  continue;
 	}
 
-        out->emplace_back(thisCell, detid.rawId(),layer,
-			   energy);
+        out->emplace_back(thisCell, detid.rawId(),layer,energy);
         auto & rh = out->back();
 	rh.setDepth(detid.plane());
 	rh.setTime(erh.time());

From ceec02921adc1f62d3e15d49b0069770b228afd6 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sun, 17 Apr 2016 16:58:35 +0200
Subject: [PATCH 10/16] use a single vector for neighbours

---
 .../ParticleFlowReco/interface/PFRecHit.h     | 36 +++++++++++-----
 DataFormats/ParticleFlowReco/src/PFRecHit.cc  | 24 +++++++----
 .../ParticleFlowReco/src/classes_def_2.xml    |  8 +++-
 .../interface/PFRecHitCaloNavigatorWithTime.h |  2 +-
 .../interface/PFRecHitNavigatorBase.h         |  2 +-
 .../src/Basic2DGenericPFlowPositionCalc.cc    | 12 +++---
 .../src/Basic2DGenericTopoClusterizer.cc      | 42 ++++++++++---------
 .../src/Basic2DGenericTopoClusterizer.h       |  2 +-
 .../src/LocalMaximumSeedFinder.cc             | 22 +++++-----
 .../PFClusterProducer/src/RBXAndHPDCleaner.cc | 11 ++---
 .../src/SpikeAndDoubleSpikeCleaner.cc         | 39 ++++++++---------
 11 files changed, 114 insertions(+), 86 deletions(-)

diff --git a/DataFormats/ParticleFlowReco/interface/PFRecHit.h b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
index 366d9a0af3f52..694e9157e56de 100644
--- a/DataFormats/ParticleFlowReco/interface/PFRecHit.h
+++ b/DataFormats/ParticleFlowReco/interface/PFRecHit.h
@@ -36,6 +36,16 @@ namespace reco {
     using RepCorners = CaloCellGeometry::RepCorners;
     using REPPointVector = RepCorners;
     using CornersVec = CaloCellGeometry::CornersVec;
+
+    struct Neighbours {
+      using Pointer = unsigned int const *;
+      Neighbours(){}
+      Neighbours(Pointer ib, unsigned int n) : b(ib), e(ib+n){}
+      Pointer b, e;
+      Pointer begin() const {return b;}
+      Pointer end() const {return e;}
+      unsigned int size() const { return e-b;}
+    };
     
     enum {
       NONE=0
@@ -64,23 +74,25 @@ namespace reco {
     void setEnergy( float energy) { energy_ = energy; }
 
 
-    void addNeighbour(short x,short y, short z,const PFRecHitRef&);
-    const PFRecHitRef getNeighbour(short x,short y, short z);
+    void addNeighbour(short x,short y, short z, unsigned int);
+    unsigned int getNeighbour(short x,short y, short z) const;
     void setTime( double time) { time_ = time; }
     void setDepth( int depth) { depth_ = depth; }
     void clearNeighbours() {
       neighbours_.clear();
+      neighbourInfos_.clear();
+      neighbours4_ = neighbours8_ = 0;
     }
 
-    const PFRecHitRefVector& neighbours4() const {
-      return neighbours4_;
+    Neighbours neighbours4() const {
+      return buildNeighbours(neighbours4_);
     }
-    const PFRecHitRefVector& neighbours8() const {
-      return neighbours8_;
+    Neighbours neighbours8() const {
+	return buildNeighbours(neighbours8_);
     }
 
-    const PFRecHitRefVector& neighbours() const {
-      return neighbours_;
+    Neighbours neighbours() const {
+      return buildNeighbours(neighbours_.size());
     }
 
     const std::vector<unsigned short>& neighbourInfos() {
@@ -139,6 +151,8 @@ namespace reco {
  
   private:
 
+    Neighbours buildNeighbours(unsigned int n) const { return  Neighbours(&neighbours_.front(),n);}
+    
     /// cell geometry
     CaloCellGeometry const * caloCell_=nullptr;
  
@@ -159,12 +173,12 @@ namespace reco {
 
   
     /// indices to existing neighbours (1 common side)
-    PFRecHitRefVector   neighbours_;
+    std::vector< unsigned int > neighbours_;
     std::vector< unsigned short >   neighbourInfos_;
 
     //Caching the neighbours4/8 per request of Lindsey
-    PFRecHitRefVector   neighbours4_;
-    PFRecHitRefVector   neighbours8_;
+    unsigned int neighbours4_ = 0;
+    unsigned int neighbours8_ = 0;
   };
 
 }
diff --git a/DataFormats/ParticleFlowReco/src/PFRecHit.cc b/DataFormats/ParticleFlowReco/src/PFRecHit.cc
index f3a80948548a6..ce1c3992d1310 100644
--- a/DataFormats/ParticleFlowReco/src/PFRecHit.cc
+++ b/DataFormats/ParticleFlowReco/src/PFRecHit.cc
@@ -1,8 +1,8 @@
 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
-
+#include<limits>
 namespace reco {
 
-void PFRecHit::addNeighbour(short x,short y,short z,const PFRecHitRef& ref) {
+void PFRecHit::addNeighbour(short x,short y,short z, unsigned int ref) {
   //bitmask interface  to accomodate more advanced naighbour finding [i.e in z as well]
   //bit 0 side for eta [0 for <=0 , 1 for >0]
   //bits 1,2,3 : abs(eta) wrt the center
@@ -28,19 +28,25 @@ void PFRecHit::addNeighbour(short x,short y,short z,const PFRecHitRef& ref) {
     bitmask = bitmask | (1<<8) ;
   bitmask = bitmask | (absz << 9);
   
-  neighbours_.push_back(ref);
-  neighbourInfos_.push_back(bitmask);
-
+   auto pos = neighbours_.size();
   if (z==0) {
-    neighbours8_.push_back(ref);
+    pos = neighbours8_++;
     //find only the 4 neighbours
     if (absx+absy==1)
-      neighbours4_.push_back(ref);
+      pos = neighbours4_++;
   }
+  neighbours_.insert(neighbours_.begin()+pos,ref);
+  neighbourInfos_.insert(neighbourInfos_.begin()+pos,bitmask);
+
+  assert( neighbours4_<5);
+  assert( neighbours8_<9);
+  assert( neighbours4_<=neighbours8_);
+  assert( neighbours8_<=neighbours_.size());
+
 }
 
 
-const PFRecHitRef PFRecHit::getNeighbour(short x,short y,short z) {
+unsigned int PFRecHit::getNeighbour(short x,short y,short z) const {
   unsigned short absx = abs(x);
   unsigned short absy = abs(y);
   unsigned short absz = abs(z);
@@ -61,7 +67,7 @@ const PFRecHitRef PFRecHit::getNeighbour(short x,short y,short z) {
     if (neighbourInfos_[i] == bitmask)
       return neighbours_[i];
   }
-  return PFRecHitRef();
+  return std::numeric_limits<unsigned int>::max();
 }
 
 }
diff --git a/DataFormats/ParticleFlowReco/src/classes_def_2.xml b/DataFormats/ParticleFlowReco/src/classes_def_2.xml
index e325ca5f965c5..99e2eed31b438 100644
--- a/DataFormats/ParticleFlowReco/src/classes_def_2.xml
+++ b/DataFormats/ParticleFlowReco/src/classes_def_2.xml
@@ -37,7 +37,7 @@
 
 
   <class name="reco::PFRecHit" ClassVersion="18">
-   <version ClassVersion="18" checksum=" 3589980878"/>
+   <version ClassVersion="18" checksum="3122773923"/>
    <version ClassVersion="17" checksum="2606786437"/>
    <version ClassVersion="10" checksum="1807687081"/>
    <version ClassVersion="11" checksum="3727193351"/>
@@ -51,7 +51,11 @@
    </ioread>
 
    <field name="caloCell_" transient="true"/>
-   
+   <field name="neighbours_" transient="true"/>
+   <field name="neighbours_info" transient="true"/>
+   <field name="neighbours4_" transient="true"/>
+   <field name="neighbours8_" transient="true"/>
+
   </class>
 
 
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h
index 52dd72d66b381..b57eddf118f90 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h
@@ -136,7 +136,7 @@ class PFRecHitCaloNavigatorWithTime : public PFRecHitNavigatorBase {
       sigma2 = _timeResolutionCalc->timeResolution2(hit.energy()) + _timeResolutionCalc->timeResolution2(found_hit->energy());
       const double deltaTime = hit.time()-found_hit->time();
       if(deltaTime*deltaTime<sigmaCut2_*sigma2) {
-	hit.addNeighbour(eta,phi,0,reco::PFRecHitRef(refProd,std::distance(hits->begin(),found_hit)));
+	hit.addNeighbour(eta,phi,0,found_hit-hits->begin());
       }
     }
 
diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h
index fc7a5167c4e22..f9da95ac1f5b8 100644
--- a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h
+++ b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h
@@ -46,7 +46,7 @@ class PFRecHitNavigatorBase {
 					return a.detId() < id;
 				      });
     if( found_hit != hits->end() && found_hit->detId() == id.rawId() ) {
-      hit.addNeighbour(eta,phi,depth,reco::PFRecHitRef(refProd,std::distance(hits->begin(),found_hit)));
+      hit.addNeighbour(eta,phi,depth,found_hit-hits->begin());
     }    
   }
 
diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
index a12a20882757f..19970373b6ba0 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
@@ -100,13 +100,13 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
   double depth = 0.0;  
   double position_norm = 0.0;
   double x(0.0),y(0.0),z(0.0);
-  const reco::PFRecHitRefVector* seedNeighbours = nullptr;
+  auto seedNeighbours = mySeed.hit->neighbours();
   switch( _posCalcNCrystals ) {
   case 5:
-    seedNeighbours = &mySeed.hit->neighbours4();
+    seedNeighbours = mySeed.hit->neighbours4();
     break;
   case 9:
-    seedNeighbours = &mySeed.hit->neighbours8();
+    seedNeighbours = mySeed.hit->neighbours8();
     break;
   default:
     break;
@@ -131,9 +131,9 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
   else {  // only seed and its neighbours
      compute(mySeed);
      // search seedNeighbours to find energy fraction in cluster (sic)
-     unInitDynArray(reco::PFRecHit const *,seedNeighbours->size(),nei);	  
-     for(auto k :seedNeighbours->refVector().keys()){ 
-      nei.push_back(&recHitCollection[k]);
+     unInitDynArray(reco::PFRecHit const *,seedNeighbours.size(),nei);	  
+     for(auto k : seedNeighbours){ 
+      nei.push_back(recHitCollection+k);
      }
      std::sort(nei.begin(),nei.end());
      struct LHitLess {
diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.cc
index a28df23c8106c..03b68889d0f33 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.cc
@@ -42,7 +42,7 @@ buildClusters(const edm::Handle<reco::PFRecHitCollection>& input,
     const int seed = idx_e.first;
     if( !rechitMask[seed] || !seedable[seed] || used[seed] ) continue;    
     temp.reset();
-    buildTopoCluster(input,rechitMask,makeRefhit(input,seed),used,temp);
+    buildTopoCluster(input,rechitMask,seed,used,temp);
     if( temp.recHitFractions().size() ) output.push_back(temp);
   }
 }
@@ -50,39 +50,41 @@ buildClusters(const edm::Handle<reco::PFRecHitCollection>& input,
 void Basic2DGenericTopoClusterizer::
 buildTopoCluster(const edm::Handle<reco::PFRecHitCollection>& input,
 		 const std::vector<bool>& rechitMask,
-		 const reco::PFRecHitRef& cell,
+		 unsigned int kcell,
 		 std::vector<bool>& used,		 
 		 reco::PFCluster& topocluster) {
-  int cell_layer = (int)cell->layer();
+  auto const & cell = (*input)[kcell];
+  int cell_layer = (int)cell.layer();
   if( cell_layer == PFLayer::HCAL_BARREL2 && 
-      std::abs(cell->positionREP().eta()) > 0.34 ) {
+      std::abs(cell.positionREP().eta()) > 0.34 ) {
       cell_layer *= 100;
     }    
   const std::pair<double,double>& thresholds =
       _thresholds.find(cell_layer)->second;
-  if( cell->energy() < thresholds.first || 
-      cell->pt2() < thresholds.second ) {
+  if( cell.energy() < thresholds.first || 
+      cell.pt2() < thresholds.second ) {
     LOGDRESSED("GenericTopoCluster::buildTopoCluster()")
-      << "RecHit " << cell->detId() << " with enegy " 
-      << cell->energy() << " GeV was rejected!." << std::endl;
+      << "RecHit " << cell.detId() << " with enegy " 
+      << cell.energy() << " GeV was rejected!." << std::endl;
     return;
   }
-
-  used[cell.key()] = true;
-  topocluster.addRecHitFraction(reco::PFRecHitFraction(cell, 1.0));
+  auto k = kcell;
+  used[k] = true;
+  auto ref = makeRefhit(input,k);
+  topocluster.addRecHitFraction(reco::PFRecHitFraction(ref, 1.0));
   
-  const reco::PFRecHitRefVector& neighbours = 
-    ( _useCornerCells ? cell->neighbours8() : cell->neighbours4() );
+  auto const & neighbours = 
+    ( _useCornerCells ? cell.neighbours8() : cell.neighbours4() );
   
-  for( const reco::PFRecHitRef nb : neighbours ) {
-    if( used[nb.key()] || !rechitMask[nb.key()] ) {
+  for( auto nb : neighbours ) {
+    if( used[nb] || !rechitMask[nb] ) {
       LOGDRESSED("GenericTopoCluster::buildTopoCluster()")
-      	<< "  RecHit " << cell->detId() << "\'s" 
-	<< " neighbor RecHit " << input->at(nb.key()).detId() 
+      	<< "  RecHit " << cell.detId() << "\'s" 
+	<< " neighbor RecHit " << input->at(nb).detId() 
 	<< " with enegy " 
-	<< input->at(nb.key()).energy() << " GeV was rejected!" 
-	<< " Reasons : " << used[nb.key()] << " (used) " 
-	<< !rechitMask[nb.key()] << " (masked)." << std::endl;
+	<< input->at(nb).energy() << " GeV was rejected!" 
+	<< " Reasons : " << used[nb] << " (used) " 
+	<< !rechitMask[nb] << " (masked)." << std::endl;
       continue;
     }
     buildTopoCluster(input,rechitMask,nb,used,topocluster);
diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.h b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.h
index 82f3b70dbc7e0..befa57a477c0b 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.h
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.h
@@ -23,7 +23,7 @@ class Basic2DGenericTopoClusterizer : public InitialClusteringStepBase {
   const bool _useCornerCells;
   void buildTopoCluster(const edm::Handle<reco::PFRecHitCollection>&,
 			const std::vector<bool>&, // masked rechits
-			const reco::PFRecHitRef&, //present rechit
+			unsigned int, //present rechit
 			std::vector<bool>&, // hit usage state
 			reco::PFCluster&); // the topocluster
   
diff --git a/RecoParticleFlow/PFClusterProducer/src/LocalMaximumSeedFinder.cc b/RecoParticleFlow/PFClusterProducer/src/LocalMaximumSeedFinder.cc
index 45724a658a5a2..467bb5b89dc0f 100644
--- a/RecoParticleFlow/PFClusterProducer/src/LocalMaximumSeedFinder.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/LocalMaximumSeedFinder.cc
@@ -7,7 +7,7 @@
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 
 namespace {
-  const reco::PFRecHitRefVector _noNeighbours;
+  const reco::PFRecHit::Neighbours  _noNeighbours(nullptr,0);
 }
 
 LocalMaximumSeedFinder::
@@ -81,36 +81,36 @@ findSeeds( const edm::Handle<reco::PFRecHitCollection>& input,
     if( !usable[idx] ) continue;
     //get the neighbours of this seed
     const auto & maybeseed = (*input)[idx];
-    const reco::PFRecHitRefVector* myNeighbours;
+    reco::PFRecHit::Neighbours  myNeighbours;
     switch( _nNeighbours ) {
     case -1:
-      myNeighbours = &maybeseed.neighbours();
+      myNeighbours = maybeseed.neighbours();
       break;
     case 0: // for HF clustering
-      myNeighbours = &_noNeighbours;
+      myNeighbours = _noNeighbours;
       break;
     case 4:
-      myNeighbours = &maybeseed.neighbours4();
+      myNeighbours = maybeseed.neighbours4();
       break;
     case 8:
-      myNeighbours = &maybeseed.neighbours8();
+      myNeighbours = maybeseed.neighbours8();
       break;
     default:
       throw cms::Exception("InvalidConfiguration")
 	<< "LocalMaximumSeedFinder only accepts nNeighbors = {-1,0,4,8}";    
     }
     seedable[idx] = true;
-    for( const reco::PFRecHitRef& neighbour : *myNeighbours ) {
-      if( !mask[neighbour.key()] ) continue;
-      if( energies[neighbour.key()] > energies[idx] ) {
+    for( auto neighbour : myNeighbours ) {
+      if( !mask[neighbour] ) continue;
+      if( energies[neighbour] > energies[idx] ) {
 //        std::cout << "how this can be?" << std::endl;
 	seedable[idx] = false;	
 	break;
       }
     }
     if( seedable[idx] ) {
-      for( const reco::PFRecHitRef& neighbour : *myNeighbours ) {
-	usable[neighbour.key()] = false;
+      for( auto neighbour : myNeighbours ) {
+	usable[neighbour] = false;
       }
     }
   }
diff --git a/RecoParticleFlow/PFClusterProducer/src/RBXAndHPDCleaner.cc b/RecoParticleFlow/PFClusterProducer/src/RBXAndHPDCleaner.cc
index 669ec1117eafb..ecd842c8b546c 100644
--- a/RecoParticleFlow/PFClusterProducer/src/RBXAndHPDCleaner.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/RBXAndHPDCleaner.cc
@@ -81,18 +81,19 @@ clean(const edm::Handle<reco::PFRecHitCollection>& input,
       totalEta2 = totalEta2W = totalPhi2 = totalPhi2W = totalEnergy2 = 1e-9;
       nSeeds = nSeeds0 = rechits.size();
       for( unsigned jh = 0; jh < rechits.size(); ++jh ) {
-	const reco::PFRecHit&  rechit = input->at(jh);
+	const reco::PFRecHit&  rechit = (*input)[jh];
 	// check if rechit is a seed
 	unsigned nN = 0 ; // neighbours over threshold
 	bool isASeed = true;
-	const reco::PFRecHitRefVector& neighbours4 = rechit.neighbours4();
-	for( const reco::PFRecHitRef& neighbour : neighbours4 ) {
-	  if( neighbour->energy() > rechit.energy() ) {	    
+	auto const & neighbours4 = rechit.neighbours4();
+	for( auto k : neighbours4 ) {
+          auto const & neighbour = (*input)[k]; 
+	  if( neighbour.energy() > rechit.energy() ) {	    
 	    --nSeeds; --nSeeds0;
 	    isASeed = false;
 	    break;
 	  } else {
-	    if( neighbour->energy() > 0.4 ) ++nN;
+	    if( neighbour.energy() > 0.4 ) ++nN;
 	  }
 	}
 	if ( isASeed && !nN ) --nSeeds0;
diff --git a/RecoParticleFlow/PFClusterProducer/src/SpikeAndDoubleSpikeCleaner.cc b/RecoParticleFlow/PFClusterProducer/src/SpikeAndDoubleSpikeCleaner.cc
index 7e7018e2e6027..765adbc6e124e 100644
--- a/RecoParticleFlow/PFClusterProducer/src/SpikeAndDoubleSpikeCleaner.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/SpikeAndDoubleSpikeCleaner.cc
@@ -106,7 +106,6 @@ clean(const edm::Handle<reco::PFRecHitCollection>& input,
   for( unsigned i = 0; i < hits.size(); ++i ) ordered_hits[i]=i;
   std::sort(ordered_hits.begin(),ordered_hits.end(),[&](unsigned i, unsigned j) { return hits[i].energy()>hits[j].energy();});
     
-
   for( const auto& idx : ordered_hits ) {
     const unsigned i = idx;
     if( !mask[i] ) continue; // don't need to re-mask things :-)
@@ -120,16 +119,17 @@ clean(const edm::Handle<reco::PFRecHitCollection>& input,
     if( rechit.energy() < clean._singleSpikeThresh ) continue;
     const double rhenergy = rechit.energy();
     // single spike cleaning
-    const reco::PFRecHitRefVector& neighbours4 = rechit.neighbours4();
+    auto const & neighbours4 = rechit.neighbours4();
     double surroundingEnergy = rechit.energy();
     double neighbourEnergy = 0.0;
     double layerEnergy = 0.0;
-    for( const reco::PFRecHitRef& neighbour : neighbours4 ) {
-      if( !mask[neighbour.key()] ) continue;
-      const double sum = neighbour->energy(); //energyUp is just rechit energy?
+    for( auto k : neighbours4 ) {
+      if( !mask[k] ) continue;
+      const auto & neighbour = hits[k];
+      const double sum = neighbour.energy(); //energyUp is just rechit energy?
       surroundingEnergy += sum;
       neighbourEnergy   += sum;
-      layerEnergy       += neighbour->energy();
+      layerEnergy       += neighbour.energy();
     }    
     //   wannaBeSeed.energyUp()/wannaBeSeed.energy() : 1.;
     // Fraction 1 is the balance between the hit and its neighbours 
@@ -159,30 +159,31 @@ clean(const edm::Handle<reco::PFRecHitCollection>& input,
       //Determine energy surrounding the seed and the most energetic neighbour
       double surroundingEnergyi = 0.0;      
       double enmax = -999.0;
-      reco::PFRecHitRef mostEnergeticNeighbour;
-      const reco::PFRecHitRefVector& neighbours4i = rechit.neighbours4();
-      for( reco::PFRecHitRef neighbour : neighbours4i ) {
-	if( !mask[neighbour.key()] ) continue;
-	const double nenergy = neighbour->energy();
+      unsigned int  mostEnergeticNeighbour=0;
+      auto const & neighbours4i = rechit.neighbours4();
+      for( auto k : neighbours4i ) {
+	if( !mask[k] ) continue;
+        const auto & neighbour = hits[k];
+	const double nenergy = neighbour.energy();
 	surroundingEnergyi += nenergy;
 	if( nenergy > enmax ) {
 	  enmax = nenergy;
-	  mostEnergeticNeighbour = neighbour;
+	  mostEnergeticNeighbour = k;
 	}
       }
       // is there an energetic neighbour
       if( enmax > 0.0 ) {
 	double surroundingEnergyj = 0.0;
-	const reco::PFRecHitRefVector& neighbours4j = 
-	  mostEnergeticNeighbour->neighbours4();
-	for( const reco::PFRecHitRef& neighbour : neighbours4j ) {
-	  //if( !mask[neighbour] &&  neighbour != i) continue; // leave out?
-	  surroundingEnergyj += neighbour->energy();
+	auto const & neighbours4j = 
+	  hits[mostEnergeticNeighbour].neighbours4();
+	for(auto k : neighbours4j ) {
+	  //if( !mask[k] &&  k != i) continue; // leave out?
+	  surroundingEnergyj += hits[k].energy();
 	}
 	// The energy surrounding the double spike candidate 
 	const double surroundingEnergyFraction = 
 	  (surroundingEnergyi+surroundingEnergyj) / 
-	  (rechit.energy()+mostEnergeticNeighbour->energy()) - 1.;
+	  (rechit.energy()+hits[mostEnergeticNeighbour].energy()) - 1.;
 	if ( surroundingEnergyFraction < clean._doubleSpikeS6S2 ) { 
 	  const double eta = rechit.positionREP().eta();
 	  const double aeta = std::abs(eta);
@@ -197,7 +198,7 @@ clean(const edm::Handle<reco::PFRecHitCollection>& input,
 		 surroundingEnergyFraction < clean._doubleSpikeS6S2/clean._fracThreshMod ) 
 		) ) {	    
 	    mask[i] = false;
-	    mask[mostEnergeticNeighbour.key()] = false;
+	    mask[mostEnergeticNeighbour] = false;
 	  }
 	}
       } // was there an energetic neighbour ? 

From f9f339a0c28211e52f3b500a167c563dcda77c37 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Sun, 17 Apr 2016 18:03:19 +0200
Subject: [PATCH 11/16] do a proper sort

---
 .../src/Basic2DGenericTopoClusterizer.cc      | 28 ++++++++-----------
 1 file changed, 11 insertions(+), 17 deletions(-)

diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.cc
index 03b68889d0f33..ddaf139806998 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericTopoClusterizer.cc
@@ -13,33 +13,27 @@
 #define LOGDRESSED(x) LogDebug(x)
 #endif
 
-namespace {
-  bool greaterByEnergy(const std::pair<unsigned,double>& a,
-		       const std::pair<unsigned,double>& b) {
-    return a.second > b.second;
-  }
-}
-
 void Basic2DGenericTopoClusterizer::
 buildClusters(const edm::Handle<reco::PFRecHitCollection>& input,
 	      const std::vector<bool>& rechitMask,
 	      const std::vector<bool>& seedable,
-	      reco::PFClusterCollection& output) {  
-  std::vector<bool> used(input->size(),false);
-  std::vector<std::pair<unsigned,double> > seeds;
+	      reco::PFClusterCollection& output) {
+  auto const & hits = *input;  
+  std::vector<bool> used(hits.size(),false);
+  std::vector<unsigned int> seeds;
   
   // get the seeds and sort them descending in energy
-  seeds.reserve(input->size());  
-  for( unsigned i = 0; i < input->size(); ++i ) {
+  seeds.reserve(hits.size());  
+  for( unsigned int i = 0; i < hits.size(); ++i ) {
     if( !rechitMask[i] || !seedable[i] || used[i] ) continue;
-    std::pair<unsigned,double> val = std::make_pair(i,input->at(i).energy());
-    auto pos = std::upper_bound(seeds.begin(),seeds.end(),val,greaterByEnergy);
-    seeds.insert(pos,val);
+    seeds.emplace_back(i);
   }
+  // maxHeap would be better
+  std::sort(seeds.begin(),seeds.end(),
+            [&](unsigned int i, unsigned int j) { return hits[i].energy()>hits[j].energy();});  
   
   reco::PFCluster temp;
-  for( const auto& idx_e : seeds ) {    
-    const int seed = idx_e.first;
+  for( auto seed : seeds ) {    
     if( !rechitMask[seed] || !seedable[seed] || used[seed] ) continue;    
     temp.reset();
     buildTopoCluster(input,rechitMask,seed,used,temp);

From cb74958faa42f535b2af1c0455a8ce0c8477700d Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Mon, 18 Apr 2016 10:37:11 +0200
Subject: [PATCH 12/16] move to floats in position computation

---
 .../src/Basic2DGenericPFlowPositionCalc.cc    | 37 ++++++++++---------
 .../src/Basic2DGenericPFlowPositionCalc.h     |  4 +-
 2 files changed, 22 insertions(+), 19 deletions(-)

diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
index 19970373b6ba0..0b1cb88d3da18 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
@@ -43,25 +43,23 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
 
   auto const recHitCollection = &(*cluster.recHitFractions()[0].recHitRef()) - cluster.recHitFractions()[0].recHitRef().key();
   auto nhits = cluster.recHitFractions().size();
-  struct LHit{ reco::PFRecHit const * hit; double energy; double fraction;};
+  struct LHit{ reco::PFRecHit const * hit; float energy; float fraction;};
   declareDynArray(LHit,nhits,hits);
   for(auto i=0U; i<nhits; ++i) {
     auto const & hf = cluster.recHitFractions()[i];
     auto k = hf.recHitRef().key();
     auto p = recHitCollection+k;
-    hits[i]= {p,(*p).energy(), hf.fraction()}; 
+    hits[i]= {p,(*p).energy(), float(hf.fraction())}; 
   }
 
-  if(_posCalcNCrystals != -1) // sorted to make neighbour search faster
-    std::sort(hits.begin(),hits.end(),[](LHit const& a, LHit const& b) { return a.hit<b.hit;});
 
-  LHit mySeed={nullptr}; 
+  LHit mySeed={nullptr};
   for( auto const & rhf : hits ) {
     const reco::PFRecHit & refhit = *rhf.hit;
-    if( refhit.detId() == cluster.seed() ) mySeed = rhf;
-    const double rh_fraction = rhf.fraction;
-    const double rh_rawenergy = rhf.energy;
-    const double rh_energy = rh_rawenergy * rh_fraction;   
+    if( refhit.detId() == cluster.seed() )  mySeed = rhf;
+    const auto rh_fraction = rhf.fraction;
+    const auto rh_rawenergy = rhf.energy;
+    const auto rh_energy = rh_rawenergy * rh_fraction;   
     if( edm::isNotFinite(rh_energy) ) {
       throw cms::Exception("PFClusterAlgo")
 	<<"rechit " << refhit.detId() << " has a NaN energy... " 
@@ -92,11 +90,12 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
       max_e_layer = refhit.layer();
     }    
   }
+  
   cluster.setEnergy(cl_energy);
   cluster.setTime(cl_time/cl_timeweight);
   cluster.setLayer(max_e_layer);
-  // calculate the position
 
+  // calculate the position
   double depth = 0.0;  
   double position_norm = 0.0;
   double x(0.0),y(0.0),z(0.0);
@@ -114,18 +113,22 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
 
   auto compute = [&](LHit const& rhf) {
     const reco::PFRecHit & refhit = *rhf.hit;  
-    const double rh_energy = rhf.energy * rhf.fraction;
-    const double norm = ( rhf.fraction < _minFractionInCalc ? 
+    const auto rh_energy = rhf.energy * rhf.fraction;
+    const auto norm = ( rhf.fraction < _minFractionInCalc ? 
 			  0.0 : 
-			  std::max(0.0,vdt::fast_log(rh_energy*_logWeightDenom)) );
-    const math::XYZPoint rhpos_xyz(refhit.position());
-    x += rhpos_xyz.X() * norm;
-    y += rhpos_xyz.Y() * norm;
-    z += rhpos_xyz.Z() * norm;
+			  std::max(0.0f,vdt::fast_logf(rh_energy*_logWeightDenom)) );
+    const auto rhpos_xyz = refhit.position();
+    x += rhpos_xyz.x() * norm;
+    y += rhpos_xyz.y() * norm;
+    z += rhpos_xyz.z() * norm;
     depth += refhit.depth()*norm;
     position_norm += norm;
   };
 
+  if(_posCalcNCrystals != -1) // sorted to make neighbour search faster (maybe)
+    std::sort(hits.begin(),hits.end(),[](LHit const& a, LHit const& b) { return a.hit<b.hit;});
+
+
   if(_posCalcNCrystals == -1)
     for( auto const & rhf : hits ) compute(rhf);
   else {  // only seed and its neighbours
diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.h b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.h
index d7d55c633da23..12248a3c07e1d 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.h
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.h
@@ -52,8 +52,8 @@ class Basic2DGenericPFlowPositionCalc : public PFCPositionCalculatorBase {
 
  private:
   const int _posCalcNCrystals;
-  const double _logWeightDenom;
-  const double _minAllowedNorm;
+  const float _logWeightDenom;
+  const float _minAllowedNorm;
   
   std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalcBarrel;
   std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalcEndcap;

From 783e0fbf06f2310ab8ad47ef189ab4754983f65a Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Mon, 18 Apr 2016 16:08:00 +0200
Subject: [PATCH 13/16] factorize conditionals and some division

---
 .../src/Basic2DGenericPFlowPositionCalc.cc    | 33 ++++++++++++-------
 .../PFClusterTools/src/LinkByRecHit.cc        |  9 ++---
 2 files changed, 26 insertions(+), 16 deletions(-)

diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
index 0b1cb88d3da18..c9a5a72013e55 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
@@ -10,6 +10,16 @@
 
 #include "vdt/vdtMath.h"
 
+
+namespace {
+  inline
+  bool isBarrel(int cell_layer){
+    return (cell_layer == PFLayer::HCAL_BARREL1 ||
+	    cell_layer == PFLayer::HCAL_BARREL2 ||
+	    cell_layer == PFLayer::ECAL_BARREL);
+  }  
+}
+
 void Basic2DGenericPFlowPositionCalc::
 calculateAndSetPosition(reco::PFCluster& cluster) {
   calculateAndSetPositionActual(cluster);
@@ -53,6 +63,7 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
   }
 
 
+  bool resGiven = bool(_timeResolutionCalcBarrel) & bool(_timeResolutionCalcEndcap);
   LHit mySeed={nullptr};
   for( auto const & rhf : hits ) {
     const reco::PFRecHit & refhit = *rhf.hit;
@@ -60,21 +71,19 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
     const auto rh_fraction = rhf.fraction;
     const auto rh_rawenergy = rhf.energy;
     const auto rh_energy = rh_rawenergy * rh_fraction;   
-    if( edm::isNotFinite(rh_energy) ) {
+#ifdef PF_DEBUG
+    if unlikely( edm::isNotFinite(rh_energy) ) {
       throw cms::Exception("PFClusterAlgo")
 	<<"rechit " << refhit.detId() << " has a NaN energy... " 
 	<< "The input of the particle flow clustering seems to be corrupted.";
     }
+#endif
     cl_energy += rh_energy;
     // If time resolution is given, calculated weighted average
-    if ( bool(_timeResolutionCalcBarrel) & bool(_timeResolutionCalcEndcap) ) {
+    if ( resGiven ) {
       double res2 = 1.e-4;
       int cell_layer = (int)refhit.layer();
-      if (cell_layer == PFLayer::HCAL_BARREL1 ||
-          cell_layer == PFLayer::HCAL_BARREL2 ||
-          cell_layer == PFLayer::ECAL_BARREL)
-        res2 = 1./_timeResolutionCalcBarrel->timeResolution2(rh_rawenergy);
-      else
+      res2 =  isBarrel(cell_layer) ? 1./_timeResolutionCalcBarrel->timeResolution2(rh_rawenergy) :
         res2 = 1./_timeResolutionCalcEndcap->timeResolution2(rh_rawenergy);
       cl_time += rh_fraction*refhit.time()*res2;
       cl_timeweight += rh_fraction*res2;
@@ -115,12 +124,12 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
     const reco::PFRecHit & refhit = *rhf.hit;  
     const auto rh_energy = rhf.energy * rhf.fraction;
     const auto norm = ( rhf.fraction < _minFractionInCalc ? 
-			  0.0 : 
+			  0.0f : 
 			  std::max(0.0f,vdt::fast_logf(rh_energy*_logWeightDenom)) );
-    const auto rhpos_xyz = refhit.position();
-    x += rhpos_xyz.x() * norm;
-    y += rhpos_xyz.y() * norm;
-    z += rhpos_xyz.z() * norm;
+    const auto rhpos_xyz = refhit.position()*norm;
+    x += rhpos_xyz.x();
+    y += rhpos_xyz.y();
+    z += rhpos_xyz.z();
     depth += refhit.depth()*norm;
     position_norm += norm;
   };
diff --git a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
index 6f319d1977509..81951bffab981 100644
--- a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
+++ b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
@@ -421,6 +421,7 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
   }
 
   // Get the rechits
+  auto zCorr = zPS/zECAL;
   const std::vector< reco::PFRecHitFraction >&  fracs = clusterECAL.recHitFractions();
   bool linkedbyrechit = false;
   //loop rechits
@@ -437,7 +438,7 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
     //getting rechit corners
     const auto &  corners = rechit_cluster.getCornersXYZ();
     
-    auto posxyz = rechit_cluster.position() * zPS/zECAL;
+    auto posxyz = rechit_cluster.position() * zCorr;
 #ifdef PFLOW_DEBUG
     if( debug ){
       std::cout << "Ecal rechit " << posxyz.X() << " "   << posxyz.Y() << std::endl;
@@ -449,10 +450,10 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
     double y[5];
     for ( unsigned jc=0; jc<4; ++jc ) {
       // corner position projected onto the preshower
-      auto cornerpos = corners[jc].basicVector() * zPS/zECAL;
+      auto cornerpos = corners[jc].basicVector() * zCorr;
       // Inflate the size by the size of the PS strips, and by 5% to include ECAL cracks.
-      x[3-jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/fabs((cornerpos.x()-posxyz.x()))*0.5*deltaX);
-      y[3-jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/fabs((cornerpos.y()-posxyz.y()))*0.5*deltaY);
+      x[3-jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/std::abs((cornerpos.x()-posxyz.x()))*0.5*deltaX);
+      y[3-jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/std::abs((cornerpos.y()-posxyz.y()))*0.5*deltaY);
       
 #ifdef PFLOW_DEBUG
       if( debug ){

From 7dd261a03872453a8c1d1d5413ed74256c5e713a Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Tue, 19 Apr 2016 10:43:21 +0200
Subject: [PATCH 14/16] vectorized (including icc compatibility)

---
 .../PFClusterTools/src/LinkByRecHit.cc        | 42 ++++++++++++++-----
 1 file changed, 31 insertions(+), 11 deletions(-)

diff --git a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
index 81951bffab981..d0f30ccb0475a 100644
--- a/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
+++ b/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc
@@ -3,6 +3,14 @@
 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
 #include "TMath.h"
 
+using  BVector2D = Basic2DVector<double>;
+using  Vector2D = Basic2DVector<double>::MathVector;
+namespace {
+  const Vector2D one2D=BVector2D(1.0,1.0).v;
+  const Vector2D fivepercent2D=BVector2D(0.05,0.05).v;
+}
+
+
 // to enable debugs
 //#define PFLOW_DEBUG
 
@@ -420,6 +428,8 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
     break;
   }
 
+  
+  auto deltaXY = BVector2D(deltaX,deltaY).v*0.5;
   // Get the rechits
   auto zCorr = zPS/zECAL;
   const std::vector< reco::PFRecHitFraction >&  fracs = clusterECAL.recHitFractions();
@@ -436,12 +446,12 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
     const reco::PFRecHit& rechit_cluster = *rh;
     
     //getting rechit corners
-    const auto &  corners = rechit_cluster.getCornersXYZ();
+    auto const & corners = rechit_cluster.getCornersXYZ();
     
-    auto posxyz = rechit_cluster.position() * zCorr;
+    auto posxy = BVector2D(rechit_cluster.position().xy()).v*zCorr;
 #ifdef PFLOW_DEBUG
     if( debug ){
-      std::cout << "Ecal rechit " << posxyz.X() << " "   << posxyz.Y() << std::endl;
+      std::cout << "Ecal rechit " << posxy.x() << " "   << posxy.y() << std::endl;
       std::cout << "PS cluster  " << xPS << " " << yPS << std::endl;
     }
 #endif
@@ -450,16 +460,25 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
     double y[5];
     for ( unsigned jc=0; jc<4; ++jc ) {
       // corner position projected onto the preshower
-      auto cornerpos = corners[jc].basicVector() * zCorr;
+      Vector2D cornerpos = BVector2D(corners[jc].basicVector().xy()).v*zCorr;
+      auto dist = (cornerpos-posxy);
+      auto adist = BVector2D(std::abs(dist[0]),std::abs(dist[1])).v; // all this beacuse icc does not support vector extension
       // Inflate the size by the size of the PS strips, and by 5% to include ECAL cracks.
-      x[3-jc] = cornerpos.x() + (cornerpos.x()-posxyz.x()) * (0.05 +1.0/std::abs((cornerpos.x()-posxyz.x()))*0.5*deltaX);
-      y[3-jc] = cornerpos.y() + (cornerpos.y()-posxyz.y()) * (0.05 +1.0/std::abs((cornerpos.y()-posxyz.y()))*0.5*deltaY);
-      
+      auto xy = cornerpos + (dist * (fivepercent2D +one2D/adist)*deltaXY);
+      /*
+      Vector2D xy(
+		  cornerpos.x() + (cornerpos.x()-posxy.x()) * (0.05 +1.0/std::abs((cornerpos.x()-posxy.x()))*deltaXY.x()),
+		  cornerpos.y() + (cornerpos.y()-posxy.y()) * (0.05 +1.0/std::abs((cornerpos.y()-posxy.y()))*deltaXY.y())
+		  );
+      */
+      x[3-jc] = xy[0];
+      y[3-jc] = xy[1];
+            
 #ifdef PFLOW_DEBUG
       if( debug ){
 	std::cout<<"corners "<<jc
-	    << " " << cornerpos.x() << " " << x[jc] 
-	    << " " << cornerpos.y() << " " << y[jc]
+	    << " " << cornerpos.x() << " " << x[3-jc] 
+	    << " " << cornerpos.y() << " " << y[3-jc]
 	    << std::endl;
       }
 #endif
@@ -485,8 +504,9 @@ LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
 #ifdef PFLOW_DEBUG
     if( debug ) std::cout << "Cluster PS and Cluster ECAL LINKED BY RECHIT" << std::endl;
 #endif
-    double dist = computeDist( xECAL/1000.,yECAL/1000.,
-			       xPS/1000.  ,yPS/1000, 
+    constexpr double scale = 1./1000.;
+    double dist = computeDist( xECAL*scale,yECAL*scale,
+			       xPS*scale  ,yPS*scale, 
 			       false);    
     return dist;
   } else { 

From 5946d0802b5aa3dcfa912863462b24c77c2476ce Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Mon, 25 Apr 2016 08:54:30 +0200
Subject: [PATCH 15/16] fix bug, remove obsolete code

---
 .../PFClusterProducer/plugins/PFCTRecHitProducer.cc            | 3 ---
 .../PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc   | 2 +-
 2 files changed, 1 insertion(+), 4 deletions(-)

diff --git a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
index 315012ba8db8f..cd2526bef00b4 100644
--- a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
+++ b/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc
@@ -858,13 +858,10 @@ PFCTRecHitProducer::createHcalRecHit( const DetId& detid,
     return 0;
   }
     
-  // double depth_correction = 0.;
   switch ( layer ) { 
   case PFLayer::HF_EM:
-    // depth_correction = position.z() > 0. ? EM_Depth_ : -EM_Depth_;
   case PFLayer::HF_HAD:
   {
-    //  depth_correction = position.z() > 0. ? HAD_Depth_ : -HAD_Depth_;
     auto zp = dynamic_cast<IdealZPrism const*>(thisCell);
     assert(zp);
     thisCell = zp->forPF();
diff --git a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
index c9a5a72013e55..c67b55cfa89f5 100644
--- a/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
+++ b/RecoParticleFlow/PFClusterProducer/src/Basic2DGenericPFlowPositionCalc.cc
@@ -84,7 +84,7 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) const {
       double res2 = 1.e-4;
       int cell_layer = (int)refhit.layer();
       res2 =  isBarrel(cell_layer) ? 1./_timeResolutionCalcBarrel->timeResolution2(rh_rawenergy) :
-        res2 = 1./_timeResolutionCalcEndcap->timeResolution2(rh_rawenergy);
+                                     1./_timeResolutionCalcEndcap->timeResolution2(rh_rawenergy);
       cl_time += rh_fraction*refhit.time()*res2;
       cl_timeweight += rh_fraction*res2;
     }

From 299d4876094913bf8fa05e92f86e69226732d206 Mon Sep 17 00:00:00 2001
From: Vincenzo Innocente <vincenzo.innocente@cern.ch>
Date: Tue, 26 Apr 2016 16:48:29 +0200
Subject: [PATCH 16/16] fix calo writers

---
 CondTools/Geometry/plugins/calowriters.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CondTools/Geometry/plugins/calowriters.cc b/CondTools/Geometry/plugins/calowriters.cc
index b6eaf19ea1364..35b825a71bfe1 100644
--- a/CondTools/Geometry/plugins/calowriters.cc
+++ b/CondTools/Geometry/plugins/calowriters.cc
@@ -75,7 +75,7 @@ CaloGeometryDBEP<HcalGeometry, CaloGeometryDBWriter>::produceAligned( const type
     
     ptr->fillDefaultNamedParameters() ;
 
-    ptr->allocateCorners( hcalTopology->ncells() ) ;
+    ptr->allocateCorners( hcalTopology->ncells() + hcalTopology->getHFSize() ) ;
 
     ptr->allocatePar(    dvec.size() ,
 			 HcalGeometry::k_NumberOfParametersPerShape ) ;