BeAnalysis.cpp

/*
 * 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), lowExpFile(0), upExpFile(0), spectra(0) {
	// TODO Auto-generated constructor stub

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

	SetCuts();
	SetChainsToDraw();

	SetNoBinsSpectra();

	noIntervals = 5;

	epsilonT = 0;
	cosThetaT = 0;
	epsilonY = 0;
	cosThetaY = 0;
	thetaAT = 0;

	kVerbose = 1;
	kRangeProportionIn = 2.0;

	kSaveFigures = 0;
	figurePath = "";
	figureFormat = "";

}

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

void BeAnalysis::OpenSimChains() {

	for (Int_t i = 0; i < 6; i++) {
		chs[i] = BeWork::OpenChain(simFiles[i].Data(),
				lowSimFile[i], upSimFile[i], "simbe", 4, "sbeam");
		Info("statesRatioFitting.cxx", "%lld events in chain \"%s\" containing simulated data", chs[i]->GetEntries(), chs[i]->GetName());
		ti[i] = BeWork::OpenTree(inputFiles[i].Data(), "sbeam", 2);
	}

}

void BeAnalysis::OpenExpChain() {

	//experimental chain
	che = BeWork::OpenChain(expFiles.Data(), lowExpFile, upExpFile, "beonly");										//original file
	Info("BeAnalysis::OpenExpChain", "%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";

	/*TCut */cBe0_14 = "fBeIM>0 && fBeIM<1.4";
	//cBe0_14 = "fBeIM>0 && fBeIM<1.6";
	//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";
	//cEpsilonT = "fTpp/fBeIM<1.";
	//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";
	//crBe0_14 = "f6BeIM>0 && f6BeIM<1.6";
	//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


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

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

	ciEpsilon = "sTpp/E_IM<0.2";
	//ciEpsilon = "sTpp/E_IM<1.";
	//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";
	//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;
	SetCMAngularCuts();
}

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(Long64_t *noExpEvents) {

	eMaxEvents = 5000000;

	if (!noExpEvents) {
		for (Int_t i = 0; i < 6; i++) {
			eEvents[i] = eMaxEvents;
		}
		Warning("BeAnalysis::SetNoExpEvents", "Default numbers of experimental events were set.");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
		eEvents[i] = noExpEvents[i];
	}

	return;

}

void BeAnalysis::SetNoSimEvents(Long64_t *noSimEvents) {

	sMaxEvents = 3000000;

	if (!noSimEvents) {
		for (Int_t i = 0; i < 6; i++) {
			sEvents[i] = sMaxEvents;
			cout << sEvents[i] << endl;
		}
		Warning("BeAnalysis::SetNoSimEvents", "Default numbers of simulated events were set.");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
		sEvents[i] = noSimEvents[i];
		cout << sEvents[i] << endl;
	}
	return;

}

void BeAnalysis::SetRatiosGStoEX(TString sRatioAl0, TString sRatioNoAl0,
		TString sRatioAl180, TString sRatioNoAl180,
		TString sRatioAl90, TString sRatioNoAl90) {

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

}

void BeAnalysis::SetCMAngularCuts() {

	TString sAngles, srAngles;

	sAngles.Form("fBeThetaCM1>%d*TMath::DegToRad() && fBeThetaCM1<%d*TMath::DegToRad()", kMinAngle, kMaxAngle);
	srAngles.Form("f6BeThetaCM1>%d*TMath::DegToRad() && f6BeThetaCM1<%d*TMath::DegToRad()", kMinAngle, kMaxAngle);

	cAngles = sAngles.Data();
	crAngles = srAngles.Data();
}

void BeAnalysis::Spectra() {

	if (!spectra) {
		return;
	}

	if (!che) {
		Error("BeAnalysis::Spectra", "Chain with experimental data was not initialized.");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
			if (!chs[i]) {
				Error("BeAnalysis::Spectra", "Chain No. %d with simulated data was not initialized.", i);
				return;
			}
	}

	TCanvas *cSpectra = new TCanvas();

	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);
	CanvasDivision(cSpectra);

	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);
		hsSpectra[j] = new TH1F(hsName.Data(), "title", noBinsSpectra, 0, 10);
//		drawCommand.Form("f6BeIM>>%s(200,0,10)", hsName.Data());
		drawCommand.Form("f6BeIM>>%s", 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);
		heSpectra[j] = new TH1F(heName.Data(), "title", noBinsSpectra, 0, 10);
//		drawCommand.Form("fBeIM>>%s(200,0,10)", heName.Data());
		drawCommand.Form("fBeIM>>%s", heName.Data());
		//			che->Draw(drawCommand.Data(), cQ && cBeWork && cAngles, "same", eEvents[j]);

		che->Draw(drawCommand.Data(), cProtons && cQ && cBeWork && cAngles, "", eEvents[j]);

////		cout << drawCommand << "\t" << eEvents[j] << endl;
////		che->Draw(drawCommand.Data(), cProtons && cQ && cBeWork && cAngles, "", eEvents[j]);
//		che->Draw("fBeIM", "", "");
//		continue;
//		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();
			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);

		cSpectra->Update();

		cout << endl;
		Info("statesRatioFitting", "MC %d", j+1);
		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) );
	}//for j

	SaveFigures(cSpectra, "Spectra", 0);

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

}

void BeAnalysis::EpsilonT() {

	if(!epsilonT) {
		Error("BeAnalysis::EpsilonT", "Energy intervals were not set.");
		return;
	}

	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, 0., 0.);
		CanvasDivision(canEpsilonT[i]);

//		cProtons.Print();
//		cQ.Print();
//		crBeE[i].Print();
//		crAngles.Print();

//		TFile fw("histos.root", "RECREATE");

		che->SetLineColor(kBlack);
		for (Int_t j = 0; j < 6; j++) {	//different files
			if (!kChains[j]) continue;

//			sRatio[j].Print();
//			cout << sEventsECuts[i][j] << endl;
//			cout << eEventsECuts[i][j] << endl;

			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(30,0,1)", hsName.Data());
			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", heName.Data());
			che->Draw(drawCommand.Data(), cProtons && cQ && cBeE[i] && cAngles, "E same", eEventsECuts[i][j]);
			heeT[i][j] = (TH1F*)gPad->FindObject(heName.Data());
//			continue;

			hseT[i][j]->Draw();
			hseT[i][j]->SetTitle("");
			hseT[i][j]->SetXTitle("\\varepsilon");
//			hseT[i][j]->GetXaxis()->SetTitleOffset(0.95);
//			hseT[i][j]->GetXaxis()->SetTitleSize(0.11);
			hseT[i][j]->GetXaxis()->CenterTitle();
			hseT[i][j]->SetYTitle("counts");
//			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", 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("BeAnalysis::EpsilonT", "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()) );
			}

//			hseT[i][j]->Write();

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

//		fw.Close();

		SaveFigures(canEpsilonT[i], "epsilonT", i);

	}//for i

}

void BeAnalysis::CosThetaTk() {

	if(!cosThetaT) {
		Error("BeAnalysis::CosThetaTk", "Energy intervals were not set.");
		return;
	}

	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);
		CanvasDivision(cThetaT[i]);

		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(30,-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(30,-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]->SetTitle("");
			hscoskT[i][j]->SetXTitle("\\cos \\theta_k");
//			hseT[i][j]->GetXaxis()->SetTitleOffset(0.95);
//			hseT[i][j]->GetXaxis()->SetTitleSize(0.11);
			hscoskT[i][j]->GetXaxis()->CenterTitle();
			hscoskT[i][j]->SetYTitle("counts");
//			hseT[i][j]->GetYaxis()->SetTitleOffset(0.7);
			hscoskT[i][j]->GetYaxis()->CenterTitle();



			//				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(30,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("BeAnalysis::CosThetaTk", "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

		SaveFigures(cThetaT[i], "cosThetakT", i);

	}//for i

}

void BeAnalysis::EpsilonY() {

	if(!epsilonY) {
		Error("BeAnalysis::EpsilonY", "Energy intervals were not set.");
		return;
	}

	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);
		CanvasDivision(canEpsilonY[i]);

		che->SetLineColor(kBlack);
		for (Int_t j = 0; j < 6; j++) {	//different files
			if (!kChains[j]) continue;
			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();
			hseY[i][j]->SetTitle("");
			hseY[i][j]->SetXTitle("\\varepsilon");
//			hseT[i][j]->GetXaxis()->SetTitleOffset(0.95);
//			hseT[i][j]->GetXaxis()->SetTitleSize(0.11);
			hseY[i][j]->GetXaxis()->CenterTitle();
			hseY[i][j]->SetYTitle("counts");
//			hseT[i][j]->GetYaxis()->SetTitleOffset(0.7);
			hseY[i][j]->GetYaxis()->CenterTitle();

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

			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("BeAnalysis::EpsilonY", "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

		SaveFigures(canEpsilonY[i], "EpsilonY", i);

	}//for i


}

void BeAnalysis::CosThetaYk() {

	if(!cosThetaY) {
		Error("BeAnalysis::CosThetaYk", "Energy intervals were not set.");
		return;
	}

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

		TH1F *hscoskY[noIntervals][6];
		TH1F *hecoskY[noIntervals][6];
		TH1F *hicoskY[noIntervals][6];

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

		cThetaY[i] = new TCanvas();

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

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

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

			cThetaY[i]->Update();

			hiName.Form("hicoskY%d_%d", i, j);
			drawCommand.Form("sCosThetaYk>>%s", hiName.Data());
			ti[j]->Draw(drawCommand.Data(), ciBeE[i] && sRatio[j], "same");
			hicoskY[i][j] = (TH1F*)gPad->FindObject(hiName.Data());
			Float_t rightmax = 1.1*hicoskY[i][j]->GetMaximum();
			Float_t scale = cThetaY[i]->GetPad(j+1)->GetUymax()/rightmax;
			hicoskY[i][j]->SetLineColor(kRed);
			hicoskY[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();

			cThetaY[i]->Update();

			if (kVerbose) {
				Info("BeAnalysis::CosThetaYk", "cut%d; case%d: %3.1f (exp)/ %3.1f (sim) = %3.3f",
						i, j, hecoskY[i][j]->Integral(0,hecoskY[i][j]->GetNbinsX()),
						hscoskY[i][j]->Integral(0,hscoskY[i][j]->GetNbinsX()),
						hecoskY[i][j]->Integral(0,hecoskY[i][j]->GetNbinsX())/hscoskY[i][j]->Integral(0,hscoskY[i][j]->GetNbinsX()) );
			}


		}//for j

		SaveFigures(cThetaY[i], "cosThetakY", i);

	}//for i

}

void BeAnalysis::ThetaAT() {

	if(!thetaAT) {
		Error("BeAnalysis::ThetaAT", "Energy intervals were not set.");
		return;
	}

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

		TH1F *hsthetaAT[noIntervals][6];
		TH1F *hethetaAT[noIntervals][6];
		TH1F *hithetaAT[noIntervals][6];

		if (!thetaAT[i]) continue;

		cThetaAT[i] = new TCanvas();

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

		che->SetLineColor(kBlack);
		for (Int_t j = 0; j < 6; j++) {	//different files
			if (!kChains[j]) continue;
			cThetaAT[i]->cd(j+1);
			hsName.Form("hsthetaAT%d_%d", i, j);
			hsthetaAT[i][j] = new TH1F(hsName.Data(), "title", 30, 0, 3.14);
			hsthetaAT[i][j]->SetLineColor(kGray+1);
			hsthetaAT[i][j]->SetFillColor(kGray+1);
			drawCommand.Form("fAThetaCM>>%s", hsthetaAT[i][j]->GetName());
			chs[j]->Draw(drawCommand.Data(), cQ && crBeE[i] && crAngles && crEpsilonT && sRatio[j], "goff", sEventsECuts[i][j]);
			heName.Form("hethetaAT%d_%d", i, j);
			hethetaAT[i][j] = new TH1F(heName.Data(), "title", 30, 0, 3.14);
			drawCommand.Form("fAThetaCM>>%s", hethetaAT[i][j]->GetName());
			che->Draw(drawCommand.Data(), cQ && cBeE[i] && cAngles && cEpsilonT, "goff", eEventsECuts[i][j]);

			hsthetaAT[i][j]->Draw();
			hsthetaAT[i][j]->SetTitle("");
			hsthetaAT[i][j]->SetXTitle("\\theta_{\\alpha} (rad)");
//			hsthetaAT[i][j]->GetXaxis()->SetTitleOffset(0.95);
			hsthetaAT[i][j]->GetXaxis()->CenterTitle();
			hsthetaAT[i][j]->SetYTitle("counts");
			//				hsthetaAT[i][j]->GetYaxis()->SetTitleOffset(1.09);
//			hsthetaAT[i][j]->GetYaxis()->SetTitleOffset(0.7);
			hsthetaAT[i][j]->GetYaxis()->CenterTitle();

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

			hiName.Form("hithetaAT%d_%d", i, j);
			drawCommand.Form("sACM.Theta()>>%s(30,0,3.14)", hiName.Data());
			ti[j]->Draw(drawCommand.Data(), ciEpsilon && ciBeE[i] && sRatio[j], "same");
			hithetaAT[i][j] = (TH1F*)gPad->FindObject(hiName.Data());
			Float_t rightmax = 1.1*hithetaAT[i][j]->GetMaximum();
			Float_t scale = cThetaAT[i]->GetPad(j+1)->GetUymax()/rightmax;
			hithetaAT[i][j]->SetLineColor(kRed);
			hithetaAT[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();

			cThetaAT[i]->Update();

			if (kVerbose) {
				Info("BeAnalysis::ThetaAT", "cut%d; case%d: %3.1f (exp)/ %3.1f (sim) = %3.3f",
						i, j, hethetaAT[i][j]->Integral(0,hethetaAT[i][j]->GetNbinsX()),
						hsthetaAT[i][j]->Integral(0,hsthetaAT[i][j]->GetNbinsX()),
						hethetaAT[i][j]->Integral(0,hethetaAT[i][j]->GetNbinsX())/hsthetaAT[i][j]->Integral(0,hsthetaAT[i][j]->GetNbinsX()) );
			}

		}//for j

		SaveFigures(cThetaAT[i], "thetaAT", i);

	}//for i

}

//void BeAnalysis::ExpEventsECuts(Long64_t **noExpEvents) {
//void BeAnalysis::ExpEventsECuts(Int_t (&noExpEvents)[5][6]) {
//void BeAnalysis::ExpEventsECuts(Int_t (*noExpEvents)[5][6]) {
void BeAnalysis::ExpEventsECuts(Long64_t noExpEvents[5][6]) {


//	eMaxEvents = 5000000;

	if (!noExpEvents) {
		for (Int_t i = 0; i < 5; i++) {
			for (Int_t j = 0; j < 6; j++) {
				eEventsECuts[i][j] = eMaxEvents;
			}
		}
		Warning("BeAnalysis::ExpEventsECuts", "Default numbers of experimental events were set.");
		return;
	}

	for (Int_t i = 0; i < 5; i++) {
		for (Int_t j = 0; j < 6; j++) {
			eEventsECuts[i][j] = (Long64_t)noExpEvents[i][j];
		}
	}
	return;


//	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(Long64_t noSimEvents[5][6]) {


//	sMaxEvents;

	if (!noSimEvents) {
		for (Int_t i = 0; i < 5; i++) {
			for (Int_t j = 0; j < 6; j++) {
				sEventsECuts[i][j] = sMaxEvents;
			}
		}
		Warning("BeAnalysis::SimEventsECuts", "Default numbers of experimental events were set.");
		return;
	}

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



//	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::SetFigures(TString figPath, TString figFormat,
		Bool_t kSave) {
	kSaveFigures = kSave;
	figureFormat = figFormat;
	figurePath = figPath;
}

void BeAnalysis::SetExpChain(TString files, Int_t minFnumber, Int_t maxFnumber) {

	lowExpFile = minFnumber;
	upExpFile = maxFnumber;
	expFiles = files;

}

void BeAnalysis::SetSimChains(TString chainAl0name, TString chainNoAl0name,
		TString chainAl180name, TString chainNoAl180name,
		TString chainAl90name, TString chainNoAl90name) {

	simFiles[0] = chainAl0name;
	simFiles[1] = chainNoAl0name;
	simFiles[2] = chainAl180name;
	simFiles[3] = chainNoAl180name;
	simFiles[4] = chainAl90name;
	simFiles[5] = chainNoAl90name;

}

void BeAnalysis::SetInputChains(TString inputAl0name, TString inputNoAl0name,
		TString inputAl180name, TString inputNoAl180name,
		TString inputAl90name, TString inputNoAl90name) {

	inputFiles[0] = inputAl0name;
	inputFiles[1] = inputNoAl0name;
	inputFiles[2] = inputAl180name;
	inputFiles[3] = inputNoAl180name;
	inputFiles[4] = inputAl90name;
	inputFiles[5] = inputNoAl90name;

}

void BeAnalysis::SetNoSimFiles(Int_t* minFnumber, Int_t* maxFnumber) {
	eMaxEvents = 5000000;

	if (!(minFnumber && maxFnumber)) {
		for (Int_t i = 0; i < 6; i++) {
			lowSimFile[i] = 0;
			upSimFile[i] = 14;
		}
		Info("BeAnalysis::SetNoSimFiles", "Default numbers of simfiles were set.");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
		lowSimFile[i] = minFnumber[i];
		upSimFile[i] = maxFnumber[i];
	}
	return;
}

void BeAnalysis::SetEpsilonTintervals(Bool_t intervals[6]) {

	epsilonT = new Bool_t[noIntervals];

	if (!intervals) {
		for (Int_t i = 0; i < 6; i++) {
			epsilonT[i] = 1;
		}
		Warning("BeAnalysis::SetEpsilonTintervals", "All intervals epsilonT will be drawn");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
//		cout << intervals[i] << endl;
		epsilonT[i] = intervals[i];
	}
	return;

}

void BeAnalysis::SetCosThetaTkIntervals(Bool_t intervals[6]) {

	cosThetaT = new Bool_t[noIntervals];

	if (!intervals) {
		for (Int_t i = 0; i < 6; i++) {
			cosThetaT[i] = 1;
		}
		Warning("BeAnalysis::SetCosThetaTkIntervals", "All intervals epsilonT will be drawn");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
//		cout << intervals[i] << endl;
		cosThetaT[i] = intervals[i];
	}
	return;

}

void BeAnalysis::SetEpsilonYintervals(Bool_t intervals[6]) {
	epsilonY = new Bool_t[noIntervals];

	if (!intervals) {
		for (Int_t i = 0; i < 6; i++) {
			epsilonY[i] = 1;
		}
		Warning("BeAnalysis::SetEpsilonYintervals", "All intervals epsilonT will be drawn");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
//		cout << intervals[i] << endl;
		epsilonY[i] = intervals[i];
	}
	return;
}

void BeAnalysis::SetCosThetaYkIntervals(Bool_t intervals[6]) {
	cosThetaY = new Bool_t[noIntervals];

	if (!intervals) {
		for (Int_t i = 0; i < 6; i++) {
			cosThetaY[i] = 1;
		}
		Warning("BeAnalysis::SetCosThetaYkIntervals", "All intervals epsilonT will be drawn");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
//		cout << intervals[i] << endl;
		cosThetaY[i] = intervals[i];
	}
	return;
}

void BeAnalysis::SetThetaATintervals(Bool_t intervals[6]) {
	thetaAT = new Bool_t[noIntervals];

	if (!intervals) {
		for (Int_t i = 0; i < 6; i++) {
			thetaAT[i] = 1;
		}
		Warning("BeAnalysis::SetThetaATintervals", "All intervals epsilonT will be drawn");
		return;
	}

	for (Int_t i = 0; i < 6; i++) {
//		cout << intervals[i] << endl;
		thetaAT[i] = intervals[i];
	}
	return;
}

void BeAnalysis::SaveFigures(TCanvas *canvas, TString variable, Int_t interval) {

	if (!kSaveFigures) return;

	canvasName.Form("%sfig%d%s:%d%d%s", figurePath.Data(), interval, variable.Data(), kMinAngle, kMaxAngle, figureFormat.Data());
	cout << canvasName << endl;
//	return;
	canvas->SaveAs(canvasName.Data());
	canvas->cd();
	canvas->Close();


}

void BeAnalysis::CanvasDivision(TCanvas* c) {
	c->Divide(2, 3, 0., 0.);
}