Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:10:22

0001 // @(#)root/mathmore:$Id$
0002 // Authors: B. List 29.4.2010
0003 
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
0007   *                                                                    *
0008   * This library is free software; you can redistribute it and/or      *
0009   * modify it under the terms of the GNU General Public License        *
0010   * as published by the Free Software Foundation; either version 2     *
0011   * of the License, or (at your option) any later version.             *
0012   *                                                                    *
0013   * This library is distributed in the hope that it will be useful,    *
0014   * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
0015   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
0016   * General Public License for more details.                           *
0017   *                                                                    *
0018   * You should have received a copy of the GNU General Public License  *
0019   * along with this library (see file COPYING); if not, write          *
0020   * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
0021   * 330, Boston, MA 02111-1307 USA, or contact the author.             *
0022   *                                                                    *
0023   **********************************************************************/
0024 
0025 // Header file for class Vavilov
0026 //
0027 // Created by: blist  at Thu Apr 29 11:19:00 2010
0028 //
0029 // Last update: Thu Apr 29 11:19:00 2010
0030 //
0031 #ifndef ROOT_Math_Vavilov
0032 #define ROOT_Math_Vavilov
0033 
0034 
0035 
0036 /**
0037    @ingroup StatFunc
0038  */
0039 
0040 
0041 
0042 namespace ROOT {
0043 namespace Math {
0044 
0045 //____________________________________________________________________________
0046 /**
0047    Base class describing a Vavilov distribution
0048 
0049    The Vavilov distribution is defined in
0050    P.V. Vavilov: Ionization losses of high-energy heavy particles,
0051    Sov. Phys. JETP 5 (1957) 749 [Zh. Eksp. Teor. Fiz. 32 (1957) 920].
0052 
0053    The probability density function of the Vavilov distribution
0054    as function of Landau's parameter is given by:
0055   \f[ p(\lambda_L; \kappa, \beta^2) =
0056   \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} \phi(s) e^{\lambda_L s} ds\f]
0057    where \f$\phi(s) = e^{C} e^{\psi(s)}\f$
0058    with  \f$ C = \kappa (1+\beta^2 \gamma )\f$
0059    and \f$\psi(s)= s \ln \kappa + (s+\beta^2 \kappa)
0060                \cdot \left ( \int \limits_{0}^{1}
0061                \frac{1 - e^{\frac{-st}{\kappa}}}{t} \,d t- \gamma \right )
0062                - \kappa \, e^{\frac{-s}{\kappa}}\f$.
0063    \f$ \gamma = 0.5772156649\dots\f$ is Euler's constant.
0064 
0065    For the class Vavilov,
0066    Pdf returns the Vavilov distribution as function of Landau's parameter
0067    \f$\lambda_L = \lambda_V/\kappa  - \ln \kappa\f$,
0068    which is the convention used in the CERNLIB routines, and in the tables
0069    by S.M. Seltzer and M.J. Berger: Energy loss stragglin of protons and mesons:
0070    Tabulation of the Vavilov distribution, pp 187-203
0071    in: National Research Council (U.S.), Committee on Nuclear Science:
0072    Studies in penetration of charged particles in matter,
0073    Nat. Akad. Sci. Publication 1133,
0074    Nucl. Sci. Series Report No. 39,
0075    Washington (Nat. Akad. Sci.) 1964, 388 pp.
0076    Available from
0077    <A HREF="http://books.google.de/books?id=kmMrAAAAYAAJ&lpg=PP9&pg=PA187#v=onepage&q&f=false">Google books</A>
0078 
0079    Therefore, for small values of \f$\kappa < 0.01\f$,
0080    pdf approaches the Landau distribution.
0081 
0082    For values \f$\kappa > 10\f$, the Gauss approximation should be used
0083    with \f$\mu\f$ and \f$\sigma\f$ given by Vavilov::Mean(kappa, beta2)
0084    and sqrt(Vavilov::Variance(kappa, beta2).
0085 
0086    The original Vavilov pdf is obtained by
0087    v.Pdf(lambdaV/kappa-log(kappa))/kappa.
0088 
0089    Two subclasses are provided:
0090    - VavilovFast uses the algorithm by
0091    A. Rotondi and P. Montagna, Fast calculation of Vavilov distribution,
0092    <A HREF="http://dx.doi.org/10.1016/0168-583X(90)90749-K">Nucl. Instr. and Meth. B47 (1990) 215-224</A>,
0093    which has been implemented in
0094    <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g115/top.html">
0095    CERNLIB (G115)</A>.
0096 
0097    - VavilovAccurate uses the algorithm by
0098    B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
0099    <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>,
0100    which has been implemented in
0101    <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g116/top.html">
0102    CERNLIB (G116)</A>.
0103 
0104    Both subclasses store coefficients needed to calculate \f$p(\lambda; \kappa, \beta^2)\f$
0105    for fixed values of \f$\kappa\f$ and \f$\beta^2\f$.
0106    Changing these values is computationally expensive.
0107 
0108    VavilovFast is about 5 times faster for the calculation of the Pdf than VavilovAccurate;
0109    initialization takes about 100 times longer than calculation of the Pdf value.
0110    For the quantile calculation, VavilovFast
0111    is 30 times faster for the initialization, and 6 times faster for
0112    subsequent calculations. Initialization for Quantile takes
0113    27 (11) times longer than subsequent calls for VavilovFast (VavilovAccurate).
0114 
0115    @ingroup StatFunc
0116 
0117    */
0118 
0119 
0120 class Vavilov {
0121 
0122 public:
0123 
0124 
0125    /**
0126      Default constructor
0127    */
0128    Vavilov();
0129 
0130    /**
0131      Destructor
0132    */
0133    virtual ~Vavilov();
0134 
0135 public:
0136 
0137    /**
0138        Evaluate the Vavilov probability density function
0139 
0140        @param x The Landau parameter \f$x = \lambda_L\f$
0141 
0142    */
0143    virtual double Pdf (double x) const = 0;
0144 
0145    /**
0146        Evaluate the Vavilov probability density function,
0147        and set kappa and beta2, if necessary
0148 
0149        @param x The Landau parameter \f$x = \lambda_L\f$
0150        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0151        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0152    */
0153    virtual double Pdf (double x, double kappa, double beta2) = 0;
0154 
0155    /**
0156        Evaluate the Vavilov cumulative probability density function
0157 
0158        @param x The Landau parameter \f$x = \lambda_L\f$
0159    */
0160    virtual double Cdf (double x) const = 0;
0161 
0162    /**
0163        Evaluate the Vavilov cumulative probability density function,
0164        and set kappa and beta2, if necessary
0165 
0166        @param x The Landau parameter \f$x = \lambda_L\f$
0167        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0168        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0169    */
0170    virtual double Cdf (double x, double kappa, double beta2) = 0;
0171 
0172    /**
0173        Evaluate the Vavilov complementary cumulative probability density function
0174 
0175        @param x The Landau parameter \f$x = \lambda_L\f$
0176    */
0177    virtual double Cdf_c (double x) const = 0;
0178 
0179    /**
0180        Evaluate the Vavilov complementary cumulative probability density function,
0181        and set kappa and beta2, if necessary
0182 
0183        @param x The Landau parameter \f$x = \lambda_L\f$
0184        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0185        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0186    */
0187    virtual double Cdf_c (double x, double kappa, double beta2) = 0;
0188 
0189    /**
0190        Evaluate the inverse of the Vavilov cumulative probability density function
0191 
0192        @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
0193    */
0194    virtual double Quantile (double z) const = 0;
0195 
0196    /**
0197        Evaluate the inverse of the Vavilov cumulative probability density function,
0198        and set kappa and beta2, if necessary
0199 
0200        @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
0201        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0202        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0203    */
0204    virtual double Quantile (double z, double kappa, double beta2) = 0;
0205 
0206    /**
0207        Evaluate the inverse of the complementary Vavilov cumulative probability density function
0208 
0209        @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
0210    */
0211    virtual double Quantile_c (double z) const = 0;
0212 
0213    /**
0214        Evaluate the inverse of the complementary Vavilov cumulative probability density function,
0215        and set kappa and beta2, if necessary
0216 
0217        @param z The argument \f$z\f$, which must be in the range \f$0 \le z \le 1\f$
0218        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0219        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0220    */
0221    virtual double Quantile_c (double z, double kappa, double beta2) = 0;
0222 
0223    /**
0224       Change \f$\kappa\f$ and \f$\beta^2\f$ and recalculate coefficients if necessary
0225 
0226        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0227        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0228    */
0229    virtual void SetKappaBeta2 (double kappa, double beta2) = 0;
0230 
0231    /**
0232       Return the minimum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
0233       is nonzero in the current approximation
0234    */
0235    virtual double GetLambdaMin() const = 0;
0236 
0237    /**
0238       Return the maximum value of \f$\lambda\f$ for which \f$p(\lambda; \kappa, \beta^2)\f$
0239       is nonzero in the current approximation
0240    */
0241    virtual double GetLambdaMax() const = 0;
0242 
0243    /**
0244       Return the current value of \f$\kappa\f$
0245    */
0246    virtual double GetKappa()     const = 0;
0247 
0248    /**
0249       Return the current value of \f$\beta^2\f$
0250    */
0251    virtual double GetBeta2()     const = 0;
0252 
0253    /**
0254       Return the value of \f$\lambda\f$ where the pdf is maximal
0255    */
0256    virtual double Mode() const;
0257 
0258    /**
0259       Return the value of \f$\lambda\f$ where the pdf is maximal function,
0260        and set kappa and beta2, if necessary
0261 
0262        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0263        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0264    */
0265    virtual double Mode(double kappa, double beta2);
0266 
0267    /**
0268       Return the theoretical mean \f$\mu = \gamma-1- \ln \kappa - \beta^2\f$,
0269       where \f$\gamma = 0.5772\dots\f$ is Euler's constant
0270    */
0271    virtual double Mean() const;
0272 
0273    /**
0274       Return the theoretical variance \f$\sigma^2 = \frac{1 - \beta^2/2}{\kappa}\f$
0275    */
0276    virtual double Variance() const;
0277 
0278    /**
0279       Return the theoretical skewness
0280       \f$\gamma_1 = \frac{1/2 - \beta^2/3}{\kappa^2 \sigma^3} \f$
0281    */
0282    virtual double Skewness() const;
0283 
0284    /**
0285       Return the theoretical kurtosis
0286       \f$\gamma_2 = \frac{1/3 - \beta^2/4}{\kappa^3 \sigma^4}\f$
0287    */
0288    virtual double Kurtosis() const;
0289 
0290    /**
0291       Return the theoretical Mean \f$\mu = \gamma-1- \ln \kappa - \beta^2\f$
0292 
0293        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0294        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0295    */
0296    static double Mean(double kappa, double beta2);
0297 
0298    /**
0299       Return the theoretical Variance \f$\sigma^2 = \frac{1 - \beta^2/2}{\kappa}\f$
0300 
0301        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0302        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0303    */
0304    static double Variance(double kappa, double beta2);
0305 
0306    /**
0307       Return the theoretical skewness
0308       \f$\gamma_1 = \frac{1/2 - \beta^2/3}{\kappa^2 \sigma^3} \f$
0309 
0310        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0311        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0312    */
0313    static double Skewness(double kappa, double beta2);
0314 
0315    /**
0316       Return the theoretical kurtosis
0317       \f$\gamma_2 = \frac{1/3 - \beta^2/4}{\kappa^3 \sigma^4}\f$
0318 
0319        @param kappa The parameter \f$\kappa\f$, which should be in the range \f$0.01 \le \kappa \le 10 \f$
0320        @param beta2 The parameter \f$\beta^2\f$, which must be in the range \f$0 \le \beta^2 \le 1 \f$
0321    */
0322    static double Kurtosis(double kappa, double beta2);
0323 
0324 
0325 };
0326 
0327 } // namespace Math
0328 } // namespace ROOT
0329 
0330 #endif /* ROOT_Math_Vavilov */