diff --git a/include/YODA/AnalysisObject.h b/include/YODA/AnalysisObject.h --- a/include/YODA/AnalysisObject.h +++ b/include/YODA/AnalysisObject.h @@ -1,308 +1,304 @@ // -*- 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_AnalysisObject_h #define YODA_AnalysisObject_h #include "YODA/Exceptions.h" #include "YODA/Utils/StringUtils.h" #include "YODA/Config/BuildConfig.h" #include #include #include -#include namespace YODA { /// AnalysisObject is the base class for histograms and scatters class AnalysisObject { public: /// Collection type for annotations, as a string-string map. typedef std::map Annotations; /// @name Creation and destruction //@{ /// Default constructor AnalysisObject() { } /// Constructor giving a type, a path and an optional title AnalysisObject(const std::string& type, const std::string& path, const std::string& title="") { setAnnotation("Type", type); setPath(path); setTitle(title); } /// Constructor giving a type, a path, another AO to copy annotation from, and an optional title AnalysisObject(const std::string& type, const std::string& path, const AnalysisObject& ao, const std::string& title="") { for (const std::string& a : ao.annotations()) setAnnotation(a, ao.annotation(a)); setAnnotation("Type", type); // might override the copied ones setPath(path); setTitle(title); } // /// Default copy constructor // AnalysisObject(const AnalysisObject& ao) { // if (ao.path().length() > 0) setPath(ao.path()); // if (ao.title().length() > 0) setTitle(ao.title()); // } /// Default destructor virtual ~AnalysisObject() { } /// Default copy assignment operator virtual AnalysisObject& operator = (const AnalysisObject& ao) { if (ao.path().length() > 0) setPath(ao.path()); if (ao.title().length() > 0) setTitle(ao.title()); return *this; } /// Make a copy on the heap, via 'new' virtual AnalysisObject* newclone() const = 0; //@} /// @name Modifiers //@{ /// Reset this analysis object virtual void reset() = 0; // variation parser - void parseVariations(){ - std::cout<< "LC DEBUG - using the analysis object parseVariations" << std:: endl; - return ; - } + void parseVariations(){ return ; } //@} ///@name Annotations //@{ /// Get all the annotation names /// @todo Change this to return the str->str map, with a separate annotationKeys, etc. std::vector annotations() const { std::vector rtn; rtn.reserve(_annotations.size()); for (const Annotations::value_type& kv : _annotations) rtn.push_back(kv.first); return rtn; } /// Check if an annotation is defined bool hasAnnotation(const std::string& name) const { return _annotations.find(name) != _annotations.end(); } /// Get an annotation by name (as a string) const std::string& annotation(const std::string& name) const { Annotations::const_iterator v = _annotations.find(name); // If not found... written this way round on purpose if (v == _annotations.end()) { std::string missing = "YODA::AnalysisObject: No annotation named " + name; throw AnnotationError(missing); } return v->second; } /// Get an annotation by name (as a string) with a default in case the annotation is not found const std::string& annotation(const std::string& name, const std::string& defaultreturn) const { Annotations::const_iterator v = _annotations.find(name); if (v != _annotations.end()) return v->second; return defaultreturn; } /// @brief Get an annotation by name (copied to another type) /// /// @note Templated on return type template const T annotation(const std::string& name) const { std::string s = annotation(name); return Utils::lexical_cast(s); } /// @brief Get an annotation by name (copied to another type) with a default in case the annotation is not found /// /// @note Templated on return type template const T annotation(const std::string& name, const T& defaultreturn) const { try { std::string s = annotation(name); return Utils::lexical_cast(s); } catch (const AnnotationError& ae) { return defaultreturn; } } /// @brief Add or set a string-valued annotation by name void setAnnotation(const std::string& name, const std::string& value) { _annotations[name] = value; } /// @brief Add or set a double-valued annotation by name /// @todo Can we cover all FP types in one function via SFINAE? void setAnnotation(const std::string& name, double value) { // Recipe from Boost docs std::stringstream ss; ss << std::setprecision(std::numeric_limits::max_digits10) << std::scientific << value; setAnnotation(name, ss.str()); } /// @brief Add or set a float-valued annotation by name /// @todo Can we cover all FP types in one function via SFINAE? void setAnnotation(const std::string& name, float value) { // Recipe from Boost docs std::stringstream ss; ss << std::setprecision(std::numeric_limits::max_digits10) << std::scientific << value; setAnnotation(name, ss.str()); } /// @brief Add or set a long-double-valued annotation by name /// @todo Can we cover all FP types in one function via SFINAE? void setAnnotation(const std::string& name, long double value) { // Recipe from Boost docs std::stringstream ss; ss << std::setprecision(std::numeric_limits::max_digits10) << std::scientific << value; setAnnotation(name, ss.str()); } /// @brief Add or set an annotation by name (templated for remaining types) /// /// @note Templated on arg type, but stored as a string. template void setAnnotation(const std::string& name, const T& value) { setAnnotation(name, Utils::lexical_cast(value)); } /// Set all annotations at once void setAnnotations(const Annotations& anns) { _annotations = anns; } /// @brief Add or set an annotation by name /// /// Note: Templated on arg type, but stored as a string. This is just a synonym for setAnnotation. template void addAnnotation(const std::string& name, const T& value) { setAnnotation(name, value); } /// Delete an annotation by name void rmAnnotation(const std::string& name) { _annotations.erase(name); } /// Delete an annotation by name void clearAnnotations() { _annotations.clear(); } //@} /// @name Standard annotations //@{ /// @brief Get the AO title. /// /// Returns a null string if undefined, rather than throwing an exception cf. the annotation("Title"). const std::string title() const { return annotation("Title", ""); } /// Set the AO title void setTitle(const std::string& title) { setAnnotation("Title", title); } /// @brief Get the AO path. /// /// Returns a null string if undefined, rather than throwing an exception cf. annotation("Path"). /// @note A leading / will be prepended if not already set. const std::string path() const { const std::string p = annotation("Path", ""); // If not set at all, return an empty string if (p.empty()) return p; // If missing a leading slash, one will be prepended return p.find("/") == 0 ? p : ("/"+p); } /// Set the AO path /// /// @note A leading / will be prepended if not already given. void setPath(const std::string& path) { const std::string p = (path.find("/") == 0) ? path : "/"+path; // if (path.length() > 0 && path.find("/") != 0) { // throw AnnotationError("Histo paths must start with a slash (/) character."); // } setAnnotation("Path", p); } /// Get the AO name -- the last part of the path. /// Returns a null string if path is undefined const std::string name() const { const std::string p = path(); const size_t lastslash = p.rfind("/"); if (lastslash == std::string::npos) return p; return p.substr(lastslash+1); } //@} public: /// @name Persistency hooks / object type info //@{ /// Get name of the analysis object type virtual std::string type() const { return annotation("Type"); } /// @brief Get the dimension of the analysis object type /// /// @note For fillable types this is the dimension of the fill space (e.g. Histo1D -> dim=1). /// For scatter types, it is the total dimension of the points (e.g. Scatter3D -> dim=3). virtual size_t dim() const = 0; //@} private: /// The annotations indexed by name std::map _annotations; }; // Convenience alias using AO = AnalysisObject; } #endif // YODA_AnalysisObject_h diff --git a/include/YODA/Point.h b/include/YODA/Point.h --- a/include/YODA/Point.h +++ b/include/YODA/Point.h @@ -1,112 +1,112 @@ // -*- 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_POINT_H #define YODA_POINT_H #include "YODA/AnalysisObject.h" namespace YODA { /// Base class for all Point*Ds, providing generic access to their numerical properties class Point { public: typedef std::pair ValuePair; /// Virtual destructor for inheritance virtual ~Point() {}; /// Space dimension of the point virtual size_t dim() = 0; //get the error map for the highest dimension virtual const std::map< std::string, std::pair> & errMap() const =0; - //Parse the annotation from the parent AO which contains any vaariaitoms + //Parse the annotation from the parent AO which contains any variations virtual void getVariationsFromParent() const =0; /// Get the point value for direction @a i virtual double val(size_t i) const = 0; /// Set the point value for direction @a i virtual void setVal(size_t i, double val) = 0; /// Get error values for direction @a i virtual const std::pair& errs(size_t i, std::string source="") const = 0; /// Set symmetric error for direction @a i virtual void setErr(size_t i, double e, std::string source="") = 0; /// Set symmetric error for direction @a i (alias) virtual void setErrs(size_t i, double e, std::string source="") { return setErr(i,e, source); } /// Set asymmetric error for direction @a i virtual void setErrs(size_t i, double eminus, double eplus, std::string source="") = 0; /// Set asymmetric error for direction @a i virtual void setErrs(size_t i, std::pair& e, std::string source="") = 0; /// Get negative error value for direction @a i virtual double errMinus(size_t i, std::string source="") const = 0; /// Set negative error for direction @a i virtual void setErrMinus(size_t i, double eminus, std::string source="") = 0; /// Get positive error value for direction @a i virtual double errPlus(size_t i, std::string source="") const = 0; /// Set positive error for direction @a i virtual void setErrPlus(size_t i, double eplus, std::string source="") = 0; /// Get average error value for direction @a i virtual double errAvg(size_t i, std::string source="") const = 0; // /// Get value minus negative error for direction @a i // double min(size_t i) const = 0; // /// Get value plus positive error for direction @a i // double max(size_t i) const = 0; //@} /// @todo Support multiple errors /// @name Combined value and error setters //@{ /// Set value and symmetric error for direction @a i virtual void set(size_t i, double val, double e, std::string source="") = 0; /// Set value and asymmetric error for direction @a i virtual void set(size_t i, double val, double eminus, double eplus, std::string source="") = 0; /// Set value and asymmetric error for direction @a i virtual void set(size_t i, double val, std::pair& e, std::string source="") = 0; //@} // @name Manipulations //@{ // /// Scaling of direction @a i // void scale(size_t i, double scale) = 0; /// @todo void transform(size_t i, FN f) = 0; //@} void setParentAO(AnalysisObject* parent){ _parentAO=parent; } AnalysisObject* getParentAO() const{ return _parentAO; } private: // pointer back to the parent AO which these points belong to. AnalysisObject* _parentAO=0; }; } #endif diff --git a/include/YODA/Point1D.h b/include/YODA/Point1D.h --- a/include/YODA/Point1D.h +++ b/include/YODA/Point1D.h @@ -1,360 +1,366 @@ // -*- 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 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 (!_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 (!_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 (!_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 (!_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 (!_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 (!_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/Point2D.h b/include/YODA/Point2D.h --- a/include/YODA/Point2D.h +++ b/include/YODA/Point2D.h @@ -1,576 +1,576 @@ // -*- 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_POINT2D_H #define YODA_POINT2D_H #include "YODA/Point.h" #include "YODA/Exceptions.h" #include "YODA/Utils/MathUtils.h" #include namespace YODA { /// A 2D data point to be contained in a Scatter2D class Point2D : public Point { public: /// @name Constructors //@{ // Default constructor Point2D() { } /// Constructor from values with optional symmetric errors Point2D(double x, double y, double ex=0.0, double ey=0.0, std::string source="") : _x(x), _y(y) { _ex = std::make_pair(ex, ex); _ey[source] = std::make_pair(ey, ey); } /// Constructor from values with explicit asymmetric errors Point2D(double x, double y, double exminus, double explus, double eyminus, double eyplus, std::string source="") : _x(x), _y(y) { _ex = std::make_pair(exminus, explus); _ey[source] = std::make_pair(eyminus, eyplus); } // /// Constructor from values with symmetric errors on x and asymmetric errors on y // Point2D(double x, double y, double ex, const std::pair& ey) // : _x(x), _y(y), _ey(ey) // { // _ex = std::make_pair(ex, ex); // } // /// Constructor from values with asymmetric errors on x and symmetric errors on y // Point2D(double x, double y, const std::pair& ex, double ey) // : _x(x), _y(y), _ex(ex) // { // _ey = std::make_pair(ey, ey); // } /// Constructor from values with asymmetric errors on both x and y Point2D(double x, double y, const std::pair& ex, const std::pair& ey, std::string source="") : _x(x), _y(y) { _ex = ex; _ey[source] = ey; } /// Copy constructor Point2D(const Point2D& p) : _x(p._x), _y(p._y) { _ex = p._ex; _ey = p._ey; this->setParentAO( p.getParentAO()); } /// Copy assignment Point2D& operator = (const Point2D& p) { _x = p._x; _y = p._y; _ex = p._ex; _ey = p._ey; this->setParentAO( p.getParentAO()); return *this; } //@} public: /// Space dimension of the point size_t dim() { return 2; } /// @name Value 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; } /// @todo Uniform "coords" accessor across all Scatters: returning fixed-size tuple? /// Get x,y value pair std::pair xy() const { return std::make_pair(_x, _y); } /// Set x and y values void setXY(double x, double y) { setX(x); setY(y); } /// Set x and y values void setXY(const std::pair& xy) { setX(xy.first); setY(xy.second); } //@} /// @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 /// @todo Remove (or extend) when multiple errors are supported /// No: doesn't need to change since (for now) we only store multiple /// errors for the highest dimentsion double xMin() const { return _x - _ex.first; } /// Get value plus positive x-error /// @todo Remove (or extend) when multiple errors are supported /// No: doesn't need to change since (for now) we only store multiple /// errors for the highest dimentsion double xMax() const { return _x + _ex.second; } //@} /// @name y error accessors //@{ /// Get y-error values const std::pair& yErrs(std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ey.count(source)) throw RangeError("yErrs has no such key: "+source); return _ey.at(source); } /// Get negative y-error value double yErrMinus(std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ey.count(source)) throw RangeError("yErrs has no such key: "+source); return _ey.at(source).first; } /// Get positive y-error value double yErrPlus(std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ey.count(source)) throw RangeError("yErrs has no such key: "+source); return _ey.at(source).second; } /// Get average y-error value double yErrAvg(std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ey.count(source)) throw RangeError("yErrs has no such key: "+source); double res=(fabs(_ey.at(source).first) + fabs(_ey.at(source).second))/2.; return res; } /// Set negative y error void setYErrMinus(double eyminus, std::string source="") { if (!_ey.count(source)) { _ey[source] = std::make_pair(0.,0.); } _ey.at(source).first = eyminus; } /// Set positive y error void setYErrPlus(double eyplus, std::string source="") { if (!_ey.count(source)) _ey[source] = std::make_pair(0.,0.); _ey.at(source).second = eyplus; } /// Set symmetric y error void setYErr(double ey, std::string source="") { setYErrMinus(ey, source ); setYErrPlus(ey, source ); } /// Set symmetric y error (alias) void setYErrs(double ey, std::string source="") { setYErr(ey, source); } /// Set asymmetric y error void setYErrs(double eyminus, double eyplus, std::string source="") { setYErrMinus(eyminus, source); setYErrPlus(eyplus, source ); } /// Set asymmetric y error void setYErrs(const std::pair& ey, std::string source="") { _ey[source] = ey; } /// Get value minus negative y-error double yMin(std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ey.count(source)) throw RangeError("yErrs has no such key: "+source); return _y - _ey.at(source).first; } /// Get value plus positive y-error double yMax(std::string source="") const { if (source!="") getVariationsFromParent(); if (!_ey.count(source)) throw RangeError("yErrs has no such key: "+source); return _y + _ey.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, std::string source="") { setY(y); setYErrs(ey, source); } /// Set y value and asymmetric error void setY(double y, double eyminus, double eyplus, std::string source="") { setY(y); setYErrs(eyminus, eyplus, source); } /// Set y value and asymmetric error void setY(double y, std::pair& ey, std::string source="") { setY(y); setYErrs(ey, 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); for (const auto &source : _ey){ setYErrs(yErrMinus()*scaley, yErrPlus()*scaley, source.first); } } /// Scaling of both axes void scaleXY(double scalex, double scaley) { scaleX(scalex); scaleY(scaley); } /// Scaling of both axes /// @deprecated Use scaleXY void scale(double scalex, double scaley) { scaleXY(scalex, scaley); } //@} /// @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(); 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; 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; - // Get parse the variations from the parent AO if it exists + // 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(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(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(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(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, 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, 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, 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, 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, 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, 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, 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, source); break; default: throw RangeError("Invalid axis int, must be in range 1..dim"); } } //@} protected: /// @name Value and error variables //@{ double _x; double _y; std::pair _ex; // a map of the errors for each source. Nominal stored under "" // to ensure backward compatibility std::map< std::string, std::pair > _ey; //@} }; /// @name Comparison operators //@{ /// Equality test of x & y characteristics only /// @todo Base on a named fuzzyEquals(a,b,tol=1e-3) unbound function inline bool operator==(const YODA::Point2D& a, const YODA::Point2D& 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; return true; } /// Equality test of x characteristics only inline bool operator != (const YODA::Point2D& a, const YODA::Point2D& b) { return !(a == b); } /// Less-than operator used to sort bins by x-ordering inline bool operator < (const YODA::Point2D& a, const YODA::Point2D& 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::Point2D& a, const YODA::Point2D& b) { if (a == b) return true; return a < b; } /// Greater-than operator used to sort bins by x-ordering inline bool operator > (const YODA::Point2D& a, const YODA::Point2D& b) { return !(a <= b); } /// Greater-than-or-equals operator used to sort bins by x-ordering inline bool operator >= (const YODA::Point2D& a, const YODA::Point2D& 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,683 +1,695 @@ // -*- 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/include/YODA/Scatter1D.h b/include/YODA/Scatter1D.h --- a/include/YODA/Scatter1D.h +++ b/include/YODA/Scatter1D.h @@ -1,334 +1,342 @@ // -*- 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_SCATTER1D_H #define YODA_SCATTER1D_H #include "YODA/AnalysisObject.h" #include "YODA/Point1D.h" #include "YODA/Utils/sortedvector.h" #include #include namespace YODA { // Forward declarations class Counter; /// A very generic data type which is just a collection of 1D data points with errors class Scatter1D : public AnalysisObject { public: /// Type of the native Point1D collection typedef Point1D Point; typedef Utils::sortedvector Points; typedef std::shared_ptr Ptr; /// @name Constructors //@{ /// Empty constructor Scatter1D(const std::string& path="", const std::string& title="") : AnalysisObject("Scatter1D", path, title) { } /// Constructor from a set of points Scatter1D(const Points& points, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter1D", path, title), _points(points) { } /// Constructor from a vector of x values with no errors Scatter1D(const std::vector& x, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter1D", path, title) { for (size_t i = 0; i < x.size(); ++i) addPoint(x[i]); } /// Constructor from vectors of x values with symmetric errors Scatter1D(const std::vector& x, const std::vector& ex, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter1D", path, title) { if (x.size() != ex.size()) throw UserError("x and ex vectors must have same length"); for (size_t i = 0; i < x.size(); ++i) addPoint(x[i], ex[i]); } /// Constructor from x values with asymmetric errors Scatter1D(const std::vector& x, const std::vector >& ex, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter1D", path, title) { if (x.size() != ex.size()) throw UserError("x and ex vectors must have same length"); for (size_t i = 0; i < x.size(); ++i) addPoint(Point1D(x[i], ex[i])); } /// Constructor from values with completely explicit asymmetric errors Scatter1D(const std::vector& x, const std::vector& exminus, const std::vector& explus, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter1D", path, title) { if (x.size() != exminus.size()) throw UserError("x and ex vectors must have same length"); if (exminus.size() != explus.size()) throw UserError("ex plus and minus vectors must have same length"); for (size_t i = 0; i < x.size(); ++i) addPoint(Point1D(x[i], exminus[i], explus[i])); } /// Copy constructor with optional new path /// @todo Also allow title setting from the constructor? Scatter1D(const Scatter1D& s1, const std::string& path="") : AnalysisObject("Scatter1D", (path.size() == 0) ? s1.path() : path, s1, s1.title()), _points(s1._points) { for ( auto &ann : annotations()){ setAnnotation(ann, annotation(ann)); } } /// Assignment operator Scatter1D& operator = (const Scatter1D& s1) { AnalysisObject::operator = (s1); //< AO treatment of paths etc. _points = s1._points; return *this; } /// Make a copy on the stack Scatter1D clone() const { return Scatter1D(*this); } /// Make a copy on the heap, via 'new' Scatter1D* newclone() const { return new Scatter1D(*this); } //@} /// Dimension of this data object size_t dim() const { return 1; } /// @name Modifiers //@{ /// Clear all points void reset() { _points.clear(); } /// Scaling of x axis void scaleX(double scalex) { for (Point1D& p : _points) p.scaleX(scalex); } //@} /////////////////////////////////////////////////// void parseVariations() ; /// Get the list of variations stored in the points const std::vector variations() const ; /// @name Point accessors //@{ /// Number of points in the scatter size_t numPoints() const { return _points.size(); } /// Get the collection of points (non-const) Points& points() { return _points; } /// Get the collection of points (const) const Points& points() const { return _points; } /// Get a reference to the point with index @a index (non-const) Point1D& point(size_t index) { if (index >= numPoints()) throw RangeError("There is no point with this index"); return _points.at(index); } /// Get a reference to the point with index @a index (const) const Point1D& point(size_t index) const { if (index >= numPoints()) throw RangeError("There is no point with this index"); return _points.at(index); } //@} /// @name Point inserters //@{ /// Insert a new point void addPoint(const Point1D& pt) { _points.insert(pt); } /// Insert a new point, defined as the x value and no errors void addPoint(double x) { - _points.insert(Point1D(x)); + Point1D thisPoint=Point1D(x); + thisPoint.setParentAO(this); + _points.insert(thisPoint); } /// Insert a new point, defined as the x value and symmetric errors void addPoint(double x, double ex) { - _points.insert(Point1D(x, ex)); + Point1D thisPoint=Point1D(x, ex); + thisPoint.setParentAO(this); + _points.insert(thisPoint); } /// Insert a new point, defined as the x value and an asymmetric error pair void addPoint(double x, const std::pair& ex) { - _points.insert(Point1D(x, ex)); + Point1D thisPoint=Point1D(x, ex); + thisPoint.setParentAO(this); + _points.insert(thisPoint); } /// Insert a new point, defined as the x value and explicit asymmetric errors void addPoint(double x, double exminus, double explus) { - _points.insert(Point1D(x, exminus, explus)); + Point1D thisPoint=Point1D(x, exminus, explus); + thisPoint.setParentAO(this); + _points.insert(thisPoint); } /// Insert a collection of new points void addPoints(const Points& pts) { for (const Point1D& pt : pts) addPoint(pt); } //@} /// @name Combining sets of scatter points //@{ /// @todo Better name? void combineWith(const Scatter1D& other) { addPoints(other.points()); } /// @todo Better name? Make this the add operation? /// @todo Convert/extend to accept a Range or generic void combineWith(const std::vector& others) { for (const Scatter1D& s : others) combineWith(s); } //@} /// Equality operator bool operator == (const Scatter1D& other) { return _points == other._points; } /// Non-equality operator bool operator != (const Scatter1D& other) { return ! operator == (other); } ////////////////////////////////// private: Points _points; bool _variationsParsed =false ; }; /// Convenience typedef typedef Scatter1D S1D; /// @name Combining scatters by merging sets of points //@{ inline Scatter1D combine(const Scatter1D& a, const Scatter1D& b) { Scatter1D rtn = a; rtn.combineWith(b); return rtn; } inline Scatter1D combine(const std::vector& scatters) { Scatter1D rtn; rtn.combineWith(scatters); return rtn; } //@} ////////////////////////////////// /// @name Conversion functions from other data types //@{ /// Make a Scatter1D representation of a Histo1D Scatter1D mkScatter(const Counter& c); /// Make a Scatter1D representation of... erm, a Scatter1D! /// @note Mainly exists to allow mkScatter to be called on any AnalysisObject type inline Scatter1D mkScatter(const Scatter1D& s) { return Scatter1D(s); } //@} ///////////////////////////////// /// @name Transforming operations on Scatter1D //@{ /// @brief Apply transformation fx(x) to all values and error positions (operates in-place on @a s) /// /// fx should be a function which takes double x -> double newx template inline void transformX(Scatter1D& s, FNX fx) { for (size_t i = 0; i < s.numPoints(); ++i) { Point1D& p = s.point(i); const double newx = fx(p.x()); const double fx_xmin = fx(p.xMin()); const double fx_xmax = fx(p.xMax()); // Deal with possible inversions of min/max ordering under the transformation const double newxmin = std::min(fx_xmin, fx_xmax); const double newxmax = std::max(fx_xmin, fx_xmax); // Set new point x values p.setX(newx); /// @todo Be careful about transforms which could switch around min and max errors, or send both in the same direction! p.setXErrMinus(newx - newxmin); p.setXErrPlus(newxmax - newx); } } //@} } #endif diff --git a/include/YODA/Scatter2D.h b/include/YODA/Scatter2D.h --- a/include/YODA/Scatter2D.h +++ b/include/YODA/Scatter2D.h @@ -1,440 +1,439 @@ // -*- 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_SCATTER2D_H #define YODA_SCATTER2D_H #include "YODA/AnalysisObject.h" #include "YODA/Point2D.h" #include "YODA/Utils/sortedvector.h" #include #include namespace YODA { // Forward declarations class Histo1D; class Profile1D; /// A very generic data type which is just a collection of 2D data points with errors class Scatter2D : public AnalysisObject { public: /// Type of the native Point2D collection typedef Point2D Point; typedef Utils::sortedvector Points; typedef std::shared_ptr Ptr; /// @name Constructors //@{ /// Empty constructor Scatter2D(const std::string& path="", const std::string& title="") : AnalysisObject("Scatter2D", path, title) { } /// Constructor from a set of points Scatter2D(const Points& points, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter2D", path, title), _points(points) { } /// Constructor from a vector of values with no errors Scatter2D(const std::vector& x, const std::vector& y, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter2D", path, title) { if (x.size() != y.size()) throw UserError("x and y vectors must have same length"); for (size_t i = 0; i < x.size(); ++i) addPoint(x[i], y[i]); } /// Constructor from vectors of values with symmetric errors on x and y Scatter2D(const std::vector& x, const std::vector& y, const std::vector& ex, const std::vector& ey, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter2D", path, title) { if (x.size() != y.size()) throw UserError("x and y vectors must have same length"); if (x.size() != ex.size()) throw UserError("x and ex vectors must have same length"); if (y.size() != ey.size()) throw UserError("y and ey vectors must have same length"); for (size_t i = 0; i < x.size(); ++i) addPoint(x[i], y[i], ex[i], ey[i]); } /// Constructor from values with asymmetric errors on both x and y Scatter2D(const std::vector& x, const std::vector& y, const std::vector >& ex, const std::vector >& ey, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter2D", path, title) { if (x.size() != y.size()) throw UserError("x and y vectors must have same length"); if (x.size() != ex.size()) throw UserError("x and ex vectors must have same length"); if (y.size() != ey.size()) throw UserError("y and ey vectors must have same length"); for (size_t i = 0; i < x.size(); ++i) addPoint(Point2D(x[i], y[i], ex[i], ey[i])); } /// Constructor from values with completely explicit asymmetric errors Scatter2D(const std::vector& x, const std::vector& y, const std::vector& exminus, const std::vector& explus, const std::vector& eyminus, const std::vector& eyplus, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter2D", path, title) { if (x.size() != y.size()) throw UserError("x and y vectors must have same length"); if (x.size() != exminus.size()) throw UserError("x and ex vectors must have same length"); if (y.size() != eyminus.size()) throw UserError("y and ey vectors must have same length"); if (exminus.size() != explus.size()) throw UserError("ex plus and minus vectors must have same length"); if (eyminus.size() != eyplus.size()) throw UserError("ey plus and minus vectors must have same length"); for (size_t i = 0; i < x.size(); ++i) addPoint(Point2D(x[i], y[i], exminus[i], explus[i], eyminus[i], eyplus[i])); } /// Copy constructor with optional new path /// @todo Also allow title setting from the constructor? Scatter2D(const Scatter2D& s2, const std::string& path="") : AnalysisObject("Scatter2D", (path.size() == 0) ? s2.path() : path, s2, s2.title()), _points(s2._points) { for ( auto &ann : annotations()){ setAnnotation(ann, annotation(ann)); } } /// Assignment operator Scatter2D& operator = (const Scatter2D& s2) { AnalysisObject::operator = (s2); //< AO treatment of paths etc. _points = s2._points; return *this; } /// Make a copy on the stack Scatter2D clone() const { return Scatter2D(*this); } /// Make a copy on the heap, via 'new' Scatter2D* newclone() const { return new Scatter2D(*this); } //@} /// Dimension of this data object size_t dim() const { return 2; } /// @name Modifiers //@{ /// Clear all points void reset() { _points.clear(); } /// Scaling of x axis void scaleX(double scalex) { for (Point2D& p : _points) p.scaleX(scalex); } /// Scaling of y axis void scaleY(double scaley) { for (Point2D& p : _points) p.scaleY(scaley); } /// Scaling of both axes void scaleXY(double scalex, double scaley) { for (Point2D& p : _points) p.scaleXY(scalex, scaley); } /// Scaling of both axes /// @deprecated Use scaleXY void scale(double scalex, double scaley) { scaleXY(scalex, scaley); } //@} /////////////////////////////////////////////////// void parseVariations() ; /// Get the list of variations stored in the points const std::vector variations() const; // Construct a covariance matrix from the error breakdown std::vector > covarianceMatrix(bool ignoreOffDiagonalTerms=false) ; /// @name Point accessors //@{ /// Number of points in the scatter size_t numPoints() const { return _points.size(); } /// Get the collection of points (non-const) Points& points() { return _points; } /// Get the collection of points (const) const Points& points() const { return _points; } /// Get a reference to the point with index @a index (non-const) Point2D& point(size_t index) { if (index >= numPoints()) throw RangeError("There is no point with this index"); return _points.at(index); } /// Get a reference to the point with index @a index (const) const Point2D& point(size_t index) const { if (index >= numPoints()) throw RangeError("There is no point with this index"); return _points.at(index); } //@} /// @name Point inserters //@{ /// Insert a new point void addPoint(const Point2D& pt) { - //pt.setParentAO(this); // how to avoid const-ness ? _points.insert(pt); } /// Insert a new point, defined as the x/y value pair and no errors void addPoint(double x, double y) { Point2D thisPoint= Point2D(x, y); + thisPoint.setParentAO(this); _points.insert(thisPoint); - thisPoint.setParentAO(this); } /// Insert a new point, defined as the x/y value pair and symmetric errors void addPoint(double x, double y, double ex, double ey) { Point2D thisPoint= Point2D(x, y, ex, ey); + thisPoint.setParentAO(this); _points.insert(thisPoint); - thisPoint.setParentAO(this); } /// Insert a new point, defined as the x/y value pair and asymmetric error pairs void addPoint(double x, double y, const std::pair& ex, const std::pair& ey) { Point2D thisPoint= Point2D(x, y, ex, ey); + thisPoint.setParentAO(this); _points.insert(thisPoint); - thisPoint.setParentAO(this); } /// Insert a new point, defined as the x/y value pair and asymmetric errors void addPoint(double x, double y, double exminus, double explus, double eyminus, double eyplus) { Point2D thisPoint=Point2D(x, y, exminus, explus, eyminus, eyplus); + thisPoint.setParentAO(this); _points.insert(thisPoint); - thisPoint.setParentAO(this); } /// Insert a collection of new points void addPoints(const Points& pts) { for (const Point2D& pt : pts) { addPoint(pt); } } //@} /// @name Combining sets of scatter points //@{ /// @todo Better name? Make this the add operation? void combineWith(const Scatter2D& other) { addPoints(other.points()); //return *this; } /// @todo Better name? /// @todo Convert/extend to accept a Range or generic void combineWith(const std::vector& others) { for (const Scatter2D& s : others) combineWith(s); //return *this; } //@} /// Equality operator bool operator == (const Scatter2D& other) { return _points == other._points; } /// Non-equality operator bool operator != (const Scatter2D& other) { return ! operator == (other); } ////////////////////////////////// private: Points _points; bool _variationsParsed =false ; }; /// Convenience typedef typedef Scatter2D S2D; /// @name Combining scatters by merging sets of points //@{ inline Scatter2D combine(const Scatter2D& a, const Scatter2D& b) { Scatter2D rtn = a; rtn.combineWith(b); return rtn; } inline Scatter2D combine(const std::vector& scatters) { Scatter2D rtn; rtn.combineWith(scatters); return rtn; } //@} ////////////////////////////////// /// @name Conversion functions from other data types //@{ /// Make a Scatter2D representation of a Histo1D /// /// Optional @c usefocus argument can be used to position the point at the bin /// focus rather than geometric midpoint. Scatter2D mkScatter(const Histo1D& h, bool usefocus=false); /// Make a Scatter2D representation of a Profile1D /// /// Optional @c usefocus argument can be used to position the point at the bin /// focus rather than geometric midpoint. Optional @c usestddev argument can /// be used to draw the y-distribution sigma rather than the standard error on /// the mean as the y-error bar size. Scatter2D mkScatter(const Profile1D& p, bool usefocus=false, bool usestddev=false); /// Make a Scatter2D representation of... erm, a Scatter2D! /// @note Mainly exists to allow mkScatter to be called on any AnalysisObject type inline Scatter2D mkScatter(const Scatter2D& s) { return Scatter2D(s); } // /// @note The usefocus arg is just for consistency and has no effect for Scatter -> Scatter // inline Scatter2D mkScatter(const Scatter2D& s, bool) { return mkScatter(s); } //@} ////////////////////////////////// /// @name Transforming operations on Scatter2D //@{ /// @brief Apply transformation fx(x) to all values and error positions (operates in-place on @a s) /// /// fx should be a function which takes double x -> double newx template inline void transformX(Scatter2D& s, FNX fx) { for (size_t i = 0; i < s.numPoints(); ++i) { Point2D& p = s.point(i); const double newx = fx(p.x()); const double fx_xmin = fx(p.xMin()); const double fx_xmax = fx(p.xMax()); // Deal with possible inversions of min/max ordering under the transformation const double newxmin = std::min(fx_xmin, fx_xmax); const double newxmax = std::max(fx_xmin, fx_xmax); // Set new point x values p.setX(newx); /// @todo Be careful about transforms which could switch around min and max errors, or send both in the same direction! p.setXErrMinus(newx - newxmin); p.setXErrPlus(newxmax - newx); } } /// @brief Apply transformation fy(y) to all values and error positions (operates in-place on @a s) /// /// fy should be a function which takes double y -> double newy template inline void transformY(Scatter2D& s, FNY fy) { for (size_t i = 0; i < s.numPoints(); ++i) { Point2D& p = s.point(i); const double newy = fy(p.y()); const double fy_ymin = fy(p.yMin()); const double fy_ymax = fy(p.yMax()); // Deal with possible inversions of min/max ordering under the transformation const double newymin = std::min(fy_ymin, fy_ymax); const double newymax = std::max(fy_ymin, fy_ymax); // Set new point y values p.setY(newy); /// @todo Be careful about transforms which could switch around min and max errors, or send both in the same direction! p.setYErrMinus(newy - newymin); p.setYErrPlus(newymax - newy); } } /// @todo Add external scale, scaleX, scaleY functions /// Exchange the x and y axes (operates in-place on @a s) inline void flip(Scatter2D& s) { for (size_t i = 0; i < s.numPoints(); ++i) { Point2D& p = s.point(i); const double newx = p.y(); const double newy = p.x(); const double newxmin = p.yMin(); const double newxmax = p.yMax(); const double newymin = p.xMin(); const double newymax = p.xMax(); p.setX(newx); p.setY(newy); /// @todo Be careful about transforms which could switch around min and max errors, or send both in the same direction! p.setXErrMinus(newx - newxmin); p.setXErrPlus(newxmax - newx); p.setYErrMinus(newy - newymin); p.setYErrPlus(newymax - newy); } } //@} } #endif diff --git a/include/YODA/Scatter3D.h b/include/YODA/Scatter3D.h --- a/include/YODA/Scatter3D.h +++ b/include/YODA/Scatter3D.h @@ -1,441 +1,449 @@ // -*- 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_SCATTER3D_H #define YODA_SCATTER3D_H #include "YODA/AnalysisObject.h" #include "YODA/Point3D.h" #include "YODA/Utils/sortedvector.h" #include #include namespace YODA { // Forward declarations class Histo2D; class Profile2D; /// A very generic data type which is just a collection of 3D data points with errors class Scatter3D : public AnalysisObject { public: /// Types of the native Point3D collection typedef Point3D Point; typedef Utils::sortedvector Points; typedef std::shared_ptr Ptr; /// @name Constructors //@{ /// Empty constructor Scatter3D(const std::string& path="", const std::string& title="") : AnalysisObject("Scatter3D", path, title) { } /// Constructor from a set of points Scatter3D(const Points& points, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter3D", path, title), _points(points) { std::sort(_points.begin(), _points.end()); } /// Constructor from vectors of values with no errors Scatter3D(const std::vector& x, const std::vector& y, const std::vector& z, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter3D", path, title) { if (x.size() != y.size() || y.size() != z.size()) { throw RangeError("There are different numbers of x, y, and z values in the provided vectors."); } const std::pair nullerr = std::make_pair(0.0, 0.0); for (size_t i = 0; i < x.size(); ++i) { addPoint(Point3D(x[i], y[i], z[i], nullerr, nullerr, nullerr)); } std::sort(_points.begin(), _points.end()); } /// Constructor from vectors of values with asymmetric errors on both x and y Scatter3D(const std::vector& x, const std::vector& y, const std::vector& z, const std::vector >& ex, const std::vector >& ey, const std::vector >& ez, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter3D", path, title) { if (x.size() != y.size() || y.size() != z.size()) { throw RangeError("There are different numbers of x, y, and z values in the provided vectors."); } if (x.size() != ex.size() || y.size() != ey.size() || z.size() != ez.size()) { throw RangeError("The sizes of the provided error vectors don't match the corresponding x, y, or z value vectors."); } for (size_t i = 0; i < x.size(); ++i) { addPoint(Point3D(x[i], y[i], z[i], ex[i], ey[i], ez[i])); } std::sort(_points.begin(), _points.end()); } /// Constructor from vectors of values with completely explicit asymmetric errors Scatter3D(const std::vector& x, const std::vector& y, const std::vector z, const std::vector& exminus, const std::vector& explus, const std::vector& eyminus, const std::vector& eyplus, const std::vector& ezminus, const std::vector& ezplus, const std::string& path="", const std::string& title="") : AnalysisObject("Scatter3D", path, title) { if(x.size() != y.size() || y.size() != z.size() || x.size() != exminus.size() || x.size() != explus.size() || y.size() != eyminus.size() || y.size() != eyplus.size() || z.size() != ezminus.size() || z.size() != ezplus.size()) throw RangeError("There are either different amounts of points on x/y/z vectors or not every of these vectors has properly defined error vectors!"); for (size_t i = 0; i < x.size(); ++i) { addPoint(Point3D(x[i], y[i], z[i], exminus[i], explus[i], eyminus[i], eyplus[i], ezminus[i], ezplus[i])); } std::sort(_points.begin(), _points.end()); } /// Copy constructor with optional new path /// @todo Also allow title setting from the constructor? Scatter3D(const Scatter3D& s3, const std::string& path="") : AnalysisObject("Scatter3D", (path.size() == 0) ? s3.path() : path, s3, s3.title()), _points(s3._points) { for ( auto &ann : annotations()){ setAnnotation(ann, annotation(ann)); } } /// Assignment operator Scatter3D& operator = (const Scatter3D& s3) { AnalysisObject::operator = (s3); //< AO treatment of paths etc. _points = s3._points; return *this; } /// Make a copy on the stack Scatter3D clone() const { return Scatter3D(*this); } /// Make a copy on the heap, via 'new' Scatter3D* newclone() const { return new Scatter3D(*this); } //@} /// Dimension of this data object size_t dim() const { return 3; } /// @name Modifiers //@{ /// Clear all points void reset() { _points.clear(); } /// Scaling of x axis void scaleX(double scalex) { for (Point3D& p : _points) p.scaleX(scalex); } /// Scaling of y axis void scaleY(double scaley) { for (Point3D& p : _points) p.scaleY(scaley); } /// Scaling of z axis void scaleZ(double scalez) { for (Point3D& p : _points) p.scaleZ(scalez); } /// Scaling of all three axes void scaleXYZ(double scalex, double scaley, double scalez) { for (Point3D& p : _points) p.scaleXYZ(scalex, scaley, scalez); } /// Scaling of all three axes /// @deprecated Use scaleXYZ void scale(double scalex, double scaley, double scalez) { scaleXYZ(scalex, scaley, scalez); } //@} /////////////////////////////////////////////////// void parseVariations() ; /// Get the list of variations stored in the points const std::vector variations() const; /// @name Point accessors //@{ /// Number of points in the scatter size_t numPoints() const { return _points.size(); } /// Get the collection of points (non-const) Points& points() { return _points; } /// Get the collection of points (const) const Points& points() const { return _points; } /// Get a reference to the point with index @a index Point3D& point(size_t index) { if (index >= numPoints()) throw RangeError("There is no point with this index"); return _points.at(index); } /// Get the point with index @a index (const version) const Point3D& point(size_t index) const { if (index >= numPoints()) throw RangeError("There is no point with such index!"); return _points.at(index); } //@} /// @name Point inserters //@{ /// Insert a new point void addPoint(const Point3D& pt) { _points.insert(pt); } /// Insert a new point, defined as the x/y/z value triplet and no errors void addPoint(double x, double y, double z) { - _points.insert(Point3D(x, y, z)); + Point3D thisPoint=Point3D(x, y, z); + thisPoint.setParentAO(this); + _points.insert(thisPoint); } /// Insert a new point, defined as the x/y/z value triplet and symmetric errors void addPoint(double x, double y, double z, double ex, double ey, double ez) { - _points.insert(Point3D(x, y, z, ex, ey, ez)); + Point3D thisPoint=Point3D(x, y, z, ex, ey, ez); + thisPoint.setParentAO(this); + _points.insert(thisPoint); } /// Insert a new point, defined as the x/y/z value triplet and asymmetric error pairs void addPoint(double x, double y, double z, const std::pair& ex, const std::pair& ey, const std::pair& ez) { - _points.insert(Point3D(x, y, z, ex, ey, ez)); + Point3D thisPoint= Point3D(x, y, z, ex, ey, ez); + thisPoint.setParentAO(this); + _points.insert(thisPoint); } /// Insert a new point, defined as the x/y/z value triplet and asymmetric errors void addPoint(double x, double y, double z, double exminus, double explus, double eyminus, double eyplus, double ezminus, double ezplus) { - _points.insert(Point3D(x, y, z, exminus, explus, eyminus, eyplus, ezminus, ezplus)); + Point3D thisPoint = Point3D(x, y, z, exminus, explus, eyminus, eyplus, ezminus, ezplus); + thisPoint.setParentAO(this); + _points.insert(thisPoint); } /// Insert a collection of new points void addPoints(const Points& pts) { for (const Point3D& pt : pts) addPoint(pt); } //@} /// @todo Better name? void combineWith(const Scatter3D& other) { addPoints(other.points()); //return *this; } /// @todo Better name? /// @todo Convert to accept a Range or generic void combineWith(const std::vector& others) { for (const Scatter3D& s : others) combineWith(s); } /// Equality operator bool operator == (const Scatter3D& other) { return _points == other._points; } /// Non-equality operator bool operator != (const Scatter3D& other) { return ! operator == (other); } ////////////////////////////////// private: Points _points; bool _variationsParsed =false ; }; /// Convenience typedef typedef Scatter3D S3D; /// @name Combining scatters by merging sets of points //@{ inline Scatter3D combine(const Scatter3D& a, const Scatter3D& b) { Scatter3D rtn = a; rtn.combineWith(b); return rtn; } inline Scatter3D combine(const std::vector& scatters) { Scatter3D rtn; rtn.combineWith(scatters); return rtn; } //@} ////////////////////////////////// /// @name Conversion functions from other data types //@{ /// Make a Scatter3D representation of a Histo2D /// /// Optional @c usefocus argument can be used to position the point at the bin /// focus rather than geometric midpoint. Scatter3D mkScatter(const Histo2D& h, bool usefocus=false); /// Make a Scatter3D representation of a Profile2D /// /// Optional @c usefocus argument can be used to position the point at the bin /// focus rather than geometric midpoint. Optional @c usestddev argument can /// be used to draw the distribution sigma rather than the standard error on /// the mean as the z-error bar size. Scatter3D mkScatter(const Profile2D& p, bool usefocus=false, bool usestddev=false); /// Make a Scatter3D representation of... erm, a Scatter3D! /// @note Mainly exists to allow mkScatter to be called on any AnalysisObject type inline Scatter3D mkScatter(const Scatter3D& s) { return Scatter3D(s); } // /// @note The usefocus arg is just for consistency and has no effect for Scatter -> Scatter //inline Scatter3D mkScatter(const Scatter3D& s, bool) { return mkScatter(s); } //@} ///////////////////////////////// /// @name Transforming operations on Scatter3D //@{ /// @brief Apply transformation fx(x) to all values and error positions (operates in-place on @a s) /// /// fx should be a function which takes double x -> double newx template inline void transformX(Scatter3D& s, FNX fx) { for (size_t i = 0; i < s.numPoints(); ++i) { Point3D& p = s.point(i); const double newx = fx(p.x()); const double fx_xmin = fx(p.xMin()); const double fx_xmax = fx(p.xMax()); // Deal with possible inversions of min/max ordering under the transformation const double newxmin = std::min(fx_xmin, fx_xmax); const double newxmax = std::max(fx_xmin, fx_xmax); // Set new point x values p.setX(newx); /// @todo Be careful about transforms which could switch around min and max errors, or send both in the same direction! p.setXErrMinus(newx - newxmin); p.setXErrPlus(newxmax - newx); } } /// @brief Apply transformation fy(y) to all values and error positions (operates in-place on @a s) /// /// fy should be a function which takes double y -> double newy template inline void transformY(Scatter3D& s, FNY fy) { for (size_t i = 0; i < s.numPoints(); ++i) { Point3D& p = s.point(i); const double newy = fy(p.y()); const double fy_ymin = fy(p.yMin()); const double fy_ymax = fy(p.yMax()); // Deal with possible inversions of min/max ordering under the transformation const double newymin = std::min(fy_ymin, fy_ymax); const double newymax = std::max(fy_ymin, fy_ymax); // Set new point y values p.setY(newy); /// @todo Be careful about transforms which could switch around min and max errors, or send both in the same direction! p.setYErrMinus(newy - newymin); p.setYErrPlus(newymax - newy); } } /// @brief Apply transformation fz(z) to all values and error positions (operates in-place on @a s) /// /// fz should be a function which takes double z -> double newz template inline void transformZ(Scatter3D& s, FNZ fz) { for (size_t i = 0; i < s.numPoints(); ++i) { Point3D& p = s.point(i); const double newz = fz(p.z()); const double fz_zmin = fz(p.zMin()); const double fz_zmax = fz(p.zMax()); // Deal with possible inversions of min/max ordering under the transformation const double newzmin = std::min(fz_zmin, fz_zmax); const double newzmax = std::max(fz_zmin, fz_zmax); // Set new point z values p.setZ(newz); /// @todo Be careful about transforms which could switch around min and max errors, or send both in the same direction! p.setZErrMinus(newz - newzmin); p.setZErrPlus(newzmax - newz); } } /// @todo Add external scale, scaleX, scaleY, scaleZ functions //@} } #endif diff --git a/src/Makefile.am b/src/Makefile.am --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,33 +1,33 @@ SUBDIRS = tinyxml yamlcpp . lib_LTLIBRARIES = libYODA.la libYODA_la_SOURCES = \ Exceptions.cc \ Reader.cc \ ReaderYODA.cc \ ReaderFLAT.cc \ ReaderAIDA.cc \ Writer.cc \ WriterYODA.cc \ WriterFLAT.cc \ WriterAIDA.cc \ - Dbn0D.cc \ + Dbn0D.cc \ Dbn1D.cc \ Counter.cc \ Histo1D.cc \ Histo2D.cc \ Profile1D.cc \ Profile2D.cc \ Scatter1D.cc \ Scatter2D.cc \ Scatter3D.cc \ Point1D.cc \ Point2D.cc \ Point3D.cc libYODA_la_LDFLAGS = -avoid-version libYODA_la_LIBADD = $(builddir)/tinyxml/libyoda-tinyxml.la $(builddir)/yamlcpp/libyoda-yaml-cpp.la libYODA_la_CPPFLAGS = $(AM_CPPFLAGS) -DTIXML_USE_STL -I$(srcdir)/yamlcpp -I$(srcdir) -DYAML_NAMESPACE=YODA_YAML EXTRA_DIST = zstr diff --git a/src/ReaderYODA.cc b/src/ReaderYODA.cc --- a/src/ReaderYODA.cc +++ b/src/ReaderYODA.cc @@ -1,526 +1,528 @@ // -*- 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/ReaderYODA.h" #include "YODA/Utils/StringUtils.h" #include "YODA/Utils/getline.h" #include "YODA/Exceptions.h" #include "YODA/Config/DummyConfig.h" #include "YODA/Counter.h" #include "YODA/Histo1D.h" #include "YODA/Histo2D.h" #include "YODA/Profile1D.h" #include "YODA/Profile2D.h" #include "YODA/Scatter1D.h" #include "YODA/Scatter2D.h" #include "YODA/Scatter3D.h" #include "yaml-cpp/yaml.h" #ifdef YAML_NAMESPACE #define YAML YAML_NAMESPACE #endif #ifdef HAVE_LIBZ #define _XOPEN_SOURCE 700 #include "zstr/zstr.hpp" #endif #include #include using namespace std; namespace YODA { /// Singleton creation function Reader& ReaderYODA::create() { static ReaderYODA _instance; return _instance; } namespace { /// Fast ASCII tokenizer, extended from FastIStringStream by Gavin Salam. class aistringstream { public: // Constructor from char* aistringstream(const char* line=0) { reset(line); } // Constructor from std::string aistringstream(const string& line) { reset(line); } // Re-init to new line as char* void reset(const char* line=0) { _next = const_cast(line); _new_next = _next; _error = false; } // Re-init to new line as std::string void reset(const string& line) { reset(line.c_str()); } // Tokenizing stream operator (forwards to specialisations) template aistringstream& operator >> (T& value) { _get(value); if (_new_next == _next) _error = true; // handy error condition behaviour! _next = _new_next; return *this; } // Allow use of operator>> in a while loop operator bool() const { return !_error; } private: void _get(double& x) { x = std::strtod(_next, &_new_next); } void _get(float& x) { x = std::strtof(_next, &_new_next); } void _get(int& i) { i = std::strtol(_next, &_new_next, 10); } // force base 10! void _get(long& i) { i = std::strtol(_next, &_new_next, 10); } // force base 10! void _get(unsigned int& i) { i = std::strtoul(_next, &_new_next, 10); } // force base 10! void _get(long unsigned int& i) { i = std::strtoul(_next, &_new_next, 10); } // force base 10! void _get(string& x) { /// @todo If _next and _new_next become null? while (std::isspace(*_next)) _next += 1; _new_next = _next; while (!std::isspace(*_new_next)) _new_next += 1; x = string(_next, _new_next-_next); } char *_next, *_new_next; bool _error; }; } void ReaderYODA::read(istream& stream_, vector& aos) { #ifdef HAVE_LIBZ // NB. zstr auto-detects if file is deflated or plain-text zstr::istream stream(stream_); #else istream& stream = stream_; #endif // Data format parsing states, representing current data type /// @todo Extension to e.g. "bar" or multi-counter or binned-value types, and new formats for extended Scatter types enum Context { NONE, //< outside any data block SCATTER1D, SCATTER2D, SCATTER3D, COUNTER, HISTO1D, HISTO2D, PROFILE1D, PROFILE2D }; /// State of the parser: line number, line, parser context, and pointer(s) to the object currently being assembled unsigned int nline = 0; string s; Context context = NONE; // AnalysisObject* aocurr = NULL; //< Generic current AO pointer vector h1binscurr; //< Current H1 bins container vector h2binscurr; //< Current H2 bins container vector p1binscurr; //< Current P1 bins container vector p2binscurr; //< Current P2 bins container vector pt1scurr; //< Current Point1Ds container vector pt2scurr; //< Current Point2Ds container vector pt3scurr; //< Current Point3Ds container Counter* cncurr = NULL; Histo1D* h1curr = NULL; Histo2D* h2curr = NULL; Profile1D* p1curr = NULL; Profile2D* p2curr = NULL; Scatter1D* s1curr = NULL; Scatter2D* s2curr = NULL; Scatter3D* s3curr = NULL; //std::vector variationscurr; string annscurr; // Loop over all lines of the input file aistringstream aiss; bool in_anns = false; string fmt = "1"; //int nfmt = 1; while (Utils::getline(stream, s)) { nline += 1; // CLEAN LINES IF NOT IN ANNOTATION MODE if (!in_anns) { // Trim the line Utils::itrim(s); // Ignore blank lines if (s.empty()) continue; // Ignore comments (whole-line only, without indent, and still allowed for compatibility on BEGIN/END lines) if (s.find("#") == 0 && s.find("BEGIN") == string::npos && s.find("END") == string::npos) continue; } // STARTING A NEW CONTEXT if (context == NONE) { // We require a BEGIN line to start a context if (s.find("BEGIN ") == string::npos) { stringstream ss; ss << "Unexpected line in YODA format parsing when BEGIN expected: '" << s << "' on line " << nline; throw ReadError(ss.str()); } // Remove leading #s from the BEGIN line if necessary while (s.find("#") == 0) s = Utils::trim(s.substr(1)); // Split into parts vector parts; istringstream iss(s); string tmp; while (iss >> tmp) parts.push_back(tmp); // Extract context from BEGIN type if (parts.size() < 2 || parts[0] != "BEGIN") { stringstream ss; ss << "Unexpected BEGIN line structure when BEGIN expected: '" << s << "' on line " << nline; throw ReadError(ss.str()); } // Second part is the context name const string ctxstr = parts[1]; // Get block path if possible const string path = (parts.size() >= 3) ? parts[2] : ""; // Set the new context and create a new AO to populate /// @todo Use the block format version for (occasional, careful) format evolution if (Utils::startswith(ctxstr, "YODA_COUNTER")) { context = COUNTER; cncurr = new Counter(path); aocurr = cncurr; } else if (Utils::startswith(ctxstr, "YODA_SCATTER1D")) { context = SCATTER1D; s1curr = new Scatter1D(path); aocurr = s1curr; } else if (Utils::startswith(ctxstr, "YODA_SCATTER2D")) { context = SCATTER2D; s2curr = new Scatter2D(path); aocurr = s2curr; } else if (Utils::startswith(ctxstr, "YODA_SCATTER3D")) { context = SCATTER3D; s3curr = new Scatter3D(path); aocurr = s3curr; } else if (Utils::startswith(ctxstr, "YODA_HISTO1D")) { context = HISTO1D; h1curr = new Histo1D(path); aocurr = h1curr; } else if (Utils::startswith(ctxstr, "YODA_HISTO2D")) { context = HISTO2D; h2curr = new Histo2D(path); aocurr = h2curr; } else if (Utils::startswith(ctxstr, "YODA_PROFILE1D")) { context = PROFILE1D; p1curr = new Profile1D(path); aocurr = p1curr; } else if (Utils::startswith(ctxstr, "YODA_PROFILE2D")) { context = PROFILE2D; p2curr = new Profile2D(path); aocurr = p2curr; } // cout << aocurr->path() << " " << nline << " " << context << endl; // Get block format version if possible (assume version=1 if none found) const size_t vpos = ctxstr.find_last_of("V"); fmt = vpos != string::npos ? ctxstr.substr(vpos+1) : "1"; // cout << fmt << endl; // From version 2 onwards, use the in_anns state from BEGIN until --- if (fmt != "1") in_anns = true; } else { //< not a BEGIN line // Throw error if a BEGIN line is found if (s.find("BEGIN ") != string::npos) ///< @todo require pos = 0 from fmt=V2 throw ReadError("Unexpected BEGIN line in YODA format parsing before ending current BEGIN..END block"); // FINISHING THE CURRENT CONTEXT // Clear/reset context and register AO /// @todo Throw error if mismatch between BEGIN (context) and END types if (s.find("END ") != string::npos) { ///< @todo require pos = 0 from fmt=V2 switch (context) { case COUNTER: break; case HISTO1D: h1curr->addBins(h1binscurr); h1binscurr.clear(); break; case HISTO2D: h2curr->addBins(h2binscurr); h2binscurr.clear(); break; case PROFILE1D: p1curr->addBins(p1binscurr); p1binscurr.clear(); break; case PROFILE2D: p2curr->addBins(p2binscurr); p2binscurr.clear(); break; case SCATTER1D: + for (auto &p : pt1scurr) { p.setParentAO(s1curr); } s1curr->addPoints(pt1scurr); pt1scurr.clear(); break; case SCATTER2D: for (auto &p : pt2scurr) { p.setParentAO(s2curr); } s2curr->addPoints(pt2scurr); pt2scurr.clear(); break; case SCATTER3D: + for (auto &p : pt3scurr) { p.setParentAO(s3curr); } s3curr->addPoints(pt3scurr); pt3scurr.clear(); break; case NONE: break; } // Set all annotations try { YAML::Node anns = YAML::Load(annscurr); // for (YAML::const_iterator it = anns.begin(); it != anns.end(); ++it) { for (const auto& it : anns) { const string key = it.first.as(); // const string val = it.second.as(); YAML::Emitter em; em << YAML::Flow << it.second; //< use single-line formatting, for lists & maps const string val = em.c_str(); // if (!(key.find("ErrorBreakdown") != string::npos)) aocurr->setAnnotation(key, val); } } catch (...) { /// @todo Is there a case for just giving up on these annotations, printing the error msg, and keep going? As an option? const string err = "Problem during annotation parsing of YAML block:\n'''\n" + annscurr + "\n'''"; // cerr << err << endl; throw ReadError(err); } annscurr.clear(); //variationscurr.clear(); in_anns = false; // Put this AO in the completed stack aos.push_back(aocurr); // Clear all current-object pointers aocurr = nullptr; cncurr = nullptr; h1curr = nullptr; h2curr = nullptr; p1curr = nullptr; p2curr = nullptr; s1curr = nullptr; s2curr = nullptr; s3curr = nullptr; context = NONE; continue; } // ANNOTATIONS PARSING if (fmt == "1") { // First convert to one-key-per-line YAML syntax const size_t ieq = s.find("="); if (ieq != string::npos) s.replace(ieq, 1, ": "); // Special-case treatment for syntax clashes const size_t icost = s.find(": *"); if (icost != string::npos) { s.replace(icost, 1, ": '*"); s += "'"; } // Store reformatted annotation const size_t ico = s.find(":"); if (ico != string::npos) { annscurr += (annscurr.empty() ? "" : "\n") + s; continue; } } else if (in_anns) { if (s == "---") { in_anns = false; } else { annscurr += (annscurr.empty() ? "" : "\n") + s; // In order to handle multi-error points in scatters, we need to know which variations are stored, if any // can't wait until we process the annotations at the end, since need to know when filling points. // This is a little inelegant though... //if (s.find("ErrorBreakdown") != string::npos) { // errorBreakdown = YAML::Load(s)["ErrorBreakdown"]; //for (const auto& it : errorBreakdown) { // const string val0 = it.first.as(); //for (const auto& it2 : it.second) { // const string val = it2.as(); //} // } // } } continue; } // DATA PARSING aiss.reset(s); // double sumw(0), sumw2(0), sumwx(0), sumwx2(0), sumwy(0), sumwy2(0), sumwz(0), sumwz2(0), sumwxy(0), sumwxz(0), sumwyz(0), n(0); switch (context) { case COUNTER: { double sumw(0), sumw2(0), n(0); aiss >> sumw >> sumw2 >> n; cncurr->setDbn(Dbn0D(n, sumw, sumw2)); } break; case HISTO1D: { string xoflow1, xoflow2; double xmin(0), xmax(0); double sumw(0), sumw2(0), sumwx(0), sumwx2(0), n(0); /// @todo Improve/factor this "bin" string-or-float parsing... esp for mixed case of 2D overflows /// @todo When outflows are treated as "infinity bins" and don't require a distinct type, string replace under/over -> -+inf if (s.find("Total") != string::npos || s.find("Underflow") != string::npos || s.find("Overflow") != string::npos) { aiss >> xoflow1 >> xoflow2; } else { aiss >> xmin >> xmax; } // The rest is the same for overflows and in-range bins aiss >> sumw >> sumw2 >> sumwx >> sumwx2 >> n; const Dbn1D dbn(n, sumw, sumw2, sumwx, sumwx2); if (xoflow1 == "Total") h1curr->setTotalDbn(dbn); else if (xoflow1 == "Underflow") h1curr->setUnderflow(dbn); else if (xoflow1 == "Overflow") h1curr->setOverflow(dbn); // else h1curr->addBin(HistoBin1D(std::make_pair(xmin,xmax), dbn)); else h1binscurr.push_back(HistoBin1D(std::make_pair(xmin,xmax), dbn)); } break; case HISTO2D: { string xoflow1, xoflow2, yoflow1, yoflow2; double xmin(0), xmax(0), ymin(0), ymax(0); double sumw(0), sumw2(0), sumwx(0), sumwx2(0), sumwy(0), sumwy2(0), sumwxy(0), n(0); /// @todo Improve/factor this "bin" string-or-float parsing... esp for mixed case of 2D overflows /// @todo When outflows are treated as "infinity bins" and don't require a distinct type, string replace under/over -> -+inf if (s.find("Total") != string::npos) { aiss >> xoflow1 >> xoflow2; // >> yoflow1 >> yoflow2; } else if (s.find("Underflow") != string::npos || s.find("Overflow") != string::npos) { throw ReadError("2D histogram overflow syntax is not yet defined / handled"); } else { aiss >> xmin >> xmax >> ymin >> ymax; } // The rest is the same for overflows and in-range bins aiss >> sumw >> sumw2 >> sumwx >> sumwx2 >> sumwy >> sumwy2 >> sumwxy >> n; const Dbn2D dbn(n, sumw, sumw2, sumwx, sumwx2, sumwy, sumwy2, sumwxy); if (xoflow1 == "Total") h2curr->setTotalDbn(dbn); // else if (xoflow1 == "Underflow") p1curr->setUnderflow(dbn); // else if (xoflow1 == "Overflow") p1curr->setOverflow(dbn); else { assert(xoflow1.empty()); // h2curr->addBin(HistoBin2D(std::make_pair(xmin,xmax), std::make_pair(ymin,ymax), dbn)); h2binscurr.push_back(HistoBin2D(std::make_pair(xmin,xmax), std::make_pair(ymin,ymax), dbn)); } } break; case PROFILE1D: { string xoflow1, xoflow2; double xmin(0), xmax(0); double sumw(0), sumw2(0), sumwx(0), sumwx2(0), sumwy(0), sumwy2(0), n(0); /// @todo Improve/factor this "bin" string-or-float parsing... esp for mixed case of 2D overflows /// @todo When outflows are treated as "infinity bins" and don't require a distinct type, string replace under/over -> -+inf if (s.find("Total") != string::npos || s.find("Underflow") != string::npos || s.find("Overflow") != string::npos) { aiss >> xoflow1 >> xoflow2; } else { aiss >> xmin >> xmax; } // The rest is the same for overflows and in-range bins aiss >> sumw >> sumw2 >> sumwx >> sumwx2 >> sumwy >> sumwy2 >> n; const double DUMMYWXY = 0; const Dbn2D dbn(n, sumw, sumw2, sumwx, sumwx2, sumwy, sumwy2, DUMMYWXY); if (xoflow1 == "Total") p1curr->setTotalDbn(dbn); else if (xoflow1 == "Underflow") p1curr->setUnderflow(dbn); else if (xoflow1 == "Overflow") p1curr->setOverflow(dbn); // else p1curr->addBin(ProfileBin1D(std::make_pair(xmin,xmax), dbn)); else p1binscurr.push_back(ProfileBin1D(std::make_pair(xmin,xmax), dbn)); } break; case PROFILE2D: { string xoflow1, xoflow2, yoflow1, yoflow2; double xmin(0), xmax(0), ymin(0), ymax(0); double sumw(0), sumw2(0), sumwx(0), sumwx2(0), sumwy(0), sumwy2(0), sumwz(0), sumwz2(0), sumwxy(0), sumwxz(0), sumwyz(0), n(0); /// @todo Improve/factor this "bin" string-or-float parsing... esp for mixed case of 2D overflows /// @todo When outflows are treated as "infinity bins" and don't require a distinct type, string replace under/over -> -+inf if (s.find("Total") != string::npos) { aiss >> xoflow1 >> xoflow2; // >> yoflow1 >> yoflow2; } else if (s.find("Underflow") != string::npos || s.find("Overflow") != string::npos) { throw ReadError("2D profile overflow syntax is not yet defined / handled"); } else { aiss >> xmin >> xmax >> ymin >> ymax; } // The rest is the same for overflows and in-range bins aiss >> sumw >> sumw2 >> sumwx >> sumwx2 >> sumwy >> sumwy2 >> sumwz >> sumwz2 >> sumwxy >> sumwxz >> sumwyz >> n; const Dbn3D dbn(n, sumw, sumw2, sumwx, sumwx2, sumwy, sumwy2, sumwz, sumwz2, sumwxy, sumwxz, sumwyz); if (xoflow1 == "Total") p2curr->setTotalDbn(dbn); // else if (xoflow1 == "Underflow") p2curr->setUnderflow(dbn); // else if (xoflow1 == "Overflow") p2curr->setOverflow(dbn); else { assert(xoflow1.empty()); // p2curr->addBin(ProfileBin2D(std::make_pair(xmin,xmax), std::make_pair(ymin,ymax), dbn)); p2binscurr.push_back(ProfileBin2D(std::make_pair(xmin,xmax), std::make_pair(ymin,ymax), dbn)); } } break; case SCATTER1D: { double x(0), exm(0), exp(0); aiss >> x >> exm >> exp; // set nominal point Point1D thispoint=Point1D(x, exm, exp); // check if we stored variations of this point //if (variationscurr.size()>0){ // // for each variation, store the alt errors. // // start at 1 since we have already filled nominal ! // for (unsigned int ivar=1; ivar> exm >> exp; // thispoint.setXErrs(exm,exp,thisvariation); // } //} pt1scurr.push_back(thispoint); } break; case SCATTER2D: { double x(0), y(0), exm(0), exp(0), eym(0), eyp(0); aiss >> x >> exm >> exp >> y >> eym >> eyp; // set nominal point Point2D thispoint=Point2D(x, y, exm, exp, eym, eyp); // check if we stored variations of this point // for each variation, store the alt errors. // start at 1 since we have already filled nominal ! pt2scurr.push_back(thispoint); } break; case SCATTER3D: { double x(0), y(0), z(0), exm(0), exp(0), eym(0), eyp(0), ezm(0), ezp(0); aiss >> x >> exm >> exp >> y >> eym >> eyp >> z >> ezm >> ezp; // set nominal point Point3D thispoint=Point3D(x, y, z, exm, exp, eym, eyp, ezm, ezp); pt3scurr.push_back(thispoint); } break; default: throw ReadError("Unknown context in YODA format parsing: how did this happen?"); } // cout << "AO CONTENT " << nline << endl; // cout << " " << xmin << " " << xmax << " " << ymin << " " << ymax << " / '" << xoflow1 << "' '" << xoflow2 << "' '" << yoflow1 << "' '" << yoflow2 << "'" << endl; // cout << " " << sumw << " " << sumw2 << " " << sumwx << " " << sumwx2 << " " << sumwy << " " << sumwy2 << " " << sumwz << " " << sumwz2 << " " << sumwxy << " " << sumwxz << " " << sumwyz << " " << n << endl; // cout << " " << x << " " << y << " " << z << " " << exm << " " << exp << " " << eym << " " << eyp << " " << ezm << " " << ezp << endl; } } } }