HGCalRecProducers/plugins/HGCalCLUEAlgo.h

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 // C/C++ headers
0024 #include <string>
0025 #include <vector>
0026 #include <limits>
0027
0028 #define DEBUG_CLUSTERS_ALPAKA 0
0029
0030 #if DEBUG_CLUSTERS_ALPAKA
0031 #include "RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h"
0032 #endif
0033
0034 template <typename TILE, typename STRATEGY>
0035 class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
0036 public:
0037   HGCalCLUEAlgoT(const edm::ParameterSet& ps)
0038       : HGCalClusteringAlgoBase(
0039             (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
0040             reco::CaloCluster::undefined),
0041         vecDeltas_(ps.getParameter<std::vector<double>>("deltac")),
0042         kappa_(ps.getParameter<double>("kappa")),
0043         ecut_(ps.getParameter<double>("ecut")),
0044         dependSensor_(ps.getParameter<bool>("dependSensor")),
0045         dEdXweights_(ps.getParameter<std::vector<double>>("dEdXweights")),
0046         thicknessCorrection_(ps.getParameter<std::vector<double>>("thicknessCorrection")),
0047         sciThicknessCorrection_(ps.getParameter<double>("sciThicknessCorrection")),
0048         deltasi_index_regemfac_(ps.getParameter<int>("deltasi_index_regemfac")),
0049         maxNumberOfThickIndices_(ps.getParameter<unsigned>("maxNumberOfThickIndices")),
0050         fcPerMip_(ps.getParameter<std::vector<double>>("fcPerMip")),
0051         fcPerEle_(ps.getParameter<double>("fcPerEle")),
0052         nonAgedNoises_(ps.getParameter<std::vector<double>>("noises")),
0053         noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
0054         thresholdW0_(ps.getParameter<std::vector<double>>("thresholdW0")),
0055         positionDeltaRho2_(ps.getParameter<double>("positionDeltaRho2")),
0056         use2x2_(ps.getParameter<bool>("use2x2")),
0057         initialized_(false) {
0058 #if DEBUG_CLUSTERS_ALPAKA
0059     moduleType_ = ps.getParameter<std::string>("type");
0060 #endif
0061   }
0062
0063   ~HGCalCLUEAlgoT() override {}
0064
0065   void getEventSetupPerAlgorithm(const edm::EventSetup& es) override;
0066
0067   void populate(const HGCRecHitCollection& hits) override;
0068
0069   // this is the method that will start the clusterisation (it is possible to invoke this method
0070   // more than once - but make sure it is with different hit collections (or else use reset)
0071
0072   void makeClusters() override;
0073
0074   // this is the method to get the cluster collection out
0075   std::vector<reco::BasicCluster> getClusters(bool) override;
0076
0077   void reset() override {
0078     clusters_v_.clear();
0079     clusters_v_.shrink_to_fit();
0080     for (auto& cl : numberOfClustersPerLayer_) {
0081       cl = 0;
0082     }
0083
0084     for (auto& cells : cells_) {
0085       cells.clear();
0086       cells.shrink_to_fit();
0087     }
0088   }
0089
0090   void computeThreshold();
0091
0092   static void fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0093     iDesc.add<std::vector<double>>("thresholdW0", {2.9, 2.9, 2.9});
0094     iDesc.add<double>("positionDeltaRho2", 1.69);
0095     iDesc.add<std::vector<double>>("deltac",
0096                                    {
0097                                        1.3,
0098                                        1.3,
0099                                        1.3,
0100                                        0.0315,  // for scintillator
0101                                    });
0102     iDesc.add<bool>("dependSensor", true);
0103     iDesc.add<double>("ecut", 3.0);
0104     iDesc.add<double>("kappa", 9.0);
0105     iDesc.addUntracked<unsigned int>("verbosity", 3);
0106     iDesc.add<std::vector<double>>("dEdXweights", {});
0107     iDesc.add<std::vector<double>>("thicknessCorrection", {});
0108     iDesc.add<double>("sciThicknessCorrection", 0.9);
0109     iDesc.add<int>("deltasi_index_regemfac", 3);
0110     iDesc.add<unsigned>("maxNumberOfThickIndices", 6);
0111     iDesc.add<std::vector<double>>("fcPerMip", {});
0112     iDesc.add<double>("fcPerEle", 0.0);
0113     iDesc.add<std::vector<double>>("noises", {});
0114     edm::ParameterSetDescription descNestedNoiseMIP;
0115     descNestedNoiseMIP.add<bool>("scaleByDose", false);
0116     descNestedNoiseMIP.add<unsigned int>("scaleByDoseAlgo", 0);
0117     descNestedNoiseMIP.add<double>("scaleByDoseFactor", 1.);
0118     descNestedNoiseMIP.add<std::string>("doseMap", "");
0119     descNestedNoiseMIP.add<std::string>("sipmMap", "");
0120     descNestedNoiseMIP.add<double>("referenceIdark", -1);
0121     descNestedNoiseMIP.add<double>("referenceXtalk", -1);
0122     descNestedNoiseMIP.add<double>("noise_MIP", 1. / 100.);
0123     iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
0124     iDesc.add<bool>("use2x2", true);  // use 2x2 or 3x3 scenario for scint density calculation
0125   }
0126
0127   /// point in the space
0128   typedef math::XYZPoint Point;
0129
0130 private:
0131   // The two parameters used to identify clusters
0132   std::vector<double> vecDeltas_;
0133   double kappa_;
0134
0135   // The hit energy cutoff
0136   double ecut_;
0137
0138   // various parameters used for calculating the noise levels for a given sensor (and whether to use
0139   // them)
0140   bool dependSensor_;
0141   std::vector<double> dEdXweights_;
0142   std::vector<double> thicknessCorrection_;
0143   double sciThicknessCorrection_;
0144   int deltasi_index_regemfac_;
0145   unsigned maxNumberOfThickIndices_;
0146   std::vector<double> fcPerMip_;
0147   double fcPerEle_;
0148   std::vector<double> nonAgedNoises_;
0149   double noiseMip_;
0150   std::vector<std::vector<double>> thresholds_;
0151   std::vector<std::vector<double>> v_sigmaNoise_;
0152   std::vector<double> thresholdW0_;
0153   double positionDeltaRho2_;
0154
0155   bool use2x2_;
0156
0157   // initialization bool
0158   bool initialized_;
0159
0160   float outlierDeltaFactor_ = 2.f;
0161
0162   struct CellsOnLayer {
0163     std::vector<DetId> detid;
0164     std::vector<float> dim1;
0165     std::vector<float> dim2;
0166
0167     std::vector<float> weight;
0168     std::vector<float> rho;
0169
0170     std::vector<float> delta;
0171     std::vector<int> nearestHigher;
0172     std::vector<int> clusterIndex;
0173     std::vector<float> sigmaNoise;
0174     std::vector<std::vector<int>> followers;
0175     std::vector<bool> isSeed;
0176     float layerDim3 = std::numeric_limits<float>::infinity();
0177
0178     void clear() {
0179       detid.clear();
0180       dim1.clear();
0181       dim2.clear();
0182       weight.clear();
0183       rho.clear();
0184       delta.clear();
0185       nearestHigher.clear();
0186       clusterIndex.clear();
0187       sigmaNoise.clear();
0188       followers.clear();
0189       isSeed.clear();
0190     }
0191
0192     void shrink_to_fit() {
0193       detid.shrink_to_fit();
0194       dim1.shrink_to_fit();
0195       dim2.shrink_to_fit();
0196       weight.shrink_to_fit();
0197       rho.shrink_to_fit();
0198       delta.shrink_to_fit();
0199       nearestHigher.shrink_to_fit();
0200       clusterIndex.shrink_to_fit();
0201       sigmaNoise.shrink_to_fit();
0202       followers.shrink_to_fit();
0203       isSeed.shrink_to_fit();
0204     }
0205   };
0206
0207   std::vector<CellsOnLayer> cells_;
0208
0209   std::vector<int> numberOfClustersPerLayer_;
0210
0211 #if DEBUG_CLUSTERS_ALPAKA
0212   std::string moduleType_;
0213 #endif
0214
0215   inline float distance2(const TILE& lt, int cell1, int cell2, int layerId) const {  // 2-d distance on the layer (x-y)
0216     return (lt.distance2(cells_[layerId].dim1[cell1],
0217                          cells_[layerId].dim2[cell1],
0218                          cells_[layerId].dim1[cell2],
0219                          cells_[layerId].dim2[cell2]));
0220   }
0221
0222   inline float distance(const TILE& lt, int cell1, int cell2, int layerId) const {  // 2-d distance on the layer (x-y)
0223     return std::sqrt(lt.distance2(cells_[layerId].dim1[cell1],
0224                                   cells_[layerId].dim2[cell1],
0225                                   cells_[layerId].dim1[cell2],
0226                                   cells_[layerId].dim2[cell2]));
0227   }
0228
0229   void prepareDataStructures(const unsigned int layerId);
0230   void calculateLocalDensity(const TILE& lt, const unsigned int layerId,
0231                              float delta);  // return max density
0232   void calculateLocalDensity(const TILE& lt, const unsigned int layerId, float delta, HGCalSiliconStrategy strategy);
0233   void calculateLocalDensity(const TILE& lt,
0234                              const unsigned int layerId,
0235                              float delta,
0236                              HGCalScintillatorStrategy strategy);
0237   void calculateDistanceToHigher(const TILE& lt, const unsigned int layerId, float delta);
0238   int findAndAssignClusters(const unsigned int layerId, float delta);
0239 };
0240
0241 // explicit template instantiation
0242 extern template class HGCalCLUEAlgoT<HGCalSiliconLayerTiles, HGCalSiliconStrategy>;
0243 extern template class HGCalCLUEAlgoT<HGCalScintillatorLayerTiles, HGCalScintillatorStrategy>;
0244 extern template class HGCalCLUEAlgoT<HFNoseLayerTiles, HGCalSiliconStrategy>;
0245
0246 using HGCalSiCLUEAlgo = HGCalCLUEAlgoT<HGCalSiliconLayerTiles, HGCalSiliconStrategy>;
0247 using HGCalSciCLUEAlgo = HGCalCLUEAlgoT<HGCalScintillatorLayerTiles, HGCalScintillatorStrategy>;
0248 using HFNoseCLUEAlgo = HGCalCLUEAlgoT<HFNoseLayerTiles, HGCalSiliconStrategy>;
0249
0250 #endif