From a0882328e0b8ebf52ffc269837b0f3745397cef1 Mon Sep 17 00:00:00 2001 From: Vratislav Chudoba Date: Wed, 5 Oct 2016 12:01:27 +0300 Subject: [PATCH] Function AculCalibration::CalibrateRawSpectra() added. --- AculData/AculCalibration.cpp | 171 +++++++++++++++++++++++++++++------ AculData/AculCalibration.h | 29 +++++- 2 files changed, 169 insertions(+), 31 deletions(-) diff --git a/AculData/AculCalibration.cpp b/AculData/AculCalibration.cpp index 23bff31..19ee6eb 100755 --- a/AculData/AculCalibration.cpp +++ b/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(atof(cA)); - fB[i][j] = static_cast(atof(cB)); + fAOld[i][j] = static_cast(atof(cA)); + fBOld[i][j] = static_cast(atof(cB)); // fMeanSigma[i][j] = static_cast(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; diff --git a/AculData/AculCalibration.h b/AculData/AculCalibration.h index 3631444..c1933de 100755 --- a/AculData/AculCalibration.h +++ b/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 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, -- 2.18.1