#ifndef __JROOT__JROOTFIT__
#define __JROOT__JROOTFIT__

#include <vector>
#include <set>
#include <limits>
#include <algorithm>

#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TF1.h"
#include "TF2.h"
#include "TF3.h"
#include "TList.h"

#include "JFit/JGandalf.hh" 
#include "JMath/JMathSupportkit.hh"
#include "JMath/JMathlib.hh"
#include "JTools/JRange.hh"


/**
 * \author mdejong
 */

namespace JROOT {}
namespace JPP { using namespace JROOT; }

namespace JROOT {

  using JFIT::JGandalf;
  using JMATH::parameter_list;
  using JMATH::poisson;
  using JMATH::getNumberOfParameters;
  using JMATH::getParameter;
  using JMATH::setParameters;


  /**
   * Type definiton of abscissa range.
   */
  typedef JTOOLS::JRange<double>  range_type;


  /**
   * Get range of given axis.
   *
   * \param  pAxis         pointer to axis
   * \return               range
   */
  inline range_type getRange(TAxis* pAxis)
  {
    if (pAxis != NULL)
      return range_type(pAxis->GetXmin(), pAxis->GetXmax());
    else
      return range_type();
  }


  /**
   * Data point for counting measurement.
   */
  struct m_count {
    /**
     * Default constructor.
     */
    m_count() :
      count(0)
    {}


    /**
     * Constructor.
     *
     * \param  count          count
     */
    m_count(const size_t count) :
      count(count)
    {}


    /**
     * Minimal value for numerical computations.
     *
     * \return                epsilon
     */
    static inline double epsilon()
    {
      return std::numeric_limits<double>::min();
    }


    /**
     * Get chi2.
     *
     * \param  z              expectation value
     * \return                chi2
     */
    inline double getRho(const double z) const
    {
      const double P = poisson(count, z > epsilon() ? z : epsilon());

      return -log(P > epsilon() ? P : epsilon());
    }


    /**
     * Get derivative of chi2.
     *
     * \param  z              expectation value
     * \return                derivative of chi2
     */
    inline double getPsi(const double z) const
    {
      return 1.0 - count/(z > epsilon() ? z : epsilon());
    }


    size_t count;             //!< count
  };


  /**
   * Data point for value with error.
   */
  struct m_value {
    /**
     * Default constructor.
     */
    m_value() :
      value(0.0),
      error(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  value          value
     * \param  error          error
     */
    m_value(const double value, 
	    const double error) :
      value(value),
      error (error)
    {}


    /**
     * Get chi2.
     *
     * \param  z              expectation value
     * \return                chi2
     */
    inline double getRho(const double z) const
    {
      const double u = (value - z) / error;

      return u*u;
    }


    /**
     * Get derivative of chi2.
     *
     * \param  z              expectation value
     * \return                derivative of chi2
     */
    inline double getPsi(const double z) const
    {
      const double u = (value - z) / error;

      return -0.5 * u / error;
    }

    
    double value;             //!< value
    double error;             //!< error
  };


  /**
   * 1D data point.
   */
  template<class T>
  struct m_1d :
    public T
  {
    /**
     * Default constructor.
     */
    m_1d() :
      T(),
      x(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  x              abscissa
     * \param  v              ordinate
     */
    m_1d(const double x, const T& v) :
      T(v),
      x(x)
    {}
    
    double x;                 //!< abscissa
  };


  /**
   * 2D data point.
   */
  template<class T>
  struct m_2d :
    public T
  {
    /**
     * Default constructor.
     */
    m_2d() :
      T(),
      x(0.0),
      y(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  x              abscissa
     * \param  y              abscissa
     * \param  v              ordinate
     */
    m_2d(const double x, const double y, const T& v) :
      T(v),
      x(x),
      y(y)
    {}
    
    double x;                 //!< abscissa
    double y;                 //!< abscissa
  };


  /**
   * 3D data point.
   */
  template<class T>
  struct m_3d :
    public T
  {
    /**
     * Default constructor.
     */
    m_3d() :
      T(),
      x(0.0),
      y(0.0),
      z(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  x              abscissa
     * \param  y              abscissa
     * \param  z              abscissa
     * \param  v              ordinate
     */
    m_3d(const double x, const double y, const double z, const T& v) :
      T(v),
      x(x),
      y(y),
      z(z)
    {}
    
    double x;                 //!< abscissa
    double y;                 //!< abscissa
    double z;                 //!< abscissa
  };


  /**
   * Template definition of data structure for set of data points.
   */
  template<class T>
  struct data_type;


  /**
   * Template specialisation of data structure for set of 1D data points.
   */
  template<class T>
  struct data_type< m_1d<T> > :
    public std::vector< m_1d<T> >
  {
    /**
     * Default constructor.
     */
    data_type()
    {}


    /**
     * Unpack constructor.
     *
     * \param  h1            1D histogram
     * \param  X             abscissa range
     */
    data_type(const TH1& h1,
	      const range_type& X = range_type())
    {
      unpack(h1, X);
    }


    /**
     * Unpack 1D-histogram.
     *
     * \param  h1            histogram
     * \param  X             abscissa range
     */
    void unpack(const TH1& h1,
		const range_type& X = range_type());
  };

  
  /**
   * Unpack 1D-histogram.
   *
   * \param  h1            histogram
   * \param  X             abscissa range
   */
  template<>
  void data_type< m_1d<m_count> >::unpack(const TH1& h1,
					  const range_type& X)
  {
    for (Int_t ix = 1; ix <= h1.GetXaxis()->GetNbins(); ++ix) {

      const double x     = h1.GetXaxis()->GetBinCenter(ix);
      const size_t count = h1.GetBinContent(ix);

      if (X(x)) {
	this->push_back(m_1d<m_count>(x, count));
      }
    }
  }

  
  /**
   * Unpack 1D-histogram.
   *
   * \param  h1            histogram
   * \param  X             abscissa range
   */
  template<>
  void data_type< m_1d<m_value> >::unpack(const TH1& h1,
					  const range_type& X)
  {
    for (Int_t ix = 1; ix <= h1.GetXaxis()->GetNbins(); ++ix) {

      const double x     = h1.GetXaxis()->GetBinCenter(ix);
      const double value = h1.GetBinContent(ix);
      const double error = h1.GetBinError  (ix);

      if (X(x) && error > 0.0) {
	this->push_back(m_1d<m_value>(x, m_value(value, error)));
      }
    }
  }


  /**
   * Template specialisation of data structure for set of 2D data points.
   */
  template<class T>
  struct data_type< m_2d<T> > :
    public std::vector< m_2d<T> >
  {
    /**
     * Default constructor.
     */
    data_type()
    {}


    /**
     * Unpack constructor.
     *
     * \param  h2            2D histogram
     * \param  X             abscissa range
     * \param  Y             abscissa range
     */
    data_type(const TH2& h2,
	      const range_type& X = range_type(),
	      const range_type& Y = range_type())
    {
      unpack(h2, X, Y);
    }


    /**
     * Unpack 2D-histogram.
     *
     * \param  h2            histogram
     * \param  X             abscissa range
     * \param  Y             abscissa range
     */
    void unpack(const TH2& h2, 
		const range_type& X = range_type(),
		const range_type& Y = range_type());
  };

  
  /**
   * Unpack 2D-histogram.
   *
   * \param  h2            histogram
   * \param  X             abscissa range
   * \param  Y             abscissa range
   */
  template<>
  void data_type< m_2d<m_count> >::unpack(const TH2& h2,
					  const range_type& X,
					  const range_type& Y)
  {
    for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
      for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {

	const double x     = h2.GetXaxis()->GetBinCenter(ix);
	const double y     = h2.GetYaxis()->GetBinCenter(iy);
	const size_t count = h2.GetBinContent(ix,iy);

	if (X(x) && Y(y)) {
	  this->push_back(m_2d<m_count>(x, y, count));
	}
      }
    }
  }

  
  /**
   * Unpack 2D-histogram.
   *
   * \param  h2            histogram
   * \param  X             abscissa range
   * \param  Y             abscissa range
   */
  template<>
  void data_type< m_2d<m_value> >::unpack(const TH2& h2,
					  const range_type& X,
					  const range_type& Y)
  {
    for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
      for (Int_t iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {

	const double x     = h2.GetXaxis()->GetBinCenter(ix);
	const double y     = h2.GetYaxis()->GetBinCenter(iy);
	const double value = h2.GetBinContent(ix,iy);
	const double error = h2.GetBinError  (ix,iy);

	if (X(x) && Y(y) && error > 0.0) {
	  this->push_back(m_2d<m_value>(x, y, m_value(value, error)));
	}
      }
    }
  }


  /**
   * Template specialisation of data structure for set of 3D data points.
   */
  template<class T>
  struct data_type< m_3d<T> > :
    public std::vector< m_3d<T> >
  {
    /**
     * Default constructor.
     */
    data_type()
    {}


    /**
     * Unpack constructor.
     *
     * \param  h3            2D histogram
     * \param  X             abscissa range
     * \param  Y             abscissa range
     * \param  Z             abscissa range
     */
    data_type(const TH3& h3,
	      const range_type& X = range_type(),
	      const range_type& Y = range_type(),
	      const range_type& Z = range_type())
    {
      unpack(h3, X, Y, Z);
    }


    /**
     * Unpack 3D-histogram.
     *
     * \param  h3            histogram
     * \param  X             abscissa range
     * \param  Y             abscissa range
     * \param  Z             abscissa range
     */
    void unpack(const TH3& h3, 
		const range_type& X = range_type(),
		const range_type& Y = range_type(),
		const range_type& Z = range_type());
  };

  
  /**
   * Unpack 3D-histogram.
   *
   * \param  h3            histogram
   * \param  X             abscissa range
   * \param  Y             abscissa range
   * \param  Z             abscissa range
   */
  template<>
  void data_type< m_3d<m_count> >::unpack(const TH3& h3,
					  const range_type& X,
					  const range_type& Y,
					  const range_type& Z)
  {
    for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
      for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {

	  const double x     = h3.GetXaxis()->GetBinCenter(ix);
	  const double y     = h3.GetYaxis()->GetBinCenter(iy);
	  const double z     = h3.GetZaxis()->GetBinCenter(iz);
	  const size_t count = h3.GetBinContent(ix,iy,iz);
	  
	  if (X(x) && Y(y) && Z(z)) {
	    this->push_back(m_3d<m_count>(x, y, z, count));
	  }
	}
      }
    }
  }

  
  /**
   * Unpack 3D-histogram.
   *
   * \param  h3            histogram
   * \param  X             abscissa range
   * \param  Y             abscissa range
   * \param  Z             abscissa range
   */
  template<>
  void data_type< m_3d<m_value> >::unpack(const TH3& h3,
					  const range_type& X,
					  const range_type& Y,
					  const range_type& Z)
  {
    for (Int_t ix = 1; ix <= h3.GetXaxis()->GetNbins(); ++ix) {
      for (Int_t iy = 1; iy <= h3.GetYaxis()->GetNbins(); ++iy) {
	for (Int_t iz = 1; iz <= h3.GetZaxis()->GetNbins(); ++iz) {

	  const double x     = h3.GetXaxis()->GetBinCenter(ix);
	  const double y     = h3.GetYaxis()->GetBinCenter(iy);
	  const double z     = h3.GetZaxis()->GetBinCenter(iz);
	  const double value = h3.GetBinContent(ix,iy,iz);
	  const double error = h3.GetBinError  (ix,iy,iz);

	  if (X(x) && Y(y) && Z(z) && error > 0.0) {
	    this->push_back(m_3d<m_value>(x, y, z, m_value(value, error)));
	  }
	}
      }
    }
  }


  /**
   * Wrapper data structure to build ROOT 1D function.
   */
  struct JF1 :
    public TF1
  {
    /**
     * Constructor.
     *
     * \param  name          name
     * \param  f1            function
     * \param  X             fit range
     */
    template<class JF1_t>
    JF1(const char* const name,
	const JF1_t& f1,
	const range_type& X  = range_type()) :
      TF1(name,
	  helper<JF1_t>(f1),
	  X.getLowerLimit(),
	  X.getUpperLimit(),
	  0)
    {
      this->SetNpx(1000);
    };


    /**
     * Helper.
     */
    template<class JF1_t>
    struct helper : 
      public JF1_t
    {
      /**
       * Constructor.
       *
       * \param  f1            function
       */
      helper(const JF1_t& f1) :
	JF1_t(f1)
      {}


      /**
       * ROOT compatible function call.
       *
       * \param  x             pointer to abscissa  values
       * \param  parameters    pointer to parameter values
       * \return               function value
       */
      double operator()(const double* x, const double* parameters)
      {
	setParameters(this, parameters);

	return this->getValue(x[0]);
      }
    };
  };


  /**
   * Wrapper data structure to build ROOT 2D function.
   */
  struct JF2 :
    public TF2
  {
    /**
     * Constructor.
     *
     * \param  name          name
     * \param  f2            function
     * \param  X             fit range
     * \param  Y             fit range
     */
    template<class JF2_t>
    JF2(const char* const name,
	const JF2_t& f2,
	const range_type& X  = range_type(),
	const range_type& Y  = range_type()) :
      TF2(name,
	  helper<JF2_t>(f2),
	  X.getLowerLimit(),
	  X.getUpperLimit(),
	  Y.getLowerLimit(),
	  Y.getUpperLimit(),
	  0)
    {
      this->SetNpx(1000);
      this->SetNpy(1000);
    };


    /**
     * Helper.
     */
    template<class JF2_t>
    struct helper :
      public JF2_t
    {
      /**
       * Constructor.
       *
       * \param  f2            function
       */
      helper(const JF2_t& f2) :
	JF2_t(f2)
      {}


      /**
       * ROOT compatible function call.
       *
       * \param  x             pointer to abscissa  values
       * \param  parameters    pointer to parameter values
       * \return               function value
       */
      double operator()(const double* x, const double* parameters)
      {
	setParameters(this, parameters);

	return this->getValue(x[0], x[1]);
      }
    };
  };


  /**
   * Wrapper data structure to build ROOT 3D function.
   */
  struct JF3 :
    public TF3
  {
    /**
     * Constructor.
     *
     * \param  name          name
     * \param  f3            function
     * \param  X             fit range
     * \param  Y             fit range
     * \param  Z             fit range
     */
    template<class JF3_t>
    JF3(const char* const name,
	const JF3_t& f3,
	const range_type& X  = range_type(),
	const range_type& Y  = range_type(),
	const range_type& Z  = range_type()) :
      TF3(name,
	  helper<JF3_t>(f3),
	  X.getLowerLimit(),
	  X.getUpperLimit(),
	  Y.getLowerLimit(),
	  Y.getUpperLimit(),
	  Z.getLowerLimit(),
	  Z.getUpperLimit(),
	  0)
    {
      this->SetNpx(300);
      this->SetNpy(300);
      this->SetNpz(300);
    };


    /**
     * Helper.
     */
    template<class JF3_t>
    struct helper :
      public JF3_t
    {
      /**
       * Constructor.
       *
       * \param  f3            function
       */
      helper(const JF3_t& f3) :
	JF3_t(f3)
      {}


      /**
       * ROOT compatible function call.
       *
       * \param  x             pointer to abscissa  values
       * \param  parameters    pointer to parameter values
       * \return               function value
       */
      double operator()(const double* x, const double* parameters)
      {
	setParameters(this, parameters);

	return this->getValue(x[0], x[1], x[2]);
      }
    };
  };


  /**
   * Auxiliary data structure for list of fixed parameters.
   */
  struct index_list :
    public std::set<size_t>
  {
    /**
     * Default constructor.
     */
    index_list()
    {}


    /**
     * Constructor.
     *
     * \param  indices       indices
     */
    index_list(const std::initializer_list<size_t>& indices) :
      std::set<size_t>(indices)
    {}


    /**
     * Conversion constructor.
     *
     * \param  parameters    parameters
     */
    template<class T>
    index_list(const std::initializer_list<double T::*>& parameters)
    {
      for (size_t i = 0; i != T::parameters.size(); ++i) {
        if (std::find(parameters.begin(), parameters.end(), T::parameters[i]) != parameters.end()) {
          this->insert(i);
        }
      }
    }
  };
  

  /**
   * Base class for result of ROOT Fit.
   */
  template<class JFs_t>
  class JRootfit_t :
    public JGandalf<JFs_t>
  {
  protected:
    /**
     * Default constructor.
     */
    JRootfit_t() :
      npx (0),
      chi2(0.0)
    {}


    size_t npx;                            //!< number of data points 
    double chi2;                           //!< chi2

  public:
    /**
     * Get function.
     *
     * \return               function.
     */
    const JFs_t& getFunction() const
    {
      return this->value;
    }


    /**
     * Get number of parameters.
     *
     * \return               number of parameters
     */
    size_t getNumberOfParameters() const
    {
      return JFs_t::parameters.size();
    }


    /**
     * Get number of free parameters.
     *
     * \return               number of free parameters
     */
    size_t getNumberOfFreeParameters() const
    {
      return this->parameters.size();
    }


    /**
     * Get number of data points.
     *
     * \return               number of data points
     */
    size_t getN() const
    {
      return npx;
    }


    /**
     * Get chi2.
     *
     * \return               chi2
     */
    double getChi2() const
    {
      return chi2;
    }


    /**
     * Get number of degrees of freedom.
     *
     * \return               number of degrees of freedom
     */
    int getNDF() const
    {
      return (int) getN() - (int) getNumberOfFreeParameters();
    }


    /**
     * Get value of parameter at given index.
     *
     * \param  i             index
     */
    double getValue(size_t i) const
    {
      return getParameter(this->value, i);
    }


    /**
     * Get error of parameter at given index.
     *
     * \param  i             index
     */
    double getError(size_t i) const
    {
      return getParameter(this->error, i);
    }
  };


  /**
   * ROOT Fit.
   */
  template<class JFs_t>
  class JRootfit :
    public JRootfit_t<JFs_t>
  {
  public:

    typedef JRootfit_t<JFs_t>  result_type;


    /**
     * Default constructor.
     */
    JRootfit()
    {}


    /**
     * Fit.
     *
     * \param  h1            histogram
     * \param  f1            start value
     * \param  type          type of data for histogram unpacking
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \return               result
     */
    template<class T>
    const result_type& operator()(const TH1&        h1, 
				  const JFs_t&      f1,
				  const T&          type,
				  const index_list& ls = index_list(),
				  const range_type& X  = range_type())
    {
      return eval(f1, ls, data_type< m_1d<T> >(h1, X));
    }


    /**
     * Fit.
     *
     * The fitted function is added to the input histogram.
     *
     * \param  h1            pointer to histogram
     * \param  f1            start value
     * \param  type          type of data for histogram unpacking
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \return               result
     */
    template<class T>
    const result_type& operator()(TH1*              h1,
				  const JFs_t&      f1,
				  const T&          type,
				  const index_list& ls = index_list(),
				  const range_type& X  = range_type())
    {
      (*this)(*h1, f1, type, ls, X);

      h1->GetListOfFunctions()->Add(new JF1("f1",
					    this->value, 
					    JTOOLS::join(X, getRange(h1->GetXaxis()))));

      return static_cast<const result_type&>(*this);
    }


    /**
     * Fit.
     *
     * \param  h2            histogram
     * \param  f2            start value
     * \param  type          type of data for histogram unpacking
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \param  Y             fit range
     * \return               result
     */
    template<class T>
    const result_type& operator()(const TH2&        h2, 
				  const JFs_t&      f2,
				  const T&          type,
				  const index_list& ls = index_list(),
				  const range_type& X  = range_type(),
				  const range_type& Y  = range_type())
    {
      return eval(f2, ls, data_type< m_2d<T> >(h2, X, Y));
    }


    /**
     * Fit.
     *
     * The fitted function is added to the input histogram.
     *
     * \param  h2            pointer to histogram
     * \param  f2            start value
     * \param  type          type of data for histogram unpacking
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \param  Y             fit range
     * \return               result
     */
    template<class T>
    const result_type& operator()(TH2*              h2,
				  const JFs_t&      f2,
				  const T&          type,
				  const index_list& ls = index_list(),
				  const range_type& X  = range_type(),
				  const range_type& Y  = range_type())
    {
      (*this)(*h2, f2, type, ls, X, Y);

      h2->GetListOfFunctions()->Add(new JF2("f2",
					    this->value, 
					    JTOOLS::join(X, getRange(h2->GetXaxis())),
					    JTOOLS::join(Y, getRange(h2->GetYaxis()))));

      return static_cast<const result_type&>(*this);
    }


    /**
     * Fit.
     *
     * \param  h3            histogram
     * \param  f3            start value
     * \param  type          type of data for histogram unpacking
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \param  Y             fit range
     * \param  Z             fit range
     * \return               result
     */
    template<class T>
    const result_type& operator()(const TH3&        h3, 
				  const JFs_t&      f3,
				  const T&          type,
				  const index_list& ls = index_list(),
				  const range_type& X  = range_type(),
				  const range_type& Y  = range_type(),
				  const range_type& Z  = range_type())
    {
      return eval(f3, ls, data_type< m_3d<T> >(h3, X, Y, Z));
    }


    /**
     * Fit.
     *
     * The fitted function is added to the input histogram.
     *
     * \param  h3            pointer to histogram
     * \param  f3            start value
     * \param  type          type of data for histogram unpacking
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \param  Y             fit range
     * \param  Z             fit range
     * \return               result
     */
    template<class T>
    const result_type& operator()(TH3*              h3,
				  const JFs_t&      f3,
				  const T&          type,
				  const index_list& ls = index_list(),
				  const range_type& X  = range_type(),
				  const range_type& Y  = range_type(),
				  const range_type& Z  = range_type())
    {
      (*this)(*h3, f3, type, ls, X, Y, Z);

      h3->GetListOfFunctions()->Add(new JF3("f3",
					    this->value, 
					    JTOOLS::join(X, getRange(h3->GetXaxis())),
					    JTOOLS::join(Y, getRange(h3->GetYaxis())),
					    JTOOLS::join(Z, getRange(h3->GetZaxis()))));

      return static_cast<const result_type&>(*this);
    }


    static JRootfit Fit;        //!< Global fit object

  private:
    size_t getNumberOfFreeParameters();          // hide method
    size_t getN();                               // hide method
    double getChi2();                            // hide method
    int    getNDF();                             // hide method


    /**
     * Evaluate fit.
     *
     * \param  fs            start value
     * \param  ls            list of fixed parameters
     * \param  data          data
     * \return               result
     */
    template<class T>
    const result_type& eval(const JFs_t&        fs,
			    const index_list&   ls,
			    const data_type<T>& data)
    {
      this->parameters.clear();

      for (size_t i = 0; i != JFs_t::parameters.size(); ++i) {
	if (ls.count(i) == 0) {
	  this->parameters.push_back(JFs_t::parameters[i]);
	}
      }

      this->value = fs;
      this->npx   = data.size();
      this->chi2  = static_cast<JGandalf<JFs_t>&>(*this)(this->fit, data.begin(), data.end()).chi2;

      return static_cast<const result_type&>(*this);
    }


    /**
     * Auxiliary data structure for fit functions.
     */
    const struct function_type {

      typedef typename JGandalf<JFs_t>::result_type  result_type;

      /**
       * Fit function.
       *
       * \param  fs             function
       * \param  mp             data point
       * \return                chi2 and gradient
       */
      template<class T>
      inline const result_type& operator()(const JFs_t& fs, const m_1d<T>& mp) const
      {
	const double y   = fs.getValue(mp.x);

	result.chi2      = mp.getRho(y);
	result.gradient  = fs.getGradient(mp.x);
	result.gradient *= mp.getPsi(y);

	return result;
      }


      /**
       * Fit function.
       *
       * \param  fs             function
       * \param  mp             data point
       * \return                chi2 and gradient
       */
      template<class T>
      inline const result_type& operator()(const JFs_t& fs, const m_2d<T>& mp) const
      {
	const double y   = fs.getValue(mp.x, mp.y);

	result.chi2      = mp.getRho(y);
	result.gradient  = fs.getGradient(mp.x, mp.y);
	result.gradient *= mp.getPsi(y);

	return result;
      }


      /**
       * Fit function.
       *
       * \param  fs             function
       * \param  mp             data point
       * \return                chi2 and gradient
       */
      template<class T>
      inline const result_type& operator()(const JFs_t& fs, const m_3d<T>& mp) const
      {
	const double y   = fs.getValue(mp.x, mp.y, mp.z);

	result.chi2      = mp.getRho(y);
	result.gradient  = fs.getGradient(mp.x, mp.y, mp.z);
	result.gradient *= mp.getPsi(y);

	return result;
      }

    private:
      mutable result_type result; 
    } fit;
  };


  /**
   * Global fit object.
   */
  template<class JFs_t>
  JRootfit<JFs_t> JRootfit<JFs_t>::Fit;


  /**
   * Global fit fuction.
   *
   * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
   *
   * \param  h1            histogram
   * \param  f1            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \return               result
   */
  template<class T, class JF1_t>
  inline JRootfit_t<JF1_t> Fit(const TH1& h1,
			       const JF1_t& f1,
			       const index_list& ls = index_list(),
			       const range_type& X  = range_type())
  {
    return JRootfit<JF1_t>::Fit(h1, f1, T(), ls, X);
  }


  /**
   * Global fit fuction.
   *
   * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
   * The fitted function is added to the input histogram.
   *
   * \param  h1            pointer to histogram
   * \param  f1            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \return               result
   */
  template<class T, class JF1_t>
  inline JRootfit_t<JF1_t> Fit(TH1* h1,
			       const JF1_t& f1,
			       const index_list& ls = index_list(),
			       const range_type& X  = range_type())
  {
    return JRootfit<JF1_t>::Fit(h1, f1, T(), ls, X);
  }


  /**
   * Global fit fuction.
   *
   * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
   *
   * \param  h2            histogram
   * \param  f2            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \param  Y             fit range
   * \return               result
   */
  template<class T, class JF2_t>
  inline JRootfit_t<JF2_t> Fit(const TH2& h2,
			       const JF2_t& f2,
			       const index_list& ls = index_list(),
			       const range_type& X  = range_type(),
			       const range_type& Y  = range_type())
  {
    return JRootfit<JF2_t>::Fit(h2, f2, T(), ls, X, Y);
  }


  /**
   * Global fit fuction.
   *
   * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
   * The fitted function is added to the input histogram.
   *
   * \param  h2            pointer to histogram
   * \param  f2            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \param  Y             fit range
   * \return               result
   */
  template<class T, class JF2_t>
  inline JRootfit_t<JF2_t> Fit(TH2* h2,
			       const JF2_t& f2,
			       const index_list& ls = index_list(),
			       const range_type& X  = range_type(),
			       const range_type& Y  = range_type())
  {
    return JRootfit<JF2_t>::Fit(h2, f2, T(), ls, X, Y);
  }


  /**
   * Global fit fuction.
   *
   * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
   *
   * \param  h3            histogram
   * \param  f3            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \param  Y             fit range
   * \param  Z             fit range
   * \return               result
   */
  template<class T, class JF3_t>
  inline JRootfit_t<JF3_t> Fit(const TH3& h3,
			       const JF3_t& f3,
			       const index_list& ls = index_list(),
			       const range_type& X  = range_type(),
			       const range_type& Y  = range_type(),
			       const range_type& Z  = range_type())
  {
    return JRootfit<JF3_t>::Fit(h3, f3, T(), ls, X, Y, Z);
  }


  /**
   * Global fit fuction.
   *
   * The template parameter <tt>T</tt> refers to the measurement and can be <tt>m_count</tt> or <tt>m_value</tt>.\n
   * The fitted function is added to the input histogram.
   *
   * \param  h3            pointer to histogram
   * \param  f3            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \param  Y             fit range
   * \param  Z             fit range
   * \return               result
   */
  template<class T, class JF3_t>
  inline JRootfit_t<JF3_t> Fit(TH3* h3,
			       const JF3_t& f3,
			       const index_list& ls = index_list(),
			       const range_type& X  = range_type(),
			       const range_type& Y  = range_type(),
			       const range_type& Z  = range_type())
  {
    return JRootfit<JF3_t>::Fit(h3, f3, T(), ls, X, Y, Z);
  }
}

#endif