diff --git a/include/YODA/Axis2D.h b/include/YODA/Axis2D.h --- a/include/YODA/Axis2D.h +++ b/include/YODA/Axis2D.h @@ -1,663 +1,666 @@ #ifndef YODA_Axis2D_h #define YODA_Axis2D_h #include "YODA/AnalysisObject.h" #include "YODA/Exceptions.h" #include "YODA/Bin.h" #include "YODA/Utils/MathUtils.h" #include "YODA/Utils/Predicates.h" #include "YODA/Utils/BinSearcher.h" #include #include namespace YODA { /// @brief 2D bin container /// /// This class handles most of the low-level operations on an axis of bins /// arranged in a 2D grid (including gaps). template class Axis2D { public: /// Typedefs //@{ /// Bin type typedef BIN2D Bin; /// A vector containing 2D bins. Not used for searching. typedef typename std::vector Bins; // Distinguishing between single edges and edge pairs (and pairs of pairs) is useful typedef std::vector Edges; typedef std::pair EdgePair1D; typedef std::pair EdgePair2D; typedef std::vector EdgePair2Ds; // Outflow distribution lists: see outflow(int, int) typedef std::vector Outflow; typedef std::vector Outflows; //@} /// @name Constructors //@{ // Empty constructor Axis2D() : _locked(false) { reset(); } /// A constructor with specified x and y axis bin cuts. Axis2D(const Edges& xedges, const Edges& yedges) : _locked(false) { addBins(xedges, yedges); reset(); } /// Constructor accepting X/Y ranges and number of bins /// on each of the axis. Both axes are divided linearly. Axis2D(size_t nbinsX, const std::pair& rangeX, size_t nbinsY, const std::pair& rangeY) : _locked(false) { addBins(linspace(nbinsX, rangeX.first, rangeX.second), linspace(nbinsY, rangeY.first, rangeY.second)); reset(); } /// Constructor accepting a list of bins Axis2D(const Bins& bins) : _locked(false) { addBins(bins); reset(); } /// State-setting constructor for persistency Axis2D(const Bins& bins, const DBN& totalDbn, const Outflows& outflows) : _dbn(totalDbn), _outflows(outflows), _locked(false) // Does this make sense? { if (_outflows.size() != 8) { throw Exception("Axis2D outflow containers must have exactly 8 elements"); } addBins(bins); } void reset() { _dbn.reset(); _outflows.assign(8, Outflow()); for (Bin& bin : _bins) bin.reset(); _locked = false; } /// Get the number of bins. size_t numBins() const { return _bins.size(); } /// Get the number of bins on the x-axis. This is only sensible for /// perfectly regular gridded bins. For irregular binnings, this is /// the number of cuts that were necessary to grid the data. size_t numBinsX() const { return _nx; } /// Get the number of bins on the y-axis. This is only sensible for /// perfectly regular gridded bins. For irregular binnings, this is /// the number of cuts that were necessary to grid the data. size_t numBinsY() const { return _ny; } //@} // /// @name Statistics accessor functions //@{ /// @brief Get the outflow by x-index and y-index (non-const version) /// /// Indices are -1 = below range, 0 = in range, +1 = above range, e.g. (+1, /// -1) is in the "bottom right" position by being greater than the greatest /// x-edge and less than the lowest y-edge. /// Outflow& outflow(int ix, int iy) { return _outflows[_outflowIndex(ix, iy)]; } /// @brief Get the outflow by x-index and y-index (const version) /// /// Indices are -1 = below range, 0 = in range, +1 = above range, e.g. (+1, /// -1) is in the "bottom right" position by being greater than the greatest /// x-edge and less than the lowest y-edge. /// const Outflow& outflow(int ix, int iy) const { return _outflows[_outflowIndex(ix, iy)]; } /// Scale each bin as if the entire x-axis had been scaled by this factor. void scaleX(double xscale) { scaleXY(xscale, 1.0); } /// Scale each bin as if the entire y-axis had been scaled by this factor. void scaleY(double yscale) { scaleXY(1.0, yscale); } /// Scale each bin as if the entire x and y-axes had been scaled by /// their respective factors. void scaleXY(double sx, double sy) { _dbn.scaleXY(sx, sy); /// @todo Reinstate when C++11 allowed in API // for (Outflow& outflow : _outflows) // for (DBN& dbn : outflow) // dbn.scaleXY(sx, sy); for (size_t io = 0; io < _outflows.size(); ++io) { Outflow& outflow = _outflows[io]; for (size_t id = 0; id < outflow.size(); ++id) { DBN& dbn = outflow[id]; dbn.scaleXY(sx, sy); } } /// @todo Reinstate when C++11 allowed in API // for (Bin& bin : _bins) // bin.scaleXY(sx, sy); for (size_t ib = 0; ib < _bins.size(); ++ib) _bins[ib].scaleXY(sx, sy); _updateAxis(_bins); } /// Rescale as if all fill weights had been different by factor @a /// scalefactor. void scaleW(double scalefactor) { _dbn.scaleW(scalefactor); /// @todo Reinstate when C++11 allowed in API // for (Outflow& outflow : _outflows) // for (DBN& dbn : outflow) // dbn.scaleW(scalefactor); for (size_t io = 0; io < _outflows.size(); ++io) { Outflow& outflow = _outflows[io]; for (size_t id = 0; id < outflow.size(); ++id) { DBN& dbn = outflow[id]; dbn.scaleW(scalefactor); } } /// @todo Reinstate when C++11 allowed in API // for (Bin& bin : _bins) // bin.scaleW(scalefactor); for (size_t ib = 0; ib < _bins.size(); ++ib) _bins[ib].scaleW(scalefactor); _updateAxis(_bins); } /// Remove the bin at the given index. If many bins need to be /// removed, prefer eraseBins(vector[size_t] &) over many calls to this, /// as recreating the binhash is comparatively expensive. void eraseBin(size_t i) { if (i >= numBins()) throw RangeError("Bin index is out of range"); // Temporarily unlock the axis during the update _bins.erase(_bins.begin() + i); _updateAxis(_bins); } /// Erase a rectangle of bins. void eraseBins(const size_t from, const size_t to) { if (from >= numBins()) throw RangeError("Initial bin index is out of range"); if (from >= numBins()) throw RangeError("Final bin index is out of range"); Bin& fromBin = bin(from); Bin& toBin = bin(to); eraseBins(std::make_pair(fromBin.xMin(), toBin.xMax()), std::make_pair(fromBin.yMin(), toBin.yMax())); } /// Erase bins in an x- and y-range. Any bins which lie entirely within the /// region are deleted. If any part of the bin lies outside this /// range, the bin remains, so this has similar behaviour to select /// tools in vector graphics GUI packages. /// /// @todo How to test this? void eraseBins(const std::pair& xrange, const std::pair& yrange) { size_t xiLow = _binSearcherX.index(xrange.first) - 1; size_t xiHigh = _binSearcherX.index(xrange.second) - 1; size_t yiLow = _binSearcherY.index(yrange.first) - 1; size_t yiHigh = _binSearcherY.index(yrange.second) - 1; /// @todo Beware the specialisation problems with vector... std::vector deleteMask(numBins(), false); for (size_t yi = yiLow; yi < yiHigh; yi++) { for (size_t xi = xiLow; xi < xiHigh; xi++) { ssize_t i = _indexes[_index(_nx, xi, yi)]; if (i == -1 || deleteMask[i]) continue; if (bin(i).fitsInside(xrange, yrange)) deleteMask[i] = true; } } // Now we just update eraseBins(deleteMask); } /// Erase using a vector, where true represents that a bin /// will be deleted, and false means it will be kept. void eraseBins(const std::vector& deleteMask) { Bins newBins; for (size_t i = 0; i < numBins(); i++) if (!deleteMask[i]) newBins.push_back(bins(i)); _update(newBins); } /// Merge together the bin range with indices from @a from to @a to, inclusive. /// Merge a series of bins, between the bins identified by indices @a from and @a to void mergeBins(size_t from, size_t to) { // Correctness checking if (from >= numBins()) throw RangeError("Initial merge index is out of range"); if (to >= numBins()) throw RangeError("Final merge index is out of range"); if (from > to) throw RangeError("Final bin must be greater than or equal to initial bin"); if (_gapInRange(from, to)) throw RangeError("Bin ranges containing gaps cannot be merged"); if (from == to) return; // nothing to be done Bin& b = bin(from); for (size_t i = from + 1; i <= to; ++i) b.merge(_bins[i]); eraseBins(from+1, to); } /// Rebin with the same rebinning factor @a n in x and y void rebin(unsigned int n) { rebinXY(n, n); } /// Rebin with separate rebinning factors @a nx, @a ny in x and y void rebinXY(unsigned int nx, unsigned int ny) { rebinX(nx); rebinY(ny); } /// Rebin in x by factor @a nx void rebinX(unsigned int nx) { /// @todo WRITE THIS! } /// Rebin in y by factor @a ny void rebinY(unsigned int ny) { /// @todo WRITE THIS! } /// Set the axis lock state /// @todo Remove? Should not be public void _setLock(bool locked) { _locked = locked; } + /// @todo Add xMins, xMaxs, xMids, xFoci, and y-versions + + /// Return the lowest-valued bin edge along the x-axis double xMin() const { return _xRange.first; } /// Return the highest-valued bin edge along the x-axis double xMax() const { return _xRange.second; } /// Return the lowest-valued bin edge along the y-axis double yMin() const { return _yRange.first; } /// Return the highest-valued bin edge along the y-axis double yMax() const { return _yRange.second; } /// Add a bin, providing its x- and y- edge ranges void addBin(EdgePair1D xrange, EdgePair1D yrange) { _checkUnlocked(); Bins newBins = _bins; newBins.push_back(Bin(xrange, yrange)); _updateAxis(newBins); } /// Add a pre-made bin void addBin(const Bin& bin) { _checkUnlocked(); Bins newBins = _bins; newBins.push_back(bin); _updateAxis(newBins); } /// Add a vector of pre-made bins void addBins(const Bins& bins) { if (bins.size() == 0) return; _checkUnlocked(); Bins newBins = _bins; /// @todo Reinstate when C++11 allowed in API // for (const Bin& b : bins) // newBins.push_back(b); for (size_t ib = 0; ib < bins.size(); ++ib) newBins.push_back(bins[ib]); _updateAxis(newBins); } /// Add a contiguous set of bins to an axis, via their list of edges void addBins(const std::vector& xedges, const std::vector& yedges) { if (xedges.size() == 0) return; if (yedges.size() == 0) return; _checkUnlocked(); Bins newBins = _bins; for (size_t xi = 0; xi < xedges.size()-1; xi++) { for (size_t yi = 0; yi < yedges.size()-1; yi++) { const std::pair xx = std::make_pair(xedges[xi], xedges[xi+1]); const std::pair yy = std::make_pair(yedges[yi], yedges[yi+1]); // std::cout << "New bin with edges: [(" << xx.first << "," << xx.second << "), " << yy.first << "," << yy.second << ")]" << std::endl; newBins.push_back(Bin(xx, yy)); } } _updateAxis(newBins); } /// Access bin by index Bin& bin(size_t i) { return _bins[i]; } /// Access bin by index (const) const Bin& bin(size_t i) const { return _bins[i]; } /// Get the bin index of the bin containing point (x, y). int binIndexAt(double x, double y) const { size_t xi = _binSearcherX.index(x) - 1; size_t yi = _binSearcherY.index(y) - 1; if (xi > _nx) return -1; if (yi > _ny) return -1; return _indexes[_index(_nx, xi, yi)]; } /// Get the bin containing point (x, y). Bin& binAt(double x, double y) { const int ret = binIndexAt(x, y); if (ret == -1) throw RangeError("No bin found!!"); return bin(ret); } /// Get the bin containing point (x, y) (const). const Bin& binAt(double x, double y) const { const int ret = binIndexAt(x, y); if (ret == -1) throw RangeError("No bin found!!"); return bin(ret); } /// Return the total distribution (non-const) DBN& totalDbn() { return _dbn; } /// Return the total distribution (const) const DBN& totalDbn() const { return _dbn; } /// Set the total distribution: CAREFUL! void setTotalDbn(const DBN& dbn) { _dbn = dbn; } /// Return the bins vector (non-const) Bins& bins() { return _bins; } /// Return the bins vector (const) const Bins& bins() const { return _bins; } /// Equality operator (on binning only) /// @todo Change as discussed below if we expose the Axis classes for direct use // (DM: Doesn't this break the semantics of equality? As it's used only // rarely, isn't there a real case for having a "binningsCompatible" or // similar method?) bool operator == (const Axis2D& other) const { if (numBins() != other.numBins()) return false; for (size_t i = 0; i < numBins(); i++) if (!(fuzzyEquals(bin(i).xMin(), other.bin(i).xMin()) && fuzzyEquals(bin(i).xMax(), other.bin(i).xMax()) && fuzzyEquals(bin(i).yMin(), other.bin(i).yMin()) && fuzzyEquals(bin(i).yMax(), other.bin(i).yMax()))) return false; return true; } /// Non-equality operator bool operator != (const Axis2D& other) const { return ! operator == (other); } /// Addition operator Axis2D& operator += (const Axis2D& toAdd) { if (*this != toAdd) { throw LogicError("YODA::Axis2D: Cannot add axes with different binnings."); } for (size_t i = 0; i < bins().size(); ++i) { bin(i) += toAdd.bin(i); } _dbn += toAdd._dbn; return *this; } /// Subtraction operator Axis2D& operator -= (const Axis2D& toSubtract) { if (*this != toSubtract) { throw LogicError("YODA::Axis2D: Cannot add axes with different binnings."); } for (size_t i = 0; i < bins().size(); ++i) { bin(i) -= toSubtract.bin(i); } _dbn -= toSubtract._dbn; return *this; } private: void _checkUnlocked(void) { // Ensure that axis is not locked if (_locked) throw LockError("Attempting to update a locked 2D axis"); } /// Detect if there is a binning gap in the given bin index range /// @todo WRITE THIS! bool _gapInRange(size_t from, size_t to) { Bin& toBin = bin(to); Bin& fromBin = bin(from); return true; } void _updateAxis(Bins& bins) { // Deal with the case that there are no bins supplied (who called that?!) if (bins.size() == 0) { _binSearcherX = Utils::BinSearcher(); _binSearcherY = Utils::BinSearcher(); _nx = 0; _ny = 0; _xRange = std::make_pair(0, 0); _yRange = std::make_pair(0, 0); } // Sort the bins std::sort(bins.begin(), bins.end()); // Create the edges std::vector xedges, yedges, xwidths, ywidths; for (const Bin& bin : bins) { xedges.push_back(bin.xMin()); xedges.push_back(bin.xMax()); xwidths.push_back(bin.xWidth()); yedges.push_back(bin.yMin()); yedges.push_back(bin.yMax()); ywidths.push_back(bin.yWidth()); } // Sort the edges and widths std::sort(xedges.begin(), xedges.end()); std::sort(yedges.begin(), yedges.end()); std::sort(xwidths.begin(), xwidths.end()); std::sort(ywidths.begin(), ywidths.end()); // Obtain the median widths as a typical scale for uniqueness comparisons const double medianxwidth = xwidths[ (xwidths.size()-1)/2 ]; const double medianywidth = ywidths[ (ywidths.size()-1)/2 ]; // Uniqueify the bin edges in the x- and y-cut vectors, with some numerical fuzziness xedges.resize(std::unique(xedges.begin(), xedges.end(), CmpFloats(1e-3, medianxwidth)) - xedges.begin()); yedges.resize(std::unique(yedges.begin(), yedges.end(), CmpFloats(1e-3, medianywidth)) - yedges.begin()); const size_t nx = xedges.size(); const size_t ny = yedges.size(); const size_t N = nx * ny; //std::cout << "Unique Axis2D edge list sizes: nx = " << nx << ", ny = " << ny << std::endl; assert(bins.size() <= (nx-1)*(ny-1) && "Input bins vector size must agree with computed number of unique bins"); // Create a sea of indices, starting with an all-gaps configuration std::vector indexes(N, -1); // Iterate through bins and find out which Utils::BinSearcher xSearcher(xedges); Utils::BinSearcher ySearcher(yedges); for (size_t i = 0; i < bins.size(); ++i) { Bin& bin = bins[i]; // std::cout << "Bin #" << i << " edges: " // << "[(" << bin.xMin() << "," << bin.xMax() << "), " // << "(" << bin.yMin() << "," << bin.yMax() << ")] " << std::endl; const size_t xiMin= xSearcher.index(bin.xMin()) - 1; const size_t xiMax= xSearcher.index(bin.xMax()) - 1; const size_t yiMin = ySearcher.index(bin.yMin()) - 1; const size_t yiMax = ySearcher.index(bin.yMax()) - 1; // std::cout << "Sub-bin range indices: x = " << xiMin << ".." << xiMax << ", y = " << yiMin << ".." << yiMax << std::endl; // Loop over sub-bins in the edge list and assign indices / detect overlaps for (size_t xi = xiMin; xi < xiMax; xi++) { for (size_t yi = yiMin; yi < yiMax; yi++) { const size_t ii = _index(nx, xi, yi); if (indexes[ii] != -1) { std::stringstream ss; ss << "Bin edges overlap! Bin #" << i << " with edges " << "[(" << bin.xMin() << "," << bin.xMax() << "), " << "(" << bin.yMin() << "," << bin.yMax() << ")] " << "overlaps bin #" << indexes[ii] << " in sub-bin #" << ii; throw RangeError(ss.str()); } indexes[ii] = i; } } } // Job's a good'n - let's change our class. _nx = nx; _ny = ny; _xRange = std::make_pair(xedges.front(), xedges.back()); _yRange = std::make_pair(yedges.front(), yedges.back()); _indexes = indexes; _bins = bins; _binSearcherX = xSearcher; _binSearcherY = ySearcher; } /// Definition of global bin ID in terms of x and y bin IDs static size_t _index(size_t nx, size_t x, size_t y) { return y * nx + x; } /// @brief Get the outflow array index by x-index and y-index /// /// Indices are -1 = below range, 0 = in range, +1 = above range, e.g. (+1, /// -1) is in the "bottom right" position by being greater than the greatest /// x-edge and less than the lowest y-edge. static size_t _outflowIndex(int ix, int iy) { if (ix == 0 || iy == 0) throw UserError("The in-range (0,0) index pair is not a valid outflow specifier"); ix += 1; iy += 1; if (ix > 2 || iy > 2) throw UserError("Outflow index out of range: valid indices are -1, 0, 1"); size_t rtn = 3*ix + iy; // uncorrected for (0,0) index offset if (rtn > 4) rtn -= 1; // offset correction (note that raw rtn == 4 is not possible) return rtn; } /// @name Data structures //@{ /// Bins vector Bins _bins; /// Total distribution DBN _dbn; // Outflows Outflows _outflows; // Binsearcher, for searching bins Utils::BinSearcher _binSearcherX; Utils::BinSearcher _binSearcherY; EdgePair1D _xRange; EdgePair1D _yRange; // Mapping from bin-searcher indices to bin indices (allowing gaps) std::vector _indexes; // Necessary for bounds checking and indexing size_t _nx; size_t _ny; /// Whether modifying bin edges is permitted bool _locked; //@} }; } #endif diff --git a/include/YODA/Histo2D.h b/include/YODA/Histo2D.h --- a/include/YODA/Histo2D.h +++ b/include/YODA/Histo2D.h @@ -1,534 +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) // #ifndef YODA_Histo2D_h #define YODA_Histo2D_h #include "YODA/AnalysisObject.h" #include "YODA/HistoBin2D.h" #include "YODA/Dbn2D.h" #include "YODA/Axis2D.h" #include "YODA/Scatter3D.h" #include "YODA/Exceptions.h" #include #include namespace YODA { // Forward declaration class Profile2D; class Scatter3D; /// Convenience typedef typedef Axis2D Histo2DAxis; /// A two-dimensional histogram. class Histo2D : public AnalysisObject { public: /// Convenience typedefs typedef Histo2DAxis Axis; typedef Axis::Bins Bins; typedef HistoBin2D Bin; typedef Axis::Outflows Outflows; typedef std::tuple FillType; typedef FillType BinType; typedef std::shared_ptr Ptr; /// @name Constructors //@{ /// Default constructor Histo2D(const std::string& path="", const std::string& title="") : AnalysisObject("Histo2D", path, title), _axis() { } /// Constructor giving range and number of bins. Histo2D(size_t nbinsX, double lowerX, double upperX, size_t nbinsY, double lowerY, double upperY, const std::string& path="", const std::string& title="") : AnalysisObject("Histo2D", path, title), _axis(nbinsX, std::make_pair(lowerX, upperX), nbinsY, std::make_pair(lowerY, upperY)) { } /// Constructor accepting the bin edges on X and Y axis. Histo2D(const std::vector& xedges, const std::vector& yedges, const std::string& path="", const std::string& title="") : AnalysisObject("Histo2D", path, title), _axis(xedges, yedges) { } /// Constructor accepting an explicit collection of bins. Histo2D(const std::vector& bins, const std::string& path="", const std::string& title="") : AnalysisObject("Histo2D", path, title), _axis(bins) { } /// Copy constructor with optional new path /// @todo Also allow title setting from the constructor? Histo2D(const Histo2D& h, const std::string& path=""); /// A constructor from a Scatter3D's binning, with optional new path /// @todo Also allow title setting from the constructor? Histo2D(const Scatter3D& s, const std::string& path=""); /// Constructor from a Profile2D's binning, with optional new path /// @todo Also allow title setting from the constructor? Histo2D(const Profile2D& h, const std::string& path=""); /// @brief State-setting constructor /// /// Mainly intended for internal persistency use. Histo2D(const std::vector& bins, const Dbn2D& totalDbn, const Outflows& outflows, const std::string& path="", const std::string& title="") : AnalysisObject("Histo2D", path, title), _axis(bins, totalDbn, outflows) { } /// Assignment operator Histo2D& operator = (const Histo2D& h2) { AnalysisObject::operator = (h2); //< AO treatment of paths etc. _axis = h2._axis; return *this; } /// Make a copy on the stack Histo2D clone() const { return Histo2D(*this); } /// Make a copy on the heap, via 'new' Histo2D* newclone() const { return new Histo2D(*this); } //@} /// Fill dimension of this data object size_t dim() const { return 2; } /// @name Modifiers //@{ /// Fill histo with weight at (x,y) virtual void fill(double x, double y, double weight=1.0, double fraction=1.0); /// virtual void fill(const FillType & xs, double weight=1.0, double fraction=1.0) { fill(std::get<0>(xs), std::get<1>(xs), weight, fraction); } /// Fill histo x-y bin i with the given weight virtual void fillBin(size_t i, double weight=1.0, double fraction=1.0); /// @brief Reset the histogram. /// /// Keep the binning but set all bin contents and related quantities to zero void reset() { _axis.reset(); } /// Rescale as if all fill weights had been different by factor @a scalefactor. void scaleW(double scalefactor) { setAnnotation("ScaledBy", annotation("ScaledBy", 1.0) * scalefactor); _axis.scaleW(scalefactor); } /// Normalize the (visible) histo "volume" to the @a normto value. /// /// If @a includeoverflows is true, the original normalisation is computed with /// the overflow bins included, so that the resulting visible normalisation can /// be less than @a normto. This is probably what you want. void normalize(double normto=1.0, bool includeoverflows=true) { const double oldintegral = integral(includeoverflows); if (oldintegral == 0) throw WeightError("Attempted to normalize a histogram with null area"); scaleW(normto / oldintegral); } /// Scale the dimensions void scaleXY(double scaleX = 1.0, double scaleY = 1.0) { _axis.scaleXY(scaleX, scaleY); } /// @brief Bin addition operator /// /// Add a bin to an axis described by its x and y ranges. void addBin(Axis::EdgePair1D xrange, Axis::EdgePair1D yrange) { _axis.addBin(xrange, yrange); } /// @brief Bin addition operator /// /// Add a bin, possibly already populated void addBin(const Bin& bin) { _axis.addBin(bin); } /// @brief Bins addition operator /// /// Add multiple bins from edge cuts without resetting void addBins(const Axis::Edges& xcuts, const Axis::Edges& ycuts) { _axis.addBins(xcuts, ycuts); } /// @brief Bins addition operator /// /// Add multiple bins without resetting void addBins(const Bins& bins) { _axis.addBins(bins); } // /// Adding bins /// @todo TODO // void addBin(const std::vector, std::pair > > coords) { // _axis.addBin(coords); // } // /// Adding bins which is not so eloquent /// @todo TODO // void addBin(double lowX, double lowY, double highX, double highY) { // _axis.addBin(lowX, lowY, highX, highY); // } // /// Merge the bins /// @todo TODO // void mergeBins(size_t from, size_t to) { // _axis.mergeBins(from, to); // } /// Rebin the whole histo by a @a factorX in the X direction and /// @a factorY in the Y direction /// @todo TODO // void rebin(size_t factorX, size_t factorY){ // _axis.rebin(factorX, factorY); // } void eraseBin(size_t index) { _axis.eraseBin(index); } //@} public: /// @name Bin accessors //@{ + /// All bin edges on this histo's x axis + /// + /// @note This only returns the finite edges, i.e. -inf and +inf are removed + /// @todo Make the +-inf stripping controllable by a default-valued bool arg + const std::vector xEdges() const { return _axis.xEdges(); } + + /// All bin edges on this histo's y axis + /// + /// @note This only returns the finite edges, i.e. -inf and +inf are removed + /// @todo Make the +-inf stripping controllable by a default-valued bool arg + const std::vector yEdges() const { return _axis.yEdges(); } + + /// @todo Add xMins, xMaxs, xMids, xFoci, and y-versions + /// Low x edge of this histo's axis double xMin() const { return _axis.xMin(); } /// High x edge of this histo's axis double xMax() const { return _axis.xMax(); } /// Low y edge of this histo's axis double yMin() const { return _axis.yMin(); } /// High y edge of this histo's axis double yMax() const { return _axis.yMax(); } /// check if binning is the same as different Histo2D bool sameBinning(const Histo2D& h2) { return _axis == h2._axis; } /// Access the bin vector (non-const version) std::vector& bins() { return _axis.bins(); } /// Access the bin vector (const version) const std::vector& bins() const { return _axis.bins(); } /// Access a bin by index (non-const version) HistoBin2D& bin(size_t index) { return _axis.bin(index); } /// Access a bin by index (const version) const HistoBin2D& bin(size_t index) const { return _axis.bin(index); } /// Access a bin index by coordinate int binIndexAt(double x, double y) { return _axis.binIndexAt(x, y); } int binIndexAt(const BinType& t) { return _axis.binIndexAt(std::get<0>(t), std::get<1>(t)); } /// Access a bin by coordinate (const version) const HistoBin2D& binAt(double x, double y) const { return _axis.binAt(x, y); } const HistoBin2D& binAt(const BinType& t) { return _axis.binAt(std::get<0>(t), std::get<1>(t)); } /// Number of bins size_t numBins() const { return _axis.numBins(); } /// Number of bins along the x axis size_t numBinsX() const { return _axis.numBinsX(); } /// Number of bins along the y axis size_t numBinsY() const { return _axis.numBinsY(); } /// Access summary distribution, including gaps and overflows (non-const version) Dbn2D& totalDbn() { return _axis.totalDbn(); } /// Access summary distribution, including gaps and overflows (const version) const Dbn2D& totalDbn() const { return _axis.totalDbn(); } /// Set summary distribution, including gaps and overflows void setTotalDbn(const Dbn2D& dbn) { _axis.setTotalDbn(dbn); } // /// @brief Access an outflow (non-const) // /// // /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow. // /// (0,0) is not a valid overflow index pair, since it is in range for both x and y. // Dbn2D& outflow(int ix, int iy) { // std::cout << "Histo2D::outflow\n"; // return _axis.outflow(ix, iy); // } // /// @brief Access an outflow (const) // /// // /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow. // /// (0,0) is not a valid overflow index pair, since it is in range for both x and y. // const Dbn2D& outflow(int ix, int iy) const { // return _axis.outflow(ix, iy); // } //@} /// @name Whole histo data //@{ /// Get the total volume of the histogram double integral(bool includeoverflows=true) const { return sumW(includeoverflows); } /// Get the number of fills (fractional fills are possible) double numEntries(bool includeoverflows=true) const; /// Get the effective number of fills double effNumEntries(bool includeoverflows=true) const; /// Get the sum of weights in histo double sumW(bool includeoverflows=true) const; /// Get the sum of squared weights in histo double sumW2(bool includeoverflows=true) const; /// Get the mean x double xMean(bool includeoverflows=true) const; /// Get the mean y double yMean(bool includeoverflows=true) const; /// Get the variance in x double xVariance(bool includeoverflows=true) const; /// Get the variance in y double yVariance(bool includeoverflows=true) const; /// Get the standard deviation in x double xStdDev(bool includeoverflows=true) const { return std::sqrt(xVariance(includeoverflows)); } /// Get the standard deviation in y double yStdDev(bool includeoverflows=true) const { return std::sqrt(yVariance(includeoverflows)); } /// Get the standard error in x double xStdErr(bool includeoverflows=true) const; /// Get the standard error in y double yStdErr(bool includeoverflows=true) const; /// Get the RMS in x double xRMS(bool includeoverflows=true) const; /// Get the RMS in y double yRMS(bool includeoverflows=true) const; //@} /// @name Adding and subtracting histograms //@{ /// @brief Add another histogram to this one /// /// @note Adding histograms will unset any ScaledBy attribute from prevous calls to scaleW or normalize. Histo2D& operator += (const Histo2D& toAdd) { if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy"); _axis += toAdd._axis; return *this; } /// @brief Subtract another histogram from this one /// /// @note Subtracting histograms will unset any ScaledBy attribute from prevous calls to scaleW or normalize. Histo2D& operator -= (const Histo2D& toSubtract) { if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy"); _axis -= toSubtract._axis; return *this; } bool operator == (const Histo2D& other) const { return _axis == other._axis; } bool operator != (const Histo2D& other) const { return ! operator == (other); } //@} // /// @name Slicing operators // //@{ // /// @brief Create a Histo2D for the bin slice parallel to the x axis at the specified y coordinate // /// // /// Note that the created histogram will not have correctly filled underflow and overflow bins. // /// @todo It's not really *at* the specified y coord: it's for the corresponding bin row. // /// @todo Change the name! // Histo2D cutterX(double atY, const std::string& path="", const std::string& title=""); // /// @brief Create a Histo2D for the bin slice parallel to the y axis at the specified x coordinate // /// // /// Note that the created histogram will not have correctly filled underflow and overflow bins. // /// @todo It's not really *at* the specified x coord: it's for the corresponding bin row. // /// @todo Change the name! // Histo2D cutterY(double atX, const std::string& path="", const std::string& title=""); // /// X-wise Profile1D creator from Histo2D // Profile1D mkProfileX(); // /// Y-wise Profile1D creator from Histo2D // Profile1D mkProfileY(); // //@} protected: /// Access a bin by coordinate (non-const version) HistoBin2D& _binAt(double x, double y) { return _axis.binAt(x, y); } private: /// @name Bin data //@{ /// Definition of bin edges and contents Axis2D _axis; //@} }; /// Convenience typedef typedef Histo2D H2D; /// @name Combining histos: global operators //@{ /// Add two histograms inline Histo2D add(const Histo2D& first, const Histo2D& second) { Histo2D tmp = first; if (first.path() != second.path()) tmp.setPath(""); tmp += second; return tmp; } /// Add two histograms inline Histo2D operator + (const Histo2D& first, const Histo2D& second) { return add(first, second); } /// Subtract two histograms inline Histo2D subtract(const Histo2D& first, const Histo2D& second) { Histo2D tmp = first; if (first.path() != second.path()) tmp.setPath(""); tmp -= second; return tmp; } /// Subtract two histograms inline Histo2D operator - (const Histo2D& first, const Histo2D& second) { return subtract(first, second); } /// @todo Add multiply(H2, H2) -> Scatter3D? /// @brief Divide two histograms, with an uncorrelated error treatment /// /// @todo Wouldn't it be nice to be able to supply a correlation matrix or function as optional arg? /// /// @note The two histos must have _exactly_ the same binning. Scatter3D divide(const Histo2D& numer, const Histo2D& denom); /// Divide two histograms, with an uncorrelated error treatment /// /// @note The two histos must have _exactly_ the same binning. inline Scatter3D operator / (const Histo2D& numer, const Histo2D& denom) { return divide(numer, denom); } /// @brief Calculate a histogrammed efficiency ratio of two histograms /// /// @note The two histos must have _exactly_ the same binning. /// /// @note An efficiency is not the same thing as a standard division of two /// histograms: the errors are treated as correlated via binomial statistics. Scatter3D efficiency(const Histo2D& accepted, const Histo2D& total); /// @brief Calculate the asymmetry (a-b)/(a+b) of two histograms /// /// @note The two histos must have _exactly_ the same binning. inline Scatter3D asymm(const Histo2D& a, const Histo2D& b) { return (a-b) / (a+b); } //@} } #endif diff --git a/include/YODA/Profile2D.h b/include/YODA/Profile2D.h --- a/include/YODA/Profile2D.h +++ b/include/YODA/Profile2D.h @@ -1,453 +1,456 @@ // -*- 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_Profile2D_h #define YODA_Profile2D_h #include "YODA/AnalysisObject.h" #include "YODA/ProfileBin2D.h" #include "YODA/Dbn3D.h" #include "YODA/Axis2D.h" #include "YODA/Scatter3D.h" #include "YODA/Exceptions.h" #include #include namespace YODA { // Forward declarations class Histo2D; class Scatter3D; /// Convenience typedef typedef Axis2D Profile2DAxis; /// A two-dimensional profile histogram. class Profile2D : public AnalysisObject { public: /// Convenience typedefs typedef Profile2DAxis Axis; typedef Axis::Bins Bins; typedef ProfileBin2D Bin; typedef Axis::Outflows Outflows; typedef std::tuple FillType; typedef std::tuple BinType; typedef std::shared_ptr Ptr; /// @name Constructors //@{ /// Default constructor Profile2D(const std::string& path="", const std::string& title="") : AnalysisObject("Profile2D", path, title), _axis() { } /// Constructor giving range and number of bins Profile2D(size_t nbinsX, double lowerX, double upperX, size_t nbinsY, double lowerY, double upperY, const std::string& path="", const std::string& title="") : AnalysisObject("Profile2D", path, title), _axis(nbinsX, std::make_pair(lowerX, upperX), nbinsY, std::make_pair(lowerY, upperY)) { } /// Constructor giving explicit bin edges in the direction of X and Y Profile2D(const std::vector& xedges, const std::vector& yedges, const std::string& path="", const std::string& title="") : AnalysisObject("Profile2D", path, title), _axis(xedges, yedges) { } /// Constructor accepting an explicit collection of bins. Profile2D(const std::vector& bins, const std::string& path="", const std::string& title="") : AnalysisObject("Profile2D", path, title), _axis(bins) { } /// A copy constructor with optional new path /// @todo Also allow title setting from the constructor? Profile2D(const Profile2D& p, const std::string& path=""); /// A constructor from a Scatter3D's binning, with optional new path /// @todo Also allow title setting from the constructor? Profile2D(const Scatter3D& s, const std::string& path=""); /// Constructor from a Histo2D's binning, with optional new path /// @todo Also allow title setting from the constructor? Profile2D(const Histo2D& h, const std::string& path=""); /// @brief State-setting constructor /// /// Mainly intended for internal persistency use. Profile2D(const std::vector& bins, const Dbn3D& totalDbn, const Outflows& outflows, const std::string& path="", const std::string& title="") : AnalysisObject("Profile2D", path, title), _axis(bins, totalDbn, outflows) { } /// Assignment operator Profile2D& operator = (const Profile2D& p2) { AnalysisObject::operator = (p2); //< AO treatment of paths etc. _axis = p2._axis; return *this; } /// Make a copy on the stack Profile2D clone() const { return Profile2D(*this); } /// Make a copy on the heap, via 'new' Profile2D* newclone() const { return new Profile2D(*this); } //@} /// Fill dimension of this data object size_t dim() const { return 2; } /// @name Modifiers //@{ /// Fill histo by value and weight virtual void fill(double x, double y, double z, double weight=1.0, double fraction=1.0); virtual void fill(const FillType & xs, double weight=1.0, double fraction=1.0) { fill(std::get<0>(xs), std::get<1>(xs), std::get<2>(xs), weight, fraction); } /// Fill histo x-y bin i with the given z value and weight virtual void fillBin(size_t i, double z, double weight=1.0, double fraction=1.0); /// @brief Reset the histogram /// /// Keep the binning but reset the statistics void reset() { _axis.reset(); } /// Rescale as if all fill weights had been different by a @a scalefactor void scaleW(double scalefactor) { /// @todo Is this ScaledBy annotation needed? setAnnotation("ScaledBy", annotation("ScaledBy", 1.0) * scalefactor); _axis.scaleW(scalefactor); } /// Rescale as if all z values had been different by factor @a scalefactor. void scaleZ(double scalefactor) { _axis.totalDbn().scaleZ(scalefactor); /// @todo Need to rescale overflows too, when they exist. // _axis.overflow().scaleZ(scalefactor); // _axis.underflow().scaleZ(scalefactor); for (size_t i = 0; i < bins().size(); ++i) bin(i).scaleZ(scalefactor); } /// @todo TODO // /// Merge together the bin range with indices from @a from to @a to, inclusive // void mergeBins(size_t from, size_t to) { // _axis.mergeBins(from, to); // } /// @todo TODO // /// Merge every group of n bins, starting from the LHS // void rebin(size_t n) { // throw "IMPLEMENT!"; // //_axis.rebin(n); // } // /// @brief Bin addition operator // /// // /// Add a bin to the axis, described by its x and y ranges. void addBin(Axis::EdgePair1D xrange, Axis::EdgePair1D yrange) { _axis.addBin(xrange, yrange); } // /// @brief Bin addition operator // /// // /// Add a bin to the axis, possibly pre-populated void addBin(const Bin& bin) { _axis.addBin(bin); } /// @brief Bins addition operator /// /// Add multiple bins from edge cuts without resetting void addBins(const Axis::Edges& xcuts, const Axis::Edges& ycuts) { _axis.addBins(xcuts, ycuts); } /// @brief Bins addition operator /// /// Add multiple bins without resetting void addBins(const Bins& bins) { _axis.addBins(bins); } /// check if binning is the same as different Profile2D bool sameBinning(const Profile2D& p2) { return _axis == p2._axis; } /// @todo TODO // /// @brief Bin addition operator // /// // /// Add a set of bins delimiting coordinates of which are contained // /// in binLimits vector. // void addBin(const std::vector& binLimits) { // _axis.addBin(binLimits); // } void eraseBin(size_t index) { _axis.eraseBin(index); } //@} /// @name Bin accessors //@{ + /// @todo Add xMins, xMaxs, xMids, xFoci, and y-versions + + /// Low x edge of this histo's axis double xMin() const { return _axis.xMin(); } /// High x edge of this histo's axis double xMax() const { return _axis.xMax(); } /// Low y edge of this histo's axis double yMin() const { return _axis.yMin(); } /// High y edge of this histo's axis double yMax() const { return _axis.yMax(); } /// Access the bin vector (non-const) std::vector& bins() { return _axis.bins(); } /// Access the bin vector (const) const std::vector& bins() const { return _axis.bins(); } /// Access a bin by index (non-const) ProfileBin2D& bin(size_t index) { return _axis.bins()[index]; } /// Access a bin by index (const) const ProfileBin2D& bin(size_t index) const { return _axis.bins()[index]; } /// Access a bin index by coordinate int binIndexAt(double x, double y) { return _axis.binIndexAt(x, y); } int binIndexAt(const BinType& t) { return _axis.binIndexAt(std::get<0>(t), std::get<1>(t)); } /// Access a bin by coordinate (const) const ProfileBin2D& binAt(double x, double y) const { return _axis.binAt(x, y); } const ProfileBin2D& binAt(const BinType& t) const { return _axis.binAt(std::get<0>(t), std::get<1>(t)); } /// Number of bins of this axis (not counting under/over flow) size_t numBins() const { return _axis.bins().size(); } /// Number of bins along the x axis size_t numBinsX() const { return _axis.numBinsX(); } /// Number of bins along the y axis size_t numBinsY() const { return _axis.numBinsY(); } /// Access summary distribution, including gaps and overflows (non-const version) Dbn3D& totalDbn() { return _axis.totalDbn(); } /// Access summary distribution, including gaps and overflows (const version) const Dbn3D& totalDbn() const { return _axis.totalDbn(); } /// Set summary distribution, including gaps and overflows void setTotalDbn(const Dbn3D& dbn) { _axis.setTotalDbn(dbn); } // /// @brief Access an outflow (non-const) // /// // /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow. // /// (0,0) is not a valid overflow index pair, since it is in range for both x and y. // Dbn3D& outflow(int ix, int iy) { // return _axis.outflow(ix, iy); // } // /// @brief Access an outflow (const) // /// // /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow. // /// (0,0) is not a valid overflow index pair, since it is in range for both x and y. // const Dbn3D& outflow(int ix, int iy) const { // return _axis.outflow(ix, iy); // } //@} /// @name Whole histo data //@{ /// Get the number of fills (fractional fills are possible) double numEntries(bool includeoverflows=true) const; /// Get the effective number of fills double effNumEntries(bool includeoverflows=true) const; /// Get sum of weights in histo double sumW(bool includeoverflows=true) const; /// Get the sum of squared weights in histo double sumW2(bool includeoverflows=true) const; /// Get the mean x double xMean(bool includeoverflows=true) const; /// Get the mean y double yMean(bool includeoverflows=true) const; /// Get the variance in x double xVariance(bool includeoverflows=true) const; /// Get the variance in y double yVariance(bool includeoverflows=true) const; /// Get the standard deviation in x double xStdDev(bool includeoverflows=true) const { return std::sqrt(xVariance(includeoverflows)); } /// Get the standard deviation in y double yStdDev(bool includeoverflows=true) const { return std::sqrt(yVariance(includeoverflows)); } /// Get the standard error on double xStdErr(bool includeoverflows=true) const; /// Get the standard error on double yStdErr(bool includeoverflows=true) const; /// Get the RMS in x double xRMS(bool includeoverflows=true) const; /// Get the RMS in y double yRMS(bool includeoverflows=true) const; //@} /// @name Adding and subtracting histograms //@{ /// Add another profile to this one Profile2D& operator += (const Profile2D& toAdd) { if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy"); _axis += toAdd._axis; return *this; } /// Subtract another profile from this one Profile2D& operator -= (const Profile2D& toSubtract) { if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy"); _axis -= toSubtract._axis; return *this; } inline bool operator == (const Profile2D& other){ return _axis == other._axis; } inline bool operator != (const Profile2D& other){ return ! operator == (other); } //@}- protected: /// Access a bin by coordinate (non-const) ProfileBin2D& _binAt(double x, double y) { return _axis.binAt(x, y); } private: /// @name Bin data //@{ /// The bins contained in this profile histogram Axis2D _axis; //@} }; /// Convenience typedef typedef Profile2D P2D; /// @name Combining profile histos: global operators //@{ /// Add two profile histograms inline Profile2D add(const Profile2D& first, const Profile2D& second) { Profile2D tmp = first; if (first.path() != second.path()) tmp.setPath(""); tmp += second; return tmp; } /// Add two profile histograms inline Profile2D operator + (const Profile2D& first, const Profile2D& second) { return add(first,second); } /// Subtract two profile histograms inline Profile2D subtract(const Profile2D& first, const Profile2D& second) { Profile2D tmp = first; if (first.path() != second.path()) tmp.setPath(""); tmp -= second; return tmp; } /// Subtract two profile histograms inline Profile2D operator - (const Profile2D& first, const Profile2D& second) { return subtract(first,second); } /// Divide two profile histograms Scatter3D divide(const Profile2D& numer, const Profile2D& denom); /// Divide two profile histograms inline Scatter3D operator / (const Profile2D& numer, const Profile2D& denom) { return divide(numer, denom); } //@} } #endif