File indexing completed on 2023-06-02 00:50:54
0001 #ifndef RecoLocalCalo_HGCalRecProducers_HGCalCLUEAlgo_h
0002 #define RecoLocalCalo_HGCalRecProducers_HGCalCLUEAlgo_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/HGCRecHit/interface/HGCRecHitCollections.h"
0010 #include "Geometry/CaloTopology/interface/HGCalTopology.h"
0011 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0012 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0013
0014 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0015 #include "DataFormats/Math/interface/Point3D.h"
0016 #include "DataFormats/Math/interface/deltaPhi.h"
0017
0018 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h"
0019 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalCLUEStrategy.h"
0020
0021 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0022
0023
0024 #include <set>
0025 #include <string>
0026 #include <vector>
0027
0028 template <typename TILE, typename STRATEGY>
0029 class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
0030 public:
0031 HGCalCLUEAlgoT(const edm::ParameterSet& ps)
0032 : HGCalClusteringAlgoBase(
0033 (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
0034 reco::CaloCluster::undefined),
0035 vecDeltas_(ps.getParameter<std::vector<double>>("deltac")),
0036 kappa_(ps.getParameter<double>("kappa")),
0037 ecut_(ps.getParameter<double>("ecut")),
0038 dependSensor_(ps.getParameter<bool>("dependSensor")),
0039 dEdXweights_(ps.getParameter<std::vector<double>>("dEdXweights")),
0040 thicknessCorrection_(ps.getParameter<std::vector<double>>("thicknessCorrection")),
0041 sciThicknessCorrection_(ps.getParameter<double>("sciThicknessCorrection")),
0042 deltasi_index_regemfac_(ps.getParameter<int>("deltasi_index_regemfac")),
0043 maxNumberOfThickIndices_(ps.getParameter<unsigned>("maxNumberOfThickIndices")),
0044 fcPerMip_(ps.getParameter<std::vector<double>>("fcPerMip")),
0045 fcPerEle_(ps.getParameter<double>("fcPerEle")),
0046 nonAgedNoises_(ps.getParameter<std::vector<double>>("noises")),
0047 noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
0048 use2x2_(ps.getParameter<bool>("use2x2")),
0049 initialized_(false) {}
0050
0051 ~HGCalCLUEAlgoT() override {}
0052
0053 void getEventSetupPerAlgorithm(const edm::EventSetup& es) override;
0054
0055 void populate(const HGCRecHitCollection& hits) override;
0056
0057
0058
0059
0060 void makeClusters() override;
0061
0062
0063 std::vector<reco::BasicCluster> getClusters(bool) override;
0064
0065 void reset() override {
0066 clusters_v_.clear();
0067 clusters_v_.shrink_to_fit();
0068 for (auto& cl : numberOfClustersPerLayer_) {
0069 cl = 0;
0070 }
0071
0072 for (auto& cells : cells_) {
0073 cells.clear();
0074 cells.shrink_to_fit();
0075 }
0076 }
0077
0078 void computeThreshold();
0079
0080 static void fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0081 iDesc.add<std::vector<double>>("thresholdW0", {2.9, 2.9, 2.9});
0082 iDesc.add<double>("positionDeltaRho2", 1.69);
0083 iDesc.add<std::vector<double>>("deltac",
0084 {
0085 1.3,
0086 1.3,
0087 1.3,
0088 0.0315,
0089 });
0090 iDesc.add<bool>("dependSensor", true);
0091 iDesc.add<double>("ecut", 3.0);
0092 iDesc.add<double>("kappa", 9.0);
0093 iDesc.addUntracked<unsigned int>("verbosity", 3);
0094 iDesc.add<std::vector<double>>("dEdXweights", {});
0095 iDesc.add<std::vector<double>>("thicknessCorrection", {});
0096 iDesc.add<double>("sciThicknessCorrection", 0.9);
0097 iDesc.add<int>("deltasi_index_regemfac", 3);
0098 iDesc.add<unsigned>("maxNumberOfThickIndices", 6);
0099 iDesc.add<std::vector<double>>("fcPerMip", {});
0100 iDesc.add<double>("fcPerEle", 0.0);
0101 iDesc.add<std::vector<double>>("noises", {});
0102 edm::ParameterSetDescription descNestedNoiseMIP;
0103 descNestedNoiseMIP.add<bool>("scaleByDose", false);
0104 descNestedNoiseMIP.add<unsigned int>("scaleByDoseAlgo", 0);
0105 descNestedNoiseMIP.add<double>("scaleByDoseFactor", 1.);
0106 descNestedNoiseMIP.add<std::string>("doseMap", "");
0107 descNestedNoiseMIP.add<std::string>("sipmMap", "");
0108 descNestedNoiseMIP.add<double>("referenceIdark", -1);
0109 descNestedNoiseMIP.add<double>("referenceXtalk", -1);
0110 descNestedNoiseMIP.add<double>("noise_MIP", 1. / 100.);
0111 iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
0112 iDesc.add<bool>("use2x2", true);
0113 }
0114
0115
0116 typedef math::XYZPoint Point;
0117
0118 private:
0119
0120 std::vector<double> vecDeltas_;
0121 double kappa_;
0122
0123
0124 double ecut_;
0125
0126
0127
0128 bool dependSensor_;
0129 std::vector<double> dEdXweights_;
0130 std::vector<double> thicknessCorrection_;
0131 double sciThicknessCorrection_;
0132 int deltasi_index_regemfac_;
0133 unsigned maxNumberOfThickIndices_;
0134 std::vector<double> fcPerMip_;
0135 double fcPerEle_;
0136 std::vector<double> nonAgedNoises_;
0137 double noiseMip_;
0138 std::vector<std::vector<double>> thresholds_;
0139 std::vector<std::vector<double>> v_sigmaNoise_;
0140
0141 bool use2x2_;
0142
0143
0144 bool initialized_;
0145
0146 float outlierDeltaFactor_ = 2.f;
0147
0148 struct CellsOnLayer {
0149 std::vector<DetId> detid;
0150 std::vector<float> dim1;
0151 std::vector<float> dim2;
0152
0153 std::vector<float> weight;
0154 std::vector<float> rho;
0155
0156 std::vector<float> delta;
0157 std::vector<int> nearestHigher;
0158 std::vector<int> clusterIndex;
0159 std::vector<float> sigmaNoise;
0160 std::vector<std::vector<int>> followers;
0161 std::vector<bool> isSeed;
0162
0163 void clear() {
0164 detid.clear();
0165 dim1.clear();
0166 dim2.clear();
0167 weight.clear();
0168 rho.clear();
0169 delta.clear();
0170 nearestHigher.clear();
0171 clusterIndex.clear();
0172 sigmaNoise.clear();
0173 followers.clear();
0174 isSeed.clear();
0175 }
0176
0177 void shrink_to_fit() {
0178 detid.shrink_to_fit();
0179 dim1.shrink_to_fit();
0180 dim2.shrink_to_fit();
0181 weight.shrink_to_fit();
0182 rho.shrink_to_fit();
0183 delta.shrink_to_fit();
0184 nearestHigher.shrink_to_fit();
0185 clusterIndex.shrink_to_fit();
0186 sigmaNoise.shrink_to_fit();
0187 followers.shrink_to_fit();
0188 isSeed.shrink_to_fit();
0189 }
0190 };
0191
0192 std::vector<CellsOnLayer> cells_;
0193
0194 std::vector<int> numberOfClustersPerLayer_;
0195
0196 inline float distance(const TILE& lt, int cell1, int cell2, int layerId) const {
0197 return std::sqrt(lt.distance2(cells_[layerId].dim1[cell1],
0198 cells_[layerId].dim2[cell1],
0199 cells_[layerId].dim1[cell2],
0200 cells_[layerId].dim2[cell2]));
0201 }
0202
0203 void prepareDataStructures(const unsigned int layerId);
0204 void calculateLocalDensity(const TILE& lt, const unsigned int layerId,
0205 float delta);
0206 void calculateLocalDensity(const TILE& lt, const unsigned int layerId, float delta, HGCalSiliconStrategy strategy);
0207 void calculateLocalDensity(const TILE& lt,
0208 const unsigned int layerId,
0209 float delta,
0210 HGCalScintillatorStrategy strategy);
0211 void calculateDistanceToHigher(const TILE& lt, const unsigned int layerId, float delta);
0212 int findAndAssignClusters(const unsigned int layerId, float delta);
0213 };
0214
0215
0216 extern template class HGCalCLUEAlgoT<HGCalSiliconLayerTiles, HGCalSiliconStrategy>;
0217 extern template class HGCalCLUEAlgoT<HGCalScintillatorLayerTiles, HGCalScintillatorStrategy>;
0218 extern template class HGCalCLUEAlgoT<HFNoseLayerTiles, HGCalSiliconStrategy>;
0219
0220 using HGCalSiCLUEAlgo = HGCalCLUEAlgoT<HGCalSiliconLayerTiles, HGCalSiliconStrategy>;
0221 using HGCalSciCLUEAlgo = HGCalCLUEAlgoT<HGCalScintillatorLayerTiles, HGCalScintillatorStrategy>;
0222 using HFNoseCLUEAlgo = HGCalCLUEAlgoT<HFNoseLayerTiles, HGCalSiliconStrategy>;
0223
0224 #endif