EIC code displayed by LXR master/​include/​YODA/​Estimate.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

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

 // -*- C++ -*-
 //
 // This file is part of YODA -- Yet more Objects for Data Analysis
 // Copyright (C) 2008-2024 The YODA collaboration (see AUTHORS for details)
 //
 #ifndef YODA_Estimate_h
 #define YODA_Estimate_h
 
 #include "YODA/Exceptions.h"
 #include "YODA/Dbn.h"
 #include "YODA/Transformation.h"
 #include "YODA/Utils/MathUtils.h"
 #include "YODA/Utils/StringUtils.h"
 #include <cmath>
 #include <map>
 #include <string>
 #include <regex>
 #include <utility>
 
 namespace YODA {
 
   /// @brief A point estimate (base class for the Estimate)
   ///
   /// This class is used internally by YODA to centralise the storage of
   /// static statistical results or estimators. Each estimate is made up of
   /// a value and an error map. The error map uses the source of uncertainty
   /// as the keys and the (possibly asymmetric) uncertainty pair as the values.
   /// The empty string is interpreted as the total uncertainty.
   class Estimate {
   public:
 
     /// @name Constructors
     //@{
 
     /// Default constructor of a new distribution.
     Estimate() { reset(); }
 
 
     /// @brief Constructor to set an Estimate with a pre-filled state.
     ///
     /// Principally designed for internal persistency use.
     Estimate(double v, std::map<std::string,std::pair<double,double>>& errors)
       : _value(v), _error(errors) { }
 
 
     /// @brief Alternative constructor to set an Estimate with value and uncertainty.
     Estimate(const double v, const std::pair<double,double>& e, const std::string& source = "")
       : _value(v) { setErr(e, source); }
 
     /// Copy constructor
     ///
     /// Sets all the parameters using the ones provided from an existing Estimate.
     Estimate(const Estimate& toCopy) {
       _value = toCopy._value;
       _error = toCopy._error;
     }
 
     /// Move constructor
     ///
     /// Sets all the parameters using the ones provided from an existing Estimate.
     Estimate(Estimate&& toMove) {
       _value = std::move(toMove._value);
       _error = std::move(toMove._error);
     }
 
     /// Copy assignment
     ///
     /// Sets all the parameters using the ones provided from an existing Estimate.
     Estimate& operator=(const Estimate& toCopy) noexcept {
       if (this != &toCopy) {
         _value = toCopy._value;
         _error = toCopy._error;
       }
       return *this;
     }
 
     /// Move assignment
     ///
     /// Sets all the parameters using the ones provided from an existing Estimate.
     Estimate& operator=(Estimate&& toMove) noexcept {
       if (this != &toMove) {
         _value = std::move(toMove._value);
         _error = std::move(toMove._error);
       }
       return *this;
     }
 
     //@}
 
 
     /// @name Operators
     //@{
 
     /// Add two Estimates
     Estimate& add(const Estimate& toAdd, const std::string& pat_uncorr="^stat|^uncor" ) {
       setVal(val() + toAdd.val());
       std::smatch match;
       const std::regex re(pat_uncorr, std::regex_constants::icase);
       for (const std::string& src : toAdd.sources()) {
         if (!hasSource(src))  setErr(toAdd.errDownUp(src), src);
         else if (std::regex_search(src, match, re)) {
           // treat as uncorrelated between AOs:
           // add in quadrature
           const auto [l_dn,l_up] = errDownUp(src);
           const auto [r_dn,r_up] = toAdd.errDownUp(src);
           const double new_dn = std::sqrt(l_dn*l_dn + r_dn*r_dn);
           const double new_up = std::sqrt(l_up*l_up + r_up*r_up);
           setErr({-new_dn,new_up}, src);
         }
         else {
           // treat as correlated between AOs: add linearly
           const auto [l_dn,l_up] = errDownUp(src);
           const auto [r_dn,r_up] = toAdd.errDownUp(src);
           setErr({l_dn+r_dn, l_up+r_up}, src);
         }
       }
       return *this;
     }
     //
     Estimate& operator+=(const Estimate& toAdd) {
       return add(toAdd);
     }
 
     /// Subtract one Estimate from another
     Estimate& subtract(const Estimate& toSubtract, const std::string& pat_uncorr="^stat|^uncor" ) {
       setVal(val() - toSubtract.val());
       std::smatch match;
       const std::regex re(pat_uncorr, std::regex_constants::icase);
       for (const std::string& src : toSubtract.sources()) {
         if (!hasSource(src))  setErr(toSubtract.errDownUp(src), src);
         else if (std::regex_search(src, match, re)) {
           // treat as uncorrelated between AOs: add in quadrature
           const auto [l_dn,l_up] = errDownUp(src);
           const auto [r_dn,r_up] = toSubtract.errDownUp(src);
           const double new_dn = std::sqrt(l_dn*l_dn + r_dn*r_dn);
           const double new_up = std::sqrt(l_up*l_up + r_up*r_up);
           setErr({-new_dn,new_up}, src);;
         }
         else {
           // treat as correlated between AOs: add linearly
           const auto [l_dn,l_up] = errDownUp(src);
           const auto [r_dn,r_up] = toSubtract.errDownUp(src);
           setErr({l_dn+r_dn, l_up+r_up}, src);
         }
       }
       return *this;
     }
     //
     Estimate& operator-=(const Estimate& toSubtract) {
       return subtract(toSubtract);
     }
 
     //@}
 
     /// @name Modifiers
     //@{
 
     /// @brief Set the central value of this estimator
     void setValue(const double value) noexcept { _value = value; }
 
     /// @brief Alias for setValue
     void setVal(const double val) noexcept { setValue(val); }
 
     /// @brief Set a signed uncertainty component.
     void setErr(const std::pair<double, double> &err, const std::string& source = "") {
       const std::string& s = Utils::toUpper(source);
       if (s == "TOTAL")
         throw UserError("Use empty string for the total uncertainty!");
       _error[source] = err;
     }
 
     /// @brief Set a symmetric uncertainty component
     ///
     /// @note The down-component will be set to -|err|
     /// and the up-component will be set to +|err|.
     void setErr(const double err, const std::string& source = "") {
       setErr({-fabs(err),fabs(err)}, source);
     }
 
     /// @brief Set both central value and uncertainty component
     void set(const double val, const std::pair<double, double> &err, const std::string source = "") {
       setVal(val);
       setErr(err, source);
     }
 
 
     /// @brief Set both central value and uncertainty component
     void set(const double val, const double err, const std::string& source = "") {
       setVal(val);
       setErr(err, source);
     }
 
     /// Reset the internal values
     void reset() noexcept {
       _value = std::numeric_limits<double>::quiet_NaN();
       _error.clear();
     }
 
 
     /// Rescale as if value and uncertainty had been different by factor @a scalefactor.
     void scale(const double scalefactor) noexcept {
       _value *= scalefactor;
       for (auto& item : _error) {
         item.second = { item.second.first * scalefactor, item.second.second * scalefactor };
       }
     }
 
     /// Generalised transformations with functors
     void scale(const Trf<1>& trf) {
       trf.transform(_value, _error);
     }
 
     /// Generalised transformations with functors
     void transform(const Trf<1>& trf) {
       scale(trf);
     }
 
     /// Replace a source label in the error breakdown
     void renameSource(const std::string& old_label, const std::string& new_label) {
       if (!hasSource(old_label)) {
         throw UserError("Error map has no such key: "+old_label);
       }
       auto entry = _error.extract(old_label);
       entry.key() = new_label;
       _error.insert(std::move(entry));
     }
 
     // Remove error source @a source
     void rmSource(const std::string& source) {
       if (hasSource(source)) {
         _error.erase(source);
       }
     }
 
     // Clear the error breakdown
     void rmErrs() {
       _error.clear();
     }
 
     //@}
 
   public:
 
 
     /// @name Estimate getters
     //@{
 
     /// The central value
     double val() const noexcept { return _value; }
 
     /// @brief The signed absolute error on the central value
     ///
     /// @note first/second element corresponds to
     /// systematic down/up variation
     std::pair<double,double> errDownUp(const std::string& source = "") const {
       const size_t count = _error.count(source);
       if (!count)
         throw RangeError("Error map has no such key: "+source);
       return _error.at(source);
     }
 
     /// @brief Convenience alias for errorDownUp(source)
     std::pair<double,double> err(const std::string& source = "") const {
       return errDownUp(source);
     }
 
     /// @brief The signed negative and positive uncertainty component
     std::pair<double,double> errNegPos(std::string source = "") const {
       return _downUp2NegPos( errDownUp(source) );
     }
 
     /// @brief The signed negative uncertainty component
     double errNeg(const std::string& source = "") const {
       return errNegPos(source).first;
     }
 
     /// @brief The signed positive uncertainty component
     double errPos(const std::string& source = "") const {
       return errNegPos(source).second;
     }
 
     /// @brief The signed error due to the systematic downward variation
     double errDown(const std::string& source = "") const {
       return errDownUp(source).first;
     }
 
     /// @brief The signed error due to the systematic upward variation
     double errUp(const std::string& source = "") const {
       return errDownUp(source).second;
     }
 
     /// @brief The unsigned error from an average of the downwards/upwards shifts
     double errAvg(const std::string& source = "") const {
       return _average(errDownUp(source));
     }
 
     /// @brief The unsigned symmetrised error from the largest component
     /// of the downwards/upwards shifts
     double errEnv(const std::string& source = "") const {
       return _envelope(errDownUp(source));
     }
 
     /// @brief The relative signed error pair with respect to the central value
     std::pair<double,double> relErrDownUp(const std::string& source = "") const {
       auto [dn,up] = errDownUp(source);
       dn = _value != 0 ? dn/fabs(_value) : std::numeric_limits<double>::quiet_NaN();
       up = _value != 0 ? up/fabs(_value) : std::numeric_limits<double>::quiet_NaN();
       return {dn,up};
     }
 
     /// @brief Convenience alias for relErrDownUp
     std::pair<double,double> relErr(const std::string& source = "") const {
       return relErrDownUp(source);
     }
 
     /// @brief The relative negative error with respect to the central value
     double relErrDown(const std::string& source = "") const {
       return relErr(source).first;
     }
 
     /// @brief The relative positive error with respect to the central value
     double relErrUp(const std::string& source = "") const {
       return relErr(source).second;
     }
 
     /// @brief The relative uncertainty from an average of the downward/upward shifts
     /// with respect to the central value
     double relErrAvg(const std::string& source = "") const {
       return _average(relErrDownUp(source));
     }
 
     /// @brief The relative unsigned symmetrised error from the largest component
     /// of the downwards/upwards shifts with respect to the central value
     double relErrEnv(const std::string& source = "") const {
       return _envelope(relErrDownUp(source));
     }
 
     /// @brief The maximal shift of the central value
     /// according to systematic variation @a source
     double valMax(const std::string& source = "") const {
       return val() + errPos(source);
     }
 
     /// @brief The minimal shift of the central value
     /// according to systematic variation @a source
     double valMin(const std::string& source = "") const {
       return val() + errNeg(source);
     }
 
     /// @brief The quadrature sum of uncertainty components
     ///
     /// @note The first element is for the negative error,
     /// the second element for the positive error.
     /// @note In case of one-sided uncertainties, a
     /// zero-deviation is used for the null component.
     /// @note if @a pat_match is supplied, only the subset
     /// of error sources matching the pattern are included.
     std::pair<double,double> quadSum(const std::string& pat_match = "") const noexcept {
       const bool check_match = pat_match != ""; std::smatch match;
       const std::regex re(pat_match, std::regex_constants::icase);
       // interpret as { neg , pos } quad sums
       std::pair<double, double> ret = { 0., 0. };
       for (const auto& item : _error) {
         if (item.first == "")  continue;
         if (check_match && !std::regex_search(item.first, match, re))  continue;
         const auto [dn,up] = _downUp2NegPos(item.second);
         ret.first += dn*dn;
         ret.second += up*up;
       }
       return { - std::sqrt(ret.first), std::sqrt(ret.second) };
     }
 
     /// The negative component of the quad-sum-based total uncertainty components
     double quadSumNeg(const std::string& pat_match = "") const noexcept {
       return quadSum(pat_match).first;
     }
 
     /// The positive component of the quad-sum-based total uncertainty components
     double quadSumPos(const std::string& pat_match = "") const noexcept {
       return quadSum(pat_match).second;
     }
 
     /// @brief The unsigned average of the two quad-sum-based total uncertainty components
     double quadSumAvg(const std::string& pat_match = "") const noexcept {
       return _average(quadSum(pat_match));
     }
 
     /// @brief The unsigned symmetrised uncertainty from the largest
     /// of the quad-sum-based total uncertainty components
     double quadSumEnv(const std::string& pat_match = "") const noexcept {
       return _envelope(quadSum(pat_match));
     }
 
     /// @brief The maximal shift of the central value
     /// from the quad-sum-based total uncertainty components
     double quadSumMax(const std::string& pat_match = "") const {
       return val() + quadSumPos(pat_match);
     }
 
     /// @brief The minimal shift of the central value
     /// from the quad-sum-based total uncertainty components
     double quadSumMin(const std::string& pat_match = "") const {
       return val() + quadSumNeg(pat_match);
     }
 
     /// @brief The total uncertainty
     ///
     /// @note If the user didn't specify a total uncertainty
     /// using the empty string, the sum in quadrature is returned.
     std::pair<double,double> totalErr(const std::string& pat_match = "") const noexcept {
       // check if the user specified the total uncertainty
       if (pat_match == "" && _error.count(""))  return _downUp2NegPos(_error.at(""));
       // otherwise return quadrature sum
       return quadSum(pat_match);
     }
 
     /// @brief The negative total uncertainty
     double totalErrNeg(const std::string& pat_match = "") const noexcept {
       return totalErr(pat_match).first;
     }
 
     /// @brief The positive total uncertainty
     double totalErrPos(const std::string& pat_match = "") const noexcept {
       return totalErr(pat_match).second;
     }
 
     /// @brief The average of the total downward/upward uncertainty shifts,
     /// which are taken to be the signed quadrature sums of the error sources
     double totalErrAvg(const std::string& pat_match = "") const {
       return _average(totalErr(pat_match));
     }
 
     /// @brief The unsigned symmetrised uncertainty from the largest of the
     /// downward/upward shifts of the quad-sum-based total uncertainty
     double totalErrEnv(const std::string& pat_match = "") const {
       return _envelope(totalErr(pat_match));
     }
 
     /// @brief The maximal shift of the central value
     /// from the total uncertainty components
     double totalErrMax(const std::string& pat_match = "") const {
       return val() + totalErrPos(pat_match);
     }
 
     /// @brief The minimal shift of the central value
     /// from the total uncertainty components
     double totalErrMin(const std::string& pat_match = "") const {
       return val() + totalErrNeg(pat_match);
     }
 
     /// @brief The unsigned average total uncertainty
     double valueErr(const std::string& pat_match = "") const {
       return fabs(totalErrAvg(pat_match));
     }
 
     /// @brief The relative negative/positive total uncertainty with respect to the central value
     std::pair<double,double> relTotalErr(const std::string& pat_match = "") const noexcept {
       auto [neg,pos] = totalErr(pat_match);
       neg = _value != 0 ? neg/fabs(_value) : std::numeric_limits<double>::quiet_NaN();
       pos = _value != 0 ? pos/fabs(_value) : std::numeric_limits<double>::quiet_NaN();
       return {neg,pos};
     }
 
     /// The relative negative total uncertainty on the central value
     double relTotalErrNeg(const std::string& pat_match = "") const noexcept {
       return relTotalErr(pat_match).first;
     }
 
     /// The relative positive total uncertainty on the central value
     double relTotalErrPos(const std::string& pat_match = "") const noexcept {
       return relTotalErr(pat_match).second;
     }
 
     /// @brief The relative average of the quad-sum-based total
     /// downward/upward uncertainty components
     double relTotalErrAvg(const std::string& pat_match = "") const noexcept {
       return _average(relTotalErr(pat_match));
     }
 
     /// @brief The unsigned symmetrised error from the largest downward/upward
     /// components of the relative total quad-sum-based uncertainty
     double relTotalErrEnv(const std::string& pat_match = "") const noexcept {
       return _envelope(relTotalErr(pat_match));
     }
 
     /// @brief The list of error source names
     std::vector<std::string> sources() const noexcept {
       std::vector<std::string> keys;
       for (const auto& item : _error)  keys.push_back(item.first);
       return keys;
     }
 
     /// @brief Returns true/false if the error map contains @a key
     bool hasSource(const std::string& key) const noexcept {
       return _error.count(key);
     }
 
     /// @brief The number of error sources in the error map
     size_t numErrs() const noexcept {
       return _error.size();
     }
 
     //@}
 
     /// @name MPI (de-)serialisation
     //@{
 
     size_t _lengthContent(bool fixed_length = false) const noexcept {
       const size_t nErrs = fixed_length? 1 : numErrs();
       return 2*(nErrs + 1);
     }
 
     std::vector<double> _serializeContent(bool fixed_length = false) const noexcept {
       std::vector<double> rtn;
       const size_t nErrs = fixed_length? 1 : numErrs();
       rtn.reserve(2*(nErrs + 1));
       rtn.push_back(_value);
       if (fixed_length) {
         rtn.push_back(1.0);
         rtn.push_back(totalErrNeg());
         rtn.push_back(totalErrPos());
         return rtn;
       }
       rtn.push_back((double)_error.size());
       for (const auto& item : _error) {
         rtn.push_back(item.second.first);
         rtn.push_back(item.second.second);
       }
       return rtn;
     }
 
     void _deserializeContent(const std::vector<double>& data, bool fixed_length = false) {
 
       if (data.size() < 2)
         throw UserError("Length of serialized data should be at least 2!");
 
       if (2*(fixed_length? 1 : data[1]) != (data.size() - 2))
         throw UserError("Expected "+std::to_string(data[1])+" error pairs!");
 
       reset();
       size_t idx = 0;
       auto itr = data.cbegin();
       const auto itrEnd = data.cend();
       while (itr != itrEnd) {
         if (!idx) {
           _value = *itr; ++itr; ++itr;
         }
         else {
           std::string name("source" + std::to_string(idx));
           const double dn = *itr; ++itr;
           const double up = *itr; ++itr;
           setErr({dn, up}, name);
         }
         ++idx;
       }
       if (numErrs() == 1)  renameSource("source1", "");
     }
 
     std::vector<std::string> serializeSources() const noexcept {
       return sources();
     }
 
     void deserializeSources(const std::vector<std::string>& data) {
 
       const size_t nErrs = numErrs();
       if (data.size() != nErrs)
         throw UserError("Expected " + std::to_string(nErrs) + " error source labels!");
 
       for (size_t i = 0; i < data.size(); ++i) {
         const std::string name("source" + std::to_string(i+1));
         if (!_error.count(name))
           throw UserError("Key names have already been updated!");
         renameSource(name, data[i]);
       }
 
     }
 
     // @}
 
   private:
 
     /// @name Utility methods
     //@{
 
     /// @brief Conversion from downward/upward shift to negative/positive error
     ///
     /// @todo Support strategies other than envelope
     std::pair<double,double> _downUp2NegPos(const std::pair<double,double> &e) const noexcept {
       const auto [dn,up] = e;
       if (dn < 0. && up < 0.) {
         // both negative
         const double env = std::min(dn, up);
         return {env, 0.};
       }
       else if (dn < 0.) {
         // dn negative, up positive
         return {dn,up};
       }
       else if (up < 0.) {
         // up negative, dn positive
         return {up, dn};
       }
       // else: both positive
       const double env = std::max(dn, up);
       return {0., env};
     }
 
     /// @brief Work out an average error from an error pair
     double _average(const std::pair<double,double> &err) const noexcept {
       return 0.5*(fabs(err.first) + fabs(err.second));
     }
 
     /// @brief Work out the largest symmetric error from an error pair
     double _envelope(const std::pair<double,double> &err) const noexcept {
       return std::max(fabs(err.first), fabs(err.second));
     }
 
     //@}
 
 
     /// @name Storage
     //@{
 
     /// @brief The central value
     double _value;
 
     /// @brief Error map for signed uncertainties.
     /// The first element of the pair is the effect on the central value
     /// from the systematic downward shift, the second element of the pair
     /// is the effect on the central value from the systematic upward shift.
     ///
     /// @note the empty string is interpreted as the total uncertainty
     std::map<std::string, std::pair<double,double>> _error;
 
     //@}
 
   public:
 
     /// @name Helper methods
     //@{
 
     std::string _toString() const noexcept {
       std::string res ="";
       res += ("value="+ std::to_string(_value));
       for (const auto& item : _error) {
         const std::string name = item.first == ""? "user-supplied total " : "";
         res += (", " + name + "error("+ item.first);
         res += (")={ "+ std::to_string(item.second.first));
         res += (", "+ std::to_string(item.second.second)+" }");
       }
       return res;
     }
 
     //@}
   };
 
 
   /// @name Global operators for Estimate objects
   //@{
 
   /// @brief Add two Estimate objects
   inline Estimate operator + (Estimate lhs, const Estimate& rhs) {
     lhs += rhs;
     return lhs;
   }
 
   /// @brief Add two Estimate objects
   inline Estimate operator + (Estimate lhs, Estimate&& rhs) {
     lhs += std::move(rhs);
     return lhs;
   }
 
   /// @brief Subtract two Estimate objects
   inline Estimate operator - (Estimate lhs, const Estimate& rhs) {
     lhs -= rhs;
     return lhs;
   }
 
   /// @brief Subtract two Estimate objects
   inline Estimate operator - (Estimate lhs, Estimate&& rhs) {
     lhs -= std::move(rhs);
     return lhs;
   }
 
   /// @brief Divide two Estimate objects
   inline Estimate divide(const Estimate& numer, const Estimate& denom,
                          const std::string& pat_uncorr="^stat|^uncor") {
 
     Estimate rtn;
     if (denom.val())  rtn.setVal(numer.val() / denom.val());
     const double newVal = rtn.val();
 
     // get error sources
     std::vector<std::string> sources = numer.sources();
     std::vector<std::string> tmp = denom.sources();
     sources.insert(std::end(sources),
                    std::make_move_iterator(std::begin(tmp)),
                    std::make_move_iterator(std::end(tmp)));
     std::sort(sources.begin(), sources.end());
     sources.erase( std::unique(sources.begin(), sources.end()), sources.end() );
 
     std::smatch match;
     const std::regex re(pat_uncorr, std::regex_constants::icase);
     for (const std::string& src : sources) {
       if (std::regex_search(src, match, re)) {
         // treat as uncorrelated between AOs:
         // add relative errors in quadrature
         double n_dn = 0.0, n_up = 0.0;
         if (numer.hasSource(src)) {
           n_dn = numer.relErrDown(src);
           n_up = numer.relErrUp(src);
         }
         double d_dn = 0.0, d_up = 0.0;
         if (denom.hasSource(src)) {
           d_dn = denom.relErrDown(src);
           d_up = denom.relErrUp(src);
         }
         const double new_dn = fabs(newVal) * std::sqrt(n_dn*n_dn + d_dn*d_dn);
         const double new_up = fabs(newVal) * std::sqrt(n_up*n_up + d_up*d_up);
         rtn.setErr({-new_dn,new_up}, src);
       }
       else {
         // treat as correlated between AOs:
         // work out correlated ratio: R+dR = (N+dN)/(D+dD)
         double n_dn = numer.val(), n_up = numer.val();
         if (numer.hasSource(src)) {
           n_dn += numer.errDown(src);
           n_up += numer.errUp(src);
         }
         double d_dn = denom.val(), d_up = denom.val();
         if (denom.hasSource(src)) {
           d_dn += denom.errDown(src);
           d_up += denom.errUp(src);
         }
         double new_dn = std::numeric_limits<double>::quiet_NaN();
         double new_up = std::numeric_limits<double>::quiet_NaN();
         if (d_dn)  new_dn = n_dn / d_dn - newVal;
         if (d_up)  new_up = n_up / d_up - newVal;
         rtn.setErr({new_dn, new_up}, src);
       }
     }
 
     return rtn;
   }
 
   /// @brief Divide two Estimate objects
   inline Estimate operator / (const Estimate& numer, const Estimate& denom) {
     return divide(numer, denom);
   }
 
   /// @brief Divide two Estimate objects
   inline Estimate operator / (Estimate&& numer, const Estimate& denom) {
     return divide(std::move(numer), denom);
   }
 
   /// @brief Divide two Estimate objects
   inline Estimate operator / (const Estimate& numer, Estimate&& denom) {
     return divide(numer, std::move(denom));
   }
 
   /// @brief Divide two Estimate objects
   inline Estimate operator / (Estimate&& numer, Estimate&& denom) {
     return divide(std::move(numer), std::move(denom));
   }
 
   /// @brief Divide two Estimate objects using binomial statistics
   inline Estimate efficiency(const Estimate& accepted, const Estimate& total,
                              const std::string& pat_uncorr="^stat|^uncor") {
 
     /// Check that the numerator is consistent with being a subset of the denominator
     ///
     /// @todo Need to check that bins are all positive? Integral could be zero due to large +ve/-ve in different bins :O
     const double acc_val = accepted.val();
     const double tot_val = total.val();
     if (acc_val > tot_val)
       throw UserError("Attempt to calculate an efficiency when the numerator is not a subset of the denominator: "
                       + Utils::toStr(acc_val) + " / " + Utils::toStr(tot_val));
 
     Estimate rtn = divide(accepted, total, pat_uncorr);
 
     // get error sources
     std::vector<std::string> sources = accepted.sources();
     std::vector<std::string> tmp = total.sources();
     sources.insert(std::end(sources),
                    std::make_move_iterator(std::begin(tmp)),
                    std::make_move_iterator(std::end(tmp)));
     std::sort(sources.begin(), sources.end());
     sources.erase( std::unique(sources.begin(), sources.end()), sources.end() );
 
     // set binomial error for uncorrelated error sources
     std::smatch match;
     const double eff = rtn.val();
     const std::regex re(pat_uncorr, std::regex_constants::icase);
     for (const std::string& src : sources) {
       double err = std::numeric_limits<double>::quiet_NaN();
       if (!tot_val) {
         rtn.setErr({-err,err}, src);
         continue;
       }
       else if (std::regex_search(src, match, re)) {
         const double acc_err = accepted.relTotalErrAvg();
         const double tot_err = total.totalErrAvg();
         err = sqrt(std::abs( ((1-2*eff)*sqr(acc_err) + sqr(eff)*sqr(tot_err)) / sqr(tot_val) ));
         rtn.setErr({-err,err}, src);
         continue;
       }
     }
 
     return rtn;
   }
 
   //@}
 
 }
 
 #endif

This page was automatically generated by the 2.3.7 LXR engine. The LXR team