BeWork.cpp

/*
 * BeWork.cpp
 *
 *  Created on: 2.2.2010
 *      Author: Vratislav
 */

//unsigned int max_number_threads = 2;

#include "BeWork.h"
//#include "BeThread.h"

ClassImp(BeWork);

#ifndef __CINT__
//boost::mutex gFillMutex;
//boost::mutex gFillMutex;
#endif  /* __CINT __ */

BeWork::BeWork()
{
	Info("BeWork::BeWork", "CsI resolution is %f for alphas and %f for protons", 100*fCsIResA, 100*fCsIResP);
	Info("BeWork::BeWork", "Si resolution is %f keV", fSiRes*1000);

	Info("BeWork::BeWork", "Creating TELoss objects ...");
	CreateTELosses();

	Info("BeWork", "End of second constructor\n\n");

}

BeWork::BeWork(const char* configfile)
{

	fConfigFile = configfile;
	ReadConfigFile();
	SetWorkDir();
	Info("BeWork::BeWork", "Work directory: \"%s\"", fWorkDir.Data());

	Info("BeWork::BeWork", "CsI resolution is %f for alphas and %f for protons", 100*fCsIResA, 100*fCsIResP);
	Info("BeWork::BeWork", "Si resolution is %f keV", fSiRes*1000);

	Info("BeWork::BeWork", "Creating TELoss objects ...");
	CreateTELosses();

//	fCalibFile = "/data2/be/parameterfiles";
	Info("BeWork", "End of constructor\n\n");
}

BeWork::BeWork(BeWork &bework) {
	//copy constructor
	fSiAlpha = new TELoss( *(bework.fSiAlpha) );
	fSiP = new TELoss( *(bework.fSiP) );

	fCsIAlpha = new TELoss( *(bework.fCsIAlpha) );
	fCsIP = new TELoss( *(bework.fCsIP) );

	fTargetAlpha = new TELoss( *(bework.fTargetAlpha) );
	fTargetP = new TELoss( *(bework.fTargetP) );
	fTargetLi = new TELoss( *(bework.fTargetLi) );

	fTargetWinAlpha = new TELoss( *(bework.fTargetWinAlpha) );
	fTargetWinP = new TELoss( *(bework.fTargetWinP) );
	fTargetWinLi = new TELoss( *(bework.fTargetWinLi) );

	//configuration and parameter files
	fConfigFile = bework.fConfigFile;
	fWorkDir = bework.fWorkDir;
	fRawFilePath = bework.fRawFilePath;


	fSiRes = bework.fSiRes;
	fCsIResP = bework.fCsIResP;
	fCsIResA = bework.fCsIResA;

	fTBeamMC = bework.fTBeamMC;
	fTBeamResMC = bework.fTBeamResMC;
	fBeamX_MC = bework.fBeamX_MC;
	fBeamY_MC = bework.fBeamY_MC;
	fBeamX_sigma_MC = bework.fBeamX_sigma_MC;
	fBeamY_sigma_MC = bework.fBeamY_sigma_MC;
//
	fT1SimPosition = bework.fT1SimPosition;
	fT2SimPosition = bework.fT2SimPosition;

	fX11_FD = bework.fX11_FD;
	fX11 = bework.fX11;
	fX11_BD = bework.fX11_BD;

	fX12_FD = bework.fX12_FD;
	fX12 = bework.fX12;
	fX12_BD = bework.fX12_BD;

	fX13_FD = bework.fX13_FD;

	fX21_FD = bework.fX21_FD;
	fX21 = bework.fX21;
	fX21_BD = bework.fX21_BD;

	fX22_FD = bework.fX22_FD;
	fX22 = bework.fX22;
	fX22_BD = bework.fX22_BD;

	fX23_FD = bework.fX23_FD;

	Info("BeWork::BeWork", "Copy constructor finished");
}

BeWork::~BeWork()
{
	Info("BeWork::~BeWork", "Destructor called");
	delete fSiAlpha;
	delete fSiP;
	delete fCsIAlpha;
	delete fCsIP;
	delete fTargetAlpha;
	delete fTargetP;
	delete fTargetLi;
	delete fTargetWinAlpha;
	delete fTargetWinP;
	delete fTargetWinLi;
	Info("BeWork::~BeWork", "Destructor finished");
}

void BeWork::ReadConfigFile() {

	std::ifstream cfile;
	//generator file opening
	if (!fConfigFile.IsNull()) {
		cfile.open(fConfigFile.Data());
		if (!cfile.is_open()) {
			Error("BeWork::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.Append(" ");
		configuration.Append(line);
	}

//	cout << configuration.Data() << endl;

	//lower and upper case are important
	ConfigDictionary CDpath(configuration.Data());
	fRawFilePath = CDpath.GetString("rawFilePath");
//	fCutFile = CDpath.GetString("cutFile");

	//numbers only
	configuration.ToLower();
	ConfigDictionary CD(configuration.Data());
	fBeOnly = CD.GetBool("beonly");
	fsRatioMax = CD.GetDouble("sratiomax");
	fsRatioMin = CD.GetDouble("sratiomin");
//	printf("%d\t%d", fBeOnly, CD.GetBool("beonly"));
	fSiRes = CD.GetDouble("sires");
	fCsIResP = CD.GetDouble("csiresp");
	fCsIBestResP = CD.GetDouble("csirespmin");
	fCsIResA = CD.GetDouble("csiresa")/1.5;
	fCsIBestResA = CD.GetDouble("csiresamin");

	fTBeamMC = CD.GetDouble("tbeam");
	fTBeamResMC = CD.GetDouble("tbeamres");
	fBeamX_MC = CD.GetDouble("beamx");
	fBeamY_MC = CD.GetDouble("beamy");
	fBeamX_sigma_MC = CD.GetDouble("beamx_sigma");
	fBeamY_sigma_MC = CD.GetDouble("beamy_sigma");


	fT1SimPosition = CD.GetDouble("t1simposition");
	fT2SimPosition = CD.GetDouble("t2simposition");

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



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

	return;

}

void BeWork::SetWorkDir() {
	fWorkDir = fConfigFile;
	fWorkDir.Remove(fConfigFile.Last('/'), fConfigFile.Length());
}

const char* BeWork::GetWorkDir() {
	return fWorkDir.Data();
}

Int_t BeWork::FillExpFile(const char* inputrawfile, const char* outputfile,
		Long64_t noevents, Option_t *opt)
{
	//convert raw data in channel units into all needed kinematical units, new file containing tree "be" filled by "beevent" is created
	//
	//inputrawfile: .root file containing raw data
	//noevents: number of events to be filled, if noevents is set to 0, whole file will be processed
	//calibfilepath: path to files with calibration data "beSi.cal", "beCsIa.cal" and "beCsIp.cal"
	//outputfile: .root file with calibrated data
	//option: config: add config file name into the outputfile name

	TString ofile;
//	if (opt
	ofile.Form("%s/%s", fWorkDir.Data(), outputfile);

	TString ifile;
	ifile.Form("%s/%s", fRawFilePath.Data(), inputrawfile);

	//read calibration parameters and thresholds

	TFile *fr = new TFile(ifile.Data(), "READ");
	if (fr->IsOpen() == 0) {
		Warning("BeWork::FillExpFile", "File %s was not opened and won't be processed", ifile.Data());
		return 0;
	}


	//RAW data tree opening
	TTree *tr = (TTree*)fr->Get("RAW");
	if (!tr) {
		Warning("BeWork::FillExpFile", "Tree \"RAW\" was not found in file %s .", ifile.Data());
		return 0;
	}
	AculRaw *eventr = new AculRaw();
	Long64_t noEntries = tr->GetEntries(/*getCondition*/);
	if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
	tr->SetBranchAddress("channels", &eventr);

	//output file containing tree creating
	TFile *fw = new TFile(ofile.Data(), "RECREATE");	//dodelat overeni vytvoreni

	BeEvent *eventw = new BeEvent(fConfigFile.Data());
	TTree *tw = new TTree("be", "strom s energetickou a jinou informaci");		//predelat jmena, ...
	tw->Bronch("BePhys", "BeEvent", &eventw, 1024000, 99);		//large buffer needed for faster processing in the course of analysis


	Info("BeWork::FillExpFile", "In file %s %d events will be processed", ifile.Data(), (int)noEntries);
	if (fBeOnly) Info("BeWork::FillExpFile", "Only events containing Be will be saved");
	else Info("BeWork::FillExpFile", "All events will be saved");
	Info("BeWork::FillExpFile", "Output file %s will be filled", ofile.Data());

	fw->cd();

	for (Int_t i = 0; i <= (Int_t)noEntries; i++) {
		printProgBar(i, noEntries);
		fr->cd();
		tr->GetEntry(i);
		eventw->FillEvent(eventr);
		fw->cd();
		if (fBeOnly) {
			if (eventw->WriteCondition()) {
				tw->Fill();
			} //zkontrolovat
		} else
			tw->Fill();
	}//for

	printf("\n\n");

	fw->cd();
	tw->AutoSave();

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

	return (Int_t)noEntries;
}

Int_t BeWork::FillBeFile(const char* inputfile, Long64_t noevents, const char* rtree, const char* rbranch,
		const char* outputfile, const char* wtree, const char* wbranch, Option_t *opt) {
	//inputfile: .root file containing calibrated data
	//noevents: number of events to be filled, if noevents is set to 0, whole file will be processed
	//rtree: name of the tree containing information with structure BeEvent class to be read
	//rbranch: name of the branch to be read
	//outputfile: .root file with data containing Be event only, which will be filled into TTree named "beonly"
	//wtree: name of the branch to be wrote
	//wbranch: name of the tree containing information with structure BePureEvent class to be wrote
	//option:

	//output file name creating
	if (!strlen(outputfile)) {
		Error("BeWork::FillBeFile", "Name of outputfile must be set");
		return 0;
	}
	TString output(outputfile);



	TFile fr(inputfile, "READ");
	if (fr.IsOpen() == 0) {
		Warning("BeWork::FillBeFile", "File %s was not opened and won't be processed", inputfile);
		return 0;
	}
	//RAW data tree opening
	TTree *tr = (TTree*)fr.Get(rtree);
	if (!tr) {
		Warning("BeWork::FillBeFile", "Tree \"be\" was not found in file %s .", inputfile);
		return 0;
	}
	BeEvent *eventr = new BeEvent(fConfigFile.Data());
	Long64_t noEntries = tr->GetEntries(/*getCondition*/);
	if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
	tr->SetBranchAddress(rbranch, &eventr);

	//tree with simulated beam coordinates opening
	TTree *trbeam = 0;
	Double_t x = 0;
	Double_t y = 0;
	Double_t z = 0;
	Double_t thetaMC = 0;
	Double_t E_IM = 0;

	TLorentzVector *beamLab = new TLorentzVector();
	TLorentzVector *alphaLab = new TLorentzVector();
	TLorentzVector *p1Lab = new TLorentzVector();
	TLorentzVector *p2Lab = new TLorentzVector();

	TLorentzVector *alphaCM = new TLorentzVector();
	TLorentzVector *p1CM = new TLorentzVector();
	TLorentzVector *p2CM = new TLorentzVector();
	TLorentzVector *beCM = new TLorentzVector();

	Double_t sTpp;
	Double_t sTapp;
	Double_t sCosThetaTk;

	TString treereadname(rtree);
	if (treereadname.Contains("sim")) {
		Info("BeWork::FillBeFile", "Tree containing information about simulated beam opening");
		trbeam = (TTree*)fr.Get("sbeam");
		if (!trbeam) {
			Warning("BeWork::FillBeFile", "Tree \"sbeam\" was not found in file %s .", inputfile);
			return 0;
		}
		trbeam->SetBranchAddress("bx", &x);
		trbeam->SetBranchAddress("by", &y);
		trbeam->SetBranchAddress("bz", &z);
		trbeam->SetBranchAddress("thetaMC", &thetaMC);
		trbeam->SetBranchAddress("E_IM", &E_IM);

		trbeam->SetBranchAddress("beamLab.", &beamLab);
		trbeam->SetBranchAddress("sALab.", &alphaLab);
		trbeam->SetBranchAddress("sP1Lab.", &p1Lab);
		trbeam->SetBranchAddress("sP2Lab.", &p2Lab);

		trbeam->SetBranchAddress("sACM.", &alphaCM);
		trbeam->SetBranchAddress("sP1CM.", &p1CM);
		trbeam->SetBranchAddress("sP2CM.", &p2CM);
		trbeam->SetBranchAddress("sBeCM.", &beCM);

		trbeam->SetBranchAddress("sTpp", &sTpp);
		trbeam->SetBranchAddress("sTapp", &sTapp);
		trbeam->SetBranchAddress("sCosThetaTk", &sCosThetaTk);
		
	}//if



	//output file containing tree creating
	TFile fw(output.Data(), "RECREATE");	//dodelat overeni vytvoreni
	BePureEvent *eventw = new BePureEvent();
	TTree *tw = new TTree(wtree, "tree containing Be events only");		//predelat jmena, ...
	tw->Bronch(wbranch, "BePureEvent", &eventw, 1024000, 99);		//large buffer needed for faster processing in the course of analysis

	//tree containing beam information
	TTree *twbeam = 0;
	if (treereadname.Contains("sim")) {
		Info("BeWork::FillBeFile", "Tree containing information about simulated beam creating");
		twbeam = new TTree("sbeam", "tree with beam position and ThetaCM read from generator (not used for simulation)");

		twbeam->Branch("bx", &x, "bx/D");
		twbeam->Branch("by", &y, "by/D");
		twbeam->Branch("bz", &z, "bz/D");

		twbeam->Branch("thetaMC", &thetaMC, "thetaMC/D");
		twbeam->Branch("E_IM", &E_IM, "E_IM/D");

		twbeam->Bronch("beamLab", "TLorentzVector", &beamLab, 1024000, 99);
		twbeam->Bronch("sALab.", "TLorentzVector", &alphaLab, 1024000, 99);
		twbeam->Bronch("sP1Lab.", "TLorentzVector", &p1Lab, 1024000, 99);
		twbeam->Bronch("sP2Lab.", "TLorentzVector", &p2Lab, 1024000, 99);


		twbeam->Bronch("sACM.", "TLorentzVector", &alphaCM, 1024000, 99);
		twbeam->Bronch("sP1CM.", "TLorentzVector", &p1CM, 1024000, 99);
		twbeam->Bronch("sP2CM.", "TLorentzVector", &p2CM, 1024000, 99);
		twbeam->Bronch("sBeCM.", "TLorentzVector", &beCM, 1024000, 99);


		twbeam->Branch("sTpp", &sTpp, "sTpp/D");
		twbeam->Branch("sTapp", &sTapp, "sTapp/D");
		twbeam->Branch("sCosThetaTk", &sCosThetaTk, "sCosThetaTk/D");
	}


	Info("BeWork::FillBeFile", "In file %s %d events will be processed", inputfile, (int)noEntries);
	Info("BeWork::FillBeFile", "Output file %s will be filled", output.Data());

	for (Int_t i = 0; i < (Int_t)noEntries; i++) {
		printProgBar(i, (Int_t)noEntries);
		fr.cd();
		tr->GetEntry(i);
		if (trbeam) trbeam->GetEntry(i);
		if ( eventr->WriteConditionPure() ) {
			eventw->FillEvent(eventr);
			fw.cd();
			tw->Fill();
			if (twbeam) twbeam->Fill();
		}//if
	}//for
	Info("BeWork::FillBeFile", "%d events containing Be were found", (int)tw->GetEntries());
	printf("\n\n");


	fw.cd();
	tw->AutoSave();
	Info("BeWork::FillBeFile", "Tree \"%s\" with Be events saving", tw->GetName());
	tw->Write();
	if (twbeam) {
		Info("BeWork::FillBeFile", "Tree \"%s\" with beam information saving", twbeam->GetName());
		twbeam->Write();
	}
	delete eventw;
	fw.Close();

	fr.cd();
	delete eventr;
	fr.Close();



	return (Int_t)noEntries;
}

Int_t BeWork::FillBeExpFile(const char *inputfile, const char* outputfile, Long64_t noevents) {
	//fill file containing tree with BePureEvent class structure
	//inputfile:
	//noevents:
	//outputfile: if empty, character sequence "Exp" in inputfile will be changed to "Be"

	TString ifile;
	ifile.Form("%s/%s", fWorkDir.Data(), inputfile);

	TString ofile;
	ofile.Form("%s/%s", fWorkDir.Data(), outputfile);

//	TString output(inputfile);
//	if (strlen(outputfile)) {
//		output = outputfile;
//	}
//	else {
//		output.ReplaceAll("Exp", "Be");
//	}

	const char* rtreename = "be";
	const char* rbranchname = "BePhys";
	const char* wtreename = "beonly";
	const char* wbranchname = "BeEvents";

	return FillBeFile(ifile.Data(), noevents, rtreename, rbranchname, ofile.Data(), wtreename, wbranchname, "");
}

Int_t BeWork::FillBeSimFile(const char *inputfile, const char* outputfile, Long64_t noevents) {
	//inputfile:
	//noevents:
	//outputfile: if empty, character sequence "Exp" in inputfile will be changed to "Be"


	TString ifile = fWorkDir + '/' + inputfile;

	TString ofile = fWorkDir + '/' + outputfile;

	const char* rtreename = "simbe";
	const char* rbranchname = "BeSimPhys";
	const char* wtreename = "simbeonly";
	const char* wbranchname = "BeSimEvents";

	return FillBeFile(ifile.Data(), noevents, rtreename, rbranchname, ofile.Data(), wtreename, wbranchname, "");
}

void BeWork::MixSimBeFiles(const Int_t infiles, const char* treename, ...) {
	// obsolete function which probably won't be developed more
	//
	//dodelat podminku "soubor" && vstupy != 0
	//	if (!i) return 0;

	va_list params;
	va_start(params, treename);

	TString file[infiles];
	Int_t events[infiles];


	for (Int_t i = 0; i < infiles; i++) {
		file[i] = va_arg(params, char*);
		events[i] = va_arg(params, Int_t);
	}
	va_end(params);

	//open output mix file
	BePureEvent *revent = new BePureEvent();
	TFile fw("mixed.root", "RECREATE");
	if (!fw.IsOpen()) {
		Error("BeWork::MixSimBeFiles", "File %s was not opened", "mixed.root");
		return;
	}
	TTree *tw = new TTree(treename, "tree of mixed simulated events");
	//!!!! make possible to mix also BeEvent trees
	tw->Bronch("BeSimMix", "BePureEvent", &revent, 1024000, 99);

	for (Int_t i = 0; i < infiles; i++) {
//		cout << file[i] << endl;
		//file opening
		TFile fr(file[i].Data());
		if (!fr.IsOpen()) {
			Warning("BeWork::MixSimBeFiles", "File \"%s\" was not open", file[i].Data());
			continue;
		}
		TTree *tr = (TTree*)fr.Get(treename);
//		tr->GetListOfBranches()->Print();
		tr->SetBranchAddress("BeSimEvents", &revent);
		Long64_t noEntries = tr->GetEntries(/*getCondition*/);
		if (events[i] != 0 && noEntries > events[i]) { noEntries = events[i]; }
//		tr->GetEntries();
		Info("BeWork::MixSimBeFiles", "%d events from file %s will be processed", (Int_t)noEntries, file[i].Data());
		for (Long64_t j = 0; j < noEntries; j++) {
			BeWork::printProgBar((Int_t)j, (Int_t)noEntries);
			tr->GetEntry(j);
			tw->Fill();
		}//for j
		fr.Close();
		if (fr.IsOpen()) {
			Warning("BeWork::MixSimBeFiles", "File \"%s\" closing error", file[i].Data());
			return;
		}
	}//for i

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

	fw.Close();

	return;
}

TTree* BeWork::OpenTree(const char* inputfile, const char* treename,
		const Color_t color, const char* friendtreename)
{
//	open file with saved TTree named "treename"
//	and return pointer to this tree

	TFile *fr = new TFile(inputfile, "READ");
	if (fr->IsOpen() == 0) {
		Warning("FillFile", "File %s was not opened and won't be processed", inputfile);
		return 0;
	}

	TTree *tr = (TTree*)fr->Get(treename);
	if (!tr) {
		Warning("OpenTree", "Tree \"%s\" was not found in file %s", treename, inputfile);
		return 0;
	}

	//add friend
	TString fr_name = friendtreename;
	if ( !fr_name.IsNull() ) {
		tr->AddFriend(fr_name.Data(), fr);
	}

	tr->SetMarkerStyle(20);
	tr->SetLineColor(color);
	tr->SetMarkerColor(color);

	gROOT->cd();

	return tr;
}

TChain* BeWork::OpenChain(const char* inputfile, int first, int last,
		const char* treename, const Color_t color, const char* friendtreename, Option_t* option)
{
	//	open series of files with identical names numbered
	//	from first to last with saved TTree named
	//	"treename" and return pointer to chain constisting
	//	all of this trees

	TString	opt = option;
	opt.ToLower();

	TChain *ch = new TChain(treename);
	Char_t	fileName[100];	//potreba delsi pole ?232?


//	if (opt.Contains("gp")) {

	for (Int_t i = first; i <= last; i++) {
		sprintf(fileName, "%s%03d.root", inputfile, i);
		if (opt.Contains("verbose")) {
			cout << fileName << endl;
		}
		ch->Add(fileName);
	}

	//adding friend chain
	TString fr_name = friendtreename;
	if ( !fr_name.IsNull() ) {
		TChain *fr_ch = new TChain(fr_name.Data());
		TString filename;

		for (Int_t i = first; i <= last; i++) {
			filename.Form("%s%03d.root", inputfile, i);
//			sprintf(fileName, "%s%03d.root", inputfile, i);
//			cout << fileName << endl;
			if (opt.Contains("verbose")) {
				cout << filename << endl;
			}
//			ch->Add(fileName);
			fr_ch->Add(filename.Data());
		}

//		ch->AddFriend(fr_name);
		ch->AddFriend(fr_ch);
	}//if



	ch->SetMarkerStyle(20);
	ch->SetLineColor(color);
	ch->SetMarkerColor(color);

	gROOT->cd();

	return ch;
}

TGeoManager* BeWork::BuildGeometry()
{
	//geometry units are cm

	//dodelat nejakou kontrolu, jestli uz nahodou takova geometrie neexistuje

	TGeoManager *geometry = new TGeoManager("SETUP", "6Be structure experimental setup");

	//materials
	TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
//	static TGeoMaterial *matSi	= new TGeoMaterial("Si", 28.085, 14, 2.332);

	//media
	TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
//	static TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);

	TGeoVolume *top = geometry->MakeBox("Top", medVacuum, 100., 100., 100.);
	geometry->SetTopVolume(top);
	top->SetLineColor(kMagenta);

	//target
	TGeoVolume *target = MakeTarget();
	target->SetName("target");
	top->AddNode(target, 1);

	//first telescope
	const Double_t z11 = fT1SimPosition;			//in cm
	const Double_t z12 = z11 + 1.;		//10.1;	//in cm		//original value
//	const Double_t z12 = z11 + 0.65;		//10.1;	//in cm
//	const Double_t z13 = z11 + 3.5;		//12.6;	//in cm		//original value
	const Double_t z13 = z11 + 2.9;		//12.6;	//in cm
	TGeoVolume *annDSD1 = MakeAnnularDetector(64, 32, fX11, fX11_FD, fX11_BD);		//backdl (from previous detector) + frontdl = 1.77
	annDSD1->SetName("T11");
	top->AddNode(annDSD1, 1, new TGeoTranslation(0, 0, z11));
	TGeoVolume *annSSD1 = MakeAnnularDetector(64, 1, fX12, fX12_FD, fX12_BD);		//backdl (from previous detector) + frontdl = 1.7 + 2.36 = 4.06 == XD12
	annSSD1->SetName("T12");
	top->AddNode(annSSD1, 1, new TGeoTranslation(0, 0, z12));
	TGeoVolume *annCsI1 = MakeCsIDetector();								//backdl (from previous detector) + frontdl = 2. + ?? = ??? ==
//	TGeoVolume *annCsI1 = MakeCsIDetectorMS("T13");
	annCsI1->SetName("T13");
	top->AddNode(annCsI1, 1, new TGeoTranslation(0., 0., z13));
//*/
	//second telescope
	const Double_t z21 = fT2SimPosition;			//in cm
	const Double_t z22 = z21 + 1.;		//31.1;	//in cm		//original value
//	const Double_t z22 = z21 + 0.65;		//31.1;	//in cm
//	const Double_t z23 = z21 + 3.5;		//33.6;	//in cm		//original value
	const Double_t z23 = z21 + 2.9;		//33.6;	//in cm

	TGeoVolume *annDSD2 = MakeAnnularDetector(64, 32, fX21, fX21_FD, fX21_BD);	//!!!!!! thickness of T21 was changed !!!!!!
	annDSD2->SetName("T21");
	top->AddNode(annDSD2, 1, new TGeoTranslation(0, 0, z21));
	TGeoVolume *annSSD2 = MakeAnnularDetector(64, 1, fX22, fX22_FD, fX22_BD);
	annSSD2->SetName("T22");
	top->AddNode(annSSD2, 1, new TGeoTranslation(0, 0, z22));
	TGeoVolume *annCsI2 = MakeCsIDetector();
//	TGeoVolume *annCsI2 = MakeCsIDetectorMS("T23");
	annCsI2->SetName("T23");
	top->AddNode(annCsI2, 1, new TGeoTranslation(0., 0., z23));


	geometry->CloseGeometry();

	Info("BeWork::BuildGeometry", "Distances in virtual setup: zT11 = %4.1f cm zT21 = %4.1f cm", fT1SimPosition, fT2SimPosition);

	return geometry;
}

TGeoVolume* BeWork::MakeAnnularDetector(Int_t nosecs, Int_t norings,
		Double_t thickness, Double_t fdeadlayer, Double_t bdeadlayer)
{
	//nosecs:	number of sectors
	//norings: number of rings
	//thickness: detector thickness in microns
	//fdeadlayer: frontside deadlayer thickness in microns
	//bdeadlayer: backside deadlayer thickness in microns



	//annular detector
	const Double_t	siThickness = thickness*MicToCm();
	const Double_t	frontdl = fdeadlayer*MicToCm();
	const Double_t	backdl = bdeadlayer*MicToCm();
	const Double_t	sSiRmin = 1.6;		//r_min of sensitive region in cm
	const Double_t	sSiRmax = 4.2;		//r_max of sensitive region in cm
	const Double_t	siRmin = 1.4;		//r_min of Si plate in cm
//	const Double_t	siRmin = 1.6;		//r_min of Si plate in cm, version of MS
	Double_t	deskThickness = 0.38;
	if (siThickness > deskThickness) { deskThickness = siThickness; }
//	const Double_t	dist = 0.01;		//in cm; distance between two neighbourgh rings or sectors			//!!!!!!!!!! original value
	const Double_t	dist = 0.00001;		//in cm; distance between two neighbourgh rings or sectors			//!!!!!!!!!! almost zero
//	const Double_t	dist = 1.;		//distance between two neighbourgh rings or sectors

	//materials
	TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
	/*static*/ TGeoMaterial *matSi	= new TGeoMaterial("Si", 28.085, 14, 2.332);

	//media
	TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
	/*static*/ TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);

//	cout << endl << endl;
//	printf("\t%3.3f\t%3.3f\n", deskThickness, siThickness);
//	cout << endl << endl;

	//detector box
//	TGeoBBox *detbox = new TGeoBBox(sSiRmax + 1., sSiRmax + 1., deskThickness);
	TGeoBBox *detbox = new TGeoBBox(sSiRmax + 1., sSiRmax + 1., siThickness/2.+0.1);
	TGeoVolume *annSSD = new TGeoVolume("annSSD", detbox, medVacuum);

	TGeoTube *senslayer = new TGeoTube(siRmin, sSiRmax, siThickness/2.);
	TGeoVolume *siDisc = new TGeoVolume("discSi", senslayer,  medSi);
	siDisc->SetLineColor(3);
	annSSD->AddNode(siDisc, 1);

	/////////////
	////rings
	/////////////
	TGeoTube *ringtube = new TGeoTube(sSiRmin, sSiRmax, backdl/2.);
	TGeoVolume *rlayer = new TGeoVolume("rlayer", ringtube, medVacuum);

	Double_t rmin = 0, rmax = 0;

	for (Int_t i = 0; i < norings; i++) {
		rmin = sSiRmin + i*(sSiRmax-sSiRmin)/norings;
		rmax = sSiRmin + (i+1)*(sSiRmax-sSiRmin)/norings - dist;
		TGeoTube *shapering = new TGeoTube(rmin, rmax, backdl/2.);
		TGeoVolume *ring = new TGeoVolume("ring", shapering, medSi);
		ring->SetLineColor(2);
		rlayer->AddNode(ring, i);
	}

	annSSD->AddNode(rlayer, 1, new TGeoTranslation(0, 0, siThickness/2.+backdl/2.));

	////////////////
	//// sectors
	////////////////
	const Double_t tubsphi = 180./nosecs;

//	cout << tubsphi << endl;

	TGeoTube *sectube = new TGeoTube(sSiRmin, sSiRmax, frontdl/2.);
	TGeoVolume *slayer = new TGeoVolume("slyer", sectube, medVacuum);

	//sensitive strip
	TGeoTubeSeg *tubesec = new TGeoTubeSeg(sSiRmin, sSiRmax, frontdl/2., -tubsphi, tubsphi);
	tubesec->SetName("tubesec");
	TGeoBBox *space = new TGeoBBox(sSiRmax, dist/2., frontdl);
	space->SetName("space");

	TGeoUnion *sub1 = new TGeoUnion(space, space, new TGeoRotation("r1", tubsphi, 0, 0), new TGeoRotation("r2", -tubsphi, 0, 0));
	TGeoCompositeShape *cs1 = new TGeoCompositeShape("cs1", sub1);
	TGeoSubtraction *sub2 = new TGeoSubtraction(tubesec, cs1, 0, 0);
	TGeoCompositeShape *ssec = new TGeoCompositeShape("ssec", sub2);
	TGeoVolume *senssec = new TGeoVolume("senssec", ssec, medSi);
	senssec->SetLineColor(2);

	for (Int_t i = 0; i < nosecs; i++) {
		slayer->AddNode(senssec, i, new TGeoRotation("secrot", 2*i*tubsphi+tubsphi, 0, 0));
	}

	TGeoRotation platerot("platerot", 90, 180, 0);
	TGeoTranslation platetrans("platetrans", 0, 0, -(siThickness/2.+frontdl/2.));
	annSSD->AddNode(slayer, 1, new TGeoCombiTrans(platetrans, platerot) );

	return annSSD;
}

TGeoVolume* BeWork::MakeCsIDetector()
{
	//all sizes are in cm

	//one CsI sector parameters
	const Double_t segX1 = 1.8;		//size of outer edge
	const Double_t segX2 = 0.6;		//size of inner edge
	const Double_t segY = 1.9;		//thickness
	const Double_t segZ = 3.0;		//length from the center to the edge
	const Double_t rIn = 3.75;		//diameter of the holel; radius = 1.875 cm
	const Double_t phisec = 180./16.;
//	const Double_t fdl = 12.;

	//materials
	TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
	TGeoMaterial *matSi	= new TGeoMaterial("Si", 28.085, 14, 2.332);

	//media
	TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
	TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);

//	cout << endl << endl;
//	printf("\t%3.3f\t%3.3f\t%3.3f\n", rIn + segZ + 1., segY/2.+0.2, rIn/2. + segZ/2.);
//	cout << endl << endl;

	TGeoBBox *detbox = new TGeoBBox(rIn + segZ + 1., rIn + segZ + 1., segY/2.+0.2);
	TGeoVolume *annDet = new TGeoVolume("annDet", detbox, medVacuum);

	TGeoTrd1 *csI = new TGeoTrd1(segX1/2., segX2/2., segY/2., segZ/2.);
	TGeoVolume *csSec = new TGeoVolume("csSec", csI, medSi);
	csSec->SetLineColor(kGray+1);

	Double_t phi = 0;



	for (Int_t i = 0; i < 16; i++) {
		phi = 2*i*phisec+phisec;
		TGeoRotation platerot("platerot", -phi, 90., 0.);
		TGeoTranslation platetrans("platetrans", (rIn/2. + segZ/2.)*TMath::Sin(phi*TMath::DegToRad()), (rIn/2. + segZ/2.)*TMath::Cos(phi*TMath::DegToRad()), 0);
		annDet->AddNode(csSec, i, new TGeoCombiTrans(platetrans, platerot));
	}

	return annDet;
}

TGeoVolume* BeWork::MakeCsIDetectorMS(TString name)
{
	//all sizes are in cm
	const Int_t nosecs = 16;
	const Double_t	siThickness = 1.9;
//	const Double_t	frontdl = 4.;
//	const Double_t	backdl = 4.;
	const Double_t	sSiRmin = 1.6;		//r_min of sensitive region in cm
	const Double_t	sSiRmax = 4.2;		//r_max of sensitive region in cm
//	const Double_t	siRmin = 1.4;		//r_min of Si plate in cm
	Double_t	deskThickness = 0.38;
	if (siThickness > deskThickness) { deskThickness = siThickness; }
	//	const Double_t	dist = 0.01;		//in cm; distance between two neighbourgh rings or sectors			//!!!!!!!!!! original value
	const Double_t	dist = 0.0001;		//in cm; distance between two neighbourgh rings or sectors			//!!!!!!!!!! almost zero
	//	const Double_t	dist = 1.;		//distance between two neighbourgh rings or sectors


	//materials
	TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
	TGeoMaterial *matSi	= new TGeoMaterial("Si", 28.085, 14, 2.332);

	//media
	TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
	TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);

	TGeoBBox *detbox = new TGeoBBox(sSiRmax + 1., sSiRmax + 1., siThickness/2.+0.1);
	TGeoVolume *annDet = new TGeoVolume("annDet", detbox, medVacuum);


	////////////////
	//// sectors
	////////////////
	const Double_t tubsphi = 180./nosecs;

	//	cout << tubsphi << endl;

	TGeoTube *sectube = new TGeoTube(sSiRmin, sSiRmax, siThickness/2.);
	TGeoVolume *slayer = new TGeoVolume(name.Data(), sectube, medVacuum);

	//sensitive strip
	TGeoTubeSeg *tubesec = new TGeoTubeSeg(sSiRmin, sSiRmax, siThickness/2., -tubsphi, tubsphi);
	tubesec->SetName("tubesec");
	TGeoBBox *space = new TGeoBBox(sSiRmax, dist/2., siThickness);
	space->SetName("space");

	TGeoUnion *sub1 = new TGeoUnion(space, space, new TGeoRotation("r1", tubsphi, 0, 0), new TGeoRotation("r2", -tubsphi, 0, 0));
	TGeoCompositeShape *cs1 = new TGeoCompositeShape("cs1", sub1);
	TGeoSubtraction *sub2 = new TGeoSubtraction(tubesec, cs1, 0, 0);
	TGeoCompositeShape *ssec = new TGeoCompositeShape("ssec", sub2);
	TGeoVolume *csSec = new TGeoVolume("csSec", ssec, medSi);
	csSec->SetLineColor(kGray+1);

	for (Int_t i = 0; i < nosecs; i++) {
		slayer->AddNode(csSec, i, new TGeoRotation("secrot", 2*i*tubsphi+tubsphi, 0, 0));
	}

	TGeoRotation platerot("platerot", 90, 180, 0);
	TGeoTranslation platetrans("platetrans", 0., 0., 0.);
	annDet->AddNode(slayer, 1, new TGeoCombiTrans(platetrans, platerot) );






//	TGeoTrd1 *csI = new TGeoTrd1(segX1/2., segX2/2., segY/2., segZ/2.);
//	TGeoVolume *csSec = new TGeoVolume("csSec", csI, medSi);
//	csSec->SetLineColor(kGray+1);

//	Double_t phi = 0;



//	for (Int_t i = 0; i < 16; i++) {
//		phi = 2*i*phisec+phisec;
//		TGeoRotation platerot("platerot", -phi, 90., 0.);
//		TGeoTranslation platetrans("platetrans", (rIn/2. + segZ/2.)*TMath::Sin(phi*TMath::DegToRad()), (rIn/2. + segZ/2.)*TMath::Cos(phi*TMath::DegToRad()), 0);
//		annDet->AddNode(csSec, i, new TGeoCombiTrans(platetrans, platerot));
//	}

	return annDet;
}

TGeoVolume* BeWork::MakeTarget()
{
	//one CsI sector parameters
	const Double_t rmin = 0.;
	const Double_t rmax = 2.;
	const Double_t dz = (T_THICKNESS+T_WIN_THICK*2)/10000;	//from mcm to cm
	const Double_t d_mat = T_THICKNESS/10000;	//from mcm to cm

	//materials
	TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0., 0., 0.);
//	TGeoMaterial *matSi	= new TGeoMaterial("Si", 28.085, 14, 2.332);

	//media
	TGeoMedium *medVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
//	TGeoMedium *medSi = new TGeoMedium("Si", 2, matSi);

	//stainless steel box for target material
	TGeoTube *tartube = new TGeoTube(rmin, rmax, dz/2.);
	TGeoVolume *target = new TGeoVolume("target", tartube, medVacuum);
	target->SetLineColor(kGray);

	//target material
	TGeoTube *tar_mat_tube = new TGeoTube(rmin, rmax, d_mat/2.);
	TGeoVolume *tar_mat = new TGeoVolume("tarmaterial", tar_mat_tube, medVacuum);
	tar_mat->SetLineColor(kBlue);
	target->AddNode(tar_mat, 1);

	return target;
}

void BeWork::FillSimFile(const char* outputfile, const Int_t noevents, const char* generator, Option_t *opt,
		Double_t beThetaMin, Double_t beThetaMax,
		UInt_t gseed/*, const char* cutfile*/)
{
	//TTree "simbe" will be filled and saved into outputfile, tree "simbe" contains two branches, "BeSimPhys" with virtually detected particles coming
	//	out of investigated reaction and the branch "BeSimKin" with whole kinematic information
	//outputfile:
	//noevents: number of simulated events, if 0 then all events from generator
	//generator: name of file with generated impulses
	//gseed: seed for beam position random generator
	//opt: root file option
	//beamX: x coordinate position of the center of beam in cm
	//beamY: y coordinate position of the center of beam in cm
	//beamXfwhm: x width of beam in cm
	//beamYfwhm: y width of beam in cm
	//beThetaMin: minimum thetaBe in CM in
	//beThetaMax: maximum thetaBe in CM in
	//beamenergy: beam energy in MeV
	//beamtheta: beam inclination in rad
	//beamphi: beam inclination in rad


	TString ofile = fWorkDir + '/' + outputfile;
//	cout << ofile.Data() << endl;

	TGeoManager *geom = BuildGeometry();

	BeEvent *simevent = new BeEvent(fConfigFile);	//fixme initialize fCutsFile

	//output file containing tree of two branches
	TFile *fw = new TFile(ofile.Data(), opt);
	if (!fw) {
		Error("BeWork::FillSimFile", "File %s was not created", ofile.Data());
		return;
	}
	TTree *tw = new TTree("simbe", "tree of simulated events");
	tw->Bronch("BeSimPhys", "BeEvent", &simevent, 1024000, 99);		//velky buffer	!!!!!!!!!!!!!!!velky problem s vetvi, prilis komprimovana informace

	//beam position information in Friend Tree
	TTree *tbeam = new TTree("sbeam", "tree with beam position");

	//uzavreno v cyklu
	Int_t noEntries = 0;
	TString gen = generator;
	if (gen.Length() != 0) {
		noEntries = CountLines(generator);
		if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
		Info("BeWork::FillSimFile", "File \"%s\" will be used as generator", generator);
		Info("BeWork::FillSimFile", "In file %s %d events will be processed", fw->GetName(), noEntries);
	}//if
	else {
		Info("BeWork::FillSimFile", "In file %s %d events will be processed", fw->GetName(), noevents);
	}
	Info("BeWork::FillSimFile", "Output file %s will be filled", fw->GetName());

	MonteCarloState(noevents, generator, tw, simevent, geom,
			tbeam, beThetaMin, beThetaMax, gseed);

//	dodelat dalsi dva generatory			------ jake?

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

	delete geom;
	delete simevent;

	return;
}

void BeWork::MonteCarloState(const Int_t noevents, const char* generator, TTree *writetree, BeEvent *besimevent, TGeoManager *geometry,
		TTree *beampos,	Double_t reactionAngleMin, Double_t reactionAngleMax, UInt_t gseed)
{
	//if generator =NULL, default internal generator will be used
	//promenne: simevent, reaction
	//			energy losses
	//beamenergy: in MeV
	//beamtheta: in rad
	//beamphi: in rad
	//reactionAngleMin: in rad
	//reactionAngleMax: in rad

	Info("BeWork::MonteCarloState", "Function was called");

	TString g(generator);
	std::ifstream gen;
	//generator file opening
	if (!g.IsNull()) {
		gen.open(g.Data());
		if (!gen.is_open()) {
			Error("BeWork::MonteCarloState", "Physical generator was not opened!");
			return;
		}//if
	}//if

	//variables for use with generator
	Double_t E, p_alpha[3], p_p1[3], p_p2[3], thetaCM;	//E is mass or exc. energy of Be, thetaCM is angle of (p,n) reaction

	//uzavreno v cyklu
	Int_t noEntries = 0;
	if (gen.is_open()) {
		noEntries = CountLines(generator);
		if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
		Info("BeWork::MonteCarloState", "%d events will be processed", noEntries);
	}//if
	else {
		noEntries = noevents;
		Info("BeWork::MonteCarloState", "%d events will be processed", noEntries);
	}


//	Int_t noEntries = CountLines(generator);
//	if (noevents != 0 && noEntries > noevents) { noEntries = noevents; }
//	Info("BeWork::FillSimFile", "In file %s %d events will be processed", fw->GetName(), noEntries);
//	Info("BeWork::FillSimFile", "Output file %s will be filled", fw->GetName());

	//beam size and position: change to function parameter
	TRandom3 ranPosition(gseed);
	TRandom3 ranEnergy(gseed+1);
	TRandom3 ranThetaCM(gseed+3);
	Double_t x = fBeamX_MC;
	Double_t y = fBeamY_MC;
	Double_t z = 0;
	Double_t beamT = 0.;
	beampos->Branch("bx", &x, "bx/D");
	beampos->Branch("by", &y, "by/D");
	beampos->Branch("bz", &z, "bz/D");
	beampos->Branch("thetaMC", &thetaCM, "thetaMC/D");
	beampos->Branch("E_IM", &E, "E_IM/D");

	Double_t stateRatio = 0.;
	beampos->Branch("sRatio", &stateRatio, "sRatio/D");

	TLorentzVector *beamLab = new TLorentzVector();
	TLorentzVector *alphaLab = new TLorentzVector();
	TLorentzVector *p1Lab = new TLorentzVector();
	TLorentzVector *p2Lab = new TLorentzVector();
	beampos->Bronch("beamLab.", "TLorentzVector", &beamLab, 1024000, 99);
	beampos->Bronch("sALab.", "TLorentzVector", &alphaLab, 1024000, 99);
	beampos->Bronch("sP1Lab.", "TLorentzVector", &p1Lab, 1024000, 99);
	beampos->Bronch("sP2Lab.", "TLorentzVector", &p2Lab, 1024000, 99);

	TLorentzVector *alphaCM = new TLorentzVector();
	TLorentzVector *p1CM = new TLorentzVector();
	TLorentzVector *p2CM = new TLorentzVector();
	TLorentzVector *beCM = new TLorentzVector();
	beampos->Bronch("sACM.", "TLorentzVector", &alphaCM, 1024000, 99);
	beampos->Bronch("sP1CM.", "TLorentzVector", &p1CM, 1024000, 99);
	beampos->Bronch("sP2CM.", "TLorentzVector", &p2CM, 1024000, 99);
	beampos->Bronch("sBeCM.", "TLorentzVector", &beCM, 1024000, 99);

	//Jacobi "T" system
	const Double_t MTx = BeEvent::ReducedMass(kPMASS, kPMASS);
	const Double_t MTy = BeEvent::ReducedMass(kPMASS + kPMASS, kAMASS);
	TVector3 kTx;
	TVector3 kTy;

	Double_t sTpp;
	Double_t sTapp;
	Double_t sCosThetaTk;
	beampos->Branch("sTpp", &sTpp, "sTpp/D");
	beampos->Branch("sTapp", &sTapp, "sTapp/D");
	beampos->Branch("sCosThetaTk", &sCosThetaTk, "sCosThetaTk/D");

	//Jacobi "Y" system
	const Double_t MYx = BeEvent::ReducedMass(kPMASS, kAMASS);
	const Double_t MYy = BeEvent::ReducedMass(kPMASS + kAMASS, kPMASS);
	TVector3 kYx;
	TVector3 kYy;

	Double_t sTap;
	Double_t sTpap;
	Double_t sCosThetaYk;
	beampos->Branch("sTap", &sTap, "sTap/D");
	beampos->Branch("sTpap", &sTpap, "sTpap/D");
	beampos->Branch("sCosThetaYk", &sCosThetaYk, "sCosThetaYk/D");

	BeReaction reaction;

	//detectors resolution
	TRandom3 *detres = new TRandom3();
	TRandom3 *proton_mix = new TRandom3();

//	const Double_t beamResolution =

	for (Int_t i = 0; i < noEntries; i++) {
//		cout << i << endl;
		printProgBar(i, noEntries);
		//particle tracking
		besimevent->Reset();
		reaction.Reset();

		//reaction position in the thickness of the target
		z = ranPosition.Uniform(-T_THICKNESS/10000./2., T_THICKNESS/10000./2.);		//in cm
		beamT = ranEnergy.Gaus(fTBeamMC*6, fTBeamResMC);
		beamT = fTargetLi->GetE(beamT, z*CmToMic());
		beamLab->SetPxPyPzE(0., 0., TMath::Sqrt(BeEvent::PC2(beamT, kLiMASS)),
				beamT + kLiMASS);
//		printf("%3.2f\t%3.2f\n", beamLab->E(), beamT + kLiMASS);
//		printf("%3.2f\t%3.2f\n", beamLab->Pz(), TMath::Sqrt(BeEvent::PC2(beamT, kLiMASS)));
//		cout << beamLab->Px() << endl
//				<< beamLab->Py() << endl;
//		cout << beamT << endl;
//		Info("BeWork::MonteCarloState", "beam z is %f mic, beam energy is %f MeV", z*CmToMic(), benergy);
//		Info("BeWork::MonteCarloState",
//				"beam z is %f mic, beam energy is %f MeV\n", -z * CmToMic(),
//				fTargetLi->GetE0(beamenergy, -z * CmToMic()) );

		//phase volume
		if (g.IsNull()) {
//			reaction.FillProcess(beamenergy, 0., 0.);			//pomala fce
			reaction.FillProcess(beamT, 0., 0.);			//pomala fce
			thetaCM = reaction.fLip.GetThetaCM();
			E = reaction.fLip.GetMa();
		}
		//generator
		else {
			if (gen.good()) {
				gen >> E >> p_alpha[0] >> p_alpha[1] >> p_alpha[2] >>
							p_p1[0] >> p_p1[1] >> p_p1[2] >>
							p_p2[0] >> p_p2[1] >> p_p2[2] >> stateRatio;
				if (E > 20) continue;		//to high energies cannot be populated with our low energy beam
				//fixme beam energy is different for different z
//				reaction.FillProcess(beamenergy, 0., 0., p_alpha, p_p1, p_p2, reactionAngleMin, reactionAngleMax);

				if (stateRatio < fsRatioMin || stateRatio > fsRatioMax) continue;	//to omit unusable events from generator

				alphaCM->SetPxPyPzE(p_alpha[0], p_alpha[1], p_alpha[2],
						BeEvent::T(
								Power(p_alpha[0], 2) + Power(p_alpha[1], 2)
										+ Power(p_alpha[2], 2), kAMASS) + kAMASS);
				p1CM->SetPxPyPzE(p_p1[0], p_p1[1], p_p1[2],
						BeEvent::T(
								Power(p_p1[0], 2) + Power(p_p1[1], 2)
								+ Power(p_p1[2], 2), kPMASS) + kPMASS);
				p2CM->SetPxPyPzE(p_p2[0], p_p2[1], p_p2[2],
						BeEvent::T(
								Power(p_p2[0], 2) + Power(p_p2[1], 2)
								+ Power(p_p2[2], 2), kPMASS) + kPMASS);
				*beCM = *p1CM + *p2CM + *alphaCM;

				//Jacobi "T" system
				kTx = (p1CM->Vect() - p2CM->Vect()) * 0.5;
				kTy =
						(4. * (p1CM->Vect() + p2CM->Vect())
								- 2. * alphaCM->Vect()) * (1. / 6.);

				sTpp = kTx.Mag2() / (2 * MTx);
				sTapp = kTy.Mag2() / (2 * MTy);
				sCosThetaTk = kTx.Dot(kTy) / (kTx.Mag() * kTy.Mag());

				//Jacobi "Y" system
				kYx = (4. * p1CM->Vect() - alphaCM->Vect()) * (1./5.);
				kYy = ((p1CM->Vect() + alphaCM->Vect())
								- 5. * p2CM->Vect()) * (1. / 6.);

				sTap = kYx.Mag2() / (2 * MYx);
				sTpap = kYy.Mag2() / (2 * MYy);
				sCosThetaYk = kYx.Dot(kYy) / (kYx.Mag() * kYy.Mag());

				thetaCM = ranThetaCM.Uniform(reactionAngleMin, reactionAngleMax);
//				reaction.FillProcess(beamT, 0., 0., p_alpha, p_p1, p_p2,
//						reactionAngleMin, reactionAngleMax);
				reaction.FillProcess(beamT, 0., 0., p_alpha, p_p1, p_p2,
						thetaCM);
				if (!gen.good()) {
					Warning("BeWork::FillSimKinFile", "End of file %s was reached at %d event ", generator, i);
					break;
				}//if
			}//if
		}//else

		TLorentzVector vAlpha(reaction.GetAlpha());
		TLorentzVector vP1;
		TLorentzVector vP2;


		if (proton_mix->Uniform() < 0.5) {
			vP1 = reaction.GetP1();
			vP2 = reaction.GetP2();
		}
		else {
			vP1 = reaction.GetP2();
			vP2 = reaction.GetP1();
		}

//		alphaT = vAlpha.E() - kAMASS;
//		p1T = vP1.E() - kPMASS;
//		p2T = vP2.E() - kPMASS;

		//beam position
		//beam size
		//Gauss
		x = ranPosition.Gaus(fBeamX_MC, fBeamX_sigma_MC);
		y = ranPosition.Gaus(fBeamY_MC, fBeamY_sigma_MC);
		//rectangle (square)
//		x = ranPosition.Uniform(-beamXfwhm, beamXfwhm);		//Gaus(beamX, beamXfwhm/2.35);
//		y = ranPosition.Uniform(-beamYfwhm, beamYfwhm);

		//reaction position in the thickness of the target
//		z = ranPosition.Uniform(-T_THICKNESS/10000./2., T_THICKNESS/10000./2.);		//in cm
//		Info("BeWork::MonteCarloState", "beam z is %f", z*CmToMic());
		//elipse (circle) is not possible
		//should be treated on level of TCut

//		Info("BeWork::MonteCarloState", "Input energy of alpha is %f", alphaT);

		*alphaLab = vAlpha;
		*p1Lab = vP1;
		*p2Lab = vP2;

		ParticleTracking(&vAlpha, 24,
				geometry, besimevent, detres, x, y, z);	//funguje pomerne rychle
		ParticleTracking(&vP1, 11,
				geometry, besimevent, detres, x, y, z);
		ParticleTracking(&vP2, 11,
				geometry, besimevent, detres, x, y, z);

		//work with energies to fill the all BeEvent variables
		besimevent->SetDeltaE();
		besimevent->SetHits();

//		fw->cd();

//		writetree->Fill();			//	pomala funkce
//		beampos->Fill();


		if (fBeOnly) {
			if (besimevent->WriteCondition()) {
				writetree->Fill();
				beampos->Fill();
			} //zkontrolovat
		} else {
			writetree->Fill();			//	pomala funkce
			beampos->Fill();
		}

	}//for all events

	//closing of generator file
	if (!g.IsNull()) {
		gen.close();
		if (gen.is_open()) {
			Warning("BeWork::FillSimKinFile", "File %s closing error\n", generator);
		}//if
	}//if

	//delete all TLorentzVectors:
	//fixme

	delete detres;
	delete proton_mix;

	delete beamLab;
	delete alphaLab;
	delete p1Lab;
	delete p2Lab;

	delete alphaCM;
	delete p1CM;
	delete p2CM;

	return;
}

void BeWork::ParticleTracking(const TLorentzVector *particle, const Int_t pID,
		TGeoManager *geom, BeEvent *event, TRandom3 *detres,
		Double_t beamX, Double_t beamY, Double_t beamZ)
{
	//particle: treated particle
	//silosses: object for energy losses calculations in Si detectors
	//csilosses: object for energy losses calculations in CsI(Tl) detectors
	//geom: geometry to be used
	//event: object of event class to be filled
	//detres:
	//beamX: in cm
	//beamY: in cm

	Double_t kinE0 = 0;

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

	if (pID == kPROTON) {	//proton
		kinE0 = particle->E() - kPMASS;
		silosses = fSiP;
		csilosses = fCsIP;
		targetlosses = fTargetP;
		target_win_losses = fTargetWinP;
	}
	else {
		if (pID == kALPHA) {	//alpha
			kinE0 = particle->E() - kAMASS;
			silosses = fSiAlpha;
			csilosses = fCsIAlpha;
			targetlosses = fTargetAlpha;
			target_win_losses = fTargetWinAlpha;
		}
		else {
			Info("BeWork::ParticleTracking",
					"I can only track either protons or alphas.");
			return;
		}
	}//else
	Double_t kinE = 0;
	const Double_t cos_theta = TMath::Cos(particle->Theta());

	TVector3 pardir(particle->Vect().Unit());	//particle direction
	geom->InitTrack(beamX, beamY, beamZ, pardir.Px(), pardir.Py(), pardir.Pz());

	TString anode;
	Double_t dist;	//in cm
	TString detector;	//mother node
	Int_t	ring, sec;
	const Int_t kNoStrip = -1;
	Double_t deltaE = 0;
	const Double_t fdlCsI = fX13_FD;	//in mcm
//	const Double_t fdlCsI = 1.;	//in mcm

	while (!geom->IsOutside()) {
		if (kinE0 == 0) {
			break; }
		geom->FindNextBoundary();
		ring = kNoStrip;
		sec = kNoStrip;
		while (geom->GetMother()) {	//inside Si detector	//XXX problem: overit
															//funguje, pouze pokud se castice nezastavi, jinak se zasekne
			geom->FindNextBoundary();
			anode = geom->FindNode()->GetName();
			if (anode.Contains("senssec") ) {	//as dead layer
				//energy losses in front dead layer of Si detector
				sec = geom->FindNode()->GetNumber();
				dist = geom->GetStep();
//				cout << dist << endl;
				kinE0 = silosses->GetE(kinE0, CmToMic()*dist);
			}//if
			if (anode.Contains("discSi") ) {
				//energy losses in sensitive layer of Si detector
				dist = geom->GetStep();
				kinE = silosses->GetE(kinE0, CmToMic()*dist);
				deltaE = kinE0 - kinE;
				kinE0 = kinE;
				detector = geom->GetMother()->GetName();	//used to identify detector
															//for write the energy
			}//if
			if (anode.Contains("ring") ) {	//as dead layer
				//energy losses in back dead layer of Si detector
				ring = geom->FindNode()->GetNumber();
				dist = geom->GetStep();
				kinE0 = silosses->GetE(kinE0, CmToMic()*dist);
			}//if
			if (anode.Contains("csSec") ) { 	//energy losses in CsI detector
//				cout << "CsI" << endl;
				sec = geom->FindNode()->GetNumber();
				ring = 0;
				//dead layer
				kinE0 = csilosses->GetE(kinE0, fdlCsI/cos_theta);
				//sensitive layer
				dist = geom->GetStep();
				kinE = csilosses->GetE(kinE0, CmToMic()*dist);
				deltaE = kinE0 - kinE;
				kinE0 = kinE;
				detector = geom->GetMother()->GetName();		//used to identify detector
																//for write the energy

//				printf("CsI: %d\t%3.2f\t%3.2f\t%s\n", sec, dist, deltaE, detector.Data());


			}//if
			if (anode.Contains("tarmaterial") ) { 	//energy losses in target
				dist = geom->GetStep();
				kinE0 = targetlosses->GetE(kinE0, CmToMic()*dist);
			}//if
			if (anode.Contains("target") ) { 	//energy losses in target windows
				dist = geom->GetStep();
				kinE0 = target_win_losses->GetE(kinE0, CmToMic()*dist);
			}//if

			geom->Step();
		}//while inside detector

		Double_t E_write;

		//write in the simulated energy deposites
		if ( event && (sec != kNoStrip) && (ring != kNoStrip) ) {
			//first telescope
			if (detector.Contains("T11")) {
				E_write = deltaE + detres->Gaus(0., fSiRes);			//calculated energy loss + detector resolution
				if (E_write >= event->calSi->fC[3][0]) event->fT11.fESec[sec/2] = E_write;
				E_write = deltaE + detres->Gaus(0., fSiRes);			//calculated energy loss + detector resolution
				if (E_write >= event->calSi->fC[5][0]) event->fT11.fERing[ring] = E_write;
			}//if T11
			if (detector.Contains("T12")) {
				E_write = deltaE + detres->Gaus(0., fSiRes);			//calculated energy loss + detector resolution
				if (E_write >= event->calSi->fC[13][sec/4]) event->fT12.fESec[sec/4] = E_write;
				event->fTau[0][sec/4] = 90.;
			}//if T12
			if (detector.Contains("T13")) {
				if (deltaE*fCsIResA>fCsIBestResA) E_write = deltaE + deltaE*detres->Gaus(0., fCsIResA/2.35); 		//pro alfy				//calculated energy loss + detector resolution
				else E_write = deltaE + deltaE*detres->Gaus(0., fCsIBestResA/2.35);
				if (E_write >= event->calCsIa->fD[15][sec])	event->fE13aSec[sec] = E_write;
				if (deltaE*fCsIResP>fCsIBestResP) E_write = deltaE + deltaE*detres->Gaus(0., fCsIResP/2.35);								//calculated energy loss + detector resolution
				else E_write = deltaE + deltaE*detres->Gaus(0., fCsIBestResP/2.35);
				if (E_write >= event->calCsIp->fD[15][sec]) event->fE13pSec[sec] = E_write;
			}//if T13

			//second telescope
			if (detector.Contains("T21")) {
				E_write = deltaE + detres->Gaus(0., fSiRes);			//calculated energy loss + detector resolution
				if (E_write >= event->calSi->fC[8][0]) event->fT21.fESec[sec/2] = E_write;
				E_write = deltaE + detres->Gaus(0., fSiRes);			//calculated energy loss + detector resolution
				if (ring < 16 && E_write >= event->calSi->fC[10][0]) event->fT21.fERing[ring] = E_write;
				if (ring > 15 && E_write >= event->calSi->fC[11][0]) event->fT21.fERing[ring] = E_write;
			}//if T21
			if (detector.Contains("T22")) {
				E_write = deltaE + detres->Gaus(0., fSiRes);			//calculated energy loss + detector resolution
				if (E_write >= event->calSi->fC[14][sec/4]) event->fT22.fESec[sec/4] = E_write;
				event->fTau[1][sec/4] = 90.;
			}//if T22
			if (detector.Contains("T23")) {
				if (deltaE*fCsIResA>fCsIBestResA) E_write = deltaE + deltaE*detres->Gaus(0., fCsIResA/2.35); 		//pro alfy				//calculated energy loss + detector resolution
				else E_write = deltaE + deltaE*detres->Gaus(0., fCsIBestResA/2.35);
				if (E_write >= event->calCsIa->fD[16][sec]) event->fE23aSec[sec] = E_write;
				if (deltaE*fCsIResP>fCsIBestResP) E_write = deltaE + deltaE*detres->Gaus(0., fCsIResP/2.35);								//calculated energy loss + detector resolution
				else E_write = deltaE + deltaE*detres->Gaus(0., fCsIBestResP/2.35);
				if (E_write >= event->calCsIp->fD[16][sec]) event->fE23pSec[sec] = E_write;
			}//if T23
		}//if writein

		geom->Step();
	}//while tracking through whole detector system

	return;
}

void BeWork::FillSimKinFile(Int_t INPUTS, Double_t beamtheta, Double_t beamphi,
		const char* outputfile, Option_t *opt) {
	//simulation of kinematic variables using BinaryReaction class;
	//make TTree object named "KIN" in file
	//INPUTS: number of events
	//outputfilename: name of root file
	//opt: file access option

	BeReaction *sReaction = new BeReaction();

	Info("BeWork:FillSimKinFile()",
			"MC simulation using phase volume generator will be done");

	if (INPUTS == 0) {
		cout << "\nAssign number of events \n";
		return;
//		cin >> INPUTS;
	}

	cout << "\n\nCalculating...\n\n";

	TFile *fw = TFile::Open(outputfile, opt);
	if (!fw) {
		Error("BeWork::FillSimKinFile", "File %s was not created", outputfile);
		return;
	}
	TTree *tree = new TTree("KIN", "two body kinematics");
	tree->Bronch("MCsimple", "BeReaction", &sReaction);

	for (Int_t i = 1; i <= INPUTS; i++) {
		printProgBar(i, INPUTS);

		sReaction->FillProcess(fTBeamMC*6, beamtheta, beamphi);
		tree->Fill();
	}

	cout << endl;

	tree->Write();

	fw->Close();
	if (fw->IsOpen()) {
		Warning("BeWork::FillSimKinFile", "File %s closing error", outputfile);
	}

	delete fw;
	delete sReaction;

}

void BeWork::FillSimKinFile(const char* generator, Int_t INPUTS,
		/*Double_t beamenergy,*/ Double_t beamtheta, Double_t beamphi,
		const char* outputfile, Option_t *opt) {
//	function variables: path to generator, INPUTS, outputfile option

	Info("BeWork::FillSimKinFile",
			"MC simulation using physical generator will be done");

//	cout << beamenergy << endl;
//	cout << beamtheta << endl;
//	cout << beamphi << endl;
	std::ifstream gen;

	gen.open(generator);
	if (!gen.is_open()) {
		Error("BeWork::FillSimKinFile", "Physical generator was not opened!");
		return;
	} //if

	//class which will be the structure of TTree branch filled into
	BeReaction *sReaction = new BeReaction();

	//checking number of inputs
	if (INPUTS == 0) {
		Info("FillSimKinFile", "O inputs was filled");
		return;
//		cout << "\nAssign number of events \n";
//		cin >> INPUTS;
	}

	cout << "\n\nCalculating...\n\n";

	//outputfile, tree and branch creation
	TFile *fw = TFile::Open(outputfile, opt); //name of the output file
	if (!fw) {
		Error("BeWork::FillSimKinFile", "File %s was not opened", outputfile);
		return;
	}
	TTree *tree = new TTree("KIN", "some suitable title"); //name of the tree, title of the three
	tree->Bronch("MCgen", "BeReaction", &sReaction); //name of the branch

	//filling the tree

	Double_t E, p_alpha[3], p_p1[3], p_p2[3], thetaCM;

	for (Int_t i = 1; i <= INPUTS; i++) {
		printProgBar(i, INPUTS); //zrychlit fci
		if (gen.good()) {
			gen >> E >> p_alpha[0] >> p_alpha[1] >> p_alpha[2] >> p_p1[0]
					>> p_p1[1] >> p_p1[2] >> p_p2[0] >> p_p2[1] >> p_p2[2]
					>> thetaCM;
//			sReaction->FillProcess(beamenergy, beamtheta, beamphi, p_alpha, p_p1, p_p2, thetaCM);
			sReaction->FillProcess(fTBeamMC*6, beamtheta, beamphi, p_alpha,
					p_p1, p_p2, 0., TMath::Pi());
			if (!gen.good()) {
				Warning("BeWork::FillSimKinFile", "End of file %s was reached",
						generator);
				tree->Fill();
				break;
			} //if
		} //if
		tree->Fill();
	} //for

	cout << endl;

	tree->Write();

	fw->Close();
	if (fw->IsOpen()) {
		Warning("BeWork::FillSimKinFile", "File %s closing error", outputfile);
	}

	gen.close();
	if (gen.is_open()) {
		Warning("BeWork::FillSimKinFile", "File %s closing error\n", generator);
	} //if

	delete fw;
	delete sReaction;

	return;
}

void BeWork::printProgBar(Int_t percent)
{
//	std::string bar;
	TString bar;

	for(int i = 0; i < 50; i++){
		if( i < (percent/2)){
//			bar.replace(i,1,"=");
			bar.Replace(i,1,"=");
		}
		else {
			if( i == (percent/2)){
//				bar.replace(i,1,">");
				bar.Replace(i,1,">");
			}
			else {
//				bar.replace(i,1," ");
				bar.Replace(i,1," ");
			}
		}
	}

	cout << "\r" "[" << bar << "] ";
	cout.width( 3 );
//	cout << percent << "%     " << flush;
	if (percent < 100) cout << percent << "%     " << flush;
	else cout << percent << "%     " << endl;

	return;
}

void BeWork::printProgBar(Int_t i, Int_t inputs)
{
//	printProgBar(0);

//	if ( (i%(inputs/100)) != 0 ) return;

	Int_t percent = i*100/inputs;
	if ( inputs >= 100 && (i%(inputs/100)) == 0 ) {
		printProgBar(percent);
	}
	if (i == inputs-1) { printProgBar(100); }
	return;
}

Int_t BeWork::CountLines(const char* file, Int_t maxlinelength) {

	char countline[maxlinelength];

	TString fname(file);
	std::ifstream rfile;
	//generator file opening
	if (!fname.IsNull()) {
		rfile.open(fname.Data());
		if (!rfile.is_open()) {
			Error("BeWork::CountLines", "File \"%s\" was not opened!", fname.Data());
			return 0;
		}//if
	}//if

	Int_t j = 0;
	while (!rfile.eof()) {
		rfile.getline(countline, maxlinelength);
		j++;
	}

	rfile.close();
	if (rfile.is_open()) {
		Warning("BeWork::CountLines", "File %s closing error\n", fname.Data());
	}//if

	return j;
}

void BeWork::CreateTELosses() {
	//energy looses in silicon
	//alpha in Si
	fSiAlpha = new TELoss();
	fSiAlpha->SetEL(SI_NELEMENTS, SI_RHO);
	fSiAlpha->AddEL(SI_1_Z, SI_1_A, SI_1_W);
	fSiAlpha->SetZP(2., 4.);		//alphas
	//100000 je urcite v poradku, 80000 staci, 50000 malo
	fSiAlpha->SetEtab(100000, 200.);	//deltaE
	//proton in Si
	fSiP = new TELoss();
	fSiP->SetEL(SI_NELEMENTS, SI_RHO);
	fSiP->AddEL(SI_1_Z, SI_1_A, SI_1_W);
	fSiP->SetZP(1., 1.);		//protons
	fSiP->SetEtab(100000, 200.);	//deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

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

	//energy looses in target
	const Double_t targetRho = T_NORMAL_RHO*(273/T_TEMPERATURE)*(T_PRESSURE/1);
	//alpha in target
	fTargetAlpha = new TELoss();
	fTargetAlpha->SetEL(T_NELEMENTS, targetRho);
	fTargetAlpha->AddEL(T_1_Z, T_1_A, T_1_W);
	fTargetAlpha->SetZP(2., 4.);		//alphas
	//100000 je urcite v poradku, 80000 staci, 50000 malo
	fTargetAlpha->SetEtab(100000, 200.);	//deltaE
	//proton in target
	fTargetP = new TELoss();
	fTargetP->SetEL(T_NELEMENTS, targetRho);
	fTargetP->AddEL(T_1_Z, T_1_A, T_1_W);
	fTargetP->SetZP(1., 1.);		//protons
	fTargetP->SetEtab(100000, 200.);	//deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

	//beam in target
	fTargetLi = new TELoss();
	fTargetLi->SetEL(T_NELEMENTS, targetRho);
	fTargetLi->AddEL(T_1_Z, T_1_A, T_1_W);
	fTargetLi->SetZP(3., 6.);		//protons
	fTargetLi->SetEtab(100000, 200.);	//deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo


	//alpha in target window
	fTargetWinAlpha = new TELoss();
	fTargetWinAlpha->SetEL(T_WIN_NELEMENTS, T_WIN_RHO);
	fTargetWinAlpha->AddEL(T_WIN_1_Z, T_WIN_1_A, T_WIN_1_W);
	fTargetWinAlpha->AddEL(T_WIN_2_Z, T_WIN_2_A, T_WIN_2_W);
	fTargetWinAlpha->AddEL(T_WIN_3_Z, T_WIN_3_A, T_WIN_3_W);
	fTargetWinAlpha->SetZP(2., 4.);		//alphas
	//100000 je urcite v poradku, 80000 staci, 50000 malo
	fTargetWinAlpha->SetEtab(100000, 200.);	//deltaE
	//proton in target window
	fTargetWinP = new TELoss();
	fTargetWinP->SetEL(T_NELEMENTS, T_WIN_RHO);
	fTargetWinP->AddEL(T_WIN_1_Z, T_WIN_1_A, T_WIN_1_W);
	fTargetWinP->AddEL(T_WIN_2_Z, T_WIN_2_A, T_WIN_2_W);
	fTargetWinP->AddEL(T_WIN_3_Z, T_WIN_3_A, T_WIN_3_W);
	fTargetWinP->SetZP(1., 1.);		//protons
	fTargetWinP->SetEtab(100000, 200.);	//deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

	//beam in target window
	fTargetWinLi = new TELoss();
	fTargetWinLi->SetEL(T_NELEMENTS, T_WIN_RHO);
	fTargetWinLi->AddEL(T_WIN_1_Z, T_WIN_1_A, T_WIN_1_W);
	fTargetWinLi->AddEL(T_WIN_2_Z, T_WIN_2_A, T_WIN_2_W);
	fTargetWinLi->AddEL(T_WIN_3_Z, T_WIN_3_A, T_WIN_3_W);
	fTargetWinLi->SetZP(3., 6.);		//protons
	fTargetWinLi->SetEtab(100000, 200.);	//deltaE		//100000 je urcite v poradku, 80000 staci, 50000 malo

	Info("BeWork::CreateTELosses", "TELoss objects were initialized.");
}