// @(#)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 VavilovFast
//
// Created by: blist  at Thu Apr 29 11:19:00 2010
//
// Last update: Thu Apr 29 11:19:00 2010
//
#ifndef ROOT_Math_VavilovFast
#define ROOT_Math_VavilovFast

/**
@ingroup StatFunc
*/

#include "Math/Vavilov.h"

namespace ROOT {
namespace Math {

//____________________________________________________________________________
/**
Class describing a Vavilov distribution.

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 VavilovFast,
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).

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.

For detailed description see
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>.

The class stores 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.

The parameter \f$\kappa\f$ must be in the range \f$0.01 \le \kappa \le 12\f$.

The parameter \f$\beta^2\f$ must be in the range \f$0 \le \beta^2 \le 1\f$.

Average times on a Pentium Core2 Duo P8400 2.26GHz:
- 9.9us per call to SetKappaBeta2 or constructor
- 0.095us per call to Pdf, Cdf
- 3.7us per first call to Quantile after SetKappaBeta2 or constructor
- 0.137us per subsequent call to Quantile

Benno List, June 2010

@ingroup StatFunc
*/

class VavilovFast: public Vavilov {

public:

/**
Initialize an object to calculate the Vavilov distribution

@param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/

VavilovFast(double kappa=1, double beta2=1);

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

public:

/**
Evaluate the Vavilov probability density function

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

/**
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 must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
double Pdf (double x, double kappa, double beta2);

/**
Evaluate the Vavilov cummulative probability density function

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

/**
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 must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
double Cdf (double x, double kappa, double beta2);

/**
Evaluate the Vavilov complementary cummulative probability density function

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

/**
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 must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
double Cdf_c (double x, double kappa, double beta2);

/**
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$
*/
double Quantile (double z) const;

/**
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 must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
double Quantile (double z, double kappa, double beta2);

/**
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$
*/
double Quantile_c (double z) const;

/**
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 must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
*/
double Quantile_c (double z, double kappa, double beta2);

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

@param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \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);

/**
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;

/**
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;

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

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

/**
Returns a static instance of class VavilovFast
*/
static VavilovFast *GetInstance();

/**
Returns a static instance of class VavilovFast,
and sets the values of kappa and beta2

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

private:
double fKappa;
double fBeta2;

double fAC[14];
double fHC[9];
double fWCM[201];
int    fItype;
int    fNpt;

static VavilovFast *fgInstance;

};

/**
The Vavilov probability density function

@param x The Landau parameter \f$x = \lambda_L\f$
@param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

@ingroup PdfFunc
*/
double vavilov_fast_pdf (double x, double kappa, double beta2);

/**
The Vavilov cummulative probability density function

@param x The Landau parameter \f$x = \lambda_L\f$
@param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

@ingroup ProbFunc
*/
double vavilov_fast_cdf (double x, double kappa, double beta2);

/**
The Vavilov complementary cummulative probability density function

@param x The Landau parameter \f$x = \lambda_L\f$
@param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

@ingroup ProbFunc
*/
double vavilov_fast_cdf_c (double x, double kappa, double beta2);

/**
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$
@param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

@ingroup QuantFunc
*/
double vavilov_fast_quantile (double z, double kappa, double beta2);

/**
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$
@param kappa The parameter \f$\kappa\f$, which must be in the range \f$0.01 \le \kappa \le 12 \f$
@param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$

@ingroup QuantFunc
*/
double vavilov_fast_quantile_c (double z, double kappa, double beta2);

} // namespace Math
} // namespace ROOT

#endif /* ROOT_Math_VavilovFast */

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