Fitting in ROOT

Welcome to the hands-on session dedicated on fitting in ROOT

Exercise 1: Gaussian fit of an histogram

We will start with an exercise where we fin an histogram with a simple function, to get familiar with the fitting options in ROOT.

  • Start creating an histogram with 50 bins in [-5,5] and fill with 1000 Gaussian distributed number
  • Fit the histogram using a Gaussian function
  • Get the value and error of the width of the gaussian
  • Retrieve the fit result and print the correlation matrix.

To solve the exercise you need first to create first the TF1 object, either using the pre-defined gaus function or by using a formula expression with "[0]*ROOT::Math::normal_pdf(x,[2],[1])". Remember that in the second case you need to set the initial function parameters (e.g. f1->SetParameters(1,0,1) ).

To get access to the TFitResult object after fitting use option "S", as shown in slide 10 of the lecture.

Use gStyle->SetOptFit(1) to display the fit result in the statistics box.

#include "TF1.h"
#include "TH1.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"

void gausFit() { 

   TH1D * h1 = new TH1D("h1","h1",50,-5,5);

   h1->FillRandom("gaus",1000);

   TF1 * f1 = new TF1("f1","gaus"); 
   
   // add also option "Q" (quite) to avoid prrinting two time the result
   TFitResultPtr r = h1->Fit(f1,"S Q");
   
   // print the result
   r->Print(); 

   // get the correlation matrix and print it
   TMatrixDSym corrMatrix = r->GetCorrelationMatrix(); 
   corrMatrix.Print(); 

   gStyle->SetOptFit(1);

  // to get the sigma of the gaussian
  std::cout  << "Gaussian sigma = " << f1->GetParameter("sigma") << "  +/- " << f1->GetParError(f1->GetParNumber("sigma")) << std::endl;

}

Pointing hand If you repeat the fit few times (without exiting ROOT) you will see that the sigma you obtain is almost always less than 1. The result is then slightly bias. Try to perform a likelihood fit (option "L"). Is the result better ?

The reason is that the lest square fit is not correct in case of low statistics. The bins with zero events are not included in the fit and this bias the result. The likelihood fit is the correct method for fitting count histograms.

Pointing hand You can notice that the sigma and the amplitude of the gaussian are quite correlated. Can you have a better parametrisation ?

You can fit using a normalised Gaussian. In this case you get a much smaller correlation. To do this create the TF1 as following:
TF1 * f1 = new TF1("f1","[0]*ROOT::Math::normal_pdf(x,[2],[1])");  
// set the parameters (needed if not using a pre-defined function as "gauss")
 f1->SetParameters(1,0,1); 

Exercise 2: Fit a double-peak histogram

We are going to fit the histogram with a more complicated function. First we create an histogram with a double gaussian peaks and then we fit it.

  • Create an histogram with 100 bins between 0 and 10.
  • Fill with 5000 gaussian number with mean 3 and width 1.5.
  • Fill with 5000 gaussian number with a mean of 8 and width of 1.
  • Create a function composed of the sum of two gaussian and fit to the histogram. Do the fit works ? What do you need to do to make the fit working ?

Before fitting you need to set sensible parameter values. You can do this by fitting first a single gaussian in the range [0,5] and then a second gaussian in the range [5,10].

If you don't set good initial parameter values, (for example if you set the amplitude of the gaussians to 1, the means of the gaussians to 5 and widths to 1), the fit will probably not converge.

#include "TF1.h"
#include "TH1.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"
#include "TStyle.h"

void doublePeakFit() { 

   TH1D * h1 = new TH1D("h1","h1",100,0,10);

   for (int i = 0; i < 5000; ++i) { 
      h1->Fill(gRandom->Gaus(3,1.5)); 
      h1->Fill(gRandom->Gaus(8,1)); 
   }

   

   // fit first a single gaussian in range [0,5]
   TFitResultPtr r1 = h1->Fit("gaus","S","",0,5);  // first fit
   TFitResultPtr r2 = h1->Fit("gaus","S","",5,10);  // first fit

   TF1 * f1 = new TF1("fitFunc","gaus(0)+gaus(3)");
   // get parameters and set in global TF1
  
   // parameters of first gaussian
   f1->SetParameter(0,r1->Parameter(0));
   f1->SetParameter(1,r1->Parameter(1));
   f1->SetParameter(2,r1->Parameter(2));
   // parameters of second gaussian
   f1->SetParameter(3,r2->Parameter(0));
   f1->SetParameter(4,r2->Parameter(1));
   f1->SetParameter(5,r2->Parameter(2));


   h1->Fit(f1); 

}

Pointing hand You can try to use a different minimiser, like the genetic minimiser. You need in this case to set parameter limits, otherwise it sill not work, To run the genetic minimiser do:

ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Genetic");

Exercise 3: Using the Fit Panel GUI

Repeat the previous exercise, but by using the Fit Panel GUI

  • Select the fit panel in the Canvas/Tools menu or by right-clicking on the histogram
  • Make the fit function (sum of two gaussian) using the Operation/Add button
  • Set the initial parameters by playing with the sliders
  • Press the Fit button

Pointing hand Explore the other functionalities of the fit panel like changing the fit method, use a different minimiser, plotting contours in the parameters.

Exercise 4: Using random numbers in ROOT

Make 2 histograms between [0,10] and 100 bins.

  • Fill the first histogram with an exponential function, with a decay length of 5 (gRandom->Exp(5)) and a gaussian with mean 5 and sigma 0.5. The number of exponential events is 90% of the total evens (use for example a total number of events of 100000) .
  • Do the same for the second histogram, but using TF1::GetRandom. For doing this create a TF1 which is the sum of a gaussian p.d.f and an exponential p.d.f., as following:
TF1 * f1 = new TF1("f1","0.9*ROOT::Math::exponential_pdf(x, 0.2)+0.1*ROOT::Math::normal_pdf(x,0.5,5)",0,10);
Plot the two histogram in the same canvas.

#include "TRandom.h"
#include "TH1.h"
#include "TF1.h"

void exampleRandom() { 

   TH1D * h1 = new TH1D("h1","h1",100,0,10);
   TH1D * h2 = new TH1D("h2","h2",100,0,10);

   int ntot = 100000;
   int n1 = 0.9*ntot; 
   int n2 = 0.1*ntot; 
   for (int i = 0; i < n1; ++i) 
      h1->Fill(gRandom->Exp(5)); 

   for (int i = 0; i < n2; ++i) 
      h1->Fill(gRandom->Gaus(5,0.5)); 

    // note that the slope parameter of the exponential pdf is the inverse of gRandom->Exp
   TF1 * f1 = new TF1("f1","0.9*ROOT::Math::exponential_pdf(x, 0.2)+0.1*ROOT::Math::normal_pdf(x,0.5,5)",0,10);
   for (int i = 0; i < ntot; ++i) 
      h2->Fill(f1->GetRandom()); 

   

   h1->Draw();
   h2->SetLineColor(kBlue);
   h2->Draw("SAME");
}

Pointing hand Why the 2 histogram are a little different ?

The reason is that TF1::GetRandom() generates the number in the interval [0,10], while gRandom->Exp() generates number between [0,+infinity]. Also the functions ROOT::Math::exponential_pdf and ROOT::Math::normal_pdf are normalised in the full range.

Pointing hand If you have still time, after having finished the exercises you can look at some of the tutorials in the ROOT distribution directory, $ROOTSYS/tutorials/fit. They are available on the Web at this location. For example look at:

Edit | Attach | Watch | Print version | History: r3 < r2 < r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r3 - 2013-02-26 - LorenzoMoneta
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback