// @(#)root/mathmore:$Id$
// Authors: B. List 29.4.2010

/**********************************************************************
*                                                                    *
* Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
*                                                                    *
* This library is free software; you can redistribute it and/or      *
* modify it under the terms of the GNU General Public License        *
* as published by the Free Software Foundation; either version 2     *
* of the License, or (at your option) any later version.             *
*                                                                    *
* This library is distributed in the hope that it will be useful,    *
* but WITHOUT ANY WARRANTY; without even the implied warranty of     *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
* General Public License for more details.                           *
*                                                                    *
* You should have received a copy of the GNU General Public License  *
* along with this library (see file COPYING); if not, write          *
* to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
* 330, Boston, MA 02111-1307 USA, or contact the author.             *
*                                                                    *
**********************************************************************/

// Header file for class Vavilov
//
// Created by: blist  at Thu Apr 29 11:19:00 2010
//
// Last update: Thu Apr 29 11:19:00 2010
//
#ifndef ROOT_Math_Vavilov
#define ROOT_Math_Vavilov

/**
@ingroup StatFunc
*/

#include <iostream>

namespace ROOT {
namespace Math {

//____________________________________________________________________________
/**
Base class describing a Vavilov distribution

The Vavilov distribution is defined in
P.V. Vavilov: Ionization losses of high-energy heavy particles,
Sov. Phys. JETP 5 (1957) 749 [Zh. Eksp. Teor. Fiz. 32 (1957) 920].

The probability density function of the Vavilov distribution
as function of Landau's parameter is given by:
\f[ p(\lambda_L; \kappa, \beta^2) =
\frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} \phi(s) e^{\lambda_L s} ds\f]
where \f$\phi(s) = e^{C} e^{\psi(s)}\f$
with  \f$C = \kappa (1+\beta^2 \gamma )\f$
and \f$\psi(s)= s \ln \kappa + (s+\beta^2 \kappa) \cdot \left ( \int \limits_{0}^{1} \frac{1 - e^{\frac{-st}{\kappa}}}{t} \,d t- \gamma \right ) - \kappa \, e^{\frac{-s}{\kappa}}\f$.
\f$\gamma = 0.5772156649\dots\f$ is Euler's constant.

For the class Vavilov,
Pdf returns the Vavilov distribution as function of Landau's parameter
\f$\lambda_L = \lambda_V/\kappa - \ln \kappa\f$,
which is the convention used in the CERNLIB routines, and in the tables
by S.M. Seltzer and M.J. Berger: Energy loss stragglin of protons and mesons:
Tabulation of the Vavilov distribution, pp 187-203
in: National Research Council (U.S.), Committee on Nuclear Science:
Studies in penetration of charged particles in matter,
Nucl. Sci. Series Report No. 39,
Washington (Nat. Akad. Sci.) 1964, 388 pp.
Available from

Therefore, for small values of \f$\kappa < 0.01\f$,
pdf approaches the Landau distribution.

For values \f$\kappa > 10\f$, the Gauss approximation should be used
with \f$\mu\f$ and \f$\sigma\f$ given by Vavilov::Mean(kappa, beta2)
and sqrt(Vavilov::Variance(kappa, beta2).

The original Vavilov pdf is obtained by
v.Pdf(lambdaV/kappa-log(kappa))/kappa.

Two subclasses are provided:
- VavilovFast uses the algorithm by
A. Rotondi and P. Montagna, Fast calculation of Vavilov distribution,
<A HREF="http://dx.doi.org/10.1016/0168-583X(90)90749-K">Nucl. Instr. and Meth. B47 (1990) 215-224</A>,
which has been implemented in
<A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/g115/top.html">
CERNLIB (G115)</A>.

- VavilovAccurate uses the algorithm by
B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
<A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>,
which has been implemented in
<A HREF="http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/g116/top.html">
CERNLIB (G116)</A>.

Both subclasses store coefficients needed to calculate \f$p(\lambda; \kappa, \beta^2)\f$
for fixed values of \f$\kappa\f$ and \f$\beta^2\f$.
Changing these values is computationally expensive.

VavilovFast is about 5 times faster for the calculation of the Pdf than VavilovAccurate;
initialization takes about 100 times longer than calculation of the Pdf value.
For the quantile calculation, VavilovFast
is 30 times faster for the initialization, and 6 times faster for
subsequent calculations. Initialization for Quantile takes
27 (11) times longer than subsequent calls for VavilovFast (VavilovAccurate).

@ingroup StatFunc

*/

class Vavilov {

public:

/**
Default constructor
*/
Vavilov();

/**
Destructor
*/
virtual ~Vavilov();

public:

/**
Evaluate the Vavilov probability density function

@param x The Landau parameter \f$x = \lambda_L\f$

*/
virtual double Pdf (double x) const = 0;

/**
Evaluate the Vavilov probability density function,
and set kappa and beta2, if necessary

@param x The Landau parameter \f$x = \lambda_L\f$
@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
virtual double Pdf (double x, double kappa, double beta2) = 0;

/**
Evaluate the Vavilov cummulative probability density function

@param x The Landau parameter \f$x = \lambda_L\f$
*/
virtual double Cdf (double x) const = 0;

/**
Evaluate the Vavilov cummulative probability density function,
and set kappa and beta2, if necessary

@param x The Landau parameter \f$x = \lambda_L\f$
@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
virtual double Cdf (double x, double kappa, double beta2) = 0;

/**
Evaluate the Vavilov complementary cummulative probability density function

@param x The Landau parameter \f$x = \lambda_L\f$
*/
virtual double Cdf_c (double x) const = 0;

/**
Evaluate the Vavilov complementary cummulative probability density function,
and set kappa and beta2, if necessary

@param x The Landau parameter \f$x = \lambda_L\f$
@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
virtual double Cdf_c (double x, double kappa, double beta2) = 0;

/**
Evaluate the inverse of the Vavilov cummulative probability density function

@param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
*/
virtual double Quantile (double z) const = 0;

/**
Evaluate the inverse of the Vavilov cummulative probability density function,
and set kappa and beta2, if necessary

@param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
virtual double Quantile (double z, double kappa, double beta2) = 0;

/**
Evaluate the inverse of the complementary Vavilov cummulative probability density function

@param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
*/
virtual double Quantile_c (double z) const = 0;

/**
Evaluate the inverse of the complementary Vavilov cummulative probability density function,
and set kappa and beta2, if necessary

@param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
virtual double Quantile_c (double z, double kappa, double beta2) = 0;

/**
Change \f$\kappa\f$ and \f$\beta^2\f$ and recalculate coefficients if necessary

@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
virtual void SetKappaBeta2 (double kappa, double beta2) = 0;

/**
Return the minimum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
is nonzero in the current approximation
*/
virtual double GetLambdaMin() const = 0;

/**
Return the maximum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
is nonzero in the current approximation
*/
virtual double GetLambdaMax() const = 0;

/**
Return the current value of \f$\kappa\f$
*/
virtual double GetKappa()     const = 0;

/**
Return the current value of \f$\beta^2\f$
*/
virtual double GetBeta2()     const = 0;

/**
Return the value of \f$\lambda\f$ where the pdf is maximal
*/
virtual double Mode() const;

/**
Return the value of \f$\lambda\f$ where the pdf is maximal function,
and set kappa and beta2, if necessary

@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
virtual double Mode(double kappa, double beta2);

/**
Return the theoretical mean \f$\mu = \gamma-1- \ln \kappa - \beta^2\f$,
where \f$\gamma = 0.5772\dots\f$ is Euler's constant
*/
virtual double Mean() const;

/**
Return the theoretical variance \f$\sigma^2 = \frac{1 - \beta^2/2}{\kappa}\f$
*/
virtual double Variance() const;

/**
Return the theoretical skewness
\f$\gamma_1 = \frac{1/2 - \beta^2/3}{\kappa^2 \sigma^3} \f$
*/
virtual double Skewness() const;

/**
Return the theoretical kurtosis
\f$\gamma_2 = \frac{1/3 - \beta^2/4}{\kappa^3 \sigma^4}\f$
*/
virtual double Kurtosis() const;

/**
Return the theoretical Mean \f$\mu = \gamma-1- \ln \kappa - \beta^2\f$

@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
static double Mean(double kappa, double beta2);

/**
Return the theoretical Variance \f$\sigma^2 = \frac{1 - \beta^2/2}{\kappa}\f$

@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
static double Variance(double kappa, double beta2);

/**
Return the theoretical skewness
\f$\gamma_1 = \frac{1/2 - \beta^2/3}{\kappa^2 \sigma^3} \f$

@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
static double Skewness(double kappa, double beta2);

/**
Return the theoretical kurtosis
\f$\gamma_2 = \frac{1/3 - \beta^2/4}{\kappa^3 \sigma^4}\f$

@param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
static double Kurtosis(double kappa, double beta2);

};

} // namespace Math
} // namespace ROOT

#endif /* ROOT_Math_Vavilov */

Vavilov.h:1
Vavilov.h:2
Vavilov.h:3
Vavilov.h:4
Vavilov.h:5
Vavilov.h:6
Vavilov.h:7
Vavilov.h:8
Vavilov.h:9
Vavilov.h:10
Vavilov.h:11
Vavilov.h:12
Vavilov.h:13
Vavilov.h:14
Vavilov.h:15
Vavilov.h:16
Vavilov.h:17
Vavilov.h:18
Vavilov.h:19
Vavilov.h:20
Vavilov.h:21
Vavilov.h:22
Vavilov.h:23
Vavilov.h:24
Vavilov.h:25
Vavilov.h:26
Vavilov.h:27
Vavilov.h:28
Vavilov.h:29
Vavilov.h:30
Vavilov.h:31
Vavilov.h:32
Vavilov.h:33
Vavilov.h:34
Vavilov.h:35
Vavilov.h:36
Vavilov.h:37
Vavilov.h:38
Vavilov.h:39
Vavilov.h:40
Vavilov.h:41
Vavilov.h:42
Vavilov.h:43
Vavilov.h:44
Vavilov.h:45
Vavilov.h:46
Vavilov.h:47
Vavilov.h:48
Vavilov.h:49
Vavilov.h:50
Vavilov.h:51
Vavilov.h:52
Vavilov.h:53
Vavilov.h:54
Vavilov.h:55
Vavilov.h:56
Vavilov.h:57
Vavilov.h:58
Vavilov.h:59
Vavilov.h:60
Vavilov.h:61
Vavilov.h:62
Vavilov.h:63
Vavilov.h:64
Vavilov.h:65
Vavilov.h:66
Vavilov.h:67
Vavilov.h:68
Vavilov.h:69
Vavilov.h:70
Vavilov.h:71
Vavilov.h:72
Vavilov.h:73
Vavilov.h:74
Vavilov.h:75
Vavilov.h:76
Vavilov.h:77
Vavilov.h:78
Vavilov.h:79
Vavilov.h:80
Vavilov.h:81
Vavilov.h:82
Vavilov.h:83
Vavilov.h:84
Vavilov.h:85
Vavilov.h:86
Vavilov.h:87
Vavilov.h:88
Vavilov.h:89
Vavilov.h:90
Vavilov.h:91
Vavilov.h:92
Vavilov.h:93
Vavilov.h:94
Vavilov.h:95
Vavilov.h:96
Vavilov.h:97
Vavilov.h:98
Vavilov.h:99
Vavilov.h:100
Vavilov.h:101
Vavilov.h:102
Vavilov.h:103
Vavilov.h:104
Vavilov.h:105
Vavilov.h:106
Vavilov.h:107
Vavilov.h:108
Vavilov.h:109
Vavilov.h:110
Vavilov.h:111
Vavilov.h:112
Vavilov.h:113
Vavilov.h:114
Vavilov.h:115
Vavilov.h:116
Vavilov.h:117
Vavilov.h:118
Vavilov.h:119
Vavilov.h:120
Vavilov.h:121
Vavilov.h:122
Vavilov.h:123
Vavilov.h:124
Vavilov.h:125
Vavilov.h:126
Vavilov.h:127
Vavilov.h:128
Vavilov.h:129
Vavilov.h:130
Vavilov.h:131
Vavilov.h:132
Vavilov.h:133
Vavilov.h:134
Vavilov.h:135
Vavilov.h:136
Vavilov.h:137
Vavilov.h:138
Vavilov.h:139
Vavilov.h:140
Vavilov.h:141
Vavilov.h:142
Vavilov.h:143
Vavilov.h:144
Vavilov.h:145
Vavilov.h:146
Vavilov.h:147
Vavilov.h:148
Vavilov.h:149
Vavilov.h:150
Vavilov.h:151
Vavilov.h:152
Vavilov.h:153
Vavilov.h:154
Vavilov.h:155
Vavilov.h:156
Vavilov.h:157
Vavilov.h:158
Vavilov.h:159
Vavilov.h:160
Vavilov.h:161
Vavilov.h:162
Vavilov.h:163
Vavilov.h:164
Vavilov.h:165
Vavilov.h:166
Vavilov.h:167
Vavilov.h:168
Vavilov.h:169
Vavilov.h:170
Vavilov.h:171
Vavilov.h:172
Vavilov.h:173
Vavilov.h:174
Vavilov.h:175
Vavilov.h:176
Vavilov.h:177
Vavilov.h:178
Vavilov.h:179
Vavilov.h:180
Vavilov.h:181
Vavilov.h:182
Vavilov.h:183
Vavilov.h:184
Vavilov.h:185
Vavilov.h:186
Vavilov.h:187
Vavilov.h:188
Vavilov.h:189
Vavilov.h:190
Vavilov.h:191
Vavilov.h:192
Vavilov.h:193
Vavilov.h:194
Vavilov.h:195
Vavilov.h:196
Vavilov.h:197
Vavilov.h:198
Vavilov.h:199
Vavilov.h:200
Vavilov.h:201
Vavilov.h:202
Vavilov.h:203
Vavilov.h:204
Vavilov.h:205
Vavilov.h:206
Vavilov.h:207
Vavilov.h:208
Vavilov.h:209
Vavilov.h:210
Vavilov.h:211
Vavilov.h:212
Vavilov.h:213
Vavilov.h:214
Vavilov.h:215
Vavilov.h:216
Vavilov.h:217
Vavilov.h:218
Vavilov.h:219
Vavilov.h:220
Vavilov.h:221
Vavilov.h:222
Vavilov.h:223
Vavilov.h:224
Vavilov.h:225
Vavilov.h:226
Vavilov.h:227
Vavilov.h:228
Vavilov.h:229
Vavilov.h:230
Vavilov.h:231
Vavilov.h:232
Vavilov.h:233
Vavilov.h:234
Vavilov.h:235
Vavilov.h:236
Vavilov.h:237
Vavilov.h:238
Vavilov.h:239
Vavilov.h:240
Vavilov.h:241
Vavilov.h:242
Vavilov.h:243
Vavilov.h:244
Vavilov.h:245
Vavilov.h:246
Vavilov.h:247
Vavilov.h:248
Vavilov.h:249
Vavilov.h:250
Vavilov.h:251
Vavilov.h:252
Vavilov.h:253
Vavilov.h:254
Vavilov.h:255
Vavilov.h:256
Vavilov.h:257
Vavilov.h:258
Vavilov.h:259
Vavilov.h:260
Vavilov.h:261
Vavilov.h:262
Vavilov.h:263
Vavilov.h:264
Vavilov.h:265
Vavilov.h:266
Vavilov.h:267
Vavilov.h:268
Vavilov.h:269
Vavilov.h:270
Vavilov.h:271
Vavilov.h:272
Vavilov.h:273
Vavilov.h:274
Vavilov.h:275
Vavilov.h:276
Vavilov.h:277
Vavilov.h:278
Vavilov.h:279
Vavilov.h:280
Vavilov.h:281
Vavilov.h:282
Vavilov.h:283
Vavilov.h:284
Vavilov.h:285
Vavilov.h:286
Vavilov.h:287
Vavilov.h:288
Vavilov.h:289
Vavilov.h:290
Vavilov.h:291
Vavilov.h:292
Vavilov.h:293
Vavilov.h:294
Vavilov.h:295
Vavilov.h:296
Vavilov.h:297
Vavilov.h:298
Vavilov.h:299
Vavilov.h:300
Vavilov.h:301
Vavilov.h:302
Vavilov.h:303
Vavilov.h:304
Vavilov.h:305
Vavilov.h:306
Vavilov.h:307
Vavilov.h:308
Vavilov.h:309
Vavilov.h:310
Vavilov.h:311
Vavilov.h:312
Vavilov.h:313
Vavilov.h:314
Vavilov.h:315
Vavilov.h:316
Vavilov.h:317
Vavilov.h:318
Vavilov.h:319
Vavilov.h:320
Vavilov.h:321
Vavilov.h:322
Vavilov.h:323
Vavilov.h:324
Vavilov.h:325
Vavilov.h:326
Vavilov.h:327
Vavilov.h:328
Vavilov.h:329
Vavilov.h:330
Vavilov.h:331