JRootfit.hh

Go to the documentation of this file.

#ifndef __JROOT__JROOTFIT__
#define __JROOT__JROOTFIT__
 
#include <ostream>
#include <iomanip>
#include <vector>
#include <set>
#include <limits>
#include <algorithm>
#include <memory>
 
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TF1.h"
#include "TF2.h"
#include "TF3.h"
#include "TGraphErrors.h"
#include "TGraph2DErrors.h"
#include "TList.h"
 
#include "JFit/JGandalf.hh" 
#include "JFit/JMEstimator.hh"
#include "JMath/JMathSupportkit.hh"
#include "JMath/JMathlib.hh"
#include "JTools/JRange.hh"
#include "JLang/JManip.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JROOT {}
namespace JPP { using namespace JROOT; }
 
namespace JROOT {
 
  using JFIT::JGandalf;
  using JFIT::JMEstimator;
  using JMATH::parameter_list;
  using JMATH::poisson;
  using JMATH::getNumberOfParameters;
  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();
  }
 
 
  /**
   * Get range of given graph.
   *
   * \param  nx            number of values
   * \param  pX            pointer to array of values
   * \return               range
   */
  inline range_type getRange(const Int_t nx, const Double_t* pX)
  {
    range_type range = range_type::DEFAULT_RANGE();
 
    for (Int_t i = 0; i != nx; ++i) {
      range.include(pX[i]);
    }
 
    return range;
  }
 
 
  /**
   * Get range of given graph.
   *
   * \param  g1            graph
   * \return               range
   */
  inline range_type getRange(const TGraph& g1)
  {
    return getRange(g1.GetN(), g1.GetX());
  }
 
 
  /**
   * Get range of given graph.
   *
   * \param  g2            graph
   * \return               range
   */
  inline range_type getRangeX(const TGraph2D& g2)
  {
    return getRange(g2.GetN(), g2.GetX());
  }
 
 
  /**
   * Get range of given graph.
   *
   * \param  g2            graph
   * \return               range
   */
  inline range_type getRangeY(const TGraph2D& g2)
  {
    return getRange(g2.GetN(), g2.GetY());
  }
 
 
  /**
   * 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 {
 
    typedef std::shared_ptr<JMEstimator>  estimator_type;
 
    /**
     * 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 = (z - value) / error;
 
      return getMEstimator()->getRho(u) / getMEstimator()->getRho(1.0);
    }
 
 
    /**
     * Get derivative of chi2.
     *
     * \param  z              expectation value
     * \return                derivative of chi2
     */
    inline double getPsi(const double z) const
    {
      const double u = (z - value) / error;
 
      return getMEstimator()->getPsi(u) * getMEstimator()->getRho(1.0);
    }
 
    
    double value;             //!< value
    double error;             //!< error
 
 
    /**
     * Get M-estimator.
     *
     * \return                M-estimator
     */
    static const estimator_type& getMEstimator()
    {
      return get_mestimator();
    }
 
 
    /**
     * Set M-estimator.
     *
     * \param  type           M-estimator type
     */
    static const void setMEstimator(int type)
    {
      get_mestimator().reset(JFIT::getMEstimator(type));
    }
 
  private:
    /**
     * Get M-estimator.
     *
     * \return                M-estimator
     */
    static estimator_type& get_mestimator()
    {
      static estimator_type estimator(new JFIT::JMEstimatorNormal());
 
      return estimator;
    }
  };
 
 
  /**
   * 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 constructor.
     *
     * \param  g1            1D graph
     * \param  X             abscissa range
     */
    data_type(const TGraphErrors& g1,
              const range_type& X = range_type())
    {
      unpack(g1, X);
    }
 
 
    /**
     * Unpack 1D-histogram.
     *
     * \param  h1            histogram
     * \param  X             abscissa range
     */
    void unpack(const TH1& h1,
                const range_type& X = range_type());
 
 
    /**
     * Unpack 1D-graph.
     *
     * \param  g1            graph
     * \param  X             abscissa range
     */
    void unpack(const TGraphErrors& g1,
                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)));
      }
    }
  }
 
  
  /**
   * Unpack 1D-graph
   *
   * \param  g1            graph
   * \param  X             abscissa range
   */
  template<>
  void data_type< m_1d<m_value> >::unpack(const TGraphErrors& g1,
                                          const range_type& X)
  {
    for (Int_t i = 0; i != g1.GetN(); ++i) {
 
      const double x     = g1.GetX()[i];
      const double value = g1.GetY()[i];
      const double error = g1.GetEY()[i];
 
      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 constructor.
     *
     * \param  g2            graph
     * \param  X             abscissa range
     * \param  Y             abscissa range
     */
    data_type(const TGraph2DErrors& g2,
              const range_type& X = range_type(),
              const range_type& Y = range_type())
    {
      unpack(g2, 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-graph.
     *
     * \param  g2            graph
     * \param  X             abscissa range
     * \param  Y             abscissa range
     */
    void unpack(const TGraph2DErrors& g2,
                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)));
        }
      }
    }
  }
 
  
  /**
   * Unpack 2D-graph
   *
   * \param  g2            graph
   * \param  X             abscissa range
   * \param  Y             abscissa range
   */
  template<>
  void data_type< m_2d<m_value> >::unpack(const TGraph2DErrors& g2,
                                          const range_type& X,
                                          const range_type& Y)
  {
    for (Int_t i = 0; i != g2.GetN(); ++i) {
 
      const double x     = g2.GetX()[i];
      const double y     = g2.GetY()[i];
      const double value = g2.GetZ()[i];
      const double error = g2.GetEZ()[i];
 
      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) :
      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,
        const range_type& Y) :
      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,
        const range_type& Y,
        const range_type& Z) :
      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)
    {
      JGandalf<JFs_t>::EPSILON_ABSOLUTE    =  true;
      JGandalf<JFs_t>::EPSILON             =  1.0e-5;
      JGandalf<JFs_t>::MAXIMUM_ITERATIONS  =  1000000;
    }
 
 
    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 (this->value)[i];
    }
 
 
    /**
     * Get error of parameter at given index.
     *
     * \param  i             index
     */
    double getError(size_t i) const
    {
      return (this->error)[i];
    }
 
 
    /**
     * Print result.
     *
     * \param  out           output stream
     */
    void print(std::ostream& out) const
    {
      using namespace std;
 
      out << "chi2/NDF " << FIXED(7,3) << this->getChi2() << "/" << this->getNDF() << endl;
 
      out << "Number of iterations " << this->numberOfIterations << endl;
 
      for (size_t i = 0; i != this->getNumberOfParameters(); ++i) {
        out << setw(2)     << i                 << ' '
            << FIXED(15,9) << this->getValue(i) << " +/- "
            << FIXED(15,9) << this->getError(i) << endl;
      }
    }
  };
 
 
  /**
   * 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  g1            graph
     * \param  f1            start value
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \return               result
     */
    const result_type& operator()(const TGraphErrors& g1, 
                                  const JFs_t&        f1,
                                  const index_list&   ls = index_list(),
                                  const range_type&   X  = range_type())
    {
      return eval(f1, ls, data_type< m_1d<m_value> >(g1, X));
    }
 
 
    /**
     * Fit.
     *
     * \param  g1            pointer to graph
     * \param  f1            start value
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \return               result
     */
    const result_type& operator()(TGraphErrors*       g1, 
                                  const JFs_t&        f1,
                                  const index_list&   ls = index_list(),
                                  const range_type&   X  = range_type())
    {
      (*this)(*g1, f1, ls, X);
 
      g1->GetListOfFunctions()->Add(new JF1("f1",
                                            this->value, 
                                            JTOOLS::join(X, getRange(*g1))));
 
      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  g2            graph
     * \param  f2            start value
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \param  Y             fit range
     * \return               result
     */
    const result_type& operator()(const TGraph2DErrors& g2, 
                                  const JFs_t&          f2,
                                  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<m_value> >(g2, X, Y));
    }
 
 
    /**
     * Fit.
     *
     * \param  g2            pointer to graph
     * \param  f2            start value
     * \param  ls            list of fixed parameters
     * \param  X             fit range
     * \param  Y             fit range
     * \return               result
     */
    const result_type& operator()(TGraph2DErrors*     g2, 
                                  const JFs_t&        f2,
                                  const index_list&   ls = index_list(),
                                  const range_type&   X  = range_type(),
                                  const range_type&   Y  = range_type())
    {
      (*this)(*g2, f2, ls, X, Y);
 
      g2->GetListOfFunctions()->Add(new JF2("f2",
                                            this->value, 
                                            JTOOLS::join(X, getRangeX(*g2)),
                                            JTOOLS::join(Y, getRangeY(*g2))));
 
      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.
   *
   * \param  g1            graph
   * \param  f1            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \return               result
   */
  template<class JF1_t>
  inline JRootfit_t<JF1_t> Fit(const TGraphErrors& g1,
                               const JF1_t& f1,
                               const index_list& ls = index_list(),
                               const range_type& X  = range_type())
  {
    return JRootfit<JF1_t>::Fit(g1, f1, ls, X);
  }
 
 
  /**
   * Global fit fuction.
   *
   * \param  g1            pointer to graph
   * \param  f1            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \return               result
   */
  template<class JF1_t>
  inline JRootfit_t<JF1_t> Fit(TGraphErrors* g1,
                               const JF1_t& f1,
                               const index_list& ls = index_list(),
                               const range_type& X  = range_type())
  {
    return JRootfit<JF1_t>::Fit(g1, f1, 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.
   *
   * \param  g2            graph
   * \param  f2            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \param  Y             fit range
   * \return               result
   */
  template<class JF2_t>
  inline JRootfit_t<JF2_t> Fit(const TGraph2DErrors& g2,
                               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(g2, f2, ls, X, Y);
  }
 
 
  /**
   * Global fit fuction.
   *
   * \param  g2            pointer to graph
   * \param  f2            start value
   * \param  ls            list of fixed parameters
   * \param  X             fit range
   * \param  Y             fit range
   * \return               result
   */
  template<class JF2_t>
  inline JRootfit_t<JF2_t> Fit(TGraph2DErrors* g2,
                               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(g2, f2, 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