From ee5da75836888f1bcfc5689c38b783ed46fb533b Mon Sep 17 00:00:00 2001 From: "Kostyleva D.A" Date: Mon, 3 Oct 2016 13:48:47 +0300 Subject: [PATCH] Try to commit all --- AculData/AculCalibration.cpp | 166 +++++++++++++++++++++++++---------- AculData/AculCalibration.h | 14 +-- TELoss/TELoss.h | 64 +++++++++++--- calibration1.cxx | 31 +++---- calibration2.cxx | 41 ++++++--- parforcal.par | 8 +- 6 files changed, 229 insertions(+), 95 deletions(-) diff --git a/AculData/AculCalibration.cpp b/AculData/AculCalibration.cpp index 98324f5..23bff31 100755 --- a/AculData/AculCalibration.cpp +++ b/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 diff --git a/AculData/AculCalibration.h b/AculData/AculCalibration.h index 00772d5..3631444 100755 --- a/AculData/AculCalibration.h +++ b/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"); diff --git a/TELoss/TELoss.h b/TELoss/TELoss.h index 80a42df..f0c5fd8 100755 --- a/TELoss/TELoss.h +++ b/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(); diff --git a/calibration1.cxx b/calibration1.cxx index 0e25c9f..9d291c6 100644 --- a/calibration1.cxx +++ b/calibration1.cxx @@ -1,21 +1,18 @@ { - 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); + } } diff --git a/calibration2.cxx b/calibration2.cxx index 4e7134d..d3c9cd0 100644 --- a/calibration2.cxx +++ b/calibration2.cxx @@ -1,33 +1,52 @@ { - 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(); + */ + } diff --git a/parforcal.par b/parforcal.par index 0d85974..68cb20b 100755 --- a/parforcal.par +++ b/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 // -- 2.18.1