Skip to content

Commit

Permalink
Energy matrix with PhotonMVA ntuplizer
Browse files Browse the repository at this point in the history
  • Loading branch information
@guitargeek
guitargeek committed Jan 12, 2019
1 parent 061cf33 commit a1a713e
Show file tree
Hide file tree
Showing 5 changed files with 144 additions and 102 deletions.
34 changes: 20 additions & 14 deletions RecoCaloTools/Navigation/interface/CaloRectangle.h
View file
Expand Up @@ -4,20 +4,11 @@
#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"

template<class T>
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;
Expand All @@ -38,6 +29,21 @@ inline auto makeCaloRectangle(int sizeX, int sizeY) {
return CaloRectangle{-sizeX, sizeX, -sizeY, sizeY};
}

template<class T>
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 T>
class CaloRectangleRange {

Expand Down
20 changes: 8 additions & 12 deletions RecoEgamma/ElectronIdentification/plugins/ElectronMVANtuplizer.cc
View file
Expand Up @@ -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_)
Expand Down Expand Up @@ -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_);
Expand Down Expand Up @@ -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<noZS::EcalClusterLazyTools> lazyTools;
if(doEnergyMatrix_) {
// Configure Lazy Tools, which will compute 5x5 quantities
Expand Down Expand Up @@ -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<float> extraVariables = variableHelper_.getAuxVariables(ele, iEvent);
vars_[iVar] = mvaVarMngr_.getValue(iVar, *ele, extraVariables);
Expand Down
58 changes: 28 additions & 30 deletions RecoEgamma/ElectronIdentification/test/testElectronMVA_cfg.py
View file
Expand Up @@ -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'
)
)
Expand Down Expand Up @@ -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)
91 changes: 54 additions & 37 deletions RecoEgamma/PhotonIdentification/plugins/PhotonMVANtuplizer.cc
View file
Expand Up @@ -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 <TTree.h>
#include <TFile.h>
Expand All @@ -64,11 +58,14 @@ class PhotonMVANtuplizer : public edm::one::EDAnalyzer<edm::one::SharedResources
// other
TTree* tree_;

//global variables
// global variables
int nEvent_, nRun_, nLumi_;
int genNpu_;
int vtxN_;

// photon variables
double pT_, eta_;
std::vector<float> energyMatrix_;

// photon genMatch variable
int matchedToGenPh_;
Expand Down Expand Up @@ -101,6 +98,8 @@ class PhotonMVANtuplizer : public edm::one::EDAnalyzer<edm::one::SharedResources
const MultiTokenT<std::vector<reco::Vertex>> vertices_;
const MultiTokenT<std::vector<PileupSummaryInfo>> pileup_;
const MultiTokenT<edm::View<reco::GenParticle>> genParticles_;
const MultiTokenT<EcalRecHitCollection> ebRecHits_;
const MultiTokenT<EcalRecHitCollection> eeRecHits_;

// to hold ID decisions and categories
std::vector<int> mvaPasses_;
Expand All @@ -115,6 +114,9 @@ class PhotonMVANtuplizer : public edm::one::EDAnalyzer<edm::one::SharedResources

const int nVars_;
std::vector<float> vars_;

const bool doEnergyMatrix_;
const CaloRectangle caloRectangle_;
};

enum PhotonMatchType {
Expand Down Expand Up @@ -155,14 +157,14 @@ namespace {

// constructor
PhotonMVANtuplizer::PhotonMVANtuplizer(const edm::ParameterSet& iConfig)
: phoMapTags_ (iConfig.getUntrackedParameter<std::vector<std::string>>("phoMVAs"))
, phoMapBranchNames_ (iConfig.getUntrackedParameter<std::vector<std::string>>("phoMVALabels"))
: phoMapTags_ (iConfig.getParameter<std::vector<std::string>>("phoMVAs"))
, phoMapBranchNames_ (iConfig.getParameter<std::vector<std::string>>("phoMVALabels"))
, nPhoMaps_ (phoMapBranchNames_.size())
, valMapTags_ (iConfig.getUntrackedParameter<std::vector<std::string>>("phoMVAValMaps"))
, valMapBranchNames_ (iConfig.getUntrackedParameter<std::vector<std::string>>("phoMVAValMapLabels"))
, valMapTags_ (iConfig.getParameter<std::vector<std::string>>("phoMVAValMaps"))
, valMapBranchNames_ (iConfig.getParameter<std::vector<std::string>>("phoMVAValMapLabels"))
, nValMaps_ (valMapBranchNames_.size())
, mvaCatTags_ (iConfig.getUntrackedParameter<std::vector<std::string>>("phoMVACats"))
, mvaCatBranchNames_ (iConfig.getUntrackedParameter<std::vector<std::string>>("phoMVACatLabels"))
, mvaCatTags_ (iConfig.getParameter<std::vector<std::string>>("phoMVACats"))
, mvaCatBranchNames_ (iConfig.getParameter<std::vector<std::string>>("phoMVACatLabels"))
, nCats_ (mvaCatBranchNames_.size())
, isMC_ (iConfig.getParameter<bool>("isMC"))
, ptThreshold_ (iConfig.getParameter<double>("ptThreshold"))
Expand All @@ -171,13 +173,17 @@ 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_)
, variableHelper_ (consumesCollector())
, mvaVarMngr_ (iConfig.getParameter<std::string>("variableDefinition"))
, nVars_ (mvaVarMngr_.getNVars())
, vars_ (nVars_)
, doEnergyMatrix_ (iConfig.getParameter<bool>("doEnergyMatrix"))
, caloRectangle_ (makeCaloRectangle(iConfig.getParameter<int>("energyMatrixSize")))
{
// phoMaps
for (auto const& tag : phoMapTags_) {
Expand Down Expand Up @@ -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]);
}
Expand Down Expand Up @@ -243,6 +251,14 @@ PhotonMVANtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe

vtxN_ = vertices->size();

// initialize cluster tools
std::unique_ptr<noZS::EcalClusterLazyTools> lazyTools;
if(doEnergyMatrix_) {
// Configure Lazy Tools, which will compute 5x5 quantities
lazyTools = std::make_unique<noZS::EcalClusterLazyTools>(
iEvent, iSetup, ebRecHits_.get(iEvent), eeRecHits_.get(iEvent));
}

// Fill with true number of pileup
if(isMC_) {
for(const auto& pu : *pileup)
Expand Down Expand Up @@ -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<float> 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<int>((*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<int>((*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();
}
Expand All @@ -325,12 +336,18 @@ PhotonMVANtuplizer::fillDescriptions(edm::ConfigurationDescriptions& description
desc.add<edm::InputTag>("verticesMiniAOD", edm::InputTag("offlineSlimmedPrimaryVertices"));
desc.add<edm::InputTag>("pileupMiniAOD", edm::InputTag("slimmedAddPileupInfo"));
desc.add<edm::InputTag>("genParticlesMiniAOD", edm::InputTag("prunedGenParticles"));
desc.addUntracked<std::vector<std::string>>("phoMVAs", {});
desc.addUntracked<std::vector<std::string>>("phoMVALabels", {});
desc.addUntracked<std::vector<std::string>>("phoMVAValMaps", {});
desc.addUntracked<std::vector<std::string>>("phoMVAValMapLabels", {});
desc.addUntracked<std::vector<std::string>>("phoMVACats", {});
desc.addUntracked<std::vector<std::string>>("phoMVACatLabels", {});
desc.add<edm::InputTag>("ebReducedRecHitCollection", edm::InputTag("reducedEcalRecHitsEB"));
desc.add<edm::InputTag>("eeReducedRecHitCollection", edm::InputTag("reducedEcalRecHitsEE"));
desc.add<edm::InputTag>("ebReducedRecHitCollectionMiniAOD", edm::InputTag("reducedEgamma","reducedEBRecHits"));
desc.add<edm::InputTag>("eeReducedRecHitCollectionMiniAOD", edm::InputTag("reducedEgamma","reducedEERecHits"));
desc.add<std::vector<std::string>>("phoMVAs", {});
desc.add<std::vector<std::string>>("phoMVALabels", {});
desc.add<std::vector<std::string>>("phoMVAValMaps", {});
desc.add<std::vector<std::string>>("phoMVAValMapLabels", {});
desc.add<std::vector<std::string>>("phoMVACats", {});
desc.add<std::vector<std::string>>("phoMVACatLabels", {});
desc.add<bool>("doEnergyMatrix", false);
desc.add<int>("energyMatrixSize", 2)->setComment("extension of crystals in each direction away from the seed");
desc.add<bool>("isMC", true);
desc.add<double>("ptThreshold", 15.0);
desc.add<double>("deltaR", 0.1);
Expand Down

0 comments on commit a1a713e

Please sign in to comment.