diff --git a/RecoEcal/EgammaClusterAlgos/BuildFile.xml b/RecoEcal/EgammaClusterAlgos/BuildFile.xml index 8ea2b99e359c2..6334a637b8a2f 100644 --- a/RecoEcal/EgammaClusterAlgos/BuildFile.xml +++ b/RecoEcal/EgammaClusterAlgos/BuildFile.xml @@ -14,6 +14,7 @@ + diff --git a/RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorDRN.h b/RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorDRN.h new file mode 100644 index 0000000000000..e0d76ef6f79dc --- /dev/null +++ b/RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorDRN.h @@ -0,0 +1,90 @@ +//-------------------------------------------------------------------------------------------------- +// +// SCEnergyCorrectorDRN +// +// Helper Class for applying regression-based energy corrections with DRN implimentation +// +// Based on RecoEcal/EgammaClusterAlgos/SCEnergyCorrectorSemiParm +// +// Author: Simon Rothman (MIT, UMN) +// +//-------------------------------------------------------------------------------------------------- + +#ifndef RecoEcal_EgammaClusterAlgos_SCEnergyCorrectorDRN_h +#define RecoEcal_EgammaClusterAlgos_SCEnergyCorrectorDRN_h + +#include "HeterogeneousCore/SonicTriton/interface/TritonData.h" + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "Geometry/Records/interface/CaloTopologyRecord.h" +#include "Geometry/CaloTopology/interface/CaloTopology.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" + +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHit.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" +#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" + +#include "CondFormats/GBRForest/interface/GBRForestD.h" +#include "CondFormats/DataRecord/interface/GBRDWrapperRcd.h" + +#include "RecoEgamma/EgammaTools/interface/EgammaBDTOutputTransformer.h" +#include "RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h" + +#include +#include +#include +#include + +class SCEnergyCorrectorDRN { +public: + SCEnergyCorrectorDRN(); + //if you want override the default on where conditions are consumed, you need to use + //the other constructor and then call setTokens approprately + SCEnergyCorrectorDRN(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc); + + static void fillPSetDescription(edm::ParameterSetDescription& desc); + static edm::ParameterSetDescription makePSetDescription(); + + template + void setTokens(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc); + + void setEventSetup(const edm::EventSetup& es); + void setEvent(const edm::Event& e); + + void makeInput(const edm::Event& iEvent, TritonInputMap& iInput, const reco::SuperClusterCollection& inputSCs) const; + TritonOutput getOutput(const TritonOutputMap& iOutput); + +private: + const CaloTopology* caloTopo_; + const CaloGeometry* caloGeom_; + edm::ESGetToken caloTopoToken_; + edm::ESGetToken caloGeomToken_; + + edm::EDGetTokenT tokenEBRecHits_; + edm::EDGetTokenT tokenEERecHits_; + edm::EDGetTokenT rhoToken_; + + edm::Handle recHitsEB_; + edm::Handle recHitsEE_; + + edm::Handle rhoHandle_; +}; + +template +void SCEnergyCorrectorDRN::setTokens(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc) { + tokenEBRecHits_ = cc.consumes(iConfig.getParameter("ecalRecHitsEB")); + tokenEERecHits_ = cc.consumes(iConfig.getParameter("ecalRecHitsEE")); + caloGeomToken_ = cc.esConsumes(); + caloTopoToken_ = cc.esConsumes(); + rhoToken_ = cc.consumes(iConfig.getParameter("rhoFastJet")); +} +#endif diff --git a/RecoEcal/EgammaClusterAlgos/src/SCEnergyCorrectorDRN.cc b/RecoEcal/EgammaClusterAlgos/src/SCEnergyCorrectorDRN.cc new file mode 100644 index 0000000000000..8bebda2d5b419 --- /dev/null +++ b/RecoEcal/EgammaClusterAlgos/src/SCEnergyCorrectorDRN.cc @@ -0,0 +1,122 @@ +#include "RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorDRN.h" + +#include "FWCore/Utilities/interface/isFinite.h" +#include "FWCore/Utilities/interface/Transition.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/Math/interface/deltaPhi.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalTools.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" +#include "RecoEgamma/EgammaTools/interface/EgammaHGCALIDParamDefaults.h" + +#include + +static const float RHO_MAX = 15.0f; +static const float X_MAX = 150.0f; +static const float X_RANGE = 300.0f; +static const float Y_MAX = 150.0f; +static const float Y_RANGE = 300.0f; +static const float Z_MAX = 330.0f; +static const float Z_RANGE = 660.0f; +static const float E_RANGE = 250.0f; + +SCEnergyCorrectorDRN::SCEnergyCorrectorDRN() : caloTopo_(nullptr), caloGeom_(nullptr) {} + +SCEnergyCorrectorDRN::SCEnergyCorrectorDRN(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc) + : SCEnergyCorrectorDRN() { + setTokens(iConfig, cc); +} + +void SCEnergyCorrectorDRN::fillPSetDescription(edm::ParameterSetDescription& desc) { + desc.add("ecalRecHitsEE", edm::InputTag("ecalRecHit", "reducedEcalRecHitsEE")); + desc.add("ecalRecHitsEB", edm::InputTag("ecalRecHit", "reducedEcalRecHitsEB")); + desc.add("rhoFastJet", edm::InputTag("fixedGridRhoAll")); +} + +edm::ParameterSetDescription SCEnergyCorrectorDRN::makePSetDescription() { + edm::ParameterSetDescription desc; + fillPSetDescription(desc); + return desc; +} + +void SCEnergyCorrectorDRN::setEventSetup(const edm::EventSetup& es) { + caloTopo_ = &es.getData(caloTopoToken_); + caloGeom_ = &es.getData(caloGeomToken_); +} + +void SCEnergyCorrectorDRN::setEvent(const edm::Event& event) { + event.getByToken(tokenEBRecHits_, recHitsEB_); + event.getByToken(tokenEERecHits_, recHitsEE_); + event.getByToken(rhoToken_, rhoHandle_); +} + +void SCEnergyCorrectorDRN::makeInput(const edm::Event& iEvent, + TritonInputMap& iInput, + const reco::SuperClusterCollection& inputSCs) const { + std::vector nHits; + nHits.reserve(inputSCs.size()); + unsigned totalHits = 0; + unsigned n; + for (const auto& inputSC : inputSCs) { + n = inputSC.hitsAndFractions().size(); + totalHits += n; + nHits.push_back(n); + } + + //set shapes + auto& input1 = iInput.at("x__0"); + input1.setShape(0, totalHits); + auto data1 = input1.allocate(); + auto& vdata1 = (*data1)[0]; + + auto& input2 = iInput.at("batch__1"); + input2.setShape(0, totalHits); + auto data2 = input2.allocate(); + auto& vdata2 = (*data2)[0]; + + auto& input3 = iInput.at("graphx__2"); + input3.setShape(0, 2 * nHits.size()); + auto data3 = input3.allocate(); + auto& vdata3 = (*data3)[0]; + + //fill + unsigned batchNum = 0; + float En, frac, x, y, z; + for (const auto& inputSC : inputSCs) { + const auto& hits = inputSC.hitsAndFractions(); + const bool isEB = hits[0].first.subdetId() == EcalBarrel; + const auto& recHitsProduct = isEB ? recHitsEB_.product() : recHitsEE_.product(); + for (const auto& hit : hits) { + En = EcalClusterTools::recHitEnergy(hit.first, recHitsProduct); + frac = hit.second; + GlobalPoint position = caloGeom_->getGeometry(hit.first)->getPosition(); + x = (position.x() + X_MAX) / X_RANGE; + y = (position.y() + Y_MAX) / Y_RANGE; + z = (position.z() + Z_MAX) / Z_RANGE; + vdata1.push_back(x); + vdata1.push_back(y); + vdata1.push_back(z); + vdata1.push_back(En * frac / E_RANGE); + //Triton does not currently support batching for pytorch GNNs + //We pass batch indices explicitely + vdata2.push_back(batchNum); + } + vdata3.push_back(*rhoHandle_ / RHO_MAX); + vdata3.push_back(0.0); + ++batchNum; + } + + // convert to server format + input1.toServer(data1); + input2.toServer(data2); + input3.toServer(data3); +} + +TritonOutput SCEnergyCorrectorDRN::getOutput(const TritonOutputMap& iOutput) { + //check the results + const auto& output1 = iOutput.begin()->second; + // convert from server format + const auto& serverout = output1.fromServer(); + + return serverout; +} diff --git a/RecoEcal/EgammaClusterProducers/BuildFile.xml b/RecoEcal/EgammaClusterProducers/BuildFile.xml index 755900bc4d8c4..6dfa027a4712c 100644 --- a/RecoEcal/EgammaClusterProducers/BuildFile.xml +++ b/RecoEcal/EgammaClusterProducers/BuildFile.xml @@ -1,4 +1,5 @@ + diff --git a/RecoEcal/EgammaClusterProducers/python/SCEnergyCorrectorDRNProducer_cfi.py b/RecoEcal/EgammaClusterProducers/python/SCEnergyCorrectorDRNProducer_cfi.py new file mode 100644 index 0000000000000..203f1dd49be78 --- /dev/null +++ b/RecoEcal/EgammaClusterProducers/python/SCEnergyCorrectorDRNProducer_cfi.py @@ -0,0 +1,26 @@ +import FWCore.ParameterSet.Config as cms + +DRNProducerEB = cms.EDProducer('SCEnergyCorrectorDRNProducer', + inputSCs = cms.InputTag('particleFlowSuperClusterECAL','particleFlowSuperClusterECALBarrel'), + Client = cms.PSet( + mode = cms.string("Async"), + modelName = cms.string("MustacheEB"), + modelConfigPath = cms.FileInPath("RecoEcal/EgammaClusterProducers/data/models/MustacheEB/config.pbtxt"), + allowedTries = cms.untracked.uint32(1), + timeout = cms.untracked.uint32(10), + ), + ) + + +DRNProducerEE = cms.EDProducer('SCEnergyCorrectorDRNProducer', + inputSCs = cms.InputTag('particleFlowSuperClusterECAL','particleFlowSuperClusterECALEndcapWithPreshower'), + Client = cms.PSet( + mode = cms.string("Async"), + modelName = cms.string('MustacheEE'), + modelConfigPath = cms.FileInPath("RecoEcal/EgammaClusterProducers/data/models/MustacheEE/config.pbtxt"), + allowedTries = cms.untracked.uint32(1), + timeout = cms.untracked.uint32(10), + ), + ) + + diff --git a/RecoEcal/EgammaClusterProducers/src/SCEnergyCorrectorDRNProducer.cc b/RecoEcal/EgammaClusterProducers/src/SCEnergyCorrectorDRNProducer.cc new file mode 100644 index 0000000000000..425e354932de5 --- /dev/null +++ b/RecoEcal/EgammaClusterProducers/src/SCEnergyCorrectorDRNProducer.cc @@ -0,0 +1,113 @@ +#include "HeterogeneousCore/SonicTriton/interface/TritonEDProducer.h" +#include "HeterogeneousCore/SonicTriton/interface/TritonData.h" + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h" +#include "DataFormats/Common/interface/ValueMap.h" + +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" + +#include "RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorDRN.h" + +#include +#include +#include +#include + +/* + * SCEnergyCorrectorDRNProducer + * + * Simple producer to generate a set of corrected superclusters with the DRN regression + * Based on RecoEcal/EgammaClusterProducers/SCEnergyCorrectorProducer by S. Harper (RAL/CERN) + * + * Author: Simon Rothman (UMN, MIT) + * + */ + +namespace { + float sigmoid(float x) { return 1.0f / (1.0f + exp(-x)); } + + float logcorrection(float x) { + static float ln2 = log(2); + return ln2 * 2 * (sigmoid(x) - 0.5); + } + + float correction(float x) { return exp(-logcorrection(x)); } +} // namespace + +class SCEnergyCorrectorDRNProducer : public TritonEDProducer<> { +public: + explicit SCEnergyCorrectorDRNProducer(const edm::ParameterSet& iConfig); + + void beginLuminosityBlock(const edm::LuminosityBlock& iLumi, const edm::EventSetup& iSetup) override; + + 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; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + SCEnergyCorrectorDRN energyCorrector_; + edm::EDGetTokenT inputSCToken_; +}; + +SCEnergyCorrectorDRNProducer::SCEnergyCorrectorDRNProducer(const edm::ParameterSet& iConfig) + : TritonEDProducer<>(iConfig, "SCEnergyCorrectorDRNProducer"), + energyCorrector_(iConfig.getParameterSet("correctorCfg"), consumesCollector()), + inputSCToken_(consumes(iConfig.getParameter("inputSCs"))) { + produces(); +} + +void SCEnergyCorrectorDRNProducer::beginLuminosityBlock(const edm::LuminosityBlock& iLumi, + const edm::EventSetup& iSetup) { + energyCorrector_.setEventSetup(iSetup); +} + +void SCEnergyCorrectorDRNProducer::acquire(edm::Event const& iEvent, edm::EventSetup const& iSetup, Input& iInput) { + const auto& inputSCs = iEvent.get(inputSCToken_); + + if (inputSCs.empty()) { + client_->setBatchSize(0); + return; + } else { + client_->setBatchSize(1); + } + + energyCorrector_.setEvent(iEvent); + energyCorrector_.makeInput(iEvent, iInput, inputSCs); +} + +void SCEnergyCorrectorDRNProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup, Output const& iOutput) { + const auto& inputSCs = iEvent.get(inputSCToken_); + if (inputSCs.empty()) + return; + + const auto& serverout = energyCorrector_.getOutput(iOutput); + + auto corrSCs = std::make_unique(); + unsigned i = 0; + for (const auto& inputSC : inputSCs) { + float corrEn = correction(serverout[0][0 + 6 * i]) * inputSC.rawEnergy(); + corrSCs->push_back(inputSC); + corrSCs->back().setEnergy(corrEn); + corrSCs->back().setCorrectedEnergy(corrEn); + ++i; + } + + auto scHandle = iEvent.put(std::move(corrSCs)); +} + +void SCEnergyCorrectorDRNProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("correctorCfg", SCEnergyCorrectorDRN::makePSetDescription()); + TritonClient::fillPSetDescription(desc); + desc.add("inputSCs", edm::InputTag("particleFlowSuperClusterECAL")); + descriptions.add("scEnergyCorrectorDRNProducer", desc); +} + +DEFINE_FWK_MODULE(SCEnergyCorrectorDRNProducer); diff --git a/RecoEcal/EgammaClusterProducers/test/BuildFile.xml b/RecoEcal/EgammaClusterProducers/test/BuildFile.xml index 2d5cdb7c4f420..669e5257fd847 100644 --- a/RecoEcal/EgammaClusterProducers/test/BuildFile.xml +++ b/RecoEcal/EgammaClusterProducers/test/BuildFile.xml @@ -8,3 +8,13 @@ + + + + + + + + + + diff --git a/RecoEcal/EgammaClusterProducers/test/DRNTest_cfg.py b/RecoEcal/EgammaClusterProducers/test/DRNTest_cfg.py new file mode 100644 index 0000000000000..98b6bd7e7d13f --- /dev/null +++ b/RecoEcal/EgammaClusterProducers/test/DRNTest_cfg.py @@ -0,0 +1,81 @@ +import FWCore.ParameterSet.Config as cms +import FWCore.ParameterSet.VarParsing as VarParsing + +from Configuration.ProcessModifiers.enableSonicTriton_cff import enableSonicTriton +process = cms.Process("Demo",enableSonicTriton) + +process.load("HeterogeneousCore.SonicTriton.TritonService_cff") +process.load("FWCore.MessageService.MessageLogger_cfi") +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000) ) +process.options.numberOfThreads = cms.untracked.uint32(4) +process.options.numberOfStreams = cms.untracked.uint32(4) +process.options.SkipEvent = cms.untracked.vstring('ProductNotFound') + +process.TritonService.verbose = False +process.TritonService.fallback.verbose = False +process.TritonService.fallback.useDocker = False +process.TritonService.fallback.useGPU = False + +process.MessageLogger.cerr.FwkReport.reportEvery = cms.untracked.int32(100) +#process.MessageLogger.suppressWarning = cms.untracked.vstring('DRNProducerEB', 'DRNProducerEE') + +process.source = cms.Source("PoolSource", + # replace 'myfile.root' with the source file you want to use + fileNames = cms.untracked.vstring( +#'root://cms-xrd-global.cern.ch///store/mc/RunIISummer20UL18RECO/DoubleElectron_Pt-1To300-gun/AODSIM/FlatPU0to70EdalIdealGT_EdalIdealGT_106X_upgrade2018_realistic_v11_L1v1_EcalIdealIC-v2/270000/4CDD9457-E14C-D84A-9BD4-3140CB6AEEB6.root' +'root://cms-xrd-global.cern.ch///store/mc/RunIISummer20UL18RECO/DoubleElectron_Pt-1To300-gun/AODSIM/FlatPU0to70EdalIdealGT_EdalIdealGT_106X_upgrade2018_realistic_v11_L1v1_EcalIdealIC-v2/270000/4CDD9457-E14C-D84A-9BD4-3140CB6AEEB6.root', +'root://cms-xrd-global.cern.ch///store/mc/RunIISummer20UL18RECO/DoubleElectron_Pt-1To300-gun/AODSIM/FlatPU0to70EdalIdealGT_EdalIdealGT_106X_upgrade2018_realistic_v11_L1v1_EcalIdealIC-v2/270000/D3A06456-7F7D-C940-9E56-0DCA06B3ECC9.root', +'root://cms-xrd-global.cern.ch///store/mc/RunIISummer20UL18RECO/DoubleElectron_Pt-1To300-gun/AODSIM/FlatPU0to70EdalIdealGT_EdalIdealGT_106X_upgrade2018_realistic_v11_L1v1_EcalIdealIC-v2/270000/477810A1-AE5C-6A49-8067-35776F0C78B6.root', +'root://cms-xrd-global.cern.ch///store/mc/RunIISummer20UL18RECO/DoubleElectron_Pt-1To300-gun/AODSIM/FlatPU0to70EdalIdealGT_EdalIdealGT_106X_upgrade2018_realistic_v11_L1v1_EcalIdealIC-v2/270000/99F8FA99-B120-BF40-BD67-7825696B9E78.root' + + ) + ) + +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") + +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, '106X_upgrade2018_realistic_v11_Ecal5', '') +#from RecoEcal.EgammaClusterProducers.SCEnergyCorrectorDRNProducer_cfi import * +from RecoEcal.EgammaClusterProducers.SCEnergyCorrectorDRNProducer_cfi import * + +process.DRNProducerEB = DRNProducerEB +process.DRNProducerEE = DRNProducerEE + +from PhysicsTools.SelectorUtils.tools.vid_id_tools import * +dataFormat = DataFormat.AOD +switchOnVIDElectronIdProducer(process, dataFormat) + +# define which IDs we want to produce +my_id_modules = ['RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Fall17_94X_V2_cff'] + +for idmod in my_id_modules: + setupAllVIDIdsInModule(process, idmod, setupVIDElectronSelection) + +process.nTuplelize = cms.EDAnalyzer('DRNTestNTuplizer', + vertexCollection = cms.InputTag('offlinePrimaryVertices'), + rhoFastJet = cms.InputTag("fixedGridRhoFastjetAll"), + pileupInfo = cms.InputTag("addPileupInfo"), + electrons = cms.InputTag("gedGsfElectrons"), + genParticles = cms.InputTag("genParticles"), + #Cut Based Id + eleLooseIdMap = cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Fall17-94X-V2-loose"), + eleMediumIdMap = cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Fall17-94X-V2-medium"), + eleTightIdMap = cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Fall17-94X-V2-tight"), + + tracks = cms.InputTag("globalTracks"), + SkipEvent = cms.untracked.vstring('ProductNotFound') + ) + +process.TFileService = cms.Service("TFileService", + fileName = cms.string("nTupleMC.root"), + closeFileFast = cms.untracked.bool(True) + ) + +process.load( "HLTrigger.Timer.FastTimerService_cfi" ) + +#process.p = cms.Path(process.DRNProducerEB*process.egmGsfElectronIDSequence*process.nTuplelize) + +#process.p = cms.Path(process.DRNProducerEB*process.DRNProducerEE*process.egmGsfElectronIDSequence*process.nTuplelize) +process.p = cms.Path(process.DRNProducerEB*process.DRNProducerEE) diff --git a/RecoEcal/EgammaClusterProducers/test/runtestRecoEcalEgammaClusterProducers.cpp b/RecoEcal/EgammaClusterProducers/test/runtestRecoEcalEgammaClusterProducers.cpp new file mode 100644 index 0000000000000..b2991bd18ae57 --- /dev/null +++ b/RecoEcal/EgammaClusterProducers/test/runtestRecoEcalEgammaClusterProducers.cpp @@ -0,0 +1,3 @@ +#include "FWCore/Utilities/interface/TestHelper.h" + +RUNTEST() diff --git a/RecoEcal/EgammaClusterProducers/test/runtests.sh b/RecoEcal/EgammaClusterProducers/test/runtests.sh new file mode 100755 index 0000000000000..b53b3c782a60e --- /dev/null +++ b/RecoEcal/EgammaClusterProducers/test/runtests.sh @@ -0,0 +1,5 @@ +#!/bin/sh + +function die { echo $1: status $2 ; exit $2; } + +cmsRun ${LOCAL_TEST_DIR}/DRNTest_cfg.py || die 'Failure using cmsRun' $?