Reference Manual

Go to the documentation of this file.
00001 #include "RecoParticleFlow/Benchmark/interface/PFTauElecRejectionBenchmark.h"
00002 
00003 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
00004 #define BOOK1D(name,title,nbinsx,lowx,highx) \
00005   h##name = db_ ? db_->book1D(#name,title,nbinsx,lowx,highx)->getTH1F() \
00006     : new TH1F(#name,title,nbinsx,lowx,highx)
00007 
00008 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
00009 #define BOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy) \
00010   h##name = db_ ? db_->book2D(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)->getTH2F() \
00011     : new TH2F(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)
00012 
00013 // all versions OK
00014 // preprocesor macro for setting axis titles
00015 #define SETAXES(name,xtitle,ytitle) \
00016   h##name->GetXaxis()->SetTitle(xtitle); h##name->GetYaxis()->SetTitle(ytitle)
00017 
00018 
00019 using namespace reco;
00020 using namespace std;
00021 
00022 class MonitorElement;
00023 
00024 PFTauElecRejectionBenchmark::PFTauElecRejectionBenchmark() : file_(0) {}
00025 
00026 PFTauElecRejectionBenchmark::~PFTauElecRejectionBenchmark() {
00027   if(file_) file_->Close();
00028 }
00029 
00030 void PFTauElecRejectionBenchmark::write() {
00031    // Store the DAQ Histograms 
00032   if (outputFile_.size() != 0) {
00033     if (db_)
00034           db_->save(outputFile_.c_str());
00035     // use bare Root if no DQM (FWLite applications)
00036     else if (file_) {
00037        file_->Write(outputFile_.c_str());
00038        cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
00039        file_->Close();
00040        }
00041   }
00042   else 
00043     cout << "No output file specified ("<<outputFile_<<"). Results will not be saved!" << endl;
00044     
00045 } 
00046 
00047 void PFTauElecRejectionBenchmark::setup(
00048                                         string Filename,
00049                                         string benchmarkLabel,
00050                                         double maxDeltaR, 
00051                                         double minRecoPt, 
00052                                         double maxRecoAbsEta, 
00053                                         double minMCPt, 
00054                                         double maxMCAbsEta, 
00055                                         string sGenMatchObjectLabel,
00056                                         bool applyEcalCrackCut,
00057                                         DQMStore * db_store) {
00058   maxDeltaR_ = maxDeltaR;
00059   benchmarkLabel_ = benchmarkLabel;
00060   outputFile_=Filename;
00061   minRecoPt_ = minRecoPt;
00062   maxRecoAbsEta_= maxRecoAbsEta;
00063   minMCPt_ = minMCPt;
00064   maxMCAbsEta_= maxMCAbsEta;
00065   sGenMatchObjectLabel_ = sGenMatchObjectLabel;
00066   applyEcalCrackCut_= applyEcalCrackCut;
00067 
00068   file_ = NULL;
00069   db_ = db_store;
00070 
00071   // print parameters
00072   cout<< "PFTauElecRejectionBenchmark Setup parameters =============================================="<<endl;
00073   cout << "Filename to write histograms " << Filename<<endl;
00074   cout << "Benchmark label name " << benchmarkLabel_<<endl;
00075   cout << "maxDeltaRMax " << maxDeltaR << endl;
00076   cout << "minRecoPt " << minRecoPt_ << endl;
00077   cout << "maxRecoAbsEta " << maxRecoAbsEta_ << endl;
00078   cout << "minMCPt " << minMCPt_ << endl;
00079   cout << "maxMCAbsEta " << maxMCAbsEta_ << endl;
00080   cout << "applyEcalCrackCut " << applyEcalCrackCut_ << endl;
00081   cout << "sGenMatchObjectLabel " << sGenMatchObjectLabel_ << endl;
00082 
00083   // Book histogram
00084 
00085   // Establish DQM Store
00086   string path = "PFTask/Benchmarks/"+ benchmarkLabel_ + "/";
00087   path += "Gen";
00088   if (db_) {
00089     db_->setCurrentFolder(path.c_str());
00090   }
00091   else {
00092     file_ = new TFile(outputFile_.c_str(), "recreate");
00093     cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. "<<endl;
00094   }
00095 
00096   // E/p
00097   BOOK1D(EoverP,"E/p",100, 0., 4.);
00098   SETAXES(EoverP, "E/p", "Entries");
00099 
00100   BOOK1D(EoverP_barrel,"E/p barrel",100, 0., 4.);
00101   SETAXES(EoverP_barrel, "E/p barrel", "Entries");
00102 
00103   BOOK1D(EoverP_endcap,"E/p endcap",100, 0., 4.);
00104   SETAXES(EoverP_endcap, "E/p endcap", "Entries");
00105 
00106   BOOK1D(EoverP_preid0,"E/p (preid=0)",100, 0., 4.);
00107   SETAXES(EoverP_preid0, "E/p", "Entries");
00108 
00109   BOOK1D(EoverP_preid1,"E/p (preid=1)",100, 0., 4.);
00110   SETAXES(EoverP_preid1, "E/p", "Entries");
00111 
00112   // H/p
00113   BOOK1D(HoverP,"H/p",100, 0., 2.);
00114   SETAXES(HoverP, "H/p", "Entries");
00115 
00116   BOOK1D(HoverP_barrel,"H/p barrel",100, 0., 2.);
00117   SETAXES(HoverP_barrel, "H/p barrel", "Entries");
00118 
00119   BOOK1D(HoverP_endcap,"H/p endcap",100, 0., 2.);
00120   SETAXES(HoverP_endcap, "H/p endcap", "Entries");
00121 
00122   BOOK1D(HoverP_preid0,"H/p (preid=0)",100, 0., 2.);
00123   SETAXES(HoverP_preid0, "H/p", "Entries");
00124 
00125   BOOK1D(HoverP_preid1,"H/p (preid=1)",100, 0., 2.);
00126   SETAXES(HoverP_preid1, "H/p", "Entries");
00127 
00128   // emfrac
00129   BOOK1D(Emfrac,"EM fraction",100, 0., 1.01);
00130   SETAXES(Emfrac, "em fraction", "Entries");
00131 
00132   BOOK1D(Emfrac_barrel,"EM fraction barrel",100, 0., 1.01);
00133   SETAXES(Emfrac_barrel, "em fraction barrel", "Entries");
00134 
00135   BOOK1D(Emfrac_endcap,"EM fraction endcap",100, 0., 1.01);
00136   SETAXES(Emfrac_endcap, "em fraction endcap", "Entries");
00137 
00138   BOOK1D(Emfrac_preid0,"EM fraction (preid=0)",100, 0., 1.01);
00139   SETAXES(Emfrac_preid0, "em fraction", "Entries");
00140 
00141   BOOK1D(Emfrac_preid1,"EM fraction (preid=1)",100, 0., 1.01);
00142   SETAXES(Emfrac_preid1, "em fraction", "Entries");
00143 
00144   // PreID
00145   BOOK1D(ElecPreID,"PFElectron PreID decision",6, 0., 1.01);
00146   SETAXES(ElecPreID, "PFElectron PreID decision", "Entries");
00147 
00148   // MVA
00149   BOOK1D(ElecMVA,"PFElectron MVA",100, -1.01, 1.01);
00150   SETAXES(ElecMVA, "PFElectron MVA", "Entries");
00151 
00152   // Discriminant
00153   BOOK1D(TauElecDiscriminant,"PFTau-Electron Discriminant",6, 0., 1.01);
00154   SETAXES(TauElecDiscriminant, "PFTau-Electron Discriminant", "Entries");
00155 
00156 
00157   // PFCand clusters
00158   BOOK1D(pfcand_deltaEta,"PFCand cluster dEta",100, 0., 0.8);
00159   SETAXES(pfcand_deltaEta, "PFCand cluster #Delta(#eta)", "Entries");
00160 
00161   BOOK1D(pfcand_deltaEta_weightE,"PFCand cluster dEta, energy weighted",100, 0., 0.8);
00162   SETAXES(pfcand_deltaEta_weightE, "PFCand cluster #Delta(#eta)", "Entries");
00163 
00164   BOOK1D(pfcand_deltaPhiOverQ,"PFCand cluster dPhi/q",100, -0.8, 0.8);
00165   SETAXES(pfcand_deltaPhiOverQ, "PFCand cluster #Delta(#phi)/q", "Entries");
00166 
00167   BOOK1D(pfcand_deltaPhiOverQ_weightE,"PFCand cluster dEta/q, energy weighted",100, -0.8, 0.8);
00168   SETAXES(pfcand_deltaPhiOverQ_weightE, "PFCand cluster #Delta(#phi)/q", "Entries");
00169 
00170 
00171   // Leading KF track
00172   BOOK1D(leadTk_pt,"leading KF track pt",100, 0., 80.);
00173   SETAXES(leadTk_pt, "leading KF track p_{T} (GeV)", "Entries");
00174 
00175   BOOK1D(leadTk_eta,"leading KF track eta",100, -4., 4.);
00176   SETAXES(leadTk_eta, "leading KF track #eta", "Entries");
00177 
00178   BOOK1D(leadTk_phi,"leading KF track phi",100, -3.2, 3.2);
00179   SETAXES(leadTk_phi, "leading KF track #phi", "Entries");
00180 
00181   // Leading Gsf track
00182   BOOK1D(leadGsfTk_pt,"leading Gsf track pt",100, 0., 80.);
00183   SETAXES(leadGsfTk_pt, "leading Gsf track p_{T} (GeV)", "Entries");
00184 
00185   BOOK1D(leadGsfTk_eta,"leading Gsf track eta",100, -4., 4.);
00186   SETAXES(leadGsfTk_eta, "leading Gsf track #eta", "Entries");
00187 
00188   BOOK1D(leadGsfTk_phi,"leading Gsf track phi",100, -3.2, 3.2);
00189   SETAXES(leadGsfTk_phi, "leading Gsf track #phi", "Entries");
00190 
00191 
00192   // H/p vs E/p
00193   BOOK2D(HoPvsEoP,"H/p vs. E/p",100, 0., 2., 100, 0., 2.);
00194   SETAXES(HoPvsEoP, "E/p", "H/p");
00195 
00196   BOOK2D(HoPvsEoP_preid0,"H/p vs. E/p (preid=0)",100, 0., 2., 100, 0., 2.);
00197   SETAXES(HoPvsEoP_preid0, "E/p", "H/p");
00198 
00199   BOOK2D(HoPvsEoP_preid1,"H/p vs. E/p (preid=0)",100, 0., 2., 100, 0., 2.);
00200   SETAXES(HoPvsEoP_preid1, "E/p", "H/p");
00201 
00202   // em fraction vs E/p
00203   BOOK2D(EmfracvsEoP,"emfrac vs. E/p",100, 0., 2., 100, 0., 1.01);
00204   SETAXES(EmfracvsEoP, "E/p", "em fraction");
00205 
00206   BOOK2D(EmfracvsEoP_preid0,"emfrac vs. E/p (preid=0)",100, 0., 2., 100, 0., 1.01);
00207   SETAXES(EmfracvsEoP_preid0, "E/p", "em fraction");
00208 
00209   BOOK2D(EmfracvsEoP_preid1,"emfrac vs. E/p (preid=0)",100, 0., 2., 100, 0., 1.01);
00210   SETAXES(EmfracvsEoP_preid1, "E/p", "em fraction");
00211 
00212 
00213 }
00214 
00215 
00216 void PFTauElecRejectionBenchmark::process(edm::Handle<edm::HepMCProduct> mcevt, edm::Handle<reco::PFTauCollection> pfTaus, 
00217                                           edm::Handle<reco::PFTauDiscriminator> pfTauIsoDiscr, 
00218                                           edm::Handle<reco::PFTauDiscriminator> pfTauElecDiscr) {
00219 
00220 
00221   // Find Gen Objects to be matched with
00222   HepMC::GenEvent * generated_event = new HepMC::GenEvent(*(mcevt->GetEvent()));
00223   _GenObjects.clear();
00224   
00225   TLorentzVector taunet;
00226   HepMC::GenEvent::particle_iterator p;
00227   for (p = generated_event->particles_begin(); p != generated_event->particles_end(); p++) {
00228     if(std::abs((*p)->pdg_id()) == 15&&(*p)->status()==2) { 
00229       bool lept_decay = false;     
00230       TLorentzVector tau((*p)->momentum().px(),(*p)->momentum().py(),(*p)->momentum().pz(),(*p)->momentum().e());
00231       HepMC::GenVertex::particle_iterator z = (*p)->end_vertex()->particles_begin(HepMC::descendants);
00232       for(; z != (*p)->end_vertex()->particles_end(HepMC::descendants); z++) {
00233         if(std::abs((*z)->pdg_id()) == 11 || std::abs((*z)->pdg_id()) == 13) lept_decay=true;
00234         if(std::abs((*z)->pdg_id()) == 16)
00235           taunet.SetPxPyPzE((*z)->momentum().px(),(*z)->momentum().py(),(*z)->momentum().pz(),(*z)->momentum().e());
00236         
00237       }
00238       if(lept_decay==false) {
00239         TLorentzVector jetMom=tau-taunet;
00240         if (sGenMatchObjectLabel_=="tau") _GenObjects.push_back(jetMom);
00241       }
00242     } else if(std::abs((*p)->pdg_id()) == 11&&(*p)->status()==1) { 
00243       TLorentzVector elec((*p)->momentum().px(),(*p)->momentum().py(),(*p)->momentum().pz(),(*p)->momentum().e());
00244       if (sGenMatchObjectLabel_=="e") _GenObjects.push_back(elec);
00245     } 
00246   }
00247  
00249 
00250   
00251   // Loop over all PFTaus
00252   math::XYZPointF myleadTkEcalPos;
00253   for (PFTauCollection::size_type iPFTau=0;iPFTau<pfTaus->size();iPFTau++) { 
00254     PFTauRef thePFTau(pfTaus,iPFTau); 
00255     if ((*pfTauIsoDiscr)[thePFTau] == 1) {
00256       if ((*thePFTau).et() > minRecoPt_ && std::abs((*thePFTau).eta()) < maxRecoAbsEta_) {
00257 
00258         // Check if track goes to Ecal crack
00259         TrackRef myleadTk;
00260         if(thePFTau->leadPFChargedHadrCand().isNonnull()){
00261           myleadTk=thePFTau->leadPFChargedHadrCand()->trackRef();
00262           myleadTkEcalPos = thePFTau->leadPFChargedHadrCand()->positionAtECALEntrance();
00263           
00264           if(myleadTk.isNonnull()){ 
00265             if (applyEcalCrackCut_ && isInEcalCrack(std::abs((double)myleadTkEcalPos.eta()))) {
00266               continue; // do nothing
00267             } else {
00268 
00269               // Match with gen object
00270               for (unsigned int i = 0; i<_GenObjects.size();i++) {
00271                 if (_GenObjects[i].Et() >= minMCPt_ && std::abs(_GenObjects[i].Eta()) < maxMCAbsEta_ ) {
00272                   TLorentzVector pftau((*thePFTau).px(),(*thePFTau).py(),(*thePFTau).pz(),(*thePFTau).energy());
00273                   double GenDeltaR = pftau.DeltaR(_GenObjects[i]);
00274                   if (GenDeltaR<maxDeltaR_) {
00275                     
00276                     hleadTk_pt->Fill((float)myleadTk->pt());
00277                     hleadTk_eta->Fill((float)myleadTk->eta());
00278                     hleadTk_phi->Fill((float)myleadTk->phi());
00279 
00280                     hEoverP->Fill((*thePFTau).ecalStripSumEOverPLead());
00281                     hHoverP->Fill((*thePFTau).hcal3x3OverPLead());
00282                     hEmfrac->Fill((*thePFTau).emFraction());
00283 
00284                     if (std::abs(myleadTk->eta())<1.5) {
00285                       hEoverP_barrel->Fill((*thePFTau).ecalStripSumEOverPLead());
00286                       hHoverP_barrel->Fill((*thePFTau).hcal3x3OverPLead());
00287                       hEmfrac_barrel->Fill((*thePFTau).emFraction());
00288                     } else if (std::abs(myleadTk->eta())>1.5 && std::abs(myleadTk->eta())<2.5) {
00289                       hEoverP_endcap->Fill((*thePFTau).ecalStripSumEOverPLead());
00290                       hHoverP_endcap->Fill((*thePFTau).hcal3x3OverPLead());
00291                       hEmfrac_endcap->Fill((*thePFTau).emFraction());
00292                     }
00293 
00294                     // if -999 fill in -1 bin!
00295                     if ((*thePFTau).electronPreIDOutput()<-1)
00296                       hElecMVA->Fill(-1);
00297                     else
00298                       hElecMVA->Fill((*thePFTau).electronPreIDOutput());
00299 
00300                     hTauElecDiscriminant->Fill((*pfTauElecDiscr)[thePFTau]);
00301 
00302                     hHoPvsEoP->Fill((*thePFTau).ecalStripSumEOverPLead(),(*thePFTau).hcal3x3OverPLead());
00303                     hEmfracvsEoP->Fill((*thePFTau).emFraction(),(*thePFTau).hcal3x3OverPLead());
00304 
00305                     if ((*thePFTau).electronPreIDDecision()==1) {
00306                       hEoverP_preid1->Fill((*thePFTau).ecalStripSumEOverPLead());
00307                       hHoverP_preid1->Fill((*thePFTau).hcal3x3OverPLead());
00308                       hEmfrac_preid1->Fill((*thePFTau).emFraction());
00309                       hHoPvsEoP_preid1->Fill((*thePFTau).ecalStripSumEOverPLead(),(*thePFTau).hcal3x3OverPLead());
00310                       hEmfracvsEoP_preid1->Fill((*thePFTau).emFraction(),(*thePFTau).hcal3x3OverPLead());
00311                     } else {
00312                       hEoverP_preid0->Fill((*thePFTau).ecalStripSumEOverPLead());
00313                       hHoverP_preid0->Fill((*thePFTau).hcal3x3OverPLead());
00314                       hEmfrac_preid0->Fill((*thePFTau).emFraction());
00315                       hHoPvsEoP_preid0->Fill((*thePFTau).ecalStripSumEOverPLead(),(*thePFTau).hcal3x3OverPLead());
00316                       hEmfracvsEoP_preid0->Fill((*thePFTau).emFraction(),(*thePFTau).hcal3x3OverPLead());
00317                     }
00318 
00319 
00320                   }
00321                 }
00322 
00323                 // Loop over all PFCands for cluster plots  
00324                 PFCandidateRefVector myPFCands=(*thePFTau).pfTauTagInfoRef()->PFCands();
00325                 for(int i=0;i<(int)myPFCands.size();i++){
00326 
00327                   math::XYZPointF candPos;
00328                   if (myPFCands[i]->particleId()==1 || myPFCands[i]->particleId()==2) // if charged hadron or electron
00329                     candPos = myPFCands[i]->positionAtECALEntrance();
00330                   else
00331                     candPos = math::XYZPointF(myPFCands[i]->px(),myPFCands[i]->py(),myPFCands[i]->pz());
00332 
00333                   //double deltaR   = ROOT::Math::VectorUtil::DeltaR(myleadTkEcalPos,candPos);
00334                   double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myleadTkEcalPos,candPos);
00335                   double deltaEta = std::abs(myleadTkEcalPos.eta()-myPFCands[i]->eta());
00336                   double deltaPhiOverQ = deltaPhi/(double)myleadTk->charge();
00337                   
00338                   hpfcand_deltaEta->Fill(deltaEta);
00339                   hpfcand_deltaEta_weightE->Fill(deltaEta*myPFCands[i]->ecalEnergy());
00340                   hpfcand_deltaPhiOverQ->Fill(deltaPhiOverQ);
00341                   hpfcand_deltaPhiOverQ_weightE->Fill(deltaPhiOverQ*myPFCands[i]->ecalEnergy());        
00342                   
00343                 }
00344 
00345               }
00346 
00347             }
00348           }
00349         }
00350 
00351       }
00352     }
00353   }
00354 }
00355 
00356 // Ecal crack  map from Egamma POG
00357 bool PFTauElecRejectionBenchmark::isInEcalCrack(double eta) const{  
00358   return (eta < 0.018 || 
00359           (eta>0.423 && eta<0.461) ||
00360           (eta>0.770 && eta<0.806) ||
00361           (eta>1.127 && eta<1.163) ||
00362           (eta>1.460 && eta<1.558));
00363 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoParticleFlow/Benchmark/src/PFTauElecRejectionBenchmark.cc
