Logo ROOT   6.12/07
Reference Guide
TestBinomial.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_fit
3 /// \notebook -js
4 /// Perform a fit to a set of data with binomial errors
5 /// like those derived from the division of two histograms.
6 /// Three different fits are performed and compared:
7 ///
8 /// - simple least square fit to the divided histogram obtained
9 /// from TH1::Divide with option b
10 /// - least square fit to the TGraphAsymmErrors obtained from
11 /// TGraphAsymmErrors::BayesDivide
12 /// - likelihood fit performed on the dividing histograms using
13 /// binomial statistics with the TBinomialEfficiency class
14 ///
15 /// The first two methods are biased while the last one is statistical correct.
16 /// Running the script passing an integer value n larger than 1, n fits are
17 /// performed and the bias are also shown.
18 /// To run the script :
19 ///
20 /// to show the bias performing 100 fits for 1000 events per "experiment"
21 ///
22 /// ~~~{.cpp}
23 /// root[0]: .x TestBinomial.C+
24 /// ~~~
25 ///
26 /// to show the bias performing 100 fits for 1000 events per "experiment"
27 ///
28 /// ~~~{.cpp}
29 /// .x TestBinomial.C+(100, 1000)
30 /// ~~~
31 ///
32 /// \macro_image
33 /// \macro_output
34 /// \macro_code
35 ///
36 /// \author Rene Brun
37 
39 #include "TVirtualFitter.h"
40 #include "TH1.h"
41 #include "TRandom3.h"
42 #include "TF1.h"
43 #include "TFitResult.h"
44 #include "TStyle.h"
45 #include "TCanvas.h"
46 #include "TLegend.h"
47 #include "TPaveStats.h"
48 #include "Math/IntegratorOptions.h"
49 #include <cassert>
50 #include <iostream>
51 
52 void TestBinomial(int nloop = 100, int nevts = 100, bool plot = false, bool debug = false, int seed = 111)
53 {
55  gStyle->SetLineWidth(2.0);
56  gStyle->SetOptStat(11);
57 
58  TObjArray hbiasNorm;
59  hbiasNorm.Add(new TH1D("h0Norm", "Bias Histogram fit",100,-5,5));
60  hbiasNorm.Add(new TH1D("h1Norm","Bias Binomial fit",100,-5,5));
61  TObjArray hbiasThreshold;
62  hbiasThreshold.Add(new TH1D("h0Threshold", "Bias Histogram fit",100,-5,5));
63  hbiasThreshold.Add(new TH1D("h1Threshold","Bias Binomial fit",100,-5,5));
64  TObjArray hbiasWidth;
65  hbiasWidth.Add(new TH1D("h0Width", "Bias Histogram fit",100,-5,5));
66  hbiasWidth.Add(new TH1D("h1Width","Bias Binomial fit",100,-5,5));
67  TH1D* hChisquared = new TH1D("hChisquared",
68  "#chi^{2} probability (Baker-Cousins)", 200, 0.0, 1.0);
69 
72 
73  // Note: in order to be able to use TH1::FillRandom() to generate
74  // pseudo-experiments, we use a trick: generate "selected"
75  // and "non-selected" samples independently. These are
76  // statistically independent and therefore can be safely
77  // added to yield the "before selection" sample.
78 
79 
80  // Define (arbitrarily?) a distribution of input events.
81  // Here: assume a x^(-2) distribution. Boundaries: [10, 100].
82 
83  Double_t xmin =10, xmax = 100;
84  TH1D* hM2D = new TH1D("hM2D", "x^(-2) denominator distribution",
85  45, xmin, xmax);
86  TH1D* hM2N = new TH1D("hM2N", "x^(-2) numerator distribution",
87  45, xmin, xmax);
88  TH1D* hM2E = new TH1D("hM2E", "x^(-2) efficiency",
89  45, xmin, xmax);
90 
91  TF1* fM2D = new TF1("fM2D", "(1-[0]/(1+exp(([1]-x)/[2])))/(x*x)",
92  xmin, xmax);
93  TF1* fM2N = new TF1("fM2N", "[0]/(1+exp(([1]-x)/[2]))/(x*x)",
94  xmin, xmax);
95  TF1* fM2Fit = new TF1("fM2Fit", "[0]/(1+exp(([1]-x)/[2]))",
96  xmin, xmax);
97  TF1* fM2Fit2 = 0;
98 
99  TRandom3 rb(seed);
100 
101  // First try: use a single set of parameters.
102  // For each try, we need to find the overall normalization
103 
104  Double_t normalization = 0.80;
105  Double_t threshold = 25.0;
106  Double_t width = 5.0;
107 
108  fM2D->SetParameter(0, normalization);
109  fM2D->SetParameter(1, threshold);
110  fM2D->SetParameter(2, width);
111  fM2N->SetParameter(0, normalization);
112  fM2N->SetParameter(1, threshold);
113  fM2N->SetParameter(2, width);
114  Double_t integralN = fM2N->Integral(xmin, xmax);
115  Double_t integralD = fM2D->Integral(xmin, xmax);
116  Double_t fracN = integralN/(integralN+integralD);
117  Int_t nevtsN = rb.Binomial(nevts, fracN);
118  Int_t nevtsD = nevts - nevtsN;
119 
120  std::cout << nevtsN << " " << nevtsD << std::endl;
121 
122  gStyle->SetOptFit(1111);
123 
124  // generate many times to see the bias
125  for (int iloop = 0; iloop < nloop; ++iloop) {
126 
127  // generate pseudo-experiments
128  hM2D->Reset();
129  hM2N->Reset();
130  hM2D->FillRandom(fM2D->GetName(), nevtsD);
131  hM2N->FillRandom(fM2N->GetName(), nevtsN);
132  hM2D->Add(hM2N);
133 
134  // construct the "efficiency" histogram
135  hM2N->Sumw2();
136  hM2E->Divide(hM2N, hM2D, 1, 1, "b");
137 
138  // Fit twice, using the same fit function.
139  // In the first (standard) fit, initialize to (arbitrary) values.
140  // In the second fit, use the results from the first fit (this
141  // makes it easier for the fit -- but the purpose here is not to
142  // show how easy or difficult it is to obtain results, but whether
143  // the CORRECT results are obtained or not!).
144 
145  fM2Fit->SetParameter(0, 0.5);
146  fM2Fit->SetParameter(1, 15.0);
147  fM2Fit->SetParameter(2, 2.0);
148  fM2Fit->SetParError(0, 0.1);
149  fM2Fit->SetParError(1, 1.0);
150  fM2Fit->SetParError(2, 0.2);
151  TH1 * hf = fM2Fit->GetHistogram();
152  // std::cout << "Function values " << std::endl;
153  // for (int i = 1; i <= hf->GetNbinsX(); ++i)
154  // std::cout << hf->GetBinContent(i) << " ";
155  // std::cout << std::endl;
156 
157  TCanvas* cEvt;
158  if (plot) {
159  cEvt = new TCanvas(Form("cEnv%d",iloop),
160  Form("plots for experiment %d", iloop),
161  1000, 600);
162  cEvt->Divide(1,2);
163  cEvt->cd(1);
164  hM2D->DrawCopy("HIST");
165  hM2N->SetLineColor(kRed);
166  hM2N->DrawCopy("HIST SAME");
167  cEvt->cd(2);
168  }
169  for (int fit = 0; fit < 2; ++fit) {
170  Int_t status = 0;
171  switch (fit) {
172  case 0:
173  {
174  // TVirtualPad * pad = gPad;
175  // new TCanvas();
176  // fM2Fit->Draw();
177  // gPad = pad;
178  TString optFit = "RN";
179  if (debug) optFit += TString("SV");
180  TFitResultPtr res = hM2E->Fit(fM2Fit, optFit);
181  if (plot) {
182  hM2E->DrawCopy("E");
183  fM2Fit->DrawCopy("SAME");
184  }
185  if (debug) res->Print();
186  status = res;
187  break;
188  }
189  case 1:
190  {
191  // if (fM2Fit2) delete fM2Fit2;
192  // fM2Fit2 = dynamic_cast<TF1*>(fM2Fit->Clone("fM2Fit2"));
193  fM2Fit2 = fM2Fit; // do not clone/copy the function
194  if (fM2Fit2->GetParameter(0) >= 1.0)
195  fM2Fit2->SetParameter(0, 0.95);
196  fM2Fit2->SetParLimits(0, 0.0, 1.0);
197 
198  // TVirtualPad * pad = gPad;
199  // new TCanvas();
200  // fM2Fit2->Draw();
201  // gPad = pad;
202 
203  TBinomialEfficiencyFitter bef(hM2N, hM2D);
204  TString optFit = "RI";
205  if (debug) optFit += TString("SV");
206  TFitResultPtr res = bef.Fit(fM2Fit2,optFit);
207  status = res;
208  if (status !=0) {
209  std::cerr << "Error performing binomial efficiency fit, result = "
210  << status << std::endl;
211  res->Print();
212  continue;
213  }
214  if (plot) {
215  fM2Fit2->SetLineColor(kRed);
216  fM2Fit2->DrawCopy("SAME");
217  }
218  if (debug) {
219  res->Print();
220  }
221  }
222  }
223 
224  if (status != 0) break;
225 
226  Double_t fnorm = fM2Fit->GetParameter(0);
227  Double_t enorm = fM2Fit->GetParError(0);
228  Double_t fthreshold = fM2Fit->GetParameter(1);
229  Double_t ethreshold = fM2Fit->GetParError(1);
230  Double_t fwidth = fM2Fit->GetParameter(2);
231  Double_t ewidth = fM2Fit->GetParError(2);
232  if (fit == 1) {
233  fnorm = fM2Fit2->GetParameter(0);
234  enorm = fM2Fit2->GetParError(0);
235  fthreshold = fM2Fit2->GetParameter(1);
236  ethreshold = fM2Fit2->GetParError(1);
237  fwidth = fM2Fit2->GetParameter(2);
238  ewidth = fM2Fit2->GetParError(2);
239  hChisquared->Fill(fM2Fit2->GetProb());
240  }
241 
242  TH1D* h = dynamic_cast<TH1D*>(hbiasNorm[fit]);
243  h->Fill((fnorm-normalization)/enorm);
244  h = dynamic_cast<TH1D*>(hbiasThreshold[fit]);
245  h->Fill((fthreshold-threshold)/ethreshold);
246  h = dynamic_cast<TH1D*>(hbiasWidth[fit]);
247  h->Fill((fwidth-width)/ewidth);
248  }
249  }
250 
251 
252  TCanvas* c1 = new TCanvas("c1",
253  "Efficiency fit biases",10,10,1000,800);
254  c1->Divide(2,2);
255 
256  TH1D *h0, *h1;
257  c1->cd(1);
258  h0 = dynamic_cast<TH1D*>(hbiasNorm[0]);
259  h0->Draw("HIST");
260  h1 = dynamic_cast<TH1D*>(hbiasNorm[1]);
261  h1->SetLineColor(kRed);
262  h1->Draw("HIST SAMES");
263  TLegend* l1 = new TLegend(0.1, 0.75, 0.5, 0.9,
264  "plateau parameter", "ndc");
265  l1->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
266  %4.2f", h0->GetMean(), h0->GetRMS()), "l");
267  l1->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
268  %4.2f", h1->GetMean(), h1->GetRMS()), "l");
269  l1->Draw();
270 
271  c1->cd(2);
272  h0 = dynamic_cast<TH1D*>(hbiasThreshold[0]);
273  h0->Draw("HIST");
274  h1 = dynamic_cast<TH1D*>(hbiasThreshold[1]);
275  h1->SetLineColor(kRed);
276  h1->Draw("HIST SAMES");
277  TLegend* l2 = new TLegend(0.1, 0.75, 0.5, 0.9,
278  "threshold parameter", "ndc");
279  l2->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
280  %4.2f", h0->GetMean(), h0->GetRMS()), "l");
281  l2->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
282  %4.2f", h1->GetMean(), h1->GetRMS()), "l");
283  l2->Draw();
284 
285  c1->cd(3);
286  h0 = dynamic_cast<TH1D*>(hbiasWidth[0]);
287  h0->Draw("HIST");
288  h1 = dynamic_cast<TH1D*>(hbiasWidth[1]);
289  h1->SetLineColor(kRed);
290  h1->Draw("HIST SAMES");
291  TLegend* l3 = new TLegend(0.1, 0.75, 0.5, 0.9, "width parameter", "ndc");
292  l3->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
293  %4.2f", h0->GetMean(), h0->GetRMS()), "l");
294  l3->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
295  %4.2f", h1->GetMean(), h1->GetRMS()), "l");
296  l3->Draw();
297 
298  c1->cd(4);
299  hChisquared->Draw("HIST");
300 }
301 
302 int main() {
303  TestBinomial();
304 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3244
An array of TObjects.
Definition: TObjArray.h:37
float xmin
Definition: THbookFile.cxx:93
Random number generator class based on M.
Definition: TRandom3.h:27
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
virtual TF1 * DrawCopy(Option_t *option="") const
Draw a copy of this function with its current attributes.
Definition: TF1.cxx:1252
return c1
Definition: legend1.C:41
Definition: Rtypes.h:59
Binomial fitter for the division of two histograms.
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:452
TH1 * h
Definition: legend2.C:5
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:688
static void SetDefaultIntegrator(const char *name)
virtual Double_t GetMean(Int_t axis=1) const
For axis = 1,2 or 3 returns the mean value of the histogram along X,Y or Z axis.
Definition: TH1.cxx:6892
Basic string class.
Definition: TString.h:125
int Int_t
Definition: RtypesCore.h:41
Double_t GetRMS(Int_t axis=1) const
Definition: TH1.h:310
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition: TH1.cxx:3016
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition: TF1.cxx:1815
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition: TF1.cxx:2402
virtual Double_t GetProb() const
Return the fit probability.
Definition: TF1.cxx:1840
int main(int argc, char **argv)
virtual void Print(Option_t *option="") const
Print result of the fit, by default chi2, parameter values and errors.
Definition: TFitResult.cxx:42
TH1F * h1
Definition: legend1.C:5
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9540
virtual Bool_t Divide(TF1 *f1, Double_t c1=1)
Performs the operation: this = this/(c1*f1) if errors are defined (see TH1::Sumw2), errors are also recalculated.
Definition: TH1.cxx:2724
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
Definition: TF1.cxx:3385
virtual void FillRandom(const char *fname, Int_t ntimes=5000)
Fill histogram following distribution in function fname.
Definition: TH1.cxx:3414
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2969
virtual void SetParError(Int_t ipar, Double_t error)
Set error for parameter number ipar.
Definition: TF1.cxx:3360
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
Definition: TFitResultPtr.h:31
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:1218
char * Form(const char *fmt,...)
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
float xmax
Definition: THbookFile.cxx:93
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:610
The Canvas class.
Definition: TCanvas.h:31
double Double_t
Definition: RtypesCore.h:55
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:359
The TH1 histogram class.
Definition: TH1.h:56
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1 if errors are defined (see TH1::Sumw2), errors are also recalculated.
Definition: TH1.cxx:777
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1153
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:496
1-Dim function class
Definition: TF1.h:211
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8276
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:1266
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:284
void Add(TObject *obj)
Definition: TObjArray.h:73
virtual TH1 * GetHistogram() const
Return a pointer to the histogram used to visualise the function.
Definition: TF1.cxx:1469
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:618
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:3688