// teste de chi2
//  ref Cowan, cap 4
//  jtmn, Set 2003

void chi2test()
{
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);
  TRandom *r3 = new TRandom3();
  r3->SetSeed(0);
  
 // simule um histograma de dados com um pico
  TCanvas *cv = new TCanvas("cv"," Hist  ", 600,600);
  Int_t Nbins = 20;
  Double_t nu_back = 1.6;
  TH1F *h  = new TH1F("hexp"," ", Nbins,0.0,20.0);
  // primeiro o background, nu_back eventos por bin em todos os bins
  for (Int_t i=1;i<=Nbins;i++){
   Double_t x = r3->Poisson(nu_back);
   h->SetBinContent(i,x);
  }
  // agora o pico
  h->SetBinContent(10,6);
  h->SetBinContent(11,5);
  // mostre o histograma
  h->Draw();
  // pegue o no. de entradas no hist.
  Int_t Nentr = h->Integral();
  cout << "Nentr = " << Nentr << endl;
  cv->Update();

  // calcule o chisq
  Double_t ExpChis = 0.0;
  for (Int_t i=1;i<=Nbins;i++){
    Double_t x=h->GetBinContent(i);
    ExpChis = ExpChis + (x-nu_back)*(x-nu_back)/nu_back;
  }
  cout << "ExpChis = " << ExpChis << endl;
 
 //faca Nexp exper. com a hipotese que so tem back
  TCanvas *cv2 = new TCanvas("cv2"," Hist Chis  ", 600,600);
  Int_t Nexp = 10000;
  Int_t NbinsChis = 60;
  TH1F *hchis  = new TH1F("hchis"," ", NbinsChis,0.0,60.0);
  for (Int_t i=0; i<Nexp; i++){
    Double_t chis = 0.0;
    for (Int_t j=1; j<Nbins; j++){
       Double_t x = r3->Poisson(nu_back);
       chis = chis + (x-nu_back)*(x-nu_back)/nu_back;
    }
    // encha hist com peso
    hchis->Fill(chis, 1.0/Double_t(Nexp));
  }
  // encontre posicao e bin que correspondem ao ExpChis
  hchis->Draw();
  TAxis *xaxis = hchis->GetXaxis();
  Int_t ilow = xaxis->FindFixBin(ExpChis);
  cout << " ilow = " << ilow << endl;
  // calcule Pval
  Double_t Pval = hchis->Integral(ilow,NbinsChis);
  cout << " Pval = " << Pval << endl;

  //vefique a mesma coisa supondo uma distribuicao chis matematica
  // NDOF = Nbins
  Double_t PvalTeo = TMath::Prob(ExpChis,Nbins);
  cout << " PvalTeo = " << PvalTeo << endl; 
  
}