Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-04-19 09:13:40

0001 // -*- C++ -*-
0002 //
0003 // This file is part of YODA -- Yet more Objects for Data Analysis
0004 // Copyright (C) 2008-2024 The YODA collaboration (see AUTHORS for details)
0005 //
0006 #ifndef YODA_Dbn_h
0007 #define YODA_Dbn_h
0008 
0009 #include "YODA/Exceptions.h"
0010 #include "YODA/Utils/DbnUtils.h"
0011 #include "YODA/Utils/MathUtils.h"
0012 #include "YODA/Utils/MetaUtils.h"
0013 #include <cmath>
0014 #include <string>
0015 #include <array>
0016 
0017 namespace YODA {
0018 
0019   /// @brief some utility methods
0020   ///
0021   /// @note Anonymous namespace to limit visibility to this file
0022   namespace {
0023 
0024     /// Helper function that checks if element i is in an (index) sequence
0025     template <size_t ... Is>
0026     constexpr bool has_match(const size_t i, std::integer_sequence<size_t, Is...>) {
0027       return ((i == Is) || ...);
0028     }
0029 
0030     /// base step to concatenate empty sequence with nothing
0031     template<typename INT>
0032     constexpr std::integer_sequence<INT>
0033     concat_sequences(std::integer_sequence<INT>){
0034         return {};
0035     }
0036 
0037     /// base step to concatenate two sequences
0038     template<typename INT, INT... s, INT... t>
0039     constexpr std::integer_sequence<INT,s...,t...>
0040     concat_sequences(std::integer_sequence<INT, s...>, std::integer_sequence<INT, t...>){
0041         return {};
0042     }
0043 
0044     /// recursion step to concatenate multuple sequences
0045     template<typename INT, INT... s, INT... t, class... R>
0046     constexpr auto
0047     concat_sequences(std::integer_sequence<INT, s...>, std::integer_sequence<INT, t...>, R...) {
0048         return concat_sequences(std::integer_sequence<INT,s...,t...>{}, R{}...);
0049     }
0050 
0051     /// @brief Helper method for AxisFilter
0052     /// provides individual axis indices if they pass
0053     template<class INT, INT a, class CHECK>
0054     constexpr auto AxisFilterSingle(std::integer_sequence<INT, a>, CHECK pass) {
0055       if constexpr (pass(a))
0056         return std::integer_sequence<INT, a>{};
0057       else
0058         return std::integer_sequence<INT>{};
0059     }
0060 
0061     /// @brief Helper method to filter axis indices
0062     /// concatenates surviving indices if they pass the check
0063     /// This is a bit of compile-time arithmetic, hence constexpr
0064     template<class INT, INT... b, class CHECK>
0065     constexpr auto AxisFilter(std::integer_sequence<INT, b...>, [[maybe_unused]] CHECK pass) {
0066       if constexpr (sizeof...(b) > 0) // non empty sequence
0067         return concat_sequences(AxisFilterSingle(std::integer_sequence<INT, b>{}, pass)...);
0068       else // empty sequence case
0069         return std::integer_sequence<INT>{};
0070     }
0071 
0072   }
0073 
0074   /// @brief A 1D distribution
0075   ///
0076   /// This class is used internally by YODA to centralise the calculation of
0077   /// statistics of unbounded, unbinned sampled distributions. Each distribution
0078   /// fill contributes a weight, \f$ w \f$, and a value, \f$ x \f$. By storing
0079   /// the total number of fills (ignoring weights), \f$ \sum w \f$, \f$ \sum w^2
0080   /// \f$, \f$ \sum wx \f$, and \f$ \sum wx^2 \f$, the Dbn can calculate the
0081   /// mean and spread (\f$ \sigma^2 \f$, \f$ \sigma \f$ and \f$ \hat{\sigma}
0082   /// \f$) of the sampled distribution. It is used to provide this information
0083   /// in bins and for the "hidden" \f$ y \f$ distribution in profile histogram
0084   /// bins.
0085   template <size_t N>
0086   class DbnBase {
0087   public:
0088 
0089     using BaseT = DbnBase<N>;
0090     using FillDim = std::integral_constant<size_t, N>;
0091     using DataSize = std::integral_constant<size_t, 1 + 2*(N+1) + (N*(N-1)/2)>;
0092 
0093     /// @name Constructors
0094     //@{
0095 
0096     /// Default constructor of a new distribution.
0097     DbnBase() { reset(); }
0098 
0099 
0100     /// @brief Constructor to set a distribution with a pre-filled state.
0101     ///
0102     /// Principally designed for internal persistency use.
0103     /// @note No cross-term possible for 0- and 1D.
0104     template<size_t dim = N, typename = std::enable_if_t<(dim < 2)>>
0105     DbnBase(const double numEntries, const std::array<double,N+1>& sumW, const std::array<double,N+1>& sumW2)
0106       : _numEntries(numEntries), _sumW(sumW), _sumW2(sumW2) { }
0107 
0108     /// @brief Constructor to set a distribution with a pre-filled state.
0109     ///
0110     /// Principally designed for internal persistency use.
0111     /// @note Includes second-order cross-terms.
0112     template<size_t dim = N, typename = std::enable_if_t<(dim >= 2)>>
0113     DbnBase(const double numEntries, const std::array<double,N+1>& sumW, const std::array<double,N+1>& sumW2,
0114                                                                          const std::array<double,N*(N-1)/2>& sumWcross)
0115       : _numEntries(numEntries), _sumW(sumW), _sumW2(sumW2), _sumWcross(sumWcross) { }
0116 
0117     /// Copy constructor
0118     ///
0119     /// Sets all the parameters using the ones provided from an existing DbnBase.
0120     DbnBase(const DbnBase&) = default;
0121 
0122     /// Move constructor
0123     ///
0124     /// Sets all the parameters using the ones provided from an existing DbnBase.
0125     DbnBase(DbnBase&&) = default;
0126 
0127     /// Copy assignment
0128     ///
0129     /// Sets all the parameters using the ones provided from an existing DbnBase.
0130     DbnBase& operator = (const DbnBase&) = default;
0131 
0132     /// Move assignment
0133     ///
0134     /// Sets all the parameters using the ones provided from an existing DbnBase.
0135     DbnBase& operator = (DbnBase&&) = default;
0136 
0137     //@}
0138 
0139 
0140     /// @name Modifiers
0141     //@{
0142 
0143 
0144     /// @brief Contribute a sample at @a val (from an array) with weight @a weight.
0145     ///
0146     /// If you use brace enclosed initializer list and call this method in this way:
0147     /// fill({1,2,3}, 2, 3)
0148     /// and get error from compiler with this note:
0149     /// "no known conversion for argument 1 from
0150     /// β€˜<brace-enclosed initializer list>’ to β€˜std::array<double, 0>’",
0151     /// it likely means that your value dimension doesn't match distribution's.
0152     void fill(const std::array<double, N> vals, const double weight = 1.0, const double fraction = 1.0) {
0153       fill(vals, weight, fraction, std::make_index_sequence<N>{});
0154     }
0155 
0156     /// @brief Templated convenience overload of fill method.
0157     ///
0158     /// Example:
0159     ///
0160     /// fill(3.0, 2.1, 3.4, 3.0) // To fill value (X=3.0,Y=2.1,Z=3.0) with weight 3.0
0161     template <typename... Args>
0162     void fill(Args&&... args) {
0163       // N + 2 - dimension + 2 default parameters
0164       static_assert(sizeof...(args) <= N + 2,
0165                     "Value's dimension doesn't match distribution's dimension.");
0166       std::array<double, sizeof...(args)> vals = {{args...}};
0167       double weight = 1.0;
0168       double fraction = 1.0;
0169       if (vals.size() > N) {
0170         weight = vals[N];
0171         if (vals.size() > N + 1)  fraction = vals[N + 1];
0172       }
0173       const double sf = fraction * weight;
0174       _numEntries += fraction;
0175       _sumW.at(0) += sf;
0176       _sumW2.at(0) += fraction*sqr(weight);
0177       for (unsigned int i = 0; i < N; ++i) {
0178         _sumW.at(i + 1) += sf * vals.at(i);
0179         _sumW2.at(i + 1) += sf * sqr(vals.at(i));
0180       }
0181 
0182       if constexpr(N >= 2) {
0183         size_t idx = 0;
0184         for (size_t i = 0; i < (N-1); ++i) {
0185           for (size_t j = i+1; j < N; ++j) {
0186             _sumWcross.at(idx++) += sf * vals.at(i) * vals.at(j);
0187           }
0188         }
0189       }
0190 
0191     }
0192 
0193 
0194     /// @brief Set a sample with array-based number of entries @a numEntries,
0195     /// sum of weights @a sumW (and corresponding moments),
0196     /// sum of squared weights @a sumW2 (and corresponding moments).
0197     void set(const double numEntries, const std::array<double,N+1>& sumW,
0198                                       const std::array<double,N+1>& sumW2,
0199                                       const std::array<double,N*(N-1)/2>& sumWcross) {
0200       _numEntries = numEntries;
0201       _sumW = sumW;
0202       _sumW2 = sumW2;
0203       _sumWcross = sumWcross;
0204     }
0205 
0206     /// @brief Set a sample with vector-based number of entries @a numEntries,
0207     /// sum of weights @a sumW (and corresponding moments),
0208     /// sum of squared weights @a sumW2 (and corresponding moments).
0209     void set(const double numEntries, const std::vector<double>& sumW,
0210                                       const std::vector<double>& sumW2,
0211                                       const std::vector<double>& sumWcross = {}) {
0212       if (!(sumW.size() <= (N + 1) || sumW2.size() <= (N + 1) || sumWcross.size() <= (N*(N-1)/2)))
0213         throw UserError("Value's dimension doesn't match distribution's dimension.");
0214       _numEntries = numEntries;
0215       std::copy_n(std::make_move_iterator(sumW.begin()),      sumW.size(),      _sumW.begin());
0216       std::copy_n(std::make_move_iterator(sumW2.begin()),     sumW2.size(),     _sumW2.begin());
0217       std::copy_n(std::make_move_iterator(sumWcross.begin()), sumWcross.size(), _sumWcross.begin());
0218     }
0219 
0220 
0221     /// @brief Set a sample with @a other
0222     void set(const DbnBase<N>& other) {
0223       if (this != &other)  *this = other;
0224     }
0225 
0226     /// @brief Set a sample with @a other
0227     void set(DbnBase<N>&& other) {
0228       if (this != &other)  *this = std::move(other);
0229     }
0230 
0231     /// Reset the internal counters.
0232     void reset() {
0233       _numEntries= 0;
0234       _sumW.fill(0);
0235       _sumW2.fill(0);
0236       _sumWcross.fill(0);
0237     }
0238 
0239 
0240     /// Rescale as if all fill weights had been different by factor @a scalefactor.
0241     void scaleW(const double scalefactor) {
0242       _sumW.at(0) *= scalefactor;
0243       _sumW2.at(0) *= sqr(scalefactor);
0244       for (size_t i = 0; i< N; ++i) {
0245         // first- and second-order moments
0246         _sumW.at(i+1) *= scalefactor;
0247         _sumW2.at(i+1) *= scalefactor;
0248       }
0249       for (size_t i = 0; i < _sumWcross.size(); ++i) {
0250         _sumWcross.at(i) *= scalefactor;
0251       }
0252     }
0253 
0254     void scale(const size_t dim, const double factor) {
0255       if (dim >= N)
0256         throw RangeError("Dimension index should be less than "+std::to_string(N));
0257       _sumW.at(dim+1) *= factor;
0258       _sumW2.at(dim+1) *= sqr(factor);
0259       size_t idx = 0;
0260       for (size_t i = 0; i < (N-1); ++i) {
0261         for (size_t j = i+1; j < N; ++j) {
0262           if (i == dim || j == dim) {
0263             _sumWcross.at(idx++) *= factor;
0264           }
0265         }
0266       } // end of double loop
0267     }
0268 
0269 
0270     //@}
0271 
0272     public:
0273 
0274 
0275     /// @name Distribution statistics
0276     //@{
0277 
0278     /// The absolute error on sumW
0279     double errW() const;
0280 
0281     /// The relative error on sumW
0282     double relErrW() const;
0283 
0284     /// Weighted mean, \f$ \bar{x} \f$, of distribution.
0285     double mean(const size_t i) const;
0286 
0287     /// Weighted variance, \f$ \sigma^2 \f$, of distribution.
0288     double variance(const size_t i) const;
0289 
0290     /// Weighted standard deviation, \f$ \sigma \f$, of distribution.
0291     double stdDev(const size_t i) const { return std::sqrt(variance(i)); }
0292 
0293     /// Weighted standard error on the mean, \f$ \sim \sigma/\sqrt{N-1} \f$, of distribution.
0294     double stdErr(const size_t i) const;
0295 
0296     /// Relative weighted standard error on the mean, \f$ \sim \sigma/\sqrt{N-1} \f$, of distribution.
0297     double relStdErr(const size_t i) const;
0298 
0299     /// Weighted RMS, \f$ \sqrt{ \sum{w x^2}/\sum{w} } \f$, of distribution.
0300     double RMS(const size_t i) const;
0301 
0302     //@}
0303 
0304 
0305     /// @name Raw distribution running sums
0306     //@{
0307 
0308     /// Number of entries (number of times @c fill was called, ignoring weights)
0309     double numEntries() const {
0310       return _numEntries;
0311     }
0312 
0313     /// Effective number of entries \f$ = (\sum w)^2 / \sum w^2 \f$
0314     double effNumEntries() const {
0315       if (_sumW2.at(0) == 0) return 0;
0316       return _sumW.at(0)*_sumW.at(0) / _sumW2.at(0);
0317     }
0318 
0319     /// The sum of weights
0320     double sumW() const {
0321       return _sumW.at(0);
0322     }
0323 
0324     /// The sum of weights squared
0325     double sumW2() const {
0326       return _sumW2.at(0);
0327     }
0328 
0329     /// The sum of x*weight
0330     double sumW(const size_t i) const {
0331       return _sumW.at(i);
0332     }
0333 
0334     /// The sum of x^2*weight
0335     double sumW2(const size_t i) const {
0336       return _sumW2.at(i);
0337     }
0338 
0339     /// The second-order cross term between axes @a k and @a l
0340     template<size_t dim = N, typename = std::enable_if_t<(dim >= 2)>>
0341     double crossTerm(const size_t A1, const size_t A2) const {
0342       if (A1 >= N || A2 >= N)  throw RangeError("Invalid axis int, must be in range 0..dim-1");
0343       if (A1 >= A2)  throw RangeError("Indices need to be different for cross term");
0344 
0345       size_t idx = 0;
0346       for (size_t i = 0; i < (N-1); ++i) {
0347         for (size_t j = i+1; j < N; ++j) {
0348           if (i == A1 && j == A2)  break;
0349           ++idx;
0350         }
0351         if (i == A1)  break;
0352       }
0353 
0354       return _sumWcross.at(idx);
0355     }
0356 
0357     size_t dim() const { return N; }
0358 
0359     //@}
0360 
0361     /// @name I/O
0362     /// @{
0363 
0364     /// @brief String representation of the DbnBase for debugging
0365     std::string toString() const {
0366       std::string res ="";
0367       res += ("numEntries="+ std::to_string(_numEntries))  ;
0368       if (sumW()) {
0369         res += (", Mean="+ std::to_string(mean(0)))  ;
0370         res += (", RMS="+ std::to_string(RMS(0)))  ;
0371       }
0372       return res;
0373     }
0374 
0375     //@}
0376 
0377     /// @name Reduce operations
0378     // @{
0379 
0380     /// @brief Reduce operation that produces a lower-dimensional DbnBase
0381     /// keeping only the specified axes in the new DbnBase.
0382     template<typename... size_t>
0383     DbnBase<sizeof...(size_t)> reduceTo(size_t... axes) const {
0384       if constexpr (N == sizeof...(axes)) { return *this; } // keep all axes
0385       else {
0386         // only select the axis moments to keep (0 is sumW, 1...N are the axis moments)
0387         std::array<double, sizeof...(axes)+1> newSumW = { _sumW[0], _sumW[axes+1]... };
0388         std::array<double, sizeof...(axes)+1> newSumW2 = { _sumW2[0], _sumW2[axes+1]... };
0389         if constexpr (sizeof...(axes) < 2) {
0390           return DbnBase<sizeof...(axes)>(_numEntries, newSumW, newSumW2);
0391         }
0392         else {
0393           constexpr auto newDim = sizeof...(axes);
0394           std::array<double, newDim*(newDim-1)/2> newSumWcross;
0395           unsigned int idx = 0;
0396           for (unsigned int i : {axes...}) {
0397             for (unsigned int j : {axes...}) {
0398               if (i < j) { newSumWcross.at(idx) = crossTerm(i, j); ++idx; }
0399             }
0400           }
0401           return DbnBase<sizeof...(axes)>(_numEntries, newSumW, newSumW2, newSumWcross);
0402         }
0403       }
0404     }
0405 
0406     /// @brief Same as above but using an std::integer_sequence
0407     template<size_t... Is>
0408     DbnBase<sizeof...(Is)> reduceTo(std::index_sequence<Is...>) const {
0409       return reduceTo(Is...);
0410     }
0411 
0412     /// @brief Reduce operation that produces a lower-dimensional DbnBase
0413     /// removing the specified axes for the new DbnBase.
0414     template<size_t... axes>
0415     auto reduce() const {
0416       constexpr auto new_axes = AxisFilter(std::make_index_sequence<N>{},
0417         [](size_t i) { return !has_match(i, std::integer_sequence<size_t, axes...>{}); }
0418       );
0419       return reduceTo(new_axes);
0420     }
0421 
0422     // @}
0423 
0424     /// @name MPI (de-)serialisation
0425     //@{
0426 
0427     size_t _lengthContent() const noexcept {
0428       return DataSize::value;
0429     }
0430 
0431     std::vector<double> _serializeContent() const noexcept {
0432       std::vector<double> rtn;
0433       rtn.reserve(DataSize::value);
0434       rtn.insert(rtn.end(), _sumW.begin(),  _sumW.end());
0435       rtn.insert(rtn.end(), _sumW2.begin(), _sumW2.end());
0436       rtn.insert(rtn.end(), _sumWcross.begin(), _sumWcross.end());
0437       rtn.push_back(_numEntries);
0438       return rtn;
0439     }
0440 
0441     void _deserializeContent(const std::vector<double>& data) {
0442 
0443       if (data.size() != DataSize::value)
0444         throw UserError("Length of serialized data should be "+std::to_string(DataSize::value)+"!");
0445 
0446       auto itr = data.cbegin();
0447       std::copy_n(itr, _sumW.size(), _sumW.begin());
0448       std::copy_n(itr+N+1, _sumW2.size(), _sumW2.begin());
0449       std::copy_n(itr+2*(N+1), _sumWcross.size(), _sumWcross.begin());
0450       _numEntries = *(itr + 2*(N+1) + (N*(N-1)/2));
0451 
0452     }
0453 
0454     // @}
0455 
0456     /// @name Operators
0457     //@{
0458 
0459     /// Add two DbnBases
0460     DbnBase& operator += (const DbnBase& d);
0461 
0462     /// Add two DbnBases
0463     DbnBase& operator += (DbnBase&& d);
0464 
0465     /// Subtract one DbnBase from another
0466     DbnBase& operator -= (const DbnBase& d);
0467 
0468     /// Subtract one DbnBase from another
0469     DbnBase& operator -= (DbnBase&& d);
0470 
0471     //@}
0472 
0473     private:
0474 
0475     /// @name Storage
0476     //@{
0477 
0478     double _numEntries;
0479 
0480     /// The 1st order weighted x moment
0481     std::array<double, N+1> _sumW; //indexed by dimension
0482 
0483     /// The 2nd order weighted x moment
0484     std::array<double, N+1> _sumW2; //indexed by dimension
0485 
0486     /// The 2nd order weighted "cross-term"
0487     std::array<double, N*(N-1)/2> _sumWcross;
0488 
0489     //@}
0490 
0491     /// @name Helper methods
0492     //@{
0493 
0494     /// @brief Expands values array
0495     template <size_t... I>
0496     void fill(const std::array<double, N> vals, const double weight, const double fraction,
0497               std::index_sequence<I...>) {
0498       fill(vals[I]..., weight, fraction);
0499     }
0500 
0501     //@}
0502   };
0503 
0504 
0505   template <size_t N>
0506   double DbnBase<N>::errW() const {
0507     return sqrt(sumW2(0));
0508   }
0509 
0510   template <size_t N>
0511   double DbnBase<N>::relErrW() const {
0512     if (effNumEntries() == 0) {
0513       return std::numeric_limits<double>::quiet_NaN();
0514     }
0515     return errW()/sumW(0);
0516   }
0517 
0518   template <size_t N>
0519   double DbnBase<N>::mean(const size_t i) const {
0520     return YODA::mean(sumW(i), sumW(0));
0521   }
0522 
0523 
0524   template <size_t N>
0525   double DbnBase<N>::variance(const size_t i) const {
0526     return YODA::variance(sumW(i), sumW(0), sumW2(i), sumW2(0));
0527   }
0528 
0529 
0530   template <size_t N>
0531   double DbnBase<N>::stdErr(const size_t i) const {
0532     return YODA::stdErr(sumW(i), sumW(0), sumW2(i), sumW2(0));
0533   }
0534 
0535 
0536   template <size_t N>
0537   double DbnBase<N>::relStdErr(const size_t i) const {
0538     if (effNumEntries() == 0) {
0539       return std::numeric_limits<double>::quiet_NaN();
0540     }
0541     return stdErr(i) / mean(i);
0542   }
0543 
0544 
0545   template <size_t N>
0546   double DbnBase<N>::RMS(const size_t i) const {
0547     return YODA::RMS(sumW2(i), sumW(0), sumW2());
0548   }
0549 
0550   template <size_t N>
0551   DbnBase<N>& DbnBase<N>::operator += (const DbnBase<N>& d) {
0552     _numEntries += d._numEntries;
0553     for (size_t i = 0; i <= N; ++i) {
0554       _sumW.at(i) += d._sumW.at(i);
0555       _sumW2.at(i) += d._sumW2.at(i);
0556     }
0557     for (size_t i = 0; i < _sumWcross.size(); ++i) {
0558       _sumWcross.at(i) += d._sumWcross.at(i);
0559     }
0560     return *this;
0561   }
0562 
0563 
0564   template <size_t N>
0565   DbnBase<N>& DbnBase<N>::operator += (DbnBase<N>&& d) {
0566     _numEntries += std::move(d._numEntries);
0567     for (size_t i = 0; i <= N; ++i) {
0568       _sumW.at(i) += std::move(d._sumW.at(i));
0569       _sumW2.at(i) += std::move(d._sumW2.at(i));
0570     }
0571     for (size_t i = 0; i < _sumWcross.size(); ++i) {
0572       _sumWcross.at(i) += std::move(d._sumWcross.at(i));
0573     }
0574     return *this;
0575   }
0576 
0577 
0578   template <size_t N>
0579   DbnBase<N>& DbnBase<N>::operator -= (const DbnBase<N>& d) {
0580     _numEntries -= d._numEntries;
0581     for (unsigned int i =0; i<= N ; ++i) {
0582       _sumW.at(i)  -= d._sumW.at(i);
0583       _sumW2.at(i) += d._sumW2.at(i);
0584     }
0585     for (size_t i = 0; i < _sumWcross.size(); ++i) {
0586       _sumWcross.at(i) -= d._sumWcross.at(i);
0587     }
0588     return *this;
0589   }
0590 
0591 
0592   template <size_t N>
0593   DbnBase<N>& DbnBase<N>::operator -= (DbnBase<N>&& d) {
0594     _numEntries -= std::move(d._numEntries);
0595     for (unsigned int i =0; i<= N ; ++i) {
0596       _sumW.at(i)  -= std::move(d._sumW.at(i));
0597       _sumW2.at(i) += std::move(d._sumW2.at(i));
0598     }
0599     for (size_t i = 0; i < _sumWcross.size(); ++i) {
0600       _sumWcross.at(i) -= std::move(d._sumWcross.at(i));
0601     }
0602     return *this;
0603   }
0604 
0605 
0606   /// Add two DbnBases
0607   template <size_t N>
0608   inline DbnBase<N> operator + (DbnBase<N> first, const DbnBase<N>& second) {
0609     first += second;
0610     return first;
0611   }
0612 
0613   /// Add two DbnBases
0614   template <size_t N>
0615   inline DbnBase<N> operator + (DbnBase<N> first, DbnBase<N>&& second) {
0616     first += std::move(second);
0617     return first;
0618   }
0619 
0620   /// Subtract one DbnBase from another one
0621   template <size_t N>
0622   inline DbnBase<N> operator - (DbnBase<N> first, const DbnBase<N>& second) {
0623     first -= second;
0624     return first;
0625   }
0626 
0627   /// Subtract one DbnBase from another one
0628   template <size_t N>
0629   inline DbnBase<N> operator - (DbnBase<N> first, DbnBase<N>&& second) {
0630     first -= std::move(second);
0631     return first;
0632   }
0633 
0634 
0635   /// @brief User-facing Dbn class inheriting from DbnBase
0636   template<size_t N>
0637   class Dbn : public DbnBase<N> {
0638     public:
0639     using BaseT = DbnBase<N>;
0640     using BaseT::BaseT;
0641     using BaseT::operator=;
0642 
0643   };
0644 
0645   /// @brief Partial template specialisation for Dbn0D
0646   template<>
0647   class Dbn<0> : public DbnBase<0> {
0648     public:
0649 
0650     using BaseT = DbnBase<0>;
0651     using BaseT::BaseT;
0652     using BaseT::operator=;
0653 
0654     Dbn() : BaseT() {}
0655 
0656     Dbn(const double numEntries, const double sumW, const double sumW2)
0657     : BaseT(numEntries, {sumW} , {sumW2}) { }
0658 
0659     Dbn(const BaseT& other) : BaseT(other) {}
0660 
0661     Dbn(BaseT&& other) : BaseT(std::move(other)) {}
0662 
0663     void fill(const std::array<double, 0>& vals, const double weight = 1.0, const double fraction = 1.0) {
0664       BaseT::fill(vals, weight, fraction);
0665     }
0666 
0667     void fill(const double weight=1.0, const double fraction=1.0) { BaseT::fill({}, weight, fraction); }
0668 
0669     void set(const double numEntries, const double sumW, const double sumW2) {
0670       BaseT::set(numEntries, {sumW}, {sumW2});
0671     }
0672 
0673   };
0674 
0675 
0676   /// @brief Partial template specialisation with CRTP for x-methods
0677   template<>
0678   class Dbn<1> : public DbnBase<1>,
0679                  public XDbnMixin<Dbn<1>> {
0680     public:
0681 
0682     using BaseT = DbnBase<1>;
0683     using BaseT::BaseT;
0684     using BaseT::operator=;
0685 
0686     Dbn() : BaseT() {}
0687 
0688     Dbn(const double numEntries, const double sumW, const double sumW2, const double sumWX, const double sumWX2)
0689     : BaseT(numEntries, {sumW, sumWX} , {sumW2, sumWX2}) { }
0690 
0691     Dbn(const BaseT& other) : BaseT(other) {}
0692 
0693     Dbn(BaseT&& other) : BaseT(std::move(other)) {}
0694 
0695     void fill(const std::array<double, 1>& vals, const double weight = 1.0, const double fraction = 1.0) {
0696       BaseT::fill(vals, weight, fraction);
0697     }
0698 
0699     void fill(const double val, const double weight=1.0, const double fraction=1.0) {
0700       BaseT::fill({val}, weight, fraction);
0701     }
0702 
0703   };
0704 
0705 
0706   /// @brief Partial template specialisation with CRTP for x- and y-methods
0707   template<>
0708   class Dbn<2> : public DbnBase<2>,
0709                  public XDbnMixin<Dbn<2>>,
0710                  public YDbnMixin<Dbn<2>> {
0711     public:
0712 
0713     using BaseT = DbnBase<2>;
0714     using BaseT::BaseT;
0715     using BaseT::operator=;
0716 
0717     Dbn() : BaseT() {}
0718 
0719     Dbn(const double numEntries, const double sumW, const double sumW2, const double sumWX, const double sumWX2,
0720                                  const double sumWY, const double sumWY2, const double sumWXY)
0721     : BaseT(numEntries, {sumW, sumWX, sumWY} , {sumW2, sumWX2, sumWY2}, {sumWXY}) { }
0722 
0723     Dbn(const BaseT& other) : BaseT(other) {}
0724 
0725     Dbn(BaseT&& other) : BaseT(std::move(other)) {}
0726 
0727     void fill(const std::array<double, 2>& vals, const double weight = 1.0, const double fraction = 1.0) {
0728       BaseT::fill(vals, weight, fraction);
0729     }
0730 
0731     void fill(const double valX, const double valY, const double weight=1.0, const double fraction=1.0) {
0732       BaseT::fill({valX, valY}, weight, fraction);
0733     }
0734 
0735 
0736   };
0737 
0738 
0739   /// @brief Partial template specialisation with CRTP for x-, y- and z-methods
0740   template<>
0741   class Dbn<3> : public DbnBase<3>,
0742                  public XDbnMixin<Dbn<3>>,
0743                  public YDbnMixin<Dbn<3>>,
0744                  public ZDbnMixin<Dbn<3>> {
0745     public:
0746 
0747     using BaseT = DbnBase<3>;
0748     using BaseT::BaseT;
0749     using BaseT::operator=;
0750 
0751     Dbn() : BaseT() {}
0752 
0753     Dbn(const double numEntries, const double sumW, const double sumW2,
0754                                  const double sumWX, const double sumWX2,
0755                                  const double sumWY, const double sumWY2,
0756                                  const double sumWZ, const double sumWZ2,
0757                                  const double sumWXY, const double sumWXZ, const double sumWYZ)
0758     : BaseT(numEntries, {sumW, sumWX, sumWY, sumWZ} , {sumW2, sumWX2, sumWY2, sumWZ2}, {sumWXY, sumWXZ, sumWYZ}) { }
0759 
0760     Dbn(const BaseT& other) : BaseT(other) {}
0761 
0762     Dbn(BaseT&& other) : BaseT(std::move(other)) {}
0763 
0764     void fill(const std::array<double, 3>& vals, const double weight = 1.0, const double fraction = 1.0) {
0765       BaseT::fill(vals, weight, fraction);
0766     }
0767 
0768     void fill(const double valX, const double valY, const double valZ, const double weight=1.0, const double fraction=1.0) {
0769       BaseT::fill({valX, valY, valZ}, weight, fraction);
0770     }
0771 
0772   };
0773 
0774   /// User-friendly aliases
0775   using Dbn0D = Dbn<0>;
0776   using Dbn1D = Dbn<1>;
0777   using Dbn2D = Dbn<2>;
0778   using Dbn3D = Dbn<3>;
0779 }
0780 
0781 #endif