FastHisto.h

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#ifndef AliceO2_TPC_FastHisto_H
#define AliceO2_TPC_FastHisto_H
 
// o2 includes
#include "MathUtils/fit.h"
#include "Framework/Logger.h"
#include <sstream>
#include <vector>
#include <iomanip>
#include "Rtypes.h"
 
namespace o2
{
namespace tpc
{
 
 
template <class T>
class FastHisto
{
 
 public:
  FastHisto(const unsigned int nBins = 20, const float xmin = 0.f, const float xmax = 2.f, const bool useUnderflow = true, const bool useOverflow = true)
    : mNBins(nBins), mXmin(xmin), mXmax(xmax), mUseUnderflow(useUnderflow), mUseOverflow(useOverflow), mBinCont(nBins + useUnderflow + useOverflow), mBinWidth((mXmax - mXmin) / mNBins){};
 
  ~FastHisto() = default;
 
  void fill(const float val, T weight = 1);
 
  int findBin(const float val) const;
 
  void reset()
  {
    mBinCount = 0;
    std::fill(mBinCont.begin(), mBinCont.end(), 0);
  }
 
  void fillBin(int index, T weight)
  {
    ++mBinCount;
    mBinCont[index] += weight;
  }
 
  // TODO function prints only integers, add to print histgramm with 0.x values
  void print(const int prec = 2) const;
 
  math_utils::StatisticsData getStatisticsData(const float low = 0.05f, const float high = 0.6f) const;
 
  T getBinContent(unsigned int index) const
  {
    if (checkBin(index)) {
      return mBinCont[index];
    } else {
      return -1;
    }
  }
 
  float getBinCenter(int bin) const
  {
    const float binWidth = getBinWidth();
    const float val = (bin - mUseUnderflow) * binWidth + mXmin;
    return val;
  }
 
  float getBinWidth() const { return mBinWidth; }
 
  unsigned int getNBins() const { return mNBins; }
 
  float getXmin() const { return mXmin; }
 
  float getXmax() const { return mXmax; }
 
  unsigned int getEntries() const { return mBinCount; }
 
  bool isUnderflowSet() const { return mUseUnderflow; }
 
  bool isOverflowSet() const { return mUseOverflow; }
 
  bool checkBin(int bin) const
  {
    unsigned int vecsize = mBinCont.size();
    if (bin >= vecsize) {
      return 0;
    } else {
      return 1;
    }
  }
 
  const FastHisto& operator+=(const FastHisto& other);
 
 private:
  unsigned int mNBins{};     
  float mXmin{};             
  float mXmax{};             
  bool mUseUnderflow = true; 
  bool mUseOverflow = true;  
  float mBinWidth{};         
  std::vector<T> mBinCont{}; 
  unsigned int mBinCount{0}; 
 
  ClassDefNV(FastHisto, 1)
};
 
//______________________________________________________________________________
template <class T>
inline void FastHisto<T>::fill(const float val, T weight)
{
  const int indexBin = findBin(val);
  if ((indexBin < 0) || (indexBin > static_cast<int>(mBinCont.size()) - 1)) { // if no underflow/overflow bin is used, but the value should be in the underflow/overflow bin return
    return;
  }
  fillBin(indexBin, weight); // fill the correct index with given weight
}
 
template <class T>
inline void FastHisto<T>::print(const int prec) const
{
  const math_utils::StatisticsData data = getStatisticsData();
  LOGP(info, "\n Entries: {}", mBinCount);
  LOGP(info, "Truncated Mean: {}", data.mCOG);
  LOGP(info, "Standard Deviation: {}", data.mStdDev);
  LOGP(info, "sum of content: {}", data.mSum);
 
  const int maxEle = *std::max_element(mBinCont.begin(), mBinCont.end()); // get the maximum element in the
                                                                          // histogram needed for printing
  // this loop prints the histogram
  // starting from the largest value in the array and go backwards to 0
  std::stringstream stream;
  stream.width(10);
 
  // getting the x axis label
  stream << "  1";
  for (int i = 2; i <= maxEle; ++i) {
    stream.width(2);
    stream << std::right << i;
  }
  LOGP(info, "{}", stream.str());
  stream.str(std::string());
  stream << "---------";
  for (int i = 0; i <= maxEle; ++i) {
    stream << "--";
  }
  LOGP(info, "{}", stream.str());
 
  const float binWidth = getBinWidth();
  const int totalBins = mBinCont.size();
 
  for (int i = 0; i < totalBins; ++i) {
    stream.str(std::string());
    if (i == 0 && mUseUnderflow) {
      stream.width(9);
      stream << std::right << "U_Flow | ";
    } else if (i == totalBins - 1 && mUseOverflow) {
      stream.width(9);
      stream << std::right << "O_Flow | ";
    } else {
      stream.width(6);
      const float xPos = mXmin + binWidth * (i - mUseUnderflow);
      stream << std::fixed << std::setprecision(prec) << xPos;
      const std::string sIn = " | ";
      stream << sIn;
    }
 
    for (int j = 1; j <= getBinContent(i); ++j) {
      const std::string sIn = "x ";
      stream << sIn;
    }
    std::string xPosString = stream.str();
    LOGP(info, "{}", xPosString);
  }
  LOGP(info, "");
}
 
template <class T>
inline math_utils::StatisticsData FastHisto<T>::getStatisticsData(const float low, const float high) const
{
  math_utils::StatisticsData data{};
  // in case something went wrong the COG is the histogram lower limit
  data.mCOG = mXmin;
  if (mBinCont.size() == 0) {
    return data;
  }
  if (low > high) {
    return data;
  }
 
  // total histogram content
  const T integral = std::accumulate(mBinCont.begin(), mBinCont.end(), 0);
  if (integral == 0) {
    return data;
  }
  const float lowerBound = integral * low;
  const float upperBound = integral * high;
  const unsigned int lastBin = mBinCont.size() - 1;
 
  float mean = 0;
  float sum = 0;
  float rms2 = 0;
 
  // add fractional values
  const float binWidth = getBinWidth();
  const int shiftBin = mUseUnderflow ? -1 : 0; // shift the bin position in case of additional bin on the left
 
  T countContent = 0;
  bool isLowerFrac = false;
  unsigned int bin = 0;
  for (bin = 0; bin <= lastBin; ++bin) {
    countContent += mBinCont[bin];
 
    if (countContent - lowerBound <= 0) { // lower truncation cut
      continue;
    } else if (countContent >= upperBound) { // upper truncation cut
      break;
    }
 
    const float xcenter = mXmin + (bin + 0.5f + shiftBin) * binWidth;
    const float tmpBinCont = isLowerFrac ? mBinCont[bin] : countContent - lowerBound; // set bincontent to countcontent for first time only
    const float tmpMean = tmpBinCont * xcenter;
    mean += tmpMean;
    sum += tmpBinCont;
    rms2 += tmpMean * xcenter;
    isLowerFrac = true;
  }
 
  if (!checkBin(bin)) {
    return data;
  }
 
  // set last bin
  // TODO move to upper loop
  const float xcenter = mXmin + (bin + 0.5f + shiftBin) * binWidth;
  const T upFrac = mBinCont[bin] - (static_cast<float>(countContent) - upperBound);
  const float tmpMean = upFrac * xcenter;
  mean += tmpMean;
  sum += upFrac;
  rms2 += tmpMean * xcenter;
 
  if (sum == 0) {
    return data;
  }
 
  mean /= sum;
  data.mCOG = mean;
  rms2 /= sum;
  data.mStdDev = std::sqrt(std::abs(rms2 - mean * mean));
  data.mSum = sum;
  return data;
}
 
template <class T>
inline int FastHisto<T>::findBin(const float val) const
{
  if (val < mXmin) {
    if (!mUseUnderflow) {
      LOGP(warning, "findBin: UNDERFLOW BIN: BIN IS NOT IN HISTOGRAM RANGE!");
      return -1; // if undeflow bin is not used BUT value is in underflow bin return -1
    } else {
      return 0; // if underflow bin is used return 0 as index
    }
  }
 
  if (val >= mXmax) {
    if (!mUseOverflow) {
      LOGP(warning, "findBin: OVERFLOW BIN: BIN IS NOT IN HISTOGRAM RANGE!");
      return -1; // if overflow bin is not used BUT value is in overflow bin return -1
    } else {
      return mBinCont.size() - 1; // if overflow bin is used return the last index in the vector
    }
  }
 
  const float binWidth = getBinWidth();
  const int bin = (val - mXmin) / binWidth + mUseUnderflow;
  return bin;
};
 
template <class T>
inline const FastHisto<T>& FastHisto<T>::operator+=(const FastHisto& other)
{
  if (other.mBinCount == 0) {
    return *this;
  }
 
  // make sure the calibration objects have the same substructure
  if (mNBins != other.mNBins || mXmin != other.mXmin || mXmax != other.mXmax || mUseUnderflow != other.mUseUnderflow || mUseOverflow != other.mUseOverflow) {
    static int errCount = 0;
    if (mBinCount && errCount < 10) {
      errCount++;
      LOGP(warning, "mBinCount {} other.mBinCount: {} mNBins {}, other.mNBins {}, mXmin {}, other.mXmin {}, mXmax {}, other.mXmax {}, mUseUnderflow {}, other.mUseUnderflow {}, mUseOverflow {}, other.mUseOverflow {}", mBinCount, other.mBinCount, mNBins, other.mNBins, mXmin, other.mXmin, mXmax, other.mXmax, mUseUnderflow, other.mUseUnderflow, mUseOverflow, other.mUseOverflow);
    }
    *this = other;
    return *this;
  }
  mBinCount += other.mBinCount;
  std::transform(mBinCont.begin(), mBinCont.end(), other.mBinCont.begin(), mBinCont.begin(), std::plus<T>());
  return *this;
}
 
} // namespace tpc
} // namespace o2
 
#endif