DataPointSetFactory.h

// -*- C++ -*-
//
// DataPointSetFactory.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_DataPointSetFactory_H
#define LWH_DataPointSetFactory_H
//
// This is the declaration of the DataPointSetFactory class.
//

#include "AIDataPointSetFactory.h"
#include "DataPointSet.h"
#include "Histogram1D.h"
#include "Tree.h"
#include <string>
#include <stdexcept>

namespace LWH {

using namespace AIDA;

class DataPointSetFactory: public IDataPointSetFactory {

public: 

  DataPointSetFactory(Tree & t)
    : tree(&t) {}

  virtual ~DataPointSetFactory() {}

  virtual IDataPointSet *
  create(const std::string & path, const std::string & title, int dim) {
    DataPointSet * dset = new DataPointSet(dim);
    dset->setTitle(title);
    if ( !tree->insert(path, dset) ) {
      delete dset;
      dset = 0;
      throw std::runtime_error("LWH could not create DataPointSet '"
                               + title + "'." );
    }
    return dset;
  }

  virtual IDataPointSet * create(const std::string & pathAndTitle, int dim) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return create(pathAndTitle, title, dim);
  }

  virtual IDataPointSet *
  createY(const std::string & path, const std::string & title,
          const std::vector<double> & y, const std::vector<double> & ey) {
    return createY(path, title, y, ey, ey);
  }

  virtual IDataPointSet *
  createY(const std::string & path, const std::string & title,
          const std::vector<double> & y, const std::vector<double> & eyp,
          const std::vector<double>  & eym) {
    IDataPointSet * dset = create(path, title, 2);
    std::vector<double> x, ex;
    for ( int i = 0, N = y.size(); i < N; ++i ) {
      dset->addPoint(DataPoint(2));
      x.push_back(i);
      ex.push_back(0);
    }
    if ( !dset->setCoordinate(0, x, ex, ex) ||
         !dset->setCoordinate(1, y, eyp, eym) )
      throw std::runtime_error("LWH could add points to DataPointSet '" +
                               title +  "'." );
    return dset;
  }
    
  virtual IDataPointSet *
  createY(const std::string & pathAndTitle, const std::vector<double> & y,
          const std::vector<double> & ey) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createY(pathAndTitle, title, y, ey);
  }

  virtual IDataPointSet *
  createY(const std::string & pathAndTitle,
          const std::vector<double>  & y, const std::vector<double>  & eyp,
          const std::vector<double>  & eym) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createY(pathAndTitle, title, y, eyp, eym);
  }

  virtual IDataPointSet *
  createX(const std::string & path, const std::string & title,
          const std::vector<double>  & x, const std::vector<double>  & ex) {
    return createX(path, title, x, ex, ex);
  }

  virtual IDataPointSet *
  createX(const std::string & path, const std::string & title,
          const std::vector<double> & x, const std::vector<double> & exp,
          const std::vector<double>  & exm) {
    IDataPointSet * dset = create(path, title, 2);
    std::vector<double> y, ey;
    for ( int i = 0, N = x.size(); i < N; ++i ) {
      dset->addPoint(DataPoint(2));
      y.push_back(i);
      ey.push_back(0);
    }
    if ( !dset->setCoordinate(0, x, exp, exm) ||
         !dset->setCoordinate(1, y, ey, ey) )
      throw std::runtime_error("LWH could add points to DataPointSet '" +
                               title +  "'." );
    return dset;
  }

  virtual IDataPointSet *
  createX(const std::string & pathAndTitle, const std::vector<double> & x,
          const std::vector<double> & ex) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createX(pathAndTitle, title, x, ex, ex);
  }

  virtual IDataPointSet *
  createX(const std::string & pathAndTitle, const std::vector<double> & x,
          const std::vector<double> & exp, const std::vector<double> & exm) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createX(pathAndTitle, title, x, exp, exm);
  }

  virtual IDataPointSet *
  createXY(const std::string & path, const std::string & title,
           const std::vector<double> & x, const std::vector<double> & y,
           const std::vector<double> & exp, const std::vector<double> & eyp,
           const std::vector<double> & exm, const std::vector<double> & eym) {
    IDataPointSet * dset = create(path, title, 2);
    for ( int i = 0, N = y.size(); i < N; ++i ) dset->addPoint(DataPoint(2));
    if ( !dset->setCoordinate(0, x, exp, exm) ||
         !dset->setCoordinate(1, y, eyp, eym) )
      throw std::runtime_error("LWH could add points to DataPointSet '" +
                               title +  "'." );
    return dset;   
  }

  virtual IDataPointSet *
  createXY(const std::string & path, const std::string & title,
           const std::vector<double> & x, const std::vector<double> & y,
           const std::vector<double> & ex, const std::vector<double> & ey) {
    return createXY(path, title, x, y, ex, ey, ex, ey);
  }

  virtual IDataPointSet *
  createXY(const std::string & pathAndTitle,
           const std::vector<double> & x, const std::vector<double> & y,
           const std::vector<double> & exp, const std::vector<double> & eyp,
           const std::vector<double> & exm, const std::vector<double> & eym) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createXY(pathAndTitle, title, x, y, exp, eyp, exm, eym);
  }

  virtual IDataPointSet *
  createXY(const std::string & pathAndTitle,
           const std::vector<double> & x, const std::vector<double> & y,
           const std::vector<double> & ex, const std::vector<double> & ey) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createXY(pathAndTitle, title, x, y, ex, ey, ex, ey);
  }

  virtual IDataPointSet *
  createXYZ(const std::string & path, const std::string & title,
            const std::vector<double> & x, const std::vector<double> & y,
            const std::vector<double> & z, const std::vector<double> & exp,
            const std::vector<double> & eyp, const std::vector<double> & ezp,
            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));
    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 '" +
                               title +  "'." );
    return dset;   
  }

  virtual IDataPointSet *
  createXYZ(const std::string & path, const std::string & title,
            const std::vector<double> & x, const std::vector<double> & y,
            const std::vector<double> & z, const std::vector<double> & ex,
            const std::vector<double> & ey, const std::vector<double> & ez) {
    return createXYZ(path, title, x, y, z, ex, ey, ez, ex, ey, ez);
  }

  virtual IDataPointSet *
  createXYZ(const std::string & pathAndTitle, const std::vector<double> & x,
            const std::vector<double> & y, const std::vector<double> & z,
            const std::vector<double> & exp, const std::vector<double> & eyp,
            const std::vector<double> & ezp, const std::vector<double> & exm,
            const std::vector<double> & eym, const std::vector<double> & ezm) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createXYZ(pathAndTitle, title, x, y, z,
                     exp, eyp, ezp, exm, eym, ezm);
  }

  virtual IDataPointSet *
  createXYZ(const std::string & pathAndTitle, const std::vector<double> & x,
            const std::vector<double> & y, const std::vector<double> & z,
            const std::vector<double> & ex, const std::vector<double> & ey,
            const std::vector<double> & ez) {
    std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
    return createXYZ(pathAndTitle, title, x, y, z, ex, ey, ez, ex, ey, ez);
  }

  virtual IDataPointSet *
  createCopy(const std::string & path, const IDataPointSet & dataPointSet) {
    IDataPointSet * dset =
      create(path, dataPointSet.title(), dataPointSet.dimension());
    for ( int i = 0, N = dataPointSet.size(); i < N; ++i )
      dset->addPoint(*dataPointSet.point(i));
    return dset;
  }

  virtual bool destroy(IDataPointSet * dataPointSet) {
    IManagedObject * mo = dynamic_cast<IManagedObject *>(dataPointSet);
    if ( !mo ) return false;
    return tree->rm(tree->findPath(*mo));
  }

  virtual IDataPointSet *
  create(const std::string & path, const IHistogram1D & hist,
         const std::string & = "") {
    IDataPointSet * dset = create(path, hist.title(), 2);
    std::vector<double> x, y, ex, ey;
    for ( int i = 2, N = hist.axis().bins() + 2; i < N; ++i ) {
      dset->addPoint(DataPoint(2));
      x.push_back(hist.binMean(i - 2));
      ex.push_back(hist.axis().binWidth(i - 2));
      y.push_back(hist.binHeight(i - 2));
      ey.push_back(hist.binError(i - 2));
    }
    if ( !dset->setCoordinate(0, x, ex, ex) ||
         !dset->setCoordinate(1, y, ey, ey) )
      throw std::runtime_error("LWH could add points to DataPointSet '" +
                               hist.title() +  "'." );
    return dset;
  }

  virtual IDataPointSet *
  create(const std::string & path, const IHistogram2D & hist,
         const std::string & = "") {
    IDataPointSet * dset = create(path, hist.title(), 3);
    
    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);
        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);
        const double binwidth = hist.xAxis().binWidth(ix - 2)*
          hist.yAxis().binWidth(iy - 2);
        z.push_back(hist.binHeight(ix - 2, iy - 2)/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;
  }


  virtual IDataPointSet * create(const std::string &, const IHistogram3D &,
                                 const std::string & = "") {
    return error<IDataPointSet>("IHistogram3D");
  }

  virtual IDataPointSet * create(const std::string &, const ICloud1D &,
                                 const std::string & = "") {
    return error<IDataPointSet>("ICloud1D");
  }

  virtual IDataPointSet * create(const std::string &, const ICloud2D &,
                                 const std::string & = "") {
    return error<IDataPointSet>("ICloud2D");
  }

  virtual IDataPointSet * create(const std::string &, const ICloud3D &,
                                 const std::string & = "") {
    return error<IDataPointSet>("ICloud3D");
  }

  virtual IDataPointSet * create(const std::string &, const IProfile1D &,
                                 const std::string & = "") {
    return error<IDataPointSet>("IProfile1D");
  }

  virtual IDataPointSet * create(const std::string &, const IProfile2D &,
                                 const std::string & = "") {
    return error<IDataPointSet>("IProfile2D");
  }

  virtual IDataPointSet * add(const std::string &,
                              const IDataPointSet &, const IDataPointSet &) {
    return error<IDataPointSet>("addition of data points");
  }

  virtual IDataPointSet * subtract(const std::string &, const IDataPointSet &,
                                   const IDataPointSet &) {
    return error<IDataPointSet>("subtraction of data points");
  }

  virtual IDataPointSet * multiply(const std::string &, const IDataPointSet &,
                                   const IDataPointSet &) {
    return error<IDataPointSet>("multiplication of data points");
  }

  virtual IDataPointSet * divide(const std::string &, const IDataPointSet &,
                                 const IDataPointSet &) {
    return error<IDataPointSet>("division of data points");
  }

  virtual IDataPointSet *
  weightedMean(const std::string &, const IDataPointSet &,
               const IDataPointSet &) {
    return error<IDataPointSet>("weighted means of data points");
  }
  
private:

  template <typename T>
  static T * error(std::string feature) {
    throw std::runtime_error("LWH cannot handle " + feature + ".");
  }

  Tree * tree;

};
}

#endif