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

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

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

	kVerbose = 1;
	kRangeProportionIn = 2.0;

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

}

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


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

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

}

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

	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;
	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];
//		cout << eEvents[i] << endl;
	}
	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;

//	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() {

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

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

////		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", 1);

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

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

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

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

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

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

			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

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

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

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

		che->SetLineColor(kBlack);
		for (Int_t j = 0; j < 6; j++) {	//different files
			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();
			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("sfAngInt_cosThetaYk.cxx", "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);


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

}

void BeAnalysis::ThetaAT() {

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

		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", 50, 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", 50, 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(50,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("sfAngInt_thetaAT.cxx", "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);
		/*if (savePictures) {
					canvasName.Form("%sfig%dthetaAT:%s%s", ppath.Data(), i, configuration.Data(), pictFormat.Data());
					cThetaAT[i]->SaveAs(canvasName.Data());
					cThetaAT[i]->cd();
					cThetaAT[i]->Close();
				}*/
	}//for i

}

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::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::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());
	canvas->SaveAs(canvasName.Data());
	canvas->cd();
	canvas->Close();


}