diff --git a/RecoCaloTools/Navigation/interface/CaloRectangle.h b/RecoCaloTools/Navigation/interface/CaloRectangle.h index 2784c5ead4a07..6f5dc7c8f1b06 100644 --- a/RecoCaloTools/Navigation/interface/CaloRectangle.h +++ b/RecoCaloTools/Navigation/interface/CaloRectangle.h @@ -4,20 +4,11 @@ #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h" #include "Geometry/CaloTopology/interface/CaloTopology.h" -template -T offsetBy(T start, CaloSubdetectorTopology const& topo, int dX, int dY) -{ - for(int x = 0; x < std::abs(dX) && start != T(0); x++) { - // east is eta in barrel - start = dX > 0 ? topo.goEast(start) : topo.goWest(start); - } - - for(int y = 0; y < std::abs(dY) && start != T(0); y++) { - // north is phi in barrel - start = dY > 0 ? topo.goNorth(start) : topo.goSouth(start); - } - return start; -} +/* + * CaloRectangle is a class to create a rectangular range of DetIds around a + * central DetId. Meant to be used for range-based loops to calculate cluster + * shape variables. + */ struct CaloRectangle { const int ixMin; @@ -38,6 +29,21 @@ inline auto makeCaloRectangle(int sizeX, int sizeY) { return CaloRectangle{-sizeX, sizeX, -sizeY, sizeY}; } +template +T offsetBy(T start, CaloSubdetectorTopology const& topo, int dX, int dY) +{ + for(int x = 0; x < std::abs(dX) && start != T(0); x++) { + // east is eta in barrel + start = dX > 0 ? topo.goEast(start) : topo.goWest(start); + } + + for(int y = 0; y < std::abs(dY) && start != T(0); y++) { + // north is phi in barrel + start = dY > 0 ? topo.goNorth(start) : topo.goSouth(start); + } + return start; +} + template class CaloRectangleRange { diff --git a/RecoEgamma/ElectronIdentification/plugins/ElectronMVANtuplizer.cc b/RecoEgamma/ElectronIdentification/plugins/ElectronMVANtuplizer.cc index fed42e5271737..a96e744789aa9 100644 --- a/RecoEgamma/ElectronIdentification/plugins/ElectronMVANtuplizer.cc +++ b/RecoEgamma/ElectronIdentification/plugins/ElectronMVANtuplizer.cc @@ -167,8 +167,8 @@ ElectronMVANtuplizer::ElectronMVANtuplizer(const edm::ParameterSet& iConfig) , vertices_ (src_, consumesCollector(), iConfig, "vertices" , "verticesMiniAOD") , pileup_ (src_, consumesCollector(), iConfig, "pileup" , "pileupMiniAOD") , genParticles_ (src_, consumesCollector(), iConfig, "genParticles", "genParticlesMiniAOD") - , ebRecHits_ (src_, consumesCollector(), iConfig, "ebReducedRecHitCollection", "ebReducedRecHitCollectionMiniAOD") - , eeRecHits_ (src_, consumesCollector(), iConfig, "eeReducedRecHitCollection", "eeReducedRecHitCollectionMiniAOD") + , ebRecHits_ (src_, consumesCollector(), iConfig, "ebReducedRecHitCollection", "ebReducedRecHitCollectionMiniAOD") + , eeRecHits_ (src_, consumesCollector(), iConfig, "eeReducedRecHitCollection", "eeReducedRecHitCollectionMiniAOD") , mvaPasses_ (nEleMaps_) , mvaValues_ (nValMaps_) , mvaCats_ (nCats_) @@ -206,13 +206,11 @@ ElectronMVANtuplizer::ElectronMVANtuplizer(const edm::ParameterSet& iConfig) tree_->Branch("ele_q",&eleQ_); tree_->Branch("ele_3q",&ele3Q_); - if (doEnergyMatrix_) tree_->Branch("ele_energyMatrix",&energyMatrix_); + if (doEnergyMatrix_) tree_->Branch("energyMatrix",&energyMatrix_); if (isMC_) tree_->Branch("matchedToGenEle", &matchedToGenEle_); - for (int i = 0; i < nVars_; ++i) { - tree_->Branch(mvaVarMngr_.getName(i).c_str(), &vars_[i]); - } + for (int i = 0; i < nVars_; ++i) tree_->Branch(mvaVarMngr_.getName(i).c_str(), &vars_[i]); tree_->Branch("ele_isEB",&eleIsEB_); tree_->Branch("ele_isEE",&eleIsEE_); @@ -252,6 +250,7 @@ ElectronMVANtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& i auto src = src_.getValidHandle(iEvent); auto vertices = vertices_.getValidHandle(iEvent); + // initialize cluster tools std::unique_ptr lazyTools; if(doEnergyMatrix_) { // Configure Lazy Tools, which will compute 5x5 quantities @@ -299,21 +298,18 @@ ElectronMVANtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& i eleIndex_ = src->size(); for(auto const& ele : src->ptrs()) { + if (ele->pt() < ptThreshold_) continue; + // Fill the energy matrix around the seed if(doEnergyMatrix_) { const auto& seed = *(ele->superCluster()->seed()); - energyMatrix_.clear(); - for(float e : lazyTools->getEnergies(seed, caloRectangle_)) energyMatrix_.push_back(e); + energyMatrix_ = lazyTools->getEnergies(seed, caloRectangle_); } // Fill various tree variable eleQ_ = ele->charge(); ele3Q_ = ele->chargeInfo().isGsfCtfScPixConsistent; - if (ele->pt() < ptThreshold_) { - continue; - } - for (int iVar = 0; iVar < nVars_; ++iVar) { std::vector extraVariables = variableHelper_.getAuxVariables(ele, iEvent); vars_[iVar] = mvaVarMngr_.getValue(iVar, *ele, extraVariables); diff --git a/RecoEgamma/ElectronIdentification/test/testElectronMVA_cfg.py b/RecoEgamma/ElectronIdentification/test/testElectronMVA_cfg.py index bc920c2056c51..665bdbe491977 100644 --- a/RecoEgamma/ElectronIdentification/test/testElectronMVA_cfg.py +++ b/RecoEgamma/ElectronIdentification/test/testElectronMVA_cfg.py @@ -15,10 +15,10 @@ outputFile = "electron_ntuple.root" -process.maxEvents = cms.PSet( input = cms.int32(-1) ) +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) process.source = cms.Source("PoolSource", - fileNames = cms.vstring( + fileNames = cms.untracked.vstring( '/store/mc/RunIIFall17MiniAOD/DYJetsToLL_M-50_TuneCP5_13TeV-madgraphMLM-pythia8/MINIAODSIM/RECOSIMstep_94X_mc2017_realistic_v10-v1/00000/0293A280-B5F3-E711-8303-3417EBE33927.root' ) ) @@ -122,35 +122,33 @@ ptThreshold = cms.double(5.0), # doEnergyMatrix = cms.bool(False), # disabled by default due to large size - energyMatrixSize = cms.int(2), # corresponding to 5x5 - """ - The energy matrix is for ecal driven electrons the n x n of raw - rec-hit energies around the seed crystal. - - The size of the electron matrix is controlled with the parameter - "energyMatrixSize", which controlls the extension of crystals in each - direction away from the seed, in other words n = 2 * energyMatrixSize + 1. - - The energy matrix gets saved as a vector but you can easily unroll it - to a two dimensional numpy array later, for example like that: - - >>> import uproot - >>> import numpy as np - >>> import matplotlib.pyplot as plt - - >>> tree = uproot.open("electron_ntuple.root")["ntuplizer/tree"] - >>> n = 5 - - >>> for a in tree.array("ele_energyMatrix"): - >>> a = a.reshape((n,n)) - >>> plt.imshow(np.log10(a)) - >>> plt.colorbar() - >>> plt.show() - """ + energyMatrixSize = cms.int32(2) # corresponding to 5x5 ) +""" +The energy matrix is for ecal driven electrons the n x n of raw +rec-hit energies around the seed crystal. -process.TFileService = cms.Service("TFileService", - fileName = cms.string( outputFile ) - ) +The size of the energy matrix is controlled with the parameter +"energyMatrixSize", which controlls the extension of crystals in each +direction away from the seed, in other words n = 2 * energyMatrixSize + 1. + +The energy matrix gets saved as a vector but you can easily unroll it +to a two dimensional numpy array later, for example like that: + +>>> import uproot +>>> import numpy as np +>>> import matplotlib.pyplot as plt + +>>> tree = uproot.open("electron_ntuple.root")["ntuplizer/tree"] +>>> n = 5 + +>>> for a in tree.array("ele_energyMatrix"): +>>> a = a.reshape((n,n)) +>>> plt.imshow(np.log10(a)) +>>> plt.colorbar() +>>> plt.show() +""" + +process.TFileService = cms.Service("TFileService", fileName = cms.string(outputFile)) process.p = cms.Path(process.egmGsfElectronIDSequence * process.ntuplizer) diff --git a/RecoEgamma/PhotonIdentification/plugins/PhotonMVANtuplizer.cc b/RecoEgamma/PhotonIdentification/plugins/PhotonMVANtuplizer.cc index 542c87ef65dd6..f418cae8e0e07 100644 --- a/RecoEgamma/PhotonIdentification/plugins/PhotonMVANtuplizer.cc +++ b/RecoEgamma/PhotonIdentification/plugins/PhotonMVANtuplizer.cc @@ -17,28 +17,22 @@ // -// user include files #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/one/EDAnalyzer.h" - #include "FWCore/Framework/interface/Event.h" - #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/Utilities/interface/InputTag.h" - #include "FWCore/ServiceRegistry/interface/Service.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" - #include "DataFormats/EgammaCandidates/interface/Photon.h" #include "DataFormats/PatCandidates/interface/Photon.h" - #include "RecoEgamma/EgammaTools/interface/MVAVariableManager.h" #include "RecoEgamma/EgammaTools/interface/MultiToken.h" - #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" #include "DataFormats/VertexReco/interface/Vertex.h" #include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h" #include #include @@ -64,11 +58,14 @@ class PhotonMVANtuplizer : public edm::one::EDAnalyzer energyMatrix_; // photon genMatch variable int matchedToGenPh_; @@ -101,6 +98,8 @@ class PhotonMVANtuplizer : public edm::one::EDAnalyzer> vertices_; const MultiTokenT> pileup_; const MultiTokenT> genParticles_; + const MultiTokenT ebRecHits_; + const MultiTokenT eeRecHits_; // to hold ID decisions and categories std::vector mvaPasses_; @@ -115,6 +114,9 @@ class PhotonMVANtuplizer : public edm::one::EDAnalyzer vars_; + + const bool doEnergyMatrix_; + const CaloRectangle caloRectangle_; }; enum PhotonMatchType { @@ -155,14 +157,14 @@ namespace { // constructor PhotonMVANtuplizer::PhotonMVANtuplizer(const edm::ParameterSet& iConfig) - : phoMapTags_ (iConfig.getUntrackedParameter>("phoMVAs")) - , phoMapBranchNames_ (iConfig.getUntrackedParameter>("phoMVALabels")) + : phoMapTags_ (iConfig.getParameter>("phoMVAs")) + , phoMapBranchNames_ (iConfig.getParameter>("phoMVALabels")) , nPhoMaps_ (phoMapBranchNames_.size()) - , valMapTags_ (iConfig.getUntrackedParameter>("phoMVAValMaps")) - , valMapBranchNames_ (iConfig.getUntrackedParameter>("phoMVAValMapLabels")) + , valMapTags_ (iConfig.getParameter>("phoMVAValMaps")) + , valMapBranchNames_ (iConfig.getParameter>("phoMVAValMapLabels")) , nValMaps_ (valMapBranchNames_.size()) - , mvaCatTags_ (iConfig.getUntrackedParameter>("phoMVACats")) - , mvaCatBranchNames_ (iConfig.getUntrackedParameter>("phoMVACatLabels")) + , mvaCatTags_ (iConfig.getParameter>("phoMVACats")) + , mvaCatBranchNames_ (iConfig.getParameter>("phoMVACatLabels")) , nCats_ (mvaCatBranchNames_.size()) , isMC_ (iConfig.getParameter("isMC")) , ptThreshold_ (iConfig.getParameter("ptThreshold")) @@ -171,6 +173,8 @@ PhotonMVANtuplizer::PhotonMVANtuplizer(const edm::ParameterSet& iConfig) , vertices_ (src_, consumesCollector(), iConfig, "vertices", "verticesMiniAOD") , pileup_ (src_, consumesCollector(), iConfig, "pileup" , "pileupMiniAOD") , genParticles_ (src_, consumesCollector(), iConfig, "genParticles", "genParticlesMiniAOD") + , ebRecHits_ (src_, consumesCollector(), iConfig, "ebReducedRecHitCollection", "ebReducedRecHitCollectionMiniAOD") + , eeRecHits_ (src_, consumesCollector(), iConfig, "eeReducedRecHitCollection", "eeReducedRecHitCollectionMiniAOD") , mvaPasses_ (nPhoMaps_) , mvaValues_ (nValMaps_) , mvaCats_ (nCats_) @@ -178,6 +182,8 @@ PhotonMVANtuplizer::PhotonMVANtuplizer(const edm::ParameterSet& iConfig) , mvaVarMngr_ (iConfig.getParameter("variableDefinition")) , nVars_ (mvaVarMngr_.getNVars()) , vars_ (nVars_) + , doEnergyMatrix_ (iConfig.getParameter("doEnergyMatrix")) + , caloRectangle_ (makeCaloRectangle(iConfig.getParameter("energyMatrixSize"))) { // phoMaps for (auto const& tag : phoMapTags_) { @@ -208,6 +214,8 @@ PhotonMVANtuplizer::PhotonMVANtuplizer(const edm::ParameterSet& iConfig) tree_->Branch("pT", &pT_); tree_->Branch("eta", &eta_); + if (doEnergyMatrix_) tree_->Branch("energyMatrix",&energyMatrix_); + for (int i = 0; i < nVars_; ++i) { tree_->Branch(mvaVarMngr_.getName(i).c_str(), &vars_[i]); } @@ -243,6 +251,14 @@ PhotonMVANtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe vtxN_ = vertices->size(); + // initialize cluster tools + std::unique_ptr lazyTools; + if(doEnergyMatrix_) { + // Configure Lazy Tools, which will compute 5x5 quantities + lazyTools = std::make_unique( + iEvent, iSetup, ebRecHits_.get(iEvent), eeRecHits_.get(iEvent)); + } + // Fill with true number of pileup if(isMC_) { for(const auto& pu : *pileup) @@ -276,36 +292,31 @@ PhotonMVANtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe for(auto const& pho : src->ptrs()) { - if (pho->pt() < ptThreshold_) { - continue; - } + if (pho->pt() < ptThreshold_) continue; + pT_ = pho->pt(); eta_ = pho->eta(); + // Fill the energy matrix around the seed + if(doEnergyMatrix_) { + const auto& seed = *(pho->superCluster()->seed()); + energyMatrix_ = lazyTools->getEnergies(seed, caloRectangle_); + } + // variables from the text file for (int iVar = 0; iVar < nVars_; ++iVar) { std::vector extraVariables = variableHelper_.getAuxVariables(pho, iEvent); vars_[iVar] = mvaVarMngr_.getValue(iVar, *pho, extraVariables); } - if (isMC_) { - matchedToGenPh_ = matchToTruth( *pho, *genParticles, deltaR_); - } + if (isMC_) matchedToGenPh_ = matchToTruth( *pho, *genParticles, deltaR_); // // Look up and save the ID decisions // - for (size_t k = 0; k < nPhoMaps_; ++k) { - mvaPasses_[k] = static_cast((*decisions[k])[pho]); - } - - for (size_t k = 0; k < nValMaps_; ++k) { - mvaValues_[k] = (*values[k])[pho]; - } - - for (size_t k = 0; k < nCats_; ++k) { - mvaCats_[k] = (*mvaCats[k])[pho]; - } + for (size_t k = 0; k < nPhoMaps_; ++k) mvaPasses_[k] = static_cast((*decisions[k])[pho]); + for (size_t k = 0; k < nValMaps_; ++k) mvaValues_[k] = (*values[k])[pho]; + for (size_t k = 0; k < nCats_ ; ++k) mvaCats_[k] = (*mvaCats[k])[pho]; tree_->Fill(); } @@ -325,12 +336,18 @@ PhotonMVANtuplizer::fillDescriptions(edm::ConfigurationDescriptions& description desc.add("verticesMiniAOD", edm::InputTag("offlineSlimmedPrimaryVertices")); desc.add("pileupMiniAOD", edm::InputTag("slimmedAddPileupInfo")); desc.add("genParticlesMiniAOD", edm::InputTag("prunedGenParticles")); - desc.addUntracked>("phoMVAs", {}); - desc.addUntracked>("phoMVALabels", {}); - desc.addUntracked>("phoMVAValMaps", {}); - desc.addUntracked>("phoMVAValMapLabels", {}); - desc.addUntracked>("phoMVACats", {}); - desc.addUntracked>("phoMVACatLabels", {}); + desc.add("ebReducedRecHitCollection", edm::InputTag("reducedEcalRecHitsEB")); + desc.add("eeReducedRecHitCollection", edm::InputTag("reducedEcalRecHitsEE")); + desc.add("ebReducedRecHitCollectionMiniAOD", edm::InputTag("reducedEgamma","reducedEBRecHits")); + desc.add("eeReducedRecHitCollectionMiniAOD", edm::InputTag("reducedEgamma","reducedEERecHits")); + desc.add>("phoMVAs", {}); + desc.add>("phoMVALabels", {}); + desc.add>("phoMVAValMaps", {}); + desc.add>("phoMVAValMapLabels", {}); + desc.add>("phoMVACats", {}); + desc.add>("phoMVACatLabels", {}); + desc.add("doEnergyMatrix", false); + desc.add("energyMatrixSize", 2)->setComment("extension of crystals in each direction away from the seed"); desc.add("isMC", true); desc.add("ptThreshold", 15.0); desc.add("deltaR", 0.1); diff --git a/RecoEgamma/PhotonIdentification/test/testPhotonMVA_cfg.py b/RecoEgamma/PhotonIdentification/test/testPhotonMVA_cfg.py index dafe8e2ebd74d..4539215b8624b 100644 --- a/RecoEgamma/PhotonIdentification/test/testPhotonMVA_cfg.py +++ b/RecoEgamma/PhotonIdentification/test/testPhotonMVA_cfg.py @@ -48,33 +48,58 @@ setupAllVIDIdsInModule(process,idmod,setupVIDPhotonSelection) process.ntuplizer = cms.EDAnalyzer('PhotonMVANtuplizer', - phoMVAs = cms.untracked.vstring( + phoMVAs = cms.vstring( ), - phoMVALabels = cms.untracked.vstring( + phoMVALabels = cms.vstring( ), - phoMVAValMaps = cms.untracked.vstring( + phoMVAValMaps = cms.vstring( "photonMVAValueMapProducer:PhotonMVAEstimatorRun2Spring16NonTrigV1Values", "photonMVAValueMapProducer:PhotonMVAEstimatorRunIIFall17v1Values", "photonMVAValueMapProducer:PhotonMVAEstimatorRunIIFall17v1p1Values", "photonMVAValueMapProducer:PhotonMVAEstimatorRunIIFall17v2Values", ), - phoMVAValMapLabels = cms.untracked.vstring( + phoMVAValMapLabels = cms.vstring( "Spring16NonTrigV1", "Fall17v1", "Fall17v1p1", "Fall17v2", ), - phoMVACats = cms.untracked.vstring( + phoMVACats = cms.vstring( "photonMVAValueMapProducer:PhotonMVAEstimatorRunIIFall17v1Categories", ), - phoMVACatLabels = cms.untracked.vstring( + phoMVACatLabels = cms.vstring( "PhoMVACats", ), variableDefinition = cms.string(mvaVariablesFile), + # + doEnergyMatrix = cms.bool(False), # disabled by default due to large size + energyMatrixSize = cms.int32(2) # corresponding to 5x5 ) +""" +The energy matrix is the n x n of raw rec-hit energies around the seed +crystal. -process.TFileService = cms.Service("TFileService", - fileName = cms.string( outputFile ) - ) +The size of the energy matrix is controlled with the parameter +"energyMatrixSize", which controlls the extension of crystals in each +direction away from the seed, in other words n = 2 * energyMatrixSize + 1. + +The energy matrix gets saved as a vector but you can easily unroll it +to a two dimensional numpy array later, for example like that: + +>>> import uproot +>>> import numpy as np +>>> import matplotlib.pyplot as plt + +>>> tree = uproot.open("photon_ntuple.root")["ntuplizer/tree"] +>>> n = 5 + +>>> for a in tree.array("ele_energyMatrix"): +>>> a = a.reshape((n,n)) +>>> plt.imshow(np.log10(a)) +>>> plt.colorbar() +>>> plt.show() +""" + +process.TFileService = cms.Service("TFileService", fileName = cms.string(outputFile)) process.p = cms.Path(process.egmPhotonIDSequence * process.ntuplizer)