No OneTemporary
Actions

Size

61 KB

Subscribers

None

View Options

	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 <<EOF

	************************************************************
	YODA CONFIGURED!

	Now build and install (to the $prefix tree) with e.g.
	make -j2 && make -j2 install
	************************************************************
	EOF
	fi
	diff --git a/include/YODA/Point1D.h b/include/YODA/Point1D.h
	--- a/include/YODA/Point1D.h
	+++ b/include/YODA/Point1D.h
	@@ -1,366 +1,369 @@
	// -- 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_POINT1D_H
	#define YODA_POINT1D_H

	#include "YODA/Point.h"
	#include "YODA/Exceptions.h"
	#include "YODA/Utils/MathUtils.h"
	#include <utility>

	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<double,double>& 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<double,double>& 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<double,double>& 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<double,double>& 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<double,double>> & 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<double,double>& 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<double,double>& 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<double,double>& 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<double,double> > _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 <utility>

	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<double,double>& ex,
	const std::pair<double,double>& ey,
	const std::pair<double,double>& 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<double,double,double> 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<double,double,double> xyz) { setX(xy.first); setY(xy.second); setZ(xy.third); }

	//@}


	/// @name x error accessors
	//@{

	/// Get x-error values
	const std::pair<double,double>& 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<double,double>& 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<double,double>& 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<double,double>& 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<double,double>& 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<double,double>& 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<double,double>& 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<double,double>& 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<double,double>& 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<double,double>> & 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<double,double>& 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<double,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 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<double,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");
	}
	}

	//@}


	protected:

	/// @name Value and error variables
	//@{

	double _x;
	double _y;
	double _z;
	std::pair<double,double> _ex;
	std::pair<double,double> _ey;
	// a map of the errors for each source. Nominal stored under ""
	// to ensure backward compatibility
	std::map< std::string, std::pair<double,double> >_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 <sstream>
	#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<std::string>();
	double eyp = variation.second["up"].as<double>();
	double eym = variation.second["dn"].as<double>();
	thispoint.setXErrs(eym,eyp,variationName);
	}
	}
	this-> _variationsParsed =true;
	}
	}

	const std::vector<std::string> Scatter1D::variations() const {
	std::vector<std::string> 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 <sstream>
	#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<double>::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<double>::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<double>::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<std::string>();
	double eyp = variation.second["up"].as<double>();
	double eym = variation.second["dn"].as<double>();
	thispoint.setZErrs(eym,eyp,variationName);
	}
	}
	this-> _variationsParsed =true;
	}
	}


	const std::vector<std::string> Scatter3D::variations() const {
	std::vector<std::string> 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 <iostream>
	#include <iomanip>
	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<std::string> 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<std::string> 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<std::string> variations= s.variations();
	+ // std::vector<std::string> 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);
	}


	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 3:12 PM (37 m, 41 s)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486871
Default Alt Text: (61 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions