EgammaTools/interface/EgammaPCAHelper.h

0001 //--------------------------------------------------------------------------------------------------
0002 //
0003 // EGammaPCAHelper
0004 //
0005 // Helper Class to compute PCA
0006 //
0007 //
0008 //--------------------------------------------------------------------------------------------------
0009 #ifndef RecoEgamma_EgammaTools_EGammaPCAHelper_h
0010 #define RecoEgamma_EgammaTools_EGammaPCAHelper_h
0011
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013
0014 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0015 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0016 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0017 #include "DataFormats/Math/interface/Vector3D.h"
0018 #include "DataFormats/ParticleFlowReco/interface/HGCalMultiCluster.h"
0019
0020 #include "RecoEgamma/EgammaTools/interface/Spot.h"
0021 #include "RecoEgamma/EgammaTools/interface/LongDeps.h"
0022 #include "RecoEgamma/EgammaTools/interface/ShowerDepth.h"
0023 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0024 #include "Math/Transform3D.h"
0025 #include <unordered_map>
0026
0027 #include "TPrincipal.h"
0028
0029 class HGCalRecHit;
0030
0031 namespace hgcal {
0032
0033   class EGammaPCAHelper {
0034   public:
0035     typedef ROOT::Math::Transform3D Transform3D;
0036     typedef ROOT::Math::Transform3D::Point Point;
0037
0038     EGammaPCAHelper();
0039
0040     // for the GsfElectrons
0041     void storeRecHits(const reco::CaloCluster &theCluster);
0042     void storeRecHits(const reco::HGCalMultiCluster &cluster);
0043
0044     const TPrincipal &pcaResult();
0045     /// to set once per event
0046     void setHitMap(const std::unordered_map<DetId, const unsigned int> *hitMap);
0047     const std::unordered_map<DetId, const unsigned int> *getHitMap() { return hitMap_; }
0048
0049     void setRecHitTools(const hgcal::RecHitTools *recHitTools);
0050     void setRecHits(edm::Handle<HGCRecHitCollection> recHitHandleEE,
0051                     edm::Handle<HGCRecHitCollection> recHitHandleFH,
0052                     edm::Handle<HGCRecHitCollection> recHitHandleBH);
0053
0054     inline void setdEdXWeights(const std::vector<double> &dEdX) { dEdXWeights_ = dEdX; }
0055
0056     void pcaInitialComputation() { computePCA(-1., false); }
0057
0058     void computePCA(float radius, bool withHalo = true);
0059     const math::XYZPoint &barycenter() const { return barycenter_; }
0060     const math::XYZVector &axis() const { return axis_; }
0061
0062     void computeShowerWidth(float radius, bool withHalo = true);
0063
0064     inline double sigmaUU() const { return checkIteration() ? sigu_ : -1.; }
0065     inline double sigmaVV() const { return checkIteration() ? sigv_ : -1.; }
0066     inline double sigmaEE() const { return checkIteration() ? sige_ : -1.; }
0067     inline double sigmaPP() const { return checkIteration() ? sigp_ : -1.; }
0068
0069     inline const TVectorD &eigenValues() const { return *pca_->GetEigenValues(); }
0070     inline const TVectorD &sigmas() const { return *pca_->GetSigmas(); }
0071     // contains maxlayer+1 values, first layer is [1]
0072     LongDeps energyPerLayer(float radius, bool withHalo = true);
0073
0074     float clusterDepthCompatibility(const LongDeps &, float &measuredDepth, float &expectedDepth, float &expectedSigma);
0075     void printHits(float radius) const;
0076     void clear();
0077
0078   private:
0079     bool checkIteration() const;
0080     void storeRecHits(const std::vector<std::pair<DetId, float>> &hf);
0081     float findZFirstLayer(const LongDeps &) const;
0082
0083     bool recHitsStored_;
0084     bool debug_;
0085
0086     //parameters
0087     std::vector<double> dEdXWeights_;
0088     std::vector<double> invThicknessCorrection_;
0089
0090     const reco::CaloCluster *theCluster_;
0091     const std::unordered_map<DetId, const unsigned int> *hitMap_;
0092     std::vector<Spot> theSpots_;
0093     int pcaIteration_;
0094     unsigned int maxlayer_;
0095
0096     // output quantities
0097     math::XYZPoint barycenter_;
0098     math::XYZVector axis_;
0099
0100     Transform3D trans_;
0101     double sigu_, sigv_, sige_, sigp_;
0102
0103     // helper
0104     std::unique_ptr<TPrincipal> pca_;
0105     const hgcal::RecHitTools *recHitTools_;
0106     ShowerDepth showerDepth_;
0107
0108     std::vector<const HGCRecHit *> hits_;
0109   };
0110
0111 }  // namespace hgcal
0112
0113 #endif