Fitting a 2D Gaussian to a Histogram

Roofit and RooStats
1

Hello,
I am fairly new to ROOT and to C++ and I am having a bit of trouble making a 2D Gaussian function to fit to a histogram. I started by trying to adapt the code from fit2.C into a Gaussian, but I am not having much luck understanding the different functions and parameters used in the code. Would it be possible for someone to explain the parameters, or suggest a way to become more familiar with them? I am attaching my attempt so far but it does not work well.
I can see that the method fills the histogram from a function, but I am not sure how to set up in a way where I can make the function be based off of the histogram.

#include "TF2.h"
#include "TH2.h"
#include "TMath.h"

// Fitting a 2-D histogram
// This tutorial illustrates :
//  - how to create a 2-d function
//  - fill a 2-d histogram randomly from this function
//  - fit the histogram
//  - display the fitted function on top of the histogram 
//
// This example can be executed via the interpreter or ACLIC
//   root > .x fit2.C
//   root > .x fit2.C++
//Author: Rene Brun
         
Double_t g2(Double_t *x, Double_t *par) {
   Double_t r1 = Double_t((x[0]-par[1])/par[2]);
   Double_t r2 = Double_t((x[1]-par[3])/par[4]);
   return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2));
}   
Double_t fun2(Double_t *x, Double_t *par) {
   Double_t *p1 = &par[0];
   Double_t *p2 = &par[5];
   Double_t *p3 = &par[10];
   Double_t result = g2(x,p1) + g2(x,p2) + g2(x,p3);
   return result;
}

void fit2() {
   const Int_t npar = 15;
   Double_t f2params[npar] = 
      {100,-3,3,-3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; //what on earth are all these params for?
   TF2 *f2 = new TF2("f2",fun2,-10,10,-10,10, npar);
   f2->SetParameters(f2params);

   //Create an histogram and fill it randomly with f2
   
   //TH2F *h2 = new TH2F("h2","from f2",40,-10,10,40,-10,10);
   //h2->SetOption("surf2z");
   Int_t nentries = 100000;
   //h2->FillRandom("f2",nentries);
  
   
   TH2F *h2 = new TH2F("h2","2DGaussian",40,-10,10,40,-10,10);
  h2->SetOption("surf2z");
  h2->GetXaxis()->SetTitle("My X Title");
  h2->GetYaxis()->SetTitle("My Y Title");
  int binx1 = h2->GetXaxis()->FindBin(-.3);
  int binx2 = h2->GetXaxis()->FindBin(.3);
  int biny1 = h2->GetYaxis()->FindBin(-.3);
  int biny2 = h2->GetYaxis()->FindBin(.3);
  TRandom3 foo(12345);
  for(int i=0;i<10000;i++)
    {
      h2->Fill(foo.Gaus(0,1),foo.Gaus(0,1));
    }
   
   //Fit h2 with original function f2
   Float_t ratio = 4*nentries/100000;
   f2params[ 0] *= ratio;
   f2params[ 5] *= ratio;
   f2params[10] *= ratio;
   f2->SetParameters(f2params);
   h2->Fit("f2");
   f2->Draw("same, SURF");
}

Sincerely,
Owen

2

Hi Owen,

I suggest to do the following:
o Take the primer of ROOT which targets exactly users which are not experienced with C++ and ROOT:https://root.cern.ch/root/htmldoc/guides/primer/ROOTPrimer.htm . For your particular case, the relevant session is root.cern.ch/root/htmldoc/guide … estimation
o Take a look to the full users guide, in particular the “Fitting Histograms” chapter: root.cern.ch/root/htmldoc/guide … histograms

If you are still lost, don’t hesitate to ask the forum, maybe posting your code.

Cheers,
Danilo

3

Dear Danilo,

Thank you very much, I now understand this macro much better. However, I am still unclear why there are so many parameters. It doesn’t look like the formula uses more than 5. I am attacking my updated and commented version of the code.

#include "TF2.h"
#include "TH2.h"
#include "TMath.h"

// Fitting a 2-D histogram
// This tutorial illustrates :
//  - how to create a 2-d function
//  - fill a 2-d histogram randomly from this function
//  - fit the histogram
//  - display the fitted function on top of the histogram 
//
// This example can be executed via the interpreter or ACLIC
//   root > .x fit2.C
//   root > .x fit2.C++
//Author: Rene Brun
         
Double_t g2(Double_t *x, Double_t *par) { //This seems to be a Gaussian function. Why does the graph look weird?
   Double_t r1 = Double_t((x[0]-par[1])/par[2]);
   Double_t r2 = Double_t((x[1]-par[3])/par[4]);
   return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2));
}   
Double_t fun2(Double_t *x, Double_t *par) { //I don't think I know what these are doing.
   Double_t *p1 = &par[0];
   Double_t *p2 = &par[5];
   Double_t *p3 = &par[10];
   //Double_t result = g2(x,p1) + g2(x,p2) + g2(x,p3);
   Double_t result = g2(x,p1) + g2(x,p2); //This is the formula, but why? Is it just adding each Gaussian, x and y? Removing the third term gets rid of the weird side bump.
   return result;
}

void fit2() {
   const Int_t npar = 15;
   Double_t f2params[npar] = 
      {100,-3,3,-3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; //what on earth are all these parameters for?
   TF2 *f2 = new TF2("f2",fun2,-10,10,-10,10, npar); //"fun2" seems to put together the Gaussian defined in g2.
   f2->SetParameters(f2params);

   //Create an histogram and fill it randomly with f2
   
   //TH2F *h2 = new TH2F("h2","from f2",40,-10,10,40,-10,10);
   //h2->SetOption("surf2z");
   Int_t nentries = 100000;
   //h2->FillRandom("f2",nentries);
  
   
   TH2F *h2 = new TH2F("h2","2DGaussian",40,-10,10,40,-10,10);
  h2->SetOption("surf2z");
  h2->GetXaxis()->SetTitle("My X Title");
  h2->GetYaxis()->SetTitle("My Y Title");
  int binx1 = h2->GetXaxis()->FindBin(-.3);
  int binx2 = h2->GetXaxis()->FindBin(.3);
  int biny1 = h2->GetYaxis()->FindBin(-.3);
  int biny2 = h2->GetYaxis()->FindBin(.3);
  TRandom3 foo(12345);
  for(int i=0;i<10000;i++)
    {
      h2->Fill(foo.Gaus(0,1),foo.Gaus(0,1));
    }
   
   //Fit h2 with original function f2
   Float_t ratio = 4*nentries/100000;
   f2params[ 0] *= ratio; //I have no idea what these are for. Scaling something to the number of entries?
   f2params[ 5] *= ratio;
   f2params[10] *= ratio;
   f2->SetParameters(f2params);
   h2->Fit("f2");
   f2->Draw("same, SURF");
}

Sincerely,
Owen

4

Hi,

the macro shows how to create a non trivial 2d function (fun2), how to fill a 2d histogram sampling it and how to fit a 2d histogram with it (for simplicity but without loss of generality, the function is still fun2).
Now, a bidimensional Gaussian is a functional form with 5 parameters, fun2 is the sum of 3 such Gaussians the number of parameters is 15.