/*
 * BeEvent.cpp
 *
 *  Created on: 9.12.2009
 *      Author: Vratislav
 */

#include "BeEvent.h"

ClassImp(BeEvent);

TRandom3 BeEvent::ranE(1);
TRandom3 BeEvent::ranChoice(1);

BeEvent::BeEvent() {
//	Info("BeEvent::BeEvent", "Default constructor called");

	fCuts = NULL;
	cutAlpha1 = NULL;
	cutProton1 = NULL;
	cutAlpha2 = NULL;
	cutProton2 = NULL;

	calSi = NULL;
	calCsIa = NULL;
	calCsIp = NULL;

//	Info("BeEvent::BeEvent", "Default constructor finished");
}

BeEvent::BeEvent(const char* configfile/*, const char* cutfile,*/ /*const TString calibfilepath*/) {
	Info("BeEvent::BeEvent", "Second constructor called");

	fConfigFile = configfile;
	ReadConfigFile();

	InitDetectors();
	Info("BeEvent::BeEvent", "Beam kinetic energy was set as %f AMeV", fTBeamRec);
	Info("BeEvent::BeEvent",
			"Distances for reconstruction: zT11 = %3.0f mm zT21 = %3.0f mm", fT1ExpPos,
			fT2ExpPos);
	Info("BeEvent::BeEvent",
			"Correction for proton calibration in T1: 0-7 --> %3.2f, 8-15 --> %3.2f",
			fCorrT1_0_7_CsI_P, fCorrT1_8_15_CsI_P);
	Info("BeEvent::BeEvent",
			"Correction for alpha calibration in T1: 0-7 --> %3.2f, 8-15 --> %3.2f",
			fCorrT1_0_7_CsI_A, fCorrT1_8_15_CsI_A);
	Info("BeEvent::BeEvent",
			"Correction for proton calibration in T2: 0-7 --> %3.2f, 8-15 --> %3.2f",
			fCorrT2_0_7_CsI_P, fCorrT2_8_15_CsI_P);
	Info("BeEvent::BeEvent",
			"Correction for alpha calibration in T2: 0-7 --> %3.2f, 8-15 --> %3.2f",
			fCorrT2_0_7_CsI_A, fCorrT2_8_15_CsI_A);
	Info("BeEvent::BeEvent", "Dead layer in CsI detector in T1: %3.2f", fX13_FD);
	Info("BeEvent::BeEvent", "Dead layer in CsI detector in T2: %3.2f", fX23_FD);

//	SetCalCorrs();	//function used for CsI calibration only

	CreateSiELosses();
	CreateCsIELosses();
	CreateTargetELosses();

	fCuts = NULL;
	cutAlpha1 = NULL;
	cutProton1 = NULL;
	cutAlpha2 = NULL;
	cutProton2 = NULL;

	ReadCuts(/*cutfile*/);

	//calibration parameters
	TString calFileSi;
//	calFileSi.Form("%s/beSi.cal", fCalibFilePath.Data());
//	calFileSi.Form("%s/beSi_lowThr.cal", fCalibFilePath.Data());
	calFileSi.Form("%s/%s", fCalibFilePath.Data(), sSiCal.Data());
	calSi = new AculCalibration(calFileSi.Data());

	TString calFileCsIa;
//	calFileCsIa.Form("%s/beCsIa.cal", fCalibFilePath.Data());
//	calFileCsIa.Form("%s/beCsIa_highThrs.cal", fCalibFilePath.Data());
	calFileCsIa.Form("%s/%s", fCalibFilePath.Data(), sCsIalphaCal.Data());
	calCsIa = new AculCalibration(calFileCsIa.Data());

	TString calFileCsIp;
//	calFileCsIp.Form("%s/beCsIp.cal", fCalibFilePath.Data());
//	calFileCsIp.Form("%s/beCsIp_highThrs.cal", fCalibFilePath.Data());
	calFileCsIp.Form("%s/%s", fCalibFilePath.Data(), sCsIprotonsCal.Data());
	calCsIp = new AculCalibration(calFileCsIp.Data());

	Info("BeEvent::BeEvent", "Second constructor finished\n\n");
}

//BeEvent::BeEvent(BeEvent &itself) {
//	fCuts;			//!
//	cutAlpha1;		//!		1st telescope
//	cutProton1;	//!
//	cutAlpha2;		//!		2nd telescope
//	cutProton2;	//!
//
//	//calibration parameters
//	AculCalibration *calSi;		//!		calibration parameters for Si detectors
//	AculCalibration *calCsIa;	//!		alpha parameters for CsI detectors
//	AculCalibration *calCsIp;	//!		proton parameters for CsI detectors
//}

BeEvent::~BeEvent() {
	//there are not dynamic class members in BeEvent
//	Info("BeEvent::~BeEvent", "Destructor called");

	delete calCsIa;
	delete calCsIp;
	delete calSi;

	delete cutAlpha1;
	delete cutProton1;
	delete cutAlpha2;
	delete cutProton2;

	if (fCuts != NULL) {
		fCuts->Close();
	}
	delete fCuts;
//	Info("BeEvent::~BeEvent", "Destructor successfully finished its work");
}

void BeEvent::ReadConfigFile() {

//	TString g(generator);
	std::ifstream cfile;
	//generator file opening
	if (!fConfigFile.IsNull()) {
		cfile.open(fConfigFile.Data());
		if (!cfile.is_open()) {
			Error("BeEvent::ReadConfigFile", "Configuration file was not opened!!!");
			return;
		}//if
	}//if


	TString line;
	TString configuration;

	while (cfile.good()) {
		//read line from file
		line.ReadLine(cfile);
		if (line.Contains("//")) line.Resize(line.Index("//"));
//		line.ToLower();
		line.Append(" ");

		configuration.Append(line);
	}
//	Info("BeEvent::ReadConfigFile", "vypis");
//	cout << configuration.Data() << endl;

	ConfigDictionary CDpath(configuration.Data());
	fCutFile = CDpath.GetString("cutFile");
	fCalibFilePath = CDpath.GetString("calibFilePath");
	sSiCal = CDpath.GetString("siCal");
	sCsIprotonsCal = CDpath.GetString("CsIprotonsCal");
	sCsIalphaCal = CDpath.GetString("CsIalphaCal");

	configuration.ToLower();
	ConfigDictionary CD(configuration.Data());

//	fBeOnly = CD.GetBool("beonly");

	//geometry
	fT1ExpPos = CD.GetDouble("t1expposition");
	fT2ExpPos = CD.GetDouble("t2expposition");

	//beam
	fTBeamRec = CD.GetDouble("tbeam");
	fBeamX = CD.GetDouble("beamx");
	fBeamY = CD.GetDouble("beamy");
	fBeamX_sigma = CD.GetDouble("beamx_sigma");
	fBeamY_sigma = CD.GetDouble("beamy_sigma");

	//detector thicknesses
	fX11_FD = CD.GetDouble("t11_frontdead");
	fX11 = CD.GetDouble("t11");
	fX11_BD = CD.GetDouble("t11_backdead");

	fX12_FD = CD.GetDouble("t12_frontdead");
	fX12 = CD.GetDouble("t12");
	fX12_BD = CD.GetDouble("t12_backdead");

	//#define	XD13F	14.	//fixme original value		probably in Si equivalent
	fX13_FD = CD.GetDouble("t13_frontdead");

	fX21_FD = CD.GetDouble("t21_frontdead");
	fX21 = CD.GetDouble("t21");
	fX21_BD = CD.GetDouble("t21_backdead");

	fX22_FD = CD.GetDouble("t22_frontdead");
	fX22 = CD.GetDouble("t22");
	fX22_BD = CD.GetDouble("t22_backdead");
	//#define	XD23F	14.		//fixme original value		probably in Si equivalent
	fX23_FD = CD.GetDouble("t23_frontdead");

	//correction for calibration parameters in CsI detectors
	//telescope 1
	fCorrT1_0_7_CsI_A = CD.GetDouble("c_t1_0_7_csi_a");
	fCorrT1_8_15_CsI_A = CD.GetDouble("c_t1_8_15_csi_a");
	fCorrT1_0_7_CsI_P = CD.GetDouble("c_t1_0_7_csi_p");
	fCorrT1_8_15_CsI_P = CD.GetDouble("c_t1_8_15_csi_p");

	//telescope 2
	fCorrT2_0_7_CsI_A = CD.GetDouble("c_t2_0_7_csi_a");
	fCorrT2_8_15_CsI_A = CD.GetDouble("c_t2_8_15_csi_a");
	fCorrT2_0_7_CsI_P = CD.GetDouble("c_t2_0_7_csi_p");
	fCorrT2_8_15_CsI_P = CD.GetDouble("c_t2_8_15_csi_p");

	//thresholds
	fRSTolerance = CD.GetDouble("rstolerance");
	fProtonThr = CD.GetDouble("protonthreshold");


	//config file closing
	cfile.close();
	if (cfile.is_open()) {
		Warning("BeWork::ReadConfigFile", "File %s closing error\n", fConfigFile.Data());
	}//if

	return;

}

void BeEvent::InitDetectors() {
	//set size and position of all detectors TXX

//	Info("BeEvent::InitDetectors", "Distances for reconstruction: zT11 = %3.0f mm zT21 = %3.0f mm", fT1ExpPos, fT2ExpPos);
	fT11.InitDetector(kTRUE, 32, 32, 16, 41, fX11, fX11_FD, fT1ExpPos);
//	fT11.InitDetector(kTRUE, 32, 32, 16, 41, X11, XD11, 106);		//86; 96; 101; 106
	fT12.InitDetector(kFALSE, NOESEC, 0, 16, 41, fX12, fX12_FD + fX11_BD, fT1ExpPos + 5.);
	fT13.InitDetector(kFALSE, NOESEC, 0, 16, 41, 1000000, fX13_FD, fT1ExpPos + 10.);

	fT21.InitDetector(kTRUE, 32, 32, 16, 41, fX21, fX21_FD, fT2ExpPos);
	fT22.InitDetector(kFALSE, NOESEC, 0, 16, 41, fX22, fX22_FD + fX11_BD, fT2ExpPos + 5.);
	fT23.InitDetector(kFALSE, NOESEC, 0, 16, 41, 1000000, fX23_FD, fT2ExpPos + 10.);

	return;
}

void BeEvent::FillEvent(AculRaw *rawevent) {
	//rawevent: event with RAW data to be converted to this type of event
	//calSi: object with calibration parameters for heavy charged particles in Si detectors
	//calA: object with calibration parameters for alpha particles in CsI detectors
	//calP: object with calibration parameters for protons in CsI detectors

//	sleep(1);

	Reset();

	fTrigger = rawevent->trigger;
//	InitDetectors();				//fixme this statement is called only once and should be located in constructor

	SetChannels(rawevent);

	//energy and time calibration
	SetSiDetEnergies(rawevent);			//ok
	SetSiDetTimes(rawevent);			//ok
	SetCsIalphaDetEnergies(rawevent);	//ok
	SetCsIprotonDetEnergies(rawevent);	//ok
	SetDeltaE();						//ok

	//searching particular hits
//	SetHitsOld();
	SetHits();

	return;
}

void BeEvent::FillEventOld(AculRaw *rawevent, AculCalibration *calSi,
		AculCalibration *calA, AculCalibration *calP) {
	//rawevent: event with RAW data to be converted to this type of event
	//calSi: object with calibration parameters for heavy charged particles in Si detectors
	//calA: object with calibration parameters for alpha particles in CsI detectors
	//calP: object with calibration parameters for protons in CsI detectors

	Reset();

	fTrigger = rawevent->trigger;
//	InitDetectors();				//fixme this statement is called only once and should be located in constructor

	SetChannels(rawevent);

	//energy and time calibration
	SetSiDetEnergiesOld(rawevent, calSi);
	SetSiDetTimesOld(rawevent, calSi);
	SetCsIalphaDetEnergiesOld(rawevent, calA);
	SetCsIprotonDetEnergiesOld(rawevent, calP);
	SetDeltaE();

	//searching particular hits
	SetHits();

	return;
}

void BeEvent::SetChannels(AculRaw *rawevent) {
	//set information in channel units for each detector obtained from RAW data file

	//sectors
	for (Int_t i = 0; i < 8; i++) {
		fT11.fChSec[2 * i] = rawevent->C3[3][2 * i + 1];
		fT11.fChSec[2 * i + 1] = rawevent->C3[3][2 * i];
		fT11.fChSec[2 * i + 16] = rawevent->C3[4][2 * i + 1];
		fT11.fChSec[2 * i + 1 + 16] = rawevent->C3[4][2 * i];

		fT21.fChSec[2 * i] = rawevent->C3[9][2 * i + 1]; //cable inversion
		fT21.fChSec[2 * i + 1] = rawevent->C3[9][2 * i];
		fT21.fChSec[2 * i + 16] = rawevent->C3[8][2 * i + 1]; //cable inversion
		fT21.fChSec[2 * i + 1 + 16] = rawevent->C3[8][2 * i];
	} //for

	//rings
	for (Int_t i = 0; i < 16; i++) {
		fT11.fChRing[i] = rawevent->C3[5][i];
		fT11.fChRing[i + 16] = rawevent->C3[6][i];
		fT12.fChSec[i] = rawevent->C3[13][i];
		fT13.fChSec[i] = rawevent->C3[15][i];
		fT21.fChRing[i] = rawevent->C3[10][i];
		fT21.fChRing[i + 16] = rawevent->C3[11][i];
		fT22.fChSec[i] = rawevent->C3[14][i];
		fT23.fChSec[i] = rawevent->C3[16][i];
	} //for

	fT11.fMultChSec = fT11.GetChMultiplicity(1);
	fT11.fMultChRing = fT11.GetChMultiplicity(0);
	fT12.fMultChSec = fT12.GetChMultiplicity(1);
	fT13.fMultChSec = fT13.GetChMultiplicity(1);

	fT21.fMultChSec = fT21.GetChMultiplicity(1);
	fT21.fMultChRing = fT21.GetChMultiplicity(0);
	fT22.fMultChSec = fT22.GetChMultiplicity(1);
	fT23.fMultChSec = fT23.GetChMultiplicity(1);

	return;
}

void BeEvent::SetSiDetEnergies(AculRaw *rawevent) {
	//set energies in MeV for each Si detector
	//rawevent: obtained from RAW data file
	//cal: object with calibration parameters for Si detectors
	//information in energy units

	Double_t energy;

	for (Int_t j = 0; j < 8; j++) {
		//T1.1 Sectors: C3[3][], C3[4][]
		energy = calSi->fA[3][2 * j + 1]
				* (rawevent->C3[3][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[3][2 * j + 1];
		if (energy > calSi->fC[3][2 * j + 1]) {
			fT11.fESec[2 * j + 1] = energy;
		}
		energy = calSi->fA[3][2 * j]
				* (rawevent->C3[3][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[3][2 * j];
		if (energy > calSi->fC[3][2 * j]) {
			fT11.fESec[2 * j] = energy;
		}
		energy = calSi->fA[4][2 * j + 1]
				* (rawevent->C3[4][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[4][2 * j + 1];
		if (energy > calSi->fC[4][2 * j + 1]) {
			fT11.fESec[2 * j + 1 + 16] = energy;
		}
		energy = calSi->fA[4][2 * j]
				* (rawevent->C3[4][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[4][2 * j];
		if (energy > calSi->fC[4][2 * j]) {
			fT11.fESec[2 * j + 16] = energy;
		}

//		//T1.2 Sectors: C3[13][], C3[14][]
		energy = calSi->fA[13][2 * j + 1]
				* (rawevent->C3[13][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[13][2 * j + 1];
		if (energy > calSi->fC[13][2 * j + 1]) {
			fT12.fESec[2 * j + 1] = energy;
		}
		energy = calSi->fA[13][2 * j]
				* (rawevent->C3[13][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[13][2 * j];
		if (energy > calSi->fC[13][2 * j]) {
			fT12.fESec[2 * j] = energy;
		}

		//T2.1 Sectors
		energy = calSi->fA[9][2 * j + 1]
				* (rawevent->C3[9][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[9][2 * j + 1];
		if (energy > calSi->fC[9][2 * j + 1]) {
			fT21.fESec[2 * j + 1] = energy;
		}
		energy = calSi->fA[9][2 * j]
				* (rawevent->C3[9][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[9][2 * j];
		if (energy > calSi->fC[9][2 * j]) {
			fT21.fESec[2 * j] = energy;
		}
		energy = calSi->fA[8][2 * j + 1]
				* (rawevent->C3[8][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[8][2 * j + 1];
		if (energy > calSi->fC[8][2 * j + 1]) {
			fT21.fESec[2 * j + 1 + 16] = energy;
		}
		energy = calSi->fA[8][2 * j]
				* (rawevent->C3[8][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[8][2 * j];
		if (energy > calSi->fC[8][2 * j]) {
			fT21.fESec[2 * j + 16] = energy;
		}

//		//T2.2 (Sectors)
		energy = calSi->fA[14][2 * j + 1]
				* (rawevent->C3[14][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[14][2 * j + 1];
		if (energy > calSi->fC[14][2 * j + 1]) {
			fT22.fESec[2 * j + 1] = energy;
		}
		energy = calSi->fA[14][2 * j]
				* (rawevent->C3[14][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[14][2 * j];
		if (energy > calSi->fC[14][2 * j]) {
			fT22.fESec[2 * j] = energy;
		}
	} //for

	for (Int_t j = 0; j < 8; j++) {
		//T1.1 Rings
		energy = calSi->fA[5][2 * j + 1]
				* (rawevent->C3[5][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[5][2 * j + 1];
		if (energy > calSi->fC[5][2 * j + 1]) {
			fT11.fERing[2 * j] = energy;
		}
		energy = calSi->fA[5][2 * j]
				* (rawevent->C3[5][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[5][2 * j];
		if (energy > calSi->fC[5][2 * j]) {
			fT11.fERing[2 * j + 1] = energy;
		}
		energy = calSi->fA[6][2 * j + 1]
				* (rawevent->C3[6][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[6][2 * j + 1];
		if (energy > calSi->fC[6][2 * j + 1]) {
			fT11.fERing[2 * j + 16] = energy;
		}
		energy = calSi->fA[6][2 * j]
				* (rawevent->C3[6][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[6][2 * j];
		if (energy > calSi->fC[6][2 * j]) {
			fT11.fERing[2 * j + 1 + 16] = energy;
		}

		energy = calSi->fA[10][2 * j + 1]
				* (rawevent->C3[10][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[10][2 * j + 1];
		if (energy > calSi->fC[10][2 * j + 1]) {
			fT21.fERing[2 * j] = energy;
		}
		energy = calSi->fA[10][2 * j]
				* (rawevent->C3[10][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[10][2 * j];
		if (energy > calSi->fC[10][2 * j]) {
			fT21.fERing[2 * j + 1] = energy;
		}
		energy = calSi->fA[11][2 * j + 1]
				* (rawevent->C3[11][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[11][2 * j + 1];
		if (energy > calSi->fC[11][2 * j + 1]) {
			fT21.fERing[2 * j + 16] = energy;
		}
		energy = calSi->fA[11][2 * j]
				* (rawevent->C3[11][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ calSi->fB[11][2 * j];
		if (energy > calSi->fC[11][2 * j]) {
			fT21.fERing[2 * j + 1 + 16] = energy;
		}
	} //for

	fT11.fMultESec = fT11.GetEMultiplicity(1);
	fT11.fMultERing = fT11.GetEMultiplicity(0);
	fT12.fMultESec = fT12.GetEMultiplicity(1);

	fT21.fMultESec = fT21.GetEMultiplicity(1);
	fT21.fMultERing = fT21.GetEMultiplicity(0);
	fT22.fMultESec = fT22.GetEMultiplicity(1);

	return;
}

void BeEvent::SetSiDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal) {
	//set energies in MeV for each Si detector
	//rawevent: obtained from RAW data file
	//cal: object with calibration parameters for Si detectors
	//information in energy units

	Double_t energy;

	for (Int_t j = 0; j < 8; j++) {
		//T1.1 Sectors: C3[3][], C3[4][]
		energy = cal->fA[3][2 * j + 1]
				* (rawevent->C3[3][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[3][2 * j + 1];
		if (energy > cal->fC[3][2 * j + 1]) {
			fT11.fESec[2 * j + 1] = energy;
		}
		energy = cal->fA[3][2 * j]
				* (rawevent->C3[3][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[3][2 * j];
		if (energy > cal->fC[3][2 * j]) {
			fT11.fESec[2 * j] = energy;
		}
		energy = cal->fA[4][2 * j + 1]
				* (rawevent->C3[4][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[4][2 * j + 1];
		if (energy > cal->fC[4][2 * j + 1]) {
			fT11.fESec[2 * j + 1 + 16] = energy;
		}
		energy = cal->fA[4][2 * j]
				* (rawevent->C3[4][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[4][2 * j];
		if (energy > cal->fC[4][2 * j]) {
			fT11.fESec[2 * j + 16] = energy;
		}

//		//T1.2 Sectors: C3[13][], C3[14][]
		energy = cal->fA[13][2 * j + 1]
				* (rawevent->C3[13][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[13][2 * j + 1];
		if (energy > cal->fC[13][2 * j + 1]) {
			fT12.fESec[2 * j + 1] = energy;
		}
		energy = cal->fA[13][2 * j]
				* (rawevent->C3[13][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[13][2 * j];
		if (energy > cal->fC[13][2 * j]) {
			fT12.fESec[2 * j] = energy;
		}

		//T2.1 Sectors
		energy = cal->fA[9][2 * j + 1]
				* (rawevent->C3[9][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[9][2 * j + 1];
		if (energy > cal->fC[9][2 * j + 1]) {
			fT21.fESec[2 * j + 1] = energy;
		}
		energy = cal->fA[9][2 * j]
				* (rawevent->C3[9][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[9][2 * j];
		if (energy > cal->fC[9][2 * j]) {
			fT21.fESec[2 * j] = energy;
		}
		energy = cal->fA[8][2 * j + 1]
				* (rawevent->C3[8][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[8][2 * j + 1];
		if (energy > cal->fC[8][2 * j + 1]) {
			fT21.fESec[2 * j + 1 + 16] = energy;
		}
		energy = cal->fA[8][2 * j]
				* (rawevent->C3[8][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[8][2 * j];
		if (energy > cal->fC[8][2 * j]) {
			fT21.fESec[2 * j + 16] = energy;
		}

//		//T2.2 (Sectors)
		energy = cal->fA[14][2 * j + 1]
				* (rawevent->C3[14][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[14][2 * j + 1];
		if (energy > cal->fC[14][2 * j + 1]) {
			fT22.fESec[2 * j + 1] = energy;
		}
		energy = cal->fA[14][2 * j]
				* (rawevent->C3[14][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[14][2 * j];
		if (energy > cal->fC[14][2 * j]) {
			fT22.fESec[2 * j] = energy;
		}
	} //for

	for (Int_t j = 0; j < 8; j++) {
		//T1.1 Rings
		energy = cal->fA[5][2 * j + 1]
				* (rawevent->C3[5][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[5][2 * j + 1];
		if (energy > cal->fC[5][2 * j + 1]) {
			fT11.fERing[2 * j] = energy;
		}
		energy = cal->fA[5][2 * j]
				* (rawevent->C3[5][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[5][2 * j];
		if (energy > cal->fC[5][2 * j]) {
			fT11.fERing[2 * j + 1] = energy;
		}
		energy = cal->fA[6][2 * j + 1]
				* (rawevent->C3[6][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[6][2 * j + 1];
		if (energy > cal->fC[6][2 * j + 1]) {
			fT11.fERing[2 * j + 16] = energy;
		}
		energy = cal->fA[6][2 * j]
				* (rawevent->C3[6][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[6][2 * j];
		if (energy > cal->fC[6][2 * j]) {
			fT11.fERing[2 * j + 1 + 16] = energy;
		}

		energy = cal->fA[10][2 * j + 1]
				* (rawevent->C3[10][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[10][2 * j + 1];
		if (energy > cal->fC[10][2 * j + 1]) {
			fT21.fERing[2 * j] = energy;
		}
		energy = cal->fA[10][2 * j]
				* (rawevent->C3[10][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[10][2 * j];
		if (energy > cal->fC[10][2 * j]) {
			fT21.fERing[2 * j + 1] = energy;
		}
		energy = cal->fA[11][2 * j + 1]
				* (rawevent->C3[11][2 * j + 1] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[11][2 * j + 1];
		if (energy > cal->fC[11][2 * j + 1]) {
			fT21.fERing[2 * j + 16] = energy;
		}
		energy = cal->fA[11][2 * j]
				* (rawevent->C3[11][2 * j] + ranE.Uniform(-0.5, 0.5))
				+ cal->fB[11][2 * j];
		if (energy > cal->fC[11][2 * j]) {
			fT21.fERing[2 * j + 1 + 16] = energy;
		}
	} //for

	fT11.fMultESec = fT11.GetEMultiplicity(1);
	fT11.fMultERing = fT11.GetEMultiplicity(0);
	fT12.fMultESec = fT12.GetEMultiplicity(1);

	fT21.fMultESec = fT21.GetEMultiplicity(1);
	fT21.fMultERing = fT21.GetEMultiplicity(0);
	fT22.fMultESec = fT22.GetEMultiplicity(1);

	return;
}

void BeEvent::SetSiDetTimes(AculRaw *rawevent) {
	//set times in ns for each Si detector
	//rawevent: obtained from RAW data file
	//cal: object with calibration parameters for Si detectors

	//information in ns

	const Double_t timeshift = 100.;

	for (Int_t j = 0; j < 16; j++) {
		//T1.2 (Sectors)
		fTau[0][j] = calSi->fA[20][j]
				* (rawevent->C3[20][j] + ranE.Uniform(-0.5, 0.5))
				- calSi->fB[20][j] + timeshift;
		//T2.2 (Sectors)
		fTau[1][j] = calSi->fA[22][j]
				* (rawevent->C3[22][j] + ranE.Uniform(-0.5, 0.5))
				- calSi->fB[22][j] + timeshift;
	} //for

	return;
}

void BeEvent::SetSiDetTimesOld(AculRaw *rawevent, AculCalibration *cal) {
	//set times in ns for each Si detector
	//rawevent: obtained from RAW data file
	//cal: object with calibration parameters for Si detectors

	//information in ns

	const Double_t timeshift = 100.;

	for (Int_t j = 0; j < 16; j++) {
		//T1.2 (Sectors)
		fTau[0][j] = cal->fA[20][j]
				* (rawevent->C3[20][j] + ranE.Uniform(-0.5, 0.5))
				- cal->fB[20][j] + timeshift;
		//T2.2 (Sectors)
		fTau[1][j] = cal->fA[22][j]
				* (rawevent->C3[22][j] + ranE.Uniform(-0.5, 0.5))
				- cal->fB[22][j] + timeshift;
	} //for

	return;
}

Double_t BeEvent::CsIEnergy(Int_t _address, Int_t _subaddress,
		AculRaw *_rawevent, AculCalibration *_cal, Double_t _x0,
		Double_t _ethr) {
	//CsI in higher channel: f2,  pol1 ax+b
	//CsI in low channel: f1,  c1*(x-p)^c2, p == piedestal position
	//
	//equations: f1(x0) = f2(x0) && f1'(x0) = f2'(x0), where fi' = df/dx
	//
	//c1*(x0-p)^c2 = a*x0+b  &&  c1*c2*(x0-p)^(c2-1) = a
	//
	//	c2 = a*(x0-p)/(a*x0+b)
	//
	//	c1 = (a*x0+b)/( (x0-p)^c2 )

	Double_t energy = 0;
	Double_t cp1 = 0;
	Double_t cp2 = 0;

	if (_rawevent->C3[_address][_subaddress]
			> _cal->fC[_address][_subaddress]) {
		if (_rawevent->C3[_address][_subaddress] < _x0) {
			cp2 = _cal->fA[_address][_subaddress]
					* (_x0 - _cal->fC[_address][_subaddress])
					/ (_cal->fA[_address][_subaddress] * _x0
							+ _cal->fB[_address][_subaddress]);
			cp1 =
					(_cal->fA[_address][_subaddress] * _x0
							+ _cal->fB[_address][_subaddress])
							/ (TMath::Power(
									_x0 - _cal->fC[_address][_subaddress], cp2));
			energy = cp1
					* TMath::Power(
							_rawevent->C3[_address][_subaddress]
									- _cal->fC[_address][_subaddress]
									+ ranE.Uniform(-0.5, 0.5), cp2);
		} else {
			energy = _cal->fA[_address][_subaddress]
					* (_rawevent->C3[_address][_subaddress]
							+ ranE.Uniform(-0.5, 0.5))
					+ _cal->fB[_address][_subaddress];
		}
	} //if
	if (energy < _ethr) {
		return 0;
	}

	return energy;
}

Double_t BeEvent::CsIcorrectedEnergy(Int_t _address, Int_t _subaddress,
		AculRaw *_rawevent, AculCalibration *_cal, Double_t _x0,
		Double_t _ethr) {
	//using this function, corrections of parameter "a" could be checked, corrections are written in DA logbook
	//
	//CsI in higher channel: f2,  pol1 ax+b
	//CsI in low channel: f1,  c1*(x-p)^c2, p == piedestal position
	//
	//equations: f1(x0) = f2(x0) && f1'(x0) = f2'(x0), where fi' = df/dx
	//
	//c1*(x0-p)^c2 = a*x0+b  &&  c1*c2*(x0-p)^(c2-1) = a
	//
	//	c2 = a*(x0-p)/(a*x0+b)
	//
	//	c1 = (a*x0+b)/( (x0-p)^c2 )

//	Double_t energy;
	Double_t energy = 0.;
	Double_t cp1;
	Double_t cp2;

	if (_rawevent->C3[_address][_subaddress]
			> _cal->fC[_address][_subaddress]) {
		if (_rawevent->C3[_address][_subaddress] < _x0) {
			cp2 = fCalCorr[_subaddress] * _cal->fA[_address][_subaddress]
					* (_x0 - _cal->fC[_address][_subaddress])
					/ (_cal->fA[_address][_subaddress] * _x0
							+ _cal->fB[_address][_subaddress]);
			cp1 =
					(fCalCorr[_subaddress] * _cal->fA[_address][_subaddress]
							* _x0 + _cal->fB[_address][_subaddress])
							/ (TMath::Power(
									_x0 - _cal->fC[_address][_subaddress], cp2));
			energy = cp1
					* TMath::Power(
							_rawevent->C3[_address][_subaddress]
									- _cal->fC[_address][_subaddress], cp2);
		} else {
			energy = fCalCorr[_subaddress] * _cal->fA[_address][_subaddress]
					* (_rawevent->C3[_address][_subaddress]
							+ ranE.Uniform(-0.5, 0.5))
					+ _cal->fB[_address][_subaddress];
		}
	} //if
	if (energy < _ethr) {
		return 0;
	}

	return energy;
}

void BeEvent::SetCsIalphaDetEnergies(AculRaw *rawevent) {
	//set energies of alpha particles in MeV for each CsI detector
	//rawevent: obtained from RAW data file
	//cal: object with calibration parameters for alpha particles in CsI detectors

	Double_t x0 = 500; //connection point [channels]

	//T13
	for (Int_t j = 0; j < 8; j++) {
//		fE13aSec[j] = CsIEnergy(15, j, rawevent, cal, x0, calCsIa->fD[15][j]);
		fE13aSec[j] = CsIEnergy(15, j, rawevent, calCsIa, x0,
				calCsIa->fD[15][j])*fCorrT1_0_7_CsI_A;
	} //for
	for (Int_t j = 8; j < 16; j++) {
//		fE13aSec[j] = CsIEnergy(15, j, rawevent, calCsIa, x0, calCsIa->fD[15][j]);
		fE13aSec[j] = CsIEnergy(15, j, rawevent, calCsIa, x0,
				calCsIa->fD[15][j])*fCorrT1_8_15_CsI_A;
	} //for

	//T23
	for (Int_t j = 0; j < 8; j++) {
//		fE23aSec[j] = CsIEnergy(16, j, rawevent, calCsIa, x0, calCsIa->fD[16][j]);
		fE23aSec[j] = CsIEnergy(16, j, rawevent, calCsIa, x0,
				calCsIa->fD[16][j])*fCorrT2_0_7_CsI_A;
	}

	for (Int_t j = 8; j < 16; j++) {
//		fE23aSec[j] = CsIEnergy(16, j, rawevent, calCsIa, x0, calCsIa->fD[16][j]);
		fE23aSec[j] = CsIEnergy(16, j, rawevent, calCsIa, x0,
				calCsIa->fD[16][j])*fCorrT2_8_15_CsI_A;
	}

	return;
}

void BeEvent::SetCsIalphaDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal) {
	//set energies of alpha particles in MeV for each CsI detector
	//rawevent: obtained from RAW data file
	//cal: object with calibration parameters for alpha particles in CsI detectors

	Double_t x0 = 500; //connection point [channels]
//	const Double_t ethreshold = 0.5;

//	const Double_t calcorrection = 1.1;	//0.90; 0.95; 1.05; 1.1

	//T13
	for (Int_t j = 0; j < 8; j++) {
//		fE13aSec[j] = CsIEnergy(15, j, rawevent, cal, x0, cal->fD[15][j]);
		fE13aSec[j] = CsIEnergy(15, j, rawevent, cal, x0,
				cal->fD[15][j])*fCorrT1_0_7_CsI_A;

//		fE13aSec[j] = CsIcorrectedEnergy(15, j, rawevent, cal, x0, ETHR);
	} //for
	for (Int_t j = 8; j < 16; j++) {
//		fE13aSec[j] = CsIEnergy(15, j, rawevent, cal, x0, cal->fD[15][j]);
		fE13aSec[j] = CsIEnergy(15, j, rawevent, cal, x0,
				cal->fD[15][j])*fCorrT1_8_15_CsI_A;

//		fE13aSec[j] = CsIEnergy(15, j, rawevent, cal, x0, ethreshold);
//		fE13aSec[j] = CsIcorrectedEnergy(15, j, rawevent, cal, x0, ETHR);
	} //for

	//T23
	for (Int_t j = 0; j < 8; j++) {
//		fE23aSec[j] = CsIEnergy(16, j, rawevent, cal, x0, cal->fD[16][j]);
		fE23aSec[j] = CsIEnergy(16, j, rawevent, cal, x0,
				cal->fD[16][j])*fCorrT2_0_7_CsI_A;

//		fE23aSec[j] = CsIEnergy(16, j, rawevent, cal, x0, ethreshold);
	}

	for (Int_t j = 8; j < 16; j++) {
//		fE23aSec[j] = CsIEnergy(16, j, rawevent, cal, x0, cal->fD[16][j]);
		fE23aSec[j] = CsIEnergy(16, j, rawevent, cal, x0,
				cal->fD[16][j])*fCorrT2_8_15_CsI_A;

//		fE23aSec[j] = CsIEnergy(16, j, rawevent, cal, x0, ethreshold);
	}

	return;
}

void BeEvent::SetCsIprotonDetEnergies(AculRaw *rawevent) {
	//set energies of protons in MeV for each CsI detector
	//rawevent: obtained from RAW data file
	//cal: object with calibration parameters for protons in CsI detectors

	Double_t x0 = 500; //connection point [channels]

	//T13
	for (Int_t j = 0; j < 8; j++) {
//		fE13pSec[j] = CsIEnergy(15, j, rawevent, calCsIp, x0, calCsIp->fD[15][j]);
		fE13pSec[j] = CsIEnergy(15, j, rawevent, calCsIp, x0,
				calCsIp->fD[15][j])*fCorrT1_0_7_CsI_P;

//		fE13pSec[j] = CsIcorrectedEnergy(15, j, rawevent, calCsIp, x0, ETHR);
	}
	for (Int_t j = 8; j < 16; j++) {
//		fE13pSec[j] = CsIEnergy(15, j, rawevent, calCsIp, x0, calCsIp->fD[15][j]);
		fE13pSec[j] = CsIEnergy(15, j, rawevent, calCsIp, x0,
				calCsIp->fD[15][j])*fCorrT1_8_15_CsI_P;

//		fE13pSec[j] = CsIcorrectedEnergy(15, j, rawevent, calCsIp, x0, ETHR);
	}

	//T23
	for (Int_t j = 0; j < 8; j++) {
//		fE23pSec[j] = CsIEnergy(16, j, rawevent, calCsIp, x0, calCsIp->fD[16][j]);
		fE23pSec[j] = CsIEnergy(16, j, rawevent, calCsIp, x0,
				calCsIp->fD[16][j])*fCorrT2_0_7_CsI_P;

//		fE23pSec[j] = CsIEnergy(16, j, rawevent, calCsIp, x0, ethreshold);
//		fE23pSec[j] = CsIcorrectedEnergy(16, j, rawevent, calCsIp, x0, ETHR);	//zkontrolovat, co je to zac

	}

	for (Int_t j = 8; j < 16; j++) {
//		fE23pSec[j] = CsIEnergy(16, j, rawevent, calCsIp, x0, calCsIp->fD[16][j]);
		fE23pSec[j] = CsIEnergy(16, j, rawevent, calCsIp, x0,
				calCsIp->fD[16][j])*fCorrT2_8_15_CsI_P;

//		fE23pSec[j] = CsIEnergy(16, j, rawevent, calCsIp, x0, ethreshold);
//		fE23pSec[j] = CsIcorrectedEnergy(16, j, rawevent, calCsIp, x0, ETHR);	//zkontrolovat, co je to zac

	}

	fMult13p = fT13.GetEMultiplicity(1);
	fMult23p = fT23.GetEMultiplicity(1);

	return;
}

void BeEvent::SetCsIprotonDetEnergiesOld(AculRaw *rawevent, AculCalibration *cal) {
	//set energies of protons in MeV for each CsI detector
	//rawevent: obtained from RAW data file
	//cal: object with calibration parameters for protons in CsI detectors

	Double_t x0 = 500; //connection point [channels]

	//T13
	for (Int_t j = 0; j < 8; j++) {
//		fE13pSec[j] = CsIEnergy(15, j, rawevent, cal, x0, cal->fD[15][j]);
		fE13pSec[j] = CsIEnergy(15, j, rawevent, cal, x0,
				cal->fD[15][j])*fCorrT1_0_7_CsI_P;

//		fE13pSec[j] = CsIcorrectedEnergy(15, j, rawevent, cal, x0, ETHR);
	}
	for (Int_t j = 8; j < 16; j++) {
//		fE13pSec[j] = CsIEnergy(15, j, rawevent, cal, x0, cal->fD[15][j]);
		fE13pSec[j] = CsIEnergy(15, j, rawevent, cal, x0,
				cal->fD[15][j])*fCorrT1_8_15_CsI_P;

//		fE13pSec[j] = CsIcorrectedEnergy(15, j, rawevent, cal, x0, ETHR);
	}

	//T23
	for (Int_t j = 0; j < 8; j++) {
//		fE23pSec[j] = CsIEnergy(16, j, rawevent, cal, x0, cal->fD[16][j]);
		fE23pSec[j] = CsIEnergy(16, j, rawevent, cal, x0,
				cal->fD[16][j])*fCorrT2_0_7_CsI_P;

//		fE23pSec[j] = CsIEnergy(16, j, rawevent, cal, x0, ethreshold);
//		fE23pSec[j] = CsIcorrectedEnergy(16, j, rawevent, cal, x0, ETHR);	//zkontrolovat, co je to zac

	}

	for (Int_t j = 8; j < 16; j++) {
//		fE23pSec[j] = CsIEnergy(16, j, rawevent, cal, x0, cal->fD[16][j]);
		fE23pSec[j] = CsIEnergy(16, j, rawevent, cal, x0,
				cal->fD[16][j])*fCorrT2_8_15_CsI_P;

//		fE23pSec[j] = CsIEnergy(16, j, rawevent, cal, x0, ethreshold);
//		fE23pSec[j] = CsIcorrectedEnergy(16, j, rawevent, cal, x0, ETHR);	//zkontrolovat, co je to zac

	}

	fMult13p = fT13.GetEMultiplicity(1);
	fMult23p = fT23.GetEMultiplicity(1);

	return;
}

Double_t BeEvent::ReactionQ(Double_t Ta1, Double_t Ta2, Double_t mac2,
		Double_t mbc2, Double_t theta) {
	//	elastic scattering;
	//
	//////////////////////////////////
	//								//
	//				   a2			//
	//				  /				//
	//			     /				//
	//	            /alpha			//
	//	a1-------->X----------		//
	//            b1\				//
	//			     \				//
	//			      b2			//
	//								//
	//								//
	//Q = (Ea1 + Eb1)-(Ea2-Eb2)
	//Q = (Ta1 + m_a*c^2 + Tb1 + m_b*c^2) - (Ta2 + m_a*c^2 + Tb2 + m_b*c^2)
	//known: masses, Ta1, Tb1==0, Ta2
	//unknown: Tb2
	//
	//pa1 + pb1 = pa2 + pb2, where (p*c)^2 = T^2 + 2*T*(m*c^2) and pb1==0
	//pb2^2 = pa1^2 + pa2^2 - 2*pa1*pa2*cos(theta)
	//
	//	Eb^2 = (Tb2 + m_b*c^2)^2 = (pb2c*c)^2 + (m_b*c^2)^2

	Double_t pa1c2 = PC2(Ta1, mac2);
	Double_t pa2c2 = PC2(Ta2, mac2);

	Double_t p_b2c2 = LawOfCosines(TMath::Power(pa1c2, 0.5),
			TMath::Power(pa2c2, 0.5), theta);

	return Ta1 + mbc2 - Ta2 - TMath::Power(p_b2c2 + TMath::Power(mbc2, 2), 0.5);
}

Double_t BeEvent::ReactionQ(Double_t Ea, Double_t Eb, Double_t Ec,
		Double_t Ed) {
	//	binary reaction;
	//
	//////////////////////////////////
	//								//
	//				   c			//
	//				  /				//
	//			     /				//
	//	            /alpha			//
	//	 a-------->X----------		//
	//             b\				//
	//			     \				//
	//			      d				//
	//								//
	//								//
	//Q = (Ea + Eb)-(Ec-Ed)
	//Q = (Ta + m_a*c^2 + Tb + m_b*c^2) - (Tc + m_c*c^2 + Td + m_d*c^2)
	//known: masses, Ta, Tb==0, Tc, Td

	return (Ea + Eb) - (Ec + Ed);
}

Double_t BeEvent::PC2(Double_t T, Double_t mc2) {
	//pc^2 = T^2 + 2*T*mc^2

	return TMath::Power(T, 2) + 2 * T * mc2;
}

Double_t BeEvent::T(Double_t pc2, Double_t mc2) {
	return TMath::Power(mc2 * mc2 + pc2, 0.5) - mc2;
}

Double_t BeEvent::LawOfCosines(Double_t pb, Double_t pc, Double_t alpha) {
	//pa^2 = pb^2 + pc^2 - 2*pb*pc*cos(alpha)

	return TMath::Power(pb, 2) + TMath::Power(pc, 2)
			- 2 * pb * pc * TMath::Cos(alpha);
}

void BeEvent::Set6LiQ() {
	if ((fT11.fMultERing != 1)) { //todo pochybne: zkontrolovat, co je pochybne
		return;
	}

	Int_t rn = fT11.FindRing();

	fSN1[0] = fT11.FindSec() / 2;
	fTheta1 = fT11.GetTheta(rn);

	fT6Li2 = fdE1[fSN1[0]] + fT12.fESec[fSN1[0]] + fE13aSec[fSN1[0]];

	fQ1 = ReactionQ(6. * 36., fT6Li2, 6000, 1000, fTheta1);

	return;
}

void BeEvent::SetDeltaE() {
	//fill variables fdE1 and fdE2, where Double_t fdEi[16]

	//deltaE1
	for (Int_t i = 0; i < NOESEC; i++) {
		if (fT11.fESec[2 * i] == 0 || fT11.fESec[2 * i + 1] == 0) {
			fdE1[i] = fT11.fESec[2 * i] + fT11.fESec[2 * i + 1];
		}
		if (fT21.fESec[2 * i] == 0 || fT21.fESec[2 * i + 1] == 0) {
			fdE2[i] = fT21.fESec[2 * i] + fT21.fESec[2 * i + 1];
		}
	}

	return;
}

void BeEvent::Reset() {

	//trigger
	fTrigger = 0;

	//detectors
	fT11.Reset();
	fT12.Reset();
	fT13.Reset();
	fT21.Reset();
	fT22.Reset();
	fT23.Reset();

	//6Be information
	fNOA = 0;
	fNOP = 0;

	f6BeIM = -100.;
	//lab
	f6BeLab.SetPxPyPzE(0., 0., 0., 0.);
	f6BeThetaLab = 0.;
	f6BePhiLab = 0;
	f6BePcLab = 0;
	fNLab.SetPxPyPzE(0., 0., 0., 0.);
	fNThetaLab = 0;
	fNPhiLab = 0;
	fNPcLab = 0;
	fP1Lab.SetPxPyPzE(0., 0., 0., 0.);
	fP1ThetaLab = 0;
	fP1PhiLab = 0;
	fP1PcLab = 0;
	fP2Lab.SetPxPyPzE(0., 0., 0., 0.);
	fP2ThetaLab = 0;
	fP2PhiLab = 0;
	fP2PcLab = 0;
	fALab.SetPxPyPzE(0., 0., 0., 0.);
	fAThetaLab = 0;
	fAPhiLab = 0;
	fAPcLab = 0;
	//CM1
	f6BeCM1.SetPxPyPzE(0., 0., 0., 0.);
	f6BeThetaCM1 = 0.;
	f6BePhiCM1 = 0;
	f6BePcCM1 = 0;
	//CM
	f6BeCM.SetPxPyPzE(0., 0., 0., 0.);
	f6BeThetaCM = 0.;
	f6BePhiCM = 0;
	f6BePcCM = 0;
	fACM.SetPxPyPzE(0., 0., 0., 0.);
	fAThetaCM = 0.;
	fAPhiCM = 0;
	fAPcCM = 0;
	fP1CM.SetPxPyPzE(0., 0., 0., 0.);
	fP1ThetaCM = 0.;
	fP1PhiCM = 0;
	fP1PcCM = 0;
	fP2CM.SetPxPyPzE(0., 0., 0., 0.);
	fP2ThetaCM = 0.;
	fP2PhiCM = 0;
	fP2PcCM = 0;
	fNCM.SetPxPyPzE(0., 0., 0., 0.);
	fNThetaCM = 0.;
	fNPhiCM = 0;
	fNPcCM = 0;

	fQLiP = kLiMASS;
	fTpp = 0.;
	fTapp = 0.;
	fCosThetaTk = 2.;
	fTap = 0.;
	fTpap = 0.;
	fCosThetaYk = 2.;

	//hits
	for (Int_t i = 0; i < NOTEL; i++) {

		fNOHits[i] = 0;
		fNOSHits[i] = 0;
		fNORHits[i] = 0;
		fNOaHits[i] = 0;
		fNOpHits[i] = 0;
		fNOtHits[i] = 0;
//		cout << "\t" << fNOtHits[i] << endl;

		for (Int_t j = 0; j < NOHITS; j++) {
			fSN1[j] = -1;
			fSN[i][j] = -1;
			fTauSN[i][j] = -1;
			fRN[i][j] = -1;
			fE1Sec[i][j] = 0.;
			fE1Ring[i][j] = 0.;
			fE2Sec[i][j] = 0.;
			fE3Sec[i][j] = 0.;
			fTheta[i][j] = -1.;
			fPhi[i][j] = -1.;
			fhTau[i][j] = 0;
			fT[i][j] = 0.;
			fP[i][j] = 0.;
			fdE1calc[i][j] = 0.;
			fID[i][j] = 0;
			fPM[j] = 0;
			fPM[NOHITS + j] = 0;
		}

	}

	fTheta1 = 0.;
	fTheta2 = 0.;

	fX11 = 0.;
	fX12 = 0.;
	fXd112 = 0.;
	fXd123 = 0.;

	fX21 = 0.;
	fX22 = 0.;
	fXd212 = 0.;
	fXd223 = 0.;

	fMult13a = 0;
	fMult23a = 0;
	fMult13p = 0;
	fMult23p = 0;

	for (Int_t i = 0; i < NOESEC; i++) {
		fn[i] = i;

		fdE1[i] = 0.;
		fdE2[i] = 0.;

		fE11aSec[i] = 0.;
		fE11pSec[i] = 0.;
		fE12aSec[i] = 0.;
		fE12pSec[i] = 0.;
		fE13aSec[i] = 0.;
		fE13pSec[i] = 0.;

		fE21aSec[i] = 0.;
		fE21pSec[i] = 0.;
		fE22aSec[i] = 0.;
		fE22pSec[i] = 0.;
		fE23aSec[i] = 0.;
		fE23pSec[i] = 0.;

		//1st telescope verification
		fE11v_0s[i] = 0;
		fE11vs[i] = 0;
		fE12v_0[i] = 0;
		fE12v[i] = 0;

		//1st telescope verification
		fE21v_0s[i] = 0;
		fE21vs[i] = 0;
		fE22v_0[i] = 0;
		fE22v[i] = 0;

		//alpha
		fE12_0a[i] = 0;
		fE12a[i] = 0;
		fE13_0a[i] = 0;

		//proton
		fE12_0p[i] = 0;
		fE12p[i] = 0;
		fE13_0p[i] = 0;

		//alpha
		fE22_0a[i] = 0;
		fE22a[i] = 0;
		fE23_0a[i] = 0;

		//proton
		fE22_0p[i] = 0;
		fE22p[i] = 0;
		fE23_0p[i] = 0;

	}

	//for BeEvent::DeltaECalibration()

	fE11v_0r = 0;
	fE11vr = 0;

	fE21v_0r = 0;
	fE21vr = 0;

	fRN2 = -1;
	fSN2 = -1;

	return;
}

void BeEvent::SetTheta1(Int_t ring) {
	fTheta1 = fT11.GetTheta(ring);
	return;
}

void BeEvent::SetTheta2(Int_t ring) {
	fTheta2 = fT21.GetTheta(ring);
	return;
}

/*
void BeEvent::DeltaEIdentification1() {

	//vypocet provedu pouze pro udalosti s nasobnosti ringu == 1
	if ((fT11.fMultERing != 1)) {
		return;
	}

	Double_t Ep = 0.;
	Double_t Ea = 0.;

	Int_t rn = fT11.FindRing(); //cislo ringu

	//function is probably not used, nevertheless:
	//potreba zavest mrtve vrstvy
	//pro kazdy sektor musim najit odpovidajici ring v prvnim detektoru a podle neho spocitat Theta1

	if (rn >= 2) {
		fRN1 = rn;
		fSN1[0] = fT11.FindSec() / 2;

		fTheta1 = fT11.GetTheta(rn);
		fX11 = X11 / TMath::Cos(fTheta1);
		fXd112 = (XD11B + XD12F) / TMath::Cos(fTheta1);
		fX12 = X12 / TMath::Cos(fTheta1);
		fXd123 = XD13F / TMath::Cos(fTheta1);

		//alphas
		if (fE13aSec[fSN1[0]] > 0.) {
			Ea = mSiAlpha.GetE0(fE13aSec[fSN1[0]], fXd123);
			fE12aSec[fSN1[0]] = mSiAlpha.GetE0(Ea, fX12) - Ea;
//			Ea = mSiAlpha.GetE0(Ea, fX22+fXd212);
		} else {
			Ea = mSiAlpha.GetE0(fT12.fESec[fSN1[0]], fXd112);
			fE11aSec[fSN1[0]] = mSiAlpha.GetE0(Ea, fX11) - Ea;
		}

		//protons
		if (fE13pSec[fSN1[0]] > 0.) {
			Ep = mSiP.GetE0(fE13pSec[fSN1[0]], fXd123);
			fE12pSec[fSN1[0]] = mSiP.GetE0(Ep, fX12) - Ep;
//			Ep = mSiP.GetE0(Ep, fX22+fXd212);
		} else {
			Ep = mSiP.GetE0(fT12.fESec[fSN1[0]], fXd112);
			fE11pSec[fSN1[0]] = mSiP.GetE0(Ep, fX11) - Ep;
		}

	} //if

	return;
}
*/

/*
void BeEvent::DeltaEIdentification2() {
	//todo tato funkce se pravdepodobne nikde nepouziva
	//nicmene, byly tady urcite nedodelky:
	//a) potreba zavest mrtve vrstvy
	//		pro kazdy sektor musim najit odpovidajici ring v prvnim detektoru a podle neho spocitat Theta1
	//b) energy loses in target material

	//vypocet provedu pouze pro udalosti s nasobnosti ringu == 1
	if ((fT21.fMultERing != 1)) {
		return;
	}

	//vsechny rozmery detektoru zapsat do jejich parametru
//	const Double_t	x0_11 = 275.;		//tloustka T11 v mcm, pouzivat 275
//	const Double_t	x0_12 = 800.;		//tloustka T12 v mcm, pouzivat 800
//	const Double_t	xd12 = 1.38 + 2.23;	//mrtva vrstva, pouzivat 1.7 + 2.36
//	const Double_t	xd23 = 14.;

	Double_t Ea = 0.;
	Double_t Ep = 0.;

	Int_t rn = fT21.FindRing(); //cislo ringu

	if (rn >= 2) {
		fRN2 = rn;
		fSN2 = fT21.FindSec() / 2;

		fTheta2 = fT21.GetTheta(rn);
		fX21 = X21 / TMath::Cos(fTheta2);
		fXd212 = (XD21B + XD22F) / TMath::Cos(fTheta2);
		fX22 = X22 / TMath::Cos(fTheta2);
		fXd223 = XD23F / TMath::Cos(fTheta2);

		//alphas
		if (fE23aSec[fSN2] > 0.) {
			Ea = mSiAlpha.GetE0(fE23aSec[fSN2], fXd223);
			fE22aSec[fSN2] = mSiAlpha.GetE0(Ea, fX22) - Ea;
//			Ea = mSiAlpha.GetE0(Ea, fX22+fXd212);
		} else {
			Ea = mSiAlpha.GetE0(fT22.fESec[fSN2], fXd212);
			fE21aSec[fSN2] = mSiAlpha.GetE0(Ea, fX21) - Ea;
		}

		//protons
		if (fE23pSec[fSN2] > 0.) {
			Ep = mSiP.GetE0(fE23pSec[fSN2], fXd223);
			fE22pSec[fSN2] = mSiP.GetE0(Ep, fX22) - Ep;
//			Ep = mSiP.GetE0(Ep, fX22+fXd212);
		} else {
			Ep = mSiP.GetE0(fT22.fESec[fSN2], fXd212);
			fE21pSec[fSN2] = mSiP.GetE0(Ep, fX21) - Ep;
		}

	} //if

	return;
}
*/

void BeEvent::DeltaECalibration1() {
	//kalibraci provadim za podminek:
	//1) mult v T11 == 1
	//2) pritomnost nenuloveho signalu v sektoru lezicim za zkoumanym sektorem

	if ((fT11.fMultERing != 1) && (fT11.fMultESec != 1)) {
		return;
	}

	Warning("BeEvent::DeltaECalibration1",
			"Detector thicknesses set incorrectly!!!");
	//should be taken from constants
	const Double_t x0_11 = 280.; //tloustka T11 v mcm, pouzivat 280
	const Double_t x0_12 = 850.; //tloustka T12 v mcm, pouzivat 900
	const Double_t xd112 = 1.7 + 2.36; //mrtva vrstva, pouzivat 1.7 + 2.36
	const Double_t xd123 = 14.; //mrtva vrstva, pouzivat XXXX

	Int_t rn = fT11.FindRing();
	if (rn >= 2) {
		//angles and sectors and rings numbers
		SetTheta1(rn);
		fRN1 = rn;
		fSN1[0] = fT11.FindSec() / 2;

		//detectors thicknessess
		//vypocitam drahu l v kremiku prvniho detektoru: l = det.d/cos(alfa)
		//	Double_t l = fT11.fThickness/TMath::Cos(fTheta1);
		fX11 = x0_11 / TMath::Cos(fTheta1);
		fX12 = x0_12 / TMath::Cos(fTheta1);
		fXd112 = xd112 / TMath::Cos(fTheta1);
		fXd123 = xd123 / TMath::Cos(fTheta1);

		///////////////
		////T11:T12////
		///////////////

		//verification of the method

		//musim udelat objekt kremikoveho materialu, dost mozna jako clen objektu annulardetector
		//z mereni znam deltaE, vypocitam E0 pred dopadem: E0 = GetE(deltaE)
		fE11v_0r = mSiAlphaV.GetE0dE(fT11.fERing[fRN1]);
		fE11v_0s[fSN1[0]] = mSiAlphaV.GetE0dE(fdE1[fSN1[0]]);
//		if ( (fE11v_0r < 0) || (fE11v_0s[fSN] < 0) ) { return; }

		// E po vystupu z prvni vrstvy kremiku (citliva oblast + mrtva vrstva) bude E = Get( E0 - deltaE, xdead/cos(alfa) )
		//fE12_0 = fE11_0s - fT11.fESec[fSN];	//prvni priblizeni
		fE11vs[fSN1[0]] = fE11v_0s[fSN1[0]]
				- mSiAlphaV.GetE(fE11v_0s[fSN1[0]], fX11);
		fE11vr = fE11v_0r - mSiAlphaV.GetE(fE11v_0r, fX11);

		fE12v_0[fSN1[0]] = mSiAlphaV.GetE(fE11v_0s[fSN1[0]] - fE11vs[fSN1[0]],
				fXd112);

		//E muzu zapsat jako vypocitanou energii v druhem kremiku
		fE12v[fSN1[0]] = fE12v_0[fSN1[0]]
				- mSiAlphaV.GetE(fE12v_0[fSN1[0]], fX12);

		///////////////
		////T12:T13////
		///////////////

		//alha
		fE12_0a[fSN1[0]] = mSiAlpha.GetE0dE(fT12.fESec[fSN1[0]]);
		fE13_0a[fSN1[0]] = mSiAlpha.GetE(fE12_0a[fSN1[0]], fX12);
		fE12a[fSN1[0]] = fE12_0a[fSN1[0]] - fE13_0a[fSN1[0]];
		fE13_0a[fSN1[0]] = mSiAlpha.GetE(fE13_0a[fSN1[0]], fXd123);

		//proton
		fE12_0p[fSN1[0]] = mSiP.GetE0dE(fT12.fESec[fSN1[0]]);
		fE13_0p[fSN1[0]] = mSiP.GetE(fE12_0p[fSN1[0]], fX12);
		fE12p[fSN1[0]] = fE12_0p[fSN1[0]] - fE13_0p[fSN1[0]];
		fE13_0p[fSN1[0]] = mSiP.GetE(fE13_0p[fSN1[0]], fXd123); //after the dead layer

	} //if

	return;

}

void BeEvent::DeltaECalibration2() {
	//kalibraci provadim za podminek:
	//1) mult v T11 == 1
	//2) pritomnost nenuloveho signalu v sektoru lezicim za zkoumanym sektorem

	if ((fT21.fMultERing != 1) && (fT21.fMultESec != 1)) {
		return;
	}

	Warning("BeEvent::DeltaECalibration2",
			"Detector thicknesses set incorrectly!!!");
	//should be taken from constants
	const Double_t x0_21 = 275.; //tloustka T21 v mcm, pouzivat 275
	const Double_t x0_22 = 800.; //tloustka T22 v mcm, pouzivat 800
	const Double_t xd212 = 1.38 + 2.23; //mrtva vrstva, pouzivat 1.7 + 2.36
	const Double_t xd223 = 14.; //mrtva vrstva, pouzivat XXXX

	Int_t rn = fT21.FindRing();
	if (rn >= 2) {
		//angles and sectors and rings numbers
		SetTheta2(rn);
		fRN2 = rn;
		fSN2 = fT21.FindSec() / 2;

		//detectors thicknessess
		//vypocitam drahu l v kremiku prvniho detektoru: l = det.d/cos(alfa)
		//	Double_t l = fT11.fThickness/TMath::Cos(fTheta1);
		fX21 = x0_21 / TMath::Cos(fTheta2);
		fX22 = x0_22 / TMath::Cos(fTheta2);
		fXd212 = xd212 / TMath::Cos(fTheta2);
		fXd223 = xd223 / TMath::Cos(fTheta2);

		///////////////
		////T11:T12////
		///////////////

		//verification of the method

		//musim udelat objekt kremikoveho materialu, dost mozna jako clen objektu annulardetector
		//z mereni znam deltaE, vypocitam E0 pred dopadem: E0 = GetE(deltaE)
		fE21v_0r = mSiAlphaV.GetE0dE(fT21.fERing[fRN2]);
		fE21v_0s[fSN2] = mSiAlphaV.GetE0dE(fdE2[fSN2]);
//		if ( (fE11v_0r < 0) || (fE11v_0s[fSN] < 0) ) { return; }

		// E po vystupu z prvni vrstvy kremiku (citliva oblast + mrtva vrstva) bude E = Get( E0 - deltaE, xdead/cos(alfa) )
		//fE12_0 = fE11_0s - fT11.fESec[fSN];	//prvni priblizeni
		fE21vs[fSN2] = fE21v_0s[fSN2] - mSiAlphaV.GetE(fE21v_0s[fSN2], fX21);
		fE21vr = fE21v_0r - mSiAlphaV.GetE(fE21v_0r, fX21);

		fE22v_0[fSN2] = mSiAlphaV.GetE(fE21v_0s[fSN2] - fE21vs[fSN2], fXd212);

		//E muzu zapsat jako vypocitanou energii v druhem kremiku
		fE22v[fSN2] = fE22v_0[fSN2] - mSiAlphaV.GetE(fE22v_0[fSN2], fX22);

		///////////////
		////T12:T13////
		///////////////

		//alha
		fE22_0a[fSN2] = mSiAlpha.GetE0dE(fT22.fESec[fSN2]);
		fE23_0a[fSN2] = mSiAlpha.GetE(fE22_0a[fSN2], fX22);
		fE22a[fSN2] = fE22_0a[fSN2] - fE23_0a[fSN2];
		fE23_0a[fSN2] = mSiAlpha.GetE(fE23_0a[fSN2], fXd223);

		//proton
		fE22_0p[fSN2] = mSiP.GetE0dE(fT22.fESec[fSN2]);
		fE23_0p[fSN2] = mSiP.GetE(fE22_0p[fSN2], fX22);
		fE22p[fSN2] = fE22_0p[fSN2] - fE23_0p[fSN2];
		fE23_0p[fSN2] = mSiP.GetE(fE23_0p[fSN2], fXd223); //after the dead layer

	} //if

	return;

}

void BeEvent::SetSecHitEnergiesSi(Int_t tel) {
	//find hits in telescope nr. tel
	//hits are determined by nonzero energy in two sector in a row of first
	//and second layer of each telescope
	//
	//time should be also taken into account, at present moment, I have no idea where


	const Double_t timemin = 50.;
	const Double_t timemax = 120.;

	Int_t noth = 0;		//number of sector hits
	Int_t tsn[NOHITS];	//time sector number

	//loop over all sectors in second detector
	for (Int_t i = 0; i < NOESEC; i++) {
		if ((fTau[tel][i] > timemin) && (fTau[tel][i] < timemax)) {
			tsn[noth] = i;
			noth++;
		} //if
	} //for


	//find hits
	Int_t nosh = 0;		//number of sector hits
	Int_t sn[noth];	//sector number


	//loop over all sectors
	//if in the sec[i] in T12 and T11 is nonzero signal ==> hit was found
	for (Int_t i = 0; i < noth; i++) {
		if (tel == 0) {
			if (fT12.fESec[tsn[i]] == 0) {
				continue;
			}
			if (fdE1[tsn[i]] == 0) {
				continue;
			}
		} //if tel1
		if (tel == 1) {
			if (fT22.fESec[tsn[i]] == 0) {
				continue;
			}
			if (fdE2[tsn[i]] == 0) {
				continue;
			}
		} //if tel2

		sn[nosh] = i;
		nosh++;
	}

//	printf("%d\t%d\n", noth-nosh, nosh);

	//determine energy of the hits
	Double_t e1s[nosh];		//energy in first detector
	Double_t e2s[nosh];		//energy in second detector
	Int_t e1sindex[nosh];

	//loop over all hits
	for (Int_t i = 0; i < nosh; i++) {
		if (tel == 0) {
			e1s[i] = fdE1[sn[i]];
			e2s[i] = fT12.fESec[sn[i]];
		}
		if (tel == 1) {
			e1s[i] = fdE2[sn[i]];
			e2s[i] = fT22.fESec[sn[i]];
		}
	}
	fNOSHits[tel] = nosh;

	//sort the hits by sector energy
	TMath::Sort(nosh, e1s, e1sindex);
	for (Int_t i = 0; i < nosh; i++) {
		fE1Sec[tel][i] = e1s[e1sindex[i]];
		fE2Sec[tel][i] = e2s[e1sindex[i]];
		fSN[tel][i] = sn[e1sindex[i]]; //?
	}

	return;
}

void BeEvent::SetSecHitEnergiesSiOld(Int_t tel) {
	//find hits in telescope nr. tel
	//hits are determined by nonzero energy in two sector in a row of first
	//and second layer of each telescope
	//
	//time should be also taken into account, at present moment, I have no idea where

	//find hits
	Int_t nosh = 0;		//number of sector hits
	Int_t sn[NOHITS];	//sector number

	//loop over all sectors
	//if in the sec[i] in T12 and T11 is nonzero signal ==> hit was found
	for (Int_t i = 0; i < NOESEC; i++) {
		if (tel == 0) {
			if (fT12.fESec[i] == 0) {
//				cout << "skjdfsdjf" << endl;
				continue;
			}
			if (fdE1[i] == 0) {
				continue;
			}
		} //if tel1
		if (tel == 1) {
			if (fT22.fESec[i] == 0) {
				continue;
			}
			if (fdE2[i] == 0) {
				continue;
			}
		} //if tel2

		sn[nosh] = i;
		nosh++;
	}

	//determine energy of the hits
	Double_t e1s[nosh];		//energy in first detector
	Double_t e2s[nosh];		//energy in second detector
	Int_t e1sindex[nosh];

	//loop over all hits
	for (Int_t i = 0; i < nosh; i++) {
		if (tel == 0) {
			e1s[i] = fdE1[sn[i]];
			e2s[i] = fT12.fESec[sn[i]];
		}
		if (tel == 1) {
			e1s[i] = fdE2[sn[i]];
			e2s[i] = fT22.fESec[sn[i]];
		}
	}
	fNOSHits[tel] = nosh;

//	printf("\n\t%d in sectors, tel %d (fE1Sec)\n", nosh, tel);
	//sort the hits by sector energy
	TMath::Sort(nosh, e1s, e1sindex);
	for (Int_t i = 0; i < nosh; i++) {
		fE1Sec[tel][i] = e1s[e1sindex[i]];
//		printf("%f\t%d\n", fE1Sec[tel][i], i);
		fE2Sec[tel][i] = e2s[e1sindex[i]];
		fSN[tel][i] = sn[e1sindex[i]]; //?
	}
//	cout << endl;

	return;
}

void BeEvent::SetRingHitEnergies(Int_t tel) {
	//find hits in rings of telescope nr. tel, match it with hits in sectors
	//by comparing energy
	//

	if (fNOSHits[tel] == 0) { return; }		//todo povolit az bude fce hotova

	//rings to read setting
	Int_t minring;
	if (tel == 0) {
		minring = 1;
	} //because first two rings are joined
	else {
		minring = 0;
	}

	Int_t norh = 0; 		//number of hits in rings
	Int_t rn[NORINGS]; 		//ring number with "ring hit"

	//loop over all rings to determine number of ring hits and save their ringsnumbers (position)
	for (Int_t i = minring; i < NORINGS; i++) {
		if (tel == 0) { //first telescope
			if (fT11.fERing[i] == 0) {
				continue;
			}
		}
		if (tel == 1) { //second telescope
			if (fT21.fERing[i] == 0) {
				continue;
			}
		}
		rn[norh] = i;
		norh++;
	}

	if (norh > NOHITS) { //it has no sense to have more hits in rings than 16
		return;
	}

	Double_t e1r[norh];			//energy in first detector
	Double_t e1rsorted[norh];			//sorted energy in first detector
	Int_t e1rindex[norh];

	//loop over all hits to read and save energies for each hit
//	printf("\t%d in rings (first search)\n", norh);
	for (Int_t i = 0; i < norh; i++) {
		if (tel == 0) { //first telescope
			e1r[i] = fT11.fERing[rn[i]];
		}
		if (tel == 1) { //second telescope
			e1r[i] = fT21.fERing[rn[i]];
		}
//		printf("%f\t%d\n", e1r[i], rn[i]);
	} //for
//	cout << endl;

	//todo potreba provest serazeni energii v rinzich
	//sort the hits by sector energy
	TMath::Sort(norh, e1r, e1rindex);
//	printf("\t%d in rings (after sort)\n", norh);
	for (Int_t i = 0; i < norh; i++) {
		e1rsorted[i] = e1r[e1rindex[i]];
//		printf("%f\t%f\n", e1r[e1rindex[i]], e1rsorted[i]);
	}
//	cout << endl;

	/////////////////////////////////////////
	//analysis of all events - new attempt
	/////////////////////////////////////////
//	const Double_t con = 0.2;
//	const Double_t con = 0.3;

	//comparing of energies in particular hits (for rings and sectors)
	//and matching them with each other according energy

	Int_t minRingHit = 0;	//maximal ring hit to be verified

	//loop over all hits in sectors
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		//loop over all hits in rings
		for (Int_t j = minRingHit; j < norh; j++) {
//			if (TMath::Abs(fE1Sec[tel][i] - e1rsorted[j]) < con) {
			if ( TMath::AreEqualAbs(fE1Sec[tel][i], e1rsorted[j], fRSTolerance) ) {
				fE1Ring[tel][i] = e1rsorted[j];
				fRN[tel][i] = rn[ e1rindex[j] ];
				minRingHit = j+1;
//				printf("\t%f\t%f\t%f\n", fE1Sec[tel][i], e1rsorted[j], fE1Sec[tel][i] - e1rsorted[j]);
				break;	//sufficient equality was found
			}
		} //for j
	} //for i


//	printf("\t%d in rings (final state)\n", fNOSHits[tel]);
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		//eliminate sector hits without corresponding ring
		if (fRN[tel][i] == -1) {
			fE1Sec[tel][i] = 0.;
			fE2Sec[tel][i] = 0.;
		}
		//eliminate repeated rings
		else {
			for (Int_t j = i + 1; j < fNOSHits[tel]; j++) {
				if (fRN[tel][i] == fRN[tel][j]) {
					if (TMath::Abs(fE1Sec[tel][i] - fE1Ring[tel][i])
							< TMath::Abs(fE1Sec[tel][j] - fE1Ring[tel][j])) {
						fE1Ring[tel][j] = 0.;
					} else {
						fE1Ring[tel][i] = 0.;
					} //else
				}
			} //for
		} //else

//		cout << fE1Ring[tel][i] << endl;

	} //for


//	cout <<"/////////////////////////"<< endl << endl;

	//v tomto miste je vse serazeno a povyrazovano z informaci o rinzich

	//	musim urcite pocet ringu se signalem - fNORHits[tel]
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		if (fE1Ring[tel][i] > 0.) {
			fNORHits[tel]++;
		}
	}

	return;
}

void BeEvent::SetHitsEnergies(Int_t tel) {
	//search hits in telescope nr. tel, fill its energies for particular layers
	//

//	SetSecHitEnergiesSi(tel); //for the first and the second detectors
	SetSecHitEnergiesSiOld(tel); //for the first and the second detectors
	SetRingHitEnergies(tel);

	//loop over all hits in sec
	//counting of generally true hits
	Int_t nohits = 0;
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		if (fE1Ring[tel][i] == 0.) {
			fE1Sec[tel][i] = 0.;
			fE2Sec[tel][i] = 0.;
		} else {
			nohits++;
		}
	}

	Double_t e1sec[fNOSHits[tel]];
	Int_t sn[fNOSHits[tel]];
	Double_t e2sec[fNOSHits[tel]];
	Double_t ering[fNOSHits[tel]];
	Int_t rn[fNOSHits[tel]];
	Int_t e1sindex[fNOSHits[tel]];

	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		e1sec[i] = fE1Sec[tel][i];
		sn[i] = fSN[tel][i];
		e2sec[i] = fE2Sec[tel][i];
		ering[i] = fE1Ring[tel][i];
		rn[i] = fRN[tel][i];
		fE1Sec[tel][i] = 0;
		fSN[tel][i] = -1;
		fE2Sec[tel][i] = 0;
		fE1Ring[tel][i] = 0;
		fRN[tel][i] = -1;
	}

	TMath::Sort(fNOSHits[tel], e1sec, e1sindex);

//	//vypis
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		fE1Sec[tel][i] = e1sec[e1sindex[i]];
		fSN[tel][i] = sn[e1sindex[i]];
		fE2Sec[tel][i] = e2sec[e1sindex[i]];
		fE1Ring[tel][i] = ering[e1sindex[i]];
		fRN[tel][i] = rn[e1sindex[i]];
	}

	fNOHits[tel] = nohits;

	return;
}

void BeEvent::SetEnergyHits(Int_t tel) {
	//search hits in telescope nr. tel, fill its energies for particular layers
	//


	//find hits in telescope nr. tel
	//hits are determined by nonzero energy in two sector in a row of first
	//and second layer of each telescope
	//
	//time should be also taken into account, at present moment, I have no idea where

	//find hits
	Int_t nosh = 0;		//number of sector hits
	Int_t sn[NOHITS];	//sector number

	//loop over all sectors
	//if in the sec[i] in T12 and T11 is nonzero signal ==> hit was found

//	cout << fNOtHits[0] << "\t" << fNOtHits[1] << endl;

	for (Int_t i = 0; i < fNOtHits[tel]; i++) {
		if (fNOtHits[tel] == 0) {
			continue;
		}
		if (tel == 0) {
			if (fT12.fESec[ fTauSN[tel][i] ] == 0) {
//				cout << "null" << endl;
//				printf("null\t%3.2f\t%d\t%3.2f\n", fT12.fESec[ fTauSN[tel][i] ], fNOtHits[tel], fTau[tel][fTauSN[tel][i]]);
				continue;
			}
			if (fdE1[ fTauSN[tel][i] ] == 0) {
//				cout << "null" << endl;
				continue;
			}
		} //if tel1
		if (tel == 1) {
			if (fT22.fESec[ fTauSN[tel][i] ] == 0) {
//				cout << "null" << endl;
				continue;
			}
			if (fdE2[ fTauSN[tel][i] ] == 0) {
//				cout << "null" << endl;
				continue;
			}
		} //if tel2

		sn[nosh] = fTauSN[tel][i];
		nosh++;
	}

//	printf("\t%d\n", nosh);
//	for (Int_t i = 0; i < nosh; i++) {
//		printf("\t\t%d\n", sn[i]);
//	}

	fNOSHits[tel] = nosh;

//	printf("%d\t%d\n", fNOSHits[0], fNOtHits[0]);

//	return;

	//determine energy of the hits
	Double_t e1s[ fNOSHits[tel] ];		//energy in first detector
	Double_t e2s[ fNOSHits[tel] ];		//energy in second detector
	Int_t e1sindex[ fNOSHits[tel] ];

	//loop over all hits
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		if (tel == 0) {
			e1s[i] = fdE1[sn[i]];
			e2s[i] = fT12.fESec[sn[i]];
//			if (e1s[i] == 0) {
//				printf("%3.2f\t%3.2f\t%d\n", e1s[i], e2s[i], sn[i]);
//			}
		}
		if (tel == 1) {
			e1s[i] = fdE2[sn[i]];
			e2s[i] = fT22.fESec[sn[i]];
		}
	}

	//	printf("\n\t%d in sectors, tel %d (fE1Sec)\n", nosh, tel);
	//sort the hits by sector energy
	TMath::Sort(nosh, e1s, e1sindex);
	for (Int_t i = 0; i < nosh; i++) {
		fE1Sec[tel][i] = e1s[e1sindex[i]];
		//		printf("%f\t%d\n", fE1Sec[tel][i], i);
		fE2Sec[tel][i] = e2s[e1sindex[i]];
		fSN[tel][i] = sn[e1sindex[i]]; //?
	}
	//	cout << endl;

//	return;



	if (fNOSHits[tel] == 0) { return; }

	//rings to read setting
	Int_t minring;
	if (tel == 0) {
		minring = 1;
	} //because first two rings are joined
	else {
		minring = 0;
	}

	Int_t norh = 0; 		//number of hits in rings
	Int_t rn[NORINGS]; 		//ring number with "ring hit"

	//loop over all rings to determine number of ring hits and save their ringsnumbers (position)
	for (Int_t i = minring; i < NORINGS; i++) {
		if (tel == 0) { //first telescope
			if (fT11.fERing[i] == 0) {
				continue;
			}
		}
		if (tel == 1) { //second telescope
			if (fT21.fERing[i] == 0) {
				continue;
			}
		}
		rn[norh] = i;
		norh++;
	}

	fNORHits[tel] = norh;

//	for (Int_t i = 0; i < fNORHits[tel]; i++) {
//		cout << rn[i] << endl;
//	}
//	cout << endl;
//	return;
	//vypis:

//	printf("%d\t%d\t%d\n", fNOtHits[tel], fNOSHits[tel], norh);
//	cout << "--------------------------------------" << endl;
//	if (fNOSHits[tel] < fNORHits[tel]) {
//		printf("%d\t%d\t%d\n", fNOtHits[tel], fNOSHits[tel], norh);
//		for (Int_t i = 0; i < fNORHits[tel]; i++) {
//			printf("\n\t\t\t%3.2f\t",fT11.fERing[ rn[i] ]);
//			if (i < fNOSHits[tel]) printf("%3.2f\t", fE1Sec[tel][i]);
//		}
//	}//if vypis


	if (fNORHits[tel] > NOHITS) { //it has no sense to have more hits in rings than 16
//		cout << norh << endl;
		return;
	}

	Double_t e1r[ fNORHits[tel] ];					//energy in first detector
	Double_t e1rsorted[ fNORHits[tel] ];			//sorted energy in first detector
	Int_t e1rindex[ fNORHits[tel] ];

	//loop over all hits to read and save energies for each hit
	for (Int_t i = 0; i < fNORHits[tel]; i++) {
		if (tel == 0) { //first telescope
			e1r[i] = fT11.fERing[rn[i]];
		}
		if (tel == 1) { //second telescope
			e1r[i] = fT21.fERing[rn[i]];
		}
		//		printf("%f\t%d\n", e1r[i], rn[i]);
	} //for

	//todo potreba provest serazeni energii v rinzich
	//sort the hits by sector energy
	TMath::Sort(fNORHits[tel], e1r, e1rindex);
	//	printf("\t%d in rings (after sort)\n", norh);
	for (Int_t i = 0; i < fNORHits[tel]; i++) {
		e1rsorted[i] = e1r[e1rindex[i]];
//				printf("%f\t%f\n", e1r[e1rindex[i]], e1rsorted[i]);
	}
//		cout << endl;

	/////////////////////////////////////////
	//analysis of all events - new attempt
	/////////////////////////////////////////

	//comparing of energies in particular hits (for rings and sectors)
	//and matching them with each other according energy

	Int_t minRingHit = 0;	//maximal ring hit to be verified


	//vypis:

//	cout << "------------------------" << endl;
//	printf("\t\t\t%d\t%d\t%d\n", fNOtHits[tel], fNOSHits[tel], fNORHits[tel]);
//
//	printf("\t\t\t%d\t%d\t%d\n", fNOtHits[tel], fNOSHits[tel], norh);
//	cout << "--------------------------------------" << endl;
////	if (fNOSHits[tel] < fNORHits[tel]) {
////		printf("%d\t%d\t%d\n", fNOtHits[tel], fNOSHits[tel], norh);
//		for (Int_t i = 0; i < fNORHits[tel]; i++) {
//			printf("\nring\t%d\tringE:\t%3.2f\t", rn[i], e1rsorted[i]);
//			if (i < fNOSHits[tel]) printf("sector: %3.2f\t", fE1Sec[tel][i]);
//		}
//	}//if vypis

//return;



//	cout << "-------" << endl;
	//loop over all hits in sectors
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		//loop over all hits in rings
		for (Int_t j = minRingHit; j < fNORHits[tel]; j++) {
			if ( TMath::AreEqualAbs(fE1Sec[tel][i], e1rsorted[j], fRSTolerance) ) {
				fE1Ring[tel][i] = e1rsorted[j];
				fRN[tel][i] = rn[ e1rindex[j] ];
				minRingHit = j+1;
//				printf("\t%f\t%f\t%f\n", fE1Sec[tel][i], fE1Ring[tel][i], fE1Sec[tel][i] - fE1Ring[tel][i]);
				break;	//sufficient equality was found
			}
		} //for j
	} //for i
//	cout << endl;







//	printf("\t%d\n", fNOSHits[tel]);
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
//		cout << i << endl;
		//eliminate sector hits without corresponding ring
		if (fRN[tel][i] == -1) {
			fE1Sec[tel][i] = 0.;
			fE2Sec[tel][i] = 0.;
		}
		//eliminate repeated rings
		else {
			for (Int_t j = i + 1; j < fNOSHits[tel]; j++) {
//				cout << "hovno" << endl;
				if (fRN[tel][i] == fRN[tel][j]) {
//					cout << "hovno" << endl;
					Warning("BeEvent::SetEnergyHits", "some strange part of code was used: check it!!!!");
					if (TMath::Abs(fE1Sec[tel][i] - fE1Ring[tel][i])
					< TMath::Abs(fE1Sec[tel][j] - fE1Ring[tel][j])) {
						fE1Ring[tel][j] = 0.;
					} else {
						fE1Ring[tel][i] = 0.;
					} //else
				}//if
			} //for
		} //else

		//		cout << fE1Ring[tel][i] << endl;

	} //for

//	cout << endl;
//
//	return;

//	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
//		printf("\t%d\t%d\t%3.2f\t%3.2f\n", fSN[tel][i], fRN[tel][i], fE1Sec[tel][i], fE1Ring[tel][i]);
//	}
//	cout << endl;


	//	cout <<"/////////////////////////"<< endl << endl;

	//v tomto miste je vse serazeno a povyrazovano z informaci o rinzich

//	cout << fNORHits[tel] << endl;
//	printf("\t\t%d\t%d\n", fNORHits[tel], fNOSHits[tel]);
//
//	//	musim urcite pocet ringu se signalem - fNORHits[tel]
//	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
//		if (fE1Ring[tel][i] > 0.) {
//			fNORHits[tel]++;
//		}
//	}
//
//	printf("\t%d\n", fNORHits[tel]);
//
//	return;





	//loop over all hits in sec
	//counting of generally true hits
	Int_t nohits = 0;
	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		if (fE1Ring[tel][i] > 0.) {
//			printf("\t%3.2f\t%3.2f\n", fE1Sec[tel][i], fE2Sec[tel][i]);
//			fE1Sec[tel][i] = 0.;
//			fE2Sec[tel][i] = 0.;
//		} else {
			nohits++;
		}//if
	}//for

	fNOHits[tel] = nohits;

//	printf("\t%d\t%d\n", fNOSHits[tel], fNOHits[tel]);
//	return;

//	/*
	//choice of proper hits
	Double_t e1secE[fNOSHits[tel]];
	Int_t secN[fNOSHits[tel]];
	Double_t e2secE[fNOSHits[tel]];
	Double_t eringE[fNOSHits[tel]];
	Int_t ringN[fNOSHits[tel]];
	Int_t e1secIndex[fNOSHits[tel]];

	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
		e1secE[i] = fE1Sec[tel][i];
		secN[i] = fSN[tel][i];
		e2secE[i] = fE2Sec[tel][i];
		eringE[i] = fE1Ring[tel][i];
		ringN[i] = fRN[tel][i];
//		printf("before:\t\t\t%d\n", fRN[tel][i]);
	}

	//small reset
	for (Int_t i = 0; i < NOHITS; i++) {
		fE1Sec[tel][i] = 0;
		fSN[tel][i] = -1;
		fE2Sec[tel][i] = 0;
		fE1Ring[tel][i] = 0;
		fRN[tel][i] = -1;
	}

//	return;
	TMath::Sort(fNOSHits[tel], e1secE, e1secIndex);

//	//vypis
//	for (Int_t i = 0; i < fNOSHits[tel]; i++) {
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
		fE1Sec[tel][i] = e1secE[e1secIndex[i]];
		fSN[tel][i] = secN[e1secIndex[i]];
		fE2Sec[tel][i] = e2secE[e1secIndex[i]];
		fE1Ring[tel][i] = eringE[e1secIndex[i]];
		fRN[tel][i] = ringN[e1secIndex[i]];
//		printf("\t%3.2f\t%3.2f\n", fE1Sec[tel][i], fE1Ring[tel][i]);
//		printf("after:\t%d\n", fRN[tel][i]);
	}
//	cout << endl;



	return;
//	*/
}

void BeEvent::SetThetas(Int_t tel) {
	//set theta angles in lab system for each hit

	if (fNOHits[tel] == 0) {
		return;
	}

	TVector3 trajectory;

	//loop over all hits to set angle theta for each of them
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
		if (tel == 0) {
//			if (fRN[tel][i] < 2) cout << fRN[tel][i] << endl;
//			if (fRN[tel][i] != 1) { //two inner rings in first telescope are joined
//			if (fRN[tel][i] != 0) { //two inner rings in first telescope are joined
//				Double_t detx = fT11.GetX(fRN[tel][i], fSN[tel][i]);
//				Double_t dety = fT11.GetY(fRN[tel][i], fSN[tel][i]);
				Double_t detx = fT11.GetX(fRN[tel][i], GetNOSecPhi(tel, fSN[tel][i]));
				Double_t dety = fT11.GetY(fRN[tel][i], GetNOSecPhi(tel, fSN[tel][i]));

//				printf("%f\t%f\t%f\n", detx, dety, fT11.GetZPosition());
//				printf("%f\t%f\t%f\n\n", detx - B_X, dety - B_Y, fT11.GetZPosition());

//				trajectory.SetXYZ(detx, dety, fT11.GetZPosition());
//				printf("%f\n", trajectory.Theta());

//				trajectory.SetXYZ(detx - B_X*10, dety - B_Y*10, fT11.GetZPosition());	//cm to mm
				trajectory.SetXYZ(detx - fBeamX*10, dety - fBeamY*10, fT11.GetZPosition());	//cm to mm
//				printf("%f\n\n", trajectory.Theta());
				fTheta[tel][i] = trajectory.Theta();
//			} //if
		} //if
		if (tel == 1) {
//			Double_t detx = fT21.GetX(fRN[tel][i], fSN[tel][i]);
//			Double_t dety = fT21.GetY(fRN[tel][i], fSN[tel][i]);
			Double_t detx = fT21.GetX(fRN[tel][i], GetNOSecPhi(tel, fSN[tel][i]));
			Double_t dety = fT21.GetY(fRN[tel][i], GetNOSecPhi(tel, fSN[tel][i]));

//			printf("%f\t%f\t%f\n", detx, dety, fT21.GetZPosition());
//			printf("%f\t%f\t%f\n\n", detx - B_X, dety - B_Y, fT21.GetZPosition());

//			trajectory.SetXYZ(detx, dety, fT21.GetZPosition());
//			printf("%f\n", trajectory.Theta());

//			trajectory.SetXYZ(detx, dety, fT21.GetZPosition());
			trajectory.SetXYZ(detx - fBeamX*10, dety - fBeamY*10, fT21.GetZPosition());		//cm to mm
//			printf("%f\n\n", trajectory.Theta());
			fTheta[tel][i] = trajectory.Theta();
		}
	} //for

	return;
}

void BeEvent::SetThetasOld(Int_t tel) {
	//set theta angles in lab system for each hit

	if (fNOHits[tel] == 0) {
		return;
	}

	//loop over all hits to set angle theta for each of them
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
		if (tel == 0) {
			if (fRN[tel][i] != 1) { //two inner rings in first telescope are joined
				fTheta[tel][i] = fT11.GetTheta(fRN[tel][i]);
			} //if
		} //if
		if (tel == 1) {
			fTheta[tel][i] = fT21.GetTheta(fRN[tel][i]);
		}
	} //for

	return;
}

void BeEvent::SetThetasNew(Int_t tel) {
	//set theta angles in lab system for each hit

	if (fNOHits[tel] == 0) {
		return;
	}

	//loop over all hits to set angle theta for each of them
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
		if (tel == 0) {
			if (fRN[tel][i] != 1) { //two inner rings in first telescope are joined
				fTheta[tel][i] = fT11.GetTheta(fRN[tel][i]); //should be redone

//				Double_t theta;
//				theta = TMath::ATan(( fT11.fInnerR + ( fRN[tel][i]+ranTheta.Uniform(1) )*
//						(fT11.fOuterR - fT11.fInnerR)/fT11.fRingsNumber )/fT11.fZ);

			} //if
		} //if
		if (tel == 1) {
			fTheta[tel][i] = fT21.GetTheta(fRN[tel][i]);
		}
	} //for

	return;
}

void BeEvent::SetPhis(Int_t tel) {
	//set phi angle for each hit

	if (fNOHits[tel] == 0) {
		return;
	}

	TVector3 trajectory;

	//loop over all right hits to set angle theta for each of them
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
		if (tel == 0) {
			//puvodni
//			fPhi[tel][i] = fT11.GetPhi( GetNOSecPhi(tel, fSN[tel][i]) );
			//nova verze

//			Double_t detx = fT11.GetX(fRN[tel][i], fSN[tel][i]);
//			Double_t dety = fT11.GetY(fRN[tel][i], fSN[tel][i]);
			Double_t detx = fT11.GetX(fRN[tel][i], GetNOSecPhi(tel, fSN[tel][i]));
			Double_t dety = fT11.GetY(fRN[tel][i], GetNOSecPhi(tel, fSN[tel][i]));

			trajectory.SetXYZ(detx - fBeamX*10, dety - fBeamY*10, fT11.GetZPosition());			//cm to mm
//			cout << fBeamY << endl;
//			cout << fT11.GetZPosition() << endl;
			fPhi[tel][i] = trajectory.Phi();

//			printf("sector: %d;\tgetnosec: %d;\tphi_wr: %3.2f\n", fSN[tel][i], GetNOSecPhi(tel, fSN[tel][i]), fPhi[tel][i]);
//			printf("sector: %d;\tx: %3.2f;\ty %3.2f;\tphi_wr: %3.2f\n", fSN[tel][i], detx, dety, fPhi[tel][i]*TMath::RadToDeg());
//			printf("sector: %d;\tphi_wr: %3.2f\n", fSN[tel][i], fPhi[tel][i]*TMath::RadToDeg());

//			fPhi[tel][i] = 2 * TMath::Pi()
//					- fT11.GetPhi(GetNOSecPhi(tel, fSN[tel][i]))
//					- TMath::Pi() / 2.;
		} //if
		if (tel == 1) {
//			fPhi[tel][i] = fT21.GetPhi( GetNOSecPhi(tel, fSN[tel][i]) );

//			Double_t detx = fT21.GetX(fRN[tel][i], fSN[tel][i]);
//			Double_t dety = fT21.GetY(fRN[tel][i], fSN[tel][i]);
			Double_t detx = fT21.GetX(fRN[tel][i], GetNOSecPhi(tel, fSN[tel][i]));
			Double_t dety = fT21.GetY(fRN[tel][i], GetNOSecPhi(tel, fSN[tel][i]));

			trajectory.SetXYZ(detx - fBeamX*10, dety - fBeamY*10, fT21.GetZPosition());			//cm to mm

			fPhi[tel][i] = trajectory.Phi();

//			fPhi[tel][i] = 2 * TMath::Pi()
//					- fT21.GetPhi(GetNOSecPhi(tel, fSN[tel][i]))
//					- TMath::Pi() / 2.;
		} //if
	} //for

	return;
}

Int_t BeEvent::GetNOSecPhi(Int_t tel, Int_t esec) {
	/////////////////////////////////////////
	// function return number of physical
	// sector in first telescope of each
	// detector, esec is number of energy
	// sector (from the second telescope
	/////////////////////////////////////////

	if (tel == 0) {
		if (fT11.fESec[2 * esec] > 0) {
			return 2 * esec;
		}
		if (fT11.fESec[2 * esec + 1] > 0) {
			return 2 * esec + 1;
		}
	}

	if (tel == 1) {
		if (fT21.fESec[2 * esec] > 0) {
			return 2 * esec;
		}
		if (fT21.fESec[2 * esec + 1] > 0) {
			return 2 * esec + 1;
		}
	}

	return -1;
}

void BeEvent::SetCalCorrs() {
	//CsI calibration corrections
	//fE23aSec:
	/*	fCalCorr[0] = 0.998;
		 fCalCorr[1] = 0.995;
		 fCalCorr[2] = 1.004;
		 fCalCorr[3] = 0.984;
		 fCalCorr[4] = 1.005;
		 fCalCorr[5] = 0.985;
		 fCalCorr[6] = 0.99;
		 fCalCorr[7] = 1.013;
		 fCalCorr[8] = 0.983;
		 fCalCorr[9] = 1.005;
		 fCalCorr[10] = 1.;
		 fCalCorr[11] = 1.01;
		 fCalCorr[12] = 1.01;
		 fCalCorr[13] = 1.;
		 fCalCorr[14] = 1.015;
		 fCalCorr[15] = 0.99;
		 //*/
	//fE23pSec:
	/*	fCalCorr[0] = 0.952;
		 fCalCorr[1] = 1.022;
		 fCalCorr[2] = 0.988;
		 fCalCorr[3] = 0.99;
		 fCalCorr[4] = 0.993;
		 fCalCorr[5] = 0.963;
		 fCalCorr[6] = 0.963;
		 fCalCorr[7] = 1.025;
		 fCalCorr[8] = 0.975;
		 fCalCorr[9] = 1.02;
		 fCalCorr[10] = 1.02;
		 fCalCorr[11] = 1.04;
		 fCalCorr[12] = 1.04;
		 fCalCorr[13] = 1.015;
		 fCalCorr[14] = 1.04;
		 fCalCorr[15] = 0.98;
		 //*/

	//fE13aSec:
	/*	fCalCorr[0] = 1.;
		 fCalCorr[1] = 0.996;
		 fCalCorr[2] = 1.005;
		 fCalCorr[3] = 1.009;
		 fCalCorr[4] = 0.989;
		 fCalCorr[5] = 1.007;
		 fCalCorr[6] = 1.;
		 fCalCorr[7] = 1.014;
		 fCalCorr[8] = 1.013;
		 fCalCorr[9] = 0.992;
		 fCalCorr[10] = 0.985;
		 fCalCorr[11] = 0.987;
		 fCalCorr[12] = 0.988;
		 fCalCorr[13] = 0.948;
		 fCalCorr[14] = 0.978;
		 fCalCorr[15] = 0.979;
		 //*/
	//silicon
	//fE13pSec
	fCalCorr[0] = 0.968;
	fCalCorr[1] = 1.018;
	fCalCorr[2] = 1.016;
	fCalCorr[3] = 1.025;
	fCalCorr[4] = 0.956;
	fCalCorr[5] = 1.022;
	fCalCorr[6] = 0.982;
	fCalCorr[7] = 1.033;
	fCalCorr[8] = 0.986;
	fCalCorr[9] = 0.926;
	fCalCorr[10] = 0.955;
	fCalCorr[11] = 0.948;
	fCalCorr[12] = 0.98;
	fCalCorr[13] = 0.905;
	fCalCorr[14] = 0.944;
	fCalCorr[15] = 0.97;
}

void BeEvent::BeParticleIdentificationNew(Int_t tel) {
	//loop over all hits
	//this method is improved BeParticleIdentification and use graphical cuts

	if (cutAlpha1 == NULL) {
		Error("BeEvent::BeParticleIdentificationNew",
				"cutAlpha1 was not initialized");
		return;
	}
	if (cutProton1 == NULL) {
		Error("BeEvent::BeParticleIdentificationNew",
				"cutProton1 was not initialized");
		return;
	}
	if (cutAlpha2 == NULL) {
		Error("BeEvent::BeParticleIdentificationNew",
				"cutAlpha2 was not initialized");
		return;
	}
	if (cutProton2 == NULL) {
		Error("BeEvent::BeParticleIdentificationNew",
				"cutProton2 was not initialized");
		return;
	}

//	static Int_t alphas1 = 0;
//	static Int_t protons1 = 0;
//	static Int_t alphas2 = 0;
//	static Int_t protons2 = 0;

//	printf("\t%d\t%d\n", fNOaHits[tel], fNOpHits[tel]);

	for (Int_t i = 0; i < fNOHits[tel]; i++) {

		//first telescope
		if (tel == 0) {
			//alphas
//			printf("\t%3.2f\t%3.2f\n", fT12.fESec[fSN[tel][i]], fE2Sec[tel][i]);
//			printf("\t%3.2f\t%3.2f\n", fdE1[fSN[tel][i]], fE1Sec[tel][i]);
//			printf("\t%3.2f\t%3.2f\n", fE13aSec[fSN[tel][i]], fE2Sec[tel][i]);
			if (cutAlpha1->IsInside(
					fT12.fESec[fSN[tel][i]] + fE13aSec[fSN[tel][i]],
					fdE1[fSN[tel][i]])) {
//				cout << "alpha1" << endl;
				fID[tel][i] = kALPHA;
				fNOaHits[tel]++;
				fE3Sec[tel][i] = fE13aSec[fSN[tel][i]];
//				printf("\t%3.2f\t%3.2f\n", fE3Sec[tel][i], fE2Sec[tel][i]);

//				alphas1++;
//				if (particle==kPROTON) { fNOpHits[tel]++; }
//				return;
			}
			//protons
			if (cutProton1->IsInside(
					fT12.fESec[fSN[tel][i]] + fE13pSec[fSN[tel][i]],
					fdE1[fSN[tel][i]])) {
//				cout << "proton 1" << endl;
				fID[tel][i] = kPROTON;
				fNOpHits[tel]++;
				fE3Sec[tel][i] = fE13pSec[fSN[tel][i]];
//				protons1++;
			}
		} //if first telescope

		//second telescope
		if (tel == 1) {
			//alphas
			if (cutAlpha2->IsInside(
					fT22.fESec[fSN[tel][i]] + fE23aSec[fSN[tel][i]],
					fdE2[fSN[tel][i]])) {
//				cout << "alpha 2" << endl;
				fID[tel][i] = kALPHA;
				fNOaHits[tel]++;
				fE3Sec[tel][i] = fE23aSec[fSN[tel][i]];
//				alphas2++;
			}
			//protons
			if (cutProton2->IsInside(
					fT22.fESec[fSN[tel][i]] + fE23pSec[fSN[tel][i]],
					fdE2[fSN[tel][i]])) {
//				cout << "proton 2" << endl;
				fID[tel][i] = kPROTON;
				fNOpHits[tel]++;
				fE3Sec[tel][i] = fE23pSec[fSN[tel][i]];
//				protons2++;
			}
		} //if second telescope

	} //for: loop over all hits

//	cout << endl;
//	cout << alphas1 << endl;
//	cout << protons1 << endl;
//	cout << alphas2 << endl;
//	cout << protons2 << endl;

	return;
}

void BeEvent::BeParticleIdentification(Int_t tel) {
	//loop over all hits
	//parameters in function BeEvent::HitIdentification(...) are sensitively set
	//for our purposes and need not to be touched

	for (Int_t i = 0; i < fNOHits[tel]; i++) {

		//first telescope
		if (tel == 0) {
			HitIdentification(tel, i, fE13aSec[fSN[tel][i]], kALPHA, -0.9, 1.2,
					-0.9, 1.2, -0.7, 0.6);
//			HitIdentification(tel, i, fE13aSec[fSN[tel][i]], kALPHA, -0.9, 1.2, -9, 2, -7, 6);
//			HitIdentification(tel, i, fE13pSec[fSN[tel][i]], kPROTON, -0.5, 0.5, -0.5, 1.03, -0.4, 0.55);
			HitIdentification(tel, i, fE13pSec[fSN[tel][i]], kPROTON, -0.5, 0.5,
					-0.0, 1.03, -0.4, 0.55);
		}
		//second telescope
		if (tel == 1) {
			HitIdentification(tel, i, fE23aSec[fSN[tel][i]], kALPHA, -0.5, 1.2,
					-0.1, 2.4, -0.1, 1.0); //ok
			HitIdentification(tel, i, fE23pSec[fSN[tel][i]], kPROTON, -0.2, 0.4,
					-0.1, 1.5, -0.1, 0.8); //ok
		}

	} //for: loop over all hits

	return;
}

void BeEvent::ParticleMarkers() {
	Int_t k = 0;
	Int_t a = 1;
	Int_t p = 1;

	for (Int_t i = 0; i < NOTEL; i++) {
//		for (Int_t j = 0; j < NOHITS; j++) {
		for (Int_t j = 0; j < fNOHits[i]; j++) {
//			cout << fNOHits[i] << endl;
			if (fID[i][j] == kALPHA) {
//				cout << i << endl;
				fPM[k] = kALPHA * 10 + a;
//				cout << fPM[k] << endl;
				k++;
				a++;
			}
			if (fID[i][j] == kPROTON) {
				fPM[k] = kPROTON * 10 + p;
//				cout << fPM[k] << endl;
				k++;
				p++;
			}
		} //for: : loop over all hits
	} //for: : loop over all telescopes

//	cout << endl;
	return;
}

void BeEvent::SetTimeHits(Int_t tel) {
//	if (fNOHits[tel] == 0) { return; }

//	Double_t time;
	const Double_t timemin = 50.;
	const Double_t timemax = 120.;

	//loop over all hits
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
//		time = fTau[tel][fSN[tel][i]];
		if ((fTau[tel][fSN[tel][i]] > timemin)
				&& (fTau[tel][fSN[tel][i]] < timemax)) {
			fTauSN[tel][i] = fSN[tel][i];
			fhTau[tel][i] = fTau[tel][fSN[tel][i]];
			fNOtHits[tel]++;
		} //if
	} //for

	return;
}

void BeEvent::SetTimeHitsNew(Int_t tel) {

	const Double_t timemin = 50.;
	const Double_t timemax = 120.;

	Int_t j = 0;

//	cout << fNOtHits[tel] << endl;

	for (Int_t i = 0; i < NOESEC; i++) {
//		cout << fTau[tel][i] << endl;
		if (fTau[tel][i] < timemin || fTau[tel][i] > timemax) {
//			cout << "continue" << endl;
			continue;
		}
		fNOtHits[tel]++;
		fTauSN[tel][j] = i;
		j++;
	}

//	cout << "\t" << fNOtHits[tel] << endl;
	return;
}

void BeEvent::HitIdentification(const Int_t tel, const Int_t hit,
		const Double_t eCsI, const Int_t particle, const Double_t amin,
		const Double_t amax, const Double_t bmin, const Double_t bmax,
		const Double_t cmin, const Double_t cmax) {
	////////////////////////////////////////////////////////////////////////////
	///
	///	Function for indetification hit #hit in telescope #tel, where
	///	eCsI is energy in third detector of the telescope, particle is
	///	particle identificator, amin and amax bounds for comparision of
	/// the dE_ex and dE_calc when eCsI==0, bmin and bmax when eCsI!=0
	///	and cmin and cmax for final cleanup of the event to be identified
	///
	////////////////////////////////////////////////////////////////////////////

//	Double_t x;
	Double_t x = 0.;
	Double_t econ;

	//case when there is no signal in CsI
//	if (eCsI==0) {
	//compare experimental and calculated dE in Tx1
	if (tel == 0) {
		x = (fT11.GetThickness() + fT12.GetDeadLayerThickness())
				/ TMath::Cos(fTheta[tel][hit]);
	}
	if (tel == 1) {
		x = (fT21.GetThickness() + fT22.GetDeadLayerThickness())
				/ TMath::Cos(fTheta[tel][hit]);
	}
	econ = CompareLosses(fE1Sec[tel][hit], fE2Sec[tel][hit], x, particle);
	//	fdE1calc[tel][hit] = econ;	//docasne
	//particle identification
	if (econ > amin && econ < amax) {
		fID[tel][hit] = particle;
		if (particle == kALPHA) {
			fNOaHits[tel]++;
		}
		if (particle == kPROTON) {
			fNOpHits[tel]++;
		}
		return; //docasne	//fixme why "docasne"
	}

//	}//if without CsI signal

	//signal in CsI
//	if (eCsI>0) {
	//compare experimental and calculated dE in Tx1
	x = (fT12.GetThickness() + fT13.GetDeadLayerThickness())
			/ TMath::Cos(fTheta[tel][hit]);
	econ = CompareLosses(fE2Sec[tel][hit], eCsI, x, particle);
	//	fdE1calc[tel][hit] = econ;	//docasne
	//particle identification
	if (econ > bmin && econ < bmax) {
		//final cleanup
		//determine dE from energy in Tx2+Tx3
		x = (fT11.GetThickness() + fT12.GetDeadLayerThickness())
				/ TMath::Cos(fTheta[tel][hit]);
		econ = CompareLosses(fE1Sec[tel][hit], fE2Sec[tel][hit] + eCsI, x,
				particle);
		fdE1calc[tel][hit] = econ; //docasne
		if (econ > cmin && econ < cmax) {
			fID[tel][hit] = particle;
			if (particle == kALPHA) {
				fNOaHits[tel]++;
			}
			if (particle == kPROTON) {
				fNOpHits[tel]++;
			}
			fE3Sec[tel][hit] = eCsI;
		} //if
	} //if

//	}//if CsI signal
//*/
	return;
}

Int_t BeEvent::NOParticlesTauE(Int_t tel, Int_t particleID) {
	Int_t noproperhits = 0;

	//loop over all hits
	for (Int_t i = 0; i < fNOHits[tel]; i++) {

//		if ( fTauSN[tel][i] == -1 )	cout << "hovno" << endl;
		if (fID[tel][i] == particleID && fTauSN[tel][i] != -1) {
//			cout << "hovno" << endl;
			noproperhits++;
		}
	}

	return noproperhits;
}

Double_t BeEvent::CompareLosses(Double_t de_ex, Double_t e, Double_t x,
		Int_t particleID) {

	if ((particleID != kALPHA) && (particleID != kPROTON)) {
		return -100.;
	}

	Double_t e0;

	//alphas
	if (particleID == kALPHA) {
		e0 = mSiAlpha.GetE0(e, x);
	}
	//protons
	if (particleID == kPROTON) {
		e0 = mSiP.GetE0(e, x);
	}

	return (e0 - e) - de_ex;
}

void BeEvent::SetHitsOld() {

	//first telescope
	SetHitsEnergies(0); //ok	search for hits in telescopes and set their energy
	SetThetas(0); //ok	set thetas of found hits
	if (fCuts) {
		BeParticleIdentificationNew(0);
	} else {
		BeParticleIdentification(0); //		try to identify particles as protons or alphas
	}
	SetPhis(0); //ok
	SetT(0); //ok
	SetP(0); //ok
	SetTimeHits(0); //ok

	//second telescope
	SetHitsEnergies(1);
	SetThetas(1);
	if (fCuts) {
		BeParticleIdentificationNew(1);
	} else {
		BeParticleIdentification(1); //		try to identify particles as protons or alphas
	}
	SetPhis(1);
	SetT(1);
	SetP(1);
	SetTimeHits(1);

	//all telescopes together
	ParticleMarkers();
//	//invariant mass
	SetProductsLab();
//	//beam energy calibration
	SetQLiP();
//	//lab --> CM_Li-p
	SetProductsCM1();
//	//lab --> CM_Be
	SetProductsCM();

	return;
}

void BeEvent::SetHits() {

	for (Int_t i = 0; i < NOTEL; i++) {
		SetTimeHitsNew(i); //ok
		SetEnergyHits(i);

		//to be checked
		SetThetas(i);
		if (fCuts) {
			BeParticleIdentificationNew(i);
		} else {
			BeParticleIdentification(i); //		try to identify particles as protons or alphas
		}
		SetPhis(i); //ok
		SetT(i); //ok
//		SetTNew(i);
		SetP(i); //ok
	}

	ParticleMarkers();
	//	//invariant mass
	SetProductsLabNew();
	//	//beam energy calibration
	SetQLiP();
	//	//lab --> CM_Li-p
	SetProductsCM1();
	//	//lab --> CM_Be
	SetProductsCM();

	return;
}

void BeEvent::SetT(Int_t tel) {
//Set kinetic energy of particle corrected to energy losses in detector's
//	dead layers and material of the target and target window.

	Double_t cos_theta = 1;

	//loop over all hits
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
//		printf("\t%d\t%d\t%d\n", tel, fNOHits[tel], fID[tel][i]);
		if (fID[tel][i] == kPROTON) {
			if (tel == 0) { //first telescope
				cos_theta = TMath::Cos(fTheta[tel][i]);
//				printf("\t%3.2f\t%3.2f\n", fTheta[tel][i]*TMath::RadToDeg(), cos_theta);
//				cout << cos_theta << endl;
				fT[tel][i] = CalcKinE(tel, i, cos_theta, fX13_FD, fX12_BD, fX12_FD,
						fX11_BD, fX11_FD);
			}
			if (tel == 1) { //second telescope
				cos_theta = TMath::Cos(fTheta[tel][i]);

				fT[tel][i] = CalcKinE(tel, i, cos_theta, fX23_FD, fX22_BD, fX22_FD,
						fX21_BD, fX21_FD);
			}
		}
		if (fID[tel][i] == kALPHA) { //first telescope
			if (tel == 0) {
				cos_theta = TMath::Cos(fTheta[tel][i]);

				fT[tel][i] = CalcKinE(tel, i, cos_theta, fX13_FD, fX12_BD, fX12_FD,
						fX11_BD, fX11_FD);

			}
			if (tel == 1) { //second telescope
				cos_theta = TMath::Cos(fTheta[tel][i]);

				fT[tel][i] = CalcKinE(tel, i, cos_theta, fX23_FD, fX22_BD, fX22_FD,
						fX21_BD, fX21_FD);

			}
		}
	}//for

//	cout << endl;

	return;
}

void BeEvent::SetTNew(Int_t tel) {
//Set kinetic energy of particle corrected to energy losses in detector's
//	dead layers and material of the target and target window.

	//loop over all hits
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
		if (tel == 0) { //first telescope
			fT[tel][i] = CalcKinE(tel, i, TMath::Cos(fTheta[tel][i]), fX13_FD,
					fX12_BD, fX12_FD, fX11_BD, fX11_FD);
		}//if 1st tel
		if (tel == 1) { //second telescope
			fT[tel][i] = CalcKinE(tel, i, TMath::Cos(fTheta[tel][i]), fX23_FD,
					fX22_BD, fX22_FD, fX21_BD, fX21_FD);
		}//if 2nd tel
	}//for

	return;
}

Double_t BeEvent::CalcKinE(const Int_t telescope, const Int_t hit,
		const Double_t cosTheta, const Double_t fDeadCsI,
		const Double_t bDeadSi2, const Double_t fDeadSi2,
		const Double_t bDeadSi1, const Double_t fDeadSi1) {
	//	cosTheta
	//	fDeadCsI:	front dead layer of CsI
	//	bDeadSi2:	back dead layer of second Si
	//	fDeadSi2:	front dead layer of second Si
	//	bDeadSi1:	back dead layer of first Si
	//	fDeadSi1:	front dead layer of first Si

	Double_t kinE = 0;

	TELoss *silosses;
	TELoss *csilosses;
	TELoss *targetlosses;
	TELoss *target_win_losses;

	if (fID[telescope][hit] == kPROTON) { //proton
//		kinE0 = particle->E() - kPMASS;
		silosses = &mSiP;
		csilosses = &mCsIP;
		targetlosses = &mH_P;
		target_win_losses = &mWin_P;
	} else {
		if (fID[telescope][hit] == kALPHA) { //alpha
//			kinE0 = particle->E() - kAMASS;
			silosses = &mSiAlpha;
			csilosses = &mCsIAlpha;
			targetlosses = &mH_Alpha;
			target_win_losses = &mWin_Alpha;
		} else {
			Info("BeWork::ParticleTracking",
					"I can only track either protons or alphas.");
			return 0;
		}
	} //else

	//correction on dead layer in CsI detector:
	kinE = csilosses->GetE0(fE3Sec[telescope][hit], (fDeadCsI) / cosTheta);
	//correction on back dead layer of the second detector:
	kinE = silosses->GetE0(kinE, bDeadSi2 / cosTheta);

	kinE = kinE + fE2Sec[telescope][hit];
	//correction on front dead layer of the second detector:
	kinE = silosses->GetE0(kinE, fDeadSi2 / cosTheta);
	//correction on back dead layer of the first detector:
	kinE = silosses->GetE0(kinE, bDeadSi1 / cosTheta);

	kinE = fE1Sec[telescope][hit] + kinE;
	//correction on front dead layer of the first detector:
	kinE = silosses->GetE0(kinE, fDeadSi1 / cosTheta);

	//losses in target window:
	kinE = target_win_losses->GetE0(kinE, T_WIN_THICK / cosTheta);

	//losses in target window:
	kinE = targetlosses->GetE0(kinE, T_THICKNESS / 2. / cosTheta);

	return kinE;
}

void BeEvent::SetP(Int_t tel) {
	for (Int_t i = 0; i < fNOHits[tel]; i++) {
		if (fID[tel][i] == kPROTON) {
			fP[tel][i] = TMath::Power(PC2(fT[tel][i], kPMASS), 0.5);
		}
		if (fID[tel][i] == kALPHA) {
			fP[tel][i] = TMath::Power(PC2(fT[tel][i], kAMASS), 0.5);
		}
	}

	return;
}

void BeEvent::SetProductsLab() {
	fNOA = NOParticlesTauE(0, kALPHA) + NOParticlesTauE(1, kALPHA);
	fNOP = NOParticlesTauE(0, kPROTON) + NOParticlesTauE(1, kPROTON);

	if (fNOtHits[0] + fNOtHits[1] != 3) {
		return;
	}

//	if (NOParticlesTauE(0, kALPHA)+NOParticlesTauE(1, kALPHA) != 1) { return; }
//	if (NOParticlesTauE(0, kPROTON)+NOParticlesTauE(1, kPROTON) != 2) { return; }
	if (fNOA != 1) {
		return;
	}
	if (fNOP != 2) {
		return;
	}

	//loop over all hits to determine number of alphas and protons
//	fce ParticleNumberTau()

	Int_t atel = 0;
	Int_t ahit = 0;

	Int_t ptel[2] = { 0, 0 };
	Int_t phit[2] = { 0, 0 };
	Int_t nop = 0; //counter of protons
	Int_t noa = 0; //counter of alphas

	//loop over all telescopes
	for (Int_t i = 0; i < NOTEL; i++) {
		//loop over all hits in each telescope
		for (Int_t j = 0; j < fNOHits[i]; j++) {
			if (fTauSN[i][j] == -1) {
				continue;
			}
			if (fID[i][j] == kPROTON) {
				//zapamatovat si ho jako proton
				ptel[nop] = i;
				phit[nop] = j;
				nop++;
			}
			if (fID[i][j] == kALPHA) {
				//zapamatovat si ho jako alfu
				atel = i;
				ahit = j;
				noa++;
			}
		} //for over all hits in each telescope
	} //for over all telescopes

	//sum(m_i*c^2) + Tcm = sqrt( ( sum(m_c*c2 + T_i) )^2 - (p_system*c)^2 )

	Double_t Elab = 0; //Elab = sum(m_c*c2 + T_i)
	Double_t PClab = 0; //PClab = p_system*c, where p_system = p1lab + p2lab + p3lab

	//impulses in ?? system
	TVector3 pa(1, 1, 1);
	TVector3 pp1(1, 1, 1);
	TVector3 pp2(1, 1, 1);
	TVector3 psum(1, 1, 1);

	pa.SetMagThetaPhi(fP[atel][ahit], fTheta[atel][ahit], fPhi[atel][ahit]);
	pp1.SetMagThetaPhi(fP[ptel[0]][phit[0]], fTheta[ptel[0]][phit[0]],
			fPhi[ptel[0]][phit[0]]);
	pp2.SetMagThetaPhi(fP[ptel[1]][phit[1]], fTheta[ptel[1]][phit[1]],
			fPhi[ptel[1]][phit[1]]);

	//alpha
	fALab.SetVect(pa);
	fALab.SetE(fT[atel][ahit] + kAMASS);
	fAThetaLab = pa.Theta();
	fAPhiLab = pa.Phi();
	fAPcLab = pa.Mag();
	//proton1
	fP1Lab.SetVect(pp1);
	fP1Lab.SetE(fT[ptel[0]][phit[0]] + kPMASS);
	fP1ThetaLab = pp1.Theta();
	fP1PhiLab = pp1.Phi();
	fP1PcLab = pp1.Mag();
	//proton2
	fP2Lab.SetVect(pp2);
	fP2Lab.SetE(fT[ptel[1]][phit[1]] + kPMASS);
	fP2ThetaLab = pp2.Theta();
	fP2PhiLab = pp2.Phi();
	fP2PcLab = pp2.Mag();

	psum = pa + pp1 + pp2;

	Elab = kAMASS + 2 * kPMASS + fT[atel][ahit] + fT[ptel[0]][phit[0]]
			+ fT[ptel[1]][phit[1]];

	PClab = psum.Mag();

	//6Be
	f6BeIM = TMath::Power((TMath::Power(Elab, 2) - TMath::Power(PClab, 2)), 0.5)
			- (kAMASS + 2 * kPMASS);
	f6BeThetaLab = psum.Theta();
	f6BePhiLab = psum.Phi();
	f6BePcLab = psum.Mag();
	f6BeLab.SetVect(psum);
	f6BeLab.SetE(Elab);

	//neutron
	const Double_t TLiA = fTBeamRec; //AMeV
	Double_t TLi = 6 * TLiA; //MeV

	Double_t pc_n2 = LawOfCosines(psum.Mag(),
			TMath::Power(PC2(TLi, kLiMASS), 0.5), psum.Theta());
	Double_t Tn = T(pc_n2, kNMASS);
	Double_t En = Tn + kNMASS;

	//neutron theta
	//	pc_be^2 = pc_li^2 + pc_n^2 - 2*pc_li*pc_n*cos(theta_n)
	//		==>
	//	cos(theta_n) = (pc_li^2 + pc_n^2 - pc_be^2)/2*pc_li*pc_n

	Double_t pc_li2 = PC2(TLi, kLiMASS);

	Double_t cos_theta_n = (pc_li2 + pc_n2 - Power(psum.Mag(), 2))
			/ (2 * Sqrt(pc_li2) * Sqrt(pc_n2));
	TVector3 pneutron;
	pneutron.SetMagThetaPhi(Sqrt(pc_n2), TMath::ACos(cos_theta_n),
			psum.Phi() + TMath::Pi());

	//neutron energy		done
	fNLab.SetVect(pneutron);
	fNLab.SetE(En);
	fNPcLab = pneutron.Mag();
	fNThetaLab = pneutron.Theta();
	fNPhiLab = pneutron.Phi();

	return;
}


void BeEvent::SetProductsLabNew() {
	fNOA = NOParticlesTauE(0, kALPHA) + NOParticlesTauE(1, kALPHA);

//	printf("%d\t%d\t%d\n", fNOA, NOParticlesTauE(0, kALPHA), NOParticlesTauE(1, kALPHA));
	fNOP = NOParticlesTauE(0, kPROTON) + NOParticlesTauE(1, kPROTON);

//	cout << fNOtHits[0] + fNOtHits[1] << endl;

//	if (fNOtHits[0] + fNOtHits[1] != 3) {
////	if (fNOHits[0] + fNOHits[1] != 3) {
//		return;
//	}

//	if (NOParticlesTauE(0, kALPHA)+NOParticlesTauE(1, kALPHA) != 1) { return; }
//	if (NOParticlesTauE(0, kPROTON)+NOParticlesTauE(1, kPROTON) != 2) { return; }
	if (fNOA != 1) {
		return;
	}
	if (fNOP != 2) {
		return;
	}

//	fTrigger = 1;

	//loop over all hits to determine number of alphas and protons
//	fce ParticleNumberTau()

	Int_t atel = 0;
	Int_t ahit = 0;

	Int_t ptel[2] = { 0, 0 };
	Int_t phit[2] = { 0, 0 };
	Int_t nop = 0; //counter of protons
	Int_t noa = 0; //counter of alphas

	//loop over all telescopes
	for (Int_t i = 0; i < NOTEL; i++) {
		//loop over all hits in each telescope
		for (Int_t j = 0; j < fNOHits[i]; j++) {
			if (fTauSN[i][j] == -1) {
				continue;
			}
			if (fID[i][j] == kPROTON) {
				//zapamatovat si ho jako proton
				ptel[nop] = i;
				phit[nop] = j;
				nop++;
			}
			if (fID[i][j] == kALPHA) {
				//zapamatovat si ho jako alfu
				atel = i;
				ahit = j;
				noa++;
			}
		} //for over all hits in each telescope
	} //for over all telescopes

	//sum(m_i*c^2) + Tcm = sqrt( ( sum(m_c*c2 + T_i) )^2 - (p_system*c)^2 )

	Double_t Elab = 0; //Elab = sum(m_c*c2 + T_i)
	Double_t PClab = 0; //PClab = p_system*c, where p_system = p1lab + p2lab + p3lab

	//impulses in ?? system
	TVector3 pa(1, 1, 1);
	TVector3 pp1(1, 1, 1);
	TVector3 pp2(1, 1, 1);
	TVector3 psum(1, 1, 1);

	pa.SetMagThetaPhi(fP[atel][ahit], fTheta[atel][ahit], fPhi[atel][ahit]);
	pp1.SetMagThetaPhi(fP[ptel[0]][phit[0]], fTheta[ptel[0]][phit[0]],
			fPhi[ptel[0]][phit[0]]);
	pp2.SetMagThetaPhi(fP[ptel[1]][phit[1]], fTheta[ptel[1]][phit[1]],
			fPhi[ptel[1]][phit[1]]);

	//alpha
	fALab.SetVect(pa);
	fALab.SetE(fT[atel][ahit] + kAMASS);
	fAThetaLab = pa.Theta();
	fAPhiLab = pa.Phi();
	fAPcLab = pa.Mag();
	//proton1
	fP1Lab.SetVect(pp1);
	fP1Lab.SetE(fT[ptel[0]][phit[0]] + kPMASS);
	fP1ThetaLab = pp1.Theta();
	fP1PhiLab = pp1.Phi();
	fP1PcLab = pp1.Mag();
	//proton2
	fP2Lab.SetVect(pp2);
	fP2Lab.SetE(fT[ptel[1]][phit[1]] + kPMASS);
	fP2ThetaLab = pp2.Theta();
	fP2PhiLab = pp2.Phi();
	fP2PcLab = pp2.Mag();

	psum = pa + pp1 + pp2;

	Elab = kAMASS + 2 * kPMASS + fT[atel][ahit] + fT[ptel[0]][phit[0]]
			+ fT[ptel[1]][phit[1]];

	PClab = psum.Mag();

	//6Be
	f6BeIM = TMath::Power((TMath::Power(Elab, 2) - TMath::Power(PClab, 2)), 0.5)
			- (kAMASS + 2 * kPMASS);
	f6BeThetaLab = psum.Theta();
	f6BePhiLab = psum.Phi();
	f6BePcLab = psum.Mag();
	f6BeLab.SetVect(psum);
	f6BeLab.SetE(Elab);

	//neutron
	const Double_t TLiA = fTBeamRec; //AMeV
	Double_t TLi = 6 * TLiA; //MeV

	Double_t pc_n2 = LawOfCosines(psum.Mag(),
			TMath::Power(PC2(TLi, kLiMASS), 0.5), psum.Theta());
	Double_t Tn = T(pc_n2, kNMASS);
	Double_t En = Tn + kNMASS;

	//neutron theta
	//	pc_be^2 = pc_li^2 + pc_n^2 - 2*pc_li*pc_n*cos(theta_n)
	//		==>
	//	cos(theta_n) = (pc_li^2 + pc_n^2 - pc_be^2)/2*pc_li*pc_n

	Double_t pc_li2 = PC2(TLi, kLiMASS);

	Double_t cos_theta_n = (pc_li2 + pc_n2 - Power(psum.Mag(), 2))
			/ (2 * Sqrt(pc_li2) * Sqrt(pc_n2));
	TVector3 pneutron;
	pneutron.SetMagThetaPhi(Sqrt(pc_n2), TMath::ACos(cos_theta_n),
			psum.Phi() + TMath::Pi());

	//neutron energy		done
	fNLab.SetVect(pneutron);
	fNLab.SetE(En);
	fNPcLab = pneutron.Mag();
	fNThetaLab = pneutron.Theta();
	fNPhiLab = pneutron.Phi();

	return;
}

void BeEvent::SetQLiP() {
	if (f6BeIM == -100.) {
		return;
	}

	const Double_t TLi = 6 * fTBeamRec;

	const Double_t EBe = f6BeLab.E();
	const Double_t En = fNLab.E();

	//6Li+p-->6Be+n
	fQLiP = ReactionQ(TLi + kLiMASS, kPMASS, EBe, En);

	return;
}

void BeEvent::SetProductsCM1() {
	//	binary reaction;
	//
	//////////////////////////////////////
	////								//
	////				 6Be			//
	////				/				//
	////			   /				//
	////	          /theta_lab		//
	//// 6Li-------->X----------		//
	////             p\					//
	////			   \				//
	////			    n				//
	////								//
	//////////////////////////////////////
	//
	//beta and gamma for Lor. trans.
	//
	//	E_lab = T_lab + m*c^2 = m*c^2/( sqrt(1-(beta_t)^2) ) and 	E_lab^2 - (p_lab*c)^2 = (m*c^2)^2
	//				lab.			CM
	//		==>
	//	E_lab^2 - (p_lab*c)^2 = E_lab^2*(1-(beta_t)^2)
	//		==>
	//	beta_t = abs(p_lab)*c)/E_lab
	//	gamma_t	= 1/sqrt(1-(beta_t)^2)
	//
	//	E_lab = m_Li*c^2 + T_labLi + m_p*c^2
	//	p_lab = p_Li ==> (p_Li*c)^2 = T_Li^2 + 2*T_Li*m_Li*c^2

	if (f6BeIM == -100.) {
		return;
	}

	const Double_t T_Li = 6 * fTBeamRec;
	const Double_t E_lab = (kLiMASS + T_Li + kPMASS);

	const Double_t beta = Sqrt(PC2(T_Li, kLiMASS)) / E_lab; //p_Li_lab/E_lab

	//lorentz transformation lab-->CM

	//6Be
	//px_lab --> px_CM
	//py_lab --> py_CM
	//pz_lab --> pz_CM
	// from (px, py, pz)_cm ==> theta_CM
	TLorentzVector pcBelab(f6BeLab);
	TLorentzRotation S(0., 0., -beta);

	f6BeCM1 = pcBelab.Transform(S);
	f6BeThetaCM1 = f6BeCM1.Theta();
	f6BePhiCM1 = f6BeCM1.Phi();
	f6BePcCM1 = f6BeCM1.P();

	return;
}

void BeEvent::SetProductsCM() {
	if (f6BeIM == -100.) {
		return;
	}

	//rotate coordinate system, new axis Z is consintent with
	// the direction of neutron in laboratory system

	TLorentzVector vlBe = f6BeLab; //vl means Lorentz vector
	TLorentzVector vlNeutron = fNLab;
	TLorentzVector vlAlpha = fALab;
	TLorentzVector vlP1 = fP1Lab;
	TLorentzVector vlP2 = fP2Lab;

	TRotation r1;
	r1.SetZAxis(-vlNeutron.Vect(), vlBe.Vect()); //opacny smer vyletu neutronu
	TLorentzRotation rot(r1);

	vlBe.Transform(rot.Inverse()); //to obtain vector beta
	vlNeutron.Transform(rot.Inverse());
	vlAlpha.Transform(rot.Inverse());
	vlP1.Transform(rot.Inverse());
	vlP2.Transform(rot.Inverse());

	//transformation from lab system to CM of 6Be
	TVector3 beta = vlBe.BoostVector();

	vlBe.Boost(-beta);
	vlNeutron.Boost(-beta);
	vlAlpha.Boost(-beta);
	vlP1.Boost(-beta);
	vlP2.Boost(-beta);


	//protons mixing
	TLorentzVector vlPa;
	TLorentzVector vlPb;

	if (ranChoice.Uniform() < 0.5) {
		vlPa = vlP1;
		vlPb = vlP2;
	} else {
		vlPa = vlP2;
		vlPb = vlP1;
	}


	f6BeCM = vlBe;
	fNCM = vlNeutron;
	fNThetaCM = vlNeutron.Theta();
	fNPhiCM = vlNeutron.Phi();
	fACM = vlAlpha;
	fAThetaCM = vlAlpha.Theta();
	fAPhiCM = vlAlpha.Phi();
	fP1CM = vlPa;
	fP1ThetaCM = vlPa.Theta();
	fP1PhiCM = vlPa.Phi();
	fP2CM = vlPb;
	fP2ThetaCM = vlPb.Theta();
	fP2PhiCM = vlPb.Phi();

	//relative kinetical energy of two protons
	TLorentzVector vlPP = vlPa + vlPb;
	TVector3 betaPP = vlPP.BoostVector();
	vlPP.Boost(-betaPP);

	////////////////////////
	//	Jacobi systems
	////////////////////////


	//"T" system
	//	(p1+p2)+alpha
	const Double_t MTx = ReducedMass(kPMASS, kPMASS);
	TVector3 kTx = (vlPa.Vect() - vlPb.Vect()) * 0.5;

	const Double_t MTy = ReducedMass(kPMASS + kPMASS, kAMASS);
	TVector3 kTy = (4. * (vlPa.Vect() + vlPb.Vect()) - 2. * vlAlpha.Vect())
			* (1. / 6.);

	fTpp = kTx.Mag2() / (2 * MTx);
	fTapp = kTy.Mag2() / (2 * MTy);
	fCosThetaTk = kTx.Dot(kTy) / (kTx.Mag() * kTy.Mag());

	//	Jacobi "Y" system
	//	(p1+alpha)+p2
	const Double_t MYx = ReducedMass(kPMASS, kAMASS);
	TVector3 kYx = (4. * vlPa.Vect() - vlAlpha.Vect()) * 0.2;

	const Double_t MYy = ReducedMass(kPMASS + kAMASS, kPMASS);
	TVector3 kYy = ((vlPa.Vect() + vlAlpha.Vect()) - 5. * vlPb.Vect())
			* (1. / 6.);

	fTap = kYx.Mag2() / (2 * MYx);
	fTpap = kYy.Mag2() / (2 * MYy);
	fCosThetaYk = kYx.Dot(kYy) / (kYx.Mag() * kYy.Mag());

	return;
}

Double_t BeEvent::ReducedMass(Double_t m1, Double_t m2) {
	return m1 * m2 / (m1 + m2);
}

void BeEvent::ReadCuts(/*const char* cutfile*/) {
	//read graphical cuts from file
	//
	//cuts must have following names:
	//protons in 1st telescope:	proton1Cut
	//alphas in 1st telescope:	alpha1Cut
	//protons in 2nd telescope:	proton2Cut
	//alphas in 2nd telescope:	alpha2Cut

//	const char *cutfile = "cuts.root";

	Info("BeEvent::ReadCuts", "Name of cutfile is: %s", fCutFile.Data());

	if (fCutFile.Length() == 0) {
		Info("BeEvent::ReadCuts", "File containing cuts was not set");
		Info("BeEvent::ReadCuts",
				"Default method to identify particles will be used");
		return;
	}

	fCuts = new TFile(fCutFile, "READ");
	if (fCuts == NULL || fCuts->IsZombie()) {
		Error("BeEvent::ReadCuts",
				"File %s containing graphical cuts was not opened",
				fCutFile.Data());
		fCuts = NULL;
		Warning("BeEvent::ReadCuts",
				"Default method to identify particles will be used");
		return;
	}
	Info("BeEvent::ReadCuts", "File %s containing graphical cuts was opened",
			fCutFile.Data());

	cutAlpha1 = (TCutG*) fCuts->Get("alpha1Cut");
//	TCutG *calpha1 = (TCutG*) fCuts->Get("alpha1Cut");
//	cutAlpha1 = new TCutG( *calpha1 );
//	cutAlpha1 = new TCutG( *(TCutG*)fCuts->Get("alpha1Cut") );		//fixme !!! redo on copys
	if (cutAlpha1 == NULL) {
		Error("BeEvent::ReadCuts", "cutAlpha1 was not loaded from file %s",
				fCutFile.Data());
	}
	Info("BeEvent::ReadCuts", "TCutG object \"%s\" was loaded as cutAlpha1",
			cutAlpha1->GetName());

	cutProton1 = (TCutG*) fCuts->Get("proton1Cut");
	if (cutProton1 == NULL) {
		Error("BeEvent::ReadCuts", "cutProton1 was not loaded from file %s",
				fCutFile.Data());
	}
	Info("BeEvent::ReadCuts", "TCutG object \"%s\" was loaded as cutProton1",
			cutProton1->GetName());

	cutAlpha2 = (TCutG*) fCuts->Get("alpha2Cut");
	if (cutAlpha2 == NULL) {
		Error("BeEvent::ReadCuts", "cutAlpha2 was not loaded from file %s",
				fCutFile.Data());
	}
	Info("BeEvent::ReadCuts", "TCutG object \"%s\" was loaded as cutAlpha2",
			cutAlpha2->GetName());

	cutProton2 = (TCutG*) fCuts->Get("proton2Cut");
	if (cutProton2 == NULL) {
		Error("BeEvent::ReadCuts", "cutProton2 was not loaded from file %s",
				fCutFile.Data());
	}
	Info("BeEvent::ReadCuts", "TCutG object \"%s\" was loaded as cutProton2",
			cutProton2->GetName());

	fCuts->Close(); //todo potential hazard
	gROOT->cd();

//	cutAlpha1->Print();
//	cout << endl;
//	cutProton1->Print();
//	cout << endl;
//	cutAlpha2->Print();
//	cout << endl;
//	cutProton2->Print();

	return;
}

void BeEvent::CreateSiELosses() {

	//verification of the method
	mSiAlphaV.SetEL(SI_NELEMENTS, SI_RHO);
	mSiAlphaV.AddEL(SI_1_Z, SI_1_A, SI_1_W);
	mSiAlphaV.SetZP(2., 4.); //alpha
	mSiAlphaV.SetEtab(100000, 200.); //deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo
	//	mSiAlphaV.SetDeltaEtab(292);		//mean thickness for all fTheta1 is 292 micron for 280 micron detector
	//	mSiAlphaV.SetDeltaEtab(276);		//mean thickness for all fTheta2 is 276 micron for 275 micron detector

	Info("BeEvent::BeEvent", "Si TELosses");

	//CsI(Tl) calibration
	mSiAlpha.SetEL(SI_NELEMENTS, SI_RHO);
	mSiAlpha.AddEL(SI_1_Z, SI_1_A, SI_1_W);
	mSiAlpha.SetZP(2., 4.); //alpha
	mSiAlpha.SetEtab(100000, 200.); //deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo
	//	mSiAlpha.SetDeltaEtab(887);		//mean thickness for all fTheta1 is 887 micron for 850 micron detector
	//	mSiAlpha.SetDeltaEtab(804);		//mean thickness for all fTheta2 is 804 micron for 800 micron detector

	mSiP.SetEL(SI_NELEMENTS, SI_RHO);
	mSiP.AddEL(SI_1_Z, SI_1_A, SI_1_W);
	mSiP.SetZP(1., 1.); //protons
	mSiP.SetEtab(120000, 200.); //deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo
	//	mSiP.SetDeltaEtab(887);		//mean thickness for all fTheta1 is 887 micron for 850 micron detector
	//	mSiP.SetDeltaEtab(804);		//mean thickness for all fTheta2 is 804 micron for 800 micron detector

	Info("BeEvent::CreateSiELosses", "Si TELosses created");

	return;
}

void BeEvent::CreateCsIELosses() {
	//energy looses in CsI(Tl)
	//alpha in CsI
	//	mCsIAlpha = new TELoss();

	mCsIAlpha.SetEL(CSI_NELEMENTS, CSI_RHO);
	mCsIAlpha.AddEL(CSI_1_Z, CSI_1_A, CSI_1_W);
	mCsIAlpha.AddEL(CSI_2_Z, CSI_2_A, CSI_2_W);
	mCsIAlpha.SetZP(2., 4.); //alphas
	//100000 je urcite v poradku, 80000 staci, 50000 malo
	mCsIAlpha.SetEtab(100000, 200.); //deltaE
	//proton in CsI
	//	fCsIP = new TELoss();
	mCsIP.SetEL(CSI_NELEMENTS, CSI_RHO);
	mCsIP.AddEL(CSI_1_Z, CSI_1_A, CSI_1_W);
	mCsIP.AddEL(CSI_2_Z, CSI_2_A, CSI_2_W);
	mCsIP.SetZP(1., 1.); //protons
	mCsIP.SetEtab(100000, 200.); //deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

	Info("BeEvent::CreateCsIELosses", "CsI TELosses created");

	return;
}

void BeEvent::CreateTargetELosses() {

	const Double_t targetRho = T_NORMAL_RHO * (273 / T_TEMPERATURE)
			* (T_PRESSURE / 1);

	mH_P.SetEL(T_NELEMENTS, targetRho);
	mH_P.AddEL(T_1_Z, T_1_A, T_1_W);
	mH_P.SetZP(1., 1.); //protons
	//	mH_P.SetEtab(100000, 150.);	//deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo
	mH_P.SetEtab(100000, 200.); //deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

	mH_Alpha.SetEL(T_NELEMENTS, targetRho);
	mH_Alpha.AddEL(T_1_Z, T_1_A, T_1_W);
	mH_Alpha.SetZP(2., 4.); //alphas
	mH_Alpha.SetEtab(100000, 200.); //deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

	mWin_P.SetEL(T_WIN_NELEMENTS, T_WIN_RHO);
	mWin_P.AddEL(T_WIN_1_Z, T_WIN_1_A, T_WIN_1_W);
	mWin_P.AddEL(T_WIN_2_Z, T_WIN_2_A, T_WIN_2_W);
	mWin_P.AddEL(T_WIN_3_Z, T_WIN_3_A, T_WIN_3_W);
	mWin_P.SetZP(1., 1.); //protons
	mWin_P.SetEtab(100000, 200.); //deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

	mWin_Alpha.SetEL(T_WIN_NELEMENTS, T_WIN_RHO);
	mWin_Alpha.AddEL(T_WIN_1_Z, T_WIN_1_A, T_WIN_1_W);
	mWin_Alpha.AddEL(T_WIN_2_Z, T_WIN_2_A, T_WIN_2_W);
	mWin_Alpha.AddEL(T_WIN_3_Z, T_WIN_3_A, T_WIN_3_W);
	mWin_Alpha.SetZP(2., 4.); //alphas
	mWin_Alpha.SetEtab(100000, 200.); //deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

	Info("BeEvent::CreateTargetELosses", "Target TELosses created");

	return;
}

Bool_t BeEvent::WriteCondition() {
	if ( /*f6BeIM != -100.*/
	f6BeIM > -50. //-100. is the indicator of no 6Be
//			&& fP1Lab.fE-938<40 && fP2Lab.fE-938<40		//protons penetrating CsI detector or deuterons
			//	( /*(fT11.fMultESec > 0) && (fT11.fMultERing > 0) &&*/ (fT12.fMultESec>0) && (fTrigger==2) )
			//	||
			//	( (fT21.fMultESec > 0) && (fT21.fMultERing > 0) && (fT22.fMultESec > 0) )
			) {
		return kTRUE;
	}

	return kFALSE;
}

Bool_t BeEvent::WriteConditionPure() {
	//write condition for pure Be events
	//it will be used in in function BeWork::FillBeExpFile(...)
	//this file will be the file which will be analysed

	if ( //f6BeIM != -100.
//	f6BeIM > -50. && fP1Lab.E() - 938 < 40 && fP2Lab.E() - 938 < 40 //protons penetrating CsI detector or deuterons
	f6BeIM > -50. && fP1Lab.E() - 938 < fProtonThr && fP2Lab.E() - 938 < fProtonThr //protons penetrating CsI detector or deuterons
	&& f6BeIM > 0 && f6BeIM < 20 //only events in interesting energy range
			//	( /*(fT11.fMultESec > 0) && (fT11.fMultERing > 0) &&*/ (fT12.fMultESec>0) && (fTrigger==2) )
			//	||
			//	( (fT21.fMultESec > 0) && (fT21.fMultERing > 0) && (fT22.fMultESec > 0) )
			) {
		return kTRUE;
	}

	return kFALSE;
}