| 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::vector<std::string> filename, std::vector<string> *histnames, std::vector<TH1F*> *hist);
|
| 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(TH1F*, 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 |
|
| 55 |
//define legend design
|
| 56 |
double legXLeft_, legXRight_;
|
| 57 |
double legYLower_, legYUpper_;
|
| 58 |
|
| 59 |
//bins
|
| 60 |
std::vector<int> Bins_;
|
| 61 |
|
| 62 |
//hists
|
| 63 |
std::vector<std::string> HistNames_, Hist2DNames_;
|
| 64 |
|
| 65 |
//validation hist
|
| 66 |
std::vector<TH1F*> hists_, refs_;
|
| 67 |
|
| 68 |
//background color of hists, if KS test fails
|
| 69 |
double ks_treshold;
|
| 70 |
int ks_color;
|
| 71 |
bool UseTitle_;
|
| 72 |
|
| 73 |
};
|
| 74 |
|
| 75 |
bool ValPlotMaker::contains(const string the_string, const string the_sub_string) const
|
| 76 |
{
|
| 77 |
return the_string.find(the_sub_string,0) != string::npos;
|
| 78 |
}
|
| 79 |
|
| 80 |
void ValPlotMaker::setHistStyle( TH1F* hist, const char* titleX, const char* titleY )
|
| 81 |
{
|
| 82 |
//hist->SetTitle( "" );
|
| 83 |
if ( hist->GetXaxis()->GetTitle()=="" )
|
| 84 |
hist->GetXaxis()->SetTitle( titleX );
|
| 85 |
hist->GetXaxis()->SetTitleSize ( 0.04 );
|
| 86 |
hist->GetXaxis()->SetTitleColor ( 1 );
|
| 87 |
hist->GetXaxis()->SetTitleOffset( 1.8 );
|
| 88 |
hist->GetXaxis()->SetTitleFont ( 62 );
|
| 89 |
hist->GetXaxis()->SetLabelSize ( 0.05 );
|
| 90 |
hist->GetXaxis()->SetLabelFont ( 62 );
|
| 91 |
hist->GetXaxis()->CenterTitle ( );
|
| 92 |
hist->GetXaxis()->SetNdivisions ( 505 );
|
| 93 |
//hist->GetYaxis()->SetTitle( titleY );
|
| 94 |
hist->GetYaxis()->SetTitleSize ( 0.07 );
|
| 95 |
hist->GetYaxis()->SetTitleColor ( 1 );
|
| 96 |
hist->GetYaxis()->SetTitleOffset( 1.3 );
|
| 97 |
hist->GetYaxis()->SetTitleFont ( 62 );
|
| 98 |
hist->GetYaxis()->SetLabelSize ( 0.05 );
|
| 99 |
hist->GetYaxis()->SetLabelFont ( 62 );
|
| 100 |
}
|
| 101 |
void ValPlotMaker::setCanvasStyle( TCanvas& canv )
|
| 102 |
{
|
| 103 |
canv.SetFillStyle ( 4000 );
|
| 104 |
canv.SetLeftMargin ( 0.20 );
|
| 105 |
canv.SetRightMargin ( 0.05 );
|
| 106 |
canv.SetBottomMargin( 0.15 );
|
| 107 |
//canv.SetTopMargin ( 0.05 );
|
| 108 |
}
|
| 109 |
void ValPlotMaker::setLegendStyle( TLegend& leg )
|
| 110 |
{
|
| 111 |
leg.SetFillStyle ( 0 );
|
| 112 |
leg.SetBorderSize( 0 );
|
| 113 |
leg.SetFillColor ( 0 );
|
| 114 |
}
|
| 115 |
|
| 116 |
void ValPlotMaker::readConfig( std::string name )
|
| 117 |
{
|
| 118 |
ConfigFile cfg( name );
|
| 119 |
input_ = bag_of_string(cfg.read<std::string>( "input_root" ));
|
| 120 |
reference_ = bag_of_string(cfg.read<std::string>( "input_reference" ));
|
| 121 |
name_data_ = cfg.read<std::string>( "name data", "data" );
|
| 122 |
name_reference_ = cfg.read<std::string>( "name reference", "reference" );
|
| 123 |
treename_ = cfg.read<std::string>( "tree_name" );
|
| 124 |
output_ = cfg.read<std::string>( "output" );
|
| 125 |
HistNames_ = bag_of_string(cfg.read<std::string>( "hists"));
|
| 126 |
ks_treshold = cfg.read<double> ( "KS treshold", 0.05 );
|
| 127 |
ks_color = cfg.read<int> ( "KS color", 5 );
|
| 128 |
gridX_ = cfg.read<int> ( "grid_x" );
|
| 129 |
gridY_ = cfg.read<int> ( "grid_y" );
|
| 130 |
legXLeft_ = cfg.read<double> ( "leg_x_left" );
|
| 131 |
legXRight_ = cfg.read<double> ( "leg_x_right" );
|
| 132 |
legYLower_ = cfg.read<double> ( "leg_y_lower" );
|
| 133 |
legYUpper_ = cfg.read<double> ( "leg_y_upper" );
|
| 134 |
cmpColor_ = bag_of<short int>( cfg.read<std::string>( "color" ));
|
| 135 |
cmpStyle_ = bag_of<short int>( cfg.read<std::string>( "style" ));
|
| 136 |
cmpWidth_ = bag_of<short int>( cfg.read<std::string>( "width" ));
|
| 137 |
UseTitle_ = cfg.read<int> ( "show title", false );
|
| 138 |
}
|
| 139 |
|
| 140 |
|
| 141 |
void ValPlotMaker::loadHistograms(vector<string> filenames, vector<string> *histnames, vector<TH1F*> *hists)
|
| 142 |
{
|
| 143 |
std::vector<std::string>::iterator filename = filenames.begin();
|
| 144 |
TFile * file_ = new TFile( (*filename).c_str() );
|
| 145 |
TH1F * dummy;
|
| 146 |
hists->clear();
|
| 147 |
for (std::vector<std::string>::iterator it=histnames->begin();
|
| 148 |
it!=histnames->end(); ) {
|
| 149 |
string dummy_name = treename_+"/"+(*it);
|
| 150 |
dummy = (TH1F*)file_->Get( dummy_name.c_str() );
|
| 151 |
if( !dummy ){
|
| 152 |
cerr << "WARNING:"
|
| 153 |
<< " Didn't find indicated hist" << " [" << (*it) << "]"
|
| 154 |
<< " in tree" << " [" << treename_ << "]"
|
| 155 |
<< " of file" << " [" << (*filename) << "]"
|
| 156 |
<< endl;
|
| 157 |
it = histnames->erase( it );
|
| 158 |
} else {
|
| 159 |
hists->push_back( dummy );
|
| 160 |
++it;
|
| 161 |
}
|
| 162 |
}
|
| 163 |
++filename;
|
| 164 |
for (; filename!=filenames.end(); ++filename){
|
| 165 |
TFile * file_ = new TFile( (*filename).c_str() );
|
| 166 |
int i=0;
|
| 167 |
std::vector<TH1F*>::iterator hist=hists->begin();
|
| 168 |
for (std::vector<std::string>::iterator it=histnames->begin();
|
| 169 |
it!=histnames->end()&&hist!=hists->end(); ++hist,++it,++i) {
|
| 170 |
string dummy_name = treename_+"/"+(*it);
|
| 171 |
dummy = (TH1F*)file_->Get( dummy_name.c_str() );
|
| 172 |
if( dummy ){
|
| 173 |
dummy->Add( (*hist) );
|
| 174 |
hist = hists->erase(hist);
|
| 175 |
hists->insert( hist, 1, dummy );
|
| 176 |
}
|
| 177 |
}
|
| 178 |
}
|
| 179 |
|
| 180 |
if (hists->size() > 0) {
|
| 181 |
filename = filenames.begin();
|
| 182 |
cout << "...loaded " << hists->size()
|
| 183 |
<<" histograms from file(s) '" << (*filename) << "'." << endl;
|
| 184 |
}
|
| 185 |
}
|
| 186 |
|
| 187 |
void ValPlotMaker::loadHistograms()
|
| 188 |
{
|
| 189 |
loadHistograms(input_, &HistNames_, &hists_);
|
| 190 |
loadHistograms(reference_, &HistNames_, &refs_);
|
| 191 |
if (hists_.size()!=refs_.size()) {
|
| 192 |
cerr << "ERROR: At least one histogram is only present in one input file!"
|
| 193 |
<< "The number of data and reference plots does not match." << endl;
|
| 194 |
exit(2);
|
| 195 |
}
|
| 196 |
|
| 197 |
}
|
| 198 |
|
| 199 |
void ValPlotMaker::drawComparison()
|
| 200 |
{
|
| 201 |
TPostScript ps(output_.c_str(), 112 );
|
| 202 |
TCanvas *c1 = new TCanvas("c1", "fit histograms", 600, 600);
|
| 203 |
setCanvasStyle( *c1 );
|
| 204 |
c1->SetGridx( 1 );
|
| 205 |
c1->SetGridy( 1 );
|
| 206 |
|
| 207 |
char * dummy_text = new char[100];
|
| 208 |
unsigned i=0;
|
| 209 |
double ks;
|
| 210 |
std::vector<std::string>::iterator it=HistNames_.begin();
|
| 211 |
std::vector<TH1F*>::iterator hst=hists_.begin();
|
| 212 |
std::vector<TH1F*>::iterator ref=refs_.begin();
|
| 213 |
|
| 214 |
for (;ref!=refs_.end()&&hst!=hists_.end()&&it!=HistNames_.end();
|
| 215 |
++ref,++hst,++it,++i){
|
| 216 |
//cout << i << "th plot: " << (*it) << endl;
|
| 217 |
|
| 218 |
TH1F& hist = *((TH1F*)(*hst));
|
| 219 |
TH1F& refer = *((TH1F*)(*ref));
|
| 220 |
//check if hist is filled
|
| 221 |
if (refer.Integral()==0.0) {
|
| 222 |
cerr << "Histogram '"<< (*it) <<"' is empty!" <<endl;
|
| 223 |
continue;
|
| 224 |
}
|
| 225 |
|
| 226 |
setHistStyle( &refer, it->c_str(), "events" );
|
| 227 |
if (UseTitle_){
|
| 228 |
//refer.SetTitle( "" );
|
| 229 |
}
|
| 230 |
else
|
| 231 |
refer.SetTitle( "" );
|
| 232 |
refer.SetLineWidth( 8 );
|
| 233 |
refer.SetLineColor( 2 );
|
| 234 |
|
| 235 |
cout << i << "th plot: " << (*it)
|
| 236 |
<< " has " << refer.Integral() << " 'reference' events and "
|
| 237 |
<< hist.Integral() << " 'data' events" << endl;
|
| 238 |
|
| 239 |
refer.Sumw2();
|
| 240 |
if ( !contains(refer.GetYaxis()->GetTitle(), "eff") ) {
|
| 241 |
refer.Scale( hist.Integral()/refer.Integral() );
|
| 242 |
} else {
|
| 243 |
if ( contains(refer.GetXaxis()->GetTitle(), "Pt") ) {
|
| 244 |
hists_.at(6)->Sumw2();
|
| 245 |
hists_.at(3)->Sumw2();
|
| 246 |
hist.Divide(hists_.at(6),hists_.at(3),1.,1.,"B");
|
| 247 |
refs_.at(6)->Sumw2();
|
| 248 |
refs_.at(3)->Sumw2();
|
| 249 |
refer.Divide(refs_.at(6),refs_.at(3),1.,1.,"B");
|
| 250 |
}
|
| 251 |
if ( contains(refer.GetXaxis()->GetTitle(), "eta") ) {
|
| 252 |
hists_.at(7)->Sumw2();
|
| 253 |
hists_.at(4)->Sumw2();
|
| 254 |
hist.Divide(hists_.at(7),hists_.at(4),1.,1.,"B");
|
| 255 |
refs_.at(7)->Sumw2();
|
| 256 |
refs_.at(4)->Sumw2();
|
| 257 |
refer.Divide(refs_.at(7),refs_.at(4),1.,1.,"B");
|
| 258 |
}
|
| 259 |
if ( contains(refer.GetXaxis()->GetTitle(), "phi") ) {
|
| 260 |
hists_.at(8)->Sumw2();
|
| 261 |
hists_.at(5)->Sumw2();
|
| 262 |
hist.Divide(hists_.at(8),hists_.at(5),1.,1.,"B");
|
| 263 |
refs_.at(8)->Sumw2();
|
| 264 |
refs_.at(5)->Sumw2();
|
| 265 |
refer.Divide(refs_.at(8),refs_.at(5),1.,1.,"B");
|
| 266 |
}
|
| 267 |
}
|
| 268 |
// Statistical test of compatibility in shape between
|
| 269 |
// THIS histogram and h2, using Kolmogorov test.
|
| 270 |
ks = refer.KolmogorovTest( &hist, "");
|
| 271 |
// Default: Ignore under- and overflow bins in comparison
|
| 272 |
// option is a character string to specify options
|
| 273 |
// "U" include Underflows in test (also for 2-dim)
|
| 274 |
// "O" include Overflows (also valid for 2-dim)
|
| 275 |
// "N" include comparison of normalizations
|
| 276 |
// "D" Put out a line of "Debug" printout
|
| 277 |
// "M" Return the Maximum Kolmogorov distance instead of prob
|
| 278 |
// "X" Run the pseudo experiments post-processor with the following procedure:
|
| 279 |
// make pseudoexperiments based on random values from the parent
|
| 280 |
// distribution and compare the KS distance of the pseudoexperiment
|
| 281 |
// to the parent distribution. Bin the KS distances in a histogram,
|
| 282 |
// and then take the integral of all the KS values above the value
|
| 283 |
// obtained from the original data to Monte Carlo distribution.
|
| 284 |
// The number of pseudo-experiments nEXPT is currently fixed at 1000.
|
| 285 |
// The function returns the integral.
|
| 286 |
// (thanks to Ben Kilminster to submit this procedure). Note that
|
| 287 |
// this option "X" is much slower.
|
| 288 |
if (ks<ks_treshold){
|
| 289 |
c1->SetFillColor(ks_color);
|
| 290 |
cout << "Deviation in plot "<< i <<" '" << (*it) << "': " << ks
|
| 291 |
<< " (Kolmogorov-Smirnov Test)" << endl;
|
| 292 |
}
|
| 293 |
else {
|
| 294 |
c1->SetFillColor(0);
|
| 295 |
}
|
| 296 |
|
| 297 |
double maximum = (refer.GetMaximum()>hist.GetMaximum() ?
|
| 298 |
1.5*refer.GetMaximum() : 1.5*hist.GetMaximum());
|
| 299 |
refer.SetMaximum(maximum);
|
| 300 |
hist.SetMarkerStyle(8);
|
| 301 |
hist.SetLineWidth(3);
|
| 302 |
refer.Draw("pe");
|
| 303 |
hist.Draw("pe,same");
|
| 304 |
TLegend* leg = new TLegend(0.5,0.67,0.95,0.87);
|
| 305 |
setLegendStyle( *leg );
|
| 306 |
sprintf(dummy_text,"KS = %2.2f",ks);
|
| 307 |
leg->SetHeader(dummy_text);
|
| 308 |
leg->AddEntry( &hist, name_data_.c_str(), "PE" );
|
| 309 |
leg->AddEntry( &refer,name_reference_.c_str(), "L" );
|
| 310 |
leg->Draw("same");
|
| 311 |
|
| 312 |
c1->Draw();
|
| 313 |
ps.NewPage();
|
| 314 |
}
|
| 315 |
|
| 316 |
ps.Close();
|
| 317 |
c1->Close();
|
| 318 |
delete c1;
|
| 319 |
}
|
| 320 |
|
| 321 |
|
| 322 |
int ValidationPlotMaker_Eff(int argc, char* argv[])
|
| 323 |
{
|
| 324 |
setNiceStyle();
|
| 325 |
gStyle->SetOptStat( 0 );
|
| 326 |
gStyle->SetOptFit ( 0 );
|
| 327 |
gStyle->SetStatColor(0);
|
| 328 |
gStyle->SetStatBorderSize(0);
|
| 329 |
gStyle->SetStatX(0.93);
|
| 330 |
gStyle->SetStatY(0.93);
|
| 331 |
gStyle->SetStatW(0.18);
|
| 332 |
gStyle->SetStatH(0.18);
|
| 333 |
|
| 334 |
//If no argument is given, use default config file name
|
| 335 |
string config_file = "bin/ValidationPlot_Eff.cfg";
|
| 336 |
if( argc>=2 )
|
| 337 |
config_file = argv[1];
|
| 338 |
|
| 339 |
//Run the Validation Plot Maker
|
| 340 |
ValPlotMaker * plots = new ValPlotMaker();
|
| 341 |
plots->readConfig( config_file.c_str() );
|
| 342 |
plots->loadHistograms();
|
| 343 |
plots->drawComparison();
|
| 344 |
delete plots;
|
| 345 |
|
| 346 |
return 0;
|
| 347 |
}
|
| 348 |
|
| 349 |
int main(int argc, char* argv[])
|
| 350 |
{
|
| 351 |
return ValidationPlotMaker_Eff(argc, argv);
|
| 352 |
}
|
