ElectronTools/test/ElectronIDValidationAnalyzer.cc

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