Reference Manual

Go to the documentation of this file.
00001 #include <TFile.h>
00002 #include <TH1.h>
00003 #include <TH2.h>
00004 #include <TStyle.h>
00005 #include <TPaveStats.h>
00006 
00007 
00008 #include <sstream>
00009 #include <cstdlib>
00010 #include <cassert>
00011 
00012 #include "Validation/RecoParticleFlow/interface/Comparator.h"
00013 #include "Validation/RecoParticleFlow/interface/TH2Analyzer.h"
00014 #include "Validation/RecoParticleFlow/interface/NicePlot.h"
00015 
00016 using namespace std; 
00017 
00018 void Comparator::SetDirs( const char* file0,
00019                           const char* dir0,
00020                           const char* file1,
00021                           const char* dir1  ) {
00022   
00023   file0_ = new TFile( file0 );
00024   if( file0_->IsZombie() ) exit(1);
00025   dir0_ = file0_->GetDirectory( dir0 );
00026   if(! dir0_ ) exit(1);
00027   
00028   file1_ = new TFile( file1 );
00029   if( file1_->IsZombie() ) exit(1);
00030   dir1_ = file1_->GetDirectory( dir1 );
00031   if(! dir1_ ) exit(1);
00032 }
00033 
00034 
00035 void Comparator::SetStyles( Style* s0, 
00036                   Style* s1,
00037                   const char* leg0,
00038                   const char* leg1) {
00039   s0_ = s0; 
00040   s1_ = s1;
00041   
00042   legend_.Clear();
00043   legend_.AddEntry( s0_, leg0, "mlf");
00044   legend_.AddEntry( s1_, leg1, "mlf");
00045 }
00046   
00047 void Comparator::DrawSlice( const char* key, 
00048                             int binxmin, int binxmax, 
00049                             Mode mode ) {
00050     
00051   static int num = 0;
00052     
00053   ostringstream out0;
00054   out0<<"h0_2d_"<<num;
00055   ostringstream out1;
00056   out1<<"h1_2d_"<<num;
00057   num++;
00058 
00059   string name0 = out0.str();
00060   string name1 = out1.str();
00061       
00062 
00063   TH1* h0 = Histo( key, 0);
00064   TH1* h1 = Histo( key, 1);
00065 
00066   TH2* h0_2d = dynamic_cast< TH2* >(h0);
00067   TH2* h1_2d = dynamic_cast< TH2* >(h1);
00068     
00069   if(h0_2d->GetNbinsY() == 1 || 
00070      h1_2d->GetNbinsY() == 1 ) {
00071     cerr<<key<<" is not 2D"<<endl;
00072     return;
00073   }
00074     
00075   TH1::AddDirectory( false );
00076 
00077   TH1D* h0_slice = h0_2d->ProjectionY(name0.c_str(),
00078                                       binxmin, binxmax, "");
00079   TH1D* h1_slice = h1_2d->ProjectionY(name1.c_str(),
00080                                       binxmin, binxmax, "");
00081   TH1::AddDirectory( true );
00082   Draw( h0_slice, h1_slice, mode);        
00083 }
00084 
00085 void Comparator::DrawMeanSlice(const char* key, const int rebinFactor, Mode mode)
00086 {
00087   TDirectory* dir = dir1_;
00088   dir->cd();
00089   TH2D *h2 = (TH2D*) dir->Get(key);
00090   TH2Analyzer TH2Ana(h2,rebinFactor);
00091   TH1D* ha=TH2Ana.Average();
00092 
00093   dir = dir0_;
00094   dir->cd();
00095   TH2D *h2b = (TH2D*) dir->Get(key);
00096   TH2Analyzer TH2Anab(h2b,rebinFactor);
00097   TH1D* hb=TH2Anab.Average();
00098 
00099   Draw(hb,ha,mode);
00100 }
00101 
00102 void Comparator::DrawSigmaSlice(const char* key, const int rebinFactor, Mode mode)
00103 {
00104   TDirectory* dir = dir1_;
00105   dir->cd();
00106   TH2D *h2 = (TH2D*) dir->Get(key);
00107   TH2Analyzer TH2Ana(h2,rebinFactor);
00108   TH1D* ha=TH2Ana.RMS();
00109 
00110   dir = dir0_;
00111   dir->cd();
00112   TH2D *h2b = (TH2D*) dir->Get(key);
00113   TH2Analyzer TH2Anab(h2b,rebinFactor);
00114   TH1D* hb=TH2Anab.RMS();
00115 
00116   Draw(hb,ha,mode);
00117 }
00118 
00119 void Comparator::DrawGaussSigmaSlice(const char* key, const int rebinFactor, Mode mode)
00120 {
00121   TDirectory* dir = dir1_;
00122   dir->cd();
00123   TH2D *h2 = (TH2D*) dir->Get(key);
00124   TH2Analyzer TH2Ana(h2,rebinFactor);
00125   TH1D* ha=TH2Ana.SigmaGauss();
00126 
00127   dir = dir0_;
00128   dir->cd();
00129   TH2D *h2b = (TH2D*) dir->Get(key);
00130   TH2Analyzer TH2Anab(h2b,rebinFactor);
00131   TH1D* hb=TH2Anab.SigmaGauss();
00132 
00133   Draw(hb,ha,mode);
00134 }
00135 
00136 Double_t fitFunction_f(Double_t *x, Double_t *par)
00137 {
00138   const Double_t value=sqrt(par[0]*par[0]+par[1]*par[1]*(x[0]-par[3])+par[2]*par[2]*(x[0]-par[3])*(x[0]-par[3]))/x[0];
00139   return value;
00140 }
00141 
00142 void Comparator::DrawGaussSigmaSlice(const char* key, const int rebinFactor, const int binxmin,
00143                                      const int binxmax, const bool cst_binning, Mode mode)
00144 {
00145   TDirectory* dir = dir1_;
00146   dir->cd();
00147   TH2D *h2 = (TH2D*) dir->Get(key);
00148   TH2Analyzer TH2Ana(h2, binxmin, binxmax, rebinFactor, cst_binning);
00149   TH1D* hrms=TH2Ana.RMS();
00150   TF1 *fitfcndgssrms3 = new TF1("fitfcndgssrms3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
00151   fitfcndgssrms3->SetNpx(500);
00152   fitfcndgssrms3->SetLineWidth(3);
00153   fitfcndgssrms3->SetLineStyle(2);
00154   fitfcndgssrms3->SetLineColor(4);
00155   hrms->Fit("fitfcndgssrms3","0R");
00156 
00157   TH1D* ha=TH2Ana.SigmaGauss();
00158   TF1 *fitfcndgsse3 = new TF1("fitfcndgsse3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
00159   fitfcndgsse3->SetNpx(500);
00160   fitfcndgsse3->SetLineWidth(3);
00161   fitfcndgsse3->SetLineStyle(1);
00162   fitfcndgsse3->SetLineColor(4);
00163   ha->Fit("fitfcndgsse3","0R");
00164 
00165   dir = dir0_;
00166   dir->cd();
00167   TH2D *h2b = (TH2D*) dir->Get(key);
00168   TH2Analyzer TH2Anab(h2b, binxmin, binxmax, rebinFactor, cst_binning);
00169   TH1D* hrmsb=TH2Anab.RMS();
00170   TF1 *fitfcndgssrmsb3 = new TF1("fitfcndgssrmsb3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
00171   fitfcndgssrmsb3->SetNpx(500);
00172   fitfcndgssrmsb3->SetLineWidth(3);
00173   fitfcndgssrmsb3->SetLineStyle(2);
00174   fitfcndgssrmsb3->SetLineColor(2);
00175   hrmsb->Fit("fitfcndgssrmsb3","0R");
00176 
00177   TH1D* hb=TH2Anab.SigmaGauss();
00178   TF1 *fitfcndgsseb3 = new TF1("fitfcndgsseb3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
00179   fitfcndgsseb3->SetNpx(500);
00180   fitfcndgsseb3->SetLineWidth(3);
00181   fitfcndgsseb3->SetLineStyle(1);
00182   fitfcndgsseb3->SetLineColor(2);
00183   hb->Fit("fitfcndgsseb3","0R");
00184 
00185   Draw(hb,ha,mode);
00186   //Draw(hrms,ha,mode);
00187   //Draw(ha,ha,mode);
00188   fitfcndgssrms3->Draw("same");
00189   fitfcndgsse3->Draw("same");
00190   fitfcndgssrmsb3->Draw("same");
00191   fitfcndgsseb3->Draw("same");
00192 }
00193 
00194 void Comparator::DrawGaussSigmaOverMeanXSlice(const char* key, const int rebinFactor, const int binxmin,
00195                                      const int binxmax, const bool cst_binning, Mode mode)
00196 {
00197 
00198   TDirectory* dir = dir1_;
00199   dir->cd();
00200   TH2D *h2 = (TH2D*) dir->Get(key);
00201   TH2Analyzer TH2Ana(h2, binxmin, binxmax, rebinFactor, cst_binning);
00202   TH1D* hrms=TH2Ana.RMS();
00203 
00204   TH1D* meanXslice = TH2Ana.MeanX();
00205   //for( int i=1; i<=meanXslice->GetNbinsX(); ++i) {
00206   //  std::cout << "meanXslice->GetBinContent(" << i << ") = "
00207   //          << meanXslice->GetBinContent(i) << std::endl;
00208   //  std::cout << "meanXslice->GetBinError(" << i << ") = "
00209   //          << meanXslice->GetBinError(i) << std::endl;
00210   //}
00211   //Draw(meanXslice,meanXslice,mode);
00212   hrms->Divide(meanXslice);
00213 
00214   TF1 *fitXfcndgssrms3 = new TF1("fitXfcndgssrms3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
00215   fitXfcndgssrms3->SetNpx(500);
00216   fitXfcndgssrms3->SetLineWidth(3);
00217   fitXfcndgssrms3->SetLineStyle(2);
00218   fitXfcndgssrms3->SetLineColor(4);
00219   hrms->Fit("fitXfcndgssrms3","0R");
00220 
00221   TH1D* ha=TH2Ana.SigmaGauss();
00222   ha->Divide(meanXslice);
00223   TF1 *fitXfcndgsse3 = new TF1("fitXfcndgsse3",fitFunction_f,ha->GetXaxis()->GetBinLowEdge(1),ha->GetXaxis()->GetBinUpEdge(ha->GetNbinsX()),4);
00224   fitXfcndgsse3->SetNpx(500);
00225   fitXfcndgsse3->SetLineWidth(3);
00226   fitXfcndgsse3->SetLineStyle(1);
00227   fitXfcndgsse3->SetLineColor(4);
00228   ha->Fit("fitXfcndgsse3","0R");
00229 
00230   dir = dir0_;
00231   dir->cd();
00232   TH2D *h2b = (TH2D*) dir->Get(key);
00233   TH2Analyzer TH2Anab(h2b, binxmin, binxmax, rebinFactor, cst_binning);
00234   TH1D* hrmsb=TH2Anab.RMS();
00235   hrmsb->Divide(meanXslice);
00236   TF1 *fitXfcndgssrmsb3 = new TF1("fitXfcndgssrmsb3",fitFunction_f,hrmsb->GetXaxis()->GetBinLowEdge(1),hrmsb->GetXaxis()->GetBinUpEdge(hrmsb->GetNbinsX()),4);
00237   fitXfcndgssrmsb3->SetNpx(500);
00238   fitXfcndgssrmsb3->SetLineWidth(3);
00239   fitXfcndgssrmsb3->SetLineStyle(2);
00240   fitXfcndgssrmsb3->SetLineColor(2);
00241   hrmsb->Fit("fitXfcndgssrmsb3","0R");
00242 
00243   TH1D* hb=TH2Anab.SigmaGauss();
00244   hb->Divide(meanXslice);
00245   TF1 *fitXfcndgsseb3 = new TF1("fitXfcndgsseb3",fitFunction_f,hb->GetXaxis()->GetBinLowEdge(1),hb->GetXaxis()->GetBinUpEdge(hb->GetNbinsX()),4);
00246   fitXfcndgsseb3->SetNpx(500);
00247   fitXfcndgsseb3->SetLineWidth(3);
00248   fitXfcndgsseb3->SetLineStyle(1);
00249   fitXfcndgsseb3->SetLineColor(2);
00250   hb->Fit("fitXfcndgsseb3","0R");
00251 
00252   Draw(hb,ha,mode);
00253   //Draw(hrms,ha,mode);
00254   //Draw(ha,ha,mode);
00255   fitXfcndgssrms3->Draw("same");
00256   fitXfcndgsse3->Draw("same");
00257   fitXfcndgssrmsb3->Draw("same");
00258   fitXfcndgsseb3->Draw("same");
00259 }
00260 
00261 void Comparator::DrawGaussSigmaOverMeanSlice(const char* key, const char* key2,
00262                                              const int rebinFactor, Mode mode)
00263 {
00264   TDirectory* dir = dir1_;
00265   dir->cd();
00266 
00267   TH2D *h2_b = (TH2D*) dir->Get(key2);
00268   TH2Analyzer TH2Ana_b(h2_b, rebinFactor);
00269   TH1D* meanslice=TH2Ana_b.Average();
00270 
00271   TH2D *h2 = (TH2D*) dir->Get(key);
00272   TH2Analyzer TH2Ana(h2, rebinFactor);
00273 
00274   TH1D* ha=TH2Ana.SigmaGauss();
00275   ha->Divide(meanslice);
00276 
00277   dir = dir0_;
00278   dir->cd();
00279   TH2D *h2b = (TH2D*) dir->Get(key);
00280   TH2Analyzer TH2Anab(h2b, rebinFactor);
00281 
00282   TH2D *h2b_b = (TH2D*) dir->Get(key2);
00283   TH2Analyzer TH2Anab_b(h2b_b, rebinFactor);
00284   TH1D* meansliceb=TH2Anab_b.Average();
00285 
00286   TH1D* hb=TH2Anab.SigmaGauss();
00287   hb->Divide(meansliceb);
00288 
00289   Draw(hb,ha,mode);
00290   //Draw(meansliceb,meanslice,mode);
00291 }
00292 
00293 void Comparator::Draw( const char* key, Mode mode) {
00294 
00295   TH1::AddDirectory( false );
00296   TH1* h0 = Histo( key, 0);
00297   TH1* h1 = (TH1*) Histo( key, 1)->Clone("h1");
00298 
00299   TH1::AddDirectory( true );
00300   Draw( h0, h1, mode);    
00301 }
00302 
00303   
00304 void Comparator::Draw( const char* key0, const char* key1, Mode mode) {
00305   TH1* h0=0;
00306   TH1* h1=0;
00307   if(mode!=EFF)
00308     {
00309       h0 = Histo( key0, 0);
00310       h1 = Histo( key1, 1);
00311     } 
00312   else
00313     {
00314       h0 = Histo( key0, 0);
00315       TH1 * h0b = Histo( key1, 0);
00316       h1 = Histo( key0, 1);
00317       TH1 * h1b = Histo( key1, 1);
00318       if(rebin_>1) {
00319         h0->Rebin( rebin_);
00320         h1->Rebin( rebin_);
00321         h0b->Rebin( rebin_);
00322         h1b->Rebin( rebin_);
00323       }
00324       if(resetAxis_) {
00325         h0->GetXaxis()->SetRangeUser( xMin_, xMax_);
00326         h1->GetXaxis()->SetRangeUser( xMin_, xMax_);
00327         h0b->GetXaxis()->SetRangeUser( xMin_, xMax_);
00328         h1b->GetXaxis()->SetRangeUser( xMin_, xMax_);
00329       }      
00330 
00331       h0b->Sumw2();
00332       h0->Sumw2();
00333       h0->Divide(h0,h0b,1.,1.,"B");
00334       h1b->Sumw2();
00335       h1->Sumw2();
00336       h1->Divide(h1,h1b,1.,1.,"B");
00337     }
00338   Draw( h0, h1, mode);
00339 }
00340 
00341 TH1* Comparator::Histo( const char* key, unsigned dirIndex) {
00342   if(dirIndex>1U) { // dirIndex >= 0, since dirIndex is unsigned
00343     cerr<<"bad dir index: "<<dirIndex<<endl;
00344     return 0;
00345   }
00346   TDirectory* dir = 0;
00347   if(dirIndex == 0) dir = dir0_;
00348   if(dirIndex == 1) dir = dir1_;
00349   assert( dir );
00350   
00351   dir->cd();
00352 
00353   TH1* h = (TH1*) dir->Get(key);
00354   if(!h)  
00355     cerr<<"no key "<<key<<" in directory "<<dir->GetName()<<endl;
00356   return h;
00357 }
00358 
00359 
00360 void Comparator::Draw( TH1* h0, TH1* h1, Mode mode ) {
00361   if( !(h0 && h1) ) { 
00362     cerr<<"invalid histo"<<endl;
00363     return;
00364   }
00365     
00366   TH1::AddDirectory( false );
00367   h0_ = (TH1*) h0->Clone( "h0_");
00368   h1_ = (TH1*) h1->Clone( "h1_");
00369   TH1::AddDirectory( true );
00370     
00371   // unsetting the title, since the title of projections
00372   // is still the title of the 2d histo
00373   // and this is better anyway
00374   h0_->SetTitle("");
00375   h1_->SetTitle("");    
00376 
00377   //h0_->SetStats(1);
00378   //h1_->SetStats(1);
00379 
00380   if(mode!=EFF)
00381     {
00382       if(rebin_>1) {
00383         h0_->Rebin( rebin_);
00384         h1_->Rebin( rebin_);
00385       }
00386       if(resetAxis_) {
00387         h0_->GetXaxis()->SetRangeUser( xMin_, xMax_);
00388         h1_->GetXaxis()->SetRangeUser( xMin_, xMax_);
00389       }
00390     }
00391 
00392   if (mode!=GRAPH)
00393   {
00394 
00395     TPaveStats *ptstats = new TPaveStats(0.7385057,0.720339,
00396                                        0.9396552,0.8792373,"brNDC");
00397     ptstats->SetName("stats");
00398     ptstats->SetBorderSize(1);
00399     ptstats->SetLineColor(2);
00400     ptstats->SetFillColor(10);
00401     ptstats->SetTextAlign(12);
00402     ptstats->SetTextColor(2);
00403     ptstats->SetOptStat(1111);
00404     ptstats->SetOptFit(0);
00405     ptstats->Draw();
00406     h0_->GetListOfFunctions()->Add(ptstats);
00407     ptstats->SetParent(h0_->GetListOfFunctions());
00408   
00409     //std::cout << "FL: h0_->GetMean() = " << h0_->GetMean() << std::endl;
00410     //std::cout << "FL: h0_->GetRMS() = " << h0_->GetRMS() << std::endl;
00411     //std::cout << "FL: h1_->GetMean() = " << h1_->GetMean() << std::endl;
00412     //std::cout << "FL: h1_->GetRMS() = " << h1_->GetRMS() << std::endl;
00413     //std::cout << "FL: test2" << std::endl;
00414     TPaveStats *ptstats2 = new TPaveStats(0.7399425,0.529661,
00415                                         0.941092,0.6885593,"brNDC");
00416     ptstats2->SetName("stats");
00417     ptstats2->SetBorderSize(1);
00418     ptstats2->SetLineColor(4);
00419     ptstats2->SetFillColor(10);
00420     ptstats2->SetTextAlign(12);
00421     ptstats2->SetTextColor(4);
00422     TText *text = ptstats2->AddText("h1_");
00423     text->SetTextSize(0.03654661);
00424   
00425     std::ostringstream oss3;
00426     oss3 << h1_->GetEntries();
00427     const std::string txt_entries="Entries = "+oss3.str();
00428     text = ptstats2->AddText(txt_entries.c_str());
00429     std::ostringstream oss;
00430     oss << h1_->GetMean();
00431     const std::string txt_mean="Mean  = "+oss.str();
00432     text = ptstats2->AddText(txt_mean.c_str());
00433     std::ostringstream oss2;
00434     oss2 << h1_->GetRMS();
00435     const std::string txt_rms="RMS  = "+oss2.str();
00436     text = ptstats2->AddText(txt_rms.c_str());
00437     ptstats2->SetOptStat(1111);
00438     ptstats2->SetOptFit(0);
00439     ptstats2->Draw();
00440     h1_->GetListOfFunctions()->Add(ptstats2);
00441     ptstats2->SetParent(h1_->GetListOfFunctions());
00442   }
00443   else
00444   {
00445     TPaveStats *ptstats = new TPaveStats(0.0,0.0,
00446                                          0.0,0.0,"brNDC");
00447     ptstats->Draw();
00448     h0_->GetListOfFunctions()->Add(ptstats);
00449     ptstats->SetParent(h0_->GetListOfFunctions());
00450   }
00451 
00452   float min=-999.;
00453   float max=+999.;
00454   switch(mode) {
00455   case SCALE:
00456     h1_->Scale( h0_->GetEntries()/h1_->GetEntries() );
00457   case NORMAL:
00458     if(s0_)
00459       Styles::FormatHisto( h0_ , s0_);
00460     if(s1_)
00461       Styles::FormatHisto( h1_ , s1_);
00462       
00463     if( h1_->GetMaximum()>h0_->GetMaximum()) {
00464       h0_->SetMaximum( h1_->GetMaximum()*1.15 );
00465     }
00466 
00467     h0_->Draw();
00468     h1_->Draw("same");
00469 
00470     break;
00471   case EFF:
00472     if(s0_)
00473       Styles::FormatHisto( h0_ , s0_);
00474     if(s1_)
00475       Styles::FormatHisto( h1_ , s1_);
00476 
00477     // rather arbitrary but useful
00478     max= h0_->GetMaximum();
00479     if ( h1_->GetMaximum() > max )
00480       max = h1_->GetMaximum();
00481     if(max > 0.8) max = 1;
00482     
00483     max *= 1.1 ;
00484 
00485     min= h0_->GetMinimum();
00486     if ( h1_->GetMinimum() < min )
00487       min = h1_->GetMinimum();
00488     if(min > 0.2) min = 0.;
00489 
00490     min *= 0.8 ;
00491     h0_->SetMaximum( max );
00492     h0_->SetMinimum( min );
00493 
00494     h0_->Draw("E");
00495     h1_->Draw("Esame");
00496     break;
00497   case GRAPH:
00498     if(s0_)
00499       Styles::FormatHisto( h0_ , s0_);
00500     if(s1_)
00501       Styles::FormatHisto( h1_ , s1_);
00502       
00503     if( h1_->GetMaximum()>h0_->GetMaximum()) {
00504       h0_->SetMaximum( h1_->GetMaximum()*1.15 );
00505     }
00506     if( h1_->GetMinimum()<h0_->GetMinimum()) {
00507       h0_->SetMinimum( h1_->GetMinimum()*1.15 );
00508     }
00509     h0_->SetMarkerStyle(21);
00510     h0_->SetMarkerColor(2);
00511     h1_->SetMarkerStyle(21);
00512     h1_->SetMarkerColor(4);
00513 
00514     h0_->Draw("E1");
00515     h1_->Draw("E1same");
00516     break;
00517   case RATIO:
00518     h0_->Sumw2();
00519     h1_->Sumw2();
00520     h0_->Divide( h1_ );
00521     if(s0_)
00522       Styles::FormatHisto( h0_ , s0_);
00523     h0_->Draw();
00524   default:
00525     break;
00526   }
00527 }
00528
CMSSW_4_4_3_patch1/src/Validation/RecoParticleFlow/src/Comparator.cc
