#include "TH1F.h"
#include <string>
void PoissonianGraph()
{
  // Calcolare l'upper limit nel caso di
  // n eventi rivelati contro un fondo atteso di b eventi

  std::string hname;
  std::string hnamec;
  TH1F* hPoiss[5];
  TH1F* hPoissC[5];
  for (int i=0; i<5; ++i)
    {
      hname  = "hPoiss"+std::to_string(i);
      hnamec = "hPoissC"+std::to_string(i);
      hPoiss[i]  = new TH1F(hname.c_str(),"Poissonian distribution",31,-0.5,30.5);
      hPoissC[i] = new TH1F(hnamec.c_str(),"Cumulative Poissonian distribution",31,-0.5,30.5);
    }
  
  
  double poisson=0.;
  double cpoisson=0.;
  double mu0;
  double bkg = 4;
  double signal = 0;
  for (int j=0; j<5; ++j)
    {
      signal = 1+j*5;
      mu0 = bkg+signal;
      std::cout<<" bkg, signal, mu_total = "<<bkg<<" "<<signal<<" "<<mu0<<std::endl;
      for (int i=0; i<31; ++i)
	{
	  poisson=TMath::PoissonI(i,mu0);
	  cpoisson=ROOT::Math::poisson_cdf(i,mu0);
	  std::cout<<"P("
		   <<i<<"|mu="
		   <<mu0<<") = "
		   <<poisson<<" P(x<"
		   <<i<<") = "
		   <<cpoisson
		   <<std::endl;
	  hPoiss[j]->SetBinContent( i+1, poisson);
	  hPoissC[j]->SetBinContent(i+1,cpoisson);
	}
    }


  hPoissC[0]->Draw("L");
  for (int i=0; i<4; ++i)
    {
      hPoiss[i] ->SetLineWidth(5);
      hPoissC[i] ->SetLineWidth(2);
      hPoiss[i] ->SetLineColor(2+i*2);
      hPoissC[i] ->SetLineColor(2+i*2);
      hPoiss[i] ->Draw("SAME");
      hPoissC[i] ->Draw("SAMEL");
    }

  
}