Skip to content

A
AculUtils
Source
Select Git revision
...
Target
Select Git revision
 
Commits (4)
Hide whitespace changes
Inline Side-by-side
Showing with 117 additions and 62 deletions
AculData/AculCalibCsI.cpp
View file @ 3a55bcf8
......@@ -125,16 +125,6 @@ void AculCalibCsI::PrintCuts() {
// return;
if (fCuts) printf("Cuts loaded from file \"%s\"\n", fCuts->GetName());
/*for (Int_t i = 0; i < NOCALFILES; i++) {
if (cTOF[i]) {
printf("TOF cut No. %d; Name: \"%s\"\n", i, cutsCol[i]->GetName());
cout << "asdjasd" << endl;
}
else {
printf("TOF cut No. %d was not loaded. Maximal number of cuts is %d\n", i, NOCALFILES);
}
}//for */
TCutG *curCut = 0;
//cutsCol
......@@ -147,19 +137,6 @@ void AculCalibCsI::PrintCuts() {
}*/
}
/* for (Int_t i = 0; i < NOCALFILES; i++) {
if (cQCD[i]) {
printf("QCD cut No. %d; Name: \"%s\"\n", i, cQCD[i]->GetName());
// cout << "kasjhd" << endl;
} else {
// cout << "bhajskd" << endl;
printf("QCD cut No. %d was not loaded. Maximal number of cuts is %d\n", i, NOCALFILES);
// cout << "vsdfjks" << endl;
}
}
printf("AculCalibCsI::PrintCuts: End of function.\n");
*/
return;
}
......@@ -256,11 +233,11 @@ void AculCalibCsI::DrawBeam(TCanvas *canvas, Int_t files, const char* variable)
con.Form("%s[5]>200", variable);
curTree->Draw(var.Data(), con.Data(), "cont");
for (Int_t j = 0; j < energyPoints; j++) {
// if ( cutsCol.At(files + j) ) {
// curCutG = (TCutG*)cutsCol.At(files + j);
// curCutG->Draw("same");
if ( cutsCol.At(files + j) ) {
curCutG = (TCutG*)cutsCol.At(files + j);
curCutG->Draw("same");
// printf("AculCalibCsI::DrawBeam: cQCD cut No. %d cannot be drawn, need to repair this function.\n", j);
// }
}
}
canvas->Update();
}
......@@ -365,14 +342,14 @@ void AculCalibCsI::GetPeakMean(const char* variable, Int_t tree, Int_t energy, T
con.Form("%s[%d]>%d && %s", variable, channel, lowRange, cutsCol.FindObject(beamcut)->GetName());
canvas->cd(i+1);
hname.Form("hfull[%d][%d]", tree, i);
hfull[tree][i] = new TH1I(hname.Data(), "title", nbins, 0, 4096);
hfull[tree][i] = new TH1I(hname.Data(), "title1", nbins, 0, 4096);
curTree->SetLineColor(1);
curTree->Draw(var.Data(), con.Data());
var.Form("%s[%d]>>hcut[%d][%d]", variable, channel, tree, i);
con.Form("%s[%d]>%d && %s[%d]<%d && %s", variable, channel, peakMin[energy][i], variable, channel, peakMax[energy][i], cutsCol.FindObject(beamcut)->GetName());
hname.Form("hcut[%d][%d]", tree, i);
hcut[tree][i] = new TH1I(hname.Data(), "title", nbins, 0, 4096);
hcut[tree][i] = new TH1I(hname.Data(), "title2", nbins, 0, 4096);
hcut[tree][i]->SetLineColor(3);
curTree->Draw(var.Data(), con.Data(), "same");
......@@ -395,7 +372,7 @@ void AculCalibCsI::Calibrate(TCanvas *canvas, Bool_t savefile, const char* filen
canvas->Divide(4,4);
// cout << alphas2.GetSize()+1 << endl;
cout << energyPoints+1 << endl;
cout <<" energyPoints+1 "<<energyPoints+1<< endl;
TString gName;
TString gTitle;
......@@ -432,7 +409,7 @@ void AculCalibCsI::FillGraph(TGraphErrors *g, Int_t npoints, Double_t *energies,
TString opt = option;
//all available energy points and (0,0)
g->SetPoint(0, 0., 0.);
g->SetPoint(0, 0., 0.); //(point number, x, y)
for (Int_t i = 0; i < npoints-1; i++) {
// g->SetPoint(i+1, energies[i], mean[i][graphNumber]);
// g->SetPointError(i+1, 0, meanRMS[i][graphNumber]);
......
AculData/AculCalibration.cpp
View file @ 3a55bcf8
......@@ -26,8 +26,8 @@ AculCalibration::AculCalibration() : fEnergy(0), fEnergyInput(0), fA(0), fB(0),
// todo: change size of fA and fB in some other place
fA.Set(32);
fB.Set(32);
fEnergy.Set(4);
fEnergyInput.Set(4);
// fEnergy.Set(4);
// fEnergyInput.Set(4);
kRaNOPEAKS = 0;
fLowerPeakRelativeHight = 0.;
......@@ -172,7 +172,6 @@ Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t s
Info("PeaksFitting", "In histogram %s was found %d peaks", hSpectrum->GetName(), peaksNumber);
return 1;
}
//should be optional output
Info("PeaksFitting", "Number of peaks in %s: %d", hSpectrum->GetName(), peaksNumber);
......@@ -259,7 +258,8 @@ Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t s
TMath::Sort(peaksNumber, peak, j, kFALSE);
fPeak.Set(peaksNumber);
for (Int_t i = 0; i < peaksNumber; i++) {
fPeak[i] = peak[j[i]];
fPeak[i] = peak[j[i]];
//printf("\tPeak peak\t%f\n", fPeak[i]);
}
if (!opt.Contains("q") || opt.Contains("v")) {
......@@ -273,13 +273,14 @@ Int_t AculCalibration::PeaksFitting(TH1* hSpectrum, Option_t* option, Double_t s
// 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")) {
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])) ) ) {
//printf("\tPeaksFitt fEnergy\t%f\n", fEnergy[i]);
if (fCalInformation /* && opt.Contains("writebad")*/) {
fCalInformation->cd();
hSpectrum->Write();
}
return 2;
//return 2;*/
}
}//for
......@@ -759,7 +760,7 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
tr->Draw(fillCommand.Data(), fillCondition.Data(), "goff");
//spectrum analysis
fitControl = PeaksFitting(hRaw, "Q", fFitMinSigma);
fitControl = PeaksFitting(hRaw, "Q V v", fFitMinSigma);
Info("CalculateCalibParameters", "Value of fitControl is: %d", fitControl); //ok
//incorrectly treated spectrum output
......@@ -782,7 +783,8 @@ Bool_t AculCalibration::CalculateCalibParameters(const char* inputfile, const ch
//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
calGraph->SetPoint(j, fPeak[j], fEnergy[j]); //calibration graph filling
printf("\tPeak\t%f and energy\t%f\n", fPeak[j], fEnergy[j]);
}//for
calGraph->Fit(calFunction, "Q", "goff", 0, 4096); //omezit hlasitost fitovani, udelat volitelne, dodelat volby rozsahu
outcalfile
......@@ -1101,13 +1103,77 @@ void AculCalibration::FindEnergyPeaks(TCanvas *c1, const char* ifile, const char
c1->cd(i+1);
PeaksFitting(hWork);
hWork->Draw();
ofile<<i<<"\t";
for(Int_t j=0; j<kRaNOPEAKS; j++) {
ofile << fPeak[j] <<"\t";
}
ofile<<endl;
}
ofile.close();
}
void AculCalibration::FindAverageEnergies(const char* ifile, const char* outfile) {
TString iFile = ifile;
TFile *fr = new TFile(iFile.Data());
if ( !fr->IsOpen() ) {
Error("FindAverageEnergies", "File %s was not opened and won't be processed.", iFile.Data());
return;
}
ofile << i <<"\t"<< fPeak[0] <<"\t"<< fPeak[1] <<"\t"<< fPeak[2] <<"\t"<< fPeak[3] <<endl;
TList *histList = fr->GetListOfKeys();
Info("FindAverageEnergies", "%d keys found in file %s.", histList->GetEntries(), fr->GetName());
//creation of output text file with average values of peak energies in MeV
TString workFile = outfile;
ofstream ofile;
ofile.open(workFile.Data());
if (!ofile.is_open()) {
Error("PeaksFitting", "Output file %s was not opened", workFile.Data());
return;
}
TH1 *hWork = 0;
Double_t hArray[histList->GetEntries()][kRaNOPEAKS];
//TString hSumName;
Double_t hSumE1 = 0.;
Double_t hAvrE1 = 0.;
Double_t hSumE2 = 0.;
Double_t hAvrE2 = 0.;
Double_t hSumE3 = 0.;
Double_t hAvrE3 = 0.;
Double_t hSumE4 = 0.;
Double_t hAvrE4 = 0.;
// c1->Clear();
// c1->Divide(6, 6);
for (Int_t i = 0; i < histList->GetEntries(); i++) {
fr->GetObject(histList->At(i)->GetName(), hWork);
PeaksFitting(hWork);
for(Int_t j = 0; j < kRaNOPEAKS; j++) {
hArray[i][j] = fPeak[j];
if(fPeak[j]==0.){
Error("FindAverageEnergies", "No peak in channel %i !", histList->GetEntries());
}
//hSumName.Form("hSumE%i",j);
}
hSumE1 += hArray[i][0];
hSumE2 += hArray[i][1];
hSumE3 += hArray[i][2];
hSumE4 += hArray[i][3];
// std::cout<<"i "<<i<<" hSumE1 "<<hSumE1<<std::endl;
}
hAvrE1 = hSumE1/histList->GetEntries();
hAvrE2 = hSumE2/histList->GetEntries();
hAvrE3 = hSumE3/histList->GetEntries();
hAvrE4 = hSumE4/histList->GetEntries();
ofile <<"Average energies are:\t"<<hAvrE1<<"\t"<<hAvrE2<<"\t"<<hAvrE3<<"\t"<<hAvrE4<<std::endl;
ofile.close();
}
}
void AculCalibration::ShowAnalyzedSpectra(const char *filename, TCanvas* fittedRawCanvas, Int_t xaxismin, Int_t xaxismax, Int_t subaddress)
{
......@@ -1334,7 +1400,9 @@ void AculCalibration::SetInputParameters() {
if ( strcmp(identificator, "noPeaks") == 0 ) {
kRaNOPEAKS = static_cast<Int_t>(atoi(parameter));
fEnergyInput.Set(kRaNOPEAKS);
fEnergyInput.Set(kRaNOPEAKS);
fEnergy.Set(kRaNOPEAKS);
fPeak.Set(kRaNOPEAKS);
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
fipr.getline(line, lineLength);
sscanf(line, "%s", parameter);
......
AculData/AculCalibration.h
View file @ 3a55bcf8
......@@ -176,10 +176,17 @@ public:
void CalibrateRawSpectra(const char* inputfile, const char* block, const char* treename, Int_t lowerchannel = 0, Int_t upperchannel = 4095, Int_t nEBins = 1000, Int_t lowersubaddress = 0, Int_t uppersubaddress = ADDRESSNUMBER-1);
void FindEnergyPeaks(TCanvas *c1, const char* ifile, const char* outfile);
// 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);
void FindEnergyPeaks(TCanvas *c1, const char* ifile, const char* outfile);
// Outputs calibrated energies for each channel in txt file
//
//
//
void FindAverageEnergies(const char* ifile, const char* outfile);
// Outputs average values of calibrated energies for the whole detector in txt file
//
//
//
Int_t PeaksFitting(TH1* hSpectrum, Option_t* option = "", Double_t sigmamin = 2);
//return values:
......
macros/calibration_CsI/calibrationSQ13Alpha.cxx
View file @ 3a55bcf8
{
// gSystem->Load("CsICalib_C.so");
gSystem->Load("../../libAculData.so");
gSystem->Load("/home/dariak/AculUtils/libAculData.so");
//CsICalib cal;
AculCalibCsI cal("parameters/SQ13Alpha.par");
......@@ -13,14 +13,17 @@
TCanvas *c2 = new TCanvas("c2", "Plain");
TCanvas *c3 = new TCanvas("c3", "Plain");
TCanvas *c4 = new TCanvas("c4", "Plain");
// TCanvas *myCanv = new TCanvas("myCanv", "Let's see what we see");
cal.PrintTrees();
cal.PrintFiles();
cal.PrintCuts();
// cal.PrintTrees();
// cal.PrintFiles();
// cal.PrintCuts();
// cal.PrintPeakRanges();
// cal.ReadParFile("ranges.dat");
// cal.PrintParameters();
// cal.PrintParameters("v");
// cal.PrintTrees();
// cal.DrawVariable("SQ13", 0, c1);
// cal.DrawVariable("SQ13", 1, c1);
......@@ -28,10 +31,10 @@
// cal.DrawVariable("TDC", 2, c2);
// cal.DrawVariable("QDC", 3, c3);
// cal.DrawBeam(c4, 4, "SQ13");
// cal.DrawBeam(myCanv, 4, "SQ13");
//return;
// cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16");
// cal.DrawVariableCut("SQ13", 0, c1, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp");
// cal.DrawVariableCut("SQ13", 0, c2, "cutsSQ13Alpha16");
// cal.DrawVariableCut("SQ13", 0, c2, "cutsSQ13Alpha16", "cutSQ13Alpha16Amp");
//return;
// cal.DrawVariableCut("SQ13", 1, c2, "cutSQ13Alpha21", "cutSQ13Alpha21Amp");
// cal.DrawVariableCut("SQ13", 2, c3, "cutSQ13Alpha26", "cutSQ13Alpha26Amp");
......@@ -50,7 +53,7 @@
// cal.SaveClbGraphs("gSQ13Alpha2points.root", "RECREATE");
// cal.SaveClbGraphs("gSQ13AlphaAlt.root", "RECREATE");
// cal.WriteClbParameters("calSQ13Alpha.clb");
cal.WriteClbParameters("calSQ13Alpha.clb");
// cal.WriteClbParameters("calSQ13Alpha2points.clb");
// cal.WriteClbParameters("calSQ13AlphaAlt.clb");
......
macros/calibration_CsI/calibrationSQ23Alpha.cxx
View file @ 3a55bcf8
{
gSystem->Load("CsICalib_C.so");
gSystem->Load("/home/dariak/AculUtils/libAculData.so");
// CsICalib cal;
CsICalib cal("parameters/SQ23Alpha.par");
AculCalibCsI cal("parameters/SQ23Alpha.par");
// cal.PrintCuts();
......@@ -13,8 +13,8 @@
TCanvas *c2 = new TCanvas("c2", "Plain");
TCanvas *c3 = new TCanvas("c3", "Plain");
TCanvas *c4 = new TCanvas("c4", "Plain");
cal.DrawBeam(c4, 4, "SQ23");
cal.PrintParameters();
// cal.DrawBeam(c4, 4, "SQ23");
return;
......