Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:37 

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)
0034       : HGCalClusteringAlgoBase(
0035             (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
0036             reco::CaloCluster::undefined),
0037         thresholdW0_(ps.getParameter<std::vector<double>>("thresholdW0")),
0038         positionDeltaRho_c_(ps.getParameter<std::vector<double>>("positionDeltaRho_c")),
0039         vecDeltas_(ps.getParameter<std::vector<double>>("deltac")),
0040         kappa_(ps.getParameter<double>("kappa")),
0041         ecut_(ps.getParameter<double>("ecut")),
0042         sigma2_(1.0),
0043         dependSensor_(ps.getParameter<bool>("dependSensor")),
0044         dEdXweights_(ps.getParameter<std::vector<double>>("dEdXweights")),
0045         thicknessCorrection_(ps.getParameter<std::vector<double>>("thicknessCorrection")),
0046         fcPerMip_(ps.getParameter<std::vector<double>>("fcPerMip")),
0047         fcPerEle_(ps.getParameter<double>("fcPerEle")),
0048         nonAgedNoises_(ps.getParameter<edm::ParameterSet>("noises").getParameter<std::vector<double>>("values")),
0049         noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
0050         initialized_(false) {}
0051 
0052   ~HGCalImagingAlgo() override {}
0053 
0054   void getEventSetupPerAlgorithm(const edm::EventSetup &es) override;
0055 
0056   void populate(const HGCRecHitCollection &hits) override;
0057   // 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
0058   // 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   // use this if you want to reuse the same cluster object but don't want to accumulate clusters (hardly useful?)
0066   void reset() override {
0067     clusters_v_.clear();
0068     clusters_v_.shrink_to_fit();
0069     layerClustersPerLayer_.clear();
0070     layerClustersPerLayer_.shrink_to_fit();
0071     for (auto &it : points_) {
0072       it.clear();
0073       it.shrink_to_fit();
0074       std::vector<KDNode>().swap(it);
0075     }
0076     for (unsigned int i = 0; i < minpos_.size(); i++) {
0077       minpos_[i][0] = 0.;
0078       minpos_[i][1] = 0.;
0079       maxpos_[i][0] = 0.;
0080       maxpos_[i][1] = 0.;
0081     }
0082   }
0083   void computeThreshold();
0084 
0085   //getDensity
0086   Density getDensity() override;
0087 
0088   static void fillPSetDescription(edm::ParameterSetDescription &iDesc) {
0089     iDesc.add<std::vector<double>>("thresholdW0", {2.9, 2.9, 2.9});
0090     iDesc.add<std::vector<double>>("positionDeltaRho_c", {1.3, 1.3, 1.3});
0091     iDesc.add<std::vector<double>>("deltac",
0092                                    {
0093                                        2.0,
0094                                        2.0,
0095                                        5.0,
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<std::vector<double>>("fcPerMip", {});
0104     iDesc.add<double>("fcPerEle", 0.0);
0105     edm::ParameterSetDescription descNestedNoises;
0106     descNestedNoises.add<std::vector<double>>("values", {});
0107     iDesc.add<edm::ParameterSetDescription>("noises", descNestedNoises);
0108     edm::ParameterSetDescription descNestedNoiseMIP;
0109     descNestedNoiseMIP.add<bool>("scaleByDose", false);
0110     descNestedNoiseMIP.add<double>("scaleByDoseFactor", 1.);
0111     iDesc.add<edm::ParameterSetDescription>("scaleByDose", descNestedNoiseMIP);
0112     descNestedNoiseMIP.add<std::string>("doseMap", "");
0113     iDesc.add<edm::ParameterSetDescription>("doseMap", descNestedNoiseMIP);
0114     descNestedNoiseMIP.add<double>("noise_MIP", 1. / 100.);
0115     iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
0116   }
0117 
0118   /// point in the space
0119   typedef math::XYZPoint Point;
0120 
0121 private:
0122   // To compute the cluster position
0123   std::vector<double> thresholdW0_;
0124   std::vector<double> positionDeltaRho_c_;
0125 
0126   // The two parameters used to identify clusters
0127   std::vector<double> vecDeltas_;
0128   double kappa_;
0129 
0130   // The hit energy cutoff
0131   double ecut_;
0132 
0133   // for energy sharing
0134   double sigma2_;  // transverse shower size
0135 
0136   // The vector of clusters
0137   std::vector<reco::BasicCluster> clusters_v_;
0138 
0139   // For keeping the density per hit
0140   Density density_;
0141 
0142   // various parameters used for calculating the noise levels for a given sensor (and whether to use them)
0143   bool dependSensor_;
0144   std::vector<double> dEdXweights_;
0145   std::vector<double> thicknessCorrection_;
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>> sigmaNoise_;
0152 
0153   // initialization bool
0154   bool initialized_;
0155 
0156   struct Hexel {
0157     double x;
0158     double y;
0159     double z;
0160     bool isHalfCell;
0161     double weight;
0162     double fraction;
0163     DetId detid;
0164     double rho;
0165     double delta;
0166     int nearestHigher;
0167     bool isBorder;
0168     bool isHalo;
0169     int clusterIndex;
0170     float sigmaNoise;
0171     float thickness;
0172     const hgcal::RecHitTools *tools;
0173 
0174     Hexel(const HGCRecHit &hit,
0175           DetId id_in,
0176           bool isHalf,
0177           float sigmaNoise_in,
0178           float thickness_in,
0179           const hgcal::RecHitTools *tools_in)
0180         : isHalfCell(isHalf),
0181           weight(0.),
0182           fraction(1.0),
0183           detid(id_in),
0184           rho(0.),
0185           delta(0.),
0186           nearestHigher(-1),
0187           isBorder(false),
0188           isHalo(false),
0189           clusterIndex(-1),
0190           sigmaNoise(sigmaNoise_in),
0191           thickness(thickness_in),
0192           tools(tools_in) {
0193       const GlobalPoint position(tools->getPosition(detid));
0194       weight = hit.energy();
0195       x = position.x();
0196       y = position.y();
0197       z = position.z();
0198     }
0199     Hexel()
0200         : x(0.),
0201           y(0.),
0202           z(0.),
0203           isHalfCell(false),
0204           weight(0.),
0205           fraction(1.0),
0206           detid(),
0207           rho(0.),
0208           delta(0.),
0209           nearestHigher(-1),
0210           isBorder(false),
0211           isHalo(false),
0212           clusterIndex(-1),
0213           sigmaNoise(0.),
0214           thickness(0.),
0215           tools(nullptr) {}
0216     bool operator>(const Hexel &rhs) const { return (rho > rhs.rho); }
0217   };
0218 
0219   typedef KDTreeLinkerAlgo<Hexel> KDTree;
0220   typedef KDTreeNodeInfo<Hexel> KDNode;
0221 
0222   std::vector<std::vector<std::vector<KDNode>>> layerClustersPerLayer_;
0223 
0224   std::vector<size_t> sort_by_delta(const std::vector<KDNode> &v) const {
0225     std::vector<size_t> idx(v.size());
0226     std::iota(std::begin(idx), std::end(idx), 0);
0227     sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) { return v[i1].data.delta > v[i2].data.delta; });
0228     return idx;
0229   }
0230 
0231   std::vector<std::vector<KDNode>> points_;  //a vector of vectors of hexels, one for each layer
0232   //@@EM todo: the number of layers should be obtained programmatically - the range is 1-n instead of 0-n-1...
0233 
0234   std::vector<std::array<float, 2>> minpos_;
0235   std::vector<std::array<float, 2>> maxpos_;
0236 
0237   //these functions should be in a helper class.
0238   inline double distance2(const Hexel &pt1, const Hexel &pt2) const {  //distance squared
0239     const double dx = pt1.x - pt2.x;
0240     const double dy = pt1.y - pt2.y;
0241     return (dx * dx + dy * dy);
0242   }  //distance squaredq
0243   inline double distance(const Hexel &pt1, const Hexel &pt2) const {  //2-d distance on the layer (x-y)
0244     return std::sqrt(distance2(pt1, pt2));
0245   }
0246   double calculateLocalDensity(std::vector<KDNode> &, KDTree &, const unsigned int) const;  //return max density
0247   double calculateDistanceToHigher(std::vector<KDNode> &) const;
0248   int findAndAssignClusters(std::vector<KDNode> &,
0249                             KDTree &,
0250                             double,
0251                             KDTreeBox<2> &,
0252                             const unsigned int,
0253                             std::vector<std::vector<KDNode>> &) const;
0254   math::XYZPoint calculatePosition(std::vector<KDNode> &) const;
0255 
0256   //For keeping the density information
0257   void setDensity(const std::vector<KDNode> &nd);
0258 
0259   // attempt to find subclusters within a given set of hexels
0260   std::vector<unsigned> findLocalMaximaInCluster(const std::vector<KDNode> &);
0261   math::XYZPoint calculatePositionWithFraction(const std::vector<KDNode> &, const std::vector<double> &);
0262   double calculateEnergyWithFraction(const std::vector<KDNode> &, const std::vector<double> &);
0263   // outputs
0264   void shareEnergy(const std::vector<KDNode> &, const std::vector<unsigned> &, std::vector<std::vector<double>> &);
0265 };
0266 
0267 #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!