Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:14:03

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/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::EDAnalyzer {
0061 public:
0062   explicit MiniAODElectronIDValidationAnalyzer(const edm::ParameterSet &);
0063   ~MiniAODElectronIDValidationAnalyzer();
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   virtual void beginJob() override;
0076   virtual void analyze(const edm::Event &, const edm::EventSetup &) override;
0077   virtual 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   edm::Service<TFileService> fs;
0140   electronTree_ = fs->make<TTree>("ElectronTree", "Electron data");
0141 
0142   electronTree_->Branch("pt", &pt_, "pt/F");
0143   electronTree_->Branch("etaSC", &etaSC_, "etaSC/F");
0144 
0145   electronTree_->Branch("dEtaIn", &dEtaIn_, "dEtaIn/F");
0146   electronTree_->Branch("dPhiIn", &dPhiIn_, "dPhiIn/F");
0147   electronTree_->Branch("hOverE", &hOverE_, "hOverE/F");
0148   electronTree_->Branch("sigmaIetaIeta", &sigmaIetaIeta_, "sigmaIetaIeta/F");
0149   electronTree_->Branch("full5x5_sigmaIetaIeta", &full5x5_sigmaIetaIeta_, "full5x5_sigmaIetaIeta/F");
0150   electronTree_->Branch("relIsoWithDBeta", &relIsoWithDBeta_, "relIsoWithDBeta/F");
0151   electronTree_->Branch("ooEmooP", &ooEmooP_, "ooEmooP/F");
0152   electronTree_->Branch("d0", &d0_, "d0/F");
0153   electronTree_->Branch("dz", &dz_, "dz/F");
0154   electronTree_->Branch("expectedMissingInnerHits", &expectedMissingInnerHits_, "expectedMissingInnerHits/I");
0155   electronTree_->Branch("passConversionVeto", &passConversionVeto_, "passConversionVeto/I");
0156 
0157   electronTree_->Branch("isTrue", &isTrue_, "isTrue/I");
0158   electronTree_->Branch("isPass", &isPass_, "isPass/I");
0159 }
0160 
0161 MiniAODElectronIDValidationAnalyzer::~MiniAODElectronIDValidationAnalyzer() {
0162   // do anything here that needs to be done at desctruction time
0163   // (e.g. close files, deallocate resources etc.)
0164 }
0165 
0166 //
0167 // member functions
0168 //
0169 
0170 // ------------ method called for each event  ------------
0171 void MiniAODElectronIDValidationAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0172   //edm::Handle<edm::ValueMap<float> > full5x5sieie;
0173   edm::Handle<edm::View<pat::Electron>> collection;
0174   //edm::Handle<edm::ValueMap<bool> > id_decisions;
0175   //iEvent.getByToken(full5x5SigmaIEtaIEtaMapToken_,full5x5sieie);
0176   iEvent.getByToken(electronCollectionToken_, collection);
0177   //iEvent.getByToken(electronIdToken_,id_decisions);
0178 
0179   // Get PV
0180   edm::Handle<reco::VertexCollection> vertices;
0181   iEvent.getByToken(vtxToken_, vertices);
0182   if (vertices->empty())
0183     return;  // skip the event if no PV found
0184   const reco::Vertex &pv = vertices->front();
0185 
0186   // Get GenParticles
0187   edm::Handle<edm::View<reco::GenParticle>> genParticles;
0188   iEvent.getByToken(genToken_, genParticles);
0189 
0190   // Get stuff for conversions
0191   edm::Handle<reco::ConversionCollection> convs;
0192   edm::Handle<reco::BeamSpot> thebs;
0193   iEvent.getByToken(convToken_, convs);
0194   iEvent.getByToken(beamToken_, thebs);
0195 
0196   for (size_t i = 0; i < collection->size(); ++i) {
0197     auto el = collection->refAt(i);
0198     pt_ = el->pt();
0199     etaSC_ = el->superCluster()->eta();
0200 
0201     // ID and matching
0202     dEtaIn_ = el->deltaEtaSuperClusterTrackAtVtx();
0203     dPhiIn_ = el->deltaPhiSuperClusterTrackAtVtx();
0204     hOverE_ = el->hcalOverEcal();
0205     sigmaIetaIeta_ = el->sigmaIetaIeta();
0206     full5x5_sigmaIetaIeta_ = el->full5x5_sigmaIetaIeta();
0207     // |1/E-1/p| = |1/E - EoverPinner/E| is computed below
0208     // The if protects against ecalEnergy == inf or zero (always
0209     // the case for electrons below 5 GeV in miniAOD)
0210     if (el->ecalEnergy() == 0) {
0211       printf("Electron energy is zero!\n");
0212       ooEmooP_ = 1e30;
0213     } else if (!std::isfinite(el->ecalEnergy())) {
0214       printf("Electron energy is not finite!\n");
0215       ooEmooP_ = 1e30;
0216     } else {
0217       ooEmooP_ = fabs(1.0 / el->ecalEnergy() - el->eSuperClusterOverP() / el->ecalEnergy());
0218     }
0219 
0220     // Isolation
0221     reco::GsfElectron::PflowIsolationVariables pfIso = el->pfIsolationVariables();
0222     // Compute isolation with delta beta correction for PU
0223     float absiso =
0224         pfIso.sumChargedHadronPt + std::max(0.0, pfIso.sumNeutralHadronEt + pfIso.sumPhotonEt - 0.5 * pfIso.sumPUPt);
0225     relIsoWithDBeta_ = absiso / pt_;
0226 
0227     // Impact parameter
0228     d0_ = (-1) * el->gsfTrack()->dxy(pv.position());
0229     dz_ = el->gsfTrack()->dz(pv.position());
0230 
0231     // Conversion rejection
0232     constexpr reco::HitPattern::HitCategory missingHitType = reco::HitPattern::MISSING_INNER_HITS;
0233     expectedMissingInnerHits_ = el->gsfTrack()->hitPattern().numberOfLostHits(missingHitType);
0234     passConversionVeto_ = false;
0235     if (thebs.isValid() && convs.isValid()) {
0236       passConversionVeto_ = !ConversionTools::hasMatchedConversion(*el, *convs, thebs->position());
0237     } else {
0238       printf("\n\nERROR!!! conversions not found!!!\n");
0239     }
0240 
0241     isTrue_ = matchToTruth(*el, genParticles);
0242     isPass_ = (el->electronID(electronIdTag_.encode()) > 0.5);
0243 
0244     electronTree_->Fill();
0245   }
0246 }
0247 
0248 // ------------ method called once each job just before starting event loop  ------------
0249 void MiniAODElectronIDValidationAnalyzer::beginJob() {}
0250 
0251 // ------------ method called once each job just after ending the event loop  ------------
0252 void MiniAODElectronIDValidationAnalyzer::endJob() {}
0253 
0254 // ------------ method called when starting to processes a run  ------------
0255 /*
0256 void 
0257 MiniAODElectronIDValidationAnalyzer::beginRun(edm::Run const&, edm::EventSetup const&)
0258 {
0259 }
0260 */
0261 
0262 // ------------ method called when ending the processing of a run  ------------
0263 /*
0264 void 
0265 MiniAODElectronIDValidationAnalyzer::endRun(edm::Run const&, edm::EventSetup const&)
0266 {
0267 }
0268 */
0269 
0270 // ------------ method called when starting to processes a luminosity block  ------------
0271 /*
0272 void 
0273 MiniAODElectronIDValidationAnalyzer::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0274 {
0275 }
0276 */
0277 
0278 // ------------ method called when ending the processing of a luminosity block  ------------
0279 /*
0280 void 
0281 MiniAODElectronIDValidationAnalyzer::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0282 {
0283 }
0284 */
0285 
0286 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0287 void MiniAODElectronIDValidationAnalyzer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0288   //The following says we do not know what parameters are allowed so do no validation
0289   // Please change this to state exactly what you do use, even if it is no parameters
0290   edm::ParameterSetDescription desc;
0291   desc.setUnknown();
0292   descriptions.addDefault(desc);
0293 }
0294 
0295 // The function that uses algorith from Josh Bendavid with
0296 // an explicit loop over gen particles.
0297 int MiniAODElectronIDValidationAnalyzer::matchToTruth(const reco::GsfElectron &el,
0298                                                       const edm::Handle<edm::View<reco::GenParticle>> &genParticles) {
0299   //
0300   // Explicit loop and geometric matching method (advised by Josh Bendavid)
0301   //
0302 
0303   // Find the closest status 1 gen electron to the reco electron
0304   double dR = 999;
0305   const reco::Candidate *closestElectron = 0;
0306   for (size_t i = 0; i < genParticles->size(); i++) {
0307     const reco::Candidate *particle = &(*genParticles)[i];
0308     // Drop everything that is not electron or not status 1
0309     if (abs(particle->pdgId()) != 11 || particle->status() != 1)
0310       continue;
0311     //
0312     double dRtmp = ROOT::Math::VectorUtil::DeltaR(el.p4(), particle->p4());
0313     if (dRtmp < dR) {
0314       dR = dRtmp;
0315       closestElectron = particle;
0316     }
0317   }
0318   // See if the closest electron (if it exists) is close enough.
0319   // If not, no match found.
0320   if (!(closestElectron != 0 && dR < 0.1)) {
0321     return UNMATCHED;
0322   }
0323 
0324   //
0325   int ancestorPID = -999;
0326   int ancestorStatus = -999;
0327   findFirstNonElectronMother(closestElectron, ancestorPID, ancestorStatus);
0328 
0329   if (ancestorPID == -999 && ancestorStatus == -999) {
0330     // No non-electron parent??? This should never happen.
0331     // Complain.
0332     printf("ElectronNtupler: ERROR! Electron does not apper to have a non-electron parent\n");
0333     return UNMATCHED;
0334   }
0335 
0336   if (abs(ancestorPID) > 50 && ancestorStatus == 2)
0337     return TRUE_NON_PROMPT_ELECTRON;
0338 
0339   if (abs(ancestorPID) == 15 && ancestorStatus == 2)
0340     return TRUE_ELECTRON_FROM_TAU;
0341 
0342   // What remains is true prompt electrons
0343   return TRUE_PROMPT_ELECTRON;
0344 }
0345 
0346 void MiniAODElectronIDValidationAnalyzer::findFirstNonElectronMother(const reco::Candidate *particle,
0347                                                                      int &ancestorPID,
0348                                                                      int &ancestorStatus) {
0349   if (particle == 0) {
0350     printf("ElectronNtupler: ERROR! null candidate pointer, this should never happen\n");
0351     return;
0352   }
0353 
0354   // Is this the first non-electron parent? If yes, return, otherwise
0355   // go deeper into recursion
0356   if (abs(particle->pdgId()) == 11) {
0357     findFirstNonElectronMother(particle->mother(0), ancestorPID, ancestorStatus);
0358   } else {
0359     ancestorPID = particle->pdgId();
0360     ancestorStatus = particle->status();
0361   }
0362 
0363   return;
0364 }
0365 
0366 //define this as a plug-in
0367 DEFINE_FWK_MODULE(MiniAODElectronIDValidationAnalyzer);