Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:25:11

0001 #ifndef ErrorPropogationTypes_h
0002 #define ErrorPropogationTypes_h
0003 
0004 #include "boost/operators.hpp"
0005 #include <cmath>
0006 
0007 class count_t
0008     : boost::addable<count_t,
0009                      boost::incrementable<count_t, boost::totally_ordered<count_t, boost::equivalent<count_t> > > > {
0010 private:
0011   unsigned long count;
0012 
0013 public:
0014   count_t() : count(0) {}
0015   count_t(unsigned long n) : count(n) {}
0016   bool operator<(const count_t& R) const { return count < R.count; }
0017   count_t operator++() {
0018     ++count;
0019     return *this;
0020   }
0021   count_t operator+=(const count_t& R) {
0022     count += R.count;
0023     return *this;
0024   }
0025   unsigned long operator()() const { return count; }
0026   unsigned long error2() const { return count; }
0027   double error() const { return sqrt(count); }
0028   double relative_error() const { return 1 / sqrt(count); }
0029 };
0030 
0031 template <class charT, class traits>
0032 inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& strm, const count_t& f) {
0033   strm << f() << "(" << f.error() << ")";
0034   return strm;
0035 }
0036 
0037 template <class T>
0038 class stats_t
0039     : boost::arithmetic1<
0040           stats_t<T>,
0041           boost::arithmetic2<
0042               stats_t<T>,
0043               T,
0044               boost::partially_ordered<
0045                   stats_t<T>,
0046                   boost::partially_ordered<
0047                       stats_t<T>,
0048                       count_t,
0049                       boost::equality_comparable<stats_t<T>, boost::equality_comparable<stats_t<T>, count_t> > > > > > {
0050 private:
0051   T value, err2;
0052 
0053 public:
0054   stats_t() : value(0), err2(1) {}
0055   stats_t(count_t c) : value(c()), err2(c.error2()) {}
0056   stats_t(T q, T e2) : value(q), err2(e2) {}
0057   static stats_t from_relative_uncertainty2(T q, T re2) { return stats_t(q, q * q * re2); }
0058 
0059   bool operator<(const stats_t& R) const { return value < R.value; }
0060   bool operator<(const count_t& R) const { return value < R(); }
0061   bool operator==(const stats_t& R) const { return value == R.value && err2 == R.err2; }
0062   bool operator==(const count_t& R) const { return value == R() && err2 == R.error2(); }
0063 
0064   stats_t operator+=(const stats_t& R) {
0065     value += R.value;
0066     err2 += R.err2;
0067     return *this;
0068   }
0069   stats_t operator-=(const stats_t& R) {
0070     value -= R.value;
0071     err2 += R.err2;
0072     return *this;
0073   }
0074   stats_t operator*=(const stats_t& R) {
0075     err2 = R.err2 * value * value + err2 * R.value * R.value;
0076     value *= R.value;
0077     return *this;
0078   }
0079   stats_t operator/=(const stats_t& R) {
0080     value /= R.value;
0081     err2 = (err2 + value * value * R.err2) / (R.value * R.value);
0082     return *this;
0083   }
0084 
0085   stats_t operator+=(const T& r) {
0086     value += r;
0087     return *this;
0088   }
0089   stats_t operator-=(const T& r) {
0090     value -= r;
0091     return *this;
0092   }
0093   stats_t operator*=(const T& r) {
0094     value *= r;
0095     err2 *= r * r;
0096     return *this;
0097   }
0098   stats_t operator/=(const T& r) {
0099     value /= r;
0100     err2 /= r * r;
0101     return *this;
0102   }
0103 
0104   stats_t inverse() const { return stats_t(1. / value, err2 / pow(value, 4)); }
0105 
0106   T operator()() const { return value; }
0107   T error2() const { return err2; }
0108   T error() const { return sqrt(err2); }
0109   T relative_error() const { return sqrt(err2) / value; }
0110   T sigmaFrom(const T& x) const { return fabs(value - x) / error(); }
0111 };
0112 
0113 template <class charT, class traits, class T>
0114 inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& strm, const stats_t<T>& f) {
0115   strm << f() << "(" << f.error() << ")";
0116   return strm;
0117 }
0118 
0119 #endif