RootHistogram.h

Go to the documentation of this file.

/***********************************************************************************\
* (c) Copyright 1998-2022 CERN for the benefit of the LHCb and ATLAS collaborations *
*                                                                                   *
* This software is distributed under the terms of the Apache version 2 licence,     *
* copied verbatim in the file "LICENSE".                                            *
*                                                                                   *
* In applying this licence, 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.                                             *
\***********************************************************************************/
#pragma once
 
#include <Gaudi/Accumulators/Histogram.h>
 
#include <type_traits>
 
namespace Gaudi::Accumulators {
 
  constexpr unsigned int NSUMS( unsigned int ND ) { return 1 + ND + ND * ( ND + 1 ) / 2; }
 
  template <typename Arithmetic, atomicity Atomicity, unsigned int ND>
  struct SigmasValueHandler {
    using InputType  = std::conditional_t<ND == 1, Arithmetic, std::array<Arithmetic, ND>>;
    using OutputType = std::array<Arithmetic, NSUMS( ND )>;
    struct OutputTypeTS : std::array<std::atomic<Arithmetic>, NSUMS( ND )> {
      using std::array<std::atomic<Arithmetic>, NSUMS( ND )>::array;
      explicit OutputTypeTS( OutputType const& other ) {
        for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { ( *this )[i] = other[i]; }
      }
      // operator= from non thread safe type
      OutputTypeTS& operator=( OutputType const& other ) {
        for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { ( *this )[i] = other[i]; }
        return *this;
      }
      // automatic conversion to non thread safe type
      operator OutputType() const {
        OutputType out;
        for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { out[i] = ( *this )[i].load( std::memory_order_relaxed ); }
        return out;
      }
    };
    using InternalType = std::conditional_t<Atomicity == atomicity::full, OutputTypeTS, OutputType>;
    static constexpr OutputType getValue( const InternalType& v ) noexcept {
      // Note automatic conversion will happen
      return v;
    };
    static OutputType exchange( InternalType& v, OutputType newv ) noexcept {
      if constexpr ( Atomicity == atomicity::full ) {
        OutputType old;
        for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { old[i] = v[i].exchange( newv[i] ); }
        return old;
      } else {
        return std::exchange( v, newv );
      }
    }
    static constexpr OutputType DefaultValue() { return InternalType{}; }
    static void                 merge( InternalType& a, OutputType b ) noexcept {
      if constexpr ( Atomicity == atomicity::full ) {
        for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { fetch_add( a[i], b[i] ); }
      } else {
        for ( unsigned int i = 0; i < ND * ( ND + 3 ); i++ ) a[i] += b[i];
      }
    }
    static void merge( InternalType& a, InputType b ) noexcept {
      // prepare values to increment internal values
      OutputType diff{};
      diff[0] = 1.;
      if constexpr ( ND == 1 ) {
        // no operator[] for b in this case
        diff[1] = b;
        diff[2] = b * b;
      } else {
        for ( unsigned int i = 0; i < ND; i++ ) diff[i + 1] += b[i];
        unsigned int n = 1 + ND;
        for ( unsigned int i = 0; i < ND; i++ ) {
          for ( unsigned int j = i; j < ND; j++ ) {
            diff[n] = b[i] * b[j];
            n++;
          }
        }
      }
      // Now increase original counter
      if constexpr ( Atomicity == atomicity::full ) {
        for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { fetch_add( a[i], diff[i] ); }
      } else {
        for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) a[i] += diff[i];
      }
    }
  };
 
  template <typename Arithmetic, atomicity Atomicity, unsigned int ND>
  struct SigmaNAccumulator
      : GenericAccumulator<std::array<Arithmetic, ND>, std::array<Arithmetic, NSUMS( ND )>, Atomicity, Identity,
                           Identity, SigmasValueHandler<Arithmetic, Atomicity, ND>> {
    std::array<Arithmetic, NSUMS( ND )> const sums() const { return this->value(); }
  };
 
  template <typename Arithmetic, atomicity Atomicity>
  struct SigmaNAccumulator<Arithmetic, Atomicity, 1>
      : GenericAccumulator<Arithmetic, std::array<Arithmetic, 3>, Atomicity, Identity, Identity,
                           SigmasValueHandler<Arithmetic, Atomicity, 1>> {
    std::array<Arithmetic, 3> const sums() const { return this->value(); }
  };
 
  template <atomicity Atomicity, typename Arithmetic, unsigned int ND>
  class RootHistogramingAccumulatorInternal
      : public HistogramingAccumulatorInternal<Atomicity, HistoInputType<Arithmetic, ND>, unsigned long,
                                               std::integral_constant<int, ND>, IntegralAccumulator> {
 
    using InputType = HistoInputType<Arithmetic, ND>;
    using Parent = HistogramingAccumulatorInternal<Atomicity, InputType, unsigned long, std::integral_constant<int, ND>,
                                                   IntegralAccumulator>;
 
    template <atomicity, typename, unsigned int>
    friend class RootHistogramingAccumulatorInternal;
 
    static_assert( ND <= 3, "Root on supports histogrmas with dimension <= 3" );
 
    struct Proxy {
      Proxy( RootHistogramingAccumulatorInternal& histo, typename InputType::ValueType& v )
          : m_histo( histo ), m_v( v ) {}
      void                                 operator++() { m_histo.update( m_v ); }
      RootHistogramingAccumulatorInternal& m_histo;
      typename InputType::ValueType        m_v;
    };
 
  public:
    using Parent::Parent;
    friend struct Proxy;
 
    [[deprecated( "Use `++h1[x]`, `++h2[{x,y}]`, etc. instead." )]] RootHistogramingAccumulatorInternal&
    operator+=( InputType v ) {
      update( v );
      return *this;
    }
    void reset() {
      m_accumulator.reset();
      Parent::reset();
    }
    template <atomicity ato>
    void mergeAndReset( RootHistogramingAccumulatorInternal<ato, Arithmetic, ND>& other ) {
      m_accumulator.mergeAndReset( other.m_accumulator );
      Parent::mergeAndReset( other );
    }
    [[nodiscard]] auto operator[]( typename InputType::ValueType v ) { return Proxy( *this, v ); }
 
    auto sums2() const { return m_accumulator.sums(); }
 
  private:
    void update( typename InputType::ValueType v ) {
      // Do not accumulate in m_accumulator if we are outside the histo range
      // We mimic here the behavior of ROOT
      if ( v.inAcceptance( this->axis() ) ) { m_accumulator += static_cast<typename InputType::InternalType>( v ); }
      ++Parent::operator[]( v );
    }
    // Accumulator for keeping squared sum of value stored in the histogram and correlation values
    // they get stored in "alphabetical" order, e.g. for ND=3 x^2, xy, xz, y^2, yz, z^2
    SigmaNAccumulator<Arithmetic, Atomicity, ND> m_accumulator;
  };
 
  namespace {
    template <typename Arithmetic>
    Arithmetic stddev( Arithmetic n, Arithmetic s, Arithmetic s2 ) {
      using Gaudi::Accumulators::sqrt;
      using std::sqrt;
      auto v = ( n > 0 ) ? ( ( s2 - s * ( s / n ) ) / n ) : Arithmetic{};
      return ( Arithmetic{ 0 } > v ) ? Arithmetic{} : sqrt( v );
    }
  } // namespace
 
  template <atomicity Atomicity, typename Arithmetic, typename ND>
  struct RootHistogramingAccumulator;
 
  template <atomicity Atomicity, typename Arithmetic>
  struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<int, 1>>
      : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1> {
    using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1>::RootHistogramingAccumulatorInternal;
    using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1>::nEntries;
    Arithmetic nEntries() const { return this->sums2()[0]; }
    Arithmetic sum() const { return this->sums2()[1]; }
    Arithmetic sum2() const { return this->sums2()[2]; }
    Arithmetic mean() const { return sum() / nEntries(); }
    Arithmetic standard_deviation() const { return stddev( nEntries(), sum(), sum2() ); }
  };
 
  template <atomicity Atomicity, typename Arithmetic>
  struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<int, 2>>
      : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2> {
    using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2>::RootHistogramingAccumulatorInternal;
    using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2>::nEntries;
    Arithmetic nEntries() const { return this->sums2()[0]; }
    Arithmetic sumx() const { return this->sums2()[1]; }
    Arithmetic sumy() const { return this->sums2()[2]; }
    Arithmetic sumx2() const { return this->sums2()[3]; }
    Arithmetic sumy2() const { return this->sums2()[5]; }
    Arithmetic sumxy() const { return this->sums2()[4]; }
    Arithmetic meanx() const { return sumx() / nEntries(); }
    Arithmetic meany() const { return sumy() / nEntries(); }
    Arithmetic standard_deviationx() const { return stddev( nEntries(), sumx(), sumx2() ); }
    Arithmetic standard_deviationy() const { return stddev( nEntries(), sumy(), sumy2() ); }
  };
 
  template <atomicity Atomicity, typename Arithmetic>
  struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<int, 3>>
      : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3> {
    using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3>::RootHistogramingAccumulatorInternal;
    using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3>::nEntries;
    Arithmetic nEntries() const { return this->sums2()[0]; }
    Arithmetic sumx() const { return this->sums2()[1]; }
    Arithmetic sumy() const { return this->sums2()[2]; }
    Arithmetic sumz() const { return this->sums2()[3]; }
    Arithmetic sumx2() const { return this->sums2()[4]; }
    Arithmetic sumy2() const { return this->sums2()[7]; }
    Arithmetic sumz2() const { return this->sums2()[9]; }
    Arithmetic sumxy() const { return this->sums2()[5]; }
    Arithmetic sumxz() const { return this->sums2()[6]; }
    Arithmetic sumyz() const { return this->sums2()[8]; }
    Arithmetic meanx() const { return sumx() / nEntries(); }
    Arithmetic meany() const { return sumy() / nEntries(); }
    Arithmetic meanz() const { return sumz() / nEntries(); }
    Arithmetic standard_deviationx() const { return stddev( nEntries(), sumx(), sumx2() ); }
    Arithmetic standard_deviationy() const { return stddev( nEntries(), sumy(), sumy2() ); }
    Arithmetic standard_deviationz() const { return stddev( nEntries(), sumz(), sumz2() ); }
  };
 
  template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type>
  class RootHistogramingCounterBase;
 
  template <atomicity Atomicity, typename Arithmetic, const char* Type>
  class RootHistogramingCounterBase<1, Atomicity, Arithmetic, Type>
      : public HistogramingCounterBaseInternal<1, Atomicity, Arithmetic, Type, RootHistogramingAccumulator,
                                               std::make_index_sequence<1>> {
  public:
    using Parent = HistogramingCounterBaseInternal<1, Atomicity, Arithmetic, Type, RootHistogramingAccumulator,
                                                   std::make_index_sequence<1>>;
    using Parent::Parent;
 
    friend void to_json( nlohmann::json& j, RootHistogramingCounterBase<1, Atomicity, Arithmetic, Type> const& h ) {
      to_json( j, static_cast<Parent const&>( h ) );
      j["nTotEntries"]        = h.nEntries();
      j["sum"]                = h.sum();
      j["mean"]               = h.mean();
      j["sum2"]               = h.sum2();
      j["standard_deviation"] = h.standard_deviation();
    }
  };
 
  template <atomicity Atomicity, typename Arithmetic, const char* Type>
  class RootHistogramingCounterBase<2, Atomicity, Arithmetic, Type>
      : public HistogramingCounterBaseInternal<2, Atomicity, Arithmetic, Type, RootHistogramingAccumulator,
                                               std::make_index_sequence<2>> {
  public:
    using Parent = HistogramingCounterBaseInternal<2, Atomicity, Arithmetic, Type, RootHistogramingAccumulator,
                                                   std::make_index_sequence<2>>;
    using Parent::Parent;
 
    friend void to_json( nlohmann::json& j, RootHistogramingCounterBase<2, Atomicity, Arithmetic, Type> const& h ) {
      to_json( j, static_cast<Parent const&>( h ) );
      j["nTotEntries"]         = h.nEntries();
      j["sumx"]                = h.sumx();
      j["sumy"]                = h.sumy();
      j["meanx"]               = h.meanx();
      j["meany"]               = h.meany();
      j["sumx2"]               = h.sumx2();
      j["sumy2"]               = h.sumy2();
      j["sumxy"]               = h.sumxy();
      j["standard_deviationx"] = h.standard_deviationx();
      j["standard_deviationy"] = h.standard_deviationy();
    }
  };
 
  template <atomicity Atomicity, typename Arithmetic, const char* Type>
  class RootHistogramingCounterBase<3, Atomicity, Arithmetic, Type>
      : public HistogramingCounterBaseInternal<3, Atomicity, Arithmetic, Type, RootHistogramingAccumulator,
                                               std::make_index_sequence<3>> {
  public:
    using Parent = HistogramingCounterBaseInternal<3, Atomicity, Arithmetic, Type, RootHistogramingAccumulator,
                                                   std::make_index_sequence<3>>;
    using Parent::Parent;
 
    friend void to_json( nlohmann::json& j, RootHistogramingCounterBase<3, Atomicity, Arithmetic, Type> const& h ) {
      to_json( j, static_cast<Parent const&>( h ) );
      j["nTotEntries"]         = h.nEntries();
      j["sumx"]                = h.sumx();
      j["sumy"]                = h.sumy();
      j["sumz"]                = h.sumz();
      j["meanx"]               = h.meanx();
      j["meany"]               = h.meany();
      j["meanz"]               = h.meanz();
      j["sumx2"]               = h.sumx2();
      j["sumy2"]               = h.sumy2();
      j["sumz2"]               = h.sumz2();
      j["sumxy"]               = h.sumxy();
      j["sumxz"]               = h.sumxz();
      j["sumyz"]               = h.sumyz();
      j["standard_deviationx"] = h.standard_deviationx();
      j["standard_deviationy"] = h.standard_deviationy();
      j["standard_deviationz"] = h.standard_deviationz();
    }
  };
 
  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
  using RootHistogram = RootHistogramingCounterBase<ND, Atomicity, Arithmetic, naming::histogramString>;
 
} // namespace Gaudi::Accumulators