/*
 * EXAMPLE: Fit to a uniform distribution
 *
 * Shows clearly that Neyman and Pearson are biased and with poor coverage.
 * The likelihood fit is not biased, and apparently it has 100% coverage.
 * Questions to students:
 *    1. Why Neyman underestimates and Person overestimates?
 *    2. Why the error of the likelihood fit seems wrong? (result doesn't fluctuate)
 *    3. Change code so that error and coverage of likelihood fit look OK
*/    

void Uniform_Binned_Fit(int N=1000) {

// Fill histogram with N random numbers, uniform in (0,100).

auto h1 = new TH1D("h1","Constant Distribution",100,0,100);
TRandom3 r(0);    //  seed=0  ->  different numbers every time
for (int i=0; i<N; ++i) h1->Fill(r.Uniform(0,100));

// Plot histogram

TCanvas* canvas = new TCanvas("canvas","canvas",10,20,700,800);
canvas->Divide(1,2);
canvas->cd(1);
gStyle->SetOptStat(10); // Only number of entries
h1->Draw("e");

// Define a constant function

auto f1 = new TF1("f1","[0]");

// Choose Minuit2 and setup TFitResultPtr pointers for each fit

//ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2"); 
TFitResultPtr result_neyman, result_pearson, result_likelihood;

// Fit histo with f1 using Neyman, Pearson and Likelihood
// (1) Neyman chi2 (default fit, using observed errors)
f1->SetLineColor(kRed);
result_neyman = h1->Fit(f1,"S"); 
// (2) P: Use Pearson chi2 (using expected errors instead of observed errors)
f1->SetLineColor(kBlue);
result_pearson = h1->Fit(f1,"S P+ ");
// (3) L: Use Loglikelihood method (default is chisquare method)
f1->SetLineColor(kGreen+1);
result_likelihood = h1->Fit(f1,"S L+");

// Summarize the three results in a TGraph

TGraphErrors* g1 = new TGraphErrors(1);
TGraphErrors* g2 = new TGraphErrors(1);
TGraphErrors* g3 = new TGraphErrors(1);

g1->SetMarkerStyle(20);
g1->SetMarkerColor(kRed);
g1->SetLineColor(kRed);
g1->SetPoint(0,1,result_neyman->Value(0));
g1->SetPointError(0,0,result_neyman->Error(0));

g2->SetMarkerStyle(20);
g2->SetMarkerColor(kBlue);
g2->SetLineColor(kBlue);
g2->SetPoint(0,2,result_pearson->Value(0));
g2->SetPointError(0,0,result_pearson->Error(0));

g3->SetMarkerStyle(20);
g3->SetMarkerColor(kGreen+1);
g3->SetLineColor(kGreen+1);
g3->SetPoint(0,3,result_likelihood->Value(0));
g3->SetPointError(0,0,result_likelihood->Error(0));

canvas->cd(2);
double trueValue = double(N)/h1->GetNbinsX();  // true result of fit
gPad->DrawFrame(0.5,trueValue*0.85,3.5,trueValue*1.1,"Summary of fit results;;Fit Bias");
g1->Draw("EP0 "); //https://root.cern.ch/root/htmldoc/guides/users-guide/Histograms.html
g2->Draw("EP0 ");
g3->Draw("EP0 ");

TLine* line = new TLine(0.5,trueValue,3.5,trueValue);
line->SetLineStyle(kDashed);
line->SetLineColor(1);
line->Draw();
gPad->Update();

cout << "Neyman     fit bias = " << result_neyman->Value(0)-trueValue<< endl;
cout << "Pearson    fit bias = " << result_pearson->Value(0)-trueValue << endl;
cout << "Likelihood fit bias = " << result_likelihood->Value(0)-trueValue << endl;

}