|
[yoda-svn] r536 - in trunk: . include/YODA srcblackhole at projects.hepforge.org blackhole at projects.hepforge.orgThu Nov 15 16:19:42 GMT 2012
Author: buckley Date: Thu Nov 15 16:19:42 2012 New Revision: 536 Log: Improving division and efficiency treatments, and allowing arbitrary f(x), f(y), and flip transformations on Scatter2D. Modified: trunk/ChangeLog trunk/include/YODA/AnalysisObject.h trunk/include/YODA/Histo1D.h trunk/include/YODA/Point2D.h trunk/include/YODA/Scatter2D.h trunk/src/Histo1D.cc trunk/src/Scatter2D.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Wed Nov 14 19:22:08 2012 (r535) +++ trunk/ChangeLog Thu Nov 15 16:19:42 2012 (r536) @@ -1,3 +1,8 @@ +2012-11-15 Andy Buckley <andy.buckley at cern.ch> + + * Improving division and efficiency treatments, and allowing + arbitrary f(x), f(y), and flip transformations on Scatter2D. + 2012-11-14 Andy Buckley <andy.buckley at cern.ch> * Converting linspace, logspace, and their usage to place the nbins argument first. Modified: trunk/include/YODA/AnalysisObject.h ============================================================================== --- trunk/include/YODA/AnalysisObject.h Wed Nov 14 19:22:08 2012 (r535) +++ trunk/include/YODA/AnalysisObject.h Thu Nov 15 16:19:42 2012 (r536) @@ -14,17 +14,6 @@ namespace YODA { - /// @name Enums - /// - /// Definitions of enums used by many/all of the classes - //@{ - - enum ErrorCombination { QUAD, BINOMIAL }; - - //@} - - - /// AnalysisObject is the base class for histograms and scatters class AnalysisObject { Modified: trunk/include/YODA/Histo1D.h ============================================================================== --- trunk/include/YODA/Histo1D.h Wed Nov 14 19:22:08 2012 (r535) +++ trunk/include/YODA/Histo1D.h Thu Nov 15 16:19:42 2012 (r536) @@ -379,16 +379,20 @@ } - /// Divide two histograms, with a specified error treatment - Scatter2D divide(const Histo1D& numer, const Histo1D& denom, - ErrorCombination erropt=QUAD); - + /// Divide two histograms, with an uncorrelated error treatment + Scatter2D divide(const Histo1D& numer, const Histo1D& denom); /// Divide two histograms, with an uncorrelated error treatment inline Scatter2D operator / (const Histo1D& numer, const Histo1D& denom) { return divide(numer, denom); } + /// @brief Calculate a histogrammed efficiency ratio of two histograms + /// + /// Note that an efficiency is not the same thing as a standard division of two + /// histograms: the errors must be treated as correlated + Scatter2D efficiency(const Histo1D& accepted, const Histo1D& total); + //@} Modified: trunk/include/YODA/Point2D.h ============================================================================== --- trunk/include/YODA/Point2D.h Wed Nov 14 19:22:08 2012 (r535) +++ trunk/include/YODA/Point2D.h Thu Nov 15 16:19:42 2012 (r536) @@ -92,7 +92,7 @@ public: - /// @name Value and error accessors + /// @name Value accessors //@{ /// Get x value @@ -107,6 +107,9 @@ /// Set y value void setY(double y) { _y = y; } + // /// Set x and y values + // void setXY(double x, double y) { setX(x); setY(y); } + //@} @@ -133,21 +136,31 @@ return (_ex.first + _ex.second)/2.0; } + /// Set negative x error + void setXErrMinus(double exminus) { + _ex.first = exminus; + } + + /// Set positive x error + void setXErrPlus(double explus) { + _ex.second = explus; + } + /// Set symmetric x error void setXErr(double ex) { - _ex.first = ex; - _ex.second = ex; + setXErrMinus(ex); + setXErrPlus(ex); } /// Set asymmetric x error - void setXErr(std::pair<double,double> ex) { - _ex = ex; + void setXErr(double exminus, double explus) { + setXErrMinus(exminus); + setXErrPlus(explus); } /// Set asymmetric x error - void setXErr(double exminus, double explus) { - _ex.first = exminus; - _ex.second = explus; + void setXErr(std::pair<double,double> ex) { + _ex = ex; } /// Get value minus negative x-error @@ -160,6 +173,18 @@ return _x + _ex.second; } + /// Set negative x-error via the xMin position + void setXMin(double xmin) { + if (xmin > x()) throw UserError("Attempted to set an xmin > x"); + setXErrMinus(x() - xmin); + } + + /// Set positive x-error via the xMax position + void setXMax(double xmax) { + if (xmax < x()) throw UserError("Attempted to set an xmax < x"); + setXErrPlus(xmax - x()); + } + //@} @@ -186,21 +211,31 @@ return (_ey.first + _ey.second)/2.0; } + /// Set negative y error + void setYErrMinus(double eyminus) { + _ey.first = eyminus; + } + + /// Set positive y error + void setYErrPlus(double eyplus) { + _ey.second = eyplus; + } + /// Set symmetric y error void setYErr(double ey) { - _ey.first = ey; - _ey.second = ey; + setYErrMinus(ey); + setYErrPlus(ey); } /// Set asymmetric y error - void setYErr(std::pair<double,double> ey) { - _ey = ey; + void setYErr(double eyminus, double eyplus) { + setYErrMinus(eyminus); + setYErrPlus(eyplus); } /// Set asymmetric y error - void setYErr(double eyminus, double eyplus) { - _ey.first = eyminus; - _ey.second = eyplus; + void setYErr(std::pair<double,double> ey) { + _ey = ey; } /// Get value minus negative y-error @@ -213,6 +248,18 @@ return _y + _ey.second; } + /// Set negative y-error via the yMin position + void setYMin(double ymin) { + if (ymin > y()) throw UserError("Attempted to set an ymin > y"); + setYErrMinus(y() - ymin); + } + + /// Set positive y-error via the yMax position + void setYMax(double ymax) { + if (ymax < y()) throw UserError("Attempted to set an ymax < y"); + setYErrPlus(ymax - y()); + } + //@} Modified: trunk/include/YODA/Scatter2D.h ============================================================================== --- trunk/include/YODA/Scatter2D.h Wed Nov 14 19:22:08 2012 (r535) +++ trunk/include/YODA/Scatter2D.h Thu Nov 15 16:19:42 2012 (r536) @@ -240,7 +240,8 @@ } /// Insert a new point, defined as the x/y value pair and asymmetric errors - Scatter2D& addPoint(double x, double y, double exminus, double explus, + Scatter2D& addPoint(double x, double y, + double exminus, double explus, double eyminus, double eyplus) { _points.insert(Point2D(x, y, exminus, explus, eyminus, eyplus)); return *this; @@ -360,6 +361,75 @@ return divide(numer, denom); } + + + //@} + + + /// @name Transforming operations on Scatter2D + //@{ + + /// @brief Apply transformation fx(x) to all values and error positions (operates in-place on @a s) + /// + /// fx should be a function which takes double x -> double newx + template<typename FNX> + inline void transformX(Scatter2D& s, FNX fx) { + for (size_t i = 0; i < s.numPoints(); ++i) { + Point2D& p = s.point(i); + const double newx = fx(p.x()); + const double fx_xmin = fx(p.xMin()); + const double fx_xmax = fx(p.xMax()); + // Deal with possible inversions of min/max ordering under the transformation + const double newxmin = std::min(fx_xmin, fx_xmax); + const double newxmax = std::max(fx_xmin, fx_xmax); + // Set new point x values + p.setX(newx); + p.setXMin(newxmin); + p.setXMax(newxmax); + } + } + + + /// @brief Apply transformation fy(y) to all values and error positions (operates in-place on @a s) + /// + /// fy should be a function which takes double y -> double newy + template<typename FNY> + inline void transformY(Scatter2D& s, FNY fy) { + for (size_t i = 0; i < s.numPoints(); ++i) { + Point2D& p = s.point(i); + const double newy = fy(p.y()); + const double fy_ymin = fy(p.yMin()); + const double fy_ymax = fy(p.yMax()); + // Deal with possible inversions of min/max ordering under the transformation + const double newymin = std::min(fy_ymin, fy_ymax); + const double newymax = std::max(fy_ymin, fy_ymax); + // Set new point y values + p.setY(newy); + p.setYMin(newymin); + p.setYMax(newymax); + } + } + + + /// @brief Exchange the x and y axes (operates in-place on @a s) + inline void flip(Scatter2D& s) { + for (size_t i = 0; i < s.numPoints(); ++i) { + Point2D& p = s.point(i); + const double newx = p.y(); + const double newy = p.x(); + const double newxmin = p.yMin(); + const double newxmax = p.yMax(); + const double newymin = p.xMin(); + const double newymax = p.xMax(); + p.setX(newx); + p.setY(newy); + p.setXMin(newxmin); + p.setXMax(newxmax); + p.setYMin(newymin); + p.setYMax(newymax); + } + } + //@} Modified: trunk/src/Histo1D.cc ============================================================================== --- trunk/src/Histo1D.cc Wed Nov 14 19:22:08 2012 (r535) +++ trunk/src/Histo1D.cc Thu Nov 15 16:19:42 2012 (r536) @@ -129,52 +129,35 @@ // Divide two histograms - Scatter2D divide(const Histo1D& numer, const Histo1D& denom, - ErrorCombination erropt) { - //if (!numer._axis.compatible(denom._axis)) ... - - Scatter2D tmp; - for (size_t i = 0; i < numer.numBins(); ++i) { - const HistoBin1D& b1 = numer.bin(i); - const HistoBin1D& b2 = denom.bin(i); - - /// @todo Can't we do this check with a method on Axis1D? - if (!fuzzyEquals(b1.midpoint(), b2.midpoint())) { - throw BinningError("Axis binnings are not equivalent"); - } + Scatter2D divide(const Histo1D& numer, const Histo1D& denom) { + return divide(mkScatter(numer), mkScatter(denom)); + } - const HistoBin1D bL = b1 + b2; - const double x = bL.focus(); - /// @todo Use exceptions instead, since this is a user-inducible error - assert(fuzzyEquals(b1.xMin(), b2.xMin())); - assert(fuzzyEquals(b1.xMax(), b2.xMax())); - const double exminus = x - b1.xMin(); - const double explus = b1.xMax() - x; - /// @todo Use exceptions instead, since this is a user-inducible error - assert(exminus >= 0); - assert(explus >= 0); - // - if (b2.height() == 0) continue; // Don't create bins with inf or nan - const double y = b1.height() / b2.height(); - double ey = -1; - - switch (erropt) { - case QUAD: - if (b1.height() == 0) continue; // Don't create bins with inf or nan - ey = y * sqrt( sqr(b1.heightErr()/b1.height()) + sqr(b2.heightErr()/b2.height()) ); - break; - case BINOMIAL: - /// @todo Check that this is correct -- isn't it a problem that this varies if the same scale - /// factor is applied to the weights on bins 1 and 2? - /// @todo I think this is only valid if the fills of b1 are a subset of the fills of b2. Throw an - /// error if Neff(b1) > Neff(b2) - if (b2.effNumEntries() == 0) continue; // Don't create bins with inf or nan - ey = sqrt( b1.effNumEntries() * (1 - b1.effNumEntries()/b2.effNumEntries()) ) / b2.effNumEntries(); - break; - } - tmp.addPoint(x, y, exminus, explus, ey, ey); + // Calculate a histogrammed efficiency ratio of two histograms + Scatter2D efficiency(const Histo1D& accepted, const Histo1D& total) { + Scatter2D tmp = divide(accepted, total); + for (size_t i = 0; i < accepted.numBins(); ++i) { + const HistoBin1D& b_acc = accepted.bin(i); + const HistoBin1D& b_tot = total.bin(i); + Point2D& point = tmp.point(i); + + // Check that the numerator is a subset of the denominator + if (b_acc.effNumEntries() > b_tot.effNumEntries() || b_acc.sumW() > b_tot.sumW()) + throw UserError("Attempt to calculate an efficiency when the numerator is not a subset of the denominator"); + + // If no entries on the denominator, set eff = 0 and move to the next bin + if (b_tot.effNumEntries() == 0) { + point.setY(0.0); + point.setYErr(0.0, 0.0); + continue; + } + // Calculate the values and errors, using numEntries rather than sumW + const double eff = b_acc.effNumEntries() / b_tot.effNumEntries(); + const double ey = sqrt( b_acc.effNumEntries() * (1 - b_acc.effNumEntries()/b_tot.effNumEntries()) ) / b_tot.effNumEntries(); + point.setY(eff); + point.setYErr(ey, ey); } return tmp; } Modified: trunk/src/Scatter2D.cc ============================================================================== --- trunk/src/Scatter2D.cc Wed Nov 14 19:22:08 2012 (r535) +++ trunk/src/Scatter2D.cc Thu Nov 15 16:19:42 2012 (r536) @@ -66,8 +66,8 @@ for (size_t i = 0; i < first.numPoints(); ++i) { const Point2D& p1 = first.point(i); const Point2D& p2 = second.point(i); - assert(fuzzyEquals(p1.xMin(), p2.xMin())); - assert(fuzzyEquals(p1.xMax(), p2.xMax())); + if (!fuzzyEquals(p1.xMin(), p2.xMin()) || !fuzzyEquals(p1.xMax(), p2.xMax())) + throw BinningError("Point x 'bins' are not equivalent"); const double x = (p1.x() + p2.x())/2.0; // const double y = p1.y() + p2.y(); @@ -86,8 +86,8 @@ for (size_t i = 0; i < first.numPoints(); ++i) { const Point2D& p1 = first.point(i); const Point2D& p2 = second.point(i); - assert(fuzzyEquals(p1.xMin(), p2.xMin())); - assert(fuzzyEquals(p1.xMax(), p2.xMax())); + if (!fuzzyEquals(p1.xMin(), p2.xMin()) || !fuzzyEquals(p1.xMax(), p2.xMax())) + throw BinningError("Point x 'bins' are not equivalent"); const double x = (p1.x() + p2.x())/2.0; // const double y = p1.y() - p2.y(); @@ -106,14 +106,16 @@ for (size_t i = 0; i < numer.numPoints(); ++i) { const Point2D& p1 = numer.point(i); const Point2D& p2 = denom.point(i); - assert(fuzzyEquals(p1.xMin(), p2.xMin())); - assert(fuzzyEquals(p1.xMax(), p2.xMax())); + if (!fuzzyEquals(p1.xMin(), p2.xMin()) || !fuzzyEquals(p1.xMax(), p2.xMax())) + throw BinningError("Point x 'bins' are not equivalent"); const double x = (p1.x() + p2.x())/2.0; // const double y = p1.y() / p2.y(); - /// @todo Deal with +/- errors separately - const double ey = y * sqrt( sqr(p1.yErrAvg()/p1.y()) + sqr(p2.yErrAvg()/p2.y()) ); - tmp.addPoint(x, p1.xErrMinus(), p1.xErrPlus(), y, ey, ey); + /// Deal with +/- errors separately, inverted for the denominator contributions + // const double eyavg = y * sqrt( sqr(p1.yErrAvg()/p1.y()) + sqr(p2.yErrAvg()/p2.y()) ); + const double eyplus = y * sqrt( sqr(p1.yErrPlus()/p1.y()) + sqr(p2.yErrMinus()/p2.y()) ); + const double eyminus = y * sqrt( sqr(p1.yErrMinus()/p1.y()) + sqr(p2.yErrPlus()/p2.y()) ); + tmp.addPoint(x, p1.xErrMinus(), p1.xErrPlus(), y, eyminus, eyplus); } assert(tmp.numPoints() == numer.numPoints()); return tmp;
More information about the yoda-svn mailing list |