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

#define Iter 10000
#define sigma 3.435188


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;
    double LimSens=LimitiPoisson(Nbkg,Nbkg,CL); //limite di sensibilità calcolato con il metodo classico
    double LimSensCls=Zech(Nbkg,Nbkg,CL);       //limite di sensibilità calcolato con il metodo Cls
    double media=0;
    int contatorePoisson=0;
    int contatoreCLs=0;
    int contatore1sPoisson=0;
    int contatore1sCLs=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",50,-5.,25.);
    //TH1F* hLimSens = new TH1F("hLimSens","Limite di sensibilità",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",30,-20.,30.,50,-20.,30.);


    double mm1s=Nbkg-sigma;
    double mp1s=Nbkg+sigma;

    printf("Misura meno una sigma %f\n",mm1s);
    printf("Misura piu' una sigma %f\n",mp1s);

    double Limmm1s=LimitiPoisson(mm1s,Nbkg,CL);
    double Limmp1s=LimitiPoisson(mp1s,Nbkg,CL);

    double CLsmm1s=Zech(mm1s,Nbkg,CL);
    double CLsmp1s=Zech(mp1s,Nbkg,CL);


    printf("Limite di sensibilita' con il metodo Classico %f\n",LimSens);
    printf("Limite classico della misura meno una sigma %f\n",Limmm1s);
    printf("Limite classico della misura piu' una sigma %f\n",Limmp1s);

    printf("Limite di sensibilita' con il metodo CLs %f\n",LimSensCls);
    printf("Limite CLs della misura meno una sigma %f\n",CLsmm1s);
    printf("Limite CLs della misura piu' una sigma %f\n",CLsmp1s);   


  	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);

        media+=LimCLs;

        if(LimPoisson<0) contatorePoisson++;
        if(LimCLs<0) contatoreCLs++;
        if(LimPoisson < Limmp1s && LimPoisson > Limmm1s) contatore1sPoisson++;
        if(LimCLs < CLsmp1s && LimCLs > CLsmm1s) contatore1sCLs++;


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

    media/=Iter;

  	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");


   // double sigma=hmisure->GetStdDev();
   // printf("Sigma %f\n",sigma);

    double mediaPoisson=hLimPoisson->GetMean();
    double sigmaPoisson=hLimPoisson->GetStdDev();
    double mediaCLs=hLimCLs->GetMean();
    double sigmaCLs=hLimCLs->GetStdDev();
    
    //printf("Media dei limiti CLs %f\n",media);

    printf("\nMedia hPoisson +1 sigma %f ",mediaPoisson+sigmaPoisson);
    printf("\nMedia hPoisson -1 sigma %f ",mediaPoisson-sigmaPoisson);
    printf("\nMedia hCLs +1 sigma %f ",mediaCLs+sigmaCLs);
    printf("\nMedia hCLs -1 sigma %f ",mediaCLs-sigmaCLs);

    printf("\nFrequenza limiti Poisson minori di 0   %f",(double) contatorePoisson/Iter);
    printf("\nFrequenza limiti CLs minori di 0   %f",(double) contatoreCLs/Iter);
    printf("\nFrequenza limiti Poisson compresi in +/- 1 sigma   %f",(double) contatore1sPoisson/Iter);
    printf("\nFrequenza limiti CLs compresi in +/- 1 sigma   %f",(double) contatore1sCLs/Iter);




  	return;
}