Commit b2df0db4 authored by Vratislav Chudoba's avatar Vratislav Chudoba

Parameters for callibration of scintillator detectors have been somehow enhanced.

parent c1566d0a
...@@ -51,10 +51,12 @@ public: ...@@ -51,10 +51,12 @@ public:
virtual const char* GetFileName(Int_t i) {return 0;}; virtual const char* GetFileName(Int_t i) {return 0;};
virtual Int_t GetNoEPoints() {return 0;} virtual Int_t GetNoEPoints() {return 0;}
virtual Double_t GetCalEnergy(Int_t i) {return 0.;}; virtual Double_t GetCalEnergy(Int_t i) {return 0.;};
virtual const char* GetCutsFileName() {return 0;}
virtual Int_t GetNoCuts() {return 0;} //todo check functions with cuts
virtual TCutG* GetCut(Int_t i) {return 0;}; // virtual const char* GetCutsFileName() {return 0;}
virtual TCutG* GetCut(const char* cutName) {return 0;}; // virtual Int_t GetNoCuts() {return 0;}
// virtual TCutG* GetCut(Int_t i) {return 0;};
virtual TCutG* GetCut(const char* cutName, Int_t treeID = 0) {return 0;};
virtual Int_t GetMinChannel(Int_t energy, Int_t crystal) {return 0;}; virtual Int_t GetMinChannel(Int_t energy, Int_t crystal) {return 0;};
virtual Int_t GetMaxChannel(Int_t energy, Int_t crystal) {return 0;}; virtual Int_t GetMaxChannel(Int_t energy, Int_t crystal) {return 0;};
......
...@@ -27,7 +27,7 @@ AculCalParsScint::~AculCalParsScint() { ...@@ -27,7 +27,7 @@ AculCalParsScint::~AculCalParsScint() {
void AculCalParsScint::Init() { void AculCalParsScint::Init() {
SetPars(); SetPars();
LoadCuts(); // LoadCuts();
} }
void AculCalParsScint::SetPars() { void AculCalParsScint::SetPars() {
...@@ -88,7 +88,7 @@ void AculCalParsScint::SetPars() { ...@@ -88,7 +88,7 @@ void AculCalParsScint::SetPars() {
continue; continue;
} }
if ( line.BeginsWith("cutFile") ) { /*if ( line.BeginsWith("cutFile") ) {
sscanf(line.Data(), "%*s %s %d", fname, &fNoCuts); sscanf(line.Data(), "%*s %s %d", fname, &fNoCuts);
fCutsFileName = fname; fCutsFileName = fname;
for (Int_t j = 0; j < fNoCuts; j++) { for (Int_t j = 0; j < fNoCuts; j++) {
...@@ -101,7 +101,7 @@ void AculCalParsScint::SetPars() { ...@@ -101,7 +101,7 @@ void AculCalParsScint::SetPars() {
fCutName.push_back(cname); fCutName.push_back(cname);
} }
continue; continue;
} }*/
if ( line.BeginsWith("detector") ) { if ( line.BeginsWith("detector") ) {
sscanf(line.Data(), "%*s %s %s", det, part); sscanf(line.Data(), "%*s %s %s", det, part);
...@@ -113,14 +113,61 @@ void AculCalParsScint::SetPars() { ...@@ -113,14 +113,61 @@ void AculCalParsScint::SetPars() {
if ( line.BeginsWith("energy") ) { if ( line.BeginsWith("energy") ) {
sscanf(line.Data(), "%*s %lf", &en); sscanf(line.Data(), "%*s %lf", &en);
fFilePars[fEnergyPoints].SetEnergy(en); fFilePars[fEnergyPoints].SetEnergy(en);
do {line.ReadLine(infile);}
while (line.BeginsWith("#") || line.BeginsWith("//"));
if ( line.BeginsWith("cutFile") ) {
Int_t noCuts;
// sscanf(line.Data(), "%*s %s %d", fname, &fNoCuts);
sscanf(line.Data(), "%*s %s %d", fname, &noCuts);
fFilePars[fEnergyPoints].SetNoCuts(noCuts);
fFilePars[fEnergyPoints].SetCutFileName(fname);
cout << line << endl;
cout << noCuts << endl;
// fCutsFileName = fname;
for (Int_t j = 0; j < noCuts; j++) {
line.ReadLine(infile);
cout << line << endl;
if ( line.BeginsWith("#") || line.BeginsWith("//") ) {
j--;
continue;
}
sscanf(line, "%s", cname);
fFilePars[fEnergyPoints].SetCutName(j, cname);
// fCutName.push_back(cname);
}
fFilePars[fEnergyPoints].LoadCuts();
line.ReadLine(infile);
cout << line << endl;
// continue;
}
// cout << fNoCrystals << endl;
for (Int_t j = 0; j < fNoCrystals; j++) {
// line.ReadLine(infile);
cout << j << endl;
if ( line.BeginsWith("#") || line.BeginsWith("//") ) {
j--;
line.ReadLine(infile);
continue;
}
sscanf(line.Data(), "%d %d %d", &j, &min, &max);
// cout << j << endl;
fFilePars.at(fEnergyPoints).SetPeakMin(j, min);
fFilePars.at(fEnergyPoints).SetPeakMax(j, max);
line.ReadLine(infile);
}
fEnergyPoints++; fEnergyPoints++;
continue; continue;
} }
// cout << "aklsjdlaksjd" << endl;
sscanf(line.Data(), "%d %d %d", &i, &min, &max); // sscanf(line.Data(), "%d %d %d", &i, &min, &max);
fFilePars.at(fEnergyPoints-1).SetPeakMin(i, min); // fFilePars.at(fEnergyPoints-1).SetPeakMin(i, min);
fFilePars.at(fEnergyPoints-1).SetPeakMax(i, max); // fFilePars.at(fEnergyPoints-1).SetPeakMax(i, max);
}//while }//while
...@@ -141,21 +188,27 @@ void AculCalParsScint::PrintParameters(const char* option) { ...@@ -141,21 +188,27 @@ void AculCalParsScint::PrintParameters(const char* option) {
} }
cout << "\tInput energies (" << fNoFiles << "):" << endl; cout << "\tInput energies (" << fNoFiles << "):" << endl;
for (Int_t i = 0; i <= (Int_t)fFilePars.size(); i++) { for (Int_t i = 0; i < (Int_t)fFilePars.size(); i++) {
cout << "\t " << fFilePars[i].GetEnergy() << " MeV" << endl; cout << "\t " << fFilePars[i].GetEnergy() << " MeV" << endl;
} }
cout << "\tInput file with " << fNoCuts << " cuts: \"" << fCutsFileName << "\"" << endl; // cout << "\tInput file with " << fNoCuts << " cuts: \"" << fCutsFileName << "\"" << endl;
if (opt.Contains("all")) { // if (opt.Contains("all")) {
for (Int_t i = 0; i < (Int_t)fCutName.size(); i++) { // for (Int_t i = 0; i < (Int_t)fCutName.size(); i++) {
cout << "\t cut: \"" << fCutName[i] << "\"" << endl; // cout << "\t cut: \"" << fCutName[i] << "\"" << endl;
} // }
} // }
if (!opt.Contains("all")) return;
cout << "\tPeak limits for particular channels and energies:" << endl;
for (Int_t k = 0; k < (Int_t)fFilePars.size(); k++) { for (Int_t k = 0; k < (Int_t)fFilePars.size(); k++) {
cout << "\t" << fFilePars[k].GetNoCuts() << " cuts from file \"" << fFilePars[k].GetCutFileName() << "\":" << endl;
for (Int_t i = 0; i < fFilePars[k].GetNoCuts(); i++) {
cout << "\t " << fFilePars[k].GetCut(i)->GetName() << endl;
}
if (!opt.Contains("all")) continue;
cout << "\tPeak limits for particular channels and energies:" << endl;
cout << "\t Set number: " << k << "; energy: " << fFilePars[k].GetEnergy() << " MeV" << endl; cout << "\t Set number: " << k << "; energy: " << fFilePars[k].GetEnergy() << " MeV" << endl;
for (Int_t i = 0; i < (Int_t)fFilePars[k].GetNoCrystals(); i++) { for (Int_t i = 0; i < (Int_t)fFilePars[k].GetNoCrystals(); i++) {
cout << "\t\t " << fFilePars[k].GetPeakMin(i) << "\t" << fFilePars[k].GetPeakMax(i) << endl; cout << "\t\t " << fFilePars[k].GetPeakMin(i) << "\t" << fFilePars[k].GetPeakMax(i) << endl;
...@@ -167,6 +220,8 @@ void AculCalParsScint::PrintParameters(const char* option) { ...@@ -167,6 +220,8 @@ void AculCalParsScint::PrintParameters(const char* option) {
void AculCalParsScint::Reset() { void AculCalParsScint::Reset() {
fNoCrystals = 0;
fParFileName = ""; fParFileName = "";
fDetName = ""; fDetName = "";
...@@ -176,10 +231,10 @@ void AculCalParsScint::Reset() { ...@@ -176,10 +231,10 @@ void AculCalParsScint::Reset() {
fFilePars.clear(); fFilePars.clear();
fEnergyPoints = 0; fEnergyPoints = 0;
fCutsFileName = ""; // fCutsFileName = "";
fNoCuts = 0; //number of cuts // fNoCuts = 0; //number of cuts
fCutName.clear(); // fCutName.clear();
return; return;
} }
...@@ -193,29 +248,40 @@ const char* AculCalParsScint::GetFileName(Int_t i) { ...@@ -193,29 +248,40 @@ const char* AculCalParsScint::GetFileName(Int_t i) {
return fFilePars[i].GetRawFileName(); return fFilePars[i].GetRawFileName();
} }
const char* AculCalParsScint::GetCutName(Int_t i) { /*const char* AculCalParsScint::GetCutName(Int_t i) {
if ( i > (Int_t)fCutName.size()-1 ) { if ( i > (Int_t)fCutName.size()-1 ) {
cerr << "\"AculCalParsScint::GetCutName\" index i cannot be higher than " << fCutName.size() - 1 << endl; cerr << "\"AculCalParsScint::GetCutName\" index i cannot be higher than " << fCutName.size() - 1 << endl;
return 0; return 0;
} }
return fCutName[i].Data(); return fCutName[i].Data();
} }*/
TCutG* AculCalParsScint::GetCut(Int_t i) { /*TCutG* AculCalParsScint::GetCut(Int_t i) {
if ( i >= (Int_t)fCuts.size() ) { if ( i >= (Int_t)fCuts.size() ) {
cerr << "\"AculCalParsScint::GetCut\" index i cannot be higher than " << fCuts.size() - 1 << endl; cerr << "\"AculCalParsScint::GetCut\" index i cannot be higher than " << fCuts.size() - 1 << endl;
return 0; return 0;
} }
return &fCuts[i]; return &fCuts[i];
} }*/
TCutG* AculCalParsScint::GetCut(const char* cutName) { TCutG* AculCalParsScint::GetCut(const char* cutName, Int_t treeID) {
const TString cName = cutName; const TString cName = cutName;
for (Int_t i = 0; i <= (Int_t)fCuts.size(); i++) { if (treeID >= (Int_t)fFilePars.size()) {
if (cName.EqualTo(fCutName[i])) { return &fCuts[i]; } cerr << "\"AculCalParsScint::GetCut\" treeID does not exist. It cannot be higher than "
<< fFilePars.size() << "." << endl;
return 0;
}
// TCutG *curCut = 0;
for (Int_t i = 0; i < fFilePars[treeID].GetNoCuts(); i++) {
// if (cName.EqualTo(fFilePars[treeID])) { return &fCuts[i]; }
if ( cName.EqualTo( fFilePars[treeID].GetCutName(i) ) ) {
return fFilePars[treeID].GetCut(i);
// return curCut;
}
} }
...@@ -266,7 +332,7 @@ Int_t AculCalParsScint::GetMaxChannel(Int_t energy, Int_t crystal) { ...@@ -266,7 +332,7 @@ Int_t AculCalParsScint::GetMaxChannel(Int_t energy, Int_t crystal) {
return fFilePars[energy].GetPeakMax(crystal); return fFilePars[energy].GetPeakMax(crystal);
} }
void AculCalParsScint::LoadCuts() { /*void AculCalParsScint::LoadCuts() {
if (fCutsFileName.Length() == 0) { if (fCutsFileName.Length() == 0) {
printf("AculCalParsScint::LoadCuts: Name of file (*.root) with cuts was not provided.\n"); printf("AculCalParsScint::LoadCuts: Name of file (*.root) with cuts was not provided.\n");
...@@ -291,4 +357,4 @@ void AculCalParsScint::LoadCuts() { ...@@ -291,4 +357,4 @@ void AculCalParsScint::LoadCuts() {
fCuts.push_back(*currentCut); fCuts.push_back(*currentCut);
} }
} }*/
...@@ -26,13 +26,13 @@ private: ...@@ -26,13 +26,13 @@ private:
vector<AculCalParsScintFile> fFilePars; vector<AculCalParsScintFile> fFilePars;
//todo following parameters should be probably moved to AculCalParsScintFile //todo following parameters should be probably moved to AculCalParsScintFile
//as far as they are related to each file separately //as far as they are related to each file separately
vector<TString> fCutName; // vector<TString> fCutName;
vector<TCutG> fCuts; // vector<TCutG> fCuts;
Int_t fEnergyPoints; Int_t fEnergyPoints;
TString fCutsFileName; // TString fCutsFileName;
Int_t fNoCuts; //number of cuts // Int_t fNoCuts; //number of cuts
public: public:
AculCalParsScint(); AculCalParsScint();
...@@ -57,15 +57,15 @@ public: ...@@ -57,15 +57,15 @@ public:
const char* GetCutName(Int_t i); const char* GetCutName(Int_t i);
Int_t GetNoEPoints() {return fEnergyPoints;} Int_t GetNoEPoints() {return fEnergyPoints;}
Double_t GetCalEnergy(Int_t i); Double_t GetCalEnergy(Int_t i);
const char* GetCutsFileName() {return fCutsFileName.Data();} // const char* GetCutsFileName() {return fCutsFileName.Data();}
Int_t GetNoCuts() {return fNoCuts;} // Int_t GetNoCuts() {return fNoCuts;}
TCutG* GetCut(Int_t i); // TCutG* GetCut(Int_t i);
TCutG* GetCut(const char* cutName); TCutG* GetCut(const char* cutName, Int_t treeID = 0);
Int_t GetMinChannel(Int_t energy, Int_t crystal); Int_t GetMinChannel(Int_t energy, Int_t crystal);
Int_t GetMaxChannel(Int_t energy, Int_t crystal); Int_t GetMaxChannel(Int_t energy, Int_t crystal);
//private: //private:
void LoadCuts(); // void LoadCuts();
protected: protected:
void SetPars(); void SetPars();
......
/*
* AculCalParsScintFile.cpp
*
* Created on: Oct 25, 2016
* Author: vratik
*/
#include "AculCalParsScintFile.h"
AculCalParsScintFile::AculCalParsScintFile() {
// TODO Auto-generated constructor stub
// fFileID = 0;
fNoCrystals = 0;
fE = 0;
Init();
}
AculCalParsScintFile::~AculCalParsScintFile() {
// TODO Auto-generated destructor stub
}
void AculCalParsScintFile::Init() {
// fNoCrystals = 0;
// fE = 0;
fPeakMin.resize(fNoCrystals);
fPeakMax.resize(fNoCrystals);
fCutName.resize(0);
}
const char* AculCalParsScintFile::GetCutName(Int_t i) {
if ( i >= (Int_t)fCutName.size() ) {
cout << "\"AculCalParsScintFile::GetCutName\" Cut number " << i
<< " is not available. Maximal i can be set as "
<< fCutName.size() - 1 << "." << endl;
}
return fCutName[i].Data();
}
TCutG* AculCalParsScintFile::GetCut(Int_t i) {
if ( i >= (Int_t)fCut.size() ) {
cerr << "\"AculCalParsScintFile::GetCut\" index i cannot be higher than " << fCut.size() - 1 << endl;
return 0;
}
return &fCut[i];
}
TCutG* AculCalParsScintFile::GetCut(const char* cutName) {
const TString cName = cutName;
for (Int_t i = 0; i <= (Int_t)fCut.size(); i++) {
if (cName.EqualTo(fCutName[i])) { return &fCut[i]; }
}
cerr << "\"AculCalParsScint::GetCut\" cut \"" << cutName << "\" was not found." << endl;
return 0;
}
Int_t AculCalParsScintFile::GetPeakMin(Int_t crystalID) {
if ( crystalID >= (Int_t)fPeakMin.size() ) {
cerr << "\"AculCalParsScintFile::GetPeakMin\" index \"crystal\" cannot be higher than "
<< fPeakMin.size() - 1 << endl;
return 0;
}
return fPeakMin[crystalID];
}
Int_t AculCalParsScintFile::GetPeakMax(Int_t crystalID) {
if ( crystalID >= (Int_t)fPeakMax.size() ) {
cerr << "\"AculCalParsScintFile::GetPeakMax\" index \"crystal\" cannot be higher than "
<< fPeakMax.size() - 1 << endl;
return 0;
}
return fPeakMax[crystalID];
}
void AculCalParsScintFile::SetRawFileName(const char* filename) {
fFileName = filename;
// fFileID++;
}
void AculCalParsScintFile::SetCutFileName(const char* filename) {
fCutsFileName = filename;
}
void AculCalParsScintFile::SetNoCuts(Int_t noCuts) {
fCutName.resize(noCuts);
}
void AculCalParsScintFile::SetCutName(Int_t i, const char* cutname) {
if ( i >= (Int_t)fCutName.size() ) {
cout << "\"AculCalParsScintFile::SetCutName\" vector for cut names is not enough large. Maximal i can be set as "
<< fCutName.size() - 1 << "." << endl;
}
fCutName[i] = cutname;
}
void AculCalParsScintFile::SetPeakMin(Int_t crystalID, Double_t limit) {
if ( crystalID>=(Int_t)fPeakMin.size() ) {
cerr << "\"AculCalParsScintFile::SetPeakMin\" maximal allowed crystal ID is " << fPeakMin.size() << "."<< endl;
return;
}
fPeakMin[crystalID] = limit;
}
void AculCalParsScintFile::SetPeakMax(Int_t crystalID, Double_t limit) {
if ( crystalID>=(Int_t)fPeakMax.size() ) {
cerr << "\"AculCalParsScintFile::SetPeakMax\" maximal allowed crystal ID is " << fPeakMax.size() << "."<< endl;
return;
}
fPeakMax[crystalID] = limit;
}
void AculCalParsScintFile::LoadCuts() {
if (fCutsFileName.Length() == 0) {
printf("AculCalParsScintFile::LoadCuts: Name of file (*.root) with cuts was not provided.\n");
printf("AculCalParsScintFile::LoadCuts: No cuts has been loaded.\n");
return;
}
TFile cutFile(fCutsFileName.Data(), "READ");
if(!cutFile.IsOpen()) {
cerr << "\"AculCalParsScintFile::LoadCuts\" File " << fCutsFileName.Data()
<< " was not open and no cuts were loaded." << endl;
return;
}
for (Int_t i = 0; i < (Int_t)fCutName.size(); i++) {
TCutG *currentCut = (TCutG*)cutFile.Get(fCutName[i]);
if (!currentCut) {
cout << "\"AculCalParsScintFile::LoadCuts\" Cut \"" << fCutName[i]
<< "\" was not found in file " << fCutsFileName << "." << endl;
continue;
}
fCut.push_back(*currentCut);
}
}
/*void AculCalParsScintFile::Reset() {
}*/
/*
* AculCalParsScintFile.h
*
* Created on: Oct 25, 2016
* Author: vratik
*/
#ifndef ACULCALIB_ACULCALPARSSCINTFILE_H_
#define ACULCALIB_ACULCALPARSSCINTFILE_H_
#include <vector>
#include <iostream>
#include "TFile.h"
#include "TCutG.h"
#include "TString.h"
using std::vector;
using std::cout;
using std::endl;
using std::cerr;
class AculCalParsScintFile {
private:
// UInt_t fFileID;
Int_t fNoCrystals;
Double_t fE;
// TCutG fCuts;
TString fFileName;
TString fCutsFileName;
vector<TString> fCutName;
vector<TCutG> fCut;
vector<Int_t> fPeakMin;
vector<Int_t> fPeakMax;
public:
AculCalParsScintFile();
virtual ~AculCalParsScintFile();
ClassDef(AculCalParsScintFile, 1);
void Init();
Int_t GetNoCrystals() {return fNoCrystals;}
const char* GetRawFileName() {return fFileName.Data();}
const char* GetCutFileName() {return fCutsFileName.Data();}
Double_t GetEnergy() {return fE;}
Int_t GetNoCuts() {return (Int_t)fCut.size();}
const char* GetCutName(Int_t i);
TCutG* GetCut(Int_t i);
TCutG* GetCut(const char* cutName);
Int_t GetPeakMin(Int_t crystalID);
Int_t GetPeakMax(Int_t crystalID);
void SetNoCrystals(Int_t nocrystals) {fNoCrystals = nocrystals;};
void SetEnergy(Double_t energy) {fE = energy;};
void SetRawFileName(const char* filename);
void SetCutFileName(const char* filename);
void SetNoCuts(Int_t noCuts);
void SetCutName(Int_t i, const char* cutname);
void SetPeakMin(Int_t crystalID, Double_t limit);
void SetPeakMax(Int_t crystalID, Double_t limit);
void LoadCuts();
// void Reset();
private:
// void LoadCuts();
};
#endif /* ACULCALIB_ACULCALPARSSCINTFILE_H_ */
/*
* AculCalParsSi.cpp
*
* Created on: Oct 26, 2016
* Author: vratik
*/
#include "AculCalParsSi.h"
AculCalParsSi::AculCalParsSi() {
// TODO Auto-generated constructor stub
Reset();
}
AculCalParsSi::AculCalParsSi(const char* parFile) {
SetParFile(parFile);
Init();
}
AculCalParsSi::~AculCalParsSi() {
// TODO Auto-generated destructor stub
}
void AculCalParsSi::Init() {
SetELosses();
SetPars();
SetCalEnergies();
}
void AculCalParsSi::Reset() {
fParFileName = "";
kRaNOPEAKS = 0;
fEnergyInput.clear();
fLowerChannel = 0.;
fUpperChannel = 0.;
fLowerPeakRelativeHight = 0.;
fUpperPeakRelativeHight = 0.;
fPeakPositionTolerance = 0.;
fFitFuncLineWidth = 1;
fFitMinSigma = 0.;
fFitPeakThreshold = 0.;
fDeadLayer = 0.;
return;
}
void AculCalParsSi::SetPars() {
if (fParFileName.Length()==0) {
cerr << "\'AculCalibration::SetInputsParameters\" File with input parameters was not set." << endl;
return;
}
const Int_t lineLength = 400;
Char_t line[lineLength];
Char_t parameter[100];
Char_t identificator[100];
ifstream fipr;
fipr.open(fParFileName.Data());
if (!fipr.is_open()) {
cerr << "\"AculCalibration::SetInputsParameters\" File with input parameters \""
<< fParFileName << "\" was not opened." << endl;
return;
}
cout << "\"AculCalibration::SetInputsParameters\" File with input parameters \""
<< fParFileName << "\" will be processed." << endl;
while (!fipr.eof()) {
fipr.getline(line, lineLength);
if (strlen(line) < 2) {
continue;
}
sscanf(line, "%s %s", parameter, identificator);
if ( strcmp(identificator, "noPeaks") == 0 ) {
kRaNOPEAKS = static_cast<Int_t>(atoi(parameter));
fEnergyInput.resize(kRaNOPEAKS);
fEtab.resize(kRaNOPEAKS);
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
fipr.getline(line, lineLength);
sscanf(line, "%s", parameter);
fEnergyInput[i] = static_cast<Double_t>(atof(parameter));
}
continue;
}//if
if ( strcmp(identificator, "lowerChannel") == 0 ) {
sscanf(line, "%s", parameter);
fLowerChannel = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "upperChannel") == 0 ) {
sscanf(line, "%s", parameter);
fUpperChannel = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "lowerPeakHight") == 0 ) {
sscanf(line, "%s", parameter);
fLowerPeakRelativeHight = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "upperPeakHight") == 0 ) {
sscanf(line, "%s", parameter);
fUpperPeakRelativeHight = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "peakPositionTolerance") == 0 ) {
sscanf(line, "%s", parameter);
fPeakPositionTolerance = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "fitFunctionLineWidth") == 0 ) {
sscanf(line, "%s", parameter);
fFitFuncLineWidth = static_cast<Width_t>(atoi(parameter));
}
if ( strcmp(identificator, "minFitSigma") == 0 ) {
sscanf(line, "%s", parameter);
fFitMinSigma = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "fitHightThreshold") == 0 ) {
sscanf(line, "%s", parameter);
fFitPeakThreshold = static_cast<Double_t>(atof(parameter));
}
if ( strcmp(identificator, "deadLayer") == 0 ) {
sscanf(line, "%s", parameter);
fDeadLayer = static_cast<Double_t>(atof(parameter));
}
}
fipr.close();
return;
}
void AculCalParsSi::PrintParameters(const char* option)
{
//print alpha source parameters
cout << "AculCalibration::PrintInputParameters:" << endl;
cout << "\tNumber of peaks: " << kRaNOPEAKS << endl;
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
cout << "\t\tfEnergyInput[" << i << "] = " << fEnergyInput[i] << " MeV" << endl;
}
cout << "\tEnergies used for calibration:" << endl;
cout << "\t(deadLayer: " << fDeadLayer << " mcm)" << endl;
for (Int_t i = 0; i < kRaNOPEAKS; i++) {
cout << "\t\tfE[" << i << "] = " << fEtab[i] << " MeV" << endl;
}
cout << "\tlowerChannel: " << fLowerChannel << "; upperChannel: " << fUpperChannel << ";" << endl;
cout << "\tlowerPeakHight: " << fLowerPeakRelativeHight << "; upperPeakHight: " << fUpperPeakRelativeHight << ";" << endl;
cout << "\tfitHightThreshold: " << fFitPeakThreshold << "; minFitSigma: " << fFitMinSigma << ";" << endl;
cout << "\tpeakPositionTolerance: " << fPeakPositionTolerance << ";" << endl;
cout << "\tfitFunctionLineWidth: " << fFitFuncLineWidth << ";" << endl;
return;
}
void AculCalParsSi::SetELosses() {
Info("AculCalibration::SetELosses", "Combination of aplha particle with silicon material only.");
fAlphaSi.SetEL(1, 2.321); // density in g/cm3
fAlphaSi.AddEL(14., 28.086, 1); //Z, mass
// mSi.SetZP(1., 1.); //protons
fAlphaSi.SetZP(2., 4.); //alphas, Z, A
fAlphaSi.SetEtab(100000, 200.); // ?, MeV calculate ranges
fAlphaSi.SetDeltaEtab(300);
}
void AculCalParsSi::SetCalEnergies() {
if (fDeadLayer<=0.) {
Warning("AculCalibration::SetCalEnergies", "Dead layer was set equal or less than 0.");
for(Int_t i = 0; i < kRaNOPEAKS; i++) {
fEtab[i] = fEnergyInput[i];
}
Info("AculCalibration::SetCalEnergies", "Energies used for calibration are the same as input file.");
return;
}
for(Int_t i = 0; i < kRaNOPEAKS; i++) {
fEtab[i] = fAlphaSi.GetE(fEnergyInput[i], fDeadLayer);
}
Info("AculCalibration::SetCalEnergies", "Energies used for calibration considering %f mcm dead layer were set.", fDeadLayer);
return;
}
/*
* AculCalParsSi.h
*
* Created on: Oct 26, 2016
* Author: vratik
*/
#ifndef ACULCALIB_ACULCALPARSSI_H_
#define ACULCALIB_ACULCALPARSSI_H_
#include "AculCalPars.h"
#include "../TELoss/TELoss.h"
class AculCalParsSi: public AculCalPars {
private:
Int_t kRaNOPEAKS;
// TArrayD fEnergy; //energy after passing through deadlayer
vector<Double_t> fEnergyInput; //incidental energy, set from .par
Double_t fLowerChannel;
Double_t fUpperChannel;
Double_t fLowerPeakRelativeHight; //pouziva se, private
Double_t fUpperPeakRelativeHight; //pouziva se, private, nastavit nenulovou prednastavenou hodnotu
Double_t fPeakPositionTolerance; //pouziva se, private
Width_t fFitFuncLineWidth; //private
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
Double_t fDeadLayer; //in mcm; it will be used for recalculation of energies outcomming from the alpha source
TELoss fAlphaSi;
public:
AculCalParsSi();
AculCalParsSi(const char* parFile);
virtual ~AculCalParsSi();
ClassDef(AculCalParsSi, 1);
void Init();
// Function which loads text file containing parameters for calibration
//
// -inputparfile: file containing information on calibration source
// -In file with the data must be preserved systematic
// -There cannot be whitelines (probably)
//
//Example of "parforcal.par" file
//................................................................
//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 //
//................................................................
void Reset();
void PrintParameters(const char* option = "");
//I hope this is selfunderstanding function
private:
void SetPars();
void SetELosses();
void SetCalEnergies();
};
#endif /* ACULCALIB_ACULCALPARSSI_H_ */
...@@ -160,7 +160,9 @@ void AculCalibScint::DrawBeam(TCanvas *canvas, Int_t file, const char* variable) ...@@ -160,7 +160,9 @@ void AculCalibScint::DrawBeam(TCanvas *canvas, Int_t file, const char* variable)
TString var; TString var;
TString con; TString con;
cout << "cuts " << fPars->GetNoCuts() << endl; // cout << "cuts " << fPars->GetNoCuts() << endl;
cout << "this function has to be checked" << endl;
return;
for (Int_t i = 0; i < noFiles; i++) { for (Int_t i = 0; i < noFiles; i++) {
canvas->cd(i+1); canvas->cd(i+1);
...@@ -173,23 +175,23 @@ void AculCalibScint::DrawBeam(TCanvas *canvas, Int_t file, const char* variable) ...@@ -173,23 +175,23 @@ void AculCalibScint::DrawBeam(TCanvas *canvas, Int_t file, const char* variable)
// curTree->Draw("QDC[0]+QDC[1]:(TDC[2]+TDC[3])/2. - (TDC[0]+TDC[1])/2.", "", "cont"); // curTree->Draw("QDC[0]+QDC[1]:(TDC[2]+TDC[3])/2. - (TDC[0]+TDC[1])/2.", "", "cont");
// cout << "aksjda\t" << energyPoints << endl; // cout << "aksjda\t" << energyPoints << endl;
for (Int_t j = 0; j < fPars->GetNoRawFiles(); j++) { for (Int_t j = 0; j < fPars->GetNoRawFiles(); j++) {
if ( fPars->GetCut(j) ) { // if ( fPars->GetCut(j) ) {
TCutG *cutTemp = fPars->GetCut(j); // TCutG *cutTemp = fPars->GetCut(j);
curCutG = new TCutG(*cutTemp); // curCutG = new TCutG(*cutTemp);
curCutG->Draw("same"); curCutG->Draw("same");
// printf("AculCalibCsI::DrawBeam: cTOF cut No. %d cannot be drawn, need to repair this function.\n", j); // printf("AculCalibCsI::DrawBeam: cTOF cut No. %d cannot be drawn, need to repair this function.\n", j);
} // }
} }
canvas->cd(noFiles+1+i); canvas->cd(noFiles+1+i);
curTree->Draw("QDC[0]:QDC[1]", "", "cont"); curTree->Draw("QDC[0]:QDC[1]", "", "cont");
for (Int_t j = 0; j < fPars->GetNoRawFiles(); j++) { for (Int_t j = 0; j < fPars->GetNoRawFiles(); j++) {
if ( fPars->GetCut(noFiles + j) ) { // if ( fPars->GetCut(noFiles + j) ) {
TCutG *cutTemp = fPars->GetCut(noFiles + j); // TCutG *cutTemp = fPars->GetCut(noFiles + j);
curCutG = new TCutG(*cutTemp); // curCutG = new TCutG(*cutTemp);
curCutG->Draw("same"); curCutG->Draw("same");
// printf("AculCalibCsI::DrawBeam: cQCD cut No. %d cannot be drawn, need to repair this function.\n", j); // printf("AculCalibCsI::DrawBeam: cQCD cut No. %d cannot be drawn, need to repair this function.\n", j);
} // }
} }
canvas->cd(2*noFiles+1+i); canvas->cd(2*noFiles+1+i);
...@@ -207,7 +209,7 @@ void AculCalibScint::DrawBeam(TCanvas *canvas, Int_t file, const char* variable) ...@@ -207,7 +209,7 @@ void AculCalibScint::DrawBeam(TCanvas *canvas, Int_t file, const char* variable)
} }
} }
void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, Int_t energy, TCanvas *canvas, const char* beamcut, const Int_t nbins, Int_t lowRange) { void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, /*Int_t energy,*/ TCanvas *canvas, const char* beamcut, const Int_t nbins, Int_t lowRange) {
// todo: change this parameter as nonconstant // todo: change this parameter as nonconstant
const Int_t noCrystals = 16; const Int_t noCrystals = 16;
...@@ -217,10 +219,22 @@ void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, Int_t energ ...@@ -217,10 +219,22 @@ void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, Int_t energ
TString con; TString con;
TString hname; TString hname;
TString canvasTitle; TString canvasTitle;
TCutG *currCut = fPars->GetCut(beamcut, treeID);
// if (!fPars->GetCut(beamcut, treeID)) {
// cout << "\"AculCalibScint::GetPeakMean\" cut with name \"" << beamcut << "\" was not found." << endl;
// } else {
// currCut =
// }
// todo manage the canvas name using cut name // todo manage the canvas name using cut name
canvasTitle.Form("variable: %s; tree: %d; cut: %s;", variable, treeID, fPars->GetCut(beamcut)->GetName()); if (currCut) {
canvasTitle.Form("variable: %s; tree: %d; cut: %s;", variable, treeID, currCut->GetName());
canvas->SetTitle(canvasTitle.Data()); canvas->SetTitle(canvasTitle.Data());
} else {
canvasTitle.Form("variable: %s; tree: %d; no cut;", variable, treeID);
canvas->SetTitle(canvasTitle.Data());
}
TTree *curTree = 0; TTree *curTree = 0;
curTree = GetTree(treeID); curTree = GetTree(treeID);
...@@ -242,7 +256,13 @@ void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, Int_t energ ...@@ -242,7 +256,13 @@ void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, Int_t energ
} }
var.Form("%s[%d]>>hfull[%d][%d]", variable, channel, treeID, i); var.Form("%s[%d]>>hfull[%d][%d]", variable, channel, treeID, i);
con.Form("%s[%d]>%d && %s", variable, channel, lowRange, fPars->GetCut(beamcut)->GetName()); if (currCut) {
con.Form("%s[%d]>%d && %s", variable, channel, lowRange, currCut->GetName());
}
else {
// cout << "\"AculCalibScint::GetPeakMean\" cut with name \"" << beamcut << "\" was not found." << endl;
con.Form("%s[%d]>%d", variable, channel, lowRange);
}
canvas->cd(i+1); canvas->cd(i+1);
hname.Form("hfull[%d][%d]", treeID, i); hname.Form("hfull[%d][%d]", treeID, i);
fHistFull[treeID].push_back( new TH1I(hname.Data(), "title", nbins, 0, 4096) ); fHistFull[treeID].push_back( new TH1I(hname.Data(), "title", nbins, 0, 4096) );
...@@ -250,10 +270,18 @@ void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, Int_t energ ...@@ -250,10 +270,18 @@ void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, Int_t energ
curTree->Draw(var.Data(), con.Data()); curTree->Draw(var.Data(), con.Data());
var.Form("%s[%d]>>hcut[%d][%d]", variable, channel, treeID, i); var.Form("%s[%d]>>hcut[%d][%d]", variable, channel, treeID, i);
if (currCut) {
con.Form("%s[%d]>%d && %s[%d]<%d && %s", con.Form("%s[%d]>%d && %s[%d]<%d && %s",
variable, channel, fPars->GetMinChannel(energy, i), variable, channel, fPars->GetMinChannel(treeID, i),
variable, channel, fPars->GetMaxChannel(energy, i), variable, channel, fPars->GetMaxChannel(treeID, i),
fPars->GetCut(beamcut)->GetName()); fPars->GetCut(beamcut, treeID)->GetName());
}
else {
con.Form("%s[%d]>%d && %s[%d]<%d",
variable, channel, fPars->GetMinChannel(treeID, i),
variable, channel, fPars->GetMaxChannel(treeID, i));
}
hname.Form("hcut[%d][%d]", treeID, i); hname.Form("hcut[%d][%d]", treeID, i);
fHistCut[treeID].push_back( new TH1I(hname.Data(), "title", nbins, 0, 4096) ); fHistCut[treeID].push_back( new TH1I(hname.Data(), "title", nbins, 0, 4096) );
fHistCut[treeID][i]->SetLineColor(3); fHistCut[treeID][i]->SetLineColor(3);
......
...@@ -29,16 +29,11 @@ private: ...@@ -29,16 +29,11 @@ private:
TTree **fTrees; TTree **fTrees;
TFile *fCutFile; TFile *fCutFile;
// TClonesArray cutsCol;
//
// TH1I *fHistFull[NOCALFILES][16];
vector< vector<TH1I*> > fHistFull; vector< vector<TH1I*> > fHistFull;
// TH1I *fHistCut[NOCALFILES][16];
vector< vector<TH1I*> > fHistCut; vector< vector<TH1I*> > fHistCut;
// //
// Double_t fMeanPeakPos[NOCALFILES][16];
vector< vector<Double_t> > fMeanPeakPos; vector< vector<Double_t> > fMeanPeakPos;
// Double_t fMeanPeakRMS[NOCALFILES][16];
vector< vector<Double_t> > fMeanPeakRMS; vector< vector<Double_t> > fMeanPeakRMS;
// //
//todo make vector of graphs //todo make vector of graphs
...@@ -64,7 +59,7 @@ public: ...@@ -64,7 +59,7 @@ public:
void DrawVariable(const char* variable, Int_t treeID, TCanvas *canvas, Int_t lowRange = 0, Int_t upRange = 4096); void DrawVariable(const char* variable, Int_t treeID, TCanvas *canvas, Int_t lowRange = 0, Int_t upRange = 4096);
void DrawBeam(TCanvas *canvas, Int_t file, const char* variable); void DrawBeam(TCanvas *canvas, Int_t file, const char* variable);
void GetPeakMean(const char* variable, Int_t treeID, Int_t energy, TCanvas *canvas, const char* beamcut, const Int_t nbins = 4096, Int_t lowRange = 0); void GetPeakMean(const char* variable, Int_t treeID, /*Int_t energy,*/ TCanvas *canvas, const char* beamcut, const Int_t nbins = 4096, Int_t lowRange = 0);
void Calibrate(TCanvas *canvas, Bool_t savefile = 0, const char* filename = "", const char* option = "READ"); void Calibrate(TCanvas *canvas, Bool_t savefile = 0, const char* filename = "", const char* option = "READ");
void WriteClbParameters(const char* filename); void WriteClbParameters(const char* filename);
......
...@@ -7,22 +7,14 @@ data/csi_13_Ealpha16.root ...@@ -7,22 +7,14 @@ data/csi_13_Ealpha16.root
data/csi_13_Ealpha21.root data/csi_13_Ealpha21.root
data/csi_13_Ealpha26.root data/csi_13_Ealpha26.root
data/csi_13_Ealpha30.root data/csi_13_Ealpha30.root
#keyword; filename; number of cuts in file
cutFile cuts/cutsSQ13Alpha.root 8
cutsSQ13Alpha16
cutSQ13Alpha21
#
cutSQ13Alpha26
cutSQ13Alpha30
#
cutSQ13Alpha16Amp
cutSQ13Alpha21Amp
cutSQ13Alpha26Amp
cutSQ13Alpha30Amp
#keyword detector particle #keyword detector particle
detector SQ13 Alpha detector SQ13 Alpha
#keyword energy [MeV/A] #keyword energy [MeV/A]
energy 66.542 energy 66.542
#keyword; filename; number of cuts in file
cutFile cuts/cutsSQ13Alpha.root 2
cutsSQ13Alpha16
cutSQ13Alpha16Amp
#channel minimum maximum #channel minimum maximum
0 650 830 0 650 830
1 920 1150 1 920 1150
...@@ -42,6 +34,11 @@ energy 66.542 ...@@ -42,6 +34,11 @@ energy 66.542
15 800 1000 15 800 1000
#CsI2, energy ~ 21 MeV/A, Alpha #CsI2, energy ~ 21 MeV/A, Alpha
energy 85.437 energy 85.437
#keyword; filename; number of cuts in file
cutFile cuts/cutsSQ13Alpha.root 2
cutSQ13Alpha21
cutSQ13Alpha21Amp
#channel minimum maximum
0 1050 1200 0 1050 1200
1 1350 1550 1 1350 1550
2 1190 1400 2 1190 1400
...@@ -60,6 +57,11 @@ energy 85.437 ...@@ -60,6 +57,11 @@ energy 85.437
15 1150 1350 15 1150 1350
#CsI2, energy ~ 26 MeV/A, Alpha #CsI2, energy ~ 26 MeV/A, Alpha
energy 105.5 energy 105.5
#keyword; filename; number of cuts in file
cutFile cuts/cutsSQ13Alpha.root 2
cutSQ13Alpha26
cutSQ13Alpha26Amp
#channel minimum maximum
0 1300 1600 0 1300 1600
1 1750 2050 1 1750 2050
2 1560 1800 2 1560 1800
...@@ -78,6 +80,11 @@ energy 105.5 ...@@ -78,6 +80,11 @@ energy 105.5
15 1500 1700 15 1500 1700
#CsI2, energy ~ 30 MeV/A, Alpha #CsI2, energy ~ 30 MeV/A, Alpha
energy 119.54 energy 119.54
#keyword; filename; number of cuts in file
cutFile cuts/cutsSQ13Alpha.root 2
cutSQ13Alpha30
cutSQ13Alpha30Amp
#channel minimum maximum
0 1450 1750 0 1450 1750
1 2200 2450 1 2200 2450
2 1900 2150 2 1900 2150
......
...@@ -34,10 +34,11 @@ void parTest() ...@@ -34,10 +34,11 @@ void parTest()
// cout << p.GetFileName(2) << endl; // cout << p.GetFileName(2) << endl;
// cout << p.GetFileName(3) << endl; // cout << p.GetFileName(3) << endl;
// cout << p.GetFileName(4) << endl; // cout << p.GetFileName(4) << endl;
// p.PrintParameters("");
p.PrintParameters("all"); p.PrintParameters("all");
// return; // return;
p.Reset(); p.Reset();
p.PrintParameters(""); // p.PrintParameters("");
// return; // return;
cout << "-----------------------" << endl; cout << "-----------------------" << endl;
...@@ -46,7 +47,7 @@ void parTest() ...@@ -46,7 +47,7 @@ void parTest()
c.SetParFile("SQ13Alpha.par"); c.SetParFile("SQ13Alpha.par");
c.Init(); c.Init();
c.PrintParameters(); c.PrintParameters();
return; // return;
// c.PrintParameters("all"); // c.PrintParameters("all");
...@@ -62,10 +63,11 @@ void parTest() ...@@ -62,10 +63,11 @@ void parTest()
// c.DrawVariable("SQ13", 2, c1); // c.DrawVariable("SQ13", 2, c1);
// c.DrawBeam(c1, 4, "SQ13"); // c.DrawBeam(c1, 4, "SQ13");
c.GetPeakMean("SQ13", 0, 0, c1, "cutSQ13Alpha16Amp", 256); // c.GetPeakMean("SQ13", 0, c1, "cutSQ13Alpha16Amp", 256);
c.GetPeakMean("SQ13", 1, 1, c2, "cutSQ13Alpha21Amp", 256); c.GetPeakMean("SQ13", 0, c1, "cutSQ13Ala16Amp", 256);
c.GetPeakMean("SQ13", 2, 2, c3, "cutSQ13Alpha26Amp", 256); c.GetPeakMean("SQ13", 1, c2, "cutSQ13Alpha21Amp", 256);
c.GetPeakMean("SQ13", 3, 3, c4, "cutSQ13Alpha30Amp", 256); c.GetPeakMean("SQ13", 2, c3, "cutSQ13Alpha26Amp", 256);
c.GetPeakMean("SQ13", 3, c4, "cutSQ13Alpha30Amp", 256);
//return; //return;
// c.PrintFiles(); // c.PrintFiles();
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment