CommonUtils.h

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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;
    }
  };
};

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);
    }
  }
};

struct AcceptCombination {
  std::function<bool(ULong64_t)> one;

  std::function<bool(ULong64_t)> two;

  bool operator()(ULong64_t entry) { return (one(entry) and two(entry)); }
};

struct AcceptAll {
  // Call operator always returns true
  bool operator()(ULong64_t /*event*/) { return true; }
};

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 and range[1] > v);
    }
    return false;
  }
};

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();
  }
};

// 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 and 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();
  }
};

// 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 and reference != nullptr) {
      float v = value->at(entry);
      float r = reference->at(entry);
      return (v - r);
    }
    return std::numeric_limits<float>::infinity();
  }
};

// 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 and reference_charge != nullptr and
        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();
  }
};

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 and theta_value != nullptr and
        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();
  }
};

// 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 and qop_error != nullptr and
        theta_value != nullptr and 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();
  }
};

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 and configuredEntries < entries) {
    entries = configuredEntries;
  }
  return entries;
}