Mathematical libraries

The ROOT mathematical libraries consist of the following components:

MathCore library

The MathCore library provides a collection of functions, C++ classes and ROOT classes for HEP numerical computing.
The MathCore is a self-consistent minimal set of tools needed for the basic numerical computing. More advanced mathematical functionalities are provided by the MathMore library. The following is included in the MathCore library:

In addition, the MathCore library contains the following ROOT classes that were originally part of libCore:

TMath

The TMath namespace provides a collection of free functions:

Elementary functions

Some of elementary mathematical functions refer to basic mathematical functions such as the square root, the power to a number of the calculus of a logarithm, while others are used for number treatment, like rounding.

Although there are some functions that are not in the standard C math library (such as Factorial), most of the functionality offered here is just a wrapper of the first ones. Nevertheless, some of them also offer some security checks or a better precision, such as the trigonometrical functions ASin(x), ACos(x) or ATan(x).

Examples

// Generate a vector with 10 random numbers.
vector<double> v(10);
std::generate(v.begin(), v.end(), rand);

// Find the minimum value of the vector (iterator version).
vector<double>::iterator it;
it = TMath::LocMin(v.begin(), v.end());
std::cout << *it << std::endl;

// The same with the old-style version.
int i;
i = TMath::LocMin(10, &v[0]);
std::cout << v[i] << std::endl;

Statistic functions operating on arrays

TMath provides functions that process arrays for calculation:

  • mean
  • median
  • geometrical mean
  • sample standard deviation (RMS)
  • the kth smallest element

Example

// Size of the array.
const int n = 100;

// Vector v with random values.
vector<double> v(n);
std::generate(v.begin(), v.end(), rand);

// Weight vector w.
vector<double> w(n);
std::fill(w.begin(), w.end, 1);
double mean;

// Calculate the mean of the vector with iterators.
mean = TMath::Mean(v.begin(), v.end());

Special and statistical functions

TMath provides special functions such as Bessel, error functions, Gamma or similar statistical mathematical functions, including probability density functions, cumulative distribution and its inverse.

The majority of the special functions and the statistical distributions are provided also as free functions in the ROOT::Math namespace.

Functions not present in the ROOT::Math name that are provided only by TMath are:

  • Special functions:
    • DiLogarithm
    • Struve
  • Statistical functions:
    • KolmogorovProb
    • Voigt function
    • LaplaceDist
    • Vavilov

ROOT::Math

The ROOT::Math namespace provides a set of function interfaces to define the basic behaviour of a mathematical function:

In addition, helper classes, wrapping the user interfaces in the ROOT::Math function interfaces are provided. With wrapper functions you can insert your own type of function in the needed function interface.

To use the self-defined functions, they must have inherited from one of the following classes:

Figure: ROOT::Math function interface structure.

One-dimensional function interfaces

This interface is used for numerical algorithms operating only on one-dimensional functions. It cannot applied to multi-dimensional functions.

ROOT::Math::IBaseFunctionOneDim
This interface provides a method to evaluate the function given a value (simple double) by implementing double operator() (const double). The defined user class only needs to reimplement the purely abstract double DoEval(double x) method, which does the work of evaluating the function at point x.

Example

Example for the implementation of a class that represents a mathematical function.

#include "Math/IFunction.h"

class MyFunction: public ROOT::Math::IBaseFunctionOneDim
{
    double DoEval(double x) const
    {
        return x*x;
    }
    ROOT::Math::IBaseFunctionOneDim* Clone() const
    {
        return new MyFunction();
    }
};

ROOT::Math::IGradientFunctionOneDim
This interface is needed by some numerical algorithms to calculate the derivatives of the function. It introduces the method double Derivative(double x), which returns the derivative of the function at point x. The class from which the user inherits must implement the abstract method double DoDerivative(double x), leaving the rest of the class untouched.

Example

Example for the implementation of a gradient one-dimensional function.

#include "Math/IFunction.h"

class MyGradientFunction: public ROOT::Math::IGradientFunctionOneDim
{
public:
    double DoEval(double x) const
    {
        return sin(x);
    }
    ROOT::Math::IBaseFunctionOneDim* Clone() const
    {
        return new MyGradientFunction();
    }
    double DoDerivative(double x) const
    {
        return -cos(x);
    }
};

Multi-dimensional function interfaces

This interface is used for numerical algorithms operating on multi-dimensional functions.

ROOT::Math::IBaseFunctionMultiDim
This interface provides the double operator() (const double*) that takes an array of doubles with all the values for the different dimensions. In this case, the user has to provide the functionality for two different functions: double DoEval(const double*) and the unsigned int NDim(). The first evaluates the function given the array representing the multiple variables. The second returns the number of dimensions of the function.

Example

Example for the implementation of a basic multi-dimensional function.

#include "Math/IFunction.h"

class MyFunction: public ROOT::Math::IBaseFunctionMultiDim
{
public:
    double DoEval(const double* x) const
    {
       return x[0] + sin(x[1]);
    }
    unsigned int NDim() const
    {
       return 2;
    }
    ROOT::Math::IBaseFunctionMultiDim* Clone() const
    {
       return new MyFunction();
    }
};

ROOT::Math::IGradientFunctionMultiDim
This interface offers the same functionality as the base function and additionally the calculation of the derivative. It only adds the double DoDerivative(double* x, uint ivar) method for the user to implement. This method must implement the derivative of the function with respect to the variable in the first parameter array at the position indicated by the second parameter.

Example

Example for the implementation of a multi-dimensional gradient function.

#include "Math/IFunction.h"

class MyGradientFunction: public ROOT::Math::IGradientFunctionMultiDim
{
public:
    double DoEval(const double* x) const
    {
        return x[0] + sin(x[1]);
    }
    unsigned int NDim() const
    {
        return 2;
    }
    ROOT::Math::IGradientFunctionMultiDim* Clone() const
    {
        return new MyGradientFunction();
    }
    double DoDerivative(const double* x, unsigned int ipar) const
    {
        if ( ipar == 0 )
            return sin(x[1]);
        else
            return x[0] + x[1] * cos(x[1]);
    }
};

Parametric function interfaces

This interface is used for fitting after evaluating multi-dimensional functions.

ROOT::Math::IParametricFunctionMultiDim
This interface describes a multi-dimensional parametric function. The user needs to provide the void SetParameters(double* p) method as well as the getter methods const double * Parameters() and uint NPar().

Example

Example for the implementation of a parametric function.

#include "Math/IFunction.h"
#include "Math/IParamFunction.h"

class MyParametricFunction: public ROOT::Math::IParametricFunctionMultiDim
{
private:
    const double* pars;
public:
    double DoEvalPar(const double* x, const double* p) const
    {
        return p[0] * x[0] + sin(x[1]) + p[1];
    }
    unsigned int NDim() const
    {
        return 2;
    }
    ROOT::Math::IParametricFunctionMultiDim* Clone() const
    {
        return new MyParametricFunction();
    }
    const double* Parameters() const
    {
        return pars;
    }
    void SetParameters(const double* p)
    {
        pars = p;
    }
    unsigned int NPar() const
    {
        return 2;
    }
};

ROOT::Math::IParametricGradFunctionMultiDim
This interface provides an interface for parametric gradient multi-dimensional functions. In addition to function evaluation, it provides the gradient with respect to the parameters, via the ParameterGradient() method. This interface is only used in case of some dedicated fitting algorithms, when is required or more efficient to provide derivatives with respect to the parameters.

Example

Example for the implementation of a parametric gradient function.

#include "Math/IFunction.h"
#include "Math/IParamFunction.h"

class MyParametricGradFunction: public ROOT::Math::IParametricGradFunctionMultiDim
{
private:
    const double* pars;
public:
    double DoEvalPar(const double* x, const double* p) const
    {
        return p[0] * x[0] + sin(x[1]) + p[1];
    }
    unsigned int NDim() const
    {
       return 2;
    }
    ROOT::Math::IParametricGradFunctionMultiDim* Clone() const
    {
        return new MyParametricGradFunction();
    }
    const double* Parameters() const
    {
        return pars;
    }
    void SetParameters(const double* p)
    {
        pars = p;
    }
    unsigned int NPar() const
    {
        return 2;
    }
    double DoParameterDerivative(const double* x,
                                 const double* p, unsigned int ipar) const
    {
        if ( ipar == 0 )
            return sin(x[1]) + p[1];
        else
           return p[0] * x[0] + x[1] * cos(x[1]) + p[1];
    }
};

Wrapper functions

To insert your own type of function in the needed function interface, helper classes, wrapping the user interface in the ROOT::Math function interfaces are provided.

There is one possible wrapper for every interface.

Interface Wrapper Description
ROOT::Math::IBaseFunctionOneDim ROOT::Math::Functor1D See Wrapping one-dimensional functions
ROOT::Math::IGradientFunctionOneDim ROOT::Math::GradFunctor1D See Wrapping one-dimensional gradient functions
ROOT::Math::IBaseFunctionMultiDim ROOT::Math::Functor See Wrapping multi-dimensional functions
ROOT::Math::IGradientFunctionMultiDim ROOT::Math::GradFunctor See Wrapping multi-dimensional gradient functions

Note the special case when wrapping TF1 objects in parametric function interfaces.

Wrapping one-dimensional functions

Use ROOT::Math::Functor1D to wrap one-dimensional functions.

ROOT::Math::Functor1D can wrap the following types:

  • A free C function of type double ()(double ).
  • Any C++ callable object implementation double operator()( double).
  • A class member function with the correct signature like double Foo::Eval(double ). In this case you pass the object pointer and a pointer to the member function (&Foo::Eval).

Example

#include "Math/Functor.h"

class MyFunction1D {
public:
    double operator()(double x) const
    {
        return x*x;
    }
    double Eval(double x) const
    {
        return x+x;
    }
};

double freeFunction1D(double x ) {
    return 2*x;
}

int main() {

    // Wrapping a free function.
    ROOT::Math::Functor1D f1(&freeFunction1D);
    MyFunction1D myf1;

    // Wrapping a function object implementing operator().
    ROOT::Math::Functor1D f2(myf1);

    // Wrapping a class member function.
    ROOT::Math::Functor1D f3(&myf1,&MyFunction1D::Eval);
    cout << f1(2) << endl;
    cout << f2(2) << endl;
    cout << f3(2) << endl;
    return 0;
}

Wrapping one-dimensional gradient functions

Use ROOT::Math::GradFunctor1D to wrap one-dimensional gradient functions.

It can be constructed in three different ways:

  • Any object implementing both double operator()( double) for the function evaluation and double Derivative(double) for the function derivative.
  • Any object implementing any member function such as Foo::XXX(double ) for the function evaluation and any other member function such as Foo::YYY(double) for the derivative.
  • Any two function objects implementing double operator()( double). One object provides the function evaluation, the other the derivative. One or both function objects can be a free C function of type double ()(double).

Wrapping multi-dimensional functions

Use the ROOT::Math::Functor to wrap multi-dimensional function objects.

It can wrap all the following types:

  • Any C++ callable object implementing double operator()( const double * ).
  • A free C function of type double ()(const double *).
  • A member function with the correct signature like Foo::Eval(const double *). In this case one pass the object pointer and a pointer to the member function (&Foo::Eval).

Example

#include "Math/Functor.h"

class MyFunction {
public:
    double operator()(const double *x) const
    {
        return x[0]+x[1];
    }
    double Eval(const double * x) const
    {
        return x[0]+x[1];
    }
};

double freeFunction(const double * x ) {
    return x[0]+x[1];
}

int main() {
    // Test directly calling the function object.
    MyFunction myf;

    // Test from a free function pointer.
    ROOT::Math::Functor f1(&freeFunction,2);

    // Test from function object.
    ROOT::Math::Functor f2(myf,2);

    // Test from a member function.
    ROOT::Math::Functor f3(&myf,&MyFunction::Eval,2);
    double x[] = {1,2};
    cout << f1(x) << endl;
    cout << f2(x) << endl;
    cout << f3(x) << endl;
    return 0;
}

Wrapping multi-dimensional gradient functions

Use ROOT::Math::GradFunctor to wrap C++ callable objects to make gradient functions.

It can be constructed in three different ways:

  • From an object implementing both double operator()( const double*) for the function evaluation and double Derivative(const double *, int icoord) for the partial derivatives.
  • From an object implementing any member function such as Foo::XXX(const double *) for the function evaluation and any member function such as Foo::XXX(const double *, int icoord) for the partial derivatives.
  • From a function object implementing double operator()( const double *) for the function evaluation and another function object implementing double operator() (const double *, int icoord) for the partial derivatives.

The function dimension is required when constructing the functor.

Wrapping TF1 objects in parametric function interfaces

Often the TF1 class is used.
Use the ROOT::Math::WrappedTF1 class, if the interface to be wrapped is one-dimensional.

The default constructor takes a TF1 reference as argument, wrapped with the interfaces of the ROOT::Math::IParametricGradFunctionOneDim class.

Example

#include "TF1.h"
#include "Math/WrappedTF1.h"
int main()
{
    TF1 f("Sin Function", "sin(x)+y",0,3);
    ROOT::Math::WrappedTF1 wf1(f);
    cout << f(1) << endl;
    cout << wf1(1) << endl;
    return 0;
}

Use the ROOT::Math::WrappedMultiTF1 class, if the interface to be wrapped is multi-dimensional.

Following the usual procedure, setting the TF1 though the constructor, wraps it into a ROOT::Math::IParametricGradFunctionMultiDim.

Example

#include "TF1.h"
#include "Math/WrappedMultiTF1.h"
int main()
{
    TF2 f("Sin Function", "sin(x) + y",0,3,0,2);
    ROOT::Math::WrappedMultiTF1 wf1(f);
    double x[] = {1,2};
    cout << f(x) << endl;
    cout << wf1(x) << endl;
    return 0;
}

Random numbers

The MathCore library provides the following classes for generating pseudo-random numbers:

  • TRandom: Using a linear congruential random generator.
  • TRandom1: Random number generator based on the Ranlux engine.
  • TRandom2: Based on the maximally equi-distributed combined Tausworthe generator by L’Ecuyer.
  • TRandom3: Based on the Mersenne and Twister pseudo-random number generator.

Note

For generating non-uniform random numbers, the UNU.RAN package (see UNU.RAN) is available.

You can work with the random number generators as follows:

Seeding the random number generators

The SetSeed() method allows to set the seed of a random generator object. When no value is given, the default seed of the generator is used. In this case, an identical sequence is generated each time the application is run. Calling SetSeed(0) generates a unique seed using either:

  • A TUUID class instance, for TRandom1, TRandom2 and TRandom3.
  • The machine clock, for TRandom. Note that in this case the resolution is about 1 s. Therefore, identical sequences are generated when the elapsed time is less than one second.

Using the random number generators

The Rndm() method generates a pseudo-random number distributed between 0 and 1.

Example

// Use the default seed (same random numbers will be generated each time).
// Generate a number in the interval ]0,1] (0 is excluded).
TRandom3 r;
r.Rndm();
double x[100];

// Generate an array of random numbers in ]0,1].
r.RndmArray(100,x);

// Construct with a user-defined seed.
TRandom3 rdm(111);

// Use 0: a unique seed is automatically generated with TUUID.
TRandom1 r1(0);
TRandom2 r2(0);
TRandom3 r3(0);

// Seed generated using machine clock (different every second).
TRandom r0(0);

Random number distributions

The TRandom class provides functions that can be used by all other derived classes to generate random variables according to predefined distributions. In the simplest cases, as in the exponential distribution, the non-uniform random number is obtained by suitable transformations. In the more complicated cases, the random variables are obtained by acceptance-rejection methods that require several random numbers.

Example

TRandom3 r;
// Generate a gaussian distributed number with:
// mu=0, sigma=1 (default values)
double x1 = r.Gaus();
double x2 = r.Gaus(10,3);
// Use mu = 10, sigma = 3;

The following table shows the various distributions that can be generated using methods of the TRandom classes. In addition, you can use TF1::GetRandom() or TH1::GetRandom() to generate random numbers distributed according to a user defined function, in a limited interval, or to a user defined histogram.

Distributions Description
Double_t Uniform(Double_t x1,Double_t x2) Uniform random numbers between x1,x2.
Double_t Gaus(Double_t mu,Double_t sigma) Gaussian random numbers. Default values: mu=0, sigma=1.
Double_t Exp(Double_t tau) Exponential random numbers with mean tau.
Double_t Landau(Double_t mean,Double_t sigma) Landau distributed random numbers. Default values: mean=0, sigma=1.
Double_t BreitWigner(Double_t mean,Double_t gamma) Breit-Wigner distributed random numbers. Default values mean=0, gamma=1.
Int_t Poisson(Double_t mean) Poisson random numbers.
Double_t PoissonD(Double_t mean) Poisson random numbers.
Int_t Binomial(Int_t ntot,Double_t prob) Binomial Random numbers
Circle(Double_t &x,Double_t &y,Double_t r) Generate a random 2D point (x,y) in a circle of radius r.
Sphere(Double_t &x,Double_t &y,Double_t &z,Double_t r) Generate a random 3D point (x,y,z) in a sphere of radius r.
Rannor(Double_t &a,Double_t &b) Generate a pair of Gaussian random numbers with mu=0 and sigma=1.

Complex numbers

The MathCore library provides with TComplex a class for complex numbers.

Numerical integration

ROOT provides algorithms for integration of one-dimensional functions, with several adaptive and non-adaptive methods and for integration of multi-dimensional function using an adaptive method or MonteCarlo Integration (GSLMCIntegrator).

ROOT::Math::VirtualIntegrator defines the most basic functionality, this is, the common methods for the numerical integrator classes of one- and multi-dimensions.

ROOT::Math::VirtualIntegratorOneDim is an abstract interface class for 1Dnumerical integration. This method must be implemented in concrete classes, so you must create the ROOT::Math::IntegratorOneDim class for integrating one-dimensional functions.

ROOT::Math::VirtualIntegratorMultiDim is an abstract interface class for multi-numerical integration. This method must be implemented in concrete classes, so you must create the ROOT::Math::IntegratorMultiDim class for integrating multi-dimensional functions.

Using ROOT::Math::IntegratorOneDim

The following code example shows how to use ROOT::Math::IntegratorOneDim.

Example

In this example different instances of the class are created using some of the available algorithms in ROOT. If no algorithm is specified, the default one is used. The default integrator together with other integration options, such as relative and absolute tolerance, can be specified using the static method of the ROOT::Math::IntegratorOneDimOptions.

#include "Math/Integrator.h"
const double ERRORLIMIT = 1E-3;

double f(double x) {
    return x;
}
double f2(const double * x) {
    return x[0] + x[1];
}

int testIntegration1D() {
    const double RESULT = 0.5;
    int status = 0;

    // Set default tolerances for all integrators.
    ROOT::Math::IntegratorOneDimOptions::SetDefaultAbsTolerance(1.E-6);
    ROOT::Math::IntegratorOneDimOptions::SetDefaultRelTolerance(1.E-6);

    ROOT::Math::Functor1D wf(&f);
    ROOT::Math::Integrator ig(ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR);
    ig.SetFunction(wf);
    double val = ig.Integral(0,1);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    ROOT::Math::Integrator ig2(ROOT::Math::IntegrationOneDim::kNONADAPTIVE);
    ig2.SetFunction(wf);
    val = ig2.Integral(0,1);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    ROOT::Math::Integrator ig3(wf, ROOT::Math::IntegrationOneDim::kADAPTIVE);
    val = ig3.Integral(0,1);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    ROOT::Math::Integrator ig4(ROOT::Math::IntegrationOneDim::kGAUSS)
    ig4.SetFunction(wf);
    val = ig4.Integral(0,1);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    ROOT::Math::Integrator ig4(ROOT::Math::IntegrationOneDim::kLEGENDRE);
    ig4.SetFunction(wf);
    val = ig4.Integral(0,1);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    return status;
}

Using ROOT::Math::IntegratorMultiDim

The following code example shows how to use ROOT::Math::IntegratorMultiDim.

Example

In this example, different instances of the class use some of the algorithms available in ROOT.

#include "Math/IntegratorMultiDim.h"
#include "Math/Functor.h"

double f2(const double * x) {
    return x[0] + x[1];
}

int testIntegrationMultiDim() {
    const double RESULT = 1.0;
    const double ERRORLIMIT = 1E-3;
    int status = 0;

    ROOT::Math::Functor wf(&f2,2);
    double a[2] = {0,0};
    double b[2] = {1,1};

    ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
    ig.SetFunction(wf);
    double val = ig.Integral(a,b);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    ROOT::Math::IntegratorMultiDim ig2(ROOT::Math::IntegrationMultiDim::kVEGAS);
    ig2.SetFunction(wf);
    val = ig2.Integral(a,b);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    ROOT::Math::IntegratorMultiDim ig3(wf,ROOT::Math::IntegrationMultiDim::kPLAIN);
    val = ig3.Integral(a,b);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    ROOT::Math::IntegratorMultiDim ig4(wf,ROOT::Math::IntegrationMultiDim::kMISER);
    val = ig4.Integral(a,b);
    std::cout << "integral result is " << val << std::endl;
    status += std::fabs(val-RESULT) > ERRORLIMIT;

    return status;
}

One-dimensional integration algorithms

You can instantiate one-dimensional integration algorithms by using the following enumeration values:

Enumeration name Integrator class
ROOT::Math::IntegratorOneDim::kGAUSS ROOT::Math::GaussianIntegrator
ROOT::Math::IntegratorOneDim::kLEGENDRE ROOT::Math:::GausLegendreIntegrator
ROOT::Math::Integration::kNONADAPTIVE ROOT::Math:::GSLIntegrator
ROOT::Math::Integration::kADAPTIVE ROOT::Math:::GSLIntegrator
ROOT::Math::Integration::kADAPTIVESINGULAR ROOT::Math:::GSLIntegrator

ROOT::Math:::GaussIntegrator

ROOT::Math:::GaussIntegrator uses the most basic Gaussian integration algorithm. It uses the 8-point and the 16-point Gaussian quadrature approximations.

Example

#include "TF1.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"

int main() {
    TF1 f("Sin Function", "sin(x)", 0, TMath::Pi());
    ROOT::Math::WrappedTF1 wf1(f);
    ROOT::Math::GaussIntegrator ig;
    ig.SetFunction(wf1, false);
    ig.SetRelTolerance(0.001);
    cout << ig.Integral(0, TMath::PiOver2()) << endl;
}

ROOT::Math::GaussLegendreIntegrator

ROOT::Math::GaussLegendreIntegrator implementes the Gauss-Legendre quadrature formulas. This sort of numerical methods requieres that you specify the number of intermediate function points used in the calculation of the integral. It automatically determines the coordinates and weights of such points before performing the integration. You can use the example above, but replacing the creation of a ROOT::Math:::GaussIntegrator object with ROOT::Math::GaussLegendreIntegrator.

ROOT::Math::GSLIntegrator

ROOT::Math::GSLIntegrator isa wrapper for the QUADPACK integrator implemented in the GSL library. It supports several integration methods that can be chosen in construction time. The default type is adaptive integration with singularity applying a Gauss-Kronrod 21-point integration rule.

Multi-dimensional integration algorithms

You can instantiate multi-dimensional integration algorithms by using the following enumeration values:

Enumeration name Integrator class
ROOT::Math::IntegratorMultiDim::kADAPTIVE ROOT::Math::AdaptiveIntegratorMultiDim
ROOT::Math::IntegratorMultiDim::kVEGAS ROOT::Math:::GSLMCIntegrator
ROOT::Math::IntegratorMultiDim::kMISER ROOT::Math:::GSLMCIntegrator
ROOT::Math::IntegratorMultiDim::kPLAIN ROOT::Math:::GSLMCIntegrator

ROOT::Math::AdaptiveIntegratorMultiDim

ROOT::Math::AdaptiveIntegratorMultiDim implements an adaptive quadrature integration method for multi-dimensional functions. It is described in the paper Genz, A.A. Malik, An adaptive algorithm for numerical integration over an N-dimensional rectangular region, J. Comput. Appl. Math. 6 (1980) 295-302.

ROOT::Math::GSLMCIntegrator

ROOT::Math::GSLMCIntegrator is a class for performing numerical integration of a multidimensional function. It uses the numerical integration algorithms of GSL, which reimplements the algorithms used in the QUADPACK, a numerical integration package written in Fortran. Plain MC, MISER and VEGAS integration algorithms are supported for integration over finite (hypercubic) ranges.

MathMore library

The MathMore library provides an advanced collection of functions and C++ classes for numerical computing. This is an extension of the functionality provided by the MathCore library. The MathMore library is implemented wrapping in C++ the GNU Scientific Library (GSL). The mathematical functions are implemented as a set of free functions in the namespace ROOT::Math.

The MathMore library includes classes and functions for:

Linear algebra packages

The linear algebra packages provide a complete environment in ROOT to perform calculations such as equation solving and eigenvalue decompositions.

There are the following linear algebra packages available:

Matrix package

The following topics are covered for the matrix package:

Matrix classes

ROOT provides the following matrix classes, among others:

  • TMatrixTBase: base class for matrices.
  • TMatrixT: template class of a general matrix. Specialized versions are available (e.g. TMatrixF for float precision).
  • TMatrixTSym: template class of a symmetric matrix. Specialized versions are available (e.g. TMatrixFSym for float precision).
  • TMatrixTSparse: template class of a general sparse matrix in the Harwell-Boeing format. Specialized versions are available (e.g. TMatrixFSparse for float precision).
  • TVectorT: template class for vectors in the linear algebra package. Specialized versions are available (e.g. TVectorF for float precision).
  • TDecompBase: Decomposition base class.
  • TDecompChol: Cholesky decomposition class.

Matrix properties

A matrix has five properties, which are all set in the constructor:

  • precision
    If the precision is float (this is single precision), use the TMatrixF class family. If the precision is double, use the TMatrixD class family.

  • type
    Possible values are: general (TMatrixD), symmetric (TMatrixDSym) or sparse (TMatrixDSparse).

  • size
    Number of rows and columns.

  • index
    Range start of row and column index. By default these start at 0.

  • sparse map
    Only relevant for a sparse matrix. It indicates where elements are unequal 0.

Accessing matrix properties

Use one of the following methods to access the information about the relevant matrix property:

*GetRowIndexArray() and *GetColIndexArray() are specific to the sparse matrix, which is implemented according to the Harwell- Boeing format. Here, besides the usual shape/size descriptors of the matrix such as fNrows, fRowLwb, fNcols and fColLwb, also a row index fRowIndex and a column index fColIndex are stored:

  • fRowIndex[0,..,fNrows]: Stores for each row the index range of the elements in the data and column array.
  • fColIndex[0,..,fNelems-1]: Stores the column number for each data element != 0.

For example, printing all non-zero elements of a matrix would look like:

Example

TMatrixDSparse a;
const Int_t *rIndex = a.GetRowIndexArray();
const Int_t *cIndex = a.GetColIndexArray();
const Double_t *pData = a.GetMatrixArray();

for (Int_t irow = 0; irow < a.getNrows(); irow++) {
    const Int_t sIndex = rIndex[irow];
    const Int_t eIndex = rIndex[irow+1];

    for (Int_t index = sIndex; index < eIndex; index++) {
        const Int_t icol = cIndex[index];
        const Double_t data = pData[index];
        printf("data(%d,%d) = %.4en", irow+a.GetfRowLwb(), icol+a.GetColLwb(), data);
    }
}

Setting matrix properties

Use one of the following methods to set a matrix property:

Creating and filling a matrix

A full list of constructors for matrices is available on the corresponding reference guide pages of TMatrixT, TMatrixTSparse and TMatrixTSym.

Use one of the following methods to fill a matrix:

  • SetMatrixArray(const Double_t*data,Option_t*option=""): copies array data. If option="F", the array fills the matrix column-wise else row-wise. This option is implemented for TMatrixD and TMatrixDSym. It is expected that the array data contains at least fNelems entries.
  • SetMatrixArray(Int_t nr,Int_t *irow,Int_t *icol,Double_t *data): only available for sparse matrices. The three arrays should each contain nr entries with row index, column index and data entry. Only the entries with non-zero data value are inserted.
  • operator(), operator[]: these operators provide the easiest way to fill a matrix but are in particular for a sparse matrix expensive. If no entry for slot (i,j) is found in the sparse index table, it is entered, which involves some memory management. Therefore, before invoking this method in a loop set the index table first through a call to the SetSparseIndex() method.
  • SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixDBase &source): the matrix to be inserted at position (row_lwb,col_lwb) can be both, dense or sparse.
  • Use(): allows inserting another matrix or data array without actually copying the data.

Example

A Hilbert matrix is created by copying an array.

TMatrixD h(5,5);
TArrayD data(25);
for (Int_t = 0; i < 25; i++) {
    const Int_t ir = i/5;
    const Int_t ic = i%5;
    data[i] = 1./(ir+ic);
}
h.SetMatrixArray(data.GetArray());

You can also assign the data array to the matrix without actually copying it.

TMatrixD h;
h.Use(5,5,data.GetArray());
h.Invert();

The array data now contains the inverted matrix.

Now a unit matrix in sparse format is created.

TMatrixDSparse unit1(5,5);
TArrayI row(5),col(5);
for (Int_t i = 0; i < 5; i++) row[i] = col[i] = i;
TArrayD data(5);
data.Reset(1.);
unit1.SetMatrixArray(5,row.GetArray(),col.GetArray(),data.GetArray());

TMatrixDSparse unit2(5,5);
unit2.SetSparseIndex(5);
unit2.SetRowIndexArray(row.GetArray());
unit2.SetColIndexArray(col.GetArray());
unit2.SetMatrixArray(data.GetArray());

Inverting a matrix

  • Use the Invert(Double_t &det=0) function to invert a matrix:
TMatrixD a(...);
a.Invert();

– or –

  • Use the appropriate constructor to invert a matrix:
TMatrixD b(kInvert,a);

Both methods are available for general and symmetric matrices.

For matrices whose size is less than or equal to 6x6, the InvertFast(Double_t &det=0) function is available. Here the Cramer algorithm is used, which is faster but less accurate.

Using decomposition classes for inverting

You can also use the following decomposition classes (see Matrix decompositions) for inverting a matrix:

Name Matrix type Comment
TDecompLU General
TDecompQRH General
TDecompSVD General Can manipulate singular matrix.
TDecompBK symmetric
TDecompChol Symmetric Matrix should also be positive definite.
TDecompSparse Sparse

If the required matrix type is general, you also can handle symmetric matrices.

Example

This example shows how to check whether the matrix is singular before attempting to invert it.

TDecompLU lu(a);
TMatrixD b;
if (!lu.Decompose()) {
    cout << "Decomposition failed, matrix singular ?" << endl;
    cout << "condition number = " << = a.GetCondition() << endl;
} else {
    lu.Invert(b);
}

Matrix operators and methods

The matrix/vector operations are classified according to BLAS (basic linear algebra subroutines) levels.

Arithmetic operations between matrices

Description Format Comment
Element C=A+B Overwrites A
Wise sum A+=B
Add (A,alpha,B)
TMatrixD(A,TMatrixD::kPlus,B)		 A = A + α B constructor
Element wise subtraction C=A-B A-=B
TMatrixD(A,TMatrixD::kMinus,B)		 Overwrites A
Constructor
Matrix multiplication C=A*B
A*=B
C.Mult(A,B)
TMatrixD(A,TMatrixD::kMult,B)
TMatrixD(A, TMatrixD(A, TMatrixD::kTransposeMult,B)
TMatrixD(A, TMatrixD::kMultTranspose,B)		 Overwrites A
 
 
Constructor of A.B
Constructor of AT .B
Constructor of A.BT
Element wise multiplication ElementMult(A,B) A(i,j)*= B(i,j)
Element wise division ElementDiv(A,B) A(i,j)/= B(i,j)

Arithmetic operations between matrices and real numbers

Description Format Comment
Element wise sum C=r+A C=A+r A+=r overwrites A
Element wise subtraction C=r-A C=A-r A-=r overwrites A
Matrix multiplication C=r*A C=A*r A*=r overwrites A

Comparison between two matrices

Format Output Description
A == B Bool_t Equal to
A != B matrix Not equal
A > B matrix Greater than
A >= B matrix Greater than or equal to
A < B matrix Smaller than
A <= B matrix Smaller than or equal to
AreCompatible(A,B) Bool_t Compare matrix properties
Compare(A,B) Bool_t Return summary of comparison
VerifyMatrixIdentity(A,B,verb, maxDev)   Check matrix identity within maxDev tolerance

Comparison between matrix and real number

Format Output Description
A == r Bool_t Equal to
A != r Bool_t Not equal
A > r Bool_t Greater than
A >= r Bool_t Greater than or equal to
A < r Bool_t Smaller than
A <= r Bool_t Smaller than or equal to
VerifyMatrixValue(A,r,verb, maxDev) Bool_t Compare matrix value with r within maxDev tolerance
A.RowNorm() Double_t Norm induced by the infinity vector norm
A.NormInf() Double_t  
A.ColNorm() Double_t Norm induced by the 1 vector norm
A.Norm1() Double_t  
A.E2Norm() Double_t Square of the Euclidean norm
A.NonZeros() Int_t  
A.Sum() Double_t Number of elements unequal zero
A.Min() Double_t  
A.Max() Double_t  
A.NormByColumn (v,"D") TMatrixD  
A.NormByRow (v,"D") TMatrixD  

Matrix views

With the following matrix view classes, you can access the matrix elements:

The classes are templated, the usual specializations are available (e.g. D instead of T for double precision).

For the matrix view classes TMatrixDRow, TMatrixDColumn and TMatrixDDiag, the necessary assignment operators are available to interact with the vector class TVectorD. The sub matrix view classes TMatrixDSub has links to the matrix classes TMatrixD and TMatrixDSym.

The next table summarizes how to access the individual matrix elements in the matrix view classes.

Format Comment
TMatrixDRow(A,i)(j) TMatrixDRow(A,i)[j] Element Aij
TMatrixDColumn(A,j)(i) TMatrixDColumn(A,j)[i] Element Aij
TMatrixDDiag(A(i) TMatrixDDiag(A[i] Element Aij
TMatrixDSub(A(i) TMatrixDSub(A,rl,rh,cl,ch)(i,j) Element Aij
Element Arl+i,cl+j

Matrix decompositions

There are the following classes available for matrix decompositions:

  • TDecompLU: Decomposes a general n x n matrix A into P A = L U.
  • TDecompBK: The Bunch-Kaufman diagonal pivoting method decomposes a real symmetric matrix A.
  • TDecompChol : The Cholesky decomposition class, which decomposes a symmetric, positive definite matrix A = U^T * U where U is a upper triangular matrix.
  • TDecompQRH: QR decomposition class.
  • TDecompSVD: Single value decomposition class.
  • TDecompSparse: Sparse symmetric decomposition class.

Matrix Eigen analysis

With the TMatrixDEigen and TMatrixDSymEigen classes, you can compute eigenvalues and eigenvectors for general dense and symmetric real matrices.

The following table lists the methods of the TMatrixDEigen and the TMatrixDSymEigen to obtain the eigenvalues and eigenvectors. TMatrixDSymEigen constructors can only be called with TMatrixDSym:

Format Output Description
eig.GetEigenVectors() TMatrixD Eigenvectors for both TMatrixDEigen and TMatrixDSymEigen.
eig.GetEigenValues() TVectorD Eigenvalues vector for TMatrixDSymEigen.
eig.GetEigenValues() TMatrixD Eigenvalues matrix for TMatrixDEigen.
eig.GetEigenValuesRe() TVectorD Real part of eigenvalues for TMatrixDEigen.
eig.GetEigenValuesIm() TVectorD Imaginary part of eigenvalues for TMatrixDEigen.

Example

The usage of the eigenvalue class is shown in this example where it is checked that the square of the singular values of a matrix c are identical to the eigenvalues of cT.c:

const TMatrixD m = THilbertMatrixD(10,10);
TDecompSVD svd(m);
TVectorD sig = svd.GetSig(); sig.Sqr();

// Symmetric matrix EigenVector algorithm.
TMatrixDSym mtm(TMatrixDBase::kAtA,m);
const TMatrixDSymEigen eigen(mtm);
const TVectorD eigenVal = eigen.GetEigenValues();
const Bool_t ok = VerifyVectorIdentity(sig,eigenVal,1,1.-e-14);

The SMatrix package

SMatrix is a C++ package for high performance vector and matrix computations. It can be used only in problems when the size of the matrices is known at compile time, like in the tracking reconstruction of HEP experiments. It is based on a C++ technique, called expression templates, to achieve an high level optimization. The C++ templates can be used to implement vector and matrix expressions in such a way that these expressions can be transformed at compile time to code equivalent to hand-optimized code in a low-level language such as FORTRAN or C.

The SMatrix has been developed initially by T. Glebe of the Max-Planck-Institut, Heidelberg, as part of the HeraB analysis framework. A subset of the original package has been now incorporated in the ROOT distribution, with the aim to provide to the LHC experiments a stand-alone and high performance matrix package for reconstruction. The API of the current package differs from the original one to conform to ROOT coding conventions.

This package contains the two following generic classes for describing matrices and vectors of arbitrary dimensions and of arbitrary type:

SVector

The template class ROOT::Math::SVector represents n-dimensional vectors for objects of arbitrary type. The class has two template parameters that define their properties at compile time:

  1. Type of the contained elements (for example float or double).
  2. Size of the vector.

Use one of the following constructors to create a vector:

  • Default constructor for a zero vector (all elements equal to zero).
  • Constructor (and assignment) from a vector expression, like v=p*q+w. Due to the expression template technique, no temporary objects are created in this operation.
  • Constructor by passing directly the elements. This is possible only for vectors up to size 10.
  • Constructor from an iterator copying the data referred by the iterator. It is possible to specify the begin and end of the iterator or the begin and the size. Note that for

Example

The namespace ROOT::Math is used.

// create an empty vector of size 3 ( v[0]=v[1]=v[2]=0).
SVector<double,3> v;
double d[3] = {1,2,3};
// create a vector from a C array.
SVector<double,3> v(d,3);

SVector<double,3> v;
v[0] = 1;                 // set the first element
v(1) = 2;                 // set the second element
*(v.begin()+3) = 3;       // set the third element

// set vector elements from a std::vector<double>::iterator
std::vector<double> w(3);
v.SetElements(w.begin(), w.end()); // v now contains three zeros

// place a sub-vector in a vector.
SVector<double,N>  v;
SVector<double,M>  w;
int ioff = 2; // offset
// M <= N otherwise a compilation error is obtained later.
// place a vector of size M starting from element ioff, v[ioff+i]=w[i]
v.Place_at(w, ioff);
// return a sub-vector of size M starting from ioff
// v[ioff]: w[i]=v[ioff+i]
w = v.Sub<SVector<double>, M>(ioff);

SMatrix

The template class ROOT::Math::SMatrix represents a matrix of arbitrary type with nrows x ncoldimension. The class has four template parameters that define their properties at compile time:

  • type of the contained elements (for example float or double)
  • number of rows
  • number of columns
  • representation type
  • ROOT::Math::MatRepStd for a general nrows x ncols matrix. This class is itself a template on the contained type T, the number of rows and the number of columns. Its data member is an array T[nrows*ncols] containing the matrix data. The data are stored in the row-major C convention.
  • ROOT::Math::MatRepSym for a symmetric matrix of size NxN. This class is a template on the contained type and on the symmetric matrix size N. It has as data member an array of type T of size N*(N+1)/2, containing the lower diagonal block of the matrix. The order follows the lower diagonal block, still in a row-major convention.

Use one of the following constructors to create a matrix:

  • Default constructor for a zero matrix (all elements equal to zero).
  • Constructor of an identity matrix.
  • Copy constructor (and assignment) for a matrix with the same representation, or from a different one when possible (for example from a symmetric to a general matrix).
  • Constructor (and assignment) from a matrix expression, like D=A*B+C. Due to the expression template technique, no temporary objects are created in this operation. In the case of an operation like A=A*B+C, a temporary object is needed and it is created automatically to store the intermediary result in order to preserve the validity of this operation.
  • Constructor from a generic STL-like iterator copying the data referred by the iterator, following its order. It is both possible, to specify the begin and end of the iterator or the begin and the size. In case of a symmetric matrix, it is required only the triangular block and the user can specify whether giving a block representing the lower (default case) or the upper diagonal part.

Example Typedef’s are used in this example to avoid the full C++ names for the matrix classes. For a general matrix, the representation has the default value ROOT::Math::MatRepStd. For a general square matrix, the number of columns can be omitted.

// Typedef definitions used in the following declarations:
typedef ROOT::Math::SMatrix<double,3> SMatrix33;
typedef ROOT::Math::SMatrix<double,2> SMatrix22;
typedef ROOT::Math::SMatrix<double,3,3,
ROOT::Math::MatRepSym<double,3>> SMatrixSym3;
typedef ROOT::Math::SVector>double,2> SVector2;
typedef ROOT::Math::SVector>double,3> SVector3;
typedef ROOT::Math::SVector>double,6> SVector6;
SMatrix33 m0; // create a zero 3x3 matrix

// Create a 3x3 identity matrix.
SMatrix33 i = ROOT::Math::SMatrixIdentity();
double a[9] = {1,2,3,4,5,6,7,8,9}; // input matrix data

// Create a matrix using the a[] data.
// This results in the 3x3 matrix:
SMatrix33 m(a,9);

Example

A symmetric matrix is filled from a std::vector.

std::vector<double> v(6);
for (int i = 0; i<6; ++i) v[i] = double(i+1);

// This creates the symmetric matrix:
SMatrixSym3 s(v.begin(),v.end())

// Create a general matrix from a symmetric matrix (the opposite does not compile)
SMatrix33 m2 = s;

SMatrix and SVector objects can be printed using the Print method or the « operator:

// m is a SMatrix or a SVector object
m.Print(std::cout);
std::cout << m << std::endl;

Minimization libraries and classes

ROOT provides several minimization libraries and classes:

TMinuit

The Minuit minimization package was originally written in Fortran by Fred James and part of PACKLIB (patch D506). It has been converted to the C++ class TMinuit, by R.Brun.

Topical manual

For TMinuit, a topical manual it available at Topical Manual - TMinuit.
It contains in-depth information about TMinuit.

Minuit2 library

The Minuit2 library is a new object-oriented implementation, written in C++, of the popular MINUIT minimization package. These new version provides basically all the functionality present in the old Fortran version, with almost equivalent numerical accuracy and computational performances.

Furthermore, it contains new functionality, like the possibility to set single side parameter limits or the FUMILI algorithm (see FUMILI minimization package), which is an optimized method for least square and log likelihood minimizations. The package has been originally developed by M. Winkler and F. James.

Topical manuals

For Minuit2, topical manuals are available at Topical Manuals - Minuit2.
They contain in-depth information about Minuit2.

FUMILI minimization package

FUMILI is used to minimize Chi-square function or to search maximum of likelihood function.

FUMILI is based on ideas, proposed by I.N. Silin. It was converted from FORTRAN to C by Sergey Yaschenko s.yaschenko@fz-juelich.de.

For detailed information on the FUMILI minimization package, see TFumili class reference.

Numerical minimization

ROOT provides algorithms for one-dimensional und multi-dimensional numerical minimizations.

One-dimensional minimization

The one-dimensional minimization algorithms are used to find the minimum of a one-dimensional minimization function. The function to minimize must be given to the class implementing the algorithm as a ROOT::Math::IBaseFunctionOneDim object.

You can apply one-dimensional minimization in the following ways:

  • ROOT::Math::BrentMinimizer1D
  • ROOT::Math::GSLMInimizer1D
  • Using the TF1 class

ROOT::Math::BrentMinimizer1D

ROOT::Math::BrentMinimizer1D implements the Brent method to minimize an one-dimensional function. You must provide an interval containing the function minimum.

Example

In this example a function is defined to minimize as a lambda function. The function to minimize must be given to the class implementing the algorithm as a ROOT::Math::IBaseFunctionOneDim object.

ROOT::Math::Functor1D func( [](double x){ return 1 + -4*x + 1*x*x; } );

ROOT::Math::BrentMinimizer1D bm;
bm.SetFunction(func, -10,10);
bm.Minimize(10,0,0);
std::cout << "Minimum: f(" << bm.XMinimum() << ") = " << bm.FValMinimum() << std::endl;

Note that when setting the function to minimize, you must provide the interval range to find the minimum. In the `Minimize call, the maximum number of function calls, the relative and absolute tolerance must be provided.

ROOT::Math::GSLMInimizer1D

ROOT::Math::GSLMInimizer1D wraps two different methods from the GNU Scientific Library (GSL). At construction time you can choose between the BRENT and the GOLDENSECTION algorithmen. The GOLDENSECTION algorithm is the simplest method but the slowest and the BRENT algorithm (default) combines the golden section with a parabolic interpolation.
You can choose the algorithm as a different enumeration in the constructor:

  • * ROOT::Math::Minim1D::kBRENT for the BRENT algorithm (default).
  • * ROOT::Math::Minim1D::kGOLDENSECTION for the GOLDENSECTION algorithm.

Example

// This creates a class with the default BRENT algorithm.
ROOT::Math::GSLMinimizer1D minBrent;

// This creates a class with the GOLDENSECTION algorithm
ROOT::Math::GSLMinimizer1D minGold(ROOT::Math::Minim1D::kGOLDENSECTION);

Using the TF1 class

You can perform a one-dimensional minimization or maximization of a function by using directly the TF1 class. The minmization implemented in TF1 uses the BrentMInimizer1D and is available with the class member functions * TF1::GetMinimum or TF1::GetMaximum to find the function minimum (* TF1::GetMinimumX) or the maximum value (TF1::GetMaximumX). You can provide the search interval for the minimum, the tolerance and the maximum iterations as optional parameters of the * TF1::GetMinimum or TF1::GetMaximum functions.

Multi-dimensional minimization

The algorithms for a multi-dimensional minimization are implemented in the ROOT::Math::Minimizer interface. They can be used the same way as it was shown for the one-dimensional mimimization.

ROOT statistics classes

ROOT provides statistics classes for computing limits and confidence levels, fitting and multi-variate analysis.

Classes for computing limits and confidence levels

  • TFeldmanCousins: calculates the confidence levels of the upper or lower limit for a Poisson process using the Feldman-Cousins method (as described in PRD V57 #7, p3873-3889). No treatment is provided in this method for the uncertainties in the signal or the background.
  • TRolke: computes the confidence intervals for the rate of a Poisson process in the presence of background and efficiency, using the profile likelihood technique for treating the uncertainties in the efficiency and background estimate. The signal is always assumed to be Poisson; background may be Poisson, Gaussian, or user-supplied. efficiency may be Binomial, Gaussian, or user-supplied. See publication at Nucl. Instrum. Meth. A551:493-503,2005.
  • TLimit: computes 95% of the confidence level limits using the likelihood ratio semi-Bayesian method (method; see e.g., T. Junk, NIM A434, p. 435-443, 1999). It takes signal background and data histograms wrapped in a TLimitDataSource as input, and runs a set of Monte Carlo experiments in order to compute the limits. If needed, inputs are fluctuated according to systematic.

Specialized classes for fitting

  • TFractionFitter: fits Monte Carlo fractions to data histogram (à la HMCMLL, R. Barlow and C. Beeston, Comp. Phys. Comm. 77 (1993) 219-228). It accounts for the both data and the statistical Monte Carlo uncertainties through a likelihood fit using Poisson statistics. However, the template (Monte Carlo) predictions are also varied within statistics, leading to additional contributions to the overall likelihood. This leads to many more fitting parameters (one per bin per template), but minimization with respect to these additional parameters is performed analytically rather than introducing them as formal fitting parameters. Some special care needs to be taken in the case of bins with zero content.
  • TMultiDimFit: implements a multi-dimensional function parametrization for multi-dimensional data by fitting them to multi-dimensional data using polynomial or Chebyshev or Legendre polynomial.
  • RooFit: toolkit for fitting and data analysis modeling.

Multi-variate analysis classes

UNU.RAN

UNU.RAN (Universal Non Uniform RAndom Number generator for generating non-uniform pseudo-random numbers) contains universal (also called automatic or black-box) algorithms that can generate random numbers from large classes of distributions:

  • continuous (in one or multi-dimensions)
  • discrete distributions
  • empirical distributions (such as histograms)

UNU.RAN is an ANSI C library licensed under GPL.

The TUnuran class is used to interface the UNURAN package.

→ Unuran tutorials

Initializing TUnuran with string API

You can initialize UNU.RAN with the string API via TUnuran::Init(), passing the distribution type and the method.

Example

TUnuran unr;
// Initialize UNU.RAN to generate normal random numbers using an "arou" method.
unr.Init("normal()","method=arou");
...

// Sample distributions N times (generate N random numbers).
for (int i = 0; i<N; ++i) double x = unr.Sample();

Using TUnuranContDist for a one-dimensional distribution

Use TUnuranContDist for creating a continuous 1-D distribution object (for example from a TF1 object providing the probability density function). You can provide additional information via TUnuranContDist::SetDomain(min,max) like the domain() for generating numbers in a restricted region.

Example

// 1D case: create a distribution from two TF1 object, pointers pdfFunc.
TUnuranContDist dist(pdfFunc);

// Initialize UNU.RAN, passing the distribution and a string.
// Define the method.
unr.Init(dist, "method=hinv");

// Sample distribution N times (generate N random numbers).
for (int i = 0; i < N; ++i) double x = unr.Sample();

Using TUnuranMultiContDist for a multi-dimensional distribution

Use TUnuranMultiContDist to create a multi-dimensional distribution that can be created from a multi-dimensional PDF (probability density function).

Example

// Multi- dimensional case from TF1 (TF2 or TF3) objects.
TUnuranMultiContDist dist(pdfFuncMulti);

// The recommended method for a multi-dimensional function is "hitro".
unr.Init(dist,"method=hitro");

// Sample distribution N times (generate N random numbers).
double x[NDIM];
for (int i = 0; i<N; ++i) unr.SampleMulti(x);

Using TUnuranDiscrDist for a discrete one-dimensional distribution

Use TUnuranDiscrDist to create a discrete one-dimensional distribution that can be initialized from a TF1 object or from a vector of probabilities.

Example

// Create a distribution from a vector of probabilities.
double pv[NSize] = {0.1,0.2,...};
TUnuranDiscrDist dist(pv,pv+NSize);

// The recommended method for a discrete distribution is "dgt".
unr.Init(dist, "method=dgt");

// Sample N times (generate N random numbers).
for (int i = 0; i < N; ++i) int k = unr.SampleDiscr();

Using TUnuranEmpDist for an empirical distribution

Use TUnuranEmpDist for creating an empirical distribution that can be initialized from a TH1 object (using the bins or from its buffer for unbinned data) or from a vector of data.

Example

// Create a distribution from a set of data.
// vdata is an std::vector containing the data.
TUnuranEmpDist dist(vdata.begin(),vdata.end());
unr.Init(dist);

// Sample N times (generate N random numbers).
for (int i = 0; i<N; ++i) double x = unr.Sample();

FOAM

FOAM is a simplified version of a multi-dimensional general purpose Monte Carlo event generator (integrator) with hyper-cubical “foam of cells”.

→ FOAM tutorials

Certain features of full version of FOAM are omitted. mFOAM is intended as an easy to use tool for Monte Carlo simulation and integration in few dimensions. It relies on the ROOT package, borrowing persistency of classes from ROOT. You can use mFOAM from the ROOT interpreter.

Examples

foam_kanwa.C is a simple example on running FOAM in interactive mode.

foam_demo.C shows the usage of FOAM in compiled mode, which is the preferred method.

foam_demopers.C demonstrates the persistency of FOAM classes.

FFTW

For computing fast Fourier transforms (FFT), ROOT uses the FFTW library. To use it, the fftw3 CMake module must be enabled.

With SetDefaultFFT() you can change the default library.

TVirtualFFT is the interface class for FFT. With TH1::FFT() you can perform a FFT for a histogram.

→ fft tutorials

Quadp

Quadp is an optimization library with linear and quadratic programming methods. It is based on the matrix package.

Example

An example of using Quadp can be found in the portfolio.C tutorial.