File indexing completed on 2025-07-09 05:00:28
0001 #ifndef RecoParticleFlow_PFClusterProducer_BarrelCLUEAlgo_h
0002 #define RecoParticleFlow_PFClusterProducer_BarrelCLUEAlgo_h
0003
0004 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalClusteringAlgoBase.h"
0005
0006 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0007
0008 #include "DataFormats/DetId/interface/DetId.h"
0009 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0010 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
0011
0012 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0013 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0014
0015 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0016 #include "DataFormats/Math/interface/Point3D.h"
0017 #include "DataFormats/Math/interface/deltaPhi.h"
0018
0019 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h"
0020
0021 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0022
0023 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0024 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0025 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0026 #include "CondFormats/HcalObjects/interface/HcalPFCuts.h"
0027
0028
0029 #include <set>
0030 #include <string>
0031 #include <vector>
0032
0033 using Density = hgcal_clustering::Density;
0034
0035 template <typename TILE>
0036 class BarrelCLUEAlgoT : public HGCalClusteringAlgoBase {
0037 public:
0038 BarrelCLUEAlgoT(const edm::ParameterSet& ps)
0039 : HGCalClusteringAlgoBase(
0040 (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
0041 reco::CaloCluster::undefined),
0042 deltac_(ps.getParameter<double>("deltac")),
0043 kappa_(ps.getParameter<double>("kappa")),
0044 fractionCutoff_(ps.getParameter<double>("fractionCutoff")),
0045 maxLayerIndex_(ps.getParameter<int>("maxLayerIndex")),
0046 outlierDeltaFactor_(ps.getParameter<double>("outlierDeltaFactor")),
0047 doSharing_(ps.getParameter<bool>("doSharing")) {}
0048 ~BarrelCLUEAlgoT() override {}
0049
0050 void getEventSetupPerAlgorithm(const edm::EventSetup& es) override;
0051 void setThresholds(edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd>,
0052 edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd>) override;
0053
0054 void populate(const HGCRecHitCollection& hits) override {};
0055 void populate(const reco::PFRecHitCollection& hits) override;
0056
0057
0058
0059 void makeClusters() override;
0060
0061
0062 std::vector<reco::BasicCluster> getClusters(bool) override;
0063
0064 void reset() override {
0065 clusters_v_.clear();
0066 clusters_v_.shrink_to_fit();
0067 for (auto& cl : numberOfClustersPerLayer_) {
0068 cl = 0;
0069 }
0070
0071 for (auto& cells : cells_) {
0072 cells.clear();
0073 cells.shrink_to_fit();
0074 }
0075 density_.clear();
0076 }
0077
0078 Density getDensity() override;
0079
0080 void computeThreshold();
0081
0082 static void fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0083 iDesc.add<double>("outlierDeltaFactor", 2.);
0084 iDesc.add<double>("kappa", 1.34);
0085 iDesc.add<int>("maxLayerIndex", 0);
0086 iDesc.add<double>("deltac", 0.0175);
0087 iDesc.add<double>("fractionCutoff", 0.);
0088 iDesc.add<bool>("doSharing", false);
0089 }
0090
0091
0092 typedef math::XYZPoint Point;
0093
0094 private:
0095
0096 edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> tok_ebThresholds_;
0097 const EcalPFRecHitThresholds* ebThresholds_;
0098
0099 edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> tok_hcalThresholds_;
0100 const HcalPFCuts* hcalThresholds_;
0101
0102 double deltac_;
0103 double kappa_;
0104 double fractionCutoff_;
0105 int maxLayerIndex_;
0106 float outlierDeltaFactor_;
0107 bool doSharing_;
0108 Density density_;
0109
0110
0111 struct BarrelCellsOnLayer {
0112 std::vector<DetId> detid;
0113 std::vector<float> eta;
0114 std::vector<float> phi;
0115 std::vector<float> r;
0116
0117 std::vector<float> weight;
0118 std::vector<float> rho;
0119
0120 std::vector<float> delta;
0121 std::vector<int> nearestHigher;
0122 std::vector<std::vector<int>> clusterIndex;
0123 std::vector<float> sigmaNoise;
0124 std::vector<std::vector<int>> followers;
0125 std::vector<bool> isSeed;
0126 std::vector<int> seedToCellIndex;
0127
0128 void clear() {
0129 detid.clear();
0130 eta.clear();
0131 phi.clear();
0132 r.clear();
0133 weight.clear();
0134 rho.clear();
0135 delta.clear();
0136 nearestHigher.clear();
0137 clusterIndex.clear();
0138 sigmaNoise.clear();
0139 followers.clear();
0140 isSeed.clear();
0141 seedToCellIndex.clear();
0142 }
0143
0144 void shrink_to_fit() {
0145 detid.shrink_to_fit();
0146 r.shrink_to_fit();
0147 eta.shrink_to_fit();
0148 phi.shrink_to_fit();
0149 weight.shrink_to_fit();
0150 rho.shrink_to_fit();
0151 delta.shrink_to_fit();
0152 nearestHigher.shrink_to_fit();
0153 clusterIndex.shrink_to_fit();
0154 sigmaNoise.shrink_to_fit();
0155 followers.shrink_to_fit();
0156 isSeed.shrink_to_fit();
0157 seedToCellIndex.clear();
0158 }
0159 };
0160
0161 std::vector<BarrelCellsOnLayer> cells_;
0162
0163 std::vector<int> numberOfClustersPerLayer_;
0164
0165 inline float distance2(int cell1, int cell2, int layerId) const {
0166 const float dphi = reco::deltaPhi(cells_[layerId].phi[cell1], cells_[layerId].phi[cell2]);
0167 const float deta = cells_[layerId].eta[cell1] - cells_[layerId].eta[cell2];
0168 return (deta * deta + dphi * dphi);
0169 }
0170
0171 inline float distance(int cell1, int cell2, int layerId) const {
0172 return std::sqrt(distance2(cell1, cell2, layerId));
0173 }
0174
0175 void prepareDataStructures(const unsigned int layerId);
0176 void calculateLocalDensity(const TILE& lt, const unsigned int layerId, float delta_c);
0177 void calculateDistanceToHigher(const TILE& lt, const unsigned int layerId, float delta_c);
0178 int findAndAssignClusters(const unsigned int layerId, float delta_c);
0179 void passSharedClusterIndex(const TILE& lt, const unsigned int layerId, float delta_c);
0180 math::XYZPoint calculatePosition(const std::vector<std::pair<int, float>>& v, const unsigned int layerId) const;
0181 void setDensity(const unsigned int layerId);
0182 };
0183
0184
0185 extern template class BarrelCLUEAlgoT<EBLayerTiles>;
0186 extern template class BarrelCLUEAlgoT<HBLayerTiles>;
0187
0188 using EBCLUEAlgo = BarrelCLUEAlgoT<EBLayerTiles>;
0189 using HBCLUEAlgo = BarrelCLUEAlgoT<HBLayerTiles>;
0190
0191 #endif