#include "TMinuit.h"
#include "TRandom.h"
#include "TH1F.h"
#include "TMath.h"
#include <iostream>
#include <stdlib.h>
#include "TGraphErrors.h"

using namespace TMath;

const double TWOPI=TMath::TwoPi();
const double RANDM=RAND_MAX;

// Neyman construction
TH2F *M_Neyman = new TH2F("M_Neyman","",20,0,20,2000,0.005,20.005);

//Signal and background model
float b = 3.00;

//---------------------------------------------------
double poisson(double nbar, double n)
//---------------------------------------------------
{
  double u = 1;
  for( int i=0 ; i<n ; i++ )
  {
    u *= nbar/(i+1);
  }
  return u*exp(-nbar);
}


//----------------------------------------------------------------------
void Neyman()
//----------------------------------------------------------------------
{

  double mu=0;
  Double_t P[40];
  Double_t Pordered[40];
  Double_t Pbest[40];
  Double_t R[40];
  Int_t idx[40];
  
  //Scan values of mu
  for (int iMu =0; iMu != 2000; iMu++)
    {
      mu=iMu*0.01;
      
      for (Int_t i=0; i!=40; i++) {
	double obs = i;
	P[i] = poisson(mu+b,obs);
	Pbest[i] = poisson(max(0.,obs-b)+b,obs);
	R[i]=P[i]/Pbest[i];
	//idx[i]=39-i;  // Use this definition of idx for standard limit 
      }
      TMath::Sort(40,R,idx); // Replace R by P for Central Confidence Interval      
      for (int i=0; i!=40; i++) {
	//cout << mu << "  " << b << "  " << i << "  " << P[i] << "  " << max(0.,i-b) << "  " << Pbest[i] << "  " << R[i] << "  " << R[idx[i]] << endl;
	Pordered[i]=P[idx[i]];
      }
      
      double a[5]={1,2,3,4,5};
      int index[5];

      float alpha=0.;
      for (int i=0; i!=40; i++) {
	if (alpha<=0.90) {
	  M_Neyman->Fill(idx[i],mu);
	  alpha+=Pordered[i];
	}
      }
      
      
    }
 
 // Compute the Feldman Cousins limit
  // Plot the result 
  
  gROOT->SetStyle("Plain");
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetPadColor(0);
  gStyle->SetLineWidth(1.);
  gStyle->SetCanvasDefH(596);
  gStyle->SetCanvasDefW(596);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetStatBorderSize(0);
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  
  c0 = new TCanvas("c0", "", 0, 0, 600, 450);
  c0->SetLeftMargin(0.12);
  c0->SetRightMargin(0.12);
  c0->SetBottomMargin(0.2);
 
  gStyle->SetPalette(0);
  M_Neyman->SetMaximum(2.);
  M_Neyman->SetMinimum(0.99);
  M_Neyman->Draw("colz");
}