diff --git a/GaudiCommonSvc/src/HistogramPersistencySvc/RootHistogramSink.cpp b/GaudiCommonSvc/src/HistogramPersistencySvc/RootHistogramSink.cpp index 9667f6dca3e0b4ff9ec901e3c98ae34bb9eb2095..ff5fe3229afaee638e5453c958d87a729b447680 100644 --- a/GaudiCommonSvc/src/HistogramPersistencySvc/RootHistogramSink.cpp +++ b/GaudiCommonSvc/src/HistogramPersistencySvc/RootHistogramSink.cpp @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2023 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -9,6 +9,7 @@ * or submit itself to any jurisdiction. * \***********************************************************************************/ +#include <Gaudi/Accumulators/StaticHistogram.h> #include <Gaudi/Histograming/Sink/Base.h> #include <Gaudi/Histograming/Sink/Utils.h> #include <GaudiKernel/Service.h> @@ -24,29 +25,29 @@ namespace Gaudi::Histograming::Sink { namespace { using namespace std::string_literals; Base::HistoBinRegistry const binRegistry = { - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<1u, Accumulators::atomicity::full, double> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<1u, Accumulators::atomicity::full, double> ) ), &saveProfileHisto<1, Accumulators::atomicity::full, double> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<1u, Accumulators::atomicity::none, double> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<1u, Accumulators::atomicity::none, double> ) ), &saveProfileHisto<1, Accumulators::atomicity::none, double> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<1u, Accumulators::atomicity::full, float> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<1u, Accumulators::atomicity::full, float> ) ), &saveProfileHisto<1, Accumulators::atomicity::full, float> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<1u, Accumulators::atomicity::none, float> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<1u, Accumulators::atomicity::none, float> ) ), &saveProfileHisto<1, Accumulators::atomicity::none, float> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<2u, Accumulators::atomicity::full, double> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<2u, Accumulators::atomicity::full, double> ) ), &saveProfileHisto<2, Accumulators::atomicity::full, double> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<2u, Accumulators::atomicity::none, double> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<2u, Accumulators::atomicity::none, double> ) ), &saveProfileHisto<2, Accumulators::atomicity::none, double> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<2u, Accumulators::atomicity::full, float> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<2u, Accumulators::atomicity::full, float> ) ), &saveProfileHisto<2, Accumulators::atomicity::full, float> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<2u, Accumulators::atomicity::none, float> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<2u, Accumulators::atomicity::none, float> ) ), &saveProfileHisto<2, Accumulators::atomicity::none, float> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<3u, Accumulators::atomicity::full, double> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<3u, Accumulators::atomicity::full, double> ) ), &saveProfileHisto<3, Accumulators::atomicity::full, double> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<3u, Accumulators::atomicity::none, double> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<3u, Accumulators::atomicity::none, double> ) ), &saveProfileHisto<3, Accumulators::atomicity::none, double> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<3u, Accumulators::atomicity::full, float> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<3u, Accumulators::atomicity::full, float> ) ), &saveProfileHisto<3, Accumulators::atomicity::full, float> }, - { std::type_index( typeid( Gaudi::Accumulators::ProfileHistogram<3u, Accumulators::atomicity::none, float> ) ), + { std::type_index( typeid( Accumulators::StaticProfileHistogram<3u, Accumulators::atomicity::none, float> ) ), &saveProfileHisto<3, Accumulators::atomicity::none, float> }, }; Base::HistoRegistry const registry = { diff --git a/GaudiExamples/TinyExperiment/src/GeneratorAlg.cpp b/GaudiExamples/TinyExperiment/src/GeneratorAlg.cpp index 7dc27d933070cef2cc5f4fe9fd4c0580cc497cd9..11d8c9f4e4b0e0f2c087597d6e3b3aac38de1c05 100644 --- a/GaudiExamples/TinyExperiment/src/GeneratorAlg.cpp +++ b/GaudiExamples/TinyExperiment/src/GeneratorAlg.cpp @@ -12,7 +12,7 @@ #include "IRandomGenSvc.h" #include "MCTrack.h" -#include <Gaudi/Accumulators/Histogram.h> +#include <Gaudi/Accumulators.h> #include <Gaudi/Functional/Transformer.h> #include <cmath> diff --git a/GaudiKernel/include/Gaudi/Accumulators/AxisAsProperty.h b/GaudiKernel/include/Gaudi/Accumulators/AxisAsProperty.h new file mode 100644 index 0000000000000000000000000000000000000000..46df385862c4b18c0d49116fba01f163c6855079 --- /dev/null +++ b/GaudiKernel/include/Gaudi/Accumulators/AxisAsProperty.h @@ -0,0 +1,81 @@ +/***********************************************************************************\ +* (c) Copyright 1998-2024 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/StaticHistogram.h> + +#include <Gaudi/Parsers/Factory.h> +#include <GaudiKernel/ToStream.h> + +#include <boost/spirit/include/qi.hpp> + +/** + * This file provides a Grammar for the type Gaudi::Accumulators::Axis + * It allows to use that type from python with a format liks : + * ( nbins, min, max, title ) + * where title can be ommited. + * note that for the parsing, min and max will be parsed as doubles before + * being converted to the actual type used as Arithmetic for Axis, so an + * automatic conversion between the 2 must be provided + */ +namespace Gaudi { + namespace Parsers { + + namespace qi = boost::spirit::qi; + + template <typename Iterator, typename Skipper, typename Arithmetic> + struct AxisGrammar : qi::grammar<Iterator, Gaudi::Accumulators::Axis<Arithmetic>(), qi::locals<char>, Skipper> { + using Axis = Gaudi::Accumulators::Axis<Arithmetic>; + struct StoreNbinsOp { + void operator()( Axis& res, unsigned int const& nBins ) const { res.setNumBins( nBins ); } + }; + struct StoreMinValueOp { + void operator()( Axis& res, Arithmetic const& minValue ) const { res.setMinValue( minValue ); } + }; + struct StoreMaxValueOp { + void operator()( Axis& res, Arithmetic const& maxValue ) const { res.setMaxValue( maxValue ); } + }; + struct StoreTitleOp { + void operator()( Axis& res, std::string const& title ) const { res.setTitle( title ); } + }; + AxisGrammar() : AxisGrammar::base_type( axis ) { + begin = enc::char_( '[' )[qi::_val = ']'] | enc::char_( '(' )[qi::_val = ')']; + end = enc::char_( qi::_r1 ); + core = qi::int_[storeNbins( qi::_val, qi::_1 )] >> "," >> qi::double_[storeMinValue( qi::_val, qi::_1 )] >> + "," >> qi::double_[storeMaxValue( qi::_val, qi::_1 )] >> + -( "," >> title[storeTitle( qi::_val, qi::_1 )] ); + axis = begin[qi::_a = qi::_1] >> core[qi::_val = qi::_1] >> end( qi::_a ); + } + qi::rule<Iterator, Axis(), qi::locals<char>, Skipper> axis; + qi::rule<Iterator, Axis(), Skipper> core; + qi::rule<Iterator, char()> begin; + qi::rule<Iterator, void( char )> end; + StringGrammar<Iterator, Skipper> title; + ph::function<StoreNbinsOp> storeNbins; + ph::function<StoreMinValueOp> storeMinValue; + ph::function<StoreMaxValueOp> storeMaxValue; + ph::function<StoreTitleOp> storeTitle; + }; + + template <typename Iterator, typename Skipper, typename Arithmetic> + struct Grammar_<Iterator, Gaudi::Accumulators::Axis<Arithmetic>, Skipper> { + using Grammar = AxisGrammar<Iterator, Skipper, Arithmetic>; + }; + + // Parse function... nothing special, but it must be done explicitely. + template <typename Arithmetic> + StatusCode parse( Gaudi::Accumulators::Axis<Arithmetic>& result, const std::string& input ) { + return parse_( result, input ); + } + + } // namespace Parsers + +} // namespace Gaudi diff --git a/GaudiKernel/include/Gaudi/Accumulators/Histogram.h b/GaudiKernel/include/Gaudi/Accumulators/Histogram.h index 3bdc9d0c00ec664fa65ab47a1a629e53d3aaa72a..58b5f3974c226faa7ba97eaf8aae7ff0faecc777 100644 --- a/GaudiKernel/include/Gaudi/Accumulators/Histogram.h +++ b/GaudiKernel/include/Gaudi/Accumulators/Histogram.h @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2022 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -10,588 +10,31 @@ \***********************************************************************************/ #pragma once -#include <Gaudi/Accumulators.h> -#include <Gaudi/MonitoringHub.h> -#include <GaudiKernel/HistoDef.h> - -#include <array> -#include <cmath> -#include <fmt/format.h> -#include <nlohmann/json.hpp> -#include <string> -#include <type_traits> -#include <utility> -#include <vector> +#include <Gaudi/Accumulators/AxisAsProperty.h> +#include <Gaudi/Accumulators/HistogramWrapper.h> +#include <Gaudi/Accumulators/StaticHistogram.h> namespace Gaudi::Accumulators { - namespace details { - template <std::size_t, typename T> - using alwaysT = T; - // get a tuple of n types the given Type, or directly the type for n = 1 - template <typename Type, unsigned int ND> - struct GetTuple; - template <typename Type, unsigned int ND> - using GetTuple_t = typename GetTuple<Type, ND>::type; - template <typename Type, unsigned int ND> - struct GetTuple { - using type = - decltype( std::tuple_cat( std::declval<std::tuple<Type>>(), std::declval<GetTuple_t<Type, ND - 1>>() ) ); - }; - template <typename Type> - struct GetTuple<Type, 1> { - using type = std::tuple<Type>; - }; - - inline void requireValidTitle( std::string_view sv ) { - if ( !sv.empty() && ( std::isspace( sv.back() ) || std::isspace( sv.front() ) ) ) { - throw GaudiException( - fmt::format( "Histogram title \'{}\' has whitespace at front or back -- please remove", sv ), - "Gaudi::Accumulators", StatusCode::FAILURE ); - } - } - } // namespace details - - /** - * A functor to extract weight, take a pair (valueTuple, weight) as input - */ - struct ExtractWeight { - template <typename Arithmetic> - constexpr decltype( auto ) operator()( const std::pair<unsigned long, Arithmetic>& v ) const noexcept { - return v.second; - } - }; - - /** - * A Product functor, take a pair (value, weight) as input - */ - struct WeightedProduct { - template <typename Arithmetic> - constexpr decltype( auto ) operator()( const std::pair<Arithmetic, Arithmetic>& v ) const noexcept { - return v.first * v.second; - } - }; - - /** - * A WeightedSquare functor, take a pair (value, weight) as input - */ - struct WeightedSquare { - template <typename Arithmetic> - constexpr decltype( auto ) operator()( const std::pair<Arithmetic, Arithmetic>& v ) const noexcept { - return v.first * v.first * v.second; - } - }; - - /** - * An Adder ValueHandler, taking weight into account and computing a count plus the sum of the weights - * In case of full atomicity, fetch_add or compare_exchange_weak are used for each element, - * that is we do not have full atomicity accross the two elements - */ - template <typename Arithmetic, atomicity Atomicity> - struct WeightedAdder { - using RegularType = std::pair<unsigned long, Arithmetic>; - using AtomicType = std::pair<std::atomic<unsigned long>, std::atomic<Arithmetic>>; - using OutputType = RegularType; - static constexpr bool isAtomic = Atomicity == atomicity::full; - using InternalType = std::conditional_t<isAtomic, AtomicType, OutputType>; - static constexpr OutputType getValue( const InternalType& v ) noexcept { - if constexpr ( isAtomic ) { - return { v.first.load( std::memory_order_relaxed ), v.second.load( std::memory_order_relaxed ) }; - } else { - return v; - } - }; - static RegularType exchange( InternalType& v, RegularType newv ) noexcept { - if constexpr ( isAtomic ) { - return { v.first.exchange( newv.first ), v.second.exchange( newv.second ) }; - } else { - return { std::exchange( v.first, newv.first ), std::exchange( v.second, newv.second ) }; - } - } - static constexpr OutputType DefaultValue() { return { 0, Arithmetic{} }; } - static void merge( InternalType& a, RegularType b ) noexcept { - if constexpr ( isAtomic ) { - fetch_add( a.first, b.first ); - fetch_add( a.second, b.second ); - } else { - a.first += b.first; - a.second += b.second; - } - }; - }; - - /** - * WeightedCountAccumulator. A WeightedCountAccumulator is an Accumulator storing the number of provided values, - * as well as the weighted version of it, aka. the sum of weights. It takes a pair (valueTuple, weight) as input - * @see Gaudi::Accumulators for detailed documentation - */ - template <atomicity Atomicity, typename Arithmetic> - struct WeightedCountAccumulator - : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, std::pair<unsigned long, Arithmetic>, Atomicity, Identity, - ExtractWeight, WeightedAdder<Arithmetic, Atomicity>> { - using Base = GenericAccumulator<std::pair<Arithmetic, Arithmetic>, std::pair<unsigned long, Arithmetic>, Atomicity, - Identity, ExtractWeight, WeightedAdder<Arithmetic, Atomicity>>; - using Base::Base; - using Base::operator+=; - /// overload of operator+= to be able to only give weight and no value - WeightedCountAccumulator operator+=( const Arithmetic weight ) { - *this += { 1ul, weight }; - return *this; - } - unsigned long nEntries() const { return this->rawValue().first; } - Arithmetic sumOfWeights() const { return this->rawValue().second; } - }; - - /** - * WeightedSumAccumulator. A WeightedSumAccumulator is an Accumulator storing a weighted sum of values. - * It takes a pair (valueTuple, weight) and basically sums the product of the last item othe 2 part of its in put pair - * : weight and value - * @see Gaudi::Accumulators for detailed documentation - */ - template <atomicity Atomicity, typename Arithmetic> - struct WeightedSumAccumulator - : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, WeightedProduct> { - using GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, - WeightedProduct>::GenericAccumulator; - Arithmetic sum() const { return this->value(); } - }; - - /** - * WeightedSquareAccumulator. A WeightedSquareAccumulator is an Accumulator storing a weighted sum of squared values. - * It basically takes a pair (value, weight) as input and sums weight*value*value - * @see Gaudi::Accumulators for detailed documentation - */ - template <atomicity Atomicity, typename Arithmetic = double> - struct WeightedSquareAccumulator - : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, WeightedSquare> { - using GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, - WeightedSquare>::GenericAccumulator; - Arithmetic sum2() const { return this->value(); }; - }; - - /** - * WeightedAveragingAccumulator. An WeightedAveragingAccumulator is an Accumulator able to compute an average - * This implementation takes a pair (value, weight) as input - * @see Gaudi::Accumulators for detailed documentation - */ - template <atomicity Atomicity, typename Arithmetic> - using WeightedAveragingAccumulator = - AveragingAccumulatorBase<Atomicity, Arithmetic, WeightedCountAccumulator, WeightedSumAccumulator>; - - /** - * WeightedSigmaAccumulator. A WeightedSigmaAccumulator is an Accumulator able to compute an average and variance/rms - * This implementation takes a pair (value, weight) as input - * @see Gaudi::Accumulators for detailed documentation - */ - template <atomicity Atomicity, typename Arithmetic> - using WeightedSigmaAccumulator = - SigmaAccumulatorBase<Atomicity, Arithmetic, WeightedAveragingAccumulator, WeightedSquareAccumulator>; - - /** - * Definition of an Histogram Axis - */ - template <typename Arithmetic> - struct Axis { - Axis( unsigned int _nBins, Arithmetic _minValue, Arithmetic _maxValue, std::string _title = {}, - std::vector<std::string> _labels = {} ) - : nBins( _nBins ) - , minValue( _minValue ) - , maxValue( _maxValue ) - , title( std::move( _title ) ) - , labels( std::move( _labels ) ) - , ratio( _nBins / ( _maxValue - _minValue ) ) { - details::requireValidTitle( title ); - for ( const auto& s : labels ) details::requireValidTitle( s ); - }; - explicit Axis( Gaudi::Histo1DDef const& def ) - : Axis( (unsigned int)def.bins(), def.lowEdge(), def.highEdge(), def.title() ){}; - /// number of bins for this Axis - unsigned int nBins; - /// min and max values on this axis - Arithmetic minValue, maxValue; - /// title of this axis - std::string title; - /// labels for the bins - std::vector<std::string> labels; - /** - * precomputed ratio to convert a value into bin number - * equal to nBins/(maxValue-minValue). Only used for floating Arithmetic - */ - Arithmetic ratio; - - /// returns the bin number for a given value, ranging from 0 (underflow) to nBins+1 (overflow) - unsigned int index( Arithmetic value ) const { - // In case we use integer as Arithmetic type, we cannot use ratio for computing indices, - // as ratios < 1.0 will simply be 0, so we have to pay the division in such a case - int idx; - if constexpr ( std::is_integral_v<Arithmetic> ) { - idx = ( ( value - minValue ) * nBins / ( maxValue - minValue ) ) + 1; - } else { - idx = std::floor( ( value - minValue ) * ratio ) + 1; - } - return idx < 0 ? 0 : ( (unsigned int)idx > nBins ? nBins + 1 : (unsigned int)idx ); - } - /// says whether the given value is within the range of the axis - bool inAcceptance( Arithmetic value ) const { return value >= minValue && value <= maxValue; } - }; - - /// automatic conversion of the Axis type to json - template <typename Arithmetic> - void to_json( nlohmann::json& j, const Axis<Arithmetic>& axis ) { - j = nlohmann::json{ { "nBins", axis.nBins }, - { "minValue", axis.minValue }, - { "maxValue", axis.maxValue }, - { "title", axis.title } }; - if ( !axis.labels.empty() ) { j["labels"] = axis.labels; } - } - - /// small class used as InputType for regular Histograms - template <typename Arithmetic, unsigned int ND, unsigned int NIndex = ND> - struct HistoInputType : std::array<Arithmetic, ND> { - // allow construction from set of values - template <class... ARGS, typename = typename std::enable_if_t<( sizeof...( ARGS ) == NIndex )>> - HistoInputType( ARGS... args ) : std::array<Arithmetic, ND>{ static_cast<Arithmetic>( args )... } {} - // The change on NIndex == 1 allow to have simpler syntax in that case, that is no tuple of one item - using ValueType = HistoInputType<Arithmetic, NIndex == 1 ? 1 : ND, NIndex>; - using AxisArithmeticType = Arithmetic; - using InternalType = std::array<Arithmetic, ND>; - unsigned int computeIndex( const std::array<Axis<Arithmetic>, NIndex>& axis ) const { - unsigned int index = 0; - for ( unsigned int j = 0; j < NIndex; j++ ) { - unsigned int dim = NIndex - j - 1; - // compute local index for a given dimension - int localIndex = axis[dim].index( ( *this )[dim] ); - // compute global index. Bins are stored in a column first manner - index *= ( axis[dim].nBins + 2 ); - index += localIndex; - } - return index; - } - bool inAcceptance( const std::array<Axis<Arithmetic>, NIndex>& axis ) const { - for ( unsigned int dim = 0; dim < NIndex; dim++ ) { - if ( !axis[dim].inAcceptance( ( *this )[dim] ) ) return false; - } - return true; - } - auto forInternalCounter() { return 1ul; } - template <typename AxisType, long unsigned NAxis> - static unsigned int computeTotNBins( std::array<AxisType, NAxis> axis ) { - unsigned int nTotBins = 1; - for ( unsigned int i = 0; i < NAxis; i++ ) { nTotBins *= ( axis[i].nBins + 2 ); } - return nTotBins; - } - }; - - /// specialization of HistoInputType for ND == 1 in order to have simpler syntax - /// that is, no tuple of one item - template <typename Arithmetic> - class HistoInputType<Arithmetic, 1> { - public: - using ValueType = HistoInputType; - using AxisArithmeticType = Arithmetic; - using InternalType = Arithmetic; - HistoInputType( Arithmetic a ) : value( a ) {} - unsigned int computeIndex( const std::array<Axis<Arithmetic>, 1>& axis ) const { return axis[0].index( value ); } - bool inAcceptance( const std::array<Axis<Arithmetic>, 1>& axis ) const { return axis[0].inAcceptance( value ); } - Arithmetic& operator[]( int ) { return value; } - operator Arithmetic() const { return value; } - auto forInternalCounter() { return 1ul; } - template <typename AxisType> - static unsigned int computeTotNBins( std::array<AxisType, 1> axis ) { - return axis[0].nBins + 2; - } - - private: - Arithmetic value; - }; - - /// small class used as InputType for weighted Histograms - template <typename Arithmetic, unsigned int ND, unsigned int NIndex = ND> - struct WeightedHistoInputType : std::pair<HistoInputType<Arithmetic, ND, NIndex>, Arithmetic> { - // The change on NIndex == 1 allow to have simpler syntax in that case, that is no tuple of one item - using ValueType = HistoInputType<Arithmetic, NIndex == 1 ? 1 : ND, NIndex>; - using AxisArithmeticType = Arithmetic; - using std::pair<HistoInputType<Arithmetic, ND, NIndex>, Arithmetic>::pair; - unsigned int computeIndex( const std::array<Axis<Arithmetic>, NIndex>& axis ) const { - return this->first.computeIndex( axis ); - } - unsigned int inAcceptance( const std::array<Axis<Arithmetic>, NIndex>& axis ) const { - return this->first.inAcceptance( axis ); - } - auto forInternalCounter() { return std::pair( this->first.forInternalCounter(), this->second ); } - template <typename AxisType, long unsigned NAxis> - static unsigned int computeTotNBins( std::array<AxisType, NAxis> axis ) { - return HistoInputType<Arithmetic, ND, NIndex>::computeTotNBins( axis ); - } - }; - - /** - * Internal Accumulator class dealing with Histograming. Templates parameters are : - * - Atomicity : none or full - * - Arithmetic : the arithmetic type used for values filled into the histogram - * - ND : the number of dimensions of the histogram. - * Note that ND is given as an integral_constant as it needs to be a type for the internal template - * business of the Counter (more precisely the Buffer class) - * - InputType : the type given as input of the Histogram. Typically (Weighted)HistoInputType - * - BaseAccumulator : the underlying accumulator used in each bin - * - * This accumulator is simply an array of BaseAccumulator, one per bin. The number of bins is - * the product of the number of bins for each dimension, each of them increased by 2 for storing - * the values under min and above max - */ - template <atomicity Atomicity, typename InputType, typename Arithmetic, typename ND, - template <atomicity Ato, typename Arith> typename BaseAccumulatorT> - class HistogramingAccumulatorInternal { - template <atomicity, typename, typename, typename, template <atomicity, typename> typename> - friend class HistogramingAccumulatorInternal; - - public: - using BaseAccumulator = BaseAccumulatorT<Atomicity, Arithmetic>; - using AxisArithmeticType = typename InputType::AxisArithmeticType; - using AxisType = Axis<AxisArithmeticType>; - template <std::size_t... Is> - HistogramingAccumulatorInternal( details::GetTuple_t<AxisType, ND::value> axis, std::index_sequence<Is...> ) - : m_axis{ std::get<Is>( axis )... } - , m_totNBins{ InputType::computeTotNBins( m_axis ) } - , m_value( new BaseAccumulator[m_totNBins] ) { - reset(); - } - template <atomicity ato> - HistogramingAccumulatorInternal( - construct_empty_t, - const HistogramingAccumulatorInternal<ato, InputType, Arithmetic, ND, BaseAccumulatorT>& other ) - : m_axis( other.m_axis ) - , m_totNBins{ InputType::computeTotNBins( m_axis ) } - , m_value( new BaseAccumulator[m_totNBins] ) { - reset(); - } - [[deprecated( "Use `++h1[x]`, `++h2[{x,y}]`, etc. instead." )]] HistogramingAccumulatorInternal& - operator+=( InputType v ) { - accumulator( v.computeIndex( m_axis ) ) += v.forInternalCounter(); - return *this; - } - void reset() { - for ( unsigned int index = 0; index < m_totNBins; index++ ) accumulator( index ).reset(); - } - template <atomicity ato> - void mergeAndReset( HistogramingAccumulatorInternal<ato, InputType, Arithmetic, ND, BaseAccumulatorT>& other ) { - assert( m_totNBins == other.m_totNBins ); - for ( unsigned int index = 0; index < m_totNBins; index++ ) { - accumulator( index ).mergeAndReset( other.accumulator( index ) ); - } - } - [[nodiscard]] auto operator[]( typename InputType::ValueType v ) { - return Buffer<BaseAccumulatorT, Atomicity, Arithmetic>{ accumulator( v.computeIndex( m_axis ) ) }; - } - - auto& axis() const { return m_axis; } - auto nBins( unsigned int i ) const { return m_axis[i].nBins; } - auto minValue( unsigned int i ) const { return m_axis[i].minValue; } - auto maxValue( unsigned int i ) const { return m_axis[i].maxValue; } - auto binValue( unsigned int i ) const { return accumulator( i ).value(); } - auto nEntries( unsigned int i ) const { return accumulator( i ).nEntries(); } - auto totNBins() const { return m_totNBins; } - - private: - BaseAccumulator& accumulator( unsigned int index ) const { - assert( index < m_totNBins ); - return m_value[index]; - } - /// set of Axis of this Histogram - std::array<AxisType, ND::value> m_axis; - /// total number of bins in this histogram, under and overflow included - unsigned int m_totNBins; - /// Histogram content - std::unique_ptr<BaseAccumulator[]> m_value; - }; - - /** - * Class implementing a regular histogram accumulator - * - * Actually only an alias to HistogramingAccumulatorInternal with proper template parameters - */ - template <atomicity Atomicity, typename Arithmetic, typename ND> - using HistogramingAccumulator = HistogramingAccumulatorInternal<Atomicity, HistoInputType<Arithmetic, ND::value>, - unsigned long, ND, IntegralAccumulator>; - - /** - * Class implementing a weighted histogram accumulator - * - * Actually only an alias to HistogramingAccumulatorInternal with proper template parameters - */ - template <atomicity Atomicity, typename Arithmetic, typename ND> - using WeightedHistogramingAccumulator = - HistogramingAccumulatorInternal<Atomicity, WeightedHistoInputType<Arithmetic, ND::value>, Arithmetic, ND, - WeightedCountAccumulator>; - - /** - * Class implementing a profile histogram accumulator - * - * Actually only an alias to HistogramingAccumulatorInternal with proper template parameters - */ - template <atomicity Atomicity, typename Arithmetic, typename ND> - using ProfileHistogramingAccumulator = - HistogramingAccumulatorInternal<Atomicity, HistoInputType<Arithmetic, ND::value + 1, ND::value>, Arithmetic, ND, - SigmaAccumulator>; - - /** - * Class implementing a weighted profile histogram accumulator - * - * Actually only an alias to HistogramingAccumulatorInternal with proper template parameters - */ - template <atomicity Atomicity, typename Arithmetic, typename ND> - using WeightedProfileHistogramingAccumulator = - HistogramingAccumulatorInternal<Atomicity, WeightedHistoInputType<Arithmetic, ND::value + 1, ND::value>, - Arithmetic, ND, WeightedSigmaAccumulator>; - - /** - * A base counter dealing with Histograms - * - * Main features of that Counter : - * - can be any number of dimensions. The dimension is its first template parameter - * - for each dimension, a triplet of values have to be given at - * construction : nbins, minValue, maxValue. These triplets have to be - * embedded into braces, as the constructor takes an array of them - * - the operator+= takes either an array of values (one per dimension) - * or a tuple<array of values, weight>. The value inside the bin - * corresponding to the given values is then increased by 1 or weight - * - the prefered syntax is to avoid operator+= and use operator[] to get a - * buffer on the bin you're updating. Syntax becomes : - * ++counter[{x,y}] or wcounter[{x,y]] += w - * - the Counter is templated by the types of the values given to - * operator+ and also by the type stored into the bins - * - the counter can be atomic or not and supports buffering. Note that - * the atomicity is classical eventual consistency. So each bin is - * atomically updated but bins are not garanted to be coherent when - * reading all of them back - * - profile histograms are also supported, operator+= takes one more - * value in the array of values in that case - * - * This base class is then aliases for the 4 standard cases of Histogram, - * WeightedHistogram, ProfileHistogram and WeightedProfileHistogram - * - * Typical usage : - * \code - * Histogram<2, double, atomicity::full> - * counter{owner, "CounterName", "HistoTitle", {{nBins1, minVal1, maxVal1}, {nBins2, minVal2, maxVal2}}}; - * ++counter[{val1, val2}]; // prefered syntax - * counter += {val1, val2}; // original syntax inherited from counters - * - * WeightedHistogram<2, double, atomicity::full> - * wcounter{owner, "CounterName", "HistoTitle", {{nBins1, minVal1, maxVal1}, {nBins2, minVal2, maxVal2}}}; - * wcounter[{val1, val2}] += w; // prefered syntax - * wcounter += {{val1, val2}, w}; // original syntax inherited from counters - * \endcode - * - * When serialized to json, this counter uses new types histogram:Histogram:<prec>, histogram:ProfileHistogram:<prec>, - * histogram:WeightedHistogram:<prec> and histrogram:WeightedProfileHistogram:<prec> - * <prec> id described in Accumulators.h and decribes here the precision of the bin data. - * All these types have the same fields, namely : - * dimension(integer), title(string), empty(bool), nEntries(integer), axis(array), bins(array) - * where : - * + axis is an array of tuples, one per dimension, with content (nBins(integer), minValue(number), - * maxValue(number), title(string)) - * + bins is an array of values - * - The length of the array is the product of (nBins+2) for all axis - * - the +2 is because the bin 0 is the one for values below minValue and bin nBins+1 is the one for values - * above maxValue bins are stored row first so we iterate first on highest dimension - * - the value is a number for non profile histograms - * - the value is of the form ( (nEntries(integer), sum(number) ), sum2(number) ) for profile histograms - * Note the pair with a pair as first entry - */ - template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type, - template <atomicity, typename, typename> typename Accumulator, typename Seq> - class HistogramingCounterBaseInternal; - template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type, - template <atomicity, typename, typename> typename Accumulator, std::size_t... NDs> - class HistogramingCounterBaseInternal<ND, Atomicity, Arithmetic, Type, Accumulator, std::index_sequence<NDs...>> - : public BufferableCounter<Atomicity, Accumulator, Arithmetic, std::integral_constant<int, ND>> { - public: - using Parent = BufferableCounter<Atomicity, Accumulator, Arithmetic, std::integral_constant<int, ND>>; - using AccumulatorType = Accumulator<Atomicity, Arithmetic, std::integral_constant<int, ND>>; - using NumberDimensions = std::integral_constant<unsigned int, ND>; - inline static const std::string typeString{ std::string{ Type } + ':' + typeid( Arithmetic ).name() }; - template <typename OWNER> - HistogramingCounterBaseInternal( OWNER* owner, std::string const& name, std::string const& title, - details::GetTuple_t<typename AccumulatorType::AxisType, ND> axis ) - : Parent( owner, name, *this, axis, std::make_index_sequence<ND>{} ), m_title( title ) { - details::requireValidTitle( m_title ); - } - template <typename OWNER> - HistogramingCounterBaseInternal( OWNER* owner, std::string const& name, std::string const& title, - details::alwaysT<NDs, typename AccumulatorType::AxisType>... allAxis ) - : HistogramingCounterBaseInternal( owner, name, title, std::make_tuple( allAxis... ) ) {} - using Parent::print; - template <typename stream> - stream& printImpl( stream& o, bool /*tableFormat*/ ) const { - o << ND << "D Histogram with config "; - for ( unsigned int i = 0; i < ND; i++ ) { - o << this->nBins( i ) << " " << this->minValue( i ) << " " << this->maxValue( i ) << " "; - } - return o; - } - std::ostream& print( std::ostream& o, bool tableFormat = false ) const override { - return printImpl( o, tableFormat ); - } - MsgStream& print( MsgStream& o, bool tableFormat = false ) const override { return printImpl( o, tableFormat ); } - friend void reset( HistogramingCounterBaseInternal& c ) { c.reset(); } - friend void mergeAndReset( HistogramingCounterBaseInternal& h, HistogramingCounterBaseInternal& o ) { - h.mergeAndReset( o ); - } - friend void to_json( nlohmann::json& j, HistogramingCounterBaseInternal const& h ) { h.to_json( j ); } - virtual void to_json( nlohmann::json& j ) const { - // get all bin values and compute total nbEntries - std::vector<typename AccumulatorType::BaseAccumulator::OutputType> bins; - bins.reserve( this->totNBins() ); - unsigned long totNEntries{ 0 }; - for ( unsigned int i = 0; i < this->totNBins(); i++ ) { - bins.push_back( this->binValue( i ) ); - totNEntries += this->nEntries( i ); - } - // build json - j = { { "type", std::string( Type ) + ":" + typeid( Arithmetic ).name() }, - { "title", m_title }, - { "dimension", ND }, - { "empty", totNEntries == 0 }, - { "nEntries", totNEntries }, - { "axis", this->axis() }, - { "bins", bins } }; - } - std::string const& title() const { return m_title; } - - protected: - std::string const m_title; - }; - template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type, - template <atomicity, typename, typename> typename Accumulator> - using HistogramingCounterBase = - HistogramingCounterBaseInternal<ND, Atomicity, Arithmetic, Type, Accumulator, std::make_index_sequence<ND>>; - - namespace naming { - inline constexpr char histogramString[] = "histogram:Histogram"; - inline constexpr char weightedHistogramString[] = "histogram:WeightedHistogram"; - inline constexpr char profilehistogramString[] = "histogram:ProfileHistogram"; - inline constexpr char weightedProfilehistogramString[] = "histogram:WeightedProfileHistogram"; - } // namespace naming - /// standard histograming counter. See HistogramingCounterBase for details - template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double> - using Histogram = - HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::histogramString, HistogramingAccumulator>; + /// standard custom histogram. See HistogramWrapper and StaticHistogram for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using Histogram = HistogramWrapper<StaticHistogram<ND, Atomicity, Arithmetic, AxisTupleType>>; - /// standard histograming counter with weight. See HistogramingCounterBase for details - template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double> - using WeightedHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::weightedHistogramString, - WeightedHistogramingAccumulator>; + /// custom histogram with weight. See HistogramWrapper and StaticWeightedHistogram for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using WeightedHistogram = HistogramWrapper<StaticWeightedHistogram<ND, Atomicity, Arithmetic, AxisTupleType>>; - /// profile histograming counter. See HistogramingCounterBase for details - template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double> - using ProfileHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::profilehistogramString, - ProfileHistogramingAccumulator>; + /// custom profile histograming. See HistogramWrapper and StaticProfileHistogram for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using ProfileHistogram = HistogramWrapper<StaticProfileHistogram<ND, Atomicity, Arithmetic, AxisTupleType>>; - /// weighted profile histograming counter. See HistogramingCounterBase for details - template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double> + /// custom weighted profile histogram. See HistogramWrapper and StaticWeightedProfileHistogram for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> using WeightedProfileHistogram = - HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::weightedProfilehistogramString, - WeightedProfileHistogramingAccumulator>; + HistogramWrapper<StaticWeightedProfileHistogram<ND, Atomicity, Arithmetic, AxisTupleType>>; } // namespace Gaudi::Accumulators diff --git a/GaudiKernel/include/Gaudi/Accumulators/HistogramArray.h b/GaudiKernel/include/Gaudi/Accumulators/HistogramArray.h index 5b06a78dd7d29adda7c0b497588089cb674bc278..e05b1d35aab78afd7749a90b705e2bb714ebeca7 100644 --- a/GaudiKernel/include/Gaudi/Accumulators/HistogramArray.h +++ b/GaudiKernel/include/Gaudi/Accumulators/HistogramArray.h @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2022 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -37,21 +37,21 @@ namespace Gaudi::Accumulators { template <typename Histo, std::size_t N> struct HistogramArrayInternal : std::array<Histo, N> { /// constructor with callables for FormatName and FormatTitle - template <typename OWNER, typename FormatName, typename FormatTitle, std::size_t... Ns, typename... Axis, + template <typename OWNER, typename FormatName, typename FormatTitle, std::size_t... Ns, typename = typename std::enable_if_t<std::is_invocable_v<FormatName, int>>, typename = typename std::enable_if_t<std::is_invocable_v<FormatTitle, int>>> HistogramArrayInternal( OWNER* owner, FormatName&& fname, FormatTitle&& ftitle, - std::integer_sequence<std::size_t, Ns...>, Axis&&... allAxis ) - : std::array<Histo, N>{ Histo{ owner, fname( Ns ), ftitle( Ns ), allAxis... }... } { + std::integer_sequence<std::size_t, Ns...>, typename Histo::AxisTupleType&& allAxis ) + : std::array<Histo, N>{ Histo{ owner, fname( Ns ), ftitle( Ns ), allAxis }... } { static_assert( sizeof...( Ns ) < 1000, "Using HistogramArray with 1000 arrays or more is prohibited. This " "would lead to very long compilation times" ); } /// constructor for strings, FormatHistDefault is used as the default callable - template <typename OWNER, std::size_t... Ns, typename... Axis> + template <typename OWNER, std::size_t... Ns> HistogramArrayInternal( OWNER* owner, std::string_view name, std::string_view title, - std::integer_sequence<std::size_t, Ns...>, Axis&&... allAxis ) + std::integer_sequence<std::size_t, Ns...>, typename Histo::AxisTupleType&& allAxis ) : std::array<Histo, N>{ - Histo{ owner, FormatHistDefault{ name }( Ns ), FormatHistDefault{ title }( Ns ), allAxis... }... } { + Histo{ owner, FormatHistDefault{ name }( Ns ), FormatHistDefault{ title }( Ns ), allAxis }... } { static_assert( sizeof...( Ns ) < 1000, "Using HistogramArray with 1000 arrays or more is prohibited. This " "would lead to very long compilation times" ); } @@ -93,17 +93,23 @@ namespace Gaudi::Accumulators { * { 21, -10.5, 10.5, "X" } }; * ++histoCustom[2][-10.0]; */ - template <typename Histo, std::size_t N, typename Seq> - struct HistogramArrayBase; - template <typename Histo, std::size_t N, std::size_t... NDs> - struct HistogramArrayBase<Histo, N, std::index_sequence<NDs...>> : details::HistogramArrayInternal<Histo, N> { + template <typename Histo, std::size_t N, + typename Seq = std::make_integer_sequence<unsigned int, std::tuple_size_v<typename Histo::AxisTupleType>>> + struct HistogramArray; + template <typename Histo, std::size_t N, unsigned int... ND> + struct HistogramArray<Histo, N, std::integer_sequence<unsigned int, ND...>> + : details::HistogramArrayInternal<Histo, N> { template <typename OWNER, typename FormatName, typename FormatTitle> - HistogramArrayBase( OWNER* owner, FormatName&& fname, FormatTitle&& ftitle, - details::alwaysT<NDs, typename Histo::AxisType>&&... allAxis ) + HistogramArray( OWNER* owner, FormatName&& fname, FormatTitle&& ftitle, typename Histo::AxisTupleType&& allAxis ) : details::HistogramArrayInternal<Histo, N>( owner, fname, ftitle, std::make_integer_sequence<std::size_t, N>{}, - allAxis... ) {} + std::forward<typename Histo::AxisTupleType>( allAxis ) ) {} + + template <unsigned int I> + using AxisType = std::tuple_element_t<I, typename Histo::AxisTupleType>; + + template <typename OWNER, typename FormatName, typename FormatTitle> + HistogramArray( OWNER* owner, FormatName&& fname, FormatTitle&& ftitle, AxisType<ND>... allAxis ) + : HistogramArray( owner, fname, ftitle, std::make_tuple( allAxis... ) ) {} }; - template <typename Histo, std::size_t N> - using HistogramArray = HistogramArrayBase<Histo, N, std::make_index_sequence<Histo::NumberDimensions::value>>; } // namespace Gaudi::Accumulators diff --git a/GaudiKernel/include/Gaudi/Accumulators/HistogramWrapper.h b/GaudiKernel/include/Gaudi/Accumulators/HistogramWrapper.h new file mode 100644 index 0000000000000000000000000000000000000000..caf7beff3e8c99882dcd9b74d3b177e6520d2a39 --- /dev/null +++ b/GaudiKernel/include/Gaudi/Accumulators/HistogramWrapper.h @@ -0,0 +1,157 @@ +/***********************************************************************************\ +* (c) Copyright 1998-2024 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/StaticHistogram.h> + +namespace Gaudi::Accumulators { + + /** + * A Wrapper of a static Histogram base class using Properties to define title and axis + * + * Practically, this is an optional of the underlying static Histogram with creation + * on demand, via the createHistogram method so one can wait that properties' values + * are known. + * By default creation will happen at initialization of the owner, but this can be + * disabled at construction time, in which case the user has to call createHistogram + * manually. + * Note that disabling automatic initialization requires using the plain + * constructor and passing title and axis. They can of course be empty and overwritten + * via Properties. But there was no easier way providing char* convert automagically + * to bool in C++ creating potential confusion between title on doNotInitialize + * + * This wrapper expects the Axis type of the Histogram to be usable in a + * Gaudi Property and thus to define the required Grammar, parse function and operator<< + * It also requires the underlying Histogram type to define : + * - AxisTupleType : a tuple type of all Axis types + * - AxisTupleArithmeticType : a tuple type of all Axis Arithmetic types + * FIXME : use concepts when available + */ + template <typename HistogramType, + typename Seq = + std::make_integer_sequence<unsigned int, std::tuple_size_v<typename HistogramType::AxisTupleType>>> + class HistogramWrapperInternal; + template <typename HistogramType, unsigned int... ND> + class HistogramWrapperInternal<HistogramType, std::integer_sequence<unsigned int, ND...>> { + public: + using AxisTupleType = typename HistogramType::AxisTupleType; + using AxisArithmeticType = typename HistogramType::AxisArithmeticType; + template <unsigned int I> + using AxisType = std::tuple_element_t<I, AxisTupleType>; + + /// constructor, only creates a set of Properties + template <typename OWNER> + HistogramWrapperInternal( OWNER* owner, std::string const& name, std::string const& title = "", + typename HistogramType::AxisTupleType axis = {}, bool doNotInitialize = false ) + : m_name{ name }, m_title{ title }, m_axis{ axis } { + // Create associated properties + owner->declareProperty( titlePropertyName(), m_title, fmt::format( "Title of histogram {}", name ) ) + ->template setOwnerType<OWNER>(); + ( owner + ->declareProperty( axisPropertyName<ND>(), std::get<ND>( m_axis ), + fmt::format( "Axis {} of histogram {}", ND, name ) ) + ->template setOwnerType<OWNER>(), + ... ); + // register creation of the Histogram at initialization time + if ( !doNotInitialize ) { + if ( owner->FSMState() >= Gaudi::StateMachine::INITIALIZED ) { + // if the owner is already initialized (e.g. the histogram is being created during `start()`) + // it is too late to register the callback + createHistogram( *owner ); + } else { + owner->registerCallBack( StateMachine::INITIALIZE, [this, owner]() { createHistogram( *owner ); } ); + } + } + } + /// constructor with const owner (i.e. outside of constructor or initialize), equivalent to + /// the StaticHistogram case (create immediately) + template <typename OWNER> + HistogramWrapperInternal( OWNER const* owner, std::string const& name, std::string const& title = "", + typename HistogramType::AxisTupleType axis = {} ) + : m_name{ name }, m_title{ title }, m_axis{ axis } { + createHistogram( *owner ); + } + /// constructor with more natural syntax for axis + template <typename OWNER> + HistogramWrapperInternal( OWNER* owner, std::string const& name, std::string const& title, AxisType<ND>... allAxis ) + : HistogramWrapperInternal( owner, name, title, std::make_tuple( allAxis... ) ) {} + + /// override of operator[] with extra checking that initialization happened + [[nodiscard]] auto operator[]( typename HistogramType::AxisTupleArithmeticType&& v ) { + if ( !m_histo ) { + throw std::logic_error( fmt::format( "Histogram {} is used before being initialized", m_name ) ); + } + return m_histo.value()[std::forward<typename HistogramType::AxisTupleArithmeticType>( v )]; + } + + /// creation of the internal histogram, from the properties + template <typename OWNER> + void createHistogram( OWNER& owner ) { + m_histo.emplace( &owner, m_name, m_title, m_axis ); + } + + friend void to_json( nlohmann::json& j, HistogramWrapperInternal const& h ) { + if ( !h.m_histo ) { + throw std::logic_error( fmt::format( "Histogram {} is converted to json before being initialized", h.m_name ) ); + } + j = h.m_histo.value(); + } + + // set directly some properties, only if histogram was not yet created + void setTitle( std::string const& title ) { + if ( m_histo ) + throw std::logic_error( + fmt::format( "Cannot modify title of histogram {} after it has been initialized", m_name ) ); + m_title = title; + } + template <unsigned int N> + void setAxis( std::tuple_element_t<N, typename HistogramType::AxisTupleType> const& axis ) { + if ( m_histo ) + throw std::logic_error( + fmt::format( "Cannot modify axis {} of histogram {} after it has been initialized", N, m_name ) ); + std::get<N>( m_axis ) = axis; + } + + void reset() { m_histo.reset(); } + + // wrapping some methods of the underlyoing histogram + auto buffer() { + if ( !m_histo ) + throw std::logic_error( fmt::format( "`buffer()` called on histogram {} before being initialized", m_name ) ); + return m_histo->buffer(); + } + + private: + std::string basePropertyName() const { + // properties are used as identifiers in python and thus cannot anything alse than _, letters and numbers + // we thus replace anything else with '_' in the property names + std::string name = m_name; + std::replace_if( + begin( name ), end( name ), []( auto& c ) { return !std::isalnum( c ); }, '_' ); + return name; + } + std::string titlePropertyName() const { return fmt::format( "{}_Title", basePropertyName() ); } + template <unsigned int N> + std::string axisPropertyName() const { + return fmt::format( "{}_Axis{}", basePropertyName(), N ); + } + + // Members of the custom histogrem + std::string m_name{}; + std::string m_title{}; + typename HistogramType::AxisTupleType m_axis{}; + std::optional<HistogramType> m_histo{}; + }; + + template <typename HistogramType> + using HistogramWrapper = HistogramWrapperInternal<HistogramType>; + +} // namespace Gaudi::Accumulators diff --git a/GaudiKernel/include/Gaudi/Accumulators/RootHistogram.h b/GaudiKernel/include/Gaudi/Accumulators/RootHistogram.h index 3d754e1bed246e57a9e6eb367b574da7e297b6ec..caffcb5c42a21e8ce737672541e7e4fb4378ede8 100644 --- a/GaudiKernel/include/Gaudi/Accumulators/RootHistogram.h +++ b/GaudiKernel/include/Gaudi/Accumulators/RootHistogram.h @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2022 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -10,336 +10,15 @@ \***********************************************************************************/ #pragma once -#include <Gaudi/Accumulators/Histogram.h> - -#include <type_traits> +#include <Gaudi/Accumulators/AxisAsProperty.h> +#include <Gaudi/Accumulators/HistogramWrapper.h> +#include <Gaudi/Accumulators/StaticRootHistogram.h> namespace Gaudi::Accumulators { - /// number of items in sums for a given dimension - /// = 1 (nb items) + ND (sums of each dimension) + ND*(ND+1)/2 (square sums) - 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 )> { - /// copy constructor from non thread safe type - 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(); } - }; - - /// specialization for ND=1 to allow for better syntax - 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(); } - }; - - /** - * Internal Accumulator class dealing with RootHistograming. - * Actually a simple extention on top of RootHistograming with an - * extra SigmaCounter embeded - */ - 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" ); - - /** - * Small procyclass allowing operator[] to work as expected on the RootHistogram - * that is to return something having an operator+= updating the histogram properly - */ - 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 ); } - - /// returns the nbentries, sums and "squared sums" of the inputs - /// Practically we have first th enumber of entries, then the simple sums of each - /// input dimension followed by all combinasions of product of 2 inputs, in "alphabetical" order, - /// e.g. for ND=3 we have sums of n, x, y, z, x^2, xy, xz, y^2, yz, z^2 - 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 { - /// helper function to compute standard_deviation - 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 - - /** - * Class implementing a root histogram accumulator - */ - 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() ); } - }; - - /** - * Extension of the standard Gaudi histogram to provide similar functionnality as ROOT - * - * The main piece of extra functionnality is that ablity to compute a mean and sigma - * of the data set hold by the histogram on top of the histogram itself - * This will then be used when filling a real Root histogram so that one gets the - * same values as if working directly with Root. - * When using pure Gaudi Histograms, and converting to ROOT, ROOT automatically - * recomputed average and sigma form the bins, but the value is not the expected one - * - * Usage is similar to HistogramingCounterBaseInternal so see the documentation there - * Serialization has the following extra fields : - * - nTotEntries, sum, sum2, mean in 1D - * - nTotEntries, sumx, sumy, sumx2, sumy2, sumxy, meanx, meany in 2D - * - nTotEntries, sumx, sumy, sumz, sumx2, sumy2, sumz2, sumxy, sumxz, sumyz, meanx, meany, meanz in 3D - * and uses same types as HistogramingCounterBaseInternal - * - */ - 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(); - } - }; - - /// Root histograming counter. See RootHistogramingCounterBase for details - template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double> - using RootHistogram = RootHistogramingCounterBase<ND, Atomicity, Arithmetic, naming::histogramString>; + /// standard custom histogram + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using RootHistogram = HistogramWrapper<StaticRootHistogram<ND, Atomicity, Arithmetic, AxisTupleType>>; } // namespace Gaudi::Accumulators diff --git a/GaudiKernel/include/Gaudi/Accumulators/StaticHistogram.h b/GaudiKernel/include/Gaudi/Accumulators/StaticHistogram.h new file mode 100644 index 0000000000000000000000000000000000000000..52894d1373ba65903a0b1660be0312e4dcc7bbbb --- /dev/null +++ b/GaudiKernel/include/Gaudi/Accumulators/StaticHistogram.h @@ -0,0 +1,736 @@ +/***********************************************************************************\ +* (c) Copyright 1998-2024 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.h> +#include <Gaudi/MonitoringHub.h> +#include <GaudiKernel/HistoDef.h> + +#include <array> +#include <cmath> +#include <fmt/format.h> +#include <nlohmann/json.hpp> +#include <string> +#include <type_traits> +#include <utility> +#include <vector> + +namespace { + // Helper class creating a "subtuple" type from a tuple type by keeping only + // the first N items. + template <typename Tuple, typename Seq> + struct SubTuple; + template <typename Tuple, size_t... I> + struct SubTuple<Tuple, std::index_sequence<I...>> { + using type = decltype( std::make_tuple( std::get<I>( std::declval<Tuple>() )... ) ); + }; + template <typename Tuple, unsigned int N> + using SubTuple_t = typename SubTuple<Tuple, std::make_index_sequence<N>>::type; + + /// helper class to create a tuple of N identical types + template <typename T, unsigned int ND, typename = std::make_integer_sequence<unsigned int, ND>> + struct make_tuple; + template <typename T, unsigned int ND, unsigned int... S> + struct make_tuple<T, ND, std::integer_sequence<unsigned int, S...>> { + template <unsigned int> + using typeMap = T; + using type = std::tuple<typeMap<S>...>; + }; + template <typename T, unsigned int ND> + using make_tuple_t = typename make_tuple<T, ND>::type; + + /// template magic converting a tuple of Axis into the tuple of corresponding Arithmetic types + template <typename AxisTupleType> + struct AxisToArithmetic; + template <typename... Axis> + struct AxisToArithmetic<std::tuple<Axis...>> { + using type = std::tuple<typename Axis::ArithmeticType...>; + }; + template <typename AxisTupleType> + using AxisToArithmetic_t = typename AxisToArithmetic<AxisTupleType>::type; + template <typename ProfArithmetic, typename AxisTupleType> + using ProfileAxisToArithmetic_t = decltype( std::tuple_cat( std::declval<AxisToArithmetic_t<AxisTupleType>>(), + std::declval<std::tuple<ProfArithmetic>>() ) ); +} // namespace + +namespace Gaudi::Accumulators { + + namespace details { + inline void requireValidTitle( std::string_view sv ) { + if ( !sv.empty() && ( std::isspace( sv.back() ) || std::isspace( sv.front() ) ) ) { + throw GaudiException( + fmt::format( "Histogram title \'{}\' has whitespace at front or back -- please remove", sv ), + "Gaudi::Accumulators", StatusCode::FAILURE ); + } + } + } // namespace details + + /** + * A functor to extract weight, take a pair (valueTuple, weight) as input + */ + struct ExtractWeight { + template <typename Arithmetic> + constexpr decltype( auto ) operator()( const std::pair<unsigned long, Arithmetic>& v ) const noexcept { + return v.second; + } + }; + + /** + * A Product functor, take a pair (value, weight) as input + */ + struct WeightedProduct { + template <typename Arithmetic> + constexpr decltype( auto ) operator()( const std::pair<Arithmetic, Arithmetic>& v ) const noexcept { + return v.first * v.second; + } + }; + + /** + * A WeightedSquare functor, take a pair (value, weight) as input + */ + struct WeightedSquare { + template <typename Arithmetic> + constexpr decltype( auto ) operator()( const std::pair<Arithmetic, Arithmetic>& v ) const noexcept { + return v.first * v.first * v.second; + } + }; + + /** + * An Adder ValueHandler, taking weight into account and computing a count plus the sum of the weights + * In case of full atomicity, fetch_add or compare_exchange_weak are used for each element, + * that is we do not have full atomicity accross the two elements + */ + template <typename Arithmetic, atomicity Atomicity> + struct WeightedAdder { + using RegularType = std::pair<unsigned long, Arithmetic>; + using AtomicType = std::pair<std::atomic<unsigned long>, std::atomic<Arithmetic>>; + using OutputType = RegularType; + static constexpr bool isAtomic = Atomicity == atomicity::full; + using InternalType = std::conditional_t<isAtomic, AtomicType, OutputType>; + static constexpr OutputType getValue( const InternalType& v ) noexcept { + if constexpr ( isAtomic ) { + return { v.first.load( std::memory_order_relaxed ), v.second.load( std::memory_order_relaxed ) }; + } else { + return v; + } + }; + static RegularType exchange( InternalType& v, RegularType newv ) noexcept { + if constexpr ( isAtomic ) { + return { v.first.exchange( newv.first ), v.second.exchange( newv.second ) }; + } else { + return { std::exchange( v.first, newv.first ), std::exchange( v.second, newv.second ) }; + } + } + static constexpr OutputType DefaultValue() { return { 0, Arithmetic{} }; } + static void merge( InternalType& a, RegularType b ) noexcept { + if constexpr ( isAtomic ) { + fetch_add( a.first, b.first ); + fetch_add( a.second, b.second ); + } else { + a.first += b.first; + a.second += b.second; + } + }; + }; + + /** + * WeightedCountAccumulator. A WeightedCountAccumulator is an Accumulator storing the number of provided values, + * as well as the weighted version of it, aka. the sum of weights. It takes a pair (valueTuple, weight) as input + * @see Gaudi::Accumulators for detailed documentation + */ + template <atomicity Atomicity, typename Arithmetic> + struct WeightedCountAccumulator + : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, std::pair<unsigned long, Arithmetic>, Atomicity, Identity, + ExtractWeight, WeightedAdder<Arithmetic, Atomicity>> { + using Base = GenericAccumulator<std::pair<Arithmetic, Arithmetic>, std::pair<unsigned long, Arithmetic>, Atomicity, + Identity, ExtractWeight, WeightedAdder<Arithmetic, Atomicity>>; + using Base::Base; + using Base::operator+=; + /// overload of operator+= to be able to only give weight and no value + WeightedCountAccumulator operator+=( const Arithmetic weight ) { + *this += { 1ul, weight }; + return *this; + } + unsigned long nEntries() const { return this->rawValue().first; } + Arithmetic sumOfWeights() const { return this->rawValue().second; } + }; + + /** + * WeightedSumAccumulator. A WeightedSumAccumulator is an Accumulator storing a weighted sum of values. + * It takes a pair (valueTuple, weight) and basically sums the product of the last item othe 2 part of its in put pair + * : weight and value + * @see Gaudi::Accumulators for detailed documentation + */ + template <atomicity Atomicity, typename Arithmetic> + struct WeightedSumAccumulator + : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, WeightedProduct> { + using GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, + WeightedProduct>::GenericAccumulator; + Arithmetic sum() const { return this->value(); } + }; + + /** + * WeightedSquareAccumulator. A WeightedSquareAccumulator is an Accumulator storing a weighted sum of squared values. + * It basically takes a pair (value, weight) as input and sums weight*value*value + * @see Gaudi::Accumulators for detailed documentation + */ + template <atomicity Atomicity, typename Arithmetic = double> + struct WeightedSquareAccumulator + : GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, WeightedSquare> { + using GenericAccumulator<std::pair<Arithmetic, Arithmetic>, Arithmetic, Atomicity, + WeightedSquare>::GenericAccumulator; + Arithmetic sum2() const { return this->value(); }; + }; + + /** + * WeightedAveragingAccumulator. An WeightedAveragingAccumulator is an Accumulator able to compute an average + * This implementation takes a pair (value, weight) as input + * @see Gaudi::Accumulators for detailed documentation + */ + template <atomicity Atomicity, typename Arithmetic> + using WeightedAveragingAccumulator = + AveragingAccumulatorBase<Atomicity, Arithmetic, WeightedCountAccumulator, WeightedSumAccumulator>; + + /** + * WeightedSigmaAccumulator. A WeightedSigmaAccumulator is an Accumulator able to compute an average and variance/rms + * This implementation takes a pair (value, weight) as input + * @see Gaudi::Accumulators for detailed documentation + */ + template <atomicity Atomicity, typename Arithmetic> + using WeightedSigmaAccumulator = + SigmaAccumulatorBase<Atomicity, Arithmetic, WeightedAveragingAccumulator, WeightedSquareAccumulator>; + + /** + * Definition of a default type of Histogram Axis + * It contains number of bins, min and max value plus a title + * and defines the basic type of Axis (non log) + * It may also contain labels for the bins + */ + template <typename Arithmetic> + class Axis { + public: + using ArithmeticType = Arithmetic; + Axis( unsigned int nBins = 0, Arithmetic minValue = Arithmetic{}, Arithmetic maxValue = Arithmetic{}, + std::string title = {}, std::vector<std::string> labels = {} ) + : m_title( std::move( title ) ) + , nBins( nBins ) + , m_minValue( minValue ) + , m_maxValue( maxValue ) + , m_labels( std::move( labels ) ) { + details::requireValidTitle( m_title ); + recomputeRatio(); + for ( const auto& s : m_labels ) details::requireValidTitle( s ); + }; + explicit Axis( Gaudi::Histo1DDef const& def ) + : Axis( (unsigned int)def.bins(), def.lowEdge(), def.highEdge(), def.title() ){}; + + /// returns the bin number for a given value, ranging from 0 (underflow) to nBins+1 (overflow) + unsigned int index( Arithmetic value ) const { + // In case we use integer as Arithmetic type, we cannot use ratio for computing indices, + // as ratios < 1.0 will simply be 0, so we have to pay the division in such a case + int idx; + if constexpr ( std::is_integral_v<Arithmetic> ) { + idx = ( ( value - m_minValue ) * nBins / ( m_maxValue - m_minValue ) ) + 1; + } else { + idx = std::floor( ( value - m_minValue ) * m_ratio ) + 1; + } + return idx < 0 ? 0 : ( (unsigned int)idx > numBins() ? numBins() + 1 : (unsigned int)idx ); + } + + friend std::ostream& operator<<( std::ostream& o, Axis const& axis ) { + // Note that we print python code here, as the generic toStream implementation uses this + // operator to generate python code. + o << "(" << axis.numBins() << ", " << axis.minValue() << ", " << axis.maxValue() << ", " + << "\"" << axis.m_title << "\", ("; + for ( auto const& label : axis.m_labels ) { o << "\"" << label << "\", "; } + return o << "))"; + } + + /// says whether the given value is within the range of the axis + bool inAcceptance( Arithmetic value ) const { return value >= m_minValue && value <= m_maxValue; } + + // accessors + unsigned int numBins() const { return nBins; }; + void setNumBins( unsigned int n ) { + nBins = n; + recomputeRatio(); + } + Arithmetic minValue() const { return m_minValue; }; + void setMinValue( Arithmetic v ) { + m_minValue = v; + recomputeRatio(); + } + Arithmetic maxValue() const { return m_maxValue; }; + void setMaxValue( Arithmetic v ) { + m_maxValue = v; + recomputeRatio(); + } + std::string const& title() const { return m_title; } + void setTitle( std::string const& t ) { m_title = t; }; + std::vector<std::string> const labels() const { return m_labels; } + + private: + /// title of this axis + std::string m_title; + + public: + /// number of bins for this Axis + /// FIXME : should be private and called m_nBins but will break backward compatibility with previous implementation. + unsigned int nBins; + + private: + /// min and max values on this axis + Arithmetic m_minValue, m_maxValue; + /** + * precomputed ratio to convert a value into bin number + * equal to nBins/(maxValue-minValue). Only used for floating Arithmetic + */ + Arithmetic m_ratio; + /// labels for the bins + std::vector<std::string> m_labels; + + void recomputeRatio() { + m_ratio = ( m_maxValue == m_minValue ) ? Arithmetic{} : nBins / ( m_maxValue - m_minValue ); + } + }; + + /// automatic conversion of the Axis type to json + template <typename Arithmetic> + void to_json( nlohmann::json& j, const Axis<Arithmetic>& axis ) { + j = nlohmann::json{ { "nBins", axis.numBins() }, + { "minValue", axis.minValue() }, + { "maxValue", axis.maxValue() }, + { "title", axis.title() } }; + if ( !axis.labels().empty() ) { j["labels"] = axis.labels(); } + } + + /** + * small class used as InputType for regular Histograms + * basically a tuple of the given values, specialized in case of a single + * entry so that the syntax is more natural. + * NIndex should be lower than number of Arithmetic types and denotes the + * number of items used as index. There can typically be one more type in + * the list for profile histograms, not use as index on an axis + * ValueType is the actual type used to fill the histogram, that is + * the ArithmeticTuple reduced to NIndex items + * + * Note : the specialization is only needed to ensure backward compatibility + * with previous implementation where there was only one Arithmetic type common + * to all axis. It should in principal be the one and only implementation + * FIXME : remove specialization when client code was adapted + */ + template <typename Arithmetic, unsigned int NIndex> + struct HistoInputType : HistoInputType<make_tuple_t<Arithmetic, NIndex>, NIndex> { + using HistoInputType<make_tuple_t<Arithmetic, NIndex>, NIndex>::HistoInputType; + }; + template <unsigned int NIndex, typename... Elements> + struct HistoInputType<std::tuple<Elements...>, NIndex> : std::tuple<Elements...> { + using InternalType = std::tuple<Elements...>; + using ValueType = HistoInputType<SubTuple_t<InternalType, NIndex>, NIndex>; + using std::tuple<Elements...>::tuple; + template <class... AxisType, typename = typename std::enable_if_t<( sizeof...( AxisType ) == NIndex )>> + unsigned int computeIndex( std::tuple<AxisType...> const& axis ) const { + return computeIndexInternal<0, std::tuple<AxisType...>>( axis ); + } + template <class... AxisType, typename = typename std::enable_if_t<( sizeof...( AxisType ) == NIndex )>> + static unsigned int computeTotNBins( std::tuple<AxisType...> const& axis ) { + return computeTotNBinsInternal<0, std::tuple<AxisType...>>( axis ); + } + auto forInternalCounter() const { return 1ul; } + template <class... AxisType, typename = typename std::enable_if_t<( sizeof...( AxisType ) == NIndex )>> + bool inAcceptance( std::tuple<AxisType...> const& axis ) const { + return inAcceptanceInternal<0, std::tuple<AxisType...>>( axis ); + } + + private: + template <int N, class Tuple> + unsigned int computeIndexInternal( Tuple const& allAxis ) const { + // compute global index. Bins are stored in a column first manner + auto const& axis = std::get<N>( allAxis ); + unsigned int localIndex = axis.index( std::get<N>( *this ) ); + if constexpr ( N + 1 == NIndex ) + return localIndex; + else + return localIndex + ( axis.numBins() + 2 ) * computeIndexInternal<N + 1, Tuple>( allAxis ); + } + template <int N, class Tuple> + static unsigned int computeTotNBinsInternal( Tuple const& allAxis ) { + auto const& axis = std::get<N>( allAxis ); + unsigned int localNBins = axis.numBins() + 2; + if constexpr ( N + 1 == NIndex ) + return localNBins; + else + return localNBins * computeTotNBinsInternal<N + 1, Tuple>( allAxis ); + } + template <int N, class Tuple> + bool inAcceptanceInternal( Tuple const& allAxis ) const { + auto const& axis = std::get<N>( allAxis ); + bool localAnswer = axis.inAcceptance( std::get<N>( *this ) ); + if constexpr ( N + 1 == NIndex ) + return localAnswer; + else + return localAnswer || inAcceptanceInternal<N + 1, Tuple>( allAxis ); + } + }; + + /** + * small class used as InputType for weighted Histograms + * only a pair of the InnerType and the weight. + * See description of HistoInputType for more details + */ + template <typename ArithmeticTuple, unsigned int NIndex, typename WArithmetic> + struct WeightedHistoInputType : std::pair<HistoInputType<ArithmeticTuple, NIndex>, WArithmetic> { + using ValueType = typename HistoInputType<ArithmeticTuple, NIndex>::ValueType; + using std::pair<HistoInputType<ArithmeticTuple, NIndex>, WArithmetic>::pair; + template <class... AxisType, typename = typename std::enable_if_t<( sizeof...( AxisType ) == NIndex )>> + unsigned int computeIndex( std::tuple<AxisType...> const& axis ) const { + return this->first.computeIndex( axis ); + } + template <class... AxisType, typename = typename std::enable_if_t<( sizeof...( AxisType ) == NIndex )>> + static unsigned int computeTotNBins( std::tuple<AxisType...> const& axis ) { + return HistoInputType<ArithmeticTuple, NIndex>::computeTotNBins( axis ); + } + auto forInternalCounter() const { return std::pair( this->first.forInternalCounter(), this->second ); } + template <class... AxisType, typename = typename std::enable_if_t<( sizeof...( AxisType ) == NIndex )>> + bool inAcceptance( std::tuple<AxisType...> const& axis ) const { + return this->first.inAcceptance( axis ); + } + }; + + /** + * Internal Accumulator class dealing with Histograming. Templates parameters are : + * - Atomicity : none or full + * - InputType : a class holding a value given as input of the Histogram, + * and able to answer questions on this value given a number of axis matching + * the type of value. + * e.g. it would hold a pair of double for a non weighted 2D histogram or + * a pair of triplet of doubles and double for the weighted 3D histogram. + * Example implementations are (Weighted)HistoInputType. + * This class must define : + * + a constructor taking a set of value to build the InputType + * + a static method `unsigned int computeTotNBins( std::tuple<AxisType...> const& )` + * able to compute the total number of bins needed with this input type and + * these axis. It will typically be the product of the number of bins for each + * dimension, potentially increased by 2 for each if underflow and overflow + * is supported + * + a type ValueType alias defining the type of the input values to give to InputType + * This type needs to implement : + * * a method + * `unsigned int computeIndex( std::tuple<AxisType...> const& ) const` + * able to compute the bin corresponding to a given value + * * a method `auto forInternalCounter() const` returning the value to be used to + * inscrease the accumulator dealing with the bin associated with the current value. + * In most simple cases, it return `Arithmetic{}` or even `1` but for weighted + * histograms, it returns a pair with the weight as second item + * * in case of usage within a RootHistogram, it should also define a method + * `bool inAcceptance( std::tuple<AxisType...> const& )` checking whether a given + * value in within the range of the accumulator + * - Arithmetic : the arithmetic type used for values stored inside the histogram + * e.g. unsigned int for regular histogram as we only count entries, or float/double + * for weighted histograms, as we store actual sums of original values + * - BaseAccumulator : the underlying accumulator used in each bin + * - AxisTupleType : the types of the axis as a tuple. Its length defines the dimension + * of the Histogram this accumulator handles. + * The constraints on the AxisType are : FIXME use concepts when available + * + that they can be copied + * + that they have a ArithmeticType type alias + * + that they have a `unsigned int numBins() const` method + * + that they have a friend operator<< using std::ostream for printing + * + that they have a friend to_json method using nlohmann library + * + that they implement whatever is needed by the computeIndex and computeTotNBins methods + * of the InputType used. Plus the inAcceptance one if Roothistograms are used + * A default Axis class is provided for most cases + * This accumulator is simply an array of BaseAccumulator, one per bin. + */ + template <atomicity Atomicity, typename InputType, typename Arithmetic, + template <atomicity Ato, typename Arith> typename BaseAccumulatorT, typename AxisTupleType> + class HistogramingAccumulatorInternal { + template <atomicity, typename, typename, template <atomicity, typename> typename, typename> + friend class HistogramingAccumulatorInternal; + + public: + using ND = std::integral_constant<unsigned int, std::tuple_size_v<AxisTupleType>>; + using BaseAccumulator = BaseAccumulatorT<Atomicity, Arithmetic>; + using AxisTupleArithmeticType = typename InputType::ValueType; + HistogramingAccumulatorInternal( AxisTupleType axis ) + : m_axis{ axis } + , m_totNBins{ InputType::computeTotNBins( m_axis ) } + , m_value( new BaseAccumulator[m_totNBins] ) { + reset(); + } + template <atomicity ato> + HistogramingAccumulatorInternal( + construct_empty_t, + const HistogramingAccumulatorInternal<ato, InputType, Arithmetic, BaseAccumulatorT, AxisTupleType>& other ) + : m_axis( other.m_axis ), m_totNBins{ other.m_totNBins }, m_value( new BaseAccumulator[m_totNBins] ) { + reset(); + } + [[deprecated( "Use `++h1[x]`, `++h2[{x,y}]`, etc. instead." )]] HistogramingAccumulatorInternal& + operator+=( InputType v ) { + accumulator( v.computeIndex( m_axis ) ) += v.forInternalCounter(); + return *this; + } + void reset() { + for ( unsigned int index = 0; index < m_totNBins; index++ ) accumulator( index ).reset(); + } + template <atomicity ato> + void mergeAndReset( + HistogramingAccumulatorInternal<ato, InputType, Arithmetic, BaseAccumulatorT, AxisTupleType>& other ) { + assert( m_totNBins == other.m_totNBins ); + for ( unsigned int index = 0; index < m_totNBins; index++ ) { + accumulator( index ).mergeAndReset( other.accumulator( index ) ); + } + } + [[nodiscard]] auto operator[]( typename InputType::ValueType v ) { + return Buffer<BaseAccumulatorT, Atomicity, Arithmetic>{ accumulator( v.computeIndex( m_axis ) ) }; + } + + template <unsigned int N> + auto& axis() const { + return std::get<N>( m_axis ); + } + auto& axis() const { return m_axis; } + auto binValue( unsigned int i ) const { return accumulator( i ).value(); } + auto nEntries( unsigned int i ) const { return accumulator( i ).nEntries(); } + auto totNBins() const { return m_totNBins; } + + // FIXME These methods are there for backwrad compatibility with previous implementation + // where all Axis had to be of type Axis<...> and were stored in an array + // Newer code should call axis<N>().foo for whatever foo is defined in that axis type + auto nBins( unsigned int i ) const { return _getAxis( i, std::integral_constant<size_t, 0>() ).numBins(); } + auto minValue( unsigned int i ) const { return _getAxis( i, std::integral_constant<size_t, 0>() ).minValue(); } + auto maxValue( unsigned int i ) const { return _getAxis( i, std::integral_constant<size_t, 0>() ).maxValue(); } + + private: + BaseAccumulator& accumulator( unsigned int index ) const { + assert( index < m_totNBins ); + return m_value[index]; + } + + // FIXME Only used for backward compatibility. should be dropped at some stage + // Can only work if all axis have same type, which is no more the case + std::tuple_element_t<0, AxisTupleType> const& _getAxis( size_t i, + typename std::tuple_size<AxisTupleType>::type ) const { + throw std::logic_error( + fmt::format( "Retrieving axis {} in Histogram of dimension {}", i, std::tuple_size_v<AxisTupleType> ) ); + } + template <size_t N, typename = std::enable_if_t<std::tuple_size_v<AxisTupleType> != N>> + auto& _getAxis( size_t i, std::integral_constant<size_t, N> ) const { + if ( i == N ) return std::get<N>( m_axis ); + return _getAxis( i, std::integral_constant<size_t, N + 1>() ); + } + + /// set of Axis of this Histogram + AxisTupleType m_axis; + /// total number of bins in this histogram, under and overflow included + unsigned int m_totNBins; + /// Histogram content + std::unique_ptr<BaseAccumulator[]> m_value; + }; + + /** + * Class implementing a regular histogram accumulator + * + * Actually only an alias to HistogramingAccumulatorInternal with proper template parameters + */ + template <atomicity Atomicity, typename Arithmetic, typename ND, typename AxisTupleType> + using HistogramingAccumulator = + HistogramingAccumulatorInternal<Atomicity, HistoInputType<AxisToArithmetic_t<AxisTupleType>, ND::value>, + unsigned long, IntegralAccumulator, AxisTupleType>; + + /** + * Class implementing a weighted histogram accumulator + * + * Actually only an alias to HistogramingAccumulatorInternal with proper template parameters + */ + template <atomicity Atomicity, typename Arithmetic, typename ND, typename AxisTupleType> + using WeightedHistogramingAccumulator = + HistogramingAccumulatorInternal<Atomicity, + WeightedHistoInputType<AxisToArithmetic_t<AxisTupleType>, ND::value, Arithmetic>, + Arithmetic, WeightedCountAccumulator, AxisTupleType>; + + /** + * Class implementing a profile histogram accumulator + * + * Actually only an alias to HistogramingAccumulatorInternal with proper template parameters + */ + template <atomicity Atomicity, typename Arithmetic, typename ND, typename AxisTupleType> + using ProfileHistogramingAccumulator = + HistogramingAccumulatorInternal<Atomicity, + HistoInputType<ProfileAxisToArithmetic_t<Arithmetic, AxisTupleType>, ND::value>, + Arithmetic, SigmaAccumulator, AxisTupleType>; + + /** + * Class implementing a weighted profile histogram accumulator + * + * Actually only an alias to HistogramingAccumulatorInternal with proper template parameters + */ + template <atomicity Atomicity, typename Arithmetic, typename ND, typename AxisTupleType> + using WeightedProfileHistogramingAccumulator = HistogramingAccumulatorInternal< + Atomicity, WeightedHistoInputType<ProfileAxisToArithmetic_t<Arithmetic, AxisTupleType>, ND::value, Arithmetic>, + Arithmetic, WeightedSigmaAccumulator, AxisTupleType>; + + /** + * A base counter dealing with Histograms + * + * Main features of that Counter : + * - can be any number of dimensions. The dimension is its first template parameter + * - for each dimension, an Axis is associated. Axis can be of any type depending + * on the underlying accumulator + * - the constructor expects one extra argument per axis, typically a tuple + * of values allowing to create the Axis objects in the back + * - the operator+= takes either an array of values (one per dimension) + * or a tuple<array of values, weight>. The value inside the bin + * corresponding to the given values is then increased by 1 or weight + * - the prefered syntax is to avoid operator+= and use operator[] to get a + * buffer on the bin you're updating. Syntax becomes : + * ++counter[{x,y}] or wcounter[{x,y}] += w + * - the Counter is templated by the types of the values given to + * operator+ and also by the type stored into the bins + * - the counter can be atomic or not and supports buffering. Note that + * the atomicity is classical eventual consistency. So each bin is + * atomically updated but bins are not garanted to be coherent when + * reading all of them back + * - profile histograms are also supported, operator+= takes one more + * value in the array of values in that case + * + * This base class is then aliases for the 4 standard cases of Histogram, + * WeightedHistogram, ProfileHistogram and WeightedProfileHistogram. + * For all predefined Histogram types, the axis type is a simple triplet + * of values nbins, minValue, maxValue plus a title. + * + * Typical usage : + * \code + * Histogram<2, double, atomicity::full> + * counter{owner, "CounterName", "HistoTitle", {{nBins1, minVal1, maxVal1}, {nBins2, minVal2, maxVal2, + * "AxisTitle"}}}; + * ++counter[{val1, val2}]; // prefered syntax + * counter += {val1, val2}; // original syntax inherited from counters + * + * WeightedHistogram<2, double, atomicity::full> + * wcounter{owner, "CounterName", "HistoTitle", {{nBins1, minVal1, maxVal1}, {nBins2, minVal2, maxVal2, + * "AxisTitle"}}}; wcounter[{val1, val2}] += w; // prefered syntax wcounter += {{val1, val2}, w}; // original + * syntax inherited from counters \endcode + * + * When serialized to json, this counter uses new types histogram:Histogram:<prec>, histogram:ProfileHistogram:<prec>, + * histogram:WeightedHistogram:<prec> and histrogram:WeightedProfileHistogram:<prec> + * <prec> id described in Accumulators.h and decribes here the precision of the bin data. + * All these types have the same fields, namely : + * dimension(integer), title(string), empty(bool), nEntries(integer), axis(array), bins(array) + * where : + * + axis is an array of tuples, one per dimension, with content (nBins(integer), minValue(number), + * maxValue(number), title(string)) for the default type of Axis + * + bins is an array of values + * - The length of the array is the product of (nBins+2) for all axis + * - the +2 is because the bin 0 is the one for values below minValue and bin nBins+1 is the one for values + * above maxValue bins are stored row first so we iterate first on highest dimension + * - the value is a number for non profile histograms + * - the value is of the form ( (nEntries(integer), sum(number) ), sum2(number) ) for profile histograms + * Note the pair with a pair as first entry + */ + template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type, + template <atomicity, typename, typename, typename> typename Accumulator, typename AxisTupleType> + class HistogramingCounterBase; + template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type, + template <atomicity, typename, typename, typename> typename Accumulator, typename... AxisTypes> + class HistogramingCounterBase<ND, Atomicity, Arithmetic, Type, Accumulator, std::tuple<AxisTypes...>> + : public BufferableCounter<Atomicity, Accumulator, Arithmetic, std::integral_constant<unsigned int, ND>, + std::tuple<AxisTypes...>> { + public: + using AxisTupleType = std::tuple<AxisTypes...>; + using NumberDimensions = std::integral_constant<unsigned int, ND>; + using Parent = BufferableCounter<Atomicity, Accumulator, Arithmetic, NumberDimensions, AxisTupleType>; + using AccumulatorType = Accumulator<Atomicity, Arithmetic, NumberDimensions, AxisTupleType>; + using AxisTupleArithmeticType = typename AccumulatorType::AxisTupleArithmeticType; + /// for backward compatibility with previous implementation, should not be used FIXME + using AxisArithmeticType = typename std::tuple_element<0, AxisTupleType>::type::ArithmeticType; + inline static const std::string typeString{ std::string{ Type } + ':' + typeid( Arithmetic ).name() }; + /// This constructor takes the axis as a tuple + template <typename OWNER> + HistogramingCounterBase( OWNER* owner, std::string const& name, std::string const& title, AxisTupleType axis ) + : Parent( owner, name, *this, axis ), m_title( title ) { + details::requireValidTitle( m_title ); + } + /// This constructor takes the axis one by one, when ND >= 2. If ND = 1, the other one can be used + template <typename OWNER> + HistogramingCounterBase( OWNER* owner, std::string const& name, std::string const& title, AxisTypes... allAxis ) + : HistogramingCounterBase( owner, name, title, std::make_tuple( allAxis... ) ) {} + using Parent::print; + template <typename stream> + stream& printImpl( stream& o, bool /*tableFormat*/ ) const { + o << ND << "D Histogram with config "; + std::apply( [&o]( auto&&... args ) { ( ( o << args << "\n" ), ... ); }, this->axis() ); + return o; + } + std::ostream& print( std::ostream& o, bool tableFormat = false ) const override { + return printImpl( o, tableFormat ); + } + MsgStream& print( MsgStream& o, bool tableFormat = false ) const override { return printImpl( o, tableFormat ); } + friend void reset( HistogramingCounterBase& c ) { c.reset(); } + friend void mergeAndReset( HistogramingCounterBase& h, HistogramingCounterBase& o ) { h.mergeAndReset( o ); } + friend void to_json( nlohmann::json& j, HistogramingCounterBase const& h ) { h.to_json( j ); } + virtual void to_json( nlohmann::json& j ) const { + // get all bin values and compute total nbEntries + std::vector<typename AccumulatorType::BaseAccumulator::OutputType> bins; + bins.reserve( this->totNBins() ); + unsigned long totNEntries{ 0 }; + for ( unsigned int i = 0; i < this->totNBins(); i++ ) { + bins.push_back( this->binValue( i ) ); + totNEntries += this->nEntries( i ); + } + // build json + j = { { "type", std::string( Type ) + ":" + typeid( Arithmetic ).name() }, + { "title", m_title }, + { "dimension", ND }, + { "empty", totNEntries == 0 }, + { "nEntries", totNEntries }, + { "axis", this->axis() }, + { "bins", bins } }; + } + std::string const& title() const { return m_title; } + + protected: + std::string const m_title; + }; + + namespace naming { + inline constexpr char histogramString[] = "histogram:Histogram"; + inline constexpr char weightedHistogramString[] = "histogram:WeightedHistogram"; + inline constexpr char profilehistogramString[] = "histogram:ProfileHistogram"; + inline constexpr char weightedProfilehistogramString[] = "histogram:WeightedProfileHistogram"; + } // namespace naming + + /// standard static histograming counter. See HistogramingCounterBase for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using StaticHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::histogramString, + HistogramingAccumulator, AxisTupleType>; + + /// standard static histograming counter with weight. See HistogramingCounterBase for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using StaticWeightedHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::weightedHistogramString, + WeightedHistogramingAccumulator, AxisTupleType>; + + /// profile static histograming counter. See HistogramingCounterBase for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using StaticProfileHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::profilehistogramString, + ProfileHistogramingAccumulator, AxisTupleType>; + + /// weighted static profile histograming counter. See HistogramingCounterBase for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using StaticWeightedProfileHistogram = + HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::weightedProfilehistogramString, + WeightedProfileHistogramingAccumulator, AxisTupleType>; + +} // namespace Gaudi::Accumulators diff --git a/GaudiKernel/include/Gaudi/Accumulators/StaticRootHistogram.h b/GaudiKernel/include/Gaudi/Accumulators/StaticRootHistogram.h new file mode 100644 index 0000000000000000000000000000000000000000..a6d14050de71182f5ab05a252684481924565eb9 --- /dev/null +++ b/GaudiKernel/include/Gaudi/Accumulators/StaticRootHistogram.h @@ -0,0 +1,361 @@ +/***********************************************************************************\ +* (c) Copyright 1998-2024 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/StaticHistogram.h> + +#include <type_traits> + +namespace { + template <typename tuple_t> + constexpr auto get_array_from_tuple( tuple_t&& tuple ) { + constexpr auto get_array = []( auto&&... x ) { return std::array{ std::forward<decltype( x )>( x )... }; }; + return std::apply( get_array, std::forward<tuple_t>( tuple ) ); + } +} // namespace + +namespace Gaudi::Accumulators { + + /// number of items in sums for a given dimension + /// = 1 (nb items) + ND (sums of each dimension) + ND*(ND+1)/2 (square sums) + 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 )> { + /// copy constructor from non thread safe type + 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(); } + }; + + /// specialization for ND=1 to allow for better syntax + 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(); } + }; + + /** + * Internal Accumulator class dealing with RootHistograming. + * Actually a simple extention on top of RootHistograming with an + * extra SigmaCounter embeded + */ + template <atomicity Atomicity, typename Arithmetic, unsigned int ND, typename AxisTupleType> + class RootHistogramingAccumulatorInternal + : public HistogramingAccumulatorInternal<Atomicity, HistoInputType<AxisToArithmetic_t<AxisTupleType>, ND>, + unsigned long, IntegralAccumulator, AxisTupleType> { + + using InputType = HistoInputType<AxisToArithmetic_t<AxisTupleType>, ND>; + using Parent = + HistogramingAccumulatorInternal<Atomicity, InputType, unsigned long, IntegralAccumulator, AxisTupleType>; + + template <atomicity, typename, unsigned int, typename> + friend class RootHistogramingAccumulatorInternal; + + static_assert( ND <= 3, "Root on supports histogrmas with dimension <= 3" ); + + /** + * Small procyclass allowing operator[] to work as expected on the RootHistogram + * that is to return something having an operator+= updating the histogram properly + */ + 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+=( typename InputType::ValueType v ) { + update( v ); + return *this; + } + void reset() { + m_accumulator.reset(); + Parent::reset(); + } + template <atomicity ato> + void mergeAndReset( RootHistogramingAccumulatorInternal<ato, Arithmetic, ND, AxisTupleType>& other ) { + m_accumulator.mergeAndReset( other.m_accumulator ); + Parent::mergeAndReset( other ); + } + [[nodiscard]] auto operator[]( typename InputType::ValueType v ) { return Proxy( *this, v ); } + + /// returns the nbentries, sums and "squared sums" of the inputs + /// Practically we have first the number of entries, then the simple sums of each + /// input dimension followed by all combinasions of product of 2 inputs, in "alphabetical" order, + /// e.g. for ND=3 we have sums of n, x, y, z, x^2, xy, xz, y^2, yz, z^2 + 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() ) ) { + if constexpr ( ND == 1 ) { + m_accumulator += std::get<0>( v ); + } else { + m_accumulator += get_array_from_tuple( 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 { + /// helper function to compute standard_deviation + 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 + + /** + * Class implementing a root histogram accumulator + */ + template <atomicity Atomicity, typename Arithmetic, typename ND, typename AxisTupleType> + struct RootHistogramingAccumulator; + + template <atomicity Atomicity, typename Arithmetic, typename AxisTupleType> + struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<unsigned int, 1>, AxisTupleType> + : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1, AxisTupleType> { + using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1, + AxisTupleType>::RootHistogramingAccumulatorInternal; + using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1, AxisTupleType>::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, typename AxisTupleType> + struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<unsigned int, 2>, AxisTupleType> + : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2, AxisTupleType> { + using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2, + AxisTupleType>::RootHistogramingAccumulatorInternal; + using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2, AxisTupleType>::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, typename AxisTupleType> + struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<unsigned int, 3>, AxisTupleType> + : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3, AxisTupleType> { + using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3, + AxisTupleType>::RootHistogramingAccumulatorInternal; + using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3, AxisTupleType>::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() ); } + }; + + /** + * Extension of the standard Gaudi histogram to provide similar functionnality as ROOT + * + * The main piece of extra functionnality is that ablity to compute a mean and sigma + * of the data set hold by the histogram on top of the histogram itself + * This will then be used when filling a real Root histogram so that one gets the + * same values as if working directly with Root. + * When using pure Gaudi Histograms, and converting to ROOT, ROOT automatically + * recomputed average and sigma form the bins, but the value is not the expected one + * + * Usage is similar to HistogramingCounterBase so see the documentation there + * Serialization has the following extra fields : + * - nTotEntries, sum, sum2, mean in 1D + * - nTotEntries, sumx, sumy, sumx2, sumy2, sumxy, meanx, meany in 2D + * - nTotEntries, sumx, sumy, sumz, sumx2, sumy2, sumz2, sumxy, sumxz, sumyz, meanx, meany, meanz in 3D + * and uses same types as HistogramingCounterBase + * + */ + template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type, typename AxisTupleType> + class RootHistogramingCounterBase; + + template <atomicity Atomicity, typename Arithmetic, const char* Type, typename AxisTupleType> + class RootHistogramingCounterBase<1, Atomicity, Arithmetic, Type, AxisTupleType> + : public HistogramingCounterBase<1, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType> { + public: + using Parent = HistogramingCounterBase<1, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType>; + using Parent::Parent; + + friend void to_json( nlohmann::json& j, + RootHistogramingCounterBase<1, Atomicity, Arithmetic, Type, AxisTupleType> 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, typename AxisTupleType> + class RootHistogramingCounterBase<2, Atomicity, Arithmetic, Type, AxisTupleType> + : public HistogramingCounterBase<2, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType> { + public: + using Parent = HistogramingCounterBase<2, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType>; + using Parent::Parent; + + friend void to_json( nlohmann::json& j, + RootHistogramingCounterBase<2, Atomicity, Arithmetic, Type, AxisTupleType> 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, typename AxisTupleType> + class RootHistogramingCounterBase<3, Atomicity, Arithmetic, Type, AxisTupleType> + : public HistogramingCounterBase<3, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType> { + public: + using Parent = HistogramingCounterBase<3, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType>; + using Parent::Parent; + + friend void to_json( nlohmann::json& j, + RootHistogramingCounterBase<3, Atomicity, Arithmetic, Type, AxisTupleType> 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(); + } + }; + + /// Root histograming counter. See RootHistogramingCounterBase for details + template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double, + typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>> + using StaticRootHistogram = + RootHistogramingCounterBase<ND, Atomicity, Arithmetic, naming::histogramString, AxisTupleType>; + +} // namespace Gaudi::Accumulators diff --git a/GaudiKernel/include/Gaudi/Histograming/Sink/Utils.h b/GaudiKernel/include/Gaudi/Histograming/Sink/Utils.h index 3cf02ccb85574d43f3a88170a70e3e53f1c77be7..f5d4a0d7bd6191054309e4780b02ab28fb16a959 100644 --- a/GaudiKernel/include/Gaudi/Histograming/Sink/Utils.h +++ b/GaudiKernel/include/Gaudi/Histograming/Sink/Utils.h @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2023 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -10,7 +10,7 @@ \***********************************************************************************/ #pragma once -#include <Gaudi/Accumulators/Histogram.h> +#include <Gaudi/Accumulators/StaticHistogram.h> #include <Gaudi/MonitoringHub.h> #include <GaudiKernel/GaudiException.h> #include <TDirectory.h> @@ -411,23 +411,23 @@ namespace Gaudi::Histograming::Sink { saveRootHisto<Traits<isProfile, ROOTHisto, N>>( file, std::move( dir ), std::move( name ), j ); } - /// Direct conversion of 1D histograms form Gaudi to ROOT format + /// Direct conversion of 1D histograms from Gaudi to ROOT format template <Accumulators::atomicity Atomicity = Accumulators::atomicity::full, typename Arithmetic = double> details::ProfileWrapper<TProfile> profileHisto1DToRoot( std::string name, Monitoring::Hub::Entity const& ent ) { // get original histogram from the Entity auto const& gaudiHisto = - *reinterpret_cast<Gaudi::Accumulators::ProfileHistogram<1, Atomicity, Arithmetic>*>( ent.id() ); + *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<1, Atomicity, Arithmetic>*>( ent.id() ); // convert to Root - auto const& gaudiAxis = gaudiHisto.axis()[0]; - details::ProfileWrapper<TProfile> histo{ name.c_str(), gaudiHisto.title().c_str(), gaudiAxis.nBins, - gaudiAxis.minValue, gaudiAxis.maxValue }; + auto const& gaudiAxis = gaudiHisto.template axis<0>(); + details::ProfileWrapper<TProfile> histo{ name.c_str(), gaudiHisto.title().c_str(), gaudiAxis.numBins(), + gaudiAxis.minValue(), gaudiAxis.maxValue() }; histo.Sumw2(); unsigned long totNEntries{ 0 }; for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) { totNEntries += gaudiHisto.nEntries( i ); } - if ( !gaudiAxis.labels.empty() ) { + if ( !gaudiAxis.labels().empty() ) { auto& axis = *histo.GetXaxis(); - for ( unsigned int i = 0; i < gaudiAxis.labels.size(); i++ ) { - axis.SetBinLabel( i + 1, gaudiAxis.labels[i].c_str() ); + for ( unsigned int i = 0; i < gaudiAxis.labels().size(); i++ ) { + axis.SetBinLabel( i + 1, gaudiAxis.labels()[i].c_str() ); } } for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) { @@ -441,29 +441,29 @@ namespace Gaudi::Histograming::Sink { return histo; } - /// Direct conversion of 2D histograms form Gaudi to ROOT format + /// Direct conversion of 2D histograms from Gaudi to ROOT format template <Accumulators::atomicity Atomicity = Accumulators::atomicity::full, typename Arithmetic = double> details::ProfileWrapper<TProfile2D> profileHisto2DToRoot( std::string name, Monitoring::Hub::Entity const& ent ) { // get original histogram from the Entity auto const& gaudiHisto = - *reinterpret_cast<Gaudi::Accumulators::ProfileHistogram<2, Atomicity, Arithmetic>*>( ent.id() ); + *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<2, Atomicity, Arithmetic>*>( ent.id() ); // convert to Root - auto const& gaudiXAxis = gaudiHisto.axis()[0]; - auto const& gaudiYAxis = gaudiHisto.axis()[1]; - details::ProfileWrapper<TProfile2D> histo{ name.c_str(), gaudiHisto.title().c_str(), gaudiXAxis.nBins, - gaudiXAxis.minValue, gaudiXAxis.maxValue, gaudiYAxis.nBins, - gaudiYAxis.minValue, gaudiYAxis.maxValue }; + auto const& gaudiXAxis = gaudiHisto.template axis<0>(); + auto const& gaudiYAxis = gaudiHisto.template axis<1>(); + details::ProfileWrapper<TProfile2D> histo{ name.c_str(), gaudiHisto.title().c_str(), gaudiXAxis.numBins(), + gaudiXAxis.minValue(), gaudiXAxis.maxValue(), gaudiYAxis.numBins(), + gaudiYAxis.minValue(), gaudiYAxis.maxValue() }; histo.Sumw2(); - if ( !gaudiXAxis.labels.empty() ) { + if ( !gaudiXAxis.labels().empty() ) { auto& axis = *histo.GetXaxis(); - for ( unsigned int i = 0; i < gaudiXAxis.labels.size(); i++ ) { - axis.SetBinLabel( i + 1, gaudiXAxis.labels[i].c_str() ); + for ( unsigned int i = 0; i < gaudiXAxis.labels().size(); i++ ) { + axis.SetBinLabel( i + 1, gaudiXAxis.labels()[i].c_str() ); } } - if ( !gaudiYAxis.labels.empty() ) { + if ( !gaudiYAxis.labels().empty() ) { auto& axis = *histo.GetYaxis(); - for ( unsigned int i = 0; i < gaudiYAxis.labels.size(); i++ ) { - axis.SetBinLabel( i + 1, gaudiYAxis.labels[i].c_str() ); + for ( unsigned int i = 0; i < gaudiYAxis.labels().size(); i++ ) { + axis.SetBinLabel( i + 1, gaudiYAxis.labels()[i].c_str() ); } } for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) { @@ -479,37 +479,37 @@ namespace Gaudi::Histograming::Sink { return histo; } - /// Direct conversion of 3D histograms form Gaudi to ROOT format + /// Direct conversion of 3D histograms from Gaudi to ROOT format template <Accumulators::atomicity Atomicity = Accumulators::atomicity::full, typename Arithmetic = double> details::ProfileWrapper<TProfile3D> profileHisto3DToRoot( std::string name, Monitoring::Hub::Entity const& ent ) { // get original histogram from the Entity auto const& gaudiHisto = - *reinterpret_cast<Gaudi::Accumulators::ProfileHistogram<3, Atomicity, Arithmetic>*>( ent.id() ); + *reinterpret_cast<Gaudi::Accumulators::StaticProfileHistogram<3, Atomicity, Arithmetic>*>( ent.id() ); // convert to Root - auto const& gaudiXAxis = gaudiHisto.axis()[0]; - auto const& gaudiYAxis = gaudiHisto.axis()[1]; - auto const& gaudiZAxis = gaudiHisto.axis()[2]; - details::ProfileWrapper<TProfile3D> histo{ name.c_str(), gaudiHisto.title().c_str(), gaudiXAxis.nBins, - gaudiXAxis.minValue, gaudiXAxis.maxValue, gaudiYAxis.nBins, - gaudiYAxis.minValue, gaudiYAxis.maxValue, gaudiZAxis.nBins, - gaudiZAxis.minValue, gaudiZAxis.maxValue }; + auto const& gaudiXAxis = gaudiHisto.template axis<0>(); + auto const& gaudiYAxis = gaudiHisto.template axis<1>(); + auto const& gaudiZAxis = gaudiHisto.template axis<2>(); + details::ProfileWrapper<TProfile3D> histo{ name.c_str(), gaudiHisto.title().c_str(), gaudiXAxis.numBins(), + gaudiXAxis.minValue(), gaudiXAxis.maxValue(), gaudiYAxis.numBins(), + gaudiYAxis.minValue(), gaudiYAxis.maxValue(), gaudiZAxis.numBins(), + gaudiZAxis.minValue(), gaudiZAxis.maxValue() }; histo.Sumw2(); - if ( !gaudiXAxis.labels.empty() ) { + if ( !gaudiXAxis.labels().empty() ) { auto& axis = *histo.GetXaxis(); - for ( unsigned int i = 0; i < gaudiXAxis.labels.size(); i++ ) { - axis.SetBinLabel( i + 1, gaudiXAxis.labels[i].c_str() ); + for ( unsigned int i = 0; i < gaudiXAxis.labels().size(); i++ ) { + axis.SetBinLabel( i + 1, gaudiXAxis.labels()[i].c_str() ); } } - if ( !gaudiYAxis.labels.empty() ) { + if ( !gaudiYAxis.labels().empty() ) { auto& axis = *histo.GetYaxis(); - for ( unsigned int i = 0; i < gaudiYAxis.labels.size(); i++ ) { - axis.SetBinLabel( i + 1, gaudiYAxis.labels[i].c_str() ); + for ( unsigned int i = 0; i < gaudiYAxis.labels().size(); i++ ) { + axis.SetBinLabel( i + 1, gaudiYAxis.labels()[i].c_str() ); } } - if ( !gaudiZAxis.labels.empty() ) { + if ( !gaudiZAxis.labels().empty() ) { auto& axis = *histo.GetZaxis(); - for ( unsigned int i = 0; i < gaudiZAxis.labels.size(); i++ ) { - axis.SetBinLabel( i + 1, gaudiZAxis.labels[i].c_str() ); + for ( unsigned int i = 0; i < gaudiZAxis.labels().size(); i++ ) { + axis.SetBinLabel( i + 1, gaudiZAxis.labels()[i].c_str() ); } } for ( unsigned int i = 0; i < gaudiHisto.totNBins(); i++ ) { diff --git a/GaudiKernel/tests/src/CounterHistoArrayUnitTest.cpp b/GaudiKernel/tests/src/CounterHistoArrayUnitTest.cpp index 33afda6c12baec995721510a5c1eb68f2a0e7f31..d82696803ae329a01600846b6f66a3796e9bd4b3 100644 --- a/GaudiKernel/tests/src/CounterHistoArrayUnitTest.cpp +++ b/GaudiKernel/tests/src/CounterHistoArrayUnitTest.cpp @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2021 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -10,7 +10,12 @@ \***********************************************************************************/ #define BOOST_TEST_DYN_LINK #define BOOST_TEST_MODULE test_CounterHistos + +#include <Gaudi/Accumulators/Histogram.h> #include <Gaudi/Accumulators/HistogramArray.h> +#include <Gaudi/Accumulators/StaticHistogram.h> + +#include <GaudiKernel/PropertyHolder.h> #include <boost/test/unit_test.hpp> @@ -24,10 +29,21 @@ namespace { MonitoringHub& monitoringHub() { return m_monitHub; } MonitoringHub m_monitHub{}; }; - struct Algo { - ServiceLocator* serviceLocator() { return &m_serviceLocator; } - std::string name() { return ""; } - ServiceLocator m_serviceLocator{}; + struct BaseAlgo : INamedInterface, IProperty { + unsigned long addRef() override { return 0; }; + unsigned long release() override { return 0; }; + void* i_cast( const InterfaceID& ) const override { return nullptr; } + std::vector<std::string> getInterfaceNames() const override { return {}; } + unsigned long refCount() const override { return 0; } + StatusCode queryInterface( const InterfaceID&, void** ) override { return StatusCode::FAILURE; }; + const std::string& name() const override { return m_name; }; + std::string m_name{}; + }; + struct Algo : PropertyHolder<BaseAlgo> { + ServiceLocator* serviceLocator() { return &m_serviceLocator; } + ServiceLocator m_serviceLocator{}; + void registerCallBack( Gaudi::StateMachine::Transition, std::function<void()> ) {} + Gaudi::StateMachine::State FSMState() const { return Gaudi::StateMachine::CONFIGURED; } }; // Little helper for using automatic nlohmann conversion mechanism @@ -38,13 +54,60 @@ namespace { } } // namespace +BOOST_AUTO_TEST_CASE( test_static_counter_histos, *boost::unit_test::tolerance( 1e-14 ) ) { + Algo algo; + + { + // testing an array of 5 1D, regular histograms, with "standard" names / titles + Gaudi::Accumulators::HistogramArray<Gaudi::Accumulators::StaticHistogram<1>, 5> histo1d{ + &algo, "SGaudiH1D-{}", "A Gaudi 1D histogram - number {}", { 21, -10.5, 10.5, "X" } }; + for ( unsigned int i = 0; i < 5; i++ ) ++histo1d[i][-10.0]; // fill the first (non-overflow) bin + for ( unsigned int i = 0; i < 5; i++ ) BOOST_TEST( toJSON( histo1d[i] ).at( "bins" )[1] == 1 ); + BOOST_TEST( toJSON( histo1d[2] ).at( "title" ) == "A Gaudi 1D histogram - number 2" ); + ++histo1d[3][-10.0]; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo1d[3] ).at( "bins" )[1] == 2 ); + } + + { + // testing an array of 7 2D, weighted histograms, with "standard" names / titles + Gaudi::Accumulators::HistogramArray<Gaudi::Accumulators::StaticWeightedHistogram<2>, 7> histo2dw{ + &algo, "SName{}", "Title {}", { 21, -10.5, 10.5, "X" }, { 21, -10.5, 10.5, "Y" } }; + for ( unsigned int i = 0; i < 7; i++ ) histo2dw[i][{ -10.0, -10.0 }] += 0.25; // fill the first (non-overflow) bin + for ( unsigned int i = 0; i < 7; i++ ) BOOST_TEST( toJSON( histo2dw[i] ).at( "bins" )[( 1 + 21 + 1 ) + 1] == 0.25 ); + for ( unsigned int i = 0; i < 7; i++ ) histo2dw[i][{ -10.0, -10.0 }] += 0.5; // fill the first (non-overflow) bin + for ( unsigned int i = 0; i < 7; i++ ) BOOST_TEST( toJSON( histo2dw[i] ).at( "bins" )[( 1 + 21 + 1 ) + 1] == 0.75 ); + } + + { + Gaudi::Accumulators::HistogramArray<Gaudi::Accumulators::StaticHistogram<1>, 5> histo1d{ + &algo, + []( int n ) { return fmt::format( "SGaudiH1D-{}-{}", n, n ^ 2 ); }, + [nb = 5]( int n ) { + return fmt::format( "Title number {} of Histogram arrays of {} histograms in total", n, nb ); + }, + { 21, -10.5, 10.5, "X" } }; + for ( unsigned int i = 0; i < 5; i++ ) ++histo1d[i][-10.0]; // fill the first (non-overflow) bin + for ( unsigned int i = 0; i < 5; i++ ) BOOST_TEST( toJSON( histo1d[i] ).at( "bins" )[1] == 1 ); + BOOST_TEST( toJSON( histo1d[3] ).at( "title" ) == "Title number 3 of Histogram arrays of 5 histograms in total" ); + ++histo1d[3][-10.0]; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo1d[3] ).at( "bins" )[1] == 2 ); + } +} + BOOST_AUTO_TEST_CASE( test_counter_histos, *boost::unit_test::tolerance( 1e-14 ) ) { Algo algo; { // testing an array of 5 1D, regular histograms, with "standard" names / titles + // Gaudi::Accumulators::HistogramArray<Gaudi::Accumulators::Histogram<1>, 5> histo1d{ Gaudi::Accumulators::HistogramArray<Gaudi::Accumulators::Histogram<1>, 5> histo1d{ - &algo, "GaudiH1D-{}", "A Gaudi 1D histogram - number {}", { 21, -10.5, 10.5, "X" } }; + &algo, "GaudiH1D-{}", "WRONG TITLE !", { 1, -1, 1, "WRONG" } }; + for ( unsigned int i = 0; i < 5; i++ ) { + algo.setProperty( fmt::format( "GaudiH1D_{}_Title", i ), fmt::format( "A Gaudi 1D histogram - number {}", i ) ) + .ignore(); + algo.setProperty( fmt::format( "GaudiH1D_{}_Axis0", i ), "( 21, -10.5, 10.5, \"X\" )" ).ignore(); + } + for ( unsigned int i = 0; i < 5; i++ ) { histo1d[i].createHistogram( algo ); } for ( unsigned int i = 0; i < 5; i++ ) ++histo1d[i][-10.0]; // fill the first (non-overflow) bin for ( unsigned int i = 0; i < 5; i++ ) BOOST_TEST( toJSON( histo1d[i] ).at( "bins" )[1] == 1 ); BOOST_TEST( toJSON( histo1d[2] ).at( "title" ) == "A Gaudi 1D histogram - number 2" ); @@ -55,7 +118,13 @@ BOOST_AUTO_TEST_CASE( test_counter_histos, *boost::unit_test::tolerance( 1e-14 ) { // testing an array of 7 2D, weighted histograms, with "standard" names / titles Gaudi::Accumulators::HistogramArray<Gaudi::Accumulators::WeightedHistogram<2>, 7> histo2dw{ - &algo, "Name{}", "Title {}", { 21, -10.5, 10.5, "X" }, { 21, -10.5, 10.5, "Y" } }; + &algo, "Name{}", "WRONG TITLE !", { 1, -1, 1, "WRONG" }, { 1, -1, 1, "WRONG" } }; + for ( unsigned int i = 0; i < 7; i++ ) { + algo.setProperty( fmt::format( "Name{}_Title", i ), fmt::format( "Title {}", i ) ).ignore(); + algo.setProperty( fmt::format( "Name{}_Axis0", i ), "( 21, -10.5, 10.5, \"X\" )" ).ignore(); + algo.setProperty( fmt::format( "Name{}_Axis1", i ), "( 21, -10.5, 10.5, \"Y\" )" ).ignore(); + } + for ( unsigned int i = 0; i < 7; i++ ) { histo2dw[i].createHistogram( algo ); } for ( unsigned int i = 0; i < 7; i++ ) histo2dw[i][{ -10.0, -10.0 }] += 0.25; // fill the first (non-overflow) bin for ( unsigned int i = 0; i < 7; i++ ) BOOST_TEST( toJSON( histo2dw[i] ).at( "bins" )[( 1 + 21 + 1 ) + 1] == 0.25 ); for ( unsigned int i = 0; i < 7; i++ ) histo2dw[i][{ -10.0, -10.0 }] += 0.5; // fill the first (non-overflow) bin @@ -65,11 +134,15 @@ BOOST_AUTO_TEST_CASE( test_counter_histos, *boost::unit_test::tolerance( 1e-14 ) { Gaudi::Accumulators::HistogramArray<Gaudi::Accumulators::Histogram<1>, 5> histo1d{ &algo, - []( int n ) { return fmt::format( "GaudiH1D-{}-{}", n, n ^ 2 ); }, + []( int n ) { return fmt::format( "GaudiH1D-{}-{}", n, n * n ); }, [nb = 5]( int n ) { return fmt::format( "Title number {} of Histogram arrays of {} histograms in total", n, nb ); }, - { 21, -10.5, 10.5, "X" } }; + { 1, -1, 1, "WRONG" } }; + for ( unsigned int i = 0; i < 5; i++ ) { + algo.setProperty( fmt::format( "GaudiH1D_{}_{}_Axis0", i, i * i ), "( 21, -10.5, 10.5, \"X\" )" ).ignore(); + } + for ( unsigned int i = 0; i < 5; i++ ) { histo1d[i].createHistogram( algo ); } for ( unsigned int i = 0; i < 5; i++ ) ++histo1d[i][-10.0]; // fill the first (non-overflow) bin for ( unsigned int i = 0; i < 5; i++ ) BOOST_TEST( toJSON( histo1d[i] ).at( "bins" )[1] == 1 ); BOOST_TEST( toJSON( histo1d[3] ).at( "title" ) == "Title number 3 of Histogram arrays of 5 histograms in total" ); diff --git a/GaudiKernel/tests/src/CounterHistosUnitTest.cpp b/GaudiKernel/tests/src/CounterHistosUnitTest.cpp index e7e1285275f6ece44c8f231979e901160bee63e8..ed9c1e7dd3bf9a31f94cf12dc6fc037fca4d620f 100644 --- a/GaudiKernel/tests/src/CounterHistosUnitTest.cpp +++ b/GaudiKernel/tests/src/CounterHistosUnitTest.cpp @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2021 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -10,14 +10,20 @@ \***********************************************************************************/ #define BOOST_TEST_DYN_LINK #define BOOST_TEST_MODULE test_CounterHistos + +#include <Gaudi/Accumulators/Histogram.h> #include <Gaudi/Accumulators/RootHistogram.h> +#include <GaudiKernel/PropertyHolder.h> +#include <GaudiKernel/StateMachine.h> + #include "LogHistogram.h" #include <boost/test/unit_test.hpp> #include <deque> #include <iostream> +#include <vector> namespace { @@ -27,10 +33,21 @@ namespace { MonitoringHub& monitoringHub() { return m_monitHub; } MonitoringHub m_monitHub{}; }; - struct Algo { - ServiceLocator* serviceLocator() { return &m_serviceLocator; } - std::string name() { return ""; } - ServiceLocator m_serviceLocator{}; + struct BaseAlgo : INamedInterface, IProperty { + unsigned long addRef() override { return 0; }; + unsigned long release() override { return 0; }; + void* i_cast( const InterfaceID& ) const override { return nullptr; } + std::vector<std::string> getInterfaceNames() const override { return {}; } + unsigned long refCount() const override { return 0; } + StatusCode queryInterface( const InterfaceID&, void** ) override { return StatusCode::FAILURE; }; + const std::string& name() const override { return m_name; }; + std::string m_name{}; + }; + struct Algo : PropertyHolder<BaseAlgo> { + ServiceLocator* serviceLocator() { return &m_serviceLocator; } + ServiceLocator m_serviceLocator{}; + void registerCallBack( Gaudi::StateMachine::Transition, std::function<void()> ) {} + Gaudi::StateMachine::State FSMState() const { return Gaudi::StateMachine::CONFIGURED; } }; struct HistSink : public Gaudi::Monitoring::Hub::Sink { virtual void registerEntity( Gaudi::Monitoring::Hub::Entity ent ) override { m_entities.push_back( ent ); } @@ -53,7 +70,8 @@ BOOST_AUTO_TEST_CASE( test_counter_histos, *boost::unit_test::tolerance( 1e-14 ) Algo algo; { - Gaudi::Accumulators::Histogram<1> histo1d{ &algo, "GaudiH1D", "A Gaudi 1D histogram", { 21, -10.5, 10.5, "X" } }; + Gaudi::Accumulators::StaticHistogram<1> histo1d{ + &algo, "GaudiH1D", "A Gaudi 1D histogram", { 21, -10.5, 10.5, "X" } }; ++histo1d[-10.0]; // fill the first (non-overflow) bin BOOST_TEST( toJSON( histo1d ).at( "bins" )[1] == 1 ); #pragma GCC diagnostic push @@ -63,7 +81,7 @@ BOOST_AUTO_TEST_CASE( test_counter_histos, *boost::unit_test::tolerance( 1e-14 ) BOOST_TEST( toJSON( histo1d ).at( "bins" )[1] == 2 ); } { - Gaudi::Accumulators::Histogram<2> histo2d{ + Gaudi::Accumulators::StaticHistogram<2> histo2d{ &algo, "GaudiH2D", "A Gaudi 2D histogram", { { 21, -10.5, 10.5, "X" }, { 21, -10.5, 10.5, "Y" } } }; ++histo2d[{ -10.0, -10.0 }]; // fill the first (non-overflow) bin BOOST_TEST( toJSON( histo2d ).at( "bins" )[( 1 + 21 + 1 ) + 1] == 1 ); @@ -74,7 +92,7 @@ BOOST_AUTO_TEST_CASE( test_counter_histos, *boost::unit_test::tolerance( 1e-14 ) BOOST_TEST( toJSON( histo2d ).at( "bins" )[( 1 + 21 + 1 ) + 1] == 2 ); } { - Gaudi::Accumulators::WeightedHistogram<1> histo1dw{ &algo, "", "", { 21, -10.5, 10.5, "X" } }; + Gaudi::Accumulators::StaticWeightedHistogram<1> histo1dw{ &algo, "", "", { 21, -10.5, 10.5, "X" } }; histo1dw[-10.0] += 0.25; // fill the first (non-overflow) bin BOOST_TEST( toJSON( histo1dw ).at( "bins" )[1] == 0.25 ); BOOST_TEST( toJSON( histo1dw ).at( "nEntries" ).get<unsigned long>() == 1 ); @@ -86,7 +104,7 @@ BOOST_AUTO_TEST_CASE( test_counter_histos, *boost::unit_test::tolerance( 1e-14 ) #pragma GCC diagnostic pop } - Gaudi::Accumulators::ProfileHistogram<1u> histo{ &algo, "GaudiP1D", "A Gaudi 1D Profile", { 10, 0, 100 } }; + Gaudi::Accumulators::StaticProfileHistogram<1u> histo{ &algo, "GaudiP1D", "A Gaudi 1D Profile", { 10, 0, 100 } }; histo[-0.5] += -0.5; for ( int i = 0; i < 10; i++ ) { histo[10.0 * double( i ) + 0.5] += double( i ); } @@ -121,7 +139,8 @@ BOOST_AUTO_TEST_CASE( test_counter_root_histos, *boost::unit_test::tolerance( 1e Algo algo; { - Gaudi::Accumulators::RootHistogram<1> histo1d{ &algo, "GaudiH1D", "A Gaudi 1D histogram", { 4, -10, 10, "X" } }; + Gaudi::Accumulators::StaticRootHistogram<1> histo1d{ + &algo, "GaudiH1D", "A Gaudi 1D histogram", { 4, -10, 10, "X" } }; ++histo1d[1]; BOOST_TEST( toJSON( histo1d ).at( "bins" )[3] == 1 ); ++histo1d[2]; @@ -135,7 +154,7 @@ BOOST_AUTO_TEST_CASE( test_counter_root_histos, *boost::unit_test::tolerance( 1e BOOST_TEST( toJSON( histo1d ).at( "standard_deviation" ) == 0.5 ); } { - Gaudi::Accumulators::RootHistogram<2> histo2d{ + Gaudi::Accumulators::StaticRootHistogram<2> histo2d{ &algo, "GaudiH2D", "A Gaudi 2D histogram", { { 4, -10, 10, "X" }, { 4, -10, 10, "Y" } } }; ++histo2d[{ 1, 1 }]; // fill the first (non-overflow) bin BOOST_TEST( toJSON( histo2d ).at( "bins" )[( 1 + 4 + 1 ) * 3 + 3] == 1 ); @@ -155,8 +174,8 @@ BOOST_AUTO_TEST_CASE( test_counter_root_histos, *boost::unit_test::tolerance( 1e BOOST_TEST( toJSON( histo2d ).at( "standard_deviationy" ) == 1 ); } { - Gaudi::Accumulators::RootHistogram<3> histo3d{ - &algo, "GaudiH2D", "A Gaudi 3D histogram", { { 4, -10, 10, "X" }, { 4, -10, 10, "Y" }, { 4, -10, 10, "Z" } } }; + Gaudi::Accumulators::StaticRootHistogram<3> histo3d{ + &algo, "GaudiH3D", "A Gaudi 3D histogram", { { 4, -10, 10, "X" }, { 4, -10, 10, "Y" }, { 4, -10, 10, "Z" } } }; ++histo3d[{ 1, 1, 1 }]; // fill the first (non-overflow) bin BOOST_TEST( toJSON( histo3d ).at( "bins" )[36 * 3 + 6 * 3 + 3] == 1 ); ++histo3d[{ 2, 3, 2 }]; // fill the first (non-overflow) bin @@ -184,8 +203,8 @@ BOOST_AUTO_TEST_CASE( test_counter_root_histos, *boost::unit_test::tolerance( 1e BOOST_AUTO_TEST_CASE( test_integer_histos ) { using namespace Gaudi::Accumulators; - Algo algo; - Gaudi::Accumulators::Histogram<1, atomicity::none, unsigned int> histo{ + Algo algo; + Gaudi::Accumulators::StaticHistogram<1, atomicity::none, unsigned int> histo{ &algo, "IntH1D", "A 1D histogram with integer content", { 10, 0, 10, "X" } }; ++histo[1]; // fill the second (non-overflow) bin ++histo[3]; // fill the fourth (non-overflow) bin @@ -197,8 +216,8 @@ BOOST_AUTO_TEST_CASE( test_integer_histos ) { BOOST_AUTO_TEST_CASE( test_integer_histos_small_ratio ) { using namespace Gaudi::Accumulators; - Algo algo; - Gaudi::Accumulators::Histogram<1, atomicity::none, unsigned int> histo{ + Algo algo; + Gaudi::Accumulators::StaticHistogram<1, atomicity::none, unsigned int> histo{ &algo, "IntH1D", "A 1D histogram with integer content", { 5, 0, 10, "X" } }; ++histo[1]; // fill the first (non-overflow) bin ++histo[3]; // fill the second (non-overflow) bin @@ -208,31 +227,58 @@ BOOST_AUTO_TEST_CASE( test_integer_histos_small_ratio ) { } enum class TestEnum { A, B, C, D }; +std::ostream& operator<<( std::ostream& o, TestEnum v ) { + switch ( v ) { + case TestEnum::A: + o << 'A'; + break; + case TestEnum::B: + o << 'B'; + break; + case TestEnum::C: + o << 'C'; + break; + case TestEnum::D: + o << 'D'; + break; + } + return o; +} namespace Gaudi::Accumulators { - template <> - struct Axis<TestEnum> { + struct EnumAxis { + using ArithmeticType = TestEnum; // helper to make the code less verbose using storage_t = std::underlying_type_t<TestEnum>; // nothing to specify in the constructor as everything is fixed in the enum - Axis() = default; + EnumAxis() = default; - unsigned int nBins = 4; + unsigned int numBins() const { return 4; } storage_t minValue = static_cast<storage_t>( TestEnum::A ), maxValue = static_cast<storage_t>( TestEnum::D ); std::string title = "TestEnum"; std::vector<std::string> labels = { "A", "B", "C", "D" }; // convert the enum value to the index in the bins, taking into account the underflow bin - unsigned int index( TestEnum val ) const { return static_cast<storage_t>( val ) + 1; } + unsigned int index( TestEnum val ) const { return static_cast<storage_t>( val ) + 1; } + friend std::ostream& operator<<( std::ostream& o, EnumAxis const& axis ) { + return o << axis.numBins() << " " << axis.minValue << " " << axis.maxValue; + } }; + void to_json( nlohmann::json& j, const EnumAxis& axis ) { + j = nlohmann::json{ { "nBins", axis.numBins() }, + { "minValue", axis.minValue }, + { "maxValue", axis.maxValue }, + { "title", axis.title } }; + j["labels"] = axis.labels; + } + StatusCode parse( EnumAxis&, const std::string& ) { return StatusCode::SUCCESS; } } // namespace Gaudi::Accumulators BOOST_AUTO_TEST_CASE( test_custom_axis ) { using namespace Gaudi::Accumulators; Algo algo; - // note that for default constructible axis, we have to specify {{}} otherwise - // it is interpreted as an empty array of axes (instead of using the constructor for a single axis) - Histogram<1, atomicity::full, TestEnum> hist{ &algo, "TestEnumHist", "TestEnum histogram", Axis<TestEnum>{} }; + StaticHistogram<1, atomicity::full, TestEnum, std::tuple<EnumAxis>> hist{ &algo, "TestEnumHist", "TestEnum histogram", + EnumAxis{} }; hist[TestEnum::A] += 1; ++hist[TestEnum::B]; @@ -255,7 +301,35 @@ BOOST_AUTO_TEST_CASE( test_custom_axis ) { BOOST_TEST( j["axis"][0]["labels"] == expected_labels ); } -BOOST_AUTO_TEST_CASE( test_custom ) { +BOOST_AUTO_TEST_CASE( test_mixed_axis ) { + using namespace Gaudi::Accumulators; + Algo algo; + + StaticHistogram<2, atomicity::full, float, std::tuple<EnumAxis, Axis<float>>> hist{ + &algo, "TestEnumHist", "TestEnum 2D histogram", EnumAxis{}, { 3, 0, 3, "Value" } }; + + hist[{ TestEnum::A, 0.5 }] += 1; + ++hist[{ TestEnum::B, 1.5 }]; + hist[{ TestEnum::C, 2.5 }] += 2; + + auto j = toJSON( hist ); + + auto bins = j["bins"]; + BOOST_TEST( bins[6 + 1] == 1 ); + BOOST_TEST( bins[6 * 2 + 2] == 1 ); + BOOST_TEST( bins[6 * 3 + 3] == 2 ); + BOOST_TEST( bins[6 * 3 + 3] == 2 ); + + BOOST_TEST( j["axis"][0]["nBins"] == 4 ); + BOOST_TEST( j["axis"][0]["title"] == "TestEnum" ); + BOOST_TEST( j["axis"][1]["nBins"] == 3 ); + BOOST_TEST( j["axis"][1]["title"] == "Value" ); + + nlohmann::json expected_labels = { "A", "B", "C", "D" }; + BOOST_TEST( j["axis"][0]["labels"] == expected_labels ); +} + +BOOST_AUTO_TEST_CASE( test_custom_input ) { using namespace Gaudi::Accumulators; Algo algo; @@ -280,7 +354,7 @@ BOOST_AUTO_TEST_CASE( test_custom ) { BOOST_TEST( j["axis"][0]["nBins"] == 4 ); } -BOOST_AUTO_TEST_CASE( test_custom_2d ) { +BOOST_AUTO_TEST_CASE( test_custom_input_2d ) { using namespace Gaudi::Accumulators; Algo algo; @@ -310,8 +384,10 @@ BOOST_AUTO_TEST_CASE( test_histos_merge_reset, *boost::unit_test::tolerance( 1e- HistSink histSink; algo.serviceLocator()->monitoringHub().addSink( &histSink ); - Histogram<1, atomicity::full, float> hist1( &algo, "TestHist1", "Test histogram 1", Axis<float>{ 10, 0., 10. } ); - Histogram<1, atomicity::full, float> hist2( &algo, "TestHist2", "Test histogram 2", Axis<float>{ 10, 0., 10. } ); + StaticHistogram<1, atomicity::full, float> hist1( &algo, "TestHist1", "Test histogram 1", + Axis<float>{ 10, 0., 10. } ); + StaticHistogram<1, atomicity::full, float> hist2( &algo, "TestHist2", "Test histogram 2", + Axis<float>{ 10, 0., 10. } ); // test all combinations of entities made (a) as part of the Histogram's constructor and (b) separately auto ent1a = histSink.m_entities[0]; @@ -346,7 +422,7 @@ BOOST_AUTO_TEST_CASE( test_2d_histos, *boost::unit_test::tolerance( 1e-14 ) ) { Algo algo; // test filling a 2D histogram with more bins in x than y // Buffer will overflow if the wrong axis' nBins is used to calculate the bin index, resulting in a double free - Histogram<2, atomicity::full, float> hist{ + StaticHistogram<2, atomicity::full, float> hist{ &algo, "Test2DHist", "Test 2D histogram", { 64, 0., 64. }, { 52, 0., 52. } }; for ( int i = 0; i < 64; ++i ) { @@ -362,7 +438,7 @@ BOOST_AUTO_TEST_CASE( test_2d_histos_unique_ptr, *boost::unit_test::tolerance( 1 using namespace Gaudi::Accumulators; Algo algo; // test constructing a 2D histogram inside a deque via emplace_back - std::deque<Histogram<2, atomicity::full, float>> histos; + std::deque<StaticHistogram<2, atomicity::full, float>> histos; histos.emplace_back( &algo, "Test2DHist", "Test 2D histogram", Axis<float>{ 10, 0., 10. }, Axis<float>{ 10, 0., 10. } ); { @@ -376,3 +452,238 @@ BOOST_AUTO_TEST_CASE( test_2d_histos_unique_ptr, *boost::unit_test::tolerance( 1 auto nEntries = j.at( "nEntries" ).get<unsigned long>(); BOOST_TEST( nEntries == 100 ); } + +BOOST_AUTO_TEST_CASE( test_custom_Histos, *boost::unit_test::tolerance( 1e-14 ) ) { + Algo algo; + + { + Gaudi::Accumulators::Histogram<1> histo1d{ &algo, "GaudiH1D" }; + algo.setProperty( "GaudiH1D_Title", "A Gaudi 1D histogram" ).ignore(); + algo.setProperty( "GaudiH1D_Axis0", "( 21, -10.5, 10.5, \"X\" )" ).ignore(); + histo1d.createHistogram( algo ); + auto jHisto = toJSON( histo1d ); + BOOST_TEST( jHisto.at( "title" ) == "A Gaudi 1D histogram" ); + auto axis = jHisto.at( "axis" )[0]; + std::cout << jHisto << std::endl; + BOOST_TEST( axis.at( "nBins" ) == 21 ); + BOOST_TEST( axis.at( "minValue" ) == -10.5 ); + BOOST_TEST( axis.at( "maxValue" ) == 10.5 ); + BOOST_TEST( axis.at( "title" ) == "X" ); + BOOST_TEST( jHisto.at( "bins" ).get<std::vector<double>>().size() == 23 ); + ++histo1d[-10.0]; // fill the first (non-overflow) bin + std::cout << toJSON( histo1d ) << std::endl; + BOOST_TEST( toJSON( histo1d ).at( "bins" )[1] == 1 ); + ++histo1d[-10.0]; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo1d ).at( "bins" )[1] == 2 ); + } + { + Gaudi::Accumulators::Histogram<2> histo2d{ &algo, "GaudiH2D" }; + algo.setProperty( "GaudiH2D_Title", "A Gaudi 2D histogram" ).ignore(); + algo.setProperty( "GaudiH2D_Axis0", "( 21, -10.5, 10.5, \"X\" )" ).ignore(); + algo.setProperty( "GaudiH2D_Axis1", "( 41, -20.5, 20.5 )" ).ignore(); + histo2d.createHistogram( algo ); + auto jHisto = toJSON( histo2d ); + BOOST_TEST( jHisto.at( "title" ) == "A Gaudi 2D histogram" ); + auto axis0 = jHisto.at( "axis" )[0]; + BOOST_TEST( axis0.at( "nBins" ) == 21 ); + BOOST_TEST( axis0.at( "minValue" ) == -10.5 ); + BOOST_TEST( axis0.at( "maxValue" ) == 10.5 ); + BOOST_TEST( axis0.at( "title" ) == "X" ); + auto axis1 = jHisto.at( "axis" )[1]; + BOOST_TEST( axis1.at( "nBins" ) == 41 ); + BOOST_TEST( axis1.at( "minValue" ) == -20.5 ); + BOOST_TEST( axis1.at( "maxValue" ) == 20.5 ); + BOOST_TEST( axis1.at( "title" ) == "" ); + BOOST_TEST( jHisto.at( "bins" ).get<std::vector<double>>().size() == 43 * 23 ); + ++histo2d[{ -10.0, -20.0 }]; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo2d ).at( "bins" )[( 1 + 21 + 1 ) + 1] == 1 ); + ++histo2d[{ -10.0, -20.0 }]; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo2d ).at( "bins" )[( 1 + 21 + 1 ) + 1] == 2 ); + } + { + Gaudi::Accumulators::WeightedHistogram<1> histo1dw{ &algo, "GaudiH1DW" }; + algo.setProperty( "GaudiH1dw_Title", "A Gaudi 1DW histogram" ).ignore(); + algo.setProperty( "GaudiH1dw_Axis0", "( 21, -10.5, 10.5, \"X\" )" ).ignore(); + histo1dw.createHistogram( algo ); + auto jHisto = toJSON( histo1dw ); + BOOST_TEST( jHisto.at( "title" ) == "A Gaudi 1DW histogram" ); + auto axis0 = jHisto.at( "axis" )[0]; + BOOST_TEST( axis0.at( "nBins" ) == 21 ); + BOOST_TEST( axis0.at( "minValue" ) == -10.5 ); + BOOST_TEST( axis0.at( "maxValue" ) == 10.5 ); + BOOST_TEST( axis0.at( "title" ) == "X" ); + histo1dw[-10.0] += 0.25; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo1dw ).at( "bins" )[1] == 0.25 ); + histo1dw[-10.0] += 0.5; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo1dw ).at( "bins" )[1] == 0.75 ); + } + { + Gaudi::Accumulators::ProfileHistogram<1u> histo{ &algo, "GaudiP1D" }; + algo.setProperty( "GaudiP1D_Title", "A Gaudi Profile histogram" ).ignore(); + algo.setProperty( "GaudiP1D_Axis0", "( 10, 0, 100 )" ).ignore(); + histo.createHistogram( algo ); + histo[-0.5] += -0.5; + for ( int i = 0; i < 10; i++ ) { histo[10.0 * double( i ) + 0.5] += double( i ); } + histo[120.0] += 120.0; + nlohmann::json j = toJSON( histo ); + auto nEntries = j.at( "nEntries" ).get<unsigned long>(); + BOOST_TEST( nEntries == 12 ); + auto bincont = j.at( "bins" ).get<std::vector<std::tuple<std::tuple<unsigned int, double>, double>>>(); + for ( size_t i = 0; i < bincont.size(); i++ ) { + auto& [tmp, sumw2] = bincont[i]; + auto& [nent, sumw] = tmp; + if ( i == 0 ) { + BOOST_TEST( nent == 1 ); + BOOST_TEST( sumw == -0.5 ); + BOOST_TEST( sumw2 == 0.5 * 0.5 ); + } else if ( i == bincont.size() - 1 ) { + BOOST_TEST( nent == 1 ); + BOOST_TEST( sumw == 120.0 ); + BOOST_TEST( sumw2 == 120.0 * 120.0 ); + } else { + auto w = i - 1; + BOOST_TEST( nent == 1 ); + BOOST_TEST( sumw == w ); + BOOST_TEST( sumw2 == w * w ); + } + } + } +} + +BOOST_AUTO_TEST_CASE( test_custom_root_histos, *boost::unit_test::tolerance( 1e-14 ) ) { + Algo algo; + + { + Gaudi::Accumulators::RootHistogram<1> histo1d{ &algo, "GaudiH1D" }; + algo.setProperty( "GaudiH1D_Title", "A Gaudi 1D histogram" ).ignore(); + algo.setProperty( "GaudiH1D_Axis0", "( 4, -10, 10, \"X\" )" ).ignore(); + histo1d.createHistogram( algo ); + ++histo1d[1]; + BOOST_TEST( toJSON( histo1d ).at( "bins" )[3] == 1 ); + ++histo1d[2]; + ++histo1d[-100]; + BOOST_TEST( toJSON( histo1d ).at( "bins" )[3] == 2 ); + BOOST_TEST( toJSON( histo1d ).at( "nEntries" ) == 3 ); + BOOST_TEST( toJSON( histo1d ).at( "nTotEntries" ) == 2 ); + BOOST_TEST( toJSON( histo1d ).at( "sum" ) == 3 ); + BOOST_TEST( toJSON( histo1d ).at( "sum2" ) == 5 ); + BOOST_TEST( toJSON( histo1d ).at( "mean" ) == 1.5 ); + BOOST_TEST( toJSON( histo1d ).at( "standard_deviation" ) == 0.5 ); + } + { + Gaudi::Accumulators::RootHistogram<2> histo2d{ &algo, "GaudiH2D" }; + algo.setProperty( "GaudiH2D_Title", "A Gaudi 2D histogram" ).ignore(); + algo.setProperty( "GaudiH2D_Axis0", "( 4, -10, 10, \"X\" )" ).ignore(); + algo.setProperty( "GaudiH2D_Axis1", "( 4, -10, 10, \"Y\" )" ).ignore(); + histo2d.createHistogram( algo ); + ++histo2d[{ 1, 1 }]; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo2d ).at( "bins" )[( 1 + 4 + 1 ) * 3 + 3] == 1 ); + ++histo2d[{ 2, 3 }]; // fill the first (non-overflow) bin + ++histo2d[{ -100, -100 }]; + BOOST_TEST( toJSON( histo2d ).at( "bins" )[( 1 + 4 + 1 ) * 3 + 3] == 2 ); + BOOST_TEST( toJSON( histo2d ).at( "nEntries" ) == 3 ); + BOOST_TEST( toJSON( histo2d ).at( "nTotEntries" ) == 2 ); + BOOST_TEST( toJSON( histo2d ).at( "sumx" ) == 3 ); + BOOST_TEST( toJSON( histo2d ).at( "sumy" ) == 4 ); + BOOST_TEST( toJSON( histo2d ).at( "sumx2" ) == 5 ); + BOOST_TEST( toJSON( histo2d ).at( "sumy2" ) == 10 ); + BOOST_TEST( toJSON( histo2d ).at( "sumxy" ) == 7 ); + BOOST_TEST( toJSON( histo2d ).at( "meanx" ) == 1.5 ); + BOOST_TEST( toJSON( histo2d ).at( "meany" ) == 2 ); + BOOST_TEST( toJSON( histo2d ).at( "standard_deviationx" ) == .5 ); + BOOST_TEST( toJSON( histo2d ).at( "standard_deviationy" ) == 1 ); + } + { + Gaudi::Accumulators::RootHistogram<3> histo3d{ &algo, "GaudiH3D" }; + algo.setProperty( "GaudiH3D_Title", "A Gaudi 3D histogram" ).ignore(); + algo.setProperty( "GaudiH3D_Axis0", "( 4, -10, 10, \"X\" )" ).ignore(); + algo.setProperty( "GaudiH3D_Axis1", "( 4, -10, 10, \"Y\" )" ).ignore(); + algo.setProperty( "GaudiH3D_Axis2", "( 4, -10, 10, \"Z\" )" ).ignore(); + histo3d.createHistogram( algo ); + ++histo3d[{ 1, 1, 1 }]; // fill the first (non-overflow) bin + BOOST_TEST( toJSON( histo3d ).at( "bins" )[36 * 3 + 6 * 3 + 3] == 1 ); + ++histo3d[{ 2, 3, 2 }]; // fill the first (non-overflow) bin + ++histo3d[{ -100, -100, -100 }]; + BOOST_TEST( toJSON( histo3d ).at( "bins" )[36 * 3 + 6 * 3 + 3] == 2 ); + BOOST_TEST( toJSON( histo3d ).at( "nEntries" ) == 3 ); + BOOST_TEST( toJSON( histo3d ).at( "nTotEntries" ) == 2 ); + BOOST_TEST( toJSON( histo3d ).at( "sumx" ) == 3 ); + BOOST_TEST( toJSON( histo3d ).at( "sumy" ) == 4 ); + BOOST_TEST( toJSON( histo3d ).at( "sumz" ) == 3 ); + BOOST_TEST( toJSON( histo3d ).at( "sumx2" ) == 5 ); + BOOST_TEST( toJSON( histo3d ).at( "sumy2" ) == 10 ); + BOOST_TEST( toJSON( histo3d ).at( "sumz2" ) == 5 ); + BOOST_TEST( toJSON( histo3d ).at( "sumxy" ) == 7 ); + BOOST_TEST( toJSON( histo3d ).at( "sumxz" ) == 5 ); + BOOST_TEST( toJSON( histo3d ).at( "sumyz" ) == 7 ); + BOOST_TEST( toJSON( histo3d ).at( "meanx" ) == 1.5 ); + BOOST_TEST( toJSON( histo3d ).at( "meany" ) == 2 ); + BOOST_TEST( toJSON( histo3d ).at( "meanz" ) == 1.5 ); + BOOST_TEST( toJSON( histo3d ).at( "standard_deviationx" ) == .5 ); + BOOST_TEST( toJSON( histo3d ).at( "standard_deviationy" ) == 1 ); + BOOST_TEST( toJSON( histo3d ).at( "standard_deviationz" ) == .5 ); + } +} + +BOOST_AUTO_TEST_CASE( test_custom_integer_histos ) { + using namespace Gaudi::Accumulators; + Algo algo; + Histogram<1, atomicity::none, unsigned int> histo{ &algo, "IntH1D" }; + algo.setProperty( "IntH1D_Title", "A 1D histogram with integer content" ).ignore(); + algo.setProperty( "IntH1D_Axis0", "( 10, 0, 10, \"X\" )" ).ignore(); + histo.createHistogram( algo ); + ++histo[1]; // fill the second (non-overflow) bin + ++histo[3]; // fill the fourth (non-overflow) bin + auto j = toJSON( histo ); + BOOST_TEST( j.at( "bins" )[0] == 0 ); + BOOST_TEST( j.at( "bins" )[2] == 1 ); + BOOST_TEST( j.at( "bins" )[4] == 1 ); +} + +BOOST_AUTO_TEST_CASE( test_custom_custom_axis ) { + using namespace Gaudi::Accumulators; + Algo algo; + Histogram<1, atomicity::full, TestEnum, std::tuple<EnumAxis>> hist{ &algo, "TestEnumHist" }; + algo.setProperty( "IntH1D_Title", "TestEnum histogram" ).ignore(); + algo.setProperty( "IntH1D_Axis0", "Junk - this is ignored" ).ignore(); + hist.createHistogram( algo ); + hist[TestEnum::A] += 1; + ++hist[TestEnum::B]; + hist[TestEnum::C] += 2; + auto j = toJSON( hist ); + auto bins = j["bins"]; + BOOST_TEST( bins[0] == 0 ); + BOOST_TEST( bins[1] == 1 ); + BOOST_TEST( bins[2] == 1 ); + BOOST_TEST( bins[3] == 2 ); + BOOST_TEST( bins[4] == 0 ); + BOOST_TEST( bins[5] == 0 ); + BOOST_TEST( j["axis"][0]["nBins"] == 4 ); + BOOST_TEST( j["axis"][0]["title"] == "TestEnum" ); + nlohmann::json expected_labels = { "A", "B", "C", "D" }; + BOOST_TEST( j["axis"][0]["labels"] == expected_labels ); +} + +BOOST_AUTO_TEST_CASE( test_custom_mixed_axis ) { + using namespace Gaudi::Accumulators; + Algo algo; + Histogram<2, atomicity::full, float, std::tuple<EnumAxis, Axis<float>>> hist{ &algo, "TestEnumHist" }; + algo.setProperty( "TestEnumHist_Title", "TestEnum 2D histogram" ).ignore(); + algo.setProperty( "TestEnumHist_Axis1", "( 3, 0, 3, \"Value\" )" ).ignore(); + hist.createHistogram( algo ); + hist[{ TestEnum::A, 0.5 }] += 1; + ++hist[{ TestEnum::B, 1.5 }]; + hist[{ TestEnum::C, 2.5 }] += 2; + auto j = toJSON( hist ); + auto bins = j["bins"]; + BOOST_TEST( bins[6 + 1] == 1 ); + BOOST_TEST( bins[6 * 2 + 2] == 1 ); + BOOST_TEST( bins[6 * 3 + 3] == 2 ); + BOOST_TEST( bins[6 * 3 + 3] == 2 ); + BOOST_TEST( j["axis"][0]["nBins"] == 4 ); + BOOST_TEST( j["axis"][0]["title"] == "TestEnum" ); + BOOST_TEST( j["axis"][1]["nBins"] == 3 ); + BOOST_TEST( j["axis"][1]["title"] == "Value" ); + nlohmann::json expected_labels = { "A", "B", "C", "D" }; + BOOST_TEST( j["axis"][0]["labels"] == expected_labels ); +} diff --git a/GaudiKernel/tests/src/LogHistogram.h b/GaudiKernel/tests/src/LogHistogram.h index a6fd493ef5cf54b22e6c97fb0f5f0eefef05cd8d..68a412096f7913281f5b66b701244358a0318b1e 100644 --- a/GaudiKernel/tests/src/LogHistogram.h +++ b/GaudiKernel/tests/src/LogHistogram.h @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2022 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -10,7 +10,7 @@ \***********************************************************************************/ #pragma once -#include <Gaudi/Accumulators/Histogram.h> +#include <Gaudi/Accumulators/StaticHistogram.h> namespace Gaudi::Accumulators { @@ -28,38 +28,50 @@ namespace Gaudi::Accumulators { */ template <typename Arithmetic, unsigned int ND> struct LogInputType : std::array<Arithmetic, ND> { + using ArithmeticType = Arithmetic; // allow construction from set of values template <class... ARGS> LogInputType( ARGS... args ) : std::array<Arithmetic, ND>{ static_cast<Arithmetic>( args )... } {} - using ValueType = LogInputType<Arithmetic, ND>; - using AxisArithmeticType = Arithmetic; - unsigned int computeIndex( const std::array<Axis<Arithmetic>, ND>& axis ) const { - unsigned int index = 0; - for ( unsigned int j = 0; j < ND; j++ ) { - unsigned int dim = ND - j - 1; - // compute local index for a given dimension - int localIndex = axis[dim].index( log( ( *this )[dim] ) ); - // compute global index. Bins are stored in a column first manner - index *= ( axis[dim].nBins + 2 ); - index += localIndex; - } - return index; + using ValueType = LogInputType<Arithmetic, ND>; + template <class... AxisType, typename = typename std::enable_if_t<( sizeof...( AxisType ) == ND )>> + unsigned int computeIndex( std::tuple<AxisType...> const& axis ) const { + return computeIndexInternal<0, std::tuple<AxisType...>>( axis ); } auto forInternalCounter() { return 1ul; } - template <typename AxisType, long unsigned NAxis> - static unsigned int computeTotNBins( std::array<AxisType, NAxis> axis ) { - unsigned int nTotBins = 1; - for ( unsigned int i = 0; i < NAxis; i++ ) { nTotBins *= ( axis[i].nBins + 2 ); } - return nTotBins; + template <class... AxisType, typename = typename std::enable_if_t<( sizeof...( AxisType ) == ND )>> + static unsigned int computeTotNBins( std::tuple<AxisType...> const& axis ) { + return computeTotNBinsInternal<0, std::tuple<AxisType...>>( axis ); + } + + private: + template <int N, class Tuple> + unsigned int computeIndexInternal( Tuple const& allAxis ) const { + // compute global index. Bins are stored in a column first manner + auto const& axis = std::get<N>( allAxis ); + unsigned int localIndex = axis.index( log( std::get<N>( *this ) ) ); + if constexpr ( N + 1 == ND ) + return localIndex; + else + return localIndex + ( axis.numBins() + 2 ) * computeIndexInternal<N + 1, Tuple>( allAxis ); + } + template <int N, class Tuple> + static unsigned int computeTotNBinsInternal( Tuple const& allAxis ) { + auto const& axis = std::get<N>( allAxis ); + unsigned int localNBins = axis.numBins() + 2; + if constexpr ( N + 1 == ND ) + return localNBins; + else + return localNBins * computeTotNBinsInternal<N + 1, Tuple>( allAxis ); } }; - template <atomicity Atomicity, typename Arithmetic, typename ND> + template <atomicity Atomicity, typename Arithmetic, typename ND, typename AxisType> using LogAccumulator = HistogramingAccumulatorInternal<Atomicity, LogInputType<Arithmetic, ND::value>, unsigned long, - ND, CountAccumulator>; + CountAccumulator, AxisType>; constexpr char logHistogramString[] = "histogram:Histogram"; template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double> - using LogHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, logHistogramString, LogAccumulator>; + using LogHistogram = HistogramingCounterBase<ND, Atomicity, Arithmetic, logHistogramString, LogAccumulator, + make_tuple_t<Axis<Arithmetic>, ND>>; } // namespace Gaudi::Accumulators diff --git a/GaudiKernel/tests/src/ProfileHistoSpeedTest.cpp b/GaudiKernel/tests/src/ProfileHistoSpeedTest.cpp index 8aec7c51f3ab332f26ea00e86ed3f5f257b1d613..bc3fc0cb9f9da41d7398b372adbe82a83fb507cb 100644 --- a/GaudiKernel/tests/src/ProfileHistoSpeedTest.cpp +++ b/GaudiKernel/tests/src/ProfileHistoSpeedTest.cpp @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2021 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -11,7 +11,7 @@ #define BOOST_TEST_DYN_LINK #define BOOST_TEST_MODULE test_ProfileHistoSpeed -#include <Gaudi/Accumulators/Histogram.h> +#include <Gaudi/Accumulators/StaticHistogram.h> #include <Gaudi/Histograming/Sink/Utils.h> #include <boost/test/unit_test.hpp> @@ -54,8 +54,8 @@ BOOST_AUTO_TEST_CASE( test_profile_histo_speed, *boost::unit_test::tolerance( 1e Algo algo; std::cout << "Elapsed time (us) :\n"; // create a large profile (1M bins) histogram and fill it with 10 million items - auto startTime = std::chrono::steady_clock::now(); - Gaudi::Accumulators::ProfileHistogram<1u> histo{ + auto startTime = std::chrono::steady_clock::now(); + Gaudi::Accumulators::StaticProfileHistogram<1u> histo{ &algo, "GaudiP1D", "A Large Gaudi 1D Profile", { 1000000, 0, 1000000 } }; for ( int i = 0; i < 10000000; i++ ) { histo[0.1 * double( i ) + 0.05] += double( i * 0.01 ); } auto filledTime = std::chrono::steady_clock::now(); @@ -113,7 +113,7 @@ BOOST_AUTO_TEST_CASE( test_profile_histo_speed, *boost::unit_test::tolerance( 1e << "\n"; // test fast conversion of 2D histos - Gaudi::Accumulators::ProfileHistogram<2u> histo2d{ + Gaudi::Accumulators::StaticProfileHistogram<2u> histo2d{ &algo, "GaudiP3D", "A Large Gaudi 2D Profile", { 1000, 0, 1000 }, { 1000, 0, 1000 } }; for ( int i = 0; i < 1000; i++ ) { for ( int j = 0; j < 1000; j++ ) { @@ -130,8 +130,8 @@ BOOST_AUTO_TEST_CASE( test_profile_histo_speed, *boost::unit_test::tolerance( 1e << "\n"; // test fast conversion of 3D histos - Gaudi::Accumulators::ProfileHistogram<3u> histo3d{ &algo, "GaudiP3D", "A Large Gaudi 3D Profile", - { 100, 0, 100 }, { 100, 0, 100 }, { 100, 0, 100 } }; + Gaudi::Accumulators::StaticProfileHistogram<3u> histo3d{ &algo, "GaudiP3D", "A Large Gaudi 3D Profile", + { 100, 0, 100 }, { 100, 0, 100 }, { 100, 0, 100 } }; for ( int i = 0; i < 100; i++ ) { for ( int j = 0; j < 100; j++ ) { for ( int k = 0; k < 100; k++ ) { diff --git a/GaudiTestSuite/options/Histograms.py b/GaudiTestSuite/options/Histograms.py index 15d75419c80cbda72c59eeadb742648e0a29539a..1a22ca2458aa7ca411b23afc584d1866c950ab00 100644 --- a/GaudiTestSuite/options/Histograms.py +++ b/GaudiTestSuite/options/Histograms.py @@ -1,5 +1,5 @@ ##################################################################################### -# (c) Copyright 1998-2019 CERN for the benefit of the LHCb and ATLAS collaborations # +# (c) Copyright 1998-2024 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". # @@ -23,15 +23,70 @@ from Configurables import ( Gaudi__TestSuite__Counter__GaudiRootHistoAlgorithm as RootCounterHistoAlg, ) -algs = [ - CounterHistoAlg("SimpleCounterHistos", OutputLevel=DEBUG), - RootCounterHistoAlg("SimpleRootCounterHistos", OutputLevel=DEBUG), -] +histoAlg = CounterHistoAlg( + "SimpleCounterHistos", + OutputLevel=DEBUG, + CustomGauss_Title="Gaussian mean=0, sigma=1, atomic", + CustomGauss_Axis0=(100, -5, 5, "X"), + CustomGaussFlat_Title="Gaussian V Flat, atomic", + CustomGaussFlat_Axis0=(50, -5, 5, "X"), + CustomGaussFlat_Axis1=(50, -5, 5, "Y"), + CustomGaussFlatGauss_Title="Gaussian V Flat V Gaussian, atomic", + CustomGaussFlatGauss_Axis0=(10, -5, 5, "X"), + CustomGaussFlatGauss_Axis1=(10, -5, 5, "Y"), + CustomGaussFlatGauss_Axis2=(10, -5, 5, "Z"), + CustomGaussW_Title="Gaussian mean=0, sigma=1, weighted", + CustomGaussW_Axis0=(100, -5, 5), + CustomGaussFlatW_Title="Gaussian V Flat, weighted", + CustomGaussFlatW_Axis0=(50, -5, 5), + CustomGaussFlatW_Axis1=(50, -5, 5), + CustomGaussFlatGaussW_Title="Gaussian V Flat V Gaussian, weighted", + CustomGaussFlatGaussW_Axis0=(10, -5, 5), + CustomGaussFlatGaussW_Axis1=(10, -5, 5), + CustomGaussFlatGaussW_Axis2=(10, -5, 5), + CustomProfGauss_Title="Profile, Gaussian mean=0, sigma=1, atomic", + CustomProfGauss_Axis0=(100, -5, 5), + CustomProfGaussFlat_Title="Profile, Gaussian V Flat, atomic", + CustomProfGaussFlat_Axis0=(50, -5, 5), + CustomProfGaussFlat_Axis1=(50, -5, 5), + CustomProfGaussFlatGauss_Title="Profile, Gaussian V Flat V Gaussian, atomic", + CustomProfGaussFlatGauss_Axis0=(10, -5, 5), + CustomProfGaussFlatGauss_Axis1=(10, -5, 5), + CustomProfGaussFlatGauss_Axis2=(10, -5, 5), + CustomProfGaussW_Title="Profile, Gaussian mean=0, sigma=1, weighted", + CustomProfGaussW_Axis0=(100, -5, 5), + CustomProfGaussFlatW_Title="Profile, Gaussian V Flat, weighted", + CustomProfGaussFlatW_Axis0=(50, -5, 5), + CustomProfGaussFlatW_Axis1=(50, -5, 5), + CustomProfGaussFlatGaussW_Title="Profile, Gaussian V Flat V Gaussian, weighted", + CustomProfGaussFlatGaussW_Axis0=(10, -5, 5), + CustomProfGaussFlatGaussW_Axis1=(10, -5, 5), + CustomProfGaussFlatGaussW_Axis2=(10, -5, 5), + CustomGaussNoInit_Title="Gaussian mean=0, sigma=1, atomic", + CustomGaussNoInit_Axis0=(100, -5, 5, "X"), +) + +rootHistoAlg = RootCounterHistoAlg( + "SimpleRootCounterHistos", + OutputLevel=DEBUG, + CustomGauss_Title="Gaussian mean=0, sigma=1, atomic", + CustomGauss_Axis0=(100, -5, 5, "X"), + CustomGaussFlat_Title="Gaussian V Flat, atomic", + CustomGaussFlat_Axis0=(50, -5, 5, "X"), + CustomGaussFlat_Axis1=(50, -5, 5, "Y"), + CustomGaussFlatGauss_Title="Gaussian V Flat V Gaussian, atomic", + CustomGaussFlatGauss_Axis0=(10, -5, 5, "X"), + CustomGaussFlatGauss_Axis1=(10, -5, 5, "Y"), + CustomGaussFlatGauss_Axis2=(10, -5, 5, "Z"), +) app = ApplicationMgr( EvtMax=50000, EvtSel="NONE", HistogramPersistency="ROOT", - TopAlg=algs, - ExtSvc=[MessageSvcSink(), RootHistoSink()], + TopAlg=[histoAlg, rootHistoAlg], + ExtSvc=[ + MessageSvcSink(TypesToSave=["counter:.*", "histogram:.*"]), + RootHistoSink(), + ], ) diff --git a/GaudiTestSuite/src/Histograms/CounterHistos.cpp b/GaudiTestSuite/src/Histograms/CounterHistos.cpp index 9cce4e1ec2ce8d27ec1257407f95568d0d5ef530..47bd305dfcd7327cddd66793114f3c8c92a9b6d1 100644 --- a/GaudiTestSuite/src/Histograms/CounterHistos.cpp +++ b/GaudiTestSuite/src/Histograms/CounterHistos.cpp @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2019 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -11,6 +11,7 @@ #include <Gaudi/Accumulators/Histogram.h> #include <Gaudi/Accumulators/RootHistogram.h> #include <Gaudi/Algorithm.h> +#include <Gaudi/FSMCallbackHolder.h> #include <GaudiKernel/RndmGenerators.h> #include "../../../GaudiKernel/tests/src/LogHistogram.h" @@ -58,7 +59,7 @@ namespace Gaudi { private: mutable Rndm::Numbers m_rand; - mutable std::deque<Gaudi::Accumulators::Histogram<1, Atomicity, Arithmetic>> m_histos; + mutable std::deque<Gaudi::Accumulators::StaticHistogram<1, Atomicity, Arithmetic>> m_histos; Gaudi::Property<unsigned int> m_nHistos{ this, "NumHistos", 20, "" }; Gaudi::Property<unsigned int> m_nTracks{ this, "NumTracks", 30, "" }; @@ -73,9 +74,14 @@ namespace Gaudi { DECLARE_COMPONENT_WITH_ID( HistoTimingAlgI, "HistoTimingAlgI" ) /// Example of algorithm using histograms accumulators. - class GaudiHistoAlgorithm : public Gaudi::Algorithm { + class GaudiHistoAlgorithm : public FSMCallbackHolder<Gaudi::Algorithm> { public: - using Gaudi::Algorithm::Algorithm; + using Base = FSMCallbackHolder<Gaudi::Algorithm>; + using Base::Base; + + StatusCode initialize() override { + return Base::initialize().andThen( [&] { m_custom_gauss_noinit.createHistogram( *this ); } ); + } StatusCode execute( const EventContext& ) const override { // some random number generators, just to provide numbers @@ -146,6 +152,22 @@ namespace Gaudi { m_log_gauss[gauss]++; m_log_gaussVflat[{ flat, gauss }]++; + // using Custom histograms, with title and axis defined in python + ++m_custom_gauss_noprop[gauss]; + ++m_custom_gauss[gauss]; + ++m_custom_gaussVflat[{ flat, gauss }]; + ++m_custom_gaussVflatVgauss[{ flat, gauss, gauss2 }]; + m_custom_gauss_w[gauss] += .5; + m_custom_gaussVflat_w[{ flat, gauss }] += .5; + m_custom_gaussVflatVgauss_w[{ flat, gauss, gauss2 }] += .5; + m_custom_prof_gauss[gauss] += gauss3; + m_custom_prof_gaussVflat[{ flat, gauss }] += gauss3; + m_custom_prof_gaussVflatVgauss[{ flat, gauss, gauss2 }] += gauss3; + m_custom_prof_gauss_w[gauss] += { gauss3, .5 }; + m_custom_prof_gaussVflat_w[{ flat, gauss }] += { gauss3, .5 }; + m_custom_prof_gaussVflatVgauss_w[{ flat, gauss, gauss2 }] += { gauss3, .5 }; + ++m_custom_gauss_noinit[gauss]; + if ( m_nCalls.nEntries() == 0 ) always() << "Filling Histograms...... Please be patient !" << endmsg; ++m_nCalls; return StatusCode::SUCCESS; @@ -157,55 +179,55 @@ namespace Gaudi { // Testing regular histograms in all dimensions and features // "default" case, that is bins containing doubles and atomic - mutable Gaudi::Accumulators::Histogram<1> m_gauss{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_gauss{ this, "Gauss", "Gaussian mean=0, sigma=1, atomic", { 100, -5, 5, "X" } }; - mutable Gaudi::Accumulators::Histogram<2> m_gaussVflat{ + mutable Gaudi::Accumulators::StaticHistogram<2> m_gaussVflat{ this, "GaussFlat", "Gaussian V Flat, atomic", { { 50, -5, 5, "X" }, { 50, -5, 5, "Y" } } }; - mutable Gaudi::Accumulators::Histogram<3> m_gaussVflatVgauss{ + mutable Gaudi::Accumulators::StaticHistogram<3> m_gaussVflatVgauss{ this, "GaussFlatGauss", "Gaussian V Flat V Gaussian, atomic", { { 10, -5, 5, "X" }, { 10, -5, 5, "Y" }, { 10, -5, 5, "Z" } } }; // non atomic versions - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::none> m_gauss_noato{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::none> m_gauss_noato{ this, "GaussNA", "Gaussian mean=0, sigma=1, non atomic", { 100, -5, 5 } }; - mutable Gaudi::Accumulators::Histogram<2, Gaudi::Accumulators::atomicity::none> m_gaussVflat_noato{ + mutable Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::none> m_gaussVflat_noato{ this, "GaussFlatNA", "Gaussian V Flat, non atomic", { { 50, -5, 5 }, { 50, -5, 5 } } }; - mutable Gaudi::Accumulators::Histogram<3, Gaudi::Accumulators::atomicity::none> m_gaussVflatVgauss_noato{ + mutable Gaudi::Accumulators::StaticHistogram<3, Gaudi::Accumulators::atomicity::none> m_gaussVflatVgauss_noato{ this, "GaussFlatGaussNA", "Gaussian V Flat V Gaussian, non atomic", { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; // using integers - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int> m_gauss_int{ + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_gauss_int{ this, "GaussInt", "Gaussian mean=0, sigma=1, integer values", { 10, -5, 5 } }; - mutable Gaudi::Accumulators::Histogram<2, Gaudi::Accumulators::atomicity::full, int> m_gaussVflat_int{ + mutable Gaudi::Accumulators::StaticHistogram<2, Gaudi::Accumulators::atomicity::full, int> m_gaussVflat_int{ this, "GaussFlatInt", "Gaussian V Flat, integer values", { { 10, -5, 5 }, { 10, -5, 5 } } }; - mutable Gaudi::Accumulators::Histogram<3, Gaudi::Accumulators::atomicity::full, int> m_gaussVflatVgauss_int{ - this, - "GaussFlatGaussInt", - "Gaussian V Flat V Gaussian, interger values", - { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; + mutable Gaudi::Accumulators::StaticHistogram<3, Gaudi::Accumulators::atomicity::full, int> + m_gaussVflatVgauss_int{ this, + "GaussFlatGaussInt", + "Gaussian V Flat V Gaussian, interger values", + { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; // weighted version, "default" case - mutable Gaudi::Accumulators::WeightedHistogram<1> m_gauss_w{ + mutable Gaudi::Accumulators::StaticWeightedHistogram<1> m_gauss_w{ this, "GaussW", "Gaussian mean=0, sigma=1, weighted", { 100, -5, 5 } }; - mutable Gaudi::Accumulators::WeightedHistogram<2> m_gaussVflat_w{ + mutable Gaudi::Accumulators::StaticWeightedHistogram<2> m_gaussVflat_w{ this, "GaussFlatW", "Gaussian V Flat, weighted", { { 50, -5, 5 }, { 50, -5, 5 } } }; - mutable Gaudi::Accumulators::WeightedHistogram<3> m_gaussVflatVgauss_w{ + mutable Gaudi::Accumulators::StaticWeightedHistogram<3> m_gaussVflatVgauss_w{ this, "GaussFlatGaussW", "Gaussian V Flat V Gaussian, weighted", { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; // "default" case, dedicated to testing buffers - mutable Gaudi::Accumulators::Histogram<1> m_gauss_buf{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_gauss_buf{ this, "GaussBuf", "Gaussian mean=0, sigma=1, buffered", { 100, -5, 5 } }; - mutable Gaudi::Accumulators::Histogram<2> m_gaussVflat_buf{ + mutable Gaudi::Accumulators::StaticHistogram<2> m_gaussVflat_buf{ this, "GaussFlatBuf", "Gaussian V Flat, buffered", { { 50, -5, 5 }, { 50, -5, 5 } } }; - mutable Gaudi::Accumulators::Histogram<3> m_gaussVflatVgauss_buf{ + mutable Gaudi::Accumulators::StaticHistogram<3> m_gaussVflatVgauss_buf{ this, "GaussFlatGaussBuf", "Gaussian V Flat V Gaussian, buffered", @@ -214,58 +236,60 @@ namespace Gaudi { // Testing profiling histograms in all dimensions and features // "default" case, that is bins containing doubles and atomic - mutable Gaudi::Accumulators::ProfileHistogram<1> m_prof_gauss{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_prof_gauss{ this, "ProfGauss", "Profile, Gaussian mean=0, sigma=1, atomic", { 100, -5, 5 } }; - mutable Gaudi::Accumulators::ProfileHistogram<2> m_prof_gaussVflat{ + mutable Gaudi::Accumulators::StaticProfileHistogram<2> m_prof_gaussVflat{ this, "ProfGaussFlat", "Profile, Gaussian V Flat, atomic", { { 50, -5, 5 }, { 50, -5, 5 } } }; - mutable Gaudi::Accumulators::ProfileHistogram<3> m_prof_gaussVflatVgauss{ + mutable Gaudi::Accumulators::StaticProfileHistogram<3> m_prof_gaussVflatVgauss{ this, "ProfGaussFlatGauss", "Profile, Gaussian V Flat V Gaussian, atomic", { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; // non atomic versions - mutable Gaudi::Accumulators::ProfileHistogram<1, Gaudi::Accumulators::atomicity::none> m_prof_gauss_noato{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::none> m_prof_gauss_noato{ this, "ProfGaussNA", "Profile, Gaussian mean=0, sigma=1, non atomic", { 100, -5, 5 } }; - mutable Gaudi::Accumulators::ProfileHistogram<2, Gaudi::Accumulators::atomicity::none> m_prof_gaussVflat_noato{ - this, "ProfGaussFlatNA", "Profile, Gaussian V Flat, non atomic", { { 50, -5, 5 }, { 50, -5, 5 } } }; - mutable Gaudi::Accumulators::ProfileHistogram<3, Gaudi::Accumulators::atomicity::none> + mutable Gaudi::Accumulators::StaticProfileHistogram<2, Gaudi::Accumulators::atomicity::none> + m_prof_gaussVflat_noato{ + this, "ProfGaussFlatNA", "Profile, Gaussian V Flat, non atomic", { { 50, -5, 5 }, { 50, -5, 5 } } }; + mutable Gaudi::Accumulators::StaticProfileHistogram<3, Gaudi::Accumulators::atomicity::none> m_prof_gaussVflatVgauss_noato{ this, "ProfGaussFlatGaussNA", "Profile, Gaussian V Flat V Gaussian, non atomic", { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; // using integers internally - mutable Gaudi::Accumulators::ProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> m_prof_gauss_int{ - this, "ProfGaussInt", "Profile, Gaussian mean=0, sigma=1, integer values", { 10, -5, 5 } }; - mutable Gaudi::Accumulators::ProfileHistogram<2, Gaudi::Accumulators::atomicity::full, int> + mutable Gaudi::Accumulators::StaticProfileHistogram<1, Gaudi::Accumulators::atomicity::full, int> + m_prof_gauss_int{ + this, "ProfGaussInt", "Profile, Gaussian mean=0, sigma=1, integer values", { 10, -5, 5 } }; + mutable Gaudi::Accumulators::StaticProfileHistogram<2, Gaudi::Accumulators::atomicity::full, int> m_prof_gaussVflat_int{ this, "ProfGaussFlatInt", "Profile, Gaussian V Flat, integer values", { { 10, -5, 5 }, { 10, -5, 5 } } }; - mutable Gaudi::Accumulators::ProfileHistogram<3, Gaudi::Accumulators::atomicity::full, int> + mutable Gaudi::Accumulators::StaticProfileHistogram<3, Gaudi::Accumulators::atomicity::full, int> m_prof_gaussVflatVgauss_int{ this, "ProfGaussFlatGaussInt", "Profile, Gaussian V Flat V Gaussian, interger values", { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; // weighted version, "default" case - mutable Gaudi::Accumulators::WeightedProfileHistogram<1> m_prof_gauss_w{ + mutable Gaudi::Accumulators::StaticWeightedProfileHistogram<1> m_prof_gauss_w{ this, "ProfGaussW", "Profile, Gaussian mean=0, sigma=1, weighted", { 100, -5, 5 } }; - mutable Gaudi::Accumulators::WeightedProfileHistogram<2> m_prof_gaussVflat_w{ + mutable Gaudi::Accumulators::StaticWeightedProfileHistogram<2> m_prof_gaussVflat_w{ this, "ProfGaussFlatW", "Profile, Gaussian V Flat, weighted", { { 50, -5, 5 }, { 50, -5, 5 } } }; - mutable Gaudi::Accumulators::WeightedProfileHistogram<3> m_prof_gaussVflatVgauss_w{ + mutable Gaudi::Accumulators::StaticWeightedProfileHistogram<3> m_prof_gaussVflatVgauss_w{ this, "ProfGaussFlatGaussW", "Profile, Gaussian V Flat V Gaussian, weighted", { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; // "default" case, dedicated to testing buffers - mutable Gaudi::Accumulators::ProfileHistogram<1> m_prof_gauss_buf{ + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_prof_gauss_buf{ this, "ProfGaussBuf", "Profile, Gaussian mean=0, sigma=1, buffered", { 100, -5, 5 } }; - mutable Gaudi::Accumulators::ProfileHistogram<2> m_prof_gaussVflat_buf{ + mutable Gaudi::Accumulators::StaticProfileHistogram<2> m_prof_gaussVflat_buf{ this, "ProfGaussFlatBuf", "Profile, Gaussian V Flat, buffered", { { 50, -5, 5 }, { 50, -5, 5 } } }; - mutable Gaudi::Accumulators::ProfileHistogram<3> m_prof_gaussVflatVgauss_buf{ + mutable Gaudi::Accumulators::StaticProfileHistogram<3> m_prof_gaussVflatVgauss_buf{ this, "ProfGaussFlatGaussBuf", "Profile, Gaussian V Flat V Gaussian, buffered", @@ -277,15 +301,41 @@ namespace Gaudi { this, "LogGaussFlat", "LogLog, Gaussian V Flat", { { 5, 0, 2 }, { 5, 0, 2 } } }; // Histogram in an absolute location - mutable Gaudi::Accumulators::Histogram<1> m_gaussAbsName{ + mutable Gaudi::Accumulators::StaticHistogram<1> m_gaussAbsName{ this, "/TopDir/SubDir/Gauss", "Gaussian mean=0, sigma=1, atomic", { 100, -5, 5, "X" } }; + + // Custom histograms, specifying title and axis in python + mutable Gaudi::Accumulators::Histogram<1> m_custom_gauss_noprop{ + this, "CustomGaussNoProp", "Gaussian mean=0, sigma=1, atomic", { 100, -5, 5, "X" } }; + mutable Gaudi::Accumulators::Histogram<1> m_custom_gauss{ this, "CustomGauss" }; + mutable Gaudi::Accumulators::Histogram<2> m_custom_gaussVflat{ this, "CustomGaussFlat" }; + mutable Gaudi::Accumulators::Histogram<3> m_custom_gaussVflatVgauss{ this, "CustomGaussFlatGauss" }; + + mutable Gaudi::Accumulators::WeightedHistogram<1> m_custom_gauss_w{ this, "CustomGaussW" }; + mutable Gaudi::Accumulators::WeightedHistogram<2> m_custom_gaussVflat_w{ this, "CustomGaussFlatW" }; + mutable Gaudi::Accumulators::WeightedHistogram<3> m_custom_gaussVflatVgauss_w{ this, "CustomGaussFlatGaussW" }; + + mutable Gaudi::Accumulators::ProfileHistogram<1> m_custom_prof_gauss{ this, "CustomProfGauss" }; + mutable Gaudi::Accumulators::ProfileHistogram<2> m_custom_prof_gaussVflat{ this, "CustomProfGaussFlat" }; + mutable Gaudi::Accumulators::ProfileHistogram<3> m_custom_prof_gaussVflatVgauss{ this, + "CustomProfGaussFlatGauss" }; + + mutable Gaudi::Accumulators::WeightedProfileHistogram<1> m_custom_prof_gauss_w{ this, "CustomProfGaussW" }; + mutable Gaudi::Accumulators::WeightedProfileHistogram<2> m_custom_prof_gaussVflat_w{ this, + "CustomProfGaussFlatW" }; + mutable Gaudi::Accumulators::WeightedProfileHistogram<3> m_custom_prof_gaussVflatVgauss_w{ + this, "CustomProfGaussFlatGaussW" }; + + // Checking DoNotInitialize property of CustomHistogram + mutable Gaudi::Accumulators::Histogram<1> m_custom_gauss_noinit{ this, "CustomGaussNoInit", "", {}, true }; }; DECLARE_COMPONENT( GaudiHistoAlgorithm ) /// Example of algorithm using root histograms accumulators. - class GaudiRootHistoAlgorithm : public Gaudi::Algorithm { + class GaudiRootHistoAlgorithm : public FSMCallbackHolder<Gaudi::Algorithm> { public: - using Gaudi::Algorithm::Algorithm; + using Base = FSMCallbackHolder<Gaudi::Algorithm>; + using Base::Base; StatusCode execute( const EventContext& ) const override { // some random number generators, just to provide numbers @@ -312,6 +362,11 @@ namespace Gaudi { ++gaussVflatVgauss_buf[{ flat, gauss, gauss2 }]; } + // custom histograms, with title and axis defined in python + ++m_custom_gauss[gauss]; + ++m_custom_gaussVflat[{ flat, gauss }]; + ++m_custom_gaussVflatVgauss[{ flat, gauss, gauss2 }]; + if ( m_nCalls.nEntries() == 0 ) always() << "Filling Histograms...... Please be patient !" << endmsg; ++m_nCalls; return StatusCode::SUCCESS; @@ -323,26 +378,34 @@ namespace Gaudi { // Testing Root histograms in all dimensions // "default" case, that is bins containing doubles and atomic - mutable Gaudi::Accumulators::RootHistogram<1> m_gauss{ + mutable Gaudi::Accumulators::StaticRootHistogram<1> m_gauss{ this, "Gauss", "Gaussian mean=0, sigma=1, atomic", { 100, -5, 5, "X" } }; - mutable Gaudi::Accumulators::RootHistogram<2> m_gaussVflat{ + mutable Gaudi::Accumulators::StaticRootHistogram<2> m_gaussVflat{ this, "GaussFlat", "Gaussian V Flat, atomic", { { 50, -5, 5, "X" }, { 50, -5, 5, "Y" } } }; - mutable Gaudi::Accumulators::RootHistogram<3> m_gaussVflatVgauss{ + mutable Gaudi::Accumulators::StaticRootHistogram<3> m_gaussVflatVgauss{ this, "GaussFlatGauss", "Gaussian V Flat V Gaussian, atomic", { { 10, -5, 5, "X" }, { 10, -5, 5, "Y" }, { 10, -5, 5, "Z" } } }; // "default" case, dedicated to testing buffers - mutable Gaudi::Accumulators::RootHistogram<1> m_gauss_buf{ + mutable Gaudi::Accumulators::StaticRootHistogram<1> m_gauss_buf{ this, "GaussBuf", "Gaussian mean=0, sigma=1, buffered", { 100, -5, 5 } }; - mutable Gaudi::Accumulators::RootHistogram<2> m_gaussVflat_buf{ + mutable Gaudi::Accumulators::StaticRootHistogram<2> m_gaussVflat_buf{ this, "GaussFlatBuf", "Gaussian V Flat, buffered", { { 50, -5, 5 }, { 50, -5, 5 } } }; - mutable Gaudi::Accumulators::RootHistogram<3> m_gaussVflatVgauss_buf{ + mutable Gaudi::Accumulators::StaticRootHistogram<3> m_gaussVflatVgauss_buf{ this, "GaussFlatGaussBuf", "Gaussian V Flat V Gaussian, buffered", { { 10, -5, 5 }, { 10, -5, 5 }, { 10, -5, 5 } } }; + + // Custom histogram cases, specifying title and axis in python + mutable Gaudi::Accumulators::RootHistogram<1> m_custom_gauss{ this, "CustomGauss", + "Gaussian mean=0, sigma=1, atomic" }; + mutable Gaudi::Accumulators::RootHistogram<2> m_custom_gaussVflat{ this, "CustomGaussFlat", + "Gaussian V Flat, atomic" }; + mutable Gaudi::Accumulators::RootHistogram<3> m_custom_gaussVflatVgauss{ this, "CustomGaussFlatGauss", + "Gaussian V Flat V Gaussian, atomic" }; }; DECLARE_COMPONENT( GaudiRootHistoAlgorithm ) } // namespace Counter diff --git a/GaudiTestSuite/src/Histograms/HistoProps.cpp b/GaudiTestSuite/src/Histograms/HistoProps.cpp index 318dee28af1abed8205f80eb22aebacc8ebae345..0e230bfdeeaf02a5d44e35547fec3e38e7354d20 100644 --- a/GaudiTestSuite/src/Histograms/HistoProps.cpp +++ b/GaudiTestSuite/src/Histograms/HistoProps.cpp @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 1998-2019 CERN for the benefit of the LHCb and ATLAS collaborations * +* (c) Copyright 1998-2024 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". * @@ -18,7 +18,7 @@ // ============================================================================ // GaudiAlg // ============================================================================ -#include "Gaudi/Accumulators/Histogram.h" +#include "Gaudi/Accumulators/StaticHistogram.h" #include "GaudiKernel/Algorithm.h" #include <fmt/format.h> // ============================================================================ @@ -32,6 +32,8 @@ namespace Gaudi { namespace TestSuite { /** @class HistoProps * simple example, which illustrates the usage of "histogram properties" + * DO NOT USE. This is legacy code, Gaudi::Accumulatos::StaticHistogram support + * properties natively. * @author Vanay BELYAEV ibelyaev@physics.syr.edu * @date 2007-09-18 */ @@ -41,9 +43,8 @@ namespace Gaudi { StatusCode initialize() override { return Algorithm::initialize().andThen( [&] { - ; - using hist_t = Gaudi::Accumulators::Histogram<1>; - using axis_t = hist_t::AccumulatorType::AxisType; + using hist_t = Gaudi::Accumulators::StaticHistogram<1>; + using axis_t = Gaudi::Accumulators::Axis<double>; m_hist1 = std::make_unique<hist_t>( this, "Histo1", "Histogram 1", axis_t{ m_hist1def.value() } ); m_hist2 = std::make_unique<hist_t>( this, "Histo2", "Histogram 2", axis_t( m_hist2def.value() ) ); } ); @@ -73,8 +74,8 @@ namespace Gaudi { Gaudi::Property<Gaudi::Histo1DDef> m_hist2def{ this, "Histo2", { "Histogram2", -5, 5, 200 }, "The parameters for the second histogram" }; - std::unique_ptr<Gaudi::Accumulators::Histogram<1>> m_hist1; - std::unique_ptr<Gaudi::Accumulators::Histogram<1>> m_hist2; + std::unique_ptr<Gaudi::Accumulators::StaticHistogram<1>> m_hist1; + std::unique_ptr<Gaudi::Accumulators::StaticHistogram<1>> m_hist2; }; } // namespace TestSuite } // end of namespace Gaudi diff --git a/GaudiTestSuite/src/testing/HistogramsTests.cpp b/GaudiTestSuite/src/testing/HistogramsTests.cpp index 5da4d13756dfe483b0af512fe76aff1c7a5e9945..62f1bd12c00d7bf29e30259fe6d6d43d474b3d83 100644 --- a/GaudiTestSuite/src/testing/HistogramsTests.cpp +++ b/GaudiTestSuite/src/testing/HistogramsTests.cpp @@ -1,5 +1,5 @@ /***********************************************************************************\ -* (c) Copyright 2021-2023 CERN for the benefit of the LHCb and ATLAS collaborations S* +* (c) Copyright 2021-2024 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". * @@ -8,23 +8,40 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \***********************************************************************************/ -#include <Gaudi/Accumulators/Histogram.h> +#include <Gaudi/Accumulators/StaticHistogram.h> #include <Gaudi/Functional/Consumer.h> #include <GaudiKernel/AlgTool.h> #include <type_traits> namespace Gaudi::Tests::Histograms::CustomAxis { enum class Category { Simple, Complex, Bad, Wrong }; -} + std::ostream& operator<<( std::ostream& o, Category v ) { + switch ( v ) { + case Category::Simple: + o << "Simple"; + break; + case Category::Complex: + o << "Complex"; + break; + case Category::Bad: + o << "Bad"; + break; + case Category::Wrong: + o << "Wrong"; + break; + } + return o; + } +} // namespace Gaudi::Tests::Histograms::CustomAxis namespace Gaudi::Accumulators { - template <> - struct Axis<Gaudi::Tests::Histograms::CustomAxis::Category> { - using Category = Gaudi::Tests::Histograms::CustomAxis::Category; + struct CustomAxis { + using Category = Gaudi::Tests::Histograms::CustomAxis::Category; + using ArithmeticType = Category; - Axis() = default; + CustomAxis() = default; /// number of bins for this Axis - unsigned int nBins = 4; + unsigned int numBins() const { return 4; } /// min and max values on this axis std::underlying_type_t<Category> minValue = 0, maxValue = 4; /// title of this axis @@ -32,13 +49,23 @@ namespace Gaudi::Accumulators { /// labels for the bins std::vector<std::string> labels{ "Simple", "Complex", "Bad", "Wrong" }; - unsigned int index( Category value ) const { return static_cast<unsigned int>( value ) + 1; } + unsigned int index( Category value ) const { return static_cast<unsigned int>( value ) + 1; } + friend std::ostream& operator<<( std::ostream& o, CustomAxis const& axis ) { + return o << axis.numBins() << " " << axis.minValue << " " << axis.maxValue; + } }; + void to_json( nlohmann::json& j, const CustomAxis& axis ) { + j = nlohmann::json{ { "nBins", axis.numBins() }, + { "minValue", axis.minValue }, + { "maxValue", axis.maxValue }, + { "title", axis.title } }; + j["labels"] = axis.labels; + } } // namespace Gaudi::Accumulators namespace Gaudi::Tests::Histograms { namespace Directories { - using MyHist_t = Gaudi::Accumulators::Histogram<1>; + using MyHist_t = Gaudi::Accumulators::StaticHistogram<1>; struct HistoGroupsTool : AlgTool { using AlgTool::AlgTool; @@ -81,7 +108,7 @@ namespace Gaudi::Tests::Histograms { namespace AxesLabels { struct HistWithLabelsAlg : Gaudi::Functional::Consumer<void()> { using Base = Gaudi::Functional::Consumer<void()>; - using MyHist_t = Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, int>; + using MyHist_t = Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, int>; using Base::Base; @@ -99,8 +126,9 @@ namespace Gaudi::Tests::Histograms { using Base = Gaudi::Functional::Consumer<void()>; using Base::Base; - mutable Gaudi::Accumulators::Histogram<1, Gaudi::Accumulators::atomicity::full, Category> m_hist{ - this, "Categories", "", Gaudi::Accumulators::Axis<Category>{} }; + mutable Gaudi::Accumulators::StaticHistogram<1, Gaudi::Accumulators::atomicity::full, Category, + std::tuple<Gaudi::Accumulators::CustomAxis>> + m_hist{ this, "Categories", "", Gaudi::Accumulators::CustomAxis{} }; void operator()() const override { ++m_hist[Category::Simple]; @@ -117,9 +145,10 @@ namespace Gaudi::Tests::Histograms { using Base = Gaudi::Functional::Consumer<void()>; using Base::Base; - mutable Gaudi::Accumulators::Histogram<1> m_h1{ this, "h1", "", { 10, 0, 10 } }; - mutable Gaudi::Accumulators::Histogram<2> m_h2{ this, "h2", "", { 10, 0, 10 }, { 10, 0, 10 } }; - mutable Gaudi::Accumulators::Histogram<3> m_h3{ this, "h3", "", { 10, 0, 10 }, { 10, 0, 10 }, { 10, 0, 10 } }; + mutable Gaudi::Accumulators::StaticHistogram<1> m_h1{ this, "h1", "", { 10, 0, 10 } }; + mutable Gaudi::Accumulators::StaticHistogram<2> m_h2{ this, "h2", "", { 10, 0, 10 }, { 10, 0, 10 } }; + mutable Gaudi::Accumulators::StaticHistogram<3> m_h3{ this, "h3", "", + { 10, 0, 10 }, { 10, 0, 10 }, { 10, 0, 10 } }; void operator()() const override { int value = 0; diff --git a/GaudiTestSuite/tests/pytest/refs/Histograms_py.yaml b/GaudiTestSuite/tests/pytest/refs/Histograms_py.yaml index 7e3066673149c5cb54c2b5784befb0ccafe0de74..460ed6c50adabb130e1aee738838c8608b9e7a01 100644 --- a/GaudiTestSuite/tests/pytest/refs/Histograms_py.yaml +++ b/GaudiTestSuite/tests/pytest/refs/Histograms_py.yaml @@ -17,8 +17,8 @@ stdout: |- RndmGenSvc INFO Using Random engine:HepRndm::Engine<CLHEP::RanluxEngine> SimpleCounterHi...SUCCESS Filling Histograms...... Please be patient ! SimpleRootCount...SUCCESS Filling Histograms...... Please be patient ! - SimpleCounterHi... INFO Number of counters : 1 - SimpleRootCount... INFO Number of counters : 1 + SimpleCounterHi... INFO Number of counters : 48 + SimpleRootCount... INFO Number of counters : 10 ApplicationMgr INFO Application Manager Stopped successfully EventLoopMgr INFO Histograms converted successfully according to request. ApplicationMgr INFO Application Manager Finalized successfully