diff --git a/CommonTools/ParticleFlow/plugins/BuildFile.xml b/CommonTools/ParticleFlow/plugins/BuildFile.xml index 6cf896a0564dc..edb79453fdcbc 100644 --- a/CommonTools/ParticleFlow/plugins/BuildFile.xml +++ b/CommonTools/ParticleFlow/plugins/BuildFile.xml @@ -14,4 +14,6 @@ + + diff --git a/CommonTools/ParticleFlow/plugins/PFCandidateRecalibrator.cc b/CommonTools/ParticleFlow/plugins/PFCandidateRecalibrator.cc new file mode 100644 index 0000000000000..f0cab8ac8c1f7 --- /dev/null +++ b/CommonTools/ParticleFlow/plugins/PFCandidateRecalibrator.cc @@ -0,0 +1,354 @@ +/// Take: +// - a PF candidate collection (which uses bugged HCAL respcorrs) +// - respCorrs values from fixed GT and bugged GT +// Produce: +// - a new PFCandidate collection containing the recalibrated PFCandidates in HF and where the neutral had pointing to problematic cells are removed +// - a second PFCandidate collection with just those discarded hadrons +// - a ValueMap that maps the old to the new, and vice-versa + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/ESWatcher.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "FWCore/Framework/interface/Event.h" + +#include "CalibFormats/HcalObjects/interface/HcalDbService.h" +#include "CalibFormats/HcalObjects/interface/HcalDbRecord.h" +#include "CalibFormats/HcalObjects/interface/HcalCalibrations.h" +#include "Geometry/CaloTopology/interface/HcalTopology.h" +#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" + +#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h" +#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" + +#include "DataFormats/Common/interface/View.h" +#include "DataFormats/Common/interface/Association.h" +#include + + +struct HEChannel { + float eta; float phi; float ratio; + HEChannel(float eta, float phi, float ratio): + eta(eta), + phi(phi), + ratio(ratio) + {} +}; +struct HFChannel { + int ieta; int iphi; int depth; float ratio; + HFChannel(int ieta, int iphi, int depth, float ratio): + ieta(ieta), + iphi(iphi), + depth(depth), + ratio(ratio) + {} +}; + + +class PFCandidateRecalibrator : public edm::stream::EDProducer<> { + public: + PFCandidateRecalibrator(const edm::ParameterSet&); + ~PFCandidateRecalibrator() override {}; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + private: + void beginRun(const edm::Run& iRun, edm::EventSetup const& iSetup) override; + void endRun(const edm::Run& iRun, edm::EventSetup const& iSetup) override; + void produce(edm::Event&, const edm::EventSetup&) override; + + edm::ESWatcher hcalDbWatcher_; + edm::ESWatcher hcalRCWatcher_; + + edm::EDGetTokenT pfcandidates_; + + std::vector badChHE_; + std::vector badChHF_; + + float shortFibreThr_; + float longFibreThr_; +}; + +PFCandidateRecalibrator::PFCandidateRecalibrator(const edm::ParameterSet &iConfig) : + pfcandidates_(consumes(iConfig.getParameter("pfcandidates"))), + shortFibreThr_(iConfig.getParameter("shortFibreThr")), + longFibreThr_(iConfig.getParameter("longFibreThr")) +{ + produces(); + produces("discarded"); + produces>(); +} + +void PFCandidateRecalibrator::beginRun(const edm::Run &iRun, const edm::EventSetup &iSetup) +{ + if (hcalDbWatcher_.check(iSetup) || hcalRCWatcher_.check(iSetup)) + { + //Get Calib Constants from current GT + edm::ESHandle gtCond; + iSetup.get().get(gtCond); + + //Get Calib Constants from bugged tag + edm::ESHandle htopo; + iSetup.get().get(htopo); + const HcalTopology* theHBHETopology = htopo.product(); + + edm::ESHandle buggedCond; + iSetup.get().get("bugged", buggedCond); + HcalRespCorrs buggedRespCorrs(*buggedCond.product()); + buggedRespCorrs.setTopo(theHBHETopology); + + //access calogeometry + edm::ESHandle calogeom; + iSetup.get().get(calogeom); + const CaloGeometry* cgeo = calogeom.product(); + const HcalGeometry* hgeom = static_cast(cgeo->getSubdetectorGeometry(DetId::Hcal,HcalForward)); + + //reset the bad channel containers + badChHE_.clear(); + badChHF_.clear(); + + //fill bad cells HE (use eta, phi) + const std::vector& cellsHE = hgeom->getValidDetIds(DetId::Detector::Hcal, HcalEndcap); + for( auto id : cellsHE) + { + float currentRespCorr = gtCond->getHcalRespCorr(id)->getValue(); + float buggedRespCorr = buggedRespCorrs.getValues(id)->getValue(); + if (buggedRespCorr == 0.) + continue; + + float ratio = currentRespCorr/buggedRespCorr; + if( std::abs(ratio - 1.f) > 0.001 ) + { + GlobalPoint pos = hgeom->getPosition(id); + badChHE_.push_back(HEChannel(pos.eta(),pos.phi(),ratio)); + } + } + + //fill bad cells HF (use ieta, iphi) + auto const& cellsHF = hgeom->getValidDetIds(DetId::Detector::Hcal, HcalForward); + for( auto id : cellsHF) + { + float currentRespCorr = gtCond->getHcalRespCorr(id)->getValue(); + float buggedRespCorr = buggedRespCorrs.getValues(id)->getValue(); + if (buggedRespCorr == 0.) + continue; + + float ratio = currentRespCorr/buggedRespCorr; + if( std::abs(ratio - 1.f) > 0.001 ) + { + HcalDetId dummyId(id); + badChHF_.push_back(HFChannel(dummyId.ieta(), dummyId.iphi(), dummyId.depth(), ratio)); + } + } + } +} + +void PFCandidateRecalibrator::endRun(const edm::Run &iRun, const edm::EventSetup &iSetup) +{ +} + +void PFCandidateRecalibrator::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) +{ + //access calogeometry + edm::ESHandle calogeom; + iSetup.get().get(calogeom); + const CaloGeometry* cgeo = calogeom.product(); + const HcalGeometry* hgeom = static_cast(cgeo->getSubdetectorGeometry(DetId::Hcal,HcalForward)); + + //access PFCandidates + edm::Handle pfcandidates; + iEvent.getByToken(pfcandidates_, pfcandidates); + + int nPfCand = pfcandidates->size(); + std::unique_ptr copy(new reco::PFCandidateCollection()); + std::unique_ptr discarded(new reco::PFCandidateCollection()); + copy->reserve(nPfCand); + std::vector oldToNew(nPfCand), newToOld, badToOld; + newToOld.reserve(nPfCand); + + LogDebug("PFCandidateRecalibrator") << "NEW EV:"; + + //loop over PFCandidates + int i = -1; + for (const reco::PFCandidate &pf : *pfcandidates) + { + ++i; + float absEta = std::abs(pf.eta()); + + //deal with HE + if( pf.particleId() == reco::PFCandidate::ParticleType::h0 && + !badChHE_.empty() && //don't touch if no miscalibration is found + absEta > 1.4 && absEta < 3.) + { + bool toKill = false; + for(auto const& badIt: badChHE_) + if( reco::deltaR2(pf.eta(), pf.phi(), badIt.eta, badIt.phi) < 0.07 ) + toKill = true; + + + if(toKill) + { + discarded->push_back(pf); + oldToNew[i] = (-discarded->size()); + badToOld.push_back(i); + continue; + } + else + { + copy->push_back(pf); + oldToNew[i] = (copy->size()); + newToOld.push_back(i); + } + } + //deal with HF + else if( (pf.particleId() == reco::PFCandidate::ParticleType::h_HF || pf.particleId() == reco::PFCandidate::ParticleType::egamma_HF) && + !badChHF_.empty() && //don't touch if no miscalibration is found + absEta >= 3.) + { + const math::XYZPointF& ecalPoint = pf.positionAtECALEntrance(); + GlobalPoint ecalGPoint(ecalPoint.X(),ecalPoint.Y(),ecalPoint.Z()); + HcalDetId closestDetId(hgeom->getClosestCell(ecalGPoint)); + + if(closestDetId.subdet() == HcalForward) + { + HcalDetId hDetId(closestDetId.subdet(),closestDetId.ieta(),closestDetId.iphi(),1); //depth1 + + //raw*calEnergy() is the same as *calEnergy() - no corrections are done for HF + float longE = pf.rawEcalEnergy() + pf.rawHcalEnergy()/2.; //depth1 + float shortE = pf.rawHcalEnergy()/2.; //depth2 + + float ecalEnergy = pf.rawEcalEnergy(); + float hcalEnergy = pf.rawHcalEnergy(); + float totEnergy = ecalEnergy + hcalEnergy; + float totEnergyOrig = totEnergy; + + bool toKill = false; + + for(auto const& badIt: badChHF_) + { + if ( hDetId.ieta() == badIt.ieta && + hDetId.iphi() == badIt.iphi ) + { + LogDebug("PFCandidateRecalibrator") << "==> orig en (tot,H,E): " + << pf.energy() << " " << pf.rawHcalEnergy() << " " << pf.rawEcalEnergy(); + if(badIt.depth == 1) //depth1 + { + longE *= badIt.ratio; + ecalEnergy = ((longE - shortE) > 0.) ? (longE - shortE) : 0.; + totEnergy = ecalEnergy + hcalEnergy; + } + else //depth2 + { + shortE *= badIt.ratio; + hcalEnergy = 2*shortE; + ecalEnergy = ((longE - shortE) > 0.) ? (longE - shortE) : 0.; + totEnergy = ecalEnergy + hcalEnergy; + } + //kill candidate if goes below thr + if((pf.particleId() == reco::PFCandidate::ParticleType::h_HF && shortE < shortFibreThr_) || + (pf.particleId() == reco::PFCandidate::ParticleType::egamma_HF && longE < longFibreThr_)) + toKill = true; + + LogDebug("PFCandidateRecalibrator") << "====> ieta,iphi,depth: " + << badIt.ieta << " " << badIt.iphi << " " << badIt.depth << " corr: " << badIt.ratio; + LogDebug("PFCandidateRecalibrator") << "====> recal en (tot,H,E): " + << totEnergy << " " << hcalEnergy << " " << ecalEnergy; + + } + } + + if(toKill) + { + discarded->push_back(pf); + oldToNew[i] = (-discarded->size()); + badToOld.push_back(i); + + LogDebug("PFCandidateRecalibrator") << "==> KILLED "; + } + else + { + copy->push_back(pf); + oldToNew[i] = (copy->size()); + newToOld.push_back(i); + + copy->back().setHcalEnergy(hcalEnergy, hcalEnergy); + copy->back().setEcalEnergy(ecalEnergy, ecalEnergy); + + float scalingFactor = totEnergy/totEnergyOrig; + math::XYZTLorentzVector recalibP4 = pf.p4() * scalingFactor; + copy->back().setP4( recalibP4 ); + + LogDebug("PFCandidateRecalibrator") << "====> stored en (tot,H,E): " + << copy->back().energy() << " " << copy->back().hcalEnergy() << " " << copy->back().ecalEnergy(); + } + } + else + { + copy->push_back(pf); + oldToNew[i] = (copy->size()); + newToOld.push_back(i); + } + } + else + { + copy->push_back(pf); + oldToNew[i] = (copy->size()); + newToOld.push_back(i); + } + } + + // Now we put things in the event + edm::OrphanHandle newpf = iEvent.put(std::move(copy)); + edm::OrphanHandle badpf = iEvent.put(std::move(discarded), "discarded"); + + std::unique_ptr> pf2pf(new edm::ValueMap()); + edm::ValueMap::Filler filler(*pf2pf); + std::vector refs; refs.reserve(nPfCand); + + // old to new + for (auto iOldToNew : oldToNew){ + if (iOldToNew > 0) { + refs.push_back(reco::PFCandidateRef(newpf, iOldToNew-1)); + } else { + refs.push_back(reco::PFCandidateRef(badpf,-iOldToNew-1)); + } + } + filler.insert(pfcandidates, refs.begin(), refs.end()); + // new good to old + refs.clear(); + for (int i : newToOld) { + refs.push_back(reco::PFCandidateRef(pfcandidates,i)); + } + filler.insert(newpf, refs.begin(), refs.end()); + // new bad to old + refs.clear(); + for (int i : badToOld) { + refs.push_back(reco::PFCandidateRef(pfcandidates,i)); + } + filler.insert(badpf, refs.begin(), refs.end()); + // done + filler.fill(); + iEvent.put(std::move(pf2pf)); +} + +void +PFCandidateRecalibrator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) +{ + edm::ParameterSetDescription desc; + + desc.add("pfcandidates", edm::InputTag("particleFlow")); + desc.add("shortFibreThr", 1.4); + desc.add("longFibreThr", 1.4); + + descriptions.add("pfCandidateRecalibrator", desc); + +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(PFCandidateRecalibrator); diff --git a/Configuration/AlCa/python/autoCond.py b/Configuration/AlCa/python/autoCond.py index dd6291b6808dd..30f037da6ae1f 100644 --- a/Configuration/AlCa/python/autoCond.py +++ b/Configuration/AlCa/python/autoCond.py @@ -10,25 +10,25 @@ # GlobalTag for MC production (p-Pb collisions) with realistic alignment and calibrations for Run1 'run1_mc_pa' : '94X_mcRun1_pA_v1', # GlobalTag for MC production with perfectly aligned and calibrated detector for Run2 - 'run2_design' : '94X_mcRun2_design_v2', + 'run2_design' : '94X_mcRun2_design_v4', # GlobalTag for MC production with pessimistic alignment and calibrations for Run2 - 'run2_mc_50ns' : '94X_mcRun2_startup_v2', + 'run2_mc_50ns' : '94X_mcRun2_startup_v3', #GlobalTag for MC production with optimistic alignment and calibrations for Run2 - 'run2_mc' : '94X_mcRun2_asymptotic_v2', + 'run2_mc' : '94X_mcRun2_asymptotic_v3', # GlobalTag for MC production (L1 Trigger Stage1) with starup-like alignment and calibrations for Run2, L1 trigger in Stage1 mode 'run2_mc_l1stage1' : '93X_mcRun2_asymptotic_v1', # GlobalTag for MC production (cosmics) with starup-like alignment and calibrations for Run2, Strip tracker in peak mode - 'run2_mc_cosmics' : '94X_mcRun2cosmics_startup_deco_v2', + 'run2_mc_cosmics' : '94X_mcRun2cosmics_startup_deco_v3', # GlobalTag for MC production (Heavy Ions collisions) with optimistic alignment and calibrations for Run2 'run2_mc_hi' : '93X_mcRun2_HeavyIon_v0', # GlobalTag for MC production (p-Pb collisions) with realistic alignment and calibrations for Run2 - 'run2_mc_pa' : '94X_mcRun2_pA_v2', + 'run2_mc_pa' : '94X_mcRun2_pA_v3', # GlobalTag for Run1 data reprocessing - 'run1_data' : '94X_dataRun2_v6', + 'run1_data' : '94X_dataRun2_v10', # GlobalTag for Run2 data reprocessing - 'run2_data' : '94X_dataRun2_v6', + 'run2_data' : '94X_dataRun2_v10', # GlobalTag for Run2 data relvals: allows customization to run with fixed L1 menu - 'run2_data_relval' : '94X_dataRun2_relval_v11', + 'run2_data_relval' : '94X_dataRun2_relval_v14', # GlobalTag for Run2 data 2016H relvals only: Prompt Conditions + fixed L1 menu (to be removed) 'run2_data_promptlike' : '94X_dataRun2_PromptLike_v9', # GlobalTag for Run1 HLT: it points to the online GT @@ -40,13 +40,13 @@ # GlobalTag for Run2 HLT for HI: it points to the online GT 'run2_hlt_hi' : '94X_dataRun2_HLTHI_frozen_v5', # GlobalTag for MC production with perfectly aligned and calibrated detector for Phase1 2017 (and 0,0,~0-centred beamspot) - 'phase1_2017_design' : '94X_mc2017_design_IdealBS_v10', + 'phase1_2017_design' : '94X_mc2017_design_IdealBS_v11', # GlobalTag for MC production with realistic conditions for Phase1 2017 detector - 'phase1_2017_realistic' : '94X_mc2017_realistic_v14', + 'phase1_2017_realistic' : '94X_mc2017_realistic_v15', # GlobalTag for MC production (cosmics) with realistic alignment and calibrations for Phase1 2017 detector, Strip tracker in DECO mode - 'phase1_2017_cosmics' : '94X_mc2017cosmics_realistic_deco_v10', + 'phase1_2017_cosmics' : '94X_mc2017cosmics_realistic_deco_v11', # GlobalTag for MC production (cosmics) with realistic alignment and calibrations for Phase1 2017 detector, Strip tracker in PEAK mode - 'phase1_2017_cosmics_peak' : '94X_mc2017cosmics_realistic_peak_v10', + 'phase1_2017_cosmics_peak' : '94X_mc2017cosmics_realistic_peak_v11', # GlobalTag for MC production with perfectly aligned and calibrated detector for full Phase1 2018 (and 0,0,0-centred beamspot) 'phase1_2018_design' : '94X_upgrade2018_design_IdealBS_v7', # GlobalTag for MC production with realistic conditions for full Phase1 2018 detector diff --git a/Configuration/PyReleaseValidation/python/relval_steps.py b/Configuration/PyReleaseValidation/python/relval_steps.py index 42e78d065981c..1b9720e1cf1c0 100644 --- a/Configuration/PyReleaseValidation/python/relval_steps.py +++ b/Configuration/PyReleaseValidation/python/relval_steps.py @@ -1960,12 +1960,15 @@ def gen2018HiMix(fragment,howMuch): '--data' : '', '--scenario' : 'pp', '--eventcontent' : 'MINIAOD,DQM', - '--datatier' : 'MINIAOD,DQMIO' + '--datatier' : 'MINIAOD,DQMIO', + '--customise_unsch' : 'PhysicsTools/PatAlgos/slimming/customizeMiniAOD_HcalFixLegacy2016.customizeAll' },stepMiniAODDefaults]) steps['REMINIAOD_data2016_HIPM'] = merge([{'--era' : 'Run2_2016_HIPM,run2_miniAOD_80XLegacy'},steps['REMINIAOD_data2016']]) -steps['REMINIAOD_data2017'] = merge([{'--era' : 'Run2_2017,run2_miniAOD_94XFall17'},steps['REMINIAOD_data2016']]) +stepReMiniAODData17 = merge([{'--era' : 'Run2_2017,run2_miniAOD_94XFall17'},steps['REMINIAOD_data2016']]) +stepReMiniAODData17 = remove(stepReMiniAODData17,'--customise_unsch') +steps['REMINIAOD_data2017'] = stepReMiniAODData17 # Not sure whether the customisations are in the dict as "--customise" or "--era" so try to # remove both. Currently premixing uses "--customise" and everything else uses "--era". diff --git a/PhysicsTools/PatAlgos/python/slimming/customizeMiniAOD_HcalFixLegacy2016.py b/PhysicsTools/PatAlgos/python/slimming/customizeMiniAOD_HcalFixLegacy2016.py new file mode 100644 index 0000000000000..c9964f8763dfd --- /dev/null +++ b/PhysicsTools/PatAlgos/python/slimming/customizeMiniAOD_HcalFixLegacy2016.py @@ -0,0 +1,105 @@ +import FWCore.ParameterSet.Config as cms + +from PhysicsTools.PatAlgos.tools.helpers import MassSearchReplaceAnyInputTagVisitor, addKeepStatement +from PhysicsTools.PatAlgos.tools.helpers import getPatAlgosToolsTask + +def loadJetMETBTag(process): + + task = getPatAlgosToolsTask(process) + + import RecoParticleFlow.PFProducer.pfLinker_cff + process.particleFlowPtrs = RecoParticleFlow.PFProducer.pfLinker_cff.particleFlowPtrs.clone() + task.add(process.particleFlowPtrs) + + process.load("CommonTools.ParticleFlow.pfNoPileUpIso_cff") + task.add(process.pfNoPileUpIsoTask) + process.load("CommonTools.ParticleFlow.ParticleSelectors.pfSortByType_cff") + task.add(process.pfSortByTypeTask) + + import RecoJets.Configuration.RecoPFJets_cff + process.ak4PFJetsCHS = RecoJets.Configuration.RecoPFJets_cff.ak4PFJetsCHS.clone() + task.add(process.ak4PFJetsCHS) + # need also the non-CHS ones as they are used to seed taus + process.ak4PFJets = RecoJets.Configuration.RecoPFJets_cff.ak4PFJets.clone() + task.add(process.ak4PFJets) + process.ak8PFJetsCHS = RecoJets.Configuration.RecoPFJets_cff.ak8PFJetsCHS.clone() + task.add(process.ak8PFJetsCHS) + + process.fixedGridRhoAll = RecoJets.Configuration.RecoPFJets_cff.fixedGridRhoAll.clone() + process.fixedGridRhoFastjetAll = RecoJets.Configuration.RecoPFJets_cff.fixedGridRhoFastjetAll.clone() + process.fixedGridRhoFastjetCentral = RecoJets.Configuration.RecoPFJets_cff.fixedGridRhoFastjetCentral.clone() + process.fixedGridRhoFastjetCentralChargedPileUp = RecoJets.Configuration.RecoPFJets_cff.fixedGridRhoFastjetCentralChargedPileUp.clone() + process.fixedGridRhoFastjetCentralNeutral = RecoJets.Configuration.RecoPFJets_cff.fixedGridRhoFastjetCentralNeutral.clone() + task.add( process.fixedGridRhoAll, + process.fixedGridRhoFastjetAll, + process.fixedGridRhoFastjetCentral, + process.fixedGridRhoFastjetCentralChargedPileUp, + process.fixedGridRhoFastjetCentralNeutral ) + + process.load("RecoJets.JetAssociationProducers.ak4JTA_cff") + task.add(process.ak4JetTracksAssociatorAtVertexPF) + + process.load('RecoBTag.Configuration.RecoBTag_cff') + task.add(process.btaggingTask) + + process.load("RecoMET.METProducers.PFMET_cfi") + task.add(process.pfMet) + + +def cleanPfCandidates(process, verbose=False): + task = getPatAlgosToolsTask(process) + + #add producer at the beginning of the schedule + process.load("CommonTools.ParticleFlow.pfCandidateRecalibrator_cfi") + task.add(process.pfCandidateRecalibrator) + + replacePFCandidates = MassSearchReplaceAnyInputTagVisitor("particleFlow", "pfCandidateRecalibrator", verbose=verbose) + replacePFTmpPtrs = MassSearchReplaceAnyInputTagVisitor("particleFlowTmpPtrs", "particleFlowPtrs", verbose=verbose) + for everywhere in [ process.producers, process.filters, process.analyzers, process.psets, process.vpsets ]: + for name,obj in everywhere.iteritems(): + if obj != process.pfCandidateRecalibrator: + replacePFCandidates.doIt(obj, name) + replacePFTmpPtrs.doIt(obj, name) + + + process.load("CommonTools.ParticleFlow.pfEGammaToCandidateRemapper_cfi") + task.add(process.pfEGammaToCandidateRemapper) + process.pfEGammaToCandidateRemapper.pf2pf = cms.InputTag("pfCandidateRecalibrator") + process.reducedEgamma.gsfElectronsPFValMap = cms.InputTag("pfEGammaToCandidateRemapper","electrons") + process.reducedEgamma.photonsPFValMap = cms.InputTag("pfEGammaToCandidateRemapper","photons") + + +def addDiscardedPFCandidates(process, inputCollection, verbose=False): + + task = getPatAlgosToolsTask(process) + + process.primaryVertexAssociationDiscardedCandidates = process.primaryVertexAssociation.clone( + particles = inputCollection, + ) + task.add(process.primaryVertexAssociationDiscardedCandidates) + + process.packedPFCandidatesDiscarded = process.packedPFCandidates.clone( + inputCollection = inputCollection, + PuppiNoLepSrc = cms.InputTag(""), + PuppiSrc = cms.InputTag(""), + secondaryVerticesForWhiteList = cms.VInputTag(), + vertexAssociator = cms.InputTag("primaryVertexAssociationDiscardedCandidates","original") + ) + task.add(process.packedPFCandidatesDiscarded) + + addKeepStatement(process, "keep patPackedCandidates_packedPFCandidates_*_*", + ["keep patPackedCandidates_packedPFCandidatesDiscarded_*_*"], + verbose=verbose) + + +def customizeAll(process, verbose=False): + + if verbose: + print "===>>> customizing the process for legacy rereco 2016" + + loadJetMETBTag(process) + + cleanPfCandidates(process, verbose) + addDiscardedPFCandidates(process, cms.InputTag("pfCandidateRecalibrator","discarded"), verbose=verbose) + + return process diff --git a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py index 30a68a9f4b642..c854dcbc466c6 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py @@ -97,3 +97,6 @@ from Configuration.Eras.Modifier_run2_miniAOD_94XFall17_cff import run2_miniAOD_94XFall17 modifyReducedEGammaRun2MiniAOD9XFall17_ = run2_miniAOD_94XFall17.makeProcessModifier(calibrateReducedEgamma) +from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy +modifyReducedEGammaRun2MiniAOD8XLegacy_ = run2_miniAOD_80XLegacy.makeProcessModifier(calibrateReducedEgamma) + diff --git a/RecoEgamma/EgammaTools/interface/ElectronEnergyCalibrator.h b/RecoEgamma/EgammaTools/interface/ElectronEnergyCalibrator.h index 350c32c4ddb03..89d9342f2a70d 100644 --- a/RecoEgamma/EgammaTools/interface/ElectronEnergyCalibrator.h +++ b/RecoEgamma/EgammaTools/interface/ElectronEnergyCalibrator.h @@ -39,7 +39,10 @@ class ElectronEnergyCalibrator //set the minimum et to apply the correction to void setMinEt(float val){minEt_=val;} - + //sets whether to use the smeared ecal energy in the combination + //note, if this is true and the E/p combination is not trained using this smeared value, this is a bug + //the E/p combination must get the ecalEnergyErr used in its training + void setUseSmearCorrEcalEnergyErrInComb(bool val){useSmearCorrEcalEnergyErrInComb_=val;} /// Correct this electron. /// StreamID is needed when used with CMSSW Random Number Generator std::array @@ -72,6 +75,7 @@ class ElectronEnergyCalibrator const EpCombinationTool *epCombinationTool_; //this is not owned TRandom *rng_; //this is not owned float minEt_; + bool useSmearCorrEcalEnergyErrInComb_; //default values to access if no correction available static const EnergyScaleCorrection::ScaleCorrection defaultScaleCorr_; diff --git a/RecoEgamma/EgammaTools/interface/EpCombinationTool.h b/RecoEgamma/EgammaTools/interface/EpCombinationTool.h index d4e6100676280..60a0bf9728612 100644 --- a/RecoEgamma/EgammaTools/interface/EpCombinationTool.h +++ b/RecoEgamma/EgammaTools/interface/EpCombinationTool.h @@ -26,6 +26,7 @@ class EpCombinationTool void setEventContent(const edm::EventSetup& iSetup); std::pair combine(const reco::GsfElectron& electron) const; + std::pair combine(const reco::GsfElectron& electron,float corrEcalEnergyErr) const; private: EgammaRegressionContainer ecalTrkEnergyRegress_; diff --git a/RecoEgamma/EgammaTools/plugins/CalibratedElectronProducers.cc b/RecoEgamma/EgammaTools/plugins/CalibratedElectronProducers.cc index 55333eb685493..6333edac5169a 100644 --- a/RecoEgamma/EgammaTools/plugins/CalibratedElectronProducers.cc +++ b/RecoEgamma/EgammaTools/plugins/CalibratedElectronProducers.cc @@ -112,6 +112,7 @@ CalibratedElectronProducerT::CalibratedElectronProducerT( const edm::Paramete produceCalibratedObjs_(conf.getParameter("produceCalibratedObjs")) { energyCorrector_.setMinEt(conf.getParameter("minEtToCalibrate")); + energyCorrector_.setUseSmearCorrEcalEnergyErrInComb(conf.getParameter("useSmearCorrEcalEnergyErrInComb")); if (conf.getParameter("semiDeterministic")) { semiDeterministicRng_.reset(new TRandom2()); @@ -137,6 +138,7 @@ void CalibratedElectronProducerT::fillDescriptions(edm::ConfigurationDescript desc.add("minEtToCalibrate",5.0); desc.add("produceCalibratedObjs",true); desc.add("semiDeterministic",true); + desc.add("useSmearCorrEcalEnergyErrInComb",false); std::vector valMapsProduced; for(auto varToStore : valMapsToStore_) valMapsProduced.push_back(EGEnergySysIndex::name(varToStore)); desc.add >("valueMapsStored",valMapsProduced)->setComment("provides to python configs the list of valuemaps stored, can not be overriden in the python config"); diff --git a/RecoEgamma/EgammaTools/plugins/EG8XObjectUpdateModifier.cc b/RecoEgamma/EgammaTools/plugins/EG8XObjectUpdateModifier.cc new file mode 100644 index 0000000000000..df65c2b13f176 --- /dev/null +++ b/RecoEgamma/EgammaTools/plugins/EG8XObjectUpdateModifier.cc @@ -0,0 +1,134 @@ +#include "CommonTools/CandAlgos/interface/ModifyObjectValueBase.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/Common/interface/Handle.h" +#include "Geometry/CaloTopology/interface/CaloTopology.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHit.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" + +#include + +//this modifier fills variables where not present in CMSSW_8X +//use case is when reading older 80X samples in newer releases, aka legacy + +class EG8XObjectUpdateModifier : public ModifyObjectValueBase { +public: + EG8XObjectUpdateModifier(const edm::ParameterSet& conf); + ~EG8XObjectUpdateModifier() override{} + + void setEvent(const edm::Event&) final; + void setEventContent(const edm::EventSetup&) final; + void setConsumes(edm::ConsumesCollector&) final; + + void modifyObject(reco::GsfElectron& ele) const final; + void modifyObject(reco::Photon& pho) const final; + + void modifyObject(pat::Electron& ele) const final{return modifyObject(static_cast(ele));} + void modifyObject(pat::Photon& pho) const final{return modifyObject(static_cast(pho));} + +private: + std::pair getSaturationInfo(const reco::SuperCluster& superClus)const; + + edm::ESHandle caloTopoHandle_; + edm::Handle ecalRecHitsEBHandle_; + edm::Handle ecalRecHitsEEHandle_; + edm::EDGetTokenT ecalRecHitsEBToken_; + edm::EDGetTokenT ecalRecHitsEEToken_; + edm::InputTag ecalRecHitsEBTag_; + edm::InputTag ecalRecHitsEETag_; + +}; + +EG8XObjectUpdateModifier::EG8XObjectUpdateModifier(const edm::ParameterSet& conf): + ModifyObjectValueBase(conf), + ecalRecHitsEBTag_(conf.getParameter("ecalRecHitsEB")), + ecalRecHitsEETag_(conf.getParameter("ecalRecHitsEE")) +{ + + +} + +void EG8XObjectUpdateModifier::setConsumes(edm::ConsumesCollector& cc) +{ + ecalRecHitsEBToken_ = cc.consumes(ecalRecHitsEBTag_); + ecalRecHitsEEToken_ = cc.consumes(ecalRecHitsEETag_); +} + +void EG8XObjectUpdateModifier::setEvent(const edm::Event& iEvent) +{ + iEvent.getByToken(ecalRecHitsEBToken_,ecalRecHitsEBHandle_); + iEvent.getByToken(ecalRecHitsEEToken_,ecalRecHitsEEHandle_); +} + +void EG8XObjectUpdateModifier::setEventContent(const edm::EventSetup& iSetup) +{ + iSetup.get().get(caloTopoHandle_); +} + + +void EG8XObjectUpdateModifier::modifyObject(reco::GsfElectron& ele)const +{ + + const reco::CaloCluster& seedClus = *(ele.superCluster()->seed()); + const EcalRecHitCollection* ecalRecHits = ele.isEB() ? &*ecalRecHitsEBHandle_ : &*ecalRecHitsEEHandle_; + const auto* caloTopo = caloTopoHandle_.product(); + + auto full5x5ShowerShapes = ele.full5x5_showerShape(); + full5x5ShowerShapes.e2x5Left = noZS::EcalClusterTools::e2x5Left (seedClus,ecalRecHits,caloTopo); + full5x5ShowerShapes.e2x5Right = noZS::EcalClusterTools::e2x5Right (seedClus,ecalRecHits,caloTopo); + full5x5ShowerShapes.e2x5Top = noZS::EcalClusterTools::e2x5Top (seedClus,ecalRecHits,caloTopo); + full5x5ShowerShapes.e2x5Bottom = noZS::EcalClusterTools::e2x5Bottom(seedClus,ecalRecHits,caloTopo); + ele.full5x5_setShowerShape(full5x5ShowerShapes); + + auto showerShapes = ele.showerShape(); + showerShapes.e2x5Left = EcalClusterTools::e2x5Left (seedClus,ecalRecHits,caloTopo); + showerShapes.e2x5Right = EcalClusterTools::e2x5Right (seedClus,ecalRecHits,caloTopo); + showerShapes.e2x5Top = EcalClusterTools::e2x5Top (seedClus,ecalRecHits,caloTopo); + showerShapes.e2x5Bottom = EcalClusterTools::e2x5Bottom(seedClus,ecalRecHits,caloTopo); + ele.setShowerShape(showerShapes); + + reco::GsfElectron::SaturationInfo eleSatInfo; + auto satInfo = getSaturationInfo(*ele.superCluster()); + eleSatInfo.nSaturatedXtals = satInfo.first; + eleSatInfo.isSeedSaturated = satInfo.second; + ele.setSaturationInfo(eleSatInfo); +} + +void EG8XObjectUpdateModifier::modifyObject(reco::Photon& pho)const +{ + reco::Photon::SaturationInfo phoSatInfo; + auto satInfo = getSaturationInfo(*pho.superCluster()); + phoSatInfo.nSaturatedXtals = satInfo.first; + phoSatInfo.isSeedSaturated = satInfo.second; + pho.setSaturationInfo(phoSatInfo); +} + + +std::pair EG8XObjectUpdateModifier::getSaturationInfo(const reco::SuperCluster& superClus)const +{ + bool isEB = superClus.seed()->seed().subdetId()==EcalBarrel; + const auto& ecalRecHits = isEB ? *ecalRecHitsEBHandle_ : *ecalRecHitsEEHandle_; + + int nrSatCrys = 0; + bool seedSaturated = false; + const auto& hitsAndFractions = superClus.seed()->hitsAndFractions(); + for(const auto& hitFractionPair : hitsAndFractions) { + auto ecalRecHitIt = ecalRecHits.find(hitFractionPair.first); + if(ecalRecHitIt != ecalRecHits.end() && + ecalRecHitIt->checkFlag(EcalRecHit::Flags::kSaturated)){ + nrSatCrys++; + if(hitFractionPair.first == superClus.seed()->seed()) seedSaturated = true; + } + } + return {nrSatCrys,seedSaturated}; +} + +DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, + EG8XObjectUpdateModifier, + "EG8XObjectUpdateModifier"); diff --git a/RecoEgamma/EgammaTools/plugins/EGEtScaleSysModifier.cc b/RecoEgamma/EgammaTools/plugins/EGEtScaleSysModifier.cc new file mode 100644 index 0000000000000..12ca932f37027 --- /dev/null +++ b/RecoEgamma/EgammaTools/plugins/EGEtScaleSysModifier.cc @@ -0,0 +1,216 @@ +#include "CommonTools/CandAlgos/interface/ModifyObjectValueBase.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "RecoEgamma/EgammaTools/interface/EpCombinationTool.h" +#include "RecoEgamma/EgammaTools/interface/EGEnergySysIndex.h" + +#include + +//in the legacy re-reco we came across an Et scale issue, specifically an inflection at 45 GeV +//it is problematic to address in the normal scale and smearing code therefore +//this patch modifies the e/gamma object to "patch in" this additional systematic +//Questions: +//1 ) Why dont we add this to the overall scale systematic? +//Well we could but the whole point of this is to get a proper template which can be evolved +//correctly. The issue is not that the current systematic doesnt cover this systematic on a +//per bin basis, it does, the problem is the sign flips at 45 GeV so its fine if you use +//only eles/phos > 45 or <45 GeV but the template cant simultaneously cover the entire spectrum +//and you'll get a nasty constrained fit. +//And if you care about these sort of things (and you probably should), you shouldnt be using +//the overall systematic anyways but its seperate parts +// +//2 ) Why dont you do it inside the standard class +//We could but we're rather hoping to solve this issue more cleanly in the future +//but we would hold up the reminiAOD too long to do so so this is just to get something +//which should be fine for 90% of analyses in the miniAOD. So this is a temporary fix +//which we hope to rip out soon (if you're reading this in 2024 because we're still doing it +//this way, all I can say is sorry, we thought it would go away soon! ) + + +class EGEtScaleSysModifier : public ModifyObjectValueBase { +public: + EGEtScaleSysModifier(const edm::ParameterSet& conf); + ~EGEtScaleSysModifier() override{} + + void setEvent(const edm::Event&) final; + void setEventContent(const edm::EventSetup&) final; + void setConsumes(edm::ConsumesCollector&) final; + + void modifyObject(pat::Electron& ele) const final; + void modifyObject(pat::Photon& pho) const final; + +private: + std::pair calCombinedMom(reco::GsfElectron& ele,const float scale, + const float smear)const; + void setEcalEnergy(reco::GsfElectron& ele,const float scale,const float smear)const; + + class UncertFuncBase { + public: + UncertFuncBase(){} + virtual ~UncertFuncBase(){} + virtual float val(const float et)const=0; + }; + + //defines two uncertaintes, one EtY + //for X("lowEt")),highEt_(conf.getParameter("highEt")), + lowEtUncert_(conf.getParameter("lowEtUncert")), + highEtUncert_(conf.getParameter("highEtUncert")), + dEt_(highEt_-lowEt_), + dUncert_(highEtUncert_ - lowEtUncert_) + { + if(highEt_<=lowEt_) throw cms::Exception("ConfigError") <<" highEt "<=highEt_) return highEtUncert_; + else{ + return (et-lowEt_)*dUncert_/dEt_+lowEtUncert_; + } + } + private: + float lowEt_; + float highEt_; + float lowEtUncert_; + float highEtUncert_; + float dEt_; + float dUncert_; + }; + + EpCombinationTool epCombTool_; + std::unique_ptr uncertFunc_; +}; + +EGEtScaleSysModifier::EGEtScaleSysModifier(const edm::ParameterSet& conf): + ModifyObjectValueBase(conf), + epCombTool_(conf.getParameter("epCombConfig")) +{ + const edm::ParameterSet funcPSet = conf.getParameter("uncertFunc"); + const std::string& funcName = funcPSet.getParameter("name"); + if(funcName == "UncertFuncV1"){ + uncertFunc_ = std::make_unique(funcPSet); + }else{ + throw cms::Exception("ConfigError") <<"Error constructing EGEtScaleSysModifier, function name "<val(ele.et()/ele.energy()*ecalEnergyPostCorr); + const float smear = getVal(ele,EGEnergySysIndex::kSmearValue); + const float corr = getVal(ele,EGEnergySysIndex::kScaleValue); + + //get the values we have to reset back to + const float oldEcalEnergy = ele.ecalEnergy(); + const float oldEcalEnergyErr = ele.ecalEnergyError(); + const auto oldP4 = ele.p4(); + const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION); + const float oldTrkMomErr = ele.trackMomentumError(); + + ele.setCorrectedEcalEnergy(ecalEnergyPreCorr); + ele.setCorrectedEcalEnergyError(ecalEnergyErrPreCorr); + + const float energyEtUncertUp = calCombinedMom(ele,corr+etUncert,smear).first; + const float energyEtUncertDn = calCombinedMom(ele,corr-etUncert,smear).first; + + //reset it back to how it was + ele.setCorrectedEcalEnergy(oldEcalEnergy); + ele.setCorrectedEcalEnergyError(oldEcalEnergyErr); + ele.correctMomentum(oldP4,oldTrkMomErr,oldP4Err); + + ele.addUserFloat("energyScaleEtUp",energyEtUncertUp); + ele.addUserFloat("energyScaleEtDown",energyEtUncertDn); + +} + +void EGEtScaleSysModifier::modifyObject(pat::Photon& pho)const +{ + auto getVal = [](const pat::Photon& pho,EGEnergySysIndex::Index valIndex){ + return pho.userFloat(EGEnergySysIndex::name(valIndex)); + }; + //so pho.energy() may be either pre or post corrections, we have no idea + //so we explicity access the pre and post correction ecal energies + //post corr for the et value for the systematic, pre corr to apply them + const float ecalEnergyPostCorr = getVal(pho,EGEnergySysIndex::kEcalPostCorr); + const float ecalEnergyPreCorr = getVal(pho,EGEnergySysIndex::kEcalPreCorr); + + //the et cut is in terms of post corr ecal energy + const float etUncert = uncertFunc_->val(pho.et()/pho.energy()*ecalEnergyPostCorr); + const float corr = getVal(pho,EGEnergySysIndex::kScaleValue); + + const float energyEtUncertUp = ecalEnergyPreCorr*(corr+etUncert); + const float energyEtUncertDn = ecalEnergyPreCorr*(corr-etUncert); + + pho.addUserFloat("energyScaleEtUp",energyEtUncertUp); + pho.addUserFloat("energyScaleEtDown",energyEtUncertDn); + +} + +std::pair EGEtScaleSysModifier::calCombinedMom(reco::GsfElectron& ele, + const float scale, + const float smear)const +{ + const float oldEcalEnergy = ele.ecalEnergy(); + const float oldEcalEnergyErr = ele.ecalEnergyError(); + const auto oldP4 = ele.p4(); + const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION); + const float oldTrkMomErr = ele.trackMomentumError(); + + setEcalEnergy(ele,scale,smear); + const auto& combinedMomentum = epCombTool_.combine(ele); + ele.setCorrectedEcalEnergy(oldEcalEnergy); + ele.setCorrectedEcalEnergyError(oldEcalEnergyErr); + ele.correctMomentum(oldP4,oldTrkMomErr,oldP4Err); + + + return combinedMomentum; +} + +void EGEtScaleSysModifier::setEcalEnergy(reco::GsfElectron& ele, + const float scale, + const float smear)const +{ + const float oldEcalEnergy = ele.ecalEnergy(); + const float oldEcalEnergyErr = ele.ecalEnergyError(); + ele.setCorrectedEcalEnergy( oldEcalEnergy*scale ); + ele.setCorrectedEcalEnergyError(std::hypot( oldEcalEnergyErr*scale, oldEcalEnergy*smear*scale ) ); +} + + + +DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, + EGEtScaleSysModifier, + "EGEtScaleSysModifier"); diff --git a/RecoEgamma/EgammaTools/python/calibratedEgammas_cff.py b/RecoEgamma/EgammaTools/python/calibratedEgammas_cff.py index af497a18c2c3c..a9458e3fb998c 100644 --- a/RecoEgamma/EgammaTools/python/calibratedEgammas_cff.py +++ b/RecoEgamma/EgammaTools/python/calibratedEgammas_cff.py @@ -1,12 +1,18 @@ import FWCore.ParameterSet.Config as cms +_correctionFile2016Legacy = "EgammaAnalysis/ElectronTools/data/ScalesSmearings/Legacy2016_07Aug2017_FineEtaR9_v3_ele_unc" +_correctionFile2017Nov17 = "EgammaAnalysis/ElectronTools/data/ScalesSmearings/Run2017_17Nov2017_v1_ele_unc" + calibratedEgammaSettings = cms.PSet(minEtToCalibrate = cms.double(5.0), semiDeterministic = cms.bool(True), - correctionFile = cms.string("EgammaAnalysis/ElectronTools/data/ScalesSmearings/Run2017_17Nov2017_v1_ele_unc"), + correctionFile = cms.string(_correctionFile2017Nov17), recHitCollectionEB = cms.InputTag('reducedEcalRecHitsEB'), recHitCollectionEE = cms.InputTag('reducedEcalRecHitsEE'), produceCalibratedObjs = cms.bool(True) ) +from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy +run2_miniAOD_80XLegacy.toModify(calibratedEgammaSettings,correctionFile = _correctionFile2016Legacy) + calibratedEgammaPatSettings = calibratedEgammaSettings.clone( recHitCollectionEB = cms.InputTag('reducedEgamma','reducedEBRecHits'), recHitCollectionEE = cms.InputTag('reducedEgamma','reducedEERecHits') @@ -40,17 +46,28 @@ ) +#a bug was found, the smear corrected ecal energy err was used in the E/p combination +#as the campaign has already started, we can not change the default nor the behaviour for the 94X Fall17 reminiAOD +#therefore we set useSmearCorrEcalEnergyErrInComb = cms.bool(True) when it should be set to false +#make no mistake, this is an incorrect configuration which will cause a scale shift at Et = 50 GeV calibratedElectrons = cms.EDProducer("CalibratedElectronProducer", - calibratedEgammaSettings, + calibratedEgammaSettings, + useSmearCorrEcalEnergyErrInComb = cms.bool(True), epCombConfig = ecalTrkCombinationRegression, src = cms.InputTag('gedGsfElectrons'), ) calibratedPatElectrons = cms.EDProducer("CalibratedPatElectronProducer", calibratedEgammaPatSettings, + useSmearCorrEcalEnergyErrInComb = cms.bool(True), epCombConfig = ecalTrkCombinationRegression, src = cms.InputTag('slimmedElectrons'), ) +#the 80XLegacy miniAOD had not started, hence we can fix it for that +run2_miniAOD_80XLegacy.toModify(calibratedElectrons,useSmearCorrEcalEnergyErrInComb = False) +run2_miniAOD_80XLegacy.toModify(calibratedPatElectrons,useSmearCorrEcalEnergyErrInComb = False) + + calibratedPhotons = cms.EDProducer("CalibratedPhotonProducer", calibratedEgammaSettings, @@ -63,3 +80,6 @@ def prefixName(prefix,name): return prefix+name[0].upper()+name[1:] + + + diff --git a/RecoEgamma/EgammaTools/python/egammaObjectModificationsInMiniAOD_cff.py b/RecoEgamma/EgammaTools/python/egammaObjectModificationsInMiniAOD_cff.py index f7c9638f4cac9..41a4d6ff3e22f 100644 --- a/RecoEgamma/EgammaTools/python/egammaObjectModificationsInMiniAOD_cff.py +++ b/RecoEgamma/EgammaTools/python/egammaObjectModificationsInMiniAOD_cff.py @@ -129,10 +129,46 @@ def setup_mva(val_pset,cat_pset,prod_name,mva_name): for valueMapName in RecoEgamma.EgammaTools.calibratedPhotonProducerTRecoPhoton_cfi.calibratedPhotonProducerTRecoPhoton.valueMapsStored: setattr(reducedEgammaEnergyScaleAndSmearingModifier.photon_config,valueMapName,cms.InputTag("reducedEgamma",prefixName("calibPho",valueMapName))) +############################################################# +# 8X to 9X modifiers (fills in variables new to 9X w.r.t 8X) +############################################################# +egamma8XObjectUpdateModifier = cms.PSet( + modifierName = cms.string('EG8XObjectUpdateModifier'), + ecalRecHitsEB = cms.InputTag("reducedEgamma","reducedEBRecHits"), + ecalRecHitsEE = cms.InputTag("reducedEgamma","reducedEERecHits"), +) + +############################################################# +# 8X legacy needs an extra Et scale systematic +# due to an inflection around 45 GeV which is handled as a +# patch on top of the standard scale and smearing systematics +############################################################# +from RecoEgamma.EgammaTools.calibratedEgammas_cff import ecalTrkCombinationRegression +egamma8XLegacyEtScaleSysModifier = cms.PSet( + modifierName = cms.string('EGEtScaleSysModifier'), + epCombConfig = ecalTrkCombinationRegression, + uncertFunc = cms.PSet( + name = cms.string("UncertFuncV1"), + lowEt = cms.double(43.5), + highEt = cms.double(46.5), + lowEtUncert = cms.double(0.002), + highEtUncert = cms.double(-0.002) + ) + ) def appendReducedEgammaEnergyScaleAndSmearingModifier(modifiers): modifiers.append(reducedEgammaEnergyScaleAndSmearingModifier) +def prependEgamma8XObjectUpdateModifier(modifiers): + modifiers.insert(0,egamma8XObjectUpdateModifier) + +def appendEgamma8XLegacyAppendableModifiers (modifiers): + modifiers.append(reducedEgammaEnergyScaleAndSmearingModifier) + modifiers.append(egamma8XLegacyEtScaleSysModifier) + from Configuration.Eras.Modifier_run2_miniAOD_94XFall17_cff import run2_miniAOD_94XFall17 run2_miniAOD_94XFall17.toModify(egamma_modifications,appendReducedEgammaEnergyScaleAndSmearingModifier) +from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy +run2_miniAOD_80XLegacy.toModify(egamma_modifications,appendEgamma8XLegacyAppendableModifiers) +run2_miniAOD_80XLegacy.toModify(egamma_modifications,prependEgamma8XObjectUpdateModifier) diff --git a/RecoEgamma/EgammaTools/src/ElectronEnergyCalibrator.cc b/RecoEgamma/EgammaTools/src/ElectronEnergyCalibrator.cc index c8d092f425634..5eeac8abea891 100644 --- a/RecoEgamma/EgammaTools/src/ElectronEnergyCalibrator.cc +++ b/RecoEgamma/EgammaTools/src/ElectronEnergyCalibrator.cc @@ -13,7 +13,8 @@ ElectronEnergyCalibrator::ElectronEnergyCalibrator(const EpCombinationTool &comb correctionRetriever_(correctionFile), epCombinationTool_(&combinator), rng_(nullptr), - minEt_(1.0) + minEt_(1.0), + useSmearCorrEcalEnergyErrInComb_(false) { } @@ -148,8 +149,8 @@ setEnergyAndSystVarations(const float scale,const float smearNrSigma,const float energyData[EGEnergySysIndex::kSmearUp] = calCombinedMom(ele,corrUp,smearUp).first; energyData[EGEnergySysIndex::kSmearDown] = calCombinedMom(ele,corrDn,smearDn).first; + const std::pair combinedMomentum = calCombinedMom(ele,corr,smear); setEcalEnergy(ele,corr,smear); - const std::pair combinedMomentum = epCombinationTool_->combine(ele); const float energyCorr = combinedMomentum.first / oldP4.t(); const math::XYZTLorentzVector newP4(oldP4.x() * energyCorr, @@ -183,10 +184,19 @@ std::pair ElectronEnergyCalibrator::calCombinedMom(reco::GsfElectro { const float oldEcalEnergy = ele.ecalEnergy(); const float oldEcalEnergyErr = ele.ecalEnergyError(); + + const auto oldP4 = ele.p4(); + const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION); + const float oldTrkMomErr = ele.trackMomentumError(); + setEcalEnergy(ele,scale,smear); - const auto& combinedMomentum = epCombinationTool_->combine(ele); + float ecalEnergyErrForComb = useSmearCorrEcalEnergyErrInComb_ ? ele.correctedEcalEnergyError() : oldEcalEnergyErr*scale; + const auto& combinedMomentum = epCombinationTool_->combine(ele,ecalEnergyErrForComb); + ele.setCorrectedEcalEnergy(oldEcalEnergy); ele.setCorrectedEcalEnergyError(oldEcalEnergyErr); + ele.correctMomentum(oldP4,oldTrkMomErr,oldP4Err); + return combinedMomentum; } diff --git a/RecoEgamma/EgammaTools/src/EnergyScaleCorrection.cc b/RecoEgamma/EgammaTools/src/EnergyScaleCorrection.cc index 0c5f594afd55a..e884a0014a552 100644 --- a/RecoEgamma/EgammaTools/src/EnergyScaleCorrection.cc +++ b/RecoEgamma/EgammaTools/src/EnergyScaleCorrection.cc @@ -8,6 +8,7 @@ #include #include + EnergyScaleCorrection::EnergyScaleCorrection(const std::string& correctionFileName, unsigned int genSeed): smearingType_(ECALELF) { @@ -359,7 +360,21 @@ EnergyScaleCorrection::CorrectionCategory::CorrectionCategory(const std::string& r9Min_ = -1; r9Max_ = 0.94; }; - + // R9 region + p1 = category.find("-R9"); + p2 = p1 + 1; + if(p1 != std::string::npos) { + p1 = category.find("_", p1); + p2 = category.find("_", p1 + 1); + r9Min_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1)); + // If there is one value, just set lower bound + if (p2 != std::string::npos) { + p1 = p2; + p2 = category.find("-", p1); + r9Max_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1)); + if(r9Max_>=1.0) r9Max_ = std::numeric_limits::max(); + } + } //------------------------------ p1 = category.find("gainEle_"); // Position of first character if(p1 != std::string::npos) { diff --git a/RecoEgamma/EgammaTools/src/EpCombinationTool.cc b/RecoEgamma/EgammaTools/src/EpCombinationTool.cc index b347e3dd8a76b..814dd6a9341c3 100644 --- a/RecoEgamma/EgammaTools/src/EpCombinationTool.cc +++ b/RecoEgamma/EgammaTools/src/EpCombinationTool.cc @@ -43,13 +43,20 @@ void EpCombinationTool::setEventContent(const edm::EventSetup& iSetup) } std::pair EpCombinationTool::combine(const reco::GsfElectron& ele)const +{ + return combine(ele,ele.correctedEcalEnergyError()); +} + +//when doing the E/p combination, its very important to ensure the ecalEnergyErr +//that the regression is trained on is used, not the actual ecalEnergyErr of the electron +//these differ when you correct the ecalEnergyErr by smearing value needed to get data/MC to agree +std::pair EpCombinationTool::combine(const reco::GsfElectron& ele,const float corrEcalEnergyErr)const { const float scRawEnergy = ele.superCluster()->rawEnergy(); const float esEnergy = ele.superCluster()->preshowerEnergy(); const float corrEcalEnergy = ele.correctedEcalEnergy(); - const float corrEcalEnergyErr = ele.correctedEcalEnergyError(); const float ecalMean = ele.correctedEcalEnergy() / (scRawEnergy+esEnergy); const float ecalSigma = corrEcalEnergyErr / corrEcalEnergy;