ROOT   Reference Guide
Searching...
No Matches
TestBinomial.C File Reference

## Detailed Description

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)
32 68
****************************************
Minimizer is Minuit2
Chi2 = 0.350714
NDf = 3
Edm = 1.27001e-06
NCalls = 81
p0 = 0.694132 +/- 0.210029
p1 = 19.3471 +/- 5.85483
p2 = 5.2245 +/- 5.11013
****************************************
Minimizer is Minuit2
Chi2 = 1.64738
NDf = 5
Edm = 1.01171e-06
NCalls = 137
p0 = 0.677329 +/- 0.135567 (limited)
p1 = 15.5028 +/- 4.11118
p2 = 4.10441 +/- 2.77044
****************************************
Minimizer is Minuit2
Chi2 = 0.928456
NDf = 6
Edm = 1.46917e-08
NCalls = 75
p0 = 0.578806 +/- 0.128412 (limited)
p1 = 19.8914 +/- 2.83678
p2 = 3.39323 +/- 1.96156
****************************************
Minimizer is Minuit2
Chi2 = 2.56551
NDf = 7
Edm = 1.19556e-08
NCalls = 109
p0 = 0.57261 +/- 0.16552 (limited)
p1 = 21.4104 +/- 5.10611
p2 = 5.58786 +/- 3.30889
****************************************
Minimizer is Minuit2
Chi2 = 7.37095
NDf = 10
Edm = 1.54804e-06
NCalls = 87
p0 = 0.667689 +/- 0.152841 (limited)
p1 = 24.0146 +/- 2.40617
p2 = 3.63777 +/- 1.74672
****************************************
Minimizer is Minuit2
Chi2 = 0.747776
NDf = 6
Edm = 2.75121e-06
NCalls = 76
p0 = 0.651877 +/- 0.155245 (limited)
p1 = 18.9865 +/- 2.68321
p2 = 2.68989 +/- 1.28756
****************************************
Minimizer is Minuit2
Chi2 = 0.433039
NDf = 4
Edm = 3.07902e-07
NCalls = 104
p0 = 0.500967 +/- 0.254652 (limited)
p1 = 19.4199 +/- 9.30614
p2 = 6.46947 +/- 8.16624
****************************************
Minimizer is Minuit2
Chi2 = 4.49197
NDf = 9
Edm = 4.04441e-06
NCalls = 113
p0 = 0.623236 +/- 0.107197 (limited)
p1 = 20.8129 +/- 1.0507
p2 = 1.22179 +/- 1.09027
****************************************
Minimizer is Minuit2
Chi2 = 2.54966
NDf = 6
Edm = 1.48922e-07
NCalls = 93
p0 = 0.382687 +/- 0.226798 (limited)
p1 = 18.0539 +/- 6.38962
p2 = 3.0464 +/- 5.58751
****************************************
Minimizer is Minuit2
Chi2 = 1.46531
NDf = 7
Edm = 2.36369e-06
NCalls = 90
p0 = 0.575891 +/- 0.187613 (limited)
p1 = 17.6126 +/- 7.68638
p2 = 2.99185 +/- 6.18408
****************************************
Minimizer is Minuit2
Chi2 = 3.31728
NDf = 6
Edm = 8.18927e-07
NCalls = 64
p0 = 0.59853 +/- 0.124912 (limited)
p1 = 17.9048 +/- 3.49062
p2 = 3.52241 +/- 2.09942
****************************************
Minimizer is Minuit2
Chi2 = 0.413866
NDf = 6
Edm = 5.34616e-09
NCalls = 82
p0 = 0.654352 +/- 0.198145 (limited)
p1 = 20.6178 +/- 3.66715
p2 = 4.03019 +/- 2.20627
****************************************
Minimizer is Minuit2
Chi2 = 1.21002
NDf = 4
Edm = 2.90565e-08
NCalls = 82
p0 = 0.651685 +/- 0.159188 (limited)
p1 = 19.3369 +/- 2.66645
p2 = 2.81842 +/- 2.20413
****************************************
Minimizer is Minuit2
Chi2 = 0.752979
NDf = 7
Edm = 5.34806e-06
NCalls = 101
p0 = 0.643852 +/- 0.170146 (limited)
p1 = 21.01 +/- 4.32394
p2 = 5.36976 +/- 3.12567
****************************************
Minimizer is Minuit2
Chi2 = 2.5161
NDf = 8
Edm = 6.68524e-10
NCalls = 61
p0 = 0.628793 +/- 0.103749 (limited)
p1 = 17.4724 +/- 2.12101
p2 = 2.59351 +/- 1.42575
****************************************
Minimizer is Minuit2
Chi2 = 1.27043
NDf = 8
Edm = 6.23845e-06
NCalls = 72
p0 = 0.561127 +/- 0.138134 (limited)
p1 = 18.1601 +/- 3.22498
p2 = 3.50338 +/- 1.84312
****************************************
Minimizer is Minuit2
Chi2 = 2.84972
NDf = 8
Edm = 1.39524e-06
NCalls = 90
p0 = 0.693419 +/- 0.285373 (limited)
p1 = 24.7834 +/- 7.96871
p2 = 5.95505 +/- 3.88805
****************************************
Minimizer is Minuit2
Chi2 = 3.67016
NDf = 5
Edm = 1.80141e-06
NCalls = 82
p0 = 0.502645 +/- 0.176121 (limited)
p1 = 20.1539 +/- 2.79976
p2 = 2.70329 +/- 1.65981
****************************************
Invalid FitResult (status = 3 )
****************************************
Minimizer is Minuit2
Chi2 = 0.936088
NDf = 1
Edm = 0.106873
NCalls = 83
p0 = 0.41816 +/- 0.669502 (limited)
p1 = 18.3049 +/- 3.25491
p2 = 3.06774 +/- 40.715
****************************************
Minimizer is Minuit2
Chi2 = 3.77203
NDf = 8
Edm = 1.56279e-05
NCalls = 84
p0 = 0.642277 +/- 0.195607 (limited)
p1 = 22.9222 +/- 4.97062
p2 = 5.05681 +/- 2.69739
****************************************
Minimizer is Minuit2
Chi2 = 0.47406
NDf = 6
Edm = 1.7448e-06
NCalls = 88
p0 = 0.635811 +/- 0.173792 (limited)
p1 = 21.0478 +/- 3.90109
p2 = 3.91609 +/- 2.85302
****************************************
Minimizer is Minuit2
Chi2 = 2.39556
NDf = 8
Edm = 4.51695e-08
NCalls = 81
p0 = 0.564008 +/- 0.188334 (limited)
p1 = 19.2305 +/- 5.64705
p2 = 4.00431 +/- 3.2294
****************************************
Minimizer is Minuit2
Chi2 = 2.71864
NDf = 3
Edm = 2.02651e-07
NCalls = 84
p0 = 0.430955 +/- 0.103386 (limited)
p1 = 15.3556 +/- 1.58614
p2 = 1.18646 +/- 1.05279
****************************************
Minimizer is Minuit2
Chi2 = 3.91497
NDf = 3
Edm = 9.71484e-06
NCalls = 172
p0 = 0.999991 +/- 0.973989 (limited)
p1 = 26.3679 +/- 4.16357
p2 = 10.2085 +/- 7.97289
****************************************
Minimizer is Minuit2
Chi2 = 3.69474
NDf = 4
Edm = 1.0286e-07
NCalls = 124
p0 = 0.566565 +/- 0.300627 (limited)
p1 = 22.9958 +/- 7.83199
p2 = 6.59909 +/- 5.25595
****************************************
Minimizer is Minuit2
Chi2 = 1.01429
NDf = 10
Edm = 2.60252e-06
NCalls = 94
p0 = 0.610818 +/- 0.175926 (limited)
p1 = 23.3058 +/- 5.77903
p2 = 5.09666 +/- 2.59915
****************************************
Minimizer is Minuit2
Chi2 = 3.63097
NDf = 7
Edm = 2.77767e-07
NCalls = 123
p0 = 0.643028 +/- 0.225348 (limited)
p1 = 21.4306 +/- 4.46994
p2 = 3.95327 +/- 2.04398
****************************************
Minimizer is Minuit2
Chi2 = 1.7669
NDf = 7
Edm = 1.75691e-07
NCalls = 98
p0 = 0.73821 +/- 0.220641 (limited)
p1 = 23.671 +/- 5.18271
p2 = 5.53649 +/- 2.72352
****************************************
Minimizer is Minuit2
Chi2 = 2.70826
NDf = 8
Edm = 5.00596e-07
NCalls = 106
p0 = 0.662591 +/- 0.352545 (limited)
p1 = 24.2783 +/- 12.3384
p2 = 7.60846 +/- 6.66009
****************************************
Minimizer is Minuit2
Chi2 = 9.74855
NDf = 4
Edm = 3.46464e-07
NCalls = 124
p0 = 0.606488 +/- 0.231103 (limited)
p1 = 22.9147 +/- 5.11615
p2 = 4.88595 +/- 2.18611
****************************************
Minimizer is Minuit2
Chi2 = 3.10119
NDf = 5
Edm = 3.71723e-07
NCalls = 234
p0 = 0.999998 +/- 0.501156 (limited)
p1 = 28.5973 +/- 4.11085
p2 = 8.15931 +/- 3.01069
****************************************
Minimizer is Minuit2
Chi2 = 2.78024
NDf = 7
Edm = 1.5102e-06
NCalls = 82
p0 = 0.65308 +/- 0.206386 (limited)
p1 = 23.0571 +/- 4.95138
p2 = 4.86356 +/- 2.47873
****************************************
Minimizer is Minuit2
Chi2 = 3.05436
NDf = 5
Edm = 1.66356e-05
NCalls = 297
p0 = 0.999547 +/- 0.735933 (limited)
p1 = 47.3138 +/- 36.6961
p2 = 43.5419 +/- 61.6394
****************************************
Minimizer is Minuit2
Chi2 = 0.139804
NDf = 2
Edm = 1.06491e-06
NCalls = 91
p0 = 0.660157 +/- 0.258707 (limited)
p1 = 19.2592 +/- 7.6936
p2 = 6.44669 +/- 4.78428
****************************************
Minimizer is Minuit2
Chi2 = 2.55321
NDf = 7
Edm = 5.89829e-06
NCalls = 89
p0 = 0.643594 +/- 0.194962 (limited)
p1 = 21.9893 +/- 5.24123
p2 = 5.90243 +/- 2.90213
****************************************
Minimizer is Minuit2
Chi2 = 2.36128
NDf = 6
Edm = 6.99629e-07
NCalls = 88
p0 = 0.560279 +/- 0.13121 (limited)
p1 = 19.8526 +/- 3.03421
p2 = 2.47661 +/- 2.23047
****************************************
Minimizer is Minuit2
Chi2 = 1.33423
NDf = 7
Edm = 2.89896e-06
NCalls = 102
p0 = 0.564562 +/- 0.203287 (limited)
p1 = 25.8413 +/- 9.79483
p2 = 8.35761 +/- 6.07925
****************************************
Minimizer is Minuit2
Chi2 = 1.74047
NDf = 7
Edm = 7.04231e-07
NCalls = 135
p0 = 0.594764 +/- 0.187954 (limited)
p1 = 25.3877 +/- 8.57244
p2 = 6.52667 +/- 4.25157
****************************************
Minimizer is Minuit2
Chi2 = 0.310648
NDf = 7
Edm = 2.45348e-06
NCalls = 169
p0 = 0.539725 +/- 0.760718 (limited)
p1 = 12.9065 +/- 72.7676
p2 = 22.8828 +/- 100.306
****************************************
Minimizer is Minuit2
Chi2 = 2.52738
NDf = 10
Edm = 1.515e-08
NCalls = 119
p0 = 0.631088 +/- 0.21992 (limited)
p1 = 26.5773 +/- 9.03434
p2 = 8.9826 +/- 4.83446
****************************************
Minimizer is Minuit2
Chi2 = 2.51398
NDf = 6
Edm = 6.94471e-07
NCalls = 84
p0 = 0.484246 +/- 0.226988 (limited)
p1 = 18.2432 +/- 5.68591
p2 = 3.93586 +/- 3.23302
****************************************
Minimizer is Minuit2
Chi2 = 0.946468
NDf = 5
Edm = 5.45857e-06
NCalls = 74
p0 = 0.584615 +/- 0.133896 (limited)
p1 = 19.3728 +/- 3.81107
p2 = 4.20495 +/- 2.53748
****************************************
Minimizer is Minuit2
Chi2 = 5.46724
NDf = 8
Edm = 1.55331e-06
NCalls = 99
p0 = 0.763745 +/- 0.233374 (limited)
p1 = 26.3444 +/- 5.03855
p2 = 6.07972 +/- 2.87905
****************************************
Minimizer is Minuit2
Chi2 = 0.716283
NDf = 7
Edm = 1.89199e-06
NCalls = 85
p0 = 0.422768 +/- 0.179544 (limited)
p1 = 17.3878 +/- 6.43713
p2 = 3.05887 +/- 4.49587
****************************************
Minimizer is Minuit2
Chi2 = 6.40375
NDf = 8
Edm = 5.62321e-07
NCalls = 104
p0 = 0.507615 +/- 0.198193 (limited)
p1 = 22.2924 +/- 5.65071
p2 = 4.8117 +/- 2.67857
****************************************
Minimizer is Minuit2
Chi2 = 3.18171
NDf = 9
Edm = 5.5664e-07
NCalls = 83
p0 = 0.694498 +/- 0.134559 (limited)
p1 = 21.5399 +/- 2.86559
p2 = 3.46906 +/- 1.42287
****************************************
Minimizer is Minuit2
Chi2 = 3.13918
NDf = 6
Edm = 4.11187e-06
NCalls = 101
p0 = 0.725096 +/- 0.145408 (limited)
p1 = 25.2372 +/- 3.82188
p2 = 4.88315 +/- 1.80529
****************************************
Minimizer is Minuit2
Chi2 = 3.15774
NDf = 6
Edm = 2.32136e-06
NCalls = 109
p0 = 0.757722 +/- 0.770168 (limited)
p1 = 26.0106 +/- 13.1967
p2 = 7.01687 +/- 5.25991
****************************************
Minimizer is Minuit2
Chi2 = 1.90008
NDf = 6
Edm = 9.90083e-06
NCalls = 99
p0 = 0.534999 +/- 0.270262 (limited)
p1 = 24.0906 +/- 7.8809
p2 = 5.63935 +/- 3.24361
****************************************
Minimizer is Minuit2
Chi2 = 2.77855
NDf = 8
Edm = 8.13609e-08
NCalls = 85
p0 = 0.702887 +/- 0.199272 (limited)
p1 = 20.9159 +/- 3.68349
p2 = 4.41504 +/- 2.04434
****************************************
Minimizer is Minuit2
Chi2 = 4.37217
NDf = 9
Edm = 2.72083e-07
NCalls = 82
p0 = 0.620959 +/- 0.15843 (limited)
p1 = 21.0013 +/- 3.72087
p2 = 4.56895 +/- 2.37229
****************************************
Minimizer is Minuit2
Chi2 = 1.04047
NDf = 4
Edm = 5.64332e-06
NCalls = 102
p0 = 0.593644 +/- 0.174313 (limited)
p1 = 21.0676 +/- 4.04434
p2 = 3.35047 +/- 2.25854
****************************************
Minimizer is Minuit2
Chi2 = 3.95319
NDf = 10
Edm = 1.32825e-07
NCalls = 96
p0 = 0.574996 +/- 0.152016 (limited)
p1 = 20.2114 +/- 5.1893
p2 = 4.26461 +/- 2.57398
****************************************
Minimizer is Minuit2
Chi2 = 1.2121
NDf = 7
Edm = 3.20149e-06
NCalls = 91
p0 = 0.498388 +/- 0.194573 (limited)
p1 = 16.3792 +/- 10.7582
p2 = 4.04759 +/- 14.7395
****************************************
Minimizer is Minuit2
Chi2 = 1.15843
NDf = 4
Edm = 3.66535e-06
NCalls = 105
p0 = 0.570211 +/- 0.236898 (limited)
p1 = 17.2578 +/- 15.7293
p2 = 9.7494 +/- 16.3952
****************************************
Minimizer is Minuit2
Chi2 = 4.96865
NDf = 8
Edm = 1.88709e-08
NCalls = 80
p0 = 0.666513 +/- 0.140936 (limited)
p1 = 21.2569 +/- 2.26557
p2 = 3.48413 +/- 2.03653
****************************************
Minimizer is Minuit2
Chi2 = 2.88079
NDf = 4
Edm = 8.70733e-06
NCalls = 96
p0 = 0.803697 +/- 0.275939 (limited)
p1 = 22.1124 +/- 5.26388
p2 = 5.42701 +/- 3.06391
****************************************
Minimizer is Minuit2
Chi2 = 3.96392
NDf = 5
Edm = 7.79729e-08
NCalls = 100
p0 = 0.694889 +/- 0.235894 (limited)
p1 = 23.9686 +/- 5.40538
p2 = 5.07476 +/- 3.36335
****************************************
Minimizer is Minuit2
Chi2 = 0.852791
NDf = 5
Edm = 8.88345e-09
NCalls = 84
p0 = 0.310283 +/- 0.164695 (limited)
p1 = 15.0842 +/- 5.86864
p2 = 2.43104 +/- 3.51103
****************************************
Minimizer is Minuit2
Chi2 = 1.32394
NDf = 6
Edm = 6.16879e-07
NCalls = 86
p0 = 0.672106 +/- 0.201669 (limited)
p1 = 23.2133 +/- 5.10201
p2 = 5.04859 +/- 2.50327
****************************************
Minimizer is Minuit2
Chi2 = 2.48447
NDf = 9
Edm = 3.99523e-06
NCalls = 92
p0 = 0.469573 +/- 0.142912 (limited)
p1 = 18.7922 +/- 3.85884
p2 = 2.98044 +/- 2.31127
****************************************
Minimizer is Minuit2
Chi2 = 4.05841
NDf = 7
Edm = 3.35768e-06
NCalls = 93
p0 = 0.667164 +/- 0.148637 (limited)
p1 = 21.9757 +/- 2.55643
p2 = 3.34362 +/- 2.015
****************************************
Minimizer is Minuit2
Chi2 = 3.14114
NDf = 6
Edm = 1.74933e-08
NCalls = 114
p0 = 0.532234 +/- 0.174899 (limited)
p1 = 20.0178 +/- 3.71281
p2 = 3.05935 +/- 2.11129
****************************************
Minimizer is Minuit2
Chi2 = 2.88405
NDf = 5
Edm = 1.25252e-05
NCalls = 77
p0 = 0.698497 +/- 0.183511 (limited)
p1 = 21.6581 +/- 4.18422
p2 = 4.04207 +/- 1.56968
****************************************
Minimizer is Minuit2
Chi2 = 0.624725
NDf = 4
Edm = 6.43892e-06
NCalls = 71
p0 = 0.380579 +/- 0.134096 (limited)
p1 = 13.9323 +/- 2.56334
p2 = 1.96554 +/- 1.64705
****************************************
Minimizer is Minuit2
Chi2 = 1.67529
NDf = 7
Edm = 4.64318e-08
NCalls = 80
p0 = 0.616278 +/- 0.1576 (limited)
p1 = 18.812 +/- 2.91455
p2 = 3.21096 +/- 1.42859
****************************************
Minimizer is Minuit2
Chi2 = 0.374196
NDf = 2
Edm = 8.80185e-06
NCalls = 80
p0 = 0.685904 +/- 0.174442 (limited)
p1 = 19.998 +/- 1.62445
p2 = 1.71179 +/- 1.75684
****************************************
Minimizer is Minuit2
Chi2 = 2.46931
NDf = 7
Edm = 3.32581e-07
NCalls = 93
p0 = 0.640812 +/- 0.135875 (limited)
p1 = 22.3459 +/- 3.96589
p2 = 4.39764 +/- 2.51185
****************************************
Minimizer is Minuit2
Chi2 = 2.49631
NDf = 3
Edm = 6.6234e-06
NCalls = 102
p0 = 0.611704 +/- 0.35871 (limited)
p1 = 20.7057 +/- 8.42105
p2 = 5.08053 +/- 3.9281
****************************************
Minimizer is Minuit2
Chi2 = 0.418946
NDf = 2
Edm = 7.44661e-07
NCalls = 53
p0 = 0.60361 +/- 0.106733 (limited)
p1 = 15.7806 +/- 2
p2 = 0.105266 +/- 2
****************************************
Minimizer is Minuit2
Chi2 = 3.54456
NDf = 5
Edm = 2.33566e-06
NCalls = 85
p0 = 0.527397 +/- 0.217626 (limited)
p1 = 17.9928 +/- 6.38912
p2 = 5.4498 +/- 3.84602
****************************************
Minimizer is Minuit2
Chi2 = 3.23312
NDf = 8
Edm = 9.46056e-07
NCalls = 70
p0 = 0.515942 +/- 0.0980029 (limited)
p1 = 17.2963 +/- 1.44299
p2 = 1.9009 +/- 1.23849
****************************************
Minimizer is Minuit2
Chi2 = 5.28211
NDf = 11
Edm = 1.87776e-06
NCalls = 71
p0 = 0.612568 +/- 0.109907 (limited)
p1 = 19.3198 +/- 2.18618
p2 = 3.37801 +/- 1.31028
****************************************
Minimizer is Minuit2
Chi2 = 2.51155
NDf = 7
Edm = 2.60724e-06
NCalls = 74
p0 = 0.647126 +/- 0.147452 (limited)
p1 = 20.8474 +/- 2.85119
p2 = 3.27715 +/- 1.66359
****************************************
Minimizer is Minuit2
Chi2 = 2.33728
NDf = 6
Edm = 7.37813e-07
NCalls = 109
p0 = 0.624689 +/- 0.141987 (limited)
p1 = 19.3681 +/- 3.57572
p2 = 3.40927 +/- 2.17685
****************************************
Minimizer is Minuit2
Chi2 = 3.75642
NDf = 8
Edm = 5.68096e-06
NCalls = 120
p0 = 0.387014 +/- 0.0702296 (limited)
p1 = 11.5811 +/- 8.22684
p2 = 0.313997 +/- 4.44032
****************************************
Minimizer is Minuit2
Chi2 = 1.68054
NDf = 6
Edm = 5.00607e-06
NCalls = 79
p0 = 0.492763 +/- 0.200749 (limited)
p1 = 18.8758 +/- 5.209
p2 = 3.88829 +/- 2.85902
****************************************
Minimizer is Minuit2
Chi2 = 1.37515
NDf = 4
Edm = 3.40318e-06
NCalls = 106
p0 = 0.759679 +/- 0.234791 (limited)
p1 = 24.2145 +/- 5.44887
p2 = 4.93415 +/- 2.65778
****************************************
Minimizer is Minuit2
Chi2 = 2.35692
NDf = 4
Edm = 1.08773e-07
NCalls = 81
p0 = 0.373315 +/- 0.215845 (limited)
p1 = 12.4572 +/- 7.57085
p2 = 4.00427 +/- 8.97995
****************************************
Minimizer is Minuit2
Chi2 = 1.93137
NDf = 7
Edm = 4.79943e-09
NCalls = 113
p0 = 0.632872 +/- 0.231351 (limited)
p1 = 27.1811 +/- 8.5188
p2 = 6.83092 +/- 3.30519
****************************************
Minimizer is Minuit2
Chi2 = 3.28016
NDf = 6
Edm = 6.04893e-06
NCalls = 99
p0 = 0.696957 +/- 0.193092 (limited)
p1 = 22.9427 +/- 3.44124
p2 = 3.63739 +/- 2.36972
****************************************
Minimizer is Minuit2
Chi2 = 2.56902
NDf = 7
Edm = 5.20615e-07
NCalls = 94
p0 = 0.68832 +/- 0.19603 (limited)
p1 = 24.2988 +/- 4.41739
p2 = 4.93562 +/- 3.1133
****************************************
Minimizer is Minuit2
Chi2 = 3.04459
NDf = 7
Edm = 7.4353e-07
NCalls = 83
p0 = 0.659997 +/- 0.249426 (limited)
p1 = 22.9747 +/- 5.95617
p2 = 5.36022 +/- 2.66473
****************************************
Minimizer is Minuit2
Chi2 = 6.27397
NDf = 8
Edm = 5.29961e-07
NCalls = 81
p0 = 0.662433 +/- 0.15714 (limited)
p1 = 21.5909 +/- 3.29623
p2 = 4.09179 +/- 1.74745
****************************************
Minimizer is Minuit2
Chi2 = 1.94311
NDf = 9
Edm = 1.89133e-07
NCalls = 106
p0 = 0.575157 +/- 0.298322 (limited)
p1 = 25.6868 +/- 11.7079
p2 = 7.54114 +/- 5.74853
****************************************
Minimizer is Minuit2
Chi2 = 2.75509
NDf = 9
Edm = 9.98946e-08
NCalls = 69
p0 = 0.593401 +/- 0.101469 (limited)
p1 = 18.8032 +/- 2.57519
p2 = 2.53183 +/- 1.7654
****************************************
Minimizer is Minuit2
Chi2 = 2.54952
NDf = 6
Edm = 4.57669e-08
NCalls = 96
p0 = 0.660242 +/- 0.259374 (limited)
p1 = 25.0771 +/- 5.98352
p2 = 5.39543 +/- 2.98249
****************************************
Minimizer is Minuit2
Chi2 = 4.62283
NDf = 10
Edm = 1.36401e-08
NCalls = 89
p0 = 0.592713 +/- 0.183874 (limited)
p1 = 20.9349 +/- 5.53213
p2 = 6.39618 +/- 4.21336
****************************************
Minimizer is Minuit2
Chi2 = 5.3127
NDf = 6
Edm = 2.53153e-06
NCalls = 119
p0 = 0.534666 +/- 0.0847463 (limited)
p1 = 12.0856 +/- 13.9766
p2 = 0.685012 +/- 8.82276
****************************************
Minimizer is Minuit2
Chi2 = 0.968364
NDf = 7
Edm = 1.15106e-06
NCalls = 74
p0 = 0.669971 +/- 0.134042 (limited)
p1 = 19.9248 +/- 3.66311
p2 = 4.7186 +/- 3.31403
****************************************
Minimizer is Minuit2
Chi2 = 2.88345
NDf = 10
Edm = 2.24225e-06
NCalls = 102
p0 = 0.629842 +/- 0.174064 (limited)
p1 = 21.9766 +/- 6.0423
p2 = 6.83408 +/- 3.9968
****************************************
Minimizer is Minuit2
Chi2 = 1.80215
NDf = 6
Edm = 2.68316e-07
NCalls = 87
p0 = 0.495591 +/- 0.211145 (limited)
p1 = 19.645 +/- 5.2052
p2 = 3.59211 +/- 2.51889
****************************************
Minimizer is Minuit2
Chi2 = 3.79465
NDf = 6
Edm = 4.35871e-08
NCalls = 69
p0 = 0.611171 +/- 0.128277 (limited)
p1 = 17.467 +/- 2.93532
p2 = 3.02848 +/- 1.64432
****************************************
Minimizer is Minuit2
Chi2 = 2.748
NDf = 7
Edm = 5.78153e-09
NCalls = 140
p0 = 1 +/- 0.733966 (limited)
p1 = 28.7334 +/- 2.38732
p2 = 6.33261 +/- 1.85434
****************************************
Minimizer is Minuit2
Chi2 = 4.29305
NDf = 9
Edm = 3.06585e-06
NCalls = 82
p0 = 0.687721 +/- 0.16847 (limited)
p1 = 22.2779 +/- 3.60769
p2 = 4.13806 +/- 1.75728
****************************************
Minimizer is Minuit2
Chi2 = 1.74986
NDf = 2
Edm = 2.69957e-06
NCalls = 256
p0 = 0.999951 +/- 0.988093 (limited)
p1 = 33.6961 +/- 9.16268
p2 = 13.7955 +/- 8.24703
****************************************
Minimizer is Minuit2
Chi2 = 4.72905
NDf = 5
Edm = 2.89445e-07
NCalls = 67
p0 = 0.417412 +/- 0.107227 (limited)
p1 = 14.5796 +/- 2.95636
p2 = 2.22066 +/- 2.12966
****************************************
Minimizer is Minuit2
Chi2 = 2.45643
NDf = 6
Edm = 2.63752e-06
NCalls = 99
p0 = 0.554001 +/- 0.226655 (limited)
p1 = 23.9121 +/- 6.69531
p2 = 6.04267 +/- 3.44117
****************************************
Minimizer is Minuit2
Chi2 = 2.91637
NDf = 4
Edm = 1.377e-05
NCalls = 84
p0 = 0.769946 +/- 0.185015 (limited)
p1 = 20.5779 +/- 2.52006
p2 = 3.50084 +/- 1.35396
****************************************
Minimizer is Minuit2
Chi2 = 1.21412
NDf = 6
Edm = 1.16053e-08
NCalls = 141
p0 = 0.632519 +/- 0.106789 (limited)
p1 = 17.4834 +/- 1.66297
p2 = 1.44082 +/- 1.10945
#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 <cassert>
#include <iostream>
void TestBinomial(int nloop = 100, int nevts = 100, bool plot = false, bool debug = false, int seed = 111)
{
TObjArray hbiasNorm;
TObjArray hbiasThreshold;
TObjArray hbiasWidth;
TH1D* hChisquared = new TH1D("hChisquared",
"#chi^{2} probability (Baker-Cousins)", 200, 0.0, 1.0);
// 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)",
TF1* fM2N = new TF1("fM2N", "[0]/(1+exp(([1]-x)/[2]))/(x*x)",
TF1* fM2Fit = new TF1("fM2Fit", "[0]/(1+exp(([1]-x)/[2]))",
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 normalization = 0.80;
Double_t threshold = 25.0;
Double_t width = 5.0;
fM2D->SetParameter(0, normalization);
fM2D->SetParameter(1, threshold);
fM2D->SetParameter(2, width);
fM2N->SetParameter(0, normalization);
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->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->SetLineColor(kBlue);
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 S";
if (debug) optFit += TString("V");
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");
bool confint = (status == 0);
if (confint) {
// compute confidence interval on fitted function
auto htemp = fM2Fit2->GetHistogram();
ROOT::Fit::FillData(xdata, fM2Fit2->GetHistogram() );
res->GetConfidenceIntervals(xdata, gr.GetEY(), 0.68, false);
gr.DrawClone("4 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-normalization)/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->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->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->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();
}
#define h(i)
Definition RSha256.hxx:106
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
@ kRed
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
include TDocParser_001 C image html pict1_TDocParser_001 png width
float xmin
float xmax
char * Form(const char *fmt,...)
R__EXTERN TStyle * gStyle
Definition TStyle.h:412
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition BinData.h:52
void GetConfidenceIntervals(unsigned int n, unsigned int stride1, unsigned int stride2, const double *x, double *ci, double cl=0.95, bool norm=false) const
get confidence intervals for an array of n points x.
static void SetDefaultIntegrator(const char *name)
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:37
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition TAttFill.h:39
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40
Binomial fitter for the division of two histograms.
The Canvas class.
Definition TCanvas.h:23
Definition TCanvas.cxx:708
1-Dim function class
Definition TF1.h:213
virtual void SetParError(Int_t ipar, Double_t error)
Set error for parameter number ipar.
Definition TF1.cxx:3486
virtual TH1 * GetHistogram() const
Return a pointer to the histogram used to visualise the function Note that this histogram is managed ...
Definition TF1.cxx:1574
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition TF1.cxx:1920
virtual Double_t GetProb() const
Return the fit probability.
Definition TF1.cxx:1945
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition TF1.cxx:2515
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
Definition TF1.cxx:3511
virtual TF1 * DrawCopy(Option_t *option="") const
Draw a copy of this function with its current attributes.
Definition TF1.cxx:1352
virtual void SetParameter(Int_t param, Double_t value)
Definition TF1.h:634
virtual Double_t GetParameter(Int_t ipar) const
Definition TF1.h:512
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
virtual void Print(Option_t *option="") const
Print result of the fit, by default chi2, parameter values and errors.
A TGraphErrors is a TGraph with error bars.
Double_t * GetEY() const
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:618
virtual void Reset(Option_t *option="")
Reset.
Definition TH1.cxx:10126
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
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:7428
virtual void FillRandom(const char *fname, Int_t ntimes=5000, TRandom *rng=nullptr)
Fill histogram following distribution in function fname.
Definition TH1.cxx:3525
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:3892
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),...
Definition TH1.cxx:822
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
Double_t GetRMS(Int_t axis=1) const
Definition TH1.h:315
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition TH1.cxx:3120
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3073
virtual Bool_t Divide(TF1 *f1, Double_t c1=1)
Performs the operation: this = this/(c1*f1) if errors are defined (see TH1::Sumw2),...
Definition TH1.cxx:2828
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition TH1.cxx:8860
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition TLegend.cxx:330
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition TLegend.cxx:423
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
An array of TObjects.
Definition TObjArray.h:37
Definition TObjArray.h:74
virtual TObject * DrawClone(Option_t *option="") const
Definition TObject.cxx:221
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Random number generator class based on M.
Definition TRandom3.h:27
Basic string class.
Definition TString.h:136
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:1589
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:1541
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
int main()
return c1
Definition legend1.C:41
TGraphErrors * gr
Definition legend1.C:25
TH1F * h1
Definition legend1.C:5