#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>

// Simple routine needed to calculate the inverse error function
// -------------------------------------------------------------
Double_t ErfInverse (Double_t x);

// Macro that computes p-value and Z-value of N observed vs B predicted
// Poisson counts
// --------------------------------------------------------------------
void Poisson_prob_fix (double B, double N) {

  int maxN = N*3/2;
  if (N<20) maxN=2*N;
  TH1D * Pois = new TH1D ("Pois", "", maxN, -0.5, maxN-1);
  TH1D * PoisGt = new TH1D ("PoisGt", "", maxN, -0.5, maxN-1);

  double sum=0.;
  double fact=1.;
  for (int i=0; i<maxN; i++) {
    if (i>1) fact*=i;
    double poisson = exp(-B)*pow(B,i)/fact;
    if (i<N) sum+= poisson;
    Pois->SetBinContent(i+1,poisson);
    if (i>=N) PoisGt->SetBinContent(i+1,poisson);
  }
  double P=1-sum;
  double Z = sqrt(2) * ErfInverse(1-2*P);

  cout << "P of observing N=" << N << " or more events if B=" 
       << B << " : P= " << 1-sum << endl;
  cout << "This corresponds to " << Z << " sigma for a Gaussian one-tailed test." << endl;

  Pois->SetLineWidth(3);
  PoisGt->SetFillColor(kBlue);
  TCanvas* T = new TCanvas ("T","Poisson distribution", 500, 500);
  T->Divide(1,2);
  T->cd(1);
  Pois->Draw();
  PoisGt->Draw("SAME");
  T->cd(2);
  T->GetPad(2)->SetLogy();
  Pois->Draw();
  PoisGt->Draw("SAME");
}

Double_t ErfInverse(Double_t y) {
// returns the inverse error function obtained
// y must be <-1<y<1

Int_t kMaxit = 50;
Double_t kEps = 1e-14;
Double_t kConst = 0.8862269254527579; // sqrt(pi)/2.0

if(TMath::Abs(y) <= kEps) return kConst*y;

// Newton iterations
Double_t erfi, derfi, Y0,Y1,DY0,DY1;
if(TMath::Abs(y) < 1.0) {
erfi = kConst*TMath::Abs(y);
Y0 = TMath::Erf(0.9*erfi);
derfi = 0.1*erfi;
for (Int_t iter=0; iter<kMaxit; iter++) {
Y1 = 1. - TMath::Erfc(erfi);
DY1 = TMath::Abs(y) - Y1;
if (TMath::Abs(DY1) < kEps) {if (y < 0) return -erfi; else return erfi;}
DY0 = Y1 - Y0;
derfi *= DY1/DY0;
Y0 = Y1;
erfi += derfi;
if(TMath::Abs(derfi/erfi) < kEps) {if (y < 0) return -erfi; else return erfi;}
}
}
 return 0; //did not converge
}