diff --git a/include/YODA/AnalysisObject.h b/include/YODA/AnalysisObject.h --- a/include/YODA/AnalysisObject.h +++ b/include/YODA/AnalysisObject.h @@ -1,301 +1,308 @@ // -*- 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 ; + } //@} ///@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,97 +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 + 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,356 +1,360 @@ // -*- 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) { } /// Copy assignment Point1D& operator = (const Point1D& p) { _x = p._x; _ex = p._ex; 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 (!_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 (!_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="") { _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,561 +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); - return (_ey.at(source).first + _ey.at(source).second)/2.0; + 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.); + 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 { - return _ey; - } + const std::map< std::string, std::pair> & errMap() const; + + // Get 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,679 +1,683 @@ // -*- 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) { } /// 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; 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 (!_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 (!_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 (!_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 (!_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 (!_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 (!_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="") { _ez[source] = ez; } /// Get value minus negative z-error double zMin( std::string source="") const { 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 (!_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,330 +1,334 @@ // -*- 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)); } /// Insert a new point, defined as the x value and symmetric errors void addPoint(double x, double ex) { _points.insert(Point1D(x, ex)); } /// 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)); } /// 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)); } /// 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,424 +1,440 @@ // -*- 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) { - _points.insert(Point2D(x, y)); + Point2D thisPoint= Point2D(x, y); + _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) { - _points.insert(Point2D(x, y, ex, ey)); + Point2D thisPoint= Point2D(x, y, ex, ey); + _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) { - _points.insert(Point2D(x, y, ex, ey)); + Point2D thisPoint= Point2D(x, y, ex, ey); + _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) { - _points.insert(Point2D(x, y, exminus, explus, eyminus, eyplus)); + Point2D thisPoint=Point2D(x, y, exminus, explus, eyminus, eyplus); + _points.insert(thisPoint); + thisPoint.setParentAO(this); } /// Insert a collection of new points void addPoints(const Points& pts) { - for (const Point2D& pt : pts) addPoint(pt); + 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,436 +1,441 @@ // -*- 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)); } /// 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)); } /// 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)); } /// 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)); } /// 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/pyext/yoda/declarations.pxd b/pyext/yoda/declarations.pxd --- a/pyext/yoda/declarations.pxd +++ b/pyext/yoda/declarations.pxd @@ -1,1433 +1,1437 @@ from libcpp.map cimport map from libcpp.pair cimport pair from libcpp.vector cimport vector from libcpp cimport bool from libcpp.string cimport string from cython.operator cimport dereference as deref cdef extern from "YODA/Config/YodaConfig.h" namespace "YODA": string version() # Import the error handling C++ routine cdef extern from "errors.hh": # Have a look in errors.cpp for implementation specifics void yodaerr "translate_yoda_error" () ctypedef map[string, string] Annotations ctypedef double (*dbl_dbl_fptr) (double) ctypedef map[string, pair[double,double]] errMap # Math utils {{{ cdef extern from "YODA/Utils/MathUtils.h" namespace "YODA": # bool isZero(double a, double tolerance) # bool fuzzyEquals(double a, double b, double tolerance) # bool fuzzyGtrEquals(double a, double b, double tolerance) # bool fuzzyLessEquals(double a, double b, double tolerance) vector[double] linspace(size_t nbins, double start, double end) vector[double] logspace(size_t nbins, double start, double end) int index_between(double&, vector[double]& binedges) double mean(vector[int]& sample) double covariance(vector[int]& sample1, vector[int]& sample2) double correlation(vector[int]& sample1, vector[int]& sample2) # }}} # Dbn0D {{{ cdef extern from "YODA/Dbn0D.h" namespace "YODA": cdef cppclass Dbn0D: Dbn0D () Dbn0D (Dbn0D) void fill(double weight, double fraction) void reset() void scaleW(double) # Raw distribution running sums unsigned long numEntries() except +yodaerr double effNumEntries() except +yodaerr double sumW() except +yodaerr double sumW2() except +yodaerr double errW() except +yodaerr double relErrW() except +yodaerr Dbn0D operator+ (Dbn0D) Dbn0D operator- (Dbn0D) # TODO: += and -= operators #}}} Dbn0D # Dbn1D {{{ cdef extern from "YODA/Dbn1D.h" namespace "YODA": cdef cppclass Dbn1D: Dbn1D () Dbn1D (Dbn1D) void fill(double val, double weight, double fraction) void reset() void scaleW(double) void scaleX(double) double errW() except +yodaerr double relErrW() except +yodaerr double xMean() except +yodaerr double xVariance() except +yodaerr double xStdDev() except +yodaerr double xStdErr() except +yodaerr double xRMS() except +yodaerr # Raw distribution running sums unsigned long numEntries() except +yodaerr double effNumEntries() except +yodaerr double sumW() except +yodaerr double sumW2() except +yodaerr double sumWX() except +yodaerr double sumWX2() except +yodaerr Dbn1D operator+ (Dbn1D) Dbn1D operator- (Dbn1D) # TODO: += and -= operators #}}} Dbn1D # Dbn2D {{{ cdef extern from "YODA/Dbn2D.h" namespace "YODA": cdef cppclass Dbn2D: Dbn2D () Dbn2D (Dbn2D) void fill(double x, double y, double weight, double fraction) except +yodaerr void reset() except +yodaerr void scaleW(double) except +yodaerr void scaleX(double) except +yodaerr void scaleY(double) except +yodaerr void scaleXY(double, double) except +yodaerr double errW() except +yodaerr double relErrW() except +yodaerr double xMean() except +yodaerr double xVariance() except +yodaerr double xStdDev() except +yodaerr double xStdErr() except +yodaerr double xRMS() except +yodaerr double yMean() except +yodaerr double yVariance() except +yodaerr double yStdDev() except +yodaerr double yStdErr() except +yodaerr double yRMS() except +yodaerr # Raw distribution running sums unsigned long numEntries() except +yodaerr double effNumEntries() except +yodaerr double sumW() except +yodaerr double sumW2() except +yodaerr double sumWX() except +yodaerr double sumWX2() except +yodaerr double sumWY() except +yodaerr double sumWY2() except +yodaerr double sumWXY() except +yodaerr # Operators void flipXY() except +yodaerr Dbn1D transformX() except +yodaerr Dbn1D transformY() except +yodaerr Dbn2D operator + (Dbn2D) Dbn2D operator - (Dbn2D) # TODO: += and -= operators #}}} Dbn2D # Dbn3D {{{ cdef extern from "YODA/Dbn3D.h" namespace "YODA": cdef cppclass Dbn3D: Dbn3D () Dbn3D (Dbn3D) void fill(double x, double y, double z, double weight, double fraction) void reset() void scaleW(double) void scaleX(double) void scaleY(double) void scaleZ(double) # void scaleXY(double, double) # void scaleYZ(double, double) # void scaleXZ(double, double) void scaleXYZ(double, double, double) double errW() except +yodaerr double relErrW() except +yodaerr double xMean() double xVariance() double xStdDev() double xStdErr() double xRMS() double yMean() double yVariance() double yStdDev() double yStdErr() double yRMS() double zMean() double zVariance() double zStdDev() double zStdErr() double zRMS() # Raw distribution running sums unsigned long numEntries() double effNumEntries() double sumW() double sumW2() double sumWX() double sumWX2() double sumWY() double sumWY2() double sumWZ() double sumWZ2() double sumWXY() double sumWXZ() double sumWYZ() double sumWXYZ() # Operators void flipXY() void flipXZ() void flipYZ() Dbn1D transformX() Dbn1D transformY() Dbn1D transformZ() Dbn3D operator + (Dbn3D) Dbn3D operator - (Dbn3D) # TODO: += and -= operators #}}} Dbn3D # Point {{{ cdef extern from "YODA/Point.h" namespace "YODA": cdef cppclass Point: int dim() except +yodaerr double val(size_t i) except +yodaerr void setVal(size_t i, double val) except +yodaerr pair[double,double] errs(size_t i) except +yodaerr pair[double,double] errs(size_t i, string source) except +yodaerr double errMinus(size_t i) except +yodaerr double errMinus(size_t i, string source) except +yodaerr void setErrMinus(size_t i, double eminus) except +yodaerr void setErrMinus(size_t i, double eminus, string source) except +yodaerr double errPlus(size_t i) except +yodaerr double errPlus(size_t i, string source) except +yodaerr void setErrPlus(size_t i, double eplus) except +yodaerr void setErrPlus(size_t i, double eplus, string source) except +yodaerr double errAvg(size_t i) except +yodaerr double errAvg(size_t i, string source) except +yodaerr void setErr(size_t i, double e) except +yodaerr void setErr(size_t i, double e, string source) except +yodaerr # void setErrs(size_t i, double e) except +yodaerr # void setErrs(size_t i, double eminus, double eplus) except +yodaerr void setErrs(size_t i, pair[double,double]& e) except +yodaerr void setErrs(size_t i, pair[double,double]& e, string source) except +yodaerr # void set(size_t i, double val, double e) except +yodaerr # void set(size_t i, double val, double eminus, double eplus) except +yodaerr void set(size_t i, double val, pair[double,double]& e) except +yodaerr void set(size_t i, double val, pair[double,double]& e, string source) except +yodaerr errMap errMap() except +yodaerr #}}} Point # Point1D {{{ cdef extern from "YODA/Point1D.h" namespace "YODA": cdef cppclass Point1D(Point): Point1D () except +yodaerr Point1D (Point1D p) except +yodaerr Point1D (double x, double exminus, double explus) except +yodaerr Point1D (double x, double exminus, double explus, string source) except +yodaerr double x() except +yodaerr void setX(double x) except +yodaerr pair[double,double] xErrs() except +yodaerr pair[double,double] xErrs(string source) except +yodaerr void setXErrs(pair[double, double]&) except +yodaerr void setXErrs(pair[double, double]&, string source) except +yodaerr double xErrAvg() except +yodaerr double xErrAvg(string source) except +yodaerr double xMin() except +yodaerr double xMin(string source) except +yodaerr double xMax() except +yodaerr double xMax(string source) except +yodaerr void scaleX(double) except +yodaerr bool operator == (Point1D) except +yodaerr bool operator != (Point1D b) except +yodaerr bool operator < (Point1D b) except +yodaerr bool operator <= (Point1D b) except +yodaerr bool operator > (Point1D b) except +yodaerr bool operator >= (Point1D b) except +yodaerr # }}} Point1D # Point2D {{{ cdef extern from "YODA/Point2D.h" namespace "YODA": cdef cppclass Point2D(Point): Point2D () except +yodaerr Point2D (Point2D p) except +yodaerr Point2D (double x, double y, double exminus, double explus, double eyminus, double eyplus) except +yodaerr Point2D (double x, double y, double exminus, double explus, double eyminus, double eyplus, string source) except +yodaerr double x() except +yodaerr double y() except +yodaerr void setX(double x) except +yodaerr void setY(double y) except +yodaerr pair[double,double] xy() except +yodaerr void setXY(pair[double,double]&) except +yodaerr pair[double,double] xErrs() except +yodaerr pair[double,double] yErrs() except +yodaerr pair[double,double] yErrs(string source) except +yodaerr void setXErrs(pair[double, double]&) except +yodaerr void setYErrs(pair[double, double]&) except +yodaerr void setYErrs(pair[double, double]&, string source) except +yodaerr double xErrAvg() except +yodaerr double yErrAvg() except +yodaerr double yErrAvg(string source) except +yodaerr double xMin() except +yodaerr double xMax() except +yodaerr double yMin() except +yodaerr double yMin(string source) except +yodaerr double yMax() except +yodaerr double yMax(string source) except +yodaerr void scaleX(double) except +yodaerr void scaleY(double) except +yodaerr void scaleXY(double, double) except +yodaerr #void scale(double, double) except +yodaerr bool operator == (Point2D) except +yodaerr bool operator != (Point2D b) except +yodaerr bool operator < (Point2D b) except +yodaerr bool operator <= (Point2D b) except +yodaerr bool operator > (Point2D b) except +yodaerr bool operator >= (Point2D b) except +yodaerr # }}} Point2D # Point3D {{{ cdef extern from "YODA/Point3D.h" namespace "YODA": cdef cppclass Point3D(Point): Point3D () except +yodaerr Point3D (Point3D& p) except +yodaerr Point3D (double x, double y, double z, double exminus, double explus, double eyminus, double eyplus, double ezminus, double ezplus) except +yodaerr Point3D (double x, double y, double z, double exminus, double explus, double eyminus, double eyplus, double ezminus, double ezplus, string source) except +yodaerr double x() except +yodaerr double y() except +yodaerr double z() except +yodaerr void setX(double x) except +yodaerr void setY(double y) except +yodaerr void setZ(double z) except +yodaerr pair[double,double] xErrs() except +yodaerr pair[double,double] yErrs() except +yodaerr pair[double,double] zErrs() except +yodaerr pair[double,double] zErrs(string source) except +yodaerr void setXErrs(pair[double, double]&) except +yodaerr void setYErrs(pair[double, double]&) except +yodaerr void setZErrs(pair[double, double]&) except +yodaerr void setZErrs(pair[double, double]&, string source) except +yodaerr double xErrAvg() double yErrAvg() double zErrAvg() double zErrAvg(string source) double xMin() except +yodaerr double xMax() except +yodaerr double yMin() except +yodaerr double yMax() except +yodaerr double zMin() except +yodaerr double zMin(string source) except +yodaerr double zMax() except +yodaerr double zMax(string source) except +yodaerr void scaleX(double) except +yodaerr void scaleY(double) except +yodaerr void scaleZ(double) except +yodaerr void scaleXYZ(double, double, double) except +yodaerr #void scale(double, double, double) except +yodaerr bool operator == (Point3D b) bool operator != (Point3D b) bool operator < (Point3D b) bool operator <= (Point3D b) bool operator > (Point3D b) bool operator >= (Point3D b) #}}} Point3D # Bin {{{ cdef extern from "YODA/Bin.h" namespace "YODA": cdef cppclass Bin: int dim() except +yodaerr unsigned long numEntries() except +yodaerr double effNumEntries() except +yodaerr double sumW() except +yodaerr double sumW2() except +yodaerr # }}} Bin #Bin1D {{{ cdef extern from "YODA/Bin1D.h" namespace "YODA": cdef cppclass Bin1D[DBN](Bin): Bin1D(pair[double, double] edges) except +yodaerr Bin1D(pair[double, double] edges, DBN dbn) except +yodaerr Bin1D(Bin1D) except +yodaerr # THIS IS A CYTHON LIMITATION... DO NOT CALL THIS Bin1D() # (DO NOT CALL THIS DO NOT CALL THIS) ### ################################################# #We're fine as long as we don't try to instantiate these from Python # void scaleW(double scale) except +yodaerr # void scaleX(double scale) except +yodaerr void reset() except +yodaerr pair[double, double] edges() except +yodaerr double xMin() except +yodaerr double xMax() except +yodaerr double xMid() except +yodaerr double xWidth() except +yodaerr double xFocus() except +yodaerr # x statistics double xMean() except +yodaerr double xVariance() except +yodaerr double xStdDev() except +yodaerr double xStdErr() except +yodaerr double xRMS() except +yodaerr # raw statistics double sumWX() except +yodaerr double sumWX2() except +yodaerr void merge (Bin1D&) except +yodaerr Bin1D operator + (Bin1D&) Bin1D operator - (Bin1D&) ctypedef Bin1D[Dbn1D] Bin1D_Dbn1D ctypedef Bin1D[Dbn2D] Bin1D_Dbn2D ctypedef Bin1D[Dbn3D] Bin1D_Dbn3D #}}} Bin1D # Bin2D {{{ cdef extern from "YODA/Bin2D.h" namespace "YODA": cdef cppclass Bin2D[DBN](Bin): Bin2D(pair[double, double] xedges, pair[double, double] yedges) except+ Bin2D(Bin2D bin) except +yodaerr # CYTHON HACK DO NOT CALL THIS IT DOES NOT EXIST Bin2D() # (DO NOT CALL DO NOT CALL) ################################################ # void scaleW(double scale) except +yodaerr # void scaleXY(double, double) except +yodaerr void reset() except +yodaerr pair[double, double] xEdges() except +yodaerr pair[double, double] yEdges() except +yodaerr double xMin() except +yodaerr double yMin() except +yodaerr double xMax() except +yodaerr double yMax() except +yodaerr double xMid() except +yodaerr double yMid() except +yodaerr double xWidth() except +yodaerr double yWidth() except +yodaerr double area() except +yodaerr double xFocus() except +yodaerr double yFocus() except +yodaerr pair[double, double] xyFocus() except +yodaerr pair[double, double] xyMid() except +yodaerr # x statistics double xMean() except +yodaerr double xVariance() except +yodaerr double xStdDev() except +yodaerr double xStdErr() except +yodaerr double xRMS() except +yodaerr double yMean() except +yodaerr double yVariance() except +yodaerr double yStdDev() except +yodaerr double yStdErr() except +yodaerr double yRMS() except +yodaerr # Raw statistics double sumWX() except +yodaerr double sumWY() except +yodaerr double sumWXY() except +yodaerr double sumWX2() except +yodaerr double sumWY2() except +yodaerr #void merge(Bin2D) except +yodaerr Bin2D operator + (Bin2D) Bin2D operator - (Bin2D) int adjacentTo(Bin2D) except +yodaerr ctypedef Bin2D[Dbn2D] Bin2D_Dbn2D ctypedef Bin2D[Dbn3D] Bin2D_Dbn3D # }}} Bin2D # HistoBin1D {{{ cdef extern from "YODA/HistoBin1D.h" namespace "YODA": cdef cppclass HistoBin1D(Bin1D_Dbn1D): HistoBin1D(double lowedge, double highedge) except +yodaerr HistoBin1D(HistoBin1D) except +yodaerr # void fill(double x, double weight, double fraction) except +yodaerr # void fillBin(double weight, double fraction) except +yodaerr double area() except +yodaerr double height() except +yodaerr double areaErr() except +yodaerr double heightErr() except +yodaerr double relErr() except +yodaerr HistoBin1D operator+(HistoBin1D) HistoBin1D operator-(HistoBin1D) #}}} HistoBin1D cdef extern from "merge.hh": void HistoBin1D_iadd_HistoBin1D "cython_iadd" (HistoBin1D*, HistoBin1D*) void HistoBin1D_isub_HistoBin1D "cython_isub" (HistoBin1D*, HistoBin1D*) # void HistoBin1D_imul_dbl "cython_imul_dbl" (HistoBin1D*, double) # void HistoBin1D_idiv_dbl "cython_idiv_dbl" (HistoBin1D*, double) HistoBin1D* HistoBin1D_add_HistoBin1D "cython_add" (HistoBin1D*, HistoBin1D*) HistoBin1D* HistoBin1D_sub_HistoBin1D "cython_sub" (HistoBin1D*, HistoBin1D*) HistoBin1D* HistoBin1D_div_HistoBin1D "cython_div" (HistoBin1D*, HistoBin1D*) # HistoBin2D {{{ cdef extern from "YODA/HistoBin2D.h" namespace "YODA": cdef cppclass HistoBin2D(Bin2D_Dbn2D): HistoBin2D(double xmin, double xmax, double ymin, double ymax) except +yodaerr HistoBin2D(HistoBin2D) except +yodaerr # void fill(double x, double y, double weight, double fraction) except +yodaerr # void fillBin(double weight, double fraction) except +yodaerr void reset() # Accessors double volume() except +yodaerr double volumeErr() except +yodaerr double height() except +yodaerr double heightErr() except +yodaerr double relErr() except +yodaerr HistoBin2D operator+(HistoBin2D) HistoBin2D operator-(HistoBin2D) #Bin2D_Dbn2D merge(HistoBin2D b) #}}} HistoBin2D # ProfileBin1D {{{ cdef extern from "YODA/ProfileBin1D.h" namespace "YODA": cdef cppclass ProfileBin1D(Bin1D_Dbn2D): ProfileBin1D(ProfileBin1D) except +yodaerr ProfileBin1D(double, double) except +yodaerr #void fill(double x, double y, double weight, double fraction) except +yodaerr #void fillBin(double y, double weight, double fraction) except +yodaerr void reset() except +yodaerr double mean() except +yodaerr double stdDev() except +yodaerr double variance() except +yodaerr double stdErr() except +yodaerr double rms() except +yodaerr double sumWY() except +yodaerr double sumWY2() except +yodaerr ProfileBin1D operator + (ProfileBin1D) ProfileBin1D operator - (ProfileBin1D) # void scaleY(double) except +yodaerr # }}} ProfileBin1D cdef extern from "merge.hh": void ProfileBin1D_iadd_ProfileBin1D "cython_iadd" (ProfileBin1D*, ProfileBin1D*) void ProfileBin1D_isub_ProfileBin1D "cython_isub" (ProfileBin1D*, ProfileBin1D*) # void ProfileBin1D_imul_dbl "cython_imul_dbl" (ProfileBin1D*, double) # void ProfileBin1D_idiv_dbl "cython_idiv_dbl" (ProfileBin1D*, double) ProfileBin1D* ProfileBin1D_add_ProfileBin1D "cython_add" (ProfileBin1D*, ProfileBin1D*) ProfileBin1D* ProfileBin1D_sub_ProfileBin1D "cython_sub" (ProfileBin1D*, ProfileBin1D*) ProfileBin1D* ProfileBin1D_div_ProfileBin1D "cython_div" (ProfileBin1D*, ProfileBin1D*) # ProfileBin2D {{{ cdef extern from "YODA/ProfileBin2D.h" namespace "YODA": cdef cppclass ProfileBin2D(Bin2D_Dbn3D): ProfileBin2D (ProfileBin2D h) except +yodaerr ProfileBin2D (double, double, double, double) except +yodaerr # void fill(double x, double y, double z, double weight, double fraction) except +yodaerr # void fillBin(double z, double weight, double fraction) except +yodaerr double mean() except +yodaerr double stdDev() except +yodaerr double variance() except +yodaerr double stdErr() except +yodaerr double rms() except +yodaerr double sumWZ() except +yodaerr double sumWZ2() except +yodaerr ProfileBin2D operator + (ProfileBin2D) ProfileBin2D operator - (ProfileBin2D) # void scaleZ(double) except +yodaerr # }}} ProfileBin2D # AnalysisObject {{{ cdef extern from "YODA/AnalysisObject.h" namespace "YODA": cdef cppclass AnalysisObject: # Constructors AnalysisObject(string type, string path, string title) except +yodaerr AnalysisObject(string type, string path, AnalysisObject ao, string title) except +yodaerr AnalysisObject() #AnalysisObject* newclone() except +yodaerr ## String used in automatic type determination string type() except +yodaerr ## Data object fill- or plot-space dimension int dim() except +yodaerr ## Annotations vector[string] annotations() except +yodaerr bool hasAnnotation(string key) except +yodaerr string annotation(string key) except +yodaerr string annotation(string key, string default) except +yodaerr void setAnnotation(string, string) except +yodaerr void rmAnnotation(string name) except +yodaerr void clearAnnotations() except +yodaerr ## Standard annotations string title() except +yodaerr void setTitle(string title) except +yodaerr string path() except +yodaerr void setPath(string title) except +yodaerr string name() except +yodaerr # }}} AnalysisObject cdef extern from "YODA/Utils/sortedvector.h" namespace "YODA::Utils": cdef cppclass sortedvector[T](vector): sortedvector(vector[T]) except +yodaerr void insert(T) except +yodaerr # TODO: forward declarations for bin-copying constructors # Counter {{{ cdef extern from "YODA/Counter.h" namespace "YODA": cdef cppclass Counter(AnalysisObject): Counter() except +yodaerr Counter(string path, string title) except +yodaerr #Counter(Dbn0D dbn, string path, string title) except +yodaerr Counter(Counter c, string path) Counter clone() except +yodaerr Counter* newclone() except +yodaerr void reset() except +yodaerr void fill(double weight, double fraction) except +yodaerr unsigned long numEntries() except +yodaerr double effNumEntries() except +yodaerr double sumW() except +yodaerr double sumW2() except +yodaerr double val() except +yodaerr double err() except +yodaerr double relErr() except +yodaerr void scaleW(double) except +yodaerr # operator += (Counter) # operator -= (Counter) Scatter1D Counter_div_Counter "divide" (const Counter&, const Counter&) except +yodaerr Scatter1D Counter_eff_Counter "efficiency" (const Counter&, const Counter&) except +yodaerr cdef extern from "merge.hh": void Counter_iadd_Counter "cython_iadd" (Counter*, Counter*) void Counter_isub_Counter "cython_isub" (Counter*, Counter*) # void Counter_imul_dbl "cython_imul_dbl" (Counter*, double) # void Counter_idiv_dbl "cython_idiv_dbl" (Counter*, double) Counter* Counter_add_Counter "cython_add" (Counter*, Counter*) Counter* Counter_sub_Counter "cython_sub" (Counter*, Counter*) #Counter* Counter_div_Counter "cython_div" (Counter*, Counter*) cdef extern from "YODA/Scatter1D.h" namespace "YODA": Scatter1D mkScatter_Counter "YODA::mkScatter" (const Counter&) except +yodaerr #}}} Counter # Scatter1D {{{ cdef extern from "YODA/Scatter1D.h" namespace "YODA": cdef cppclass Scatter1D(AnalysisObject): Scatter1D() except +yodaerr Scatter1D(string path, string title) except +yodaerr Scatter1D(sortedvector[Point1D], string path, string title) except +yodaerr Scatter1D(vector[double], vector[double], vector[pair[double, double]], vector[pair[double, double]]) except +yodaerr Scatter1D(Scatter1D p, string path) Scatter1D clone() except +yodaerr Scatter1D* newclone() except +yodaerr void reset() except +yodaerr size_t numPoints() except +yodaerr # TODO: have to ignore exception handling on ref-returning methods until Cython bug is fixed vector[Point1D]& points() #except +yodaerr Point1D& point(size_t index) #except +yodaerr void addPoint(const Point1D&) #except +yodaerr void addPoint(double) #except +yodaerr void addPoint(double, const pair[double, double]&) #except +yodaerr void addPoints(const sortedvector[Point1D]&) #except +yodaerr void combineWith(const Scatter1D&) #except +yodaerr void combineWith(const vector[Scatter1D]&) #except +yodaerr void scaleX(double) except +yodaerr + void parseVariations() except +yodaerr + vector[string] variations() except +yodaerr vector[vector[double]] covarianceMatrix(bool) except +yodaerr void Scatter1D_transformX "YODA::transformX" (Scatter1D&, dbl_dbl_fptr) #}}} Scatter1D # cdef extern from "merge.hh": # Scatter2D* Scatter2D_add_Scatter2D "cython_add" (Scatter2D*, Scatter2D*) # Scatter2D* Scatter2D_sub_Scatter2D "cython_sub" (Scatter2D*, Scatter2D*) cdef extern from "YODA/Scatter1D.h" namespace "YODA": Scatter1D mkScatter_Scatter1D "YODA::mkScatter" (const Scatter1D&) except +yodaerr # Scatter2D {{{ cdef extern from "YODA/Scatter2D.h" namespace "YODA": cdef cppclass Scatter2D(AnalysisObject): Scatter2D() except +yodaerr Scatter2D(string path, string title) except +yodaerr Scatter2D(sortedvector[Point2D], string path, string title) except +yodaerr Scatter2D(vector[double], vector[double], vector[pair[double, double]], vector[pair[double, double]]) except +yodaerr Scatter2D(Scatter2D p, string path) Scatter2D clone() except +yodaerr Scatter2D* newclone() except +yodaerr void reset() except +yodaerr size_t numPoints() except +yodaerr # TODO: have to ignore exception handling on ref-returning methods until Cython bug is fixed vector[Point2D]& points() #except +yodaerr Point2D& point(size_t index) #except +yodaerr void addPoint(const Point2D&) #except +yodaerr void addPoint(double, double) #except +yodaerr void addPoint(double, double, const pair[double, double]&, const pair[double, double]&) #except +yodaerr void addPoints(const sortedvector[Point2D]&) #except +yodaerr void combineWith(const Scatter2D&) #except +yodaerr void combineWith(const vector[Scatter2D]&) #except +yodaerr void scaleX(double) except +yodaerr void scaleY(double) except +yodaerr void scaleXY(double, double) except +yodaerr #void scale(double, double) except +yodaerr + void parseVariations() except +yodaerr vector[string] variations() except +yodaerr vector[vector[double]] covarianceMatrix(bool) except +yodaerr void Scatter2D_transformX "YODA::transformX" (Scatter2D&, dbl_dbl_fptr) void Scatter2D_transformY "YODA::transformY" (Scatter2D&, dbl_dbl_fptr) #}}} Scatter2D # cdef extern from "merge.hh": # Scatter2D* Scatter2D_add_Scatter2D "cython_add" (Scatter2D*, Scatter2D*) # Scatter2D* Scatter2D_sub_Scatter2D "cython_sub" (Scatter2D*, Scatter2D*) cdef extern from "YODA/Scatter2D.h" namespace "YODA": Scatter2D mkScatter_Scatter2D "YODA::mkScatter" (const Scatter2D&) except +yodaerr # Scatter3D {{{ cdef extern from "YODA/Scatter3D.h" namespace "YODA": cdef cppclass Scatter3D(AnalysisObject): Scatter3D() except +yodaerr Scatter3D(string path, string title) except +yodaerr Scatter3D(sortedvector[Point3D], string path, string title) except +yodaerr Scatter3D(vector[double], vector[double], vector[pair[double, double]], vector[pair[double, double]], vector[pair[double, double]]) except +yodaerr Scatter3D(Scatter3D p, string path) Scatter3D clone() except +yodaerr Scatter3D* newclone() except +yodaerr void reset() except +yodaerr size_t numPoints() except +yodaerr # TODO: have to ignore exception handling on ref-returning methods until Cython bug is fixed sortedvector[Point3D]& points() #except +yodaerr Point3D& point(size_t index) #except +yodaerr void addPoint(const Point3D&) #except +yodaerr void addPoint(double, double, double) #except +yodaerr void addPoint(double, double, double, const pair[double, double]&, const pair[double, double]&, const pair[double, double]&) #except +yodaerr void addPoints(const sortedvector[Point3D]&) #except +yodaerr void combineWith(const Scatter3D&) #except +yodaerr void combineWith(const vector[Scatter3D]&) #except +yodaerr void scaleX(double) except +yodaerr void scaleY(double) except +yodaerr void scaleZ(double) except +yodaerr void scaleXYZ(double, double, double) except +yodaerr #void scale(double, double, double) except +yodaerr + void parseVariations() except +yodaerr vector[string] variations() except +yodaerr vector[vector[double]] covarianceMatrix() except +yodaerr vector[vector[double]] covarianceMatrix(bool) except +yodaerr void Scatter3D_transformX "YODA::transformX" (Scatter3D&, dbl_dbl_fptr) void Scatter3D_transformY "YODA::transformY" (Scatter3D&, dbl_dbl_fptr) void Scatter3D_transformZ "YODA::transformZ" (Scatter3D&, dbl_dbl_fptr) #}}} Scatter3D # cdef extern from "merge.hh": # Scatter3D* Scatter3D_add_Scatter3D "cython_add" (Scatter3D*, Scatter3D*) # Scatter3D* Scatter3D_sub_Scatter3D "cython_sub" (Scatter3D*, Scatter3D*) cdef extern from "YODA/Scatter3D.h" namespace "YODA": Scatter3D mkScatter_Scatter3D "YODA::mkScatter" (const Scatter3D&) except +yodaerr # Histo1D#{{{ cdef extern from "YODA/Histo1D.h" namespace "YODA": cdef cppclass Histo1D(AnalysisObject): Histo1D() except +yodaerr Histo1D(string path, string title) except +yodaerr Histo1D(size_t nbins, double lower, double upper, string path, string title) except +yodaerr Histo1D(vector[double] binedges, string path, string title) except +yodaerr Histo1D(vector[Bin] bins, string path, string title) except +yodaerr Histo1D(Histo1D h, string path) except +yodaerr #Histo1D(Profile1D p, string path) #Histo1D(Scatter2D p, string path) Histo1D clone() except +yodaerr Histo1D* newclone() except +yodaerr void reset() except +yodaerr void fill(double x, double weight, double fraction) except +yodaerr void fillBin(size_t i, double weight, double fraction) except +yodaerr void scaleW(double s) except +yodaerr void normalize(double normto, bool includeoverflows) except +yodaerr void mergeBins(size_t, size_t) except +yodaerr void rebinBy(unsigned int n, size_t begin, size_t end) except +yodaerr void rebinTo(vector[double] edges) except +yodaerr void addBin(double, double) except +yodaerr void addBins(vector[double] edges) except +yodaerr void eraseBin(size_t index) except +yodaerr vector[double] xEdges() except +yodaerr double xMin() except +yodaerr double xMax() except +yodaerr size_t numBins() except +yodaerr vector[HistoBin1D]& bins() int binIndexAt(double x) except +yodaerr const HistoBin1D& bin(size_t ix) const HistoBin1D& binAt(double x) except +yodaerr # TODO: Some Cython mapping problem? Dbn1D& totalDbn() Dbn1D& underflow() Dbn1D& overflow() # Whole histo data double integral(bool) double integralTo(int, bool) double integralRange(int, int) unsigned long numEntries(bool) double effNumEntries(bool) double sumW(bool) double sumW2(bool) double xMean(bool) double xVariance(bool) double xStdDev(bool) double xStdErr(bool) double xRMS(bool) # operator == (Histo1D) # operator != (Histo1D) operator + (Histo1D) operator - (Histo1D) operator / (Histo1D) Scatter2D Histo1D_toIntegral "toIntegralHisto" (const Histo1D& h, bool includeunderflow) except +yodaerr Scatter2D Histo1D_toIntegralEff "toIntegralEfficiencyHisto" (const Histo1D& h, bool includeunderflow, bool includeoverflow) except +yodaerr Scatter2D Histo1D_div_Histo1D "divide" (const Histo1D&, const Histo1D&) except +yodaerr Scatter2D Histo1D_eff_Histo1D "efficiency" (const Histo1D&, const Histo1D&) except +yodaerr cdef extern from "merge.hh": void Histo1D_iadd_Histo1D "cython_iadd" (Histo1D*, Histo1D*) void Histo1D_isub_Histo1D "cython_isub" (Histo1D*, Histo1D*) # void Histo1D_imul_dbl "cython_imul_dbl" (Histo1D*, double) # void Histo1D_idiv_dbl "cython_idiv_dbl" (Histo1D*, double) Histo1D* Histo1D_add_Histo1D "cython_add" (Histo1D*, Histo1D*) Histo1D* Histo1D_sub_Histo1D "cython_sub" (Histo1D*, Histo1D*) Histo1D* Histo1D_div_Histo1D "cython_div" (Histo1D*, Histo1D*) cdef extern from "YODA/Scatter2D.h" namespace "YODA": Scatter2D mkScatter_Histo1D "YODA::mkScatter" (const Histo1D&, bool) except +yodaerr #}}} Histo1D # Histo2D {{{ cdef extern from "YODA/Histo2D.h" namespace "YODA": cdef cppclass Histo2D(AnalysisObject): Histo2D() except +yodaerr Histo2D(string path, string title) except +yodaerr Histo2D(size_t nBinsX, double lowerX, double upperX, size_t nBinsY, double lowerY, double upperY, string path, string title) except +yodaerr Histo2D(vector[double] xedges, vector[double] yedges, string path, string title) except +yodaerr Histo2D(Histo2D, string path) #Histo2D(Profile1D p, string path) #Histo2D(Scatter2D p, string path) Histo2D clone() except +yodaerr Histo2D* newclone() except +yodaerr # TODO: add missing functions and enable refs + exceptions when Cython allows void reset() except +yodaerr void fill(double x, double y, double weight, double fraction) except +yodaerr void fillBin(size_t i, double weight, double fraction) except +yodaerr void normalize(double normto, bool includeoverflows) except +yodaerr void scaleW(double scalefactor) except +yodaerr void scaleXY(double, double) # void mergeBins(size_t, size_t) except +yodaerr # void rebin(unsigned int n) except +yodaerr size_t numBins() except +yodaerr size_t numBinsX() except +yodaerr size_t numBinsY() except +yodaerr vector[HistoBin2D]& bins() #except +yodaerr int binIndexAt(double x, double y) except +yodaerr const HistoBin2D& bin(size_t ix) #except +yodaerr const HistoBin2D& binAt(double x, double y) #except +yodaerr void addBin(const pair[double, double]&, const pair[double, double]&) void addBins(const vector[HistoBin2D]&) void addBin(double, double) except +yodaerr void addBins(const vector[double]& edges) except +yodaerr # void eraseBin(size_t index) except +yodaerr vector[double] xEdges() except +yodaerr vector[double] yEdges() except +yodaerr double xMin() except +yodaerr double xMax() except +yodaerr double yMin() except +yodaerr double yMax() except +yodaerr # Dbn2D& outflow(int, int) #except +yodaerr # Whole histo data Dbn2D& totalDbn() #except +yodaerr double integral(bool) unsigned long numEntries(bool) double effNumEntries(bool) double sumW(bool) double sumW2(bool) double xMean(bool) double yMean(bool) double xVariance(bool) double yVariance(bool) double xStdDev(bool) double yStdDev(bool) double xStdErr(bool) double yStdErr(bool) double xRMS(bool) double yRMS(bool) # operator == (Histo2D) # operator != (Histo2D) operator + (Histo2D) operator - (Histo2D) operator / (Histo2D) Scatter3D Histo2D_div_Histo2D "divide" (const Histo2D&, const Histo2D&) except +yodaerr Scatter3D Histo2D_eff_Histo2D "efficiency" (const Histo2D&, const Histo2D&) except +yodaerr cdef extern from "merge.hh": void Histo2D_iadd_Histo2D "cython_iadd" (Histo2D*, Histo2D*) void Histo2D_isub_Histo2D "cython_isub" (Histo2D*, Histo2D*) # void Histo2D_imul_dbl "cython_imul_dbl" (Histo2D*, double) # void Histo2D_idiv_dbl "cython_idiv_dbl" (Histo2D*, double) Histo2D* Histo2D_add_Histo2D "cython_add" (Histo2D*, Histo2D*) Histo2D* Histo2D_sub_Histo2D "cython_sub" (Histo2D*, Histo2D*) Histo2D* Histo2D_div_Histo2D "cython_div" (Histo2D*, Histo2D*) cdef extern from "YODA/Scatter3D.h" namespace "YODA": Scatter3D mkScatter_Histo2D "YODA::mkScatter" (const Histo2D&, bool) except +yodaerr # Histo2D }}} # Profile1D {{{ cdef extern from "YODA/Profile1D.h" namespace "YODA": cdef cppclass Profile1D(AnalysisObject): Profile1D() except +yodaerr Profile1D(string path, string title) except +yodaerr Profile1D(size_t nxbins, double xlower, double xupper, string path, string title) except +yodaerr Profile1D(vector[double] xbinedges, string path, string title) except +yodaerr Profile1D(Profile1D p, string path) except +yodaerr Profile1D(Scatter2D s, string path) except +yodaerr #Profile1D(Histo1D p, string path) Profile1D clone() except +yodaerr Profile1D* newclone() except +yodaerr void reset() except +yodaerr void fill(double x, double y, double weight, double fraction) except +yodaerr void fillBin(size_t i, double y, double weight, double fraction) except +yodaerr void scaleW(double s) except +yodaerr void scaleY(double s) except +yodaerr void mergeBins(size_t, size_t) except +yodaerr void rebinBy(unsigned int n, size_t begin, size_t end) except +yodaerr void rebinTo(vector[double] edges) except +yodaerr void addBin(double, double) except +yodaerr void addBins(vector[double] edges) except +yodaerr # TODO: void eraseBin(size_t index) except +yodaerr vector[double] xEdges() except +yodaerr double xMin() except +yodaerr double xMax() except +yodaerr size_t numBins() except +yodaerr vector[ProfileBin1D] bins() #except +yodaerr int binIndexAt(double x) except +yodaerr const ProfileBin1D& bin(size_t ix) #except +yodaerr const ProfileBin1D& binAt(double x) #except +yodaerr # The trick here is to treat these not as references. # I suppose when you think about it, it makes sense Dbn2D& totalDbn() Dbn2D& underflow() Dbn2D& overflow() unsigned long numEntries(bool) double effNumEntries(bool) double sumW(bool) double sumW2(bool) double xMean(bool) double xVariance(bool) double xStdDev(bool) double xStdErr(bool) double xRMS(bool) operator + (Profile1D) operator - (Profile1D) operator / (Profile1D) Scatter2D Profile1D_div_Profile1D "divide" (const Profile1D&, const Profile1D&) except +yodaerr cdef extern from "merge.hh": void Profile1D_iadd_Profile1D "cython_iadd" (Profile1D*, Profile1D*) void Profile1D_isub_Profile1D "cython_isub" (Profile1D*, Profile1D*) # void Profile1D_imul_dbl "cython_imul_dbl" (Profile1D*, double) # void Profile1D_idiv_dbl "cython_idiv_dbl" (Profile1D*, double) Profile1D* Profile1D_add_Profile1D "cython_add" (Profile1D*, Profile1D*) Profile1D* Profile1D_sub_Profile1D "cython_sub" (Profile1D*, Profile1D*) Profile1D* Profile1D_div_Profile1D "cython_div" (Profile1D*, Profile1D*) cdef extern from "YODA/Scatter2D.h" namespace "YODA": Scatter2D mkScatter_Profile1D "YODA::mkScatter" (const Profile1D&, bool, bool) except +yodaerr #}}} Profile1D # Profile2D {{{ cdef extern from "YODA/Profile2D.h" namespace "YODA": cdef cppclass Profile2D(AnalysisObject): Profile2D() except +yodaerr Profile2D(string path, string title) except +yodaerr Profile2D(size_t nbinsX, double lowerX, double upperX, size_t nbinsY, double lowerY, double upperY, string path, string title) except +yodaerr Profile2D(vector[double] xedges, vector[double] yedges, string path, string title) except +yodaerr Profile2D(Profile2D p, string path) except +yodaerr #Profile2D(Scatter3D s, string path) except +yodaerr #Profile2D(Histo2D p, string path) Profile2D clone() except +yodaerr Profile2D* newclone() except +yodaerr # TODO: add missing functions and enable refs + exceptions when Cython allows void reset() except +yodaerr void fill(double x, double y, double z, double weight, double fraction) except +yodaerr void fillBin(size_t i, double z, double weight, double fraction) except +yodaerr void scaleW(double s) except +yodaerr void scaleXY(double, double) # void mergeBins(size_t, size_t) except +yodaerr # void rebin(unsigned int n) except +yodaerr size_t numBins() except +yodaerr size_t numBinsX() except +yodaerr size_t numBinsY() except +yodaerr vector[ProfileBin2D]& bins() #except +yodaerr int binIndexAt(double x, double y) except +yodaerr const ProfileBin2D& bin(size_t ix) #except +yodaerr const ProfileBin2D& binAt(double x, y) #except +yodaerr void addBin(const pair[double, double]&, const pair[double, double]&) except +yodaerr void addBins(const vector[double]&, const vector[double]&) except +yodaerr # void eraseBin(size_t index) except +yodaerr vector[double] xEdges() except +yodaerr vector[double] yEdges() except +yodaerr double xMin() except +yodaerr double xMax() except +yodaerr double yMin() except +yodaerr double yMax() except +yodaerr # Dbn3D& outflow(int, int) #except +yodaerr # Whole histo data Dbn3D& totalDbn() #except +yodaerr unsigned long numEntries(bool) double effNumEntries(bool) double sumW(bool) double sumW2(bool) double xMean(bool) double yMean(bool) double xVariance(bool) double yVariance(bool) double xStdDev(bool) double yStdDev(bool) double xStdErr(bool) double yStdErr(bool) double xRMS(bool) double yRMS(bool) operator + (Profile2D) operator - (Profile2D) operator / (Profile2D) Scatter3D Profile2D_div_Profile2D "divide" (const Profile2D&, const Profile2D&) except +yodaerr cdef extern from "merge.hh": void Profile2D_iadd_Profile2D "cython_iadd" (Profile2D*, Profile2D*) void Profile2D_isub_Profile2D "cython_isub" (Profile2D*, Profile2D*) # void Profile2D_imul_dbl "cython_imul_dbl" (Profile2D*, double) # void Profile2D_idiv_dbl "cython_idiv_dbl" (Profile2D*, double) Profile2D* Profile2D_add_Profile2D "cython_add" (Profile2D*, Profile2D*) Profile2D* Profile2D_sub_Profile2D "cython_sub" (Profile2D*, Profile2D*) Profile2D* Profile2D_div_Profile2D "cython_div" (Profile2D*, Profile2D*) cdef extern from "YODA/Scatter3D.h" namespace "YODA": Scatter3D mkScatter_Profile2D "YODA::mkScatter" (const Profile2D&, bool, bool) except +yodaerr #}}} Profile2D # Streams {{{ cdef extern from "" namespace "std": cdef cppclass istringstream: istringstream() string& str(string&) cdef cppclass ostringstream: ostringstream() string& str() cdef extern from "YODA/IO.h" namespace "YODA": void IO_read_from_file "YODA::read" (string&, vector[AnalysisObject*]&) except +yodaerr cdef extern from "YODA/Reader.h" namespace "YODA": cdef cppclass Reader: void read(istringstream&, vector[AnalysisObject*]&) except +yodaerr void read_from_file "YODA::Reader::read" (string&, vector[AnalysisObject*]&) except +yodaerr cdef extern from "YODA/ReaderYODA.h" namespace "YODA": Reader& ReaderYODA_create "YODA::ReaderYODA::create" () cdef extern from "YODA/ReaderFLAT.h" namespace "YODA": Reader& ReaderFLAT_create "YODA::ReaderFLAT::create" () cdef extern from "YODA/ReaderAIDA.h" namespace "YODA": Reader& ReaderAIDA_create "YODA::ReaderAIDA::create" () cdef extern from "YODA/Reader.h" namespace "YODA": Reader& Reader_create "YODA::mkReader" (string& filename) cdef extern from "YODA/IO.h" namespace "YODA": void IO_write_to_file "YODA::write" (string&, vector[AnalysisObject*]&) except +yodaerr cdef extern from "YODA/Writer.h" namespace "YODA": cdef cppclass Writer: void write(ostringstream&, vector[AnalysisObject*]&) except +yodaerr void write_to_file "YODA::Writer::write" (string&, vector[AnalysisObject*]&) except +yodaerr cdef extern from "YODA/WriterYODA.h" namespace "YODA": Writer& WriterYODA_create "YODA::WriterYODA::create" () cdef extern from "YODA/WriterFLAT.h" namespace "YODA": Writer& WriterFLAT_create "YODA::WriterFLAT::create" () cdef extern from "YODA/WriterAIDA.h" namespace "YODA": Writer& WriterAIDA_create "YODA::WriterAIDA::create" () cdef extern from "YODA/Reader.h" namespace "YODA": Writer& Writer_create "YODA::mkWriter" (string& filename) # Streams }}} # Axis1D {{{ cdef extern from "YODA/Axis1D.h" namespace "YODA": cdef cppclass Axis1D[BIN1D, DBN]: Axis1D() except +yodaerr Axis1D(vector[double] binedges) except +yodaerr Axis1D(size_t, double, double) except +yodaerr Axis1D(vector[BIN1D] bins) except +yodaerr void addBin(double, double) except +yodaerr size_t numBins() except +yodaerr vector[BIN1D]& bins() double xMin() except +yodaerr double xMax() except +yodaerr vector[double] xEdges() except +yodaerr long getBinIndex(double) void reset() DBN& totalDbn() DBN& underflow() DBN& overflow() void eraseBin(size_t index) except +yodaerr void mergeBins(size_t, size_t) except +yodaerr # Axis1D }}} # Axis2D {{{ cdef extern from "YODA/Axis2D.h" namespace "YODA": cdef cppclass Axis2D[BIN2D, DBN]: Axis2D() except +yodaerr Axis2D(vector[double], vector[double]) except +yodaerr Axis2D(size_t, pair[double, double], size_t, pair[double, double]) except +yodaerr Axis2D(vector[BIN2D] bins) except +yodaerr void addBin(pair[double, double], pair[double, double]) except +yodaerr size_t numBins() except +yodaerr vector[BIN2D]& bins() double xMin() except +yodaerr double xMax() except +yodaerr double yMin() except +yodaerr double yMax() except +yodaerr long getBinIndex(double, double) void reset() DBN& totalDbn() # TODO: reinstate DBN& outflow(int, int) void eraseBin(size_t index) except +yodaerr void mergeBins(size_t, size_t) except +yodaerr # Axis2D }}} diff --git a/pyext/yoda/include/Point1D.pyx b/pyext/yoda/include/Point1D.pyx --- a/pyext/yoda/include/Point1D.pyx +++ b/pyext/yoda/include/Point1D.pyx @@ -1,89 +1,121 @@ cimport util cdef class Point1D(Point): """ A 1D point with errors, used by the Scatter1D class. """ cdef c.Point1D* p1ptr(self) except NULL: return self.ptr() def __init__(self, x=0, xerrs=0, source=""): if source==None: source="" cutil.set_owned_ptr(self, new c.Point1D()) self.setX(x) self.setXErrs(xerrs, source) def copy(self): return cutil.new_owned_cls(Point1D, new c.Point1D(deref(self.p1ptr()))) # TODO: add clone() as mapping to (not yet existing) C++ newclone()? def setXErrs(self, val, source): if source==None: source="" self.p1ptr().setXErrs(util.read_symmetric(val)) # property x: # """x coordinate""" # def __get__(self): # return self.p1ptr().x() # def __set__(self, x): # self.p1ptr().setX(x) # property xErrs: # """The x errors""" # def __get__(self): # return util.read_error_pair(self.p1ptr().xErrs()) # def __set__(self, val): # self.p1ptr().setXErrs(util.read_symmetric(val)) def x(self): """The x value""" return self.p1ptr().x() def setX(self, x): """Set the x value""" self.p1ptr().setX(x) def xErrs(self): """The x errors""" return util.read_error_pair(self.p1ptr().xErrs()) - def setXErrs(self, val): - """Set the x errors""" - def __set__(self, val): - self.p1ptr().setXErrs(util.read_symmetric(val)) - + + def xErrsFromSource(self, source): + """The y errors""" + if isinstance(source, str): + source = source.encode('utf-8') + return util.read_error_pair(self.p1ptr().xErrs(source)) + + def setXErrs(self, *es): + """(int, float) -> None + (int, [float, float]) -> None + (int, float, float) -> None + Set asymmetric errors on axis i""" + source = None + es = list(es) + if type(es[-1]) is str: + source = es[-1] + es = es[:-1] + else: + pass + errs = es + if source is None: + source = "" + if len(errs) == 1: + if not hasattr(errs[0], "__iter__"): + self.setErr(1,errs[0], source) + return + errs = errs[0] + # assert len(errs) == 2: + if isinstance(source, str): + source = source.encode('utf-8') + self.pptr().setErrs(1, tuple(errs), source) + + def setYErrs(self, val, source): + if source is None: + source = "" + self.p1ptr().setXErrs(util.read_symmetric(val)) + #@property def xMin(self): """The minimum x position, i.e. lowest error""" return self.p1ptr().xMin() #@property def xMax(self): """The maximum x position, i.e. highest error""" return self.p1ptr().xMax() def xErrAvg(self): return self.p1ptr().xErrAvg() def scaleX(self, a): """(float) -> None Scale the x values and errors by factor a.""" self.p1ptr().scaleX(a) def __repr__(self): return '' % self.x() def __richcmp__(Point1D self, Point1D other, int op): if op == 0: return deref(self.p1ptr()) < deref(other.p1ptr()) elif op == 1: return deref(self.p1ptr()) <= deref(other.p1ptr()) elif op == 2: return deref(self.p1ptr()) == deref(other.p1ptr()) elif op == 3: return deref(self.p1ptr()) != deref(other.p1ptr()) elif op == 4: return deref(self.p1ptr()) > deref(other.p1ptr()) elif op == 5: return deref(self.p1ptr()) >= deref(other.p1ptr()) diff --git a/pyext/yoda/include/Point2D.pyx b/pyext/yoda/include/Point2D.pyx --- a/pyext/yoda/include/Point2D.pyx +++ b/pyext/yoda/include/Point2D.pyx @@ -1,168 +1,175 @@ cimport util cdef class Point2D(Point): """ A 2D point with errors, used by the Scatter2D class. """ cdef c.Point2D* p2ptr(self) except NULL: return self.ptr() def __init__(self, x=0, y=0, xerrs=0, yerrs=0, source=""): if source==None: source="" cutil.set_owned_ptr(self, new c.Point2D()) self.setX(x) self.setY(y) self.setXErrs(xerrs) self.setYErrs(yerrs, source) def copy(self): return cutil.new_owned_cls(Point2D, new c.Point2D(deref(self.p2ptr()))) # TODO: add clone() as mapping to (not yet existing) C++ newclone()? def x(self): """The x value""" return self.p2ptr().x() def setX(self, x): """Set the x value""" self.p2ptr().setX(x) def xErrs(self): """The x errors""" return util.read_error_pair(self.p2ptr().xErrs()) def setXErrs(self, val): """Set the x errors""" self.p2ptr().setXErrs(util.read_symmetric(val)) def xMin(self): """The minimum x position, i.e. lowest error""" return self.p2ptr().xMin() def xMax(self): """The maximum x position, i.e. highest error""" return self.p2ptr().xMax() def xErrAvg(self): return self.p2ptr().xErrAvg() def y(self): """The y value""" return self.p2ptr().y() def setY(self, y): """Set the y value""" self.p2ptr().setY(y) def yErrs(self): """The y errors""" return util.read_error_pair(self.p2ptr().yErrs()) + + def yErrsFromSource(self, source): + """The y errors""" + if isinstance(source, str): + source = source.encode('utf-8') + return util.read_error_pair(self.p2ptr().yErrs(source)) # def setYErrs(self, val): # """Set the y errors""" # self.p2ptr().setYErrs(util.read_symmetric(val)) def setYErrs(self, *es): """(int, float) -> None (int, [float, float]) -> None (int, float, float) -> None Set asymmetric errors on axis i""" source = None es = list(es) if type(es[-1]) is str: source = es[-1] es = es[:-1] else: pass errs = es if source is None: source = "" if len(errs) == 1: if not hasattr(errs[0], "__iter__"): self.setErr(2,errs[0], source) return errs = errs[0] # assert len(errs) == 2: if isinstance(source, str): source = source.encode('utf-8') self.pptr().setErrs(2, tuple(errs), source) + def setYErrs(self, val, source): if source is None: source = "" self.p2ptr().setYErrs(util.read_symmetric(val)) def yMin(self): """The minimum y position, i.e. lowest error""" return self.p2ptr().yMin() def yMax(self): """The maximum y position, i.e. highest error""" return self.p2ptr().yMax() def yErrAvg(self): return self.p2ptr().yErrAvg() # property xy: # """x and y coordinates as a tuple""" # def __get__(self): # return util.XY(self.x, self.y) # def __set__(self, val): # self.x, self.y = val # # TODO: How does this fit into the multi-error API? Still useful, but just reports first errs... how to get _all_ +- err pairs? # # LC: I think it's Ok to leave this like this, for most users the nominal is what they want anyway, # # and for those who want multi-errs, they can set using a method eg setErrs(dim,(ed,eu),source) and access using errs(dim,(ed,eu),source) # property yErrs: # """The y errors""" # def __get__(self): # return util.read_error_pair(self.p2ptr().yErrs()) # def __set__(self, val): # self.p2ptr().setYErrs(util.read_symmetric(val)) def scaleX(self, a): """(float) -> None Scale the x values and errors by factor a.""" self.p2ptr().scaleX(a) def scaleY(self, a): """(float) -> None Scale the y values and errors by factor a.""" self.p2ptr().scaleY(a) def scaleXY(self, x=1.0, y=1.0): """ (float=1, float=1) -> None Scale the point coordinates by the given factors. """ self.p2ptr().scaleXY(x, y) # TODO: remove def scale(self, x=1.0, y=1.0): """ (float=1, float=1) -> None DEPRECATED! Use scaleXY Scale the point coordinates by the given factors. """ self.p2ptr().scaleXY(x, y) def __repr__(self): return '' % (self.x(), self.y()) def __richcmp__(Point2D self, Point2D other, int op): if op == 0: return deref(self.p2ptr()) < deref(other.p2ptr()) elif op == 1: return deref(self.p2ptr()) <= deref(other.p2ptr()) elif op == 2: return deref(self.p2ptr()) == deref(other.p2ptr()) elif op == 3: return deref(self.p2ptr()) != deref(other.p2ptr()) elif op == 4: return deref(self.p2ptr()) > deref(other.p2ptr()) elif op == 5: return deref(self.p2ptr()) >= deref(other.p2ptr()) diff --git a/pyext/yoda/include/Point3D.pyx b/pyext/yoda/include/Point3D.pyx --- a/pyext/yoda/include/Point3D.pyx +++ b/pyext/yoda/include/Point3D.pyx @@ -1,222 +1,228 @@ cimport util cdef class Point3D(Point): """ A 3D point with errors, used by the Scatter3D class. """ cdef c.Point3D* p3ptr(self) except NULL: return self.ptr() def __init__(self, x=0, y=0, z=0, xerrs=0, yerrs=0, zerrs=0, source=""): if source==None: source="" cutil.set_owned_ptr(self, new c.Point3D()) # TODO: need shortcuts self.setX(x) self.setY(y) self.setZ(z) self.setXErrs(xerrs) self.setYErrs(yerrs) self.setZErrs(zerrs, source) def copy(self): return cutil.new_owned_cls(Point3D, new c.Point3D(deref(self.p3ptr()))) def x(self): """The x value""" return self.p3ptr().x() def setX(self, x): """Set the x value""" self.p3ptr().setX(x) def xErrs(self): """The x errors""" return util.read_error_pair(self.p3ptr().xErrs()) def setXErrs(self, val): """Set the x errors""" self.p3ptr().setXErrs(util.read_symmetric(val)) def xMin(self): """The minimum x position, i.e. lowest error""" return self.p3ptr().xMin() def xMax(self): """The maximum x position, i.e. highest error""" return self.p3ptr().xMax() def xErrAvg(self): return self.p3ptr().xErrAvg() def y(self): """The y value""" return self.p3ptr().y() def setY(self, y): """Set the y value""" self.p3ptr().setY(y) def yErrs(self): """The y errors""" return util.read_error_pair(self.p3ptr().yErrs()) def setYErrs(self, val): """Set the y errors""" self.p3ptr().setYErrs(util.read_symmetric(val)) def yMin(self): """The minimum y position, i.e. lowest error""" return self.p3ptr().yMin() def yMax(self): """The maximum y position, i.e. highest error""" return self.p3ptr().yMax() def yErrAvg(self): return self.p3ptr().yErrAvg() def z(self): """The z value""" return self.p3ptr().z() def setZ(self, z): """Set the z value""" self.p3ptr().setZ(z) def zErrs(self): """The z errors""" return util.read_error_pair(self.p3ptr().zErrs()) + + def zErrsFromSource(self, source): + """The z errors""" + if isinstance(source, str): + source = source.encode('utf-8') + return util.read_error_pair(self.p3ptr().zErrs(source)) # def setZErrs(self, val): # """Set the z errors""" # self.p3ptr().setZErrs(util.read_symmetric(val)) def setZErrs(self, *es): """(int, float) -> None (int, [float, float]) -> None (int, float, float) -> None Set asymmetric errors on axis i""" source = None es = list(es) if type(es[-1]) is str: source = es[-1] es = es[:-1] else: pass errs = es if source is None: source = "" if len(errs) == 1: if not hasattr(errs[0], "__iter__"): self.setErr(2,errs[0], source) return errs = errs[0] # assert len(errs) == 2: if isinstance(source, str): source = source.encode('utf-8') self.pptr().setErrs(2, tuple(errs), source) def setZErrs(self, val, source): if source is None: source = "" self.p3ptr().setZErrs(util.read_symmetric(val)) def zMin(self): """The minimum z position, i.e. lowest error""" return self.p3ptr().zMin() def zMax(self): """The maximum z position, i.e. highest error""" return self.p3ptr().zMax() def zErrAvg(self): return self.p3ptr().zErrAvg() # property xyz: # def __get__(self): # return util.XYZ(self.x, self.y, self.z) # def __set__(self, val): # self.x, self.y, self.z = val # # TODO: How does this fit into the multi-error API? Still useful, but just reports first errs... how to get _all_ +- err pairs? # # LC: This is fine because preserntly only the highest dimension supports multi-errors # property xErrs: # def __get__(self): # return util.read_error_pair(self.p3ptr().xErrs()) # def __set__(self, val): # self.p3ptr().setXErrs(util.read_symmetric(val)) # # TODO: How does this fit into the multi-error API? Still useful, but just reports first errs... how to get _all_ +- err pairs? # # LC: This is fine because preserntly only the highest dimension supports multi-errors # property yErrs: # def __get__(self): # return util.read_error_pair(self.p3ptr().yErrs()) # def __set__(self, val): # self.p3ptr().setYErrs(util.read_symmetric(val)) # # TODO: How does this fit into the multi-error API? Still useful, but just reports first errs... how to get _all_ +- err pairs? # # LC: I think it's Ok to leave this like this, for most users the nominal is what they want anyway, # # and for those who want multi-errs, they can set using a method eg setErrs(dim,(ed,eu),source) and access using errs(dim,(ed,eu),source) # property zErrs: # def __get__(self): # return util.read_error_pair(self.p3ptr().zErrs()) # def __set__(self, val): # self.p3ptr().setZErrs(util.read_symmetric(val)) def scaleX(self, ax): """ (float) -> None Scale the x point coordinates by the given factor. """ self.p3ptr().scaleX(ax) def scaleY(self, ay): """ (float) -> None Scale the y point coordinates by the given factor. """ self.p3ptr().scaleY(ay) def scaleZ(self, az): """ (float) -> None Scale the z point coordinates by the given factor. """ self.p3ptr().scaleZ(az) def scaleXYZ(self, ax=1.0, ay=1.0, az=1.0): """ (float=1.0, float=1.0, float=1.0) -> None Scale the point coordinates by the given factors. """ self.p3ptr().scaleXYZ(ax, ay, az) # TODO: remove def scaleXYZ(self, ax=1.0, ay=1.0, az=1.0): """ (double=1.0, double=1.0, double=1.0) -> None DEPRECATED: USE scaleXYZ Scale the point coordinates by the given factors. """ self.scaleXYZ(ax, ay, az) # TODO: transformX,Y,Z def __repr__(self): return '' % (self.x(), self.y(), self.z()) def __richcmp__(Point3D self, Point3D other, int op): if op == 0: return deref(self.p3ptr()) < deref(other.p3ptr()) elif op == 1: return deref(self.p3ptr()) <= deref(other.p3ptr()) elif op == 2: return deref(self.p3ptr()) == deref(other.p3ptr()) elif op == 3: return deref(self.p3ptr()) != deref(other.p3ptr()) elif op == 4: return deref(self.p3ptr()) > deref(other.p3ptr()) elif op == 5: return deref(self.p3ptr()) >= deref(other.p3ptr()) diff --git a/pyext/yoda/include/Scatter2D.pyx b/pyext/yoda/include/Scatter2D.pyx --- a/pyext/yoda/include/Scatter2D.pyx +++ b/pyext/yoda/include/Scatter2D.pyx @@ -1,288 +1,294 @@ cimport util cdef class Scatter2D(AnalysisObject): """ 2D scatter plot, i.e. a collection of Point2D objects with positions and errors. Constructor calling idioms: Scatter2D(path="", title="") Create a new empty scatter, with optional path and title. Scatter2D(points, path="", title=""): Create a new empty scatter from an iterable of points, with optional path and title. TODO: more documentation! """ cdef inline c.Scatter2D* s2ptr(self) except NULL: return self.ptr() def __init__(self, *args, **kwargs): util.try_loop([self.__init_2, self.__init_3], *args, **kwargs) def __init_2(self, path="", title=""): path = path.encode('utf-8') title = title.encode('utf-8') cutil.set_owned_ptr(self, new c.Scatter2D(path, title)) def __init_3(self, points, path="", title=""): self.__init_2(path, title) self.addPoints(points) def clone(self): """() -> Scatter2D. Clone this Scatter2D.""" return cutil.new_owned_cls(Scatter2D, self.s2ptr().newclone()) def __repr__(self): return "<%s '%s' %d points>" % (self.__class__.__name__, self.path(), len(self.points())) #@property def numPoints(self): """() -> int Number of points in this scatter.""" return self.s2ptr().numPoints() def __len__(self): return self.numPoints() #@property def points(self): """Access the ordered list of points.""" return [self.point(i) for i in xrange(self.numPoints())] def point(self, size_t i): """Access the i'th point.""" return cutil.new_borrowed_cls(Point2D, &self.s2ptr().point(i), self) def __getitem__(self, py_ix): cdef size_t i = cutil.pythonic_index(py_ix, self.numPoints()) return cutil.new_borrowed_cls(Point2D, &self.s2ptr().point(i), self) def addPoint(self, *args, **kwargs): """Add a new point. Provide either a single yoda.Point2D object, or the four args: x, y, xerrs=0, yerrs=0. """ try: self.__addPoint_point(*args, **kwargs) except TypeError: self.__addPoint_explicit(*args, **kwargs) def __addPoint_explicit(self, x, y, xerrs=0, yerrs=0): self.__addPoint_point(Point2D(x, y, xerrs, yerrs)) def __addPoint_point(self, Point2D p): self.s2ptr().addPoint(p.p2ptr()[0]) def addPoints(self, iterable): """Add several new points.""" for row in iterable: try: self.addPoint(*row) except TypeError: self.addPoint(row) def combineWith(self, others): """Try to add points from other Scatter2Ds into this one.""" cdef Scatter2D other try: # Can we type it as a Scatter2D? other = others except TypeError: # Could be an iterable... for other in others: self.s2ptr().combineWith(deref(other.s2ptr())) else: self.s2ptr().combineWith(deref(other.s2ptr())) def mkScatter(self): """None -> Scatter2D. Make a new Scatter2D. Exists to allow mkScatter calls on any AnalysisObject, even if it already is a scatter.""" cdef c.Scatter2D s2 = c.mkScatter_Scatter2D(deref(self.s2ptr())) return cutil.new_owned_cls(Scatter2D, s2.newclone()) def scaleX(self, a): """(float) -> None Scale the x values and errors of the points in this scatter by factor a.""" self.s2ptr().scaleX(a) def scaleY(self, a): """(float) -> None Scale the y values and errors of the points in this scatter by factor a.""" self.s2ptr().scaleY(a) def scaleXY(self, ax=1.0, ay=1.0): """(float=1, float=1) -> None Scale the values and errors of the points in this scatter by factors ax, ay.""" self.s2ptr().scaleXY(ax, ay) # TODO: remove def scale(self, ax=1.0, ay=1.0): """(float=1, float=1) -> None DEPRECATED: USE scaleXY Scale the values and errors of the points in this scatter by factors ax, ay.""" self.scaleXY(ax, ay) def transformX(self, f): """(fn) -> None Transform the x values and errors of the points in this scatter by function f.""" import ctypes try: callback = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double)(f) except: raise RuntimeError("Callback is not of type (double) -> double") fptr = (ctypes.addressof(callback))[0] c.Scatter2D_transformX(deref(self.s2ptr()), fptr) def transformY(self, f): """(fn) -> None Transform the y values and errors of the points in this scatter by function f.""" import ctypes try: callback = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double)(f) except: raise RuntimeError("Callback is not of type (double) -> double") fptr = (ctypes.addressof(callback))[0] c.Scatter2D_transformY(deref(self.s2ptr()), fptr) + + def parseVariations(self): + """None -> None + Parse the YAML which contains the variations stored in the poins of the Scatter. + Only needs to be done once!""" + return self.s2ptr().parseVariations() def variations(self): """None -> vector[string] Get the list of variations stored in the poins of the Scatter""" return self.s2ptr().variations() def _mknp(self, xs): try: import numpy return numpy.array(xs) except ImportError: return xs def covarianceMatrix(self, *ignoreOffDiagonalTerms): """bool -> vector[vector[float]] Construct the covariance matrix""" return self._mknp(self.s2ptr().covarianceMatrix(ignoreOffDiagonalTerms)) # # TODO: remove? # def __add__(Scatter2D self, Scatter2D other): # return cutil.new_owned_cls(Scatter2D, c.Scatter2D_add_Scatter2D(self.s2ptr(), other.s2ptr())) # # TODO: remove? # def __sub__(Scatter2D self, Scatter2D other): # return cutil.new_owned_cls(Scatter2D, c.Scatter2D_sub_Scatter2D(self.s2ptr(), other.s2ptr())) def hasValidErrorBreakdown(self): """ Check if the AO's error breakdown is not empty and has no bins withh 0 uncertainty """ counter = -1 for p in self.points(): counter += 1 binErrs = p.errMap() if len(binErrs) < 2: return False binTotal = [0.,0.] for sys, err in binErrs.iteritems(): binTotal[0] = (binTotal[0]**2 + err[0]**2)**0.5 binTotal[1] = (binTotal[1]**2 + err[1]**2)**0.5 if binTotal[0] == 0 and binTotal[1] == 0: return False return True def correlationMatrix(self): """ `covMatrix` numpy matrix Convert a covariance matrix to a correlation matrix (ie normalise entry in i,j by uncertainty of bin i * uncertainty in bin j) """ covMatrix = self.covarianceMatrix() nbins = len(covMatrix) corr = [[0 for i in range(nbins)] for j in range (nbins)] for i in range(nbins): sigma_i = covMatrix[i][i] for j in range(nbins): sigma_j = covMatrix[j][j] corr[i][j] = covMatrix[i][j] / (sigma_i * sigma_j)**0.5 return self._mknp(corr) def xVals(self): return self._mknp([p.x() for p in self.points()]) def xMins(self): """All x low values.""" # TODO: add extra dimensionality for multiple errors? return self._mknp([p.xMin() for p in self.points()]) def xMaxs(self): """All x high values.""" # TODO: add extra dimensionality for multiple errors? return self._mknp([p.xMax() for p in self.points()]) def xErrs(self): """All x error pairs""" # TODO: add extra dimensionality for multiple errors? return self._mknp([p.xErrs() for p in self.points()]) def xErrAvgs(self): """All x average errors""" # TODO: add extra dimensionality for multiple errors? return self._mknp([p.xAvgErr() for p in self.points()]) def xMin(self): """Lowest x value.""" # TODO: add extra dimensionality for multiple errors? return min(self.xMins()) def xMax(self): """Highest x value.""" # TODO: add extra dimensionality for multiple errors? return max(self.xMaxs()) def yVals(self): return self._mknp([p.y() for p in self.points()]) def yMins(self): """All y low values.""" # TODO: add extra dimensionality for multiple errors? return self._mknp([p.yMin() for p in self.points()]) def yMaxs(self): """All y high values.""" # TODO: add extra dimensionality for multiple errors? return self._mknp([p.yMax() for p in self.points()]) def yErrs(self): """All y error pairs""" # TODO: add extra dimensionality for multiple errors? return self._mknp([p.yErrs() for p in self.points()]) def yErrAvgs(self): """All y average errors""" # TODO: add extra dimensionality for multiple errors? return self._mknp([p.yAvgErr() for p in self.points()]) def yMin(self): """Lowest x value.""" # TODO: add extra dimensionality for multiple errors? return min(self.yMins()) def yMax(self): """Highest y value.""" # TODO: add extra dimensionality for multiple errors? return max(self.yMaxs()) ## Convenience alias S2D = Scatter2D diff --git a/src/Makefile.am b/src/Makefile.am --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,30 +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 \ - Dbn1D.cc \ + Dbn0D.cc \ + Dbn1D.cc \ Counter.cc \ Histo1D.cc \ Histo2D.cc \ Profile1D.cc \ Profile2D.cc \ Scatter1D.cc \ Scatter2D.cc \ - Scatter3D.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,548 +1,526 @@ // -*- 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; - YAML::Node errorBreakdown; // 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: s1curr->addPoints(pt1scurr); pt1scurr.clear(); break; case SCATTER2D: + for (auto &p : pt2scurr) { p.setParentAO(s2curr); } s2curr->addPoints(pt2scurr); pt2scurr.clear(); break; case SCATTER3D: 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(); - // - // The Variations annotation is just a placeholder to help collect the right columns - // Don't want to be saving it to the actual AO, since the method variations() - // provides the info that's needed without needing to keep the annotation up to date - // LC: Not using the column representation of multiple errors right now. - //if (!(key.find("Variations") != string::npos)) aocurr->setAnnotation(key, val); - if (!(key.find("ErrorBreakdown") != string::npos)) aocurr->setAnnotation(key, val); + // 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"]; + //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); - int thisPointIndex =pt2scurr.size(); // check if we stored variations of this point // for each variation, store the alt errors. // start at 1 since we have already filled nominal ! - if (errorBreakdown.size()) { - YAML::Node variations = errorBreakdown[thisPointIndex]; - for (const auto& variation : variations) { - const string variationName = variation.first.as(); - double eyp = variation.second["up"].as(); - double eym = variation.second["dn"].as(); - thispoint.setYErrs(eym,eyp,variationName); - } - - // LC: Not using the column representation of multiple errors right now. - //for (unsigned int ivar=1; ivar> eym >> eyp; - // thispoint.setYErrs(eym,eyp,thisvariation); - //} - } 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; } } } } + diff --git a/src/Scatter1D.cc b/src/Scatter1D.cc --- a/src/Scatter1D.cc +++ b/src/Scatter1D.cc @@ -1,30 +1,57 @@ #include "YODA/Scatter1D.h" #include "YODA/Counter.h" #include +#include "yaml-cpp/yaml.h" +#ifdef YAML_NAMESPACE +#define YAML YAML_NAMESPACE +#endif namespace YODA { /// Make a Scatter1D representation of a Histo1D Scatter1D mkScatter(const Counter& c) { Scatter1D rtn; for (const std::string& a : c.annotations()) rtn.setAnnotation(a, c.annotation(a)); rtn.setAnnotation("Type", c.type()); // might override the copied ones rtn.addPoint(c.val(), c.err()); return rtn; } + + void Scatter1D::parseVariations() { + if (this-> _variationsParsed) { return;} + if (!(this->hasAnnotation("ErrorBreakdown"))) { return; } + YAML::Node errorBreakdown; + errorBreakdown = YAML::Load(this->annotation("ErrorBreakdown")); + + if (errorBreakdown.size()) { + for (unsigned int thisPointIndex=0 ; thisPointIndex< this->numPoints() ; ++thisPointIndex){ + Point1D &thispoint = this->_points[thisPointIndex]; + YAML::Node variations = errorBreakdown[thisPointIndex]; + for (const auto& variation : variations) { + const std::string variationName = variation.first.as(); + double eyp = variation.second["up"].as(); + double eym = variation.second["dn"].as(); + thispoint.setXErrs(eym,eyp,variationName); + } + } + this-> _variationsParsed =true; + } + } + const std::vector Scatter1D::variations() const { std::vector vecvariations; for (auto &point : this->_points){ for (auto &it : point.errMap()){ //if the variation is not already in the vector, add it ! if (std::find(vecvariations.begin(), vecvariations.end(), it.first) == vecvariations.end()){ vecvariations.push_back(it.first); } } } return vecvariations; } + } diff --git a/src/Scatter2D.cc b/src/Scatter2D.cc --- a/src/Scatter2D.cc +++ b/src/Scatter2D.cc @@ -1,137 +1,166 @@ #include "YODA/Scatter2D.h" #include "YODA/Histo1D.h" #include "YODA/Profile1D.h" #include +#include "yaml-cpp/yaml.h" +#ifdef YAML_NAMESPACE +#define YAML YAML_NAMESPACE +#endif namespace YODA { /// Make a Scatter2D representation of a Histo1D Scatter2D mkScatter(const Histo1D& h, bool usefocus) { Scatter2D rtn; for (const std::string& a : h.annotations()) rtn.setAnnotation(a, h.annotation(a)); rtn.setAnnotation("Type", h.type()); // might override the copied ones for (const HistoBin1D& b : h.bins()) { const double x = usefocus ? b.xFocus() : b.xMid(); const double ex_m = x - b.xMin(); const double ex_p = b.xMax() - x; double y; try { y = b.height(); } catch (const Exception&) { // LowStatsError or WeightError y = std::numeric_limits::quiet_NaN(); } double ey; try { ey = b.heightErr(); } catch (const Exception&) { // LowStatsError or WeightError ey = std::numeric_limits::quiet_NaN(); } - const Point2D pt(x, y, ex_m, ex_p, ey, ey); + Point2D pt(x, y, ex_m, ex_p, ey, ey); + pt.setParentAO(&rtn); rtn.addPoint(pt); } assert(h.numBins() == rtn.numPoints()); return rtn; } /// Make a Scatter2D representation of a Profile1D Scatter2D mkScatter(const Profile1D& p, bool usefocus, bool usestddev) { Scatter2D rtn; for (const std::string& a : p.annotations()) rtn.setAnnotation(a, p.annotation(a)); rtn.setAnnotation("Type", p.type()); for (const ProfileBin1D& b : p.bins()) { const double x = usefocus ? b.xFocus() : b.xMid(); const double ex_m = x - b.xMin(); const double ex_p = b.xMax() - x; double y; try { y = b.mean(); } catch (const Exception&) { // LowStatsError or WeightError y = std::numeric_limits::quiet_NaN(); } double ey; try { ey = usestddev ? b.stdDev() : b.stdErr(); ///< Control y-error scheme via usestddev arg } catch (const Exception&) { // LowStatsError or WeightError ey = std::numeric_limits::quiet_NaN(); } - const Point2D pt(x, y, ex_m, ex_p, ey, ey); + //const Point2D pt(x, y, ex_m, ex_p, ey, ey); + Point2D pt(x, y, ex_m, ex_p, ey, ey); + pt.setParentAO(&rtn); rtn.addPoint(pt); } assert(p.numBins() == rtn.numPoints()); return rtn; } + // retrieve variations from annoation, parse them as YAML, and update the points + void Scatter2D::parseVariations() { + if (this-> _variationsParsed) { return; } + if (!(this->hasAnnotation("ErrorBreakdown"))) { return; } + YAML::Node errorBreakdown; + errorBreakdown = YAML::Load(this->annotation("ErrorBreakdown")); + + if (errorBreakdown.size()) { + for (unsigned int thisPointIndex=0 ; thisPointIndex< this->numPoints() ; ++thisPointIndex){ + Point2D &thispoint = this->_points[thisPointIndex]; + YAML::Node variations = errorBreakdown[thisPointIndex]; + for (const auto& variation : variations) { + const std::string variationName = variation.first.as(); + double eyp = variation.second["up"].as(); + double eym = variation.second["dn"].as(); + thispoint.setYErrs(eym,eyp,variationName); + } + } + this-> _variationsParsed =true; + } + } + const std::vector Scatter2D::variations() const { std::vector vecVariations; for (auto &point : this->_points){ for (auto &it : point.errMap()){ //if the variation is not already in the vector, add it ! if (std::find(vecVariations.begin(), vecVariations.end(), it.first) == vecVariations.end()){ vecVariations.push_back(it.first); } } } return vecVariations; } - + std::vector > Scatter2D::covarianceMatrix( bool ignoreOffDiagonalTerms) { int nPoints= this->numPoints(); //double covM[nPoints][nPoints] ={}; std::vector > covM; // initialose cov matrix to be the right shape! for (int i=0; i row; row.resize(nPoints); covM.push_back(row); } // case where only have nominal, ie total uncertainty, labelled "" (empty string) if (this->variations().size()==1) { for (int i=0; i_points[i].yErrs().first+this->_points[i].yErrs().second)/2),2); if (covM[i][i]==0 ) covM[i][i]=1; } return covM; } //more interesting case where we actually have some uncertainty breakdown! auto systList= this->variations(); for (auto sname : systList){ if (sname.length()==0) continue; std::vector< double> systErrs; systErrs.resize(nPoints); for (int i=0; i_points[i]; - try { - auto variations=point.errMap().at(sname); - systErrs[i]=(fabs(variations.first)+fabs(variations.second))*0.5 ;//up/dn are symmetrized since this method can't handle asymmetric errors - } catch (const std::exception e) { // Missing bin. - systErrs[i]=0.0; - } + try { + auto variations=point.errMap().at(sname); + systErrs[i]=(fabs(variations.first)+fabs(variations.second))*0.5 ;//up/dn are symmetrized since this method can't handle asymmetric errors + } catch (const std::exception e) { // Missing bin. + systErrs[i]=0.0; + } } if (ignoreOffDiagonalTerms || sname.find("stat") != std::string::npos || sname.find("uncor") != std::string::npos){ for (int i=0; i +#include "yaml-cpp/yaml.h" +#ifdef YAML_NAMESPACE +#define YAML YAML_NAMESPACE +#endif namespace YODA { Scatter3D mkScatter(const Histo2D& h, bool usefocus) { Scatter3D rtn; for (const std::string& a : h.annotations()) rtn.setAnnotation(a, h.annotation(a)); rtn.setAnnotation("Type", h.type()); for (size_t i = 0; i < h.numBins(); ++i) { const HistoBin2D& b = h.bin(i); /// SAME FOR ALL 2D BINS double x = b.xMid(); if (usefocus) { try { x = b.xFocus(); } catch (const LowStatsError& lse) { x = b.xMid(); } } const double exminus = x - b.xMin(); const double explus = b.xMax() - x; double y = b.yMid(); if (usefocus) { try { y = b.yFocus(); } catch (const LowStatsError& lse) { y = b.yMid(); } } const double eyminus = y - b.yMin(); const double eyplus = b.yMax() - y; /// END SAME FOR ALL 2D BINS const double z = b.height(); const double ez = b.heightErr(); rtn.addPoint(x, y, z, exminus, explus, eyminus, eyplus, ez, ez); } return rtn; } Scatter3D mkScatter(const Profile2D& h, bool usefocus, bool usestddev) { Scatter3D rtn; for (const std::string& a : h.annotations()) rtn.setAnnotation(a, h.annotation(a)); rtn.setAnnotation("Type", h.type()); for (size_t i = 0; i < h.numBins(); ++i) { const ProfileBin2D& b = h.bin(i); /// SAME FOR ALL 2D BINS double x = b.xMid(); if (usefocus) { try { x = b.xFocus(); } catch (const LowStatsError& lse) { x = b.xMid(); } } const double exminus = x - b.xMin(); const double explus = b.xMax() - x; double y = b.yMid(); if (usefocus) { try { y = b.yFocus(); } catch (const LowStatsError& lse) { y = b.yMid(); } } const double eyminus = y - b.yMin(); const double eyplus = b.yMax() - y; /// END SAME FOR ALL 2D BINS double z; try { z = b.mean(); } catch (const LowStatsError& lse) { z = std::numeric_limits::quiet_NaN(); } double ez; try { ez = usestddev ? b.stdDev() : b.stdErr(); ///< Control z-error scheme via usestddev arg } catch (const LowStatsError& lse) { ez = std::numeric_limits::quiet_NaN(); } rtn.addPoint(x, y, z, exminus, explus, eyminus, eyplus, ez, ez); } return rtn; } + void Scatter3D::parseVariations() { + if (this-> _variationsParsed) { return; } + if (!(this->hasAnnotation("ErrorBreakdown"))) { return;} + YAML::Node errorBreakdown; + errorBreakdown = YAML::Load(this->annotation("ErrorBreakdown")); + if (errorBreakdown.size()) { + for (unsigned int thisPointIndex=0 ; thisPointIndex< this->numPoints() ; ++thisPointIndex){ + Point3D &thispoint = this->_points[thisPointIndex]; + YAML::Node variations = errorBreakdown[thisPointIndex]; + for (const auto& variation : variations) { + const std::string variationName = variation.first.as(); + double eyp = variation.second["up"].as(); + double eym = variation.second["dn"].as(); + thispoint.setZErrs(eym,eyp,variationName); + } + } + this-> _variationsParsed =true; + } + } + const std::vector Scatter3D::variations() const { std::vector vecvariations; for (auto &point : this->_points){ for (auto &it : point.errMap()){ //if the variation is not already in the vector, add it ! if (std::find(vecvariations.begin(), vecvariations.end(), it.first) == vecvariations.end()){ vecvariations.push_back(it.first); } } } return vecvariations; } }