Merge pull request cms-sw#24407 from gpetruc/PF_102X

10_2_X backport of PF HEM mitigation

Configuration/Eras/python/Modifier_pf_badHcalMitigation_cff.py

@@ -0,0 +1,3 @@
+import FWCore.ParameterSet.Config as cms
+
+pf_badHcalMitigation =  cms.Modifier()

Configuration/StandardSequences/python/Eras.py

@@ -50,7 +50,7 @@ def __init__(self):
                            'phase2_hgcal', 'phase2_muon', 'phase2_timing', 'phase2_hgcalV9',
                            'phase2_timing_layer','phase2_timing_layer_new','phase2_hcal',
                            'trackingLowPU', 'trackingPhase1', 'ctpps_2016', 'trackingPhase2PU140','highBetaStar_2018',
-                           'tracker_apv_vfp30_2016', 'run2_miniAOD_80XLegacy','run2_miniAOD_94XFall17', 'run2_nanoAOD_92X',
+                           'tracker_apv_vfp30_2016', 'pf_badHcalMitigation', 'run2_miniAOD_80XLegacy','run2_miniAOD_94XFall17', 'run2_nanoAOD_92X',
                            'run2_nanoAOD_94XMiniAODv1', 'run2_nanoAOD_94XMiniAODv2', 'run2_nanoAOD_94X2016',
                            'hcalHardcodeConditions', 'hcalSkipPacker']
         internalUseModChains = ['run2_2017_noTrackingModifier']

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"/>

Geometry/CaloTopology/src/CaloTowerConstituentsMap.cc

@@ -69,7 +69,12 @@ void CaloTowerConstituentsMap::sort() {
 }
 
 std::vector<DetId> CaloTowerConstituentsMap::constituentsOf(const CaloTowerDetId& id) const {
+#ifdef EDM_ML_DEBUG
+  std::cout << "Get constituent of " << std::hex << id.rawId() << std::dec
+	    << " ID " << id << " ieta " << id.ieta() << std::endl;
+#endif
   std::vector<DetId> items;
+  if (id.ieta() == 0) return items;
 
   // build reverse map if needed
   if(!m_reverseItems.load(std::memory_order_acquire)) {

HLTrigger/Configuration/python/customizeHLTforCMSSW.py

@@ -186,6 +186,24 @@ def customiseForEcalTestPR22254thresholdC(process):
     return process
 
 
+def customiseFor24212(process):
+    for pfName in "hltParticleFlow", "hltParticleFlowForTaus", "hltParticleFlowReg":
+        pf = getattr(process,pfName,None)
+        if pf: # Treatment of tracks in region of bad HCal
+            pf.goodTrackDeadHcal_ptErrRel = cms.double(0.2) # trackRef->ptError()/trackRef->pt() < X
+            pf.goodTrackDeadHcal_chi2n = cms.double(5)      # trackRef->normalizedChi2() < X
+            pf.goodTrackDeadHcal_layers = cms.uint32(4)     # trackRef->hitPattern().trackerLayersWithMeasurement() >= X
+            pf.goodTrackDeadHcal_validFr = cms.double(0.5)  # trackRef->validFraction() > X
+            pf.goodTrackDeadHcal_dxy = cms.double(0.5)      # [cm] abs(trackRef->dxy(primaryVertex_.position())) < X
+            pf.goodPixelTrackDeadHcal_minEta = cms.double(2.3)   # abs(trackRef->eta()) > X
+            pf.goodPixelTrackDeadHcal_maxPt  = cms.double(50.)   # trackRef->ptError()/trackRef->pt() < X
+            pf.goodPixelTrackDeadHcal_ptErrRel = cms.double(1.0) # trackRef->ptError()/trackRef->pt() < X
+            pf.goodPixelTrackDeadHcal_chi2n = cms.double(2)      # trackRef->normalizedChi2() < X
+            pf.goodPixelTrackDeadHcal_maxLost3Hit = cms.int32(0) # max missing outer hits for a track with 3 valid pixel layers (can set to -1 to reject all these tracks)
+            pf.goodPixelTrackDeadHcal_maxLost4Hit = cms.int32(1) # max missing outer hits for a track with >= 4 valid pixel layers
+            pf.goodPixelTrackDeadHcal_dxy = cms.double(0.02)     # [cm] abs(trackRef->dxy(primaryVertex_.position())) < X
+            pf.goodPixelTrackDeadHcal_dz  = cms.double(0.05)     # [cm] abs(trackRef->dz(primaryVertex_.position())) < X
+    return process
 
 
 
@@ -194,5 +212,6 @@ def customizeHLTforCMSSW(process, menuType="GRun"):
 
     # add call to action function in proper order: newest last!
     # process = customiseFor12718(process)
+    process = customiseFor24212(process)
 
     return process

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

@@ -608,6 +608,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
    {
@@ -616,6 +619,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
@@ -683,6 +689,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
    {
@@ -691,6 +700,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

@@ -10,7 +10,9 @@
 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
 #include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
+#include "Geometry/CaloTopology/interface/HcalTopology.h"
 
+class HcalChannelQuality;
 
 class EgammaHadTower {
  public:
@@ -26,12 +28,17 @@ 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_;
+  const HcalTopology * hcalTopology_;
 };
 
 bool ClusterGreaterThan(const reco::CaloClusterPtr& c1, const reco::CaloClusterPtr& c2) ;