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 <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   // this is the method that will start the clusterisation (it is possible to invoke this method
0058   // more than once - but make sure it is with different hit collections (or else use reset)
0059
0060   void makeClusters() override;
0061
0062   // this is the method to get the cluster collection out
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,  // for scintillator
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);  // use 2x2 or 3x3 scenario for scint density calculation
0113   }
0114
0115   /// point in the space
0116   typedef math::XYZPoint Point;
0117
0118 private:
0119   // The two parameters used to identify clusters
0120   std::vector<double> vecDeltas_;
0121   double kappa_;
0122
0123   // The hit energy cutoff
0124   double ecut_;
0125
0126   // various parameters used for calculating the noise levels for a given sensor (and whether to use
0127   // them)
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   // initialization bool
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 {  // 2-d distance on the layer (x-y)
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);  // return max density
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 // explicit template instantiation
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