# 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,
Nucl. Sci. Series Report No. 39,
Washington (Nat. Akad. Sci.) 1964, 388 pp.
Available from

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 = 5.0000000000000001E-4, double epsilon = 1.0000000000000001E-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::Vavilov ROOT::Math::Vavilov::Vavilov() ROOT::Math::Vavilov ROOT::Math::Vavilov::Vavilov(const ROOT::Math::Vavilov&) ROOT::Math::VavilovAccurate VavilovAccurate(const ROOT::Math::VavilovAccurate&) ROOT::Math::VavilovAccurate VavilovAccurate(double kappa = 1, double beta2 = 1, double epsilonPM = 5.0000000000000001E-4, double epsilon = 1.0000000000000001E-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, double(ROOT::Math::VavilovAccurate::*)(double)const f) const

## Data Members

private:
 static ROOT::Math::VavilovAccurate::(anonymous) MAXTERMS 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 static ROOT::Math::VavilovAccurate::(anonymous) kNquantMax

## Function documentation

VavilovAccurate(double kappa = 1, double beta2 = 1, double epsilonPM = 5.0000000000000001E-4, double epsilon = 1.0000000000000001E-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 = 5.0000000000000001E-4, double epsilon = 1.0000000000000001E-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



      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, double(ROOT::Math::VavilovAccurate::*)(double)const f) const