Adding shower shape variables

RecoEgamma/EgammaHLTProducers/plugins/BuildFile.xml

@@ -4,6 +4,7 @@
 <use name="Geometry/CaloGeometry"/>
 <use name="RecoEcal/EgammaCoreTools"/>
 <use name="RecoEgamma/EgammaIsolationAlgos"/>
+<use name="RecoEgamma/EgammaTools"/>
 <use name="MagneticField/Engine"/>
 <use name="MagneticField/Records"/>
 <use name="DataFormats/DetId"/>

RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTHGCalIDVarProducer.cc

@@ -0,0 +1,139 @@
+
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
+#include "DataFormats/EgammaReco/interface/SuperCluster.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
+
+#include "DataFormats/EgammaReco/interface/SuperCluster.h"
+#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
+
+#include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
+#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
+
+#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
+
+#include "RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h"
+#include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h"
+
+class EgammaHLTHGCalIDVarProducer : public edm::stream::EDProducer<> {
+public:
+  explicit EgammaHLTHGCalIDVarProducer(const edm::ParameterSet&);
+  ~EgammaHLTHGCalIDVarProducer() override;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+  void produce(edm::Event&, const edm::EventSetup&) override;
+
+  class PCAAssocMap {
+  public:
+    PCAAssocMap(double HGCalShowerShapeHelper::ShowerWidths::*var, const std::string& name) : var_(var), name_(name) {}
+
+    std::unique_ptr<reco::RecoEcalCandidateIsolationMap>& initMap(
+        const edm::Handle<reco::RecoEcalCandidateCollection>& candHandle) {
+      assocMap_ = std::make_unique<reco::RecoEcalCandidateIsolationMap>(candHandle);
+      return assocMap_;
+    }
+
+    void insert(reco::RecoEcalCandidateRef& ref, const HGCalShowerShapeHelper::ShowerWidths& showerWidths) {
+      assocMap_->insert(ref, showerWidths.*var_);
+    }
+
+    std::unique_ptr<reco::RecoEcalCandidateIsolationMap>& operator()() { return assocMap_; }
+    const std::string& name() const { return name_; }
+
+  private:
+    double HGCalShowerShapeHelper::ShowerWidths::*var_;
+    std::string name_;
+    std::unique_ptr<reco::RecoEcalCandidateIsolationMap> assocMap_;
+  };
+
+private:
+  // ----------member data ---------------------------
+  float rCylinder_;
+  float hOverECone_;
+  std::vector<PCAAssocMap> pcaAssocMaps_;
+  const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandidateToken_;
+  const edm::EDGetTokenT<reco::PFRecHitCollection> hgcalRecHitToken_;
+  const edm::EDGetTokenT<reco::CaloClusterCollection> layerClusterToken_;
+};
+
+EgammaHLTHGCalIDVarProducer::EgammaHLTHGCalIDVarProducer(const edm::ParameterSet& config)
+    : rCylinder_(config.getParameter<double>("rCylinder")),
+      hOverECone_(config.getParameter<double>("hOverECone")),
+      recoEcalCandidateToken_(
+          consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
+      hgcalRecHitToken_(consumes<reco::PFRecHitCollection>(config.getParameter<edm::InputTag>("hgcalRecHits"))),
+      layerClusterToken_(consumes<reco::CaloClusterCollection>(config.getParameter<edm::InputTag>("layerClusters"))) {
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2xx, "sigma2xx"));
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2yy, "sigma2yy"));
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2zz, "sigma2zz"));
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2xy, "sigma2xy"));
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2yz, "sigma2yz"));
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2zx, "sigma2zx"));
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2uu, "sigma2uu"));
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2vv, "sigma2vv"));
+  pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2ww, "sigma2ww"));
+
+  produces<reco::RecoEcalCandidateIsolationMap>("rVar");
+  produces<reco::RecoEcalCandidateIsolationMap>("hForHOverE");
+  for (auto& var : pcaAssocMaps_) {
+    produces<reco::RecoEcalCandidateIsolationMap>(var.name());
+  }
+}
+
+EgammaHLTHGCalIDVarProducer::~EgammaHLTHGCalIDVarProducer() {}
+
+void EgammaHLTHGCalIDVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("recoEcalCandidateProducer", edm::InputTag("hltL1SeededRecoEcalCandidate"));
+  desc.add<edm::InputTag>("hgcalRecHits", edm::InputTag("hgcalRecHits"));
+  desc.add<edm::InputTag>("layerClusters", edm::InputTag("layerClusters"));
+  desc.add<double>("rCylinder", 2.8);
+  desc.add<double>("hOverECone", 0.15);
+  descriptions.add(("hltEgammaHLTHGCalIDVarProducer"), desc);
+}
+
+void EgammaHLTHGCalIDVarProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandidateToken_);
+  auto hgcalRecHitHandle = iEvent.getHandle(hgcalRecHitToken_);
+  auto layerClustersHandle = iEvent.getHandle(layerClusterToken_);
+
+  HGCalShowerShapeHelper ssCalc;
+  ssCalc.initPerEvent(iSetup, *hgcalRecHitHandle);
+
+  auto rVarMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
+  auto hForHoverEMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
+  for (auto& pcaMap : pcaAssocMaps_) {
+    pcaMap.initMap(recoEcalCandHandle);
+  }
+
+  for (size_t candNr = 0; candNr < recoEcalCandHandle->size(); candNr++) {
+    reco::RecoEcalCandidateRef candRef(recoEcalCandHandle, candNr);
+    ssCalc.initPerObject(candRef->superCluster()->hitsAndFractions());
+    rVarMap->insert(candRef, ssCalc.getRvar(rCylinder_, candRef->superCluster()->energy()));
+
+    float hForHoverE = HGCalClusterTools::hadEnergyInCone(
+        candRef->superCluster()->eta(), candRef->superCluster()->phi(), *layerClustersHandle, 0., hOverECone_, 0., 0.);
+    hForHoverEMap->insert(candRef, hForHoverE);
+    auto pcaWidths = ssCalc.getPCAWidths(rCylinder_);
+    for (auto& pcaMap : pcaAssocMaps_) {
+      pcaMap.insert(candRef, pcaWidths);
+    }
+  }
+  iEvent.put(std::move(rVarMap), "rVar");
+  iEvent.put(std::move(hForHoverEMap), "hForHOverE");
+  for (auto& pcaMap : pcaAssocMaps_) {
+    iEvent.put(std::move(pcaMap()), pcaMap.name());
+  }
+}
+
+DEFINE_FWK_MODULE(EgammaHLTHGCalIDVarProducer);

RecoEgamma/EgammaTools/interface/HGCalClusterTools.h

@@ -0,0 +1,44 @@
+#ifndef RecoEgamma_EgammaTools_HGCalClusterTools_h
+#define RecoEgamma_EgammaTools_HGCalClusterTools_h
+
+#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
+#include <vector>
+
+class HGCalClusterTools {
+public:
+  enum class EType { ET, ENERGY };
+
+  static float energyInCone(const float eta,
+                            const float phi,
+                            const std::vector<reco::CaloCluster>& layerClusters,
+                            const float minDR,
+                            const float maxDR,
+                            const float minEt,
+                            const float minEnergy,
+                            const std::vector<DetId::Detector>& subDets,
+                            const HGCalClusterTools::EType& eType = EType::ENERGY);
+
+  static float hadEnergyInCone(const float eta,
+                               const float phi,
+                               const std::vector<reco::CaloCluster>& layerClusters,
+                               const float minDR,
+                               const float maxDR,
+                               const float minEt,
+                               const float minEnergy,
+                               const HGCalClusterTools::EType& eType = EType::ENERGY) {
+    return energyInCone(
+        eta, phi, layerClusters, minDR, maxDR, minEt, minEnergy, {DetId::HGCalHSi, DetId::HGCalHSc}, eType);
+  }
+  static float emEnergyInCone(const float eta,
+                              const float phi,
+                              const std::vector<reco::CaloCluster>& layerClusters,
+                              const float minDR,
+                              const float maxDR,
+                              const float minEt,
+                              const float minEnergy,
+                              const HGCalClusterTools::EType& eType = EType::ENERGY) {
+    return energyInCone(eta, phi, layerClusters, minDR, maxDR, minEt, minEnergy, {DetId::HGCalEE}, eType);
+  }
+};
+
+#endif

RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h

@@ -0,0 +1,131 @@
+#ifndef RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h
+#define RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h
+
+// system include files
+#include <memory>
+
+// user include files
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
+#include "DataFormats/FWLite/interface/ESHandle.h"
+#include "DataFormats/ForwardDetId/interface/HGCEEDetId.h"
+#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
+#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
+#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.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 "Geometry/CaloTopology/interface/HGCalTopology.h"
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+#include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
+
+#include <CLHEP/Vector/LorentzVector.h>
+
+#include <algorithm>
+#include <cstdlib>
+#include <iostream>
+#include <map>
+#include <string>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include <TMatrixD.h>
+#include <TVectorD.h>
+
+#include <Math/Point3D.h>
+#include <Math/Point3Dfwd.h>
+
+class HGCalShowerShapeHelper {
+private:
+  // 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.
+  // If a different set of preselections (E, ET, etc.) is required for a given object, then reinitialize using initPerObject(...).
+
+  // In principle should consider the HGCalHSi and HGCalHSc hits (leakage) also.
+  // Can have subdetector dependent thresholds and layer selection.
+  // To be implemented.
+
+  double minHitE_;
+  double minHitET_;
+  double minHitET2_;
+  int minLayer_;
+  int maxLayer_;
+  int nLayer_;
+  DetId::Detector subDet_;
+
+  hgcal::RecHitTools recHitTools_;
+
+  std::unordered_map<uint32_t, const reco::PFRecHit *> pfRecHitPtrMap_;
+  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_;
+
+  void setPFRecHitPtrMap(const std::vector<reco::PFRecHit> &recHits);
+
+  void setFilteredHitsAndFractions(const std::vector<std::pair<DetId, float> > &hitsAndFracs);
+
+public:
+  static constexpr double kLDWaferCellSize_ = 0.698;
+  static constexpr double kHDWaferCellSize_ = 0.465;
+
+  void setLayerWiseStuff();
+
+  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) {}
+  };
+
+  void initPerEvent(const edm::EventSetup &iSetup, const std::vector<reco::PFRecHit> &recHits);
+
+  void initPerObject(const std::vector<std::pair<DetId, float> > &hitsAndFracs,
+                     double minHitE = 0,
+                     double minHitET = 0,
+                     int minLayer = 1,
+                     int maxLayer = 28,
+                     DetId::Detector subDet = DetId::HGCalEE);
+
+  const double getCellSize(DetId detId);
+
+  // Compute Rvar in a cylinder around the layer centroids
+  const double getRvar(double cylinderR, double energyNorm, bool useFractions = true, bool useCellSize = true);
+
+  // Compute PCA widths around the layer centroids
+  const ShowerWidths getPCAWidths(double cylinderR, bool useFractions = false);
+};
+
+#endif

RecoEgamma/EgammaTools/src/HGCalClusterTools.cc

@@ -0,0 +1,56 @@
+#include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h"
+
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "DataFormats/Math/interface/deltaR.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+
+float HGCalClusterTools::energyInCone(const float eta,
+                                      const float phi,
+                                      const std::vector<reco::CaloCluster>& layerClusters,
+                                      const float minDR,
+                                      const float maxDR,
+                                      const float minEt,
+                                      const float minEnergy,
+                                      const std::vector<DetId::Detector>& subDets,
+                                      const HGCalClusterTools::EType& eType) {
+  float hadValue = 0.;
+
+  const float minDR2 = minDR * minDR;
+  const float maxDR2 = maxDR * maxDR;
+
+  for (auto& clus : layerClusters) {
+    if (clus.energy() < minEnergy) {
+      continue;
+    }
+
+    if (std::find(subDets.begin(), subDets.end(), clus.seed().det()) == subDets.end()) {
+      continue;
+    }
+
+    float clusEt = clus.energy() * std::sin(clus.position().theta());
+    if (clusEt < minEt) {
+      continue;
+    }
+
+    //this is a prefilter on the clusters before we calculuate
+    //the expensive eta() of the cluster
+    float dPhi = reco::deltaPhi(phi, clus.phi());
+    if (dPhi > maxDR) {
+      continue;
+    }
+
+    float dR2 = reco::deltaR2(eta, phi, clus.eta(), clus.phi());
+    if (dR2 < minDR2 || dR2 > maxDR2) {
+      continue;
+    }
+    switch (eType) {
+      case EType::ET:
+        hadValue += clusEt;
+        break;
+      case EType::ENERGY:
+        hadValue += clus.energy();
+        break;
+    }
+  }
+  return hadValue;
+}