CommonUtils.h

// This file is part of the Acts project.
//
// Copyright (C) 2019-2021 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
#pragma once
 
#include <functional>
#include <limits>
 
#include <TColor.h>
#include <TDirectory.h>
#include <TH1F.h>
#include <TString.h>
 
template <typename hist_t>
void setHistStyle(hist_t* hist, short color = 1) {
  hist->GetXaxis()->SetTitleSize(0.04);
  hist->GetYaxis()->SetTitleSize(0.04);
  hist->GetXaxis()->SetLabelSize(0.04);
  hist->GetYaxis()->SetLabelSize(0.04);
  hist->GetXaxis()->SetTitleOffset(1.0);
  hist->GetYaxis()->SetTitleOffset(1.0);
  hist->GetXaxis()->SetNdivisions(505);
  hist->SetMarkerStyle(20);
  hist->SetMarkerSize(0.8);
  hist->SetLineWidth(2);
  hist->SetTitle("");
  hist->SetLineColor(color);
  hist->SetMarkerColor(color);
}
 
template <typename eff_t>
void setEffStyle(eff_t* eff, short color = 1) {
  eff->SetMarkerStyle(20);
  eff->SetMarkerSize(0.8);
  eff->SetLineWidth(2);
  eff->SetLineColor(color);
  eff->SetMarkerColor(color);
}
 
template <typename hist_t>
void adaptColorPalette(hist_t* h, float rmin, float rmax, float rgood,
                       float rwindow, int n) {
  // min - max is the range of the axis
  float rel_good = (rgood - rmin) / (rmax - rmin);
  float rel_window = rwindow / (rmax - rmin);
 
  // Stops are
  const int number = 5;
  double red[number] = {0., 0., 0., 1., 1.};
  double green[number] = {0., 1., 1., 1., 0.};
  double blue[number] = {1., 1., 0., 0., 0.};
  double stops[number] = {0., rel_good - rel_window, rel_good,
                          rel_good + rel_window, 1.};
  h->SetContour(n);
 
  TColor::CreateGradientColorTable(number, stops, red, green, blue, n);
}
 
template <typename eff_t>
void adaptEffRange(eff_t* eff, float minScale = 1, float maxScale = 1.1) {
  gPad->Update();
  auto ymin = gPad->GetUymin();
  auto ymax = gPad->GetUymax();
  auto graph = eff->GetPaintedGraph();
  graph->SetMinimum(ymin * minScale);
  graph->SetMaximum(ymax * maxScale);
  gPad->Modified();
  gPad->Update();
}
 
struct ResidualPullHandle {
  std::string tag = "";
 
  std::string residualStr = "";
  std::string residualUnit = "";
 
  std::string errorStr = "";
 
  std::string rangeDrawStr = "";
  std::string rangeMaxStr = "";
  std::string rangeCutStr = "";
 
  std::array<float, 2> range = {0., 0.};
 
  std::function<float(ULong64_t)> value;
 
  std::function<float(ULong64_t)> error;
 
  std::function<bool(ULong64_t)> accept;
 
  TH1F* rangeHist = nullptr;
 
  TH1F* residualHist = nullptr;
 
  TH1F* pullHist = nullptr;
 
  ULong64_t accepted = 0;
 
  void fill(unsigned int entry) {
    if (accept(entry)) {
      // Access the value, error
      float v = value(entry);
      residualHist->Fill(v);
      pullHist->Fill(v / error(entry));
      // Count the accessor
      ++accepted;
    }
  };
  void fill(unsigned int entry) {…}
};
struct ResidualPullHandle {…};
 
struct SingleHandle {
  std::string tag = "";
 
  std::string label = "";
 
  // Range draw string
  std::string rangeDrawStr = "";
 
  unsigned int bins = 1;
 
  std::array<float, 2> range = {0., 0.};
 
  std::function<float(ULong64_t)> value;
 
  std::function<bool(ULong64_t)> accept;
 
  TH1F* hist = nullptr;
 
  void fill(unsigned int entry) {
    if (accept(entry)) {
      // Access the value, error
      float v = value(entry);
      hist->Fill(v);
    }
  }
  void fill(unsigned int entry) {…}
};
struct SingleHandle {…};
 
struct AcceptCombination {
  std::function<bool(ULong64_t)> one;
 
  std::function<bool(ULong64_t)> two;
 
  bool operator()(ULong64_t entry) { return (one(entry) && two(entry)); }
};
struct AcceptCombination {…};
 
struct AcceptAll {
  // Call operator always returns true
  bool operator()(ULong64_t /*event*/) { return true; }
};
struct AcceptAll {…};
 
struct AcceptRange {
  std::vector<float>* value = nullptr;
 
  std::array<float, 2> range = {0., 0.};
 
  bool operator()(ULong64_t entry) {
    if (value != nullptr) {
      float v = value->at(entry);
      return (range[0] <= v && range[1] > v);
    }
    return false;
  }
  bool operator()(ULong64_t entry) {…}
};
struct AcceptRange {…};
 
template <typename primitive_t>
struct DirectAccessor {
  std::vector<primitive_t>* value = nullptr;
 
  primitive_t operator()(ULong64_t entry) {
    if (value) {
      primitive_t v = value->at(entry);
      return v;
    }
    return std::numeric_limits<primitive_t>::max();
  }
  primitive_t operator()(ULong64_t entry) {…}
};
struct DirectAccessor {…};
 
// Division accessor
template <typename primitive_one_t, typename primitive_two_t>
struct DivisionAccessor {
  std::vector<primitive_one_t>* one = nullptr;
 
  std::vector<primitive_two_t>* two = nullptr;
 
  primitive_one_t operator()(ULong64_t entry) {
    if (one && two) {
      primitive_one_t vo = one->at(entry);
      primitive_two_t vt = two->at(entry);
      return vo / vt;
    }
    return std::numeric_limits<primitive_one_t>::max();
  }
  primitive_one_t operator()(ULong64_t entry) {…}
};
struct DivisionAccessor {…};
 
// This is a residual type accessor
struct ResidualAccessor {
  std::vector<float>* value = nullptr;
 
  std::vector<float>* reference = nullptr;
 
  float operator()(ULong64_t entry) {
    if (value != nullptr && reference != nullptr) {
      float v = value->at(entry);
      float r = reference->at(entry);
      return (v - r);
    }
    return std::numeric_limits<float>::infinity();
  }
  float operator()(ULong64_t entry) {…}
};
struct ResidualAccessor {…};
 
// This is a  dedicated qop residual accessor
struct QopResidualAccessor {
  std::vector<float>* qop_value = nullptr;
 
  std::vector<int>* reference_charge = nullptr;
 
  std::vector<float>* reference_p = nullptr;
 
  float operator()(ULong64_t entry) {
    if (qop_value != nullptr && reference_charge != nullptr &&
        reference_p != nullptr) {
      float v = qop_value->at(entry);
      float q_true = reference_charge->at(entry);
      float p_true = reference_p->at(entry);
      return (v - q_true / p_true);
    }
    return std::numeric_limits<float>::infinity();
  }
  float operator()(ULong64_t entry) {…}
};
struct QopResidualAccessor {…};
 
struct PtResidualAccessor {
  std::vector<float>* qop_value = nullptr;
 
  std::vector<float>* theta_value = nullptr;
 
  std::vector<float>* reference_pt = nullptr;
 
  float operator()(ULong64_t entry) {
    if (qop_value != nullptr && theta_value != nullptr &&
        reference_pt != nullptr) {
      float p = 1. / std::abs(qop_value->at(entry));
      float theta = theta_value->at(entry);
      float pt_true = reference_pt->at(entry);
      return (p * sin(theta) - pt_true);
    }
    return std::numeric_limits<float>::infinity();
  }
  float operator()(ULong64_t entry) {…}
};
struct PtResidualAccessor {…};
 
// This is a dedicated pT error accessor
struct PtErrorAccessor {
  std::vector<float>* qop_value = nullptr;
  std::vector<float>* qop_error = nullptr;
 
  std::vector<float>* theta_value = nullptr;
  std::vector<float>* theta_error = nullptr;
 
  float operator()(ULong64_t entry) {
    if (qop_value != nullptr && qop_error != nullptr &&
        theta_value != nullptr && theta_error != nullptr) {
      float qop_v = qop_value->at(entry);
      float qop_e = qop_error->at(entry);
      float theta_v = theta_value->at(entry);
      float theta_e = theta_error->at(entry);
      return std::cos(theta_v) / qop_v * theta_e -
             std::sin(theta_v) / (qop_v * qop_v) * qop_e;
    }
    return std::numeric_limits<float>::infinity();
  }
  float operator()(ULong64_t entry) {…}
};
struct PtErrorAccessor {…};
 
template <typename dir_t, typename tree_t>
void estimateResiudalRange(ResidualPullHandle& handle, dir_t& directory,
                           tree_t& tree, unsigned long peakEntries,
                           unsigned int hBarcode) {
  // Change into the Directory
  directory.cd();
  TString rangeHist = handle.rangeDrawStr;
  rangeHist += ">>";
  // Hist name snipped
  TString rangeHN = "hrg_";
  rangeHN += hBarcode;
  // Full histogram
  rangeHist += rangeHN;
  rangeHist += handle.rangeMaxStr;
 
  // Do the drawing
  tree.Draw(rangeHist.Data(), handle.rangeCutStr.c_str(), "", peakEntries);
  handle.rangeHist = dynamic_cast<TH1F*>(gDirectory->Get(rangeHN.Data()));
  if (handle.rangeHist != nullptr) {
    float rms = handle.rangeHist->GetRMS();
    handle.range = {-rms, rms};
  }
}
 
template <typename dir_t, typename tree_t>
void estimateIntegerRange(SingleHandle& handle, dir_t& directory, tree_t& tree,
                          unsigned long peakEntries, unsigned int startBins,
                          unsigned int addBins, unsigned int hBarcode) {
  // Change into the Directory
  directory.cd();
  TString rangeHist = handle.rangeDrawStr;
  rangeHist += ">>";
  // Hist name snipped
  TString rangeHN = "hrg_";
  rangeHN += hBarcode;
  // Full histogram
  rangeHist += rangeHN;
  rangeHist += "(";
  rangeHist += startBins;
  rangeHist += ",-0.5,";
  rangeHist += static_cast<float>(startBins - 0.5);
  rangeHist += ")";
 
  unsigned int nBins = startBins;
  // Do the drawing
  tree.Draw(rangeHist.Data(), "", "", peakEntries);
  auto rhist = dynamic_cast<TH1F*>(gDirectory->Get(rangeHN.Data()));
  if (rhist != nullptr) {
    for (unsigned int ib = 1; ib <= startBins; ++ib) {
      if (rhist->GetBinContent(ib) > 0.) {
        nBins = ib;
      }
    }
    handle.bins = (nBins + addBins);
    handle.range = {-0.5, static_cast<float>(handle.bins - 0.5)};
    return;
  }
  handle.bins = (startBins);
  handle.range = {-0.5, static_cast<float>(handle.bins - 0.5)};
}
 
void bookHistograms(ResidualPullHandle& handle, float pullRange,
                    unsigned int hBins, unsigned int hBarcode) {
  // Residual histogram
  TString rName = std::string("res_") + handle.tag;
  rName += hBarcode;
  handle.residualHist =
      new TH1F(rName.Data(), handle.tag.c_str(), hBins,
               pullRange * handle.range[0], pullRange * handle.range[1]);
  std::string xAxisTitle =
      handle.residualStr + std::string(" ") + handle.residualUnit;
  handle.residualHist->GetXaxis()->SetTitle(xAxisTitle.c_str());
  handle.residualHist->GetYaxis()->SetTitle("Entries");
 
  // Pull histogram
  TString pName = std::string("pull_") + handle.tag;
  pName += hBarcode;
  handle.pullHist =
      new TH1F(pName.Data(), (std::string("pull ") + handle.tag).c_str(), hBins,
               -pullRange, pullRange);
  xAxisTitle = std::string("(") + handle.residualStr + std::string(")/") +
               handle.errorStr;
  handle.pullHist->GetXaxis()->SetTitle(xAxisTitle.c_str());
  handle.pullHist->GetYaxis()->SetTitle("Entries");
}
 
template <typename tree_t>
unsigned long estimateEntries(const tree_t& tree,
                              unsigned long configuredEntries) {
  unsigned long entries = static_cast<unsigned long>(tree.GetEntries());
  if (configuredEntries > 0 && configuredEntries < entries) {
    entries = configuredEntries;
  }
  return entries;
}