Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
errorOfHistograms.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file errorOfHistograms.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 /*
10  * errorOfHistograms.C
11  *
12  * Created on: 24 Jan 2017
13  * Author: jhrdinka
14  */
15 
16 #include <iostream>
17 #include "TFile.h"
18 #include "TH1F.h"
19 #include "TROOT.h"
20 
21 // This root script compares either two 1D or two 2D histograms with each other,
22 // by plotting them in different colors into the same canvas and displaying the
23 // relative error beneath. The histograms to be compared can either be directly
24 // in the file or in a directory.
25 
26 void
28  std::string hist1Name,
29  std::string legendName1,
30  std::string inFile2,
31  std::string hist2Name,
32  std::string legendName2,
33  std::string dirName1 = "",
34  std::string dirName2 = "")
35 {
36  std::cout << "Opening file: " << inFile1 << std::endl;
37  TFile inputFile1(inFile1.c_str());
38  std::cout << "Opening file: " << inFile2 << std::endl;
39  TFile inputFile2(inFile2.c_str());
40  std::cout << "Comparing Histograms: " << hist1Name << " & " << hist2Name
41  << std::endl;
42 
43  TH1F* h1 = nullptr;
44  TH1F* h2 = nullptr;
45 
46  if (!dirName1.empty() && !dirName2.empty()) {
47  TDirectory* dir1 = inputFile1.GetDirectory(dirName1.c_str());
48  TDirectory* dir2 = inputFile2.GetDirectory(dirName2.c_str());
49  h1 = (TH1F*)dir1->Get(hist1Name.c_str());
50  h2 = (TH1F*)dir2->Get(hist2Name.c_str());
51  } else {
52  h1 = (TH1F*)inputFile1.Get(hist1Name.c_str());
53  h2 = (TH1F*)inputFile2.Get(hist2Name.c_str());
54  }
55 
56  h1->Sumw2();
57  h2->Sumw2();
58 
59  int bins = h1->GetXaxis()->GetNbins();
60  float min = h1->GetXaxis()->GetXmin();
61  float max = h1->GetXaxis()->GetXmax();
62  std::cout << "creating new histogram with min: " << min << ", max: " << max
63  << ", nBins: " << bins << std::endl;
64  TH1F* error = new TH1F("error", "relative error", bins, min, max);
65 
66  if (!(bins == h2->GetXaxis()->GetNbins())) {
67  std::cout << "Number of bins of the two histograms not equal : return"
68  << std::endl;
69  return;
70  }
71  for (int i = 0; i < bins; i++) {
72  error->Fill(h1->GetBinCenter(i),
73  (h1->GetBinContent(i) - h2->GetBinContent(i)));
74  }
75  TCanvas* c1 = new TCanvas();
76  gStyle->SetOptStat(0);
77  c1->Divide(1, 2);
78  c1->cd(1);
79  h1->SetMarkerColor(1);
80  h1->SetLineColor(1);
81  h1->Draw("");
82  h2->SetMarkerColor(2);
83  h2->SetLineColor(2);
84  h2->Draw("same");
85  TLegend* leg = new TLegend(0.72, 0.696, 0.99, 0.936);
86  leg->AddEntry(h1, legendName1.c_str());
87  leg->AddEntry(h2, legendName2.c_str());
88  leg->Draw();
89  h1->SetDirectory(0);
90  h2->SetDirectory(0);
91  c1->cd(2);
92  error->Divide(h2);
93  error->Draw("");
94  error->SetDirectory(0);
95  inputFile1.Close();
96  inputFile2.Close();
97 }
98 
99 void
101  std::string hist1Name,
102  std::string legendName1,
103  std::string inFile2,
104  std::string hist2Name,
105  std::string legendName2,
106  std::string dirName1 = "",
107  std::string dirName2 = "")
108 {
109  std::cout << "Opening file: " << inFile1 << std::endl;
110  TFile inputFile1(inFile1.c_str());
111  std::cout << "Opening file: " << inFile2 << std::endl;
112  TFile inputFile2(inFile2.c_str());
113  std::cout << "Comparing Histograms: " << hist1Name << " & " << hist2Name
114  << std::endl;
115 
116  TH2F* h1 = nullptr;
117  TH2F* h2 = nullptr;
118 
119  if (!dirName1.empty() && !dirName2.empty()) {
120  TDirectory* dir1 = inputFile1.GetDirectory(dirName1.c_str());
121  TDirectory* dir2 = inputFile2.GetDirectory(dirName2.c_str());
122  h1 = (TH2F*)dir1->Get(hist1Name.c_str());
123  h2 = (TH2F*)dir2->Get(hist2Name.c_str());
124  } else {
125  h1 = (TH2F*)inputFile1.Get(hist1Name.c_str());
126  h2 = (TH2F*)inputFile2.Get(hist2Name.c_str());
127  }
128  h1->Sumw2();
129  h2->Sumw2();
130 
131  int bins1 = h1->GetXaxis()->GetNbins();
132  int bins2 = h1->GetYaxis()->GetNbins();
133  float min1 = h1->GetXaxis()->GetXmin();
134  float max1 = h1->GetXaxis()->GetXmax();
135  float min2 = h1->GetYaxis()->GetXmin();
136  float max2 = h1->GetYaxis()->GetXmax();
137  std::cout << "creating new histogram with min1: " << min1
138  << ", max1: " << max1 << ", nBins1: " << bins1 << ", min2: " << min2
139  << ", max2: " << max2 << ", nBins2: " << bins2 << std::endl;
140 
141  std::string title = "relative error";
142  if (hist1Name == hist2Name)
143  title += " of " + hist1Name + " in " + legendName1 + " and " + legendName2;
144  else
145  title += " of " + hist1Name + " in " + legendName1 + " and " + hist1Name
146  + legendName2;
147  TH2F* error
148  = new TH2F("error", title.c_str(), bins1, min1, max1, bins2, min2, max2);
149 
150  if (!(bins1 == h2->GetXaxis()->GetNbins())
151  && !(bins2 == h2->GetYaxis()->GetNbins())) {
152  std::cout << "Number of bins of the two histograms not equal : return"
153  << std::endl;
154  return;
155  }
156  for (int i = 0; i < bins1; i++) {
157  for (int j = 0; j < bins2; j++) {
158  error->Fill(h1->GetXaxis()->GetBinCenter(i),
159  h1->GetYaxis()->GetBinCenter(j),
160  (h1->GetBinContent(i, j) - h2->GetBinContent(i, j)));
161  }
162  }
163  TCanvas* c1 = new TCanvas();
164  gStyle->SetOptStat(0);
165  error->Divide(h2);
166  error->Draw("colZ");
167 
168  error->SetDirectory(0);
169 
170  inputFile1.Close();
171  inputFile2.Close();
172 }