File indexing completed on 2021-04-27 00:16:56
0001 #ifndef RecoLocalCalo_HGCalRecProducers_HGCalImagingAlgo_h
0002 #define RecoLocalCalo_HGCalRecProducers_HGCalImagingAlgo_h
0003
0004 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalClusteringAlgoBase.h"
0005
0006 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0007
0008 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0009 #include "Geometry/CaloTopology/interface/HGCalTopology.h"
0010 #include "DataFormats/DetId/interface/DetId.h"
0011 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0012 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0013 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0014 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0015 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0016 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0017
0018 #include "DataFormats/Math/interface/Point3D.h"
0019 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0020
0021 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0022 #include "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.h"
0023
0024
0025 #include <string>
0026 #include <vector>
0027 #include <set>
0028
0029 using Density = hgcal_clustering::Density;
0030
0031 class HGCalImagingAlgo : public HGCalClusteringAlgoBase {
0032 public:
0033 HGCalImagingAlgo(const edm::ParameterSet &ps, edm::ConsumesCollector iC)
0034 : HGCalClusteringAlgoBase(
0035 (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
0036 reco::CaloCluster::undefined,
0037 iC),
0038 thresholdW0_(ps.getParameter<std::vector<double>>("thresholdW0")),
0039 positionDeltaRho_c_(ps.getParameter<std::vector<double>>("positionDeltaRho_c")),
0040 vecDeltas_(ps.getParameter<std::vector<double>>("deltac")),
0041 kappa_(ps.getParameter<double>("kappa")),
0042 ecut_(ps.getParameter<double>("ecut")),
0043 sigma2_(1.0),
0044 dependSensor_(ps.getParameter<bool>("dependSensor")),
0045 dEdXweights_(ps.getParameter<std::vector<double>>("dEdXweights")),
0046 thicknessCorrection_(ps.getParameter<std::vector<double>>("thicknessCorrection")),
0047 fcPerMip_(ps.getParameter<std::vector<double>>("fcPerMip")),
0048 fcPerEle_(ps.getParameter<double>("fcPerEle")),
0049 nonAgedNoises_(ps.getParameter<edm::ParameterSet>("noises").getParameter<std::vector<double>>("values")),
0050 noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
0051 initialized_(false) {}
0052
0053 ~HGCalImagingAlgo() override {}
0054
0055 void getEventSetupPerAlgorithm(const edm::EventSetup &es) override;
0056
0057 void populate(const HGCRecHitCollection &hits) override;
0058
0059
0060
0061 void makeClusters() override;
0062
0063
0064 std::vector<reco::BasicCluster> getClusters(bool) override;
0065
0066
0067 void reset() override {
0068 clusters_v_.clear();
0069 clusters_v_.shrink_to_fit();
0070 layerClustersPerLayer_.clear();
0071 layerClustersPerLayer_.shrink_to_fit();
0072 for (auto &it : points_) {
0073 it.clear();
0074 it.shrink_to_fit();
0075 std::vector<KDNode>().swap(it);
0076 }
0077 for (unsigned int i = 0; i < minpos_.size(); i++) {
0078 minpos_[i][0] = 0.;
0079 minpos_[i][1] = 0.;
0080 maxpos_[i][0] = 0.;
0081 maxpos_[i][1] = 0.;
0082 }
0083 }
0084 void computeThreshold();
0085
0086
0087 Density getDensity() override;
0088
0089 static void fillPSetDescription(edm::ParameterSetDescription &iDesc) {
0090 iDesc.add<std::vector<double>>("thresholdW0", {2.9, 2.9, 2.9});
0091 iDesc.add<std::vector<double>>("positionDeltaRho_c", {1.3, 1.3, 1.3});
0092 iDesc.add<std::vector<double>>("deltac",
0093 {
0094 2.0,
0095 2.0,
0096 5.0,
0097 });
0098 iDesc.add<bool>("dependSensor", true);
0099 iDesc.add<double>("ecut", 3.0);
0100 iDesc.add<double>("kappa", 9.0);
0101 iDesc.addUntracked<unsigned int>("verbosity", 3);
0102 iDesc.add<std::vector<double>>("dEdXweights", {});
0103 iDesc.add<std::vector<double>>("thicknessCorrection", {});
0104 iDesc.add<std::vector<double>>("fcPerMip", {});
0105 iDesc.add<double>("fcPerEle", 0.0);
0106 edm::ParameterSetDescription descNestedNoises;
0107 descNestedNoises.add<std::vector<double>>("values", {});
0108 iDesc.add<edm::ParameterSetDescription>("noises", descNestedNoises);
0109 edm::ParameterSetDescription descNestedNoiseMIP;
0110 descNestedNoiseMIP.add<bool>("scaleByDose", false);
0111 descNestedNoiseMIP.add<double>("scaleByDoseFactor", 1.);
0112 iDesc.add<edm::ParameterSetDescription>("scaleByDose", descNestedNoiseMIP);
0113 descNestedNoiseMIP.add<std::string>("doseMap", "");
0114 iDesc.add<edm::ParameterSetDescription>("doseMap", descNestedNoiseMIP);
0115 descNestedNoiseMIP.add<double>("noise_MIP", 1. / 100.);
0116 iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
0117 }
0118
0119
0120 typedef math::XYZPoint Point;
0121
0122 private:
0123
0124 std::vector<double> thresholdW0_;
0125 std::vector<double> positionDeltaRho_c_;
0126
0127
0128 std::vector<double> vecDeltas_;
0129 double kappa_;
0130
0131
0132 double ecut_;
0133
0134
0135 double sigma2_;
0136
0137
0138 std::vector<reco::BasicCluster> clusters_v_;
0139
0140
0141 Density density_;
0142
0143
0144 bool dependSensor_;
0145 std::vector<double> dEdXweights_;
0146 std::vector<double> thicknessCorrection_;
0147 std::vector<double> fcPerMip_;
0148 double fcPerEle_;
0149 std::vector<double> nonAgedNoises_;
0150 double noiseMip_;
0151 std::vector<std::vector<double>> thresholds_;
0152 std::vector<std::vector<double>> sigmaNoise_;
0153
0154
0155 bool initialized_;
0156
0157 struct Hexel {
0158 double x;
0159 double y;
0160 double z;
0161 bool isHalfCell;
0162 double weight;
0163 double fraction;
0164 DetId detid;
0165 double rho;
0166 double delta;
0167 int nearestHigher;
0168 bool isBorder;
0169 bool isHalo;
0170 int clusterIndex;
0171 float sigmaNoise;
0172 float thickness;
0173 const hgcal::RecHitTools *tools;
0174
0175 Hexel(const HGCRecHit &hit,
0176 DetId id_in,
0177 bool isHalf,
0178 float sigmaNoise_in,
0179 float thickness_in,
0180 const hgcal::RecHitTools *tools_in)
0181 : isHalfCell(isHalf),
0182 weight(0.),
0183 fraction(1.0),
0184 detid(id_in),
0185 rho(0.),
0186 delta(0.),
0187 nearestHigher(-1),
0188 isBorder(false),
0189 isHalo(false),
0190 clusterIndex(-1),
0191 sigmaNoise(sigmaNoise_in),
0192 thickness(thickness_in),
0193 tools(tools_in) {
0194 const GlobalPoint position(tools->getPosition(detid));
0195 weight = hit.energy();
0196 x = position.x();
0197 y = position.y();
0198 z = position.z();
0199 }
0200 Hexel()
0201 : x(0.),
0202 y(0.),
0203 z(0.),
0204 isHalfCell(false),
0205 weight(0.),
0206 fraction(1.0),
0207 detid(),
0208 rho(0.),
0209 delta(0.),
0210 nearestHigher(-1),
0211 isBorder(false),
0212 isHalo(false),
0213 clusterIndex(-1),
0214 sigmaNoise(0.),
0215 thickness(0.),
0216 tools(nullptr) {}
0217 bool operator>(const Hexel &rhs) const { return (rho > rhs.rho); }
0218 };
0219
0220 typedef KDTreeLinkerAlgo<Hexel> KDTree;
0221 typedef KDTreeNodeInfo<Hexel> KDNode;
0222
0223 std::vector<std::vector<std::vector<KDNode>>> layerClustersPerLayer_;
0224
0225 std::vector<size_t> sort_by_delta(const std::vector<KDNode> &v) const {
0226 std::vector<size_t> idx(v.size());
0227 std::iota(std::begin(idx), std::end(idx), 0);
0228 sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) { return v[i1].data.delta > v[i2].data.delta; });
0229 return idx;
0230 }
0231
0232 std::vector<std::vector<KDNode>> points_;
0233
0234
0235 std::vector<std::array<float, 2>> minpos_;
0236 std::vector<std::array<float, 2>> maxpos_;
0237
0238
0239 inline double distance2(const Hexel &pt1, const Hexel &pt2) const {
0240 const double dx = pt1.x - pt2.x;
0241 const double dy = pt1.y - pt2.y;
0242 return (dx * dx + dy * dy);
0243 }
0244 inline double distance(const Hexel &pt1, const Hexel &pt2) const {
0245 return std::sqrt(distance2(pt1, pt2));
0246 }
0247 double calculateLocalDensity(std::vector<KDNode> &, KDTree &, const unsigned int) const;
0248 double calculateDistanceToHigher(std::vector<KDNode> &) const;
0249 int findAndAssignClusters(std::vector<KDNode> &,
0250 KDTree &,
0251 double,
0252 KDTreeBox<2> &,
0253 const unsigned int,
0254 std::vector<std::vector<KDNode>> &) const;
0255 math::XYZPoint calculatePosition(std::vector<KDNode> &) const;
0256
0257
0258 void setDensity(const std::vector<KDNode> &nd);
0259
0260
0261 std::vector<unsigned> findLocalMaximaInCluster(const std::vector<KDNode> &);
0262 math::XYZPoint calculatePositionWithFraction(const std::vector<KDNode> &, const std::vector<double> &);
0263 double calculateEnergyWithFraction(const std::vector<KDNode> &, const std::vector<double> &);
0264
0265 void shareEnergy(const std::vector<KDNode> &, const std::vector<unsigned> &, std::vector<std::vector<double>> &);
0266 };
0267
0268 #endif