PATJetIDAnalyzer/bin/ValidationPlotMaker.cc

//System
#include <memory>
#include <string>
#include <fstream>
#include <iostream>
//Root
#include <TF1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
#include <TTree.h>
#include <TPostScript.h>
#include <TMath.h>
#include <TMatrix.h>
#include <TFormula.h>
#include <TAxis.h>
#include <TAttLine.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TPaveText.h>
//User
#include "ConfigFile.h"
#include "NiceStyle.cc"

using namespace std;

class ValPlotMaker{
public:
  ValPlotMaker(){};
  ~ValPlotMaker(){};
  
  void readConfig( std::string );
  void loadHistograms();
  template <class T> void loadHistograms(std::vector<std::string> filename, vector<string> *histnames, std::vector<T*> *hists);
  void drawComparison();
  
protected:   
  bool contains(const string the_string, const string the_sub_string) const;
  //Style Stuff
  void setHistStyle(TH1*, const char*, const char*);
  void setCanvasStyle( TCanvas& canv );
  void setLegendStyle( TLegend& leg );
  
private:
  //define input/output
  std::string treename_, output_, name_data_, name_reference_;
  std::vector<std::string> input_, reference_;
  
  //define histogram design
  std::vector<short int> cmpColor_, cmpStyle_, cmpWidth_; 
  
  //define canvas design
  int gridX_, gridY_;
  bool scale_;
  
  //define legend design
  double legXLeft_,  legXRight_;
  double legYLower_, legYUpper_;
  
  //bins
  std::vector<int> Bins_;
  
  //hists
  std::vector<std::string> HistNames_, Hist2DNames_;
  
  //validation hist
  std::vector<TH1F*> hists_, refs_;
  std::vector<TH2F*> hists2D_, refs2D_;
  
  //background color of hists, if KS test fails
  double ks_treshold;
  int    ks_color;
  bool   UseTitle_;

};

bool ValPlotMaker::contains(const string the_string, const string the_sub_string) const
{
  return the_string.find(the_sub_string,0) != string::npos;
}

void ValPlotMaker::setHistStyle( TH1* hist, const char* titleX, const char* titleY )
{
  //hist->SetTitle( "" );
  if ( hist->GetXaxis()->GetTitle()=="" ) 
    hist->GetXaxis()->SetTitle( titleX );
  hist->GetXaxis()->SetTitleSize  ( 0.04 );
  hist->GetXaxis()->SetTitleColor (    1 );
  hist->GetXaxis()->SetTitleOffset(  1.8 );
  hist->GetXaxis()->SetTitleFont  (   62 );
  hist->GetXaxis()->SetLabelSize  ( 0.05 );
  hist->GetXaxis()->SetLabelFont  (   62 );
  hist->GetXaxis()->CenterTitle   (      );
  hist->GetXaxis()->SetNdivisions (  505 );
  //hist->GetYaxis()->SetTitle( titleY );
  hist->GetYaxis()->SetTitleSize  ( 0.07 );
  hist->GetYaxis()->SetTitleColor (    1 );
  hist->GetYaxis()->SetTitleOffset(  1.3 );
  hist->GetYaxis()->SetTitleFont  (   62 );
  hist->GetYaxis()->SetLabelSize  ( 0.05 );
  hist->GetYaxis()->SetLabelFont  (   62 );
}
void ValPlotMaker::setCanvasStyle( TCanvas& canv )
{
  canv.SetFillStyle   ( 4000 );
  canv.SetLeftMargin  ( 0.20 );
  canv.SetRightMargin ( 0.05 );
  canv.SetBottomMargin( 0.15 );
  //canv.SetTopMargin   ( 0.05 );
}
void ValPlotMaker::setLegendStyle( TLegend& leg )
{
  leg.SetFillStyle ( 0 );
  leg.SetBorderSize( 0 );
  leg.SetFillColor ( 0 ); 
}

void ValPlotMaker::readConfig( std::string name )
{
  ConfigFile cfg( name );  
  input_       = bag_of_string(cfg.read<std::string>( "input_root" ));
  reference_   = bag_of_string(cfg.read<std::string>( "input_reference" ));
  name_data_      = cfg.read<std::string>( "name data", "data" );
  name_reference_ = cfg.read<std::string>( "name reference", "reference" );
  treename_    = cfg.read<std::string>( "tree_name" );
  output_      = cfg.read<std::string>( "output"   );
  HistNames_   = bag_of_string(cfg.read<std::string>( "hists"));
  Hist2DNames_ = bag_of_string(cfg.read<std::string>( "2D hists"));
  ks_treshold  = cfg.read<double> ( "KS treshold", 0.05  );
  ks_color     = cfg.read<int>    ( "KS color", 5     );
  gridX_       = cfg.read<int>    ( "grid_x"     );
  gridY_       = cfg.read<int>    ( "grid_y"     );
  legXLeft_    = cfg.read<double> ( "leg_x_left"  );
  legXRight_   = cfg.read<double> ( "leg_x_right" );
  legYLower_   = cfg.read<double> ( "leg_y_lower" );
  legYUpper_   = cfg.read<double> ( "leg_y_upper" );
  cmpColor_    = bag_of<short int>( cfg.read<std::string>( "color" ));
  cmpStyle_    = bag_of<short int>( cfg.read<std::string>( "style" ));
  cmpWidth_    = bag_of<short int>( cfg.read<std::string>( "width" ));
  UseTitle_    = cfg.read<int>    ( "show title", false  );
  scale_       = cfg.read<bool>   ( "scale", false       );
}


template <class T> void ValPlotMaker::loadHistograms(vector<string> filenames, vector<string> *histnames, vector<T*> *hists)
{
  std::vector<std::string>::iterator filename = filenames.begin();
  TFile * file_ = new TFile( (*filename).c_str() );
  T * dummy;
  hists->clear();
  for (std::vector<std::string>::iterator it=histnames->begin(); 
       it!=histnames->end(); ) {
    string dummy_name = treename_+"/"+(*it);
    dummy = (T*)file_->Get( dummy_name.c_str() );
    if( !dummy ){
      cerr << "WARNING:" 
           << " Didn't find indicated hist" << " [" << (*it) << "]" 
           << " in tree" << " [" << treename_ << "]" 
           << " of file" << " [" << (*filename) << "]" 
           << endl;
      it = histnames->erase( it ); 
    } else {
      hists->push_back( dummy );
      ++it;
    }
  }
  ++filename;
  for (; filename!=filenames.end(); ++filename){
    TFile * file_ = new TFile( (*filename).c_str() );
    int i=0;
    typename std::vector<T*>::iterator hist=hists->begin();
    for (std::vector<std::string>::iterator it=histnames->begin(); 
         it!=histnames->end()&&hist!=hists->end(); ++hist,++it,++i) {
      string dummy_name = treename_+"/"+(*it);
      dummy = (T*)file_->Get( dummy_name.c_str() );
      if( dummy ){
        dummy->Add( (*hist) );
        hist = hists->erase(hist);
        hists->insert( hist, 1, dummy );
      }
    }
  }
  
  if (hists->size() > 0) {
    filename = filenames.begin();
    cout << "...loaded " << hists->size() 
         <<" histograms from file(s) '" << (*filename) << "'." << endl;
  }     
}

void ValPlotMaker::loadHistograms()
{
  loadHistograms<TH1F>(input_,     &HistNames_, &hists_);
  loadHistograms<TH1F>(reference_, &HistNames_, &refs_);
  if (hists_.size()!=refs_.size()) {
    cerr << "ERROR: At least one histogram is only present in one input file!"
         << "The number of data and reference plots does not match." << endl;
    exit(2);     
  }
  
  loadHistograms<TH2F>(input_,     &Hist2DNames_, &hists2D_);
  loadHistograms<TH2F>(reference_, &Hist2DNames_, &refs2D_);
  if (hists2D_.size()!=refs2D_.size()) {
    cerr << "ERROR: At least one 2D-histogram is only present in one input file!"
         << "The number of data and reference plots does not match." << endl;
    exit(3);     
  }
  
}

void ValPlotMaker::drawComparison()
{
  TPostScript ps(output_.c_str(), 112 );
  TCanvas *c1 = new TCanvas("c1", "fit histograms", 600, 600);
  setCanvasStyle( *c1 );
  c1->SetGridx( 1 );
  c1->SetGridy( 1 );
  
  char * dummy_text = new char[100];
  unsigned i=0;
  double ks;
  std::vector<std::string>::iterator it=HistNames_.begin();
  std::vector<TH1F*>::iterator hst=hists_.begin();
  std::vector<TH1F*>::iterator ref=refs_.begin();
  
  for (;ref!=refs_.end()&&hst!=hists_.end()&&it!=HistNames_.end(); 
       ++ref,++hst,++it,++i){
    //cout << i << "th plot: " << (*it) << endl;
    
    TH1F& hist   = *((TH1F*)(*hst));
    TH1F& refer  = *((TH1F*)(*ref));
    //check if hist is filled
    if (refer.Integral()==0.0) {
      cerr << "Histogram '"<< (*it) <<"' is empty!" <<endl;
      continue;
    }   
    
    setHistStyle( &refer, it->c_str(), "events"  );
    if (UseTitle_){
      //refer.SetTitle( "" );
    }   
    else
      refer.SetTitle( "" );
    refer.SetLineWidth( 8 );
    refer.SetLineColor( 2 );
    
    cout << i << "th plot: " << (*it) 
         << " has " << refer.Integral() << " 'reference' events and "
         <<  hist.Integral() << " 'data' events" << endl;
    
    refer.Sumw2();
    if ( scale_ && !contains(refer.GetYaxis()->GetTitle(), "eff") ) {
      refer.Scale( hist.Integral()/refer.Integral() );
    }
//  Statistical test of compatibility in shape between
//  THIS histogram and h2, using Kolmogorov test.
    ks = refer.KolmogorovTest( &hist, "");
//     Default: Ignore under- and overflow bins in comparison
//     option is a character string to specify options
//         "U" include Underflows in test  (also for 2-dim)
//         "O" include Overflows     (also valid for 2-dim)
//         "N" include comparison of normalizations
//         "D" Put out a line of "Debug" printout
//         "M" Return the Maximum Kolmogorov distance instead of prob
//         "X" Run the pseudo experiments post-processor with the following procedure:
//             make pseudoexperiments based on random values from the parent
//             distribution and compare the KS distance of the pseudoexperiment
//             to the parent distribution. Bin the KS distances in a histogram,
//             and then take the integral of all the KS values above the value
//             obtained from the original data to Monte Carlo distribution.
//             The number of pseudo-experiments nEXPT is currently fixed at 1000.
//             The function returns the integral.
//             (thanks to Ben Kilminster to submit this procedure). Note that
//             this option "X" is much slower.
    if (ks<ks_treshold){
      c1->SetFillColor(ks_color);
      cout << "Deviation in plot "<< i <<" '" << (*it) << "': " << ks 
           << " (Kolmogorov-Smirnov Test)" << endl;
    }   
    else {
      c1->SetFillColor(0);
    }   
    
    //Show under- and over-flow bins
    refer.SetBinContent(1,refer.GetBinContent(0)+refer.GetBinContent(1));
    hist.SetBinContent(1,hist.GetBinContent(0)+hist.GetBinContent(1));
    unsigned int mb = refer.GetNbinsX();
    refer.SetBinContent(mb-1,refer.GetBinContent(mb)+refer.GetBinContent(mb-1));
    mb = hist.GetNbinsX();
    hist.SetBinContent(mb-1,hist.GetBinContent(mb)+hist.GetBinContent(mb-1));
    
    double maximum = (refer.GetMaximum()>hist.GetMaximum() ?
                      1.5*refer.GetMaximum() : 1.5*hist.GetMaximum());
    double minimum = (refer.GetMinimum()<hist.GetMinimum() ?
                      refer.GetMinimum() : hist.GetMinimum());
    refer.SetMaximum(maximum);
    refer.SetMinimum(minimum);
    for (unsigned int b=0; b<mb; ++b)
      if (refer.GetBinContent(b)<=0.)
        refer.SetBinError(b, sqrt(-1.*refer.GetBinContent(b)));
    
    
    hist.SetMarkerStyle(8);
    hist.SetLineWidth(3);
    refer.Draw("h");
    hist.Draw("pe,same");
    TLegend* leg = new TLegend(0.5,0.67,0.95,0.87);
    setLegendStyle( *leg );
    sprintf(dummy_text,"KS = %2.2f",ks);
    leg->SetHeader(dummy_text);
    leg->AddEntry( &hist, name_data_.c_str(), "PE" );
    leg->AddEntry( &refer,name_reference_.c_str(), "L" );
    leg->Draw("same");
    
    c1->Draw();
    ps.NewPage();
  }  
  
  
  //2D-stuff 
  
  i=0;
  it=Hist2DNames_.begin();
  std::vector<TH2F*>::iterator hst2D=hists2D_.begin();
  std::vector<TH2F*>::iterator ref2D=refs2D_.begin();
  
  for (;ref2D!=refs2D_.end()&&hst2D!=hists2D_.end()&&it!=Hist2DNames_.end(); 
       ++ref2D,++hst2D,++it,++i){
    cout << i << "th 2D-plot: " << (*it) << endl;
    
    TH1F& hist2D   = *((TH1F*)(*hst2D));
    TH1F& refer2D  = *((TH1F*)(*ref2D));
    TH2F * diff2D  = new TH2F( *((TH2F*)(*hst2D)) );
    setHistStyle( &refer2D, it->c_str(), "events"  );

    for (int binx=0; binx<=refer2D.GetNbinsX();++binx) {
    for (int biny=0; biny<=refer2D.GetNbinsY();++biny) {
    //  if (refer2D.GetBinContent(b)<=0.)
    //    refer2D.SetBinContent(b, fabs(refer2D.GetBinContent(b)));
    //  if (hist2D.GetBinContent(b)<=0.)
    //    hist2D.SetBinContent(b, fabs(hist2D.GetBinContent(b)));
        diff2D->SetBinContent(binx,biny,
                refer2D.GetBinContent(binx,biny)-hist2D.GetBinContent(binx,biny) );
    }}
    if (UseTitle_){
      //refer2D.SetTitle( "" );
    }   
    else
      refer2D.SetTitle( "" );
    //refer2D.SetLineWidth( 8 );
    refer2D.SetLineColor( 2 );
    refer2D.Sumw2();
    if ( !contains(refer2D.GetYaxis()->GetTitle(), "eff") ) {
      refer2D.Scale( hist2D.Integral()/refer2D.Integral() );
    }
    ks = refer2D.KolmogorovTest( &hist2D, "");
    if (ks<ks_treshold){
      c1->SetFillColor(ks_color);
      cout << "Deviation in plot "<< i <<" '" << (*it) << "': " << ks 
           << " (Kolmogorov-Smirnov Test)" << endl;
    }   
    else {
      c1->SetFillColor(0);
    }   

    double maximum = (refer2D.GetMaximum()>hist2D.GetMaximum() ?
                      1.5*refer2D.GetMaximum() : 1.5*hist2D.GetMaximum());
    double minimum = (refer2D.GetMinimum()<hist2D.GetMinimum() ?
                      refer2D.GetMinimum() : hist2D.GetMinimum());
    //refer2D.SetMaximum(maximum);
    //refer2D.SetMinimum(minimum);
    
    refer2D.Draw("BOX");
    hist2D.Draw("BOX,same");
    
    TLegend* leg = new TLegend(0.7,0.75,0.95,0.95);
    setLegendStyle( *leg );
    sprintf(dummy_text,"KS = %2.2f",ks);
    leg->SetHeader(dummy_text);
    leg->AddEntry( &hist2D, "data", "F" );
    leg->AddEntry( &refer2D,"reference", "F" );
    leg->Draw("same");
    c1->Draw();
    ps.NewPage();


    diff2D->Draw("BOX");
    
    TLegend* dleg = new TLegend(0.7,0.8,0.95,0.95);
    setLegendStyle( *dleg );
    sprintf(dummy_text,"Difference:",ks);
    dleg->SetHeader(dummy_text);
    dleg->AddEntry( diff2D,"reference - data", "F" );
    dleg->Draw("same");
    c1->Draw();
    ps.NewPage();
    delete diff2D;
  }
  
  ps.Close();
  c1->Close();
  delete c1;
}


int ValidationPlotMaker(int argc, char* argv[])
{
  setNiceStyle();
  gStyle->SetOptStat( 0 );
  gStyle->SetOptFit ( 0 );
  gStyle->SetStatColor(0);
  gStyle->SetStatBorderSize(0);
  gStyle->SetStatX(0.93);
  gStyle->SetStatY(0.93);
  gStyle->SetStatW(0.18);
  gStyle->SetStatH(0.18);
  
  //If no argument is given, use default config file name
  string config_file = "bin/ValidationPlot.cfg";
  if( argc>=2 )
    config_file = argv[1];
  
  //Run the Validation Plot Maker
  ValPlotMaker * plots = new ValPlotMaker();
  plots->readConfig( config_file.c_str() );
  plots->loadHistograms();
  plots->drawComparison();
  delete plots;
  
  return 0;
}

int main(int argc, char* argv[])
{
  return ValidationPlotMaker(argc, argv);
}
Revision:	1.1.1.2 (vendor branch)
Committed:	Mon Apr 14 12:45:47 2008 UTC (17 years ago) by auterman
Content type:	text/plain
Branch:	tex, JetID, Demo
CVS Tags:	start
Changes since 1.1.1.1:	+283 -115 lines
Log Message:	A JetID skeleton
#	User	Rev	Content
1	auterman	1.1	//System
2			#include <memory>
3			#include <string>
4			#include <fstream>
5			#include <iostream>
6			//Root
7			#include <TF1.h>
8			#include <TH1F.h>
9			#include <TH2F.h>
10			#include <TFile.h>
11			#include <TTree.h>
12			#include <TPostScript.h>
13			#include <TMath.h>
14			#include <TMatrix.h>
15			#include <TFormula.h>
16			#include <TAxis.h>
17			#include <TAttLine.h>
18			#include <TCanvas.h>
19			#include <TLegend.h>
20			#include <TPaveText.h>
21			//User
22			#include "ConfigFile.h"
23			#include "NiceStyle.cc"
24
25			using namespace std;
26
27			class ValPlotMaker{
28	auterman	1.1.1.2	public:
29	auterman	1.1	ValPlotMaker(){};
30			~ValPlotMaker(){};
31	auterman	1.1.1.2
32	auterman	1.1	void readConfig( std::string );
33			void loadHistograms();
34	auterman	1.1.1.2	template <class T> void loadHistograms(std::vector<std::string> filename, vector<string> histnames, std::vector<T> *hists);
35	auterman	1.1	void drawComparison();
36	auterman	1.1.1.2
37	auterman	1.1	protected:
38	auterman	1.1.1.2	bool contains(const string the_string, const string the_sub_string) const;
39	auterman	1.1	//Style Stuff
40			void setHistStyle(TH1, const char, const char*);
41			void setCanvasStyle( TCanvas& canv );
42			void setLegendStyle( TLegend& leg );
43	auterman	1.1.1.2
44			private:
45	auterman	1.1	//define input/output
46	auterman	1.1.1.2	std::string treename_, output_, name_data_, name_reference_;
47			std::vector<std::string> input_, reference_;
48
49	auterman	1.1	//define histogram design
50			std::vector<short int> cmpColor_, cmpStyle_, cmpWidth_;
51	auterman	1.1.1.2
52	auterman	1.1	//define canvas design
53			int gridX_, gridY_;
54	auterman	1.1.1.2	bool scale_;
55
56	auterman	1.1	//define legend design
57			double legXLeft_, legXRight_;
58			double legYLower_, legYUpper_;
59	auterman	1.1.1.2
60	auterman	1.1	//bins
61			std::vector<int> Bins_;
62	auterman	1.1.1.2
63	auterman	1.1	//hists
64	auterman	1.1.1.2	std::vector<std::string> HistNames_, Hist2DNames_;
65
66	auterman	1.1	//validation hist
67			std::vector<TH1F*> hists_, refs_;
68	auterman	1.1.1.2	std::vector<TH2F*> hists2D_, refs2D_;
69
70			//background color of hists, if KS test fails
71			double ks_treshold;
72			int ks_color;
73			bool UseTitle_;
74
75	auterman	1.1	};
76
77	auterman	1.1.1.2	bool ValPlotMaker::contains(const string the_string, const string the_sub_string) const
78			{
79			return the_string.find(the_sub_string,0) != string::npos;
80			}
81
82	auterman	1.1	void ValPlotMaker::setHistStyle( TH1* hist, const char* titleX, const char* titleY )
83			{
84	auterman	1.1.1.2	//hist->SetTitle( "" );
85			if ( hist->GetXaxis()->GetTitle()=="" )
86			hist->GetXaxis()->SetTitle( titleX );
87			hist->GetXaxis()->SetTitleSize ( 0.04 );
88	auterman	1.1	hist->GetXaxis()->SetTitleColor ( 1 );
89	auterman	1.1.1.2	hist->GetXaxis()->SetTitleOffset( 1.8 );
90	auterman	1.1	hist->GetXaxis()->SetTitleFont ( 62 );
91			hist->GetXaxis()->SetLabelSize ( 0.05 );
92			hist->GetXaxis()->SetLabelFont ( 62 );
93			hist->GetXaxis()->CenterTitle ( );
94			hist->GetXaxis()->SetNdivisions ( 505 );
95	auterman	1.1.1.2	//hist->GetYaxis()->SetTitle( titleY );
96	auterman	1.1	hist->GetYaxis()->SetTitleSize ( 0.07 );
97			hist->GetYaxis()->SetTitleColor ( 1 );
98			hist->GetYaxis()->SetTitleOffset( 1.3 );
99			hist->GetYaxis()->SetTitleFont ( 62 );
100			hist->GetYaxis()->SetLabelSize ( 0.05 );
101			hist->GetYaxis()->SetLabelFont ( 62 );
102			}
103			void ValPlotMaker::setCanvasStyle( TCanvas& canv )
104			{
105			canv.SetFillStyle ( 4000 );
106			canv.SetLeftMargin ( 0.20 );
107			canv.SetRightMargin ( 0.05 );
108			canv.SetBottomMargin( 0.15 );
109	auterman	1.1.1.2	//canv.SetTopMargin ( 0.05 );
110	auterman	1.1	}
111			void ValPlotMaker::setLegendStyle( TLegend& leg )
112			{
113			leg.SetFillStyle ( 0 );
114			leg.SetBorderSize( 0 );
115			leg.SetFillColor ( 0 );
116			}
117
118			void ValPlotMaker::readConfig( std::string name )
119			{
120			ConfigFile cfg( name );
121	auterman	1.1.1.2	input_ = bag_of_string(cfg.read<std::string>( "input_root" ));
122			reference_ = bag_of_string(cfg.read<std::string>( "input_reference" ));
123			name_data_ = cfg.read<std::string>( "name data", "data" );
124			name_reference_ = cfg.read<std::string>( "name reference", "reference" );
125			treename_ = cfg.read<std::string>( "tree_name" );
126			output_ = cfg.read<std::string>( "output" );
127			HistNames_ = bag_of_string(cfg.read<std::string>( "hists"));
128			Hist2DNames_ = bag_of_string(cfg.read<std::string>( "2D hists"));
129			ks_treshold = cfg.read<double> ( "KS treshold", 0.05 );
130			ks_color = cfg.read<int> ( "KS color", 5 );
131			gridX_ = cfg.read<int> ( "grid_x" );
132			gridY_ = cfg.read<int> ( "grid_y" );
133			legXLeft_ = cfg.read<double> ( "leg_x_left" );
134			legXRight_ = cfg.read<double> ( "leg_x_right" );
135			legYLower_ = cfg.read<double> ( "leg_y_lower" );
136			legYUpper_ = cfg.read<double> ( "leg_y_upper" );
137			cmpColor_ = bag_of<short int>( cfg.read<std::string>( "color" ));
138			cmpStyle_ = bag_of<short int>( cfg.read<std::string>( "style" ));
139			cmpWidth_ = bag_of<short int>( cfg.read<std::string>( "width" ));
140			UseTitle_ = cfg.read<int> ( "show title", false );
141			scale_ = cfg.read<bool> ( "scale", false );
142	auterman	1.1	}
143
144
145	auterman	1.1.1.2	template <class T> void ValPlotMaker::loadHistograms(vector<string> filenames, vector<string> histnames, vector<T> *hists)
146	auterman	1.1	{
147	auterman	1.1.1.2	std::vector<std::string>::iterator filename = filenames.begin();
148			TFile * file_ = new TFile( (*filename).c_str() );
149			T * dummy;
150	auterman	1.1	hists->clear();
151	auterman	1.1.1.2	for (std::vector<std::string>::iterator it=histnames->begin();
152			it!=histnames->end(); ) {
153			string dummy_name = treename_+"/"+(*it);
154			dummy = (T*)file_->Get( dummy_name.c_str() );
155			if( !dummy ){
156			cerr << "WARNING:"
157			<< " Didn't find indicated hist" << " [" << (*it) << "]"
158			<< " in tree" << " [" << treename_ << "]"
159			<< " of file" << " [" << (*filename) << "]"
160			<< endl;
161			it = histnames->erase( it );
162			} else {
163			hists->push_back( dummy );
164			++it;
165			}
166			}
167			++filename;
168			for (; filename!=filenames.end(); ++filename){
169			TFile * file_ = new TFile( (*filename).c_str() );
170			int i=0;
171			typename std::vector<T*>::iterator hist=hists->begin();
172			for (std::vector<std::string>::iterator it=histnames->begin();
173			it!=histnames->end()&&hist!=hists->end(); ++hist,++it,++i) {
174	auterman	1.1	string dummy_name = treename_+"/"+(*it);
175	auterman	1.1.1.2	dummy = (T*)file_->Get( dummy_name.c_str() );
176			if( dummy ){
177			dummy->Add( (*hist) );
178			hist = hists->erase(hist);
179			hists->insert( hist, 1, dummy );
180	auterman	1.1	}
181	auterman	1.1.1.2	}
182	auterman	1.1	}
183	auterman	1.1.1.2
184			if (hists->size() > 0) {
185			filename = filenames.begin();
186			cout << "...loaded " << hists->size()
187			<<" histograms from file(s) '" << (*filename) << "'." << endl;
188			}
189	auterman	1.1	}
190
191			void ValPlotMaker::loadHistograms()
192			{
193	auterman	1.1.1.2	loadHistograms<TH1F>(input_, &HistNames_, &hists_);
194			loadHistograms<TH1F>(reference_, &HistNames_, &refs_);
195	auterman	1.1	if (hists_.size()!=refs_.size()) {
196			cerr << "ERROR: At least one histogram is only present in one input file!"
197			<< "The number of data and reference plots does not match." << endl;
198	auterman	1.1.1.2	exit(2);
199			}
200
201			loadHistograms<TH2F>(input_, &Hist2DNames_, &hists2D_);
202			loadHistograms<TH2F>(reference_, &Hist2DNames_, &refs2D_);
203			if (hists2D_.size()!=refs2D_.size()) {
204			cerr << "ERROR: At least one 2D-histogram is only present in one input file!"
205			<< "The number of data and reference plots does not match." << endl;
206	auterman	1.1	exit(3);
207			}
208	auterman	1.1.1.2
209	auterman	1.1	}
210
211			void ValPlotMaker::drawComparison()
212			{
213			TPostScript ps(output_.c_str(), 112 );
214	auterman	1.1.1.2	TCanvas *c1 = new TCanvas("c1", "fit histograms", 600, 600);
215			setCanvasStyle( *c1 );
216			c1->SetGridx( 1 );
217			c1->SetGridy( 1 );
218
219	auterman	1.1	char * dummy_text = new char[100];
220			unsigned i=0;
221			double ks;
222			std::vector<std::string>::iterator it=HistNames_.begin();
223			std::vector<TH1F*>::iterator hst=hists_.begin();
224			std::vector<TH1F*>::iterator ref=refs_.begin();
225
226			for (;ref!=refs_.end()&&hst!=hists_.end()&&it!=HistNames_.end();
227	auterman	1.1.1.2	++ref,++hst,++it,++i){
228			//cout << i << "th plot: " << (*it) << endl;
229
230			TH1F& hist = ((TH1F)(*hst));
231			TH1F& refer = ((TH1F)(*ref));
232			//check if hist is filled
233			if (refer.Integral()==0.0) {
234			cerr << "Histogram '"<< (*it) <<"' is empty!" <<endl;
235			continue;
236			}
237
238			setHistStyle( &refer, it->c_str(), "events" );
239			if (UseTitle_){
240			//refer.SetTitle( "" );
241			}
242			else
243	auterman	1.1	refer.SetTitle( "" );
244	auterman	1.1.1.2	refer.SetLineWidth( 8 );
245			refer.SetLineColor( 2 );
246
247			cout << i << "th plot: " << (*it)
248			<< " has " << refer.Integral() << " 'reference' events and "
249			<< hist.Integral() << " 'data' events" << endl;
250
251			refer.Sumw2();
252			if ( scale_ && !contains(refer.GetYaxis()->GetTitle(), "eff") ) {
253	auterman	1.1	refer.Scale( hist.Integral()/refer.Integral() );
254	auterman	1.1.1.2	}
255	auterman	1.1	// Statistical test of compatibility in shape between
256			// THIS histogram and h2, using Kolmogorov test.
257	auterman	1.1.1.2	ks = refer.KolmogorovTest( &hist, "");
258	auterman	1.1	// Default: Ignore under- and overflow bins in comparison
259			// option is a character string to specify options
260			// "U" include Underflows in test (also for 2-dim)
261			// "O" include Overflows (also valid for 2-dim)
262			// "N" include comparison of normalizations
263			// "D" Put out a line of "Debug" printout
264			// "M" Return the Maximum Kolmogorov distance instead of prob
265			// "X" Run the pseudo experiments post-processor with the following procedure:
266			// make pseudoexperiments based on random values from the parent
267			// distribution and compare the KS distance of the pseudoexperiment
268			// to the parent distribution. Bin the KS distances in a histogram,
269			// and then take the integral of all the KS values above the value
270			// obtained from the original data to Monte Carlo distribution.
271			// The number of pseudo-experiments nEXPT is currently fixed at 1000.
272			// The function returns the integral.
273			// (thanks to Ben Kilminster to submit this procedure). Note that
274			// this option "X" is much slower.
275	auterman	1.1.1.2	if (ks<ks_treshold){
276			c1->SetFillColor(ks_color);
277			cout << "Deviation in plot "<< i <<" '" << (*it) << "': " << ks
278			<< " (Kolmogorov-Smirnov Test)" << endl;
279			}
280			else {
281			c1->SetFillColor(0);
282			}
283
284			//Show under- and over-flow bins
285			refer.SetBinContent(1,refer.GetBinContent(0)+refer.GetBinContent(1));
286			hist.SetBinContent(1,hist.GetBinContent(0)+hist.GetBinContent(1));
287			unsigned int mb = refer.GetNbinsX();
288			refer.SetBinContent(mb-1,refer.GetBinContent(mb)+refer.GetBinContent(mb-1));
289			mb = hist.GetNbinsX();
290			hist.SetBinContent(mb-1,hist.GetBinContent(mb)+hist.GetBinContent(mb-1));
291
292			double maximum = (refer.GetMaximum()>hist.GetMaximum() ?
293			1.5refer.GetMaximum() : 1.5hist.GetMaximum());
294			double minimum = (refer.GetMinimum()<hist.GetMinimum() ?
295			refer.GetMinimum() : hist.GetMinimum());
296			refer.SetMaximum(maximum);
297			refer.SetMinimum(minimum);
298			for (unsigned int b=0; b<mb; ++b)
299			if (refer.GetBinContent(b)<=0.)
300			refer.SetBinError(b, sqrt(-1.*refer.GetBinContent(b)));
301
302
303			hist.SetMarkerStyle(8);
304			hist.SetLineWidth(3);
305			refer.Draw("h");
306			hist.Draw("pe,same");
307			TLegend* leg = new TLegend(0.5,0.67,0.95,0.87);
308			setLegendStyle( *leg );
309			sprintf(dummy_text,"KS = %2.2f",ks);
310			leg->SetHeader(dummy_text);
311			leg->AddEntry( &hist, name_data_.c_str(), "PE" );
312			leg->AddEntry( &refer,name_reference_.c_str(), "L" );
313			leg->Draw("same");
314
315			c1->Draw();
316			ps.NewPage();
317			}
318
319
320
321			//2D-stuff
322
323			i=0;
324			it=Hist2DNames_.begin();
325			std::vector<TH2F*>::iterator hst2D=hists2D_.begin();
326			std::vector<TH2F*>::iterator ref2D=refs2D_.begin();
327
328			for (;ref2D!=refs2D_.end()&&hst2D!=hists2D_.end()&&it!=Hist2DNames_.end();
329			++ref2D,++hst2D,++it,++i){
330			cout << i << "th 2D-plot: " << (*it) << endl;
331
332			TH1F& hist2D = ((TH1F)(*hst2D));
333			TH1F& refer2D = ((TH1F)(*ref2D));
334			TH2F * diff2D = new TH2F( ((TH2F)(*hst2D)) );
335			setHistStyle( &refer2D, it->c_str(), "events" );
336
337			for (int binx=0; binx<=refer2D.GetNbinsX();++binx) {
338			for (int biny=0; biny<=refer2D.GetNbinsY();++biny) {
339			// if (refer2D.GetBinContent(b)<=0.)
340			// refer2D.SetBinContent(b, fabs(refer2D.GetBinContent(b)));
341			// if (hist2D.GetBinContent(b)<=0.)
342			// hist2D.SetBinContent(b, fabs(hist2D.GetBinContent(b)));
343			diff2D->SetBinContent(binx,biny,
344			refer2D.GetBinContent(binx,biny)-hist2D.GetBinContent(binx,biny) );
345			}}
346			if (UseTitle_){
347			//refer2D.SetTitle( "" );
348			}
349			else
350			refer2D.SetTitle( "" );
351			//refer2D.SetLineWidth( 8 );
352			refer2D.SetLineColor( 2 );
353			refer2D.Sumw2();
354			if ( !contains(refer2D.GetYaxis()->GetTitle(), "eff") ) {
355			refer2D.Scale( hist2D.Integral()/refer2D.Integral() );
356			}
357			ks = refer2D.KolmogorovTest( &hist2D, "");
358			if (ks<ks_treshold){
359			c1->SetFillColor(ks_color);
360			cout << "Deviation in plot "<< i <<" '" << (*it) << "': " << ks
361			<< " (Kolmogorov-Smirnov Test)" << endl;
362			}
363			else {
364			c1->SetFillColor(0);
365			}
366
367			double maximum = (refer2D.GetMaximum()>hist2D.GetMaximum() ?
368			1.5refer2D.GetMaximum() : 1.5hist2D.GetMaximum());
369			double minimum = (refer2D.GetMinimum()<hist2D.GetMinimum() ?
370			refer2D.GetMinimum() : hist2D.GetMinimum());
371			//refer2D.SetMaximum(maximum);
372			//refer2D.SetMinimum(minimum);
373
374			refer2D.Draw("BOX");
375			hist2D.Draw("BOX,same");
376
377			TLegend* leg = new TLegend(0.7,0.75,0.95,0.95);
378			setLegendStyle( *leg );
379			sprintf(dummy_text,"KS = %2.2f",ks);
380			leg->SetHeader(dummy_text);
381			leg->AddEntry( &hist2D, "data", "F" );
382			leg->AddEntry( &refer2D,"reference", "F" );
383			leg->Draw("same");
384			c1->Draw();
385			ps.NewPage();
386	auterman	1.1
387	auterman	1.1.1.2
388			diff2D->Draw("BOX");
389
390			TLegend* dleg = new TLegend(0.7,0.8,0.95,0.95);
391			setLegendStyle( *dleg );
392			sprintf(dummy_text,"Difference:",ks);
393			dleg->SetHeader(dummy_text);
394			dleg->AddEntry( diff2D,"reference - data", "F" );
395			dleg->Draw("same");
396			c1->Draw();
397			ps.NewPage();
398			delete diff2D;
399			}
400
401	auterman	1.1	ps.Close();
402	auterman	1.1.1.2	c1->Close();
403			delete c1;
404	auterman	1.1	}
405
406	auterman	1.1.1.2
407			int ValidationPlotMaker(int argc, char* argv[])
408	auterman	1.1	{
409			setNiceStyle();
410			gStyle->SetOptStat( 0 );
411			gStyle->SetOptFit ( 0 );
412	auterman	1.1.1.2	gStyle->SetStatColor(0);
413	auterman	1.1	gStyle->SetStatBorderSize(0);
414			gStyle->SetStatX(0.93);
415			gStyle->SetStatY(0.93);
416			gStyle->SetStatW(0.18);
417			gStyle->SetStatH(0.18);
418	auterman	1.1.1.2
419			//If no argument is given, use default config file name
420			string config_file = "bin/ValidationPlot.cfg";
421			if( argc>=2 )
422			config_file = argv[1];
423
424			//Run the Validation Plot Maker
425	auterman	1.1	ValPlotMaker * plots = new ValPlotMaker();
426	auterman	1.1.1.2	plots->readConfig( config_file.c_str() );
427	auterman	1.1	plots->loadHistograms();
428			plots->drawComparison();
429			delete plots;
430	auterman	1.1.1.2
431	auterman	1.1	return 0;
432			}
433	auterman	1.1.1.2
434			int main(int argc, char* argv[])
435			{
436			return ValidationPlotMaker(argc, argv);
437			}