#include "TH1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "plot_tools.h"

void gauss_dist()
{
    plot_nice();
    gStyle->SetPalette(1); 

    // We need to create an object that is a random number generator
    TRandom* random = new TRandom();


    TH1F *h1  = new TH1F("h1","Gaussian distribution",2001,-0.5,2000.5);
    TH1F *h_dist  = new TH1F("h_dist","Gaussian distribution",2001,-0.5,2000.5);
    

    int n_throw = 500000; // number of events thrown

    float n_background = 1000;

    float n_on;
    float n_off;

    float mean = n_background;
    float sigma = sqrt(mean);
    float n_sig = 50;

	float mean_on  = mean + n_sig;
	float sigma_on = sqrt(mean_on);
	for(int i=0;i<n_throw;i++)
	{
	    n_on  = random->Gaus(mean_on, sigma_on);
	    n_off = random->Gaus(mean, sigma);
	    
	    h1->Fill(n_off);
	    h_dist->Fill(n_on-n_off);	    
	}
	// now find 90/95% C.L upper limit
        float sum=0;
	float sum_prev=0;
	float n5=0;
	float n10=0;
	float n90=0;
	float n95=0;
	for(int i=0;i<1001;i++)
	{
	    sum_prev=sum;
	    sum+=h_dist->GetBinContent(i);
	    printf("i=%d sum=%f\n",i,sum);
	    if((sum/n_throw)>0.05 && (sum_prev/n_throw)<0.05) 
	    {
		n5=i-500;
	    }
	    if((sum/n_throw)>0.10 && (sum_prev/n_throw)<0.10) 
	    {
		n10=i-500;
	    }
	    if((sum/n_throw)>0.90 && (sum_prev/n_throw)<0.90) 
	    {
		n90=i-500;
	    }
	    if((sum/n_throw)>0.95 && (sum_prev/n_throw)<0.95) 
	    {
		n95=i-500;
		i=1001;
	    }
	}
	//h_dist->Reset();
    }
    TCanvas *c1 = new TCanvas("c1","c1",200,300,500,500);
    c1->Divide(2,1);
    h1->SetTitle("");
    h1->GetXaxis()->SetTitle("Events");
    h1->GetYaxis()->SetTitle("rel. probability");
    h1->Draw();
    h_dist->SetTitle("");
    h_dist->GetXaxis()->SetTitle("Events");
    h_dist->GetYaxis()->SetTitle("rel. probability");
    h_dist->Draw();
    
}