PhotonIdentification/plugins/PhotonMVANtuplizer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoEgamma/PhotonIdentification
0004 // Class:      PhotonMVANtuplizer
0005 //
0006 /**\class PhotonMVANtuplizer PhotonMVANtuplizer.cc RecoEgamma/PhotonIdentification/plugins/PhotonMVANtuplizer.cc
0007
0008  Description: Ntuplizer to use for testing photon MVA IDs.
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Jonas REMBSER
0015 //         Created:  Thu, 22 Mar 2018 14:54:24 GMT
0016 //
0017 //
0018
0019 #include "FWCore/Framework/interface/Frameworkfwd.h"
0020 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024 #include "FWCore/ServiceRegistry/interface/Service.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0027 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0028 #include "DataFormats/PatCandidates/interface/Photon.h"
0029 #include "RecoEgamma/EgammaTools/interface/MVAVariableHelper.h"
0030 #include "RecoEgamma/EgammaTools/interface/MVAVariableManager.h"
0031 #include "FWCore/Utilities/interface/EDGetToken.h"
0032 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0033 #include "DataFormats/VertexReco/interface/Vertex.h"
0034 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0035 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
0036
0037 #include <TTree.h>
0038 #include <TFile.h>
0039 #include <Math/VectorUtil.h>
0040
0041 //
0042 // class declaration
0043 //
0044
0045 class PhotonMVANtuplizer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0046 public:
0047   explicit PhotonMVANtuplizer(const edm::ParameterSet&);
0048
0049   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 private:
0052   void analyze(const edm::Event&, const edm::EventSetup&) override;
0053
0054   // ----------member data ---------------------------
0055
0056   // other
0057   TTree* tree_;
0058
0059   // global variables
0060   int nEvent_, nRun_, nLumi_;
0061   int genNpu_;
0062   int vtxN_;
0063
0064   // photon variables
0065   double pT_, eta_;
0066   std::vector<float> energyMatrix_;
0067
0068   // photon genMatch variable
0069   int matchedToGenPh_;
0070   int matchedGenIdx_;
0071
0072   // ID decisions objects
0073   const std::vector<std::string> phoMapTags_;
0074   std::vector<edm::EDGetTokenT<edm::ValueMap<bool>>> phoMapTokens_;
0075   const std::vector<std::string> phoMapBranchNames_;
0076   const size_t nPhoMaps_;
0077
0078   // MVA values and categories (optional)
0079   const std::vector<std::string> valMapTags_;
0080   std::vector<edm::EDGetTokenT<edm::ValueMap<float>>> valMapTokens_;
0081   const std::vector<std::string> valMapBranchNames_;
0082   const size_t nValMaps_;
0083
0084   const std::vector<std::string> mvaCatTags_;
0085   std::vector<edm::EDGetTokenT<edm::ValueMap<int>>> mvaCatTokens_;
0086   const std::vector<std::string> mvaCatBranchNames_;
0087   const size_t nCats_;
0088
0089   // config
0090   const bool isMC_;
0091   const double ptThreshold_;
0092   const double deltaR_;
0093
0094   // Tokens
0095   const edm::EDGetTokenT<edm::View<reco::Photon>> src_;
0096   const edm::EDGetTokenT<std::vector<reco::Vertex>> vertices_;
0097   const edm::EDGetTokenT<std::vector<PileupSummaryInfo>> pileup_;
0098   const edm::EDGetTokenT<edm::View<reco::GenParticle>> genParticles_;
0099   const edm::EDGetTokenT<EcalRecHitCollection> ebRecHits_;
0100   const edm::EDGetTokenT<EcalRecHitCollection> eeRecHits_;
0101
0102   const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_;
0103
0104   // to hold ID decisions and categories
0105   std::vector<int> mvaPasses_;
0106   std::vector<float> mvaValues_;
0107   std::vector<int> mvaCats_;
0108
0109   // To get the auxiliary MVA variables
0110   const MVAVariableHelper variableHelper_;
0111
0112   // To manage the variables which are parsed from the text file
0113   MVAVariableManager<reco::Photon> mvaVarMngr_;
0114
0115   const int nVars_;
0116   std::vector<float> vars_;
0117
0118   const bool doEnergyMatrix_;
0119   const int energyMatrixSize_;
0120 };
0121
0122 enum PhotonMatchType {
0123   FAKE_PHOTON,
0124   TRUE_PROMPT_PHOTON,
0125   TRUE_NON_PROMPT_PHOTON,
0126 };
0127
0128 namespace {
0129
0130   int matchToTruth(const reco::Photon& ph, const edm::View<reco::GenParticle>& genParticles, double deltaR) {
0131     // Find the closest status 1 gen photon to the reco photon
0132     double dR = 999;
0133     reco::GenParticle const* closestPhoton = &genParticles[0];
0134     for (auto& particle : genParticles) {
0135       // Drop everything that is not photon or not status 1
0136       if (abs(particle.pdgId()) != 22 || particle.status() != 1)
0137         continue;
0138
0139       double dRtmp = ROOT::Math::VectorUtil::DeltaR(ph.p4(), particle.p4());
0140       if (dRtmp < dR) {
0141         dR = dRtmp;
0142         closestPhoton = &particle;
0143       }
0144     }
0145     // See if the closest photon (if it exists) is close enough.
0146     // If not, no match found.
0147     if (dR < deltaR) {
0148       if (closestPhoton->isPromptFinalState())
0149         return TRUE_PROMPT_PHOTON;
0150       else
0151         return TRUE_NON_PROMPT_PHOTON;
0152     }
0153     return FAKE_PHOTON;
0154   }
0155
0156 };  // namespace
0157
0158 // constructor
0159 PhotonMVANtuplizer::PhotonMVANtuplizer(const edm::ParameterSet& iConfig)
0160     : phoMapTags_(iConfig.getParameter<std::vector<std::string>>("phoMVAs")),
0161       phoMapBranchNames_(iConfig.getParameter<std::vector<std::string>>("phoMVALabels")),
0162       nPhoMaps_(phoMapBranchNames_.size()),
0163       valMapTags_(iConfig.getParameter<std::vector<std::string>>("phoMVAValMaps")),
0164       valMapBranchNames_(iConfig.getParameter<std::vector<std::string>>("phoMVAValMapLabels")),
0165       nValMaps_(valMapBranchNames_.size()),
0166       mvaCatTags_(iConfig.getParameter<std::vector<std::string>>("phoMVACats")),
0167       mvaCatBranchNames_(iConfig.getParameter<std::vector<std::string>>("phoMVACatLabels")),
0168       nCats_(mvaCatBranchNames_.size()),
0169       isMC_(iConfig.getParameter<bool>("isMC")),
0170       ptThreshold_(iConfig.getParameter<double>("ptThreshold")),
0171       deltaR_(iConfig.getParameter<double>("deltaR")),
0172       src_(consumes<edm::View<reco::Photon>>(iConfig.getParameter<edm::InputTag>("src"))),
0173       vertices_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("vertices"))),
0174       pileup_(consumes<std::vector<PileupSummaryInfo>>(iConfig.getParameter<edm::InputTag>("pileup"))),
0175       genParticles_(consumes<edm::View<reco::GenParticle>>(iConfig.getParameter<edm::InputTag>("genParticles"))),
0176       ebRecHits_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ebReducedRecHitCollection"))),
0177       eeRecHits_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("eeReducedRecHitCollection"))),
0178       ecalClusterToolsESGetTokens_{consumesCollector()},
0179       mvaPasses_(nPhoMaps_),
0180       mvaValues_(nValMaps_),
0181       mvaCats_(nCats_),
0182       variableHelper_(consumesCollector()),
0183       mvaVarMngr_(iConfig.getParameter<std::string>("variableDefinition"), MVAVariableHelper::indexMap()),
0184       nVars_(mvaVarMngr_.getNVars()),
0185       vars_(nVars_),
0186       doEnergyMatrix_(iConfig.getParameter<bool>("doEnergyMatrix")),
0187       energyMatrixSize_(iConfig.getParameter<int>("energyMatrixSize")) {
0188   // phoMaps
0189   for (auto const& tag : phoMapTags_) {
0190     phoMapTokens_.push_back(consumes<edm::ValueMap<bool>>(edm::InputTag(tag)));
0191   }
0192   // valMaps
0193   for (auto const& tag : valMapTags_) {
0194     valMapTokens_.push_back(consumes<edm::ValueMap<float>>(edm::InputTag(tag)));
0195   }
0196   // categories
0197   for (auto const& tag : mvaCatTags_) {
0198     mvaCatTokens_.push_back(consumes<edm::ValueMap<int>>(edm::InputTag(tag)));
0199   }
0200
0201   // Book tree
0202   usesResource(TFileService::kSharedResource);
0203   edm::Service<TFileService> fs;
0204   tree_ = fs->make<TTree>("tree", "tree");
0205
0206   tree_->Branch("nEvent", &nEvent_);
0207   tree_->Branch("nRun", &nRun_);
0208   tree_->Branch("nLumi", &nLumi_);
0209   if (isMC_) {
0210     tree_->Branch("genNpu", &genNpu_);
0211     tree_->Branch("matchedToGenPh", &matchedToGenPh_);
0212   }
0213   tree_->Branch("vtxN", &vtxN_);
0214   tree_->Branch("pT", &pT_);
0215   tree_->Branch("eta", &eta_);
0216
0217   if (doEnergyMatrix_)
0218     tree_->Branch("energyMatrix", &energyMatrix_);
0219
0220   for (int i = 0; i < nVars_; ++i) {
0221     tree_->Branch(mvaVarMngr_.getName(i).c_str(), &vars_[i]);
0222   }
0223
0224   // IDs
0225   for (size_t k = 0; k < nValMaps_; ++k) {
0226     tree_->Branch(valMapBranchNames_[k].c_str(), &mvaValues_[k]);
0227   }
0228
0229   for (size_t k = 0; k < nPhoMaps_; ++k) {
0230     tree_->Branch(phoMapBranchNames_[k].c_str(), &mvaPasses_[k]);
0231   }
0232
0233   for (size_t k = 0; k < nCats_; ++k) {
0234     tree_->Branch(mvaCatBranchNames_[k].c_str(), &mvaCats_[k]);
0235   }
0236 }
0237
0238 // ------------ method called for each event  ------------
0239 void PhotonMVANtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0240   // Fill global event info
0241   nEvent_ = iEvent.id().event();
0242   nRun_ = iEvent.id().run();
0243   nLumi_ = iEvent.luminosityBlock();
0244
0245   // Get Handles
0246   auto src = iEvent.getHandle(src_);
0247   auto vertices = iEvent.getHandle(vertices_);
0248   auto pileup = iEvent.getHandle(pileup_);
0249   auto genParticles = iEvent.getHandle(genParticles_);
0250
0251   vtxN_ = vertices->size();
0252
0253   // initialize cluster tools
0254   std::unique_ptr<noZS::EcalClusterLazyTools> lazyTools;
0255   if (doEnergyMatrix_) {
0256     // Configure Lazy Tools, which will compute 5x5 quantities
0257     lazyTools = std::make_unique<noZS::EcalClusterLazyTools>(
0258         iEvent, ecalClusterToolsESGetTokens_.get(iSetup), ebRecHits_, eeRecHits_);
0259   }
0260
0261   // Fill with true number of pileup
0262   if (isMC_) {
0263     for (const auto& pu : *pileup) {
0264       int bx = pu.getBunchCrossing();
0265       if (bx == 0) {
0266         genNpu_ = pu.getPU_NumInteractions();
0267         break;
0268       }
0269     }
0270   }
0271
0272   // Get MVA decisions
0273   edm::Handle<edm::ValueMap<bool>> decisions[nPhoMaps_];
0274   for (size_t k = 0; k < nPhoMaps_; ++k) {
0275     iEvent.getByToken(phoMapTokens_[k], decisions[k]);
0276   }
0277
0278   // Get MVA values
0279   edm::Handle<edm::ValueMap<float>> values[nValMaps_];
0280   for (size_t k = 0; k < nValMaps_; ++k) {
0281     iEvent.getByToken(valMapTokens_[k], values[k]);
0282   }
0283
0284   // Get MVA categories
0285   edm::Handle<edm::ValueMap<int>> mvaCats[nCats_];
0286   for (size_t k = 0; k < nCats_; ++k) {
0287     iEvent.getByToken(mvaCatTokens_[k], mvaCats[k]);
0288   }
0289
0290   std::vector<float> extraVariables = variableHelper_.getAuxVariables(iEvent);
0291
0292   for (auto const& pho : src->ptrs()) {
0293     if (pho->pt() < ptThreshold_)
0294       continue;
0295
0296     pT_ = pho->pt();
0297     eta_ = pho->eta();
0298
0299     // Fill the energy matrix around the seed
0300     if (doEnergyMatrix_) {
0301       const auto& seed = *(pho->superCluster()->seed());
0302       energyMatrix_ = lazyTools->energyMatrix(seed, energyMatrixSize_);
0303     }
0304
0305     // variables from the text file
0306     for (int iVar = 0; iVar < nVars_; ++iVar) {
0307       vars_[iVar] = mvaVarMngr_.getValue(iVar, *pho, extraVariables);
0308     }
0309
0310     if (isMC_)
0311       matchedToGenPh_ = matchToTruth(*pho, *genParticles, deltaR_);
0312
0313     //
0314     // Look up and save the ID decisions
0315     //
0316     for (size_t k = 0; k < nPhoMaps_; ++k)
0317       mvaPasses_[k] = static_cast<int>((*decisions[k])[pho]);
0318     for (size_t k = 0; k < nValMaps_; ++k)
0319       mvaValues_[k] = (*values[k])[pho];
0320     for (size_t k = 0; k < nCats_; ++k)
0321       mvaCats_[k] = (*mvaCats[k])[pho];
0322
0323     tree_->Fill();
0324   }
0325 }
0326
0327 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0328 void PhotonMVANtuplizer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0329   edm::ParameterSetDescription desc;
0330   desc.add<edm::InputTag>("src", edm::InputTag("slimmedPhotons"));
0331   desc.add<edm::InputTag>("vertices", edm::InputTag("offlineSlimmedPrimaryVertices"));
0332   desc.add<edm::InputTag>("pileup", edm::InputTag("slimmedAddPileupInfo"));
0333   desc.add<edm::InputTag>("genParticles", edm::InputTag("prunedGenParticles"));
0334   desc.add<edm::InputTag>("ebReducedRecHitCollection", edm::InputTag("reducedEgamma", "reducedEBRecHits"));
0335   desc.add<edm::InputTag>("eeReducedRecHitCollection", edm::InputTag("reducedEgamma", "reducedEERecHits"));
0336   desc.add<std::vector<std::string>>("phoMVAs", {});
0337   desc.add<std::vector<std::string>>("phoMVALabels", {});
0338   desc.add<std::vector<std::string>>("phoMVAValMaps", {});
0339   desc.add<std::vector<std::string>>("phoMVAValMapLabels", {});
0340   desc.add<std::vector<std::string>>("phoMVACats", {});
0341   desc.add<std::vector<std::string>>("phoMVACatLabels", {});
0342   desc.add<bool>("doEnergyMatrix", false);
0343   desc.add<int>("energyMatrixSize", 2)->setComment("extension of crystals in each direction away from the seed");
0344   desc.add<bool>("isMC", true);
0345   desc.add<double>("ptThreshold", 15.0);
0346   desc.add<double>("deltaR", 0.1);
0347   desc.add<std::string>("variableDefinition");
0348   descriptions.addDefault(desc);
0349 }
0350
0351 //define this as a plug-in
0352 DEFINE_FWK_MODULE(PhotonMVANtuplizer);