////////////////////////////////////////////////////////// // // // AculCalibration // // // // //Some description of this very useful class, //its properties and a short example how to //use it. This text may be partly used in //PhD thesis. // //Description of the detector itself. // // // ////////////////////////////////////////////////////////// #include "AculCalibration.h" ClassImp(AculCalibration); AculCalibration::AculCalibration() { //default constructor fCurrentHStack = NULL; fCurrentHistList.IsOwner(); kRaNOPEAKS = 0; fLowerPeakRelativeHight = 0.; fUpperPeakRelativeHight = 0.; fPeakPositionTolerance = 0.; fFitFuncLineWidth = 1; fFitMinSigma = 0.; fFitPeakThreshold = 0.; for(Int_t i = 0; i < DEFAULTNOPEAKS; i++) { fEnergy[i] = 0; } fCalInformation = 0; Reset(); } AculCalibration::AculCalibration(const char* parfile) { //constructor which fills fA, fB, fC, fD from file parfile fCurrentHStack = NULL; fCurrentHistList.IsOwner(); kRaNOPEAKS = 0; fLowerPeakRelativeHight = 0.; fUpperPeakRelativeHight = 0.; fPeakPositionTolerance = 0.; fFitFuncLineWidth = 1; fFitMinSigma = 0.; fFitPeakThreshold = 0.; for(Int_t i = 0; i < DEFAULTNOPEAKS; i++) { fEnergy[i] = 0; } fCalInformation = 0; SetCalibrationParameters(parfile); } AculCalibration::~AculCalibration() { DeleteStacks(); // delete fCalInformation; // fCalInformation->Close(); } Int_t AculCalibration::SearchPeaks(const TH1 *hin, Double_t sigma, Option_t *option, Double_t threshold, const Int_t searchedpeaks) { //Function searching peaks in inputed TH1 spectrum and selects the peaks in the histogram. // // hin: // sigma: // option: // threshold: // searchedpeaks: TSpectrum sc; //by default for 100 peaks Int_t nopeaks = sc.Search(hin, sigma, "goff", threshold); TString opt = option; opt.ToLower(); const Double_t tStep = 0.05; while ( nopeaks > searchedpeaks && threshold <= 1) { threshold = threshold + tStep; nopeaks = sc.Search(hin, sigma, "goff", threshold); } if (!nopeaks) { return 0; } if (opt.Contains("goff")) { return nopeaks; } TPolyMarker *pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker"); if (pm) { hin->GetListOfFunctions()->Remove(pm); delete pm; } pm = new TPolyMarker(nopeaks, sc.GetPositionX(), sc.GetPositionY()); hin->GetListOfFunctions()->Add(pm); pm->SetMarkerStyle(23); pm->SetMarkerColor(kRed); pm->SetMarkerSize(1.3); return nopeaks; } Int_t AculCalibration::PeaksFitting(TH1I* hSpectrum, Option_t* option, Double_t sigmamin) { //function searching peaks in inputed spectrum // // hSpectrum: // option: posible option are // GP: explanation needed // WRITEBAD: explanation needed // Q: displays on the histogram how the spectrum were fitted and doesn't writes out // the numbers of channels in which are peaks // V: displays on the histogram how the spectrum were fitted and writes out the numbers // of channels in which are peaks // sigmamin: if (!hSpectrum) return 1; Int_t dimension = hSpectrum->GetDimension(); if (dimension > 1) { Error("PeaksFitting", "Only implemented for 1-d histograms"); return 1; } TString opt = option; opt.ToLower(); if (!kRaNOPEAKS) { Error("PeaksFitting", "kRaNOPEAKS is set to zero; calibration spectrum must be set"); return 1; } Int_t peaksNumber = SearchPeaks(hSpectrum, sigmamin, "", fFitPeakThreshold, kRaNOPEAKS); if (peaksNumber != kRaNOPEAKS) { Info("PeaksFitting", "In histogram %s was found %d peaks", hSpectrum->GetName(), peaksNumber); return 1; } //predelat jako volitelny vypis Info("PeaksFitting", "Number of peaks: %d", peaksNumber); Double_t peak[peaksNumber]; //pracovni pole pro zapis piku v neusporadanem poradi Double_t *peakPosition; Double_t *peakHight; for (Int_t i = 0; i < peaksNumber; i++) { TList *functions = hSpectrum->GetListOfFunctions(); TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker"); peakPosition = pm->GetX(); peakHight = pm->GetY(); Double_t currentHight = peakHight[i]; //posouva se smerem nahoru a urcuje hrubou hranici piku, ktera je urcena 1/10 jeho vysky Int_t j = 0; Int_t fitMin = 0; while ( currentHight > (peakHight[i]*fUpperPeakRelativeHight) ) { j++; fitMin = static_cast(peakPosition[i]) + j; currentHight = hSpectrum->GetBinContent(fitMin); } //totez, ale urcujeme dolni hranici piku currentHight = peakHight[i]; j = 0; Int_t fitMax = 0; while ( currentHight > (peakHight[i]*fLowerPeakRelativeHight) ) { j++; fitMax = static_cast(peakPosition[i]) - j; currentHight = hSpectrum->GetBinContent(fitMax); } //fitting if (opt.Contains("gp")) { Info("PeaksFitting", "Option containing gp"); char fncname[20]; sprintf(fncname, "gaus_aux_%d", i); TF1 *gausAux = new TF1(fncname, "gaus", fitMin - 10, fitMax + 10); //pomocny gaus hSpectrum->Fit(fncname, "0 Q", "", fitMin - 15, fitMax + 15); //prvotni fitovani sprintf(fncname, "auto_gp_%d", i); TF1 *fitAuto = new TF1(fncname, "gaus(0) + pol0(3)", fitMin - 15, fitMax + 15); //fce pro automaticke fitovani fitAuto->SetParameter(0, gausAux->GetParameter(0)); //nastavovani parametru fitovaci fce fitAuto->SetParameter(1, gausAux->GetParameter(1)); fitAuto->SetParameter(2, gausAux->GetParameter(2)); hSpectrum->Fit(fncname, "0 R Q +", "", fitMin - 15, fitMax + 15); //dodelat zapis vsech fci hSpectrum->GetFunction(fncname)->ResetBit(TF1::kNotDraw); peak[i] = fitAuto->GetParameter(1); //zapis asi pozice v kanalech do pomocneho pole if (opt.Contains("V")) { Info("PeaksFitting", "Peak position is\t %4.2f \tresolution is \t %2.1f %%", fitAuto->GetParameter(1), 235*(fitAuto->GetParameter(2))/(fitAuto->GetParameter(1))); } } else { char fncname[20]; sprintf(fncname, "auto_g%d", i); TF1 *fitAuto = new TF1(fncname, "gaus", fitMin, fitMax); //fce pro automaticke fitovani fitAuto->SetLineWidth(fFitFuncLineWidth); hSpectrum->Fit(fncname, "+ 0 R Q", ""/*, fitMin - 1, fitMax + 1*/); // hSpectrum->GetFunction(fncname)->ResetBit(TF1::kNotDraw); hSpectrum->GetFunction(fncname)->InvertBit(TF1::kNotDraw); peak[i] = fitAuto->GetParameter(1); //zapis asi pozice v kanalech do pomocneho pole if (opt.Contains("v")) { Info("PeaksFitting", "Peak position is\t%4.2f\tresolution is \t%2.1f %%", fitAuto->GetParameter(1), 235*(fitAuto->GetParameter(2))/(fitAuto->GetParameter(1))); } }//else //end of fitting }//for over all analyzed peaks //peaks sorting Int_t j[peaksNumber]; TMath::Sort(peaksNumber, peak, j, kFALSE); for (Int_t i = 0; i < 4; i++) { fPeak[i] = peak[j[i]]; } if (!opt.Contains("q") || opt.Contains("v")) { Info("PeaksFitting", "Control output:"); for (Int_t i = 0; i < peaksNumber; i++) { printf("\tPeak position is\t%f\n", fPeak[i]); } } // provest kontrolu pomerne polohy piku, // 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")) { fCalInformation->cd(); hSpectrum->Write(); } return 2; } }//for return 0; } Bool_t AculCalibration::SetInputParameters(const char* inputparfile) { // Function which loads the data file for calibration // // -inputparfile: file containing information about calibration source // -In file with the data must be preserved systematic // -There can't be pause between the lines // //For example //................................................................ //223Ra // //4 nopeaks //number of peaks //4.415 E1 //in MeV //5.153 E2 //in MeV //5.683 E3 //in MeV //7.419 E4 //in MeV //100 lowerchannel //in channels //4096 upperchannel //in channels //0.1 lowerpeakhight //in relative units //0.1 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 // //................................................................ const Int_t lineLength = 400; Char_t line[lineLength]; Char_t parameter[100]; Char_t identificator[100]; ifstream fipr; fipr.open(inputparfile); if (!fipr.is_open()) { Error("SetInputsParameters", "File with input parameters was not opened"); return kFALSE; } while (!fipr.eof()) { fipr.getline(line, lineLength); if (strlen(line) == 0) { continue; } sscanf(line, "%s %s", parameter, identificator); if ( strcmp(identificator, "nopeaks") == 0 ) { kRaNOPEAKS = static_cast(atoi(parameter)); for (Int_t i = 0; i < kRaNOPEAKS; i++) { fipr.getline(line, lineLength); sscanf(line, "%s", parameter); fEnergy[i] = static_cast(atof(parameter)); } }//if if ( strcmp(identificator, "lowerpeakhight") == 0 ) { sscanf(line, "%s", parameter); fLowerPeakRelativeHight = static_cast(atof(parameter)); } if ( strcmp(identificator, "upperpeakhight") == 0 ) { sscanf(line, "%s", parameter); fUpperPeakRelativeHight = static_cast(atof(parameter)); } if ( strcmp(identificator, "peakpositiontolerance") == 0 ) { sscanf(line, "%s", parameter); fPeakPositionTolerance = static_cast(atof(parameter)); } if ( strcmp(identificator, "fitfunctionlinewidth") == 0 ) { sscanf(line, "%s", parameter); fFitFuncLineWidth = static_cast(atoi(parameter)); } if ( strcmp(identificator, "minfitsigma") == 0 ) { sscanf(line, "%s", parameter); fFitMinSigma = static_cast(atof(parameter)); } if ( strcmp(identificator, "fithightthreshold") == 0 ) { sscanf(line, "%s", parameter); fFitPeakThreshold = static_cast(atof(parameter)); } } fipr.close(); return kTRUE; } void AculCalibration::PrintInputParameters() { //print alpha source parameters cout << "Number of peaks: " << kRaNOPEAKS << endl << endl; for (Int_t i = 0; i < kRaNOPEAKS; i++) { cout << "fEnergy[" << i << "] = " << fEnergy[i] << endl; } return; } Bool_t AculCalibration::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 const Int_t lineLength = 200; char line[lineLength]; // Int_t crate; char crate[100]; Int_t i, j; char cA[40], cB[40], cC[40], cD[40]; //open file with calibration parameters ifstream calFileR; calFileR.open(calparfile); if( !calFileR.is_open() ) { Error("SetCalibrationParameters", "File %s with calibration data was not opened", calparfile); return kFALSE; } Reset(); //read calibration parameters from file while (!calFileR.eof()) { calFileR.getline(line, lineLength); if ( line[0] != '*' && line[0] != '#' && line[0] != '%') { //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); fC[i][j] = atof(cC); fD[i][j] = atof(cD); } } calFileR.close(); return kTRUE; } void AculCalibration::PrintCalibrationParameters(const Int_t blockmin, const Int_t blockmax) { //Print calibration parameters fA, fB, fC, fD // // blockmin: minimum block data are displayed for // blockmax: maximum block data are displayed for 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; } } } void AculCalibration::ShowRawSpectra(const char* filename, const Int_t block, TCanvas* rawCanvas, Int_t xaxismin, Int_t xaxismax, /*TObjArray* histList,*/ const Int_t subaddress) { //Displays the spectrum from a file, divides the canvas into a sufficient number of pads and displays spectrums of each //block subaddress on the suitable pads. // // filename: input .root file containing spectra to be showed // block: block which will be drawn // rawCanvas: canvas on which one you will see the spectrum // xaxismin: minimal channel, which will be displayed // xaxismax: maximal channel, which will be displayed // subaddress: Char_t address[40]; Char_t histName[40]; Char_t fillCommand[40]; Char_t fillCondition[40]; if (!rawCanvas) { //rawCanvas = new TCanvas("RawSpectra", "Raw spectra in channels", 1); cout << "You have to assign TCanvas for raw spectra drawing" << endl; return; } rawCanvas->Clear(); // cout << "hovno" << endl; rawCanvas->SetFillColor(10); // cout << "hovno" << endl; TFile *fr = new TFile(filename); if (fr->IsOpen() == 0) { cout << endl << "File " << filename << " was not opened and won't be processed" << endl << endl; return; } TH1I *hRead = 0; TTree *tr = (TTree*)fr->Get("RAW"); // cout << "hovno" << endl; if (subaddress > 15) { rawCanvas->Divide(4, 4); rawCanvas->SetFillColor(10); // cout << "hovno" << endl; for (Int_t i = 0; i < 16; i++) { cout << i << endl; rawCanvas->cd(i+1); hRead = new TH1I("name", "title", 4096, 0, 4095); sprintf(address, "C3[%d][%d]", block, i); sprintf(histName, "H3[%d][%d]", block, i); hRead->SetName(histName); sprintf(fillCommand, "%s >> %s", address, hRead->GetName()); sprintf(fillCondition, "%s > 0", address); tr->Draw(fillCommand, fillCondition, ""); if (hRead) { hRead->SetDirectory(0); // cout << hRead->GetEntries() << endl; // if (fHRawList) { // fHRawList->Add(hRead); fHRawList.Add(hRead); // } hRead->SetAxisRange(xaxismin, xaxismax); } }//for } else { fr->cd(); hRead = new TH1I("name", "title", 4096, 0, 4095); sprintf(address, "C3[%d][%d]", block, subaddress); sprintf(histName, "H3[%d][%d]", block, subaddress); hRead->SetName(histName); sprintf(fillCommand, "%s >> %s", address, hRead->GetName()); sprintf(fillCondition, "%s > 0", address); // cout << fillCommand << setw(20) << fillCondition << endl; tr->Draw(fillCommand, fillCondition, "goff"); if (hRead) { hRead->SetDirectory(0); // if (fHRawList) { // fHRawList->Add(hRead); // } fHRawList.Add(hRead); hRead->Draw(); hRead->SetAxisRange(xaxismin, xaxismax); } }//else fr->Close(); rawCanvas->Update(); return; } void AculCalibration::ShowSpectra(const char* filename, TCanvas* rawCanvas, Option_t *option, Int_t xaxismin, Int_t xaxismax, const Int_t subaddress) { //filename: input .root file with saved filled histograms to be showed //rawCanvas: canvas on which one you will see the spectrum //option: THStack options //xaxismin: Minimum channel, which will be displayed //xaxismax: Maximum channel, which will be displayed //subaddress: TString opt = option; opt.ToLower(); if (!rawCanvas) { Error("ShowRawSpectra", "You have to assign TCanvas for raw spectra drawing"); return; } rawCanvas->Clear(); TFile fr(filename); if (fr.IsOpen() == 0) { Error("ShowRawSpectra", "File %s was not opened and won't be processed", filename); return; } TList *histList; histList = fr.GetListOfKeys(); Int_t listEntries = histList->GetEntries(); TH1 *hDraw = 0; DeleteStacks(); if (subaddress >= listEntries) { fCurrentHStack = new THStack(); for (Int_t i = 0; i < listEntries; i++) { //zkontrolovat hranice Info("ShowRawSpectra", "Histogram with spectrum of subaddress %d is loading", i); fr.GetObject(histList->At(i)->GetName(), hDraw); if (hDraw) { hDraw->SetDirectory(0); //hDraw->SetAxisRange(xaxismax, xaxismin); //nefunguje fCurrentHistList.Add(hDraw); fCurrentHStack->Add(hDraw); } }//for if ( !fCurrentHStack->GetHists()->IsEmpty() ) { Info("ShowRawSpectra", "Histogram stack drawing"); fCurrentHStack->Draw(opt.Data()); } }//if all subaddresses else { //zkontrolovat fr.GetObject(histList->At(subaddress)->GetName(), hDraw); if (hDraw) { hDraw->SetAxisRange(xaxismin, xaxismax, "X"); hDraw->Draw(); hDraw->SetDirectory(0); fCurrentHistList.Add(hDraw); } }//else fr.Close(); rawCanvas->Update(); return; } void AculCalibration::FillRawSpectraFile(const char* rawdatafile, const char* block, const char* treename, TCanvas* rawCanvas, Option_t *option, Int_t xaxismin, Int_t xaxismax) { //filename: input .root file containing spectra to be showed //block: //rawCanvas: //xaxismin: //xaxismax: //variables to be became function parameter TString opt(option); opt.ToLower(); if (!rawCanvas) { Error("ShowRawSpectra", "You have to assign TCanvas for raw spectra drawing"); return; } rawCanvas->Clear(); TFile fr(rawdatafile); if (fr.IsOpen() == 0) { Error("ShowRawSpectra", "File %s was not opened and won't be processed", rawdatafile); return; } TTree *tr = (TTree*)fr.Get(treename); char outputfile[300]; sprintf(outputfile, "%s[]Raw.root", block); TFile fw(outputfile, opt.Data()); if (fw.IsOpen() == 0) { Error("CalculateCalibParameters", "Output file %s was not created.", outputfile); return; } if (fw.IsWritable() == 0) { Error("CalculateCalibParameters", "Output file %s is not writable. Set option to \"RECREATE\".", outputfile); return; } char address[40]; char histName[40]; char histTitle[40]; char fillCommand[40]; char fillCondition[40]; fw.cd(); TH1I *hRead = 0; for (Int_t i = 0; i < 16; i++) { //zkontrolovat hranice cout << i << endl; //predelat na info hRead = new TH1I("name", "title", 4096, 0, 4095); sprintf(address, "%s[%d]", block, i); sprintf(histName, "%s[%d]", block, i); sprintf(histTitle, "%s : %s", rawdatafile, histName); hRead->SetName(histName); hRead->SetTitle(histTitle); sprintf(fillCommand, "%s >> %s", address, hRead->GetName()); sprintf(fillCondition, "%s > 0", address); tr->Draw(fillCommand, fillCondition, "goff"); //prozkoumat goff hRead->Write(); }//for fw.Close(); fr.cd(); delete tr; fr.Close(); rawCanvas->Update(); return; } Bool_t AculCalibration::CalculateCalibParameters(const char* inputfilename, const Int_t block, Int_t lowerchannel, Int_t upperchannel, Int_t lowersubaddress, Int_t uppersubaddress) { //function for calculation of calibrate parameters for DAQ system based on "Black Windows" program // // inputfile: root file with calibration spectra // block: block to be calibrated as number // lowerchannel: minimal channel from which the spectrum will be analysed // upperchannel: maximal channel up to which the spectrum will be analysed // lowersubaddress: block subaddress // uppersubaddress: block subbaddress if ( (block > BLOCKSNUMBER) || (block < 1) ) { Error("CalculateCalibParameters", "Possible block values have to be in range 1 - %d", BLOCKSNUMBER - 1); return kFALSE; } if ( (uppersubaddress - lowersubaddress) >= ADDRESSNUMBER ) { Error("CalculateCalibParameters", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1); return kFALSE; } //promenne potrebne pro fitovani TF1 *calFunction = new TF1("calib", "pol1", 0, 1000); //predelat jako lokalni promennou fce TGraph *calGraph = new TGraph(kRaNOPEAKS, fPeak, fEnergy); //lokalni promenna, dohodit pocet vstupu pomoci parametru //auxiliary variables, particularly for text parameter fields Char_t outputfilename[40]; Char_t address[40]; Char_t histName[40]; Char_t fillCommand[40]; Char_t fillCondition[40]; Int_t fitControl = 0; Char_t histTitle[40]; //creation of the output text file if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "C3[%d][].cal", block); } else { if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "C3[%d][%d].cal", block, lowersubaddress); } else { sprintf(outputfilename, "C3[%d][%d-%d].cal", block, lowersubaddress, uppersubaddress); } } // fFileName = outputfilename; //doubtful // cout << "hovno1" << endl; // cout << outputfilename << endl; // fOutCalFile.open(outputfilename); ofstream outcalfile; outcalfile.open(outputfilename); // cout << "hovno2" << endl; // if (!fOutCalFile.is_open()) { if (!outcalfile.is_open()) { Error("CalculateCalibParameters", "Output file %s was not opened", outputfilename); return kFALSE; } //creation of the output root file if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "C3[%d][].root", block); } else { if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "C3[%d][%d].root", block, lowersubaddress); } else { sprintf(outputfilename, "C3[%d][%d-%d].root", block, lowersubaddress, uppersubaddress); } } fCalInformation = new TFile(outputfilename, "RECREATE"); if (fCalInformation->IsOpen() == 0) { Error("CalculateCalibParameters", "File %s was not opened and won't be processed", outputfilename); return kFALSE; } //input file with raw data opening TFile *fr = new TFile(inputfilename); if (fr->IsOpen() == 0) { Error("CalculateCalibParameters", "File %s was not opened and won't be processed", inputfilename); return kFALSE; } TTree *tr = (TTree*)fr->Get("RAW"); if (!tr) { Error("CalculateCalibParameters", "Tree \"RAW\" was not found in file %s", inputfilename); return kFALSE; } //raw data histogram filling TH1I *hRaw = 0; for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) { Info("\n\nCalculateCalibParameters", "Calculating calibration parametres for address C3[%d][%d].", block, i); // cout << endl // << endl // << "Calculating calibration parametres for address C3[" << block << "][" << i << "]." << endl; //TH1I object preparing hRaw = new TH1I("name", "title", 4096, 0, 4095); sprintf(address, "C3[%d][%d]", block, i); sprintf(histName, "H3[%d][%d]", block, i); hRaw->SetName(histName); sprintf(fillCommand, "%s >> %s", address, hRaw->GetName()); sprintf(fillCondition, "%s > 0", address); //filling from the .root raw data file and content arrangement tr->Draw(fillCommand, fillCondition, "goff"); if (lowerchannel != 0) { for (Int_t i = 0; i < lowerchannel; i++) { hRaw->SetBinContent(i, 0); } } if (upperchannel != 0) { for (Int_t i = upperchannel; i < 4095; i++) { hRaw->SetBinContent(i, 0); } } //spectrum analysis fitControl = PeaksFitting(hRaw, "", fFitMinSigma); // cout << "Value of fitControl is: " << fitControl << endl; Info("CalculateCalibParameters", "Value of fitControl is: %d", fitControl); //incorrectly treated spectrum output if (fitControl != 0) { // fOutCalFile << setw(39) << fitControl << endl; outcalfile << setw(39) << fitControl << endl; fCalInformation->cd(); hRaw->SetLineColor(2); //red hRaw->Write(); continue; } //correctly treated spectrum saving if (fCalInformation->IsOpen()) { fCalInformation->cd(); hRaw->SetLineColor(3); //green hRaw->Write(); } //calibration parameters calculation for (Int_t j = 0; j < kRaNOPEAKS; j++) { //delat podle poctu zkoumanych piku calGraph->SetPoint(j, fPeak[j], fEnergy[j]); //calibration graph filling // cout << "vypis:\t fPeak = " << fPeak[i] << "\tfEnergy = " << fEnergy[i] << endl; }//for calGraph->Fit(calFunction, "", "goff", 0, 4096); //omezit hlasitost fitovani, udelat volitelne // fOutCalFile << std::right outcalfile << std::right << setw(2) << "3" << setw(4) << block << setw(4) << i << setw(12) << setprecision(4) << calFunction->GetParameter(1) << setw(12) << setprecision(4) << calFunction->GetParameter(0) << setw(5) << fitControl // << setw(10) << calFunction->GetParameter(1)*2.35*(fSigma[0] + fSigma[1] + fSigma[2] + fSigma[3])/4. << endl; fA[block][i] = calFunction->GetParameter(1); fB[block][i] = calFunction->GetParameter(0); }//for //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! //urcite by se to dalo hodit do samostatne fce //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! //outputfile with calibrated spectra if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "C3[%d][]E.root", block); } else { if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "C3[%d][%d]E.root", block, lowersubaddress); } else { sprintf(outputfilename, "C3[%d][%d-%d]E.root", block, lowersubaddress, uppersubaddress); } } TFile *fw = new TFile(outputfilename, "RECREATE"); if (fw->IsOpen() == 0) { Error("\nCalculateCalibParameters", "File %s was not created and won't be processed\n\n", outputfilename); // cout << endl << "File " << outputfilename << " was not created and won't be processed" << endl << endl; return kFALSE; } AculRaw *eventr = new AculRaw(); Long64_t noEntries = tr->GetEntries(/*getCondition*/); tr->SetBranchAddress("channels", &eventr); TH1F *hEnergy = 0; TRandom3 ranGen(1); for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) { cout << "Calibration spectrum from address C3[" << block << "][" << i << "]." << endl; sprintf(histName, "C3[%d][%d]", block, i); sprintf(histTitle, "%s: %s", inputfilename, histName); hEnergy = new TH1F(histName, histTitle, 10000, 0., 10.); for (Int_t j = 0; j < noEntries; j++) { tr->GetEntry(j); if ((eventr->C3[block][i]) > 0) { hEnergy->Fill( (eventr->C3[block][i] + ranGen.Uniform(-0.5, 0.5) )*fA[block][i] + fB[block][i] ); } } fw->cd(); hEnergy->Write(); }//for fw->Close(); ///////////////////////////////////////////////////////// fr->Close(); fCalInformation->Close(); delete fCalInformation; //pokusne // fOutCalFile.close(); outcalfile.close(); return kTRUE; } Bool_t AculCalibration::CalculateCalibParameters(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) { //function is not completely ready to use // // //function for calculation of calibrate parameters for DAQ system based on "Go4" // // inputfile: root file with calibration spectra // block: block name to be calibrated // lowerchannel: minimal channel from which the spectrum will be analysed // upperchannel: maximal channel up to which the spectrum will be analysed // lowersubaddress: block subaddress // uppersubaddress: block subbaddress if (kRaNOPEAKS == 0) { Error("CalculateCalibParameters", "Alpha source parameters was not red"); return 0; } if ( (address > BLOCKSNUMBER) || (address < 1) ) { Error("CalculateCalibParameters", "Possible address values have to be in range 1 - %d", BLOCKSNUMBER - 1); return kFALSE; } //muzu nechat if ( (uppersubaddress - lowersubaddress) >= ADDRESSNUMBER ) { Error("CalculateCalibParameters", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1); return 0; } //probrat podle potreby //auxiliary variables, particularly for text parameter fields char outputfilename[300]; //creation of the output text file if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "%s[].cal", block); } else { if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "%s[%d].cal", block, lowersubaddress); } else { sprintf(outputfilename, "%s[%d-%d].cal", block, lowersubaddress, uppersubaddress); } }//if ofstream outcalfile; outcalfile.open(outputfilename); if (!outcalfile.is_open()) { Error("CalculateCalibParameters", "Output file %s was not opened", outputfilename); return 0; }//if //predelat podle nazvu bloku //creation of the output root file if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "%s[].root", block); } else { if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "%s[%d].root", block, lowersubaddress); } else { sprintf(outputfilename, "%s[%d-%d].root", block, lowersubaddress, uppersubaddress); } } fCalInformation = new TFile(outputfilename, "RECREATE"); if ( !fCalInformation->IsOpen() ) { Error("CalculateCalibParameters", "File %s was not opened and won't be processed", outputfilename); return 0; } //nechat //input file with raw data opening TFile *fr = new TFile(inputfile); if ( !fr->IsOpen() ) { Error("CalculateCalibParameters", "File %s was not opened and won't be processed", inputfile); return 0; } TTree *tr = (TTree*)fr->Get(treename); if (!tr) { Error("CalculateCalibParameters", "Tree %s was not found in file %s", treename, inputfile); return 0; } //promenne potrebne pro fitovani: presunout nize //pohlidat delete TF1 *calFunction = new TF1("calib", "pol1", 0, 1000); //predelat jako lokalni promennou fce TGraph *calGraph = new TGraph(kRaNOPEAKS, fPeak, fEnergy); //lokalni promenna, dohodit pocet vstupu pomoci parametru Char_t detectorChannel[100]; Char_t histName[100]; Char_t fillCommand[512]; Char_t fillCondition[200]; Int_t fitControl = 0; //predelat nazvy histogramu //zrusit cyklus, napsat jako fci //raw data histogram filling TH1I *hRaw = 0; for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) { printf("\n\n"); Info("CalculateCalibParameters", "Calculating calibration parametres for detector channel %s[%d].", block, i); //TH1I object preparing hRaw = new TH1I("name", "title", 4096, 0, 4095); //nastavovat hranice histogramu podle parametru fce sprintf(detectorChannel, "%s[%d]", block, i); sprintf(histName, "Hist%s[%d]", block, i); hRaw->SetName(histName); sprintf(fillCommand, "%s >> %s", detectorChannel, hRaw->GetName()); sprintf(fillCondition, "%s > 0", detectorChannel); //filling from the .root raw data file and content arrangement tr->Draw(fillCommand, fillCondition, "goff"); if (lowerchannel != 0) { //zbytecna cast // for (Int_t i = 0; i < lowerchannel; i++) { // hRaw->SetBinContent(i, 0); // } //muzu napsat jednoduseji: } // vytvaret histogram s pozadovanym rozsahem if (upperchannel != 0) { //zbytecna cast // for (Int_t i = upperchannel; i < 4095; i++) { // taky ovsem nemusim, chtelo by to jeste zvazit hRaw->SetBinContent(i, 0); // } // } // //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 hRaw->Write(); continue; }//if //correctly treated spectrum saving if (fCalInformation->IsOpen()) { //ok fCalInformation->cd(); hRaw->SetLineColor(3); //green 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); }//for //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! //zapis histogramu v MeV do souboru .root // //urcite by se to dalo hodit do samostatne fce //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Char_t histTitle[40]; //outputfile with calibrated spectra if ( (lowersubaddress == 0) && (uppersubaddress == 15) ) { sprintf(outputfilename, "%s[]E.root", block); } else { if (lowersubaddress == uppersubaddress) { sprintf(outputfilename, "%s[%d]E.root", block, lowersubaddress); } else { sprintf(outputfilename, "%s[%d-%d]E.root", block, lowersubaddress, uppersubaddress); } } TFile *fw = new TFile(outputfilename, "RECREATE"); if (fw->IsOpen() == 0) { Error("CalculateCalibParameters", "File %s was not created and won't be processed\n\n", outputfilename); return 1; } TH1F *hEnergy = 0; TRandom3 ranGen(1); for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) { Info("CalculateCalibParameters", "Calibration spectrum from address %s[%d]", block, i); sprintf(histName, "Hist%s[%d]E", block, i); sprintf(histTitle, "%s: %s", inputfile, histName); hEnergy = new TH1F(histName, histTitle, 10000, 0., 10.); sprintf(detectorChannel, "%s[%d]", block, i); hEnergy->SetName(histName); sprintf(fillCommand, "%f*(%s+%f)+%f >> %s", fA[address][i], detectorChannel, ranGen.Uniform(-0.5, 0.5), fB[address][i], hEnergy->GetName()); sprintf(fillCondition, "%s > 0", detectorChannel); //filling from the .root raw data file and content arrangement tr->Draw(fillCommand, fillCondition, "goff"); fw->cd(); hEnergy->Write(); }//for fw->Close(); fr->Close(); fCalInformation->Close(); delete fCalInformation; //pokusne outcalfile.close(); return 1; } Int_t AculCalibration::CalibrateBlock(const Char_t* inputfilename, const Int_t block, const Char_t* outputfilename, Int_t lowersubaddress, Int_t uppersubaddress) { //probably some of the obsolete functions, maybe does not do anything important //input data root file opening TFile *fr = new TFile(inputfilename); if (fr->IsOpen() == 0) { cout << endl << "File " << inputfilename << " was not opened and won't be processed" << endl << endl; return 1; } TTree *tr = (TTree*)fr->Get("RAW"); if (!tr) { cout << "Tree \"RAW\" was not found in file " << inputfilename << "." << endl; return 1; } //outputfile with calibrated spectra TFile *fw = new TFile(outputfilename, "RECREATE"); if (fw->IsOpen() == 0) { cout << endl << "File " << outputfilename << " was not created." << endl << endl; return 1; } //check for the calibration parameters for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) { if ( fA[block][i] == 0 && fB[block][i] == 0) { Warning("CalibrateBlock", "Calibration parameters are empty"); return 2; } } //auxiliary variables for histogram names Char_t histName[40]; Char_t histTitle[40]; //beginning of the tree reading AculRaw *eventr = new AculRaw(); Long64_t noEntries = tr->GetEntries(/*getCondition*/); tr->SetBranchAddress("channels", &eventr); TH1F *hEnergy = 0; TRandom3 ranGen(1); for (Int_t i = lowersubaddress; i <= uppersubaddress; i++) { cout << "Calibration spectrum from address C3[" << block << "][" << i << "]." << endl; sprintf(histName, "C3[%d][%d]", block, i); sprintf(histTitle, "%s: %s", inputfilename, histName); hEnergy = new TH1F(histName, histTitle, 10000, 0., 10.); for (Int_t j = 0; j < noEntries; j++) { tr->GetEntry(j); if ((eventr->C3[block][i]) > 0) { //v tomto miste se potom muze predelat jako bud vyplnovani histogramu nebo noveho stromu hEnergy->Fill( (eventr->C3[block][i] + ranGen.Uniform(-0.5, 0.5) )*fA[block][i] + fB[block][i] ); } } fw->cd(); hEnergy->Write(); }//for fw->Close(); return 0; } void AculCalibration::ShowAnalyzedSpectra(const char *filename, TCanvas* fittedRawCanvas, Int_t xaxismin, Int_t xaxismax, Int_t subaddress) { //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 . //Selects the peaks in the histogram and displays on the histogram how the spectrum were fitted. // // filename: file .root containing analysed spectra // fittedRawCanvas: canvas on which one you will see the spectrum // xaxismin: Minimum channel, which will be displayed // xaxismax: Maximum channel, which will be displayed // subaddress: if ( subaddress > ADDRESSNUMBER ) { Error("ShowAnalyzedSpectra", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1); return; } if (!fittedRawCanvas) { Warning("ShowAnalyzedSpectra", "You have to assign TCanvas for fitted raw spectra drawing"); return; } TFile *fr = new TFile(filename, "READ"); if (!fr->IsOpen()) { cout << "File " << filename << " was not opened" << endl; return; } TList *histList; histList = fr->GetListOfKeys(); Int_t listEntries = histList->GetEntries(); TH1I *hDraw = 0; fittedRawCanvas->Clear(); // fittedRawCanvas->SetFillColor(10); if ( (listEntries > 1) && (listEntries <= 8) ) { fittedRawCanvas->Divide(2, 4); fittedRawCanvas->SetFillColor(10); } if ( (listEntries > 8) && (listEntries <= 16) ) { fittedRawCanvas->Divide(4, 4); fittedRawCanvas->SetFillColor(10); } if (subaddress >= listEntries) { for (Int_t i = 0; i < listEntries; i++) { fittedRawCanvas->cd(i+1); fr->GetObject(histList->At(i)->GetName(), hDraw); if (hDraw) { hDraw->SetAxisRange(xaxismin, xaxismax, "X"); hDraw->Draw(); hDraw->SetDirectory(0); // if (fHAnalyzedList) { // fHAnalyzedList->Add(hDraw); // } fHAnalyzedList.Add(hDraw); } }//for } else { fr->GetObject(histList->At(subaddress)->GetName(), hDraw); if (hDraw) { hDraw->SetAxisRange(xaxismin, xaxismax, "X"); hDraw->Draw(); hDraw->SetDirectory(0); fHAnalyzedList.Add(hDraw); } } fr->Close(); fittedRawCanvas->Update(); return; } void AculCalibration::ShowEnergySpectra(const char *filename, TCanvas* energyCanvas, const Int_t subaddress, Option_t* option, Double_t xaxismin, Double_t xaxismax) { //Displays the spectrum of the selected subbaddress block in MeV // // filename: file .root containing calibrated spectra in MeV // energyCanvas: : canvas on which one you will see the spectrum // subaddress: block subaddress which will be drawn // option: sum ,+ ,c // xaxismin: Minimum channel, which will be displayed // xaxismax: Maximum channel, which will be displayed if ( subaddress > ADDRESSNUMBER ) { Error("ShowEnergySpectra", "Possible subaddress values have to be in range 0 - %d", ADDRESSNUMBER - 1); return; } if (!energyCanvas) { Warning("ShowEnergySpectra", "You have to assign TCanvas for fitted raw spectra drawing"); return; } TString opt = option; opt.ToLower(); TFile *fr = new TFile(filename, "READ"); if (!fr->IsOpen()) { cout << "File " << filename << " was not opened" << endl; return; } TList *histList; histList = fr->GetListOfKeys(); Int_t listEntries = histList->GetEntries(); TH1F *hDraw = 0; energyCanvas->Clear(); energyCanvas->SetFillColor(10); if ( (listEntries > 1) && (listEntries <= 8) ) { energyCanvas->Divide(2, 4); energyCanvas->SetFillColor(10); } if ( (listEntries > 8) && (listEntries <= 16) ) { energyCanvas->Divide(4, 4); energyCanvas->SetFillColor(10); } if (subaddress >= listEntries) { if (opt.Contains("sum")) { energyCanvas->cd(0); for (Int_t i = 0; i < listEntries; i++) { fr->GetObject(histList->At(i)->GetName(), hDraw); if (hDraw) { hDraw->SetDirectory(0); //hDraw->SetAxisRange(xaxismin, xaxismax, "X"); if (opt.Contains("c")) { hDraw->SetLineColor(i+1); } if (opt.Contains("c") && opt.Contains("+")) { hDraw->SetFillColor(i+1); } // fHEnergyStack->Add(hDraw); fHEnergyStack.Add(hDraw); } } // if (fHEnergyStack) { // if (opt.Contains("+")) { fHEnergyStack->Draw(); } // else { fHEnergyStack->Draw("nostack"); } // } if (opt.Contains("+")) { fHEnergyStack.Draw(); } else { fHEnergyStack.Draw("nostack"); } } else { for (Int_t i = 0; i < listEntries; i++) { energyCanvas->cd(i+1); fr->GetObject(histList->At(i)->GetName(), hDraw); if (hDraw) { hDraw->SetAxisRange(xaxismin, xaxismax, "X"); hDraw->Draw(); hDraw->SetDirectory(0); // if (fHEnergyList) { // fHEnergyList->Add(hDraw); // } fHEnergyList.Add(hDraw); } }//for }//else }//if else { fr->GetObject(histList->At(subaddress)->GetName(), hDraw); energyCanvas->cd(0); if (hDraw) { hDraw->SetAxisRange(xaxismin, xaxismax, "X"); hDraw->Draw(); hDraw->SetDirectory(0); // if (fHEnergyList) { // fHEnergyList->Add(hDraw); // } fHEnergyList.Add(hDraw); } }//else fr->Close(); // fFileName = filename; // fFileName.Resize(fFileName.Length() - 6); // fFileName.Append(".cal", 4); energyCanvas->Update(); return; } Bool_t AculCalibration::AddCalFileToList(const char* calfilelist) { //this function does not work at the moment //some problem with TString object fFileName TString fl = calfilelist; fl.ToLower(); ofstream fw; fw.open(fl.Data(), ofstream::app); if (!fw.is_open()) { cout << "File " << fl.Data() << " was not opened" << endl; return kFALSE; } // fw << fFileName.Data() << endl; fw.close(); return kTRUE; } void AculCalibration::ClearHistograms(Option_t* option) { //clear THStack and TObjArray members //this function will be removed as soon as possible TString opt = option; opt.ToLower(); // fHRawList->Add(hRead); fHRawList.Clear(); fHAnalyzedList.Clear(); fHEnergyList.Clear(); fHEnergyStack.Clear(); return; } void AculCalibration::MakeCalibrationFile(Char_t* calibrationfile, Char_t *calfilelist) { //calibrationfile: file with calibration parameters to be created //calfilelist: file containing list of existing text files with calibration parameters ifstream calListR; calListR.open(calfilelist); if( !calListR.is_open() ) { cout << "File with list of calibration files was not opened" << endl; return; } //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; } } const Int_t lineLength = 100; char line[lineLength]; char filename[50]; // ifstream calFileR; Int_t crate, i, j, id; // int crate, i, j, id; char cA[40], cB[40], cSigma[40]; while (!calListR.eof()) { calListR.getline(line, lineLength); // cout << line << endl; sscanf(line, "%s", filename); cout << filename << endl; ifstream calFileR; calFileR.open(filename); if (calFileR.is_open()) { cout << filename << " processing" << endl; calFileR.seekg(0); while (!calFileR.eof()) { // cout << " in inner while" << endl; calFileR.getline(line, lineLength); 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)); // fMeanSigma[i][j] = static_cast(atof(cSigma)); cout << fA[i][j] << "\t" << fB[i][j] << endl; }//if }//while calFileR.close(); cout << "calFileR was closed" << endl << endl; }//if }//while calListR.close(); ofstream CalibFileW; CalibFileW.open(calibrationfile); if (!CalibFileW.is_open()) { cout << "Calibration file was not opened" << endl; return; } for (Int_t i = 0; i < BLOCKSNUMBER; i++) { for (Int_t j = 0; j < ADDRESSNUMBER; j++) { if (fA[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(10) << fMeanSigma[i][j] << endl; } } } CalibFileW.close(); return; } void AculCalibration::DeleteStacks(Option_t* option) { if (fCurrentHStack) { delete fCurrentHStack; fCurrentHStack = NULL; } fCurrentHistList.Delete(); return; } void AculCalibration::Reset() { //reset calibration parameters fA, fB, 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; fC[i][j] = 0; fD[i][j] = 0; // fMeanSigma[i][j] = 0; } } return; }