JBell.hh

Go to the documentation of this file.

#ifndef __JMATH__JBELL__
#define __JMATH__JBELL__
 
#include <cmath>
#include <limits>
 
#include "JMath/JGauss.hh"
 
/**
 * \author mdejong
 */
 
namespace JMATH {}
namespace JPP { using namespace JMATH; }
 
namespace JMATH {
 
  /**
   * Type definition for Bell model.
   */
  typedef JGauss_t  JBell_t;
 
 
  /**
   * Bell function object.\n
   * This function object correspods to a Gauss with a different power (<tt>alpha</tt>) for the exponential term. 
   *
   * Evaluates function, derivative and gradient values.
   */
  struct JBell :
    public JBell_t
  {
    /**
     * Type definition of fit parameter.
     */
    typedef double JBell_t::*parameter_type;
 
 
    /**
     * Default constructor.
     */
    JBell() :
      JBell_t()
    {}
 
 
    /**
     * Copy constructor.
     *
     * \param  bell          bell
     * \param  alpha         alpha
     */
    JBell(const JBell_t& bell,
          const double   alpha    = 2.0) :
      JBell_t(bell),
      alpha(alpha),
      V(2.0 * pow(2.0,1.0/alpha) * tgamma(1.0/alpha) / alpha)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  mean          mean
     * \param  sigma         sigma
     * \param  signal        signal
     * \param  background    background
     * \param  alpha         alpha
     */
    JBell(const double mean,
          const double sigma,
          const double signal     = 1.0,
          const double background = 0.0,
          const double alpha      = 2.0) :
      JBell(JBell_t(mean, sigma, signal, background), alpha)
    {}
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      const double u  = (x - mean) / sigma;
 
      return signal * get(u) + background;
    }
 
 
    /**
     * Derivative value.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      const double u  = (x - mean) / sigma;
 
      return signal * get(u) * copysign(_f(u, alpha - 1.0), -u) * alpha / sigma;
    }
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double operator()(const double x) const
    {
      return getValue(x);
    }
 
 
    /**
     * Get gradient.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    const JBell_t& getGradient(const double x) const
    {
      const double w  = 1.0 / sigma;
      const double u  = (x - mean) * w;
      const double f0 =  get(u);
      const double fs =  signal * f0;
 
      gradient.mean       =  fs * copysign(_f(u, alpha - 1.0), u) * alpha * w;      // d(f)/d(mean)
      gradient.sigma      =  fs *        -(_f(u, alpha) * alpha + 1.0)    * w;      // d(f)/d(sigma)
      gradient.signal     =  f0;                                                    // d(f)/d(signal)
      gradient.background =  1.0;                                                   // d(f)/d(background)
 
      return gradient;
    }
 
  private:
    double alpha;            //!< power
    double V;                //!< normalisation constant
    
    /**
     * Power.
     *
     * \param  x             abscissa value
     * \param  a             power
     * \return               ordinate value
     */
    static inline double _f(const double x, const double a)
    {
      return -0.5 * pow(fabs(x), a);
    }
 
 
    /**
     * Get ordinate value.
     *
     * \param  u             abscissa value
     * \return               ordinate value
     */
    inline double get(const double u) const
    {
      const double v = _f(u, alpha);
 
      return exp(v) / (sigma * V);
    }
 
    mutable JBell_t gradient;
  };
}
 
#endif