Logo ROOT  
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
42namespace ROOT {
43namespace 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
120class Vavilov {
121
122public:
123
124
125 /**
126 Default constructor
127 */
128 Vavilov();
129
130 /**
131 Destructor
132 */
133 virtual ~Vavilov();
134
135public:
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
virtual double Pdf(double x) const =0
Evaluate the Vavilov probability density function.
virtual ~Vavilov()
Destructor.
Definition: Vavilov.cxx:50
virtual double Cdf(double x) const =0
Evaluate the Vavilov cumulative probability density function.
virtual double Quantile_c(double z, double kappa, double beta2)=0
Evaluate the inverse of the complementary Vavilov cumulative probability density function,...
virtual double Pdf(double x, double kappa, double beta2)=0
Evaluate the Vavilov probability density function, and set kappa and beta2, if necessary.
virtual double Cdf(double x, double kappa, double beta2)=0
Evaluate the Vavilov cumulative probability density function, and set kappa and beta2,...
virtual double GetKappa() const =0
Return the current value of .
virtual double GetLambdaMax() const =0
Return the maximum value of for which is nonzero in the current approximation.
virtual double Quantile(double z) const =0
Evaluate the inverse of the Vavilov cumulative probability density function.
virtual double GetLambdaMin() const =0
Return the minimum value of for which is nonzero in the current approximation.
virtual double Quantile(double z, double kappa, double beta2)=0
Evaluate the inverse of the Vavilov cumulative probability density function, and set kappa and beta2,...
virtual double Cdf_c(double x) const =0
Evaluate the Vavilov complementary cumulative probability density function.
virtual double Variance() const
Return the theoretical variance .
Definition: Vavilov.cxx:88
virtual double GetBeta2() const =0
Return the current value of .
virtual void SetKappaBeta2(double kappa, double beta2)=0
Change and and recalculate coefficients if necessary.
virtual double Quantile_c(double z) const =0
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
virtual double Mode() const
Return the value of where the pdf is maximal.
Definition: Vavilov.cxx:56
virtual double Mean() const
Return the theoretical mean , where is Euler's constant.
Definition: Vavilov.cxx:80
virtual double Skewness() const
Return the theoretical skewness .
Definition: Vavilov.cxx:96
virtual double Kurtosis() const
Return the theoretical kurtosis .
Definition: Vavilov.cxx:105
virtual double Cdf_c(double x, double kappa, double beta2)=0
Evaluate the Vavilov complementary cumulative probability density function, and set kappa and beta2,...
Vavilov()
Default constructor.
Definition: Vavilov.cxx:46
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
VSD Structures.
Definition: StringConv.hxx:21