Fitting in ROOT is based on the TMinuit class. See in particular
[Basic Concepts of Minuit](http://root.cern.ch/doc/master/classTMinuit.html#TMinuit:Basic Concepts of Minuit). The user has the choice of using directly the TMinuit class, as illustrated in the last section below or use specialized functions provided in the histogram class TH1 or the graph classes TGraph and TGraphErrors.
The Minuit package acts on a multiparameter Fortran function FCN. In the ROOT implementation, the function FCN is defined via the TMinuit::SetFCN member function when an Histogram.Fit command is invoked. The value of FCN will in general depend on one or more variable parameters.
To take a simple example, in case of ROOT histograms (classes TH1C,TH1S,TH1F,TH1D) the Fit function defines the Minuit fitting function as being H1FitChisquare. H1FitChisquare calculates the chisquare between the user fitting function (gaussian, polynomial, user defined, etc) and the data for given values of the parameters. It is the task of MINUIT to find those values of the parameters which give the lowest value of chisquare.
Fitting 1-Dim histograms
Fitting histograms is done via TH1::Fit. The name of the fitted function (the model) is passed as first parameter. This name may be one of the ROOT pre-defined function names or a user-defined function.
Algorithm to compute the errors per bin
If TH1::Sumw2 has been called for the histogram, the error computed is the square root of the sum of the square of weigths for this bin.
If not, the error computed is the square root of the bin content.
With pre-defined functions
The following functions are automatically created by ROOT when any fitting function is invoked (eg, TH1::FIT):
- gaus. A gaussian with 3 parameters: f(x) = p0exp(-0.5((x-p1)/p2)^2)).
- expo. An exponential with 2 parameters: f(x) = exp(p0+p1*x).
- polN. a polynomial of degree N: f(x) = p0 + p1*x + p2*x^2 +…
As an example, the following command will fit histogram object hist with a gaussian. By default, the fitting function object will be added to the histogram structure and it will be drawn in the current pad. For pre-defined functions, there is no need to set initial values for the parameters. ROOT will do it automatically for you.
Root > hist.Fit("gaus");
With user-defined functions
You can create a TF1 object. The TF1 function may be created out of existing known expressions (see TFormula). For example, the following command creates a function called "myfit" with 3 parameters in the range between 0 and 2.
Root > myfit->SetParName(0,"c0");
Root > myfit->SetParName(1,"c1");
Root > myfit->SetParName(2,"slope");
Root > myfit->SetParameter(0, 1);
Root > myfit->SetParameter(1, 0.05);
Root > myfit->SetParameter(2, 0.2);
We are now ready to fit:
Root > hist->Fit("myfit");
You can also create your own fitting function. This function must have 2 parameters:
- Double_t *v: a pointer to the variables array. This array must be a 1-dim array with v[0] = x in case of a 1-dim histogram, v[0] =x, v[1] = y for a 2-dim histogram, etc.
- Double_t *par: a pointer to the parameters array. par will contain the current values of parameters when it is called by the FCN function.
The following macro fitexample.C illustrates how fit a 1-dim histogram stored in a ROOT file with a user-defined function.
//_____________________macro fitexample.C___________________________
Double_t fitf(Double_t *v, Double_t *par)
{
Double_t arg = 0;
if (par[2] != 0) arg = (v[0] - par[1])/par[2];
Double_t fitval = par[0]*TMath::Exp(-0.5*arg*arg);
return fitval;
}
void fitexample()
{
TFile *f = new TFile("hsimple.root");
TH1F *hpx = (TH1F*)f->Get("hpx");
//create a function with 3 parameters in the range [-3,3]
TF1 *func = new TF1("fit",fitf,-3,3,3);
func->SetParameters(500,hpx->GetMean(),hpx->GetRMS());
func->SetParNames("Constant","Mean_value","Sigma");
hpx->Fit("fit");
}
Fitting a sub-range of the histogram bins
By default,TH1::Fit will fit the function on the defined histogram range. You can specify the option “r” in the second parameter to restrict the fit to the range specified in the TF1object. For more complete examples, see A simple fitting example and Fitting histogram sub-ranges.
Fitting several functions on the same histogram
The example Fitting histogram sub-ranges also illustrates how to fit several functions to the same histogram. By default, a Fit operation deletes the previously fitted function in the histogram data structure. You can specify the option “+” in the second parameter of TH1::Fit to add the newly fitted function to the existing list of fitted functions for this histogram. Note that the fitted function(s) are saved with the histogram when the histogram is written to a ROOT file.
Fitting non-equidistant data points
To fit non-equidistant data points, you can use the services of the classes TGraph or TGraphErrors. Examples of creation of such objects are given in:
- <a href="http://root.cern.ch/doc/master/classtutorials/graphs/graph.C.html">A simple Graph example</a>.
- <a href="http://root.cern.ch/doc/master/classtutorials/graphs/gerrors.C.html">Example of a graph with error bars</a>.
To fit such objects, use TGraph::Fit. For example, to fit the graph called "gr in the above examples, with a polynomial of degree 4, you do:
Root > gr->Fit("pol4");```
Like for histograms above, you can use the interactive Fitpanel option in the graph context menu.
Fitting in general using the class TMinuit
The following example (part of the tutorials: Using Minuit to fit non-equidistant points)) macro illustrates how to fit in general by using the services of the TMinuit class only.
//
// Example of a program to fit non-equidistant data points
// =======================================================
//
// The fitting function fcn is a simple chisquare function
// The data consists of 5 data points (arrays x,y,z) + the errors in errorsz
// More details on the various functions or parameters for these functions
// can be obtained in an interactive ROOT session with:
// Root > TMinuit *minuit = new TMinuit(10);
// Root > minuit->mnhelp("*",0) to see the list of possible keywords
// Root > minuit->mnhelp("SET",0) explains most parameters
//
Float_t z[5],x[5],y[5],errorz[5];
//______________________________________________________________________________
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
const Int_t nbins = 5;
Int_t i;
//calculate chisquare
Double_t chisq = 0;
Double_t delta;
for (i=0;i<nbins; i++) {
delta = (z[i]-func(x[i],y[i],par))/errorz[i];
chisq += delta*delta;
}
f = chisq;
}
//______________________________________________________________________________
Double_t func(float x,float y,Double_t *par)
{
Double_t value=( (par[0]*par[0])/(x*x)-1)/ ( par[1]+par[2]*y-par[3]*y*y);
return value;
}
//______________________________________________________________________________
void Ifit()
{
// The z values
z[0]=1;
z[1]=0.96;
z[2]=0.89;
z[3]=0.85;
z[4]=0.78;
// The errors on z values
Float_t error = 0.01;
errorz[0]=error;
errorz[1]=error;
errorz[2]=error;
errorz[3]=error;
errorz[4]=error;
// the x values
x[0]=1.5751;
x[1]=1.5825;
x[2]=1.6069;
x[3]=1.6339;
x[4]=1.6706;
// the y values
y[0]=1.0642;
y[1]=0.97685;
y[2]=1.13168;
y[3]=1.128654;
y[4]=1.44016;
// initialize TMinuit with a maximum of 5 params
TMinuit *gMinuit = new TMinuit(5);
gMinuit->SetFCN(fcn);
Double_t arglist[10];
Int_t ierflg = 0;
arglist[0] = 1;
gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
// Set starting values and step sizes for parameters
static Double_t vstart[4] = {3, 1 , 0.1 , 0.01};
static Double_t step[4] = {0.1 , 0.1 , 0.01 , 0.001};
gMinuit->mnparm(0, "a1", vstart[0], step[0], 0,0,ierflg);
gMinuit->mnparm(1, "a2", vstart[1], step[1], 0,0,ierflg);
gMinuit->mnparm(2, "a3", vstart[2], step[2], 0,0,ierflg);
gMinuit->mnparm(3, "a4", vstart[3], step[3], 0,0,ierflg);
// Now ready for minimization step
arglist[0] = 500;
arglist[1] = 1.;
gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
// Print results
Double_t amin,edm,errdef;
Int_t nvpar,nparx,icstat;
gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
gMinuit->mnprin(3,amin);
}
Access to the Fit parameters and results
Once an histogram has been fitted, you can access the fitted function parameters via the “Get” functions of the TF1 and TFormula classes. Examples:
Root > TF1 *fit = hist->GetFunction(function_name);
Root > Double_t chi2 = fit->GetChisquare();
Root > Double_t p1 = fit->GetParameter(1);
Root > Double_t e1 = fit->GetParError(1);
Note that above p1 refers to the second parameter.
You can enable the display of the fit parameters in the current canvas by selecting the corresponding item in the Options canvas tool bar menu. You can also set the display of the fit parameters by default in a ROOT session via:
Root > gStyle->SetOptFit</a>();
Fitting with one or more parameters between some bounds
To set bounds for one parameter, use TF1::SetParLimits:
Root > func->SetParLimits(0, -1, 1);
where func is the pointer to the function to be fitted. If you only have the function name, you can get the pointer to this function with:
Root > gROOT->GetFunction(func_name);