Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:12:11

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 /*
0010  * errorOfHistograms.C
0011  *
0012  *  Created on: 24 Jan 2017
0013  *      Author: jhrdinka
0014  */
0015 
0016 #include <iostream>
0017 #include "TFile.h"
0018 #include "TH1F.h"
0019 #include "TROOT.h"
0020 
0021 // This root script compares either two 1D or two 2D histograms with each other,
0022 // by plotting them in different colors into the same canvas and displaying the
0023 // relative error beneath. The histograms to be compared can either be directly
0024 // in the file or in a directory.
0025 
0026 void
0027 errorOf1DHistograms(std::string inFile1,
0028                     std::string hist1Name,
0029                     std::string legendName1,
0030                     std::string inFile2,
0031                     std::string hist2Name,
0032                     std::string legendName2,
0033                     std::string dirName1 = "",
0034                     std::string dirName2 = "")
0035 {
0036   std::cout << "Opening file: " << inFile1 << std::endl;
0037   TFile inputFile1(inFile1.c_str());
0038   std::cout << "Opening file: " << inFile2 << std::endl;
0039   TFile inputFile2(inFile2.c_str());
0040   std::cout << "Comparing Histograms: " << hist1Name << " & " << hist2Name
0041             << std::endl;
0042 
0043   TH1F* h1 = nullptr;
0044   TH1F* h2 = nullptr;
0045 
0046   if (!dirName1.empty() && !dirName2.empty()) {
0047     TDirectory* dir1 = inputFile1.GetDirectory(dirName1.c_str());
0048     TDirectory* dir2 = inputFile2.GetDirectory(dirName2.c_str());
0049     h1               = (TH1F*)dir1->Get(hist1Name.c_str());
0050     h2               = (TH1F*)dir2->Get(hist2Name.c_str());
0051   } else {
0052     h1 = (TH1F*)inputFile1.Get(hist1Name.c_str());
0053     h2 = (TH1F*)inputFile2.Get(hist2Name.c_str());
0054   }
0055 
0056   h1->Sumw2();
0057   h2->Sumw2();
0058 
0059   int   bins = h1->GetXaxis()->GetNbins();
0060   float min  = h1->GetXaxis()->GetXmin();
0061   float max  = h1->GetXaxis()->GetXmax();
0062   std::cout << "creating new histogram with min: " << min << ", max: " << max
0063             << ", nBins: " << bins << std::endl;
0064   TH1F* error = new TH1F("error", "relative error", bins, min, max);
0065 
0066   if (!(bins == h2->GetXaxis()->GetNbins())) {
0067     std::cout << "Number of bins of the two histograms not equal : return"
0068               << std::endl;
0069     return;
0070   }
0071   for (int i = 0; i < bins; i++) {
0072     error->Fill(h1->GetBinCenter(i),
0073                 (h1->GetBinContent(i) - h2->GetBinContent(i)));
0074   }
0075   TCanvas* c1 = new TCanvas();
0076   gStyle->SetOptStat(0);
0077   c1->Divide(1, 2);
0078   c1->cd(1);
0079   h1->SetMarkerColor(1);
0080   h1->SetLineColor(1);
0081   h1->Draw("");
0082   h2->SetMarkerColor(2);
0083   h2->SetLineColor(2);
0084   h2->Draw("same");
0085   TLegend* leg = new TLegend(0.72, 0.696, 0.99, 0.936);
0086   leg->AddEntry(h1, legendName1.c_str());
0087   leg->AddEntry(h2, legendName2.c_str());
0088   leg->Draw();
0089   h1->SetDirectory(0);
0090   h2->SetDirectory(0);
0091   c1->cd(2);
0092   error->Divide(h2);
0093   error->Draw("");
0094   error->SetDirectory(0);
0095   inputFile1.Close();
0096   inputFile2.Close();
0097 }
0098 
0099 void
0100 errorOf2DHistograms(std::string inFile1,
0101                     std::string hist1Name,
0102                     std::string legendName1,
0103                     std::string inFile2,
0104                     std::string hist2Name,
0105                     std::string legendName2,
0106                     std::string dirName1 = "",
0107                     std::string dirName2 = "")
0108 {
0109   std::cout << "Opening file: " << inFile1 << std::endl;
0110   TFile inputFile1(inFile1.c_str());
0111   std::cout << "Opening file: " << inFile2 << std::endl;
0112   TFile inputFile2(inFile2.c_str());
0113   std::cout << "Comparing Histograms: " << hist1Name << " & " << hist2Name
0114             << std::endl;
0115 
0116   TH2F* h1 = nullptr;
0117   TH2F* h2 = nullptr;
0118 
0119   if (!dirName1.empty() && !dirName2.empty()) {
0120     TDirectory* dir1 = inputFile1.GetDirectory(dirName1.c_str());
0121     TDirectory* dir2 = inputFile2.GetDirectory(dirName2.c_str());
0122     h1               = (TH2F*)dir1->Get(hist1Name.c_str());
0123     h2               = (TH2F*)dir2->Get(hist2Name.c_str());
0124   } else {
0125     h1 = (TH2F*)inputFile1.Get(hist1Name.c_str());
0126     h2 = (TH2F*)inputFile2.Get(hist2Name.c_str());
0127   }
0128   h1->Sumw2();
0129   h2->Sumw2();
0130 
0131   int   bins1 = h1->GetXaxis()->GetNbins();
0132   int   bins2 = h1->GetYaxis()->GetNbins();
0133   float min1  = h1->GetXaxis()->GetXmin();
0134   float max1  = h1->GetXaxis()->GetXmax();
0135   float min2  = h1->GetYaxis()->GetXmin();
0136   float max2  = h1->GetYaxis()->GetXmax();
0137   std::cout << "creating new histogram with min1: " << min1
0138             << ", max1: " << max1 << ", nBins1: " << bins1 << ", min2: " << min2
0139             << ", max2: " << max2 << ", nBins2: " << bins2 << std::endl;
0140 
0141   std::string title = "relative error";
0142   if (hist1Name == hist2Name)
0143     title += " of " + hist1Name + " in " + legendName1 + " and " + legendName2;
0144   else
0145     title += " of " + hist1Name + " in " + legendName1 + " and " + hist1Name
0146         + legendName2;
0147   TH2F* error
0148       = new TH2F("error", title.c_str(), bins1, min1, max1, bins2, min2, max2);
0149 
0150   if (!(bins1 == h2->GetXaxis()->GetNbins())
0151       && !(bins2 == h2->GetYaxis()->GetNbins())) {
0152     std::cout << "Number of bins of the two histograms not equal : return"
0153               << std::endl;
0154     return;
0155   }
0156   for (int i = 0; i < bins1; i++) {
0157     for (int j = 0; j < bins2; j++) {
0158       error->Fill(h1->GetXaxis()->GetBinCenter(i),
0159                   h1->GetYaxis()->GetBinCenter(j),
0160                   (h1->GetBinContent(i, j) - h2->GetBinContent(i, j)));
0161     }
0162   }
0163   TCanvas* c1 = new TCanvas();
0164   gStyle->SetOptStat(0);
0165   error->Divide(h2);
0166   error->Draw("colZ");
0167 
0168   error->SetDirectory(0);
0169 
0170   inputFile1.Close();
0171   inputFile2.Close();
0172 }