Back to home page

Project CMSSW displayed by LXR
 
 

    

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 // C/C++ headers
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   // this is the method that will start the clusterisation (it is possible to invoke this method more than once - but make sure it is with
0059   // different hit collections (or else use reset)
0060 
0061   void makeClusters() override;
0062 
0063   // this is the method to get the cluster collection out
0064   std::vector<reco::BasicCluster> getClusters(bool) override;
0065 
0066   // use this if you want to reuse the same cluster object but don't want to accumulate clusters (hardly useful?)
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   //getDensity
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   /// point in the space
0120   typedef math::XYZPoint Point;
0121 
0122 private:
0123   // To compute the cluster position
0124   std::vector<double> thresholdW0_;
0125   std::vector<double> positionDeltaRho_c_;
0126 
0127   // The two parameters used to identify clusters
0128   std::vector<double> vecDeltas_;
0129   double kappa_;
0130 
0131   // The hit energy cutoff
0132   double ecut_;
0133 
0134   // for energy sharing
0135   double sigma2_;  // transverse shower size
0136 
0137   // The vector of clusters
0138   std::vector<reco::BasicCluster> clusters_v_;
0139 
0140   // For keeping the density per hit
0141   Density density_;
0142 
0143   // various parameters used for calculating the noise levels for a given sensor (and whether to use them)
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   // initialization bool
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_;  //a vector of vectors of hexels, one for each layer
0233   //@@EM todo: the number of layers should be obtained programmatically - the range is 1-n instead of 0-n-1...
0234 
0235   std::vector<std::array<float, 2>> minpos_;
0236   std::vector<std::array<float, 2>> maxpos_;
0237 
0238   //these functions should be in a helper class.
0239   inline double distance2(const Hexel &pt1, const Hexel &pt2) const {  //distance squared
0240     const double dx = pt1.x - pt2.x;
0241     const double dy = pt1.y - pt2.y;
0242     return (dx * dx + dy * dy);
0243   }                                                                   //distance squaredq
0244   inline double distance(const Hexel &pt1, const Hexel &pt2) const {  //2-d distance on the layer (x-y)
0245     return std::sqrt(distance2(pt1, pt2));
0246   }
0247   double calculateLocalDensity(std::vector<KDNode> &, KDTree &, const unsigned int) const;  //return max density
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   //For keeping the density information
0258   void setDensity(const std::vector<KDNode> &nd);
0259 
0260   // attempt to find subclusters within a given set of hexels
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   // outputs
0265   void shareEnergy(const std::vector<KDNode> &, const std::vector<unsigned> &, std::vector<std::vector<double>> &);
0266 };
0267 
0268 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!