Logo ROOT   6.10/09
Reference Guide
VavilovAccurate.h
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: B. List 29.4.2010
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24 
25 // Header file for class VavilovAccurate
26 //
27 // Created by: blist at Thu Apr 29 11:19:00 2010
28 //
29 // Last update: Thu Apr 29 11:19:00 2010
30 //
31 #ifndef ROOT_Math_VavilovAccurate
32 #define ROOT_Math_VavilovAccurate
33 
34 
35 #include "Math/Vavilov.h"
36 
37 namespace ROOT {
38 namespace Math {
39 
40 //____________________________________________________________________________
41 /**
42  Class describing a Vavilov distribution.
43 
44  The probability density function of the Vavilov distribution
45  as function of Landau's parameter is given by:
46  \f[ p(\lambda_L; \kappa, \beta^2) =
47  \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} \phi(s) e^{\lambda_L s} ds\f]
48  where \f$\phi(s) = e^{C} e^{\psi(s)}\f$
49  with \f$ C = \kappa (1+\beta^2 \gamma )\f$
50  and \f$\psi(s)= s \ln \kappa + (s+\beta^2 \kappa)
51  \cdot \left ( \int \limits_{0}^{1}
52  \frac{1 - e^{\frac{-st}{\kappa}}}{t} \,d t- \gamma \right )
53  - \kappa \, e^{\frac{-s}{\kappa}}\f$.
54  \f$ \gamma = 0.5772156649\dots\f$ is Euler's constant.
55 
56  For the class VavilovAccurate,
57  Pdf returns the Vavilov distribution as function of Landau's parameter
58  \f$\lambda_L = \lambda_V/\kappa - \ln \kappa\f$,
59  which is the convention used in the CERNLIB routines, and in the tables
60  by S.M. Seltzer and M.J. Berger: Energy loss stragglin of protons and mesons:
61  Tabulation of the Vavilov distribution, pp 187-203
62  in: National Research Council (U.S.), Committee on Nuclear Science:
63  Studies in penetration of charged particles in matter,
64  Nat. Akad. Sci. Publication 1133,
65  Nucl. Sci. Series Report No. 39,
66  Washington (Nat. Akad. Sci.) 1964, 388 pp.
67  Available from
68  <A HREF="http://books.google.de/books?id=kmMrAAAAYAAJ&lpg=PP9&pg=PA187#v=onepage&q&f=false">Google books</A>
69 
70  Therefore, for small values of \f$\kappa < 0.01\f$,
71  pdf approaches the Landau distribution.
72 
73  For values \f$\kappa > 10\f$, the Gauss approximation should be used
74  with \f$\mu\f$ and \f$\sigma\f$ given by Vavilov::mean(kappa, beta2)
75  and sqrt(Vavilov::variance(kappa, beta2).
76 
77  The original Vavilov pdf is obtained by
78  v.Pdf(lambdaV/kappa-log(kappa))/kappa.
79 
80  For detailed description see
81  B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
82  <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>,
83  which has been implemented in
84  <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g116/top.html">
85  CERNLIB (G116)</A>.
86 
87  The class stores coefficients needed to calculate \f$p(\lambda; \kappa, \beta^2)\f$
88  for fixed values of \f$\kappa\f$ and \f$\beta^2\f$.
89  Changing these values is computationally expensive.
90 
91  The parameter \f$\kappa\f$ should be in the range \f$0.01 \le \kappa \le 10\f$.
92  In contrast to the CERNLIB implementation, all values of \f$\kappa \ge 0.001\f$ may be used,
93  but may result in slower running and/or inaccurate results.
94 
95  The parameter \f$\beta^2\f$ must be in the range \f$0 \le \beta^2 \le 1\f$.
96 
97  Two parameters which are fixed in the CERNLIB implementation may be set by the user:
98  - epsilonPM corresponds to \f$\epsilon^+ = \epsilon^-\f$ in Eqs. (2.1) and (2.2) of Schorr's paper.
99  epsilonPM gives an estimate on the integral of the cumulative distribution function
100  outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
101  where the approximation is valid.
102  Thus, it determines the support of the approximation used here (called $T_0 - T_1$ in the paper).
103  Schorr recommends \f$\epsilon^+ = \epsilon^- = 5\cdot 10^{-4}\f$.
104  The code from CERNLIB has been extended such that also smaller values are possible.
105 
106  - epsilon corresponds to \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper.
107  It determines the accuracy of the series expansion.
108  Schorr recommends \f$\epsilon = 10^{-5}\f$.
109 
110  For the quantile calculation, the algorithm given by Schorr is not used,
111  because it turns out to be very slow and still inaccurate.
112  Instead, an initial estimate is calculated based on a pre-calculated table,
113  which is subsequently improved by Newton iterations.
114 
115  While the CERNLIB implementation calculates at most 156 terms in the series expansion
116  for the pdf and cdf calculation, this class calculates up to 500 terms, depending
117  on the values of epsilonPM and epsilon.
118 
119  Average times on a Pentium Core2 Duo P8400 2.26GHz:
120  - 38us per call to SetKappaBeta2 or constructor
121  - 0.49us per call to Pdf, Cdf
122  - 8.2us per first call to Quantile after SetKappaBeta2 or constructor
123  - 0.83us per subsequent call to Quantile
124 
125  Benno List, June 2010
126 
127  @ingroup StatFunc
128  */
129 
130 
131 class VavilovAccurate: public Vavilov {
132 
133 public:
134 
135 
136  /**
137  Initialize an object to calculate the Vavilov distribution
138 
139  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
140  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
141  @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 cumulative distribution function
142  outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
143  where the approximation is valid.
144  @param epsilon: \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper; determines the accuracy of the series expansion.
145  */
146 
147  VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5);
148 
149  /**
150  Destructor
151  */
152  virtual ~VavilovAccurate();
153 
154 
155 public:
156 
157  /**
158  Evaluate the Vavilov probability density function
159 
160  @param x The Landau parameter \f$x = \lambda_L\f$
161  */
162  double Pdf (double x) const;
163 
164  /**
165  Evaluate the Vavilov probability density function,
166  and set kappa and beta2, if necessary
167 
168  @param x The Landau parameter \f$x = \lambda_L\f$
169  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
170  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
171  */
172  double Pdf (double x, double kappa, double beta2);
173 
174  /**
175  Evaluate the Vavilov cumulative probability density function
176 
177  @param x The Landau parameter \f$x = \lambda_L\f$
178  */
179  double Cdf (double x) const;
180 
181  /**
182  Evaluate the Vavilov cumulative probability density function,
183  and set kappa and beta2, if necessary
184 
185  @param x The Landau parameter \f$x = \lambda_L\f$
186  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
187  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
188  */
189  double Cdf (double x, double kappa, double beta2);
190 
191  /**
192  Evaluate the Vavilov complementary cumulative probability density function
193 
194  @param x The Landau parameter \f$x = \lambda_L\f$
195  */
196  double Cdf_c (double x) const;
197 
198  /**
199  Evaluate the Vavilov complementary cumulative probability density function,
200  and set kappa and beta2, if necessary
201 
202  @param x The Landau parameter \f$x = \lambda_L\f$
203  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
204  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
205  */
206  double Cdf_c (double x, double kappa, double beta2);
207 
208  /**
209  Evaluate the inverse of the Vavilov cumulative probability density function
210 
211  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
212  */
213  double Quantile (double z) const;
214 
215  /**
216  Evaluate the inverse of the Vavilov cumulative probability density function,
217  and set kappa and beta2, if necessary
218 
219  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
220  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
221  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
222  */
223  double Quantile (double z, double kappa, double beta2);
224 
225  /**
226  Evaluate the inverse of the complementary Vavilov cumulative probability density function
227 
228  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
229  */
230  double Quantile_c (double z) const;
231 
232  /**
233  Evaluate the inverse of the complementary Vavilov cumulative probability density function,
234  and set kappa and beta2, if necessary
235 
236  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
237  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
238  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
239  */
240  double Quantile_c (double z, double kappa, double beta2);
241 
242  /**
243  Change \f$\kappa\f$ and \f$\beta^2\f$ and recalculate coefficients if necessary
244 
245  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
246  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
247  */
248  virtual void SetKappaBeta2 (double kappa, double beta2);
249 
250 
251  /**
252  (Re)Initialize the object
253 
254  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
255  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
256  @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 cumulative distribution function
257  outside the range \f$\lambda_{min} \le \lambda \le \lambda_{max}\f$
258  where the approximation is valid.
259  @param epsilon \f$\epsilon\f$ in Eq. (4.10) of Schorr's paper; determines the accuracy of the series expansion.
260  */
261  void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5);
262 
263 
264  /**
265  Return the minimum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
266  is nonzero in the current approximation
267  */
268  virtual double GetLambdaMin() const;
269 
270  /**
271  Return the maximum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
272  is nonzero in the current approximation
273  */
274  virtual double GetLambdaMax() const;
275 
276  /**
277  Return the current value of \f$\kappa\f$
278  */
279  virtual double GetKappa() const;
280 
281  /**
282  Return the current value of \f$\beta^2\f$
283  */
284  virtual double GetBeta2() const;
285 
286  /**
287  Return the value of \f$\lambda\f$ where the pdf is maximal
288  */
289  virtual double Mode() const;
290 
291  /**
292  Return the value of \f$\lambda\f$ where the pdf is maximal function,
293  and set kappa and beta2, if necessary
294 
295  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
296  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
297  */
298  virtual double Mode(double kappa, double beta2);
299 
300  /**
301  Return the current value of \f$\epsilon^+ = \epsilon^-\f$
302  */
303 
304  double GetEpsilonPM() const;
305 
306  /**
307  Return the current value of \f$\epsilon\f$
308  */
309  double GetEpsilon() const;
310 
311  /**
312  Return the number of terms used in the series expansion
313  */
314  double GetNTerms() const;
315 
316  /**
317  Returns a static instance of class VavilovFast
318  */
319  static VavilovAccurate *GetInstance();
320 
321  /**
322  Returns a static instance of class VavilovFast,
323  and sets the values of kappa and beta2
324 
325  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
326  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
327  */
328  static VavilovAccurate *GetInstance(double kappa, double beta2);
329 
330 
331 private:
332  enum{MAXTERMS=500};
334  double fKappa, fBeta2;
336 
337  mutable bool fQuantileInit;
338  mutable int fNQuant;
339  enum{kNquantMax=32};
340  mutable double fQuant[kNquantMax];
341  mutable double fLambda[kNquantMax];
342 
343  void InitQuantile() const;
344 
346 
347  double G116f1 (double x) const;
348  double G116f2 (double x) const;
349 
350  int Rzero (double a, double b, double& x0,
351  double eps, int mxf, double (VavilovAccurate::*f)(double)const) const;
352  static double E1plLog (double x); // Calculates log(|x|)+E_1(x)
353 
354 };
355 
356  /**
357  The Vavilov probability density function
358 
359  @param x The Landau parameter \f$x = \lambda_L\f$
360  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
361  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
362 
363  @ingroup PdfFunc
364  */
365 double vavilov_accurate_pdf (double x, double kappa, double beta2);
366 
367  /**
368  The Vavilov cumulative probability density function
369 
370  @param x The Landau parameter \f$x = \lambda_L\f$
371  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
372  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
373 
374  @ingroup ProbFunc
375  */
376 double vavilov_accurate_cdf (double x, double kappa, double beta2);
377 
378  /**
379  The Vavilov complementary cumulative probability density function
380 
381  @param x The Landau parameter \f$x = \lambda_L\f$
382  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
383  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
384 
385  @ingroup ProbFunc
386  */
387 double vavilov_accurate_cdf_c (double x, double kappa, double beta2);
388 
389  /**
390  The inverse of the Vavilov cumulative probability density function
391 
392  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
393  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
394  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
395 
396  @ingroup QuantFunc
397  */
398 double vavilov_accurate_quantile (double z, double kappa, double beta2);
399 
400  /**
401  The inverse of the complementary Vavilov cumulative probability density function
402 
403  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
404  @param kappa The parameter \f$\kappa\f$, which must be in the range \f$\kappa \ge 0.001 \f$
405  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
406 
407  @ingroup QuantFunc
408  */
409 double vavilov_accurate_quantile_c (double z, double kappa, double beta2);
410 
411 } // namespace Math
412 } // namespace ROOT
413 
414 #endif /* ROOT_Math_VavilovAccurate */
Base class describing a Vavilov distribution.
Definition: Vavilov.h:120
virtual double GetKappa() const
Return the current value of .
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
TArc * a
Definition: textangle.C:12
int Rzero(double a, double b, double &x0, double eps, int mxf, double(VavilovAccurate::*f)(double) const) const
double G116f2(double x) const
double Cdf(double x) const
Evaluate the Vavilov cumulative probability density function.
double GetNTerms() const
Return the number of terms used in the series expansion.
double vavilov_accurate_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
double GetEpsilonPM() const
Return the current value of .
Double_t x[n]
Definition: legend1.C:17
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
double GetEpsilon() const
Return the current value of .
void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5)
(Re)Initialize the object
static VavilovAccurate * GetInstance()
Returns a static instance of class VavilovFast.
double Quantile(double z) const
Evaluate the inverse of the Vavilov cumulative probability density function.
double Cdf_c(double x) const
Evaluate the Vavilov complementary cumulative probability density function.
virtual double Mode() const
Return the value of where the pdf is maximal.
static VavilovAccurate * fgInstance
VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5)
Initialize an object to calculate the Vavilov distribution.
REAL epsilon
Definition: triangle.c:617
constexpr Double_t E()
Definition: TMath.h:74
double f(double x)
double vavilov_accurate_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
Class describing a Vavilov distribution.
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
double G116f1(double x) const
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
virtual double GetBeta2() const
Return the current value of .
static double E1plLog(double x)
double Pdf(double x) const
Evaluate the Vavilov probability density function.
double fLambda[kNquantMax]
double Quantile_c(double z) const
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
virtual ~VavilovAccurate()
Destructor.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
double vavilov_accurate_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function. ...
double vavilov_accurate_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cumulative probability density function.
double vavilov_accurate_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.