diff --git a/CalibCalorimetry/EcalTPGTools/BuildFile.xml b/CalibCalorimetry/EcalTPGTools/BuildFile.xml index 3653c6a1d060d..462039b5defc7 100644 --- a/CalibCalorimetry/EcalTPGTools/BuildFile.xml +++ b/CalibCalorimetry/EcalTPGTools/BuildFile.xml @@ -6,6 +6,9 @@ + + + diff --git a/CalibCalorimetry/EcalTPGTools/interface/EcalReadoutTools.h b/CalibCalorimetry/EcalTPGTools/interface/EcalReadoutTools.h new file mode 100644 index 0000000000000..79529c8fb1cc3 --- /dev/null +++ b/CalibCalorimetry/EcalTPGTools/interface/EcalReadoutTools.h @@ -0,0 +1,28 @@ +#ifndef CalibCalorimetry_EcalTPGTools_EcalReadoutTools_H +#define CalibCalorimetry_EcalTPGTools_EcalReadoutTools_H + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h" +#include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h" +#include "Geometry/EcalMapping/interface/EcalMappingRcd.h" + +class EcalReadoutTools { + + private: + const EcalTrigTowerConstituentsMap * triggerTowerMap_; + const EcalElectronicsMapping* elecMap_; + + public: + EcalReadoutTools(const edm::Event &iEvent, const edm::EventSetup &iSetup); + EcalReadoutTools(const EcalReadoutTools&) = delete; + EcalReadoutTools& operator=(const EcalReadoutTools&) = delete; + + EcalTrigTowerDetId readOutUnitOf(const EBDetId& xtalId) const; + EcalScDetId readOutUnitOf(const EEDetId& xtalId) const; +}; + +#endif diff --git a/CalibCalorimetry/EcalTPGTools/src/EcalReadoutTools.cc b/CalibCalorimetry/EcalTPGTools/src/EcalReadoutTools.cc new file mode 100644 index 0000000000000..0c932666bfe34 --- /dev/null +++ b/CalibCalorimetry/EcalTPGTools/src/EcalReadoutTools.cc @@ -0,0 +1,26 @@ +#include "CalibCalorimetry/EcalTPGTools/interface/EcalReadoutTools.h" + +EcalReadoutTools::EcalReadoutTools(const edm::Event &iEvent, const edm::EventSetup &iSetup) { + + edm::ESHandle hTriggerTowerMap; + iSetup.get().get(hTriggerTowerMap); + triggerTowerMap_ = hTriggerTowerMap.product(); + + edm::ESHandle< EcalElectronicsMapping > ecalmapping; + iSetup.get< EcalMappingRcd >().get(ecalmapping); + elecMap_ = ecalmapping.product(); + +} + +EcalTrigTowerDetId EcalReadoutTools::readOutUnitOf(const EBDetId& xtalId) const{ + return triggerTowerMap_->towerOf(xtalId); +} + +EcalScDetId EcalReadoutTools::readOutUnitOf(const EEDetId& xtalId) const{ + const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId); + int iDCC= EcalElecId.dccId(); + int iDccChan = EcalElecId.towerId(); + const bool ignoreSingle = true; + const std::vector id = elecMap_->getEcalScDetId(iDCC, iDccChan, ignoreSingle); + return id.size()>0?id[0]:EcalScDetId(); +} diff --git a/Configuration/AlCa/python/autoCond.py b/Configuration/AlCa/python/autoCond.py index d46827133a8af..97b9a44314f84 100644 --- a/Configuration/AlCa/python/autoCond.py +++ b/Configuration/AlCa/python/autoCond.py @@ -2,60 +2,60 @@ ### NEW KEYS ### # GlobalTag for MC production with perfectly aligned and calibrated detector for Run1 - 'run1_design' : '92X_mcRun1_design_v1', + 'run1_design' : '92X_mcRun1_design_v2', # GlobalTag for MC production (pp collisions) with realistic alignment and calibrations for Run1 - 'run1_mc' : '92X_mcRun1_realistic_v1', + 'run1_mc' : '92X_mcRun1_realistic_v2', # GlobalTag for MC production (Heavy Ions collisions) with realistic alignment and calibrations for Run1 - 'run1_mc_hi' : '92X_mcRun1_HeavyIon_v1', + 'run1_mc_hi' : '92X_mcRun1_HeavyIon_v2', # GlobalTag for MC production (p-Pb collisions) with realistic alignment and calibrations for Run1 - 'run1_mc_pa' : '92X_mcRun1_pA_v1', + 'run1_mc_pa' : '92X_mcRun1_pA_v2', # GlobalTag for MC production with perfectly aligned and calibrated detector for Run2 - 'run2_design' : '92X_mcRun2_design_v1', + 'run2_design' : '92X_mcRun2_design_v2', # GlobalTag for MC production with pessimistic alignment and calibrations for Run2 - 'run2_mc_50ns' : '92X_mcRun2_startup_v1', + 'run2_mc_50ns' : '92X_mcRun2_startup_v2', #GlobalTag for MC production with optimistic alignment and calibrations for Run2 - 'run2_mc' : '92X_mcRun2_asymptotic_v1', + 'run2_mc' : '92X_mcRun2_asymptotic_v2', # GlobalTag for MC production (cosmics) with starup-like alignment and calibrations for Run2, Strip tracker in peak mode - 'run2_mc_cosmics' : '92X_mcRun2cosmics_startup_deco_v1', + 'run2_mc_cosmics' : '92X_mcRun2cosmics_startup_deco_v2', # GlobalTag for MC production (Heavy Ions collisions) with optimistic alignment and calibrations for Run2 - 'run2_mc_hi' : '92X_mcRun2_HeavyIon_v1', + 'run2_mc_hi' : '92X_mcRun2_HeavyIon_v2', # GlobalTag for MC production (p-Pb collisions) with realistic alignment and calibrations for Run2 - 'run2_mc_pa' : '92X_mcRun2_pA_v1', + 'run2_mc_pa' : '92X_mcRun2_pA_v2', # GlobalTag for Run1 data reprocessing - 'run1_data' : '92X_dataRun2_v1', + 'run1_data' : '92X_dataRun2_v2', # GlobalTag for Run2 data reprocessing - 'run2_data' : '92X_dataRun2_v1', + 'run2_data' : '92X_dataRun2_v2', # GlobalTag for Run2 data relvals: allows customization to run with fixed L1 menu - 'run2_data_relval' : '92X_dataRun2_relval_v1', + 'run2_data_relval' : '92X_dataRun2_relval_v2', # GlobalTag for Run2 data 2016H relvals only: Prompt Conditions + fixed L1 menu (to be removed) - 'run2_data_promptlike' : '92X_dataRun2_PromptLike_v0', + 'run2_data_promptlike' : '92X_dataRun2_PromptLike_v1', # GlobalTag for Run1 HLT: it points to the online GT - 'run1_hlt' : '92X_dataRun2_HLT_frozen_v1', + 'run1_hlt' : '92X_dataRun2_HLT_frozen_v2', # GlobalTag for Run2 HLT: it points to the online GT - 'run2_hlt' : '92X_dataRun2_HLT_frozen_v1', + 'run2_hlt' : '92X_dataRun2_HLT_frozen_v2', # GlobalTag for Run2 HLT RelVals: customizations to run with fixed L1 Menu - 'run2_hlt_relval' : '92X_dataRun2_HLT_relval_v3', + 'run2_hlt_relval' : '92X_dataRun2_HLT_relval_v4', # GlobalTag for Run2 HLT for HI: it points to the online GT - 'run2_hlt_hi' : '92X_dataRun2_HLTHI_frozen_v2', + 'run2_hlt_hi' : '92X_dataRun2_HLTHI_frozen_v3', # GlobalTag for MC production with perfectly aligned and calibrated detector for Phase1 2017 (and 0,0,~0-centred beamspot) - 'phase1_2017_design' : '92X_upgrade2017_design_IdealBS_v5', + 'phase1_2017_design' : '92X_upgrade2017_design_IdealBS_v6', # GlobalTag for MC production with realistic conditions for Phase1 2017 detector - 'phase1_2017_realistic' : '92X_upgrade2017_realistic_v5', + 'phase1_2017_realistic' : '92X_upgrade2017_realistic_v6', # GlobalTag for MC production (cosmics) with realistic alignment and calibrations for Phase1 2017 detector, Strip tracker in DECO mode - 'phase1_2017_cosmics' : '92X_upgrade2017cosmics_realistic_deco_v5', + 'phase1_2017_cosmics' : '92X_upgrade2017cosmics_realistic_deco_v6', # GlobalTag for MC production (cosmics) with realistic alignment and calibrations for Phase1 2017 detector, Strip tracker in PEAK mode - 'phase1_2017_cosmics_peak' : '92X_upgrade2017cosmics_realistic_peak_v5', + 'phase1_2017_cosmics_peak' : '92X_upgrade2017cosmics_realistic_peak_v6', # GlobalTag for MC production with perfectly aligned and calibrated detector for full Phase1 2018 (and 0,0,0-centred beamspot) - 'phase1_2018_design' : '92X_upgrade2018_design_IdealBS_v4', + 'phase1_2018_design' : '92X_upgrade2018_design_IdealBS_v5', # GlobalTag for MC production with realistic conditions for full Phase1 2018 detector - 'phase1_2018_realistic' : '92X_upgrade2018_realistic_v5', + 'phase1_2018_realistic' : '92X_upgrade2018_realistic_v6', # GlobalTag for MC production (cosmics) with realistic conditions for full Phase1 2018 detector, Strip tracker in DECO mode - 'phase1_2018_cosmics' : '92X_upgrade2018cosmics_realistic_deco_v5', + 'phase1_2018_cosmics' : '92X_upgrade2018cosmics_realistic_deco_v6', # GlobalTag for MC production with perfectly aligned and calibrated detector for Phase1 2019 'phase1_2019_design' : 'DES19_70_V2', # placeholder (GT not meant for standard RelVal) # GlobalTag for MC production with realistic conditions for Phase2 2023 - 'phase2_realistic' : '92X_upgrade2023_realistic_v0' + 'phase2_realistic' : '92X_upgrade2023_realistic_v1' } aliases = { diff --git a/RecoParticleFlow/PFClusterProducer/BuildFile.xml b/RecoParticleFlow/PFClusterProducer/BuildFile.xml index 4368b7d1d8212..e4617d7fed9ea 100644 --- a/RecoParticleFlow/PFClusterProducer/BuildFile.xml +++ b/RecoParticleFlow/PFClusterProducer/BuildFile.xml @@ -10,6 +10,7 @@ + diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFClusterEMEnergyCorrector.h b/RecoParticleFlow/PFClusterProducer/interface/PFClusterEMEnergyCorrector.h index bd5c30f1a9011..38cd633cafffa 100644 --- a/RecoParticleFlow/PFClusterProducer/interface/PFClusterEMEnergyCorrector.h +++ b/RecoParticleFlow/PFClusterProducer/interface/PFClusterEMEnergyCorrector.h @@ -5,14 +5,23 @@ #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" #include "DataFormats/ParticleFlowReco/interface/PFCluster.h" #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h" -#include "DataFormats/VertexReco/interface/VertexFwd.h" #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h" #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" -#include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h" #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/EcalDigi/interface/EBSrFlag.h" +#include "DataFormats/EcalDigi/interface/EESrFlag.h" +#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h" + +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h" +#include "CalibCalorimetry/EcalTPGTools/interface/EcalReadoutTools.h" +#include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h" + +#include "CondFormats/DataRecord/interface/GBRDWrapperRcd.h" +#include "CondFormats/EgammaObjects/interface/GBRForestD.h" class PFClusterEMEnergyCorrector { public: @@ -23,25 +32,54 @@ class PFClusterEMEnergyCorrector { void correctEnergies(const edm::Event &evt, const edm::EventSetup &es, const reco::PFCluster::EEtoPSAssociation &assoc, reco::PFClusterCollection& cs); private: - bool _applyCrackCorrections; - bool _applyMVACorrections; - double _maxPtForMVAEvaluation; - - bool autoDetectBunchSpacing_; - int bunchSpacingManual_; - + double maxPtForMVAEvaluation_; + + edm::EDGetTokenT ebSrFlagToken_; + edm::EDGetTokenT eeSrFlagToken_; + + //required for reading SR flags + edm::EDGetTokenT recHitsEB_; + edm::EDGetTokenT recHitsEE_; edm::EDGetTokenT bunchSpacing_; - edm::EDGetTokenT _recHitsEB; - edm::EDGetTokenT _recHitsEE; - - std::vector _condnames_mean_50ns; - std::vector _condnames_sigma_50ns; - std::vector _condnames_mean_25ns; - std::vector _condnames_sigma_25ns; + std::vector condnames_mean_; + std::vector condnames_sigma_; + + std::vector condnames_mean_25ns_; + std::vector condnames_sigma_25ns_; + std::vector condnames_mean_50ns_; + std::vector condnames_sigma_50ns_; + + bool srfAwareCorrection_; + bool applyCrackCorrections_; + bool applyMVACorrections_; + + bool autoDetectBunchSpacing_; + int bunchSpacingManual_; + + std::unique_ptr calibrator_; + void getAssociatedPSEnergy(const size_t clusIdx, const reco::PFCluster::EEtoPSAssociation &assoc, float& e1, float& e2); + + double meanlimlowEB_; + double meanlimhighEB_; + double meanoffsetEB_; + double meanscaleEB_; - std::unique_ptr _calibrator; + double meanlimlowEE_; + double meanlimhighEE_; + double meanoffsetEE_; + double meanscaleEE_; + double sigmalimlowEB_; + double sigmalimhighEB_; + double sigmaoffsetEB_; + double sigmascaleEB_; + + double sigmalimlowEE_; + double sigmalimhighEE_; + double sigmaoffsetEE_; + double sigmascaleEE_; + }; #endif diff --git a/RecoParticleFlow/PFClusterProducer/plugins/CorrectedECALPFClusterProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/CorrectedECALPFClusterProducer.cc index ac56278062bfd..7149e8c3f9092 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/CorrectedECALPFClusterProducer.cc +++ b/RecoParticleFlow/PFClusterProducer/plugins/CorrectedECALPFClusterProducer.cc @@ -131,13 +131,16 @@ void CorrectedECALPFClusterProducer::fillDescriptions(edm::ConfigurationDescript edm::ParameterSetDescription psd0; psd0.add("applyCrackCorrections",false); psd0.add("applyMVACorrections",false); + psd0.add("srfAwareCorrection",false); + psd0.add("autoDetectBunchSpacing",true); + psd0.add("bunchSpacing",25); psd0.add("maxPtForMVAEvaluation",-99.); psd0.add("algoName","PFClusterEMEnergyCorrector"); psd0.add("recHitsEBLabel",edm::InputTag("ecalRecHit","EcalRecHitsEB")); psd0.add("recHitsEELabel",edm::InputTag("ecalRecHit","EcalRecHitsEE")); psd0.add("verticesLabel",edm::InputTag("offlinePrimaryVertices")); - psd0.add("autoDetectBunchSpacing",true); - psd0.add("bunchSpacing",25); + psd0.add("ebSrFlagLabel",edm::InputTag("ecalDigis")); + psd0.add("eeSrFlagLabel",edm::InputTag("ecalDigis")); desc.add("energyCorrector",psd0); } desc.add("inputECAL",edm::InputTag("particleFlowClusterECALUncorrected")); diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterECAL_cff.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterECAL_cff.py new file mode 100644 index 0000000000000..5c876a52a36ad --- /dev/null +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterECAL_cff.py @@ -0,0 +1,10 @@ +import FWCore.ParameterSet.Config as cms + +from RecoParticleFlow.PFClusterProducer.particleFlowClusterECAL_cfi import particleFlowClusterECAL +particleFlowClusterECAL.energyCorrector.applyMVACorrections = True +particleFlowClusterECAL.energyCorrector.maxPtForMVAEvaluation = 90. + +from Configuration.Eras.Modifier_run2_ECAL_2017_cff import run2_ECAL_2017 +run2_ECAL_2017.toModify(particleFlowClusterECAL, + energyCorrector = dict(srfAwareCorrection = True, maxPtForMVAEvaluation = 300.)) + diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterECAL_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterECAL_cfi.py deleted file mode 100644 index 86bc88cab54c0..0000000000000 --- a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterECAL_cfi.py +++ /dev/null @@ -1,23 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -#### PF CLUSTER ECAL #### - -#energy corrector for corrected cluster producer -_emEnergyCorrector = cms.PSet( - algoName = cms.string("PFClusterEMEnergyCorrector"), - applyCrackCorrections = cms.bool(False), - applyMVACorrections = cms.bool(True), - maxPtForMVAEvaluation = cms.double(90.), - recHitsEBLabel = cms.InputTag('ecalRecHit', 'EcalRecHitsEB'), - recHitsEELabel = cms.InputTag('ecalRecHit', 'EcalRecHitsEE'), - autoDetectBunchSpacing = cms.bool(True), -) - -particleFlowClusterECAL = cms.EDProducer( - "CorrectedECALPFClusterProducer", - inputECAL = cms.InputTag("particleFlowClusterECALUncorrected"), - inputPS = cms.InputTag("particleFlowClusterPS"), - minimumPSEnergy = cms.double(0.0), - energyCorrector = _emEnergyCorrector - ) - diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterOOTECAL_cff.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterOOTECAL_cff.py index 216a04444c7d6..6d92581546b32 100644 --- a/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterOOTECAL_cff.py +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowClusterOOTECAL_cff.py @@ -1,5 +1,5 @@ import FWCore.ParameterSet.Config as cms -from RecoParticleFlow.PFClusterProducer.particleFlowClusterECAL_cfi import * +from RecoParticleFlow.PFClusterProducer.particleFlowClusterECAL_cff import * particleFlowClusterOOTECAL = particleFlowClusterECAL.clone() particleFlowClusterOOTECAL.inputECAL = cms.InputTag("particleFlowClusterOOTECALUncorrected") diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py index c241333fc685b..a2c3f2a6cf7a3 100644 --- a/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowCluster_cff.py @@ -11,8 +11,7 @@ from RecoParticleFlow.PFClusterProducer.particleFlowRecHitPS_cfi import * from RecoParticleFlow.PFClusterProducer.particleFlowClusterECALUncorrected_cfi import * -from RecoParticleFlow.PFClusterProducer.particleFlowClusterECAL_cfi import * - +from RecoParticleFlow.PFClusterProducer.particleFlowClusterECAL_cff import * from RecoParticleFlow.PFClusterProducer.particleFlowClusterHBHE_cfi import * from RecoParticleFlow.PFClusterProducer.particleFlowClusterHBHETimeSelected_cfi import * @@ -22,7 +21,6 @@ from RecoParticleFlow.PFClusterProducer.particleFlowClusterPS_cfi import * particleFlowClusterECALSequence = cms.Sequence(particleFlowClusterECAL) - pfClusteringECAL = cms.Sequence(particleFlowRecHitECAL* particleFlowClusterECALUncorrected * particleFlowClusterECALSequence) diff --git a/RecoParticleFlow/PFClusterProducer/src/PFClusterEMEnergyCorrector.cc b/RecoParticleFlow/PFClusterProducer/src/PFClusterEMEnergyCorrector.cc index cff62a09fe800..b87026e75a562 100644 --- a/RecoParticleFlow/PFClusterProducer/src/PFClusterEMEnergyCorrector.cc +++ b/RecoParticleFlow/PFClusterProducer/src/PFClusterEMEnergyCorrector.cc @@ -1,8 +1,5 @@ #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterEMEnergyCorrector.h" -#include "DataFormats/VertexReco/interface/Vertex.h" -#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h" -#include "CondFormats/DataRecord/interface/GBRDWrapperRcd.h" -#include "CondFormats/EgammaObjects/interface/GBRForestD.h" + #include "vdt/vdtMath.h" #include @@ -11,261 +8,424 @@ namespace { bool sortByKey(const EEPSPair& a, const EEPSPair& b) { return a.first < b.first; } + + double getOffset(const double lo, const double hi) { return lo + 0.5*(hi-lo); } + double getScale(const double lo, const double hi) { return 0.5*(hi-lo); } } PFClusterEMEnergyCorrector::PFClusterEMEnergyCorrector(const edm::ParameterSet& conf, edm::ConsumesCollector &&cc) : - _calibrator(new PFEnergyCalibration) { + calibrator_(new PFEnergyCalibration) { - _applyCrackCorrections = conf.getParameter("applyCrackCorrections"); - _applyMVACorrections = conf.getParameter("applyMVACorrections"); - _maxPtForMVAEvaluation = conf.getParameter("maxPtForMVAEvaluation"); - + applyCrackCorrections_ = conf.getParameter("applyCrackCorrections"); + applyMVACorrections_ = conf.getParameter("applyMVACorrections"); + srfAwareCorrection_ = conf.getParameter("srfAwareCorrection"); + + maxPtForMVAEvaluation_ = conf.getParameter("maxPtForMVAEvaluation"); + + if (applyMVACorrections_) { + + meanlimlowEB_ = -0.336; + meanlimhighEB_ = 0.916; + meanoffsetEB_ = getOffset(meanlimlowEB_, meanlimhighEB_); + meanscaleEB_ = getScale(meanlimlowEB_, meanlimhighEB_); + + meanlimlowEE_ = -0.336; + meanlimhighEE_ = 0.916; + meanoffsetEE_ = getOffset(meanlimlowEE_, meanlimhighEE_); + meanscaleEE_ = getScale(meanlimlowEE_, meanlimhighEE_); + + sigmalimlowEB_ = 0.001; + sigmalimhighEB_ = 0.4; + sigmaoffsetEB_ = getOffset(sigmalimlowEB_, sigmalimhighEB_); + sigmascaleEB_ = getScale(sigmalimlowEB_, sigmalimhighEB_); + + sigmalimlowEE_ = 0.001; + sigmalimhighEE_ = 0.4; + sigmaoffsetEE_ = getOffset(sigmalimlowEE_, sigmalimhighEE_); + sigmascaleEE_ = getScale(sigmalimlowEE_, sigmalimhighEE_); - if (_applyMVACorrections) { - _recHitsEB = cc.consumes(conf.getParameter("recHitsEBLabel")); - _recHitsEE = cc.consumes(conf.getParameter("recHitsEELabel")); + recHitsEB_ = cc.consumes(conf.getParameter("recHitsEBLabel")); + recHitsEE_ = cc.consumes(conf.getParameter("recHitsEELabel")); + autoDetectBunchSpacing_ = conf.getParameter("autoDetectBunchSpacing"); + + if (autoDetectBunchSpacing_) { + bunchSpacing_ = cc.consumes(edm::InputTag("bunchSpacingProducer")); + bunchSpacingManual_ = 0; + } else { + bunchSpacingManual_ = conf.getParameter("bunchSpacing"); + } + + condnames_mean_25ns_.insert(condnames_mean_25ns_.end(), { + "ecalPFClusterCorV2_EB_pfSize1_mean_25ns", + "ecalPFClusterCorV2_EB_pfSize2_mean_25ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin1_mean_25ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin2_mean_25ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin3_mean_25ns", + "ecalPFClusterCorV2_EE_pfSize1_mean_25ns", + "ecalPFClusterCorV2_EE_pfSize2_mean_25ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin1_mean_25ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin2_mean_25ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin3_mean_25ns"}); + condnames_sigma_25ns_.insert(condnames_sigma_25ns_.end(), { + "ecalPFClusterCorV2_EB_pfSize1_sigma_25ns", + "ecalPFClusterCorV2_EB_pfSize2_sigma_25ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin1_sigma_25ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin2_sigma_25ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin3_sigma_25ns", + "ecalPFClusterCorV2_EE_pfSize1_sigma_25ns", + "ecalPFClusterCorV2_EE_pfSize2_sigma_25ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin1_sigma_25ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin2_sigma_25ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin3_sigma_25ns"}); + condnames_mean_50ns_.insert(condnames_mean_50ns_.end(), { + "ecalPFClusterCorV2_EB_pfSize1_mean_50ns", + "ecalPFClusterCorV2_EB_pfSize2_mean_50ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin1_mean_50ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin2_mean_50ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin3_mean_50ns", + "ecalPFClusterCorV2_EE_pfSize1_mean_50ns", + "ecalPFClusterCorV2_EE_pfSize2_mean_50ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin1_mean_50ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin2_mean_50ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin3_mean_50ns"}); + condnames_sigma_50ns_.insert(condnames_sigma_50ns_.end(), { + "ecalPFClusterCorV2_EB_pfSize1_sigma_50ns", + "ecalPFClusterCorV2_EB_pfSize2_sigma_50ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin1_sigma_50ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin2_sigma_50ns", + "ecalPFClusterCorV2_EB_pfSize3_ptbin3_sigma_50ns", + "ecalPFClusterCorV2_EE_pfSize1_sigma_50ns", + "ecalPFClusterCorV2_EE_pfSize2_sigma_50ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin1_sigma_50ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin2_sigma_50ns", + "ecalPFClusterCorV2_EE_pfSize3_ptbin3_sigma_50ns"}); + + if (srfAwareCorrection_) { + + sigmalimlowEE_ = 0.001; + sigmalimhighEE_ = 0.1; + sigmaoffsetEE_ = getOffset(sigmalimlowEE_, sigmalimhighEE_); + sigmascaleEE_ = getScale(sigmalimlowEE_, sigmalimhighEE_); + + ebSrFlagToken_ = cc.consumes(conf.getParameter("ebSrFlagLabel")); + eeSrFlagToken_ = cc.consumes(conf.getParameter("eeSrFlagLabel")); + + condnames_mean_.insert(condnames_mean_.end(), { + "ecalPFClusterCor2017V2_EB_ZS_mean_25ns", + "ecalPFClusterCor2017V2_EB_Full_ptbin1_mean_25ns", + "ecalPFClusterCor2017V2_EB_Full_ptbin2_mean_25ns", + "ecalPFClusterCor2017V2_EB_Full_ptbin3_mean_25ns", + "ecalPFClusterCor2017V2_EE_ZS_mean_25ns", + "ecalPFClusterCor2017V2_EE_Full_ptbin1_mean_25ns", + "ecalPFClusterCor2017V2_EE_Full_ptbin2_mean_25ns", + "ecalPFClusterCor2017V2_EE_Full_ptbin3_mean_25ns"}); + + condnames_sigma_.insert(condnames_sigma_.end(), { + "ecalPFClusterCor2017V2_EB_ZS_sigma_25ns", + "ecalPFClusterCor2017V2_EB_Full_ptbin1_sigma_25ns", + "ecalPFClusterCor2017V2_EB_Full_ptbin2_sigma_25ns", + "ecalPFClusterCor2017V2_EB_Full_ptbin3_sigma_25ns", + "ecalPFClusterCor2017V2_EE_ZS_sigma_25ns", + "ecalPFClusterCor2017V2_EE_Full_ptbin1_sigma_25ns", + "ecalPFClusterCor2017V2_EE_Full_ptbin2_sigma_25ns", + "ecalPFClusterCor2017V2_EE_Full_ptbin3_sigma_25ns"}); + } + } - autoDetectBunchSpacing_ = conf.getParameter("autoDetectBunchSpacing"); +} + +void PFClusterEMEnergyCorrector::correctEnergies(const edm::Event &evt, + const edm::EventSetup &es, + const reco::PFCluster::EEtoPSAssociation &assoc, + reco::PFClusterCollection& cs) { + // First deal with pre-MVA corrections + // Kept for backward compatibility (and for HLT) + if (!applyMVACorrections_) { + for (unsigned int idx = 0; idxenergyEm(cluster,ePS1,ePS2,applyCrackCorrections_); + cluster.setCorrectedEnergy(correctedEnergy); + } + return; + } + + // Common objects for SRF-aware and old style corrections + EcalClusterLazyTools lazyTool(evt, es, recHitsEB_, recHitsEE_); + EcalReadoutTools readoutTool(evt, es); + + if (!srfAwareCorrection_) { + + int bunchspacing = 450; if (autoDetectBunchSpacing_) { - bunchSpacing_ = cc.consumes(edm::InputTag("bunchSpacingProducer")); - bunchSpacingManual_ = 0; + edm::Handle bunchSpacingH; + evt.getByToken(bunchSpacing_,bunchSpacingH); + bunchspacing = *bunchSpacingH; } else { - bunchSpacingManual_ = conf.getParameter("bunchSpacing"); + bunchspacing = bunchSpacingManual_; } - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EB_pfSize1_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EB_pfSize2_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin1_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin2_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin3_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EE_pfSize1_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EE_pfSize2_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin1_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin2_mean_50ns"); - _condnames_mean_50ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin3_mean_50ns"); + const std::vector& condnames_mean = (bunchspacing == 25) ? condnames_mean_25ns_ : condnames_mean_50ns_; + const std::vector& condnames_sigma = (bunchspacing == 25) ? condnames_sigma_25ns_ : condnames_sigma_50ns_; + const unsigned int ncor = condnames_mean.size(); + std::vector > forestH_mean(ncor); + std::vector > forestH_sigma(ncor); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EB_pfSize1_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EB_pfSize2_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin1_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin2_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin3_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EE_pfSize1_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EE_pfSize2_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin1_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin2_sigma_50ns"); - _condnames_sigma_50ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin3_sigma_50ns"); - - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EB_pfSize1_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EB_pfSize2_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin1_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin2_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin3_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EE_pfSize1_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EE_pfSize2_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin1_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin2_mean_25ns"); - _condnames_mean_25ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin3_mean_25ns"); + for (unsigned int icor=0; icor().get(condnames_mean[icor],forestH_mean[icor]); + es.get().get(condnames_sigma[icor],forestH_sigma[icor]); + } - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EB_pfSize1_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EB_pfSize2_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin1_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin2_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EB_pfSize3_ptbin3_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EE_pfSize1_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EE_pfSize2_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin1_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin2_sigma_25ns"); - _condnames_sigma_25ns.push_back("ecalPFClusterCorV2_EE_pfSize3_ptbin3_sigma_25ns"); - } - - - -} - - -void PFClusterEMEnergyCorrector::correctEnergies(const edm::Event &evt, const edm::EventSetup &es, const reco::PFCluster::EEtoPSAssociation &assoc, reco::PFClusterCollection& cs) { - - //legacy corrections - if (!_applyMVACorrections) { - for (unsigned int idx = 0; idx eval; + for (unsigned int idx = 0; idx0. && pt>maxPtForMVAEvaluation_) { + evale *= maxPtForMVAEvaluation_/pt; + } + double invE = (e == 0.) ? 0. : 1./e; //guard against dividing by 0. + int size = lazyTool.n5x5(cluster); + + int ietaix=0; + int iphiiy=0; + if (iseb) { + EBDetId ebseed(cluster.seed()); + ietaix = ebseed.ieta(); + iphiiy = ebseed.iphi(); + } else { + EEDetId eeseed(cluster.seed()); + ietaix = eeseed.ix(); + iphiiy = eeseed.iy(); + } + + //find index of corrections (0-4 for EB, 5-9 for EE, depending on cluster size and raw pt) + int coridx = std::min(size,3)-1; + if (coridx==2) { + if (pt>4.5) { + coridx += 1; + } + if (pt>18.) { + coridx += 1; + } + } + if (!iseb) { + coridx += 5; + } - //compute preshower energies for endcap clusters - double ePS1=0, ePS2=0; - if(!iseb) { - auto ee_key_val = std::make_pair(idx,edm::Ptr()); - const auto clustops = std::equal_range(assoc.begin(), - assoc.end(), - ee_key_val, - sortByKey); - for( auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) { - edm::Ptr psclus(i_ps->second); - switch( psclus->layer() ) { - case PFLayer::PS1: - ePS1 += psclus->energy(); - break; - case PFLayer::PS2: - ePS2 += psclus->energy(); - break; - default: - break; - } - } + const GBRForestD &meanforest = *forestH_mean[coridx].product(); + const GBRForestD &sigmaforest = *forestH_sigma[coridx].product(); + + //fill array for forest evaluation + eval[0] = evale; + eval[1] = ietaix; + eval[2] = iphiiy; + if (!iseb) { + eval[3] = ePS1*invE; + eval[4] = ePS2*invE; } - double correctedEnergy = _calibrator->energyEm(cluster,ePS1,ePS2,_applyCrackCorrections); - cluster.setCorrectedEnergy(correctedEnergy); + //these are the actual BDT responses + double rawmean = meanforest.GetResponse(eval.data()); + double rawsigma = sigmaforest.GetResponse(eval.data()); + + //apply transformation to limited output range (matching the training) + double mean = iseb ? meanoffsetEB_ + meanscaleEB_*vdt::fast_sin(rawmean) : meanoffsetEE_ + meanscaleEE_*vdt::fast_sin(rawmean); + double sigma = iseb ? sigmaoffsetEB_ + sigmascaleEB_*vdt::fast_sin(rawsigma) : sigmaoffsetEE_ + sigmascaleEE_*vdt::fast_sin(rawsigma); + + //regression target is ln(Etrue/Eraw) + //so corrected energy is ecor=exp(mean)*e, uncertainty is exp(mean)*eraw*sigma=ecor*sigma + double ecor = vdt::fast_exp(mean)*e; + double sigmacor = sigma*ecor; + + cluster.setCorrectedEnergy(ecor); + cluster.setCorrectedEnergyUncertainty(sigmacor); } return; } + + // Selective Readout Flags + edm::Handle ebSrFlags; + edm::Handle eeSrFlags; + evt.getByToken(ebSrFlagToken_, ebSrFlags); + evt.getByToken(eeSrFlagToken_, eeSrFlags); + if (!ebSrFlags.isValid() || !eeSrFlags.isValid()) + throw cms::Exception("PFClusterEMEnergyCorrector") << "This version of PFCluster corrections requires the SrFlagCollection information to proceed!\n"; - int bunchspacing = 450; - - if (autoDetectBunchSpacing_) { - edm::Handle bunchSpacingH; - evt.getByToken(bunchSpacing_,bunchSpacingH); - bunchspacing = *bunchSpacingH; - } - else { - bunchspacing = bunchSpacingManual_; - } - - const std::vector condnames_mean = (bunchspacing == 25) ? _condnames_mean_25ns : _condnames_mean_50ns; - const std::vector condnames_sigma = (bunchspacing == 25) ? _condnames_sigma_25ns : _condnames_sigma_50ns; - - const unsigned int ncor = condnames_mean.size(); + const unsigned int ncor = condnames_mean_.size(); std::vector > forestH_mean(ncor); std::vector > forestH_sigma(ncor); for (unsigned int icor=0; icor().get(condnames_mean[icor],forestH_mean[icor]); - es.get().get(condnames_sigma[icor],forestH_sigma[icor]); + es.get().get(condnames_mean_[icor],forestH_mean[icor]); + es.get().get(condnames_sigma_[icor],forestH_sigma[icor]); } - std::array eval; + std::array evalEB; + std::array evalEE; - EcalClusterLazyTools lazyTool(evt, es, _recHitsEB, _recHitsEE); - - //magic numbers for MINUIT-like transformation of BDT output onto limited range - //(These should be stored inside the conditions object in the future as well) - const double meanlimlow = -0.336; - const double meanlimhigh = 0.916; - const double meanoffset = meanlimlow + 0.5*(meanlimhigh-meanlimlow); - const double meanscale = 0.5*(meanlimhigh-meanlimlow); - - const double sigmalimlow = 0.001; - const double sigmalimhigh = 0.4; - const double sigmaoffset = sigmalimlow + 0.5*(sigmalimhigh-sigmalimlow); - const double sigmascale = 0.5*(sigmalimhigh-sigmalimlow); - for (unsigned int idx = 0; idx0. && pt>_maxPtForMVAEvaluation) { - evale *= _maxPtForMVAEvaluation/pt; - } - - double invE = (e == 0.) ? 0. : 1./e; //guard against dividing by 0. - - int size = lazyTool.n5x5(cluster); - - bool iseb = cluster.layer() == PFLayer::ECAL_BARREL; - - //find index of corrections (0-4 for EB, 5-9 for EE, depending on cluster size and raw pt) - int coridx = std::min(size,3)-1; - if (coridx==2) { - if (pt>4.5) { - coridx += 1; - } - if (pt>18.) { - coridx += 1; - } - } - if (!iseb) { - coridx += 5; + if (maxPtForMVAEvaluation_>0. && pt>maxPtForMVAEvaluation_) { + evale *= maxPtForMVAEvaluation_/pt; } + double invE = (e == 0.) ? 0. : 1./e; //guard against dividing by 0. + int size = lazyTool.n5x5(cluster); + int reducedHits = size; + if(size>=3) + reducedHits = 3; - const GBRForestD &meanforest = *forestH_mean[coridx].product(); - const GBRForestD &sigmaforest = *forestH_sigma[coridx].product(); - - //find seed crystal indices int ietaix=0; int iphiiy=0; if (iseb) { EBDetId ebseed(cluster.seed()); ietaix = ebseed.ieta(); iphiiy = ebseed.iphi(); - } - else { + } else { EEDetId eeseed(cluster.seed()); ietaix = eeseed.ix(); iphiiy = eeseed.iy(); } + // Hardcoded number are positions of modules boundaries of ECAL + int signeta = (ietaix > 0) ? 1 : -1; + int ietamod20 = (std::abs(ietaix) < 26) ? ietaix - signeta : (ietaix - 26*signeta) % 20; + int iphimod20 = (iphiiy-1) % 20; - //compute preshower energies for endcap clusters - double ePS1=0, ePS2=0; - if(!iseb) { - auto ee_key_val = std::make_pair(idx,edm::Ptr()); - const auto clustops = std::equal_range(assoc.begin(), - assoc.end(), - ee_key_val, - sortByKey); - for( auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) { - edm::Ptr psclus(i_ps->second); - switch( psclus->layer() ) { - case PFLayer::PS1: - ePS1 += psclus->energy(); - break; - case PFLayer::PS2: - ePS2 += psclus->energy(); - break; - default: - break; - } - } + int clusFlag = 0; + if (iseb){ + auto ecalUnit = readoutTool.readOutUnitOf(static_cast(cluster.seed())); + EBSrFlagCollection::const_iterator srf = ebSrFlags->find(ecalUnit); + if (srf != ebSrFlags->end()) + clusFlag = srf->value(); + else + clusFlag = 3; + } else { + auto ecalUnit = readoutTool.readOutUnitOf(static_cast(cluster.seed())); + EESrFlagCollection::const_iterator srf = eeSrFlags->find(ecalUnit); + if (srf != eeSrFlags->end()) + clusFlag = srf->value(); + else + clusFlag = 3; } - + + //find index of corrections (0-3 for EB, 4-7 for EE, depending on cluster size and raw pt) + int coridx = 0; + int regind = 0; + if (!iseb) regind = 4; + if (clusFlag==1) + coridx = 0 + regind; + else { + if (pt<2.5) coridx = 1 + regind; + else if (pt>=2.5 && pt<6.) coridx = 2 + regind; + else if (pt>=6.) coridx = 3 + regind; + } + if (clusFlag!=1 && clusFlag!=3) { + edm::LogWarning("PFClusterEMEnergyCorrector") << "We can only correct regions readout in ZS (flag 1) or FULL readout (flag 3). Flag " << clusFlag << " is not recognized." + << "\n" << "Assuming FULL readout and continuing"; + } + + const GBRForestD &meanforest = *forestH_mean[coridx].product(); + const GBRForestD &sigmaforest = *forestH_sigma[coridx].product(); + //fill array for forest evaluation - eval[0] = evale; - eval[1] = ietaix; - eval[2] = iphiiy; - if (!iseb) { - eval[3] = ePS1*invE; - eval[4] = ePS2*invE; + if (iseb) { + evalEB[0] = evale; + evalEB[1] = ietaix; + evalEB[2] = iphiiy; + evalEB[3] = ietamod20; + evalEB[4] = iphimod20; + evalEB[5] = reducedHits; + } else { + evalEE[0] = evale; + evalEE[1] = ietaix; + evalEE[2] = iphiiy; + evalEE[3] = (ePS1+ePS2)*invE; + evalEE[4] = reducedHits; } - + //these are the actual BDT responses - double rawmean = meanforest.GetResponse(eval.data()); - double rawsigma = sigmaforest.GetResponse(eval.data()); - + double rawmean = 1; + double rawsigma = 0; + + if(iseb){ + rawmean = meanforest.GetResponse(evalEB.data()); + rawsigma = sigmaforest.GetResponse(evalEB.data()); + } else { + rawmean = meanforest.GetResponse(evalEE.data()); + rawsigma = sigmaforest.GetResponse(evalEE.data()); + } + //apply transformation to limited output range (matching the training) - double mean = meanoffset + meanscale*vdt::fast_sin(rawmean); - double sigma = sigmaoffset + sigmascale*vdt::fast_sin(rawsigma); + //the training was done with different transformations for EB and EE (width only) + //makes a the code a bit more cumbersome, but it is not a problem per se + double mean = iseb ? meanoffsetEB_ + meanscaleEB_*vdt::fast_sin(rawmean) : meanoffsetEE_ + meanscaleEE_*vdt::fast_sin(rawmean); + double sigma = iseb ? sigmaoffsetEB_ + sigmascaleEB_*vdt::fast_sin(rawsigma) : sigmaoffsetEE_ + sigmascaleEE_*vdt::fast_sin(rawsigma); //regression target is ln(Etrue/Eraw) //so corrected energy is ecor=exp(mean)*e, uncertainty is exp(mean)*eraw*sigma=ecor*sigma - double ecor = exp(mean)*e; + double ecor = iseb ? vdt::fast_exp(mean)*e : vdt::fast_exp(mean)*(e+ePS1+ePS2); double sigmacor = sigma*ecor; + + LogDebug("PFClusterEMEnergyCorrector") << "ieta : iphi : ietamod20 : iphimod20 : size : reducedHits = " + << ietaix << " " << iphiiy << " " + << ietamod20 << " " << iphimod20 << " " + << size << " " << reducedHits + << "\n" << "isEB : eraw : ePS1 : ePS2 : (eps1+eps2)/raw : Flag = " + << iseb << " " << evale << " " << ePS1 << " " << ePS2 << " " << (ePS1+ePS2)/evale << " " << clusFlag + << "\n" << "response : correction = " + << exp(mean) << " " << ecor; cluster.setCorrectedEnergy(ecor); cluster.setCorrectedEnergyUncertainty(sigmacor); - } } +void PFClusterEMEnergyCorrector::getAssociatedPSEnergy(const size_t clusIdx, const reco::PFCluster::EEtoPSAssociation &assoc, float& e1, float& e2) { + e1 = 0; + e2 = 0; + auto ee_key_val = std::make_pair(clusIdx,edm::Ptr()); + const auto clustops = std::equal_range(assoc.begin(), + assoc.end(), + ee_key_val, + sortByKey); + for( auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) { + edm::Ptr psclus(i_ps->second); + switch( psclus->layer() ) { + case PFLayer::PS1: + e1 += psclus->energy(); + break; + case PFLayer::PS2: + e2 += psclus->energy(); + break; + default: + break; + } + } +} +