diff --git a/documentation/users-guide/FittingHistograms.md b/documentation/users-guide/FittingHistograms.md
index 6b39d45faed25fa293f2248c99a83e6e90127ad8..989636ba454ac33ab18a4f16c9a3299837536856 100644
--- a/documentation/users-guide/FittingHistograms.md
+++ b/documentation/users-guide/FittingHistograms.md
@@ -9,33 +9,25 @@ needs to be drawn in a pad before the Fit Panel is invoked. The method
## The Fit Method
+The Fit method is implemented in ROOT for the histogram classes `TH1`,
+the sparce histogram classes, `THnSparse`, the graph classes, `TGraph`,
+`TGraph2D` and `TMultiGraph` for fitting a collection of Graphs with the same function.
+
+### The TH1::Fit Method
To fit a histogram programmatically, you can use the `TH1::Fit`
-method. Here is the signature of `TH1::Fit` and an explanation of the
+method. Here is the signatures of `TH1::Fit` and an explanation of the
parameters:
``` {.cpp}
- void Fit(const char *fname, Option_t *option, Option_t *goption,
+ TFitResultPtr Fit(TF1 *function, Option_t *option, Option_t *goption,
Axis_t xxmin, Axis_t xxmax)
```
+- `function` a pointer to the fitted function (the fit model) object.
+ One can also use the function name. This name may be one of ROOT pre-defined
+ function names or a user-defined function. See the next paragraph for the list of pre-defined functions.
-- `*fname: `The name of the fitted function (the model) is passed as
- the first parameter. This name may be one of ROOT pre-defined
- function names or a user-defined function. The functions below are
- predefined, and can be used with the `TH1::Fit` method:
-
-- "`gaus`" Gaussian function with 3 parameters:
- `f(x) = p0*exp(-0.5*((x-p1)/p2)^2)`
-
-- "`expo`"An Exponential with 2 parameters: `f(x) = exp(p0+p1*x)`
-
-- "`pol`*`N`*" A polynomial of degree *N*:
- `f(x) = p0 + p1*x + p2*x2 +...`
-
-- "`landau`" Landau function with mean and sigma. This function has
- been adaptedfrom the `CERNLIB` routine `G110 denlan`.
-
-- `*option:`The second parameter is the fitting option. Here is the
+- `*option:` The second parameter is the fitting option. Here is the
list of fitting options:
- "`W`" Set all weights to 1 for non empty bins; ignore error bars
@@ -46,17 +38,21 @@ parameters:
- "`I`" Use integral of function in bin instead of value at bin
center
-- "`L`" Use log likelihood method (default is chi-square method)
+- "`L`" Use log likelihood method (default is chi-square method). To be used when
+the histogram represents counts
-- "`U`" Use a user specified fitting algorithm
+- "`WL`" Weighted log likelihood method. To be used when the histogram has been filled with
+ weights different than 1.
- "`Q`" Quiet mode (minimum printing)
- "`V`" Verbose mode (default is between Q and V)
+- "`S`" The result of the fit is returned in the `TFitResultPtr`.
+
- "`E`" Perform better errors estimation using the Minos technique
-- "`M`" Improve fit results
+- "`M`" Improve fit results, by using the *IMPROVE* algorithm of TMinuit.
- "`R`" Use the range specified in the function range
@@ -70,17 +66,15 @@ parameters:
one is kept)
- "`B`"Use this option when you want to fix one or more parameters
- and the fitting function is like `polN`, `expo`, `landau`, `gaus`.
+and the fitting function is a predefined one, like `polN`, `expo`, `landau`, `gaus`.
+Note that in case of pre-defined functions some default initial values and limits are set.
-- "`LL`"An improved Log Likelihood fit in case of very low
- statistics and when bincontentsare not integers. Do not use this
- option if bin contents are large (greater than 100).
- "`C`"In case of linear fitting, don't calculate the chisquare
(saves time).
-- "`F`"If fitting a `polN`, switch to `Minuit` fitter (by default,
- `polN` functions are fitted by the linear fitter).
+- "`F`"If fitting a linear function (e.g. `polN`), switch to use the default minimizer (e.g. `Minuit`). By default,
+ `polN` functions are fitted by the linear fitter.
- `*goption: `The third parameter is the graphics option that is the
same as in the **`TH1`**`::Draw` (see the chapter Draw Options).
@@ -88,10 +82,29 @@ parameters:
- `xxmin`, `xxmax:`Thee fourth and fifth parameters specify the
range over which to apply the fit.
-By default, the fitting function object is added to the histogram and
+By default, the fitted function object is added to the histogram and
is drawn in the current pad.
-## Fit with a Predefined Function
+
+### The TGraph::Fit Method
+
+The signature for fitting a TGraph is exactly the same as for the `TH1`. Only some options apply only for fitting histograms,
+these are the options `L`, `WL` and `I`.
+These options apply instead only for `TGraph::Fit`:
+
+* `TGraph` specific *options*
+
+- "`EX0`" When fitting a `TGraphErrors` or a `TgraphAsymErrors` the errors on the coordinates are not used in the fit
+
+- "`ROB`" in case of linear fitting use the Robust fitting. Compute the LTS regression coefficients (robust (resistant) regression),
+ using the default fraction of good points.
+- "`ROB=0.x`" as above, but compute the LTS regression coefficients, using 0.x as a fraction of good points.
+
+## The `TF1` function class
+
+Here we will show how to create the `TF1` class that is used for fitting histograms and graphs.
+
+### Fit with a Predefined Function
To fit a histogram with a predefined function, simply pass the name of
@@ -102,10 +115,30 @@ this line fits histogram object `hist` with a Gaussian.
root[] hist.Fit("gaus");
```
-The initial parameter values for pre-defined functions are set
-automatically.
+The initial parameter values (and eventual limits) for pre-defined functions are set
+automatically. For overriding the default limits values use the fit option `B`.
+
+The list of pre-defined functions that can be used with the `Fit` method is the following:
-## Fit with a User-Defined Function
+- "`gaus`" Gaussian function with 3 parameters:
+ `f(x) = p0*exp(-0.5*((x-p1)/p2)^2)`
+
+- "`expo`"An Exponential with 2 parameters: `f(x) = exp(p0+p1*x)`
+
+- "`pol`*`N`*" A polynomial of degree *N*, where N is a number between 0 and 9:
+ `f(x) = p0 + p1*x + p2*x2 +...`
+
+- "`chebyshev`*`N`*" A Chebyshev polynomial of degree *N*, where N is a number between 0 and 9:
+ `f(x) = p0 + p1*x + p2*(2*x2-1) +...`
+
+- "`landau`" Landau function with mean and sigma. This function has
+been adapted from the `CERNLIB` routine `G110 denlan` (see `TMath::Landau`).
+
+- "`gausn` Normalized form of the gaussian function with 3 parameters
+`f(x) = p0*exp(-0.5*((x-p1)/p2)^2)/(p2 *sqrt(2PI))`
+
+
+### Creating User-Defined Functions (TF1)
You can create a **`TF1`** object and use it in the call the
@@ -117,9 +150,10 @@ the **`TF1`** object. There are three ways to create a **`TF1`**.
- Same as first one, with parameters
-- Using a function that you have defined
+- Using a function that you have defined. This can be a free function or
+ a functor object or a particular member function of a class.
-### Creating a TF1 with a Formula
+#### Creating a TF1 with a Formula
Let's look at the first case. Here we call the **`TF1`** constructor
@@ -136,7 +170,7 @@ You can also use a **`TF1`** object in the constructor of another
root[] TF1 *f2 = new TF1("f2","f1*2",0,10)
```
-### Creating a TF1 with Parameters
+#### Creating a TF1 with Parameters
The second way to construct a **`TF1`** is to add parameters to the
@@ -169,11 +203,11 @@ This sets parameter 0 to 10 and parameter 1 to 5. We can now draw the
root[] f1->Draw()
```
-### Creating a TF1 with a User Function
+#### Creating a TF1 with a User Function
The third way to build a **`TF1`** is to define a function yourself
-and then give its name to the constructor. A function for a **`TF1`**
+and then pass the function pointer to the constructor. A function for a **`TF1`**
constructor needs to have this exact signature:
``` {.cpp}
@@ -232,8 +266,48 @@ Now we use the function:
}
```
-## Fixing and Setting Parameters' Bounds
+You can create a TF1 also from a C++ function object (functor) with parameters
+A TF1 can be created from any C++ class implementing this member function:
+
+```{.cpp}
+double operator()(double *x, double *p)
+```
+The advantage of the function object is that it can have a state and reference therefore what-ever other object
+the user needs, without using globals. This is an example to define first the function object
+
+``` {.cpp}
+class MyFunctionObject {
+ public:
+ // use constructor to customize your function object
+ MyFunctionObject(......) { ......}
+
+ double operator() (double *x, double *p) {
+ // function implementation using class data members
+ }
+};
+```
+
+and then use it to create the `TF1`:
+
+```{.cpp}
+ MyFunctionObject fobj(....); // create the function object
+ TF1 * f = new TF1("f",fobj,xmin,xmax,npar); // create TF1 class with n-parameters and range [xmin,xmax]
+```
+
+If using C++11, one can create a `TF1` also from a C++ `lambda` function:
+
+```{.cpp}
+// create TF1 class with 2 parameters and range [xmin,xmax] using a lambda
+TF1 * f = new TF1("f",[](double*x,double*p){return p[0] + p[1]*x[0];},xmin,xmax,2);
+```
+
+## Configuring the Fit
+
+We will show here some configuration actions that can or must be done
+when fitting histogram or graph using the `Fit` method.
+
+### Fixing and Setting Parameters' Bounds
Parameters must be initialized before invoking the `Fit` method. The
setting of the parameter initial values is automatic for the
@@ -278,8 +352,7 @@ With this setup, parameters 0`->`2 can vary freely, parameter 3 has
boundaries [-10, 4] with initial value -1.5, and parameter 4 is fixed
to 0.
-## Fitting Sub Ranges
-
+### Fitting Sub Ranges
By default, `TH1::Fit` will fit the function on the defined histogram
range. You can specify the option "`R`" in the second parameter of
@@ -301,140 +374,8 @@ root[] hist->Fit("f1","","",-2,2)
See macros `$ROOTSYS/tutorials/fit/myfit.C` and `multifit.C` as more
completed examples.
-## The Fit Panel
-
-
-
-
-To display the Fit Panel right click on a histogram to pop up the
-context menu, and then select the menu entry Fit Panel.
-
-The new Fit Panel GUI is available in ROOT v5.14. Its goal is to
-replace the old Fit Panel and to provide more user friendly way for
-performing, exploring and comparing fits.
-
-By design, this user interface is planned to contain two tabs:
-"General" and "Minimization". Currently, the "General" tab provides
-user interface elements for setting the fit function, fit method and
-different fit, draw, print options.
-
-The new fit panel is a modeless dialog, i.e. when opened, it does not
-prevent users from interacting with other windows. Its first prototype
-is a singleton application. When the Fit Panel is activated, users can
-select an object for fitting in the usual way, i.e. by left-mouse
-click on it. If the selected object is suitable for fitting, the fit
-panel is connected with this object and users can perform fits by
-setting different parameters and options.
-
-### Function Choice and Settings
-
-*‘Predefined' combo box* - contains a list of predefined functions in
-ROOT. You have a choice of several polynomials, a Gaussian, a Landau,
-and an Exponential function. The default one is Gaussian.
-
-*‘Operation' radio button group* defines the selected operational mode
-between functions:
-
-*Nop* - no operation (default);
-
-*Add* - addition;
-
-*Conv* - convolution (will be implemented in the future).
-
-Users can enter the function expression into the text entry field
-below the ‘Predefined' combo box. The entered string is checked after
-the Enter key was pressed and an error message shows up, if the
-function string is not accepted.
-
-‘*Set Parameters*' button opens a dialog for parameters settings,
-which will be explaned later.
-
-### Fitter Settings
-
-
-*‘Method' combo box* currently provides only two fit model choices:
-Chi-square and Binned Likelihood. The default one is Chi-square. The
-Binned Likelihood is recomended for bins with low statistics.
-
-*‘Linear Fit' check button* sets the use of Linear fitter when is
-selected. Otherwise the minimization is done by Minuit, i.e. fit
-option "`F`" is applied. The Linear fitter can be selected only for
-functions linears in parameters (for example - `polN)`.
-
-*‘Robust' number entry* sets the robust value when fitting graphs.
-
-*‘No Chi-square' check button* switch On/Off the fit option "`C`" -
-do not calculate Chi-square (for Linear fitter).
-
-*‘Integral' check button* switch On/Off the option "`I`" - use
-integral of function instead of value in bin center.
-
-*‘Best Errors'* sets On/Off the option "`E`" - better errors
-estimation by using Minos technique.
-
-*‘All weights = 1'* sets On/Off the option "`W`"- all weights set to 1
-excluding empty bins; error bars ignored.
-
-*‘Empty bins, weights=1'* sets On/Off the option "`WW`" - all weights
-equal to 1 including empty bins; error bars ignored.
-
-*‘Use range'* sets On/Off the option "`R`" - fit only data within the
-specified function range. Sliders settings are used if this option is
-set to On. Users can change the function range values by pressing the
-left mouse button near to the left/right slider edges. It is possible
-to change both values simultaneously by pressing the left mouse button
-near to the slider center and moving it to a new position.
-
-*‘Improve fit results'* sets On/Off the option "`M`"- after minimum is
-found, search for a new one.
-
-*‘Add to list'* sets On/Off the option "`+`"- add function to the list
-without deleting the previous one. When fitting a histogram, the
-function is attached to the histogram's list of functions. By default,
-the previously fitted function is deleted and replaced with the most
-recent one, so the list only contains one function. Setting this
-option to On will add the newly fitted function to the existing list
-of functions for the histogram. Note that the fitted functions are
-saved with the histogram when it is written to a ROOT file. By
-default, the function is drawn on the pad displaying the histogram.
-
-### Draw Options
-
-
-*‘SAME'* sets On/Off function drawing on the same pad. When a fit is
-executed, the image of the function is drawn on the current pad.
-
-*‘No drawing'* sets On/Off the option "`0`"- do not draw the fit
-results.
-
-*‘Do not store/draw'* sets On/Off option "`N`"- do not store the
-function and do not draw it.
-
-### Print Options
-
-
-This set of options specifies the amount of feedback printed on the
-root command line after performed fits.
-
-*‘Verbose'* - prints fit results after each iteration.
-
-*‘Quiet'* - no fit information is printed.
-
-*‘Default'* - between Verbose and Quiet.
-
-### Command Buttons
-
-
-*Fit button* - performs a fit taking different option settings via the
-Fit Panel interface.
-
-*Reset* - sets the GUI elements and related fit settings to the
-default ones.
-
-*Close* - closes the Fit panel window.
-
-## Fitting Multiple Sub Ranges
+### Fitting Multiple Sub Ranges
The script for this example is `$ROOTSYS/tutorials/fit/multifit.C`. It
@@ -493,7 +434,7 @@ the "+" sign is explained below:
h->Fit(total,"R+");
```
-## Adding Functions to the List
+### Adding Functions to the List
The example `$ROOTSYS/tutorials/fit/multifit.C` also illustrates how
@@ -510,8 +451,8 @@ root[] hist->Fit("f1","+","",-2,2)
Note that the fitted function(s) are saved with the histogram when it
is written to a ROOT file.
-## Combining Functions
+## Example of fit: Combining Functions
You can combine functions to fit a histogram with their sum as it is
illustrated in the macro `FitDemo.C`
@@ -592,7 +533,12 @@ function:
-## Associated Function
+## Result of the fit
+
+Here we will show how to obtain the result of the fit (fitted function, parameter values, errors
+and eventually the covariance and correlation matrix).
+
+### Associated Function
One or more objects (typically a **`TF1`**\*) can be added to the list
@@ -604,7 +550,7 @@ of functions (`fFunctions`) associated to each histogram. A call to
TF1 *myfunc = h->GetFunction("myfunc");
```
-## Access to the Fit Parameters and Results
+### Access to the Fit Parameters and Results
If the histogram (or graph) is made persistent, the list of associated
@@ -621,7 +567,10 @@ root[] Double_t p1 = fit->GetParameter(0);
root[] Double_t e1 = fit->GetParError(0);
```
-## Associated Errors
+Using the fit option `S` one can access the full result of the fit including the covariance and correlation matrix.
+See later the paragraph `TFitResult`.
+
+### Associated Errors
By default, for each bin, the sum of weights is computed at fill time.
@@ -636,10 +585,13 @@ weights)`; otherwise, the error is set equal to the
```
Empty bins are excluded in the fit when using the Chi-square fit method.
-When fitting the histogram with the low statistics, it is recommended to
-use the Log-Likelihood method (option ‘`L`' or "`LL`").
+When fitting an histogram representing counts (i.e with Poisson statistics) it is recommended to
+use the Log-Likelihood method (option ‘`L`' or "`WL`"), particularly in case of low statistics.
+When the histogram has been filled with weights different than one, a weighted likelihood method can be used
+and the errors retrieved from the fit are corrected following a procedure described in paragraph 8.5.2 of the book,
+*F. James, Statistical Methods in Experimental Physics, 2nd Edition*.
-## Fit Statistics
+### Fit Statistics
You can change the statistics box to display the fit parameters with the
@@ -658,7 +610,595 @@ errors, use:
gStyle->SetOptFit(1011);
```
-## The Minimization Package
+## The Fit Panel
+
+
+
+
+To display the Fit Panel right click on a histogram to pop up the
+context menu, and then select the menu entry Fit Panel.
+
+The new Fit Panel GUI is available in ROOT v5.14. Its goal is to
+replace the old Fit Panel and to provide more user friendly way for
+performing, exploring and comparing fits.
+
+By design, this user interface is planned to contain two tabs:
+"General" and "Minimization". Currently, the "General" tab provides
+user interface elements for setting the fit function, fit method and
+different fit, draw, print options.
+The "Minimization tab" provides the option to set the Minimizer to use in the fit and
+its specific options.
+
+The new fit panel is a modeless dialog, i.e. when opened, it does not
+prevent users from interacting with other windows. Its first prototype
+is a singleton application. When the Fit Panel is activated, users can
+select an object for fitting in the usual way, i.e. by left-mouse
+click on it. If the selected object is suitable for fitting, the fit
+panel is connected with this object and users can perform fits by
+setting different parameters and options.
+
+### Function Choice and Settings
+
+
+*‘Predefined' combo box* - contains a list of predefined functions in
+ROOT. You have a choice of several polynomials, a Gaussian, a Landau,
+and an Exponential function. The default one is Gaussian.
+
+*‘Operation' radio button group* defines the selected operational mode
+between functions:
+
+*Nop* - no operation (default);
+
+*Add* - addition;
+
+*Conv* - convolution (will be implemented in the future).
+
+Users can enter the function expression into the text entry field
+below the ‘Predefined' combo box. The entered string is checked after
+the Enter key was pressed and an error message shows up, if the
+function string is not accepted.
+
+‘*Set Parameters*' button opens a dialog for parameters settings,
+which will be explaned later.
+
+### Fitter Settings
+
+
+*‘Method' combo box* currently provides only two fit model choices:
+Chi-square and Binned Likelihood. The default one is Chi-square. The
+Binned Likelihood is recomended for bins with low statistics.
+
+*‘Linear Fit' check button* sets the use of Linear fitter when is
+selected. Otherwise the minimization is done by Minuit, i.e. fit
+option "`F`" is applied. The Linear fitter can be selected only for
+functions linears in parameters (for example - `polN)`.
+
+*‘Robust' number entry* sets the robust value when fitting graphs.
+
+*‘No Chi-square' check button* switch On/Off the fit option "`C`" -
+do not calculate Chi-square (for Linear fitter).
+
+*‘Integral' check button* switch On/Off the option "`I`" - use
+integral of function instead of value in bin center.
+
+*‘Best Errors'* sets On/Off the option "`E`" - better errors
+estimation by using Minos technique.
+
+*‘All weights = 1'* sets On/Off the option "`W`"- all weights set to 1
+excluding empty bins; error bars ignored.
+
+*‘Empty bins, weights=1'* sets On/Off the option "`WW`" - all weights
+equal to 1 including empty bins; error bars ignored.
+
+*‘Use range'* sets On/Off the option "`R`" - fit only data within the
+specified function range. Sliders settings are used if this option is
+set to On. Users can change the function range values by pressing the
+left mouse button near to the left/right slider edges. It is possible
+to change both values simultaneously by pressing the left mouse button
+near to the slider center and moving it to a new position.
+
+*‘Improve fit results'* sets On/Off the option "`M`"- after minimum is
+found, search for a new one.
+
+*‘Add to list'* sets On/Off the option "`+`"- add function to the list
+without deleting the previous one. When fitting a histogram, the
+function is attached to the histogram's list of functions. By default,
+the previously fitted function is deleted and replaced with the most
+recent one, so the list only contains one function. Setting this
+option to On will add the newly fitted function to the existing list
+of functions for the histogram. Note that the fitted functions are
+saved with the histogram when it is written to a ROOT file. By
+default, the function is drawn on the pad displaying the histogram.
+
+### Draw Options
+
+
+*‘SAME'* sets On/Off function drawing on the same pad. When a fit is
+executed, the image of the function is drawn on the current pad.
+
+*‘No drawing'* sets On/Off the option "`0`"- do not draw the fit
+results.
+
+*‘Do not store/draw'* sets On/Off option "`N`"- do not store the
+function and do not draw it.
+
+### Advances Options
+
+The advance option button is enabled only after having performed the fit and provides
+additional drawing options that can be used after having done the fit. These new drawing tools,
+which can be selected by the "Advanced Drawing Tool" panel that pops up when clicking the "Advanced" button, are:
+
+* *Contour*: to plot the confidence contour of two chosen parameters. One can select the number of points to draw the contour
+(more points might require more time to compute it), the parameters and the desired confidence level .
+
+* *Scan* : to plot a scan of the minimization function (likelihood or chi-squared) around the minimum as function of the chosen parameter.
+
+* *Conf Interval* : to plot the confidence interval of the fitted function as a filled coloured band around its central value.
+ One can select the desired confidence level for the band to be plotted.
+
+### Print Options
+
+
+This set of options specifies the amount of feedback printed on the
+root command line after performed fits.
+
+*‘Verbose'* - prints fit results after each iteration.
+
+*‘Quiet'* - no fit information is printed.
+
+*‘Default'* - between Verbose and Quiet.
+
+### Command Buttons
+
+
+*Fit button* - performs a fit taking different option settings via the
+Fit Panel interface.
+
+*Reset* - sets the GUI elements and related fit settings to the
+default ones.
+
+*Close* - closes the Fit panel window.
+
+### Minimization Options
+
+With this tab one can select specific options for minimization. These include
+
+* The minimizer library ( *Minuit*, *Minuit2*, *Fumili*, *GSL*, *Genetics* )
+* The method (algorithm) for minimization. For example for Minuit one can choose between (*Migrad*, *Simplex* or *Scan*)
+* Error definition
+* Minimization tolerance
+* Number of iterations/function calls
+* Print Level: (*Default*, *Verbose* or *Quiet*).
+
+
+
+## New ROOT::Fit classes
+
+The fitting of the data objects in ROOT, histograms, graphs and tree is performed via some common classes,
+which are defined in the `ROOT::Fit` namespace.
+These classes can be classified in the following groups:
+
+* User classes driving the fit: `ROOT::Fit::Fitter` for executing the fit, `ROOT::Fit::FitConfig` for configuring the fit,
+ `ROOT::Fit::ParameterSettings` to define the properties of the fit parameters (initial
+ values, bounds, etc..), `ROOT::Fit::FitResult` for storing the result of the fit.
+* Data classes containing the data sets used in the fitting. These classes are the`ROOT::Fit::BinData`for describing bin data sets,
+ thus data points containing both coordinates and a corresponding value/weight
+ with optionally an error on the value or the coordinate and the `ROOT::Fit::UnBinData` for un-binned data sets,
+ which consisst only of a vector of coordinate values. The coordinate values can be
+ one-dimensional (i.e. one entry per event) or multi-dimensional (N entries per event).
+* Function classes defining the type of fit (the objective function used for fitting):
+ - `ROOT::Fit::Chi2FCN` for chi2 (least-square fits),
+ - `ROOT::Fit::PoissonLikelihoodFCN` for binned likelihood fits of histograms,
+ - `ROOT::Fit::LogLikelihoodFCN` for generic un-binned likelihood fits.
+ These classes are templated on the type of function interface they implement (see later). User convenient typedefs are also provided.
+ They derive from the common generic interface multi-dimensional for function evaluation, `ROOT::Math::IBaseFunctionMultiDim`.
+
+In addition the fitter classes make uses of the generic interfaces for parametric function evaluations, `ROOT::Math::IParametricFunctionMultiDim`
+to define the fit model function and use the `ROOT::Math::Minimizer` interface to perform the minimization of the objective function.
+More information about the function interface and the multi-dimensional minimization in ROOT is given in the Mathematical Library chapter.
+
+Here we present a detailed description of the `ROOT::Fit` classes and how to use them.
+Using these classes instead of the interface provided directly in the ROOT data objects, like `TH1::Fit` allow are more fine control
+to configure and customise the fits. For example, using these classes a combined fit of several histograms can be performed.
+
+To understand how these class work, let's go through a simple example, such as fitting an histogram.
+
+When fitting an histogram, instead of using `TH1::Fit` we will show in the following hot wo use the `ROOT::Fit` classes.
+We will show how to perform the following different type of fits with the histogram data:
+* a least square fit using the observed errors (Neyman chi-squared);
+* a least square fit using the expected errors from the function (Pearson chi-squared);
+* a binned likelihood fit;
+* an extended unbinned likelihood fits, if the histogram has been set to store in the buffer the original data used to fill it.
+
+
+Let's go through all the steps required for performing these fits using the `ROOT::Fit::Fitter` class.
+These steps are:
+1. Create the input fit data object.
+2. Create the input model function.
+3. Configure the fit.
+4. Perform the data fitting.
+5. Examine the result.
+
+### Creating the input fit data
+
+We have two types of input data, binned data (class `ROOT::Fit::BinData`) used for least square (chi-square) fits of histograms or `TGraph` objects
+or un-binned data (class `ROOT::Fit::UnBinData`) used for
+fitting vectors of data points (e.g. from a `TTree`).
+
+#### Using Binned data
+
+Let's suppose we have an histogram, represented as a `TH1` type object (it can be one or multi-dimensional). The following shows how to create and
+fill a `ROOT:Fit::BinData` object.
+
+``` {.cpp}
+ ROOT::Fit::DataOptions opt;
+ opt.fIntegral = true;
+ ROOT::Fit::BinData data(opt);
+ // fill the bin data using the histogram
+ // we can do this using the following helper function from the Hist library
+ TH1 * h1 = (TH1*) gDirectory->Get("myHistogram");
+ ROOT::Fit::FillData(data, h1);
+```
+
+In this code example, we have used the utility function of the *Hist* library, `ROOT::Fit::FillData` to fill the `BinData` object.
+The `ROOT::Fit::FillData` is defined in the headerfile `HFitInterface.h` and it has a signature for all different ROOT objects,
+like `TH1`, `THnBase`, `TGraph`, `TGraph2D` and `TMultiGraph`
+It is possible to specify, when creating the `BinData` object, the data range we want to use and some fitting options we want to apply
+to fill in the object and later when fitting.
+The fit data options are controlled by the ``ROOT::Fit::DataOptions`` class, the range by the ``ROOT::Fit::DataRange`` class.
+
+Here is an example how to specify the input option to use the integral of the function value in the bin instead of using the function value
+evaluated at the bin center, when doing the fit and to use a
+range beween the 'xmin' and 'xmax' values.
+
+``` {.cpp}
+ ROOT::Fit::DataOptions opt;
+ opt.fIntegral = true;
+ ROOT::Fit::DataRange range(xmin,xmax);
+ ROOT::Fit::BinData data(opt,range);
+ // fill the bin data using the histogram
+ // we can do this using the following helper function from the Hist library
+ TH1 * h1 = (TH1*) gDirectory->Get("myHistogram");
+ ROOT::Fit::FillData(data, h1);
+```
+
+The list of possible fit options available is the following:
+``` {.cpp}
+ ROOT::Fit::DataOptions opt;
+ opt.fIntegral = true; // use integral of bin content instead of bin center (default is false).
+ opt.fBinVolume = true; // normalize data by the bin volume (default is false).
+ // This is for fitting density functions in histograms with variable bin sizes.
+ opt.fUseRange =true; // use the function range when creating the fit data (default is false).
+ opt.fExpErrors = true; // use the expected errors estimated from the function values
+ // assuming Poisson statistics and not the observed errors (default is false).
+ opt.fUseEmpty = true; // use empty bins when fitting (default is false). If fExpErrors
+ // is not set an arbitrary error = 1 is assigned to those bins.
+ opt.fErrors1 = true; // Set all measured errors to 1 (default is false).
+ opt.fCoordErrors = false; // When available coordinate errors are not used in the fit
+ // (default is true: the errors are used when they are available,
+ // e.g. fitting a TGraphErrors).
+ opt.fAsymErrors = false; // When available asymmetric errors are considered in the fit
+ // (default is true, the asymmetric errors are used when they are available,
+ // e.g. fitting a TGraphAsymmErrors).
+```
+
+The `ROOT::Fit::DataRange` class supports defining multiple rectangular ranges in each dimension, and supports n-dimension.
+The function `DataRange::AddRange(icoord,xmin,xmax)` adds a range in
+the coordinate `icoord` with lower value `xmin` and upper value `xmax`:
+
+
+``` {.cpp}
+ ROOT::Fit::DataRange range;
+ range.AddRange(icoordinate, xmin, xmax);
+```
+
+#### Using Un-Binned data
+
+The unbinned data sets are represented with the `ROOT::Fit::UnBinData` class.
+For creating un-binned data sets, a `ROOT::Fit::UnBinData` object, one has two possibilities:
+1. Copy the data inside `ROOT::Fit::UnBinData`. One can create an empty `UnBinData` object, iterate on the data and add the data point one by one, or directly create the `UnBinData`
+ object from a data iterator. In this case an input `ROOT::Fit::DataRange` object is passed in order to copy the data according to the given range.
+ 2. Use `ROOT::Fit::UnBinData` as a wrapper to an external data storage. In this case the `UnBinData` object is created from an iterator or pointers to the data and the data are not copied
+ inside. In this case the data cannot be selected according to a specified range. All the data points will be included in the fit.
+
+The `ROOT::Fit::UnBinData` class supports also weighted data. In addition to the data points (coordinates), which can be of arbitrary `k` dimensions, the class can be constructed from a vector of
+weights. This is an example of taking data from an histogram buffer of a `TH1` object:
+
+``` {.cpp}
+ double * buffer = histogram->GetBuffer();
+ // number of entry is first entry in the buffer
+ int n = buffer[0];
+ // when creating the data object it is important to create with the size of the data
+ ROOT::Fit::UnBinData data(n);
+ for (int i = 0; i < n; ++i)
+ data.add(buffer[2*i+1]); // the buffer of 1D histogram contains nevt,x1,w1,x2,w2,......
+```
+
+Instead in this example we will create a 2-dim `UnBinData` object with the contents from a ROOT `TTree`
+
+``` {.cpp}
+ TFile * file = TFile::Open("hsimple.root");
+ TTree *ntuple = 0; file->GetObject("ntuple",ntuple);
+ // select from the tree the data we want to use for fitting
+ // we use TTree::Draw for this
+ int nevt = ntuple->Draw("px:py","","goff");
+ double * x = ntuple->GetV1();
+ double * y = ntuple->GetV2();
+ ROOT::Fit::UnBinData data(nevt, x, y );
+```
+
+### Creating the Fit model
+
+In order to fit a data sets we need a model to describe our data, e.g. a probability density function describing our observed data or
+an hypothetical function describing the relation between the independent variables **`X`** and the single dependent variable `Y`.
+We can have an arbitrary number `k` of independent variables. For example, when fitting a `k`-dimensional histogram,
+the independent variables **`X`** are the bin center coordinates and `Y` is the bin weight.
+
+The model function needs to be expressed as function of some unknown parameters. The fitting will find the best parameter value to describe
+the observed data.
+
+We can use the ROOT `TF1` class, the parametric function class, to describe the model function. However the `ROOT::Fit::Fitter` class, to be independent of the ROOT *`Hist`* library,
+takes as input a more general parametric function object, the interface (abstract) class `ROOT::Math::IParametricFunctionMultiDim`, which describe a generic one or multi-dimensional function
+with parameters. This interface extends the abstract class `ROOT::Math::IBaseFunctionMultiDim`, with methods to set/retrieve parameter values and to evaluate the function given the
+independent vector of values **`X`** and vector of parameters `P`.
+More information about the different `ROOT::Math` function interfaces is available in the Mathematical Library chapter.
+
+An end-user can convert a `TF1` object in a `ROOT::Math::IParametricFunctionMultiDim`, using the wrapper class `ROOT::Math::WrapperMultiTF1`:
+
+``` {.cpp}
+ TF1 * f1 = new TF1("f1","gaus");
+ ROOT::Math::WrappedMultiTF1 fitFunction(f1, f1->GetNdim() );
+ ROOT::Fit::Fitter fitter;
+ fitter.SetFunction( fitFunction, false);
+```
+
+When creating the wrapper, the parameter values stored in `TF1` will be copied in the `ROOT::Math::WrappedMultiTF1` object.
+The function object representing the model function is given to the `ROOT::Fitter` class using the `Fitter::SetFunction` method.
+
+The user has also the possibility to provide a function object, which implements the derivatives of the function with respect
+to the parameters.
+This information might be useful for some types of fits. In this case he needs to provide the function object as a class deriving from the
+`ROOT::Math::IParametricGradFunctionMultiDim` interface.
+Note that the wrapper class `ROOT::Math::WrappedMultiTF1` implements also the gradient interface, using internally `TF1::GradientPar`,
+which is based on numerical differentiation, apart for the case of linear functions (i.e. when `TF1::IsLinear()` is `true`).
+The parameter derivatives of the model function can be useful to some minimization algorithms, such as Fumili.
+However, in general is better to leave the minimization algorithm (e.g. Minuit) to compute the needed derivatives using its own customised
+numerical differentiation algorithm.
+In order to not provide to the fitter the parameter derivatives, we explicitly passed in `Fitter::SetFunction` a `false` value.
+
+### Fit Configuration
+
+The configuration of the fit is done via the `ROOT::Fit::FitConfig` class and its contained `ROOT::Fit::ParameterSettings` class.
+These are the possible allowed fit configurations:
+
+ - setting the initial values of the parameters;
+ - setting the parameter step sizes;
+ - setting eventual parameter bounds;
+ - setting the minimizer library and the particular algorithm to use;
+ - setting different minimization options (print level, tolerance, max iterations, etc...)
+ - setting the type of parameter errors to compute (parabolic error, Minos errors, re-normalize errors using fitted chi2 values)
+
+
+The initial parameter values can be set directly in the input model function object.
+However, for setting parameter bounds and step sizes to values different than the automatically computed ones, one needs to use the `ROOT::Fit::ParameterSetting` class.
+This example code will set the lower/upper bounds for the first parameter and a lower bound for the second parameter
+
+``` {.cpp}
+ fitter.SetFunction( fitFunction, false);
+ fitter.Config().ParSettings(0).SetLimits(0,1.E6);
+ fitter.Config().ParSettings(2).SetLowerLimit(0);
+```
+
+Note that a `ROOT::Fit::ParameterSettings` objects exists for each fit parameter and it created by the `ROOT::Fit::FitConfig` class, after the model function has been set in the Fitter.
+Only when the function is set, the number of parameter is known and
+automatically the `FitConfig` creates the corresponding `ParameterSetting` objects.
+
+When fitting, different minimizer can be used. The can be implemented in different libraries and loaded ar run time by the plug-in manager system of ROOT.
+Each different minimizer (e.g. *Minuit, Minuit2, Fumili,* etc.) consists of a different implementation of the `ROOT::Math::Minimizer` interface.
+Within the same minimizer, thus within the same class implementing the `Minimizer` interface, different algorithms can exist.
+For example in the case of Minuit, we have *Migrad, Simplex* or *Minimize*. The minimizer and its corresponding algorithm, when available,
+can be set by using the function `FitConfig::SetMinimizer("minimizerName")` or by using directly the `ROOT:Math::MinimizerOptions` class.
+
+If the requested minimizer is not available in ROOT, the default one is used. The default minimizer type and algorithm can be specified by using the
+static function `ROOT::Math::MinimizerOptions::SetDefaultMinimizer("minimizerName")`
+
+
+### Minimizer Libraries and Algorithms
+
+The list of available minimizer libraries currently available in ROOT, with their corresponding available algorithms is the following one.
+Some minimizers (e.g. *Minuit*) contain several algorithms that the user can
+choose. Others are based on a single algorithm (e.g. *Fumili*)
+
+
+ - **`Minuit`** (library *libMinuit*). Old version of Minuit, based on the `TMinuit` class. The list of possible algorithms are:
+ - *`Migrad`* (default one)
+ - *`Simplex`*
+ - *`Minimize`* (it is a combination of Migrad and Simplex)
+ - *`MigradImproved`*
+ - *`Scan`*
+ - *`Seek`*
+
+
+ - **`Minuit2`** (library *libMinuit2*). New C++ version of Minuit. The list of possible algorithm is :
+ - *`Migrad`* (default)
+ - *`Simplex`*
+ - *`Minimize`*
+ - *`Scan`*
+ - *`Fumili`* . This is the same algorithm of `TFumili`, but implemented in the Minuit2 library.
+
+ - **`Fumili`**. Implement a dedicated minimization algorithm for least-square and likelihood fits. It has requirements on the type of method function to be used.
+ No specific algorithm exists
+
+ - **`GSLMultiMin`** (library *libMathMore*). Minimizer based on the Multidimensional Minimization routines of the Gnu Scientific Library (GSL). The list of available algorithms is
+ - *`BFGS2`* (default) : second version of the vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm;
+ - *`BFGS`* : old version of the vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm;
+ - *`ConjugateFR`* : Fletcher-Reeves conjugate gradient algorithm;
+ - *`ConjugatePR`* : Polak-Ribiere conjugate gradient algorithm;
+ - *`SteepestDescent`*: steepest descent algorithm;
+
+ - **`GSLMultiFit`** (library *libMathMore*). Minimizer based on the Non-Linear Least-Square routines of GSL. This minimizer can be used only for least-square fits.
+
+ - **`GSLSimAn`** (library *libMathMore*). Minimizer based on simulated annealing.
+
+ - **`Genetic`** (library *libGenetic*). Genetic minimizer based on an algorithm implemented in the *TMVA* package.
+
+Each minimizer can be configured using the `ROOT::Math::MinimizerOptions` class. The list of possible option that can be set are:
+* Minimizer type (`MinimizerOptions::SetMinimizerType(const char *)`) .
+* Minimizer algorithm (`MinimizerOptions::SetMinimizerAlgorithm(const char *)`).
+* Print Level (`MinimizerOptions::SetPrintLevel(int )`) to set the verbose printing level (default is 0).
+* Tolerance (`MinimizerOptions::SetTolerance(double )`) tolerance used to control the iterations.
+* Maximum number of function calls (`MinimizerOptions::SetMaxFunctionCalls(int )`).
+* Maximum number of iterations (`MinimizerOptions::SetMaxIterations(int )`). Note that this is not used by *Minuit*
+* FCN Upper value for Error Definition (`MinimizerOptions::SetMaxIterations(int )`). Value in the minimization function used to compute the parameter errors.
+ The default is to get the uncertainties at the 68% CL is a value of 1 for a chi-squared function minimization and 0.5 for a log-likelihood function.
+ * Strategy (`MinimizerOptions::SetStrategy(int )`), minimization strategy used. For each minimization strategy *Minuit* uses different configuration parameters
+ (e.g. different requirements in computing derivatives, computing full Hessian (strategy = 2) or an approximate version. The default is a value of 1. In this case the full Hessian matrix
+ is computed only after the minimization.
+* Precision (`MinimizerOptions::SetTolerance(double )`). Precision value in the evaluation of the minimization function. Default is numerical double precision.
+
+Note that not all the options are implemented by all the minimizers.
+For example in *Minuit* is possible to set the maximum number of function calls, but not the maximum number of iterations. The Strategy and the Precision options apply instead only for *Minuit* (and
+*Minuit2*).
+
+The class supports alo setting different default values for the options, by using the static functions `MinimizerOptions::SetDefault...` (e.g. `MinimizerOptions::SetDefaultPrintLevel(int )`).
+The list of the current option values can be inspected by using `MinimizerOptions::Print`.
+```{.cpp}
+ROOT::Math::MinimizerOptions() opt;
+// print the default minimizer option values
+opt.Print();
+```
+
+In addition it is possible to provide extra options which might apply for a particular minimizer `MinimizerOptions::SetExtraOptions(const IOptions & )`.
+See the documentation of the particular minimizer to use for the list of possible additional options available.
+
+
+### Performing the Fit
+
+Here we have now all the required input ingredients for the fit, the data and the function to fit.
+Depending on these we have now several different way to perform the fit, using the corresponding methods of the
+`ROOT::Fit::Fitter` class and depending on the type of input data.
+
+#### Available fit methods
+
+* **Least-square fit**: `Fitter::LeastSquare(const BinData & )` or `Fitter::Fit(const Bindata &)`. It requires the user to pass a `BinData` object. It should be used when the data values follow a
+ Gaussian distribution. This fit method is implemented using the class `ROOT::Fit::Chi2FCN`.
+* **Binned Likelihood fit** : `Fitter::LikelihoodFit(const Bindata & )`. The user needs to pass a `BinData` object. It should be used when the data values follow a Poisson or a multinomial
+distribution. The Poisson case (extended fit) is the default and in this case the function normalization is also fit to the data. The Multi-nominal case can be selected by passing the optional
+*extended* boolean flag as *false*. This method is implemented by the class `ROOT::Fit:::PoissonLikelihoodFCN`.
+* **Un-Binned likelihood fit**: `Fitter::LikelihoodFit(const UnBindata &)`. The user needs to pass an `UnBinData` object. By default the fit is not extended (i.e. the normalization is not fitted to the
+data). As above the user can select an extended likelihood fit by passing the optional
+*extended* boolean flag as *true*. This methos is implemented using the class `LogLikelihoodFCN`
+* **Linear Fit**: A linear fit can be selected (no iterative minimization is needed in this case, but using linear algebra algorithms from the *Matrix* library), if the model function is linear in the
+ parameters.
+
+#### Customised Fit methods
+
+Above we described the pre-defined methods used for fitting. A user can also implement its own fitting methods, thus its version of the chi-square or likelihood function he wants to minimize.
+In this cas, the user does not really need to build as input a `ROOT::Fit` data set and model function as we described before. He can implements its own version of the method function using on its own
+data set objects and functions.
+
+In this case `ROOT::Fit::Fitter::SetFCN` is used to set the method function and `ROOT::Fit::FitFCN` is used for fitting. The method function can be passed also in `ROOT::Fit::FitFCN`, but in this
+case a previously defined fitting configuration is used.
+
+The possible type of method functions that can be bassed in `ROOT::Fit::Fitter::SetFCN` are:
+
+* A generic functor object implementing `operator()(const double * p)` where **`p`** is the parameter vectors. In this case one needs to pass the number of parameters,
+ the function object and optionally a vector of initial parameter values. Other optional parameter include the size of the data sets and a flag specifying if it is a chi2 (least-square fit).
+ In the last two parameters are given, the `chi2/ndf` can be computed after fitting the data.
+``` {.cpp}
+template <class Function>
+bool Fitter::SetFCN(unsigned int npar, Function & f,
+ const double * initialParameters = 0,
+ unsigned int dataSize=0, bool isChi2Fit = false)
+```
+
+* A function object implementing the `ROOT::Math::IBaseFunctionMultiDim` interface:
+``` {.cpp}
+bool Fitter::SetFCN(const ROOT::Math::IBaseFunctionMultiDim & f,
+ const double * initialParameters = 0,
+ unsigned int dataSize=0, bool isChi2Fit = false)
+```
+
+* A function object implementing the `ROOT::Math::FitMethodFunction` interface. This is an interface class extending
+the `ROOT::Math::IBaseFunctionMultiDim` with some extra functionality which can be used when fitting.
+This extra functionality is required by dedicated fitting algorithms like *Fumili* or *GSLMultiFit*.
+``` {.cpp}
+bool Fitter::SetFCN(const ROOT::Math::FitMethodFunction & f,
+ const double * initialParameters = 0, unsigned int dataSize=0)
+```
+
+* A old-Minuit like FCN interface (i.e. a free function with the signature `fcn(int &npar, double *gin, double &f, double *u, int flag)`.
+``` {.cpp}
+typedef void(* MinuitFCN)(int &npar, double *gin, double &f, double *u, int flag)
+bool Fitter::SetFCN(MinuitFCN fcn, int npar,
+ const double * initialParameters = 0,
+ unsigned int dataSize=0, bool isChi2Fit = false)
+```
+
+### Fit Result
+
+The result of the fit is contained in the `ROOT::Fit::Result` object. A reference to the result object is obtained with the function
+`Fitter::Result()`.
+The `ROOT::Fit::FitResult` class provides an API for retrieving parameter values, errors, covariance and correlation matrix from the fit,
+minimum chi2/likelihood values, etc...
+
+A `FitResult::Print` method is also available to print the result of the fit.
+
+The class has a self-explanatory API so, see its reference documentation for the possible information available after the fit.
+
+One extra functionality offered by `ROOT::Fit::FitResult` is the possibility to compute the confidence intervals of the function after the fit.
+The function `ROOT::Fit::FitResult::GetConfidenceInterval` given an input data sets (e.g. a `BinData` object) and a confidence level value (e.g. 68%)
+computes the lower/upper band values of the model function at the given data points.
+
+### TFitResult
+
+`TFitResult` is a class deriving from `ROOT::Fit::Result` and providing in addition some convenient methods to return a
+covariance or correlation matrix as a `TMatrixDSym` object. In addition `TFitResult` derives from `TNamed` and can be conveniently
+stored in a file.
+
+When fitting an histogram ( a `TH1` object) or a graph (a `TGraph` object) it is possible to return a `TFitResult` via the `TFitResultPtr` object,
+which behaves as a smart pointer to a `TFitResult`.
+`TFitResultPtr` is the return object by `TH1::Fit` or `TGraph::Fit`.
+By default the TFitResultPtr contains only the status of the fit and can be obtained by an automatic conversion of the TFitResultPtr to an integer.
+
+If the fit option *`S`* is instead used, `TFitResultPtr` contains the `TFitResult` and behaves as a smart
+pointer to it. This is an example:
+
+``` {.cpp}
+int fitStatus = hist->Fit(myFunction); // TFitResultPtr contains only the fit status
+
+TFitResultPtr r = hist->Fit(myFunction,"S"); // TFitResultPtr contains the TFitResult
+TMatrixDSym cov = r->GetCovarianceMatrix(); // to access the covariance matrix
+Double_t chi2 = r->Chi2(); // to retrieve the fit chi2
+Double_t par0 = r->Parameter(0); // retrieve the value for the parameter 0
+Double_t err0 = r->ParError(0); // retrieve the error for the parameter 0
+r->Print("V"); // print full information of fit including covariance matrix
+r->Write(); // store the result in a file
+```
+
+## The Minimization packages
+
+As explained before various minimization packages can be used when fitting in ROOT.
+We have seen before how to configure the `Fitter` class to use different minimization packages
+and different minimization options.
+When using the `Fit` method the minimization package (and its options) can be selected using the
+static methods of the `ROOT::Math::MinimizerOptions` class.
+For example to select `Minuit2` instead of `Minuit` for fitting an histogram do:
+
+``` {.cpp}
+ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
+// fit the histogram histo with the gaussian pre-defined function
+histo->Fit("gaus");
+```
+
+In the following we will give some brief description of the minimization packages.
+The packages all implement the `ROOT::Math::Minimizer` interface which can be use for
+finding the minimum of a multi-dimensional function.
+The interface is documented in the Mathematical Library Chapter.
+
+In addition packages like Minuit or Minuit2 provide their own interfaces.
+
+## MINUIT (Old TMInuit Version)
This package was originally written in FORTRAN by Fred James and part
@@ -691,8 +1231,6 @@ conversion of the original FORTRAN version. The main changes are:
- The ROOT static function `Printf` is provided to replace all
format statements and to print on currently defined output file
-- The derived class **`TMinuitOld`** contains obsolete routines from
- the FORTRAN based version
- The functions `SetObjectFit/GetObjectFit` can be used inside the
`FCN` function to set/get a referenced object instead of using
@@ -939,6 +1477,52 @@ following:
- Starting too far from the solution - the function may have
unphysical local minima, especially at infinity in some variables.
+
+
+## Minuit2 Package
+
+
+`Minuit2` is a new object-oriented implementation, written in C++, of
+the popular `MINUIT` minimization package. Compared with the
+**`TMinuit`** class, which is a direct conversion from FORTRAN to C++,
+`Minuit2` is a complete redesign and re-implementation of the package.
+This new version provides all the functionality present in the old
+FORTRAN version, with almost equivalent numerical accuracy and
+computational performances.
+Furthermore, it contains some fixes and small improvements and this new functionality:
+* The possibility to set single side parameter limits
+* the FUMILI algorithm (see the next paragraph "FUMILI Minimization Package"),
+which is an optimized method for least square and log
+likelihood minimizations.
+
+Minuit2 has been originally developed by M.
+Winkler and F. James in the SEAL project. More information can be found
+on the [MINUIT Web Site](MINUIT Web Site) and in particular at the
+following documentation page at
+<http://www.cern.ch/minuit/doc/doc.html>.
+
+A detailed User Guide for Minuit2 exists, describing the API of the internal classes.
+ROOT uses `Minuit2` for fitting via the `Minuit2Minimizer` class which implements
+the `ROOT::Math::Minimizer` interface.
+
+`Minuit2` is also distributed as an independent package of ROOT and can be built
+without any other dependency on the ROOT libraries.
+
+Examples on how to use the `Minuit2` and `Fumili2` plug-ins are provided
+in the tutorials' directory `$ROOTSYS/tutorials/fit`:
+`minuit2FitBench.C`, `minuit2FitBench2D.C` and `minuit2GausFit.C`.
+More information on the classes and functions present in `Minuit2` is
+available at
+[online reference documentation](online reference documentation).
+
+Useful information on MINUIT and minimization in general is provided in the
+following documents:
+
+F. James, *Minuit Tutorial on Function Minimization* (
+<http://seal.cern.ch/documents/minuit/mntutorial.pdf>); F. James, *The
+Interpretation of Errors in Minuit* (
+<http://seal.cern.ch/documents/minuit/mnerror.pdf>);
+
## FUMILI Minimization Package
@@ -1031,6 +1615,8 @@ similar step formulae are used in FUMILI for negative logarithm of
the likelihood function with the same idea - linearization of function
argument.
+
+
## Neural Networks
diff --git a/documentation/users-guide/MathLibraries.md b/documentation/users-guide/MathLibraries.md
index fc94d6471c4dae4753bbd59f43a7ad7caaac0d2e..3859e6d58f87461d3e15757fcde5223173436f14 100644
--- a/documentation/users-guide/MathLibraries.md
+++ b/documentation/users-guide/MathLibraries.md
@@ -13,30 +13,252 @@ and vector with fixed sizes and their operation has been developed
+
+## MathCore Library
+
+`MathCore` provides a collection of functions and C++ classes for
+numerical computing. This library includes only the basic mathematical
+functions and algorithms and not all the functionality required by the
+physics community. A more advanced mathematical functionality is
+provided by the `MathMore` library. The current set of included classes,
+which are provided in the `ROOT::Math` namespace are:
+
+
+- Basic special functions like the gamma, beta and error function.
+
+- Mathematical functions used in statistics, such as the probability
+ density functions and the cumulative distributions functions (lower
+ and upper integral of the pdf's).
+
+- Generic function classes and interfaces
+ for evaluating one-dimensional (`ROOT::Math::IBaseFunctiononeDim`) and multi-dimensional functions
+ (`ROOT::Math::IBaseFunctionMultiDim`) and parametric function interfaces for evaluating functions with parameters in one
+ (`ROOT::Math::IParametricFunctionOneDim`) or multi dimensions (`ROOT::Math::IParametricFunctionMultiDim`).
+ A set of user convenient wrapper classes, such as `ROOT::Math::Functor` is provided for wrapping user-classes in the needed interface,
+ required to use the algorithms of the `ROOT` Mathematical libraries.
+
+- Numerical algorithms interfaces and in same cases default implementations for:
+ - numerical integration;
+ - numerical differentiation;
+ - one dimensional root-finding;
+ - one-dimensional minimization;
+ - multi-dimensional minimization (only the `ROOT::Math::Minimizer` interface)
+
+- Fitting classes: set of classes for fitting generic data sets. These classes are provided in the namespace `ROOT::Fit`.
+ They are describing separately in the Fitting chapter.
+
+The sets described above is independednt of ROOT libraries and can be built as a set of standalone classes.
+In addition `MathCore` provides the following classes (depending on ROOT *libCore* library):
+
+- `TMath`: namespace with mathematical functions and basic function algorithms.
+- `TComplex`: class for complex numbers.
+- Random classes: the base class `TRandom` and the derived classes `TRandom1`, `TRandom2` and `TRandom3`, implementing the pseudo-random number generators.
+
+A detailed description for all `MathCore` classes is available in the Doxygen
+[online reference documentation](online reference documentation).
+
+## 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 current set, provided in
+the `ROOT::Math` namespace
+include:
+
+- Special mathematical functions (like Bessel functions, Legendre polynomials, etc.. )
+
+- Additional mathematical functions used in statistics such as probability density
+functions, cumulative distributions functions and their inverse which are not in `MathCore` but present in
+the `GSL` library.
+
+- Numerical algorithms for one dimensional functions based on
+implementation of the GNU Scientific Library (GSL):
+
+- Numerical integration classes implementing the interface **`ROOT::Math::Integrator`**
+ which is based on the Adaptive integration algorithms of QUADPACK
+
+- Numerical differentiation via **`ROOT::Math::GSLDerivator`**
+
+- Root finder implementing the **`ROOT::Math::RootFinder`** interface, using different
+ solver algorithms from GSL
+
+- one-dimensional Minimization implementing the interface**`ROOT::Math::IMinimizer1D`**
+
+- Interpolation via **`ROOT::Math::Interpolation`**. All the GSL
+ interpolation types are supported
+
+- Function approximation based on Chebyshev polynomials via the class
+ **`ROOT::Math::Chebyshev`**
+
+- Random number generators and distributions based on GSL using the `ROOT::Math::Random<Engine_type>` class.
+
+- Polynomial evaluation and root solvers
+
+The mathematical functions are implemented as a set of free functions in
+the namespace **`ROOT::Math`**. The naming used for the special
+functions is the same proposed for the C++ standard (see C++ standard
+extension [proposal document](proposal document)).The `MathMore` library
+is implemented wrapping in C++ the GNU Scientific Library ( <GSL>).
+Building `MathMore` requires a version of GSL larger or equal 1.8. The
+source code of `MathMore` is distributed under the GNU General Public
+License.
+
+`MathMore` (and its ROOT Cling dictionary) can be built within ROOT
+whenever a GSL library is found in the system. The GSL library and
+header file location can be specified in the ROOT configure script, by
+doing:
+
+```
+./configure --with-gsl-incdir=... --with-gsl-libdir=...
+```
+
+`MathMore` can be built also a stand-alone library (without requiring
+ROOT) downloding the tar file from the Web at this link. In this case
+the library will not contain the dictionary information and therefore
+cannot be used interactively
+
+More information on the classes and functions present in `MathMore` is
+available in the
+[online reference documentation](online reference documentation).
+
+
## TMath
-In the namespace, **`TMath`** a collection of free functions is provided
-for the following functionality:
+In the namespace, **`TMath`**, a collection of free functions is provided for the following functionality:
- numerical constants (like `pi`, `e`, `h`, etc.);
-- elementary and trigonometric functions;
+- trigonometric and elementary mathematical functions;
-- functions to find `min` and `max` of arrays;
+- functions to work with arrays and collections (e.g. functions to find `min` and `max` of arrays);
-- statistic functions to find mean and `rms` of arrays of data;
+- statistic functions to work on array of data (e.g. mean and `RMS` of arrays);
- algorithms for binary search/hashing sorting;
- special mathematical functions like `Bessel`, `Erf`, `Gamma`, etc.;
- statistical functions, like common probability and cumulative
- (quantile) distributions
+(quantile) distributions
+
+- geometrical functions.
For more details, see the reference documentation of **`TMath`** at
<http://root.cern.ch/root/htmldoc/TMath.html>.
+
+### Numerical Constants
+
+`TMath` offers a wide range of constants in the form of inline functions. Notice that they are not defined as C/C++ preprocessor macros. This set of functions includes one or more definitions for the following constants:
+
+* Pi.
+* Base of natural logarithm.
+* Velocity of light.
+* Gravitational constant (G).
+* Standard acceleration of gravity (g).
+* Standard acceleration of Gravity.
+* Plank's contant.
+* Boltzmann's and Steffan-Boltzmann's constants.
+* Avogadro's number.
+* Universal gas constant.
+* Molecular weight of dry air.
+* Dry air gas constant.
+* Euler-Mascheroni Constant.
+* Elementary charge.
+
+### Elementary Functions
+
+A set of miscellaneous elementary mathematical functions is provided along with a set of basic trigonometrical functions. Some of this functions refer to basic mathematical functions like 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 (like `Factorial`), most of the functionality offered here is just a wrapper of the first ones. Nevertheless, some of the them also offer some security checks or a better precision, like the trigonometrical functions `ASin(x)`, `ACos(x)` or `ATan(x)`.
+
+```{.cpp}
+ // Generate a vector with 10 random numbers
+ vector<double> v(10);
+ std::generate(v.begin(), v.end(), rand);
+
+ // Find the minumum 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;
+ ```
+
+Another example of these functions can be found in $ROOTSYS/tutorials/permute.C.
+
+### Statistic Functions Operating on Arrays.
+
+This set of functions processes arrays to calculate:
+
+* Mean.
+* Median.
+* Geometrical mean.
+* Sample Standard Deviation (*RMS*).
+* The kth smallest element.
+
+These functions, as the array algorithms, have two different interfaces. An old-style one where the size of the array is passed as a first argument followed by a pointer to the array itself
+and a modern C++-like interface that receives two iterators to it.
+
+```{.cpp}
+ // 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());
+
+ // old-style
+ mean = TMath::Mean(n, &v[0]);
+
+ // Calculate the mean with a weight vector
+ // with iterators
+ mean = TMath::Mean(v.begin(), v.end(), w.begin());
+
+ // old-style
+ mean = TMath::Mean(n, &v[0], &w[0]);
+ ```
+
+### Special and Statistical Functions.
+
+`TMath` also provides special functions like Bessel, Error functions, Gamma or similar plus statistical mathematical functions, including probability density functions, cumulative distribution and their inverse.
+
+The majority of the special functions and the statitical distributions are provided also as free functions in the `ROOT::Math` namespace.
+See one of the next paragraph for the complete description of the functions provided in `ROOT::Math`.
+The user is encourage to use those versions of the algorithms rather than the ones in TMath.
+
+Functions not present in `ROOT::Math` and provided only by `TMath` are:
+
+* Special functions:
+ * DiLogarithm
+ * Struve
+
+* Statistical functions:
+ * KolmogorovProb
+ * Voigt function
+ * LaplaceDist
+ * Vavilov
+
+The example tutorial `GammaFun.C` and `mathBeta.C` in `$ROOTSYS/tutorials` shows an example of use of the `ROOT::Math` special functions
+
+
+
## Random Numbers
@@ -370,1291 +592,2029 @@ Here are the CPU times obtained using the four random classes on an
| `UNURAN` | | | | |
+--------------------+---------------+----------------+----------------+----------------+
-## MathCore Library
-
-
-`MathCore` provides a collection of functions and C++ classes for
-numerical computing. This library includes only the basic mathematical
-functions and algorithms and not all the functionality required by the
-physics community. A more advanced mathematical functionality is
-provided by the `MathMore` library. The current set included classes
-are:
-- Basic special functions like the gamma, beta and error function.
+## Mathematical Functions
-- Mathematical functions used in statistics, such as the probability
- density functions and the cumulative distributions functions (lower
- and upper integral of the pdf's).
-- `GenVector`: physics and geometry vectors for 3 and 4 dimensions
- with their transformations (rotations and boost).
-
-- Generic (`ROOT::Math::IFunction`) and parametric
- (**`ROOT::Math::IParamFunction`**) function interfaces for one and
- multi dimensions.
-
-A detailed description for all `MathCore` classes is available in the
-[online reference documentation](online reference documentation). The
-`MathCore` library presented in the ROOT distribution contains the Cling
-dictionary for I/O and interactive usage. For the template classes, the
-dictionary is provided for some of the possible types, such as those
-based on double and Double32\_t. For the I/O or interactive use of other
-types, the dictionary must be first generated. An example on how to
-generate the required dictionary is provided in the tutorial
-`mathcoreVectorFloatIO.C` (in `$ROOTSYS/tutorials/math`). `MathCore` can
-also be built as an independent package using `configure/make`. In this
-case the library will not contain the dictionary information and cannot
-be used interactively in ROOT.
+The mathematical functions are present in both `MathCore` and `MathMore`
+libraries. All mathematical functions are implemented as free functions
+in the namespace **`ROOT::Math`**. The most used functions are in the
+`MathCore` library while the others are in the `MathMore` library. The
+functions in `MathMore` are all using the implementation of the GNU
+Scientific Library (GSL). The naming of the special functions is the
+same defined in the C++
+[Technical Report on Standard Library extensions](Technical Report on
+Standard Library extensions).
+The special functions are defined in the header file `Math/SpecFunc.h`.
-## Generic Vectors for 2, 3 and 4 Dimensions (GenVector)
+### Special Functions in MathCore
-`GenVector` is a package intended to represent vectors and their
-operations and transformations, such as rotations and Lorentz
-transformations, in 3 and 4 dimensions. The 3D space is used to describe
-the geometry vectors and points, while the 4D space-time is used for
-physics vectors representing relativistic particles. These 3D and 4D
-vectors are different from vectors of the linear algebra package, which
-describe generic N-dimensional vectors. Similar functionality is
-currently provided by the CLHEP <Vector> and <Geometry> packages and the
-ROOT Physics vector classes (See "Physics Vectors"). It also re-uses
-concepts and ideas from the CMS
-[Common Vector package](Common Vector package). In contrast to CLHEP or
-the ROOT physics libraries, `GenVector` provides class templates for
-modeling the vectors. The user can control how the vector is internally
-represented. This is expressed by a choice of coordinate system, which
-is supplied as a template parameter when the vector is constructed.
-Furthermore, each coordinate system is itself a template, so that the
-user can specify the underlying scalar type.
+- `ROOT::Math::beta(double x,double y) - `evaluates the beta function:
+ $$B(x,y) = \frac{\Gamma(x) \Gamma(y)}{\Gamma(x+y)}$$
-The `GenVector` classes do not inherit from **`TObject`**, therefore
-cannot be used as in the case of the physics vector classes in ROOT
-collections.
+- `double ROOT::Math::erf(double x)` - evaluates the error function
+ encountered in integrating the normal
+ distribution:
+ $$erf(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt$$
-In addition, to optimize performances, no virtual destructors are
-provided. In the following paragraphs, the main characteristics of
-`GenVector` are described. A more detailed description of all the
-`GenVector` classes is available also at
-<http://seal.cern.ch/documents/mathlib/GenVector.pdf>
+- `double ROOT::Math::erfc(double x)` - evaluates the complementary
+ error function:
+ $$erfc(x) = 1 - erf(x) = \frac{2}{\sqrt{\pi}} \int_{x}^{\infty} e^{-t^2} dt$$
-### Main Characteristics
+- `double ROOT::Math::tgamma(double x)` - calculates the gamma
+ function:
+ $$\Gamma(x) = \int_{0}^{\infty} t^{x-1} e^{-t} dt$$
+### Special Functions in MathMore
-#### Optimal Runtime Performances
-We try to minimize any overhead in the run-time performance. We have
-deliberately avoided the use of any virtual function and even virtual
-destructors in the classes. In addition, as much as possible functions
-are defined as inline. For this reason, we have chosen to use template
-classes to implement the `GenVector` concepts instead of abstract or
-base classes and virtual functions. It is then recommended to avoid
-using the `GenVector` classes polymorphically and developing classes
-inheriting from them.
+- `double ROOT::Math::assoc_legendre(unsigned l,unsigned m,double x) -`computes
+ the associated Legendre polynomials (with `m>=0`, `l>=m` and
+ `|x|<1)`:
+ $$P_{l}^{m}(x) = (1-x^2)^{m/2} \frac{d^m}{dx^m} P_{l}(x)$$
-#### Points and Vector Concept
+- `double ROOT::Math::comp_ellint_1(double k)` - calculates the
+ complete elliptic integral of the first kind (with $0 \le k^2 \le 1$:
+ $$
+ K(k) = F(k, \pi / 2) = \int_{0}^{\pi /2} \frac{d \theta}{\sqrt{1 - k^2 \sin^2{\theta}}}
+ $$
-Mathematically vectors and points are two distinct concepts. They have
-different transformations, as vectors only rotate while points rotate
-and translate. You can add two vectors but not two points and the
-difference between two points is a vector. We then distinguish for the 3
-dimensional case, between points and vectors, modeling them with
-different classes:
+- `double ROOT::Math::comp_ellint_2(double k)` - calculates the
+ complete elliptic integral of the second kind (with $0 \le k^2 \le 1$):
+ $$
+ E(k) = E(k , \pi / 2) = \int_{0}^{\pi /2} \sqrt{1 - k^2 \sin^2{\theta}} d \theta
+ $$
-- `ROOT::Math::`**`DisplacementVector2D`** and
- `ROOT::Math::`**`DisplacementVector3D`** template classes describing
- 2 and 3 component direction and magnitude vectors, not rooted at any
- particular point;
+- `double ROOT::Math::comp_ellint_3(double n,double k)` - calculates
+ the complete elliptic integral of the third kind (with $0 \le k^2 \le 1$):
+ $$
+ \Pi (n, k, \pi / 2) = \int_{0}^{\pi /2} \frac{d \theta}{(1 - n \sin^2{\theta})\sqrt{1 - k^2 \sin^2{\theta}}}
+ $$
-- `ROOT::Math::`**`PositionVector2D`** and
- `ROOT::Math::`**`PositionVector3D`** template classes modeling the
- points in 2 and 3 dimensions.
+- `double ROOT::Math::conf_hyperg(double a,double b,double z)` -
+ calculates the confluent hyper-geometric functions of the first
+ kind:
+ $$
+ _{1}F_{1}(a;b;z) = \frac{\Gamma(b)}{\Gamma(a)} \sum_{n=0}^{\infty} \frac{\Gamma(a+n)}{\Gamma(b+n)} \frac{z^n}{n!}
+ $$
-For the 4D space-time vectors, we use the same class to model them,
-`ROOT::Math::`**`LorentzVector`**, since we have recognized a limited
-need for modeling the functionality of a 4D point.
+- `double ROOT::Math::conf_hypergU(double a,double b,double z)` -
+ calculates the confluent hyper-geometric
+ functions of the second kind, known also as Kummer function of the second type. It is
+ related to the confluent hyper-geometric function of the first kind:
+ $$
+ U(a,b,z) = \frac{ \pi}{ \sin{\pi b } } \left[ \frac{ _{1}F_{1}(a,b,z) } {\Gamma(a-b+1) } - \frac{ z^{1-b} { _{1}F_{1}}(a-b+1,2-b,z)}{\Gamma(a)} \right]
+ $$
-#### Generic Coordinate System
+- `double ROOT::Math::cyl_bessel_i(double nu,double x)` - calculates
+ the modified Bessel function of the first kind, also called regular
+ modified (cylindrical) Bessel function:
+ $$
+ I_{\nu} (x) = i^{-\nu} J_{\nu}(ix) = \sum_{k=0}^{\infty} \frac{(\frac{1}{2}x)^{\nu + 2k}}{k! \Gamma(\nu + k + 1)}
+ $$
-The vector classes are based on a generic type of coordinate system,
-expressed as a template parameter of the class. Various classes exist to
-describe the various coordinates systems:
+- `double ROOT::Math::cyl_bessel_j(double nu,double x)` - calculates
+ the (cylindrical) Bessel function of the first kind, also called
+ regular (cylindrical) Bessel function:
+ $$
+ J_{\nu} (x) = \sum_{k=0}^{\infty} \frac{(-1)^k(\frac{1}{2}x)^{\nu + 2k}}{k! \Gamma(\nu + k + 1)}
+ $$
-2D coordinate system classes:
+- `double ROOT::Math::cyl_bessel_k(double nu,double x)` - calculates
+ the modified Bessel function of the second kind, also called
+ irregular modified (cylindrical) Bessel function for $x > 0$, $v > 0$:
+ $$
+ K_{\nu} (x) = \frac{\pi}{2} i^{\nu + 1} (J_{\nu} (ix) + iN(ix)) = \left\{ \begin{array}{cl} \frac{\pi}{2} \frac{I_{-\nu}(x) - I_{\nu}(x)}{\sin{\nu \pi}} & \mbox{for non-integral $\nu$} \\ \frac{\pi}{2} \lim{\mu \to \nu} \frac{I_{-\mu}(x) - I_{\mu}(x)}{\sin{\mu \pi}} & \mbox{for integral $\nu$} \end{array} \right.
+ $$
-- **`ROOT::Math::Cartesian2D`**, based on (`x,y`);
+- `double ROOT::Math::cyl_neumann(double nu,double x)` - calculates
+ the (cylindrical) Bessel function of the second kind, also called
+ irregular (cylindrical) Bessel function or (cylindrical) Neumann
+ function:
+ $$
+ N_{\nu} (x) = Y_{\nu} (x) = \left\{ \begin{array}{cl} \frac{J_{\nu} \cos{\nu \pi}-J_{-\nu}(x)}{\sin{\nu \pi}} & \mbox{for non-integral $\nu$} \\ \lim{\mu \to \nu} \frac{J_{\mu} \cos{\mu \pi}-J_{-\mu}(x)}{\sin{\mu \pi}} & \mbox{for integral $\nu$} \end{array} \right.
+ $$
-- **`ROOT::Math::Polar2D`**, based on (`r,phi`);
+- `double ROOT::Math::ellint_1(double k,double phi)` - calculates
+ incomplete elliptic integral of the first kind (with $0 \le k^2 \le 1$):
+ $$
+ K(k) = F(k, \pi / 2) = \int_{0}^{\pi /2} \frac{d \theta}{\sqrt{1 - k^2 \sin^2{\theta}}}
+ $$
-3D coordinate system classes:
+- `double ROOT::Math::ellint_2(double k,double phi)` - calculates
+ the complete elliptic integral of the second kind (with $0 \le k^2 \le 1$):
+ $$
+ E(k) = E(k , \pi / 2) = \int_{0}^{\pi /2} \sqrt{1 - k^2 \sin^2{\theta}} d \theta
+ $$
-- **`ROOT::Math::Cartesian3D`**, based on (`x,y,z`);
+- `double ROOT::Math::ellint_3(double n,double k,double phi)` - calculates
+ the complete elliptic integral of the third kind (with $0 \le k^2 \le 1$):
+ $$
+ \Pi (n, k, \pi / 2) = \int_{0}^{\pi /2} \frac{d \theta}{(1 - n \sin^2{\theta})\sqrt{1 - k^2 \sin^2{\theta}}}
+ $$
-- **`ROOT::Math::Polar3D`**, based on (`r,theta,phi`);
+- `double ROOT::Math::expint(double x)` - calculates the exponential
+ integral:
+ $$
+ Ei(x) = - \int_{-x}^{\infty} \frac{e^{-t}}{t} dt
+ $$
-- **`ROOT::Math::Cylindrical3D`**, based on (`rho,z,phi`)
+- `double ROOT::Math::hyperg(double a,double b,double c,double x)` -
+ calculates Gauss' hyper-geometric function:
+ $$
+ _{2}F_{1}(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a) \Gamma(b)} \sum_{n=0}^{\infty} \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)} \frac{x^n}{n!}
+ $$
-- **`ROOT::Math::CylindricalEta3D`**, based on (`rho,eta,phi`), where
- `eta` is the pseudo-rapidity;
+- `double ROOT::Math::legendre(unsigned l,double x)` - calculates
+ the Legendre polynomials for $l \ge 0$, $|x| \le 1$ in the Rodrigues
+ representation:
+ $$
+ P_{l}(x) = \frac{1}{2^l l!} \frac{d^l}{dx^l} (x^2 - 1)^l
+ $$
-4D coordinate system classes:
+- `double ROOT::Math::riemann_zeta(double x)` - calculates the
+ Riemann zeta function:
+ $$
+ \zeta (x) = \left\{ \begin{array}{cl} \sum_{k=1}^{\infty}k^{-x} & \mbox{for $x > 1$} \\ 2^x \pi^{x-1} \sin{(\frac{1}{2}\pi x)} \Gamma(1-x) \zeta (1-x) & \mbox{for $x < 1$} \end{array} \right.
+ $$
-- **`ROOT::Math::PxPyPzE4D`**, based on based on (`px,py,pz,E`);
+- `double ROOT::Math::sph_bessel(unsigned n,double x)` - calculates
+ the spherical Bessel functions of the first kind (also called
+ regular spherical Bessel functions):
+ $$
+ j_{n}(x) = \sqrt{\frac{\pi}{2x}} J_{n+1/2}(x)
+ $$
-- **`ROOT::Math::PxPyPzM4D`**, based on based on (`px,py,pz,M`);
+- `double ROOT::Math::sph_neumann(unsigned n,double x)` - calculates
+ the spherical Bessel functions of the second kind (also called
+ irregular spherical Bessel functions or spherical Neumann
+ functions):
+ $$
+ n_n(x) = y_n(x) = \sqrt{\frac{\pi}{2x}} N_{n+1/2}(x)
+ $$
-- **`ROOT::Math::PtEtaPhiE4D`**, based on based on (`pt,eta,phi,E`);
+### Probability Density Functions (PDF)
-- **`ROOT::Math::PtEtaPhiM4D`**, based on based on (`pt,eta,phi,M`);
-Users can define the vectors according to the coordinate type, which is
-the most efficient for their use. Transformations between the various
-coordinate systems are available through copy constructors or the
-assignment (=) operator. For maximum flexibility and minimize memory
-allocation, the coordinate system classes are templated on the scalar
-type. To avoid exposing templated parameter to the users, typedefs are
-defined for all types of vectors based on doubles. See in the examples
-for all the possible types of vector classes, which can be constructed
-by users with the available coordinate system types.
+Probability density functions of various distributions. All the
+functions, apart from the discrete ones, have the extra location
+parameter `x0`, which by default is zero. For example, in the case of a
+gaussian `pdf`, `x0` is the `mean`, `mu`, of the distribution. All the
+probability density functions are defined in the header file
+`Math/DistFunc.h` and are part of the `MathCore` libraries. The
+definition of these functions is documented in the
+[reference doc for statistical functions](reference doc for statistical functions):
-#### Coordinate System Tag
+``` {.cpp}
+double ROOT::Math::beta_pdf(double x,double a, double b);
+double ROOT::Math::binomial_pdf(unsigned int k,double p,unsigned int n);
+double ROOT::Math::breitwigner_pdf(double x,double gamma,double x0=0);
+double ROOT::Math::cauchy_pdf(double x,double b=1,double x0=0);
+double ROOT::Math::chisquared_pdf(double x,double r,double x0=0);
+double ROOT::Math::exponential_pdf(double x,double lambda,double x0=0);
+double ROOT::Math::fdistribution_pdf(double x,double n,double m,double x0=0);
+double ROOT::Math::gamma_pdf(double x,double alpha,double theta,double x0=0);
+double ROOT::Math::gaussian_pdf(double x,double sigma,double x0=0);
+double ROOT::Math::landau_pdf(double x,double s,double x0=0);
+double ROOT::Math::lognormal_pdf(double x,double m,double s,double x0=0);
+double ROOT::Math::normal_pdf(double x,double sigma,double x0=0);
+double ROOT::Math::poisson_pdf(unsigned int n,double mu);
+double ROOT::Math::tdistribution_pdf(double x,double r,double x0=0);
+double ROOT::Math::uniform_pdf(double x,double a,double b,double x0=0);
+```
-The 2D and 3D points and vector classes can be associated to a tag
-defining the coordinate system. This can be used to distinguish between
-vectors of different coordinate systems like global or local vectors.
-The coordinate system tag is a template parameter of the
-**`ROOT::Math::`**`DisplacementVector3D` and
-`ROOT::Math::PositionVector3D` (and also for 2D classes). A default tag
-exists for users who do not need this functionality,
-`ROOT::Math::DefaultCoordinateSystemTag`.
+### Cumulative Distribution Functions (CDF)
-#### Transformations
-The transformations are modeled using simple (non-template) classes,
-using double as the scalar type to avoid too large numerical errors. The
-transformations are grouped in rotations (in 3 dimensions), Lorentz
-transformations and Poincare transformations, which are
-translation`/`rotation combinations. Each group has several members
-which may model physically equivalent transformations but with different
-internal representations. Transformation classes can operate on all type
-of vectors by using the operator `() `or the operator `*` and the
-transformations can be combined via the operator `*`. The available
-transformations are:
+For all the probability density functions, we have the corresponding
+cumulative distribution functions and their complements. The functions
+with extension `_cdf` calculate the lower tail integral of the
+probability density function:
-- 3D rotation classes
+$$
+D(x) = \int_{-\infty}^{x} p(x') dx'
+$$
-- rotation described by a 3x3 matrix (**`ROOT::Math::Rotation3D`**)
+while those with the `cdf_c` extension calculate the upper tail of the
+probability density function, so-called in statistics the survival
+function. For example, the function:
-- rotation described by Euler angles (**`ROOT::Math::EulerAngles`**)
+``` {.cpp}
+double ROOT::Math::gaussian_cdf(double x,double sigma,double x0=0);
+```
+evaluates the lower tail of the Gaussian distribution:
+$$
+D(x) = \int_{-\infty}^{x} {1 \over \sqrt{2 \pi \sigma^2}} e^{-(x'-x_0)^2 / 2\sigma^2} dx'
+$$
-- rotation described by a direction axis and an angle
- (**`ROOT::Math::AxisAngle`**)
+while the function:
-- rotation described by a quaternion (**`ROOT::Math::Quaternion`**)
+``` {.cpp}
+double ROOT::Math::gaussian_cdf_c(double x, double sigma, double x0=0);
+```
+evaluates the upper tail of the Gaussian distribution:
+$$
+D(x) = \int_{x}^{+\infty} {1 \over \sqrt{2 \pi \sigma^2}} e^{-(x'-x_0)^2 / 2\sigma^2} dx'
+$$
-- optimized rotation around `x` (**`ROOT::Math::RotationX`**), `y`
- (**`ROOT::Math::RotationY`**) and `z` (**`ROOT::Math::RotationZ`**)
- and described by just one angle.
+The cumulative distributions functions are defined in the header file
+`Math/ProbFunc.h`. The majority of the CDF's are present in the
+`MathCore`, apart from the `chisquared`, `fdistribution`, `gamma` and
+`tdistribution`, which are in the `MathMore` library.
-- 3D transformation: we describe the transformations defined as a
-composition between a rotation and a translation using the class
-**`ROOT::Math::Transform3D`**. It is important to note that
-transformations act differently on vectors and points. The vectors only
-rotate, therefore when applying a transformation (rotation +
-translation) on a vector, only the rotation operates while the
-translation has no effect. The **`Transform3D`** class interface is
-similar to the one used in the CLHEP Geometry package (class
-<HepGeom::Transform3D>).
+#### Inverse of the Cumulative Distribution Functions(Quantiles)
-- Lorentz rotation:
+For almost all the cumulative distribution functions (`_cdf`) and their
+complements (`_cdf_c`) present in the library, we provide the inverse
+functions. The inverse of the cumulative distribution function is called
+in statistics quantile function. The functions with the extension
+`_quantile` calculate the inverse of the cumulative distribution
+function (lower tail integral of the probability density function),
+while those with the *`quantile_c`* extension calculate the inverse of
+the complement of the cumulative distribution (upper tail integral). All
+the inverse distributions are in the MathMore library and are defined in
+the header file `Math/ProbFuncInv.h`.
-- generic Lorentz rotation described by a `4x4` matrix containing a 3D
- rotation part and a boost part (class
- **`ROOT::Math::LorentzRotation`**)
+The following picture illustrates the available statistical functions
+(PDF, CDF and quantiles) in the case of the normal distribution.
-- a pure boost in an arbitrary direction and described by a 4x4
- symmetric matrix or 10 numbers (class **`ROOT::Math::Boost`**)
+
-- boost along the axis:` x `(**`ROOT::Math::BoostX`**),
- `y `(**`ROOT::Math::BoostY`**) and `z `(**`ROOT::Math::BoostZ`**).
+## Numerical Algorithms
-#### Minimal Vector Classes Interface
+ROOT provides C++ classes implementing numerical algorithms to solve a wide set of problem, like:
+
+* Evaluation of function derivatives.
+* Evaluation of integrals.
+* Finding the roots of a function
+* Finding the minimum/maximum of a function
+
+In order to use these algorithm the user needs to provide a function.
+ROOT provides a common way of specifying them via some interfaces
+
+## ROOT::Math Function interfaces
+
+To get a consistency in the mathematical methods within ROOT, there exists a set of interfaces to define the basic behaviour of a mathematical function.
+In order to use the classes presented in this chapter, the mathematical functions defined by the user must inherit from any of the classes seen in the figure:
+
+
+
+
+### One-dimensional Function Interfaces
+
+These interfaces are used for numerical algorithms operating only on one-dimensional functions and cannot be applied to multi-dimensional functions.
+For this case the users needs to define a function object which evaluates in one dimension, and the object will have to derivate from the following:
+
+* `ROOT::Math::IBaseFunctionOneDim`: This class is the most basic function. Provides a method to evaluate the function given a value (simple double) by implementing
+`double operator() (const double )`. The user class defined only needs to reimplement the pure abstract method `double DoEval(double x)`,
+that will do the work of evaluating the function at point x.
+
+Example on how to create a class that represents a mathematical function. The user only has to override two methods from `IBaseFunctionOneDim`:
+
+```{.cpp}
+#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`: Some of the numerical algorithm will need to calculate the derivatives of the function. In these cases, the user will have to provide the neccesary code for
+ this to happen. The interface defined in `IGradientFunctionOneDim` introduced the method `double Derivative(double x)` that will return the derivative of the function at the point `x`. The class
+ inherit by the user will have to implement the abstract method `double DoDerivative(double x)`, leaving the rest of the class untouched.
+
+ Example for implementing a gradient one-dimensional function:
+
+```{.cpp}
+#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
+
+
+The most generic case of a multidimensional function has similar approach. Some examples will be shown next. It is important to notice, that one dimensional functions can be also implemented through
+the interfaces that will be presented here. Nevertheless, the user needs to implement those following the indications of the previous chapter, for algorithm woring exclusivly on one-dimensional
+functions. For algorithms working on both one-dimensional and multi-dimensional functions they should instead use this interface.
+
+* `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 `unsigned int NDim()`. The first ones evaluates the function given the array that represents
+ the multiple variables. The second returns the number of dimensions of the function.
+
+ Example of implementing a basic multi-dimensional function:
+
+```{.cpp}
+#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 plus the calcualtion of the derivative.
+It only adds the `double Derivative(double* x, uint ivar)` method for the user to implement. This method must implement the derivative of the function with respect to the variable indicated with the
+second parameter.
+
+Example of implementing a multi-dimensional gradient function
+
+```{.cpp}
+#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
+
+These interfaces, for evaluating multi-dimensional functions are used for fitting. These interfaces are defined in the header file
+`Math/IParamFunction.h`.
+See also the documentation of the `ROOT::Fit` classes in the Fitting chaper for more information.
+
+* **`ROOT::Math::IParametricFunctionMultiDim`**: Describes a multi dimensional parametric function. Similarly to the one dimensional version, the user needs to provide the
+method `void SetParameters(double* p)` as well as the getter methods `const double * Parameters()` and `uint NPar()`.
+Example of creating a parametric function:
+
+```{.cpp}
+#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`**:
+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 method `ParameterGradient()`. 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. Here is an example:
+
+```{.cpp}
+#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 facilitate the user to insert their own type of function in the needed function interface, helper classes, wrapping the user interface in the
+`ROOT::Math` function interfaces are provided.
+this will avoid the user to re-implement dedicated funcition classes, following the code example shown in the previous paragraphs.
+
+There is one possible wrapper for every interface explained in the previous section.
+The following table indicates the wrapper for the most basic ones:
+
+
+
+
+| **Interface**| **Function Wrapper** |
+|------------------------------------------|------------------------|
+| `ROOT::Math::IBaseFunctionOneDim` | `ROOT::Math::Functor1D` |
+| `ROOT::Math::IGradientFunctionOneDim` | `ROOT::Math::GradFunctor1D` |
+| `ROOT::Math::IBaseFunctionMultiDim` | `ROOT::Math::Functor` |
+| `ROOT::Math::IGradientFunctionMultiDim` | `ROOT::Math::GradFunctor` |
+
+
+
+Thee functor wrapper are defined in the header file `Math/Functor.h`.
+
+#### Wrapping One Dimensional Functions
+
+The `ROOT::Math::Functor1D` is used to wrap one-dimensional functions It can wrap all the following types:
+* A free C function of type `double ()(double )`.
+* Any C++ callable object implemention `double operator()( double )`.
+* A class member function with the correct signature like `double Foo::Eval(double )`. In this case one pass the object pointer and a pointer to the member function `(&Foo::Eval)`.
+
+Example:
+
+```{.cpp}
+#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
+
+The `ROOT::Math::GradFunctor1D` class is used 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 like `Foo::XXX(double )` for the function evaluation and any other member function like `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 object can be a free C function of
+type `double ()(double )`.
+
+#### Wrapping Multi-dimensional Functions
+
+The class `ROOT::Math::Functor` is used to wrap in a very simple and convenient way 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)`.
+
+The function dimension is required when constructing the functor.
+
+Example of using `Functor`:
+```{.cpp}
+#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
+
+The class `ROOT::Math::GradFunctor` is used to wrap in a very C++ callable object to make gradient functions. It can be constructed in three different way:
+* 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 like `Foo::XXX(const double *)` for the function evaluation and any member function like `Foo::XXX(const double *, int icoord)` for the partial derivatives.
+* From an 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.
+
+#### Special case: Wrapping TF1 objects in Parametric Function interfaces
+
+In many cases, the user works with the `TF1` class. The mathematical library in ROOT provides some solutions to wrap these into the interfaces needed by other methods.
+If the desired interface to wrap is one-dimensional, the class to use is `ROOT::Math::WrappedTF1`.
+The default constructor takes a `TF1` reference as an argument, that will be wrapped with the interfaces of a `ROOT::Math::IParametricGradFunctionOneDim`.
+Example:
+```{.cpp}
+#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;
+}
+```
+
+For a TF1 defining a multidimensional function or in case we need to wrap in a multi-dimensional function interface, the class to use is `ROOT::Math::WrappedMultiTF1`.
+Following the usual procedure, setting the `TF1` though the constructor, will wrap it into a `ROOT::Math::IParametricGradFunctionMultiDim`.
+Example:
+
+```{.cpp}
+#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;
+}
+```
+
+
+## Numerical Integration
+
+The algorithms provided by ROOT for numerical integration are implemented following the hierarchy shown in the next image.
+`ROOT::Math::VirtualIntegrator` defines the most basic functionality while the specific behaviours for one or multiple dimensions are implemented in
+`ROOT::Math::VirtualIntegratorOneDim` and `ROOT::Math::VirtualIntegratorMultiDim`.
+These interfaces define the integrator functionality with abstract methods to set the function, to compute the integral or to set the integration tolerance.
+These methods must be implemented in the concrete classes existing for the different integration algorithms.
+The user cannot create directly these virtual integrator interfaces. He needs to create the
+`ROOT::Math::IntegratorOneDim` class for integrating one-dimensional functions and `ROOT::Math::IntegratorMultiDim` for multi-dimensional functions.
+Through the ROOT Plug-In Manager, the user can initialize `ROOT::Math::IntegratorOneDim` or `ROOT::Math::IntegratorMultiDim` with
+any of the concrete integration classes without dealing with them directly.
+These two classes provide the same interface as in `VirtualIntegratorOneDim` and `VirtualIntegratorMultiDim`, but with the possibility to choose in the constructor,
+which method will be used to perform the integration.
+
+The method to set the function to be integrated, must be of the function interface type described before.
+`ROOT::Math::IBaseFunctionOneDimFunction` is used for `ROOT::Math::IBaseFunctionMultiDim` and
+The only difference between the `ROOT::Math::IntegratorOneDim` and `ROOT::Math::IntegratorMultiDim` resides
+in the dimensionality of that function and some specific that will be seen afterwards for the one dimensional one.
+
+
+
+The rest of the classes shown above in the diagram are the specialized classes provided. Each one implements a different method that will be explained in detail. It is important
+to notice that the two grayed classes (the one which name starts by GSL) are part of the *MathMore* library.
+We will later show in more detail the differences between the implementations.
+
+
+### Integration of One-dimensional Functions
+
+#### Using `ROOT::Math::IntegratorOneDim`
+
+Here is a code example on how to use the `ROOT::Math::IntegratorOneDim` class
+(note that the class is defined in the header file `Math/Integrator.h`). In this example we create
+different instance of the class 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`
+
+```{.cpp}
+#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;
+}
+```
+
+### One-dimensional Integration Algorithms
+
+Here we provide a brief description of the different integration algorithms, which are also
+implemented as separate classes. The algorithms can be instantiated 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
+
+It uses the most basic Gaussian integration algorithm, it uses the 8-point and the 16-point Gaussian
+quadrature approximations. It is derived from the `DGAUSS` rountine of the *CERNLIB* by S. Kolbig.
+This class
+Here is an example of using directly the `GaussIntegrator` class
+
+```{.cpp}
+#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;
+
+ return 0;
+}
+```
+#### ROOT::Math::GaussLegendreIntegrator
+
+This class implementes the Gauss-Legendre quadrature formulas. This sort of numerical methods requieres that the user specifies the number of intermediate function points
+used in the calculation of the integral. It will automatically determine the coordinates and weights of such points before performing the integration.
+We can use the example above, but replacing the creation of a `ROOT::Math::GaussIntegrator` object with `ROOT::Math::GaussLegendreIntegrator`.
+
+#### ROOT::Math::GSLIntegrator
+
+This is a 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. For a detail description of the GSL methods visit the GSL user guide
+This class implements the best algorithms for numerical integration for one dimensional functions. We encourage the use it as the main option, bearing in mind that it uses code from the
+GSL library, wich is provided in the *MathMore* library of ROOT.
+
+The interface to use is the same as above. We have now the possibility to specify a different integration algorithm in the constructor of the `ROOT::Math::GSLIntegrator` class.
+```{.cpp}
+ROOT::Math::GSLIntegrator ig(ROOT::Math::Integration::kADAPTIVE, ROOT::Math::Integration::kGAUSS51); // create the adaptive integrator with the 51 point rule
+ig.SetRelTolerance(1.E-6); // set relative tolerance
+ig.SetAbsTolerance(1.E-6); // set absoulte tolerance
+```
+
+The algorithm is controlled by the given absolute and relative tolerance. The iterations are continued until the following condition is satisfied
+$$
+absErr <= max ( epsAbs, epsRel * Integral)
+$$
+Where *absErr* is an estimate of the absolute error (it can be retrieved with `GSLIntegrator::Error()`) and *Integral* is the estimate of the function integral
+(it can be obtained with `GSLIntegrator::Result()`)
+
+The possible integration algorithm types to use with the GSLIntegrator are the following. More information is provided in the `GSL` users documentation.
+* `ROOT::Math::Integration::kNONADAPTIVE` : based on `gsl_integration_qng`. It is a non-adaptive procedure which uses fixed Gauss-Kronrod-Patterson abscissae
+to sample the integrand at a maximum of 87 points. It is provided for fast integration of smooth functions.
+* `ROOT::Math::Integration::kADAPTIVE`: based on `gsl_integration_qag`. It is an adaptiva Gauss-Kronrod integration algorithm, the integration region is divided into subintervals, and on each
+iteration the subinterval with the largest estimated error is bisected. It is possible to specify the integration rule as an extra enumeration parameter. The possible rules are
+ * `Integration::kGAUSS15` : 15 points Gauss-Konrod rule (value = 1)
+ * `Integration::kGAUSS21` : 21 points Gauss-Konrod rule (value = 2)
+ * `Integration::kGAUSS31` : 31 points Gauss-Konrod rule (value = 3)
+ * `Integration::kGAUSS41` : 41 points Gauss-Konrod rule (value = 4)
+ * `Integration::kGAUSS51` : 51 points Gauss-Konrod rule (value = 5)
+ * `Integration::kGAUSS61` : 61 points Gauss-Konrod rule (value = 6)
+ The higher-order rules give better accuracy for smooth functions, while lower-order rules save time when the function contains local difficulties, such as discontinuities. If no integration rule
+ is passed, the 31 points rule is used as default.
+
+* `ROOT::Math::Integration::kADAPTIVESINGULAR`: based on `gsl_integration_qags`. It is an integration type which can be used in the case of the presence of singularities.It uses the
+ Gauss-Kronrod 21-point integration rule. This is the default algorithm
+
+Note that when using the common `ROOT::Math::IntegratorOneDIm` class the enumeration type defining the algorithm must be defined in the namespace `ROOT::Math::IntegrationOneDim` (to distinguish from
+the multi-dimensional case) and the rule enumeration (or its corresponding integer) can be passed in the constructor of the `ROOT::Math::IntegratorOneDIm`.
+
+### Multi-dimensional Integration
+
+The multi-dimensional integration algorithm should be applied to functions with dimension larger than one.
+Adaptive multi-dimensional integration works for low function dimension, while MC integration can be applied to higher dimensions.
+
+#### Using `ROOT::Math::IntegratorMultiDim`
+
+Here is a code example on how to use the `ROOT::Math::IntegratorOneDim` class
+(note that the class is defined in the header file `Math/Integrator.h`). In this example we create
+different instance of the class using some of the available algorithms in ROOT.
+
+```{.cpp}
+#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;
+}
+```
+
+#### Multi-dimensions Integration Algorithms
+
+Here is the types, that can be specified as enumeration and the corresponding classes
+
+| **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` |
+
+The control parameters for the integration algorithms can be specified using the
+`ROOT::Math::IntegratorMultiDimOptions` class. Static methods are provided to change the default values.
+It is possible to print the list of default control parameters using the `ROOT::Math::IntegratorMultiDimOptions::Print` function.
+Example:
+```{.cpp}
+ROOT::Math::IntegratorMultiDimOptions opt;
+opt.Print();
+ Integrator Type : ADAPTIVE
+ Absolute tolerance : 1e-09
+ Relative tolerance : 1e-09
+ Workspace size : 100000
+ (max) function calls : 100000
+```
+Depending on the algorithm, some of the control parameters might have no effect.
+
+#### `ROOT::Math::AdaptiveIntegratorMultiDim`
+
+This class implements an adaptive quadrature integration method for multi dimensional functions. It is described in this 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*.
+It is part of the *MathCore* library.
+The user can control the relative and absolute tolerance and the maximum allowed number of function evaluation.
-We have tried to keep the interface to a minimal level by:
-- Avoiding methods that provide the same functionality but use
- different names (like `getX()` and `x()`).
+#### `ROOT::Math::GSLMCIntegrator`
-- Minimizing the number of setter methods, avoiding methods, which can
- be ambiguous and can set the vector classes in an inconsistent
- state. We provide only methods which set all the coordinates at the
- same time or set only the coordinates on which the vector is based,
- for example `SetX()` for a Cartesian vector. We then enforce the use
- of transformations as rotations or translations (additions) for
- modifying the vector contents.
+It 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.
+For a detail description of the GSL methods visit the GSL users guide.
+Specific configuration options (documented in the GSL user guide) for the `ROOT::Math::GSLMCIntegration` can be set directly in the class, or when using it via the `ROOT::Math::IntegratorMultiDim`
+interface, can be defined using the `ROOT::Math::IntegratorMultiDimOptions`.
-- The majority of the functionality, which is present in the CLHEP
- package, involving operations on two vectors, is moved in separated
- helper functions (see `ROOT::Math::VectorUtil`). This has the
- advantage that the basic interface will remain more stable with time
- while additional functions can be added easily.
-#### Naming Convention
-As part of ROOT, the `GenVector` package adheres to the prescribed ROOT
-naming convention, with some (approved) exceptions, as described here:
+## Function Derivation
-- Every class and function is in the **`ROOT::Math`** namespace.
+There are in ROOT only two classes to perform numerical derivation. One of them is in the MathCore library while the other is in the MathMore wrapping an integration function from the GSL library.
+* RichardsonDerivator: Implements the Richardson method for numerical integration. It can calculate up to the third derivative of a function.
+* GSLDerivator of *MathMore* based on GSL.
-- Member function names start with upper-case letter, apart some
- exceptions (see the next section about CLHEP compatibility).
+## Numerical Minimization
-#### Compatibility with CLHEP Vector Classes
+The algorithms provided by ROOT for numerical integration are implemented following the hierarchy shown in the next image. The left branch of classes are used for one dimensional minimization, while
+the right one is used for multidimensional minimization. In the case of multidimensional minimization we have also the classes `TMinuitMinimizer` implemented using `TMinuit`, `TFumiliMinimizer`
+implemented using `TFumili` for least square or likelihood minimizations.
+We encourage the use of the GSL algorithms for one dimensional minimization and `Minuit2` (or the old version`Minuit`) for multi dimensional minimization.
-- For backward compatibility with CLHEP the vector classes can be
- constructed from a CLHEP `HepVector` or **`HepLorentzVector`**, by
- using a template constructor, which requires only that the classes
- implement the accessors` x()`, `y()`, and `z()` (and `t()` for the
- 4D).
+
-- We provide vector member function with the same naming convention as
- CLHEP for the most used functions like `x()`, `y()` and `z()`.
-#### Connection to Linear Algebra Package
+### One-Dimensional Minimization
-In some use cases, like in track reconstruction, it is needed to use the
-content of the vector and rotation classes in conjunction with linear
-algebra operations. We prefer to avoid any direct dependency to any
-linear algebra package. However, we provide some hooks to convert to and
-from linear algebra classes. The vector and the transformation classes
-have methods which allow to get and set their data members (like
-`SetCoordinates` and `GetCoordinates`) passing either a generic iterator
-or a pointer to a contiguous set of data, like a C array. This allows an
-easy connection with the linear algebra package, which in turn, allows
-creation of matrices using C arrays (like the ROOT **`TMatrix`**
-classes) or iterators (`SMatrix` classes). Multiplication between linear
-algebra matrices and `GenVector` vectors is possible by using the
-template free functions `ROOT::Math::VectorUtil::Mult`. This function
-works for any linear algebra matrix, which implements the operator
-(`i,j`) and with first matrix element at `i=j=0`.
+These algorithms are for finding 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.
+The algorithms supported are only bracketing algorithm which do not use derivatives information.
-### Example: 3D Vector Classes
+Two classes exist. One in the *MathCore* library implementing the Brent method (not using the derivatives)
+and one in the *MathMore* library implementing several different methods, using in some case the derivatives.
-To avoid exposing template parameter to the users, typedef's are defined
-for all types of vectors based on double's and float's. To use them, one
-must include the header file `Math/Vector3D.h`. The following typedef's,
-defined in the header file `Math/Vector3Dfwd.h`, are available for the
-different instantiations of the template class
-`ROOT::Math::`**`DisplacementVector3D`**:
+#### `ROOT::Math::BrentMinimizer1D`
-- `ROOT::Math::`**`XYZVector`** vector based on `x,y,z` coordinates
- (Cartesian) in double precision
+This class implements the Brent method to minimize one-dimensional function.
+An interval containing the function minimum must be provided.
+Here is an example where we define the function to minimize as a *lambda* function
+(requires C++11). The function to minimize must be given to the class implementing the algorithm as a
+`ROOT::Math::IBaseFunctionOneDim` object.
-- `ROOT::Math::`**`XYZVectorF`** vector based on `x,y,z` coordinates
- (Cartesian) in float precision
+```{.cpp}
+ 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);
+ cout << "Minimum: f(" << bm.XMinimum() << ") = " <<bm.FValMinimum() << endl;
+```
-- `ROOT::Math::`**`Polar3DVector`** vector based on `r,theta,phi`
- coordinates (polar) in double precision
+Note that when setting the function to minimize, one needs to 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::`**`Polar3DVectorF`** vector based on `r,theta,phi`
- coordinates (polar) in float precision
+#### `ROOT::Math::GSLMInimizer1D`
-- `ROOT::Math::`**`RhoZPhiVector`** vector based on `rho,z,phi`
- coordinates (cylindrical) in double precision
+This class wraps two different methods from the GSL.
+The algorithms which can be choosen at construction time are *GOLDENSECTION*, which is the simplest method
+but the slowest and *BRENT* (the default one) which combines the golden section with a parabolic interpolation.
+The algorithm can be chosen as a different enumeration in the constructor:
+* `ROOT::Math::Minim1D::kBRENT` for the Brent algorithm (default)
+* `ROOT::Math::Minim1D::kGOLDENSECTION` for the golden section algorithm
-- `ROOT::Math::`**`RhoZPhiVectorF`** vector based on `rho,z,phi`
- coordinates (cylindrical) in float precision
+```{.cpp}
+ // this makes class with the default Brent algorithm
+ ROOT::Math::GSLMinimizer1D minBrent;
+ // this make the class with the Golden Section algorithm
+ ROOT::Math::GSLMinimizer1D minGold(ROOT::Math::Minim1D::kGOLDENSECTION);
+```
-- `ROOT::Math::`**`RhoEtaPhiVector`** vector based on `rho,eta,phi`
- coordinates (cylindrical using `eta` instead of `z`) in double
- precision
+The interface to set the function and to minimize is the same as in the case of the `BrentMinimizer1D`.
-- `ROOT::Math::`**`RhoEtaPhiVectorF`** vector based on `rho,eta,phi`
- coordinates (cylindrical using `eta` instead of `z`) in float
- precision
+#### Using the TF1 class
-#### Constructors and Assignment
+It is possible to perform the one-dimensional minimization/maximization of a function by using directly the function class in ROOT, `TF1` of the *Hist* library.
+The minmization is implemented in `TF1` using the BrentMInimizer1D and available with the class member functions
+* `TF1::GetMinimum`/`TF1::GetMaximum` to find the function minimum/maximum value
+* `TF1::GetMinimumX`/`TF1::GetMaximumX` to find the x value corresponding at the function minimum.
-The following declarations are available:
+The interval to search for the minimum (the default is the `TF1` range), tolerance and maximum iterations can be provided as optional parameters of the
+`TF1::GetMinimum/Maximum` functions.
-``` {.cpp}
-XYZVector v1; //an empty vector (x=0, y=0, z=0)
-XYZVector v2(1,2,3); //vector with x=1, y=2, z=3;
-Polar3DVector v3(1,PI/2,PI); //vector with r=1, theta=PI/2, phi=PI
-RhoEtaPHiVector v4(1,2, PI); //vector with rho=1, eta=2, phi=PI
-```
-Note that each vector type is constructed by passing its coordinate
-representation, so a `XYZVector(1,2,3)` is different from a
-`Polar3DVector(1,2,3)`. In addition, the vector classes can be
-constructed by any vector, which implements the accessors `x()`, `y()`
-and `z()`. This can be another 3D vector based on a different coordinate
-system type. It can be even any vector of a different package, like the
-CLHEP **`HepThreeVector`** that implements the required signature.
+### Multi-Dimensional Minimization
-``` {.cpp}
- XYZVector v1(1,2,3);
- RhoEtaPhiVector r2(v1);
- CLHEP::HepThreeVector q(1,2,3);
- XYZVector v3(q);
-```
+All the algorithms for multi-dimensional minimization are implementing the `ROOT::Math::Minimizer`
+interface and they can be used in the same way and one can switch between minimizer at run-time.
+The minimizer concrete class can be in different ROOT libraries and they can be instantiate using the ROOT
+plug-in manager.
+More information on multi-dimensional minimization is provided in the Fitting Histogram chapter.
-#### Coordinate Accessors
+## ROOT Finder Algorithms
-All coordinate accessors are available through the class
-`ROOT::Math::`**`DisplacementVector3D`**:
+The function must be given to the class implementing the algorithm as a
+`ROOT::Math::IBaseFunctionOneDim` object.
+Some of the algorithm requires the derivatives of the function.
+In that case a `ROOT::Math::IGradientFunctionOneDim` object must be provided.
-``` {.cpp}
- //returns cartesian components for the cartesian vector v1
- v1.X(); v1.Y(); v1.Z();
- //returns cylindrical components for the cartesian vector v1
- v1.Rho(); v1.Eta(); v1.Phi();
- //returns cartesian components for the cylindrical vector r2
- r2.X(); r2.Y(); r2.Z()
-```
-In addition, all the 3 coordinates of the vector can be retrieved with
-the `GetCoordinates` method:
-``` {.cpp}
- double d[3];
- //fill d array with (x,y,z) components of v1
- v1.GetCoordinates(d);
- //fill d array with (r,eta,phi) components of r2
- r2.GetCoordinates(d);
- std::vector vc(3);
- //fill std::vector with (x,y,z) components of v1
- v1.GetCoordinates(vc.begin(),vc.end());
-```
-See the reference documentation of
-`ROOT::Math::`**`DisplacementVector3D`** for more details on all the
-coordinate accessors.
-#### Setter Methods
-One can set only all the three coordinates via:
+## Generic Vectors for 2, 3 and 4 Dimensions (GenVector)
-``` {.cpp}
- v1.SetCoordinates(c1,c2,c3); // (x,y,z) for a XYZVector
- r2.SetCoordinates(c1,c2,c3); // r,theta,phi for a Polar3DVector
- r2.SetXYZ(x,y,z); // 3 cartesian components for Polar3DVector
-```
-Single coordinate setter methods are available for the basic vector
-coordinates, like `SetX()` for a **`XYZVector`** or `SetR()` for a polar
-vector. Attempting to do a `SetX()` on a polar vector will not compile.
+`GenVector` is a package intended to represent vectors and their
+operations and transformations, such as rotations and Lorentz
+transformations, in 3 and 4 dimensions. The 3D space is used to describe
+the geometry vectors and points, while the 4D space-time is used for
+physics vectors representing relativistic particles. These 3D and 4D
+vectors are different from vectors of the linear algebra package, which
+describe generic N-dimensional vectors. Similar functionality is
+currently provided by the CLHEP <Vector> and <Geometry> packages and the
+ROOT Physics vector classes (See "Physics Vectors"). It also re-uses
+concepts and ideas from the CMS
+[Common Vector package](Common Vector package). In contrast to CLHEP or
+the ROOT physics libraries, `GenVector` provides class templates for
+modeling the vectors. The user can control how the vector is internally
+represented. This is expressed by a choice of coordinate system, which
+is supplied as a template parameter when the vector is constructed.
+Furthermore, each coordinate system is itself a template, so that the
+user can specify the underlying scalar type.
-``` {.cpp}
- XYZVector v1;
- v1.SetX(1); //OK setting x for a Cartesian vector
- Polar3DVector v2;
- v2.SetX(1); //ERROR: cannot set X for a Polar vector.
- //Method will not compile
- v2.SetR(1); //OK setting r for a Polar vector
-```
+The `GenVector` classes do not inherit from **`TObject`**, therefore
+cannot be used as in the case of the physics vector classes in ROOT
+collections.
-In addition, there are setter methods from C arrays or iterator
+In addition, to optimize performances, no virtual destructors are
+provided. In the following paragraphs, the main characteristics of
+`GenVector` are described. A more detailed description of all the
+`GenVector` classes is available also at
+<http://seal.cern.ch/documents/mathlib/GenVector.pdf>
-``` {.cpp}
- double d[3] = {1.,2.,3.};
- XYZVector v;
- // set (x,y,z) components of v using values from d
- v.SetCoordinates(d);
-```
+### Main Characteristics
-or, for example, from an `std::vector` using the iterator
-``` {.cpp}
- std::vector w(3);
- // set (x,y,z) components of v using values from w
- v.SetCoordinates(w.begin(),w.end());
-```
+#### Optimal Runtime Performances
-#### Arithmetic Operations
+We try to minimize any overhead in the run-time performance. We have
+deliberately avoided the use of any virtual function and even virtual
+destructors in the classes. In addition, as much as possible functions
+are defined as inline. For this reason, we have chosen to use template
+classes to implement the `GenVector` concepts instead of abstract or
+base classes and virtual functions. It is then recommended to avoid
+using the `GenVector` classes polymorphically and developing classes
+inheriting from them.
-The following operations are possible between vector classes, even of
-different coordinate system types: (`v1,v2` are any type of
-**`ROOT::Math::DisplacementVector3D`** classes, `v3` is the same type
-of `v1`; `a` is a scalar value)
+#### Points and Vector Concept
-``` {.cpp}
- v1 += v2;
- v1 -= v2;
- v1 = - v2;
- v1 *= a;
- v1 /= a;
- v2 = a * v1;
- v2 = v1 / a;
- v2 = v1 * a;
- v3 = v1 + v2;
- v3 = v1 - v2;
-```
+Mathematically vectors and points are two distinct concepts. They have
+different transformations, as vectors only rotate while points rotate
+and translate. You can add two vectors but not two points and the
+difference between two points is a vector. We then distinguish for the 3
+dimensional case, between points and vectors, modeling them with
+different classes:
-#### Comparison
+- `ROOT::Math::`**`DisplacementVector2D`** and
+ `ROOT::Math::`**`DisplacementVector3D`** template classes describing
+ 2 and 3 component direction and magnitude vectors, not rooted at any
+ particular point;
-For `v1` and `v2` of the same type (same coordinate system and same
-scalar type):
+- `ROOT::Math::`**`PositionVector2D`** and
+ `ROOT::Math::`**`PositionVector3D`** template classes modeling the
+ points in 2 and 3 dimensions.
-``` {.cpp}
- v1 == v2;
- v1 != v2;
-```
+For the 4D space-time vectors, we use the same class to model them,
+`ROOT::Math::`**`LorentzVector`**, since we have recognized a limited
+need for modeling the functionality of a 4D point.
-#### Dot and Cross Product
+#### Generic Coordinate System
-We support the dot and cross products, through the `Dot()` and `Cross()`
-method, with any vector (`q`) implementing `x()`, `y()` and `z()`.
+The vector classes are based on a generic type of coordinate system,
+expressed as a template parameter of the class. Various classes exist to
+describe the various coordinates systems:
-``` {.cpp}
- XYZVector v1(x,y,z);
- double s = v1.Dot(q);
- XYZVector v2 = v1.Cross(q);
-```
+2D coordinate system classes:
-Note that the multiplication between two vectors using the operator `*`
-is not supported because it is ambiguous.
+- **`ROOT::Math::Cartesian2D`**, based on (`x,y`);
-#### Other Methods
+- **`ROOT::Math::Polar2D`**, based on (`r,phi`);
-``` {.cpp}
- XYZVector u = v1.Unit(); //return unit vector parallel to v1
-```
+3D coordinate system classes:
-### Example: 3D Point Classes
+- **`ROOT::Math::Cartesian3D`**, based on (`x,y,z`);
+- **`ROOT::Math::Polar3D`**, based on (`r,theta,phi`);
-To use all possible types of 3D points one must include the header file
-`Math/Point3D.h`. The following typedef's defined in the header file
-`Math/Point3Dfwd.h`, are available for different instantiations of the
-template class **`ROOT::Math`**`::`**`PositionVector3D`**:
+- **`ROOT::Math::Cylindrical3D`**, based on (`rho,z,phi`)
-- `ROOT::Math::`**`XYZPoint`** point based on `x`, `y`, `z`
- coordinates (Cartesian) in double precision
+- **`ROOT::Math::CylindricalEta3D`**, based on (`rho,eta,phi`), where
+ `eta` is the pseudo-rapidity;
-- `ROOT::Math::`**`XYZPointF`** point based on `x`, `y`, `z`
- coordinates (Cartesian) in float precision
+4D coordinate system classes:
-- `ROOT::Math::`**`Polar3DPoint`** point based on `r`, `theta`, `phi`
- coordinates (polar) in double precision
+- **`ROOT::Math::PxPyPzE4D`**, based on based on (`px,py,pz,E`);
-- `ROOT::Math::`**`Polar3DPointF`** point based on `r`, `theta`, `phi`
- coordinates (polar) in float precision
+- **`ROOT::Math::PxPyPzM4D`**, based on based on (`px,py,pz,M`);
-- `ROOT::Math::`**`RhoZPhiPoint`** point based on `rho`, `z`, `phi`
- coordinates (cylindrical using `z`) in double precision
+- **`ROOT::Math::PtEtaPhiE4D`**, based on based on (`pt,eta,phi,E`);
-- `ROOT::Math::`**`RhoZPhiPointF`** point based on `rho`, `z`, `phi`
- coordinates (cylindrical using `z`) in float precision
+- **`ROOT::Math::PtEtaPhiM4D`**, based on based on (`pt,eta,phi,M`);
-- `ROOT::Math::`**`RhoEtaPhiPoint`** point based on `rho`, `eta`,
- `phi` coordinates (cylindrical using eta instead of `z`) in double
- precision
+Users can define the vectors according to the coordinate type, which is
+the most efficient for their use. Transformations between the various
+coordinate systems are available through copy constructors or the
+assignment (=) operator. For maximum flexibility and minimize memory
+allocation, the coordinate system classes are templated on the scalar
+type. To avoid exposing templated parameter to the users, typedefs are
+defined for all types of vectors based on doubles. See in the examples
+for all the possible types of vector classes, which can be constructed
+by users with the available coordinate system types.
-- `ROOT::Math::`**`RhoEtaPhiPointF`** point based on `rho`, `eta`,
- `phi` coordinates (cylindrical using eta instead of `z`) in float
- precision
+#### Coordinate System Tag
-#### Constructors and Assignment
+The 2D and 3D points and vector classes can be associated to a tag
+defining the coordinate system. This can be used to distinguish between
+vectors of different coordinate systems like global or local vectors.
+The coordinate system tag is a template parameter of the
+**`ROOT::Math::`**`DisplacementVector3D` and
+`ROOT::Math::PositionVector3D` (and also for 2D classes). A default tag
+exists for users who do not need this functionality,
+`ROOT::Math::DefaultCoordinateSystemTag`.
-The following declarations are available:
+#### Transformations
-``` {.cpp}
- XYZPoint p1; //an empty vector (x=0, y=0, z=0)
- XYZPoint p2(1,2,3); //
-```
+The transformations are modeled using simple (non-template) classes,
+using double as the scalar type to avoid too large numerical errors. The
+transformations are grouped in rotations (in 3 dimensions), Lorentz
+transformations and Poincare transformations, which are
+translation`/`rotation combinations. Each group has several members
+which may model physically equivalent transformations but with different
+internal representations. Transformation classes can operate on all type
+of vectors by using the operator `() `or the operator `*` and the
+transformations can be combined via the operator `*`. The available
+transformations are:
-Note that each point type is constructed by passing its coordinate
-representation, so a `XYZPoint(1,2,3)` is different from a
-`Polar3DPoint(1,2,3)`. In addition the point classes can be constructed
-by any vector, which implements the accessors `x()`, `y()` and `z()`.
-This can be another 3D point based on a different coordinate system type
-or even any vector of a different package, like the CLHEP
-**`HepThreePoint`** that implements the required signatures.
+- 3D rotation classes
-``` {.cpp}
- XYZPoint p1(1,2,3);
- RhoEtaPHiPoint r2(v1);
- CLHEP::HepThreePoint q(1,2,3);
- XYZPoint p3(q);
-```
+- rotation described by a 3x3 matrix (**`ROOT::Math::Rotation3D`**)
-#### Coordinate Accessors and Setter Methods
+- rotation described by Euler angles (**`ROOT::Math::EulerAngles`**)
-For the points classes we have the same getter and setter methods as for
-the vector classes. See "Example: 3D Vector Classes".
+- rotation described by a direction axis and an angle
+ (**`ROOT::Math::AxisAngle`**)
-#### Point-Vector Operations
+- rotation described by a quaternion (**`ROOT::Math::Quaternion`**)
-The following operations are possible between points and vector classes:
-(`p1`, `p2` and `p3` are instantiations of the
-`ROOT::Math::`**`PositionVector3D`** objects with `p1` and `p3` of the
-same type; `v1` and `v2` are `ROOT::Math::`**`DisplacementVector3D`**
-objects).
+- optimized rotation around `x` (**`ROOT::Math::RotationX`**), `y`
+ (**`ROOT::Math::RotationY`**) and `z` (**`ROOT::Math::RotationZ`**)
+ and described by just one angle.
-``` {.cpp}
- p1 += v1;
- p1 -= v1;
- p3 = p1 + v1; // p1 and p3 are the same type
- p3 = v1 + p1; // p3 is based on the same coordinate system as v1
- p3 = p1 - v1;
- p3 = v1 - p1;
- v2 = p1 - p2; // difference between points returns a vector v2
- // based on the same coordinate system as p1
-```
+- 3D transformation: we describe the transformations defined as a
+composition between a rotation and a translation using the class
+**`ROOT::Math::Transform3D`**. It is important to note that
+transformations act differently on vectors and points. The vectors only
+rotate, therefore when applying a transformation (rotation +
+translation) on a vector, only the rotation operates while the
+translation has no effect. The **`Transform3D`** class interface is
+similar to the one used in the CLHEP Geometry package (class
+<HepGeom::Transform3D>).
-Note that the addition between two points is **NOT** possible and the
-difference between points returns a vector.
+- Lorentz rotation:
-#### Other Operations
+- generic Lorentz rotation described by a `4x4` matrix containing a 3D
+ rotation part and a boost part (class
+ **`ROOT::Math::LorentzRotation`**)
-Exactly as for the 3D Vectors, the following operations are allowed:
+- a pure boost in an arbitrary direction and described by a 4x4
+ symmetric matrix or 10 numbers (class **`ROOT::Math::Boost`**)
-- comparison of points
+- boost along the axis:` x `(**`ROOT::Math::BoostX`**),
+ `y `(**`ROOT::Math::BoostY`**) and `z `(**`ROOT::Math::BoostZ`**).
-- scaling and division of points with a scalar
+#### Minimal Vector Classes Interface
-- dot and cross product with any type of vector
+We have tried to keep the interface to a minimal level by:
-### Example: LorentzVector Classes
+- Avoiding methods that provide the same functionality but use
+ different names (like `getX()` and `x()`).
+
+- Minimizing the number of setter methods, avoiding methods, which can
+ be ambiguous and can set the vector classes in an inconsistent
+ state. We provide only methods which set all the coordinates at the
+ same time or set only the coordinates on which the vector is based,
+ for example `SetX()` for a Cartesian vector. We then enforce the use
+ of transformations as rotations or translations (additions) for
+ modifying the vector contents.
+- The majority of the functionality, which is present in the CLHEP
+ package, involving operations on two vectors, is moved in separated
+ helper functions (see `ROOT::Math::VectorUtil`). This has the
+ advantage that the basic interface will remain more stable with time
+ while additional functions can be added easily.
-As in the 3D case, typedef's are defined for user convenience. and can
-be used by including the header file `Math/Vector4D.h`. The following
-typedef's, defined in the header file `Math/Vector4Dfwd.h`, are
-available for the different instantiations of the template class
-**`ROOT::Math::LorentzVector`**:
+#### Naming Convention
-- `ROOT::Math::`**`XYZTVector`** vector based on `x`, `y`, `z`, `t`
- coordinates (Cartesian) in double precision
+As part of ROOT, the `GenVector` package adheres to the prescribed ROOT
+naming convention, with some (approved) exceptions, as described here:
-- `ROOT::Math::`**`XYZTVectorF`** vector based on `x`, `y`, `z`, `t`
- coordinates (Cartesian) in float precision
+- Every class and function is in the **`ROOT::Math`** namespace.
-- `ROOT::Math::`**`PtEtaPhiEVector`** vector based on `pt(rho)`,
- `eta`, `phi` and `E(t)` coordinates in double precision
+- Member function names start with upper-case letter, apart some
+ exceptions (see the next section about CLHEP compatibility).
-- `ROOT::Math::`**`PtEtaPhiMVector`** vector based on `pt(rho)`,
- `eta`, `phi` and `M(t)` coordinates in double precision
+#### Compatibility with CLHEP Vector Classes
-- `ROOT::Math::`**`PxPyPzMVector`** vector based on `px`, `py`, `pz`
- and `M(mass)` coordinates in double precision
+- For backward compatibility with CLHEP the vector classes can be
+ constructed from a CLHEP `HepVector` or **`HepLorentzVector`**, by
+ using a template constructor, which requires only that the classes
+ implement the accessors` x()`, `y()`, and `z()` (and `t()` for the
+ 4D).
-The metric used for all the LorentzVector is (`-,-,-,+`) .
+- We provide vector member function with the same naming convention as
+ CLHEP for the most used functions like `x()`, `y()` and `z()`.
+
+#### Connection to Linear Algebra Package
+
+In some use cases, like in track reconstruction, it is needed to use the
+content of the vector and rotation classes in conjunction with linear
+algebra operations. We prefer to avoid any direct dependency to any
+linear algebra package. However, we provide some hooks to convert to and
+from linear algebra classes. The vector and the transformation classes
+have methods which allow to get and set their data members (like
+`SetCoordinates` and `GetCoordinates`) passing either a generic iterator
+or a pointer to a contiguous set of data, like a C array. This allows an
+easy connection with the linear algebra package, which in turn, allows
+creation of matrices using C arrays (like the ROOT **`TMatrix`**
+classes) or iterators (`SMatrix` classes). Multiplication between linear
+algebra matrices and `GenVector` vectors is possible by using the
+template free functions `ROOT::Math::VectorUtil::Mult`. This function
+works for any linear algebra matrix, which implements the operator
+(`i,j`) and with first matrix element at `i=j=0`.
+
+### Example: 3D Vector Classes
+
+
+To avoid exposing template parameter to the users, typedef's are defined
+for all types of vectors based on double's and float's. To use them, one
+must include the header file `Math/Vector3D.h`. The following typedef's,
+defined in the header file `Math/Vector3Dfwd.h`, are available for the
+different instantiations of the template class
+`ROOT::Math::`**`DisplacementVector3D`**:
+
+- `ROOT::Math::`**`XYZVector`** vector based on `x,y,z` coordinates
+ (Cartesian) in double precision
+
+- `ROOT::Math::`**`XYZVectorF`** vector based on `x,y,z` coordinates
+ (Cartesian) in float precision
+
+- `ROOT::Math::`**`Polar3DVector`** vector based on `r,theta,phi`
+ coordinates (polar) in double precision
+
+- `ROOT::Math::`**`Polar3DVectorF`** vector based on `r,theta,phi`
+ coordinates (polar) in float precision
+
+- `ROOT::Math::`**`RhoZPhiVector`** vector based on `rho,z,phi`
+ coordinates (cylindrical) in double precision
+
+- `ROOT::Math::`**`RhoZPhiVectorF`** vector based on `rho,z,phi`
+ coordinates (cylindrical) in float precision
+
+- `ROOT::Math::`**`RhoEtaPhiVector`** vector based on `rho,eta,phi`
+ coordinates (cylindrical using `eta` instead of `z`) in double
+ precision
+
+- `ROOT::Math::`**`RhoEtaPhiVectorF`** vector based on `rho,eta,phi`
+ coordinates (cylindrical using `eta` instead of `z`) in float
+ precision
#### Constructors and Assignment
The following declarations are available:
``` {.cpp}
- // create an empty vector (x=0, y=0, z=0, t=0)
- XYZTVector v1;
- // vector with x=1, y=2, z=3, t=4
- XYZTVector v2(1,2,3,4);
- // vector with pt=1, eta=2, phi=PI, E=5
- PtEtaPhiEVector v3(1,2,PI,5);
+XYZVector v1; //an empty vector (x=0, y=0, z=0)
+XYZVector v2(1,2,3); //vector with x=1, y=2, z=3;
+Polar3DVector v3(1,PI/2,PI); //vector with r=1, theta=PI/2, phi=PI
+RhoEtaPHiVector v4(1,2, PI); //vector with rho=1, eta=2, phi=PI
```
-Note that each type of vector is constructed by passing its coordinate
-representation, so a **`XYZTVector`**`(1,2,3,4)` is different from a
-`PtEtaPhiEVector(1,2,3,4)`. In addition, the Vector classes can be
-constructed by any vector, which implements the accessors `x()`, `y()`,
-`z()` and `t()`.
-
-This can be another `ROOT::Math::`**`LorentzVector`** based on a
-different coordinate system or any vector of a different package, like
-the CLHEP **`HepLorentzVector`** that implements the required signature.
+Note that each vector type is constructed by passing its coordinate
+representation, so a `XYZVector(1,2,3)` is different from a
+`Polar3DVector(1,2,3)`. In addition, the vector classes can be
+constructed by any vector, which implements the accessors `x()`, `y()`
+and `z()`. This can be another 3D vector based on a different coordinate
+system type. It can be even any vector of a different package, like the
+CLHEP **`HepThreeVector`** that implements the required signature.
``` {.cpp}
- XYZTVector v1(1,2,3,4);
- PtEtaPhiEVector v2(v1);
- CLHEP::HepLorentzVector q(1,2,3,4);
- XYZTVector v3(q);
+ XYZVector v1(1,2,3);
+ RhoEtaPhiVector r2(v1);
+ CLHEP::HepThreeVector q(1,2,3);
+ XYZVector v3(q);
```
#### Coordinate Accessors
-All the same coordinate accessors are available through the interface of
-`ROOT::Math::`**`LorentzVector`**. For example:
+All coordinate accessors are available through the class
+`ROOT::Math::`**`DisplacementVector3D`**:
``` {.cpp}
//returns cartesian components for the cartesian vector v1
- v1.X(); v1.X(); v1.Z(); v1.T();
- //returns cartesian components for the cylindrical vector v2
- v2.Px(); v2.Py(); v2.Pz(); v2.E();
- //returns other components for the cartesian vector v1
- v1.Pt(); v1.Eta(); v1.Phi(); v1.M()
+ v1.X(); v1.Y(); v1.Z();
+ //returns cylindrical components for the cartesian vector v1
+ v1.Rho(); v1.Eta(); v1.Phi();
+ //returns cartesian components for the cylindrical vector r2
+ r2.X(); r2.Y(); r2.Z()
```
-In addition, all 4 vector coordinates can be retrieved with the
-`GetCoordinates` method:
+In addition, all the 3 coordinates of the vector can be retrieved with
+the `GetCoordinates` method:
``` {.cpp}
- double d[4];
- //fill d array with (x,y,z,t) components of v1
+ double d[3];
+ //fill d array with (x,y,z) components of v1
v1.GetCoordinates(d);
- //fill d array with (pt,eta,phi,e) components of v2
- v2.GetCoordinates(d);
- std::vector w(4);
- //fill std::vector with (x,y,z,t)
- v1.GetCoordinates(w.begin(),w.end());
- //components of v1
+ //fill d array with (r,eta,phi) components of r2
+ r2.GetCoordinates(d);
+ std::vector vc(3);
+ //fill std::vector with (x,y,z) components of v1
+ v1.GetCoordinates(vc.begin(),vc.end());
```
-To get information on all the coordinate accessors see the
-`ROOT::Math::`**`LorentzVector`** reference documentation.
+See the reference documentation of
+`ROOT::Math::`**`DisplacementVector3D`** for more details on all the
+coordinate accessors.
#### Setter Methods
One can set only all the three coordinates via:
``` {.cpp}
- //sets the (x,y,z,t) for a XYZTVector
- v1.SetCoordinates(c1,c2,c3,c4);
- //sets pt,eta,phi,e for a PtEtaPhiEVector
- v2.SetCoordinates(c1,c2,c3,c4);
- //sets cartesian components for PtEtaPhiEVector
- v2.SetXYZ(x,y,z,t);
+ v1.SetCoordinates(c1,c2,c3); // (x,y,z) for a XYZVector
+ r2.SetCoordinates(c1,c2,c3); // r,theta,phi for a Polar3DVector
+ r2.SetXYZ(x,y,z); // 3 cartesian components for Polar3DVector
```
Single coordinate setter methods are available for the basic vector
-coordinates, like `SetX()` for a `XYZTVector` or `SetPt()` for a
-**`PtEtaPhiEVector`**. Attempting to do a `SetX()` on a non-Cartesian
-vector will not compile.
+coordinates, like `SetX()` for a **`XYZVector`** or `SetR()` for a polar
+vector. Attempting to do a `SetX()` on a polar vector will not compile.
``` {.cpp}
- XYZTVector v1;
- v1.SetX(1); //OK setting x for a cartesian vector
- PtEtaPhiEVector v2;
- v2.SetX(1); //ERROR: cannot set X for a non-cartesian
- //vector. Method will not compile.
- v2.SetR(1) // OK setting Pt for a PtEtaPhiEVector vector
-```
-
-In addition, there are setter methods from C arrays or iterators.
-
+ XYZVector v1;
+ v1.SetX(1); //OK setting x for a Cartesian vector
+ Polar3DVector v2;
+ v2.SetX(1); //ERROR: cannot set X for a Polar vector.
+ //Method will not compile
+ v2.SetR(1); //OK setting r for a Polar vector
+```
+
+In addition, there are setter methods from C arrays or iterator
+
``` {.cpp}
- double d[4] = {1.,2.,3.,4.};
- XYZTVector v;
- //set (x,y,z,t) components of v using values from d
+ double d[3] = {1.,2.,3.};
+ XYZVector v;
+ // set (x,y,z) components of v using values from d
v.SetCoordinates(d);
```
-or for example from an `std::vector `using the iterators
+or, for example, from an `std::vector` using the iterator
``` {.cpp}
- std::vector w(4);
- //set (x,y,z,t) components of v using values from w
+ std::vector w(3);
+ // set (x,y,z) components of v using values from w
v.SetCoordinates(w.begin(),w.end());
```
#### Arithmetic Operations
-The following operations are possible between Lorentz vectors classes,
-even of different coordinate system types: (`v` and` w` are two Lorentz
-vector of the same type, `q `is a generic Lorentz vector implementing
-`x()`, `y()`, `z()` and `t()`, and `a` is a generic scalar type: double,
-float, int, etc.) .
+The following operations are possible between vector classes, even of
+different coordinate system types: (`v1,v2` are any type of
+**`ROOT::Math::DisplacementVector3D`** classes, `v3` is the same type
+of `v1`; `a` is a scalar value)
``` {.cpp}
- v += q;
- v -= q;
- v = -q;
- v *= a;
- v /= a;
- w = v + q;
- w = v - q;
- w = v * a;
- w = a * v;
- w = v / a;
+ v1 += v2;
+ v1 -= v2;
+ v1 = - v2;
+ v1 *= a;
+ v1 /= a;
+ v2 = a * v1;
+ v2 = v1 / a;
+ v2 = v1 * a;
+ v3 = v1 + v2;
+ v3 = v1 - v2;
```
#### Comparison
+For `v1` and `v2` of the same type (same coordinate system and same
+scalar type):
+
``` {.cpp}
- v == w;
- v != w;
+ v1 == v2;
+ v1 != v2;
```
-#### Other Methods
+#### Dot and Cross Product
+
+We support the dot and cross products, through the `Dot()` and `Cross()`
+method, with any vector (`q`) implementing `x()`, `y()` and `z()`.
``` {.cpp}
- a = v.Dot(q); //dot product in metric(+,+,+,-) of 2 LorentzVectors
- XYZVector s = v.Vect() //return the spatial components (x,y,z)
- v.Beta(); //return beta and gamma value (vector must
- v.Gamma() // be time-like otherwise result is meaningless)
- XYZVector b = v.BoostToCM(); //return boost vector which will bring
- //the Vector in its mas frame (P=0)
+ XYZVector v1(x,y,z);
+ double s = v1.Dot(q);
+ XYZVector v2 = v1.Cross(q);
```
-### Example: Vector Transformations
+Note that the multiplication between two vectors using the operator `*`
+is not supported because it is ambiguous.
+#### Other Methods
-Transformation classes are grouped in rotations (in three dimensions),
-Lorentz transformations and Poincarre transformations, which are
-translation`/`rotation combinations. Each group has several members
-which may model physically equivalent transformations but with
-different internal representations. All the classes are non-template
-and use double precision as the scalar type. The following types of
-transformation classes are defined:
+``` {.cpp}
+ XYZVector u = v1.Unit(); //return unit vector parallel to v1
+```
-3D rotations:
+### Example: 3D Point Classes
-- `ROOT::Math::`**`Rotation3D`**, rotation described by a 3x3 matrix
- of doubles
-- `ROOT::Math::`**`EulerAngles`** rotation described by the three
- Euler angles (`phi`, `theta` and `psi`) following the
- [GoldStein definition](GoldStein definition).
+To use all possible types of 3D points one must include the header file
+`Math/Point3D.h`. The following typedef's defined in the header file
+`Math/Point3Dfwd.h`, are available for different instantiations of the
+template class **`ROOT::Math`**`::`**`PositionVector3D`**:
-- **`ROOT::Math::RotationZYX`** rotation described by three angles
- defining a rotation first along the `Z` axis, then along the rotated
- `Y'` axis and then along the rotated `X''` axis.
+- `ROOT::Math::`**`XYZPoint`** point based on `x`, `y`, `z`
+ coordinates (Cartesian) in double precision
-- `ROOT::Math::`**`AxisAngle`**, rotation described by a vector (axis)
- and an angle
+- `ROOT::Math::`**`XYZPointF`** point based on `x`, `y`, `z`
+ coordinates (Cartesian) in float precision
-- `ROOT::Math::`**`Quaternion`**, rotation described by a quaternion
- (4 numbers)
+- `ROOT::Math::`**`Polar3DPoint`** point based on `r`, `theta`, `phi`
+ coordinates (polar) in double precision
-- `ROOT::Math::`**`RotationX`**, specialized rotation along the X axis
+- `ROOT::Math::`**`Polar3DPointF`** point based on `r`, `theta`, `phi`
+ coordinates (polar) in float precision
-- `ROOT::Math::`**`RotationY`**, specialized rotation along the Y axis
+- `ROOT::Math::`**`RhoZPhiPoint`** point based on `rho`, `z`, `phi`
+ coordinates (cylindrical using `z`) in double precision
-- `ROOT::Math::`**`RotationZ`**, specialized rotation along the Z axis
+- `ROOT::Math::`**`RhoZPhiPointF`** point based on `rho`, `z`, `phi`
+ coordinates (cylindrical using `z`) in float precision
-3D transformations (rotations + translations)
+- `ROOT::Math::`**`RhoEtaPhiPoint`** point based on `rho`, `eta`,
+ `phi` coordinates (cylindrical using eta instead of `z`) in double
+ precision
-- `ROOT::Math::`**`Transform3D`**, (rotations and then translation)
- described by a 3x4 matrix (12 double numbers)
+- `ROOT::Math::`**`RhoEtaPhiPointF`** point based on `rho`, `eta`,
+ `phi` coordinates (cylindrical using eta instead of `z`) in float
+ precision
-- **`ROOT::Math::Translation3D`** (only translation) described by a 3D
- Vector
+#### Constructors and Assignment
-Lorentz rotations and boosts
+The following declarations are available:
-- `ROOT::Math::`**`LorentzRotation`**, 4D rotation (3D rotation plus a
- boost) described by a 4x4 matrix
+``` {.cpp}
+ XYZPoint p1; //an empty vector (x=0, y=0, z=0)
+ XYZPoint p2(1,2,3); //
+```
-- `ROOT::Math::`**`Boost`**, a Lorentz boost in an arbitrary direction
- and described by a 4x4 symmetrix matrix (10 numbers)
+Note that each point type is constructed by passing its coordinate
+representation, so a `XYZPoint(1,2,3)` is different from a
+`Polar3DPoint(1,2,3)`. In addition the point classes can be constructed
+by any vector, which implements the accessors `x()`, `y()` and `z()`.
+This can be another 3D point based on a different coordinate system type
+or even any vector of a different package, like the CLHEP
+**`HepThreePoint`** that implements the required signatures.
-- `ROOT::Math::`**`BoostX`**, a boost in the X axis direction
+``` {.cpp}
+ XYZPoint p1(1,2,3);
+ RhoEtaPHiPoint r2(v1);
+ CLHEP::HepThreePoint q(1,2,3);
+ XYZPoint p3(q);
+```
-- `ROOT::Math::`**`BoostY`**, a boost in the Y axis direction
+#### Coordinate Accessors and Setter Methods
-- `ROOT::Math::`**`BoostZ`**, a boost in the Z axis direction
+For the points classes we have the same getter and setter methods as for
+the vector classes. See "Example: 3D Vector Classes".
-#### Constructors
+#### Point-Vector Operations
-All rotations and transformations are default constructible (giving the
-identity transformation). All rotations are constructible taking a
-number of scalar arguments matching the number (and order of
-components).
+The following operations are possible between points and vector classes:
+(`p1`, `p2` and `p3` are instantiations of the
+`ROOT::Math::`**`PositionVector3D`** objects with `p1` and `p3` of the
+same type; `v1` and `v2` are `ROOT::Math::`**`DisplacementVector3D`**
+objects).
``` {.cpp}
- Rotation3D rI; //a summy rotation (Identity matrix)
- RotationX rX(PI); //a RotationX with an angle PI
- EulerAngles rE(phi,theta,psi); //an Euler rotation with phi,
- //theta,psi angles
- XYZVector u(ux,uy,uz);
- AxisAngle rA(u,delta); //a rotation based on direction u,
- //angle delta
+ p1 += v1;
+ p1 -= v1;
+ p3 = p1 + v1; // p1 and p3 are the same type
+ p3 = v1 + p1; // p3 is based on the same coordinate system as v1
+ p3 = p1 - v1;
+ p3 = v1 - p1;
+ v2 = p1 - p2; // difference between points returns a vector v2
+ // based on the same coordinate system as p1
```
-In addition, all rotations and transformations (other than the axial
-rotations) and transformations are constructible from (`begin`,`end`)
-iterators or from pointers behave like iterators.
+Note that the addition between two points is **NOT** possible and the
+difference between points returns a vector.
-``` {.cpp}
- double data[9];
- //create a rotation from a rotation matrix
- Rotation3D r(data,data+9);
- std::vector w(12);
- //create Transform3D from std::vector content
- Transform3D t(w.begin(),w.end());
-```
+#### Other Operations
-All rotations, except the axial rotations, are constructible and
-assigned from any other type of rotation (including the axial):
+Exactly as for the 3D Vectors, the following operations are allowed:
-``` {.cpp}
- //create a rotation 3D from a rotation along X axis of angle PI
- Rotation3D r(ROOT::Math::RotationX(PI));
+- comparison of points
- //construct an Euler rotation from A Rotation3D
- EulerAngles r2(r);
+- scaling and division of points with a scalar
- //assign an Axis rotation from an Euler Rotation
- AxisAngle r3; r3 = r2;
-```
+- dot and cross product with any type of vector
-**`Transform3D`** (rotation + translation) can be constructed from a
-rotation and a translation vector:
+### Example: LorentzVector Classes
-``` {.cpp}
- Rotation3D r;
- XYZVector v;
- Transform3D t1(r,v); //construct from rotation and then
- //translation
- Transform3D t2(v,r); //construct inverse from first translation
- //then rotation
- Transform3D t3(r); //construct from only a rotation
- //(zero translation)
- Transform3D t4(v); //construct from only translation
- //(identity rotation)
-```
-#### Operations
+As in the 3D case, typedef's are defined for user convenience. and can
+be used by including the header file `Math/Vector4D.h`. The following
+typedef's, defined in the header file `Math/Vector4Dfwd.h`, are
+available for the different instantiations of the template class
+**`ROOT::Math::LorentzVector`**:
-All transformations can be applied to vector and points using the
-operator \* or using the operator()
+- `ROOT::Math::`**`XYZTVector`** vector based on `x`, `y`, `z`, `t`
+ coordinates (Cartesian) in double precision
-``` {.cpp}
- XYZVector v1(...);
- Rotation3D r(...);
- XYZVector v2 = r*v1; //rotate vector v1 using r
- v2 = r(v1); //equivalent
-```
+- `ROOT::Math::`**`XYZTVectorF`** vector based on `x`, `y`, `z`, `t`
+ coordinates (Cartesian) in float precision
-Transformations can be combined using the operator `*`. Rotation,
-translation and **`Transform3D`** classes can be all combined with the
-operator \*. The result of a combination of a rotation and a translation
-will be a **`Transform3D`** class. Note that the rotations are not
-commutative, the order is then important.
+- `ROOT::Math::`**`PtEtaPhiEVector`** vector based on `pt(rho)`,
+ `eta`, `phi` and `E(t)` coordinates in double precision
-``` {.cpp}
- Rotation3D r1(...);
- Rotation3D r2(...);
- Rotation3D r3 = r2*r1; //a combine rotation r3 by
- //applying first r1 then r2
-```
+- `ROOT::Math::`**`PtEtaPhiMVector`** vector based on `pt(rho)`,
+ `eta`, `phi` and `M(t)` coordinates in double precision
-We can combine rotations of different types, like `Rotation3D` with any
-other type of rotations. The product of two different axial rotations
-returns a `Rotation3D`:
+- `ROOT::Math::`**`PxPyPzMVector`** vector based on `px`, `py`, `pz`
+ and `M(mass)` coordinates in double precision
+
+The metric used for all the LorentzVector is (`-,-,-,+`) .
+
+#### Constructors and Assignment
+
+The following declarations are available:
``` {.cpp}
- RotationX rx(1.);
- RotationY ry(2.);
- Rotation3D r = ry * rx; //rotation along X and then Y axis
+ // create an empty vector (x=0, y=0, z=0, t=0)
+ XYZTVector v1;
+ // vector with x=1, y=2, z=3, t=4
+ XYZTVector v2(1,2,3,4);
+ // vector with pt=1, eta=2, phi=PI, E=5
+ PtEtaPhiEVector v3(1,2,PI,5);
```
-It is also possible to invert all the transformation or return their
-inverse:
+Note that each type of vector is constructed by passing its coordinate
+representation, so a **`XYZTVector`**`(1,2,3,4)` is different from a
+`PtEtaPhiEVector(1,2,3,4)`. In addition, the Vector classes can be
+constructed by any vector, which implements the accessors `x()`, `y()`,
+`z()` and `t()`.
+
+This can be another `ROOT::Math::`**`LorentzVector`** based on a
+different coordinate system or any vector of a different package, like
+the CLHEP **`HepLorentzVector`** that implements the required signature.
``` {.cpp}
- Rotation3D r1(...);
- r1.Invert(); //invert the rotation modifying its content
- Rotation3D r2 =r1.Inverse(); //return the inverse in a new
- //rotation class
+ XYZTVector v1(1,2,3,4);
+ PtEtaPhiEVector v2(v1);
+ CLHEP::HepLorentzVector q(1,2,3,4);
+ XYZTVector v3(q);
```
-We have used rotation as examples, but all these operations can be
-applied to all the transformation classes.
-
-#### Set/GetComponents Methods
+#### Coordinate Accessors
-Common methods to all transformations are `Get` and `SetComponents`.
-They can be used to retrieve all the scalar values on which the
-transformation is based.
+All the same coordinate accessors are available through the interface of
+`ROOT::Math::`**`LorentzVector`**. For example:
``` {.cpp}
- RotationX rx;
- rx.SetComponents(1.); //set agle of the X rotation
- double d[9] = {........};
- Rotation3D r;
- r.SetComponents(d,d+9); //set 9 components of 3D rotation
- double d[16];
- LorentzRotation lr;
- lr.GetComponents(d,d+16); //get 16 components of a LorentzRotation
- TMatrixD(3,4) m;
- Transform3D t;
- t.GetComponens(m); //fill 3x4 matrix with components of t
+ //returns cartesian components for the cartesian vector v1
+ v1.X(); v1.X(); v1.Z(); v1.T();
+ //returns cartesian components for the cylindrical vector v2
+ v2.Px(); v2.Py(); v2.Pz(); v2.E();
+ //returns other components for the cartesian vector v1
+ v1.Pt(); v1.Eta(); v1.Phi(); v1.M()
```
-The` GetComponents` and `SetComponents` methods can be used with a
-signature based iterators or by using any foreign matrix which
-implements the `operator(i,j)` or a different signatures depending on
-the transformation type. For more details on all methods see the
-reference documentation of any specific transformation class.
+In addition, all 4 vector coordinates can be retrieved with the
+`GetCoordinates` method:
-### Example with External Packages
+``` {.cpp}
+ double d[4];
+ //fill d array with (x,y,z,t) components of v1
+ v1.GetCoordinates(d);
+ //fill d array with (pt,eta,phi,e) components of v2
+ v2.GetCoordinates(d);
+ std::vector w(4);
+ //fill std::vector with (x,y,z,t)
+ v1.GetCoordinates(w.begin(),w.end());
+ //components of v1
+```
+To get information on all the coordinate accessors see the
+`ROOT::Math::`**`LorentzVector`** reference documentation.
-#### Connection to Linear Algebra Classes
+#### Setter Methods
-It is possible to use the vector and rotation classes together with the
-linear algebra classes and to set and get the contents of any 3D or 4D
-vector from a linear algebra vector class which implements an iterator
-or something which behaves like an iterator. For example a pointer to a
-C array (double`*`) behaves like an iterator. It is then assumed that
-the coordinates, like (`x,y,z`) will be stored contiguously.
+One can set only all the three coordinates via:
``` {.cpp}
- TVectorD r2(N); //ROOT Linear Algebra Vector containing
- //many vectors
- XYZVector v2;
- //construct vector from x=r[INDEX], y=r[INDEX+1], z=r[INDEX+2]
- v2.SetCoordinates(&r2[INDEX],&r2[index]+3);
+ //sets the (x,y,z,t) for a XYZTVector
+ v1.SetCoordinates(c1,c2,c3,c4);
+ //sets pt,eta,phi,e for a PtEtaPhiEVector
+ v2.SetCoordinates(c1,c2,c3,c4);
+ //sets cartesian components for PtEtaPhiEVector
+ v2.SetXYZ(x,y,z,t);
```
-To fill a linear algebra vector from a 3D or 4D vector, with
-`GetCoordinates()` one can get the internal coordinate data.
+Single coordinate setter methods are available for the basic vector
+coordinates, like `SetX()` for a `XYZTVector` or `SetPt()` for a
+**`PtEtaPhiEVector`**. Attempting to do a `SetX()` on a non-Cartesian
+vector will not compile.
``` {.cpp}
- HepVector c(3); //CLHEP Linear algebra vector
- //fill HepVector c with c[0]=x, c[1]=y, c[2]=z
- v2.GetCoordinates(&c[0],&c[index]+3)
+ XYZTVector v1;
+ v1.SetX(1); //OK setting x for a cartesian vector
+ PtEtaPhiEVector v2;
+ v2.SetX(1); //ERROR: cannot set X for a non-cartesian
+ //vector. Method will not compile.
+ v2.SetR(1) // OK setting Pt for a PtEtaPhiEVector vector
```
-or using **`TVectorD`**:
+In addition, there are setter methods from C arrays or iterators.
``` {.cpp}
- double *data[3];
- v2.GetCoordinates(data,data+3);
- TVectorD r1(3,data); //create a new Linear Algebra vector
- //copying the data
+ double d[4] = {1.,2.,3.,4.};
+ XYZTVector v;
+ //set (x,y,z,t) components of v using values from d
+ v.SetCoordinates(d);
```
-In the case of transformations, constructor and method to set`/`get
-components exist with linear algebra matrices. The requisite is that the
-matrix data are stored, for example in the case of a Lorentz rotation,
-from (`0,0`) thru (`3,3`)
+or for example from an `std::vector `using the iterators
``` {.cpp}
- TMatrixD(4,4) m;
- LorentzRotation r(m); //create Lorentz r
+ std::vector w(4);
+ //set (x,y,z,t) components of v using values from w
+ v.SetCoordinates(w.begin(),w.end());
```
-#### Connection to Other Vector Classes
+#### Arithmetic Operations
-The 3D and 4D vectors of the `GenVector` package can be constructed and
-assigned from any vector which satisfies the following requisites:
+The following operations are possible between Lorentz vectors classes,
+even of different coordinate system types: (`v` and` w` are two Lorentz
+vector of the same type, `q `is a generic Lorentz vector implementing
+`x()`, `y()`, `z()` and `t()`, and `a` is a generic scalar type: double,
+float, int, etc.) .
-- for 3D vectors implementing the `x()`, `y()` and `z()` methods
+``` {.cpp}
+ v += q;
+ v -= q;
+ v = -q;
+ v *= a;
+ v /= a;
+ w = v + q;
+ w = v - q;
+ w = v * a;
+ w = a * v;
+ w = v / a;
+```
-- for Lorentz vectors implementing the `x()`, `y()`, `z()` and `t()`
- methods.
+#### Comparison
``` {.cpp}
- CLHEP::Hep3Vector hv;
- XYZVector v1(hv); //create 3D vector from
- //CLHEP 3D Vector
- HepGeom::Point3D hp;
- XYZPoint p1(hp); //create a 3D p
+ v == w;
+ v != w;
```
-## MathMore Library
+#### Other Methods
+``` {.cpp}
+ a = v.Dot(q); //dot product in metric(+,+,+,-) of 2 LorentzVectors
+ XYZVector s = v.Vect() //return the spatial components (x,y,z)
+ v.Beta(); //return beta and gamma value (vector must
+ v.Gamma() // be time-like otherwise result is meaningless)
+ XYZVector b = v.BoostToCM(); //return boost vector which will bring
+ //the Vector in its mas frame (P=0)
+```
-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 current set
-includes:
+### Example: Vector Transformations
-- Special functions (see Special Functions in MathMore)
-- Mathematical functions used in statistics such as probability density
-functions, cumulative distributions functions and their inverse.
+Transformation classes are grouped in rotations (in three dimensions),
+Lorentz transformations and Poincarre transformations, which are
+translation`/`rotation combinations. Each group has several members
+which may model physically equivalent transformations but with
+different internal representations. All the classes are non-template
+and use double precision as the scalar type. The following types of
+transformation classes are defined:
-- Numerical algorithms for one dimensional functions based on
-implementation of the GNU Scientific Library (GSL):
+3D rotations:
-- Numerical integration using the class **`ROOT::Math::Integrator`**
- which is based on the Adaptive integration algorithms of QUADPACK
+- `ROOT::Math::`**`Rotation3D`**, rotation described by a 3x3 matrix
+ of doubles
-- Numerical differentiation via **`ROOT::Math::Derivator`**
+- `ROOT::Math::`**`EulerAngles`** rotation described by the three
+ Euler angles (`phi`, `theta` and `psi`) following the
+ [GoldStein definition](GoldStein definition).
-- Root finder via **`ROOT::Math::RootFinder`** which uses different
- solver algorithms from GSL
+- **`ROOT::Math::RotationZYX`** rotation described by three angles
+ defining a rotation first along the `Z` axis, then along the rotated
+ `Y'` axis and then along the rotated `X''` axis.
-- Minimization via **`ROOT::Math::Minimizer1D`**
+- `ROOT::Math::`**`AxisAngle`**, rotation described by a vector (axis)
+ and an angle
-- Interpolation via **`ROOT::Math::Interpolation`**. All the GSL
- interpolation types are supported
+- `ROOT::Math::`**`Quaternion`**, rotation described by a quaternion
+ (4 numbers)
-- Function approximation based on Chebyshev polynomials via the class
- **`ROOT::Math::Chebyshev`**
+- `ROOT::Math::`**`RotationX`**, specialized rotation along the X axis
-- Random number generators and distributions
+- `ROOT::Math::`**`RotationY`**, specialized rotation along the Y axis
-- Polynomial evaluation and root solvers
+- `ROOT::Math::`**`RotationZ`**, specialized rotation along the Z axis
-The mathematical functions are implemented as a set of free functions in
-the namespace **`ROOT::Math`**. The naming used for the special
-functions is the same proposed for the C++ standard (see C++ standard
-extension [proposal document](proposal document)).The `MathCore` library
-is implemented wrapping in C++ the GNU Scientific Library ( <GSL>).
-Building `MathMore` requires a version of GSL larger or equal 1.8. The
-source code of `MathMore` is distributed under the GNU General Public
-License.
+3D transformations (rotations + translations)
-`MathMore` (and its ROOT Cling dictionary) can be built within ROOT
-whenever a GSL library is found in the system. The GSL library and
-header file location can be specified in the ROOT configure script, by
-doing:
+- `ROOT::Math::`**`Transform3D`**, (rotations and then translation)
+ described by a 3x4 matrix (12 double numbers)
-```
-./configure --with-gsl-incdir=... --with-gsl-libdir=...
-```
+- **`ROOT::Math::Translation3D`** (only translation) described by a 3D
+ Vector
-`MathMore` can be built also a stand-alone library (without requiring
-ROOT) downloding the tar file from the Web at this link. In this case
-the library will not contain the dictionary information and therefore
-cannot be used interactively
+Lorentz rotations and boosts
-More information on the classes and functions present in `MathMore` is
-available in the
-[online reference documentation](online reference documentation).
+- `ROOT::Math::`**`LorentzRotation`**, 4D rotation (3D rotation plus a
+ boost) described by a 4x4 matrix
-## Mathematical Functions
+- `ROOT::Math::`**`Boost`**, a Lorentz boost in an arbitrary direction
+ and described by a 4x4 symmetrix matrix (10 numbers)
+- `ROOT::Math::`**`BoostX`**, a boost in the X axis direction
-The mathematical functions are present in both `MathCore` and `MathMore`
-libraries. All mathematical functions are implemented as free functions
-in the namespace **`ROOT::Math`**. The most used functions are in the
-`MathCore` library while the others are in the `MathMore` library. The
-functions in `MathMore` are all using the implementation of the GNU
-Scientific Library (GSL). The naming of the special functions is the
-same defined in the C++
-[Technical Report on Standard Library extensions](Technical Report on
-Standard Library extensions).
-The special functions are defined in the header file `Math/SpecFunc.h`.
+- `ROOT::Math::`**`BoostY`**, a boost in the Y axis direction
-### Special Functions in MathCore
+- `ROOT::Math::`**`BoostZ`**, a boost in the Z axis direction
+#### Constructors
-- `ROOT::Math::beta(double x,double y) - `evaluates the beta function:
- $$B(x,y) = \frac{\Gamma(x) \Gamma(y)}{\Gamma(x+y)}$$
+All rotations and transformations are default constructible (giving the
+identity transformation). All rotations are constructible taking a
+number of scalar arguments matching the number (and order of
+components).
-- `double ROOT::Math::erf(double x)` - evaluates the error function
- encountered in integrating the normal
- distribution:
- $$erf(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt$$
+``` {.cpp}
+ Rotation3D rI; //a summy rotation (Identity matrix)
+ RotationX rX(PI); //a RotationX with an angle PI
+ EulerAngles rE(phi,theta,psi); //an Euler rotation with phi,
+ //theta,psi angles
+ XYZVector u(ux,uy,uz);
+ AxisAngle rA(u,delta); //a rotation based on direction u,
+ //angle delta
+```
-- `double ROOT::Math::erfc(double x)` - evaluates the complementary
- error function:
- $$erfc(x) = 1 - erf(x) = \frac{2}{\sqrt{\pi}} \int_{x}^{\infty} e^{-t^2} dt$$
+In addition, all rotations and transformations (other than the axial
+rotations) and transformations are constructible from (`begin`,`end`)
+iterators or from pointers behave like iterators.
-- `double ROOT::Math::tgamma(double x)` - calculates the gamma
- function:
- $$\Gamma(x) = \int_{0}^{\infty} t^{x-1} e^{-t} dt$$
+``` {.cpp}
+ double data[9];
+ //create a rotation from a rotation matrix
+ Rotation3D r(data,data+9);
+ std::vector w(12);
+ //create Transform3D from std::vector content
+ Transform3D t(w.begin(),w.end());
+```
-### Special Functions in MathMore
+All rotations, except the axial rotations, are constructible and
+assigned from any other type of rotation (including the axial):
+``` {.cpp}
+ //create a rotation 3D from a rotation along X axis of angle PI
+ Rotation3D r(ROOT::Math::RotationX(PI));
-- `double ROOT::Math::assoc_legendre(unsigned l,unsigned m,double x) -`computes
- the associated Legendre polynomials (with `m>=0`, `l>=m` and
- `|x|<1)`:
- $$P_{l}^{m}(x) = (1-x^2)^{m/2} \frac{d^m}{dx^m} P_{l}(x)$$
+ //construct an Euler rotation from A Rotation3D
+ EulerAngles r2(r);
-- `double ROOT::Math::comp_ellint_1(double k)` - calculates the
- complete elliptic integral of the first kind (with $0 \le k^2 \le 1$:
- $$
- K(k) = F(k, \pi / 2) = \int_{0}^{\pi /2} \frac{d \theta}{\sqrt{1 - k^2 \sin^2{\theta}}}
- $$
+ //assign an Axis rotation from an Euler Rotation
+ AxisAngle r3; r3 = r2;
+```
-- `double ROOT::Math::comp_ellint_2(double k)` - calculates the
- complete elliptic integral of the second kind (with $0 \le k^2 \le 1$):
- $$
- E(k) = E(k , \pi / 2) = \int_{0}^{\pi /2} \sqrt{1 - k^2 \sin^2{\theta}} d \theta
- $$
+**`Transform3D`** (rotation + translation) can be constructed from a
+rotation and a translation vector:
-- `double ROOT::Math::comp_ellint_3(double n,double k)` - calculates
- the complete elliptic integral of the third kind (with $0 \le k^2 \le 1$):
- $$
- \Pi (n, k, \pi / 2) = \int_{0}^{\pi /2} \frac{d \theta}{(1 - n \sin^2{\theta})\sqrt{1 - k^2 \sin^2{\theta}}}
- $$
+``` {.cpp}
+ Rotation3D r;
+ XYZVector v;
+ Transform3D t1(r,v); //construct from rotation and then
+ //translation
+ Transform3D t2(v,r); //construct inverse from first translation
+ //then rotation
+ Transform3D t3(r); //construct from only a rotation
+ //(zero translation)
+ Transform3D t4(v); //construct from only translation
+ //(identity rotation)
+```
-- `double ROOT::Math::conf_hyperg(double a,double b,double z)` -
- calculates the confluent hyper-geometric functions of the first
- kind:
- $$
- _{1}F_{1}(a;b;z) = \frac{\Gamma(b)}{\Gamma(a)} \sum_{n=0}^{\infty} \frac{\Gamma(a+n)}{\Gamma(b+n)} \frac{z^n}{n!}
- $$
+#### Operations
-- `double ROOT::Math::conf_hypergU(double a,double b,double z)` -
- calculates the confluent hyper-geometric
- functions of the second kind, known also as Kummer function of the second type. It is
- related to the confluent hyper-geometric function of the first kind:
- $$
- U(a,b,z) = \frac{ \pi}{ \sin{\pi b } } \left[ \frac{ _{1}F_{1}(a,b,z) } {\Gamma(a-b+1) } - \frac{ z^{1-b} { _{1}F_{1}}(a-b+1,2-b,z)}{\Gamma(a)} \right]
- $$
+All transformations can be applied to vector and points using the
+operator \* or using the operator()
-- `double ROOT::Math::cyl_bessel_i(double nu,double x)` - calculates
- the modified Bessel function of the first kind, also called regular
- modified (cylindrical) Bessel function:
- $$
- I_{\nu} (x) = i^{-\nu} J_{\nu}(ix) = \sum_{k=0}^{\infty} \frac{(\frac{1}{2}x)^{\nu + 2k}}{k! \Gamma(\nu + k + 1)}
- $$
+``` {.cpp}
+ XYZVector v1(...);
+ Rotation3D r(...);
+ XYZVector v2 = r*v1; //rotate vector v1 using r
+ v2 = r(v1); //equivalent
+```
-- `double ROOT::Math::cyl_bessel_j(double nu,double x)` - calculates
- the (cylindrical) Bessel function of the first kind, also called
- regular (cylindrical) Bessel function:
- $$
- J_{\nu} (x) = \sum_{k=0}^{\infty} \frac{(-1)^k(\frac{1}{2}x)^{\nu + 2k}}{k! \Gamma(\nu + k + 1)}
- $$
+Transformations can be combined using the operator `*`. Rotation,
+translation and **`Transform3D`** classes can be all combined with the
+operator \*. The result of a combination of a rotation and a translation
+will be a **`Transform3D`** class. Note that the rotations are not
+commutative, the order is then important.
-- `double ROOT::Math::cyl_bessel_k(double nu,double x)` - calculates
- the modified Bessel function of the second kind, also called
- irregular modified (cylindrical) Bessel function for $x > 0$, $v > 0$:
- $$
- K_{\nu} (x) = \frac{\pi}{2} i^{\nu + 1} (J_{\nu} (ix) + iN(ix)) = \left\{ \begin{array}{cl} \frac{\pi}{2} \frac{I_{-\nu}(x) - I_{\nu}(x)}{\sin{\nu \pi}} & \mbox{for non-integral $\nu$} \\ \frac{\pi}{2} \lim{\mu \to \nu} \frac{I_{-\mu}(x) - I_{\mu}(x)}{\sin{\mu \pi}} & \mbox{for integral $\nu$} \end{array} \right.
- $$
+``` {.cpp}
+ Rotation3D r1(...);
+ Rotation3D r2(...);
+ Rotation3D r3 = r2*r1; //a combine rotation r3 by
+ //applying first r1 then r2
+```
-- `double ROOT::Math::cyl_neumann(double nu,double x)` - calculates
- the (cylindrical) Bessel function of the second kind, also called
- irregular (cylindrical) Bessel function or (cylindrical) Neumann
- function:
- $$
- N_{\nu} (x) = Y_{\nu} (x) = \left\{ \begin{array}{cl} \frac{J_{\nu} \cos{\nu \pi}-J_{-\nu}(x)}{\sin{\nu \pi}} & \mbox{for non-integral $\nu$} \\ \lim{\mu \to \nu} \frac{J_{\mu} \cos{\mu \pi}-J_{-\mu}(x)}{\sin{\mu \pi}} & \mbox{for integral $\nu$} \end{array} \right.
- $$
+We can combine rotations of different types, like `Rotation3D` with any
+other type of rotations. The product of two different axial rotations
+returns a `Rotation3D`:
-- `double ROOT::Math::ellint_1(double k,double phi)` - calculates
- incomplete elliptic integral of the first kind (with $0 \le k^2 \le 1$):
- $$
- K(k) = F(k, \pi / 2) = \int_{0}^{\pi /2} \frac{d \theta}{\sqrt{1 - k^2 \sin^2{\theta}}}
- $$
+``` {.cpp}
+ RotationX rx(1.);
+ RotationY ry(2.);
+ Rotation3D r = ry * rx; //rotation along X and then Y axis
+```
-- `double ROOT::Math::ellint_2(double k,double phi)` - calculates
- the complete elliptic integral of the second kind (with $0 \le k^2 \le 1$):
- $$
- E(k) = E(k , \pi / 2) = \int_{0}^{\pi /2} \sqrt{1 - k^2 \sin^2{\theta}} d \theta
- $$
+It is also possible to invert all the transformation or return their
+inverse:
-- `double ROOT::Math::ellint_3(double n,double k,double phi)` - calculates
- the complete elliptic integral of the third kind (with $0 \le k^2 \le 1$):
- $$
- \Pi (n, k, \pi / 2) = \int_{0}^{\pi /2} \frac{d \theta}{(1 - n \sin^2{\theta})\sqrt{1 - k^2 \sin^2{\theta}}}
- $$
+``` {.cpp}
+ Rotation3D r1(...);
+ r1.Invert(); //invert the rotation modifying its content
+ Rotation3D r2 =r1.Inverse(); //return the inverse in a new
+ //rotation class
+```
-- `double ROOT::Math::expint(double x)` - calculates the exponential
- integral:
- $$
- Ei(x) = - \int_{-x}^{\infty} \frac{e^{-t}}{t} dt
- $$
+We have used rotation as examples, but all these operations can be
+applied to all the transformation classes.
-- `double ROOT::Math::hyperg(double a,double b,double c,double x)` -
- calculates Gauss' hyper-geometric function:
- $$
- _{2}F_{1}(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a) \Gamma(b)} \sum_{n=0}^{\infty} \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)} \frac{x^n}{n!}
- $$
+#### Set/GetComponents Methods
-- `double ROOT::Math::legendre(unsigned l,double x)` - calculates
- the Legendre polynomials for $l \ge 0$, $|x| \le 1$ in the Rodrigues
- representation:
- $$
- P_{l}(x) = \frac{1}{2^l l!} \frac{d^l}{dx^l} (x^2 - 1)^l
- $$
+Common methods to all transformations are `Get` and `SetComponents`.
+They can be used to retrieve all the scalar values on which the
+transformation is based.
-- `double ROOT::Math::riemann_zeta(double x)` - calculates the
- Riemann zeta function:
- $$
- \zeta (x) = \left\{ \begin{array}{cl} \sum_{k=1}^{\infty}k^{-x} & \mbox{for $x > 1$} \\ 2^x \pi^{x-1} \sin{(\frac{1}{2}\pi x)} \Gamma(1-x) \zeta (1-x) & \mbox{for $x < 1$} \end{array} \right.
- $$
+``` {.cpp}
+ RotationX rx;
+ rx.SetComponents(1.); //set agle of the X rotation
+ double d[9] = {........};
+ Rotation3D r;
+ r.SetComponents(d,d+9); //set 9 components of 3D rotation
+ double d[16];
+ LorentzRotation lr;
+ lr.GetComponents(d,d+16); //get 16 components of a LorentzRotation
+ TMatrixD(3,4) m;
+ Transform3D t;
+ t.GetComponens(m); //fill 3x4 matrix with components of t
+```
-- `double ROOT::Math::sph_bessel(unsigned n,double x)` - calculates
- the spherical Bessel functions of the first kind (also called
- regular spherical Bessel functions):
- $$
- j_{n}(x) = \sqrt{\frac{\pi}{2x}} J_{n+1/2}(x)
- $$
+The` GetComponents` and `SetComponents` methods can be used with a
+signature based iterators or by using any foreign matrix which
+implements the `operator(i,j)` or a different signatures depending on
+the transformation type. For more details on all methods see the
+reference documentation of any specific transformation class.
-- `double ROOT::Math::sph_neumann(unsigned n,double x)` - calculates
- the spherical Bessel functions of the second kind (also called
- irregular spherical Bessel functions or spherical Neumann
- functions):
- $$
- n_n(x) = y_n(x) = \sqrt{\frac{\pi}{2x}} N_{n+1/2}(x)
- $$
+### Example with External Packages
-### Probability Density Functions (PDF)
+#### Connection to Linear Algebra Classes
-Probability density functions of various distributions. All the
-functions, apart from the discrete ones, have the extra location
-parameter `x0`, which by default is zero. For example, in the case of a
-gaussian `pdf`, `x0` is the `mean`, `mu`, of the distribution. All the
-probability density functions are defined in the header file
-`Math/DistFunc.h` and are part of the `MathCore` libraries. The
-definition of these functions is documented in the
-[reference doc for statistical functions](reference doc for statistical functions):
+It is possible to use the vector and rotation classes together with the
+linear algebra classes and to set and get the contents of any 3D or 4D
+vector from a linear algebra vector class which implements an iterator
+or something which behaves like an iterator. For example a pointer to a
+C array (double`*`) behaves like an iterator. It is then assumed that
+the coordinates, like (`x,y,z`) will be stored contiguously.
``` {.cpp}
-double ROOT::Math::beta_pdf(double x,double a, double b);
-double ROOT::Math::binomial_pdf(unsigned int k,double p,unsigned int n);
-double ROOT::Math::breitwigner_pdf(double x,double gamma,double x0=0);
-double ROOT::Math::cauchy_pdf(double x,double b=1,double x0=0);
-double ROOT::Math::chisquared_pdf(double x,double r,double x0=0);
-double ROOT::Math::exponential_pdf(double x,double lambda,double x0=0);
-double ROOT::Math::fdistribution_pdf(double x,double n,double m,double x0=0);
-double ROOT::Math::gamma_pdf(double x,double alpha,double theta,double x0=0);
-double ROOT::Math::gaussian_pdf(double x,double sigma,double x0=0);
-double ROOT::Math::landau_pdf(double x,double s,double x0=0);
-double ROOT::Math::lognormal_pdf(double x,double m,double s,double x0=0);
-double ROOT::Math::normal_pdf(double x,double sigma,double x0=0);
-double ROOT::Math::poisson_pdf(unsigned int n,double mu);
-double ROOT::Math::tdistribution_pdf(double x,double r,double x0=0);
-double ROOT::Math::uniform_pdf(double x,double a,double b,double x0=0);
+ TVectorD r2(N); //ROOT Linear Algebra Vector containing
+ //many vectors
+ XYZVector v2;
+ //construct vector from x=r[INDEX], y=r[INDEX+1], z=r[INDEX+2]
+ v2.SetCoordinates(&r2[INDEX],&r2[index]+3);
```
-### Cumulative Distribution Functions (CDF)
+To fill a linear algebra vector from a 3D or 4D vector, with
+`GetCoordinates()` one can get the internal coordinate data.
+``` {.cpp}
+ HepVector c(3); //CLHEP Linear algebra vector
+ //fill HepVector c with c[0]=x, c[1]=y, c[2]=z
+ v2.GetCoordinates(&c[0],&c[index]+3)
+```
-For all the probability density functions, we have the corresponding
-cumulative distribution functions and their complements. The functions
-with extension `_cdf` calculate the lower tail integral of the
-probability density function:
+or using **`TVectorD`**:
-$$
-D(x) = \int_{-\infty}^{x} p(x') dx'
-$$
+``` {.cpp}
+ double *data[3];
+ v2.GetCoordinates(data,data+3);
+ TVectorD r1(3,data); //create a new Linear Algebra vector
+ //copying the data
+```
-while those with the `cdf_c` extension calculate the upper tail of the
-probability density function, so-called in statistics the survival
-function. For example, the function:
+In the case of transformations, constructor and method to set`/`get
+components exist with linear algebra matrices. The requisite is that the
+matrix data are stored, for example in the case of a Lorentz rotation,
+from (`0,0`) thru (`3,3`)
``` {.cpp}
-double ROOT::Math::gaussian_cdf(double x,double sigma,double x0=0);
+ TMatrixD(4,4) m;
+ LorentzRotation r(m); //create Lorentz r
```
-evaluates the lower tail of the Gaussian distribution:
-$$
-D(x) = \int_{-\infty}^{x} {1 \over \sqrt{2 \pi \sigma^2}} e^{-(x'-x_0)^2 / 2\sigma^2} dx'
-$$
-while the function:
+#### Connection to Other Vector Classes
+
+The 3D and 4D vectors of the `GenVector` package can be constructed and
+assigned from any vector which satisfies the following requisites:
+
+- for 3D vectors implementing the `x()`, `y()` and `z()` methods
+
+- for Lorentz vectors implementing the `x()`, `y()`, `z()` and `t()`
+ methods.
``` {.cpp}
-double ROOT::Math::gaussian_cdf_c(double x, double sigma, double x0=0);
+ CLHEP::Hep3Vector hv;
+ XYZVector v1(hv); //create 3D vector from
+ //CLHEP 3D Vector
+ HepGeom::Point3D hp;
+ XYZPoint p1(hp); //create a 3D p
```
-evaluates the upper tail of the Gaussian distribution:
-$$
-D(x) = \int_{x}^{+\infty} {1 \over \sqrt{2 \pi \sigma^2}} e^{-(x'-x_0)^2 / 2\sigma^2} dx'
-$$
-The cumulative distributions functions are defined in the header file
-`Math/ProbFunc.h`. The majority of the CDF's are present in the
-`MathCore`, apart from the `chisquared`, `fdistribution`, `gamma` and
-`tdistribution`, which are in the `MathMore` library.
-
-#### Inverse of the Cumulative Distribution Functions(Quantiles)
-For almost all the cumulative distribution functions (`_cdf`) and their
-complements (`_cdf_c`) present in the library, we provide the inverse
-functions. The inverse of the cumulative distribution function is called
-in statistics quantile function. The functions with the extension
-`_quantile` calculate the inverse of the cumulative distribution
-function (lower tail integral of the probability density function),
-while those with the *`quantile_c`* extension calculate the inverse of
-the complement of the cumulative distribution (upper tail integral). All
-the inverse distributions are in the MathMore library and are defined in
-the header file `Math/ProbFuncInv.h`.
-The following picture illustrates the available statistical functions
-(PDF, CDF and quantiles) in the case of the normal distribution.
-
## Linear Algebra: SMatrix Package
@@ -2135,70 +3095,6 @@ and `SVector` for for **`Double_t`**, **`Float_t`** and **`Double32_t`**
up to dimension 7. This allows the possibility to store them in a ROOT
file.
-## Minuit2 Package
-
-
-`Minuit2` is a new object-oriented implementation, written in C++, of
-the popular `MINUIT` minimization package. Compared with the
-**`TMinuit`** class, which is a direct conversion from FORTRAN to C++,
-`Minuit2` is a complete redesign and re-implementation of the package.
-This new version provides 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" in "Fitting Histograms"
-chapter), which is an optimized method for least square and log
-likelihood minimizations. Minuit2 has been originally developed by M.
-Winkler and F. James in the SEAL project. More information can be found
-on the [MINUIT Web Site](MINUIT Web Site) and in particular at the
-following documentation page at
-<http://www.cern.ch/minuit/doc/doc.html>.
-
-The API has been then changed in this new version to follow the ROOT
-coding convention (function names starting with capital letters) and the
-classes have been moved inside the namespace **`ROOT::Minuit2`**. In
-addition, the ROOT distribution contains classes needed to integrate
-`Minuit2` in the ROOT framework, like **`TFitterMinuit`** and
-**`TFitterFumili`**. `Minuit2` can be used in ROOT as another fitter
-plug-in. For example for using it in histogram fitting, one only needs
-to do:
-
-``` {.cpp}
- TVirtualFitter::SetDefaultFitter("Minuit2"); //or Fumili2 for the
- // FUMILI algorithmhistogram->Fit();
-```
-
-For minimization problem, providing an FCN function to minimize, one can
-do:
-
-``` {.cpp}
- TVirtualFitter::SetDefaultFitter("Minuit2");
- TVirtualFitter * minuit2 = TVirtualFitter::Fitter(0,2);
-```
-
-Then set the parameters, the FCN and minimize using the
-**`TVirtualFitter`** methods: `SetParameter`, `SetFCN` and
-`ExecuteCommand`. The FCN function can also be given to Minuit2 as an
-instance of a class implementing the **`ROOT::Minuit2::FCNBase`**
-interface. In this case one must use directly the **`TFitterMinuit`**
-class via the method `SetMinuitFCN`.
-
-Examples on how to use the `Minuit2` and `Fumili2` plug-ins are provided
-in the tutorials' directory `$ROOTSYS/tutorials/fit`:
-`minuit2FitBench.C`, `minuit2FitBench2D.C` and `minuit2GausFit.C`. More
-information on the classes and functions present in `Minuit2` is
-available at
-[online reference documentation](online reference documentation). In
-addition, the C++ MINUIT User Guide provides all the information needed
-for using directly the package without **`TVirtualFitter`** interface
-(see <http://seal.cern.ch/documents/minuit/mnusersguide.pdf>). Useful
-information on MINUIT and minimization in general is provided in the
-following documents:
-
-F. James, *Minuit Tutorial on Function Minimization* (
-<http://seal.cern.ch/documents/minuit/mntutorial.pdf>); F. James, *The
-Interpretation of Errors in Minuit* (
-<http://seal.cern.ch/documents/minuit/mnerror.pdf>);
## ROOT Statistics Classes
diff --git a/documentation/users-guide/pictures/Integration.png b/documentation/users-guide/pictures/Integration.png
new file mode 100644
index 0000000000000000000000000000000000000000..4cb019a1cc4a4b6e457a332eb8dc8bf4fb7e3663
Binary files /dev/null and b/documentation/users-guide/pictures/Integration.png differ
diff --git a/documentation/users-guide/pictures/Minimization.png b/documentation/users-guide/pictures/Minimization.png
new file mode 100644
index 0000000000000000000000000000000000000000..d3b4a0ac37b20ac621f06b66981357489d82b356
Binary files /dev/null and b/documentation/users-guide/pictures/Minimization.png differ
diff --git a/documentation/users-guide/pictures/function-hierarchy.png b/documentation/users-guide/pictures/function-hierarchy.png
new file mode 100644
index 0000000000000000000000000000000000000000..f1a8c5f482c6b2021d5f6e8705fecfff532101db
Binary files /dev/null and b/documentation/users-guide/pictures/function-hierarchy.png differ