Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
fitFraction.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_fit
3/// \notebook
4/// FractionFitter example À la [HMCMLL](https://cds.cern.ch/record/2296378/files/hbook.pdf),
5/// see R. Barlow and C. Beeston, Comp. Phys. Comm. 77 (1993) 219-228,
6/// \see https://indico.in2p3.fr/event/2635/contributions/25070/
7/// \note An alternative interface is described in https://root.cern.ch/doc/master/rf709__BarlowBeeston_8C_source.html
8///
9/// \macro_image
10/// \macro_output
11/// \macro_code
12///
13/// \author
14
15#include "TH1F.h"
16#include "TF1.h"
17#include "TFractionFitter.h"
18#include "TRandom.h"
19#include "TCanvas.h"
20#include "TStyle.h"
21#include "TLatex.h"
22#include <iostream>
23
24void fitFraction()
25{
26 // parameters and functions to generate the data
27 Int_t Ndata = 1000;
28 Int_t N0 = 1000;
29 Int_t N1 = 1000;
30 Int_t N2 = 1000;
31
32 Int_t nBins = 40;
33
34 Double_t trueP0 = .01;
35 Double_t trueP1 = .3;
37
38 // contribution 0
39 TF1 f0("f0", "[0]*(1-cos(x))/TMath::Pi()", 0., TMath::Pi());
40 f0.SetParameter(0, 1.);
41 f0.SetLineColor(2);
42 Double_t int0 = f0.Integral(0., TMath::Pi());
43
44 // contribution 1
45 TF1 f1("f1", "[0]*(1-cos(x)*cos(x))*2./TMath::Pi()", 0., TMath::Pi());
46 f1.SetParameter(0, 1.);
49
50 // contribution 2
51 TF1 f2("f2", "[0]*(1+cos(x))/TMath::Pi()", 0., TMath::Pi());
52 f2.SetParameter(0, 1.);
53 f2.SetLineColor(4);
54 Double_t int2 = f2.Integral(0., TMath::Pi());
55
56 // generate data
57 TH1F data("data", "Data angle distribution", nBins, 0, TMath::Pi()); // data histogram
58 data.SetXTitle("x");
59 data.SetMarkerStyle(20);
60 data.SetMarkerSize(.7);
61 data.SetMinimum(0);
63 htruemc0.SetLineColor(2);
65 htruemc1.SetLineColor(3);
67 htruemc2.SetLineColor(4);
68 Double_t p, x;
69 for (Int_t i = 0; i < Ndata; i++) {
70 p = gRandom->Uniform();
71 if (p < trueP0) {
72 x = f0.GetRandom();
73 htruemc0.Fill(x);
74 } else if (p < trueP0 + trueP1) {
75 x = f1.GetRandom();
76 htruemc1.Fill(x);
77 } else {
78 x = f2.GetRandom();
79 htruemc2.Fill(x);
80 }
81 data.Fill(x);
82 }
83
84 // generate MC samples
85 TH1F mc0("mc0", "MC sample 0 angle distribution", nBins, 0, TMath::Pi()); // first MC histogram
86 mc0.SetXTitle("x");
87 mc0.SetLineColor(2);
88 mc0.SetMarkerColor(2);
89 mc0.SetMarkerStyle(24);
90 mc0.SetMarkerSize(.7);
91 for (Int_t i = 0; i < N0; i++) {
92 mc0.Fill(f0.GetRandom());
93 }
94
95 TH1F mc1("mc1", "MC sample 1 angle distribution", nBins, 0, TMath::Pi()); // second MC histogram
96 mc1.SetXTitle("x");
97 mc1.SetLineColor(3);
98 mc1.SetMarkerColor(3);
99 mc1.SetMarkerStyle(24);
100 mc1.SetMarkerSize(.7);
101 for (Int_t i = 0; i < N1; i++) {
102 mc1.Fill(f1.GetRandom());
103 }
104
105 TH1F mc2("mc2", "MC sample 2 angle distribution", nBins, 0, TMath::Pi()); // third MC histogram
106 mc2.SetXTitle("x");
107 mc2.SetLineColor(4);
108 mc2.SetMarkerColor(4);
109 mc2.SetMarkerStyle(24);
110 mc2.SetMarkerSize(.7);
111 for (Int_t i = 0; i < N2; i++) {
112 mc2.Fill(f2.GetRandom());
113 }
114
115 // FractionFitter
116 TObjArray mc(3); // MC histograms are put in this array
117 mc.Add(&mc0);
118 mc.Add(&mc1);
119 mc.Add(&mc2);
120 TFractionFitter fit(&data, &mc); // initialise
121 fit.Constrain(0, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1
122 fit.Constrain(1, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1
123 fit.Constrain(2, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1
124 // fit.SetRangeX(1,15); // use only the first 15 bins in the fit
125 Int_t status = fit.Fit(); // perform the fit
126 std::cout << "Status: " << status << std::endl;
127
128 // Display
129 gStyle->SetOptStat(0);
130 TCanvas c("c", "FractionFitter example", 700, 700);
131 c.Divide(2, 2);
132
133 c.cd(1);
134 f0.DrawClone();
135 f0.GetHistogram()->SetTitle("Original MC distributions");
136 f1.DrawClone("same");
137 f2.DrawClone("same");
138
139 c.cd(2);
140 data.SetTitle("Data distribution with true contributions");
141 data.DrawClone("EP");
142 htruemc0.Draw("same");
143 htruemc1.Draw("same");
144 htruemc2.Draw("same");
145
146 c.cd(3);
147 mc0.SetTitle("MC generated samples with fit predictions");
148 mc0.Draw("PE");
149 mc1.Draw("PEsame");
150 mc2.Draw("PEsame");
151 if (status == 0) { // check on fit status
152 auto mcp0 = (TH1F *)fit.GetMCPrediction(0);
153 mcp0->SetLineColor(2);
154 mcp0->Draw("same");
155 auto mcp1 = (TH1F *)fit.GetMCPrediction(1);
156 mcp1->SetLineColor(3);
157 mcp1->Draw("same");
158 auto mcp2 = (TH1F *)fit.GetMCPrediction(2);
159 mcp2->SetLineColor(4);
160 mcp2->Draw("same");
161 }
162
163 c.cd(4);
165 TLatex l;
166 l.SetTextSize(.035);
167 Char_t text[200];
168 if (status == 0) { // check on fit status
169 auto result = (TH1F *)fit.GetPlot();
170 fit.GetResult(0, p0, errP0);
171 printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 0, trueP0, p0, errP0);
172 fit.GetResult(1, p1, errP1);
173 printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 1, trueP1, p1, errP1);
174 fit.GetResult(2, p2, errP2);
175 printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 2, trueP2, p2, errP2);
176 data.SetTitle("Data distribution with fitted contributions");
177 data.DrawClone("Ep");
178 result->Draw("same");
179 f0.SetParameter(0, Ndata * p0 / int0 * data.GetBinWidth(1));
180 f0.SetLineStyle(2);
181 f0.DrawClone("same");
182 f1.SetParameter(0, Ndata * p1 / int1 * data.GetBinWidth(1));
183 f1.SetLineStyle(2);
184 f1.DrawClone("same");
185 f2.SetParameter(0, Ndata * p2 / int2 * data.GetBinWidth(1));
186 f2.SetLineStyle(2);
187 f2.DrawClone("same");
188 sprintf(text, "%d: true %.2f, estimated %.2f +/- %.2f\n", 0, trueP0, p0, errP0);
189 l.DrawTextNDC(.45, .30, text);
190 sprintf(text, "%d: true %.2f, estimated %.2f +/- %.2f\n", 1, trueP1, p1, errP1);
191 l.DrawTextNDC(.45, .25, text);
192 sprintf(text, "%d: true %.2f, estimated %.2f +/- %.2f\n", 2, trueP2, p2, errP2);
193 l.DrawTextNDC(.45, .20, text);
194 }
195
196 auto cnew = c.DrawClone();
197 // Cleanup
198 // delete cnew;
199}
#define c(i)
Definition RSha256.hxx:101
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
char Char_t
Character 1 byte (char)
Definition RtypesCore.h:51
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
winID h TVirtualViewer3D TVirtualGLPainter p
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 result
Option_t Option_t TPoint TPoint const char text
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
R__EXTERN TStyle * gStyle
Definition TStyle.h:442
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:44
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:42
The Canvas class.
Definition TCanvas.h:23
1-Dim function class
Definition TF1.h:182
virtual Double_t GetRandom(TRandom *rng=nullptr, Option_t *opt=nullptr)
Return a random number following this function shape.
Definition TF1.cxx:2218
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition TF1.cxx:2557
virtual void SetParameter(Int_t param, Double_t value)
Definition TF1.h:623
Fits MC fractions to data histogram.
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:879
To draw Mathematical Formula.
Definition TLatex.h:18
An array of TObjects.
Definition TObjArray.h:31
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad with: gROOT->SetSelectedPad(c1).
Definition TObject.cxx:319
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition TRandom.cxx:682
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:1642
Double_t x[n]
Definition legend1.C:17
TF1 * f1
Definition legend1.C:11
constexpr Double_t Pi()
Definition TMath.h:38
TLine l
Definition textangle.C:4