HGCalShowerShapeHelper.h

Go to the documentation of this file.

#ifndef RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h
#define RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h

// system include files
#include <map>
#include <memory>
#include <utility>
#include <vector>

// external include files
#include <Math/Vector3Dfwd.h>

// CMSSW include files
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/ForwardDetId/interface/HGCEEDetId.h"
#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
#include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/StreamID.h"
#include "FWCore/Utilities/interface/Transition.h"
#include "Geometry/CaloTopology/interface/HGCalTopology.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"

class HGCalShowerShapeHelper {
  // Good to filter/compute/store this stuff beforehand as they are common to the shower shape variables.
  // No point in filtering, computing layer-wise centroids, etc. for each variable again and again.
  // Once intitialized, one can the calculate different variables one after another for a given object.
  // This is all handled by ShowerShapeCalc class which caches the layer-wise centroids and other
  // heavy variables for an object + set of cuts
  // It was changed to this approach so that we could use this in constant functions

  // In principle should consider the HGCalHSi and HGCalHSc hits (leakage) also.
  // Can have subdetector dependent thresholds and layer selection.
  // To be implemented.

public:
  static const double kLDWaferCellSize_;
  static const double kHDWaferCellSize_;

  struct ShowerWidths {
    double sigma2xx;
    double sigma2yy;
    double sigma2zz;

    double sigma2xy;
    double sigma2yz;
    double sigma2zx;

    double sigma2uu;
    double sigma2vv;
    double sigma2ww;

    ShowerWidths()
        : sigma2xx(0.0),
          sigma2yy(0.0),
          sigma2zz(0.0),
          sigma2xy(0.0),
          sigma2yz(0.0),
          sigma2zx(0.0),
          sigma2uu(0.0),
          sigma2vv(0.0),
          sigma2ww(0.0) {}
  };

  class ShowerShapeCalc {
  public:
    ShowerShapeCalc(std::shared_ptr<const hgcal::RecHitTools> recHitTools,
                    std::shared_ptr<const std::unordered_map<uint32_t, const reco::PFRecHit *> > pfRecHitPtrMap,
                    const std::vector<std::pair<DetId, float> > &hitsAndFracs,
                    const double rawEnergy,
                    const double minHitE = 0,
                    const double minHitET = 0,
                    const int minLayer = 1,
                    const int maxLayer = -1,
                    DetId::Detector subDet = DetId::HGCalEE);

    double getCellSize(DetId detId) const;

    // Compute Rvar in a cylinder around the layer centroids
    double getRvar(double cylinderR, bool useFractions = true, bool useCellSize = true) const;

    // Compute PCA widths around the layer centroids
    ShowerWidths getPCAWidths(double cylinderR, bool useFractions = false) const;

    std::vector<double> getEnergyHighestHits(unsigned int nrHits, bool useFractions = true) const;

  private:
    void setFilteredHitsAndFractions(const std::vector<std::pair<DetId, float> > &hitsAndFracs);
    void setLayerWiseInfo();

    std::shared_ptr<const hgcal::RecHitTools> recHitTools_;
    std::shared_ptr<const std::unordered_map<uint32_t, const reco::PFRecHit *> > pfRecHitPtrMap_;
    double rawEnergy_;

    double minHitE_;
    double minHitET_;
    double minHitET2_;
    int minLayer_;
    int maxLayer_;
    int nLayer_;
    DetId::Detector subDet_;

    std::vector<std::pair<DetId, float> > hitsAndFracs_;
    std::vector<double> hitEnergies_;
    std::vector<double> hitEnergiesWithFracs_;

    ROOT::Math::XYZVector centroid_;
    std::vector<double> layerEnergies_;
    std::vector<ROOT::Math::XYZVector> layerCentroids_;
  };

  HGCalShowerShapeHelper();
  HGCalShowerShapeHelper(edm::ConsumesCollector &&sumes);
  ~HGCalShowerShapeHelper() = default;
  HGCalShowerShapeHelper(const HGCalShowerShapeHelper &rhs) = delete;
  HGCalShowerShapeHelper(const HGCalShowerShapeHelper &&rhs) = delete;
  HGCalShowerShapeHelper &operator=(const HGCalShowerShapeHelper &rhs) = delete;
  HGCalShowerShapeHelper &operator=(const HGCalShowerShapeHelper &&rhs) = delete;

  template <edm::Transition tr = edm::Transition::Event>
  void setTokens(edm::ConsumesCollector consumesCollector) {
    caloGeometryToken_ = consumesCollector.esConsumes<CaloGeometry, CaloGeometryRecord, tr>();
  }

  void initPerSetup(const edm::EventSetup &iSetup);
  void initPerEvent(const std::vector<reco::PFRecHit> &recHits);
  void initPerEvent(const edm::EventSetup &iSetup, const std::vector<reco::PFRecHit> &recHits);

  HGCalShowerShapeHelper::ShowerShapeCalc createCalc(const std::vector<std::pair<DetId, float> > &hitsAndFracs,
                                                     double rawEnergy,
                                                     double minHitE = 0,
                                                     double minHitET = 0,
                                                     int minLayer = 1,
                                                     int maxLayer = -1,
                                                     DetId::Detector subDet = DetId::HGCalEE) const;
  HGCalShowerShapeHelper::ShowerShapeCalc createCalc(const reco::SuperCluster &sc,
                                                     double minHitE = 0,
                                                     double minHitET = 0,
                                                     int minLayer = 1,
                                                     int maxLayer = -1,
                                                     DetId::Detector subDet = DetId::HGCalEE) const {
    return createCalc(sc.hitsAndFractions(), sc.rawEnergy(), minHitE, minHitET, minLayer, maxLayer, subDet);
  }

private:
  void setPFRecHitPtrMap(const std::vector<reco::PFRecHit> &recHits);

  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
  std::shared_ptr<hgcal::RecHitTools> recHitTools_;
  std::shared_ptr<std::unordered_map<uint32_t, const reco::PFRecHit *> > pfRecHitPtrMap_;
};

#endif