HistogramFactory.h

// -*- C++ -*-
//
// HistogramFactory.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef LWH_HistogramFactory_H
#define LWH_HistogramFactory_H
//
// This is the declaration of the HistogramFactory class.
//
 
#include "AIHistogramFactory.h"
#include "Histogram1D.h"
#include "Histogram2D.h"
#include "Tree.h"
#include <string>
#include <stdexcept>
 
namespace LWH {
 
using namespace AIDA;
 
class HistogramFactory: public IHistogramFactory {
 
public:
 
  HistogramFactory(Tree & t)
    : tree(&t) {}
 
  virtual ~HistogramFactory() {}
 
  bool destroy(IBaseHistogram * hist) {
    IManagedObject * mo = dynamic_cast<IManagedObject *>(hist);
    if ( !mo ) return false;
    return tree->rm(tree->findPath(*mo));
  }
 
  ICloud1D * createCloud1D(const std::string &, const std::string &,
                           int = -1, const std::string & = "") {
    return error<ICloud1D>("ICloud1D");
  }
 
  ICloud1D * createCloud1D(const std::string &) {
    return error<ICloud1D>("ICloud1D");
  }
 
  ICloud1D * createCopy(const std::string &, const ICloud1D &) {
    return error<ICloud1D>("ICloud1D");
  }
 
  ICloud2D * createCloud2D(const std::string &, const std::string &, int = -1,
                           const std::string & = "") {
    return error<ICloud2D>("ICloud2D");
  }
 
 
  ICloud2D * createCloud2D(const std::string &) {
    return error<ICloud2D>("ICloud2D");
  }
 
  ICloud2D * createCopy(const std::string &, const ICloud2D &) {
    return error<ICloud2D>("ICloud2D");
  }
 
  ICloud3D * createCloud3D(const std::string &, const std::string &, int = -1,
                           const std::string & = "") {
    return error<ICloud3D>("ICloud3D");
  }
 
  ICloud3D * createCloud3D(const std::string &) {
    return error<ICloud3D>("ICloud3D");
  }
 
  ICloud3D * createCopy(const std::string &, const ICloud3D &) {
    return error<ICloud3D>("ICloud3D");
  }
 
  IHistogram1D *
  createHistogram1D(const std::string & path, const std::string & title,
                    int nBins, double lowerEdge, double upperEdge,
                    const std::string & = "") {
    Histogram1D * hist = new Histogram1D(nBins, lowerEdge, upperEdge);
    hist->setTitle(title);
    if ( !tree->insert(path, hist) ) {
      delete hist;
      hist = 0;
      throw std::runtime_error("LWH could not create histogram '"
                               + title + "'." );
    }
    return hist;
  }
 
  IHistogram1D *
  createHistogram1D(const std::string & pathAndTitle,
                    int nBins, double lowerEdge, double upperEdge) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createHistogram1D(pathAndTitle, title, nBins, lowerEdge, upperEdge);
  }
 
 
  IHistogram1D *
  createHistogram1D(const std::string & path, const std::string & title,
                    const std::vector<double> & binEdges,
                    const std::string & = "") {
    Histogram1D * hist = new Histogram1D(binEdges);
    hist->setTitle(title);
    if ( !tree->insert(path, hist) ) {
      delete hist;
      hist = 0;
      throw std::runtime_error("LWH could not create histogram '"
                               + title + "'." );
    }
    return hist;
  }
 
  IHistogram1D *
  createCopy(const std::string & path, const IHistogram1D & hist) {
    Histogram1D * h = new Histogram1D(dynamic_cast<const Histogram1D &>(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;
  }
 
  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;
  }
 
  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);
  }
 
  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;
  }
 
  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;
  }
 
  IHistogram3D * createHistogram3D(const std::string &, const std::string &,
                                   int, double, double, int, double, double,
                                   int, double, double,
                                   const std::string & = "") {
    return error<IHistogram3D>("IHistogram3D");
  }
 
  IHistogram3D * createHistogram3D(const std::string &, int, double, double,
                                   int, double, double, int, double, double) {
    return error<IHistogram3D>("IHistogram3D");
  }
 
  IHistogram3D * createHistogram3D(const std::string &, const std::string &,
                                   const std::vector<double> &,
                                   const std::vector<double> &,
                                   const std::vector<double> &,
                                   const std::string & = "") {
    return error<IHistogram3D>("IHistogram3D");
  }
 
  IHistogram3D * createCopy(const std::string &, const IHistogram3D &) {
    return error<IHistogram3D>("IHistogram3D");
  }
 
  IProfile1D * createProfile1D(const std::string &, const std::string &,
                               int, double, double, const std::string & = "") {
    return error<IProfile1D>("IProfile1D");
  }
 
  IProfile1D * createProfile1D(const std::string &, const std::string &,
                               int, double, double, double, double,
                               const std::string & = "") {
    return error<IProfile1D>("IProfile1D");
  }
 
  IProfile1D * createProfile1D(const std::string &, const std::string &,
                               const std::vector<double> &,
                               const std::string & = "") {
    return error<IProfile1D>("IProfile1D");
  }
 
  IProfile1D * createProfile1D(const std::string &, const std::string &,
                               const std::vector<double> &, double, double,
                               const std::string & = "") {
    return error<IProfile1D>("IProfile1D");
  }
 
  IProfile1D * createProfile1D(const std::string &, int, double, double) {
    return error<IProfile1D>("IProfile1D");
  }
 
  IProfile1D * createProfile1D(const std::string &,
                               int, double, double, double, double) {
    return error<IProfile1D>("IProfile1D");
  }
 
  IProfile1D * createCopy(const std::string &, const IProfile1D &) {
    return error<IProfile1D>("IProfile1D");
  }
 
  IProfile2D * createProfile2D(const std::string &, const std::string &,
                               int, double, double, int, double, double,
                               const std::string & = "") {
    return error<IProfile2D>("IProfile2D");
  }
 
  IProfile2D * createProfile2D(const std::string &, const std::string &,
                               int, double, double, int,
                               double, double, double, double,
                               const std::string & = "") {
    return error<IProfile2D>("IProfile2D");
  }
 
  IProfile2D * createProfile2D(const std::string &, const std::string &,
                               const std::vector<double> &,
                               const std::vector<double> &,
                               const std::string & = "") {
    return error<IProfile2D>("IProfile2D");
  }
 
  IProfile2D * createProfile2D(const std::string &, const std::string &,
                               const std::vector<double> &,
                               const std::vector<double> &,
                               double, double, const std::string & = "") {
    return error<IProfile2D>("IProfile2D");
  }
 
  IProfile2D * createProfile2D(const std::string &, int, double, double,
                               int, double, double) {
    return error<IProfile2D>("IProfile2D");
  }
 
  IProfile2D * createProfile2D(const std::string &, int, double, double,
                               int, double, double, double, double) {
    return error<IProfile2D>("IProfile2D");
  }
 
  IProfile2D * createCopy(const std::string &, const IProfile2D &) {
    return error<IProfile2D>("IProfile2D");
  }
 
  Histogram1D * add(const std::string & path,
                     const Histogram1D & hist1, const Histogram1D & hist2) {
    if ( !checkBins(hist1, hist2) ) return 0;
    Histogram1D * h = new Histogram1D(hist1);
    h->setTitle(path.substr(path.rfind('/') + 1));
    h->add(hist2);
    if ( !tree->insert(path, h) ) return 0;
    return h;
  }
    
  IHistogram1D * add(const std::string & path,
                     const IHistogram1D & hist1, const IHistogram1D & hist2) {
    return add(path, dynamic_cast<const Histogram1D &>(hist1),
               dynamic_cast<const Histogram1D &>(hist2));
  }
    
  Histogram1D * subtract(const std::string & path,
                    const Histogram1D & h1, const Histogram1D & h2) {
    if ( !checkBins(h1, h2) ) return 0;
    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] += h2.sumw2[i];
    }
    if ( !tree->insert(path, h) ) return 0;
    return h;
  }
 
  IHistogram1D * subtract(const std::string & path, const IHistogram1D & hist1,
                          const IHistogram1D & hist2) {
    return subtract(path, dynamic_cast<const Histogram1D &>(hist1),
                    dynamic_cast<const Histogram1D &>(hist2));
  }
 
  Histogram1D * multiply(const std::string & path,
                    const Histogram1D & h1, const Histogram1D & h2) {
    if ( !checkBins(h1, h2) ) return 0;
    Histogram1D * h = new Histogram1D(h1);
    h->setTitle(path.substr(path.rfind('/') + 1));
    for ( int i = 0; i < h->ax->bins() + 2; ++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];
    }
    if ( !tree->insert(path, h) ) return 0;
    return h;
  }
 
  IHistogram1D * multiply(const std::string & path, const IHistogram1D & hist1,
                          const IHistogram1D & hist2) {
    return multiply(path, dynamic_cast<const Histogram1D &>(hist1),
                    dynamic_cast<const Histogram1D &>(hist2));
  }
 
  Histogram1D * divide(const std::string & path,
                  const Histogram1D & h1, const Histogram1D & h2) {
    if ( !checkBins(h1, h2) ) return 0;
    Histogram1D * h = new Histogram1D(h1);
    h->setTitle(path.substr(path.rfind('/') + 1));
    for ( int i = 0; i < h->ax->bins() + 2; ++i ) {
      if ( h2.sum[i] == 0 || h2.sumw[i] == 0.0 ) {
        h->sum[i] = 0;
        h->sumw[i] = h->sumw2[i] = 0.0;
        continue;
      }
      h->sumw[i] /= h2.sumw[i];
      h->sumw2[i] = h1.sumw2[i]/(h2.sumw[i]*h2.sumw[i]) +
        h1.sumw[i]*h1.sumw[i]*h2.sumw2[i]/
        (h2.sumw[i]*h2.sumw[i]*h2.sumw[i]*h2.sumw[i]);
    }
    if ( !tree->insert(path, h) ) return 0;
    return h;
  }
 
  IHistogram1D * divide(const std::string & path, const IHistogram1D & hist1,
                        const IHistogram1D & hist2) {
    return divide(path, dynamic_cast<const Histogram1D &>(hist1),
                    dynamic_cast<const Histogram1D &>(hist2));
  }
 
  inline bool _neq(double a, double b, double eps = 1e-5) const {
    using std::abs;
    if ( a == 0 && b == 0 ) return false;
    if ( abs(a-b) < eps*(abs(a) + abs(b)) ) return false;
    return true;
  }
 
  bool checkBins(const Histogram1D & h1, const Histogram1D & h2) const {
    if ( _neq(h1.ax->upperEdge(), h2.ax->upperEdge()) ||
         _neq(h1.ax->lowerEdge(), h2.ax->lowerEdge()) ||
         _neq(h1.ax->bins(), h2.ax->bins()) ) return false;
    if ( h1.fax && h2.fax ) return true;
    for ( int i = 0; i < h1.ax->bins(); ++i ) {
      if ( _neq(h1.ax->binUpperEdge(i), h2.ax->binUpperEdge(i)) ||
           _neq(h1.ax->binLowerEdge(i), h2.ax->binLowerEdge(i)) ) return false;
    }
    return true;
  }
 
  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;
    if ( h1.xfax && h2.xfax && h1.yfax && h2.yfax ) return true;
    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;
  }
 
  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));
  }
 
  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;    
  }
 
  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->sumw2[ix][iy] += h2.sumw2[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;
  }
 
  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));
  }
 
  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));
  }
 
  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->sum[ix][iy] *= h2.sum[ix][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;
  }
 
  Histogram2D * divide(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 ) {
      if ( h2.sum[ix][iy] == 0 || h2.sumw[ix][iy] == 0.0 ) {
        h->sum[ix][iy] = 0;
        h->sumw[ix][iy] = h->sumw2[ix][iy] = 0.0;
        continue;
      }
      h->sumw[ix][iy] /= h2.sumw[ix][iy];
      h->sumw2[ix][iy] = h1.sumw2[ix][iy]/(h2.sumw[ix][iy]*h2.sumw[ix][iy]) +
        h1.sumw[ix][iy]*h1.sumw[ix][iy]*h2.sumw2[ix][iy]/
        (h2.sumw[ix][iy]*h2.sumw[ix][iy]*h2.sumw[ix][iy]*h2.sumw[ix][iy]);
    }
    if ( !tree->insert(path, h) ) {
      delete h;
      return 0;
    }
    return h;
  }
 
 
  IHistogram2D * divide(const std::string & path,
                        const IHistogram2D & h1, const IHistogram2D & h2) {
    return divide(path, dynamic_cast<const Histogram2D &>(h1),
                  dynamic_cast<const Histogram2D &>(h2));
  }
 
  IHistogram3D * add(const std::string &,
                     const IHistogram3D &, const IHistogram3D &) {
    return error<IHistogram3D>("3D histograms");
  }
 
  IHistogram3D * subtract(const std::string &,
                          const IHistogram3D &, const IHistogram3D &) {
    return error<IHistogram3D>("3D histograms");
  }
 
  IHistogram3D * multiply(const std::string &,
                          const IHistogram3D &, const IHistogram3D &) {
    return error<IHistogram3D>("3D histograms");
  }
 
  IHistogram3D * divide(const std::string &,
                        const IHistogram3D &, const IHistogram3D &) {
    return error<IHistogram3D>("3D histograms");
  }
 
  IHistogram1D * projectionX(const std::string & path, const IHistogram2D & h) {
    return projectionX(path, dynamic_cast<const Histogram2D &>(h));
  }
 
  Histogram1D * projectionX(const std::string & path, const Histogram2D & h) {
    return sliceX(path, h, 0, h.yax->bins() - 1);
  }
 
  IHistogram1D * projectionY(const std::string & path, const IHistogram2D & h) {
    return projectionY(path, dynamic_cast<const Histogram2D &>(h));
  }
 
  Histogram1D * projectionY(const std::string & path, const Histogram2D & h) {
    return sliceY(path, h, 0, h.xax->bins() - 1);
  }
 
  IHistogram1D *
  sliceX(const std::string & path, const IHistogram2D & h, int i) {
    return sliceX(path, dynamic_cast<const Histogram2D &>(h), i, i);
  }
 
  Histogram1D *
  sliceX(const std::string & path, const Histogram2D & h, int i) {
    return sliceX(path, h, i, i);
  }
 
  IHistogram1D *
  sliceY(const std::string & path, const IHistogram2D & h, int i) {
    return sliceY(path, dynamic_cast<const Histogram2D &>(h), i, i);
  }
 
  Histogram1D * sliceY(const std::string & path, const Histogram2D & h, int i) {
    return sliceY(path, h, i, i);
  }
 
  IHistogram1D *
  sliceX(const std::string & path, const IHistogram2D & h, int il, int iu) {
    return sliceX(path, dynamic_cast<const Histogram2D &>(h), il, iu);
  }
 
  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;
  }
 
  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;
  }
 
  IHistogram2D * projectionXY(const std::string &, const IHistogram3D &) {
    return error<IHistogram2D>("2D histograms");
  }
 
  IHistogram2D * projectionXZ(const std::string &, const IHistogram3D &) {
    return error<IHistogram2D>("2D histograms");
  }
 
  IHistogram2D * projectionYZ(const std::string &, const IHistogram3D &) {
    return error<IHistogram2D>("2D histograms");
  }
 
  IHistogram2D * sliceXY(const std::string &, const IHistogram3D &, int, int) {
    return error<IHistogram2D>("2D histograms");
  }
 
  IHistogram2D * sliceXZ(const std::string &, const IHistogram3D &, int, int) {
    return error<IHistogram2D>("2D histograms");
  }
 
  IHistogram2D * sliceYZ(const std::string &, const IHistogram3D &, int, int) {
    return error<IHistogram2D>("2D histograms");
  }
 
 
private:
 
  template <typename T>
  static T * error(std::string feature) {
    throw std::runtime_error("LWH cannot handle " + feature + ".");
  }
 
  Tree * tree;
 
};
 
}
 
#endif /* LWH_HistogramFactory_H */