HGCal ID 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,130 @@
+#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) {}
+
+    void initMap(const edm::Handle<reco::RecoEcalCandidateCollection>& candHandle) {
+      assocMap_ = std::make_unique<reco::RecoEcalCandidateIsolationMap>(candHandle);
+    }
+
+    void insert(reco::RecoEcalCandidateRef& ref, const HGCalShowerShapeHelper::ShowerWidths& showerWidths) {
+      assocMap_->insert(ref, showerWidths.*var_);
+    }
+
+    std::unique_ptr<reco::RecoEcalCandidateIsolationMap> releaseMap() { return std::move(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_;
+  HGCalShowerShapeHelper ssCalc_;
+};
+
+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"))),
+      ssCalc_(consumesCollector()) {
+  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_);
+  const auto& hgcalRecHits = iEvent.get(hgcalRecHitToken_);
+  const auto& layerClusters = iEvent.get(layerClusterToken_);
+
+  ssCalc_.initPerEvent(iSetup, hgcalRecHits);
+
+  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(), layerClusters, 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(pcaMap.releaseMap(), pcaMap.name());
+  }
+}
+
+DEFINE_FWK_MODULE(EgammaHLTHGCalIDVarProducer);

RecoEgamma/EgammaHLTProducers/plugins/HLTHGCalLayerClusterIsolationProducer.cc

@@ -0,0 +1,197 @@
+#include <iostream>
+#include <vector>
+#include <memory>
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/Utilities/interface/Exception.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include "DataFormats/Math/interface/deltaR.h"
+
+#include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
+
+#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/RecoEcalCandidate.h"
+#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
+
+#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
+#include "DataFormats/RecoCandidate/interface/RecoChargedCandidateIsolation.h"
+
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
+
+#include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h"
+
+template <typename T1>
+class HLTHGCalLayerClusterIsolationProducer : public edm::stream::EDProducer<> {
+  typedef std::vector<T1> T1Collection;
+  typedef edm::Ref<T1Collection> T1Ref;
+  typedef edm::AssociationMap<edm::OneToValue<std::vector<T1>, float>> T1IsolationMap;
+
+public:
+  explicit HLTHGCalLayerClusterIsolationProducer(const edm::ParameterSet&);
+  ~HLTHGCalLayerClusterIsolationProducer() override = default;
+
+  void produce(edm::Event&, const edm::EventSetup&) override;
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  edm::EDGetTokenT<T1Collection> recoCandidateProducer_;
+  const edm::EDGetTokenT<reco::CaloClusterCollection> layerClusterProducer_;
+  const edm::EDGetTokenT<double> rhoProducer_;
+
+  const double drMax_;
+  const double drVetoEM_;
+  const double drVetoHad_;
+  const double minEnergyEM_;
+  const double minEnergyHad_;
+  const double minEtEM_;
+  const double minEtHad_;
+  const bool useEt_;
+  const bool doRhoCorrection_;
+  const double rhoMax_;
+  const double rhoScale_;
+  const std::vector<double> effectiveAreas_;
+};
+
+template <typename T1>
+HLTHGCalLayerClusterIsolationProducer<T1>::HLTHGCalLayerClusterIsolationProducer(const edm::ParameterSet& config)
+    : layerClusterProducer_(
+          consumes<reco::CaloClusterCollection>(config.getParameter<edm::InputTag>("layerClusterProducer"))),
+      rhoProducer_(consumes<double>(config.getParameter<edm::InputTag>("rhoProducer"))),
+      drMax_(config.getParameter<double>("drMax")),
+      drVetoEM_(config.getParameter<double>("drVetoEM")),
+      drVetoHad_(config.getParameter<double>("drVetoHad")),
+      minEnergyEM_(config.getParameter<double>("minEnergyEM")),
+      minEnergyHad_(config.getParameter<double>("minEnergyHad")),
+      minEtEM_(config.getParameter<double>("minEtEM")),
+      minEtHad_(config.getParameter<double>("minEtHad")),
+      useEt_(config.getParameter<bool>("useEt")),
+      doRhoCorrection_(config.getParameter<bool>("doRhoCorrection")),
+      rhoMax_(config.getParameter<double>("rhoMax")),
+      rhoScale_(config.getParameter<double>("rhoScale")),
+      effectiveAreas_(config.getParameter<std::vector<double>>("effectiveAreas")) {
+  if (doRhoCorrection_) {
+    if (effectiveAreas_.size() != 2)
+      throw cms::Exception("IncompatibleVects")
+          << "effectiveAreas should have two elements for em and had components. \n";
+  }
+
+  std::string recoCandidateProducerName = "recoCandidateProducer";
+  if ((typeid(HLTHGCalLayerClusterIsolationProducer<T1>) ==
+       typeid(HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate>)))
+    recoCandidateProducerName = "recoEcalCandidateProducer";
+
+  recoCandidateProducer_ = consumes<T1Collection>(config.getParameter<edm::InputTag>(recoCandidateProducerName));
+  produces<T1IsolationMap>();
+  produces<T1IsolationMap>("em");
+  produces<T1IsolationMap>("had");
+}
+
+template <typename T1>
+void HLTHGCalLayerClusterIsolationProducer<T1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  std::string recoCandidateProducerName = "recoCandidateProducer";
+  if ((typeid(HLTHGCalLayerClusterIsolationProducer<T1>) ==
+       typeid(HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate>)))
+    recoCandidateProducerName = "recoEcalCandidateProducer";
+
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>(recoCandidateProducerName, edm::InputTag("hltL1SeededRecoEcalCandidatePF"));
+  desc.add<edm::InputTag>("layerClusterProducer", edm::InputTag("hltParticleFlowClusterECAL"));
+  desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAllCalo"));
+  desc.add<bool>("doRhoCorrection", false);
+  desc.add<bool>("useEt", false);
+  desc.add<double>("rhoMax", 9.9999999E7);
+  desc.add<double>("rhoScale", 1.0);
+  desc.add<double>("drMax", 0.3);
+  desc.add<double>("drVetoEM", 0.0);
+  desc.add<double>("drVetoHad", 0.0);
+  desc.add<double>("minEnergyEM", 0.0);
+  desc.add<double>("minEnergyHad", 0.0);
+  desc.add<double>("minEtEM", 0.0);
+  desc.add<double>("minEtHad", 0.0);
+  desc.add<std::vector<double>>("effectiveAreas", {0.0, 0.0});  // for em and had components
+  descriptions.add(defaultModuleLabel<HLTHGCalLayerClusterIsolationProducer<T1>>(), desc);
+}
+
+template <typename T1>
+void HLTHGCalLayerClusterIsolationProducer<T1>::produce(edm::Event& iEvent, const edm::EventSetup&) {
+  edm::Handle<double> rhoHandle;
+  double rho = 0.0;
+  if (doRhoCorrection_) {
+    iEvent.getByToken(rhoProducer_, rhoHandle);
+    rho = *(rhoHandle.product());
+  }
+
+  rho = std::min(rho, rhoMax_);
+  rho = rho * rhoScale_;
+
+  edm::Handle<T1Collection> recoCandHandle;
+  edm::Handle<reco::CaloClusterCollection> clusterHandle;
+
+  iEvent.getByToken(recoCandidateProducer_, recoCandHandle);
+  iEvent.getByToken(layerClusterProducer_, clusterHandle);
+
+  const std::vector<reco::CaloCluster> layerClusters = *(clusterHandle.product());
+
+  T1IsolationMap recoCandMap(recoCandHandle);
+  T1IsolationMap recoCandMapEm(recoCandHandle);
+  T1IsolationMap recoCandMapHad(recoCandHandle);
+
+  for (unsigned int iReco = 0; iReco < recoCandHandle->size(); iReco++) {
+    T1Ref candRef(recoCandHandle, iReco);
+
+    float sumEm =
+        HGCalClusterTools::emEnergyInCone(candRef->eta(),
+                                          candRef->phi(),
+                                          layerClusters,
+                                          drVetoEM_,
+                                          drMax_,
+                                          minEtEM_,
+                                          minEnergyEM_,
+                                          useEt_ ? HGCalClusterTools::EType::ET : HGCalClusterTools::EType::ENERGY);
+
+    float sumHad =
+        HGCalClusterTools::hadEnergyInCone(candRef->eta(),
+                                           candRef->phi(),
+                                           layerClusters,
+                                           drVetoHad_,
+                                           drMax_,
+                                           minEtHad_,
+                                           minEnergyHad_,
+                                           useEt_ ? HGCalClusterTools::EType::ET : HGCalClusterTools::EType::ENERGY);
+
+    if (doRhoCorrection_) {
+      sumEm = sumEm - rho * effectiveAreas_.at(0);
+      sumHad = sumHad - rho * effectiveAreas_.at(1);
+    }
+
+    float sum = sumEm + sumHad;
+
+    recoCandMap.insert(candRef, sum);
+    recoCandMapEm.insert(candRef, sumEm);
+    recoCandMapHad.insert(candRef, sumHad);
+  }
+
+  iEvent.put(std::make_unique<T1IsolationMap>(recoCandMap));
+  iEvent.put(std::make_unique<T1IsolationMap>(recoCandMapEm), "em");
+  iEvent.put(std::make_unique<T1IsolationMap>(recoCandMapHad), "had");
+}
+
+typedef HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate> EgammaHLTHGCalLayerClusterIsolationProducer;
+typedef HLTHGCalLayerClusterIsolationProducer<reco::RecoChargedCandidate> MuonHLTHGCalLayerClusterIsolationProducer;
+
+DEFINE_FWK_MODULE(EgammaHLTHGCalLayerClusterIsolationProducer);
+DEFINE_FWK_MODULE(MuonHLTHGCalLayerClusterIsolationProducer);