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]                          *
 *                                                                              *
 *    This software is governed by the CeCILL-B license under French law and    *
 *    abiding by the  rules of distribution of free  software.  You can use,    *
 *    modify  and/ or  redistribute  the  software under  the  terms of  the    *
 *    CeCILL-B license as circulated by CEA, CNRS and INRIA at the following    *
 *    URL \"http://www.cecill.info\".                                           *
 *                                                                              *
 *    As a counterpart to the access  to the source code and rights to copy,    *
 *    modify  and redistribute granted  by the  license, users  are provided    *
 *    only with a limited warranty  and the software's author, the holder of    *
 *    the economic  rights, and the  successive licensors have  only limited    *
 *    liability.                                                                *
 *                                                                              *
 *    In this respect, the user's attention is drawn to the risks associated    *
 *    with loading,  using, modifying  and/or developing or  reproducing the    *
 *    software by the user in light of its specific status of free software,    *
 *    that  may mean that  it is  complicated to  manipulate, and  that also    *
 *    therefore  means that it  is reserved  for developers  and experienced    *
 *    professionals having in-depth  computer knowledge. Users are therefore    *
 *    encouraged  to load  and test  the software's  suitability  as regards    *
 *    their  requirements  in  conditions  enabling the  security  of  their    *
 *    systems  and/or data to  be ensured  and, more  generally, to  use and    *
 *    operate it in the same conditions as regards security.                    *
 *                                                                              *
 *    The fact that  you are presently reading this means  that you have had    *
 *    knowledge of the CeCILL-B license and that you accept its terms.          *
 ********************************************************************************/
 
#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() = default;
 
    bool is_asym_errors() const { return fAsymError; }
 
    void set_info_tag(const tkstring &_tag);
    const tkstring &get_info_tag() {return finfo_tag;}
    const tkstring &get_info_tag_const() const {return finfo_tag;}
 
    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);
    double get_value(const tkstring &_unit_name) {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);
    double get_error(const tkstring &_unit_name) {return get_error(gunits->get_key(_unit_name));}
 
    double get_error_low() const;
    double get_error_low(const tkunit_manager::units_keys &_unit);
    double get_error_low(const tkstring &_unit_name) {return get_error_low(gunits->get_key(_unit_name));}
 
    double get_error_high() const;
    double get_error_high(const tkunit_manager::units_keys &_unit);
    double get_error_high(const tkstring &_unit_name) {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(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());
            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;
        if(is_asym_errors()) {
            if(_multiplier<0) {
                low_e1 = get_error_high()*_multiplier;
                high_e1 = get_error_low()*_multiplier;
            }
            else {
                low_e1 = get_error_low()*_multiplier;
                high_e1 = get_error_high()*_multiplier;
            }
        }
        else {
            low_e1 = get_error()*_multiplier;
            high_e1 = get_error()*_multiplier;
        }
        if(is_asym_errors()) result.set_error(low_e1,high_e1);
        else result.set_error(low_e1);
 
        if(_multiplier<0) {
            if(result.get_info_tag().contains("LT")) result.set_info_tag(result.get_info_tag().copy().replace_all("LT","GT"));
            if(result.get_info_tag().contains("LE")) result.set_info_tag(result.get_info_tag().copy().replace_all("LE","GE"));
            if(result.get_info_tag().contains("GT")) result.set_info_tag(result.get_info_tag().copy().replace_all("GT","LT"));
        }
        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
    ClassDef(tkmeasure,0);
#endif
 
};
}
 
#endif