Logo ROOT  
Reference Guide
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Loading...
Searching...
No Matches
FitHistoInFile.C File Reference

Detailed Description

View in nbviewer Open in SWAN

Example for fitting a signal + background model to a histogram found in a file.

This example can be executed as: root > .x FitHistoInFile.C

/// Based on FittingDemo.C by Rene Brun
#include "TH1.h"
#include "TMath.h"
#include "TF1.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TStyle.h"
// Function parameters are passed as an array to TF1. Here, we
// define the position of each parameter in this array.
// Note: N_PAR will give us the total number of parameters. Make
// sure it is always the last entry!
enum ParIndex_t {
Bkg0=0, Bkg1=1, Bkg2,
SigScale, SigSigma, SigMean,
N_PAR};
// Use this map to (re-)name parameters for the plot
const std::map<ParIndex_t,std::string> parNames{
{Bkg0, "Bkg0"}, {Bkg1, "Bkg1"}, {Bkg2, "Bkg2"},
{SigScale, "Gauss scale"}, {SigSigma, "Gauss #sigma"}, {SigMean, "Gauss #mu"}
};
// Quadratic background function
double background(double *x, double *par) {
return par[Bkg0] + par[Bkg1]*x[0] + par[Bkg2]*x[0]*x[0];
}
// Gauss Peak function
double signal(double *x, double *par) {
return par[SigScale]*TMath::Gaus(x[0], par[SigMean], par[SigSigma], true);
}
// Sum of background and peak function. We pass x and the fit parameters
// down to the signal and background functions.
double fitFunction(double *x, double *par) {
return background(x, par) + signal(x, par);
}
// Fit "fitFunction" to the histogram, and draw results on the canvas `c1`.
void FitRoutine(TCanvas* c1, TH1* histo, float fitxmin, float fitxmax, TString filename){
c1->cd();
// create a TF1 with the range from 0 to 3 and N_PAR parameters (six by default)
TF1 fitFcn("fitFcn",fitFunction,fitxmin,fitxmax,N_PAR);
fitFcn.SetNpx(500);
fitFcn.SetLineWidth(2);
fitFcn.SetLineColor(kBlue);
// Assign the names from the map "parNames". Optional, but makes a nicer plot.
for (auto& idx_name : parNames) {
fitFcn.SetParName(idx_name.first, idx_name.second.c_str());
}
// Fit. First set ok-ish starting values for the parameters
fitFcn.SetParameters(30,0,0,50.,0.1,1.);
histo->GetXaxis()->SetRange(2,40);
histo->Fit("fitFcn","VR+","ep");
// improve the picture:
// Draw signal and background functions separately
TF1 backFcn("backFcn",background,fitxmin,fitxmax,N_PAR);
backFcn.SetLineColor(kRed);
TF1 signalFcn("signalFcn",signal,fitxmin,fitxmax,N_PAR);
signalFcn.SetLineColor(kBlue);
signalFcn.SetNpx(500);
// Retrieve fit parameters, and copy them to the signal and background functions
double par[N_PAR];
fitFcn.GetParameters(par);
backFcn.SetParameters(par);
backFcn.DrawCopy("same");
signalFcn.SetParameters(par);
signalFcn.SetLineColor(kGreen);
signalFcn.DrawCopy("same");
const double binwidth = histo->GetXaxis()->GetBinWidth(1);
const double integral = signalFcn.Integral(0.,3.);
std::cout << "number of signal events = " << integral/binwidth << " " << binwidth<< std::endl;
// draw the legend
TLegend legend(0.15,0.7,0.28,0.85);
legend.SetTextFont(72);
legend.SetTextSize(0.03);
legend.AddEntry(histo,"Data","lpe");
legend.AddEntry(&backFcn,"Bgd","l");
legend.AddEntry(&signalFcn,"Sig","l");
legend.AddEntry(&fitFcn,"Sig+Bgd","l");
legend.DrawClone();
histo->Draw("esame");
c1->SaveAs(filename);
}
// Create a file with example data
void CreateRootFile(){
// The data in array form
const int nBins = 60;
double data[nBins] = { 6, 1,10,12, 6,13,23,22,15,21,
23,26,36,25,27,35,40,44,66,81,
75,57,43,37,36,31,35,36,43,32,
40,37,38,33,36,44,42,37,32,32,
43,44,35,33,33,39,29,41,32,44,
26,39,29,35,32,21,21,15,25,15};
TFile* f = new TFile("exampleRootFile.root","RECREATE");
TH1D *histo = new TH1D("histo", "Gauss Peak on Quadratic Background;x;Events/0.05",60,0,3);
for(int i=0; i < nBins; i++) histo->SetBinContent(i+1,data[i]);
f->Write();
f->Close();
}
// Show how to fit a histogram that's in a file.
void FitHistoInFile() {
gStyle->SetOptFit(1111);
// fit range from fitxmin to fitxmax
float fitxmin=0.2;
float fitxmax=2.7;
TCanvas *c1 = new TCanvas("c1","Fitting Demo of Histogram in File",10,10,700,500);
CreateRootFile();
TFile* f = new TFile("exampleRootFile.root");
TH1D* histo= nullptr;
f->GetObject("histo",histo);
if (!histo){
std::cout << "histo not found" << std::endl;
return;
}
histo->SetMarkerStyle(21);
histo->SetMarkerSize(0.8);
// now call the fit routine
FitRoutine(c1,histo, fitxmin, fitxmax,"FitHistoInFile.pdf");
}
#define f(i)
Definition RSha256.hxx:104
@ kRed
Definition Rtypes.h:66
@ kGreen
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition TAttMarker.h:45
virtual void SetRange(Int_t first=0, Int_t last=0)
Set the viewing range for the axis using bin numbers.
Definition TAxis.cxx:1052
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition TAxis.cxx:540
The Canvas class.
Definition TCanvas.h:23
1-Dim function class
Definition TF1.h:214
A ROOT file is composed of a header, followed by consecutive data records (TKey instances) with a wel...
Definition TFile.h:53
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:620
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
TAxis * GetXaxis()
Definition TH1.h:322
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition TH1.cxx:3901
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3067
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition TH1.cxx:9139
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
Basic string class.
Definition TString.h:139
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition TStyle.cxx:1636
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition TStyle.cxx:1589
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculates a gaussian function with mean and sigma.
Definition TMath.cxx:471
Minuit2Minimizer: Minimize with max-calls 1780 convergence for edm < 0.01 strategy 1
Number of iterations 17
----------> Iteration 0
FVAL = 1529.09554298 Edm = 1995.56841581 Nfcn = 25
Error matrix change = 1
Parameters : p0 = 30 p1 = 0 p2 = 0 p3 = 50 p4 = 0.1 p5 = 1
----------> Iteration 1
FVAL = 417.624308317 Edm = 480.13293748 Nfcn = 41
Error matrix change = 0.640972
Parameters : p0 = 21.0028 p1 = -4.96702 p2 = -1.39455 p3 = 27.1085 p4 = 0.12243 p5 = 0.995333
----------> Iteration 2
FVAL = 88.8885076595 Edm = 16.2936179231 Nfcn = 59
Error matrix change = 0.367602
Parameters : p0 = 22.2062 p1 = 1.01384 p2 = 3.74492 p3 = 15.4684 p4 = 0.133674 p5 = 0.994635
----------> Iteration 3
FVAL = 63.8241467861 Edm = 3.17758579488 Nfcn = 73
Error matrix change = 0.34157
Parameters : p0 = 17.9274 p1 = 1.52758 p2 = 4.93171 p3 = 12.2791 p4 = 0.132646 p5 = 0.994209
----------> Iteration 4
FVAL = 55.1442862092 Edm = 0.22822646638 Nfcn = 89
Error matrix change = 0.423967
Parameters : p0 = 14.5751 p1 = 2.94924 p2 = 5.87871 p3 = 14.8197 p4 = 0.131275 p5 = 0.991903
----------> Iteration 5
FVAL = 53.4289073945 Edm = 1.11050817173 Nfcn = 104
Error matrix change = 0.499462
Parameters : p0 = 13.8128 p1 = 4.31214 p2 = 4.77436 p3 = 15.2637 p4 = 0.129147 p5 = 0.984581
----------> Iteration 6
FVAL = 48.4756174099 Edm = 0.26351167941 Nfcn = 118
Error matrix change = 0.560446
Parameters : p0 = 12.0901 p1 = 10.4687 p2 = 2.42177 p3 = 13.7425 p4 = 0.124228 p5 = 0.962145
----------> Iteration 7
FVAL = 46.6349583124 Edm = 1.56781379254 Nfcn = 133
Error matrix change = 0.587229
Parameters : p0 = 11.2362 p1 = 14.7787 p2 = 0.198283 p3 = 11.597 p4 = 0.120788 p5 = 0.958233
----------> Iteration 8
FVAL = 32.3981722197 Edm = 6.55363294533 Nfcn = 151
Error matrix change = 0.634164
Parameters : p0 = 1.02818 p1 = 46.3715 p2 = -15.2613 p3 = 9.25378 p4 = 0.0928931 p5 = 0.946441
----------> Iteration 9
FVAL = 21.4555626905 Edm = 0.277914347616 Nfcn = 166
Error matrix change = 0.325569
Parameters : p0 = -2.66809 p1 = 58.5453 p2 = -20.8916 p3 = 8.43353 p4 = 0.0835661 p5 = 0.983301
----------> Iteration 10
FVAL = 21.2576224712 Edm = 0.0985943659181 Nfcn = 180
Error matrix change = 0.235324
Parameters : p0 = -2.531 p1 = 57.2304 p2 = -20.1351 p3 = 8.91543 p4 = 0.0842051 p5 = 0.98436
----------> Iteration 11
FVAL = 21.0418662198 Edm = 0.163200465587 Nfcn = 194
Error matrix change = 0.121684
Parameters : p0 = -2.93549 p1 = 57.7085 p2 = -20.1524 p3 = 8.78327 p4 = 0.0819441 p5 = 0.980776
----------> Iteration 12
FVAL = 20.1793241317 Edm = 0.148586112923 Nfcn = 211
Error matrix change = 0.0902235
Parameters : p0 = -2.66094 p1 = 58.3868 p2 = -20.6874 p3 = 8.45469 p4 = 0.0722423 p5 = 0.985085
----------> Iteration 13
FVAL = 20.0356153356 Edm = 0.0027813824137 Nfcn = 224
Error matrix change = 0.0782338
Parameters : p0 = -1.6371 p1 = 55.5779 p2 = -19.3667 p3 = 8.43758 p4 = 0.0727503 p5 = 0.984327
----------> Iteration 14
FVAL = 20.033124412 Edm = 3.26053080905e-05 Nfcn = 237
Error matrix change = 0.0415848
Parameters : p0 = -1.72039 p1 = 55.7224 p2 = -19.4085 p3 = 8.41476 p4 = 0.0723863 p5 = 0.984483
----------> Iteration 15
FVAL = 20.0330977486 Edm = 1.41976040596e-07 Nfcn = 251
Error matrix change = 0.0223374
Parameters : p0 = -1.72294 p1 = 55.7228 p2 = -19.4083 p3 = 8.41522 p4 = 0.0724347 p5 = 0.984485
----------> Iteration 16
FVAL = 20.0330977486 Edm = 1.41976040596e-07 Nfcn = 251
Error matrix change = 0.0223374
Parameters : p0 = -1.72294 p1 = 55.7228 p2 = -19.4083 p3 = 8.41522 p4 = 0.0724347 p5 = 0.984485
Minuit2Minimizer : Valid minimum - status = 0
FVAL = 20.0330977485550434
Edm = 1.41976040595893912e-07
Nfcn = 251
Bkg0 = -1.72294 +/- 3.20827
Bkg1 = 55.7228 +/- 8.6473
Bkg2 = -19.4083 +/- 4.14403
Gauss scale = 8.41522 +/- 1.33872
Gauss #sigma = 0.0724347 +/- 0.0116855
Gauss #mu = 0.984485 +/- 0.0109757
Covariance Matrix:
Bkg0 Bkg1 Bkg2 Gauss scaleGauss #sigma Gauss #mu
Bkg0 10.293 -25.273 11.233 1.4782 0.0095082 0.0011726
Bkg1 -25.273 74.776 -35.237 -5.827 -0.037636 0.00042072
Bkg2 11.233 -35.237 17.173 2.8018 0.018113 -0.00087613
Gauss scale 1.4782 -5.827 2.8018 1.7922 0.008954 -0.00050823
Gauss #sigma 0.0095082 -0.037636 0.018113 0.008954 0.00013655 -1.1754e-05
Gauss #mu 0.0011726 0.00042072 -0.00087613 -0.00050823 -1.1754e-05 0.00012047
Correlation Matrix:
Bkg0 Bkg1 Bkg2 Gauss scaleGauss #sigma Gauss #mu
Bkg0 1 -0.91099 0.84489 0.34417 0.25362 0.033301
Bkg1 -0.91099 1 -0.98333 -0.50335 -0.37246 0.0044329
Bkg2 0.84489 -0.98333 1 0.50504 0.37404 -0.019263
Gauss scale 0.34417 -0.50335 0.50504 1 0.57237 -0.034589
Gauss #sigma 0.25362 -0.37246 0.37404 0.57237 1 -0.091642
Gauss #mu 0.033301 0.0044329 -0.019263 -0.034589 -0.091642 1
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 20.0331
NDf = 30
Edm = 1.41976e-07
NCalls = 251
Bkg0 = -1.72294 +/- 3.20827
Bkg1 = 55.7228 +/- 8.6473
Bkg2 = -19.4083 +/- 4.14403
Gauss scale = 8.41522 +/- 1.33872
Gauss #sigma = 0.0724347 +/- 0.0116855
Gauss #mu = 0.984485 +/- 0.0109757
number of signal events = 168.304 0.05
Author
Author E. von Toerne

Definition in file FitHistoInFile.C.