Comparator.cc

Go to the documentation of this file.

#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <TStyle.h>
#include <TPaveStats.h>


#include <sstream>
#include <cstdlib>
#include <cassert>

#include "Validation/RecoParticleFlow/interface/Comparator.h"
#include "Validation/RecoParticleFlow/interface/TH2Analyzer.h"
#include "Validation/RecoParticleFlow/interface/NicePlot.h"

using namespace std; 

void Comparator::SetDirs( const char* file0,
              const char* dir0,
              const char* file1,
              const char* dir1  ) {
  
  file0_ = new TFile( file0 );
  if( file0_->IsZombie() ) exit(1);
  dir0_ = file0_->GetDirectory( dir0 );
  if(! dir0_ ) exit(1);
  
  file1_ = new TFile( file1 );
  if( file1_->IsZombie() ) exit(1);
  dir1_ = file1_->GetDirectory( dir1 );
  if(! dir1_ ) exit(1);
}


void Comparator::SetStyles( Style* s0, 
          Style* s1,
          const char* leg0,
          const char* leg1) {
  s0_ = s0; 
  s1_ = s1;
  
  legend_.Clear();
  legend_.AddEntry( s0_, leg0, "mlf");
  legend_.AddEntry( s1_, leg1, "mlf");
}
  
void Comparator::DrawSlice( const char* key, 
                int binxmin, int binxmax, 
                Mode mode ) {
    
  static int num = 0;
    
  ostringstream out0;
  out0<<"h0_2d_"<<num;
  ostringstream out1;
  out1<<"h1_2d_"<<num;
  num++;

  string name0 = out0.str();
  string name1 = out1.str();
      

  TH1* h0 = Histo( key, 0);
  TH1* h1 = Histo( key, 1);

  TH2* h0_2d = dynamic_cast< TH2* >(h0);
  TH2* h1_2d = dynamic_cast< TH2* >(h1);
    
  if(h0_2d->GetNbinsY() == 1 || 
     h1_2d->GetNbinsY() == 1 ) {
    cerr<<key<<" is not 2D"<<endl;
    return;
  }
    
  TH1::AddDirectory( false );

  TH1D* h0_slice = h0_2d->ProjectionY(name0.c_str(),
                      binxmin, binxmax, "");
  TH1D* h1_slice = h1_2d->ProjectionY(name1.c_str(),
                      binxmin, binxmax, "");
  TH1::AddDirectory( true );
  Draw( h0_slice, h1_slice, mode);        
}

void Comparator::DrawMeanSlice(const char* key, const int rebinFactor, Mode mode)
{
  TDirectory* dir = dir1_;
  dir->cd();
  TH2D *h2 = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Ana(h2,rebinFactor);
  TH1D* ha=TH2Ana.Average();

  dir = dir0_;
  dir->cd();
  TH2D *h2b = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Anab(h2b,rebinFactor);
  TH1D* hb=TH2Anab.Average();

  Draw(hb,ha,mode);
}

void Comparator::DrawSigmaSlice(const char* key, const int rebinFactor, Mode mode)
{
  TDirectory* dir = dir1_;
  dir->cd();
  TH2D *h2 = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Ana(h2,rebinFactor);
  TH1D* ha=TH2Ana.RMS();

  dir = dir0_;
  dir->cd();
  TH2D *h2b = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Anab(h2b,rebinFactor);
  TH1D* hb=TH2Anab.RMS();

  Draw(hb,ha,mode);
}

void Comparator::DrawGaussSigmaSlice(const char* key, const int rebinFactor, Mode mode)
{
  TDirectory* dir = dir1_;
  dir->cd();
  TH2D *h2 = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Ana(h2,rebinFactor);
  TH1D* ha=TH2Ana.SigmaGauss();

  dir = dir0_;
  dir->cd();
  TH2D *h2b = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Anab(h2b,rebinFactor);
  TH1D* hb=TH2Anab.SigmaGauss();

  Draw(hb,ha,mode);
}

Double_t fitFunction_f(Double_t *x, Double_t *par)
{
  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];
  return value;
}

void Comparator::DrawGaussSigmaSlice(const char* key, const int rebinFactor, const int binxmin,
                     const int binxmax, const bool cst_binning, Mode mode)
{
  TDirectory* dir = dir1_;
  dir->cd();
  TH2D *h2 = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Ana(h2, binxmin, binxmax, rebinFactor, cst_binning);
  TH1D* hrms=TH2Ana.RMS();
  TF1 *fitfcndgssrms3 = new TF1("fitfcndgssrms3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
  fitfcndgssrms3->SetNpx(500);
  fitfcndgssrms3->SetLineWidth(3);
  fitfcndgssrms3->SetLineStyle(2);
  fitfcndgssrms3->SetLineColor(4);
  hrms->Fit(fitfcndgssrms3,"0R");

  TH1D* ha=TH2Ana.SigmaGauss();
  TF1 *fitfcndgsse3 = new TF1("fitfcndgsse3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
  fitfcndgsse3->SetNpx(500);
  fitfcndgsse3->SetLineWidth(3);
  fitfcndgsse3->SetLineStyle(1);
  fitfcndgsse3->SetLineColor(4);
  ha->Fit(fitfcndgsse3,"0R");

  dir = dir0_;
  dir->cd();
  TH2D *h2b = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Anab(h2b, binxmin, binxmax, rebinFactor, cst_binning);
  TH1D* hrmsb=TH2Anab.RMS();
  TF1 *fitfcndgssrmsb3 = new TF1("fitfcndgssrmsb3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
  fitfcndgssrmsb3->SetNpx(500);
  fitfcndgssrmsb3->SetLineWidth(3);
  fitfcndgssrmsb3->SetLineStyle(2);
  fitfcndgssrmsb3->SetLineColor(2);
  hrmsb->Fit(fitfcndgssrmsb3,"0R");

  TH1D* hb=TH2Anab.SigmaGauss();
  TF1 *fitfcndgsseb3 = new TF1("fitfcndgsseb3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
  fitfcndgsseb3->SetNpx(500);
  fitfcndgsseb3->SetLineWidth(3);
  fitfcndgsseb3->SetLineStyle(1);
  fitfcndgsseb3->SetLineColor(2);
  hb->Fit(fitfcndgsseb3,"0R");

  Draw(hb,ha,mode);
  //Draw(hrms,ha,mode);
  //Draw(ha,ha,mode);
  fitfcndgssrms3->Draw("same");
  fitfcndgsse3->Draw("same");
  fitfcndgssrmsb3->Draw("same");
  fitfcndgsseb3->Draw("same");
}

void Comparator::DrawGaussSigmaOverMeanXSlice(const char* key, const int rebinFactor, const int binxmin,
                     const int binxmax, const bool cst_binning, Mode mode)
{

  TDirectory* dir = dir1_;
  dir->cd();
  TH2D *h2 = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Ana(h2, binxmin, binxmax, rebinFactor, cst_binning);
  TH1D* hrms=TH2Ana.RMS();

  TH1D* meanXslice = TH2Ana.MeanX();
  //for( int i=1; i<=meanXslice->GetNbinsX(); ++i) {
  //  std::cout << "meanXslice->GetBinContent(" << i << ") = "
  //          << meanXslice->GetBinContent(i) << std::endl;
  //  std::cout << "meanXslice->GetBinError(" << i << ") = "
  //          << meanXslice->GetBinError(i) << std::endl;
  //}
  //Draw(meanXslice,meanXslice,mode);
  hrms->Divide(meanXslice);

  TF1 *fitXfcndgssrms3 = new TF1("fitXfcndgssrms3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
  fitXfcndgssrms3->SetNpx(500);
  fitXfcndgssrms3->SetLineWidth(3);
  fitXfcndgssrms3->SetLineStyle(2);
  fitXfcndgssrms3->SetLineColor(4);
  hrms->Fit(fitXfcndgssrms3,"0R");

  TH1D* ha=TH2Ana.SigmaGauss();
  ha->Divide(meanXslice);
  TF1 *fitXfcndgsse3 = new TF1("fitXfcndgsse3",fitFunction_f,ha->GetXaxis()->GetBinLowEdge(1),ha->GetXaxis()->GetBinUpEdge(ha->GetNbinsX()),4);
  fitXfcndgsse3->SetNpx(500);
  fitXfcndgsse3->SetLineWidth(3);
  fitXfcndgsse3->SetLineStyle(1);
  fitXfcndgsse3->SetLineColor(4);
  ha->Fit(fitXfcndgsse3,"0R");

  dir = dir0_;
  dir->cd();
  TH2D *h2b = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Anab(h2b, binxmin, binxmax, rebinFactor, cst_binning);
  TH1D* hrmsb=TH2Anab.RMS();
  hrmsb->Divide(meanXslice);
  TF1 *fitXfcndgssrmsb3 = new TF1("fitXfcndgssrmsb3",fitFunction_f,hrmsb->GetXaxis()->GetBinLowEdge(1),hrmsb->GetXaxis()->GetBinUpEdge(hrmsb->GetNbinsX()),4);
  fitXfcndgssrmsb3->SetNpx(500);
  fitXfcndgssrmsb3->SetLineWidth(3);
  fitXfcndgssrmsb3->SetLineStyle(2);
  fitXfcndgssrmsb3->SetLineColor(2);
  hrmsb->Fit(fitXfcndgssrmsb3,"0R");

  TH1D* hb=TH2Anab.SigmaGauss();
  hb->Divide(meanXslice);
  TF1 *fitXfcndgsseb3 = new TF1("fitXfcndgsseb3",fitFunction_f,hb->GetXaxis()->GetBinLowEdge(1),hb->GetXaxis()->GetBinUpEdge(hb->GetNbinsX()),4);
  fitXfcndgsseb3->SetNpx(500);
  fitXfcndgsseb3->SetLineWidth(3);
  fitXfcndgsseb3->SetLineStyle(1);
  fitXfcndgsseb3->SetLineColor(2);
  hb->Fit(fitXfcndgsseb3,"0R");

  Draw(hb,ha,mode);
  //Draw(hrms,ha,mode);
  //Draw(ha,ha,mode);
  fitXfcndgssrms3->Draw("same");
  fitXfcndgsse3->Draw("same");
  fitXfcndgssrmsb3->Draw("same");
  fitXfcndgsseb3->Draw("same");
}

void Comparator::DrawGaussSigmaOverMeanSlice(const char* key, const char* key2,
                         const int rebinFactor, Mode mode)
{
  TDirectory* dir = dir1_;
  dir->cd();

  TH2D *h2_b = (TH2D*) dir->Get(key2);
  TH2Analyzer TH2Ana_b(h2_b, rebinFactor);
  TH1D* meanslice=TH2Ana_b.Average();

  TH2D *h2 = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Ana(h2, rebinFactor);

  TH1D* ha=TH2Ana.SigmaGauss();
  ha->Divide(meanslice);

  dir = dir0_;
  dir->cd();
  TH2D *h2b = (TH2D*) dir->Get(key);
  TH2Analyzer TH2Anab(h2b, rebinFactor);

  TH2D *h2b_b = (TH2D*) dir->Get(key2);
  TH2Analyzer TH2Anab_b(h2b_b, rebinFactor);
  TH1D* meansliceb=TH2Anab_b.Average();

  TH1D* hb=TH2Anab.SigmaGauss();
  hb->Divide(meansliceb);

  Draw(hb,ha,mode);
  //Draw(meansliceb,meanslice,mode);
}

void Comparator::Draw( const char* key, Mode mode) {

  TH1::AddDirectory( false );
  TH1* h0 = Histo( key, 0);
  TH1* h1 = (TH1*) Histo( key, 1)->Clone("h1");

  TH1::AddDirectory( true );
  Draw( h0, h1, mode);    
}

  
void Comparator::Draw( const char* key0, const char* key1, Mode mode) {
  TH1* h0=0;
  TH1* h1=0;
  if(mode!=EFF)
    {
      h0 = Histo( key0, 0);
      h1 = Histo( key1, 1);
    } 
  else
    {
      h0 = Histo( key0, 0);
      TH1 * h0b = Histo( key1, 0);
      h1 = Histo( key0, 1);
      TH1 * h1b = Histo( key1, 1);
      if(rebin_>1) {
    h0->Rebin( rebin_);
    h1->Rebin( rebin_);
    h0b->Rebin( rebin_);
    h1b->Rebin( rebin_);
      }
      if(resetAxis_) {
    h0->GetXaxis()->SetRangeUser( xMin_, xMax_);
    h1->GetXaxis()->SetRangeUser( xMin_, xMax_);
    h0b->GetXaxis()->SetRangeUser( xMin_, xMax_);
    h1b->GetXaxis()->SetRangeUser( xMin_, xMax_);
      }      

      h0b->Sumw2();
      h0->Sumw2();
      h0->Divide(h0,h0b,1.,1.,"B");
      h1b->Sumw2();
      h1->Sumw2();
      h1->Divide(h1,h1b,1.,1.,"B");
    }
  Draw( h0, h1, mode);
}

TH1* Comparator::Histo( const char* key, unsigned dirIndex) {
  if(dirIndex>1U) { // dirIndex >= 0, since dirIndex is unsigned
    cerr<<"bad dir index: "<<dirIndex<<endl;
    return 0;
  }
  TDirectory* dir = 0;
  if(dirIndex == 0) dir = dir0_;
  if(dirIndex == 1) dir = dir1_;
  assert( dir );
  
  dir->cd();

  TH1* h = (TH1*) dir->Get(key);
  if(!h)  
    cerr<<"no key "<<key<<" in directory "<<dir->GetName()<<endl;
  return h;
}


void Comparator::Draw( TH1* h0, TH1* h1, Mode mode ) {
  if( !(h0 && h1) ) { 
    cerr<<"invalid histo"<<endl;
    return;
  }
    
  TH1::AddDirectory( false );
  h0_ = (TH1*) h0->Clone( "h0_");
  h1_ = (TH1*) h1->Clone( "h1_");
  TH1::AddDirectory( true );
    
  // unsetting the title, since the title of projections
  // is still the title of the 2d histo
  // and this is better anyway
  h0_->SetTitle("");
  h1_->SetTitle("");    

  //h0_->SetStats(1);
  //h1_->SetStats(1);

  if(mode!=EFF)
    {
      if(rebin_>1) {
    h0_->Rebin( rebin_);
    h1_->Rebin( rebin_);
      }
      if(resetAxis_) {
    h0_->GetXaxis()->SetRangeUser( xMin_, xMax_);
    h1_->GetXaxis()->SetRangeUser( xMin_, xMax_);
      }
    }

  if (mode!=GRAPH)
  {

    TPaveStats *ptstats = new TPaveStats(0.7385057,0.720339,
                       0.9396552,0.8792373,"brNDC");
    ptstats->SetName("stats");
    ptstats->SetBorderSize(1);
    ptstats->SetLineColor(2);
    ptstats->SetFillColor(10);
    ptstats->SetTextAlign(12);
    ptstats->SetTextColor(2);
    ptstats->SetOptStat(1111);
    ptstats->SetOptFit(0);
    ptstats->Draw();
    h0_->GetListOfFunctions()->Add(ptstats);
    ptstats->SetParent(h0_->GetListOfFunctions());
  
    //std::cout << "FL: h0_->GetMean() = " << h0_->GetMean() << std::endl;
    //std::cout << "FL: h0_->GetRMS() = " << h0_->GetRMS() << std::endl;
    //std::cout << "FL: h1_->GetMean() = " << h1_->GetMean() << std::endl;
    //std::cout << "FL: h1_->GetRMS() = " << h1_->GetRMS() << std::endl;
    //std::cout << "FL: test2" << std::endl;
    TPaveStats *ptstats2 = new TPaveStats(0.7399425,0.529661,
                    0.941092,0.6885593,"brNDC");
    ptstats2->SetName("stats");
    ptstats2->SetBorderSize(1);
    ptstats2->SetLineColor(4);
    ptstats2->SetFillColor(10);
    ptstats2->SetTextAlign(12);
    ptstats2->SetTextColor(4);
    TText *text = ptstats2->AddText("h1_");
    text->SetTextSize(0.03654661);
  
    std::ostringstream oss3;
    oss3 << h1_->GetEntries();
    const std::string txt_entries="Entries = "+oss3.str();
    text = ptstats2->AddText(txt_entries.c_str());
    std::ostringstream oss;
    oss << h1_->GetMean();
    const std::string txt_mean="Mean  = "+oss.str();
    text = ptstats2->AddText(txt_mean.c_str());
    std::ostringstream oss2;
    oss2 << h1_->GetRMS();
    const std::string txt_rms="RMS  = "+oss2.str();
    text = ptstats2->AddText(txt_rms.c_str());
    ptstats2->SetOptStat(1111);
    ptstats2->SetOptFit(0);
    ptstats2->Draw();
    h1_->GetListOfFunctions()->Add(ptstats2);
    ptstats2->SetParent(h1_->GetListOfFunctions());
  }
  else
  {
    TPaveStats *ptstats = new TPaveStats(0.0,0.0,
                     0.0,0.0,"brNDC");
    ptstats->Draw();
    h0_->GetListOfFunctions()->Add(ptstats);
    ptstats->SetParent(h0_->GetListOfFunctions());
  }

  float min=-999.;
  float max=+999.;
  switch(mode) {
  case SCALE:
    h1_->Scale( h0_->GetEntries()/h1_->GetEntries() );
  case NORMAL:
    if(s0_)
      Styles::FormatHisto( h0_ , s0_);
    if(s1_)
      Styles::FormatHisto( h1_ , s1_);
      
    if( h1_->GetMaximum()>h0_->GetMaximum()) {
      h0_->SetMaximum( h1_->GetMaximum()*1.15 );
    }

    h0_->Draw();
    h1_->Draw("same");

    break;
  case EFF:
    if(s0_)
      Styles::FormatHisto( h0_ , s0_);
    if(s1_)
      Styles::FormatHisto( h1_ , s1_);

    // rather arbitrary but useful
    max= h0_->GetMaximum();
    if ( h1_->GetMaximum() > max )
      max = h1_->GetMaximum();
    if(max > 0.8) max = 1;
    
    max *= 1.1 ;

    min= h0_->GetMinimum();
    if ( h1_->GetMinimum() < min )
      min = h1_->GetMinimum();
    if(min > 0.2) min = 0.;

    min *= 0.8 ;
    h0_->SetMaximum( max );
    h0_->SetMinimum( min );

    h0_->Draw("E");
    h1_->Draw("Esame");
    break;
  case GRAPH:
    if(s0_)
      Styles::FormatHisto( h0_ , s0_);
    if(s1_)
      Styles::FormatHisto( h1_ , s1_);
      
    if( h1_->GetMaximum()>h0_->GetMaximum()) {
      h0_->SetMaximum( h1_->GetMaximum()*1.15 );
    }
    if( h1_->GetMinimum()<h0_->GetMinimum()) {
      h0_->SetMinimum( h1_->GetMinimum()*1.15 );
    }
    h0_->SetMarkerStyle(21);
    h0_->SetMarkerColor(2);
    h1_->SetMarkerStyle(21);
    h1_->SetMarkerColor(4);

    h0_->Draw("E1");
    h1_->Draw("E1same");
    break;
  case RATIO:
    h0_->Sumw2();
    h1_->Sumw2();
    h0_->Divide( h1_ );
    if(s0_)
      Styles::FormatHisto( h0_ , s0_);
    h0_->Draw();
  default:
    break;
  }
}