#define SimPipeline_cxx
#include "SimPipeline.h"

/*************************************************
Default constructor
*************************************************/
SimPipeline::SimPipeline()
{
	luminosity = 1;
}

/*************************************************
Default destructor
*************************************************/
SimPipeline::~SimPipeline()
{

}

/*************************************************
This adds a background sample to the collection.
*************************************************/
void SimPipeline::AddBackground(Sample* bg)
{
	backgrounds.add(bg);

	luminosity = bg->GetLuminosity();			//Hopefully all your backgrounds have the same luminosity!
}

/*************************************************
This adds a signal sample to the collection.
*************************************************/
void SimPipeline::AddSignal(Sample* sig)
{
	signals.add(sig);
}

/*************************************************
This computes the signal limits for all the 
backgrounds.
*************************************************/
void SimPipeline::ComputeSignalLimits()
{
	//Int_t n_background = 0;			//Total number of background events

	//First, we add the tables for the backgrounds
	Int_t metbins = backgrounds[0]->GetMETBins();		//Hopefully all the backgrounds are binned the same way
	Int_t htbins = backgrounds[0]->GetHTBins();

	//cout<<"Good so far....."<<endl;
	//Int_t a;
	//cin>>a;

	total_background[0].ResizeTo(1, metbins);			//Set dimensions for total background tables
	total_background[1].ResizeTo(htbins, metbins);
	total_background[2].ResizeTo(htbins, metbins);
	total_background[3].ResizeTo(htbins, metbins);
	for (Int_t i = 0; i < backgrounds.size(); i++)		//Sum tables
	{
		total_background[0] += backgrounds[i]->table[0];
		total_background[1] += backgrounds[i]->table[1];
		total_background[2] += backgrounds[i]->table[2];
		total_background[3] += backgrounds[i]->table[3];

		//n_background += backgrounds[i].GetNumOfEvents();
	}

	//Double_t sl_sys = 0.5*n_background;

	//Int_t metbins = 5;
	//Int_t htbins = 5;

	Double_t b;

	//Now we build the signal limit tables
	signal_limits[0].ResizeTo(1, metbins);			//Set dimensions for signal limit tables
	for (Int_t row = 0; row < 1; row++)
	{
		for (Int_t col = 0; col < metbins; col++)
		{
			b = luminosity*total_background[0][row][col];
			if (b < 10)									//Poissonian case
			{
				signal_limits[0][row][col] = sqrt( pow(CalcSLStat(b), 2) + pow(SL_SYS*b, 2) );
				//cout<<"Table 0, row "<<row<<", col "<<col<<": background "<<b/4<<", statistical "<<CalcSLStat(b)/4<<", systematic "<<SL_SYS*b/4<<endl;
				//Int_t x;
				//cin>>x;
			}
			else										//Gaussian limit
			{
				signal_limits[0][row][col] = sqrt( b + pow(0.5*b, 2) );
				//cout<<"Table 0, row "<<row<<", col "<<col<<": background "<<b/4<<", statistical "<<sqrt(b)/4<<", systematic "<<SL_SYS*b/4<<endl;
				//Int_t x;
				//cin>>x;
			}
		}
	}
	signal_limits[1].ResizeTo(htbins, metbins);
	signal_limits[2].ResizeTo(htbins, metbins);
	signal_limits[3].ResizeTo(htbins, metbins);
	for (Int_t i = 1; i < 4; i++)
	{
		for (Int_t row = 0; row < htbins; row++)
		{
			for (Int_t col = 0; col < metbins; col++)
			{
				b =	luminosity*total_background[i][row][col];
				if (b < 10)								//Poissonian case
				{
					signal_limits[i][row][col] = sqrt( pow(CalcSLStat(b), 2) + pow(SL_SYS*b, 2) );
					//cout<<"Table "<<i<<", row "<<row<<", col "<<col<<": background "<<b/4<<", statistical "<<CalcSLStat(b)/4<<", systematic "<<SL_SYS*b/4<<endl;
					//Int_t x;
					//cin>>x;
				}
				else									//Gaussian limit
				{
					signal_limits[i][row][col] = sqrt( b + pow(0.5*b, 2) );
					//cout<<"Table "<<i<<", row "<<row<<", col "<<col<<": background "<<b/4<<", statistical "<<sqrt(b)/4<<", systematic "<<SL_SYS*b/4<<endl;
					//Int_t x;
					//cin>>x;
				}	
			}
		}
	}
}

/*************************************************
This uses the signal limit tables to determine 
whether the signal samples are excluded.
*************************************************/
void SimPipeline::GetSignalExclusions()
{
	Double_t chi_square;
	Double_t significance;
	Double_t dof;

	for (Int_t j = 0; j < signals.size(); j++)			//Loop through each signal sample
	{
		chi_square = 0;
		dof = 0;
		for (Int_t i = 0; i < 4; i++)					//Loop through each table
		{
			for (Int_t row = 0; row < total_background[i].GetNrows(); row++)
			{
				for (Int_t col = 0; col < total_background[i].GetNcols(); col++)
				{
					if (signal_limits[i][row][col] != 0)
					{
						significance = (luminosity*signals[j]->table[i][row][col])/signal_limits[i][row][col];
						if (significance > 0.5)
						{
							cout<<"Table "<<i<<", row "<<row<<", col "<<col<<", significance "<<significance<<endl;
							dof++;
							chi_square += pow(significance, 2);
						}
					}
				}
			}
		}
		cout<<"Degrees of freedom: "<<dof<<endl;
		if (dof != 0)
		{
			chi_square = chi_square;
			cout<<"Degrees of freedom: "<<dof<<endl;
			cout<<"Minimum chi^2 (non-normalized) for 95%: "<<TMath::ChisquareQuantile(0.95, dof)<<endl;
			cout<<"Actual chi^2 (non-normalized): "<<chi_square<<endl;
			if (TMath::ChisquareQuantile(0.95, dof) < chi_square)
				cout<<"Signal sample "<<j<<" is excluded."<<endl;
			else
				cout<<"Signal sample "<<j<<" is NOT excluded (chi^2 too small)."<<endl;
		}
		else
			cout<<"Signal sample "<<j<<" is NOT excluded (no good bins)."<<endl;
	}
}

/*************************************************
This takes a background cross section and 
calculates the statistical signal limit.
*************************************************/
Double_t SimPipeline::CalcSLStat(Double_t b)
{
	Double_t sl_stat = 0;		//Statistical signal limit

	for (Int_t nm = 0; nm < 41; nm++)
	{
		if ( (CalcMu(nm) - b) > 0 )
			sl_stat += ( (CalcMu(nm) - b)*exp(-b)*pow(b, nm)/(TMath::Factorial(nm)) );
	}
	
	return sl_stat;
}

/*************************************************
This returns mu_excl for a given N_m.
*************************************************/
Double_t SimPipeline::CalcMu(Int_t nm)
{
	if (nm == 0) return 1.840;
	if (nm == 1) return 3.290;
	if (nm == 2) return 4.630;
	if (nm == 3) return 5.910;
	if (nm == 4) return 7.150;
	if (nm == 5) return 8.370;
	if (nm == 6) return 9.570;
	if (nm == 7) return 10.760;
	if (nm == 8) return 11.930;
	if (nm == 9) return 13.090;
	if (nm == 10) return 14.250;
	if (nm == 11) return 15.400;
	if (nm == 12) return 16.540;
	if (nm == 13) return 17.680;
	if (nm == 14) return 18.810;
	if (nm == 15) return 19.940;
	if (nm == 16) return 21.060;
	if (nm == 17) return 22.180;
	if (nm == 18) return 23.300;
	if (nm == 19) return 24.410;
	if (nm == 20) return 25.520;
	if (nm == 21) return 26.630;
	if (nm == 22) return 27.740;
	if (nm == 23) return 28.840;
	if (nm == 24) return 29.940;
	if (nm == 25) return 31.040;
	if (nm == 26) return 32.140;
	if (nm == 27) return 33.230;
	if (nm == 28) return 34.330;
	if (nm == 29) return 35.420;
	if (nm == 30) return 36.510;
	if (nm == 31) return 37.600;
	if (nm == 32) return 38.690;
	if (nm == 33) return 39.770;
	if (nm == 34) return 40.860;
	if (nm == 35) return 41.940;
	if (nm == 36) return 43.020;
	if (nm == 37) return 44.100;
	if (nm == 38) return 45.190;
	if (nm == 39) return 46.260;
	if (nm == 40) return 47.340;

	return 0;
}