From b2df0db46a276cef4a910d0fb01a8ee0e606b350 Mon Sep 17 00:00:00 2001 From: Vratislav Chudoba Date: Wed, 2 Nov 2016 10:11:43 +0300 Subject: [PATCH] Parameters for callibration of scintillator detectors have been somehow enhanced. --- AculCalib/AculCalPars.h | 10 +- AculCalib/AculCalParsScint.cpp | 120 +++++++++++++---- AculCalib/AculCalParsScint.h | 18 +-- AculCalib/AculCalParsScintFile.cpp | 150 ++++++++++++++++++++++ AculCalib/AculCalParsScintFile.h | 74 +++++++++++ AculCalib/AculCalParsSi.cpp | 200 +++++++++++++++++++++++++++++ AculCalib/AculCalParsSi.h | 74 +++++++++++ AculCalib/AculCalibScint.cpp | 62 ++++++--- AculCalib/AculCalibScint.h | 9 +- macros/myMacros/SQ13Alpha.par | 69 +++++----- macros/myMacros/parTest.cxx | 14 +- 11 files changed, 699 insertions(+), 101 deletions(-) create mode 100644 AculCalib/AculCalParsScintFile.cpp create mode 100644 AculCalib/AculCalParsScintFile.h create mode 100644 AculCalib/AculCalParsSi.cpp create mode 100644 AculCalib/AculCalParsSi.h diff --git a/AculCalib/AculCalPars.h b/AculCalib/AculCalPars.h index f98b7eb..2923434 100644 --- a/AculCalib/AculCalPars.h +++ b/AculCalib/AculCalPars.h @@ -51,10 +51,12 @@ public: virtual const char* GetFileName(Int_t i) {return 0;}; virtual Int_t GetNoEPoints() {return 0;} virtual Double_t GetCalEnergy(Int_t i) {return 0.;}; - virtual const char* GetCutsFileName() {return 0;} - virtual Int_t GetNoCuts() {return 0;} - virtual TCutG* GetCut(Int_t i) {return 0;}; - virtual TCutG* GetCut(const char* cutName) {return 0;}; + + //todo check functions with cuts +// virtual const char* GetCutsFileName() {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 GetMaxChannel(Int_t energy, Int_t crystal) {return 0;}; diff --git a/AculCalib/AculCalParsScint.cpp b/AculCalib/AculCalParsScint.cpp index 2465b6a..c537b50 100644 --- a/AculCalib/AculCalParsScint.cpp +++ b/AculCalib/AculCalParsScint.cpp @@ -27,7 +27,7 @@ AculCalParsScint::~AculCalParsScint() { void AculCalParsScint::Init() { SetPars(); - LoadCuts(); +// LoadCuts(); } void AculCalParsScint::SetPars() { @@ -88,7 +88,7 @@ void AculCalParsScint::SetPars() { continue; } - if ( line.BeginsWith("cutFile") ) { + /*if ( line.BeginsWith("cutFile") ) { sscanf(line.Data(), "%*s %s %d", fname, &fNoCuts); fCutsFileName = fname; for (Int_t j = 0; j < fNoCuts; j++) { @@ -101,7 +101,7 @@ void AculCalParsScint::SetPars() { fCutName.push_back(cname); } continue; - } + }*/ if ( line.BeginsWith("detector") ) { sscanf(line.Data(), "%*s %s %s", det, part); @@ -113,14 +113,61 @@ void AculCalParsScint::SetPars() { if ( line.BeginsWith("energy") ) { sscanf(line.Data(), "%*s %lf", &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++; continue; } +// cout << "aklsjdlaksjd" << endl; - sscanf(line.Data(), "%d %d %d", &i, &min, &max); - fFilePars.at(fEnergyPoints-1).SetPeakMin(i, min); - fFilePars.at(fEnergyPoints-1).SetPeakMax(i, max); +// sscanf(line.Data(), "%d %d %d", &i, &min, &max); +// fFilePars.at(fEnergyPoints-1).SetPeakMin(i, min); +// fFilePars.at(fEnergyPoints-1).SetPeakMax(i, max); }//while @@ -141,21 +188,27 @@ void AculCalParsScint::PrintParameters(const char* option) { } 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 << "\tInput file with " << fNoCuts << " cuts: \"" << fCutsFileName << "\"" << endl; - if (opt.Contains("all")) { - for (Int_t i = 0; i < (Int_t)fCutName.size(); i++) { - cout << "\t cut: \"" << fCutName[i] << "\"" << endl; - } - } +// cout << "\tInput file with " << fNoCuts << " cuts: \"" << fCutsFileName << "\"" << endl; +// if (opt.Contains("all")) { +// for (Int_t i = 0; i < (Int_t)fCutName.size(); i++) { +// 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++) { + 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; 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; @@ -167,6 +220,8 @@ void AculCalParsScint::PrintParameters(const char* option) { void AculCalParsScint::Reset() { + fNoCrystals = 0; + fParFileName = ""; fDetName = ""; @@ -176,10 +231,10 @@ void AculCalParsScint::Reset() { fFilePars.clear(); fEnergyPoints = 0; - fCutsFileName = ""; - fNoCuts = 0; //number of cuts +// fCutsFileName = ""; +// fNoCuts = 0; //number of cuts - fCutName.clear(); +// fCutName.clear(); return; } @@ -193,29 +248,40 @@ const char* AculCalParsScint::GetFileName(Int_t i) { 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 ) { cerr << "\"AculCalParsScint::GetCutName\" index i cannot be higher than " << fCutName.size() - 1 << endl; return 0; } return fCutName[i].Data(); -} +}*/ -TCutG* AculCalParsScint::GetCut(Int_t i) { +/*TCutG* AculCalParsScint::GetCut(Int_t i) { if ( i >= (Int_t)fCuts.size() ) { cerr << "\"AculCalParsScint::GetCut\" index i cannot be higher than " << fCuts.size() - 1 << endl; return 0; } return &fCuts[i]; -} +}*/ -TCutG* AculCalParsScint::GetCut(const char* cutName) { +TCutG* AculCalParsScint::GetCut(const char* cutName, Int_t treeID) { const TString cName = cutName; - for (Int_t i = 0; i <= (Int_t)fCuts.size(); i++) { - if (cName.EqualTo(fCutName[i])) { return &fCuts[i]; } + if (treeID >= (Int_t)fFilePars.size()) { + 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) { return fFilePars[energy].GetPeakMax(crystal); } -void AculCalParsScint::LoadCuts() { +/*void AculCalParsScint::LoadCuts() { if (fCutsFileName.Length() == 0) { printf("AculCalParsScint::LoadCuts: Name of file (*.root) with cuts was not provided.\n"); @@ -291,4 +357,4 @@ void AculCalParsScint::LoadCuts() { fCuts.push_back(*currentCut); } -} +}*/ diff --git a/AculCalib/AculCalParsScint.h b/AculCalib/AculCalParsScint.h index 5980e04..90b7b5e 100644 --- a/AculCalib/AculCalParsScint.h +++ b/AculCalib/AculCalParsScint.h @@ -26,13 +26,13 @@ private: vector fFilePars; //todo following parameters should be probably moved to AculCalParsScintFile //as far as they are related to each file separately - vector fCutName; - vector fCuts; +// vector fCutName; +// vector fCuts; Int_t fEnergyPoints; - TString fCutsFileName; - Int_t fNoCuts; //number of cuts +// TString fCutsFileName; +// Int_t fNoCuts; //number of cuts public: AculCalParsScint(); @@ -57,15 +57,15 @@ public: const char* GetCutName(Int_t i); Int_t GetNoEPoints() {return fEnergyPoints;} Double_t GetCalEnergy(Int_t i); - const char* GetCutsFileName() {return fCutsFileName.Data();} - Int_t GetNoCuts() {return fNoCuts;} - TCutG* GetCut(Int_t i); - TCutG* GetCut(const char* cutName); +// const char* GetCutsFileName() {return fCutsFileName.Data();} +// Int_t GetNoCuts() {return fNoCuts;} +// TCutG* GetCut(Int_t i); + TCutG* GetCut(const char* cutName, Int_t treeID = 0); Int_t GetMinChannel(Int_t energy, Int_t crystal); Int_t GetMaxChannel(Int_t energy, Int_t crystal); //private: - void LoadCuts(); +// void LoadCuts(); protected: void SetPars(); diff --git a/AculCalib/AculCalParsScintFile.cpp b/AculCalib/AculCalParsScintFile.cpp new file mode 100644 index 0000000..cce8173 --- /dev/null +++ b/AculCalib/AculCalParsScintFile.cpp @@ -0,0 +1,150 @@ +/* + * 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() { + +}*/ diff --git a/AculCalib/AculCalParsScintFile.h b/AculCalib/AculCalParsScintFile.h new file mode 100644 index 0000000..db5ad8d --- /dev/null +++ b/AculCalib/AculCalParsScintFile.h @@ -0,0 +1,74 @@ +/* + * AculCalParsScintFile.h + * + * Created on: Oct 25, 2016 + * Author: vratik + */ + +#ifndef ACULCALIB_ACULCALPARSSCINTFILE_H_ +#define ACULCALIB_ACULCALPARSSCINTFILE_H_ + +#include +#include + +#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 fCutName; + vector fCut; + + + vector fPeakMin; + vector 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_ */ diff --git a/AculCalib/AculCalParsSi.cpp b/AculCalib/AculCalParsSi.cpp new file mode 100644 index 0000000..c4d3f66 --- /dev/null +++ b/AculCalib/AculCalParsSi.cpp @@ -0,0 +1,200 @@ +/* + * 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(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(atof(parameter)); + } + continue; + }//if + + if ( strcmp(identificator, "lowerChannel") == 0 ) { + sscanf(line, "%s", parameter); + fLowerChannel = static_cast(atof(parameter)); + } + + if ( strcmp(identificator, "upperChannel") == 0 ) { + sscanf(line, "%s", parameter); + fUpperChannel = static_cast(atof(parameter)); + } + + 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)); + } + + if ( strcmp(identificator, "deadLayer") == 0 ) { + sscanf(line, "%s", parameter); + fDeadLayer = static_cast(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; +} diff --git a/AculCalib/AculCalParsSi.h b/AculCalib/AculCalParsSi.h new file mode 100644 index 0000000..c4a0e6a --- /dev/null +++ b/AculCalib/AculCalParsSi.h @@ -0,0 +1,74 @@ +/* + * 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 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_ */ diff --git a/AculCalib/AculCalibScint.cpp b/AculCalib/AculCalibScint.cpp index e48f604..0bed6fa 100644 --- a/AculCalib/AculCalibScint.cpp +++ b/AculCalib/AculCalibScint.cpp @@ -160,7 +160,9 @@ void AculCalibScint::DrawBeam(TCanvas *canvas, Int_t file, const char* variable) TString var; 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++) { canvas->cd(i+1); @@ -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"); // cout << "aksjda\t" << energyPoints << endl; for (Int_t j = 0; j < fPars->GetNoRawFiles(); j++) { - if ( fPars->GetCut(j) ) { - TCutG *cutTemp = fPars->GetCut(j); - curCutG = new TCutG(*cutTemp); +// if ( fPars->GetCut(j) ) { + // TCutG *cutTemp = fPars->GetCut(j); + // curCutG = new TCutG(*cutTemp); curCutG->Draw("same"); // printf("AculCalibCsI::DrawBeam: cTOF cut No. %d cannot be drawn, need to repair this function.\n", j); - } +// } } canvas->cd(noFiles+1+i); curTree->Draw("QDC[0]:QDC[1]", "", "cont"); for (Int_t j = 0; j < fPars->GetNoRawFiles(); j++) { - if ( fPars->GetCut(noFiles + j) ) { - TCutG *cutTemp = fPars->GetCut(noFiles + j); - curCutG = new TCutG(*cutTemp); +// if ( fPars->GetCut(noFiles + j) ) { +// TCutG *cutTemp = fPars->GetCut(noFiles + j); +// curCutG = new TCutG(*cutTemp); curCutG->Draw("same"); // printf("AculCalibCsI::DrawBeam: cQCD cut No. %d cannot be drawn, need to repair this function.\n", j); - } +// } } canvas->cd(2*noFiles+1+i); @@ -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 const Int_t noCrystals = 16; @@ -217,10 +219,22 @@ void AculCalibScint::GetPeakMean(const char* variable, Int_t treeID, Int_t energ TString con; TString hname; 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 - canvasTitle.Form("variable: %s; tree: %d; cut: %s;", variable, treeID, fPars->GetCut(beamcut)->GetName()); - canvas->SetTitle(canvasTitle.Data()); + if (currCut) { + canvasTitle.Form("variable: %s; tree: %d; cut: %s;", variable, treeID, currCut->GetName()); + canvas->SetTitle(canvasTitle.Data()); + } else { + canvasTitle.Form("variable: %s; tree: %d; no cut;", variable, treeID); + canvas->SetTitle(canvasTitle.Data()); + } TTree *curTree = 0; curTree = GetTree(treeID); @@ -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); - 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); hname.Form("hfull[%d][%d]", treeID, i); 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 curTree->Draw(var.Data(), con.Data()); var.Form("%s[%d]>>hcut[%d][%d]", variable, channel, treeID, i); - con.Form("%s[%d]>%d && %s[%d]<%d && %s", - variable, channel, fPars->GetMinChannel(energy, i), - variable, channel, fPars->GetMaxChannel(energy, i), - fPars->GetCut(beamcut)->GetName()); + if (currCut) { + con.Form("%s[%d]>%d && %s[%d]<%d && %s", + variable, channel, fPars->GetMinChannel(treeID, i), + variable, channel, fPars->GetMaxChannel(treeID, i), + 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); fHistCut[treeID].push_back( new TH1I(hname.Data(), "title", nbins, 0, 4096) ); fHistCut[treeID][i]->SetLineColor(3); diff --git a/AculCalib/AculCalibScint.h b/AculCalib/AculCalibScint.h index a6225fd..580e0eb 100644 --- a/AculCalib/AculCalibScint.h +++ b/AculCalib/AculCalibScint.h @@ -29,16 +29,11 @@ private: TTree **fTrees; TFile *fCutFile; -// TClonesArray cutsCol; -// -// TH1I *fHistFull[NOCALFILES][16]; + vector< vector > fHistFull; -// TH1I *fHistCut[NOCALFILES][16]; vector< vector > fHistCut; // -// Double_t fMeanPeakPos[NOCALFILES][16]; vector< vector > fMeanPeakPos; -// Double_t fMeanPeakRMS[NOCALFILES][16]; vector< vector > fMeanPeakRMS; // //todo make vector of graphs @@ -64,7 +59,7 @@ public: 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 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 WriteClbParameters(const char* filename); diff --git a/macros/myMacros/SQ13Alpha.par b/macros/myMacros/SQ13Alpha.par index 96df96a..c11a380 100644 --- a/macros/myMacros/SQ13Alpha.par +++ b/macros/myMacros/SQ13Alpha.par @@ -7,59 +7,61 @@ data/csi_13_Ealpha16.root data/csi_13_Ealpha21.root data/csi_13_Ealpha26.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 detector SQ13 Alpha #keyword energy [MeV/A] energy 66.542 +#keyword; filename; number of cuts in file +cutFile cuts/cutsSQ13Alpha.root 2 +cutsSQ13Alpha16 +cutSQ13Alpha16Amp #channel minimum maximum -0 650 830 -1 920 1150 -2 850 1000 +0 650 830 +1 920 1150 +2 850 1000 3 1200 1400 -4 960 1100 -5 700 850 -6 650 800 +4 960 1100 +5 700 850 +6 650 800 7 1040 1240 -8 850 1050 -9 650 850 -10 950 1200 -11 900 1100 -12 740 860 +8 850 1050 +9 650 850 +10 950 1200 +11 900 1100 +12 740 860 13 1140 1380 -14 420 500 -15 800 1000 +14 420 500 +15 800 1000 #CsI2, energy ~ 21 MeV/A, Alpha energy 85.437 +#keyword; filename; number of cuts in file +cutFile cuts/cutsSQ13Alpha.root 2 +cutSQ13Alpha21 +cutSQ13Alpha21Amp +#channel minimum maximum 0 1050 1200 1 1350 1550 2 1190 1400 3 1650 1950 4 1300 1500 -5 950 1100 -6 900 1050 +5 950 1100 +6 900 1050 7 1450 1700 8 1150 1400 -9 900 1100 +9 900 1100 10 1350 1600 11 1250 1500 -12 950 1200 +12 950 1200 13 1600 1860 -14 550 640 +14 550 640 15 1150 1350 #CsI2, energy ~ 26 MeV/A, Alpha energy 105.5 +#keyword; filename; number of cuts in file +cutFile cuts/cutsSQ13Alpha.root 2 +cutSQ13Alpha26 +cutSQ13Alpha26Amp +#channel minimum maximum 0 1300 1600 1 1750 2050 2 1560 1800 @@ -74,10 +76,15 @@ energy 105.5 11 1650 1900 12 1250 1600 13 2150 2450 -14 690 900 +14 690 900 15 1500 1700 #CsI2, energy ~ 30 MeV/A, Alpha energy 119.54 +#keyword; filename; number of cuts in file +cutFile cuts/cutsSQ13Alpha.root 2 +cutSQ13Alpha30 +cutSQ13Alpha30Amp +#channel minimum maximum 0 1450 1750 1 2200 2450 2 1900 2150 diff --git a/macros/myMacros/parTest.cxx b/macros/myMacros/parTest.cxx index f87c60e..f82aa3a 100644 --- a/macros/myMacros/parTest.cxx +++ b/macros/myMacros/parTest.cxx @@ -34,10 +34,11 @@ void parTest() // cout << p.GetFileName(2) << endl; // cout << p.GetFileName(3) << endl; // cout << p.GetFileName(4) << endl; +// p.PrintParameters(""); p.PrintParameters("all"); // return; p.Reset(); - p.PrintParameters(""); +// p.PrintParameters(""); // return; cout << "-----------------------" << endl; @@ -46,7 +47,7 @@ void parTest() c.SetParFile("SQ13Alpha.par"); c.Init(); c.PrintParameters(); - return; +// return; // c.PrintParameters("all"); @@ -62,10 +63,11 @@ void parTest() // c.DrawVariable("SQ13", 2, c1); // c.DrawBeam(c1, 4, "SQ13"); - c.GetPeakMean("SQ13", 0, 0, c1, "cutSQ13Alpha16Amp", 256); - c.GetPeakMean("SQ13", 1, 1, c2, "cutSQ13Alpha21Amp", 256); - c.GetPeakMean("SQ13", 2, 2, c3, "cutSQ13Alpha26Amp", 256); - c.GetPeakMean("SQ13", 3, 3, c4, "cutSQ13Alpha30Amp", 256); +// c.GetPeakMean("SQ13", 0, c1, "cutSQ13Alpha16Amp", 256); + c.GetPeakMean("SQ13", 0, c1, "cutSQ13Ala16Amp", 256); + c.GetPeakMean("SQ13", 1, c2, "cutSQ13Alpha21Amp", 256); + c.GetPeakMean("SQ13", 2, c3, "cutSQ13Alpha26Amp", 256); + c.GetPeakMean("SQ13", 3, c4, "cutSQ13Alpha30Amp", 256); //return; // c.PrintFiles(); -- 2.18.1