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; }
// Avogadro constant (Avogadro's Number)
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[0] corresponding to the largest arr1 value and
/// arr2[Narr] the smallest.
/// Author: Adrian Bevan (bevan@slac.stanford.edu)
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[3],const T v2[3], T out[3]);
Float_t Normalize(Float_t v[3]); // Normalize a vector
/// Normalize a vector v in place.
/// Returns the norm of the original vector.
Double_t Normalize(Double_t v[3]); // 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[3],const T v2[3],T out[3]);
// Calculate a normal vector of a plane
template <typename T> T *Normal2Plane(const T v1[3],const T v2[3],const T v3[3], T normal[3]);
/* ************************ */
/* * Polynomial Functions * */
/* ************************ */
Bool_t RootsCubic(const Double_t coef[4],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[3], b == coef[2], c == coef[1], d == coef[0]
///coef[3] 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.
/// About F distribution:
/// 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.
/// About Student distribution:
/// 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 [2].
/// The algoritm was developed by J. Humlicek [1].
/// This code is based on fortran code presented by R. J. Wells [2].
/// 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.
/// [1] J. Humlicek, JQSRT, 21, 437 (1982).
/// [2] 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, [0], [1])", -5, 5);
fc->SetParameters(0, 1);
fc->Draw();