Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
VavilovFast.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 VavilovFast
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_VavilovFast
32#define ROOT_Math_VavilovFast
33
34
35/**
36 @ingroup StatFunc
37 */
38
39
40#include "Math/Vavilov.h"
41
42namespace ROOT {
43namespace Math {
44
45//____________________________________________________________________________
46/**
47 Class describing a Vavilov distribution.
48
49 The probability density function of the Vavilov distribution
50 as function of Landau's parameter is given by:
51 \f[ p(\lambda_L; \kappa, \beta^2) =
52 \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} \phi(s) e^{\lambda_L s} ds\f]
53 where \f$\phi(s) = e^{C} e^{\psi(s)}\f$
54 with \f$ C = \kappa (1+\beta^2 \gamma )\f$
55 and \f$\psi(s)= s \ln \kappa + (s+\beta^2 \kappa)
56 \cdot \left ( \int \limits_{0}^{1}
57 \frac{1 - e^{\frac{-st}{\kappa}}}{t} \,d t- \gamma \right )
58 - \kappa \, e^{\frac{-s}{\kappa}}\f$.
59 \f$ \gamma = 0.5772156649\dots\f$ is Euler's constant.
60
61 For the class VavilovFast,
62 Pdf returns the Vavilov distribution as function of Landau's parameter
63 \f$\lambda_L = \lambda_V/\kappa - \ln \kappa\f$,
64 which is the convention used in the CERNLIB routines, and in the tables
65 by S.M. Seltzer and M.J. Berger: Energy loss stragglin of protons and mesons:
66 Tabulation of the Vavilov distribution, pp 187-203
67 in: National Research Council (U.S.), Committee on Nuclear Science:
68 Studies in penetration of charged particles in matter,
69 Nat. Akad. Sci. Publication 1133,
70 Nucl. Sci. Series Report No. 39,
71 Washington (Nat. Akad. Sci.) 1964, 388 pp.
72 Available from
73 <A HREF="http://books.google.de/books?id=kmMrAAAAYAAJ&lpg=PP9&pg=PA187#v=onepage&q&f=false">Google books</A>
74
75 Therefore, for small values of \f$\kappa < 0.01\f$,
76 pdf approaches the Landau distribution.
77
78 For values \f$\kappa > 10\f$, the Gauss approximation should be used
79 with \f$\mu\f$ and \f$\sigma\f$ given by Vavilov::mean(kappa, beta2)
80 and sqrt(Vavilov::variance(kappa, beta2).
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 For detailed description see
90 A. Rotondi and P. Montagna, Fast calculation of Vavilov distribution,
91 <A HREF="http://dx.doi.org/10.1016/0168-583X(90)90749-K">Nucl. Instr. and Meth. B47 (1990) 215-224</A>,
92 which has been implemented in
93 <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g115/top.html">
94 CERNLIB (G115)</A>.
95
96 The class stores coefficients needed to calculate \f$p(\lambda; \kappa, \beta^2)\f$
97 for fixed values of \f$\kappa\f$ and \f$\beta^2\f$.
98 Changing these values is computationally expensive.
99
100 The parameter \f$\kappa\f$ must be in the range \f$0.01 \le \kappa \le 12\f$.
101
102 The parameter \f$\beta^2\f$ must be in the range \f$0 \le \beta^2 \le 1\f$.
103
104 Average times on a Pentium Core2 Duo P8400 2.26GHz:
105 - 9.9us per call to SetKappaBeta2 or constructor
106 - 0.095us per call to Pdf, Cdf
107 - 3.7us per first call to Quantile after SetKappaBeta2 or constructor
108 - 0.137us per subsequent call to Quantile
109
110 Benno List, June 2010
111
112 @ingroup StatFunc
113 */
114
115
116class VavilovFast: public Vavilov {
117
118public:
119
120
121 /**
122 Initialize an object to calculate the Vavilov distribution
123
124 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
125 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
126 */
127
128 VavilovFast(double kappa=1, double beta2=1);
129
130
131 /**
132 Destructor
133 */
134 ~VavilovFast() override;
135
136
137public:
138
139 /**
140 Evaluate the Vavilov probability density function
141
142 @param x The Landau parameter \f$x = \lambda_L\f$
143 */
144 double Pdf (double x) const override;
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 must be in the range \f$0.01 \le \kappa \le 12 \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 double Pdf (double x, double kappa, double beta2) override;
155
156 /**
157 Evaluate the Vavilov cumulative probability density function
158
159 @param x The Landau parameter \f$x = \lambda_L\f$
160 */
161 double Cdf (double x) const override;
162
163 /**
164 Evaluate the Vavilov cumulative 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 must be in the range \f$0.01 \le \kappa \le 12 \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 double Cdf (double x, double kappa, double beta2) override;
172
173 /**
174 Evaluate the Vavilov complementary cumulative probability density function
175
176 @param x The Landau parameter \f$x = \lambda_L\f$
177 */
178 double Cdf_c (double x) const override;
179
180 /**
181 Evaluate the Vavilov complementary cumulative 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 must be in the range \f$0.01 \le \kappa \le 12 \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 double Cdf_c (double x, double kappa, double beta2) override;
189
190 /**
191 Evaluate the inverse of the Vavilov cumulative 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 double Quantile (double z) const override;
196
197 /**
198 Evaluate the inverse of the Vavilov cumulative 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 must be in the range \f$0.01 \le \kappa \le 12 \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 double Quantile (double z, double kappa, double beta2) override;
206
207 /**
208 Evaluate the inverse of the complementary Vavilov cumulative 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 double Quantile_c (double z) const override;
213
214 /**
215 Evaluate the inverse of the complementary Vavilov cumulative 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 must be in the range \f$0.01 \le \kappa \le 12 \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 double Quantile_c (double z, double kappa, double beta2) override;
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 must be in the range \f$0.01 \le \kappa \le 12 \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 void SetKappaBeta2 (double kappa, double beta2) override;
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 double GetLambdaMin() const override;
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 double GetLambdaMax() const override;
243
244 /**
245 Return the current value of \f$\kappa\f$
246 */
247 double GetKappa() const override;
248
249 /**
250 Return the current value of \f$\beta^2\f$
251 */
252 double GetBeta2() const override;
253
254 /**
255 Returns a static instance of class VavilovFast
256 */
257 static VavilovFast *GetInstance();
258
259 /**
260 Returns a static instance of class VavilovFast,
261 and sets the values of kappa and beta2
262
263 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \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 static VavilovFast *GetInstance(double kappa, double beta2);
267
268
269private:
270 double fKappa;
271 double fBeta2;
272
273 double fAC[14];
274 double fHC[9];
275 double fWCM[201];
277 int fNpt;
278
280
281};
282
283 /**
284 The Vavilov probability density function
285
286 @param x The Landau parameter \f$x = \lambda_L\f$
287 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
288 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
289
290 @ingroup PdfFunc
291 */
292double vavilov_fast_pdf (double x, double kappa, double beta2);
293
294 /**
295 The Vavilov cumulative probability density function
296
297 @param x The Landau parameter \f$x = \lambda_L\f$
298 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
299 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
300
301 @ingroup ProbFunc
302 */
303double vavilov_fast_cdf (double x, double kappa, double beta2);
304
305 /**
306 The Vavilov complementary cumulative probability density function
307
308 @param x The Landau parameter \f$x = \lambda_L\f$
309 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
310 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
311
312 @ingroup ProbFunc
313 */
314double vavilov_fast_cdf_c (double x, double kappa, double beta2);
315
316 /**
317 The inverse of the Vavilov cumulative probability density function
318
319 @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
320 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \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 @ingroup QuantFunc
324 */
325double vavilov_fast_quantile (double z, double kappa, double beta2);
326
327 /**
328 The inverse of the complementary Vavilov cumulative probability density function
329
330 @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
331 @param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
332 @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
333
334 @ingroup QuantFunc
335 */
336double vavilov_fast_quantile_c (double z, double kappa, double beta2);
337
338} // namespace Math
339} // namespace ROOT
340
341#endif /* ROOT_Math_VavilovFast */
@ Pdf
Definition RooChi2Var.h:1
Class describing a Vavilov distribution.
double GetBeta2() const override
Return the current value of .
~VavilovFast() override
Destructor.
double GetLambdaMax() const override
Return the maximum value of for which is nonzero in the current approximation.
static VavilovFast * fgInstance
double GetKappa() const override
Return the current value of .
void SetKappaBeta2(double kappa, double beta2) override
Change and and recalculate coefficients if necessary.
double Quantile_c(double z) const override
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
double Cdf_c(double x) const override
Evaluate the Vavilov complementary cumulative probability density function.
double Quantile(double z) const override
Evaluate the inverse of the Vavilov cumulative probability density function.
double Cdf(double x) const override
Evaluate the Vavilov cumulative probability density function.
double GetLambdaMin() const override
Return the minimum value of for which is nonzero in the current approximation.
static VavilovFast * GetInstance()
Returns a static instance of class VavilovFast.
Base class describing a Vavilov distribution.
Definition Vavilov.h:120
double vavilov_fast_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
double vavilov_fast_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.
double vavilov_fast_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
double vavilov_fast_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cumulative probability density function.
double vavilov_fast_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function.
Double_t x[n]
Definition legend1.C:17
Namespace for new Math classes and functions.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.