EgammaPCAHelper.h

Go to the documentation of this file.

//--------------------------------------------------------------------------------------------------
//
// EGammaPCAHelper
//
// Helper Class to compute PCA
//
//
//--------------------------------------------------------------------------------------------------
#ifndef RecoEgamma_EgammaTools_EGammaPCAHelper_h
#define RecoEgamma_EgammaTools_EGammaPCAHelper_h
 
#include "FWCore/Framework/interface/Frameworkfwd.h"
 
#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
#include "DataFormats/Math/interface/Vector3D.h"
#include "DataFormats/ParticleFlowReco/interface/HGCalMultiCluster.h"
 
#include "RecoEgamma/EgammaTools/interface/Spot.h"
#include "RecoEgamma/EgammaTools/interface/LongDeps.h"
#include "RecoEgamma/EgammaTools/interface/ShowerDepth.h"
#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
#include "Math/Transform3D.h"
#include <map>
 
#include "TPrincipal.h"
 
class HGCalRecHit;
 
namespace hgcal {
 
  class EGammaPCAHelper {
  public:
    typedef ROOT::Math::Transform3D Transform3D;
    typedef ROOT::Math::Transform3D::Point Point;
 
    EGammaPCAHelper();
    ~EGammaPCAHelper();
 
    // for the GsfElectrons
    void storeRecHits(const reco::CaloCluster &theCluster);
    void storeRecHits(const reco::HGCalMultiCluster &cluster);
 
    const TPrincipal &pcaResult();
    void setHitMap(std::map<DetId, const HGCRecHit *> *hitMap);
    void fillHitMap(const HGCRecHitCollection &HGCEERecHits,
                    const HGCRecHitCollection &HGCFHRecHits,
                    const HGCRecHitCollection &HGCBHRecHits);
 
    std::map<DetId, const HGCRecHit *> *getHitMap() { return hitMap_; }
 
    void setRecHitTools(const hgcal::RecHitTools *recHitTools);
 
    inline void setdEdXWeights(const std::vector<double> &dEdX) { dEdXWeights_ = dEdX; }
 
    void pcaInitialComputation() { computePCA(-1., false); }
 
    void computePCA(float radius, bool withHalo = true);
    const math::XYZPoint &barycenter() const { return barycenter_; }
    const math::XYZVector &axis() const { return axis_; }
 
    void computeShowerWidth(float radius, bool withHalo = true);
 
    inline double sigmaUU() const { return checkIteration() ? sigu_ : -1.; }
    inline double sigmaVV() const { return checkIteration() ? sigv_ : -1.; }
    inline double sigmaEE() const { return checkIteration() ? sige_ : -1.; }
    inline double sigmaPP() const { return checkIteration() ? sigp_ : -1.; }
 
    inline const TVectorD &eigenValues() const { return *pca_->GetEigenValues(); }
    inline const TVectorD &sigmas() const { return *pca_->GetSigmas(); }
    // contains maxlayer+1 values, first layer is [1]
    LongDeps energyPerLayer(float radius, bool withHalo = true);
 
    float clusterDepthCompatibility(const LongDeps &, float &measuredDepth, float &expectedDepth, float &expectedSigma);
    void printHits(float radius) const;
    void clear();
 
  private:
    bool checkIteration() const;
    void storeRecHits(const std::vector<std::pair<DetId, float>> &hf);
    float findZFirstLayer(const LongDeps &) const;
 
    bool recHitsStored_;
    bool debug_;
 
    //parameters
    std::vector<double> dEdXWeights_;
    std::vector<double> invThicknessCorrection_;
 
    int hitMapOrigin_;  // 0 not initialized; 1 set from outside ; 2 set from inside
    const reco::CaloCluster *theCluster_;
    std::map<DetId, const HGCRecHit *> *hitMap_;
    std::vector<Spot> theSpots_;
    int pcaIteration_;
    unsigned int maxlayer_;
 
    // output quantities
    math::XYZPoint barycenter_;
    math::XYZVector axis_;
 
    Transform3D trans_;
    double sigu_, sigv_, sige_, sigp_;
 
    // helper
    std::unique_ptr<TPrincipal> pca_;
    const hgcal::RecHitTools *recHitTools_;
    ShowerDepth showerDepth_;
  };
 
}  // namespace hgcal
 
#endif