ROOT logo
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 doubleCdf(double x) const
virtual doubleCdf(double x, double kappa, double beta2)
virtual doubleCdf_c(double x) const
virtual doubleCdf_c(double x, double kappa, double beta2)
virtual doubleGetBeta2() const
doubleGetEpsilon() const
doubleGetEpsilonPM() const
static ROOT::Math::VavilovAccurate*GetInstance()
static ROOT::Math::VavilovAccurate*GetInstance(double kappa, double beta2)
virtual doubleGetKappa() const
virtual doubleGetLambdaMax() const
virtual doubleGetLambdaMin() const
doubleGetNTerms() const
virtual doubleROOT::Math::Vavilov::Kurtosis() const
static doubleROOT::Math::Vavilov::Kurtosis(double kappa, double beta2)
virtual doubleROOT::Math::Vavilov::Mean() const
static doubleROOT::Math::Vavilov::Mean(double kappa, double beta2)
virtual doubleMode() const
virtual doubleMode(double kappa, double beta2)
ROOT::Math::VavilovAccurate&operator=(const ROOT::Math::VavilovAccurate&)
virtual doublePdf(double x) const
virtual doublePdf(double x, double kappa, double beta2)
virtual doubleQuantile(double z) const
virtual doubleQuantile(double z, double kappa, double beta2)
virtual doubleQuantile_c(double z) const
virtual doubleQuantile_c(double z, double kappa, double beta2)
voidSet(double kappa, double beta2, double epsilonPM = 5E-4, double epsilon = 1E-5)
virtual voidSetKappaBeta2(double kappa, double beta2)
virtual doubleROOT::Math::Vavilov::Skewness() const
static doubleROOT::Math::Vavilov::Skewness(double kappa, double beta2)
virtual doubleROOT::Math::Vavilov::Variance() const
static doubleROOT::Math::Vavilov::Variance(double kappa, double beta2)
ROOT::Math::VavilovAccurateVavilovAccurate(const ROOT::Math::VavilovAccurate&)
ROOT::Math::VavilovAccurateVavilovAccurate(double kappa = 1, double beta2 = 1, double epsilonPM = 5E-4, double epsilon = 1E-5)
private:
static doubleE1plLog(double x)
doubleG116f1(double x) const
doubleG116f2(double x) const
voidInitQuantile() const
intRzero(double a, double b, double& x0, double eps, int mxf, void f) const

Data Members

public:
enum { MAXTERMS
kNquantMax
};
private:
doublefA_cdf[501]
doublefA_pdf[501]
doublefB_cdf[501]
doublefB_pdf[501]
doublefBeta2
doublefEpsilon
doublefEpsilonPM
doublefH[8]
doublefKappa
doublefLambda[32]
intfNQuant
doublefOmega
doublefQuant[32]
boolfQuantileInit
doublefT
doublefT0
doublefT1
doublefX0
static ROOT::Math::VavilovAccurate*fgInstance

Class Charts

Inheritance Inherited Members Includes Libraries
Class Charts

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