diff --git a/Configuration/ProcessModifiers/python/photonDRN_cff.py b/Configuration/ProcessModifiers/python/photonDRN_cff.py new file mode 100644 index 0000000000000..438677ec800c3 --- /dev/null +++ b/Configuration/ProcessModifiers/python/photonDRN_cff.py @@ -0,0 +1,9 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.ProcessModifiers.enableSonicTriton_cff import enableSonicTriton + +#behind-the-scenes modifier that only turns on the DRN photon regression +_photonDRN = cms.Modifier() + +#modifier to enable DRN energy regression for photons +#requires also having enableSonicTriton +photonDRN = cms.ModifierChain(_photonDRN, enableSonicTriton) diff --git a/Configuration/PyReleaseValidation/README.md b/Configuration/PyReleaseValidation/README.md index 2d597155457e6..f6b3aa1a361aa 100644 --- a/Configuration/PyReleaseValidation/README.md +++ b/Configuration/PyReleaseValidation/README.md @@ -53,6 +53,7 @@ The offsets currently in use are: * 0.17: Run-3 deep core seeding for JetCore iteration * 0.21: Production-like sequence * 0.24: 0 Tesla (Run-2, Run-3) +* 0.31: Photon energy corrections with DRN architecture * 0.61: `phase2_ecal_devel` era * 0.91: Track DNN modifier * 0.97: Premixing stage1 diff --git a/Configuration/PyReleaseValidation/python/relval_2017.py b/Configuration/PyReleaseValidation/python/relval_2017.py index 1686cbef71e4e..0a6d5192ad573 100644 --- a/Configuration/PyReleaseValidation/python/relval_2017.py +++ b/Configuration/PyReleaseValidation/python/relval_2017.py @@ -18,6 +18,7 @@ # (TTbar trackingOnly, trackingRun2, trackingOnlyRun2, trackingLowPU, pixelTrackingOnly) # (TTbar PU with JME NanoAOD, disable for now due to Run-3 Nano-Prompt preparation) # 2018 (ele guns 10, 35, 1000; pho guns 10, 35; mu guns 1, 10, 100, 1000, QCD 3TeV, QCD Flat) +# (pho guns 10, 35 with photonDRN enabled) # 2018 (ZMM, TTbar, ZEE, MinBias, TTbar PU, ZEE PU, TTbar design) # (TTbar trackingOnly, pixelTrackingOnly) # (HE collapse: TTbar, TTbar PU, TTbar design) @@ -47,6 +48,7 @@ 10024.1,10024.2,10024.3,10024.4,10024.5, #10224.15, 10801.0,10802.0,10803.0,10804.0,10805.0,10806.0,10807.0,10808.0,10809.0,10859.0,10871.0, + 10804.31, 10805.31, 10842.0,10824.0,10825.0,10826.0,10823.0,11024.0,11025.0,11224.0, 10824.1,10824.5, 10824.6,11024.6,11224.6, diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index f079efe18c78a..37567986655d5 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -490,6 +490,32 @@ def condition(self, fragment, stepList, key, hasHarvest): '--procModifiers': 'mlpf' } +# photonDRN workflows +class UpgradeWorkflow_photonDRN(UpgradeWorkflow): + def setup_(self, step, stepName, stepDict, k, properties): + if 'Reco' in step: + stepDict[stepName][k] = merge([self.step3, stepDict[step][k]]) + def condition(self, fragment, stepList, key, hasHarvest): + return '2018' in key and "SingleGamma" in fragment + +upgradeWFs['photonDRN'] = UpgradeWorkflow_photonDRN( + steps = [ + 'Reco', + 'RecoNano', + 'RecoFakeHLT' + ], + PU = [ + 'Reco', + 'RecoNano', + 'RecoFakeHLT' + ], + suffix = '_photonDRN', + offset = 0.31, +) +upgradeWFs['photonDRN'].step3 = { + '--procModifiers': 'enableSonicTriton,photonDRN' +} + # Patatrack workflows: # - 2018 conditions, TTbar # - 2018 conditions, Z->mumu, diff --git a/PhysicsTools/PatAlgos/python/slimming/patElectronDRNCorrector_cfi.py b/PhysicsTools/PatAlgos/python/slimming/patElectronDRNCorrector_cfi.py new file mode 100644 index 0000000000000..94fe00cb1b93a --- /dev/null +++ b/PhysicsTools/PatAlgos/python/slimming/patElectronDRNCorrector_cfi.py @@ -0,0 +1,15 @@ +import FWCore.ParameterSet.Config as cms + +from RecoEgamma.EgammaTools.patElectronDRNCorrectionProducer_cfi import patElectronDRNCorrectionProducer + +patElectronsDRN = patElectronDRNCorrectionProducer.clone( + particleSource = 'selectedPatElectrons', + rhoName = 'fixedGridRhoFastjetAll', + Client = patElectronDRNCorrectionProducer.Client.clone( + mode = 'Async', + allowedTries = 1, + modelName = 'electronObjectEnsemble', + modelConfigPath = 'RecoEgamma/EgammaElectronProducers/data/models/electronObjectEnsemble/config.pbtxt', + timeout = 10 + ) + ) diff --git a/PhysicsTools/PatAlgos/python/slimming/patPhotonDRNCorrector_cfi.py b/PhysicsTools/PatAlgos/python/slimming/patPhotonDRNCorrector_cfi.py new file mode 100644 index 0000000000000..68ee77d41f1e1 --- /dev/null +++ b/PhysicsTools/PatAlgos/python/slimming/patPhotonDRNCorrector_cfi.py @@ -0,0 +1,14 @@ +import FWCore.ParameterSet.Config as cms +from RecoEgamma.EgammaTools.patPhotonDRNCorrectionProducer_cfi import patPhotonDRNCorrectionProducer + +patPhotonsDRN = patPhotonDRNCorrectionProducer.clone( + particleSource = 'selectedPatPhotons', + rhoName = 'fixedGridRhoFastjetAll', + Client = patPhotonDRNCorrectionProducer.Client.clone( + mode = 'Async', + allowedTries = 1, + modelName = 'photonObjectEnsemble', + modelConfigPath = 'RecoEgamma/EgammaPhotonProducers/data/models/photonObjectEnsemble/config.pbtxt', + timeout = 10 + ) + ) diff --git a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py index 306efa71e690e..91edf0e01dc43 100644 --- a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py +++ b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py @@ -102,3 +102,7 @@ _phase2_timing_slimmingTask = cms.Task(slimmingTask.copy(), offlineSlimmedPrimaryVertices4D) phase2_timing.toReplaceWith(slimmingTask,_phase2_timing_slimmingTask) + +from PhysicsTools.PatAlgos.slimming.patPhotonDRNCorrector_cfi import patPhotonsDRN +from Configuration.ProcessModifiers.photonDRN_cff import _photonDRN +_photonDRN.toReplaceWith(slimmingTask, cms.Task(slimmingTask.copy(), patPhotonsDRN)) diff --git a/RecoEgamma/EgammaTools/BuildFile.xml b/RecoEgamma/EgammaTools/BuildFile.xml index 33d90ff89a012..8a5d985ab6994 100644 --- a/RecoEgamma/EgammaTools/BuildFile.xml +++ b/RecoEgamma/EgammaTools/BuildFile.xml @@ -17,6 +17,7 @@ + diff --git a/RecoEgamma/EgammaTools/plugins/DRNCorrectionProducerT.cc b/RecoEgamma/EgammaTools/plugins/DRNCorrectionProducerT.cc new file mode 100644 index 0000000000000..9efcb4bb5e7c7 --- /dev/null +++ b/RecoEgamma/EgammaTools/plugins/DRNCorrectionProducerT.cc @@ -0,0 +1,426 @@ +#include "HeterogeneousCore/SonicTriton/interface/TritonEDProducer.h" +#include "HeterogeneousCore/SonicTriton/interface/TritonData.h" + +#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 "DataFormats/PatCandidates/interface/Photon.h" +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/EgammaCandidates/interface/Photon.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h" +#include "DataFormats/EgammaCandidates/interface/PhotonCore.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" +#include "DataFormats/EcalDetId/interface/ESDetId.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHit.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" + +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h" +#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h" +#include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h" +#include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/CaloTopology/interface/CaloTopology.h" + +#include "CondFormats/EcalObjects/interface/EcalPedestals.h" +#include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h" + +#include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h" + +#include + +/* + * DRNCorrectionProducerT + * + * Producer template generate a ValueMap of corrected energies and resolutions + * ValueMap contains std::pair of corrected energy, resolution + * + * Author: Simon Rothman (MIT) + * Written 2022 + * + */ + +namespace { + float sigmoid(float x) { return 1.0f / (1.0f + std::exp(-x)); } + + float logcorrection(float x) { + static float ln2 = std::log(2); + return ln2 * 2 * (sigmoid(x) - 0.5); + } + + //correction factor is transformed by sigmoid and "logratioflip target" + float correction(float x) { return std::exp(-logcorrection(x)); } + + inline float rescale(float x, float min, float range) { return (x - min) / range; } + + //resolution is transformed by softplus function + float resolution(float x) { return std::log(1 + std::exp(x)); } + + const float RHO_MIN = 0.0f; + const float RHO_RANGE = 13.0f; + + const float HOE_MIN = 0.0f; + const float HOE_RANGE = 0.05f; + + const float XY_MIN = -150.0f; + const float XY_RANGE = 300.0f; + + const float Z_MIN = -330.0f; + const float Z_RANGE = 660.0f; + + const float NOISE_MIN = 0.9f; + const float NOISE_RANGE = 3.0f; + + const float ECAL_MIN = 0.0f; + const float ECAL_RANGE = 250.0f; + + const float ES_MIN = 0.0f; + const float ES_RANGE = 0.1f; + +} // namespace + +template +class DRNCorrectionProducerT : public TritonEDProducer<> { +public: + explicit DRNCorrectionProducerT(const edm::ParameterSet& iConfig); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + void acquire(edm::Event const& iEvent, edm::EventSetup const& iSetup, Input& input) override; + void produce(edm::Event& iEvent, const edm::EventSetup& iSetup, Output const& iOutput) override; + +private: + edm::EDGetTokenT> particleToken_; + + edm::EDGetTokenT rhoToken_; + + edm::EDGetTokenT EBRecHitsToken_, EERecHitsToken_, ESRecHitsToken_; + + edm::ESGetToken pedToken_; + edm::ESGetToken geomToken_; + + size_t nPart_, nValidPart_; + + bool isEB(const T& part); + bool isEE(const T& part); + bool skip(const T& part); +}; + +template +DRNCorrectionProducerT::DRNCorrectionProducerT(const edm::ParameterSet& iConfig) + : TritonEDProducer<>(iConfig), + particleToken_(consumes(iConfig.getParameter("particleSource"))), + rhoToken_(consumes(iConfig.getParameter("rhoName"))), + EBRecHitsToken_(consumes(iConfig.getParameter("reducedEcalRecHitsEB"))), + EERecHitsToken_(consumes(iConfig.getParameter("reducedEcalRecHitsEE"))), + ESRecHitsToken_(consumes(iConfig.getParameter("reducedEcalRecHitsES"))), + pedToken_(esConsumes()), + geomToken_(esConsumes()) { + produces>>(); +} + +template +bool DRNCorrectionProducerT::isEB(const T& part) { + return part.superCluster()->seed()->hitsAndFractions().at(0).first.subdetId() == EcalBarrel; +} + +template +bool DRNCorrectionProducerT::isEE(const T& part) { + return part.superCluster()->seed()->hitsAndFractions().at(0).first.subdetId() == EcalEndcap; +} + +template +bool DRNCorrectionProducerT::skip(const T& part) { + /* + * Separated out from acquire() and produce() to ensure that skipping check is identical in both + * N.B. in MiniAOD there are sometimes particles with no RecHits + * We can not apply our regression to these, so we skip them + */ + return (!isEB(part) && !isEE(part)) || part.superCluster()->hitsAndFractions().empty(); +} + +template +void DRNCorrectionProducerT::acquire(edm::Event const& iEvent, edm::EventSetup const& iSetup, Input& iInput) { + /* + * Get products from event and event setup + */ + const auto& particles_ = iEvent.getHandle(particleToken_); + float rho = iEvent.get(rhoToken_); + + const auto& ped = &iSetup.getData(pedToken_); + const auto& geo = &iSetup.getData(geomToken_); + + const CaloSubdetectorGeometry* ecalEBGeom = + static_cast(geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel)); + const CaloSubdetectorGeometry* ecalEEGeom = + static_cast(geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap)); + const CaloSubdetectorGeometry* ecalESGeom = + static_cast(geo->getSubdetectorGeometry(DetId::Ecal, EcalPreshower)); + + const auto& recHitsEB = iEvent.get(EBRecHitsToken_); + const auto& recHitsEE = iEvent.get(EERecHitsToken_); + const auto& recHitsES = iEvent.get(ESRecHitsToken_); + + nPart_ = particles_->size(); + + if (nPart_ == 0) { + client_->setBatchSize(0); + return; + } else { + client_->setBatchSize(1); + } + + /* + * Determine how many particles, how many RecHits there are in each subdetector + */ + unsigned nHitsECAL = 0, nHitsES = 0; + nValidPart_ = 0; + for (auto& part : *particles_) { + const reco::SuperClusterRef& sc = part.superCluster(); + + if (skip(part)) + continue; + + nHitsECAL += sc->hitsAndFractions().size(); + + for (auto iES = sc->preshowerClustersBegin(); iES != sc->preshowerClustersEnd(); ++iES) { + nHitsES += (*iES)->hitsAndFractions().size(); + } + + ++nValidPart_; + } + + /* + * Allocate DRN inputs ({SB} is one of ECAL, ES): + * x{SB}: (x, y, z, energy, [noise]) continuous-valued inputs per RecHit + * f{SB}: (flagVal) integer denoting RecHit flag values + * gainECAL: (gain) integer in (0, 1, 2) denoting gain value + * gx: (rho, H/E) additional high-level features. + * batch{SB}: graph models require explicitely passing the particle index for each vertex + */ + auto& inputxECAL = iInput.at("xECAL"); + inputxECAL.setShape(0, nHitsECAL); + auto dataxECAL = inputxECAL.allocate(); + auto& vdataxECAL = (*dataxECAL)[0]; + + auto& inputfECAL = iInput.at("fECAL"); + inputfECAL.setShape(0, nHitsECAL); + auto datafECAL = inputfECAL.allocate(); + auto& vdatafECAL = (*datafECAL)[0]; + + auto& inputGainECAL = iInput.at("gainECAL"); + inputGainECAL.setShape(0, nHitsECAL); + auto dataGainECAL = inputGainECAL.allocate(); + auto& vdataGainECAL = (*dataGainECAL)[0]; + + auto& inputGx = iInput.at("gx"); + inputGx.setShape(0, nValidPart_); + auto dataGx = inputGx.allocate(); + auto& vdataGx = (*dataGx)[0]; + + auto& inputBatchECAL = iInput.at("batchECAL"); + inputBatchECAL.setShape(0, nHitsECAL); + auto dataBatchECAL = inputBatchECAL.allocate(); + auto& vdataBatchECAL = (*dataBatchECAL)[0]; + + auto& inputxES = iInput.at("xES"); + inputxES.setShape(0, nHitsES); + auto dataxES = inputxES.allocate(); + auto& vdataxES = (*dataxES)[0]; + + auto& inputfES = iInput.at("fES"); + inputfES.setShape(0, nHitsES); + auto datafES = inputfES.allocate(); + auto& vdatafES = (*datafES)[0]; + + auto& inputBatchES = iInput.at("batchES"); + inputBatchES.setShape(0, nHitsES); + auto dataBatchES = inputBatchES.allocate(); + auto& vdataBatchES = (*dataBatchES)[0]; + + /* + * Fill input tensors by iterating over particles... + */ + int64_t partNum = 0; + std::shared_ptr geom; + for (auto& part : *particles_) { + const reco::SuperClusterRef& sc = part.superCluster(); + + if (skip(part)) + continue; + + std::vector> hitsAndFractions = sc->hitsAndFractions(); + EcalRecHitCollection::const_iterator hit; + + //iterate over ECAL hits... + for (const auto& detitr : hitsAndFractions) { + DetId id = detitr.first.rawId(); + if (isEB(part)) { + geom = ecalEBGeom->getGeometry(id); + hit = recHitsEB.find(detitr.first); + } else { + geom = ecalEEGeom->getGeometry(id); + hit = recHitsEE.find(detitr.first); + } + + //fill xECAL + auto pos = geom->getPosition(); + vdataxECAL.push_back(rescale(pos.x(), XY_MIN, XY_RANGE)); + vdataxECAL.push_back(rescale(pos.y(), XY_MIN, XY_RANGE)); + vdataxECAL.push_back(rescale(pos.z(), Z_MIN, Z_RANGE)); + vdataxECAL.push_back(rescale(hit->energy() * detitr.second, ECAL_MIN, ECAL_RANGE)); + vdataxECAL.push_back(rescale(ped->find(detitr.first)->rms(1), NOISE_MIN, NOISE_RANGE)); + + //fill fECAL + int64_t flagVal = 0; + if (hit->checkFlag(EcalRecHit::kGood)) + flagVal += 1; + if (hit->checkFlag(EcalRecHit::kOutOfTime)) + flagVal += 2; + if (hit->checkFlag(EcalRecHit::kPoorCalib)) + flagVal += 4; + + vdatafECAL.push_back(flagVal); + + //fill gain + int64_t gainVal = 0; + if (hit->checkFlag(EcalRecHit::kHasSwitchToGain6)) + gainVal = 1; + else if (hit->checkFlag(EcalRecHit::kHasSwitchToGain1)) + gainVal = 0; + else + gainVal = 2; + + vdataGainECAL.push_back(gainVal); + + //fill batch number + vdataBatchECAL.push_back(partNum); + } //end iterate over ECAL hits + + //iterate over ES clusters... + for (auto iES = sc->preshowerClustersBegin(); iES != sc->preshowerClustersEnd(); ++iES) { + for (const auto& ESitr : (*iES)->hitsAndFractions()) { //iterate over ES hits + hit = recHitsES.find(ESitr.first); + geom = ecalESGeom->getGeometry(ESitr.first); + auto& pos = geom->getPosition(); + + //fill xES + vdataxES.push_back(rescale(pos.x(), XY_MIN, XY_RANGE)); + vdataxES.push_back(rescale(pos.y(), XY_MIN, XY_RANGE)); + vdataxES.push_back(rescale(pos.z(), Z_MIN, Z_RANGE)); + vdataxES.push_back(rescale(hit->energy(), ES_MIN, ES_RANGE)); + + //fill fES + int64_t flagVal = 0; + if (hit->checkFlag(EcalRecHit::kESGood)) + flagVal += 1; + + vdatafES.push_back(flagVal); + + //fill batchES + vdataBatchES.push_back(partNum); + } //end iterate over ES hits + } //end iterate over ES clusters + + //fill gx + vdataGx.push_back(rescale(rho, RHO_MIN, RHO_RANGE)); + vdataGx.push_back(rescale(part.hadronicOverEm(), HOE_MIN, HOE_RANGE)); + + //increment particle number + ++partNum; + } // end iterate over particles + + /* + * Convert input tensors to server data format + */ + inputxECAL.toServer(dataxECAL); + inputfECAL.toServer(datafECAL); + inputGainECAL.toServer(dataGainECAL); + inputBatchECAL.toServer(dataBatchECAL); + + inputGx.toServer(dataGx); + + inputxES.toServer(dataxES); + inputfES.toServer(datafES); + inputBatchES.toServer(dataBatchES); +} + +template +void DRNCorrectionProducerT::produce(edm::Event& iEvent, const edm::EventSetup& iSetup, Output const& iOutput) { + const auto& particles_ = iEvent.getHandle(particleToken_); + + std::vector> corrections; + corrections.reserve(nPart_); + + //if there are no particles, the fromServer() call will fail + //but we can just put() an empty valueMap + if (nPart_) { + const auto& muOut = iOutput.at("mu").fromServer(); + const auto& sigmaOut = iOutput.at("sigma").fromServer(); + + unsigned i = 0; + float mu, sigma, Epred, sigmaPred, rawE; + for (unsigned iPart = 0; iPart < nPart_; ++iPart) { + const auto& part = particles_->at(iPart); + if (!skip(part)) { + mu = correction(muOut[0][0 + 6 * i]); + sigma = resolution(sigmaOut[0][0 + 5 * i]); + ++i; + + rawE = part.superCluster()->rawEnergy(); + Epred = mu * rawE; + sigmaPred = sigma * rawE; + corrections.emplace_back(Epred, sigmaPred); + } else { + corrections.emplace_back(-1.0f, -1.0f); + } + } + } + + //fill + auto out = std::make_unique>>(); + edm::ValueMap>::Filler filler(*out); + filler.insert(particles_, corrections.begin(), corrections.end()); + filler.fill(); + + iEvent.put(std::move(out)); +} + +template +void DRNCorrectionProducerT::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + TritonClient::fillPSetDescription(desc); + desc.add("particleSource"); + desc.add("rhoName"); + desc.add("reducedEcalRecHitsEB", edm::InputTag("reducedEcalRecHitsEB")); + desc.add("reducedEcalRecHitsEE", edm::InputTag("reducedEcalRecHitsEE")); + desc.add("reducedEcalRecHitsES", edm::InputTag("reducedEcalRecHitsES")); + descriptions.addWithDefaultLabel(desc); +} + +//reco:: template instances are supported +//uncomment the lines below to enable them + +//using GsfElectronDRNCorrectionProducer = DRNCorrectionProducerT; +//using GedPhotonDRNCorrectionProducer = DRNCorrectionProducerT; +//DEFINE_FWK_MODULE(GedPhotonDRNCorrectionProducer); +//DEFINE_FWK_MODULE(GsfElectronDRNCorrectionProducer); + +using PatElectronDRNCorrectionProducer = DRNCorrectionProducerT; +using PatPhotonDRNCorrectionProducer = DRNCorrectionProducerT; + +DEFINE_FWK_MODULE(PatPhotonDRNCorrectionProducer); +DEFINE_FWK_MODULE(PatElectronDRNCorrectionProducer); diff --git a/RecoEgamma/EgammaTools/plugins/EGRegressionModifierDRN.cc b/RecoEgamma/EgammaTools/plugins/EGRegressionModifierDRN.cc new file mode 100644 index 0000000000000..a6f54a847e7b2 --- /dev/null +++ b/RecoEgamma/EgammaTools/plugins/EGRegressionModifierDRN.cc @@ -0,0 +1,218 @@ +#include "CommonTools/CandAlgos/interface/ModifyObjectValueBase.h" + +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Utilities/interface/EDGetToken.h" + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/ESHandle.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 "DataFormats/Common/interface/ValueMap.h" + +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/EgammaCandidates/interface/Photon.h" + +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" + +#include "RecoEgamma/EgammaTools/interface/EgammaRegressionContainer.h" +#include "RecoEgamma/EgammaTools/interface/EpCombinationTool.h" +#include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h" + +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" + +/* + * EGRegressionModifierDRN + * + * Object modifier to apply DRN regression. + * Designed to be a drop-in replacement for EGRegressionModifierVX + * + * Requires the appropriate DRNCorrectionProducerX(s) to also be in the path + * You can specify which of reco::GsfElectron, reco::Photon, pat::Electron, pat::Photon + * to apply corrections to in the config + * + */ + +class EGRegressionModifierDRN : public ModifyObjectValueBase { +public: + EGRegressionModifierDRN(const edm::ParameterSet& conf, edm::ConsumesCollector& cc); + + void setEvent(const edm::Event&) final; + void setEventContent(const edm::EventSetup&) final; + + void modifyObject(reco::GsfElectron&) const final; + void modifyObject(reco::Photon&) const final; + + void modifyObject(pat::Electron&) const final; + void modifyObject(pat::Photon&) const final; + +private: + template + struct partVars { + edm::InputTag source; + edm::EDGetTokenT> token; + const edm::View* particles; + + edm::InputTag correctionsSource; + edm::EDGetTokenT>> correctionsToken; + const edm::ValueMap>* corrections; + + bool userFloat; + std::string energyFloat, resFloat; + + unsigned i; + + partVars(const edm::ParameterSet& config, edm::ConsumesCollector& cc) { + source = config.getParameter("source"); + token = cc.consumes(source); + + correctionsSource = config.getParameter("correctionsSource"); + correctionsToken = cc.consumes(correctionsSource); + + if (config.exists("energyFloat")) { + userFloat = true; + energyFloat = config.getParameter("energyFloat"); + resFloat = config.getParameter("resFloat"); + } else { + userFloat = false; + } + + i = 0; + } + + const std::pair getCorrection(T& part); + + const void doUserFloat(T& part, const std::pair& correction) const { + part.addUserFloat(energyFloat, correction.first); + part.addUserFloat(resFloat, correction.second); + } + }; + + std::unique_ptr> patPhotons_; + std::unique_ptr> patElectrons_; + std::unique_ptr> gedPhotons_; + std::unique_ptr> gsfElectrons_; +}; + +EGRegressionModifierDRN::EGRegressionModifierDRN(const edm::ParameterSet& conf, edm::ConsumesCollector& cc) + : ModifyObjectValueBase(conf) { + if (conf.exists("patPhotons")) { + patPhotons_ = std::make_unique>(conf.getParameterSet("patPhotons"), cc); + } + + if (conf.exists("gedPhotons")) { + gedPhotons_ = std::make_unique>(conf.getParameterSet("gedPhotons"), cc); + } + + if (conf.exists("patElectrons")) { + patElectrons_ = std::make_unique>(conf.getParameterSet("patElectrons"), cc); + } + + if (conf.exists("gsfElectrons")) { + gsfElectrons_ = std::make_unique>(conf.getParameterSet("gsfElectrons"), cc); + } +} + +void EGRegressionModifierDRN::setEvent(const edm::Event& evt) { + if (patElectrons_) { + patElectrons_->particles = &evt.get(patElectrons_->token); + patElectrons_->corrections = &evt.get(patElectrons_->correctionsToken); + patElectrons_->i = 0; + } + + if (patPhotons_) { + patPhotons_->particles = &evt.get(patPhotons_->token); + patPhotons_->corrections = &evt.get(patPhotons_->correctionsToken); + patPhotons_->i = 0; + } + + if (gsfElectrons_) { + gsfElectrons_->particles = &evt.get(gsfElectrons_->token); + gsfElectrons_->corrections = &evt.get(gsfElectrons_->correctionsToken); + gsfElectrons_->i = 0; + } + + if (gedPhotons_) { + gedPhotons_->particles = &evt.get(gedPhotons_->token); + gedPhotons_->corrections = &evt.get(gedPhotons_->correctionsToken); + gedPhotons_->i = 0; + } +} + +void EGRegressionModifierDRN::setEventContent(const edm::EventSetup& iSetup) {} + +void EGRegressionModifierDRN::modifyObject(reco::GsfElectron& ele) const { + if (!gsfElectrons_) + return; + + const std::pair& correction = gsfElectrons_->getCorrection(ele); + + if (correction.first > 0 && correction.second > 0) { + ele.setCorrectedEcalEnergy(correction.first, true); + ele.setCorrectedEcalEnergyError(correction.second); + } + + throw cms::Exception("EGRegressionModifierDRN") + << "Electron energy corrections not fully implemented yet:" << std::endl + << "Still need E/p combination" << std::endl + << "Do not enable DRN for electrons" << std::endl; +} + +void EGRegressionModifierDRN::modifyObject(pat::Electron& ele) const { + if (!patElectrons_) + return; + + const std::pair& correction = patElectrons_->getCorrection(ele); + + if (patElectrons_->userFloat) { + patElectrons_->doUserFloat(ele, correction); + } else if (correction.first > 0 && correction.second > 0) { + ele.setCorrectedEcalEnergy(correction.first, true); + ele.setCorrectedEcalEnergyError(correction.second); + } + + throw cms::Exception("EGRegressionModifierDRN") + << "Electron energy corrections not fully implemented yet:" << std::endl + << "Still need E/p combination" << std::endl + << "Do not enable DRN for electrons" << std::endl; +} + +void EGRegressionModifierDRN::modifyObject(pat::Photon& pho) const { + if (!patPhotons_) + return; + + const std::pair& correction = patPhotons_->getCorrection(pho); + + if (patPhotons_->userFloat) { + patPhotons_->doUserFloat(pho, correction); + } else if (correction.first > 0 && correction.second > 0) { + pho.setCorrectedEnergy(pat::Photon::P4type::regression2, correction.first, correction.second, true); + } +} + +void EGRegressionModifierDRN::modifyObject(reco::Photon& pho) const { + if (!gedPhotons_) + return; + + const std::pair& correction = gedPhotons_->getCorrection(pho); + + if (correction.first > 0 && correction.second > 0) { + pho.setCorrectedEnergy(reco::Photon::P4type::regression2, correction.first, correction.second, true); + } +}; + +template +const std::pair EGRegressionModifierDRN::partVars::getCorrection(T& part) { + edm::Ptr ptr = particles->ptrAt(i++); + + std::pair correction = (*corrections)[ptr]; + + return correction; +} + +DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, EGRegressionModifierDRN, "EGRegressionModifierDRN"); diff --git a/RecoEgamma/EgammaTools/python/egammaObjectModificationsInMiniAOD_cff.py b/RecoEgamma/EgammaTools/python/egammaObjectModificationsInMiniAOD_cff.py index 351cb21abf91c..848b9a16c182b 100644 --- a/RecoEgamma/EgammaTools/python/egammaObjectModificationsInMiniAOD_cff.py +++ b/RecoEgamma/EgammaTools/python/egammaObjectModificationsInMiniAOD_cff.py @@ -145,6 +145,16 @@ def setup_mva(val_pset,cat_pset,prod_name,mva_name): ) ) +photonDRNModifier = cms.PSet( + modifierName = cms.string("EGRegressionModifierDRN"), + patPhotons = cms.PSet( + source = cms.InputTag("selectedPatPhotons"), + correctionsSource = cms.InputTag('patPhotonsDRN'), + energyFloat = cms.string("energyDRN"), + resFloat = cms.string("resolutionDRN") + ) + ) + def appendReducedEgammaEnergyScaleAndSmearingModifier(modifiers): modifiers.append(reducedEgammaEnergyScaleAndSmearingModifier) @@ -158,6 +168,9 @@ def appendEgamma8XLegacyAppendableModifiers (modifiers): def appendEgammaHIPhotonIsolationModifier(modifiers): modifiers.append(egammaHIPhotonIsolationModifier) +def appendPhotonDRNModifier(modifiers): + modifiers.append(photonDRNModifier) + from Configuration.Eras.Modifier_run2_miniAOD_94XFall17_cff import run2_miniAOD_94XFall17 from Configuration.ProcessModifiers.run2_miniAOD_UL_cff import run2_miniAOD_UL (run2_miniAOD_94XFall17 | run2_miniAOD_UL).toModify(egamma_modifications,appendReducedEgammaEnergyScaleAndSmearingModifier) @@ -170,3 +183,6 @@ def appendEgammaHIPhotonIsolationModifier(modifiers): from Configuration.ProcessModifiers.pp_on_AA_cff import pp_on_AA pp_on_AA.toModify(egamma_modifications, appendEgammaHIPhotonIsolationModifier) + +from Configuration.ProcessModifiers.photonDRN_cff import _photonDRN +_photonDRN.toModify(egamma_modifications, appendPhotonDRNModifier)