#include "Zech.C"
#include "LimitiPoisson.C"
#include "FeldmanCousins.C"
#include "TRandom.h"

#define Iter 10000


void LimiteAtteso(double Nfondo = 12)
{
	double Nbkg = Nfondo;
  	double CL = 0.9;
  	double dato=0;
  	TRandom3 *R = new TRandom3();
  	double segnale=0;
  	double LimPoisson=0;
  	double LimCLs=0;


  	TH1F* hmisure = new TH1F("hmisure","Misure",100,-6.,25.);
  	TH1F* hsegnale = new TH1F("hsegnale","Segnale",100,-20.,20.);
  	TH1F* hLimPoisson = new TH1F("hLimPoisson","Limite Classico",100,-15.,25.);
  	TH1F* hLimCLs = new TH1F("hLimCLs","Limite CLs-Zech",100,-5.,25.);

  	TH2F* hCLsVsStd = new TH2F("hCLsVsStd","CLs vs Standard Poisson",50,-20.,30.,50,-20.,30.);
  	TH2F* hCLsVsSign = new TH2F("hCLsVsSign","CLs vs Segnale",50,-20.,30.,50,-20.,30.);
  	TH2F* hStdVsSign = new TH2F("hStdVsSign","Standard Poisson vs Segnale",50,-20.,30.,50,-20.,30.);

  	for(int i=0;i<Iter;i++)
  	{
  		dato=R->Poisson(Nbkg);
  		hmisure->Fill(dato);

  		segnale = dato - Nbkg;
  		hsegnale->Fill(segnale);

  		LimPoisson=LimitiPoisson(dato,Nbkg,CL);
  		LimCLs=Zech(dato,Nbkg,CL);
  		hLimPoisson->Fill(LimPoisson);
  		hLimCLs->Fill(LimCLs);

  		hCLsVsStd->Fill(LimCLs,LimPoisson);
  		hStdVsSign->Fill(segnale,LimPoisson);
  		hCLsVsSign->Fill(segnale,LimCLs);
  	}

  	TCanvas* c = new TCanvas("c","c",1100,800);
    c->Divide(2,2);
    c->cd(1);
   	hmisure->Draw();
    c->cd(2);
    hsegnale->Draw();
    c->cd(3);
    hLimPoisson->Draw();
    c->cd(4);
    hLimCLs->Draw();

    TCanvas* c1 = new TCanvas("c1","c1",800,600);
    c1->Divide(2,2);
    c1->cd(1);
    hCLsVsStd->Draw("colz");
    c1->cd(2);
    hStdVsSign->Draw("colz");
    c1->cd(3);
    hCLsVsSign->Draw("colz");

  	return;
}