PatExamples/plugins/PatTauAnalyzer.cc

0001 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0002 #include "DataFormats/PatCandidates/interface/Tau.h"
0003 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ServiceRegistry/interface/Service.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011
0012 #include <TH1.h>
0013 #include <TMath.h>
0014
0015 #include <string>
0016
0017 class PatTauAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0018 public:
0019   explicit PatTauAnalyzer(const edm::ParameterSet&);
0020   ~PatTauAnalyzer() override;
0021
0022   //--- methods inherited from EDAnalyzer base-class
0023   void beginJob() override;
0024   void analyze(const edm::Event&, const edm::EventSetup&) override;
0025   void endJob() override;
0026
0027 private:
0028   //--- configuration parameters
0029   edm::InputTag src_;
0030   edm::EDGetTokenT<pat::TauCollection> srcToken_;
0031
0032   bool requireGenTauMatch_;
0033
0034   std::string discrByLeadTrack_;
0035   std::string discrByIso_;
0036   std::string discrByTaNC_;
0037
0038   //--- generator level histograms
0039   TH1* hGenTauEnergy_;
0040   TH1* hGenTauPt_;
0041   TH1* hGenTauEta_;
0042   TH1* hGenTauPhi_;
0043
0044   //--- reconstruction level histograms
0045   TH1* hTauJetEnergy_;
0046   TH1* hTauJetPt_;
0047   TH1* hTauJetEta_;
0048   TH1* hTauJetPhi_;
0049
0050   TH1* hNumTauJets_;
0051
0052   TH1* hTauLeadTrackPt_;
0053
0054   TH1* hTauNumSigConeTracks_;
0055   TH1* hTauNumIsoConeTracks_;
0056
0057   TH1* hTauDiscrByIso_;
0058   TH1* hTauDiscrByTaNC_;
0059   TH1* hTauDiscrAgainstElectrons_;
0060   TH1* hTauDiscrAgainstMuons_;
0061
0062   TH1* hTauJetEnergyIsoPassed_;
0063   TH1* hTauJetPtIsoPassed_;
0064   TH1* hTauJetEtaIsoPassed_;
0065   TH1* hTauJetPhiIsoPassed_;
0066
0067   TH1* hTauJetEnergyTaNCpassed_;
0068   TH1* hTauJetPtTaNCpassed_;
0069   TH1* hTauJetEtaTaNCpassed_;
0070   TH1* hTauJetPhiTaNCpassed_;
0071 };
0072
0073 const reco::GenParticle* getGenTau(const pat::Tau& patTau) {
0074   std::vector<reco::GenParticleRef> associatedGenParticles = patTau.genParticleRefs();
0075   for (std::vector<reco::GenParticleRef>::const_iterator it = associatedGenParticles.begin();
0076        it != associatedGenParticles.end();
0077        ++it) {
0078     if (it->isAvailable()) {
0079       const reco::GenParticleRef& genParticle = (*it);
0080       if (genParticle->pdgId() == -15 || genParticle->pdgId() == +15)
0081         return genParticle.get();
0082     }
0083   }
0084
0085   return nullptr;
0086 }
0087
0088 PatTauAnalyzer::PatTauAnalyzer(const edm::ParameterSet& cfg) {
0089   //std::cout << "<PatTauAnalyzer::PatTauAnalyzer>:" << std::endl;
0090   usesResource(TFileService::kSharedResource);
0091
0092   //--- read name of pat::Tau collection
0093   src_ = cfg.getParameter<edm::InputTag>("src");
0094   //std::cout << " src = " << src_ << std::endl;
0095   srcToken_ = consumes<pat::TauCollection>(src_);
0096
0097   //--- fill histograms for all tau-jet candidates or for "real" taus only ?
0098   requireGenTauMatch_ = cfg.getParameter<bool>("requireGenTauMatch");
0099   //std::cout << " requireGenTauMatch = " << requireGenTauMatch_ << std::endl;
0100
0101   //--- read names of tau id. discriminators
0102   discrByLeadTrack_ = cfg.getParameter<std::string>("discrByLeadTrack");
0103   //std::cout << " discrByLeadTrack = " << discrByLeadTrack_ << std::endl;
0104
0105   discrByIso_ = cfg.getParameter<std::string>("discrByIso");
0106   //std::cout << " discrByIso = " << discrByIso_ << std::endl;
0107
0108   discrByTaNC_ = cfg.getParameter<std::string>("discrByTaNC");
0109   //std::cout << " discrByTaNC = " << discrByTaNC_ << std::endl;
0110 }
0111
0112 PatTauAnalyzer::~PatTauAnalyzer() {
0113   //std::cout << "<PatTauAnalyzer::~PatTauAnalyzer>:" << std::endl;
0114
0115   //--- clean-up memory;
0116   //    delete all histograms
0117   /*
0118   deletion of histograms taken care of by TFileService;
0119   do not delete them here (if the histograms are deleted here,
0120   they will not appear in the ROOT file written by TFileService)
0121
0122   delete hGenTauEnergy_;
0123   delete hGenTauPt_;
0124   delete hGenTauEta_;
0125   delete hGenTauPhi_;
0126   delete hTauJetEnergy_;
0127   delete hTauJetPt_;
0128   delete hTauJetEta_;
0129   delete hTauJetPhi_;
0130   delete hNumTauJets_;
0131   delete hTauLeadTrackPt_;
0132   delete hTauNumSigConeTracks_;
0133   delete hTauNumIsoConeTracks_;
0134   delete hTauDiscrByIso_;
0135   delete hTauDiscrByTaNC_;
0136   delete hTauDiscrAgainstElectrons_;
0137   delete hTauDiscrAgainstMuons_;
0138   delete hTauJetEnergyIsoPassed_;
0139   delete hTauJetPtIsoPassed_;
0140   delete hTauJetEtaIsoPassed_;
0141   delete hTauJetPhiIsoPassed_;
0142   delete hTauJetEnergyTaNCpassed_;
0143   delete hTauJetPtTaNCpassed_;
0144   delete hTauJetEtaTaNCpassed_;
0145   delete hTauJetPhiTaNCpassed_;
0146  */
0147 }
0148
0149 void PatTauAnalyzer::beginJob() {
0150   //--- retrieve handle to auxiliary service
0151   //    used for storing histograms into ROOT file
0152   edm::Service<TFileService> fs;
0153
0154   //--- book generator level histograms
0155   hGenTauEnergy_ = fs->make<TH1F>("GenTauEnergy", "GenTauEnergy", 30, 0., 150.);
0156   hGenTauPt_ = fs->make<TH1F>("GenTauPt", "GenTauPt", 30, 0., 150.);
0157   hGenTauEta_ = fs->make<TH1F>("GenTauEta", "GenTauEta", 24, -3., +3.);
0158   hGenTauPhi_ = fs->make<TH1F>("GenTauPhi", "GenTauPhi", 18, -TMath::Pi(), +TMath::Pi());
0159
0160   //--- book reconstruction level histograms
0161   //    for tau-jet Energy, Pt, Eta, Phi
0162   hTauJetEnergy_ = fs->make<TH1F>("TauJetEnergy", "TauJetEnergy", 30, 0., 150.);
0163   hTauJetPt_ = fs->make<TH1F>("TauJetPt", "TauJetPt", 30, 0., 150.);
0164   //
0165   // TO-DO: add histograms for eta and phi of the tau-jet candidate
0166   //
0167   // NOTE:
0168   //  1.) please use
0169   //       "TauJetEta" and "TauJetPhi"
0170   //      for the names of the histograms and choose the exact same binning
0171   //      as is used for the histograms
0172   //       "TauJetEtaIsoPassed" and "TauJetPhiIsoPassed"
0173   //      below
0174   //
0175   //  2.) please check the histograms
0176   //       hTauJetEta_ and hTauJetPt_
0177   //      have already been defined in PatTauAnalyzer.h
0178   //
0179   //hTauJetEta_ =...
0180   //hTauJetPt_ =...
0181
0182   //... for number of tau-jet candidates
0183   hNumTauJets_ = fs->make<TH1F>("NumTauJets", "NumTauJets", 10, -0.5, 9.5);
0184
0185   //... for Pt of highest Pt track within signal cone tau-jet...
0186   hTauLeadTrackPt_ = fs->make<TH1F>("TauLeadTrackPt", "TauLeadTrackPt", 40, 0., 100.);
0187
0188   //... for total number of tracks within signal/isolation cones
0189   hTauNumSigConeTracks_ = fs->make<TH1F>("TauNumSigConeTracks", "TauNumSigConeTracks", 10, -0.5, 9.5);
0190   hTauNumIsoConeTracks_ = fs->make<TH1F>("TauNumIsoConeTracks", "TauNumIsoConeTracks", 20, -0.5, 19.5);
0191
0192   //... for values of tau id. discriminators based on track isolation cut/
0193   //    neural network-based tau id.
0194   hTauDiscrByIso_ = fs->make<TH1F>("TauDiscrByIso", "TauDiscrByIso", 103, -0.015, 1.015);
0195   hTauDiscrByTaNC_ = fs->make<TH1F>("TauDiscrByTaNC", "TauDiscrByTaNC", 103, -0.015, 1.015);
0196
0197   //... for values of tau id. discriminators against (unidentified) electrons and muons
0198   hTauDiscrAgainstElectrons_ =
0199       fs->make<TH1F>("TauDiscrAgainstElectrons", "TauDiscrAgainstElectrons", 103, -0.015, 1.015);
0200   hTauDiscrAgainstMuons_ = fs->make<TH1F>("TauDiscrAgainstMuons", "TauDiscrAgainstMuons", 103, -0.015, 1.015);
0201
0202   //... for Energy, Pt, Eta, Phi of tau-jets passing the discriminatorByIsolation selection
0203   hTauJetEnergyIsoPassed_ = fs->make<TH1F>("TauJetEnergyIsoPassed", "TauJetEnergyIsoPassed", 30, 0., 150.);
0204   hTauJetPtIsoPassed_ = fs->make<TH1F>("TauJetPtIsoPassed", "TauJetPtIsoPassed", 30, 0., 150.);
0205   hTauJetEtaIsoPassed_ = fs->make<TH1F>("TauJetEtaIsoPassed", "TauJetEtaIsoPassed", 24, -3., +3.);
0206   hTauJetPhiIsoPassed_ = fs->make<TH1F>("TauJetPhiIsoPassed", "TauJetPhiIsoPassed", 18, -TMath::Pi(), +TMath::Pi());
0207
0208   //... for Energy, Pt, Eta, Phi of tau-jets passing the discriminatorByTaNC ("Tau Neural Classifier") selection
0209   hTauJetEnergyTaNCpassed_ = fs->make<TH1F>("TauJetEnergyTaNCpassed", "TauJetEnergyTaNCpassed", 30, 0., 150.);
0210   hTauJetPtTaNCpassed_ = fs->make<TH1F>("TauJetPtTaNCpassed", "TauJetPtTaNCpassed", 30, 0., 150.);
0211   hTauJetEtaTaNCpassed_ = fs->make<TH1F>("TauJetEtaTaNCpassed", "TauJetEtaTaNCpassed", 24, -3., +3.);
0212   hTauJetPhiTaNCpassed_ = fs->make<TH1F>("TauJetPhiTaNCpassed", "TauJetPhiTaNCpassed", 18, -TMath::Pi(), +TMath::Pi());
0213 }
0214
0215 void PatTauAnalyzer::analyze(const edm::Event& evt, const edm::EventSetup& es) {
0216   //std::cout << "<PatTauAnalyzer::analyze>:" << std::endl;
0217
0218   edm::Handle<pat::TauCollection> patTaus;
0219   evt.getByToken(srcToken_, patTaus);
0220
0221   hNumTauJets_->Fill(patTaus->size());
0222
0223   for (pat::TauCollection::const_iterator patTau = patTaus->begin(); patTau != patTaus->end(); ++patTau) {
0224     //--- skip fake taus in case configuration parameters set to do so...
0225     const reco::GenParticle* genTau = getGenTau(*patTau);
0226     if (requireGenTauMatch_ && !genTau)
0227       continue;
0228
0229     //--- fill generator level histograms
0230     if (genTau) {
0231       hGenTauEnergy_->Fill(genTau->energy());
0232       hGenTauPt_->Fill(genTau->pt());
0233       hGenTauEta_->Fill(genTau->eta());
0234       hGenTauPhi_->Fill(genTau->phi());
0235     }
0236
0237     //--- fill reconstruction level histograms
0238     //    for Pt of highest Pt track within signal cone tau-jet...
0239     hTauJetEnergy_->Fill(patTau->energy());
0240     hTauJetPt_->Fill(patTau->pt());
0241     //
0242     // TO-DO:
0243     //  1.) fill histograms
0244     //       hTauJetEta_ and hTauJetPhi_
0245     //      with the pseudo-rapidity and azimuthal angle
0246     //      of the tau-jet candidate respectively
0247     //  hTauJetEta_->...
0248     //  hTauJetPhi_->...
0249     //
0250     //  2.) fill histogram
0251     //       hTauLeadTrackPt_
0252     //      with the transverse momentum of the highest Pt ("leading") track within the tau-jet
0253     //
0254     // NOTE:
0255     //  1.) please have a look at
0256     //       http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/DataFormats/Candidate/interface/Particle.h?revision=1.28&view=markup
0257     //      to find the methods for accessing eta and phi of the tau-jet
0258     //
0259     //  2.) please have a look at
0260     //       http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/DataFormats/PatCandidates/interface/Tau.h?revision=1.25&view=markup
0261     //      to find the method for accessing the leading track
0262     //
0263     //  3.) the method pat::Tau::leadTrack returns a reference (reco::TrackRef) to a reco::Track object
0264     //      this reference can be null (in case no high Pt track has been reconstructed within the tau-jet),
0265     //      so a check if the leadTrack exists is needed before dereferencing the reco::TrackRef via operator->
0266     //
0267     //  if ( patTau->leadTrack().isAvailable() ) hTauLeadTrackPt_->Fill(patTau->leadTrack()->pt());
0268
0269     //... for total number of tracks within signal/isolation cones
0270     hTauNumSigConeTracks_->Fill(patTau->signalTracks().size());
0271     hTauNumIsoConeTracks_->Fill(patTau->isolationTracks().size());
0272
0273     //... for values of tau id. discriminators based on track isolation cut/
0274     //    neural network-based tau id.
0275     //    (combine with requirement of at least one "leading" track of Pt > 5. GeV
0276     //     within the signal cone of the tau-jet)
0277     float discrByIso = (patTau->tauID(discrByLeadTrack_.data()) > 0.5) ? patTau->tauID(discrByIso_.data()) : 0.;
0278     hTauDiscrByIso_->Fill(discrByIso);
0279     float discrByTaNC = (patTau->tauID(discrByLeadTrack_.data()) > 0.5) ? patTau->tauID(discrByTaNC_.data()) : 0.;
0280     hTauDiscrByTaNC_->Fill(discrByTaNC);
0281
0282     //... for values of tau id. discriminators against (unidentified) electrons and muons
0283     //
0284     // TO-DO: fill histogram
0285     //         hTauDiscrAgainstElectrons_
0286     //        with the value of the discriminatorAgainstElectronsLoose
0287     //
0288     // NOTE:
0289     //  1.) please have a look at
0290     //       http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/DataFormats/PatCandidates/interface/Tau.h?revision=1.25&view=markup
0291     //      to find the method for accessing the tau id. information
0292     //
0293     //  2.) please have a look at
0294     //       http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/PhysicsTools/PatAlgos/python/tools/tauTools.py?revision=1.43&view=markup
0295     //      and convince yourself that the string "againstElectronLoose" needs to be passed as argument
0296     //      of the pat::Tau::tauID method
0297     //
0298     //  hTauDiscrAgainstElectrons_->Fill...
0299     hTauDiscrAgainstMuons_->Fill(patTau->tauID("againstMuonLoose"));
0300
0301     //... for Energy, Pt, Eta, Phi of tau-jets passing the discriminatorByIsolation selection
0302     if (discrByIso > 0.5) {
0303       hTauJetEnergyIsoPassed_->Fill(patTau->energy());
0304       hTauJetPtIsoPassed_->Fill(patTau->pt());
0305       hTauJetEtaIsoPassed_->Fill(patTau->eta());
0306       hTauJetPhiIsoPassed_->Fill(patTau->phi());
0307     }
0308
0309     //... for Energy, Pt, Eta, Phi of tau-jets passing the discriminatorByTaNC ("Tau Neural Classifier") selection
0310     if (discrByTaNC > 0.5) {
0311       hTauJetEnergyTaNCpassed_->Fill(patTau->energy());
0312       hTauJetPtTaNCpassed_->Fill(patTau->pt());
0313       hTauJetEtaTaNCpassed_->Fill(patTau->eta());
0314       hTauJetPhiTaNCpassed_->Fill(patTau->phi());
0315     }
0316   }
0317 }
0318
0319 void PatTauAnalyzer::endJob() {
0320   //--- nothing to be done yet...
0321 }
0322
0323 #include "FWCore/Framework/interface/MakerMacros.h"
0324
0325 DEFINE_FWK_MODULE(PatTauAnalyzer);