HGCal3DClustering.h

Go to the documentation of this file.

#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 "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.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),
        es(0),
        clusterTools(std::make_unique<hgcal::ClusterTools>(conf, sumes)),
        caloGeomToken_(sumes.esConsumes<CaloGeometry, CaloGeometryRecord>()) {}

  HGCal3DClustering(const edm::ParameterSet& conf, edm::ConsumesCollector& sumes)
      : HGCal3DClustering(conf,
                          sumes,
                          conf.getParameter<std::vector<double>>("multiclusterRadii"),
                          conf.getParameter<unsigned>("minClusters")) {}

  void getEvent(const edm::Event& ev) { clusterTools->getEvent(ev); }
  void getEventSetup(const edm::EventSetup& es) {
    clusterTools->getEventSetup(es);
    const CaloGeometry& geom = es.getData(caloGeomToken_);
    rhtools_.setGeometry(geom);
    maxlayer = rhtools_.lastLayerBH();
    points.clear();
    minpos.clear();
    maxpos.clear();
    zees.clear();
    points.resize(2 * (maxlayer + 1));
    minpos.resize(2 * (maxlayer + 1), {{0.0f, 0.0f}});
    maxpos.resize(2 * (maxlayer + 1), {{0.0f, 0.0f}});
    zees.resize(2 * (maxlayer + 1), 0.);
  }

  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
  unsigned int 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> KDTree;
  typedef KDTreeNodeInfo<ClusterRef> 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;                            
  std::vector<float> zees;                           
  std::unique_ptr<hgcal::ClusterTools> clusterTools; 
  hgcal::RecHitTools rhtools_;                       
  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
};

#endif