## From \$ROOTSYS/tutorials/fit/TestBinomial.C

```// Perform a fit to a set of data with binomial errors
// like those derived from the division of two histograms.
// Three different fits are performed and compared:
//
//   -  simple least square fit to the divided histogram obtained
//      from TH1::Divide with option b
//   -  least square fit to the TGraphAsymmErrors obtained from
//      TGraphAsymmErrors::BayesDivide
//   -  likelihood fit performed on the dividing histograms using
//      binomial statistics with the TBinomialEfficiency class
//
// The first two methods are biased while the last one  is statistical correct.
// Running the script passing an integer value n larger than 1, n fits are
// performed and the bias are also shown.
// To run the script :
//
//  to show the bias performing 100 fits for 1000 events per "experiment"
//  root[0]: .x TestBinomial.C+
//
//  to show the bias performing 100 fits for 1000 events per "experiment"
//           .x TestBinomial.C+(100, 1000)
//
//
#include "TBinomialEfficiencyFitter.h"
#include "TVirtualFitter.h"
#include "TH1.h"
#include "TRandom3.h"
#include "TF1.h"
#include "TFitResult.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "Math/IntegratorOptions.h"
#include <cassert>
#include <iostream>

void TestBinomial(int nloop = 100, int nevts = 100, bool plot=false, bool debug = false, int seed = 111)
{
gStyle->SetMarkerStyle(20);
gStyle->SetLineWidth(2.0);
gStyle->SetOptStat(11);

TObjArray hbiasNorm;
TObjArray hbiasThreshold;
TObjArray hbiasWidth;
TH1D* hChisquared = new TH1D("hChisquared",
"#chi^{2} probability (Baker-Cousins)", 200, 0.0, 1.0);

TVirtualFitter::SetDefaultFitter("Minuit2");
ROOT::Math::IntegratorOneDimOptions::SetDefaultIntegrator("Gauss");

// Note: in order to be able to use TH1::FillRandom() to generate
//       pseudo-experiments, we use a trick: generate "selected"
//       and "non-selected" samples independently. These are
//       statistically independent and therefore can be safely
//       added to yield the "before selection" sample.

// Define (arbitrarily?) a distribution of input events.
// Here: assume a x^(-2) distribution. Boundaries: [10, 100].

Double_t xmin =10, xmax = 100;
TH1D* hM2D = new TH1D("hM2D", "x^(-2) denominator distribution",
45, xmin, xmax);
TH1D* hM2N = new TH1D("hM2N", "x^(-2) numerator distribution",
45, xmin, xmax);
TH1D* hM2E = new TH1D("hM2E", "x^(-2) efficiency",
45, xmin, xmax);

TF1*  fM2D = new TF1("fM2D", "(1-[0]/(1+exp(([1]-x)/[2])))/(x*x)",
xmin, xmax);
TF1*  fM2N = new TF1("fM2N", "[0]/(1+exp(([1]-x)/[2]))/(x*x)",
xmin, xmax);
TF1*  fM2Fit = new TF1("fM2Fit", "[0]/(1+exp(([1]-x)/[2]))",
xmin, xmax);
TF1*  fM2Fit2 = 0;

TRandom3 rb(seed);

// First try: use a single set of parameters.
// For each try, we need to find the overall normalization

Double_t norm = 0.80;
Double_t threshold = 25.0;
Double_t width = 5.0;

fM2D->SetParameter(0, norm);
fM2D->SetParameter(1, threshold);
fM2D->SetParameter(2, width);
fM2N->SetParameter(0, norm);
fM2N->SetParameter(1, threshold);
fM2N->SetParameter(2, width);
Double_t integralN = fM2N->Integral(xmin, xmax);
Double_t integralD = fM2D->Integral(xmin, xmax);
Double_t fracN = integralN/(integralN+integralD);
Int_t nevtsN = rb.Binomial(nevts, fracN);
Int_t nevtsD = nevts - nevtsN;

std::cout << nevtsN << "  " << nevtsD << std::endl;

gStyle->SetOptFit(1111);

// generate many times to see the bias
for (int iloop = 0; iloop < nloop; ++iloop) {

// generate pseudo-experiments
hM2D->Reset();
hM2N->Reset();
hM2D->FillRandom(fM2D->GetName(), nevtsD);
hM2N->FillRandom(fM2N->GetName(), nevtsN);

// construct the "efficiency" histogram
hM2N->Sumw2();
hM2E->Divide(hM2N, hM2D, 1, 1, "b");

// Fit twice, using the same fit function.
// In the first (standard) fit, initialize to (arbitrary) values.
// In the second fit, use the results from the first fit (this
// makes it easier for the fit -- but the purpose here is not to
// show how easy or difficult it is to obtain results, but whether
// the CORRECT results are obtained or not!).

fM2Fit->SetParameter(0, 0.5);
fM2Fit->SetParameter(1, 15.0);
fM2Fit->SetParameter(2, 2.0);
fM2Fit->SetParError(0, 0.1);
fM2Fit->SetParError(1, 1.0);
fM2Fit->SetParError(2, 0.2);
TH1 * hf = fM2Fit->GetHistogram();
// std::cout << "Function values " << std::endl;
// for (int i = 1; i <= hf->GetNbinsX(); ++i)
//    std::cout << hf->GetBinContent(i) << "  ";
// std::cout << std::endl;

TCanvas* cEvt;
if (plot) {
cEvt = new TCanvas(Form("cEnv%d",iloop),
Form("plots for experiment %d", iloop),
1000, 600);
cEvt->Divide(1,2);
cEvt->cd(1);
hM2D->DrawCopy("HIST");
hM2N->SetLineColor(kRed);
hM2N->DrawCopy("HIST SAME");
cEvt->cd(2);
}
for (int fit = 0; fit < 2; ++fit) {
Int_t status = 0;
switch (fit) {
case 0:
{
// new TCanvas();
// fM2Fit->Draw();
TString optFit = "RN";
if (debug) optFit += TString("SV");
TFitResultPtr res = hM2E->Fit(fM2Fit, optFit);
if (plot) {
hM2E->DrawCopy("E");
fM2Fit->DrawCopy("SAME");
}
if (debug) res->Print();
status = res;
break;
}
case 1:
{
// if (fM2Fit2) delete fM2Fit2;
// fM2Fit2 = dynamic_cast<TF1*>(fM2Fit->Clone("fM2Fit2"));
fM2Fit2 = fM2Fit; // do not clone/copy the function
if (fM2Fit2->GetParameter(0) >= 1.0)
fM2Fit2->SetParameter(0, 0.95);
fM2Fit2->SetParLimits(0, 0.0, 1.0);

// new TCanvas();
// fM2Fit2->Draw();

TBinomialEfficiencyFitter bef(hM2N, hM2D);
TString optFit = "RI";
if (debug) optFit += TString("SV");
TFitResultPtr res = bef.Fit(fM2Fit2,optFit);
status = res;
if (status !=0) {
std::cerr << "Error performing binomial efficiency fit, result = "
<< status << std::endl;
res->Print();
continue;
}
if (plot) {
fM2Fit2->SetLineColor(kRed);
fM2Fit2->DrawCopy("SAME");
}
if (debug) {
res->Print();
}
}
}

if (status != 0) break;

Double_t fnorm = fM2Fit->GetParameter(0);
Double_t enorm = fM2Fit->GetParError(0);
Double_t fthreshold = fM2Fit->GetParameter(1);
Double_t ethreshold = fM2Fit->GetParError(1);
Double_t fwidth = fM2Fit->GetParameter(2);
Double_t ewidth = fM2Fit->GetParError(2);
if (fit == 1) {
fnorm = fM2Fit2->GetParameter(0);
enorm = fM2Fit2->GetParError(0);
fthreshold = fM2Fit2->GetParameter(1);
ethreshold = fM2Fit2->GetParError(1);
fwidth = fM2Fit2->GetParameter(2);
ewidth = fM2Fit2->GetParError(2);
hChisquared->Fill(fM2Fit2->GetProb());
}

TH1D* h = dynamic_cast<TH1D*>(hbiasNorm[fit]);
h->Fill((fnorm-norm)/enorm);
h = dynamic_cast<TH1D*>(hbiasThreshold[fit]);
h->Fill((fthreshold-threshold)/ethreshold);
h = dynamic_cast<TH1D*>(hbiasWidth[fit]);
h->Fill((fwidth-width)/ewidth);
}
}

TCanvas* c1 = new TCanvas("c1",
"Efficiency fit biases",10,10,1000,800);
c1->Divide(2,2);

TH1D *h0, *h1;
c1->cd(1);
h0 = dynamic_cast<TH1D*>(hbiasNorm[0]);
h0->Draw("HIST");
h1 = dynamic_cast<TH1D*>(hbiasNorm[1]);
h1->SetLineColor(kRed);
h1->Draw("HIST SAMES");
TLegend* l1 = new TLegend(0.1, 0.75, 0.5, 0.9,
"plateau parameter", "ndc");
l1->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
%4.2f", h0->GetMean(), h0->GetRMS()), "l");
l1->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
%4.2f", h1->GetMean(), h1->GetRMS()), "l");
l1->Draw();

c1->cd(2);
h0 = dynamic_cast<TH1D*>(hbiasThreshold[0]);
h0->Draw("HIST");
h1 = dynamic_cast<TH1D*>(hbiasThreshold[1]);
h1->SetLineColor(kRed);
h1->Draw("HIST SAMES");
TLegend* l2 = new TLegend(0.1, 0.75, 0.5, 0.9,
"threshold parameter", "ndc");
l2->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
%4.2f", h0->GetMean(), h0->GetRMS()), "l");
l2->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
%4.2f", h1->GetMean(), h1->GetRMS()), "l");
l2->Draw();

c1->cd(3);
h0 = dynamic_cast<TH1D*>(hbiasWidth[0]);
h0->Draw("HIST");
h1 = dynamic_cast<TH1D*>(hbiasWidth[1]);
h1->SetLineColor(kRed);
h1->Draw("HIST SAMES");
TLegend* l3 = new TLegend(0.1, 0.75, 0.5, 0.9, "width parameter", "ndc");
l3->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
%4.2f", h0->GetMean(), h0->GetRMS()), "l");
l3->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
%4.2f", h1->GetMean(), h1->GetRMS()), "l");
l3->Draw();

c1->cd(4);
hChisquared->Draw("HIST");
}

int main() {
TestBinomial();
}
```
TestBinomial.C:1
TestBinomial.C:2
TestBinomial.C:3
TestBinomial.C:4
TestBinomial.C:5
TestBinomial.C:6
TestBinomial.C:7
TestBinomial.C:8
TestBinomial.C:9
TestBinomial.C:10
TestBinomial.C:11
TestBinomial.C:12
TestBinomial.C:13
TestBinomial.C:14
TestBinomial.C:15
TestBinomial.C:16
TestBinomial.C:17
TestBinomial.C:18
TestBinomial.C:19
TestBinomial.C:20
TestBinomial.C:21
TestBinomial.C:22
TestBinomial.C:23
TestBinomial.C:24
TestBinomial.C:25
TestBinomial.C:26
TestBinomial.C:27
TestBinomial.C:28
TestBinomial.C:29
TestBinomial.C:30
TestBinomial.C:31
TestBinomial.C:32
TestBinomial.C:33
TestBinomial.C:34
TestBinomial.C:35
TestBinomial.C:36
TestBinomial.C:37
TestBinomial.C:38
TestBinomial.C:39
TestBinomial.C:40
TestBinomial.C:41
TestBinomial.C:42
TestBinomial.C:43
TestBinomial.C:44
TestBinomial.C:45
TestBinomial.C:46
TestBinomial.C:47
TestBinomial.C:48
TestBinomial.C:49
TestBinomial.C:50
TestBinomial.C:51
TestBinomial.C:52
TestBinomial.C:53
TestBinomial.C:54
TestBinomial.C:55
TestBinomial.C:56
TestBinomial.C:57
TestBinomial.C:58
TestBinomial.C:59
TestBinomial.C:60
TestBinomial.C:61
TestBinomial.C:62
TestBinomial.C:63
TestBinomial.C:64
TestBinomial.C:65
TestBinomial.C:66
TestBinomial.C:67
TestBinomial.C:68
TestBinomial.C:69
TestBinomial.C:70
TestBinomial.C:71
TestBinomial.C:72
TestBinomial.C:73
TestBinomial.C:74
TestBinomial.C:75
TestBinomial.C:76
TestBinomial.C:77
TestBinomial.C:78
TestBinomial.C:79
TestBinomial.C:80
TestBinomial.C:81
TestBinomial.C:82
TestBinomial.C:83
TestBinomial.C:84
TestBinomial.C:85
TestBinomial.C:86
TestBinomial.C:87
TestBinomial.C:88
TestBinomial.C:89
TestBinomial.C:90
TestBinomial.C:91
TestBinomial.C:92
TestBinomial.C:93
TestBinomial.C:94
TestBinomial.C:95
TestBinomial.C:96
TestBinomial.C:97
TestBinomial.C:98
TestBinomial.C:99
TestBinomial.C:100
TestBinomial.C:101
TestBinomial.C:102
TestBinomial.C:103
TestBinomial.C:104
TestBinomial.C:105
TestBinomial.C:106
TestBinomial.C:107
TestBinomial.C:108
TestBinomial.C:109
TestBinomial.C:110
TestBinomial.C:111
TestBinomial.C:112
TestBinomial.C:113
TestBinomial.C:114
TestBinomial.C:115
TestBinomial.C:116
TestBinomial.C:117
TestBinomial.C:118
TestBinomial.C:119
TestBinomial.C:120
TestBinomial.C:121
TestBinomial.C:122
TestBinomial.C:123
TestBinomial.C:124
TestBinomial.C:125
TestBinomial.C:126
TestBinomial.C:127
TestBinomial.C:128
TestBinomial.C:129
TestBinomial.C:130
TestBinomial.C:131
TestBinomial.C:132
TestBinomial.C:133
TestBinomial.C:134
TestBinomial.C:135
TestBinomial.C:136
TestBinomial.C:137
TestBinomial.C:138
TestBinomial.C:139
TestBinomial.C:140
TestBinomial.C:141
TestBinomial.C:142
TestBinomial.C:143
TestBinomial.C:144
TestBinomial.C:145
TestBinomial.C:146
TestBinomial.C:147
TestBinomial.C:148
TestBinomial.C:149
TestBinomial.C:150
TestBinomial.C:151
TestBinomial.C:152
TestBinomial.C:153
TestBinomial.C:154
TestBinomial.C:155
TestBinomial.C:156
TestBinomial.C:157
TestBinomial.C:158
TestBinomial.C:159
TestBinomial.C:160
TestBinomial.C:161
TestBinomial.C:162
TestBinomial.C:163
TestBinomial.C:164
TestBinomial.C:165
TestBinomial.C:166
TestBinomial.C:167
TestBinomial.C:168
TestBinomial.C:169
TestBinomial.C:170
TestBinomial.C:171
TestBinomial.C:172
TestBinomial.C:173
TestBinomial.C:174
TestBinomial.C:175
TestBinomial.C:176
TestBinomial.C:177
TestBinomial.C:178
TestBinomial.C:179
TestBinomial.C:180
TestBinomial.C:181
TestBinomial.C:182
TestBinomial.C:183
TestBinomial.C:184
TestBinomial.C:185
TestBinomial.C:186
TestBinomial.C:187
TestBinomial.C:188
TestBinomial.C:189
TestBinomial.C:190
TestBinomial.C:191
TestBinomial.C:192
TestBinomial.C:193
TestBinomial.C:194
TestBinomial.C:195
TestBinomial.C:196
TestBinomial.C:197
TestBinomial.C:198
TestBinomial.C:199
TestBinomial.C:200
TestBinomial.C:201
TestBinomial.C:202
TestBinomial.C:203
TestBinomial.C:204
TestBinomial.C:205
TestBinomial.C:206
TestBinomial.C:207
TestBinomial.C:208
TestBinomial.C:209
TestBinomial.C:210
TestBinomial.C:211
TestBinomial.C:212
TestBinomial.C:213
TestBinomial.C:214
TestBinomial.C:215
TestBinomial.C:216
TestBinomial.C:217
TestBinomial.C:218
TestBinomial.C:219
TestBinomial.C:220
TestBinomial.C:221
TestBinomial.C:222
TestBinomial.C:223
TestBinomial.C:224
TestBinomial.C:225
TestBinomial.C:226
TestBinomial.C:227
TestBinomial.C:228
TestBinomial.C:229
TestBinomial.C:230
TestBinomial.C:231
TestBinomial.C:232
TestBinomial.C:233
TestBinomial.C:234
TestBinomial.C:235
TestBinomial.C:236
TestBinomial.C:237
TestBinomial.C:238
TestBinomial.C:239
TestBinomial.C:240
TestBinomial.C:241
TestBinomial.C:242
TestBinomial.C:243
TestBinomial.C:244
TestBinomial.C:245
TestBinomial.C:246
TestBinomial.C:247
TestBinomial.C:248
TestBinomial.C:249
TestBinomial.C:250
TestBinomial.C:251
TestBinomial.C:252
TestBinomial.C:253
TestBinomial.C:254
TestBinomial.C:255
TestBinomial.C:256
TestBinomial.C:257
TestBinomial.C:258
TestBinomial.C:259
TestBinomial.C:260
TestBinomial.C:261
TestBinomial.C:262
TestBinomial.C:263
TestBinomial.C:264
TestBinomial.C:265
TestBinomial.C:266
TestBinomial.C:267
TestBinomial.C:268
TestBinomial.C:269
TestBinomial.C:270
TestBinomial.C:271
TestBinomial.C:272
TestBinomial.C:273
TestBinomial.C:274
TestBinomial.C:275
TestBinomial.C:276
TestBinomial.C:277
TestBinomial.C:278
TestBinomial.C:279
TestBinomial.C:280
TestBinomial.C:281
TestBinomial.C:282
TestBinomial.C:283
TestBinomial.C:284
TestBinomial.C:285
TestBinomial.C:286
TestBinomial.C:287
TestBinomial.C:288
TestBinomial.C:289
TestBinomial.C:290
TestBinomial.C:291