Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:57:28 

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