1 Fitting Histograms

To fit a histogram you can use the Fit Panel on a visible histogram using the GUI, or you can use the TH1::Fit method. The Fit Panel, which is limited, is best for prototyping. The histogram needs to be drawn in a pad before the Fit Panel is available. The TH1::Fit method is more powerful and used in scripts and programs.

The Fit Panel

To display the Fit Panel right click on a histogram to bring up the context menu, then select the menu option: FitPanel.

The first sets of buttons are the predefined functions of ROOT that can be used to fit the histogram. You have a choice of several polynomials, a gaussian, a landau, and an exponential function. You can also define a function and call it "user". It will be linked to the user button on this panel.

You have the option to specify Quiet or Verbose. This is the amount of feedback printed on the root command line on the result of the fit.

When a fit is executed the image of the function is drawn on the current pad. By default the image of the histogram is replaced with the image of the function. Select Same Picture to see the function drawn and the histogram on the same picture.

Select W: Set all weights to 1, to set all errors to 1.

Select E: Compute best errors to use the Minos technique to compute best errors.

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. You can select + : Add to list of functions to 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. Select N: Do not store/draw function to avoid adding the function to the histogram and to avoid drawing it.

Select 0: Do not draw function to avoid drawing the result of the fit.

Select L: Log Likelihood to use loglikelihood method (default is chisquare method).

The slider at the bottom of the panel allows you to set a range for the fit. Drag the edges of the slider towards the center to narrow the range. Draw the entire range to change the beginning and end.

To returns to the original setting, you need press Defaults.

To apply the fit, press the Fit button.

The 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 parameters:

void Fit(const char *fname , Option_t *option , Option_t *goption, Axis_t xxmin, Axis_t xxmax)

*fname:The name of the fitted function (the model) is passed as the first parameter. This name may be one of the of ROOT's pre-defined function names or a user-defined function.

The following functions are predefined, and can be used with the TH1::Fit method.

*option: The second parameter is the fitting option. Here is the list of fitting options:

*goption: The third parameter is the graphics option (goption), it is the same as in the TH1::Draw (see Draw Options above) .

xxmin, xxmax: The 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 is drawn in the current pad.

Fit with a Predefined Function

To fit a histogram with a predefined function, simply pass the name of the function in the first parameter of TH1::Fit. For example, this line fits histogram object hist with a gaussian.

root[] hist.Fit("gaus");

For pre-defined functions, there is no need to set initial values for the parameters, it is done automatically.

Fit with a User- Defined Function

You can create a TF1 object and use it in the call the TH1::Fit. The parameter in to the Fit method is the NAME of the TF1 object.

There are three ways to create a TF1.

  1. Using C++ like expression using x with a fixed set of operators and functions defined in TFormula.

  2. Same as #1, with parameters

  3. Using a function that you have defined

Creating a TF1 with a Formula

Let's look at the first case. Here we call the TF1 constructor by giving it the formula: sin(x)/x.

root[] TF1 *f1 = new TF1("f1", "sin(x)/x", 0,10)

You can also use a TF1 object in the constructor of another TF1.

root[] TF1 *f2 = new TF1("f2", "f1 * 2", 0,10)

Creating a TF1 with Parameters

The second way to construct a TF1 is to add parameters to the expression. For example, this TF1 has 2 parameters:

root[] TF1 *f1 = new TF1("f1","[0]*x*sin([1]*x)",-3,3);

The parameter index is enclosed in square brackets. To set the initial parameters explicitly you can use the SetParameter method.

root[] f1->SetParameter(0,10);

This sets parameter 0 to 10. You can also use SetParameters to set multiple parameters at once.

root[] f1->SetParameters(10,5);

This sets parameter 0 to 10 and parameter 1 to 5.

We can now draw the TF1:

root[] f1->Draw()


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 constructor needs to have this exact signature:

Double_t fitf(Double_t *x, Double_t *par)

The two parameters are:

The following script $ROOTSYS/tutorials/myfit.C illustrates how to fit a 1D histogram with a user-defined function. First we declare the function.

// define a function with 3 parameters

Double_t fitf(Double_t *x, Double_t *par)

{

Double_t arg = 0;

if (par[2]) arg = (x[0] - par[1])/par[2];

Double_t fitval = par[0]*TMath::Exp(-0.5*arg*arg);

return fitval;

}


Now we use the function:

// this function used fitf to fit a histogram

void fitexample()

{

// open a file and get a histogram

TFile *f = new TFile("hsimple.root");

TF1 *hpx = (TF1*)f->Get("hpx");


// create a TF1 object using the function defined above.

// The last 3 specifies the number of parameters

// for the function.

TF1 *func = new TF1 "fit",fitf,-3,3,3);


// set the parameters to the mean and RMS of the histogram

func->SetParameters(500,hpx->GetMean(),hpx->GetRMS());

// give the parameters meaningful names

func->SetParNames ("Constant","Mean_value","Sigma");


// call TH1::Fit with the name of the TF1 object

hpx->Fit ("fit");
}

Fixing and Setting Bounds for Parameters

Parameters must be initialized before invoking the Fit method. The setting of the parameter initial values is automatic for the predefined functions: poln, exp, gaus. You can disable the automatic computation by specifying the "B" option when calling the Fit method.

When a functions is not predefined, the fit paramters must be initialized to some value as close as possible to the expected values before calling the fit function.

To set bounds for one parameter, use TF1::SetParLimits:

func->SetParLimits(0, -1, 1);

When the lower and upper limits are equal, the parameter is fixed. This statement fixes parameter 4 at 10.

func->SetParameter(4,10)

func->SetParLimits(4,77,77);

However, to fix a parameter to 0, one must call the FixParameter function:

func->SetParameter(4,0)

func->FixParameter(4,0);

Note that you are not forced to fix the limits for all parameters. For example, if you fit a function with 6 parameters, you can:

func->SetParameters(0,3.1,1.e-6,-1.5,0,100);

func->SetParLimits(3,-10,-4);

func->FixParameter(4,0);

With this setup, parameters 0->2 can vary freely, parameter 3 has boundaries [-10,-4] with initial value -8, and parameter 4 is fixed to 0.

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 TH1::Fit to restrict the fit to the range specified in the TF1 constructor. In this example, the fit will be limited to -3 to 3, the range specified in the TF1 constructor.

root[] TF1 *f1 = new TF1("f1","[0]*x*sin([1]*x)",-3,3);

root[] hist->Fit("f1", "R");

You can also specify a range in the call to TH1::Fit:

root[] hist->Fit("f1","","",-2,2)

For more complete examples, see $ROOTSYS/tutorials/myfit.C and $ROOTSYS/tutorials/multifit.C.

Example: Fitting Multiple Sub Ranges

The script for this example is in $ROOTSYS/tutorials/multifit.C. It shows how to use several gaussian functions with different parameters on separate sub ranges of the same histogram.

To use a gaussian, or any other ROOT built in function, on a sub range you need to define a new TF1. Each is 'derived' from the canned function gaus.



// Create 4 TF1 objects, one for each subrange

g1 = new TF1("m1","gaus",85,95);

g2 = new TF1("m2","gaus",98,108);

g3 = new TF1("m3","gaus",110,121);

// The total is the sum of the three, each has three //parameters.

total = new TF1("mstotal","gaus(0)+gaus(3)+gaus(6)",85,125);

Here we fill a histogram with bins defined in the array x (see $ROOTSYS/tutorials/multifit.C).

// Create a histogram and set it's contents

h = new TH1F("g1",

"Example of several fits in subranges",np,85,134);

h->SetMaximum(7);

for (int i=0;i<np;i++) {

h->SetBinContent(i+1,x[i]);

}

// Define the parameter array for the total function

Double_t par[9];

When fitting simple functions, such as a gaussian, the initial values of the parameters are automatically computed by ROOT. In the more complicated case of the sum of 3 gaussian functions, the initial values of parameters must be set. In this particular case, the initial values are taken from the result of the individual fits.

The use of the "+" sign is explained below.

//fit each function and add it to the list of functions

h->Fit(g1,"R");

h->Fit(g2,"R+");

h->Fit(g3,"R+");

// Get the parameters from the fit

g1->GetParameters(&par[0]);

g2->GetParameters(&par[3]);

g3->GetParameters(&par[6]);

// Use the parameters on the sum

total->SetParameters(par);

h->Fit(total,"R+");

Adding Functions to The List

The example $ROOTSYS/tutorials/multifit.C also illustrates how to fit several functions on the same histogram. By default a Fit command deletes the previously fitted function in the histogram object. You can specify the option "+" in the second parameter to add the newly fitted function to the existing list of functions for the histogram.

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

You can combine functions to fit a histogram with their sum. Here is an example, the code is in $ROOTSYS/tutorials/FitDemo.C. We have a function that is the combination of a background and lorenzian peak. Each function contributes 3 parameters.

y(E) = a1 + a2E + a3E2 + AP (G / 2 p)/( (E-m)2 + (G/2)2)

background lorenzianPeak

par[0] = a1 par[0] = AP

par[1] = a2 par[1] = G

par[2] = a3 par[2] = m


The combination function (fitFunction) has six parameters:

fitFunction = background (x, par ) + lorenzianPeak (x, &par[3])

par[0] = a1

par[1] = a2

par[2] = a3

par[3] = Ap

par[4] = G

par[5] = m

This script creates a histogram and fits the combination of the two functions. First we define the two functions and the combination function:

// Quadratic background function

Double_t background(Double_t *x, Double_t *par) {

return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];

}



// Lorenzian Peak function

Double_t lorentzianPeak(Double_t *x, Double_t *par) {

return (0.5*par[0]*par[1]/TMath::Pi()) /

TMath::Max( 1.e-10,

(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]

);

}


// Sum of background and peak function

Double_t fitFunction(Double_t *x, Double_t *par) {

return background(x,par) + lorentzianPeak(x,&par[3]);

}


// ? continued on the next page void FittingDemo() {

// Bevington Exercise by Peter Malzacher,

// modified by Rene Brun


const int nBins = 60;

Stat_t data[nBins] = { 6, 1,10,12, 6,13,23,22,15,21,

23,26,36,25,27,35,40,44,66,81,

75,57,48,45,46,41,35,36,53,32,

40,37,38,31,36,44,42,37,32,32,

43,44,35,33,33,39,29,41,32,44,

26,39,29,35,32,21,21,15,25,15};

TH1F *histo = new TH1F("example_9_1",

"Lorentzian Peak on Quadratic Background",60,0,3);

for(int i=0; i < nBins; i++) {

// we use these methods to explicitly set the content

// and error instead of using the fill method.

histo->SetBinContent(i+1,data[i]);

histo->SetBinError(i+1,TMath::Sqrt(data[i]));

}

// create a TF1 with the range from 0 to 3

// and 6 parameters

TF1 *fitFcn = new TF1("fitFcn",fitFunction,0,3,6);

// first try without starting values for the parameters

// This defaults to 1 for each param.

histo->Fit("fitFcn");

// this results in an ok fit for the polynomial function

// however the non-linear part (lorenzian) does not

// respond well.


// second try: set start values for some parameters

fitFcn->SetParameter(4,0.2); // width

fitFcn->SetParameter(5,1); // peak


histo->Fit("fitFcn","V+");

//? continued on next page


// improve the picture:

TF1 *backFcn = new TF1("backFcn",background,0,3,3);

backFcn->SetLineColor(3);

TF1 *signalFcn = new TF1("signalFcn",lorentzianPeak,0,3,3);

signalFcn->SetLineColor(4);

Double_t par[6];

// writes the fit results into the par array

fitFcn->GetParameters(par);

backFcn->SetParameters(par);

backFcn->Draw("same");


signalFcn->SetParameters(&par[3]);

signalFcn->Draw("same");

}



This is the result:

For another example see: http://root.cern.ch/root/html/examples/backsig.C.html

Associated Function

One or more objects (typically a TF1*) can be added to the list of functions (fFunctions) associated to each histogram. A call to TH1::Fit adds the fitted function to this list. Given a histogram h, one can retrieve the associated function with:

TF1 *myfunc = h->GetFunction("myfunc");


Access to the Fit Parameters and Results

If the histogram (or graph) is made persistent, the list of associated functions is also persistent. Retrieve a pointer to the function with the TH1::GetFunction()method. Then you can retrieve the fit parameters from the function (TF1) with calls such as:

root[] TF1 *fit = hist->GetFunction(function_name);

root[] Double_t chi2 = fit->GetChisquare();

// value of the first parameter

root[] Double_t p1 = fit->GetParameter(0);

// errro of the first parameter

root[] Double_t e1 = fit->GetParError(0);


Associated Errors

By default, for each bin, the sum of weights is computed at fill time. One can also call TH1::Sumw2 to force the storage and computation of the sum of the square of weights per bin. If Sumw2 has been called, the error per bin is computed as the sqrt(sum of squares of weights), otherwise the error is set equal to the sqrt (bin content). To return the error for a given bin number, do:

Double_t error = h->GetBinError(bin);

Fit Statistics

You can change the statistics box to display the fit parameters with the TH1::SetOptFit(mode) method. This mode has four digits.

Mode = pcev (default = 0111)

For example:

gStyle->SetOptFit(1011);

This prints the fit probability, parameter names/values, and errors.