//////////////////////////////////////////////////////////
//									//
//	AculCalibration					//
//									//
//
//Some description of this very useful class,
//its properties and a short example how to
//use it. This text may be partly used in
//PhD thesis.
//
//Description of the detector itself.
//
//														//
//////////////////////////////////////////////////////////

#include "AculCalibration.h"

ClassImp(AculCalibration);

AculCalibration::AculCalibration()
{
	//default constructor

	fCurrentHStack = NULL;
	fCurrentHistList.IsOwner();

	kRaNOPEAKS = 0;
	fLowerPeakRelativeHight = 0.;
	fUpperPeakRelativeHight = 0.;
	fPeakPositionTolerance = 0.;
	fFitFuncLineWidth = 1;
	fFitMinSigma = 0.;
	fFitPeakThreshold = 0.;

	for(Int_t i = 0; i < DEFAULTNOPEAKS; i++) {
		fEnergy[i] = 0;
	}
	fCalInformation = 0;

	Reset();

}

AculCalibration::AculCalibration(const char* parfile)
{
	//constructor which fills fA, fB, fC, fD from file parfile

	fCurrentHStack = NULL;
	fCurrentHistList.IsOwner();

	kRaNOPEAKS = 0;
	fLowerPeakRelativeHight = 0.;
	fUpperPeakRelativeHight = 0.;
	fPeakPositionTolerance = 0.;
	fFitFuncLineWidth = 1;
	fFitMinSigma = 0.;
	fFitPeakThreshold = 0.;

	for(Int_t i = 0; i < DEFAULTNOPEAKS; i++) {
		fEnergy[i] = 0;
	}

	fCalInformation = 0;

	SetCalibrationParameters(parfile);

}

AculCalibration::~AculCalibration()
{

	DeleteStacks();
//	delete fCalInformation;
//	fCalInformation->Close();

}

Int_t AculCalibration::SearchPeaks(const TH1 *hin, Double_t sigma, Option_t *option, const Int_t searchedpeaks)
{
	//Function searching peaks in inputed TH1 spectrum and selects the peaks in the histogram.
	//
	//  hin:
	//  sigma:
	//  option:
	//  threshold:
	//  searchedpeaks:

	TSpectrum sc;	//by default for 100 peaks
	Int_t nopeaks = sc.Search(hin, sigma, "goff", fFitPeakThreshold);

	TString opt = option;
	opt.ToLower();

	const Double_t tStep = 0.05;

	while ( nopeaks > searchedpeaks && fFitPeakThreshold <= 1) {
		fFitPeakThreshold = fFitPeakThreshold + tStep;
		nopeaks = sc.Search(hin, sigma, "goff", fFitPeakThreshold);
	}

	if (!nopeaks) {
		return 0;
	}

	if (opt.Contains("goff")) {
		return nopeaks;
	}

	TPolyMarker *pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
	if (pm) {
		hin->GetListOfFunctions()->Remove(pm);
		delete pm;
	}
	pm = new TPolyMarker(nopeaks, sc.GetPositionX(), sc.GetPositionY());
	hin->GetListOfFunctions()->Add(pm);
	pm->SetMarkerStyle(23);
	pm->SetMarkerColor(kRed);
	pm->SetMarkerSize(1.3);

	return nopeaks;
}

Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t sigmamin)
{

	if (!hSpectrum) return 1;
	Int_t dimension = hSpectrum->GetDimension();
	if (dimension > 1) {
		Error("PeaksFitting", "Only implemented for 1-d histograms");
		return 1;
	}

	TString	opt = option;
	opt.ToLower();

	if (!kRaNOPEAKS) {
		Error("PeaksFitting", "kRaNOPEAKS is set to zero; calibration spectrum must be set");
		return 1;
	}

	Int_t peaksNumber =	SearchPeaks(hSpectrum, sigmamin, "", kRaNOPEAKS);

	if (peaksNumber != kRaNOPEAKS) {
		Info("PeaksFitting", "In histogram %s was found %d peaks", hSpectrum->GetName(), peaksNumber);
		return 1;
	}

	//should be optional output
	Info("PeaksFitting", "Number of peaks in %s: %d", hSpectrum->GetName(), peaksNumber);

	//working array for peaks, there are founded in accidental order
	Double_t peak[peaksNumber];
	Double_t *peakPosition;
	Double_t *peakHight;

	TList *functions = hSpectrum->GetListOfFunctions();
	TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");

	peakPosition = pm->GetX();
	peakHight = pm->GetY();

	for (Int_t i = 0; i < peaksNumber; i++) {

		Double_t fitMin = 0;
		Double_t fitMax = 0;
		Double_t fitStep = hSpectrum->GetXaxis()->GetBinWidth(0);
//		cout << fitStep << endl;
//		cout << fLowerPeakRelativeHight << "\t" << fUpperPeakRelativeHight << endl;

		//fitting range:
		//shift a range of fit and search for raw boarder of peak determined by fUpperPeakRelativeHight
		//maximum
		Int_t j = 0;
		Double_t currentHight = peakHight[i];
		while ( currentHight > (peakHight[i]*fUpperPeakRelativeHight) ) {
			j++;
			fitMax = static_cast<Double_t>(peakPosition[i]) + j*fitStep;
			currentHight = hSpectrum->GetBinContent(hSpectrum->GetXaxis()->FindBin(fitMax));
		}

		//minimum
		j = 0;
		currentHight = peakHight[i];
		while ( currentHight > (peakHight[i]*fLowerPeakRelativeHight) ) {
			j++;
			fitMin = static_cast<Double_t>(peakPosition[i]) - j*fitStep;
			currentHight = hSpectrum->GetBinContent(hSpectrum->GetXaxis()->FindBin(fitMin));
		}

		//fitting
		if (opt.Contains("gp")) {
			Info("PeaksFitting", "Option containing gp");
			char fncname[20];
			sprintf(fncname, "gaus_aux_%d", i);
			TF1 *gausAux = new TF1(fncname, "gaus", fitMin - 10, fitMax + 10);		//pomocny gaus
			hSpectrum->Fit(fncname, "0 Q", "", fitMin - 15, fitMax + 15);				//prvotni fitovani

			sprintf(fncname, "auto_gp_%d", i);
			TF1 *fitAuto = new TF1(fncname, "gaus(0) + pol0(3)", fitMin - 15, fitMax + 15);		//fce pro automaticke fitovani
			fitAuto->SetParameter(0, gausAux->GetParameter(0));		//nastavovani parametru fitovaci fce
			fitAuto->SetParameter(1, gausAux->GetParameter(1));
			fitAuto->SetParameter(2, gausAux->GetParameter(2));

			hSpectrum->Fit(fncname, "0 R Q +", "", fitMin - 15, fitMax + 15); //dodelat zapis vsech fci
			hSpectrum->GetFunction(fncname)->ResetBit(TF1::kNotDraw);
			peak[i] = fitAuto->GetParameter(1);			//zapis asi pozice v kanalech do pomocneho pole
			if (opt.Contains("V")) {
				Info("PeaksFitting", "Peak position is\t %4.2f \tresolution is \t %2.1f %%", fitAuto->GetParameter(1), 235*(fitAuto->GetParameter(2))/(fitAuto->GetParameter(1)));
			}
		}
		else {
			char fncname[20];
			sprintf(fncname, "auto_g%d", i);
			TF1 *fitAuto = new TF1(fncname, "gaus", fitMin, fitMax);		//fce pro automaticke fitovani
//			cout << fitMin << "\t" << fitMax << endl;
//			fitAuto->SetParameter(2, fitMax-fitMin);
			fitAuto->SetLineWidth(fFitFuncLineWidth);
			hSpectrum->Fit(fncname, "+ 0 R Q", ""/*, fitMin - 1, fitMax + 1*/);
//			hSpectrum->GetFunction(fncname)->ResetBit(TF1::kNotDraw);
			hSpectrum->GetFunction(fncname)->InvertBit(TF1::kNotDraw);
			peak[i] = fitAuto->GetParameter(1);			//zapis asi pozice v kanalech do pomocneho pole
			if (opt.Contains("v")) {
				Info("PeaksFitting", "Peak position is\t%4.2f\tresolution is \t%2.1f %%", fitAuto->GetParameter(1), 235*(fitAuto->GetParameter(2))/(fitAuto->GetParameter(1)));
			}
		}//else
		//end of fitting
	}//for over all analyzed peaks

	//peaks sorting
	Int_t j[peaksNumber];
	TMath::Sort(peaksNumber, peak, j, kFALSE);
	for (Int_t i = 0; i < 4; i++) {
		fPeak[i] = peak[j[i]];
	}

	if (!opt.Contains("q") || opt.Contains("v")) {
		Info("PeaksFitting", "Control output:");
		for (Int_t i = 0; i < peaksNumber; i++) {
			printf("\tPeak position is\t%f\n", fPeak[i]);
		}
	}

	//	provest kontrolu pomerne polohy piku,
	//	jestli jsou spatne, provest urcita opatreni,
	//	napr. zapis daneho histogramu do souboru,
	//	zapis do souboru s chybama, vypis na obrazovku, ...
	for (Int_t i = 0; i < peaksNumber; i++) {
		if ( !( (((1-fPeakPositionTolerance)*(fEnergy[0]/fEnergy[i])) < (fPeak[0]/fPeak[i])) && (((1+fPeakPositionTolerance)*(fEnergy[0]/fEnergy[i])) > (fPeak[0]/fPeak[i])) ) ) {
			if (fCalInformation && opt.Contains("writebad")) {
				fCalInformation->cd();
				hSpectrum->Write();
			}
			return 2;
		}
	}//for

	return 0;
}


Bool_t AculCalibration::EnergyPositions(const char* inputfile, const char* block,
			const Int_t address, const char* treename, Int_t lowerchannel,
			Int_t upperchannel, Int_t lowersubaddress, Int_t uppersubaddress)
{
	TString iFile = inputfile;
	TFile fr(iFile.Data());



	return 1;
}

Bool_t AculCalibration::SetInputParameters(const char* inputparfile)
{
	const Int_t lineLength = 400;
	Char_t	line[lineLength];
	Char_t	parameter[100];
	Char_t	identificator[100];


	ifstream fipr;
	fipr.open(inputparfile);
	if (!fipr.is_open()) {
		Error("SetInputsParameters", "File with input parameters was not opened");
		return kFALSE;
	}

	while (!fipr.eof()) {

		fipr.getline(line, lineLength);
		if (strlen(line) == 0) {
			continue;
		}
		sscanf(line, "%s %s", parameter, identificator);
		if ( strcmp(identificator, "nopeaks") == 0 ) {
			kRaNOPEAKS = static_cast<Int_t>(atoi(parameter));
			for (Int_t i = 0; i < kRaNOPEAKS; i++) {
				fipr.getline(line, lineLength);
				sscanf(line, "%s", parameter);
				fEnergy[i] = static_cast<Double_t>(atof(parameter));
			}
		}//if

		if ( strcmp(identificator, "lowerchannel") == 0 ) {
			sscanf(line, "%s", parameter);
			fLowerChannel = static_cast<Double_t>(atof(parameter));
		}

		if ( strcmp(identificator, "upperchannel") == 0 ) {
			sscanf(line, "%s", parameter);
			fUpperChannel = static_cast<Double_t>(atof(parameter));
		}

		if ( strcmp(identificator, "lowerpeakhight") == 0 ) {
			sscanf(line, "%s", parameter);
			fLowerPeakRelativeHight = static_cast<Double_t>(atof(parameter));
		}

		if ( strcmp(identificator, "upperpeakhight") == 0 ) {
			sscanf(line, "%s", parameter);
			fUpperPeakRelativeHight = static_cast<Double_t>(atof(parameter));
		}

		if ( strcmp(identificator, "peakpositiontolerance") == 0 ) {
			sscanf(line, "%s", parameter);
			fPeakPositionTolerance = static_cast<Double_t>(atof(parameter));
		}

		if ( strcmp(identificator, "fitfunctionlinewidth") == 0 ) {
			sscanf(line, "%s", parameter);
			fFitFuncLineWidth = static_cast<Width_t>(atoi(parameter));
		}

		if ( strcmp(identificator, "minfitsigma") == 0 ) {
			sscanf(line, "%s", parameter);
			fFitMinSigma = static_cast<Double_t>(atof(parameter));
		}

		if ( strcmp(identificator, "fithightthreshold") == 0 ) {
			sscanf(line, "%s", parameter);
			fFitPeakThreshold = static_cast<Double_t>(atof(parameter));
		}

	}


	fipr.close();

	return kTRUE;

}

void AculCalibration::PrintInputParameters()
{
	//print alpha source parameters

	cout << "Number of peaks: " << kRaNOPEAKS << endl
		<< endl;

	for (Int_t i = 0; i < kRaNOPEAKS; i++) {
		cout << "fEnergy[" << i << "] = " << fEnergy[i] << endl;
	}


	return;

}

Bool_t AculCalibration::SetCalibrationParameters(const char* calparfile)
{

	const Int_t lineLength = 200;
	char line[lineLength];
//	Int_t	crate;
	char	crate[100];
	Int_t	i, j;
	char	cA[40], cB[40], cC[40], cD[40];

	//open file with calibration parameters
	ifstream calFileR;
	calFileR.open(calparfile);


	if( !calFileR.is_open() ) {
		Error("SetCalibrationParameters", "File %s with calibration data was not opened", calparfile);
		return kFALSE;
	}

	Reset();

	//read calibration parameters from file
	while (!calFileR.eof()) {
		calFileR.getline(line, lineLength);
		if ( line[0] != '*' && line[0] != '#' && line[0] != '%' && (line[0] != '/' && line[1] != '/') ) {		//possible comment characters
			sscanf(line, "%s %d %d %s %s %s %s", crate, &i, &j, cA, cB, cC, cD);
			fA[i][j] = atof(cA);
			fB[i][j] = atof(cB);
			fC[i][j] = atof(cC);
			fD[i][j] = atof(cD);
		}

	}
	calFileR.close();
	return kTRUE;
}

void AculCalibration::PrintCalibrationParameters(const Int_t blockmin, const Int_t blockmax)
{
	for (Int_t i = blockmin; i <= blockmax; i++) {
		for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
			cout << "C3[" << i << "][" << j << "]" << setw(10) << fA[i][j] << setw(10) << fB[i][j] << setw(10) << fC[i][j] << setw(10) << fD[i][j] << endl;
		}
	}
}

void AculCalibration::ShowRawSpectra(const char* filename, const Int_t block, TCanvas* rawCanvas, Int_t xaxismin, Int_t xaxismax, /*TObjArray* histList,*/ const Int_t subaddress)
{
	//Displays the spectrum from a file, divides the canvas into a sufficient number of pads and displays spectrums of each  
	//block subaddress on the suitable pads.
	//
        //  filename: input .root file containing spectra to be showed
	//  block: block which will be drawn
	//  rawCanvas: canvas on which one you will see the spectrum
	//  xaxismin: minimal channel, which will be displayed
	//  xaxismax: maximal channel, which will be displayed
	//  subaddress:

	Char_t address[40];
	Char_t histName[40];
	Char_t fillCommand[40];
	Char_t fillCondition[40];


	if (!rawCanvas) {
		//rawCanvas = new TCanvas("RawSpectra", "Raw spectra in channels", 1);
		cout << "You have to assign TCanvas for raw spectra drawing" << endl;
		return;
	}

	rawCanvas->Clear();

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

	rawCanvas->SetFillColor(10);

//	cout << "hovno" << endl;
	TFile *fr = new TFile(filename);
	if (fr->IsOpen() == 0) {
		cout << endl << "File " << filename << " was not opened and won't be processed" << endl << endl;
		return;
	}
	TH1I *hRead = 0;
	TTree *tr = (TTree*)fr->Get("RAW");
//	cout << "hovno" << endl;


	if (subaddress > 15) {
		rawCanvas->Divide(4, 4);
		rawCanvas->SetFillColor(10);
//		cout << "hovno" << endl;
		for (Int_t i = 0; i < 16; i++) {
			cout << i << endl;
			rawCanvas->cd(i+1);
			hRead = new TH1I("name", "title", 4096, 0, 4095);
			sprintf(address, "C3[%d][%d]", block, i);
			sprintf(histName, "H3[%d][%d]", block, i);
			hRead->SetName(histName);
			sprintf(fillCommand, "%s >> %s", address, hRead->GetName());
			sprintf(fillCondition, "%s > 0", address);
			tr->Draw(fillCommand, fillCondition, "");
			if (hRead) {
				hRead->SetDirectory(0);
//				cout << hRead->GetEntries() << endl;
//				if (fHRawList) {
//					fHRawList->Add(hRead);
					fHRawList.Add(hRead);
//				}
				hRead->SetAxisRange(xaxismin, xaxismax);
			}
		}//for
	}
	else {
		fr->cd();
		hRead = new TH1I("name", "title", 4096, 0, 4095);
		sprintf(address, "C3[%d][%d]", block, subaddress);
		sprintf(histName, "H3[%d][%d]", block, subaddress);
		hRead->SetName(histName);
		sprintf(fillCommand, "%s >> %s", address, hRead->GetName());
		sprintf(fillCondition, "%s > 0", address);
//		cout << fillCommand << setw(20) << fillCondition << endl;
		tr->Draw(fillCommand, fillCondition, "goff");
		if (hRead) {
			hRead->SetDirectory(0);
//			if (fHRawList) {
//				fHRawList->Add(hRead);
//			}
			fHRawList.Add(hRead);
			hRead->Draw();
			hRead->SetAxisRange(xaxismin, xaxismax);
		}
	}//else

	fr->Close();

	rawCanvas->Update();

	return;

}

void AculCalibration::ShowSpectra(const char* filename, TCanvas* rawCanvas, Option_t *option, Int_t xaxismin, Int_t xaxismax, const Int_t subaddress)
{
	//filename: input .root file with saved filled histograms to be showed
	//rawCanvas: canvas on which one you will see the spectrum
	//option: THStack options
	//xaxismin: Minimum channel, which will be displayed
	//xaxismax: Maximum channel, which will be displayed
	//subaddress:
 
	TString opt = option;
	opt.ToLower();

	if (!rawCanvas) {
		Error("ShowRawSpectra", "You have to assign TCanvas for raw spectra drawing");
		return;
	}
	rawCanvas->Clear();

	TFile fr(filename);
	if (fr.IsOpen() == 0) {
		Error("ShowRawSpectra", "File %s was not opened and won't be processed", filename);
		return;
	}

	TList *histList;
	histList = fr.GetListOfKeys();
	Int_t listEntries = histList->GetEntries();
	TH1 *hDraw = 0;
	DeleteStacks();

	if (subaddress >= listEntries) {
		fCurrentHStack = new THStack();
		for (Int_t i = 0; i < listEntries; i++) {	//zkontrolovat hranice
			Info("ShowRawSpectra", "Histogram with spectrum of subaddress %d is loading", i);
			fr.GetObject(histList->At(i)->GetName(), hDraw);
			if (hDraw) {
				hDraw->SetDirectory(0);
				//hDraw->SetAxisRange(xaxismax, xaxismin);	//nefunguje
				fCurrentHistList.Add(hDraw);
				fCurrentHStack->Add(hDraw);
			}
		}//for
		if ( !fCurrentHStack->GetHists()->IsEmpty() ) {
			Info("ShowRawSpectra", "Histogram stack drawing");
			fCurrentHStack->Draw(opt.Data());
		}
	}//if all subaddresses
	else {
		//zkontrolovat
		fr.GetObject(histList->At(subaddress)->GetName(), hDraw);
		if (hDraw) {
			hDraw->SetAxisRange(xaxismin, xaxismax, "X");
			hDraw->Draw();
			hDraw->SetDirectory(0);
			fCurrentHistList.Add(hDraw);
		}
	}//else

	fr.Close();
	rawCanvas->Update();
	return;
}

void AculCalibration::FillRawSpectraFile(const char* rawdatafile, const char* block, const char* treename, TCanvas* rawCanvas, Option_t *option, Int_t xaxismin, Int_t xaxismax)
{
	//filename: input .root file containing spectra to be showed
	//block:
	//rawCanvas:
	//xaxismin:
	//xaxismax:

	//variables to be became function parameter
	TString opt(option);
	opt.ToLower();

	if (!rawCanvas) {
		Error("ShowRawSpectra", "You have to assign TCanvas for raw spectra drawing");
		return;
	}
	rawCanvas->Clear();

	TFile fr(rawdatafile);
	if (fr.IsOpen() == 0) {
		Error("ShowRawSpectra", "File %s was not opened and won't be processed", rawdatafile);
		return;
	}
	TTree *tr = (TTree*)fr.Get(treename);

	char outputfile[300];
	sprintf(outputfile, "%s[]Raw.root", block);
	TFile fw(outputfile, opt.Data());
	if (fw.IsOpen() == 0) {
		Error("CalculateCalibParameters", "Output file %s was not created.", outputfile);
		return;
	}
	if (fw.IsWritable() == 0) {
		Error("CalculateCalibParameters", "Output file %s is not writable. Set option to \"RECREATE\".", outputfile);
		return;
	}

	char address[40];
	char histName[40];
	char histTitle[40];
	char fillCommand[40];
	char fillCondition[40];

	fw.cd();
	TH1I *hRead = 0;

	for (Int_t i = 0; i < 16; i++) {	//zkontrolovat hranice
		cout << i << endl;	//predelat na info
		hRead = new TH1I("name", "title", 4096, 0, 4095);
		sprintf(address, "%s[%d]", block, i);
		sprintf(histName, "%s[%d]", block, i);
		sprintf(histTitle, "%s : %s", rawdatafile, histName);
		hRead->SetName(histName);
		hRead->SetTitle(histTitle);
		sprintf(fillCommand, "%s >> %s", address, hRead->GetName());
		sprintf(fillCondition, "%s > 0", address);
		tr->Draw(fillCommand, fillCondition, "goff");		//prozkoumat goff
		hRead->Write();
	}//for

	fw.Close();
	fr.cd();
	delete tr;
	fr.Close();

	rawCanvas->Update();

	return;
}

Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const char* block, const Int_t address,
		const char* treename, Int_t lowerchannel, Int_t upperchannel, Int_t nEBins, Int_t lowersubaddress,
		Int_t uppersubaddress)
{
	if (kRaNOPEAKS == 0) {
		Error("CalculateCalibParameters", "Alpha source parameters was not read");
		return 0;
	}

	if ( (address > BLOCKSNUMBER) || (address < 1) ) {
		Error("CalculateCalibParameters", "Possible address values have to be in range 1 - %d", BLOCKSNUMBER - 1);
		return kFALSE;
	}

	//muzu nechat
	if ( (uppersubaddress - lowersubaddress) >= ADDRESSNUMBER ) {
		Error("CalculateCalibParameters", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1);
		return 0;
	}

	//auxiliary variables, particularly for text parameter fields
	TString oFileName;

	//creation of the output text file
	if ( (lowersubaddress == 0) && (uppersubaddress == ADDRESSNUMBER-1) ) { oFileName.Form("%s[].cal", block); }
	else {
		if (lowersubaddress == uppersubaddress) { oFileName.Form("%s[%d].cal", block, lowersubaddress); }
		else { oFileName.Form("%s[%d-%d].cal", block, lowersubaddress, uppersubaddress); }
	}//if

	ofstream outcalfile;
	outcalfile.open(oFileName.Data());
	if (!outcalfile.is_open()) {
		Error("CalculateCalibParameters", "Output file %s was not opened", oFileName.Data());
		return 0;
	}//if

	//creation of the output root file
	if ( (lowersubaddress == 0) && (uppersubaddress == ADDRESSNUMBER-1) ) { oFileName.Form("%s[].root", block); }
	else {
		if (lowersubaddress == uppersubaddress) { oFileName.Form("%s[%d].root", block, lowersubaddress); }
		else { oFileName.Form("%s[%d-%d].root", block, lowersubaddress, uppersubaddress); }
	}
	fCalInformation = new TFile(oFileName.Data(), "RECREATE");
	if ( !fCalInformation->IsOpen() ) {
		Error("CalculateCalibParameters", "File %s was not opened and won't be processed", oFileName.Data());
		return 0;
	}

	//input file with raw data opening
	TString iFileName = inputfile;
	TFile *fr = new TFile(iFileName.Data());
	if ( !fr->IsOpen() ) {
		Error("CalculateCalibParameters", "File %s was not opened and won't be processed", iFileName.Data());
		return 0;
	}
	TTree *tr = (TTree*)fr->Get(treename);
	if (!tr) {
		Error("CalculateCalibParameters", "Tree %s was not found in file %s", treename, iFileName.Data());
		return 0;
	}


	//promenne potrebne pro fitovani: presunout nize
	//pohlidat delete
	TF1 	*calFunction = new TF1("calib", "pol1", 0, 1000);	//predelat jako lokalni promennou fce (nebo snad tridy?)
	TGraph 	*calGraph = new TGraph(kRaNOPEAKS, fPeak, fEnergy); //lokalni promenna, dohodit pocet vstupu pomoci parametru

	TString	detectorChannel;
	TString	histName;
	TString histTitle;
	TString fillCommand;
	TString	fillCondition;
	Int_t	fitControl = 0;

	//predelat nazvy histogramu
	//zrusit cyklus, napsat jako fci
	//raw data histogram filling
	TH1I *hRaw = 0;
	TH1F *hEnergy = 0;

	TRandom3 ranGen(1);

	//outputfile with calibrated spectra
	if ( (lowersubaddress == 0) && (uppersubaddress == ADDRESSNUMBER-1) ) { oFileName.Form("%s[]E.root", block); }
	else {
		if (lowersubaddress == uppersubaddress) { oFileName.Form("%s[%d]E.root", block, lowersubaddress); }
		else { oFileName.Form("%s[%d-%d]E.root", block, lowersubaddress, uppersubaddress); }

	}
	TFile *fw = new TFile(oFileName.Data(), "RECREATE");
	if (fw->IsOpen() == 0) {
		Error("CalculateCalibParameters", "File %s was not created and won't be processed\n\n", oFileName.Data());
		return 1;
	}

	for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) {
		printf("\n\n");
		Info("CalculateCalibParameters", "Calculating calibration parameters for detector channel %s[%d].", block, i);
		//TH1I object preparing
		hRaw = new TH1I("name", "title", 4096, 0, 4095);					//nastavovat hranice histogramu podle parametru fce
		detectorChannel.Form("%s[%d]", block, i);
		histName.Form("Hist%s[%d]", block, i);
		hRaw->SetName(histName.Data());
		fillCommand.Form("%s >> %s", detectorChannel.Data(), hRaw->GetName());
		fillCondition.Form("%s > %d && %s < %d",
				detectorChannel.Data(), lowerchannel, detectorChannel.Data(), upperchannel);
		//filling from the .root raw data file and content arrangement
		tr->Draw(fillCommand.Data(), fillCondition.Data(), "goff");

		//spectrum analysis
		fitControl = PeaksFitting(hRaw, "Q", fFitMinSigma);
		Info("CalculateCalibParameters", "Value of fitControl is: %d", fitControl);		//ok

		//incorrectly treated spectrum output
		if (fitControl != 0 && fCalInformation->IsOpen()) {															//ok
			outcalfile << setw(39) << fitControl << endl;
			fCalInformation->cd();
			hRaw->SetLineColor(2);	//red
			fCalInformation->cd();
			hRaw->Write();
			continue;
		}//if

		//correctly treated spectrum saving
		if (fCalInformation->IsOpen()) {
			fCalInformation->cd();
			hRaw->SetLineColor(kGreen+3);	//green
			hRaw->SetFillColor(kGreen+1);
			hRaw->Write();
		}//if

		//calibration parameters calculation																//ok
		for (Int_t j = 0; j < kRaNOPEAKS; j++) {		//delat podle poctu zkoumanych piku
			calGraph->SetPoint(j, fPeak[j], fEnergy[j]);  //calibration graph filling
		}//for
		calGraph->Fit(calFunction, "Q", "goff", 0, 4096);  //omezit hlasitost fitovani, udelat volitelne, dodelat volby rozsahu
		outcalfile
			<< block << "\t"
			<< address << "\t"
			<< i << "\t"
			<< setprecision(4) << calFunction->GetParameter(1) << "\t"
			<< setprecision(4) << calFunction->GetParameter(0) << "\t\t"
			<< fitControl
			<< endl;
		fA[address][i] = calFunction->GetParameter(1);
		fB[address][i] = calFunction->GetParameter(0);


		//calibration of raw spectra using obtained parameters
		Info("CalculateCalibParameters", "Energy spectrum from address %s[%d] calculating", block, i);
		histName.Form("%sE", hRaw->GetName());
		histTitle.Form("%s: %s", iFileName.Data(), histName.Data());
		hEnergy = new TH1F(histName.Data(), histTitle.Data(), nEBins, 0., 10.);
		//			detectorChannel.Form("%s[%d]", block, i);
		//			hEnergy->SetName(histName.Data());

		for (Int_t j = lowerchannel; j < upperchannel; j++) {
			Int_t binCont = (Int_t)hRaw->GetBinContent(j);
			//				cout << j << ":\t" << hRaw->GetBinContent(j) << endl;
			for (Int_t k = 0; k < binCont; k++) {
				hEnergy->Fill( fA[address][i]*( j+ranGen.Uniform(-0.5, 0.5) ) + fB[address][i] );
//				cout << j << ":\t" << fA[address][i]*( j+ranGen.Uniform(-0.5, 0.5) ) + fB[address][i] << endl;
//				cout << j << ":\t" << fA[address][i] << endl;
			}
		}

		fw->cd();
		hEnergy->Write();

	}//for

	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	//zapis histogramu v MeV do souboru .root
	//
	//urcite by se to dalo hodit do samostatne fce
	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	/*Char_t	histTitle[40];

	//outputfile with calibrated spectra
	if ( (lowersubaddress == 0) && (uppersubaddress == ADDRESSNUMBER-1) ) { oFileName.Form("%s[]E.root", block); }
	else {
		if (lowersubaddress == uppersubaddress) { oFileName.Form("%s[%d]E.root", block, lowersubaddress); }
		else { oFileName.Form("%s[%d-%d]E.root", block, lowersubaddress, uppersubaddress); }

	}
	TFile *fw = new TFile(oFileName.Data(), "RECREATE");
	if (fw->IsOpen() == 0) {
		Error("CalculateCalibParameters", "File %s was not created and won't be processed\n\n", oFileName.Data());
		return 1;
	}

	TH1F *hEnergy = 0;

	TRandom3 ranGen(1);

	Long64_t nEntries = tr->GetEntries();

	for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) {
		Info("CalculateCalibParameters", "Calibration spectrum from address %s[%d]", block, i);
		histName.Form("Hist%s[%d]E", block, i);
		sprintf(histTitle, "%s: %s", iFileName.Data(), histName.Data());
		hEnergy = new TH1F(histName.Data(), histTitle, 10000, 0., 10.);
		detectorChannel.Form("%s[%d]", block, i);
		hEnergy->SetName(histName.Data());
		fillCommand.Form("%f*(%s+%f)+%f >> %s", fA[address][i], detectorChannel.Data(), ranGen.Uniform(-0.5, 0.5), fB[address][i], hEnergy->GetName());
//		fillCommand.Form("%f*(%s+%f)+%f >> %s", fA[address][i], detectorChannel.Data(), 0., fB[address][i], hEnergy->GetName());
		cout << fillCommand.Data() << endl;
		fillCondition.Form("%s > %f && %s < %f", detectorChannel.Data(), fLowerChannel, detectorChannel.Data(), fUpperChannel);
		//filling from the .root raw data file and content arrangement

		for (Int_t j = 0; j < nEntries; j++) {

		}

		tr->Draw(fillCommand.Data(), fillCondition.Data(), "goff");


		fw->cd();
		hEnergy->Write();
	}//for

	fw->Close();*/

	fw->Close();
	fr->Close();
	fCalInformation->Close();
	outcalfile.close();

	return 1;

}



void AculCalibration::ShowAnalyzedSpectra(const char *filename, TCanvas* fittedRawCanvas, Int_t xaxismin, Int_t xaxismax, Int_t subaddress)
{

	if ( subaddress > ADDRESSNUMBER ) {
		Error("ShowAnalyzedSpectra", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1);
		return;
	}

	if (!fittedRawCanvas) {
		Warning("ShowAnalyzedSpectra", "You have to assign TCanvas for fitted raw spectra drawing");
		return;
	}


	TFile *fr = new TFile(filename, "READ");
	if (!fr->IsOpen()) {
		cout << "File " << filename << " was not opened" << endl;
		return;
	}

	TList *histList;
	histList = fr->GetListOfKeys();
	Int_t listEntries = histList->GetEntries();
	TH1I *hDraw = 0;

	fittedRawCanvas->Clear();
//	fittedRawCanvas->SetFillColor(10);

	if ( (listEntries > 1) && (listEntries <= 8) ) {
		fittedRawCanvas->Divide(2, 4);
		fittedRawCanvas->SetFillColor(10);
	}
	if ( (listEntries > 8) && (listEntries <= 16) ) {
		fittedRawCanvas->Divide(4, 4);
		fittedRawCanvas->SetFillColor(10);
	}

	if (subaddress >= listEntries) {
		for (Int_t i = 0; i < listEntries; i++) {
			fittedRawCanvas->cd(i+1);
			fr->GetObject(histList->At(i)->GetName(), hDraw);
			if (hDraw) {
				hDraw->SetAxisRange(xaxismin, xaxismax, "X");
				hDraw->Draw();
				hDraw->SetDirectory(0);
//				if (fHAnalyzedList) {
//					fHAnalyzedList->Add(hDraw);
//				}
				fHAnalyzedList.Add(hDraw);
			}
		}//for
	}
	else {
		fr->GetObject(histList->At(subaddress)->GetName(), hDraw);
		if (hDraw) {
			hDraw->SetAxisRange(xaxismin, xaxismax, "X");
			hDraw->Draw();
			hDraw->SetDirectory(0);
			fHAnalyzedList.Add(hDraw);
		}
	}

	fr->Close();

	fittedRawCanvas->Update();

	return;

}

void AculCalibration::ShowEnergySpectra(const char *filename, TCanvas* energyCanvas, const Int_t subaddress, Option_t* option, Double_t xaxismin, Double_t xaxismax)
{
	if ( subaddress > ADDRESSNUMBER ) {
		Error("ShowEnergySpectra", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1);
		return;
	}

	if (!energyCanvas) {
		Warning("ShowEnergySpectra", "You have to assign TCanvas for fitted raw spectra drawing");
		return;
	}

	TString	opt = option;
	opt.ToLower();

	TFile *fr = new TFile(filename, "READ");
	if (!fr->IsOpen()) {
		cout << "File " << filename << " was not opened" << endl;
		return;
	}

	TList *histList;
	histList = fr->GetListOfKeys();
	Int_t listEntries = histList->GetEntries();
	TH1F *hDraw = 0;

	energyCanvas->Clear();
	energyCanvas->SetFillColor(10);
	if ( (listEntries > 1) && (listEntries <= 8) ) {
		energyCanvas->Divide(2, 4);
		energyCanvas->SetFillColor(10);
	}
	if ( (listEntries > 8) && (listEntries <= 16) ) {
		energyCanvas->Divide(4, 4);
		energyCanvas->SetFillColor(10);
	}


	if (subaddress >= listEntries) {
		if (opt.Contains("sum")) {
			energyCanvas->cd(0);
			for (Int_t i = 0; i < listEntries; i++) {
				fr->GetObject(histList->At(i)->GetName(), hDraw);
				if (hDraw) {
					hDraw->SetDirectory(0);
					//hDraw->SetAxisRange(xaxismin, xaxismax, "X");
					if (opt.Contains("c")) { hDraw->SetLineColor(i+1); }
					if (opt.Contains("c") && opt.Contains("+")) { hDraw->SetFillColor(i+1); }
//					fHEnergyStack->Add(hDraw);
					fHEnergyStack.Add(hDraw);
				}
			}
//			if (fHEnergyStack) {
//				if (opt.Contains("+")) { fHEnergyStack->Draw(); }
//				else { fHEnergyStack->Draw("nostack"); }
//			}
			if (opt.Contains("+")) { fHEnergyStack.Draw(); }
			else { fHEnergyStack.Draw("nostack"); }
		}
		else {
			for (Int_t i = 0; i < listEntries; i++) {
				energyCanvas->cd(i+1);
				fr->GetObject(histList->At(i)->GetName(), hDraw);
				if (hDraw) {
					hDraw->SetAxisRange(xaxismin, xaxismax, "X");
					hDraw->Draw();
					hDraw->SetDirectory(0);
//					if (fHEnergyList) {
//						fHEnergyList->Add(hDraw);
//					}
					fHEnergyList.Add(hDraw);
				}
			}//for
		}//else
	}//if
	else {
		fr->GetObject(histList->At(subaddress)->GetName(), hDraw);
		energyCanvas->cd(0);
		if (hDraw) {
			hDraw->SetAxisRange(xaxismin, xaxismax, "X");
			hDraw->Draw();
			hDraw->SetDirectory(0);
//			if (fHEnergyList) {
//				fHEnergyList->Add(hDraw);
//			}
			fHEnergyList.Add(hDraw);
		}
	}//else

	fr->Close();

//	fFileName = filename;
//	fFileName.Resize(fFileName.Length() - 6);
//	fFileName.Append(".cal", 4);

	energyCanvas->Update();

	return;

}

Bool_t AculCalibration::AddCalFileToList(const char* calfilelist)
{
	//this function does not work at the moment
	//some problem with TString object fFileName

	TString fl = calfilelist;
	fl.ToLower();

	ofstream fw;
	fw.open(fl.Data(), ofstream::app);
	if (!fw.is_open()) {
		cout << "File " << fl.Data() << " was not opened" << endl;
		return kFALSE;
	}

//	fw << fFileName.Data() << endl;

	fw.close();

	return kTRUE;

}

void AculCalibration::ClearHistograms(Option_t* option)
{
	//clear THStack and TObjArray members
	//this function will be removed as soon as possible

	TString opt = option;
	opt.ToLower();

//	fHRawList->Add(hRead);
	fHRawList.Clear();
	fHAnalyzedList.Clear();
	fHEnergyList.Clear();
	fHEnergyStack.Clear();

	return;

}

void AculCalibration::MakeCalibrationFile(Char_t* calibrationfile, Char_t *calfilelist)
{
	//calibrationfile: file with calibration parameters to be created
	//calfilelist: file containing list of existing text files with calibration parameters

	ifstream calListR;
	calListR.open(calfilelist);
	if( !calListR.is_open() ) {
		cout << "File with list of calibration files was not opened" << endl;
		return;
	}

	//asi fce Reset()
	for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
		for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
			fA[i][j] = 0;
			fB[i][j] = 0;
		}
	}

	const Int_t lineLength = 100;
	char line[lineLength];
	char filename[50];
//	ifstream calFileR;
	Int_t	crate, i, j, id;
//	int	crate, i, j, id;
	char	cA[40], cB[40], cSigma[40];

	while (!calListR.eof()) {
		calListR.getline(line, lineLength);
//		cout << line << endl;
		sscanf(line, "%s", filename);
		cout << filename << endl;
		ifstream calFileR;
		calFileR.open(filename);
		if (calFileR.is_open()) {
			cout << filename << " processing" << endl;
			calFileR.seekg(0);
			while (!calFileR.eof()) {
//				cout << " in inner while" << endl;
				calFileR.getline(line, lineLength);
				sscanf(line, "%d %d %d %s %s %d %s", &crate, &i, &j, cA, cB, &id, cSigma);
//				cout << line << endl;
				if (id == 0) {
					fA[i][j] = static_cast<Double_t>(atof(cA));
					fB[i][j] = static_cast<Double_t>(atof(cB));
	//				fMeanSigma[i][j] = static_cast<Double_t>(atof(cSigma));
					cout << fA[i][j] << "\t" << fB[i][j] << endl;
				}//if
			}//while
			calFileR.close();
			cout << "calFileR was closed" << endl << endl;
		}//if
	}//while

	calListR.close();

	ofstream CalibFileW;
	CalibFileW.open(calibrationfile);
	if (!CalibFileW.is_open()) {
		cout << "Calibration file was not opened" << endl;
		return;
	}

	for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
		for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
			if (fA[i][j] != 0) {
				CalibFileW << std::right
					<< setw(2) << "3"
					<< setw(4) << i
					<< setw(4) << j
					<< setw(12) << fA[i][j]
					<< setw(12) << fB[i][j]
	//				<< setw(10) << fMeanSigma[i][j]
					                     << endl;
			}
		}
	}

	CalibFileW.close();

	return;

}

void AculCalibration::DeleteStacks(Option_t* option) {

	if (fCurrentHStack) {
		delete fCurrentHStack;
		fCurrentHStack = NULL;
	}

	fCurrentHistList.Delete();

	return;
}

void AculCalibration::Reset()
{
	//reset calibration parameters fA, fB, fC, fD to zero

	for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
		for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
			fA[i][j] = 0;
			fB[i][j] = 0;
			fC[i][j] = 0;
			fD[i][j] = 0;
	//		fMeanSigma[i][j] = 0;
		}
	}

	return;
}