Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-07 23:58:21

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