#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "../src/TNeEvent.h"
#include "../../AculUtils/TELoss/TELoss.h"

using std::cout;
using std::endl;

void fillTree(Int_t nofile = 0, const Int_t noevents = 0) {

	TString inFile;
	inFile.Form("~/data/exp1804/h5_14_%02d.root", nofile);
//	inFile.Form("~/data/exp1804/calib/si_20_03.root");

	TString outFile;
	outFile.Form("~/data/exp1804/h5_14_%02d_calib.root", nofile);
//	outFile.Form("~/data/exp1804/calib/si_20_03_calib.root");

	cout << "Input file: " << inFile << endl;
	cout << "Output file: " << outFile << endl;

	TFile *fr = new TFile(inFile);
	TTree *tr = (TTree*)fr->Get("AnalysisxTree");

	TNeEvent *revent = new TNeEvent();

	tr->SetBranchAddress("NeEvent.", &revent);

	TFile *fw = new TFile(outFile, "RECREATE");
	TTree *tw = new TTree("cal", "Calibrated information");

	Float_t SQ20E[16];
	Float_t SQ20Ecorr[16];
	Float_t SQ20Esum;
	Float_t SQLXE[32];
	Float_t SQLXEsum;
	Float_t SQLYE[16];
	Float_t SQLYEsum;

	Int_t SQLXmult;
	Int_t SQLYmult;
	Int_t SQ20EcorrMult;

	Float_t SQRXE[32];
	Float_t SQRXEsum;
	Int_t SQRXmult;
	Float_t TOF, dEbeam;
	Float_t x1p, y1p, x2p, y2p, xt, yt;

	tw->Branch("SQ20E",SQ20E,"SQ20E[16]/F");
	tw->Branch("SQ20Ecorr",SQ20Ecorr,"SQ20Ecorr[16]/F");
	tw->Branch("SQ20Esum",&SQ20Esum,"SQ20Esum/F");
	tw->Branch("SQLXE",SQLXE,"SQLXE[32]/F");
	tw->Branch("SQLXEsum",&SQLXEsum,"SQLXE/F");
	tw->Branch("SQLYE",SQLYE,"SQLYE[16]/F");
	tw->Branch("SQLYEsum",&SQLYEsum,"SQLYEsum/F");

	tw->Branch("SQLXmult",&SQLXmult,"SQLXmult/I");
	tw->Branch("SQLYmult",&SQLYmult,"SQLYmult/I");
	tw->Branch("SQ20EcorrMult",&SQ20EcorrMult,"SQ20EcorrMult/I");

	tw->Branch("SQRXE",SQRXE,"SQRXE[32]/F");
	tw->Branch("SQRXEsum",&SQRXEsum,"SQRXE/F");
	tw->Branch("SQRXmult",&SQRXmult,"SQRXmult/I");

	tw->Branch("TOF", &TOF, "TOF/F");
	tw->Branch("dEbeam", &dEbeam, "dEbeam/F");

	tw->Branch("x1p", &x1p, "x1p/F");
	tw->Branch("y1p", &y1p, "y1p/F");
	tw->Branch("x2p", &x2p, "x2p/F");
	tw->Branch("y2p", &y2p, "y2p/F");

	tw->Branch("xt", &xt, "xt/F");
	tw->Branch("yt", &yt, "yt/F");


	Long64_t nevents;
	if (noevents == 0) nevents = tr->GetEntries();
	else nevents = noevents;
	if (nevents > tr->GetEntries()) nevents = tr->GetEntries();

//	TNeDet16 *pSQX_L_EC = new TNeDet16("SQX_L_EC");
//	TNeDet16 pSQX_L_EC("../SQX_L_EC");
	TNeDet16 pSQX_L_EC("./SQX_L");
	pSQX_L_EC.ReadData();

	TNeDet16 pSQX_L_TC("./SQX_L_time");
	pSQX_L_TC.ReadData();

	TNeDet16 pSQY_L_EC("./SQY_L");
	pSQY_L_EC.ReadData();

	TNeDet16 pSQX_R_EC("../SQX_R_EC");
	pSQX_R_EC.ReadData();

	TNeDet16 pSQ20_EC("./sq20_58");
	pSQ20_EC.ReadData();

//	for (Int_t i = 0; i < 32; i++) {
//		cout << pSQX_L_EC.Energy(1, i) << endl;
//	}

	Float_t energy = 0;

	cout << nevents << " entries will be treated." << endl;

//	TFile fThickness("thicknessPreliminary.root", "OPEN");
	TFile fThickness("thicknessFinal.root", "OPEN");
	fr->cd();
	TH2F *hThickness = new TH2F(*(TH2F*)fThickness.Get("hTh"));
	hThickness->Draw("col");

//	std::cout << std::setprecision(1) << std::fixed;


	const Int_t kSQL_X_strips = 32;
	const Int_t kSQL_Y_strips = 16;
	const Int_t kSQL_20_strips = 16;

	const Double_t kSQLY_energy_thr = 1.;
	const Double_t kSQLX_energy_thr = 1.;
	const Double_t kSQL20_energy_thr = 1.2;

	const Double_t kSQ20_norm_thickness = 20.;

//	for (Int_t yi = 0; yi < kSQL_Y_strips; yi++) {
//		cout << "y bin: " << yi+1 << "\t\t";
//		for (Int_t xi = 0; xi < kSQL_X_strips; xi++) {
//			cout << hThickness->GetBinContent(xi+1, yi+1) << "\t";
//			if (xi == kSQL_X_strips-1) cout << endl;
//		}
//	}


	fw->cd();


	TELoss fAlphaSi;
	fAlphaSi.SetEL(1, 2.321); // density in g/cm3
	fAlphaSi.AddEL(14., 28.086, 1);  //Z, mass
	//	mSi.SetZP(1., 1.);		//protons
	fAlphaSi.SetZP(2., 3.);		//3He, Z, A
	fAlphaSi.SetEtab(100000, 100.);	// ?, MeV calculate ranges
	fAlphaSi.SetDeltaEtab(300);

	//////////////////////////////////
	//event processing
	//////////////////////////////////

	for (Int_t eventNo = 0; eventNo < nevents; eventNo++) {
		tr->GetEvent(eventNo);
//		cout << revent->SQX_L[0] << endl;
//		SQLXE[0] = revent->SQX_L[0];

		SQ20Esum = 0.;
		SQLXEsum = 0.;
		SQLYEsum = 0.;
		SQRXEsum = 0.;
		SQRXmult = 0;

		SQLXmult = 0;
		SQLYmult = 0;
		SQ20EcorrMult = 0;


		for (Int_t i = 0; i < 32; i++) {
//			cout << pSQX_L_EC.Energy(1, i) << endl;
			SQLXE[i] = 0;
			SQRXE[i] = 0;

			energy = pSQX_L_EC.Energy(revent->SQX_L[i], i);

//			cout << energy << "\t" << pSQX_L_EC.a[i][0] << "\t" << pSQX_L_EC.a[i][1] << "\t" << pSQX_L_EC.a[i][2] << endl;
//			cout << energy << endl;

//			if (i == 5 && i == 6) continue;

			if (energy > kSQLX_energy_thr) {
				SQLXE[i] = energy;
				SQLXEsum += SQLXE[i];
				SQLXmult++;
			}

			energy = pSQX_R_EC.Energy(revent->SQX_R[i], i);

			if (energy > 1.0) {
				SQRXE[i] = energy;
				SQRXEsum += SQRXE[i];
				SQRXmult++;
			}
		}

		for (Int_t i = 0; i < 16; i++) {
			//			cout << pSQX_L_EC.Energy(1, i) << endl;

			SQLYE[i] = 0.;
			SQ20E[i] = 0.;
			SQ20Ecorr[i] = 0.;
			//				SQRXE[i] = 0;

			energy = pSQ20_EC.Energy(revent->SQ20[i], i);

//			if (energy > 1.) {
				SQ20E[i] = energy;
				SQ20Esum += SQ20E[i];
//				cout << energy << endl;
//				cout << i << "\t" << energy << "\t" << pSQ20_EC.a[i][0] << "\t" << pSQ20_EC.a[i][1] << "\t" << pSQ20_EC.a[i][2] << endl;
//			}


			if (i==0) continue;

			energy = pSQY_L_EC.Energy(revent->SQY_L[i], i);

			//				cout << energy << "\t" << pSQY_L_EC.a[i][0] << "\t" << pSQY_L_EC.a[i][1] << "\t" << pSQY_L_EC.a[i][2] << endl;
			//			cout << energy << endl;

			if (energy > kSQLY_energy_thr) {
				SQLYE[i] = energy;
				SQLYEsum += SQLYE[i];
				SQLYmult++;
			}

			//				energy = pSQX_R_EC.Energy(revent->SQX_R[i], i);
			//
			//				if (energy > 1.0) {
			//					SQRXE[i] = energy;
			//					SQRXEsum += SQRXE[i];
			//					SQRXmult++;
			//				}
		}//for 16

//		cout << endl;

		///////////////////////////////////////////////
		//correction for thickness inhomogeneity
		///////////////////////////////////////////////

		Double_t currThickness;

		if (SQLXmult == 1 && SQLYmult == 1) {

			for (Int_t yi = 0; yi < kSQL_Y_strips; yi++) {
				//			cout << "y bin: " << yi+1 << "\t\t";
				if (SQLYE[yi]>kSQLY_energy_thr) {

					for (Int_t xi = 0; xi < kSQL_X_strips; xi++) {

						if (SQLXE[xi] > kSQLX_energy_thr) {
							currThickness = hThickness->GetBinContent(xi+1, yi+1);
//							if (currThickness<30.) {
//								Int_t probableThinStrip = xi-2/2;
//								SQ20Ecorr[0] =
//							}

							for (Int_t x20 = 0; x20 < kSQL_20_strips; x20++) {
								if (SQ20E[x20] > kSQL20_energy_thr && currThickness<30.) {
//									SQ20Ecorr[x20] = SQ20E[x20] + 1.;

//									currThickness = kSQ20_norm_thickness + (currThickness - kSQ20_norm_thickness)*0.5;

									SQ20Ecorr[x20] = (kSQ20_norm_thickness/currThickness)*SQ20E[x20];
									SQ20EcorrMult++;

//									cout << xi << "\t" << yi << "\t" << x20
//											<< "\t" << currThickness
//											<< "\t" << SQ20E[x20]  << "\t" << SQ20Ecorr[x20] << "\t" << endl;
								}
							}//for*/
						}



						//					cout << hThickness->GetBinContent(xi+1, yi+1) << "\t";
						//					if (xi == kSQL_X_strips-1) cout << endl;
					}
				}
			}

		}


		///////////////////////////////////////////////
		//beam diagnostics
		///////////////////////////////////////////////

		dEbeam = 0.;
		TOF = 0.;
		if (revent->tF5[0]>0 && revent->tF5[1]>0 && revent->tF5[2]>0 && revent->tF5[3]>0
				&& revent->tF3[0]>0 && revent->tF3[1]>0 && revent->tF3[2]>0 && revent->tF3[3]>0
				&& revent->F5[0]>0 && revent->F5[1]>0 && revent->F5[2]>0 && revent->F5[3]>0
//				&& ( (revent->tF5[0]+revent->tF5[1]+revent->tF5[2]+revent->tF5[3]) - (revent->tF3[0]+revent->tF3[1]+revent->tF3[2]+revent->tF3[3]) )/4.*0.125+89.165<200
//				&& ( (tF5[0]+tF5[1]+tF5[2]+tF5[3]) - (tF3[0]+tF3[1]+tF3[2]+tF3[3]) )/4.*0.125+89.165>100
//				&& (F5[0]+F5[1]+F5[2]+F5[3])/4. < 2500
		) {
			dEbeam = (revent->F5[0]+revent->F5[1]+revent->F5[2]+revent->F5[3])/4.;
			TOF = ( (revent->tF5[0]+revent->tF5[1]+revent->tF5[2]+revent->tF5[3]) - (revent->tF3[0]+revent->tF3[1]+revent->tF3[2]+revent->tF3[3]) )/4.*0.125+89.165;
		}

		const Float_t l12 = 546.;
		const Float_t lt = 270.;

		x1p = -100.;
		y1p = -100.;
		x2p = -100.;
		y2p = -100.;
		xt = -100.;
		yt = -100.;

		if (revent->x1[0]<1000 && revent->y1[0]<1000 && revent->nx1==1 && revent->ny1==1) {
			x1p = revent->x1[0]*1.25-20.;
			y1p = revent->y1[0]*1.25-20.;
		}

		if (revent->x2[0]<1000 && revent->y2[0]<1000 && revent->nx2==1 && revent->ny2==1) {
			x2p = revent->x2[0]*1.25-20.;
			y2p = revent->y2[0]*1.25-20.;
		}

		if (x1p > -50. && y1p > -50. && x2p > -50. && y2p > -50.) {
			xt = x1p - (x1p - x2p)*(l12/lt);
			yt = y1p - (y1p - y2p)*(l12/lt);
		}


//		cout << endl;

		tw->Fill();

//		cout << eventNo << "\t" << eventNo%100 << endl;
		if (eventNo%100000 == 0 && eventNo !=0) {
			cout << eventNo << " events processed." << endl;
		}
	}

	cout << nevents << " entries processed." << endl;

	fw->cd();
	tw->Write();
	fw->Close();

}


void fillChain(const Int_t from = 0, const Int_t to = 9, const Int_t noevents = 0) {
	for (Int_t i = from; i <= to; i++) {
		fillTree(i, noevents);
	}
}