|
|
[Rivet-svn] r2705 - in trunk: include include/LWH include/Rivet src/Core
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Tue Oct 5 17:38:09 BST 2010
Author: lonnblad
Date: Tue Oct 5 17:38:08 2010
New Revision: 2705
Log:
Added 2D histograms to LWH.
Added:
trunk/include/LWH/AIHistogram2D.h
trunk/include/LWH/Histogram2D.h
Modified:
trunk/include/LWH/AIBaseHistogram.h
trunk/include/LWH/AIDataPointSet.h
trunk/include/LWH/AIHistogramFactory.h
trunk/include/LWH/DataPointSet.h
trunk/include/LWH/DataPointSetFactory.h
trunk/include/LWH/HistogramFactory.h
trunk/include/Makefile.am
trunk/include/Rivet/Analysis.hh
trunk/include/Rivet/RivetAIDA.fhh
trunk/include/Rivet/RivetAIDA.hh
trunk/src/Core/Analysis.cc
Modified: trunk/include/LWH/AIBaseHistogram.h
==============================================================================
--- trunk/include/LWH/AIBaseHistogram.h Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/LWH/AIBaseHistogram.h Tue Oct 5 17:38:08 2010 (r2705)
@@ -87,11 +87,32 @@
return true;
}
+ /**
+ * Get the z-axis title.
+ * @return The title.
+ *
+ */
+ std::string ztitle() const {
+ return theZTitle;
+ }
+
+ /**
+ * Set the z-axis title.
+ * @param title The new title.
+ * @return false If the title cannot be set.
+ *
+ */
+ bool setZTitle(const std::string & ztitle) {
+ theZTitle = ztitle;
+ return true;
+ }
+
protected:
std::string theTitle;
std::string theXTitle;
std::string theYTitle;
+ std::string theZTitle;
};
Modified: trunk/include/LWH/AIDataPointSet.h
==============================================================================
--- trunk/include/LWH/AIDataPointSet.h Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/LWH/AIDataPointSet.h Tue Oct 5 17:38:08 2010 (r2705)
@@ -113,11 +113,60 @@
return true;
}
+ /**
+ * Get the z-axis title.
+ * @return The title.
+ *
+ */
+ std::string ztitle() const {
+ return theZTitle;
+ }
+
+ /**
+ * Set the z-axis title.
+ * @param title The new title.
+ * @return false If the title cannot be set.
+ *
+ */
+ bool setZTitle(const std::string & ztitle) {
+ theZTitle = ztitle;
+ return true;
+ }
+
+ /**
+ * Get the z-axis title.
+ * @return The title.
+ *
+ */
+ std::string title(int i) const {
+ switch ( i ) {
+ case 0: return theXTitle;
+ case 1: return theYTitle;
+ case 2: return theZTitle;
+ default: return "";
+ }
+ }
+
+ /**
+ * Set the <i>-axis title.
+ * @param title The new title.
+ * @return false If the title cannot be set.
+ *
+ */
+ bool setTitle(int i, const std::string & title) {
+ if ( i == 0 ) theXTitle = title;
+ else if ( i == 1 ) theYTitle = title;
+ else if ( i == 2 ) theZTitle = title;
+ else return false;
+ return true;
+ }
+
protected:
std::string theTitle;
std::string theXTitle;
std::string theYTitle;
+ std::string theZTitle;
};
Added: trunk/include/LWH/AIHistogram2D.h
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/include/LWH/AIHistogram2D.h Tue Oct 5 17:38:08 2010 (r2705)
@@ -0,0 +1,191 @@
+// -*- C++ -*-
+#ifndef LWH_AIHistogram2D_H
+#define LWH_AIHistogram2D_H
+
+#ifndef LWH_USING_AIDA
+
+/** @cond DONT_DOCUMENT_STRIPPED_DOWN_AIDA_INTERFACES */
+
+#include "AIHistogram1D.h"
+
+namespace AIDA {
+
+class IAxis;
+
+class IHistogram2D : virtual public IHistogram {
+
+public:
+
+ virtual ~IHistogram2D() { /* nop */; }
+
+ /**
+ * Fill the IHistogram2D with a couple of values and the
+ * corresponding weight.
+ * @param x The x value to be filled in.
+ * @param y The y value to be filled in.
+ * @param weight The corresponding weight (by default 1).
+ * @return false If the weight is <0 or >1 (?).
+ *
+ */
+ virtual bool fill(double x, double y, double weight = 1.) = 0;
+
+ /**
+ * The weighted mean along the x axis of a given bin.
+ * @param indexX The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param indexY The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The mean of the corresponding bin along the x axis.
+ *
+ */
+ virtual double binMeanX(int indexX, int indexY) const = 0;
+
+ /**
+ * The weighted mean along the y axis of a given bin.
+ * @param indexX The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param indexY The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The mean of the corresponding bin along the y axis.
+ *
+ */
+ virtual double binMeanY(int indexX, int indexY) const = 0;
+
+ /**
+ * Number of entries in the corresponding bin (ie the number of times fill was called for this bin).
+ * @param indexX The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param indexY The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The number of entries in the corresponding bin.
+ *
+ */
+ virtual int binEntries(int indexX, int indexY) const = 0;
+
+ /**
+ * Sum of all the entries of the bins along a given x bin.
+ * This is equivalent to <tt>projectionX().binEntries(index)</tt>.
+ * @param index The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The number of entries in the corresponding set of bins.
+ *
+ */
+ virtual int binEntriesX(int index) const = 0;
+
+ /**
+ * Sum of all the entries of the bins along a given y bin.
+ * This is equivalent to <tt>projectionY().binEntries(index)</tt>.
+ * @param index The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The number of entries in the corresponding set of bins.
+ *
+ */
+ virtual int binEntriesY(int index) const = 0;
+
+ /**
+ * Total height of a give bin (ie the sum of the weights in this bin).
+ * @param indexX The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param indexY The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The height of the corresponding bin.
+ *
+ */
+ virtual double binHeight(int indexX, int indexY) const = 0;
+
+ /**
+ * Sum of all the heights of the bins along a given x bin.
+ * This is equivalent to <tt>projectionX().binHeight(index)</tt>.
+ * @param index The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The sum of the heights in the corresponding set of bins.
+ *
+ */
+ virtual double binHeightX(int index) const = 0;
+
+ /**
+ * Sum of all the heights of the bins along a given y bin.
+ * This is equivalent to <tt>projectionY().binHeight(index)</tt>.
+ * @param index The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The sum of the heights in the corresponding set of bins.
+ *
+ */
+ virtual double binHeightY(int index) const = 0;
+
+ /**
+ * The error of a given bin.
+ * @param indexX The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param indexY The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The error on the corresponding bin.
+ *
+ */
+ virtual double binError(int indexX, int indexY) const = 0;
+
+ /**
+ * The mean of the IHistogram2D along the x axis.
+ * @return The mean of the IHistogram2D along the x axis.
+ *
+ */
+ virtual double meanX() const = 0;
+
+ /**
+ * The mean of the IHistogram2D along the y axis.
+ * @return The mean of the IHistogram2D along the y axis.
+ *
+ */
+ virtual double meanY() const = 0;
+
+ /**
+ * The RMS of the IHistogram2D along the x axis.
+ * @return The RMS if the IHistogram2D along the x axis.
+ *
+ */
+ virtual double rmsX() const = 0;
+
+ /**
+ * The RMS of the IHistogram2D along the y axis.
+ * @return The RMS if the IHistogram2D along the y axis.
+ *
+ */
+ virtual double rmsY() const = 0;
+
+ /**
+ * Get the x axis of the IHistogram2D.
+ * @return The x coordinate IAxis.
+ *
+ */
+ virtual const IAxis & xAxis() const = 0;
+
+ /**
+ * Get the y axis of the IHistogram2D.
+ * @return The y coordinate IAxis.
+ *
+ */
+ virtual const IAxis & yAxis() const = 0;
+
+ /**
+ * Get the bin number corresponding to a given coordinate along the x axis.
+ * This is a convenience method, equivalent to <tt>xAxis().coordToIndex(coord)</tt>.
+ * @see IAxis#coordToIndex(double)
+ * @param coord The coordinalte along the x axis.
+ * @return The corresponding bin number.
+ *
+ */
+ virtual int coordToIndexX(double coord) const = 0;
+
+ /**
+ * Get the bin number corresponding to a given coordinate along the y axis.
+ * This is a convenience method, equivalent to <tt>yAxis().coordToIndex(coord)</tt>.
+ * @see IAxis#coordToIndex(double)
+ * @param coord The coordinalte along the y axis.
+ * @return The corresponding bin number.
+ *
+ */
+ virtual int coordToIndexY(double coord) const = 0;
+
+ /**
+ * Add to this IHistogram2D the contents of another IHistogram2D.
+ * @param hist The IHistogram2D to be added to this IHistogram2D.
+ * @return false If the IHistogram2Ds binnings are incompatible.
+ *
+ */
+ virtual bool add(const IHistogram2D & hist) = 0;
+}; // class
+} // namespace AIDA
+/** @endcond */
+
+#else
+#include "AIDA/IHistogram2D.h"
+#endif
+
+#endif /* LWH_AIHistogram2D_H */
+
Modified: trunk/include/LWH/AIHistogramFactory.h
==============================================================================
--- trunk/include/LWH/AIHistogramFactory.h Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/LWH/AIHistogramFactory.h Tue Oct 5 17:38:08 2010 (r2705)
@@ -67,7 +67,26 @@
// virtual IDataPointSet * divide(const std::string &, const IProfile1D &,
// const IProfile1D &) = 0;
-
+ virtual IHistogram2D *
+ createHistogram2D(const std::string & path, const std::string & title,
+ int nx, double xlo, double xup,
+ int ny, double ylo, double yup,
+ const std::string & = "") = 0;
+
+ virtual IHistogram2D *
+ createHistogram2D(const std::string & pathAndTitle,
+ int nx, double xlo, double xup,
+ int ny, double ylo, double yup) = 0;
+
+ virtual IHistogram2D *
+ createHistogram2D(const std::string & path, const std::string & title,
+ const std::vector<double> & xedges,
+ const std::vector<double> & yedges,
+ const std::string & = "") = 0;
+
+ virtual IHistogram2D *
+ createCopy(const std::string & path, const IHistogram2D & hist) = 0;
+
};
}
Modified: trunk/include/LWH/DataPointSet.h
==============================================================================
--- trunk/include/LWH/DataPointSet.h Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/LWH/DataPointSet.h Tue Oct 5 17:38:08 2010 (r2705)
@@ -332,8 +332,9 @@
<< "\"\n title=\"" << encodeForXML(title())
<< "\" path=\"" << path
<< "\" dimension=\"" << dimension() << "\">\n";
- os << " <dimension dim=\"0\" title=\"" << encodeForXML(xtitle()) << "\" />\n";
- os << " <dimension dim=\"1\" title=\"" << encodeForXML(ytitle()) << "\" />\n";
+ for ( int d = 0; d < dimension(); ++d )
+ os << " <dimension dim=\"" << d << "\" title=\""
+ << encodeForXML(title(d)) << "\" />\n";
for ( int i = 0, N = size(); i < N; ++i ) {
os << " <dataPoint>\n";
for ( int j = 0, M = dimension(); j < M; ++j )
Modified: trunk/include/LWH/DataPointSetFactory.h
==============================================================================
--- trunk/include/LWH/DataPointSetFactory.h Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/LWH/DataPointSetFactory.h Tue Oct 5 17:38:08 2010 (r2705)
@@ -8,6 +8,7 @@
#include "AIDataPointSetFactory.h"
#include "DataPointSet.h"
#include "Histogram1D.h"
+#include "Histogram2D.h"
#include "Profile1D.h"
#include "Tree.h"
#include <string>
@@ -58,7 +59,7 @@
if ( !tree->insert(path, dset) ) {
delete dset;
dset = 0;
- throw std::runtime_error("LWH could not create DataPointSet '"
+ throw std::runtime_error("LWH couldn't create DataPointSet '"
+ title + "'." );
}
return dset;
@@ -399,11 +400,11 @@
const std::vector<double> & exm, const std::vector<double> & eym,
const std::vector<double> & ezm) {
IDataPointSet * dset = create(path, title, 3);
- for ( int i = 0, N = y.size(); i < N; ++i ) dset->addPoint(DataPoint(3));
+ for ( int i = 0, N = z.size(); i < N; ++i ) dset->addPoint(DataPoint(3));
if ( !dset->setCoordinate(0, x, exp, exm) ||
!dset->setCoordinate(1, y, eyp, eym) ||
!dset->setCoordinate(2, z, ezp, ezm) )
- throw std::runtime_error("LWH could add points to DataPointSet '" +
+ throw std::runtime_error("LWH could not add points to DataPointSet '" +
title + "'." );
return dset;
}
@@ -562,6 +563,56 @@
/**
+ * Create an IDataPointSet from an IHistogram2D.
+ * @param path The path of the IDataPointSet. The path can either
+ * be a relative or full path.
+ * ("/folder1/folder2/dataName" and "../folder/dataName"
+ * are valid paths). All the directories in the path
+ * must exist. The characther `/` cannot be used in
+ * names; it is only used to delimit directories within
+ * paths.
+ * @param hist The IHistogram2D from which the data is taken.
+ * @param options Options, currently not specified
+ * @return The newly created IDataPointSet.
+ */
+ virtual IDataPointSet *
+ create(const std::string & path, const IHistogram2D & hist,
+ const std::string & = "") {
+ IDataPointSet * dset = create(path, hist.title(), 3);
+ //std::cout << hist.xtitle() << " & " << hist.ytitle() << std::endl;
+
+ dset->setXTitle(hist.xtitle());
+ dset->setYTitle(hist.ytitle());
+ std::vector<double> x, y, z, ex, ey, ez;
+ for ( int ix = 2, Nx = hist.xAxis().bins() + 2; ix < Nx; ++ix )
+ for ( int iy = 2, Ny = hist.yAxis().bins() + 2; iy < Ny; ++iy ) {
+ dset->addPoint(DataPoint(3));
+ //x.push_back(hist.binMean(i - 2)); // < "Dynamic" version
+ // Shouldn't IAxis have a binCentre(size_t binId) method? (According to Java AIDA v3.3.0 API)
+ x.push_back((hist.xAxis().binLowerEdge(ix - 2) +
+ hist.xAxis().binUpperEdge(ix - 2))/2.0);
+ ex.push_back(hist.xAxis().binWidth(ix - 2)/2.0);
+ /// @todo This is not really the height or error: width needs to be included...
+ y.push_back((hist.yAxis().binLowerEdge(iy - 2) +
+ hist.yAxis().binUpperEdge(iy - 2))/2.0);
+ ey.push_back(hist.yAxis().binWidth(iy - 2)/2.0);
+ /// @todo This is not really the height or error: width needs to be included...
+ const double binwidth = hist.xAxis().binWidth(ix - 2)*
+ hist.yAxis().binWidth(iy - 2);
+ z.push_back(hist.binHeight(ix - 2, iy - 2)/binwidth);
+ //ey.push_back(hist.binError(i - 2)/2.0/binwidth);
+ ez.push_back(hist.binError(ix - 2, iy - 2)/binwidth);
+ }
+ if ( !dset->setCoordinate(0, x, ex, ex) ||
+ !dset->setCoordinate(1, y, ey, ey) ||
+ !dset->setCoordinate(2, z, ez, ez) )
+ throw std::runtime_error("LWH could not add points to DataPointSet '" +
+ hist.title() + "'." );
+ return dset;
+ }
+
+
+ /**
* Create an IDataPointSet from an IProfile1D.
* @param path The path of the IDataPointSet. The path can either
* be a relative or full path.
@@ -601,14 +652,6 @@
/**
- * LWH cannot handle a IHistogram2D.
- */
- virtual IDataPointSet * create(const std::string &, const IHistogram2D &,
- const std::string & = "") {
- return error<IDataPointSet>("IHistogram2D");
- }
-
- /**
* LWH cannot handle a IHistogram3D.
*/
virtual IDataPointSet * create(const std::string &, const IHistogram3D &,
Added: trunk/include/LWH/Histogram2D.h
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/include/LWH/Histogram2D.h Tue Oct 5 17:38:08 2010 (r2705)
@@ -0,0 +1,873 @@
+// -*- C++ -*-
+#ifndef LWH_Histogram2D_H
+#define LWH_Histogram2D_H
+//
+// This is the declaration of the Histogram1D class.
+//
+
+#include "AIHistogram2D.h"
+#include "ManagedObject.h"
+#include "Axis.h"
+#include "VariAxis.h"
+#include <vector>
+#include <stdexcept>
+
+#include <iostream>
+#ifdef HAVE_ROOT
+ #include "TH2D.h"
+#endif
+
+
+
+namespace LWH {
+
+ using namespace AIDA;
+
+
+ /**
+ * User level interface to 1D Histogram.
+ */
+ class Histogram2D: public IHistogram2D, public ManagedObject {
+
+ public:
+
+ /** HistFactory is a friend. */
+ friend class HistogramFactory;
+
+ public:
+
+ /**
+ * Standard constructor.
+ */
+ Histogram2D(int nx, double lox, double upx,
+ int ny, double loy, double upy)
+ : xfax(new Axis(nx, lox, upx)), xvax(0), yfax(new Axis(ny, loy, upy)),
+ sum(nx + 2, std::vector<int>(ny + 2)),
+ sumw(nx + 2, std::vector<double>(ny + 2)),
+ sumw2(nx + 2, std::vector<double>(ny + 2)),
+ sumxw(nx + 2, std::vector<double>(ny + 2)),
+ sumx2w(nx + 2, std::vector<double>(ny + 2)),
+ sumyw(nx + 2, std::vector<double>(ny + 2)),
+ sumy2w(nx + 2, std::vector<double>(ny + 2)) {
+ xax = xfax;
+ yax = yfax;
+ }
+
+ /**
+ * Standard constructor for variable bin width.
+ */
+ Histogram2D(const std::vector<double> & xedges,
+ const std::vector<double> & yedges)
+ : xfax(0), xvax(new VariAxis(xedges)),
+ yfax(0), yvax(new VariAxis(xedges)),
+ sum(xedges.size() + 1, std::vector<int>(yedges.size() + 1)),
+ sumw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
+ sumw2(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
+ sumxw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
+ sumx2w(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
+ sumyw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
+ sumy2w(xedges.size() + 1, std::vector<double>(yedges.size() + 1)) {
+ xax = xvax;
+ yax = yvax;
+ }
+
+ /**
+ * Copy constructor.
+ */
+ Histogram2D(const Histogram2D & h)
+ : IBaseHistogram(h), IHistogram(h), IHistogram2D(h), ManagedObject(h),
+ xfax(0), xvax(0), yfax(0), yvax(0),
+ sum(h.sum), sumw(h.sumw), sumw2(h.sumw2),
+ sumxw(h.sumxw), sumx2w(h.sumx2w) ,
+ sumyw(h.sumyw), sumy2w(h.sumy2w){
+ const VariAxis * hxvax = dynamic_cast<const VariAxis *>(h.xax);
+ if ( hxvax ) xax = xvax = new VariAxis(*hxvax);
+ else xax = xfax = new Axis(dynamic_cast<const Axis &>(*h.xax));
+ const VariAxis * hyvax = dynamic_cast<const VariAxis *>(h.yax);
+ if ( hyvax ) yax = yvax = new VariAxis(*hyvax);
+ else yax = yfax = new Axis(dynamic_cast<const Axis &>(*h.yax));
+ }
+
+ /// Destructor.
+ virtual ~Histogram2D() {
+ delete xax;
+ delete yax;
+ }
+
+ /**
+ * Get the Histogram's title.
+ * @return The Histogram's title.
+ */
+ // std::string title() const {
+ // return theTitle;
+ // }
+
+ /**
+ * Get the Histogram's name.
+ * @return The Histogram's name
+ */
+ std::string name() const {
+ return title();
+ }
+
+ /**
+ * Set the histogram title.
+ * @param title The title.
+ * @return false If title cannot be changed.
+ */
+ // bool setTitle(const std::string & title) {
+ // theTitle = title;
+ // return true;
+ // }
+
+ /**
+ * Not implemented in LWH. will throw an exception.
+ */
+ IAnnotation & annotation() {
+ throw std::runtime_error("LWH cannot handle annotations");
+ return *anno;
+ }
+
+ /**
+ * Not implemented in LWH. will throw an exception.
+ */
+ const IAnnotation & annotation() const {
+ throw std::runtime_error("LWH cannot handle annotations");
+ return *anno;
+ }
+
+ /**
+ * Get the Histogram's dimension.
+ * @return The Histogram's dimension.
+ */
+ int dimension() const {
+ return 2;
+ }
+
+ /**
+ * Reset the Histogram; as if just created.
+ * @return false If something goes wrong.
+ */
+ bool reset() {
+ const int nx = xax->bins() + 2;
+ const int ny = yax->bins() + 2;
+ sum = std::vector< std::vector<int> >(nx, std::vector<int>(ny));
+ sumw = std::vector< std::vector<double> >(nx, std::vector<double>(ny));
+ sumw2 = sumw;
+ sumxw = sumw;
+ sumx2w = sumw;
+ sumyw = sumw;
+ sumy2w = sumw;
+ return true;
+ }
+
+ /**
+ * Get the number of in-range entries in the Histogram.
+ * @return The number of in-range entries.
+ *
+ */
+ int entries() const {
+ int si = 0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy ) si += sum[ix][iy];
+ return si;
+ }
+
+ /**
+ * Sum of the entries in all the IHistogram's bins,
+ * i.e in-range bins, UNDERFLOW and OVERFLOW.
+ * This is equivalent to the number of times the
+ * method fill was invoked.
+ * @return The sum of all the entries.
+ */
+ int allEntries() const {
+ return entries() + extraEntries();
+ }
+
+ /**
+ * Number of entries in the UNDERFLOW and OVERFLOW bins.
+ * @return The number of entries outside the range of the IHistogram.
+ */
+ int extraEntries() const {
+ int esum = sum[0][0] + sum[1][0] + sum[0][1] + sum[1][1];
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ esum += sum[ix][0] + sum[ix][1];
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy )
+ esum += sum[0][iy] + sum[1][iy];
+ return esum;
+ }
+
+ /**
+ * Number of equivalent entries,
+ * i.e. <tt>SUM[ weight ] ^ 2 / SUM[ weight^2 ]</tt>
+ * @return The number of equivalent entries.
+ */
+ double equivalentBinEntries() const {
+ double sw = 0.0;
+ double sw2 = 0.0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
+ sw += sumw[ix][iy];
+ sw2 += sumw2[ix][iy];
+ }
+ return sw2/(sw*sw);
+ }
+
+ /**
+ * Sum of in-range bin heights in the IHistogram,
+ * UNDERFLOW and OVERFLOW bins are excluded.
+ * @return The sum of the in-range bins heights.
+ *
+ */
+ double sumBinHeights() const {
+ double sw = 0.0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy ) sw += sumw[ix][iy];
+ return sw;
+ }
+
+ /**
+ * Sum of the heights of all the IHistogram's bins,
+ * i.e in-range bins, UNDERFLOW and OVERFLOW.
+ * @return The sum of all the bins heights.
+ */
+ double sumAllBinHeights() const {
+ return sumBinHeights() + sumExtraBinHeights();
+ }
+
+ /**
+ * Sum of heights in the UNDERFLOW and OVERFLOW bins.
+ * @return The sum of the heights of the out-of-range bins.
+ */
+ double sumExtraBinHeights() const {
+ int esum = sumw[0][0] + sumw[1][0] + sumw[0][1] + sumw[1][1];
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ esum += sumw[ix][0] + sumw[ix][1];
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy )
+ esum += sumw[0][iy] + sumw[1][iy];
+ return esum;
+ }
+
+ /**
+ * Minimum height of the in-range bins,
+ * i.e. not considering the UNDERFLOW and OVERFLOW bins.
+ * @return The minimum height among the in-range bins.
+ */
+ double minBinHeight() const {
+ double minw = sumw[2][2];
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy )
+ minw = std::min(minw, sumw[ix][iy]);
+ return minw;
+ }
+
+ /**
+ * Maximum height of the in-range bins,
+ * i.e. not considering the UNDERFLOW and OVERFLOW bins.
+ * @return The maximum height among the in-range bins.
+ */
+ double maxBinHeight() const{
+ double maxw = sumw[2][2];
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy )
+ maxw = std::max(maxw, sumw[ix][iy]);
+ return maxw;
+ }
+
+ /**
+ * Fill the IHistogram1D with a value and the
+ * corresponding weight.
+ * @param x The value to be filled in.
+ * @param weight The corresponding weight (by default 1).
+ * @return false If the weight is <0 or >1 (?).
+ */
+ bool fill(double x, double y, double weight = 1.) {
+ int ix = xax->coordToIndex(x) + 2;
+ int iy = yax->coordToIndex(y) + 2;
+ ++sum[ix][iy];
+ sumw[ix][iy] += weight;
+ sumxw[ix][iy] += x*weight;
+ sumx2w[ix][iy] += x*x*weight;
+ sumyw[ix][iy] += y*weight;
+ sumy2w[ix][iy] += y*y*weight;
+ sumw2[ix][iy] += weight*weight;
+ return weight >= 0 && weight <= 1;
+ }
+
+ /**
+ * The weighted mean along the x-axis of a bin.
+ * @param xindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param yindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The mean of the corresponding bin.
+ */
+ double binMeanX(int xindex, int yindex) const {
+ int ix = xindex + 2;
+ int iy = yindex + 2;
+ return sumw[ix][iy] != 0.0? sumxw[ix][iy]/sumw[ix][iy]:
+ ( xvax? xvax->binMidPoint(xindex): xfax->binMidPoint(xindex) );
+ };
+
+ /**
+ * The weighted mean along the y-axis of a bin.
+ * @param xindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param yindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The mean of the corresponding bin.
+ */
+ double binMeanY(int xindex, int yindex) const {
+ int ix = xindex + 2;
+ int iy = yindex + 2;
+ return sumw[ix][iy] != 0.0? sumyw[ix][iy]/sumw[ix][iy]:
+ ( yvax? yvax->binMidPoint(yindex): xfax->binMidPoint(yindex) );
+ };
+
+ /**
+ * The weighted x-RMS of a bin.
+ * @param xindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param yindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The RMS of the corresponding bin.
+ */
+ double binRmsX(int xindex, int yindex) const {
+ int ix = xindex + 2;
+ int iy = yindex + 2;
+ return sumw[ix][iy] == 0.0 || sum[ix][iy] < 2? xax->binWidth(xindex):
+ std::sqrt(std::max(sumw[ix][iy]*sumx2w[ix][iy] -
+ sumxw[ix][iy]*sumxw[ix][iy], 0.0))/sumw[ix][iy];
+ };
+
+ /**
+ * The weighted y-RMS of a bin.
+ * @param xindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @param yindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The RMS of the corresponding bin.
+ */
+ double binRmsY(int xindex, int yindex) const {
+ int ix = xindex + 2;
+ int iy = yindex + 2;
+ return sumw[ix][iy] == 0.0 || sum[ix][iy] < 2? yax->binWidth(yindex):
+ std::sqrt(std::max(sumw[ix][iy]*sumy2w[ix][iy] -
+ sumyw[ix][iy]*sumyw[ix][iy], 0.0))/sumw[ix][iy];
+ };
+
+ /**
+ * Number of entries in the corresponding bin (ie the number of
+ * times fill was called for this bin).
+ * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The number of entries in the corresponding bin.
+ */
+ int binEntries(int xindex, int yindex) const {
+ return sum[xindex + 2][yindex + 2];
+ }
+
+ /**
+ * Sum of all the entries of the bins along a given x bin.
+ * This is equivalent to <tt>projectionX().binEntries(index)</tt>.
+ * @param index The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The number of entries in the corresponding set of bins.
+ *
+ */
+ virtual int binEntriesX(int index) const {
+ int ret = 0;
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy )
+ ret += sum[index + 2][iy];
+ return ret;
+ }
+
+ /**
+ * Sum of all the entries of the bins along a given y bin.
+ * This is equivalent to <tt>projectionY().binEntries(index)</tt>.
+ * @param index The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The number of entries in the corresponding set of bins.
+ *
+ */
+ virtual int binEntriesY(int index) const {
+ int ret = 0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ ret += sum[ix][index + 2];
+ return ret;
+ }
+
+ /**
+ * Total height of the corresponding bin (ie the sum of the weights
+ * in this bin).
+ * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The height of the corresponding bin.
+ */
+ double binHeight(int xindex, int yindex) const {
+ /// @todo While this is compatible with the reference AIDA
+ /// implementation, it is not the bin height!
+ return sumw[xindex + 2][yindex + 2];
+ }
+
+ /**
+ * Sum of all the heights of the bins along a given x bin.
+ * This is equivalent to <tt>projectionX().binHeight(index)</tt>.
+ * @param index The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The sum of the heights in the corresponding set of bins.
+ *
+ */
+ virtual double binHeightX(int index) const {
+ double ret = 0;
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy )
+ ret += sumw[index + 2][iy];
+ return ret;
+ }
+
+ /**
+ * Sum of all the heights of the bins along a given y bin.
+ * This is equivalent to <tt>projectionY().binHeight(index)</tt>.
+ * @param index The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The sum of the heights in the corresponding set of bins.
+ *
+ */
+ virtual double binHeightY(int index) const {
+ double ret = 0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ ret += sumw[ix][index + 2];
+ return ret;
+ }
+
+ /**
+ * The error of a given bin.
+ * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
+ * @return The error on the corresponding bin.
+ *
+ */
+ double binError(int xindex, int yindex) const {
+ return std::sqrt(sumw2[xindex + 2][yindex + 2]);
+ }
+
+ /**
+ * The mean of the IHistogram2D along the x axis.
+ * @return The mean of the IHistogram2D along the x axis.
+ *
+ */
+ double meanX() const {
+ double s = 0.0;
+ double sx = 0.0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
+ s += sumw[ix][iy];
+ sx += sumxw[ix][iy];
+ }
+ return s != 0.0? sx/s: 0.0;
+ }
+
+ /**
+ * The mean of the IHistogram2D along the y axis.
+ * @return The mean of the IHistogram2D along the y axis.
+ *
+ */
+ double meanY() const {
+ double s = 0.0;
+ double sy = 0.0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
+ s += sumw[ix][iy];
+ sy += sumyw[ix][iy];
+ }
+ return s != 0.0? sy/s: 0.0;
+ }
+
+ /**
+ * The RMS of the IHistogram2D along the x axis.
+ * @return The RMS if the IHistogram2D along the x axis.
+ *
+ */
+ double rmsX() const {
+ double s = 0.0;
+ double sx = 0.0;
+ double sx2 = 0.0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
+ s += sumw[ix][iy];
+ sx += sumxw[ix][iy];
+ sx2 += sumx2w[ix][iy];
+ }
+ return s != 0.0? std::sqrt(std::max(s*sx2 - sx*sx, 0.0))/s:
+ xax->upperEdge() - xax->lowerEdge();
+ }
+
+ /**
+ * The RMS of the IHistogram2D along the x axis.
+ * @return The RMS if the IHistogram2D along the x axis.
+ *
+ */
+ double rmsY() const {
+ double s = 0.0;
+ double sy = 0.0;
+ double sy2 = 0.0;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
+ s += sumw[ix][iy];
+ sy += sumyw[ix][iy];
+ sy2 += sumy2w[ix][iy];
+ }
+ return s != 0.0? std::sqrt(std::max(s*sy2 - sy*sy, 0.0))/s:
+ yax->upperEdge() - yax->lowerEdge();
+ }
+
+ /** The weights. */
+ double getSumW(int xindex, int yindex) const {
+ return sumw[xindex + 2][yindex + 2];
+ }
+
+ /** The squared weights. */
+ double getSumW2(int xindex, int yindex) const {
+ return sumw2[xindex + 2][yindex + 2];
+ }
+
+ /** The weighted x-values. */
+ double getSumXW(int xindex, int yindex) const {
+ return sumxw[xindex + 2][yindex + 2];
+ }
+
+ /** The weighted x-square-values. */
+ double getSumX2W(int xindex, int yindex) const {
+ return sumx2w[xindex + 2][yindex + 2];
+ }
+
+ /** The weighted x-values. */
+ double getSumYW(int xindex, int yindex) const {
+ return sumyw[xindex + 2][yindex + 2];
+ }
+
+ /** The weighted x-square-values. */
+ double getSumY2W(int xindex, int yindex) const {
+ return sumy2w[xindex + 2][yindex + 2];
+ }
+
+ /**
+ * Get the x axis of the IHistogram2D.
+ * @return The x coordinate IAxis.
+ */
+ const IAxis & xAxis() const {
+ return *xax;
+ }
+
+ /**
+ * Get the y axis of the IHistogram2D.
+ * @return The y coordinate IAxis.
+ */
+ const IAxis & yAxis() const {
+ return *yax;
+ }
+
+ /**
+ * Get the bin number corresponding to a given coordinate along the
+ * x axis. This is a convenience method, equivalent to
+ * <tt>axis().coordToIndex(coord)</tt>.
+ * @param coord The coordinalte along the x axis.
+ * @return The corresponding bin number.
+ */
+ int coordToIndexX(double coord) const {
+ return xax->coordToIndex(coord);
+ }
+
+ /**
+ * Get the bin number corresponding to a given coordinate along the
+ * y axis. This is a convenience method, equivalent to
+ * <tt>axis().coordToIndex(coord)</tt>.
+ * @param coord The coordinalte along the y axis.
+ * @return The corresponding bin number.
+ */
+ int coordToIndexY(double coord) const {
+ return yax->coordToIndex(coord);
+ }
+
+ /**
+ * Add to this Histogram2D the contents of another IHistogram2D.
+ * @param h The Histogram2D to be added to this IHistogram2D.
+ * @return false If the IHistogram1Ds binnings are incompatible.
+ */
+ bool add(const Histogram2D & h) {
+ if ( xax->upperEdge() != h.xax->upperEdge() ||
+ xax->lowerEdge() != h.xax->lowerEdge() ||
+ xax->bins() != h.xax->bins() ) return false;
+ if ( yax->upperEdge() != h.yax->upperEdge() ||
+ yax->lowerEdge() != h.yax->lowerEdge() ||
+ yax->bins() != h.yax->bins() ) return false;
+ for ( int ix = 0; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
+ sum[ix][iy] += h.sum[ix][iy];
+ sumw[ix][iy] += h.sumw[ix][iy];
+ sumxw[ix][iy] += h.sumxw[ix][iy];
+ sumx2w[ix][iy] += h.sumx2w[ix][iy];
+ sumyw[ix][iy] += h.sumyw[ix][iy];
+ sumy2w[ix][iy] += h.sumy2w[ix][iy];
+ sumw2[ix][iy] += h.sumw2[ix][iy];
+ }
+ return true;
+ }
+
+ /**
+ * Add to this IHistogram1D the contents of another IHistogram1D.
+ * @param hist The IHistogram1D to be added to this IHistogram1D.
+ * @return false If the IHistogram1Ds binnings are incompatible.
+ */
+ bool add(const IHistogram2D & hist) {
+ return add(dynamic_cast<const Histogram2D &>(hist));
+ }
+
+ /**
+ * Scale the contents of this histogram with the given factor.
+ * @param s the scaling factor to use.
+ */
+ bool scale(double s) {
+ for ( int ix = 0; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
+ sumw[ix][iy] *= s;
+ sumxw[ix][iy] *= s;
+ sumx2w[ix][iy] *= s;
+ sumyw[ix][iy] *= s;
+ sumy2w[ix][iy] *= s;
+ sumw2[ix][iy] *= s*s;
+ }
+ return true;
+ }
+
+ /**
+ * Scale the given histogram so that the integral over all bins
+ * (including overflow) gives \a intg. This function also corrects
+ * for the bin-widths, which means that it should only be run once
+ * for each histogram. Further rescaling must be done with the
+ * scale(double) function.
+ */
+ void normalize(double intg) {
+ double oldintg = sumAllBinHeights();
+ if ( oldintg == 0.0 ) return;
+ for ( int ix = 0; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
+ double fac = intg/oldintg;
+ if ( ix >= 2 && iy >= 2 )
+ fac /= (xax->binUpperEdge(ix - 2) - xax->binLowerEdge(ix - 2))*
+ (yax->binUpperEdge(iy - 2) - yax->binLowerEdge(iy - 2));
+ sumw[ix][iy] *= fac;
+ sumxw[ix][iy] *= fac;
+ sumx2w[ix][iy] *= fac;
+ sumyw[ix][iy] *= fac;
+ sumy2w[ix][iy] *= fac;
+ sumw2[ix][iy] *= fac*fac;
+ }
+ }
+
+ /**
+ * Return the integral over the histogram bins assuming it has been
+ * normalize()d.
+ */
+ // double integral() const {
+ // double intg = sumw[0] + sumw[1];
+ // for ( int i = 2; i < ax->bins() + 2; ++i )
+
+ // is this right? Leave out bin width factor?
+
+ // intg += sumw[ix][iy]*(ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
+ // return intg;
+ // }
+
+ /**
+ * Not implemented in LWH.
+ * @return null pointer always.
+ */
+ void * cast(const std::string &) const {
+ return 0;
+ }
+
+ /**
+ * Write out the histogram in the AIDA xml format.
+ */
+ bool writeXML(std::ostream & os, std::string path, std::string name) {
+ //std::cout << "Writing out histogram " << name << " in AIDA file format!" << std::endl;
+ os << " <histogram2d name=\"" << encodeForXML(name)
+ << "\"\n title=\"" << encodeForXML(title())
+ << "\" path=\"" << path
+ << "\">\n <axis max=\"" << xax->upperEdge()
+ << "\" numberOfBins=\"" << xax->bins()
+ << "\" min=\"" << xax->lowerEdge()
+ << "\" direction=\"x\"";
+ if ( xvax ) {
+ os << ">\n";
+ for ( int i = 0, N = xax->bins() - 1; i < N; ++i )
+ os << " <binBorder value=\"" << xax->binUpperEdge(i) << "\"/>\n";
+ os << " </axis>\n";
+ } else {
+ os << "/>\n";
+ }
+ os << " <axis max=\"" << yax->upperEdge()
+ << "\" numberOfBins=\"" << yax->bins()
+ << "\" min=\"" << yax->lowerEdge()
+ << "\" direction=\"y\"";
+ if ( yvax ) {
+ os << ">\n";
+ for ( int i = 0, N = yax->bins() - 1; i < N; ++i )
+ os << " <binBorder value=\"" << yax->binUpperEdge(i) << "\"/>\n";
+ os << " </axis>\n";
+ } else {
+ os << "/>\n";
+ }
+ os << " <statistics entries=\"" << entries()
+ << "\">\n <statistic mean=\"" << meanX()
+ << "\" direction=\"x\"\n rms=\"" << rmsX()
+ << "\"/>\n </statistics>\n <statistic mean=\"" << meanY()
+ << "\" direction=\"y\"\n rms=\"" << rmsY()
+ << "\"/>\n </statistics>\n <data2d>\n";
+ for ( int ix = 0; ix < xax->bins() + 2; ++ix )
+ for ( int iy = 0; iy < yax->bins() + 2; ++iy )
+ if ( sum[ix][iy] ) {
+ os << " <bin2d binNumX=\"";
+ if ( ix == 0 ) os << "UNDERFLOW";
+ else if ( ix == 1 ) os << "OVERFLOW";
+ else os << ix - 2;
+ os << "\" binNumY=\"";
+ if ( iy == 0 ) os << "UNDERFLOW";
+ else if ( iy == 1 ) os << "OVERFLOW";
+ else os << iy - 2;
+ os << "\" entries=\"" << sum[ix][iy]
+ << "\" height=\"" << sumw[ix][iy]
+ << "\"\n error=\"" << std::sqrt(sumw2[ix][iy])
+ << "\" error2=\"" << sumw2[ix][iy]
+ << "\"\n weightedMeanX=\"" << binMeanX(ix - 2, iy - 2)
+ << "\" weightedRmsX=\"" << binRmsX(ix - 2, iy - 2)
+ << "\"\n weightedMeanY=\"" << binMeanY(ix - 2, iy - 2)
+ << "\" weightedRmsY=\"" << binRmsY(ix - 2, iy - 2)
+ << "\"/>\n";
+ }
+ os << " </data2d>\n </histogram2d>" << std::endl;
+ return true;
+ }
+
+
+ /**
+ * Write out the histogram in a flat text file suitable for
+ * eg. gnuplot to read. The coloums are layed out as 'x w w2 n'.
+ */
+ bool writeFLAT(std::ostream & os, std::string path, std::string name) {
+ os << "#2D " << path << "/" << name << " " << xax->lowerEdge()
+ << " " << xax->bins() << " " << xax->upperEdge() << " "
+ << yax->lowerEdge() << " " << yax->bins() << " " << yax->upperEdge()
+ << " \"" << title() << "\"" << std::endl;
+ for ( int ix = 2; ix < xax->bins() + 2; ++ix ) {
+ for ( int iy = 2; iy < yax->bins() + 2; ++iy )
+ os << 0.5*(xax->binLowerEdge(ix - 2)+xax->binUpperEdge(ix - 2)) << " "
+ << 0.5*(yax->binLowerEdge(iy - 2)+yax->binUpperEdge(iy - 2))
+ << " " << sumw[ix][iy] << " " << sqrt(sumw2[ix][iy])
+ << " " << sum[ix][iy] << std::endl;
+ os << std::endl;
+ }
+ os << std::endl;
+ return true;
+ }
+
+
+
+ #ifdef HAVE_ROOT
+ /**
+ * Write out the histogram in Root file format.
+ */
+ //bool writeROOT(std::ostream & os, std::string path, std::string name) {
+ bool writeROOT(TFile* file, std::string path, std::string name) {
+
+ //std::cout << "Writing out histogram " << name.c_str() << " in ROOT file format" << std::endl;
+
+ TH1D* hist1d;
+ int nbins;
+ if (!vax || vax->isFixedBinning() ) {//equidistant binning (easier case)
+ nbins = ax->bins();
+ hist1d = new TH1D(name.c_str(), title().c_str(), nbins, ax->lowerEdge(), ax->upperEdge());
+ }
+ else {
+ nbins = vax->bins();
+ double* bins = new double[nbins+1];
+ for (int i=0; i<nbins; ++i) {
+ bins[ix][iy] = vax->binEdges(i).first;
+ }
+ bins[nbins] = vax->binEdges(nbins-1).second; //take last bin right border
+ hist1d = new TH1D(name.c_str(), title().c_str(), nbins, bins);
+ delete bins;
+ }
+
+
+ double entries = 0;
+ for ( int i = 0; i < nbins + 2; ++i ) {
+ if ( sum[ix][iy] ) {
+ //i==0: underflow->RootBin(0), i==1: overflow->RootBin(NBins+1)
+ entries = entries + sum[ix][iy];
+ int j=i;
+ if (i==0) j=0; //underflow
+ else if (i==1) j=nbins+1; //overflow
+ if (i>=2) j=i-1; //normal bin entries
+ hist1d->SetBinContent(j, sumw[ix][iy]);
+ hist1d->SetBinError(j, sqrt(sumw2[ix][iy]));
+ //hist1d->Fill(binMean(i), sumw[ix][iy]);
+ }
+ }
+
+ hist1d->Sumw2();
+ hist1d->SetEntries(entries);
+
+ std::string DirName; //remove preceding slash from directory name, else ROOT error
+ for (unsigned int i=1; i<path.size(); ++i) DirName += path[i];
+ if (!file->Get(DirName.c_str())) file->mkdir(DirName.c_str());
+ file->cd(DirName.c_str());
+ hist1d->Write();
+
+ delete hist1d;
+
+ return true;
+ }
+
+ #endif
+
+
+
+ private:
+
+ /** The title */
+ // std::string theTitle;
+
+ /** The axis. */
+ IAxis * xax;
+
+ /** Pointer (possibly null) to a axis with fixed bin width. */
+ Axis * xfax;
+
+ /** Pointer (possibly null) to a axis with fixed bin width. */
+ VariAxis * xvax;
+
+ /** The axis. */
+ IAxis * yax;
+
+ /** Pointer (possibly null) to a axis with fixed bin width. */
+ Axis * yfax;
+
+ /** Pointer (possibly null) to a axis with fixed bin width. */
+ VariAxis * yvax;
+
+ /** The counts. */
+ std::vector< std::vector<int> > sum;
+
+ /** The weights. */
+ std::vector< std::vector<double> > sumw;
+
+ /** The squared weights. */
+ std::vector< std::vector<double> > sumw2;
+
+ /** The weighted x-values. */
+ std::vector< std::vector<double> > sumxw;
+
+ /** The weighted x-square-values. */
+ std::vector< std::vector<double> > sumx2w;
+
+ /** The weighted y-values. */
+ std::vector< std::vector<double> > sumyw;
+
+ /** The weighted y-square-values. */
+ std::vector< std::vector<double> > sumy2w;
+
+ /** dummy pointer to non-existen annotation. */
+ IAnnotation * anno;
+
+ };
+
+}
+
+#endif /* LWH_Histogram1D_H */
Modified: trunk/include/LWH/HistogramFactory.h
==============================================================================
--- trunk/include/LWH/HistogramFactory.h Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/LWH/HistogramFactory.h Tue Oct 5 17:38:08 2010 (r2705)
@@ -7,6 +7,7 @@
#include "AIHistogramFactory.h"
#include "Histogram1D.h"
+#include "Histogram2D.h"
#include "DataPointSet.h"
#include "Profile1D.h"
#include "Tree.h"
@@ -220,35 +221,65 @@
/**
* LWH cannot create a IHistogram2D.
*/
- IHistogram2D * createHistogram2D(const std::string &, const std::string &,
- int, double, double, int, double, double,
- const std::string & = "") {
- return error<IHistogram2D>("IHistogram2D");
+ IHistogram2D *
+ createHistogram2D(const std::string & path, const std::string & title,
+ int nx, double xlo, double xup,
+ int ny, double ylo, double yup,
+ const std::string & = "") {
+ Histogram2D * hist = new Histogram2D(nx, xlo, xup, ny, ylo, yup);
+ hist->setTitle(title);
+ if ( !tree->insert(path, hist) ) {
+ delete hist;
+ hist = 0;
+ throw std::runtime_error("LWH could not create histogram '"
+ + title + "'." );
+ }
+ return hist;
}
/**
* LWH cannot create a IHistogram2D.
*/
- IHistogram2D * createHistogram2D(const std::string &,
- int, double, double, int, double, double) {
- return error<IHistogram2D>("IHistogram2D");
+ IHistogram2D * createHistogram2D(const std::string & pathAndTitle,
+ int nx, double xlo, double xup,
+ int ny, double ylo, double yup) {
+ std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
+ return createHistogram2D(pathAndTitle, title, nx, xlo, xup, ny, ylo, yup);
}
/**
* LWH cannot create a IHistogram2D.
*/
- IHistogram2D * createHistogram2D(const std::string &, const std::string &,
- const std::vector<double> &,
- const std::vector<double> &,
- const std::string & = "") {
- return error<IHistogram2D>("IHistogram2D");
+ IHistogram2D *
+ createHistogram2D(const std::string & path, const std::string & title,
+ const std::vector<double> & xedges,
+ const std::vector<double> & yedges,
+ const std::string & = "") {
+ Histogram2D * hist = new Histogram2D(xedges, yedges);
+ hist->setTitle(title);
+ if ( !tree->insert(path, hist) ) {
+ delete hist;
+ hist = 0;
+ throw std::runtime_error("LWH could not create histogram '"
+ + title + "'." );
+ }
+ return hist;
}
/**
* LWH cannot create a copy of an IHistogram2D.
*/
- IHistogram2D * createCopy(const std::string &, const IHistogram2D &) {
- return error<IHistogram2D>("IHistogram2D");
+ IHistogram2D *
+ createCopy(const std::string & path, const IHistogram2D & hist) {
+ Histogram2D * h = new Histogram2D(dynamic_cast<const Histogram2D &>(hist));
+ h->setTitle(path.substr(path.rfind('/') + 1));
+ if ( !tree->insert(path, h) ) {
+ delete h;
+ h = 0;
+ throw std::runtime_error("LWH could not create a copy of histogram '"
+ + hist.title() + "'." );
+ }
+ return h;
}
/**
@@ -584,7 +615,6 @@
Histogram1D * h = new Histogram1D(h1);
h->setTitle(path.substr(path.rfind('/') + 1));
for ( int i = 0; i < h->ax->bins() + 2; ++i ) {
- h->sum[i] *= h2.sum[i];
h->sumw[i] *= h2.sumw[i];
h->sumw2[i] += h1.sumw[i]*h1.sumw[i]*h2.sumw2[i] +
h2.sumw[i]*h2.sumw[i]*h1.sumw2[i];
@@ -771,35 +801,178 @@
}
/**
+ * Check if two histograms have the same bins.
+ */
+ bool checkBins(const Histogram2D & h1, const Histogram2D & h2) const {
+ if (_neq( h1.xax->upperEdge(), h2.xax->upperEdge()) ||
+ _neq( h1.xax->lowerEdge(), h2.xax->lowerEdge()) ||
+ h1.xax->bins() != h2.xax->bins() ) return false;
+ if (_neq( h1.yax->upperEdge(), h2.yax->upperEdge()) ||
+ _neq( h1.yax->lowerEdge(), h2.yax->lowerEdge()) ||
+ h1.yax->bins() != h2.yax->bins() ) return false;
+ for ( int i = 0; i < h1.xax->bins(); ++i ) {
+ if ( _neq(h1.xax->binUpperEdge(i), h2.xax->binUpperEdge(i)) ||
+ _neq(h1.xax->binLowerEdge(i), h2.xax->binLowerEdge(i)) )
+ return false;
+ }
+ for ( int i = 0; i < h1.yax->bins(); ++i ) {
+ if ( _neq(h1.yax->binUpperEdge(i), h2.yax->binUpperEdge(i)) ||
+ _neq(h1.yax->binLowerEdge(i), h2.yax->binLowerEdge(i)) )
+ return false;
+ }
+ return true;
+ }
+
+ /**
* LWH cannot create an IHistogram2D by adding two IHistogram2D.
*/
- IHistogram2D * add(const std::string &,
- const IHistogram2D &, const IHistogram2D &) {
- return error<IHistogram2D>("2D histograms");
+ IHistogram2D * add(const std::string & path,
+ const IHistogram2D & hist1, const IHistogram2D & hist2) {
+ return add(path, dynamic_cast<const Histogram2D &>(hist1),
+ dynamic_cast<const Histogram2D &>(hist2));
+ }
+
+ /**
+ * LWH cannot create an IHistogram2D by adding two IHistogram2D.
+ */
+ Histogram2D * add(const std::string & path,
+ const Histogram2D & h1, const Histogram2D & h2) {
+ if ( !checkBins(h1, h2) ) return 0;
+ Histogram2D * h = new Histogram2D(h1);
+ h->setTitle(path.substr(path.rfind('/') + 1));
+ h->add(h2);
+ if ( !tree->insert(path, h) ) {
+ delete h;
+ return 0;
+ }
+ return h;
}
/**
* LWH cannot create an IHistogram2D by subtracting two IHistogram2D.
*/
- IHistogram2D * subtract(const std::string &,
- const IHistogram2D &, const IHistogram2D &) {
- return error<IHistogram2D>("2D histograms");
+ Histogram2D * subtract(const std::string & path,
+ const Histogram2D & h1, const Histogram2D & h2) {
+ if ( !checkBins(h1, h2) ) {
+ //std::cout << "!!!!!!!" << std::endl;
+ return 0;
+ }
+ Histogram2D * h = new Histogram2D(h1);
+ h->setTitle(path.substr(path.rfind('/') + 1));
+ for ( int ix = 0; ix < h->xax->bins() + 2; ++ix )
+ for ( int iy = 0; iy < h->yax->bins() + 2; ++iy ) {
+ h->sum[ix][iy] += h2.sum[ix][iy];
+ h->sumw[ix][iy] -= h2.sumw[ix][iy];
+ h->sumw[ix][iy] -= h2.sumw[ix][iy];
+ h->sumxw[ix][iy] -= h2.sumxw[ix][iy];
+ h->sumx2w[ix][iy] -= h2.sumx2w[ix][iy];
+ h->sumyw[ix][iy] -= h2.sumyw[ix][iy];
+ h->sumy2w[ix][iy] -= h2.sumy2w[ix][iy];
+ }
+ if ( !tree->insert(path, h) ) {
+ //std::cout << "&&&&&&&" << std::endl;
+ delete h;
+ return 0;
+ }
+ return h;
+ }
+
+ /**
+ * LWH cannot create an IHistogram2D by subtracting two IHistogram2D.
+ */
+ IHistogram2D * subtract(const std::string & path,
+ const IHistogram2D & h1, const IHistogram2D & h2) {
+ return subtract(path, dynamic_cast<const Histogram2D &>(h1),
+ dynamic_cast<const Histogram2D &>(h2));
}
/**
* LWH cannot create an IHistogram2D by multiplying two IHistogram2D.
*/
- IHistogram2D * multiply(const std::string &,
- const IHistogram2D &, const IHistogram2D &) {
- return error<IHistogram2D>("2D histograms");
+ IHistogram2D * multiply(const std::string & path,
+ const IHistogram2D & h1, const IHistogram2D & h2) {
+ return multiply(path, dynamic_cast<const Histogram2D &>(h1),
+ dynamic_cast<const Histogram2D &>(h2));
+ }
+
+ /**
+ * LWH cannot create an IHistogram2D by multiplying two IHistogram2D.
+ */
+ Histogram2D * multiply(const std::string & path,
+ const Histogram2D & h1, const Histogram2D & h2) {
+ if ( !checkBins(h1, h2) ) return 0;
+ Histogram2D * h = new Histogram2D(h1);
+ h->setTitle(path.substr(path.rfind('/') + 1));
+ for ( int ix = 0; ix < h->xax->bins() + 2; ++ix )
+ for ( int iy = 0; iy < h->yax->bins() + 2; ++iy ) {
+ h->sumw[ix][iy] *= h2.sumw[ix][iy];
+ h->sumw2[ix][iy] += h1.sumw[ix][iy]*h1.sumw[ix][iy]*h2.sumw2[ix][iy] +
+ h2.sumw[ix][iy]*h2.sumw[ix][iy]*h1.sumw2[ix][iy];
+ }
+ if ( !tree->insert(path, h) ) {
+ delete h;
+ return 0;
+ }
+ return h;
}
/**
* LWH cannot create an IHistogram2D by dividing two IHistogram2D.
*/
- IHistogram2D * divide(const std::string &,
- const IHistogram2D &, const IHistogram2D &) {
- return error<IHistogram2D>("2D histograms");
+ DataPointSet * divide(const std::string & path,
+ const Histogram2D & h1, const Histogram2D & h2) {
+ if ( !checkBins(h1,h2) ) return 0;
+ DataPointSet * h = new DataPointSet(3);
+ h->setTitle(path.substr(path.rfind('/') + 1));
+ for (int ix = 0; ix < h1.xax->bins(); ++ix) {
+ const double xbinwidth = h1.xax->binWidth(ix);
+ const double xbincentre =
+ ( h1.xax->binLowerEdge(ix) + h1.xax->binUpperEdge(ix) ) / 2.0;
+ for (int iy = 0; iy < h1.yax->bins(); ++iy) {
+ const double ybinwidth = h1.yax->binWidth(iy);
+ const double ybincentre =
+ ( h1.yax->binLowerEdge(iy) + h1.yax->binUpperEdge(iy) ) / 2.0;
+ IDataPoint* point = h->addPoint();
+ IMeasurement* x = point->coordinate(0);
+ x->setValue(xbincentre);
+ x->setErrorPlus(xbinwidth/2.0);
+ x->setErrorMinus(xbinwidth/2.0);
+ IMeasurement* y = point->coordinate(1);
+ y->setValue(ybincentre);
+ y->setErrorPlus(ybinwidth/2.0);
+ y->setErrorMinus(ybinwidth/2.0);
+
+ double zval(0), zerr(0);
+ if ( h1.binHeight(ix, iy) == 0 || h2.binHeight(ix, iy) == 0 ) {
+ /// @todo Bad way of handling div by zero!
+ zval = 0.0;
+ zerr = 0.0;
+ } else {
+ zval = h1.binHeight(ix, iy) / h2.binHeight(ix, iy);
+ double z1relerr = h1.binError(ix, iy)/h1.binHeight(ix,iy);
+ double z2relerr = h2.binError(ix, iy)/h2.binHeight(ix,iy);
+ zerr = zval * sqrt(pow(z1relerr, 2) + pow(z2relerr, 2));
+ }
+ IMeasurement* z = point->coordinate(2);
+ z->setValue(zval);
+ z->setErrorPlus(zerr);
+ z->setErrorMinus(zerr);
+ }
+ }
+ if ( !tree->insert(path, h) ) {
+ delete h;
+ return 0;
+ }
+ return h;
+ }
+
+ /**
+ * LWH cannot create an IHistogram2D by dividing two IHistogram2D.
+ */
+ IDataPointSet * divide(const std::string & path,
+ const IHistogram2D & h1, const IHistogram2D & h2) {
+ return divide(path, dynamic_cast<const Histogram1D &>(h1),
+ dynamic_cast<const Histogram1D &>(h2));
}
/**
@@ -838,48 +1011,145 @@
* LWH cannot create an IHistogram1D by projecting an IHistogram2D
* along its x axis.
*/
- IHistogram1D * projectionX(const std::string &, const IHistogram2D &) {
- return error<IHistogram1D>("2D histograms");
+ IHistogram1D * projectionX(const std::string & path, const IHistogram2D & h) {
+ return projectionX(path, dynamic_cast<const Histogram2D &>(h));
+ }
+
+ /**
+ * LWH cannot create an IHistogram1D by projecting an IHistogram2D
+ * along its x axis.
+ */
+ Histogram1D * projectionX(const std::string & path, const Histogram2D & h) {
+ return sliceX(path, h, h.xax->bins() - 1);
+ }
+
+ /**
+ * LWH cannot create an IHistogram1D by projecting an IHistogram2D
+ * along its y axis.
+ */
+ IHistogram1D * projectionY(const std::string & path, const IHistogram2D & h) {
+ return projectionY(path, dynamic_cast<const Histogram2D &>(h));
}
/**
* LWH cannot create an IHistogram1D by projecting an IHistogram2D
* along its y axis.
*/
- IHistogram1D * projectionY(const std::string &, const IHistogram2D &) {
- return error<IHistogram1D>("2D histograms");
+ Histogram1D * projectionY(const std::string & path, const Histogram2D & h) {
+ return sliceY(path, h, h.yax->bins() - 1);
}
/**
* LWH cannot create an IHistogram1D by slicing an IHistogram2D
* parallel to the y axis at a given bin.
*/
- IHistogram1D * sliceX(const std::string &, const IHistogram2D &, int) {
- return error<IHistogram1D>("2D histograms");
+ IHistogram1D *
+ sliceX(const std::string & path, const IHistogram2D & h, int i) {
+ return sliceX(path, dynamic_cast<const Histogram2D &>(h), i, i);
+ }
+
+ /**
+ * LWH cannot create an IHistogram1D by slicing an IHistogram2D
+ * parallel to the y axis at a given bin.
+ */
+ Histogram1D *
+ sliceX(const std::string & path, const Histogram2D & h, int i) {
+ return sliceX(path, h, i, i);
+ }
+
+ /**
+ * LWH cannot create an IHistogram1D by slicing an IHistogram2D
+ * parallel to the x axis at a given bin.
+ */
+ IHistogram1D *
+ sliceY(const std::string & path, const IHistogram2D & h, int i) {
+ return sliceY(path, dynamic_cast<const Histogram2D &>(h), i, i);
}
/**
* LWH cannot create an IHistogram1D by slicing an IHistogram2D
* parallel to the x axis at a given bin.
*/
- IHistogram1D * sliceY(const std::string &, const IHistogram2D &, int) {
- return error<IHistogram1D>("2D histograms");
+ Histogram1D * sliceY(const std::string & path, const Histogram2D & h, int i) {
+ return sliceY(path, h, i, i);
}
/**
* LWH cannot create an IHistogram1D by slicing an IHistogram2D
* parallel to the y axis between two bins (inclusive).
*/
- IHistogram1D * sliceX(const std::string &, const IHistogram2D &, int, int) {
- return error<IHistogram1D>("2D histograms");
+ IHistogram1D *
+ sliceX(const std::string & path, const IHistogram2D & h, int il, int iu) {
+ return sliceX(path, dynamic_cast<const Histogram2D &>(h), il, iu);
+ }
+
+ /**
+ * LWH cannot create an IHistogram1D by slicing an IHistogram2D
+ * parallel to the y axis between two bins (inclusive).
+ */
+ Histogram1D *
+ sliceX(const std::string & path, const Histogram2D & h2, int il, int iu) {
+ Histogram1D * h1;
+ if ( h2.xfax )
+ h1 = new Histogram1D(h2.xfax->bins(), h2.xfax->lowerEdge(),
+ h2.xfax->upperEdge());
+ else {
+ std::vector<double> edges(h2.xax->bins() + 1);
+ edges.push_back(h2.xax->lowerEdge());
+ for ( int i = 0; i < h2.xax->bins(); ++i )
+ edges.push_back(h2.xax->binLowerEdge(i));
+ h1 = new Histogram1D(edges);
+ }
+ for ( int ix = 0; ix < h2.xax->bins() + 2; ++ix )
+ for ( int iy = il + 2; iy <= iu + 2; ++iy ) {
+ h1->sum[ix] += h2.sum[ix][iy];
+ h1->sumw[ix] += h2.sumw[ix][iy];
+ h1->sumw2[ix] += h2.sumw2[ix][iy];
+ h1->sumxw[ix] += h2.sumxw[ix][iy];
+ h1->sumx2w[ix] += h2.sumx2w[ix][iy];
+ }
+ if ( !tree->insert(path, h1) ) {
+ delete h1;
+ return 0;
+ }
+ return h1;
}
/**
* LWH cannot create an IHistogram1D by slicing an IHistogram2D
* parallel to the x axis between two bins (inclusive).
*/
- IHistogram1D * sliceY(const std::string &, const IHistogram2D &, int, int) {
- return error<IHistogram1D>("2D histograms");
+ IHistogram1D *
+ sliceY(const std::string & path, const IHistogram2D & h, int il, int iu) {
+ return sliceY(path, dynamic_cast<const Histogram2D &>(h), il, iu);
+ }
+
+ Histogram1D *
+ sliceY(const std::string & path, const Histogram2D & h2, int il, int iu) {
+ Histogram1D * h1;
+ if ( h2.yfax )
+ h1 = new Histogram1D(h2.yfax->bins(), h2.yfax->lowerEdge(),
+ h2.yfax->upperEdge());
+ else {
+ std::vector<double> edges(h2.yax->bins() + 1);
+ edges.push_back(h2.yax->lowerEdge());
+ for ( int i = 0; i < h2.yax->bins(); ++i )
+ edges.push_back(h2.yax->binLowerEdge(i));
+ h1 = new Histogram1D(edges);
+ }
+ for ( int iy = 0; iy < h2.yax->bins() + 2; ++iy )
+ for ( int ix = il + 2; ix <= iu + 2; ++ix ) {
+ h1->sum[iy] += h2.sum[ix][iy];
+ h1->sumw[iy] += h2.sumw[ix][iy];
+ h1->sumw2[iy] += h2.sumw2[ix][iy];
+ h1->sumxw[iy] += h2.sumyw[ix][iy];
+ h1->sumx2w[iy] += h2.sumy2w[ix][iy];
+ }
+ if ( !tree->insert(path, h1) ) {
+ delete h1;
+ return 0;
+ }
+ return h1;
}
/**
Modified: trunk/include/Makefile.am
==============================================================================
--- trunk/include/Makefile.am Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/Makefile.am Tue Oct 5 17:38:08 2010 (r2705)
@@ -10,6 +10,7 @@
LWH/AIManagedObject.h \
LWH/AIBaseHistogram.h \
LWH/AIHistogram1D.h \
+ LWH/AIHistogram2D.h \
LWH/AIProfile1D.h \
LWH/AIDataPointSet.h \
LWH/AIMeasurement.h \
@@ -22,6 +23,7 @@
LWH/DataPointSetFactory.h \
LWH/ManagedObject.h \
LWH/Histogram1D.h \
+ LWH/Histogram2D.h \
LWH/Profile1D.h \
LWH/DataPointSet.h \
LWH/DataPoint.h \
Modified: trunk/include/Rivet/Analysis.hh
==============================================================================
--- trunk/include/Rivet/Analysis.hh Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/Rivet/Analysis.hh Tue Oct 5 17:38:08 2010 (r2705)
@@ -218,6 +218,22 @@
/// @warning The old histogram will be deleted, and its pointer set to zero.
void scale(AIDA::IHistogram1D*& histo, double scale);
+ /// Normalize the given histogram, @a histo. After this call the
+ /// histogram will have been transformed to a DataPointSet with the
+ /// same name and path. It has the same effect as
+ /// @c scale(histo, norm/sumOfWeights).
+ /// @param histo The histogram to be normalised.
+ /// @param norm The new area of the histogram.
+ /// @warning The old histogram will be deleted, and its pointer set to zero.
+ void normalize(AIDA::IHistogram2D*& histo, double norm=1.0);
+
+ /// Multiplicatively scale the given histogram, @a histo. After this call the
+ /// histogram will have been transformed to a DataPointSet with the same name and path.
+ /// @param histo The histogram to be scaled.
+ /// @param scale The factor used to multiply the histogram bin heights.
+ /// @warning The old histogram will be deleted, and its pointer set to zero.
+ void scale(AIDA::IHistogram2D*& histo, double scale);
+
/// Set the cross section from the generator
Analysis& setCrossSection(double xs);
@@ -301,6 +317,33 @@
const std::vector<double>& binedges, const std::string& title="",
const std::string& xtitle="", const std::string& ytitle="");
+ /// Book a 2D histogram with @a nxbins and @a nybins uniformly
+ /// distributed across the ranges @a xlower - @a xupper and @a
+ /// ylower - @a yupper respectively along the x- and y-axis.
+ /// (NB. this returns a pointer rather than a reference since it
+ /// will have to be stored in the analysis class - there's no
+ /// point in forcing users to explicitly get the pointer from a
+ /// reference before they can use it!)
+ AIDA::IHistogram2D*
+ bookHistogram2D(const std::string& name,
+ size_t nxbins, double xlower, double xupper,
+ size_t nybins, double ylower, double yupper,
+ const std::string& title="", const std::string& xtitle="",
+ const std::string& ytitle="", const std::string& ztitle="");
+
+ /// Book a 2D histogram with non-uniform bins defined by the
+ /// vectorx of bin edges @a xbinedges and @a ybinedges.
+ /// (NB. this returns a pointer rather than a reference since it
+ /// will have to be stored in the analysis class - there's no
+ /// point in forcing users to explicitly get the pointer from a
+ /// reference before they can use it!)
+ AIDA::IHistogram2D*
+ bookHistogram2D(const std::string& name,
+ const std::vector<double>& xbinedges,
+ const std::vector<double>& ybinedges,
+ const std::string& title="", const std::string& xtitle="",
+ const std::string& ytitle="", const std::string& ztitle="");
+
/// Book a 1D histogram based on the name in the corresponding AIDA
/// file. The binnings will be obtained by reading the bundled AIDA data
/// record file with the same filename as the analysis' name() property.
Modified: trunk/include/Rivet/RivetAIDA.fhh
==============================================================================
--- trunk/include/Rivet/RivetAIDA.fhh Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/Rivet/RivetAIDA.fhh Tue Oct 5 17:38:08 2010 (r2705)
@@ -9,6 +9,7 @@
class IAnalysisFactory;
class IHistogramFactory;
class IHistogram1D;
+ class IHistogram2D;
class IProfile1D;
class IDataPointSetFactory;
class IDataPointSet;
Modified: trunk/include/Rivet/RivetAIDA.hh
==============================================================================
--- trunk/include/Rivet/RivetAIDA.hh Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/include/Rivet/RivetAIDA.hh Tue Oct 5 17:38:08 2010 (r2705)
@@ -10,6 +10,7 @@
#include "LWH/AIAnalysisFactory.h"
#include "LWH/AIHistogramFactory.h"
#include "LWH/AIHistogram1D.h"
+#include "LWH/AIHistogram2D.h"
#include "LWH/AIProfile1D.h"
#include "LWH/AITreeFactory.h"
#include "LWH/AIDataPointSetFactory.h"
Modified: trunk/src/Core/Analysis.cc
==============================================================================
--- trunk/src/Core/Analysis.cc Tue Sep 28 18:12:56 2010 (r2704)
+++ trunk/src/Core/Analysis.cc Tue Oct 5 17:38:08 2010 (r2705)
@@ -378,6 +378,44 @@
return hist;
}
+ IHistogram2D*
+ Analysis::bookHistogram2D(const string& hname,
+ size_t nxbins, double xlower, double xupper,
+ size_t nybins, double ylower, double yupper,
+ const string& title, const string& xtitle,
+ const string& ytitle, const string& ztitle) {
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IHistogram2D* hist =
+ histogramFactory().createHistogram2D(path, title, nxbins, xlower, xupper,
+ nybins, ylower, yupper);
+ getLog() << Log::TRACE << "Made 2D histogram "
+ << hname << " for " << name() << endl;
+ hist->setXTitle(xtitle);
+ hist->setYTitle(ytitle);
+ hist->setZTitle(ztitle);
+ return hist;
+ }
+
+
+ IHistogram2D*
+ Analysis::bookHistogram2D(const string& hname,
+ const vector<double>& xbinedges,
+ const vector<double>& ybinedges,
+ const string& title, const string& xtitle,
+ const string& ytitle, const string& ztitle) {
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IHistogram2D* hist =
+ histogramFactory().createHistogram2D(path, title, xbinedges, ybinedges);
+ getLog() << Log::TRACE << "Made 2D histogram " << hname << " for "
+ << name() << endl;
+ hist->setXTitle(xtitle);
+ hist->setYTitle(ytitle);
+ hist->setZTitle(ztitle);
+ return hist;
+ }
+
/////////////////
@@ -576,4 +614,85 @@
}
+ void Analysis::normalize(AIDA::IHistogram2D*& histo, double norm) {
+ if (!histo) {
+ getLog() << Log::ERROR << "Failed to normalise histo=NULL in analysis "
+ << name() << " (norm=" << norm << ")" << endl;
+ return;
+ }
+ const string hpath = tree().findPath(dynamic_cast<const AIDA::IManagedObject&>(*histo));
+ getLog() << Log::TRACE << "Normalizing histo " << hpath << " to " << norm << endl;
+
+ double oldintg = 0.0;
+ int nxBins = histo->xAxis().bins();
+ int nyBins = histo->yAxis().bins();
+ for (int ixBin = 0; ixBin != nxBins; ++ixBin)
+ for (int iyBin = 0; iyBin != nyBins; ++iyBin) {
+ // Leaving out factor of binWidth because AIDA's "height"
+ // already includes a width factor.
+ oldintg += histo->binHeight(ixBin, iyBin); // * histo->axis().binWidth(iBin);
+ }
+ if (oldintg == 0.0) {
+ getLog() << Log::WARN << "Histo " << hpath
+ << " has null integral during normalisation" << endl;
+ return;
+ }
+
+ // Scale by the normalisation factor.
+ scale(histo, norm/oldintg);
+ }
+
+
+ void Analysis::scale(AIDA::IHistogram2D*& histo, double scale) {
+ if (!histo) {
+ getLog() << Log::ERROR << "Failed to scale histo=NULL in analysis "
+ << name() << " (scale=" << scale << ")" << endl;
+ return;
+ }
+ const string hpath =
+ tree().findPath(dynamic_cast<const AIDA::IManagedObject&>(*histo));
+ getLog() << Log::TRACE << "Scaling histo " << hpath << endl;
+
+ vector<double> x, y, z, ex, ey, ez;
+ for (size_t ix = 0, Nx = histo->xAxis().bins(); ix < Nx; ++ix)
+ for (size_t iy = 0, Ny = histo->yAxis().bins(); iy < Ny; ++iy) {
+ x.push_back(0.5 * (histo->xAxis().binLowerEdge(ix) +
+ histo->xAxis().binUpperEdge(ix)));
+ ex.push_back(histo->xAxis().binWidth(ix)*0.5);
+ y.push_back(0.5 * (histo->yAxis().binLowerEdge(iy) +
+ histo->yAxis().binUpperEdge(iy)));
+ ey.push_back(histo->yAxis().binWidth(iy)*0.5);
+
+ // "Bin height" is a misnomer in the AIDA spec: width is neglected.
+ // We'd like to do this: y.push_back(histo->binHeight(i) * scale);
+ z.push_back(histo->binHeight(ix, iy)*scale/
+ (histo->xAxis().binWidth(ix)*histo->yAxis().binWidth(iy)));
+ // "Bin error" is a misnomer in the AIDA spec: width is neglected.
+ // We'd like to do this: ey.push_back(histo->binError(i) * scale);
+ ez.push_back(histo->binError(ix, iy)*scale/
+ (histo->xAxis().binWidth(ix)*histo->yAxis().binWidth(iy)));
+ }
+
+ string title = histo->title();
+ string xtitle = histo->xtitle();
+ string ytitle = histo->ytitle();
+ string ztitle = histo->ztitle();
+
+ tree().mkdir("/tmpnormalize");
+ tree().mv(hpath, "/tmpnormalize");
+
+ AIDA::IDataPointSet* dps =
+ datapointsetFactory().createXYZ(hpath, title, x, y, z, ex, ey, ez);
+ dps->setXTitle(xtitle);
+ dps->setYTitle(ytitle);
+ dps->setZTitle(ztitle);
+
+ tree().rm(tree().findPath(dynamic_cast<AIDA::IManagedObject&>(*histo)));
+ tree().rmdir("/tmpnormalize");
+
+ // Set histo pointer to null - it can no longer be used.
+ histo = 0;
+ }
+
+
}
More information about the Rivet-svn
mailing list
|