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
0064   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0065
0066   enum ElectronMatchType {
0067     UNMATCHED = 0,
0068     TRUE_PROMPT_ELECTRON,
0069     TRUE_ELECTRON_FROM_TAU,
0070     TRUE_NON_PROMPT_ELECTRON
0071   };  // The last does not include tau parents
0072
0073 private:
0074   void analyze(const edm::Event &, const edm::EventSetup &) override;
0075
0076   int matchToTruth(const reco::GsfElectron &el, const edm::Handle<edm::View<reco::GenParticle>> &genParticles);
0077   void findFirstNonElectronMother(const reco::Candidate *particle, int &ancestorPID, int &ancestorStatus);
0078
0079   // ----------member data ---------------------------
0080   edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0081   edm::EDGetTokenT<edm::View<reco::GenParticle>> genToken_;
0082   edm::EDGetTokenT<reco::ConversionCollection> convToken_;
0083   edm::EDGetTokenT<reco::BeamSpot> beamToken_;
0084   edm::EDGetTokenT<edm::ValueMap<float>> full5x5SigmaIEtaIEtaMapToken_;
0085   edm::EDGetTokenT<edm::View<pat::Electron>> electronCollectionToken_;
0086   edm::InputTag electronIdTag_;
0087   edm::EDGetTokenT<edm::ValueMap<bool>> electronIdToken_;
0088
0089   TTree *electronTree_;
0090   Float_t pt_;
0091   Float_t etaSC_;
0092   // All ID variables
0093   Float_t dEtaIn_;
0094   Float_t dPhiIn_;
0095   Float_t hOverE_;
0096   Float_t sigmaIetaIeta_;
0097   Float_t full5x5_sigmaIetaIeta_;
0098   Float_t relIsoWithDBeta_;
0099   Float_t ooEmooP_;
0100   Float_t d0_;
0101   Float_t dz_;
0102   Int_t expectedMissingInnerHits_;
0103   Int_t passConversionVeto_;
0104
0105   Int_t isTrue_;
0106   Int_t isPass_;
0107 };
0108
0109 //
0110 // constants, enums and typedefs
0111 //
0112
0113 //
0114 // static data member definitions
0115 //
0116
0117 //
0118 // constructors and destructor
0119 //
0120 MiniAODElectronIDValidationAnalyzer::MiniAODElectronIDValidationAnalyzer(const edm::ParameterSet &iConfig)
0121     : vtxToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0122       genToken_(consumes<edm::View<reco::GenParticle>>(iConfig.getParameter<edm::InputTag>("genparticles"))),
0123       convToken_(consumes<reco::ConversionCollection>(iConfig.getParameter<edm::InputTag>("convcollection"))),
0124       beamToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamspot"))),
0125       full5x5SigmaIEtaIEtaMapToken_(
0126           consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("full5x5SigmaIEtaIEtaMap"))),
0127       electronCollectionToken_(consumes<edm::View<pat::Electron>>(iConfig.getParameter<edm::InputTag>("electrons"))),
0128       electronIdTag_(iConfig.getParameter<edm::InputTag>("electronIDs")),
0129       electronIdToken_(consumes<edm::ValueMap<bool>>(iConfig.getParameter<edm::InputTag>("electronIDs"))) {
0130   //now do what ever initialization is needed
0131   usesResource(TFileService::kSharedResource);
0132   edm::Service<TFileService> fs;
0133   electronTree_ = fs->make<TTree>("ElectronTree", "Electron data");
0134
0135   electronTree_->Branch("pt", &pt_, "pt/F");
0136   electronTree_->Branch("etaSC", &etaSC_, "etaSC/F");
0137
0138   electronTree_->Branch("dEtaIn", &dEtaIn_, "dEtaIn/F");
0139   electronTree_->Branch("dPhiIn", &dPhiIn_, "dPhiIn/F");
0140   electronTree_->Branch("hOverE", &hOverE_, "hOverE/F");
0141   electronTree_->Branch("sigmaIetaIeta", &sigmaIetaIeta_, "sigmaIetaIeta/F");
0142   electronTree_->Branch("full5x5_sigmaIetaIeta", &full5x5_sigmaIetaIeta_, "full5x5_sigmaIetaIeta/F");
0143   electronTree_->Branch("relIsoWithDBeta", &relIsoWithDBeta_, "relIsoWithDBeta/F");
0144   electronTree_->Branch("ooEmooP", &ooEmooP_, "ooEmooP/F");
0145   electronTree_->Branch("d0", &d0_, "d0/F");
0146   electronTree_->Branch("dz", &dz_, "dz/F");
0147   electronTree_->Branch("expectedMissingInnerHits", &expectedMissingInnerHits_, "expectedMissingInnerHits/I");
0148   electronTree_->Branch("passConversionVeto", &passConversionVeto_, "passConversionVeto/I");
0149
0150   electronTree_->Branch("isTrue", &isTrue_, "isTrue/I");
0151   electronTree_->Branch("isPass", &isPass_, "isPass/I");
0152 }
0153
0154 //
0155 // member functions
0156 //
0157
0158 // ------------ method called for each event  ------------
0159 void MiniAODElectronIDValidationAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0160   //edm::Handle<edm::ValueMap<float> > full5x5sieie;
0161   edm::Handle<edm::View<pat::Electron>> collection;
0162   //edm::Handle<edm::ValueMap<bool> > id_decisions;
0163   //iEvent.getByToken(full5x5SigmaIEtaIEtaMapToken_,full5x5sieie);
0164   iEvent.getByToken(electronCollectionToken_, collection);
0165   //iEvent.getByToken(electronIdToken_,id_decisions);
0166
0167   // Get PV
0168   edm::Handle<reco::VertexCollection> vertices;
0169   iEvent.getByToken(vtxToken_, vertices);
0170   if (vertices->empty())
0171     return;  // skip the event if no PV found
0172   const reco::Vertex &pv = vertices->front();
0173
0174   // Get GenParticles
0175   edm::Handle<edm::View<reco::GenParticle>> genParticles;
0176   iEvent.getByToken(genToken_, genParticles);
0177
0178   // Get stuff for conversions
0179   edm::Handle<reco::ConversionCollection> convs;
0180   edm::Handle<reco::BeamSpot> thebs;
0181   iEvent.getByToken(convToken_, convs);
0182   iEvent.getByToken(beamToken_, thebs);
0183
0184   for (size_t i = 0; i < collection->size(); ++i) {
0185     auto el = collection->refAt(i);
0186     pt_ = el->pt();
0187     etaSC_ = el->superCluster()->eta();
0188
0189     // ID and matching
0190     dEtaIn_ = el->deltaEtaSuperClusterTrackAtVtx();
0191     dPhiIn_ = el->deltaPhiSuperClusterTrackAtVtx();
0192     hOverE_ = el->hcalOverEcal();
0193     sigmaIetaIeta_ = el->sigmaIetaIeta();
0194     full5x5_sigmaIetaIeta_ = el->full5x5_sigmaIetaIeta();
0195     // |1/E-1/p| = |1/E - EoverPinner/E| is computed below
0196     // The if protects against ecalEnergy == inf or zero (always
0197     // the case for electrons below 5 GeV in miniAOD)
0198     if (el->ecalEnergy() == 0) {
0199       printf("Electron energy is zero!\n");
0200       ooEmooP_ = 1e30;
0201     } else if (!std::isfinite(el->ecalEnergy())) {
0202       printf("Electron energy is not finite!\n");
0203       ooEmooP_ = 1e30;
0204     } else {
0205       ooEmooP_ = fabs(1.0 / el->ecalEnergy() - el->eSuperClusterOverP() / el->ecalEnergy());
0206     }
0207
0208     // Isolation
0209     reco::GsfElectron::PflowIsolationVariables pfIso = el->pfIsolationVariables();
0210     // Compute isolation with delta beta correction for PU
0211     float absiso =
0212         pfIso.sumChargedHadronPt + std::max(0.0, pfIso.sumNeutralHadronEt + pfIso.sumPhotonEt - 0.5 * pfIso.sumPUPt);
0213     relIsoWithDBeta_ = absiso / pt_;
0214
0215     // Impact parameter
0216     d0_ = (-1) * el->gsfTrack()->dxy(pv.position());
0217     dz_ = el->gsfTrack()->dz(pv.position());
0218
0219     // Conversion rejection
0220     constexpr reco::HitPattern::HitCategory missingHitType = reco::HitPattern::MISSING_INNER_HITS;
0221     expectedMissingInnerHits_ = el->gsfTrack()->hitPattern().numberOfLostHits(missingHitType);
0222     passConversionVeto_ = false;
0223     if (thebs.isValid() && convs.isValid()) {
0224       passConversionVeto_ = !ConversionTools::hasMatchedConversion(*el, *convs, thebs->position());
0225     } else {
0226       printf("\n\nERROR!!! conversions not found!!!\n");
0227     }
0228
0229     isTrue_ = matchToTruth(*el, genParticles);
0230     isPass_ = (el->electronID(electronIdTag_.encode()) > 0.5);
0231
0232     electronTree_->Fill();
0233   }
0234 }
0235
0236 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0237 void MiniAODElectronIDValidationAnalyzer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0238   //The following says we do not know what parameters are allowed so do no validation
0239   // Please change this to state exactly what you do use, even if it is no parameters
0240   edm::ParameterSetDescription desc;
0241   desc.setUnknown();
0242   descriptions.addDefault(desc);
0243 }
0244
0245 // The function that uses algorith from Josh Bendavid with
0246 // an explicit loop over gen particles.
0247 int MiniAODElectronIDValidationAnalyzer::matchToTruth(const reco::GsfElectron &el,
0248                                                       const edm::Handle<edm::View<reco::GenParticle>> &genParticles) {
0249   //
0250   // Explicit loop and geometric matching method (advised by Josh Bendavid)
0251   //
0252
0253   // Find the closest status 1 gen electron to the reco electron
0254   double dR = 999;
0255   const reco::Candidate *closestElectron = 0;
0256   for (size_t i = 0; i < genParticles->size(); i++) {
0257     const reco::Candidate *particle = &(*genParticles)[i];
0258     // Drop everything that is not electron or not status 1
0259     if (abs(particle->pdgId()) != 11 || particle->status() != 1)
0260       continue;
0261     //
0262     double dRtmp = ROOT::Math::VectorUtil::DeltaR(el.p4(), particle->p4());
0263     if (dRtmp < dR) {
0264       dR = dRtmp;
0265       closestElectron = particle;
0266     }
0267   }
0268   // See if the closest electron (if it exists) is close enough.
0269   // If not, no match found.
0270   if (!(closestElectron != 0 && dR < 0.1)) {
0271     return UNMATCHED;
0272   }
0273
0274   //
0275   int ancestorPID = -999;
0276   int ancestorStatus = -999;
0277   findFirstNonElectronMother(closestElectron, ancestorPID, ancestorStatus);
0278
0279   if (ancestorPID == -999 && ancestorStatus == -999) {
0280     // No non-electron parent??? This should never happen.
0281     // Complain.
0282     printf("ElectronNtupler: ERROR! Electron does not apper to have a non-electron parent\n");
0283     return UNMATCHED;
0284   }
0285
0286   if (abs(ancestorPID) > 50 && ancestorStatus == 2)
0287     return TRUE_NON_PROMPT_ELECTRON;
0288
0289   if (abs(ancestorPID) == 15 && ancestorStatus == 2)
0290     return TRUE_ELECTRON_FROM_TAU;
0291
0292   // What remains is true prompt electrons
0293   return TRUE_PROMPT_ELECTRON;
0294 }
0295
0296 void MiniAODElectronIDValidationAnalyzer::findFirstNonElectronMother(const reco::Candidate *particle,
0297                                                                      int &ancestorPID,
0298                                                                      int &ancestorStatus) {
0299   if (particle == 0) {
0300     printf("ElectronNtupler: ERROR! null candidate pointer, this should never happen\n");
0301     return;
0302   }
0303
0304   // Is this the first non-electron parent? If yes, return, otherwise
0305   // go deeper into recursion
0306   if (abs(particle->pdgId()) == 11) {
0307     findFirstNonElectronMother(particle->mother(0), ancestorPID, ancestorStatus);
0308   } else {
0309     ancestorPID = particle->pdgId();
0310     ancestorStatus = particle->status();
0311   }
0312
0313   return;
0314 }
0315
0316 //define this as a plug-in
0317 DEFINE_FWK_MODULE(MiniAODElectronIDValidationAnalyzer);