#include <iostream>

#include <algorithm>
#include <fstream>
#include <iterator>

//#include <array>

#include <math.h>
#include <stdio.h>    
#include <stdlib.h>

#include "determine_offset.h"

using namespace std;



int Q2_2p4_EB_4p2[] = {47048, 47049, 47050, 47051, 47052};
int Q2_6p5_EB_5p2[] = {47354, 47355, 47356};
int Q2_4p4_EB_3p7[] = {47575, 47576, 47577, 47578, 47579, 47580, 47581};
int Q2_5p4_EB_4p7[] = {47679, 47680, 47681, 47682, 47683, 47684, 47685}; 


vector<int> Q2_2p4_EB_4p2_vec (Q2_2p4_EB_4p2, Q2_2p4_EB_4p2 + sizeof(Q2_2p4_EB_4p2)/sizeof(Q2_2p4_EB_4p2[0]));
vector<int> Q2_6p5_EB_5p2_vec (Q2_6p5_EB_5p2, Q2_6p5_EB_5p2 + sizeof(Q2_6p5_EB_5p2)/sizeof(Q2_6p5_EB_5p2[0]));
vector<int> Q2_4p4_EB_3p7_vec (Q2_4p4_EB_3p7, Q2_4p4_EB_3p7 + sizeof(Q2_4p4_EB_3p7)/sizeof(Q2_4p4_EB_3p7[0]));
vector<int> Q2_5p4_EB_4p7_vec (Q2_5p4_EB_4p7, Q2_5p4_EB_4p7 + sizeof(Q2_5p4_EB_4p7)/sizeof(Q2_5p4_EB_4p7[0])); 


/*--------------------------------------------------*/
// Main function

int main() {


	Init();

 	TCut sos_e_selection = "scer_npe > 1 && ssshsum > 0.8 && ssshsum < 1.2";
 
    Add_Root_File(Q2_2p4_EB_4p2);
  	Add_Root_File(Q2_6p5_EB_5p2);
  	Add_Root_File(Q2_4p4_EB_3p7);
  	Add_Root_File(Q2_5p4_EB_4p7);
 	
 	globchain_total->Add(globchain1);
 	globchain_total->Add(globchain2);
 	globchain_total->Add(globchain3);
 	globchain_total->Add(globchain4);
 
  	Coin_time_loop(Q2_2p4_EB_4p2_vec);
  	Coin_time_loop(Q2_6p5_EB_5p2_vec);
  	Coin_time_loop(Q2_4p4_EB_3p7_vec);
  	Coin_time_loop(Q2_5p4_EB_4p7_vec);
 
 	Write_overall();

	return 0;
}


/*--------------------------------------------------*/
/// Initialization
void Init() {

	globchain1 		=  new TChain("h9500");
	globchain2 		=  new TChain("h9500");
	globchain3 		=  new TChain("h9500");
	globchain4 		=  new TChain("h9500");
	globchain_total =  new TChain("h9500");



	/*--------------------------------------------------*/
	/// Output Director

	file_out = new TFile("plot_out.root", "RECREATE");

	test_dir               = file_out->mkdir("test_dir");
	plot_dir_Q2_2p4_EB_4p2 = file_out->mkdir("Q2_2p4_EB_4p2");
	plot_dir_Q2_6p5_EB_5p2 = file_out->mkdir("Q2_6p5_EB_5p2");
	plot_dir_Q2_4p4_EB_3p7 = file_out->mkdir("Q2_4p4_EB_3p7");
	plot_dir_Q2_5p4_EB_4p7 = file_out->mkdir("Q2_5p4_EB_4p7");



	/*--------------------------------------------------*/
	/// Combining Plot for each kinematic settings

	Q2_2p4_EB_4p2_combine  = new TH2F("Q2_2p4_EB_4p2_combine", "", 441, -11, 11, 501, -1.5, 1);
	Q2_6p5_EB_5p2_combine  = new TH2F("Q2_6p5_EB_5p2_combine", "", 441, -11, 11, 501, -1.5, 1);
	Q2_4p4_EB_3p7_combine  = new TH2F("Q2_4p4_EB_3p7_combine", "", 441, -11, 11, 501, -1.5, 1);
	Q2_5p4_EB_4p7_combine  = new TH2F("Q2_5p4_EB_4p7_combine", "", 441, -11, 11, 501, -1.5, 1);
	total_hbeta_cointime   = new TH2F("Overall", "", 441, -11, 11, 501, -1.5, 1);


	beta_cut                   = 0.1;
	global_centered_hsbeta_cut = -0.15;
	gradi                      = 9;	


}

/*--------------------------------------------------*/
// Adding all the root file togather

void Add_Root_File(Int_t* grp) {

//	cout << "  asdasddd"  << Q2_2p4_EB_4p2 << endl;

//	cout << "  asdas"  << grp << endl;

	TString path("../data/");
    TString filename;  
//	cout << grp[3]  << "  " << sizeof(grp) << endl;

	for(int i=0; i < sizeof(grp); i++) {
		//cout << i << "     "<< grp[i] << endl;
	    filename.Form("coin%i.root",grp[i]);
		filename = path + filename;
		// cout << filename << endl;
		

		ifstream my_file(filename);

		if (my_file.good()) {

			if (grp == Q2_2p4_EB_4p2) {
				globchain1->Add(filename);
			} else if (grp == Q2_6p5_EB_5p2) {
				globchain2->Add(filename);
			} else if (grp == Q2_4p4_EB_3p7) {
				globchain3->Add(filename);
			} else if (grp == Q2_5p4_EB_4p7) {
				globchain4->Add(filename);
			} else {
				cout << "*************************" << endl;
				cout << "There is something wrong!" << endl;
				cout << "*************************" << endl;
			}
		}
	}

}






/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// Computing the HMS Cerenkov (h_cer) efficiency

void Cer_plot(TChain*chain, TString plot_title, TCut data_cut, Float_t threshhold, Float_t bin_min, Float_t bin_max, Int_t bin_num) {


	TH1F *myhist = new TH1F("myhist", plot_title, bin_num, bin_min, bin_max);

 	TCanvas *c5 = new TCanvas();
 	chain->Draw("hcer_npe>>myhist", data_cut); 
	myhist->Draw("hist");

	c5->SetLogy();

	Int_t bin_thre = (threshhold - bin_min) / ( (fabs(bin_min)+fabs(bin_max))/(float) bin_num) ;
	Int_t lower = myhist->Integral(0, bin_thre);
	Int_t upper = myhist->Integral(bin_thre+1, -1);

	float eff = (float) lower / (float) (upper + lower);

	cout << threshhold << "    "<< bin_thre << "     "<< lower << "    " << upper << "     " << eff << endl; 
	
}


/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// Computing the HMS Beta efficiency

void Beta_plot(TChain*chain, TString plot_title, TCut data_cut, Float_t threshhold, Float_t bin_min, Float_t bin_max, Int_t bin_num) {

	TH1F *myhist = new TH1F("myhist", plot_title, bin_num, bin_min, bin_max);

 	TCanvas *c5 = new TCanvas();
 	chain->Draw("hsbeta>>myhist", data_cut); 
	myhist->Draw("hist");

	c5->SetLogy();

	Int_t bin_thre = (threshhold - bin_min) / ( (fabs(bin_min)+fabs(bin_max))/(float) bin_num) ;
	Int_t lower = myhist->Integral(0, bin_thre);
	Int_t upper = myhist->Integral(bin_thre+1, -1);

	float eff = (float) upper / (float) (upper + lower);

	cout << threshhold << "    "<< bin_thre << "     "<< lower << "    " << upper << "     " << eff << endl; 

}


/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// Computing the HMS Aerogel Detector Efficiency


void Aero_plot(TChain*chain, TString plot_title, TCut data_cut, Float_t threshhold, Float_t bin_min, Float_t bin_max, Int_t bin_num) {


	TH1F *myhist = new TH1F("myhist", plot_title, bin_num, bin_min, bin_max);

 	TCanvas *c5 = new TCanvas();
 	chain->Draw("haero_su>>myhist", data_cut); 
	myhist->Draw("hist");

	c5->SetLogy();

	Int_t bin_thre = (threshhold - bin_min) / ( (fabs(bin_min)+fabs(bin_max))/(float) bin_num) ;

//	Int_t bin_thre = myhist->GetBin(10);

	Int_t lower = myhist->Integral(0, bin_thre);
	Int_t upper = myhist->Integral(bin_thre+1, -1);

	float eff = (float) lower / (float) (upper + lower);

	cout << threshhold << "    "<< bin_thre << "     "<< lower << "    " << upper << "     " << eff << endl; 
	
}




/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// Some Checking Functions


// void file_open(TString filename) {
// 
// //	cout << filename << endl;
// 
// 
// //	TFile*rootfile  = new TFile(filename, "READ");
// //	globchain->Add(filename);
// //	globchain->Add("../data/coin47048.root");
// 
// 
// 	//TTree *data_tree = (TTree*)globchain->Get("h9500");	
// 	//data_tree-> Draw("hcer_npe");	
// 
// 
// //	globchain->Draw("hcer_npe");
// 	
// }






void check_plots(TChain*chain, TString parameter, TString log) {
	
 	TCanvas *c1 = new TCanvas();

	if(log == "log") {
		c1->SetLogy();
	}

 	chain->Draw(parameter); 
	
}


/*--------------------------------------------------*/


void plot_cointime_finebin(TChain*chain, TString plot_title, Float_t bin_min, Float_t bin_max, Int_t bin_num) {


	TH1F *myhist = new TH1F("myhist", plot_title, bin_num, bin_min, bin_max);

 	TCanvas *c5 = new TCanvas();
 	chain->Draw("cointime>>myhist"); 
	myhist->Draw("hist");

	c5->SetLogy();

	c5->Update();

}

void plot_beta_finebin(TChain*chain, TString plot_title, Float_t bin_min, Float_t bin_max, Int_t bin_num) {


	TH1F *myhist = new TH1F("myhist", plot_title, bin_num, bin_min, bin_max);

 	TCanvas *c5 = new TCanvas();
 	chain->Draw("hsbeta>>myhist"); 
	myhist->Draw("hist");

	c5->SetLogy();

	c5->Update();

}



void plot_beta_vs_cointime_finebin(TChain*chain, TString plot_title, Float_t x_bin_min, Float_t x_bin_max, Int_t x_bin_num, Float_t y_bin_min, Float_t y_bin_max, Int_t y_bin_num) {

 	TCanvas *c5 = new TCanvas();

	TH2F *myhist = new TH2F(plot_title, plot_title, x_bin_num, x_bin_min, x_bin_max, y_bin_num, y_bin_min, y_bin_max);
 	chain->Draw("hsbeta:cointime>>"+plot_title); 
	myhist->Draw("hist");
	c5->Update();

}


Int_t plot_beta_vs_cointime_finebin_cut(TChain*chain, TCut cut, TString plot_title, Float_t x_bin_min, Float_t x_bin_max, Int_t x_bin_num, Float_t y_bin_min, Float_t y_bin_max, Int_t y_bin_num) {

	TH2F *myhist = new TH2F("myhist_2D", plot_title, x_bin_num, x_bin_min, x_bin_max, y_bin_num, y_bin_min, y_bin_max);
 	chain->Draw("hsbeta:cointime>>myhist_2D", cut, "goff");

	return myhist->GetEntries();

}

void plot_beta_vs_cointime_finebin_pretty(TChain*chain, TString plot_title, Float_t x_bin_min, Float_t x_bin_max, Int_t x_bin_num, Float_t y_bin_min, Float_t y_bin_max, Int_t y_bin_num) {

 	TCanvas *c5 = new TCanvas();

	TH2F *myhist = new TH2F(plot_title, plot_title, x_bin_num, x_bin_min, x_bin_max, y_bin_num, y_bin_min, y_bin_max);
 	chain->Draw("hsbeta:cointime>>"+plot_title); 
	myhist->Draw("hist");

	c5->Update();


//	DrawGoodBox();
//	DrawRandomBox();

	c5->Update();

}



/*--------------------------------------------------*/
/// Coin time loop

//void Coin_time_loop(Int_t* grp) {

void Coin_time_loop(vector<int> grp) {
	
 	for(int i=0; i < grp.size(); i++) {
		cout << grp[i] << "    ";
 		Run_by_Run(grp[i], -18, 4, 441, -0.5, 2, 501 );		
	}

}



/*--------------------------------------------------*/
/// Correction at Run by Run basis

void Run_by_Run(Int_t run_num, Float_t x_bin_min, Float_t x_bin_max, Int_t x_bin_num, Float_t y_bin_min, Float_t y_bin_max, Int_t y_bin_num) {


	TString plot_title;
	TString path("../../data/");
	TString filename;

	plot_title.Form("%i", run_num);
	filename.Form("coin%i.root",run_num);
	filename = path + filename;

	TFile *file = TFile::Open(filename);
	TTree *tree = (TTree*)file->Get("h9500");

	TCanvas *c1 = new TCanvas();

	TH2F *myhist_2D = new TH2F(plot_title, plot_title, x_bin_num, x_bin_min, x_bin_max, y_bin_num, y_bin_min, y_bin_max);
	TH2F *get_center = new TH2F("get_center", "get center", x_bin_num, x_bin_min, x_bin_max, y_bin_num, y_bin_min, y_bin_max);

	tree->Draw("hsbeta:cointime>>" + plot_title, "hsbeta > 0.5 && hsbeta < 1.9 && cointime < 3.9 && cointime > -17.9");
//	tree->Draw("hsbeta:cointime>>" + plot_title, "hsbeta > 0.85 && hsbeta < 1.5 && cointime > -9 && cointime < -4");

	tree->Draw("hsbeta:cointime>>get_center", "hsbeta > 0.5 && hsbeta < 1.5 && cointime > -10 && cointime < -4", "goff");


 	float coin_time_mean = get_center->GetMean(1);
 	float hbeta_mean = get_center->GetMean(2);

	TString plotname;
	plotname.Form("plot_%i.eps",run_num);
	TString ps_plotname = "plot/" + plotname;


	/*--------------------------------------------------*/
	/// Re-adjust the center
	
// 	cout << myhist_2D->GetXaxis()->GetXmax() + ceil(fabs(coin_time_mean)) << "  " 
// 		 << myhist_2D->GetXaxis()->GetXmin() + ceil(fabs(coin_time_mean)) << "  " 
// 		 << myhist_2D->GetYaxis()->GetXmax() - 1  << "  " 
// 		 << myhist_2D->GetYaxis()->GetXmin() - 1  << "  " 
// 		 << endl; 
	
//	cout << "center " << coin_time_mean << "    " << hbeta_mean << endl; 


	float x_bin_max_cen;
	float x_bin_min_cen;


	if (coin_time_mean > -7) {
		x_bin_max_cen = myhist_2D->GetXaxis()->GetXmax() + ceil(fabs(coin_time_mean));
		x_bin_min_cen = myhist_2D->GetXaxis()->GetXmin() + ceil(fabs(coin_time_mean));
	} else {
		x_bin_max_cen = myhist_2D->GetXaxis()->GetXmax() + floor(fabs(coin_time_mean));
		x_bin_min_cen = myhist_2D->GetXaxis()->GetXmin() + floor(fabs(coin_time_mean));
	}

	float y_bin_max_cen = myhist_2D->GetYaxis()->GetXmax() - 1;
	float y_bin_min_cen = myhist_2D->GetYaxis()->GetXmin() - 1;


	TH2F *myhist_2D_cen = new TH2F(plot_title+"_centered", plot_title, x_bin_num, x_bin_min_cen, x_bin_max_cen, y_bin_num, y_bin_min_cen, y_bin_max_cen);

// 	cout << coin_time_mean << "  " << plot_title << "  " << x_bin_num << "  " << x_bin_min_cen << "  " 
// 		<< x_bin_max_cen << "  " << y_bin_num << "    " << y_bin_min_cen << "   "<< y_bin_max_cen << endl;

	for (int i = 0; i < myhist_2D->GetXaxis()->GetNbins(); i++) {
		for (int j = 0; j < myhist_2D->GetYaxis()->GetNbins(); j++) {
	 		myhist_2D_cen -> Fill(myhist_2D->GetXaxis()->GetBinCenter(i) + fabs(coin_time_mean), 
 					myhist_2D->GetYaxis()->GetBinCenter(j) - fabs(hbeta_mean), 
 					myhist_2D->GetBinContent(i,j));
    	}
   	}


	////*--------------------------------------------------*/
	/// It is very important to reset the StatBox after filling the histogram
	myhist_2D_cen->ResetStats();
    myhist_2D_cen->Draw("hist");

 	float coin_time_mean_cen = myhist_2D_cen->GetMean(1);
  	float hbeta_mean_cen = myhist_2D_cen->GetMean(2);


// 	cout << "before recenter : " << myhist_2D->GetEntries()  << "    " 
// 		 << myhist_2D_cen->GetEntries()<< endl;





// 	Int_t good_events_beta_cut = myhist_2D->Integral(myhist_2D->GetXaxis()->FindBin(coin_time_mean - 1), 
//  			myhist_2D->GetXaxis()->FindBin(coin_time_mean + 1), 
//  			myhist_2D->GetYaxis()->FindBin(hbeta_mean-beta_cut), 
//  			myhist_2D->GetYaxis()->FindBin(1.6)); 
//
// 	Int_t good_events_beta_cut_cen = myhist_2D_cen->Integral(
// 	    		myhist_2D_cen->GetXaxis()->FindBin(coin_time_mean_cen - 1), 
//  			myhist_2D_cen->GetXaxis()->FindBin(coin_time_mean_cen + 1), 
//  			myhist_2D_cen->GetYaxis()->FindBin(hbeta_mean_cen-beta_cut), 
//  			myhist_2D_cen->GetYaxis()->FindBin(1.1)); 
// 
// 	cout << "Compare : "<< good_events_beta_cut << "     " << good_events_beta_cut_cen << endl;
 


//	Calculate_efficiency(myhist_2D);


//	cout << "center " << coin_time_mean << "    " << hbeta_mean << endl; 

	Calculate_efficiency(myhist_2D_cen, global_centered_hsbeta_cut);

 	DrawGoodBox(coin_time_mean_cen, global_centered_hsbeta_cut);
 	DrawRandomBox(coin_time_mean_cen, global_centered_hsbeta_cut);


//    myhist_2D->Draw("hist");

 	c1->Print(ps_plotname);

	test_dir->cd();
 	c1->Write(plotname);

	//myhist_2D_cen->Write();


	Histogram_output(myhist_2D_cen);

//	exit(0);

	delete get_center;

}


/*--------------------------------------------------*/
/// Calculate the efficiency

Int_t Calculate_efficiency(TH2F* hist_target) {

	float x_mean = hist_target->GetMean(1);
	float y_mean = hist_target->GetMean(2);

	float y_max  = hist_target->GetYaxis()->GetXmax();
	float y_min  = hist_target->GetYaxis()->GetXmin();


//	float y_min  = hbeta_cut;

//	cout <<  y_max << "    " << y_min << endl;


	/*--------------------------------------------------*/
	Int_t good_events_beta_cut = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean-1), 
			hist_target->GetXaxis()->FindBin(x_mean + 1), 
			hist_target->GetYaxis()->FindBin(y_mean - beta_cut), 
			hist_target->GetYaxis()->FindBin(y_max)); 

	Int_t random_late_beta_cut = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean + 3), 
			hist_target->GetXaxis()->FindBin(x_mean + 9), 
			hist_target->GetYaxis()->FindBin(y_mean - beta_cut), 
			hist_target->GetYaxis()->FindBin(y_max));

	Int_t random_early_beta_cut = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean -9), 
			hist_target->GetXaxis()->FindBin(x_mean - 3 ), 
			hist_target->GetYaxis()->FindBin(y_mean - beta_cut), 
			hist_target->GetYaxis()->FindBin(y_max)); 



	/*--------------------------------------------------*/
	Int_t good_events_total = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean - 1), 
			hist_target->GetXaxis()->FindBin(x_mean + 1), 
			hist_target->GetYaxis()->FindBin(y_min), 
			hist_target->GetYaxis()->FindBin(y_max)); 

	Int_t random_late_total = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean +3), 
			hist_target->GetXaxis()->FindBin(x_mean + 9), 
			hist_target->GetYaxis()->FindBin(y_min), 
			hist_target->GetYaxis()->FindBin(y_max));

	Int_t random_early_total = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean - 9), 
			hist_target->GetXaxis()->FindBin(x_mean - 3), 
			hist_target->GetYaxis()->FindBin(y_min), 
			hist_target->GetYaxis()->FindBin(y_max)); 

	/*--------------------------------------------------*/
	Int_t total = hist_target->GetEntries();
	
	float good_event_random_delta_cut_subtracted = (float) good_events_beta_cut - 
		((float)random_early_beta_cut + (float)random_late_beta_cut)/2.;

	float good_event_random_subtracted = (float) good_events_total - 
		((float)random_early_total + (float)random_late_total)/2.;

	float beta_eff = good_event_random_delta_cut_subtracted / good_event_random_subtracted;

	float beta_eff_error = beta_eff * sqrt( pow(sqrt(good_event_random_subtracted)/good_event_random_subtracted, 2)  + 
	pow(sqrt(good_event_random_subtracted)/good_event_random_subtracted, 2) );




	/// Find events in the tail (side cut region)

	Int_t count = 0;

	
	Float_t top;
	Float_t down;

// 	for (int j = 0; j < hist_target->GetYaxis()->GetNbins(); j++) {
// 
// 		for (int i = 0; i < hist_target->GetXaxis()->GetNbins(); i++) {
// 	 		
// // 			if ( 
// //  					myhist_2D_cen->GetYaxis()->GetBinCenter(j)  
// //  					 myhist_2D_cen->GetXaxis()->GetBinCenter(i)
// // 
// // 				) {
// 
// // 			if ( ){ 
// 
// 				top  = Calculate_miss_x (x_mean+1, global_centered_hsbeta_cut, hist_target->GetYaxis()->GetBinCenter(j));
// 				down = Calculate_miss_x (x_mean-1, global_centered_hsbeta_cut, hist_target->GetYaxis()->GetBinCenter(j));
// 
// 	
// 				cout << "top and down" << top << "     " << down << endl;			
// 
// 
// //				}
// 
// 			}
// 
//     	}

//   	}



	
// 	cout << good_events_beta_cut  << "   " << good_events_total    << "   " 
// 		 << random_early_beta_cut << "   " << random_late_beta_cut << "   " 
// 		 << random_early_total    << "   " << random_late_total    << "   " 
// 		 << total 				  << "   " << beta_eff 			   << endl; 
// // 

	/*--------------------------------------------------*/
	// Int_t total = hist_target->GetEntries();
	
// 	float good_event_random_delta_cut_subtracted = (float) good_events_beta_cut - 
// 		((float)random_early_beta_cut + (float)random_late_beta_cut)/2.;
// 
// 	float good_event_random_subtracted = (float) good_events_total - 
// 		((float)random_early_total + (float)random_late_total)/2.;
// 
// 	float beta_eff = good_event_random_delta_cut_subtracted / good_event_random_subtracted;
// 
// 	float beta_eff_error = beta_eff * sqrt( pow(sqrt(good_event_random_subtracted)/good_event_random_subtracted, 2)  + 
// 	pow(sqrt(good_event_random_subtracted)/good_event_random_subtracted, 2) );

	printf ("%5f    %5f \n", beta_eff, beta_eff_error);

	return 0;

}




/*--------------------------------------------------*/
/// Calculate the efficiency

Int_t Calculate_efficiency(TH2F* hist_target, Float_t hsbeta_limit) {

	float x_mean = hist_target->GetMean(1);
	float y_mean = hist_target->GetMean(2);

	float y_max  = hist_target->GetYaxis()->GetXmax();
	float y_min  = hist_target->GetYaxis()->GetXmin();


//	float y_min  = hbeta_cut;

//	cout <<  y_max << "    " << y_min << endl;


	/*--------------------------------------------------*/
	Int_t good_events_beta_cut = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean-1), 
			hist_target->GetXaxis()->FindBin(x_mean + 1), 
			hist_target->GetYaxis()->FindBin(hsbeta_limit), 
			hist_target->GetYaxis()->FindBin(y_max)); 

	Int_t random_late_beta_cut = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean + 3), 
			hist_target->GetXaxis()->FindBin(x_mean + 9), 
			hist_target->GetYaxis()->FindBin(hsbeta_limit), 
			hist_target->GetYaxis()->FindBin(y_max));

	Int_t random_early_beta_cut = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean -9), 
			hist_target->GetXaxis()->FindBin(x_mean - 3 ), 
			hist_target->GetYaxis()->FindBin(hsbeta_limit), 
			hist_target->GetYaxis()->FindBin(y_max)); 



// 	/*--------------------------------------------------*/
// 	Int_t good_events_total = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean - 1), 
// 			hist_target->GetXaxis()->FindBin(x_mean + 1), 
// 			hist_target->GetYaxis()->FindBin(y_min), 
// 			hist_target->GetYaxis()->FindBin(y_max)); 
// 
// 	Int_t random_late_total = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean +3), 
// 			hist_target->GetXaxis()->FindBin(x_mean + 9), 
// 			hist_target->GetYaxis()->FindBin(y_min), 
// 			hist_target->GetYaxis()->FindBin(y_max));
// 
// 	Int_t random_early_total = hist_target->Integral(hist_target->GetXaxis()->FindBin(x_mean - 9), 
// 			hist_target->GetXaxis()->FindBin(x_mean - 3), 
// 			hist_target->GetYaxis()->FindBin(y_min), 
// 			hist_target->GetYaxis()->FindBin(y_max)); 
// 
// 
// 
// 
// 	/*--------------------------------------------------*/
// 	// Int_t total = hist_target->GetEntries();
// 	
// 	float good_event_random_delta_cut_subtracted = (float) good_events_beta_cut - 
// 		((float)random_early_beta_cut + (float)random_late_beta_cut)/2.;
// 
// 	float good_event_random_subtracted = (float) good_events_total - 
// 		((float)random_early_total + (float)random_late_total)/2.;
// 
// 	float beta_eff = good_event_random_delta_cut_subtracted / good_event_random_subtracted;
// 
// 	float beta_eff_error = beta_eff * sqrt( pow(sqrt(good_event_random_subtracted)/good_event_random_subtracted, 2)  + 
// 	pow(sqrt(good_event_random_subtracted)/good_event_random_subtracted, 2) );

	
// 	cout << good_events_beta_cut  << "   " << good_events_total    << "   " 
// 		 << random_early_beta_cut << "   " << random_late_beta_cut << "   " 
// 		 << random_early_total    << "   " << random_late_total    << "   " 
// 		 << total 				  << "   " << beta_eff 			   << endl; 
// // 


	/// Find events in the tail (side cut region)

	Int_t tail_count_good  = 0;
	Int_t tail_count_early = 0;
	Int_t tail_count_late  = 0;

	Int_t blank  = 0;
	
	Float_t top_good;
	Float_t down_good;

	Float_t top_early;
	Float_t down_early;

	Float_t top_late;
	Float_t down_late;



	for (int j = 0; j < hist_target->GetYaxis()->GetNbins(); j++) {
		for (int i = 0; i < hist_target->GetXaxis()->GetNbins(); i++) {
				top_good  = Calculate_miss_x (x_mean+1, global_centered_hsbeta_cut, hist_target->GetYaxis()->GetBinCenter(j));
				down_good = Calculate_miss_x (x_mean-1, global_centered_hsbeta_cut, hist_target->GetYaxis()->GetBinCenter(j));

				top_early  = Calculate_miss_x (x_mean-3, global_centered_hsbeta_cut, hist_target->GetYaxis()->GetBinCenter(j));
				down_early = Calculate_miss_x (x_mean-9, global_centered_hsbeta_cut, hist_target->GetYaxis()->GetBinCenter(j));

				top_late  = Calculate_miss_x (x_mean+9, global_centered_hsbeta_cut, hist_target->GetYaxis()->GetBinCenter(j));
				down_late = Calculate_miss_x (x_mean+3, global_centered_hsbeta_cut, hist_target->GetYaxis()->GetBinCenter(j));


				if ( hist_target->GetYaxis()->GetBinCenter(j) > -0.5 && 
					 hist_target->GetYaxis()->GetBinCenter(j) < global_centered_hsbeta_cut) {

					if ( hist_target->GetXaxis()->GetBinCenter(i) < top_good  &&  
						 hist_target->GetXaxis()->GetBinCenter(i) > down_good  ) {

//   						 cout << "top and down check "  <<  hist_target->GetYaxis()->GetBinCenter(j) <<  "  "
// 							 << hist_target->GetXaxis()->GetBinCenter(i) << "  " << top_good << "  " << down_good
// 							 << "  "  << hist_target->GetBinContent(i, j) << endl;			
// 
						 tail_count_good = tail_count_good + hist_target->GetBinContent(i, j);

					} else if (hist_target->GetXaxis()->GetBinCenter(i) < top_early  &&  
						 hist_target->GetXaxis()->GetBinCenter(i) > down_early  ) {

//   						 cout << "top and down check "  <<  hist_target->GetYaxis()->GetBinCenter(j) <<  "  "
// 							 << hist_target->GetXaxis()->GetBinCenter(i) << "  " << top_early << "  " << down_early
// 							 << "  "  << hist_target->GetBinContent(i, j) << endl;			
// 
						 tail_count_early = tail_count_early + hist_target->GetBinContent(i, j);

					} else if (hist_target->GetXaxis()->GetBinCenter(i) < top_late  &&  
						 hist_target->GetXaxis()->GetBinCenter(i) > down_late  ) {

//   						 cout << "top and down check "  <<  hist_target->GetYaxis()->GetBinCenter(j) <<  "  "
// 							 << hist_target->GetXaxis()->GetBinCenter(i) << "  " << top_late << "  " << down_late
// 							 << "  "  << hist_target->GetBinContent(i, j) << endl;			

						 tail_count_late = tail_count_late+ hist_target->GetBinContent(i, j);

					} else {
						blank = blank + hist_target->GetBinContent(i, j);		
					}

			}
		}
    }






// 	cout << "sdsssssssss: " << tail_count_good <<  "   " << tail_count_early << "   "
// 		 << tail_count_late << endl;


//	exit(0);



	/*--------------------------------------------------*/
	Int_t good_events_total  = good_events_beta_cut + tail_count_good;

	Int_t random_early_total = random_late_beta_cut + tail_count_late;

	Int_t random_late_total  = random_early_beta_cut + tail_count_early;


	float good_event_random_delta_cut_subtracted = (float) good_events_beta_cut - 
		((float)random_early_beta_cut + (float)random_late_beta_cut)/6.;

	float good_event_random_subtracted = (float) good_events_total - 
		((float)random_early_total + (float)random_late_total)/6.;

 	float beta_eff = good_event_random_delta_cut_subtracted / good_event_random_subtracted;
 
 	float beta_eff_error = beta_eff * sqrt( pow(sqrt(good_event_random_subtracted)/good_event_random_subtracted, 2)  + 
 	pow(sqrt(good_event_random_subtracted)/good_event_random_subtracted, 2) );

	
	cout  << " %%%  " << good_events_total + random_early_total + random_late_total << "   "<< good_events_total << "   "<<  random_late_total  << "   "<< random_early_total << "   " << blank << endl;


	printf ("%5f    %5f \n", beta_eff, beta_eff_error);

	return 0;

}





/*--------------------------------------------------*/
/// Plot boxes


void DrawGoodBox(float x_center, float y_center) {

//	cout << x_center - 1 << endl; 

	TLine *line1 = new TLine();
 	line1->SetLineWidth(1.4);
 	line1->SetLineColor(2);
	line1->DrawLine(x_center-1, y_center, x_center-1, y_center+1.15);
	line1->DrawLine(x_center+1, y_center, x_center+1, y_center+1.15);
	line1->DrawLine(x_center-1, y_center, x_center+1, y_center);

	line1->DrawLine(x_center-1, y_center, Calculate_tail_x (x_center-1, y_center ), -0.6);
	line1->DrawLine(x_center+1, y_center, Calculate_tail_x (x_center+1, y_center ), -0.6);


}



void DrawRandomBox(float x_center, float y_center) {

	TLine *line1 = new TLine();
 	line1->SetLineStyle(2);
 	line1->SetLineWidth(1.4);
	line1->DrawLine(x_center+3, y_center, x_center+3, y_center+1.15);
	line1->DrawLine(x_center+9, y_center, x_center+9, y_center+0.49);
	line1->DrawLine(x_center+3, y_center, x_center+9, y_center);

	line1->DrawLine(x_center+3, y_center, Calculate_tail_x (x_center+3, y_center ), -0.6);
	line1->DrawLine(x_center+9, y_center, Calculate_tail_x (x_center+9, y_center ), -0.6);


	TLine *line2 = new TLine();
 	line2->SetLineStyle(2);
 	line2->SetLineWidth(1.4);
	line2->DrawLine(x_center-9, y_center, x_center-9, y_center+1.15);
	line2->DrawLine(x_center-3, y_center, x_center-3, y_center+1.15);
	line2->DrawLine(x_center-9, y_center, x_center-3, y_center);

	line2->DrawLine(x_center-3, y_center, Calculate_tail_x (x_center-3, y_center), -0.6);
	line2->DrawLine(x_center-9, y_center, Calculate_tail_x (x_center-9, y_center), -0.6);

}




void Histogram_output(TH2F *hist) {


	TString run_num = hist->GetName();

	run_num.ReplaceAll("_centered", "");

	// cout << atoi(run_num) << endl;


	//hist->Write();




	bool exist1 = Check_exists(Q2_2p4_EB_4p2_vec, atoi(run_num));
	bool exist2 = Check_exists(Q2_6p5_EB_5p2_vec, atoi(run_num));
	bool exist3 = Check_exists(Q2_4p4_EB_3p7_vec, atoi(run_num));
	bool exist4 = Check_exists(Q2_5p4_EB_4p7_vec, atoi(run_num));

//	cout << exist1 << exist2 << exist3 << exist4 << endl;


	if (exist1) {

		plot_dir_Q2_2p4_EB_4p2->cd();
		Q2_2p4_EB_4p2_combine->Add(hist);

	} else if (exist2) {

		plot_dir_Q2_6p5_EB_5p2->cd();
		Q2_6p5_EB_5p2_combine->Add(hist);

	} else if (exist3) {

		plot_dir_Q2_4p4_EB_3p7->cd();
		Q2_4p4_EB_3p7_combine->Add(hist);

	} else if (exist4) {

		plot_dir_Q2_5p4_EB_4p7->cd();
		Q2_5p4_EB_4p7_combine->Add(hist);

	} 

	hist->Write();

}


bool Check_exists(vector<int> vectora, int number) {

	bool exists = find(vectora.begin(), vectora.end(), number)  != vectora.end();

	return exists;

}


void Write_overall() {

	file_out->cd();
	Q2_2p4_EB_4p2_combine->Write();
	Q2_6p5_EB_5p2_combine->Write();
	Q2_4p4_EB_3p7_combine->Write();
	Q2_5p4_EB_4p7_combine->Write();

	total_hbeta_cointime->Add(Q2_2p4_EB_4p2_combine);
	total_hbeta_cointime->Add(Q2_6p5_EB_5p2_combine);
	total_hbeta_cointime->Add(Q2_4p4_EB_3p7_combine);
	total_hbeta_cointime->Add(Q2_5p4_EB_4p7_combine);

	total_hbeta_cointime->Write();

}

Float_t Calculate_tail_x (Float_t xx, Float_t yy) {
	
	Float_t xx_tail;

	xx_tail = xx - (yy - (-0.6)) * gradi;

//	cout << "Calculation Check: "  << xx_tail << "    " << xx << "    " << float (yy - (-0.6)) * gradi << endl;

	return xx_tail;


}






Float_t Calculate_miss_x (Float_t x_boundary, Float_t y_boudary, Float_t y_pos) {


	Float_t x_pos;

	x_pos = x_boundary - (y_boudary - ( y_pos )) * gradi;

	return x_pos;

}



// int Calculate_tail_y (Float_t xx, Float_t yy) {
// 
// 	Float_t yy_tail;
// 
// 	
// 	return yy_tail;
// 
// }