ThePEG/Utilities/UnitIO.h

0001 // -*- C++ -*-
0002 //
0003 // UnitIO.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef ThePEG_UnitIO_H
0010 #define ThePEG_UnitIO_H
0011 // This is the declaration of the IUnit and OUnit classes and
0012 // associated templated functions.
0013
0014 #include <complex>
0015 #include <iomanip>
0016 #include <sstream>
0017 #include <cstdlib>
0018 #include <cmath>
0019
0020 namespace ThePEG {
0021
0022 /**
0023  * The OUnit< class is used to
0024  * facilitate output of unitful numbers to a
0025  * persistent stream. An Energy can hence be written like
0026  * this:<BR> <code>os
0027  * << ounit(x, GeV);</code><BR> Also containers of unitful
0028  * numbers may be written like this, as well as LorentzVector and
0029  * ThreeVector.
0030  *
0031  * @see PersistentOStream
0032  * @see PersistentIStream
0033  *
0034  */
0035 template <typename T, typename UT>
0036 struct OUnit {
0037
0038   /** Constructor given an object to be written assuming the given
0039    *  unit. */
0040   OUnit(const T & t, const UT & u): theX(t), theUnit(u) {}
0041
0042   /** Copy constructor */
0043   OUnit(const OUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}
0044
0045   /** Reference to the object to be written. */
0046   const T & theX;
0047
0048   /** The unit assumed when writing the object. */
0049   const UT & theUnit;
0050 };
0051
0052 /**
0053  * The IUnit class is used to facilitate input of unitful numbers from
0054  * and to a persistent stream. An Energy can hence be read like
0055  * this:<BR> <code>is >> iunit(x, GeV);</code><BR> Also containers of
0056  * unitful numbers may be read like this, as well as LorentzVector and
0057  * ThreeVector.
0058  *
0059  * @see PersistentOStream
0060  * @see PersistentIStream
0061  *
0062  */
0063 template <typename T, typename UT>
0064 struct IUnit {
0065
0066   /** Constructor given an object to be read assuming the given
0067    *  unit. */
0068   IUnit(T & t, const UT & u): theX(t), theUnit(u) {}
0069
0070   /** Copy constructor */
0071   IUnit(const IUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}
0072
0073   /** Reference to the object to be read. */
0074   T & theX;
0075
0076   /** The unit assumed when reading the object. */
0077   const UT & theUnit;
0078
0079 };
0080
0081 /** Helper function creating a OUnit object given an object and a
0082  *  unit. */
0083 template <typename T, typename UT>
0084 inline OUnit<T,UT> ounit(const T & t, const UT & ut) {
0085   return OUnit<T,UT>(t, ut);
0086 }
0087
0088 /** Helper function creating a IUnit object given an object and a
0089  *  unit. */
0090 template <typename T, typename UT>
0091 inline IUnit<T,UT> iunit(T & t, const UT & ut) {
0092   return IUnit<T,UT>(t, ut);
0093 }
0094
0095 /** Helper function writing out an object with a given unit to an
0096  *  output stream. */
0097 template <typename OStream, typename T, typename UT>
0098 void ounitstream(OStream & os, const T & t, UT & u) {
0099   os << t/u;
0100 }
0101
0102 /** Helper function reading an object with a given unit from an
0103  *  input stream. */
0104 template <typename IStream, typename T, typename UT>
0105 void iunitstream(IStream & is, T & t, UT & u) {
0106   double d;
0107   is >> d;
0108   t = d*u;;
0109 }
0110
0111 /** Helper function reading a complex object with a given unit from an
0112  *  input stream. */
0113 template <typename IStream, typename T, typename UT>
0114 void iunitstream(IStream & is, std::complex<T> & t, UT & u) {
0115   std::complex<double> d;
0116   is >> d;
0117   t = d*u;;
0118 }
0119
0120 /** Output an OUnit object to a stream. */
0121 template <typename OStream, typename T, typename UT>
0122 OStream & operator<<(OStream & os, const OUnit<T,UT> & u) {
0123   ounitstream(os, u.theX, u.theUnit);
0124   return os;
0125 }
0126
0127 /** Input an IUnit object from a stream. */
0128 template <typename IStream, typename T, typename UT>
0129 IStream & operator>>(IStream & is, const IUnit<T,UT> & u) {
0130   iunitstream(is, u.theX, u.theUnit);
0131   return is;
0132 }
0133
0134 /**
0135  * OUnitErr is used to write out unitful numbers with an error
0136  * estimate on a standard ostream. using the helper function ouniterr
0137  * an energy <code>e</code> with an error estimate <code>de</code> can
0138  * be written out as eg. <code>cout << ouniterr(e, de,
0139  * GeV);</code>. The result will be presented in scientific format
0140  * (with the exponent divisible by three) with the relevant number of
0141  * significant digits with a single digit in parenthesis indicating
0142  * the error in the least significant digit,
0143  * eg. <code>1.23(2)e+03</code>.
0144  */
0145 template <typename T, typename UT>
0146 struct OUnitErr {
0147
0148   /** Constructor given an object to be written assuming the given
0149    *  unit. */
0150   OUnitErr(const T & t, const T & dt, const UT & u): x(t/u), dx(dt/u) {}
0151
0152   /** The number to be written. */
0153   double x;
0154
0155   /** The estimated error of the number to be written. */
0156   double dx;
0157
0158 };
0159
0160 /** Helper function creating a OUnitErr object. */
0161 template <typename T, typename UT>
0162 inline OUnitErr<T,UT> ouniterr(const T & t, const T & dt, const UT & ut) {
0163   return OUnitErr<T,UT>(t, dt, ut);
0164 }
0165
0166 /** Helper function creating a OUnitErr object. */
0167 inline OUnitErr<double,double> ouniterr(double t, double dt) {
0168   return OUnitErr<double,double>(t, dt, 1.0);
0169 }
0170
0171 /** Output an OUnitErr object to a stream. */
0172 template <typename OStream, typename T, typename UT>
0173 OStream & operator<<(OStream & os, const OUnitErr<T,UT> & u) {
0174   if ( ! isfinite(u.x) ) return os << u.x;
0175   if ( ! isfinite(u.dx) ) {
0176     ostringstream out;
0177     out << u.x << '(' << u.dx << ')';
0178     return os << out.str();
0179   }
0180   double dx = min(u.dx, abs(u.x));
0181   if ( dx <= 0.0 ) return os << u.x;
0182   double x = abs(u.x);
0183   ostringstream osse;
0184   osse << std::scientific << setprecision(0) << dx;
0185   string sse = osse.str();
0186   string::size_type ee = sse.find('e');
0187   long m = static_cast<long>(round(abs(x)/std::pow(10.0,std::atoi(sse.substr(ee + 1).c_str()))));
0188   int powx = m <= 0? os.precision(): int(log10(double(m)));
0189   if ( m <= 0 || powx > os.precision() ) sse[0]='0';
0190   ostringstream oss;
0191   oss << std::scientific << setprecision(powx) << x;
0192   string ss = oss.str();
0193   string::size_type e = ss.find('e');
0194   ostringstream out;
0195   int pp = std::atoi(ss.substr(e + 1).c_str());
0196   if ( pp%3 == 0 )
0197     out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
0198   else if ( (pp - 1)%3 == 0 ) {
0199     ostringstream oss;
0200     oss << std::scientific << setprecision(powx) << x/10.0;
0201     string ss = oss.str();
0202     string::size_type e = ss.find('e');
0203     if ( powx == 0 )
0204       out << ss.substr(0, e) << "0(" << sse[0] << "0)" << ss.substr(e);
0205     else if ( powx == 1 )
0206       out << ss.substr(0, ss.find('.'))
0207       << ss.substr(ss.find('.') + 1, e - ss.find('.') - 1)
0208       << "(" << sse[0] << ")" << ss.substr(e);
0209     else {
0210       swap(ss[ss.find('.')], ss[ss.find('.') + 1]);
0211       out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
0212     }
0213   }
0214   else {
0215     ostringstream oss;
0216     oss << std::scientific << setprecision(powx) << x*10.0;
0217     string ss = oss.str();
0218     string::size_type e = ss.find('e');
0219     if ( powx == 0 )
0220       out << "0." << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
0221     else {
0222       swap(ss[ss.find('.')], ss[ss.find('.') - 1]);
0223       out << ss.substr(0, ss.find('.')) << "0" << ss.substr(ss.find('.'), e)
0224       << "(" << sse[0] << ")" << ss.substr(e);
0225     }
0226   }
0227   string res = out.str();
0228   if ( u.x < 0.0 )
0229     res = "-" + res;
0230   return os << res;
0231 }
0232
0233 /**
0234  * The IUnitErr class is used to facilitate input of unitful numbers
0235  * with error estimates written out using the OUnitErr class.
0236  *
0237  */
0238 template <typename T, typename UT>
0239 struct IUnitErr {
0240
0241   /** Constructor given an object to be read assuming the given
0242    *  unit. */
0243   IUnitErr(T & t, T & dt, const UT & u): x(t), dx(dt), ut(u) {}
0244
0245   /** Reference to the object to be read. */
0246   T & x;
0247
0248   /** The estimated error of the number to be read. */
0249   T & dx;
0250
0251   /** The unit assumed when reading the object. */
0252   UT ut;
0253
0254 };
0255
0256 /** Helper function creating a IUnitErr object. */
0257 template <typename T, typename UT>
0258 inline IUnitErr<T,UT> iuniterr(T & t, T & dt, const UT & ut) {
0259   return IUnitErr<T,UT>(t, dt, ut);
0260 }
0261
0262 /** Helper function creating a OUnitErr object. */
0263 inline IUnitErr<double,double> iuniterr(double & t, double & dt) {
0264   return IUnitErr<double,double>(t, dt, 1.0);
0265 }
0266
0267 /** Input an IUnit object from a stream. */
0268 template <typename IStream, typename T, typename UT>
0269 IStream & operator>>(IStream & is, const IUnitErr<T,UT> & u) {
0270   string s;
0271   double x = 0.0;
0272   double dx = 0.0;
0273   double ex = 1.0;
0274   is >> s;
0275   string::size_type open = s.find('(');
0276   string::size_type close = s.find(')');
0277   string se = "0";
0278   string sp = "1";
0279   double pe = 1.0;
0280   if ( open != string::npos && close != string::npos ) {
0281     se = s.substr(open + 1);
0282     sp += s.substr(close + 1);
0283     string::size_type dot = s.find('.');
0284     if ( dot != string::npos && dot < open ) pe = std::pow(10.0, 1.0 - (open - dot));
0285   }
0286
0287   istringstream(s) >> x;
0288   istringstream(se) >> dx;
0289   istringstream(sp) >> ex;
0290
0291   u.x = x*ex*u.ut;
0292   u.dx = dx*ex*pe*u.ut;
0293
0294   return is;
0295 }
0296
0297 }
0298
0299 #endif /* ThePEG_UnitIO_H */