/*
 * BeAnalysis.cpp
 *
 *  Created on: Jul 19, 2017
 *      Author: vratik
 */

#include "BeAnalysis.h"

#include "BeWork.h"
#include "TCanvas.h"
#include "TGaxis.h"

BeAnalysis::BeAnalysis() : che(0), spectra(1) {
	// TODO Auto-generated constructor stub

	for (Int_t i = 0; i < 6; i++) {
		chs[i] = NULL;
		ti[i] = NULL;
	}

	SetCuts();
	SetChainsToDraw();

	noIntervals = 5;
	epsilonT = new Bool_t[noIntervals];
	for (Int_t i = 0; i < noIntervals; i++) {
		epsilonT[i] = 0;
	}
	epsilonT[0] = 1;
	epsilonT[2] = 1;


	cosThetaT = new Bool_t[noIntervals];
	for (Int_t i = 0; i < noIntervals; i++) {
		cosThetaT[i] = 0;
	}
	cosThetaT[0] = 1;
	cosThetaT[2] = 1;
	cosThetaT[4] = 1;

	epsilonY = new Bool_t[noIntervals];
	for (Int_t i = 0; i < noIntervals; i++) {
		epsilonY[i] = 0;
	}
	epsilonY[0] = 1;
	epsilonY[2] = 1;


	kVerbose = 1;
	kRangeProportionIn = 2.0;

}

BeAnalysis::~BeAnalysis() {
	// TODO Auto-generated destructor stub
}

void BeAnalysis::OpenSimChains() {

	const Int_t lowSimFile[6] = {0, 0, 0, 0, 0, 0};
	const Int_t upSimFile[6] = {14, 14, 14, 14, 14, 14};

	TString chainAl0name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_0_50-85_";
	TString inputTreeAl0name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_0_50-85_000.root";
	//isotropic, 0 degrees
	TString chainNoAl0name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_0_50-85_";
	TString inputTreeNoAl0name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_0_50-85_000.root";
	//aligned, 180 degrees
	TString chainAl180name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_180_50-85_";
	TString inputTreeAl180name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_180_50-85_000.root";
	//isotropic, 180 degrees
	TString chainNoAl180name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_180_50-85_";
	TString inputTreeNoAl180name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_180_50-85_000.root";
	//isotropic, 90 degrees
	TString chainAl90name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_90_50-85_";
	TString inputTreeAl90name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_al_90_50-85_000.root";
	//isotropic, 90 degrees
	TString chainNoAl90name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_90_50-85_";
	TString inputTreeNoAl90name = "../../../be/rootdata/correlations/v5_6075/Sim_mix_br_noal_90_50-85_000.root";

//	TChain *chsAl0 = BeWork::OpenChain(chainAl0name.Data(),
	chs[0] = BeWork::OpenChain(chainAl0name.Data(),
			lowSimFile[0], upSimFile[0], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[0]->GetEntries(), chs[0]->GetName());
	ti[0] = BeWork::OpenTree(inputTreeAl0name.Data(), "sbeam", 2);
	//	tiAl0->SetLineWidth(2);

//	TChain *chsNoAl0 = BeWork::OpenChain(chainNoAl0name.Data(),
	chs[1] = BeWork::OpenChain(chainNoAl0name.Data(),
			lowSimFile[1], upSimFile[1], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[1]->GetEntries(), chs[1]->GetName());
	ti[1] = BeWork::OpenTree(inputTreeNoAl0name.Data(), "sbeam", 2);
	//	tiNoAl0->SetLineWidth(2);
	//
//	TChain *chsAl180 = BeWork::OpenChain(chainAl180name.Data(),
	chs[2] = BeWork::OpenChain(chainAl180name.Data(),
			lowSimFile[2], upSimFile[2], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[2]->GetEntries(), chs[2]->GetName());
	ti[2] = BeWork::OpenTree(inputTreeAl180name.Data(), "sbeam", 2);
	//	tiAl180->SetLineWidth(2);

//	TChain *chsNoAl180 = BeWork::OpenChain(chainNoAl180name.Data(),
	chs[3] = BeWork::OpenChain(chainNoAl180name.Data(),
			lowSimFile[3], upSimFile[3], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[3]->GetEntries(), chs[3]->GetName());
	ti[3] = BeWork::OpenTree(inputTreeNoAl180name.Data(), "sbeam", 2);
	//	tiNoAl180->SetLineWidth(2);

//	TChain *chsAl90 = BeWork::OpenChain(chainAl90name.Data(),
	chs[4] = BeWork::OpenChain(chainAl90name.Data(),
			lowSimFile[4], upSimFile[4], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[4]->GetEntries(), chs[4]->GetName());
	ti[4] = BeWork::OpenTree(inputTreeAl90name.Data(), "sbeam", 2);
	//	tiAl90->SetLineWidth(2);

//	TChain *chsNoAl90 = BeWork::OpenChain(chainNoAl90name.Data(),
	chs[5] = BeWork::OpenChain(chainNoAl90name.Data(),
			lowSimFile[5], upSimFile[5], "simbe", 4, "sbeam");
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[5]->GetEntries(), chs[5]->GetName());
	ti[5] = BeWork::OpenTree(inputTreeNoAl90name.Data(), "sbeam", 2);
	//	tiNoAl90->SetLineWidth(2);

}

void BeAnalysis::OpenExpChain() {

	Int_t lowExpFile = 0;
	Int_t upExpFile = 40;

	//experimental chain
	che = BeWork::OpenChain("../../../be/rootdata/correlations/Be.",
			lowExpFile, upExpFile, "beonly");										//original file
	Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing experimental data",
			che->GetEntries(), che->GetName());

}

void BeAnalysis::SetCuts() {
	cBe20 = "fBeIM>0 && fBeIM<20";
	cBe3 = "fBeIM>0 && fBeIM<3";
	cBeWork = "fBeIM>0 && fBeIM<10";

	cBe0_14 = "fBeIM>0 && fBeIM<1.4";
	//	TCut cBe0_14 = "fBeIM>0 && fBeIM<1.6";
	//	TCut cBe0_14 = "fBeIM>0 && fBeIM<1.20";
	cBe14_19 = "fBeIM>1.4 && fBeIM<1.9";
	cBe19_25 = "fBeIM>1.9 && fBeIM<2.5";
	cBe25_31 = "fBeIM>2.5 && fBeIM<3.1";
	cBe31_37 = "fBeIM>3.1 && fBeIM<3.7";
	//array:
	cBeE[0] = cBe0_14;
	cBeE[1] = cBe14_19;
	cBeE[2] = cBe19_25;
	cBeE[3] = cBe25_31;
	cBeE[4] = cBe31_37;

	cEpsilonT = "fTpp/fBeIM<0.2";
	//	TCut cEpsilonT = "fTpp/fBeIM<1.";
	//	TCut cEpsilonY = "fTap/fBeIM<0.5";
	cEpsilonY = "fTap/fBeIM>0.7";

	//raw files
	//energy cuts
	crBe20 = "f6BeIM>0 && f6BeIM<20";
	crBe3 = "f6BeIM>0 && f6BeIM<3";
	crBeWork = "f6BeIM>0 && f6BeIM<10";

	crBe0_14 = "f6BeIM>0 && f6BeIM<1.4";
	//	TCut crBe0_14 = "f6BeIM>0 && f6BeIM<1.6";
	//	TCut crBe0_14 = "f6BeIM>0.5 && f6BeIM<1.20";
	crBe14_19 = "f6BeIM>1.4 && f6BeIM<1.9";
	crBe19_25 = "f6BeIM>1.9 && f6BeIM<2.5";
	crBe25_31 = "f6BeIM>2.5 && f6BeIM<3.1";
	crBe31_37 = "f6BeIM>3.1 && f6BeIM<3.7";
	//array:
	crBeE[0] = crBe0_14;
	crBeE[1] = crBe14_19;
	crBeE[2] = crBe19_25;
	crBeE[3] = crBe25_31;
	crBeE[4] = crBe31_37;
	//angular cuts


	//	TCut crEpsilonT = "fTpp/f6BeIM>0.8";
	crEpsilonT = "fTpp/f6BeIM<0.2";
	//	TCut crEpsilonT = "fTpp/f6BeIM<1.";
	//	TCut crEpsilonY = "fTap/f6BeIM<0.5";
	crEpsilonY = "fTap/f6BeIM>0.7";

	//simulation input
	ciBe0_14 = "E_IM>0 && E_IM<1.4";
	//	TCut ciBe0_14 = "E_IM>0 && E_IM<1.6";
	//	TCut ciBe0_14 = "E_IM>0.5 && E_IM<1.20";
	ciBe14_19 = "E_IM>1.4 && E_IM<1.9";
	ciBe19_25 = "E_IM>1.9 && E_IM<2.5";
	ciBe25_31 = "E_IM>2.5 && E_IM<3.1";
	ciBe31_37 = "E_IM>3.1 && E_IM<3.7";
	//array:
	ciBeE[5] = ciBe0_14;
	ciBeE[5] = ciBe14_19;
	ciBeE[5] = ciBe19_25;
	ciBeE[5] = ciBe25_31;
	ciBeE[5] = ciBe31_37;

	ciEpsilon = "sTpp/E_IM<0.2";
	//	TCut ciEpsilon = "sTpp/E_IM<1.";
	//	TCut ciEpsilonY = "sTap/E_IM<0.5";
	ciEpsilonY = "sTap/E_IM>0.7";

	cQ = "TMath::Abs(fQLiP)<10";
	cProtons = "fP1Lab.fE-938.272<50 && fP2Lab.fE-938.272<50";
	//	TCut cProtons = "fP1Lab.fE-938.272<34 && fP2Lab.fE-938.272<34 && fP1Lab.fE-938.272>1 && fP2Lab.fE-938.272>1";
}

void BeAnalysis::SetCMAngularRange(Int_t minAngle, Int_t maxAngle) {
	kMinAngle = minAngle;
	kMaxAngle = maxAngle;
}

void BeAnalysis::SetChainsToDraw(Bool_t ch0, Bool_t ch1, Bool_t ch2, Bool_t ch3, Bool_t ch4,
		Bool_t ch5) {

	kChains[0] = ch0;
	kChains[1] = ch1;
	kChains[2] = ch2;
	kChains[3] = ch3;
	kChains[4] = ch4;
	kChains[5] = ch5;

}

void BeAnalysis::SetRangeProportion(Float_t rangeProportion, Bool_t autoRange) {
	kRangeProportion = rangeProportion;
	kAutoRange = autoRange;
}

void BeAnalysis::SetNoExpEvents() {

	eMaxEvents = 5000000;

	eEventsAl0 = eMaxEvents;
	eEventsNoAl0 = eMaxEvents;
	eEventsAl180 = eMaxEvents;
	eEventsNoAl180 = eMaxEvents;
	eEventsAl90 = eMaxEvents;
	eEventsNoAl90 = eMaxEvents;

	eEvents[0] = eEventsAl0;
	eEvents[1] = eEventsNoAl0;
	eEvents[2] = eEventsAl180;
	eEvents[3] = eEventsNoAl180;
	eEvents[4] = eEventsAl90;
	eEvents[5] = eEventsNoAl90;

}

void BeAnalysis::SetNoSimEvents() {

	sMaxEvents = 3000000;

	sEventsAl0 = sMaxEvents;
	sEventsNoAl0 = 2700000;
	sEventsAl180 = sMaxEvents;
	sEventsNoAl180 = 2700000;
	sEventsAl90 = 2690000;
	sEventsNoAl90 = 2710000;

	sEvents[0] = sEventsAl0;
	sEvents[1] = sEventsNoAl0;
	sEvents[2] = sEventsAl180;
	sEvents[3] = sEventsNoAl180;
	sEvents[4] = sEventsAl90;
	sEvents[5] = sEventsNoAl90;

}

void BeAnalysis::SetRatiosGStoEX() {

	sRatioAl0 = "sRatio>0.051 && sRatio<0.101";
	//	TCut sRatioAl0 = "sRatio>0.060 && sRatio<0.110";
	sRatioNoAl0 = "sRatio>0.050 && sRatio<0.100";
	sRatioAl180 = "sRatio>0.047 && sRatio<0.097";
	sRatioNoAl180 = "sRatio>0.048 && sRatio<0.098";
	sRatioAl90 = "sRatio>0.051 && sRatio<0.101";
	sRatioNoAl90 = "sRatio>0.051 && sRatio<0.101";

	sRatio[0] = sRatioAl0;
	sRatio[1] = sRatioNoAl0;
	sRatio[2] = sRatioAl180;
	sRatio[3] = sRatioNoAl180;
	sRatio[4] = sRatioAl90;
	sRatio[5] = sRatioNoAl90;

}

void BeAnalysis::SetCMAngularCuts() {
	cAngles = "fBeThetaCM1>60*TMath::DegToRad() && fBeThetaCM1<75*TMath::DegToRad()";
	crAngles = "f6BeThetaCM1>60*TMath::DegToRad() && f6BeThetaCM1<75*TMath::DegToRad()";
}

void BeAnalysis::Spectra() {

	if (!spectra) {
		return;
	}

	if (!che) return;

	TCanvas *cSpectra = new TCanvas();
//	TString canvasTitle;

	TH1F *hsSpectra[6];
	TH1F *heSpectra[6];
	TH1F *hSdiff[6];

	canvasTitle.Form("energy spectra fitting;\t(%d,%d) degrees", kMinAngle, kMaxAngle);
	cSpectra->SetTitle(canvasTitle.Data());
	cSpectra->Divide(2, 3);

	che->SetLineColor(kGreen);

	for (Int_t j = 0; j < 6; j++) {	//different files
		if (!kChains[j]) continue;
		cout << kChains[j] << endl;
		cSpectra->cd(j+1);
		chs[j]->SetLineColor(kBlue);
		chs[j]->SetFillColor(kWhite);
		hsName.Form("hsSpectra%d", j);
		drawCommand.Form("f6BeIM>>%s(200,0,10)", hsName.Data());
		//			chs[j]->Draw(drawCommand.Data(), cQ && crBeWork && crAngles && sRatio[j], "", sEvents[j]);
		chs[j]->Draw(drawCommand.Data(), cProtons && cQ && crBeWork && crAngles && sRatio[j], "", sEvents[j]);
		hsSpectra[j] = (TH1F*)gPad->FindObject(hsName.Data());
		heName.Form("heSpectra%d", j);
		drawCommand.Form("fBeIM>>%s(200,0,10)", heName.Data());
		//			che->Draw(drawCommand.Data(), cQ && cBeWork && cAngles, "same", eEvents[j]);
		che->Draw(drawCommand.Data(), cProtons && cQ && cBeWork && cAngles, "", eEvents[j]);
		heSpectra[j] = (TH1F*)gPad->FindObject(heName.Data());

		//			heSpectra[j]->Draw("");
		heSpectra[j]->Draw("E");
		heSpectra[j]->SetLineColor(1);
		heSpectra[j]->SetTitle("");
		heSpectra[j]->SetXTitle("E_{T} [MeV]");
		heSpectra[j]->GetXaxis()->SetTitleOffset(0.95);
		heSpectra[j]->GetXaxis()->CenterTitle();
		heSpectra[j]->SetYTitle("counts");
		heSpectra[j]->GetYaxis()->SetTitleOffset(0.7);
		heSpectra[j]->GetYaxis()->CenterTitle();
		hsSpectra[j]->SetLineColor(kGray+1);
		hsSpectra[j]->SetFillColor(kGray+1);
		hsSpectra[j]->Draw("same");
		heSpectra[j]->Draw("same E");

		if (kAutoRange) {
			Float_t leftMaxMC = kRangeProportion*hsSpectra[j]->GetMaximum();
			Float_t leftMaxE = kRangeProportion*heSpectra[j]->GetMaximum();
			//				hsSpectra[j]->GetYaxis()->SetRangeUser(0, leftMaxE);
			leftMaxMC > leftMaxE ? heSpectra[j]->GetYaxis()->SetRangeUser(0, leftMaxMC) : heSpectra[j]->GetYaxis()->SetRangeUser(0, leftMaxE);
		}


		hSdiff[j] = new TH1F(*heSpectra[j]);
		hSdiff[j]->Add(hsSpectra[j], -1);
		hSdiff[j]->SetLineColor(kRed);
		//			hSdiff[j]->Draw("same");



		cSpectra->Update();

		cout << endl;
		Info("statesRatioFitting", "MC %d", j+1);
		//			Info("statesRatioFitting", "ground state: %3.1f (exp)/ %3.1f (sim) = %3.3f", heSpectra[j]->Integral(0,20), hsSpectra[j]->Integral(0,20), heSpectra[j]->Integral(0,20)/hsSpectra[j]->Integral(0,20) );
		//			Info("statesRatioFitting", "left slope: %3.1f (exp)/ %3.1f (sim) = %3.3f\n", heSpectra[j]->Integral(20,32), hsSpectra[j]->Integral(20,32), heSpectra[j]->Integral(20,32)/hsSpectra[j]->Integral(20,32) );
		//0-2 MeV, 2-3.1 MeV
		//			Info("statesRatioFitting", "ground state: %3.1f (exp)/ %3.1f (sim) = %3.3f", heSpectra[j]->Integral(0,40), hsSpectra[j]->Integral(0,40), heSpectra[j]->Integral(0,40)/hsSpectra[j]->Integral(0,40) );
		//			Info("statesRatioFitting", "left slope: %3.1f (exp)/ %3.1f (sim) = %3.3f\n", heSpectra[j]->Integral(40,64), hsSpectra[j]->Integral(40,64), heSpectra[j]->Integral(40,64)/hsSpectra[j]->Integral(40,64) );
		//0-2 MeV, 2.5-3.5 MeV
		Info("statesRatioFitting", "ground state: %3.1f (sim)/ %3.1f (exp) = %3.3f", hsSpectra[j]->Integral(0,40), heSpectra[j]->Integral(0,40), hsSpectra[j]->Integral(0,40)/heSpectra[j]->Integral(0,40) );
		Info("statesRatioFitting", "left slope: %3.1f (sim)/ %3.1f (exp) = %3.3f\n", hsSpectra[j]->Integral(50,70), heSpectra[j]->Integral(50,70), hsSpectra[j]->Integral(50,70)/heSpectra[j]->Integral(50,70) );
		//0-1.3 MeV, 2-3.1 MeV
		//			Info("statesRatioFitting", "ground state: %3.1f (exp)/ %3.1f (sim) = %3.3f", heSpectra[j]->Integral(0,26), hsSpectra[j]->Integral(0,26), heSpectra[j]->Integral(0,26)/hsSpectra[j]->Integral(0,26) );
		//			Info("statesRatioFitting", "left slope: %3.1f (exp)/ %3.1f (sim) = %3.3f\n", heSpectra[j]->Integral(40,64), hsSpectra[j]->Integral(40,64), heSpectra[j]->Integral(40,64)/hsSpectra[j]->Integral(40,64) );
	}//for j

	/*if (savePictures) {
			canvasName.Form("%sfigSpectra:%s%s", ppath.Data(), configuration.Data(), pictFormat.Data());
			cout << canvasName.Data() << endl;
			cSpectra->SaveAs(canvasName.Data());
			cSpectra->cd();
			cSpectra->Close();
		}//if save pictures*/


	Info("sfAngInt_spectra.cxx", "Finished.");

}

void BeAnalysis::EpsilonT() {


	TCanvas *canEpsilonT[noIntervals];
	for (Int_t i = 0; i<noIntervals; i++) {	//energy intervals (canvases)

		TH1F *hseT[noIntervals][6];
		TH1F *heeT[noIntervals][6];
		TH1F *hieT[noIntervals][6];

		if (!epsilonT[i]) continue;
//		{

			canEpsilonT[i] = new TCanvas();

			canvasTitle.Form("epsilon in %s T-system;\t(%d,%d) degrees", (const char*)cBeE[i], kMinAngle, kMaxAngle);
			canEpsilonT[i]->SetTitle(canvasTitle.Data());
			canEpsilonT[i]->Divide(2, 3);

			che->SetLineColor(kBlack);
			for (Int_t j = 0; j < 6; j++) {	//different files
				if (!kChains[j]) continue;
				canEpsilonT[i]->cd(j+1);
				chs[j]->SetLineColor(kGray+1);
				chs[j]->SetFillColor(kGray+1);
				hsName.Form("hseT%d_%d", i, j);
				//				drawCommand.Form("fTpp/f6BeIM>>%s(50,0,1)", hsName.Data());
				drawCommand.Form("fTpp/f6BeIM>>%s(30,0,1)", hsName.Data());
				//				chs[j]->Draw(drawCommand.Data(), cQ && crBeE[i] && crAngles && sRatio[j], "", sEventsECuts[i][j]);
				chs[j]->Draw(drawCommand.Data(), cProtons && cQ && crBeE[i] && crAngles && sRatio[j], "", sEventsECuts[i][j]);
				hseT[i][j] = (TH1F*)gPad->FindObject(hsName.Data());
				heName.Form("heeT%d_%d", i, j);
				//				drawCommand.Form("fTpp/fBeIM>>%s(50,0,1)", heName.Data());
				drawCommand.Form("fTpp/fBeIM>>%s", heName.Data());
				//				che->Draw(drawCommand.Data(), cQ && cBeE[i] && cAngles, "E same", eEventsECuts[i][j]);
				che->Draw(drawCommand.Data(), cProtons && cQ && cBeE[i] && cAngles, "E same", eEventsECuts[i][j]);
				heeT[i][j] = (TH1F*)gPad->FindObject(heName.Data());

				hseT[i][j]->Draw();
				hseT[i][j]->SetTitle("");
				hseT[i][j]->SetXTitle("\\varepsilon");
				hseT[i][j]->GetXaxis()->SetTitleOffset(0.95);
				hseT[i][j]->GetXaxis()->CenterTitle();
				hseT[i][j]->SetYTitle("counts");
				//				hseT[i][j]->GetYaxis()->SetTitleOffset(1.09);
				hseT[i][j]->GetYaxis()->SetTitleOffset(0.7);
				hseT[i][j]->GetYaxis()->CenterTitle();
				heeT[i][j]->Draw("E same");
				if (kAutoRange) {
					Float_t leftMaxMC = kRangeProportion*hseT[i][j]->GetMaximum();
					Float_t leftMaxE = kRangeProportion*heeT[i][j]->GetMaximum();
					hseT[i][j]->GetYaxis()->SetRangeUser(0, leftMaxE);
					leftMaxMC > leftMaxE ? hseT[i][j]->GetYaxis()->SetRangeUser(0, leftMaxMC) : hseT[i][j]->GetYaxis()->SetRangeUser(0, leftMaxE);
				}
				else {
					hseT[i][j]->GetYaxis()->SetRangeUser(0, eTRange[i][j]);
				}
				canEpsilonT[i]->Update();


				hiName.Form("hieT%d_%d", i, j);
				//				drawCommand.Form("sTpp/E_IM>>%s(50,0,1)", hiName.Data());
				drawCommand.Form("sTpp/E_IM>>%s", hiName.Data());
				ti[j]->Draw(drawCommand.Data(), ciBeE[i] && sRatio[j], "same");
				hieT[i][j] = (TH1F*)gPad->FindObject(hiName.Data());
				Float_t rightmax = kRangeProportionIn*hieT[i][j]->GetMaximum();
				Float_t scale = canEpsilonT[i]->GetPad(j+1)->GetUymax()/rightmax;
				hieT[i][j]->SetLineColor(kRed);
				hieT[i][j]->Scale(scale);

				//draw an axis on the right side
				TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
						gPad->GetUxmax(), gPad->GetUymax(),0,rightmax,510,"+L");
				axis->SetLineColor(kRed);
				axis->SetLabelColor(kRed);
				axis->Draw();

				canEpsilonT[i]->Update();

				if (kVerbose) {
					Info("sfAngInt_epsilonT.cxx", "cut%d; case%d: %3.1f (exp)/ %3.1f (sim) = %3.3f",
							i, j, heeT[i][j]->Integral(0,heeT[i][j]->GetNbinsX()),
							hseT[i][j]->Integral(0,hseT[i][j]->GetNbinsX()),
							heeT[i][j]->Integral(0,heeT[i][j]->GetNbinsX())/hseT[i][j]->Integral(0,hseT[i][j]->GetNbinsX()) );
				}

			}//for j
			if (kVerbose) cout << endl;

			/*if (savePictures) {
				canvasName.Form("%sfig%depsilonT:%s%s", ppath.Data(), i, configuration.Data(), pictFormat.Data());
				cout << canvasName.Data() << endl;
				canEpsilonT[i]->SaveAs(canvasName.Data());
				canEpsilonT[i]->cd();
				canEpsilonT[i]->Close();
			}//if save pictures*/
//		}//if
	}//for i

}

void BeAnalysis::CosThetaTk() {

	TCanvas *cThetaT[noIntervals];
	for (Int_t i = 0; i<noIntervals; i++) {	//energy intervals (canvases)

		TH1F *hscoskT[noIntervals][6];
		TH1F *hecoskT[noIntervals][6];
		TH1F *hicoskT[noIntervals][6];

		if (!cosThetaT[i]) continue;
		//		{

		cThetaT[i] = new TCanvas();

		canvasTitle.Form("cos thetak in %s T-system;\t(%d,%d) degrees", (const char*)cBeE[i], kMinAngle, kMaxAngle);
		cThetaT[i]->SetTitle(canvasTitle.Data());
		cThetaT[i]->Divide(2, 3);

		che->SetLineColor(kBlack);
		for (Int_t j = 0; j < 6; j++) {	//different files
			if (!kChains[j]) continue;
			cThetaT[i]->cd(j+1);
			chs[j]->SetLineColor(kGray+1);
			chs[j]->SetFillColor(kGray+1);
			hsName.Form("hscoskT%d_%d", i, j);
			drawCommand.Form("fCosThetaTk>>%s(50,-1,1)", hsName.Data());
			chs[j]->Draw(drawCommand.Data(), cQ && crBeE[i] && crAngles && sRatio[j], "", sEventsECuts[i][j]);
			hscoskT[i][j] = (TH1F*)gPad->FindObject(hsName.Data());
			heName.Form("hecoskT%d_%d", i, j);
			drawCommand.Form("fCosThetaTk>>%s(50,-1,1)", heName.Data());
			che->Draw(drawCommand.Data(), cQ && cBeE[i] && cAngles, "E same", eEventsECuts[i][j]);
			hecoskT[i][j] = (TH1F*)gPad->FindObject(heName.Data());

			hscoskT[i][j]->Draw();
			//				hscoskT[i][j]->GetYaxis()->SetRangeUser(0, coskTRange[i][j]);
			hecoskT[i][j]->Draw("E same");
			if (kAutoRange) {
				Float_t leftMaxMC = kRangeProportion*hecoskT[i][j]->GetMaximum();
				Float_t leftMaxE = kRangeProportion*hscoskT[i][j]->GetMaximum();
				hscoskT[i][j]->GetYaxis()->SetRangeUser(0, leftMaxE);
				leftMaxMC > leftMaxE ? hscoskT[i][j]->GetYaxis()->SetRangeUser(0, leftMaxMC) : hscoskT[i][j]->GetYaxis()->SetRangeUser(0, leftMaxE);
			}
			else {
				hscoskT[i][j]->GetYaxis()->SetRangeUser(0, eTRange[i][j]);
			}
			cThetaT[i]->Update();


			hiName.Form("hicoskT%d_%d", i, j);
			drawCommand.Form("sCosThetaTk>>%s(50,0,1)", hiName.Data());
			ti[j]->Draw(drawCommand.Data(), ciBeE[i] && sRatio[j], "same");
			hicoskT[i][j] = (TH1F*)gPad->FindObject(hiName.Data());
			Float_t rightmax = 1.1*hicoskT[i][j]->GetMaximum();
			Float_t scale = cThetaT[i]->GetPad(j+1)->GetUymax()/rightmax;
			hicoskT[i][j]->SetLineColor(kRed);
			hicoskT[i][j]->Scale(scale);

			//draw an axis on the right side
			TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
					gPad->GetUxmax(), gPad->GetUymax(),0,rightmax,510,"+L");
			axis->SetLineColor(kRed);
			axis->SetLabelColor(kRed);
			axis->Draw();

			cThetaT[i]->Update();

			if (kVerbose) {
				Info("sfAngInt_cosThetaTk.cxx", "cut%d; case%d: %3.1f (exp)/ %3.1f (sim) = %3.3f",
						i, j, hecoskT[i][j]->Integral(0,hecoskT[i][j]->GetNbinsX()),
						hscoskT[i][j]->Integral(0,hscoskT[i][j]->GetNbinsX()),
						hecoskT[i][j]->Integral(0,hecoskT[i][j]->GetNbinsX())/hscoskT[i][j]->Integral(0,hscoskT[i][j]->GetNbinsX()) );
			}

		}//for j

		/*if (savePictures) {
			canvasName.Form("%sfig%dcosThetakT:%s%s", ppath.Data(), i, configuration.Data(), pictFormat.Data());
			cThetaT[i]->SaveAs(canvasName.Data());
			cThetaT[i]->cd();
			cThetaT[i]->Close();
		}*/
		//		}//if
	}//for i

}

void BeAnalysis::EpsilonY() {

	TCanvas *canEpsilonY[noIntervals];
	for (Int_t i = 0; i<noIntervals; i++) {	//energy intervals (canvases)

		TH1F *hseY[noIntervals][6];
		TH1F *heeY[noIntervals][6];
		TH1F *hieY[noIntervals][6];

		if (!epsilonY[i]) continue;
		//		{

		canEpsilonY[i] = new TCanvas();

		canvasTitle.Form("epsilon in %s Y-system;\t(%d,%d) degrees", (const char*)cBeE[i], kMinAngle, kMaxAngle);
		canEpsilonY[i]->SetTitle(canvasTitle.Data());
		canEpsilonY[i]->Divide(2, 3);

		che->SetLineColor(kBlack);
		for (Int_t j = 0; j < 6; j++) {	//different files
			canEpsilonY[i]->cd(j+1);
			chs[j]->SetLineColor(kGray+1);
			chs[j]->SetFillColor(kGray+1);
			hsName.Form("hseY%d_%d", i, j);
			drawCommand.Form("fTap/f6BeIM>>%s(50,0,1)", hsName.Data());
			chs[j]->Draw(drawCommand.Data(), cQ && crBeE[i] && crAngles && sRatio[j], "", sEvents[j]);
			hseY[i][j] = (TH1F*)gPad->FindObject(hsName.Data());
			heName.Form("heeY%d_%d", i, j);
			drawCommand.Form("fTap/fBeIM>>%s(50,0,1)", heName.Data());
			che->Draw(drawCommand.Data(), cQ && cBeE[i] && cAngles, "E same", eEvents[j]);
			heeY[i][j] = (TH1F*)gPad->FindObject(heName.Data());

			hseY[i][j]->Draw();
			heeY[i][j]->Draw("E same");
			if (kAutoRange) {
				Float_t leftMaxMC = kRangeProportion*heeY[i][j]->GetMaximum();
				Float_t leftMaxE = kRangeProportion*hseY[i][j]->GetMaximum();
				hseY[i][j]->GetYaxis()->SetRangeUser(0, leftMaxE);
				leftMaxMC > leftMaxE ? hseY[i][j]->GetYaxis()->SetRangeUser(0, leftMaxMC) : hseY[i][j]->GetYaxis()->SetRangeUser(0, leftMaxE);
			}
			else {
				hseY[i][j]->GetYaxis()->SetRangeUser(0, eTRange[i][j]);
			}
			canEpsilonY[i]->Update();



			//				ti[j]->Draw("sTap/E_IM", ciBeE[i] && sRatio[j], "same", iEvents[j]);

			hiName.Form("hieY%d_%d", i, j);
			drawCommand.Form("sTap/E_IM>>%s(50,0,1)", hiName.Data());
			ti[j]->Draw(drawCommand.Data(), ciBeE[i] && sRatio[j], "same");
			hieY[i][j] = (TH1F*)gPad->FindObject(hiName.Data());
			Float_t rightmax = 1.1*hieY[i][j]->GetMaximum();
			Float_t scale = canEpsilonY[i]->GetPad(j+1)->GetUymax()/rightmax;
			hieY[i][j]->SetLineColor(kRed);
			hieY[i][j]->Scale(scale);

			//draw an axis on the right side
			TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
					gPad->GetUxmax(), gPad->GetUymax(),0,rightmax,510,"+L");
			axis->SetLineColor(kRed);
			axis->SetLabelColor(kRed);
			axis->Draw();



			canEpsilonY[i]->Update();

			if (kVerbose) {
				Info("sfAngInt_epsilonY.cxx", "cut%d; case%d: %3.1f (exp)/ %3.1f (sim) = %3.3f",
						i, j, heeY[i][j]->Integral(0,heeY[i][j]->GetNbinsX()),
						hseY[i][j]->Integral(0,hseY[i][j]->GetNbinsX()),
						heeY[i][j]->Integral(0,heeY[i][j]->GetNbinsX())/hseY[i][j]->Integral(0,hseY[i][j]->GetNbinsX()) );
			}

		}//for j

		/*if (savePictures) {
			canvasName.Form("%sfig%dEpsilonY:%s%s", ppath.Data(), i, configuration.Data(), pictFormat.Data());
			cout << canvasName.Data() << endl;
			canEpsilonY[i]->SaveAs(canvasName.Data());
			canEpsilonY[i]->cd();
			canEpsilonY[i]->Close();
		}//if save pictures*/
		//		}//if
	}//for i


}

void BeAnalysis::CosThetaYk() {
}

void BeAnalysis::ThetaP1T() {
}

void BeAnalysis::ExpEventsECuts() {
	eEventsECuts[0][0] = 3950000;
	eEventsECuts[0][1] = 3950000;
	eEventsECuts[0][2] = 3800000;
	eEventsECuts[0][3] = 3850000;
	eEventsECuts[0][4] = 3950000;
	eEventsECuts[0][5] = 3900000;

	eEventsECuts[1][0] = eMaxEvents;
	eEventsECuts[1][1] = eMaxEvents;
	eEventsECuts[1][2] = 4100000;
	eEventsECuts[1][3] = eMaxEvents;
	eEventsECuts[1][4] = eMaxEvents;
	eEventsECuts[1][5] = eMaxEvents;

	eEventsECuts[2][0] = 3050000;
	eEventsECuts[2][1] = 3150000;
	eEventsECuts[2][2] = 3200000;
	eEventsECuts[2][3] = 3200000;
	eEventsECuts[2][4] = 3250000;
	eEventsECuts[2][5] = 3350000;

	eEventsECuts[3][0] = 3850000;
	eEventsECuts[3][1] = 3950000;
	eEventsECuts[3][2] = 3900000;
	eEventsECuts[3][3] = 3900000;
	eEventsECuts[3][4] = 4000000;
	eEventsECuts[3][5] = 4000000;

	eEventsECuts[4][0]	= 4100000;
	eEventsECuts[4][1]	= 4100000;
	eEventsECuts[4][2]	= 4000000;
	eEventsECuts[4][3]	= 4100000;
	eEventsECuts[4][4]	= eMaxEvents;
	eEventsECuts[4][5]	= eMaxEvents;
}

void BeAnalysis::SimEventsECuts() {

	for (Int_t i = 0; i < 5; i++) {
		for (Int_t j = 0; j < 6; j++) {
			sEventsECuts[i][j] = sMaxEvents;
		}
	}

	sEventsECuts[1][1] = 2630000;
	sEventsECuts[1][4] = 2660000;
	sEventsECuts[1][5] = 2660000;

}

void BeAnalysis::EpsilonTRange() {

	for (Int_t i = 0; i < 6; i++) {
		eTRange[0][i] = 300;
		eTRange[1][i] = 380;
		eTRange[2][i] = 300;
		eTRange[3][i] = 500;
		eTRange[4][i] = 500;
	}


}

void BeAnalysis::ThetaAT() {
}