Dataformats and inputs only for PF HEM 15-16 mitigation

DataFormats/CaloRecHit/interface/CaloCluster.h

@@ -36,6 +36,8 @@ namespace reco {
 
     // super-cluster flags
     enum SCFlags { cleanOnly = 0, common = 100, uncleanOnly = 200 };
+    // hcal cluster flags (used for pf)
+    enum HCalFlags { badHcalMarker = 1 };
 
    //FIXME:  
    //temporary fix... to be removed before 310 final

DataFormats/EgammaCandidates/interface/GsfElectron.h

@@ -396,6 +396,7 @@ class GsfElectron : public RecoCandidate
       std::vector<CaloTowerDetId> hcalTowersBehindClusters ; //
       float hcalDepth1OverEcalBc ; // hcal over ecal seed cluster energy using 1st hcal depth (using hcal towers behind clusters)
       float hcalDepth2OverEcalBc ; // hcal over ecal seed cluster energy using 2nd hcal depth (using hcal towers behind clusters)
+      bool  invalidHcal ;          // set to true if the hcal energy estimate is not valid (e.g. the corresponding tower was off or masked)
       float sigmaIetaIphi;
       float eMax;
       float e2nd;
@@ -415,6 +416,7 @@ class GsfElectron : public RecoCandidate
 	     r9(-std::numeric_limits<float>::max()),
        hcalDepth1OverEcal(0.), hcalDepth2OverEcal(0.),
        hcalDepth1OverEcalBc(0.), hcalDepth2OverEcalBc(0.),
+       invalidHcal(false),
        sigmaIetaIphi(0.f),
        eMax(0.f),
        e2nd(0.f),
@@ -444,6 +446,7 @@ class GsfElectron : public RecoCandidate
     float hcalDepth1OverEcalBc() const { return showerShape_.hcalDepth1OverEcalBc ; }
     float hcalDepth2OverEcalBc() const { return showerShape_.hcalDepth2OverEcalBc ; }
     float hcalOverEcalBc() const { return hcalDepth1OverEcalBc() + hcalDepth2OverEcalBc() ; } 
+    float hcalOverEcalValid() const { return !showerShape_.invalidHcal; } 
     float eLeft() const { return showerShape_.eLeft; }
     float eRight() const { return showerShape_.eRight; }
     float eTop() const { return showerShape_.eTop; }
@@ -464,6 +467,7 @@ class GsfElectron : public RecoCandidate
     float full5x5_hcalDepth1OverEcalBc() const { return full5x5_showerShape_.hcalDepth1OverEcalBc ; }
     float full5x5_hcalDepth2OverEcalBc() const { return full5x5_showerShape_.hcalDepth2OverEcalBc ; }
     float full5x5_hcalOverEcalBc() const { return full5x5_hcalDepth1OverEcalBc() + full5x5_hcalDepth2OverEcalBc() ; }
+    float full5x5_hcalOverEcalValid() const { return !full5x5_showerShape_.invalidHcal; } 
     float full5x5_e2x5Left() const { return full5x5_showerShape_.e2x5Left; }
     float full5x5_e2x5Right() const { return full5x5_showerShape_.e2x5Right; }
     float full5x5_e2x5Top() const { return full5x5_showerShape_.e2x5Top; }

DataFormats/EgammaCandidates/interface/Photon.h

@@ -150,6 +150,7 @@ namespace reco {
       float hcalDepth1OverEcalBc;
       float hcalDepth2OverEcalBc;
       std::vector<CaloTowerDetId> hcalTowersBehindClusters;
+      bool  invalidHcal;
       float effSigmaRR;
       float sigmaIetaIphi;
       float sigmaIphiIphi;
@@ -177,6 +178,7 @@ namespace reco {
 	  hcalDepth2OverEcal(0),
 	  hcalDepth1OverEcalBc(0),
           hcalDepth2OverEcalBc(0),
+          invalidHcal(false),
           effSigmaRR(std::numeric_limits<float>::max()),
           sigmaIetaIphi(std::numeric_limits<float>::max()),
           sigmaIphiIphi(std::numeric_limits<float>::max()),
@@ -206,6 +208,8 @@ namespace reco {
     float hadronicDepth1OverEm() const {return  showerShapeBlock_.hcalDepth1OverEcal  ;}
     /// the  hadronic release in depth2 over electromagnetic fraction
     float hadronicDepth2OverEm() const {return  showerShapeBlock_.hcalDepth2OverEcal  ;}
+    /// returns false if hadronicOverEm is not reliably estimated (e.g. because hcal was off or masked)
+    float hadronicOverEmValid() const {return !showerShapeBlock_.invalidHcal  ;}
 
     /// the ration of hadronic energy in towers behind the BCs in the SC  and the SC energy 
     float hadTowOverEm() const {return   showerShapeBlock_.hcalDepth1OverEcalBc + showerShapeBlock_.hcalDepth2OverEcalBc  ;}
@@ -214,6 +218,8 @@ namespace reco {
     /// the ration of hadronic energy in towers depth2 behind the BCs in the SC  and the SC energy 
     float hadTowDepth2OverEm() const {return  showerShapeBlock_.hcalDepth2OverEcalBc  ;}
     const std::vector<CaloTowerDetId> & hcalTowersBehindClusters() const { return showerShapeBlock_.hcalTowersBehindClusters ; }
+    /// returns false if hadTowOverEm is not reliably estimated (e.g. because hcal was off or masked)
+    float hadTowOverEmValid() const {return !showerShapeBlock_.invalidHcal  ;}
 
     ///  Shower shape variables
     float e1x5()            const {return showerShapeBlock_.e1x5;}

DataFormats/EgammaCandidates/src/classes_def.xml

@@ -38,7 +38,8 @@
   <class name="reco::Photon::FiducialFlags" ClassVersion="10">
    <version ClassVersion="10" checksum="2148735796"/>
   </class>
-  <class name="reco::Photon::ShowerShape" ClassVersion="12">
+  <class name="reco::Photon::ShowerShape" ClassVersion="13">
+   <version ClassVersion="13" checksum="641067911"/>
    <version ClassVersion="12" checksum="1586705784"/>
    <version ClassVersion="11" checksum="2494542444"/>
   </class>
@@ -161,7 +162,8 @@
   <class name="reco::GsfElectron::FiducialFlags" ClassVersion="10">
    <version ClassVersion="10" checksum="2564471186"/>
   </class>
-  <class name="reco::GsfElectron::ShowerShape" ClassVersion="13">
+  <class name="reco::GsfElectron::ShowerShape" ClassVersion="14">
+   <version ClassVersion="14" checksum="1877043619"/>
    <version ClassVersion="13" checksum="1955217576"/>
    <version ClassVersion="12" checksum="1486422105"/>
    <version ClassVersion="11" checksum="3451347438"/>

RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h

@@ -26,7 +26,7 @@ class ElectronHcalHelper {
     double hOverEConeSize;
 
     // strategy
-    bool useTowers;
+    bool useTowers, checkHcalStatus;
 
     // specific parameters if use towers
     edm::EDGetTokenT<CaloTowerCollection> hcalTowers; 
@@ -52,7 +52,11 @@ class ElectronHcalHelper {
   std::vector<CaloTowerDetId> hcalTowersBehindClusters( const reco::SuperCluster & sc ) ;
   double hcalESumDepth1BehindClusters( const std::vector<CaloTowerDetId> & towers ) ;
   double hcalESumDepth2BehindClusters( const std::vector<CaloTowerDetId> & towers ) ;
-
+
+  // forward EgammaHadTower methods, if checkHcalStatus is enabled, using towers and H/E 
+  // otherwise, return true
+  bool hasActiveHcal( const reco::SuperCluster & sc ) ;
+
  private:
 
   const Configuration cfg_ ;

RecoEgamma/EgammaElectronAlgos/src/ElectronHcalHelper.cc

@@ -107,6 +107,16 @@ double ElectronHcalHelper::hcalESumDepth2( const SuperCluster & sc ,const std::v
    { return hcalIso_->getHcalESumDepth2(&sc) ; }
  }
 
+bool ElectronHcalHelper::hasActiveHcal( const reco::SuperCluster & sc ) 
+ {
+     if (cfg_.checkHcalStatus && cfg_.hOverEConeSize != 0 && cfg_.useTowers) {
+         return hadTower_->hasActiveHcal( sc );
+     } else {
+         return true;
+     }
+ }
+
+
 ElectronHcalHelper::~ElectronHcalHelper()
  {
   if (cfg_.hOverEConeSize==0)

RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc

@@ -609,6 +609,9 @@ void GsfElectronAlgo::calculateShowerShape( const reco::SuperClusterRef & theClu
     showerShape.hcalTowersBehindClusters = generalData_->hcalHelperPflow->hcalTowersBehindClusters(*theClus) ;
     showerShape.hcalDepth1OverEcalBc = generalData_->hcalHelperPflow->hcalESumDepth1BehindClusters(showerShape.hcalTowersBehindClusters)/theClus->energy() ;
     showerShape.hcalDepth2OverEcalBc = generalData_->hcalHelperPflow->hcalESumDepth2BehindClusters(showerShape.hcalTowersBehindClusters)/theClus->energy() ;
+    showerShape.invalidHcal = (showerShape.hcalDepth1OverEcalBc == 0 && 
+                               showerShape.hcalDepth2OverEcalBc == 0 &&
+                               !generalData_->hcalHelperPflow->hasActiveHcal(*theClus));
    }
   else
    {
@@ -617,6 +620,9 @@ void GsfElectronAlgo::calculateShowerShape( const reco::SuperClusterRef & theClu
     showerShape.hcalTowersBehindClusters = generalData_->hcalHelper->hcalTowersBehindClusters(*theClus) ;
     showerShape.hcalDepth1OverEcalBc = generalData_->hcalHelper->hcalESumDepth1BehindClusters(showerShape.hcalTowersBehindClusters)/theClus->energy() ;
     showerShape.hcalDepth2OverEcalBc = generalData_->hcalHelper->hcalESumDepth2BehindClusters(showerShape.hcalTowersBehindClusters)/theClus->energy() ;
+    showerShape.invalidHcal = (showerShape.hcalDepth1OverEcalBc == 0 && 
+                               showerShape.hcalDepth2OverEcalBc == 0 &&
+                               !generalData_->hcalHelper->hasActiveHcal(*theClus));
    }
 
   // extra shower shapes
@@ -684,6 +690,9 @@ void GsfElectronAlgo::calculateShowerShape_full5x5( const reco::SuperClusterRef
     showerShape.hcalTowersBehindClusters = generalData_->hcalHelperPflow->hcalTowersBehindClusters(*theClus) ;
     showerShape.hcalDepth1OverEcalBc = generalData_->hcalHelperPflow->hcalESumDepth1BehindClusters(showerShape.hcalTowersBehindClusters)/showerShape.e5x5 ;
     showerShape.hcalDepth2OverEcalBc = generalData_->hcalHelperPflow->hcalESumDepth2BehindClusters(showerShape.hcalTowersBehindClusters)/showerShape.e5x5 ;
+    showerShape.invalidHcal = (showerShape.hcalDepth1OverEcalBc == 0 && 
+                               showerShape.hcalDepth2OverEcalBc == 0 &&
+                               !generalData_->hcalHelperPflow->hasActiveHcal(*theClus));
    }
   else
    {
@@ -692,6 +701,9 @@ void GsfElectronAlgo::calculateShowerShape_full5x5( const reco::SuperClusterRef
     showerShape.hcalTowersBehindClusters = generalData_->hcalHelper->hcalTowersBehindClusters(*theClus) ;
     showerShape.hcalDepth1OverEcalBc = generalData_->hcalHelper->hcalESumDepth1BehindClusters(showerShape.hcalTowersBehindClusters)/showerShape.e5x5 ;
     showerShape.hcalDepth2OverEcalBc = generalData_->hcalHelper->hcalESumDepth2BehindClusters(showerShape.hcalTowersBehindClusters)/showerShape.e5x5 ;
+    showerShape.invalidHcal = (showerShape.hcalDepth1OverEcalBc == 0 && 
+                               showerShape.hcalDepth2OverEcalBc == 0 &&
+                               !generalData_->hcalHelper->hasActiveHcal(*theClus));
    }
 
   // extra shower shapes

RecoEgamma/EgammaElectronProducers/plugins/GsfElectronBaseProducer.cc

@@ -283,6 +283,7 @@ GsfElectronBaseProducer::GsfElectronBaseProducer( const edm::ParameterSet& cfg,
   if (hcalCfg_.hOverEConeSize>0)
    {
     hcalCfg_.useTowers = true ;
+    hcalCfg_.checkHcalStatus = cfg.getParameter<bool>("checkHcalStatus") ;
     hcalCfg_.hcalTowers = 
       consumes<CaloTowerCollection>(cfg.getParameter<edm::InputTag>("hcalTowers")) ;
     hcalCfg_.hOverEPtMin = cfg.getParameter<double>("hOverEPtMin") ;
@@ -291,6 +292,7 @@ GsfElectronBaseProducer::GsfElectronBaseProducer( const edm::ParameterSet& cfg,
   if (hcalCfgPflow_.hOverEConeSize>0)
    {
     hcalCfgPflow_.useTowers = true ;
+    hcalCfgPflow_.checkHcalStatus = cfg.getParameter<bool>("checkHcalStatus") ;
     hcalCfgPflow_.hcalTowers = 
       consumes<CaloTowerCollection>(cfg.getParameter<edm::InputTag>("hcalTowers")) ;
     hcalCfgPflow_.hOverEPtMin = cfg.getParameter<double>("hOverEPtMinPflow") ;
@@ -334,7 +336,7 @@ GsfElectronBaseProducer::GsfElectronBaseProducer( const edm::ParameterSet& cfg,
   EcalClusterFunctionBaseClass * superClusterErrorFunction = nullptr ;
   std::string superClusterErrorFunctionName
    = cfg.getParameter<std::string>("superClusterErrorFunction") ;
-  if (superClusterErrorFunctionName!="")
+  if (!superClusterErrorFunctionName.empty())
    {
     superClusterErrorFunction
      = EcalClusterFunctionFactory::get()->create(superClusterErrorFunctionName,cfg) ;
@@ -347,7 +349,7 @@ GsfElectronBaseProducer::GsfElectronBaseProducer( const edm::ParameterSet& cfg,
   EcalClusterFunctionBaseClass * crackCorrectionFunction = nullptr ;
   std::string crackCorrectionFunctionName
    = cfg.getParameter<std::string>("crackCorrectionFunction") ;
-  if (crackCorrectionFunctionName!="")
+  if (!crackCorrectionFunctionName.empty())
    {
     crackCorrectionFunction
      = EcalClusterFunctionFactory::get()->create(crackCorrectionFunctionName,cfg) ;

RecoEgamma/EgammaElectronProducers/python/gedGsfElectrons_cfi.py

@@ -15,6 +15,7 @@
     barrelRecHitCollectionTag = cms.InputTag("ecalRecHit","EcalRecHitsEB"),
     endcapRecHitCollectionTag = cms.InputTag("ecalRecHit","EcalRecHitsEE"),
     hcalTowers = cms.InputTag("towerMaker"),
+    checkHcalStatus = cms.bool(True),
     pfMvaTag =  cms.InputTag(""),
     seedsTag = cms.InputTag("ecalDrivenElectronSeeds"),
     beamSpotTag = cms.InputTag("offlineBeamSpot"),

RecoEgamma/EgammaElectronProducers/python/gsfElectrons_cfi.py

@@ -18,6 +18,7 @@
     barrelRecHitCollectionTag = cms.InputTag("ecalRecHit","EcalRecHitsEB"),
     endcapRecHitCollectionTag = cms.InputTag("ecalRecHit","EcalRecHitsEE"),
     hcalTowers = cms.InputTag("towerMaker"),
+    checkHcalStatus = cms.bool(True),
     pfMvaTag =  cms.InputTag(""),
     seedsTag = cms.InputTag("ecalDrivenElectronSeeds"),
     beamSpotTag = cms.InputTag("offlineBeamSpot"),
@@ -181,6 +182,7 @@
     pflowGsfElectronsTag = cms.InputTag("pfElectronTranslator:pf"),
     gsfElectronCoresTag = cms.InputTag("gsfElectronCores"),
     hcalTowers = cms.InputTag("towerMaker"),
+    checkHcalStatus = cms.bool(True),
     barrelRecHitCollectionTag = cms.InputTag("ecalRecHit","EcalRecHitsEB"),
     endcapRecHitCollectionTag = cms.InputTag("ecalRecHit","EcalRecHitsEE"),
     pfMvaTag =  cms.InputTag("pfElectronTranslator:pf"),

RecoEgamma/EgammaIsolationAlgos/BuildFile.xml

@@ -8,8 +8,10 @@
 <use   name="DataFormats/RecoCandidate"/>
 <use   name="DataFormats/ParticleFlowReco"/>
 <use   name="DataFormats/ParticleFlowCandidate"/>
+<use   name="DataFormats/HcalDetId"/>
 <use   name="DataFormats/PatCandidates"/>
 <use   name="CondFormats/EcalObjects"/>
+<use   name="CondFormats/HcalObjects"/>
 <use   name="CondFormats/DataRecord"/>
 <use   name="RecoLocalCalo/EcalRecAlgos"/>
 <use   name="CommonTools/Statistics"/>

RecoEgamma/EgammaIsolationAlgos/interface/EgammaHadTower.h

@@ -11,6 +11,7 @@
 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
 #include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
 
+class HcalChannelQuality;
 
 class EgammaHadTower {
  public:
@@ -26,12 +27,16 @@ class EgammaHadTower {
   std::vector<CaloTowerDetId> towersOf(const reco::SuperCluster& sc) const ;
   CaloTowerDetId  towerOf(const reco::CaloCluster& cluster) const ;
   void setTowerCollection(const CaloTowerCollection* towercollection);
+  bool hasActiveHcal( const reco::SuperCluster & sc ) const ;
+  bool hasActiveHcal( const std::vector<CaloTowerDetId> & towers ) const ;
+
 
  private:
   const CaloTowerConstituentsMap * towerMap_;
   HoeMode mode_;
   const CaloTowerCollection * towerCollection_;
   unsigned int NMaxClusters_;
+  const HcalChannelQuality * hcalQuality_;
 };
 
 bool ClusterGreaterThan(const reco::CaloClusterPtr& c1, const reco::CaloClusterPtr& c2) ;

RecoEgamma/EgammaIsolationAlgos/src/EgammaHadTower.cc

@@ -3,6 +3,12 @@
 #include "FWCore/Framework/interface/ESHandle.h"
 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
 #include "Geometry/Records/interface/CaloGeometryRecord.h"
+#include "CondFormats/HcalObjects/interface/HcalChannelStatus.h"
+#include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
+#include "CondFormats/HcalObjects/interface/HcalCondObjectContainer.h"
+#include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
+#include "DataFormats/HcalDetId/interface/HcalDetId.h"
+
 
 #include <algorithm>
 #include <iostream>
@@ -12,6 +18,10 @@ EgammaHadTower::EgammaHadTower(const edm::EventSetup &es,HoeMode mode):mode_(mod
   es.get<CaloGeometryRecord>().get(ctmaph);
   towerMap_ = &(*ctmaph);
   NMaxClusters_ = 4;
+
+  edm::ESHandle<HcalChannelQuality> hQuality;
+  es.get<HcalChannelQualityRcd>().get("withTopo",hQuality);
+  hcalQuality_ = hQuality.product();
 }
 
 CaloTowerDetId  EgammaHadTower::towerOf(const reco::CaloCluster& cluster) const {
@@ -100,6 +110,41 @@ double EgammaHadTower::getDepth2HcalESum(const std::vector<CaloTowerDetId> & tow
   return esum;
 }
 
+bool EgammaHadTower::hasActiveHcal( const std::vector<CaloTowerDetId> & towers ) const {
+  bool debug = false; // change this to true to get debug output
+  bool active = false;
+  int statusMask = ((1<<HcalChannelStatus::HcalCellOff) | (1<<HcalChannelStatus::HcalCellMask) | (1<<HcalChannelStatus::HcalCellDead));
+  if (debug) std::cout << "DEBUG: hasActiveHcal called with " << towers.size() << " detids. First tower detid ieta " << towers.front().ieta() << " iphi " << towers.front().iphi() << std::endl;
+  for (auto towerid : towers) {
+      unsigned int ngood = 0, nbad = 0;
+      for (DetId id : towerMap_->constituentsOf(towerid)) {
+          if (id.det() != DetId::Hcal) {
+              continue;
+          }
+          HcalDetId hid(id);
+          if (debug) std::cout << "      hcal constituent on subdet " << hid.subdet() << ", ieta " << hid.ieta() << " iphi " << hid.iphi() << ", depth " << hid.depth() << std::endl;
+          if (hid.subdet() != HcalBarrel && hid.subdet() != HcalEndcap) continue;
+          const auto * values = hcalQuality_->getValues(id,/*throwOnFail=*/false);
+          if (!values) {
+              if (debug) std::cout << "          MISSING IN CONDITIONS" << std::endl;
+              continue;
+          }
+          int status = values->getValue();
+          if (status & statusMask) {
+              if (debug) std::cout << "          BAD!" << std::endl;
+              nbad++;
+          } else {
+              ngood++;
+          }
+      }
+      if (debug) std::cout << "    overall ngood " << ngood << " nbad " << nbad << std::endl;
+      if (nbad == 0 || (ngood > 0 && nbad < ngood)) {
+          active = true;
+      }
+  }
+  return active;
+}
+
 double EgammaHadTower::getDepth1HcalESum( const reco::SuperCluster& sc ) const {
   return getDepth1HcalESum(towersOf(sc)) ;
 }
@@ -112,6 +157,10 @@ void EgammaHadTower::setTowerCollection(const CaloTowerCollection* towerCollecti
   towerCollection_ = towerCollection;
 }
 
+bool EgammaHadTower::hasActiveHcal( const reco::SuperCluster & sc ) const {
+    return hasActiveHcal(towersOf(sc)) ;
+}
+
 bool ClusterGreaterThan(const reco::CaloClusterPtr& c1, const reco::CaloClusterPtr& c2)  {
   return (*c1 > *c2);
 }

RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonProducer.h

@@ -160,5 +160,7 @@ class GEDPhotonProducer : public edm::stream::EDProducer<> {
   PhotonEnergyCorrector* thePhotonEnergyCorrector_;
   std::string  candidateP4type_;
 
+  bool checkHcalStatus_;
+
 };
 #endif

RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py

@@ -87,6 +87,7 @@
     RecHitSeverityToBeExcludedEB = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded,
     RecHitFlagToBeExcludedEE = multi5x5BasicClustersCleaned.RecHitFlagToBeExcluded,
     RecHitSeverityToBeExcludedEE = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded,
+    checkHcalStatus = cms.bool(True),
 )
 
 

RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py

@@ -78,6 +78,7 @@
     RecHitSeverityToBeExcludedEB = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded,
     RecHitFlagToBeExcludedEE = multi5x5BasicClustersCleaned.RecHitFlagToBeExcluded,
     RecHitSeverityToBeExcludedEE = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded,
+    checkHcalStatus = cms.bool(True),
 )
 
 photonsFromMultiCl = photons.clone(