# Fitting a 2D Gaussian to a Histogram

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

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

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

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.