EgammaTools/interface/HGCalShowerShapeHelper.h

0001 #ifndef RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h
0002 #define RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h
0003
0004 // system include files
0005 #include <map>
0006 #include <memory>
0007 #include <utility>
0008 #include <vector>
0009
0010 // external include files
0011 #include <Math/Vector3Dfwd.h>
0012
0013 // CMSSW include files
0014 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0015 #include "DataFormats/ForwardDetId/interface/HGCEEDetId.h"
0016 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0017 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0020 #include "FWCore/Framework/interface/ConsumesCollector.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 #include "FWCore/Framework/interface/stream/EDProducer.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Utilities/interface/StreamID.h"
0026 #include "FWCore/Utilities/interface/Transition.h"
0027 #include "Geometry/CaloTopology/interface/HGCalTopology.h"
0028 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0029 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0030
0031 class HGCalShowerShapeHelper {
0032   // Good to filter/compute/store this stuff beforehand as they are common to the shower shape variables.
0033   // No point in filtering, computing layer-wise centroids, etc. for each variable again and again.
0034   // Once intitialized, one can the calculate different variables one after another for a given object.
0035   // This is all handled by ShowerShapeCalc class which caches the layer-wise centroids and other
0036   // heavy variables for an object + set of cuts
0037   // It was changed to this approach so that we could use this in constant functions
0038
0039   // In principle should consider the HGCalHSi and HGCalHSc hits (leakage) also.
0040   // Can have subdetector dependent thresholds and layer selection.
0041   // To be implemented.
0042
0043 public:
0044   static const double kLDWaferCellSize_;
0045   static const double kHDWaferCellSize_;
0046
0047   struct ShowerWidths {
0048     double sigma2xx;
0049     double sigma2yy;
0050     double sigma2zz;
0051
0052     double sigma2xy;
0053     double sigma2yz;
0054     double sigma2zx;
0055
0056     double sigma2uu;
0057     double sigma2vv;
0058     double sigma2ww;
0059
0060     ShowerWidths()
0061         : sigma2xx(0.0),
0062           sigma2yy(0.0),
0063           sigma2zz(0.0),
0064           sigma2xy(0.0),
0065           sigma2yz(0.0),
0066           sigma2zx(0.0),
0067           sigma2uu(0.0),
0068           sigma2vv(0.0),
0069           sigma2ww(0.0) {}
0070   };
0071
0072   class ShowerShapeCalc {
0073   public:
0074     ShowerShapeCalc(std::shared_ptr<const hgcal::RecHitTools> recHitTools,
0075                     std::shared_ptr<const std::unordered_map<uint32_t, const reco::PFRecHit *> > pfRecHitPtrMap,
0076                     const std::vector<std::pair<DetId, float> > &hitsAndFracs,
0077                     const double rawEnergy,
0078                     const double minHitE = 0,
0079                     const double minHitET = 0,
0080                     const int minLayer = 1,
0081                     const int maxLayer = -1,
0082                     DetId::Detector subDet = DetId::HGCalEE);
0083
0084     double getCellSize(DetId detId) const;
0085
0086     // Compute Rvar in a cylinder around the layer centroids
0087     double getRvar(double cylinderR, bool useFractions = true, bool useCellSize = true) const;
0088
0089     // Compute PCA widths around the layer centroids
0090     ShowerWidths getPCAWidths(double cylinderR, bool useFractions = false) const;
0091
0092     std::vector<double> getEnergyHighestHits(unsigned int nrHits, bool useFractions = true) const;
0093
0094   private:
0095     void setFilteredHitsAndFractions(const std::vector<std::pair<DetId, float> > &hitsAndFracs);
0096     void setLayerWiseInfo();
0097
0098     std::shared_ptr<const hgcal::RecHitTools> recHitTools_;
0099     std::shared_ptr<const std::unordered_map<uint32_t, const reco::PFRecHit *> > pfRecHitPtrMap_;
0100     double rawEnergy_;
0101
0102     double minHitE_;
0103     double minHitET_;
0104     double minHitET2_;
0105     int minLayer_;
0106     int maxLayer_;
0107     int nLayer_;
0108     DetId::Detector subDet_;
0109
0110     std::vector<std::pair<DetId, float> > hitsAndFracs_;
0111     std::vector<double> hitEnergies_;
0112     std::vector<double> hitEnergiesWithFracs_;
0113
0114     ROOT::Math::XYZVector centroid_;
0115     std::vector<double> layerEnergies_;
0116     std::vector<ROOT::Math::XYZVector> layerCentroids_;
0117   };
0118
0119   HGCalShowerShapeHelper();
0120   HGCalShowerShapeHelper(edm::ConsumesCollector &&sumes);
0121   ~HGCalShowerShapeHelper() = default;
0122   HGCalShowerShapeHelper(const HGCalShowerShapeHelper &rhs) = delete;
0123   HGCalShowerShapeHelper(const HGCalShowerShapeHelper &&rhs) = delete;
0124   HGCalShowerShapeHelper &operator=(const HGCalShowerShapeHelper &rhs) = delete;
0125   HGCalShowerShapeHelper &operator=(const HGCalShowerShapeHelper &&rhs) = delete;
0126
0127   template <edm::Transition tr = edm::Transition::Event>
0128   void setTokens(edm::ConsumesCollector consumesCollector) {
0129     caloGeometryToken_ = consumesCollector.esConsumes<CaloGeometry, CaloGeometryRecord, tr>();
0130   }
0131
0132   void initPerSetup(const edm::EventSetup &iSetup);
0133   void initPerEvent(const std::vector<reco::PFRecHit> &recHits);
0134   void initPerEvent(const edm::EventSetup &iSetup, const std::vector<reco::PFRecHit> &recHits);
0135
0136   HGCalShowerShapeHelper::ShowerShapeCalc createCalc(const std::vector<std::pair<DetId, float> > &hitsAndFracs,
0137                                                      double rawEnergy,
0138                                                      double minHitE = 0,
0139                                                      double minHitET = 0,
0140                                                      int minLayer = 1,
0141                                                      int maxLayer = -1,
0142                                                      DetId::Detector subDet = DetId::HGCalEE) const;
0143   HGCalShowerShapeHelper::ShowerShapeCalc createCalc(const reco::SuperCluster &sc,
0144                                                      double minHitE = 0,
0145                                                      double minHitET = 0,
0146                                                      int minLayer = 1,
0147                                                      int maxLayer = -1,
0148                                                      DetId::Detector subDet = DetId::HGCalEE) const {
0149     return createCalc(sc.hitsAndFractions(), sc.rawEnergy(), minHitE, minHitET, minLayer, maxLayer, subDet);
0150   }
0151
0152 private:
0153   void setPFRecHitPtrMap(const std::vector<reco::PFRecHit> &recHits);
0154
0155   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0156   std::shared_ptr<hgcal::RecHitTools> recHitTools_;
0157   std::shared_ptr<std::unordered_map<uint32_t, const reco::PFRecHit *> > pfRecHitPtrMap_;
0158 };
0159
0160 #endif