Function AculCalibration::CalibrateRawSpectra() added.

AculData/AculCalibration.cpp

@@ -17,7 +17,7 @@
 
 ClassImp(AculCalibration);
 
-AculCalibration::AculCalibration()
+AculCalibration::AculCalibration() : fA(0), fB(0)
 {
 	//default constructor
 
@@ -43,7 +43,7 @@ AculCalibration::AculCalibration()
 
 AculCalibration::AculCalibration(const char* parfile)
 {
-	//constructor which fills fA, fB, fC, fD from file parfile
+	//constructor which fills fAOld, fBOld, fC, fD from file parfile
 
 	fCurrentHStack = NULL;
 	fCurrentHistList.IsOwner();
@@ -393,8 +393,8 @@ Bool_t AculCalibration::SetCalibrationParameters(const char* calparfile)
 		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);
+			fAOld[i][j] = atof(cA);
+			fBOld[i][j] = atof(cB);
 			fC[i][j] = atof(cC);
 			fD[i][j] = atof(cD);
 		}
@@ -408,7 +408,7 @@ void AculCalibration::PrintCalibrationParameters(const Int_t blockmin, const Int
 {
 	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;
+			cout << "C3[" << i << "][" << j << "]" << setw(10) << fAOld[i][j] << setw(10) << fBOld[i][j] << setw(10) << fC[i][j] << setw(10) << fD[i][j] << endl;
 		}
 	}
 }
@@ -785,8 +785,8 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
 			<< setprecision(4) << calFunction->GetParameter(0) << "\t\t"
 			<< fitControl
 			<< endl;
-		fA[address][i] = calFunction->GetParameter(1);
-		fB[address][i] = calFunction->GetParameter(0);
+		fAOld[address][i] = calFunction->GetParameter(1);
+		fBOld[address][i] = calFunction->GetParameter(0);
 
 
 		//calibration of raw spectra using obtained parameters
@@ -801,9 +801,9 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
 			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;
+				hEnergy->Fill( fAOld[address][i]*( j+ranGen.Uniform(-0.5, 0.5) ) + fBOld[address][i] );
+//				cout << j << ":\t" << fAOld[address][i]*( j+ranGen.Uniform(-0.5, 0.5) ) + fBOld[address][i] << endl;
+//				cout << j << ":\t" << fAOld[address][i] << endl;
 			}
 		}
 
@@ -821,6 +821,125 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
 
 }
 
+void AculCalibration::CalibrateRawSpectra() {
+
+	//function parameters:
+	const char* iFileName = "clb01_0001.root";
+	const char* treeName = "AnalysisxTree";
+	const char* branchName = "LiEvent.SQ22[32]";
+	const Int_t address = 22;
+	const Int_t lowerElement = 0;
+	const Int_t upperElement = 32;
+	const Int_t lowerChannel = 100, upperChannel = 4096;
+	const Int_t nEBins = 1000;
+	//optional:
+	Long64_t nentries = 0;
+
+
+	//function itself:
+	TString iFile = iFileName;
+	TFile fr( iFile.Data() );
+	if ( !fr .IsOpen() ) {
+		Error("CalibrateRawSpectra", "File %s was not opened and won't be processed", iFile.Data());
+		return;
+	}
+
+	TString tName = treeName;
+	TTree *tr = (TTree*)fr.Get(tName.Data());
+	if (!tr) {
+		Error("CalibrateRawSpectra", "Tree %s was not found in file %s", tName.Data(), iFile.Data());
+		return;
+	}
+	tr->SetMakeClass(1);
+
+
+	//!!!!!!!!!!!!!!!!!!!!!!!!!!!
+	UShort_t variable[32];
+	//!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+	TString bName = branchName;
+	tr->SetBranchAddress(bName.Data(), variable);
+
+	if (nentries == 0) nentries = tr->GetEntries();
+	Info("CalibrateRawSpectra", "%lld entries from tree %s will be read.", nentries, tr->GetName());
+
+	//make histogram to fill
+	const Int_t noElements = upperElement - lowerElement +1;	//number of treated detector elements
+	TH1F *hEnergy[noElements];
+	for (Int_t i = 0; i < noElements; i++) {
+		hEnergy[i] = 0;
+	}
+	TString	histName;
+	TString histTitle;
+	for (Int_t i = 0; i < noElements; i++) {
+		histName.Form("%sE%d", bName.Data(), i+lowerElement);
+		histTitle.Form("%s: %s", iFile.Data(), bName.Data());
+		hEnergy[i] = new TH1F(histName.Data(), histTitle.Data(), nEBins, 0., 10.);
+	}
+
+	TRandom3 ranGen(1);
+
+	for (Long64_t j = 0; j < nentries; j++) {
+		tr->GetEntry(j);
+
+		for (Int_t i = lowerElement; i <= upperElement; 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");
+
+			//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());
+			//			detectorChannel.Form("%s[%d]", block, i);
+			//			hEnergy->SetName(histName.Data());
+
+			if (variable[i] > lowerChannel && variable[i] < upperChannel) {
+				hEnergy[i-lowerElement]->Fill( fAOld[address][i]*( variable[i]+ranGen.Uniform(-0.5, 0.5) ) + fBOld[address][i] );
+			}
+
+
+
+		}//for subaddresses
+
+	}//for entries
+
+	fr.Close();
+
+	TString oFileName;
+	//outputfile with calibrated spectra
+	if ( (lowerElement == 0) && (upperElement == ADDRESSNUMBER-1) ) { oFileName.Form("%sE.root", bName.Data()); }
+	else {
+		if (lowerElement == upperElement) { oFileName.Form("%s[%d]E.root", bName.Data(), lowerElement); }
+		else { oFileName.Form("%s[%d-%d]E.root", bName.Data(), lowerElement, upperElement); }
+
+	}
+
+	TFile fw(oFileName.Data(), "RECREATE");
+	if (fw.IsOpen() == 0) {
+		Error("CalculateCalibParameters", "File %s was not created and won't be processed\n\n", oFileName.Data());
+		return;
+	}
+
+	fw.cd();
+	for (Int_t i = 0; i < noElements; i++) {
+		hEnergy[i]->Write();
+	}
+
+	fw.Close();
+
+	return;
+}
+
 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) {
@@ -915,8 +1034,8 @@ void AculCalibration::CalibrateRawSpectra(const char* inputfile, const char* blo
 //		<< setprecision(4) << calFunction->GetParameter(0) << "\t\t"
 //		<< fitControl
 //		<< endl;
-//		fA[address][i] = calFunction->GetParameter(1);
-//		fB[address][i] = calFunction->GetParameter(0);
+//		fAOld[address][i] = calFunction->GetParameter(1);
+//		fBOld[address][i] = calFunction->GetParameter(0);
 
 
 		//calibration of raw spectra using obtained parameters
@@ -931,9 +1050,9 @@ void AculCalibration::CalibrateRawSpectra(const char* inputfile, const char* blo
 			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;
+				hEnergy->Fill( fAOld[address][i]*( j+ranGen.Uniform(-0.5, 0.5) ) + fBOld[address][i] );
+				//				cout << j << ":\t" << fAOld[address][i]*( j+ranGen.Uniform(-0.5, 0.5) ) + fBOld[address][i] << endl;
+				//				cout << j << ":\t" << fAOld[address][i] << endl;
 			}
 		}
 
@@ -1176,8 +1295,8 @@ void AculCalibration::MakeCalibrationFile(Char_t* calibrationfile, Char_t *calfi
 	//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;
+			fAOld[i][j] = 0;
+			fBOld[i][j] = 0;
 		}
 	}
 
@@ -1205,10 +1324,10 @@ void AculCalibration::MakeCalibrationFile(Char_t* calibrationfile, Char_t *calfi
 				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));
+					fAOld[i][j] = static_cast<Double_t>(atof(cA));
+					fBOld[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;
+					cout << fAOld[i][j] << "\t" << fBOld[i][j] << endl;
 				}//if
 			}//while
 			calFileR.close();
@@ -1227,13 +1346,13 @@ void AculCalibration::MakeCalibrationFile(Char_t* calibrationfile, Char_t *calfi
 
 	for (Int_t i = 0; i < BLOCKSNUMBER; i++) {
 		for (Int_t j = 0; j < ADDRESSNUMBER; j++) {
-			if (fA[i][j] != 0) {
+			if (fAOld[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(12) << fAOld[i][j]
+					<< setw(12) << fBOld[i][j]
 	//				<< setw(10) << fMeanSigma[i][j]
 					                     << endl;
 			}
@@ -1260,12 +1379,12 @@ void AculCalibration::DeleteStacks(Option_t* option) {
 
 void AculCalibration::Reset()
 {
-	//reset calibration parameters fA, fB, fC, fD to zero
+	//reset calibration parameters fAOld, fBOld, 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;
+			fAOld[i][j] = 0;
+			fBOld[i][j] = 0;
 			fC[i][j] = 0;
 			fD[i][j] = 0;
 	//		fMeanSigma[i][j] = 0;

AculData/AculCalibration.h

@@ -75,12 +75,14 @@ public:
 	Double_t	fFitMinSigma;				//pouziva se, private
 	Double_t	fFitPeakThreshold;			//pouziva se, private, prozkoumat, k cemu vlastne slouzi ve fci ShowPeaks, popremyslet o vhodnem prednastaveni v konstruktoru
 
-	//tyto promenne jsou smyslem tridy, mely by tedy byt snad jako jedine verejne
-	Double_t	fA[BLOCKSNUMBER][ADDRESSNUMBER];	//kalibracni parametry, f(x) = fA*x + fB
-	Double_t	fB[BLOCKSNUMBER][ADDRESSNUMBER];	//kalibracni parametry, f(x) = fA*x + fB
+	//these variables are the main for the class
+	Double_t	fAOld[BLOCKSNUMBER][ADDRESSNUMBER];	//calibration parameter, f(x) = fAOld*x + fBOld
+	Double_t	fBOld[BLOCKSNUMBER][ADDRESSNUMBER];	//calibration parameter, f(x) = fAOld*x + fBOld
 	Double_t	fC[BLOCKSNUMBER][ADDRESSNUMBER];	//treti kalibracni parametr, jine zavislosti nez pol1
 	Double_t	fD[BLOCKSNUMBER][ADDRESSNUMBER];	//ctvrty kalibracni parametr
 
+	TArrayD fA;
+	TArrayD fB;
 
 
 
@@ -126,7 +128,7 @@ public:
 	Bool_t	SetCalibrationParameters(const char* calparfile);
 	//Loads the file with calibration parameters
 	//
-	//  calparfile: file containing calibration parameters in format: crate number, address, subaddress, fA, fB, fC, fD
+	//  calparfile: file containing calibration parameters in format: crate number, address, subaddress, fAOld, fBOld, fC, fD
 	//
 	//	allowed comment characters: *, #, %, //
 
@@ -134,7 +136,7 @@ public:
 	//I hope this is selfunderstanding function
 
 	void	PrintCalibrationParameters(const Int_t blockmin = 1, const Int_t blockmax = BLOCKSNUMBER - 1);
-	//Print calibration parameters fA, fB, fC, fD which are currently saved in object AculCalibration
+	//Print calibration parameters fAOld, fBOld, fC, fD which are currently saved in object AculCalibration
 	//
 	//  blockmin: minimum  block data are displayed for
 	//  blockmax: maximum  block data are displayed for
@@ -154,6 +156,23 @@ public:
 	//  lowersubaddress: block subaddress
 	//  uppersubaddress: block subbaddress
 
+	void CalibrateRawSpectra();
+
+	UShort_t SomeFunction(UShort_t x, int len);
+
+	template<typename T> T SomeFunction(T x, int len)
+	{
+		/*T max = x[0];
+		for(int i = 1; i < len; i++)
+			if(max < x[i])
+				max = x[i];*/
+
+//		return max;
+		return x;
+	}
+
+
+
 	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,