UnitIO.h

// -*- C++ -*-
//
// UnitIO.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_UnitIO_H
#define ThePEG_UnitIO_H
// This is the declaration of the IUnit and OUnit classes and
// associated templated functions.
 
#include <complex>
#include <iomanip>
#include <sstream>
#include <cstdlib>
#include <cmath>
 
namespace ThePEG {
 
template <typename T, typename UT>
struct OUnit {
 
  OUnit(const T & t, const UT & u): theX(t), theUnit(u) {}
 
  OUnit(const OUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}
 
  const T & theX;
 
  const UT & theUnit;
};
 
template <typename T, typename UT>
struct IUnit {
 
  IUnit(T & t, const UT & u): theX(t), theUnit(u) {}
 
  IUnit(const IUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}
 
  T & theX;
 
  const UT & theUnit;
 
};
 
template <typename T, typename UT>
inline OUnit<T,UT> ounit(const T & t, const UT & ut) {
  return OUnit<T,UT>(t, ut);
}
 
template <typename T, typename UT>
inline IUnit<T,UT> iunit(T & t, const UT & ut) {
  return IUnit<T,UT>(t, ut);
}
 
template <typename OStream, typename T, typename UT>
void ounitstream(OStream & os, const T & t, UT & u) {
  os << t/u;
}
 
template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, T & t, UT & u) {
  double d;
  is >> d;
  t = d*u;;
}
 
template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, std::complex<T> & t, UT & u) {
  std::complex<double> d;
  is >> d;
  t = d*u;;
}
 
template <typename OStream, typename T, typename UT>
OStream & operator<<(OStream & os, const OUnit<T,UT> & u) {
  ounitstream(os, u.theX, u.theUnit);
  return os;
}
 
template <typename IStream, typename T, typename UT>
IStream & operator>>(IStream & is, const IUnit<T,UT> & u) {
  iunitstream(is, u.theX, u.theUnit);
  return is;
}
 
template <typename T, typename UT>
struct OUnitErr {
 
  OUnitErr(const T & t, const T & dt, const UT & u): x(t/u), dx(dt/u) {}
 
  double x;
 
  double dx;
  
};
 
template <typename T, typename UT>
inline OUnitErr<T,UT> ouniterr(const T & t, const T & dt, const UT & ut) {
  return OUnitErr<T,UT>(t, dt, ut);
}
 
inline OUnitErr<double,double> ouniterr(double t, double dt) {
  return OUnitErr<double,double>(t, dt, 1.0);
}
 
template <typename OStream, typename T, typename UT>
OStream & operator<<(OStream & os, const OUnitErr<T,UT> & u) {
  if ( ! isfinite(u.x) ) return os << u.x;
  if ( ! isfinite(u.dx) ) {
    ostringstream out;
    out << u.x << '(' << u.dx << ')';
    return os << out.str();
  }
  double dx = min(u.dx, abs(u.x));
  if ( dx <= 0.0 ) return os << u.x;
  double x = abs(u.x);
  ostringstream osse;
  osse << std::scientific << setprecision(0) << dx;
  string sse = osse.str();
  string::size_type ee = sse.find('e');
  long m = static_cast<long>(round(abs(x)/std::pow(10.0,std::atoi(sse.substr(ee + 1).c_str()))));
  int powx = m <= 0? os.precision(): int(log10(double(m)));
  if ( m <= 0 || powx > os.precision() ) sse[0]='0';  
  ostringstream oss;
  oss << std::scientific << setprecision(powx) << x;
  string ss = oss.str();
  string::size_type e = ss.find('e');
  ostringstream out;
  int pp = std::atoi(ss.substr(e + 1).c_str());
  if ( pp%3 == 0 )
    out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
  else if ( (pp - 1)%3 == 0 ) {
    ostringstream oss;
    oss << std::scientific << setprecision(powx) << x/10.0;
    string ss = oss.str();
    string::size_type e = ss.find('e');
    if ( powx == 0 )
      out << ss.substr(0, e) << "0(" << sse[0] << "0)" << ss.substr(e);
    else if ( powx == 1 )
      out << ss.substr(0, ss.find('.'))
          << ss.substr(ss.find('.') + 1, e - ss.find('.') - 1)
          << "(" << sse[0] << ")" << ss.substr(e);
    else {
      swap(ss[ss.find('.')], ss[ss.find('.') + 1]);
      out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
    }
  }
  else {
    ostringstream oss;
    oss << std::scientific << setprecision(powx) << x*10.0;
    string ss = oss.str();
    string::size_type e = ss.find('e');
    if ( powx == 0 )
      out << "0." << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
    else {
      swap(ss[ss.find('.')], ss[ss.find('.') - 1]);
      out << ss.substr(0, ss.find('.')) << "0" << ss.substr(ss.find('.'), e)
          << "(" << sse[0] << ")" << ss.substr(e);
    }
  }
  string res = out.str();
  if ( u.x < 0.0 )
    res = "-" + res;
  return os << res;
}
 
template <typename T, typename UT>
struct IUnitErr {
 
  IUnitErr(T & t, T & dt, const UT & u): x(t), dx(dt), ut(u) {}
 
  T & x;
 
  T & dx;
  
  UT ut;
 
};
 
template <typename T, typename UT>
inline IUnitErr<T,UT> iuniterr(T & t, T & dt, const UT & ut) {
  return IUnitErr<T,UT>(t, dt, ut);
}
 
inline IUnitErr<double,double> iuniterr(double & t, double & dt) {
  return IUnitErr<double,double>(t, dt, 1.0);
}
 
template <typename IStream, typename T, typename UT>
IStream & operator>>(IStream & is, const IUnitErr<T,UT> & u) {
  string s;
  double x = 0.0;
  double dx = 0.0;
  double ex = 1.0;
  is >> s;
  string::size_type open = s.find('(');
  string::size_type close = s.find(')');
  string se = "0";
  string sp = "1";
  double pe = 1.0;
  if ( open != string::npos && close != string::npos ) {
    se = s.substr(open + 1);
    sp += s.substr(close + 1);
    string::size_type dot = s.find('.');
    if ( dot != string::npos && dot < open ) pe = std::pow(10.0, 1.0 - (open - dot));
  }
 
  istringstream(s) >> x;
  istringstream(se) >> dx;
  istringstream(sp) >> ex;
 
  u.x = x*ex*u.ut;
  u.dx = dx*ex*pe*u.ut;
 
  return is;
}
 
}
 
#endif /* ThePEG_UnitIO_H */