# TMath

Encapsulate most frequently used Math functions.
NB. The basic functions Min, Max, Abs and Sign are defined in TMathBase.

namespace TMath

## namespace

   /* ************************* */
/* * Fundamental constants * */
/* ************************* */

inline Double_t Pi()       { return 3.14159265358979323846; }
inline Double_t TwoPi()    { return 2.0 * Pi(); }
inline Double_t PiOver2()  { return Pi() / 2.0; }
inline Double_t PiOver4()  { return Pi() / 4.0; }
inline Double_t InvPi()    { return 1.0 / Pi(); }
inline Double_t RadToDeg() { return 180.0 / Pi(); }
inline Double_t DegToRad() { return Pi() / 180.0; }
inline Double_t Sqrt2()    { return 1.4142135623730950488016887242097; }

// e (base of natural log)
inline Double_t E()        { return 2.71828182845904523536; }

// natural log of 10 (to convert log to ln)
inline Double_t Ln10()     { return 2.30258509299404568402; }

// base-10 log of e  (to convert ln to log)
inline Double_t LogE()     { return 0.43429448190325182765; }

// velocity of light
inline Double_t C()        { return 2.99792458e8; }        // m s^-1
inline Double_t Ccgs()     { return 100.0 * C(); }         // cm s^-1
inline Double_t CUncertainty() { return 0.0; }             // exact

// gravitational constant
inline Double_t G()        { return 6.673e-11; }           // m^3 kg^-1 s^-2
inline Double_t Gcgs()     { return G() / 1000.0; }        // cm^3 g^-1 s^-2
inline Double_t GUncertainty() { return 0.010e-11; }

// G over h-bar C
inline Double_t GhbarC()   { return 6.707e-39; }           // (GeV/c^2)^-2
inline Double_t GhbarCUncertainty() { return 0.010e-39; }

// standard acceleration of gravity
inline Double_t Gn()       { return 9.80665; }             // m s^-2
inline Double_t GnUncertainty() { return 0.0; }            // exact

// Planck's constant
inline Double_t H()        { return 6.62606876e-34; }      // J s
inline Double_t Hcgs()     { return 1.0e7 * H(); }         // erg s
inline Double_t HUncertainty() { return 0.00000052e-34; }

// h-bar (h over 2 pi)
inline Double_t Hbar()     { return 1.054571596e-34; }     // J s
inline Double_t Hbarcgs()  { return 1.0e7 * Hbar(); }      // erg s
inline Double_t HbarUncertainty() { return 0.000000082e-34; }

// hc (h * c)
inline Double_t HC()       { return H() * C(); }           // J m
inline Double_t HCcgs()    { return Hcgs() * Ccgs(); }     // erg cm

// Boltzmann's constant
inline Double_t K()        { return 1.3806503e-23; }       // J K^-1
inline Double_t Kcgs()     { return 1.0e7 * K(); }         // erg K^-1
inline Double_t KUncertainty() { return 0.0000024e-23; }

// Stefan-Boltzmann constant
inline Double_t Sigma()    { return 5.6704e-8; }           // W m^-2 K^-4
inline Double_t SigmaUncertainty() { return 0.000040e-8; }

inline Double_t Na()       { return 6.02214199e+23; }      // mol^-1
inline Double_t NaUncertainty() { return 0.00000047e+23; }

// universal gas constant (Na * K)
// http://scienceworld.wolfram.com/physics/UniversalGasConstant.html
inline Double_t R()        { return K() * Na(); }          // J K^-1 mol^-1
inline Double_t RUncertainty() { return R()*((KUncertainty()/K()) + (NaUncertainty()/Na())); }

// Molecular weight of dry air
// 1976 US Standard Atmosphere,
// also see http://atmos.nmsu.edu/jsdap/encyclopediawork.html
inline Double_t MWair()    { return 28.9644; }             // kg kmol^-1 (or gm mol^-1)

// Dry Air Gas Constant (R / MWair)
// http://atmos.nmsu.edu/education_and_outreach/encyclopedia/gas_constant.htm
inline Double_t Rgair()    { return (1000.0 * R()) / MWair(); }  // J kg^-1 K^-1

// Euler-Mascheroni Constant
inline Double_t EulerGamma() { return 0.577215664901532860606512090082402431042; }

// Elementary charge
inline Double_t Qe()       { return 1.602176462e-19; }     // C
inline Double_t QeUncertainty() { return 0.000000063e-19; }

/* ************************** */
/* * Mathematical Functions * */
/* ************************** */

/* ***************************** */
/* * Trigonometrical Functions * */
/* ***************************** */
inline Double_t Sin(Double_t);
inline Double_t Cos(Double_t);
inline Double_t Tan(Double_t);
inline Double_t SinH(Double_t);
inline Double_t CosH(Double_t);
inline Double_t TanH(Double_t);
inline Double_t ASin(Double_t);
inline Double_t ACos(Double_t);
inline Double_t ATan(Double_t);
inline Double_t ATan2(Double_t, Double_t);
Double_t ASinH(Double_t);
Double_t ACosH(Double_t);
Double_t ATanH(Double_t);
Double_t Hypot(Double_t x, Double_t y);

/* ************************ */
/* * Elementary Functions * */
/* ************************ */
inline Double_t Ceil(Double_t x);
inline Int_t    CeilNint(Double_t x);
inline Double_t Floor(Double_t x);
inline Int_t    FloorNint(Double_t x);
template<typename T>
inline Int_t    Nint(T x);

inline Double_t Sq(Double_t x);
inline Double_t Sqrt(Double_t x);
inline Double_t Exp(Double_t x);
inline Double_t Ldexp(Double_t x, Int_t exp);
Double_t Factorial(Int_t i);/// Compute factorial(n).
inline LongDouble_t Power(LongDouble_t x, LongDouble_t y);
inline LongDouble_t Power(LongDouble_t x, Long64_t y);
inline LongDouble_t Power(Long64_t x, Long64_t y);
inline Double_t Power(Double_t x, Double_t y);
inline Double_t Power(Double_t x, Int_t y);
inline Double_t Log(Double_t x);
Double_t Log2(Double_t x);
inline Double_t Log10(Double_t x);
inline Int_t    Finite(Double_t x);
inline Int_t    IsNaN(Double_t x);

inline Double_t QuietNaN();
inline Double_t SignalingNaN();
inline Double_t Infinity();

template <typename T>
struct Limits {
inline static T Min();
inline static T Max();
inline static T Epsilon();
};

// Some integer math
Long_t   Hypot(Long_t x, Long_t y);     // sqrt(px*px + py*py)

// Comparing floating points
inline Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon) {
//return kTRUE if absolute difference between af and bf is less than epsilon
return TMath::Abs(af-bf) < epsilon;
}
inline Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec) {
//return kTRUE if relative difference between af and bf is less than relPrec
return TMath::Abs(af-bf) <= 0.5*relPrec*(TMath::Abs(af)+TMath::Abs(bf));
}

/* ******************** */
/* * Array Algorithms * */
/* ******************** */

// Min, Max of an array
template <typename T> T MinElement(Long64_t n, const T *a);
template <typename T> T MaxElement(Long64_t n, const T *a);

// Locate Min, Max element number in an array
template <typename T> Long64_t  LocMin(Long64_t n, const T *a);
template <typename Iterator> Iterator LocMin(Iterator first, Iterator last);
template <typename T> Long64_t  LocMax(Long64_t n, const T *a);
template <typename Iterator> Iterator LocMax(Iterator first, Iterator last);

// Binary search
template <typename T> Long64_t BinarySearch(Long64_t n, const T  *array, T value);
template <typename T> Long64_t BinarySearch(Long64_t n, const T **array, T value);
template <typename Iterator, typename Element> Iterator BinarySearch(Iterator first, Iterator last, Element value);

// Hashing
ULong_t Hash(const void *txt, Int_t ntxt);
/// Calculates hash index from any char string.
/// Based on precalculated table of 256 specially selected numbers.
/// These numbers are selected in such a way, that for string
/// length == 4 (integer number) the hash is unambigous, i.e.
/// from hash value we can recalculate input (no degeneration).
/// The quality of hash method is good enough, that
/// "random" numbers made as R = Hash(1), Hash(2), ...Hash(N)
/// tested by <R>, <R*R>, <Ri*Ri+1> gives the same result
/// as for libc rand().
/// For string:  i = TMath::Hash(string,nstring);
/// For int:     i = TMath::Hash(&intword,sizeof(int));
/// For pointer: i = TMath::Hash(&pointer,sizeof(void*));
///              V.Perev
/// This function is kept for back compatibility. The code previously in this function
/// has been moved to the static function TString::Hash

ULong_t Hash(const char *str);

// Sorting
template <typename Element, typename Index>
void Sort(Index n, const Element* a, Index* index, Bool_t down=kTRUE);
template <typename Iterator, typename IndexIterator>
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE);

void BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2);
/// Bubble sort variant to obtain the order of an array's elements into
/// an index in order to do more useful things than the standard built
/// in functions.
/// *arr1 is unchanged;
/// *arr2 is the array of indicies corresponding to the decending value
/// of arr1 with arr2 corresponding to the largest arr1 value and
/// arr2[Narr] the smallest.

void BubbleLow (Int_t Narr, Double_t *arr1, Int_t *arr2);
/// Opposite ordering of the array arr2[] to that of BubbleHigh.

Bool_t   Permute(Int_t n, Int_t *a); // Find permutations
/// Simple recursive algorithm to find the permutations of
/// n natural numbers, not necessarily all distinct
/// adapted from CERNLIB routine PERMU.
/// The input array has to be initialised with a non descending
/// sequence. The method returns kFALSE when
/// all combinations are exhausted.

/* ************************* */
/* * Geometrical Functions * */
/* ************************* */

//Sample quantiles
void      Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob,
Bool_t isSorted=kTRUE, Int_t *index = 0, Int_t type=7);
///Computes sample quantiles, corresponding to the given probabilities
///Parameters:
///  x -the data sample
///  n - its size
///  quantiles - computed quantiles are returned in there
///  prob - probabilities where to compute quantiles
///  nprob - size of prob array
///  isSorted - is the input array x sorted?
///  NOTE, that when the input is not sorted, an array of integers of size n needs
///        to be allocated. It can be passed by the user in parameter index,
///        or, if not passed, it will be allocated inside the function
///  type - method to compute (from 1 to 9). Following types are provided:
///  Discontinuous:
///    type=1 - inverse of the empirical distribution function
///    type=2 - like type 1, but with averaging at discontinuities
///    type=3 - SAS definition: nearest even order statistic
///  Piecwise linear continuous:
///    In this case, sample quantiles can be obtained by linear interpolation
///    between the k-th order statistic and p(k).
///    type=4 - linear interpolation of empirical cdf, p(k)=k/n;
///    type=5 - a very popular definition, p(k) = (k-0.5)/n;
///    type=6 - used by Minitab and SPSS, p(k) = k/(n+1);
///    type=7 - used by S-Plus and R, p(k) = (k-1)/(n-1);
///    type=8 - resulting sample quantiles are approximately median unbiased
///             regardless of the distribution of x. p(k) = (k-1/3)/(n+1/3);
///    type=9 - resulting sample quantiles are approximately unbiased, when
///             the sample comes from Normal distribution. p(k)=(k-3/8)/(n+1/4);
///    default type = 7
/// References:
/// 1) Hyndman, R.J and Fan, Y, (1996) "Sample quantiles in statistical packages"
///                                     American Statistician, 50, 361-365
/// 2) R Project documentation for the function quantile of package {stats}

// IsInside
template <typename T> Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y);

// Calculate the Cross Product of two vectors
template <typename T> T *Cross(const T v1,const T v2, T out);

Float_t   Normalize(Float_t v);  // Normalize a vector
/// Normalize a vector v in place.
/// Returns the norm of the original vector.

Double_t  Normalize(Double_t v); // Normalize a vector
/// Normalize a vector v in place.
/// Returns the norm of the original vector.
/// This implementation (thanks Kevin Lynch <krlynch@bu.edu>) is protected
/// against possible overflows.

//Calculate the Normalized Cross Product of two vectors
template <typename T> inline T NormCross(const T v1,const T v2,T out);

// Calculate a normal vector of a plane
template <typename T> T *Normal2Plane(const T v1,const T v2,const T v3, T normal);

/* ************************ */
/* * Polynomial Functions * */
/* ************************ */

Bool_t    RootsCubic(const Double_t coef,Double_t &a, Double_t &b, Double_t &c);
/// Calculates roots of polynomial of 3rd order a*x^3 + b*x^2 + c*x + d, where
/// a == coef, b == coef, c == coef, d == coef
///coef must be different from 0
/// If the boolean returned by the method is false:
///    ==> there are 3 real roots a,b,c
/// If the boolean returned by the method is true:
///    ==> there is one real root a and 2 complex conjugates roots (b+i*c,b-i*c)
/// Author: Francois-Xavier Gentit

/* *********************** */
/* * Statistic Functions * */
/* *********************** */

Double_t Binomial(Int_t n,Int_t k);  // Calculate the binomial coefficient n over k
Double_t BinomialI(Double_t p, Int_t n, Int_t k);
/// Suppose an event occurs with probability _p_ per trial
/// Then the probability P of its occuring _k_ or more times
/// in _n_ trials is termed a cumulative binomial probability
/// the formula is P = sum_from_j=k_to_n(TMath::Binomial(n, j)*
/// *TMath::Power(p, j)*TMath::Power(1-p, n-j)
/// For _n_ larger than 12 BetaIncomplete is a much better way
/// to evaluate the sum than would be the straightforward sum calculation
/// for _n_ smaller than 12 either method is acceptable
/// ("Numerical Recipes")
///     --implementation by Anna Kreshuk

Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1);
/// Calculate a Breit Wigner function with mean and gamma.

Double_t CauchyDist(Double_t x, Double_t t=0, Double_t s=1);
/// Computes the density of Cauchy distribution at point x
/// by default, standard Cauchy distribution is used (t=0, s=1)
///    t is the location parameter
///    s is the scale parameter
/// The Cauchy distribution, also called Lorentzian distribution,
/// is a continuous distribution describing resonance behavior
/// The mean and standard deviation of the Cauchy distribution are undefined.
/// The practical meaning of this is that collecting 1,000 data points gives
/// no more accurate an estimate of the mean and standard deviation than
/// does a single point.
/// The formula was taken from "Engineering Statistics Handbook" on site
/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3663.htm
/// Implementation by Anna Kreshuk.

Double_t ChisquareQuantile(Double_t p, Double_t ndf);
/// Evaluate the quantiles of the chi-squared probability distribution function.
/// Algorithm AS 91   Appl. Statist. (1975) Vol.24, P.35
/// implemented by Anna Kreshuk.
/// Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
/// Parameters:
///   p   - the probability value, at which the quantile is computed
///   ndf - number of degrees of freedom

Double_t FDist(Double_t F, Double_t N, Double_t M);
/// Computes the density function of F-distribution
/// (probability function, integral of density, is computed in FDistI).
/// Parameters N and M stand for degrees of freedom of chi-squares
/// mentioned above parameter F is the actual variable x of the
/// density function p(x) and the point at which the density function
/// is calculated.
/// F-distribution arises in testing whether two random samples
/// have the same variance. It is the ratio of two chi-square
/// distributions, with N and M degrees of freedom respectively,
/// where each chi-square is first divided by it's number of degrees
/// of freedom.
/// Implementation by Anna Kreshuk.

Double_t FDistI(Double_t F, Double_t N, Double_t M);
/// Calculates the cumulative distribution function of F-distribution,
/// this function occurs in the statistical test of whether two observed
/// samples have the same variance. For this test a certain statistic F,
/// the ratio of observed dispersion of the first sample to that of the
/// second sample, is calculated. N and M stand for numbers of degrees
/// of freedom in the samples 1-FDistI() is the significance level at
/// which the hypothesis "1 has smaller variance than 2" can be rejected.
/// A small numerical value of 1 - FDistI() implies a very significant
/// rejection, in turn implying high confidence in the hypothesis
/// "1 has variance greater than 2".
/// Implementation by Anna Kreshuk.

Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE);
/// Calculate a gaussian function with mean and sigma.
/// If norm=kTRUE (default is kFALSE) the result is divided
/// by sqrt(2*Pi)*sigma.

Double_t KolmogorovProb(Double_t z);/// Calculates the Kolmogorov distribution function,
Double_t KolmogorovTest(Int_t na, const Double_t *a, Int_t nb, const Double_t *b, Option_t *option);
///  Statistical test whether two one-dimensional sets of points are compatible
///  with coming from the same parent distribution, using the Kolmogorov test.
///  That is, it is used to compare two experimental distributions of unbinned data.
///  Input:
///  a,b: One-dimensional arrays of length na, nb, respectively.
///       The elements of a and b must be given in ascending order.
///  option is a character string to specify options
///         "D" Put out a line of "Debug" printout
///         "M" Return the Maximum Kolmogorov distance instead of prob
///  Output:
/// The returned value prob is a calculated confidence level which gives a
/// statistical test for compatibility of a and b.
/// Values of prob close to zero are taken as indicating a small probability
/// of compatibility. For two point sets drawn randomly from the same parent
/// distribution, the value of prob should be uniformly distributed between
/// zero and one.
///   in case of error the function return -1
///   If the 2 sets have a different number of points, the minimum of
///   the two sets is used.
/// Method:
/// The Kolmogorov test is used. The test statistic is the maximum deviation
/// between the two integrated distribution functions, multiplied by the
/// normalizing factor (rdmax*sqrt(na*nb/(na+nb)).
///  Code adapted by Rene Brun from CERNLIB routine TKOLMO (Fred James)
///   (W.T. Eadie, D. Drijard, F.E. James, M. Roos and B. Sadoulet,
///      Statistical Methods in Experimental Physics, (North-Holland,
///      Amsterdam 1971) 269-271)
///  Method Improvement by Jason A Detwiler (JADetwiler@lbl.gov)
///  -----------------------------------------------------------
///   The nuts-and-bolts of the TMath::KolmogorovTest() algorithm is a for-loop
///   over the two sorted arrays a and b representing empirical distribution
///   functions. The for-loop handles 3 cases: when the next points to be
///   evaluated satisfy a>b, a<b, or a=b:
///      for (Int_t i=0;i<na+nb;i++) {
///         if (a[ia-1] < b[ib-1]) {
///            rdiff -= sa;
///            ia++;
///            if (ia > na) {ok = kTRUE; break;}
///         } else if (a[ia-1] > b[ib-1]) {
///            rdiff += sb;
///            ib++;
///            if (ib > nb) {ok = kTRUE; break;}
///         } else {
///            rdiff += sb - sa;
///            ia++;
///            ib++;
///            if (ia > na) {ok = kTRUE; break;}
///            if (ib > nb) {ok = kTRUE; break;}
///        }
///         rdmax = TMath::Max(rdmax,TMath::Abs(rdiff));
///      }
///   For the last case, a=b, the algorithm advances each array by one index in an
///   attempt to move through the equality. However, this is incorrect when one or
///   the other of a or b (or both) have a repeated value, call it x. For the KS
///   statistic to be computed properly, rdiff needs to be calculated after all of
///   the a and b at x have been tallied (this is due to the definition of the
///   empirical distribution function; another way to convince yourself that the
///   old CERNLIB method is wrong is that it implies that the function defined as the
///   difference between a and b is multi-valued at x -- besides being ugly, this
///   would invalidate Kolmogorov's theorem).
///   The solution is to just add while-loops into the equality-case handling to
///   perform the tally:
///         } else {
///            double x = a[ia-1];
///            while(a[ia-1] == x && ia <= na) {
///              rdiff -= sa;
///              ia++;
///            }
///            while(b[ib-1] == x && ib <= nb) {
///              rdiff += sb;
///              ib++;
///            }
///            if (ia > na) {ok = kTRUE; break;}
///            if (ib > nb) {ok = kTRUE; break;}
///         }
///  NOTE1
///  A good description of the Kolmogorov test can be seen at:
///    http://www.itl.nist.gov/div898/handbook/eda/section3/eda35g.htm

Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE);
/// The LANDAU function.
/// mu is a location parameter and correspond approximatly to the most probable value
/// and sigma is a scale parameter (not the sigma of the full distribution which is not defined)
/// Note that for mu=0 and sigma=1 (default values) the exact location of the maximum of the distribution
/// (most proble value) is at x = -0.22278
/// This function has been adapted from the CERNLIB routine G110 denlan.
/// If norm=kTRUE (default is kFALSE) the result is divided by sigma

Double_t LandauI(Double_t x);
///Returns the value of the Landau distribution function at point x.
///The algorithm was taken from the Cernlib function dislan(G110)
///Reference: K.S.Kolbig and B.Schorr, "A program package for the Landau
///distribution", Computer Phys.Comm., 31(1984), 97-111

Double_t LaplaceDist(Double_t x, Double_t alpha=0, Double_t beta=1);
/// Computes the probability density function of Laplace distribution
/// at point x, with location parameter alpha and shape parameter beta.
/// By default, alpha=0, beta=1
/// This distribution is known under different names, most common is
/// double exponential distribution, but it also appears as
/// the two-tailed exponential or the bilateral exponential distribution

Double_t LaplaceDistI(Double_t x, Double_t alpha=0, Double_t beta=1);
/// Computes the distribution function of Laplace distribution
/// at point x, with location parameter alpha and shape parameter beta.
/// By default, alpha=0, beta=1
/// This distribution is known under different names, most common is
/// double exponential distribution, but it also appears as
/// the two-tailed exponential or the bilateral exponential distribution

Double_t LogNormal(Double_t x, Double_t sigma, Double_t theta=0, Double_t m=1);
/// Computes the density of LogNormal distribution at point x.
/// Variable X has lognormal distribution if Y=Ln(X) has normal distribution
/// sigma is the shape parameter
/// theta is the location parameter
/// m is the scale parameter
/// The formula was taken from "Engineering Statistics Handbook" on site
/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm
/// Implementation using ROOT::Math::lognormal_pdf

Double_t NormQuantile(Double_t p);
/// Computes quantiles for standard normal distribution N(0, 1)
/// at probability p
/// ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3, 477-484.

Double_t Poisson(Double_t x, Double_t par);
/// compute the Poisson distribution function for (x,par)
/// The Poisson PDF is implemented by means of Euler's Gamma-function
/// (for the factorial), so for any x integer argument it is correct.
/// BUT for non-integer x values, it IS NOT equal to the Poisson distribution.
/// see TMath::PoissonI to get a non-smooth function.
/// Note that for large values of par, it is better to call
///     TMath::Gaus(x,par,sqrt(par),kTRUE)

Double_t PoissonI(Double_t x, Double_t par);
/// compute the Poisson distribution function for (x,par)
/// This is a non-smooth function.
/// This function is equivalent to ROOT::Math::poisson_pdf

Double_t Prob(Double_t chi2,Int_t ndf);
/// Computation of the probability for a certain Chi-squared (chi2)
/// and number of degrees of freedom (ndf).
/// Calculations are based on the incomplete gamma function P(a,x),
/// where a=ndf/2 and x=chi2/2.
/// P(a,x) represents the probability that the observed Chi-squared
/// for a correct model should be less than the value chi2.
/// The returned probability corresponds to 1-P(a,x),
/// which denotes the probability that an observed Chi-squared exceeds
/// the value chi2 by chance, even for a correct model.
///--- NvE 14-nov-1998 UU-SAP Utrecht

Double_t Student(Double_t T, Double_t ndf);
/// Computes density function for Student's t- distribution
/// (the probability function (integral of density) is computed in StudentI).
/// First parameter stands for x - the actual variable of the
/// density function p(x) and the point at which the density is calculated.
/// Second parameter stands for number of degrees of freedom.
/// Student's t-distribution is used for many significance tests, for example,
/// for the Student's t-tests for the statistical significance of difference
/// between two sample means and for confidence intervals for the difference
/// between two population means.
/// Example: suppose we have a random sample of size n drawn from normal
/// distribution with mean Mu and st.deviation Sigma. Then the variable
///   t = (sample_mean - Mu)/(sample_deviation / sqrt(n))
/// has Student's t-distribution with n-1 degrees of freedom.
/// NOTE that this function's second argument is number of degrees of freedom,
/// not the sample size.
/// As the number of degrees of freedom grows, t-distribution approaches
/// Normal(0,1) distribution.
/// Implementation by Anna Kreshuk.

Double_t StudentI(Double_t T, Double_t ndf);
/// Calculates the cumulative distribution function of Student's
/// t-distribution second parameter stands for number of degrees of freedom,
/// not for the number of samples
/// if x has Student's t-distribution, the function returns the probability of
/// x being less than T.
/// Implementation by Anna Kreshuk.

Double_t StudentQuantile(Double_t p, Double_t ndf, Bool_t lower_tail=kTRUE);
/// Computes quantiles of the Student's t-distribution
/// 1st argument is the probability, at which the quantile is computed
/// 2nd argument - the number of degrees of freedom of the
/// Student distribution
/// When the 3rd argument lower_tail is kTRUE (default)-
/// the algorithm returns such x0, that
///   P(x < x0)=p
/// upper tail (lower_tail is kFALSE)- the algorithm returns such x0, that
///   P(x > x0)=p
/// the algorithm was taken from
///   G.W.Hill, "Algorithm 396, Student's t-quantiles"
///             "Communications of the ACM", 13(10), October 1970

Double_t Vavilov(Double_t x, Double_t kappa, Double_t beta2);
///Returns the value of the Vavilov density function
///Parameters: 1st - the point were the density function is evaluated
///            2nd - value of kappa (distribution parameter)
///            3rd - value of beta2 (distribution parameter)
///The algorithm was taken from the CernLib function vavden(G115)
///Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
///Nucl.Instr. and Meth. B47(1990), 215-224
///Accuracy: quote from the reference above:
///"The resuls of our code have been compared with the values of the Vavilov
///density function computed numerically in an accurate way: our approximation
///shows a difference of less than 3% around the peak of the density function, slowly
///increasing going towards the extreme tails to the right and to the left"
///Begin_Html

Double_t VavilovI(Double_t x, Double_t kappa, Double_t beta2);
///Returns the value of the Vavilov distribution function
///Parameters: 1st - the point were the density function is evaluated
///            2nd - value of kappa (distribution parameter)
///            3rd - value of beta2 (distribution parameter)
///The algorithm was taken from the CernLib function vavden(G115)
///Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
///Nucl.Instr. and Meth. B47(1990), 215-224
///Accuracy: quote from the reference above:
///"The resuls of our code have been compared with the values of the Vavilov
///density function computed numerically in an accurate way: our approximation
///shows a difference of less than 3% around the peak of the density function, slowly
///increasing going towards the extreme tails to the right and to the left"

Double_t Voigt(Double_t x, Double_t sigma, Double_t lg, Int_t r = 4);
/// Computation of Voigt function (normalised).
/// Voigt is a convolution of
/// gauss(xx) = 1/(sqrt(2*pi)*sigma) * exp(xx*xx/(2*sigma*sigma)
/// and
/// lorentz(xx) = (1/pi) * (lg/2) / (xx*xx + lg*lg/4)
/// functions.
/// The Voigt function is known to be the real part of Faddeeva function also
/// called complex error function .
/// The algoritm was developed by J. Humlicek .
/// This code is based on fortran code presented by R. J. Wells .
/// Translated and adapted by Miha D. Puc
/// To calculate the Faddeeva function with relative error less than 10^(-r).
/// r can be set by the the user subject to the constraints 2 <= r <= 5.
///  J. Humlicek, JQSRT, 21, 437 (1982).
///  R.J. Wells "Rapid Approximation to the Voigt/Faddeeva Function and its
/// Derivatives" JQSRT 62 (1999), pp 29-48.
/// http://www-atm.physics.ox.ac.uk/user/wells/voigt.html

/* ************************** */
/* * Statistics over arrays * */
/* ************************** */

//Mean, Geometric Mean, Median, RMS(sigma)

template <typename T> Double_t Mean(Long64_t n, const T *a, const Double_t *w=0);
template <typename Iterator> Double_t Mean(Iterator first, Iterator last);
template <typename Iterator, typename WeightIterator> Double_t Mean(Iterator first, Iterator last, WeightIterator wfirst);

template <typename T> Double_t GeomMean(Long64_t n, const T *a);
template <typename Iterator> Double_t GeomMean(Iterator first, Iterator last);

template <typename T> Double_t RMS(Long64_t n, const T *a, const Double_t *w=0);
template <typename Iterator> Double_t RMS(Iterator first, Iterator last);
template <typename Iterator, typename WeightIterator> Double_t RMS(Iterator first, Iterator last, WeightIterator wfirst);

template <typename T> Double_t StdDev(Long64_t n, const T *a, const Double_t * w = 0) { return RMS<T>(n,a,w); }
template <typename Iterator> Double_t StdDev(Iterator first, Iterator last) { return RMS<Iterator>(first,last); }
template <typename Iterator, typename WeightIterator> Double_t StdDev(Iterator first, Iterator last, WeightIterator wfirst) { return RMS<Iterator,WeightIterator>(first,last,wfirst); }

template <typename T> Double_t Median(Long64_t n, const T *a,  const Double_t *w=0, Long64_t *work=0);

//k-th order statistic
template <class Element, typename Size> Element KOrdStat(Size n, const Element *a, Size k, Size *work = 0);

/* ******************* */
/* * Special Functions */
/* ******************* */

Double_t Beta(Double_t p, Double_t q);
/// Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).

Double_t BetaCf(Double_t x, Double_t a, Double_t b);
/// Continued fraction evaluation by modified Lentz's method
/// used in calculation of incomplete Beta function.

Double_t BetaDist(Double_t x, Double_t p, Double_t q);
/// Computes the probability density function of the Beta distribution
/// (the distribution function is computed in BetaDistI).
/// The first argument is the point, where the function will be
/// computed, second and third are the function parameters.
/// Since the Beta distribution is bounded on both sides, it's often
/// used to represent processes with natural lower and upper limits.

Double_t BetaDistI(Double_t x, Double_t p, Double_t q);
/// Computes the distribution function of the Beta distribution.
/// The first argument is the point, where the function will be
/// computed, second and third are the function parameters.
/// Since the Beta distribution is bounded on both sides, it's often
/// used to represent processes with natural lower and upper limits.

Double_t BetaIncomplete(Double_t x, Double_t a, Double_t b);/// Calculates the incomplete Beta-function.

// Bessel functions
Double_t BesselI(Int_t n,Double_t x);  // integer order modified Bessel function I_n(x)
Double_t BesselK(Int_t n,Double_t x);  // integer order modified Bessel function K_n(x)
Double_t BesselI0(Double_t x);         // modified Bessel function I_0(x)
Double_t BesselK0(Double_t x);         // modified Bessel function K_0(x)
Double_t BesselI1(Double_t x);         // modified Bessel function I_1(x)
Double_t BesselK1(Double_t x);         // modified Bessel function K_1(x)
Double_t BesselJ0(Double_t x);         // Bessel function J0(x) for any real x
Double_t BesselJ1(Double_t x);         // Bessel function J1(x) for any real x
Double_t BesselY0(Double_t x);         // Bessel function Y0(x) for positive x
Double_t BesselY1(Double_t x);         // Bessel function Y1(x) for positive x
Double_t StruveH0(Double_t x);         // Struve functions of order 0
Double_t StruveH1(Double_t x);         // Struve functions of order 1
Double_t StruveL0(Double_t x);         // Modified Struve functions of order 0
Double_t StruveL1(Double_t x);         // Modified Struve functions of order 1

Double_t DiLog(Double_t x);/// The DiLogarithm function
Double_t Erf(Double_t x);
/// Computation of the error function erf(x).
/// Erf(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between 0 and x

Double_t ErfInverse(Double_t x);
/// returns  the inverse error function
/// x must be  <-1<x<1

Double_t Erfc(Double_t x);
/// Compute the complementary error function erfc(x).
/// Erfc(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between x and infinity

Double_t ErfcInverse(Double_t x);
// returns  the inverse of the complementary error function
// x must be  0<x<2
// implement using  the quantile of the normal distribution
// instead of ErfInverse for better numerical precision for large x
// erfc-1(x) = - 1/sqrt(2) * normal_quantile( 0.5 * x)

Double_t Freq(Double_t x);
/// Computation of the normal frequency function freq(x).
/// Freq(x) = (1/sqrt(2pi)) Integral(exp(-t^2/2))dt between -infinity and x.

Double_t Gamma(Double_t z);/// Computation of gamma(z) for all z.
Double_t Gamma(Double_t a,Double_t x);
/// Computation of the normalized lower incomplete gamma function P(a,x) as defined in the
/// Handbook of Mathematical Functions by Abramowitz and Stegun, formula 6.5.1 on page 260 .
/// Its normalization is such that TMath::Gamma(a,+infinity) = 1 .

Double_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1);
/// Computes the density function of Gamma distribution at point x.
///   gamma - shape parameter
///   mu    - location parameter
///   beta  - scale parameter
/// The definition can be found in "Engineering Statistics Handbook" on site
/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
/// use now implementation in ROOT::Math::gamma_pdf

Double_t LnGamma(Double_t z);
/// Computation of ln[gamma(z)] for all z.
/// C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
/// The accuracy of the result is better than 2e-10.
///--- Nve 14-nov-1998 UU-SAP Utrecht

   Double_t GamCf(Double_t a,Double_t x);/// Computation of the incomplete gamma function P(a,x) via its continued fraction representation.
Double_t GamSer(Double_t a,Double_t x);/// Computation of the incomplete gamma function P(a,x) via its series representation.
Double_t VavilovDenEval(Double_t rlam, Double_t *AC, Double_t *HC, Int_t itype);///Internal function, called by Vavilov and VavilovSet
void VavilovSet(Double_t rkappa, Double_t beta2, Bool_t mode, Double_t *WCM, Double_t *AC, Double_t *HC, Int_t &itype, Int_t &npt);///Internal function, called by Vavilov and VavilovI


## code

TF1* fc = new TF1("fc", "TMath::CauchyDist(x, , )", -5, 5);
fc->SetParameters(0, 1);
fc->Draw();


## example

updated: 2018-05-14 17:18:22