Getting details of a gaussian fit

Hello, I need details of a gaussian fit:

Least squares, Maximum Likelihood, fit goodness, residual plot. etc…

For the fit goodness, by default, ROOT print the X^2 in the stat box, …but, how to get the other needed information?

@moneta , surely you can help me
Thank you


Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided


Hi,

You can get the majority of information from the FitResult class, that is returned from TH1::Fit when using the option S. See ROOT: ROOT::Fit::FitResult Class Reference

Lorenzo

Thank you @moneta.
I looked the page but I could only find

Chi2();
Ndf();

I can’t find functions for the other parameters that I need (Least squares, Maximum Likelihood, residual plot. etc…)

you see…there aren’t hit searching these parameters

Hi,
Are you doing a binned maximum likelihood fit or a least square fit of an histogram ?
If doing a maximum likelihood fit, MinFcnValue() will return the maximum likelihood value, while Chi2() will return the corresponding chi2 (least-square value)… If doing a least square fit the two function will return the same value.

For the residual plot, look at TRatioPlot, ROOT: TRatioPlot Class Reference

Lorenzo

Hello @moneta…I will explain…

This is the fit

we used this image for the submission of an article and the referee asked us these parameters…he wrotes:

It would be useful to include the details of the fit: Least squares, Maximum Likelihood, fit goodness, residual plot. etc…

Therefore, being the request of a referee, I’ve to get all these details (least squares, maximum likelihood, fit goodness and residual plot)

I don’t know…I will show you the macro…I fitted in this way

TF1 *fitgaus = new TF1("fitgaus", "gaus", 0,4000);
			auto result= hene4subsum->Fit("fitgaus","S");

and this is the full code for the graph

TCanvas *c19 = new TCanvas("c19",canvtitle,1280,1280);
		t->SetLineColor(kBlue);
		TCut cutsubdet = TString::Format("subdet==%d ", subdetnumb).Data();
		TCut cutxhmin = TString::Format("xh>%d ", xhmin).Data();
		TCut cutxhmax = TString::Format("xh<%d ", xhmax).Data();
		TString hene4substringsum = TString::Format("Calo_EnDep[%d] >> htemp(200., 0., 4000.)", c);
		t->Draw(hene4substringsum, cut && cut3sub && cut4sub && cutsubdet && cutxhmin && cutxhmax);
	 	gPad->Modified();
		gPad->Update(); // make sure it's really (re)drawn
		//t->GetHistogram()->SetTitle(cdenenamesum);
		t->GetHistogram()->SetTitle("");
		TH1F *hene4subsum  = (TH1F*)gPad->GetPrimitive("htemp"); 
		hene4subsum ->GetXaxis()->SetTitle(cdeneXnameab); 
		hene4subsum ->GetYaxis()->SetTitle(cdeneYnameab);
		hene4subsum ->GetYaxis()->SetTitleSize(c_YTitleSize);
		hene4subsum ->GetYaxis()->SetTitleFont(c_YTitleFont);
		hene4subsum ->GetYaxis()->SetTitleOffset(c_YTitleOffset);
		hene4subsum ->GetYaxis()->SetLabelFont(c_YLabelFont); 
		hene4subsum ->GetYaxis()->SetLabelSize(c_YLabelSize);
		hene4subsum ->GetXaxis()->SetTitleSize(c_XTitleSize);
		hene4subsum ->GetXaxis()->SetTitleFont(c_XTitleFont);
		hene4subsum ->GetXaxis()->SetTitleOffset(c_XTitleOffset);
		hene4subsum ->GetXaxis()->SetLabelFont(c_XLabelFont); 
		hene4subsum ->GetXaxis()->SetLabelSize(c_XLabelSize);
		hene4subsum ->SetFillColorAlpha(kBlue, heneabfillcolor);
		hene4subsum ->Draw("E");
		TF1 *fitgaus = new TF1("fitgaus", "gaus", 0,4000);
		auto result= hene4subsum->Fit("fitgaus","S");
		double chi2 = result->Chi2(); 
		double ndf = result->Ndf(); 
		hene4subsum ->SetName(heneabname);
		gPad->Modified();
		gPad->Update(); // make sure it's really (re)drawn
		TLegend* legene4subsum  = new TLegend(0.15, 0.85, .30, .75);
		legene4subsum ->SetHeader("Legend");
		legene4subsum ->SetNColumns(1);
		legene4subsum ->AddEntry(hene4subsum , "Data", "l");
		legene4subsum->AddEntry(fitgaus, "Fit", "l"); 
		legene4subsum ->Draw(); 
		gPad->Update();
		gPad->Update();
		c19->Print(myploteneout4subsum);
		std::cout << "Chi2  : " << chi2 << std::endl;
		std::cout << "Ndf  " << ndf << std::endl;
		results << "Chi2  : " << chi2 << endl;
		results << "Nfd " << ndf<<endl;
		delete c19;	

here the macro…it makes several plots, but this fit is the last one (lines 1356-1402)
calocal.cpp (68.4 KB)

Hi,

Thank you , I understand it now.
You are doing a default fit of an histogram. Then the fit is a Least squares fit using as error for the data point the square root of the observed events.
The fit goodness (chi2 and ndf) you obtain from the fit result, as you are doing in the macro.
The problem with the least square fits are the bins with very small or zero entries. The fit method of ROOT exclude bins with zero entries, but this creates a small bias.
A more correct approach for fitting histogram with counts is to use the binned maximum likelihood method (option L). When doing this you can compute for goodness of fit a chi2 from 2 * maximum likelihood value at the minimum. Since the empty bins are used in the fit, the number of degree of freedom will include now also the empty bins. However, be aware that if you have a large number of empty bins the obtained chi2 value will deviate from following a chi2 distribution.
I attach an example on fitting with the two methods and on how to make a residual plot

Cheers

Lorenzo

exampleFitHistogram.C (1.1 KB)

Dear @moneta thank you.

  1. I tried to fit by Maximum likelihood,
	TF1 *fitgaus = new TF1("fitgaus", "gaus", 0,4000);
			auto result= hene4subsum->Fit("fitgaus","S L");

but
a. you see, the function doesnt’ fit well data in this way

b.
double chi2 = 2* result->MinFcnValue();
gives chi2=1620.75
instead, you see…in the stat box chi2=1932

image

Adding these lines

	auto rp2 = new TRatioPlot(hene4subsum);
   			rp2->Draw();

to make the residual plot I get this error, but in this case, the main plot loses the “E” option

even if I write
rp2->Draw("E");

TCanvas *c19 = new TCanvas("c19",canvtitle,1280,1280);
 			t->SetLineColor(kBlue);
 			TCut cutsubdet = TString::Format("subdet==%d ", subdetnumb).Data();
 			TCut cutxhmin = TString::Format("xh>%d ", xhmin).Data();
 			TCut cutxhmax = TString::Format("xh<%d ", xhmax).Data();
 			TString hene4substringsum = TString::Format("Calo_EnDep[%d] >> htemp(200., 0., 4000.)", c);
			t->Draw(hene4substringsum, cut && cut3sub && cut4sub && cutsubdet && cutxhmin && cutxhmax);
		 	gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
			//t->GetHistogram()->SetTitle(cdenenamesum);
			t->GetHistogram()->SetTitle("");
			TH1F *hene4subsum  = (TH1F*)gPad->GetPrimitive("htemp"); 
			hene4subsum ->GetXaxis()->SetTitle(cdeneXnameab); 
			hene4subsum ->GetYaxis()->SetTitle(cdeneYnameab);
			hene4subsum ->GetYaxis()->SetTitleSize(c_YTitleSize);
			hene4subsum ->GetYaxis()->SetTitleFont(c_YTitleFont);
			hene4subsum ->GetYaxis()->SetTitleOffset(c_YTitleOffset);
			hene4subsum ->GetYaxis()->SetLabelFont(c_YLabelFont); 
			hene4subsum ->GetYaxis()->SetLabelSize(c_YLabelSize);
			hene4subsum ->GetXaxis()->SetTitleSize(c_XTitleSize);
			hene4subsum ->GetXaxis()->SetTitleFont(c_XTitleFont);
			hene4subsum ->GetXaxis()->SetTitleOffset(c_XTitleOffset);
			hene4subsum ->GetXaxis()->SetLabelFont(c_XLabelFont); 
			hene4subsum ->GetXaxis()->SetLabelSize(c_XLabelSize);
			hene4subsum ->SetFillColorAlpha(kBlue, heneabfillcolor);
			hene4subsum ->SetStats(0);
			hene4subsum ->Draw("E");
			TF1 *fitgaus = new TF1("fitgaus", "gaus", 0,4000);
			auto result= hene4subsum->Fit("fitgaus","S L");
			double chi2 = 2* result->MinFcnValue(); //
			double ndf = result->Ndf(); 
   		//	double chi2 = result->Chi2();  // Least squares fit (i.e. without L option)
			hene4subsum ->SetName(heneabname);
   			gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			TLegend* legene4subsum  = new TLegend(0.15, 0.85, .30, .75);
   			legene4subsum ->SetHeader("Legend");
			legene4subsum ->SetNColumns(1);
			legene4subsum ->AddEntry(hene4subsum , "Data", "l");
			legene4subsum->AddEntry(fitgaus, "Fit", "l");  
			legene4subsum ->Draw(); 
			gPad->Update();
			gPad->Update();
			auto rp2 = new TRatioPlot(hene4subsum);
   			rp2->Draw("E");
   			gPad->Update();
			gPad->Update();
			c19->Print(myploteneout4subsum);
			std::cout << "Chi2  : " << chi2 << std::endl;
  			std::cout << "Ndf  " << ndf << std::endl;
   			results << "Chi2  : " << chi2 << endl;
   			results << "Nfd " << ndf<<endl;
			delete c19;

Lastly, in your macro there is the statbox, but, if I plot with the statbox (i.e. I don’t add

hene4subsum ->SetStats(0);

I get this error

Error in <TRatioPlot::BuildLowerPlot>: h1 does not have a fit function

@moneta In your example macro, if you set “gStyle->SetOptFit();”, in the “statistics box” of the “Maximum Likelihood Fit” you will see “chi2 / ndf” which is completely different from “2 * FCN”. The value that is displayed is the same as returned by “result2->Chi2()”.

Hi,
@Wile_E_Coyote Good point! This should be fixed, since the displayed chi2 is the Least Square value and not the 2 * likelihood ratio value, called sometimes Baker-Cousins chi2 (see their paper ).

I will make a PR fixing this

Lorenzo

Hi @moneta

This replies also to my problem 1b in my previous message.
Instead, is there a way to solve:

The problem 1a (i.e. the maximum likelihood doesn’t fit well) as you can see here

  1. Adding the residual plot it loses the “E” drawing and the legend option as you can see here

1 ) The fit is bad because you are having a large non gaussian tail. You should consider modifying your fit function to include the tails, or you can fit in a restricted range, but it will be a biased fit.
When using the Least squares fit (default case in ROOT) the tails are under-estimated because bins with zero entries are not considered.

  1. For fixing legend you should draw the legend in the upper pad using gPad=ratioPlot->GetUpperPad() before drawing the legend.
    For fixing the histogram drawing option use ratioPlot->SetH1DrawOpt("E")

Lorenzo

HEllo @moneta thank you.
I restricted the range and now the fit looks like good

The fit results are
Chi2 = 220.492
Nfd =57

it doesn’t look like a high value…then I think it’s a good result

Just two last questions please,

  1. if I want the statbox, the macro crashes with this error

Error in <TRatioPlot::BuildLowerPlot>: h1 does not have a fit function

  1. Is there a way to define the y-axis range of the residual plot? (now it’s [-0.8;0.8], but I would enlarge it
    I wrote
rp-SetMinimum(-1);
			rp->SetMaximum(1);

but it doewn’t works

TCanvas *c19 = new TCanvas("c19",canvtitle,1280,1280);
 			t->SetLineColor(kBlue);
 			TCut cutsubdet = TString::Format("subdet==%d ", subdetnumb).Data();
 			TCut cutxhmin = TString::Format("xh>%d ", xhmin).Data();
 			TCut cutxhmax = TString::Format("xh<%d ", xhmax).Data();
 			TString hene4substringsum = TString::Format("Calo_EnDep[%d] >> htemp(60., 2800., 4000.)", c);
			t->Draw(hene4substringsum, cut && cut3sub && cut4sub && cutsubdet && cutxhmin && cutxhmax);
		 	gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
			//t->GetHistogram()->SetTitle(cdenenamesum);
			t->GetHistogram()->SetTitle("");
			TH1F *hene4subsum  = (TH1F*)gPad->GetPrimitive("htemp"); 
			hene4subsum ->GetXaxis()->SetTitle(cdeneXnameab); 
			hene4subsum ->GetYaxis()->SetTitle(cdeneYnameab);
			hene4subsum ->GetYaxis()->SetTitleSize(c_YTitleSize);
			hene4subsum ->GetYaxis()->SetTitleFont(c_YTitleFont);
			hene4subsum ->GetYaxis()->SetTitleOffset(c_YTitleOffset);
			hene4subsum ->GetYaxis()->SetLabelFont(c_YLabelFont); 
			hene4subsum ->GetYaxis()->SetLabelSize(c_YLabelSize);
			hene4subsum ->GetXaxis()->SetTitleSize(c_XTitleSize);
			hene4subsum ->GetXaxis()->SetTitleFont(c_XTitleFont);
			hene4subsum ->GetXaxis()->SetTitleOffset(c_XTitleOffset);
			hene4subsum ->GetXaxis()->SetLabelFont(c_XLabelFont); 
			hene4subsum ->GetXaxis()->SetLabelSize(c_XLabelSize);
			hene4subsum ->SetFillColorAlpha(kBlue, heneabfillcolor);
			//hene4subsum ->SetStats(0);
			hene4subsum ->Draw("E");
			TF1 *fitgaus = new TF1("fitgaus", "gaus", 2800,4000);
			auto result= hene4subsum->Fit("fitgaus","S L");
			//auto result= hene4subsum->Fit("fitgaus","S");
			double chi2 = 2* result->MinFcnValue(); //
			double ndf = result->Ndf(); 
   		//	double chi2 = result->Chi2();  // Least squares fit (i.e. without L option)
			hene4subsum ->SetName(heneabname);
   			gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
		//	gPad=ratioPlot->GetUpperPad();
			gPad->Update();
			auto rp = new TRatioPlot(hene4subsum);
			rp->SetH1DrawOpt("E");
   			rp->Draw();
   			gPad->Update();
   			gPad=rp->GetUpperPad();
   			gPad->Update(); // make sure it's really (re)drawn
   			TLegend* legene4subsum  = new TLegend(0.15, 0.85, .30, .75);
   			legene4subsum ->SetHeader("Legend");
			legene4subsum ->SetNColumns(1);
			legene4subsum ->AddEntry(hene4subsum , "Data", "l");
			legene4subsum->AddEntry(fitgaus, "Fit", "l");  
			legene4subsum ->Draw(); 
			gPad->Update();
			gPad->Update();
			c19->Print(myploteneout4subsum);
			std::cout << "Chi2  : " << chi2 << std::endl;
  			std::cout << "Ndf  " << ndf << std::endl;
   			results << "Chi2  : " << chi2 << endl;
   			results << "Nfd " << ndf<<endl;
			delete c19;

Hi,

You can do rp->GetLowerRefYaxis()->SetRange(...) , there might be examples in the TRatioPlot reference documentation.

Hello @moneta, I tried

auto rp = new TRatioPlot(hene4subsum);
			rp->GetLowerRefYaxis()->SetRange(-1.,1.);
			rp->SetH1DrawOpt("E");
   			rp->Draw();

but I get this error

Error in <TRatioPlot::GetLowerRefGraph>: Lower pad does not have primitives

and I get the same error using

rp->GetLowerRefYaxis()->SetTitle("Ratio");

As shown in this example GetLowerRefYaxis needs the lower plot. The lower plot is produce bye Draw . So call GetLowerRefYaxis after Draw.

Thank you @couet it worked

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.