Merge pull request #19026 from clelange/PR18236Backport91X

PR #18236 backport to 91X

RecoLocalCalo/HGCalRecAlgos/interface/ClusterTools.h

@@ -3,6 +3,7 @@
 
 #include <array>
 #include <cmath>
+#include <numeric>
 
 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
@@ -23,25 +24,39 @@ namespace edm {
 }
 
 namespace hgcal {
-  class ClusterTools {    
+  class ClusterTools {
   public:
+    ClusterTools();
     ClusterTools(const edm::ParameterSet&, edm::ConsumesCollector&);
     ~ClusterTools() {}
 
     void getEvent(const edm::Event&);
     void getEventSetup(const edm::EventSetup&);
-    
+
     float getClusterHadronFraction(const reco::CaloCluster&) const;
 
-    math::XYZPoint getMultiClusterPosition(const reco::HGCalMultiCluster&, double vz = 0.) const;
-
-
+    math::XYZPoint getMultiClusterPosition(const reco::HGCalMultiCluster&) const;
+
+    int getLayer(const DetId) const;
+
     double getMultiClusterEnergy(const reco::HGCalMultiCluster&) const;
 
+    // only for EE
+    bool getWidths(const reco::CaloCluster & clus,double & sigmaetaeta, double & sigmaphiphi, double & sigmaetaetalog, double & sigmaphiphilog ) const;
   private:
+
+    std::vector<size_t> sort_by_z(const reco::HGCalMultiCluster&v) const {
+      std::vector<size_t> idx(v.size());
+      std::iota (std::begin(idx), std::end(idx), 0);
+      sort(idx.begin(), idx.end(),
+	   [&v](size_t i1, size_t i2) {return v.clusters()[i1]->z() < v.clusters()[i2]->z();});
+      return idx;
+    }
+
     RecHitTools rhtools_;
     const edm::EDGetTokenT<HGCRecHitCollection> eetok, fhtok, bhtok;
     const HGCRecHitCollection *eerh_, *fhrh_, *bhrh_;
+    static const int lastLayerEE = 28;
   };
 }
 

RecoLocalCalo/HGCalRecAlgos/interface/HGCal3DClustering.h

@@ -0,0 +1,94 @@
+#ifndef RecoLocalCalo_HGCalRecAlgos_HGCal3DClustering
+#define RecoLocalCalo_HGCalRecAlgos_HGCal3DClustering
+
+
+#include "DataFormats/Math/interface/Point3D.h"
+#include "DataFormats/EgammaReco/interface/BasicCluster.h"
+#include "DataFormats/ParticleFlowReco/interface/HGCalMultiCluster.h"
+
+#include <vector>
+#include <array>
+
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/ClusterTools.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/HGCalImagingAlgo.h"
+
+#include "KDTreeLinkerAlgoT.h"
+
+class HGCal3DClustering
+{
+public:
+
+ HGCal3DClustering() : radii({0.,0.,0.}), minClusters(0), clusterTools(nullptr)
+  {
+  }
+
+ HGCal3DClustering(const edm::ParameterSet& conf, edm::ConsumesCollector& sumes, const std::vector<double>& radii_in, uint32_t min_clusters) :
+  radii(radii_in),
+  minClusters(min_clusters),
+  points(2*(maxlayer+1)),
+  minpos(2*(maxlayer+1),{ {0.0f,0.0f} }),
+  maxpos(2*(maxlayer+1),{ {0.0f,0.0f} }),
+  es(0),
+  zees(2*(maxlayer+1),0.),
+  clusterTools(std::make_unique<hgcal::ClusterTools>(conf,sumes))
+  {
+  }
+
+  void getEvent(const edm::Event& ev) { clusterTools->getEvent(ev); }
+  void getEventSetup(const edm::EventSetup& es)
+  {
+    clusterTools->getEventSetup(es);
+    rhtools_.getEventSetup(es);
+  }
+
+  typedef std::vector<reco::BasicCluster> ClusterCollection;
+
+  std::vector<reco::HGCalMultiCluster> makeClusters(const reco::HGCalMultiCluster::ClusterCollection &);
+
+private:
+
+  void organizeByLayer(const reco::HGCalMultiCluster::ClusterCollection &);
+  void reset(){
+    for( auto& it: points)
+      {
+        it.clear();
+        std::vector<KDNode>().swap(it);
+      }
+    std::fill(zees.begin(), zees.end(), 0.);
+    for(unsigned int i = 0; i < minpos.size(); i++)
+      {
+	minpos[i][0]=0.;minpos[i][1]=0.;
+	maxpos[i][0]=0.;maxpos[i][1]=0.;
+      }
+  }
+  void layerIntersection(std::array<double,3> &to, const std::array<double,3> &from) const;
+
+  //max number of layers
+  static const unsigned int maxlayer = HGCalImagingAlgo::maxlayer;
+
+  std::vector<double> radii;
+  uint32_t minClusters;
+  struct ClusterRef {
+    int ind;
+    float z;
+    ClusterRef(int ind_i, float z_i): ind(ind_i),z(z_i){}
+    ClusterRef(): ind(-1),z(0.){}
+  };
+
+  typedef KDTreeLinkerAlgo<ClusterRef,2> KDTree;
+  typedef KDTreeNodeInfoT<ClusterRef,2> KDNode;
+  std::vector< std::vector<KDNode> > points;
+  std::vector<std::array<float,2> > minpos;
+  std::vector<std::array<float,2> > maxpos;
+  std::vector<size_t> es; /*!< vector to contain sorted indices of all clusters. */
+  std::vector<float> zees; /*!< vector to contain z position of each layer. */
+  std::unique_ptr<hgcal::ClusterTools> clusterTools; /*!< instance of tools to simplify cluster access. */
+  hgcal::RecHitTools rhtools_; /*!< instance of tools to access RecHit information. */
+  static const unsigned int lastLayerEE = 28;
+  static const unsigned int lastLayerFH = 40;
+  static const unsigned int lastLayerBH = 52;
+
+};
+
+#endif

RecoLocalCalo/HGCalRecAlgos/interface/HGCalDepthPreClusterer.h

@@ -10,17 +10,18 @@
 
 #include "RecoLocalCalo/HGCalRecAlgos/interface/ClusterTools.h"
 
-class HGCalDepthPreClusterer 
+class HGCalDepthPreClusterer
 {
 public:
-  
- HGCalDepthPreClusterer() : radius(0.), minClusters(0.), clusterTools(nullptr)
+
+ HGCalDepthPreClusterer() : radii({0.,0.,0.,}), minClusters(0), realSpaceCone(false), clusterTools(nullptr)
     {
   }
-  
- HGCalDepthPreClusterer(const edm::ParameterSet& conf, edm::ConsumesCollector& sumes, double radius_in, uint32_t min_clusters) : 
-  radius(radius_in),
+
+ HGCalDepthPreClusterer(const edm::ParameterSet& conf, edm::ConsumesCollector& sumes, std::vector<float> radii_in, uint32_t min_clusters, bool real_space_cone) :
+  radii(radii_in),
   minClusters(min_clusters),
+  realSpaceCone(real_space_cone),
   clusterTools(std::make_unique<hgcal::ClusterTools>(conf,sumes)) {
   }
 
@@ -30,14 +31,19 @@ class HGCalDepthPreClusterer
   typedef std::vector<reco::BasicCluster> ClusterCollection;
   //  typedef std::vector<reco::BasicCluster> MultiCluster;
 
-  std::vector<reco::HGCalMultiCluster> makePreClusters(const reco::HGCalMultiCluster::ClusterCollection &) const;    
+  std::vector<reco::HGCalMultiCluster> makePreClusters(const reco::HGCalMultiCluster::ClusterCollection &) const;
 
-private:  
-  float radius;
+private:
+  std::vector<float> radii;
   uint32_t minClusters;
-
+  bool realSpaceCone; /*!< flag to use cartesian space clustering. */
+
   std::unique_ptr<hgcal::ClusterTools> clusterTools;
 
+  static const unsigned int lastLayerEE = 28;
+  static const unsigned int lastLayerFH = 40;
+  static const unsigned int lastLayerBH = 52;
+
 };
 
 #endif