#include <iostream>

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

//#include <array>

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

#include <sstream>

#include "determine_efficiency.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}; 
 

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


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();
// 
// 	Read_Offset();
// 
// 
// // 	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();
// 
// 
// //	cout << Get_Offset(47048) << endl;
// 
// 
// 
// //	loop_run(47048);


	Init();
 	Read_Offset();

	fstream run_file;
	run_file.open("../file_list_test.lst", std::fstream::in);

    Int_t a;

	vector<int> run_vec;

    while (run_file >> a) {
        printf("%i \n", a);
		run_vec.push_back(a);
    }


	

	Coin_time_loop(run_vec);



	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");

	beta_coin_dir= file_out->mkdir("beta_coin");
	aerogel_dir = file_out->mkdir("aerogel");
	cer_dir = file_out->mkdir("cer");
	cal_dir = file_out->mkdir("cal");





	/*--------------------------------------------------*/
	/// 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);



	/*--------------------------------------------------*/
	/// Important Parameter definition for the 
	///
	///	beta_cut: this beta cut is allied if the cointime is assumed to be independent wrt beta 
	///
	/// global_centered_hsbeta_cut: this parameter gives the hsbeta limit of when the cointime becomes 
	///								beta independent. For events have hsbeta less than this value, 
	///                             the coincidence time has to be calculated.  
	///	
	/// gradi: since the cut on the proton coincidence time is beta dependent, 
	
	beta_cut                   = 0.1;
	global_centered_hsbeta_cut = -0.15;
	gradi                      = 1./9.;	




	offset_file_in.open("offset.dat", ios::in);

	plot_put_dir = "plot/";

	efficiency_file = fopen( "detector_efficiency.dat", "w" );

	//efficiency_file << "!Run   gbeta   err   gcer   err   gaero   err" << endl;

	fprintf (efficiency_file, "!Run   beta   err    cer    err    aero   err    cal    err\n");	

}

/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// The loop that actually ceing used. This loop takes into account that cointime is dependent 
/// of beta.

void loop_run(int run_num) {

	Event *event = new Event();

	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");

	//tree->SetBranchAddress("event", &event);
 
	TH2F* test_beta_vs_cointime  = new TH2F(plot_title, "", 441, -11, 11, 501, -1.5, 1);

	int events_num = tree->GetEntries();
	cout << events_num << endl;

	//TTree *new_tree = tree->CopyTree("scer_npe > 5 && ssshsum > 0.8 && ssshsum < 1.2");


	TFile* dummy_file = new TFile("dummy_file", "RECREATE");
	
	dummy_file->cd();

	TTree * new_tree = tree->CopyTree("abs(haero_su) <= 2.5 && cointime > -11 && cointime < 4 && hcer_npe < 2.0");


	float hcer_npe;
	float hsbeta;
	float cointime;
	float haero_su;

	new_tree->SetBranchAddress("hcer_npe", &hcer_npe);
	new_tree->SetBranchAddress("hsbeta",   &hsbeta);
	new_tree->SetBranchAddress("cointime", &cointime);
	new_tree->SetBranchAddress("haero_su", &haero_su);

	//tree->SetBranchAddress("scer_npe", &temp_num);

	float run_offset = Get_Offset(run_num);

	for (int i=0; i < new_tree->GetEntries(); i++) {

		new_tree->GetEntry(i);
		// tree->GetEntry(i);

		//cout << temp_num  << "    " << hsbeta   <<  "     " << cointime << endl;	

		Float_t cointime_corrected = cointime - run_offset;
		Float_t hsbeta_corrected   = hsbeta - 1;




// 		Float_t cointime_upper_boudary = Calculate_miss_x( 1, global_centered_hsbeta_cut, hsbeta_corrected);
// 		Float_t cointime_lower_boudary = Calculate_miss_x(-1, global_centered_hsbeta_cut, hsbeta_corrected);
// 
// 
// 
// 		if( (abs(cointime_corrected) < 1 && hsbeta_corrected > -0.15)) {
// 	
// 
// // 		if( (abs(cointime_corrected) < 1 && hsbeta_corrected > -0.15) || 
// // 			(hsbeta_corrected <= -0.15 && hsbeta_corrected > -0.6 && 
// // 			 cointime_corrected < cointime_upper_boudary && cointime_corrected> cointime_lower_boudary) ) {
// // 		
// //			cout << "boundary sati: " << hsbeta_corrected << "    " << cointime_corrected << endl;
// //			test_beta_vs_cointime->Fill(cointime_corrected, hsbeta_corrected);
// 
// 
// 		}
// 

//		if (hcer_npe <= 0.5){


			if (abs(haero_su) < 50.0 && hcer_npe < 2.0 )
 			  test_beta_vs_cointime->Fill(cointime_corrected, hsbeta_corrected);

//		}

	}

	TCanvas *c1 = new TCanvas();

	test_beta_vs_cointime->Draw("hist");

	DrawGoodBox(0, global_centered_hsbeta_cut);
	DrawRandomBox(0, global_centered_hsbeta_cut);
	
	c1->Print(plot_put_dir + "beta_cointime_"+ plot_title + ".eps");

	
	beta_coin_dir->cd();
	c1->Write(plot_title);
	

 	Calculate_efficiency(test_beta_vs_cointime, global_centered_hsbeta_cut);

	// cout <<  new_tree->GetEntries() << endl;





	/*--------------------------------------------------*/
	/// Define Cer Cut 
	TString cer_cut_str;
	cer_cut_str.Form("cointime > %f && cointime < %f &&  abs(haero_su) < 2.5", run_offset-1, run_offset+1.);
	TCut cer_cut = (TCut) cer_cut_str;


	/*--------------------------------------------------*/
	/// Define Aero Cut 
	TString aero_cut_str;
	aero_cut_str.Form("cointime > %f && cointime < %f && hcer_npe < 2.0 && haero_su > -50", run_offset-1, run_offset+1.);
	TCut aero_cut = (TCut) aero_cut_str;


	/*--------------------------------------------------*/
	/// Define Cal Cut 
	TString cal_cut_str;
	cal_cut_str.Form("cointime > %f && cointime < %f && hcer_npe < 2.0 && abs(haero_su) < 50", run_offset-1, run_offset+1.);
	TCut cal_cut = (TCut) cal_cut_str;


	Cer_plot(new_tree, "cer_" + plot_title, cer_cut, 2.0, 0, 4, 40);

	Aero_plot(new_tree, "aero_" + plot_title, aero_cut, 50, -60, 140, 200);

	Cal_plot(new_tree, "cal_" + plot_title, cal_cut, 0.5, 0, 0.6, 60);

	fprintf (efficiency_file, "%i  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f  %.3f\n", run_num, gbeta_eff, gbeta_eff_err, gcer_eff, gcer_eff_err, gaero_eff, gaero_eff_err, cal_eff, cal_eff_err);	



	delete event;
	delete c1;
	delete test_beta_vs_cointime;
	delete new_tree;
	delete file;
	delete dummy_file;


}




/*--------------------------------------------------*/
// Read Offset



void Read_Offset() {

	string test_str;

	int run_num_tmp;
	float offset_tmp;
	float offset_err_tmp;

	if (offset_file_in.is_open()) {

		//while(!offset_file_in.eof()){
		//}

		while(getline(offset_file_in, test_str)){

			istringstream ss(test_str);
			ss >> run_num_tmp >> offset_tmp >> offset_err_tmp;
			//cout <<  run_num_tmp << "  "  << offset_tmp  << endl;
			offset_map[run_num_tmp] = offset_tmp;

		}

	}

}


float Get_Offset(int run_num_tmp) {

	return offset_map.find(run_num_tmp)->second;

}







/*--------------------------------------------------*/
// 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; 
	

	delete myhist;
	delete c5;

}

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

void Cer_plot(TTree*target_tree, 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(plot_title, "",  bin_num, bin_min, bin_max);

 	TCanvas *c5 = new TCanvas();
 	target_tree->Draw("hcer_npe>>" + plot_title, 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; 

	c5->Print(plot_put_dir + plot_title+".eps");

	cer_dir->cd();
	c5->Write(plot_title);


	float eff_error =  eff * sqrt( pow(sqrt(lower)/lower, 2) + pow(sqrt(upper+lower)/(upper+lower), 2) );

	gcer_eff = eff;
	gcer_eff_err = eff_error;	


	delete myhist;
	delete c5;

}

/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// 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; 


	delete myhist;
	delete c5;

}


/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// 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; 
	
}

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


void Aero_plot(TTree*target_tree, 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(plot_title, "", bin_num, bin_min, bin_max);

 	TCanvas *c5 = new TCanvas();
 	target_tree->Draw("haero_su>>" + plot_title, 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; 
	
	c5->Print(plot_put_dir + plot_title+".eps");


	float eff_error =  eff * sqrt( pow(sqrt(lower)/lower, 2) + pow(sqrt(upper+lower)/(upper+lower), 2) );

	gaero_eff = eff;
	gaero_eff_err = eff_error;

	aerogel_dir->cd();
	c5->Write(plot_title);


	delete myhist;
	delete c5;


}



/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// Computing the HMS Calorimeter Efficiency


void Cal_plot(TTree*target_tree, 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(plot_title, "", bin_num, bin_min, bin_max);

 	TCanvas *c5 = new TCanvas();
 	target_tree->Draw("hsshtrk>>" + plot_title, 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; 
	
	c5->Print(plot_put_dir + plot_title+".eps");


	float eff_error =  eff * sqrt( pow(sqrt(lower)/lower, 2) + pow(sqrt(upper+lower)/(upper+lower), 2) );

	cal_eff = eff;
	cal_eff_err = eff_error;

	cal_dir->cd();
	c5->Write(plot_title);


	delete myhist;
	delete c5;


}



/*--------------------------------------------------*/
/*--------------------------------------------------*/
/// 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 );		
		loop_run(grp[i]);

	}

}



/*--------------------------------------------------*/
/// 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 abs(haero_su) <= 2.5  && cointime > -10 && cointime < -4", "goff");

	tree->Draw("hsbeta:cointime>>" + plot_title, "abs(haero_su) <= 2.5 && cointime < 3.9 && cointime > -17.9");
	tree->Draw("hsbeta:cointime>>get_center", "abs(haero_su) <= 2.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;
	delete file;
	delete myhist_2D;
	delete myhist_2D_cen;
	delete tree;

}


/*--------------------------------------------------*/
/// 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 -5), 
			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-5, 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.6 && 
					 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)/4.;

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

 	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);



	gbeta_eff = beta_eff;

	gbeta_eff_err = 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-5, y_center, x_center-5, y_center+1.15);
	line2->DrawLine(x_center-3, y_center, x_center-3, y_center+1.15);
	line2->DrawLine(x_center-5, 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-5, y_center, Calculate_tail_x (x_center-5, 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;
// 
// }
//
//
//
//