Logo ROOT  
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
37namespace ROOT {
38namespace 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
132
133public:
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
155public:
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 */
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
331private:
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 */
365double 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 */
376double 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 */
387double 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 */
398double 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 */
409double vavilov_accurate_quantile_c (double z, double kappa, double beta2);
410
411} // namespace Math
412} // namespace ROOT
413
414#endif /* ROOT_Math_VavilovAccurate */
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
Class describing a Vavilov distribution.
double Cdf_c(double x) const
Evaluate the Vavilov complementary cumulative probability density function.
double GetEpsilon() const
Return the current value of .
static VavilovAccurate * GetInstance()
Returns a static instance of class VavilovFast.
double G116f1(double x) const
int Rzero(double a, double b, double &x0, double eps, int mxf, double(VavilovAccurate::*f)(double) const) const
static VavilovAccurate * fgInstance
virtual double GetBeta2() const
Return the current value of .
double fLambda[kNquantMax]
static double E1plLog(double x)
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
double Quantile_c(double z) const
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
virtual double Mode() const
Return the value of where the pdf is maximal.
double Pdf(double x) const
Evaluate the Vavilov probability density function.
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
double GetEpsilonPM() const
Return the current value of .
virtual double GetKappa() const
Return the current value of .
double Quantile(double z) const
Evaluate the inverse of the Vavilov cumulative probability density function.
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
double G116f2(double x) const
void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5)
(Re)Initialize the object
VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5)
Initialize an object to calculate the Vavilov distribution.
double GetNTerms() const
Return the number of terms used in the series expansion.
virtual ~VavilovAccurate()
Destructor.
double Cdf(double x) const
Evaluate the Vavilov cumulative probability density function.
Base class describing a Vavilov distribution.
Definition: Vavilov.h:120
double vavilov_accurate_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
double vavilov_accurate_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
double vavilov_accurate_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.
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_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:96
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:618