// *********************************************************************
// 					Test of TLimit
//	compiled with : g++ -o prog prog.cpp `root-config --cflags --libs`
// *********************************************************************

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>

#include "TRandom.h"
#include "TMath.h"
#include "TH1F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraph2D.h"
#include "TLimitDataSource.h"
#include "TConfidenceLevel.h"
#include "TLimit.h"

using namespace std;

// *********************************************************************

#define XMIN_RANGE -7
#define XMAX_RANGE 7

#define NBBINS 50

#define NEVENTS_BKG 100000
#define RESCALE_BKG 10
#define NEVENTS_SIG 150

// *********************************************************************

int main()
{

	// Init stuff

		TRandom* gen = new TRandom();
		gen->SetSeed();

		TH1F* signal = new TH1F("signal","Signal histo"			   ,NBBINS,XMIN_RANGE,XMAX_RANGE);
		TH1F* background = new TH1F("background","Background histo",NBBINS,XMIN_RANGE,XMAX_RANGE);
		TH1F* data = new TH1F("data","Data histo"				   ,NBBINS,XMIN_RANGE,XMAX_RANGE);

		// background->Sumw2();
		// signal->Sumw2();

	// Generate gaussians for signal and bkg

		for (int i = 0 ; i < NEVENTS_BKG ; i++)
			background->Fill(gen->Gaus(0,3));

		for (int i = 0 ; i < NEVENTS_SIG ; i++)
			signal->Fill(gen->Gaus(-4,0.3));

		background->Scale(1.0 /(float) RESCALE_BKG);

	// Compute data as signal + bkg (or just bkg)

		data->Add(signal,background);
		//data->Add(background);

	// Plot histos

		TFile* output = new TFile("histos.root","RECREATE");
		TCanvas* canvas = new TCanvas("canvas","canvas",800,600);
		canvas->Divide(2,2);

		canvas->cd(1);
		background->Draw();
		canvas->cd(2);
		signal->Draw();
		canvas->cd(4);
		data->Draw();

		canvas->Write();
		output->Close();

	// Computing TConfidence Level to TLimit

		TLimitDataSource* mydatasource = new TLimitDataSource(signal,background,data);
		TConfidenceLevel* myconfidence = TLimit::ComputeLimit(mydatasource,500000);

	// Display outputs

		cout << "  Clb   : " << myconfidence->CLb()  					<< endl;
		cout << "  Clsb   : " << myconfidence->CLsb() 					<< endl;
		cout << "  CLs    : " << myconfidence->CLs()  					<< endl;
		cout << endl;
		cout << "  ExpCLs_b  : " << myconfidence->GetExpectedCLs_b()  		<< endl;
		cout << "  ExpCLs_b +1s  : " << myconfidence->GetExpectedCLs_b(1)  	<< endl;
		cout << "  ExpCLs_b -1s  : " << myconfidence->GetExpectedCLs_b(-1) 	<< endl;
		cout << endl;
		cout << "  ExpClb_b   : " << myconfidence->GetExpectedCLb_b()  	<< endl;
		cout << "  ExpClb_sb   : " << myconfidence->GetExpectedCLb_sb() << endl;

	// Draw -2ln(Q)

		TFile* output2 = new TFile("-2ln_Q.root","RECREATE");
		TCanvas* canvas2 = new TCanvas("canvas2","canvas2",800,600);

			myconfidence->Draw();

		canvas2->Write();
		output2->Close();

	// Exit

	delete myconfidence;
	delete mydatasource;

	return 0;

}