Class describing a Vavilov distribution.
The probability density function of the Vavilov distribution as function of Landau's parameter is given by:
\[ 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\]
where \(\phi(s) = e^{C} e^{\psi(s)}\) with \( C = \kappa (1+\beta^2 \gamma )\) and \(\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}}\). \( \gamma = 0.5772156649\dots\) is Euler's constant.
For the class VavilovAccurate, Pdf returns the Vavilov distribution as function of Landau's parameter \(\lambda_L = \lambda_V/\kappa - \ln \kappa\), 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 Google books
Therefore, for small values of \(\kappa < 0.01\), pdf approaches the Landau distribution.
For values \(\kappa > 10\), the Gauss approximation should be used with \(\mu\) and \(\sigma\) 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, Computer Phys. Comm. 7 (1974) 215-224, which has been implemented in CERNLIB (G116).
The class stores coefficients needed to calculate \(p(\lambda; \kappa, \beta^2)\) for fixed values of \(\kappa\) and \(\beta^2\). Changing these values is computationally expensive.
The parameter \(\kappa\) should be in the range \(0.01 \le \kappa \le 10\). In contrast to the CERNLIB implementation, all values of \(\kappa \ge 0.001\) may be used, but may result in slower running and/or inaccurate results.
The parameter \(\beta^2\) must be in the range \(0 \le \beta^2 \le 1\).
Two parameters which are fixed in the CERNLIB implementation may be set by the user:
- epsilonPM corresponds to \(\epsilon^+ = \epsilon^-\) in Eqs. (2.1) and (2.2) of Schorr's paper. epsilonPM gives an estimate on the integral of the cumulative distribution function outside the range \(\lambda_{min} \le \lambda \le \lambda_{max}\) 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 \(\epsilon^+ = \epsilon^- = 5\cdot 10^{-4}\). The code from CERNLIB has been extended such that also smaller values are possible.
- epsilon corresponds to \(\epsilon\) in Eq. (4.10) of Schorr's paper. It determines the accuracy of the series expansion. Schorr recommends \(\epsilon = 10^{-5}\).
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 pre-calculated 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
Definition at line 131 of file VavilovAccurate.h.