HistogramFactory.h

// -*- C++ -*-
//
// HistogramFactory.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2007 Leif Lonnblad
//
// ThePEG is licenced under version 2 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 "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 &, const std::string &,
                                   int, double, double, int, double, double,
                                   const std::string & = "") {
    return error<IHistogram2D>("IHistogram2D");
  }

  IHistogram2D * createHistogram2D(const std::string &,
                                   int, double, double, int, double, double) {
    return error<IHistogram2D>("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 * createCopy(const std::string &, const IHistogram2D &) {
    return error<IHistogram2D>("IHistogram2D");
  }

  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->sum[i] *= h2.sum[i];
      h->sumw[i] *= h2.sumw[i];
      h->sumxw[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->sum[i] /= h2.sum[i];
      h->sumw[i] /= h2.sumw[i];
      h->sumxw[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));
  }

  bool checkBins(const Histogram1D & h1, const Histogram1D & h2) const {
    if ( h1.ax->upperEdge() != h2.ax->upperEdge() ||
         h1.ax->lowerEdge() != h2.ax->lowerEdge() ||
         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 ( h1.ax->binUpperEdge(i) != h2.ax->binUpperEdge(i) ||
           h1.ax->binLowerEdge(i) != h2.ax->binLowerEdge(i) ) return false;
    }
    return true;
  }

  IHistogram2D * add(const std::string &,
                     const IHistogram2D &, const IHistogram2D &) {
    return error<IHistogram2D>("2D histograms");
  }

  IHistogram2D * subtract(const std::string &,
                          const IHistogram2D &, const IHistogram2D &) {
    return error<IHistogram2D>("2D histograms");
  }

  IHistogram2D * multiply(const std::string &,
                          const IHistogram2D &, const IHistogram2D &) {
    return error<IHistogram2D>("2D histograms");
  }

  IHistogram2D * divide(const std::string &,
                        const IHistogram2D &, const IHistogram2D &) {
    return error<IHistogram2D>("2D histograms");
  }

  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 &, const IHistogram2D &) {
    return error<IHistogram1D>("2D histograms");
  }

  IHistogram1D * projectionY(const std::string &, const IHistogram2D &) {
    return error<IHistogram1D>("2D histograms");
  }

  IHistogram1D * sliceX(const std::string &, const IHistogram2D &, int) {
    return error<IHistogram1D>("2D histograms");
  }

  IHistogram1D * sliceY(const std::string &, const IHistogram2D &, int) {
    return error<IHistogram1D>("2D histograms");
  }

  IHistogram1D * sliceX(const std::string &, const IHistogram2D &, int, int) {
    return error<IHistogram1D>("2D histograms");
  }

  IHistogram1D * sliceY(const std::string &, const IHistogram2D &, int, int) {
    return error<IHistogram1D>("2D histograms");
  }

  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 */

---