Logo ROOT   6.10/09
Reference Guide
Vavilov.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 Vavilov
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_Vavilov
32 #define ROOT_Math_Vavilov
33 
34 
35 
36 /**
37  @ingroup StatFunc
38  */
39 
40 
41 
42 namespace ROOT {
43 namespace Math {
44 
45 //____________________________________________________________________________
46 /**
47  Base class describing a Vavilov distribution
48 
49  The Vavilov distribution is defined in
50  P.V. Vavilov: Ionization losses of high-energy heavy particles,
51  Sov. Phys. JETP 5 (1957) 749 [Zh. Eksp. Teor. Fiz. 32 (1957) 920].
52 
53  The probability density function of the Vavilov distribution
54  as function of Landau's parameter is given by:
55  \f[ p(\lambda_L; \kappa, \beta^2) =
56  \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} \phi(s) e^{\lambda_L s} ds\f]
57  where \f$\phi(s) = e^{C} e^{\psi(s)}\f$
58  with \f$ C = \kappa (1+\beta^2 \gamma )\f$
59  and \f$\psi(s)= s \ln \kappa + (s+\beta^2 \kappa)
60  \cdot \left ( \int \limits_{0}^{1}
61  \frac{1 - e^{\frac{-st}{\kappa}}}{t} \,d t- \gamma \right )
62  - \kappa \, e^{\frac{-s}{\kappa}}\f$.
63  \f$ \gamma = 0.5772156649\dots\f$ is Euler's constant.
64 
65  For the class Vavilov,
66  Pdf returns the Vavilov distribution as function of Landau's parameter
67  \f$\lambda_L = \lambda_V/\kappa - \ln \kappa\f$,
68  which is the convention used in the CERNLIB routines, and in the tables
69  by S.M. Seltzer and M.J. Berger: Energy loss stragglin of protons and mesons:
70  Tabulation of the Vavilov distribution, pp 187-203
71  in: National Research Council (U.S.), Committee on Nuclear Science:
72  Studies in penetration of charged particles in matter,
73  Nat. Akad. Sci. Publication 1133,
74  Nucl. Sci. Series Report No. 39,
75  Washington (Nat. Akad. Sci.) 1964, 388 pp.
76  Available from
77  <A HREF="http://books.google.de/books?id=kmMrAAAAYAAJ&lpg=PP9&pg=PA187#v=onepage&q&f=false">Google books</A>
78 
79  Therefore, for small values of \f$\kappa < 0.01\f$,
80  pdf approaches the Landau distribution.
81 
82  For values \f$\kappa > 10\f$, the Gauss approximation should be used
83  with \f$\mu\f$ and \f$\sigma\f$ given by Vavilov::Mean(kappa, beta2)
84  and sqrt(Vavilov::Variance(kappa, beta2).
85 
86  The original Vavilov pdf is obtained by
87  v.Pdf(lambdaV/kappa-log(kappa))/kappa.
88 
89  Two subclasses are provided:
90  - VavilovFast uses the algorithm by
91  A. Rotondi and P. Montagna, Fast calculation of Vavilov distribution,
92  <A HREF="http://dx.doi.org/10.1016/0168-583X(90)90749-K">Nucl. Instr. and Meth. B47 (1990) 215-224</A>,
93  which has been implemented in
94  <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g115/top.html">
95  CERNLIB (G115)</A>.
96 
97  - VavilovAccurate uses the algorithm by
98  B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
99  <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>,
100  which has been implemented in
101  <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g116/top.html">
102  CERNLIB (G116)</A>.
103 
104  Both subclasses store coefficients needed to calculate \f$p(\lambda; \kappa, \beta^2)\f$
105  for fixed values of \f$\kappa\f$ and \f$\beta^2\f$.
106  Changing these values is computationally expensive.
107 
108  VavilovFast is about 5 times faster for the calculation of the Pdf than VavilovAccurate;
109  initialization takes about 100 times longer than calculation of the Pdf value.
110  For the quantile calculation, VavilovFast
111  is 30 times faster for the initialization, and 6 times faster for
112  subsequent calculations. Initialization for Quantile takes
113  27 (11) times longer than subsequent calls for VavilovFast (VavilovAccurate).
114 
115  @ingroup StatFunc
116 
117  */
118 
119 
120 class Vavilov {
121 
122 public:
123 
124 
125  /**
126  Default constructor
127  */
128  Vavilov();
129 
130  /**
131  Destructor
132  */
133  virtual ~Vavilov();
134 
135 public:
136 
137  /**
138  Evaluate the Vavilov probability density function
139 
140  @param x The Landau parameter \f$x = \lambda_L\f$
141 
142  */
143  virtual double Pdf (double x) const = 0;
144 
145  /**
146  Evaluate the Vavilov probability density function,
147  and set kappa and beta2, if necessary
148 
149  @param x The Landau parameter \f$x = \lambda_L\f$
150  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
151  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
152  */
153  virtual double Pdf (double x, double kappa, double beta2) = 0;
154 
155  /**
156  Evaluate the Vavilov cumulative probability density function
157 
158  @param x The Landau parameter \f$x = \lambda_L\f$
159  */
160  virtual double Cdf (double x) const = 0;
161 
162  /**
163  Evaluate the Vavilov cumulative probability density function,
164  and set kappa and beta2, if necessary
165 
166  @param x The Landau parameter \f$x = \lambda_L\f$
167  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
168  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
169  */
170  virtual double Cdf (double x, double kappa, double beta2) = 0;
171 
172  /**
173  Evaluate the Vavilov complementary cumulative probability density function
174 
175  @param x The Landau parameter \f$x = \lambda_L\f$
176  */
177  virtual double Cdf_c (double x) const = 0;
178 
179  /**
180  Evaluate the Vavilov complementary cumulative probability density function,
181  and set kappa and beta2, if necessary
182 
183  @param x The Landau parameter \f$x = \lambda_L\f$
184  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
185  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
186  */
187  virtual double Cdf_c (double x, double kappa, double beta2) = 0;
188 
189  /**
190  Evaluate the inverse of the Vavilov cumulative probability density function
191 
192  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
193  */
194  virtual double Quantile (double z) const = 0;
195 
196  /**
197  Evaluate the inverse of the Vavilov cumulative probability density function,
198  and set kappa and beta2, if necessary
199 
200  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
201  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
202  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
203  */
204  virtual double Quantile (double z, double kappa, double beta2) = 0;
205 
206  /**
207  Evaluate the inverse of the complementary Vavilov cumulative probability density function
208 
209  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
210  */
211  virtual double Quantile_c (double z) const = 0;
212 
213  /**
214  Evaluate the inverse of the complementary Vavilov cumulative probability density function,
215  and set kappa and beta2, if necessary
216 
217  @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
218  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
219  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
220  */
221  virtual double Quantile_c (double z, double kappa, double beta2) = 0;
222 
223  /**
224  Change \f$\kappa\f$ and \f$\beta^2\f$ and recalculate coefficients if necessary
225 
226  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
227  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
228  */
229  virtual void SetKappaBeta2 (double kappa, double beta2) = 0;
230 
231  /**
232  Return the minimum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
233  is nonzero in the current approximation
234  */
235  virtual double GetLambdaMin() const = 0;
236 
237  /**
238  Return the maximum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
239  is nonzero in the current approximation
240  */
241  virtual double GetLambdaMax() const = 0;
242 
243  /**
244  Return the current value of \f$\kappa\f$
245  */
246  virtual double GetKappa() const = 0;
247 
248  /**
249  Return the current value of \f$\beta^2\f$
250  */
251  virtual double GetBeta2() const = 0;
252 
253  /**
254  Return the value of \f$\lambda\f$ where the pdf is maximal
255  */
256  virtual double Mode() const;
257 
258  /**
259  Return the value of \f$\lambda\f$ where the pdf is maximal function,
260  and set kappa and beta2, if necessary
261 
262  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
263  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
264  */
265  virtual double Mode(double kappa, double beta2);
266 
267  /**
268  Return the theoretical mean \f$\mu = \gamma-1- \ln \kappa - \beta^2\f$,
269  where \f$\gamma = 0.5772\dots\f$ is Euler's constant
270  */
271  virtual double Mean() const;
272 
273  /**
274  Return the theoretical variance \f$\sigma^2 = \frac{1 - \beta^2/2}{\kappa}\f$
275  */
276  virtual double Variance() const;
277 
278  /**
279  Return the theoretical skewness
280  \f$\gamma_1 = \frac{1/2 - \beta^2/3}{\kappa^2 \sigma^3} \f$
281  */
282  virtual double Skewness() const;
283 
284  /**
285  Return the theoretical kurtosis
286  \f$\gamma_2 = \frac{1/3 - \beta^2/4}{\kappa^3 \sigma^4}\f$
287  */
288  virtual double Kurtosis() const;
289 
290  /**
291  Return the theoretical Mean \f$\mu = \gamma-1- \ln \kappa - \beta^2\f$
292 
293  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
294  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
295  */
296  static double Mean(double kappa, double beta2);
297 
298  /**
299  Return the theoretical Variance \f$\sigma^2 = \frac{1 - \beta^2/2}{\kappa}\f$
300 
301  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
302  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
303  */
304  static double Variance(double kappa, double beta2);
305 
306  /**
307  Return the theoretical skewness
308  \f$\gamma_1 = \frac{1/2 - \beta^2/3}{\kappa^2 \sigma^3} \f$
309 
310  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
311  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
312  */
313  static double Skewness(double kappa, double beta2);
314 
315  /**
316  Return the theoretical kurtosis
317  \f$\gamma_2 = \frac{1/3 - \beta^2/4}{\kappa^3 \sigma^4}\f$
318 
319  @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
320  @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
321  */
322  static double Kurtosis(double kappa, double beta2);
323 
324 
325 };
326 
327 } // namespace Math
328 } // namespace ROOT
329 
330 #endif /* ROOT_Math_Vavilov */
Base class describing a Vavilov distribution.
Definition: Vavilov.h:120
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
virtual double Variance() const
Return the theoretical variance .
Definition: Vavilov.cxx:88
virtual double Cdf_c(double x) const =0
Evaluate the Vavilov complementary cumulative probability density function.
virtual double Pdf(double x) const =0
Evaluate the Vavilov probability density function.
virtual double GetLambdaMin() const =0
Return the minimum value of for which is nonzero in the current approximation.
Double_t x[n]
Definition: legend1.C:17
Vavilov()
Default constructor.
Definition: Vavilov.cxx:46
virtual double Kurtosis() const
Return the theoretical kurtosis .
Definition: Vavilov.cxx:105
virtual ~Vavilov()
Destructor.
Definition: Vavilov.cxx:50
virtual double Cdf(double x) const =0
Evaluate the Vavilov cumulative probability density function.
virtual double Quantile(double z) const =0
Evaluate the inverse of the Vavilov cumulative probability density function.
virtual double Skewness() const
Return the theoretical skewness .
Definition: Vavilov.cxx:96
virtual double Mode() const
Return the value of where the pdf is maximal.
Definition: Vavilov.cxx:56
virtual double GetBeta2() const =0
Return the current value of .
virtual double Quantile_c(double z) const =0
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
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 GetLambdaMax() const =0
Return the maximum value of for which is nonzero in the current approximation.
virtual void SetKappaBeta2(double kappa, double beta2)=0
Change and and recalculate coefficients if necessary.
virtual double GetKappa() const =0
Return the current value of .
virtual double Mean() const
Return the theoretical mean , where is Euler&#39;s constant.
Definition: Vavilov.cxx:80