JGauss.hh

Go to the documentation of this file.

#ifndef __JMATH__JGAUSS__
#define __JMATH__JGAUSS__
 
#include <istream>
#include <ostream>
#include <iomanip>
#include <cmath>
#include <limits>
 
#include "JLang/JEquals.hh"
#include "JLang/JManip.hh"
#include "JMath/JMath.hh"
#include "JMath/JConstants.hh"
 
/**
 * \author mdejong
 */
 
namespace JMATH {}
namespace JPP { using namespace JMATH; }
 
namespace JMATH {
 
  using JLANG::JEquals;
  using JLANG::JValueOutOfRange;
 
  /**
   * Gauss model.
   */
  struct JGauss_t :
    public JMath<JGauss_t>,
    public JEquals<JGauss_t>
  {
    /**
     * Default constructor.
     */
    JGauss_t() :
      mean      (0.0),
      sigma     (0.0),
      signal    (0.0),
      background(0.0)
    {}
 
 
    /**
     * Constructor.
     * 
     * Note that an error will be thrown for a sigma smaller or equal to zero.
     *
     * \param  mean          mean
     * \param  sigma         sigma
     * \param  signal        signal
     * \param  background    background
     */
    JGauss_t(const double mean,
             const double sigma,
             const double signal,
             const double background) :
      mean      (mean),
      sigma     (sigma),
      signal    (signal),
      background(background)
    {}
 
 
    /**
     * Equality.
     *
     * \param  gauss         gauss
     * \param  eps           numerical precision
     * \return               true if gauss's identical; else false
     */
    bool equals(const JGauss_t& gauss,
                const double    eps = std::numeric_limits<double>::min()) const
    {
      return (fabs(mean       - gauss.mean)       <=  eps  &&
              fabs(sigma      - gauss.sigma)      <=  eps  &&
              fabs(signal     - gauss.signal)     <=  eps  &&
              fabs(background - gauss.background) <=  eps);
    }
 
 
    /**
     * Add gauss.
     *
     * \param  gauss         gauss
     * \return               this gauss
     */
    JGauss_t& add(const JGauss_t& gauss)
    {
      mean       += gauss.mean;
      sigma      += gauss.sigma;
      signal     += gauss.signal;
      background += gauss.background;
 
      return *this;
    }
 
 
    /**
     * Subtract gauss.
     *
     * \param  gauss         gauss
     * \return               this gauss
     */
    JGauss_t& sub(const JGauss_t& gauss)
    {
      mean       -= gauss.mean;
      sigma      -= gauss.sigma;
      signal     -= gauss.signal;
      background -= gauss.background;
 
      return *this;
    }
 
 
    /**
     * Scale gauss.
     *
     * \param  factor        multiplication factor
     * \return               this gauss
     */
    JGauss_t& mul(const double factor)
    {
      mean       *= factor;
      sigma      *= factor;
      signal     *= factor;
      background *= factor;
 
      return *this;
    }
 
 
    /**
     * Write Gauss to input stream.
     *
     * \param  in            input stream
     * \param  gauss         gauss
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JGauss_t& gauss) 
    {
      return in >> gauss.mean >> gauss.sigma >> gauss.signal >> gauss.background;
    }
 
 
    /**
     * Write Gauss to output stream.
     *
     * \param  out           output stream
     * \param  gauss         gauss
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JGauss_t& gauss) 
    {
      using namespace std;
 
      return out << FIXED(7,3) << gauss.mean   << ' '
                 << FIXED(7,3) << gauss.sigma  << ' '
                 << FIXED(9,3) << gauss.signal << ' '
                 << FIXED(9,3) << gauss.background;
    }
 
    double mean;
    double sigma;
    double signal;
    double background;
  };
 
 
  /**
   * Gauss function object.
   *
   * Evaluates function, derivative and gradient values.
   */
  struct JGauss :
    public JGauss_t
  {
    /**
     * Type definition of fit parameter.
     */
    typedef double JGauss_t::*parameter_type;
 
 
    /**
     * Default constructor.
     */
    JGauss() :
      JGauss_t()
    {}
 
 
    /**
     * Copy constructor.
     *
     * \param  gauss         gauss
     */
    JGauss(const JGauss_t& gauss) :
      JGauss_t(gauss)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  mean          mean
     * \param  sigma         sigma
     * \param  signal        signal
     * \param  background    background
     */
    JGauss(const double mean,
           const double sigma,
           const double signal     = 1.0,
           const double background = 0.0) :
      JGauss_t(mean, sigma, signal, background)
    {}
 
 
    /**
     * Function value.
     * 
     * Note that an error will be thrown if sigma <= 0.0.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double getValue(const double x) const
    {
      if (sigma <= 0.0) { THROW(JValueOutOfRange, "JGauss::getValue: Invalid sigma " << sigma); }
      
      const double u  = (x - mean) / sigma;
 
      return signal * get(u) + background;
    }
 
 
    /**
     * Derivative value.
     * 
     * Note that an error will be thrown if sigma <= 0.0.
     *
     * \param  x             abscissa value
     * \return               derivative value
     */
    double getDerivative(const double x) const
    {
      if (sigma <= 0.0) { THROW(JValueOutOfRange, "JGauss::getDerivative: Invalid sigma " << sigma); }
 
      const double u  = (x - mean) / sigma;
 
      return signal * get(u) * -u / sigma;
    }
 
 
    /**
     * Integral value.
     * 
     * Note that an error will be thrown if sigma <= 0.0.
     * Note that zero will be returned if xmin is equal to xmax.
     * Note that the integration limits can be set to infinity using std::numeric_limits<double>::inifinity().
     * 
     * \param  xmin          lower limit
     * \param  xmax          upper limit 
     */
    double getIntegral(const double xmin, const double xmax) const
    {
      if (sigma <= 0.0) { THROW(JValueOutOfRange, "JGauss::getIntegral: Invalid sigma " << sigma); }
 
      double zmin = xmin;
      double zmax = xmax;
 
      if (xmax < xmin) {
 
        zmin = xmax;
        zmax = xmin;
      }
        
      if (zmin < zmax) {
 
        const double umin = (zmin - mean) / sigma / sqrt(2.0);
        const double umax = (zmax - mean) / sigma / sqrt(2.0);
        
        double I;
 
        if        (std::isinf(zmin) && std::isinf(zmax)) {
 
          I = 1.0;
        
        } else if (std::isinf(zmin)) {
 
          I = 0.5 * (1 + erf(umax));
 
        } else if (std::isinf(zmax)) {
 
          I = 0.5 * erfc(umin);
 
        } else {
 
          I = 0.5 * (erfc(umin) - erfc(umax));
        }
 
        return signal * I + background;
      
      } else {
 
        return 0.0;
      }
    }
 
 
    /**
     * Function value.
     *
     * \param  x             abscissa value
     * \return               function value
     */
    double operator()(const double x) const
    {
      return getValue(x);
    }
 
 
    /**
     * Get gradient.
     * 
     * Note that an error will be thrown if sigma <= 0.0.
     *
     * \param  x             abscissa value
     * \return               gradient
     */
    const JGauss_t& getGradient(const double x) const
    {
      if (sigma <= 0.0) { THROW(JValueOutOfRange, "JGauss::getGradient: Invalid sigma " << sigma); }
 
      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 * (u)                   * w;             // d(f)/d(mean)
      gradient.sigma      =  fs * (u + 1.0) * (u - 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:
    /**
     * Get ordinate value.
     * 
     * Note that an error will be thrown if sigma <= 0.0.
     *
     * \param  u             abscissa value
     * \return               ordinate value
     */
    inline double get(const double u) const
    {
      if (sigma <= 0.0) { THROW(JValueOutOfRange, "JGauss::get: Invalid sigma " << sigma); }
 
      return exp(-0.5*u*u) / (sqrt(2.0*PI) * sigma);
    }
 
    mutable JGauss_t gradient;
  };
}
 
#endif