Logo ROOT   6.10/09
Reference Guide
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 
42 namespace ROOT {
43 namespace 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 
116 class VavilovFast: public Vavilov {
117 
118 public:
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  virtual ~VavilovFast();
135 
136 
137 public:
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;
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);
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;
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);
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;
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);
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;
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);
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;
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);
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  virtual void SetKappaBeta2 (double kappa, double beta2);
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;
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;
243 
244  /**
245  Return the current value of \f$\kappa\f$
246  */
247  virtual double GetKappa() const;
248 
249  /**
250  Return the current value of \f$\beta^2\f$
251  */
252  virtual double GetBeta2() const;
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 
269 private:
270  double fKappa;
271  double fBeta2;
272 
273  double fAC[14];
274  double fHC[9];
275  double fWCM[201];
276  int fItype;
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  */
292 double 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  */
303 double 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  */
314 double 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  */
325 double 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  */
336 double vavilov_fast_quantile_c (double z, double kappa, double beta2);
337 
338 } // namespace Math
339 } // namespace ROOT
340 
341 #endif /* ROOT_Math_VavilovFast */
double vavilov_fast_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
Base class describing a Vavilov distribution.
Definition: Vavilov.h:120
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
virtual ~VavilovFast()
Destructor.
Definition: VavilovFast.cxx:57
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
virtual double GetKappa() const
Return the current value of .
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
Definition: VavilovFast.cxx:62
double Quantile(double z) const
Evaluate the inverse of the Vavilov cumulative probability density function.
double Cdf(double x) const
Evaluate the Vavilov cumulative probability density function.
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
VavilovFast(double kappa=1, double beta2=1)
Initialize an object to calculate the Vavilov distribution.
Definition: VavilovFast.cxx:51
Class describing a Vavilov distribution.
Definition: VavilovFast.h:116
double Cdf_c(double x) const
Evaluate the Vavilov complementary cumulative probability density function.
Double_t x[n]
Definition: legend1.C:17
double vavilov_fast_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function. ...
static VavilovFast * GetInstance()
Returns a static instance of class VavilovFast.
double vavilov_fast_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.
double vavilov_fast_quantile(double z, double kappa, double beta2)
The inverse of the 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 GetBeta2() const
Return the current value of .
static VavilovFast * fgInstance
Definition: VavilovFast.h:279
double Quantile_c(double z) const
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
double Pdf(double x) const
Evaluate the Vavilov probability density function.
double vavilov_fast_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.