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();
  void loadHistograms(std::string const filename, std::vector<TH1F*> *hists);
  void drawComparison();

protected:   
  //Style Stuff
  void setHistStyle(TH1*, const char*, const char*);
  void setCanvasStyle( TCanvas& canv );
  void setLegendStyle( TLegend& leg );

 private:
  //define input/output
  std::string treename_, input_, reference_, output_;

  //define histogram design
  std::vector<short int> cmpColor_, cmpStyle_, cmpWidth_; 

  //define canvas design
  int gridX_, gridY_;

  //define legend design
  double legXLeft_,  legXRight_;
  double legYLower_, legYUpper_;

  //bins
  std::vector<int> Bins_;

  //hists
  std::vector<std::string> HistNames_;

  //validation hist
  std::vector<TH1F*> hists_, refs_;
};

void ValPlotMaker::setHistStyle( TH1* hist, const char* titleX, const char* titleY )
{
  hist->SetTitle( "" );
  hist->GetXaxis()->SetTitle( titleX );
  hist->GetXaxis()->SetTitleSize  ( 0.06 );
  hist->GetXaxis()->SetTitleColor (    1 );
  hist->GetXaxis()->SetTitleOffset(  1.0 );
  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_     = cfg.read<std::string>( "input_root" );
  reference_ = cfg.read<std::string>( "input_reference" );
  treename_  = cfg.read<std::string>( "tree_name" );
  output_    = cfg.read<std::string>( "output"   );
  Bins_      = bag_of<int>(  cfg.read<std::string>( "bins","0 1 2"));
  HistNames_ = bag_of_string(cfg.read<std::string>( "hists"));

  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" ));
}


void ValPlotMaker::loadHistograms(string const filename, vector<TH1F*> *hists)
{
  TFile * file_ = new TFile( filename.c_str() );
  TH1F *dummy;
  hists->clear();
  for (std::vector<std::string>::iterator it=HistNames_.begin(); 
       it!=HistNames_.end(); ) {
    
      string dummy_name = treename_+"/"+(*it);
      dummy = (TH1F*)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;
      }
  }
  cout << "...loaded " << hists->size() 
       <<" histograms from file '" << filename << "'." << endl;
}

void ValPlotMaker::loadHistograms()
{
  loadHistograms(input_, &hists_);
  loadHistograms(reference_, &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(3);     
  }
}

void ValPlotMaker::drawComparison()
{
  TPostScript ps(output_.c_str(), 112 );
  TCanvas *canv = new TCanvas("canv", "fit histograms", 600, 600);
  setCanvasStyle( *canv );
  canv->SetGridx( 1 );
  canv->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));
      setHistStyle( &refer, it->c_str(), "events"  );
      refer.SetTitle( "" );
      refer.SetLineWidth( 8 );
      refer.SetLineColor( 2 );
      refer.Sumw2();
      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<0.05){
        canv->SetFillColor(5);
        cout << "Deviation in plot '" << (*it) << "': " << ks 
             << " (Kolmogorov-Smirnov Test)" << endl;
      } 
      else {
        canv->SetFillColor(0);
      } 
        
      double maximum = (refer.GetMaximum()>hist.GetMaximum() ?
                        1.5*refer.GetMaximum() : 1.5*hist.GetMaximum());
      refer.SetMaximum(maximum);
      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, "data", "PE" );
      leg->AddEntry( &refer,"reference", "L" );
      leg->Draw("same");
      
      canv->Draw();
      ps.NewPage();
  }     

  ps.Close();
  canv->Close();
  delete canv;
}

int main(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( argc<2 ){
    cerr << "ERROR:" 
         << " Missing argument (the cfg file)!" << endl;
    return 1;
  }

  ValPlotMaker * plots = new ValPlotMaker();
  plots->readConfig( argv[1] );
  plots->loadHistograms();
  
  plots->drawComparison();
  delete plots;

  cout << "done." << endl; 
  return 0;
}
Revision:	1.1.1.1 (vendor branch)
Committed:	Mon Feb 25 15:54:04 2008 UTC (17 years, 2 months ago) by auterman
Content type:	text/plain
Branch:	tex, JetID, Demo
Changes since 1.1:	+0 -0 lines
Log Message:	PAT Jet Analyzer
#	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			public:
29			ValPlotMaker(){};
30			~ValPlotMaker(){};
31
32			void readConfig( std::string );
33			void loadHistograms();
34			void loadHistograms(std::string const filename, std::vector<TH1F> hists);
35			void drawComparison();
36
37			protected:
38			//Style Stuff
39			void setHistStyle(TH1, const char, const char*);
40			void setCanvasStyle( TCanvas& canv );
41			void setLegendStyle( TLegend& leg );
42
43			private:
44			//define input/output
45			std::string treename_, input_, reference_, output_;
46
47			//define histogram design
48			std::vector<short int> cmpColor_, cmpStyle_, cmpWidth_;
49
50			//define canvas design
51			int gridX_, gridY_;
52
53			//define legend design
54			double legXLeft_, legXRight_;
55			double legYLower_, legYUpper_;
56
57			//bins
58			std::vector<int> Bins_;
59
60			//hists
61			std::vector<std::string> HistNames_;
62
63			//validation hist
64			std::vector<TH1F*> hists_, refs_;
65			};
66
67			void ValPlotMaker::setHistStyle( TH1* hist, const char* titleX, const char* titleY )
68			{
69			hist->SetTitle( "" );
70			hist->GetXaxis()->SetTitle( titleX );
71			hist->GetXaxis()->SetTitleSize ( 0.06 );
72			hist->GetXaxis()->SetTitleColor ( 1 );
73			hist->GetXaxis()->SetTitleOffset( 1.0 );
74			hist->GetXaxis()->SetTitleFont ( 62 );
75			hist->GetXaxis()->SetLabelSize ( 0.05 );
76			hist->GetXaxis()->SetLabelFont ( 62 );
77			hist->GetXaxis()->CenterTitle ( );
78			hist->GetXaxis()->SetNdivisions ( 505 );
79			hist->GetYaxis()->SetTitle( titleY );
80			hist->GetYaxis()->SetTitleSize ( 0.07 );
81			hist->GetYaxis()->SetTitleColor ( 1 );
82			hist->GetYaxis()->SetTitleOffset( 1.3 );
83			hist->GetYaxis()->SetTitleFont ( 62 );
84			hist->GetYaxis()->SetLabelSize ( 0.05 );
85			hist->GetYaxis()->SetLabelFont ( 62 );
86			}
87			void ValPlotMaker::setCanvasStyle( TCanvas& canv )
88			{
89			canv.SetFillStyle ( 4000 );
90			canv.SetLeftMargin ( 0.20 );
91			canv.SetRightMargin ( 0.05 );
92			canv.SetBottomMargin( 0.15 );
93			canv.SetTopMargin ( 0.05 );
94			}
95			void ValPlotMaker::setLegendStyle( TLegend& leg )
96			{
97			leg.SetFillStyle ( 0 );
98			leg.SetBorderSize( 0 );
99			leg.SetFillColor ( 0 );
100			}
101
102			void ValPlotMaker::readConfig( std::string name )
103			{
104			ConfigFile cfg( name );
105			input_ = cfg.read<std::string>( "input_root" );
106			reference_ = cfg.read<std::string>( "input_reference" );
107			treename_ = cfg.read<std::string>( "tree_name" );
108			output_ = cfg.read<std::string>( "output" );
109			Bins_ = bag_of<int>( cfg.read<std::string>( "bins","0 1 2"));
110			HistNames_ = bag_of_string(cfg.read<std::string>( "hists"));
111
112			gridX_ = cfg.read<int> ( "grid_x" );
113			gridY_ = cfg.read<int> ( "grid_y" );
114			legXLeft_ = cfg.read<double> ( "leg_x_left" );
115			legXRight_= cfg.read<double> ( "leg_x_right" );
116			legYLower_= cfg.read<double> ( "leg_y_lower" );
117			legYUpper_= cfg.read<double> ( "leg_y_upper" );
118			cmpColor_ = bag_of<short int>( cfg.read<std::string>( "color" ));
119			cmpStyle_ = bag_of<short int>( cfg.read<std::string>( "style" ));
120			cmpWidth_ = bag_of<short int>( cfg.read<std::string>( "width" ));
121			}
122
123
124			void ValPlotMaker::loadHistograms(string const filename, vector<TH1F> hists)
125			{
126			TFile * file_ = new TFile( filename.c_str() );
127			TH1F *dummy;
128			hists->clear();
129			for (std::vector<std::string>::iterator it=HistNames_.begin();
130			it!=HistNames_.end(); ) {
131
132			string dummy_name = treename_+"/"+(*it);
133			dummy = (TH1F*)file_->Get( dummy_name.c_str() );
134			if( !dummy ){
135			cerr << "WARNING:"
136			<< " Didn't find indicated hist" << " [" << (*it) << "]"
137			<< " in tree" << " [" << treename_ << "]"
138			<< " of file" << " [" << filename << "]"
139			<< endl;
140			it = HistNames_.erase( it );
141			} else {
142			hists->push_back( dummy );
143			++it;
144			}
145			}
146			cout << "...loaded " << hists->size()
147			<<" histograms from file '" << filename << "'." << endl;
148			}
149
150			void ValPlotMaker::loadHistograms()
151			{
152			loadHistograms(input_, &hists_);
153			loadHistograms(reference_, &refs_);
154			if (hists_.size()!=refs_.size()) {
155			cerr << "ERROR: At least one histogram is only present in one input file!"
156			<< "The number of data and reference plots does not match." << endl;
157			exit(3);
158			}
159			}
160
161			void ValPlotMaker::drawComparison()
162			{
163			TPostScript ps(output_.c_str(), 112 );
164			TCanvas *canv = new TCanvas("canv", "fit histograms", 600, 600);
165			setCanvasStyle( *canv );
166			canv->SetGridx( 1 );
167			canv->SetGridy( 1 );
168
169			char * dummy_text = new char[100];
170			unsigned i=0;
171			double ks;
172			std::vector<std::string>::iterator it=HistNames_.begin();
173			std::vector<TH1F*>::iterator hst=hists_.begin();
174			std::vector<TH1F*>::iterator ref=refs_.begin();
175
176			for (;ref!=refs_.end()&&hst!=hists_.end()&&it!=HistNames_.end();
177			++ref,++hst,++it,++i){
178			//cout << i << "th plot: " << (*it) << endl;
179
180			TH1F& hist = ((TH1F)(*hst));
181			TH1F& refer = ((TH1F)(*ref));
182			setHistStyle( &refer, it->c_str(), "events" );
183			refer.SetTitle( "" );
184			refer.SetLineWidth( 8 );
185			refer.SetLineColor( 2 );
186			refer.Sumw2();
187			refer.Scale( hist.Integral()/refer.Integral() );
188			// Statistical test of compatibility in shape between
189			// THIS histogram and h2, using Kolmogorov test.
190			ks = refer.KolmogorovTest( &hist, "");
191			// Default: Ignore under- and overflow bins in comparison
192			// option is a character string to specify options
193			// "U" include Underflows in test (also for 2-dim)
194			// "O" include Overflows (also valid for 2-dim)
195			// "N" include comparison of normalizations
196			// "D" Put out a line of "Debug" printout
197			// "M" Return the Maximum Kolmogorov distance instead of prob
198			// "X" Run the pseudo experiments post-processor with the following procedure:
199			// make pseudoexperiments based on random values from the parent
200			// distribution and compare the KS distance of the pseudoexperiment
201			// to the parent distribution. Bin the KS distances in a histogram,
202			// and then take the integral of all the KS values above the value
203			// obtained from the original data to Monte Carlo distribution.
204			// The number of pseudo-experiments nEXPT is currently fixed at 1000.
205			// The function returns the integral.
206			// (thanks to Ben Kilminster to submit this procedure). Note that
207			// this option "X" is much slower.
208			if (ks<0.05){
209			canv->SetFillColor(5);
210			cout << "Deviation in plot '" << (*it) << "': " << ks
211			<< " (Kolmogorov-Smirnov Test)" << endl;
212			}
213			else {
214			canv->SetFillColor(0);
215			}
216
217			double maximum = (refer.GetMaximum()>hist.GetMaximum() ?
218			1.5refer.GetMaximum() : 1.5hist.GetMaximum());
219			refer.SetMaximum(maximum);
220			hist.SetMarkerStyle(8);
221			hist.SetLineWidth(3);
222			refer.Draw("h");
223			hist.Draw("pe,same");
224			TLegend* leg = new TLegend(0.5,0.67,0.95,0.87);
225			setLegendStyle( *leg );
226			sprintf(dummy_text,"KS = %2.2f",ks);
227			leg->SetHeader(dummy_text);
228			leg->AddEntry( &hist, "data", "PE" );
229			leg->AddEntry( &refer,"reference", "L" );
230			leg->Draw("same");
231
232			canv->Draw();
233			ps.NewPage();
234			}
235
236			ps.Close();
237			canv->Close();
238			delete canv;
239			}
240
241			int main(int argc, char* argv[])
242			{
243			setNiceStyle();
244			gStyle->SetOptStat( 0 );
245			gStyle->SetOptFit ( 0 );
246
247			//gStyle->SetStatColor(0);
248			gStyle->SetStatBorderSize(0);
249			gStyle->SetStatX(0.93);
250			gStyle->SetStatY(0.93);
251			gStyle->SetStatW(0.18);
252			gStyle->SetStatH(0.18);
253
254			if( argc<2 ){
255			cerr << "ERROR:"
256			<< " Missing argument (the cfg file)!" << endl;
257			return 1;
258			}
259
260			ValPlotMaker * plots = new ValPlotMaker();
261			plots->readConfig( argv[1] );
262			plots->loadHistograms();
263
264			plots->drawComparison();
265			delete plots;
266
267			cout << "done." << endl;
268			return 0;
269			}