Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-04-26 09:01:08 

0001 /***********************************************************************************\
0002 * (c) Copyright 1998-2024 CERN for the benefit of the LHCb and ATLAS collaborations *
0003 *                                                                                   *
0004 * This software is distributed under the terms of the Apache version 2 licence,     *
0005 * copied verbatim in the file "LICENSE".                                            *
0006 *                                                                                   *
0007 * In applying this licence, CERN does not waive the privileges and immunities       *
0008 * granted to it by virtue of its status as an Intergovernmental Organization        *
0009 * or submit itself to any jurisdiction.                                             *
0010 \***********************************************************************************/
0011 #pragma once
0012 
0013 #include <Gaudi/Accumulators/StaticHistogram.h>
0014 
0015 #include <type_traits>
0016 
0017 namespace {
0018   template <typename tuple_t>
0019   constexpr auto get_array_from_tuple( tuple_t&& tuple ) {
0020     constexpr auto get_array = []( auto&&... x ) { return std::array{ std::forward<decltype( x )>( x )... }; };
0021     return std::apply( get_array, std::forward<tuple_t>( tuple ) );
0022   }
0023 } // namespace
0024 
0025 namespace Gaudi::Accumulators {
0026 
0027   /// number of items in sums for a given dimension
0028   /// = 1 (nb items) + ND (sums of each dimension) + ND*(ND+1)/2 (square sums)
0029   constexpr unsigned int NSUMS( unsigned int ND ) { return 1 + ND + ND * ( ND + 1 ) / 2; }
0030 
0031   template <typename Arithmetic, atomicity Atomicity, unsigned int ND>
0032   struct SigmasValueHandler {
0033     using InputType  = std::conditional_t<ND == 1, Arithmetic, std::array<Arithmetic, ND>>;
0034     using OutputType = std::array<Arithmetic, NSUMS( ND )>;
0035     struct OutputTypeTS : std::array<std::atomic<Arithmetic>, NSUMS( ND )> {
0036       /// copy constructor from non thread safe type
0037       using std::array<std::atomic<Arithmetic>, NSUMS( ND )>::array;
0038       explicit OutputTypeTS( OutputType const& other ) {
0039         for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { ( *this )[i] = other[i]; }
0040       }
0041       // operator= from non thread safe type
0042       OutputTypeTS& operator=( OutputType const& other ) {
0043         for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { ( *this )[i] = other[i]; }
0044         return *this;
0045       }
0046       // automatic conversion to non thread safe type
0047       operator OutputType() const {
0048         OutputType out;
0049         for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { out[i] = ( *this )[i].load( std::memory_order_relaxed ); }
0050         return out;
0051       }
0052     };
0053     using InternalType = std::conditional_t<Atomicity == atomicity::full, OutputTypeTS, OutputType>;
0054     static constexpr OutputType getValue( const InternalType& v ) noexcept {
0055       // Note automatic conversion will happen
0056       return v;
0057     };
0058     static OutputType exchange( InternalType& v, OutputType newv ) noexcept {
0059       if constexpr ( Atomicity == atomicity::full ) {
0060         OutputType old;
0061         for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { old[i] = v[i].exchange( newv[i] ); }
0062         return old;
0063       } else {
0064         return std::exchange( v, newv );
0065       }
0066     }
0067     static constexpr OutputType DefaultValue() { return InternalType{}; }
0068     static void                 merge( InternalType& a, OutputType b ) noexcept {
0069       if constexpr ( Atomicity == atomicity::full ) {
0070         for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { fetch_add( a[i], b[i] ); }
0071       } else {
0072         for ( unsigned int i = 0; i < ND * ( ND + 3 ); i++ ) a[i] += b[i];
0073       }
0074     }
0075     static void merge( InternalType& a, InputType b ) noexcept {
0076       // prepare values to increment internal values
0077       OutputType diff{};
0078       diff[0] = 1.;
0079       if constexpr ( ND == 1 ) {
0080         // no operator[] for b in this case
0081         diff[1] = b;
0082         diff[2] = b * b;
0083       } else {
0084         for ( unsigned int i = 0; i < ND; i++ ) diff[i + 1] += b[i];
0085         unsigned int n = 1 + ND;
0086         for ( unsigned int i = 0; i < ND; i++ ) {
0087           for ( unsigned int j = i; j < ND; j++ ) {
0088             diff[n] = b[i] * b[j];
0089             n++;
0090           }
0091         }
0092       }
0093       // Now increase original counter
0094       if constexpr ( Atomicity == atomicity::full ) {
0095         for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) { fetch_add( a[i], diff[i] ); }
0096       } else {
0097         for ( unsigned int i = 0; i < NSUMS( ND ); i++ ) a[i] += diff[i];
0098       }
0099     }
0100   };
0101 
0102   template <typename Arithmetic, atomicity Atomicity, unsigned int ND>
0103   struct SigmaNAccumulator
0104       : GenericAccumulator<std::array<Arithmetic, ND>, std::array<Arithmetic, NSUMS( ND )>, Atomicity, Identity,
0105                            Identity, SigmasValueHandler<Arithmetic, Atomicity, ND>> {
0106     std::array<Arithmetic, NSUMS( ND )> const sums() const { return this->value(); }
0107   };
0108 
0109   /// specialization for ND=1 to allow for better syntax
0110   template <typename Arithmetic, atomicity Atomicity>
0111   struct SigmaNAccumulator<Arithmetic, Atomicity, 1>
0112       : GenericAccumulator<Arithmetic, std::array<Arithmetic, 3>, Atomicity, Identity, Identity,
0113                            SigmasValueHandler<Arithmetic, Atomicity, 1>> {
0114     std::array<Arithmetic, 3> const sums() const { return this->value(); }
0115   };
0116 
0117   /**
0118    * Internal Accumulator class dealing with RootHistograming.
0119    * Actually a simple extention on top of RootHistograming with an
0120    * extra SigmaCounter embeded
0121    */
0122   template <atomicity Atomicity, typename Arithmetic, unsigned int ND, typename AxisTupleType>
0123   class RootHistogramingAccumulatorInternal
0124       : public HistogramingAccumulatorInternal<Atomicity, HistoInputType<AxisToArithmetic_t<AxisTupleType>, ND>,
0125                                                unsigned long, IntegralAccumulator, AxisTupleType> {
0126 
0127     using InputType = HistoInputType<AxisToArithmetic_t<AxisTupleType>, ND>;
0128     using Parent =
0129         HistogramingAccumulatorInternal<Atomicity, InputType, unsigned long, IntegralAccumulator, AxisTupleType>;
0130 
0131     template <atomicity, typename, unsigned int, typename>
0132     friend class RootHistogramingAccumulatorInternal;
0133 
0134     static_assert( ND <= 3, "Root on supports histogrmas with dimension <= 3" );
0135 
0136     /**
0137      * Small procyclass allowing operator[] to work as expected on the RootHistogram
0138      * that is to return something having an operator+= updating the histogram properly
0139      */
0140     struct Proxy {
0141       Proxy( RootHistogramingAccumulatorInternal& histo, typename InputType::ValueType& v )
0142           : m_histo( histo ), m_v( v ) {}
0143       void                                 operator++() { m_histo.update( m_v ); }
0144       RootHistogramingAccumulatorInternal& m_histo;
0145       typename InputType::ValueType        m_v;
0146     };
0147 
0148   public:
0149     using Parent::Parent;
0150     friend struct Proxy;
0151 
0152     [[deprecated( "Use `++h1[x]`, `++h2[{x,y}]`, etc. instead." )]] RootHistogramingAccumulatorInternal&
0153     operator+=( typename InputType::ValueType v ) {
0154       update( v );
0155       return *this;
0156     }
0157     void reset() {
0158       m_accumulator.reset();
0159       Parent::reset();
0160     }
0161     template <atomicity ato>
0162     void mergeAndReset( RootHistogramingAccumulatorInternal<ato, Arithmetic, ND, AxisTupleType>& other ) {
0163       m_accumulator.mergeAndReset( other.m_accumulator );
0164       Parent::mergeAndReset( other );
0165     }
0166     [[nodiscard]] auto operator[]( typename InputType::ValueType v ) { return Proxy( *this, v ); }
0167 
0168     /// returns the nbentries, sums and "squared sums" of the inputs
0169     /// Practically we have first the number of entries, then the simple sums of each
0170     /// input dimension followed by all combinasions of product of 2 inputs, in "alphabetical" order,
0171     /// e.g. for ND=3 we have sums of n, x, y, z, x^2, xy, xz, y^2, yz, z^2
0172     auto sums2() const { return m_accumulator.sums(); }
0173 
0174   private:
0175     void update( typename InputType::ValueType v ) {
0176       // Do not accumulate in m_accumulator if we are outside the histo range
0177       // We mimic here the behavior of ROOT
0178       if ( v.inAcceptance( this->axis() ) ) {
0179         if constexpr ( ND == 1 ) {
0180           m_accumulator += std::get<0>( v );
0181         } else {
0182           m_accumulator += get_array_from_tuple( static_cast<typename InputType::InternalType>( v ) );
0183         }
0184       }
0185       ++Parent::operator[]( v );
0186     }
0187     // Accumulator for keeping squared sum of value stored in the histogram and correlation values
0188     // they get stored in "alphabetical" order, e.g. for ND=3 x^2, xy, xz, y^2, yz, z^2
0189     SigmaNAccumulator<Arithmetic, Atomicity, ND> m_accumulator;
0190   };
0191 
0192   namespace {
0193     /// helper function to compute standard_deviation
0194     template <typename Arithmetic>
0195     Arithmetic stddev( Arithmetic n, Arithmetic s, Arithmetic s2 ) {
0196       using Gaudi::Accumulators::sqrt;
0197       using std::sqrt;
0198       auto v = ( n > 0 ) ? ( ( s2 - s * ( s / n ) ) / n ) : Arithmetic{};
0199       return ( Arithmetic{ 0 } > v ) ? Arithmetic{} : sqrt( v );
0200     }
0201   } // namespace
0202 
0203   /**
0204    * Class implementing a root histogram accumulator
0205    */
0206   template <atomicity Atomicity, typename Arithmetic, typename ND, typename AxisTupleType>
0207   struct RootHistogramingAccumulator;
0208 
0209   template <atomicity Atomicity, typename Arithmetic, typename AxisTupleType>
0210   struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<unsigned int, 1>, AxisTupleType>
0211       : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1, AxisTupleType> {
0212     using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1,
0213                                               AxisTupleType>::RootHistogramingAccumulatorInternal;
0214     using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 1, AxisTupleType>::nEntries;
0215     Arithmetic nEntries() const { return this->sums2()[0]; }
0216     Arithmetic sum() const { return this->sums2()[1]; }
0217     Arithmetic sum2() const { return this->sums2()[2]; }
0218     Arithmetic mean() const { return sum() / nEntries(); }
0219     Arithmetic standard_deviation() const { return stddev( nEntries(), sum(), sum2() ); }
0220   };
0221 
0222   template <atomicity Atomicity, typename Arithmetic, typename AxisTupleType>
0223   struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<unsigned int, 2>, AxisTupleType>
0224       : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2, AxisTupleType> {
0225     using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2,
0226                                               AxisTupleType>::RootHistogramingAccumulatorInternal;
0227     using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 2, AxisTupleType>::nEntries;
0228     Arithmetic nEntries() const { return this->sums2()[0]; }
0229     Arithmetic sumx() const { return this->sums2()[1]; }
0230     Arithmetic sumy() const { return this->sums2()[2]; }
0231     Arithmetic sumx2() const { return this->sums2()[3]; }
0232     Arithmetic sumy2() const { return this->sums2()[5]; }
0233     Arithmetic sumxy() const { return this->sums2()[4]; }
0234     Arithmetic meanx() const { return sumx() / nEntries(); }
0235     Arithmetic meany() const { return sumy() / nEntries(); }
0236     Arithmetic standard_deviationx() const { return stddev( nEntries(), sumx(), sumx2() ); }
0237     Arithmetic standard_deviationy() const { return stddev( nEntries(), sumy(), sumy2() ); }
0238   };
0239 
0240   template <atomicity Atomicity, typename Arithmetic, typename AxisTupleType>
0241   struct RootHistogramingAccumulator<Atomicity, Arithmetic, std::integral_constant<unsigned int, 3>, AxisTupleType>
0242       : RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3, AxisTupleType> {
0243     using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3,
0244                                               AxisTupleType>::RootHistogramingAccumulatorInternal;
0245     using RootHistogramingAccumulatorInternal<Atomicity, Arithmetic, 3, AxisTupleType>::nEntries;
0246     Arithmetic nEntries() const { return this->sums2()[0]; }
0247     Arithmetic sumx() const { return this->sums2()[1]; }
0248     Arithmetic sumy() const { return this->sums2()[2]; }
0249     Arithmetic sumz() const { return this->sums2()[3]; }
0250     Arithmetic sumx2() const { return this->sums2()[4]; }
0251     Arithmetic sumy2() const { return this->sums2()[7]; }
0252     Arithmetic sumz2() const { return this->sums2()[9]; }
0253     Arithmetic sumxy() const { return this->sums2()[5]; }
0254     Arithmetic sumxz() const { return this->sums2()[6]; }
0255     Arithmetic sumyz() const { return this->sums2()[8]; }
0256     Arithmetic meanx() const { return sumx() / nEntries(); }
0257     Arithmetic meany() const { return sumy() / nEntries(); }
0258     Arithmetic meanz() const { return sumz() / nEntries(); }
0259     Arithmetic standard_deviationx() const { return stddev( nEntries(), sumx(), sumx2() ); }
0260     Arithmetic standard_deviationy() const { return stddev( nEntries(), sumy(), sumy2() ); }
0261     Arithmetic standard_deviationz() const { return stddev( nEntries(), sumz(), sumz2() ); }
0262   };
0263 
0264   /**
0265    * Extension of the standard Gaudi histogram to provide similar functionnality as ROOT
0266    *
0267    * The main piece of extra functionnality is that ablity to compute a mean and sigma
0268    * of the data set hold by the histogram on top of the histogram itself
0269    * This will then be used when filling a real Root histogram so that one gets the
0270    * same values as if working directly with Root.
0271    * When using pure Gaudi Histograms, and converting to ROOT, ROOT automatically
0272    * recomputed average and sigma form the bins, but the value is not the expected one
0273    *
0274    * Usage is similar to HistogramingCounterBase so see the documentation there
0275    * Serialization has the following extra fields :
0276    *   - nTotEntries, sum, sum2, mean in 1D
0277    *   - nTotEntries, sumx, sumy, sumx2, sumy2, sumxy, meanx, meany in 2D
0278    *   - nTotEntries, sumx, sumy, sumz, sumx2, sumy2, sumz2, sumxy, sumxz, sumyz, meanx, meany, meanz in 3D
0279    * and uses same types as HistogramingCounterBase
0280    *
0281    */
0282   template <unsigned int ND, atomicity Atomicity, typename Arithmetic, const char* Type, typename AxisTupleType>
0283   class RootHistogramingCounterBase;
0284 
0285   template <atomicity Atomicity, typename Arithmetic, const char* Type, typename AxisTupleType>
0286   class RootHistogramingCounterBase<1, Atomicity, Arithmetic, Type, AxisTupleType>
0287       : public HistogramingCounterBase<1, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType> {
0288   public:
0289     using Parent = HistogramingCounterBase<1, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType>;
0290     using Parent::Parent;
0291 
0292     friend void to_json( nlohmann::json&                                                                   j,
0293                          RootHistogramingCounterBase<1, Atomicity, Arithmetic, Type, AxisTupleType> const& h ) {
0294       to_json( j, static_cast<Parent const&>( h ) );
0295       j["nTotEntries"]        = h.nEntries();
0296       j["sum"]                = h.sum();
0297       j["mean"]               = h.mean();
0298       j["sum2"]               = h.sum2();
0299       j["standard_deviation"] = h.standard_deviation();
0300     }
0301   };
0302 
0303   template <atomicity Atomicity, typename Arithmetic, const char* Type, typename AxisTupleType>
0304   class RootHistogramingCounterBase<2, Atomicity, Arithmetic, Type, AxisTupleType>
0305       : public HistogramingCounterBase<2, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType> {
0306   public:
0307     using Parent = HistogramingCounterBase<2, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType>;
0308     using Parent::Parent;
0309 
0310     friend void to_json( nlohmann::json&                                                                   j,
0311                          RootHistogramingCounterBase<2, Atomicity, Arithmetic, Type, AxisTupleType> const& h ) {
0312       to_json( j, static_cast<Parent const&>( h ) );
0313       j["nTotEntries"]         = h.nEntries();
0314       j["sumx"]                = h.sumx();
0315       j["sumy"]                = h.sumy();
0316       j["meanx"]               = h.meanx();
0317       j["meany"]               = h.meany();
0318       j["sumx2"]               = h.sumx2();
0319       j["sumy2"]               = h.sumy2();
0320       j["sumxy"]               = h.sumxy();
0321       j["standard_deviationx"] = h.standard_deviationx();
0322       j["standard_deviationy"] = h.standard_deviationy();
0323     }
0324   };
0325 
0326   template <atomicity Atomicity, typename Arithmetic, const char* Type, typename AxisTupleType>
0327   class RootHistogramingCounterBase<3, Atomicity, Arithmetic, Type, AxisTupleType>
0328       : public HistogramingCounterBase<3, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType> {
0329   public:
0330     using Parent = HistogramingCounterBase<3, Atomicity, Arithmetic, Type, RootHistogramingAccumulator, AxisTupleType>;
0331     using Parent::Parent;
0332 
0333     friend void to_json( nlohmann::json&                                                                   j,
0334                          RootHistogramingCounterBase<3, Atomicity, Arithmetic, Type, AxisTupleType> const& h ) {
0335       to_json( j, static_cast<Parent const&>( h ) );
0336       j["nTotEntries"]         = h.nEntries();
0337       j["sumx"]                = h.sumx();
0338       j["sumy"]                = h.sumy();
0339       j["sumz"]                = h.sumz();
0340       j["meanx"]               = h.meanx();
0341       j["meany"]               = h.meany();
0342       j["meanz"]               = h.meanz();
0343       j["sumx2"]               = h.sumx2();
0344       j["sumy2"]               = h.sumy2();
0345       j["sumz2"]               = h.sumz2();
0346       j["sumxy"]               = h.sumxy();
0347       j["sumxz"]               = h.sumxz();
0348       j["sumyz"]               = h.sumyz();
0349       j["standard_deviationx"] = h.standard_deviationx();
0350       j["standard_deviationy"] = h.standard_deviationy();
0351       j["standard_deviationz"] = h.standard_deviationz();
0352     }
0353   };
0354 
0355   /// Root histograming counter. See RootHistogramingCounterBase for details
0356   template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double,
0357             typename AxisTupleType = make_tuple_t<Axis<Arithmetic>, ND>>
0358   using StaticRootHistogram =
0359       RootHistogramingCounterBase<ND, Atomicity, Arithmetic, naming::histogramString, AxisTupleType>;
0360 
0361 } // namespace Gaudi::Accumulators
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!