From 8e18c810a2c58741273cffda8ed25fa910f9fce5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?V=C3=ADt=20Ku=C4=8Dera?= <vit.kucera@cern.ch>
Date: Wed, 22 Nov 2023 13:04:11 +0100
Subject: [PATCH] Improve y-range calculation and plotting

---
 exec/utilitiesPlot.h       | 63 ++++++++++++++++++++++++++++++++++++++
 exec/utilitiesValidation.h | 16 +++++-----
 2 files changed, 72 insertions(+), 7 deletions(-)

diff --git a/exec/utilitiesPlot.h b/exec/utilitiesPlot.h
index 91f8fdc3..0aa9995f 100644
--- a/exec/utilitiesPlot.h
+++ b/exec/utilitiesPlot.h
@@ -16,11 +16,19 @@
 
 void SetCanvas(TCanvas* can, int nPadsX, int nPadsY)
 {
+  if (!can) {
+    Fatal("SetCanvas", "Invalid canvas");
+    return;
+  }
   can->Divide(nPadsX, nPadsY, 0.005, 0.005);
 }
 
 void SetPad(TVirtualPad* pad, bool logScale)
 {
+  if (!pad) {
+    Fatal("SetPad", "Invalid pad");
+    return;
+  }
   pad->SetBottomMargin(0.11);
   pad->SetLeftMargin(0.1);
   pad->SetTopMargin(0.08);
@@ -30,8 +38,59 @@ void SetPad(TVirtualPad* pad, bool logScale)
   }
 }
 
+/// @brief Get maximum and minimum of a histograms with non-negative bin contents and errors
+/// @param his  histogram
+/// @param yMin  variable to set minimum
+/// @param yMax  variable to set maximum
+/// @param onlyPositive  make sure that yMin, yMax are positive
+/// @param includeErrors  consider heights of error bars
+void GetYRange(TH1* his, Float_t& yMin, Float_t& yMax, bool onlyPositive = true, bool includeErrors = true)
+{
+  if (!his) {
+    Fatal("GetYRange", "Invalid histogram");
+    return;
+  }
+  yMin = his->GetMinimum(0);
+  yMax = his->GetMaximum();
+  if (onlyPositive) {
+    if (yMin <= 0.) {
+      yMin = 1.;
+    }
+    if (yMax <= 0.) {
+      yMax = 1.;
+    }
+  }
+  if (includeErrors) {
+    Float_t cont, err, yLow;
+    for (Int_t i = 0; i < his->GetNbinsX(); i++) {
+      cont = his->GetBinContent(i + 1);
+      if (cont <= 0.) {
+        continue;
+      }
+      err = his->GetBinError(i + 1);
+      yLow = cont - err;
+      if (onlyPositive && yLow <= 0.) {
+        yLow = cont;
+      }
+      yMin = std::min(yMin, yLow);
+      yMax = std::max(yMax, cont + err);
+    }
+  }
+}
+
+/// @brief Set plotting properties of a histogram
+/// @param his  histogram
+/// @param yMin  minimum y value to display
+/// @param yMax  maximum y value to display
+/// @param marginLow  margin to keep below yMin (expressed as a fraction of the full y-axis plotting range)
+/// @param marginHigh  margin to keep above yMax (expressed as a fraction of the full y-axis plotting range)
+/// @note The full y-axis plotting range is calculated from yMin, yMax, marginLow, marginHigh, logScale.
 void SetHistogram(TH1* his, Float_t yMin, Float_t yMax, Float_t marginLow, Float_t marginHigh, bool& logScale)
 {
+  if (!his) {
+    Fatal("SetHistogram", "Invalid histogram");
+    return;
+  }
   Float_t textsize = 0.05;
   his->GetYaxis()->SetTitleSize(textsize);
   his->GetXaxis()->SetTitleSize(textsize);
@@ -51,6 +110,10 @@ void SetHistogram(TH1* his, Float_t yMin, Float_t yMax, Float_t marginLow, Float
 
 void SetHistogramStyle(TH1* his, Int_t colour = 1, Int_t markerStyle = 1, Float_t markerSize = 1, Float_t lineWidth = 1)
 {
+  if (!his) {
+    Fatal("SetHistogramStyle", "Invalid histogram");
+    return;
+  }
   his->SetLineColor(colour);
   his->SetLineWidth(lineWidth);
   his->SetMarkerColor(colour);
diff --git a/exec/utilitiesValidation.h b/exec/utilitiesValidation.h
index 538099c5..662ad0f7 100644
--- a/exec/utilitiesValidation.h
+++ b/exec/utilitiesValidation.h
@@ -161,12 +161,15 @@ Int_t MakePlots(const VecSpecVecSpec& vecSpecVecSpec,
       hO2->SetLineWidth(1);
       hAli->SetTitle(Form(";%s;Entries", labelAxis.Data()));
       hAli->GetYaxis()->SetMaxDigits(3);
-      yMin = TMath::Min(hO2->GetMinimum(0), hAli->GetMinimum(0));
-      yMax = TMath::Max(hO2->GetMaximum(), hAli->GetMaximum());
+      Float_t yMinO2, yMaxO2;
+      GetYRange(hO2, yMinO2, yMaxO2, logScaleH, false);
+      GetYRange(hAli, yMin, yMax, logScaleH, false);
+      yMin = TMath::Min(yMinO2, yMin);
+      yMax = TMath::Max(yMaxO2, yMax);
       SetHistogram(hAli, yMin, yMax, marginLow, marginHigh, logScaleH);
       SetPad(padH, logScaleH);
-      hAli->Draw();
-      hO2->Draw("same");
+      hAli->Draw("hist");
+      hO2->Draw("hist same");
       TLegend* legend = new TLegend(0.2, 0.92, 0.82, 1.0);
       legend->SetNColumns(2);
       legend->SetBorderSize(0);
@@ -180,11 +183,10 @@ Int_t MakePlots(const VecSpecVecSpec& vecSpecVecSpec,
         hRatio = reinterpret_cast<TH1F*>(hO2->Clone(Form("hRatio%d", index)));
         hRatio->Divide(hAli);
         hRatio->SetTitle(Form("Entries ratio: %g;%s;O^{2}/Ali", static_cast<float>(nO2) / static_cast<float>(nAli), labelAxis.Data()));
-        yMin = hRatio->GetMinimum(0);
-        yMax = hRatio->GetMaximum();
+        GetYRange(hRatio, yMin, yMax, logScaleR, false);
         SetHistogram(hRatio, yMin, yMax, marginRLow, marginRHigh, logScaleR);
         SetPad(padR, logScaleR);
-        hRatio->Draw();
+        hRatio->Draw("hist");
       }
     }
     canHis->SaveAs(Form("comparison_histos_%s.pdf", nameSpec.Data()));