diff --git a/include/YODA/Scatter2D.h b/include/YODA/Scatter2D.h --- a/include/YODA/Scatter2D.h +++ b/include/YODA/Scatter2D.h @@ -1,421 +1,424 @@ // -*- 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); } //@} /////////////////////////////////////////////////// /// 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) { _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)); } /// 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)); } /// 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)); } /// 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)); } /// Insert a collection of new points void addPoints(const Points& pts) { for (const Point2D& pt : pts) addPoint(pt); } //@} /// @name Combining sets of scatter points //@{ /// @todo Better name? Make this the add operation? void combineWith(const Scatter2D& other) { addPoints(other.points()); //return *this; } /// @todo Better name? /// @todo Convert/extend to accept a Range or generic void combineWith(const std::vector& others) { for (const Scatter2D& s : others) combineWith(s); //return *this; } //@} /// Equality operator bool operator == (const Scatter2D& other) { return _points == other._points; } /// Non-equality operator bool operator != (const Scatter2D& other) { return ! operator == (other); } ////////////////////////////////// private: Points _points; }; /// 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/pyext/yoda/declarations.pxd b/pyext/yoda/declarations.pxd --- a/pyext/yoda/declarations.pxd +++ b/pyext/yoda/declarations.pxd @@ -1,1426 +1,1433 @@ 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 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 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 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/Point.pyx b/pyext/yoda/include/Point.pyx --- a/pyext/yoda/include/Point.pyx +++ b/pyext/yoda/include/Point.pyx @@ -1,133 +1,150 @@ cimport util cdef class Point(util.Base): """ A generic point with errors, used by the Scatter classes. """ cdef c.Point* pptr(self) except NULL: return self.ptr() def __dealloc__(self): cdef c.Point *p = self.pptr() if self._deallocate: del p # def __init__(self): # cutil.set_owned_ptr(self, new c.Point()) # def copy(self): # return cutil.new_owned_cls(Point, new c.Point(deref(self.pptr()))) # TODO: add clone() as mapping to (not yet existing) C++ newclone()? @property def dim(self): """None -> int Space dimension of the point (should match containing Scatter)""" return self.pptr().dim() def val(self, i): """int -> float Value on axis i""" return self.pptr().val(i) def setVal(self, i, val): """(int, float) -> None Value on axis i""" self.pptr().setVal(i, val) def errs(self, i, source=""): """int -> float Errors on axis i""" if source is None: source = "" + if isinstance(source, str): + source = source.encode('utf-8') return util.read_error_pair(self.pptr().errs(i,source)) def setErr(self, i, e, source=""): """(int, float) -> None Set symmetric errors on axis i""" if source is None: source = "" - print "LC DEBUG setErr ", e, source + if isinstance(source, str): + source = source.encode('utf-8') self.pptr().setErr(i, e, source) def setErrs(self, i, *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(i,errs[0], source) return errs=errs[0] # assert len(errs) == 2: + if isinstance(source, str): + source = source.encode('utf-8') self.pptr().setErrs(i, tuple(errs), source) def errMinus(self, i, source=""): """int -> float Minus error on axis i""" if source is None: source = "" + if isinstance(source, str): + source = source.encode('utf-8') return self.pptr().errMinus(i ,source) def setErrMinus(self, i, e, source=""): """(int, float) -> None Set minus error on axis i""" if source is None: source = "" + if isinstance(source, str): + source = source.encode('utf-8') self.pptr().setErrMinus(i, e, source) def errPlus(self, i, source=""): """int -> float Plus error on axis i""" if source is None: source = "" + if isinstance(source, str): + source = source.encode('utf-8') return self.pptr().errPlus(i, source) def setErrPlus(self, i, e, source=""): """(int, float) -> None Set plus error on axis i""" if source is None: source = "" + if isinstance(source, str): + source = source.encode('utf-8') self.pptr().setErrPlus(i, e, source) def errAvg(self, i, source=""): """int -> float Average error on axis i""" if source is None: source = "" + if isinstance(source, str): + source = source.encode('utf-8') return self.pptr().errAvg(i, source) def set(self, i, val, *es, source=""): """(int, float, float) -> None (int, float, [float, float]) -> None (int, float, float, float) -> None Set value and errors on axis i""" errs = es if source is None: source = "" if len(es) == 1: if hasattr(es[0], "__iter__"): errs = [es[0], es[0]] else: errs = es[0] # assert len(errs) == 2: + if isinstance(source, str): + source = source.encode('utf-8') self.pptr().set(i, val, errs, source) def errMap(self): """None -> {string: [float,float]} error map of this point""" return self.pptr().errMap() # def __repr__(self): # return '' % self.x 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,141 +1,165 @@ 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.x = x self.y = y self.xErrs = 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 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==None: source="" self.p2ptr().setYErrs(util.read_symmetric(val)) property x: """x coordinate""" def __get__(self): return self.p2ptr().x() def __set__(self, x): self.p2ptr().setX(x) property y: """y coordinate""" def __get__(self): return self.p2ptr().y() def __set__(self, y): self.p2ptr().setY(y) 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: This is fine because preserntly only the highest dimension supports multi-errors property xErrs: """The x errors""" def __get__(self): return util.read_error_pair(self.p2ptr().xErrs()) def __set__(self, val): self.p2ptr().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: 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)) @property def xMin(self): """The minimum x position, i.e. lowest error""" return self.p2ptr().xMin() @property def xMax(self): """The maximum x position, i.e. highest error""" return self.p2ptr().xMax() @property def yMin(self): """The minimum y position, i.e. lowest error""" return self.p2ptr().yMin() @property def yMax(self): """The maximum y position, i.e. highest error""" return self.p2ptr().yMax() property xErrAvg: def __get__(self): return self.p2ptr().xErrAvg() property yErrAvg: def __get__(self): return self.p2ptr().yErrAvg() 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/Scatter2D.pyx b/pyext/yoda/include/Scatter2D.pyx --- a/pyext/yoda/include/Scatter2D.pyx +++ b/pyext/yoda/include/Scatter2D.pyx @@ -1,224 +1,267 @@ 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.s2ptr().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 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): + """ + Construct the covariance matrix""" + return self._mknp(self.s2ptr().covarianceMatrix(False)) + + 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() + 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 - def _mknp(self, xs): - try: - import numpy - return numpy.array(xs) - except ImportError: - return xs + return self._mknp(corr) + def xVals(self): return self._mknp([p.x for p in self.points]) def xMins(self): """All x low values.""" return self._mknp([p.xMin for p in self.points]) def xMaxs(self): """All x high values.""" return self._mknp([p.xMax for p in self.points]) # TODO: xErrs def xMin(self): """Lowest x value.""" return min(self.xMins()) def xMax(self): """Highest x value.""" return max(self.xMaxs()) def yVals(self): return self._mknp([p.y for p in self.points]) def yMins(self): """All y low values.""" return self._mknp([p.yMin for p in self.points]) def yMaxs(self): """All y high values.""" return self._mknp([p.yMax for p in self.points]) # TODO: yErrs def yMin(self): """Lowest x value.""" return min(self.yMins()) def yMax(self): """Highest y value.""" return max(self.yMaxs()) ## Convenience alias S2D = Scatter2D diff --git a/src/ReaderYODA.cc b/src/ReaderYODA.cc --- a/src/ReaderYODA.cc +++ b/src/ReaderYODA.cc @@ -1,551 +1,548 @@ // -*- 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; + //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: 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 - if (!(key.find("Variations") != string::npos)) aocurr->setAnnotation(key, val); + // 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); } } 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(); + //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("Variations") != string::npos) { - YAML::Node anns = YAML::Load(s); - for (const auto& it : anns) { - assert(it.second.IsSequence()); - for (const auto& it2 : it.second) { - const string val = it2.as(); - //const string val=""; - variationscurr.push_back(val); - } - } + 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); - } - } + //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 - 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> eym >> eyp; - thispoint.setYErrs(eym,eyp,thisvariation); + // 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); - // 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> ezm >> ezp; - thispoint.setZErrs(ezm,ezp,thisvariation); - } - } 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/Scatter2D.cc b/src/Scatter2D.cc --- a/src/Scatter2D.cc +++ b/src/Scatter2D.cc @@ -1,85 +1,133 @@ #include "YODA/Scatter2D.h" #include "YODA/Histo1D.h" #include "YODA/Profile1D.h" #include 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); 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); rtn.addPoint(pt); } assert(p.numBins() == rtn.numPoints()); return rtn; } const std::vector Scatter2D::variations() const { - std::vector vecvariations; + 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); + if (std::find(vecVariations.begin(), vecVariations.end(), it.first) == vecVariations.end()){ + vecVariations.push_back(it.first); } } } - return vecvariations; + 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]; + 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 + } + if (ignoreOffDiagonalTerms || sname.find("stat") != std::string::npos || sname.find("uncor") != std::string::npos){ + for (int i=0; i #include using namespace std; namespace YODA { /// Singleton creation function Writer& WriterYODA::create() { static WriterYODA _instance; _instance.setPrecision(6); return _instance; } // Format version: // - V1/empty = make-plots annotations style // - V2 = YAML annotations static const int YODA_FORMAT_VERSION = 2; // Version-formatting helper function inline string _iotypestr(const string& baseiotype) { ostringstream os; os << "YODA_" << Utils::toUpper(baseiotype) << "_V" << YODA_FORMAT_VERSION; return os.str(); } void WriterYODA::_writeAnnotations(std::ostream& os, const AnalysisObject& ao) { os << scientific << setprecision(_precision); for (const string& a : ao.annotations()) { if (a.empty()) continue; /// @todo Write out floating point annotations as scientific notation os << a << ": " << ao.annotation(a) << "\n"; } os << "---\n"; } void WriterYODA::writeCounter(std::ostream& os, const Counter& c) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("COUNTER") << " " << c.path() << "\n"; _writeAnnotations(os, c); os << "# sumW\t sumW2\t numEntries\n"; os << c.sumW() << "\t" << c.sumW2() << "\t" << c.numEntries() << "\n"; os << "END " << _iotypestr("COUNTER") << "\n"; os.flags(oldflags); } void WriterYODA::writeHisto1D(std::ostream& os, const Histo1D& h) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("HISTO1D") << " " << h.path() << "\n"; _writeAnnotations(os, h); try { //if ( h.totalDbn().effNumEntries() > 0 ) { os << "# Mean: " << h.xMean() << "\n"; os << "# Area: " << h.integral() << "\n"; } catch (LowStatsError& e) { // } os << "# ID\t ID\t sumw\t sumw2\t sumwx\t sumwx2\t numEntries\n"; os << "Total \tTotal \t"; os << h.totalDbn().sumW() << "\t" << h.totalDbn().sumW2() << "\t"; os << h.totalDbn().sumWX() << "\t" << h.totalDbn().sumWX2() << "\t"; os << h.totalDbn().numEntries() << "\n"; os << "Underflow\tUnderflow\t"; os << h.underflow().sumW() << "\t" << h.underflow().sumW2() << "\t"; os << h.underflow().sumWX() << "\t" << h.underflow().sumWX2() << "\t"; os << h.underflow().numEntries() << "\n"; os << "Overflow\tOverflow\t"; os << h.overflow().sumW() << "\t" << h.overflow().sumW2() << "\t"; os << h.overflow().sumWX() << "\t" << h.overflow().sumWX2() << "\t"; os << h.overflow().numEntries() << "\n"; os << "# xlow\t xhigh\t sumw\t sumw2\t sumwx\t sumwx2\t numEntries\n"; for (const HistoBin1D& b : h.bins()) { os << b.xMin() << "\t" << b.xMax() << "\t"; os << b.sumW() << "\t" << b.sumW2() << "\t"; os << b.sumWX() << "\t" << b.sumWX2() << "\t"; os << b.numEntries() << "\n"; } os << "END " << _iotypestr("HISTO1D") << "\n"; os.flags(oldflags); } void WriterYODA::writeHisto2D(std::ostream& os, const Histo2D& h) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("HISTO2D") << " " << h.path() << "\n"; _writeAnnotations(os, h); try { //if ( h.totalDbn().numEntries() > 0 ) os << "# Mean: (" << h.xMean() << ", " << h.yMean() << ")\n"; os << "# Volume: " << h.integral() << "\n"; } catch (LowStatsError& e) { // } os << "# ID\t ID\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t sumwxy\t numEntries\n"; // Total distribution const Dbn2D& td = h.totalDbn(); os << "Total \tTotal \t"; os << td.sumW() << "\t" << td.sumW2() << "\t"; os << td.sumWX() << "\t" << td.sumWX2() << "\t"; os << td.sumWY() << "\t" << td.sumWY2() << "\t"; os << td.sumWXY() << "\t"; os << td.numEntries() << "\n"; // Outflows /// @todo Disabled for now, reinstate with a *full* set of outflow info to allow marginalisation os << "# 2D outflow persistency not currently supported until API is stable\n"; // for (int ix = -1; ix <= 1; ++ix) { // for (int iy = -1; iy <= 1; ++iy) { // if (ix == 0 && iy == 0) continue; // os << "Outflow\t" << ix << ":" << iy << "\t"; // const Dbn2D& d = h.outflow(ix, iy); // os << d.sumW() << "\t" << d.sumW2() << "\t"; // os << d.sumWX() << "\t" << d.sumWX2() << "\t"; // os << d.sumWY() << "\t" << d.sumWY2() << "\t"; // os << d.sumWXY() << "\t"; // os << d.numEntries() << "\n"; // } // } // Bins os << "# xlow\t xhigh\t ylow\t yhigh\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t sumwxy\t numEntries\n"; for (const HistoBin2D& b : h.bins()) { os << b.xMin() << "\t" << b.xMax() << "\t"; os << b.yMin() << "\t" << b.yMax() << "\t"; os << b.sumW() << "\t" << b.sumW2() << "\t"; os << b.sumWX() << "\t" << b.sumWX2() << "\t"; os << b.sumWY() << "\t" << b.sumWY2() << "\t"; os << b.sumWXY() << "\t"; os << b.numEntries() << "\n"; } os << "END " << _iotypestr("HISTO2D") << "\n"; os.flags(oldflags); } void WriterYODA::writeProfile1D(std::ostream& os, const Profile1D& p) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("PROFILE1D") << " " << p.path() << "\n"; _writeAnnotations(os, p); os << "# ID\t ID\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t numEntries\n"; os << "Total \tTotal \t"; os << p.totalDbn().sumW() << "\t" << p.totalDbn().sumW2() << "\t"; os << p.totalDbn().sumWX() << "\t" << p.totalDbn().sumWX2() << "\t"; os << p.totalDbn().sumWY() << "\t" << p.totalDbn().sumWY2() << "\t"; os << p.totalDbn().numEntries() << "\n"; os << "Underflow\tUnderflow\t"; os << p.underflow().sumW() << "\t" << p.underflow().sumW2() << "\t"; os << p.underflow().sumWX() << "\t" << p.underflow().sumWX2() << "\t"; os << p.underflow().sumWY() << "\t" << p.underflow().sumWY2() << "\t"; os << p.underflow().numEntries() << "\n"; os << "Overflow\tOverflow\t"; os << p.overflow().sumW() << "\t" << p.overflow().sumW2() << "\t"; os << p.overflow().sumWX() << "\t" << p.overflow().sumWX2() << "\t"; os << p.overflow().sumWY() << "\t" << p.overflow().sumWY2() << "\t"; os << p.overflow().numEntries() << "\n"; os << "# xlow\t xhigh\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t numEntries\n"; for (const ProfileBin1D& b : p.bins()) { os << b.xMin() << "\t" << b.xMax() << "\t"; os << b.sumW() << "\t" << b.sumW2() << "\t"; os << b.sumWX() << "\t" << b.sumWX2() << "\t"; os << b.sumWY() << "\t" << b.sumWY2() << "\t"; os << b.numEntries() << "\n"; } os << "END " << _iotypestr("PROFILE1D") << "\n"; os.flags(oldflags); } void WriterYODA::writeProfile2D(std::ostream& os, const Profile2D& p) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("PROFILE2D") << " " << p.path() << "\n"; _writeAnnotations(os, p); os << "# sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t sumwz\t sumwz2\t sumwxy\t numEntries\n"; // Total distribution const Dbn3D& td = p.totalDbn(); os << "Total \tTotal \t"; os << td.sumW() << "\t" << td.sumW2() << "\t"; os << td.sumWX() << "\t" << td.sumWX2() << "\t"; os << td.sumWY() << "\t" << td.sumWY2() << "\t"; os << td.sumWZ() << "\t" << td.sumWZ2() << "\t"; os << td.sumWXY() << "\t"; // << td.sumWXZ() << "\t" << td.sumWYZ() << "\t"; os << td.numEntries() << "\n"; // Outflows /// @todo Disabled for now, reinstate with a *full* set of outflow info to allow marginalisation os << "# 2D outflow persistency not currently supported until API is stable\n"; // for (int ix = -1; ix <= 1; ++ix) { // for (int iy = -1; iy <= 1; ++iy) { // if (ix == 0 && iy == 0) continue; // os << "Outflow\t" << ix << ":" << iy << "\t"; // const Dbn3D& d = p.outflow(ix, iy); // os << d.sumW() << "\t" << d.sumW2() << "\t"; // os << d.sumWX() << "\t" << d.sumWX2() << "\t"; // os << d.sumWY() << "\t" << d.sumWY2() << "\t"; // os << d.sumWZ() << "\t" << d.sumWZ2() << "\t"; // os << d.sumWXY() << "\t"; // << d.sumWXZ() << "\t" << d.sumWYZ() << "\t"; // os << d.numEntries() << "\n"; // } // } // Bins os << "# xlow\t xhigh\t ylow\t yhigh\t sumw\t sumw2\t sumwx\t sumwx2\t sumwy\t sumwy2\t sumwz\t sumwz2\t sumwxy\t numEntries\n"; for (const ProfileBin2D& b : p.bins()) { os << b.xMin() << "\t" << b.xMax() << "\t"; os << b.yMin() << "\t" << b.yMax() << "\t"; os << b.sumW() << "\t" << b.sumW2() << "\t"; os << b.sumWX() << "\t" << b.sumWX2() << "\t"; os << b.sumWY() << "\t" << b.sumWY2() << "\t"; os << b.sumWZ() << "\t" << b.sumWZ2() << "\t"; os << b.sumWXY() << "\t"; // << b.sumWXZ() << "\t" << b.sumWYZ() << "\t"; os << b.numEntries() << "\n"; } os << "END " << _iotypestr("PROFILE2D") << "\n"; os.flags(oldflags); } void WriterYODA::writeScatter1D(std::ostream& os, const Scatter1D& s) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("SCATTER1D") << " " << s.path() << "\n"; //first write the Variations, a dummy annotation which //contains the additional columns which will be written out //for sytematic variations YAML::Emitter out; out << YAML::Flow ; out << s.variations(); - os << "Variations" << ": " << out.c_str() << "\n"; + //os << "Variations" << ": " << out.c_str() << "\n"; // then write the regular annotations _writeAnnotations(os, s); std::vector variations= s.variations(); //write headers std::string headers="# xval\t "; for (const auto &source : variations){ headers+=" xerr-"+source+"\t xerr+"+source+"\t"; } os << headers << "\n"; //write points for (const Point1D& pt : s.points()) { // fill central value os << pt.x(); // fill errors for variations. The first should always be "" which is nominal. // Assumes here that all points in the Scatter have the same // variations... if not a range error will get thrown from // the point when the user tries to access a variation it // doesn't have... @todo maybe better way to do this? for (const auto &source : variations){ os << "\t" << pt.xErrMinus(source) << "\t" << pt.xErrPlus(source) ; } os << "\n"; } os << "END " << _iotypestr("SCATTER1D") << "\n"; os << flush; os.flags(oldflags); } void WriterYODA::writeScatter2D(std::ostream& os, const Scatter2D& s) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("SCATTER2D") << " " << s.path() << "\n"; //first write the Variations, a dummy annotation which //contains the additional columns which will be written out //for sytematic variations YAML::Emitter out; - out << YAML::Flow ; - out << s.variations(); - os << "Variations" << ": " << out.c_str() << "\n"; + out << YAML::Flow << YAML::BeginMap; + int counter=0; + std::vector variations= s.variations(); + //write ErrBreakdown Annotation + for (const Point2D& pt : s.points()) { + out << YAML::Key << counter; + out << YAML::Value << YAML::BeginMap; + for (const auto &source : variations){ + if (source.length()==0) continue; + out << YAML::Key << source; + out << YAML::Value << YAML::BeginMap; + out << YAML::Key << "dn" << YAML::Value << pt.yErrMinus(source); + out << YAML::Key << "up" << YAML::Value << pt.yErrPlus(source); + out << YAML::EndMap; + } + out << YAML::EndMap; + } + out << YAML::EndMap; + os << "ErrorBreakdown" << ": " << out.c_str() << "\n"; // then write the regular annotations _writeAnnotations(os, s); - std::vector variations= s.variations(); //write headers /// @todo Change ordering to {vals} {errs} {errs} ... - std::string headers="# xval\t xerr-\t xerr+\t yval\t"; - for (const auto &source : variations){ - headers+=" yerr-"+source+"\t yerr+"+source+"\t"; - } + std::string headers="# xval\t xerr-\t xerr+\t yval\t yerr-\t yerr+\t"; + //for (const auto &source : variations){ + // headers+=" yerr-"+source+"\t yerr+"+source+"\t"; + //} os << headers << "\n"; //write points for (const Point2D& pt : s.points()) { /// @todo Change ordering to {vals} {errs} {errs} ... // fill central value os << pt.x() << "\t" << pt.xErrMinus() << "\t" << pt.xErrPlus() << "\t"; os << pt.y(); // fill errors for variations. The first should always be "" which is nominal. // Assumes here that all points in the Scatter have the same // variations... if not a range error will get thrown from // the point when the user tries to access a variation it // doesn't have... @todo maybe better way to do this? - for (const auto &source : variations){ - os << "\t" << pt.yErrMinus(source) << "\t" << pt.yErrPlus(source) ; - } + //for (const auto &source : variations){ + os << "\t" << pt.yErrMinus() << "\t" << pt.yErrPlus() ; + // } os << "\n"; } - os << "END " << _iotypestr("SCATTER2D") << "\n"; + os << "END " << _iotypestr("SCATTER2D") << "\n\n"; os << flush; os.flags(oldflags); } void WriterYODA::writeScatter3D(std::ostream& os, const Scatter3D& s) { ios_base::fmtflags oldflags = os.flags(); os << scientific << showpoint << setprecision(_precision); os << "BEGIN " << _iotypestr("SCATTER3D") << " " << s.path() << "\n"; //first write the Variations, a dummy annotation which //contains the additional columns which will be written out //for sytematic variations YAML::Emitter out; out << YAML::Flow ; out << s.variations(); - os << "Variations" << ": " << out.c_str() << "\n"; // then write the regular annotations _writeAnnotations(os, s); std::vector variations= s.variations(); //write headers /// @todo Change ordering to {vals} {errs} {errs} ... std::string headers="# xval\t xerr-\t xerr+\t yval\t yerr-\t yerr+\t zval\t "; for (const auto &source : variations){ headers+=" zerr-"+source+"\t zerr+"+source+"\t"; } os << headers << "\n"; //write points for (const Point3D& pt : s.points()) { /// @todo Change ordering to {vals} {errs} {errs} ... // fill central value os << pt.x() << "\t" << pt.xErrMinus() << "\t" << pt.xErrPlus() << "\t"; os << pt.y() << "\t" << pt.yErrMinus() << "\t" << pt.yErrPlus() << "\t"; os << pt.z(); // fill errors for variations. The first should always be "" which is nominal. // Assumes here that all points in the Scatter have the same // variations... if not a range error will get thrown from // the point when the user tries to access a variation it // doesn't have... @todo maybe better way to do this? for (const auto &source : variations){ os << "\t" << pt.zErrMinus(source) << "\t" << pt.zErrPlus(source) ; } os << "\n"; } os << "END " << _iotypestr("SCATTER3D") << "\n"; os << flush; os.flags(oldflags); } }