include/YODA/BinnedEstimate.h

0001 // -*- C++ -*-
0002 //
0003 // This file is part of YODA -- Yet more Objects for Data Analysis
0004 // Copyright (C) 2008-2025 The YODA collaboration (see AUTHORS for details)
0005 //
0006 #ifndef YODA_BinnedEstimate_h
0007 #define YODA_BinnedEstimate_h
0008
0009 #include "YODA/AnalysisObject.h"
0010 #include "YODA/BinnedStorage.h"
0011 #include "YODA/Estimate.h"
0012 #include "YODA/Scatter.h"
0013 #include "YODA/Transformation.h"
0014 #include "YODA/Utils/BinnedUtils.h"
0015 #include <regex>
0016
0017 #ifdef HAVE_HDF5
0018 #include "YODA/Utils/H5Utils.h"
0019 #endif
0020
0021 namespace YODA {
0022
0023
0024   /// @note Anonymous namespace to limit visibility to this file
0025   namespace {
0026     /// @brief Nullify bin uncertainty for discrete axes
0027     template<size_t axisN, typename BinT>
0028     std::pair<double, double> nullifyIfDisc(const BinT& /* b */, const double /* val */, std::false_type, const double null = 0.0) {
0029       return { null, null };
0030     }
0031
0032     /// @brief Work out bin uncertainty with respect to
0033     /// central value for continuous axes
0034     template<size_t axisN, typename BinT>
0035     std::pair<double, double> nullifyIfDisc(const BinT& b, const double val, std::true_type, const double /* null */ = 0.0) {
0036       return { val - b.template min<axisN>(), b.template max<axisN>() - val };
0037     }
0038
0039     /// @brief Unify bin interval for discrete axes
0040     template<size_t axisN, typename BinT>
0041     double unifyIfDisc(const BinT& /* b */, std::false_type, const double unity = 1.0) {
0042       return unity;
0043     }
0044
0045     /// @brief Get bin width for continuous axes
0046     template<size_t axisN, typename BinT>
0047     double unifyIfDisc(const BinT& b, std::true_type, const double /* unity */ = 1.0) {
0048       return b.template width<axisN>();
0049     }
0050
0051     /// @brief Return bin index for non-integral discrete axes
0052     template<size_t axisN, typename BinT>
0053     double coordPicker(const BinT& b, std::false_type, std::false_type) {
0054       return b.index();
0055     }
0056
0057     /// @brief Return bin edge for integral discrete axes
0058     template<size_t axisN, typename BinT>
0059     double coordPicker(const BinT& b, std::true_type, std::false_type) {
0060       return b.template edge<axisN>();
0061     }
0062     /// @brief Return bin mid point for continuous axes
0063     template<size_t axisN, typename BinT>
0064     double coordPicker(const BinT& b, std::true_type, std::true_type) {
0065       return b.template mid<axisN>();
0066     }
0067   }
0068
0069
0070   /// Forward declaration
0071   template <typename... AxisT>
0072   class BinnedEstimate;
0073
0074   /// @brief EstimateStorage convenience class based on BinnedStorage.
0075   ///
0076   /// @note We use this abstraction layer to implement the bulk once and only once.
0077   /// The user-facing BinnedEstimate type will inherit all its methods from this
0078   /// base class along with a few axis-specifc mixins.
0079   template<typename... AxisT>
0080   class EstimateStorage
0081       : public BinnedStorage<Estimate, AxisT...>,
0082         public AnalysisObject {
0083   public:
0084
0085     using BaseT = BinnedStorage<Estimate, AxisT...>;
0086     using BinningT = typename BaseT::BinningT;
0087     using BinT = typename BaseT::BinT;
0088     using BinType = BinT;
0089     using AnalysisObject::operator =;
0090
0091     /// @name Constructors
0092     /// @{
0093
0094     /// @brief Nullary constructor for unique pointers etc.
0095     EstimateStorage(const std::string& path = "", const std::string& title = "")
0096         : BaseT(), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0097
0098     /// @brief Constructor giving explicit bin edges as rvalue reference.
0099     ///
0100     /// Discrete axes have as many edges as bins.
0101     /// Continuous axes have bins.size()+1 edges, the last one
0102     /// being the upper bound of the last bin.
0103     EstimateStorage(std::vector<AxisT>&&... binsEdges, const std::string& path = "", const std::string& title = "")
0104         : BaseT(Axis<AxisT>(std::move(binsEdges))...), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0105
0106     /// @brief Constructor giving explicit bin edges as lvalue reference.
0107     ///
0108     /// Discrete axes have as many edges as bins.
0109     /// Continuous axes have bins.size()+1 edges, the last one
0110     /// being the upper bound of the last bin.
0111     EstimateStorage(const std::vector<AxisT>&... binsEdges, const std::string& path = "", const std::string& title = "")
0112         : BaseT(Axis<AxisT>(binsEdges)...), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0113
0114     /// @brief Constructor giving explicit bin edges as initializer list.
0115     ///
0116     /// Discrete axes have as many edges as bins.
0117     /// Continuous axes have bins.size()+1 edges, the last one
0118     /// being the upper bound of the last bin.
0119     EstimateStorage(std::initializer_list<AxisT>&&... binsEdges, const std::string& path = "", const std::string& title = "")
0120         : BaseT(Axis<AxisT>(std::move(binsEdges))...), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0121
0122     /// @brief Constructor giving range and number of bins.
0123     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT...>>
0124     EstimateStorage(const std::vector<size_t>& nBins,
0125                     const std::vector<std::pair<EdgeT, EdgeT>>& limitsLowUp,
0126                     const std::string& path = "", const std::string& title = "")
0127         //: BaseT( BinningT(nBins, limitsLowUp, std::make_index_sequence<sizeof...(AxisT)>{}) ),
0128         : BaseT( _mkBinning(nBins, limitsLowUp, std::make_index_sequence<sizeof...(AxisT)>{}) ),
0129           AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0130
0131     /// @brief Constructor given a binning type
0132     EstimateStorage(const BinningT& binning, const std::string& path = "", const std::string& title = "")
0133         : BaseT(binning), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0134
0135     /// @brief Constructor given an rvalue BinningT
0136     EstimateStorage(BinningT&& binning, const std::string& path = "", const std::string& title = "")
0137         : BaseT(std::move(binning)), AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0138
0139     /// @brief Constructor given a scatter
0140     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT...>>
0141     EstimateStorage(const ScatterND<sizeof...(AxisT)+1>& s, const std::string& path = "", const std::string& title = "")
0142          : BaseT(_mkBinning(s, std::make_index_sequence<sizeof...(AxisT)>{})),
0143            AnalysisObject(mkTypeString<AxisT...>(), path, title) { }
0144
0145     /// @brief Copy constructor
0146     ///
0147     /// @todo Also allow title setting from the constructor?
0148     EstimateStorage(const EstimateStorage& other, const std::string& path = "") : BaseT(other),
0149         AnalysisObject(mkTypeString<AxisT...>(), path!=""? path : other.path(), other, other.title()) { }
0150
0151     /// @brief Move constructor
0152     ///
0153     /// @todo Also allow title setting from the constructor?
0154     EstimateStorage(EstimateStorage&& other, const std::string& path = "") : BaseT(std::move(other)),
0155         AnalysisObject(mkTypeString<AxisT...>(), path!=""? path : other.path(), other, other.title()) { }
0156
0157     /// @brief Make a copy on the stack
0158     EstimateStorage clone() const noexcept {
0159       return EstimateStorage(*this);
0160     }
0161
0162     /// @brief Make a copy on the heap
0163     EstimateStorage* newclone() const noexcept {
0164       return new EstimateStorage(*this);
0165     }
0166
0167     /// @}
0168
0169
0170     /// @name I/O
0171     /// @{
0172
0173   private:
0174
0175     /// @brief Render information about this AO (private implementation)
0176     template<size_t... Is>
0177     void _renderYODA_aux(std::ostream& os, const int width, std::index_sequence<Is...>) const noexcept {
0178
0179       // print bin edges
0180       BaseT::_binning._renderYODA(os);
0181
0182       // Assemble union of error sources, as it's not guaranteed
0183       // that every bin has the same error breakdown
0184       const std::vector<std::string> labels = this->sources();
0185       if (labels.size()) {
0186         os << "ErrorLabels: [";
0187         for (size_t i = 0; i < labels.size(); ++i) {
0188           const std::string& src = labels[i];
0189           if (i)  os << ", ";
0190           os << std::quoted(src);
0191         }
0192         os << "]\n";
0193       }
0194
0195       // column header: content types
0196       os << std::setw(width) << std::left << "# value" << "\t";
0197       const int errwidth = std::max(int(std::to_string(labels.size()).size()+7), width); // "errDn(" + src + ")"
0198       for (size_t i = 0; i < labels.size(); ++i) {
0199         const std::string& src = labels[i];
0200         if (src.empty()) {
0201           os << std::setw(errwidth) << std::left << "totalDn" << "\t"
0202              << std::setw(errwidth) << std::left << "totalUp" << "\t";
0203         }
0204         else {
0205           os << std::setw(errwidth) << std::left << ("errDn(" + std::to_string(i+1) + ")") << "\t"
0206              << std::setw(errwidth) << std::left << ("errUp(" + std::to_string(i+1) + ")") << "\t";
0207         }
0208       }
0209       os << "\n";
0210
0211       for (const auto& b : BaseT::bins(true, true)) {
0212         os << std::setw(width) << std::left << b.val() << "\t"; // print value
0213         // print systs if available
0214         for (const std::string& src : labels) {
0215           if (!b.hasSource(src)) {
0216             os << std::setw(errwidth) << std::left << "---" << "\t"
0217                << std::setw(errwidth) << std::left << "---" << "\t";
0218             continue;
0219           }
0220           const auto& err = b.err(src);
0221           os << std::setw(errwidth) << std::left << err.first << "\t"
0222              << std::setw(errwidth) << std::left << err.second << "\t";
0223         }
0224         os << "\n";
0225       }
0226     }
0227
0228   public:
0229
0230     /// @brief Render information about this AO
0231     void _renderYODA(std::ostream& os, const int width = 13) const noexcept {
0232       _renderYODA_aux(os, width, std::make_index_sequence<sizeof...(AxisT)>{});
0233     }
0234
0235     /// @brief Render scatter-like information about this AO
0236     void _renderFLAT(std::ostream& os, const int width = 13) const noexcept {
0237       const ScatterND<sizeof...(AxisT)+1> tmp = mkScatter();
0238       tmp._renderYODA(os, width);
0239     }
0240
0241     #ifdef HAVE_HDF5
0242
0243     /// @brief Extract error labels of this AO
0244     void _extractLabels(std::vector<std::string>& labels,
0245                         std::vector<size_t>& labelsizes) const noexcept {
0246       size_t nlabels = 0;
0247       for (const auto& b : BaseT::bins(true, true)) {
0248         std::vector<std::string> keys = b.sources();
0249         nlabels += keys.size()+1; //< +1 for length info itself
0250         labels.insert(std::end(labels),
0251                       std::make_move_iterator(std::begin(keys)),
0252                       std::make_move_iterator(std::end(keys)));
0253       }
0254       std::sort(labels.begin(), labels.end());
0255       labels.erase( std::unique(labels.begin(), labels.end()), labels.end() );
0256       labelsizes.emplace_back(std::move(nlabels));
0257     }
0258
0259     /// @brief Extract axes edges of this AO into map of @a edges
0260     void _extractEdges(std::map<std::string, EdgeHandlerBasePtr>& edges,
0261                        const std::vector<std::string>& labels) const noexcept {
0262
0263       using lenT = EdgeHandler<size_t>;
0264       using lenPtr = EdgeHandlerPtr<size_t>;
0265       const std::string lengthID("sizeinfo");
0266       lenPtr nedges = std::static_pointer_cast<lenT>(edges.find(lengthID)->second);
0267
0268       auto extractEdges = [&edges, &binning = BaseT::_binning, &nedges](auto I) {
0269
0270         using thisEdgeT = typename BinningT::template getEdgeT<I>;
0271         using thisHandlerT = EdgeHandler<thisEdgeT>;
0272         using thisHandlerPtr = EdgeHandlerPtr<thisEdgeT>;
0273
0274         const std::string edgeID = std::string("edges_") + TypeID<thisEdgeT>::name();
0275         std::vector<thisEdgeT> tmp = binning.template edges<I>();
0276         nedges->extend({ tmp.size() });
0277
0278         auto itr = edges.find(edgeID);
0279         if (itr != edges.cend()) {
0280           thisHandlerPtr edgeset = std::static_pointer_cast<thisHandlerT>(itr->second);
0281           edgeset->extend(std::move(tmp));
0282         }
0283         else {
0284           edges[edgeID] = std::make_shared<thisHandlerT>(std::move(tmp));
0285         }
0286       };
0287       MetaUtils::staticFor<sizeof...(AxisT)>(extractEdges);
0288
0289       std::vector<size_t> masks = BaseT::_binning.maskedBins();
0290       nedges->extend({ masks.size() });
0291       nedges->extend(std::move(masks));
0292
0293       const auto& itr = labels.cbegin();
0294       const auto& itrEnd = labels.cend();
0295       for (const auto& b : BaseT::bins(true, true)) {
0296         vector<string> keys = b.sources();
0297         nedges->extend({ keys.size() });
0298         for (const string& k : keys) {
0299           size_t dist = std::find(itr, itrEnd, k) - itr;
0300           nedges->extend({ dist });
0301         }
0302       }
0303     }
0304
0305     #endif
0306
0307     /// @}
0308
0309
0310     /// @name Transformations
0311     /// @{
0312
0313     /// @brief Rescale as if all fill weights had been different by factor @a scalefactor.
0314     void scale(const double scalefactor) noexcept {
0315       setAnnotation("ScaledBy", annotation<double>("ScaledBy", 1.0) * scalefactor);
0316       for (auto& bin : BaseT::_bins) {
0317         bin.scale(scalefactor);
0318       }
0319     }
0320
0321
0322     /// @brief Merge every group of @a n bins, from start to end inclusive
0323     ///
0324     /// If the number of bins is not a multiple of @a n, the last @a m < @a n
0325     /// bins on the RHS will also be merged, as the closest possible approach to
0326     /// factor @n rebinning everywhere.
0327     ///
0328     /// @note Only visible bins are being rebinned. Underflow (index = 0) and
0329     /// overflow (index = numBins() + 1) are not included.
0330     ///
0331     /// @note Only supported for continuous axes
0332     template <size_t axisN, typename = std::enable_if_t<BinningT::template is_CAxis<axisN>::value>>
0333     void rebinBy(unsigned int n, size_t begin=1, size_t end=UINT_MAX) {
0334       if (n < 1)  throw UserError("Rebinning requested in groups of 0!");
0335       if (!begin) throw UserError("Visible bins start with index 1!");
0336       if (end > BaseT::numBinsAt(axisN)+1)  end = BaseT::numBinsAt(axisN) + 1;
0337       for (size_t m = begin; m < end; ++m) {
0338         if (m > BaseT::numBinsAt(axisN)+1) break; // nothing to be done
0339         const size_t myend = (m+n-1 < BaseT::numBinsAt(axisN)+1) ? m+n-1 : BaseT::numBinsAt(axisN);
0340         if (myend > m) {
0341           BaseT::template mergeBins<axisN>({m, myend});
0342           end -= myend-m; //< reduce upper index by the number of removed bins
0343         }
0344       }
0345     }
0346
0347     /// @brief Overloaded alias for rebinBy
0348     ///
0349     /// @note Only supported for continuous axes
0350     template <size_t axisN, typename = std::enable_if_t<BinningT::template is_CAxis<axisN>::value>>
0351     void rebin(unsigned int n, size_t begin=1, size_t end=UINT_MAX) {
0352       rebinBy<axisN>(n, begin, end);
0353     }
0354
0355     /// @brief Rebin to the given list of bin edges
0356     ///
0357     /// @note Only supported for continuous axes
0358     template <size_t axisN, typename = std::enable_if_t<BinningT::template is_CAxis<axisN>::value>>
0359     void rebinTo(const std::vector<typename BinningT::template getAxisT<axisN>::EdgeT>& newedges) {
0360       if (newedges.size() < 2)
0361         throw UserError("Requested rebinning to an edge list which defines no bins");
0362       using thisAxisT = typename BinningT::template getAxisT<axisN>;
0363       using thisEdgeT = typename thisAxisT::EdgeT;
0364       // get list of shared edges
0365       thisAxisT& oldAxis = BaseT::_binning.template axis<axisN>();
0366       const thisAxisT newAxis(std::move(newedges));
0367       const std::vector<thisEdgeT> eshared = oldAxis.sharedEdges(newAxis);
0368       if (eshared.size() != newAxis.edges().size())
0369         throw BinningError("Requested rebinning to incompatible edges");
0370       // loop over new lower bin edges (= first bin index of merge range)
0371       for (size_t begin = 0; begin < eshared.size() - 1; ++begin) {
0372         // find index of upper edge along old axis
0373         // (subtracting 1 gives index of last bin to be merged)
0374         size_t end = oldAxis.index(eshared[begin+1]) - 1;
0375         // if the current edge is the last visible edge before the overflow
0376         // merge the remaining bins into the overflow
0377         if (begin == newAxis.numBins(true)-1)  end = oldAxis.numBins(true)-1;
0378         // merge this range
0379         if (end > begin)  BaseT::template mergeBins<axisN>({begin, end});
0380         if (eshared.size() == oldAxis.edges().size())  break; // we're done
0381       }
0382     }
0383
0384     /// @brief Overloaded alias for rebinTo
0385     ///
0386     /// @note Only supported for continuous axes
0387     template <size_t axisN, typename = std::enable_if_t<BinningT::template is_CAxis<axisN>::value>>
0388     void rebin(const std::vector<typename BinningT::template getAxisT<axisN>::EdgeT>& newedges) {
0389       rebinTo<axisN>(std::move(newedges));
0390     }
0391
0392     /// @brief Reset the EstimateStorage
0393     ///
0394     /// Keep the binning but set all bin contents and related quantities to zero
0395     void reset() noexcept { BaseT::clearBins(); }
0396
0397     /// @brief Copy assignment
0398     EstimateStorage& operator = (const EstimateStorage& est) noexcept {
0399       if (this != &est) {
0400         AnalysisObject::operator = (est);
0401         BaseT::operator = (est);
0402       }
0403       return *this;
0404     }
0405
0406     /// Move assignment
0407     EstimateStorage& operator = (EstimateStorage&& est) noexcept {
0408       if (this != &est) {
0409         AnalysisObject::operator = (est);
0410         BaseT::operator = (std::move(est));
0411       }
0412       return *this;
0413     }
0414
0415     /// @brief Add two EstimateStorages
0416     ///
0417     /// @note Adding EstimateStorages will unset any ScaledBy
0418     /// attribute from previous calls to scale or normalize.
0419     EstimateStorage& add(const EstimateStorage& est,
0420                          const std::string& pat_uncorr="^stat|^uncor" ) {
0421       if (*this != est)
0422         throw BinningError("Arithmetic operation requires compatible binning!");
0423       if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
0424       for (size_t i = 0; i< BaseT::numBins(true, true); ++i) {
0425         BaseT::bin(i).add(est.bin(i), pat_uncorr);
0426       }
0427       BaseT::maskBins(est.maskedBins(), true);
0428       return *this;
0429     }
0430     //
0431     EstimateStorage& operator += (const EstimateStorage& est) {
0432       return add(est);
0433     }
0434
0435     /// @brief Add two (rvalue) EstimateStorages
0436     ///
0437     /// @note Adding EstimateStorages will unset any ScaledBy
0438     /// attribute from previous calls to scale or normalize.
0439     EstimateStorage& add(EstimateStorage&& est,
0440                          const std::string& pat_uncorr="^stat|^uncor" ) {
0441       if (*this != est)
0442         throw BinningError("Arithmetic operation requires compatible binning!");
0443       if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
0444       for (size_t i = 0; i< BaseT::numBins(true, true); ++i) {
0445         BaseT::bin(i).add(std::move(est.bin(i)), pat_uncorr);
0446       }
0447       BaseT::maskBins(est.maskedBins(), true);
0448       return *this;
0449     }
0450     //
0451     EstimateStorage& operator += (EstimateStorage&& est) {
0452       return add(std::move(est));
0453     }
0454
0455
0456     /// @brief Subtract one EstimateStorages from another one
0457     ///
0458     /// @note Subtracting EstimateStorages will unset any ScaledBy
0459     /// attribute from previous calls to scale or normalize.
0460     EstimateStorage& subtract(const EstimateStorage& est,
0461                               const std::string& pat_uncorr="^stat|^uncor" ) {
0462       if (*this != est)
0463         throw BinningError("Arithmetic operation requires compatible binning!");
0464       if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
0465       for (size_t i = 0; i< BaseT::numBins(true, true); ++i) {
0466         BaseT::bin(i).subtract(est.bin(i), pat_uncorr);
0467       }
0468       BaseT::maskBins(est.maskedBins(), true);
0469       return *this;
0470     }
0471     //
0472     EstimateStorage& operator -= (const EstimateStorage& est) {
0473       return subtract(est);
0474     }
0475
0476     /// @brief Subtract one (rvalue) EstimateStorages from another one
0477     EstimateStorage& subtract(EstimateStorage&& est,
0478                               const std::string& pat_uncorr="^stat|^uncor" ) {
0479       if (*this != est)
0480         throw BinningError("Arithmetic operation requires compatible binning!");
0481       if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
0482       for (size_t i = 0; i< BaseT::numBins(true, true); ++i) {
0483         BaseT::bin(i) -= std::move(est.bin(i), pat_uncorr);
0484       }
0485       BaseT::maskBins(est.maskedBins(), true);
0486       return *this;
0487     }
0488     //
0489     EstimateStorage& operator -= (EstimateStorage&& est) {
0490       return subtract(std::move(est));
0491     }
0492
0493     /// @}
0494
0495
0496     /// @name Binning info.
0497     /// @{
0498
0499     /// @brief Total dimension of the object (number of axes + estimate)
0500     size_t dim() const noexcept { return sizeof...(AxisT) + 1; }
0501
0502     /// @brief Returns the axis configuration
0503     std::string _config() const noexcept { return mkAxisConfig<AxisT...>(); }
0504
0505     /// @}
0506
0507
0508     /// @name Whole EstimateStorage data
0509     /// @{
0510
0511     /// @brief Get list of central values
0512     std::vector<double> vals(const bool includeOverflows=false,
0513                              const bool includeMaskedBins=false) const {
0514       std::vector<double> rtn;
0515       rtn.reserve(BaseT::numBins(includeOverflows, includeMaskedBins));
0516       for (const auto& b : BaseT::bins(includeOverflows, includeMaskedBins)) {
0517         rtn.push_back(b.val());
0518       }
0519       return rtn;
0520     }
0521
0522     /// @brief Get list of error sources
0523     std::vector<std::string> sources() const {
0524       // Assemble union of error sources, as it's not guaranteed
0525       // that every bin has the same error breakdown
0526       std::vector<std::string> rtn;
0527
0528       for (const auto& b : BaseT::bins(true)) {
0529         std::vector<std::string> keys = b.sources();
0530         rtn.insert(std::end(rtn),
0531                    std::make_move_iterator(std::begin(keys)),
0532                    std::make_move_iterator(std::end(keys)));
0533       }
0534       std::sort(rtn.begin(), rtn.end());
0535       rtn.erase( std::unique(rtn.begin(), rtn.end()), rtn.end() );
0536
0537       return rtn;
0538     }
0539
0540     /// @brief Calculate the volume underneath the EstimateStorage
0541     ///
0542     /// @note overflow bins have infinite bin width
0543     double areaUnderCurve(const bool includeBinVol=true,
0544                           const bool includeOverflows=false,
0545                           const bool includeMaskedBins=false) const {
0546       double ret = 0.;
0547       for (const auto& b : BaseT::bins(includeOverflows, includeMaskedBins)) {
0548         const double val = fabs(b.val());
0549         const double vol = includeBinVol? b.dVol() : 1.0;
0550         if (std::isfinite(vol))  ret += val*vol; // aggregate bin volume
0551       }
0552       return ret;
0553     }
0554
0555     /// @brief Convenient alias for areaUnderCurve()
0556     double auc(const bool includeBinVol=true,
0557                const bool includeOverflows=false,
0558                const bool includeMaskedBins=false) const {
0559       return areaUnderCurve(includeBinVol, includeOverflows, includeMaskedBins);
0560     }
0561
0562
0563
0564     /// @brief Construct a covariance matrix from the error breakdown
0565     std::vector<std::vector<double>> covarianceMatrix(const bool ignoreOffDiagonalTerms=false,
0566                                                       const bool includeOverflows=false,
0567                                                       const bool includeMaskedBins=false,
0568                                                       const std::string& pat_uncorr="^stat|^uncor") const {
0569       const size_t nBins = BaseT::numBins(includeOverflows,includeMaskedBins);
0570       std::vector<std::vector<double> > covM(nBins);
0571
0572       // initialise cov matrix to be the right shape
0573       for (size_t i = 0; i < nBins; ++i) {
0574         covM[i] = std::vector<double>(nBins, 0.0);
0575       }
0576
0577       const std::vector<std::string> error_sources = sources();
0578
0579       // nominal-only case, i.e. total uncertainty, labelled as empty string
0580       if (error_sources.size() == 1 && error_sources[0] == "") {
0581         size_t i = 0;
0582         for (const auto& b : BaseT::bins(includeOverflows,includeMaskedBins)) {
0583           covM[i][i] = b.hasSource("")? sqr(0.5*(b.err("").first+b.err("").second)) : 1.0;
0584           ++i;
0585         }
0586         return covM;
0587       }
0588
0589       // more interesting case where we actually have some uncertainty breakdown!
0590       std::smatch match;
0591       const std::regex re(pat_uncorr, std::regex_constants::icase);
0592       for (const std::string& sname : error_sources) {
0593         if (sname == "")  continue;
0594         std::vector<double> systErrs(nBins, 0.0);
0595         size_t k = 0;
0596         for (const auto& b : BaseT::bins(includeOverflows,includeMaskedBins)) {
0597           if (b.hasSource(sname)) {
0598             const auto& errs = b.err(sname); // dn-up pair
0599             systErrs[k] = 0.5 *( fabs(errs.first)+fabs(errs.second));
0600           }
0601           ++k;
0602         }
0603         const bool skipOffDiag = (ignoreOffDiagonalTerms
0604                                   || std::regex_search(sname, match, re));
0605         for (size_t i = 0; i < nBins; ++i) {
0606           for (size_t j = 0; j < nBins; ++j) {
0607             if (skipOffDiag && i != j)  continue;
0608             covM[i][j] += systErrs[i]*systErrs[j];
0609           }
0610         }
0611       }
0612
0613       return covM;
0614     }
0615
0616     /// @}
0617
0618     /// @name Type reductions
0619     /// @{
0620
0621     /// @brief Produce a ScatterND from a EstimateStorage
0622     auto mkScatter(const std::string& path="", const std::string& pat_match = "",
0623                                                const bool includeOverflows=false,
0624                                                const bool includeMaskedBins=false) const {
0625
0626       constexpr size_t N = sizeof...(AxisT);
0627
0628       ScatterND<N+1> rtn;
0629       for (const std::string& a : annotations()) {
0630         if (a != "Type")  rtn.setAnnotation(a, annotation(a));
0631       }
0632       rtn.setAnnotation("Path", path);
0633
0634       const bool incOF = all_CAxes<AxisT...>::value && includeOverflows;
0635       for (const auto& b : BaseT::bins(incOF, includeMaskedBins)) {
0636
0637         // Create the vector of coordinates
0638         Utils::ndarray<double, N+1> vals;
0639         // first fill bin centres, use bin index if axis non-arithmetic
0640         auto indexIfDiscrete = [&vals, &b](auto I) {
0641           using isContinuous = typename BinningT::template is_CAxis<I>;
0642           using isArithmetic = typename BinningT::template is_Arithmetic<I>;
0643           vals[I] = coordPicker<I>(b, std::integral_constant<bool, isArithmetic::value>{},
0644                                       std::integral_constant<bool, isContinuous::value>{});
0645         };
0646         MetaUtils::staticFor<BinningT::Dimension::value>(indexIfDiscrete);
0647         // then fill bin content
0648         vals[N] = b.val();
0649
0650         // Create the vector of error pairs, use 0 if axis not continuous
0651         Utils::ndarray<std::pair<double,double>, N+1> errs;
0652         auto nullifyDiscrete = [&errs, &vals, &b](auto I) {
0653           using isContinuous = typename BinningT::template is_CAxis<I>;
0654           errs[I] = nullifyIfDisc<I>(b, vals[I], std::integral_constant<bool, isContinuous::value>{});
0655         };
0656         MetaUtils::staticFor<BinningT::Dimension::value>(nullifyDiscrete);
0657         const std::pair<double,double> tot = b.totalErr(pat_match);
0658         errs[N] = {fabs(tot.first), tot.second};
0659
0660         // Add the PointND
0661         rtn.addPoint( PointND<N+1>(vals, errs) );
0662       }
0663
0664       // Decorate output scatter with the discrete edges
0665       const BinningT& binning = BaseT::_binning;
0666       auto decorateEdges = [&rtn, &binning](auto I) {
0667         using isContinuous = typename BinningT::template is_CAxis<I>;
0668         if constexpr( !isContinuous::value ) {
0669           const auto& axis  = binning.template axis<I>();
0670           if (axis.numBins()) {
0671             std::stringstream ss;
0672             axis._renderYODA(ss);
0673             rtn.setAnnotation("EdgesA" + std::to_string(I+1), ss.str());
0674           }
0675         }
0676       };
0677       MetaUtils::staticFor<BinningT::Dimension::value>(decorateEdges);
0678
0679       return rtn;
0680     }
0681
0682     /// @brief Split into vector of BinnedEstimates along @a axisN
0683     ///
0684     /// The binning dimension of the returned objects are reduced by one unit.
0685     /// @note Requires at least two binning dimensions.
0686     template<size_t axisN, typename = std::enable_if_t< (axisN < sizeof...(AxisT) && sizeof...(AxisT)>=2) >>
0687     auto mkEstimates(const std::string& path="", const bool includeOverflows=false) const {
0688
0689       // Need to provide a prescription for how to add the two bin contents
0690       auto how2add = [](auto& pivot, const BinType& toCopy) { pivot = toCopy; };
0691       auto rtn = BaseT::template mkBinnedSlices<axisN, BinnedEstimate>(how2add, includeOverflows);
0692       for (const std::string& a : annotations()) {
0693         if (a == "Type")  continue;
0694         for (size_t i = 0; i < rtn.size(); ++i) {
0695           rtn[i].setAnnotation(a, annotation(a));
0696         }
0697       }
0698       for (size_t i = 0; i < rtn.size(); ++i) {
0699         rtn[i].setAnnotation("Path", path);
0700       }
0701
0702       return rtn;
0703     }
0704
0705
0706     /// @brief Method returns clone of the estimate with streamlined error source
0707     AnalysisObject* mkInert(const std::string& path = "",
0708                             const std::string& source = "") const noexcept {
0709       EstimateStorage* rtn = newclone();
0710       rtn->setPath(path);
0711       for (auto& b : rtn->bins(true, true)) {
0712         if (b.numErrs() == 1) {
0713           try {
0714             b.renameSource("", source);
0715           }
0716           catch (YODA::UserError& e) { }
0717         }
0718       }
0719       return rtn;
0720     }
0721
0722     /// @}
0723
0724     /// @name MPI (de-)serialisation
0725     /// @{
0726
0727     size_t lengthContent(bool fixed_length = false) const noexcept {
0728       size_t rtn = 0;
0729       for (const auto& bin : BaseT::bins(true, true)) {
0730         rtn += bin._lengthContent(fixed_length);
0731       }
0732       return rtn;
0733     }
0734
0735     std::vector<double> serializeContent(bool fixed_length = false) const noexcept {
0736       std::vector<double> rtn;
0737       const size_t nBins = BaseT::numBins(true, true);
0738       rtn.reserve(nBins * 4);
0739       for (size_t i = 0; i < nBins; ++i) {
0740         const auto& b = BaseT::bin(i);
0741         std::vector<double> bdata = b._serializeContent(fixed_length);
0742         rtn.insert(std::end(rtn),
0743                    std::make_move_iterator(std::begin(bdata)),
0744                    std::make_move_iterator(std::end(bdata)));
0745       }
0746       return rtn;
0747     }
0748
0749
0750     void deserializeContent(const std::vector<double>& data) {
0751
0752       const size_t nBins = BaseT::numBins(true, true);
0753       const size_t minLen = 2*nBins;
0754       if (data.size() < minLen)
0755         throw UserError("Length of serialized data should be at least " + std::to_string(minLen)+"!");
0756
0757       size_t i = 0;
0758       auto itr = data.cbegin();
0759       const auto itrEnd = data.cend();
0760       const bool fixedLen = data.size() == 2*minLen;
0761       while (itr != itrEnd) {
0762         // for estimates, the first element represents the central,
0763         // the subsequent value represents the number of error pairs
0764         const size_t nErrs = fixedLen? 1 : (*(itr + 1) + 0.5); // add 0.5 to avoid rounding issues
0765         auto last = itr + 2*(nErrs+1); // last element + 1
0766         BaseT::bin(i)._deserializeContent(std::vector<double>{itr, last}, fixedLen);
0767         // update for next iteration
0768         itr = last;
0769         ++i;
0770       }
0771     }
0772
0773     /// @}
0774
0775   private:
0776
0777     /// @brief Helper function to create a BinningT from
0778     /// a given set @a nBins within a range @a limitsLowUp
0779     template<size_t... Is>
0780     BinningT _mkBinning(const std::vector<size_t>& nBins,
0781                         const std::vector<std::pair<double, double>>& limitsLowUp,
0782                         std::index_sequence<Is...>) const {
0783       return BinningT({((void)Is, Axis<AxisT>(nBins[Is], limitsLowUp[Is].first, limitsLowUp[Is].second))...});
0784     }
0785
0786     /// @brief Helper function to create a BinningT from a scatter @a s
0787     template<size_t... Is>
0788     BinningT _mkBinning(const ScatterND<sizeof...(AxisT)+1>& s, std::index_sequence<Is...>) const {
0789       return BinningT({((void)Is, Axis<AxisT>(s.edges(Is)))...});
0790     }
0791
0792   };
0793
0794
0795
0796   /// @brief User-facing BinnedEstimate class in arbitrary dimensions
0797   template <typename... AxisT>
0798   class BinnedEstimate : public EstimateStorage<AxisT...> {
0799   public:
0800     using EstimateT = BinnedEstimate<AxisT...>;
0801     using BaseT = EstimateStorage<AxisT...>;
0802     using BinType = typename BaseT::BinT;
0803     using Ptr = std::shared_ptr<EstimateT>;
0804
0805     /// @brief Inherit constructors.
0806     using BaseT::BaseT;
0807
0808     BinnedEstimate(const EstimateT&) = default;
0809     BinnedEstimate(EstimateT&&) = default;
0810     BinnedEstimate& operator =(const EstimateT&) = default;
0811     BinnedEstimate& operator =(EstimateT&&) = default;
0812     using AnalysisObject::operator =;
0813
0814     /// @brief Copy constructor (needed for clone functions).
0815     ///
0816     /// @note Compiler won't generate this constructor automatically.
0817     BinnedEstimate(const BaseT& other) : BaseT(other) {}
0818     //
0819     BinnedEstimate(const EstimateT& other, const std::string& path) : BaseT(other, path) {}
0820
0821     /// @brief Move constructor.
0822     BinnedEstimate(BaseT&& other) : BaseT(std::move(other)) {}
0823     //
0824     BinnedEstimate(EstimateT&& other, const std::string& path) : BaseT(std::move(other), path) {}
0825
0826   };
0827
0828
0829   /// @brief Specialisation of the BinnedEstimate for a 1D histogram
0830   template <typename AxisT>
0831   class BinnedEstimate<AxisT>
0832       : public EstimateStorage<AxisT>,
0833         public XAxisMixin<BinnedEstimate<AxisT>, AxisT> {
0834   public:
0835     using EstimateT = BinnedEstimate<AxisT>;
0836     using BaseT = EstimateStorage<AxisT>;
0837     using BinType = typename BaseT::BinT;
0838     using Ptr = std::shared_ptr<EstimateT>;
0839
0840     /// @brief Inherit constructors.
0841     using BaseT::BaseT;
0842
0843     BinnedEstimate(const EstimateT&) = default;
0844     BinnedEstimate(EstimateT&&) = default;
0845     BinnedEstimate& operator =(const EstimateT&) = default;
0846     BinnedEstimate& operator =(EstimateT&&) = default;
0847     using AnalysisObject::operator =;
0848
0849     /// @brief Copy constructor (needed for clone functions).
0850     ///
0851     /// @note Compiler won't generate this constructor automatically.
0852     BinnedEstimate(const BaseT& other) : BaseT(other) {}
0853     //
0854     BinnedEstimate(const EstimateT& other, const std::string& path) : BaseT(other, path) {}
0855
0856     /// @brief Move constructor.
0857     BinnedEstimate(BaseT&& other) : BaseT(std::move(other)) {}
0858     //
0859     BinnedEstimate(EstimateT&& other, const std::string& path) : BaseT(std::move(other), path) {}
0860
0861
0862     BinnedEstimate(std::vector<AxisT>&& edges, const std::string& path="", const std::string& title="")
0863               : BaseT(std::move(edges), path, title) {}
0864
0865     BinnedEstimate(const std::vector<AxisT>& edges, const std::string& path="", const std::string& title="")
0866               : BaseT(edges, path, title) {}
0867
0868     /// @brief Constructor with auto-setup of evenly spaced axes.
0869     ///
0870     /// The constructor argument uses double rather than EdgeT to
0871     /// allow for auto-conversion of int to double.
0872     ///
0873     /// @note This constructor is only supported when all axes are continuous.
0874     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT>>
0875     BinnedEstimate(size_t nbins, double lower, double upper, const std::string& path = "", const std::string& title = "")
0876               : BaseT({nbins}, {{lower, upper}}, path, title) {}
0877
0878     /// @brief Make a copy on the stack
0879     EstimateT clone() const noexcept {
0880       return EstimateT(*this);
0881     }
0882
0883     /// @brief Make a copy on the heap
0884     EstimateT* newclone() const noexcept {
0885       return new EstimateT(*this);
0886     }
0887
0888     /// @brief Find bin index for given coordinates
0889     size_t indexAt(const AxisT xCoord) const noexcept {
0890       return BaseT::binAt( {xCoord} ).index();
0891     }
0892
0893     /// @brief Mask/Unmask bin at given set of coordinates
0894     void maskBinAt(const AxisT xCoord, const bool status = true) noexcept {
0895       return BaseT::maskBin({xCoord}, status);
0896     }
0897
0898   };
0899
0900
0901
0902   /// @brief Specialisation of the BinnedEstimate for a 2D BinnedEstimate
0903   template <typename AxisT1, typename AxisT2>
0904   class BinnedEstimate<AxisT1, AxisT2>
0905       : public EstimateStorage<AxisT1, AxisT2>,
0906         public XAxisMixin<BinnedEstimate<AxisT1, AxisT2>, AxisT1>,
0907         public YAxisMixin<BinnedEstimate<AxisT1, AxisT2>, AxisT2> {
0908   public:
0909     using EstimateT = BinnedEstimate<AxisT1, AxisT2>;
0910     using BaseT = EstimateStorage<AxisT1, AxisT2>;
0911     using BinType = typename BaseT::BinT;
0912     using Ptr = std::shared_ptr<EstimateT>;
0913
0914     /// @brief Inherit constructors.
0915     using BaseT::BaseT;
0916
0917     BinnedEstimate(const EstimateT&) = default;
0918     BinnedEstimate(EstimateT&&) = default;
0919     BinnedEstimate& operator =(const EstimateT&) = default;
0920     BinnedEstimate& operator =(EstimateT&&) = default;
0921     using AnalysisObject::operator =;
0922
0923     /// @brief Copy constructor (needed for clone functions).
0924     ///
0925     /// @note Compiler won't generate this constructor automatically.
0926     BinnedEstimate(const BaseT& other) : BaseT(std::move(other)) {}
0927     //
0928     BinnedEstimate(const EstimateT& other, const std::string& path) : BaseT(other, path) {}
0929
0930     /// @brief Move constructor.
0931     BinnedEstimate(BaseT&& other) : BaseT(std::move(other)) {}
0932     //
0933     BinnedEstimate(EstimateT&& other, const std::string& path) : BaseT(std::move(other), path) {}
0934
0935     BinnedEstimate(std::vector<AxisT1>&& xEdges, std::vector<AxisT2>&& yEdges,
0936                    const std::string& path="", const std::string& title="")
0937               : BaseT(std::move(xEdges), std::move(yEdges), path, title) {}
0938
0939     BinnedEstimate(const std::vector<AxisT1>& xEdges, const std::vector<AxisT2>& yEdges,
0940                    const std::string& path="", const std::string& title="")
0941               : BaseT(xEdges, yEdges, path, title) {}
0942
0943     /// @brief Constructor with auto-setup of evenly spaced axes.
0944     ///
0945     /// The constructor argument uses double rather than EdgeT to
0946     /// allow for auto-conversion of int to double.
0947     ///
0948     /// @note This constructor is only supported when all axes are continuous.
0949     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT1, AxisT2>>
0950     BinnedEstimate(size_t nbinsX, double lowerX, double upperX,
0951                    size_t nbinsY, double lowerY, double upperY,
0952                    const std::string& path = "", const std::string& title = "")
0953               : BaseT({nbinsX, nbinsY}, {{lowerX, upperX}, {lowerY, upperY}}, path, title) {}
0954
0955     /// @brief Make a copy on the stack
0956     EstimateT clone() const noexcept {
0957       return EstimateT(*this);
0958     }
0959
0960     /// @brief Make a copy on the heap
0961     EstimateT* newclone() const noexcept {
0962       return new EstimateT(*this);
0963     }
0964
0965     /// @brief Bin access using global index
0966     BinType& bin(const size_t index) noexcept {
0967       return BaseT::bin(index);
0968     }
0969
0970     /// @brief Bin access using global index (const version)
0971     const BinType& bin(const size_t index) const noexcept {
0972       return BaseT::bin(index);
0973     }
0974
0975     /// @brief Bin access using local indices
0976     BinType& bin(const size_t localX, const size_t localY) noexcept {
0977       return BaseT::bin( {localX, localY} );
0978     }
0979
0980     /// @brief Bin access using local indices (const version)
0981     const BinType& bin(const size_t localX, const size_t localY) const noexcept {
0982       return BaseT::bin( {localX, localY} );
0983     }
0984
0985     /// @brief Bin access using coordinates
0986     BinType& binAt(const AxisT1 xCoord, const AxisT2 yCoord) noexcept {
0987       return BaseT::binAt( {xCoord, yCoord} );
0988     }
0989
0990     /// @brief Bin access using coordinates (const version)
0991     const BinType& binAt(const AxisT1 xCoord, const AxisT2 yCoord) const noexcept {
0992       return BaseT::binAt( {xCoord, yCoord} );
0993     }
0994
0995     /// @brief Find bin index for given coordinates
0996     size_t indexAt(const AxisT1 xCoord, const AxisT2 yCoord) const noexcept {
0997       return BaseT::binAt( {xCoord, yCoord} ).index();
0998     }
0999
1000     /// @brief Mask/Unmask bin at given set of coordinates
1001     void maskBinAt(const AxisT1 xCoord, const AxisT2 yCoord, const bool status = true) noexcept {
1002       return BaseT::maskBin({xCoord, yCoord}, status);
1003     }
1004
1005   };
1006
1007
1008   /// @brief Specialisation of the BinnedEstimate for a 3D BinnedEstimate
1009   template <typename AxisT1, typename AxisT2, typename AxisT3>
1010   class BinnedEstimate<AxisT1, AxisT2, AxisT3>
1011       : public EstimateStorage<AxisT1, AxisT2, AxisT3>,
1012         public XAxisMixin<BinnedEstimate<AxisT1, AxisT2, AxisT3>, AxisT1>,
1013         public YAxisMixin<BinnedEstimate<AxisT1, AxisT2, AxisT3>, AxisT2>,
1014         public ZAxisMixin<BinnedEstimate<AxisT1, AxisT2, AxisT3>, AxisT3> {
1015   public:
1016     using EstimateT = BinnedEstimate<AxisT1, AxisT2, AxisT3>;
1017     using BaseT = EstimateStorage<AxisT1, AxisT2, AxisT3>;
1018     using BinType = typename BaseT::BinT;
1019     using Ptr = std::shared_ptr<EstimateT>;
1020
1021     /// @brief Inherit constructors.
1022     using BaseT::BaseT;
1023
1024     BinnedEstimate(const EstimateT&) = default;
1025     BinnedEstimate(EstimateT&&) = default;
1026     BinnedEstimate& operator =(const EstimateT&) = default;
1027     BinnedEstimate& operator =(EstimateT&&) = default;
1028     using AnalysisObject::operator =;
1029
1030     /// @brief Copy constructor (needed for clone functions).
1031     ///
1032     /// @note Compiler won't generate this constructor automatically.
1033     BinnedEstimate(const BaseT& other) : BaseT(other) {}
1034     //
1035     BinnedEstimate(const EstimateT& other, const std::string& path) : BaseT(other, path) {}
1036
1037     /// @brief Move constructor.
1038     BinnedEstimate(BaseT&& other) : BaseT(std::move(other)) {}
1039     //
1040     BinnedEstimate(EstimateT&& other, const std::string& path) : BaseT(std::move(other), path) {}
1041
1042     /// @brief Constructor with auto-setup of evenly spaced axes.
1043     ///
1044     /// The constructor argument uses double rather than EdgeT to
1045     /// allow for auto-conversion of int to double.
1046     ///
1047     /// @note This constructor is only supported when all axes are continuous.
1048     template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT1, AxisT2, AxisT3>>
1049     BinnedEstimate(size_t nbinsX, double lowerX, double upperX,
1050                    size_t nbinsY, double lowerY, double upperY,
1051                    size_t nbinsZ, double lowerZ, double upperZ,
1052                    const std::string& path = "", const std::string& title = "")
1053               : BaseT({nbinsX, nbinsY, nbinsZ}, {{lowerX, upperX}, {lowerY, upperY},
1054                       {lowerZ, upperZ}}, path, title) {}
1055
1056     /// @brief Make a copy on the stack
1057     EstimateT clone() const noexcept {
1058       return EstimateT(*this);
1059     }
1060
1061     /// @brief Make a copy on the heap
1062     EstimateT* newclone() const noexcept {
1063       return new EstimateT(*this);
1064     }
1065
1066     /// @brief Bin access using global index
1067     BinType& bin(const size_t index) noexcept {
1068       return BaseT::bin(index);
1069     }
1070
1071     /// @brief Bin access using global index (const version)
1072     const BinType& bin(const size_t index) const noexcept {
1073       return BaseT::bin(index);
1074     }
1075
1076     /// @brief Bin access using local indices
1077     BinType& bin(const size_t localX, const size_t localY, const size_t localZ) noexcept {
1078       return BaseT::bin( {localX, localY, localZ} );
1079     }
1080
1081     /// @brief Bin access using local indices
1082     const BinType& bin(const size_t localX, const size_t localY, const size_t localZ) const noexcept {
1083       return BaseT::bin( {localX, localY, localZ} );
1084     }
1085
1086     /// @brief Bin access using coordinates
1087     BinType& binAt(const AxisT1 xCoord, const AxisT2 yCoord, const AxisT3 zCoord) noexcept {
1088       return BaseT::binAt( {xCoord, yCoord, zCoord} );
1089     }
1090
1091     /// @brief Bin access using coordinates (const version)
1092     const BinType& binAt(const AxisT1 xCoord, const AxisT2 yCoord, const AxisT3 zCoord) const noexcept {
1093       return BaseT::binAt( {xCoord, yCoord, zCoord} );
1094     }
1095
1096     /// @brief Find bin index for given coordinates
1097     size_t indexAt(const AxisT1 xCoord, const AxisT2 yCoord, const AxisT3 zCoord) const noexcept {
1098       return BaseT::binAt( {xCoord, yCoord, zCoord} ).index();
1099     }
1100
1101     /// @brief Mask/Unmask bin at given set of coordinates
1102     void maskBinAt(const AxisT1 xCoord, const AxisT2 yCoord, const AxisT3 zCoord, const bool status = true) noexcept {
1103       return BaseT::maskBin({xCoord, yCoord, zCoord}, status);
1104     }
1105
1106   };
1107
1108   /// @name Combining BinnedEstimates: global operators
1109   /// @{
1110
1111   /// @brief Add two BinnedEstimates
1112   template <typename... AxisT>
1113   inline BinnedEstimate<AxisT...>
1114   operator + (BinnedEstimate<AxisT...> first, const BinnedEstimate<AxisT...>& second) {
1115     first += second;
1116     return first;
1117   }
1118
1119
1120   /// @brief Add two BinnedEstimates
1121   template <typename... AxisT>
1122   inline BinnedEstimate<AxisT...>
1123   operator + (BinnedEstimate<AxisT...> first, BinnedEstimate<AxisT...>&& second) {
1124     first += std::move(second);
1125     return first;
1126   }
1127
1128
1129   /// @brief Subtract two BinnedEstimates
1130   template <typename... AxisT>
1131   inline BinnedEstimate< AxisT...>
1132   operator - (BinnedEstimate<AxisT...> first, const BinnedEstimate<AxisT...>& second) {
1133     first -= second;
1134     return first;
1135   }
1136
1137   /// @brief Subtract two BinnedEstimates
1138   template <typename... AxisT>
1139   inline BinnedEstimate< AxisT...>
1140   operator - (BinnedEstimate<AxisT...> first, BinnedEstimate<AxisT...>&& second) {
1141     first -= std::move(second);
1142     return first;
1143   }
1144
1145   /// @brief Divide two BinnedEstimates
1146   template <typename... AxisT>
1147   inline BinnedEstimate<AxisT...>
1148   divide(const BinnedEstimate<AxisT...>& numer, const BinnedEstimate<AxisT...>& denom,
1149          const std::string& pat_uncorr="^stat|^uncor" ) {
1150     if (numer != denom) {
1151       throw BinningError("Arithmetic operation requires compatible binning!");
1152     }
1153
1154     BinnedEstimate<AxisT...> rtn(numer.binning());
1155     if (numer.path() == denom.path())  rtn.setPath(numer.path());
1156     if (rtn.hasAnnotation("ScaledBy")) rtn.rmAnnotation("ScaledBy");
1157
1158     for (const auto& b_num : numer.bins(true, true)) {
1159       const size_t idx = b_num.index();
1160       rtn.bin(idx) = divide(b_num, denom.bin(idx), pat_uncorr);
1161     }
1162     rtn.maskBins(denom.maskedBins(), true);
1163
1164     return rtn;
1165   }
1166
1167   /// @brief Divide two BinnedEstimates
1168   template <typename... AxisT>
1169   inline BinnedEstimate<AxisT...>
1170   operator / (const BinnedEstimate<AxisT...>& numer, const BinnedEstimate<AxisT...>& denom) {
1171     return divide(numer, denom);
1172   }
1173
1174   /// @brief Divide two BinnedEstimates
1175   template <typename... AxisT>
1176   inline BinnedEstimate<AxisT...>
1177   operator / (BinnedEstimate<AxisT...>&& numer, const BinnedEstimate<AxisT...>& denom) {
1178     return divide(std::move(numer), denom);
1179   }
1180
1181   /// @brief Divide two BinnedEstimates
1182   template <typename... AxisT>
1183   inline BinnedEstimate<AxisT...>
1184   operator / (const BinnedEstimate<AxisT...>& numer, BinnedEstimate<AxisT...>&& denom) {
1185     return divide(numer, std::move(denom));
1186   }
1187
1188   /// @brief Divide two BinnedEstimates
1189   template <typename... AxisT>
1190   inline BinnedEstimate<AxisT...>
1191   operator / (BinnedEstimate<AxisT...>&& numer, BinnedEstimate<AxisT...>&& denom) {
1192     return divide(std::move(numer), std::move(denom));
1193   }
1194
1195
1196   /// @brief Calculate a binned efficiency ratio of two BinnedEstimate objects
1197   ///
1198   /// @note An efficiency is not the same thing as a standard division of two
1199   /// BinnedEstimate objects: the errors are treated as correlated via binomial statistics.
1200   template <typename... AxisT>
1201   inline BinnedEstimate<AxisT...>
1202   efficiency(const BinnedEstimate<AxisT...>& accepted, const BinnedEstimate<AxisT...>& total,
1203              const std::string& pat_uncorr="^stat|^uncor" ) {
1204
1205     if (accepted != total) {
1206       throw BinningError("Arithmetic operation requires compatible binning!");
1207     }
1208
1209     BinnedEstimate<AxisT...> rtn(accepted.binning());
1210
1211     for (const auto& b_acc : accepted.bins(true, true)) {
1212       Estimate est;
1213       const size_t idx = b_acc.index();
1214       try {
1215         est = efficiency(b_acc, total.bin(idx), pat_uncorr);
1216       } catch (const UserError& e) {
1217         //
1218       }
1219       rtn.bin(idx).set(est);
1220     }
1221     return rtn;
1222   }
1223
1224
1225   /// @brief Calculate the asymmetry (a-b)/(a+b) of two BinnedDbn objects
1226   template <typename... AxisT>
1227   inline BinnedEstimate<AxisT...>
1228   asymm(const BinnedEstimate<AxisT...>& a,
1229         const BinnedEstimate<AxisT...>& b,
1230         const std::string& pat_uncorr="^stat|^uncor" ) {
1231     return divde(subtract(a-b, pat_uncorr), add(a+b, pat_uncorr), pat_uncorr);
1232   }
1233
1234
1235   /// @}
1236
1237
1238   /// @name Generalised transformations
1239   /// @{
1240
1241   template <typename... AxisT>
1242   inline void transform(BinnedEstimate<AxisT...>& est, const Trf<1>& fn) {
1243     for (auto& b : est.bins(true, true)) {
1244       b.transform(fn);
1245     }
1246   }
1247
1248   template <typename... AxisT, typename FN>
1249   inline void transform(BinnedEstimate<AxisT...>& est, const FN& fn) {
1250     transform(est, Trf<1>(fn));
1251   }
1252
1253   /// @}
1254
1255
1256   /// @name Convenience aliases
1257   /// @{
1258
1259   /// Define dimension-specific short-hands (Cython sugar)
1260   template<typename A1>
1261   using BinnedEstimate1D = BinnedEstimate<A1>;
1262
1263   template <typename A1, typename A2>
1264   using BinnedEstimate2D = BinnedEstimate<A1, A2>;
1265
1266   template <typename A1, typename A2, typename A3>
1267   using BinnedEstimate3D = BinnedEstimate<A1, A2, A3>;
1268
1269   /// Anonymous namespace to limit visibility
1270   namespace {
1271     template <class T>
1272     struct EstimateMaker;
1273
1274     template<size_t... Is>
1275     struct EstimateMaker<std::index_sequence<Is...>> {
1276       using type = BinnedEstimate< std::decay_t<decltype((void)Is, std::declval<double>())>... >;
1277     };
1278   }
1279
1280   template<size_t N>
1281   using EstimateND = typename EstimateMaker<std::make_index_sequence<N>>::type;
1282
1283
1284   // User-friendly names (continuous axes only)
1285   using Estimate1D = BinnedEstimate<double>;
1286   using Estimate2D = BinnedEstimate<double,double>;
1287   using Estimate3D = BinnedEstimate<double,double,double>;
1288
1289   /// @}
1290
1291 }
1292
1293 #endif