Reference Manual

Go to the documentation of this file.
00001 #include "RecoParticleFlow/Benchmark/interface/GenericBenchmark.h"
00002 #include "RecoParticleFlow/Benchmark/interface/BenchmarkTree.h"
00003 #include "DQMServices/Core/interface/MonitorElement.h"
00004 
00005 // CMSSW_2_X_X
00006 #include "DQMServices/Core/interface/DQMStore.h"
00007 
00008 #include <TROOT.h>
00009 #include <TFile.h>
00010 
00011 //Colin: it seems a bit strange to use the preprocessor for that kind of 
00012 //thing. Looks like all these macros could be replaced by plain functions.
00013 
00014 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
00015 #define BOOK1D(name,title,nbinsx,lowx,highx)                            \
00016   h##name = DQM ? DQM->book1D(#name,title,nbinsx,lowx,highx)->getTH1F() \
00017     : new TH1F(#name,title,nbinsx,lowx,highx)
00018 
00019 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
00020 #define BOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)          \
00021   h##name = DQM ? DQM->book2D(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)->getTH2F() \
00022     : new TH2F(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)
00023 
00024 
00025 // all versions OK
00026 // preprocesor macro for setting axis titles
00027 
00028 #define SETAXES(name,xtitle,ytitle)                                     \
00029   h##name->GetXaxis()->SetTitle(xtitle); h##name->GetYaxis()->SetTitle(ytitle)
00030 
00031 #define ET (PlotAgainstReco_)?"reconstructed E_{T} [GeV]":"generated E_{T} [GeV]"
00032 #define ETA (PlotAgainstReco_)?"reconstructed #eta":"generated #eta"
00033 #define PHI (PlotAgainstReco_)?"reconstructed #phi (rad)":"generated #phi (rad)"
00034 
00035 using namespace std;
00036 
00037 GenericBenchmark::GenericBenchmark() {}
00038 
00039 GenericBenchmark::~GenericBenchmark() {}
00040 
00041 void GenericBenchmark::setup(DQMStore *DQM, bool PlotAgainstReco_, float minDeltaEt, float maxDeltaEt,
00042                              float minDeltaPhi, float maxDeltaPhi, bool doMetPlots) {
00043 
00044   //std::cout << "minDeltaPhi = " << minDeltaPhi << std::endl;
00045 
00046   // CMSSW_2_X_X
00047   // use bare Root if no DQM (FWLite applications)
00048   //if (!DQM)
00049   //{
00050   //  file_ = new TFile("pfmetBenchmark.root", "recreate");
00051   //  cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. "<<endl;
00052   //}
00053   // Book Histograms
00054 
00055   //std::cout << "FL : pwd = ";
00056   //gDirectory->pwd();
00057   //std::cout << std::endl;
00058 
00059   int nbinsEt = 1000;
00060   float minEt = 0;
00061   float maxEt = 1000;
00062 
00063   //float minDeltaEt = -100;
00064   //float maxDeltaEt = 50;
00065 
00066   int nbinsEta = 200;
00067   float minEta = -5;
00068   float maxEta = 5;
00069 
00070   int nbinsDeltaEta = 1000;
00071   float minDeltaEta = -0.5;
00072   float maxDeltaEta = 0.5;
00073 
00074   int nbinsDeltaPhi = 1000;
00075   //float minDeltaPhi = -0.5;
00076   //float maxDeltaPhi = 0.5;
00077 
00078 
00079   // delta et quantities
00080   BOOK1D(DeltaEt,"#DeltaE_{T}", nbinsEt, minDeltaEt, maxDeltaEt);
00081   BOOK1D(DeltaEx,"#DeltaE_{X}", nbinsEt, minDeltaEt, maxDeltaEt);
00082   BOOK1D(DeltaEy,"#DeltaE_{Y}", nbinsEt, minDeltaEt, maxDeltaEt);
00083   BOOK2D(DeltaEtvsEt,"#DeltaE_{T} vs E_{T}",
00084          nbinsEt, minEt, maxEt,
00085          nbinsEt,minDeltaEt, maxDeltaEt);
00086   BOOK2D(DeltaEtOverEtvsEt,"#DeltaE_{T}/E_{T} vsE_{T}",
00087          nbinsEt, minEt, maxEt,
00088          nbinsEt,-1,1);
00089   BOOK2D(DeltaEtvsEta,"#DeltaE_{T} vs #eta",
00090          nbinsEta, minEta, maxEta,
00091          nbinsEt,minDeltaEt, maxDeltaEt);
00092   BOOK2D(DeltaEtOverEtvsEta,"#DeltaE_{T}/E_{T} vs #eta",
00093          nbinsEta, minEta, maxEta,
00094          100,-1,1);
00095   BOOK2D(DeltaEtvsPhi,"#DeltaE_{T} vs #phi",
00096          200,-M_PI,M_PI,
00097          nbinsEt,minDeltaEt, maxDeltaEt);
00098   BOOK2D(DeltaEtOverEtvsPhi,"#DeltaE_{T}/E_{T} vs #Phi",
00099          200,-M_PI,M_PI,
00100          100,-1,1);
00101   BOOK2D(DeltaEtvsDeltaR,"#DeltaE_{T} vs #DeltaR",
00102          100,0,1,
00103          nbinsEt,minDeltaEt, maxDeltaEt);
00104   BOOK2D(DeltaEtOverEtvsDeltaR,"#DeltaE_{T}/E_{T} vs #DeltaR",
00105          100,0,1,
00106          100,-1,1);
00107 
00108   // delta eta quantities
00109   BOOK1D(DeltaEta,"#Delta#eta",nbinsDeltaEta,minDeltaEta,maxDeltaEta);
00110   BOOK2D(DeltaEtavsEt,"#Delta#eta vs E_{T}",
00111          nbinsEt, minEt, maxEt,
00112          nbinsDeltaEta,minDeltaEta,maxDeltaEta);
00113   BOOK2D(DeltaEtavsEta,"#Delta#eta vs #eta",
00114          nbinsEta, minEta, maxEta,
00115          nbinsDeltaEta,minDeltaEta,maxDeltaEta);
00116 
00117   // delta phi quantities
00118   BOOK1D(DeltaPhi,"#Delta#phi",nbinsDeltaPhi,minDeltaPhi,maxDeltaPhi);
00119   BOOK2D(DeltaPhivsEt,"#Delta#phi vs E_{T}",
00120          nbinsEt, minEt, maxEt,
00121          nbinsDeltaPhi,minDeltaPhi,maxDeltaPhi);
00122   BOOK2D(DeltaPhivsEta,"#Delta#phi vs #eta",
00123          nbinsEta, minEta, maxEta,
00124          nbinsDeltaPhi,minDeltaPhi,maxDeltaPhi);
00125 
00126   // delta R quantities
00127   BOOK1D(DeltaR,"#DeltaR",100,0,1);
00128   BOOK2D(DeltaRvsEt,"#DeltaR vs E_{T}",
00129          nbinsEt, minEt, maxEt,
00130          100,0,1);
00131   BOOK2D(DeltaRvsEta,"#DeltaR vs #eta",
00132          nbinsEta, minEta, maxEta,
00133          100,0,1);
00134 
00135   BOOK1D(NRec,"Number of reconstructed objects",20,0,20);
00136 
00137   // seen and gen distributions, for efficiency computation
00138   BOOK1D(EtaSeen,"seen #eta",100,-5,5);
00139   BOOK1D(PhiSeen,"seen #phi",100,-3.5,3.5);
00140   BOOK1D(EtSeen,"seen E_{T}",nbinsEt, minEt, maxEt);
00141   BOOK2D(EtvsEtaSeen,"seen E_{T} vs eta",100,-5,5,200, 0, 200);
00142   BOOK2D(EtvsPhiSeen,"seen E_{T} vs seen #phi",100,-3.5,3.5,200, 0, 200);  
00143 
00144   BOOK1D(PhiRec,"Rec #phi",100,-3.5,3.5);
00145   BOOK1D(EtRec,"Rec E_{T}",nbinsEt, minEt, maxEt);
00146   BOOK1D(ExRec,"Rec E_{X}",nbinsEt, -maxEt, maxEt);
00147   BOOK1D(EyRec,"Rec E_{Y}",nbinsEt, -maxEt, maxEt);
00148 
00149   BOOK2D(EtRecvsEt,"Rec E_{T} vs E_{T}",
00150          nbinsEt, minEt, maxEt,
00151          nbinsEt, minEt, maxEt);
00152   BOOK2D(EtRecOverTrueEtvsTrueEt,"Rec E_{T} / E_{T} vs E_{T}",
00153          nbinsEt, minEt, maxEt,
00154          1000, 0., 100.);
00155 
00156   BOOK1D(EtaGen,"generated #eta",100,-5,5);
00157   BOOK1D(PhiGen,"generated #phi",100,-3.5,3.5);
00158   BOOK1D(EtGen,"generated E_{T}",nbinsEt, minEt, maxEt);
00159   BOOK2D(EtvsEtaGen,"generated E_{T} vs generated #eta",100,-5,5,200, 0, 200);  
00160   BOOK2D(EtvsPhiGen,"generated E_{T} vs generated #phi",100,-3.5,3.5, 200, 0, 200);  
00161 
00162   BOOK1D(NGen,"Number of generated objects",20,0,20);
00163 
00164   if (doMetPlots)
00165   {
00166     BOOK1D(SumEt,"SumEt", 1000, 0., 3000.);
00167     BOOK1D(TrueSumEt,"TrueSumEt", 1000, 0., 3000.);
00168     BOOK2D(DeltaSetvsSet,"#DeltaSEt vs trueSEt",
00169            3000, 0., 3000.,
00170            1000,-1000., 1000.);
00171     BOOK2D(DeltaMexvsSet,"#DeltaMEX vs trueSEt",
00172            3000, 0., 3000.,
00173            1000,-400., 400.);
00174     BOOK2D(DeltaSetOverSetvsSet,"#DeltaSetOverSet vs trueSet",
00175            3000, 0., 3000.,
00176            1000,-1., 1.);
00177     BOOK2D(RecSetvsTrueSet,"Set vs trueSet",
00178            3000, 0., 3000.,
00179            3000,0., 3000.);
00180     BOOK2D(RecSetOverTrueSetvsTrueSet,"Set/trueSet vs trueSet",
00181            3000, 0., 3000.,
00182            500,0., 2.);
00183     BOOK2D(TrueMexvsTrueSet,"trueMex vs trueSet",
00184            3000,0., 3000.,
00185            nbinsEt, -maxEt, maxEt);
00186   }
00187   // number of truth particles found within given cone radius of reco
00188   //BOOK2D(NumInConeVsConeSize,NumInConeVsConeSize,100,0,1,25,0,25);
00189 
00190   // Set Axis Titles
00191  
00192   // delta et quantities
00193   SETAXES(DeltaEt,"#DeltaE_{T} [GeV]","");
00194   SETAXES(DeltaEx,"#DeltaE_{X} [GeV]","");
00195   SETAXES(DeltaEy,"#DeltaE_{Y} [GeV]","");
00196   SETAXES(DeltaEtvsEt,ET,"#DeltaE_{T} [GeV]");
00197   SETAXES(DeltaEtOverEtvsEt,ET,"#DeltaE_{T}/E_{T}");
00198   SETAXES(DeltaEtvsEta,ETA,"#DeltaE_{T} [GeV]");
00199   SETAXES(DeltaEtOverEtvsEta,ETA,"#DeltaE_{T}/E_{T}");
00200   SETAXES(DeltaEtvsPhi,PHI,"#DeltaE_{T} [GeV]");
00201   SETAXES(DeltaEtOverEtvsPhi,PHI,"#DeltaE_{T}/E_{T}");
00202   SETAXES(DeltaEtvsDeltaR,"#DeltaR","#DeltaE_{T} [GeV]");
00203   SETAXES(DeltaEtOverEtvsDeltaR,"#DeltaR","#DeltaE_{T}/E_{T}");
00204 
00205   // delta eta quantities
00206   SETAXES(DeltaEta,"#Delta#eta","Events");
00207   SETAXES(DeltaEtavsEt,ET,"#Delta#eta");
00208   SETAXES(DeltaEtavsEta,ETA,"#Delta#eta");
00209 
00210   // delta phi quantities
00211   SETAXES(DeltaPhi,"#Delta#phi [rad]","Events");
00212   SETAXES(DeltaPhivsEt,ET,"#Delta#phi [rad]");
00213   SETAXES(DeltaPhivsEta,ETA,"#Delta#phi [rad]");
00214 
00215   // delta R quantities
00216   SETAXES(DeltaR,"#DeltaR","Events");
00217   SETAXES(DeltaRvsEt,ET,"#DeltaR");
00218   SETAXES(DeltaRvsEta,ETA,"#DeltaR");
00219 
00220   SETAXES(NRec,"Number of Rec Objects","");
00221   
00222   SETAXES(EtaSeen,"seen #eta","");
00223   SETAXES(PhiSeen,"seen #phi [rad]","");
00224   SETAXES(EtSeen,"seen E_{T} [GeV]","");
00225   SETAXES(EtvsEtaSeen,"seen #eta","seen E_{T}");
00226   SETAXES(EtvsPhiSeen,"seen #phi [rad]","seen E_{T}");
00227 
00228   SETAXES(PhiRec,"#phi [rad]","");
00229   SETAXES(EtRec,"E_{T} [GeV]","");
00230   SETAXES(ExRec,"E_{X} [GeV]","");
00231   SETAXES(EyRec,"E_{Y} [GeV]","");
00232   SETAXES(EtRecvsEt,ET,"Rec E_{T} [GeV]");
00233   SETAXES(EtRecOverTrueEtvsTrueEt,ET,"Rec E_{T} / E_{T} [GeV]");
00234 
00235   SETAXES(EtaGen,"generated #eta","");
00236   SETAXES(PhiGen,"generated #phi [rad]","");
00237   SETAXES(EtGen,"generated E_{T} [GeV]","");
00238   SETAXES(EtvsPhiGen,"generated #phi [rad]","generated E_{T} [GeV]");
00239   SETAXES(EtvsEtaGen,"generated #eta","generated E_{T} [GeV]");
00240 
00241   SETAXES(NGen,"Number of Gen Objects","");
00242   
00243   if (doMetPlots)
00244   {
00245     SETAXES(SumEt,"SumEt [GeV]","");
00246     SETAXES(TrueSumEt,"TrueSumEt [GeV]","");
00247     SETAXES(DeltaSetvsSet,"TrueSumEt","#DeltaSumEt [GeV]");
00248     SETAXES(DeltaMexvsSet,"TrueSumEt","#DeltaMEX [GeV]");
00249     SETAXES(DeltaSetOverSetvsSet,"TrueSumEt","#DeltaSumEt/trueSumEt");
00250     SETAXES(RecSetvsTrueSet,"TrueSumEt","SumEt");
00251     SETAXES(RecSetOverTrueSetvsTrueSet,"TrueSumEt","SumEt/trueSumEt");
00252     SETAXES(TrueMexvsTrueSet,"TrueSumEt","TrueMEX");
00253   }
00254 
00255   TDirectory* oldpwd = gDirectory;
00256 
00257 
00258   TIter next( gROOT->GetListOfFiles() );
00259 
00260   const bool debug=false;
00261 
00262   while ( TFile *file = (TFile *)next() )
00263   {
00264     if (debug) cout<<"file "<<file->GetName()<<endl;
00265   }
00266   if (DQM)
00267   {
00268     cout<<"DQM subdir"<<endl;
00269     cout<< DQM->pwd().c_str()<<endl;
00270 
00271     DQM->cd( DQM->pwd() );
00272   }
00273 
00274   if (debug)
00275   {
00276     cout<<"current dir"<<endl;
00277     gDirectory->pwd();
00278   }
00279   
00280   oldpwd->cd();
00281   //gDirectory->pwd();
00282 
00283   doMetPlots_=doMetPlots;
00284 
00285   //   tree_ = new BenchmarkTree("genericBenchmark", "Generic Benchmark TTree");
00286 }
00287 
00288 bool GenericBenchmark::accepted(const reco::Candidate* particle,
00289                                 double ptCut,
00290                                 double minEtaCut,
00291                                 double maxEtaCut ) const {
00292  
00293   //skip reconstructed PFJets with p_t < recPt_cut
00294   if (particle->pt() < ptCut and ptCut != -1.)
00295     return false;
00296 
00297   if (fabs(particle->eta())>maxEtaCut and maxEtaCut > 0)
00298     return false;
00299   if (fabs(particle->eta())<minEtaCut and minEtaCut > 0)
00300     return false;
00301 
00302   //accepted!
00303   return true;
00304  
00305 }
00306 
00307 
00308 void GenericBenchmark::fillHistos( const reco::Candidate* genParticle,
00309                                    const reco::Candidate* recParticle,
00310                                    double deltaR_cut,
00311                                    bool plotAgainstReco ) {
00312 
00313   // get the quantities to place on the denominator and/or divide by
00314   double et = genParticle->et();
00315   double eta = genParticle->eta();
00316   double phi = genParticle->phi();
00317   //std::cout << "FL : et = " << et << std::endl;
00318   //std::cout << "FL : eta = " << eta << std::endl;
00319   //std::cout << "FL : phi = " << phi << std::endl;
00320   //std::cout << "FL : rec et = " << recParticle->et() << std::endl;
00321   //std::cout << "FL : rec eta = " << recParticle->eta() << std::endl;
00322   //std::cout << "FL : rec phi = " <<recParticle-> phi() << std::endl;
00323 
00324   if (plotAgainstReco) { 
00325     et = recParticle->et(); 
00326     eta = recParticle->eta(); 
00327     phi = recParticle->phi(); 
00328   }
00329   
00330     
00331   // get the delta quantities
00332   double deltaEt = algo_->deltaEt(recParticle,genParticle);
00333   double deltaR = algo_->deltaR(recParticle,genParticle);
00334   double deltaEta = algo_->deltaEta(recParticle,genParticle);
00335   double deltaPhi = algo_->deltaPhi(recParticle,genParticle);
00336    
00337   //std::cout << "FL :deltaR_cut = " << deltaR_cut << std::endl;
00338   //std::cout << "FL :deltaR = " << deltaR << std::endl;
00339 
00340   if (deltaR>deltaR_cut && deltaR_cut != -1.)
00341     return;  
00342 
00343   hDeltaEt->Fill(deltaEt);
00344   hDeltaEx->Fill(recParticle->px()-genParticle->px());
00345   hDeltaEy->Fill(recParticle->py()-genParticle->py());
00346   hDeltaEtvsEt->Fill(et,deltaEt);
00347   hDeltaEtOverEtvsEt->Fill(et,deltaEt/et);
00348   hDeltaEtvsEta->Fill(eta,deltaEt);
00349   hDeltaEtOverEtvsEta->Fill(eta,deltaEt/et);
00350   hDeltaEtvsPhi->Fill(phi,deltaEt);
00351   hDeltaEtOverEtvsPhi->Fill(phi,deltaEt/et);
00352   hDeltaEtvsDeltaR->Fill(deltaR,deltaEt);
00353   hDeltaEtOverEtvsDeltaR->Fill(deltaR,deltaEt/et);
00354     
00355   hDeltaEta->Fill(deltaEta);
00356   hDeltaEtavsEt->Fill(et,deltaEta);
00357   hDeltaEtavsEta->Fill(eta,deltaEta);
00358     
00359   hDeltaPhi->Fill(deltaPhi);
00360   hDeltaPhivsEt->Fill(et,deltaPhi);
00361   hDeltaPhivsEta->Fill(eta,deltaPhi);
00362 
00363   hDeltaR->Fill(deltaR);
00364   hDeltaRvsEt->Fill(et,deltaR);
00365   hDeltaRvsEta->Fill(eta,deltaR);
00366 
00367   BenchmarkTreeEntry entry;
00368   entry.deltaEt = deltaEt;
00369   entry.deltaEta = deltaEta;
00370   entry.et = et;
00371   entry.eta = eta;
00372 
00373   if (doMetPlots_)
00374   {
00375     const reco::MET* met1=static_cast<const reco::MET*>(genParticle);
00376     const reco::MET* met2=static_cast<const reco::MET*>(recParticle);
00377     if (met1!=NULL && met2!=NULL)
00378     {
00379       //std::cout << "FL: met1.sumEt() = " << (*met1).sumEt() << std::endl;
00380       hTrueSumEt->Fill((*met1).sumEt());
00381       hSumEt->Fill((*met2).sumEt());
00382       hDeltaSetvsSet->Fill((*met1).sumEt(),(*met2).sumEt()-(*met1).sumEt());
00383       hDeltaMexvsSet->Fill((*met1).sumEt(),recParticle->px()-genParticle->px());
00384       hDeltaMexvsSet->Fill((*met1).sumEt(),recParticle->py()-genParticle->py());
00385       if ((*met1).sumEt()>0.01) hDeltaSetOverSetvsSet->Fill((*met1).sumEt(),((*met2).sumEt()-(*met1).sumEt())/(*met1).sumEt());
00386       hRecSetvsTrueSet->Fill((*met1).sumEt(),(*met2).sumEt());
00387       hRecSetOverTrueSetvsTrueSet->Fill((*met1).sumEt(),(*met2).sumEt()/((*met1).sumEt()));
00388       hTrueMexvsTrueSet->Fill((*met1).sumEt(),(*met1).px());
00389       hTrueMexvsTrueSet->Fill((*met1).sumEt(),(*met1).py());
00390     }
00391     else
00392     {
00393       std::cout << "Warning : reco::MET* == NULL" << std::endl;
00394     }
00395   }
00396 
00397   //     tree_->Fill(entry);
00398 
00399 }
00400 
00401 void GenericBenchmark::write(std::string Filename) {
00402 
00403   if (Filename.size() != 0 && file_)
00404     file_->Write(Filename.c_str());
00405 
00406 }
00407 
00408 void GenericBenchmark::setfile(TFile *file)
00409 {
00410   file_=file;
00411 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/RecoParticleFlow/Benchmark/src/GenericBenchmark.cc
