ElectronTools/test/MiniAODElectronIDValidationAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    TestElectronID/MiniAODElectronIDValidationAnalyzer
0004 // Class:      MiniAODElectronIDValidationAnalyzer
0005 //
0006 /**\class MiniAODElectronIDValidationAnalyzer MiniAODElectronIDValidationAnalyzer.cc TestElectronID/MiniAODElectronIDValidationAnalyzer/plugins/MiniAODElectronIDValidationAnalyzer.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Ilya Kravchenko
0015 //         Created:  Thu, 14 Aug 2014 08:27:41 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030
0031 #include "DataFormats/Common/interface/ValueMap.h"
0032 #include "DataFormats/Common/interface/View.h"
0033
0034 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0035 #include "DataFormats/VertexReco/interface/Vertex.h"
0036 #include "DataFormats/PatCandidates/interface/Electron.h"
0037
0038 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0039 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0040 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0041 #include "DataFormats/PatCandidates/interface/Electron.h"
0042
0043 #include "DataFormats/Candidate/interface/Candidate.h"
0044 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0045
0046 #include "DataFormats/EgammaCandidates/interface/ConversionFwd.h"
0047 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0048 #include "CommonTools/Egamma/interface/ConversionTools.h"
0049
0050 #include "FWCore/ServiceRegistry/interface/Service.h"
0051 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0052
0053 #include "TTree.h"
0054 #include "Math/VectorUtil.h"
0055
0056 //
0057 // class declaration
0058 //
0059
0060 class MiniAODElectronIDValidationAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0061 public:
0062   explicit MiniAODElectronIDValidationAnalyzer(const edm::ParameterSet &);
0063   ~MiniAODElectronIDValidationAnalyzer() override;
0064
0065   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0066
0067   enum ElectronMatchType {
0068     UNMATCHED = 0,
0069     TRUE_PROMPT_ELECTRON,
0070     TRUE_ELECTRON_FROM_TAU,
0071     TRUE_NON_PROMPT_ELECTRON
0072   };  // The last does not include tau parents
0073
0074 private:
0075   void beginJob() override;
0076   void analyze(const edm::Event &, const edm::EventSetup &) override;
0077   void endJob() override;
0078
0079   //virtual void beginRun(edm::Run const&, edm::EventSetup const&) override;
0080   //virtual void endRun(edm::Run const&, edm::EventSetup const&) override;
0081   //virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
0082   //virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override;
0083
0084   int matchToTruth(const reco::GsfElectron &el, const edm::Handle<edm::View<reco::GenParticle>> &genParticles);
0085   void findFirstNonElectronMother(const reco::Candidate *particle, int &ancestorPID, int &ancestorStatus);
0086
0087   // ----------member data ---------------------------
0088   edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0089   edm::EDGetTokenT<edm::View<reco::GenParticle>> genToken_;
0090   edm::EDGetTokenT<reco::ConversionCollection> convToken_;
0091   edm::EDGetTokenT<reco::BeamSpot> beamToken_;
0092   edm::EDGetTokenT<edm::ValueMap<float>> full5x5SigmaIEtaIEtaMapToken_;
0093   edm::EDGetTokenT<edm::View<pat::Electron>> electronCollectionToken_;
0094   edm::InputTag electronIdTag_;
0095   edm::EDGetTokenT<edm::ValueMap<bool>> electronIdToken_;
0096
0097   TTree *electronTree_;
0098   Float_t pt_;
0099   Float_t etaSC_;
0100   // All ID variables
0101   Float_t dEtaIn_;
0102   Float_t dPhiIn_;
0103   Float_t hOverE_;
0104   Float_t sigmaIetaIeta_;
0105   Float_t full5x5_sigmaIetaIeta_;
0106   Float_t relIsoWithDBeta_;
0107   Float_t ooEmooP_;
0108   Float_t d0_;
0109   Float_t dz_;
0110   Int_t expectedMissingInnerHits_;
0111   Int_t passConversionVeto_;
0112
0113   Int_t isTrue_;
0114   Int_t isPass_;
0115 };
0116
0117 //
0118 // constants, enums and typedefs
0119 //
0120
0121 //
0122 // static data member definitions
0123 //
0124
0125 //
0126 // constructors and destructor
0127 //
0128 MiniAODElectronIDValidationAnalyzer::MiniAODElectronIDValidationAnalyzer(const edm::ParameterSet &iConfig)
0129     : vtxToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0130       genToken_(consumes<edm::View<reco::GenParticle>>(iConfig.getParameter<edm::InputTag>("genparticles"))),
0131       convToken_(consumes<reco::ConversionCollection>(iConfig.getParameter<edm::InputTag>("convcollection"))),
0132       beamToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamspot"))),
0133       full5x5SigmaIEtaIEtaMapToken_(
0134           consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("full5x5SigmaIEtaIEtaMap"))),
0135       electronCollectionToken_(consumes<edm::View<pat::Electron>>(iConfig.getParameter<edm::InputTag>("electrons"))),
0136       electronIdTag_(iConfig.getParameter<edm::InputTag>("electronIDs")),
0137       electronIdToken_(consumes<edm::ValueMap<bool>>(iConfig.getParameter<edm::InputTag>("electronIDs"))) {
0138   //now do what ever initialization is needed
0139   usesResource(TFileService::kSharedResource);
0140   edm::Service<TFileService> fs;
0141   electronTree_ = fs->make<TTree>("ElectronTree", "Electron data");
0142
0143   electronTree_->Branch("pt", &pt_, "pt/F");
0144   electronTree_->Branch("etaSC", &etaSC_, "etaSC/F");
0145
0146   electronTree_->Branch("dEtaIn", &dEtaIn_, "dEtaIn/F");
0147   electronTree_->Branch("dPhiIn", &dPhiIn_, "dPhiIn/F");
0148   electronTree_->Branch("hOverE", &hOverE_, "hOverE/F");
0149   electronTree_->Branch("sigmaIetaIeta", &sigmaIetaIeta_, "sigmaIetaIeta/F");
0150   electronTree_->Branch("full5x5_sigmaIetaIeta", &full5x5_sigmaIetaIeta_, "full5x5_sigmaIetaIeta/F");
0151   electronTree_->Branch("relIsoWithDBeta", &relIsoWithDBeta_, "relIsoWithDBeta/F");
0152   electronTree_->Branch("ooEmooP", &ooEmooP_, "ooEmooP/F");
0153   electronTree_->Branch("d0", &d0_, "d0/F");
0154   electronTree_->Branch("dz", &dz_, "dz/F");
0155   electronTree_->Branch("expectedMissingInnerHits", &expectedMissingInnerHits_, "expectedMissingInnerHits/I");
0156   electronTree_->Branch("passConversionVeto", &passConversionVeto_, "passConversionVeto/I");
0157
0158   electronTree_->Branch("isTrue", &isTrue_, "isTrue/I");
0159   electronTree_->Branch("isPass", &isPass_, "isPass/I");
0160 }
0161
0162 MiniAODElectronIDValidationAnalyzer::~MiniAODElectronIDValidationAnalyzer() {
0163   // do anything here that needs to be done at desctruction time
0164   // (e.g. close files, deallocate resources etc.)
0165 }
0166
0167 //
0168 // member functions
0169 //
0170
0171 // ------------ method called for each event  ------------
0172 void MiniAODElectronIDValidationAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0173   //edm::Handle<edm::ValueMap<float> > full5x5sieie;
0174   edm::Handle<edm::View<pat::Electron>> collection;
0175   //edm::Handle<edm::ValueMap<bool> > id_decisions;
0176   //iEvent.getByToken(full5x5SigmaIEtaIEtaMapToken_,full5x5sieie);
0177   iEvent.getByToken(electronCollectionToken_, collection);
0178   //iEvent.getByToken(electronIdToken_,id_decisions);
0179
0180   // Get PV
0181   edm::Handle<reco::VertexCollection> vertices;
0182   iEvent.getByToken(vtxToken_, vertices);
0183   if (vertices->empty())
0184     return;  // skip the event if no PV found
0185   const reco::Vertex &pv = vertices->front();
0186
0187   // Get GenParticles
0188   edm::Handle<edm::View<reco::GenParticle>> genParticles;
0189   iEvent.getByToken(genToken_, genParticles);
0190
0191   // Get stuff for conversions
0192   edm::Handle<reco::ConversionCollection> convs;
0193   edm::Handle<reco::BeamSpot> thebs;
0194   iEvent.getByToken(convToken_, convs);
0195   iEvent.getByToken(beamToken_, thebs);
0196
0197   for (size_t i = 0; i < collection->size(); ++i) {
0198     auto el = collection->refAt(i);
0199     pt_ = el->pt();
0200     etaSC_ = el->superCluster()->eta();
0201
0202     // ID and matching
0203     dEtaIn_ = el->deltaEtaSuperClusterTrackAtVtx();
0204     dPhiIn_ = el->deltaPhiSuperClusterTrackAtVtx();
0205     hOverE_ = el->hcalOverEcal();
0206     sigmaIetaIeta_ = el->sigmaIetaIeta();
0207     full5x5_sigmaIetaIeta_ = el->full5x5_sigmaIetaIeta();
0208     // |1/E-1/p| = |1/E - EoverPinner/E| is computed below
0209     // The if protects against ecalEnergy == inf or zero (always
0210     // the case for electrons below 5 GeV in miniAOD)
0211     if (el->ecalEnergy() == 0) {
0212       printf("Electron energy is zero!\n");
0213       ooEmooP_ = 1e30;
0214     } else if (!std::isfinite(el->ecalEnergy())) {
0215       printf("Electron energy is not finite!\n");
0216       ooEmooP_ = 1e30;
0217     } else {
0218       ooEmooP_ = fabs(1.0 / el->ecalEnergy() - el->eSuperClusterOverP() / el->ecalEnergy());
0219     }
0220
0221     // Isolation
0222     reco::GsfElectron::PflowIsolationVariables pfIso = el->pfIsolationVariables();
0223     // Compute isolation with delta beta correction for PU
0224     float absiso =
0225         pfIso.sumChargedHadronPt + std::max(0.0, pfIso.sumNeutralHadronEt + pfIso.sumPhotonEt - 0.5 * pfIso.sumPUPt);
0226     relIsoWithDBeta_ = absiso / pt_;
0227
0228     // Impact parameter
0229     d0_ = (-1) * el->gsfTrack()->dxy(pv.position());
0230     dz_ = el->gsfTrack()->dz(pv.position());
0231
0232     // Conversion rejection
0233     constexpr reco::HitPattern::HitCategory missingHitType = reco::HitPattern::MISSING_INNER_HITS;
0234     expectedMissingInnerHits_ = el->gsfTrack()->hitPattern().numberOfLostHits(missingHitType);
0235     passConversionVeto_ = false;
0236     if (thebs.isValid() && convs.isValid()) {
0237       passConversionVeto_ = !ConversionTools::hasMatchedConversion(*el, *convs, thebs->position());
0238     } else {
0239       printf("\n\nERROR!!! conversions not found!!!\n");
0240     }
0241
0242     isTrue_ = matchToTruth(*el, genParticles);
0243     isPass_ = (el->electronID(electronIdTag_.encode()) > 0.5);
0244
0245     electronTree_->Fill();
0246   }
0247 }
0248
0249 // ------------ method called once each job just before starting event loop  ------------
0250 void MiniAODElectronIDValidationAnalyzer::beginJob() {}
0251
0252 // ------------ method called once each job just after ending the event loop  ------------
0253 void MiniAODElectronIDValidationAnalyzer::endJob() {}
0254
0255 // ------------ method called when starting to processes a run  ------------
0256 /*
0257 void
0258 MiniAODElectronIDValidationAnalyzer::beginRun(edm::Run const&, edm::EventSetup const&)
0259 {
0260 }
0261 */
0262
0263 // ------------ method called when ending the processing of a run  ------------
0264 /*
0265 void
0266 MiniAODElectronIDValidationAnalyzer::endRun(edm::Run const&, edm::EventSetup const&)
0267 {
0268 }
0269 */
0270
0271 // ------------ method called when starting to processes a luminosity block  ------------
0272 /*
0273 void
0274 MiniAODElectronIDValidationAnalyzer::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0275 {
0276 }
0277 */
0278
0279 // ------------ method called when ending the processing of a luminosity block  ------------
0280 /*
0281 void
0282 MiniAODElectronIDValidationAnalyzer::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0283 {
0284 }
0285 */
0286
0287 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0288 void MiniAODElectronIDValidationAnalyzer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0289   //The following says we do not know what parameters are allowed so do no validation
0290   // Please change this to state exactly what you do use, even if it is no parameters
0291   edm::ParameterSetDescription desc;
0292   desc.setUnknown();
0293   descriptions.addDefault(desc);
0294 }
0295
0296 // The function that uses algorith from Josh Bendavid with
0297 // an explicit loop over gen particles.
0298 int MiniAODElectronIDValidationAnalyzer::matchToTruth(const reco::GsfElectron &el,
0299                                                       const edm::Handle<edm::View<reco::GenParticle>> &genParticles) {
0300   //
0301   // Explicit loop and geometric matching method (advised by Josh Bendavid)
0302   //
0303
0304   // Find the closest status 1 gen electron to the reco electron
0305   double dR = 999;
0306   const reco::Candidate *closestElectron = 0;
0307   for (size_t i = 0; i < genParticles->size(); i++) {
0308     const reco::Candidate *particle = &(*genParticles)[i];
0309     // Drop everything that is not electron or not status 1
0310     if (abs(particle->pdgId()) != 11 || particle->status() != 1)
0311       continue;
0312     //
0313     double dRtmp = ROOT::Math::VectorUtil::DeltaR(el.p4(), particle->p4());
0314     if (dRtmp < dR) {
0315       dR = dRtmp;
0316       closestElectron = particle;
0317     }
0318   }
0319   // See if the closest electron (if it exists) is close enough.
0320   // If not, no match found.
0321   if (!(closestElectron != 0 && dR < 0.1)) {
0322     return UNMATCHED;
0323   }
0324
0325   //
0326   int ancestorPID = -999;
0327   int ancestorStatus = -999;
0328   findFirstNonElectronMother(closestElectron, ancestorPID, ancestorStatus);
0329
0330   if (ancestorPID == -999 && ancestorStatus == -999) {
0331     // No non-electron parent??? This should never happen.
0332     // Complain.
0333     printf("ElectronNtupler: ERROR! Electron does not apper to have a non-electron parent\n");
0334     return UNMATCHED;
0335   }
0336
0337   if (abs(ancestorPID) > 50 && ancestorStatus == 2)
0338     return TRUE_NON_PROMPT_ELECTRON;
0339
0340   if (abs(ancestorPID) == 15 && ancestorStatus == 2)
0341     return TRUE_ELECTRON_FROM_TAU;
0342
0343   // What remains is true prompt electrons
0344   return TRUE_PROMPT_ELECTRON;
0345 }
0346
0347 void MiniAODElectronIDValidationAnalyzer::findFirstNonElectronMother(const reco::Candidate *particle,
0348                                                                      int &ancestorPID,
0349                                                                      int &ancestorStatus) {
0350   if (particle == 0) {
0351     printf("ElectronNtupler: ERROR! null candidate pointer, this should never happen\n");
0352     return;
0353   }
0354
0355   // Is this the first non-electron parent? If yes, return, otherwise
0356   // go deeper into recursion
0357   if (abs(particle->pdgId()) == 11) {
0358     findFirstNonElectronMother(particle->mother(0), ancestorPID, ancestorStatus);
0359   } else {
0360     ancestorPID = particle->pdgId();
0361     ancestorStatus = particle->status();
0362   }
0363
0364   return;
0365 }
0366
0367 //define this as a plug-in
0368 DEFINE_FWK_MODULE(MiniAODElectronIDValidationAnalyzer);