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.
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
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 ) 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.
gPad=ratioPlot->GetUpperPad() before drawing the legend.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,
Error in <TRatioPlot::BuildLowerPlot>: h1 does not have a fit function
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.