diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -1,204 +1,204 @@ ## Process this file with autoconf to produce a configure script. AC_PREREQ(2.59) AC_INIT([YODA],[1.8.0],[yoda@projects.hepforge.org],[YODA]) ## Check and block installation into the src/build dir if test "$prefix" = "$PWD"; then AC_MSG_ERROR([Installation into the build directory is not supported: use a different --prefix argument]) fi ## Force default prefix to have a path value rather than NONE if test "$prefix" = "NONE"; then prefix=/usr/local fi AC_CONFIG_SRCDIR([src/Counter.cc]) AM_INIT_AUTOMAKE([subdir-objects -Wall dist-bzip2 1.10]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) m4_ifdef([AM_PROG_AR], [AM_PROG_AR]) AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_HEADERS([include/YODA/Config/DummyConfig.h include/YODA/Config/YodaConfig.h include/YODA/Config/BuildConfig.h]) AC_DEFINE_UNQUOTED(YODA_VERSION, "$PACKAGE_VERSION", "YODA version string") AC_DEFINE_UNQUOTED(YODA_NAME, "$PACKAGE_NAME", "YODA name string") AC_DEFINE_UNQUOTED(YODA_STRING, "$PACKAGE_STRING", "YODA name and version string") AC_DEFINE_UNQUOTED(YODA_TARNAME, "$PACKAGE_TARNAME", "YODA short name string") AC_DEFINE_UNQUOTED(YODA_BUGREPORT, "$PACKAGE_BUGREPORT", "YODA contact email address") ## OS X AC_CEDAR_OSX ## Set default compiler flags if test "x$CXXFLAGS" = "x"; then CXXFLAGS="-O3"; fi ## Compiler setup AC_LANG(C++) AC_PROG_CXX AX_CXX_COMPILE_STDCXX([11], [noext], [mandatory]) AC_PROG_INSTALL AC_PROG_LN_S AC_DISABLE_STATIC AC_LIBTOOL_DLOPEN AC_PROG_LIBTOOL ## Work out library suffix for the build LIB_SUFFIX=\\\"$shrext_cmds\\\" AC_SUBST([LIB_SUFFIX]) ## Set default build flags AC_CEDAR_CHECKCXXFLAG([-pedantic], [AM_CXXFLAGS="$AM_CXXFLAGS -pedantic"]) AC_CEDAR_CHECKCXXFLAG([-Wall], [AM_CXXFLAGS="$AM_CXXFLAGS -Wall -Wno-format"]) dnl AC_CEDAR_CHECKCXXFLAG([-std=c++98], [AM_CXXFLAGS="$AM_CXXFLAGS -std=c++98"]) dnl AC_CEDAR_CHECKCXXFLAG([-Wno-unused-variable], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-unused-variable"]) ## Debug flag (default=none) AC_ARG_ENABLE([debug], [AC_HELP_STRING(--enable-debug, [build with debugging symbols @<:@default=no@:>@])], [], [enable_debug=no]) if test x$enable_debug = xyes; then [AM_CXXFLAGS="$AM_CXXFLAGS -g"] fi ## Optional zlib support for gzip-compressed data streams/files AX_CHECK_ZLIB ## Optional ROOT compatibility AC_ARG_ENABLE([root], [AC_HELP_STRING(--disable-root, [don't try to build YODA interface to PyROOT (needs root-config) @<:@default=yes@:>@])], [], [enable_root=yes]) if test "x$enable_root" = "xyes"; then AC_PATH_PROG(ROOTCONFIG, [root-config]) if test "x$ROOTCONFIG" = "x"; then AC_MSG_WARN([root-config not found -- not building extra ROOT compatibility tools]) enable_root=no; else AC_MSG_CHECKING([ROOT version]) ROOT_VERSION=`$ROOTCONFIG --version` ROOT_MAJOR_VERSION=`echo $ROOT_VERSION | cut -d. -f1` ROOT_MINOR_VERSION=`echo $ROOT_VERSION | cut -d. -f2 | cut -d/ -f1` ROOT_MICRO_VERSION=`echo $ROOT_VERSION | cut -d. -f2 | cut -d/ -f2` AC_MSG_RESULT([$ROOT_VERSION ($ROOT_MAJOR_VERSION,$ROOT_MINOR_VERSION,$ROOT_MICRO_VERSION)]) if test "$ROOT_MAJOR_VERSION" -lt 5; then enable_root=no; AC_MSG_WARN([ROOT major version is < 5 -- not building extra ROOT compatibility tools]) elif test "$ROOT_MAJOR_VERSION" -eq 5; then if test "$ROOT_MINOR_VERSION" -lt 33; then enable_root=no; AC_MSG_WARN([ROOT version is less than 5.33 -- not building extra ROOT compatibility tools]) elif test "$ROOT_MINOR_VERSION" -eq 33 && test "$ROOT_MICRO_VERSION" -lt 2; then enable_root=no; AC_MSG_WARN([ROOT version is less than 5.33/02 -- not building extra ROOT compatibility tools]) fi fi # TODO: Test for existence of TPython, instance_from_void API, etc. #AM_CXXFLAGS="$AM_CXXFLAGS -Wno-long-long" ROOT_CXXFLAGS=`$ROOTCONFIG --cflags` ROOT_LDFLAGS=`$ROOTCONFIG --ldflags` ROOT_LIBS=`$ROOTCONFIG --libs` AC_SUBST(ROOT_CXXFLAGS) AC_SUBST(ROOT_LDFLAGS) AC_SUBST(ROOT_LIBS) fi fi AM_CONDITIONAL(ENABLE_ROOT, [test x$enable_root = xyes]) if test x$enable_root = xyes; then AC_MSG_NOTICE([Building extra ROOT compatibility tools]) else AC_MSG_NOTICE([Not building extra ROOT compatibility tools]) fi ## Python extension AC_ARG_ENABLE(pyext, [AC_HELP_STRING(--disable-pyext, [don't build Python module (default=build)])], [], [enable_pyext=yes]) ## Basic Python checks if test x$enable_pyext = xyes; then AX_PYTHON_DEVEL([>= '2.7.3']) AC_SUBST(PYTHON_VERSION) YODA_PYTHONPATH=`$PYTHON -c "from __future__ import print_function; import distutils.sysconfig; print(distutils.sysconfig.get_python_lib(prefix='$prefix', plat_specific=True));"` AC_SUBST(YODA_PYTHONPATH) if test -z "$PYTHON"; then AC_MSG_ERROR([Can't build Python extension since python can't be found]) enable_pyext=no fi if test -z "$PYTHON_CPPFLAGS"; then AC_MSG_ERROR([Can't build Python extension since Python.h header file cannot be found]) enable_pyext=no fi fi AM_CONDITIONAL(ENABLE_PYEXT, [test x$enable_pyext == xyes]) ## Cython checks if test x$enable_pyext == xyes; then AM_CHECK_CYTHON([0.24], [:], [:]) if test x$CYTHON_FOUND = xyes; then AC_MSG_NOTICE([Cython >= 0.24 found: Python extension source can be rebuilt (for developers)]) # Force rebuild since we have a sufficient Cython - touch pyext/rivet/core.pyx + test -f $(top_srcdir)/pyext/yoda/core.pyx && touch $(top_srcdir)/yoda/rivet/core.pyx fi - AC_CHECK_FILE([pyext/yoda/core.cpp], + AC_CHECK_FILE([$(top_srcdir)/pyext/yoda/core.cpp], [], [if test "x$CYTHON_FOUND" != "xyes"; then AC_MSG_ERROR([Cython is required for --enable-pyext, no pre-built core.cpp was found.]) fi]) ## Set extra Python extension build flags (to cope with Cython output code oddities) PYEXT_CXXFLAGS=$CXXFLAGS AC_CEDAR_CHECKCXXFLAG([-Wno-unused-but-set-variable], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-unused-but-set-variable"]) AC_CEDAR_CHECKCXXFLAG([-Wno-sign-compare], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-sign-compare"]) AC_CEDAR_CHECKCXXFLAG([-Wno-strict-prototypes], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-strict-prototypes"]) AC_SUBST(PYEXT_CXXFLAGS) AC_MSG_NOTICE([All Python build checks successful: 'yoda' Python extension will be built]) fi AM_CONDITIONAL(WITH_CYTHON, [test x$CYTHON_FOUND = xyes]) ## Extend and substitute the default build flags after lib testing AM_CPPFLAGS="-I\$(top_srcdir)/include -I\$(top_builddir)/include" AC_SUBST(AM_CPPFLAGS) AC_SUBST(AM_CXXFLAGS) dnl dnl setup.py puts its build artifacts into a labelled path dnl this helps the test scripts to find them locally instead of dnl having to install first dnl YODA_SETUP_PY_PATH=$(${PYTHON} -c 'from __future__ import print_function; import distutils.util as u, sys; vi=sys.version_info; print("lib.%s-%s.%s" % (u.get_platform(),vi.major, vi.minor))') AC_SUBST(YODA_SETUP_PY_PATH) ## Build Doxygen if possible AC_PATH_PROG(DOXYGEN, doxygen) AM_CONDITIONAL(WITH_DOXYGEN, test "$DOXYGEN") ## Build file output AC_EMPTY_SUBST AC_CONFIG_FILES([Makefile Doxyfile]) AC_CONFIG_FILES([include/Makefile include/YODA/Makefile]) AC_CONFIG_FILES([src/Makefile src/tinyxml/Makefile src/yamlcpp/Makefile ]) AC_CONFIG_FILES([tests/Makefile]) AC_CONFIG_FILES([pyext/Makefile pyext/setup.py pyext/yoda/Makefile ]) AC_CONFIG_FILES([bin/Makefile bin/yoda-config]) AC_CONFIG_FILES([yodaenv.sh yoda.pc]) AC_OUTPUT if test x$enable_pyext == xyes; then cat < namespace YODA { /// A 1D data point to be contained in a Scatter1D class Point1D : public Point { public: /// @name Constructors //@{ // Default constructor Point1D() { } /// Constructor from values with optional symmetric errors Point1D(double x, double ex=0.0, std::string source="") : _x(x) { _ex[source] = std::make_pair(ex, ex); } /// Constructor from values with explicit asymmetric errors Point1D(double x, double exminus, double explus, std::string source="") : _x(x) { _ex[source] = std::make_pair(exminus, explus); } /// Constructor from values with asymmetric errors Point1D(double x, const std::pair& ex, std::string source="") : _x(x) { _ex[source] = ex; } /// Copy constructor Point1D(const Point1D& p) : _x(p._x), _ex(p._ex) { this->setParentAO( p.getParentAO()); } /// Copy assignment Point1D& operator = (const Point1D& p) { _x = p._x; _ex = p._ex; this->setParentAO( p.getParentAO()); return *this; } //@} public: /// Space dimension of the point size_t dim() { return 1; } /// @name Value accessors //@{ /// Get x value double x() const { return _x; } /// Set x value void setX(double x) { _x = x; } /// @todo Uniform "coords" accessor across all Scatters: returning fixed-size tuple? //@} /// @name x error accessors //@{ /// Get x-error values const std::pair& xErrs( std::string source="") const { + if (source!="") getVariationsFromParent(); if (!_ex.count(source)) throw RangeError("xErrs has no such key: "+source); return _ex.at(source); } /// Get negative x-error value double xErrMinus( std::string source="") const { + if (source!="") getVariationsFromParent(); if (!_ex.count(source)) throw RangeError("xErrs has no such key: "+source); return _ex.at(source).first; } /// Get positive x-error value double xErrPlus( std::string source="") const { + if (source!="") getVariationsFromParent(); if (!_ex.count(source)) throw RangeError("xErrs has no such key: "+source); return _ex.at(source).second; } /// Get average x-error value double xErrAvg( std::string source="") const { + if (source!="") getVariationsFromParent(); if (!_ex.count(source)) throw RangeError("xErrs has no such key: "+source); return (_ex.at(source).first + _ex.at(source).second)/2.0; } /// Set negative x error void setXErrMinus(double exminus, std::string source="") { - if (source!="") getVariationsFromParent(); if (!_ex.count(source)) _ex[source] = std::make_pair(0.,0.); _ex.at(source).first = exminus; } /// Set positive x error void setXErrPlus(double explus, std::string source="") { - if (source!="") getVariationsFromParent(); if (!_ex.count(source)) _ex[source] = std::make_pair(0.,0.); _ex.at(source).second = explus; } /// Set symmetric x error void setXErr(double ex, std::string source="") { setXErrMinus(ex, source); setXErrPlus(ex, source); } /// Set symmetric x error (alias) void setXErrs(double ex, std::string source="") { setXErr(ex, source); } /// Set asymmetric x error void setXErrs(double exminus, double explus, std::string source="") { setXErrMinus(exminus, source); setXErrPlus(explus, source); } /// Set asymmetric x error void setXErrs(const std::pair& ex, std::string source="") { - if (source!="") getVariationsFromParent(); _ex[source] = ex; } /// Get value minus negative x-error double xMin(std::string source="") const { + if (source!="") getVariationsFromParent(); if (!_ex.count(source)) throw RangeError("xErrs has no such key: "+source); return _x - _ex.at(source).first; } /// Get value plus positive x-error double xMax(std::string source="") const { + if (source!="") getVariationsFromParent(); if (!_ex.count(source)) throw RangeError("xErrs has no such key: "+source); return _x + _ex.at(source).second; } //@} /// @name Combined x value and error setters //@{ /// Set x value and symmetric error void setX(double x, double ex, std::string source="") { setX(x); setXErr(ex, source); } /// Set x value and asymmetric error void setX(double x, double exminus, double explus, std::string source="") { setX(x); setXErrs(exminus, explus, source); } /// Set x value and asymmetric error void setX(double x, std::pair& ex, std::string source="") { setX(x); setXErrs(ex, source); } //@} // @name Manipulations //@{ /// Scaling of x axis void scaleX(double scalex) { setX(x()*scalex); for (const auto &source : _ex){ setXErrs(xErrMinus()*scalex, xErrPlus()*scalex, source.first); } } //@} /// @name Integer axis accessor equivalents //@{ /// Get the point value for direction @a i double val(size_t i) const { if (i == 0 || i > 1) throw RangeError("Invalid axis int, must be in range 1..dim"); return x(); } /// Set the point value for direction @a i void setVal(size_t i, double val) { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setX(val); } /// Get error map for direction @a i const std::map< std::string, std::pair> & errMap() const { getVariationsFromParent(); return _ex; } // Parse the variations from the parent AO if it exists void getVariationsFromParent() const; /// Get error values for direction @a i const std::pair& errs(size_t i, std::string source="") const { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); return xErrs(source); } /// Get negative error value for direction @a i double errMinus(size_t i, std::string source="") const { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); return xErrMinus(source); } /// Get positive error value for direction @a i double errPlus(size_t i, std::string source="") const { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); return xErrPlus(source); } /// Get average error value for direction @a i double errAvg(size_t i, std::string source="") const { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); return xErrAvg(source); } /// Set negative error for direction @a i void setErrMinus(size_t i, double eminus, std::string source="") { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setXErrMinus(eminus, source); } /// Set positive error for direction @a i void setErrPlus(size_t i, double eplus, std::string source="") { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setXErrPlus(eplus, source); } /// Set symmetric error for direction @a i void setErr(size_t i, double e, std::string source="") { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setXErr(e, source); } /// Set asymmetric error for direction @a i void setErrs(size_t i, double eminus, double eplus, std::string source="") { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setXErrs(eminus, eplus, source); } /// Set asymmetric error for direction @a i void setErrs(size_t i, std::pair& e, std::string source="") { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setXErrs(e, source); } /// Set value and symmetric error for direction @a i void set(size_t i, double val, double e, std::string source="") { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setX(val, e, source); } /// Set value and asymmetric error for direction @a i void set(size_t i, double val, double eminus, double eplus, std::string source="") { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setX(val, eminus, eplus, source); } /// Set value and asymmetric error for direction @a i void set(size_t i, double val, std::pair& e, std::string source="") { if (i != 1) throw RangeError("Invalid axis int, must be in range 1..dim"); setX(val, e, source); } //@} protected: /// @name Value and error variables //@{ double _x; // a map of the errors for each source. Nominal stored under "" // to ensure backward compatibility std::map< std::string, std::pair > _ex; //@} }; /// @name Comparison operators //@{ /// Equality test of x characteristics only /// @todo Base on a named fuzzyEquals(a,b,tol=1e-3) unbound function inline bool operator==(const YODA::Point1D& a, const YODA::Point1D& b) { if (!YODA::fuzzyEquals(a.x(), b.x()) || !YODA::fuzzyEquals(a.xErrMinus(), b.xErrMinus()) || !YODA::fuzzyEquals(a.xErrPlus(), b.xErrPlus()) ) return false; return true; } /// Equality test of x characteristics only inline bool operator != (const YODA::Point1D& a, const YODA::Point1D& b) { return !(a == b); } /// Less-than operator used to sort bins by x-ordering inline bool operator < (const YODA::Point1D& a, const YODA::Point1D& b) { if (!YODA::fuzzyEquals(a.x(), b.x())) { return a.x() < b.x(); } if (!YODA::fuzzyEquals(a.xErrMinus(), b.xErrMinus())) { return a.xErrMinus() < b.xErrMinus(); } if (!YODA::fuzzyEquals(a.xErrPlus(), b.xErrPlus())) { return a.xErrPlus() < b.xErrPlus(); } return false; } /// Less-than-or-equals operator used to sort bins by x-ordering inline bool operator <= (const YODA::Point1D& a, const YODA::Point1D& b) { if (a == b) return true; return a < b; } /// Greater-than operator used to sort bins by x-ordering inline bool operator > (const YODA::Point1D& a, const YODA::Point1D& b) { return !(a <= b); } /// Greater-than-or-equals operator used to sort bins by x-ordering inline bool operator >= (const YODA::Point1D& a, const YODA::Point1D& b) { return !(a < b); } //@} } #endif diff --git a/include/YODA/Point3D.h b/include/YODA/Point3D.h --- a/include/YODA/Point3D.h +++ b/include/YODA/Point3D.h @@ -1,695 +1,692 @@ // -*- C++ -*- // // This file is part of YODA -- Yet more Objects for Data Analysis // Copyright (C) 2008-2018 The YODA collaboration (see AUTHORS for details) // #ifndef YODA_POINT3D_H #define YODA_POINT3D_H #include "YODA/Point.h" #include "YODA/Exceptions.h" #include "YODA/Utils/MathUtils.h" #include namespace YODA { /// A 3D data point to be contained in a Scatter3D class Point3D : public Point { public: /// @name Constructors //@{ // Default constructor Point3D() { } /// Constructor from values with optional symmetric errors Point3D(double x, double y, double z, double ex=0.0, double ey=0.0, double ez=0.0, std::string source="") : _x(x), _y(y), _z(z) { _ex = std::make_pair(ex, ex); _ey = std::make_pair(ey, ey); _ez[source] = std::make_pair(ez, ez); } /// Constructor from values with explicit asymmetric errors Point3D(double x, double y, double z, double exminus, double explus, double eyminus, double eyplus, double ezminus, double ezplus, std::string source="") : _x(x), _y(y), _z(z) { _ex = std::make_pair(exminus, explus); _ey = std::make_pair(eyminus, eyplus); _ez[source] = std::make_pair(ezminus, ezplus); } /// Constructor from asymmetric errors given as vectors Point3D(double x, double y, double z, const std::pair& ex, const std::pair& ey, const std::pair& ez, std::string source="") : _x(x), _y(y), _z(z), _ex(ex), _ey(ey) { _ez[source] = ez; } /// Copy constructor Point3D(const Point3D& p) : _x(p._x), _y(p._y), _z(p._z), _ex(p._ex), _ey(p._ey), _ez(p._ez) { this->setParentAO( p.getParentAO()); } /// Copy assignment Point3D& operator = (const Point3D& p) { _x = p._x; _y = p._y; _z = p._z; _ex = p._ex; _ey = p._ey; _ez = p._ez; this->setParentAO( p.getParentAO()); return *this; } //@} public: /// Space dimension of the point size_t dim() { return 3; } /// @name Value and error accessors //@{ /// Get x value double x() const { return _x; } /// Set x value void setX(double x) { _x = x; } /// Get y value double y() const { return _y; } /// Set y value void setY(double y) { _y = y; } /// Get z value double z() const { return _z;} /// Set z value void setZ(double z) { _z = z;} /// @todo Uniform "coords" accessor across all Scatters: returning fixed-size tuple? // /// Get x,y,z value tuple // triple xyz() const { return std::triple(_x, _y, _z); } /// Set x, y and z values void setXYZ(double x, double y, double z) { setX(x); setY(y); setZ(z); } // /// Set x and y values // void setXY(triple xyz) { setX(xy.first); setY(xy.second); setZ(xy.third); } //@} /// @name x error accessors //@{ /// Get x-error values const std::pair& xErrs() const { return _ex; } /// Get negative x-error value double xErrMinus() const { return _ex.first; } /// Get positive x-error value double xErrPlus() const { return _ex.second; } /// Get average x-error value double xErrAvg() const { return (_ex.first + _ex.second)/2.0; } /// Set negative x error void setXErrMinus(double exminus) { _ex.first = exminus; } /// Set positive x error void setXErrPlus(double explus) { _ex.second = explus; } /// Set symmetric x error void setXErr(double ex) { setXErrMinus(ex); setXErrPlus(ex); } /// Set symmetric x error (alias) void setXErrs(double ex) { setXErr(ex); } /// Set asymmetric x error void setXErrs(double exminus, double explus) { setXErrMinus(exminus); setXErrPlus(explus); } /// Set asymmetric x error void setXErrs(const std::pair& ex) { _ex = ex; } /// Get value minus negative x-error double xMin() const { return _x - _ex.first; } /// Get value plus positive x-error double xMax() const { return _x + _ex.second; } //@} /// @name y error accessors //@{ /// Get y-error values const std::pair& yErrs() const { return _ey; } /// Get negative y-error value double yErrMinus() const { return _ey.first; } /// Get positive y-error value double yErrPlus() const { return _ey.second; } /// Get average y-error value double yErrAvg() const { return (_ey.first + _ey.second)/2.0; } /// Set negative y error void setYErrMinus(double eyminus) { _ey.first = eyminus; } /// Set positive y error void setYErrPlus(double eyplus) { _ey.second = eyplus; } /// Set symmetric y error void setYErr(double ey) { setYErrMinus(ey); setYErrPlus(ey); } /// Set symmetric y error (alias) void setYErrs(double ey) { setYErr(ey); } /// Set asymmetric y error void setYErrs(double eyminus, double eyplus) { setYErrMinus(eyminus); setYErrPlus(eyplus); } /// Set asymmetric y error void setYErrs(const std::pair& ey) { _ey = ey; } /// Get value minus negative y-error double yMin() const { return _y - _ey.first; } /// Get value plus positive y-error double yMax() const { return _y + _ey.second; } //@} /// @name z error accessors //@{ /// Get z-error values const std::pair& zErrs( std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ez.count(source)) throw RangeError("zErrs has no such key: "+source); return _ez.at(source); } /// Get negative z-error value double zErrMinus( std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ez.count(source)) throw RangeError("zErrs has no such key: "+source); return _ez.at(source).first; } /// Get positive z-error value double zErrPlus( std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ez.count(source)) throw RangeError("zErrs has no such key: "+source); return _ez.at(source).second; } /// Get average z-error value double zErrAvg( std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ez.count(source)) throw RangeError("zErrs has no such key: "+source); return (_ez.at(source).first + _ez.at(source).second)/2.0; } /// Set negative z error void setZErrMinus(double ezminus, std::string source="") { - if (source!="") getVariationsFromParent(); if (!_ez.count(source)) _ez[source] = std::make_pair(0.,0.); _ez.at(source).first = ezminus; } /// Set positive z error void setZErrPlus(double ezplus, std::string source="") { - if (source!="") getVariationsFromParent(); if (!_ez.count(source)) _ez[source] = std::make_pair(0.,0.); _ez.at(source).second = ezplus; } /// Set symmetric z error void setZErr(double ez, std::string source="") { setZErrMinus(ez, source); setZErrPlus(ez, source); } /// Set symmetric z error (alias) void setZErrs(double ez, std::string source="") { setZErr(ez, source); } /// Set asymmetric z error void setZErrs(double ezminus, double ezplus, std::string source="") { setZErrMinus(ezminus, source); setZErrPlus(ezplus, source); } /// Set asymmetric z error void setZErrs(const std::pair& ez, std::string source="") { - if (source!="") getVariationsFromParent(); _ez[source] = ez; } /// Get value minus negative z-error double zMin( std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ez.count(source)) throw RangeError("zErrs has no such key: "+source); return _z - _ez.at(source).first; } /// Get value plus positive z-error double zMax( std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ez.count(source)) throw RangeError("zErrs has no such key: "+source); return _z + _ez.at(source).second; } //@} /// @name Combined x/y value and error setters //@{ /// Set x value and symmetric error void setX(double x, double ex) { setX(x); setXErrs(ex); } /// Set x value and asymmetric error void setX(double x, double exminus, double explus) { setX(x); setXErrs(exminus, explus); } /// Set x value and asymmetric error void setX(double x, std::pair& ex) { setX(x); setXErrs(ex); } /// Set y value and symmetric error void setY(double y, double ey) { setY(y); setYErrs(ey); } /// Set y value and asymmetric error void setY(double y, double eyminus, double eyplus) { setY(y); setYErrs(eyminus, eyplus); } /// Set y value and asymmetric error void setY(double y, std::pair& ey) { setY(y); setYErrs(ey); } /// Set z value and symmetric error void setZ(double z, double ez, std::string source="") { setZ(z); setZErrs(ez, source); } /// Set z value and asymmetric error void setZ(double z, double ezminus, double ezplus, std::string source="") { setZ(z); setZErrs(ezminus, ezplus, source); } /// Set z value and asymmetric error void setZ(double z, std::pair& ez, std::string source="") { setZ(z); setZErrs(ez, source); } //@} // @name Manipulations //@{ /// Scaling of x axis void scaleX(double scalex) { setX(x()*scalex); setXErrs(xErrMinus()*scalex, xErrPlus()*scalex); } /// Scaling of y axis void scaleY(double scaley) { setY(y()*scaley); setYErrs(yErrMinus()*scaley, yErrPlus()*scaley); } /// Scaling of z axis void scaleZ(double scalez) { setZ(z()*scalez); for (const auto &source : _ez){ setZErrs(zErrMinus()*scalez, zErrPlus()*scalez, source.first); } } /// Scaling of all three axes void scaleXYZ(double scalex, double scaley, double scalez) { scaleX(scalex); scaleY(scaley); scaleZ(scalez); } /// Scaling of both axes /// @deprecated Use scaleXYZ void scale(double scalex, double scaley, double scalez) { scaleXYZ(scalex, scaley, scalez); } //@} /// @name Integer axis accessor equivalents //@{ /// Get the point value for direction @a i double val(size_t i) const { switch (i) { case 1: return x(); case 2: return y(); case 3: return z(); default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set the point value for direction @a i void setVal(size_t i, double val) { switch (i) { case 1: setX(val); break; case 2: setY(val); break; case 3: setZ(val); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Get error map for direction @a i const std::map< std::string, std::pair> & errMap() const { getVariationsFromParent(); return _ez; } // Parse the variations from the parent AO if it exists void getVariationsFromParent() const; /// Get error values for direction @a i const std::pair& errs(size_t i, std::string source="") const { switch (i) { case 1: return xErrs(); case 2: return yErrs(); case 3: return zErrs(source); default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Get negative error value for direction @a i double errMinus(size_t i, std::string source="") const { switch (i) { case 1: return xErrMinus(); case 2: return yErrMinus(); case 3: return zErrMinus(source); default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Get positive error value for direction @a i double errPlus(size_t i, std::string source="") const { switch (i) { case 1: return xErrPlus(); case 2: return yErrPlus(); case 3: return zErrPlus(source); default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Get average error value for direction @a i double errAvg(size_t i, std::string source="") const { switch (i) { case 1: return xErrAvg(); case 2: return yErrAvg(); case 3: return zErrAvg(source); default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set negative error for direction @a i void setErrMinus(size_t i, double eminus, std::string source="") { switch (i) { case 1: setXErrMinus(eminus); break; case 2: setYErrMinus(eminus); break; case 3: setZErrMinus(eminus, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set positive error for direction @a i void setErrPlus(size_t i, double eplus, std::string source="") { switch (i) { case 1: setXErrPlus(eplus); break; case 2: setYErrPlus(eplus); break; case 3: setZErrPlus(eplus, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set symmetric error for direction @a i void setErr(size_t i, double e, std::string source="") { switch (i) { case 1: setXErrs(e); break; case 2: setYErrs(e); break; case 3: setZErrs(e, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set asymmetric error for direction @a i void setErrs(size_t i, double eminus, double eplus, std::string source="") { switch (i) { case 1: setXErrs(eminus, eplus); break; case 2: setYErrs(eminus, eplus); break; case 3: setZErrs(eminus, eplus, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set asymmetric error for direction @a i void setErrs(size_t i, std::pair& e, std::string source="") { switch (i) { case 1: setXErrs(e); break; case 2: setYErrs(e); break; case 3: setZErrs(e, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set value and symmetric error for direction @a i void set(size_t i, double val, double e, std::string source="") { switch (i) { case 1: setX(val, e); break; case 2: setY(val, e); break; case 3: setZ(val, e, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set value and asymmetric error for direction @a i void set(size_t i, double val, double eminus, double eplus, std::string source="") { switch (i) { case 1: setX(val, eminus, eplus); break; case 2: setY(val, eminus, eplus); break; case 3: setZ(val, eminus, eplus, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } /// Set value and asymmetric error for direction @a i void set(size_t i, double val, std::pair& e, std::string source="") { switch (i) { case 1: setX(val, e); break; case 2: setY(val, e); break; case 3: setZ(val, e, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } //@} protected: /// @name Value and error variables //@{ double _x; double _y; double _z; std::pair _ex; std::pair _ey; // a map of the errors for each source. Nominal stored under "" // to ensure backward compatibility std::map< std::string, std::pair >_ez; //@} }; /// @name Comparison operators //@{ /// Equality test of x, y & z characteristics only /// @todo Base on a named fuzzyEquals(a,b,tol=1e-3) unbound function inline bool operator==(const Point3D& a, const YODA::Point3D& b) { if (!YODA::fuzzyEquals(a.x(), b.x()) || !YODA::fuzzyEquals(a.xErrMinus(), b.xErrMinus()) || !YODA::fuzzyEquals(a.xErrPlus(), b.xErrPlus()) ) return false; if (!YODA::fuzzyEquals(a.y(), b.y()) || !YODA::fuzzyEquals(a.yErrMinus(), b.yErrMinus()) || !YODA::fuzzyEquals(a.yErrPlus(), b.yErrPlus()) ) return false; if (!YODA::fuzzyEquals(a.z(), b.z()) || !YODA::fuzzyEquals(a.zErrMinus(), b.zErrMinus()) || !YODA::fuzzyEquals(a.zErrPlus(), b.zErrPlus()) ) return false; return true; const bool same_val = fuzzyEquals(a.x(), b.x()) && fuzzyEquals(a.y(), b.y()); const bool same_eminus = fuzzyEquals(a.xErrMinus(), b.xErrMinus()) && fuzzyEquals(a.yErrMinus(), b.yErrMinus()); const bool same_eplus = fuzzyEquals(a.xErrPlus(), b.xErrPlus()) && fuzzyEquals(a.yErrPlus(), b.yErrPlus()); return same_val && same_eminus && same_eplus; } /// Inequality operator inline bool operator != (const Point3D& a, const YODA::Point3D& b) { return !(a == b); } /// Less-than operator used to sort bins by x-first ordering inline bool operator < (const Point3D& a, const YODA::Point3D& b) { if (! fuzzyEquals(a.x(), b.x())) { return a.x() < b.x(); } if (!fuzzyEquals(a.y(), b.y())) { return a.y() < b.y(); } if (! fuzzyEquals(a.xErrMinus(), b.xErrMinus())) { return a.xErrMinus() < b.xErrMinus(); } if (!fuzzyEquals(a.yErrMinus(), b.yErrMinus())) { return a.yErrMinus() < b.yErrMinus(); } if (! fuzzyEquals(a.xErrPlus(), b.xErrPlus())) { return a.xErrPlus() < b.xErrPlus(); } if (!fuzzyEquals(a.yErrPlus(), b.yErrPlus())) { return a.yErrPlus() < b.yErrPlus(); } return false; } /// Less-than-or-equals operator inline bool operator <= (const Point3D& a, const YODA::Point3D& b) { if (a == b) return true; return a < b; } /// Greater-than operator inline bool operator > (const Point3D& a, const YODA::Point3D& b) { return !(a <= b); } /// Greater-than-or-equals operator inline bool operator >= (const Point3D& a, const YODA::Point3D& b) { return !(a < b); } //@} } #endif diff --git a/src/Scatter1D.cc b/src/Scatter1D.cc --- a/src/Scatter1D.cc +++ b/src/Scatter1D.cc @@ -1,57 +1,59 @@ #include "YODA/Scatter1D.h" #include "YODA/Counter.h" #include #include "yaml-cpp/yaml.h" #ifdef YAML_NAMESPACE #define YAML YAML_NAMESPACE #endif namespace YODA { /// Make a Scatter1D representation of a Histo1D Scatter1D mkScatter(const Counter& c) { Scatter1D rtn; for (const std::string& a : c.annotations()) rtn.setAnnotation(a, c.annotation(a)); rtn.setAnnotation("Type", c.type()); // might override the copied ones - rtn.addPoint(c.val(), c.err()); + Point1D pt(c.val(), c.err()); + pt.setParentAO(&rtn); + rtn.addPoint(pt); return rtn; } void Scatter1D::parseVariations() { if (this-> _variationsParsed) { return;} if (!(this->hasAnnotation("ErrorBreakdown"))) { return; } YAML::Node errorBreakdown; errorBreakdown = YAML::Load(this->annotation("ErrorBreakdown")); if (errorBreakdown.size()) { for (unsigned int thisPointIndex=0 ; thisPointIndex< this->numPoints() ; ++thisPointIndex){ Point1D &thispoint = this->_points[thisPointIndex]; YAML::Node variations = errorBreakdown[thisPointIndex]; for (const auto& variation : variations) { const std::string variationName = variation.first.as(); double eyp = variation.second["up"].as(); double eym = variation.second["dn"].as(); thispoint.setXErrs(eym,eyp,variationName); } } this-> _variationsParsed =true; } } const std::vector Scatter1D::variations() const { std::vector vecvariations; for (auto &point : this->_points){ for (auto &it : point.errMap()){ //if the variation is not already in the vector, add it ! if (std::find(vecvariations.begin(), vecvariations.end(), it.first) == vecvariations.end()){ vecvariations.push_back(it.first); } } } return vecvariations; } } diff --git a/src/Scatter3D.cc b/src/Scatter3D.cc --- a/src/Scatter3D.cc +++ b/src/Scatter3D.cc @@ -1,154 +1,156 @@ #include "YODA/Scatter3D.h" #include "YODA/Histo2D.h" #include "YODA/Profile2D.h" #include "YODA/Exceptions.h" #include #include "yaml-cpp/yaml.h" #ifdef YAML_NAMESPACE #define YAML YAML_NAMESPACE #endif namespace YODA { Scatter3D mkScatter(const Histo2D& h, bool usefocus, bool binareadiv) { Scatter3D rtn; for (const std::string& a : h.annotations()) rtn.setAnnotation(a, h.annotation(a)); rtn.setAnnotation("Type", h.type()); for (size_t i = 0; i < h.numBins(); ++i) { const HistoBin2D& b = h.bin(i); /// SAME FOR ALL 2D BINS double x = b.xMid(); if (usefocus) { try { x = b.xFocus(); } catch (const LowStatsError& lse) { x = b.xMid(); } } const double exminus = x - b.xMin(); const double explus = b.xMax() - x; double y = b.yMid(); if (usefocus) { try { y = b.yFocus(); } catch (const LowStatsError& lse) { y = b.yMid(); } } const double eyminus = y - b.yMin(); const double eyplus = b.yMax() - y; /// END SAME FOR ALL 2D BINS double z; try { z = b.sumW(); } catch (const Exception&) { // LowStatsError or WeightError z = std::numeric_limits::quiet_NaN(); } if (binareadiv) z /= b.xWidth()*b.yWidth(); const double ez = b.relErr() * z; - /// @todo Set up parent link cf. Scatter2D - rtn.addPoint(x, y, z, exminus, explus, eyminus, eyplus, ez, ez); + + Point3D pt(x, y, z, exminus, explus, eyminus, eyplus, ez, ez); + pt.setParentAO(&rtn); + rtn.addPoint(pt); } assert(h.numBins() == rtn.numPoints()); return rtn; } Scatter3D mkScatter(const Profile2D& h, bool usefocus, bool usestddev) { Scatter3D rtn; for (const std::string& a : h.annotations()) rtn.setAnnotation(a, h.annotation(a)); rtn.setAnnotation("Type", h.type()); for (size_t i = 0; i < h.numBins(); ++i) { const ProfileBin2D& b = h.bin(i); /// SAME FOR ALL 2D BINS double x = b.xMid(); if (usefocus) { try { x = b.xFocus(); } catch (const LowStatsError& lse) { x = b.xMid(); } } const double exminus = x - b.xMin(); const double explus = b.xMax() - x; double y = b.yMid(); if (usefocus) { try { y = b.yFocus(); } catch (const LowStatsError& lse) { y = b.yMid(); } } const double eyminus = y - b.yMin(); const double eyplus = b.yMax() - y; /// END SAME FOR ALL 2D BINS double z; try { z = b.mean(); } catch (const LowStatsError& lse) { z = std::numeric_limits::quiet_NaN(); } double ez; try { ez = usestddev ? b.stdDev() : b.stdErr(); ///< Control z-error scheme via usestddev arg } catch (const LowStatsError& lse) { ez = std::numeric_limits::quiet_NaN(); } rtn.addPoint(x, y, z, exminus, explus, eyminus, eyplus, ez, ez); } return rtn; } void Scatter3D::parseVariations() { if (this-> _variationsParsed) { return; } if (!(this->hasAnnotation("ErrorBreakdown"))) { return;} YAML::Node errorBreakdown; errorBreakdown = YAML::Load(this->annotation("ErrorBreakdown")); if (errorBreakdown.size()) { for (unsigned int thisPointIndex=0 ; thisPointIndex< this->numPoints() ; ++thisPointIndex){ Point3D &thispoint = this->_points[thisPointIndex]; YAML::Node variations = errorBreakdown[thisPointIndex]; for (const auto& variation : variations) { const std::string variationName = variation.first.as(); double eyp = variation.second["up"].as(); double eym = variation.second["dn"].as(); thispoint.setZErrs(eym,eyp,variationName); } } this-> _variationsParsed =true; } } const std::vector Scatter3D::variations() const { std::vector vecvariations; for (auto &point : this->_points){ for (auto &it : point.errMap()){ //if the variation is not already in the vector, add it ! if (std::find(vecvariations.begin(), vecvariations.end(), it.first) == vecvariations.end()){ vecvariations.push_back(it.first); } } } return vecvariations; } } diff --git a/src/WriterYODA.cc b/src/WriterYODA.cc --- a/src/WriterYODA.cc +++ b/src/WriterYODA.cc @@ -1,402 +1,328 @@ // -*- C++ -*- // // This file is part of YODA -- Yet more Objects for Data Analysis // Copyright (C) 2008-2018 The YODA collaboration (see AUTHORS for details) // #include "YODA/WriterYODA.h" #include "yaml-cpp/yaml.h" #ifdef YAML_NAMESPACE #define YAML YAML_NAMESPACE #endif #include #include using namespace std; namespace YODA { /// Singleton creation function Writer& WriterYODA::create() { static WriterYODA _instance; _instance.setPrecision(6); return _instance; } // Format version: // - V1/empty = make-plots annotations style // - V2 = YAML annotations static const int YODA_FORMAT_VERSION = 2; // Version-formatting helper function inline string _iotypestr(const string& baseiotype) { ostringstream os; os << "YODA_" << Utils::toUpper(baseiotype) << "_V" << YODA_FORMAT_VERSION; return os.str(); } void WriterYODA::_writeAnnotations(std::ostream& os, const AnalysisObject& ao) { os << scientific << setprecision(_precision); for (const string& a : ao.annotations()) { if (a.empty()) continue; /// @todo Write out floating point annotations as scientific notation string ann = ao.annotation(a); // remove stpurious line returns at the end of a string so that we don't // end up with two line returns. ann.erase(std::remove(ann.begin(), ann.end(), '\n'), ann.end()); os << a << ": " << ann << "\n"; } os << "---\n"; } void WriterYODA::writeCounter(std::ostream& os, const Counter& c) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("COUNTER") << " " << c.path() << "\n"; _writeAnnotations(os, c); os << "# sumW\t sumW2\t numEntries\n"; os << c.sumW() << "\t" << c.sumW2() << "\t" << c.numEntries() << "\n"; os << "END " << _iotypestr("COUNTER") << "\n\n"; os.flags(oldflags); } void WriterYODA::writeHisto1D(std::ostream& os, const Histo1D& h) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("HISTO1D") << " " << h.path() << "\n"; _writeAnnotations(os, h); try { //if ( h.totalDbn().effNumEntries() > 0 ) { os << "# Mean: " << h.xMean() << "\n"; os << "# Area: " << h.integral() << "\n"; } catch (LowStatsError& e) { // } os << "# ID\t ID\t sumw\t sumw2\t sumwx\t sumwx2\t numEntries\n"; os << "Total \tTotal \t"; os << h.totalDbn().sumW() << "\t" << h.totalDbn().sumW2() << "\t"; os << h.totalDbn().sumWX() << "\t" << h.totalDbn().sumWX2() << "\t"; os << h.totalDbn().numEntries() << "\n"; os << "Underflow\tUnderflow\t"; os << h.underflow().sumW() << "\t" << h.underflow().sumW2() << "\t"; os << h.underflow().sumWX() << "\t" << h.underflow().sumWX2() << "\t"; os << h.underflow().numEntries() << "\n"; os << "Overflow\tOverflow\t"; os << h.overflow().sumW() << "\t" << h.overflow().sumW2() << "\t"; os << h.overflow().sumWX() << "\t" << h.overflow().sumWX2() << "\t"; os << h.overflow().numEntries() << "\n"; os << "# xlow\t xhigh\t sumw\t sumw2\t sumwx\t sumwx2\t numEntries\n"; for (const HistoBin1D& b : h.bins()) { os << b.xMin() << "\t" << b.xMax() << "\t"; os << b.sumW() << "\t" << b.sumW2() << "\t"; os << b.sumWX() << "\t" << b.sumWX2() << "\t"; os << b.numEntries() << "\n"; } os << "END " << _iotypestr("HISTO1D") << "\n\n"; os.flags(oldflags); } void WriterYODA::writeHisto2D(std::ostream& os, const Histo2D& h) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("HISTO2D") << " " << h.path() << "\n"; _writeAnnotations(os, h); try { //if ( h.totalDbn().numEntries() > 0 ) os << "# Mean: (" << h.xMean() << ", " << h.yMean() << ")\n"; os << "# Volume: " << h.integral() << "\n"; } catch (LowStatsError& e) { // } os << "# ID\t ID\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t sumwxy\t numEntries\n"; // Total distribution const Dbn2D& td = h.totalDbn(); os << "Total \tTotal \t"; os << td.sumW() << "\t" << td.sumW2() << "\t"; os << td.sumWX() << "\t" << td.sumWX2() << "\t"; os << td.sumWY() << "\t" << td.sumWY2() << "\t"; os << td.sumWXY() << "\t"; os << td.numEntries() << "\n"; // Outflows /// @todo Disabled for now, reinstate with a *full* set of outflow info to allow marginalisation os << "# 2D outflow persistency not currently supported until API is stable\n"; // for (int ix = -1; ix <= 1; ++ix) { // for (int iy = -1; iy <= 1; ++iy) { // if (ix == 0 && iy == 0) continue; // os << "Outflow\t" << ix << ":" << iy << "\t"; // const Dbn2D& d = h.outflow(ix, iy); // os << d.sumW() << "\t" << d.sumW2() << "\t"; // os << d.sumWX() << "\t" << d.sumWX2() << "\t"; // os << d.sumWY() << "\t" << d.sumWY2() << "\t"; // os << d.sumWXY() << "\t"; // os << d.numEntries() << "\n"; // } // } // Bins os << "# xlow\t xhigh\t ylow\t yhigh\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t sumwxy\t numEntries\n"; for (const HistoBin2D& b : h.bins()) { os << b.xMin() << "\t" << b.xMax() << "\t"; os << b.yMin() << "\t" << b.yMax() << "\t"; os << b.sumW() << "\t" << b.sumW2() << "\t"; os << b.sumWX() << "\t" << b.sumWX2() << "\t"; os << b.sumWY() << "\t" << b.sumWY2() << "\t"; os << b.sumWXY() << "\t"; os << b.numEntries() << "\n"; } os << "END " << _iotypestr("HISTO2D") << "\n\n"; os.flags(oldflags); } void WriterYODA::writeProfile1D(std::ostream& os, const Profile1D& p) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("PROFILE1D") << " " << p.path() << "\n"; _writeAnnotations(os, p); os << "# ID\t ID\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t numEntries\n"; os << "Total \tTotal \t"; os << p.totalDbn().sumW() << "\t" << p.totalDbn().sumW2() << "\t"; os << p.totalDbn().sumWX() << "\t" << p.totalDbn().sumWX2() << "\t"; os << p.totalDbn().sumWY() << "\t" << p.totalDbn().sumWY2() << "\t"; os << p.totalDbn().numEntries() << "\n"; os << "Underflow\tUnderflow\t"; os << p.underflow().sumW() << "\t" << p.underflow().sumW2() << "\t"; os << p.underflow().sumWX() << "\t" << p.underflow().sumWX2() << "\t"; os << p.underflow().sumWY() << "\t" << p.underflow().sumWY2() << "\t"; os << p.underflow().numEntries() << "\n"; os << "Overflow\tOverflow\t"; os << p.overflow().sumW() << "\t" << p.overflow().sumW2() << "\t"; os << p.overflow().sumWX() << "\t" << p.overflow().sumWX2() << "\t"; os << p.overflow().sumWY() << "\t" << p.overflow().sumWY2() << "\t"; os << p.overflow().numEntries() << "\n"; os << "# xlow\t xhigh\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t numEntries\n"; for (const ProfileBin1D& b : p.bins()) { os << b.xMin() << "\t" << b.xMax() << "\t"; os << b.sumW() << "\t" << b.sumW2() << "\t"; os << b.sumWX() << "\t" << b.sumWX2() << "\t"; os << b.sumWY() << "\t" << b.sumWY2() << "\t"; os << b.numEntries() << "\n"; } os << "END " << _iotypestr("PROFILE1D") << "\n\n"; os.flags(oldflags); } void WriterYODA::writeProfile2D(std::ostream& os, const Profile2D& p) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("PROFILE2D") << " " << p.path() << "\n"; _writeAnnotations(os, p); os << "# sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t sumwz\t sumwz2\t sumwxy\t numEntries\n"; // Total distribution const Dbn3D& td = p.totalDbn(); os << "Total \tTotal \t"; os << td.sumW() << "\t" << td.sumW2() << "\t"; os << td.sumWX() << "\t" << td.sumWX2() << "\t"; os << td.sumWY() << "\t" << td.sumWY2() << "\t"; os << td.sumWZ() << "\t" << td.sumWZ2() << "\t"; os << td.sumWXY() << "\t"; // << td.sumWXZ() << "\t" << td.sumWYZ() << "\t"; os << td.numEntries() << "\n"; // Outflows /// @todo Disabled for now, reinstate with a *full* set of outflow info to allow marginalisation os << "# 2D outflow persistency not currently supported until API is stable\n"; // for (int ix = -1; ix <= 1; ++ix) { // for (int iy = -1; iy <= 1; ++iy) { // if (ix == 0 && iy == 0) continue; // os << "Outflow\t" << ix << ":" << iy << "\t"; // const Dbn3D& d = p.outflow(ix, iy); // os << d.sumW() << "\t" << d.sumW2() << "\t"; // os << d.sumWX() << "\t" << d.sumWX2() << "\t"; // os << d.sumWY() << "\t" << d.sumWY2() << "\t"; // os << d.sumWZ() << "\t" << d.sumWZ2() << "\t"; // os << d.sumWXY() << "\t"; // << d.sumWXZ() << "\t" << d.sumWYZ() << "\t"; // os << d.numEntries() << "\n"; // } // } // Bins os << "# xlow\t xhigh\t ylow\t yhigh\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t sumwz\t sumwz2\t sumwxy\t numEntries\n"; for (const ProfileBin2D& b : p.bins()) { os << b.xMin() << "\t" << b.xMax() << "\t"; os << b.yMin() << "\t" << b.yMax() << "\t"; os << b.sumW() << "\t" << b.sumW2() << "\t"; os << b.sumWX() << "\t" << b.sumWX2() << "\t"; os << b.sumWY() << "\t" << b.sumWY2() << "\t"; os << b.sumWZ() << "\t" << b.sumWZ2() << "\t"; os << b.sumWXY() << "\t"; // << b.sumWXZ() << "\t" << b.sumWYZ() << "\t"; os << b.numEntries() << "\n"; } os << "END " << _iotypestr("PROFILE2D") << "\n\n"; os.flags(oldflags); } void WriterYODA::writeScatter1D(std::ostream& os, const Scatter1D& s) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("SCATTER1D") << " " << s.path() << "\n"; - //first write the Variations, a dummy annotation which - //contains the additional columns which will be written out - //for sytematic variations - YAML::Emitter out; - out << YAML::Flow ; - out << s.variations(); - //os << "Variations" << ": " << out.c_str() << "\n"; - // then write the regular annotations _writeAnnotations(os, s); - std::vector variations= s.variations(); - //write headers - std::string headers="# xval\t "; - for (const auto &source : variations){ - headers+=" xerr-"+source+"\t xerr+"+source+"\t"; - } + std::string headers="# xval\t xerr-\t xerr+\t"; os << headers << "\n"; //write points for (const Point1D& pt : s.points()) { // fill central value - os << pt.x(); - // fill errors for variations. The first should always be "" which is nominal. - // Assumes here that all points in the Scatter have the same - // variations... if not a range error will get thrown from - // the point when the user tries to access a variation it - // doesn't have... @todo maybe better way to do this? - for (const auto &source : variations){ - os << "\t" << pt.xErrMinus(source) << "\t" << pt.xErrPlus(source) ; - } + os << pt.x() << "\t" << pt.xErrMinus() << "\t" << pt.xErrPlus() ; os << "\n"; } os << "END " << _iotypestr("SCATTER1D") << "\n\n"; os << flush; os.flags(oldflags); } void WriterYODA::writeScatter2D(std::ostream& os, const Scatter2D& s) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); - os << "BEGIN " << _iotypestr("SCATTER2D") << " " << s.path() << "\n"; - //first write the Variations, a dummy annotation which - //contains the additional columns which will be written out - //for sytematic variations - YAML::Emitter out; - out << YAML::Flow << YAML::BeginMap; - int counter=0; - std::vector variations= s.variations(); - //write ErrBreakdown Annotation - for (const Point2D& pt : s.points()) { - out << YAML::Key << counter; - out << YAML::Value << YAML::BeginMap; - for (const auto &source : variations){ - if (source.length()==0) continue; - out << YAML::Key << source; - out << YAML::Value << YAML::BeginMap; - out << YAML::Key << "dn" << YAML::Value << pt.yErrMinus(source); - out << YAML::Key << "up" << YAML::Value << pt.yErrPlus(source); - out << YAML::EndMap; - } - out << YAML::EndMap; - } - out << YAML::EndMap; - os << "ErrorBreakdown" << ": " << out.c_str() << "\n"; - // then write the regular annotations + // write annotations _writeAnnotations(os, s); //write headers /// @todo Change ordering to {vals} {errs} {errs} ... std::string headers="# xval\t xerr-\t xerr+\t yval\t yerr-\t yerr+\t"; - //for (const auto &source : variations){ - // headers+=" yerr-"+source+"\t yerr+"+source+"\t"; - //} os << headers << "\n"; //write points for (const Point2D& pt : s.points()) { /// @todo Change ordering to {vals} {errs} {errs} ... // fill central value os << pt.x() << "\t" << pt.xErrMinus() << "\t" << pt.xErrPlus() << "\t"; - os << pt.y(); - // fill errors for variations. The first should always be "" which is nominal. - // Assumes here that all points in the Scatter have the same - // variations... if not a range error will get thrown from - // the point when the user tries to access a variation it - // doesn't have... @todo maybe better way to do this? - //for (const auto &source : variations){ - os << "\t" << pt.yErrMinus() << "\t" << pt.yErrPlus() ; - // } + os << pt.y() << "\t" << pt.yErrMinus() << "\t" << pt.yErrPlus() ; os << "\n"; } os << "END " << _iotypestr("SCATTER2D") << "\n\n"; os << flush; os.flags(oldflags); } void WriterYODA::writeScatter3D(std::ostream& os, const Scatter3D& s) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); - os << "BEGIN " << _iotypestr("SCATTER3D") << " " << s.path() << "\n"; - //first write the Variations, a dummy annotation which - //contains the additional columns which will be written out - //for sytematic variations - YAML::Emitter out; - out << YAML::Flow ; - out << s.variations(); - // then write the regular annotations + // write the annotations _writeAnnotations(os, s); - std::vector variations= s.variations(); + // std::vector variations= s.variations(); //write headers /// @todo Change ordering to {vals} {errs} {errs} ... std::string headers="# xval\t xerr-\t xerr+\t yval\t yerr-\t yerr+\t zval\t "; - for (const auto &source : variations){ - headers+=" zerr-"+source+"\t zerr+"+source+"\t"; - } os << headers << "\n"; //write points for (const Point3D& pt : s.points()) { /// @todo Change ordering to {vals} {errs} {errs} ... // fill central value os << pt.x() << "\t" << pt.xErrMinus() << "\t" << pt.xErrPlus() << "\t"; os << pt.y() << "\t" << pt.yErrMinus() << "\t" << pt.yErrPlus() << "\t"; - os << pt.z(); - // fill errors for variations. The first should always be "" which is nominal. - // Assumes here that all points in the Scatter have the same - // variations... if not a range error will get thrown from - // the point when the user tries to access a variation it - // doesn't have... @todo maybe better way to do this? - for (const auto &source : variations){ - os << "\t" << pt.zErrMinus(source) << "\t" << pt.zErrPlus(source) ; - } + os << pt.z() << "\t" << pt.zErrMinus() << "\t" << pt.zErrPlus() ; os << "\n"; } os << "END " << _iotypestr("SCATTER3D") << "\n\n"; os << flush; os.flags(oldflags); } }