CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Comparator.cc
Go to the documentation of this file.
1 #include <TFile.h>
2 #include <TH1.h>
3 #include <TH2.h>
4 #include <TStyle.h>
5 #include <TPaveStats.h>
6 
7 
8 #include <sstream>
9 #include <cstdlib>
10 #include <cassert>
11 
15 
16 using namespace std;
17 
18 void Comparator::SetDirs( const char* file0,
19  const char* dir0,
20  const char* file1,
21  const char* dir1 ) {
22 
23  file0_ = new TFile( file0 );
24  if( file0_->IsZombie() ) exit(1);
25  dir0_ = file0_->GetDirectory( dir0 );
26  if(! dir0_ ) exit(1);
27 
28  file1_ = new TFile( file1 );
29  if( file1_->IsZombie() ) exit(1);
30  dir1_ = file1_->GetDirectory( dir1 );
31  if(! dir1_ ) exit(1);
32 }
33 
34 
36  Style* s1,
37  const char* leg0,
38  const char* leg1) {
39  s0_ = s0;
40  s1_ = s1;
41 
42  legend_.Clear();
43  legend_.AddEntry( s0_, leg0, "mlf");
44  legend_.AddEntry( s1_, leg1, "mlf");
45 }
46 
47 void Comparator::DrawSlice( const char* key,
48  int binxmin, int binxmax,
49  Mode mode ) {
50 
51  static int num = 0;
52 
53  ostringstream out0;
54  out0<<"h0_2d_"<<num;
55  ostringstream out1;
56  out1<<"h1_2d_"<<num;
57  num++;
58 
59  string name0 = out0.str();
60  string name1 = out1.str();
61 
62 
63  TH1* h0 = Histo( key, 0);
64  TH1* h1 = Histo( key, 1);
65 
66  TH2* h0_2d = dynamic_cast< TH2* >(h0);
67  TH2* h1_2d = dynamic_cast< TH2* >(h1);
68 
69  if(h0_2d->GetNbinsY() == 1 ||
70  h1_2d->GetNbinsY() == 1 ) {
71  cerr<<key<<" is not 2D"<<endl;
72  return;
73  }
74 
75  TH1::AddDirectory( false );
76 
77  TH1D* h0_slice = h0_2d->ProjectionY(name0.c_str(),
78  binxmin, binxmax, "");
79  TH1D* h1_slice = h1_2d->ProjectionY(name1.c_str(),
80  binxmin, binxmax, "");
81  TH1::AddDirectory( true );
82  Draw( h0_slice, h1_slice, mode);
83 }
84 
85 void Comparator::DrawMeanSlice(const char* key, const int rebinFactor, Mode mode)
86 {
87  TDirectory* dir = dir1_;
88  dir->cd();
89  TH2D *h2 = (TH2D*) dir->Get(key);
90  TH2Analyzer TH2Ana(h2,rebinFactor);
91  TH1D* ha=TH2Ana.Average();
92 
93  dir = dir0_;
94  dir->cd();
95  TH2D *h2b = (TH2D*) dir->Get(key);
96  TH2Analyzer TH2Anab(h2b,rebinFactor);
97  TH1D* hb=TH2Anab.Average();
98 
99  Draw(hb,ha,mode);
100 }
101 
102 void Comparator::DrawSigmaSlice(const char* key, const int rebinFactor, Mode mode)
103 {
104  TDirectory* dir = dir1_;
105  dir->cd();
106  TH2D *h2 = (TH2D*) dir->Get(key);
107  TH2Analyzer TH2Ana(h2,rebinFactor);
108  TH1D* ha=TH2Ana.RMS();
109 
110  dir = dir0_;
111  dir->cd();
112  TH2D *h2b = (TH2D*) dir->Get(key);
113  TH2Analyzer TH2Anab(h2b,rebinFactor);
114  TH1D* hb=TH2Anab.RMS();
115 
116  Draw(hb,ha,mode);
117 }
118 
119 void Comparator::DrawGaussSigmaSlice(const char* key, const int rebinFactor, Mode mode)
120 {
121  TDirectory* dir = dir1_;
122  dir->cd();
123  TH2D *h2 = (TH2D*) dir->Get(key);
124  TH2Analyzer TH2Ana(h2,rebinFactor);
125  TH1D* ha=TH2Ana.SigmaGauss();
126 
127  dir = dir0_;
128  dir->cd();
129  TH2D *h2b = (TH2D*) dir->Get(key);
130  TH2Analyzer TH2Anab(h2b,rebinFactor);
131  TH1D* hb=TH2Anab.SigmaGauss();
132 
133  Draw(hb,ha,mode);
134 }
135 
136 Double_t fitFunction_f(Double_t *x, Double_t *par)
137 {
138  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];
139  return value;
140 }
141 
142 void Comparator::DrawGaussSigmaSlice(const char* key, const int rebinFactor, const int binxmin,
143  const int binxmax, const bool cst_binning, Mode mode)
144 {
145  TDirectory* dir = dir1_;
146  dir->cd();
147  TH2D *h2 = (TH2D*) dir->Get(key);
148  TH2Analyzer TH2Ana(h2, binxmin, binxmax, rebinFactor, cst_binning);
149  TH1D* hrms=TH2Ana.RMS();
150  TF1 *fitfcndgssrms3 = new TF1("fitfcndgssrms3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
151  fitfcndgssrms3->SetNpx(500);
152  fitfcndgssrms3->SetLineWidth(3);
153  fitfcndgssrms3->SetLineStyle(2);
154  fitfcndgssrms3->SetLineColor(4);
155  hrms->Fit("fitfcndgssrms3","0R");
156 
157  TH1D* ha=TH2Ana.SigmaGauss();
158  TF1 *fitfcndgsse3 = new TF1("fitfcndgsse3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
159  fitfcndgsse3->SetNpx(500);
160  fitfcndgsse3->SetLineWidth(3);
161  fitfcndgsse3->SetLineStyle(1);
162  fitfcndgsse3->SetLineColor(4);
163  ha->Fit("fitfcndgsse3","0R");
164 
165  dir = dir0_;
166  dir->cd();
167  TH2D *h2b = (TH2D*) dir->Get(key);
168  TH2Analyzer TH2Anab(h2b, binxmin, binxmax, rebinFactor, cst_binning);
169  TH1D* hrmsb=TH2Anab.RMS();
170  TF1 *fitfcndgssrmsb3 = new TF1("fitfcndgssrmsb3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
171  fitfcndgssrmsb3->SetNpx(500);
172  fitfcndgssrmsb3->SetLineWidth(3);
173  fitfcndgssrmsb3->SetLineStyle(2);
174  fitfcndgssrmsb3->SetLineColor(2);
175  hrmsb->Fit("fitfcndgssrmsb3","0R");
176 
177  TH1D* hb=TH2Anab.SigmaGauss();
178  TF1 *fitfcndgsseb3 = new TF1("fitfcndgsseb3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
179  fitfcndgsseb3->SetNpx(500);
180  fitfcndgsseb3->SetLineWidth(3);
181  fitfcndgsseb3->SetLineStyle(1);
182  fitfcndgsseb3->SetLineColor(2);
183  hb->Fit("fitfcndgsseb3","0R");
184 
185  Draw(hb,ha,mode);
186  //Draw(hrms,ha,mode);
187  //Draw(ha,ha,mode);
188  fitfcndgssrms3->Draw("same");
189  fitfcndgsse3->Draw("same");
190  fitfcndgssrmsb3->Draw("same");
191  fitfcndgsseb3->Draw("same");
192 }
193 
194 void Comparator::DrawGaussSigmaOverMeanXSlice(const char* key, const int rebinFactor, const int binxmin,
195  const int binxmax, const bool cst_binning, Mode mode)
196 {
197 
198  TDirectory* dir = dir1_;
199  dir->cd();
200  TH2D *h2 = (TH2D*) dir->Get(key);
201  TH2Analyzer TH2Ana(h2, binxmin, binxmax, rebinFactor, cst_binning);
202  TH1D* hrms=TH2Ana.RMS();
203 
204  TH1D* meanXslice = TH2Ana.MeanX();
205  //for( int i=1; i<=meanXslice->GetNbinsX(); ++i) {
206  // std::cout << "meanXslice->GetBinContent(" << i << ") = "
207  // << meanXslice->GetBinContent(i) << std::endl;
208  // std::cout << "meanXslice->GetBinError(" << i << ") = "
209  // << meanXslice->GetBinError(i) << std::endl;
210  //}
211  //Draw(meanXslice,meanXslice,mode);
212  hrms->Divide(meanXslice);
213 
214  TF1 *fitXfcndgssrms3 = new TF1("fitXfcndgssrms3",fitFunction_f,hrms->GetXaxis()->GetBinLowEdge(1),hrms->GetXaxis()->GetBinUpEdge(hrms->GetNbinsX()),4);
215  fitXfcndgssrms3->SetNpx(500);
216  fitXfcndgssrms3->SetLineWidth(3);
217  fitXfcndgssrms3->SetLineStyle(2);
218  fitXfcndgssrms3->SetLineColor(4);
219  hrms->Fit("fitXfcndgssrms3","0R");
220 
221  TH1D* ha=TH2Ana.SigmaGauss();
222  ha->Divide(meanXslice);
223  TF1 *fitXfcndgsse3 = new TF1("fitXfcndgsse3",fitFunction_f,ha->GetXaxis()->GetBinLowEdge(1),ha->GetXaxis()->GetBinUpEdge(ha->GetNbinsX()),4);
224  fitXfcndgsse3->SetNpx(500);
225  fitXfcndgsse3->SetLineWidth(3);
226  fitXfcndgsse3->SetLineStyle(1);
227  fitXfcndgsse3->SetLineColor(4);
228  ha->Fit("fitXfcndgsse3","0R");
229 
230  dir = dir0_;
231  dir->cd();
232  TH2D *h2b = (TH2D*) dir->Get(key);
233  TH2Analyzer TH2Anab(h2b, binxmin, binxmax, rebinFactor, cst_binning);
234  TH1D* hrmsb=TH2Anab.RMS();
235  hrmsb->Divide(meanXslice);
236  TF1 *fitXfcndgssrmsb3 = new TF1("fitXfcndgssrmsb3",fitFunction_f,hrmsb->GetXaxis()->GetBinLowEdge(1),hrmsb->GetXaxis()->GetBinUpEdge(hrmsb->GetNbinsX()),4);
237  fitXfcndgssrmsb3->SetNpx(500);
238  fitXfcndgssrmsb3->SetLineWidth(3);
239  fitXfcndgssrmsb3->SetLineStyle(2);
240  fitXfcndgssrmsb3->SetLineColor(2);
241  hrmsb->Fit("fitXfcndgssrmsb3","0R");
242 
243  TH1D* hb=TH2Anab.SigmaGauss();
244  hb->Divide(meanXslice);
245  TF1 *fitXfcndgsseb3 = new TF1("fitXfcndgsseb3",fitFunction_f,hb->GetXaxis()->GetBinLowEdge(1),hb->GetXaxis()->GetBinUpEdge(hb->GetNbinsX()),4);
246  fitXfcndgsseb3->SetNpx(500);
247  fitXfcndgsseb3->SetLineWidth(3);
248  fitXfcndgsseb3->SetLineStyle(1);
249  fitXfcndgsseb3->SetLineColor(2);
250  hb->Fit("fitXfcndgsseb3","0R");
251 
252  Draw(hb,ha,mode);
253  //Draw(hrms,ha,mode);
254  //Draw(ha,ha,mode);
255  fitXfcndgssrms3->Draw("same");
256  fitXfcndgsse3->Draw("same");
257  fitXfcndgssrmsb3->Draw("same");
258  fitXfcndgsseb3->Draw("same");
259 }
260 
261 void Comparator::DrawGaussSigmaOverMeanSlice(const char* key, const char* key2,
262  const int rebinFactor, Mode mode)
263 {
264  TDirectory* dir = dir1_;
265  dir->cd();
266 
267  TH2D *h2_b = (TH2D*) dir->Get(key2);
268  TH2Analyzer TH2Ana_b(h2_b, rebinFactor);
269  TH1D* meanslice=TH2Ana_b.Average();
270 
271  TH2D *h2 = (TH2D*) dir->Get(key);
272  TH2Analyzer TH2Ana(h2, rebinFactor);
273 
274  TH1D* ha=TH2Ana.SigmaGauss();
275  ha->Divide(meanslice);
276 
277  dir = dir0_;
278  dir->cd();
279  TH2D *h2b = (TH2D*) dir->Get(key);
280  TH2Analyzer TH2Anab(h2b, rebinFactor);
281 
282  TH2D *h2b_b = (TH2D*) dir->Get(key2);
283  TH2Analyzer TH2Anab_b(h2b_b, rebinFactor);
284  TH1D* meansliceb=TH2Anab_b.Average();
285 
286  TH1D* hb=TH2Anab.SigmaGauss();
287  hb->Divide(meansliceb);
288 
289  Draw(hb,ha,mode);
290  //Draw(meansliceb,meanslice,mode);
291 }
292 
293 void Comparator::Draw( const char* key, Mode mode) {
294 
295  TH1::AddDirectory( false );
296  TH1* h0 = Histo( key, 0);
297  TH1* h1 = (TH1*) Histo( key, 1)->Clone("h1");
298 
299  TH1::AddDirectory( true );
300  Draw( h0, h1, mode);
301 }
302 
303 
304 void Comparator::Draw( const char* key0, const char* key1, Mode mode) {
305  TH1* h0=0;
306  TH1* h1=0;
307  if(mode!=EFF)
308  {
309  h0 = Histo( key0, 0);
310  h1 = Histo( key1, 1);
311  }
312  else
313  {
314  h0 = Histo( key0, 0);
315  TH1 * h0b = Histo( key1, 0);
316  h1 = Histo( key0, 1);
317  TH1 * h1b = Histo( key1, 1);
318  if(rebin_>1) {
319  h0->Rebin( rebin_);
320  h1->Rebin( rebin_);
321  h0b->Rebin( rebin_);
322  h1b->Rebin( rebin_);
323  }
324  if(resetAxis_) {
325  h0->GetXaxis()->SetRangeUser( xMin_, xMax_);
326  h1->GetXaxis()->SetRangeUser( xMin_, xMax_);
327  h0b->GetXaxis()->SetRangeUser( xMin_, xMax_);
328  h1b->GetXaxis()->SetRangeUser( xMin_, xMax_);
329  }
330 
331  h0b->Sumw2();
332  h0->Sumw2();
333  h0->Divide(h0,h0b,1.,1.,"B");
334  h1b->Sumw2();
335  h1->Sumw2();
336  h1->Divide(h1,h1b,1.,1.,"B");
337  }
338  Draw( h0, h1, mode);
339 }
340 
341 TH1* Comparator::Histo( const char* key, unsigned dirIndex) {
342  if(dirIndex>1U) { // dirIndex >= 0, since dirIndex is unsigned
343  cerr<<"bad dir index: "<<dirIndex<<endl;
344  return 0;
345  }
346  TDirectory* dir = 0;
347  if(dirIndex == 0) dir = dir0_;
348  if(dirIndex == 1) dir = dir1_;
349  assert( dir );
350 
351  dir->cd();
352 
353  TH1* h = (TH1*) dir->Get(key);
354  if(!h)
355  cerr<<"no key "<<key<<" in directory "<<dir->GetName()<<endl;
356  return h;
357 }
358 
359 
360 void Comparator::Draw( TH1* h0, TH1* h1, Mode mode ) {
361  if( !(h0 && h1) ) {
362  cerr<<"invalid histo"<<endl;
363  return;
364  }
365 
366  TH1::AddDirectory( false );
367  h0_ = (TH1*) h0->Clone( "h0_");
368  h1_ = (TH1*) h1->Clone( "h1_");
369  TH1::AddDirectory( true );
370 
371  // unsetting the title, since the title of projections
372  // is still the title of the 2d histo
373  // and this is better anyway
374  h0_->SetTitle("");
375  h1_->SetTitle("");
376 
377  //h0_->SetStats(1);
378  //h1_->SetStats(1);
379 
380  if(mode!=EFF)
381  {
382  if(rebin_>1) {
383  h0_->Rebin( rebin_);
384  h1_->Rebin( rebin_);
385  }
386  if(resetAxis_) {
387  h0_->GetXaxis()->SetRangeUser( xMin_, xMax_);
388  h1_->GetXaxis()->SetRangeUser( xMin_, xMax_);
389  }
390  }
391 
392  if (mode!=GRAPH)
393  {
394 
395  TPaveStats *ptstats = new TPaveStats(0.7385057,0.720339,
396  0.9396552,0.8792373,"brNDC");
397  ptstats->SetName("stats");
398  ptstats->SetBorderSize(1);
399  ptstats->SetLineColor(2);
400  ptstats->SetFillColor(10);
401  ptstats->SetTextAlign(12);
402  ptstats->SetTextColor(2);
403  ptstats->SetOptStat(1111);
404  ptstats->SetOptFit(0);
405  ptstats->Draw();
406  h0_->GetListOfFunctions()->Add(ptstats);
407  ptstats->SetParent(h0_->GetListOfFunctions());
408 
409  //std::cout << "FL: h0_->GetMean() = " << h0_->GetMean() << std::endl;
410  //std::cout << "FL: h0_->GetRMS() = " << h0_->GetRMS() << std::endl;
411  //std::cout << "FL: h1_->GetMean() = " << h1_->GetMean() << std::endl;
412  //std::cout << "FL: h1_->GetRMS() = " << h1_->GetRMS() << std::endl;
413  //std::cout << "FL: test2" << std::endl;
414  TPaveStats *ptstats2 = new TPaveStats(0.7399425,0.529661,
415  0.941092,0.6885593,"brNDC");
416  ptstats2->SetName("stats");
417  ptstats2->SetBorderSize(1);
418  ptstats2->SetLineColor(4);
419  ptstats2->SetFillColor(10);
420  ptstats2->SetTextAlign(12);
421  ptstats2->SetTextColor(4);
422  TText *text = ptstats2->AddText("h1_");
423  text->SetTextSize(0.03654661);
424 
425  std::ostringstream oss3;
426  oss3 << h1_->GetEntries();
427  const std::string txt_entries="Entries = "+oss3.str();
428  text = ptstats2->AddText(txt_entries.c_str());
429  std::ostringstream oss;
430  oss << h1_->GetMean();
431  const std::string txt_mean="Mean = "+oss.str();
432  text = ptstats2->AddText(txt_mean.c_str());
433  std::ostringstream oss2;
434  oss2 << h1_->GetRMS();
435  const std::string txt_rms="RMS = "+oss2.str();
436  text = ptstats2->AddText(txt_rms.c_str());
437  ptstats2->SetOptStat(1111);
438  ptstats2->SetOptFit(0);
439  ptstats2->Draw();
440  h1_->GetListOfFunctions()->Add(ptstats2);
441  ptstats2->SetParent(h1_->GetListOfFunctions());
442  }
443  else
444  {
445  TPaveStats *ptstats = new TPaveStats(0.0,0.0,
446  0.0,0.0,"brNDC");
447  ptstats->Draw();
448  h0_->GetListOfFunctions()->Add(ptstats);
449  ptstats->SetParent(h0_->GetListOfFunctions());
450  }
451 
452  float min=-999.;
453  float max=+999.;
454  switch(mode) {
455  case SCALE:
456  h1_->Scale( h0_->GetEntries()/h1_->GetEntries() );
457  case NORMAL:
458  if(s0_)
459  Styles::FormatHisto( h0_ , s0_);
460  if(s1_)
461  Styles::FormatHisto( h1_ , s1_);
462 
463  if( h1_->GetMaximum()>h0_->GetMaximum()) {
464  h0_->SetMaximum( h1_->GetMaximum()*1.15 );
465  }
466 
467  h0_->Draw();
468  h1_->Draw("same");
469 
470  break;
471  case EFF:
472  if(s0_)
473  Styles::FormatHisto( h0_ , s0_);
474  if(s1_)
475  Styles::FormatHisto( h1_ , s1_);
476 
477  // rather arbitrary but useful
478  max= h0_->GetMaximum();
479  if ( h1_->GetMaximum() > max )
480  max = h1_->GetMaximum();
481  if(max > 0.8) max = 1;
482 
483  max *= 1.1 ;
484 
485  min= h0_->GetMinimum();
486  if ( h1_->GetMinimum() < min )
487  min = h1_->GetMinimum();
488  if(min > 0.2) min = 0.;
489 
490  min *= 0.8 ;
491  h0_->SetMaximum( max );
492  h0_->SetMinimum( min );
493 
494  h0_->Draw("E");
495  h1_->Draw("Esame");
496  break;
497  case GRAPH:
498  if(s0_)
499  Styles::FormatHisto( h0_ , s0_);
500  if(s1_)
501  Styles::FormatHisto( h1_ , s1_);
502 
503  if( h1_->GetMaximum()>h0_->GetMaximum()) {
504  h0_->SetMaximum( h1_->GetMaximum()*1.15 );
505  }
506  if( h1_->GetMinimum()<h0_->GetMinimum()) {
507  h0_->SetMinimum( h1_->GetMinimum()*1.15 );
508  }
509  h0_->SetMarkerStyle(21);
510  h0_->SetMarkerColor(2);
511  h1_->SetMarkerStyle(21);
512  h1_->SetMarkerColor(4);
513 
514  h0_->Draw("E1");
515  h1_->Draw("E1same");
516  break;
517  case RATIO:
518  h0_->Sumw2();
519  h1_->Sumw2();
520  h0_->Divide( h1_ );
521  if(s0_)
522  Styles::FormatHisto( h0_ , s0_);
523  h0_->Draw();
524  default:
525  break;
526  }
527 }
528 
TH1D * SigmaGauss()
Definition: TH2Analyzer.h:56
void SetStyles(Style *s0, Style *s1, const char *leg0, const char *leg1)
Definition: Comparator.cc:35
static void FormatHisto(TH1 *h, const Style *s)
Definition: NicePlot.cc:58
void DrawSigmaSlice(const char *key, const int rebinFactor, Mode mode)
Definition: Comparator.cc:102
void Draw(const char *key, Mode mode)
Definition: Comparator.cc:293
void DrawGaussSigmaOverMeanXSlice(const char *key, const int rebinFactor, const int binxmin, const int binxmax, const bool cst_binning, Mode mode)
Definition: Comparator.cc:194
void DrawMeanSlice(const char *key, const int rebinFactor, Mode mode)
Definition: Comparator.cc:85
void SetDirs(const char *file0, const char *dir0, const char *file1, const char *dir1)
set the 2 files, and the directory within each file, in which the histograms will be compared ...
Definition: Comparator.cc:18
const T & max(const T &a, const T &b)
T sqrt(T t)
Definition: SSEVec.h:48
#define RATIO(A, B)
tuple text
Definition: runonSM.py:42
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
TH1 * Histo(const char *key, unsigned dirIndex)
Definition: Comparator.cc:341
Double_t fitFunction_f(Double_t *x, Double_t *par)
Definition: Comparator.cc:136
TH1D * MeanX()
Definition: TH2Analyzer.h:57
void DrawSlice(const char *key, int binxmin, int binxmax, Mode mode)
Definition: Comparator.cc:47
TH1D * RMS()
Definition: TH2Analyzer.h:55
void DrawGaussSigmaOverMeanSlice(const char *key, const char *key2, const int rebinFactor, Mode mode)
Definition: Comparator.cc:261
list key
Definition: combine.py:13
void DrawGaussSigmaSlice(const char *key, const int rebinFactor, Mode mode)
Definition: Comparator.cc:119
dbl *** dir
Definition: mlp_gen.cc:35
Definition: DDAxes.h:10
TH1D * Average()
Definition: TH2Analyzer.h:54
Definition: NicePlot.h:7