Quick Links: ROOT Homepage Class Index Class Hierarchy

Sections: class description function members data members class charts

ROOT » MATH » MATHMORE » ROOT::Math::VavilovAccurate

class ROOT::Math::VavilovAccurate: public ROOT::Math::Vavilov


   Class describing a Vavilov distribution.

   The probability density function of the Vavilov distribution
   as function of Landau's parameter is given by:
  \f[ p(\lambda_L; \kappa, \beta^2) =
  \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} \phi(s) e^{\lambda_L s} ds\f]
   where \f$\phi(s) = e^{C} e^{\psi(s)}\f$
   with  \f$ C = \kappa (1+\beta^2 \gamma )\f$
   and \f$\psi(s)= s \ln \kappa + (s+\beta^2 \kappa)
               \cdot \left ( \int \limits_{0}^{1}
               \frac{1 - e^{\frac{-st}{\kappa}}}{t} \,d t- \gamma \right )
               - \kappa \, e^{\frac{-s}{\kappa}}\f$.
   \f$ \gamma = 0.5772156649\dots\f$ is Euler's constant.

   For the class VavilovAccurate,
   Pdf returns the Vavilov distribution as function of Landau's parameter
   \f$\lambda_L = \lambda_V/\kappa  - \ln \kappa\f$,
   which is the convention used in the CERNLIB routines, and in the tables
   by S.M. Seltzer and M.J. Berger: Energy loss stragglin of protons and mesons:
   Tabulation of the Vavilov distribution, pp 187-203
   in: National Research Council (U.S.), Committee on Nuclear Science:
   Studies in penetration of charged particles in matter,
   Nat. Akad. Sci. Publication 1133,
   Nucl. Sci. Series Report No. 39,
   Washington (Nat. Akad. Sci.) 1964, 388 pp.
   Available from
   <A HREF="http://books.google.de/books?id=kmMrAAAAYAAJ&lpg=PP9&pg=PA187#v=onepage&q&f=false">Google books</A>

   Therefore, for small values of \f$\kappa < 0.01\f$,
   pdf approaches the Landau distribution.

   For values \f$\kappa > 10\f$, the Gauss approximation should be used
   with \f$\mu\f$ and \f$\sigma\f$ given by Vavilov::mean(kappa, beta2)
   and sqrt(Vavilov::variance(kappa, beta2).

   The original Vavilov pdf is obtained by
   v.Pdf(lambdaV/kappa-log(kappa))/kappa.

   For detailed description see
   B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
   <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>,
   which has been implemented in
   <A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/g116/top.html">
   CERNLIB (G116)</A>.

   The class stores coefficients needed to calculate \f$p(\lambda; \kappa, \beta^2)\f$
   for fixed values of \f$\kappa\f$ and \f$\beta^2\f$.
   Changing these values is computationally expensive.

   The parameter \f$\kappa\f$ should be in the range \f$0.01 \le \kappa \le 10\f$.
   In contrast to the CERNLIB implementation, all values of \f$\kappa \ge 0.001\f$ may be used,
   but may result in slower running and/or inaccurate results.

   The parameter \f$\beta^2\f$ must be in the range \f$0 \le \beta^2 \le 1\f$.

   Two parameters which are fixed in the CERNLIB implementation may be set by the user:
   - epsilonPM corresponds to \f$\epsilon^+ = \epsilon^-\f$ in Eqs. (2.1) and (2.2) of Schorr's paper.
   epsilonPM gives an estimate on the integral of the cummulative distribution function
   outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
   where the approximation is valid.
   Thus, it determines the support of the approximation used here (called $T_0 - T_1$ in the paper).
   Schorr recommends  \f$\epsilon^+ = \epsilon^- = 5\cdot 10^{-4}\f$.
   The code from CERNLIB has been extended such that also smaller values are possible.

   - epsilon corresponds to \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper.
   It determines the accuracy of the series expansion.
   Schorr recommends  \f$\epsilon = 10^{-5}\f$.

   For the quantile calculation, the algorithm given by Schorr is not used,
   because it turns out to be very slow and still inaccurate.
   Instead, an initial estimate is calculated based on a precalculated table,
   which is subsequently improved by Newton iterations.

   While the CERNLIB implementation calculates at most 156 terms in the series expansion
   for the pdf and cdf calculation, this class calculates up to 500 terms, depending
   on the values of epsilonPM and epsilon.

   Average times on a Pentium Core2 Duo P8400 2.26GHz:
   - 38us per call to SetKappaBeta2 or constructor
   - 0.49us per call to Pdf, Cdf
   - 8.2us per first call to Quantile after SetKappaBeta2 or constructor
   - 0.83us per subsequent call to Quantile

   Benno List, June 2010

   @ingroup StatFunc

Function Members (Methods)

public:

virtual	~VavilovAccurate()
virtual double	Cdf(double x) const
virtual double	Cdf(double x, double kappa, double beta2)
virtual double	Cdf_c(double x) const
virtual double	Cdf_c(double x, double kappa, double beta2)
virtual double	GetBeta2() const
double	GetEpsilon() const
double	GetEpsilonPM() const
static ROOT::Math::VavilovAccurate*	GetInstance()
static ROOT::Math::VavilovAccurate*	GetInstance(double kappa, double beta2)
virtual double	GetKappa() const
virtual double	GetLambdaMax() const
virtual double	GetLambdaMin() const
double	GetNTerms() const
virtual double	ROOT::Math::Vavilov::Kurtosis() const
static double	ROOT::Math::Vavilov::Kurtosis(double kappa, double beta2)
virtual double	ROOT::Math::Vavilov::Mean() const
static double	ROOT::Math::Vavilov::Mean(double kappa, double beta2)
virtual double	Mode() const
virtual double	Mode(double kappa, double beta2)
ROOT::Math::VavilovAccurate&	operator=(const ROOT::Math::VavilovAccurate&)
virtual double	Pdf(double x) const
virtual double	Pdf(double x, double kappa, double beta2)
virtual double	Quantile(double z) const
virtual double	Quantile(double z, double kappa, double beta2)
virtual double	Quantile_c(double z) const
virtual double	Quantile_c(double z, double kappa, double beta2)
void	Set(double kappa, double beta2, double epsilonPM = 5E-4, double epsilon = 1E-5)
virtual void	SetKappaBeta2(double kappa, double beta2)
virtual double	ROOT::Math::Vavilov::Skewness() const
static double	ROOT::Math::Vavilov::Skewness(double kappa, double beta2)
virtual double	ROOT::Math::Vavilov::Variance() const
static double	ROOT::Math::Vavilov::Variance(double kappa, double beta2)
ROOT::Math::VavilovAccurate	VavilovAccurate(const ROOT::Math::VavilovAccurate&)
ROOT::Math::VavilovAccurate	VavilovAccurate(double kappa = 1, double beta2 = 1, double epsilonPM = 5E-4, double epsilon = 1E-5)

private:

static double	E1plLog(double x)
double	G116f1(double x) const
double	G116f2(double x) const
void	InitQuantile() const
int	Rzero(double a, double b, double& x0, double eps, int mxf, void f) const

Data Members

public:

enum {	MAXTERMS
	kNquantMax
};

private:

double	fA_cdf[501]
double	fA_pdf[501]
double	fB_cdf[501]
double	fB_pdf[501]
double	fBeta2
double	fEpsilon
double	fEpsilonPM
double	fH[8]
double	fKappa
double	fLambda[32]
int	fNQuant
double	fOmega
double	fQuant[32]
bool	fQuantileInit
double	fT
double	fT0
double	fT1
double	fX0
static ROOT::Math::VavilovAccurate*	fgInstance

Class Charts

Inheritance Inherited Members Includes Libraries

Function documentation

VavilovAccurate(double kappa = 1, double beta2 = 1, double epsilonPM = 5E-4, double epsilon = 1E-5)

      Initialize an object to calculate the Vavilov distribution

       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
       @param epsilonPM: \f$\epsilon^+ = \epsilon^-\f$ in Eqs. (2.1) and (2.2) of Schorr's paper; gives an estimate on the integral of the cummulative distribution function
              outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
              where the approximation is valid.
       @param epsilon: \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper; determines the accuracy of the series expansion.

virtual ~VavilovAccurate()

     Destructor

double Pdf(double x) const

       Evaluate the Vavilov probability density function

       @param x The Landau parameter \f$x = \lambda_L\f$

double Pdf(double x, double kappa, double beta2)

       Evaluate the Vavilov probability density function,
       and set kappa and beta2, if necessary

       @param x The Landau parameter \f$x = \lambda_L\f$
       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

double Cdf(double x) const

       Evaluate the Vavilov cummulative probability density function

       @param x The Landau parameter \f$x = \lambda_L\f$

double Cdf(double x, double kappa, double beta2)

       Evaluate the Vavilov cummulative probability density function,
       and set kappa and beta2, if necessary

       @param x The Landau parameter \f$x = \lambda_L\f$
       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

double Cdf_c(double x) const

       Evaluate the Vavilov complementary cummulative probability density function

       @param x The Landau parameter \f$x = \lambda_L\f$

double Cdf_c(double x, double kappa, double beta2)

       Evaluate the Vavilov complementary cummulative probability density function,
       and set kappa and beta2, if necessary

       @param x The Landau parameter \f$x = \lambda_L\f$
       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

double Quantile(double z) const

       Evaluate the inverse of the Vavilov cummulative probability density function

       @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$

double Quantile(double z, double kappa, double beta2)

       Evaluate the inverse of the Vavilov cummulative probability density function,
       and set kappa and beta2, if necessary

       @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

double Quantile_c(double z) const

       Evaluate the inverse of the complementary Vavilov cummulative probability density function

       @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$

double Quantile_c(double z, double kappa, double beta2)

       Evaluate the inverse of the complementary Vavilov cummulative probability density function,
       and set kappa and beta2, if necessary

       @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

void SetKappaBeta2(double kappa, double beta2)

      Change \f$\kappa\f$ and \f$\beta^2\f$ and recalculate coefficients if necessary

       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

void Set(double kappa, double beta2, double epsilonPM = 5E-4, double epsilon = 1E-5)

      (Re)Initialize the object

       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
       @param epsilonPM \f$\epsilon^+ = \epsilon^-\f$ in Eqs. (2.1) and (2.2) of Schorr's paper; gives an estimate on the integral of the cummulative distribution function
              outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
              where the approximation is valid.
       @param epsilon \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper; determines the accuracy of the series expansion.

double GetLambdaMin() const

      Return the minimum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
      is nonzero in the current approximation

double GetLambdaMax() const

      Return the maximum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
      is nonzero in the current approximation

double GetKappa() const

      Return the current value of \f$\kappa\f$

double GetBeta2() const

      Return the current value of \f$\beta^2\f$

double Mode() const

      Return the value of \f$\lambda\f$ where the pdf is maximal

double Mode(double kappa, double beta2)

      Return the value of \f$\lambda\f$ where the pdf is maximal function,
       and set kappa and beta2, if necessary

       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

double GetEpsilonPM() const

      Return the current value of \f$\epsilon^+ = \epsilon^-\f$

double GetEpsilon() const

      Return the current value of \f$\epsilon\f$

double GetNTerms() const

      Return the number of terms used in the series expansion

VavilovAccurate * GetInstance()

      Returns a static instance of class VavilovFast

VavilovAccurate * GetInstance(double kappa, double beta2)

      Returns a static instance of class VavilovFast,
      and sets the values of kappa and beta2

       @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
       @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

void InitQuantile() const

double G116f1(double x) const

double G116f2(double x) const

int Rzero(double a, double b, double& x0, double eps, int mxf, void f) const