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