AculData/AculCalibCsI.cpp

@@ -125,16 +125,6 @@ void AculCalibCsI::PrintCuts() {
 //	return;
 	if (fCuts) printf("Cuts loaded from file \"%s\"\n", fCuts->GetName());
 
-	/*for (Int_t i = 0; i < NOCALFILES; i++) {
-		if (cTOF[i]) {
-			printf("TOF cut No. %d; Name: \"%s\"\n", i, cutsCol[i]->GetName());
-			cout << "asdjasd" << endl;
-		}
-		else {
-			printf("TOF cut No. %d was not loaded. Maximal number of cuts is %d\n", i, NOCALFILES);
-		}
-	}//for */
-
 	TCutG *curCut = 0;
 	//cutsCol
 
@@ -147,19 +137,6 @@ void AculCalibCsI::PrintCuts() {
 		}*/
 	}
 
-/*	for (Int_t i = 0; i < NOCALFILES; i++) {
-		if (cQCD[i]) {
-			printf("QCD cut No. %d; Name: \"%s\"\n", i, cQCD[i]->GetName());
-//			cout << "kasjhd" << endl;
-		} else {
-//			cout << "bhajskd" << endl;
-			printf("QCD cut No. %d was not loaded. Maximal number of cuts is %d\n", i, NOCALFILES);
-//			cout << "vsdfjks" << endl;
-		}
-	}
-
-	printf("AculCalibCsI::PrintCuts: End of function.\n");
-*/
 	return;
 }
 
@@ -256,11 +233,11 @@ void AculCalibCsI::DrawBeam(TCanvas *canvas, Int_t files, const char* variable)
 		con.Form("%s[5]>200", variable);
 		curTree->Draw(var.Data(), con.Data(), "cont");
 		for (Int_t j = 0; j < energyPoints; j++) {		
-//			if ( cutsCol.At(files + j) ) {
-//				curCutG = (TCutG*)cutsCol.At(files + j);
-//				curCutG->Draw("same");
+			if ( cutsCol.At(files + j) ) {
+				curCutG = (TCutG*)cutsCol.At(files + j);
+				curCutG->Draw("same");
 //				printf("AculCalibCsI::DrawBeam: cQCD cut No. %d cannot be drawn, need to repair this function.\n", j);
-//			}
+			}
 		}
 		canvas->Update();
 	}
@@ -365,14 +342,14 @@ void AculCalibCsI::GetPeakMean(const char* variable, Int_t tree, Int_t energy, T
 		con.Form("%s[%d]>%d && %s", variable, channel, lowRange, cutsCol.FindObject(beamcut)->GetName());
 		canvas->cd(i+1);
 		hname.Form("hfull[%d][%d]", tree, i);
-		hfull[tree][i] = new TH1I(hname.Data(), "title", nbins, 0, 4096);
+		hfull[tree][i] = new TH1I(hname.Data(), "title1", nbins, 0, 4096);
 		curTree->SetLineColor(1);
 		curTree->Draw(var.Data(), con.Data());
 
 		var.Form("%s[%d]>>hcut[%d][%d]", variable, channel, tree, i);
 		con.Form("%s[%d]>%d && %s[%d]<%d && %s", variable, channel, peakMin[energy][i], variable, channel, peakMax[energy][i], cutsCol.FindObject(beamcut)->GetName());
 		hname.Form("hcut[%d][%d]", tree, i);
-		hcut[tree][i] = new TH1I(hname.Data(), "title", nbins, 0, 4096);
+		hcut[tree][i] = new TH1I(hname.Data(), "title2", nbins, 0, 4096);
 		hcut[tree][i]->SetLineColor(3);
 		curTree->Draw(var.Data(), con.Data(), "same");
 
@@ -395,7 +372,7 @@ void AculCalibCsI::Calibrate(TCanvas *canvas, Bool_t savefile, const char* filen
 	canvas->Divide(4,4);
 
 //	cout << alphas2.GetSize()+1 << endl;
-	cout << energyPoints+1 << endl;
+	cout <<" energyPoints+1 "<<energyPoints+1<< endl;
 
 	TString gName;
 	TString gTitle;
@@ -432,7 +409,7 @@ void AculCalibCsI::FillGraph(TGraphErrors *g, Int_t npoints, Double_t *energies,
 	TString opt = option;
 
 	//all available energy points and (0,0)
-	g->SetPoint(0, 0., 0.);
+	g->SetPoint(0, 0., 0.);	//(point number, x, y)
 	for (Int_t i = 0; i < npoints-1; i++) {
 //		g->SetPoint(i+1, energies[i], mean[i][graphNumber]);
 //		g->SetPointError(i+1, 0, meanRMS[i][graphNumber]);

AculData/AculCalibration.cpp

@@ -26,8 +26,8 @@ AculCalibration::AculCalibration() : fEnergy(0), fEnergyInput(0), fA(0), fB(0),
 //	todo: change size of fA and fB in some other place
 	fA.Set(32);
 	fB.Set(32);
-	fEnergy.Set(4);
-	fEnergyInput.Set(4);
+//	fEnergy.Set(4);
+//	fEnergyInput.Set(4);
 
 	kRaNOPEAKS = 0;
 	fLowerPeakRelativeHight = 0.;
@@ -172,7 +172,6 @@ Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t s
 		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);
 
@@ -259,7 +258,8 @@ Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t s
 	TMath::Sort(peaksNumber, peak, j, kFALSE);
 	fPeak.Set(peaksNumber);
 	for (Int_t i = 0; i < peaksNumber; i++) {
-		fPeak[i] = peak[j[i]];
+		fPeak[i] = peak[j[i]];
+		//printf("\tPeak peak\t%f\n", fPeak[i]);
 	}
 
 	if (!opt.Contains("q") || opt.Contains("v")) {
@@ -273,13 +273,14 @@ Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t s
 	//	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")) {
+	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])) ) ) {
+			//printf("\tPeaksFitt fEnergy\t%f\n", fEnergy[i]);		
+			if (fCalInformation /* && opt.Contains("writebad")*/) {
 				fCalInformation->cd();
 				hSpectrum->Write();
 			}
-			return 2;
+			//return 2;*/
 		}
 	}//for
 
@@ -759,7 +760,7 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
 		tr->Draw(fillCommand.Data(), fillCondition.Data(), "goff");
 
 		//spectrum analysis
-		fitControl = PeaksFitting(hRaw, "Q", fFitMinSigma);
+		fitControl = PeaksFitting(hRaw, "Q V v", fFitMinSigma);
 		Info("CalculateCalibParameters", "Value of fitControl is: %d", fitControl);		//ok
 
 		//incorrectly treated spectrum output
@@ -782,7 +783,8 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
 
 		//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
+			calGraph->SetPoint(j, fPeak[j], fEnergy[j]);  //calibration graph filling
+			printf("\tPeak\t%f and energy\t%f\n", fPeak[j], fEnergy[j]);
 		}//for
 		calGraph->Fit(calFunction, "Q", "goff", 0, 4096);  //omezit hlasitost fitovani, udelat volitelne, dodelat volby rozsahu
 		outcalfile
@@ -1101,13 +1103,77 @@ void AculCalibration::FindEnergyPeaks(TCanvas *c1, const char* ifile, const char
 		c1->cd(i+1);
 		PeaksFitting(hWork);
 		hWork->Draw();
+		ofile<<i<<"\t";
+		for(Int_t j=0; j<kRaNOPEAKS; j++) {
+			ofile << fPeak[j] <<"\t";
+		}
+		ofile<<endl;
+
+	}
+
+	ofile.close();
+}
+
+void AculCalibration::FindAverageEnergies(const char* ifile, const char* outfile) {
+
+	TString iFile = ifile;
+	TFile *fr = new TFile(iFile.Data());
+	if ( !fr->IsOpen() ) {
+		Error("FindAverageEnergies", "File %s was not opened and won't be processed.", iFile.Data());
+		return;
+	}
 
-		ofile << i <<"\t"<< fPeak[0] <<"\t"<< fPeak[1] <<"\t"<< fPeak[2] <<"\t"<< fPeak[3] <<endl;
+	TList *histList = fr->GetListOfKeys();
+	Info("FindAverageEnergies", "%d keys found in file %s.", histList->GetEntries(), fr->GetName());
 
+	//creation of output text file with average values of peak energies in MeV
+	TString workFile = outfile;
+	ofstream ofile;
+	ofile.open(workFile.Data());
+	if (!ofile.is_open()) {
+		Error("PeaksFitting", "Output file %s was not opened", workFile.Data());
+		return;
 	}
 
+	TH1 *hWork = 0;
+	Double_t hArray[histList->GetEntries()][kRaNOPEAKS];
+	//TString hSumName;
+	Double_t hSumE1 = 0.;
+	Double_t hAvrE1 = 0.;
+	Double_t hSumE2 = 0.;
+	Double_t hAvrE2 = 0.;
+	Double_t hSumE3 = 0.;
+	Double_t hAvrE3 = 0.;
+	Double_t hSumE4 = 0.;
+	Double_t hAvrE4 = 0.;
+//	c1->Clear();
+//	c1->Divide(6, 6);
+
+	for (Int_t i = 0; i < histList->GetEntries(); i++) {
+		fr->GetObject(histList->At(i)->GetName(), hWork);
+		PeaksFitting(hWork);
+		for(Int_t j = 0; j < kRaNOPEAKS; j++) {
+			hArray[i][j] = fPeak[j];
+			if(fPeak[j]==0.){
+				Error("FindAverageEnergies", "No peak in channel %i !", histList->GetEntries());			
+			}
+			//hSumName.Form("hSumE%i",j);
+		}
+
+		hSumE1 += hArray[i][0];
+		hSumE2 += hArray[i][1];
+		hSumE3 += hArray[i][2];
+		hSumE4 += hArray[i][3];
+//		std::cout<<"i "<<i<<" hSumE1 "<<hSumE1<<std::endl;
+	}
+	hAvrE1 = hSumE1/histList->GetEntries();
+	hAvrE2 = hSumE2/histList->GetEntries();
+	hAvrE3 = hSumE3/histList->GetEntries();
+	hAvrE4 = hSumE4/histList->GetEntries();
+	ofile <<"Average energies are:\t"<<hAvrE1<<"\t"<<hAvrE2<<"\t"<<hAvrE3<<"\t"<<hAvrE4<<std::endl;
 	ofile.close();
-}
+}
+
 
 void AculCalibration::ShowAnalyzedSpectra(const char *filename, TCanvas* fittedRawCanvas, Int_t xaxismin, Int_t xaxismax, Int_t subaddress)
 {
@@ -1334,7 +1400,9 @@ void AculCalibration::SetInputParameters() {
 
 		if ( strcmp(identificator, "noPeaks") == 0 ) {
 			kRaNOPEAKS = static_cast<Int_t>(atoi(parameter));
-			fEnergyInput.Set(kRaNOPEAKS);
+			fEnergyInput.Set(kRaNOPEAKS);
+			fEnergy.Set(kRaNOPEAKS);
+			fPeak.Set(kRaNOPEAKS);
 			for (Int_t i = 0; i < kRaNOPEAKS; i++) {
 				fipr.getline(line, lineLength);
 				sscanf(line, "%s", parameter);

AculData/AculCalibration.h

@@ -176,10 +176,17 @@ public:
 
 	void CalibrateRawSpectra(const char* inputfile, const char* block, const char* treename, Int_t lowerchannel = 0, Int_t upperchannel = 4095, Int_t nEBins = 1000, Int_t lowersubaddress = 0, Int_t uppersubaddress = ADDRESSNUMBER-1);
 
-	void FindEnergyPeaks(TCanvas *c1, const char* ifile, const char* outfile);
-//	const char* inputfile, const char* block,
-//			const Int_t address, const Int_t address, const char* treename, Int_t lowerchannel = 0,
-//			Int_t upperchannel = 4095, Int_t lowersubaddress = 0, Int_t uppersubaddress = ADDRESSNUMBER-1);
+	void FindEnergyPeaks(TCanvas *c1, const char* ifile, const char* outfile); 
+	// Outputs calibrated energies for each channel in txt file 
+	//
+	// 
+	//  
+
+	void FindAverageEnergies(const char* ifile, const char* outfile);
+	// Outputs average values of calibrated energies for the whole detector in txt file
+	//
+	// 
+	//  
 
 	Int_t	PeaksFitting(TH1* hSpectrum, Option_t* option = "", Double_t sigmamin = 2);
 	//return values:

macros/calibration_CsI/calibrationSQ13Alpha.cxx

 {
 //	gSystem->Load("CsICalib_C.so");
-	gSystem->Load("../../libAculData.so");
+	gSystem->Load("/home/dariak/AculUtils/libAculData.so");
 
 	//CsICalib cal;
 	AculCalibCsI cal("parameters/SQ13Alpha.par");
@@ -13,14 +13,17 @@
 	TCanvas *c2 = new TCanvas("c2", "Plain");
 	TCanvas *c3 = new TCanvas("c3", "Plain");
 	TCanvas *c4 = new TCanvas("c4", "Plain");
+//	TCanvas *myCanv = new TCanvas("myCanv", "Let's see what we see");
 
-	cal.PrintTrees();
-	cal.PrintFiles();
-	cal.PrintCuts();
+//	cal.PrintTrees();
+//	cal.PrintFiles();
+//	cal.PrintCuts();
+//	cal.PrintPeakRanges();
 
 //	cal.ReadParFile("ranges.dat");
-//	cal.PrintParameters();
-
+//	cal.PrintParameters("v");
+//	cal.PrintTrees();
+//	cal.DrawVariable("SQ13", 0, c1);
 
 
 //	cal.DrawVariable("SQ13", 1, c1);
@@ -28,10 +31,10 @@
 //	cal.DrawVariable("TDC", 2, c2);
 //	cal.DrawVariable("QDC", 3, c3);
 
-//	cal.DrawBeam(c4, 4, "SQ13");
+//	cal.DrawBeam(myCanv, 4, "SQ13");
 //return;
-//	cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16");
-//	cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp");
+//	cal.DrawVariableCut("SQ13", 0, c2, "cutsSQ13Alpha16");
+//	cal.DrawVariableCut("SQ13", 0, c2, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp");
 //return;
 //	cal.DrawVariableCut("SQ13", 1, c2, "cutSQ13Alpha21", "cutSQ13Alpha21Amp");
 //	cal.DrawVariableCut("SQ13", 2, c3, "cutSQ13Alpha26", "cutSQ13Alpha26Amp");
@@ -50,7 +53,7 @@
 //	cal.SaveClbGraphs("gSQ13Alpha2points.root", "RECREATE");
 //	cal.SaveClbGraphs("gSQ13AlphaAlt.root", "RECREATE");
 
-//	cal.WriteClbParameters("calSQ13Alpha.clb");
+	cal.WriteClbParameters("calSQ13Alpha.clb");
 //	cal.WriteClbParameters("calSQ13Alpha2points.clb");
 //	cal.WriteClbParameters("calSQ13AlphaAlt.clb");
 

macros/calibration_CsI/calibrationSQ23Alpha.cxx

 {
-	gSystem->Load("CsICalib_C.so");
+	gSystem->Load("/home/dariak/AculUtils/libAculData.so");
 
 //	CsICalib cal;
-	CsICalib cal("parameters/SQ23Alpha.par");
+	AculCalibCsI cal("parameters/SQ23Alpha.par");
 
 
 //	cal.PrintCuts();
@@ -13,8 +13,8 @@
 	TCanvas *c2 = new TCanvas("c2", "Plain");
 	TCanvas *c3 = new TCanvas("c3", "Plain");
 	TCanvas *c4 = new TCanvas("c4", "Plain");
-
-	cal.DrawBeam(c4, 4, "SQ23");
+	cal.PrintParameters();
+//	cal.DrawBeam(c4, 4, "SQ23");
 
 return;