CMS 3D CMS Logo

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