Try to commit all

AculData/AculCalibration.cpp

@@ -258,17 +258,17 @@ Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t s
 }
 
 
-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::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)
 {
@@ -545,7 +545,6 @@ void AculCalibration::ShowSpectra(const char* filename, TCanvas* rawCanvas, Opti
 			fr.GetObject(histList->At(i)->GetName(), hDraw);
 			if (hDraw) {
 				hDraw->SetDirectory(0);
-				//hDraw->SetAxisRange(xaxismax, xaxismin);	//nefunguje
 				fCurrentHistList.Add(hDraw);
 				fCurrentHStack->Add(hDraw);
 			}
@@ -813,13 +812,46 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
 
 	}//for
 
-	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-	//zapis histogramu v MeV do souboru .root
-	//
-	//urcite by se to dalo hodit do samostatne fce
-	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+	fw->Close();
+	fr->Close();
+	fCalInformation->Close();
+	outcalfile.close();
+
+	return 1;
+
+}
+
+void AculCalibration::CalibrateRawSpectra(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) {
+
+
+	//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;
+	}
+	TTree *tr = (TTree*)fr->Get(treename);
+	if (!tr) {
+		Error("CalculateCalibParameters", "Tree %s was not found in file %s", treename, iFileName.Data());
+		return;
+	}
+
+
+	TH1I *hRaw = 0;
+	TH1F *hEnergy = 0;
 
-	/*Char_t	histTitle[40];
+	TRandom3 ranGen(1);
+
+
+	TString	detectorChannel;
+	TString	histName;
+	TString histTitle;
+	TString fillCommand;
+	TString	fillCondition;
+	TString oFileName;
 
 	//outputfile with calibrated spectra
 	if ( (lowersubaddress == 0) && (uppersubaddress == ADDRESSNUMBER-1) ) { oFileName.Form("%s[]E.root", block); }
@@ -831,50 +863,86 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
 	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;
+		return;
 	}
 
-	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.);
+		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);
-		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);
+		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");
 
-		for (Int_t j = 0; j < nEntries; j++) {
+		//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
 
-		tr->Draw(fillCommand.Data(), fillCondition.Data(), "goff");
+		//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);
 
-		fw->cd();
-		hEnergy->Write();
-	}//for
 
-	fw->Close();*/
+		//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());
 
-	fw->Close();
-	fr->Close();
-	fCalInformation->Close();
-	outcalfile.close();
+		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;
+			}
+		}
 
-	return 1;
+		fw->cd();
+		hEnergy->Write();
 
-}
+	}//for
 
+}
 
 
 void AculCalibration::ShowAnalyzedSpectra(const char *filename, TCanvas* fittedRawCanvas, Int_t xaxismin, Int_t xaxismax, Int_t subaddress)
@@ -1048,6 +1116,10 @@ void AculCalibration::ShowEnergySpectra(const char *filename, TCanvas* energyCan
 
 }
 
+void AculCalibration::DivideCanvas(TCanvas *c1, Int_t inputs) {}
+
+
+
 Bool_t AculCalibration::AddCalFileToList(const char* calfilelist)
 {
 	//this function does not work at the moment

AculData/AculCalibration.h

@@ -154,9 +154,11 @@ public:
 	//  lowersubaddress: block subaddress
 	//  uppersubaddress: block subbaddress
 
-	Bool_t	EnergyPositions(const char* inputfile, const char* block,
-			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 CalibrateRawSpectra(const char* inputfile, const char* block, const Int_t address, const char* treename, Int_t lowerchannel = 0, Int_t upperchannel = 4095, Int_t nEBins = 1000, Int_t lowersubaddress = 0, Int_t uppersubaddress = ADDRESSNUMBER-1);
+
+//	Bool_t	EnergyPositions(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);
 
 	Int_t	PeaksFitting(TH1* hSpectrum, Option_t* option = "", Double_t sigmamin = 2);
 	//
@@ -167,7 +169,9 @@ public:
 	void	FillRawSpectraFile(const char* rawdatafile, const char* block, const char* treename, TCanvas* rawCanvas = NULL, Option_t *option = "", Int_t xaxismin = 0, Int_t xaxismax = 4096);
 
 	void	ShowRawSpectra(const char* filename, const Int_t block, TCanvas* rawCanvas = NULL, Int_t xaxismin = 0, Int_t xaxismax = 4096, /*TObjArray* histList = NULL,*/ const Int_t subaddress = 16);
-	void	ShowSpectra(const char* filename, TCanvas* rawCanvas = NULL, Option_t *option = "", Int_t xaxismin = 0, Int_t xaxismax = 4096, /*TObjArray* histList = NULL,*/ const Int_t subaddress = 16);
+	//this function is written for old format of raw data
+	//do not use this function!!!!!!!!!!!
+	void	ShowSpectra(const char* filename, TCanvas* rawCanvas = NULL, Option_t *option = "", Int_t xaxismin = 0, Int_t xaxismax = 4096, /*TObjArray* histList = NULL,*/ const Int_t subaddress = ADDRESSNUMBER);
 	void	ShowAnalyzedSpectra(const char* filename, TCanvas* fittedRawCanvas = NULL, Int_t xaxismin = 0, Int_t xaxismax = 4096, Int_t subaddress = 16);
 	//This function displays analyzed spectrum from a file, divides the canvas into a sufficient number of pads and displays
 	//spectrums of each  block subadress on the suitable pads or displays one  selected spectrum .
@@ -191,7 +195,7 @@ public:
 
 	void DivideCanvas(TCanvas *c1, Int_t inputs);
 
-//	void
+
 
 	//dodelat funkce TTree* Get...(...) pro raw, anal i E spectra
 	Bool_t	AddCalFileToList(const char* calfilelist = "CalFileList.log");

TELoss/TELoss.h

@@ -28,7 +28,11 @@ private:
 	double a, b;		//linear interpolation coeficients, y = a*x + b
 	int bsearch(int ntab, double *xtab, double x);
 	double aitken(int ntab, double *xtab, double *ytab, double x);
-	double aitken3(int ntab, double *xtab, double *ytab, double x);
+	double aitken3(int ntab, double *xtab, double *ytab, double x);
+	//====================================================================
+	//== Interpolation by 4 points
+	//====================================================================
+
 	Double_t linear(int ntab, Double_t *xtab, Double_t *ytab, double x);
 	Double_t linear(int ntab, Double_t *xtab, Double_t *ytab, Int_t i0, double x);
 
@@ -38,20 +42,58 @@ public:
 	virtual ~TELoss(void);
 	void SetEL(int _mel, double _den);
 	void SetDen(double _den) {den = _den;}
-	int  AddEL(double _zel, double _ael, double _wel=1.);
+	int  AddEL(double _zel, double _ael, double _wel=1.);
+	//===============================================================
+	//== Add one element
+	//===============================================================
 	void SetZP(double _zp, double _ap);
-	void SetEtab(int _ne, double e2);	//from 0 MeV to e2 MeV, _ne elements
-	void SetDeltaEtab(double r);		//pro kazdou tlostku bude zvlastni tabulka, nebo taky ne
+	void SetEtab(int _ne, double e2);	//from 0 MeV to e2 MeV, _ne elements
+	//===============================================================
+	//== Set list of energies and calculate R(E)
+	//===============================================================
+
+	void SetDeltaEtab(double r);		//pro kazdou tlostku bude zvlastni tabulka, nebo taky ne
+	//===============================================================
+	//== Calculate and set list of dE(E)
+	//===============================================================
+
 	double GetZ() { return zp; };
 	double GetA() { return ap; };
-	double GetE(double e0, double r);
-	double GetE0dE(double de, double r);	//de in MeV, r in microns
+	double GetE(double e0, double r);
+	//==================================================================
+	//== Calculate new energy using linear interpolation
+	//==================================================================
+
+	double GetE0dE(double de, double r);	//de in MeV, r in microns
+	//==================================================================
+	//== Calculate new energy from deltaE using linear interpolation
+	//==================================================================
+
 	double GetE0dE(double de);				//de in MeV
-	double GetE0(double e, double r);
-	double GetEold(double e0, double r);
-	double GetE0old(double e, double r);
-	double GetR(double e0, double e);
-	double GetRold(double e0, double e);
+	double GetE0(double e, double r);
+	//==================================================================
+	//== Calculate new energy using linear interpolation
+	//==================================================================
+
+	double GetEold(double e0, double r);
+	//==================================================================
+	//== Calculate new energy using aitken interpolation
+	//==================================================================
+
+	double GetE0old(double e, double r);
+	//==================================================================
+	//== Calculate new energy
+	//==================================================================
+
+	double GetR(double e0, double e);
+	//==================================================================
+	//== Calculate new energy using linear interpolation
+	//==================================================================
+
+	double GetRold(double e0, double e);
+	//==================================================================
+	//== Calculate new energy
+	//==================================================================
 	void PrintRE();
 	void PrintdEE();
 	void PrintREtoFile();

calibration1.cxx

 {
-	gSystem->Load("/home/dariak/Utilities/libAculData.so");
+	gSystem->Load("/home/dariak/AculUtils/libAculData.so");
 
 	AculCalibration cal;
-	cal.SetInputParameters();
-	cal.CalculateCalibParameters("clb01_0001.root", "SQ22", 22, "AnalysisxTree", 50, 700);
-
-	TCanvas *c1 = new TCanvas();
-
-//	cal->ShowRawSpectra("clb01_0001.root", 22, c1);
-//	cal.ShowSpectra("SQ22[].root", c1, "", 300, 800, 10);
-
-
-//	TFile filerootE("SQ22[]E.root");
-//	TFile filecal("SQ22[].cal");
-
-//	TDirectory dir("dir","dash");
-//	dir.Add("SQ22[].root");
-
-
+	cal.SetInputParameters(); //from .par
+	cal.CalculateCalibParameters("clb04_0001.root", "SQ22", 22, "AnalysisxTree", 100, 4095);
+//	    CalculateCalibParameters(const char* inputfile, const char* block, const Int_t address, const char* treename, Int_t lowerchannel = 0, Int_t upperchannel = 4095, Int_t nEBins = 1000, Int_t lowersubaddress = 0, Int_t uppersubaddress = ADDRESSNUMBER-1);
+
+	TFile fr("SQ22[]E.root");
+	cout << fr.GetListOfKeys()->GetEntries() << " histograms" << endl;
+	TList *histList = fr.GetListOfKeys();
+	TH1 *hWork = 0;
+	for (Int_t i = 0; i < 32; i++) 
+	{
+		fr.GetObject(histList->At(i)->GetName(), hWork);		
+		cal->PeaksFitting(hWork);
+	}
 }

calibration2.cxx

 {
-	gSystem->Load("/home/dariak/Utilities/libAculData.so");
+	gSystem->Load("/home/dariak/AculUtils/libAculData.so");
 
 	AculCalibration cal;
 	cal->SetInputParameters();
 
+	cal->SetCalibrationParameters("./run1/SQ22[].cal"); //takes calibration parameters 
+	cal->PrintCalibrationParameters(22, 22);
 
 	TCanvas *c1 = new TCanvas();
+	cal->CalibrateRawSpectra("clb03_0001.root", "SQ22", 22, "AnalysisxTree", 100, 800, 500); //takes data from raw file and calibrates it with obtained calibration parameters
 
 	TFile fr("SQ22[]E.root");
-//	TH1F *h1 = (TH1F*)fr.Get("HistSQ22[10]E");
-//	h1->Draw();
-
 	cout << fr.GetListOfKeys()->GetEntries() << " histograms" << endl;
 
 	TList *histList = fr.GetListOfKeys();
-	const Int_t noHist = histList->GetEntries();
-
-	const Double_t minE = 3;
-	const Double_t maxE = 8;
-	Int_t minBin, maxBin;
-
 	TH1 *hWork = 0;
 	c1->Clear();
 	c1->Divide(6, 6);
-	for (Int_t i = 0; i < 32; i++) {
+
+	for (Int_t i = 0; i < 32; i++) 
+	{
 		fr.GetObject(histList->At(i)->GetName(), hWork);		
 		c1->cd(i+1);
 		cal->PeaksFitting(hWork);
 		hWork->Draw();
 	}
 
+/*	c1->Divide(2, 2);
+
+	c1->cd(1);
+	fr.GetObject(histList->At(7)->GetName(), hWork);
+	cal->PeaksFitting(hWork);
+	hWork->Draw();
+
+	fr.GetObject(histList->At(8)->GetName(), hWork);
+	c1->cd(2);
+	cal->PeaksFitting(hWork);
+	hWork->Draw();
+
+	fr.GetObject(histList->At(23)->GetName(), hWork);
+	c1->cd(3);
+	cal->PeaksFitting(hWork);
+	hWork->Draw();
+
+	fr.GetObject(histList->At(24)->GetName(), hWork);
+	c1->cd(4);
+	cal->PeaksFitting(hWork);
+	hWork->Draw();
+	*/	
+
 }

parforcal.par

@@ -6,10 +6,10 @@
 6.002	E3	//in MeV
 7.687	E4	//in MeV
 100	lowerchannel	//in channels
-4096	upperchannel	//in channels
-0.35	lowerpeakhight	//in relative units
+800	upperchannel	//in channels
+0.5	lowerpeakhight	//in relative units
 0.5	upperpeakhight	//in relative units
 0.1	peakpositiontolerance	//in relative units
 2	fitfunctionlinewidth		//integer 1 - 10
-5	minfitsigma	//minimal sigma of the peaks to be fitted
-0.4	fithightthreshold	//
+2	minfitsigma	//minimal sigma of the peaks to be fitted
+0.3	fithightthreshold	//