diff --git a/Configuration/Eras/python/Modifier_pf_badHcalMitigation_cff.py b/Configuration/Eras/python/Modifier_pf_badHcalMitigation_cff.py new file mode 100644 index 0000000000000..3a29156e0ba85 --- /dev/null +++ b/Configuration/Eras/python/Modifier_pf_badHcalMitigation_cff.py @@ -0,0 +1,3 @@ +import FWCore.ParameterSet.Config as cms + +pf_badHcalMitigation = cms.Modifier() diff --git a/Configuration/StandardSequences/python/Eras.py b/Configuration/StandardSequences/python/Eras.py index 4a70e75741dff..617e3743f4bda 100644 --- a/Configuration/StandardSequences/python/Eras.py +++ b/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'] diff --git a/DataFormats/CaloRecHit/interface/CaloCluster.h b/DataFormats/CaloRecHit/interface/CaloCluster.h index 6ec3df3845110..002d3b4b76c6e 100644 --- a/DataFormats/CaloRecHit/interface/CaloCluster.h +++ b/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 diff --git a/DataFormats/EgammaCandidates/interface/GsfElectron.h b/DataFormats/EgammaCandidates/interface/GsfElectron.h index 9cd1140d5cee9..0071ddafe308d 100644 --- a/DataFormats/EgammaCandidates/interface/GsfElectron.h +++ b/DataFormats/EgammaCandidates/interface/GsfElectron.h @@ -396,6 +396,7 @@ class GsfElectron : public RecoCandidate std::vector 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::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; } diff --git a/DataFormats/EgammaCandidates/interface/Photon.h b/DataFormats/EgammaCandidates/interface/Photon.h index 5e915eee8c10c..932ca599398be 100644 --- a/DataFormats/EgammaCandidates/interface/Photon.h +++ b/DataFormats/EgammaCandidates/interface/Photon.h @@ -150,6 +150,7 @@ namespace reco { float hcalDepth1OverEcalBc; float hcalDepth2OverEcalBc; std::vector 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::max()), sigmaIetaIphi(std::numeric_limits::max()), sigmaIphiIphi(std::numeric_limits::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 & 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;} diff --git a/DataFormats/EgammaCandidates/src/classes_def.xml b/DataFormats/EgammaCandidates/src/classes_def.xml index a1f16258f969a..b6ccb80436975 100644 --- a/DataFormats/EgammaCandidates/src/classes_def.xml +++ b/DataFormats/EgammaCandidates/src/classes_def.xml @@ -38,7 +38,8 @@ - + + @@ -161,7 +162,8 @@ - + + diff --git a/Geometry/CaloTopology/src/CaloTowerConstituentsMap.cc b/Geometry/CaloTopology/src/CaloTowerConstituentsMap.cc index 6f943500948a0..2e32c4bb7558e 100644 --- a/Geometry/CaloTopology/src/CaloTowerConstituentsMap.cc +++ b/Geometry/CaloTopology/src/CaloTowerConstituentsMap.cc @@ -69,7 +69,12 @@ void CaloTowerConstituentsMap::sort() { } std::vector 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 items; + if (id.ieta() == 0) return items; // build reverse map if needed if(!m_reverseItems.load(std::memory_order_acquire)) { diff --git a/HLTrigger/Configuration/python/customizeHLTforCMSSW.py b/HLTrigger/Configuration/python/customizeHLTforCMSSW.py index 7a71d4c06a3d0..f70aac4f6c0b9 100644 --- a/HLTrigger/Configuration/python/customizeHLTforCMSSW.py +++ b/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 diff --git a/RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h b/RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h index cea421af1c456..1868b848b327f 100644 --- a/RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h +++ b/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 hcalTowers; @@ -52,7 +52,11 @@ class ElectronHcalHelper { std::vector hcalTowersBehindClusters( const reco::SuperCluster & sc ) ; double hcalESumDepth1BehindClusters( const std::vector & towers ) ; double hcalESumDepth2BehindClusters( const std::vector & 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_ ; diff --git a/RecoEgamma/EgammaElectronAlgos/src/ElectronHcalHelper.cc b/RecoEgamma/EgammaElectronAlgos/src/ElectronHcalHelper.cc index 984d3b61ff23d..f67535f8452e5 100644 --- a/RecoEgamma/EgammaElectronAlgos/src/ElectronHcalHelper.cc +++ b/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) diff --git a/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc b/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc index 0159246b94345..8ac4a316778a7 100644 --- a/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc +++ b/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 diff --git a/RecoEgamma/EgammaElectronProducers/plugins/GsfElectronBaseProducer.cc b/RecoEgamma/EgammaElectronProducers/plugins/GsfElectronBaseProducer.cc index a7470c1da1ce3..29afeaa25e07c 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/GsfElectronBaseProducer.cc +++ b/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("checkHcalStatus") ; hcalCfg_.hcalTowers = consumes(cfg.getParameter("hcalTowers")) ; hcalCfg_.hOverEPtMin = cfg.getParameter("hOverEPtMin") ; @@ -291,6 +292,7 @@ GsfElectronBaseProducer::GsfElectronBaseProducer( const edm::ParameterSet& cfg, if (hcalCfgPflow_.hOverEConeSize>0) { hcalCfgPflow_.useTowers = true ; + hcalCfgPflow_.checkHcalStatus = cfg.getParameter("checkHcalStatus") ; hcalCfgPflow_.hcalTowers = consumes(cfg.getParameter("hcalTowers")) ; hcalCfgPflow_.hOverEPtMin = cfg.getParameter("hOverEPtMinPflow") ; @@ -334,7 +336,7 @@ GsfElectronBaseProducer::GsfElectronBaseProducer( const edm::ParameterSet& cfg, EcalClusterFunctionBaseClass * superClusterErrorFunction = nullptr ; std::string superClusterErrorFunctionName = cfg.getParameter("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("crackCorrectionFunction") ; - if (crackCorrectionFunctionName!="") + if (!crackCorrectionFunctionName.empty()) { crackCorrectionFunction = EcalClusterFunctionFactory::get()->create(crackCorrectionFunctionName,cfg) ; diff --git a/RecoEgamma/EgammaElectronProducers/python/gedGsfElectrons_cfi.py b/RecoEgamma/EgammaElectronProducers/python/gedGsfElectrons_cfi.py index 46c6868e7cd47..887f31e9a5b18 100644 --- a/RecoEgamma/EgammaElectronProducers/python/gedGsfElectrons_cfi.py +++ b/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"), diff --git a/RecoEgamma/EgammaElectronProducers/python/gsfElectrons_cfi.py b/RecoEgamma/EgammaElectronProducers/python/gsfElectrons_cfi.py index f6b80f0ce9b05..00050cc0963f7 100644 --- a/RecoEgamma/EgammaElectronProducers/python/gsfElectrons_cfi.py +++ b/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"), diff --git a/RecoEgamma/EgammaIsolationAlgos/BuildFile.xml b/RecoEgamma/EgammaIsolationAlgos/BuildFile.xml index 254ee40b15074..7f8da9c416aad 100644 --- a/RecoEgamma/EgammaIsolationAlgos/BuildFile.xml +++ b/RecoEgamma/EgammaIsolationAlgos/BuildFile.xml @@ -8,8 +8,10 @@ + + diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/EgammaHadTower.h b/RecoEgamma/EgammaIsolationAlgos/interface/EgammaHadTower.h index 4157952a9483a..d6e181fe9b135 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/EgammaHadTower.h +++ b/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 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 & 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) ; diff --git a/RecoEgamma/EgammaIsolationAlgos/src/EgammaHadTower.cc b/RecoEgamma/EgammaIsolationAlgos/src/EgammaHadTower.cc index 60eb9d92eb01b..d58e702e16110 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/EgammaHadTower.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/EgammaHadTower.cc @@ -3,15 +3,34 @@ #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 "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/Records/interface/HcalRecNumberingRecord.h" #include #include +//#define EDM_ML_DEBUG + EgammaHadTower::EgammaHadTower(const edm::EventSetup &es,HoeMode mode):mode_(mode) { edm::ESHandle ctmaph; es.get().get(ctmaph); towerMap_ = &(*ctmaph); NMaxClusters_ = 4; + + edm::ESHandle hQuality; + es.get().get("withTopo",hQuality); + hcalQuality_ = hQuality.product(); + + edm::ESHandle hcalTopology; + es.get().get( hcalTopology ); + hcalTopology_ = hcalTopology.product(); + } CaloTowerDetId EgammaHadTower::towerOf(const reco::CaloCluster& cluster) const { @@ -52,6 +71,9 @@ std::vector EgammaHadTower::towersOf(const reco::SuperCluster& for ( unsigned iclus =0 ; iclus ::const_iterator itcheck=find(towers.begin(),towers.end(),id); if( itcheck == towers.end() ) { towers.push_back(id); @@ -100,6 +122,58 @@ double EgammaHadTower::getDepth2HcalESum(const std::vector & tow return esum; } +bool EgammaHadTower::hasActiveHcal( const std::vector & towers ) const { + + bool active = false; + int statusMask = ((1<constituentsOf(towerid)) { + if (id.det() != DetId::Hcal) { + continue; + } + HcalDetId hid(id); + if (hid.subdet() != HcalBarrel && hid.subdet() != HcalEndcap) continue; +#ifdef EDM_ML_DEBUG + std::cout << "EgammaHadTower DetId " << std::hex << id.rawId() + << " hid.rawId " << hid.rawId() << std::dec + << " sub " << hid.subdet() << " ieta " << hid.ieta() + << " iphi " << hid.iphi() + << " depth " << hid.depth() << std::endl; +#endif + // Sunanda's fix for 2017 Plan1 + // and removed protection + int status = hcalQuality_->getValues((DetId)(hcalTopology_->idFront(HcalDetId(id))),/*throwOnFail=*/true)->getValue(); + +#ifdef EDM_ML_DEBUG + std::cout << "channels status = " << std::hex << status << std::dec + << " int value = " << status << std::endl; +#endif + + if (status & statusMask) { +#ifdef EDM_ML_DEBUG + std::cout << " BAD!" << std::endl; +#endif + nbad++; + } else { + ngood++; + } + } +#ifdef EDM_ML_DEBUG + std::cout << " overall ngood " << ngood << " nbad " << nbad << "\n"; +#endif + 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 +186,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); } diff --git a/RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonProducer.h b/RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonProducer.h index 7bd9ddc1c7ae6..69c5a5fb19c08 100644 --- a/RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonProducer.h +++ b/RecoEgamma/EgammaPhotonProducers/interface/GEDPhotonProducer.h @@ -160,5 +160,7 @@ class GEDPhotonProducer : public edm::stream::EDProducer<> { PhotonEnergyCorrector* thePhotonEnergyCorrector_; std::string candidateP4type_; + bool checkHcalStatus_; + }; #endif diff --git a/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py index c984055a69a6d..42a637378eaac 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/gedPhotons_cfi.py @@ -87,6 +87,7 @@ RecHitSeverityToBeExcludedEB = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded, RecHitFlagToBeExcludedEE = multi5x5BasicClustersCleaned.RecHitFlagToBeExcluded, RecHitSeverityToBeExcludedEE = cleanedHybridSuperClusters.RecHitSeverityToBeExcluded, + checkHcalStatus = cms.bool(True), ) diff --git a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py index 75b237bad14cc..7912ab28a8487 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/photons_cfi.py +++ b/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( diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index dc535aedabc88..50ee3b8bc9fba 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -227,6 +227,9 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config) : thePhotonIsolationCalculator_=nullptr; thePhotonMIPHaloTagger_=nullptr; } + + checkHcalStatus_ = conf_.getParameter("checkHcalStatus"); + // Register the product produces< reco::PhotonCollection >(photonCollection_); if (not pfEgammaCandidates_.isUninitialized()) @@ -549,6 +552,7 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, std::vector TowersBehindClus; float hcalDepth1OverEcalBc,hcalDepth2OverEcalBc; hcalDepth1OverEcalBc=hcalDepth2OverEcalBc=0.f; + bool invalidHcal = false; if (not hcalTowers_.isUninitialized()) { const CaloTowerCollection* hcalTowersColl = hcalTowersHandle.product(); @@ -562,6 +566,10 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, TowersBehindClus = towerIsoBehindClus.towersOf(*scRef); hcalDepth1OverEcalBc = towerIsoBehindClus.getDepth1HcalESum(TowersBehindClus)/scRef->energy(); hcalDepth2OverEcalBc = towerIsoBehindClus.getDepth2HcalESum(TowersBehindClus)/scRef->energy(); + + if (checkHcalStatus_ && hcalDepth1OverEcalBc == 0 && hcalDepth2OverEcalBc == 0) { + invalidHcal = !towerIsoBehindClus.hasActiveHcal(TowersBehindClus); + } } // std::cout << " GEDPhotonProducer calculation of HoE with towers in a cone " << HoE1 << " " << HoE2 << std::endl; @@ -637,6 +645,7 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, showerShape.hcalDepth1OverEcalBc = hcalDepth1OverEcalBc; showerShape.hcalDepth2OverEcalBc = hcalDepth2OverEcalBc; showerShape.hcalTowersBehindClusters = TowersBehindClus; + showerShape.invalidHcal = invalidHcal; /// fill extra shower shapes const float spp = (!edm::isFinite(locCov[2]) ? 0. : sqrt(locCov[2])); const float sep = locCov[1]; diff --git a/RecoJets/JetProducers/python/caloJetsForTrk_cff.py b/RecoJets/JetProducers/python/caloJetsForTrk_cff.py index 96a649c847b62..fad402d974e77 100644 --- a/RecoJets/JetProducers/python/caloJetsForTrk_cff.py +++ b/RecoJets/JetProducers/python/caloJetsForTrk_cff.py @@ -15,6 +15,9 @@ caloJetsForTrkTask = cms.Task(caloTowerForTrk,ak4CaloJetsForTrk) caloJetsForTrk = cms.Sequence(caloJetsForTrkTask) +from Configuration.Eras.Modifier_pf_badHcalMitigation_cff import pf_badHcalMitigation +pf_badHcalMitigation.toModify( caloTowerForTrk, missingHcalRescaleFactorForEcal = 1.0 ) + from Configuration.Eras.Modifier_phase2_hcal_cff import phase2_hcal phase2_hcal.toModify( caloTowerForTrk, hbheInput = cms.InputTag("hbhereco") ) diff --git a/RecoLocalCalo/CaloTowersCreator/interface/CaloTowersCreationAlgo.h b/RecoLocalCalo/CaloTowersCreator/interface/CaloTowersCreationAlgo.h index a7b2047cd206f..2534ea23a0b2b 100644 --- a/RecoLocalCalo/CaloTowersCreator/interface/CaloTowersCreationAlgo.h +++ b/RecoLocalCalo/CaloTowersCreator/interface/CaloTowersCreationAlgo.h @@ -186,6 +186,7 @@ class CaloTowersCreationAlgo { void setUseRejectedRecoveredHcalHits(bool flag) {useRejectedRecoveredHcalHits = flag; }; void setUseRejectedRecoveredEcalHits(bool flag) {useRejectedRecoveredEcalHits = flag; }; + void setMissingHcalRescaleFactorForEcal(float factor) {missingHcalRescaleFactorForEcal = factor; }; //------------------------------------------------------------------------------------------------------------------- @@ -315,6 +316,8 @@ class CaloTowersCreationAlgo { unsigned int useRejectedRecoveredHcalHits; unsigned int useRejectedRecoveredEcalHits; + // if Hcal is missing, fudge it scaling Ecal by this factor (if > 0) + float missingHcalRescaleFactorForEcal; /// only affects energy and ET calculation. HO is still recorded in the tower bool theHOIsUsed; @@ -340,7 +343,7 @@ class CaloTowersCreationAlgo { // Number of channels in the tower that were not used in RecHit production (dead/off,...). // These channels are added to the other "bad" channels found in the recHit collection. - typedef std::map HcalDropChMap; + typedef std::map> HcalDropChMap; HcalDropChMap hcalDropChMap; // Number of bad Ecal channel in each tower diff --git a/RecoLocalCalo/CaloTowersCreator/python/calotowermaker_cfi.py b/RecoLocalCalo/CaloTowersCreator/python/calotowermaker_cfi.py index 154c70c23cc06..06e3db9e680eb 100644 --- a/RecoLocalCalo/CaloTowersCreator/python/calotowermaker_cfi.py +++ b/RecoLocalCalo/CaloTowersCreator/python/calotowermaker_cfi.py @@ -138,6 +138,9 @@ UseRejectedRecoveredHcalHits = cms.bool(True), UseRejectedRecoveredEcalHits = cms.bool(False), +# If Hcal is masked, and Ecal is present, pretend Hcal = (this factor) * Ecal + missingHcalRescaleFactorForEcal = cms.double(0), + # flag to allow/disallow missing inputs AllowMissingInputs = cms.bool(False), diff --git a/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreationAlgo.cc b/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreationAlgo.cc index d027cc52ea354..192e93c8e0e8f 100644 --- a/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreationAlgo.cc +++ b/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreationAlgo.cc @@ -81,6 +81,7 @@ CaloTowersCreationAlgo::CaloTowersCreationAlgo() theHcalAcceptSeverityLevelForRejectedHit(0), useRejectedRecoveredHcalHits(0), useRejectedRecoveredEcalHits(0), + missingHcalRescaleFactorForEcal(0.), theHOIsUsed(true), // (for momentum reconstruction algorithm) theMomConstrMethod(0), @@ -187,6 +188,7 @@ CaloTowersCreationAlgo::CaloTowersCreationAlgo(double EBthreshold, double EEthre theHcalAcceptSeverityLevelForRejectedHit(0), useRejectedRecoveredHcalHits(0), useRejectedRecoveredEcalHits(0), + missingHcalRescaleFactorForEcal(0.), theHOIsUsed(useHO), // (momentum reconstruction algorithm) theMomConstrMethod(momConstrMethod), @@ -300,6 +302,7 @@ CaloTowersCreationAlgo::CaloTowersCreationAlgo(double EBthreshold, double EEthre theHcalAcceptSeverityLevelForRejectedHit(0), useRejectedRecoveredHcalHits(0), useRejectedRecoveredEcalHits(0), + missingHcalRescaleFactorForEcal(0.), theHOIsUsed(useHO), // (momentum reconstruction algorithm) theMomConstrMethod(momConstrMethod), @@ -939,6 +942,14 @@ void CaloTowersCreationAlgo::convert(const CaloTowerDetId& id, const MetaTower& if(metaContains.empty()) return; + if (missingHcalRescaleFactorForEcal > 0 && E_had == 0 && E_em > 0) { + auto match = hcalDropChMap.find(id); + if (match != hcalDropChMap.end() && match->second.second) { + E_had = missingHcalRescaleFactorForEcal * E_em; + E += E_had; + } + } + double E_had_tot = (theHOIsUsed && id.ietaAbs()<=theTowerTopology->lastHORing())? E_had+E_outer : E_had; @@ -1107,7 +1118,7 @@ void CaloTowersCreationAlgo::convert(const CaloTowerDetId& id, const MetaTower& // now add dead/off/... channels not used in RecHit reconstruction for HCAL HcalDropChMap::iterator dropChItr = hcalDropChMap.find(id); - if (dropChItr != hcalDropChMap.end()) numBadHcalChan += dropChItr->second; + if (dropChItr != hcalDropChMap.end()) numBadHcalChan += dropChItr->second.first; // for ECAL the number of all bad channels is obtained here ----------------------- @@ -1618,7 +1629,7 @@ void CaloTowersCreationAlgo::makeHcalDropChMap() { CaloTowerDetId twrId = theTowerConstituentsMap->towerOf(id); - hcalDropChMap[twrId] +=1; + hcalDropChMap[twrId].first +=1; HcalDetId hid(*it); @@ -1629,11 +1640,34 @@ void CaloTowersCreationAlgo::makeHcalDropChMap() { bool merge = theHcalTopology->mergedDepth29(hid); if (merge) { CaloTowerDetId twrId29(twrId.ieta()+twrId.zside(), twrId.iphi()); - hcalDropChMap[twrId29] +=1; + hcalDropChMap[twrId29].first +=1; } } } - + } + // now I know how many bad channels, but I also need to know if there's any good ones + if (missingHcalRescaleFactorForEcal > 0) { + for (auto & pair : hcalDropChMap) { + if (pair.second.first == 0) continue; // unexpected, but just in case + int ngood = 0, nbad = 0; + for (DetId id : theTowerConstituentsMap->constituentsOf(pair.first)) { + if (id.det() != DetId::Hcal) continue; + HcalDetId hid(id); + if (hid.subdet() != HcalBarrel && hid.subdet() != HcalEndcap) continue; + const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue(); + if (dbStatusFlag == 0 || ! theHcalSevLvlComputer->dropChannel(dbStatusFlag)) { + ngood += 1; + } else { + nbad += 1; // recount, since pair.second.first may include HO + } + } + if (nbad > 0 && nbad >= ngood) { + //uncomment for debug (may be useful to tune the criteria above) + //CaloTowerDetId id(pair.first); + //std::cout << "CaloTower at ieta = " << id.ieta() << ", iphi " << id.iphi() << ": set Hcal as not efficient (ngood =" << ngood << ", nbad = " << nbad << ")" << std::endl; + pair.second.second = true; + } + } } } diff --git a/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreator.cc b/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreator.cc index 960e195a9f136..7ec0a9e9e81f7 100644 --- a/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreator.cc +++ b/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreator.cc @@ -8,6 +8,7 @@ #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h" #include "CommonTools/Utils/interface/StringToEnumValue.h" +//#define EDM_ML_DEBUG CaloTowersCreator::CaloTowersCreator(const edm::ParameterSet& conf) : algo_(conf.getParameter("EBThreshold"), @@ -92,6 +93,7 @@ CaloTowersCreator::CaloTowersCreator(const edm::ParameterSet& conf) : { + algo_.setMissingHcalRescaleFactorForEcal(conf.getParameter("missingHcalRescaleFactorForEcal")); // register for data access tok_hbhe_ = consumes(conf.getParameter("hbheInput")); @@ -125,13 +127,13 @@ CaloTowersCreator::CaloTowersCreator(const edm::ParameterSet& conf) : if (eScales_.instanceLabel.empty()) produces(); else produces(eScales_.instanceLabel); - /* +#ifdef EDM_ML_DEBUG std::cout << "VI Producer " << (useRejectedHitsOnly_ ? "use rejectOnly " : " ") << (allowMissingInputs_ ? "allowMissing " : " " ) << nLabels << ' ' << severitynames.size() << std::endl; - */ +#endif } void CaloTowersCreator::produce(edm::Event& e, const edm::EventSetup& c) { @@ -199,7 +201,7 @@ void CaloTowersCreator::produce(edm::Event& e, const edm::EventSetup& c) { algo_.setUseRejectedRecoveredHcalHits(useRejectedRecoveredHcalHits_); algo_.setUseRejectedRecoveredEcalHits(useRejectedRecoveredEcalHits_); - /* +#ifdef EDM_ML_DEBUG std::cout << "VI Produce: " << (useRejectedHitsOnly_ ? "use rejectOnly " : " ") << (allowMissingInputs_ ? "allowMissing " : " " ) @@ -208,7 +210,7 @@ void CaloTowersCreator::produce(edm::Event& e, const edm::EventSetup& c) { << ' ' << theEcalSeveritiesToBeExcluded_.size() << ' ' << theEcalSeveritiesToBeUsedInBadTowers_.size() << std::endl; - */ +#endif algo_.begin(); // clear the internal buffer @@ -283,12 +285,18 @@ void CaloTowersCreator::produce(edm::Event& e, const edm::EventSetup& c) { // Step C: Process algo_.finish(*prod); - /* +#ifdef EDM_ML_DEBUG int totc=0; float totE=0; reco::LeafCandidate::LorentzVector totP4; - for (auto const & tw : (*prod) ) { totc += tw.constituents().size(); totE+=tw.energy(); totP4+=tw.p4();} + for (auto const & tw : (*prod) ) { + totc += tw.constituents().size(); + totE+=tw.energy(); + totP4+=tw.p4(); + std::cout << "CaloTowerCreator: " << tw.id() << " with E " << tw.energy() + << " and " << tw.constituents().size() << " constituents\n"; + } std::cout << "VI " << (*prod).size() << " " << totc << " " << totE << " " << totP4 << std::endl; - */ +#endif // Step D: Put into the event if (eScales_.instanceLabel.empty()) e.put(std::move(prod)); @@ -341,6 +349,7 @@ void CaloTowersCreator::fillDescriptions(edm::ConfigurationDescriptions& descrip desc.add("UseRejectedHitsOnly", false); desc.add("UseRejectedRecoveredHcalHits", true); desc.add("UseRejectedRecoveredEcalHits", false); + desc.add("missingHcalRescaleFactorForEcal", 0.0); desc.add("AllowMissingInputs", false); desc.add >("HBGrid", {-1.0, 1.0, 10.0, 100.0, 1000.0}); desc.add >("EEWeights", {1.0, 1.0, 1.0, 1.0, 1.0}); diff --git a/RecoParticleFlow/PFClusterProducer/plugins/PFBadHcalPseudoClusterProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/PFBadHcalPseudoClusterProducer.cc new file mode 100644 index 0000000000000..6d599b24502c8 --- /dev/null +++ b/RecoParticleFlow/PFClusterProducer/plugins/PFBadHcalPseudoClusterProducer.cc @@ -0,0 +1,173 @@ +// system include files +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/ParticleFlowReco/interface/PFCluster.h" +#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h" +#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" + +#include "DataFormats/HcalDetId/interface/HcalDetId.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 "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/CaloTopology/interface/HcalTopology.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" + +#include "FWCore/Framework/interface/ESHandle.h" + + +using namespace std; +using namespace edm; + +// +// class declaration +// + +class PFBadHcalPseudoClusterProducer : public edm::stream::EDProducer<> { + public: + explicit PFBadHcalPseudoClusterProducer(const edm::ParameterSet&); + ~PFBadHcalPseudoClusterProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions & descriptions); + + private: + virtual void init(const EventSetup& c) ; + void produce(edm::Event&, const edm::EventSetup&) override; + + bool enabled_; + bool debug_; + + edm::ESHandle hQuality_; + edm::ESHandle hGeom_; + unsigned long long cacheId_quality_, cacheId_geom_; + std::vector badAreasC_; + std::vector badAreasRH_; +}; + + +PFBadHcalPseudoClusterProducer::PFBadHcalPseudoClusterProducer(const edm::ParameterSet& ps) : + enabled_(ps.getParameter("enable")), + debug_(ps.getUntrackedParameter("debug",false)), + cacheId_quality_(0), cacheId_geom_(0) +{ + produces>(); + produces>("hits"); +} + + +PFBadHcalPseudoClusterProducer::~PFBadHcalPseudoClusterProducer() +{ +} + +void PFBadHcalPseudoClusterProducer::init(const EventSetup& iSetup) +{ + badAreasC_.clear(); + badAreasRH_.clear(); + + edm::ESHandle hQuality_; + iSetup.get().get("withTopo",hQuality_); + const HcalChannelQuality & chanquality = *hQuality_; + + edm::ESHandle hGeom_; + iSetup.get().get(hGeom_); + const CaloGeometry& caloGeom = *hGeom_; + const CaloSubdetectorGeometry *hbGeom = caloGeom.getSubdetectorGeometry(DetId::Hcal, HcalBarrel); + const CaloSubdetectorGeometry *heGeom = caloGeom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap); + + int statusMask = ((1<, int> good, bads; + std::map, std::pair> minDepths; + for (const DetId & i : chanquality.getAllChannels()) { + if (i.det()!=DetId::Hcal) continue; // not an hcal cell + HcalDetId id = HcalDetId(i); + if (id.subdet() != HcalBarrel && id.subdet() != HcalEndcap) continue; // we don't deal with HO and HF here + int status = chanquality.getValues(id)->getValue(); + auto tower = std::make_pair(id.ieta(), id.iphi()); + if (status & statusMask) { + bads[tower]++; + if (debug_) std::cout << "Channel " << i() << " (subdet " << id.subdet() << ", zside " << id.zside() << ", ieta " << id.ieta() << ", iphi " << id.iphi() << " depth " << id.depth() << " has status " << status << " masked " << (status & statusMask) << std::endl; + } else { + good[tower]++; + } + auto & minD = minDepths[tower]; + if (minD.second == HcalEmpty || minD.first > id.depth()) { + minD.first = id.depth(); + minD.second = id.subdet(); + } + } + + const float dummyEnergy = 1e-5; // non-zero, but small (even if it shouldn't ever be used) + for (const auto & rec : bads) { + int ieta = rec.first.first, iphi = rec.first.second, nbad = rec.second, ngood = good[rec.first]; + auto minDepth = minDepths[rec.first]; + bool barrel = minDepth.second == HcalBarrel; + HcalDetId id(minDepth.second, ieta, iphi, minDepth.first); + bool isBad = (nbad > 0 && nbad >= ngood); + if (debug_) std::cout << "At ieta " << id.ieta() << ", iphi " << id.iphi() << " I have " << nbad << " bad depths, " << ngood << " good depths. First depth is in " << (barrel ? "HB" : "HE") << " depth " << minDepth.first <<"; " << (isBad ? " MARK BAD": " ignore") << std::endl; + if (!isBad) continue; + PFLayer::Layer layer = (barrel ? PFLayer::HCAL_BARREL1 : PFLayer::HCAL_ENDCAP); + // make a PF RecHit + std::shared_ptr thisCell = (barrel ? hbGeom : heGeom)->getGeometry(id); + const GlobalPoint & pos = thisCell->getPosition(); + badAreasRH_.emplace_back( thisCell, id(), layer, dummyEnergy ); + // make a PF cluster (but can't add the rechit, as for that I need an edm::Ref) + badAreasC_.emplace_back( layer, dummyEnergy, pos.x(), pos.y(), pos.z() ); + badAreasC_.back().setFlags(reco::CaloCluster::badHcalMarker); + } + + cacheId_quality_ = iSetup.get().cacheIdentifier(); + cacheId_geom_ = iSetup.get().cacheIdentifier(); +} + +// ------------ method called to produce the data ------------ + void +PFBadHcalPseudoClusterProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + if (enabled_) { + if (cacheId_quality_ != iSetup.get().cacheIdentifier() || + cacheId_geom_ != iSetup.get().cacheIdentifier()) { + init(iSetup); + } + } + + + auto outRH = std::make_unique>(badAreasRH_); + auto rhHandle = iEvent.put(std::move(outRH), "hits"); + + auto outC = std::make_unique>(badAreasC_); + + // now we go set references + for (unsigned int i = 0, n = rhHandle->size(); i < n; ++i) { + (*outC)[i].addRecHitFraction( reco::PFRecHitFraction(reco::PFRecHitRef(rhHandle,i), 1.0) ); + } + + iEvent.put(std::move(outC)); +} + +void PFBadHcalPseudoClusterProducer::fillDescriptions(edm::ConfigurationDescriptions & descriptions) { + edm::ParameterSetDescription desc; + desc.add("enable", false)->setComment("activate the module (if false, it doesn't check the DB and produces an empty collection)"); + desc.addUntracked("debug", false); + descriptions.add("particleFlowBadHcalPseudoCluster", desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(PFBadHcalPseudoClusterProducer); + diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowBadHcalPseudoCluster_cff.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowBadHcalPseudoCluster_cff.py new file mode 100644 index 0000000000000..028b5381922ea --- /dev/null +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowBadHcalPseudoCluster_cff.py @@ -0,0 +1,6 @@ +from RecoParticleFlow.PFClusterProducer.particleFlowBadHcalPseudoCluster_cfi import * + +# OFF by default, turned on via modifier +from Configuration.Eras.Modifier_pf_badHcalMitigation_cff import pf_badHcalMitigation +pf_badHcalMitigation.toModify(particleFlowBadHcalPseudoCluster, enable = True) + diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py index bde692c9f6483..b855fe46a3548 100644 --- a/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py @@ -19,6 +19,7 @@ from RecoParticleFlow.PFClusterProducer.particleFlowClusterHCAL_cfi import * from RecoParticleFlow.PFClusterProducer.particleFlowClusterHO_cfi import * from RecoParticleFlow.PFClusterProducer.particleFlowClusterPS_cfi import * +from RecoParticleFlow.PFClusterProducer.particleFlowBadHcalPseudoCluster_cff import * particleFlowClusterECALTask = cms.Task(particleFlowClusterECAL) particleFlowClusterECALSequence = cms.Sequence(particleFlowClusterECALTask) @@ -40,6 +41,7 @@ pfClusteringHO = cms.Sequence(pfClusteringHOTask) particleFlowClusterWithoutHOTask = cms.Sequence( + particleFlowBadHcalPseudoCluster, pfClusteringPSTask, pfClusteringECALTask, pfClusteringHBHEHFTask @@ -47,6 +49,7 @@ particleFlowClusterWithoutHO = cms.Sequence(particleFlowClusterWithoutHOTask) particleFlowClusterTask = cms.Task( + particleFlowBadHcalPseudoCluster, pfClusteringPSTask, pfClusteringECALTask, pfClusteringHBHEHFTask, diff --git a/RecoParticleFlow/PFProducer/interface/PFAlgo.h b/RecoParticleFlow/PFProducer/interface/PFAlgo.h index eef463198d665..35d05f40cd2b6 100644 --- a/RecoParticleFlow/PFProducer/interface/PFAlgo.h +++ b/RecoParticleFlow/PFProducer/interface/PFAlgo.h @@ -63,7 +63,7 @@ class PFAlgo { void setAlgo( int algo ) {algo_ = algo;} void setPFMuonAlgo(PFMuonAlgo* algo) {pfmu_ =algo;} void setMuonHandle(const edm::Handle&); - void setDebug( bool debug ) {debug_ = debug; connector_.setDebug(debug_);} + void setDebug( bool debug ) {debug_ = debug; connector_.setDebug(debug_); if (pfegamma_) pfegamma_->setDebug(debug); } void setParameters(double nSigmaECAL, double nSigmaHCAL, @@ -85,7 +85,9 @@ class PFAlgo { void setPFMuonAndFakeParameters(const edm::ParameterSet& pset); - + + void setBadHcalTrackParams(const edm::ParameterSet& pset); + PFMuonAlgo* getPFMuonAlgo(); void setPFEleParameters(double mvaEleCut, @@ -125,12 +127,14 @@ class PFAlgo { unsigned int ele_missinghits, bool useProtectionsForJetMET, const edm::ParameterSet& ele_protectionsForJetMET, + const edm::ParameterSet& ele_protectionsForBadHcal, double ph_MinEt, double ph_combIso, double ph_HoE, double ph_sietaieta_eb, double ph_sietaieta_ee, - const edm::ParameterSet& ph_protectionsForJetMET); + const edm::ParameterSet& ph_protectionsForJetMET, + const edm::ParameterSet& ph_protectionsForBadHcal); void setEGammaCollections(const edm::View & pfEgammaCandidates, @@ -396,6 +400,23 @@ class PFAlgo { double ptError_; std::vector factors45_; + /// Variables for track cleaning in bad HCal areas + float goodTrackDeadHcal_ptErrRel_; + float goodTrackDeadHcal_chi2n_; + int goodTrackDeadHcal_layers_; + float goodTrackDeadHcal_validFr_; + float goodTrackDeadHcal_dxy_; + + float goodPixelTrackDeadHcal_minEta_; + float goodPixelTrackDeadHcal_maxPt_; + float goodPixelTrackDeadHcal_ptErrRel_; + float goodPixelTrackDeadHcal_chi2n_; + int goodPixelTrackDeadHcal_maxLost3Hit_; + int goodPixelTrackDeadHcal_maxLost4Hit_; + float goodPixelTrackDeadHcal_dxy_; + float goodPixelTrackDeadHcal_dz_; + + // Parameters for post HF cleaning bool postHFCleaning_; bool postMuonCleaning_; diff --git a/RecoParticleFlow/PFProducer/interface/PFEGammaFilters.h b/RecoParticleFlow/PFProducer/interface/PFEGammaFilters.h index 1aac6ea98a9d3..e59cd7938b85a 100644 --- a/RecoParticleFlow/PFProducer/interface/PFEGammaFilters.h +++ b/RecoParticleFlow/PFProducer/interface/PFEGammaFilters.h @@ -23,6 +23,7 @@ class PFEGammaFilters { float ph_sietaieta_eb, float ph_sietaieta_ee, const edm::ParameterSet& ph_protectionsForJetMET, + const edm::ParameterSet& ph_protectionsForBadHcal, float ele_iso_pt, float ele_iso_mva_eb, float ele_iso_mva_ee, @@ -31,7 +32,8 @@ class PFEGammaFilters { float ele_noniso_mva, unsigned int ele_missinghits, const std::string& ele_iso_path_mvaWeightFile, - const edm::ParameterSet& ele_protectionsForJetMET + const edm::ParameterSet& ele_protectionsForJetMET, + const edm::ParameterSet& ele_protectionsForBadHcal ); ~PFEGammaFilters(){}; @@ -50,6 +52,7 @@ class PFEGammaFilters { bool isPhotonSafeForJetMET(const reco::Photon &, const reco::PFCandidate &); + void setDebug( bool debug ) { debug_ = debug; } private: @@ -78,10 +81,21 @@ class PFEGammaFilters { ele_maxEleHcalEOverEcalE, ele_maxEcalEOverPRes, ele_maxEeleOverPoutRes, ele_maxHcalEOverP, ele_maxHcalEOverEcalE, ele_maxEcalEOverP_1, ele_maxEcalEOverP_2, ele_maxEeleOverPout, ele_maxDPhiIN; - - + + // dead hcal selections (electrons) + std::array badHcal_full5x5_sigmaIetaIeta_; + std::array badHcal_eInvPInv_; + std::array badHcal_dEta_; + std::array badHcal_dPhi_; + bool badHcal_eleEnable_; + static void readEBEEParams_(const edm::ParameterSet &pset, const std::string &name, std::array & out) ; + + // dead hcal selections (photons) + float badHcal_phoTrkSolidConeIso_offs_, badHcal_phoTrkSolidConeIso_slope_; + bool badHcal_phoEnable_; + // Event variables - + bool debug_; }; #endif diff --git a/RecoParticleFlow/PFProducer/plugins/PFEGammaProducer.cc b/RecoParticleFlow/PFProducer/plugins/PFEGammaProducer.cc index 2401d69ffe8d8..66bf6b09551b0 100644 --- a/RecoParticleFlow/PFProducer/plugins/PFEGammaProducer.cc +++ b/RecoParticleFlow/PFProducer/plugins/PFEGammaProducer.cc @@ -340,6 +340,7 @@ PFEGammaProducer::produce(edm::Event& iEvent, case reco::PFBlockElement::HFEM: case reco::PFBlockElement::HFHAD: case reco::PFBlockElement::HCAL: + if (elements[0].clusterRef()->flags() & reco::CaloCluster::badHcalMarker) continue; hcalBlockRefs.push_back( blockref ); singleEcalOrHcal = true; break; diff --git a/RecoParticleFlow/PFProducer/plugins/PFProducer.cc b/RecoParticleFlow/PFProducer/plugins/PFProducer.cc index 3dee18a5c1ba8..3d2c1287966b1 100644 --- a/RecoParticleFlow/PFProducer/plugins/PFProducer.cc +++ b/RecoParticleFlow/PFProducer/plugins/PFProducer.cc @@ -207,8 +207,7 @@ PFProducer::PFProducer(const edm::ParameterSet& iConfig) { double ph_MinEt(0.0), ph_combIso(0.0), ph_HoE(0.0), ph_sietaieta_eb(0.0),ph_sietaieta_ee(0.0); string ele_iso_mvaWeightFile(""), ele_iso_path_mvaWeightFile(""); - edm::ParameterSet ele_protectionsForJetMET,ph_protectionsForJetMET; - + edm::ParameterSet ele_protectionsForJetMET,ele_protectionsForBadHcal,ph_protectionsForJetMET,ph_protectionsForBadHcal; // Reading new EGamma ubiased collections and value maps if(use_EGammaFilters_) { ele_iso_mvaWeightFile = iConfig.getParameter("isolatedElectronID_mvaWeightFile"); @@ -232,8 +231,12 @@ PFProducer::PFProducer(const edm::ParameterSet& iConfig) { iConfig.getParameter("useProtectionsForJetMET"); ele_protectionsForJetMET = iConfig.getParameter("electron_protectionsForJetMET"); + ele_protectionsForBadHcal = + iConfig.getParameter("electron_protectionsForBadHcal"); ph_protectionsForJetMET = iConfig.getParameter("photon_protectionsForJetMET"); + ph_protectionsForBadHcal = + iConfig.getParameter("photon_protectionsForBadHcal"); } //Secondary tracks and displaced vertices parameters @@ -330,12 +333,14 @@ PFProducer::PFProducer(const edm::ParameterSet& iConfig) { ele_missinghits, useProtectionsForJetMET, ele_protectionsForJetMET, + ele_protectionsForBadHcal, ph_MinEt, ph_combIso, ph_HoE, ph_sietaieta_eb, ph_sietaieta_ee, - ph_protectionsForJetMET); + ph_protectionsForJetMET, + ph_protectionsForBadHcal); //Secondary tracks and displaced vertices parameters @@ -353,6 +358,7 @@ PFProducer::PFProducer(const edm::ParameterSet& iConfig) { // Set muon and fake track parameters pfAlgo_->setPFMuonAndFakeParameters(iConfig); + pfAlgo_->setBadHcalTrackParams(iConfig); //Post cleaning of the HF postHFCleaning_ diff --git a/RecoParticleFlow/PFProducer/python/particleFlowBlock_cfi.py b/RecoParticleFlow/PFProducer/python/particleFlowBlock_cfi.py index 0a64ece7de28a..8e1cb6f247642 100644 --- a/RecoParticleFlow/PFProducer/python/particleFlowBlock_cfi.py +++ b/RecoParticleFlow/PFProducer/python/particleFlowBlock_cfi.py @@ -60,6 +60,8 @@ BCtoPFCMap = cms.InputTag('particleFlowSuperClusterECAL:PFClusterAssociationEBEE') ), cms.PSet( importerName = cms.string("GenericClusterImporter"), source = cms.InputTag("particleFlowClusterHCAL") ), + cms.PSet( importerName = cms.string("GenericClusterImporter"), + source = cms.InputTag("particleFlowBadHcalPseudoCluster") ), cms.PSet( importerName = cms.string("GenericClusterImporter"), source = cms.InputTag("particleFlowClusterHO") ), cms.PSet( importerName = cms.string("GenericClusterImporter"), diff --git a/RecoParticleFlow/PFProducer/python/particleFlow_cfi.py b/RecoParticleFlow/PFProducer/python/particleFlow_cfi.py index 5093cbab42e19..3780799dec45e 100644 --- a/RecoParticleFlow/PFProducer/python/particleFlow_cfi.py +++ b/RecoParticleFlow/PFProducer/python/particleFlow_cfi.py @@ -70,12 +70,24 @@ maxEeleOverPout = cms.double(0.2), maxDPhiIN = cms.double(0.1) ), + electron_protectionsForBadHcal = cms.PSet( + enableProtections = cms.bool(False), + full5x5_sigmaIetaIeta = cms.vdouble(0.0106, 0.0387), # EB, EE; 94Xv2 cut-based medium id + eInvPInv = cms.vdouble(0.184, 0.0721), + dEta = cms.vdouble(0.0032*2, 0.00632*2), # relax factor 2 to be safer against misalignment + dPhi = cms.vdouble(0.0547, 0.0394), + ), # New photon selection cuts for CMSSW_700 photon_MinEt = cms.double(10.), photon_combIso = cms.double(10.), photon_HoE = cms.double(0.05), photon_SigmaiEtaiEta_barrel = cms.double(0.0125), photon_SigmaiEtaiEta_endcap = cms.double(0.034), + photon_protectionsForBadHcal = cms.PSet( + enableProtections = cms.bool(False), + solidConeTrkIsoOffset = cms.double(10.), + solidConeTrkIsoSlope = cms.double(0.3), + ), # sumPtTrackIso, sumPtTrackIsoSlope photon_protectionsForJetMET = cms.PSet( @@ -186,6 +198,23 @@ # Factors to be applied in the four and fifth steps to the pt error factors_45 = cms.vdouble(10.,100.), + # Treatment of tracks in region of bad HCal + goodTrackDeadHcal_ptErrRel = cms.double(0.2), # trackRef->ptError()/trackRef->pt() < X + goodTrackDeadHcal_chi2n = cms.double(5), # trackRef->normalizedChi2() < X + goodTrackDeadHcal_layers = cms.uint32(4), # trackRef->hitPattern().trackerLayersWithMeasurement() >= X + goodTrackDeadHcal_validFr = cms.double(0.5), # trackRef->validFraction() > X + goodTrackDeadHcal_dxy = cms.double(0.5), # [cm] abs(trackRef->dxy(primaryVertex_.position())) < X + + goodPixelTrackDeadHcal_minEta = cms.double(2.3), # abs(trackRef->eta()) > X + goodPixelTrackDeadHcal_maxPt = cms.double(50.), # trackRef->ptError()/trackRef->pt() < X + goodPixelTrackDeadHcal_ptErrRel = cms.double(1.0), # trackRef->ptError()/trackRef->pt() < X + goodPixelTrackDeadHcal_chi2n = cms.double(2), # trackRef->normalizedChi2() < X + 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) + goodPixelTrackDeadHcal_maxLost4Hit = cms.int32(1), # max missing outer hits for a track with >= 4 valid pixel layers + goodPixelTrackDeadHcal_dxy = cms.double(0.02), # [cm] abs(trackRef->dxy(primaryVertex_.position())) < X + goodPixelTrackDeadHcal_dz = cms.double(0.05), # [cm] abs(trackRef->dz(primaryVertex_.position())) < X + + # Post HF cleaning postHFCleaning = cms.bool(False), # Clean only objects with pt larger than this value @@ -244,3 +273,8 @@ +from Configuration.Eras.Modifier_pf_badHcalMitigation_cff import pf_badHcalMitigation +pf_badHcalMitigation.toModify(particleFlowTmp, + electron_protectionsForBadHcal = dict(enableProtections = True), + photon_protectionsForBadHcal = dict(enableProtections = True)) + diff --git a/RecoParticleFlow/PFProducer/src/PFAlgo.cc b/RecoParticleFlow/PFProducer/src/PFAlgo.cc index e14de09b0e4d2..186550cc48d16 100644 --- a/RecoParticleFlow/PFProducer/src/PFAlgo.cc +++ b/RecoParticleFlow/PFProducer/src/PFAlgo.cc @@ -220,12 +220,14 @@ void PFAlgo::setEGammaParameters(bool use_EGammaFilters, unsigned int ele_missinghits, bool useProtectionsForJetMET, const edm::ParameterSet& ele_protectionsForJetMET, + const edm::ParameterSet& ele_protectionsForBadHcal, double ph_MinEt, double ph_combIso, double ph_HoE, double ph_sietaieta_eb, double ph_sietaieta_ee, - const edm::ParameterSet& ph_protectionsForJetMET + const edm::ParameterSet& ph_protectionsForJetMET, + const edm::ParameterSet& ph_protectionsForBadHcal ) { @@ -252,6 +254,7 @@ void PFAlgo::setEGammaParameters(bool use_EGammaFilters, ph_sietaieta_eb, ph_sietaieta_ee, ph_protectionsForJetMET, + ph_protectionsForBadHcal, ele_iso_pt, ele_iso_mva_barrel, ele_iso_mva_endcap, @@ -260,7 +263,8 @@ void PFAlgo::setEGammaParameters(bool use_EGammaFilters, ele_noniso_mva, ele_missinghits, ele_iso_path_mvaWeightFile, - ele_protectionsForJetMET); + ele_protectionsForJetMET, + ele_protectionsForBadHcal); return; } @@ -326,6 +330,24 @@ PFAlgo::setPFMuonAndFakeParameters(const edm::ParameterSet& pset) } +void +PFAlgo::setBadHcalTrackParams(const edm::ParameterSet& pset) +{ + goodTrackDeadHcal_ptErrRel_ = pset.getParameter("goodTrackDeadHcal_ptErrRel"); + goodTrackDeadHcal_chi2n_ = pset.getParameter("goodTrackDeadHcal_chi2n"); + goodTrackDeadHcal_layers_ = pset.getParameter("goodTrackDeadHcal_layers"); + goodTrackDeadHcal_validFr_ = pset.getParameter("goodTrackDeadHcal_validFr"); + goodTrackDeadHcal_dxy_ = pset.getParameter("goodTrackDeadHcal_dxy"); + + goodPixelTrackDeadHcal_minEta_ = pset.getParameter("goodPixelTrackDeadHcal_minEta"); + goodPixelTrackDeadHcal_maxPt_ = pset.getParameter("goodPixelTrackDeadHcal_maxPt"); + goodPixelTrackDeadHcal_ptErrRel_ = pset.getParameter("goodPixelTrackDeadHcal_ptErrRel"); + goodPixelTrackDeadHcal_chi2n_ = pset.getParameter("goodPixelTrackDeadHcal_chi2n"); + goodPixelTrackDeadHcal_maxLost3Hit_ = pset.getParameter("goodPixelTrackDeadHcal_maxLost3Hit"); + goodPixelTrackDeadHcal_maxLost4Hit_ = pset.getParameter("goodPixelTrackDeadHcal_maxLost4Hit"); + goodPixelTrackDeadHcal_dxy_ = pset.getParameter("goodPixelTrackDeadHcal_dxy"); + goodPixelTrackDeadHcal_dz_ = pset.getParameter("goodPixelTrackDeadHcal_dz"); +} void PFAlgo::setMuonHandle(const edm::Handle& muons) { @@ -632,7 +654,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, if(useEGammaFilters_) { // const edm::ValueMap & myGedElectronValMap(*valueMapGedElectrons_); - bool egmLocalDebug = false; + bool egmLocalDebug = debug_; bool egmLocalBlockDebug = false; unsigned int negmcandidates = pfEgammaCandidates_->size(); @@ -666,7 +688,11 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, cout << "** Good Electron, pt " << gedEleRef->pt() << " eta, phi " << gedEleRef->eta() << ", " << gedEleRef->phi() << " charge " << gedEleRef->charge() - << " isPrimary " << isPrimaryElectron << endl; + << " isPrimary " << isPrimaryElectron + << " isoDr03 " << (gedEleRef->dr03TkSumPt() + gedEleRef->dr03EcalRecHitSumEt() + gedEleRef->dr03HcalTowerSumEt()) + << " mva_isolated " << gedEleRef->mva_Isolated() + << " mva_e_pi " << gedEleRef->mva_e_pi() + << endl; } @@ -728,7 +754,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, for (PFCandidate::ElementsInBlocks::const_iterator ieb = theElements.begin(); iebsecond] = false; - if(egmLocalBlockDebug) + if(egmLocalBlockDebug||(debug_&&egmLocalDebug)) cout << " Elements used " << ieb->second << endl; } @@ -737,6 +763,8 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, const PFCandidate::ElementsInBlocks& extraTracks = myPFElectron.egammaExtraRef()->extraNonConvTracks(); for (PFCandidate::ElementsInBlocks::const_iterator itrk = extraTracks.begin(); itrksecond << endl; active[itrk->second] = false; } } @@ -780,7 +808,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, for (PFCandidate::ElementsInBlocks::const_iterator ieb = theElements.begin(); iebsecond] = false; - if(egmLocalBlockDebug) + if(egmLocalBlockDebug||(debug_&&egmLocalDebug)) cout << " Elements used " << ieb->second << endl; } pfCandidates_->push_back(myPFPhoton); @@ -850,7 +878,8 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, vector hfEmIs; vector hfHadIs; - + + vector deadArea(elements.size(), false); for(unsigned iEle=0; iEleflags() & reco::CaloCluster::badHcalMarker) { + if(debug_) cout<<"HCAL DEAD AREA: remember and skip."<second]) { hasDeadHcal = true; it = hcalElems.erase(it); } // std::multimap::erase returns iterator to next + // else ++it; + //} + //--- option (2) -- + bool hasDeadHcal = false; + if (!hcalElems.empty() && deadArea[hcalElems.begin()->second]) { + hasDeadHcal = true; + hcalElems.clear(); + } + //--- option (3) -- + //bool hasDeadHcal = true; + //for (auto it = hcalElems.begin(), ed = hcalElems.end(); it != ed; /*NOTE NO ++it HERE */ ) { + // if (deadArea[it->second]) { it = hcalElems.erase(it); } // std::multimap::erase returns iterator to next + // else { hasDeadHcal = false; } + //} + + // for tracks with bad Hcal, check the quality + bool goodTrackDeadHcal = false; + if (hasDeadHcal) { + goodTrackDeadHcal = ( trackRef->ptError() < goodTrackDeadHcal_ptErrRel_ * trackRef->pt() && + trackRef->normalizedChi2() < goodTrackDeadHcal_chi2n_ && + trackRef->hitPattern().trackerLayersWithMeasurement() >= goodTrackDeadHcal_layers_ && + trackRef->validFraction() > goodTrackDeadHcal_validFr_ && + std::abs(trackRef->dxy(primaryVertex_.position())) < goodTrackDeadHcal_dxy_ ); + // now we add an extra block for tracks at high |eta| + if (!goodTrackDeadHcal && + std::abs(trackRef->eta()) > goodPixelTrackDeadHcal_minEta_ && // high eta + trackRef->hitPattern().pixelLayersWithMeasurement() >= 3 && // pixel track + trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS) == 0 && + trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS) == 0 && + trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS) <= ( + trackRef->hitPattern().pixelLayersWithMeasurement() > 3 ? + goodPixelTrackDeadHcal_maxLost4Hit_ : goodPixelTrackDeadHcal_maxLost3Hit_ ) && + trackRef->normalizedChi2() < goodPixelTrackDeadHcal_chi2n_ && // tighter cut + std::abs(trackRef->dxy(primaryVertex_.position())) < goodPixelTrackDeadHcal_dxy_ && + std::abs(trackRef->dz(primaryVertex_.position())) < goodPixelTrackDeadHcal_dz_ && + trackRef->ptError() < goodPixelTrackDeadHcal_ptErrRel_*trackRef->pt() && // sanity + trackRef->pt() < goodPixelTrackDeadHcal_maxPt_) { // sanity + goodTrackDeadHcal = true; + // FIXME: may decide to do something to the track pT + } + //if (!goodTrackDeadHcal && trackRef->hitPattern().trackerLayersWithMeasurement() == 4 && trackRef->validFraction() == 1 + if (debug_) cout << " track pt " << trackRef->pt() << " +- " << trackRef->ptError() + << " layers valid " << trackRef->hitPattern().trackerLayersWithMeasurement() + << ", lost " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS) + << ", lost outer " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS) + << ", lost inner " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS) + << "(valid fraction " << trackRef->validFraction() << ")" + << " chi2/ndf " << trackRef->normalizedChi2() + << " |dxy| " << std::abs(trackRef->dxy(primaryVertex_.position())) << " +- " << trackRef->dxyError() + << " |dz| " << std::abs(trackRef->dz(primaryVertex_.position())) << " +- " << trackRef->dzError() + << (goodTrackDeadHcal ? " passes " : " fails ") << "quality cuts" << std::endl; + } + // When a track has no HCAL cluster linked, but another track is linked to the same // ECAL cluster and an HCAL cluster, link the track to the HCAL cluster for // later analysis - if ( hcalElems.empty() && !ecalElems.empty() ) { + if ( hcalElems.empty() && !ecalElems.empty() && !hasDeadHcal) { // debug_ = true; unsigned ntt = 1; unsigned index = ecalElems.begin()->second; @@ -963,12 +1076,12 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL ); - // std::cout << "The ECAL is linked to " << sortedTracks.size() << std::endl; + if (debug_ ) std::cout << "The closest ECAL cluster is linked to " << sortedTracks.size() << " tracks, with distance = " << ecalElems.begin()->first << std::endl; // Loop over all tracks for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) { unsigned jTrack = ie->second; - // std::cout << "Track " << jTrack << std::endl; + //std::cout << "Track " << jTrack << std::endl; // Track must be active if ( !active[jTrack] ) continue; //std::cout << "Active " << std::endl; @@ -985,7 +1098,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL ); if ( sortedECAL.begin()->second != index ) continue; - //std::cout << "With closest ECAL identical " << std::endl; + if (debug_ ) std::cout << " track " << jTrack << " with closest ECAL identical " << std::endl; // Check if this track is also linked to an HCAL @@ -995,7 +1108,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL ); if ( sortedHCAL.empty() ) continue; - //std::cout << "With an HCAL cluster " << sortedHCAL.begin()->first << std::endl; + if (debug_ ) std::cout << " and with an HCAL cluster " << sortedHCAL.begin()->second << std::endl; ntt++; // In that case establish a link with the first track @@ -1078,7 +1191,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, if( hcalElems.empty() ) { if ( debug_ ) - std::cout << "Now deals with tracks linked to no HCAL clusters" << std::endl; + std::cout << "Now deals with tracks linked to no HCAL clusters. Was HCal active? " << (!hasDeadHcal) << std::endl; // vector elementIndices; // elementIndices.push_back(iTrack); @@ -1106,19 +1219,40 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, if ( !ecalElems.empty() ) { unsigned thisEcal = ecalElems.begin()->second; reco::PFClusterRef clusterRef = elements[thisEcal].clusterRef(); - deficit -= clusterRef->energy(); - resol = neutralHadronEnergyResolution(trackMomentum, - clusterRef->positionREP().Eta()); - resol *= trackMomentum; + bool useCluster = true; + if (hasDeadHcal) { + std::multimap sortedTracks; + block.associatedElements( thisEcal, linkData, + sortedTracks, + reco::PFBlockElement::TRACK, + reco::PFBlock::LINKTEST_ALL ); + useCluster = (sortedTracks.begin()->second == iTrack); + } + if (useCluster) { + deficit -= clusterRef->energy(); + resol = neutralHadronEnergyResolution(trackMomentum, + clusterRef->positionREP().Eta()); + resol *= trackMomentum; + } } bool isPrimary = isFromSecInt(elements[iTrack], "primary"); - if ( deficit > nSigmaTRACK_*resol && !isPrimary) { + if ( deficit > nSigmaTRACK_*resol && !isPrimary && !goodTrackDeadHcal) { rejectFake = true; active[iTrack] = false; if ( debug_ ) std::cout << elements[iTrack] << std::endl + << " deficit " << deficit << " > " << nSigmaTRACK_ << " x " << resol + << " track pt " << trackRef->pt() << " +- " << trackRef->ptError() + << " layers valid " << trackRef->hitPattern().trackerLayersWithMeasurement() + << ", lost " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS) + << ", lost outer " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS) + << ", lost inner " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS) + << "(valid fraction " << trackRef->validFraction() <<")" + << " chi2/ndf " << trackRef->normalizedChi2() + << " |dxy| " << std::abs(trackRef->dxy(primaryVertex_.position())) << " +- " << trackRef->dxyError() + << " |dz| " << std::abs(trackRef->dz(primaryVertex_.position())) << " +- " << trackRef->dzError() << "is probably a fake (1) --> lock the track" << std::endl; } @@ -1137,7 +1271,8 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, if ( rejectTracks_Step45_ && ecalElems.empty() && trackMomentum > 30. && Dpt > 0.5 && - ( PFTrackAlgoTools::step45(trackRef->algo())) ) { + ( PFTrackAlgoTools::step45(trackRef->algo()) + && !goodTrackDeadHcal) ) { double dptRel = Dpt/trackRef->pt()*100; bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all"); @@ -1203,7 +1338,26 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL ); - + if (debug_) cout << "the closest ECAL cluster, id " << thisEcal << ", has " << sortedTracks.size() << " associated tracks, now processing them. " << endl; + + if (hasDeadHcal && !sortedTracks.empty()) { + // GP: only allow the ecal cluster closest to this track to be considered + if (debug_) cout << " the closest track to ECAL " << thisEcal << " is " << sortedTracks.begin()->second << " (distance " << sortedTracks.begin()->first << ")" << endl; + if (sortedTracks.begin()->second != iTrack) { + if (debug_) cout << " the closest track to ECAL " << thisEcal << " is " << sortedTracks.begin()->second << " which is not the one being processed. Will skip ECAL linking for this track" << endl; + (*pfCandidates_)[tmpi[0]].setEcalEnergy( 0., 0. ); + (*pfCandidates_)[tmpi[0]].setHcalEnergy( 0., 0. ); + (*pfCandidates_)[tmpi[0]].setHoEnergy( 0., 0. ); + (*pfCandidates_)[tmpi[0]].setPs1Energy( 0 ); + (*pfCandidates_)[tmpi[0]].setPs2Energy( 0 ); + (*pfCandidates_)[tmpi[0]].addElementInBlock( blockref, kTrack[0] ); + continue; + } else { + if (debug_) cout << " the closest track to ECAL " << thisEcal << " is " << sortedTracks.begin()->second << " which is the one being processed. Will skip ECAL linking for all other track" << endl; + sortedTracks.clear(); + } + } + for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) { unsigned jTrack = ie->second; @@ -1216,6 +1370,9 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, // Check if the ECAL cluster closest to this track is // indeed the current ECAL cluster. Only tracks not // linked to an HCAL cluster are considered here + // (GP: this is because if there's a jTrack linked + // to the same Ecal cluster and that also has an Hcal link + // we would have also linked iTrack to that Hcal above) std::multimap sortedECAL; block.associatedElements( jTrack, linkData, sortedECAL, @@ -1225,23 +1382,31 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, // Check if this track is a muon bool thatIsAMuon = PFMuonAlgo::isMuon(elements[jTrack]); + if (debug_) { + cout << " found track " << jTrack << (thatIsAMuon ? " (muon) " : " (non-muon)") << ", with distance = " << sortedECAL.begin()->first << std::endl; + } // Check if this track is not a fake bool rejectFake = false; reco::TrackRef trackRef = elements[jTrack].trackRef(); if ( !thatIsAMuon && trackRef->ptError() > ptError_) { + // GP: FIXME this selection should be adapted depending on hasDeadHcal + // but we don't know if jTrack is linked to a dead Hcal or not double deficit = trackMomentum + trackRef->p() - clusterRef->energy(); double resol = nSigmaTRACK_*neutralHadronEnergyResolution(trackMomentum+trackRef->p(), clusterRef->positionREP().Eta()); resol *= (trackMomentum+trackRef->p()); - if ( deficit > nSigmaTRACK_*resol ) { + if ( deficit > nSigmaTRACK_*resol && !goodTrackDeadHcal) { rejectFake = true; kTrack.push_back(jTrack); active[jTrack] = false; if ( debug_ ) std::cout << elements[jTrack] << std::endl << "is probably a fake (2) --> lock the track" + << "(trackMomentum = " << trackMomentum << ", clusterEnergy = " << clusterRef->energy() << + ", deficit = " << deficit << " > " << nSigmaTRACK_ << " x " << resol << + " assuming neutralHadronEnergyResolution " << neutralHadronEnergyResolution(trackMomentum+trackRef->p(), clusterRef->positionREP().Eta()) << ") " << std::endl; } } @@ -1358,7 +1523,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, if ( totalEcal > 0.) slopeEcal = calibEcal/totalEcal; if ( debug_ ) - std::cout << "The total calibrated energy so far amounts to = " << calibEcal << std::endl; + std::cout << "The total calibrated energy so far amounts to = " << calibEcal << " (slope = " << slopeEcal << ")" << std::endl; // Stop the loop when adding more ECAL clusters ruins the compatibility @@ -1501,7 +1666,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, } } - } + } // end if( hcalElems.empty() ) // In case several HCAL elements are linked to this track, @@ -2262,7 +2427,7 @@ void PFAlgo::processBlock( const reco::PFBlockRef& blockref, Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError); } } - /* + if(debug_){ cout<<"\tBefore Cleaning "< Calo Energy = "< Calo Energy- total charged momentum = " - < hadron "< check if the excess is photonic or hadronic"< check if the excess is photonic or hadronic"<("maxEcalEOverP_2")), ele_maxEeleOverPout(ele_protectionsForJetMET.getParameter("maxEeleOverPout")), - ele_maxDPhiIN(ele_protectionsForJetMET.getParameter("maxDPhiIN")) + ele_maxDPhiIN(ele_protectionsForJetMET.getParameter("maxDPhiIN")), + badHcal_eleEnable_(ele_protectionsForBadHcal.getParameter("enableProtections")), + badHcal_phoTrkSolidConeIso_offs_(ph_protectionsForBadHcal.getParameter("solidConeTrkIsoOffset")), + badHcal_phoTrkSolidConeIso_slope_(ph_protectionsForBadHcal.getParameter("solidConeTrkIsoSlope")), + badHcal_phoEnable_(ph_protectionsForBadHcal.getParameter("enableProtections")), + debug_(false) { + readEBEEParams_(ele_protectionsForBadHcal, "full5x5_sigmaIetaIeta", badHcal_full5x5_sigmaIetaIeta_); + readEBEEParams_(ele_protectionsForBadHcal, "eInvPInv", badHcal_eInvPInv_); + readEBEEParams_(ele_protectionsForBadHcal, "dEta", badHcal_dEta_); + readEBEEParams_(ele_protectionsForBadHcal, "dPhi", badHcal_dPhi_); } +void PFEGammaFilters::readEBEEParams_(const edm::ParameterSet &pset, const std::string &name, std::array & out) { + const auto & vals = pset.getParameter>(name); + if (vals.size() != 2) throw cms::Exception("Configuration") << "Parameter " << name << " does not contain exactly 2 values (EB, EE)\n"; + out[0] = vals[0]; + out[1] = vals[1]; +} + + bool PFEGammaFilters::passPhotonSelection(const reco::Photon & photon) { // First simple selection, same as the Run1 to be improved in CMSSW_710 // Photon ET if(photon.pt() < ph_Et_ ) return false; -// std::cout<< "Cuts " << ph_combIso_ << " H/E " << ph_loose_hoe_ -// << " SigmaiEtaiEta_EB " << ph_sietaieta_eb_ -// << " SigmaiEtaiEta_EE " << ph_sietaieta_ee_ << std::endl; + bool validHoverE = photon.hadTowOverEmValid(); + if (debug_) std::cout<< "PFEGammaFilters:: photon pt " << photon.pt() + << " eta, phi " << photon.eta() << ", " << photon.phi() + << " isoDr03 " << (photon.trkSumPtHollowConeDR03()+photon.ecalRecHitSumEtConeDR03()+photon.hcalTowerSumEtConeDR03()) << " (cut: " << ph_combIso_ << ")" + << " H/E " << photon.hadTowOverEm() << " (valid? " << validHoverE << ", cut: " << ph_loose_hoe_ << ")" + << " s(ieie) " << photon.sigmaIetaIeta() << " (cut: " << (photon.isEB() ? ph_sietaieta_eb_ : ph_sietaieta_ee_) << ")" + << " isoTrkDr03Solid " << (photon.trkSumPtSolidConeDR03()) << " (cut: " << ( + validHoverE || !badHcal_phoEnable_ ? + -1 : + badHcal_phoTrkSolidConeIso_offs_ + badHcal_phoTrkSolidConeIso_slope_*photon.pt()) << ")" + << std::endl; if (photon.hadTowOverEm() >ph_loose_hoe_ ) return false; //Isolation variables in 0.3 cone combined if(photon.trkSumPtHollowConeDR03()+photon.ecalRecHitSumEtConeDR03()+photon.hcalTowerSumEtConeDR03() > ph_combIso_) - return false; + return false; + + //patch for bad hcal + if (!validHoverE && badHcal_phoEnable_ && photon.trkSumPtSolidConeDR03() > badHcal_phoTrkSolidConeIso_offs_ + badHcal_phoTrkSolidConeIso_slope_*photon.pt()) { + return false; + } if(photon.isEB()) { if(photon.sigmaIetaIeta() > ph_sietaieta_eb_) @@ -88,7 +120,22 @@ bool PFEGammaFilters::passElectronSelection(const reco::GsfElectron & electron, const reco::PFCandidate & pfcand, const int & nVtx) { // First simple selection, same as the Run1 to be improved in CMSSW_710 - + + bool validHoverE = electron.hcalOverEcalValid(); + if (debug_) std::cout << "PFEGammaFilters:: Electron pt " << electron.pt() + << " eta, phi " << electron.eta() << ", " << electron.phi() + << " charge " << electron.charge() + << " isoDr03 " << (electron.dr03TkSumPt() + electron.dr03EcalRecHitSumEt() + electron.dr03HcalTowerSumEt()) + << " mva_isolated " << electron.mva_Isolated() + << " mva_e_pi " << electron.mva_e_pi() + << " H/E_valid " << validHoverE + << " s(ieie) " << electron.full5x5_sigmaIetaIeta() + << " H/E " << electron.hcalOverEcal() + << " 1/e-1/p " << (1.0-electron.eSuperClusterOverP())/electron.ecalEnergy() + << " deta " << std::abs(electron.deltaEtaSeedClusterTrackAtVtx()) + << " dphi " << std::abs(electron.deltaPhiSuperClusterTrackAtVtx()) + << endl; + bool passEleSelection = false; // Electron ET @@ -111,7 +158,17 @@ bool PFEGammaFilters::passElectronSelection(const reco::GsfElectron & electron, // cout << " My OLD MVA " << pfcand.mva_e_pi() << " MyNEW MVA " << electron.mva() << endl; if(electron.mva_e_pi() > ele_noniso_mva_) { - passEleSelection = true; + if (validHoverE || !badHcal_eleEnable_) { + passEleSelection = true; + } else { + bool EE = (std::abs(electron.eta()) > 1.485); // for prefer consistency with above than with E/gamma for now + if ((electron.full5x5_sigmaIetaIeta() < badHcal_full5x5_sigmaIetaIeta_[EE]) && + (std::abs(1.0-electron.eSuperClusterOverP())/electron.ecalEnergy() < badHcal_eInvPInv_[EE]) && + (std::abs(electron.deltaEtaSeedClusterTrackAtVtx()) < badHcal_dEta_[EE]) && // looser in case of misalignment + (std::abs(electron.deltaPhiSuperClusterTrackAtVtx()) < badHcal_dPhi_[EE])) { + passEleSelection = true; + } + } } return passEleSelection;