tkmeasure.h

/********************************************************************************
 *   Copyright (c) : Université de Lyon 1, CNRS/IN2P3, UMR5822,                 *
 *                   IP2I, F-69622 Villeurbanne Cedex, France                   *
 *                   Normandie Université, ENSICAEN, UNICAEN, CNRS/IN2P3,       *
 *                   LPC Caen, F-14000 Caen, France                             *
 *   Contibutor(s) :                                                            *
 *      Jérémie Dudouet jeremie.dudouet@cnrs.fr [2020]                          *
 *      Diego Gruyer    diego.gruyer@cnrs.fr    [2020]                          *
 *                                                                              *
 *   Licensed under the MIT License <http://opensource.org/licenses/MIT>.       *
 *   SPDX-License-Identifier: MIT                                               *
 ********************************************************************************/
 
#ifndef tkmeasure_H
#define tkmeasure_H
 
#include <cmath>
#include <iomanip>
 
#include "tkn_config.h"
#include "tklog.h"
#include "tkproperty.h"
#include "tkunit.h"
 
#ifdef HAS_ROOT
#include "TClass.h"
#endif
 
using namespace std;
 
namespace tkn {
 
typedef std::map<tkstring, pair<tkstring, tkstring>> properties; // map properties, value, unit

class tkmeasure : public tkproperty {
protected:
   tkunit fValue{};     // the value of the measure
   tkunit fError{};     // the error on this measure
   tkunit fErrorLow{};  // if define, assym low error value
   tkunit fErrorHigh{}; // if define, assym high error value
 
   tkstring ftype = ""; // type of measured data
 
   bool fAsymError = false;
 
   tkstring finfo_tag{}; // LT(<), GT(>), LE(<=), GE(>=), AP(~), CA (calculated), SY (from systematics), ? (uncertain, see comments), STABLE
 
public:
   tkmeasure() { ; }
   tkmeasure(double _value, const tkstring &_unit_name, double _err = 0) : tkmeasure(_value, gunits->get_key(_unit_name), _err) {}
   tkmeasure(double _value, const tkstring &_unit_name, double _errlow, double _errhigh) : tkmeasure(_value, gunits->get_key(_unit_name), _errlow, _errhigh) {}
   tkmeasure(double _value, const tkunit_manager::units_keys &_unit, double _err = 0);
   tkmeasure(double _value, const tkunit_manager::units_keys &_unit, double _errlow, double _errhigh);
   virtual ~tkmeasure() override = default;

   bool is_asym_errors() const { return fAsymError; }

   void set_info_tag(const tkstring &_tag);
   const tkstring &get_info_tag() const { return finfo_tag; }
   const tkstring &get_info_tag_const() const { return finfo_tag; }

   bool has_alternative() const { return finfo_tag.contains(";ALT="); }
   tkstring get_alternative() const
   {
      const tkstring marker = ";ALT=";
      const auto position = finfo_tag.find(marker);
      if (position == tkstring::npos) return "";
      const auto begin = position + marker.size();
      const auto end = finfo_tag.find(';', begin);
      return finfo_tag.substr(begin, end == tkstring::npos ? tkstring::npos : end - begin);
   }

   bool has_incompatible_alternative() const
   { return finfo_tag.contains(";ALT_INCOMPATIBLE=1"); }

   void set_type(const tkstring &_type) { ftype = _type; }
   const tkstring &get_type() const { return ftype; }

   double get_value() const { return fValue.get_value(); }
   double get_value(const tkunit_manager::units_keys &_unit) const;
   double get_value(const tkstring &_unit_name) const { return get_value(gunits->get_key(_unit_name)); }

   void set(const double &_val, const tkstring &_unit = "");
   void set(const double &_val, const tkunit_manager::units_keys &_unit);

   virtual void set_error(const double &_err)
   {
      fError.set_value(_err);
      fErrorLow.clear();
      fErrorHigh.clear();
      fAsymError = false;
   }

   virtual void set_error(const double &_err_low, const double &_err_high)
   {
      fErrorLow.set_value(_err_low);
      fErrorHigh.set_value(_err_high);
      fError.clear();
      fAsymError = true;
   }

   double get_error() const;
   double get_error(const tkunit_manager::units_keys &_unit) const;
   double get_error(const tkstring &_unit_name) const { return get_error(gunits->get_key(_unit_name)); }

   double get_error_low() const;
   double get_error_low(const tkunit_manager::units_keys &_unit) const;
   double get_error_low(const tkstring &_unit_name) const { return get_error_low(gunits->get_key(_unit_name)); }

   double get_error_high() const;
   double get_error_high(const tkunit_manager::units_keys &_unit) const;
   double get_error_high(const tkstring &_unit_name) const { return get_error_high(gunits->get_key(_unit_name)); }

   tkstring get_unit() const { return fValue.get_unit_name(); }
   tkunit_manager::units_keys get_unit_key() const { return fValue.get_unit_key(); }

   bool set_unit(const tkstring &_unit_name) { return set_unit(gunits->get_key(_unit_name)); }
   bool set_unit(const tkunit_manager::units_keys &_unit);

   friend std::ostream &operator<<(std::ostream &os, const tkmeasure *_measure)
   {
 
      int save_precision = os.precision();
      os << left << _measure->get_type();
 
      // check energy offset
      tkstring offset = "";
      if (_measure->get_info_tag_const().contains(";OFF=")) {
         int idx = _measure->get_info_tag_const().index(";OFF=");
         offset = _measure->get_info_tag_const().substr(idx + 5, 1) + "+";
      }
 
      if (_measure->is_info(kLimit)) {
         if (_measure->get_info_tag_const().contains(";ERR=LT"))
            os << " < ";
         else if (_measure->get_info_tag_const().contains(";ERR=LE"))
            os << " <= ";
         else if (_measure->get_info_tag_const().contains(";ERR=GT"))
            os << " > ";
         else if (_measure->get_info_tag_const().contains(";ERR=GE"))
            os << " >= ";
      } else if (_measure->is_info(kAbout))
         os << " ~ ";
      else
         os << " = ";
 
      if (!_measure->get_info_tag_const().empty() && _measure->get_info_tag_const().contains("STABLE")) {
         os << left << std::setw(7) << "STABLE" << std::setw(9) << "";
      } else {
         if (_measure->is_asym_errors() && _measure->get_error_low() > 0 && _measure->get_error_high() > 0) {
            int exp_error_low = _measure->get_exp_error_precision(_measure->get_error_low());
            int exp_error_high = _measure->get_exp_error_precision(_measure->get_error_high());
            int exp_error = max(max(exp_error_low, exp_error_high), _measure->get_exp_error_precision(_measure->get_value()));
 
            if (_measure->get_value() != 0 && (_measure->get_value() > 1e5 || _measure->get_value() < 1e-5))
               os << scientific;
            else
               os << fixed;
 
            os << left << setprecision(exp_error) << offset << std::setw(7) << _measure->get_value();
            os << " (-" << std::setw(5) << _measure->get_error_low() << " ; +" << std::setw(5) << _measure->get_error_high() << ") " << _measure->get_unit();
         } else {
            int exp_error = max(_measure->get_exp_error_precision(_measure->get_value()), _measure->get_exp_error_precision(_measure->get_error()));
            if (_measure->get_value() != 0 && (abs(_measure->get_value()) > 1e5 || abs(_measure->get_value()) < 1e-5))
               os << scientific;
            else
               os << fixed;
 
            if (_measure->get_error() > 0.) {
               os << left << setprecision(exp_error) << offset << std::setw(7) << _measure->get_value();
               os << " (" << std::setw(5) << _measure->get_error() << ") " << _measure->get_unit();
            } else if (_measure->get_error() == 0.) {
               os << left << setprecision(_measure->get_exp_error_precision(_measure->get_value())) << offset << std::setw(7) << _measure->get_value();
               os << _measure->get_unit();
            } else {
               os << left << setprecision(_measure->get_exp_error_precision(_measure->get_value())) << offset << std::setw(7) << _measure->get_value();
               if (_measure->get_value() == 0.)
                  os << "  " << std::setw(5) << "     " << "  " << _measure->get_unit();
               else
                  os << fixed << left << setprecision(exp_error) << " (" << std::setw(5) << "  ?  " << ") " << _measure->get_unit();
            }
         }
      }
 
      if (!_measure->is_info(kKnown) || _measure->is_converted()) {
         os << " [" << _measure->get_info_str(false) << "]";
      }
 
      os.precision(save_precision);
 
      return os;
   }

   void print(bool with_return = true)
   {
      cout << (this);
      if (with_return) cout << endl;
   }
 
   tkmeasure operator+(const tkmeasure &_measure) const
   {
      tkmeasure result = *this;
      tkmeasure measure_to_add = _measure;
      if (this->get_unit() != measure_to_add.get_unit()) {
         if (!measure_to_add.set_unit(get_unit())) {
            glog << warning_o << "trying to add two different types of units... ignored" << do_endl;
            return result;
         }
      }
      result.set(this->get_value() + measure_to_add.get_value());
      double low_e1, high_e1, low_e2, high_e2;
      if (is_asym_errors()) {
         low_e1 = get_error_low();
         high_e1 = get_error_high();
      } else {
         low_e1 = get_error();
         high_e1 = get_error();
      }
      if (measure_to_add.is_asym_errors()) {
         low_e2 = measure_to_add.get_error_low();
         high_e2 = measure_to_add.get_error_high();
      } else {
         low_e2 = measure_to_add.get_error();
         high_e2 = measure_to_add.get_error();
      }
      if (low_e1 < 0. || high_e1 < 0. || low_e2 < 0. || high_e2 < 0.) {
         // Unknown uncertainty is not a numeric standard deviation.
         result.set_error(-1.);
      } else if (is_asym_errors() || measure_to_add.is_asym_errors()) {
         result.set_error(sqrt(low_e1 * low_e1 + low_e2 * low_e2), sqrt(high_e1 * high_e1 + high_e2 * high_e2));
      } else
         result.set_error(sqrt(low_e1 * low_e1 + low_e2 * low_e2));
 
      if (get_info_const() != kKnown || measure_to_add.get_info() != kKnown) {
         if (get_info_const() == kKnown && measure_to_add.get_info() != kKnown)
            result.set_info(measure_to_add.get_info());
         else if (get_info_const() != kKnown && measure_to_add.get_info() == kKnown)
            result.set_info(get_info_const());
         if (get_info_const() != kKnown && measure_to_add.get_info() != kKnown) result.set_info(kUncertain);
 
         result.set_info_tag(get_info_tag_const() + measure_to_add.get_info_tag());
         if (measure_to_add.is_converted()) result.set_converted();
      }
      return result;
   }
   tkmeasure operator*(double _multiplier) const
   {
      tkmeasure result = *this;
      result.set(get_value() * _multiplier);
      double low_e1, high_e1;
      const double scale = std::abs(_multiplier);
      if (is_asym_errors()) {
         if (_multiplier < 0) {
            low_e1 = get_error_high() * scale;
            high_e1 = get_error_low() * scale;
         } else {
            low_e1 = get_error_low() * scale;
            high_e1 = get_error_high() * scale;
         }
      } else {
         low_e1 = get_error() * scale;
         high_e1 = get_error() * scale;
      }
      if (is_asym_errors())
         result.set_error(low_e1, high_e1);
      else
         result.set_error(low_e1);
 
      if (_multiplier < 0) {
         const tkstring tag = result.get_info_tag();
         if (tag.contains(";ERR=LT"))
            result.set_info_tag(tag.copy().replace_all(";ERR=LT", ";ERR=GT"));
         else if (tag.contains(";ERR=LE"))
            result.set_info_tag(tag.copy().replace_all(";ERR=LE", ";ERR=GE"));
         else if (tag.contains(";ERR=GT"))
            result.set_info_tag(tag.copy().replace_all(";ERR=GT", ";ERR=LT"));
         else if (tag.contains(";ERR=GE"))
            result.set_info_tag(tag.copy().replace_all(";ERR=GE", ";ERR=LE"));
      }
      return result;
   }
   friend tkmeasure operator*(double _multiplier, const tkmeasure &obj) { return obj * _multiplier; }
   tkmeasure operator-() const
   {
      tkmeasure result = *this * -1.;
      return result;
   }
   tkmeasure operator-(const tkmeasure &_measure) const
   {
      tkmeasure measure_to_add = -_measure;
      tkmeasure result = *this + measure_to_add;
      return result;
   }

   int get_exp_error_precision(const double &_err) const
   {
      double precision_float_error = tkstring::get_precision(tkstring::Form("%g", _err));
      if (_err < 0.) precision_float_error = tkstring::get_precision(tkstring::Form("%g", get_value()));
      int precision_int_err = (int)precision_float_error;
      if (precision_float_error < 1) precision_int_err = (int)(-(1. / precision_float_error + 0.5));
 
      int exp_error = (int)std::floor(std::log10(std::fabs(precision_int_err)));
      if (precision_int_err >= 0) exp_error = 0;
      return exp_error;
   }
 
#ifdef HAS_ROOT
   ClassDefOverride(tkmeasure, 0);
#endif
};
} // namespace tkn
 
#endif