Benchmark/src/PFTauElecRejectionBenchmark.cc

0001 #include "RecoParticleFlow/Benchmark/interface/PFTauElecRejectionBenchmark.h"
0002
0003 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
0004 #define BOOK1D(name, title, nbinsx, lowx, highx) \
0005   h##name =                                      \
0006       db_ ? db_->book1D(#name, title, nbinsx, lowx, highx)->getTH1F() : new TH1F(#name, title, nbinsx, lowx, highx)
0007
0008 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
0009 #define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)                            \
0010   h##name = db_ ? db_->book2D(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)->getTH2F() \
0011                 : new TH2F(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
0012
0013 // all versions OK
0014 // preprocesor macro for setting axis titles
0015 #define SETAXES(name, xtitle, ytitle)    \
0016   h##name->GetXaxis()->SetTitle(xtitle); \
0017   h##name->GetYaxis()->SetTitle(ytitle)
0018
0019 using namespace reco;
0020 using namespace std;
0021
0022 PFTauElecRejectionBenchmark::PFTauElecRejectionBenchmark() : file_(nullptr) {}
0023
0024 PFTauElecRejectionBenchmark::~PFTauElecRejectionBenchmark() {
0025   if (file_)
0026     file_->Close();
0027 }
0028
0029 void PFTauElecRejectionBenchmark::write() {
0030   // Store the DAQ Histograms
0031   if (!outputFile_.empty()) {
0032     if (db_)
0033       db_->save(outputFile_);
0034     // use bare Root if no DQM (FWLite applications)
0035     else if (file_) {
0036       file_->Write(outputFile_.c_str());
0037       cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
0038       file_->Close();
0039     }
0040   } else
0041     cout << "No output file specified (" << outputFile_ << "). Results will not be saved!" << endl;
0042 }
0043
0044 void PFTauElecRejectionBenchmark::setup(string Filename,
0045                                         string benchmarkLabel,
0046                                         double maxDeltaR,
0047                                         double minRecoPt,
0048                                         double maxRecoAbsEta,
0049                                         double minMCPt,
0050                                         double maxMCAbsEta,
0051                                         string sGenMatchObjectLabel,
0052                                         bool applyEcalCrackCut,
0053                                         DQMStore* db_store) {
0054   maxDeltaR_ = maxDeltaR;
0055   benchmarkLabel_ = benchmarkLabel;
0056   outputFile_ = Filename;
0057   minRecoPt_ = minRecoPt;
0058   maxRecoAbsEta_ = maxRecoAbsEta;
0059   minMCPt_ = minMCPt;
0060   maxMCAbsEta_ = maxMCAbsEta;
0061   sGenMatchObjectLabel_ = sGenMatchObjectLabel;
0062   applyEcalCrackCut_ = applyEcalCrackCut;
0063
0064   file_ = nullptr;
0065   db_ = db_store;
0066
0067   // print parameters
0068   cout << "PFTauElecRejectionBenchmark Setup parameters ==============================================" << endl;
0069   cout << "Filename to write histograms " << Filename << endl;
0070   cout << "Benchmark label name " << benchmarkLabel_ << endl;
0071   cout << "maxDeltaRMax " << maxDeltaR << endl;
0072   cout << "minRecoPt " << minRecoPt_ << endl;
0073   cout << "maxRecoAbsEta " << maxRecoAbsEta_ << endl;
0074   cout << "minMCPt " << minMCPt_ << endl;
0075   cout << "maxMCAbsEta " << maxMCAbsEta_ << endl;
0076   cout << "applyEcalCrackCut " << applyEcalCrackCut_ << endl;
0077   cout << "sGenMatchObjectLabel " << sGenMatchObjectLabel_ << endl;
0078
0079   // Book histogram
0080
0081   // Establish DQM Store
0082   string path = "PFTask/Benchmarks/" + benchmarkLabel_ + "/";
0083   path += "Gen";
0084   if (db_) {
0085     db_->setCurrentFolder(path);
0086   } else {
0087     file_ = new TFile(outputFile_.c_str(), "recreate");
0088     cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. " << endl;
0089   }
0090
0091   // E/p
0092   BOOK1D(EoverP, "E/p", 100, 0., 4.);
0093   SETAXES(EoverP, "E/p", "Entries");
0094
0095   BOOK1D(EoverP_barrel, "E/p barrel", 100, 0., 4.);
0096   SETAXES(EoverP_barrel, "E/p barrel", "Entries");
0097
0098   BOOK1D(EoverP_endcap, "E/p endcap", 100, 0., 4.);
0099   SETAXES(EoverP_endcap, "E/p endcap", "Entries");
0100
0101   BOOK1D(EoverP_preid0, "E/p (preid=0)", 100, 0., 4.);
0102   SETAXES(EoverP_preid0, "E/p", "Entries");
0103
0104   BOOK1D(EoverP_preid1, "E/p (preid=1)", 100, 0., 4.);
0105   SETAXES(EoverP_preid1, "E/p", "Entries");
0106
0107   // H/p
0108   BOOK1D(HoverP, "H/p", 100, 0., 2.);
0109   SETAXES(HoverP, "H/p", "Entries");
0110
0111   BOOK1D(HoverP_barrel, "H/p barrel", 100, 0., 2.);
0112   SETAXES(HoverP_barrel, "H/p barrel", "Entries");
0113
0114   BOOK1D(HoverP_endcap, "H/p endcap", 100, 0., 2.);
0115   SETAXES(HoverP_endcap, "H/p endcap", "Entries");
0116
0117   BOOK1D(HoverP_preid0, "H/p (preid=0)", 100, 0., 2.);
0118   SETAXES(HoverP_preid0, "H/p", "Entries");
0119
0120   BOOK1D(HoverP_preid1, "H/p (preid=1)", 100, 0., 2.);
0121   SETAXES(HoverP_preid1, "H/p", "Entries");
0122
0123   // emfrac
0124   BOOK1D(Emfrac, "EM fraction", 100, 0., 1.01);
0125   SETAXES(Emfrac, "em fraction", "Entries");
0126
0127   BOOK1D(Emfrac_barrel, "EM fraction barrel", 100, 0., 1.01);
0128   SETAXES(Emfrac_barrel, "em fraction barrel", "Entries");
0129
0130   BOOK1D(Emfrac_endcap, "EM fraction endcap", 100, 0., 1.01);
0131   SETAXES(Emfrac_endcap, "em fraction endcap", "Entries");
0132
0133   BOOK1D(Emfrac_preid0, "EM fraction (preid=0)", 100, 0., 1.01);
0134   SETAXES(Emfrac_preid0, "em fraction", "Entries");
0135
0136   BOOK1D(Emfrac_preid1, "EM fraction (preid=1)", 100, 0., 1.01);
0137   SETAXES(Emfrac_preid1, "em fraction", "Entries");
0138
0139   // PreID
0140   BOOK1D(ElecPreID, "PFElectron PreID decision", 6, 0., 1.01);
0141   SETAXES(ElecPreID, "PFElectron PreID decision", "Entries");
0142
0143   // MVA
0144   BOOK1D(ElecMVA, "PFElectron MVA", 100, -1.01, 1.01);
0145   SETAXES(ElecMVA, "PFElectron MVA", "Entries");
0146
0147   // Discriminant
0148   BOOK1D(TauElecDiscriminant, "PFTau-Electron Discriminant", 6, 0., 1.01);
0149   SETAXES(TauElecDiscriminant, "PFTau-Electron Discriminant", "Entries");
0150
0151   // PFCand clusters
0152   BOOK1D(pfcand_deltaEta, "PFCand cluster dEta", 100, 0., 0.8);
0153   SETAXES(pfcand_deltaEta, "PFCand cluster #Delta(#eta)", "Entries");
0154
0155   BOOK1D(pfcand_deltaEta_weightE, "PFCand cluster dEta, energy weighted", 100, 0., 0.8);
0156   SETAXES(pfcand_deltaEta_weightE, "PFCand cluster #Delta(#eta)", "Entries");
0157
0158   BOOK1D(pfcand_deltaPhiOverQ, "PFCand cluster dPhi/q", 100, -0.8, 0.8);
0159   SETAXES(pfcand_deltaPhiOverQ, "PFCand cluster #Delta(#phi)/q", "Entries");
0160
0161   BOOK1D(pfcand_deltaPhiOverQ_weightE, "PFCand cluster dEta/q, energy weighted", 100, -0.8, 0.8);
0162   SETAXES(pfcand_deltaPhiOverQ_weightE, "PFCand cluster #Delta(#phi)/q", "Entries");
0163
0164   // Leading KF track
0165   BOOK1D(leadTk_pt, "leading KF track pt", 100, 0., 80.);
0166   SETAXES(leadTk_pt, "leading KF track p_{T} (GeV)", "Entries");
0167
0168   BOOK1D(leadTk_eta, "leading KF track eta", 100, -4., 4.);
0169   SETAXES(leadTk_eta, "leading KF track #eta", "Entries");
0170
0171   BOOK1D(leadTk_phi, "leading KF track phi", 100, -3.2, 3.2);
0172   SETAXES(leadTk_phi, "leading KF track #phi", "Entries");
0173
0174   // Leading Gsf track
0175   BOOK1D(leadGsfTk_pt, "leading Gsf track pt", 100, 0., 80.);
0176   SETAXES(leadGsfTk_pt, "leading Gsf track p_{T} (GeV)", "Entries");
0177
0178   BOOK1D(leadGsfTk_eta, "leading Gsf track eta", 100, -4., 4.);
0179   SETAXES(leadGsfTk_eta, "leading Gsf track #eta", "Entries");
0180
0181   BOOK1D(leadGsfTk_phi, "leading Gsf track phi", 100, -3.2, 3.2);
0182   SETAXES(leadGsfTk_phi, "leading Gsf track #phi", "Entries");
0183
0184   // H/p vs E/p
0185   BOOK2D(HoPvsEoP, "H/p vs. E/p", 100, 0., 2., 100, 0., 2.);
0186   SETAXES(HoPvsEoP, "E/p", "H/p");
0187
0188   BOOK2D(HoPvsEoP_preid0, "H/p vs. E/p (preid=0)", 100, 0., 2., 100, 0., 2.);
0189   SETAXES(HoPvsEoP_preid0, "E/p", "H/p");
0190
0191   BOOK2D(HoPvsEoP_preid1, "H/p vs. E/p (preid=0)", 100, 0., 2., 100, 0., 2.);
0192   SETAXES(HoPvsEoP_preid1, "E/p", "H/p");
0193
0194   // em fraction vs E/p
0195   BOOK2D(EmfracvsEoP, "emfrac vs. E/p", 100, 0., 2., 100, 0., 1.01);
0196   SETAXES(EmfracvsEoP, "E/p", "em fraction");
0197
0198   BOOK2D(EmfracvsEoP_preid0, "emfrac vs. E/p (preid=0)", 100, 0., 2., 100, 0., 1.01);
0199   SETAXES(EmfracvsEoP_preid0, "E/p", "em fraction");
0200
0201   BOOK2D(EmfracvsEoP_preid1, "emfrac vs. E/p (preid=0)", 100, 0., 2., 100, 0., 1.01);
0202   SETAXES(EmfracvsEoP_preid1, "E/p", "em fraction");
0203 }
0204
0205 void PFTauElecRejectionBenchmark::process(edm::Handle<edm::HepMCProduct> mcevt,
0206                                           edm::Handle<reco::PFTauCollection> pfTaus,
0207                                           edm::Handle<reco::PFTauDiscriminator> pfTauIsoDiscr,
0208                                           edm::Handle<reco::PFTauDiscriminator> pfTauElecDiscr) {
0209   // Find Gen Objects to be matched with
0210   HepMC::GenEvent* generated_event = new HepMC::GenEvent(*(mcevt->GetEvent()));
0211   _GenObjects.clear();
0212
0213   TLorentzVector taunet;
0214   HepMC::GenEvent::particle_iterator p;
0215   for (p = generated_event->particles_begin(); p != generated_event->particles_end(); p++) {
0216     if (std::abs((*p)->pdg_id()) == 15 && (*p)->status() == 2) {
0217       bool lept_decay = false;
0218       TLorentzVector tau((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
0219       HepMC::GenVertex::particle_iterator z = (*p)->end_vertex()->particles_begin(HepMC::descendants);
0220       for (; z != (*p)->end_vertex()->particles_end(HepMC::descendants); z++) {
0221         if (std::abs((*z)->pdg_id()) == 11 || std::abs((*z)->pdg_id()) == 13)
0222           lept_decay = true;
0223         if (std::abs((*z)->pdg_id()) == 16)
0224           taunet.SetPxPyPzE((*z)->momentum().px(), (*z)->momentum().py(), (*z)->momentum().pz(), (*z)->momentum().e());
0225       }
0226       if (lept_decay == false) {
0227         TLorentzVector jetMom = tau - taunet;
0228         if (sGenMatchObjectLabel_ == "tau")
0229           _GenObjects.push_back(jetMom);
0230       }
0231     } else if (std::abs((*p)->pdg_id()) == 11 && (*p)->status() == 1) {
0232       TLorentzVector elec((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
0233       if (sGenMatchObjectLabel_ == "e")
0234         _GenObjects.push_back(elec);
0235     }
0236   }
0237
0238   ////////////
0239
0240   // Loop over all PFTaus
0241   math::XYZPointF myleadTkEcalPos;
0242   for (PFTauCollection::size_type iPFTau = 0; iPFTau < pfTaus->size(); iPFTau++) {
0243     PFTauRef thePFTau(pfTaus, iPFTau);
0244     if ((*pfTauIsoDiscr)[thePFTau] == 1) {
0245       if ((*thePFTau).et() > minRecoPt_ && std::abs((*thePFTau).eta()) < maxRecoAbsEta_) {
0246         // Check if track goes to Ecal crack
0247         reco::TrackRef myleadTk;
0248         const reco::PFCandidatePtr& pflch = thePFTau->leadPFChargedHadrCand();
0249         if (pflch.isNonnull()) {
0250           myleadTk = pflch->trackRef();
0251           myleadTkEcalPos = pflch->positionAtECALEntrance();
0252
0253           if (myleadTk.isNonnull()) {
0254             if (applyEcalCrackCut_ && isInEcalCrack(std::abs((double)myleadTkEcalPos.eta()))) {
0255               continue;  // do nothing
0256             } else {
0257               // Match with gen object
0258               for (unsigned int i = 0; i < _GenObjects.size(); i++) {
0259                 if (_GenObjects[i].Et() >= minMCPt_ && std::abs(_GenObjects[i].Eta()) < maxMCAbsEta_) {
0260                   TLorentzVector pftau((*thePFTau).px(), (*thePFTau).py(), (*thePFTau).pz(), (*thePFTau).energy());
0261                   double GenDeltaR = pftau.DeltaR(_GenObjects[i]);
0262                   if (GenDeltaR < maxDeltaR_) {
0263                     hleadTk_pt->Fill((float)myleadTk->pt());
0264                     hleadTk_eta->Fill((float)myleadTk->eta());
0265                     hleadTk_phi->Fill((float)myleadTk->phi());
0266
0267                     hEoverP->Fill((*thePFTau).ecalStripSumEOverPLead());
0268                     hHoverP->Fill((*thePFTau).hcal3x3OverPLead());
0269                     hEmfrac->Fill((*thePFTau).emFraction());
0270
0271                     if (std::abs(myleadTk->eta()) < 1.5) {
0272                       hEoverP_barrel->Fill((*thePFTau).ecalStripSumEOverPLead());
0273                       hHoverP_barrel->Fill((*thePFTau).hcal3x3OverPLead());
0274                       hEmfrac_barrel->Fill((*thePFTau).emFraction());
0275                     } else if (std::abs(myleadTk->eta()) > 1.5 && std::abs(myleadTk->eta()) < 2.5) {
0276                       hEoverP_endcap->Fill((*thePFTau).ecalStripSumEOverPLead());
0277                       hHoverP_endcap->Fill((*thePFTau).hcal3x3OverPLead());
0278                       hEmfrac_endcap->Fill((*thePFTau).emFraction());
0279                     }
0280
0281                     // if -999 fill in -1 bin!
0282                     if ((*thePFTau).electronPreIDOutput() < -1)
0283                       hElecMVA->Fill(-1);
0284                     else
0285                       hElecMVA->Fill((*thePFTau).electronPreIDOutput());
0286
0287                     hTauElecDiscriminant->Fill((*pfTauElecDiscr)[thePFTau]);
0288
0289                     hHoPvsEoP->Fill((*thePFTau).ecalStripSumEOverPLead(), (*thePFTau).hcal3x3OverPLead());
0290                     hEmfracvsEoP->Fill((*thePFTau).emFraction(), (*thePFTau).hcal3x3OverPLead());
0291
0292                     if ((*thePFTau).electronPreIDDecision() == 1) {
0293                       hEoverP_preid1->Fill((*thePFTau).ecalStripSumEOverPLead());
0294                       hHoverP_preid1->Fill((*thePFTau).hcal3x3OverPLead());
0295                       hEmfrac_preid1->Fill((*thePFTau).emFraction());
0296                       hHoPvsEoP_preid1->Fill((*thePFTau).ecalStripSumEOverPLead(), (*thePFTau).hcal3x3OverPLead());
0297                       hEmfracvsEoP_preid1->Fill((*thePFTau).emFraction(), (*thePFTau).hcal3x3OverPLead());
0298                     } else {
0299                       hEoverP_preid0->Fill((*thePFTau).ecalStripSumEOverPLead());
0300                       hHoverP_preid0->Fill((*thePFTau).hcal3x3OverPLead());
0301                       hEmfrac_preid0->Fill((*thePFTau).emFraction());
0302                       hHoPvsEoP_preid0->Fill((*thePFTau).ecalStripSumEOverPLead(), (*thePFTau).hcal3x3OverPLead());
0303                       hEmfracvsEoP_preid0->Fill((*thePFTau).emFraction(), (*thePFTau).hcal3x3OverPLead());
0304                     }
0305                   }
0306                 }
0307
0308                 // Loop over all PFCands for cluster plots
0309                 std::vector<CandidatePtr> myPFCands = (*thePFTau).pfTauTagInfoRef()->PFCands();
0310                 for (int i = 0; i < (int)myPFCands.size(); i++) {
0311                   const reco::PFCandidate* pfCand = dynamic_cast<const reco::PFCandidate*>(myPFCands[i].get());
0312                   if (pfCand == nullptr)
0313                     continue;
0314
0315                   math::XYZPointF candPos;
0316                   if (std::abs(pfCand->pdgId()) == 211 ||
0317                       std::abs(pfCand->pdgId()) == 11)  // if charged hadron or electron
0318                     candPos = pfCand->positionAtECALEntrance();
0319                   else
0320                     candPos = math::XYZPointF(pfCand->px(), pfCand->py(), pfCand->pz());
0321
0322                   //double deltaR   = ROOT::Math::VectorUtil::DeltaR(myleadTkEcalPos,candPos);
0323                   double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myleadTkEcalPos, candPos);
0324                   double deltaEta = std::abs(myleadTkEcalPos.eta() - myPFCands[i]->eta());
0325                   double deltaPhiOverQ = deltaPhi / (double)myleadTk->charge();
0326
0327                   hpfcand_deltaEta->Fill(deltaEta);
0328                   hpfcand_deltaEta_weightE->Fill(deltaEta * pfCand->ecalEnergy());
0329                   hpfcand_deltaPhiOverQ->Fill(deltaPhiOverQ);
0330                   hpfcand_deltaPhiOverQ_weightE->Fill(deltaPhiOverQ * pfCand->ecalEnergy());
0331                 }
0332               }
0333             }
0334           }
0335         }
0336       }
0337     }
0338   }
0339 }
0340
0341 // Ecal crack  map from Egamma POG
0342 bool PFTauElecRejectionBenchmark::isInEcalCrack(double eta) const {
0343   return (eta < 0.018 || (eta > 0.423 && eta < 0.461) || (eta > 0.770 && eta < 0.806) || (eta > 1.127 && eta < 1.163) ||
0344           (eta > 1.460 && eta < 1.558));
0345 }