#include <iostream>

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

//#include <array>

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

#include "determine_offset.h"

using namespace std;



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
//
// For Omega

int main() {


//  
//  	off->Coin_time_loop(Q2_2p4_EB_4p2_vec);
//  	off->Coin_time_loop(Q2_6p5_EB_5p2_vec);
//  	off->Coin_time_loop(Q2_4p4_EB_3p7_vec);
//  	off->Coin_time_loop(Q2_5p4_EB_4p7_vec);
 
	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);
    }


 	Offset *off = new Offset;
 	off->Init();
	

  	off->Coin_time_loop(run_vec);


// 	Write_overall();

	off->Cleanup();




	return 0;
}







/*--------------------------------------------------*/
/// Initialization
void Offset::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");

	offset_out_dir = file_out->mkdir("offset");


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

	fit_fun = new TF1("my_function", "gaus");
	fit_fun_improve = new TF1("my_function_imp", "gaus");

	offset_file = fopen( "offset.dat", "w" );

}



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

//void Coin_time_loop(Int_t* grp) {

void Offset::Coin_time_loop(vector<int> grp) {
	
 	for(int i=0; i < grp.size(); i++) {
		cout << grp[i] << "    " << endl;
 		Run_by_Run(grp[i], -8, -4.5, 60);		
// 		Run_by_Run(grp[i], -20, 20, 100);		
	}

}



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

void Offset::Run_by_Run(Int_t run_num, Float_t x_bin_min, Float_t x_bin_max, Int_t x_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();

 	TH1F *myhist = new TH1F(plot_title, plot_title, x_bin_num, x_bin_min, x_bin_max);

 	tree->Draw("cointime >> " + plot_title, "hsbeta > 0.1 && hsbeta < 1.5 && hcer_npe < 0.5 && abs(haero_su) < 50 && ssshtrk > 0.70 && scer_npe > 0.50 && hsbeta > 0.85 && hsbeta < 1.15", "goff");

	myhist->Draw("hist");

	Float_t fun_mean = myhist->GetMean();
	fit_fun->SetRange(fun_mean-0.4, fun_mean+0.45);

	// myhist->Fit("my_function", "QR");
	myhist->Fit("my_function", "R");


	fit_fun_improve->SetRange(fit_fun->GetParameter("Mean")-0.41, fit_fun->GetParameter("Mean")+0.43);

	myhist->Fit("my_function_imp", "R");

	fit_fun_improve->Draw("same");
	



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

  	c1->Print(ps_plotname);

	cout << "Bin width: " << myhist->GetBinWidth(10) << endl;

	cout << "Mean: " << myhist->GetMean() << "    "<< fit_fun_improve->GetParameter("Mean") << endl;
	cout << "Mean Error: " << myhist->GetMeanError() << "    "<< fit_fun_improve->GetParError(1) << endl;

	fprintf (offset_file, " %i  %.5f  %.5f\n", run_num, fit_fun_improve->GetParameter("Mean"), fit_fun_improve->GetParError(1));	

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


}



/*--------------------------------------------------*/
// Clean

void Offset::Cleanup() {
	
}