diff --git a/include/YODA/Counter.h b/include/YODA/Counter.h
--- a/include/YODA/Counter.h
+++ b/include/YODA/Counter.h
@@ -1,279 +1,279 @@
// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2016 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_Counter_h
#define YODA_Counter_h
#include "YODA/AnalysisObject.h"
#include "YODA/Dbn0D.h"
#include "YODA/Scatter1D.h"
#include "YODA/Exceptions.h"
#include <vector>
#include <string>
#include <map>
namespace YODA {
/// A weighted counter.
class Counter : public AnalysisObject {
public:
/// @name Constructors
//@{
/// Default constructor
Counter(const std::string& path="", const std::string& title="")
: AnalysisObject("Counter", path, title)
{ }
/// @brief Constructor accepting an explicit Dbn0D.
///
/// Intended both for internal persistency and user use.
Counter(const Dbn0D& dbn,
const std::string& path="", const std::string& title="")
: AnalysisObject("Counter", path, title),
_dbn(dbn)
{ }
/// @brief Constructor accepting a double (treated as the weight of a single fill).
///
/// Intended for user convenience only, so Counter can be treated as a number.
Counter(double w,
const std::string& path="", const std::string& title="")
: AnalysisObject("Counter", path, title)
{
_dbn.fill(w);
}
/// Copy constructor with optional new path
/// @todo Don't copy the path?
Counter(const Counter& c, const std::string& path="");
/// Assignment operator
Counter& operator = (const Counter& c) {
setPath(c.path());
setTitle(c.title());
_dbn = c._dbn;
return *this;
}
/// Make a copy on the stack
Counter clone() const {
return Counter(*this);
}
/// Make a copy on the heap, via 'new'
Counter* newclone() const {
return new Counter(*this);
}
//@}
/// Fill dimension of this data object
size_t dim() const { return 0; }
/// @name Modifiers
//@{
/// Fill histo by value and weight
- void fill(double weight=1.0, double fraction=1.0) {
+ virtual void fill(double weight=1.0, double fraction=1.0) {
_dbn.fill(weight, fraction);
}
/// @brief Reset the histogram.
///
/// Keep the binning but set all bin contents and related quantities to zero
virtual void reset() {
_dbn.reset();
}
/// Rescale as if all fill weights had been different by factor @a scalefactor.
void scaleW(double scalefactor) {
setAnnotation("ScaledBy", annotation<double>("ScaledBy", 1.0) * scalefactor);
_dbn.scaleW(scalefactor);
}
//@}
/// @name Data access
//@{
/// Get the number of fills
unsigned long numEntries() const { return _dbn.numEntries(); }
/// Get the effective number of fills
double effNumEntries() const { return _dbn.effNumEntries(); }
/// Get the sum of weights
double sumW() const { return _dbn.sumW(); }
/// Get the sum of squared weights
double sumW2() const { return _dbn.sumW2(); }
/// Get the value
double val() const { return sumW(); }
/// Get the uncertainty on the value
/// @todo Implement on Dbn0D and feed through to this and Dbn1D, 2D, etc.
// double err() const { return _dbn.err(); }
double err() const { return sqrt(sumW2()); }
/// Get the relative uncertainty on the value
/// @todo Implement on Dbn0D and feed through to this and Dbn1D, 2D, etc.
// double err() const { return _dbn.err(); }
double relErr() const {
/// @todo Throw excp if sumW2 is 0?
return sumW2() != 0 ? err()/sumW() : 0;
}
//@}
/// @name Internal state access and modification (mainly for persistency use)
//@{
/// Get the internal distribution object
const Dbn0D& dbn() const {
return _dbn;
}
/// Set the internal distribution object: CAREFUL!
void setDbn(const Dbn0D& dbn) {
_dbn = dbn;
}
// /// Set the whole object state
// void setState(const Dbn0D& dbn, const AnalysisObject::Annotations& anns=AnalysisObject::Annotations()) {
// setDbn(dbn);
// setAnnotations(anns);
// }
//@}
/// @name Adding and subtracting counters
//@{
/// Add another counter to this
Counter& operator += (const Counter& toAdd) {
_dbn += toAdd._dbn;
return *this;
}
/// Subtract another counter from this
Counter& operator -= (const Counter& toSubtract) {
_dbn -= toSubtract._dbn;
return *this;
}
/// Increment as if by a fill of weight = 1
/// @note This is post-increment only, i.e. cn++ not ++cn
Counter& operator ++ () {
*this += 1;
return *this;
}
/// Increment as if by a fill of weight = -1
/// @note This is post-decrement only, i.e. cn-- not --cn
Counter& operator -- () {
*this -= 1;
return *this;
}
/// Scale by a double (syntactic sugar for @c scaleW(s))
Counter& operator *= (double s) {
scaleW(s);
return *this;
}
/// Inverse-scale by a double (syntactic sugar for @c scaleW(1/s))
Counter& operator /= (double s) {
scaleW(1/s);
return *this;
}
//@}
private:
/// @name Data
//@{
/// Contained 0D distribution
Dbn0D _dbn;
//@}
};
/// @name Combining counters: global operators
//@{
/// Add two counters
inline Counter add(const Counter& first, const Counter& second) {
Counter tmp = first;
tmp += second;
return tmp;
}
/// Add two counters
inline Counter operator + (const Counter& first, const Counter& second) {
return add(first, second);
}
/// Subtract two counters
inline Counter subtract(const Counter& first, const Counter& second) {
Counter tmp = first;
tmp -= second;
return tmp;
}
/// Subtract two counters
inline Counter operator - (const Counter& first, const Counter& second) {
return subtract(first, second);
}
/// Divide two counters, with an uncorrelated error treatment
/// @todo Or just return a Point1D?
Scatter1D divide(const Counter& numer, const Counter& denom);
/// Divide two counters, with an uncorrelated error treatment
/// @todo Or just return a Point1D?
inline Scatter1D operator / (const Counter& numer, const Counter& denom) {
return divide(numer, denom);
}
/// @todo Add divide functions/operators on pointers
/// @brief Calculate an efficiency ratio of two counters
///
/// Note that an efficiency is not the same thing as a standard division of two
/// histograms: the errors must be treated as correlated
///
/// @todo Or just return a Point1D?
Scatter1D efficiency(const Counter& accepted, const Counter& total);
//@}
}
#endif
diff --git a/include/YODA/Histo1D.h b/include/YODA/Histo1D.h
--- a/include/YODA/Histo1D.h
+++ b/include/YODA/Histo1D.h
@@ -1,519 +1,519 @@
// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2016 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_Histo1D_h
#define YODA_Histo1D_h
#include "YODA/AnalysisObject.h"
#include "YODA/HistoBin1D.h"
#include "YODA/Dbn1D.h"
#include "YODA/Scatter2D.h"
#include "YODA/Axis1D.h"
#include "YODA/Exceptions.h"
#include <vector>
#include <string>
#include <map>
namespace YODA {
/// Convenience typedef
typedef Axis1D<HistoBin1D, Dbn1D> Histo1DAxis;
/// A one-dimensional histogram.
class Histo1D : public AnalysisObject {
public:
/// Convenience typedefs
typedef Histo1DAxis Axis;
typedef Axis::Bins Bins;
typedef HistoBin1D Bin;
/// @name Constructors
//@{
/// Default constructor
Histo1D(const std::string& path="", const std::string& title="")
: AnalysisObject("Histo1D", path, title),
_axis()
{ }
/// Constructor giving range and number of bins.
Histo1D(size_t nbins, double lower, double upper,
const std::string& path="", const std::string& title="")
: AnalysisObject("Histo1D", path, title),
_axis(nbins, lower, upper)
{ }
/// @brief Constructor giving explicit bin edges.
///
/// For n bins, binedges.size() == n+1, the last one being the upper bound
/// of the last bin
Histo1D(const std::vector<double>& binedges,
const std::string& path="", const std::string& title="")
: AnalysisObject("Histo1D", path, title),
_axis(binedges)
{ }
/// Constructor accepting an explicit collection of bins.
Histo1D(const std::vector<Bin>& bins,
const std::string& path="", const std::string& title="")
: AnalysisObject("Histo1D", path, title),
_axis(bins)
{ }
/// Copy constructor with optional new path
/// @todo Also allow title setting from the constructor?
Histo1D(const Histo1D& h, const std::string& path="");
/// Constructor from a Scatter2D's binning, with optional new path
/// @todo Also allow title setting from the constructor?
Histo1D(const Scatter2D& s, const std::string& path="");
/// Constructor from a Profile1D's binning, with optional new path
/// @todo Also allow title setting from the constructor?
Histo1D(const Profile1D& p, const std::string& path="");
/// @brief State-setting constructor
///
/// Intended principally for internal persistency use.
Histo1D(const std::vector<HistoBin1D>& bins,
const Dbn1D& dbn_tot, const Dbn1D& dbn_uflow, const Dbn1D& dbn_oflow,
const std::string& path="", const std::string& title="")
: AnalysisObject("Histo1D", path, title),
_axis(bins, dbn_tot, dbn_uflow, dbn_oflow)
{ }
/// Assignment operator
Histo1D& operator = (const Histo1D& h1) {
AnalysisObject::operator = (h1); //< AO treatment of paths etc.
_axis = h1._axis;
return *this;
}
/// Make a copy on the stack
Histo1D clone() const {
return Histo1D(*this);
}
/// Make a copy on the heap, via 'new'
Histo1D* newclone() const {
return new Histo1D(*this);
}
//@}
/// Fill dimension of this data object
size_t dim() const { return 1; }
/// @name Modifiers
//@{
/// @brief Reset the histogram.
///
/// Keep the binning but set all bin contents and related quantities to zero
virtual void reset() {
_axis.reset();
}
/// Fill histo by value and weight, optionally as a fractional fill
virtual void fill(double x, double weight=1.0, double fraction=1.0);
/// Fill histo bin i with the given weight, optionally as a fractional fill
- void fillBin(size_t i, double weight=1.0, double fraction=1.0);
+ virtual void fillBin(size_t i, double weight=1.0, double fraction=1.0);
/// Rescale as if all fill weights had been different by factor @a scalefactor.
void scaleW(double scalefactor) {
setAnnotation("ScaledBy", annotation<double>("ScaledBy", 1.0) * scalefactor);
_axis.scaleW(scalefactor);
}
/// Normalize the (visible) histo area to the @a normto value.
///
/// If @a includeoverflows is true, the original normalisation is computed with
/// the overflow bins included, so that the resulting visible normalisation can
/// be less than @a normto. This is probably what you want.
void normalize(double normto=1.0, bool includeoverflows=true) {
const double oldintegral = integral(includeoverflows);
if (oldintegral == 0) throw WeightError("Attempted to normalize a histogram with null area");
/// @todo Check that this is the desired behaviour
scaleW(normto / oldintegral);
}
/// Merge together the bin range with indices from @a from to @a to, inclusive
void mergeBins(size_t from, size_t to) {
_axis.mergeBins(from, to);
}
/// Merge every group of n bins, starting from the LHS
void rebinBy(unsigned int n, size_t begin=0, size_t end=UINT_MAX) {
_axis.rebinBy(n, begin, end);
}
/// Overloaded alias for rebinBy
void rebin(unsigned int n, size_t begin=0, size_t end=UINT_MAX) {
rebinBy(n, begin, end);
}
/// Rebin to the given list of bin edges
void rebinTo(const std::vector<double>& newedges) {
_axis.rebinTo(newedges);
}
/// Overloaded alias for rebinTo
void rebin(const std::vector<double>& newedges) {
rebinTo(newedges);
}
//@}
public:
/// @name Bin accessors
//@{
/// Number of bins on this axis (not counting under/overflow)
size_t numBins() const { return bins().size(); }
/// Low edge of this histo's axis
double xMin() const { return _axis.xMin(); }
/// High edge of this histo's axis
double xMax() const { return _axis.xMax(); }
/// All bin edges on this histo's axis
///
/// @note This only returns the finite edges, i.e. -inf and +inf are removed
/// @todo Make the +-inf stripping controllable by a default-valued bool arg
const std::vector<double> xEdges() const { return _axis.xEdges(); }
/// Access the bin vector
std::vector<YODA::HistoBin1D>& bins() { return _axis.bins(); }
/// Access the bin vector (const version)
const std::vector<YODA::HistoBin1D>& bins() const { return _axis.bins(); }
/// Access a bin by index (non-const version)
HistoBin1D& bin(size_t index) { return _axis.bins()[index]; }
/// Access a bin by index (const version)
const HistoBin1D& bin(size_t index) const { return _axis.bins()[index]; }
/// Access a bin index by coordinate
/// @todo Convert to ssize_t?
int binIndexAt(double x) {
return _axis.binIndexAt(x);
}
/// Access a bin by coordinate (const version)
const HistoBin1D& binAt(double x) const { return _axis.binAt(x); }
/// Access summary distribution, including gaps and overflows (non-const version)
Dbn1D& totalDbn() { return _axis.totalDbn(); }
/// Access summary distribution, including gaps and overflows (const version)
const Dbn1D& totalDbn() const { return _axis.totalDbn(); }
/// Set summary distribution, mainly for persistency: CAREFUL!
void setTotalDbn(const Dbn1D& dbn) { _axis.setTotalDbn(dbn); }
/// Access underflow (non-const version)
Dbn1D& underflow() { return _axis.underflow(); }
/// Access underflow (const version)
const Dbn1D& underflow() const { return _axis.underflow(); }
/// Set underflow distribution, mainly for persistency: CAREFUL!
void setUnderflow(const Dbn1D& dbn) { _axis.setUnderflow(dbn); }
/// Access overflow (non-const version)
Dbn1D& overflow() { return _axis.overflow(); }
/// Access overflow (const version)
const Dbn1D& overflow() const { return _axis.overflow(); }
/// Set overflow distribution, mainly for persistency: CAREFUL!
void setOverflow(const Dbn1D& dbn) { _axis.setOverflow(dbn); }
/// Add a new bin specifying its lower and upper bound
void addBin(double from, double to) { _axis.addBin(from, to); }
/// Add new bins by specifying a vector of edges
void addBins(std::vector<double> edges) { _axis.addBins(edges); }
// /// Add new bins specifying a beginning and end of each of them
// void addBins(std::vector<std::pair<double,double> > edges) {
// _axis.addBins(edges);
// }
/// Add a new bin, perhaps already populated: CAREFUL!
void addBin(const HistoBin1D& b) { _axis.addBin(b); }
/// @brief Bins addition operator
///
/// Add multiple bins without resetting
void addBins(const Bins& bins) {
_axis.addBins(bins);
}
/// Remove a bin
void eraseBin(size_t index) { _axis.eraseBin(index); }
//@}
/// @name Whole histo data
//@{
/// Get the total area (sumW) of the histogram
double integral(bool includeoverflows=true) const { return sumW(includeoverflows); }
/// @brief Get the integrated area of the histogram between bins @a binindex1 and @a binindex2.
///
/// @note The area of bin @a binindex2 _is_ included in the returned
/// value. To include the underflow and overflow areas, you should add them
/// explicitly with the underflow() and overflow() methods.
///
/// @todo Allow int bin index args for type compatibility with binIndexAt()?
double integralRange(size_t binindex1, size_t binindex2) const {
assert(binindex2 >= binindex1);
if (binindex1 >= numBins()) throw RangeError("binindex1 is out of range");
if (binindex2 >= numBins()) throw RangeError("binindex2 is out of range");
double rtn = 0;
for (size_t i = binindex1; i < binindex2; ++i) {
rtn += bin(i).sumW();
}
return rtn;
}
/// @brief Get the integrated area of the histogram up to bin @a binindex.
///
/// @note The area of bin @a binindex _is_ included in the returned
/// value. To not include the underflow, set includeunderflow=false.
///
/// @todo Allow int bin index args for type compatibility with binIndexAt()?
double integralTo(size_t binindex, bool includeunderflow=true) const {
double rtn = includeunderflow ? underflow().sumW() : 0;
rtn += integralRange(0, binindex);
return rtn;
}
/// Get the number of fills
unsigned long numEntries(bool includeoverflows=true) const;
/// Get the effective number of fills
double effNumEntries(bool includeoverflows=true) const;
/// Get sum of weights in histo
double sumW(bool includeoverflows=true) const;
/// Get sum of squared weights in histo
double sumW2(bool includeoverflows=true) const;
/// Get the mean in x
double xMean(bool includeoverflows=true) const;
/// Get the variance in x
double xVariance(bool includeoverflows=true) const;
/// Get the standard deviation in x
double xStdDev(bool includeoverflows=true) const {
if (includeoverflows) return _axis.totalDbn().xStdDev();
return std::sqrt(xVariance(includeoverflows));
}
/// Get the standard error in x
double xStdErr(bool includeoverflows=true) const;
/// Get the RMS in x
double xRMS(bool includeoverflows=true) const;
//@}
/// @name Adding and subtracting histograms
//@{
/// @brief Add another histogram to this one
///
/// @note Adding histograms will unset any ScaledBy attribute from prevous calls to scaleW or normalize.
Histo1D& operator += (const Histo1D& toAdd) {
if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
_axis += toAdd._axis;
return *this;
// if (!hasAnnotation("ScaledBy") && !toAdd.hasAnnotation("ScaledBy")) {
// _axis += toAdd._axis;
// } else {
// // Undo scaling of both histograms
// double scaledBy = annotation<double>("ScaledBy", 1.0);
// _axis.scaleW(1.0/scaledBy);
// double toAddScaledBy = toAdd.annotation<double>("ScaledBy", 1.0);
// Axis1D<HistoBin1D, Dbn1D> toAddAxis = toAdd._axis;
// toAddAxis.scaleW(1.0/toAddScaledBy);
// _axis += toAddAxis;
// // Re-apply combined scaling
// double newScaledBy = scaledBy*toAddScaledBy/(scaledBy+toAddScaledBy);
// _axis.scaleW(newScaledBy);
// setAnnotation("ScaledBy", newScaledBy);
// }
/// @todo What about if one histo sets ScaledBy, and the other doesn't?!? Aaaargh
// return *this;
}
/// @brief Subtract another histogram from this one
///
/// @note Subtracting histograms will unset any ScaledBy attribute from prevous calls to scaleW or normalize.
Histo1D& operator -= (const Histo1D& toSubtract) {
if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
_axis -= toSubtract._axis;
return *this;
}
//@}
protected:
/// Access a bin by coordinate (non-const version)
HistoBin1D& binAt(double x) { return _axis.binAt(x); }
private:
/// @name Bin data
//@{
/// Definition of bin edges and contents
Axis1D<HistoBin1D, Dbn1D> _axis;
//@}
};
/// @name Combining histos: global operators
//@{
/// Add two histograms
inline Histo1D add(const Histo1D& first, const Histo1D& second) {
Histo1D tmp = first;
if (first.path() != second.path()) tmp.setPath("");
tmp += second;
return tmp;
}
/// Add two histograms
inline Histo1D operator + (const Histo1D& first, const Histo1D& second) {
return add(first, second);
}
/// Subtract two histograms
inline Histo1D subtract(const Histo1D& first, const Histo1D& second) {
Histo1D tmp = first;
if (first.path() != second.path()) tmp.setPath("");
tmp -= second;
return tmp;
}
/// Subtract two histograms
inline Histo1D operator - (const Histo1D& first, const Histo1D& second) {
return subtract(first, second);
}
/// @todo Add multiply(H1, H1) -> Scatter2D?
/// Divide two histograms, with an uncorrelated error treatment
///
/// @todo Wouldn't it be nice to be able to supply a correlation matrix or function as optional arg?
///
/// @note The two histos must have _exactly_ the same binning.
Scatter2D divide(const Histo1D& numer, const Histo1D& denom);
/// Divide two histograms, with an uncorrelated error treatment
///
/// @note The two histos must have _exactly_ the same binning.
inline Scatter2D operator / (const Histo1D& numer, const Histo1D& denom) {
return divide(numer, denom);
}
/// @todo Add functions/operators on pointers
/// @brief Calculate a histogrammed efficiency ratio of two histograms
///
/// @note The two histos must have _exactly_ the same binning.
///
/// @note An efficiency is not the same thing as a standard division of two
/// histograms: the errors are treated as correlated via binomial statistics.
Scatter2D efficiency(const Histo1D& accepted, const Histo1D& total);
/// @brief Calculate the asymmetry (a-b)/(a+b) of two histograms
///
/// @note The two histos must have _exactly_ the same binning.
inline Scatter2D asymm(const Histo1D& a, const Histo1D& b) {
return (a-b) / (a+b);
}
/// @brief Convert a Histo1D to a Scatter2D representing the integral of the histogram
///
/// @note The integral histo errors are calculated as sqrt(binvalue), as if they
/// are uncorrelated. This is not in general true for integral histograms, so if you
/// need accurate errors you should explicitly monitor bin-to-bin correlations.
///
/// The includeunderflow param chooses whether the underflow bin is included
/// in the integral numbers as an offset.
///
/// @todo Rename/alias as mkIntegral
Scatter2D toIntegralHisto(const Histo1D& h, bool includeunderflow=true);
/// @brief Convert a Histo1D to a Scatter2D where each bin is a fraction of the total
///
/// @note This sounds weird: let's explain a bit more! Sometimes we want to
/// take a histo h, make an integral histogram H from it, and then divide H by
/// the total integral of h, such that every bin in H represents the
/// cumulative efficiency of that bin as a fraction of the total. I.e. an
/// integral histo, scaled by 1/total_integral and with binomial errors.
///
/// The includeunderflow param behaves as for toIntegral, and applies to both
/// the initial integration and the integral used for the scaling. The
/// includeoverflow param applies only to obtaining the scaling factor.
///
/// @todo Rename/alias as mkIntegralEff
Scatter2D toIntegralEfficiencyHisto(const Histo1D& h, bool includeunderflow=true, bool includeoverflow=true);
//@}
}
#endif
diff --git a/include/YODA/Histo2D.h b/include/YODA/Histo2D.h
--- a/include/YODA/Histo2D.h
+++ b/include/YODA/Histo2D.h
@@ -1,512 +1,512 @@
// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2016 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_Histo2D_h
#define YODA_Histo2D_h
#include "YODA/AnalysisObject.h"
#include "YODA/HistoBin2D.h"
#include "YODA/Dbn2D.h"
#include "YODA/Axis2D.h"
#include "YODA/Scatter3D.h"
#include "YODA/Exceptions.h"
#include <vector>
namespace YODA {
// Forward declaration
class Profile2D;
class Scatter3D;
/// Convenience typedef
typedef Axis2D<HistoBin2D, Dbn2D> Histo2DAxis;
/// A two-dimensional histogram.
class Histo2D : public AnalysisObject {
public:
/// Convenience typedefs
typedef Histo2DAxis Axis;
typedef Axis::Bins Bins;
typedef HistoBin2D Bin;
typedef Axis::Outflows Outflows;
/// @name Constructors
//@{
/// Default constructor
Histo2D(const std::string& path="", const std::string& title="")
: AnalysisObject("Histo2D", path, title),
_axis()
{ }
/// Constructor giving range and number of bins.
Histo2D(size_t nbinsX, double lowerX, double upperX,
size_t nbinsY, double lowerY, double upperY,
const std::string& path="", const std::string& title="")
: AnalysisObject("Histo2D", path, title),
_axis(nbinsX, std::make_pair(lowerX, upperX), nbinsY, std::make_pair(lowerY, upperY))
{ }
/// Constructor accepting the bin edges on X and Y axis.
Histo2D(const std::vector<double>& xedges, const std::vector<double>& yedges,
const std::string& path="", const std::string& title="")
: AnalysisObject("Histo2D", path, title),
_axis(xedges, yedges)
{ }
/// Constructor accepting an explicit collection of bins.
Histo2D(const std::vector<Bin>& bins,
const std::string& path="", const std::string& title="")
: AnalysisObject("Histo2D", path, title),
_axis(bins)
{ }
/// Copy constructor with optional new path
/// @todo Also allow title setting from the constructor?
Histo2D(const Histo2D& h, const std::string& path="");
/// A constructor from a Scatter3D's binning, with optional new path
/// @todo Also allow title setting from the constructor?
Histo2D(const Scatter3D& s, const std::string& path="");
/// Constructor from a Profile2D's binning, with optional new path
/// @todo Also allow title setting from the constructor?
Histo2D(const Profile2D& h, const std::string& path="");
/// @brief State-setting constructor
///
/// Mainly intended for internal persistency use.
Histo2D(const std::vector<HistoBin2D>& bins,
const Dbn2D& totalDbn,
const Outflows& outflows,
const std::string& path="", const std::string& title="")
: AnalysisObject("Histo2D", path, title),
_axis(bins, totalDbn, outflows)
{ }
/// Assignment operator
Histo2D& operator = (const Histo2D& h2) {
AnalysisObject::operator = (h2); //< AO treatment of paths etc.
_axis = h2._axis;
return *this;
}
/// Make a copy on the stack
Histo2D clone() const {
return Histo2D(*this);
}
/// Make a copy on the heap, via 'new'
Histo2D* newclone() const {
return new Histo2D(*this);
}
//@}
/// Fill dimension of this data object
size_t dim() const { return 2; }
/// @name Modifiers
//@{
/// Fill histo with weight at (x,y)
- void fill(double x, double y, double weight=1.0, double fraction=1.0);
+ virtual void fill(double x, double y, double weight=1.0, double fraction=1.0);
/// Fill histo x-y bin i with the given weight
- void fillBin(size_t i, double weight=1.0, double fraction=1.0);
+ virtual void fillBin(size_t i, double weight=1.0, double fraction=1.0);
/// @brief Reset the histogram.
///
/// Keep the binning but set all bin contents and related quantities to zero
void reset() {
_axis.reset();
}
/// Rescale as if all fill weights had been different by factor @a scalefactor.
void scaleW(double scalefactor) {
setAnnotation("ScaledBy", annotation<double>("ScaledBy", 1.0) * scalefactor);
_axis.scaleW(scalefactor);
}
/// Normalize the (visible) histo "volume" to the @a normto value.
///
/// If @a includeoverflows is true, the original normalisation is computed with
/// the overflow bins included, so that the resulting visible normalisation can
/// be less than @a normto. This is probably what you want.
void normalize(double normto=1.0, bool includeoverflows=true) {
const double oldintegral = integral(includeoverflows);
if (oldintegral == 0) throw WeightError("Attempted to normalize a histogram with null area");
/// @todo Check that this is the desired behaviour
scaleW(normto / oldintegral);
}
/// Scale the dimensions
void scaleXY(double scaleX = 1.0, double scaleY = 1.0) {
_axis.scaleXY(scaleX, scaleY);
}
/// @brief Bin addition operator
///
/// Add a bin to an axis described by its x and y ranges.
void addBin(Axis::EdgePair1D xrange, Axis::EdgePair1D yrange) {
_axis.addBin(xrange, yrange);
}
/// @brief Bin addition operator
///
/// Add a bin, possibly already populated
void addBin(const Bin& bin) {
_axis.addBin(bin);
}
/// @brief Bins addition operator
///
/// Add multiple bins from edge cuts without resetting
void addBins(const Axis::Edges& xcuts, const Axis::Edges& ycuts) {
_axis.addBins(xcuts, ycuts);
}
/// @brief Bins addition operator
///
/// Add multiple bins without resetting
void addBins(const Bins& bins) {
_axis.addBins(bins);
}
// /// Adding bins
/// @todo TODO
// void addBin(const std::vector<std::pair<std::pair<double,double>, std::pair<double,double> > > coords) {
// _axis.addBin(coords);
// }
// /// Adding bins which is not so eloquent
/// @todo TODO
// void addBin(double lowX, double lowY, double highX, double highY) {
// _axis.addBin(lowX, lowY, highX, highY);
// }
// /// Merge the bins
/// @todo TODO
// void mergeBins(size_t from, size_t to) {
// _axis.mergeBins(from, to);
// }
/// Rebin the whole histo by a @a factorX in the X direction and
/// @a factorY in the Y direction
/// @todo TODO
// void rebin(size_t factorX, size_t factorY){
// _axis.rebin(factorX, factorY);
// }
void eraseBin(size_t index) {
_axis.eraseBin(index);
}
//@}
public:
/// @name Bin accessors
//@{
/// Low x edge of this histo's axis
double xMin() const { return _axis.xMin(); }
/// High x edge of this histo's axis
double xMax() const { return _axis.xMax(); }
/// Low y edge of this histo's axis
double yMin() const { return _axis.yMin(); }
/// High y edge of this histo's axis
double yMax() const { return _axis.yMax(); }
/// Access the bin vector (non-const version)
std::vector<YODA::HistoBin2D>& bins() { return _axis.bins(); }
/// Access the bin vector (const version)
const std::vector<YODA::HistoBin2D>& bins() const { return _axis.bins(); }
/// Access a bin by index (non-const version)
HistoBin2D& bin(size_t index) { return _axis.bin(index); }
/// Access a bin by index (const version)
const HistoBin2D& bin(size_t index) const { return _axis.bin(index); }
/// Access a bin index by coordinate
int binIndexAt(double x, double y) { return _axis.binIndexAt(x, y); }
/// Access a bin by coordinate (const version)
const HistoBin2D& binAt(double x, double y) const { return _axis.binAt(x, y); }
/// Number of bins
size_t numBins() const { return _axis.numBins(); }
/// Number of bins along the x axis
size_t numBinsX() const { return _axis.numBinsX(); }
/// Number of bins along the y axis
size_t numBinsY() const { return _axis.numBinsY(); }
/// Access summary distribution, including gaps and overflows (non-const version)
Dbn2D& totalDbn() { return _axis.totalDbn(); }
/// Access summary distribution, including gaps and overflows (const version)
const Dbn2D& totalDbn() const { return _axis.totalDbn(); }
/// Set summary distribution, including gaps and overflows
void setTotalDbn(const Dbn2D& dbn) { _axis.setTotalDbn(dbn); }
// /// @brief Access an outflow (non-const)
// ///
// /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow.
// /// (0,0) is not a valid overflow index pair, since it is in range for both x and y.
// Dbn2D& outflow(int ix, int iy) {
// std::cout << "Histo2D::outflow\n";
// return _axis.outflow(ix, iy);
// }
// /// @brief Access an outflow (const)
// ///
// /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow.
// /// (0,0) is not a valid overflow index pair, since it is in range for both x and y.
// const Dbn2D& outflow(int ix, int iy) const {
// return _axis.outflow(ix, iy);
// }
//@}
/// @name Whole histo data
//@{
/// Get the total volume of the histogram
double integral(bool includeoverflows=true) const { return sumW(includeoverflows); }
/// Get the number of fills
unsigned long numEntries(bool includeoverflows=true) const;
/// Get the effective number of fills
double effNumEntries(bool includeoverflows=true) const;
/// Get the sum of weights in histo
double sumW(bool includeoverflows=true) const;
/// Get the sum of squared weights in histo
double sumW2(bool includeoverflows=true) const;
/// Get the mean x
double xMean(bool includeoverflows=true) const;
/// Get the mean y
double yMean(bool includeoverflows=true) const;
/// Get the variance in x
double xVariance(bool includeoverflows=true) const;
/// Get the variance in y
double yVariance(bool includeoverflows=true) const;
/// Get the standard deviation in x
double xStdDev(bool includeoverflows=true) const {
return std::sqrt(xVariance(includeoverflows));
}
/// Get the standard deviation in y
double yStdDev(bool includeoverflows=true) const {
return std::sqrt(yVariance(includeoverflows));
}
/// Get the standard error in x
double xStdErr(bool includeoverflows=true) const;
/// Get the standard error in y
double yStdErr(bool includeoverflows=true) const;
/// Get the RMS in x
double xRMS(bool includeoverflows=true) const;
/// Get the RMS in y
double yRMS(bool includeoverflows=true) const;
//@}
/// @name Adding and subtracting histograms
//@{
/// @brief Add another histogram to this one
///
/// @note Adding histograms will unset any ScaledBy attribute from prevous calls to scaleW or normalize.
Histo2D& operator += (const Histo2D& toAdd) {
if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
_axis += toAdd._axis;
return *this;
}
/// @brief Subtract another histogram from this one
///
/// @note Subtracting histograms will unset any ScaledBy attribute from prevous calls to scaleW or normalize.
Histo2D& operator -= (const Histo2D& toSubtract) {
if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
_axis -= toSubtract._axis;
return *this;
}
bool operator == (const Histo2D& other) const {
return _axis == other._axis;
}
bool operator != (const Histo2D& other) const {
return ! operator == (other);
}
//@}
// /// @name Slicing operators
// //@{
// /// @brief Create a Histo2D for the bin slice parallel to the x axis at the specified y coordinate
// ///
// /// Note that the created histogram will not have correctly filled underflow and overflow bins.
// /// @todo It's not really *at* the specified y coord: it's for the corresponding bin row.
// /// @todo Change the name!
// Histo2D cutterX(double atY, const std::string& path="", const std::string& title="");
// /// @brief Create a Histo2D for the bin slice parallel to the y axis at the specified x coordinate
// ///
// /// Note that the created histogram will not have correctly filled underflow and overflow bins.
// /// @todo It's not really *at* the specified x coord: it's for the corresponding bin row.
// /// @todo Change the name!
// Histo2D cutterY(double atX, const std::string& path="", const std::string& title="");
// /// X-wise Profile1D creator from Histo2D
// Profile1D mkProfileX();
// /// Y-wise Profile1D creator from Histo2D
// Profile1D mkProfileY();
// //@}
protected:
/// Access a bin by coordinate (non-const version)
HistoBin2D& binAt(double x, double y) { return _axis.binAt(x, y); }
private:
/// @name Bin data
//@{
/// Definition of bin edges and contents
Axis2D<HistoBin2D, Dbn2D> _axis;
//@}
};
/// @name Combining histos: global operators
//@{
/// Add two histograms
inline Histo2D add(const Histo2D& first, const Histo2D& second) {
Histo2D tmp = first;
if (first.path() != second.path()) tmp.setPath("");
tmp += second;
return tmp;
}
/// Add two histograms
inline Histo2D operator + (const Histo2D& first, const Histo2D& second) {
return add(first, second);
}
/// Subtract two histograms
inline Histo2D subtract(const Histo2D& first, const Histo2D& second) {
Histo2D tmp = first;
if (first.path() != second.path()) tmp.setPath("");
tmp -= second;
return tmp;
}
/// Subtract two histograms
inline Histo2D operator - (const Histo2D& first, const Histo2D& second) {
return subtract(first, second);
}
/// @todo Add multiply(H2, H2) -> Scatter3D?
/// @brief Divide two histograms, with an uncorrelated error treatment
///
/// @todo Wouldn't it be nice to be able to supply a correlation matrix or function as optional arg?
///
/// @note The two histos must have _exactly_ the same binning.
Scatter3D divide(const Histo2D& numer, const Histo2D& denom);
/// Divide two histograms, with an uncorrelated error treatment
///
/// @note The two histos must have _exactly_ the same binning.
inline Scatter3D operator / (const Histo2D& numer, const Histo2D& denom) {
return divide(numer, denom);
}
/// @brief Calculate a histogrammed efficiency ratio of two histograms
///
/// @note The two histos must have _exactly_ the same binning.
///
/// @note An efficiency is not the same thing as a standard division of two
/// histograms: the errors are treated as correlated via binomial statistics.
Scatter3D efficiency(const Histo2D& accepted, const Histo2D& total);
/// @brief Calculate the asymmetry (a-b)/(a+b) of two histograms
///
/// @note The two histos must have _exactly_ the same binning.
inline Scatter3D asymm(const Histo2D& a, const Histo2D& b) {
return (a-b) / (a+b);
}
//@}
}
#endif
diff --git a/include/YODA/Profile1D.h b/include/YODA/Profile1D.h
--- a/include/YODA/Profile1D.h
+++ b/include/YODA/Profile1D.h
@@ -1,402 +1,402 @@
// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2016 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_Profile1D_h
#define YODA_Profile1D_h
#include "YODA/AnalysisObject.h"
#include "YODA/ProfileBin1D.h"
#include "YODA/Scatter2D.h"
#include "YODA/Dbn2D.h"
#include "YODA/Axis1D.h"
#include "YODA/Exceptions.h"
#include <vector>
#include <string>
#include <map>
namespace YODA {
// Forward declarations
class Histo1D;
class Scatter2D;
/// Convenience typedef
typedef Axis1D<ProfileBin1D, Dbn2D> Profile1DAxis;
/// A one-dimensional profile histogram.
class Profile1D : public AnalysisObject {
public:
/// Convenience typedefs
typedef Profile1DAxis Axis;
typedef Axis::Bins Bins;
typedef ProfileBin1D Bin;
/// @name Constructors
//@{
/// Default constructor
Profile1D(const std::string& path="", const std::string& title="")
: AnalysisObject("Profile1D", path, title),
_axis()
{ }
/// Constructor giving range and number of bins
Profile1D(size_t nxbins, double xlower, double xupper,
const std::string& path="", const std::string& title="")
: AnalysisObject("Profile1D", path, title),
_axis(nxbins, xlower, xupper)
{ }
/// Constructor giving explicit bin edges
///
/// For n bins, binedges.size() == n+1, the last one being the upper bound
/// of the last bin
Profile1D(const std::vector<double>& xbinedges,
const std::string& path="", const std::string& title="")
: AnalysisObject("Profile1D", path, title),
_axis(xbinedges)
{ }
/// Copy constructor with optional new path
/// @todo Also allow title setting from the constructor?
Profile1D(const Profile1D& p, const std::string& path="");
/// Constructor from a Scatter2D's binning, with optional new path
/// @todo Also allow title setting from the constructor?
Profile1D(const Scatter2D& s, const std::string& path="");
/// Constructor from a Histo1D's binning, with optional new path
/// @todo Also allow title setting from the constructor?
Profile1D(const Histo1D& h, const std::string& path="");
/// @brief State-setting constructor.
///
/// Intended principally for internal persistency use.
Profile1D(const std::vector<ProfileBin1D>& bins,
const Dbn2D& dbn_tot, const Dbn2D& dbn_uflow, const Dbn2D& dbn_oflow,
const std::string& path="", const std::string& title="")
: AnalysisObject("Profile1D", path, title),
_axis(bins, dbn_tot, dbn_uflow, dbn_oflow)
{ }
/// Assignment operator
Profile1D& operator = (const Profile1D& p1) {
AnalysisObject::operator = (p1); //< AO treatment of paths etc.
_axis = p1._axis;
return *this;
}
/// Make a copy on the stack
Profile1D clone() const {
return Profile1D(*this);
}
/// Make a copy on the heap, via 'new'
Profile1D* newclone() const {
return new Profile1D(*this);
}
//@}
/// Fill dimension of this data object
size_t dim() const { return 1; }
/// @name Modifiers
//@{
/// Fill histo by value and weight
- void fill(double x, double y, double weight=1.0, double fraction=1.0);
+ virtual void fill(double x, double y, double weight=1.0, double fraction=1.0);
/// Fill histo x bin i with the given y value and weight
- void fillBin(size_t i, double y, double weight=1.0, double fraction=1.0);
+ virtual void fillBin(size_t i, double y, double weight=1.0, double fraction=1.0);
/// @brief Reset the histogram
///
/// Keep the binning but set all bin contents and related quantities to zero
void reset() {
_axis.reset();
}
/// Rescale as if all fill weights had been different by factor @a scalefactor.
void scaleW(double scalefactor) {
_axis.scaleW(scalefactor);
}
/// Rescale as if all y values had been different by factor @a scalefactor.
void scaleY(double scalefactor) {
_axis.totalDbn().scaleY(scalefactor);
_axis.overflow().scaleY(scalefactor);
_axis.underflow().scaleY(scalefactor);
for (size_t i = 0; i < bins().size(); ++i)
bin(i).scaleY(scalefactor);
}
/// Merge together the bin range with indices from @a from to @a to, inclusive
void mergeBins(size_t from, size_t to) {
_axis.mergeBins(from, to);
}
/// Merge every group of n bins, starting from the LHS
void rebinBy(unsigned int n, size_t begin=0, size_t end=UINT_MAX) {
_axis.rebinBy(n, begin, end);
}
/// Overloaded alias for rebinBy
void rebin(unsigned int n, size_t begin=0, size_t end=UINT_MAX) {
rebinBy(n, begin, end);
}
/// Rebin to the given list of bin edges
void rebinTo(const std::vector<double>& newedges) {
_axis.rebinTo(newedges);
}
/// Overloaded alias for rebinTo
void rebin(const std::vector<double>& newedges) {
rebinTo(newedges);
}
/// Bin addition operator
void addBin(double xlow, double xhigh) {
_axis.addBin(xlow, xhigh);
}
/// Bin addition operator
void addBins(const std::vector<double> binedges) {
_axis.addBins(binedges);
}
// /// Bin addition operator
// void addBins(const std::vector<std::pair<double,double> > edges) {
// _axis.addBins(edges);
// }
/// Add a new bin, perhaps already populated: CAREFUL!
void addBin(const ProfileBin1D& b) { _axis.addBin(b); }
/// @brief Bins addition operator
///
/// Add multiple bins without resetting
void addBins(const Bins& bins) {
_axis.addBins(bins);
}
//@}
/// @name Bin accessors
//@{
/// Number of bins on this axis (not counting under/overflow)
size_t numBins() const { return bins().size(); }
/// Low edge of this histo's axis
double xMin() const { return _axis.xMin(); }
/// High edge of this histo's axis
double xMax() const { return _axis.xMax(); }
/// All bin edges on this histo's axis
///
/// @note This only returns the finite edges, i.e. -inf and +inf are removed
/// @todo Make the +-inf stripping controllable by a default-valued bool arg
const std::vector<double> xEdges() const { return _axis.xEdges(); }
/// Access the bin vector
std::vector<YODA::ProfileBin1D>& bins() { return _axis.bins(); }
/// Access the bin vector
const std::vector<YODA::ProfileBin1D>& bins() const { return _axis.bins(); }
/// Access a bin by index (non-const version)
ProfileBin1D& bin(size_t index) { return _axis.bins()[index]; }
/// Access a bin by index (const version)
const ProfileBin1D& bin(size_t index) const { return _axis.bins()[index]; }
/// Access a bin index by x-coordinate.
int binIndexAt(double x) { return _axis.binIndexAt(x); }
/// Access a bin by x-coordinate (const version)
const ProfileBin1D& binAt(double x) const { return _axis.binAt(x); }
/// Access summary distribution, including gaps and overflows (non-const version)
Dbn2D& totalDbn() { return _axis.totalDbn(); }
/// Access summary distribution, including gaps and overflows (const version)
const Dbn2D& totalDbn() const { return _axis.totalDbn(); }
/// Set summary distribution, mainly for persistency: CAREFUL!
void setTotalDbn(const Dbn2D& dbn) { _axis.setTotalDbn(dbn); }
/// Access underflow (non-const version)
Dbn2D& underflow() { return _axis.underflow(); }
/// Access underflow (const version)
const Dbn2D& underflow() const { return _axis.underflow(); }
/// Set underflow distribution, mainly for persistency: CAREFUL!
void setUnderflow(const Dbn2D& dbn) { _axis.setUnderflow(dbn); }
/// Access overflow (non-const version)
Dbn2D& overflow() { return _axis.overflow(); }
/// Access overflow (const version)
const Dbn2D& overflow() const { return _axis.overflow(); }
/// Set overflow distribution, mainly for persistency: CAREFUL!
void setOverflow(const Dbn2D& dbn) { _axis.setOverflow(dbn); }
//@}
/// @name Whole histo data
//@{
/// @todo Add integrals? Or are they too ambiguous to make a core function?
/// Get the number of fills
unsigned long numEntries(bool includeoverflows=true) const;
/// Get the effective number of fills
double effNumEntries(bool includeoverflows=true) const;
/// Get sum of weights in histo.
double sumW(bool includeoverflows=true) const;
/// Get sum of squared weights in histo.
double sumW2(bool includeoverflows=true) const;
/// Get the mean x
double xMean(bool includeoverflows=true) const;
/// Get the variance in x
double xVariance(bool includeoverflows=true) const;
/// Get the standard deviation in x
double xStdDev(bool includeoverflows=true) const {
return std::sqrt(xVariance(includeoverflows));
}
/// Get the standard error on <x>
double xStdErr(bool includeoverflows=true) const;
/// Get the RMS in x
double xRMS(bool includeoverflows=true) const;
//@}
/// @name Adding and subtracting histograms
//@{
/// Add another profile to this one
Profile1D& operator += (const Profile1D& toAdd) {
if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
_axis += toAdd._axis;
return *this;
}
/// Subtract another profile from this one
Profile1D& operator -= (const Profile1D& toSubtract) {
if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
_axis -= toSubtract._axis;
return *this;
}
inline bool operator == (const Profile1D& other){
return _axis == other._axis;
}
inline bool operator != (const Profile1D& other){
return ! operator == (other);
}
//@}
protected:
/// Access a bin by x-coordinate (non-const version)
ProfileBin1D& binAt(double x) { return _axis.binAt(x); }
private:
/// @name Bin data
//@{
/// The bins contained in this profile histogram
Axis1D<ProfileBin1D, Dbn2D> _axis;
//@}
};
/// @name Combining profile histos: global operators
//@{
/// Add two profile histograms
inline Profile1D add(const Profile1D& first, const Profile1D& second) {
Profile1D tmp = first;
if (first.path() != second.path()) tmp.setPath("");
tmp += second;
return tmp;
}
/// Add two profile histograms
inline Profile1D operator + (const Profile1D& first, const Profile1D& second) {
return add(first, second);
}
/// Subtract two profile histograms
inline Profile1D subtract(const Profile1D& first, const Profile1D& second) {
Profile1D tmp = first;
if (first.path() != second.path()) tmp.setPath("");
tmp -= second;
return tmp;
}
/// Subtract two profile histograms
inline Profile1D operator - (const Profile1D& first, const Profile1D& second) {
return subtract(first, second);
}
/// Divide two profile histograms
Scatter2D divide(const Profile1D& numer, const Profile1D& denom);
/// Divide two profile histograms
inline Scatter2D operator / (const Profile1D& numer, const Profile1D& denom) {
return divide(numer, denom);
}
//@}
}
#endif
diff --git a/include/YODA/Profile2D.h b/include/YODA/Profile2D.h
--- a/include/YODA/Profile2D.h
+++ b/include/YODA/Profile2D.h
@@ -1,433 +1,433 @@
// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2016 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_Profile2D_h
#define YODA_Profile2D_h
#include "YODA/AnalysisObject.h"
#include "YODA/ProfileBin2D.h"
#include "YODA/Dbn3D.h"
#include "YODA/Axis2D.h"
#include "YODA/Scatter3D.h"
#include "YODA/Exceptions.h"
#include <vector>
namespace YODA {
// Forward declarations
class Histo2D;
class Scatter3D;
/// Convenience typedef
typedef Axis2D<ProfileBin2D, Dbn3D> Profile2DAxis;
/// A two-dimensional profile histogram.
class Profile2D : public AnalysisObject {
public:
/// Convenience typedefs
typedef Profile2DAxis Axis;
typedef Axis::Bins Bins;
typedef ProfileBin2D Bin;
typedef Axis::Outflows Outflows;
/// @name Constructors
//@{
/// Default constructor
Profile2D(const std::string& path="", const std::string& title="")
: AnalysisObject("Profile2D", path, title),
_axis()
{ }
/// Constructor giving range and number of bins
Profile2D(size_t nbinsX, double lowerX, double upperX,
size_t nbinsY, double lowerY, double upperY,
const std::string& path="", const std::string& title="")
: AnalysisObject("Profile2D", path, title),
_axis(nbinsX, std::make_pair(lowerX, upperX), nbinsY, std::make_pair(lowerY, upperY))
{ }
/// Constructor giving explicit bin edges in the direction of X and Y
Profile2D(const std::vector<double>& xedges, const std::vector<double>& yedges,
const std::string& path="", const std::string& title="")
: AnalysisObject("Profile2D", path, title),
_axis(xedges, yedges)
{ }
/// Constructor accepting an explicit collection of bins.
Profile2D(const std::vector<Bin>& bins,
const std::string& path="", const std::string& title="")
: AnalysisObject("Profile2D", path, title),
_axis(bins)
{ }
/// A copy constructor with optional new path
/// @todo Also allow title setting from the constructor?
Profile2D(const Profile2D& p, const std::string& path="");
/// A constructor from a Scatter3D's binning, with optional new path
/// @todo Also allow title setting from the constructor?
Profile2D(const Scatter3D& s, const std::string& path="");
/// Constructor from a Histo2D's binning, with optional new path
/// @todo Also allow title setting from the constructor?
Profile2D(const Histo2D& h, const std::string& path="");
/// @brief State-setting constructor
///
/// Mainly intended for internal persistency use.
Profile2D(const std::vector<ProfileBin2D>& bins,
const Dbn3D& totalDbn,
const Outflows& outflows,
const std::string& path="", const std::string& title="")
: AnalysisObject("Profile2D", path, title),
_axis(bins, totalDbn, outflows)
{ }
/// Assignment operator
Profile2D& operator = (const Profile2D& p2) {
AnalysisObject::operator = (p2); //< AO treatment of paths etc.
_axis = p2._axis;
return *this;
}
/// Make a copy on the stack
Profile2D clone() const {
return Profile2D(*this);
}
/// Make a copy on the heap, via 'new'
Profile2D* newclone() const {
return new Profile2D(*this);
}
//@}
/// Fill dimension of this data object
size_t dim() const { return 2; }
/// @name Modifiers
//@{
/// Fill histo by value and weight
- void fill(double x, double y, double z, double weight=1.0, double fraction=1.0);
+ virtual void fill(double x, double y, double z, double weight=1.0, double fraction=1.0);
/// Fill histo x-y bin i with the given z value and weight
- void fillBin(size_t i, double z, double weight=1.0, double fraction=1.0);
+ virtual void fillBin(size_t i, double z, double weight=1.0, double fraction=1.0);
/// @brief Reset the histogram
///
/// Keep the binning but reset the statistics
void reset() {
_axis.reset();
}
/// Rescale as if all fill weights had been different by a @a scalefactor
void scaleW(double scalefactor) {
/// @todo Is this ScaledBy annotation needed?
setAnnotation("ScaledBy", annotation<double>("ScaledBy", 1.0) * scalefactor);
_axis.scaleW(scalefactor);
}
/// Rescale as if all z values had been different by factor @a scalefactor.
void scaleZ(double scalefactor) {
_axis.totalDbn().scaleZ(scalefactor);
/// @todo Need to rescale overflows too, when they exist.
// _axis.overflow().scaleZ(scalefactor);
// _axis.underflow().scaleZ(scalefactor);
for (size_t i = 0; i < bins().size(); ++i)
bin(i).scaleZ(scalefactor);
}
/// @todo TODO
// /// Merge together the bin range with indices from @a from to @a to, inclusive
// void mergeBins(size_t from, size_t to) {
// _axis.mergeBins(from, to);
// }
/// @todo TODO
// /// Merge every group of n bins, starting from the LHS
// void rebin(size_t n) {
// throw "IMPLEMENT!";
// //_axis.rebin(n);
// }
// /// @brief Bin addition operator
// ///
// /// Add a bin to the axis, described by its x and y ranges.
void addBin(Axis::EdgePair1D xrange, Axis::EdgePair1D yrange) {
_axis.addBin(xrange, yrange);
}
// /// @brief Bin addition operator
// ///
// /// Add a bin to the axis, possibly pre-populated
void addBin(const Bin& bin) {
_axis.addBin(bin);
}
/// @brief Bins addition operator
///
/// Add multiple bins from edge cuts without resetting
void addBins(const Axis::Edges& xcuts, const Axis::Edges& ycuts) {
_axis.addBins(xcuts, ycuts);
}
/// @brief Bins addition operator
///
/// Add multiple bins without resetting
void addBins(const Bins& bins) {
_axis.addBins(bins);
}
/// @todo TODO
// /// @brief Bin addition operator
// ///
// /// Add a set of bins delimiting coordinates of which are contained
// /// in binLimits vector.
// void addBin(const std::vector<Segment>& binLimits) {
// _axis.addBin(binLimits);
// }
void eraseBin(size_t index) {
_axis.eraseBin(index);
}
//@}
/// @name Bin accessors
//@{
/// Low x edge of this histo's axis
double xMin() const { return _axis.xMin(); }
/// High x edge of this histo's axis
double xMax() const { return _axis.xMax(); }
/// Low y edge of this histo's axis
double yMin() const { return _axis.yMin(); }
/// High y edge of this histo's axis
double yMax() const { return _axis.yMax(); }
/// Access the bin vector (non-const)
std::vector<YODA::ProfileBin2D>& bins() { return _axis.bins(); }
/// Access the bin vector (const)
const std::vector<YODA::ProfileBin2D>& bins() const { return _axis.bins(); }
/// Access a bin by index (non-const)
ProfileBin2D& bin(size_t index) { return _axis.bins()[index]; }
/// Access a bin by index (const)
const ProfileBin2D& bin(size_t index) const { return _axis.bins()[index]; }
/// Access a bin index by coordinate
int binIndexAt(double x, double y) { return _axis.binIndexAt(x, y); }
/// Access a bin by coordinate (const)
const ProfileBin2D& binAt(double x, double y) const { return _axis.binAt(x, y); }
/// Number of bins of this axis (not counting under/over flow)
size_t numBins() const { return _axis.bins().size(); }
/// Number of bins along the x axis
size_t numBinsX() const { return _axis.numBinsX(); }
/// Number of bins along the y axis
size_t numBinsY() const { return _axis.numBinsY(); }
/// Access summary distribution, including gaps and overflows (non-const version)
Dbn3D& totalDbn() { return _axis.totalDbn(); }
/// Access summary distribution, including gaps and overflows (const version)
const Dbn3D& totalDbn() const { return _axis.totalDbn(); }
/// Set summary distribution, including gaps and overflows
void setTotalDbn(const Dbn3D& dbn) { _axis.setTotalDbn(dbn); }
// /// @brief Access an outflow (non-const)
// ///
// /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow.
// /// (0,0) is not a valid overflow index pair, since it is in range for both x and y.
// Dbn3D& outflow(int ix, int iy) {
// return _axis.outflow(ix, iy);
// }
// /// @brief Access an outflow (const)
// ///
// /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow.
// /// (0,0) is not a valid overflow index pair, since it is in range for both x and y.
// const Dbn3D& outflow(int ix, int iy) const {
// return _axis.outflow(ix, iy);
// }
//@}
/// @name Whole histo data
//@{
/// Get the number of fills
unsigned long numEntries(bool includeoverflows=true) const;
/// Get the effective number of fills
double effNumEntries(bool includeoverflows=true) const;
/// Get sum of weights in histo
double sumW(bool includeoverflows=true) const;
/// Get the sum of squared weights in histo
double sumW2(bool includeoverflows=true) const;
/// Get the mean x
double xMean(bool includeoverflows=true) const;
/// Get the mean y
double yMean(bool includeoverflows=true) const;
/// Get the variance in x
double xVariance(bool includeoverflows=true) const;
/// Get the variance in y
double yVariance(bool includeoverflows=true) const;
/// Get the standard deviation in x
double xStdDev(bool includeoverflows=true) const {
return std::sqrt(xVariance(includeoverflows));
}
/// Get the standard deviation in y
double yStdDev(bool includeoverflows=true) const {
return std::sqrt(yVariance(includeoverflows));
}
/// Get the standard error on <x>
double xStdErr(bool includeoverflows=true) const;
/// Get the standard error on <y>
double yStdErr(bool includeoverflows=true) const;
/// Get the RMS in x
double xRMS(bool includeoverflows=true) const;
/// Get the RMS in y
double yRMS(bool includeoverflows=true) const;
//@}
/// @name Adding and subtracting histograms
//@{
/// Add another profile to this one
Profile2D& operator += (const Profile2D& toAdd) {
if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
_axis += toAdd._axis;
return *this;
}
/// Subtract another profile from this one
Profile2D& operator -= (const Profile2D& toSubtract) {
if (hasAnnotation("ScaledBy")) rmAnnotation("ScaledBy");
_axis -= toSubtract._axis;
return *this;
}
inline bool operator == (const Profile2D& other){
return _axis == other._axis;
}
inline bool operator != (const Profile2D& other){
return ! operator == (other);
}
//@}-
protected:
/// Access a bin by coordinate (non-const)
ProfileBin2D& binAt(double x, double y) { return _axis.binAt(x, y); }
private:
/// @name Bin data
//@{
/// The bins contained in this profile histogram
Axis2D<ProfileBin2D, Dbn3D> _axis;
//@}
};
/// @name Combining profile histos: global operators
//@{
/// Add two profile histograms
inline Profile2D add(const Profile2D& first, const Profile2D& second) {
Profile2D tmp = first;
if (first.path() != second.path()) tmp.setPath("");
tmp += second;
return tmp;
}
/// Add two profile histograms
inline Profile2D operator + (const Profile2D& first, const Profile2D& second) {
return add(first,second);
}
/// Subtract two profile histograms
inline Profile2D subtract(const Profile2D& first, const Profile2D& second) {
Profile2D tmp = first;
if (first.path() != second.path()) tmp.setPath("");
tmp -= second;
return tmp;
}
/// Subtract two profile histograms
inline Profile2D operator - (const Profile2D& first, const Profile2D& second) {
return subtract(first,second);
}
/// Divide two profile histograms
Scatter3D divide(const Profile2D& numer, const Profile2D& denom);
/// Divide two profile histograms
inline Scatter3D operator / (const Profile2D& numer, const Profile2D& denom) {
return divide(numer, denom);
}
//@}
}
#endif