Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
FitHistoInFile.C
Go to the documentation of this file.
1/// \file
2/// \notebook
3/// \brief Example for fitting a signal + background model to a histogram found in a file.
4///
5/// This example can be executed as:
6/// root > .x FitHistoInFile.C
7///
8/// \macro_image
9/// \macro_code
10/// \macro_output
11/// \author Author E. von Toerne
12/// Based on FittingDemo.C by Rene Brun
13
14#include "TH1.h"
15#include "TMath.h"
16#include "TF1.h"
17#include "TLegend.h"
18#include "TCanvas.h"
19#include "TFile.h"
20#include "TStyle.h"
21
22// Function parameters are passed as an array to TF1. Here, we
23// define the position of each parameter in this array.
24// Note: N_PAR will give us the total number of parameters. Make
25// sure it is always the last entry!
26enum ParIndex_t {
27 Bkg0=0, Bkg1=1, Bkg2,
28 SigScale, SigSigma, SigMean,
29 N_PAR};
30// Use this map to (re-)name parameters for the plot
31const std::map<ParIndex_t,std::string> parNames{
32 {Bkg0, "Bkg0"}, {Bkg1, "Bkg1"}, {Bkg2, "Bkg2"},
33 {SigScale, "Gauss scale"}, {SigSigma, "Gauss #sigma"}, {SigMean, "Gauss #mu"}
34};
35
36
37// Quadratic background function
38double background(double *x, double *par) {
39 return par[Bkg0] + par[Bkg1]*x[0] + par[Bkg2]*x[0]*x[0];
40}
41
42// Gauss Peak function
43double signal(double *x, double *par) {
44 return par[SigScale]*TMath::Gaus(x[0], par[SigMean], par[SigSigma], true);
45}
46
47// Sum of background and peak function. We pass x and the fit parameters
48// down to the signal and background functions.
49double fitFunction(double *x, double *par) {
50 return background(x, par) + signal(x, par);
51}
52
53// Fit "fitFunction" to the histogram, and draw results on the canvas `c1`.
54void FitRoutine(TCanvas* c1, TH1* histo, float fitxmin, float fitxmax, TString filename){
55 c1->cd();
56 // create a TF1 with the range from 0 to 3 and N_PAR parameters (six by default)
57 TF1 fitFcn("fitFcn",fitFunction,fitxmin,fitxmax,N_PAR);
58 fitFcn.SetNpx(500);
59 fitFcn.SetLineWidth(2);
60 fitFcn.SetLineColor(kBlue);
61
62 // Assign the names from the map "parNames". Optional, but makes a nicer plot.
63 for (auto& idx_name : parNames) {
64 fitFcn.SetParName(idx_name.first, idx_name.second.c_str());
65 }
66
67 // Fit. First set ok-ish starting values for the parameters
68 fitFcn.SetParameters(30,0,0,50.,0.1,1.);
69 histo->GetXaxis()->SetRange(2,40);
70 histo->Fit("fitFcn","VR+","ep");
71
72 // improve the picture:
73 // Draw signal and background functions separately
74 TF1 backFcn("backFcn",background,fitxmin,fitxmax,N_PAR);
75 backFcn.SetLineColor(kRed);
76 TF1 signalFcn("signalFcn",signal,fitxmin,fitxmax,N_PAR);
77 signalFcn.SetLineColor(kBlue);
78 signalFcn.SetNpx(500);
79
80 // Retrieve fit parameters, and copy them to the signal and background functions
81 double par[N_PAR];
82 fitFcn.GetParameters(par);
83
84 backFcn.SetParameters(par);
85 backFcn.DrawCopy("same");
86
87 signalFcn.SetParameters(par);
88 signalFcn.SetLineColor(kGreen);
89 signalFcn.DrawCopy("same");
90
91 const double binwidth = histo->GetXaxis()->GetBinWidth(1);
92 const double integral = signalFcn.Integral(0.,3.);
93 std::cout << "number of signal events = " << integral/binwidth << " " << binwidth<< std::endl;
94
95 // draw the legend
96 TLegend legend(0.15,0.7,0.28,0.85);
97 legend.SetTextFont(72);
98 legend.SetTextSize(0.03);
99 legend.AddEntry(histo,"Data","lpe");
100 legend.AddEntry(&backFcn,"Bgd","l");
101 legend.AddEntry(&signalFcn,"Sig","l");
102 legend.AddEntry(&fitFcn,"Sig+Bgd","l");
103 legend.DrawClone();
104 histo->Draw("esame");
105 c1->SaveAs(filename);
106}
107
108// Create a file with example data
109void CreateRootFile(){
110 // The data in array form
111 const int nBins = 60;
112 double data[nBins] = { 6, 1,10,12, 6,13,23,22,15,21,
113 23,26,36,25,27,35,40,44,66,81,
114 75,57,43,37,36,31,35,36,43,32,
115 40,37,38,33,36,44,42,37,32,32,
116 43,44,35,33,33,39,29,41,32,44,
117 26,39,29,35,32,21,21,15,25,15};
118
119 TFile* f = new TFile("exampleRootFile.root","RECREATE");
120 TH1D *histo = new TH1D("histo", "Gauss Peak on Quadratic Background;x;Events/0.05",60,0,3);
121 for(int i=0; i < nBins; i++) histo->SetBinContent(i+1,data[i]);
122 f->Write();
123 f->Close();
124}
125
126// Show how to fit a histogram that's in a file.
127void FitHistoInFile() {
128 gStyle->SetOptFit(1111);
129 gStyle->SetOptStat(0);
130
131 // fit range from fitxmin to fitxmax
132 float fitxmin=0.2;
133 float fitxmax=2.7;
134
135 TCanvas *c1 = new TCanvas("c1","Fitting Demo of Histogram in File",10,10,700,500);
136 CreateRootFile();
137 TFile* f = new TFile("exampleRootFile.root");
138 TH1D* histo= nullptr;
139 f->GetObject("histo",histo);
140 if (!histo){
141 std::cout << "histo not found" << std::endl;
142 return;
143 }
144 histo->SetMarkerStyle(21);
145 histo->SetMarkerSize(0.8);
146 // now call the fit routine
147 FitRoutine(c1,histo, fitxmin, fitxmax,"FitHistoInFile.pdf");
148}
#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 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
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
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:233
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:664
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
TAxis * GetXaxis()
Definition TH1.h:324
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:3894
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3062
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:9186
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