|
|
[yoda-svn] r253 - in trunk: . include/YODA src
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Tue Aug 16 15:04:58 BST 2011
Author: buckley
Date: Tue Aug 16 15:04:57 2011
New Revision: 253
Log:
Tidying (Doxy commenting and code style), removing redundant functions, adding inlines. Hopefully no clashes this time!
Modified:
trunk/TODO
trunk/include/YODA/Axis2D.h
trunk/include/YODA/Dbn1D.h
trunk/include/YODA/Dbn2D.h
trunk/include/YODA/Histo2D.h
trunk/include/YODA/HistoBin1D.h
trunk/include/YODA/HistoBin2D.h
trunk/include/YODA/Point2D.h
trunk/include/YODA/Point3D.h
trunk/include/YODA/ProfileBin1D.h
trunk/include/YODA/Scatter2D.h
trunk/include/YODA/Scatter3D.h
trunk/src/Histo1D.cc
trunk/src/Histo2D.cc
trunk/src/Scatter3D.cc
Modified: trunk/TODO
==============================================================================
--- trunk/TODO Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/TODO Tue Aug 16 15:04:57 2011 (r253)
@@ -6,14 +6,18 @@
* Check if the weight2 treatment in Dbn1D is really correct: we're currently
multiplying it by sign(w). Hmm. (AB)
+* Rebinning: merges of n adjacent bins and global rebinning by integer factor
+ (on widths or on bin groups?) (AB for 1D, MICHAL for 2D (done))
+
+* Add copy constructors for Dbn1D/2D, Scatter2D, Histo/ProfileBin1D, Point2D/3D, HistoBin2D.
+ Rule of three: also need explicit assignment operator= for all classes
+ with copy constructors?
+
* Conversion functions to build Histo1D and Profile1D objects for slicings and
marginalisations along both X and Y directions of 2D histos. Throw an
exception if binnings are not complete grids. (MICHAL)
(done -- to be reviewed)
-* Rebinning: merges of n adjacent bins and global rebinning by integer factor
- (on widths or on bin groups?) (AB for 1D, MICHAL for 2D (done))
-
* Implement hierarchy-cascading state-setting constructor for Histo2D.(done) (MICHAL)
* Use state-setting constructors for Histo and Profile persistency via the
Modified: trunk/include/YODA/Axis2D.h
==============================================================================
--- trunk/include/YODA/Axis2D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/Axis2D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -21,819 +21,820 @@
const double LARGENUM = numeric_limits<double>::max();
const double SMALLNUM = numeric_limits<double>::min();
-
namespace YODA {
- /// @todo Use the same 2-space indentation scheme as elsewhere
- /// @brief 2D bin container and provider
- /// This class handles almost all boiler-plate operations
- /// on 2D bins (like creating axis, adding, searching, testing).
- template <typename BIN>
- class Axis2D {
- public:
-
- /// A collection of helpful typedefs
- typedef BIN Bin;
- typedef typename std::vector<BIN> Bins;
-
- /// When edge is added to the collection it must
- /// obey the following format. size_t specifies the bin this
- /// edge is a member of, a pair contains a beginning and end of the edge.
- typedef typename std::pair<size_t, std::pair<double,double> > Edge;
-
- /// A type being a basic substructure of _binHashSparse. It contains a indicator
- /// specifying the major coordinate and a collection of edges sharing the same major
- /// coordinate.
- typedef typename std::pair<double, std::vector<Edge> > EdgeCollection;
-
- /// A simple point in 2D @todo Should Point2D be used?
- typedef typename std::pair<double, double> Point;
-
- /// Segment, having a beginning and end.
- typedef typename std::pair<Point, Point> Segment;
-
- public:
- /// @name Constructors:
- //@{
-
- /// @brief Empty constructor
- /// Only added because it is required by SWIG.
- /// It doesn't make much sense to use it.
- Axis2D()
- {
- vector<Segment> edges;
- _mkAxis(edges);
- }
-
- /// Constructor provided with a vector of bin delimiters
- Axis2D(const vector<Segment>& binLimits)
- {
- _mkAxis(binLimits);
- }
-
- ///Most standard constructor, should be self-explanatory
- Axis2D(size_t nbinsX, double lowerX, double upperX, size_t nbinsY, double lowerY, double upperY)
- {
- vector<Segment> binLimits;
- double coeffX = (upperX - lowerX)/(double)nbinsX;
- double coeffY = (upperY - lowerX)/(double)nbinsY;
-
- for(double i=lowerX; i < upperX; i+=coeffX) {
- for(double j=lowerY; j < upperY; j+=coeffY) {
- binLimits.push_back(make_pair(make_pair(i, j),
- make_pair((double)(i+coeffX), (double)(j+coeffY))));
- }
- }
- _mkAxis(binLimits);
- }
+ /// @brief 2D bin container and provider
+ /// This class handles almost all boiler-plate operations
+ /// on 2D bins (like creating axis, adding, searching, testing).
+ template <typename BIN>
+ class Axis2D {
+ public:
+
+ // A few helpful typedefs
+ typedef BIN Bin;
+ typedef typename std::vector<BIN> Bins;
+
+ /// When edge is added to the collection it must
+ /// obey the following format. size_t specifies the bin this
+ /// edge is a member of, a pair contains a beginning and end of the edge.
+ typedef typename std::pair<size_t, std::pair<double,double> > Edge;
+
+ /// A type being a basic substructure of _binHashSparse. It contains a indicator
+ /// specifying the major coordinate and a collection of edges sharing the same major
+ /// coordinate.
+ typedef typename std::pair<double, std::vector<Edge> > EdgeCollection;
+
+ /// A simple point in 2D @todo Should Point2D be used?
+ typedef typename std::pair<double, double> Point;
+
+ /// Segment, having a beginning and end.
+ typedef typename std::pair<Point, Point> Segment;
+
+
+ public:
+
+ /// @name Constructors:
+ //@{
+
+ /// @brief Empty constructor
+ /// Only added because it is required by SWIG.
+ /// It doesn't make much sense to use it.
+ /// @todo Really "required"? Eliminate the requirement in the SWIG mapping.
+ Axis2D() {
+ vector<Segment> edges;
+ _mkAxis(edges);
+ }
- /// @brief A bin constructor
- /// Creates an Axis2D from existing bins
- Axis2D(Bins bins)
- {
- vector<Segment> binLimits;
- for(size_t i=0; i < bins.size(); ++i) binLimits.push_back(make_pair(make_pair(bins[i].xMin(), bins[i].yMin()),
- make_pair(bins[i].xMax(), bins[i].yMax())));
- _mkAxis(binLimits);
- for(size_t i=0; i < _bins.size(); ++i) _bins[i] = bins[i];
- }
- //@}
-
- /// @name Addition operators:
- //@{
-
- /// @brief Bin addition operator
- /// This operator is provided a vector of limiting
- /// points in the format required by _mkAxis().
- /// It should be noted that there is nothing special about
- /// the initiation stage of Axis2D, and the edges can be added
- /// online if they meet all the requirements of non-degeneracy.
- /// No merging is supported, and I don't think it should before the support
- /// for merging for '+' operator (codewise it should be the same thing).
- void addBin(const vector<Segment>& binLimits)
- {
- _mkAxis(binLimits);
- }
-
- /// @brief Bin addition operator
- /// This operator is supplied with whe extreamal coordinates of just
- /// one bin. It then launches the standard bin addition procedure.
- void addBin(double lowX, double lowY, double highX, double highY)
- {
- vector<Segment> coords;
- coords.push_back(make_pair(make_pair(lowX, lowY), make_pair(highX, highY)));
-
- addBin(coords);
- }
- //@}
-
- /// @name Some helper functions:
- //@{
-
- /// @brief Bin merger
- /// Try to merge a certain amount of bins
- void mergeBins(size_t from, size_t to) {
- HistoBin2D& start = bin(from);
- HistoBin2D& end = bin(to);
- HistoBin2D temp = start;
- start.isReal = false;
-
- if(start.midpoint().first > end.midpoint().first ||
- start.midpoint().second > end.midpoint().second)
- throw GridError("The start/end bins are wrongly defined.");
-
-
- for(size_t y = (*_binHashSparse.first._cache.lower_bound(start.yMin())).second; y <= (*_binHashSparse.first._cache.lower_bound(end.yMin())).second; ++y)
- {
- for(size_t x = 0; x < _binHashSparse.first[y].second.size(); ++x) {
- if((_binHashSparse.first[y].second[x].second.first > start.xMin() ||
- fuzzyEquals(_binHashSparse.first[y].second[x].second.first, start.xMin())) &&
- (_binHashSparse.first[y].second[x].second.second < end.xMax() ||
- fuzzyEquals(_binHashSparse.first[y].second[x].second.second, end.xMax())) &&
- bin(_binHashSparse.first[y].second[x].first).isReal)
- {
- temp += bin(_binHashSparse.first[y].second[x].first);
- bin(_binHashSparse.first[y].second[x].first).isReal = false;
- }
- }
- }
- _addEdge(temp.edges(), _binHashSparse, false);
- _bins.push_back(temp);
- _binHashSparse.first.regenCache();
- _binHashSparse.second.regenCache();
- _regenDelimiters();
- }
- /// @brief Checks if our bins form a grid.
- /// This function uses a neat property of _binHashSparse.
- /// If it is containing a set of edges forming a grid without
- /// gaps in the middle it will have the same number of edges in the
- /// inner subcaches and half of this amount in the outer (grid boundary)
- /// subcaches. This makes isGriddy() a very, very fast function.
- /// @todo Is the name appropriate?
- bool isGriddy() const
- {
-
- /// Check if the number of edges parallel to X axis
- /// is proper in every subcache.
- unsigned int sizeX = _binHashSparse.first[0].second.size();
- for(unsigned int i=1; i < _binHashSparse.first.size(); i++) {
- if(i == _binHashSparse.first.size() - 1) {
- if(_binHashSparse.first[i].second.size() != sizeX) {
- return false;
- }
- }
- else if(_binHashSparse.first[i].second.size() != 2*sizeX) {
- return false;
- }
- }
+ /// Constructor provided with a vector of bin delimiters
+ Axis2D(const vector<Segment>& binLimits) {
+ _mkAxis(binLimits);
+ }
- /// Do the same for edges parallel to Y axis.
- unsigned int sizeY = _binHashSparse.second[0].second.size();
- for(unsigned int i=1; i < _binHashSparse.second.size(); i++) {
- if(i!= _binHashSparse.second.size() - 1) {
- if(2*sizeY != _binHashSparse.second[i].second.size()) {
- return false;
- }
- }
- else if(_binHashSparse.second[i].second.size() != sizeY) return -1;
- }
- /// If everything is proper, announce it.
- return true;
+ /// Most standard constructor, should be self-explanatory
+ Axis2D(size_t nbinsX, double lowerX, double upperX, size_t nbinsY, double lowerY, double upperY) {
+ vector<Segment> binLimits;
+ double coeffX = (upperX - lowerX)/(double)nbinsX;
+ double coeffY = (upperY - lowerX)/(double)nbinsY;
+ for (double i = lowerX; i < upperX; i += coeffX) {
+ for(double j = lowerY; j < upperY; j += coeffY) {
+ binLimits.push_back(make_pair(make_pair(i, j), make_pair((double)(i+coeffX), (double)(j+coeffY))));
}
+ }
+ _mkAxis(binLimits);
+ }
- /// Return a total number of bins in a Histo
- unsigned int numBinsTotal() const
- {
- return _bins.size();
- }
- /// Get inf(X) (non-const version)
- double lowEdgeX()
- {
- return _lowEdgeX;
- }
+ /// @brief A bin constructor
+ /// Creates an Axis2D from existing bins
+ Axis2D(Bins bins) {
+ vector<Segment> binLimits;
+ for (size_t i = 0; i < bins.size(); ++i) {
+ binLimits.push_back(make_pair(make_pair(bins[i].xMin(), bins[i].yMin()),
+ make_pair(bins[i].xMax(), bins[i].yMax())));
+ }
+ _mkAxis(binLimits);
+ for (size_t i = 0; i < _bins.size(); ++i) {
+ _bins[i] = bins[i];
+ }
+ }
- /// Get sup(X) (non-const version)
- double highEdgeX()
- {
- return _highEdgeX;
- }
+ //@}
- /// Get inf(Y) (non-const version)
- double lowEdgeY()
- {
- return _lowEdgeY;
- }
+ /// @name Addition operators:
+ //@{
- /// Get sup(Y) (non-const version)
- double highEdgeY()
- {
- return _highEdgeY;
- }
+ /// @brief Bin addition operator
+ /// This operator is provided a vector of limiting
+ /// points in the format required by _mkAxis().
+ /// It should be noted that there is nothing special about
+ /// the initiation stage of Axis2D, and the edges can be added
+ /// online if they meet all the requirements of non-degeneracy.
+ /// No merging is supported, and I don't think it should before the support
+ /// for merging for '+' operator (codewise it should be the same thing).
+ void addBin(const vector<Segment>& binLimits) {
+ _mkAxis(binLimits);
+ }
- /// Get inf(X) (const version)
- const double lowEdgeX() const
- {
- return _lowEdgeX;
- }
+ /// @brief Bin addition operator
+ /// This operator is supplied with whe extreamal coordinates of just
+ /// one bin. It then launches the standard bin addition procedure.
+ void addBin(double lowX, double lowY, double highX, double highY) {
+ vector<Segment> coords;
+ coords.push_back(make_pair(make_pair(lowX, lowY), make_pair(highX, highY)));
+ addBin(coords);
+ }
- /// Get sup(X) (const version)
- const double highEdgeX() const
- {
- return _highEdgeX;
- }
+ //@}
- /// Get inf(Y)
- const double lowEdgeY() const
- {
- return _lowEdgeY;
- }
- ///Get sup(Y)
- const double highEdgeY() const
- {
- return _highEdgeY;
- }
+ /// @name Helper functions
+ //@{
- /// Get the bins from an Axis (non-const version)
- Bins& bins()
- {
- return _bins;
- }
+ /// @brief Bin merging
+ /// Try to merge a certain amount of bins
+ void mergeBins(size_t from, size_t to) {
+ HistoBin2D& start = bin(from);
+ HistoBin2D& end = bin(to);
+ HistoBin2D temp = start;
+ start.isReal = false;
+
+ if (start.midpoint().first > end.midpoint().first ||
+ start.midpoint().second > end.midpoint().second)
+ throw GridError("The start/end bins are wrongly defined.");
+
+ /// @todo Tidy!
+ for (size_t y = (*_binHashSparse.first._cache.lower_bound(start.yMin())).second;
+ y <= (*_binHashSparse.first._cache.lower_bound(end.yMin())).second; ++y) {
+ for (size_t x = 0; x < _binHashSparse.first[y].second.size(); ++x) {
+ /// @todo Tidy!
+ if ((_binHashSparse.first[y].second[x].second.first > start.xMin() ||
+ fuzzyEquals(_binHashSparse.first[y].second[x].second.first, start.xMin())) &&
+ (_binHashSparse.first[y].second[x].second.second < end.xMax() ||
+ fuzzyEquals(_binHashSparse.first[y].second[x].second.second, end.xMax())) &&
+ bin(_binHashSparse.first[y].second[x].first).isReal)
+ {
+ temp += bin(_binHashSparse.first[y].second[x].first);
+ bin(_binHashSparse.first[y].second[x].first).isReal = false;
+ }
+ }
+ }
+ _addEdge(temp.edges(), _binHashSparse, false);
+ _bins.push_back(temp);
+
+ _binHashSparse.first.regenCache();
+ _binHashSparse.second.regenCache();
+ _regenDelimiters();
+ }
- /// Get the bins from an Axis (const version)
- const Bins& bins() const
- {
- return _bins;
- }
- /// Get a bin with a given index (non-const version)
- BIN& bin(size_t index)
- {
- if(index >= _bins.size()) throw RangeError("Bin index out of range.");
- return _bins[index];
+ /// @brief Checks if our bins form a grid.
+ /// This function uses a neat property of _binHashSparse.
+ /// If it is containing a set of edges forming a grid without
+ /// gaps in the middle it will have the same number of edges in the
+ /// inner subcaches and half of this amount in the outer (grid boundary)
+ /// subcaches. This makes isGriddy() a very, very fast function.
+ /// @todo Is the name appropriate? Should it be private?
+ bool isGriddy() const {
+
+ /// Check if the number of edges parallel to X axis
+ /// is proper in every subcache.
+ size_t sizeX = _binHashSparse.first[0].second.size();
+ for (size_t i = 1; i < _binHashSparse.first.size(); ++i) {
+ if (i == _binHashSparse.first.size() - 1) {
+ if (_binHashSparse.first[i].second.size() != sizeX) {
+ return false;
+ }
}
-
- /// Get a bin with a given index (const version)
- const BIN& bin(size_t index) const
- {
- if(index >= _bins.size()) throw RangeError("Bin index out of range.");
- return _bins[index];
+ else if (_binHashSparse.first[i].second.size() != 2*sizeX) {
+ return false;
}
+ }
- /// Get a bin at given coordinates (non-const version)
- BIN& binByCoord(double x, double y)
- {
- int ret = _findBinIndex(x, y);
- if(ret != -1) return bin(ret);
- else throw RangeError("No bin found!!");
+ /// Do the same for edges parallel to Y axis.
+ size_t sizeY = _binHashSparse.second[0].second.size();
+ for (size_t i=1; i < _binHashSparse.second.size(); ++i) {
+ if (i != _binHashSparse.second.size() - 1) {
+ if (2*sizeY != _binHashSparse.second[i].second.size()) {
+ return false;
+ }
}
-
- /// Get a bin at given coordinates (const version)
- const BIN& binByCoord(double x, double y) const
- {
- int ret = _findBinIndex(x, y);
- if(ret != -1) return bin(ret);
- else throw RangeError("No bin found!!");
+ else if (_binHashSparse.second[i].second.size() != sizeY) {
+ return -1;
}
+ }
- /// Get a bin at given coordinates (non-const version)
- BIN& binByCoord(pair<double, double>& coords)
- {
- int ret = _findBinIndex(coords.first, coords.second);
- if(ret != -1) return bin(ret);
- else throw RangeError("No bin found!!");
- }
+ /// If everything is proper, announce it.
+ return true;
+ }
- /// Get a bin at given coordinates (const version)
- const BIN& binByCoord(pair<double, double>& coords) const
- {
- int ret = _findBinIndex(coords.first, coords.second);
- if(ret != -1) return bin(ret);
- else throw RangeError("No bin found!!");
- }
- /// Get a total distribution (non-const version)
- Dbn2D& totalDbn()
- {
- return _dbn;
- }
+ /// Get the total number of bins
+ unsigned int numBinsTotal() const {
+ return _bins.size();
+ }
- /// Get a total distribution (const version)
- const Dbn2D& totalDbn() const
- {
- return _dbn;
- }
- /// Get the overflow distribution (non-const version)
- Dbn2D& overflow()
- {
- return _overflow;
- }
+ /// Get the value of the lowest x-edge on the axis
+ double lowEdgeX() const {
+ return _lowEdgeX;
+ }
- /// Get the overflow distribution (const version)
- const Dbn2D& overflow() const
- {
- return _overflow;
- }
- /// Get the underflow distribution (non-const version)
- Dbn2D& underflow()
- {
- return _underflow;
- }
+ /// Get the value of the highest x-edge on the axis
+ double highEdgeX() const {
+ return _highEdgeX;
+ }
- /// Get the underflow distribution (const version)
- const Dbn2D& underflow() const
- {
- return _underflow;
- }
- /// Get bin index from external classes (const version)
- int getBinIndex(double coordX, double coordY) const
- {
- return _findBinIndex(coordX, coordY);
- }
+ /// Get the value of the lowest y-edge on the axis
+ double lowEdgeY() const {
+ return _lowEdgeY;
+ }
- /// Resetts the axis statistics ('fill history')
- void reset()
- {
- _dbn.reset();
- _underflow.reset();
- _overflow.reset();
- for (size_t i=0; i<_bins.size(); i++) _bins[i].reset();
- }
- /// @brief Axis scaler
- /// Scales the axis with a given scale. If no scale is given, assumes
- /// identity transform.
- void scaleXY(double scaleX = 1.0, double scaleY = 1.0)
- {
- // Two loops are put on purpose, just to protect
- // against improper _binHashSparse
- for (unsigned int i=0; i < _binHashSparse.first.size(); i++) {
- _binHashSparse.first[i].first *= scaleY;
- for (unsigned int j=0; j < _binHashSparse.first[i].second.size(); j++){
- _binHashSparse.first[i].second[j].second.first *=scaleX;
- _binHashSparse.first[i].second[j].second.second *=scaleX;
- }
- }
- for (unsigned int i=0; i < _binHashSparse.second.size(); i++) {
- _binHashSparse.second[i].first *= scaleX;
- for (unsigned int j=0; j < _binHashSparse.second[i].second.size(); j++){
- _binHashSparse.second[i].second[j].second.first *=scaleY;
- _binHashSparse.second[i].second[j].second.second *=scaleY;
- }
- }
+ /// Get the value of the highest y-edge on the axis
+ double highEdgeY() const {
+ return _highEdgeY;
+ }
- /// Regenerate the bin edges cache.
- _binHashSparse.first.regenCache();
- _binHashSparse.second.regenCache();
-
- /// Now, as we have the map rescaled, we need to update the bins
- for (size_t i = 0; i < _bins.size(); ++i) {
- _bins[i].scaleXY(scaleX, scaleY);
- }
- _dbn.scaleXY(scaleX, scaleY);
- _underflow.scaleXY(scaleX, scaleY);
- _overflow.scaleXY(scaleX, scaleY);
- // Making sure that we have correct boundaries set after rescaling
- _regenDelimiters();
- }
+ /// Get the bins (non-const version)
+ Bins& bins() {
+ return _bins;
+ }
+ /// Get the bins (const version)
+ const Bins& bins() const {
+ return _bins;
+ }
- /// Scales the bin weights
- void scaleW(double scalefactor) {
- _dbn.scaleW(scalefactor);
- _underflow.scaleW(scalefactor);
- _overflow.scaleW(scalefactor);
- for (size_t i=0; i<_bins.size(); i++) _bins[i].scaleW(scalefactor);
- }
- //@}
- /// @name Operators:
- //@{
+ /// Get the bin with a given index (non-const version)
+ BIN& bin(size_t index) {
+ if (index >= _bins.size()) throw RangeError("Bin index out of range.");
+ return _bins[index];
+ }
- /// Equality operator
- bool operator == (const Axis2D& other) const
- {
- return _binHashSparse == other._binHashSparse;
- }
+ /// Get the bin with a given index (const version)
+ const BIN& bin(size_t index) const {
+ if (index >= _bins.size()) throw RangeError("Bin index out of range.");
+ return _bins[index];
+ }
- /// Non-equality operator
- bool operator != (const Axis2D& other) const
- {
- return ! operator == (other);
- }
- /// @brief Addition operator
- /// At this stage it is only possible to add two histograms with
- /// the same binnings. Compatible but not equal binning to come soon.
- Axis2D<BIN>& operator += (const Axis2D<BIN>& toAdd)
- {
- if (*this != toAdd) {
- throw LogicError("YODA::Histo1D: Cannot add axes with different binnings.");
- }
- for (unsigned int i=0; i < bins().size(); i++) bins().at(i) += toAdd.bins().at(i);
+ /// Get a bin at given coordinates (non-const version)
+ BIN& binByCoord(double x, double y) {
+ const int ret = _findBinIndex(x, y);
+ if (ret != -1) return bin(ret);
+ else throw RangeError("No bin found!!");
+ }
+
+ /// Get a bin at given coordinates (const version)
+ const BIN& binByCoord(double x, double y) const {
+ const int ret = _findBinIndex(x, y);
+ if (ret != -1) return bin(ret);
+ else throw RangeError("No bin found!!");
+ }
+
+ /// Get a bin at given coordinates (non-const version)
+ BIN& binByCoord(pair<double, double>& coords) {
+ return binByCoord(coords.first, coords.second);
+ }
+
+ /// Get a bin at given coordinates (const version)
+ const BIN& binByCoord(pair<double, double>& coords) const {
+ return binByCoord(coords.first, coords.second);
+ }
+
+
+ /// Get a total distribution (non-const version)
+ Dbn2D& totalDbn() {
+ return _dbn;
+ }
+
+ /// Get a total distribution (const version)
+ const Dbn2D& totalDbn() const {
+ return _dbn;
+ }
+
- _dbn += toAdd._dbn;
- _underflow += toAdd._underflow;
- _overflow += toAdd._overflow;
- return *this;
+ // /// Get the overflow distribution (non-const version)
+ // Dbn2D& overflow()
+ // {
+ // return _overflow;
+ // }
+
+ // /// Get the overflow distribution (const version)
+ // const Dbn2D& overflow() const
+ // {
+ // return _overflow;
+ // }
+
+ // /// Get the underflow distribution (non-const version)
+ // Dbn2D& underflow()
+ // {
+ // return _underflow;
+ // }
+
+ // /// Get the underflow distribution (const version)
+ // const Dbn2D& underflow() const
+ // {
+ // return _underflow;
+ // }
+
+
+ /// Get bin index from external classes
+ int getBinIndex(double coordX, double coordY) const
+ {
+ return _findBinIndex(coordX, coordY);
+ }
+
+ /// Reset the axis statistics
+ void reset()
+ {
+ _dbn.reset();
+ // _underflow.reset();
+ // _overflow.reset();
+ for (size_t i = 0; i < _bins.size(); ++i) {
+ _bins[i].reset();
+ }
+ }
+
+
+ /// @brief Axis scaler
+ /// Scales the axis with a given scale. If no scale is given, assumes
+ /// identity transform.
+ void scaleXY(double scaleX = 1.0, double scaleY = 1.0) {
+ // Two loops are put on purpose, just to protect
+ // against improper _binHashSparse
+ for (size_t i = 0; i < _binHashSparse.first.size(); ++i) {
+ _binHashSparse.first[i].first *= scaleY;
+ for (size_t j = 0; j < _binHashSparse.first[i].second.size(); ++j) {
+ _binHashSparse.first[i].second[j].second.first *= scaleX;
+ _binHashSparse.first[i].second[j].second.second *= scaleX;
+ }
+ }
+ for (size_t i = 0; i < _binHashSparse.second.size(); ++i) {
+ _binHashSparse.second[i].first *= scaleX;
+ for (size_t j = 0; j < _binHashSparse.second[i].second.size(); ++j){
+ _binHashSparse.second[i].second[j].second.first *= scaleY;
+ _binHashSparse.second[i].second[j].second.second *= scaleY;
+ }
+ }
+
+ /// Regenerate the bin edges cache.
+ _binHashSparse.first.regenCache();
+ _binHashSparse.second.regenCache();
+
+ /// Now, as we have the map rescaled, we need to update the bins
+ for (size_t i = 0; i < _bins.size(); ++i) {
+ _bins[i].scaleXY(scaleX, scaleY);
+ }
+ _dbn.scaleXY(scaleX, scaleY);
+ // _underflow.scaleXY(scaleX, scaleY);
+ // _overflow.scaleXY(scaleX, scaleY);
+
+ // Making sure that we have correct boundaries set after rescaling
+ _regenDelimiters();
+ }
+
+
+ /// Scales the bin weights
+ void scaleW(double scalefactor) {
+ _dbn.scaleW(scalefactor);
+ // _underflow.scaleW(scalefactor);
+ // _overflow.scaleW(scalefactor);
+ for (size_t i = 0; i < _bins.size(); ++i) {
+ _bins[i].scaleW(scalefactor);
+ }
+ }
+ //@}
+
+
+ /// @name Operators
+ //@{
+
+ /// Equality operator
+ bool operator == (const Axis2D& other) const {
+ return _binHashSparse == other._binHashSparse;
+ }
+
+ /// Non-equality operator
+ bool operator != (const Axis2D& other) const {
+ return ! operator == (other);
+ }
+
+ /// @brief Addition operator
+ /// At this stage it is only possible to add histograms with the same binnings.
+ /// @todo Compatible but not equal binning to come soon.
+ Axis2D<BIN>& operator += (const Axis2D<BIN>& toAdd) {
+ if (*this != toAdd) {
+ throw LogicError("YODA::Histo1D: Cannot add axes with different binnings.");
+ }
+ for (size_t i = 0; i < bins().size(); ++i) {
+ bins().at(i) += toAdd.bins().at(i);
+ }
+ _dbn += toAdd._dbn;
+ // _underflow += toAdd._underflow;
+ // _overflow += toAdd._overflow;
+ return *this;
+ }
+
+
+ /// Subtraction operator
+ ///
+ /// At this stage it is only possible to subtract histograms with the same binnings.
+ /// @todo Compatible but not equal binning to come soon.
+ Axis2D<BIN>& operator -= (const Axis2D<BIN>& toSubtract) {
+ if (*this != toSubtract) {
+ throw LogicError("YODA::Histo1D: Cannot add axes with different binnings.");
+ }
+ for (size_t i = 0; i < bins().size(); ++i) {
+ bins().at(i) -= toSubtract.bins().at(i);
+ }
+ _dbn -= toSubtract._dbn;
+ // _underflow -= toSubtract._underflow;
+ // _overflow -= toSubtract._overflow;
+ return *this;
+ }
+
+ //@}
+
+
+ private:
+
+ /// @brief Segment validator function
+ ///
+ /// This a 'dispatcher' function. It checks if the segment in question
+ /// is vertical or horizontal and launches the appropriate function
+ /// searching for cuts in the prope direction. Since it operates on
+ /// a vector of segments it is prepared to act on arbitrarly large sets of edges,
+ /// in practice usually being four sides of a rectangular bin.
+ ///
+ /// Notice that it will never be checked, in the current state, if there is a cut
+ /// in edges in the edgeset. This imposes the requirement of provide the program
+ /// with non-degenerate bins, until checking in implemented. However, it doesn't seem
+ /// to be needed, as edges are not generated by a user.
+ ///
+ /// This is also a perfect place to parallelize, if required.
+ bool _validateEdge(vector<Segment>& edges) {
+ /// Setting the return variable. True means that no cuts were detected.
+ bool ret = true;
+
+ /// Looping over all the edges provided
+ for (size_t i = 0; i < edges.size(); ++i) {
+ /// If the X coordinate of the starting point is the same
+ /// as X coordinate of the ending one, checks if there are cuts
+ /// on this vertical segment.
+ if (fuzzyEquals(edges[i].first.first, edges[i].second.first)) {
+ ret = _findCutsY(edges[i]);
+ }
+
+ /// Check if the segment is horizontal and is it cutting any bin that already exists
+ else if(fuzzyEquals(edges[i].first.second, edges[i].second.second)) {
+ ret = _findCutsX(edges[i]);
+ }
+
+ /// This is a check that discards the bin if it is not a rectangle
+ /// composed of vertical and horizontal segments.
+ else {
+ ret = false;
+ }
+
+ /// If a cut was detected, return false immediately.
+ if (!ret) return false;
+ }
+ /// If no cuts were detected in any of the edges, tell the launching function about this
+ return true;
+ }
+
+
+ /// @brief A binary search function
+ ///
+ /// This is conceptually the same implementation as in STL
+ /// but because it returns the index of an element in a vector,
+ /// it is easier to use in our case than the STL implementation
+ /// that returns a pointer at the element.
+ size_t _binaryS(Utils::cachedvector<EdgeCollection>& toSearch,
+ double value, size_t lower, size_t higher) {
+ // Just a check if such degenerate situation happens
+ if (lower == higher) return lower;
+
+ // Choose a midpoint that will be our pivot
+ size_t where = (higher+lower)/2;
+
+ // Launch the same procedure on half of the range above the pivot
+ if (value >= toSearch[where].first) {
+ if (where == toSearch.size() - 1) return where;
+ if (value <= toSearch[where+1].first) return where;
+ return _binaryS(toSearch, value, where, higher);
+ }
+
+ // This is not a redundant check, because of the nature of int division.
+ if (where == 0) return where;
+
+ // Check if the value is somewhere in-between an element at the position
+ // in question and an element at a lower position. If so, return an index
+ // to the current positon.
+ if(value >= toSearch[where-1].first) return where;
+
+ // If none of the above occurs, the value must be smaller that the element
+ // at the current position. In such case, launch the search on half of the
+ // interval below the current position.
+ return _binaryS(toSearch, value, lower, where);
+ }
+
+
+ /// @brief Function that finds cuts of horizontal edges.
+ ///
+ /// A specialised function that tries to exploit the fact
+ /// that edges can cut one another only on right angles
+ /// to the highest extent. The inner workings are explained
+ /// in the comments placed in the function body.
+ bool _findCutsX(Segment& edge) {
+ // Look up the limits of search in the _binHashSparse
+ // structure. We are not interested in the vertical edges
+ // that are before the beginning of our segment or after its end.
+ size_t i = _binaryS(_binHashSparse.second, edge.first.first, 0, _binHashSparse.second.size());
+ size_t end = _binaryS(_binHashSparse.second, edge.second.first, 0, _binHashSparse.second.size());
+
+ for (; i < end; ++i) {
+ // Scroll through all the vertical segments with a given X coordinate
+ // and see if any of those fulfills the cutting requirement. If it does,
+ // announce it.
+ for (size_t j = 0; j < _binHashSparse.second[i].second.size(); ++j) {
+ // Note that we are not taking into account the edges touching the
+ // segment in question. That's because sides of a bin touch
+ if (_binHashSparse.second[i].second[j].second.first < edge.first.second &&
+ _binHashSparse.second[i].second[j].second.second > edge.first.second &&
+ !fuzzyEquals(_binHashSparse.second[i].second[j].second.first, edge.first.second) &&
+ !fuzzyEquals(_binHashSparse.second[i].second[j].second.second, edge.first.second)) {
+ return false;
+ }
}
+ }
+ /// If none of the existing edges is cutting this edge, announce it
+ return true;
+ }
- /// Subtraction operator
- Axis2D<BIN>& operator -= (const Axis2D<BIN>& toSubtract)
- {
- if (*this != toSubtract) {
- throw LogicError("YODA::Histo1D: Cannot add axes with different binnings.");
- }
- for (unsigned int i=0; i < bins().size(); i++) bins().at(i) -= toSubtract.bins().at(i);
- _dbn -= toSubtract._dbn;
- _underflow -= toSubtract._underflow;
- _overflow -= toSubtract._overflow;
- return *this;
- }
- //@}
-
-
- private:
-
- /// @brief Segment validator function
- /// This a 'dispatcher' function. It checks if the segment in question
- /// is vertical or horizontal and launches the appropriate function
- /// searching for cuts in the prope direction. Since it operates on
- /// a vector of segments it is prepared to act on arbitrarly large sets of edges,
- /// in practice usually being four sides of a rectangular bin.
- ///
- /// Notice that it will never be checked, in the current state, if there is a cut
- /// in edges in the edgeset. This imposes the requirement of provide the program
- /// with non-degenerate bins, until checking in implemented. However, it doesn't seem
- /// to be needed, as edges are not generated by a user.
- ///
- /// This is also a perfect place to paralellize the program, if required.
- bool _validateEdge(vector<Segment>& edges)
- {
- /// Setting the return variable. True means that no cuts were detected.
- bool ret = true;
-
- /// Looping over all the edges provided
- for(unsigned int i=0; i < edges.size(); i++) {
- /// If the X coordinate of the starting point is the same
- /// as X coordinate of the ending one, checks if there are cuts
- /// on this vertical segment.
- if(fuzzyEquals(edges[i].first.first, edges[i].second.first)) ret = _findCutsY(edges[i]);
-
- /// Check if the segment is horizontal and is it cutting any bin that already exists
- else if(fuzzyEquals(edges[i].first.second, edges[i].second.second)) ret = _findCutsX(edges[i]);
-
- /// This is a check that discards the bin if it is not a rectangle
- /// composed of vertical and horizontal segments.
- else ret = false;
-
- /// If a cut was detected, say it. There is no point in checking other edges
- /// in the set.
- if(!ret) return false;
- }
- /// If no cuts were detected in any of the edges, tell the launching function about this
- return true;
+ /// @brief Function that finds cuts of vertical edges
+ ///
+ /// For a detailed description, see the documentation for _findCutsX().
+ bool _findCutsY(Segment& edge) {
+ size_t i = _binaryS(_binHashSparse.first, edge.first.second, 0, _binHashSparse.first.size());
+ size_t end = _binaryS(_binHashSparse.first, edge.second.second, 0, _binHashSparse.first.size());
+ for (; i < end; ++i) {
+ for (size_t j = 0; j < _binHashSparse.first[i].second.size(); ++j) {
+ if (_binHashSparse.first[i].second[j].second.first < edge.first.first &&
+ _binHashSparse.first[i].second[j].second.second > edge.first.first &&
+ !fuzzyEquals(_binHashSparse.first[i].second[j].second.first, edge.first.first) &&
+ !fuzzyEquals(_binHashSparse.first[i].second[j].second.second, edge.first.first)) {
+ return false;
+ }
}
+ }
+ return true;
+ }
+
+
+ /// @brief Function executed when a set of edges is dropped.
+ ///
+ /// It does not have any information about which edge in the set
+ /// had failed the check. If this is needed such additional information
+ /// can be readily implemented.
+ void _dropEdge(vector<Segment>& edges) {
+ std::cerr << "A set of edges was dropped." << endl;
+ }
- /// @brief A binary search function
- /// This is conceptually the same implementation as in STL
- /// but because it returns the index of an element in a vector,
- /// it is easier to use in our case than the STL implementation
- /// that returns a pointer at the element.
- size_t _binaryS(Utils::cachedvector<EdgeCollection>& toSearch,
- double value, size_t lower, size_t higher)
- {
- /// Just a check if such degenerate situation happens
- if(lower == higher) return lower;
-
- /// Choose a midpoint that will be our pivot
- size_t where = (higher+lower)/2;
-
- /// Launch the same procedure on half of the range above the pivot
- if(value >= toSearch[where].first) {
- if(where == toSearch.size() - 1) return where;
- if(value <= toSearch[where+1].first) return where;
- return _binaryS(toSearch, value, where, higher);
- }
- /// This is not a redundant check, because
- /// of the nature of int division.
- if (where == 0) return where;
-
- /// Check if the value is somewhere inbetween
- /// an element at the position in question and
- /// an element at a lower position. If so, return
- /// an index to the current positon.
- if(value >= toSearch[where-1].first) return where;
-
- /// If none of the above occurs, the value must
- /// be smaller that the element at the current
- /// position. In such case, launch the search on
- /// half of the interval below the current position.
- return _binaryS(toSearch, value, lower, where);
- }
-
- /// @brief Function that finds cuts of horizontal edges.
- /// A specialised function that tries to exploit the fact
- /// that edges can cut one another only on right angles
- /// to the highest extent. The inner workings are explained
- /// in the comments placed in the function body.
- bool _findCutsX(Segment& edge) {
- /// Look up the limits of search in the _binHashSparse
- /// structure. We are not interested in the vertical edges
- /// that are before the beginning of our segment or after its end.
- size_t i = _binaryS(_binHashSparse.second, edge.first.first, 0, _binHashSparse.second.size());
- size_t end = _binaryS(_binHashSparse.second, edge.second.first, 0, _binHashSparse.second.size());
-
- for(; i < end; i++) {
-
- /// Scroll through all the vertical segments with a given X coordinate
- /// and see if any of those fulfills the cutting requirement. If it does,
- /// announce it.
- for(unsigned int j = 0; j < _binHashSparse.second[i].second.size(); j++) {
- /// Note that we are not taking into account the edges touching the
- /// segment in question. That's because sides of a bin touch
- if(_binHashSparse.second[i].second[j].second.first < edge.first.second &&
- _binHashSparse.second[i].second[j].second.second > edge.first.second &&
- !fuzzyEquals(_binHashSparse.second[i].second[j].second.first, edge.first.second) &&
- !fuzzyEquals(_binHashSparse.second[i].second[j].second.second, edge.first.second)) {
- return false;
- }
- }
- }
- /// If none of the existing edges is cutting this edge, announce it
- return true;
- }
-
- /// @brief Function that finds cuts of vertical edges
- /// For a detailed descpription, please look into
- /// documentation for _findCutsX().
- bool _findCutsY(Segment& edge) {
- size_t i = _binaryS(_binHashSparse.first, edge.first.second, 0, _binHashSparse.first.size());
- size_t end = _binaryS(_binHashSparse.first, edge.second.second, 0, _binHashSparse.first.size());
- for(; i < end; i++) {
-
- for(unsigned int j = 0; j < _binHashSparse.first[i].second.size(); j++) {
- if(_binHashSparse.first[i].second[j].second.first < edge.first.first &&
- _binHashSparse.first[i].second[j].second.second > edge.first.first &&
- !fuzzyEquals(_binHashSparse.first[i].second[j].second.first, edge.first.first) &&
- !fuzzyEquals(_binHashSparse.first[i].second[j].second.second, edge.first.first)) {
- return false;
- }
- }
- }
- return true;
- }
-
- /// @brief Function executed when a set of edges is dropped.
- /// It does not have any information about which edge in the set
- /// had failed the check. If this is needed such additional information
- /// can be readily implemented.
- void _dropEdge(vector<Segment>& edges) {
- std::cerr << "A set of edges was dropped." << endl;
- }
-
- /// @brief Bin adder.
- /// It contains all the commands that need to executed
- /// to properly add a bin. Specifially edges are added to
- /// the edge cache (_binHashSparse) and a bin is created from
- /// those edges.
- void _addEdge(vector<Segment> edges, pair<Utils::cachedvector<EdgeCollection>,
- Utils::cachedvector<EdgeCollection> >& binHash, bool addBin = true) {
- /// Check if there was no mistake made when adding segments to a vector.
- if(edges.size() != 4) throw Exception("The segments supplied don't describe a full bin!");
-
- /// This is the part in charge of adding each of the segments
- /// to the edge cache. Segments are assumed to be validated
- /// beforehand.
- for(unsigned int j=0; j < edges.size(); j++) {
- /// Association made for convinience.
- Segment edge = edges[j];
-
- /// Do the following if the edge is vertical.
- /// Those two cases need to be distinguished
- /// because of the way in which edge cache is structured.
- if(edge.first.first == edge.second.first) {
- /// See if our edge has the same X coordinate as any other
- /// edge that is currently in the cache.
-
- /// Keeps the status of the search
- bool found = false;
-
- /// There is only a certain set of X coordinates that we need to sweep
- /// to check if our segment has the same X coordinate. Find them.
- size_t i = _binaryS(binHash.second, edge.first.first, 0, binHash.second.size())-1;
- if(i < 0) i = 0;
- size_t end = i+3;
-
- /// For the coordinates in range, check if one of them is an X coordinate of
- /// the sement.
- for(; i < binHash.second.size() && i < end ; i++) {
- /// If this is the case, do what is needed to be done.
- if(fuzzyEquals(binHash.second[i].first, edge.first.first)) {
- binHash.second[i].second.push_back(make_pair(_bins.size(),make_pair(edge.first.second, edge.second.second)));
- found = true;
- break;
- }
- }
-
- /// If no edge with the same X coordinate exist, create
- /// a new subhash at the X coordinate of a segment.
- if(!found) {
- vector<Edge> temp;
- temp.push_back(make_pair(_bins.size(), make_pair(edge.first.second, edge.second.second)));
- binHash.second.push_back(make_pair(edge.first.first,temp));
- sort(binHash.second.begin(), binHash.second.end());
- }
- }
-
- /// See the vertical case for description of a horizontal one
- else if(edge.first.second == edge.second.second) {
- bool found = false;
- size_t i = _binaryS(binHash.first, edge.first.second, 0, binHash.first.size())-1;
- if(i < 0) i = 0;
- size_t end = i+3;
- for(; i < binHash.first.size() && i < end; i++) {
- if(fuzzyEquals(binHash.first[i].first, edge.first.second)) {
- binHash.first[i].second.push_back(make_pair(_bins.size(),make_pair(edge.first.first, edge.second.first)));
- found = true;
- }
- }
- if(!found) {
- vector<Edge> temp;
- temp.push_back(make_pair(_bins.size(), make_pair(edge.first.first, edge.second.first)));
- binHash.first.push_back(make_pair(edge.second.second, temp));
- sort(binHash.first.begin(), binHash.first.end());
- }
- }
+ /// @brief Bin adder
+ ///
+ /// It contains all the commands that need to executed
+ /// to properly add a bin. Specifially edges are added to
+ /// the edge cache (_binHashSparse) and a bin is created from
+ /// those edges.
+ void _addEdge(vector<Segment> edges, pair<Utils::cachedvector<EdgeCollection>,
+ Utils::cachedvector<EdgeCollection> >& binHash, bool addBin = true) {
+ // Check if there was no mistake made when adding segments to a vector.
+ if (edges.size() != 4) throw Exception("The segments supplied don't describe a full bin!");
+
+ // This is the part in charge of adding each of the segments
+ // to the edge cache. Segments are assumed to be validated
+ // beforehand.
+ for (size_t j = 0; j < edges.size(); ++j) {
+ // Association made for convinience.
+ Segment edge = edges[j];
+
+ // Do the following if the edge is vertical.
+ // Those two cases need to be distinguished
+ // because of the way in which edge cache is structured.
+ if (edge.first.first == edge.second.first) {
+ // See if our edge has the same X coordinate as any other
+ // edge that is currently in the cache.
+
+ // Keeps the status of the search
+ bool found = false;
+
+ // There is only a certain set of X coordinates that we need to sweep
+ // to check if our segment has the same X coordinate. Find them.
+ size_t i = _binaryS(binHash.second, edge.first.first, 0, binHash.second.size()) - 1;
+ if (i < 0) i = 0;
+ size_t end = i + 3;
+
+ // For the coordinates in range, check if one of them is an X
+ // coordinate of the segment.
+ for (; i < binHash.second.size() && i < end ; ++i) {
+ /// If this is the case, do what is needed to be done.
+ if (fuzzyEquals(binHash.second[i].first, edge.first.first)) {
+ binHash.second[i].second.push_back(make_pair(_bins.size(),make_pair(edge.first.second, edge.second.second)));
+ found = true;
+ break;
}
- /// Now, create a bin with the edges provided
- if(addBin) _bins.push_back(BIN(edges));
+ }
+
+ // If no edge with the same X coordinate exist, create
+ // a new subhash at the X coordinate of a segment.
+ if (!found) {
+ vector<Edge> temp;
+ temp.push_back(make_pair(_bins.size(), make_pair(edge.first.second, edge.second.second)));
+ binHash.second.push_back(make_pair(edge.first.first,temp));
+ sort(binHash.second.begin(), binHash.second.end());
+ }
}
- /// @brief Orientation fixer
- /// Check if the orientation of an edge is proper
- /// for the rest of the algorithm to work on, and if it is not
- /// fix it.
- void _fixOrientation(Segment& edge) {
- if(fuzzyEquals(edge.first.first, edge.second.first)) {
- if(edge.first.second > edge.second.second) {
- double temp = edge.second.second;
- edge.second.second = edge.first.second;
- edge.first.second = temp;
- }
- }
- else if(edge.first.first > edge.second.first) {
- double temp = edge.first.first;
- edge.first.first = edge.second.first;
- edge.second.first = temp;
+ // See the vertical case for description of a horizontal one
+ else if (edge.first.second == edge.second.second) {
+ bool found = false;
+ size_t i = _binaryS(binHash.first, edge.first.second, 0, binHash.first.size())-1;
+ if (i < 0) i = 0;
+ size_t end = i+3;
+ for (; i < binHash.first.size() && i < end; ++i) {
+ if (fuzzyEquals(binHash.first[i].first, edge.first.second)) {
+ binHash.first[i].second.push_back(make_pair(_bins.size(),make_pair(edge.first.first, edge.second.first)));
+ found = true;
}
+ }
+ if (!found) {
+ vector<Edge> temp;
+ temp.push_back(make_pair(_bins.size(), make_pair(edge.first.first, edge.second.first)));
+ binHash.first.push_back(make_pair(edge.second.second, temp));
+ sort(binHash.first.begin(), binHash.first.end());
+ }
}
+ }
- /// @brief Axis creator
- /// The top-level function taking part in the process of
- /// adding edges. Creating an axis is the same operation
- /// for it as adding new bins so it can be as well used to
- /// add some custom bins.
- ///
- /// It accepts two extremal points of a rectangle
- /// (top-right and bottom-left) as input.
- void _mkAxis(const vector<Segment>& binLimits) {
-
- /// For each of the rectangles
- for(unsigned int i=0; i < binLimits.size(); i++) {
- /// Produce the segments that a rectangle is composed of
- Segment edge1 =
- make_pair(binLimits[i].first,
- make_pair(binLimits[i].first.first, binLimits[i].second.second));
- Segment edge2 =
- make_pair(make_pair(binLimits[i].first.first, binLimits[i].second.second),
- binLimits[i].second);
- Segment edge3 =
- make_pair(make_pair(binLimits[i].second.first, binLimits[i].first.second),
- binLimits[i].second);
- Segment edge4 =
- make_pair(binLimits[i].first,
- make_pair(binLimits[i].second.first, binLimits[i].first.second));
-
- /// Check if they are made properly
- _fixOrientation(edge1); _fixOrientation(edge2);
- _fixOrientation(edge3); _fixOrientation(edge4);
-
- /// Add all the segments to a vector
- vector<Segment> edges;
- edges.push_back(edge1); edges.push_back(edge2);
- edges.push_back(edge3); edges.push_back(edge4);
-
- /// And check if a bin is a proper one, if it is, add it.
- if(_validateEdge(edges)) _addEdge(edges, _binHashSparse);
- else _dropEdge(edges);
+ // Now, create a bin with the edges provided
+ if (addBin) _bins.push_back(BIN(edges));
+ }
- }
- /// Setting all the caches
- _binHashSparse.first.regenCache();
- _binHashSparse.second.regenCache();
- _regenDelimiters();
- }
-
- /// @brief Plot extrema (re)generator.
- /// Since scrolling through every bin is an expensive
- /// operation to do every time we need the limits of
- /// the plot, there are caches set up. This function
- /// regenerates them. It should be run after any change is made
- /// to bin layout.
- void _regenDelimiters() {
- double highEdgeX = SMALLNUM;
- double highEdgeY = SMALLNUM;
- double lowEdgeX = LARGENUM;
- double lowEdgeY = LARGENUM;
-
- /// Scroll through the bins and set the delimiters.
- for(unsigned int i=0; i < _bins.size(); i++) {
- if(_bins[i].xMin() < lowEdgeX) lowEdgeX = _bins[i].xMin();
- if(_bins[i].xMax() > highEdgeX) highEdgeX = _bins[i].xMax();
- if(_bins[i].yMin() < lowEdgeY) lowEdgeY = _bins[i].yMin();
- if(_bins[i].yMax() > highEdgeY) highEdgeY = _bins[i].yMax();
- }
-
- cout << "LowX: " << lowEdgeX << " LowY: " << lowEdgeY << " HighX: " << highEdgeX << " HighY: " << highEdgeY << endl;
- _lowEdgeX = lowEdgeX;
- _highEdgeX = highEdgeX;
- _lowEdgeY = lowEdgeY;
- _highEdgeY = highEdgeY;
- }
-
- /// @brief Bin index finder
- /// Looks through all the bins to see which
- /// one contains the point of interest.
- int _findBinIndex(double coordX, double coordY) const
- {
- for(size_t i=0; i < _bins.size(); i++) {
- if(_bins[i].xMin() <= coordX && _bins[i].xMax() >= coordX &&
- _bins[i].yMin() <= coordY && _bins[i].yMax() >= coordY &&
- _bins[i].isReal) return i;
- }
- return -1;
- }
-
- private:
-
- /// Bins contained in this histogram
- Bins _bins;
-
- /// Underflow distribution
- Dbn2D _underflow;
-
- /// Overflow distribution
- Dbn2D _overflow;
-
- /// The total distribution
- Dbn2D _dbn;
-
- /// @brief Bin hash structure
- /// First in pair is holding the horizontal edges indexed by first.first
- /// which is an y coordinate. The last pair specifies x coordinates (begin, end) of
- /// the horizontal edge.
- /// Analogous for the second member of the pair.
- ///
- /// For the fullest description, see typedefs at the beginning
- /// of this file.
- std::pair<Utils::cachedvector<EdgeCollection>,
- Utils::cachedvector<EdgeCollection> >
- _binHashSparse;
-
- /// Low/High edges:
- double _highEdgeX, _highEdgeY, _lowEdgeX, _lowEdgeY;
-
- };
-
- /// Addition operator
- template <typename BIN>
- Axis2D<BIN> operator + (const Axis2D<BIN>& first, const Axis2D<BIN>& second)
- {
- Axis2D<BIN> tmp = first;
- tmp += second;
- return tmp;
+ /// @brief Orientation fixer
+ ///
+ /// Check if the orientation of an edge is proper
+ /// for the rest of the algorithm to work on, and if it is not
+ /// fix it.
+ void _fixOrientation(Segment& edge) {
+ if (fuzzyEquals(edge.first.first, edge.second.first)) {
+ if (edge.first.second > edge.second.second) {
+ double temp = edge.second.second;
+ edge.second.second = edge.first.second;
+ edge.first.second = temp;
+ }
+ }
+ else if (edge.first.first > edge.second.first) {
+ double temp = edge.first.first;
+ edge.first.first = edge.second.first;
+ edge.second.first = temp;
+ }
}
- /// Subtraction operator
- template <typename BIN>
- Axis2D<BIN> operator - (const Axis2D<BIN>& first, const Axis2D<BIN>& second)
+
+ /// @brief Axis creator
+ ///
+ /// The top-level function taking part in the process of adding
+ /// edges. Creating an axis is the same operation for it as adding new bins
+ /// so it can be as well used to add some custom bins.
+ ///
+ /// It accepts two extremal points of a rectangle (top-right and
+ /// bottom-left) as input.
+ void _mkAxis(const vector<Segment>& binLimits) {
+
+ // For each of the rectangles
+ for (size_t i = 0; i < binLimits.size(); ++i) {
+ // Produce the segments that a rectangle is composed of
+ Segment edge1 =
+ make_pair(binLimits[i].first,
+ make_pair(binLimits[i].first.first, binLimits[i].second.second));
+ Segment edge2 =
+ make_pair(make_pair(binLimits[i].first.first, binLimits[i].second.second),
+ binLimits[i].second);
+ Segment edge3 =
+ make_pair(make_pair(binLimits[i].second.first, binLimits[i].first.second),
+ binLimits[i].second);
+ Segment edge4 =
+ make_pair(binLimits[i].first,
+ make_pair(binLimits[i].second.first, binLimits[i].first.second));
+
+ // Check if they are made properly
+ _fixOrientation(edge1); _fixOrientation(edge2);
+ _fixOrientation(edge3); _fixOrientation(edge4);
+
+ // Add all the segments to a vector
+ vector<Segment> edges;
+ edges.push_back(edge1); edges.push_back(edge2);
+ edges.push_back(edge3); edges.push_back(edge4);
+
+ // And check if a bin is a proper one, if it is, add it.
+ if (_validateEdge(edges)) _addEdge(edges, _binHashSparse);
+ else _dropEdge(edges);
+ }
+
+ // Setting all the caches
+ _binHashSparse.first.regenCache();
+ _binHashSparse.second.regenCache();
+ _regenDelimiters();
+ }
+
+
+ /// @brief Plot extrema (re)generator.
+ ///
+ /// Since scrolling through every bin is an expensive operation to do every
+ /// time we need the limits of the plot, there are caches set up. This
+ /// function regenerates them. It should be run after any change is made to
+ /// bin layout.
+ void _regenDelimiters() {
+ double highEdgeX = SMALLNUM;
+ double highEdgeY = SMALLNUM;
+ double lowEdgeX = LARGENUM;
+ double lowEdgeY = LARGENUM;
+
+ // Scroll through the bins and set the delimiters.
+ for (size_t i = 0; i < _bins.size(); ++i) {
+ if (_bins[i].xMin() < lowEdgeX) lowEdgeX = _bins[i].xMin();
+ if (_bins[i].xMax() > highEdgeX) highEdgeX = _bins[i].xMax();
+ if (_bins[i].yMin() < lowEdgeY) lowEdgeY = _bins[i].yMin();
+ if (_bins[i].yMax() > highEdgeY) highEdgeY = _bins[i].yMax();
+ }
+
+ // cout << "LowX: " << lowEdgeX << " LowY: " << lowEdgeY << " HighX: " << highEdgeX << " HighY: " << highEdgeY << endl;
+ _lowEdgeX = lowEdgeX;
+ _highEdgeX = highEdgeX;
+ _lowEdgeY = lowEdgeY;
+ _highEdgeY = highEdgeY;
+ }
+
+ /// @brief Bin index finder
+ ///
+ /// Looks through all the bins to see which one contains the point of
+ /// interest.
+ int _findBinIndex(double coordX, double coordY) const
{
- Axis2D<BIN> tmp = first;
- tmp -= second;
- return tmp;
+ for (size_t i=0; i < _bins.size(); ++i) {
+ if(_bins[i].xMin() <= coordX && _bins[i].xMax() >= coordX &&
+ _bins[i].yMin() <= coordY && _bins[i].yMax() >= coordY &&
+ _bins[i].isReal) return i;
+ }
+ return -1;
}
+
+
+ private:
+
+ /// Bins contained in this histogram
+ Bins _bins;
+
+ /// Underflow distribution
+ // /// @todo Need many (binned?) outflows instead. This can't work at the moment.
+ // Dbn2D _underflow;
+ // Dbn2D _overflow;
+
+ /// The total distribution
+ Dbn2D _dbn;
+
+ /// @brief Bin hash structure
+ /// First in pair is holding the horizontal edges indexed by first.first
+ /// which is an y coordinate. The last pair specifies x coordinates (begin, end) of
+ /// the horizontal edge.
+ /// Analogous for the second member of the pair.
+ ///
+ /// For the fullest description, see typedefs at the beginning
+ /// of this file.
+ std::pair<Utils::cachedvector<EdgeCollection>, Utils::cachedvector<EdgeCollection> > _binHashSparse;
+
+ /// Low/high edges
+ double _highEdgeX, _highEdgeY, _lowEdgeX, _lowEdgeY;
+
+ };
+
+
+ /// @name Operators
+ //@{
+
+ /// Addition operator
+ template <typename BIN>
+ inline Axis2D<BIN> operator + (const Axis2D<BIN>& first, const Axis2D<BIN>& second) {
+ Axis2D<BIN> tmp = first;
+ tmp += second;
+ return tmp;
+ }
+
+
+ /// Subtraction operator
+ template <typename BIN>
+ inline Axis2D<BIN> operator - (const Axis2D<BIN>& first, const Axis2D<BIN>& second)
+ {
+ Axis2D<BIN> tmp = first;
+ tmp -= second;
+ return tmp;
+ }
+
+ //@}
+
+
}
#endif
Modified: trunk/include/YODA/Dbn1D.h
==============================================================================
--- trunk/include/YODA/Dbn1D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/Dbn1D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -42,6 +42,9 @@
_sumWX2 = sumWX2;
}
+
+ /// @todo Add copy constructor
+
//@}
Modified: trunk/include/YODA/Dbn2D.h
==============================================================================
--- trunk/include/YODA/Dbn2D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/Dbn2D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -14,6 +14,7 @@
/// @name Constructors
//@{
+ /// Default constructor
Dbn2D() {
reset();
}
@@ -28,6 +29,9 @@
_sumWXY(sumWXY)
{ }
+
+ /// @todo Add copy constructor
+
//@}
Modified: trunk/include/YODA/Histo2D.h
==============================================================================
--- trunk/include/YODA/Histo2D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/Histo2D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -18,6 +18,7 @@
namespace YODA {
+
// Forward declaration
class Scatter3D;
@@ -57,49 +58,51 @@
_axis(binedges)
{ }
+
/// Copy constructor with optional new path
- Histo2D(const Histo2D& h, const std::string& path="");
+ Histo2D(const Histo2D& h, const std::string& path="")
+ : AnalysisObject("Histo2D", (path.size() == 0) ? h.path() : path, h, h.title())
+ {
+ _axis = h._axis;
+ }
- /// @brief Bin constructor
- /// A constructor that accepts a vector of bins and produces a meaningful Histo2D out of those
- Histo2D(const std::vector<HistoBin2D> bins, const std::string& path="", const std::string& title="")
- : AnalysisObject("Histo2D", path, title), _axis(bins)
- { }
-
- //@}
+ /// @todo Add binning constructors from Scatter3D (and Profile2D when it exists)
- public:
+ /// @brief State-setting constructor
+ /// Mainly intended for internal persistency use.
+ Histo2D(const std::vector<HistoBin2D> bins,
- /// @name Persistency hooks
- //@{
+ const std::string& path="", const std::string& title="")
+ : AnalysisObject("Histo2D", path, title), _axis(bins)
+ { }
- /// Get name of the analysis object type, for persisting
- std::string _aotype() const { return "Histo2D"; }
+ //@}
- /// Set the state of the histo object, for unpersisting
- /// @todo Need to set annotations (do that on AO), all-histo Dbns, and dbns for every bin. Delegate!
- // void _setstate() = 0;
- //@}
+ public:
/// @name Modifiers
//@{
- /// Fill histo by value and weight
+ /// @brief Fill histo by value and weight
+ ///
+ /// It relies on Axis2D for searching procedures and complains loudly if no
+ /// bin was found.
int fill(double x, double y, double weight=1.0);
/// @brief Reset the histogram.
/// Keep the binning but set all bin contents and related quantities to zero
- virtual void reset() {
+ void reset() {
_axis.reset();
}
- ///Check if is a grid:
+ /// Check if this binning is a grid
+ /// @todo Really expose this on the public API?
int isGriddy() {
- return _axis.isGriddy();
+ return _axis.isGriddy();
}
/// Rescale as if all fill weights had been different by factor @a scalefactor.
@@ -122,32 +125,36 @@
void addBin(double lowX, double lowY, double highX, double highY) {
_axis.addBin(lowX, lowY, highX, highY);
}
-
+
/// Merge the bins
void mergeBins(size_t from, size_t to) {
_axis.mergeBins(from, to);
}
+
//@}
+
public:
/// @name Bin accessors
//@{
- /// Low edges of this histo's axis
+ /// Low x edge of this histo's axis
double lowEdgeX() const {
return _axis.lowEdgeX();
}
+ /// Low y edge of this histo's axis
double lowEdgeY() const {
return _axis.lowEdgeY();
}
- /// High edges of this histo's axis
+ /// High x edge of this histo's axis
double highEdgeX() const {
return _axis.highEdgeX();
}
+ /// High y edge of this histo's axis
double highEdgeY() const {
return _axis.highEdgeY();
}
@@ -193,25 +200,25 @@
return _axis.getBinIndex(coordX, coordY);
}
- /// Underflow (const version)
- const Dbn2D& underflow() const {
- return _axis.underflow();
- }
-
- /// Return underflow (non-const version)
- Dbn2D& underflow() {
- return _axis.underflow();
- }
-
- /// Return overflow (const version)
- const Dbn2D& overflow() const {
- return _axis.overflow();
- }
-
- /// Return overflow (non-const version)
- Dbn2D& overflow() {
- return _axis.overflow();
- }
+ // /// Underflow (const version)
+ // const Dbn2D& underflow() const {
+ // return _axis.underflow();
+ // }
+
+ // /// Return underflow (non-const version)
+ // Dbn2D& underflow() {
+ // return _axis.underflow();
+ // }
+
+ // /// Return overflow (const version)
+ // const Dbn2D& overflow() const {
+ // return _axis.overflow();
+ // }
+
+ // /// Return overflow (non-const version)
+ // Dbn2D& overflow() {
+ // return _axis.overflow();
+ // }
/// Return a total number of bins in Histo(non-const version)
unsigned int numBinsTotal() {
@@ -224,6 +231,7 @@
//@}
+
public:
/// @name Whole histo data
@@ -235,24 +243,30 @@
return sumW(includeoverflows);
}
- /// Get sum of weights in histo
+ /// Get the sum of weights in histo
double sumW(bool includeoverflows=true) const;
- /// Get sum of squared weights in histo
+ /// Get the sum of squared weights in histo
double sumW2(bool includeoverflows=true) const;
- /// Get the mean
+ /// Get the mean x
double xMean(bool includeoverflows=true) const;
+
+ /// Get the mean y
double yMean(bool includeoverflows=true) const;
- /// Get the variance
+ /// 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
+ /// 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));
}
@@ -292,16 +306,17 @@
//@{
/// @brief Create a Histo1D 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 Need to check that there is a continuous row for this y
/// @todo Change the name!
Histo1D cutterX(double atY, const std::string& path="", const std::string& title="") {
- if(!_axis.isGriddy()) throw GridError("I cannot cut a Histo2D that is not a grid!");
+ if (!_axis.isGriddy()) throw GridError("I cannot cut a Histo2D that is not a grid!");
if (atY < lowEdgeY() || atY > highEdgeY()) throw RangeError("Y is outside the grid");
vector<HistoBin1D> tempBins;
-
+
/// @todo Make all Bin1D constructions happen in loop, to reduce code duplication
const HistoBin2D& first = binByCoord(lowEdgeX(), atY);
const Dbn1D dbn(first.numEntries(), first.sumW(), first.sumW2(), first.sumWX(), first.sumWX2());
@@ -319,17 +334,18 @@
}
- /// Create a Histo1D for the bin slice parallel to the y axis at the specified x coordinate
+ /// @brief Create a Histo1D 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 Need to check that there is a continuous column for this x
/// @todo Change the name!
Histo1D cutterY(double atX, const std::string& path="", const std::string& title="") {
- if(!_axis.isGriddy()) throw GridError("I cannot cut a Histo2D that is not a grid!");
+ if (!_axis.isGriddy()) throw GridError("I cannot cut a Histo2D that is not a grid!");
if (atX < lowEdgeX() || atX > highEdgeX()) throw RangeError("X is outside the grid");
vector<HistoBin1D> tempBins;
-
+
/// @todo Make all Bin1D constructions happen in loop, to reduce code duplication
const HistoBin2D& first = binByCoord(atX, lowEdgeY());
const Dbn1D dbn(first.numEntries(), first.sumW(), first.sumW2(), first.sumWX(), first.sumWX2());
@@ -347,7 +363,7 @@
}
- /// @todo Create x-wise and y-wise conversions to Profile1D
+ /// @todo Create x-wise and y-wise conversions to Profile1D -- ignore outflows for now, but mark as such
//@}
@@ -384,6 +400,12 @@
return tmp;
}
+ /// @todo Multiply histograms?
+
+ /// @brief Divide two histograms
+ ///
+ /// This uses the midpoint instead of the focus
+ /// @todo Set the output to be the focus of the new distribution
Scatter3D operator / (const Histo2D& numer, const Histo2D& denom);
//@}
Modified: trunk/include/YODA/HistoBin1D.h
==============================================================================
--- trunk/include/YODA/HistoBin1D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/HistoBin1D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -37,12 +37,16 @@
: Bin1D(edges)
{ }
+
/// @brief Init a bin with all the components of a fill history.
/// Mainly intended for internal persistency use.
HistoBin1D(double lowedge, double highedge, const Dbn1D& dbn)
: Bin1D(lowedge, highedge, dbn)
{ }
+
+ /// @todo Add copy constructor
+
//@}
Modified: trunk/include/YODA/HistoBin2D.h
==============================================================================
--- trunk/include/YODA/HistoBin2D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/HistoBin2D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -30,6 +30,9 @@
: Bin2D(edges)
{ }
+
+ /// @todo Add copy constructor
+
//@}
Modified: trunk/include/YODA/Point2D.h
==============================================================================
--- trunk/include/YODA/Point2D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/Point2D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -27,7 +27,7 @@
Point2D() { }
- /// Values with optional symmetric errors
+ /// Constructor from values with optional symmetric errors
Point2D(double x, double y, double ex=0.0, double ey=0.0)
: _x(x), _y(y)
{
@@ -36,7 +36,7 @@
}
- /// Values with explicit asymmetric errors
+ /// Constructor from values with explicit asymmetric errors
Point2D(double x, double y,
double exminus,
double explus,
@@ -49,7 +49,7 @@
}
- /// Values with symmetric errors on x and asymmetric errors on y
+ /// Constructor from values with symmetric errors on x and asymmetric errors on y
Point2D(double x, double y, double ex, const std::pair<double,double>& ey)
: _x(x), _y(y), _ey(ey)
{
@@ -57,7 +57,7 @@
}
- /// Values with asymmetric errors on x and symmetric errors on y
+ /// Constructor from values with asymmetric errors on x and symmetric errors on y
Point2D(double x, double y, const std::pair<double,double>& ex, double ey)
: _x(x), _y(y), _ex(ex)
{
@@ -65,11 +65,13 @@
}
- /// Values with asymmetric errors on both x and y
+ /// Constructor from values with asymmetric errors on both x and y
Point2D(double x, double y, const std::pair<double,double>& ex, const std::pair<double,double>& ey)
: _x(x), _y(y), _ex(ex), _ey(ey)
- {
- }
+ { }
+
+
+ /// @todo Add copy constructor
//@}
@@ -92,16 +94,18 @@
void setY(double y) { _y = y; }
//@}
-
+
+
/// Scaling
void scale(double scaleX, double scaleY) {
setX(x()*scaleX);
setY(y()*scaleY);
-
+
setXErr(xErrMinus()*scaleX, xErrPlus()*scaleX);
setYErr(yErrMinus()*scaleY, yErrPlus()*scaleY);
}
+
/// @name x error accessors
//@{
Modified: trunk/include/YODA/Point3D.h
==============================================================================
--- trunk/include/YODA/Point3D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/Point3D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -27,7 +27,7 @@
Point3D() { }
- /// Values with optional symmetric errors
+ /// Constructor from values with optional symmetric errors
Point3D(double x, double y, double z, double ex=0.0, double ey=0.0, double ez=0.0)
: _x(x), _y(y), _z(z)
{
@@ -37,7 +37,7 @@
}
- /// Values with explicit asymmetric errors
+ /// Constructor from values with explicit asymmetric errors
Point3D(double x, double y, double z,
double exminus,
double explus,
@@ -51,10 +51,10 @@
_ey = std::make_pair(eyminus, eyplus);
_ez = std::make_pair(ezminus, ezplus);
}
-
- /// Assymetric Errors given as vectors:
- Point3D(const double& x,
- const double& y,
+
+ /// Constructor from asymmetric errors given as vectors
+ Point3D(const double& x,
+ const double& y,
const double& z,
const std::pair<double,double>& ex,
const std::pair<double,double>& ey,
@@ -63,6 +63,8 @@
}
+ /// @todo Add copy constructor
+
//@}
@@ -281,7 +283,7 @@
/// Equality operator
inline bool operator==(const Point3D& a, const YODA::Point3D& b) {
const bool same_val = fuzzyEquals(a.x(), b.x()) && fuzzyEquals(a.y(), b.y());
- const bool same_eminus = fuzzyEquals(a.xErrMinus(), b.xErrMinus()) &&
+ const bool same_eminus = fuzzyEquals(a.xErrMinus(), b.xErrMinus()) &&
fuzzyEquals(a.yErrMinus(), b.yErrMinus());
const bool same_eplus = fuzzyEquals(a.xErrPlus(), b.xErrPlus()) &&
fuzzyEquals(a.yErrPlus(), b.yErrPlus());
Modified: trunk/include/YODA/ProfileBin1D.h
==============================================================================
--- trunk/include/YODA/ProfileBin1D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/ProfileBin1D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -48,6 +48,8 @@
_ydbn(dbny)
{ }
+ /// @todo Add copy constructor
+
//@}
Modified: trunk/include/YODA/Scatter2D.h
==============================================================================
--- trunk/include/YODA/Scatter2D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/Scatter2D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -32,20 +32,24 @@
/// @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)
{ }
+
/// @todo Add constructor from generic container/Range
- /// Values with no errors
+
+ /// Constructor from a vector of values with no errors
Scatter2D(const std::vector<double>& x, const std::vector<double>& y,
const std::string& path="", const std::string& title="")
: AnalysisObject("Scatter2D", path, title)
@@ -56,7 +60,8 @@
}
}
- /// Values with symmetric errors on x and y
+
+ /// Constructor from vectors of values with symmetric errors on x and y
Scatter2D(const std::vector<double>& x, const std::vector<double>& y,
const std::vector<double>& ex, const std::vector<double>& ey,
const std::string& path="", const std::string& title="")
@@ -68,7 +73,8 @@
}
}
- /// Values with symmetric errors on x and asymmetric errors on y
+
+ /// Constructor from values with symmetric errors on x and asymmetric errors on y
Scatter2D(const std::vector<double>& x, const std::vector<double>& y,
const std::vector<double>& ex, const std::vector<std::pair<double,double> >& ey,
const std::string& path="", const std::string& title="")
@@ -80,7 +86,8 @@
}
}
- /// Values with asymmetric errors on x and symmetric errors on y
+
+ /// Constructor from values with asymmetric errors on x and symmetric errors on y
Scatter2D(const std::vector<double>& x, const std::vector<double>& y,
const std::vector<std::pair<double,double> >& ex, const std::vector<double>& ey,
const std::string& path="", const std::string& title="")
@@ -92,7 +99,8 @@
}
}
- /// Values with asymmetric errors on both x and y
+
+ /// Constructor from values with asymmetric errors on both x and y
Scatter2D(const std::vector<double>& x, const std::vector<double>& y,
const std::vector<std::pair<double,double> >& ex, const std::vector<std::pair<double,double> >& ey,
const std::string& path="", const std::string& title="")
@@ -104,7 +112,8 @@
}
}
- /// Values with completely explicit asymmetric errors
+
+ /// Constructor from values with completely explicit asymmetric errors
Scatter2D(const std::vector<double>& x, const std::vector<double>& y,
const std::vector<double>& exminus,
const std::vector<double>& explus,
@@ -121,46 +130,56 @@
}
}
+
+ /// @todo Add copy constructor
+
//@}
+ /// @name Modifiers
+ //@{
+
/// Clear all points
void reset() {
_points.clear();
}
-
- /// @name Persistency hooks
- //@{
-
- /// Set the state of the profile object, for unpersisting
- /// @todo Need to set annotations (do that on AO), all-histo Dbns, and dbns for every bin. Delegate!
- // void _setstate() = 0;
+ /// Scaling
+ void scale(double scaleX, double scaleY) {
+ for (size_t i = 0; i < _points.size(); ++i) {
+ _points[i].scale(scaleX, scaleY);
+ }
+ }
//@}
///////////////////////////////////////////////////
+
/// @name Point accessors
//@{
+ /// Number of points in the scatter
size_t numPoints() const {
return _points.size();
}
+ /// Get the collection of points
const Points& points() const {
return _points;
}
+ /// Get a reference to the point with index @a index
Point2D& point(size_t index) {
assert(index < numPoints());
return _points.at(index);
}
+ /// Get the point with index @a index (const version)
const Point2D& point(size_t index) const {
assert(index < numPoints());
return _points.at(index);
@@ -169,45 +188,53 @@
//@}
- /// @name Point adders
+ /// @name Point inserters
//@{
+ /// Insert a new point
Scatter2D& addPoint(const Point2D& pt) {
_points.insert(pt);
return *this;
}
+ /// Insert a new point, defined as the x/y value pair
Scatter2D& addPoint(double x, double y) {
_points.insert(Point2D(x, y));
return *this;
}
+ /// Insert a new point, defined as the x/y value pair and symmetric errors
Scatter2D& addPoint(double x, double y, double ex, double ey) {
_points.insert(Point2D(x, y, ex, ey));
return *this;
}
+ /// Insert a new point, defined as the x/y value pair and mixed errors
Scatter2D& addPoint(double x, double y, std::pair<double,double> ex, double ey) {
_points.insert(Point2D(x, y, ex, ey));
return *this;
}
+ /// Insert a new point, defined as the x/y value pair and mixed errors
Scatter2D& addPoint(double x, double y, double ex, std::pair<double,double> ey) {
_points.insert(Point2D(x, y, ex, ey));
return *this;
}
+ /// Insert a new point, defined as the x/y value pair and asymmetric errors
Scatter2D& addPoint(double x, double y, std::pair<double,double> ex, std::pair<double,double> ey) {
_points.insert(Point2D(x, y, ex, ey));
return *this;
}
+ /// Insert a new point, defined as the x/y value pair and asymmetric errors
Scatter2D& addPoint(double x, double exminus, double explus,
double y, double eyminus, double eyplus) {
_points.insert(Point2D(x, exminus, explus, y, eyminus, eyplus));
return *this;
}
+ /// Insert a collection of new points
Scatter2D& addPoints(Points pts) {
foreach (const Point2D& pt, pts) {
addPoint(pt);
@@ -215,21 +242,18 @@
return *this;
}
- /// Scaling
- void scale(double scaleX, double scaleY) {
- for(unsigned int i=0; i < _points.size(); i++) _points[i].scale(scaleX, scaleY);
- }
-
//@}
+ /// @name Combining sets of scatter points
+ //@{
+
/// @todo Better name?
Scatter2D& combineWith(const Scatter2D& other) {
addPoints(other.points());
return *this;
}
-
/// @todo Better name?
/// @todo Convert to accept a Range or generic
Scatter2D& combineWith(const std::vector<Scatter2D>& others) {
@@ -239,6 +263,8 @@
return *this;
}
+ //@}
+
private:
@@ -249,6 +275,8 @@
};
+ /// @name Combining scatters by merging sets of points
+ //@{
inline Scatter2D combine(const Scatter2D& a, const Scatter2D& b) {
Scatter2D rtn = a;
@@ -266,6 +294,8 @@
return rtn;
}
+ //@}
+
//////////////////////////////////
Modified: trunk/include/YODA/Scatter3D.h
==============================================================================
--- trunk/include/YODA/Scatter3D.h Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/include/YODA/Scatter3D.h Tue Aug 16 15:04:57 2011 (r253)
@@ -14,12 +14,11 @@
#include <set>
#include <string>
#include <utility>
+#include <cassert>
namespace YODA {
-
-
/// A very generic data type which is just a collection of 3D data points with errors
class Scatter3D : public AnalysisObject {
public:
@@ -31,11 +30,13 @@
/// @name Constructors
//@{
+ /// Empty constructor
Scatter3D(const std::string& path="", const std::string& title="")
: AnalysisObject("Scatter3D", path, title)
{ }
+ /// Constructor from a set of points
Scatter3D(const Points& points,
const std::string& path="", const std::string& title="")
: AnalysisObject("Scatter3D", path, title),
@@ -43,7 +44,7 @@
{ }
- /// Values with asymmetric errors on both x and y
+ /// Constructor from vectors of values with asymmetric errors on both x and y
Scatter3D(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z,
const std::vector<std::pair<double,double> >& ex, const std::vector<std::pair<double,double> >& ey, const std::vector<std::pair<double,double> >& ez,
const std::string& path="", const std::string& title="")
@@ -55,7 +56,8 @@
}
}
- /// Values with no errors:
+
+ /// Constructor from vectors of values with no errors
Scatter3D(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z,
const std::string& path="", const std::string& title="") {
std::vector<std::pair<double,double> > null;
@@ -65,7 +67,7 @@
Scatter3D(x, y, z, null, null, null, path, title);
}
- /// Values with completely explicit asymmetric errors
+ /// Constructor from vectors of values with completely explicit asymmetric errors
Scatter3D(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double> z,
const std::vector<double>& exminus,
const std::vector<double>& explus,
@@ -84,11 +86,21 @@
addPoint(Point3D(x[i], exminus[i], explus[i], y[i], eyminus[i], eyplus[i], z[i], ezminus[i], ezplus[i]));
}
}
+
+
/// Copy constructor
- Scatter3D(const Scatter3D&, const std::string& path="");
+ Scatter3D(const Scatter3D& s, const std::string& path="")
+ : AnalysisObject("Scatter3D", (path.size() == 0) ? s.path() : path, s, s.title())
+ {
+ _points = s._points;
+ }
+
//@}
+ /// @name Modifiers
+ //@{
+
/// Clear all points
void reset() {
_points.clear();
@@ -96,21 +108,14 @@
/// Scale
void scale(double scaleX, double scaleY, double scaleZ) {
- for(unsigned int i=0; i < _points.size(); i++) _points[i].scale(scaleX, scaleY, scaleZ);
+ for (size_t i = 0; i < _points.size(); ++i) {
+ _points[i].scale(scaleX, scaleY, scaleZ);
+ }
}
- /// @name Persistency hooks
- //@{
-
- /// Set the state of the profile object, for unpersisting
- /// @todo Need to set annotations (do that on AO), all-histo Dbns, and dbns for every bin. Delegate!
- // void _setstate() = 0;
-
//@}
- ///////////////////////////////////////////////////
-
/// @name Point accessors
//@{
@@ -138,7 +143,7 @@
//@}
- /// @name Point adders
+ /// @name Point inserters
//@{
Scatter3D& addPoint(const Point3D& pt) {
@@ -201,6 +206,9 @@
+ /// @name Combining scatters by merging sets of points
+ //@{
+
inline Scatter3D combine(const Scatter3D& a, const Scatter3D& b) {
Scatter3D rtn = a;
rtn.combineWith(b);
@@ -216,6 +224,8 @@
return rtn;
}
+ //@}
+
//////////////////////////////////
@@ -225,7 +235,10 @@
/// Make a Scatter3D representation of a Histo2D
/// @todo Implement!
- // Scatter3D mkScatter(const Histo2D& h);
+ inline Scatter3D mkScatter(const Histo2D& h) {
+ assert(0 && "Not implemented!");
+ return Scatter3D();
+ }
/// Make a Scatter3D representation of a Profile2D
/// @todo Implement!
Modified: trunk/src/Histo1D.cc
==============================================================================
--- trunk/src/Histo1D.cc Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/src/Histo1D.cc Tue Aug 16 15:04:57 2011 (r253)
@@ -14,7 +14,7 @@
namespace YODA {
- typedef vector<HistoBin1D> Bins;
+ // typedef vector<HistoBin1D> Bins;
void Histo1D::fill(double x, double weight) {
Modified: trunk/src/Histo2D.cc
==============================================================================
--- trunk/src/Histo2D.cc Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/src/Histo2D.cc Tue Aug 16 15:04:57 2011 (r253)
@@ -13,30 +13,29 @@
namespace YODA {
- typedef vector<HistoBin2D> Bins;
+ // typedef vector<HistoBin2D> Bins;
- /// @brief Fill function
- /// It relies on Axis2D for searching procedures and
- /// complies loudly if no bin was found.
int Histo2D::fill(double x, double y, double weight) {
- /// Fill the normal bins
+ // Fill the normal bins
int index = findBinIndex(x, y);
- if(index != -1) {
+ if (index != -1) {
HistoBin2D& b = bin(index);
-
- /// Fill the underflow and overflow nicely
+
+ // Fill the total and outflow distributions
_axis.totalDbn().fill(x, y, weight);
-
+ /// @todo Fill the underflow and overflow nicely
+
b.fill(x, y, weight);
}
- else if (x < _axis.lowEdgeX()) { _axis.underflow().fill(x, y, weight); }
- else if (x >= _axis.highEdgeX()) { _axis.overflow().fill(x, y, weight); }
+ // else if (x < _axis.lowEdgeX()) { _axis.underflow().fill(x, y, weight); }
+ // else if (x >= _axis.highEdgeX()) { _axis.overflow().fill(x, y, weight); }
else throw GridError("You are trying to fill an empty space on a grid!");
return index;
}
+
double Histo2D::sumW(bool includeoverflows) const {
if (includeoverflows) return _axis.totalDbn().sumW();
double sumw = 0;
@@ -61,20 +60,21 @@
if (includeoverflows) return _axis.totalDbn().xMean();
double sumwx = 0;
double sumw = 0;
- for (size_t i = 0; i < bins().size(); i++) {
- sumwx += bins().at(i).sumWX();
- sumw += bins().at(i).sumW();
+ foreach (const Bin2D& b, bins()) {
+ sumwx += b.sumWX();
+ sumw += b.sumW();
}
return sumwx/sumw;
}
+
double Histo2D::yMean(bool includeoverflows) const {
if (includeoverflows) return _axis.totalDbn().yMean();
double sumwy = 0;
double sumw = 0;
- for (size_t i = 0; i < bins().size(); i++) {
- sumwy += bins().at(i).sumWY();
- sumw += bins().at(i).sumW();
+ foreach (const Bin2D& b, bins()) {
+ sumwy += b.sumWY();
+ sumw += b.sumW();
}
return sumwy/sumw;
}
@@ -91,6 +91,7 @@
return sigma2/sumW();
}
+
double Histo2D::yVariance(bool includeoverflows) const {
if (includeoverflows) return _axis.totalDbn().yVariance();
double sigma2 = 0;
@@ -102,21 +103,7 @@
return sigma2/sumW();
}
- ////////////////////////////////////////
- /// Copy constructor with optional new path
- Histo2D::Histo2D(const Histo2D& h, const std::string& path)
- : AnalysisObject("Histo2D", (path.size() == 0) ? h.path() : path, h, h.title())
- {
- _axis = h._axis;
- }
-
-
- ///////////////////////////////////////
-
-
- /// @brief Divide two histograms
- /// This uses the midpoint instead of the focus
Scatter3D operator / (const Histo2D& numer, const Histo2D& denom) {
if(numer != denom) throw GridError("The two Histos are not the same!");
Scatter3D tmp;
@@ -126,7 +113,7 @@
const HistoBin2D& bL = b1 + b2;
assert(fuzzyEquals(b1.midpoint().first, b2.midpoint().first));
assert(fuzzyEquals(b1.midpoint().second, b2.midpoint().second));
-
+
const double x = bL.focus().first;
const double y = bL.focus().second;
@@ -136,7 +123,7 @@
const double eyminus = y - b1.yMin();
const double eyplus = bL.yMax() - y;
//cout << b1.xMin() << " " << b1.xMax() << " " << b1.yMin() << " " << b1.yMax() << " EMinus: " << exminus << " " << eyminus << " Focus: " << x << " " << y << endl;
-
+
const double z = b1.height() / b2.height();
const double ez = z * sqrt( sqr(b1.heightErr()/b1.height()) + sqr(b2.heightErr()/b2.height()) );
tmp.addPoint(x, exminus, explus, y, eyminus, eyplus, z, ez, ez);
@@ -145,4 +132,5 @@
return tmp;
}
+
}
Modified: trunk/src/Scatter3D.cc
==============================================================================
--- trunk/src/Scatter3D.cc Tue Aug 16 12:59:40 2011 (r252)
+++ trunk/src/Scatter3D.cc Tue Aug 16 15:04:57 2011 (r253)
@@ -4,37 +4,29 @@
namespace YODA {
- /// Make a Scatter3D representation of a Histo2D
-/* Scatter3D mkScatter(const Histo2D& h) {
- Scatter3D rtn;
- rtn.setAnnotations(h.annotations());
- rtn.setAnnotation("Type", h.type());
- foreach (const HistoBin2D& b, h.bins()) {
- const double x = b.focus().first;
- const double ex_m = b.focus().first - b.lowEdgeX();
- const double ex_p = b.highEdgeX() - b.focus().first;
-
- const double y = b.focus().second;
- const double ey_m = b.focus().second - b.lowEdgeY();
- const double ey_p = b.highEdgeY() - b.focus().second;
-
- const double z = b.height();
- const double ez = b.heightErr();
- const Point3D pt(x, ex_m, ex_p, y, ey_m, ey_p, z, ez, ez);
- rtn.addPoint(pt);
- }
- //assert(h.numBins() == rtn.numPoints());
- return rtn;
- }
-*/
- ////////////////////////////////////////
+ // /// Make a Scatter3D representation of a Histo2D
+ // Scatter3D mkScatter(const Histo2D& h) {
+ // Scatter3D rtn;
+ // rtn.setAnnotations(h.annotations());
+ // rtn.setAnnotation("Type", h.type());
+ // foreach (const HistoBin2D& b, h.bins()) {
+ // const double x = b.focus().first;
+ // const double ex_m = b.focus().first - b.lowEdgeX();
+ // const double ex_p = b.highEdgeX() - b.focus().first;
+
+ // const double y = b.focus().second;
+ // const double ey_m = b.focus().second - b.lowEdgeY();
+ // const double ey_p = b.highEdgeY() - b.focus().second;
+
+ // const double z = b.height();
+ // const double ez = b.heightErr();
+ // const Point3D pt(x, ex_m, ex_p, y, ey_m, ey_p, z, ez, ez);
+ // rtn.addPoint(pt);
+ // }
+ // //assert(h.numBins() == rtn.numPoints());
+ // return rtn;
+ // }
- /// A copy constructor:
- Scatter3D::Scatter3D(const Scatter3D& s, const std::string& path)
- : AnalysisObject("Scatter3D", (path.size() == 0) ? s.path() : path, s, s.title())
- {
- _points = s._points;
- }
/// Subtract two scatters
Scatter3D operator + (const Scatter3D& first, const Scatter3D& second) {
More information about the yoda-svn
mailing list
|