Back to home page

Project CMSSW displayed by LXR

 
 

    


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 // C/C++ headers
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   // this is the method that will start the clusterisation (it is possible to invoke this method
0057   // more than once - but make sure it is with different hit collections (or else use reset)
0058 
0059   void makeClusters() override;
0060 
0061   // this is the method to get the cluster collection out
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   /// point in the space
0092   typedef math::XYZPoint Point;
0093 
0094 private:
0095   // To get ECAL sigmaNoise
0096   edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> tok_ebThresholds_;
0097   const EcalPFRecHitThresholds* ebThresholds_;
0098   //To get HCAL sigmaNoise
0099   edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> tok_hcalThresholds_;
0100   const HcalPFCuts* hcalThresholds_;
0101   // The two parameters used to identify clusters
0102   double deltac_;
0103   double kappa_;
0104   double fractionCutoff_;
0105   int maxLayerIndex_;
0106   float outlierDeltaFactor_;
0107   bool doSharing_;
0108   Density density_;
0109   // For keeping the density per hit
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 {  // distance squared
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 {  // 2-d distance on the layer (x-y)
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 // explicit template instantiation
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