| 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 |
template <class T> void loadHistograms(std::vector<std::string> filename, vector<string> *histnames, std::vector<T*> *hists);
|
| 35 |
void drawComparison();
|
| 36 |
|
| 37 |
protected:
|
| 38 |
bool contains(const string the_string, const string the_sub_string) const;
|
| 39 |
//Style Stuff
|
| 40 |
void setHistStyle(TH1*, const char*, const char*);
|
| 41 |
void setCanvasStyle( TCanvas& canv );
|
| 42 |
void setLegendStyle( TLegend& leg );
|
| 43 |
|
| 44 |
private:
|
| 45 |
//define input/output
|
| 46 |
std::string treename_, output_, name_data_, name_reference_;
|
| 47 |
std::vector<std::string> input_, reference_;
|
| 48 |
|
| 49 |
//define histogram design
|
| 50 |
std::vector<short int> cmpColor_, cmpStyle_, cmpWidth_;
|
| 51 |
|
| 52 |
//define canvas design
|
| 53 |
int gridX_, gridY_;
|
| 54 |
bool scale_;
|
| 55 |
|
| 56 |
//define legend design
|
| 57 |
double legXLeft_, legXRight_;
|
| 58 |
double legYLower_, legYUpper_;
|
| 59 |
|
| 60 |
//bins
|
| 61 |
std::vector<int> Bins_;
|
| 62 |
|
| 63 |
//hists
|
| 64 |
std::vector<std::string> HistNames_, Hist2DNames_;
|
| 65 |
|
| 66 |
//validation hist
|
| 67 |
std::vector<TH1F*> hists_, refs_;
|
| 68 |
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 |
};
|
| 76 |
|
| 77 |
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 |
void ValPlotMaker::setHistStyle( TH1* hist, const char* titleX, const char* titleY )
|
| 83 |
{
|
| 84 |
//hist->SetTitle( "" );
|
| 85 |
if ( hist->GetXaxis()->GetTitle()=="" )
|
| 86 |
hist->GetXaxis()->SetTitle( titleX );
|
| 87 |
hist->GetXaxis()->SetTitleSize ( 0.04 );
|
| 88 |
hist->GetXaxis()->SetTitleColor ( 1 );
|
| 89 |
hist->GetXaxis()->SetTitleOffset( 1.8 );
|
| 90 |
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 |
//hist->GetYaxis()->SetTitle( titleY );
|
| 96 |
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 |
//canv.SetTopMargin ( 0.05 );
|
| 110 |
}
|
| 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 |
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 |
}
|
| 143 |
|
| 144 |
|
| 145 |
template <class T> void ValPlotMaker::loadHistograms(vector<string> filenames, vector<string> *histnames, vector<T*> *hists)
|
| 146 |
{
|
| 147 |
std::vector<std::string>::iterator filename = filenames.begin();
|
| 148 |
TFile * file_ = new TFile( (*filename).c_str() );
|
| 149 |
T * dummy;
|
| 150 |
hists->clear();
|
| 151 |
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 |
string dummy_name = treename_+"/"+(*it);
|
| 175 |
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 |
}
|
| 181 |
}
|
| 182 |
}
|
| 183 |
|
| 184 |
if (hists->size() > 0) {
|
| 185 |
filename = filenames.begin();
|
| 186 |
cout << "...loaded " << hists->size()
|
| 187 |
<<" histograms from file(s) '" << (*filename) << "'." << endl;
|
| 188 |
}
|
| 189 |
}
|
| 190 |
|
| 191 |
void ValPlotMaker::loadHistograms()
|
| 192 |
{
|
| 193 |
loadHistograms<TH1F>(input_, &HistNames_, &hists_);
|
| 194 |
loadHistograms<TH1F>(reference_, &HistNames_, &refs_);
|
| 195 |
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 |
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 |
exit(3);
|
| 207 |
}
|
| 208 |
|
| 209 |
}
|
| 210 |
|
| 211 |
void ValPlotMaker::drawComparison()
|
| 212 |
{
|
| 213 |
TPostScript ps(output_.c_str(), 112 );
|
| 214 |
TCanvas *c1 = new TCanvas("c1", "fit histograms", 600, 600);
|
| 215 |
setCanvasStyle( *c1 );
|
| 216 |
c1->SetGridx( 1 );
|
| 217 |
c1->SetGridy( 1 );
|
| 218 |
|
| 219 |
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 |
++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 |
refer.SetTitle( "" );
|
| 244 |
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 |
refer.Scale( hist.Integral()/refer.Integral() );
|
| 254 |
}
|
| 255 |
// Statistical test of compatibility in shape between
|
| 256 |
// THIS histogram and h2, using Kolmogorov test.
|
| 257 |
ks = refer.KolmogorovTest( &hist, "");
|
| 258 |
// 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 |
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.5*refer.GetMaximum() : 1.5*hist.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.5*refer2D.GetMaximum() : 1.5*hist2D.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 |
|
| 387 |
|
| 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 |
ps.Close();
|
| 402 |
c1->Close();
|
| 403 |
delete c1;
|
| 404 |
}
|
| 405 |
|
| 406 |
|
| 407 |
int ValidationPlotMaker(int argc, char* argv[])
|
| 408 |
{
|
| 409 |
setNiceStyle();
|
| 410 |
gStyle->SetOptStat( 0 );
|
| 411 |
gStyle->SetOptFit ( 0 );
|
| 412 |
gStyle->SetStatColor(0);
|
| 413 |
gStyle->SetStatBorderSize(0);
|
| 414 |
gStyle->SetStatX(0.93);
|
| 415 |
gStyle->SetStatY(0.93);
|
| 416 |
gStyle->SetStatW(0.18);
|
| 417 |
gStyle->SetStatH(0.18);
|
| 418 |
|
| 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 |
ValPlotMaker * plots = new ValPlotMaker();
|
| 426 |
plots->readConfig( config_file.c_str() );
|
| 427 |
plots->loadHistograms();
|
| 428 |
plots->drawComparison();
|
| 429 |
delete plots;
|
| 430 |
|
| 431 |
return 0;
|
| 432 |
}
|
| 433 |
|
| 434 |
int main(int argc, char* argv[])
|
| 435 |
{
|
| 436 |
return ValidationPlotMaker(argc, argv);
|
| 437 |
}
|
