Fitting in ROOT

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;

}

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.

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); 

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); 

}

• 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

#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");
}

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.