CMSSW Reference Manual

00001 #include "RecoParticleFlow/Benchmark/interface/GenericBenchmark.h"
00002 #include "DQMServices/Core/interface/MonitorElement.h"
00003 
00004 // CMSSW_2_X_X
00005 #include "DQMServices/Core/interface/DQMStore.h"
00006 
00007 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
00008 #define BOOK1D(name,title,nbinsx,lowx,highx) \
00009   h##name = DQM ? DQM->book1D(#name,title,nbinsx,lowx,highx)->getTH1F() \
00010     : new TH1F(#name,title,nbinsx,lowx,highx)
00011 
00012 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
00013 #define BOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy) \
00014   h##name = DQM ? DQM->book2D(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)->getTH2F() \
00015     : new TH2F(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)
00016 
00017 
00018 
00019 
00020 
00021 
00022 // all versions OK
00023 // preprocesor macro for setting axis titles
00024 
00025 #define SETAXES(name,xtitle,ytitle) \
00026   h##name->GetXaxis()->SetTitle(xtitle); h##name->GetYaxis()->SetTitle(ytitle)
00027 
00028 #define ET (PlotAgainstReco_)?"reconstructed E_{T}":"generated E_{T}"
00029 #define ETA (PlotAgainstReco_)?"reconstructed #eta":"generated #eta"
00030 #define PHI (PlotAgainstReco_)?"reconstructed #phi":"generated #phi"
00031 
00032 
00033 
00034 GenericBenchmark::GenericBenchmark() {}
00035 
00036 GenericBenchmark::~GenericBenchmark() {}
00037 
00038 void GenericBenchmark::setup(DQMStore *DQM, bool PlotAgainstReco_) { 
00039 
00040   // CMSSW_2_X_X
00041   // use bare Root if no DQM (FWLite applications)
00042   if (!DQM) file_ = new TFile();
00043 
00044   // Book Histograms
00045 
00046   // delta et quantities
00047   BOOK1D(DeltaEt,"#DeltaE_{T}",1000,-60,40);
00048   BOOK2D(DeltaEtvsEt,"#DeltaE_{T} vs E_{T}",1000,0,1000,1000,-100,100);
00049   BOOK2D(DeltaEtOverEtvsEt,"#DeltaE_{T}/E_{T} vsE_{T}",1000,0,1000,100,-1,1);
00050   BOOK2D(DeltaEtvsEta,"#DeltaE_{T} vs #eta",200,-5,5,1000,-100,100);
00051   BOOK2D(DeltaEtOverEtvsEta,"#DeltaE_{T}/E_{T} vs #eta",200,-5,5,100,-1,1);
00052   BOOK2D(DeltaEtvsPhi,"#DeltaE_{T} vs #phi",200,-M_PI,M_PI,1000,-100,100);
00053   BOOK2D(DeltaEtOverEtvsPhi,"#DeltaE_{T}/E_{T} vs #Phi",200,-M_PI,M_PI,100,-1,1);
00054   BOOK2D(DeltaEtvsDeltaR,"#DeltaE_{T} vs #DeltaR",100,0,1,1000,-100,100);
00055   BOOK2D(DeltaEtOverEtvsDeltaR,"#DeltaE_{T}/E_{T} vs #DeltaR",100,0,1,100,-1,1);
00056 
00057   // delta eta quantities
00058   BOOK1D(DeltaEta,"#Delta#eta",100,-0.2,0.2);
00059   BOOK2D(DeltaEtavsEt,"#Delta#eta vs E_{T}",1000,0,1000,100,-3,3);
00060   BOOK2D(DeltaEtaOverEtavsEt,"#Delta#eta/#eta vs E_(T}",1000,0,1000,100,-1,1); // ms: propose remove
00061   BOOK2D(DeltaEtavsEta,"#Delta#eta vs #eta",200,-5,5,100,-3,3);
00062   BOOK2D(DeltaEtaOverEtavsEta,"EDelta#eta/#eta vs #eta",200,-5,5,100,-1,1); // ms: propose remove
00063   BOOK2D(DeltaEtavsPhi,"#Delta#eta vs #phi",200,-M_PI,M_PI,200,-M_PI,M_PI); // ms: propose remove
00064   BOOK2D(DeltaEtaOverEtavsPhi,"#Delta#eta/#eta vs #phi",200,-M_PI,M_PI,100,-1,1); // ms: propose remove
00065 
00066   // delta phi quantities
00067   BOOK1D(DeltaPhi,"#Delta#phi",100,-0.2,0.2);
00068   BOOK2D(DeltaPhivsEt,"#Delta#phi vs E_{T}",1000,0,1000,100,-M_PI_2,M_PI_2);
00069   BOOK2D(DeltaPhiOverPhivsEt,"#Delta#phi/#phi vs E_{T}",1000,0,1000,100,-1,1); // ms: propose remove
00070   BOOK2D(DeltaPhivsEta,"#Delta#phi vs #eta",200,-5,5,100,-M_PI_2,M_PI_2);
00071   BOOK2D(DeltaPhiOverPhivsEta,"#Delta#phi/#phi vs #eta",200,-5,5,100,-1,1); // ms: propose remove
00072   BOOK2D(DeltaPhivsPhi,"#Delta#phi vs #phi",200,-M_PI,M_PI,200,-M_PI,M_PI); // ms: propose remove
00073   BOOK2D(DeltaPhiOverPhivsPhi,"#Delta#phi/#phi vs #phi",200,-M_PI,M_PI,100,-1,1); // ms: propose remove
00074 
00075   // delta R quantities
00076   BOOK1D(DeltaR,"#DeltaR",100,0,1);
00077   BOOK2D(DeltaRvsEt,"#DeltaR vs E_{T}",1000,0,1000,100,0,1);
00078   BOOK2D(DeltaRvsEta,"#DeltaR vs #eta",200,-5,5,100,0,1);
00079   BOOK2D(DeltaRvsPhi,"#DeltaR vs #phi",200,-M_PI,M_PI,100,0,1); // ms: propose remove
00080 
00081   // number of truth particles found within given cone radius of reco
00082   //BOOK2D(NumInConeVsConeSize,NumInConeVsConeSize,100,0,1,25,0,25);
00083 
00084   // Set Axis Titles
00085  
00086   // delta et quantities
00087   SETAXES(DeltaEt,"#DeltaE_{T}","Events");
00088   SETAXES(DeltaEtvsEt,ET,"#DeltaE_{T}");
00089   SETAXES(DeltaEtOverEtvsEt,ET,"#DeltaE_{T}/E_{T}");
00090   SETAXES(DeltaEtvsEta,ETA,"#DeltaE_{T}");
00091   SETAXES(DeltaEtOverEtvsEta,ETA,"#DeltaE_{T}/E_{T}");
00092   SETAXES(DeltaEtvsPhi,PHI,"#DeltaE_{T}");
00093   SETAXES(DeltaEtOverEtvsPhi,PHI,"#DeltaE_{T}/E_{T}");
00094   SETAXES(DeltaEtvsDeltaR,"#DeltaR","#DeltaE_{T}");
00095   SETAXES(DeltaEtOverEtvsDeltaR,"#DeltaR","#DeltaE_{T}/E_{T}");
00096 
00097   // delta eta quantities
00098   SETAXES(DeltaEta,"#Delta#eta","Events");
00099   SETAXES(DeltaEtavsEt,ET,"#Delta#eta");
00100   SETAXES(DeltaEtavsEta,ETA,"#Delta#eta");
00101   SETAXES(DeltaEtaOverEtavsEt,ET,"#Delta#eta/#eta");
00102   SETAXES(DeltaEtaOverEtavsEta,ETA,"#Delta#eta/#eta");
00103   SETAXES(DeltaEtavsPhi,PHI,"#Delta#eta");
00104   SETAXES(DeltaEtaOverEtavsPhi,PHI,"#Delta#eta/#eta");
00105 
00106   // delta phi quantities
00107   SETAXES(DeltaPhi,"#Delta#phi","Events");
00108   SETAXES(DeltaPhivsEt,ET,"#Delta#phi");
00109   SETAXES(DeltaPhivsEta,ETA,"#Delta#phi");
00110   SETAXES(DeltaPhiOverPhivsEt,ET,"#Delta#phi/#phi");
00111   SETAXES(DeltaPhiOverPhivsEta,ETA,"#Delta#phi/#phi");
00112   SETAXES(DeltaPhivsPhi,PHI,"#Delta#phi");
00113   SETAXES(DeltaPhiOverPhivsPhi,PHI,"#Delta#phi/#phi");
00114 
00115   // delta R quantities
00116   SETAXES(DeltaR,"#DeltaR","Events");
00117   SETAXES(DeltaRvsEt,ET,"#DeltaR");
00118   SETAXES(DeltaRvsEta,ETA,"#DeltaR");
00119   SETAXES(DeltaRvsPhi,PHI,"#DeltaR");
00120   
00121 }
00122 
00123 
00124 void GenericBenchmark::fill(const edm::View<reco::Candidate> *RecoCollection, const edm::View<reco::Candidate> *GenCollection, bool PlotAgainstReco, bool onlyTwoJets, double recPt_cut, double maxEta_cut, double deltaR_cut) {
00125 
00126   // loop over reco particles
00127   for (unsigned int i = 0; i < RecoCollection->size(); i++) {
00128     
00129     // generate histograms comparing the reco and truth candidate (truth = closest in delta-R)
00130     const reco::Candidate *particle = &(*RecoCollection)[i];
00131     const reco::Candidate *gen_particle = algo_->matchByDeltaR(particle,GenCollection);
00132 
00133    if  ((particle!=NULL) &&  (gen_particle!=NULL)){
00134       //   std::cout << RecoCollection->size() << " " << GenCollection->size() <<particle  <<gen_particle << std::endl;   
00135 
00136     // Count the number of jets with a larger energy
00137     unsigned highJets = 0;
00138     for(unsigned j=0; j<RecoCollection->size(); j++) { 
00139       const reco::Candidate *otherParticle = &(*RecoCollection)[j];
00140       if ( j != i && otherParticle->pt() > particle->pt() ) highJets++;
00141     }
00142     if ( onlyTwoJets && highJets > 1 ) continue;
00143                 
00144     //skip reconstructed PFJets with p_t < recPt_cut
00145     if (particle->pt() < recPt_cut and recPt_cut != -1.)
00146       continue;
00147     //skip if PFJet within eta>maxEta_cut
00148     if (fabs(particle->eta())>maxEta_cut and maxEta_cut != -1.)
00149       continue;
00150 
00151     // get the quantities to place on the denominator and/or divide by
00152     double et, eta, phi;
00153     if (PlotAgainstReco) { et = particle->et(); eta = particle->eta(); phi = particle->phi(); }
00154     else { et = gen_particle->et(); eta = gen_particle->eta(); phi = gen_particle->phi(); }
00155     
00156     // get the delta quantities
00157     double deltaEt = algo_->deltaEt(particle,gen_particle);
00158     double deltaR = algo_->deltaR(particle,gen_particle);
00159     double deltaEta = algo_->deltaEta(particle,gen_particle);
00160     double deltaPhi = algo_->deltaPhi(particle,gen_particle);
00161    
00162     //TODO implement variable Cut:
00163      if (fabs(deltaR)>deltaR_cut and deltaR_cut != -1.)
00164        continue;
00165 
00166     // fill histograms
00167     hDeltaEt->Fill(deltaEt);
00168     hDeltaEtvsEt->Fill(et,deltaEt);
00169     hDeltaEtOverEtvsEt->Fill(et,deltaEt/et);
00170     hDeltaEtvsEta->Fill(eta,deltaEt);
00171     hDeltaEtOverEtvsEta->Fill(eta,deltaEt/et);
00172     hDeltaEtvsPhi->Fill(phi,deltaEt);
00173     hDeltaEtOverEtvsPhi->Fill(phi,deltaEt/et);
00174     hDeltaEtvsDeltaR->Fill(deltaR,deltaEt);
00175     hDeltaEtOverEtvsDeltaR->Fill(deltaR,deltaEt/et);
00176     
00177     hDeltaEta->Fill(deltaEta);
00178     hDeltaEtavsEt->Fill(et,deltaEta/eta);
00179     hDeltaEtaOverEtavsEt->Fill(et,deltaEta/eta);
00180     hDeltaEtavsEta->Fill(eta,deltaEta);
00181     hDeltaEtaOverEtavsEta->Fill(eta,deltaEta/eta);
00182     hDeltaEtavsPhi->Fill(phi,deltaEta);
00183     hDeltaEtaOverEtavsPhi->Fill(phi,deltaEta/eta);
00184     
00185     hDeltaPhi->Fill(deltaPhi);
00186     hDeltaPhivsEt->Fill(et,deltaPhi);
00187     hDeltaPhiOverPhivsEt->Fill(et,deltaPhi/phi);
00188     hDeltaPhivsEta->Fill(eta,deltaPhi);
00189     hDeltaPhiOverPhivsEta->Fill(eta,deltaPhi/phi);
00190     hDeltaPhivsPhi->Fill(phi,deltaPhi);
00191     hDeltaPhiOverPhivsPhi->Fill(phi,deltaPhi/phi);
00192 
00193     hDeltaR->Fill(deltaR);
00194     hDeltaRvsEt->Fill(et,deltaR);
00195     hDeltaRvsEta->Fill(eta,deltaR);
00196    }
00197   }
00198 
00199 }
00200 
00201 void GenericBenchmark::write(std::string Filename) {
00202 
00203   if (Filename.size() != 0 && file_)
00204     file_->Write(Filename.c_str());
00205 
00206 }
GenericBenchmark.cc
