#include <iostream>
#include <sstream>

#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TCanvas.h"
#include "TString.h"
#include "TStyle.h"
#include "TChain.h"
#include "TH2.h"
#include "TH1.h"
#include "TF1.h"
#include "TTree.h"
#include "TKey.h"
#include "Riostream.h"
#include "TMath.h"
#include "TRandom.h"
#include "TRandom3.h"
#include "TVirtualFitter.h"

#include <stdio.h>
#include <math.h>

void Chi2vsLik (int mumax=100, int Nbins=100, int Maxrep=100) {

  TH1D * data = new TH1D ("data", "", Nbins, 0., (double)Nbins);
  TH1D * Pearson = new TH1D ("Pearson", "", mumax, 0.5, 0.5+(double)mumax);
  TH1D * Neyman = new TH1D ("Neyman", "", mumax, 0.5, 0.5+(double)mumax);
  TH1D * Likelihood = new TH1D ("Likelihood", "", mumax, 0.5, 0.5+(double)mumax);

  for (int i=1; i<=mumax; i++) {
    double mu = (double)i;

    double bigsumP=0;
    double bigsumN=0;
    double bigsumL=0;
    for (int Nrep=0; Nrep<Maxrep; Nrep++) {
      
      data->Reset();
      for (int j=0; j<Nbins; j++) {
	double k = gRandom->Poisson(mu);
	data->SetBinContent(j+1,k);
      }
      double sum2=0;
      double suminv=0;
      double sum=0;
      for (int j=0; j<Nbins; j++) {
	double k=data->GetBinContent(j+1);
	sum2+=k*k/Nbins;
	suminv+=1/k;
	sum+=k/Nbins;
      }
      bigsumP+=sqrt(sum2)/mu;
      bigsumN+=Nbins/suminv/mu;
      bigsumL+=sum/mu;
    }
    Pearson->SetBinContent(i,bigsumP/Maxrep);
    Neyman->SetBinContent(i,bigsumN/Maxrep);
    Likelihood->SetBinContent(i,bigsumL/Maxrep);
  }


  TCanvas * C = new TCanvas ("C","",500,500);
  C->cd();
  Pearson->SetMinimum(0);
  Pearson->SetLineWidth(3);
  Neyman->SetLineWidth(3);
  Likelihood->SetLineWidth(3);
  Pearson->SetLineColor(kRed);
  Neyman->SetLineColor(kBlue);
  Likelihood->SetLineColor(kBlack);
  Pearson->Draw();
  Neyman->Draw("SAME");
  Likelihood->Draw("SAME");

}